手术的最优化分配（4）——随机规划模型及样本平均近似(SAA)方法代码实现

随机规划SAA方法的原理

经过这段时间艰难的边学边讲，我们终于推进到了论文的第五部分：随机规划模型，首先我们来看一下文中给出的数学模型。

之前的部分中，我们求解的是手术安排的确定性模型，但在实际的操作中，手术时长有很强的不确定性。如果不考虑手术时长的变动，则会给医院和患者都带来很大的风险。

在上述随机规划模型中，是否开放j号手术室x_j和是否将i手术安排到j手术室中y_i_j是第一阶段决策；j号手术室加班时间o_j(ω)是第二阶段的追索（补偿）决策。想了解关于随机规划和追索问题的基础内容请点击：“随机规划及其Recourse(追索/补偿）问题 (qq.com)”

我们实际求解上述随机规划的方式与文章中类似，为样本平均近似法(Simple Average Approximation，SAA)，将不确定问题看作k个确定性问题进行联合求解。即求出第一阶段的决策变量，并同时求出n个不同场景下的补偿变量o_j。

下面，给出三种不同的思路来考虑手术时长的不确定性。注意！三种思路中仅有一种是正确的，请各位读者思考：

(1)已知每台手术可能时长的上下限，假设手术时长为均匀分布。在十台手术情况下，假设此时为k=3个场景，则抽取3组d值如下表所示。

d_max	322	280	189	297	257	276	274	285	321	275
d_min	61	105	87	114	126	77	169	155	100	61
d1	210	151	148	167	149	245	235	233	229	180
d2	133	199	119	283	136	87	235	250	258	199
d3	120	123	187	167	242	163	214	207	235	85

此时，将d1、d2、d3中每列纵向求和并取均值d_mean,将均值直接代入确定性模型中求解。

d_mean=[154, 157, 151, 205, 175, 165, 228, 230, 240, 154]

    此种思路是错误的，因为其实我们并未考虑到不确定性。可以设想一个极端例子，假如只有两台手术均服从[10, 20]之间的均匀分布，第一个场景中手术时长均取为10，第二个场景中手术时长均取20，我们应该得到两个很极端的解，但采用这种思路只会得到一个平均的解。而且由于d_mean是对列求均值，其实我们并没有分别考虑到每种情况。

    (2)同样已知每台手术时长的上下限，抽取k=3个场景:

d_max	322	280	189	297	257	276	274	285	321	275
d_min	61	105	87	114	126	77	169	155	100	61
d1	210	151	148	167	149	245	235	233	229	180
d2	133	199	119	283	136	87	235	250	258	199
d3	120	123	187	167	242	163	214	207	235	85

    此时，将d1、d2、d3分别带入确定性模型中，得到三套x、y、o的取值，然后取平均数。

    此种思路同样是错误的，很简单，如果把所有的决策变量取均值，那0-1变量x，y的取值会出现小数，我们还是不知道应该如何分配手术室。

    （3）假设四间手术室，则分别设三种情况下的加班时长为o_1_1、o_1_2......o_3_4。如下表所示：

	X1=	X2=	X3=	X4=
y1
y2
y3
y4
y5
y6
y7
y8
y9
y10
加班时长	o11=	o12=	o13=	o14=
	o21=	o22=	o23=	o24=
	o33=	o32=	o33=	o34=

    上表中彩色部分为决策变量取值。蓝色部分为决策变量x，代表手术室是否开放；红色部分y_i_j代表手术i是否安排在手术室j中；紫色部分o_k_j代表第k种情况下第j号手术室的加班时长。

    此时，将目标函数变为

    即我们得出一套一阶段变量x、y，以及3种情况下分别的二阶段变量o。此时将约束条件

    变为：

d_1_1(第一种情况下第一台手术的时长)  * y_1_1(是否将手术1分配到手术室1中) + d_1_2 * y_2_1 + d_1_3 * y_3_1 + ... + d_1_10 * y_10_1 - o_1_1(第一种情况下1号手术室的加班时长) ≤ T * x_1 (1号手术室正常工作时长)

d_1_1 * y_1_2 + d_1_2  *  y_2_2+...+d_1_10 * y_10_2 - o_1_2 ≤ T * x_2

d_1_1 *  y_1_3.....

d_1_1  *  y_1_4...

d_2_1 * y_1_1 + d_2_2 * y_2_1 +...+ d_2_10  *  y_10_1 - o_2_1  ≤ T * x_1

..

..

..

d_3_1 * y_1_4 +...  + d_3_10 * y_10_4 - o_3_4  ≤ T * x_4

    共4间手术室x3种场景下12条约束。此处有一个小技巧，首先我们需要确定什么是一阶段决策变量，什么是二阶段追索变量。我们在建模时一阶段决策变量只有一种取值，二阶段变量会考虑到k种场景下的k种取值。

Python+DoCplex代码实现

 代码如下：

from docplex.mp.model import Model
from icecream import ic
import numpy as np
​
​
​
model = Model(name="The Stochastic OR Allocation")
def generate_decision_variables(n, m, k):
    '''
    依次生成0-1决策变量x[j] ,y[i][j] 以及整数变量o[j]
    n:手术数目    m:手术室数目    k:此时场景数    return x_j y_i_j o_k_j
    '''
    x = np.array(model.binary_var_list(m, name="x"))
    o = np.array(model.integer_var_list([(kk, j) for kk in range(k) for j in range(m)], lb=0, name='o')).reshape(k, m)
    y = np.array(model.binary_var_list([(i, j) for i in range(n) for j in range(m)], name='y')).reshape(n, m)
    return x, y, o
​
​
if __name__ == "__main__":
    # 输入所有超参数
    n = 10        #手术(块)数目(i)
    m = 4         #手术室数目(j)
    k = 3  # 场景数
    c_f = 1         #手术室正常开放一天的固定消耗
    c_v = 0.0333         #单个手术室加班一分钟的消耗
    #手术块所需时长（min）
    np.random.seed(1111110)
    D_max = np.tile(np.random.randint(180, 360, size=n), k).reshape(k, n)
    D_min = np.tile(np.random.randint(60, 180, size=n), k).reshape(k, n)
    ic(D_min[0], D_max[0])
    # 与确定性模型取值相同，用于测试解的质量
    # D_min = np.tile([142, 276, 9, 211, 117, 223, 244, 333, 352, 94], k).reshape(k, n)
    # D_max = np.tile([143, 277, 10, 212, 118, 224, 245, 334, 353, 95], k).reshape(k, n)
    # 假如为正态分布
    # D = [142, 276, 9, 211, 117, 223, 244, 333, 352, 94] * k
    # d = np.abs(np.random.normal(D, scale=20). reshape(k, m)).astype(int)
    # 假如为在dmin-dmax之间的均匀分布
    d = np.random.randint(D_min, D_max)
    ic(d)
    T = 480           #手术室正常开放时长（min）
    #生成决策变量
    x, y, o = generate_decision_variables(n, m, k)
    # ic(x, y, o)
    # 目标函数
    model.minimize(c_f * np.sum(x) + c_v * np.sum(o) / k)
    # ic(c_f * np.sum(x) + c_v * np.sum(o))
    # 添加约束函数
    for Yij,Xj in zip(y.reshape(n*m), np.tile(x, n)):
        model.add_constraint(Yij <= Xj)
        # ic(Yij <= Xj)
    for i in range(n):
        model.add_constraint(np.sum(y, axis=1)[i] == 1)
        # ic(np.sum(y, axis=1)[i] == 1)
    for kk in range(k):
        for j in range(m):
            model.add_constraint(np.inner(y[:, j], d[kk]) - o[kk][j] <= x[j] * T)
            ic(kk, np.inner(y[:, j], d[kk]) - o[kk][j] <= x[j] * T)
    # 求解
    model.print_information()
    sol = model.solve()
    print(sol.solve_details)
    model.print_solution()