Matlab应用笔记--模拟退火

注：本篇随笔依据《Matlab在数学建模上的应用》中第6章介绍来写，主要介绍模拟退火思想及其Matlab实现
（博客以及Matlab小白，若有不当欢迎指出）

模拟退火简介

模拟退火（SA）是一种通用概率算法，用来在一个大的搜寻空间内寻找问题的最优解。
优点：可以有效解决NP难问题，避免陷入局部最优。对初值没有强依赖关系。编程工作量小，易于实现。统计上可以保证找到全局最优解。能够处理任意程度的非线性、不连续性、随机性的目标函数。
缺点：找到最优解要耗费非常多的时间。对某个具体问题的求解需要更困难的参数调整。

模拟退火的基础原理

基本思想：

模拟固体退火的过程。将固体加到足够高的温度，再缓慢冷却到低温。在高温时，固体内部粒子无序（内能大），若降温过程足够缓慢，那么冷却中任一温度时固体都能达到热平衡，而冷却到低温时将达到这一低温下的内能最小状态。
在退火过程中，任一温度都能达到热平衡是个重要步骤，而达到平衡过程需要考虑接受新的状态。设当前状态$i$生成新状态$j$，若新状态的内能小于状态i的内能（$E_{j}<E_{i}$），则接受新状态j新的当前状态；否则以概率$e^{\frac{-(E_{j}-E_{i})}{kT}}$接受状态j，其中$k$为Boltzmann常数。
在将SA算法运用到最优化问题时，可以将温度$T$当做控制参数，函数值$f$视为内能$E$，而固体在某温度T时的一个状态对应一个解$x_{i}$。

其他一些参数的说明

退火过程由冷却进度表控制（一组初始参数），其核心是尽量使系统达到准平衡，以使算法在有限的时间内逼近最优解。
冷却进度表包括：
（1）初温$T_{0}$：冷却开始的温度。
（2）温度$T$的衰减函数：计算退火过程中的温度值的函数。
（3）最终温度$T_{f}$：停止准则。
（4）Markov链的长度$L_{k}$：任一温度$T$的迭代次数。

算法的基本步骤

其它注意点

1.状态表达

将实际问题的解用合适的数学形式表达出来。
常见有：
0-1编码（解决背包问题和指派问题等）
自然数编码（TSP问题、调度问题等）
实数编码（连续函数优化等）

2.新解的产生

要能尽量遍及空间各个区域，以便跳出局部最优解。

3.收敛的一般性条件

（1）初温够高
（2）热平衡时间够长
（3）降温够慢
（4）终止温度够低

4.参数的选择

（1）初温$T_{0}$
T0过大，计算量大；T0过小，容易陷入局部最优
（2）T的衰减函数
常用的衰减函数：$T_{k+1}=\alpha T_{k}, k=0,1,2,...$
其中$\alpha$是一个常数可以取为$0.5\sim 0.99$。
小的衰减量致使迭代次数的增加，使算法能搜索更大范围的解空间，返回更好的最终解。
同时由于在$T_{k}$值上已经到达准平衡，则在$T_{k+1}$时只需要少量的变换就可到达准平衡。这样就能选取较短长度的Markov链来减少算法时间。
（3）Markov链长度
T的衰减函数已经选定的前提下，$L_{k}$应能使在$T$的每一取值上达到准平衡。
对于简单的情况，可令$L_{k}=100n$，$n$为问题的规模。
（4）终止温度$T_{f}$
一般$T_{f}$设为一个足够小的整数，比如$0.01\sim 5$（粗糙经验）。

例题以及函数模版

模拟退火

模拟退火解决TSP问题的一般步骤

（1）确定解空间和初始解
解空间：城市的所有排列构成的集合。
初始解：随机生成一个排列（TSP对初始解的依赖低）。
（2）目标函数
某排列下的路径长度。
（3）新解产生
1.二变换法：任选序号$u,v$（设$u<v<n$），交换$u$和$v$之间的访问顺序。
2.三变换法：任选序号$u,v,w$（设$u\leqslant v<w$），将$u$和$v$之间的路径插到$w$之后访问。
（4）目标函数差
计算变换前的解和变换后目标函数的差值
（5）Metropolis接受准则
若$\Delta \leqslant 0$则接受新解
否则以$e^{\frac{-\Delta }{T}}$的概率接受新解。

经典TSP问题，数据来源于“TSPLIB95”中的“berlin52”

-----------
（52座城市）

解决的函数模版

clear;clc;
 
a = 0.99;   %温度衰减函数的参数
t0 = 97;tf = 3;t = t0;
Markov_length = 10000;  %Markov链长度
 
%读取tsp文件的数据
fid = fopen('data\TSP_Pro\berlin52.tsp','r');
FormatString = repmat('%f',1,3);
coordinates = cell2mat(textscan(fid,FormatString,52,'HeaderLines',6));
 
coordinates(:,1) = [];      %去掉数据中城市的序号
amount = size(coordinates,1);   %城市的数目
%通过向量化的方法计算距离矩阵
coor_x_tmp1 = coordinates(:,1)*ones(1,amount);
coor_x_tmp2 = coor_x_tmp1';
coor_y_tmp1 = coordinates(:,2)*ones(1,amount);
coor_y_tmp2 = coor_y_tmp1';
%dist_matrix(i,j)即代表了i号城市与j号城市间的距离
dist_matrix = sqrt((coor_x_tmp1-coor_x_tmp2).^2+...
    (coor_y_tmp1-coor_y_tmp2).^2);
 
%sol_new是每次产生的新解；sol_current是当前解；sol_best是冷却中最好的解
%E_new是新解的回路距离；E_current是当前解的回路距离；E_best是最优解的回路距离
sol_new = 1:amount; %产生初始解
sol_current = sol_new;
sol_best = sol_new;
E_current = inf;
E_best = inf;
 
while t>=tf
    for r = 1:Markov_length     %Markov链长度
        %随机扰动产生新解
        if(rand < 0.5)  %随机决定是进行两交换还是三交换
            %两交换
            ind1 = 0;ind2 = 0;
            while(ind1 == ind2)
                ind1 = ceil(rand.*amount);
                ind2 = ceil(rand.*amount);
            end
            tmp1 = sol_new(ind1);
            sol_new(ind1) = sol_new(ind2);
            sol_new(ind2) = tmp1;
        else
            %三交换
            ind1 = 0;ind2 = 0;ind3 = 0;
            while(abs(ind1-ind2)<2||abs(ind1-ind3)<2 ...
                  ||abs(ind2-ind3)<2)
                ind1 = ceil(rand.*amount);
                ind2 = ceil(rand.*amount);
                ind3 = ceil(rand.*amount);
            end
            tmp1 = ind1;tmp2 = ind2;tmp3 = ind3;
            
            if (ind1 < ind2)&&(ind2 < ind3)      %确保ind1 < ind2 < ind3
            elseif (ind1 < ind3)&&(ind3 < ind2)
                ind2 = tmp3;ind3 = tmp2;
            elseif (ind2 < ind1)&&(ind1 < ind3)
                ind1 = tmp2;ind2 = tmp1;
            elseif (ind2 < ind3)&&(ind3 < ind1)
                ind1 = tmp2;ind2 = tmp3;ind3 = tmp1;
            elseif (ind3 < ind1)&&(ind1 < ind2)
                ind1 = tmp3;ind2 = tmp1;ind3 = tmp2;
            else
                ind1 = tmp3;ind2 = tmp2;ind3 = tmp1;
            end
            
            tmplist1 = sol_new((ind1+1):(ind2-1));
            sol_new((ind1+1):(ind1+ind3-ind2+1)) = ...
                sol_new((ind2):(ind3));
            sol_new((ind1+ind3-ind2+2):ind3) = ...
                tmplist1;
        end
        
        %计算目标函数值（即内能）
        E_new = 0;
        for i = 1:(amount-1)
            E_new = E_new+...
                dist_matrix(sol_new(i),sol_new(i+1));
        end
        E_new = E_new+...      %再算上从最后一个城市到第一个城市的距离
            dist_matrix(sol_new(amount),sol_new(1));
        
        %考虑接受新解
        if E_new<E_current
            E_current = E_new;
            sol_current = sol_new;
            if E_new<E_best
                %把冷却过程中最好的解保留下来
                E_best = E_new;
                sol_best = sol_new;
            end
        else
            %若新解的目标函数值小于当前解的，
            %则仅以一定概率接受新解
            if rand<exp(-(E_new-E_current)./t)
                E_current = E_new;
                sol_current = sol_new;
            else
                sol_new = sol_current;
            end
        end
    end
    t = t.*a;       %控制参数t（温度）减少为原来的a倍
end
 
disp('最优解为：')
disp(sol_best)
disp('最短距离')
disp(E_best)

实际测试结果

实机跑：

实际最优结果：

途径城市顺序：
1 49 32 45 19 41 8 9 10 43 33 51 11 52 14 13 47 26 27 28 12 25 4 6 15 5 24 48 38 37 40
39 36 35 34 44 46 16 29 50 20 23 30 2 7 42 21 17 3 18 31 22
最短路径长度

相关资料

TSP标准测试算例的网站 http://comopt.ifi.uni-heidelberg.de/software/TSPLIB95/tsp/
其中当前求得的最优解可查询 http://comopt.ifi.uni-heidelberg.de/software/TSPLIB95/STSP.html
使用方法 https://blog.csdn.net/junzhepan/article/details/8498707
按行读取文本文件数据 https://www.cnblogs.com/sinux/p/3456577.html

01背包问题

clear;clc;
 
a = 0.95;
k = [5;10;13;4;3;11;13;10;8;16;7;4];
k = -k;
d = [2;5;18;3;2;5;10;4;11;7;14;6];
restriction = 46;
num = 12;
sol_new = ones(1,num);  %生成初始解
E_curr = inf;
E_best = inf;
%E_curr是当前解对应的目标函数值
%E_new是新解的目标函数值
%E_best是最优解
sol_curr = sol_new;
sol_best = sol_new;
t0 = 97;tf = 3;t = t0;
p = 1;
 
while t >= tf
    for r = 1:100
        %产生随机扰动
        tmp = ceil(rand.*num);
        sol_new(1,tmp) = ~sol_new(1,tmp);
        %检查是否满足约束
        while 1
            q = (sol_new*d <= restriction);
            if ~q
                p = ~p; %实现交错逆转头尾的首个1
                tmp = find(sol_new == 1);
                if p
                    sol_new(1,tmp) = 0;
                else
                    sol_new(1,tmp(end)) = 0;
                end
            else
                break
            end
        end
        
        %计算背包中的物品价值
        E_new = sol_new*k;
        if E_new < E_curr
            E_curr = E_new;
            sol_curr = sol_new;
            if E_new < E_best
                %把冷却过程中最好的解保存下来
                E_best = E_new;
                sol_best = sol_new;
            end
        else
            if rand < exp(-(E_new-E_curr)./t)
                E_curr = E_new;
                sol_curr = sol_new;
            else
                sol_new = sol_curr;
            end
        end
    end
    t = t.*a;
end
 
disp('最优解为：')
sol_best
disp('物品总价值等于：')
val = -E_best;
disp(val)
disp('背包中物品重量是：')
disp(sol_best*d)

实际运行结果

动态解决01背包问题
https://www.cnblogs.com/kkbill/p/12081172.html

关于ASA程序包和ASA对MATLAB的接口程序ASAMIN的使用

（由于大都是外网链接，所以加载会很慢，若加载错误可刷新尝试一下）

（1）在Matlab中配置mex（若有配置，则忽略）

本人使用2018的Matlab，可能当初安装时没有安装MinGW-w64 C/C++ 编译器，在Matlab命令行窗口输入mex -setup后报错“未找到支持的编译器”，即尚未配置mex。
我的解决方法（如果不成功，可以看看下载链接里评论区其他人的做法，我之前也失败了好几次）：
于是按照提示指引，注册账号（已有则不用），用浏览器打开 https://ww2.mathworks.cn/matlabcentral/fileexchange/52848-matlab-support-for-mingw-w64-c-c-compiler
点击Download下载mingw（mingw.mlpkginstall），然后放到桌面上。
在matlab已打开的情况下，双击打开mingw然后按指示操作，Matlab自动帮你下载安装好。

（下载好的）

（2）下载ASA程序包和ASAMIN程序包

ASA程序包的下载地址 https://www.ingber.com/#ASA
找到ASA.zip 下载解压

ASAMIN程序的下载地址 http://ssakata.sdf.org/software/

从解压后的ASA程序包中复制asa.c、asa.h、asa_usr_asa.h到ASAMIN所在的文件夹。

（3）编译生成mex文件，并执行该文件

将Matlab的工作路径转到ASAMIN所在的文件夹。
在Matlab命令行窗口中输入 mex asamin.c asa.c -DUSER_ACCEPTANCE_TEST=TRUE -DUSER_ASA_OUT=TRUE -DDBL_MIN=2.2250738585072014e-308
若报错，可参考以下两个链接中的处理
https://blog.csdn.net/qq_32439305/article/details/103201515
https://blog.csdn.net/summerdream_/article/details/86422024
在Matlab命令行窗口中输入asatest运行该文件，生成asatest1.log和asatest2.log两个结果文件。
这两个结果文件的内容可以与ASA文件下的asa_test_asa文件进行对比，以验证ASAMIN的安装。

（4）设置搜索路径

移动ASAMIN的mex文件和asamin.m文件到Matlab的搜索目录下。
（我个人是在Matlab的toolbox目录下新建asamin文件夹，把那两个文件拖进去，然后再设置路径把新建的asamin文件添加进去）

（5）使用方法

至此ASAMIN的安装就完成了。
关于查看及调整参数列表可以在ASAMIN文件夹下的asamin.m查看。
ASAMIN的函数说明如下：

[fstar, xstar, grad, hessian, state] = ...
  asamin('minimize', func, xinit, xmin, xmax, xtype,parm1,parm2,...);

其它使用ASA的方法

ASAMIN是方便Matlab使用ASA，若想用别的方法调用ASA，还有别的方法。
这里只简单谈书上说的使用Cygwin对ASA编译。
关于Cygwin的安装教程 https://www.douban.com/note/723078223/

这里编译ASA中自带的测试例子asa_usr_cst.c文件，ASA的参数可以以ASA-Makefile为模版进行设置，asa_usr_cst.c中包含用户定义的目标函数及约束条件，其它详细的使用方法可以查看ASA程序包中的说明文档。
将ASA移动到Cygwin的home下的用户文件夹内，在Cygwin的lib文件下找到cygwin1.dll将其拷贝到ASA同级目录下

把运行Cygwin（Windows 批处理文件 (.bat)），输入如下三个指令。

可以看到home下用户文件内生成的程序以及程序运行生成的两个文件（红框内）。
可以文本打开两个绿框内的两个文件的内容分别进行对比结果。