运筹学（Operational
Research）是一门应用于管理有组织系统的科学，最早的朴素思想在中国的古文献中多有记载，比如耳熟能详的田忌赛马的故事。运筹的一般思想是：在各项资源条件优先的情况下，如何确定一个方案，使得预期目标最优；或者为了达到预期目标，确定资源消耗最小的方案。在二次世界大战之后，组织和企业的活动规模更大，信息系统化空前完备（想象一下水晶报表的诞生多么让人兴奋），加之各类数学算法模型层出不穷，研究如何做好决策的运筹学也有了极大的发展。
运筹学方向很多，比如线性规划、非线性规划、整数规划、目标规划、动态规划、排队论、对策论等。笔者偷个懒，找一些在整数规划体系下的例子，让大家感受一下在日常企业中这些方法的应用。
1. 一个简化问题
公司有 4 条生产线，每条生产线的月产量分别为 0.56, 3.11, 3.04,
2.11。近期因为经济不景气，需要将月产量总和控制在 5
以内，但出于总成本摊销的考虑，又要保证产出尽可能的大，那么哪几条产线需要被关闭。我们盲猜结果：3
和 4
需要被关闭。当然这个问题手指头可以掰过来，超过十个手指头怎么办？
我们先约定一些参数的数学表达：
每条生产线的月产量为 \(w_i\)
产线是否开启用 \(x_i\) 表示，要么是
0 要么是 1
那这个问题就可以用统一框架结构化表述该问题。
设置决策变量 \(x_{i} \in \{0, 1\}, i =
1,2,3,4\)
目标：最大化 \(\sum x_i w_i\)
约束条件：\(\sum x_i w_i \le
5\)
求解 \(x_i\)
我们尝试使用 R 包 ompr
来实现上述过程（Optimization Modeling Package）：
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 library( ompr) library( ompr.roi) library( ROI.plugin.glpk) max_capacity <- 5 n <- 4 weights <- c ( 0.56 , 3.11 , 3.04 , 2.11 ) m <- MIPModel( ) |> add_variable( x[ i] , i = 1 : n, type = "binary" ) |> set_objective( sum_over( weights[ i] * x[ i] , i = 1 : n) , "max" ) |> add_constraint( sum_over( weights[ i] * x[ i] , i = 1 : n) <= max_capacity) |> solve_model( with_ROI( solver = "glpk" ) ) |> get_solution( x[ i] ) > m variable i value 1 x 1 1 2 x 2 1 3 x 3 0 4 x 4 0
果然和我们猜想一样，3 和 4 需要被关闭。
P.S.
ompr
这个包非常优秀，将设置决策变量、最优化目标、约束条件，以及求解四个步骤，很清晰地通过管道操作来表达。具体细节各位看官可以参考官方文档。
2. 排班
假设某培训机构，周一到周日所需要的老师数为3782, 3718, 2580, 2912,
4112, 7862, 7007，设为 \(b\) ，且要保证老师每周有连续两天的休息，有两个问题需要回答：
机构需要准备多少老师
这些老师该如何排班
首先生成所有的排班可能 \(A_{ij}\) ，表示第\(i\) 天的第 \(j\) 个班型需要的人数，同时假设第 \(i\) 天的工作人数为 \(x_{i}\) ，那么满足
\[
\min z= \sum_{i} x_i \\
\]
\[
s.t. \sum_i\sum_j x_i A_{i,j} \le b_i
\]
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 b <- c ( 3782 , 3718 , 2580 , 2912 , 4112 , 7862 , 7007 ) A <- matrix( 1 , nrow = 7 , ncol = 7 ) for ( i in 1 : 7 ) { A[ i, i %% 7 + 1 ] <- 0 ; A[ i, ( i + 1 ) %% 7 + 1 ] <- 0 ; } m <- MIPModel( ) |> add_variable( x[ i] , i = 1 : 7 , type = "integer" , lb = 0 ) |> set_objective( sum_over( x[ i] , i = 1 : 7 ) , "min" ) |> add_constraint( sum_over( A[ i, j] * x[ i] , j = 1 : 7 ) >= b[ i] , i = 1 : 7 ) |> solve_model( with_ROI( solver = "glpk" ) ) |> get_solution( x[ i] ) > m variable i value 1 x 1 757 2 x 2 744 3 x 3 516 4 x 4 583 5 x 5 823 6 x 6 1573 7 x 7 1402
排班的班型和结果如下：
1
0
0
1
1
1
1
757
1
1
0
0
1
1
1
744
1
1
1
0
0
1
1
516
1
1
1
1
0
0
1
583
1
1
1
1
1
0
0
823
0
1
1
1
1
1
0
1573
0
0
1
1
1
1
1
1402
看一下差异有多大
1 2 3 4 > m$ value * A |> colSums( ) [ 1 ] 3785 3720 2580 2915 4115 7865 7010 > b[ 1 ] 3782 3718 2580 2912 4112 7862 7007
3. 装箱问题
有一堆物品长度分别为 \(w_i\) ，若干长度为 \(C\) 的箱子（\(w_i
\le C\) ）。问：最少需要几个箱子才能把所有物品全部装进箱子？
最坏的情况是一个箱子装一个物品，最好的情况是物品恰好可以完美的放入在箱子内，那么箱子数为\(\sum
w_i/C\) ，所以说取值应该在这两个值之间。抽象一下这个问题，我们希望最小化箱子数量，设
\(y_i = 1\) ，表示第 \(i\) 个箱子被使用\(x_{ij} = 1\) ，表示 第 \(j\) 个物品被放入了第 \(i\) 个箱子中
那么优化的目标是：
\[
\min z = \sum_{i=1}^{n} y_i
\] 约束条件是：
\[
\begin{aligned}
&\sum_{j=1}^n w_i x_{i j} \leq C y_i, \quad i=1,2, \ldots, n\\
&\sum_{i=1}^n x_{i j}=1, \quad j=1,2, \ldots, n\\
&y_i \in\{0,1\}, \quad i=1,2, \ldots, n\\
&x_{i j} \in\{0,1\}, \quad i, j=1,2, \ldots, n
\end{aligned}
\] 用一个例子来说明：
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 max_bins <- 10 bin_size <- 3 n <- 10 weights <- runif( n, max = bin_size) MIPModel( ) |> add_variable( y[ i] , i = 1 : max_bins, type = "binary" ) |> add_variable( x[ i, j] , i = 1 : max_bins, j = 1 : n, type = "binary" ) |> set_objective( sum_over( y[ i] , i = 1 : max_bins) , "min" ) |> add_constraint( sum_over( weights[ j] * x[ i, j] , j = 1 : n) <= y[ i] * bin_size, i = 1 : max_bins) |> add_constraint( sum_over( x[ i, j] , i = 1 : max_bins) == 1 , j = 1 : n) |> solve_model( with_ROI( solver = "glpk" , verbose = TRUE ) ) |> get_solution( x[ i, j] ) |> filter( value > 0 ) |> arrange( i)
4. 指派问题
指派的标准问题是，有 \(n\) 个人和
\(n\) 件事，已知第 \(i\) 个人做第 \(j\) 个事的费用为 \(C_{ij}\) ，要求确定人和事之间的关系，使得这
\(n\) 件事情的总费用最小。
数学模型可以写成：
\[
\min z=\sum\limits_{i=1}^n\sum\limits_{j=1}^n c_{ij}x_{ij}
\]
\[
\mathrm{s.t}
\begin{cases}
\sum\limits_{i=1}^n x_{ij}=1,\quad i=1,2,\cdots,n \\
\sum\limits_{j=1}^n x_{ij}=1,\quad j=1,2,\cdots,n \\
x_{ij} \in \{0,1\}, \quad i,j=1,2,\cdots,n \end{cases}
\]
模拟一个示例做演示：
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 set.seed( 1 ) cost <- rnorm( 16 , 4 , 1 ) |> round ( 2 ) |> matrix( 4 , 4 ) m <- MIPModel( ) |> add_variable( x[ i, j] , i = 1 : 4 , j = 1 : 4 , type = "binary" ) |> set_objective( sum_over( cost[ i, j] * x[ i, j] , i = 1 : 4 , j = 1 : 4 ) , "min" ) |> add_constraint( sum_over( x[ i, j] , i = 1 : 4 ) == 1 , j = 1 : 4 ) |> add_constraint( sum_over( x[ i, j] , j = 1 : 4 ) == 1 , i = 1 : 4 ) |> solve_model( with_ROI( solver = "glpk" ) ) |> get_solution( x[ i, j] ) > cost [ , 1 ] [ , 2 ] [ , 3 ] [ , 4 ] [ 1 , ] 3.37 4.33 4.58 3.38 [ 2 , ] 4.18 3.18 3.69 1.79 [ 3 , ] 3.16 4.49 5.51 5.12 [ 4 , ] 5.60 4.74 4.39 3.96 > z <- matrix( 0 , 4 , 4 ) > z[ cbind( m$ i, m$ j) ] <- m$ value> z [ , 1 ] [ , 2 ] [ , 3 ] [ , 4 ] [ 1 , ] 0 1 0 0 [ 2 , ] 0 0 0 1 [ 3 , ] 1 0 0 0 [ 4 , ] 0 0 1 0
5. 最优传输问题
最优传输比指派问题稍复杂一点点：假如有 6 个生产点和 8
个销售点，他们之间的存在距离成本矩阵。且约束条件有：
A1-A6 的产量是：60,55,51,43,41,52
B1-B8 的销量是：35,37,22,32,41,32,43,38
1 2 3 4 5 6 7 8 costs<- matrix( c ( 6 , 4 , 5 , 7 , 2 , 5 , 2 , 9 , 2 , 6 , 3 , 5 , 6 , 5 , 1 , 7 , 9 , 2 , 7 , 3 , 9 , 3 , 5 , 2 , 4 , 8 , 7 , 9 , 7 , 8 , 2 , 5 , 4 , 2 , 2 , 1 , 5 , 8 , 3 , 7 , 6 , 4 , 9 , 2 , 3 , 1 , 5 , 3 ) , nrow = 6 ) n_row <- nrow( costs) n_col <- ncol( costs) row.rhs <- c ( 60 , 55 , 51 , 43 , 41 , 52 ) col.rhs <- c ( 35 , 37 , 22 , 32 , 41 , 32 , 43 , 38 )
构建传输模型，假设生产点的产量为 \(P_i\) , 销售点的预期销量为 \(S_j\) ，生产点到销售点运费矩阵为 \(C_{ij}\) ，不考虑对于每次传输的限制，销售点的需求必须满足：
\[
\begin{gathered}
min \sum_i \sum_j C_{ij} X_{ij} \\
s.t. \\
\sum_{i}\sum_j X_{ij} \le P_i;i = 1,2...,M \\
\sum_{i}\sum_j X_{ij} = S_j;j = 1,2,...,N \\
\end{gathered}
\]
构建混合整数规划模型：
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 tm <- MIPModel( ) |> add_variable( x[ i, j] , i = 1 : n_row, j = 1 : n_col, type = "integer" , lb = 0 ) |> set_objective( sum_over( costs[ i, j] * x[ i, j] , i = 1 : n_row, j = 1 : n_col) , "min" ) |> add_constraint( sum_over( x[ i, j] , j = 1 : n_col) <= row.rhs[ i] , i = 1 : n_row) |> add_constraint( sum_over( x[ i, j] , i = 1 : n_row) == col.rhs[ j] , j = 1 : n_col) |> solve_model( with_ROI( solver = "glpk" ) ) |> get_solution( x[ i, j] ) m <- matrix( 0 , n_row, n_col) index <- cbind( tm$ i, tm$ j) m[ index] <- tm$ value > m [ , 1 ] [ , 2 ] [ , 3 ] [ , 4 ] [ , 5 ] [ , 6 ] [ , 7 ] [ , 8 ] [ 1 , ] 0 19 0 0 41 0 0 0 [ 2 , ] 1 0 0 32 0 0 0 0 [ 3 , ] 0 11 0 0 0 0 40 0 [ 4 , ] 0 0 0 0 0 5 0 38 [ 5 , ] 34 7 0 0 0 0 0 0 [ 6 , ] 0 0 22 0 0 27 3 0
6. 分货优化
最优传输问题是一个相对简化的模型，我们将这个问题更一般化。假设有 M
个大仓，N 个门店。第 \(i\)
个大仓的分货限额为 \(C_i=1,2,3,...,M\) ，需要将 \(C_i\) 分配给各个门店，第 \(i\) 个门店分得 \(X_{ij},
i=1,2,3,...,M,j=1,2,...,N\) ，使得
\[
min(OBJ)
\]
约束条件为
\[
\sum_{j}^{N} X_{i, j} \leq C_{i}, i=1,2, \ldots, M
\]
这里的 OBJ
是某种可以衡量分货策略好坏的函数，分货策略越好，该函数值越小。举个简单的例子，OBJ
可以是该分货策略的预期缺货量，那么显然，缺货量越小，说明分货策略越好。当然，在实际生产决策中，衡量分货策略好坏往往是从多个维度来考量，诸如现货率、补货次数、补货周期等核心优化指标均要纳入考虑。
6.1. 传输优先级
构建最简化的模型，仅仅考虑一个约束条件，传输的优先级（该情况即为最优传输模型）。\(j\) 个门店的库存为 0，完全达成门店要求。
此时抽象的公式为：
\[
min \sum_{i}\sum_{j} P_{ij} X_{ij}
\] 设门店的目标要求为 \(V_j\) ，约束条件为：
\[
\begin{gathered}
\sum_{j}^{M} X_{i, j} \leq C_{i}, i=1,2, \ldots, M \\
\sum_{i}^{N} X_{i, j} = V_{j}, j=1,2, \ldots, N
\end{gathered}
\]
构建模型，及结果是：
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 ## 大仓总量 wh_capacity <- c(100, 300) ## 门店目标库存水平 store_target <- c(100, 80, 120) ## 大仓优先级 priority <- matrix(c(0.9, 0.1, 0.9, 0.1, 0.05, 0.95), 2, 3) ## 构建模型 tm1 <- MIPModel() |> # 添加长度为 i * j 变量组 add_variable(x[i,j], i = 1:2, j = 1:3, type = "integer", lb = 0) |> ## 设置目标函数 set_objective(sum_over(priority[i,j] * x[i,j], i = 1:2, j = 1:3), "min") |> ## 增加约束 add_constraint(sum_over(x[i,j], j = 1:3) <= wh_capacity[i], i = 1:2) |> add_constraint(sum_over(x[i,j], i = 1:2) == store_target[j], j = 1:3) |> solve_model(with_ROI(solver = "glpk")) |> get_solution(x[i,j]) > tm1 variable i j value 1 x 1 1 0 4 x 2 1 100 2 x 1 2 0 5 x 2 2 80 3 x 1 3 100 6 x 2 3 20
6.2. 考虑实际库存
上面的抽象模型显然太过于理想化，事实上门店显然是有库存的。
假设目标库存车辆数是 \(Tar_j \ge
0\) ，当前库存数车辆数为 \(Inv_j \ge
0\) ，\(X_{ij} \ge 0\) 为从 \(i\) 位置到 \(j\) 站点的调度量，\(S_j \ge 0\) 为站点 \(j\)
取得调度量之后的未满足量。那么则有这样的限制条件：
\[
Inv_j + \sum_j^N X_{ij} + S_j \ge Tar_j
\]
显然我们希望未满足量最小：
\[
min \sum_j^N S_j
\] 也就是说我们需要的目标是两个同时优化：
\[
\min (\sum_i^M \sum_j^N P_{ij} X_{ij} + \sum_{j} S_{j})
\]
约束条件是：
\[
\begin{gathered}
\sum_{j} X_{ij} \leq C_{i}, i=1,2, \ldots, M \\
I n v_{j}+\sum_{i} X_{ij}+S_{j} \geq Tar_{j}, j=1,2, \ldots, N
\end{gathered}
\]
这个混合整数规划问题变成了更复杂的样子：
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 store_inv <- c ( 50 , 20 , 30 ) store_mip <- MIPModel( ) |> add_variable( x[ i, j] , i = 1 : 2 , j = 1 : 3 , type = "integer" , lb = 0 ) |> add_variable( b[ j] , j = 1 : 3 , type = 'integer' , lb = 0 ) |> set_objective( sum_over( priority[ i, j] * x[ i, j] , i = 1 : 2 , j = 1 : 3 ) + sum_over( 1000 * b[ j] , j = 1 : 3 ) , "min" ) |> add_constraint( sum_over( x[ i, j] , j = 1 : 3 ) <= wh_capacity[ i] , i = 1 : 2 ) |> add_constraint( sum_over( x[ i, j] , i = 1 : 2 ) + store_inv[ j] + b[ j] >= store_target[ j] , j = 1 : 3 ) |> solve_model( with_ROI( solver = "glpk" ) ) |> get_solution( x[ i, j] ) > store_mip variable i j value 1 x 1 1 0 4 x 2 1 50 2 x 1 2 0 5 x 2 2 60 3 x 1 3 90 6 x 2 3 0
6.3. 最小分货单元
因为分货的车每次只能装载 6 件商品，因此我们还有一个条件：分货量必须为
6 的倍数。可以考虑引入一个辅助变量 \(O_{ij}\) ，引入 mlc 之后约束条件增加
\[
\begin{gathered}
mlc \times O_{ij} = X_{ij}\
\end{gathered}
\] 模型又发生了一些改变：
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 mlc <- 6 store_mip <- MIPModel( ) |> add_variable( x[ i, j] , i = 1 : 2 , j = 1 : 3 , type = "integer" , lb = 0 ) |> add_variable( b[ j] , j = 1 : 3 , type = 'integer' , lb = 0 ) |> add_variable( o[ i, j] , i = 1 : 2 , j = 1 : 3 , type = 'integer' , lb = 0 ) |> set_objective( sum_over( priority[ i, j] * x[ i, j] , i = 1 : 2 , j = 1 : 3 ) + sum_over( b[ j] , j = 1 : 3 ) , "min" ) |> add_constraint( sum_over( x[ i, j] , j = 1 : 3 ) <= wh_capacity[ i] , i = 1 : 2 ) |> add_constraint( sum_over( x[ i, j] , i = 1 : 2 ) + store_inv[ j] + b[ j] >= store_target[ j] , j = 1 : 3 ) |> add_constraint( mlc * o[ i, j] == x[ i, j] , i = 1 : 2 , j = 1 : 3 ) |> solve_model( with_ROI( solver = "glpk" ) ) |> get_solution( x[ i, j] ) > store_mip variable i j value 1 x 1 1 0 4 x 2 1 54 2 x 1 2 0 5 x 2 2 60 3 x 1 3 90 6 x 2 3 0
当然，以上示例均为演示而构造，实际生产过程会远比这些示例复杂，有些明确需要使用其他方法来进行拟合或逼近。