智能优化算法之模拟退火（Simulated Annealing，SA）-附源码

获取更多资讯，赶快关注上面的公众号吧！

模拟退火（Simulated Annealing，SA）

本章将介绍模拟退火算法，它是受冶金中退火过程的启发而提出的一种元启发式优化算法。首先回顾退火算法的发展和应用，然后描述物理退火的过程及其与退火算法的映射关系，并详细描述算法的步骤，最后给出SA算法的伪代码和java代码实现。

灵感

SA由Kirkpatrick等人于1983年提出，其揭示了如何利用模拟固体退火过程的模型求解优化问题，在优化问题中，最小化目标函数对应于固体中的能量状态。

在物理退火过程中，固体首先被加热，然后慢慢冷却，直到到达最规则的晶格排列而没有晶格缺陷。加热会改变物质的物理甚至是化学特性，以增加其延展性并降低其硬度。固体中的粒子具有与最稳定状态下最小能量排列相对应的几何构型，物理退火是通过熔化一种物质，然后缓慢降低其温度来实现固体的低能排列的过程，值得注意的是，在冷却阶段，物质的温度必须缓慢下降，否则，生成的固体就会凝固成存在结构缺陷的晶体。

当应用SA时，物质的不同状态表示优化问题的不同解，该物质的能量相当于待优化的适应度函数，原子的移动引入新的解。在SA中，引入更优解的原理移动将被接受，引入更劣解的非改进移动以概率接受，此概率取决于接受函数。下表列出了模拟退火与优化算法之间的对应关系。

优化算法	模拟退火
决策变量	物质原子位置
解	物质状态
旧解	物质当前状态
新解	物质新的状态
最优解	最优状态
适应度函数	物质能量
初始解	随机位置
选择	接受函数
生成新解的过程	原子移动

SA首先生成一个初始解，也就是物质的当前状态，通过适当的机制在当前状态的邻域中生成新的状态，并评估其适应度值，如果新的状态优于当前状态，则接受新的状态，否则以一定概率接受，此概率与系统温度有关，该算法在温度逐渐降低的同时，在每个温度下生成一定数量的新状态，不断尝试一部分解，直达每个温度下都达到热平衡准则，此时再次降低温度。随着温度的降低，选择非改进原子移动的概率也随之降低。整个生成新解的过程随着系统冷却逐渐重复，直到满足终止条件。SA的算法流程图如下所示。

图1 SA流程图

生成初始状态

SA生成的每个优化问题的解对应于物质原子的一种排列。系统状态包括在N维空间中的N个决策变量，可以表示为一个大小为$1\times N$的数组，如下：

$$\text{State}=X=\left(x_1,x_2,\dots ,x_i,\dots ,x_N\right)\text{\hspace{1em}(1)}$$

其中$X$为优化问题的一个解，$x_i$为解中$X$的第$i$个决策变量，$N$为决策变量的个数。对于连续问题和离散问题，决策变量值$\left(x_1,x_2,\dots ,x_i,\dots ,x_N\right)$可以表示为实数或一组预定义的值。初始时随机生成一个初始化状态，该状态就是系统的当前状态。

生成新的状态

在退火过程中，原子移动到新的位置，以降低系统能量，实现稳定状态。在当前状态基础上生成新的系统可能状态，有多种确定性的和随机性的机制可以根据给定的解生成邻域解。最常用的一种机制就是随机行走，即：

$$x_i^{\prime }=x_i+Rnd(-\epsilon ,\epsilon ),i=1,2,\dots ,N\text{\hspace{1em}(2)}$$

其中$x_i$第$i$个决策变量的新值，$Rnd(a,b)$为$[a,b]$范围内的随机值，$\epsilon$为一个较小的数。

所有决策变量的新值在生成新解 之前先进行评估，然后新解按下式生成：
$$X^{(new)}=\left(x_1^{\prime },x_2^{\prime },\dots ,x_i^{\prime },\dots ,x_N^{\prime }\right)\text{\hspace{1em}(3)}$$

接受函数

接受函数用于确定是否使用新生成的解替换当前解，新解代表物质的一种可能的排列，其接受与否取决于旧解和新解的适应度值。当新解适应度值优于旧解时，新解替换旧解，否则新解一定概率替换旧解，该概率根据新旧解适应度值之间的差异计算得到。针对最小化问题，根据如下的接受概率函数来接受新解：
$$P\left(X,X^{(new)}\right)=\left\{\begin{array}{ll}1& \text{ if }F\left(X^{(new)}\right)<F(X)\\ e^{\frac{-\Delta F}{\lambda }}& \text{ Otherwise }\end{array}\right\}\text{\hspace{1em}(4)}$$

$$\Delta F=\left|F\left(X^{(new)}\right)-F(X)\right|\text{\hspace{1em}(5)}$$

其中$X$为旧解，$X^{(new)}$为新生成的解，$P\left(X,X^{(new)}\right)$大于等于$Rand\in [0,1]$为使用$X^{(new)}$替换$X$的概率，$F(X)$为解$X$的适应度值，$\lambda$为控制参数，对应于物理退火中的温度。如果$P\left(X,X^{(new)}\right)$大于等于$Rand\in [0,1]$，则使用$X^{(new)}$替换$X$，否则就不接受新解。式（4）和式（5）定义的接受函数为玻尔兹曼分布，当适应度值差异较大时，概率将变小，因此适应度值差异较小时更容易被接受。同理，$\lambda$较大时非改进改变更容易被接受，即当$\lambda$较大时选择压力较小。参数$\lambda$对于算法正确收敛起着重要的作用，算法初期$\lambda$应该较大，以避免陷入局部最优，随着算法逐渐进行其也应该逐渐下降。

热平衡

当在每一温度下发生大量的邻域运动时，系统在每一温度下达到热平衡，可以证明，在热平衡时，系统状态的概率分布服从玻尔兹曼分布。SA通过尝试在每个温度下的多个邻域移动来进行，换言之，对于每个温度$\lambda$，在温度降低之前，需要生成多个新的状态（不管是接受还是拒绝），这些新的状态需要通过接受函数进行测试，至于需要生成多少个新的状态，是用户定义的算法参数，通常表示为$\beta$。当生成了$\beta$个新解并由接受函数测试后，就满足了热平衡。

温度下降

在测试多个新的状态后，系统温度逐渐下降：

$$\underset{t\to +\infty }{\lim }{\lambda }_t=0,t>0\text{\hspace{1em}(6)}$$
其中$t$为算法迭代次数。

两种常见的降低$\lambda$的方法分别是线性法和几何法。

线性法使用下式在每次迭代中更改$\lambda$：
$${\lambda }_t={\lambda }_0-\alpha \times t\text{\hspace{1em}(7)}$$

$$\alpha =\frac{{\lambda }_0-{\lambda }_T}{T}\text{\hspace{1em}(8)}$$

其中${\lambda }_0$为初始温度，${\lambda }_t$为第$t$次迭代时的温度，$T$为总迭代次数，$\alpha$为冷却因子。

系统冷却的几何法如下：
$${\lambda }_t={\lambda }_0\times {\alpha }^t,0<\alpha <1\text{\hspace{1em}(9)}$$

几何法的优点在于其不需要指定算法的最大迭代次数。算法的终止条件可以是迭代$T$的最大次数，也可以是其他终止条件，如运行时间。在这种情况下，T的值是未知的，SA算法将继续执行，直到满足停止条件(例如，运行时间)。

终止条件

终止条件决定何时终止SA算法。选择合适的终止准则对算法的正确收敛有重要作用。迭代次数、连续迭代之间目标函数的增量改进和运行时间是SA算法实现中常用的终止条件。

伪代码和Java实现

开始

  输入算法参数和初始数据；

  生成初始解$X$并进行评估；

  While(终止条件未满足)

    For $j$=1 to $\beta$

      生成新解$X^{(new)}$，并进行评估；

      If 新解$(X^{(new)})$优于旧解$(X)$

        令$X=X^{(new)}$

      Otherwise

        计算$P\left(X,X^{(new)}\right)$，生成[0,1]内的随机数$Rand$;

        If $P\left(X,X^{(new)}\right)$>$Rand$

          令$X=X^{(new)}$

        End if

      End if

    Next $j$

    降低温度

  End while

  输出解$X$

End

下面结合代码具体说一下算法执行流程，java源码关注公众号后回复“模拟退火”即可获取。

1. 初始参数设置。主要设置最大迭代次数 maxGen，求解问题维度 dim，初始温度T0，温度下降率alpha，热平衡beta和目标函数 objectiveFun 等。

ObjectiveFun objectiveFun = new ObjectiveFun();
SANormal sa = SANormal.builder().maxGen(400).dim(2).T0(10).beta(200).objectiveFun(objectiveFun).build();

2. 生成初始解。

private Substance genInitSolution() {
	// TODO Auto-generated method stub
	Substance substance = new Substance(new Position(dim, objectiveFun.getRange()));
	substance.setEnergy(objectiveFun.getObjValue(substance.getPosition().getPositionCode()));
	return substance;
}

3. 生成新的解。采用随机游走的方式对旧解进行扰动。

private Substance genNewSolution(Substance oldSubstance) {
	// TODO Auto-generated method stub
	Substance newSubstance = new Substance(new Position(dim, objectiveFun.getRange()));
	double[] newPositionCode = IntStream.range(0, this.dim)
			.mapToDouble(ind -> oldSubstance.getPosition().getPositionCode()[ind]
					+ 0.1 * objectiveFun.getRange().getScale()[ind] * (random.nextDouble() - 0.5))
			.toArray();
	newSubstance.getPosition().setPositionCode(newPositionCode);
	newSubstance.setEnergy(this.objectiveFun.getObjValue(newSubstance.getPosition().getPositionCode()));
	return newSubstance;
}

4. 接受新解与否。

// 针对最大化问题
if (newSubstance.getEnergy() > bestSubstance.getEnergy()) {
	substance = newSubstance;
} else {
	double deldaF = Math.abs(newSubstance.getEnergy() - bestSubstance.getEnergy());
	double p = Math.pow(Math.E, -deldaF / T);
	double rand = random.nextDouble();
	if (rand <= p) {
		substance = newSubstance;
	}
}

5. 更新迭代器。迭代次数加1，温度降低。

public void incrementIter() {
  iterator++;
  T = T * alpha;
}

6. 循环。

while (iterator < maxGen) {
	for (int j = 0; j < beta; j++) {
		Substance newSubstance = genNewSolution(substance);
		// 针对最大化问题
		if (newSubstance.getEnergy() > bestSubstance.getEnergy()) {
			substance = newSubstance;
		} else {
			double deldaF = Math.abs(newSubstance.getEnergy() - bestSubstance.getEnergy());
			double p = Math.pow(Math.E, -deldaF / T);
			double rand = random.nextDouble();
			if (rand <= p) {
				substance = newSubstance;
			}
		}
	}
	if (substance.getEnergy() > bestSubstance.getEnergy()) {
		bestSubstance.getPosition().setPositionCode(substance.getPosition().getPositionCode());
		bestSubstance.setEnergy(substance.getEnergy());
	}
	incrementIter();
	System.out.println("**********第" + iterator + "代最优解：" + bestSubstance + "**********");
	T = T * alpha;
}

测试中使用的目标函数（最大化）方程和图像分别如下：

exp(-(x-4)^2-(y-4)^2)+exp(-(x+4)^2-(y-4)^2)+2*exp(-x^2-(y+4)^2)+2*exp(-x^2-y^2)

实验结果如下: