簡述

算法設(shè)計(jì)課這周的作業(yè)：
趕緊寫了先，不然搞不完了。

文章目錄

• • 簡述
  • 算法理論部分
  • • 變量簡單分析
    • 從狀態(tài)轉(zhuǎn)移概率到狀態(tài)概率
    • 推導(dǎo)
    • 理解當(dāng)溫度收斂到接近0的時(shí)候，收斂到結(jié)果
    • 理論部分的后記
  • python實(shí)現(xiàn)
  • • python實(shí)現(xiàn)模擬退火解極值
  • C++實(shí)現(xiàn)解函數(shù)值問題
  • TSP問題求解
  • • 代碼詳細(xì)解釋
    • • 定義了兩個(gè)宏，主要是為了加速運(yùn)算
      • 定義全局變量
      • 函數(shù)解釋
      • 給一組數(shù)據(jù)和結(jié)果

算法理論部分

• 用粒子的排列或相應(yīng)的能量表示物體所處的狀態(tài)，在溫度T下，物體(系統(tǒng))所處的狀態(tài)具有一定的隨機(jī)性。主流趨勢(shì)是系統(tǒng)向能量較低的狀態(tài)發(fā)展，但粒子的不規(guī)則熱運(yùn)動(dòng)妨礙系統(tǒng)準(zhǔn)確落入低能狀態(tài)。

簡單來說就是，在溫度T下，

有兩種可能

• $f(x_i)\ge f(x_j)$， 那沒得說有更好的肯定選更好的。
• $f(x_i)<f(x_j)$，這個(gè)時(shí)候我們也有一定的概率選這個(gè)。

概率為

$$P(X_{i\to j})=e^{\frac{f(x_i)?f(x_j)}{KT}}$$

其中K為玻爾茲曼常數(shù)（在這個(gè)算法上，我們經(jīng)常使用數(shù)值1代替這個(gè)）

• $K>0$，一般設(shè)置為$K=1$
• $T>0$，我們?cè)谶f減的過程中，終止的T，一般認(rèn)為是$T=1$

變量簡單分析

不難看出，隨著T下降，這個(gè)狀態(tài)轉(zhuǎn)移的概率是下降的。

• 原因：$f(x_i)?f(x_j)<0$

物理上解釋：

• 因?yàn)橹暗倪@個(gè)轉(zhuǎn)移概率，認(rèn)為是，在高溫情況下，分子可能會(huì)產(chǎn)生較為不穩(wěn)定的隨機(jī)運(yùn)動(dòng)。且溫度越高，這個(gè)分子不規(guī)則運(yùn)動(dòng)的可能是更大的。這是在模擬這個(gè)過程。

• 所以，我們稱這個(gè)算法為，模擬退火算法

從狀態(tài)轉(zhuǎn)移概率到狀態(tài)概率

這個(gè)分析過程，其實(shí)是類似于推理馬爾可夫鏈的過程。

• 之前給給出的$P(X_{i\to j})$其實(shí)是條件概率。

我們假設(shè)第n個(gè)狀態(tài)為$X_i$，那么現(xiàn)在就是考慮下一個(gè)狀態(tài)的概率。

$$P(s_{n+1}=X_j|s_n=X_i)=P(X_{i\to j})=e^{\frac{f(x_i)?f(x_j)}{KT}}$$

我們用$s_n$，表示第n個(gè)狀態(tài)。

用條件概率公式得到

$$P(s_{n+1}=X_j|s_n=X_i)=\frac{P(s_{n+1}=X_j,s_n=X_i)}{P(s_n=X_i)}$$

$$P(s_{n+1}=X_j,s_n=X_i)=P(s_n=X_i)?P(s_{n+1}=X_j|s_n=X_i)$$
$$P(s_{n+1}=X_j)=\sum _{i}P(s_n=X_i)?P(s_{n+1}=X_j|s_n=X_i)$$

而解這個(gè)過程，就用到了以前解馬爾科夫過程的方法，這里我需要假設(shè)一個(gè)均衡態(tài)，在這種狀態(tài)下，經(jīng)過任意次狀態(tài)轉(zhuǎn)移之后，整個(gè)模型的概率保持一致。

得到結(jié)果是

$$P_i(T)=\frac{e^{\frac{?f(x_i)}{KT}}}{Z_T}$$

$Z_T$只是一個(gè)歸一化的因子而已，就是把所有的這樣的分子的數(shù)值求個(gè)和。使得概率和為一而已。

這個(gè)分布稱之為 Boltzmann分布。

推導(dǎo)

$$P_i(T)?P_j(T)=\frac{e^{\frac{?f(i)}{KT}}}{Z_T}(1?e^{?\frac{f(j)?f(i)}{KT}})$$

理解當(dāng)溫度收斂到接近0的時(shí)候，收斂到結(jié)果

其實(shí)只需要理解下面這個(gè)函數(shù)的導(dǎo)數(shù)就可以了

$$e^{?\frac{E(i)}{KT}}$$
關(guān)于T求導(dǎo)。

$$e^{?\frac{E(i)}{KT}}?\frac{E(i)}{k{T}^2}$$

那么，這就是這個(gè)變量關(guān)于T的變化導(dǎo)數(shù)。再考慮關(guān)于函數(shù)值$E(i)$

這個(gè)函數(shù)是關(guān)于$E(i)$的單調(diào)減函數(shù)。$E(i)>0$的情況下。

所以說，函數(shù)值稍微小的數(shù)值所對(duì)應(yīng)的狀態(tài)概率，受到T的減小而導(dǎo)致的減小的幅度，其實(shí)是較為小的。
那么當(dāng)T趨于0的時(shí)候，就可以得到，狀態(tài)概率，會(huì)集中在函數(shù)值較為小的數(shù)值點(diǎn)上（數(shù)值小的點(diǎn)的概率大）

也就是說，當(dāng)模擬退火，溫度越低，越有可能收斂到正確解。而且，這是收斂的。

理論部分的后記

這里，我們證明了，當(dāng)溫度越小，越會(huì)收斂到正確解。這只是在理論上的證明。但是我們都知道，當(dāng)T趨于0，但是特別小的時(shí)候，作為分母時(shí)，這個(gè)精度就會(huì)變得非常低了（計(jì)算機(jī)上是離散的）。

• 所以，為了簡單，我們一般設(shè)置T到1就截止了。

python實(shí)現(xiàn)

先嘗試解決下面的這個(gè)圖

import numpy as np import matplotlib.pyplot as pltf = lambda x: (x - 1) ** 2 + x + 100 * np.sin(x)**2 x = np.linspace(-10., 10, 101) plt.plot(x, f(x)) plt.show()

python實(shí)現(xiàn)模擬退火解極值

• 先在這個(gè)區(qū)間上隨機(jī)生成一個(gè)起始點(diǎn)
• 每次在給定點(diǎn)的一個(gè)鄰近區(qū)域（自己設(shè)置），找到一個(gè)新的點(diǎn)。
• 然后這兩個(gè)點(diǎn)做比較。通過概率模擬來確定是否發(fā)生位置遷移。
• 直到溫度下降到一定的數(shù)值。

import numpy as np import matplotlib.pyplot as pltf = lambda x: (x - 1) ** 2 + x + 100 * np.sin(x) ** 2 x = np.linspace(-10., 10, 101) plt.plot(x, f(x))T = 1000 # 起始溫度 alpha = 0.9 # T_{k+1} = alpha * T_k方式更新溫度 limitedT = 1. # 最小值的T iterTime = 1000 # 每個(gè)溫度下迭代的次數(shù) sigmaX = 3 # 每次的搜索半徑 K = 1 # 系數(shù)K X = 20 * np.random.rand() - 10 Y = f(X) while T > limitedT:for i in range(iterTime):xnew = X + sigmaX * (2 * np.random.rand() - 1)if xnew <= 10. and xnew >= -10.:fnew = f(xnew)if fnew < Y:Y = fnewX = xnewelse:res = Y-fnewp = np.exp(res / (K * T))if np.random.rand() < p:X = xnewY = fnewT *= alpha print(X, Y) plt.plot(X, Y, '*r') plt.show()

基本上準(zhǔn)確~

C++實(shí)現(xiàn)解函數(shù)值問題

對(duì)于上面一個(gè)問題的C++解法
但是C++的精度似乎沒有Python的好。這個(gè)可以適當(dāng)調(diào)節(jié)參數(shù)來獲得更高的精度，反正C++的速度也快很多。

#include <iostream> #include <cmath> #include <ctime> using namespace std; #define FUN(x) ((x-1)*(x-1) + x + 100 * sin(x) * sin(x)) #define RAND() ( (double)rand() / double(RAND_MAX))int main() {srand((unsigned)time(NULL));double T = 1000; // 起始溫度double alpha = 0.99; // T_{ k + 1 } = alpha * T_k方式更新溫度double limitedT = 1; // 最小值的Tint iterTime = 1000; // 每個(gè)溫度下迭代的次數(shù)double sigmaX = 3; // 每次的搜索半徑double K = 1; // 系數(shù)Kdouble X = 20 * RAND() - 10;double Y = FUN(X);double xnew, ynew, rest, p;while (T > limitedT) {for (int i = 0; i < iterTime; ++i) {xnew = X + sigmaX * (2 * RAND() - 1);if (xnew >= -10. && xnew <= 10) {ynew = FUN(xnew);if (ynew <= Y) {X = xnew;Y = ynew;}else {rest = Y-ynew;p = exp(rest / (K*T));if (RAND() < p) {X = xnew;Y = ynew;}}}}T *= alpha;}cout << X << " " << Y << endl;system("pause"); }

TSP問題求解

---

代碼詳細(xì)解釋

---

定義了兩個(gè)宏，主要是為了加速運(yùn)算

• 第一個(gè)是RAND(b,e) 生成在$[b,e)$這個(gè)區(qū)間上的整數(shù)
• 第二個(gè)是RANDFLOAT() 是為了生成(0,1)之間的數(shù)。

---

定義全局變量

double **Mat; int *Path, *tempPath; double Value, tempValue; int N = 0;

• Mat是來存儲(chǔ)距離矩陣的
• Path存儲(chǔ)路徑
• Value存儲(chǔ)路勁所對(duì)應(yīng)的函數(shù)。這里考慮的是具體的數(shù)值
• N表示有多少個(gè)點(diǎn)
• 所有前面加了temp都表示臨時(shí)變量。

---

函數(shù)解釋

• double CalValue(int *p) 函數(shù)用于給輸入的路勁下，求對(duì)應(yīng)的value值。
• void refresh() 將temp的數(shù)組和具體數(shù)值恢復(fù)。為了下一次的考慮
• void change() 覆蓋掉原來的路勁。進(jìn)行存儲(chǔ)。
• void initialPath() 初始化路徑。并算出初始值。

#include <iostream> #include <cmath> #include <ctime> #include <fstream> using namespace std; #define RAND(b, e) (rand() % (e-b) + b) // 左閉右開 #define RANDFLOAT() ((double)rand() / double(RAND_MAX)) // 0-1浮點(diǎn)數(shù)double **Mat; int *Path, *tempPath; double Value, tempValue; int N = 0;double CalValue(int *p) {double t = 0;for (int i = 1; i < N; ++i) {t += Mat[p[i - 1]][p[i]];}t += Mat[p[N - 1]][0];return t; }// 將temp重置為path void refresh() {for (int i = 0; i < N; ++i) {tempPath[i] = Path[i];}tempValue = Value; }// 覆蓋修改Path void change() {for (int i = 0; i < N; ++i) {Path[i] = tempPath[i];}Value = tempValue; }void initialPath() {for (int i = 0; i < N; ++i) {Path[i] = i;}// Path[i]表示路上第i個(gè)點(diǎn)的標(biāo)記為Path[i]// Path[0] = 0int tx,t;for (int i = 1; i < N-1; ++i) {tx = RAND(i, N);if (tx != i) { // swapt = Path[i];Path[i] = Path[tx];Path[tx] = t;}}Value = CalValue(Path); }int main() {srand((unsigned)time(NULL));ifstream cin("data.txt");cin >> N;// initializeMat = new double *[N];for (int i = 0; i < N; ++i)Mat[i] = new double[N];for (int i = 0; i < N; ++i) {for (int j = 0; j < N; ++j) {cin >> Mat[i][j];}}Path = new int[N];tempPath = new int[N];double T = 1000; // 起始溫度double alpha = 0.99; // T_{ k + 1 } = alpha * T_k方式更新溫度double limitedT = 1e-9; // 最小值的Tint iterTime = 1000; // 每個(gè)溫度下迭代的次數(shù)double K = 1; // 系數(shù)Kdouble p = 0;initialPath();int tx, ty, t;while (T >= limitedT) {for (int i = 0; i < iterTime; ++i) {// 任意交換兩點(diǎn)的順序refresh();tx = RAND(1, N);ty = RAND(1, N);if (tx != ty) {t = tempPath[tx];tempPath[tx] = tempPath[ty];tempPath[ty] = t;tempValue = CalValue(tempPath);if (tempValue <= Value) {change();}else {p = exp((Value-tempValue) / (K*T));if (RANDFLOAT() < p) { change(); }}}}T *= alpha;}cout << "Value:" << Value << endl;for (int i = 0; i < N; ++i) {cout << Path[i] << " -> ";}cout << 0 << endl;// releasedelete[]tempPath;delete[] Path;for (int i = 0; i < N; ++i)delete[] Mat[i];delete[] Mat;system("pause"); }

給一組數(shù)據(jù)和結(jié)果

10 0 58 82 89 17 50 26 48 70 19 58 0 74 46 70 2 70 49 87 60 82 74 0 58 76 98 37 97 34 67 89 46 58 0 15 17 28 69 46 79 17 70 76 15 0 98 60 69 97 89 50 2 98 17 98 0 81 14 43 47 26 70 37 28 60 81 0 43 73 56 48 49 97 69 69 14 43 0 39 0 70 87 34 46 97 43 73 39 0 53 19 60 67 79 89 47 56 0 53 0

輸出：

Value:244 0 -> 4 -> 3 -> 6 -> 2 -> 8 -> 5 -> 1 -> 7 -> 9 -> 0

總結(jié)

以上是生活随笔為你收集整理的模拟退火算法理论+Python解决函数极值+C++实现解决TSP问题的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。