模拟退火算法来源于固体退火原理,将固体加温至充分高,再让其徐徐冷却,加温时,固体内部粒子随温升变为无序状,内能增大,而徐徐冷却时粒子渐趋有序,在每个温度都达到平衡态,最后在常温时达到基态,内能减为最小。根据Metropolis准则,粒子在温度T时趋于平衡的概率为e-ΔE/(kT),其中E为温度T时的内能,ΔE为其改变量,k为Boltzmann常数。用固体退火模拟组合优化问题,将内能E模拟为目标函数值f,温度T演化成控制参数t,即得到解组合优化问题的模拟退火算法:由初始解i和控制参数初值t开始,对当前解重复“产生新解→计算目标函数差→接受或舍弃”的迭代,并逐步衰减t值,算法终止时的当前解即为所得近似最优解,这是基于蒙特卡罗迭代求解法的一种启发式随机搜索过程。退火过程由冷却进度表(Cooling Schedule)控制,包括控制参数的初值t及其衰减因子Δt、每个t值时的迭代次数L和停止条件S。
模拟退火算法可以分解为解空间、目标函数和初始解三部分。
(1) 初始化:初始温度T(充分大),初始解状态S(是算法迭代的起点), 每个T值的迭代次数L
(2) 对k=1,……,L做第(3)至第6步:
(3) 产生新解S′
(4) 计算增量Δt′=C(S′)-C(S),其中C(S)为评价函数
(5) 若Δt′<0则接受S′作为新的当前解,否则以概率exp(-Δt′/T)接受S′作为新的当前解.
(6) 如果满足终止条件则输出当前解作为最优解,结束程序。 终止条件通常取为连续若干个新解都没有被接受时终止算法。
(7) T逐渐减少,且T->0,然后转第2步。
模拟退火的算法流程图如下:
模拟退火算法新解的产生和接受可分为如下四个步骤:
第一步是由一个产生函数从当前解产生一个位于解空间的新解;为便于后续的计算和接受,减少算法耗时,通常选择由当前新解经过简单地变换即可产生新解的方法,如对构成新解的全部或部分元素进行置换、互换等,注意到产生新解的变换方法决定了当前新解的邻域结构,因而对冷却进度表的选取有一定的影响。
第二步是计算与新解所对应的目标函数差。因为目标函数差仅由变换部分产生,所以目标函数差的计算最好按增量计算。事实表明,对大多数应用而言,这是计算目标函数差的最快方法。
第三步是判断新解是否被接受,判断的依据是一个接受准则,最常用的接受准则是Metropo1is准则: 若Δt′<0则接受S′作为新的当前解S,否则以概率exp(-Δt′/T)接受S′作为新的当前解S。
第四步是当新解被确定接受时,用新解代替当前解,这只需将当前解中对应于产生新解时的变换部分予以实现,同时修正目标函数值即可。此时,当前解实现了一次迭代。可在此基础上开始下一轮试验。而当新解被判定为舍弃时,则在原当前解的基础上继续下一轮试验。 模拟退火算法与初始值无关,算法求得的解与初始解状态S(是算法迭代的起点)无关;模拟退火算法具有渐近收敛性,已在理论上被证明是一种以概率l 收敛于全局最优解的全局优化算法;模拟退火算法具有并行性
如果你对退火的物理意义还是晕晕的,没关系我们还有更为简单的理解方式。想象一下如果我们现在有下面这样一个函数,现在想求函数的(全局)最优解。如果采用Greedy策略,那么从A点开始试探,如果函数值继续减少,那么试探过程就会继续。而当到达点B时,显然我们的探求过程就结束了(因为无论朝哪个方向努力,结果只会越来越大)。最终我们只能找打一个局部最后解B。
模拟退火其实也是一种Greedy算法,但是它的搜索过程引入了随机因素。模拟退火算法以一定的概率来接受一个比当前解要差的解,因此有可能会跳出这个局部的最优解,达到全局的最优解。以上图为例,模拟退火算法在搜索到局部最优解B后,会以一定的概率接受向右继续移动。也许经过几次这样的不是局部最优的移动后会到达B 和C之间的峰点,于是就跳出了局部最小值B。
根据Metropolis准则,粒子在温度T时趋于平衡的概率为exp(-ΔE/(kT)),其中E为温度T时的内能,ΔE为其改变数,k为Boltzmann常数。Metropolis准则常表示为
Metropolis准则表明,在温度为T时,出现能量差为dE的降温的概率为P(dE),表示为:P(dE) = exp( dE/(kT) )。其中k是一个常数,exp表示自然指数,且dE<0。所以P和T正相关。这条公式就表示:温度越高,出现一次能量差为dE的降温的概率就越大;温度越低,则出现降温的概率就越小。又由于dE总是小于0(因为退火的过程是温度逐渐下降的过程),因此dE/kT < 0 ,所以P(dE)的函数取值范围是(0,1) 。随着温度T的降低,P(dE)会逐渐降低。 我们将一次向较差解的移动看做一次温度跳变过程,我们以概率P(dE)来接受这样的移动。也就是说,在用固体退火模拟组合优化问题,将内能E模拟为目标函数值 f,温度T演化成控制参数 t,即得到解组合优化问题的模拟退火演算法:由初始解 i 和控制参数初值 t 开始,对当前解重复“产生新解→计算目标函数差→接受或丢弃”的迭代,并逐步衰减 t 值,算法终止时的当前解即为所得近似最优解,这是基于蒙特卡罗迭代求解法的一种启发式随机搜索过程。退火过程由冷却进度表(Cooling Schedule)控制,包括控制参数的初值 t 及其衰减因子Δt 、每个 t 值时的迭代次数L和停止条件S。
总结起来就是:
关于普通Greedy算法与模拟退火,有一个有趣的比喻:
作为模拟退火算法应用,讨论货郎担问题(Travelling Salesman Problem,简记为TSP):设有n个城市,用数码1,…,n代表。城市i和城市j之间的距离为d(i,j) i, j=1,…,n.TSP问题是要找遍访每个域市恰好一次的一条回路,且其路径总长度为最短.。
求解TSP的模拟退火算法模型可描述如下:
解空间 解空间S是遍访每个城市恰好一次的所有回路,是{1,……,n}的所有循环排列的集合,S中的成员记为(w1,w2 ,……,wn),并记wn+1= w1。初始解可选为(1,……,n)
目标函数 此时的目标函数即为访问所有城市的路径总长度或称为代价函数:
我们要求此代价函数的最小值。
新解的产生 随机产生1和n之间的两相异数k和m,若k<m,则将
(w1, w2 ,…,wk , wk+1 ,…,wm ,…,wn)
变为:
(w1, w2 ,…,wm , wm-1 ,…,wk+1 , wk ,…,wn).
如果是k>m,则将
(w1, w2 ,…,wk , wk+1 ,…,wm ,…,wn)
变为:
(wm, wm-1 ,…,w1 , wm+1 ,…,wk-1 ,wn , wn-1 ,…,wk).
上述变换方法可简单说成是“逆转中间或者逆转两端”。
也可以采用其他的变换方法,有些变换有独特的优越性,有时也将它们交替使用,得到一种更好方法。
代价函数差 设将(w1, w2 ,……,wn)变换为(u1, u2 ,……,un), 则代价函数差为:
根据上述分析,可写出用模拟退火算法求解TSP问题的伪程序:
Procedure TSPSA:
begin
init-of-T; { T为初始温度}
S={1,……,n}; {S为初始值}
termination=false;
while termination=false
begin
for i=1 to L do
begin
generate(S′form S); { 从当前回路S产生新回路S′}
Δt:=f(S′))-f(S);{f(S)为路径总长}
IF(Δt<0) OR (EXP(-Δt/T)>Random-of-[0,1])
S=S′;
IF the-halt-condition-is-TRUE THEN
termination=true;
End;
T_lower;
End;
End
下面给出C++实现参考源码:
1 /*
2 模拟退火算法解决TSP问题
3 输入格式(tsp.in):
4 第1行:1个整数N,表示城市的数量
5 第2..N+1行:每行有2个空格分开的整数x,y,第i+1行的x,y表示城市i的坐标
6 */
7 #include <iostream>
8 #include <string.h>
9 #include <stdlib.h>
10 #include <algorithm>
11 #include <stdio.h>
12 #include <time.h>
13 #include <math.h>
14
15 #define N 30 //城市数量
16 #define T 3000 //初始温度
17 #define EPS 1e-8 //终止温度
18 #define DELTA 0.98 //温度衰减率
19
20 #define LIMIT 1000 //概率选择上限
21 #define OLOOP 20 //外循环次数
22 #define ILOOP 100 //内循环次数
23
24 using namespace std;
25
26 //定义路线结构体
27 struct Path
28 {
29 int citys[N];
30 double len;
31 };
32
33 //定义城市点坐标
34 struct Point
35 {
36 double x, y;
37 };
38
39 Path bestPath; //记录最优路径
40 Point p[N]; //每个城市的坐标
41 double w[N][N]; //两两城市之间路径长度
42 int nCase; //测试次数
43
44 double dist(Point A, Point B)
45 {
46 return sqrt((A.x - B.x) * (A.x - B.x) + (A.y - B.y) * (A.y - B.y));
47 }
48
49 void GetDist(Point p[], int n)
50 {
51 for(int i = 0; i < n; i++)
52 for(int j = i + 1; j < n; j++)
53 w[i][j] = w[j][i] = dist(p[i], p[j]);
54 }
55
56 void Input(Point p[], int &n)
57 {
58 scanf("%d", &n);
59 for(int i = 0; i < n; i++)
60 scanf("%lf %lf", &p[i].x, &p[i].y);
61 }
62
63 void Init(int n)
64 {
65 nCase = 0;
66 bestPath.len = 0;
67 for(int i = 0; i < n; i++)
68 {
69 bestPath.citys[i] = i;
70 if(i != n - 1)
71 {
72 printf("%d--->", i);
73 bestPath.len += w[i][i + 1];
74 }
75 else
76 printf("%d\n", i);
77 }
78 printf("\nInit path length is : %.3lf\n", bestPath.len);
79 printf("-----------------------------------\n\n");
80 }
81
82 void Print(Path t, int n)
83 {
84 printf("Path is : ");
85 for(int i = 0; i < n; i++)
86 {
87 if(i != n - 1)
88 printf("%d-->", t.citys[i]);
89 else
90 printf("%d\n", t.citys[i]);
91 }
92 printf("\nThe path length is : %.3lf\n", t.len);
93 printf("-----------------------------------\n\n");
94 }
95
96 Path GetNext(Path p, int n)
97 {
98 Path ans = p;
99 int x = (int)(n * (rand() / (RAND_MAX + 1.0)));
100 int y = (int)(n * (rand() / (RAND_MAX + 1.0)));
101 while(x == y)
102 {
103 x = (int)(n * (rand() / (RAND_MAX + 1.0)));
104 y = (int)(n * (rand() / (RAND_MAX + 1.0)));
105 }
106 swap(ans.citys[x], ans.citys[y]);
107 ans.len = 0;
108 for(int i = 0; i < n - 1; i++)
109 ans.len += w[ans.citys[i]][ans.citys[i + 1]];
110 cout << "nCase = " << nCase << endl;
111 Print(ans, n);
112 nCase++;
113 return ans;
114 }
115
116 void SA(int n)
117 {
118 double t = T;
119 srand((unsigned)(time(NULL)));
120 Path curPath = bestPath;
121 Path newPath = bestPath;
122 int P_L = 0;
123 int P_F = 0;
124 while(1) //外循环,主要更新参数t,模拟退火过程
125 {
126 for(int i = 0; i < ILOOP; i++) //内循环,寻找在一定温度下的最优值
127 {
128 newPath = GetNext(curPath, n);
129 double dE = newPath.len - curPath.len;
130 if(dE < 0) //如果找到更优值,直接更新
131 {
132 curPath = newPath;
133 P_L = 0;
134 P_F = 0;
135 }
136 else
137 {
138 double rd = rand() / (RAND_MAX + 1.0);
139 //如果找到比当前更差的解,以一定概率接受该解,并且这个概率会越来越小
140 if(exp(dE / t) > rd && exp(dE / t) < 1)
141 curPath = newPath;
142 P_L++;
143 }
144 if(P_L > LIMIT)
145 {
146 P_F++;
147 break;
148 }
149 }
150 if(curPath.len < bestPath.len)
151 bestPath = curPath;
152 if(P_F > OLOOP || t < EPS)
153 break;
154 t *= DELTA;
155 }
156 }
157
158 int main(int argc, const char * argv[]) {
159
160 freopen("TSP.data", "r", stdin);
161 int n;
162 Input(p, n);
163 GetDist(p, n);
164 Init(n);
165 SA(n);
166 Print(bestPath, n);
167 printf("Total test times is : %d\n", nCase);
168 return 0;
169 }
TSP.data的数据格式如下,第一行的数字表示一个有多少座城市,第2至最后一行,每行有两个数字表示,城市的坐标(平面直角坐标系)。例如: 6 20 80 16 84 23 66 62 90 11 9 35 28
注意由于是基于蒙特卡洛的方法,所以上面代码每次得出的结果并不完全一致。你可以通过增加迭代的次数来获得一个更优的结果。
我们这里需要说明的是,在之前的文章里,我们用求最小值的例子来解释模拟退火的执行:如果新一轮的计算结果更前一轮之结果更小,那么我们就接受它,否则就以一个概率来拒绝或接受它,而这个拒绝的概率会随着温度的降低(也即是迭代次数的增加)而变大(也就是接受的概率会越来越小)。
但现在我们面对一个TSP问题,我们如何定义或者说如何获取下一轮将要被考察的哈密尔顿路径呢?在一元函数最小值的例子中,下一轮就是指向左或者向右移动一小段距离。而在TSP问题中,我们可以采用的方式其实是很多的。上面代码中GetNext()函数所采用的方式是随机交换两个城市在路径中的顺序。例如当前路径为 A->B->C->D->A,那么下一次路径就可能是A->D->C->B->A,即交换B和D。
1 public class Tour{
2 ... ...
3 // Creates a random individual
4 public void generateIndividual() {
5 // Loop through all our destination cities and add them to our tour
6 for (int cityIndex = 0; cityIndex < TourManager.numberOfCities(); cityIndex++) {
7 setCity(cityIndex, TourManager.getCity(cityIndex));
8 }
9 // Randomly reorder the tour
10 Collections.shuffle(tour);
11 }
12 ... ...
13 }
可见把上一轮路径做一个随机的重排(这显然也是一种策略)。
我们对上述问题提出一种新的策略:
首先,我们需要创建一个城市类,它可以用来为旅行推销员的不同目的地建模。
1 /*
2 * City.java
3 * Models a city
4 */
5
6 package sa;
7
8 public class City {
9 int x;
10 int y;
11
12 // Constructs a randomly placed city
13 public City(){
14 this.x = (int)(Math.random()*200);
15 this.y = (int)(Math.random()*200);
16 }
17
18 // Constructs a city at chosen x, y location
19 public City(int x, int y){
20 this.x = x;
21 this.y = y;
22 }
23
24 // Gets city's x coordinate
25 public int getX(){
26 return this.x;
27 }
28
29 // Gets city's y coordinate
30 public int getY(){
31 return this.y;
32 }
33
34 // Gets the distance to given city
35 public double distanceTo(City city){
36 int xDistance = Math.abs(getX() - city.getX());
37 int yDistance = Math.abs(getY() - city.getY());
38 double distance = Math.sqrt( (xDistance*xDistance) + (yDistance*yDistance) );
39
40 return distance;
41 }
42
43 @Override
44 public String toString(){
45 return getX()+", "+getY();
46 }
47 }
接下来让我们创建一个可以跟踪城市的类:
1 /*
2 * TourManager.java
3 * Holds the cities of a tour
4 */
5
6 package sa;
7
8 import java.util.ArrayList;
9
10 public class TourManager {
11
12 // Holds our cities
13 private static ArrayList destinationCities = new ArrayList<City>();
14
15 // Adds a destination city
16 public static void addCity(City city) {
17 destinationCities.add(city);
18 }
19
20 // Get a city
21 public static City getCity(int index){
22 return (City)destinationCities.get(index);
23 }
24
25 // Get the number of destination cities
26 public static int numberOfCities(){
27 return destinationCities.size();
28 }
29
30 }
现在来创建一个可以模拟旅行推销员之旅:
1 /*
2 * Tour.java
3 * Stores a candidate tour through all cities
4 */
5
6 package sa;
7
8 import java.util.ArrayList;
9 import java.util.Collections;
10
11 public class Tour{
12
13 // Holds our tour of cities
14 private ArrayList tour = new ArrayList<City>();
15 // Cache
16 private int distance = 0;
17
18 // Constructs a blank tour
19 public Tour(){
20 for (int i = 0; i < TourManager.numberOfCities(); i++) {
21 tour.add(null);
22 }
23 }
24
25 // Constructs a tour from another tour
26 public Tour(ArrayList tour){
27 this.tour = (ArrayList) tour.clone();
28 }
29
30 // Returns tour information
31 public ArrayList getTour(){
32 return tour;
33 }
34
35 // Creates a random individual
36 public void generateIndividual() {
37 // Loop through all our destination cities and add them to our tour
38 for (int cityIndex = 0; cityIndex < TourManager.numberOfCities(); cityIndex++) {
39 setCity(cityIndex, TourManager.getCity(cityIndex));
40 }
41 // Randomly reorder the tour
42 Collections.shuffle(tour);
43 }
44
45 // Gets a city from the tour
46 public City getCity(int tourPosition) {
47 return (City)tour.get(tourPosition);
48 }
49
50 // Sets a city in a certain position within a tour
51 public void setCity(int tourPosition, City city) {
52 tour.set(tourPosition, city);
53 // If the tours been altered we need to reset the fitness and distance
54 distance = 0;
55 }
56
57 // Gets the total distance of the tour
58 public int getDistance(){
59 if (distance == 0) {
60 int tourDistance = 0;
61 // Loop through our tour's cities
62 for (int cityIndex=0; cityIndex < tourSize(); cityIndex++) {
63 // Get city we're traveling from
64 City fromCity = getCity(cityIndex);
65 // City we're traveling to
66 City destinationCity;
67 // Check we're not on our tour's last city, if we are set our
68 // tour's final destination city to our starting city
69 if(cityIndex+1 < tourSize()){
70 destinationCity = getCity(cityIndex+1);
71 }
72 else{
73 destinationCity = getCity(0);
74 }
75 // Get the distance between the two cities
76 tourDistance += fromCity.distanceTo(destinationCity);
77 }
78 distance = tourDistance;
79 }
80 return distance;
81 }
82
83 // Get number of cities on our tour
84 public int tourSize() {
85 return tour.size();
86 }
87
88 @Override
89 public String toString() {
90 String geneString = "|";
91 for (int i = 0; i < tourSize(); i++) {
92 geneString += getCity(i)+"|";
93 }
94 return geneString;
95 }
96 }
最后,让我们创建模拟退火算法:
1 package sa;
2
3 public class SimulatedAnnealing {
4
5 // Calculate the acceptance probability
6 public static double acceptanceProbability(int energy, int newEnergy, double temperature) {
7 // If the new solution is better, accept it
8 if (newEnergy < energy) {
9 return 1.0;
10 }
11 // If the new solution is worse, calculate an acceptance probability
12 return Math.exp((energy - newEnergy) / temperature);
13 }
14
15 public static void main(String[] args) {
16 // Create and add our cities
17 City city = new City(60, 200);
18 TourManager.addCity(city);
19 City city2 = new City(180, 200);
20 TourManager.addCity(city2);
21 City city3 = new City(80, 180);
22 TourManager.addCity(city3);
23 City city4 = new City(140, 180);
24 TourManager.addCity(city4);
25 City city5 = new City(20, 160);
26 TourManager.addCity(city5);
27 City city6 = new City(100, 160);
28 TourManager.addCity(city6);
29 City city7 = new City(200, 160);
30 TourManager.addCity(city7);
31 City city8 = new City(140, 140);
32 TourManager.addCity(city8);
33 City city9 = new City(40, 120);
34 TourManager.addCity(city9);
35 City city10 = new City(100, 120);
36 TourManager.addCity(city10);
37 City city11 = new City(180, 100);
38 TourManager.addCity(city11);
39 City city12 = new City(60, 80);
40 TourManager.addCity(city12);
41 City city13 = new City(120, 80);
42 TourManager.addCity(city13);
43 City city14 = new City(180, 60);
44 TourManager.addCity(city14);
45 City city15 = new City(20, 40);
46 TourManager.addCity(city15);
47 City city16 = new City(100, 40);
48 TourManager.addCity(city16);
49 City city17 = new City(200, 40);
50 TourManager.addCity(city17);
51 City city18 = new City(20, 20);
52 TourManager.addCity(city18);
53 City city19 = new City(60, 20);
54 TourManager.addCity(city19);
55 City city20 = new City(160, 20);
56 TourManager.addCity(city20);
57
58 // Set initial temp
59 double temp = 10000;
60
61 // Cooling rate
62 double coolingRate = 0.003;
63
64 // Initialize intial solution
65 Tour currentSolution = new Tour();
66 currentSolution.generateIndividual();
67
68 System.out.println("Initial solution distance: " + currentSolution.getDistance());
69
70 // Set as current best
71 Tour best = new Tour(currentSolution.getTour());
72
73 // Loop until system has cooled
74 while (temp > 1) {
75 // Create new neighbour tour
76 Tour newSolution = new Tour(currentSolution.getTour());
77
78 // Get a random positions in the tour
79 int tourPos1 = (int) (newSolution.tourSize() * Math.random());
80 int tourPos2 = (int) (newSolution.tourSize() * Math.random());
81
82 // Get the cities at selected positions in the tour
83 City citySwap1 = newSolution.getCity(tourPos1);
84 City citySwap2 = newSolution.getCity(tourPos2);
85
86 // Swap them
87 newSolution.setCity(tourPos2, citySwap1);
88 newSolution.setCity(tourPos1, citySwap2);
89
90 // Get energy of solutions
91 int currentEnergy = currentSolution.getDistance();
92 int neighbourEnergy = newSolution.getDistance();
93
94 // Decide if we should accept the neighbour
95 if (acceptanceProbability(currentEnergy, neighbourEnergy, temp) > Math.random()) {
96 currentSolution = new Tour(newSolution.getTour());
97 }
98
99 // Keep track of the best solution found
100 if (currentSolution.getDistance() < best.getDistance()) {
101 best = new Tour(currentSolution.getTour());
102 }
103
104 // Cool system
105 temp *= 1-coolingRate;
106 }
107
108 System.out.println("Final solution distance: " + best.getDistance());
109 System.out.println("Tour: " + best);
110 }
111 }
结果如下:
Initial solution distance: 1966
Final solution distance: 911
Tour: |180, 200|200, 160|140, 140|180, 100|180, 60|200, 40|160, 20|120, 80|100, 40|60, 20|20, 20|20, 40|60, 80|100, 120|40, 120|20, 160|60, 200|80, 180|100, 160|140, 180|
在这个例子中,我们能够超过我们初始随机生成路径的一半以上。很大程度上说明,当应用到某些类型的优化问题时,这个相对简单的算法是多么方便。
模拟退火算法的应用很广泛,可以求解NP完全问题,但其参数难以控制,其主要问题有以下三点:
(1) 温度T的初始值设置问题。 温度T的初始值设置是影响模拟退火算法全局搜索性能的重要因素之一、初始温度高,则搜索到全局最优解的可能性大,但因此要花费大量的计算时间;反之,则可节约计算时间,但全局搜索性能可能受到影响。实际应用过程中,初始温度一般需要依据实验结果进行若干次调整。 (2) 退火速度问题。 模拟退火算法的全局搜索性能也与退火速度密切相关。一般来说,同一温度下的“充分”搜索(退火)是相当必要的,但这需要计算时间。实际应用中,要针对具体问题的性质和特征设置合理的退火平衡条件。 (3) 温度管理问题。 温度管理问题也是模拟退火算法难以处理的问题之一。实际应用中,由于必须考虑计算复杂度的切实可行性等问题,常采用如下所示的降温方式: T(t+1)=k×T(t) 式中k为正的略小于1.00的常数,t为降温的次数。
,其中
,输入
,求
的最小值。 本题可以用经典的二分法求解,这种方法比较简单,就不说了。主要来说模拟退火做法。---HDU 2899