模拟退火算法原理及C++实例

1.1算法简介
模拟退火算法得益于材料的统计力学研究成果。统计力学表明材料中粒子的不同结构对应于粒子的不同能量水平。在高温条件下,粒子的能量较高,可以自由运动和重新排序。在低温条件下,粒子的能力较低。如果聪高温开始,非常缓慢的降温(这个过程被称为退火),粒子就可以在每个温度下达到热平衡。当系统完全被冷却时,最终形成处于低能状态的结晶。
如果粒子的能量定义材料的状态,Metropolis算法用一个简单的数学模型描述了退火过程。假设材料在状态 i i i之下的能量为 E ( i ) E(i) E(i),那么材料在温度 T T T时从状态 i i i进入状态 j j j就遵循如下规律:
(1)如果 E ( j ) ≤ E ( i ) E(j)\leq E(i) E(j)E(i),接受该状态被转换。
(2)如果 E ( j ) > E ( i ) E(j)>E(i) E(j)>E(i),则状态以如下概率接受:
e E ( i ) − E ( j ) K T e^{\frac{E(i)-E(j)}{ KT}} eKTE(i)E(j)
其中 K K K是物理学中的玻尔兹曼常数, T T T是材料温度。
在某一特定温度下,进行了充分的转换之后,材料将达到热平衡。这是材料处于状态 i i i的概率满足玻尔兹曼分布:
P T ( x = i ) = e − E ( i ) K T ∑ j ∈ S e − E ( i ) K T P_{T}(x=i)=\frac{e^{-\frac{E(i)}{KT}}}{\sum\limits_{j\in S}e^{-\frac{E(i)}{KT}}} PT(x=i)=jSeKTE(i)eKTE(i)
其中 x x x表示材料当前状态的随机变量, S S S表示状态空间集合。
显然
lim ⁡ T → ∞ e − E ( i ) K T ∑ j ∈ S e − E ( i ) K T = 1 ∣ S ∣ \lim\limits_{T\to\infty}\frac{e^{-\frac{E(i)}{KT}}}{\sum\limits_{j\in S}e^{-\frac{E(i)}{KT}}}=\frac{1}{|S|} TlimjSeKTE(i)eKTE(i)=S1
其中 ∣ S ∣ |S| S表示集合S中状态的数量。当当前温度很高趋于无穷时,那么当前的能量可近似达到一个最大值,那么它只能沿着小能量方向变化,这时遵循规律(1),状态直接被接受。这表明所有状态在高温下具有相同的概率。而当温度下降时,
lim ⁡ T → 0 e − E ( i ) − E m i n K T ∑ j ∈ S e − E ( i ) − E m i n K T = lim ⁡ T → 0 e − E ( i ) − E m i n K T ∑ j ∈ S m i n e − E ( i ) − E m i n K T + ∑ j ∉ S m i n e − E ( i ) − E m i n K T = lim ⁡ T → 0 e − E ( i ) − E m i n K T ∑ j ∉ S e − E ( i ) − E m i n K T = { 1 ∣ S m i n ∣ , 若 i ∈ S m i n 0 , 其 他 \begin{aligned}\lim\limits_{T\to0}\frac{e^{-\frac{E(i)-E_{min}}{KT}}}{\sum\limits_{j\in S}e^{-\frac{E(i)-E_{min}}{KT}}}&=\lim\limits_{T\to0}\frac{e^{-\frac{E(i)-E_{min}}{KT}}}{\sum\limits_{j\in S_{min}}e^{-\frac{E(i)-E_{min}}{KT}}+\sum\limits_{j\notin S_{min}}e^{-\frac{E(i)-E_{min}}{KT}}}\\ &=\lim\limits_{T\to0}\frac{e^{-\frac{E(i)-E_{min}}{KT}}}{\sum\limits_{j\notin S}e^{-\frac{E(i)-E_{min}}{KT}}} \\&=\begin{cases} \frac{1}{|S_{min}|},\qquad 若i\in S_{min}\\ 0, \qquad\qquad 其他 \end{cases} \end{aligned} T0limjSeKTE(i)EmineKTE(i)Emin=T0limjSmineKTE(i)Emin+j/SmineKTE(i)EmineKTE(i)Emin=T0limj/SeKTE(i)EmineKTE(i)Emin={Smin1iSmin0
其中 E m i n = min ⁡ j ∈ S E ( j ) E_{min}= \min \limits_{j\in S}E(j) Emin=jSminE(j) S m i n = { i ∣ E ( i ) = E m i n } S_{min}=\{i|E(i)=E_{min}\} Smin={iE(i)=Emin}
因为 E m i n E_{min} Emin为最低温度下的能量最小值,那么其他任何状态下的值都会比其更高,所以遵循规律(2)。上式表明当温度降低至很低时,材料会以很大概率进入最小能量状态。

理论运用:
假定我们要解决的问题是寻找最小值的优化问题。将物理学中模拟退火的思想应用于优化问题就可以得到模拟退火寻优方法。
考虑这样一个组合优化问题:优化函数为: f : x → R + f:x \to R^+ f:xR+,其中 x ∈ S x \in S xS,他表示优化问题的一个可行解, R + = { y ∣ y ∈ R , y > 0 } R^+=\{y|y\in R,y>0\} R+={yyR,y>0} S S S表示函数的定义域。 N ( x ) ⫅ S N(x)\subseteqq S N(x)S表示 x x x的一个领域集合。
首先给定一个初始温度 T 0 T_{0} T0和该优化问题的一个初始解 x 0 x{0} x0,并由 x ( 0 ) x(0) x(0)生成下一个解 x ′ ∈ N ( x ( 0 ) ) x^{'} \in N(x(0)) xN(x(0)),是否作为一个新解 x ( 1 ) x(1) x(1)依赖于下面概率:
P ( x ( 0 ) → x ′ ) = { 1 若 f ( x ′ ) &lt; f ( x ( 0 ) ) e − f ( x ′ ) − f ( x ( 0 ) ) T 0 其 他 P(x(0)\to x^{&#x27;})=\begin{cases} 1 \qquad\qquad\qquad\qquad 若f(x^{&#x27;})&lt;f(x(0))\\ e^{-\frac{f(x^{&#x27;})-f(x(0))}{T_0}} \qquad\quad其他 \end{cases} P(x(0)x)=1f(x)<f(x(0))eT0f(x)f(x(0))
换句话说,如果生成的解 x ′ x^{&#x27;} x的函数值比前一个解的函数值更小,则接受 x ( 1 ) = x ′ x(1)=x^{&#x27;} x(1)=x作为一个新解。否则以概率 e − f ( x ′ ) − f ( x ( 0 ) ) T 0 e^{-\frac{f(x^{&#x27;})-f(x(0))}{T_0}} eT0f(x)f(x(0))接受 x ′ x^{&#x27;} x作为一个新解。
泛泛地说,对于某一个温度 T i T_i Ti和该优化问题的一个解 x ( k ) x(k) x(k),可以生成 x ′ x^{&#x27;} x。接受 x ′ x^{&#x27;} x作为下一个新解 x ( k + 1 ) x(k+1) x(k+1)的概率为:
(1) P ( x ( 0 ) → x ′ ) = { 1 若 f ( x ′ ) &lt; f ( x ( 0 ) ) e − f ( x ′ ) − f ( x ( 0 ) ) T 0 其 他 P(x(0)\to x^{&#x27;})=\begin{cases} 1 \qquad\qquad\qquad\qquad 若f(x^{&#x27;})&lt;f(x(0))\\ e^{-\frac{f(x^{&#x27;})-f(x(0))}{T_0}} \qquad\quad其他 \end{cases}\tag 1 P(x(0)x)=1f(x)<f(x(0))eT0f(x)f(x(0))(1)
在温度 T i T_{i} Ti下,经过很多次的转移之后,降低温度 T i T_{i} Ti,得到 T I + 1 &lt; T i T_{I+1}&lt;T_{i} TI+1<Ti。在 T i + 1 T_{i+1} Ti+1下重复上述过程。因此整个优化过程就是不断寻找新解和缓慢降温的交替过程。最终的解是对该问题寻优的结果。
我们注意到,在每个 T i T_{i} Ti下,所得到的一个新状态 x ( k + 1 ) x(k+1) x(k+1)完全依赖于前一个状态 x ( k ) x(k) x(k),可以和前面的状态 x ( 0 ) , . . . , x ( k − 1 ) x(0),…,x(k-1) x(0),...,x(k1)无关,因此是一个马尔可夫过程。使用马尔可夫过程对上述模拟退火的步骤进行分析,结果表明:从任何一个状态 x ( k ) x(k) x(k)生成 x ′ x^{&#x27;} x的概率,在 N ( x ( k ) ) N(x(k)) N(x(k))中是均匀分布的,且新状态 x ′ x^{&#x27;} x被接受的概率满足式(1),那么经过有限次转换,在温度 T i T_{i} Ti下平衡态 x i x_{i} xi的分布由下式给出:
P i ( T i ) = e − f ( x i ) T ∑ j ∈ S e − f ( x i ) T i P_{i}(T_{i})=\frac{e^{-\frac{f(x_{i})}{T}}}{\sum\limits_{j\in S}e^{-\frac{f(x_i)}{T_{i}}}} Pi(Ti)=jSeTif(xi)eTf(xi)
当温度 T T T降为0时, x i x_{i} xi的分布为:
P i ∗ = { 1 ∣ S m i n ∣ 若 x i ∈ S m i n 0 其 他 P^*_{i}=\begin{cases} \frac{1}{|S_{min}|}\qquad\qquad 若x_{i}\in S_{min}\\ 0 \qquad\qquad\qquad 其他 \end{cases} Pi={Smin1xiSmin0
并且
∑ x i ∈ S m i n P i ∗ = 1 \sum\limits_{x_{i}\in S_{min}}P^*_i=1 xiSminPi=1
这说明如果温度如果温度下降十分缓慢,而在每个温度都有足够多次的状态状态转移,使之在每一个温度下达到热平衡,则全局最优解将以概率1被找到。因此可以说模拟退火算法可以找到全局最优解。
在模拟退火算法中,应注意以下问题:
(1)理论上,降温过程要足够缓慢,要使得在每一温度下达到热平衡。但在计算机实现中,如果降温速度过缓,所得到的解性能会较为令人满意,但是算法会太慢,相对于简单的搜索算法不具有明显的优势。如果降温速度过快,很可能最终得不到全局最优解。因此使用时要综合考虑解的性能和算法速度,在两者之间采取一种折中。
(2)要确定在每一温度下状态转换的结束准则。实际操作可以考虑当连续 m m m次的转换过程没有使状态发生变化时结束该温度下的状态转换。最终温度的确定可以提前定位一个较小值 T e T_{e} Te,或连续几个温度下转换过程没有使状态发生变化算法就结束。
(3)选择初始温度和确定某个可行解的领域的方法也要恰当。

1.2 应用举例
例 已知地方100个目标的经、纬度如表所示。
《模拟退火算法原理及C++实例》
《模拟退火算法原理及C++实例》
我方有一个基地,经度和维度为(70,40)。假设我方飞机的速度为1000公里/小时。我方派一架飞机从基地出发,侦察完敌方所有目标,再返回原来基地。在敌方每一目标点的侦察时间不计,求该架飞机所花费的时间(假设我方飞机巡航时间充分长)。
c++求解代码见链接戳这里.

点赞