手术的最优化分配（1）——确定性模型及Python调用Docplex求解

本文讲Brian T.Denton在2010年发表于Operations Research上的文章《Optimal Allocation of Surgery Blocks to Operating Rooms Under Uncertainty》,DOI号为“10.1287/opre.1090.0791”。笔者将标题译为《不确定条件下手术（块）在手术室中的最优化分配》。文中标下划线处为翻译内容，其余为讲解。

文章关键词：优化、随机规划、手术。

摘要

众所周知，一篇论文最重要的部分就是摘要。该论文摘要包括三个部分，分别为：问题介绍、文章内容和本文结论。

问题介绍：手术在不同手术室之间的分配是一个很有挑战的组合优化问题，手术时长的高度不确定性更大大增加了分配问题的复杂度。

内容：本文提出了一个随机规划模型，用以在特定某天内将手术分配到不同手术室中。目标函数包括开放手术室的固定费用和可变的加班费用。本文描绘了两种模型：

①两阶段随机规划：一阶段——二元决策变量；二阶段——追索问题。想了解追索问题可以点击“随机规划及其Recourse(追索/补偿）问题 (qq.com)”

②：对应的鲁棒优化模型：目标为在手术时长的不确定集上最小化最大成本。简单来说，鲁棒优化就是针对最坏情况的优化。

分别描述了1、数学模型，2、最优解约束，3、解法，还有4、易于实施的启发式算法。本文还进行了基于真实医疗数据的数值实验，用于比较两个模型的优劣，并且展示了模型求解现实问题的潜力。

结论：基于数值实验，我们发现启发式算法在所有情况下表现均匀，速度快并且易于实施。

我们也发现了鲁棒模型解的质量和启发式算法一样好，比随机优化模型收敛更快，且更有利于在追索问题中限制最差的结果。

确定性模型

读完摘要部分我们已经对文章框架有了大概的了解，接下来我们直接跳过背景介绍和综述部分进入正题：确定性手术室分配模型（The Deterministic OperationRoom Allocation Problem,简称DORA)

手术分配问题中包括两个重要决策：

①当日开放多少手术室(Operation Room,简称OR);

②手术在手术室之间如何分配。

手术室开放决策①和手术分配决策②的目标为最小化开放手术室以及手术室加班的总成本。决策①和决策②是对抗/矛盾关系，因为：

• 少开放手术室会导致更多的加班成本；

• 增加手术室可以降低加班时长，但会增加手术室成本。

 我们将开放手术室的成本看作一个固定值，该成本基于①支持手术室开放必须的人员；②上下游资源（例如护士和麻醉病房）。手术室开放的固定成本假设手术室只能开放一整天或是不开。在经典的可拓展装箱问题模型中手术室(集装箱）的数目是固定的，与之不同，我们把开放手术室的数目当作决策变量。

在DORA模型中，主要有以下变量：

i	手术编号，i = 1, …, n
j	手术室编号, j=1, …, m
di	手术i所需时间(按分钟计算）
T	手术室的正常开放时长(即医护人员正常上班时长)
cf	开放手术室的固定成本
cv	单个手术室超时开放单位时间的成本(可以理解为一个手术室加班一分钟的成本)
xj	0-1变量，代表手术室OR是否开放
yij	0-1变量，代表手术i是否安排到手术室j中
oj	决策变量，代表手术室j的加班时长(按分钟计算)

根据以上参数，确定性手术分配数学模型如下：

请大家记住，该模型仅针对择期手术患者，换而言之我们实现已经决定了哪些手术要在今天完成，现在的任务是将手术分别安排在手术室中。

接下来，我们用自然语言解读一下数学模型：

目标函数(1)为求总成本的最小值,c^f * x_j为手术室开放的固定成本乘以手术室开放数目，两者相乘再对手术室编号j求和 即为所有手术室开放的固定成本；c^v * o_j为加班成本乘以加班时长，两者相乘对手术室编号j求和即为所有手术室加班的成本。(此处和下文的c^f代表f是c的上标，x_j代表j为x的下标，是latex语法。因为微信的文章编辑器不支持上下标，大家可以理解就好)总成本=固定成本+加班成本

约束函数(2)的意义为：只有手术室j开放时，才能将手术i安排在手术室j中。因为x，j均为0-1变量，所以此处可以理解为只有x_j=1时，y_i_j才能为一。

假设i=10，j=4，即十台手术四个手术室时，只有x_4=1 即4号手术室开放时，y_10_4才有可能等于1，即才可能将10号手术安排在4号手术室中。

约束函数(3)为，对于任意手术i，y_i_j对手术室j求和结果均为1，即每个手术i都必须安排在手术室中，可以理解为所有手术都必须在当日完成。

约束函数(4)，不等式左边的值为该手术室今日手术的总时长，即该手术室安排的手术y_i_j乘以手术时长d_i；不等式右侧的值为该手术室的正常开放时长T乘以开放手术室的决策变量x_j，加上该手术室的加班时长o_j。(4)约束意义为，对于每一个手术室j，手术时长之和要小于该手术室的上半时长+加班时长。

一般来说，仅当x_j==1时该约束有意义，x_j==0时手术室不开放，该约束等价为约束(2).

约束(5)则很易于理解，所有x，y均为0-1变量，即取值只能为0或1.加班时长o_j必须>=0，不能为负值。

为了不失一般性，模型中c_f<c_v * T，即每单位的手术室正常开放消耗会小于每单位的加班消耗。否则我们会得到很荒谬的解，类似于单个手术室大量加班，其他手术室不予开放。

Docplex库介绍

以下划线内容引用自“https://zhuanlan.zhihu.com/p/69400759”

首先需要了解IBM ILOG CPLEX Optimization Studio还提供两个优化器：CPLEX和CP

• CPLEX：用于Mathematical Programming（数学规划）--多基于线性规划理论

• CP：Constraint Programming (基于约束的规划) --多基于图论的理论

IBM ILOG CPLEX Optimization Studio 有两个Python API库：

• CPLEX Python API

• DOCplex

以上两个库的差别如下

• CPLEX Python API：

  • 是遗留传统的Python API.

  • Callable Library的应该原生C语言的

  • 目前仅用于CPLEX优化器，即Mathematical Programming（数学规划）领域。

  • 同时这个API和 其他C++, Java, and .NET APIs一样，都是C Callable Library在其他各领域语言中提供的API。

  • 只能基于本地的IBM ILOG CPLEX Optimization Studio软件及其优化引擎求解优化问题

  • 由于这只是一个领域语言的API接口库，因此使用该API库，还需要依赖原始软件或优化引擎的安装，就还需要安装IBM ILOG CPLEX Optimization Studio软件，提供软件和C Callable Library库和优化引擎。

DOCplex：又叫DOcplex Python Modeling API 或 Decision Optimization CPLEX Modeling for Python

• 原生Python的优化建模API库

• 和Python的Numpy，pandas库兼容

• 同时支持CPLEX和CP优化引擎

• 目前可以对接本地安装的IBM ILOG CPLEX Optimization Studio软件及其优化引擎，或者支持对接云端DOCloud（ IBM® Decision Optimization on Cloud service ）

• 对接云端的DOCloud，则需要注册DOCloud的账号，并获得相应的API Key

  个人体验上讲，Docplex库使用比cplex库友好得多，而且对np的支持很大程度上简化了建模难度（后面会讲到）。DOcplex库的编程思路与cplex类似，均为

①输入决策变量；

②输入目标函数；

③输入约束方程。

python建模及求解

了解了手术分配的确定模型后，我们自然想到了将模型在计算机中表达并使用真实数据进行求解。

首先给大家展示整个代码↓，然后进行逐行的讲解，请记住：代码和上述的确定性模型是严格一一对应的，模型中每一个方程都顺序对应着一个代码块。

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 = [152, 74, 218, 127, 183, 175, 224, 266, 461, 310]    #手术块所需时长（min）    # np.random.seed(2021)    # d = np.random.randint(0, 480, 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(np.sum(d)//T + 1):        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:                model.add_constraint(y[i*j_+j] <= np.sum(y[(j-1)*j_+j-1::j_][:i-j+1]))    #             ic(i,j,y[i*j_+j])    #             ic(np.sum(y[(j-1)*j_+j-1::j_][:i-j+1]))

    #上下限    # 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行为引用该程序所需的库；10-20行根据输入的i，j值生成了决策函数x、y和o；26-35行输入了此次所用的求解参数；39行调用10行函数生成决策变量；43行是模型的目标函数；47-62行中的四个循环分别实现了上一章模型中(1,2,3,4)四个约束条件；65-88行是为了提高求解效率，加入与否不影响最优解的值，我们会在下篇文章中进行讲解；92-95行是模型的求解和输出部分，将该模型运行输出↓如下：

Model: The Deterministic OR Allocation - number of variables: 48   - binary=44, integer=4, continuous=0 - number of constraints: 90   - linear=90 - parameters: defaults - objective: minimize - problem type is: MILPobjective: 12.991  x_0=1  x_1=1  x_2=1  x_3=1  o_0=157  o_1=55  o_2=4  o_3=54  y_0_0=1  y_1_1=1  y_2_2=1  y_3_3=1  y_4_3=1  y_5_0=1  y_6_3=1  y_7_2=1  y_8_1=1  y_9_0=1status  = integer optimal solutiontime    = 0.015 s.problem = MILPgap     = 0%

可以看到模型中当i=1,...,10;j=1,...,4时模型共有48个决策变量(第2行),90个约束条件(第4行),该问题为混合整数线性规划(Mixed Integer Linear Programming，MILP）问题(第8行)。

第9行输出为目标函数最优值12.991，10-27行为所有决策变量的取值，求解时间为0.015s(第29行），笔者PC运行求解时间一般在0.02s以下。

接下来我们详细地讲解代码：

• 1-7行

from docplex.mp.model import Modelfrom icecream import icimport numpy as npfrom multiprocessing import Process

model = Model(name="The Deterministic OR Allocation")

1-4行为导入该模型需要的库，Model为docplex库中数学建模模块，ic是一个方便调试的输出库，功能与python自带print函数类似，但在输出时会附上变量名称，代码中的34,44,50,54,58,62行注释（其他ic也一样）均可变为可执行代码输出，可以清楚地看到我们给模型中增加了哪些变量和约束条件；numpy为一个矩阵运算库，docplex库与numpy库完美兼容，后续也会讲到；Process为多线程库，可以同时调用cpu多个核用以提高计算效率，在该代码中并未实际用到。

第七行为实例化(创建)一个docplex模型，命名为“The Deterministic OR Allocation”。

• 10-20行

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      ic(generate_decision_variables(3,2))

以上代码我们创建了一个函数，输入手术数目i和手术室数目j即可生成0-1决策变量x_j(分别代表1..j手术室是否开放),0-1决策变量y_i_j(代表手术i是否安排在手术室j中），以及整型决策变量o_j(代表手术室j的加班时长)。

上面这个代码块中，7行是生成单下标0-1变量，8行是生成单下标整数型变量，9行是生成双下标0-1变量。

13行是为了讲解临时添加，即当3台手术在2间手术室之间分配时，变量输出如下：

ic| generate_decision_variables(3,2): ([docplex.mp.Var(type=B,name='x_0'), docplex.mp.Var(type=B,name='x_1')],                                       [docplex.mp.Var(type=B,name='y_0_0'),                                        docplex.mp.Var(type=B,name='y_0_1'),                                        docplex.mp.Var(type=B,name='y_1_0'),                                        docplex.mp.Var(type=B,name='y_1_1'),                                        docplex.mp.Var(type=B,name='y_2_0'),                                        docplex.mp.Var(type=B,name='y_2_1')],                                       [docplex.mp.Var(type=I,name='o_0'), docplex.mp.Var(type=I,name='o_1')])

由上述输出可知，3台手术2间手术室对应着2个B(Binary)二元变量x,分别对于两间手术室是否开放，两个I(integer)变量o对应两个手术室的加班时长，3×2个变量y对应手术的分配。

• 23-39行：输入超参数

if __name__ == "__main__":
    # 输入所有超参数    i_ = 10         #手术(块)数目    j_ = 4          #手术室数目    c_f = 1         #手术室正常开放一天的固定消耗    c_v = 0.0333         #单个手术室加班一分钟的消耗    d = [152, 74, 218, 127, 183, 175, 224, 266, 461, 310]    #手术块所需时长（min）    # np.random.seed(2021)    # d = np.random.randint(0, 480, size=i_)    # ic(d)    T = 480           #手术室正常开放时长（min）

    #生成决策变量    x, y, o = generate_decision_variables(i_, j_)

第1行标志着主程序的开始,其余行均为可以随意更改的参数，c_v=0.0333,c_f=1取值意味着手术室加班30分钟成本相当于正常开放一整天。

8行和9-11行只需要运行两者之一，另一个需转变为注释，8为手动输入的手术时长，因为手术时长列表d需要与手术数目i_相对应，所以当改变i_时应把8转变为注释，把9-11变为代码执行。

其中9为随机数种子，确保10生成的结果可复现；

10为生成包含i_个位于[0,480]区间内随机元素内的列表；

11行为输出此时的d取值；

16行调用之前设定的参数生成决策变量。

• 43-44行：目标函数（1）

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)))

np.sum表示对数组求和，n.sum(1,2,3)值为1+2+3=6。

k*np.array()为用常数k分别乘()内的每一个数，3*np.array([1, 2, 3]) = [1*3, 2*3, 3*3]=[3, 6, 9]。

输出如下

ic| np.sum(c_f * np.array(x))+ np.sum(c_v * np.array(o)): docplex.mp.LinearExpr(x_0+x_1+x_2+x_3+0.033o_0+0.033o_1+0.033o_2+0.033o_3)

由输出可知，此时我们的目标函数为取下式的最小值，与数学模型中(1)完美对应。

1*(x_0+x_1+x_2+x_3)+0.0333*(o_0+o_1+o_2+o_3)

这里插一句题外话，docplex库与numpy库的兼容性好得让人吃惊，nump.sum()函数是对()里的数组求和，()中的对象应该有明确的取值，此处的模型还未求解，x_0、x_1、x_2、,o_0、o_1、o_2均为未知数，此处竟然可以运行并输出。不仅如此，本文中用到的其他np函数都完美兼容docplex，这给模型的搭建过程降低了很多难度，我们可以明确地看到在哪一步中添加了什么样的约束条件。

• 48-50行：约束条件（2）

    for Yij,Xj in zip(y, x * i_):        model.add_constraint(Yij <= Xj)        ic(Yij, Xj)

此处有一个建模小技巧，当我们看到对于任意(i,j)时，就应该想到(i,j)共有40种不同取值，此处应该有40个约束条件。对于其他约束条件也类似，抓住这一点可以快速理清思路。

因为y有40个变量，x仅有10个，所以我们将x拓展到与y同样的维度添加约束，x，y输出如下：

ic| Yij: docplex.mp.Var(type=B,name='y_0_0')    Xj: docplex.mp.Var(type=B,name='x_0')ic| Yij: docplex.mp.Var(type=B,name='y_0_1')    Xj: docplex.mp.Var(type=B,name='x_1')ic| Yij: docplex.mp.Var(type=B,name='y_0_2')    Xj: docplex.mp.Var(type=B,name='x_2')ic| Yij: docplex.mp.Var(type=B,name='y_0_3')    Xj: docplex.mp.Var(type=B,name='x_3')ic| Yij: docplex.mp.Var(type=B,name='y_1_0')    Xj: docplex.mp.Var(type=B,name='x_0')ic| Yij: docplex.mp.Var(type=B,name='y_1_1')    Xj: docplex.mp.Var(type=B,name='x_1')ic| Yij: docplex.mp.Var(type=B,name='y_1_2')    Xj: docplex.mp.Var(type=B,name='x_2')ic| Yij: docplex.mp.Var(type=B,name='y_1_3')    Xj: docplex.mp.Var(type=B,name='x_3')ic| Yij: docplex.mp.Var(type=B,name='y_2_0')

因为输出很长，此处就不全部展示了，大家可以看到xj与yij的j取值完全相同，这一段函数完美对应了模型中的约束(2)

• 51-53行：约束条件(3)

    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_]))

因为(3)条件为对于任意i，所以应该有10个约束条件，y有i个手术，j个手术室共计i*j个变量。对于任意i，y_i_j对j求和均为1,所以取从i到i+j的y进行切片后相加等于1，输出如下。

ic| np.sum(y[j_ * i:j_ * i + j_]): docplex.mp.LinearExpr(y_0_0+y_0_1+y_0_2+y_0_3)ic| np.sum(y[j_ * i:j_ * i + j_]): docplex.mp.LinearExpr(y_1_0+y_1_1+y_1_2+y_1_3)ic| np.sum(y[j_ * i:j_ * i + j_]): docplex.mp.LinearExpr(y_2_0+y_2_1+y_2_2+y_2_3)ic| np.sum(y[j_ * i:j_ * i + j_]): docplex.mp.LinearExpr(y_3_0+y_3_1+y_3_2+y_3_3)ic| np.sum(y[j_ * i:j_ * i + j_]): docplex.mp.LinearExpr(y_4_0+y_4_1+y_4_2+y_4_3)ic| np.sum(y[j_ * i:j_ * i + j_]): docplex.mp.LinearExpr(y_5_0+y_5_1+y_5_2+y_5_3)ic| np.sum(y[j_ * i:j_ * i + j_]): docplex.mp.LinearExpr(y_6_0+y_6_1+y_6_2+y_6_3)ic| np.sum(y[j_ * i:j_ * i + j_]): docplex.mp.LinearExpr(y_7_0+y_7_1+y_7_2+y_7_3)ic| np.sum(y[j_ * i:j_ * i + j_]): docplex.mp.LinearExpr(y_8_0+y_8_1+y_8_2+y_8_3)ic| np.sum(y[j_ * i:j_ * i + j_]): docplex.mp.LinearExpr(y_9_0+y_9_1+y_9_2+y_9_3)

y的10台手术分别对4个手术室进行求和，对应了原模型中的约束(3)

• 55-57行：约束条件（4）

    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])

(4)条件为对于所有j，所以共有四个约束条件，对于每一个手术室j，手术时长求和<开放时长+加班时长。np.inner()为求两个数组的内积，例如

np.inner([1, 2, 3], [4, 5, 6]) = 1 × 2 + 2 × 5 + 3 × 6

y[j::j_]为从第j个开始，以步长j_对y进行切片，例如当j=1,则y[j::j_]包括了y_1_1 + y_2_1 + y_3_1....

上述函数输出如下：

ic| np.inner(d, y[j::j_]) <= T * x[j] + o[j]: docplex.mp.LinearConstraint[](152y_0_0+74y_1_0+218y_2_0+127y_3_0+183y_4_0+175y_5_0+224y_6_0+266y_7_0+461y_8_0+310y_9_0,LE,480x_0+o_0)ic| np.inner(d, y[j::j_]) <= T * x[j] + o[j]: docplex.mp.LinearConstraint[](152y_0_1+74y_1_1+218y_2_1+127y_3_1+183y_4_1+175y_5_1+224y_6_1+266y_7_1+461y_8_1+310y_9_1,LE,480x_1+o_1)ic| np.inner(d, y[j::j_]) <= T * x[j] + o[j]: docplex.mp.LinearConstraint[](152y_0_2+74y_1_2+218y_2_2+127y_3_2+183y_4_2+175y_5_2+224y_6_2+266y_7_2+461y_8_2+310y_9_2,LE,480x_2+o_2)ic| np.inner(d, y[j::j_]) <= T * x[j] + o[j]: docplex.mp.LinearConstraint[](152y_0_3+74y_1_3+218y_2_3+127y_3_3+183y_4_3+175y_5_3+224y_6_3+266y_7_3+461y_8_3+310y_9_3,LE,480x_3+o_3)

其中"LE"意为≤,所以可以清楚地看出，代码对应了模型中的约束条件(4)

• 59-61行：约束条件（5）

    for j in range(j_):        model.add_constraint(o[j] >= 0)        ic(o[j] >= 0)

这段代码则很好理解，原约束(5)中的x，y取值范围已经在生成决策变量时规定好，所以这里只需添加手术时长o_j的非负约束，代码输出如下：

ic| o[j] >= 0: docplex.mp.LinearConstraint[](o_0,GE,0)ic| o[j] >= 0: docplex.mp.LinearConstraint[](o_1,GE,0)ic| o[j] >= 0: docplex.mp.LinearConstraint[](o_2,GE,0)ic| o[j] >= 0: docplex.mp.LinearConstraint[](o_3,GE,0)

“GE”代表≥，所以这里对应了原模型中约束(5)的后半部分

• 91-94行：模型求解与输出

    model.print_information()    sol = model.solve()    model.print_solution()    print(sol.solve_details)

此段第1行为输出模型基本信息例如约束和变量个数；第2行为求解模型；第3行为输出所有决策变量x,y,o的最终取值；第4行为输出求解的时间和是否最优解。

结论

因为代码调试过程中和原数学模型做了严格的对应，所以整段代码不会有较大的问题。本文跳过了64-87行代码的讲解，而这一部分代码主要用于加速解的收敛。

在原代码10台手术4个手术室的条件下，不加64-87代码的求解时间为0.078秒，而增加之后在最优解不变的情况下，求解时间缩短为0.015s，大大缩短了求解时间。

而在手术数目i=100，手术室j=40的条件下，不加64-87代码求解时间为2.984s，增加了之后在最优解不变的条件下，求解时间缩短为0.282秒，提升明显。

笔者尝试运行1000台手术，400个手术室的情况，64-87代码运行的条件下，模型共有400800个变量，721898个线性约束，求解时间为306.468 s.总体来说表现优秀。

根据上述例子可以看出64-87行代码提升求解效率非常明显，其中81-87行是决策变量取值的上下限，由数学推导得出；而64-78行打破“对称性破缺”，原理是在不改变最优解的情况下增加约束函数，大大收缩解空间，提高求解速率。具体原理会在下一篇文章中讲解。

p.s. 在线性模型求解的过程中，约束函数越多目标函数收敛得越快，因为每一个约束函数都在超平面上缩小了解空间。可以类比现实中求解多元方程组，如果方程足够多，我们只需要做简单的加减就可以把一些未知数变为确定值。

不足

这段代码在模型的求解速率之外，为了方便调试和讲解使用了大量的循环来添加决策变量和约束条件，如果想要进一步提高代码运行速率应该对循环函数合并，以减少循环次数。