模拟退火算法来源于固体退火原理，将固体加温至充分高，再让其徐徐冷却，加温时，固体内部粒子随温升变为无序状，内能增大，而徐徐冷却时粒子渐趋有序，在每个温度都达到平衡态，最后在常温时达到基态，内能减为最小。根据 Metropolis 准则，粒子在温度T时趋于平衡的概率为 e-ΔE/(kT) ，其中 E 为温度 T 时的内能， ΔE 为其改变量， k 为 Boltzmann 常数。用固体退火模拟组合优化问题，将内能E模拟为目标函数值 f ，温度 T 演化成控制参数 t ，即得到解组合优化问题的模拟退火算法：由初始解 i 和控制参数初值 t 开始，对当前解重复“产生新解计算目标函数差接受或舍弃”的迭代，并逐步衰减t值，算法终止时的当前解即为所得近似最优解，这是基于蒙特卡罗迭代求解法的一种启发式随机搜索过程。退火过程由冷却进度表 (Cooling Schedule) 控制，包括控制参数的初值 t 及其衰减因子 Δt 、每个 t 值时的迭代次数 L 和停止条件 S 。
模拟退火算法的模型
模拟退火算法可以分解为解空间、目标函数和初始解三部分。
模拟退火的基本思想:
(1) 初始化：初始温度 T (充分大)，初始解状态 S (是算法迭代的起点)， 每个 T 值的迭代次数 L
(2) 对 k=1，……，L 做第 3 至第 6 步：
(3) 产生新解 S'
(4) 计算增量 Δt′=C(S′)-C(S) ，其中 C(S) 为评价函数
(5) 若 Δt′<0 则接受 S′ 作为新的当前解，否则以概率 exp(-Δt′/T) 接受 S′ 作为新的当前解.
(6) 如果满足终止条件则输出当前解作为最优解，结束程序。 终止条件通常取为连续若干个新解都没有被接受时终止算法。
(7) T 逐渐减少，且 T->0 ，然后转第 2 步。
模拟退火的算法流程图如下:
模拟退火算法新解的产生和接受可分为如下四个步骤：
第一步是由一个产生函数从当前解产生一个位于解空间的新解；为便于后续的计算和接受，减少算法耗时，通常选择由当前新解经过简单地变换即可产生新解的方法，如对构成新解的全部或部分元素进行置换、互换等，注意到产生新解的变换方法决定了当前新解的邻域结构，因而对冷却进度表的选取有一定的影响。
第二步是计算与新解所对应的目标函数差。因为目标函数差仅由变换部分产生，所以目标函数差的计算最好按增量计算。事实表明，对大多数应用而言，这是计算目标函数差的最快方法。
第三步是判断新解是否被接受,判断的依据是一个接受准则，最常用的接受准则是 Metropolis 准则: 若 Δt′<0 则接受 S′ 作为新的当前解 S ，否则以概率 exp(-Δt′/T) 接受 S′ 作为新的当前解 S 。
第四步是当新解被确定接受时，用新解代替当前解，这只需将当前解中对应于产生新解时的变换部分予以实现，同时修正目标函数值即可。此时，当前解实现了一次迭代。可在此基础上开始下一轮试验。而当新解被判定为舍弃时，则在原当前解的基础上继续下一轮试验。模拟退火算法与初始值无关，算法求得的解与初始解状态 S (是算法迭代的起点)无关；模拟退火算法具有渐近收敛性，已在理论上被证明是一种以概率 l 收敛于全局最优解的全局优化算法；模拟退火算法具有并行性
如果你对退火的物理意义还是晕晕的，没关系我们还有更为简单的理解方式。想象一下如果我们现在有下面这样一个函数，现在想求函数的(全局)最优解。如果采用 Greedy 策略，那么从 A 点开始试探，如果函数值继续减少，那么试探过程就会继续。而当到达点B时，显然我们的探求过程就结束了(因为无论朝哪个方向努力，结果只会越来越大)。最终我们只能找打一个局部最后解 B 。
模拟退火其实也是一种 Greedy 算法，但是它的搜索过程引入了随机因素。模拟退火算法以一定的概率来接受一个比当前解要差的解，因此有可能会跳出这个局部的最优解，达到全局的最优解。以上图为例，模拟退火算法在搜索到局部最优解 B 后，会以一定的概率接受向右继续移动。也许经过几次这样的不是局部最优的移动后会到达 B 和 C 之间的峰点，于是就跳出了局部最小值B。
根据Metropolis准则，粒子在温度T时趋于平衡的概率为 exp(-ΔE/(kT)) ，其中 E 为温度 T 时的内能，ΔE为其改变数, k 为 Boltzmann 常数。 Metropolis 准则常表示为
Metropolis 准则表明，在温度为 T 时，出现能量差为 dE 的降温的概率为 P(dE)，表示为： P(dE) = exp(dE/(kT)) 。其中 k 是一个常数， exp 表示自然指数，且 dE<0 。所以 P 和 T 正相关。这条公式就表示：温度越高，出现一次能量差为 dE 的降温的概率就越大；温度越低，则出现降温的概率就越小。又由于 dE 总是小于 0 (因为退火的过程是温度逐渐下降的过程)，因此 dE/kT < 0 ，所以 P(dE) 的函数取值范围是 (0,1) 。随着温度 T 的降低， P(dE) 会逐渐降低。
我们将一次向较差解的移动看做一次温度跳变过程，我们以概率 P(dE) 来接受这样的移动。也就是说，在用固体退火模拟组合优化问题，将内能E模拟为目标函数值 f，温度T演化成控制参数 t，即得到解组合优化问题的模拟退火演算法：由初始解 i 和控制参数初值 t 开始，对当前解重复“产生新解计算目标函数差接受或丢弃”的迭代，并逐步衰减 t 值，算法终止时的当前解即为所得近似最优解，这是基于蒙特卡罗迭代求解法的一种启发式随机搜索过程。退火过程由冷却进度表 (Cooling Schedule) 控制，包括控制参数的初值 t 及其衰减因子 Δt 、每个 t 值时的迭代次数 L 和停止条件 S 。
总结起来就是：
若 f( Y(i+1) ) <= f( Y(i) ) (即移动后得到更优解)，则总是接受该移动；
若 f( Y(i+1) ) > f( Y(i) ) (即移动后的解比当前解要差)，则以一定的概率接受移动，而且这个概率随着时间推移逐渐降低(逐渐降低才能趋向稳定)相当于上图中，从 B 移向 BC 之间的小波峰时，每次右移(即接受一个更糟糕值)的概率在逐渐降低。如果这个坡特别长，那么很有可能最终我们并不会翻过这个坡。如果它不太长，这很有可能会翻过它，这取决于衰减 t 值的设定。
关于普通 Greedy 算法与模拟退火，有一个有趣的比喻：
普通 Greedy 算法：兔子朝着比现在低的地方跳去。它找到了不远处的最低的山谷。但是这座山谷不一定最低的。这就是普通 Greedy 算法，它不能保证局部最优值就是全局最优值。
模拟退火：兔子喝醉了。它随机地跳了很长时间。这期间，它可能走向低处，也可能踏入平地。但是，它渐渐清醒了并朝最低的方向跳去。这就是模拟退火。
模拟退火算法的简单应用
作为模拟退火算法应用，讨论货郎担问题 ( Travelling Salesman Problem ，简记为 TSP )：设有 n 个城市，用数码1,…,n 代表。城市i和城市j之间的距离为 d(i，j) i, j=1,…,n ．TSP 问题是要找遍访每个域市恰好一次的一条回路，且其路径总长度为最短.。
求解 TSP 的模拟退火算法模型可描述如下：
解空间 解空间S是遍访每个城市恰好一次的所有回路，是 {1，……，n} 的所有循环排列的集合， S 中的成员记为(w1,w2 ,……，wn) ，并记 wn+1= w1 。初始解可选为 (1，……，n)
目标函数 此时的目标函数即为访问所有城市的路径总长度或称为代价函数：
我们要求此代价函数的最小值。
新解的产生 随机产生 1 和 n 之间的两相异数 k 和 m ，若 k
(w1, w2 ,…，wk , wk+1 ,…，wm ,…，wn)
变为：
(w1, w2 ,…，wm , wm-1 ,…，wk+1 , wk ,…，wn) .
如果是 k>m ，则将
(w1, w2 ,…，wk , wk+1 ,…，wm ,…，wn)
变为：
(wm, wm-1 ,…，w1 , wm+1 ,…，wk-1 ,wn , wn-1 ,…，wk) .
上述变换方法可简单说成是“逆转中间或者逆转两端”。
也可以采用其他的变换方法，有些变换有独特的优越性，有时也将它们交替使用，得到一种更好方法。
代价函数差 设将 (w1, w2 ,……，wn) 变换为 (u1, u2 ,……，un) , 则代价函数差为：
根据上述分析，可写出用模拟退火算法求解 TSP 问题的伪程序：
Procedure TSPSA:
begin
init-of-T; { T为初始温度}
S={1，……，n}; {S为初始值}
termination=false;
while termination=false
begin
for i=1 to L do
begin
generate(S′form S); { 从当前回路S产生新回路S′}
Δt:=f(S′))-f(S);{f(S)为路径总长}
IF(Δt<0) OR (EXP(-Δt/T)>Random-of-[0,1])
S=S′;
IF the-halt-condition-is-TRUE THEN
termination=true;
End;
T_lower;
End;
End
下面给出C++实现参考源码:
/*
模拟退火算法解决TSP问题
输入格式(tsp.in)：
第1行:1个整数N，表示城市的数量
第2..N+1行：每行有2个空格分开的整数x,y，第i+1行的x,y表示城市i的坐标
*/
#include
#include
#include
#include
#include
#include
#include
#define N 30//城市数量
#define T 3000//初始温度
#define EPS 1e-8//终止温度
#define DELTA 0.98//温度衰减率
#define LIMIT 1000//概率选择上限
#define OLOOP 20//外循环次数
#define ILOOP 100//内循环次数
using namespace std;
//定义路线结构体
struct Path
{
int citys[N];
double len;
};
//定义城市点坐标
struct Point
{
double x, y;
};
Path bestPath; //记录最优路径
Point p[N]; //每个城市的坐标
double w[N][N]; //两两城市之间路径长度
int nCase; //测试次数
double dist(Point A, Point B)
{
return sqrt((A.x - B.x) * (A.x - B.x) + (A.y - B.y) * (A.y - B.y));
}
void GetDist(Point p[], int n)
{
for(int i = 0; i < n; i++)
for(int j = i + 1; j < n; j++)
w[i][j] = w[j][i] = dist(p[i], p[j]);
}
void Input(Point p[], int &n)
{
scanf("%d", &n);
for(int i = 0; i < n; i++)
scanf("%lf %lf", &p[i].x, &p[i].y);
}
void Init(int n)
{
nCase = 0;
bestPath.len = 0;
for(int i = 0; i < n; i++)
{
bestPath.citys[i] = i;
if(i != n - 1)
{
printf("%d--->", i);
bestPath.len += w[i][i + 1];
}
else
printf("%d\n", i);
}
printf("\nInit path length is : %.3lf\n", bestPath.len);
printf("-----------------------------------\n\n");
}
void Print(Path t, int n)
{
printf("Path is : ");
for(int i = 0; i < n; i++)
{
if(i != n - 1)
printf("%d-->", t.citys[i]);
else
printf("%d\n", t.citys[i]);
}
printf("\nThe path length is : %.3lf\n", t.len);
printf("-----------------------------------\n\n");
}
Path GetNext(Path p, int n)
{
Path ans = p;
int x = (int)(n * (rand() / (RAND_MAX + 1.0)));
int y = (int)(n * (rand() / (RAND_MAX + 1.0)));
while(x == y)
{
x = (int)(n * (rand() / (RAND_MAX + 1.0)));
y = (int)(n * (rand() / (RAND_MAX + 1.0)));
}
swap(ans.citys[x], ans.citys[y]);
ans.len = 0;
for(int i = 0; i < n - 1; i++)
ans.len += w[ans.citys[i]][ans.citys[i + 1]];
cout << "nCase = " << nCase << endl;
Print(ans, n);
nCase++;
return ans;
}
void SA(int n)
{
double t = T;
srand((unsigned)(time(NULL)));
Path curPath = bestPath;
Path newPath = bestPath;
int P_L = 0;
int P_F = 0;
while(1) //外循环，主要更新参数t，模拟退火过程
{
for(int i = 0; i < ILOOP; i++) //内循环，寻找在一定温度下的最优值
{
newPath = GetNext(curPath, n);
double dE = newPath.len - curPath.len;
if(dE < 0) //如果找到更优值，直接更新
{
curPath = newPath;
P_L = 0;
P_F = 0;
}
else
{
double rd = rand() / (RAND_MAX + 1.0);
//如果找到比当前更差的解，以一定概率接受该解，并且这个概率会越来越小
if(exp(dE / t) > rd && exp(dE / t) < 1)
curPath = newPath;
P_L++;
}
if(P_L > LIMIT)
{
P_F++;
break;
}
}
if(curPath.len < bestPath.len)
bestPath = curPath;
if(P_F > OLOOP || t < EPS)
break;
t *= DELTA;
}
}
int main(int argc, const char * argv[]) {
freopen("TSP.data", "r", stdin);
int n;
Input(p, n);
GetDist(p, n);
Init(n);
SA(n);
Print(bestPath, n);
printf("Total test times is : %d\n", nCase);
return 0;
}
TSP.data 的数据格式如下，第一行的数字表示一个有多少座城市，第 2 至最后一行，每行有两个数字表示，城市的坐标(平面直角坐标系)。例如：
6
20 80
16 84
23 66
62 90
11 9
35 28
注意由于是基于蒙特卡洛的方法，所以上面代码每次得出的结果并不完全一致。你可以通过增加迭代的次数来获得一个更优的结果。
我们这里需要说明的是，在之前的文章里，我们用求最小值的例子来解释模拟退火的执行：如果新一轮的计算结果更前一轮之结果更小，那么我们就接受它，否则就以一个概率来拒绝或接受它，而这个拒绝的概率会随着温度的降低(也即是迭代次数的增加)而变大(也就是接受的概率会越来越小)。
但现在我们面对一个 TSP 问题，我们如何定义或者说如何获取下一轮将要被考察的哈密尔顿路径呢？在一元函数最小值的例子中，下一轮就是指向左或者向右移动一小段距离。而在 TSP 问题中，我们可以采用的方式其实是很多的。上面代码中 GetNext() 函数所采用的方式是随机交换两个城市在路径中的顺序。例如当前路径为 A->B->C->D->A ，那么下一次路径就可能是 A->D->C->B->A ，即交换 B 和 D 。
public class Tour{
... ...
// Creates a random individual
public void generateIndividual() {
// Loop through all our destination cities and add them to our tour
for (int cityIndex = 0; cityIndex < TourManager.numberOfCities(); cityIndex++) {
setCity(cityIndex, TourManager.getCity(cityIndex));
}
// Randomly reorder the tour
Collections.shuffle(tour);
}
... ...
}
可见把上一轮路径做一个随机的重排(这显然也是一种策略)。
我们对上述问题提出一种新的策略:
首先，我们需要创建一个城市类，它可以用来为旅行推销员的不同目的地建模。
/*
* City.java
* Models a city
*/
package sa;
public class City {
int x;
int y;
// Constructs a randomly placed city
public City(){
this.x = (int)(Math.random()*200);
this.y = (int)(Math.random()*200);
}
// Constructs a city at chosen x, y location
public City(int x, int y){
this.x = x;
this.y = y;
}
// Gets city's x coordinate
public int getX(){
return this.x;
}
// Gets city's y coordinate
public int getY(){
return this.y;
}
// Gets the distance to given city
public double distanceTo(City city){
int xDistance = Math.abs(getX() - city.getX());
int yDistance = Math.abs(getY() - city.getY());
double distance = Math.sqrt( (xDistance*xDistance) + (yDistance*yDistance) );
return distance;
}
@Override
public String toString(){
return getX()+", "+getY();
}
}
接下来让我们创建一个可以跟踪城市的类：
/*
* TourManager.java
* Holds the cities of a tour
*/
package sa;
import java.util.ArrayList;
public class TourManager {
// Holds our cities
private static ArrayList destinationCities = new ArrayList();
// Adds a destination city
public static void addCity(City city) {
destinationCities.add(city);
}
// Get a city
public static City getCity(int index){
return (City)destinationCities.get(index);
}
// Get the number of destination cities
public static int numberOfCities(){
return destinationCities.size();
}
}
现在来创建一个可以模拟旅行推销员之旅：
/*
* Tour.java
* Stores a candidate tour through all cities
*/
package sa;
import java.util.ArrayList;
import java.util.Collections;
public class Tour{
// Holds our tour of cities
private ArrayList tour = new ArrayList();
// Cache
private int distance = 0;
// Constructs a blank tour
public Tour(){
for (int i = 0; i < TourManager.numberOfCities(); i++) {
tour.add(null);
}
}
// Constructs a tour from another tour
public Tour(ArrayList tour){
this.tour = (ArrayList) tour.clone();
}
// Returns tour information
public ArrayList getTour(){
return tour;
}
// Creates a random individual
public void generateIndividual() {
// Loop through all our destination cities and add them to our tour
for (int cityIndex = 0; cityIndex < TourManager.numberOfCities(); cityIndex++) {
setCity(cityIndex, TourManager.getCity(cityIndex));
}
// Randomly reorder the tour
Collections.shuffle(tour);
}
// Gets a city from the tour
public City getCity(int tourPosition) {
return (City)tour.get(tourPosition);
}
// Sets a city in a certain position within a tour
public void setCity(int tourPosition, City city) {
tour.set(tourPosition, city);
// If the tours been altered we need to reset the fitness and distance
distance = 0;
}
// Gets the total distance of the tour
public int getDistance(){
if (distance == 0) {
int tourDistance = 0;
// Loop through our tour's cities
for (int cityIndex=0; cityIndex < tourSize(); cityIndex++) {
// Get city we're traveling from
City fromCity = getCity(cityIndex);
// City we're traveling to
City destinationCity;
// Check we're not on our tour's last city, if we are set our
// tour's final destination city to our starting city
if(cityIndex+1 < tourSize()){
destinationCity = getCity(cityIndex+1);
}
else{
destinationCity = getCity(0);
}
// Get the distance between the two cities
tourDistance += fromCity.distanceTo(destinationCity);
}
distance = tourDistance;
}
return distance;
}
// Get number of cities on our tour
public int tourSize() {
return tour.size();
}
@Override
public String toString() {
String geneString = "|";
for (int i = 0; i < tourSize(); i++) {
geneString += getCity(i)+"|";
}
return geneString;
}
}
最后，让我们创建模拟退火算法：
package sa;
public class SimulatedAnnealing {
// Calculate the acceptance probability
public static double acceptanceProbability(int energy, int newEnergy, double temperature) {
// If the new solution is better, accept it
if (newEnergy < energy) {
return 1.0;
}
// If the new solution is worse, calculate an acceptance probability
return Math.exp((energy - newEnergy) / temperature);
}
public static void main(String[] args) {
// Create and add our cities
City city = new City(60, 200);
TourManager.addCity(city);
City city2 = new City(180, 200);
TourManager.addCity(city2);
City city3 = new City(80, 180);
TourManager.addCity(city3);
City city4 = new City(140, 180);
TourManager.addCity(city4);
City city5 = new City(20, 160);
TourManager.addCity(city5);
City city6 = new City(100, 160);
TourManager.addCity(city6);
City city7 = new City(200, 160);
TourManager.addCity(city7);
City city8 = new City(140, 140);
TourManager.addCity(city8);
City city9 = new City(40, 120);
TourManager.addCity(city9);
City city10 = new City(100, 120);
TourManager.addCity(city10);
City city11 = new City(180, 100);
TourManager.addCity(city11);
City city12 = new City(60, 80);
TourManager.addCity(city12);
City city13 = new City(120, 80);
TourManager.addCity(city13);
City city14 = new City(180, 60);
TourManager.addCity(city14);
City city15 = new City(20, 40);
TourManager.addCity(city15);
City city16 = new City(100, 40);
TourManager.addCity(city16);
City city17 = new City(200, 40);
TourManager.addCity(city17);
City city18 = new City(20, 20);
TourManager.addCity(city18);
City city19 = new City(60, 20);
TourManager.addCity(city19);
City city20 = new City(160, 20);
TourManager.addCity(city20);
// Set initial temp
double temp = 10000;
// Cooling rate
double coolingRate = 0.003;
// Initialize intial solution
Tour currentSolution = new Tour();
currentSolution.generateIndividual();
System.out.println("Initial solution distance: " + currentSolution.getDistance());
// Set as current best
Tour best = new Tour(currentSolution.getTour());
// Loop until system has cooled
while (temp > 1) {
// Create new neighbour tour
Tour newSolution = new Tour(currentSolution.getTour());
// Get a random positions in the tour
int tourPos1 = (int) (newSolution.tourSize() * Math.random());
int tourPos2 = (int) (newSolution.tourSize() * Math.random());
// Get the cities at selected positions in the tour
City citySwap1 = newSolution.getCity(tourPos1);
City citySwap2 = newSolution.getCity(tourPos2);
// Swap them
newSolution.setCity(tourPos2, citySwap1);
newSolution.setCity(tourPos1, citySwap2);
// Get energy of solutions
int currentEnergy = currentSolution.getDistance();
int neighbourEnergy = newSolution.getDistance();
// Decide if we should accept the neighbour
if (acceptanceProbability(currentEnergy, neighbourEnergy, temp) > Math.random()) {
currentSolution = new Tour(newSolution.getTour());
}
// Keep track of the best solution found
if (currentSolution.getDistance() < best.getDistance()) {
best = new Tour(currentSolution.getTour());
}
// Cool system
temp *= 1-coolingRate;
}
System.out.println("Final solution distance: " + best.getDistance());
System.out.println("Tour: " + best);
}
}
结果如下：
Initial solution distance: 1966
Final solution distance: 911
Tour: |180, 200|200, 160|140, 140|180, 100|180, 60|200, 40|160, 20|120, 80|100, 40|60, 20|20, 20|20, 40|60, 80|100, 120|40, 120|20, 160|60, 200|80, 180|100, 160|140, 180|
在这个例子中，我们能够超过我们初始随机生成路径的一半以上。很大程度上说明，当应用到某些类型的优化问题时，这个相对简单的算法是多么方便。
模拟退火算法的参数控制问题
模拟退火算法的应用很广泛，可以求解 NP 完全问题，但其参数难以控制，其主要问题有以下三点：
(1) 温度 T 的初始值设置问题。
温度T的初始值设置是影响模拟退火算法全局搜索性能的重要因素之一、初始温度高，则搜索到全局最优解的可能性大，但因此要花费大量的计算时间；反之，则可节约计算时间，但全局搜索性能可能受到影响。实际应用过程中，初始温度一般需要依据实验结果进行若干次调整。
(2) 退火速度问题。
模拟退火算法的全局搜索性能也与退火速度密切相关。一般来说，同一温度下的“充分”搜索(退火)是相当必要的，但这需要计算时间。实际应用中，要针对具体问题的性质和特征设置合理的退火平衡条件。
(3) 温度管理问题。
温度管理问题也是模拟退火算法难以处理的问题之一。实际应用中，由于必须考虑计算复杂度的切实可行性等问题，常采用如下所示的降温方式：T(t+1)＝k×T(t)
式中 k 为正的略小于 1.00 的常数， t 为降温的次数。
例题推荐
给定 n 个质点，求重心，这 n 个质点的重心满足Σ(重心到点i的距离)*g[i]最小。—BZOJ 3680 参考题解请看这里
给 n 个点，找出一个点，使这个点到其他所有点的距离之和最小，也就是求费马点。—POJ 2420
给定三维空间的 n 点，找出一个半径最小的球把这些点全部包围住。—POJ 2069
平面上给定 n 条线段，找出一个点，使这个点到这 n 条线段的距离和最小。参考源码在这里
地图中有 N 个陷阱，给出他们的坐标，求一个点，使得这个点到所有陷阱的最小距离最大。—POJ 1379
求一个椭球面上的一个点到原点的最短距离。—HDU 5017
找出一个点使得这个店到 n 个点的最长距离最短，即求最小覆盖圆的半径。—HDU 3932
给一个矩阵的长宽，再给 n 个点，求矩阵区域内某个点到各个点的最小距离的最大值，输出所求点的坐标。—HDU 1109
给定 n 个点的一个多边形，一个圆的半径，判断圆是否可以放在多边形里。—HDU 3644
给定n个点的坐标和它 x 和 y 方向的分速度，要求在任意时刻两两点之间距离最大值中的最小值。—HDU 4717
参考文献
6234
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
