poj 2420 A Star not a Tree? 费马点的求解 模拟退火

缘起

学习模拟退火算法. poj 2420 A Star not a Tree?

分析

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
每组数据中给n个点(n<=100)，求平面中一个点使得这个点到n个点的距离之和最小。这个点数学上叫费马点

【输入】
第一行是N，然后是N行，每行给出(x,y), 所有的坐标都是整数，并且在[0, 10000]

【输出】
输出最短距离之和（保留整数部分即可）

【样例输入】
4
0 0
0 10000
10000 10000
10000 0

【样例输出】
28284

【限制】
1s

模拟退火算法（simulate anneal, SA）是一种通用的概率算法. 用于在一个较大的状态空间中寻找最优解. SA 是1983年左右发明的算法. SA是解决TSP这种完全NP问题的有效方法之一.

讲到SA，就不得不提及爬山算法

如上图，我们画了一座山（~QAQ~）. 一只喵酱起始在C点. 它想去山的最高处（即B点）. 然后它发现往右爬山就越爬越高.尝到甜头之后， 它就不断的向右爬啊，爬啊，爬到A点的时候, 它发现不论怎么移动，周遭都会降低山的高度. 遂它下了一个错误的结论——我它喵的已经在最高处啦~ 咳咳咳，井底之喵~ 这就是爬山算法。

但是整个过程如果多一瓶酒结局就可能完全不一样了. 喵喵喝醉了. 当然它也不是那么的不清醒. 如果移动之后高度上升了它还是知道自己正在往上爬的. 它当然会义无反顾的向上爬啦~ 但是，它喝醉了. 有的时候它发现自己在下降的话, 例如从A爬到E的话, 它也会接受的——因为它醉了嘛~ 但是醉意总会消散的——随着时间的推移, 它逐渐清醒了. 例如在D处它慢慢清醒了, 则它渐渐的不接受会让它状态变差的爬行方向了. 则野喵酱将成功的到达真正的最高处B了. 这就是模拟退火算法.

显然, SA和爬山算法的核心基础都是贪心——欣然接受对自己有利的状态转移（喵喵发现下一次移动的状态会让自己的函数值（即山的高度）更大）. 但是不同之处在于， 模拟退火算法会以一定概率（此概率随着时间推移而降低）接受这种让自己的函数值降低的状态转移（即喵喵喝醉了）. 这种接受就很有可能让喵喵跳出（即经过E越过D）爬山算法仅仅找到局部最优而非全局最优解A而找到真正的全局最优解B.

于是我们不难明白为什么这叫做模拟退火. 因为和金属的退火（退火就是一块金属加热到高温然后置于空气中慢慢冷却）十分相似. 温度高的时候金属分子十分活跃. 温度渐渐降低, 分子就慢慢动能降低了. 这里的活跃就是上面说的接受状态变差的概率——温度高的时候分子比较活跃（状态转移比较活跃）就越有可能接受这种变差，而后面温度渐渐降低，分子慢慢不活跃的时候, 就越不可能接受这种变差了.

所以一个自然的问题就是上图，喵酱从左至右接受状态变差（即山的高度降低）的概率越来越低，如果没到E和D之间的低谷概率就已经很低了——导致喵酱很有可能最终过不去这个低谷，就可能致使模拟退火失败——依旧找到的是A这个错误的最优解. 而我们也敏锐的注意到了——这完全取决于我们让这个概率衰减的速率——不能让这个衰减速率太快, 太快了就过不了低谷了.

我们可以给出模拟退火的算法框架

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
T=初始化温度
cur = 初始化状态
while(T>最低温度)
{
   	nxt=cur转移后的的状态
    double delta = E(nxt)-E(cur); // E是状态函数
    if(delta>0) // 借用爬山问题背景, 我们希望E越大越好
    {
        cur = nxt; // 欣然转移
        ans = E(nxt) // 更优解
    }
    else // delta<=0
    {
        bool change = exp(delta/T)>random(0,1); // delta越大（例如0最大,其实表示虽然状态变差了,但其实没有变得太糟）, 则change越可能为真,即越可能接受这种变差,T越大(即退火早期),则接受变差的概率也越大. 这和金属退火是一样的.
        if(change)
        {
            cur = nxt; // ans不变, ans对cur说: "你去开拓嘛,没有更优的,我不变."
        }
    }
    T*=alpha; // alpha是温度衰减因子,一般为了防止概率衰减太快, alpha取0.98或者0.99.若alpha过大，则搜索到全局最优解的概率可能会较高，但搜索的过程也就较长。若alpha过小，则搜索的过程会很快，但最终可能只会搜索到一个局部最优值而非全局最优解
}

不难看出，SA的要素是

  1. 状态如何转移, 这个视实际问题而定.
  2. 状态函数是什么，一般就取实际问题中的目标函数

模拟退火算法的优缺点

  • 迭代搜索效率高，并且可以并行化；
  • 算法中有一定概率接受比当前解较差的解，因此一定程度上可以跳出局部最优, 优于爬山算法；
  • 算法求得的解与初始解状态S无关，因此有一定的鲁棒性；
  • 最重要的，具有渐近收敛性，模拟退火已在理论上被证明是一种依概率收敛于全局最优解的全局优化算法。

好了，来看看本题该怎么做.

本题给的状态空间太大了——10000*10000. 如果枚举，一定T. 所以要用模拟退火. 但是我按照上面的框架写了一版

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
#include "stdafx.h"
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <time.h>
#include <algorithm>
using namespace std;
#define LOCAL
#define SQUARE(x) ((x)*(x))
const int maxn = 105, dir[4][2]={{1,0},{-1,0},{0,1},{0,-1}}; // 四个转移方向
const double eps = 1e-8;
double ans;
int n;

inline double kk()
{
	return (rand()%10000+1.0)/10000.0;
}

struct Point
{
	double x,y;
}p[maxn];

double dis(Point &a, Point &b)
{
	return sqrt(SQUARE(a.x-b.x)+SQUARE(a.y-b.y));
}

double E(Point &fermat) // 状态函数
{
	double ans = 0;
	for (int i = 0; i<n; i++)
	{
		ans+=dis(p[i], fermat);
	}
	return ans;
}

void transform(Point &nxt, Point &cur, int i, double t) // 定义如何转移状态——当前状态cur按照dir[i]方向转移为nxt状态
{
	nxt.x = cur.x + dir[i][0]*t;
	nxt.y = cur.y + dir[i][1]*t;
}

void sa()
{
	double t = 100000, alpha = 0.9999, cure, nxte, delta; // t是初始温度, alpha是温度的衰减速率, cure是当前状态的函数值, nxte是转移后的状态的函数值
	Point cur = p[0], nxt; // 初始状态为p[0], nxt是转移后的状态
	ans = cure = E(cur);
	while(t>eps)
	{
		int i =  rand()%4; // 随机选择方向
		transform(nxt, cur, i, t);
		double nxte = E(nxt);
		if (nxte+eps<cure)
		{
			cur = nxt; // 状态转移
			ans = min(ans, nxte);
		}
		else
		{
			bool change = false; // 要不要转移
			delta = cure-nxte; // 变差了, delta小于0
			double z = kk();
			if (exp(delta/t)>z)
			{
				change = true;
			}
			if (change)
			{
				cur = nxt; // 转移状态, 但是不优化ans(你去开拓嘛~ 没有更优的, 我ans不变)
			}
		}
		t*=alpha;
	}
}

int main()
{
#ifdef LOCAL
	freopen("d:\\data.in", "r", stdin);
//	freopen("d:\\my.out", "w", stdout);
#endif
	srand((int)(time(0)));
	scanf("%d", &n);
	for (int i = 0; i<n; i++)
	{
		scanf("%lf%lf", &p[i].x, &p[i].y);
	}
	ans = 0x3f3f3f3f;
	sa(); // 跑模拟退火算法
	printf("%.0lf\n", ans);
	return 0;
}

但是很遗憾, 效率太低了, 所以找到网上的题解（思想依旧是上面的思想哈~） 其实大部分OIer都把模拟退火叫做“随机化的步长递减贪心”. 随机化体现在探索方向上——无规律可循的

算法很简单, 就是每次探索的步长固定. 然后沿着上下左右四个方向（甚至8个方向）探索. 如果能找到更好的, 则下一轮还是这个步长. 否则，将步长/2, 直至步长几乎为0即可.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
//#include "stdafx.h"
#include <stdio.h>
#include <math.h>
//#define LOCAL
#define SQUARE(x) ((x)*(x))
const int maxn = 105, dir[4][2] = {{0,1},{0,-1},{1,0},{-1,0}};
const double inf = 0x3f3f3f3f, eps = 1e-8;
int n;
struct Point
{
	double x,y;
}p[maxn], cur, nxt;
double ans;

void transform(int i, double step) // i方向step改变当前状态
{
	nxt.x = cur.x+dir[i][0]*step;
	nxt.y = cur.y+dir[i][1]*step;
}

double dis(Point &a, Point &b)
{
	return sqrt(SQUARE(a.x-b.x)+SQUARE(a.y-b.y));
}

double e(Point &fermat) // 状态函数
{
	double ans = 0;
	for (int i = 0; i<n; i++)
	{
		ans+= dis(p[i], fermat);
	}
	return ans;
}

void sa() // 模拟退火
{
	double step = 100; // 初始步长
	cur = p[0], ans = e(p[0]);
	while(step>eps)
	{
		bool ok = true; // 此轮有没有找到更优的?
		while(ok)
		{
			ok = false; // 默认没找到
			for (int i = 0; i<4; i++) // 沿着四个方向走四步
			{
				transform(i, step); // 状态转移
				double nxte = e(nxt); // 计算当前状态和转移后的状态
				if (ans>eps+nxte) // 如果此轮（一轮指的是走4步）找到更优的了
				{
					ans = nxte;
					cur = nxt; // 状态转移
					ok = true; // 将继续使用step这个步长而不缩小
				}
			}
		}
		step/=2;
	}
}

int main()
{
#ifdef LOCAL
	freopen("d:\\data.in", "r", stdin);
//	freopen("d:\\my.out", "w", stdout);
#endif
	scanf("%d", &n);
	for (int i = 0; i<n; i++)
	{
		scanf("%lf%lf", &p[i].x, &p[i].y);
	}
	sa();
	printf("%.0f\n", ans+eps); // 加上eps是为了防止丢失精度,WA到哭~ 而且这里要用 f 不能使用 lf, 因为lf为双精度输出，有15-16位有效数字，f为单精度输出，有6-7位有效数字。而这里的eps是1e-8
	return 0;
}

ac情况

Status Accepted
Memory 364kB
Length 1378
Lang G++
Submitted 2019-09-27 22:27:37
Shared
RemoteRunId 20905451

通过本题我们知道了——模拟退火适宜什么类型的问题: 搜索空间极大的最优解问题, 几乎不能通过枚举来搞定.