现代优化算法(模拟退火和遗传算法-matlab)

一、概述

优化算法
是一种用于在给定问题的解空间中寻找最优解的方法。
现代优化算法
强调全局搜索能力和局部搜索能力的结合，利用概率论、统计学和生物学等学科的理论。

参考资料

现代优化算法 （一）：模拟退火算法 及应用举例
现代优化算法 （二）： 遗传算法 及应用举例
改进的遗传算法
现代优化算法
现代智能优化算法

现代优化算法概述

现代优化算法是上世纪 80 年代初兴起的启发式算法,包括:

• 禁忌搜索（tabu search）
• 模拟退火（simulated annealing）
• 遗传算法（genetic algorithms）
• 人工神经网络（neural networks）

现代优化算法的特点：

• 具有全局优化性能
• 通用性强
• 适合于并行处理

NP-hard 问题

这些算法在理论和实际应用方面得到了较大的发展。无论这些算法是怎样产生的，它们有一个共同的目标——求 NP-hard 组合优化问题的全局最优解。尽管这些目标，但 NP-hard 理论限制它们只能以启发式的算法去求解实际问题。

启发式算法

启发式算法包含的算法很多，例如解决复杂优化问题的蚁群算法（Ant Colony Algorithms）。有些启发式算法是根据实际问题而产生的，如解空间分解、解空间的限制等；另一类算法是集成算法，这些算法是诸多启发式算法的合成。

组合优化问题

现代优化算法解决组合优化问题，如 TSP（Traveling Salesman Problem）问题，QAP（Quadratic Assignment Problem）问题，JSP（Job-shop Scheduling Problem）问题等效果很好。

旅行商问题:TSP（Traveling Salesman Problem）

组合优化中的经典问题之一。其定义是：给定一组城市以及这些城市之间的距离，要求找到一条经过每个城市一次且仅一次的最短路径，并返回起始城市。TSP 广泛应用于物流配送、电子电路设计和基因排序等领域。

二次分配问题:QAP（Quadratic Assignment Problem）

一个复杂的组合优化问题，常用于设施布局设计。其定义是：给定 n 个设施和 n 个位置，以及设施之间的流量和位置之间的距离，目标是找到一种设施到位置的分配方式，使得总的流量-距离乘积最小。QAP 广泛应用于工厂布局设计、键盘布局设计和计算机芯片布局等。

作业车间调度问题:JSP（Job-shop Scheduling Problem）

生产调度中的经典问题。其定义是：给定一组作业，每个作业由一系列特定顺序的任务组成，每个任务需要在特定的机器上处理，目标是找到一种调度方案，使得所有作业的完成时间最短。JSP 广泛应用于制造业中的生产调度、项目管理和软件开发等领域。

二、相关模型

1.模拟退火算法

参考资料

算法简介

模拟退火算法得益于材料统计力学的研究成果。统计力学表明材料中粒子的不同结构对应于粒子的不同能量水平。

• 高温条件下：粒子的能量较高，可以自由运动和重新排列。
• 低温条件下：粒子能量较低。

如果从高温开始，缓慢降低温度（这一过程被称为退火），粒子可以在每个温度下达到热平衡。当系统完全被冷却时，最终形成处于低能状态的晶体。

Metropolis 算法

如果用粒子的能量定义材料的状态，Metropolis 算法用一个简单的数学模型描述了退火过程。

假设材料在状态 $i$ 之下的能量为 $E(i)$，那么材料在温度 $T$时从状态 $i$ 进入状态 $j$ 就遵循如下规律：

1. 若 $E(j)\le E(i)$，则接受该状态被转换。
2. 若 $E(j)>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$ 表示状态空间集合。

显然：

$$\underset{T\to \infty }{\lim }\frac{e^{-\frac{E(i)}{KT}}}{{\sum }_{j\in S}{e}^{-\frac{E(j)}{KT}}}=\frac{1}{|S|}$$

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

当温度下降时，有：

$$\underset{T\to 0}{\lim }\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}}}=\{\begin{array}{ll}\frac{1}{|S_{\min }|},& \text{若 }i\in S_{\min }\\ 0,& \text{其他}\end{array}$$
其中：
$$E_{\min }=\underset{j\in S}{\min }E(j)$$
且：
S min ⁡ = { i ∣ E ( i ) = E min ⁡ } S_{\min} = \{i | E(i) = E_{\min}\} Smin​={i∣E(i)=Emin​}
上式表明当温度降至很低时，材料会以很大概率进入最小能量状态。

优化问题中的模拟退火

假设要解决的问题是一个寻找最小值的优化问题。将物理学中模拟退火的思想应用于优化问题就可以得到模拟退火优化方法。
考虑这样一个组合优化问题：优化函数为 $$f:x\to {\mathbb{R}}^{+}$$其中 $x\in S$，它表示优化问题的一个可行解 R + = { y ∣ y ∈ R , y ≥ 0 } \mathbb{R}^+ = \{y | y \in \mathbb{R}, y \geq 0\} R+={y∣y∈R,y≥0}$S$ 表示函数的定义域。$$N(x)\subseteq S$$ 表示$x$的一个邻域集合。
首先给定一个初始温度 $T_0$ 和该优化问题的一个初始解 $x(0)$，并由 $x(0)$ 生成下一个解 $x^{\prime }\in N(x(0))$ ，是否接受 $x^{\prime }$ 作为一个新解 $x(1)$ 依赖于下面概率：
$$P(x(0)\to x^{\prime })=\{\begin{array}{ll}1,& \text{若 }f(x^{\prime })<f(x(0))\\ e^{-\frac{f(x^{\prime })-f(x(0))}{T_0}},& \text{其他}\end{array}$$
换句话说，如果生成的解 $x^{\prime }$ 的函数值比前一个解的函数值更小，则接受 $x(1)=x^{\prime }$ 作为一个新解。否则以概率 $$e^{-\frac{f(x^{\prime })-f(x(0))}{T_0}}$$ 接受 $x^{\prime }$ 作为一个新解。
在 $T_{i+1}$ 下重复上述过程。因此整个优化过程就是不断寻找新解和缓慢降温的交替过程。最终的解是对该问题寻优的结果。
注意到在每个 $T_i$ 下，所得到的一个新状态 $x(k+1)$ 完全依赖于前一个状态 $x(k)$，和前面的状态 $x(0),\dots ,x(k-1)$ 无关，因此这是一个马尔可夫过程。使用马尔可夫过程的上述模拟退火的步骤进行分析，结果表明从任何一个状态 $x(k)$ 生成 $x^{\prime }$ 的概率，在 $N(x(k))$ 中是均匀分布的，且新状态 $x^{\prime }$ 被接受的概率满足下式 ，那么经过有限次的转换，在温度 $T_i$ 下的平衡态 $x_i$ 的分布由下式给出：
$$P_i(T_i)=\frac{e^{-\frac{f(x_i)}{T_i}}}{{\sum }_{j\in S}{e}^{-\frac{f(x_j)}{T_i}}}$$
当温度 $T$ 降为 0 时，$x_i$ 的分布为：
$$P_i^{\ast }=\{\begin{array}{ll}\frac{1}{|S_{\min }|},& \text{若 }{x}_i\in S_{\min }\\ 0,& \text{其他}\end{array}$$
并且：
$$\sum _{x_i\in S_{\min }}{P}_i^{\ast }=1$$
这说明如果温度下降得十分缓慢，而在每个温度都有足够多次的状态转移，使之在每一个温度下达到热平衡则全局最优解将以概率 1 被找到。因此可以说模拟退火算法可以找到全局最优解。

注意:和马尔可夫过程一样,无后效性,即当前过程之于以前的过程有关,与以后结果无关

模拟退火算法中的注意事项

1. 理论上，降温过程要足够缓慢，要使得在每一温度下达到热平衡。但在计算机实现中，如果降温速度过缓，所得到的解的性能会较为令人满意，但是算法会太慢，相对于简单的搜索算法不具有明显优势。如果降温速度过快，很可能最终得不到全局最优解。因此使用时要综合考虑解的性能和算法速度，在两者之间采取一种折衷。
2. 要确定在每一温度下状态转换的结束准则。实际操作可以考虑当连续 $m$ 次的转换过程中没有使状态发生变化时结束该温度下的状态转换。最终温度的确定可以提前定为一个较小的值 $T_e$
   或连续几个温度下转换过程中没有使状态发生变化算法就结束。
3. 选择初始温度和确定某个可行解的邻域的方法也要恰当。

应用举例

例:已知 100 个目标的经度、纬度如表 12.1 所示。

53.7121   15.3046   51.1758    0.0322   46.3253   28.2753   30.3313    6.9348
56.5432   21.4188   10.8198   16.2529   22.7891   23.1045   10.1584   12.4819
20.1050   15.4562   1.9451    0.2057    26.4951   22.1221   31.4847    8.9640
26.2418   18.1760   44.0356   13.5401   28.9836   25.9879   38.4722   20.1731
28.2694   29.0011   32.1910    5.8699   36.4863   29.7284   0.9718   28.1477
8.9586   24.6635    16.5618   23.6143   10.5597   15.1178   50.2111   10.2944
8.1519    9.5325    22.1075   18.5569   0.1215   18.8726    48.2077   16.8889
31.9499   17.6309   0.7732    0.4656    47.4134   23.7783   41.8671    3.5667
43.5474    3.9061   53.3524   26.7256   30.8165   13.4595   27.7133    5.0706
23.9222    7.6306   51.9612   22.8511   12.7938   15.7307   4.9568    8.3669
21.5051   24.0909   15.2548   27.2111   6.2070    5.1442    49.2430   16.7044
17.1168   20.0354   34.1688   22.7571   9.4402    3.9200    11.5812   14.5677
52.1181    0.4088   9.5559   11.4219    24.4509    6.5634   26.7213   28.5667
37.5848   16.8474   35.6619    9.9333   24.4654    3.1644   0.7775    6.9576
14.4703   13.6368   19.8660   15.1224   3.1616    4.2428    18.5245   14.3598
58.6849   27.1485   39.5168   16.9371   56.5089   13.7090   52.5211   15.7957
38.4300    8.4648   51.8181   23.0159   8.9983   23.6440    50.1156   23.7816
13.7909    1.9510   34.0574   23.3960   23.0624    8.4319   19.9857    5.7902
40.8801   14.2978   58.8289   14.5229   18.6635    6.7436   52.8423   27.2880
39.9494   29.5114   47.5099   24.0664   10.1121   27.2662   28.7812   27.6659
8.0831   27.6705    9.1556   14.1304    53.7989    0.2199   33.6490    0.3980
1.3496   16.8359    49.9816    6.0828   19.3635   17.6622   36.9545   23.0265
15.7320   19.5697   11.5118   17.3884   44.0398   16.2635   39.7139   28.4203
6.9909   23.1804    38.3392   19.9950   24.6543   19.6057   36.9980   24.3992
4.1591    3.1853    40.1400   20.3030   23.9876    9.4030   41.1084   27.7149

我方有一个基地，经度和纬度为 $(70,40)$。假设我方飞机的速度为 1000 公里/小时。我方派一架飞机从基地出发，侦察完所有目标，再返回原来的基地。在每一目标点的侦察时间不计，求该架飞机所花费的时间（假设我方飞机巡航时间可以充分长）。
这是一个旅行商问题。给我方基地编号为 1，目标依次编号为 2, 3, …, 101，最后我方基地再重复编号为 102（这样便于程序中计算）。距离矩阵 $D=(d_{ij}{)}_{102\times 102}$，其中 $d_{ij}$ 表示 $i,j$ 两点的距离， $i,j=1,2,...,102$，这里 $D$ 为实对称矩阵。则问题是求一个从点 1 出发，走遍所有中间点，到达点 102 的一个最短路径。
上面问题中给定的是地理坐标（经度和纬度），必须求两点间的实际距离。设 $A,B$ 两点的地理坐标分别为 $(x_1,y_1)$，$(x_2,y_2)$，过 $A,B$ 两点的大圆的劣弧即为两点的实际距离。以地心为坐标原点 $O$，以赤道平面为 $XOY$ 平面，以 0 度经线圈所在的平面为 $XOZ$ 平面建立三维直角坐标系。则 $A,B$ 两点的直角坐标分别为：
$$A(R\cos x_1\cos y_1,R\sin x_1\cos y_1,R\sin y_1)$$
$$B(R\cos x_2\cos y_2,R\sin x_2\cos y_2,R\sin y_2)$$
其中 $R=6370\text{km}$ 为地球半径。

模拟退火算法的详细步骤

1. 解空间

解空间 $S$ 可表为 { 1 , 2 , ⋯   , 101 , 102 } \{1, 2, \cdots, 101, 102\} {1,2,⋯,101,102} 的所有固定起点和终点的循环排列集合，即
S = { ( π 1 , ⋯   , π 102 ) ∣ π 1 = 1 , ( π 2 , ⋯   , π 101 ) 为 { 2 , 3 , ⋯   , 101 } 的循环排列 , π 102 = 102 } S = \{(\pi_1, \cdots, \pi_{102}) | \pi_1 = 1, (\pi_2, \cdots, \pi_{101}) 为 \{2, 3, \cdots, 101\} 的循环排列, \pi_{102} = 102\} S={(π1​,⋯,π102​)∣π1​=1,(π2​,⋯,π101​)为{2,3,⋯,101}的循环排列,π102​=102}
其中每一个循环排列表示侦察 100 个目标的一个回路，${\pi }_i=j$ 表示在第 $i-1$ 次侦察目标 $j$。初始解可选为 { 1 , 2 , ⋯   , 102 } \{1, 2, \cdots, 102\} {1,2,⋯,102}，我们先使用 Monte Carlo（蒙特卡洛）方法来得到一个较好的初始解。

2. 目标函数

目标函数（或称代价函数）为侦察所有目标的路径长度。要求
$$\min f({\pi }_1,{\pi }_2,\cdots ,{\pi }_{102})=\sum _{i=1}^{101}{d}_{{\pi }_i,{\pi }_{i+1}}$$
而一次迭代由下列三步构成。

3. 新解的产生

设上一步迭代的解为
$${\pi }_1,{\pi }_2,\cdots ,{\pi }_{u-1},{\pi }_u,{\pi }_{u+1},\cdots ,{\pi }_v,{\pi }_{v+1},\cdots ,{\pi }_w,{\pi }_{w+1},\cdots ,{\pi }_{102}$$

i. 2 交换法

任选序号 $u,v$，交换 $u$ 与 $v$ 之间的顺序，变成逆序，此时的新路径为
$${\pi }_1,\cdots ,{\pi }_{u-1},{\pi }_v,{\pi }_{v-1},\cdots ,{\pi }_u,{\pi }_{u+1},\cdots ,{\pi }_w,\cdots ,{\pi }_{102}$$

ii. 3 交换法

任选序号 $u,v,w$，将 $u$ 和 $v$ 之间的路径插到 $w$ 之后，对应的新路径为
$${\pi }_1,\cdots ,{\pi }_{u-1},{\pi }_{v+1},\cdots ,{\pi }_w,{\pi }_u,\cdots ,{\pi }_v,{\pi }_{w+1},\cdots ,{\pi }_{102}$$

4. 代价函数差

对于 2 交换法，路径差可表示为
$$\Delta f=(d_{{\pi }_{u-1},{\pi }_u}+d_{{\pi }_v,{\pi }_{v+1}})-(d_{{\pi }_{u-1},{\pi }_v}+d_{{\pi }_u,{\pi }_{v+1}})$$

5. 降温

利用选定的降温系数 $\alpha$ 进行降温，取新的温度 $T$ 为 $\alpha T$（这里 $T$ 为上一步迭代的温度），这里选定 $\alpha =0.999$。

6. 结束条件

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

为什么要进行 2,3 交换法?

用于产生新解的两种邻域搜索方法。它们的主要目的是通过在当前解的基础上进行局部调整，探索解空间中的其他可能解，从而逐步接近全局最优解。

2 交换法（2-opt）

2 交换法是最简单的邻域搜索方法之一，其基本操作是选择路径中的两个点，将这两个点之间的路径段反转。这种方法的作用和优点包括：

1. 减少路径长度：通过反转路径中的一段，可以减少路径的总长度，特别是当反转后的路径更短时。
2. 简单且高效：2 交换法的操作简单且计算量小，适合快速生成新解。
3. 避免局部最优：通过小范围的调整，可以避免陷入局部最优解，有助于进一步探索解空间。

3 交换法（3-opt）

3 交换法相比 2 交换法更复杂一些，其操作是选择路径中的三个点，然后重新排列这些点之间的路径段。具体的操作包括将其中一段路径插入到另一段路径之后。3 交换法的作用和优点包括：

1. 更大范围的调整：3 交换法可以进行更大范围的路径调整，相比 2 交换法能更有效地避开局部最优解。
2. 增加多样性：通过更复杂的路径调整，可以增加新解的多样性，有助于在更广阔的解空间中进行搜索。
3. 提高解的质量：由于 3 交换法涉及更多路径段的重新排列，因此有可能找到比 2 交换法更优的解，从而提高解的整体质量。
   2 交换法和 3 交换法在模拟退火算法中的作用是通过不同程度的路径调整，产生新的可能解，以此来搜索解空间中的最优解。可以在保证算法效率的同时，提高找到全局最优解的概率。通过邻域搜索方法的不断调整和降温过程的迭代，模拟退火算法能够逐步逼近全局最优解。

clc, clear, close all
sj0=load('data12_1.txt');
x=sj0(:,[1:2:8]); x=x(:);
y=sj0(:,[2:2:8]); y=y(:);
sj=[x y]; d1=[70,40]; 
xy=[d1;sj;d1]; sj=xy*pi/180; %角度化成弧度
d=zeros(102); %距离矩阵d初始化
for i=1:101
   for j=i+1:102
       d(i,j)=6370*acos(cos(sj(i,1)-sj(j,1))*cos(sj(i,2))*...
           cos(sj(j,2))+sin(sj(i,2))*sin(sj(j,2)));
   end
end
d=d+d';
path=[];long=inf; %巡航路径及长度初始化
for j=1:1000  %求较好的初始解
    path0=[1 1+randperm(100),102]; temp=0;
    for i=1:101
        temp=temp+d(path0(i),path0(i+1));
    end
    if temp<long
        path=path0; long=temp;
    end
end
e=0.1^30;L=20000;at=0.999;T=1;
for k=1:L  %退火过程
    c=2+floor(100*rand(1,2));  %产生新解
    c=sort(c); c1=c(1);c2=c(2);
    %计算代价函数值的增量
    df=d(path(c1-1),path(c2))+d(path(c1),path(c2+1))-...
        d(path(c1-1),path(c1))-d(path(c2),path(c2+1));
    if df<0 %接受准则
        path=[path(1:c1-1),path(c2:-1:c1),path(c2+1:102)]; long=long+df;
    elseif exp(-df/T)>=rand
        path=[path(1:c1-1),path(c2:-1:c1),path(c2+1:102)]; long=long+df;
    end
    T=T*at;
    if T<e
        break;
    end
end
path, long % 输出巡航路径及路径长度
xx=xy(path,1); yy=xy(path,2);
plot(xx,yy,'-*') %画出巡航路径

2.遗传算法

参考资料

现代优化算法 （二）： 遗传算法 及应用举例

算法简介

遗传算法（Genetic Algorithms，简称 GA）是一种基于自然选择原理和自然遗传机制的搜索（寻优）算法，它是模拟自然界中的生命进化机制，在人工系统中实现特定目标的优化。遗传算法的实质是通过群体搜索技术，根据适者生存的原则逐代进化，最终得到最优解或准最优解。
它必须做以下操作：初始群体的产生，求每一个体的适应度，根据适者生存的原则选择优良个体，被选出的优良个体两两配对，通过随机交叉其染色体的基因并随机变异某些染色体的基因生成下一代群体，按此方法使群体逐代进化，直到满足进化终止条件。
其实现方法如下：

1. 根据具体问题确定可行解域，确定一种编码方法，能用数值串或字符串表示可行解域的每一解。
2. 对每一解应有一个度量好坏的依据，它用一个函数表示，叫做适应度函数，一般由目标函数构成。
3. 确定进化参数群体规模 $M$，交叉概率 $p_c$、变异概率 $p_m$、进化终止条件。

为便于计算，一般来说，每一代群体的个体数目都取相等。群体规模越大，越容易找到最优解，但由于受到计算机的运算能力的限制，群体规模越大，计算所需要的时间也相应地增加。进化终止条件指的是当进化到什么时候结束，它可以设定到某一代进化结束，也可以根据找出近似最优解是否满足精度要求来确定。下表列出了生物遗传概念在遗传算法中的对应关系。

生物遗传概念	遗传算法中的作用
适者生存	算法停止时，最优目标值的可行解有最大的可能被保留住
个体	可行解
染色体	可行解的编码
基因	可行解中每一分量的特征
适应性	适应度函数值
种群	根据适应度函数值选取的一组可行解
交配	通过交配原则产生一组新的可行解的过程
变异	编码的某一分量发生变化的过程

遗传算法的详细步骤

用遗传算法继续探究应用举例中的问题。
求解的遗传算法的参数设定如下：
种群大小 $M=50$；最大代数 $G=100$；
交叉率 $p_c=1$，交叉概率为 1 能保证种群的充分进化；共有 $\frac{M}{2}$ 对相互配对的个体组；此时，被交叉的个体数目为 $M{p}_c=M$。
变异率 $p_m=0.1$，一般而言，变异发生的可能性较小。

编码策略

采用十进制编码，用随机数列 ${\omega }_1,{\omega }_2,\cdots ,{\omega }_{102}$ 作为染色体，其中 $0\le {\omega }_i\le 1$ $(i=2,3,\cdots ,101)$，${\omega }_1=0$，${\omega }_{102}=1$；每一个随机序列都和种群中的一个个体相对应，例如，目标问题的一个染色体为
$$[0.23,0.82,0.45,0.74,0.87,0.11,0.56,0.69,0.78]$$
其中编码位置 $i$ 代表目标，位置 $i$ 的随机数表示目标在巡回中的顺序，将这些随机数按升序排列得到如下巡回
$$6-1-3-7-8-4-9-2-5$$

初始种群

先利用经典的近似算法（改良圈算法）求得一个较好的初始种群。
对于随机产生的初始圈
$$C={\pi }_1,\cdots ,{\pi }_{u-1},{\pi }_u,{\pi }_{u+1},\cdots ,{\pi }_v,{\pi }_{v+1},\cdots ,{\pi }_{102}$$$$2\le u\le 101,2\le {\pi }_u\le 101$$
交换 $u$ 与 $v$ 之间的顺序，此时的新路径为
$${\pi }_1,\cdots ,{\pi }_{u-1},{\pi }_v,{\pi }_{v-1},\cdots ,{\pi }_u,{\pi }_{u+1},\cdots ,{\pi }_{102}$$
记
$$\Delta f=(d_{{\pi }_{u-1},{\pi }_u}+d_{{\pi }_v,{\pi }_{v+1}})-(d_{{\pi }_{u-1},{\pi }_v}+d_{{\pi }_u,{\pi }_{v+1}})$$
若 $\Delta f<0$，则以新路径修改旧路径，直到不能修改为止，就得到一个比较好的可行解。
直到产生 $M$ 个可行解，并把这 $M$ 个可行解转换成染色体编码。

目标函数

目标函数为侦察所有目标的路径长度，适应度函数就取为目标函数。我们要求
$$\min f({\pi }_1,{\pi }_2,\cdots ,{\pi }_{102})=\sum _{i=1}^{101}{d}_{{\pi }_i,{\pi }_{i+1}}$$

交叉操作

交叉操作采用单点交叉。设计如下，对于选定的两个父代个体 $f_1={\omega }_1{\omega }_2\cdots {\omega }_{102}$，$f_2={\omega }_1^{\prime }{\omega }_2^{\prime }\cdots {\omega }_{102}^{\prime }$，我们随机地选取第 $t$ 个基因处为交叉点，则经过交叉运算后得到的子代个体为 $s_1$ 和 $s_2$，$s_1$ 的基因由 $f_1$ 的前 $t$ 个基因和 $f_2$ 的后 $(102-t)$ 个基因构成， 的基因由 $f_2$ 的前 $t$ 个基因和 $f_1$ 的后 $(102-t)$ 个基因构成。
例如：
$$f_1=[0,0.14,0.25,0.27,\mid 0.29,0.54,\cdots ,0.19,1]$$$$f_2=[0,0.23,0.44,0.56,\mid 0.74,0.21,\cdots ,0.24,1]$$
设交叉点为第四个基因处，则
$$s_1=[0,0.14,0.25,0.27,\mid 0.74,0.21,\cdots ,0.24,1]$$$$s_2=[0,0.23,0.44,0.56,\mid 0.29,0.54,\cdots ,0.19,1]$$

变异操作

变异也是实现群体多样性的一种手段，同时也是全局寻优的保证。具体设计如下，按照给定的变异率，对选定变异的个体，随机地取三个整数 $u,v,w$，满足 $1<u<v<w<102$，把 $u,v$ 之间（包括 $u$ 和 $v$）的基因段插到 $w$ 后面。
变异操作（3 交换法）示意图：
$${\pi }_1\cdots {\pi }_{u-1}{\pi }_u{\pi }_{u+1}\cdots {\pi }_v{\pi }_{v+1}\cdots {\pi }_w{\pi }_{w+1}\cdots {\pi }_{102}$$

选择

采用确定性的选择策略，也就是说在父代种群和子代种群中选择目标函数数值最小的 $M$ 个个体进化到下一代，这样可以保证父代的优良特性被保留下来。
另外，选择策略也可选用轮盘轮选择、随机遍历抽样、锦标赛选择等方式。

结束准则

遗传算法的结束准则有以下三种：

1. 种群中个体的进化代数超过预定值；
2. 当达到目标函数的最优目标值时，可采用控制偏差的方式实现终止；
3. 在连续 $m$ 次进化后，种群适应度并没有变化，即可停止。
   

clc,clear, close all
sj0=load('data12_1.txt');
x=sj0(:,1:2:8); x=x(:);
y=sj0(:,2:2:8); y=y(:);
sj=[x y]; d1=[70,40]; 
xy=[d1;sj;d1]; sj=xy*pi/180;  %单位化成弧度
d=zeros(102); %距离矩阵d的初始值
for i=1:101
  for j=i+1:102
      d(i,j)=6370*acos(cos(sj(i,1)-sj(j,1))*cos(sj(i,2))*...
          cos(sj(j,2))+sin(sj(i,2))*sin(sj(j,2)));
  end
end
d=d+d'; w=50; g=100; %w为种群的个数，g为进化的代数
for k=1:w  %通过改良圈算法选取初始种群
    c=randperm(100); %产生1，...，100的一个全排列  
    c1=[1,c+1,102]; %生成初始解
    for t=1:102 %该层循环是修改圈 
        flag=0; %修改圈退出标志
    for m=1:100
      for n=m+2:101
        if d(c1(m),c1(n))+d(c1(m+1),c1(n+1))<...
                d(c1(m),c1(m+1))+d(c1(n),c1(n+1))
           c1(m+1:n)=c1(n:-1:m+1);  flag=1; %修改圈
        end
      end
    end
   if flag==0
      J(k,c1)=1:102; break %记录下较好的解并退出当前层循环
   end
   end
end
J(:,1)=0; J=J/102; %把整数序列转换成[0,1]区间上的实数，即转换成染色体编码
for k=1:g  %该层循环进行遗传算法的操作 
    A=J; %交配产生子代A的初始染色体
    c=randperm(w); %产生下面交叉操作的染色体对 
    for i=1:2:w  
        F=2+floor(100*rand(1)); %产生交叉操作的地址
        temp=A(c(i),[F:102]); %中间变量的保存值
        A(c(i),[F:102])=A(c(i+1),[F:102]); %交叉操作
        A(c(i+1),F:102)=temp;  
    end
    by=[];  %为了防止下面产生空地址，这里先初始化
while ~length(by)
    by=find(rand(1,w)<0.1); %产生变异操作的地址
end
B=A(by,:); %产生变异操作的初始染色体
for j=1:length(by)
   bw=sort(2+floor(100*rand(1,3)));  %产生变异操作的3个地址
   B(j,:)=B(j,[1:bw(1)-1,bw(2)+1:bw(3),bw(1):bw(2),bw(3)+1:102]); %交换位置
end
   G=[J;A;B]; %父代和子代种群合在一起
   [SG,ind1]=sort(G,2); %把染色体翻译成1，...,102的序列ind1
   num=size(G,1); long=zeros(1,num); %路径长度的初始值
   for j=1:num
       for i=1:101
           long(j)=long(j)+d(ind1(j,i),ind1(j,i+1)); %计算每条路径长度
       end
   end
     [slong,ind2]=sort(long); %对路径长度按照从小到大排序
     J=G(ind2(1:w),:); %精选前w个较短的路径对应的染色体
end
path=ind1(ind2(1),:), flong=slong(1)  %解的路径及路径长度
xx=xy(path,1);yy=xy(path,2);
plot(xx,yy,'-o') %画出路径

3.改进的遗传算法

参考资料

改进的遗传算法

引言

无人机航路规划问题实际上是一个组合优化问题，是优化理论中的 $NP-hard$问题。因为其解空间不连续，解邻域表达困难，所以难以用通常的算法求解。遗传算法作为现代优化算法之一，其主要特点是对非线性极值问题能以概率 1 跳出局部最优解，找到全局最优解。而遗传算法这种跳出局部最优寻找全局最优特性都基于算法中的交叉和变异。在传统遗传算法的结构中，变异操作在交叉操作基础上进行，强调的是交叉作用，认为变异只是一个生物学背景机制。
在具体交叉操作中，人们通常采用单点交叉（段交叉）、多点交叉与均匀交叉，其中单点交叉是指随机地在基因序列中选择一个断点，然后交换双亲断点右端的所有染色体。在变异操作中，变异算子一般是用 Guassian 分布的随机变量来实现。近年来，也有学者尝试用 Cauchy 分布的随机序列来实现变异，希望通过 Cauchy 分布宽大的两翼特性实现更大范围的变异，以利于找到全局最优解。Rudolph G 从理论上分析了采用 Cauchy 分布随机变异进化算法的局部收敛性。
Chellapilla K 进一步把二者结合起来，采用两种分布的线性叠加，但仿真结果显示，算法改进效果并不十分明显。文献将生物进化看成是随机性加上反馈，并指出其中的随机性主要是由系统的内在因素所引起，而不是由外部环境的随机扰动所造成。而混沌系统在其混沌域中表现为随机性，它是确定系统内部随机性的反映，不同于外在的随机特性。
根据以上特点对基于求解航路规划的遗传算法进行改进，首先将变异操作从交叉操作中分离出来，使其成为独立的并列于交叉的寻优操作，在具体遗传操作中，混沌与遗传操作联系在一起，在交叉操作中，以“门当户对”原则进行个体的配对，利用混沌序列确定交叉点，实行强度最弱的单点交叉，以确保算法收敛精度，前弱和避免寻优抖振问题；在变异操作中，利用混沌序列对染色体中多个基因进行变异以避免算法早熟。

算法改进

与标准的遗传算法相比，本节做了如下的两点改进：

交叉操作

本节的交叉操作采用改进型交叉。具体设计如下，首先以“门当户对”原则，对父代个体进行配对，即对父代以适应度函数（目标函数）值进行排序，目标函数值小的与小的配对，目标函数值大的与大的配对。然后利用混沌序列确定交叉点的位置，最后对确定的交叉项进行交叉。
例如: $({\Omega }_1,{\Omega }_2)$ 配对，它们的染色体分别是
$${\Omega }_1={\omega }_1^1{\omega }_2^1\cdots {\omega }_{102}^1,{\Omega }_2={\omega }_1^2{\omega }_2^2\cdots {\omega }_{102}^2$$
采用 Logistic 混沌序列 $x(n+1)=4x(n)(1-x(n))$ 产生一个 2 到 101 之间的正整数，具体步骤如下：

• 取一个 $(0,1)$ 区间上的随机数作为初始值，然后利用 $x(n+1)=4x(n)(1-x(n))$ 迭代一次产生 1 个 $(0,1)$ 区间上的混沌值，保存以上混沌值作为产生下一代交叉项的混沌迭代初值，再把这个值分别乘以 100 并加上 2，最后向下取整即可。
• 假如这个数为 $33$，那么以此做为交叉点对 $({\Omega }_1,{\Omega }_2)$ 染色体中相应的基因对进行单点交叉，得到新的染色体 $({\Omega }_1^{\prime },{\Omega }_2^{\prime })$
  $${\Omega }_1^{\prime }={\omega }_1^1{\omega }_2^1{\omega }_3^1{\omega }_4^1\cdots {\omega }_{33}^1{\omega }_{34}^2{\omega }_{35}^2\cdots {\omega }_{102}^2$$
  $${\Omega }_2^{\prime }={\omega }_1^2{\omega }_2^2{\omega }_3^2{\omega }_4^2\cdots {\omega }_{33}^2{\omega }_{34}^1{\omega }_{35}^1\cdots {\omega }_{102}^1$$
  很明显这种单点交叉对原来的解改动很小，这可以削弱避免遗传算法在组合优化应用中产生的寻优抖振问题，可以提高算法收敛精度。

变异操作

变异也是实现群体多样性的一种手段，是跳出局部最优、全局寻优的重要保证。这里变异算子设计如下，首先根据给定的变异率（本节选为 0.02），随机地取两个在 2 到 101 之间的整数，对这两个数对应位置的基因进行变异，变异时利用混沌序列把这两个位置的基因换成新的基因值，从而得到新的染色体。

性能比较

在仿真实验中，本节对航路规划问题分别利用单点交叉和换位变异结合的遗传算法、多点交叉和移位变异结合的遗传算法和本节中提出的改进算法进行求解比较。表 12.3 是各种算法种群规模（$M=50$）和迭代次数（$G=100$）相同的情况下连续 20 次求解的平均值（千米）、算法平均运算时间（秒）。

算法性能比较表

指标	单点交叉算法	多点交叉算法	改进的算法
平均航路距离（千米）	41572	40416	39849
算法执行时间（秒）	0.28	0.31	0.14

tic  %计时开始
clc,clear, close all
sj0=load('data12_1.txt'); 
x=sj0(:,1:2:8); x=x(:);
y=sj0(:,2:2:8); y=y(:);
sj=[x y]; d1=[70,40]; 
xy=[d1;sj;d1]; sj=xy*pi/180;  %单位化成弧度
d=zeros(102); %距离矩阵d的初始值
for i=1:101
  for j=i+1:102
      d(i,j)=6370*acos(cos(sj(i,1)-sj(j,1))*cos(sj(i,2))*...
             cos(sj(j,2))+sin(sj(i,2))*sin(sj(j,2)));
  end
end
d=d+d'; w=50; g=100; %w为种群的个数，g为进化的代数
rand('state',sum(clock)); %初始化随机数发生器
for k=1:w  %通过改良圈算法选取初始种群
    c=randperm(100); %产生1，...，100的一个全排列  
    c1=[1,c+1,102]; %生成初始解
    for t=1:102 %该层循环是修改圈 
        flag=0; %修改圈退出标志
    for m=1:100
      for n=m+2:101
        if d(c1(m),c1(n))+d(c1(m+1),c1(n+1))<...
                d(c1(m),c1(m+1))+d(c1(n),c1(n+1))
           c1(m+1:n)=c1(n:-1:m+1);  flag=1; %修改圈
        end
      end
    end
   if flag==0
      J(k,c1)=1:102; break %记录下较好的解并退出当前层循环
   end
   end
end
J(:,1)=0; J=J/102; %把整数序列转换成[0,1]区间上的实数，即转换成染色体编码
for k=1:g  %该层循环进行遗传算法的操作 
    A=J; %交配产生子代A的初始染色体
    for i=1:2:w
        ch1(1)=rand; %混沌序列的初始值
        for j=2:50
            ch1(j)=4*ch1(j-1)*(1-ch1(j-1)); %产生混沌序列
        end
        ch1=2+floor(100*ch1); %产生交叉操作的地址
        temp=A(i,ch1); %中间变量的保存值
        A(i,ch1)=A(i+1,ch1); %交叉操作
        A(i+1,ch1)=temp;
    end
    by=[];  %为了防止下面产生空地址，这里先初始化
while ~length(by)
    by=find(rand(1,w)<0.1); %产生变异操作的地址
end
num1=length(by); B=J(by,:); %产生变异操作的初始染色体
ch2=rand;  %产生混沌序列的初始值
for t=2:2*num1 
       ch2(t)=4*ch2(t-1)*(1-ch2(t-1)); %产生混沌序列
end
for j=1:num1
   bw=sort(2+floor(100*rand(1,2)));  %产生变异操作的2个地址
   B(j,bw)=ch2([j,j+1]); %bw处的两个基因发生了变异
end
   G=[J;A;B]; %父代和子代种群合在一起
   [SG,ind1]=sort(G,2); %把染色体翻译成1，...,102的序列ind1
   num2=size(G,1); long=zeros(1,num2); %路径长度的初始值
   for j=1:num2
       for i=1:101
           long(j)=long(j)+d(ind1(j,i),ind1(j,i+1)); %计算每条路径长度
       end
   end
     [slong,ind2]=sort(long); %对路径长度按照从小到大排序
     J=G(ind2(1:w),:); %精选前w个较短的路径对应的染色体
end
path=ind1(ind2(1),:), flong=slong(1)  %解的路径及路径长度
toc  %计时结束
xx=xy(path,1);yy=xy(path,2);
plot(xx,yy,'-o') %画出路径

4.Matlab 遗传算法函数

参考资料

matlab 遗传算法（GA）详解（一）算法入门
matlab 官方介绍

使用规则

在每一步中，遗传算法随机地从当前种群中选择若干个体作为父辈，并且使用它们产生下一代的子种群。在连续若干代之后，种群朝着优化解的方向进化。可以用遗传算法来求解各种不适宜于用标准算法求解的优化问题，包括目标函数数不连续、不可微、随机或高度非线性的问题。
遗传算法在每一步使用下列三类规则从当前种群来创建下一代：

1. 选择规则 (Selection Rules): 选择对下一代种群有贡献的个体（称为父辈）。
2. 交叉规则 (Crossover Rules): 将两个父辈结合起来构成下一代的子辈种群。
3. 变异规则 (Mutation Rules): 施加随机变化用父辈个体来构成子辈。

2. 使用方式

1. 以命令行方式调用遗传算法函数 $ga$。
2. 通过用户图形界面使用遗传算法工具。

遗传算法函数

遗传算法函数 ga 调用格式如下：

[x, fval] = ga(fun, nvars, A, b, Aeq, beq, LB, UB, nonlcon, options)

返回值 x 是最终值到达的点，它为所求函数的局部极小点，这里 x 为行向量，fval 为目标函数的极小值。

示例

求下列问题的解：
$$\begin{gathered} \max f(x)=2x_1+3x_1^2+3x_2+x_2^2+x_3, \\ \text{s.t.}\{\begin{array}{ll}& x_1+2x_1^2+x_2+2x_2^2+x_3\le 10,\\ & x_1+x_1^2+x_2+x_2^2-x_3\le 50,\\ & 2x_1+x_1^2+2x_2+x_3\le 40,\\ & x_1^2+x_3=2,\\ & x_1+2x_2\ge 1,\\ & x_1\ge 0,x_2,x_3\text{不约束}.\end{array} \end{gathered}$$

clc, clear
a=[-1 -2 0;-1 0 0];b=[-1;0]; %两个线性约束条件
[x,y]=ga(@obj,3,a,b,[],[],[],[],@constr);  %%调用遗传算法函数
x, y=-y

%再单独创建一个.m后缀的文件命名为obj.m
function y=obj(x);   %定义目标函数，其中x为行向量
c1 = [2 3 1]; c2 = [3 1 0];
y = c1*x' + c2*x'.^2; 
y=-y;   %ga中目标函数极小化
end

%再单独创建一个.m后缀的文件命名为constr.m
function [f,g]=constr(x);  %定义非线性约束函数
f=[x(1)+2*x(1)^2+x(2)+2*x(2)^2+x(3)-10
   x(1)+x(1)^2+x(2)+x(2)^2-x(3)-50
   2*x(1)+x(1)^2+2*x(2)+x(3)-40];
g=x(1)^2+x(3)-2;
end

可以如上将对应两个函数分别建立为文件再调用,也可以如下放在同一个文件夹

注意：遗传算法程序的运行结果每一次都是不一样的，要运行多次，找一个最好的结果。

Matlab 遗传算法工具

遗传算法与直接搜索概述

Matlab 中遗传算法与直接搜索（Genetic Algorithm and Direct Search，简称 GADS）工具箱是一系列函数的集合，它们扩展了优化工具箱和 Matlab 数值计算环境。遗传算法与直接搜索工具箱包含了要使用遗传算法和直接搜索算法来求解优化问题的一些例程。这些算法使我们能够求解那些标准优化工具箱范围之外的各种优化问题。
工具箱函数可以通过图形界面或 Matlab 命令行来访问，它们是用 Matlab 语言编写的，对用户开放，因此可以查看算法，修改源代码（例如 >>edit ga）或生成用户函数。
遗传算法与直接搜索工具箱有助于求解那些不易用传统方法解决的问题，譬如旅行商问题等。
遗传算法与直接搜索工具箱有一个精心设计的用户图形界面，可直观、方便、快速地求解优化问题。

功能特点

遗传算法与直接搜索工具箱的功能特点如下：

1. 用户图形界面和命令行函数可用来快速地描述问题、设置算法选项以及监控进程。
2. 具有多个选项的遗传算法工具可用于问题创建、适应度计算、选择、交叉和变异。
3. 直接搜索工具实现了一种模式搜索方法，其选项可用于定义网格尺寸、表决方法和搜索方法。
4. 遗传算法与直接搜索工具箱函数可与 Matlab 的优化工具箱或其它的 Matlab 程序结合使用。
5. 支持自动的 M 代码生成。

2.2.3 通过 GUI 使用遗传算法

遗传算法工具有一个用户图形界面 GUI，它使我们可以使用遗传算法而不必在命令行方式工作。遗传算法的用户图形界面集成在优化工具箱里，打开遗传算法用户图形界面，可键入以下命令：>>optimtool （根据 Matlab 的不同版本, 也可用命令：>>gatool）

2.2.4 直接搜索工具

直接搜索的命令为：

[x,fval] = patternsearch(@fun, x0, A, b, Aeq, beq, LB, UB, @nonlcon, options)

直接搜索的用户图形界面也集成到了 optimtool 中。

function ex12_42023
a=[-1 -2 0;-1 0 0];b=[-1;0]; %两个线性约束条件
[x,y]=patternsearch(@obj,rand(1,3),a,b,[],[],[],[],@constr);  %%直接搜索算法
x, y=-y

function y=obj(x);   %定义目标函数，其中x为行向量
c1 = [2 3 1]; c2 = [3 1 0];
y = c1*x' + c2*x'.^2; 
y=-y;   %ga中目标函数极小化

function [f,g]=constr(x);  %定义非线性约束函数
f=[x(1)+2*x(1)^2+x(2)+2*x(2)^2+x(3)-10
   x(1)+x(1)^2+x(2)+x(2)^2-x(3)-50
   2*x(1)+x(1)^2+2*x(2)+x(3)-40];
g=x(1)^2+x(3)-2;

三、相关经验

1. matlab 中的代码删除分号可以将结果显示出来,这样可以帮助我们理解别人代码中运行步骤
2. 启发式算法程序的运行结果每一次都是不一样的，可以运行多次，找一个最好的结果。