模拟退火算法之scikit Opt库使用

Table of Contents

一、模拟退火算法概念

1. 简介

模拟退火算法（Simulate Annealing Arithmetic，SAA）**是一种通用概率演算法，用来在一个大的搜寻空间内找寻命题的最优解。**它是基于Monte-Carlo迭代求解策略的一种随机寻优算法。

模拟退火算法思想源于固体退火原理：

• 固体退火（物理退火）：将材料加热后再经特定速率冷却，目的是增大晶粒的体积，并且减少晶格中的缺陷。材料中的原子原来会停留在使内能有局部最小值的位置，加热使能量变大，原子会离开原来位置，而随机在其他位置中移动。退火冷却时速度较慢，使得原子有较多可能可以找到内能比原先更低的位置。

• 模拟退火：其原理也和固体退火的原理近似。模拟退火算法从某一较高初温出发，伴随温度参数的不断下降,结合概率突跳特性在解空间中随机寻找目标函数的全局最优解，即在局部最优解能概率性地跳出并最终趋于全局最优。

2. 基本思想

1. 爬山算法

   介绍模拟退火前，先介绍爬山算法。爬山算法是一种简单的贪心搜索算法，该算法每次从当前解的临近解空间中选择一个最优解作为当前解，直到达到一个局部最优解。

   ​ 爬山算法实现很简单，其主要缺点是会陷入局部最优解，而不一定能搜索到全局最优解。如图1所示：假设C点为当前解，爬山算法搜索到A点这个局部最优解就会停止搜索，因为在A点无论向那个方向小幅度移动都不能得到更优的解。

   

2. 模拟退火算法

   爬山法是完完全全的贪心法，这种贪心是很鼠目寸光的，只把眼光放在局部最优解上，因此只能搜索到局部的最优值。模拟退火其实也是一种贪心算法，只不过与爬山法不同的是，模拟退火算法在搜索过程引入了随机因素。模拟退火算法以一定的概率来接受一个比当前解要差的解，因此有可能会跳出这个局部的最优解，达到全局的最优解。

二、模拟退火算法流程

模拟退火算法新解的产生和接受可分为如下四个步骤：

1. 第一步是由一个产生函数从当前解产生一个位于解空间的新解；为便于后续的计算和接受，减少算法耗时，通常选择由当前新解经过简单地变换即可产生新解的方法，如对构成新解的全部或部分元素进行置换、互换等，注意到产生新解的变换方法决定了当前新解的邻域结构，因而对冷却进度表的选取有一定的影响。

2. 第二步是计算与新解所对应的目标函数差。因为目标函数差仅由变换部分产生，所以目标函数差的计算最好按增量计算。事实表明，对大多数应用而言，这是计算目标函数差的最快方法。

3. 第三步是判断新解是否被接受，判断的依据是一个接受准则，最常用的接受准则是Metropolis准则: 若Δt′<0则接受S′作为新的当前解S，否则以概率exp(-Δt′/T)接受S′作为新的当前解S。

4. 第四步是当新解被确定接受时，用新解代替当前解，这只需将当前解中对应于产生新解时的变换部分予以实现，同时修正目标函数值即可。此时，当前解实现了一次迭代。可在此基础上开始下一轮试验。而当新解被判定为舍弃时，则在原当前解的基础上继续下一轮试验。

三、模拟退火算法工具箱

1. 快速开始

1. 定义问题

   1	demo_func = lambda x: x[0] ** 2 + (x[1] - 0.05) ** 2 + x[2] ** 2

2. 运行模拟退火算法

   1 2 3 4 5

   from sko.SA import SA sa = SA(func=demo_func, x0=[1, 1, 1], T_max=1, T_min=1e-9, L=300, max_stay_counter=150) best_x, best_y = sa.run() print('best_x:', best_x, 'best_y', best_y)

3. 画图

   1 2 3 4 5

   import matplotlib.pyplot as plt import pandas as pd plt.plot(pd.DataFrame(sa.best_y_history).cummin(axis=0)) plt.show()

   

2. 算法参数

模拟退火算法参数控制

• 温度T的初始值设置

  温度T的初始值设置是影响模拟退火算法全局搜索性能的重要因素之一、初始温度高，则搜索到全局最优解的可能性大，但因此要花费大量的计算时间；反之，则可节约计算时间，但全局搜索性能可能受到影响。实际应用过程中，初始温度一般需要依据实验结果进行若干次调整。

• 退火速度问题

  模拟退火算法的全局搜索性能也与退火速度密切相关。一般来说，同一温度下的“充分”搜索(退火)是相当必要的，但这需要计算时间。实际应用中，要针对具体问题的性质和特征设置合理的退火平衡条件。

  这个参数一般设置为0.99或0.98，越小降温越快，如果设置过小，就如同淬火一下，算法来不及充分寻优。

  更多需要参考一些论文

3. 算法进阶

三种形式的模拟退火算法具体参考https://scikit-opt.github.io/scikit-opt/#/zh/more_sa

模拟退火有三种具体形式

fast:

1
2
3
4
5
6
7
8

u ~ Uniform(0, 1, size = d)
y = sgn(u - 0.5) * T * ((1 + 1/T)**abs(2*u - 1) - 1.0)

xc = y * (upper - lower)
x_new = x_old + xc

c = n * exp(-n * quench)
T_new = T0 * exp(-c * k**quench)

cauchy:

1
2
3
4
5

u ~ Uniform(-pi/2, pi/2, size=d)
xc = learn_rate * T * tan(u)
x_new = x_old + xc

T_new = T0 / (1 + k)

boltzmann:

1
2
3
4
5

std = minimum(sqrt(T) * ones(d), (upper - lower) / (3*learn_rate))
y ~ Normal(0, std, size = d)
x_new = x_old + learn_rate * y

T_new = T0 / log(1 + k)

问题分析

通过运行发现他的文档提供的代码不对，如下

1
2
3
4
5
6
7
8

from sko.SA import SAFast

demo_func = lambda x: x[0] ** 2 + (x[1] - 8) ** 2 + x[2] ** 2

sa_fast = SAFast(func=demo_func, x0=[1, 1, 1], T_max=1, T_min=1e-9, q=0.99, L=300, max_stay_counter=150,
                 lb=[-1, 1, -1], ub=[2, 3, 4])
sa_fast.run()
print('Fast Simulated Annealing with bounds: best_x is ', sa_fast.best_x, 'best_y is ', sa_fast.best_y)

结果：

可以看到其中解超出原本该在的范围，他指定了lb、ub以及q，发现其实是没有起到作用的，因为 SAFast 类中没有对应的变量。只有lower、upper、queue，没有对应的参数。

通过测试相关值也可以知道还是默认值

可以通过改变参数设置来解决他的问题，但是在源码中关于范围的设置是有点问题的，lower、upper 用于计算时使用的数值类型，不是 list，会发生报错，可以继承类来实现，也可以直接修改代码解决。

**其实在源码中的lower、upper并没有对新值的产生起到约束作用，只是用于更新新值。**如果希望对变量范围进行约束，需要继承该类实现，可参考模拟退火算法处理有约束问题。

四、总结

模拟退火算法优缺点

• 优点：计算过程简单，通用，鲁棒性强，适用于并行处理，可用于求解复杂的非线性优化问题

• 缺点：收敛速度慢，执行时间长，算法性能与初始值有关及参数敏感等缺点

模拟退火算法则更适合处理无约束优化问题，采用一些简单的处理方法也可以适用于简单的约束优化问题。具体参考[Python数模笔记-模拟退火算法（2）约束条件的处理]