智能优化算法（四）：模拟退火算法

1.算法简介

\quad 模拟退火算法（Simulated Annealing，SA）的思想最早由N. Metropolis等人在1953年提出，1983年，S. Kirkpatrick等人成功地将退火思想引入到组合优化领域，该算法基于Monte-Carlo迭代求解策略，是一种随机寻优算法，其核心思想是物理中固体物质的退火过程与组合优化问题之间的相似性。
\quad 模拟退火算法是一种基于概率的随机优化算法，它尽管也是一种贪心算法，但其搜索过程引入了随机因素，使得算法在一定程度上能够接受较差的解，从而有可能跳出局部最优解，最终达到全局最优解。模拟退火算法通过以一定概率接受一个比当前解更差的解，避免陷入局部最优解。
\quad 在退火算法中，算法在找到局部最优解E后，可能会以一定的概率接受到F的移动，经过几次这样的非最优移动后，可能会到达G点，并进一步上升到H点，从而跳出局部最大值E。模拟退火算法来源于固体退火的原理，固体退火过程是将固体加热至高温，然后逐渐冷却，在加热过程中，固体内部的粒子变得无序，内能增大；在逐渐冷却的过程中，粒子趋于有序，并在每个温度下达到平衡态，最终在常温下达到基态，内能最小。

2.算法推导

2.1.Bolzman方程

\quad 设热力学系统S中有n个状态（有限且离散的），其中状态i的能量为$E_i$ ，在温度$T_k$下，经一段时间达到热平衡，此时处在状态i的概率为：
$$P_i\left(T_k\right)=C_k\exp \left(\begin{array}{c}-E_i\\ T_k\end{array}\right)$$
\quad 那么对于系数$C_k$我们应该如何确定？确定方法如下所示：
$$\begin{array}{rl}& \sum _{j=1}^n{P}_j(T_k)=1\Rightarrow \sum _{j=1}^n{C}_k\cdot \exp \left(\begin{array}{c}-E_j\\ T_k\end{array}\right)=1\\ & \Rightarrow C_k\cdot \sum _{j=1}^n\exp \left(\begin{array}{c}-E_j\\ T_k\end{array}\right)=1⟹C_k=\frac{1}{{\sum }_{j=1}^n\exp \left(\begin{array}{c}-E_j\\ T_k\end{array}\right)}\end{array}$$
\quad 通过求$C_k$， 获得Bolzman方程：
$$P_i\left(T_k\right)=\frac{\exp \left(\begin{array}{c}-E_i\\ T_k\end{array}\right)}{{\sum }_{j=1}^n\exp \left(\begin{array}{c}-E_j\\ T_k\end{array}\right)}$$

2.2.Metropolis准则

\quad 本质思想为以概率接受新状态：
\quad 若在温度$T_k$ ， 当前状态$i$ →新状态$j$，若$E_j$<$E_i$， 则接受$j$为当前状态.
\quad 否则， 若概率 $p=exp[-(Ej-Ei)/Tk]$大于[0,1)区间的随机数， 则仍接受状态 $j$为当前状态； 若不成立则保留状态$i$为当前状态.

2.3.算法要素

\quad $\bullet$初始温度$T_0$
\quad ①均匀抽样一组状态，以各状态目标值的方差为初温。
\quad ② 随机产生一组状态，确定两两状态间的最大目标值差，根据差值，利用一定的函数确定初温。
\quad ③ 利用经验公式。
实验表明：初温越大，获得高质量解的机率越大，但花费较多的计算时间！！
\quad $\bullet$降温函数
\quad ①$T_{k+1}=T_k\cdot r$ \quad $r\in (\mathrm{0.950.99})$
\quad ②$T_{k+1}=T_k-\Delta T$
\quad $\bullet$热平衡条件
\quad Metropolis抽样稳定准则:
\quad ①设置内循环迭代次数
\quad ②检验目标函数的均值是否稳定
\quad ③连续若干步的目标值变化较小
\quad 停止准则：
\quad ①设置外循环迭代次数
\quad ②设置终止温度的阈值
\quad ③算法搜索到的最优值连续若干步保持不变

2.4.算法流程

3.算法实现

3.1.算法MTALAB实现

%%%%%%%%%%%%%%%%%%%%%%模拟退火算法解决 TSP 问题%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%初始化%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all; %清除所有变量
close all; %清图
clc; %清屏
tic
%成绩坐标点
C=[1304 2312;3639 1315;4177 2244;3712 1399;3488 1535;3326 1556;...
3238 1229;4196 1044;4312 790;4386 570;3007 1970;2562 1756;...
2788 1491;2381 1676;1332 695;3715 1678;3918 2179;4061 2370;...
3780 2212;3676 2578;4029 2838;4263 2931;3429 1908;3507 2376;...
3394 2643;3439 3201;2935 3240;3140 3550;2545 2357;2778 2826;...
2370 2975]; %31 个省会城市坐标
n=size(C,1); %TSP 问题的规模,即城市数目
T=100*n; %初始温度
L=100; %马可夫链长度
K=0.8; %衰减参数
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%城市坐标结构体%%%%%%%%%%%%%%%%%%%%%%%%%%
city=struct([]);
for i=1:n
    city(i).x=C(i,1);
    city(i).y=C(i,2);
end
l=1; %统计迭代次数
len(l)=func3(city,n); %每次迭代后的路线长度
figure(1);
while T>0.001 %停止迭代温度
%%%%%%%%%%%%%%%%多次迭代扰动，温度降低之前多次实验%%%%%%%%%%%%%%%
    for i=1:L
        %%%%%%%%%%%%%%%%%%%计算原路线总距离%%%%%%%%%%%%%%%%%%%%%%%%%
        len1=func3(city,n);
        %%%%%%%%%%%%%%%%%%%%%%%%%产生随机扰动%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%随机置换两个不同的城市的坐标%%%%%%%%%%%%%%%%%
        p1=floor(1+n*rand());
        p2=floor(1+n*rand());
        while p1==p2
            p1=floor(1+n*rand());
            p2=floor(1+n*rand());
        end
        tmp_city=city;
        tmp=tmp_city(p1);
        tmp_city(p1)=tmp_city(p2);
        tmp_city(p2)=tmp;
        %%%%%%%%%%%%%%%%%%%%%%%%计算新路线总距离%%%%%%%%%%%%%%%%%%%%
        len2=func3(tmp_city,n);
        %%%%%%%%%%%%%%%%%%新老距离的差值，相当于能量%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%新路线好于旧路线，用新路线代替旧路线%%%%%%%%%%%%%%
        delta_e=len2-len1;
        if delta_e<0
            city=tmp_city;
        else
        %%%%%%%%%%%%%%%%%%以概率选择是否接受新解%%%%%%%%%%%%%%%%%
        if exp(-delta_e/T)>rand()
            city=tmp_city;
        end
        end
    end
    l=l+1;
    %%%%%%%%%%%%%%%%%%%%%%%%%计算新路线距离%%%%%%%%%%%%%%%%%%%%%%%%%%
    len(l)=func3(city,n);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%温度不断下降%%%%%%%%%%%%%%%%%%%%%%%%%%
    T=T*K;
    for i=1:n-1
        plot([city(i).x,city(i+1).x],[city(i).y,city(i+1).y],'bo-');
        hold on;
    end
    plot([city(n).x,city(1).x],[city(n).y,city(1).y],'ro-');
    
    title(['优化最短距离(K=0.8):',num2str(len(l))]);
    hold off;
    pause(0.005);
end
%print(gcf,'C:\Users\Zeng Zhong Yan\Desktop\SA_result','-dpng','-r600')

figure(2);
plot(len)
xlabel('迭代次数')
ylabel('目标函数值')
title('适应度进化曲线')
print(gcf,'C:\Users\Zeng Zhong Yan\Desktop\SA','-dpng','-r600')

3.2.算法PYTHON实现

import numpy as np
import matplotlib.pyplot as plt
import math
import random
import time

# Function to calculate the total distance of the tour
def calculate_total_distance(city, n):
    total_distance = 0
    for i in range(n - 1):
        total_distance += np.linalg.norm(np.array([city[i]['x'], city[i]['y']]) - np.array([city[i + 1]['x'], city[i + 1]['y']]))
    total_distance += np.linalg.norm(np.array([city[n - 1]['x'], city[n - 1]['y']]) - np.array([city[0]['x'], city[0]['y']]))
    return total_distance

# Initialize coordinates
C = np.array([
    [1304, 2312], [3639, 1315], [4177, 2244], [3712, 1399], [3488, 1535], [3326, 1556],
    [3238, 1229], [4196, 1044], [4312, 790], [4386, 570], [3007, 1970], [2562, 1756],
    [2788, 1491], [2381, 1676], [1332, 695], [3715, 1678], [3918, 2179], [4061, 2370],
    [3780, 2212], [3676, 2578], [4029, 2838], [4263, 2931], [3429, 1908], [3507, 2376],
    [3394, 2643], [3439, 3201], [2935, 3240], [3140, 3550], [2545, 2357], [2778, 2826],
    [2370, 2975]
])

n = len(C)  # Number of cities
T = 100 * n  # Initial temperature
L = 100  # Length of the Markov chain
K = 0.8  # Decay parameter

# City coordinates as a list of dictionaries
city = [{'x': C[i, 0], 'y': C[i, 1]} for i in range(n)]

# Initialize variables for iteration
iteration = 1
len_list = [calculate_total_distance(city, n)]

plt.figure(1)
start_time = time.time()

while T > 0.001:  # Stop iteration when temperature is low enough
    for _ in range(L):
        len1 = calculate_total_distance(city, n)
        p1, p2 = random.sample(range(n), 2)
        tmp_city = city.copy()
        tmp_city[p1], tmp_city[p2] = tmp_city[p2], tmp_city[p1]
        len2 = calculate_total_distance(tmp_city, n)
        delta_e = len2 - len1
        if delta_e < 0 or math.exp(-delta_e / T) > random.random():
            city = tmp_city

    iteration += 1
    len_list.append(calculate_total_distance(city, n))
    T *= K

    for i in range(n - 1):
        plt.plot([city[i]['x'], city[i + 1]['x']], [city[i]['y'], city[i + 1]['y']], 'bo-')
    plt.plot([city[n - 1]['x'], city[0]['x']], [city[n - 1]['y'], city[0]['y']], 'ro-')

    plt.title(f'Optimized Shortest Distance (K=0.8): {len_list[-1]}')
    plt.pause(0.005)
    plt.clf()

print(f"Execution Time: {time.time() - start_time:.2f} seconds")

# Save the final result plot
plt.figure(1)
for i in range(n - 1):
    plt.plot([city[i]['x'], city[i + 1]['x']], [city[i]['y'], city[i + 1]['y']], 'bo-')
plt.plot([city[n - 1]['x'], city[0]['x']], [city[n - 1]['y'], city[0]['y']], 'ro-')
plt.title(f'Optimized Shortest Distance (K=0.8): {len_list[-1]}')
plt.savefig('SA_result.png', dpi=600)

# Plot the evolution of the objective function value
plt.figure(2)
plt.plot(len_list)
plt.xlabel('Iteration')
plt.ylabel('Objective Function Value')
plt.title('Fitness Evolution Curve')
plt.savefig('SA.png', dpi=600)
plt.show()

4.参考文献

[1]https://blog.csdn.net/bingokunkun/article/details/118583729
[2]https://www.cnblogs.com/dion-90/articles/8469351.html
[3]https://chatgpt.com/