現代優化演算法探究 模擬退火演算法
相較於完全貪心的單點爬山演算法，模擬退火用概率接受新解的方式，對貪心法容易陷入區域性最優解的缺陷進行了改進，是一種常用的在大的搜尋空間中逼近全域性最優解的元啟發式方法。這裡先大致描述演算法。然後先用一個簡單的求函式最小值例子，解釋我寫的c++模擬退火演算法模板的基本使用，之後用其解決15節點的TSP問題，並與動態規劃得到的全域性最優解進行比較，解釋模擬退火演算法解決一般問題的方法和效果。
演算法描述：
隨機從一個合法解s開始，從s擴展出足夠的狀態，對於新的狀態，若新狀態優於原始狀態，則接受新解，否則，以一定的概率接受新解（這保證了演算法有機會跳出區域性最優解），而該概率隨演算法執行時間下降，類似於粒子退火的過程（這也是合理的，因為隨著演算法的進行，當前解就逐步逼近於最優解，那麼接受新解的概率自然要下降）。
其實模擬退火就是不停的基於當前最優解，試圖進行改良，並且以概率接受不是那麼好的解來避免陷入區域性最優解，並不複雜。
演算法虛擬碼：
模擬退火模板題： 求函式(x - 2)*(x + 3)*(x + 8)*(x - 9)的最小值
分析：
直接使用模板即可。在定義問題時，關鍵要素有如下幾個：
1. 狀態的定義
2. 對應於一個狀態的效用值的計算
3. 產生新解的方式
然後還要注意的是，引數的調整，比如時間到溫度如何轉化，迭代次數等等。引數調整也是所有這些非經典演算法都要考慮的問題。
程式碼：
#include<iostream> #include<algorithm> #include<map> #include<set> #include<queue> #include<sstream> #include<cmath> #include<iterator> #include<bitset> #include<stdio.h> #include<unordered_set> #include<ctime> using namespace std; #define _for(i,a,b) for(int i=(a);i<(b);++i) #define _rep(i,a,b) for(int i=(a);i<=(b);++i) typedef long long LL; const int INF = 1 << 30; const int maxn = 100005; const double eps = 1e-6; template< typename T = double> class State { //狀態的定義 public: T x, val; T getval(State<T> & s) { //效用值得計算 auto x = s.val; return (x - 2)*(x + 3)*(x + 8)*(x - 9); } T getval(T & x) { return (x - 2)*(x + 3)*(x + 8)*(x - 9); } State(T x) :x(x), val(getval(x)) {} State() {}; State(const State<T> & rhs) :x(rhs.x), val(rhs.val) {} }; template<typename T = double> class SA { //Simulated Annealing private: static double schedule(int t,double tem) { return tem*temrate; } static State<T> randomNe(State<T> & s) { //產生新解的方式 double y = s.x + 2 * ((double)(rand() % 1000)) / 1000 - 1; return State<T>(y); } static const int maxtry = 100; public: static const double temrate; static const double inittem; State<T> minans(State<T> s) { State<T> current = s, ans = s; double tem = SA::inittem; for (int t = 1;; t++) { tem=schedule(t,tem); if (tem <= eps)return ans; for (int i = 0; i<maxtry; ++i) { State<T> ne = randomNe(current); double detaE = current.val - ne.val; if (detaE > 0)ans = current = ne; else if (rand()<exp(detaE / tem) * 0x7fff) current = ne; } //cout << current.val << endl; } } }; const double SA<double>::inittem = 1000; const double SA<double>::temrate = 0.99; int main() { // freopen("C:\\Users\\admin\\Desktop\\in.txt", "r", stdin); // freopen("C:\\Users\\admin\\Desktop\\out.txt", "w", stdout); srand((unsigned)time(NULL)); double s0 = rand() % 100 / 10; SA<double> solver; State<double> ans = solver.minans(State<double>(s0)); cout << ans.x << " " << ans.val << endl; return 0; }
最後結果可以大致穩定在x=6.5左右，某次執行結果為x=6.492 y=-1549.72.
在求解過程中，當前解的演變過程如下：
可以看到隨著迭代的次數增加，當前解逐步趨向於最小值，過程類似於粒子退火，初始時能量大，無序性大，最後趨於穩定。
該函式的圖形如下：
可以看出得出的結果接近於全域性最優值。
15節點的TSP問題：
問題描述：
有15個城市, 給出每個城市的座標，每兩個城市之間相互有路徑可走，需要消耗一定的費用。問一個商人，想從某個起點出發經過每個城市一次且僅僅一次最後回到起點所需的最小費用是多少。本題直接把兩點的距離作為費用。
資料：
53.7121 15.3046 51.1758 0.0322 46.3253 28.2753 30.3313 6.9348
56.5432 21.4188 10.8198 16.2529 22.7891 23.1045 10.1584 12.4819
20.1050 15.4562 1.9451 0.2057 26.4951 22.1221 31.4847 8.9640
26.2418 18.1760 44.0356 13.5401 28.9836 25.9879
從上到下，從左往右，每兩個數代表一個城市的x 和y 座標。
分析：
仍然使用模擬退火演算法。同樣需要考慮三個問題：狀態的定義，效用值得計算，產生新解的方式。 具體看程式碼。
程式碼：
#include<iostream>
#include<cstdio>
#include<string>
#include<cstring>
#include<vector>
#include<stack>
#include<algorithm>
#include<map>
#include<set>
#include<queue>
#include<sstream>
#include<cmath>
#include<iterator>
#include<bitset>
#include<stdio.h>
#include<unordered_set>
#include<ctime>
#include<assert.h>
using namespace std;
#define _for(i,a,b) for(int i=(a);i<(b);++i)
#define _rep(i,a,b) for(int i=(a);i<=(b);++i)
typedef long long LL;
const int INF = 1 << 30;
const int maxn = 15;
struct Node {
double x, y;
};
Node a[maxn];
int order[maxn];
double dist[maxn][maxn];
template< typename T = double>
class State { //狀態定義
public:
T val;
int path[maxn];
T getval() { //效用值定義
double x = 0;
for (int i = 1; i <= maxn; ++i) x += dist[path[i%maxn]][path[i - 1]];
return val=x;
}
//double getval(double & x) { return (x - 2)*(x + 3)*(x + 8)*(x - 9); }
State() {
memcpy(path, order, sizeof(path));
random_shuffle(path, path + maxn);
val = getval();
};
State(const State<T> & rhs) :val(rhs.val) { memcpy(path, rhs.path, sizeof(path)); }
};
template<typename T = double>
class SA { //Simulated Annealing
private:
static double schedule(int t) {
return 100- 0.05*t;
}
static State<T> randomNe(State<T> & s) { //產生新解
State<T> ne = s;
int head = max(1, rand() % maxn), tail = max(1, rand() % maxn);
if (head > tail)swap(head, tail);
reverse(ne.path+head, ne.path + tail);
ne.val=ne.getval();
return ne;
}
static const int maxtry = 1000;
public:
State<T> minans(State<T> & s) {
State<T> current = s, ans = s;
for (int t = 1;; t++) {
double tem = schedule(t);
if (tem <= 0)return ans;
for (int i = 0; i<maxtry; ++i) {
State<T> ne = randomNe(current);
double detaE = current.val - ne.val;
if (detaE > 0)ans = current = ne;
else if (rand()<exp(detaE / tem) * 0x7fff) current = ne;
}
//printf("%lf\n", ans.val);
}
}
};
int main()
{
//freopen("C:\\Users\\admin\\Desktop\\in.txt", "r", stdin);
//freopen("C:\\Users\\admin\\Desktop\\out.txt", "w", stdout);
srand((unsigned)time(NULL));
for (int i = 0; i < maxn; ++i) {
scanf("%lf%lf", &a[i].x, &a[i].y);
}
for (int i = 0; i < maxn; ++i) order[i] = i;
for (int i = 0; i < maxn; ++i)
for (int j = 0; j < maxn; ++j)
dist[i][j] = sqrt((a[i].x - a[j].x)*(a[i].x - a[j].x) + (a[i].y - a[j].y)*(a[i].y - a[j].y));
SA<double> solver;
State<double> s;
State<double> ans = solver.minans(s);
printf("%lf\n", ans.val);
for (int i = 0; i < maxn; ++i) {
printf("%d ", ans.path[i]);
}
return 0;
}結果能穩定在170上下，某次執行結果如下：
161.241317 （最小總費用）
2 14 10 6 12 8 5 7 9 3 11 13 1 0 4 （路徑，城市下標從0開始）
解的演變過程如下：
動態規劃驗證：
由於這裡的重點只在模擬退火演算法，這裡只給出動態規劃的狀態轉移方程，不做更多的解釋了。
dp(S,i) 表示當前處在城市i上，已經走過bitset集合S中的城市，所耗費的最小費用。
所以狀態轉移方程如下：
其中 dist(i,j) 表示城市i和j之間的距離(或是說費用），初始條件dp(1,0)=0，時間複雜度大致為
程式碼：
#include<iostream>
#include<cstdio>
#include<string>
#include<cstring>
#include<vector>
#include<stack>
#include<algorithm>
#include<map>
#include<set>
#include<queue>
#include<sstream>
#include<cmath>
#include<iterator>
#include<bitset>
#include<stdio.h>
#include<unordered_set>
#include<ctime>
#include<assert.h>
using namespace std;
#define _for(i,a,b) for(int i=(a);i<(b);++i)
#define _rep(i,a,b) for(int i=(a);i<=(b);++i)
typedef long long LL;
const int INF = 1 << 30;
const int maxn = 15;
struct Node {
double x, y;
};
Node a[maxn];
int order[maxn];
double dist[maxn][maxn];
double d[1 << maxn][maxn];
int p[1 << maxn][maxn];
void print(int S, int u, double & tot) {
if (u == 0) {
cout << u << " ";
//cout << endl << tot;
return;
}
tot += dist[u][p[(S)][u]];
cout << u << " ";
print(S ^ (1 << u), p[S][u], tot);
}
int main()
{
//freopen("C:\\Users\\admin\\Desktop\\in.txt", "r", stdin);
//freopen("C:\\Users\\admin\\Desktop\\out.txt", "w", stdout);
srand((unsigned)time(NULL));
for (int i = 0; i < maxn; ++i) {
scanf("%lf%lf", &a[i].x, &a[i].y);
}
for (int i = 0; i < maxn; ++i) order[i] = i;
for (int i = 0; i < maxn; ++i)
for (int j = 0; j < maxn; ++j)
dist[i][j] = sqrt((a[i].x - a[j].x)*(a[i].x - a[j].x) + (a[i].y - a[j].y)*(a[i].y - a[j].y));
int m = 1 << maxn;
for (int i = 0; i < m; ++i)
for (int j = 0; j < maxn; ++j)
d[i][j] = INF;
d[1][0] = 0;
for (int S = 1; S < m; ++S) {
if (!S & 1) continue;
for (int i = 0; i < maxn; ++i) {
if (S&(1 << i)) {
for (int j = 0; j < maxn; ++j)if (S&(1 << j)) {
if (d[S][i] > d[S ^ (1 << i)][j] + dist[i][j]) {
d[S][i] = min(d[S][i], d[S ^ (1 << i)][j] + dist[i][j]);
p[S][i] = j;
}
}
}
}
}
double ans = INF;
int ind;
for (int i = 0; i < maxn; ++i) {
if (ans > d[m - 1][i] + dist[i][0]) {
ans = min(ans, d[m - 1][i] + dist[i][0]);
ind = i;
}
}
cout << ans << endl;
double tot = 0;
print((1 << maxn) - 1, ind, tot);
return 0;
}所得結果如下：
161.241
1 13 11 3 9 7 5 8 12 6 10 14 2 4 0
與以上模擬退火演算法結果一致。
現在將maxn調至50，模擬退火演算法很快給出結果：
300.458339
4 37 33 2 30 27 43 0 23 48 1 32 31 13 11 17 3 35 36 8 25 12 34 28 15 45 18 16 14 10 6 40 41 21 44 38 47 5 22 7 49 24 46 9 29 42 39 26 19 20
而因動態規劃的時空複雜度過高，此時動態規劃無法求解。
結論：
模擬退火演算法能有效求解一般優化問題，並且能通過調整引數在較短的時間內得到很好的結果。是一個在找不到多項式複雜度演算法的情況下求解一般較大規模優化問題的較好選擇。