手术的最优化分配（2）——打破对称性的约束(Symmetry-Breaking Constraints)

​上篇文章中的代码和测试用数据出现了问题，需要做一些更正：

（1）测试用的手术时长太长了，如果每个手术室都要加班那么怎么排班总消耗都是一样的，无法验证解的最优性。故将手术时长数据d改用为[142, 276, 9, 211, 117, 223, 244, 333, 352, 94]。另外，手术时长的选取最好满足以下两种情况之一：①所有手术室开放，部分手术室加班；②部分手术室开放，个别手术室少量加班。

（2）70-72行代码有误，应把“y[j::j_]”改为“y[j * j_:]”

这篇文章主要讲解用于提高求解速度的打破对称性约束(Symmetry-Breaking Constraints)，首先我们来了解一下概念的起源，以下划线内容参考“https://swarma.org/?p=28890”。

Symmetry Breaking(对称性破缺)是一个物理学的概念，在物理学中，一个作用于系统的小扰动使系统跨过临界点，通过决定去向分叉的哪个分支来决定系统的命运，这种现象叫做对称性破缺。

图一：一个小球位于中央山丘的山峰处（C）。这是一种不稳定平衡：一个很小的扰动会使它落到左边（L）或右边（R）稳定点。尽管山丘是对称的，没有理由让球落在哪一侧，但观察到的最终状态仍然是不对称的，它总会落到某一侧。

而在优化问题中也常常出现对称性，在上篇文章中求解的手术规划问题中，四个手术室之间的安排可以两两互换，也就是说有4!=4x3x2x1=24个目标函数取值相同的最优解，而最终选用哪一种方案对我们并无区别。

手术室安排问题中有一个解如下：

	手术室0	手术室1	手术室2	手术室3
手术0				1
手术1	1
手术2		1
手术3	1
手术4		1
手术5			1
手术6				1
手术7			1
手术8		1
手术9				1
加班	7		76

因为每个手术室的固定开放时间都相同，我们可以将任意两个手术室计划互换而不改变目标函数的值。这种对称性大大增加了解空间的范围，在求解过程中我们浪费了许多时间来搜索对称性解。在上表安排中，我们自然而然想到了将解固定为下表中的形式

	手术室0	手术室1	手术室2	手术室3
手术0	1
手术1		1
手术2			1
手术3		1
手术4			1
手术5				1
手术6	1
手术7				1
手术8			1
手术9	1
加班		7		76

我们想要将前几台手术在不改变成本的情况下安排在前几个手术室中，即在表中我们想要左上角红色的大部分为1，通过这种方式我们可以将4！个最优解变为单个最优解，且目标函数取值不变，这就是打破对称(Symmetry-Breaking Constraints)的基本思想。

此时整个代码如下，74-79为了与论文中形式一样做了修改，原理与之前基本相同

from docplex.mp.model import Modelfrom icecream import icimport numpy as npfrom multiprocessing import Process​​model = Model(name="The Deterministic OR Allocation")​​def generate_decision_variables(i, j):    '''    依次生成0-1决策变量x[j] ,y[i][j] 以及整数变量o[j]    i:手术数目    j:手术室数目    '''    x = model.binary_var_list(j, name="x")    o = model.integer_var_list(j, name='o')    varlist = [(ii, jj) for ii in range(i) for jj in range(j)]    y = model.binary_var_list(varlist, name='y')    return x, y, o​​if __name__ == "__main__":​​    # 输入所有超参数    i_ = 10         #手术(块)数目    j_ = 4         #手术室数目    c_f = 1         #手术室正常开放一天的固定消耗    c_v = 0.0333         #单个手术室加班一分钟的消耗    # d = [142, 276, 9, 211, 117, 223, 244, 333, 352, 94]    #手术块所需时长（min）    np.random.seed(1111110)    d = np.random.randint(0, 360, size=i_)    ic(d)    T = 480           #手术室正常开放时长（min）​​    #生成决策变量    x, y, o = generate_decision_variables(i_, j_)​​    #目标函数    model.minimize(np.sum(c_f * np.array(x))+ np.sum(c_v * np.array(o)))    # ic(np.sum(c_f * np.array(x))+ np.sum(c_v * np.array(o)))​​    #添加约束函数    for Yij,Xj in zip(y, x * i_):        model.add_constraint(Yij <= Xj)        # ic(Yij, Xj)​    for i in range(i_):        model.add_constraint(np.sum(y[j_ * i:j_ * i + j_]) == 1)        # ic(np.sum(y[j_ * i:j_ * i + j_]))​    for j in range(j_):        model.add_constraint(np.inner(d, y[j::j_]) <= T * x[j] + o[j])        # ic(np.inner(d, y[j::j_]) <= T * x[j] + o[j])​    for j in range(j_):        model.add_constraint(o[j] >= 0)        # ic(o[j] >= 0)​​    # 对称性破缺    for j in range(j_ - 1):        model.add_constraint(x[j] >= x[j+1])    #     ic(x[j] >= x[j+1])​    for j in range(j_):                         #已修改        model.add_constraint(np.sum(y[j * j_:][:j+1]) == 1)    #     ic(np.sum(y[j * j_:][:j+1]))​    for i in range(i_):        for j in range( j_):            if i>=j and i*j>=1 and i < j_:                model.add_constraint(np.sum(y[i * j_ + j : i * j_ + j_]) <= np.sum(y[(i - 1) * j_ + j - 1: : -j_][:j_ - j]))                # ic(i, j, y[i * j_ + j : i * j_ + j_])                # ic(y[(i - 1) * j_ + j - 1: : -j_][:j_ - j])​​    # #上下限    # L = (np.sum(d)) / (T * (1 + c_f / (c_v * T)))    # U = 2 * L    # model.add_constraint(np.sum(x) <= U)    # if j_ > L:    #     model.add_constraint(np.sum(x) >= L)    # ic(L, U)​​    # 求解    model.print_information()    sol = model.solve()    model.print_solution()    print(sol.solve_details)

接下来我们依次添加约束来收缩解空间（注意：以下所有的讨论都基于解为最优解的情况）：

（1）我们希望手术室依次开放，即如果4间手术室仅开放3间，应是1、2、3号。即当手术室编号j=1...m时，添加约束

对应代码中的65-68行。

    for j in range(j_ - 1):        model.add_constraint(x[j] >= x[j+1])    #     ic(x[j] >= x[j+1])

（2）不失一般性，我们假设手术台数m应大于手术室数目n。此时，手术0应安排在手术室0中，手术1应安排在手术室0或1中，手术2应安排在手术室0、1、2三者中，手术3应......。添加以下约束

代码为

    for j in range(j_):                         #已修改        model.add_constraint(np.sum(y[j * j_:][:j+1]) == 1)    #     ic(np.sum(y[j * j_:][:j+1]))

这条约束只能约束前m台手术，例如10台手术4个手术室情况下，如果前4台手术都安排在了第1手术室中，则后面手术的安排依然会带来对称性的解。所以我们考虑了下面这条约束。

（3）在上段例中，如果前4台手术已经被安排在第一手术室中，则第5台手术有两种可能①也被安排在第一手术室中；②被安排在其他手术室中。在情况②中，我们希望保证手术5被安排在2号手术室而不是3、4号中。

故我们添加了以下约束

简而言之为在针对顶部4x4的方阵中下图中每个颜色的方框，左侧框内之和都大于右侧(手绘比较丑，大家可以理解就好)

此处代码如下：

    for i in range(i_):        for j in range( j_):            if i>=j and i*j>=1 and i < j_:                model.add_constraint(np.sum(y[i * j_ + j : i * j_ + j_]) <= np.sum(y[(i - 1) * j_ + j - 1: : -j_][:j_ - j]))                # ic(i, j, y[i * j_ + j : i * j_ + j_])                # ic(y[(i - 1) * j_ + j - 1: : -j_][:j_ - j])

添加三个对称性破缺约束后，解的输出为：

objective: 6.764  x_0=1  x_1=1  x_2=1  x_3=1  o_1=7  o_3=76  y_0_0=1  y_1_1=1  y_2_2=1  y_3_1=1  y_4_2=1  y_5_3=1  y_6_0=1  y_7_3=1  y_8_2=1  y_9_0=1status  = integer optimal solutiontime    = 0.047 s.problem = MILPgap     = 0%

可见在目标函数不变的情况下，缩短了求解时间，将解在表格中表示如下

	手术室0	手术室1	手术室2	手术室3
手术0	1
手术1		1
手术2			1
手术3		1
手术4			1
手术5				1
手术6	1
手术7				1
手术8			1
手术9	1
加班		7		76

完美达成了我们想要求得的解，以上就是打破对称性破缺的约束，下篇文章会谈谈手术室开放个数的上下界，以及随机​规划模型。