最近学习了模拟退火智能算法,顺便学习了一下matlab,顺带一提matlba真香,对矩阵的操作会比python的numpy操作要更加方便,这里我们是以TSP问题为例子,因为比较好理解。
模拟退火介绍
模拟退火总的来说还是一种优化算法,他模拟的是淬火冶炼的一个过程,通过升温增强分子的热运动,然后再慢慢降温,使其达到稳定的状态。
初始解
通常是以一个随机解作为初始解. 并保证理论上能够生成解空间中任意的解,也可以是一个经挑选过的较好的解,初始解不宜“太好”, 否则很难从这个解的邻域跳出,针对问题去分析。
扰动邻解
邻解生成函数应尽可能保证产生的侯选解能够遍布解空间,邻域应尽可能的小,能够在少量循环步中允分探测.,但每次的改变不应该引起太大的变化。
初始温度
- 初始温度应该设置的尽可能的高, 以确保最终解不受初始解
影响. 但过高又会增加计算时间. - 均匀抽样一组状态,以各状态目标值的方差为初温.
- 如果能确定邻解间目标函数(COST 函数) 的最大差值, 就可以确定出初始温度 T 0 T_0 T0 , 以使初始接受概率 P = e − ∣ Δ C ∣ m a x / T P=e^{-|\Delta C|_{max}/T} P=e−∣ΔC∣max/T足够大. ∣ Δ C ∣ m a x |\Delta C|_{max} ∣ΔC∣max可由随机产生一组状态的最大目标值差来替代.
- 在正式开始退火算法前, 可进行一个升温过程确定初始温度:逐渐增加温度, 直到所有的尝试尝试运动都被接受, 将此时的温度设置为初始温度.
- 由经验给出, 或通过尝试找到较好的初始温度.
等温步数 - 等温步数也称Metropolis 抽样稳定准则, 用于决定在各温度下产生候选解的数目. 通常取决于解空间和邻域的大小.
- 等温过程是为了让系统达到平衡, 因此可通过检验目标函数的均值是否稳定(或连续若干步的目标值变化较小) 来确定等温步数.
- 等温步数受温度的影响. 高温时, 等温步数可以较小, 温度较小时, 等温步数要大. 随着温度的降低, 增加等温步数.
- 有时为了考虑方便, 也可直接按一定的步数抽样。
**Metropolis法则 **
Metropolis法则是SA接受新解的概率。
P ( x = > x ′ ) = { 1 , f ( x ′ ) < f ( x ) e − f ( x ′ ) − f ( x ) T , f ( x ′ ) > f ( x ) P(x=>x')= \begin{cases} 1 &\text{$,f(x')<f(x)$} \\ e^{-\frac{f(x')-f(x)}{T}} & \text{$,f(x')>f(x)$}\\ \end{cases} P(x=>x′)={1e−Tf(x′)−f(x),f(x′)<f(x),f(x′)>f(x)
x x x是表示当前解, x ′ x' x′是新解,其实这也是模拟退火区别于贪心的一点,我们在更新新解的时候对于不满足条件的情况,我们也有一定的概率来进行选取,从而可以使得退火模拟可以跳出局部最优解。
降温公式
经典模拟退火算法的降温方式:
T ( t ) = T 0 l o g ( 1 + t ) T(t)=\frac{T_0}{log(1+t)} T(t)=log(1+t)T0
快速模拟退火算法的降温方式:
T ( t ) = T 0 1 + t T(t)=\frac{T_0}{1+t} T(t)=1+tT0
常用的模拟退火算法的降温方式还有(通常 0.8 < α < 0.99 0.8<\alpha<0.99 0.8<α<0.99)
T ( t + Δ t ) = α T ( t ) T(t+\Delta t)=\alpha T(t) T(t+Δt)=αT(t)
终止条件自己设定阈值即可。
花费函数COST
这个基本就是分析题目,一般设置为我们需要求解的目标函数最值。
TSP问题
已知中国34 个省会城市(包括直辖市) 的经纬度, 要求从北京出发, 游遍34 个城市, 最后回到北京. 用模拟退火算法求最短路径。
1.main.m
主函数
clc
clear;
%%
load('china.mat');
plotcities(province, border, city);%画出地图
cityaccount=length(city);%城市数量
dis=distancematrix(city)%距离矩阵
route =randperm(cityaccount);%路径格式
temperature=1000;%初始化温度
cooling_rate=0.95;%温度下降比率
item=1;%用来控制降温的循环记录次数
distance=totaldistance(route,dis);%计算路径总长度
temperature_iterations = 1;
% This is a flag used to plot the current route after 200 iterations
plot_iterations = 1;
plotroute(city, route, distance,temperature);%画出路线
%%
while temperature>1.0 %循环条件
temp_route=change(route,'reverse');%产生扰动。分子序列变化
% fprintf("%d\n",temp_route(1));
temp_distance=totaldistance(temp_route,dis);%产生变化后的长度
dist=temp_distance-distance;%两个路径的差距
if(dist<0)||(rand < exp(-dist/(temperature)))
route=temp_route;
distance=temp_distance;
item=item+1;
temperature_iterations=temperature_iterations+1;
plot_iterations=plot_iterations+1;
end
if temperature_iterations>=10
temperature=cooling_rate*temperature;
temperature_iterations=0;
end
if plot_iterations >= 20
plotroute(city, route, distance,temperature);%画出路线
plot_iterations = 0;
end
% fprintf("it=%d",item);
end
2.distance.m
计算两点之间的距离
function d = distance(lat1, long1, lat2, long2, R)
% DISTANCE
% d = DISTANCE(lat1, long1, lat2, long2, R) compute distance between points
% on sphere with radians R.
%
% Latitude/Longitude Distance Calculation:
% http://www.mathforum.com/library/drmath/view/51711.html
y1 = lat1/180*pi; x1 = long1/180*pi;
y2 = lat2/180*pi; x2 = long2/180*pi;
dy = y1-y2; dx = x1-x2;
d = 2*R*asin(sqrt(sin(dy/2)^2+sin(dx/2)^2*cos(y1)*cos(y2)));
end
3.totaldistance.m
一条route路线的总路程,也是我们需要优化的目标函数
function totaldis = totaldistance(route,dis)%传入距离矩阵和当前路线
%TOTALDISTANCE 此处显示有关此函数的摘要
% 此处显示详细说明
% fprintf("%d\n",route(1));
totaldis=dis(route(end),route(1));
% totaldis=dis(route(end),route(1));
for k=1:length(route)-1
totaldis=totaldis+dis(route(k),route(k+1));%直接加两个点之间的距离
end
4.distancematrix.m
任意两点之间的距离矩阵
function dis = distancematrix(city)
%DISTANCEMATRIX 此处显示有关此函数的摘要
% 此处显示详细说明
cityacount=length(city);
R=6378.137;%地球半径用于求两个城市的球面距离
for i = 1:cityacount
for j = i+1:cityacount
dis(i,j)=distance(city(i).lat,city(i).long,...
city(j).lat,city(j).long,R);%distance函数原来是设计来计算球面上距离的
dis(j,i)=dis(i,j);%对称,无向图
end
end
end
5.change.m
产生分子扰动,求当前解的邻解
function route = change(pre_route,method)
%CHANGE 此处显示有关此函数的摘要
% 此处显示详细说明
route=pre_route;
cityaccount=length(route);
%随机取数,相当于把0-34个城市映射到0-1的等概率抽取上,再取整
city1=ceil(cityaccount*rand); % [1, 2, ..., n-1, n]
city2=ceil(cityaccount*rand); % 1<=city1, city2<=n
switch method
case 'reverse' %[1 2 3 4 5 6] -> [1 5 4 3 2 6]
cmin = min(city1,city2);
cmax = max(city1,city2);
route(cmin:cmax) = route(cmax:-1:cmin);%反转某一段
case 'swap' %[1 2 3 4 5 6] -> [1 5 3 4 2 6]
route([city1, city2]) = route([city2, city1]);
end
6.plotcities.m
画出city
function h = plotcities(province, border, city)
% PLOTCITIES
% h = PLOTCITIES(province, border, city) draw the map of China, and return
% the route handle.
global h;
% draw the map of China
plot(province.long, province.lat, 'color', [0.7,0.7,0.7])
hold on
plot(border.long , border.lat , 'color', [0.5,0.5,0.5], 'linewidth', 1.5);
% plot a NaN route, and global the handle h.
h = plot(NaN, NaN, 'b-', 'linewidth', 1);
% plot cities as green dots
plot([city(2:end).long], [city(2:end).lat], 'o', 'markersize', 3, ...
'MarkerEdgeColor','b','MarkerFaceColor','g');
% plot Beijing as a red pentagram
plot([city(1).long],[city(1).lat],'p','markersize',5, ...
'MarkerEdgeColor','r','MarkerFaceColor','g');
axis([70 140 15 55]);
7.plotroute.m
画出路线
function plotroute(city, route, current_distance, temperature)
% PLOTROUTE
% PLOTROUTE(city, route, current_distance, temperature) plots the route and
% display current temperautre and distance.
global h;
cycle = route([1:end, 1]);
% update route
set(h,'Xdata',[city(cycle).long],'Ydata',[city(cycle).lat]);
% display current temperature and total distance
xlabel(sprintf('T = %6.1f Total Distance = %6.1f', ...
temperature, current_distance));
drawnow
一开始的路径
最后运行出来的最优化路径
有需要的朋友可以下载相关的数据集下载链接