优化算法----模拟退火算法（Simulated Annealing,SA）

简介

\quad\quad 模拟退火算法是80年代初发展起来的一种基于Monte Carlo迭代求解策略的随机性寻优算法。其思想最早由Metropolis等人于1953年提出,由Krikpatrick等人于1983年引入组合优化领域，目前已在工程中得到了广泛的应用。其克服了贪心（Greedy）算法等局部寻优算法容易陷入局部最优解的缺陷与对初值的依赖性，其优势在于可以以一定的概率选择当前解邻域中较差的解，可以跳出局部最优，从而实现全局寻优的目的。

\quad\quad 在热力学上，退火(annealing)现象指物体逐渐降温物理过程。温度越低， 物体的能量状态越低，其分子结构组成越有规律，从而使得物体在完全结晶状态时，系统的能量状态达到最低。但如果在某一刻骤然降温（热力学中称为淬炼（quenching）），会导致非能量最低的非完全晶体状态。

与贪心算法的比较

\quad\quad 如下图所示，以目标函数取最小值为例。若初始解定为A点，根据贪心算法，我们可以一步步的减小目标函数值，从而到达B点。此时算法会结束。但使用模拟退火算法，其重要的一个特点就是可以以一定的概率跳出局部最优点，如B点，接受一个在BC之间的解，进而通过逐次的跳变，可以取得D点，实现全局最优。

物理退火过程

\quad\quad 模拟退火算法的思想来源于对固体退火降温过程的模拟。即将固体加温至充分高的，再给其逐渐降温冷却的过程。具体为：
（1）升温。将固体温度升高，从而使固体中粒子的运动状态处于完全无序的状态。
（2）等温。固定在当前的温度，不再升温。对于与周围环境交换热量而温度不变的封闭系统，系统会自发地总是朝着自由能减小的方向进行，从而系统达到在当前温度下平衡的新状态。
（3）逐渐冷却。对新的平衡状态进行冷却处理，使得系统中粒子的热运动减弱并逐渐趋于有序，系统的总能量逐渐降低，从而得到低能量的晶体结构。

算法流程

\quad\quad 在组合优化的实际应用领域，可以将固体内能模拟为最小化的目标函数值 $f$, 将温度 $T$ 模拟为控制参数。可以假设初始将温度升高为一个定值，在该温度下随机产生一个初始解，在当前温度下，算法持续进行“产生新解------计算目标函数差------判断是否接受新解------接受或舍弃” 的迭代过程，从而达到在该温度下的热平衡状态。然后继续减小控制参数 $T$ 的取值，重复上述在一定温度下达到热平衡状态的过程。当控制参数逐渐减小并趋于零时，系统整体趋于平衡状态，从而达到整体最优解。 算法流程如下图

Metropolis准则

Metropolis准则是一种有效的重点抽样法。其思想为：系统从一个当前能量状态变化到新的能量状态。若新的能量状态小于当前能量状态，则以概率1接受新的能量状态（取得局部最优）；若新的能量状态大于当前能量状态（全局搜索最优点），则以一定的概率接受或舍弃新的能力状态。
数学性表示： 设 $E_{Tnew}$ 表示当前温度 $T$ 下新的状态能量，$E_{Tcur}$ 表示当前温度 $T$ 下当前状态能量，则接受新状态的概率 $p$ 可表示为：
$$p=\{\begin{array}{l}1,if{E}_{Tnew}<E_{Tcur}\\ exp(-\frac{E_{Tnew}-E_{Tcur}}{T}),if{E}_{Tnew}\ge E_{Tcur}\end{array}(1)$$从方程（1）可得，当新的能量状态大于当前状态能量时，温度越高，新状态的接受概率越高；当前状态的能量与新状态的能量差越高，新状态的接受概率越高。
注： 温度 $T$ 为Metropolis准则的一个重要控制参数，这个参数的大小控制了随机过程向局部或全局最优解移动的快慢。为避免所谓热力学中淬炼（quenching）的操作，控制参数 $T$ 的值必须缓慢衰减。

重要参数分析

1. 状态产生函数
设计的状态产生函数应该尽可能的保证所产生的候选解遍布全部解空间。一般的状态产生函数由两部分产生，即产生候选解的方式和产生候选解的概率分布。候选解的产生方式由问题的性质决定，通常在当前状态的邻域结构内以一定的概率产生。
2. 初始控制参数温度 $T$
常用的确定初温参数的方法：
（1）随机生成一组可行解，以该解所对应的目标函数值的方差作为初温；
（2）可利用经验或试验确定。
3. 退温函数
常用的退温函数有：
（1）$T_{n+1}=k{T}_n$, $n$ 为当前循环次数，$0<k<1$ 为一个非常接近1的常数。
（2）$T_{n+1}=T_n-\Delta T$, 其中 $\Delta T$ 为自定义温度下降值。
4. 内循环的循环次数（Markov链长度的选取）
（1）每一个温度下迭代相同次数，次数与问题有关。
（2）随着温度的下降，增加迭代次数。
5. 算法终止准则
（1）设置终止温度
（2）设置最大迭代次数
（3）算法搜索到的最优值连续若干步保持不变

Matlab代码演示

【示例】旅行商问题（TSP问题）。假设有一个旅行商人要拜访全国31个省会城市，他需要选择所要走的路径，路径的限制是每个城市只能拜访一次，而且最后要回到原来出发的城市。路径的选择要求是：所选路径的路程为所有路径之中的最小值。
全国31个省会城市的坐标为 [1304 2312；3639 1315；4177 2244；3712 1399；3488 1535；3326 1556；3238 1229；4196 1004；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 2367；3394 2643；3439 3201；2935 3240；3140 3550；2545 2357；2778 2826；2370 2975]。

代码详情

%%################################ 1. 主要函数：模拟退火算法实现 ###################################################
function [obj,order,city] = SA(TIn,TdIn,LIn,city,order,k)
%% Input parameters
%SA Simulated Annealing
%TIn: Current Temperature
%TdIn:Final Temperature
%LIn: Number of initial iterations of the inner loop  
%city: Current City-coordinate
%order: Current city-order
%k:Temperature Attenuation Coefficient
%% Output parameters
%obj: Objective function value(total distance between 31 cities)
%order: New city-order
%city: New City-coordinate
%% Data preparation
iter = 1;%Number of outer loop iterations
T = TIn;
Td = TdIn;
distcur = dist(city);
obj(iter) = distcur;
L = LIn;
%% Main code
% ##################Outer loop############################################
while T>Td
    %% Setting cycle times for inner cycle
    if T<100
        L = 100;
    elseif T<10
        L = 200;
    end
    %% inner cycle
    for i = 1:L
        [tmp_city,order] = new_solution(city,order);
        distnew = dist(tmp_city);
        delta = distnew - distcur;
        if delta < 0
            city = tmp_city;
        else
            p = exp(-delta/T);
            if p>rand
                city = tmp_city;
            end
        end
        distcur = dist(city);
    end
    %% Update the information of the outer loop
    iter = iter + 1;%Oteration number updating
    obj(iter) = distcur;%Objective function value updating
    T = k*T;%Decay temperature
end
end
%%################################ 2. 新解产生函数 ##############################################################
function [tmp_Codnt,order] = new_solution(Codnt,Order)
%% input parameter
%Codnt:a struct defines current city coordinate
%Order:the order of the current solution
%% output parameter
%tmp_Codnt:a struct defines new city coordinate
%order:the new city order
%%
n = size(Codnt,2);
tmp_Codnt = Codnt;
order = Order;
%%
%###################随机产生两个城市编号#######################
p1 = unidrnd(n);
p2 = unidrnd(n);
while p1==p2
    p1 = unidrnd(n);
    p2 = unidrnd(n);
end
%%
%###################对换两个城市位置形成新解#######################
tmp = tmp_Codnt(p1);
tmp_Codnt(p1) = tmp_Codnt(p2);
tmp_Codnt(p2) = tmp;
tmp = order(p1);
order(p1)=order(p2);
order(p2)=tmp;
end
%%######################## 3. 目标函数计算函数 ###############################################################
function distvalue = dist(Codnt)
%% input parameter
%Codnt: a struct defines city coordinate 
%% output parameter
% distvalue:Euclidean distance based on current solution
%%
n = size(Codnt,2);
distvalue = 0;
for i=1:n-1
    distvalue = distvalue + sqrt((Codnt(i).x-Codnt(i+1).x)^2+(Codnt(i).y-Codnt(i+1).y)^2);
end
distvalue = distvalue + sqrt((Codnt(n).x-Codnt(1).x)^2+(Codnt(n).y-Codnt(1).y)^2);
end
%%############################ 4. 图形绘制函数 ################################################################
function Fgplot(city,order)
%% iuput parameters
%city: city coordinate
%order: city order of the current solution
n = size(city,2);
for i=1:n-1
    plot([city(i).x,city(i+1).x],[city(i).y,city(i+1).y],'bo-');
    hold on;
    text(city(i).x+30,city(i).y+10,num2str(order(i)));
end
plot([city(n).x,city(1).x],[city(n).y,city(1).y],'bo-');
text(city(n).x+10,city(n).y+10,num2str(order(n)));
xlabel('城市位置横坐标');
ylabel('城市位置纵坐标');
end
%%##################################### main.m (参数输入) #######################################################
clear all;
close all;
clc;
T = 300;%Initial Temperature
Td = 0.005;%Terminal Temperature
L = 50;%Number of initial iterations of the inner loop 
k = 0.995;%Temperature attenuation coefficient
CityCoordinate = [1304 2312;3639 1315;4177 2244;3712 1399;3488 1535;3326 1556;3238 1229;4196 1004;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 2367;3394 2643;3439 3201;2935 3240;3140 3550;2545 2357;2778 2826;2370 2975];
nCity = size(CityCoordinate,1);
city = struct([]);
for i = 1:nCity
    city(i).x = CityCoordinate(i,1);
    city(i).y = CityCoordinate(i,2);
end
tic
figure(1)
order = 1:nCity;
Fgplot(city,order);
title('初始解路径')
[obj,order,City] = SA(T,Td,L,city,order,k);
figure(2)
m = length(obj);
plot(obj)
xlabel('迭代次数');
ylabel('城市间总距离');
title('总距离演变曲线')
figure(3)
Fgplot(City,order);
title('最优解路径')
disp('1.新的路径顺序为：')
disp(order)
disp('2.初始目标距离为：')
disp(obj(1))
disp('3.最优目标距离值为：')
disp(obj(m))
toc
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147

图形展示

参考文献

[1] http://cighao.com/2015/12/03/introduction-of-SA/
[2] https://www.cnblogs.com/ranjiewen/p/6084052.html
[3] 公众号----交通与优化之模拟退火
[4] S. Kirkpatrick, C. D. Gelatt, Jr., M. P. Vecchi, 1983. Optimization by Simulated Annealing. Science, 220, No.4598.
[5] 包子阳，余继周. 智能优化算法及其MATLAB实例[M]. 电子工业出版社，2016.
[6] python 优化包scikit-opt