进化算法（遗传算法、粒子群算法、差分算法）—— Rosenbrock函数优化

进化算法

进化算法是一类算法簇，这类算法模拟的是自然界当中生物种群的进化机制。进化算法在解决最优化问题或者机器学习问题的时候，通常会维护一个解的集合X={x1,x2,…,xN}X=\{x_1,x_2,\dots,x_N\}X={x1​,x2​,…,xN​}，称之为种群，$N$称为种群规模。算法运行的机制，是在每一轮迭代当中，对种群进行变异、交叉等操作，形成新的个体添加到种群当中，然后对种群进行筛选，使得种群的规模始终维持在$N$. 筛选的规则是根据种群个体的适应度$f(x_i)$设计的，理想的筛选规则应当保留适应度高的个体，淘汰适应度低的个体。

1. 三种算法的标准式表述

我们这里介绍的进化算法，主要是三种：遗传算法、粒子群算法和标准差分算法。

首先我们先定义基本的问题背景：在一个环境$\text{env}$（可能是一个最优化问题、一个强化学习场景、一个实际的工程模拟）当中，每一个解都表示为向量$x_i\in {\mathbb{R}}^d$. 存在一个适应度函数$f(x_i)$能够评价每一个解的适应性，函数评分越高，说明解越适应于该环境。

我们接下来分别介绍三种算法解决该问题的流程。

1.1 遗传算法

遗传算法模拟的是自然界中生物种群的繁衍以及自然选择机制。一个遗传算法的流程可以表示为：

1. 初始化种群$X$，设置进化代序为$t=0$
2. 计算种群的适应度，并且利用选择算子筛选进行遗传的个体$X_c=\text{select}(X)$
3. 对进行遗传的个体进行交叉算子操作$X_m=\text{crossover}(X_c)$
4. 对全体种群进行变异操作$X^{\prime }=\text{mutation}(X_m)$
5. 更新种群$X=\text{Reproduction}(X^{\prime })$
6. 检查是否满足终止条件，若不满足，更新进化代序$t=t+1$并返回第2步

需要注意的是，遗传算法同样是一个算法簇，而不是某一具体算法。对于具体的实现，筛选、交叉、变异算子在算法流程中的顺序，以及它们对当前种群的更新机制会有一定的差异。

在这个流程当中，最为关键的是三个遗传算子：$\text{select}(\cdot ),\text{crossover}(\cdot ),\text{mutation}(\cdot )$，还有一个非关键的更新算子$\text{Reproduction}$.

1.1.1 比例选择算子

假设每一个个体的适应度$f(x_i)$都是非负的，那么进行变换

$$p_i=\frac{f(x_i)}{{\sum }_{j=1}^{N}f(x_j)}$$

在选择过程中，每个个体被选中加入交叉种群$X_c$的概率是$p_i$. 实现这种选择的方法，可以采用轮盘抽奖的算法。假设我们在面积为$1$的抽奖圆盘上，为每一个个体分配面积为$p_i$的扇面，然后转动圆盘，当圆盘停止时，指针指向第$i$个个体的扇面的概率恰好是$p_i$. 如图：

从计算机程序实现角度来说，我们可以规定一个起始坐标$q_0=0$，然后用累加概率的方式找出剩下的各个分界点坐标

$$q_i=\sum _{j=1}^i{p}_i$$

随机取一个$[q_0,q_N]$之间的点$\lambda$，如果$\lambda \in [q_{i-1},q_i]$，那就意味着指针落在了$x_i$对应的圆盘扇面内，$x_i$被选中加入$X_c$。把这个过程重复$N$次，就会选择出期望为$N{p}_c$个的种群个体。

这种采样过程产生的交叉种群数$N_c$是不固定的，而且比原始种群数$N$要少，需要一定的机制，例如将交叉运算后的个体添加回原种群当中，才能保证种群的数目不会随着迭代次数增加而减少。

1.1.2 其他选择算子

最优选择算子：这种算子通常指定一个数量$c<N$，在每一轮迭代当中，它用当前种群里适应度最高的$c$个个体，替换掉适应度最低的$c$个个体，然后直接用得到的新种群作为$X_c$。更新后的种群继续进行交叉、变异算子的操作。

最优保存算子：这种算子不对种群进行筛选，$X_c=X$. 只是记录下当前适应度最高的个体。最后在进行交叉、变异结束之后的种群$X^{\prime }$中，将最优个体插入其中。使得最优个体保持到下一轮迭代。

生存数目算子：计算每一个个体适应度与平均适应度的比值，并且下取整：

$$N_i=\lfloor \frac{Nf(x_i)}{{\sum }_{j=1}^{N}f(x_j)}\rfloor$$

将原种群当中的每一个个体$x_i$，复制$N_i$次，加入到交叉种群$X_c$当中。按照适应度对于$X_c$排序，截取其前$N$个个体。

1.1.3 交叉算子

一般来说，交叉算子总是规定了交叉概率$P_c$，并且希望交叉运算产生的子代数目期望为$N_c{P}_c$. 交叉算子有一个重要的问题是，交叉的两个父本如何选择。这个问题没有固定解法。有的算法当中，会将交叉种群$X_c$打乱顺序，之后按顺序两两匹配，对于每一个匹配都以$P_c$的概率进行交叉或不交叉；有的算法会随机抽取$n$轮，每轮抽取一对父本以$P_c$的概率进行交叉或不交叉；有的算法会随机产生一个$(0,1)$之间的随机数$\lambda$，如果$\lambda <P_c$，就对整个$N_c$进行打乱、按顺序两两匹配、交叉运算，否则就完全不进行交叉运算。

交叉算子最常用的是单点交叉算子，即随机地在个体向量上取一个点位$\text{min}$. 对于配对的两个向量$x_i,x_j$，将它们从点位$\text{min}$开始直到结束的片段交换：

$$\begin{gathered} {\overline{x}}_{i,n}=x_{j,n},n\ge \text{min} \\ {\overline{x}}_{j,n}=x_{i,n},n\ge \text{min} \end{gathered}$$

其次还有双点位交叉，即选择两个点位$\text{min},\text{max}$，交换父本在这两个点位之间的片段。

有的时候，还会使用凸组合交叉算子，即随机采样一个$a\in (0,1)$，令

$$\begin{gathered} {\overline{x}}_i=a{x}_i+(1-a)x_j \\ {\overline{x}}_j=(1-a)x_i+a{x}_j \end{gathered}$$

这种算子适用于解空间是高维凸域的情形。

最后一个问题是，交叉算子产生的子代和父代如何进行取舍。有的算法是将子代添加到种群$X_c$中组合形成$X^{\prime }$，然后通过更新算子$\text{Reproduction}(X^{\prime })$进行取舍（例如选择出适应度前$N$个个体）。有的算法会比较子代与父代的适应度，然后选择保留适应度强的子代。

1.1.4 变异算子

变异算子同样指定一个变异概率$P_m$，并且试图使变异的个体数期望为$N_m{P}_m$. 而实现的方法，可以是对于$X_m$的每一个个体都进行二分决策，以$P_m$的概率进行突变，也可以是对整个$X_m$进行二分决策，突变或者不突变。

对于每一个突变的个体，如果是二进制串（二进制的向量），可以直接进行位反转。

对于一般的向量，还可以随机选取一个片段进行倒置（顺序颠倒），或者选择两个点位互换数值。

和交叉算子类似，变异算子也可以采用线性组合的方式进行变异：随机选取一个$v\in {\mathbb{R}}^d$和一个$\lambda \in (0,1)$，令

$${\hat{x}}_i=x_i+\lambda v$$

1.1.5 更新算子

更新算子的实现很简单，如果选择、交叉、变异算子最终产生的种群$X^{\prime }$规模是$N$，可以直接将其更新到$X$上. 如果数目超过$N$，那么可以选择保留适应度最高的$N$个个体。（也就是需要再次遍历并且进行适应度评价）

1.2 粒子群算法

粒子群算法和遗传算法的区别在于，粒子群算法没有交叉、变异的相关概念，而是会在每一轮迭代当中，更新每一个个体的值。在粒子群算法当中，种群当中的个体称为粒子，由两个$d$维的向量$x_i\in {\mathbb{R}}^d,v_i\in {\mathbb{R}}^d$表示，$x_i$称为粒子的位置，$v_i$称为粒子的速度。一个种群可以表示为$E=(X,V)$，其中$X=(x_1,x_2,\dots ,x_N)$，$V=(v_1,v_2,\dots ,v_N)$. 环境中的适应性函数仅仅取决于位置$f(x_i)$.

粒子群算法当中，每一个粒子都需要维护两个变量$p_i\in \mathbb{R},P_i\in {\mathbb{R}}^d$，其中$p_i$是在以往迭代轮数当中，粒子$i$所取得的最大适应度，而$P_i$是取得该适应度时的位置。

同时在每一轮迭代当中，还要找出两个变量$s\in \mathbb{R},S\in {\mathbb{R}}^d$，分别是当前种群当中，适应度的最大值，以及适应度最大的粒子所处位置。

标准粒子群算法的参数包括：

1. 最大进化代数$G$
2. 学习因子$c_1,c_2$
3. 粒子群规模$N$
4. 可行集：$D\subset {\mathbb{R}}^d$

标准粒子群算法的流程可以描述为：

1 初始化

随机生成（或按照一定分布）粒子群$X^{(t=1)},V^{(t=1)}$，根据初始化条件计算$p_i^{(0)},P_i^{(0)}$，迭代轮数置为$t=1$

2 计算适应度

计算每一个粒子的适应度$f(x_i)$，若$f(x_i)>p_i^{(t-1)}$，则更新$p_i^{(t)}=f(x_i),P_i^{(t)}=x_i$，否则$p_i^{(t)}=p_i^{(t-1)},P_i^{(t)}=P_i^{(t-1)}$

计算出当前种群最优：$s^{(t)},S^{(t)}$

3 更新粒子

按照如下公式更新粒子的速度和位置：

$$\begin{array}{rl}& v_i^{(t+1)}=\omega (t)v_i^{(t)}+c_1{r}_1(P_i^{(t)}-x_i^{(t)})+c_2{r}_2(S^{(t)}-x_i^{(t)})\\ & x_i^{(t+1)}=x_i^{(t)}+v_i^{(t)}\end{array}$$

其中$\omega (t)$是随迭代轮数自适应的惯性因子，而$r_1$和$r_2$是两个$(0,1)$间的随机数.

4 检查结束条件

如果满足结束条件或$t>G$，则算法终止，否则返回到2，并且设置$t=t+1$.

粒子群算法当中的位置和速度更新公式

$$\begin{array}{rl}& v_i^{(t+1)}=\omega (t)v_i^{(t)}+c_1{r}_1(P_i^{(t)}-x_i^{(t)})+c_2{r}_2(S^{(t)}-x_i^{(t)})\\ & x_i^{(t+1)}=x_i^{(t)}+v_i^{(t)}\end{array}$$

模拟的是单个智能体的决策过程。每个智能体都有回到其历史最优值的趋势$P_i^{(t)}-x_i^{(t)}$，前往种群最优值的趋势$S^{(t)}-x_i^{(t)}$，智能体通过一个随机权值组合$(c_1{r}_1,c_2{r}_2)$来将这两个趋势结合起来。同时，智能体由于先前的速度，在物理空间中存在惯性$\omega (t)v_i^{(t)}$. 这三个部分组合在一起，就构成了对速度的改变量。而位置的改变，就是速度，这是符合物理直觉的。

1.2.1 局部粒子群算法

为了增加更多的随机性，避免算法陷入局部最优值，我们还可以对标准粒子群算法进行修改，得到局部粒子群算法。局部粒子群算法相比较于标准粒子群算法的改动在于：当前迭代轮次当中，每个粒子都计算一个局部适应度最优值$s_i^{(t)}$和局部最优位置$S_i^{(t)}$. 局部$s_i^{(t)}$和$S_i^{(t)}$的选择范围，通常是编号与$i$相近的个体。例如我们可以取

$$\begin{gathered} s_i^{(t)}=\text{Best}(f(x_{i-1}^{(t)}),f(x_i^{(t)}),f(x_{i+1}^{(t)})) \\ {S}_i^{(t)}=\text{BestX}(f(x_{i-1}^{(t)}),f(x_i^{(t)}),f(x_{i+1}^{(t)})) \end{gathered}$$

相应的，粒子的更新算法也要修改为

$$\begin{array}{rl}& v_i^{(t+1)}=\omega (t)v_i^{(t)}+c_1{r}_1(P_i^{(t)}-x_i^{(t)})+c_2{r}_2(S_i^{(t)}-x_i^{(t)})\\ & x_i^{(t+1)}=x_i^{(t)}+v_i^{(t)}\end{array}$$

1.3 差分算法

差分算法类似于遗传算法，具有选择、交叉、变异三个流程。

标准的差分算法可以表述为：

1 初始化

随机初始化一个由$N$个$d$维行向量组成的集合：X(0)={xi(0)∈Rd∣i=1,2,…,N}X^{(0)}=\{x_i^{(0)}\in\R^d|i=1,2,\dots,N\}X(0)={xi(0)​∈Rd∣i=1,2,…,N}，称为种群

初始化变异矩阵$H^{(0)}=0\in {\mathbb{R}}^{N\times d}$

初始化子代矩阵$V^{(0)}=0\in {\mathbb{R}}^{N\times d}$

初始化进化代数$t=0$

2 差分变异

对于每一个个体$x_i^{(t)}$，随机选择三个各不相同的个体$x_{p_1}^{(t)},x_{p_2}^{(t)},x_{p_3}^{(t)},(i\ne p_1\ne p_2\ne p_3)$. 计算出变异向量：

$$h_i^{(t+1)}=x_{p_1}^{(t)}+F(x_{p_2}^{(t)}-x_{p_3}^{(t)})$$

如果不存在局部最优问题，也就是环境当中的局部最优适应度等价于整体最优适应度，那么可以把$x_{p_1}^{(t)}$替换为当前种群中适应度最高的个体$x_b^t$，即

$$h_i^{(t+1)}=x_b^{(t)}+F(x_{p_2}^{(t)}-x_{p_3}^{(t)})$$

这样会使得收敛更快

3 交叉操作

对于每一个个体$i,j$，随机采样一个$\lambda \in (0,1)$，计算子代单点位值：

$$v_{ij}^{(t+1)}=\{\begin{array}{ll}{h}_{ij}^{(t+1)},& \lambda <c\\ x_{ij}(t),& \lambda \ge c\end{array}$$

4 选择操作

对于每一个个体编号$i$，判断子代向量$v_i^{(t+1)}$和父代$x_i^{(t)}$的适应度，保留适应度更大的一个：

$$x_i^{(t+1)}=\{\begin{array}{ll}{x}_i^{(t)},& f(x_i^{(t)})\ge f(v_i^{(t+1)})\\ v_i^{(t+1)},& f(x_i^{(t)})<f(v_i^{(t+1)})\end{array}$$

2. 案例：Rosenbrock函数的最优化

Rosenbrock函数的定义式为：

$$\{\begin{array}{ll}f(x_1,x_2)=100(x_1^2-x_2^2)+(1-x_1{)}^2& \\ -2.048\le x_i\le 2.048,& i=1,2\end{array}$$

我们想要求解如下最优化问题

$$\begin{array}{rlr}\underset{x_1,x_2}{\max }& f(x_1,x_2)& \\ \text{s.t.}& -2.048\le x_i\le 2.048,& i=1,2\end{array}$$

2.1 解的编码

Rosenbrock最优化的可行解空间是$[-2.048,2.048{]}^2\subset {\mathbb{R}}^2$，我们可以把这个空间离散化，在$x_1,x_2$轴上，分别以长度$0.004$为单位，划分成离散的区间片段。以区间$[-2.048,-2.044)$编号为$0$，沿着实数轴的正向递增地编号为$0,1,\dots ,1023$.

在两个轴上分别取一个区间，编号组合为$(i,j)$，它们的交叉部分就是一个可行解空间上的矩形区域。把$i,j$都转换成10位的二进制串，再首位拼接为20位的二进制串$x$，就把这个二进制串视为区间$[0.004i-2.048,0.004i-2.044]\times [0.004j-2.048,0.004j-2.044]$内任何一个点的编码，而反过来，用点$(0.004i-2.048,0.004j-2.048)$作为$x$的解码值。

把编码和解码的流程写作算子$\text{encode}(\cdot )$，$\text{decode}(\cdot )$

2.2 应用遗传算法

参数：

1. 种群规模：$N=80$
2. 迭代次数：$G=100$
3. 交叉概率：$p_c=0.6$
4. 变异概率：$p_m=0.1$

流程：

1. 随机生成$N$个20位的二进制串E={x1,…,xN}E=\{x_1,\dots,x_N\}E={x1​,…,xN​}
2. 对每一个$x_i\in E$，解码为浮点数$(u_i,v_i)=\text{decode}(x_i)$
3. 计算每一个解码值的Rosenbrock值$f_i=f(u_i,v_i)$，并作为对应个体的适应度
4. 按适应度对$E$进行排序，记录下$x_b=x_N$，删除排在$x_N$以后的个体（如果有的话）
5. 选择交叉种群$E_c=\text{Select}(E)$
6. 进行交叉计算$E=\text{Crossover}(E_c)$
7. 把$E$当中的第$N$个个体置为$x_m$（这是利用了最优保存策略）
8. 进行变异计算$E=\text{Mutation}(E_m)$
9. 把$E$当中的第$N$个个体置为$x_m$
10. 判断是否达到最大迭代次数$G$，若达到则终止迭代，否则返回第2步
11. 取终止迭代时，群体中的最优个体为算法的结果

2.2.1 选择算子

$E_c=\text{Select}(E)$是这样计算的：

1. 创建一个空集合$E_c$
2. 计算每个个体的适应度，与平均适应度的比值（下取整）
   $$\begin{array}{rl}& \overline{f}=\sum _{i=1}^N{f}_i\\ & q_i=\text{floor}(f_i/\overline{f})\end{array}$$
3. 对于每个个体$x_i\in E$，复制$q_i$倍，加入到$E_c$当中
4. 返回$E_c$

2.2.2 交叉算子

$E_m=\text{Crossover}(E_c)$的计算规则是这样的：

1. 创建一个空集合$E_m$
2. 随机选取一个交叉点位$n\le 20$
3. 随机打乱$E_c$当中个体的排序，记打乱后的个体是$x_i^{\prime }$
4. 对于每一个$i=1,3,\dots ,N-1orN-2$，进行概率为$p_c$的二项决策，决定是否进行交叉
5. 如果决定交叉，那么就将$x_i^{\prime }$和$x_{i+1}^{\prime }$在$n$到20位之间的二进制片段交换，得到${\overline{x}}_i,{\overline{x}}_{i+1}$
6. 若不交换，则${\overline{x}}_i,{\overline{x}}_{i+1}$就是$x_i^{\prime }$和$x_{i+1}^{\prime }$
7. 把${\overline{x}}_i,{\overline{x}}_{i+1}$加入到$E_m$
8. 返回$E_m$

2.2.3 变异算子

$E=\text{Mutation}(E_m)$是这样计算的：

1. 创建空集合$E^{\prime }$
2. 对于每一个$i=1,2,\dots ,N$，进行概率为$p_m$的二项决策，决定是否进行变异
3. 如果决定进行变异，则将个体${\overline{x}}_i$的每一位取反，得到${\hat{x}}_i$，否则${\hat{x}}_i={\overline{x}}_i$
4. 把${\hat{x}}_i$加入到$E^{\prime }$
5. 返回$E=E^{\prime }$

2.3 应用局部粒子群算法

参数：

1. 粒子群规模：$N=50$
2. 学习因子：$c_1=1.4,c_2=1.7$
3. 自适应惯性因子：$\omega (k)$.
   $$\omega (k)=\frac{k}{G}({\omega }_{\text{max}}-{\omega }_{\text{min}})$$
4. 边界：$V_{\text{min}}=-1,V_{\text{max}}=1,{\omega }_{\text{max}}=0.1,{\omega }_{\text{min}}=0.9$

算法流程：

1. 初始化粒子群$\{(x_1,v_1),\dots ,(x_N,v_N)\}$，其中$x_i$和$v_i$都是20位的二进制串. 按照当前的粒子群分布，初始化历史最优值$p_i,P_i$，初始化粒子局部最优值最优值$s_i,S_i$，设置迭代次数$k=1$
2. 对每个粒子更新速度$v_i=\omega (k)v_i+c_1{r}_1(P_i-x_i)+c_2{r}_2(S_i-x_i)$
   ，$r_1$和$r_2$是随机采样的$(0,1)$之间的值。对于超过区间$[V_{\text{min}},V_{\text{max}}]$的值，取边界值。
3. 对每个粒子更新位置$x_i=x_i+v_i$
4. 对每一个粒子$x_i$，计算适应度$f_i=f(\text{decode}(x_i))$
5. 如果$f_i>p_i$，则更新$p_i=f_i,P_i=X_i$
6. 计算每个粒子的局部最优值$s_i,S_i$
7. 如果$k>G$，终止。否则更新$k=k+1$并且回到2
8. 取终止迭代时，粒子群中的最优粒子为算法的结果

2.4 应用差分进化算法

参数：

1. 种群规模：$N=30$
2. 迭代次数：$G=50$
3. 变异因子：$F=1.2$
4. 交叉因子：$c=0.9$

流程：

1. 随机初始化种群E={x1,x2,…,xN}E=\{x_1,x_2,\dots,x_N\}E={x1​,x2​,…,xN​}，计算最优值$x_b$和最优适应度$f_b=f(x_b)$，置迭代次数为$t=0$
2. 对于每一个$x_i$，随机选择$i\ne r_1\ne r_2\ne r_3$，计算差分变异：

$$h_i=x_{r_1}+F(x_{r_2}-x_{r_3})$$

3. 对于每一个$i,j$，随机生成一个$\lambda \in (0,1)$，计算交叉子代矩阵

$$v_{ij}=\{\begin{array}{ll}{h}_{ij},& \lambda <c\\ x_{ij},& \lambda \ge c\end{array}$$

4. 对于每一个$i$，进行选择

$$x_i=\{\begin{array}{ll}{x}_i,& f(x_i)\ge f(v_i)\\ v_i,& f(x_i)<f(v_i)\end{array}$$

5. 判断是否得到最大迭代次数，若是则终止迭代，否则置$t=t+1$
6. 取终止迭代时群体中的最优个体作为算法的结果