模拟退火

模拟退火

爬山算法（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