模拟退火算法求解大规模0-1规划

前言

较小规模的0-1规划问题可以用Lingo求解（可以参考《2011年国赛B题论文求解复现》），在数据量较大时，Lingo会因为复杂度过高无法求解。这个时候需要使用模拟退火的算法。
下面用来实践模拟退火的例子取自2011年国赛B题的某篇获奖论文。关于这篇论文的分析，可以参考《2011年国赛B题论文研读》。

一、模拟退火算法的原理【1】

模拟退火求解0-1规划有几个要素：
一是解空间，二是目标函数，三是接受准则。
假设你在一个山谷里，要寻找最低点，这里海拔就是目标函数($F$)。一个可行的策略是，在你的周围用脚试探，如果发现比现在所处的点低的点，你就移动过去（我们把这种策略叫“趋优”），直到在你周围找不到比现在所处的点更低的点。这种方法的弊端是，你可能会落入局部最优出不来。
模拟退火算法就是为了避免落入局部最优的。思想就是，如果你试探的点比你现在所处的点还高，也仍允许你有一定机会移动到那个地方，也就是为了能试探更多的点，放弃了“趋优”的原则。什么时候可以移动，这由接受准则($P$)决定，接受准则简单理解就是可移动的概率。
模拟退火算法里还有几个概念：
温度、温度衰减系数、马尔可夫链长。
在你寻找最优解的开始，你一定是希望多做尝试，所以即便你试探的点比你现在所处的点高很多，你也应该宽容。而到尝试后期，你对不满足“趋优”原则的尝试应该减少宽容，因为你已经趋近最优解了。
概括一下这个规律：开始活跃，后来稳定。
这个规律和固体退火过程类似。“将固体加温至充分高，再让其徐徐冷却，加温时，固体内部粒子随温升变为无序状，内能增大，而徐徐冷却时粒子渐趋有序，在每个温度都达到平衡态，最后在常温时达到基态，内能减为最小【2】”。所以这个算法被成为模拟退火算法。
而温度(T)就是控制这个活越程度的因素。我们来看接受准则的表达：
$$p=0,if尝试的点不满足约束条件$$
$$p=1,if尝试的点不满足约束条件$$
$$p=exp(-\Delta F/T),else$$
$\Delta F$是F的变化量。从这里就能看出温度越低，else条件下可移动的概率$p$越低。
在同一个温度下，要尝试迭代若干轮，这个轮数，就是马尔可夫链长($L_k$)。而在每迭代$L_k$后，温度T就衰减到原来的α倍（$0<\alpha <1$），α为温度衰减系数。
初始时的温度记为$T_0$，终止温度记为$T_f$，当$T<T_f$时，算法结束。

二、代码（结合问题背景）

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