1 模拟退火物理背景介绍

金属退火是将金属加热到一定温度，保持足够时间，然后以适宜速度冷却(通常是缓慢冷却，有时是控制冷却)的一种金属热处理工艺。而模拟退火算法来源于固体退火原理，将固体加温至充分高，再让其徐徐冷却，加温时，固体内部粒子随温升变为无序状，内能增大，而徐徐冷却时粒子渐趋有序，在每个温度都达到平衡态，最后在常温时达到基态，内能减为最小。

如上图，处在低温状态时，固体中分子具有的内能很低，在原本的位置上做小范围的振动。若是将固体加热到一定温度，分子内能将会增加，热运动加剧，分子排列的无序度增加。此时再将温度缓缓降低，在每个温度都达到平衡态（即准静态过程），分子具有的能量逐渐降低，最终回归到有序排列的状态，分子内能也跟着降到最低。

2 模拟退火基本原理

2.1 局部寻优与全局寻优

模拟退火算法(Simulated Annealing, SA)最早的思想是由N. Metropolis等人于1953年提出。模拟退火算法是基于Monte-Carlo 迭代求解策略的一种随机寻优算法，其出发点是基于物理中固体物质的退火过程与一般组合优化问题之间的相似性。介绍模拟退火前，有必要先介绍一下爬山算法。

爬山算法是一种简单的贪心搜索算法，该算法每次从当前解的临近解空间中选择一个最优解作为当前解，直到达到一个局部最优解。爬山算法实现很简单，其主要缺点是会陷入局部最优解，而不一定能搜索到全局最优解。如上图所示：我们要求取定义域内该曲线的最大值。假设C点为当前解，爬山算法搜索到A点这个局部最优解就会停止搜索，因为在A点无论向哪个方向小幅度移动都不能得到更优的解。而模拟退火算法则是从某一较高初温出发，伴随温度参数的不断下降，结合一定的概率突跳特性在解空间中随机寻找目标函数的全局最优解，即在局部最优解能概率性地跳出并最终趋于全局最优。下图关于求解定义域范围内函数最大值的动图可以帮助我们更好地理解本模型：

这里的“一定的概率”的计算参考了金属冶炼的退火过程，这也是模拟退火算法名称的由来。将温度T当作控制参数，目标函数值f视为内能E，而固体在某温度T时的一个状态对应一个解，然后算法试图随着控制参数T的降低，使目标函数f(内能E)也逐渐降低，直至趋于全局最小值（退火中低温时的最低能量状态），就像金属退火过程一样。

2.2 Metropolis准则

上文提到，模拟退火算法可以结合一定的概率突跳特性在解空间中随机寻找目标函数的全局最优解，即在局部最优解能概率性地跳出并最终趋于全局最优。具体是怎么实现”突跳”的呢？设$E_t$为当前能量值，$E_{t+1}$为下一时刻的能量值，需要求解的是全范围内的能量最小值，接受概率为$P$,那么：

若下一时刻的能量值相较当前能量值更小，则完全接受该点的坐标与能量值，$P = 1$
若下一时刻的能量值更大，则以一定的概率接受这个新点的坐标与能量值，接受的概率为 $e^{\frac{-(E_{t+1}-E_t)}{kT}}$ ，表达式为：
$P=\left{ \begin{aligned} 1\qquad\qquad\quad,E_{t+1} < E_t
e^{\frac{-(E_{t+1}-E_t)}{kT}}, E_{t+1} >= E_t \end{aligned} \right.$ 其中k是接受概率系数，一般取$k = 1$。

Metropolis准则使得在目标点有可能陷入局部最优解B、C、E时，及时跳出，从而可能寻找到更优的全局最优解。

2.3 随机种子选取

在每个新温度or新的迭代次数下，我们究竟是如何完成”跳变”的呢？答案是通过一个随机生成的数字进行前/后移动。随机种子选取准则有两种方法：

2.3.1 Gauss分布(0-1分布)

采用高斯分布，其概率密度函数为

通常μ称为均值或者密度的peak、mode， $σ^2$称为方差。如果一个随机变量$X$服从这一分布，则记作

以均值为1，方差为0.38的高斯分布曲线为例，如下图所示：

服从Gauss分布的数是模型中的各个参数。它们的初次均值为initial(也是一个在定义域内随机生成的数)，之后均值为当前值对应的参数，均方差则为参数上下限差值的1/3。这些参数落在均值附近的概率最大，落在边缘的概率最小。不可避免地，在本项目中定义域是有限的，它们位于[lower, upper]区间，而不是遵从(-∞, +∞)分布。因此制定逻辑：若本次生成的随机数落在参数区间之外，则返回重新生成一个随机数。代码块如下：

#include <iostream>
#include <random>
#include <chrono>
#include <ctime>

double normalDist(double initial, double stdDev, double lower, double upper)
{
    if(initial < lower) return lower;
    if(initial > upper) return upper;
    double res = 0.0;
    unsigned seed = std::chrono::system_clock::now().time_since_epoch.count();
    std::default_system_engine gen(seed);
    std::normal_distribution<double> dist_x(initial, stdDev);
    res = dist(gen);
    if(res < lower || res > upper) return normalDist(initial, stdDev, lower, upper);
    return res;
}

2.3.2 均匀分布

假设x服从[a,b]上的均匀分布，则x的概率密度函数如下： $f(x)=\left{ \begin{aligned} 0\quad,xb\\ \frac{1}{a+b},a≤x≤b \end{aligned} \right.$ 以1000个输入点为例，进行范围为[-1,1]的均匀分布，结果如下所示： ![image](https://github.com/user-attachments/assets/1515cd56-507b-4329-b139-8c720e90a85c)

由于均匀分布不可能出现位于区间之外的随机数，因此不必考虑边界条件。代码块如下：

#include <iostream>
#include <random>
#include <chrono>
#include <ctime>

double realDist(double lower, double upper)
{
    unsigned seed = std::chrono::system_clock::now().time_since_epoch.count();
    std::default_random_engine gen(seed);
    std::uniform_real_distribution<double> dist_x(lower, upper);
    return dist_x(gen);
}

2.4 基于随机种子的参数跳变方法

2.4.1 Fast

2.4.2 Cauchy

2.4.3 BoltzMann

2.5 退火策略

当前使用最广泛的退火策略为快速退火策略，其表达式为： $currT = currT * T_k$ 其中currT为当前温度，Tk为降温系数，Tk∈(0,1)，Tk越大降温速度就越慢。此外，其他可行的降温策略如下： $currT = currT / (T_k + 1)$ $(T_k > 0)$ $currT = currT / log_{10}(T_k + 10)$(T_k > 0)$ $currT = currT * e^{-T_k}$(T_k > 0)$

2.6 评价函数计算方法

2.6.1 RMSE

当前评价计算结果正确性的标准为RMSE(均方根误差)，其计算表达式为： $RMSE = \sqrt{\frac{SSE}{W}} = \sqrt{\frac{1}{W}\sum_{i=1}^{N}{w_iv_i^2}}$ 其中，

$SSE$: 加权平方和
$W$: 总体的总权重（如基于Array Height, Step Height, Cu_Thickness 和Surface Height进行评价函数综合计算，那么W=每个对应观测值不为无理数的参数的$w_i$之和）
$N$: 待比对计算结果的数量（排除了观测值为无理数的情况）
$w_i$: 第i个观测值的权重，可以为≥0的任何值
$v_i$: 第i个仿真值与第i个观测值的误差 RMSE取值范围$[0, +\infin)$，RMSE越小则意味着仿真计算结果与参考结果数值越接近，被认为是越好的计算结果。
2.6.2 CorrCoef

带权重的Pearson相关系数（Pearson Correlation Coefficient），单个参数以及全局的计算方式如下： $CorrCoef = \frac{\sum_{i=1}^{N} w_i(x_i - \bar{x})(y_i - \bar{y})}{\sqrt{\sum_{i=1}^N{} w_i(x_i - \bar{x})^2}\sqrt{\sum_{i=1}^{N} w_i(y_i - \bar{y})^2}}$ 其中，
$N$: 待比对计算结果的数量（排除了观测值为无理数的情况）
$x_i$: 第i个仿真值
$\bar{x}$: 仿真值的加权平均值（排除了对应观测值为无理数的情况）
$y_i$: 第i个观测值
$\bar{y}$: 观测值的加权平均值（排除了观测值为无理数的情况）
$w_i$: 第i个观测值的权重，可以为≥0的任何值

该公式与普通的Pearson相关系数公式类似，不同之处在于每个数据点都乘以了一个权重因子$w_i$。如果所有数据点的权重都相等，则带权重的Pearson相关系数就等于普通的Pearson相关系数。 CorrCoef取值范围$[-1, +1]$，取值为1表示两个变量之间存在完全正向线性关系，取值为-1表示存在完全负向线性关系，而取值为0表示两个变量之间没有线性关系。

2.6.3 R2

带权重的R平方（R-squared）用于衡量一个线性回归模型对数据的拟合程度，考虑了每个数据点的权重。与普通的R平方类似。单个参数以及全局计算带权重的R平方的公式如下： $R^2= 1 - \frac{\sum_{i=1}^N w_i(y_i - x_i)^2}{\sum_{i=1}^N w_i(y_i - \bar{y})^2}$ 其中，

$N$: 待比对计算结果的数量（排除了观测值为无理数的情况）
$x_i$: 第i个仿真值
$y_i$: 第i个观测值
$\bar{y}$: 观测值的加权平均值
$w_i$: 第i个观测值的权重，可以为≥0的任何值公式中的分子表示模型预测值与实际观测值之间的加权残差平方和，分母表示观测值的加权总平方和。带权重的R平方通过比较模型预测值的误差和观测值的误差来评估模型的拟合效果，其中每个误差项都乘以相应的权重因子。带权重的R平方取值范围$(-\infty, 1]$，越接近1，表示模型能够很好地解释观测数据的变异性，而越接近0则表示模型的解释能力较差。
3 模型参数含义介绍

本小节主要介绍在Runset.rst中，用户可调各参数的含义。
T——模拟退火过程当前温度
Tmax——模拟退火过程初始温度
Tmin——模拟退火过程最小温度，若当前温度低于该温度，则完成退火全过程
SA_L——当前温度下执行退火操作的迭代次数，该值可由用户决定
Tk——退火系数，该值也可由用户决定
Learn_Rate——模拟退火学习率。学习率越大，校准参数在退火初期跳动能力越强，参数波动范围就越大
max_stay_counter——保持统一计算结果的最大次数。若max_stay_counter次计算结果都为同一值，则结束计算，完成退火全过程
Param_init——模拟退火初始参数，由随机数产生算法生成，不可由用户决定

模拟退火算法的全过程流图如下： Metropolis判断流图如下图所示：

模拟退火算法