CIKM 2021 | UCL、华为、复旦提出超大规模建模系统Grassland，助力华为供应链快速决策

基于数学规划方法进行商业决策的流程图，虚线框内的部分即为代数建模系统

对于一组决策变量，一个完整的数学规划模型由目标f(x)和约束两部分组成，即. 在实际商业场景中，目标f(x)对应于希望达成的商业目标（例如，最大化订单满足率，最小化成本），而约束对应于需要遵守的商业规则和现实限制（例如，仓库的出库量与入库量之差不能大于库存量）。在大规模场景中，虽然目标往往只有一行表达式，但约束的数量却可以数以百万计，从而成为建模过程的主要瓶颈。我们接下来以经典的网络流模型为例介绍这一问题。

如图所示，可以将有向图想象成水管网络，源点 1 和汇点 4 想象成水管网络的供水口和出水口，对于任意一条有向边(i,j)，有流量，代表从顶点 i 流到 j 的水流量。网络流问题需要满足一个自然的约束，即对于所有中间顶点，其总流出量应当等于总流入量。而对于源点和汇点，其总流出量和总流入量应当等于一个定值，即网络的总流量（简单起见令其为 1）。我们将其写作

i=1 时，找到所有满足和的元素，分别为 (1,2), (1,3) 和空集，因此有

i=2 时，找到所有满足和的元素，分别为 (2,4) 和(1,2), (3,2)，因此有

i=3 时，找到所有满足和的元素，分别为 (3,2), (3,4) 和(1,3)，因此有

i=4 时，找到所有满足和的元素，分别为空集和(2,4),(3,4)，因此有

然而容易看出，以上过程需要对顶点集合V中的每一个元素进行迭代，每次迭代时需要遍历边集E，本质是一个二重循环，时间复杂度是O(|V||E|) 的！如果我们有数以百万计的顶点和边，这种生成约束的方法就完全行不通了。虽然针对这里的简单例子可以有一些加速方案（例如对边集E做预先索引），但通用性差，需要根据问题定制，难以根本性地解决复杂商业场景建模效率低下的问题。

我们使用红色和蓝色分别标记出第一个求和表达式（总流出量）和第二个求和表达式（总流入量）在四个约束中所对应的变量项x。我们会发现，它们恰好都在四个约束中分别出现了 5 次，且完整对应于边集数据E。这是一个巧合吗？并非如此！以总流出量为例，求和符号的条件是，而作为约束的下标，会遍历所有可能的取值。因此当 i 的遍历完成后，E也会被不多不少恰好遍历一次，使得四个约束里一定会出现 5 个与E一一对应的变量项x。

2. 无视表达式后的约束下标，直接对第一个和第二个求和表达式中的分别进行遍历，得到总计 10 个变量项。同时，每个变量项对应的约束下标值进行额外标记（这里使用红色标出）。

以上步骤其实很类似于一个 “种草” 的过程。如果我们将线性约束的约束矩阵视作一片土地（Land），那么我们其实是在第 1、2 步批量生成了所有要种在这片土地上的草苗（Grass），然后在第 3、4 步将这些草苗一棵棵地种在了约束矩阵A的适当位置。这也是本系统被命名为 “Grassland” 的由来。

不难看出，如果我们将第三步哈希表归类的单次时间复杂度视为O(1)，那么整个约束建模过程的时间复杂度是的！也就是说，我们仅仅是将边集数据E扫描了常数次，就完成了约束的建模过程。同时值得注意的是，这种方案的通用性极高。只要约束是线性的（绝大多数可大规模求解的模型都要求约束为线性），且求和表达式的条件都写作 “下标∈数据” 的格式，就可以使用这一方法进行约束建模，在华为供应链不同场景下数十种超大规模复杂约束下均表现出色。论文中详细讨论了上述建模过程和直接按照公式建模的等价性。

为例。假设我们有两个 CPU 核，可以互不干扰地并行处理子任务，那么是否有一种方法可以将建模任务平均分配到两个 CPU 上呢？从公式的字面意义来看，这件事情并不难。我们只需将顶点集合平分为，并让两个 CPU 核分别处理，生成和，最后合并即可，我们将这个过程称为“对约束空间的切分”。然而，正如我们前面的分析所示，大规模场景下按照公式字面意义进行建模是完全不实际的，即使能将任务平分到多个 CPU 上也不例外。

import pulp
V = [0, 1, 2, 3]
E = [(0, 1), (0, 2), (1, 3), (2, 1), (2, 3)]
S = [1, 0, 0, -1]
prob = pulp.LpProblem()
x = pulp.LpVariable.dicts('x', (range(len(V)), range(len(V))))
for v in V:    
    flow_in = pulp.lpSum(x[i][j] for i, j in E if i == v)    
    flow_out = pulp.lpSum(x[i][j] for i, j in E if j == v)    
    prob += flow_in - flow_out == S[v]
print(prob)

import grassland as gl
max_V = 4
model = gl.Model()
x = gl.Variable(name='x', shape=(max_V, max_V))
i = gl.IndexPlaceholder(range=(0, max_V))
E_in = gl.Constant(shape=(max_V, max_V))
E_out = gl.Constant(shape=(max_V, max_V))
S_ = gl.Constant(shape=(max_V,))
flow_out = gl.cond_sum(index_placeholders(i, 'j'),cond=E_out,  expr=gl.Expression(x[i, 'j']))
flow_in = gl.cond_sum(index_placeholders=('j', i), cond=E_in,expr=gl.Expression(x[i, 'j']))
model.add_constraint(gl.constraint(expr=flow_in - flow_out, sense='=', rhs=S,axes=[i], parallel_axes[i]) )
V = gl.COO([0, 1, 2, 3])
E = gl.COO([(0, 1), (0, 2), (1, 3), (2, 1), (2, 3)])
E_T = gl.COO([(1, 0), (1, 2), (2, 0), (3, 1), (3, 2)])
S = gl.COO([1, 0, 0, -1])
model.generate_lp(feed_dict={E_out: E, E_in: E_T, S_: S},     num_workers=2, filename='network_flow.lp')

接下来，构建总流入量 flow_in 和总流出量 flow_out 表达式，这一部分无需对约束下标 i 的所有可能取值进行循环遍历，而是使用 Grassland 中的计算图节点 gl.cond_sum 构造出表达式和。gl.cond_sum 包含三个参数：下标 index_placeholders、数据 cond 和表达式，即

下标 index_placeholders：传入一个 Python 元组，依次列出与数据对应的下标。其中，约束中出现的下标（例如前文中使用红色标识的下标 i，在论文中称为全局下标，global index placeholder）使用 gl.IndexPlaceholder()定义，而其他在求和符号中局部使用的下标（例如网络流约束中的下标 j，在论文中称为局部下标，local index placeholder）使用 python 字符串定义。例如，可以传入描述的下标部分（红色标出）。

数据 cond：传入一个 gl.Constant 对象，作为计算图的占位符以供后续传入实际数据。值得注意的是，Grassland 使用稀疏张量（Sparse Tensor）来表示数据。例如，E = [(0, 1), (0, 2), (1, 3), (2, 1), (2, 3)] 可以被视为一个 4×4 的布尔类型二维张量（矩阵）的稀疏表示形式。特别值得注意的是，稀疏张量数据默认按照第一个维度进行排序，因此对于总流出量，我们将其改写为，将数据的第二个维度和第一个维度交换，并按照交换后的第一个维度（原第二个维度）排序即可。

表达式 expr：传入一个 gl.Expression 对象，可以是经由其他计算节点处理后的表达式，也可以是一个叶子节点表达式，由变量和变量系数（默认为 1）表示。例如，可以传入叶子结点表达式来描述的表达式部分（红色标出）。

然后，构建约束对象 gl.constraint 并使用 model.add_constraint 方法添加到模型。约束需要指定约束公式的左表达式 expr、符号 sense（，或=）和右边的常数项 rhs。除此以外，还需指定约束下标 axes，列出所有在约束中使用的下标维度。例如，针对网络流约束的而指定 axes=[i]。在大规模排产模型中，约束往往至少有三维，即时刻 t - 工厂 p - 产品编号 i，此时需要指定 axes=[t, p, i]。为了进行并行化建模，还需指定用于进行数据切分的下标维度 parallel_axes，一般指定约束的第一个下标维度即可。

2. 高效求解：如何提升大规模商业数学规划模型的求解效率

那么，我们是不是只要对决策周期 T 进行简单的切分就行呢？例如，对于 28 天的排产周期，我们将其分成 4 份，首先排产 1-7 天，然后排产 8-14 天，以此类推。这样每次排产时求解器只需要求解原先 1/4 规模的模型，如果我们根据经验性的求解时间复杂度（N为决策变量数量）推断，分解后的求解时间会是直接求解的，不可谓不快。如果分解的块数更多，速度提升甚至会更加显著。我们将这种方式称之为“滚动时域求解”（Rolling Horizon，RH）。然而，这样做有一个致命的缺点，就是模型的决策将非常“短视”，每次决策都只关注于当前小的决策片段的局部最优性，而无法顾及全局。例如，若在排产周期的第 28 天对某物品有一个大的订单需求，且生产该物品需要 10 天，那么第 1-7、8-14、15-21 天的决策中都不会考虑这一需求，从而不进行足量的生产，而等到第 22-28 天开工生产已为时晚矣。分块的数量越多，短视的现象就越严重，从而严重影响到模型决策的最优性。

Grassland 在华为 78 周供需模拟场景下的端到端性能测试结果。蓝色部分为建模时间，红色部分为求解时间。

• 李锡涵，伦敦大学学院计算机系博士研究生，主要研究方向为学习优化，在 AAMAS、CIKM、COLING 发表学术论文。TensorFlow技术专家，开源技术手册《简明的 TensorFlow 2》（https://tf.wiki）作者。

• 韩雄威，研究生毕业于华中科技大学管理科学与工程专业，现任华为诺亚方舟实验室高级研究员，供应链高级算法科学家，主要研究方向为运筹学、学习优化等技术，致力将其运用于供应链领域的组合优化问题。

• 周之烁，本硕毕业于复旦大学大数据学院。硕士研究生期间师从江如俊副教授，研究方向为最优化算法、非凸优化。