模拟退火
模拟退火
爬山算法（Hill Climbing）
介绍模拟退火前，先介绍爬山算法。爬山算法是一种简单的贪心搜索算法，该算法每次从当前的解空中选择一个最优解作为当前解，直到达到一个局部最优解。
爬山算法实现很简单，其主要缺点是会陷入局部最优解，而不一定能搜索到全局最优解。如下图所示：假设C点为当前解，爬山算法搜索到A点这个局部最优解就会停止搜索，因为在A点无论向那个方向小幅度移动都不能得到更优的解。
模拟退火算法
爬山法是完完全全的贪心法，每次都鼠目寸光的选择一个当前最优解，因此只能搜索到局部的最优值。模拟退火其实也是一种贪心算法，但是它的搜索过程引入了随机因素。模拟退火算法以一定的概率来接受一个比当前解要差的解，因此有可能会跳出这个局部的最优解，达到全局的最优解。以上图为列，模拟退火算法在搜索到局部最优解A后，会以一定的概率接受到E的移动。也许经过几次这样的不是局部最优的移动后会到达D点，于是就跳出了局部最大值A。
模拟退火算法描述：
若(即移动后得到更优解)，则总是接受该移动
若(即移动后的解比当前解要差)，则以一定的概率接受移动，而且这个概率随着时间推移逐渐降低（逐渐降低才能去趋向稳定）
这里的"一定概率"的计算参考了金属冶炼的退火过程，这也是模拟退火算法名称的由来。
根据热力学原理，在温度为T时，出现能量差为dE的降温的概率为P(dE)，表示为：
P(dE)=exp(dE/(KT))
其中k是一个常数，exp表示自然指数，且dE<0。这条公式说白了就是：温度越高，出现一次能量差为dE的降温的概率就越大。温度越低，则出现降温的概率就越小。又由于dE总是小于0（否则就不叫退火了），因此dE/kT<0,所以P(dE)的函数取值范围是（0,1）。
随着温度T的降低，P(dE)会逐渐降低。
我们将一次向较差解的移动看作一次温度跳变过程，我们以概率P(dE)来接受这样的移动。
关于爬山算法与模拟退火，有一个有趣的比喻：
爬山算法：兔子朝着比现在高的地方跳去。它找到了不远处的最高峰。但是这座山不一定是珠穆朗玛峰。这就是爬山算法，它不能保证局部最优值就是全局最优值。
模拟退火：兔子喝醉了。它随机地跳了很长时间。这期间，它可能走向高处，也可能踏入平地。但是，它渐渐清醒了并朝最高方向跳去。这就是模拟退火。
模拟退火算法的关键步骤
模拟退火算法的基本思想是：由初始解和控制参数初值开始，对当前解重复"产生新解->计算目标函数差->接受或舍弃"的迭代，并逐步衰减控制参数值，算法终止时，当前解的值即为近似最优解，这便是基于蒙特卡罗迭代求解法的一种启发式随机搜索过程。
用模拟退火解决TSP问题
import java.io.BufferedReader; import java.io.FileInputStream; import java.io.IOException; import java.io.InputStreamReader; import java.util.Random; /** * an example of Simulated Annealing Agorithm to solve TSP * 算法主要思路： * 1、邻域变换函数：随机交换当前的两个元素。 * 2、评估函数：以当前方案的路程总和作为评估值。 * 数据集来源：tsplib 上的数据eil51.tsp去除头部。 */ public class SimulatedAnnealingTSP { private int cityNum; // 城市数量，编码长度 private int N;// 每个温度迭代步长 private int T;// 降温次数 private float a;// 降温系数 private float t0;// 初始温度， 尽可能大，最后结果与初始温度无关 private int[][] distance; // 距离矩阵 private int bestT;// 最佳出现代数 private int[] Ghh;// 初始路径编码，当前编码（相对于temGhh 是上次编码） private int GhhEvaluation;//评估值 private int[] bestGh;// 最好的路径编码 private int bestEvaluation;//最好的路径编码的评估值 private int[] tempGhh;// 存放临时编码 private int tempEvaluation;//评估值 private Random random; /** * constructor * * @param cn * 城市数量 * @param t * 降温次数 * @param n * 每个温度迭代步长 * @param tt * 初始温度 * @param aa * 降温系数 * **/ public SimulatedAnnealingTSP(int cn, int n, int t, float tt, float aa) { cityNum = cn; N = n; T = t; t0 = tt; a = aa; } @SuppressWarnings("resource")// 给编译器一条指令，告诉它对被批注的代码元素内部的某些警告保持静默 /** * 读取数据集，初始化cityNum, distance[][], * @param filename 数据文件名，该文件存储所有城市节点坐标数据 * @throws IOException */ private void init(String filename) throws IOException { // 读取数据 int[] x; int[] y; String strbuff; BufferedReader data = new BufferedReader(new InputStreamReader( new FileInputStream(filename))); distance = new int[cityNum][cityNum]; x = new int[cityNum]; y = new int[cityNum]; for (int i = 0; i < cityNum; i++) { // 读取一行数据，数据格式1 6734 1453 (编号 x坐标 y坐标) strbuff = data.readLine(); // 字符分割 String[] strcol = strbuff.split(" "); x[i] = Integer.valueOf(strcol[1]);// x坐标 y[i] = Integer.valueOf(strcol[2]);// y坐标 } // 计算距离矩阵 // 针对具体问题，距离计算方法也不一样，此处用的是att48作为案例，它有48个城市，距离计算方法为伪欧氏距离，最优值为10628 for (int i = 0; i < cityNum - 1; i++) { distance[i][i] = 0; // 对角线为0 for (int j = i + 1; j < cityNum; j++) { double rij = Math .sqrt(((x[i] - x[j]) * (x[i] - x[j]) + (y[i] - y[j]) * (y[i] - y[j]))); //小数部分进位， 1.1-> 2 1.7->2 int tij = (int) Math.round(rij);// 四舍五入 if (tij < rij) { distance[i][j] = tij + 1; distance[j][i] = distance[i][j]; } else { distance[i][j] = tij; distance[j][i] = distance[i][j]; } } } distance[cityNum - 1][cityNum - 1] = 0; Ghh = new int[cityNum]; bestGh = new int[cityNum]; bestEvaluation = Integer.MAX_VALUE;//设置初始值位最大 tempGhh = new int[cityNum]; tempEvaluation = Integer.MAX_VALUE;//设置初始值位最大 bestT = 0;//迭代歩数 random = new Random(System.currentTimeMillis()); System.out.println("城市数："+ cityNum+","+ "每个温度迭代步长（内层迭代上限）："+N+","+ "降温次数（外层迭代上限）："+T+","+ "降温系数："+a+","+ "初始温度："+t0); } /** * 初始化编码Ghh * 随机产生一组城市号 0到cityNum的序列 * 可以用任意方式指定 */ void initGroup() { int i, j; Ghh[0] = random.nextInt(65535) % cityNum; for (i = 1; i < cityNum;)// 编码长度 { Ghh[i] = random.nextInt(65535) % cityNum; for (j = 0; j < i; j++) { if (Ghh[i] == Ghh[j]) { break; } } if (j == i) { i++; } } } // 复制编码体，复制编码Gha到Ghb public void copyGh(int[] Gha, int[] Ghb) { for (int i = 0; i < cityNum; i++) { Ghb[i] = Gha[i]; } } /** * 估计函数 * @param chr 路径编码，起始城市,城市1,城市2...城市n * @return 路径总长度 */ public int evaluate(int[] chr) { // 0123 int len = 0; // 计算路径总长度 for (int i = 1; i < cityNum; i++) { len += distance[chr[i - 1]][chr[i]]; } // 加上从最后一个城市回到出发城市的路程 len += distance[chr[cityNum - 1]][chr[0]]; return len; } /** * 邻域交换，得到当前编码Ghh的邻域编码tempGhh * 随机交换Ghh的两个编码位 */ public void Linju(int[] Gh, int[] tempGh) { int i, temp; int ran1, ran2; //copy for (i = 0; i < cityNum; i++) { tempGh[i] = Gh[i]; } //随机生成两个不同的编码位下标 ran1 = random.nextInt(65535) % cityNum; ran2 = random.nextInt(65535) % cityNum; while (ran1 == ran2) { ran2 = random.nextInt(65535) % cityNum; } //交换两个编码位 temp = tempGh[ran1]; tempGh[ran1] = tempGh[ran2]; tempGh[ran2] = temp; } public void solve() { // 初始化编码Ghh initGroup(); copyGh(Ghh, bestGh);// 复制当前编码Ghh到最好编码bestGh bestEvaluation = evaluate(Ghh); GhhEvaluation = bestEvaluation; int k = 0;// 降温次数 int n = 0;// 迭代步数 float t = t0; float r = 0.0f; while (k < T) {//T位降温次数上限 n = 0; while (n < N) { Linju(Ghh, tempGhh);// 得到当前编码Ghh的邻域编码tempGhh tempEvaluation = evaluate(tempGhh);//评估当前路劲，结果为总的路程 if (tempEvaluation < bestEvaluation)//if the temp one is better { copyGh(tempGhh, bestGh);//copy tempGhh to bestGh bestT = k; bestEvaluation = tempEvaluation; } //根据Metropolis准则判断是否接受当前解 r = random.nextFloat();// 返回值属于[0,1] if (tempEvaluation < GhhEvaluation || Math.exp((GhhEvaluation - tempEvaluation) / t) > r) {// t = current temperature copyGh(tempGhh, Ghh); GhhEvaluation = tempEvaluation; } n++; } t = a * t; k++; } System.out.println("此次运行的最佳长度出现代数："); System.out.println(bestT); System.out.println("此次运行的最佳长度:"); System.out.println(bestEvaluation); System.out.println("此次运行的最佳路径："); for (int i = 0; i < cityNum; i++) { System.out.print(bestGh[i] + ","); if (i % 10 == 0 && i != 0) { System.out.println(); } } } /** * @param args * @throws IOException */ public static void main(String[] args) throws IOException { System.out.println("Start...."); SimulatedAnnealingTSP sa = new SimulatedAnnealingTSP(51, 8000, 1500, 10000.0f, 0.902f); sa.init("eil51.tsp"); sa.solve(); } }
bestEvaluation
bestEvaluation
