手术的最优化分配（3）——手术室开放数目上下限及手术分配的启发式算法

在上期“手术的最优化分配（2）——打破对称性的约束(Symmetry-Breaking Constraints) (qq.com)”的代码中，没有对82-86行代码进行讲解，今天我们先来谈谈这几行代码背后的原理。

手术室开放数目的下限

最优的手术室开放数目的下限为：

证明过程如下：

手术室分配的前置条件为医院已经决定了哪些手术需要在今天完成，我们的任务为将手术分配在手术室中。

在之前的文章中，我们谈到c^f为手术室正常开放一整天的固定成本，c^v为手术室加班一分钟的成本。在我们用以计算的例中，取值为

c^f = 1, c^v = 0.0333，

此时 

c^f / c^v = 30

即手术室加班30min的成本相当于正常开放一天。因此

（1）在手术室未全部开放的情况下：任何一个手术室的总开放时长不应超过T+c^f/c^v,否则此时可以通过开放一间新的手术室来降低成本。

（2）在手术室全开放的情况下，不可能所有手术室的工作时长都超过T+c^f/c^v，否则此时加班成本已经超过了第二天正常上班的成本。这种情况下手术的安排首先就出现了问题。

所以手术室开放的数目的下限为

        分子部分为所有手术的总时长，分母部分为单个手术室开放的最大时长。L即为手术室的最小开放数目。

手术室开放数目的上限

最优的手术室开放数目的上限为：

    我们可以发现，手术室开放数目的上限恰好为其下限的二倍。这是因为如果同时有两间手术室的手术时长小于T+c^f/c^v，那么可以将其合并在一间手术室完成来降低成本。此时虽然有可能产生加班成本，但减少了一间手术室开放所带来的固定成本。故手术室开放数目的上限为下限的二倍。

一个简单的启发式算法

此处我们用到的基本算法为最长处理时间法（Longest processing time，LPT），以下划线内容翻译自“https://riot.ieor.berkeley.edu/Applications/Scheduling/algorithms.html”。

LPT规则按照处理时间对作业进行排序，无论何时释放一台机器，此时就绪的最长作业就开始处理。是一种启发式的算法，用于求解调度问题。它首先安排时间最长的作业，这样就不会有大型作业在计划末尾“突出”。

在我们之前的例子中，共有4间手术室，10台手术，手术时长为[142, 276, 9, 211, 117, 223, 244, 333, 352, 94]，假如四间手术室全部开放，则安排的顺序如下

首先将手术时长从大到小排序，为[352, 333, 276, 244, 223, 211, 142, 117, 94, 9]；然后依次将手术安排到总时长最小（即最先释放）的手术室中，此时得到的解为加班时长为88min，总成本为6.93；略大于我们数学模型求得的加班时长83min，总成本6.764。可见我们此时得到的解虽然为次优解，但表现还算不错。

但是在实际的手术分配问题中，我们并不知道应该开放几间手术室，因此我们结合推导得出的手术室上下限，得出了启发式的算法如下：

原理非常简单，我们从下限L到上限U中依次遍历手术室的开放数目，然后使用LPT法求出此时的成本z，在所有z之中保存一个Zmin作为输出的最后解。

上述的启发式算法是一个经典的可拓展装箱问题的解法，Dell’Olmo在1998年的论文《

A 13/12

approximation algorithm for bin packing with extendable bins

》中证明了此解法求出的解≤最优解的1.0833倍（13/12倍)。

接下来使用代码实现启发式算法：

import numpy as np
from icecream import ic
​
​
def LPT(i, d):
    '''
    使用LPT法计算此时的手术分配
    i:此时开放的手术室数目
    d:手术块时长
    return: 此时的成本z
    '''
    d = sorted(d)[::-1]
    assign = np.zeros((n+1, i))
    for nn in range(n):
        for mm in range(i):
            assign[-1] = np.sum(assign[:-1], axis=0)
            if assign[-1][mm] == np.min(assign[-1]):
                assign[nn][mm] = d[nn]
                break
    return assign
​
​
def Bound(T, d, c_f, c_v):
    '''
    计算手术室开放数目的上下限
    T: 手术室正常开放时长
    d: 每个手术所需的时长
    c_f: 手术室开放一天的固定成本
    c_v: 单个手术室加班一分钟的成本
    return: 手术室开放下限L，手术室开放上限U
    '''
    L = int((np.sum(d)) / (T * (1 + c_f / (c_v * T))))
    U = int(2 * L)
    return L, U
​
​
def cost(assign, i):
    z = i*c_f + np.sum(np.maximum(assign[-1] - 480, 0)) * c_v
    return z
​
if __name__ == "__main__":
​
    n = 10
    m = 4
    c_f = 1
    c_v = 0.0333
    T = 480
    d = [142, 276, 9, 211, 117, 223, 244, 333, 352, 94]
    # np.random.seed(2021)
    # d = np.random.randint(30, 360, size=n)
    ic(d)
​
    L, U = Bound(T, d, c_f, c_v)
    # ic(L,U)
​
    z = {}
    for i in range(L, min(U, m) + 1):
        assi = LPT(i, d)
        z[i] = cost(assi, i)
        # ic(assi)
        # ic(cost(assi, i))
    ic(z)
​

代码部分不进行详细讲解了，毕竟启发式算法的原理非常简单，每个人都可以用自己的风格编写出来。