模拟退火

参考: http://www.cnblogs.com/heaad/archive/2010/12/20/1911614.html

爬山算法

只能接受下一个解比当前解好的情况

模拟退火

以一定概率接受比当前解差的解，此概率随时间经过逐渐变小

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
/*
* J(y)：在状态y时的评价函数值
* Y(i)：表示当前状态
* Y(i+1)：表示新的状态
* r： 用于控制降温的快慢
* T： 系统的温度，系统初始应该要处于一个高温的状态
* T_min ：温度的下限，若温度T达到T_min，则停止搜索
*/
while( T > T_min )
{
　　dE = J( Y(i+1) ) - J( Y(i) ) ; 

　　if ( dE >=0 ) //表达移动后得到更优解，则总是接受移动
 Y(i+1) = Y(i) ; //接受从Y(i)到Y(i+1)的移动
　　else
　　{
      // 函数exp( dE/T )的取值范围是(0,1) ，dE/T越大，则exp( dE/T )也
      if ( exp( dE/T ) > random( 0 , 1 ) )
      Y(i+1) = Y(i) ; //接受从Y(i)到Y(i+1)的移动
　　}
　　T = r * T ; //降温退火 ，0<r<1 。r越大，降温越慢；r越小，降温越快
　　/*
　　* 若r过大，则搜索到全局最优解的可能会较高，但搜索的过程也就较长。若r过小，则搜索的过程会很快，但最终可能会达到一个局部最优值
　　*/
　　i ++ ;
}

参考：https://zh.wikipedia.org/wiki/%E6%A8%A1%E6%8B%9F%E9%80%80%E7%81%AB

演算步骤

初始化

生成一个可行的解作为当前解输入迭代过程，并定义一个足够大的数值作为初始温度。

迭代过程

迭代过程是模拟退火算法的核心步骤，分为新解的产生和接受新解两部分：

  1. 由一个产生函数从当前解产生一个位于解空间的新解；为便于后续的计算和接受，减少算法耗时，通常选择由当前新解经过简单地变换即可产生新解的方法，如对构成新解的全部或部分元素进行置换、互换等，注意到产生新解的变换方法决定了当前新解的邻域结构，因而对冷却进度表的选取有一定的影响。
  2. 计算与新解所对应的目标函数差。因为目标函数差仅由变换部分产生，所以目标函数差的计算最好按增量计算。事实表明，对大多数应用而言，这是计算目标函数差的最快方法。
  3. 判断新解是否被接受，判断的依据是一个接受准则，最常用的接受准则是Metropolis准则：若Δt′<0则接受S′作为新的当前解S，否则以概率exp（-Δt′/T）接受S′作为新的当前解S。
  4. 当新解被确定接受时，用新解代替当前解，这只需将当前解中对应于产生新解时的变换部分予以实现，同时修正目标函数值即可。此时，当前解实现了一次迭代。可在此基础上开始下一轮试验。而当新解被判定为舍弃时，则在原当前解的基础上继续下一轮试验。

模拟退火算法与初始值无关，算法求得的解与初始解状态S（是算法迭代的起点）无关；模拟退火算法具有渐近收敛性，已在理论上被证明是一种以概率1收敛于全局最优解的全局优化算法；模拟退火算法具有并行性。

停止准则

迭代过程的停止准则：温度T降至某最低值时，完成给定数量迭代中无法接受新解，停止迭代，接受当前寻找的最优解为最终解。

退火方案

在某个温度状态T下，当一定数量的迭代操作完成后，降低温度T，在新的温度状态下执行下一个批次的迭代操作。

TSP问题示例代码

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
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
import numpy as np
import math
import random
from copy import deepcopy
import matplotlib.pyplot as plt


class Tsp(object):
    def __init__(self, city_num, times, steps, init_temperature, simulated_k):
        self.city_num = city_num
        self.distance = np.mat(np.zeros((city_num, city_num)), dtype=int)
        self.init_temperature = init_temperature
        self.times = times
        self.steps = steps
        self.simulated_k = simulated_k

        self.x = list(range(self.city_num))
        self.y = list(range(self.city_num))
        self.now_path = list(range(self.city_num))
        self.new_path = list(range(self.city_num))
        self.best_path = list(range(self.city_num))
        self.now_value = 0
        self.new_value = 0
        self.best_value = -1
        self.path = [int(item)-1 for item in "1 8 38 31 44 18 7 28 6 37 19 27 17 43 30 36 46 33 20 47 21 32 39 48 5  42 24 10 45 35 4 26 2 29 34 41 16 22 3 23 14 25 13 11 12 15 40 9".split()]

        self.best_time = -1

    def read(self, path):
        print("Reading data...")
        # x = [0 for i in range(self.city_num)]
        # y = [0 for i in range(self.city_num)]
        with open(path, "r") as file:
            lines = file.readlines()
            self.x = [int(line.split()[1]) for line in lines]
            self.y = [int(line.split()[2]) for line in lines]
        x = self.x
        y = self.y
        for i in range(0, self.city_num):
            for j in range(0, self.city_num):
                self.distance[i, j] = round(math.sqrt((x[i]-x[j])*(x[i]-x[j]) + (y[i]-y[j])*(y[i]-y[j])))
        print(self.distance)
        print("Read data done.")

    def init_path(self):
        random.shuffle(self.now_path)

    def value_of_path(self, path):
        value = 0
        for (index, i) in enumerate(path):
            value += self.distance[path[index-1], i]
        return value

    def get_neighbor(self):
        self.new_path = deepcopy(self.now_path)
        i_1, i_2 = random.randint(0, self.city_num-1), random.randint(0, self.city_num-1)
        self.new_path[i_1], self.new_path[i_2] = self.new_path[i_2], self.new_path[i_1]

    def solve(self):
        self.init_path()
        self.now_value = self.value_of_path(self.now_path)
        print(self.now_value)
        temperature = self.init_temperature
        k = 0
        # while k < self.times:
        for k in range(self.times):
            print("k: ", k, "  now best: ", self.best_value)
            for n in range(self.steps):
                self.get_neighbor()
                self.new_value = self.value_of_path(self.new_path)
                if self.new_value < self.best_value or self.best_value == -1:
                    self.best_value = self.new_value
                    self.best_path = deepcopy(self.new_path)
                    self.best_time = k
                    print(self.best_value)
                # random_value = random.random()
                if self.new_value < self.now_value or math.exp((self.now_value - self.new_value)/temperature)>random.random():
                    self.now_path = deepcopy(self.new_path)
                    self.now_value = self.new_value
            temperature *= self.simulated_k

        print("best: ", self.best_value)
        # self.draw_best()
        # self.draw()
        # self.show()

    def test(self):
        # best = 33551(ceil)  33522(round)
        print(self.value_of_path(self.path))

    def draw_path(self, path):
        x = [self.x[i] for i in path]
        y = [self.y[i] for i in path]
        x.append(self.x[path[0]])
        y.append(self.y[path[0]])
        plt.plot(x, y, "-o")
        # plt.show()

    def draw_best(self):
        self.draw_path(self.path)

    def draw(self):
        self.draw_path(self.best_path)

    def show(self):
        plt.show()


if __name__ == "__main__":
    tsp = Tsp(48, 1000, 1000, 10000, 0.992)
    tsp.read("att48.tsp")
    tsp.solve()
    tsp.draw_best()
    # tsp.init_path()
    # tsp.draw_path(tsp.now_path)
    # tsp.test()
    tsp.draw()
    tsp.show()

运行结果

已知结果取round的情况下最优解是33522。

我以代码中的参数迭代1000次的结果是34384。

可视化如下：

图中蓝色路线为最优解，橙色路线为SA跑出来的解。