1 模拟退火算法（Simulated Annealing Algorithm）介绍
模拟退火算法是一种通用概率演算法，用来在一个大的搜索空间内寻找命题的最优解，它是基于Monte-Carlo迭代求解策略的一种随机寻优算法。
模拟退火算法来源于固体退火原理。
物理退火: 材料中的原子原来会停留在使内能有局部最小值的位置，加热使能量变大，原子会离开原来位置，而随机在其他位置中移动。退火冷却时速度较慢，使得原子有较多可能可以找到内能比原先更低的位置。
模拟火: 其原理也和固体退火的原理近似。模拟退火算法从某一较高初温出发，伴随温度参数的不断下降,结合概率突跳特性在解空间中随机寻找目标函数的全局最优解，即在局部最优解能概率性地跳出并最终趋于全局最优。
2 模拟退火算法（Simulated Annealing Algorithm）描述
模拟退火其实也是一种贪心算法，只不过与Local Search不同的是，模拟退火算法在搜索过程引入了随机因素。模拟退火算法以一定的概率来接受一个比当前解要差的解，因此有可能会跳出这个局部的最优解，达到全局的最优解。
从下图来说，模拟退火算法在搜索到局部最优解B后，会以一定的概率接受向右的移动。也许经过几次这样的不是局部最优的移动后会到达BC之间的峰点D，这样一来便跳出了局部最优解B，继续往右移动就有可能获得全局最优解C。如下图：
模拟退火算法示例
关于普通Greedy算法与模拟退火，这里也有一个有趣的比喻：
普通贪心算法：兔子朝着比现在低的地方跳去。它找到了不远处的最低的山谷。但是这座山谷不一定最低的。
模拟退火：兔子喝醉了。它随机地跳了很长时间。这期间，它可能走向低处，也可能踏入平地。但是，它渐渐清醒了并朝最低的方向跳去。
如此一来，大家对模拟退火算法有了一定的认识，但是这还是不够的。对比上面两种算法，对于模拟退火算法我们提到了一个很important的概念--一定的概率，关于这个一定的概率是如何计算的。这里还是参考了固体的物理退火过程。
2.1 模拟退火算法概率的确定
根据热力学的原理，在温度为T时，出现能量差为dE的降温的概率为P(dE)，表示为：
P(dE) = exp( dE/(kT) )
其中k是一个常数，exp表示自然指数，且dE<0(温度总是降低的)。这条公式指明了：
1) 温度越高，出现一次能量差为dE的降温的概率就越大。
2) 温度越低，则出现降温的概率就越小。又由于dE总是小于0（不然怎么叫退火），因此dE/kT < 0 ，exp(dE/kT)取值是(0,1),那么P(dE)的函数取值范围是(0,1) 。
随着温度T的降低，P(dE)会逐渐降低。我们将一次向较差解的移动看做一次温度跳变过程，我们以概率P(dE)来接受这样的移动。也就是说，在用固体退火模拟组合优化问题，将内能E模拟为目标函数值 f，温度T演化成控制参数 t，即得到解组合优化问题的模拟退火演算法：
由初始解 i 和控制参数初值 t 开始，对当前解重复“产生新解→计算目标函数差→接受或丢弃”的迭代，并逐步衰减 t 值，算法终止时的当前解即为所得近似最优解。
因此我们归结起来就是以下几点：
1) 若f( Y(i+1) ) <= f( Y(i) ) (即移动后得到更优解)，则总是接受该移动。
2) 若f( Y(i+1) ) > f( Y(i) ) (即移动后的解比当前解要差)，则以一定的概率接受移动，而且这个概率随着时间推移逐渐降低（逐渐降低才能趋向稳定）。
相当于上图中，从B移向BC之间的小波峰D时，每次右移(即接受一个更糟糕值)的概率在逐渐降低。如果这个坡特别长，那么很有可能最终我们并不会翻过这个坡。如果它不太长，这很有可能会翻过它，这取决于衰减 t 值的设定。
2.1 模拟退火算法伪代码
模拟退火算法伪代码
2.2 使用模拟退火算法解决旅行商问题
TSP是经典的NP完全问题。精确的解决TSP的算法的时间复杂度是O(2^N), 其中N是节点的个数 。而使用模拟退火算法则可以快速地获得一条近似最优路径。大体的思路如下：
1) 产生一条新的遍历路径P(i+1)，计算路径P(i+1)的长度L( P(i+1) )。
2) 若L(P(i+1)) < L(P(i))，则接受P(i+1)为新的路径，否则以模拟退火的那个概率接受P(i+1) ，然后降温。
3) 重复步骤1，2直到满足退出条件。
好了多说无益，下面大家一起看代码吧。
代码是以中国31个城市为例跑的。
/*
* 使用模拟退火算法(SA)求解TSP问题(以中国TSP问题为例)
* 参考自《Matlab 智能算法30个案例分析》
*/
#include<stdio.h>
#include<stdlib.h>
#include<string.h>
#include<time.h>
#include<math.h>
#define T0 50000.0 // 初始温度
#define T_end (1e-8)
#define q 0.98 // 退火系数
#define L 1000 // 每个温度时的迭代次数，即链长
#define N 31 // 城市数量
int city_list[N]; // 用于存放一个解
// 中国31个城市坐标
double city_pos[N][2] =
{
{1304,2312},{3639,1315},{4177,2244},{3712,1399},
{3488,1535},{3326,1556},{3238,1229},{4196,1004},
{4312,790},{4386,570},{3007,1970},{2562,1756},
{2788,1491},{2381,1676},{1332,695},
{3715,1678},{3918,2179},{4061,2370},
{3780,2212},{3676,2578},{4029,2838},
{4263,2931},{3429,1908},{3507,2367},
{3394,2643},{3439,3201},{2935,3240},
{3140,3550},{2545,2357},{2778,2826},
{2370,2975}};
//函数声明
double distance(double *,double *); // 计算两个城市距离
double path_len(int *); // 计算路径长度
void init(); //初始化函数
void create_new(); // 产生新解
// 距离函数
double distance(double * city1,double * city2)
{
double x1 = *city1;
double y1 = *(city1+1);
double x2 = *(city2);
double y2 = *(city2+1);
double dis = sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
return dis;
}
// 计算路径长度
double path_len(int * arr)
{
double path = 0; // 初始化路径长度
int index = *arr; // 定位到第一个数字(城市序号)
for(int i=0;i<N-1;i++)
{
int index1 = *(arr+i);
int index2 = *(arr+i+1);
double dis = distance(city_pos[index1-1],
city_pos[index2-1]);
path += dis;
}
int last_index = *(arr+N-1); // 最后一个城市序号
int first_index = *arr; // 第一个城市序号
double last_dis = distance(city_pos[last_index-1],
city_pos[first_index-1]);
path = path + last_dis;
return path; // 返回总的路径长度
}
// 初始化函数
void init()
{
for(int i=0;i<N;i++)
city_list[i] = i+1; // 初始化一个解
}
// 产生一个新解
// 此处采用随机交叉两个位置的方式产生新的解
void create_new()
{
double r1 = ((double)rand())/(RAND_MAX+1.0);
double r2 = ((double)rand())/(RAND_MAX+1.0);
int pos1 = (int)(N*r1); //第一个交叉点的位置
int pos2 = (int)(N*r2);
int temp = city_list[pos1];
city_list[pos1] = city_list[pos2];
city_list[pos2] = temp; // 交换两个点
}
// 主函数
int main(void)
{
srand((unsigned)time(NULL)); //初始化随机数种子
time_t start,finish;
start = clock(); // 程序运行开始计时
double T;
int count = 0; // 记录降温次数
T = T0; //初始温度
init(); //初始化一个解
int city_list_copy[N]; // 用于保存原始解
double f1,f2,df; //f1为初始解目标函数值，
//f2为新解目标函数值，df为二者差值
double r; // 0-1之间的随机数，用来决定是否接受新解
while(T > T_end) // 当温度低于结束温度时，退火结束
{
for(int i=0;i<L;i++)
{
// 复制数组
memcpy(city_list_copy,city_list,N*sizeof(int));
create_new(); // 产生新解
f1 = path_len(city_list_copy);
f2 = path_len(city_list);
df = f2 - f1;
// 以下是Metropolis准则
if(df >= 0)
{
r = ((double)rand())/(RAND_MAX);
if(exp(-df/T) <= r) // 保留原来的解
{
memcpy(city_list,city_list_copy,N*sizeof(int));
}
}
}
T *= q; // 降温
count++;
}
finish = clock(); // 退火过程结束
double duration = ((double)(finish-start))/CLOCKS_PER_SEC; // 计算时间
printf("模拟退火算法，初始温度T0=%.2f,降温系数q=%.2f,每个温度迭代%d次,共降温%d次，得到的TSP最优路径为:\n",T0,q,L,count);
for(int i=0;i<N-1;i++) // 输出最优路径
{
printf("%d--->",city_list[i]);
}
printf("%d\n",city_list[N-1]);
double len = path_len(city_list); // 最优路径长度
printf("最优路径长度为:%lf\n",len);
printf("程序运行耗时:%lf秒.\n",duration);
return 0;
}
网友评论