DSA(direct simulated annealing，DSA)一种强大的随机搜索方法
author: san
email: visualsan@yahoo.cn
address: nuaa,cepe,202，2012.12
摘要：本文将探讨DSA优化算法在数值优化方面的应用，首先简单对比DSA和模拟退火算法，然后就DSA优化算法开发通用软件包，最后进行数值算例分析。
关键词：DSA,数值优化，模拟退火
1.DSA算法和模拟退火
2.DSA算法流程
3.DSA优化算法软件包
4.数值算例
5.结论
1.DSA算法和模拟退火
模拟退火算法是模拟金属高温退火过程的一种随机搜索算法。‘退火’的定义如下：将金属构件加热到高于或低于临界点，保持一定时间，随后缓慢冷却，从而获得接近平衡状态的组织与性能的金属热处理工艺。
模拟退火算法的概念：“模拟退火”算法是源于对热力学中退火过程的模拟，在某一给定初温下，通过缓慢下降温度参数，使算法能够在多项式时间内给出一个近似最优解。退火与冶金学上的‘退火’相似，而与冶金学的淬火有很大区别，前者是温度缓慢下降，后者是温度迅速下降。
算法流程图如下：
图1.SA算法流程图
其中T0为初始温度，一般取较高的数值；S0为随机选取的初始点；C(S)为目标函数值；K为在某一温度下的迭代次数，L为最大迭代次数，当K和L相等时，温度做一次下降处理，然后重新开始迭代；S’为随机产生的计算点；C(S’)为S’的目标函数值；t’为新产生的点的目标函数值和当前最佳点的目标函数值的差值，当t’小于0时，表示新的点比原来的点要好，算法接受新的点，若t’小于0，则S’有一点概念被接受，接受概念服从波尔兹曼概率分布，并且和温度有关，随着温度降低，S’被接受的概率降低，这表现为实际退火过程中温度较低后原子能量降低，跳动的概率降低；pa为接受概率，pr为产生的0和1之间的随机数；apha为温度降低指数，可取0.9到0.99之间的数。
2.DSA算法流程图
DSA为直接搜索模拟退火算法，是模拟退火算法的改进版本，该算法由Ali于2002年提出。该算法在两个方面区别于SA。首先在SA中，算法只维持一个当前最优点，而在DSA过程中，算法维持一个工作点集合，所以在SA中，算法只在一个点附近搜索，这使得SA可能会陷入局部收敛，而在DSA中，算法在一组工作点集合附近搜索，从而能有效地跳出局部最优点；其次在DSA优化过程中，一直更新和保持当前最佳工作点，这使得在优化过程中DSA拥有记忆功能，当一个新的工作点被接受时，它会取代工作组中的最差工作点。和SA相同的是，在DSA优化过程中，需要随机产生新工作点。
DSA算法流程图如图2所示：
图2.DSA算法流程图
3.DSA算法软件包开发
用面向对象语言C++开发DSA优化算法，主要类对象有
DSA：管理整个DSA算法流程
DSAConfigure：DSA优化过程中当前最优点集合管理类
DSAParmSet：DSA算法参数配置，包括初始温度T0以及一些迭代参数
DSAItem：对象类，DSA优化的中的一个工作点对象，类DSAItem为抽象类，定义了一些接口函数，具体类需要派生并实现接口。
DSAItemFactory：对象生产工厂，用于新的工作点的创建和解空间搜索。同样DSAItemFactory也是抽象类，具体类需要派生
DSAObserver：迭代观察器，用于输出每次迭代的计算结果，同样需要从该类派生出具体类
DSA类结构如下：
View Code
/*
直接搜索模拟退火算法（direct simulated annealing,DSA)。
测试例子里详细阐述了使用步骤，测试函数包括多变量函数优化、模式匹配、经典的GA测试函数、
TSP问题和一个九宫格匹配游戏。出了DSA外，程序还包含了通用的模拟退火算法（SA)和DSA的
并行版本。关于并行版本不建议使用，除非求解的问题的评估函数是耗时计算而且你所在的计算机
是多核计算机，因为对于普通的问题线程的创建和通信的时间可能会比你评估的时间还多，
关于DSA的算法原理可以自己查阅文献。
这是一个免费的程序，使用时不要删除作者的这段信息就可以了。求值结果请使用者自己斟酌，
作者概不负责！
author: SAN
data: 2011.7.25
address: 南京航空航天大学、能源与动力学院、202实验室、SAN
email: visualsan@yahoo.cn
------------------------------------------------------------------------------------*/
#include "define.h"
#include "DSAParmSet.h"
#include "DSAConfigure.h"
namespace DSA_OPT{
class DSA_API DSA
{
friend class DSAItem;
friend class DSAItemFactory;
friend class DSAResult;
friend class DSAObserver;
public:
//create rand number
static double GetRandVal(double min_=0,double max_=1);//获得随机数
DSA(DSAItemFactory*,DSAResult*r=0,DSAObserver*po=0);
virtual ~DSA();
//get item creater factory
DSAItemFactory* GetItemCreater();
//get observer
DSAObserver* GetObserver(){return m_pObserver;}
//get the length of marklow
long GetLK(){return m_LK;}
//get the actual times of marklow iteration
long GetLKActual(){return m_LK_ACTUAL;}
//reset 重置系统
virtual void Reset();
//SetRandSeed: set seed for rand function,when d=0 ,we use current time as seed
//初始化随机数种子(可选），d=0表示以当前时间为随机数种子，每次计算均不同
//对于需要重现结果的计算，建议自定义一个随机数种子。
void SetRandSeed(unsigned long d=0);
//get current temperature
double GetCurrentTemp();//获得当前温度
//get item configure
DSAConfigure& GetConfigure(){return m_Configure;}
//1.SET BASE PARM FOR DSA 1.设置基本信息
DSAParmSet m_InitInfo;
//2.INIT WORKSPACE 2.初始化
virtual void InitWorkspace();
//3.MAIN LOOP OF DSA 3.主循环
virtual void MainLoop();
//4.GET RESULT OF DSA 4.获得结果
DSAResult* GetResult(){return m_pIResult;}
//5.coerce stop 5.强制停止
void StopIter(){m_Stop=1;}
protected:
int m_bInit;
DSAConfigure m_Configure;
DSAItemFactory*m_pIItemCreater;
DSAObserver*m_pObserver;
DSAResult*m_pIResult;
double m_apha;//当前温度衰减系数
double m_T; //当前温度
int m_iter;//当前迭代次数
bool m_Stop;//强制停止迭代
//marklow链长度和实际计算的长度
int m_LK;//第k次计算的marklow链长度
int m_LK_ACTUAL;//第k次计算的marklow链实际使用长度
int m_LK_ACTUAL_pre;////第k-1次计算的marklow链实际使用长度
void CalMarklow();//计算marklow长度
//获得apha,初始状态apha=(m_apha_min+m_apha_max)/2
//以后apha由公式计算获得
void CalApha();
double GetApha(){return m_apha;}
//温度递减
void CoolTemp();
//终止条件
bool IsOver();
};
}
主要迭代步骤分如下：
1.设置基本信息
DSAParmSet m_InitInfo;
2.初始化
virtual void InitWorkspace();
3.主循环
virtual void MainLoop();
4.获得结果
DSAResult* GetResult(){return m_pIResult;}
其中主循环MainLoop执行DSA优化步骤如图2所示，主循环中主要调用了CalMarklow用于计算marklow链长度、CalApha用于计算冷却系数、CoolTemp用于冷却温度。主循环代码如下：
View Code
//主循环：visualsan@yahoo.cn
void DSA::MainLoop()
{
if(m_Stop)
return;
dbg("MainLoop\n");
int j=0;
int stop_flag=0;
int mk_stop=0;
m_iter=0;
m_pIItemCreater->UserOnStartRun(this);
do
{
m_pIItemCreater->UserTempControl(m_T,this);
j=0;
m_iter++;
mk_stop=0;
if(m_iter!=1)
{
m_LK_ACTUAL_pre=m_LK_ACTUAL;
}
dbg("SA iter=%d\n",m_iter);
CalMarklow();
//marklow loop
do
{
if(j==0)
{
dbg("marklow loop start\n");
}
j++;
dbg("marklow loop =%d\n",j);
DSAItem*p=m_pIItemCreater->CreateItem_CG_RG();
assert(p);
//计算目标函数
p->CalObj();
//统计函数求值次数
if(m_pIResult)
m_pIResult->m_TOTAL_COUNT++;
double pa=p->Acceptabiility();
//产生一个随机数
double rd=GetRandVal();
dbg("random=%f\n",rd);
if(pa>rd)//accept
{
dbg("pa>rd,accept\n");
//replace worst
m_Configure.ReplaceHigh(p);
//check if new obj is current best
if( p->ObjFunc()<m_Configure.GetLowItem()->ObjFunc() )
{
m_Configure.UpdateDataHighLow();
dbg("new obj is current best\n");
//set m_LK_ACTUAL=j
m_LK_ACTUAL=j;
//check if is over
stop_flag=m_pIItemCreater->UserStopCriterion()||m_Stop;
//stop marklow loop
mk_stop=1;
}
else
{
m_Configure.UpdateDataHighLow();
dbg("new obj is not current best\n");
if(j==m_LK)
{
m_LK_ACTUAL=j;
stop_flag=m_pIItemCreater->UserStopCriterion()||m_Stop;
//stop marklow loop
mk_stop=1;
}
}
}
else//reject
{
dbg("pa<rd,reject\n");
//delete object
delete p;
if(j==m_LK)
{
stop_flag=m_pIItemCreater->UserStopCriterion()||m_Stop;
//stop marklow loop
mk_stop=1;
m_LK_ACTUAL=j;
}
}
//如果有观察器，那么输出marklow循环信息
if(m_pObserver)
{
m_pObserver->RunOnMarkowLoop(j, //Markow内部第几次迭代
m_T, //当前温度
m_Configure.GetLowItem(), //当前最佳个体
m_Configure.GetHighItem()); //获取最差个体
}
} while (mk_stop==0);
//如果有观察器，则输出SA循环信息
if(m_pObserver)
{
m_pObserver->RunOnDSALoop(
m_iter, //退火过程第几次迭代
m_T, //当前温度
&m_Configure,//当前集合
m_Configure.GetLowItem(), //当前最佳个体
m_Configure.GetHighItem()); //获取最差个体
}
//保存结果
if(m_pIResult)
{
DSAResult::ITEM it;
it.bst=m_Configure.GetLowItem()->Copy();
it.wrst=m_Configure.GetHighItem()->Copy();
it.func_val_best=it.bst->ObjFunc();
it.func_val_worst=it.wrst->ObjFunc();
it.m_LK=m_LK;
it.m_LK_ACTUAL=m_LK_ACTUAL;
it.m_T=m_T;
m_pIResult->InsertItemBefore(it);
m_pIResult->m_ResultSet.push_back(it);
m_pIResult->InsertItemAfter(it);
}
if(stop_flag==0)
{
if(m_iter==1)
{
m_LK_ACTUAL_pre=m_LK_ACTUAL;
}
//计算冷却系数
CalApha();
//冷却
CoolTemp();
}
} while (stop_flag==0);
//退火迭代次数
if(m_pIResult)
m_pIResult->m_DSA_COUNT=m_iter;
if(m_pIResult)
{
time_t ltime;
m_pIResult->m_end_time=time( <ime );
m_pIResult->m_cpu_time=
(double)(clock()-m_pIResult->m_cpu_time)/CLOCKS_PER_SEC;
DSAResult::ITEM it;
it.bst=m_Configure.GetLowItem()->Copy();
it.wrst=m_Configure.GetHighItem()->Copy();
it.func_val_best=it.bst->ObjFunc();
it.func_val_worst=it.wrst->ObjFunc();
it.m_LK=m_LK;
it.m_LK_ACTUAL=m_LK_ACTUAL;
it.m_T=m_T;
m_pIResult->InsertItemBefore(it);
m_pIResult->m_ResultSet.push_back(it);
m_pIResult->InsertItemAfter(it);
m_pIResult->m_final=it;
}
m_pIItemCreater->UserOnEndRun(this);
}
类DSAConfigure主要用于保存迭代过程中的最优点集合，其结构如下
View Code
/*
直接搜索模拟退火算法（direct simulated annealing,DSA)。
测试例子里详细阐述了使用步骤，测试函数包括多变量函数优化、模式匹配、经典的GA测试函数、
TSP问题和一个九宫格匹配游戏。出了DSA外，程序还包含了通用的模拟退火算法（SA)和DSA的
并行版本。关于并行版本不建议使用，除非求解的问题的评估函数是耗时计算而且你所在的计算机
是多核计算机，因为对于普通的问题线程的创建和通信的时间可能会比你评估的时间还多，
关于DSA的算法原理可以自己查阅文献。
这是一个免费的程序，使用时不要删除作者的这段信息就可以了。求值结果请使用者自己斟酌，
作者概不负责！
author: SAN
data: 2011.7.25
address: 南京航空航天大学、能源与动力学院、202实验室、SAN
email: visualsan@yahoo.cn
------------------------------------------------------------------------------------*/
#include "define.h"
namespace DSA_OPT{
class DSA;
class DSAItem;
class DSA_API DSAConfigure
{
friend class DSAItem;
friend class DSA;
public:
DSAConfigure(DSA*);
virtual ~DSAConfigure();
//get current worst item pointer 当前目标函数值最大的个体
DSAItem* GetHighItem();
//get current best item pointer 当前目标函数值最大的个体
DSAItem* GetLowItem();
//replace the current worst item with new item 替换当前最佳个体
void ReplaceHigh(DSAItem*);
//get n random items from item set 随机获得n个集合中的个体
ItemVector GetRandItem(int n);
//get best n items from item set 获得最好的n个
ItemVector GetBestItem(int n);
//get random one depending on the scores,the lower has higher probability 根据目标函数值的大小选取一个，目标函数值小的选中概率大
DSAItem* GetItemWithFuncValue();
//get pointer to DSA algorithm
DSA* GetDSA(){return m_pDSA;}
private:
//init dsa algorithm
void Init();
//reset
void Reset();
//update the best item and worst item
void UpdateDataHighLow();
//pointer to DSA algorithm class
DSA*m_pDSA;
//m_ItemMap is a vector that keep a vector of the current best items,
//the size of m_ItemMap is determined by the size of problem.
//size of determine=7*(n+1),where n is the size of problem
ItemVector m_ItemMap;
//a pointer to item that has highest scores,always a
//higer scores item is worse than a lower one
DSAItem*m_pHigh;
//a pointer to item that has the lowest scores.
DSAItem*m_pLow;
//m_N is size of m_ItemMap,m_N=7*(n+1),where n is the size of problem
_Size m_N;
//compute acceptable probability of a item
double Factor(DSAItem*);
};
}
类DSAItem用于一个工作点的描述，工作点可以是一组数值，也可以是一组用于匹配的字符串，或者是一个对象。类DSAItem是一个抽象类，该类定义了一些接口，任何从该类派生的具体类都可以用于DSA算法的一个工作点。类DSAItem的接口如下：
//interface that must defined for a item 必须重写的函数
//how to cal object function 计算目标函数
virtual void CalObj()=0;
//how to cpoy an item 拷贝函数
virtual DSAItem* Copy()=0;
//a string that used to print item result 以字符串的形式获得变量值，用于结果显示
virtual _String GetVariableInString()=0;
接口CalObj用于计算目标函数；Copy定义拷贝自身；GetVariableInString接口用于获得一串字符串以描述当前的工作点的值。以函数优化为例，一个二元函数f(x1,x2)有两个优化变量x1和x2，所以工作点定义如下：
View Code
/*
直接搜索模拟退火算法（direct simulated annealing,DSA)。
测试例子里详细阐述了使用步骤，测试函数包括多变量函数优化、模式匹配、经典的GA测试函数、
TSP问题和一个九宫格匹配游戏。出了DSA外，程序还包含了通用的模拟退火算法（SA)和DSA的
并行版本。关于并行版本不建议使用，除非求解的问题的评估函数是耗时计算而且你所在的计算机
是多核计算机，因为对于普通的问题线程的创建和通信的时间可能会比你评估的时间还多，
关于DSA的算法原理可以自己查阅文献。
这是一个免费的程序，使用时不要删除作者的这段信息就可以了。求值结果请使用者自己斟酌，
作者概不负责！
author: SAN
data: 2011.7.25
address: 南京航空航天大学、能源与动力学院、202实验室、SAN
email: visualsan@yahoo.cn
------------------------------------------------------------------------------------*/
class fx:public DSAItem
{
public:
fx(DSAConfigure*p):DSAItem(p){}
//DEFINE COPY FUNCTION
virtual DSAItem* Copy()//拷贝函数
{
fx*f=new fx(m_pConfigure);
f->x1=x1;
f->x2=x2;
f->m_objfuncval=m_objfuncval;
return f;
}
//CAL OBJ FUNC
virtual void CalObj()
{//计算目标函数
m_objfuncval =x1*x2.......
}
//GET STRING OF VARIABLES 以字符串的形式获得变量值，用于结果显示
virtual _String GetVariableInString()
{
static char tmp[200];
sprintf(tmp,"x1=%f,x2=%f",x1,x2);
return tmp;
}
//变量值
double x1,x2;
}
DSA算法需要产生新的工作点，分两种情况：（1）初始化是随机产生工作点。（2）迭代过程中需要从最优点集合中进行搜索并产生新的工作点。迭代选取一般模式是：以概率p从最优的集合中随机选取一个点，概率1-p选取当前最好的点。选取后的点需要进行随机搜索，随机搜索的模式根据变量的类型的不同而不同。对于离散变量，则一般是轮盘随机选取；对于连续变量的话，则为加上一个随机数，数值大小根据搜索区间大小由用户自己设定。
连续变量搜索的原则是继承最好的同时进行新的尝试。继承最好表示从最优点集合中随机选取一个点进行搜索，而且有一定概率选择最好的点；新的尝试表示在选取点附近进行搜索。如何搜索是DSA算法的关键，而且不同的搜索方法得到的结果和计算效率是不一样的。搜索方法有很多，用户可以自行设置，管理对象创建的类是DSAItemFactory。该类有两个接口需要用户实现：初始化创建接口CreateItem_RAND()和随机搜索接口CreateItem_CG_RG。DSAItemFactory类如下所示
class DSA_API DSAItemFactory
{
friend class DSA;
public:
DSAItemFactory(){m_pIRandOperator=0;}
virtual ~DSAItemFactory(){if(m_pIRandOperator)delete m_pIRandOperator;}
/**
有两种方法实现随机操作：
1.重载函数CreateItem_CG_RG
2.重载类IRandOperator定义随机操作，实现IItem* RandOperator(Configure&c)
然后用函数void SetRandOper(IRandOperator*p)设置重载操作。
第二种方法能实现重载操作的动态载入，通过派生不同的重载操作可以在运行时
切换不同的操作。
*/
virtual DSAItem* CreateItem_CG_RG();
virtual DSAItem* CreateItem_RAND()=0;
//用户定义的停止准则，当系统满足用户定义停止准则
//时，算法终止，默认该停止准则为DSA算法定义的停止准则
virtual bool UserStopCriterion();
//用户对温度进行操作,默认不做任何操作
virtual void UserTempControl(double&T,DSA*){}
//通知用户开始优化
virtual void UserInit(DSA*){};
virtual void UserOnStartRun(DSA*){};
virtual void UserOnEndRun(DSA*){};
//设置随机操作方法，默认为NULL。
void SetRandOper(IRandOperator*p);
IRandOperator*GetRandOperator();
protected:
DSA*m_pDSA;
IRandOperator*m_pIRandOperator;
};
以优化函数f(x1,x2)为例，其创建类如下：
//DEFINE ITEM CREATER FACTORY
class cx:public DSAItemFactory
{
public:
virtual DSAItem* CreateItem_CG_RG()
{
DSAConfigure &c=m_pDSA->GetConfigure();
fx*f=(fx*)c.GetRandItem(1)[0];
f=(fx*)f->Copy();
fx*f2=(fx*)c.GetLowItem();
//RAND CHANGE NEW ITEM 搜索方法很多，在这里新的x为最好的点和随机点的
//平均再加上一个很小的随机数
f->x1=(f->x1+DSA::GetRandVal(-0.2,0.2)+f2->x1)/2.0;
f->x2=(f->x2+DSA::GetRandVal(-0.2,0.2)+f2->x2)/2.0;
return f;
}
virtual DSAItem* CreateItem_RAND()
{
fx*f=new fx(&m_pDSA->GetConfigure());
f->x1=DSA::GetRandVal(-1,1);
f->x2=DSA::GetRandVal(-1,1);
return f;
}
};
DSAObserver：输出迭代过程，需要迭代观察时请派生DSAObserver，重载接口函数实现迭代过程信息输出。DSAObserver结构如下：
View Code
#include "define.h"
namespace DSA_OPT{
class DSAItem;
class DSAConfigure;
class DSA_API DSAObserver
{
public:
friend class DSA;
//run RunOnMarkowLoop in every Markow iteration
virtual void RunOnMarkowLoop
(int index, //index of Markow iteration Markow内部第几次迭代
double T, //current temperature 当前温度
const DSAItem*bst, //current best 当前最佳个体
const DSAItem*wst){}; //current worst 当前最差个体
//run RunOnDSALoop in every dsa iteration
virtual void RunOnDSALoop
(int index, //index of DSA iteration 退火过程第几次迭代
double T, //current temperature 当前温度
const DSAConfigure*p, //current best items configure 当前集合
const DSAItem*bst, //current best item 当前最佳个体
const DSAItem*wst){}; //current worst item 当前最差个体
virtual~DSAObserver(){}
};
}
类DSAParmSet用于配置参数，DSA优化过程主要参数如下：
#include "define.h"
/*
class DSAParmSet is used to keep the base parameters of the dsa algorithm。
*/
namespace DSA_OPT{
class DSAParmSet
{
public:
DSAParmSet();
/*m_description:
a description string of the problem
*/
_String m_description;
/*
m_T0:
initial temperature of the dsa algorithm，default value is 500
*/
_Double m_T0;
/*
m_n:
number of the variables of the dsa probem,default 2
*/
_Size m_n;
/*
m_pRG:
the probability that create a new item with RG method,deafult 0.85
*/
_Double m_pRG;
/*m_apha_min,m_apha_max:
the range of cold factor,default [0.97,0.997]
*/
_Double m_apha_min,m_apha_max;
/*
termination criterion，deafult m_eps_1=0.008，m_eps_=0.008
*/
_Double m_eps_1,m_eps_2;
};
}
4.数值算例
(a)
优化函数 min f(x,y)=(x-2.0)*(x-2.0)+(y-2.0)*(y-2.0)
函数f(x,y)的最小值为f=0,(x,y)=(2,2)。使用DSA求解，首先是工作点类。该函数有两个参数，所以工作点类有两个成员变量x,y；搜索方式为：保留最优点，同时进行随机搜索。工作点类为：
View Code
class fx:public IItem
{
public:
fx(Configure*c):IItem(c){}
//define copy function
virtual IItem* Copy()//拷贝函数
{
fx*f=new fx(m_pConfigure);
f->x=x;
f->y=y;
f->m_objfuncval=m_objfuncval;
return f;
}
//define the objective value
virtual void CalObj()
{
m_objfuncval = (x-2.0)*(x-2.0)+(y-2.0)*(y-2.0);
}
//以字符串的形式获得变量值，用于结果显示
//return string which can be displayed during running
virtual std::string GetVariableInString()
{
static char tmp[200];
sprintf(tmp,"x=%f,y=%f",x,y);
return tmp;
}
//x,y data
double x,y;
}
创建工厂类为：
View Code
//define a creater whcih can producted rand value
class xc:public IItemCreater
{
public:
virtual IItem* CreateItem_CG_RG()
{
Configure &c=m_pDSA->GetConfigure();
fx*f=(fx*)c.GetRandItem(1)[0];
f=(fx*)f->Copy();
f->x=f->x+DSA::GetRandVal(-0.01,0.01);
f->y=f->y+DSA::GetRandVal(-0.01,0.01);
return f;
}
virtual IItem* CreateItem_RAND()
{
fx*f=new fx(&m_pDSA->GetConfigure());
f->x=DSA::GetRandVal(-2,2);
f->y=DSA::GetRandVal(-2,2);
return f;
}
};
主函数：
View Code
//solve
double range[]={-10.0,10.0};
fx::xc *x=new fx::xc;
DSA d(x);
d.InitWorkspace();
d.MainLoop();
cout<<"cost "<<d.GetResult()->GetCostTime()<<endl;
cout<<"(x,y)="<<d.GetConfigure().GetLowItem()->GetVariableInString().c_str()<<endl;
cout<<"f(x,y)="<<d.GetConfigure().GetLowItem()->ObjFunc()<<endl;
计算结果：
（b）min y = 21.5 + x1*sin(4*PI*x1) + x2*sin(20*PI*x2)
x1=[-3.0,12.1]
x2=[4.1,5.8]
设计变量x1,x2带有范围限制，所以采用罚函数的形式进行限制。首先用matlab计算，计算代码如下：
View Code
banana = @(x)21.5 + x(1)*sin(4*pi*x(1)) + x(2)*sin(20*pi*(2));
[x,fval,exitflag] = fminsearch(banana,[5,5])
=》
x1=4.8763 x2=5.1244 y=16.6244
工作点函数：
View Code
#define PI 3.1415926
//定义目标函数fx
class fx:public IItem
{
public:
fx(Configure*c):IItem(c){}
virtual IItem* Copy()//拷贝函数
{
fx*f=new fx(m_pConfigure);
f-> x1=x1;
f-> x2=x2;
f->m_objfuncval=m_objfuncval;
return f;
}
//计算目标函数
virtual void CalObj()
{
double a=0;
a=21.5 + x1*sin(4.0*PI*x1) + x2*sin(20.0*PI*x2);
if(x1<-3.0)
a=a+1000.0*(-3.0-x1);
else
if(x1>12.1)
a=a+1000.0*(x1-12.1);
if(x2<4.1)
a=a+1000.0*(4.1-x2);
else
if(x2>5.8)
a=a+1000.0*(x2-5.8);
m_objfuncval= a;
}
//以字符串的形式获得变量值，用于结果显示
virtual std::string GetVariableInString()
{
static char tmp[300];
sprintf(tmp,"x1=%f，x2=%f",
x1, x2);
return tmp;
}
//数据项
double x1,x2;
}
创建工厂：
View Code
//创建函数
class xc:public IItemCreater
{
public:
virtual IItem* CreateItem_CG_RG()
{
Configure &c=m_pDSA->GetConfigure();
fx*f=(fx*)c.GetDSA()->GetItemCreater()
->CreateItem_RAND(),*f2=(fx*)c.GetLowItem();
//精英保留策略
double p;
f->x1=f->x1*0.3+0.7*f2->x1;
f->x2=f->x2*0.3+0.7*f2->x2;
//随机突变策略
p=DSA::GetRandVal();
if(p<0.01)
f->x1+=DSA::GetRandVal(-0.1,0.1);
p=DSA::GetRandVal();
if(p<0.01)
f->x2+=DSA::GetRandVal(-0.1,0.1);
return f;
}
virtual IItem* CreateItem_RAND()
{
fx*f=new fx(&m_pDSA->GetConfigure());
f->x1=DSA::GetRandVal(-50,50);
f->x2=DSA::GetRandVal(-50,50);
return f;
}
bool UserStopCriterion()
{
return m_pDSA->GetCurrentTemp()<1e-50;
}
};
观察类：
View Code
class ob:public IObserver
{
public:
ob(){lp=0;}
int lp;
virtual void RunOnSALoop
(int index, //退火过程第几次迭代
double T, //当前温度
const Configure*p, //当前集合
const IItem*bst, //当前最佳个体
const IItem*wst) //获取最差个体
{
if(lp%300==0)
{
lp=0;
// fx*f=(fx*)((IItem*)bst);
// cout<<"SA ITER:"<<index<<",T="<<T<<",x1="<<f->x1<<",x2="<<f->x2
// <<":"<<((IItem*)bst)->ObjFunc()<<endl;
}lp++;
}
};
主函数：
View Code
fx::xc *x=new fx::xc;
fx::ob* o=new fx::ob;
DSA d(x,0,o);
d.m_InitInfo.m_n=2;//问题变量个数
d.m_InitInfo.m_T0=500;//初始化温度
d.m_InitInfo.m_eps_2=0.00001;
d.m_InitInfo.m_eps_1=0.00001;
d.InitWorkspace();
d.MainLoop();
cout<<"----------------------------result---------------------\n";
cout<<"cost "<<d.GetResult()->GetCostTime()<<endl;
cout<<"(x,y)="<<d.GetConfigure().GetLowItem()->GetVariableInString().c_str()<<endl;
cout<<"f(x,y)="<<d.GetConfigure().GetLowItem()->ObjFunc()<<endl;
完整代码：
View Code
#pragma warning(disable:4786)
#include "test.h"
#include <iostream>
#include "../src/AConfigure.h"
#include <math.h>
#include <algorithm>
#include "../src/DSA_PARALLEL.h"
#include "../src/Result.h"
#include "../src/DSA_PARALLEL.h"
#include <complex>
using namespace std;
//21.5 + v[0]*sin(4*PI*v[0]) + v[1]*sin(20*PI*v[1]);
void test::test_func3()
{
cout << "This program finds the maximum value in the function\n";
cout << " y = 21.5 + x1*sin(4*PI*x1) + x2*sin(20*PI*x2)\n";
cout << "with the constraints\n";
cout << " x1=[-3.0,12.1]\n";
cout << " x2=[4.1,5.8]\n";
cout << "\n\n";
#define PI 3.1415926
//定义目标函数fx
class fx:public IItem
{
public:
fx(Configure*c):IItem(c){}
virtual IItem* Copy()//拷贝函数
{
fx*f=new fx(m_pConfigure);
f-> x1=x1;
f-> x2=x2;
f->m_objfuncval=m_objfuncval;
return f;
}
//计算目标函数
virtual void CalObj()
{
double a=0;
a=21.5 + x1*sin(4.0*PI*x1) + x2*sin(20.0*PI*x2);
if(x1<-3.0)
a=a+1000.0*(-3.0-x1);
else
if(x1>12.1)
a=a+1000.0*(x1-12.1);
if(x2<4.1)
a=a+1000.0*(4.1-x2);
else
if(x2>5.8)
a=a+1000.0*(x2-5.8);
m_objfuncval= a;
}
//以字符串的形式获得变量值，用于结果显示
virtual std::string GetVariableInString()
{
static char tmp[300];
sprintf(tmp,"x1=%f，x2=%f",
x1, x2);
return tmp;
}
//数据项
double x1,x2;
//创建函数
class xc:public IItemCreater
{
public:
virtual IItem* CreateItem_CG_RG()
{
Configure &c=m_pDSA->GetConfigure();
fx*f=(fx*)c.GetDSA()->GetItemCreater()
->CreateItem_RAND(),*f2=(fx*)c.GetLowItem();
double p;
f->x1=f->x1*0.3+0.7*f2->x1;
f->x2=f->x2*0.3+0.7*f2->x2;
p=DSA::GetRandVal();
if(p<0.01)
f->x1+=DSA::GetRandVal(-0.1,0.1);
p=DSA::GetRandVal();
if(p<0.01)
f->x2+=DSA::GetRandVal(-0.1,0.1);
return f;
}
virtual IItem* CreateItem_RAND()
{
fx*f=new fx(&m_pDSA->GetConfigure());
f->x1=DSA::GetRandVal(-50,50);
f->x2=DSA::GetRandVal(-50,50);
return f;
}
bool UserStopCriterion()
{
return m_pDSA->GetCurrentTemp()<1e-50;
}
};
class ob:public IObserver
{
public:
ob(){lp=0;}
int lp;
virtual void RunOnSALoop
(int index, //退火过程第几次迭代
double T, //当前温度
const Configure*p, //当前集合
const IItem*bst, //当前最佳个体
const IItem*wst) //获取最差个体
{
if(lp%300==0)
{
lp=0;
// fx*f=(fx*)((IItem*)bst);
// cout<<"SA ITER:"<<index<<",T="<<T<<",x1="<<f->x1<<",x2="<<f->x2
// <<":"<<((IItem*)bst)->ObjFunc()<<endl;
}lp++;
}
};
};
fx::xc *x=new fx::xc;
fx::ob* o=new fx::ob;
DSA d(x,0,o);
d.m_InitInfo.m_n=2;//问题变量个数
d.m_InitInfo.m_T0=500;//初始化温度
d.m_InitInfo.m_eps_2=0.00001;
d.m_InitInfo.m_eps_1=0.00001;
d.InitWorkspace();
d.MainLoop();
cout<<"----------------------------result---------------------\n";
cout<<"cost "<<d.GetResult()->GetCostTime()<<endl;
cout<<"(x,y)="<<d.GetConfigure().GetLowItem()->GetVariableInString().c_str()<<endl;
cout<<"f(x,y)="<<d.GetConfigure().GetLowItem()->ObjFunc()<<endl;
}
计算结果：
(c)模式匹配，模式字符串为：
//模式匹配例子
char sed[]={' ','D','S','A'};
#define sz 50
#define h 15
char tag[]=
" DDDDDDDDDD SSSSSSSS AAAAA "
" DDDDDDDDDDD SSSSSSSSSS AAAAAAA "
" DDDDDDDDDDDD SSSSS SSSS AAAAAAA "
" DDDD DDDDD SSSS SSSSS AAAAAAA "
" DDDD DDDDD SSSSS AAAA AAAA "
" DDDD DDDD SSSSSSSS AAAA AAAA "
" DDDD DDDD SSSSSSSSS AAAA AAA "
" DDDD DDDD SSSSSSSS AAAAA AAAA "
" DDDD DDDD SSSSSS AAAA AAAA "
" DDDD DDDDD S SSSS AAAAAAAAAAA "
" DDDD DDDDD SSSS SSSS SAAAAAAAAAAAA "
" DDDDDDDDDDDD SSSSS SSSSS SAAAAAAAAAAAA "
" DDDDDDDDDDD SSSSSSSSSSS SAAA AAAAA "
" DDDDDDDDDD SSSSSSSSS SSAAA AAAAA "
" DDDDDDD SSSS SSAA AAAA ";
通过DSA算法进行匹配。该问题的工作点为由'D','S','A',' '组成的字符串，其工作点类为：
View Code
//定义目标函数fx
class fx:public IItem
{
public:
fx(Configure*c):IItem(c){}
virtual IItem* Copy()//拷贝函数
{
fx*f=new fx(m_pConfigure);
strcpy(f->pt,pt);
f->m_objfuncval=m_objfuncval;
return f;
}
virtual void CalObj()
{
double a=0;
for (int i=0;i<sz*h;i++)
{
if(tag[i]==pt[i])
a++;
}
m_objfuncval= (double)h*sz/a;
}
//以字符串的形式获得变量值，用于结果显示
virtual std::string GetVariableInString()
{
return pt;
}
//数据项
char pt[h*sz+1];
创建工厂类重写了终止准则，代码如下：
View Code
//创建函数
class xc:public IItemCreater
{
public:
virtual IItem* CreateItem_CG_RG()
{
Configure &c=m_pDSA->GetConfigure();
fx*f=(fx*)c.GetRandItem(1)[0],*f2=(fx*)c.GetLowItem();
f=(fx*)f->Copy();
int r;
double p;
for (int i=0;i<sz*h;i++)
{
r=rand()%4;
p=DSA::GetRandVal();
if(p<0.8)
f->pt[i]=f2->pt[i];
p=DSA::GetRandVal();
if(p<0.01)
f->pt[i]=sed[r];
}
return f;
}
virtual IItem* CreateItem_RAND()
{
fx*f=new fx(&m_pDSA->GetConfigure());
int r;
for (int i=0;i<sz*h;i++)
{
r=rand()%4;
f->pt[i]=sed[r];
}
f->pt[i]='\0';
return f;
}
//重写停止准则:要求绝对匹配
bool UserStopCriterion()
{
return
(m_pDSA->GetConfigure().GetLowItem()->ObjFunc()<1.06);
}
};
观察对象时时输出迭代结果：
View Code
class ob:public IObserver
{
public:
ob(){lp=0;}
int lp;
virtual void RunOnSALoop
(int index, //退火过程第几次迭代
double T, //当前温度
const Configure*p, //当前集合
const IItem*bst, //当前最佳个体
const IItem*wst) //获取最差个体
{
if(lp%80==0)
{ double v=((IItem*)bst)->ObjFunc();
cout<<v<<"---"<<((IItem*)wst)->ObjFunc()<<endl;
fx*f=(fx*)bst;
std::string s;
for (int i=0;i<h*sz;i++)
{
if( (i%sz)==0 && i!=0 )
s+="\n";
s+=f->pt[i];
}
cout<<s.c_str()<<endl;
lp=0;
}
lp++;
}
};
迭代过程：
（d）DeJongsFunction
DeJong's first function is: f1(x1,x2,x3) = x1*x1 + x2*x2 + x3*x3
where each x is in the range [-5.12,5.12]
DeJong's second function is:f2(x1,x2) = 100 * (x1*x1 - x2)^2 + (1 - x1)^2
where each x is in the range [-2.048,2.048
// DeJong's third function is:
//
// f3(x1,x2,x3,x4,x5) =
// 25 + floor(x1) + floor(x2) + floor(x3) + floor(x4) + floor(x5)
//
// where each x is in the range [-5.12,5.12]
// DeJong's fourth function is:
//
// 30
// ___
// f4(xi) = \ { i * xi^4 + Gauss(0,1) }
// /__
// i=1
//
// where each x is in the range [-1.28,1.28]
代码：
View Code
/*
直接搜索模拟退火算法（direct simulated annealing,DSA)。
测试例子里详细阐述了使用步骤，测试函数包括多变量函数优化、模式匹配、经典的GA测试函数、
TSP问题和一个九宫格匹配游戏。出了DSA外，程序还包含了通用的模拟退火算法（SA)和DSA的
并行版本。关于并行版本不建议使用，除非求解的问题的评估函数是耗时计算而且你所在的计算机
是多核计算机，因为对于普通的问题线程的创建和通信的时间可能会比你评估的时间还多，
关于DSA的算法原理可以自己查阅文献。
这是一个免费的程序，使用时不要删除作者的这段信息就可以了。求值结果请使用者自己斟酌，
作者概不负责！
author: SAN
data: 2011.7.25
address: 南京航空航天大学、能源与动力学院、202实验室、SAN
email: visualsan@yahoo.cn
------------------------------------------------------------------------------*/
#pragma warning(disable:4786)
#include "test.h"
#include <iostream>
#include "../src/AConfigure.h"
#include <math.h>
#include <algorithm>
#include "../src/DSA_PARALLEL.h"
#include "../src/Result.h"
#include "../src/DSA_PARALLEL.h"
using namespace std;
double
Gauss(double mean, double variance){
for(;;) {
double u1 = DSA::GetRandVal();
double u2 = DSA::GetRandVal();
double v1 = 2 * u1 - 1;
double v2 = 2 * u2 - 1;
double w = (v1 * v1) + (v2 * v2);
// cout<<w<<endl;
if (w <= 1) {
double y = sqrt( (-2 * log(w)) / w);
double x1 = v1 * y;
// double x2 = v2 * y; // we don't use this one
return(x1 * sqrt(variance) + mean);
}
}
}
// DeJong's first function is:
//
// f1(x1,x2,x3) = x1*x1 + x2*x2 + x3*x3
//
// where each x is in the range [-5.12,5.12]
float
DeJong1(double *x)
{
double x1=x[0];
double x2=x[1];
double x3=x[2];
float value=0;
value += x1* x1;
value += x2 * x2;
value += x3 * x3;
for (int i=0;i<3;i++)
{
if(fabs(x[i])-5.12>0)
value=value+1000.0*(fabs(x[i])-5.12);
}
return(value);
}
// DeJong's second function is:
//
// f2(x1,x2) = 100 * (x1*x1 - x2)^2 + (1 - x1)^2
//
// where each x is in the range [-2.048,2.048]
float
DeJong2(double*x)
{
double x1=x[0];
double x2=x[1];
float value=100.0;
value *= x1 * x1 - x2;
value *= x1 * x1 - x2;
value += (1 - x1)*(1 - x1);
for (int i=0;i<2;i++)
{
if(fabs(x[i])-2.048>0)
value=value+1000.0*(fabs(x[i])-2.048);
}
return(value);
}
// DeJong's third function is:
//
// f3(x1,x2,x3,x4,x5) =
// 25 + floor(x1) + floor(x2) + floor(x3) + floor(x4) + floor(x5)
//
// where each x is in the range [-5.12,5.12]
float
DeJong3(double *x)
{
double x1=x[0];
double x2=x[1];
double x3=x[2];
double x4=x[3];
double x5=x[4];
float value=25.0;
value += floor(x1);
value += floor(x2);
value += floor(x3);
value += floor(x4);
value += floor(x5);
for (int i=0;i<5;i++)
{
if(fabs(x[i])-5.12>0)
value=value+1000.0*(fabs(x[i])-5.12);
}
return(value);
}
// DeJong's fourth function is:
//
// 30
// ___
// f4(xi) = \ { i * xi^4 + Gauss(0,1) }
// /__
// i=1
//
// where each x is in the range [-1.28,1.28]
float
DeJong4(double *x)
{
float value = 0;
int i=0;
for( i=0; i<30; i++){
float v = x[i];
v *= v; // xi^2
v *= v; // xi^4
v *= i;
v += Gauss(0,1);
value += v;
}
for ( i=0;i<30;i++)
{
if(fabs(x[i])-1.28>0)
value=value+1000.0*(fabs(x[i])-1.28);
}
return(value);
}
void test::test_DeJongsFunction()
{
//定义目标函数fx
class fx:public IItem
{
public:
float (*pFunc)(double*);
fx(Configure*c):IItem(c){x=0;}
virtual~fx(){delete []x;}
virtual IItem* Copy()//拷贝函数
{
fx*f=new fx(m_pConfigure);
f->pFunc=pFunc;
f->len=len;
if(f->x)
delete [] f->x;
f->x= new double [len];
for (int i=0;i<len;i++)
{
f->x[i] = x[i];
}
f->m_objfuncval=m_objfuncval;
return f;
}
virtual void CalObj()
{
/*double a,b,c;
a=x*x+y*y;
b=pow(a,0.25);
c=pow(sin(50.0*pow(a,0.1)),2.0)+1;
m_objfuncval = b*c;*/
m_objfuncval = pFunc(x);
}
//以字符串的形式获得变量值，用于结果显示
virtual std::string GetVariableInString()
{
static char tmp[200];
std::string str;
for (int i=0;i<len;i++)
{
sprintf(tmp,"x%d=%f,",i+1,x[i]);
str=str+tmp;
}
return str;
}
//数据项
double *x;
int len;
//创建函数
class xc:public IItemCreater
{
public:
float (*pFunc)(double*);
int len;
virtual IItem* CreateItem_CG_RG()
{
Configure &c=m_pDSA->GetConfigure();
fx*f=(fx*)c.GetRandItem(1)[0];
f=(fx*)f->Copy();
fx*f2=(fx*)c.GetLowItem();
for (int i=0;i<f->len;i++)
{
f->x[i]=(0.3*f->x[i]+0.7*f2->x[i]);
if(DSA::GetRandVal()<0.01)
f->x[i]+=DSA::GetRandVal(-2,2);
}
/*
//when test second function,use this method
for (int i=0;i<f->len;i++)
{
f->x[i]=(0.3*f->x[i]+0.7*f2->x[i]);
if(DSA::GetRandVal()<0.01)
f->x[i]+=DSA::GetRandVal(-0.1,0.1);
}
*/
return f;
}
virtual IItem* CreateItem_RAND()
{
fx*f=new fx(&m_pDSA->GetConfigure());
f->pFunc=pFunc;
f->x=new double[len];
f->len=len;
for (int i=0;i<len;i++)
{
f->x[i]=DSA::GetRandVal(-10,10);
}
return f;
}
bool UserStopCriterion()
{
return m_pDSA->GetCurrentTemp()<1e-10;
}
};
class ob:public IObserver
{
public:
ob(){lp=0;}
int lp;
virtual void RunOnSALoop
(int index, //退火过程第几次迭代
double T, //当前温度
const Configure*p, //当前集合
const IItem*bst, //当前最佳个体
const IItem*wst) //获取最差个体
{
// cout<<"T="<<T<<","<<((IItem*)bst)->ObjFunc()<<",";
// cout<<((IItem*)wst)->ObjFunc()<<endl;
// cout<<((IItem*)bst)->GetVariableInString().c_str()<<endl;
}
};
};
// DeJong's first function is:
//
// f1(x1,x2,x3) = x1*x1 + x2*x2 + x3*x3
//
// where each x is in the range [-5.12,5.12]
{
fx::xc *x=new fx::xc;
x->pFunc=DeJong1;
x->len=3;
DSA d(x);
d.m_InitInfo.m_n=3;//问题变量个数
d.m_InitInfo.m_T0=500;//初始化温度
d.m_InitInfo.m_eps_2=0.01;
d.InitWorkspace();
d.MainLoop();
cout<<"f1(x1,x2,x3) = x1*x1 + x2*x2 + x3*x3\n";
cout<<"where each x is in the range [-5.12,5.12]\n";
cout<<"------------------------------------------------\n";
cout<<d.GetConfigure().GetLowItem()->GetVariableInString().c_str()<<endl;
cout<<"val="<<d.GetConfigure().GetLowItem()->ObjFunc()<<endl;
cout<<"------------------------------------------------\n";
}
// DeJong's second function is:
//
// f2(x1,x2) = 100 * (x1*x1 - x2)^2 + (1 - x1)^2
//
// where each x is in the range [-2.048,2.048]
{
fx::xc *x=new fx::xc;
x->pFunc=DeJong2;
x->len=2;
DSA d(x);
d.m_InitInfo.m_n=2;//问题变量个数
d.m_InitInfo.m_T0=5000;//初始化温度
d.m_InitInfo.m_eps_2=0.000001;
d.InitWorkspace();
d.MainLoop();
cout<<"DeJong's second function is:\n";
cout<<"f2(x1,x2) = 100 * (x1*x1 - x2)^2 + (1 - x1)^2\n";
cout<<"where each x is in the range [-2.048,2.048]\n";
cout<<"------------------------------------------------\n";
cout<<d.GetConfigure().GetLowItem()->GetVariableInString().c_str()<<endl;
cout<<"val="<<d.GetConfigure().GetLowItem()->ObjFunc()<<endl;
cout<<"------------------------------------------------\n";
}
// DeJong's third function is:
//
// f3(x1,x2,x3,x4,x5) =
// 25 + floor(x1) + floor(x2) + floor(x3) + floor(x4) + floor(x5)
//
// where each x is in the range [-5.12,5.12]
{
fx::xc *x=new fx::xc;
x->pFunc=DeJong3;
x->len=5;
DSA d(x);
d.m_InitInfo.m_n=5;//问题变量个数
d.m_InitInfo.m_T0=500;//初始化温度
d.m_InitInfo.m_eps_2=0.01;
d.InitWorkspace();
d.MainLoop();
cout<<"DeJong's third function is:\n";
cout<<"f3(x1,x2,x3,x4,x5) = 25 + floor(x1) + floor(x2) + floor(x3) + floor(x4) + floor(x5)\n";
cout<<"where each x is in the range [-5.12,5.12]\n";
cout<<"------------------------------------------------\n";
cout<<d.GetConfigure().GetLowItem()->GetVariableInString().c_str()<<endl;
cout<<"val="<<d.GetConfigure().GetLowItem()->ObjFunc()<<endl;
cout<<"------------------------------------------------\n";
}
// DeJong's fourth function is:
//
// 30
// ___
// f4(xi) = \ { i * xi^4 + Gauss(0,1) }
// /__
// i=1
//
// where each x is in the range [-1.28,1.28]
{
char hd[]=
"DeJong's fourth function is:\n"
" \n"
" 30\n"
" ___\n"
" f4(xi) = \ { i * xi^4 + Gauss(0,1) }\n"
" /__\n"
" i=1\n"
" \n"
" where each x is in the range [-1.28,1.28]\n";
fx::xc *x=new fx::xc;
fx::ob*o=new fx::ob;
x->pFunc=DeJong4;
x->len=30;
DSA d(x,0,o);
d.m_InitInfo.m_n=30;//问题变量个数
d.m_InitInfo.m_T0=500;//初始化温度
d.m_InitInfo.m_eps_2=0.1;
d.InitWorkspace();
d.MainLoop();
cout<<hd;
cout<<"------------------------------------------------\n";
cout<<d.GetConfigure().GetLowItem()->GetVariableInString().c_str()<<endl;
cout<<"val="<<d.GetConfigure().GetLowItem()->ObjFunc()<<endl;
cout<<"------------------------------------------------\n";
}
}
计算结果：
(e) tsp问题（Traveling Salesman Problem）
测试数据：
View Code
1 16.47 96.1
2 16.47 94.44
3 20.09 92.54
4 22.39 93.37
5 25.23 97.24
6 22 96.05
7 20.47 97.02
8 17.2 96.29
9 16.3 97.38
10 14.05 98.12
11 16.53 97.38
12 21.52 95.59
13 19.41 97.13
14 20.09 92.55
tsp问题代码
View Code
struct CITY
{
CITY(std::string n,double x_,double y_)
{
name=n;
x=x_;
y=y_;
}
std::string name;
CITY*pre,*next;
double x,y;//城市坐标
double dis(CITY*c)
{
return sqrt( (x-c->x)*(x-c->x) + (y-c->y)*(y-c->y) );
}
};
std::vector<CITY> ct;
int MAX_CITY=0;
void test::test_tsp()
{
//read data
char tmp[2000];
cout<<"Input File name:";
cin>>tmp;
cout<<endl;
FILE*f=fopen(tmp,"r");
if(!f)
{
cout<<"cant open:"<<tmp<<endl;
return;
}
char a[200],b[200],c[200];
fscanf(f,"%s",a);
fscanf(f,"%s",b);int var=atoi(b);
cout<<a<<"="<<b<<endl;
fscanf(f,"%s",a);
fscanf(f,"%s",b);double t0=atof(b);
cout<<a<<"="<<b<<endl;
while (feof(f)==0)
{
fscanf(f,"%s",a);
fscanf(f,"%s",b);
fscanf(f,"%s",c);
printf("%10s %10s %10s\n",a,b,c);
ct.push_back( CITY(a,atof(b),atof(c)));
}
MAX_CITY=ct.size();
fclose(f);
//定义目标函数fx
class fx:public IItem
{
public:
fx(Configure*c):IItem(c){}
virtual IItem* Copy()//拷贝函数
{
fx*f=new fx(m_pConfigure);
f->m_index=m_index;
f->m_objfuncval=m_objfuncval;
return f;
}
//计算目标函数
virtual void CalObj()
{
double a=0;
//统计城市距离
for (int i=0;i<m_index.size()-1;i++)
{
// printf("m_index[i]=%d,m_index[i+1]=%d\n",
// m_index[i],m_index[i+1]);
a=a+ct[ m_index[i] ].dis(&ct[ m_index[i+1] ]);
}
a=a+ct[ m_index[0] ].dis(&ct[ m_index[m_index.size()-1] ]);
m_objfuncval= a;
}
//以字符串的形式获得变量值，用于结果显示
virtual std::string GetVariableInString()
{
static char tmp[300];
std::string str;
for (int i=0;i<m_index.size();i++)
{
sprintf(tmp,"%d----",m_index[i]);
str+=tmp;
}
return str;
}
//数据项
std::vector<int> m_index;//城市顺序
//创建函数
class xc:public IItemCreater
{
public:
xc(){m_max_div=4;};
int m_max_div;//1-10最大段数
virtual IItem* CreateItem_CG_RG()
{
Configure &c=m_pDSA->GetConfigure();
fx*f;
if(DSA::GetRandVal()<0.3)
f=(fx*)c.GetItemWithFuncValue()->Copy();
else
f=(fx*)c.GetLowItem()->Copy();
int i=0;
//for (int i=0;i<f->m_index.size();i++)
{
//cout<<f->m_index[i]<<endl;
}
int r;
//随机产生m_max_div不相等的个整数
std::vector<int> d;
std::vector<int> val;
for ( i=1;i<f->m_index.size();i++)
{
d.push_back(i);
}
int k=m_max_div;
for (i=0;i<k;i++)
{
r=(rand()%d.size());
val.push_back( d[ r ] );
d.erase(d.begin()+r);
}
// cout<<"val==\n";
// for (i=0;i<val.size();i++)
{
// cout<<val[i]<<endl;
}
// cout<<"-------\n";
//随机连接
std::vector< std::vector<int> > lst;
d.clear();
std::sort(val.begin(),val.end());
for (i=0;i<val[0];i++)
{
d.push_back(f->m_index[i]);
}
lst.push_back(d);
for (i=1;i<val.size();i++)
{
d.clear();
for(int j=val[i-1];j<val[i];j++)
d.push_back( f->m_index[ j ] );
lst.push_back(d);
}
d.clear();
for (i=val[val.size()-1];i<f->m_index.size();i++)
{
d.push_back(f->m_index[i]);
}
lst.push_back(d);
f->m_index.clear();
while (lst.size()!=0)
{
int n=rand()%lst.size();
d=lst[n];
lst.erase(lst.begin()+n);
for (i=0;i<d.size();i++)
{
f->m_index.push_back(d[i]);
}
}
// for (i=0;i<f->m_index.size();i++)
{
// cout<<f->m_index[i]<<endl;
}
return f;
}
virtual IItem* CreateItem_RAND()
{
fx*f=new fx(&m_pDSA->GetConfigure());
std::vector<int> d;
for (int i=0;i<MAX_CITY;i++)
{
d.push_back(i);
}
int r;
//随机打乱排序
while( d.size()!=0 )
{
f->m_index.push_back( d[ r=(rand()%d.size()) ] );
d.erase(d.begin()+r);
}
return f;
}
};
class ob:public IObserver
{
public:
ob(){lp=0;}
int lp;
virtual void RunOnSALoop
(int index, //退火过程第几次迭代
double T, //当前温度
const Configure*p, //当前集合
const IItem*bst, //当前最佳个体
const IItem*wst) //获取最差个体
{
lp++;
if(lp%3000==0)
{
lp=0;
cout<<((IItem*)bst)->ObjFunc()*100.0<<"%"<<endl;
}
// cout<<((IItem*)bst)->GetVariableInString().c_str()<<endl;
}
};
};
fx::xc *x=new fx::xc;
fx::ob* o=new fx::ob;
DSA d(x,0,o);
if(var<=1)
var=1;
if(t0<=100)
t0=100;
d.m_InitInfo.m_n=var;//问题变量个数
d.m_InitInfo.m_T0=t0;//初始化温度
d.InitWorkspace();
d.MainLoop();
cout<<d.GetConfigure().GetLowItem()->GetVariableInString().c_str()<<endl;
cout<<"time cost:"<<d.GetResult()->GetCostTime()<<" sec"<<endl;
f=fopen("final_bst_wst_val.txt","w");
FILE*f2=fopen("final_result.txt","w");
for (int i=0;i<d.GetResult()->GetResultCount();i++)
{
fprintf(f,"%f %f\n",d.GetResult()->GetResultItem(i).func_val_best,
d.GetResult()->GetResultItem(i).func_val_worst);
fprintf(f2,"------------------iter%d--------------\n",i+1);
fprintf(f2,"%s\n",d.GetResult()->GetResultItem(i).bst->GetVariableInString().c_str());
}
fclose(f);
fclose(f2);
getchar();
}
计算结果：
来自tsplib的bier127的测试结果
5.结论
通过一系列函数的测试，得出：DSA算法收敛效率高，适用于离散变量和连续变量混合的优化问题。
1674
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
