模拟退火，参数难调，~~急需欧皇~~
——《AK杀》

前言

今天（2019/6/19），$Payphone-X$想要学习最小生成树，结果翻到了模拟退火……

于是乎

就有了今天这篇博客

何为模拟退火

先来看看百度百科对于模拟退火的解释

模拟退火的出发点是基于物理中固体物质的退火过程与一般组合优化问题之间的相似性。模拟退火算法是一种通用的优化算法，其物理退火过程由加温过程、等温过程、冷却过程这三部分组成。

——百度百科

嗯……

看起来有点玄学的气息……

还是来看看通俗的解释吧。

简单的说，模拟退火是一种随机化搜索算法。~~（所有除了机惨的算法里需要用到rand与random_shuffle的算法~~

它可以求得一段多峰函数的极值~~（也可能求不出来~~

也就是说，当我们面对一个方案数极大（甚至是无穷大）的题目时，可以使用模拟退火

~~是一个比较看脸的算法~~

爬山算法

在讲模拟退火之前，我们先来看看它的基础算法——爬山算法

爬山算法是一种基于贪心的暴力搜索，可以求得一个单峰函数的极值。

它的原理是每一次寻找一个比当前更优的解，一旦找不出，就结束。

就像这样：

当爬到顶峰，即无法继续向上爬时，即停止搜索。就像这样：

嗯……

现在为止，一切好像都很正常

但……如果整张图是这样的

就自闭了……

而模拟退火会在寻找到一个非局部最优解时，赋予一个接受他的概率。从而继续向下寻找，直到找到全局最优解。
（ps.即在爬山时，模拟退火会有一定概率下降当前高度去寻找是否可以有别的路到达更高的高度。）

算法原理

要想了解模拟退火，就必须了解金属的退火原理

简单来说，就是将金属加热到一定温度，保持足够时间，然后以适宜速度冷却。

而OI上，就是每次随机出一个新解。如果这个解更优，就接受它，否则以一个与温度和与最优解的差相关的概率接受它。

设这个新的解与最优解的差为 $\Delta E$ ，温度为 $T$ ， $k$ 为一个随机数，那么这个概率为 $\displaystyle e^{\frac{\Delta E}{kT}}$

算法实现

降温

模拟退火时有3个参数，分别为初始温度$T_0$，降温系数$Δ$，终止温度$T_k$。

其中，$T_0$的值会很大，$Δ$是一个略小于$1$的正数，$T_k$是一个略大于$0$的数

我们先把初始温度$T_0$的值赋给$T$，每一次降温时让$T = T * Δ$，直到$T ≤ T_k$为止

伪代码如下

double start = 120, s = 0.96 , finish = 15;
double tem = start;
while(tem > finish){
    tem = tem * s;
}

如果转换成图片，大概是这样的

可以看出，随着温度的不断降低，解也不断稳定下来，并集中在最优解附近。

others

不多说，直接看表。

首先，我们在程序开始时需要srand(一个常数)来随机化，这个常数可以是任意的常数，比如$23333$ , $1212121$ ~~或古老东方的某个神秘的八位信仰质数~~

一遍退火一般不能直接得到最优解，所以要在不TLE的情况下多跑几遍。

模拟退火的时间复杂度比较玄学，一般降温系数 $Δ$ 与 $1$ 的差减少一个数量级, 耗时大约多 $10$ 倍；但$T_0$ 和 $T_k$ 如果变化一个数量级, 耗时却不会变化很大。

最近有的出题人卡常数随机化，针对此种状况，我们可以 srand((unsigned)time(NULL));

调参

众所周知，调参是一门~~极其玄学~~的艺术

不信可以看我的提交记录：

在调参过程中，不同的$Δ$，$T_0$，$T_k$，$srand$常数甚至退火的次数都会影响最后的答案。

下面，我们就来探讨一下如何调参

Q：答案不是最优的怎么办？

A：以下几种方法：调大 $Δ$ ，调大 $T_0$ ，调小 $T_k$ ，以及多跑几遍。其中，如果您的程序跑的比较快，可以尝试多跑几遍或调大$Δ$，调大 $T_0$ 与调小 $T_k$对运行时间影响较小，但可能效果不会很明显。

Q：还是跑不出最优解怎么办？

A：可以尝试更换随机数种子，或者 srand(rand())srand(rand());总有可能跑过。

Q：那我提交了两页都没有$AC$，咋办啊？

A：~~那您还是写正解吧~~

例题

现在，我们以这道题为例来讲解模拟退火的应用。

题目描述

如图：有$n$个重物，每个重物系在一条足够长的绳子上。

每条绳子自上而下穿过桌面上的洞，然后系在一起，图中$X$处就是公共的绳结。

假设绳子是完全弹性的（不会造成能量损失），桌子足够高（因而重物不会垂到地上），且忽略所有的摩擦。

问绳结$X$最终平衡于何处。

（注意：桌面上的洞都比绳结$X$小得多，所以即使某个重物特别重，绳结$X$也不可能穿过桌面上的洞掉下来，最多是卡在某个洞口处。）

输入格式：

第一行为一个正整数$n$，表示重物和洞的数目。

接下来的$n$行，每行是3个整数：$X_i.Y_i.W_i$，分别表示第$i$个洞的坐标以及第 $i$个重物的重量。

输出格式：

输出两个浮点数（保留小数点后三位）分别表示处于最终平衡状态时绳结$X$的横坐标和纵坐标。两个数以一个空格隔开。

Input & Output ‘s examples

Input ‘s eg

3
0 0 1
0 2 1
1 1 1

Output ‘s eg

0.577 1.000

数据范围和约定

对于$100\%$数据，$1≤n≤1000$，$-10000≤x,y≤10000, 0< w ≤ 1000$

分析

恐怕这是OI上为数不多的物理题了……~~（看不懂题的话先补物理去吧）~~

根据初中物理的知识，当系统处于平衡状态时，系统的总能量最小。~~（有人说初中不学这玩意，但老宋确实讲了，而且创客转OI的同学们应该觉得这是基础~~

又因为物体静止,不用考虑动能与弹性势能，则此时系统的总能量是等于各个物体的重力势能之和。

由公式$E_重 = mgh$可得，物体的重力势能取决于物体的质量与高度，

由于每个物体的质量是固定的（题目给出），我们能调整的，只有物体的高度。

为了让系统的重力势能尽可能小，我们需要尽量每一个物体的高度更小，即使每一个$h$更小，也就是让桌子下面的绳子长度更长。

也就是说，我们要让在桌子上边的绳子长度尽量的短。

当桌子上边的绳子长度最短时，就达到了平衡状态。

所以我们用模拟退火跑出桌上绳子的最短长度，然后找出此时绳结的坐标即可。

附件[AC标程]

/*#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#include<iomanip>*/
#include<bits/stdcpp.h>
#define N 1005
#define debug(x){\
	printf("----------------------------\n");\
	printf("DEBUG: There is no bug /n",x);\
}
#define dai 0.996		//降温系数
#define C 4				//模拟退火的执行次数

using namespace std;

int n;

struct Things{
	int x;	//物体横坐标
	int y;	//物体纵坐标
	int w;	//物体重量
}t[N];     

double ansx , ansy , answ;	//答案

double energy(double x , double y){		//公式模拟
	double r = 0;
	double dx , dy;
	for(int i = 1; i <= n; i ++){
		dx = x - t[i].x;
		dy = y - t[i].y;
		r = r + sqrt(dx * dx + dy * dy) * t[i].w;
	}
	return r;
}

void SA(){		//模拟退火
	double tem = 3000;		//初始温度调高
	while(tem > 1e-15){
		double ex = ansx + (rand() * 2 - RAND_MAX) * tem;
		double ey = ansy + (rand() * 2 - RAND_MAX) * tem;
		double ew = energy(ex , ey);
		double ac = ew - answ;
		if(ac < 0){			//如果答案比当前答案更优，则接受
			ansx = ex;
			ansy = ey;
			answ = ew;
		}
		else if(exp (-(ew - answ) / tem) * RAND_MAX > rand()){
			ansx = ex;
			ansy = ey;		//否则以一定概率接受
		}
		tem = dai * tem;
	}
}

void work(){				//多跑几遍才能跑出最优解
	for(int i = 1; i <= C; i ++){
        SA();
    }
}


int main(){
	cin >> n;
	for(int i = 1; i <= n; i ++){
		cin >> t[i].x >> t[i].y >> t[i].w;
		ansx = ansx + t[i].x;
		ansy = ansy + t[i].y;
	}
	ansx = ansx / n;		//取平均值作为初始答案
	ansy = ansy / n;
	answ = energy(ansx , ansy);
	work();
	cout << fixed << setprecision(3) << ansx <<" "<< ansy << endl;
	return 0;
}