模拟退火算法是用来求解最优化问题的算法。比如著名的TSP问题,函数最大值最小值问题等等。接下来将以如下几个方面来详细介绍模拟退火算法。
Contents
1. 模拟退火算法认识
2. 模拟退火算法描述
3. 费马点问题求解
4. 最小包含球问题求解
5. 函数最值问题求解
6. TSP问题求解
1. 模拟退火算法认识
爬山算法也是一个用来求解最优化问题的算法,每次都向着当前上升最快的方向往上爬,但是初始化不同可能
会得到不同的局部最优值,模拟退火算法就可能跳出这种局部最优解的限制。模拟退火算法是模拟热力学系统
中的退火过程。在退火过程中是将目标函数作为能量函数。大致过程如下
初始高温 => 温度缓慢下降=> 终止在低温 (这时能量函数达到最小,目标函数最小)
在热力学中的退火过程大致是变温物体缓慢降温而达到分子之间能量最低的状态。设热力学系统S中有有限个且
离散的n个状态,状态的能量为,在温度下,经过一段时间达到热平衡,这时处于状态的概率为
模拟退火算法也是贪心算法,但是在其过程中引入了随机因素,以一定的概率接受一个比当前解要差的解,并且
这个概率随着时间的推移而逐渐降低。
2. 模拟退火算法描述
模拟退火算法来源于固体退火原理,将固体加温至充分高,再让其徐徐冷却,加温时,固体内部粒子随温升变为无序状,内能增大,而徐徐冷却时粒子渐趋有序,在每个温度都达到平衡态,最后在常温时达到基态,内能减为最小。
模拟退火算法(Simulated Annealing,SA)最早由Kirkpatrick等应用于组合优化领域,它是基于Monte-Carlo迭代求解策略的一种随机寻优算法,其出发点是基于物理中固体物质的退火过程与一般组合优化问题之间的相似性。模拟退火算法从某一较高初温出发,伴随温度参数的不断下降,结合概率突跳特性在解空间中随机寻找目标函数的全局最优解,即在局部最优解能概率性地跳出并最终趋于全局最优。模拟退火算法是一种通用的优化算法,理论上算法具有概率的全局优化性能,目前已在工程中得到了广泛应用,诸如VLSI、生产调度、控制工程、机器学习、神经网络、信号处理等领域。
模拟退火算法是通过赋予搜索过程一种时变且最终趋于零的概率突跳性,从而可有效避免陷入局部极小并最终趋于全局最优的串行结构的优化算法。
模拟退火算法的模型
模拟退火算法可以分解为解空间、目标函数和初始解三部分。
模拟退火的基本思想:
(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步。
模拟退火算法流程图
若,即移动后得到更优的解,那么总是接受改移动。
若,即移动后得到更差的解,则以一定的概率接受该移动,并且这个概率随时间推移
逐渐降低。这个概率表示为
由于是退火过程,所以dE < 0,这个公式说明了温度越高出现一次能量差为dE的降温概率就越大,温度越低,
出现降温的概率就越小,由于dE总是小于0,所以P(dE)取值在0到1之间。伪码如下
3. 费马点问题求解
题目:http://poj.org/problem?id=2420
题意:给n个点,找出一个点,使这个点到其他所有点的距离之和最小,也就是求费马点。
代码:
- #include <iostream>
- #include <string.h>
- #include <stdlib.h>
- #include <stdio.h>
- #include <time.h>
- #include <math.h>
-
- #define N 1005
- #define eps 1e-8 //搜索停止条件阀值
- #define INF 1e99
- #define delta 0.98 //温度下降速度
- #define T 100 //初始温度
-
- using namespace std;
-
- int dx[4] = {0, 0, -1, 1};
- int dy[4] = {-1, 1, 0, 0};
-
- struct Point
- {
- double x, y;
- };
-
- Point p[N];
-
- double dist(Point A, Point B)
- {
- return sqrt((A.x - B.x) * (A.x - B.x) + (A.y - B.y) * (A.y - B.y));
- }
-
- double GetSum(Point p[], int n, Point t)
- {
- double ans = 0;
- while(n--)
- ans += dist(p[n], t);
- return ans;
- }
-
-
- double Search(Point p[], int n)
- {
- Point s = p[0];
- double t = T;
- double ans = INF;
- while(t > eps)
- {
- bool flag = 1;
- while(flag)
- {
- flag = 0;
- for(int i = 0; i < 4; i++)
- {
- Point z;
- z.x = s.x + dx[i] * t;
- z.y = s.y + dy[i] * t;
- double tp = GetSum(p, n, z);
- if(ans > tp)
- {
- ans = tp;
- s = z;
- flag = 1;
- }
- }
- }
- t *= delta;
- }
- return ans;
- }
-
- int main()
- {
- int n;
- while(scanf("%d", &n) != EOF)
- {
- for(int i = 0; i < n; i++)
- scanf("%lf %lf", &p[i].x, &p[i].y);
- printf("%.0lf\n", Search(p, n));
- }
- return 0;
- }
题目:平面上给定n条线段,找出一个点,使这个点到这n条线段的距离和最小。
代码:
- #include <iostream>
- #include <string.h>
- #include <stdlib.h>
- #include <stdio.h>
- #include <time.h>
- #include <math.h>
-
- #define N 1005
- #define eps 1e-8 //搜索停止条件阀值
- #define INF 1e99
- #define delta 0.98 //温度下降速度
- #define T 100 //初始温度
-
- using namespace std;
-
- int dx[4] = {0, 0, -1, 1};
- int dy[4] = {-1, 1, 0, 0};
-
- struct Point
- {
- double x, y;
- };
-
- Point s[N], t[N];
-
- double cross(Point A, Point B, Point C)
- {
- return (B.x - A.x) * (C.y - A.y) - (B.y - A.y) * (C.x - A.x);
- }
-
- double multi(Point A, Point B, Point C)
- {
- return (B.x - A.x) * (C.x - A.x) + (B.y - A.y) * (C.y - A.y);
- }
-
- double dist(Point A, Point B)
- {
- return sqrt((A.x - B.x) * (A.x - B.x) + (A.y - B.y) * (A.y - B.y));
- }
-
- double GetDist(Point A, Point B, Point C)
- {
- if(dist(A, B) < eps) return dist(B, C);
- if(multi(A, B, C) < -eps) return dist(A, C);
- if(multi(B, A, C) < -eps) return dist(B, C);
- return fabs(cross(A, B, C) / dist(A, B));
- }
-
- double GetSum(Point s[], Point t[], int n, Point o)
- {
- double ans = 0;
- while(n--)
- ans += GetDist(s[n], t[n], o);
- return ans;
- }
-
- double Search(Point s[], Point t[], int n, Point &o)
- {
- o = s[0];
- double tem = T;
- double ans = INF;
- while(tem > eps)
- {
- bool flag = 1;
- while(flag)
- {
- flag = 0;
- for(int i = 0; i < 4; i++)
- {
- Point z;
- z.x = o.x + dx[i] * tem;
- z.y = o.y + dy[i] * tem;
- double tp = GetSum(s, t, n, z);
- if(ans > tp)
- {
- ans = tp;
- o = z;
- flag = 1;
- }
- }
- }
- tem *= delta;
- }
- return ans;
- }
-
- int main()
- {
- int n;
- while(scanf("%d", &n) != EOF)
- {
- for(int i = 0; i < n; i++)
- scanf("%lf %lf %lf %lf", &s[i].x, &s[i].y, &t[i].x, &t[i].y);
- Point o;
- double ans = Search(s, t, n, o);
- printf("%.1lf %.1lf %.1lf\n", o.x, o.y, ans);
- }
- return 0;
- }
4. 最小包含球问题求解
题目:http://poj.org/problem?id=2069
题意:给定三维空间的n点,找出一个半径最小的球把这些点全部包围住。
代码:
- #include <iostream>
- #include <string.h>
- #include <stdio.h>
- #include <math.h>
-
- #define N 55
- #define eps 1e-7
- #define T 100
- #define delta 0.98
- #define INF 1e99
-
- using namespace std;
-
- struct Point
- {
- double x, y, z;
- };
-
- Point p[N];
-
- double dist(Point A, Point B)
- {
- return sqrt((A.x - B.x) * (A.x - B.x) + (A.y - B.y) * (A.y - B.y) + (A.z - B.z) * (A.z - B.z));
- }
-
- double Search(Point p[], int n)
- {
- Point s = p[0];
- double t = T;
- double ans = INF;
- while(t > eps)
- {
- int k = 0;
- for(int i = 0; i < n; i++)
- if(dist(s, p[i]) > dist(s, p[k]))
- k = i;
- double d = dist(s, p[k]);
- ans = min(ans, d);
- s.x += (p[k].x - s.x) / d * t;
- s.y += (p[k].y - s.y) / d * t;
- s.z += (p[k].z - s.z) / d * t;
- t *= delta;
- }
- return ans;
- }
-
- int main()
- {
- int n;
- while(cin >> n && n)
- {
- for(int i = 0; i < n; i++)
- cin >> p[i].x >> p[i].y >> p[i].z;
- double ans = Search(p, n);
- printf("%.5lf\n", ans);
- }
- return 0;
- }
5. 函数最值问题求解
题目:http://acm.hdu.edu.cn/showproblem.php?pid=2899
题意:给出方程,其中,输入,求的最小值。
分析:本题可以用经典的二分法求解,这种方法比较简单,就不说了。主要来说模拟退火做法。
代码:
- #include <iostream>
- #include <string.h>
- #include <stdlib.h>
- #include <algorithm>
- #include <stdio.h>
- #include <time.h>
- #include <math.h>
-
- #define ITERS 10
- #define T 100
- #define eps 1e-8
- #define delta 0.98
- #define INF 1e99
-
- using namespace std;
-
- double x[ITERS];
-
- int Judge(double x, double y)
- {
- if(fabs(x - y) < eps) return 0;
- return x < y ? -1 : 1;
- }
-
- double Random()
- {
- double x = rand() * 1.0 / RAND_MAX;
- if(rand() & 1) x *= -1;
- return x;
- }
-
- double F(double x, double y)
- {
- return 6 * x * x * x * x * x * x * x + 8 * x * x * x * x * x * x + 7 * x * x * x + 5 * x * x - y * x;
- }
-
- void Init()
- {
- for(int i = 0; i < ITERS; i++)
- x[i] = fabs(Random()) * 100;
- }
-
- double SA(double y)
- {
- double ans = INF;
- double t = T;
- while(t > eps)
- {
- for(int i = 0; i < ITERS; i++)
- {
- double tmp = F(x[i], y);
- for(int j = 0; j < ITERS; j++)
- {
- double _x = x[i] + Random() * t;
- if(Judge(_x, 0) >= 0 && Judge(_x, 100) <= 0)
- {
- double f = F(_x, y);
- if(tmp > f)
- x[i] = _x;
- }
- }
- }
- t *= delta;
- }
- for(int i = 0; i < ITERS; i++)
- ans = min(ans, F(x[i], y));
- return ans;
- }
-
- int main()
- {
- int t;
- scanf("%d", &t);
- while(t--)
- {
- double y;
- scanf("%lf", &y);
- srand(time(NULL));
- Init();
- printf("%.4lf\n", SA(y));
- }
- return 0;
- }
6. TSP问题求解
TSP问题是一个NP问题,但是可以求近似解,通过模拟退火算法实现,代码如下
代码:
- #include <iostream>
- #include <string.h>
- #include <stdlib.h>
- #include <algorithm>
- #include <stdio.h>
- #include <time.h>
- #include <math.h>
-
- #define N 30 //城市数量
- #define T 3000 //初始温度
- #define EPS 1e-8 //终止温度
- #define DELTA 0.98 //温度衰减率
- #define LIMIT 10000 //概率选择上限
- #define OLOOP 1000 //外循环次数
- #define ILOOP 15000 //内循环次数
-
- using namespace std;
-
-
- struct Path
- {
- int citys[N];
- double len;
- };
-
-
- struct Point
- {
- double x, y;
- };
-
- Path path;
- Point p[N];
- double w[N][N];
- 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(Point 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]);
- }
-
- void Input(Point p[], int &n)
- {
- scanf("%d", &n);
- for(int i = 0; i < n; i++)
- scanf("%lf %lf", &p[i].x, &p[i].y);
- }
-
- void Init(int n)
- {
- nCase = 0;
- path.len = 0;
- for(int i = 0; i < n; i++)
- {
- path.citys[i] = i;
- if(i != n - 1)
- {
- printf("%d--->", i);
- path.len += w[i][i + 1];
- }
- else
- printf("%d\n", i);
- }
- printf("\nInit path length is : %.3lf\n", path.len);
- }
-
- void Print(Path t, int n)
- {
- printf("Path is : ");
- 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);
- }
-
- Path GetNext(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(time(NULL));
- Path curPath = path;
- Path newPath = path;
- int P_L = 0;
- int P_F = 0;
- while(1)
- {
- for(int i = 0; i < ILOOP; i++)
- {
- newPath = GetNext(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)
- curPath = newPath;
- P_L++;
- }
- if(P_L > LIMIT)
- {
- P_F++;
- break;
- }
- }
- if(curPath.len < newPath.len)
- path = curPath;
- if(P_F > OLOOP || t < EPS)
- break;
- t *= DELTA;
- }
- }
-
- int main()
- {
- freopen("TSP.data", "r", stdin);
- int n;
- Input(p, n);
- GetDist(p, n);
- Init(n);
- SA(n);
- Print(path, n);
- printf("Total test times is : %d\n", nCase);
- return 0;
- }
数据:TSP.data
- 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