🌞

最优化:Python3 模拟退火算法

背景

这个TP的背景是这样的:有一个商务旅行团需要去n个城市,怎样安排他们的行程可以使行程变得更短?行程的起点和终点都在巴黎。为了解决这个问题,可以使用两种算法:第一种就是最简单的求解n个城市n!种行程的各个距离,然后选取最短的行程;第二种使用la méthode du recuit simulé 模拟退火算法

关于模拟退火算法

模拟退火算法的思想借鉴于固体的退火原理,当固体的温度很高的时候,内能比较大,内部粒子处于快速无序运动。当温度慢慢降低的过程中,固体的内能减小,慢慢趋于有序。最终,当固体处于常温时,内能达到最小。此时,粒子最为稳定。

模拟退火算法从某一较高的温度出发,这个温度称为初始温度,伴随着温度参数的不断下降,算法中的解趋于稳定。但是,可能这样的稳定解是一个局部最优解。此时,模拟退火算法中会以一定的概率跳出这样的局部最优解,以寻找目标函数的全局最优解。

导入所需的模块

1
2
3
4
import numpy as np
import time
import itertools
import matplotlib.pyplot as plt

单个路程的总距离计算

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
# 两个输入

# city_list : 需要去的城市列表

# city_dict : 城市的坐标字典

# 一个输出

# 路程的总距离

    distance = 0 # 初始化距离

    distance+=np.sqrt(city_dict[city_list[0]]['x']**2+city_dict[city_list[0]]['y']**2)  # 就是求两个坐标点之间的距离,先求第一个(巴黎的坐标为(0,0))

    # 循环求并累加

    for i in range(len(city_list)-1): 
        x1=city_dict[city_list[i]]['x']
        y1=city_dict[city_list[i]]['y']
        x2=city_dict[city_list[i+1]]['x']
        y2=city_dict[city_list[i+1]]['y']
        distance=distance+np.sqrt((x1-x2)**2+(y1-y2)**2)
    # 求最后一个(巴黎的坐标为(0,0))

    distance+=np.sqrt(city_dict[city_list[len(city_list)-1]]['x']**2+city_dict[city_list[len(city_list)-1]]['y']**2)
    return distance  # 返回总距离

随机打乱原有路程

1
2
3
4
5
# 一个输入:list_in 当前旅行团的行程列表

# 一个输出:尝试路程列表

# np.random.permutation:输入一个数或者数组,生成一个随机序列,对多维数组来说是多维随机打乱

# np.asarray:从已有的数组创建数组

    return np.random.permutation(np.asarray(list_in).tolist())

计算退火新的温度(降温)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
# 输入:h>0 计算参数,浮点数。h越小,算法越有风险;越大算法收敛时间越长。

# 输入:k 算法参数,整数。

# 输入:ind 算法的迭代次数,整数。

# 输入:Temp 当前算法的温度

# 输出:新的参数k

# 输出:新的温度Temp

    while ind <= np.exp((k-1)*h) or ind>np.exp(k*h):
         k+=1
         Temp=1.0/k

    return k,Temp

绘图函数

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
    chemin=['Paris']+chemin+['Paris']
    plt.plot([dico[ville]['x'] for ville in chemin],[dico[ville]['y'] for ville in chemin],'bo-') # 其他城市用蓝色表示

    plt.plot(dico['Paris']['x'],dico['Paris']['y'],'ro') # 巴黎用红色表示

    
    plt.xlim([-500.0,500.0]) # 坐标范围

    plt.ylim([-800.0,300.0])
    plt.title("Trajet")
    ax=plt.gca()
    ax.set_aspect('equal')

def plot_temp(Temp_list):
    plt.figure()
    plt.plot(Temp_list)
    plt.xlabel('$n$')
    plt.ylabel('$T$')
    plt.title(u'Profil de température')
    plt.grid()

基本参数

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
# 城市列表,字典

dico=dict()
dico['Lille']      ={'x':52.0, 'y':197.0}
dico['Orléans']    ={'x':-33.0, 'y':-105.0}
dico['Lyon']       ={'x':185.0, 'y':-343.0}
dico['Paris']      ={'x':0.0, 'y':0.0}
dico['Marseille']  ={'x':225.0, 'y':-617.0}
dico['Strasbourg'] ={'x':403.0, 'y':-31.0}
dico['Rennes']     ={'x':-300.0, 'y':-88.0}
dico['Metz']       ={'x':285.0, 'y':30.0}
dico['Bordeaux']   ={'x':-213.0, 'y':-448.0}
dico['Perpignan']  ={'x':40.0, 'y':-688.0}
dico['Cherbourg']  ={'x':-289.0, 'y':86.0}

parcours=np.random.permutation(['Marseille',
                                'Lyon',
                                'Rennes',
                                'Lille',
                                'Perpignan',
                                'Cherbourg']).tolist()
# 城市越多,传统的算法将花费更多时间。</pre>

用传统解决办法

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
# 计算所有的可能路线,如果是5个城市就是5!种情况;如果是10个城市就是10!种情况。

t1=time.time()   #计时开始


permutations=list()
for i in itertools.permutations(parcours):  #排列顺序可变,元素不重复:每一项用长度为r的元组表示

    permutations.append(i)      #所有的可能情况,一种种情况往上叠。

print('Nombre de trajets étudiés : ', len(permutations))

cost=list()  
for i in range(0,len(permutations)-1):
    parcours_chaque=list(permutations[i])
    parcours_chaque=['Paris']+parcours_chaque+['Paris']
    print(parcours_chaque)
    print(cost_function(parcours_chaque,dico)) #传到上方的cost_function算法里

    cost.append(cost_function(parcours_chaque,dico)) #所有的可能情况,一种种情况往上叠。

t2=time.time() # 计时结束


# 打印结果

cost=np.asarray(cost)
print('Trajet le plus court :')
for ind, item in zip(range(1,len(parcours)+1),permutations[cost.argmin()]):
    print(ind,' ',item)
print('Distance totale : ', cost.min())
print('Temps de calcul : ',t2-t1)

用模拟退火法

 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
candidate=parcours
# 先设置参数

itermax=10000
hpar=3.0
kpar=1
Temp=1.0/kpar
Temp_list=[Temp]

t1=time.time()   #计时开始

for ind in range(itermax):
    #计算新的行程尝试(就是根据现有行程随机生成一个新的)

    new_candidate=compute_new_candidate(candidate)
    #计算新的行程和旧的行程的路程差距

    delta=cost_function(new_candidate,dico)-cost_function(candidate,dico)
    #如果新的候选路线更长,则仍可以一定的概率接受它。

    if delta &lt;= 0:
        candidate = new_candidate
    else:
        u=np.random.uniform(0,1) # 产生一个0到1之间的随机数

        if np.exp(-delta/Temp) >= u :
            candidate = new_candidate           
    # 降温

    kpar, Temp = compute_Temp(hpar,kpar,ind+2,Temp)
    Temp_list.append(Temp)
t2=time.time()

# 打印结果

print('Trajet le plus court :')
for ind, item in zip(range(1,len(parcours)+1),candidate):
    print(ind,' ',item)
print('Distance totale : ', cost_function(candidate,dico))
print('Temps de calcul : ',t2-t1)

# 绘制图表

plot_chemin(parcours,dico)
plot_chemin(candidate,dico)

plot_temp(Temp_list)
plt.show()

参考资料

https://blog.csdn.net/google19890102/article/details/45395257

updatedupdated2020-01-272020-01-27