前言
较小规模的0-1规划问题可以用Lingo求解(可以参考《2011年国赛B题论文求解复现》),在数据量较大时,Lingo会因为复杂度过高无法求解。这个时候需要使用模拟退火的算法。
下面用来实践模拟退火的例子取自2011年国赛B题的某篇获奖论文。关于这篇论文的分析,可以参考《2011年国赛B题论文研读》。
一、模拟退火算法的原理【1】
模拟退火求解0-1规划有几个要素:
一是解空间,二是目标函数,三是接受准则。
假设你在一个山谷里,要寻找最低点,这里海拔就是目标函数(
F
F
F)。一个可行的策略是,在你的周围用脚试探,如果发现比现在所处的点低的点,你就移动过去(我们把这种策略叫“趋优”),直到在你周围找不到比现在所处的点更低的点。这种方法的弊端是,你可能会落入局部最优出不来。
模拟退火算法就是为了避免落入局部最优的。思想就是,如果你试探的点比你现在所处的点还高,也仍允许你有一定机会移动到那个地方,也就是为了能试探更多的点,放弃了“趋优”的原则。什么时候可以移动,这由接受准则(
P
P
P)决定,接受准则简单理解就是可移动的概率。
模拟退火算法里还有几个概念:
温度、温度衰减系数、马尔可夫链长。
在你寻找最优解的开始,你一定是希望多做尝试,所以即便你试探的点比你现在所处的点高很多,你也应该宽容。而到尝试后期,你对不满足“趋优”原则的尝试应该减少宽容,因为你已经趋近最优解了。
概括一下这个规律:开始活跃,后来稳定。
这个规律和固体退火过程类似。“将固体加温至充分高,再让其徐徐冷却,加温时,固体内部粒子随温升变为无序状,内能增大,而徐徐冷却时粒子渐趋有序,在每个温度都达到平衡态,最后在常温时达到基态,内能减为最小【2】”。所以这个算法被成为模拟退火算法。
而温度(T)就是控制这个活越程度的因素。我们来看接受准则的表达:
p
=
0
,
i
f
尝
试
的
点
不
满
足
约
束
条
件
p=0,if尝试的点不满足约束条件
p=0,if尝试的点不满足约束条件
p
=
1
,
i
f
尝
试
的
点
不
满
足
约
束
条
件
p=1,if尝试的点不满足约束条件
p=1,if尝试的点不满足约束条件
p
=
e
x
p
(
−
Δ
F
/
T
)
,
e
l
s
e
p=exp(-ΔF/T),else
p=exp(−ΔF/T),else
Δ
F
ΔF
ΔF是F的变化量。从这里就能看出温度越低,else条件下可移动的概率
p
p
p越低。
在同一个温度下,要尝试迭代若干轮,这个轮数,就是马尔可夫链长(
L
k
L_k
Lk)。而在每迭代
L
k
L_k
Lk后,温度T就衰减到原来的α倍(
0
<
α
<
1
0<α<1
0<α<1),α为温度衰减系数。
初始时的温度记为
T
0
T_0
T0,终止温度记为
T
f
T_f
Tf,当
T
<
T
f
T<T_f
T<Tf时,算法结束。
二、代码(结合问题背景)
1、预处理
最终得到表示无向图的二维数组time和表示发案率的列表w。
import xlrd
import numpy as np
INFINITY = 99999.0
# 打开文件
data = xlrd.open_workbook('cumcm2011B附件2_全市六区交通网路和平台设置的数据表.xls')
#print("sheets:" + str(data.sheet_names())) #打印sheets名
table1 = data.sheet_by_name('全市交通路口节点数据')
loc = [] #节点位置列表
w = [] #发案率列表
# 为使索引值和节点编号对应,填入一个占位数据
loc.append([0, 0]) # 为使索引值和节点编号对应,填入一个占位数据
w.append(0.0) # 同上
#读入A区各节点的坐标数据
for k in range(1, 93):
x = table1.cell(k, 1).value #注意,excel里的计数是从0开始的,和图像化界面里显示的不同
y = table1.cell(k, 2).value
loc.append([x, y])
w.append(table1.cell(k, 4).value)
table2 = data.sheet_by_name('全市交通路口的路线')
distance = np.zeros((93, 93), dtype=np.float)
# 初始化一个二维数组,用来记录无向图
for i in range(1, 93):
for j in range(1, 93):
if i == j:
distance[i, j] = 0.0
else:
distance[i, j] = INFINITY
# 得到距离矩阵
for k in range(1, table2.nrows):
i = int(table2.cell(k, 0).value)
j = int(table2.cell(k, 1).value)
# print(i, j)
if i <= 92 and j <= 92:
distance[i, j] = np.sqrt((loc[i][0]-loc[j][0])**2 + (loc[i][1]-loc[j][1])**2)
distance[j, i] = np.sqrt((loc[i][0]-loc[j][0])**2 + (loc[i][1]-loc[j][1])**2)
# floyd算法求最短距离
for k in range(1, 93):
for i in range(1, 93):
for j in range(1, 93):
if distance[i, k] + distance[k, j] < distance[i, j]:
distance[i, j] = distance[i, k] + distance[k, j]
# 考虑车速60km/h,把距离矩阵变成时间矩阵
time = np.zeros((93, 93), dtype=np.float)
for i in range(1, 93):
for j in range(1, 93):
time[i, j] = distance[i, j]/10
2、规划模型表示
原模型:
目标函数计算函数:
def F (t, w, x):
numerator = 0.0
denominator = 0.0
for i in range(1, 21):
for j in range(1, 93):
numerator += w[j]*t[i, j]*x[i, j]
pass
for j in range(1, 93):
denominator += w[j]
return numerator/denominator
约束1判断函数:
def F1 (x):
judge = True
for j in range(1, 93):
f = 0.0
for i in range(1, 21):
f += x[i, j]
if f < 1:
judge = False
break
return judge
约束2判断函数,注意这里为了不陷入局部最优或者无解情况,弱化了约束条件:
def F2(t, x):
judge = True
flag = 0
a = 20 # 发现原约束太强,这里用来弱化约束, 这样做也有利于找到最优解
for i in range(1, 21):
for j in range(1, 93):
if t[i, j]*x[i, j] > 3+a:
judge = False
flag = 1
break
if flag == 1:
break
return judge
3、模拟退火
# 模拟退火参数设定
T0 = 100000
Tf = 0.0001
alpha = 0.995 #温度衰减系数
Lk = 100 #马尔可夫链长
# 初始化
T = T0 # 初始化温度
x = np.zeros((21, 93), dtype=np.int) # 初始化当前可行解
f = INFINITY # 初始化当前可行解目标函数值
x_best = np.zeros((21, 93), dtype=np.int) # 初始化最优解
f_best = INFINITY #初始化 最优目标函数值
# 设置最初解,保证每个节点都有对应的平台
for j in range(1, 93):
id = np.random.randint(1, 21)
x[id, j] = 1
print("开始迭代:")
iter = 0
while T >= Tf:
print("iter = "+str(iter)+", T = "+str(T)+", f = "+str(f)+", f_best = "+str(f_best)) # 打印迭代信息,用于监控退火过程
iter += 1
for k in range(Lk): # 在每个温度下,搜索Lk次
# 尝试更新解
x_test = x.copy()
for j in np.random.randint(1, 93, 46): # 随机抽取一半平台做更新
for i in range(1, 21):
x_test[i, j] = 0
id = np.random.randint(1, 21)
x_test[id, j] = 1
# 判断是否更新解
P = -1
p = 0.0 # 这两个都用于表示概率
if not (F1(x_test) and F2(time, x_test)):
P = 0
elif F1(x_test) and F2(time, x_test) and F(time, w, x_test) < f:
P = 1
else:
p = np.exp((F(time, w, x_test) - f)/T)
if P == 1 or np.random.rand(1) < p: # 更新解
x = x_test
f = F(time, w, x_test)
if f < f_best:
x_best = x.copy()
f_best = f
T = alpha*T # 降温
print("迭代结束")
print("f_best = "+str(f_best))
4、结果打印
for i in range(1, 21):
print(str(i)+"号平台管辖:")
for j in range(1, 93):
if x_best[i, j] == 1:
print(str(j), end=" ")
print("\n")
三、结果
最终计算出来的结果使目标函数达到7.14585741526687,与原论文的1.013相差仍比较大。后续仍有待进一步调试。欢迎大家给出建议。
参考资料
【1】梁国宏,张生,黄辉,等.一种改进的模拟退火算法求解0-1背包问题[J].广西民族大学学报(自然科学版),2007,13(3):91-93.
【2】百度百科https://baike.baidu.com/item/模拟退火算法/355508?fromtitle=%E6%A8%A1%E6%8B%9F%E9%80%80%E7%81%AB&fromid=8664695&fr=aladdin