12现代优化算法

现代优化算法是 20世纪 80年代初兴起的启发式算法，包括禁忌搜索（Tabu Search）、模拟退火（Simuated Annealing）、遗传算法（Cenetic Algorithnas）、人工神经网络（Neural Networks）。
它们主要用于解决那些传统优化方法（如梯度下降法、牛顿法、单纯形法等）难以处理或效率低下的复杂问题。这些问题通常具有高维度、非线性、非凸性、多峰性、不连续性、存在大量局部最优解，或者目标函数难以求导等特点。

12.1 模拟退火算法

模拟退火 (Simulated Annealing)。它来源于冶金学中对金属进行“退火”处理的过程。

1. 加热 (High Temperature)：将金属加热到很高的温度，此时金属内部的粒子（原子）能量非常高，可以自由、无序地运动。
2. 缓慢冷却 (Slow Cooling)：然后，非常缓慢地降低温度。在这个过程中，粒子能量逐渐降低，它们会自发地寻找能量最低、最稳定的位置，最终形成规整、坚固的晶体结构。
3. 结果 (Low Energy State)：最终得到的金属，其内部处于一种能量最低的、最稳定的状态。

---

在算法中，问题的解 (Solution)对应粒子的状态 (State)；解的优劣程度 (Cost/Objective Function)对应能量 (Energy)。一个“好”的解，其“能量”就低，算法目标就是找到能量最低的解；算法的控制参数 (Control Parameter)对应 温度 (Temperature)。

12.1.1 算法介绍

1 状态转移的规则 (Metropolis准则)

假设材料在状态 $j$ 之下的能量为 $E(i)$ , 那么材料在温度 $T$ 时从状态 $i$ 进人状态 $j$ 就遵循如下规律:
（1）若 $E(j) \leq E(i)$ 则接受该状态被转换；
（2）若 $E(j) \leq E(i)$ 则状态转换以下面概率被接受：

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

其中：$K$ 为物理学中的玻耳兹曼常数，$T$ 为材料温度。

算法理解：
（1）接受“更好的解”：
若新解 $j$ 优于或与当前解 $i$ 相当（即 $E(j) \le E(i)$），总接受这个新解。
（2）有条件接受“更差的解”（跳出局部最优解）：
若新解 $j$ 劣于当前解 $i$ （即 $E(j) > E(i)$），不直接拒绝，以特定概率接受。
这个概率是 $P = e^{-\frac{\Delta E}{kT}}$，其中 $\Delta E = E(j) - E(i)$ 是新旧解的能量差（一个正数）。

关于 $P$ ：

• 与能量差 $\Delta E$ 的关系：能量差 $\Delta E$ 越大（新解差得越多），指数的绝对值就越大，概率 $P$ 就越小。这意味着，我们更倾向于接受一个“只差一点点”的解，而不愿意接受一个“差得离谱”的解。
• 与温度 $T$ 的关系：
  • 当温度 $T$ 很高时，分母 $kT$ 很大，$\frac{\Delta E}{kT}$ 接近于0，所以概率 $P$ 接近于 $e^0=1$。这意味着在高温时，算法非常“大胆”，几乎愿意接受任何新解，无论好坏。这对应了物理退火中粒子自由运动的阶段，目的是充分探索整个解空间。
  • 当温度 $T$ 很低时，分母 $kT$ 很小，$\frac{\Delta E}{kT}$ 是一个很大的正数，所以概率 $P$ 接近于0。这意味着在低温时，算法变得非常“保守”，几乎不再接受更差的解，只会向更好的方向移动。这对应了物理退火中粒子能量降低、逐渐稳定下来的阶段，目的是收敛到最优解。

2 热平衡状态与玻耳兹曼分布

1. 提议概率 $g(i, j)$ ：从当前状态 $i$，根据某种规则生成一个候选的新状态 $j$。
2. 接受概率 $A(i, j)$ ：根据Metropolis准则，计算一个接受概率来决定是否接受这个新状态 $j$。

   $$$$

   $$$$

   状态转移概率 $P(i \to j)$ 的分段定义：

• 如果 $E(j) \leq E(i)$，则 $P(i \to j) = 1$
• 如果 $E(j) > E(i)$，则 $P(i \to j) = e^{\frac{E(i) - E(j)}{kT}}$

现仅考虑$E(i) > E(j)$ 的情况：
在热平衡状态下，任意两个状态 $i$ 和 $j$ 的转移应该是平衡的，根据状态转移的原则，平衡条件可以写为：

$$P_T(i) \cdot P(i \to j) = P_T(j) \cdot P(j \to i) \implies \frac{P_T(i)}{P_T(j)} = \frac{P(j \to i)}{P(i \to j)}$$

得到：

$$\frac{P_T(i)}{P_T(j)} = \frac{1}{e^{\frac{E(i) - E(j)}{kT}}} = e^{-\frac{E(i) - E(j)}{kT}} = e^{\frac{E(j) - E(i)}{kT}}$$

要得到具体的概率分布，需满足归一化条件：

$$\sum_{j \in S} P_T(X = j) = 1$$

假设 $P_T(X = i) = C e^{-\frac{E(i)}{kT}}$，其中 $C$ 是一个常数，代入归一化条件得到：

$$\sum_{j \in S} C e^{-\frac{E(j)}{kT}} = 1$$

从而得到：

$$C = \frac{1}{\sum\limits_{j \in S} e^{-\frac{E(j)}{kT}}}$$

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

$$P_T(X = i) = \frac{e^{-\frac{E(i)}{kT}}}{\sum\limits_{j \in S} e^{-\frac{E(j)}{kT}}}$$

其中：$X$ 为材料当前状态的随机变量; $S$ 为状态空间集合。

算法理解：
如果我们让温度 $T$ 固定不变，并让算法运行足够长的时间，系统会达到一个“稳定状态”（热平衡）。在这个稳定状态下，算法并不是停留在某一个解上，而是在不同的解之间来回跳转。
玻耳兹曼分布描述了在温度 $T$ 下，算法停留在任意一个解 $i$ 的概率，它的核心是 $e^{-\frac{E(i)}{kT}}$ 。

• 能量 $E(i)$ 越低，指数的负值越小， $e^{-\frac{E(i)}{kT}}$ 的值就越大。
• 能量 $E(i)$ 越高，指数的负值越大，$e^{-\frac{E(i)}{kT}}$ 的值就越小。
  结论是：在任何给定的温度下，系统总是更倾向于处于能量较低的状态（更好的解）。

3 极端温度下的系统行为

这部分是从数学上证明，模拟退火算法的“高温探索”和“低温收敛”这两个行为是合理的。
（1）当温度 $T \to \infty$ (高温情况)
当温度 $T$ 趋向于无穷大时，$\frac{E(i)}{kT}$ 的值会趋向于零。这样，指数项会变为：

$$e^{-\frac{E(i)}{kT}} \to e^0 = 1$$

同样地，所有状态的指数项也会趋向于1：

$$e^{-\frac{E(j)}{kT}} \to 1 \quad \text{对于所有 } j \in S$$

因此：

$$\sum\limits_{j \in S} e^{-\frac{E(j)}{kT}} \to \sum\limits_{j \in S} 1 = |S|$$

其中：$|S|$ 为集合 $S$ 中状态的数量。
代入原来的概率表达式得到：

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

算法理解：
在极高的温度下，所有解（无论好坏）被访问到的概率都变得几乎一样。这正是我们想要的“全局探索”阶段。在算法开始时，设置一个很高的初始温度，让算法可以不受能量（解的优劣）的束缚，在整个解空间中自由地、随机地漫游，避免过早地陷入某个局部最优区域。

（2） 当温度 $T \to 0$ (低温情况)
在数学上，我们无法直接从 $\frac{0}{0}$ 得出极限值，因此将分子和分母上同乘 $e^{\frac{E_{\min}}{kT}}$，转化为：

$$\lim_{T \to 0} \frac{e^{-\frac{E(i) - E_{\min}}{kT}}}{\sum\limits_{j \in S} e^{-\frac{E(j) - E_{\min}}{kT}}}$$

$$= \lim_{T \to 0} \frac{e^{-\frac{E(i) - E_{\min}}{kT}}} {\sum\limits_{j \in S_{\min}} e^{-\frac{E(j) - E_{\min}}{kT}} + \sum\limits_{j \notin S_{\min}} e^{-\frac{E(j) - E_{\min}}{kT}}}$$

$$= \lim_{T \to 0} \frac{e^{-\frac{E(i) - E_{\min}}{kT}}}{\sum\limits_{j \in S_{\min}} e^{-\frac{E(j) - E_{\min}}{kT}}} = \begin{cases} \frac{1}{|S_{\min}|}, & i \in S_{\min} \\\\ 0, & else\end{cases}$$

其中：$E_{\min} = \min E(j)$ 且 $S_{\min} = \{ i \mid E(i) = E_{\min} \}$.
当温度 $T$ 趋近于0时，能量的微小差异都会被急剧放大。

• 对于能量最低的解（$i \in S_{\min}$），$E(i) - E_{\min} = 0$，所以分子中的项 $e^{-\frac{E(i)-E_{\min}}{kT}}$ 变成了 $e^0=1$。
• 对于任何能量不是最低的解（$j \notin S_{\min}$），$E(j) - E_{\min}$ 是一个正数。当 $T \to 0$ 时，$\frac{E(j) - E_{\min}}{kT}$ 趋向于正无穷大，所以 $e^{-\frac{E(j) - E_{\min}}{kT}}$ 趋近于0。
  算法理解：
  当温度降到极低时，只有那些能量最低的解（全局最优解）出现的概率不为零，而其他所有解的概率都趋近于0。如果有多个全局最优解，它们会平分这个概率。这正是“收敛和利用”阶段，随着温度的缓慢降低，算法接受“差解”的概率越来越小，最终基本只接受“好解”，从而稳定地收敛到它所找到的全局最优解上。

4 优化问题的应用

考虑这样一个组合优化问题：优化函数为 $f: x \to \mathbb{R}^+$，其中 $x \in S$，它表示优化问题的一个可行解，$\mathbb{R}^+ = \{ y \;|\; y \in \mathbb{R}, y \geq 0 \}$，$S$ 表示函数的定义域。$N(x) \subseteq S$ 表示 $x$ 的一个邻域集合。
先给定一个初始温度 $T_0$ 和该优化问题的一个初始解 $x(0)$，并由 $x(0)$ 生成下一个解 $x' \in N(x(0))$，是否接受 $x'$ 作为一个新解 $x(1)$ 依赖于如下概率：

$$P(x(0) \to x') = \begin{cases} 1, & f(x') < f(x(0)) \\ \frac{e^{\frac{f(x') - f(x(0))}{T_0}}}{Z}, & else\\ \end{cases}$$

在温度 $T_i$ 下，经过很多次的转移后降低温度 $T_i$ 得到 $T_{i+1}$ ＜ $T_i$ ，在 $T_{i+1}$ 下重复上述过程。因此整个优化过程就是不断寻找新解和缓慢降温的交替过程，最终的解是对该问题寻优的结果。

马尔可夫过程：

Wiki：马尔可夫性质是概率论中的一个概念。当一个随机过程在给定现在状态及所有过去状态情况下，其未来状态的条件概率分布仅依赖于当前状态；换句话说，在给定现在状态时，它与过去状态（即该过程的历史路径）是条件独立的，那么此随机过程即具有马尔可夫性质。具有马尔可夫性质的过程通常称之为"马尔可夫过程")。马尔可夫过程不具备记忆特质（memorylessness）。

注意到在每个 $T_i$ 下所得到的一个新状态 $x(k+1)$ 完全依赖于前一个状态 $x(k)$ ，和前面的状态 $x(0),x(1),…,x(k-1)$ 无关，因此这是一个马尔可夫过程。
因此任何一个状态从 $x(k)$ 生成 $x(k+1)$ 的概率，在 $N(x_k)$ 中是均匀分布的，且新状态被接受的概率满足上式，那么经过有限次转换，在温度 $T_i$ 下的平衡态 $x_i$ 的分布由下式给出:

$$P_i(T_i) = \frac{e^{-\frac{f(x_i)}{T_i}}}{\sum\limits_{j \in S} e^{-\frac{f(x_j)}{T_j}}}$$

当温度 $T$ 降为 0 时，$x_i$ 分布为：

$$P_i^* = \begin{cases} \frac{1}{|S_{min}|}, & x_i \in S_{min}\\ \\ 0, & else\\ \end{cases}$$

且

$$\sum_{x_i \in S_{min}} P_i^* = 1$$

这说明如果温度下降十分缓慢，而在每个温度都有足够多次的状态转移，使之在每一个温度下达到热平衡，则全局最优解将以概率 1 被找到。因此可以说模拟退火算法可以找到全局最优解。

注意事项：

1. 理论上，降温过程要足够缓慢，要使得在每一温度下达到热平衡。计算机实现中用时要综合考虑解的性能和算法速度，在两者之间折中。
2. 要确定在每一温度下状态转换的结束准则。实际操作可以考虑当连续 m 次的转换过程没有使状态发生变化时结束该温度下的状态转换。最终温度的确定可以提前定为一个较小的值 $T_e$ 或连续几个温度下转换过程没有使状态发生变化算法就结束。
3. 选择初始温度和确定某个可行解的邻域的方法也要恰当。

12.1.2 旅行商问题

基地经度和纬度为 (70 ’ 40) ，假设飞机速度为 1000 km/h 。派一架飞机从基地出发，侦察完所有目标，再返回原来基地。已知 100 个目标的经纬度，在每个目标点的侦察时间不计，求飞机所花费的时间（假设飞机巡航时间可以充分长）。
首先进行编号，将基地编号为 0 和 101 (返航)，目标点分别标号 1~100 。利用距离矩阵 $D(d_{ij})_{102*102}$ 表示，注意 $D$ 为实对称矩阵。问题将转化为求 0 到 101 的最短路径。
需将经纬度坐标转化为两点间的实际距离。设 $A,B$ 两点的地理坐标分别为$A(x_1, y_1), B(x_2, y_2)$ ，则两点间的实际距离应是过两点大圆的劣弧长。

Wiki：优弧是指圆心位于弧与弦连接成的封闭图形之内；劣弧是指圆心位于弧与弦连接成的封闭图形之外。

以地心为坐标原点 $O$ ，以赤道平面为 $XOY$ 平面，以0度经线圈所在的平面为 $XOZ$ 平面建立三维直角坐标系，则两点直角坐标分别为：

$$\begin{cases} 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)\\ \end{cases}$$

因此：

$$d = R\;arccos\left(\frac{OA \cdot OB}{|OA|\cdot|OB|}\right)$$

化简得到：

$$d = R\;arccos \,[\,cos\,(x_1-x_2)\,cos\,y_1\,cos\,y_2\,+\,sin\,y_1\,sin\,y_2]$$

其中：地球半径 $R=6370$ km。

1 解空间

有固定起点为0，终点为101，中间是1~100的循环排列集合，表示为：

$$S = \{(\pi_0, \pi_1, ..., \pi_{101})\;|\;\pi_0 = 0, (\pi_1, ..., \pi_{100})\text{是{1,..., 100}的循环排列}, \pi_{101} = 101 \}$$

可先使用蒙特卡洛方法求得一个较好的初始解。

2 目标函数

目标函数（或称代价函数）为侦察所有目标的路径长度：

$$min \quad f(\pi_0, \pi_1, ..., \pi_{101}) = \sum\limits_{i=0}^{100} d_{\pi_i\pi_{i+1}}$$

3 迭代

（1）新解的产生
设迭代的解为：

$$\pi_0 ... \pi_{u-1} \pi_{u} \pi_{u+1}...\pi_{v-1} \pi_{v} \pi_{v+1}...\pi_{w-1}\pi_{w}\pi_{w+1} ... \pi_{100}$$

• 2变换法
  任选序号 $u, v$ ，交换 $u$ 与 $v$ 之间的顺序，变成逆序。

  $$\pi_0 ... \pi_{u-1} \pi_{v} \pi_{v-1}...\pi_{u+1} \pi_{u} \pi_{v+1}...\pi_{w-1}\pi_{w}\pi_{w+1} ... \pi_{100}$$

• 3变换法
  任选序号 $u, v$ 和 $w$ ，将 $u$ 与 $v$ 之间的路径插到 $w$ 后，变成逆序。

  $$\pi_0 ... \pi_{u-1}\pi_{v+1}...\pi_{w-1}\pi_{w} \pi_{u} \pi_{u+1}...\pi_{v-1} \pi_{v} \pi_{w+1} ... \pi_{100}$$

  （2）代价函数差
• 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}} )$$

  （3）接受准则

  $$P = \begin{cases} 1, & \Delta f < 0, \\ \exp\left(-\frac{\Delta f}{T}\right), & \Delta f \geq 0. \end{cases}$$

4 降温

利用选定的降温系数 $\alpha$ 进行降温，取新的温度 $T$ 为 $\alpha T$ ，这里选定 $\alpha = 0.999$ 。

5 结束

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

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);

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(c1), 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, '-*') % 画出巡航路径

12.2 遗传算法

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

实现方法

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

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

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

求解的遗传算法的参数设定如下：种群大小 $M=50$ ，最大代数 $G=1000$ ，交叉率 $p_c=1$ (保证种群的充分进化)，变异率 $p=0.1$ (一般而言变异发生的可能性较小)。

12.2 改进的遗传算法

遗传算法作为现代优化算法之一，其主要特点是对非线性极值问题能以概率1跳出局部最优解，找到全局最优解。而遗传算法这种跳出局部最优、寻找全局最优的特性都基于算法中的交叉和变异。

在传统遗传算法的结构中，变异操作在交叉操作基础上进行，强调的是交叉作用，认为变异只是一个生物学背景机制。在具体交叉操作中，人们通常采用单点交叉（段交叉）、多点交叉与均匀交叉 其中单点交叉是猪随机地在基因序列中选操一个断点，然后交换双亲上断点右端的所有染色体。

首先以"门当户对"原则对父代个体进行配对，即对父代以适应度函数（目标函数）值进行排序，目标函数值小的与小的配对，目标函数值大的与大的配对。然后利用混沌序列确定交叉点的位置，最后对确定的交叉项进行交叉。