模拟退火算法（代码）

╰半夏微凉° 2022-11-30 04:21 213阅读 0赞

在实际日常中，人们会经常遇到如下问题：在某个给定的定义域![X][]内，求函数![f(x)][f_x]对应的最优值。此处以最小值问题举例（最大值问题可以等价转化成最小值问题），形式化为：

![\\min\_\{x \\in X\}f(x)][min_x _in X_f_x]

如果*![X][]*是离散有限取值，那么可以通过穷取法获得问题的最优解；如果*![X][]*连续，但![f(x)][f_x]是凸的，那可以通过梯度下降等方法获得最优解；如果*![X][]*连续且![f(x)][f_x]非凸，虽说根据已有的近似求解法能够找到问题解，可解是否是最优的还有待考量，很多时候若初始值选择的不好，非常容易陷入局部最优值。

![20200823164828511.png][]

随着日常业务场景的复杂化，第三种问题经常遇见。如何有效地避免局部最优的困扰？模拟退火算法应运而生。其实模拟退火也算是**启发式算法**的一种，具体学习的是冶金学中金属加热-冷却的过程。由S.Kirkpatrick, C.D.Gelatt和M.P.Vecchi在1983年所发明的，V.Čern在1985年也独立发明此演算法。

不过模拟退火算法到底是如何模拟金属退火的原理？主要是将热力学的理论套用到统计学上，将搜寻空间内每一点想象成空气内的分子；分子的能量，就是它本身的动能；而搜寻空间内的每一点，也像空气分子一样带有“能量”，以表示该点对命题的合适程度。演算法先以搜寻空间内一个任意点作起始：每一步先选择一个“邻居”，然后再计算从现有位置到达“邻居”的概率。若概率大于给定的阈值，则跳转到“邻居”；若概率较小，则停留在原位置不动。

**一、模拟退火算法基本思想**

模拟退火是启发示算法的一种，也是一种贪心算法，但是它的搜索过程引入了随机因素。在迭代更新可行解时，以一定的概率来接受一个比当前解要差的解，因此有可能会跳出这个局部的最优解，达到全局的最优解。以下图为例，假定初始解为左边蓝色点A，模拟退火算法会快速搜索到局部最优解B，但在搜索到局部最优解后，不是就此结束，而是会以一定的概率接受到左边的移动。也许经过几次这样的不是局部最优的移动后会到达全局最优点D，于是就跳出了局部最小值。

![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L0xpbnV4X2Jpbg_size_16_color_FFFFFF_t_70][]

根据热力学的原理，在温度为![T][]时，出现能量差为d*E*的降温的概率为![p(dE)][p_dE]，表示为：

![p(dE)=exp(\\frac\{dE\}\{kT\})][p_dE_exp_frac_dE_kT]

其中![k][]是波尔兹曼常数，值为![k=1.3806488\\times10^\{-23\}][k_1.3806488_times10_-23]，exp表示自然指数，且![dE<0][dE_0]。因此![dE/kT < 0][dE_kT _ 0]，所以![p(dE)][p_dE]函数的取值范围是(0, 1)。满足概率密度函数的定义。其实这条公式更直观意思就是：温度越高，出现一次能量差为![p(dE)][p_dE]的降温的概率就越大；温度越低，则出现降温的概率就越小。

在实际问题中，这里的“一定的概率”的计算参考了金属冶炼的退火过程。假定当前可行解为![x][]，迭代更新后的解为![x\_\{new\}][x_new]，那么对应的“能量差”定义为：

![\\Delta f = f(x\_\{new\}) - f(x)][Delta f _ f_x_new_ - f_x]

其对应的“一定概率”为：

最小值优化：![p(\\Delta f)=exp(-\\frac\{\\Delta f\}\{kT\})][p_Delta f_exp_-_frac_Delta f_kT]

最大值优化：![p(\\Delta f)=exp(\\frac\{\\Delta f\}\{kT\})][p_Delta f_exp_frac_Delta f_kT]

注：在实际问题中，可以设定![k=1][k_1]，即将参数*![k][]*与*![T][]*合并。

**二、模拟退火算法描述**

1.  初始化：初始温度*![T][]*（充分大），温度下限![T\_\{min\}][T_min]（充分小），初始解状态*![x][]*（是算法迭代的起点），每个*![T][]*值的迭代次数*![L][]*；
2.  对*![l][]*=1, 2, ..., *![L][]*做第3至第6步；
3.  产生新解*![x\_\{new\}][x_new]*: ![x\_\{new\}=x+\\Delta x][x_new_x_Delta x]；![\\Delta x][Delta x]为![\[d\_\{min\}, d\_\{max\}\]][d_min_ d_max]之间的随机数。
4.  利计算增量![\\Delta f = f(x\_\{new\}) - f(x)][Delta f _ f_x_new_ - f_x]，其中![f(x)][f_x]为优化目标；
5.  若![\\Delta f<0][Delta f_0](若寻找最大值，![\\Delta f>0][Delta f_0 1])则接受*![x\_\{new\}][x_new]*作为新的当前解，否则以概率![exp(-\\Delta f / (kT))][exp_-_Delta f _ _kT]接受*![x\_\{new\}][x_new]*作为新的当前解；
6.  如果满足终止条件则输出当前解作为最优解，结束程序。(终止条件通常取为连续若干个新解都没有被接受时终止算法。)；
7.  *T*逐渐减少，且![T>T\_\{min\}][T_T_min]，然后转第2步。

![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L0xpbnV4X2Jpbg_size_16_color_FFFFFF_t_70 1][]

在每个温度*![T][]*下迭代*![L][]*次，通过不断改变*x*来寻找当前温度下的最优值，然后降低温度继续寻找，直到温度达到最低，即选择概率![p][]接近于0。

**注意：**生成新的*![x][]*后，要判断是否在定义域内，对于超出的![x][]值要抛弃。

**三、模拟退火算法的优缺点**

模拟退火算法的应用很广泛，可以高效地求解NP完全问题，如货郎担问题(Travelling Salesman Problem，简记为TSP)、最大截问题(Max Cut Problem)、0-1背包问题(Zero One Knapsack Problem)、图着色问题(Graph Colouring Problem)等等，但其参数难以控制，不能保证一次就收敛到最优值，一般需要多次尝试才能获得（大部分情况下还是会陷入局部最优值）。观察模拟退火算法的过程，发现其主要存在如下三个参数问题：

**(1) ****温度*****T*****的初始值设置问题**

温度*T*的初始值设置是影响模拟退火算法全局搜索性能的重要因素之一、初始温度高，则搜索到全局最优解的可能性大，但因此要花费大量的计算时间；反之，则可节约计算时间，但全局搜索性能可能受到影响。

**(2) ****退火速度问题，即每个***T***值的迭代次数**

模拟退火算法的全局搜索性能也与退火速度密切相关。一般来说，同一温度下的“充分”搜索是相当必要的，但这也需要计算时间。循环次数增加必定带来计算开销的增大。

**(3) ****温度管理问题**

温度管理问题也是模拟退火算法难以处理的问题之一。实际应用中，由于必须考虑计算复杂度的切实可行性等问题，常采用如下所示的降温方式：

![T=\\alpha \\times T , \\alpha \\in (0, 1)][T_alpha _times T _ _alpha _in _0_ 1]

注：为了保证较大的搜索空间，*α*一般取接近于1的值，如0.95、0.9，代码中设置为0.98。

算法原理转自：[https://blog.csdn.net/huahua19891221/article/details/81737053][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;
    }

[X]: https://private.codecogs.com/gif.latex?X
[f_x]: https://private.codecogs.com/gif.latex?f%28x%29
[min_x _in X_f_x]: https://private.codecogs.com/gif.latex?%5Cmin_%7Bx%20%5Cin%20X%7Df%28x%29
[20200823164828511.png]: /images/20221124/fa3548bb0cdc4a1288671b1971e11c04.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L0xpbnV4X2Jpbg_size_16_color_FFFFFF_t_70]: /images/20221124/f43ef0654552455d90025df20d3ac450.png
[T]: https://private.codecogs.com/gif.latex?T
[p_dE]: https://private.codecogs.com/gif.latex?p%28dE%29
[p_dE_exp_frac_dE_kT]: https://private.codecogs.com/gif.latex?p%28dE%29%3Dexp%28%5Cfrac%7BdE%7D%7BkT%7D%29
[k]: https://private.codecogs.com/gif.latex?k
[k_1.3806488_times10_-23]: https://private.codecogs.com/gif.latex?k%3D1.3806488%5Ctimes10%5E%7B-23%7D
[dE_0]: https://private.codecogs.com/gif.latex?dE%3C0
[dE_kT _ 0]: https://private.codecogs.com/gif.latex?dE/kT%20%3C%200
[x]: https://private.codecogs.com/gif.latex?x
[x_new]: https://private.codecogs.com/gif.latex?x_%7Bnew%7D
[Delta f _ f_x_new_ - f_x]: https://private.codecogs.com/gif.latex?%5CDelta%20f%20%3D%20f%28x_%7Bnew%7D%29%20-%20f%28x%29
[p_Delta f_exp_-_frac_Delta f_kT]: https://private.codecogs.com/gif.latex?p%28%5CDelta%20f%29%3Dexp%28-%5Cfrac%7B%5CDelta%20f%7D%7BkT%7D%29
[p_Delta f_exp_frac_Delta f_kT]: https://private.codecogs.com/gif.latex?p%28%5CDelta%20f%29%3Dexp%28%5Cfrac%7B%5CDelta%20f%7D%7BkT%7D%29
[k_1]: https://private.codecogs.com/gif.latex?k%3D1
[T_min]: https://private.codecogs.com/gif.latex?T_%7Bmin%7D
[L]: https://private.codecogs.com/gif.latex?L
[l]: https://private.codecogs.com/gif.latex?l
[x_new_x_Delta x]: https://private.codecogs.com/gif.latex?x_%7Bnew%7D%3Dx&plus;%5CDelta%20x
[Delta x]: https://private.codecogs.com/gif.latex?%5CDelta%20x
[d_min_ d_max]: https://private.codecogs.com/gif.latex?%5Bd_%7Bmin%7D%2C%20d_%7Bmax%7D%5D
[Delta f_0]: https://private.codecogs.com/gif.latex?%5CDelta%20f%3C0
[Delta f_0 1]: https://private.codecogs.com/gif.latex?%5CDelta%20f%3E0
[exp_-_Delta f _ _kT]: https://private.codecogs.com/gif.latex?exp%28-%5CDelta%20f%20/%20%28kT%29%29
[T_T_min]: https://private.codecogs.com/gif.latex?T%3ET_%7Bmin%7D
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L0xpbnV4X2Jpbg_size_16_color_FFFFFF_t_70 1]: /images/20221124/08941d686d504b3d8751588336303a18.png
[p]: https://private.codecogs.com/gif.latex?p
[T_alpha _times T _ _alpha _in _0_ 1]: https://private.codecogs.com/gif.latex?T%3D%5Calpha%20%5Ctimes%20T%20%2C%20%5Calpha%20%5Cin%20%280%2C%201%29
[https_blog.csdn.net_huahua19891221_article_details_81737053]: https://blog.csdn.net/huahua19891221/article/details/81737053