1. 应用场景

列生成算法是不断添加求解变量的算法，可参考论文。上一节我们介绍了约束条件存在分块结构的情形下可以使用单纯型法的判断规则进行列生成，下面我们介绍另一种场景，set-partitionning问题的列生成算法。
具体来说，场景如下：有$n$个0-1变量$z_1...z_n$，每个$z_i$带着很多中间变量$x_{i,j}$用来进行约束，这是个变量很多，约束也很多的模型。我们首先对问题模型进行转化，将所有$z_i$的组合枚举出来，在枚举的过程中，就把约束条件都过了一遍，只有满足所有约束条件的组合才会保留下来。最后的枚举结果分别对应$y_1...y_m$。
举个VRP问题，每一列是一个可能路径（列生成指的就是生成闭环路径的规则），主问题的约束条件则是除了depot外的每个点访问过一次。
再举个航空排班的例子。假设机组有x个人需要执飞y个航班，每一列是一个可能的部分航班计划，而主问题的约束条件则是每个航班执飞过一次。

不难看出，$m$是$n$的指数形式，使用set-partitioning模型之后，变量的数量大的可怕。好消息是，原问题的约束条件比较紧，只有少数组合成立，可以通过一定的规则产生，这样m的数量比较小，问题就比较好解决。

2. 算法步骤

问题的列生成算法对应的模型是：
$\min {\Sigma }_j{y}_j$
s.t. $Ay\ge b^T$
$A$是一个矩阵，$b=[b_1,...,b_m]$和$y=[y_1,...,y_n]$是向量；并且有重要的一点：set-covering的变量$y_j$对应的列$a_j=[a_{1j},...,a_{mj}{]}^T$生成的规则，可以通过线性约束条件 $B\ast a_j\le D$得到。这使得我们可以用线性规划的方法来求解子问题~
注意由于列生成算法求解的是松弛线性规划问题，因此对整数规划模型需要结合branch and bound算法进行。

2.1 限制主问题

首先用启发式算法找出一部分$A$，比如说选出了$k$列。然后我们的线性规划问题就变成了：
min $y_1+y_2+...+y_k$
s.t. $a_{i,1}\ast y_1+a_{i,2}\ast y_2+...+a_{i,k}\ast y_k\ge b_i$，$i=1...n$
相比原来的模型，相当于把$y_{k+1}$~$y_m$强制限制为非基变量了，称为限制主问题（Restricted Master Problem，RMP）。上面的限制主问题求解完成后，我们想使用单纯型法进行基变量的转换，看看$y_{k+1}$~$y_m$中，是否有可以转入基变量的列。检验数$\sigma =c_N-\pi \ast a_j$，并且$\pi =c_B{B}^{-1}$可以由求解器给出。我们要找出非基变量中最小的负数$\sigma$，将其转入基变量。正如前面所述，我们这里使用线性规划来找这个新的列。

2.2 子问题

注意$c=[1,...,1]$，因此 min $c_N-\pi a_j$等价于
min $1-\pi \ast a$
s.t. $B\ast a\le D$ (列生成的规则)
称为子问题（sub problem）。如果目标函数最优值＜0，就将新生成的列y_k+1转入基变量，生成新的限制主问题进行求解。如此往复，直至子问题的目标函数≥0。

3. 例子

3.1 问题描述

来看网上流传很多的cutting stock problem的例子，问题如下：
我们有一些纸筒，每个纸筒16m。顾客需要裁剪的长度为：25个3m，20个6m，15个7m。要求在满足顾客需求的情况下，裁剪的纸筒数最小。

3.2 数学模型和初始可行解

我们可以用启发式算法找到一个upper bound的初始解：
5个方案1：5个3m
10个方案2：2个6m
8个方案3：2个7m
总计23筒。也就是说，我们用23筒是肯定可以满足要求的，这算问题的一个上界。下面我们探索一下其他的切筒方式，看能不能给出下界。
用set-covering对问题进行建模：设$P$是所有可行的裁剪方案的集合（其中前3个可以设置为我们前面的3个裁剪方案），里面方案的总数为$n$(我们并不需要确切的知道这个值是多少，只需要知道它很大)。令$a_{ij}$表示第$j$种方案里类别$i$的个数，$y_j$表示第$j$种方案的选择个数，原问题可以变为：
$\min y_1+...+y_n$
$a_{1,1}{y}_1+a_{1,2}{y}_2+...+a_{1,n}{y}_n\ge 25$
$a_{2,1}{y}_1+a_{2,2}{y}_2+...+a_{2,n}{y}_n\ge 20$
$a_{3,1}{y}_1+a_{3,2}{y}_2+...+a_{3,n}{y}_n\ge 15$
其中
$a_{1,1}=5$，$a_{2,1}=0$，$a_{3,1}=0$
$a_{1,2}=0$，$a_{2,2}=2$，$a_{3,2}=0$
$a_{1,3}=0$，$a_{2,3}=0$，$a_{3,3}=2$
对应最初的3种方案。
$a_{i,j}$满足条件：$3\ast a_{1,j}+6\ast a_{2,j}+7\ast a_{3,j}\le 16$（前面所说的列生成的规则）
注意我们在使用列生成算法的时候，将整数变量松弛为了连续变量。

3.3 第一轮

问题对应的初始解为$y_1=5$，$y_2=10$，$y_3=8$，$y_4=...y_n=0$
限制主问题如下：
$\min y_1+y_2+y_3$
$5y_1+0y_2+0y_3\ge 25$
$0y_1+2y_2+0y_3\ge 20$
$0y_1+0y_2+2y_3\ge 15$

下面是pymprog的sensitivity输出：

对应$\pi =(0.2,0.5,0.5)$，我们把新转入的变量记作$a_4=(a_{1,4},a_{2,4},a_{3,4})$，则子问题为：
min $1-0.2\ast a_{1,4}-0.5\ast a_{2,4}-0.5\ast a_{3,4}$
s.t. $3\ast a_{1,4}+6\ast a_{2,4}+7\ast a_{3,4}\le 16$ (列生成的规则)
$a_{i,j}\in Z$
求解结果为$a_4=(1,2,0)$，检验数$\sigma =1-\pi \ast a_4=-0.2\le 0$。此时问题目标为22.5。
将$y_4$进基，把这个结果添加到主问题当中去，开始第二轮迭代。

3.5 第二轮

添加了$y_4$后，限制主问题为：
$\min y_1+y_2+y_3+y_4$
$5y_1+0y_2+0y_3+1y_4\ge 25$
$0y_1+2y_2+0y_3+2y_4\ge 20$
$0y_1+0y_2+2y_3+0y_4\ge 15$

直观上来看，是添加了一个新列，所以称为列生成算法。上面的主问题求解结果为：

对应$\pi =(0.2,0.4,0.5)$，我们把新转入的变量记作$a_5=(a_{1,5},a_{2,5},a_{3,5})$，则子问题为：
min $1-0.2\ast a_{1,5}-0.4\ast a_{2,5}-0.5\ast a_{3,5}$
s.t. $3\ast a_{1,5}+6\ast a_{2,5}+7\ast a_{3,5}\le 16$ (列生成的规则)
$a_{i,j}\in Z$
求解结果为$a_5=(3,0,1)$，检验数$\sigma =1-\pi \ast a_5=-0.1<0$。此时问题目标为20.5。
将$y_5$进基，把这个结果添加到主问题当中去，开始第三轮迭代。

3.6 第三轮

添加了$y_5$后，限制主问题为：
$\min y_1+y_2+y_3+y_4+y_5$
$5y_1+0y_2+0y_3+1y_4+3y_5\ge 25$
$0y_1+2y_2+0y_3+2y_4+0y_5\ge 20$
$0y_1+0y_2+2y_3+0y_4+1y_5\ge 15$

上面的主问题求解结果为：

对应$\pi =(0.166667,0.416667,0.5)$，我们把新转入的变量记作$a_6=(a_{1,6},a_{2,6},a_{3,6})$，则子问题为：
min $1-0.166667\ast a_{1,6}-0.416667\ast a_{2,6}-0.5\ast a_{3,6}$
s.t. $3\ast a_{1,6}+6\ast a_{2,6}+7\ast a_{3,6}\le 16$ (列生成的规则)
$a_{i,j}\in Z$
求解结果为$a_6=(1,1,1)$，检验数$\sigma =1-\pi \ast a_6=-0.083334<0$。此时问题目标为20。
将$y_6$进基，把这个结果添加到主问题当中去，开始第三轮迭代。

3.7 第四轮

添加了$y_6$后，限制主问题为：
$\min y_1+y_2+y_3+y_4+y_5+y_6$
$5y_1+0y_2+0y_3+1y_4+3y_5+y_6\ge 25$
$0y_1+2y_2+0y_3+2y_4+0y_5+y_6\ge 20$
$0y_1+0y_2+2y_3+0y_4+1y_5+y_6\ge 15$

上面的主问题求解结果为：

对应$\pi =(0.2,0.4,0.4)$，我们把新转入的变量记作$a_7=(a_{1,7},a_{2,7},a_{3,7})$，则子问题为：
min $1-0.2\ast a_{1,7}-0.4\ast a_{2,7}-0.4\ast a_{3,7}$
s.t. $3\ast a_{1,7}+6\ast a_{2,7}+7\ast a_{3,7}\le 16$ (列生成的规则)
$a_{i,j}\in Z$
求解结果为$a_7=(3,1,0)$，检验数$\sigma =1-\pi \ast a_6=0$，结束迭代，此时问题目标为19。

3.8 结果

使用下列代码求出最优解：

from pymprog import *
p = model('example')
y = p.var('y',6,int) # variables
p.minimize(sum(y), 'master')
r1 = 5*y[0]+ y[3]+3*y[4]+y[5]>= 25 
r2 = 2*y[1]+2*y[3]+y[5] >= 20 
r3 = 2*y[2] +y[4]+y[5]>= 15 
p.solve()
for i in range(6):
    print(y[i].primal)

求解结果为：$y_1=1,y_4=3,y_5=1,y_6=14$，总计需要19个纸筒。
也可以将最后一轮的$y$取整得到一个可行解：$y_1=2,y_4=3,y_6=15$，需要20个纸筒。
这里给出第四轮的限制主问题和子问题的代码供参考：

from pymprog import *
p = model('master problem')
y = p.var('y',6) # variables
p.minimize(sum(y), 'master')
r1 = 5*y[0]+ y[3]+3*y[4]+y[5]>= 25 
r2 = 2*y[1]+2*y[3]+y[5] >= 20 
r3 = 2*y[2] +y[4]+y[5]>= 15 
p.solve()
p.sensitivity()

from pymprog import *
s = model('sub problem')
a = s.var('a',3) # variables
for i in range(3):
    a[i].kind = int
s.minimize(1-(0.2*a[0] + 0.4*a[1] +0.4*a[2]), 'sub')
3*a[0] +6*a[1] +7*a[2]<= 16 # mountain bike limit
s.solve()
print(s.vobj())
for i in range(3):
    print(a[i].primal)

然后下面给出直接求最优解的代码。pymprog用了一天都没有求出来，因此改用cplex进行求解：

import cplex
from cplex.exceptions import CplexError
my_obj = [0.0,0.0,0.0,1.0]*23
my_ub = [5.0,2.0,2.0,1.0]*23
my_lb = [0.0,0.0,0.0,0.0]*23
my_ctype = "IIII"*23
my_colnames = []
my_rhs = [25,20,15]+[0,16]*23
my_sense = "GGG"+"LL"*23
for i in range(23):
    my_colnames = my_colnames + ["x0,"+str(i),"x1,"+str(i),"x2,"+str(i),"y"+str(i)]

def populatebyrow(prob):
    prob.objective.set_sense(prob.objective.sense.minimize)
    prob.variables.add(obj=my_obj, lb=my_lb, ub=my_ub, types=my_ctype,names=my_colnames)
    row0 = [1.0,0.0,0.0,0.0]*23
    row1 = [0.0,1.0,0.0,0.0]*23
    row2 = [0.0,0.0,1.0,0.0]*23    
    rows = [[my_colnames, row0],
            [my_colnames, row1],
            [my_colnames, row2]]
    for i in range(23):
        row3 = [0.0,0.0,0.0,0.0]*23
        row3[i*4:i*4+4]=[1.0,1.0,1.0,-100.0]
        row4 = [0.0,0.0,0.0,0.0]*23
        row4[i*4:i*4+4]=[3.0,6.0,7.0,0.0]
        rows = rows+[[my_colnames, row3],
            [my_colnames, row4]]
    prob.linear_constraints.add(lin_expr=rows, senses=my_sense,rhs=my_rhs)

try:
    my_prob = cplex.Cplex()
    handle = populatebyrow(my_prob)
    my_prob.solve()
    
except CplexError as exc:
    print(exc)

print("Solution status = ", my_prob.solution.status[my_prob.solution.get_status()])
print("Solution value  = ", my_prob.solution.get_objective_value())
x = my_prob.solution.get_values()
for i in range(23):
    if x[i*4+3]>0:
        print(x[i*4:i*4+3])

结果如下：

CPXPARAM_Read_DataCheck                          1
Tried aggregator 1 time.
MIP Presolve modified 23 coefficients.
Reduced MIP has 49 rows, 92 columns, and 230 nonzeros.
Reduced MIP has 23 binaries, 69 generals, 0 SOSs, and 0 indicators.
Presolve time = 0.00 sec. (0.12 ticks)
Found incumbent of value 22.000000 after 0.00 sec. (0.48 ticks)
Probing time = 0.00 sec. (0.01 ticks)
Tried aggregator 1 time.
Reduced MIP has 49 rows, 92 columns, and 230 nonzeros.
Reduced MIP has 23 binaries, 69 generals, 0 SOSs, and 0 indicators.
Presolve time = 0.00 sec. (0.21 ticks)
Probing time = 0.00 sec. (0.01 ticks)
MIP emphasis: balance optimality and feasibility.
MIP search method: dynamic search.
Parallel mode: deterministic, using up to 4 threads.
Root relaxation solution time = 0.00 sec. (0.13 ticks)

        Nodes                                         Cuts/
   Node  Left     Objective  IInf  Best Integer    Best Bound    ItCnt     Gap

*     0+    0                           22.0000        0.0000           100.00%
      0     0        6.6667    38       22.0000        6.6667       30   69.70%
      0     0       10.2791    39       22.0000      Cuts: 62      143   53.28%
      0     0       18.0972    40       22.0000      Cuts: 80      213   17.74%
*     0+    0                           21.0000       18.0972            13.82%
      0     0       18.4688    35       21.0000      Cuts: 38      305   12.05%
      0     0       18.7500    19       21.0000      Cuts: 41      395   10.71%
      0     0       18.7500    12       21.0000      Cuts: 18      444   10.71%
      0     0       18.8667    10       21.0000      Cuts: 12      496   10.16%
*     0+    0                           19.0000       18.8667             0.70%
      0     0        cutoff             19.0000       18.8667      496    0.70%
Elapsed time = 0.18 sec. (7.56 ticks, tree = 0.01 MB, solutions = 3)

Implied bound cuts applied:  22
Mixed integer rounding cuts applied:  35
Zero-half cuts applied:  5
Gomory fractional cuts applied:  7

Root node processing (before b&c):
  Real time             =    0.19 sec. (7.57 ticks)
Parallel b&c, 4 threads:
  Real time             =    0.00 sec. (0.00 ticks)
  Sync time (average)   =    0.00 sec.
  Wait time (average)   =    0.00 sec.
                          ------------
Total (root+branch&cut) =    0.19 sec. (7.57 ticks)
Solution status =  MIP_optimal
Solution value  =  19.0
[1.0, 1.0, 1.0]
[1.0, 1.0, 1.0]
[1.0, 2.0, 0.0]
[1.0, 1.0, 1.0]
[1.0, 1.0, 1.0]
[1.0, 1.0, 1.0]
[1.0, 1.0, 1.0]
[1.0, 1.0, 1.0]
[1.0, 1.0, 1.0]
[1.0, 1.0, 1.0]
[1.0, 2.0, 0.0]
[1.0, 1.0, 1.0]
[1.0, 2.0, 0.0]
[1.0, 1.0, 1.0]
[5.0, 0.0, 0.0]
[1.0, 1.0, 1.0]
[1.0, 1.0, 1.0]
[3.0, 0.0, 1.0]
[1.0, 1.0, 1.0]

最终是19个。

前言

之前碰到了扔鸡蛋问题（给你2个鸡蛋，在100层楼上扔，要求想出一个好的策略，去测出哪一层楼开始鸡蛋就会碎掉），一直摸不着头脑。后来才知道可以使用“动态规划”这种思想（或者叫算法范式）去解决这个问题。

但是看了一些鸡蛋问题和动态规划的文章，依然只是流于形式，并不能理解其中的精髓。想想或许是鸡蛋问题对我现在而言难了一些，所以只好找了一些其它的动态规划问题，从简单入手，先去理解“动态规划”的思想精髓，再反过来去思考“鸡蛋问题”。

其中一个经典问题是01背包问题，这是我之前一直想搞懂的一个问题。看了一篇文章，就大致上理解了这个问题的解决办法和思路。所以我就不班门弄斧了，直接推荐这篇文章：动态规划之01背包问题（最易理解的讲解）。

比起一开始就搬概念和公式，这篇文章用了一个“填表格”的方式去让读者理解背包问题的解决过程，这对初学者更友好。推荐大家也试着去画出那个表格，然后自己填一填，不需要编程，只需要一张纸和笔和你的脑子，就能按照文中的思路把表格填完，然后就解决了背包问题。算法编程实际只是对你脑子的这种思维过程的一个模拟。

我想，之所以这种所谓“思想”、“范式”之类的东西难以理解，也许是因为它经过了双重的“抽象”吧。比如，以下是动态规划思想的描述：

动态规划算法与分治法类似，其基本思想也是将待求解问题分解成若干个子问题，先求解子问题，然后从这些子问题的解得到原问题的解。与分治法不同的是，适合于用动态规划求解的问题，经分解得到子问题往往不是互相独立的。若用分治法来解这类问题，则分解得到的子问题数目太多，有些子问题被重复计算了很多次。如果我们能够保存已解决的子问题的答案，而在需要时再找出已求得的答案，这样就可以避免大量的重复计算，节省时间。我们可以用一个表来记录所有已解的子问题的答案。不管该子问题以后是否被用到，只要它被计算过，就将其结果填入表中。这就是动态规划法的基本思路。具体的动态规划算法多种多样，但它们具有相同的填表格式…

有时候，真的想吐槽这些东西就是“正确的废话”。你不能说这些描述有什么错，但是真遇到问题的时候，光看这些描述你又完全没办法想出解决问题的方法。

解决某个具体问题，我们会把问题进行抽象，例如把搜索问题抽象成树问题。而这些所谓的“思想”，为了归纳这一类问题的共性，又成了对这类问题的抽象。即，成了“抽象之抽象”，双重抽象。

普适的代价即抽象

不知道别人如何，反正我对于“抽象”这类反人类大脑结构的东西实在无法理解。只能重新“化抽象为形象”，才能去“识”得它。就算是概念、思想这类看不见摸不着的东西，还是得转换成我脑海一个具象化的东西才行。

理解了01背包问题之后，我开始找到第2个动态规划的问题，即数塔问题。企图通过解决这些问题，一窥“动态规划”的冰山一角。

下面进入正题。

问题

如图，有一个数塔，从顶部出发，在每一个节点可以选择向左或者向右走，一直走到底层。要求找出一条路径，使得所走路径上的节点数字之和最大。

（图片来自百度）

分析

1.穷举

1）思路来源

这是在没有去看正确的动态规划方法之前，我自己想出的思路。

虽然不是“动态规划”，但是这个走歪的方向，也许更能让我们理解究竟怎样才算是“动态规划”。所以这个思路也是有价值的，可以说一说。

我一开始的理解，“动态规划”就是“使用已经计算好的信息，避免重复的计算”。这个想法的来源来自于这篇文章：

如何理解动态规划？

文章用了一张图来解释什么是动态规划：

我用当时的个人理解来复述一下：要计算A -> Q -> Z这条线路的最优选择，实际上并不需要我们把 A -> H -> Q -> Z、 A -> I -> Q -> Z 这两条路径都计算出来。 我们在计算A -> Q 的时候， 实际上已经把A -> H ->Q 、A -> I -> Q 这部分计算过了，因此已经得出了A -> Q的最优线路。计算A -> Q -> Z的时候，只需要利用这个结果即可。A -> Q是A -> Q -> Z的一个部分（注意：这个说法不完全正确，所以才会导致我得出了一个不完美的思路），因此能够“使用已经计算好的信息，避免重复的计算”。

当时我的理解不能说十分准确，认为动态规划就是：一个子问题如果属于一个大问题的一部分，就可以利用子问题的结果信息去计算大问题，避免重复计算。这个说法不准确，我后面会纠正这种想法。

但此时我们就顺着这种错误想法，看看这种想法会产生怎样的解题思路。

2）思维过程

毫无疑问，我们可以把问题图中的数塔看成一个大数塔，它其实可以看作包含了很多小数塔。比如只有两层的 9、12、15：

或者三层的9、12、15、10、6、8：

总之一个大数塔是可以拆分成小数塔的。那么，为了找到一个大数塔哪条路最优，是不是可以对小数塔找到最优？然后每一层每一层我们都最优，最后就能得到大数塔的最优了呢？

咋听起来很诱人，似乎找到了很厉害很巧妙的办法。但可惜，在这里行不通。这实际上就是“贪心算法”（em…又一种算法思想范式…）的本质：只要每一步我们都找最优，找到最后我们就能找到最后的最优。

还是引用推荐文章的一张图：

去计算A -> Z的最优线路，我们先一步一步来。先算A到“H or I”的最优线路，然后往下算“H or I”到“P or Q”的最优线路，最后算“P or Q”到Z的最优线路。算到最后，我们就得到了A到Z的最优线路。

这种情况下，我们才能用“贪心法”。看上去似乎跟“动态规划”的图不太一样，我把这张图改改，把“H or I”、“P or Q”这两个节点展开：

可以看到，任意一层的任意节点，都与下一层的所有节点有连接。第一层的A，与第二层的H、I都是连接的。第二层的H、I，与三层的P、Q都是连接的。第三层的P、Q，与第四层的Z都是连接的。

大概只有问题的推理形式具有这样的图特性，我们才能使用“贪心”法去解决吧。

我来再来看看我们的数塔，如果我把数塔横着倒放，变成这样：

看起来是不是跟“贪心法”的结构图不一样？是不是更像“动态规划”的那张结构图：

这也就说明了，为什么数塔问题可以用“动态规划”去解决。

3）解决办法（当然，最后想歪了，得出了穷举法。。。）

既然没办法通过一步步计算最优，来得出最后的最优，那我就妥协一些。每一次计算某一级子树塔的时候，我就用低一层级的子数塔的全部路线，然后再跟这一层的全部节点作一个连接，得出新的路线。

这样一步步不断计算出路线、以及这些路线的节点数值之和，最后把所有路线全部计算出来，作一个排序，就能知道哪条路线节点数值之和最大，不是就能得出问题答案了么？

想出这个办法的时候，我并不能确定这到底是不是“动态规划”，毕竟这时候对“动态规划”的理解还太浅了。但是又感觉这似乎不是“穷举”，因为我并没有把从根节点到末尾节点的路线全部列举一遍（至少从程序上我不会），我是“一步一步”计算的路线，并且我也把大数塔拆成来小数塔解决，相当于“把大问题拆分成相似的子问题”。从这些特征来看，这似乎也不算是“太笨的方法”，怎么会是“穷举”呢？(￣▽￣)"

再根据这个思路写写程序，发现我居然可以用递归！？这么一来似乎就更自信了。因为之前曾经看过，“递归”是“动态规划”过程中使用的一种手段（当然，这句话也不准确。。。后面也会纠正它），这些特征表面，我使用的方法就算是“动态规划”了吧？

怀着将信将疑的态度，试着去查找了别人的文章，发现想了那么多层，最后自己还是没绕出“穷举”这个坑。/(ㄒoㄒ)/~~

4）分析纠正

前面说了几点“动态规划”特征的个人理解，很遗憾的是，那几点理解都不是严格地准确。这里我们来分析、纠正一下它们：

<1> 一个子问题如果属于一个大问题的一部分，就可以利用子问题的结果信息去计算大问题，避免重复计算

其实这并算是“动态规划”独有的特征。算法优化的本质就是为了“避免重复计算”。那么问题来了：都是在“避免重复计算”，不同的方法又有什么不同？它们避免重复计算的部分不一样吗？

而且，不是只有动态规划问题才能拆成子问题，例如贪心类型问题，把一段长的路拆成一步步的小路，把整体最优拆成计算局部最优，不也是可以大问题拆成小问题么？

<2>“递归”是“动态规划”过程中使用的一种手段

同理，并不是只有动态规划才会用到递归。“递归”说来，实际上是属于编程层面的东西，“函数自己调用自己”，是一种代码文字的书写、组织、表达方式，一种手段。而动态规划是一种思想，是可以脱离编程、程序而存在的，你可以在自己脑海中去运行这种思想。

2.动态规划

1）理论分析

正经的动态规划方法我参考学习了这篇文章：几个经典的动态规划问题

事实上，动态规划方法最重要的，是找到一个当前状态与前一状态对应的关系。因为在动态规划结构的问题中，当前状态会依赖于上一个状态（或认为是，当前的“可选”与之前的“选择”有关）。

如果拿数塔来说的话，就是我这一层能够选择哪些岔路，跟我之前走过、选择的路是有关系的。比如下图，我走了第一层的“9”，然后走了第二层的“15”，到了第三层我就只能走“6”、“8”这两种选择了，而不可能走“10”这个节点。这就说明我第三层的可选项跟我第一、二层的选择是有关系的。

如果是贪心结构的问题，比如下图，那当前的选择，就跟我之前的选择没有什么关系。比如我走到第二层，将要考虑如何选择第三层的路，不管我现在是处于H还是I，我的可选项都是一样（都是P、Q），我可以走到第三层的任意节点，因为H、I与第三层所有节点都是通的。每一步都随你瞎走，反正最后给你的选项都是一样的，至于路线是不是最优，那就是另说了。

如果把这种“可选项”理解成一种状态的话，那么动态规划结构的问题，每一层都是有很多种状态的，且不同的状态取决于上一层是的哪种状态。而贪心结构的问题，每一层的状态都只有一种。

明白了动态规划问题存在着“状态”，并且状态与上个状态是存在着某种关系之后，我们去找出这种关系，将这种关系总结、归纳成“状态转换方程”，然后就是去编程实现了。动态规划问题的基本思路就是如此（当然这还是“正确的废话”，就算告诉你要去找状态转换方程，你也未必就能找出๑乛◡乛๑）

直接说吧，数塔问题的状态转换方程是这个：

dp(i, j） = Max( dp(i+1, j), dp(i+1, j+1 ) )+ f [ i ][ j ]

还需要把它翻译成人话：

为了便于说明，咱们先约定一下节点的表示方法吧：[ i, j ] ，它指的是第 i 层的第 j 个节点。比如下图[ 1, 1 ] 节点就是第1层、第1个节点，即“9”。

dp(i, j)函数的值，是[ i, j ] 这个节点，到最底层之后，所有可能路线的节点数值之和的最大值。还拿上图来说，比如dp(5, 1)其实就是[5, 1]这个节点（第5层、第1个）到最底层（第5层，也是它自己所在这一层）所有可能路线只有一条，即它自身“19”，所以dp[5, 1] = 19。

又比如dp(4, 2)，即[4, 2]这个节点，它到最底层（还是第5层）的所有路线有两条：7 -> 18、10 -> 18。数值和分别是25、28，28更大，所以最大值为28。因此dp(4, 2) = 28。

所有情况依次类推。

然后我们来解释一下这个状态转移方程式在说什么吧，dp(i, j） = Max( dp(i+1, j), dp(i+1, j+1 ) )+ f [ i ][ j ]，意思是dp(i , j)这个的值，等于dp(i+1, j)、dp(i+1, j+1 ) 两者之间取大的，然后再加上 f [ i ][ j ]这个。

f [ i ][ j ]的含义是节点[i, j]的数值，比如上图[1, 1]节点，值是9，所以f[1][1 ]= 9。

再翻译得更通俗，就是你想要知道[i, j]这个节点到最底层的最大数值和，首先要知道[i+1, j]和[i+1, j+1]这两节点到最底层的最大数值和。知道之后，我们取大的，然后加上[i, j]节点的值，即f[i][j]，就是我们要求的结果。

我们要求第1层、第1个节点到最底层的最大数值和，本质上就是要求dp(1, 1)。

2）实际操作

即使把公式翻译得通俗了，可能依然很绕口。我们就直接拿题目来说，一步步地来说明我们每一步怎么做。

想要知道根节点“9”到最底层的最大数值和，首先得知道第二层的“12”和“15”到最底层的最大数值和。我假设它们分别是M和N，假如“15”这条线路的数值和更大，即Max(M, N) = N，那么“9”到最底层的最大数值和，就是N + 9, 即下一层级的数值和 加上它自己。

然后我们反推一步，按照第一步的类似原理，我们要知道“12”和“15”到最底层的最大数值和，就得知道第三层的“10”、”6“、”8“到最底层的最大数值和。

其中，连接“12”这个节点的是“10”、“6”。我假设这两条支线的数值和最大分别是O、P，假如”10“这条线路的更大，即Max(O、P) = O，那么“12”这条线到最底层的最大数值和就是 O + 12。

其中，连接“15”这个节点的是“6”、“8”。我假设这两条支线的数值和最大分别是P、Q（注意到，P就是上段话里的那个P，因为节点都是“6”，所以这个值是一样的），假如”6“这条线路的更大，即Max(P，Q) = P，那么“15”这条线到最底层的最大数值和就是 P + 12。

然后我们继续，往下一层推理、不断推理…

最后推理到了最底层，就没法继续往下推了，已经到了一个结束点。此时节点到最底层的数值和最大值，就是它自身。比如“19”，它到最底层线路的数值和最大就是它自己——19了。

现在我们明白了，要想知道第一层的信息，就必须知道第二层。 要想知道第二层的信息，又必须知道第三层…依次类推，只有最底层的信息，我们不需要由别的层来推导，而是可以直接得出。因此，最底层就是一个计算的起始点。

若从实现来说，确实也需要用到递归。最底层既是计算的起始点，同时又是递归的终止点、回溯点。

实现

说了那么多，我们就用程序来实现、验证吧。

1）数塔生成器

先来实现一个数塔生成器吧：

function createNumberTower (level) {
    let tower = new Array();
    for (let i = 1; i < level + 1 ; i++) {
        let arr = new Array(i)
        for (let j = 0; j < arr.length; j++) {
            arr[j] = parseInt(Math.random() * 20)
        }
        tower.push(arr)
    }
    return tower;
}

实际上就是生成一个二维数组，只是第二维的长度变化是有规律的，依次是1、2、3、4…这样自增。例如题目图里的数塔用数组表示就是[[9], [12, 15], [10, 6, 8], [2, 18, 9, 5], [19, 7, 10, 4, 16]]

参数是level，即要生成的数塔的层数。

这里我把每个节点设置为0~20的随机数，要想生成别的数，可以自行修改。

arr[j] = parseInt(Math.random() * 20)

2) 穷举法

虽然这个方法的效率低，但是既然已经思考出来了，那还是把它也一块实现了吧。穷举法的结果一般比较准确，正好可以用来验证我们的动态规划法的结果是否正确。

function findPath(tower) {
    let level = tower.length - 1;
    let prevPaths = tower.length === 2 ? [{path: [0], value: tower[0][0]}] : findPath(tower.slice(0, level))
    let newPaths = new Array();
    for (let path of prevPaths) {
        let lastNode = path.path[path.path.length - 1];
        newPaths.push({ 
            path: Array.from([...path.path,lastNode]),
            value: path.value + tower[level][lastNode]
        })
        newPaths.push({ 
            path: Array.from([...path.path,lastNode + 1]),
            value: path.value + tower[level][lastNode + 1]
        })

    }
    return newPaths;
}

结果将是一个数组，数组的每个元素是一条路径。里面包含两个字段信息：path、value。

path也是一个数组，用来表示路径每一层的节点索引。那题目图的数塔举例子，例如某条路径的path值是[0, 1, 2, 3, 4]，那就说明第一层我走的是第1个元素（个数 = 索引 + 1嘛），第二层走的也是第2个元素，第三层走的是第3个元素，第四层走的是第4个元素，第五层走的是第5个元素。即9 -> 15 -> 8 -> 5 -> 16这一条路线。

value就是这一条路线所有节点的数值之和了。穷举了所用路线之后，按照value由大到小排序，就能知道哪一条路线数值之和最大了。

3）动态规划法

出乎意料，动态规划法的代码极其简单优雅，不过数行。

let fullTower = createNumberTower(30); // 生成数塔
function dp(i, j) {
    if (i === fullTower.length - 1) { // 起始点
        return fullTower[i][j];
    } else {
        return Math.max(dp(i+1, j), dp(i+1, j+1)) + fullTower[i][j];
    }
}

注意：为了编程方便，i，j 我用的是数组索引，即初始值是0。在这里dp(0 , 0)表示的才是对第1层第1个节点求解，而不是上文分析说明里的dp(1, 1)了。

4）验证

我们来验证一下，先把原题目求解一下。这里就不需要我们的数塔生成器了，直接把题目数塔写成一个二维数组

let fullTower = [[9], [12, 15], [10, 6, 8], [2, 18, 9, 5], [19, 7, 10, 4, 16]]

然后是穷举法：

findPath(fullTower).sort((prev, next) => next.value - prev.value)
输出：
 [{ path: [ 0, 0, 0, 1, 2 ], value: 59 },
  { path: [ 0, 1, 1, 1, 2 ], value: 58 },
  { path: [ 0, 0, 0, 1, 1 ], value: 56 },
  { path: [ 0, 0, 1, 1, 2 ], value: 55 },
  { path: [ 0, 1, 1, 1, 1 ], value: 55 },
  { path: [ 0, 1, 2, 3, 4 ], value: 53 },
  { path: [ 0, 0, 0, 0, 0 ], value: 52 },
  { path: [ 0, 0, 1, 1, 1 ], value: 52 },
  { path: [ 0, 1, 2, 2, 2 ], value: 51 },
  { path: [ 0, 1, 1, 2, 2 ], value: 49 },
  { path: [ 0, 0, 1, 2, 2 ], value: 46 },
  { path: [ 0, 1, 2, 2, 3 ], value: 45 },
  { path: [ 0, 1, 1, 2, 3 ], value: 43 },
  { path: [ 0, 1, 2, 3, 3 ], value: 41 },
  { path: [ 0, 0, 0, 0, 1 ], value: 40 },
  { path: [ 0, 0, 1, 2, 3 ], value: 40 }]

可见，数值和最大的路径是[ 0, 0, 0, 1, 2 ]这种路线，也就是数塔图的9 - > 12 -> 10 -> 18 -> 10这种走法。该路线的和为59。

我们用动态规划来验证一下：

dp(0, 0)
输出：59

对比一下，发现与穷举法结果一致。证明我们的动态规划函数也是对的。为了简单性，这个dp函数里，我就不像穷举法一样，去列出具体的路径节点了。

我们再用数塔函数，随意生成一个数塔来试试。来生成一个层数更多的数塔吧，比如10层的：

let fullTower = createNumberTower(10);
输出：[
  [ 16 ],
  [ 13, 0 ],
  [ 2, 14, 1 ],
  [ 10, 5, 17, 15 ],     
  [ 14, 9, 10, 14, 16 ], 
  [ 1, 7, 2, 18, 5, 19 ],
  [
    11,  0, 8, 5,        
     4, 11, 5
  ],
  [
    13, 13, 11,  3,
     2,  3,  9, 16
  ],
  [
     3, 18,  3, 17, 13,
    11, 18, 10, 18
  ],
  [
    17,  2, 9,  7, 8,
     2, 13, 8, 11, 5
  ],
  [
     4, 5, 15,  0,  6,
     9, 0, 12, 14, 19,
    19
  ],
  [
     7,  3, 17,  3, 3,
    17, 13, 10, 11, 6,
     5,  2
  ],
  [
    14,  0, 5, 2, 9, 17,
    10, 18, 6, 4, 6,  4,
    19
  ],
  [
    18,  9, 5, 12, 15,  8,
    13,  0, 3,  2, 14, 17,
    11, 15
  ],
  [
     1, 6,  4, 18,  8, 3,
     0, 1, 12, 12, 15, 6,
    11, 7,  4
  ],
  [
     7,  5,  4,  0,  1, 19,
    10,  6, 18,  5, 16, 12,
     2, 13, 10, 11
  ],
  [
     1, 19, 11,  0, 15, 4, 15,
    18,  0,  9, 14, 18, 4, 11,
     1,  5,  7
  ],
  [
    10, 8,  0,  8,  5, 11,  2,
     1, 5, 16, 19, 16, 15, 14,
    14, 4, 16, 15
  ],
  [
    3, 19,  2,  4, 10,  2, 16,
    1, 12, 10,  6,  5, 16, 15,
    8,  7, 18, 10, 15
  ],
  [
    2,  1, 10,  4, 12, 19, 11,
    7, 16, 13,  0, 18, 14, 17,
    0,  1, 13, 15,  5,  5
  ],
  [
    16,  3,  1,  0,  1,  2,  5,
     4,  8, 10, 17,  2, 18, 19,
     6, 12,  8, 13, 10,  6, 16
  ],
  [
     4, 15,  9,  6, 10, 17, 2,
    11,  7, 17, 16,  2,  5, 2,
    10, 15,  9, 15,  0, 10, 7,
     8
  ],
  [
    18, 18, 11, 11,  4,  1, 4, 13,
    18, 13,  7, 18, 10, 15, 7,  5,
    12, 18, 10, 15, 12,  0, 9
  ],
  [
     8,  8, 11, 18,  6,  4, 17, 15,
     7,  2,  7, 12, 12, 11,  9,  7,
    17, 15,  5,  3, 12, 10,  5,  7
  ],
  [
    14,  2,  5, 16,  1,  6,  7,  6,
    10, 15,  6,  1,  6,  3, 16,  4,
     8,  0, 16,  5, 11, 10,  9, 16,
     1
  ],
  [
    13,  0, 13,  3, 19, 11, 3, 9,
     5,  5, 11, 11, 15,  7, 9, 3,
    11,  2,  0,  2, 11, 19, 5, 1,
    18, 12
  ],
  [
    10, 10, 16, 12, 14,  3, 19, 11,
     9,  0, 10,  6,  7,  0,  6, 14,
     5,  1, 11, 10,  6, 13,  5,  5,
     3, 14,  5
  ],
  [
     3,  5,  4, 14, 13, 15,  1,  8,
    18, 14,  5, 16, 11,  5, 10,  6,
    15, 10, 10,  6
  ],
  [
    10, 18,  7, 17,  1, 10, 17,  9, 3,
    17, 16, 10,  3,  1, 17,  1, 12, 6,
     9, 19, 13, 18, 14,  0,  5,  5, 0,
    19, 17
  ],
  [
     0, 17, 18, 16,  2,  5,  9, 13, 15,
     0, 18, 13,  8,  0,  0, 15, 16,  4,
    12, 10, 12,  8, 15, 16,  2,  9,  1,
     4,  7, 15
  ]
]

穷举法的结果是：

findPath(fullTower).sort((prev, next) => next.value - prev.value)
输出：[
  {
    path: [
      0, 1, 2, 3, 4,
      4, 5, 6, 7, 8
    ],
    value: 126
  },
  {
    path: [
      0, 1, 1, 1, 2,
      3, 3, 4, 5, 5
    ],
    value: 123
  },
  {
    path: [
      0, 1, 1, 2, 3,
      4, 5, 6, 7, 8
    ],
    value: 122
  } ....]

动态规划法给出的结果是：

dp(0, 0)
输出：126

可见也是一致的。多次尝试，结果都是一致，这样就验证了我们想法的正确性。

作者：王源，徐思坤，陈贤邦

OR Paper Weekly 栏目将会从运筹学顶级期刊上选择一部分有趣的文章，对这些文章的主要研究内容进行一个概述/点评。OR Paper Weekly 的特点是 不做大而全的照搬，也未必都只选择优质的文章，而是精选一部分有趣的文章。辅之以科普/点评/吐槽的方式，让大家随时了解最新的科研动态。欢迎大家一起来 欣赏优质文章，学习脑洞文章，鄙视灌水文章。

精选论文 (一)

论文题目：Machine-Learning–Based Column Selection for Column Generation

期刊：Transportation Science

发表年份：2020年

作者：Mouad Morabit, Guy Desaulniers, Andrea Lodia

原文链接：

https://pubsonline.informs.org/doi/abs/10.1287/trsc.2021.1045

摘要：列生成算法是常见的求解大规模整数规划问题的方法。本文提出了一种基于机器学习来加速列生成求解的新方法。该方法的特点是 在每一步列生成算法的迭代中，通过机器学习训练出一个模型来生成新的列（选择一个变量的子集添加到主问题构成新的一列）。该方法通过选择最合适的列来降低主问题的计算时间。我们在带有时间窗约束的车辆路径规划问题和组环问题上验证了该方法。实验结果表明该方法可以降低30%的计算时间。

文章亮点/点评：传统的列生成算法是通过极大化 reduced cost 的子问题来生成新的列的。从全局的视角来看 这样添加列的方式会比较的greedy，因此传统列生成算法的效率会受到较大影响。本文通过机器学习来生成新的列，将机器学习和传统的运筹优化算法相结合，以基于数据的机器学习模型来解决传统运筹学中一些难以解决的问题已经成为了目前的研究趋势。

精选论文 (二)

论文题目：Metaheuristics “In the Large”

期刊：European Journal of Operational Research

发表年份：2021年

作者：Jerry Swan, Steven Adriaensen,  Alexander E. I. Brownlee c, Kevin Hammond d, Colin G. Johnson, et al.

原文链接：

https://www.sciencedirect.com/science/article/pii/S0377221721004707

摘要：在过去的几十年里涌现出大量关于元启发式算法（Metaheuristics）的各类改进，可以说元启发式算法早已成为了优化领域的显学。为了避免元启发式算法缺乏系统性和缺乏可重复性的问题，我们迫切需要强大的计算设备来支持新方法的开发，分析和比较。为此，我们提出了 Metaheuristics “In the Large” 项目，以实现我们的愿景。该项目的愿景是：实现真正可扩展的算法模版，无需对算法进行修改即可重用。为特定领域的知识（domain knowledge）的融入提供通用的白盒问题描述，以及可以远程访问的框架，组建和问题。相信这一举措会提高算法的可重复性，加速该领域的进一步发展。我们认为通过这种基础施设的建设，可以促使该领域朝着更水平的科学研究加速前进。在本文中我们描述了我们的期望，并且展示了所有元启发采用的通用协议，以帮助进一步释放该领域的潜力，简化对元启发式算法设计空间的探索。

文章亮点/点评：诚如本篇论文所述，元启发式算法在过去的几十年中出（guan）现（shui）了大量的论文。各种各样的生物在元启发式算法中粉墨登场，例如蛙跳算法、猫群算法、蟑螂算法，天牛算法，帝王蝶算法，头脑风暴算法，黏菌觅食算法，植物胞群算法，树种优化算法。由此可见，动物都不够用的了，元启发式算法已经将触手伸到了昆虫，植物和微生物。元启发式算法有着非常好的应用场景，例如在多目标优化，数据驱动优化和黑箱优化等场景下可以和数学优化算法形成互补。希望通过Metaheuristics “In the Large” 项目能促使元启发式算法的研究进一步的向着健康的方向发展。

精选论文 (三)

论文题目：Operations Research: Topics, Impact, and Trends from 1952–2019

期刊：Operations Research

发表年份：2021年

作者：Mouad Morabit, Guy Desaulniers, Andrea Lodia

原文链接：

https://pubsonline.informs.org/doi/abs/10.1287/opre.2021.2139

摘要：本文回顾了 Operations Research 期刊近68年来的出版成果，揭示了其文章、作者的变化及其随时间的影响，以及这些变化对研究人员和从业者的影响。Operations Research 自1952年创立到2019年，总共接收了5440篇文章。

通过对论文作者的数据进行分析可得以下结论：1美国、加拿大和英国分列发文的前三位，其中美国以4058篇名列第一，并且美国的发文量超过了第2名到第10名的总和。中国（大陆）以125篇名列第7名。2 库存问题是北美、亚洲和中东国家最热门的研究问题，而欧洲国家则专注于调度问题。

通过对研究内容的数据进行分析可得以下结论：1 数学规划是最常见的研究方法，2 库存问题是研究最多的问题，3近年来与定价有关的研究正在显著增加。

通过对研究趋势的数据进行分析可得以下结论：1 近十年来动态规划是最常用的方法，2 近十年来定价是研究最多的问题。由此可知动态规划和定价问题在未来依然会是研究热点。

文章亮点/点评：

关于本文的详细解读 可以参考之前我们的一篇文章：【主编推荐】哪个国家在 Operations Research 上发文最多？（数据解析OR创刊68年来发文趋势）

精选论文 (四)

论文题目：From Predictive to Prescriptive Analytics

期刊：Management Science

发表年份：2018年

作者：Dimitris Bertsimas, Nathan Kallusb

原文链接：

https://pubsonline.informs.org/doi/abs/10.1287/mnsc.2018.3253

摘要：我们构造了一个结合了机器学习、运筹与管理科学以及一些特殊的方法的框架，通过这个框架，我们可以利用数据为运筹与管理科学中的问题推荐（Prescribe）最优决策。与其他数据驱动优化的研究不同的是，我们考虑的数据不仅包含了对成本和收益有直接影响的观测量，还包含了问题相关的主要辅助变量的观测值。本文考虑的主要问题是一个基于不完美样本的条件随机优化问题，其中问题相关的联合概率分布未知。我们展示了我们的方法可广泛应用于非常多的决策问题，证明了这些方法在计算上易解并且在一般条件下渐进最优，这些结果即使在不满足独立同分布的数据以及删失数据下仍然成立。我们将这套方法拓展至某些决策变量（例如价格）与不确定性有未知因果关系的场景。我们定义了指导性系数（Coefficient of Prescriptiveness）P，用以从运营的角度度量数据的指导性内容以及策略的有效性。我们以一个年均发货量超过十亿单的大型媒体公司分销部门所面对的库存管理问题为例来展示这套方法。我们利用内部数据以及从IMDb、烂番茄和谷歌上收集得到的公共数据对他们的运营决策所作出的指导性建议提升了他们运营的基准表现。具体而言，利用我们的方法，我们收集到的数据在指导性系数的度量下为他们的运营能力带来了88%的提升。

文章亮点/点评：从描述性分析（Descriptive analysis）中的收集、整理、统计，到预测性分析（Predictive analysis）中更为复杂的挖掘、聚类、预测，我们对数据的了解越来越深。现在我们不仅需要理解数据，更要利用数据来做出更智能的决策，这就是本文所描述的指导性分析（Prescriptive Analysis）。数据驱动决策是目前绝大多数商业模式迈向智能与高效的关键，而实现这样的跨越需要背后完备的概率统计和优化理论的支撑。本文从概述利用数据推荐最优决策的框架出发，通过研究一类特殊的随机优化问题，以及它在经典的库存管理问题中的应用，展示了数据驱动决策在运营效率提升上的强大之处。虽然 Prescriptive Analysis一词有新瓶装旧酒之嫌，但无论是为了了解数据驱动决策的框架还是为了学习它在具体运营问题中的应用，本文都是一篇值得一读的经典。

精选论文 (五)

论文题目：Actor-Critic–Like Stochastic Adaptive Search for Continuous Simulation Optimization

期刊：Operation Research

发表年份：2021年

作者：Qi Zhang, Jiaqiao Hua

原文链接：

https://pubsonline.informs.org/doi/abs/10.1287/opre.2021.2214

摘要：我们为求解一类带Lipschitz连续性质的仿真优化问题提出了一个随机搜索算法。该算法从一个参数化的解空间概率分布中抽取候选解，并通过一个基于所谓的收敛球方法的异步学习过程衡量抽取的样本点（候选解）的表现。这个算法最大的特点在于它在完整地保留了此前仿真的信息的同时采用了一个近似架构来利用对目标函数的知识来提升解的质量。算法的每一步同时调整参数化的解空间概率分布以及目标函数的近似估计，这与强化学习中的Actor-Critic架构非常类似。我们提出了算法表现的有限时间概率边界，并展示了在每次迭代中仅有一个仿真观测值的情况下的全局收敛性。我们的实证研究表明该算法很有前景，并有可能在效率以及可靠性上超越目前的一些搜索方法。

文章亮点/点评：本文所提出的随机搜索算法，在传统随机优化的思路，即近似目标函数的基础上，赋予解空间一个先验并用随机抽样的方法去学习解空间的分布，来达到提升算法表现的目的。文章算法的思路非常清晰，并且有一定的启发性。与强化学习中Actor-Critic算法的类比非常抓人眼球。这也从侧面说明了一篇顶刊文章不仅需要方法论上的创新，还需要讲好一个故事，有一个抓人眼球的点。

精选论文 (六)

论文题目：Feature-driven Economic Improvement for Network-constrained Unit Commitment: A Closed-loop Predict-and-optimize Framework

期刊：IEEE Transactions on Power Systems

发表年份：2021年

作者：Xianbang Chen, Yafei Yang, Yikui Liu, Lei Wu

原文链接：

https://www.researchgate.net/profile/Xianbang-Chen/publication/356189174_Feature-driven_Economic_Improvement_for_Network-constrained_Unit_Commitment_A_Closed-loop_Predict-and-optimize_Framework/links/61906a2961f098772093c4d3/Feature-driven-Economic-Improvement-for-Network-constrained-Unit-Commitment-A-Closed-loop-Predict-and-optimize-Framework.pdf

Abstract:

机组组合(unit commitment, UC)是电力市场结算的重要环节。一般来说，UC会在开环型预测-再优化(Open-loop predict-then-optimize, O-PO)框架下执行：先以统计精度为导向预测可再生能源出力，再以系统运行成本最小为目标优化UC，预测与优化环节呈开环状态。但由于UC模型含有大量非线性约束，因此统计上更精确的预测值未必会诱导出经济性更好的UC决策。针对该问题，本文提出了一种闭环型预测-优化(Closed-loop predict-and-optimize, C-PO)框架以提高UC经济性。C-PO以UC经济性为导向进行预测，同时考虑了UC模型的约束与优化目标，使得预测与优化环节呈闭环状态。为提高C-PO的实用性，设计了拉格朗日分解算法，以加速C-PO预测器的训练过程。最后，基于IEEE RTS 24-Bus系统和5655-Bus系统的算例得到两个结论：C-PO框架能够有效提高UC经济性，并具有较好的实用性。

文章亮点/点评：

1.   生成以UC经济性为导向的可再生能源预测值，C-PO能够有效提高UC经济性。

2.   无需替换已有的预测工具，C-PO与当前普遍使用的O-PO框架兼容。

3.   拉格朗日分解算法能有效加速预测器的训练，C-PO具备应用于大型系统的潜力。