模拟退火学习报告（P3878 [TJOI2010]分金币）

P3878 [TJOI2010]分金币

题目描述

现在有n枚金币，它们可能会有不同的价值，现在要把它们分成两部分，要求这两部分金币数目之差不超过1，问这样分成的两部分金币的价值之差最小是多少？

题目分析

分析发现，我们可以每次交换在两个部分中的两个元素，然后求两部分差最小的情况.

但是观察$n$的范围不大于$30$，那么在数据最强的情况下，搜索的复杂度可以达到$O(2^n)$的水平，显然是承受不了的.

又因为这是一个最优解问题，所以我们考虑模拟退火算法.

模拟退火算法

模拟退火算法来源于固体退火原理，将固体加温至充分高，再让其徐徐冷却，加温时，固体内部粒子随温升变为无序状，内能增大，而徐徐冷却时粒子渐趋有 序，在每个温度都达到平衡态，最后在常温时达到基态，内能减为最小。根据Metropolis准则，粒子在温度T时趋于平衡的概率为 e(-ΔE/(kT))，其中E为温度T时的内能，ΔE为其改变量，k为Boltzmann常数。用固体退火模拟组合优化问题，将内能E模拟为目标函数值 f，温度T演化成控制参数t，即得到解组合优化问题的模拟退火算法：由初始解i和控制参数初值t开始，对当前解重复“产生新解→计算目标函数差→接受或舍弃”的迭代，并逐步衰减t值，算法终止时的当前解即为所得近似最优解，这是基于蒙特卡罗迭代求解法的一种启发式随机搜索过程。退火过程由冷却进度表 (Cooling Schedule)控制，包括控制参数的初值t及其衰减因子Δt、每个t值时的迭代次数L和停止条件S。

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

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

我们还可以用一个形象的比喻来描述模拟退火算法：兔子喝醉了。它随机地跳了很长时间. 这期间，它可能走向高处，也可能踏入平地. 但是，它渐渐清醒了并朝最高方向跳去. 这就是模拟退火.

由上述可以发现，模拟退火大致有这些步骤：

1.设置初始温度T，初始符合条件的答案；
2.通过某种神奇的方式，找到另一个符合条件的新状态；
3.分别将两个状态的答案计算出来，并作差得到$\Delta E$；
4.根据题目要求，贪心的决定是否更换答案（有些题目是改变状态）. 即：选择最优解；如果无法替换答案，则根据一定概率替换答案. 即运用到上述的平衡概率$exp(\Delta E/T)$（表示$e$的$\Delta E/T$次方）随机的决定是否替换.
每一次操作后，进行降温操作。即：将温度$T$乘上某一个系数，一般是$0.985-0.999$随具体题目随缘定.

我们会发现，新状态离最优解越近，温度越高，越容易被接纳

伪代码实现如下:

esp=1e-15;//终止温度
T=初始温度;
while(T>esp)
{
    now=从当前最优状态随机更新的一个状态；
    delta=calc(now)-calc(ans);
    if(delta与题目要求的满足更优)ans=now;
    else if(exp(delta/T)*RAND_MAX>rand())ans=now;//PS:这里的delta前面可能要加'-'；也可能不是对最优解的改变而是对状态的改变（限于最优解记录状态的时候）
    T*=t0;//t0一般在0.985-0.999之间，根据具体题目时间，随缘调试。。。
}

关于上述delta的符号问题：

自然对数的底$e$的幂必须为负数，所以有时候delta要加负号，有时候不用.

---

以上，我们可以假定状态是当前的差值，每次更新状态时就交换一组元素即可.

用ans记录最优解，所以ans的值不改变，改变的是状态.

程序实现

#include<bits/stdc++.h>
#define esp 1e-8//escape温度也是随缘
#define dt 0.998
#define ll long long
using namespace std;
ll tot1,tot2,a[40],ans;
int n,len1,len2;
void SA(){
	double T=1200;
	while(T>esp){
		int x=(int )(rand()%len1+1),y=(int )(rand()%len2+len1+1);
		ll delta=abs((tot1+a[y]-a[x])-(tot2+a[x]-a[y]))-ans;//随机新状态
		if(delta<0){
			tot1=tot1+a[y]-a[x];
			tot2=tot2+a[x]-a[y];
			ans=abs(tot1-tot2);
			swap(a[x],a[y]);
		}
		else if(exp((-delta)/T)*RAND_MAX>rand()){//RAND_MAX表示rand函数能够取得到的最大值，这是对是否更新状态的估值；exp表示自然对数的底的k次方
			tot1=tot1+a[y]-a[x];
			tot2=tot2+a[x]-a[y];
//			ans=abs(tot1-tot2);此行有误
			swap(a[x],a[y]);//改变的是状态，不是用来记录最优解的ans 
		}
		T*=dt;
	}
}
int main(){
//	freopen("make.in","r",stdin);
//	freopen("1.out","w",stdout);
	srand((int )time(0));//随机时间种子
	int T;
	scanf("%d",&T);
	for(int i=1;i<=T;i++){
		len1=0,len2=0;
		memset(a,0,sizeof a);
		tot1=0,tot2=0,ans=0;//没初始化就完了
		scanf("%d",&n);
		len1=n/2;
		len2=n-len1;
		for(int i=1;i<=n;i++){
			scanf("%lld",&a[i]);
			if(i<=len1)tot1+=a[i];
			else if(i>len1)tot2+=a[i];
		}
		ans=abs(tot1-tot2);
		if(n==1){printf("%lld\n",a[1]);continue;}
		for(int i=1;i<=20;i++)SA();//SA是模拟退火英文Simulate Anneal的缩写
		printf("%lld\n",ans);
	}
	return 0;
} 

感谢算法支持：

SA算法介绍本家

rand()和srand()的用法