现代优化算法

启发式算法包括了：禁忌搜索，模拟退火，遗传算法，人工神经网络。目标：求NP-hard组合优化问题的全局最优解。解决复杂优化问题的蚁群算法，根据实际问题产生的解空间分解、解空间限制，还有对多个启发式算法的合成的集成算法。

主要用于解决组合优化问题：TSP,QAP,JSP

1 模拟退化算法

1.1 算法简介

该算法主要思想来自于材料的统计力学，高温条件下粒子的能量较高，可以自由运动和重新排列。而低温条件下，粒子能量较低。如果从高温开始，缓慢降温（退火），粒子就可以在每个温度下达到热平衡。系统完全冷却时，最终形成处于低能状态的晶体。

退火算法用下列数学模型描述退火过程。

假设材料在状态$i$之下的能量为$E(i)$，则材料在温度$T$时从状态$i$进入状态$j$遵循以下规律：
（1）如果$E(j) \le E(i)$，接受该状态被转换

（2）如果$E(j) \gt E(i)$，则状态转换以如下概率被接受：

$$ e^{\frac{E(i)-E(j)}{KT}} $$

其中$K$为玻尔兹曼常数，$T$是材料温度

在某个特定的温度下，经过充分转换之后，材料达到热平衡，这时候材料处于状态$i$的概率满足玻尔兹曼分布：

$$ P_{T}(x=i)=\frac{e^-{\frac{E(i)}{KT}}}{\sum_{j \in S}e^-{\frac{E(j)}{KT}}} $$

其中$x$表示材料当前状态的随机变量，$S$表示状态空间集合。

显然有：

$$ \lim_{T \to \infty}\frac{e^-{\frac{E(i)}{KT}}}{\sum_{j \in S}e^-{\frac{E(j)}{KT}}}=\frac{1}{|S|} $$

其中$|S|$表示集合$S$中状态的数量。这表明所有状态在高温下具有相同的概率。而当温度下降时，

$$ \lim_{T \to 0}\frac{e^-{\frac{E(i)-E_{min}}{KT}}}{\sum_{j \in S}e^-{\frac{E(j)-E_{min}}{KT}}}=\lim_{T \to 0}\frac{e^-{\frac{E(i)-E_{min}}{KT}}}{\sum_{j \in S_{min}}e^-{\frac{E(j)-E_{min}}{KT}}+\sum_{j \notin S_{min}}e^-{\frac{E(j)-E_{min}}{KT}}}\\ =\lim_{T \to 0}\frac{e^-{\frac{E(i)-E_{min}}{KT}}}{\sum_{j \in S_{min}}e^-{\frac{E(j)-E_{min}}{KT}}}=\begin{cases} \frac{1}{|S_{min}|} \quad 若i \in S_{min}\\ 0 \quad 其它 \end{cases} $$

其中$E_{min}=\min_{j \in S}E(j)$且$S_{min}=\left\{i|E(i)=E_{min}\right\}$

上式表明当温度降至很低时，材料会以很大概率进入最小能量状态。

1.2 模拟退火寻优方法：

解是不断由退火生成的，但是是以概率来确定的。并且是一个马尔可夫过程，转移概率服从均匀分布。不断地进行状态转移，直到温度降低为0。而根据理论分析，模拟退火算法可以找到全局最优解。

要注意的问题有：
（1），理论上，降温过程要足够缓慢，使得在每个温度下达到热平衡，才能不至于跳过任何一个可能的解。不过会导致算法比较慢，正常选用的时候都需要折衷。

（2），要确定在每一温度下状态的结束准则。实际操作可以考虑当连续$m$次的转换过程没有使状态发生改变的时候，结束该温度下的状态转换。 最终温度的确定可以提前定为一个较小的值$T_{e}$，或连续几个温度下转换过程没有使状态发生变化算法就结束。

（3），选择初始温度和确定某个可行解的邻域的方法也要恰当。

1.3 举个例子

例子：已经知道敌方的100个目标的经纬度（具体数据可以找我索取）

我方有一个基地，经度和纬度为（70，40）。假设我方飞机的速度为 1000 公里/小时。我方派一架飞机从基地出发，侦察完敌方所有目标，再返回原来的基地。在敌方每一目标点的侦察时间不计，求该架飞机所花费的时间（假设我方飞机巡航时间可以充分长）。

这是一个旅行商问题。我们依次给基地编号为 1，敌方目标依次编号为 2， 3，…，101， 最后我方基地再重复编号为 102（这样便于程序中计算）。 距离矩阵$D=(d_{ij})_{102*102}$，其中$d_{ij}$表示$i,j$两点的距离，$i,j=1,2,\cdots,102$，这里$D$为实对称矩阵。则问题变成了从点1出发，走遍中间的所有点，到达点120的一个最短路径。

但是由于给定的是地理坐标，我们要转化为求两点之间的实际距离。设$A,B$两点的地理坐标为$(x_1,y_1),(x_2,y_2)$，过$A,B$两点的 大圆的劣弧长即为两点的实际距离。以地心为坐标原点$O$,以赤道平面为$XOY$平面，以0度经线圈所在的平面为$XOZ$平面建立三维直角坐标系。则$A,B$两点的直角坐标分别为：

$$ A(Rcosx_1cosy_1,Rsinx_1cosy_1,Rsiny_1)\\ B(Rcosx_2cosy_2,Rsinx_2cosy_2,Rsiny_2)\\ $$

其中$R=6370$为地球半径。

$A,B$两点的实际距离

$$ d=Rarccos(\frac{\vec{OA}{\cdot}\vec{OB}}{|\vec{OA}|{\cdot}|\vec{OB}|}) $$

化简得到：

$$ d=Rarccos[cos(x_1-x_2)cosy_1cosy_2+siny_1siny_2]。 $$

1.4 模拟退火算法描述如下：

（1）解空间

解空间$S$可表为$\left\{1,2,\cdots,101,102\right\}$的所有固定起点和重点的循环排列集合，即$S=\left\{\pi_1,\cdots,\pi_{102}|\pi_1=1,(\pi_2,\cdots,\pi_{101})\right\}$为$\left\{2,3,\cdots,101\right\}$的循环排列，$\pi_{102}=102$,其中每一个循环表示侦察100个目标的一个回路，$\pi_i=j$表示在第$i$次侦察$j$点，初始解可选为(1,2,...,102)

（2）目标函数

此时的目标函数为侦察所有目标的路径长度或者称为代价函数，我们要求：

$$ \min f(\pi_1,\pi_2,\cdots,\pi_{102})=\sum^101_{i=1}d_{\pi_i\pi_{i+1}}\\ $$

而一次迭代由下列三步构成；

（3）新解的产生

2变换法：

任选序号$u,v(u \lt v)$交换$u$与$v$之间的顺序，此时的新路径为：

$$ \pi_1\cdots\pi_{u-1}\pi_v\pi_{v+1}\cdots\pi_{u+1}\pi_u\pi_{v+1}\cdots\pi_{102}\\ $$

（4）代价函数差

这里采用2变换法得到的路径差表示为：

$$ \Delta f=(d_{\pi_{u-1}\pi_v}+d_{\pi_u\pi_{v+1}})-(d_{\pi_{u-1}\pi_u}+d_{\pi_v\pi_{v+1}}) $$

（5）接受准则

$$ P=\begin{cases} 1\quad \Delta f \lt0\\ exp(-\Delta f/T) \quad \Delta f \ge 0\\ \end{cases} $$

如果$\Delta f \lt 0$则接受新的路径，否则以概率$exp(-\Delta f/T)$接受新的路径，即若$exp(-\Delta f/T)$大于0到1之间的随机数则接受。

（6）降温

利用选定的降温系数$\alpha$进行降温即：$T \gets \alpha{T}$,得到新的温度，这里我们取$\alpha=0.999$

（7）结束条件

用选定的终止温度$e=10^{-30}$，判断退火过程是否结束。若$T \lt e$,算法结束，输出当前状态。

1.5 Matlab算法程序

matlab程序代码如下：

clc,clear;
load sj.txt %加载敌方 100 个目标的数据， 数据按照表格中的位置保存在纯文本文件 sj.txt 中,为25*8的矩阵格式
x=sj(:,1:2:8);x=x(:);  %取1到8列中步长为2的列数，并且转换为列向量。（原本为25*4的矩阵格式）
y=sj(:,2:2:8);y=y(:);  %取2到8列中步长为2的列数,并且转换为列向量。
sj=[x y];              %转换为矩阵100*2的经纬度矩阵
d1=[70,40];            %我方基地的经纬度
sj=[d1;sj;d1];         %加入我方基地，从基地绕一圈再回到基地。所以102*2矩阵
sj=sj*pi/180;          %转换为弧度制
%距离矩阵 d
d=zeros(102);          %给出一个基本的零矩阵初始化以存储距离
for i=1:101             %迭代行数
 for j=i+1:102          %迭代列数
 temp=cos(sj(i,1)-sj(j,1))*cos(sj(i,2))*cos(sj(j,2))+sin(sj(i,2))*sin(sj(j,2));   %两点间距离公式（三维）
 d(i,j)=6370*acos(temp);   %地球半径6370
 end
end
d=d+d';             %搞出实对称矩阵，两两基地之间的距离
S0=[];Sum=inf;      %初始化走过的结点，初始化路径长度
rand('state',sum(clock));   %根据当前时间设定初始状态随机数器
for j=1:1000
 S=[1 1+randperm(100),102];  %s是一维数组，第一位是1，中间是2到101的随机序列，最后一位是102
 temp=0;                     %临时变量设为0
 for i=1:101             
 temp=temp+d(S(i),S(i+1));   %临时变量是距离，不断遍历这些飞机飞过的点
 end
 if temp<Sum                 %在循环过程中，如果临时距离小于Sum，则下面两个语句执行，也就是找到飞过这些路线的最短距离。           
 S0=S;Sum=temp;
 end
end
e=0.1^30;L=20000;at=0.999;T=1;   %退火参数
%退火过程
for k=1:L   %我tm迭代个2w次
 %产生新解
c=2+floor(100*rand(1,2));  %产生102之间的两个数字——经纬度
c=sort(c);                 %进行升序排列
c1=c(1);c2=c(2);           %取出经纬度
 %计算代价函数（也就是路径差）
 df=d(S0(c1-1),S0(c2))+d(S0(c1),S0(c2+1))-d(S0(c1-1),S0(c1))-d(S0(c2),S0(c2+1));
 %接受准则
 if df<0            %代价函数差接受路径
 S0=[S0(1:c1-1),S0(c2:-1:c1),S0(c2+1:102)];
 Sum=Sum+df;        %更正路径长度
 elseif exp(-df/T)>rand(1)  %按照概率接受准则接受新路径
 S0=[S0(1:c1-1),S0(c2:-1:c1),S0(c2+1:102)];   %接受路径
 Sum=Sum+df;        %更正路径长度
 end
 T=T*at;     %退火计算得到新的温度
 if T<e      %直到冷却
 break;
 end
end
 % 输出巡航路径及路径长度
S0,Sum

本质上是通过用随机序列爆算出最小长度，然后再通过退火算法的2代换，用随机序列再次计算代价函数差，然后根据温度退火决定后面的2代换次数，最后冷却输出最小的距离长度。