文章目录
将 benders 分解、SAA、动态规划 方法结合在一起,就产生了随机对偶动态规划方法(stochastic dual dynamic programming, SDDP)。这两年在论文里看到了这个方法,一般用在运输或能源规划问题中,在库存管理问题中目前用的还不多。
1. SDDP 的假设与优点
应用该方法的两个重要假设是:
- 各阶段(stage)的随机性要相互独立
- 随机性具有有限个实现值 (finite realizations)
该方法的优点:
- 每次运行的算法复杂度只与单个阶段的情景数有关,而与总情景数无关
该方法一般是从大的情景树中随机抽取一定数量的情景,因此单个阶段的情景数量就是随机抽取的情景数量
2. 一般的多阶段随机规划模型与 SDDP 模型
2.1 一般的多阶段随机规划模型
一般的多阶段规划模型为:
stage 1 1 1 模型
z = min x 1 c 1 T x 1 + E [ Q 2 ( x 1 , ω 2 ) ] s.t. A 1 x 1 = b 1 x 1 ≥ 0 \begin{aligned} z=&\min_{x_1}~~ &&c_1^Tx_1+\mathbb{E}[Q_2(x_1, \omega_2)]\\ &\text{s.t.} &&A_1x_1=b_1\\ & &&x_1\geq 0 \end{aligned} z=x1min s.t.c1Tx1+E[Q2(x1,ω2)]A1x1=b1x1≥0
其中, ω 2 \omega_2 ω2 表示第二阶段随机性的一个具体实现值(第一阶段是不存在随机性的,即第一阶段的决策 x 1 x_1 x1 需要在随机性发生前做出),同理, ω t \omega_t ωt 表示第 t t t 阶段随机性的一个具体实现值,这个随机性可以表现在 A t A_t At, b t b_t bt , E t E_t Et 或 c t c_t ct 上。
stage t 模型:
Q t ( x t − 1 , w t ) = min x t c t T x t + E [ Q t + 1 ( x t , ω t + 1 ) ] s.t. A t x t = b t − E t x t − 1 x t ≥ 0 \begin{aligned} Q_t(x_{t-1}, w_t)=&\min_{x_t}~~ &&c_t^Tx_t+\mathbb{E}[Q_{t+1}(x_{t}, \omega_{t+1})]\\ &\text{s.t.} &&A_tx_t=b_t-E_{t}x_{t-1}\\ & &&x_t\geq 0 \end{aligned} Qt(xt−1,wt)=xtmin s.t.ctTxt+E[Qt+1(xt,ωt+1)]Atxt=bt−Etxt−1xt≥0
最后 stage T 的模型:
Q
T
(
x
T
−
1
,
w
T
)
=
min
x
T
c
T
T
x
T
s.t.
A
T
x
T
=
b
T
−
E
T
x
T
−
1
x
T
≥
0
\begin{aligned} Q_T(x_{T-1}, w_T)=&\min_{x_T}~~ &&c_T^Tx_T\\ &\text{s.t.} &&A_Tx_T=b_T-E_{T}x_{T-1}\\ & &&x_T\geq 0 \end{aligned}
QT(xT−1,wT)=xTmin s.t.cTTxTATxT=bT−ETxT−1xT≥0
2.2 一般的 SDDP 模型
SDDP 方法在运用时,有几个特点:
- 每次迭代,有一个 forward 与 backward 过程
- 在 forward 中,得出每个阶段的解 x t k x_t^k xtk
- 在 backward 中,添加切平面约束条件到每个阶段的线性规划模型中。
- 不断循环上面的过程,直到统计意义的上界与下界比较接近为止,或者直到迭代达到一定次数为止
经我个人判断,这个 backward 其实也基本等价于上一个迭代过程中的 forward,区别在于考虑的样本量不一样:forward 时只针对抽取的样本计算,而backward时,每一个 stage 都对该 stage 所有的样本计算。
对应的 SDDP 模型:
第 k + 1 k+1 k+1 次迭代时,stage 1 的SDDP 模型
z = min x 1 k + 1 , θ 2 k + 1 c 1 T x 1 k + 1 + θ 2 k + 1 s.t. A 1 x 1 k + 1 = b 1 θ 2 k + 1 + π ‾ 2 , m T E 2 x 1 k + 1 ≥ g ‾ 2 , m m = 1 , 2 , … , k θ 2 k + 1 ≥ θ ‾ x 1 k + 1 ≥ 0 (1) \begin{aligned} z=&\min_{x^{k+1}_1, \theta_2^{k+1}}~~ &&c_1^Tx^{k+1}_1+\theta_2^{k+1}\\ &\text{s.t.} &&A_1x^{k+1}_1=b_1\\ &&& \theta_2^{k+1}+\overline{\pi}^T_{2,m}E_2x^{k+1}_1\geq \overline{g}_{2,m}\quad m=1,2,\dots,k\tag{1}\\ &&& \theta_2^{k+1}\geq \underline{\theta}\\ & &&x^{k+1}_1\geq 0 \end{aligned} z=x1k+1,θ2k+1min s.t.c1Tx1k+1+θ2k+1A1x1k+1=b1θ2k+1+π2,mTE2x1k+1≥g2,mm=1,2,…,kθ2k+1≥θx1k+1≥0(1)
其中, θ 2 k + 1 \theta_2^{k+1} θ2k+1 为 stage 2模型的近似期望最优值,而 θ ‾ \underline{\theta} θ 为它的一个给定的下界(可以自己定一个比较小的值)。
注意:在第 k + 1 k+1 k+1 次迭代时,前 k k k 次的求解结果是已知的。因此,此时已经知道了第 k k k 次迭代时的求解结果 x 1 k x_1^k x1k。
约束条件(1)为切平面条件。 m m m 表示第 m m m 个切平面,在第 k + 1 k+1 k+1 次迭代时,有 k k k 个切平面。这些切平面方程是前 k k k 次迭代时,每次迭代添加的。因此,在第 1 次迭代的 forward 求解时,这个切平面约束条件还没有。
注意:要想得到 stage t t t 第 k k k 个切平面,需要在第 k k k 次迭代时的第 stage t + 1 t+1 t+1 的规划模型求出。
这个切平面条件也可以认为是一个梯度下降的约束条件,梯度为 − π ‾ 2 , k T E 2 -\overline{\pi}^T_{2,k}E_2 −π2,kTE2。其中, π ‾ 2 , k T \overline{\pi}^T_{2,k} π2,kT 是一个期望值: π ‾ 2 , k T = E ( π 2 , k ( ω 2 ) ) \overline{\pi}^T_{2,k}=\mathbb{E}(\pi_{2,k}(\omega_2)) π2,kT=E(π2,k(ω2)),即为所有可能情景下的期望值,同理 g ‾ 2 , m \overline{g}_{2,m} g2,m 也是一个期望值。而 π 2 ( ω 2 ) \pi_2(\omega_2) π2(ω2) 是其中一个情景 ω 2 \omega_2 ω2 第 2 阶段线性规划模型的对偶模型的解。
这个切平面条件也等价于 benders 分解中的 optimality cut,即在第
k
+
1
k+1
k+1 次迭代时,第二阶段的最优解大于等于它对偶问题的目标函数值:
θ
2
≥
π
2
T
(
b
2
−
E
2
x
1
k
+
1
)
=
π
2
T
b
2
−
π
2
T
E
2
x
1
k
+
π
2
T
E
2
x
1
k
−
π
2
T
E
2
x
1
k
+
1
=
Q
2
k
−
π
2
T
E
2
(
x
1
k
+
1
−
x
1
k
)
\begin{aligned} \theta_2&\geq \pi_2^T(b_2-E_2x_1^{k+1})\\ &=\pi_2^Tb_2-\pi_2^TE_2x_1^k+\pi_2^TE_2x_1^k-\pi_2^TE_2x_1^{k+1}\\ &=Q_2^k-\pi_2^TE_2(x_1^{k+1}-x_1^k) \end{aligned}
θ2≥π2T(b2−E2x1k+1)=π2Tb2−π2TE2x1k+π2TE2x1k−π2TE2x1k+1=Q2k−π2TE2(x1k+1−x1k)
上面的第一个不等式中,右边其实是对偶问题的目标函数值(个人感觉在求出第二阶段的线性规划后,可以通过求解器得到约束条件的对偶值 π 2 T \pi_2^T π2T,然后可以写出对偶问题的目标函数,最后将 π 2 T \pi_2^T π2T 带入到这个目标函数中就能立刻写出这个切平面条件了)。
在第 k k k 次迭代时中,stage 2 线性规划的近似模型为:
Q 2 ( x 1 k , ω 2 ) = min x 2 k , θ 3 k c 2 T x 2 k + θ 3 k s.t. A 2 x 2 k = b 2 − E 2 x 1 k [ π 2 , k ( ω 2 ) ] θ 3 k + π ‾ 3 , m T E 2 x 2 k ≥ g ‾ 3 , m m = 1 , 2 , … , k − 1 [ π 2 , k ( ω 2 ) ] θ 3 k ≥ θ ‾ x 2 k ≥ 0 \begin{aligned} \mathcal{Q}_2(x^k_{1}, \omega_2)=&\min_{x_2^k,\theta_3^k}~~ &&c_2^Tx^k_2+\theta_3^k\\ &\text{s.t.} &&A_2x_2^k=b_2-E_{2}x_{1}^k\qquad~~~~~~ [\pi_{2,k}(\omega_2)] \\ &&& \theta_3^{k}+\overline{\pi}^T_{3,m}E_2x^{k}_2\geq \overline{g}_{3,m}\quad m=1,2,\dots,k-1\qquad~~~~~~ [\pi_{2,k}(\omega_2)]\\ &&& \theta_3^k\geq \underline{\theta}\\ & &&x_2^k\geq 0 \end{aligned} Q2(x1k,ω2)=x2k,θ3kmin s.t.c2Tx2k+θ3kA2x2k=b2−E2x1k [π2,k(ω2)]θ3k+π3,mTE2x2k≥g3,mm=1,2,…,k−1 [π2,k(ω2)]θ3k≥θx2k≥0
其中, Q 2 ( x 1 k , ω 2 ) \mathcal{Q}_2(x^k_{1}, \omega_2) Q2(x1k,ω2) 为 Q 2 ( x 1 k , ω 2 ) Q_2(x^k_1, \omega_2) Q2(x1k,ω2) 的近似, θ 3 k \theta^k_3 θ3k 是 E [ Q 3 ( x 2 k , ω 3 ) ] \mathbb{E}[Q_3(x^k_2, \omega_3)] E[Q3(x2k,ω3)] 的近似。
Q
2
(
x
1
k
,
ω
2
)
Q_2(x^k_1, \omega_2)
Q2(x1k,ω2) 中约束条件
A
2
x
2
k
=
b
2
−
E
2
x
1
k
A_2x_2^k=b_2-E_{2}x_{1}^k
A2x2k=b2−E2x1k 对应的对偶求解变量为
π
2
k
(
ω
2
)
\pi_2^k(\omega_2)
π2k(ω2),可以证明,
Q
2
(
x
1
k
,
ω
2
)
Q_2(x^k_{1}, \omega_2)
Q2(x1k,ω2) 的次梯度为:
∂
Q
2
(
x
1
k
,
ω
2
)
=
−
E
2
D
2
(
x
1
k
,
ω
2
)
\partial Q_2(x^k_{1}, \omega_2)=-E_2\mathcal{D}_2(x^k_{1}, \omega_2)
∂Q2(x1k,ω2)=−E2D2(x1k,ω2)
其中,
D
2
(
x
1
k
,
ω
2
)
\mathcal{D}_2(x^k_{1}, \omega_2)
D2(x1k,ω2) 为
Q
2
(
x
1
k
,
ω
2
)
Q_2(x^k_{1}, \omega_2)
Q2(x1k,ω2)对偶问题的最优解(其实就是
π
2
k
(
ω
2
)
\pi_2^k(\omega_2)
π2k(ω2))。根据次梯度的定义,可以得到一个切平面约束条件:
θ
2
k
+
1
−
Q
2
(
x
1
k
,
ω
2
)
≥
−
E
2
π
2
k
(
ω
2
)
(
x
1
k
+
1
−
x
1
k
)
\theta_2^{k+1}-Q_2(x^k_1, \omega_2)\geq -E_2\pi_2^k(\omega_2)(x_1^{k+1}-x_1^k)
θ2k+1−Q2(x1k,ω2)≥−E2π2k(ω2)(x1k+1−x1k)
实际计算时,用 Q 2 ( x 1 k , ω 2 ) \mathcal{Q}_2(x^k_{1}, \omega_2) Q2(x1k,ω2) 近似替代 Q 2 ( x 1 k , ω 2 ) Q_2(x^k_1, \omega_2) Q2(x1k,ω2),并且使用 Q 2 ( x 1 k , ω 2 ) \mathcal{Q}_2(x^k_{1}, \omega_2) Q2(x1k,ω2) 对偶模型得到的 π 2 k ( ω 2 ) \pi_2^k(\omega_2) π2k(ω2)。由于每一个随机性 ω 2 \omega_2 ω2 都有一个上面的切平面条件,求出它们的期望,就得到最终的切平面条件。
因此,为了求出第 k + 1 k+1 k+1 次迭代时 stage 1 的切平面条件,需要在第 k k k 次迭代时,通过 stage 2 的规划模型,分别求出 Q 2 ( x 1 k , ω 2 ) \mathcal{Q}_2(x^k_1, \omega_2) Q2(x1k,ω2) 以及 π 2 , k ( ω 2 ) \pi_{2,k}(\omega_2) π2,k(ω2)。进而,取期望,得到切平面条件。
切平面条件中,截距
g
‾
2
,
k
\overline{g}_{2,k}
g2,k 为:
g
‾
2
,
k
=
E
[
Q
2
(
x
1
k
,
ω
2
)
]
+
π
‾
2
,
k
T
E
2
x
1
k
\overline{g}_{2,k}=\mathbb{E}[\mathcal{Q}_2(x_1^k, \omega_2)]+\overline{\pi}^T_{2,k}E_2x_1^k
g2,k=E[Q2(x1k,ω2)]+π2,kTE2x1k
切平面条件(等价于 benders 分解中的 optimality cut)也等价于梯度下降约束条件为:
θ 2 k + 1 ≥ E [ Q 2 ( x 1 k , ω 2 ) ] − π ‾ 2 , m T E 2 ( x 1 k + 1 − x 1 k ) m = 1 , 2 , … , k \theta_2^{k+1}\geq \mathbb{E}[\mathcal{Q}_2(x^k_1, \omega_2)]- \overline{\pi}^T_{2,m}E_2(x^{k+1}_1-x_1^k)\qquad m=1,2,\dots,k θ2k+1≥E[Q2(x1k,ω2)]−π2,mTE2(x1k+1−x1k)m=1,2,…,k
同理,第
k
+
1
k+1
k+1 次迭代时,stage t 的SDDP 模型:
Q
t
(
x
t
−
1
k
+
1
,
ω
t
)
=
min
x
t
k
+
1
,
θ
t
+
1
k
+
1
c
t
T
x
t
k
+
1
+
θ
t
+
1
k
+
1
s.t.
A
t
x
t
k
+
1
=
b
t
−
E
t
x
t
−
1
k
+
1
θ
t
+
1
k
+
1
+
π
‾
t
+
1
,
m
T
E
t
+
1
x
t
k
≥
g
‾
t
+
1
,
m
m
=
1
,
2
,
…
,
k
θ
t
+
1
k
≥
θ
‾
x
t
k
+
1
≥
0
\begin{aligned} \mathcal{Q}_t(x^{k+1}_{t-1}, \omega_{t})=&\min_{x_t^{k+1},\theta_{t+1}^{k+1}}~~ &&c_t^Tx^{k+1}_t+\theta_{t+1}^{k+1}\\ &\text{s.t.} &&A_tx_t^{k+1}=b_t-E_{t}x_{t-1}^{k+1}\\ &&& \theta_{t+1}^{k+1}+\overline{\pi}^{T}_{t+1,m}E_{t+1}x^{k}_t\geq \overline{g}_{t+1,m}\quad m=1,2,\dots,k\\ &&& \theta_{t+1}^k\geq \underline{\theta}\\ & &&x_t^{k+1}\geq 0 \end{aligned}
Qt(xt−1k+1,ωt)=xtk+1,θt+1k+1min s.t.ctTxtk+1+θt+1k+1Atxtk+1=bt−Etxt−1k+1θt+1k+1+πt+1,mTEt+1xtk≥gt+1,mm=1,2,…,kθt+1k≥θxtk+1≥0
它的第 k k k 个切平面约束,是通过在第 k k k 次迭代时,第 stage t + 1 t+1 t+1 的规划模型分别求出 Q t + 1 ( x t k , ω t + 1 ) \mathcal{Q}_{t+1}(x^{k}_{t}, \omega_{t+1}) Qt+1(xtk,ωt+1) 及 π t + 1 k ( w t + 1 ) \pi_{t+1}^k(w_{t+1}) πt+1k(wt+1),对它们取期望得出。
规划模型
Q
t
+
1
(
x
t
k
,
ω
t
+
1
)
=
min
x
t
+
1
k
,
θ
t
+
2
k
c
t
+
1
T
x
t
+
1
k
+
θ
t
+
2
k
s.t.
A
t
+
1
x
t
+
1
k
=
b
t
+
1
−
E
t
+
1
x
t
k
[
π
t
+
1
,
k
(
ω
t
+
1
)
]
θ
t
+
2
k
+
π
‾
t
+
2
,
m
T
E
t
+
2
x
t
k
≥
g
‾
t
+
1
,
m
m
=
1
,
2
,
…
,
k
−
1
[
π
t
+
1
,
k
(
ω
t
+
1
)
]
θ
t
+
2
k
≥
θ
‾
x
t
+
1
k
≥
0
\begin{aligned} \mathcal{Q}_{t+1}(x^{k}_{t}, \omega_{t+1})=&\min_{x_{t+1}^{k},\theta_{t+2}^{k}}~~ &&c_{t+1}^Tx^{k}_{t+1}+\theta_{t+2}^{k}\\ &\text{s.t.} &&A_{t+1}x_{t+1}^{k}=b_{t+1}-E_{t+1}x_{t}^{k}\qquad~~ [\pi_{t+1,k}(\omega_{t+1})]\\ &&& \theta_{t+2}^{k}+\overline{\pi}^T_{t+2,m}E_{t+2}x^{k}_t\geq \overline{g}_{t+1,m}\quad m=1,2,\dots,k-1\qquad~~ [\pi_{t+1,k}(\omega_{t+1})]\\ &&& \theta_{t+2}^k\geq \underline{\theta}\\ & &&x_{t+1}^{k}\geq 0 \end{aligned}
Qt+1(xtk,ωt+1)=xt+1k,θt+2kmin s.t.ct+1Txt+1k+θt+2kAt+1xt+1k=bt+1−Et+1xtk [πt+1,k(ωt+1)]θt+2k+πt+2,mTEt+2xtk≥gt+1,mm=1,2,…,k−1 [πt+1,k(ωt+1)]θt+2k≥θxt+1k≥0
第
k
+
1
k+1
k+1 次迭代时,stage
T
T
T 阶 SDDP 模型:
Q
T
(
x
T
−
1
k
+
1
,
ω
t
)
=
min
x
T
k
+
1
c
T
T
x
T
k
+
1
s.t.
A
t
x
T
k
+
1
=
b
T
−
E
T
x
T
−
1
k
+
1
[
π
T
,
k
+
1
]
x
T
k
+
1
≥
0
\begin{aligned} \mathcal{Q}_T(x^{k+1}_{T-1}, \omega_{t})=&\min_{x_T^{k+1}}~~ &&c_T^Tx^{k+1}_T\\ &\text{s.t.} &&A_tx_T^{k+1}=b_T-E_{T}x_{T-1}^{k+1}\qquad[\pi_{T,k+1}]\\ & &&x_T^{k+1}\geq 0 \end{aligned}
QT(xT−1k+1,ωt)=xTk+1min s.t.cTTxTk+1AtxTk+1=bT−ETxT−1k+1[πT,k+1]xTk+1≥0
最后一 stage T T T 就没有切平面约束了,直接求解。当然,这个模型可以得到一个对偶变量 π T , k + 1 \pi_{T,k+1} πT,k+1 以及 Q T ( x T − 1 k + 1 , ω t ) \mathcal{Q}_T(x^{k+1}_{T-1}, \omega_{t}) QT(xT−1k+1,ωt),可以用来构造第 k + 2 k+2 k+2 次迭代时第 stage T − 1 T-1 T−1 的切平面约束条件。
需要注意的是
- SDDP 最关键的步骤是添加切平面条件,而切平面条件与对偶模型的构建与求解非常相关,极易出错
- 第 k + 1 k+1 k+1 次迭代时,stage t t t 时,切平面条件中对偶变量 π \pi π 不仅包括 stage t + 1 t+1 t+1原线性规划模型中的约束条件对应的对偶变量值,还包括之前迭代时 stage t + 1 t+1 t+1 每次添加的切平面条件对应的对偶变量值
- 求解变量有非零的上下界时,最好把上下界放到约束条件里,这样求解器就不会遗漏一些对偶变量了
- 添加切平面条件时,原问题目标函数中的常数值,以及上一个 stage 的求解变量不能遗忘
3. 报童模型的 SDDP 模型
下面以一个 2 阶段的报童模型为例,构造它的 SDDP 求解模型。首先用最传统的报童模型,缺货能补,不考虑产品的销售价格,也不考虑最后一阶段的剩余库存残值。
报童模型中的参数:
- 单位进货成本 c c c
- 单位库存持有成本 h h h
- 单位延迟补货成本(backorder cost) b b b
- t t t 阶段的订货量 q t q_t qt
- 初始库存量是 I 0 I_0 I0, t t t 阶段期末库存量 I t I_t It
如果是单阶段报童模型的话,最优订货点为:
q ∗ = b − c b + h q^\ast=\frac{b-c}{b+h} q∗=b+hb−c
下面构建 2 个阶段(3 个stage)的SDDP模型
k+1 次迭代时,stage 1 的模型:
z = min q 1 k + 1 , θ 2 k + 1 c 1 q 1 k + 1 + θ 2 k + 1 s.t. θ 2 k + 1 ≥ E [ F 2 ( q 1 k , d 1 s ) ] m = 1 , 2 , … , k θ 2 k + 1 ≥ θ ‾ q 1 k + 1 ≥ 0 \begin{aligned} z=&\min_{q_1^{k+1}, \theta_2^{k+1}}~~ &&c_1q_1^{k+1}+\theta_2^{k+1}\\ &\text{s.t.}&& \theta_2^{k+1}\geq \mathbb{E}[\mathcal{F}_2(q^k_{1}, d_1^s)]\qquad m=1,2,\dots,k\\ & &&\theta^{k+1}_2\geq \underline{\theta}\\ & &&q_1^{k+1}\geq 0 \end{aligned} z=q1k+1,θ2k+1min s.t.c1q1k+1+θ2k+1θ2k+1≥E[F2(q1k,d1s)]m=1,2,…,kθ2k+1≥θq1k+1≥0
(stage 1 的模型不用写它的对偶问题了,而接下来的 stage 模型都要涉及到约束条件的对偶)
上面的公式中,用 F \mathcal{F} F 表示子问题对偶问题的目标函数值,这一项是切平面条件。
k+1 次迭代时,stage 2 的模型:
Q
2
(
q
1
k
+
1
,
s
,
d
1
s
)
=
min
I
1
k
+
1
,
B
1
k
+
1
,
q
2
k
+
1
h
I
1
k
+
1
+
b
B
1
k
+
1
+
c
q
2
k
+
1
+
θ
3
k
+
1
s.t.
I
1
k
+
1
−
B
1
k
+
1
=
I
0
+
q
1
k
+
1
−
d
1
s
[
π
2
,
1
s
]
θ
2
k
+
1
≥
E
[
F
3
(
q
2
k
,
d
2
s
)
]
m
=
1
,
2
,
…
,
k
[
π
2
,
2
s
]
θ
2
k
+
1
≥
θ
‾
I
1
k
+
1
≥
0
,
B
1
k
+
1
≥
0
,
q
2
k
+
1
≥
0
\begin{aligned} \mathcal{Q}_2(q^{k+1,s}_1, d_1^s)=&\min_{I_1^{k+1}, B_1^{k+1},q_2^{k+1}}~~ &&hI_1^{k+1}+bB_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\qquad [\pi_{2,1}^s]\\ &&&\theta_2^{k+1}\geq \mathbb{E}[\mathcal{F}_3(q^k_{2}, d_2^s)]\qquad m=1,2,\dots,k\qquad [\pi_{2,2}^s]\\ &&& \theta^{k+1}_2\geq \underline{\theta}\\ & &&I_1^{k+1}\geq 0,B_1^{k+1}\geq 0,q_2^{k+1}\geq 0 \end{aligned}
Q2(q1k+1,s,d1s)=I1k+1,B1k+1,q2k+1min s.t.hI1k+1+bB1k+1+cq2k+1+θ3k+1I1k+1−B1k+1=I0+q1k+1−d1s[π2,1s]θ2k+1≥E[F3(q2k,d2s)]m=1,2,…,k[π2,2s]θ2k+1≥θI1k+1≥0,B1k+1≥0,q2k+1≥0
可以据此写出它对偶问题的目标函数,即
F
2
=
π
2
,
1
s
(
I
0
+
q
1
k
+
1
−
d
1
s
)
−
∑
π
2
,
2
s
π
3
s
d
2
s
\mathcal{F}_2=\pi_{2,1}^s(I_0+q_1^{k+1}-d_1^s)-\sum\pi_{2,2}^s\pi_3^sd_2^s
F2=π2,1s(I0+q1k+1−d1s)−∑π2,2sπ3sd2s
d 1 s d_1^s d1s 表示在情景 s s s 的第一阶段需求值,需求注意的是:求解得到的变量 I 1 k + 1 I_1^{k+1} I1k+1, B 1 k + 1 B_1^{k+1} B1k+1 , q 2 k + 1 q_2^{k+1} q2k+1 其实也是跟情景 s s s 相关的。
- 上面表达式的第二项不能丢,否则出错。因为之前迭代时添加的切平面条件也有对应的对偶值
k+1 次迭代时,stage 3 的模型:
Q 3 ( q 2 k + 1 , s , d 1 s ) = min I 2 k + 1 , B 2 k + 1 h I 2 k + 1 + b B 2 k + 1 s.t. I 2 k + 1 − B 2 k + 1 = I 1 k + 1 + q 2 k + 1 − d 2 s [ π 3 s ] I 2 k + 1 ≥ 0 , B 2 k + 1 ≥ 0 \begin{aligned} \mathcal{Q}_3(q^{k+1,s}_2, d_1^s)=&\min_{I_2^{k+1}, B_2^{k+1}}~~ &&hI_2^{k+1}+bB_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\qquad [\pi_3^s]\\ & &&I_2^{k+1}\geq 0,B_2^{k+1}\geq 0 \end{aligned} Q3(q2k+1,s,d1s)=I2k+1,B2k+1min s.t.hI2k+1+bB2k+1I2k+1−B2k+1=I1k+1+q2k+1−d2s[π3s]I2k+1≥0,B2k+1≥0
可以据此写出它对偶问题的目标函数,即
F
3
=
π
3
s
(
I
1
k
+
1
+
q
2
k
+
1
−
d
2
s
)
\mathcal{F}_3=\pi_3^s(I_1^{k+1}+q_2^{k+1}-d_2^s)
F3=π3s(I1k+1+q2k+1−d2s)
4. 求解单阶段报童模型的一个例子
我用这个方法试了试求解单阶段的报童,求解结果不错。但应用到多阶段时编程要复杂多了,有几点需要注意。
- 约束条件要都是线性表达式,含有 max 或 min 的表达式要处理下
- 若有整数求解变量,更加复杂,似乎要用到 sddip 方法
单阶段报童模型其实也是按照一个两阶段模型来求解:第一阶段的决策变量是订货量(需求已知前),而第二阶段的决策变量是销售量(需求已知后)。
当然,报童模型的最终求解变量其实是订货量。
k+1 次迭代时,stage 1 的模型:
z = min q 1 k + 1 , θ 2 k + 1 c 1 q 1 k + 1 + θ 2 k + 1 s.t. θ 2 k + 1 ≥ E [ Q 2 ( q 1 k , d 1 s ) ] − π 2 k + 1 , s ( q 1 k + 1 − q 1 k ) m = 1 , 2 , … , k θ 2 k + 1 ≥ θ ‾ q 1 k + 1 ≥ 0 \begin{aligned} z=&\min_{q_1^{k+1}, \theta_2^{k+1}}~~ &&c_1q_1^{k+1}+\theta_2^{k+1}\\ &\text{s.t.}&& \theta_2^{k+1}\geq \mathbb{E}[\mathcal{Q}_2(q^k_{1}, d_1^s)]-\pi_2^{k+1,s}(q^{k+1}_1-q^k_1)\qquad m=1,2,\dots,k\\ & &&\theta^{k+1}_2\geq \underline{\theta}\\ & &&q_1^{k+1}\geq 0 \end{aligned} z=q1k+1,θ2k+1min s.t.c1q1k+1+θ2k+1θ2k+1≥E[Q2(q1k,d1s)]−π2k+1,s(q1k+1−q1k)m=1,2,…,kθ2k+1≥θq1k+1≥0
k+1 次迭代时,stage 2 的模型(此时的决策变量是销量):
Q
2
(
q
1
k
+
1
,
s
,
d
1
s
)
=
min
y
1
k
+
1
−
p
y
1
k
+
1
−
γ
(
q
1
k
+
1
,
s
−
y
1
k
+
1
,
s
)
s.t.
−
y
1
k
+
1
≥
−
q
1
k
+
1
,
s
−
y
1
k
+
1
≥
−
d
2
s
y
1
k
+
1
≥
0
\begin{aligned} \mathcal{Q}_2(q^{k+1,s}_1, d_1^s)=&\min_{y_1^{k+1}}~~ &&-py_1^{k+1}-\gamma (q_1^{k+1,s}-y_1^{k+1,s})\\ &\text{s.t.} &&-y_1^{k+1}\geq -q_1^{k+1,s}\\ &&&-y_1^{k+1}\geq -d_2^s\\ & &&y_1^{k+1}\geq 0 \end{aligned}
Q2(q1k+1,s,d1s)=y1k+1min s.t.−py1k+1−γ(q1k+1,s−y1k+1,s)−y1k+1≥−q1k+1,s−y1k+1≥−d2sy1k+1≥0
也可以写成(此时的决策变量是库存量):
Q
2
(
q
1
k
+
1
,
s
,
d
1
s
)
=
min
I
1
k
+
1
−
γ
I
1
k
+
1
−
p
(
q
1
k
+
1
,
s
−
I
1
k
+
1
,
s
)
s.t.
I
1
k
+
1
≥
q
1
k
+
1
,
s
−
d
1
s
−
I
1
k
+
1
≥
−
q
1
k
+
1
,
s
I
1
k
+
1
≥
0
\begin{aligned} \mathcal{Q}_2(q^{k+1,s}_1, d_1^s)=&\min_{I_1^{k+1}}~~ &&-\gamma I_1^{k+1}-p (q_1^{k+1,s}-I_1^{k+1,s})\\ &\text{s.t.} &&I_1^{k+1}\geq q_1^{k+1,s}-d_1^s\\ &&&-I_1^{k+1}\geq -q_1^{k+1,s} \\ &&&I_1^{k+1}\geq 0 \\ \end{aligned}
Q2(q1k+1,s,d1s)=I1k+1min s.t.−γI1k+1−p(q1k+1,s−I1k+1,s)I1k+1≥q1k+1,s−d1s−I1k+1≥−q1k+1,sI1k+1≥0
之前漏写了一个约束条件,一直求不对割平面约束条件,上面模型中的约束条件一个都不能少。
它们的对偶模型分别为:
max
π
1
k
+
1
,
π
2
k
+
1
−
π
1
k
+
1
q
1
k
+
1
−
π
2
k
+
1
d
1
s
−
γ
q
1
k
+
1
s.t.
−
π
1
k
+
1
−
π
2
k
+
1
≤
−
p
+
γ
−
π
1
k
+
1
≥
0
,
−
π
2
k
+
1
≥
0
\begin{aligned} &\max_{\pi_1^{k+1},\pi_2^{k+1}}~~ &&-\pi_1^{k+1}q_1^{k+1}-\pi^{k+1}_2d_1^s-\gamma q_1^{k+1}\\ &\text{s.t.} &&-\pi_1^{k+1}-\pi_2^{k+1}\leq -p+\gamma\\ &&&-\pi_1^{k+1}\geq0, -\pi^{k+1}_2\geq 0 \end{aligned}
π1k+1,π2k+1max s.t.−π1k+1q1k+1−π2k+1d1s−γq1k+1−π1k+1−π2k+1≤−p+γ−π1k+1≥0,−π2k+1≥0
以及:
max
π
1
k
+
1
,
π
2
k
+
1
π
1
k
+
1
(
q
1
k
+
1
−
d
1
s
)
−
π
2
k
+
1
q
1
k
+
1
−
p
q
1
k
+
1
s.t.
π
1
k
+
1
−
π
2
k
+
1
≤
p
−
γ
π
1
k
+
1
≥
0
,
π
2
k
+
1
≥
0
\begin{aligned} &\max_{\pi_1^{k+1},\pi_2^{k+1}}~~ &&\pi_1^{k+1}(q_1^{k+1}-d_1^s)-\pi^{k+1}_2q_1^{k+1}-p q_1^{k+1}\\ &\text{s.t.} &&\pi_1^{k+1}-\pi_2^{k+1}\leq p-\gamma\\ &&&\pi_1^{k+1}\geq0, \pi^{k+1}_2\geq 0 \end{aligned}
π1k+1,π2k+1max s.t.π1k+1(q1k+1−d1s)−π2k+1q1k+1−pq1k+1π1k+1−π2k+1≤p−γπ1k+1≥0,π2k+1≥0
几点总结:
- 对偶变量 π \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)