模拟退火算法

一，问题描述

在实际日常中，人们会经常遇到如下问题：在某个给定的定义域X内，求函数 $f(x)$ 对应的最优值。此处以最小值问题举例（最大值问题可以等价转化成最小值问题），形式化为：${\min }_{x\in C}f(x)$

概述：基于迭代跳转到正解的进阶搜索算法
条件：问题的解连续性很强，绘制成函数是连续的，包含可能多峰或者多谷，退火的误差就在求解的谷不是问题空间内的最低谷或峰

要素：
1，温度 $T$，可以转化为控制参数，常见的意义包括步长，可以划定问题解的范围
2，初始温度： $T$₀，要求这个控制参数包括所有的状态空间，避免漏解
3，最值温度： 最终降温的结局低温也就是搜索终止的条件

如果X是离散有限取值，那么可以通过穷取法获得问题的最优解；如果X连续，但$f(x)$是凸的，那可以通过梯度下降等方法获得最优解；如果X连续且$f(x)$非凸，虽说根据已有的近似求解法能够找到问题解，可解是否是最优的还有待考量，很多时候若初始值选择的不好，非常容易陷入局部最优值。

所以为了更有效的跳转最优解，模拟退火应运而生

• 不过模拟退火算法到底是如何模拟金属退火的原理？
  1，主要是将热力学的理论套用到统计学上，将搜寻空间内每一点想象成空气内的分子；
  2，分子的能量，就是它本身的动能；
  3，搜寻空间内的每一点，也像空气分子一样带有“能量”，以表示该点对命题的合适程度。
  4，演算法先以搜寻空间内一个任意点作起始：每一步先选择一个“邻居”，然后再计算从现有位置到达“邻居”的概率。
  5，若概率大于给定的阈值，则跳转到“邻居”；若概率较小，则停留在原位置不动。
  6，随着温度的降低，跳跃越来越不随机，最优解也越来越稳定
  （gif全程需10s，引用自wiki）

二，算法实现

1，温度

• 理解角度：温度，步长或者搜索范围
  1，起始（温度）步长是极大的，可能是真的定义域
  2，终止温度是人为认定的，可以认为是我们接受的误差范围
  3，衰减系数：

取值范围： $(0,1)$
常见取值： 0.9，0.99，0.999
衰减形式： 指数衰减

• 一次模拟退火中温度 t 需要覆盖数据范围的3倍，我觉得是因为这符合正态分布 $3\delta$ 覆盖的绝大部分(大概$96\%$) 数据范围，取得经验值一般不会出错

2，跳转

• 根据一个公理： 计算机科学中，判定总比求解简单
• 每次在邻域内随机一个点，算出其价值函数： $f(new-point)$
• $f(new-point)$ 与现有的价值比较大小，我们取得其变化量 $\delta$

1，根据热力学的原理，在温度为$T$时，出现能量差为 $\delta$ 的降温的概率为 $p(\delta )$，表示为$p(\delta )=\exp (\frac{\delta }{kT})$

基于热力学原理：跳转的概率就是（一个经验值）： $-\exp (\delta /t)$

2，随机化算法收敛到局部最优解的概率是极低的，但是多次随机，获得全局最优的概率会大幅提升

因为是随机化的算法，整个计时函数，有时间的话就多跑跑，最优解还是很好搞的

三，例题

1，二维费马点

在二维平面上有 n 个点，第 i 个点的坐标为 (xi,yi)，请你找出一个点，使得该点到这 n 个点的距离之和最小，该点可以选择在平面中的任意位置，甚至与这 n 个点的位置重合。

#include<iostream>
#include<cmath>
#include<ctime>
#include<cstring>
#include<algorithm>

#define x first
#define y second
using namespace std;

typedef pair<double ,double >pdd;
double ans;
int n;
pdd q[109];

inline int read ()
{
    int x; bool flag = 1; char ch = getchar();
    while(ch>'9'&&ch<'0')if(ch=='-')flag = 0,ch= getchar();
    while(ch<='9'&&ch>='0')x=(x<<1)+(x<<3)+ch-'0',ch = getchar ();
    if(!flag)return ~(x-1);
    else return x;
}

inline double  randx (double l , double r)
{
    return (double)rand() / RAND_MAX * (r - l ) + l;
}

inline double dis(pdd a, pdd b)
{
    return sqrt((a.x - b.x)*(a.x - b.x ) + (a.y - b.y ) * (a.y - b.y ));
}

double clac (pdd k)
{
    double res=0;
    for(int i= 1;i<=n;i ++)res+=dis(k,q[i]);
    ans =min(ans,res);
    return res;
}
void simulate_anneal()
{
    pdd dx (randx(0,10000),randx (0,10000));
    for(double  t = 1e4; t >1e-4 ; t*=0.99)
    {
        pdd dt  (randx (dx.x -t,dx.x + t) , randx (dx.y-t,dx.y+t));
        double dlata= clac(dt)- clac(dx);
        if (exp(-dlata/t) > randx(0, 1) ) dx =dt ;
    }
}

signed main(void)
{
    ans = 1e8+9;
    n=read();
    for(int i = 1 ; i <=n; i ++ )cin>>q[i].x>>q[i].y;
    for(int i= 1; i<=100;i++)simulate_anneal();
    printf ("%.0lf",ans);
}

2，NOIP 2022 方差

3，AHOI/JSOI 2014 保龄球

https://www.luogu.com.cn/problem/P4044

这是一个完全的最优解搜寻问题，实际性很强
这是一个借助交换的调整搜寻最优解的例子，交换在贪心和SA中都是很常见的微扰跳转最优的例子
你很想DP，但是DP 很难，那很方便先搜一搜，那随机化还是比较OK的

#include <iostream>
#include <cstring>
#include <algorithm>
#include <cmath>
#include <ctime>

#define x first
#define y second

using namespace std;

typedef pair<int, int> PII;
const int N = 55;

int n, m;
PII q[N];
int ans;

int calc()
{
    int res = 0;
    for (int i = 0; i < m; i ++ )
    {
        res += q[i].x + q[i].y;
        if (i < n)
        {
            if (q[i].x == 10) res += q[i + 1].x + q[i + 1].y;
            else if (q[i].x + q[i].y == 10)
                res += q[i + 1].x;
        }
    }
    ans = max(ans, res);
    return res;
}

void simulate_anneal()
{
    for (double t = 1e4; t > 1e-4; t *= 0.99)
    {
        int a = rand() % m, b = rand() % m;
        int x = calc();
        swap(q[a], q[b]);
        if (n + (q[n - 1].x == 10) == m)
        {
            int y = calc();
            int delta = y - x;
            if (exp(delta / t) < (double)rand() / RAND_MAX)
                swap(q[a], q[b]);
        }
        else swap(q[a], q[b]);
    }
}

int main()
{
    cin >> n;
    for (int i = 0; i < n; i ++ ) cin >> q[i].x >> q[i].y;
    if (q[n - 1].x == 10) m = n + 1, cin >> q[n].x >> q[n].y;
    else m = n;

    for (int i = 0; i < 100; i ++ ) simulate_anneal();

    cout << ans << endl;
    return 0;
}