随机对偶动态规划 SDDP，报童模型的一个 python 例子

将 benders 分解、SAA、动态规划 方法结合在一起，就产生了随机对偶动态规划方法（stochastic dual dynamic programming, SDDP）。这两年在论文里看到了这个方法，一般用在运输或能源规划问题中，在库存管理问题中目前用的还不多。

1. SDDP 的假设与优点

应用该方法的两个重要假设是：

• 各阶段（stage）的随机性要相互独立
• 随机性具有有限个实现值 (finite realizations)

该方法的优点：

• 每次运行的算法复杂度只与单个阶段的情景数有关，而与总情景数无关

该方法一般是从大的情景树中随机抽取一定数量的情景，因此单个阶段的情景数量就是随机抽取的情景数量

2. 一般的多阶段随机规划模型与 SDDP 模型

2.1 一般的多阶段随机规划模型

一般的多阶段规划模型为：

stage $1$ 模型

$$\begin{array}{rlrl}z=& \underset{x_1}{\min }& & c_1^T{x}_1+\mathbb{E}[Q_2(x_1,{\omega }_2)]\\ & \text{s.t.}& & A_1{x}_1=b_1\\ & & & x_1\ge 0\end{array}$$

其中， ${\omega }_2$ 表示第二阶段随机性的一个具体实现值（第一阶段是不存在随机性的，即第一阶段的决策 $x_1$ 需要在随机性发生前做出），同理，${\omega }_t$ 表示第 $t$ 阶段随机性的一个具体实现值，这个随机性可以表现在 $A_t$，$b_t$ ， $E_t$ 或 $c_t$ 上。

stage t 模型：

$$\begin{array}{rlrl}{Q}_t(x_{t-1},w_t)=& \underset{x_t}{\min }& & c_t^T{x}_t+\mathbb{E}[Q_{t+1}(x_t,{\omega }_{t+1})]\\ & \text{s.t.}& & A_t{x}_t=b_t-E_t{x}_{t-1}\\ & & & x_t\ge 0\end{array}$$

最后 stage T 的模型：
$$\begin{array}{rlrl}{Q}_T(x_{T-1},w_T)=& \underset{x_T}{\min }& & c_T^T{x}_T\\ & \text{s.t.}& & A_T{x}_T=b_T-E_T{x}_{T-1}\\ & & & x_T\ge 0\end{array}$$

2.2 一般的 SDDP 模型

SDDP 方法在运用时，有几个特点：

• 每次迭代，有一个 forward 与 backward 过程
• 在 forward 中，得出每个阶段的解 $x_t^k$
• 在 backward 中，添加切平面约束条件到每个阶段的线性规划模型中。
• 不断循环上面的过程，直到统计意义的上界与下界比较接近为止，或者直到迭代达到一定次数为止

经我个人判断，这个 backward 其实也基本等价于上一个迭代过程中的 forward，区别在于考虑的样本量不一样：forward 时只针对抽取的样本计算，而backward时，每一个 stage 都对该 stage 所有的样本计算。

对应的 SDDP 模型：

第 $k+1$ 次迭代时，stage 1 的SDDP 模型

$$\begin{array}{rlrl}z=& \underset{x_1^{k+1},{\theta }_2^{k+1}}{\min }& & c_1^T{x}_1^{k+1}+{\theta }_2^{k+1}\\ & \text{s.t.}& & A_1{x}_1^{k+1}=b_1\\ & & & {\theta }_2^{k+1}+{\overline{\pi }}_{2,m}^T{E}_2{x}_1^{k+1}\ge {\overline{g}}_{2,m}m=1,2,\dots ,k\\ & & & {\theta }_2^{k+1}\ge \underline{\theta }\\ & & & x_1^{k+1}\ge 0\end{array}\text{\hspace{1em}(1)}$$

其中，${\theta }_2^{k+1}$ 为 stage 2模型的近似期望最优值，而 $\underline{\theta }$ 为它的一个给定的下界（可以自己定一个比较小的值）。

注意：在第 $k+1$ 次迭代时，前 $k$ 次的求解结果是已知的。因此，此时已经知道了第 $k$ 次迭代时的求解结果 $x_1^k$。

约束条件（1）为切平面条件。$m$ 表示第 $m$ 个切平面，在第 $k+1$ 次迭代时，有 $k$ 个切平面。这些切平面方程是前 $k$ 次迭代时，每次迭代添加的。因此，在第 1 次迭代的 forward 求解时，这个切平面约束条件还没有。

注意：要想得到 stage $t$ 第 $k$ 个切平面，需要在第 $k$ 次迭代时的第 stage $t+1$ 的规划模型求出。

这个切平面条件也可以认为是一个梯度下降的约束条件，梯度为 $-{\overline{\pi }}_{2,k}^T{E}_2$。其中，${\overline{\pi }}_{2,k}^T$ 是一个期望值：${\overline{\pi }}_{2,k}^T=\mathbb{E}({\pi }_{2,k}({\omega }_2))$，即为所有可能情景下的期望值，同理 ${\overline{g}}_{2,m}$ 也是一个期望值。而 ${\pi }_2({\omega }_2)$ 是其中一个情景 ${\omega }_2$ 第 2 阶段线性规划模型的对偶模型的解。

这个切平面条件也等价于 benders 分解中的 optimality cut，即在第 $k+1$ 次迭代时，第二阶段的最优解大于等于它对偶问题的目标函数值：
$$\begin{array}{rl}{\theta }_2& \ge {\pi }_2^T(b_2-E_2{x}_1^{k+1})\\ & ={\pi }_2^T{b}_2-{\pi }_2^T{E}_2{x}_1^k+{\pi }_2^T{E}_2{x}_1^k-{\pi }_2^T{E}_2{x}_1^{k+1}\\ & =Q_2^k-{\pi }_2^T{E}_2(x_1^{k+1}-x_1^k)\end{array}$$

上面的第一个不等式中，右边其实是对偶问题的目标函数值（个人感觉在求出第二阶段的线性规划后，可以通过求解器得到约束条件的对偶值 ${\pi }_2^T$，然后可以写出对偶问题的目标函数，最后将 ${\pi }_2^T$ 带入到这个目标函数中就能立刻写出这个切平面条件了）。

在第 $k$ 次迭代时中，stage 2 线性规划的近似模型为：

$$\begin{array}{rlrl}{\mathcal{Q}}_2(x_1^k,{\omega }_2)=& \underset{x_2^k,{\theta }_3^k}{\min }& & c_2^T{x}_2^k+{\theta }_3^k\\ & \text{s.t.}& & A_2{x}_2^k=b_2-E_2{x}_1^k[{\pi }_{2,k}({\omega }_2)]\\ & & & {\theta }_3^k+{\overline{\pi }}_{3,m}^T{E}_2{x}_2^k\ge {\overline{g}}_{3,m}m=1,2,\dots ,k-1[{\pi }_{2,k}({\omega }_2)]\\ & & & {\theta }_3^k\ge \underline{\theta }\\ & & & x_2^k\ge 0\end{array}$$

其中，${\mathcal{Q}}_2(x_1^k,{\omega }_2)$ 为 $Q_2(x_1^k,{\omega }_2)$ 的近似，${\theta }_3^k$ 是 $\mathbb{E}[Q_3(x_2^k,{\omega }_3)]$ 的近似。

$Q_2(x_1^k,{\omega }_2)$ 中约束条件 $A_2{x}_2^k=b_2-E_2{x}_1^k$ 对应的对偶求解变量为 ${\pi }_2^k({\omega }_2)$，可以证明，$Q_2(x_1^k,{\omega }_2)$ 的次梯度为：
$$\partial Q_2(x_1^k,{\omega }_2)=-E_2{\mathcal{D}}_2(x_1^k,{\omega }_2)$$

其中，${\mathcal{D}}_2(x_1^k,{\omega }_2)$ 为 $Q_2(x_1^k,{\omega }_2)$对偶问题的最优解（其实就是 ${\pi }_2^k({\omega }_2)$）。根据次梯度的定义，可以得到一个切平面约束条件：
$${\theta }_2^{k+1}-Q_2(x_1^k,{\omega }_2)\ge -E_2{\pi }_2^k({\omega }_2)(x_1^{k+1}-x_1^k)$$

实际计算时，用 ${\mathcal{Q}}_2(x_1^k,{\omega }_2)$ 近似替代 $Q_2(x_1^k,{\omega }_2)$，并且使用 ${\mathcal{Q}}_2(x_1^k,{\omega }_2)$ 对偶模型得到的 ${\pi }_2^k({\omega }_2)$。由于每一个随机性 ${\omega }_2$ 都有一个上面的切平面条件，求出它们的期望，就得到最终的切平面条件。

因此，为了求出第 $k+1$ 次迭代时 stage 1 的切平面条件，需要在第 $k$ 次迭代时，通过 stage 2 的规划模型，分别求出 ${\mathcal{Q}}_2(x_1^k,{\omega }_2)$ 以及 ${\pi }_{2,k}({\omega }_2)$。进而，取期望，得到切平面条件。

切平面条件中，截距 ${\overline{g}}_{2,k}$ 为：
$${\overline{g}}_{2,k}=\mathbb{E}[{\mathcal{Q}}_2(x_1^k,{\omega }_2)]+{\overline{\pi }}_{2,k}^T{E}_2{x}_1^k$$

切平面条件（等价于 benders 分解中的 optimality cut）也等价于梯度下降约束条件为：

$${\theta }_2^{k+1}\ge \mathbb{E}[{\mathcal{Q}}_2(x_1^k,{\omega }_2)]-{\overline{\pi }}_{2,m}^T{E}_2(x_1^{k+1}-x_1^k)m=1,2,\dots ,k$$

同理，第 $k+1$ 次迭代时，stage t 的SDDP 模型：
$$\begin{array}{rlrl}{\mathcal{Q}}_t(x_{t-1}^{k+1},{\omega }_t)=& \underset{x_t^{k+1},{\theta }_{t+1}^{k+1}}{\min }& & c_t^T{x}_t^{k+1}+{\theta }_{t+1}^{k+1}\\ & \text{s.t.}& & A_t{x}_t^{k+1}=b_t-E_t{x}_{t-1}^{k+1}\\ & & & {\theta }_{t+1}^{k+1}+{\overline{\pi }}_{t+1,m}^T{E}_{t+1}{x}_t^k\ge {\overline{g}}_{t+1,m}m=1,2,\dots ,k\\ & & & {\theta }_{t+1}^k\ge \underline{\theta }\\ & & & x_t^{k+1}\ge 0\end{array}$$

它的第 $k$ 个切平面约束，是通过在第 $k$ 次迭代时，第 stage $t+1$ 的规划模型分别求出 ${\mathcal{Q}}_{t+1}(x_t^k,{\omega }_{t+1})$ 及 ${\pi }_{t+1}^k(w_{t+1})$，对它们取期望得出。

规划模型
$$\begin{array}{rlrl}{\mathcal{Q}}_{t+1}(x_t^k,{\omega }_{t+1})=& \underset{x_{t+1}^k,{\theta }_{t+2}^k}{\min }& & c_{t+1}^T{x}_{t+1}^k+{\theta }_{t+2}^k\\ & \text{s.t.}& & A_{t+1}{x}_{t+1}^k=b_{t+1}-E_{t+1}{x}_t^k[{\pi }_{t+1,k}({\omega }_{t+1})]\\ & & & {\theta }_{t+2}^k+{\overline{\pi }}_{t+2,m}^T{E}_{t+2}{x}_t^k\ge {\overline{g}}_{t+1,m}m=1,2,\dots ,k-1[{\pi }_{t+1,k}({\omega }_{t+1})]\\ & & & {\theta }_{t+2}^k\ge \underline{\theta }\\ & & & x_{t+1}^k\ge 0\end{array}$$

第 $k+1$ 次迭代时，stage $T$ 阶 SDDP 模型：
$$\begin{array}{rlrl}{\mathcal{Q}}_T(x_{T-1}^{k+1},{\omega }_t)=& \underset{x_T^{k+1}}{\min }& & c_T^T{x}_T^{k+1}\\ & \text{s.t.}& & A_t{x}_T^{k+1}=b_T-E_T{x}_{T-1}^{k+1}[{\pi }_{T,k+1}]\\ & & & x_T^{k+1}\ge 0\end{array}$$

最后一 stage $T$ 就没有切平面约束了，直接求解。当然，这个模型可以得到一个对偶变量 ${\pi }_{T,k+1}$ 以及 ${\mathcal{Q}}_T(x_{T-1}^{k+1},{\omega }_t)$，可以用来构造第 $k+2$ 次迭代时第 stage $T-1$ 的切平面约束条件。

需要注意的是

• SDDP 最关键的步骤是添加切平面条件，而切平面条件与对偶模型的构建与求解非常相关，极易出错
• 第 $k+1$ 次迭代时，stage $t$ 时，切平面条件中对偶变量 $\pi$ 不仅包括 stage $t+1$原线性规划模型中的约束条件对应的对偶变量值，还包括之前迭代时 stage $t+1$ 每次添加的切平面条件对应的对偶变量值
• 求解变量有非零的上下界时，最好把上下界放到约束条件里，这样求解器就不会遗漏一些对偶变量了
• 添加切平面条件时，原问题目标函数中的常数值，以及上一个 stage 的求解变量不能遗忘

3. 报童模型的 SDDP 模型

下面以一个 2 阶段的报童模型为例，构造它的 SDDP 求解模型。首先用最传统的报童模型，缺货能补，不考虑产品的销售价格，也不考虑最后一阶段的剩余库存残值。

报童模型中的参数：

• 单位进货成本 $c$
• 单位库存持有成本 $h$
• 单位延迟补货成本（backorder cost）$b$
• $t$ 阶段的订货量 $q_t$
• 初始库存量是 $I_0$，$t$ 阶段期末库存量 $I_t$

如果是单阶段报童模型的话，最优订货点为：

$$q^{\ast }=\frac{b-c}{b+h}$$

下面构建 2 个阶段（3 个stage）的SDDP模型

k+1 次迭代时，stage 1 的模型：

$$\begin{array}{rlrl}z=& \underset{q_1^{k+1},{\theta }_2^{k+1}}{\min }& & c_1{q}_1^{k+1}+{\theta }_2^{k+1}\\ & \text{s.t.}& & {\theta }_2^{k+1}\ge \mathbb{E}[{\mathcal{F}}_2(q_1^k,d_1^s)]m=1,2,\dots ,k\\ & & & {\theta }_2^{k+1}\ge \underline{\theta }\\ & & & q_1^{k+1}\ge 0\end{array}$$

(stage 1 的模型不用写它的对偶问题了，而接下来的 stage 模型都要涉及到约束条件的对偶)

上面的公式中，用 $\mathcal{F}$ 表示子问题对偶问题的目标函数值，这一项是切平面条件。

k+1 次迭代时，stage 2 的模型：

$$\begin{array}{rlrl}{\mathcal{Q}}_2(q_1^{k+1,s},d_1^s)=& \underset{I_1^{k+1},B_1^{k+1},q_2^{k+1}}{\min }& & h{I}_1^{k+1}+b{B}_1^{k+1}+c{q}_2^{k+1}+{\theta }_3^{k+1}\\ & & & \\ & \text{s.t.}& & I_1^{k+1}-B_1^{k+1}=I_0+q_1^{k+1}-d_1^s[{\pi }_{2,1}^s]\\ & & & {\theta }_2^{k+1}\ge \mathbb{E}[{\mathcal{F}}_3(q_2^k,d_2^s)]m=1,2,\dots ,k[{\pi }_{2,2}^s]\\ & & & {\theta }_2^{k+1}\ge \underline{\theta }\\ & & & I_1^{k+1}\ge 0,B_1^{k+1}\ge 0,q_2^{k+1}\ge 0\end{array}$$
可以据此写出它对偶问题的目标函数，即
$${\mathcal{F}}_2={\pi }_{2,1}^s(I_0+q_1^{k+1}-d_1^s)-\sum {\pi }_{2,2}^s{\pi }_3^s{d}_2^s$$

$d_1^s$ 表示在情景 $s$ 的第一阶段需求值，需求注意的是：求解得到的变量 $I_1^{k+1}$, $B_1^{k+1}$ , $q_2^{k+1}$ 其实也是跟情景 $s$ 相关的。

• 上面表达式的第二项不能丢，否则出错。因为之前迭代时添加的切平面条件也有对应的对偶值

k+1 次迭代时，stage 3 的模型：

$$\begin{array}{rlrl}{\mathcal{Q}}_3(q_2^{k+1,s},d_1^s)=& \underset{I_2^{k+1},B_2^{k+1}}{\min }& & h{I}_2^{k+1}+b{B}_2^{k+1}\\ & \text{s.t.}& & I_2^{k+1}-B_2^{k+1}=I_1^{k+1}+q_2^{k+1}-d_2^s[{\pi }_3^s]\\ & & & I_2^{k+1}\ge 0,B_2^{k+1}\ge 0\end{array}$$

可以据此写出它对偶问题的目标函数，即
$${\mathcal{F}}_3={\pi }_3^s(I_1^{k+1}+q_2^{k+1}-d_2^s)$$

4. 求解单阶段报童模型的一个例子

我用这个方法试了试求解单阶段的报童，求解结果不错。但应用到多阶段时编程要复杂多了，有几点需要注意。

• 约束条件要都是线性表达式，含有 max 或 min 的表达式要处理下
• 若有整数求解变量，更加复杂，似乎要用到 sddip 方法

单阶段报童模型其实也是按照一个两阶段模型来求解：第一阶段的决策变量是订货量（需求已知前），而第二阶段的决策变量是销售量（需求已知后）。

当然，报童模型的最终求解变量其实是订货量。

k+1 次迭代时，stage 1 的模型：

$$\begin{array}{rlrl}z=& \underset{q_1^{k+1},{\theta }_2^{k+1}}{\min }& & c_1{q}_1^{k+1}+{\theta }_2^{k+1}\\ & \text{s.t.}& & {\theta }_2^{k+1}\ge \mathbb{E}[{\mathcal{Q}}_2(q_1^k,d_1^s)]-{\pi }_2^{k+1,s}(q_1^{k+1}-q_1^k)m=1,2,\dots ,k\\ & & & {\theta }_2^{k+1}\ge \underline{\theta }\\ & & & q_1^{k+1}\ge 0\end{array}$$

k+1 次迭代时，stage 2 的模型（此时的决策变量是销量）：
$$\begin{array}{rlrl}{\mathcal{Q}}_2(q_1^{k+1,s},d_1^s)=& \underset{y_1^{k+1}}{\min }& & -p{y}_1^{k+1}-\gamma (q_1^{k+1,s}-y_1^{k+1,s})\\ & \text{s.t.}& & -y_1^{k+1}\ge -q_1^{k+1,s}\\ & & & -y_1^{k+1}\ge -d_2^s\\ & & & y_1^{k+1}\ge 0\end{array}$$

也可以写成（此时的决策变量是库存量）：
$$\begin{array}{rlrl}{\mathcal{Q}}_2(q_1^{k+1,s},d_1^s)=& \underset{I_1^{k+1}}{\min }& & -\gamma I_1^{k+1}-p(q_1^{k+1,s}-I_1^{k+1,s})\\ & \text{s.t.}& & I_1^{k+1}\ge q_1^{k+1,s}-d_1^s\\ & & & -I_1^{k+1}\ge -q_1^{k+1,s}\\ & & & I_1^{k+1}\ge 0\end{array}$$

之前漏写了一个约束条件，一直求不对割平面约束条件，上面模型中的约束条件一个都不能少。

它们的对偶模型分别为：
$$\begin{array}{rlrl}& \underset{{\pi }_1^{k+1},{\pi }_2^{k+1}}{\max }& & -{\pi }_1^{k+1}{q}_1^{k+1}-{\pi }_2^{k+1}{d}_1^s-\gamma q_1^{k+1}\\ & \text{s.t.}& & -{\pi }_1^{k+1}-{\pi }_2^{k+1}\le -p+\gamma \\ & & & -{\pi }_1^{k+1}\ge 0,-{\pi }_2^{k+1}\ge 0\end{array}$$

以及：
$$\begin{array}{rlrl}& \underset{{\pi }_1^{k+1},{\pi }_2^{k+1}}{\max }& & {\pi }_1^{k+1}(q_1^{k+1}-d_1^s)-{\pi }_2^{k+1}{q}_1^{k+1}-p{q}_1^{k+1}\\ & \text{s.t.}& & {\pi }_1^{k+1}-{\pi }_2^{k+1}\le p-\gamma \\ & & & {\pi }_1^{k+1}\ge 0,{\pi }_2^{k+1}\ge 0\end{array}$$

几点总结：

• 对偶变量 $\pi$ 的求解非常重要，最好直接将对偶模型写出来，然后添加类似 benders 的 optimality cut

用 SDDP 求解单阶段报童模型的一个例子：

# -*- coding: utf-8 -*-
"""
Created on Sat Jan  8 10:57:33 2022

@author: zhen chen
@Email: chen.zhen5526@gmail.com

MIT Licence.

Python version: 3.8


Description: test sddp in a single period news vender problem
    
"""


import numpy as np
import scipy.stats as st
from gurobipy import *
import math
import random


# generate latin hypercube samples 
def generate_sample(sample_num, trunQuantile, mu):
    samples = [0 for i in range(sample_num)]
    for i in range(sample_num):
        rand_p = np.random.uniform(trunQuantile*i/sample_num, trunQuantile*(i+1)/sample_num)
        samples[i] = st.poisson.ppf(rand_p, mu)
    return samples


price  =  6
vari_cost = 2
sal_value = 1
mean_demands = 10
sample_num = 100
trunQuantile = 0.9999 # affective to the final ordering quantity

samples = generate_sample(sample_num, trunQuantile, mean_demands)
   
ini_Q = 6
ini_Obj = float("inf") # largest float number

## compute seconde stage objective in each scenario
N = len(samples)

k = 1
feasible_cut_index = []
Q = ini_Q
T= [1, 0]
obj = [0.0  for i in range(N)]
g = [0.0 for i in range(N)]
pi1 = [0.0  for i in range(N)]
Dpi2 = [0.0 for i in range(N)]

m2 = Model()
nita = m2.addVar(lb = -GRB.INFINITY, vtype = GRB.CONTINUOUS) # objective of the second stage
nita_value = -3000
m2.addConstr(nita >= nita_value)
x = m2.addVar(vtype = GRB.CONTINUOUS)
this_obj = vari_cost*x + nita
m2.setObjective(this_obj, GRB.MINIMIZE)

while True:
    last_nita = nita_value
    last_master_obj = nita_value + vari_cost*Q
    for i in range(N):
        m = Model()    
        y = m.addVar(vtype = GRB.CONTINUOUS)
        this_s_obj = -price*y - sal_value*(Q - y) # or sal_value*(Q - y) , max(0, Q-samples[i])
        m.setObjective(this_s_obj, GRB.MINIMIZE)
        
        m.addConstr(-y >= -Q)
        m.addConstr(-y >= -samples[i])
        m.optimize()       
        y_value = y.x
        
        obj[i] = m.objVal
        pi = m.getAttr(GRB.Attr.Pi)
        g[i] = sum([a*b for a,b in zip(pi, T)])
        pi1[i] = -pi[0]
        Dpi2[i] = -pi[1] * samples[i]
        print()    
   
    avg_obj = sum(obj)/N    
    M = round(0.5*N)
    obj2 = random.choices(obj, k=M)
    avg_obj2 = sum(obj2)/N 
    var_obj2 = sum([(obj2[i] - avg_obj2)**2 for i in range(M)]) / (M-1)
    std_obj2 = math.sqrt(var_obj2)
    avg_g = sum(g)/N
    avg_pi1 = sum(pi1)/N
    avg_Dpi2 = sum(Dpi2)/N
    m2.addConstr(nita >= avg_obj - avg_g*(x-Q)) # add cut
    #m2.addConstr(nita >= avg_pi1*x + avg_Dpi2) # same as the above constraint
    m2.update()
    m2.optimize()
    
    Q = x.x
    nita_value = nita.x
    master_obj = m2.objVal
    print()
    
    k = k + 1
    # if abs(nita_value - avg_obj) < 1e-2 and k > 4: 
    #     break

    lb = avg_obj2 - 1.96*std_obj2/math.sqrt(M)
    ub = avg_obj2 + 1.96*std_obj2/math.sqrt(M)    
    # print(nita_value)
    # print(lb)
    # print(ub)  
    # print(Q)
    # if nita_value >= lb and nita_value <= ub:
    #     break
    # if abs(ub-lb) < 0.8:
    #     break
    if abs(last_master_obj - master_obj) < 1e-2 and k > 4:  # 这个终止条件最好
        break


print('iteration steps are %d' % k)    
print('ordering quantity is %.2f' % Q)
print('expected profit is %.2f' % master_obj)