模拟退火算法研究

模拟退火算法的物理背景：

　　　　　　　　　　　　

模拟退火算法的核心思想与热力学的原理极为相似，在高温下，液体的大量分子彼此之间进行相对移动，如果液体慢慢冷却，原子的可动性就会消失，原子自行排列成行，形成一个纯净的晶体． 如果温度冷却迅速，那它不一定能达到这个状态，这一个过程的本质在于，慢慢冷却，使原子在丧失可动性之前重新排列，达到低能状态的必要条件．

模拟退火算法的描述：

　　　　　　　　　　　　　　　　　　　　　

模拟退火算法实质是一种贪心算法，但和贪心算法不同的是，模拟退火算法在搜索过程中加入了随机的因素，它在搜索的过程中，会以一定的概率来接受比当前解要更差的解，因此有可能会跳出这个局部最优解，找到全局的最优值．

如图，假设从A点出发，搜索过程中找到了局部最优解B,模拟退火算法会以一定的概率接受比B更差的解，从而找到全局最优解C

模拟退火算法流程图

 　　　　　　　　　　　　　　

　模拟退火算法主要分三个步骤：　　

　　１．在当前解的基础上产生一个新 解

　　２．利用Metropolis准则(即以e(- ∆E/T)的概率接受最优解)，判 断新解是否被接受． 　

　　３．当新解确定被接受时，用新解 代替当前解，在此基础上进行 下一轮实验，当新解被判定舍 弃时，在原解的基础上进行下 一轮实验

模拟退火算法伪代码：

模拟退火算法解决实际问题：

１．TSP问题：

旅行商问题(Traveling Salesman Problem,TSP),是数学领域一个著名问题之一，假设有一个旅行商人要拜访n个城市，他必须选择所要走的路径，路径的限制是每个城市只能拜访一次，而且最后要回到原来出发的城市，路径的选择目标是要求的路径的路程为所有路径中的最小值 .

执行代码

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <algorithm>
#include <iostream>
using namespace std;

int N;               //城市的数量
const double T=3000; //初始温度
const double k=0.98; //温度下降比例
const double mint=1e-8; //温度的极值

const int Oloop=1000;  //温度控制的最小值为t*0.98^20
const int Iloop=15000; //寻找一定温度下的最大值
const int lim=10000;  //概率选择次数的最大值

double dist[100][100]; //任意两座城市之间的距离

struct Path
{
    int city[100];
    double len;  //路径的总长度
};


struct Point    //城市点的坐标
{
    double x;
    double y;
};

Path bestpath;   //记录最优路径
Path currpath;   //记录当前的路径
Point p[100];     //各点的坐标

void input()
{
    scanf("%d",&N);
    for(int i=0;i<N;i++)
    {
        scanf("%lf%lf",&p[i].x,&p[i].y);
    }
}

//计算任意两点之间的距离
void compu()
{
    for(int i=0;i<N;i++)
    {
        for(int j=i+1;j<N;j++)
        {
            dist[i][j]=dist[j][i]=sqrt((p[i].x-p[j].x)*(p[i].x-p[j].x)+(p[i].y-p[j].y)*(p[i].y-p[j].y));
        }
    }
}

//初始化最优路径
void init()
{
    bestpath.len=0;
    bestpath.city[0]=0;
    for(int i=1;i<N;i++)
    {
        bestpath.city[i]=i;
        bestpath.len=bestpath.len+dist[i-1][i];
    }
    bestpath.len=bestpath.len+dist[0][N-1];
}

//随机选择目的优化路径，此处有两种方式，一是随机选取两个节点交换，二是随机排列
Path choose(Path newpath)
{
    int city1,city2;
    city1=(int)rand()%N;
    city2=(int)rand()%N;
    while(city1==city2)
    {
        city1=(int)rand()%N;
        city2=(int)rand()%N;
    }
    swap(newpath.city[city1],newpath.city[city2]);
    newpath.len=0;
    for(int i=1;i<N;i++)
    {
        newpath.len=newpath.len+dist[newpath.city[i-1]][newpath.city[i]];
    }
    newpath.len=newpath.len+dist[newpath.city[0]][newpath.city[N-1]];
    return newpath;
}

//模拟退火算法核心代码
Path sa()
{
    srand((unsigned)(time(NULL)));//时刻更新随机模板
    Path newpath=bestpath;
    Path currpath=bestpath;
    int pl=0; //以一定概率接受此路径的次数
    int pf=0;
    double t2=T;
    double t;
    while(true){
        for(int i=0;i<Iloop;i++)
        {
            newpath=choose(currpath);
            t=newpath.len-currpath.len;
            if(t<0)
            {
                currpath=newpath;
                pl=0;
                pf=0;

            }
            else
            {
                double ran=rand()/(RAND_MAX+1.0);
                double exx=-t*1.0/t2*1.0;
                double ex=exp(exx);
                if(ran<ex)
                {
                    currpath=newpath;
                }
                pl++;
            }
           if(pl>lim){
                pf++;
                break;
           }
        }
        if(currpath.len<bestpath.len)
        {
            bestpath=currpath;
        }
        if(pf>Oloop||t2<mint)
        {
            break;
        }
        t2=t2*k;
    }
}

void Print()
{
    for(int i=0;i<N;i++)
    {
        printf("%d ",bestpath.city[i]);
    }
    printf("length：%lf
",bestpath.len);
}

int main()
{
    for(int i=0;i<50;i++){
        freopen("in.txt","r",stdin);
        input();
        compu();
        init();
        sa();
        Print();
    }
    return 0;
}

验证数据：

27
41 94
37 84
53 67
25 62
7  64
2  99
68 58
71 44
54 62
83 69
64 60
18 54
22 60
83 46
91 38
25 38
24 42
58 69
71 71
74 78
87 76
18 40
13 40
82  7
62 32
58 35
45 21

执行过程

答案大概在401左右

２．求函数f(x,y)=5cos(xy)+xy+y^3的最小值，其中x的取值范围是(-5,5),y的取值范围是(-5,5)

执行代码：

进行过程中，容易出现陷入局部最优的问题

解决方法有：

#include <stdio.h>
#include <math.h>
#include <time.h>
#include <stdlib.h>
#include <algorithm>
#include <iostream>

using namespace std;

const double MAXx=5.0;          //x的最大值
const double MINx=-5.0;         //x的最小值
const double MAXy=5.0;          //y的最大值
const double MINy=-5.0;         //y的最小值
const double T=100;             //初始温度
const double S=0.5;             //每次的变化量
const double K=0.99;            //温度递减率
const double ex=1e-9;
double Bestx;
double Besty;

double functionf(double x1,double y1)       //所要求解的函数
{
    return 5*cos(x1*y1)+x1*y1+y1*y1*y1;
}

void sa()
{
    srand((unsigned)(time(NULL)));
    double currx,curry;
    double newx,newy;
    double Bestfun,currfun,newfun;
    double t=T;
    double de;
    Bestx=MINx+(MAXx-MINx)*rand()/RAND_MAX;
    Besty=MINy+(MAXx-MINy)*rand()/RAND_MAX;
    Bestfun=functionf(Bestx,Besty);
    currx=Bestx;
    curry=Besty;
    while(t>0.00001)
    {   int P=0;
        for(int i=0;i<10;i++)
        {
            int p=0;
            while(p==0)
            {
                newx=currx+S*(MINx+(MAXx-MINx)*rand()/RAND_MAX);
                if(newx>=MINx&&newx<=MAXx)
                {
                    p=1;
                }
            }
            de=functionf(currx,curry)-functionf(newx,newy);
            double dexp=exp(de/t);
            if(de>0)
            {
                currx=newx;
                //P++;
            }
            else
            {
                if(dexp>(rand()/RAND_MAX+1))
                {
                    currx=newx;
                    //P++;
                }
            }
            for(int j=0;j<10;j++)
            {
                p=0;
                while(p==0){
                    newy=curry+S*(MINy+(MAXy-MINy)*rand()/RAND_MAX);
                    if(newy>=MINy&&newy<=MAXy)
                    {
                        p=1;
                    }
                }
                de=functionf(currx,curry)-functionf(newx,newy);
                double dexp=exp(de/t);
                if(de>0)
                {
                    curry=newy;
                }
                else
                {
                    if(dexp>rand()/RAND_MAX+1)
                    {
                        curry=newy;
                    }
                }
            }



//            double y2=curry;
//            for(int j=0;j<10;j++)
//            {
//                p=0;
//                while(p==0){
//                    newy=y2+S*(MINy+(MAXy-MINy)*rand()/RAND_MAX);
//                    if(newy>=MINy&&newy<=MAXy)
//                    {
//                        p=1;
//                    }
//                }
//                de=functionf(currx,y2)-functionf(newx,newy);
//                double dexp=exp(de/t);
//                if(de>0)
//                {
//                    y2=newy;
//                }
//                else
//                {
//                    if(dexp>rand()/RAND_MAX+1)
//                    {
//                        y2=newy;
//                    }
//                }
//            }
//
//            if(functionf(currx,y2)<functionf(currx,curry))
//            {
//                curry=y2;
//            }
//            if(P>300){
//               P=0;
//                break;
//            }


        }
        if(functionf(currx,curry)<functionf(Bestx,Besty))
        {
                Bestx=currx;
                Besty=curry;
        }
        t=t*K;
    }
}

int main()
{
    for(int i=1;i<205;i++)
    {
        sa();
        printf("%f ",functionf(Bestx,Besty));
        if(i%6==0)
        printf("
");
    }
    printf("
");
}

 执行过程

答案在-151.1左右

总结：

模拟退火算法计算过程简单，引入概率跳出局部最优解，全局性较强． 但在执行过程中，算法的执行时间较长，相关参数的变化对算法性能影响大，例如，如果在某一温度下充分搜索，很容易陷入到局部最优解中，并且收敛速度慢，但如果不充分搜索，那么很难的到全局最优解． 所以，约束条件的选择，参数控制都非常关键，还有待进一步研究．