在实际日常中,人们会经常遇到如下问题:在某个给定的定义域内,求函数对应的最优值。此处以最小值问题举例(最大值问题可以等价转化成最小值问题),形式化为:
如果是离散有限取值,那么可以通过穷取法获得问题的最优解;如果连续,但是凸的,那可以通过梯度下降等方法获得最优解;如果连续且非凸,虽说根据已有的近似求解法能够找到问题解,可解是否是最优的还有待考量,很多时候若初始值选择的不好,非常容易陷入局部最优值。
随着日常业务场景的复杂化,第三种问题经常遇见。如何有效地避免局部最优的困扰?模拟退火算法应运而生。其实模拟退火也算是启发式算法的一种,具体学习的是冶金学中金属加热-冷却的过程。由S.Kirkpatrick, C.D.Gelatt和M.P.Vecchi在1983年所发明的,V.Čern在1985年也独立发明此演算法。
不过模拟退火算法到底是如何模拟金属退火的原理?主要是将热力学的理论套用到统计学上,将搜寻空间内每一点想象成空气内的分子;分子的能量,就是它本身的动能;而搜寻空间内的每一点,也像空气分子一样带有“能量”,以表示该点对命题的合适程度。演算法先以搜寻空间内一个任意点作起始:每一步先选择一个“邻居”,然后再计算从现有位置到达“邻居”的概率。若概率大于给定的阈值,则跳转到“邻居”;若概率较小,则停留在原位置不动。
一、模拟退火算法基本思想
模拟退火是启发示算法的一种,也是一种贪心算法,但是它的搜索过程引入了随机因素。在迭代更新可行解时,以一定的概率来接受一个比当前解要差的解,因此有可能会跳出这个局部的最优解,达到全局的最优解。以下图为例,假定初始解为左边蓝色点A,模拟退火算法会快速搜索到局部最优解B,但在搜索到局部最优解后,不是就此结束,而是会以一定的概率接受到左边的移动。也许经过几次这样的不是局部最优的移动后会到达全局最优点D,于是就跳出了局部最小值。
根据热力学的原理,在温度为时,出现能量差为dE的降温的概率为,表示为:
其中是波尔兹曼常数,值为,exp表示自然指数,且。因此,所以函数的取值范围是(0, 1)。满足概率密度函数的定义。其实这条公式更直观意思就是:温度越高,出现一次能量差为的降温的概率就越大;温度越低,则出现降温的概率就越小。
在实际问题中,这里的“一定的概率”的计算参考了金属冶炼的退火过程。假定当前可行解为,迭代更新后的解为,那么对应的“能量差”定义为:
其对应的“一定概率”为:
最小值优化:
最大值优化:
注:在实际问题中,可以设定,即将参数与合并。
二、模拟退火算法描述
- 初始化:初始温度(充分大),温度下限(充分小),初始解状态(是算法迭代的起点),每个值的迭代次数;
- 对=1, 2, ..., 做第3至第6步;
- 产生新解: ;为之间的随机数。
- 利计算增量,其中为优化目标;
- 若(若寻找最大值,)则接受作为新的当前解,否则以概率接受作为新的当前解;
- 如果满足终止条件则输出当前解作为最优解,结束程序。(终止条件通常取为连续若干个新解都没有被接受时终止算法。);
- T逐渐减少,且,然后转第2步。
在每个温度下迭代次,通过不断改变x来寻找当前温度下的最优值,然后降低温度继续寻找,直到温度达到最低,即选择概率接近于0。
注意:生成新的后,要判断是否在定义域内,对于超出的值要抛弃。
模拟退火算法的应用很广泛,可以高效地求解NP完全问题,如货郎担问题(Travelling Salesman Problem,简记为TSP)、最大截问题(Max Cut Problem)、0-1背包问题(Zero One Knapsack Problem)、图着色问题(Graph Colouring Problem)等等,但其参数难以控制,不能保证一次就收敛到最优值,一般需要多次尝试才能获得(大部分情况下还是会陷入局部最优值)。观察模拟退火算法的过程,发现其主要存在如下三个参数问题:
(1) 温度T的初始值设置问题
温度T的初始值设置是影响模拟退火算法全局搜索性能的重要因素之一、初始温度高,则搜索到全局最优解的可能性大,但因此要花费大量的计算时间;反之,则可节约计算时间,但全局搜索性能可能受到影响。
(2) 退火速度问题,即每个T值的迭代次数
模拟退火算法的全局搜索性能也与退火速度密切相关。一般来说,同一温度下的“充分”搜索是相当必要的,但这也需要计算时间。循环次数增加必定带来计算开销的增大。
(3) 温度管理问题
温度管理问题也是模拟退火算法难以处理的问题之一。实际应用中,由于必须考虑计算复杂度的切实可行性等问题,常采用如下所示的降温方式:
注:为了保证较大的搜索空间,α一般取接近于1的值,如0.95、0.9,代码中设置为0.98。
算法原理转自:https://blog.csdn.net/huahua19891221/article/details/81737053。
四、案例
TSP问题即旅行商问题,经典的TSP可以描述为:一个商品推销员要去若干个城市推销商品,该推销员从一个城市出发,需要经过所有城市后,回到出发地。应如何选择行进路线,以使总的行程最短。从图论的角度来看,该问题实质是在一个带权完全无向图中,找一个权值最小的哈密尔顿回路。
分析:不同路线构求得的解构成了解空间,算法以一定的概率从解空间中迭代出极值,不同每次放入路线更新方式采用随机调整路线。
五、代码
#include <iostream>
#include <string.h>
#include <stdlib.h>
#include <algorithm>
#include <stdio.h>
#include <time.h>
#include <math.h>
#include <fstream>
#include <string>
#include <vector>
using namespace std;
#define T 3000 //初始温度
#define EPS 1e-8 //终止温度
#define DELTA 0.98 //温度衰减率
#define LIMIT 1000 //概率选择上限
#define OLOOP 20 //外循环次数
#define ILOOP 100 //内循环次数
#define cN 6 //城市数量
using namespace std;
//定义路线结构体
struct Path
{
int citys[cN];
double len;
};
//定义城市点坐标
struct Point
{
double x, y;
};
typedef struct Point pt;
Path bestPath; //记录最优路径
double w[cN][cN]; //两两城市之间路径长度
int nCase; //测试次数
double dist(Point A, Point B)
{
return sqrt((A.x - B.x) * (A.x - B.x) + (A.y - B.y) * (A.y - B.y));
}
void GetDist(std::vector<pt> p, int n)
{
for (int i = 0; i < n; i++)
for (int j = i + 1; j < n; j++)
w[i][j] = w[j][i] = dist(p[i], p[j]);//A到 B的距离等于B到A的距离。
}
void Init(int n)
{
nCase = 0;
bestPath.len = 0;
for (int i = 0; i < n; i++)
{
bestPath.citys[i] = i;//对条线设置编号
if (i != n - 1)
{
printf("%d--->", i);
bestPath.len += w[i][i + 1];
}
else
printf("%d\n", i);
}
printf("\n path length is : %.4f\n", bestPath.len);
}
void Print(Path t, int n)
{
printf("current path : ");
for (int i = 0; i < n; i++)
{
if (i != n - 1)
printf("%d-->", t.citys[i]);
else
printf("%d\n", t.citys[i]);
}
printf("\nThe path length is : %.3lf\n", t.len);
printf("-----------------------------------\n\n");
}
/*-----------------------------
生成新路线
------------------------------*/
Path GetcNext(Path p, int n)
{
Path ans = p;
int x = (int)(n * (rand() / (RAND_MAX + 1.0)));
int y = (int)(n * (rand() / (RAND_MAX + 1.0)));
while (x == y)
{
x = (int)(n * (rand() / (RAND_MAX + 1.0)));
y = (int)(n * (rand() / (RAND_MAX + 1.0)));
}
swap(ans.citys[x], ans.citys[y]);
ans.len = 0;
for (int i = 0; i < n - 1; i++)
ans.len += w[ans.citys[i]][ans.citys[i + 1]];
cout << "nCase = " << nCase << endl;
Print(ans, n);
nCase++;
return ans;
}
void SA(int n)
{
double t = T;
srand((unsigned)(time(NULL)));
Path curPath = bestPath;
Path newPath = bestPath;
int P_L = 0;
int P_F = 0;
while (1) //外循环,主要更新参数t,模拟退火过程
{
for (int i = 0; i < ILOOP; i++) //内循环,寻找在一定温度下的最优值
{
newPath = GetcNext(curPath, n);
double dE = newPath.len - curPath.len;
if (dE < 0) //当前值小于前一个值(求最小值问题),此时更新接受新值。
{
curPath = newPath;
P_L = 0;
P_F = 0;
}
else//以一定概率接受新值,可能会跳出局部最优解
{
double rd = rand() / (RAND_MAX + 1.0);
if (exp(dE / t) > rd && exp(dE / t) < 1)//概率大于随机生成的0-1之间的值时,接受新值。
curPath = newPath;
P_L++;
}
if (P_L > LIMIT)//终止条件
{
P_F++;
break;
}
}
if (curPath.len < bestPath.len)
bestPath = curPath;
if (P_F > OLOOP || t < EPS)
break;
t *= DELTA;//以一定速率更新温度
}
}
int main(int argc, const char * argv[]) {
ifstream f_data("C:/Users/91324/Desktop/Tsp.data");
std::vector<pt> city_pt;
city_pt.clear();
string data = "";
if (f_data.fail())//文件打开失败:返回0
{
std::cout << "Open file faild!" << std::endl;
return 0;
}
getline(f_data, data);
int n = atoi(data.c_str());
while (getline(f_data, data))
{
int pos=data.find(" ");
pt c_pt;
string x_str = data.substr(0, pos);
c_pt.x = atof(x_str.c_str());
string y_str = data.substr(pos+1, pos);
c_pt.y = atof(y_str.c_str());
city_pt.push_back(c_pt);
}
GetDist(city_pt, n);
Init(n);
SA(n);
Print(bestPath, n);
printf("Total test times is : %d\n", nCase);
system("pause");
return 0;
}