模擬退火算法SA原理及python、java、php、c++語言代碼實現TSP旅行商問題，智能優化算法，隨機尋優算法，全局最短路徑

模擬退火算法SA原理及python、java、php、c++語言代碼實現TSP旅行商問題，智能優化算法，隨機尋優算法，全局最短路徑

模擬退火算法(Simulated Annealing，SA)最早的思想是由N. Metropolis等人於1953年提出。1983 年,S. Kirkpatrick 等成功地將退火思想引入到組合優化領域。
來源於固體退火原理，將固體加溫至充分高，再讓其徐徐冷卻，加溫時，固體內部粒子隨溫升變為無序狀，內能增大，而徐徐冷卻時粒子漸趨有序，在每個溫度都達到平衡態，最后在常溫時達到基態，內能減為最小。

它是基於Monte-Carlo(蒙特卡洛)迭代求解策略的一種隨機尋優算法，其出發點是基於物理中固體物質的退火過程與一般組合優化問題之間的相似性。

模擬退火算法從某一較高初溫出發，伴隨溫度參數的不斷下降,結合概率突跳特性在解空間中隨機尋找目標函數的全局最優解，即在局部最優解能概率性地跳出並最終趨於全局最優。模擬退火算法是一種通用的優化算法，理論上算法具有概率的全局優化性能,目前已在工程中得到了廣泛應用，諸如VLSI、生產調度、控制工程、機器學習、神經網絡、信號處理等領域。

模擬退火算法是通過賦予搜索過程一種時變且最終趨於零的概率突跳性，從而可有效避免陷入局部極小並最終趨於全局最優的串行結構的優化算法。

這里的“一定的概率”的計算參考了金屬冶煉的退火過程，這也是模擬退火算法名稱的由來。
　　根據熱力學的原理，在溫度為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)來接受這樣的移動。

二、什么是智能優化算法
啟發式算法，是一種具有全局優化性能、通用性強且適用於並行處理的算法。這種算法一般具有嚴密的理論依據，而不是單純憑專家經驗，理論上可以在一定的時間內找到最優解或近似最優解。

常用的智能優化算法
遺傳算法（Genetic Algorithm, GA)
模擬退火算法（Simulated Annealing, SA)
禁忌搜索算法（Tabu Search, TS)
神經網絡 （Neural Network）
蟻群算法（Ant Colony Optimization，ACO）
爬山算法(局部最優解法) ( Hill Climbing )

三、TSP問題描述

旅行商問題描述的是一個旅行商要在N之間游走來販賣貨物，而路費的成本與距離成正比，因此他想在N個城市之間找到僅一條最短路徑，即能去所有城市一次且僅一次，最后回到出發點。
這個簡單的問題被證明為NP問題，對於N個城市的問題計算量為O(N!),對於N=40時的計算量就已是目前全世界最強的計算機都難以解決。因此必須尋找其它的可行的算法（TSP問題算法目前有：爬山、禁忌、蟻群、退火、遺傳算法），模擬退火算法就是其中一種。

模擬退火算法其特點是在開始搜索階段解的質量提高比較緩慢，但是到了迭代后期，它的解的質量提高明顯，所以如果在求解過程中，對迭代步數限制比較嚴格的話，模擬退火算法在有限的迭代步數內很難得到高質量的解。總體而言模擬退火算法比較適合用於有充足計算資源的問題求解。

根據Metropolis准則，粒子在溫度T時趨於平衡的概率為exp(-ΔE/(kT))，其中E為溫度T時的內能，ΔE為其改變數,k為Boltzmann常數。
若f( Y(i+1) ) <= f( Y(i) ) (即移動后得到更優解)，則總是接受該移動；
若f( Y(i+1) ) > f( Y(i) ) (即移動后的解比當前解要差)，則以一定的概率接受移動，而且這個概率隨着時間推移逐漸降低（逐漸降低才能趨向穩定）相當於上圖中，從B移向BC之間的小波峰時，每次右移（即接受一個更糟糕值）的概率在逐漸降低。如果這個坡特別長，那么很有可能最終我們並不會翻過這個坡。如果它不太長，這很有可能會翻過它，這取決於衰減 t 值的設定。

 

 

====================

模擬退火算法 偽代碼
s:=s0;e:=E(s)//設定目前狀態為s0，其能量E(s0)
k:=0//評估次數k
whilek<kmaxande>emax//若還有時間（評估次數k還不到kmax）且結果還不夠好（能量e不夠低）則：
sn:=neighbour(s)//隨機選取一臨近狀態sn
en:=Esn)//sn的能量為E(sn)
ifrandom()<P(e,en,temp(k/kmax))then//決定是否移至臨近狀態sn
s:=sn;e:=en//移至臨近狀態sn
k:=k+1//評估完成，次數k加一
returns//回轉狀態s
===============

模擬退火算法python有庫可以調用，實現起來很簡單

Python 代碼實現編輯
step1：在GitHub上下載常用的 scikit-opt [2] 庫。
step2：設立目標函數並執行模擬退火算法。

def demo_func(x):
    x1, x2, x3 = x
    return x1 ** 2 + (x2 - 0.05) ** 2 + x3 ** 2
 
from sko.SA import SA

sa = SA(func=demo_func, x0=[1, 1, 1])
x_star, y_star = sa.fit()

step2(2)：如果是最短路徑問題（TSP）
sa_tsp = SA_TSP(func=demo_func, x0=range(num_points))
best_points, best_distance = sa_tsp.fit()

====================

TSP問題的SA算法PHP實現

<?php
class SimulatedAnnealing
{
    private $map;//地圖，按照矩陣的方式
    private $n;//地圖規模
    private $T;//初始溫度T
    private $L;//每個溫度下達到降溫狀態的迭代次數
    private $l;//降溫常數小於卻接近於1的常數。λ通常的取值在0.8至0.99之間。在每一溫度下，實驗足夠多的轉移次數
    private $ord_temp;//降溫終止條件，相當於降溫到常溫
    private $path;//輸出結果
    
    public function __construct($_map,$_n,$_T,$_L)
    {
        $this->map = $_map;
        $this->n = $_n;
        $this->T = $_T;
        $this->L = $_L;
        $this->l = 0.95;
        $this->ord_temp = 1;
    }
    
    function randomFloat($min = 0, $max = 1) 
    {
        return $min + mt_rand() / mt_getrandmax() * ($max - $min);
    }
    
    public function output()
    {
        foreach($this->path as $key => $value)
        {
            echo $value."->";
        }
    }
    
    //第二個元素向右隨機移動1至n-2位
    public function new_S($_S)
    {
        $temp_S = $_S;
        $shift = rand(1,$this->n-2);
        
        $ts = $temp_S[1];
        $temp_S[1] = $temp_S[1+$shift];
        $temp_S[1+$shift] = $ts;
        
        return $temp_S;
    }
    
    public function cost($_S)
    {
        $_cost = 0;
        for($i=0;$i<$this->n-1;$i++)
        {
            $_cost += $this->map[$_S[$i]][$_S[$i+1]];
        }
        $_cost += $this->map[$_S[$i]][0];
        return $_cost;
    }
    
    public function calculate()
    {
        //初始解
        $S = array();
        
        for($i=0;$i<$this->n;$i++)
        {
            $S[$i] = $i;
        }
        $S[] = 0;
        $this->path = $S;
        //進入降溫過程，初始溫度為T
        $t = $this->T;
        while($t > $this->ord_temp)
        {
            for($i=0;$i<$this->L;$i++)
            { 
                //產生新解
                $S1 = $this->new_S($S);
                //判斷新解的價值
                $K = $this->cost($S1) - $this->cost($S);
                if($K < 0)
                {
                    $S = $S1;
                    if($this->cost($S) < $this->cost($this->path))
                        $this->path = $S;
                }
                else
                {
                    //不好的解根據Metropolis准則接受
                    if($this->randomFloat(0,1) < exp(-$K/$t))
                    {
                        $S = $S1;
                    }
                }
            }
            //這里可以按照如果連續幾次降溫能量不變作為終止循序條件
            //本例按照降溫到常溫終止，不計算是否能量不變，需要初始溫度給定足夠好
            $t = $t * $this->l;
        }
        //輸出結果
        $this->output();
        echo '<br>The min cost is '.$this->cost($this->path);
    }
}
?>

調用過程：

$map = array(
        array(0,3,6,7),
        array(5,0,2,3),
        array(6,4,0,2),
        array(3,7,5,0),
    );
    $n = 4;
    $T = 20;
    $L = 100;
    $SA = new SimulatedAnnealing($map,$n,$T,$L);
    $SA->calculate();

結果輸出：

0->1->2->3->0->
The min cost is 10

=========================

TSP問題JAVA實現代碼：

package noah;
 
import java.io.BufferedReader;
import java.io.FileInputStream;
import java.io.IOException;
import java.io.InputStreamReader;
import java.util.Random;
 
public class SA {
 
    private int cityNum; // 城市數量，編碼長度
    private int N;// 每個溫度迭代步長
    private int T;// 降溫次數
    private float a;// 降溫系數
    private float t0;// 初始溫度
 
    private int[][] distance; // 距離矩陣
    private int bestT;// 最佳出現代數
 
    private int[] Ghh;// 初始路徑編碼
    private int GhhEvaluation;
    private int[] bestGh;// 最好的路徑編碼
    private int bestEvaluation;
    private int[] tempGhh;// 存放臨時編碼
    private int tempEvaluation;
 
    private Random random;
 
    public SA() {
 
    }
 
    /**
     * constructor of GA
     * 
     * @param cn
     *            城市數量
     * @param t
     *            降溫次數
     * @param n
     *            每個溫度迭代步長
     * @param tt
     *            初始溫度
     * @param aa
     *            降溫系數
     * 
     **/
    public SA(int cn, int n, int t, float tt, float aa) {
        cityNum = cn;
        N = n;
        T = t;
        t0 = tt;
        a = aa;
    }
 
    // 給編譯器一條指令，告訴它對被批注的代碼元素內部的某些警告保持靜默
    @SuppressWarnings("resource")
    /**
     * 初始化Tabu算法類
     * @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
            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])) / 10.0);
                // 四舍五入，取整
                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
    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];
        }
    }
 
    public int evaluate(int[] chr) {
        // 0123
        int len = 0;
        // 編碼，起始城市,城市1,城市2...城市n
        for (int i = 1; i < cityNum; i++) {
            len += distance[chr[i - 1]][chr[i]];
        }
        // 城市n,起始城市
        len += distance[chr[cityNum - 1]][chr[0]];
        return len;
    }
 
    // 鄰域交換，得到鄰居
    public void Linju(int[] Gh, int[] tempGh) {
        int i, temp;
        int ran1, ran2;
 
        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) {
            n = 0;
            while (n < N) {
                Linju(Ghh, tempGhh);// 得到當前編碼Ghh的鄰域編碼tempGhh
                tempEvaluation = evaluate(tempGhh);
                if (tempEvaluation < bestEvaluation)
                {
                    copyGh(tempGhh, bestGh);
                    bestT = k;
                    bestEvaluation = tempEvaluation;
                }
                r = random.nextFloat();
                if (tempEvaluation < GhhEvaluation
                        || Math.exp((GhhEvaluation - tempEvaluation) / t) > r) {
                    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....");
        SA sa = new SA(48, 40, 400, 250.0f, 0.98f);
        sa.init("c://data.txt");
        sa.solve();
    }
}

==================
C++實現的代碼：

#include <iostream>
#include <string.h>
#include <stdlib.h>
#include <algorithm>
#include <stdio.h>
#include <time.h>
#include <math.h>

#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;
}