鲁棒优化| C&CG算法求解两阶段鲁棒优化：全网最完整、最 ...

rustum · 发表于 2022-7-20 09:52

作者：刘兴禄，清华大学，清华伯克利深圳学院博士在读，E-mail: hsinglul@163.com,
初次完稿日期：2021.09.28
审校：张瑞三，四川大学，硕士在读，E-mail:zrssum3@stu.scu.edu.cn，知乎ID：MunSum3
需要加入读者群的读者，可以关注公众号或者联系知乎编辑,感谢您的支持！
<hr/>

特别感谢Zeng Bo老师非常耐心地解答我的疑惑，您的解答对我来讲非常重要。
本文非常详细、完整地介绍了求解Two-stage Robust Optimization Model的一种精确算法：Column and Constraint Generation Algorithm。该算法由Zeng Bo, Zhao Long于2013年首次提出，相关论文发表在期刊《Operations Research Letters》上，论文题目为《Solving two-stage robust optimization problems using a column-and-constraint generation method》。
本文对该论文进行了非常细致的解读，包括：

对鲁棒优化问题进行了分类；
对两阶段鲁棒优化模型进行了详细解析；
对Bender-dual算法进行了讲解；
对Column-and-constraint generation method进行了非常详细的解读、完整的推导和完整的解释；
对论文中的案例进行了相当详细的推导；
对Column-and-constraint generation method进行了完整的代码复现（Python调用Gurobi实现）。

图1-CCG经典文章

两阶段鲁棒优化问题（Two-stage Robust Optimization）

鲁棒优化是一类考虑参数不确定的数学规划问题，是运筹学中比较高阶的方法。近几年，关于鲁棒优化的研究越来越火热。常见的鲁棒优化问题包括基本的鲁棒优化、多阶段鲁棒优化、分布式鲁棒优化等。
鲁棒优化旨在处理优化问题中参数不确定的挑战，致力于优化最坏情况(worst case)下的解，或者最坏参数分布下的解(对应分布式鲁棒优化)，从而使最终得到的解具有非常强的抵抗不确定性或者风险的能力。
处理参数不确定的优化问题，还有一类方法，叫作随机规划，关于二者的区别，见往期推文。
两阶段鲁棒优化问题，是多阶段鲁棒优化的一个特例(multi-stage robust optimization)。在两阶段鲁棒优化问题中，一般包含两个不同层级的决策变量，又称为第一阶段的决策(first-stage decision)和第二阶段的决策(second-stage decision)。

战略层决策和战术层决策；
战术层决策和操作层决策；
...

一个典型例子，就是选址-配送规划问题。选址决策属于战略层决策，配送问题，可以视为战术层或者作业层决策。
两阶段鲁棒优化问题描述如下：

第一阶段的决策相关的参数是确定的，并且第一阶段的决策需要首先做出；
第二阶段的决策的相关参数有一些是不确定的。第二阶段的决策需要在第一阶段的决策确定之后，等到第二阶段参数的揭示后，再做出相应的第二阶段决策。
两阶段鲁棒优化的目标为：对两阶段决策进行联合优化，同时考虑第二阶段参数的不确定性。也就是优化在第二阶段参数取到最坏情况下的两个阶段决策对应的总的目标值。

下面我们用标准的数学语言来描述两阶段鲁棒优化问题。两阶段鲁棒优化问题的数学模型如下所示(参考自文献[1])

其中:

为第一阶段的决策变量;
为第二阶段的决策变量;
$\mathbf{c}$ 为第一阶段的决策变量相关的参数，是确定的;
为第二阶段的决策变量相关的参数，是不确定的，且其不确定参数的取值范围由不确定集刻画；具体的例子见往期推文；
目标函数中，第一项 $\min_{\mathbf{y}} \,\mathbf{c^Ty}$ 为第一阶段决策对应的目标函数。该部分决策需要首先确定；
目标函数中，第二项 $\max_{u \in \mathcal{U}} \min_{ \mathbf{x}\in {\mathbf{F}(\mathbf{y}, u )}} \mathbf{b}^T \mathbf{x}$ 为第二阶段决策对应的目标函数。该部分决策需要在第一阶段决策确定之后再做出；该部分目标函数旨在找到worst case;也就是在最坏情况下，第二阶段的目标值；
整个目标函数，就是在优化在最坏情况下的总目标函数，使得最坏的情况最好，使得解具有非常好的鲁棒性。
约束(2)只跟第一阶段的决策有关；
约束(3)中跟第一阶段，第二阶段的决策都有关，且包含不确定参数。

以上就是两阶段鲁棒优化的一个简要介绍。目前还比较晦涩的数学符号，我们将在下文进行更通俗的解读。

注意，如果的取值可能性非常少，我们就可以通过穷举的取值，显式地将每个可能的对应的约束都加进模型，然后将两阶段鲁棒优化模型等价为一个可直接求解的一阶段数学规划模型，从而可以达到直接求解上述Two-stage robust optimization模型的目的。因为最坏情况，一定是对应着某一种的取值，而我们已经穷举了所有的取值可能，且将其显式地写进了模型，因此原两阶段鲁棒优化模型可以直接被解决。但是通常，的取值非常多（组合数级别或者无穷多），穷举往往是低效甚至不可行的，因此，需要利用一些其他方法来求解Two-stage robust optimization model。本文将要介绍的Column and Constraint Generation算法就是其中一种（此外，文献[1]中提到的Benders-dual algorithm也是一种）。

两阶段鲁棒优化问题的详细解读

参考文献[1]中提到：

Two-stage RO 的计算非常困难，即使是一个简单的Two-stage RO问题也可能是NP-Hard的。为了克服计算负担，目前主要流行两种求解策略。第一种是使用 approximation algorithms，它假设第二阶段的决策是不确定性的简单函数，如仿射函数。第二类算法是利用Benders分解法求出精确解，即利用第二阶段决策问题的对偶解逐步构造第一阶段决策的价值函数。因此，我们称它们为Benders-dual cutting plane algorithms算法。

文献[1]中关于两阶段鲁棒优化问题模型的描述如下：

图2-文献[1]中关于两阶段鲁棒优化问题模型的描述

该论文探讨的问题，是基于如下的假设：

决策变量类型假设：first-stage and second-stage are linear optimization
参数不确定性(uncertainty): is either a finite discrete set or a polyhedron. 换句话说，改论文考虑的不确定集为是一个有限的离散集合或者一个多面体集合，比如
a finite discrete set：不确定参数的可能取值为一个有限的可选集，例如： $u \in \mathcal{U} = { \text{scenario}_1, \text{scenario}_2, \text{scenario}_3}$ .
a polyhedron: 不确定参数的可能取值的范围为一个多面体，例如 $u \in \mathcal{U} = { 1 \leqslant u_1 \leqslant 2, 3 \leqslant u_2 \leqslant 5}$ ，该不确定集为一个矩形，学术上一般称其为盒子不确定集(box uncertainty set)。

根据上图，两阶段鲁棒优化中，第一阶段的决策需要首先做出，其优化方向为。
但是第二阶段的决策，却是一个双层问题(bi-level)。具体来讲，目标函数第二项的符号为 $\\ \max_{u \in \mathcal{U}} \min_{ \mathbf{x}\in {\mathbf{F}(\mathbf{y}, u )}} \mathbf{b}^T \mathbf{x}$
这个符号需要详细解释一下：

$\mathbf{x}\in {\mathbf{F}(\mathbf{y}, u )}$ 的意思: 表示given 不确定参数的一个取值，以及given first-stage决策变量的取值，我们决策，使得 $\mathbf{b}^T \mathbf{x}$ 最小。并且满足如下约束 $\\ \begin{aligned} \mathbf{F}(\mathbf{y}, u ) = \{ \mathbf{x} \in \mathbf{S}_{\mathbf{x}} : \mathbf{Gx \geqslant h - Ey - M}u \} \end{aligned}$
with 。

其实就是说，和分别是 $n\times 1$ 维和 $m\times 1$ 维的非负连续决策变量而已。
我们将上图中的模型展开，写成一个比较容易看明白的形式：

只不过，论文截图中 $\mathbf{y \in S_y}$ ，就等价于是下面模型的可行解 $\\ \begin{aligned} \min_{\mathbf{y}} \quad & \mathbf{c^Ty} +\max_{u \in \mathcal{U}} \min_{ \mathbf{x}\in {\mathbf{F}(\mathbf{y}, u )}} \mathbf{b}^T \mathbf{x} \\ s.t.\quad & \mathbf{Ay \geqslant d}, \\ & \mathbf{S_y} \subseteq \mathbb{R}_{+}^{n} \end{aligned}$
其中 $u, \mathbf{x}$ 是一个给定的值。

上述Two-stage Robust Optimization模型，一般是无法直接调用优化求解器进行求解的，需要设计对应的算法进行求解。下面就来介绍求解该模型的算法。
Two-stage RO和Benders-dual cutting plane method

类似于使用标准的Benders Decomposition算法求解确定性问题，Two-stage RO也是可以使用Benders Decomposition算法来求解的。文献中将该算法称之为Benders-dual cutting plane method。
本文仅对Benders-dual cutting plane method做简要介绍。小编其实早就已经写好了Benders-dual cutting plane method的完整笔记和代码，只不过想留着后面再分享，尽请期待！
Benders-dual cutting plane method求解Two-stage RO问题可以类比Benders Decomposition求解确定性问题的情形。我们可以将Two-stage RO分解为：

一个主问题(Master Problem)，主要包含第一阶段的决策及只跟有关的约束、子问题返回的割，评估第二阶段目标函数取值的辅助变量（类似于Benders Decomposition算法中的 $q$ ）;
一个子问题(Subproblem)，主要包含第二阶段的决策和不确定变量，旨在给出第二阶段目标函数值的一个界限（类似于标准的Benders Decomposition算法）。

Benders-dual cutting plane method中的初始主问题和Benders Decomposition算法类似，为

$\\ \begin{aligned} \min_{\mathbf{y}} \quad & \mathbf{c^Ty} +\eta && (1) \\ s.t.\quad & \mathbf{Ay \geqslant d}, && (2) \\ & \mathbf{S_y} \subseteq \mathbb{R}_{+}^{n}, && (3) \\ &\eta \geqslant 0. && (4) \end{aligned}$
初始主问题并没有任何Benders cut的加入。

下面我们来大致推导Benders-dual cutting plane method中的割平面。首先来看second-stage的部分，也就是

$\\ \begin{aligned} \max_{u \in \mathcal{U}} \min_{ \mathbf{x}\in {\mathbf{F}(\mathbf{y}, u )}} \quad & \mathbf{b}^T \mathbf{x} \\ s.t.\quad& \mathbf{Gx \geqslant h - Ey - M}u , \qquad u \in \mathcal{U} \\ & \mathbf{S_x} \subseteq \mathbb{R}_{+}^{m} \end{aligned}$

这部分原文是这么描述的

Consider the case where the second-stage decision problem is a linear programming (LP) problem in . We first take the relatively complete recourse assumption that this LP is feasible for any given and . Let be its dual variables. Then, we obtain its dual problem, which is a maximization problem and can be merged with the maximization over . As a result, we have the following problem, which yields the subproblem in the Benders-dual method.

由于上式中，外层是 $\max_{u \in \mathcal{U}}$ ，外层的决策变量是不确定变量，因此对于内层， $\min_{ \mathbf{x}\in {\mathbf{F}(\mathbf{y}, u )}}$ 而言，和都是已知(given)的，也就是，inner level的模型其实是 given , given ,

$\\ \begin{aligned} \min_{ \mathbf{x}\in {\mathbf{F}(\mathbf{\bar{y}}, \bar{u} )}} \quad & \mathbf{b}^T \mathbf{x} \\ s.t.\quad& \mathbf{Gx \geqslant h - E\bar{y} - M}\bar{u} , \qquad \rightarrow \mathbf{\pi} \\ & \mathbf{S_x} \subseteq \mathbb{R}_{+}^{m} \end{aligned}$
其中是求解first-stage的模型得到的解，而 $\bar{u}$ 是在second-stage的outer-level中已经fix的不确定的量的取值。另外， 是一个列向量，并不是标量。因此， $h - E\bar{y} - \mathbf{M}\bar{u}$ 是已知常数列向量，是约束的右端常数项。

由于上述inner level的模型是线性规划，所以我们直接写出上述模型的对偶问题，为

$\\ \begin{aligned} \max_{ \pi} \quad & (\mathbf{h - E\bar{y} - M}\bar{u} )^T \pi \\ s.t.\quad& \mathbf{G^T \pi \leqslant b} , \\ & \pi \geqslant \mathbf{0}. \end{aligned}$
注意，是一个列向量，并不是标量。以及 $\mathbf{0}$ 是一个元素全部为0的列向量。

如果我们把second-stage的outer-level和inner level的模型写在一起，就等价于。于是我们就将第二阶段bi-level（双层）优化模型转化成了一个单层优化模型，也就是把

$\\ \begin{aligned} \max_{u} \max_{ \pi} \quad & (\mathbf{h - E\bar{y} - M}u )^T \pi \\ s.t.\quad& \mathbf{G^T \pi \leqslant b} , \\ & \pi \geqslant \mathbf{0} \\ & u \in \mathcal{U}. \end{aligned}$

等价为

$\\ \begin{aligned} \max_{ \pi, u} \quad & (\mathbf{h - E\bar{y} - M}u )^T \pi \\ s.t.\quad& \mathbf{G^T \pi \leqslant b} , \\ & \pi \geqslant \mathbf{0} \\ & u \in \mathcal{U}. \end{aligned}$

注意，由于是将等价成了，因此outer-level的不确定变量也就纳入到等价后的一阶段模型中来了。其中是从第一阶段得到的fix的值。

由于上述模型中，已经被固定，因此就为最优解中第二阶段的目标函数提供了一个下界。

为什么就为最优解中第二阶段的目标函数提供了一个下界呢？
因为，如果没有被固定，则 $\max_{ \pi, u} \, (\mathbf{h - Ey - M}u )^T \pi \geqslant \max_{ \pi, u} \, (\mathbf{h - E\bar{y} - M}u )^T \pi$ 一定成立，因此就为第二阶段的目标函数提供了一个下界。

我们将上述模型简略地写为

$\\ \begin{aligned} \mathbf{SP_1}: \mathcal{Q}(\mathbf{y}) = \max_{ \pi, u} \{ (\mathbf{h - E\bar{y} - M}u )^T \pi : \mathbf{G^T \pi \leqslant b} ,u \in \mathcal{U}, \pi \geqslant \mathbf{0}.\} \end{aligned}$

这相当于把上面的数学规划模型，写成了一个函数形式，实际上是等价的，只是另一个简便，节省地方的写法。这个部分，也被称为是recourse problem，可以理解为，第一阶段决策确定了，第二阶段等到不确定参数揭示了，我们做第二阶段的决策，来对整个问题进行追索，秋后算账。
写好了子问题的模型，那接下来的重要工作当然是：求解。

如何求解呢？：
注意，是一个 bilinear optimization problem （因为 $u\cdot \pi$ 是二次项，两个都是决策变量），根据文献[1]中的参考文献，可知：

可以用一些启发式算法求解 ;这种方法得到的不一定是的最优解，因此算法收敛速度也许会慢一些；
另外，如果具有特殊的结构，可以构造精确性算法进行精确求解; 比如，如果是0-1变量，可以进行等价线性化，变成MIP进行求解；
当然，根据笔者的经验，还可以使用Gurobi，COPT等求解器来求解 ;
最后，还可以使用文献[1]介绍的方法：使用KKT条件来求解。本推文也是使用了KKT条件来求解的，见下文的超级详细的介绍；

综上，大概可以用四种可以来求解：
1. 启发式算法
2. 具有特殊的结构，可以构造精确性算法进行精确求解;
3. 使用Gurobi，COPT等求解器来求解；
4. 使用KKT条件来求解。
这部分原文是这么描述的

Note that the resulting problem in (2) is a bilinear optimization problem. Several solution strategies have been developed, either in a heuristic fashion or for instances with specially-structured .

求解了之后，又能有什么用呢？
答案是：为主问题提供割平面，从而改进全局界限，改进当前解。
假设在Benders-dual算法中，对于给定的第步迭代中第一阶段（Master Problem）得到的 $\mathbf{y}_k^{}$ ，我们都能得到一个相应 $\mathcal{Q}(\mathbf{y}_k^{})$ 的第二阶段的解 $(u_k^*\pi_k^*)$ ，我们就可以生成一个下面形式的割平面：

上面已经介绍过， $\max_{ \pi, u} \, (\mathbf{h - E\bar{y} - M}u )^T \pi$ 为最优解中第二阶段的目标函数的取值提供了一个下界

因此，割平面的形式如下：

$\begin{aligned} \eta \geqslant \mathbf{(h - Ey - M}u_k^* )^T \pi_k^* \end{aligned}$ 其中，表示最优解中第二阶段的目标函数的取值。这个割平面，其实就可以看做是Benders decomposition中的optimality cut。（比较直观的理解请参考我们往期讲解Benders decomposition的推文中最后的Cut 5和 Cut 7）

这个割平面可以被加入到主问题中，主问题就变化为

$\\ \begin{aligned} \mathbf{MP_1}: \quad &\min_{\mathbf{y}, \eta}\quad \mathbf{c^T y + \eta} \\ s.t.\quad& \mathbf{Ay \geqslant d} \\ & \eta \geqslant \mathbf{(h - Ey - M}u_l^* )^T \pi_l^*, \qquad \forall l \leqslant k \\ & \mathbf{y \in S_y}, \\ & \eta \in \mathbb{R}. \end{aligned}$
注意，第2条约束中， $\forall \ l \leqslant k$ 指的是在前次的迭代中，每一次产生的割平面都被加入到了 中。并且相当于一个辅助变量，是为了识别worst-case下的second-stage产生的成本的(也叫作recourse)。

上述模型，我们可以计算到目前为止得到的当前最优解 $(\mathbf{y_{k+1}^{}}, \eta_{k+1}^{})$ ，这也就是第 $k+1$ 次迭代的时候二者的取值。

注意到，随着割平面 $\eta \geqslant \mathbf{(h - Ey - M}u_l^* )^T \pi_l^* \,\,\, \forall l \leqslant k$ 的不断添加，中的值会不断增加(至少是不减)，因此的目标函数会不断上升。下面会讲到，的目标函数是全局最优解的一个下界。

与传统的Benders Decomposition算法中相同，Benders-dual算法迭代过程中，也可以计算任意一步迭代过程中，全局的上界和下界的变化，以此来判断当前解的质量，以及是否达到了最优解。
文献[1]中原文对Benders-dual method下的主问题的表述如下：

图3-原文对Benders-dual method下的主问题的表述

注意：文中提到了一些非常关键的信息，即Benders-dual method迭代过程中，全局上界和下界的计算。根据原文：

为原始问题目标值 $\text{Obj} = {\min_{\mathbf{y}} \mathbf{c^Ty} +\max_{u \in \mathcal{U}} \min_{ \mathbf{x}\in {\mathbf{F}(\mathbf{y}, u )}} \mathbf{b}^T \mathbf{x} | \mathbf{Ay} \geqslant \mathbf{d}, \mathbf{y} \in \mathbf{S}{\mathbf{y}}}$ （即原始问题）提供了一个上界；
而 $\mathbf{MP_1}: \,\, \mathbf{c^Ty_{k+1}^{*}} + \eta_{k+1}^{*}$ 为原始问题提供了一个下界。
在Benders-dual method算法迭代的过程中，上界和下界会得到不断改进。并且，循环地求解主问题，得到固定的，并将传给子问题，更新子问题，然后求解，根据的解，构造割平面 $\eta \geqslant \mathbf{(h - Ey - M}u_k^* )^T \pi_k^*$ ，并且将其添加到中。
同时，更新全局的上界和下界。循环这个过程，最终上下界就可以在有限步之内收敛。

总结一句，Benders-dual算法迭代过程中，任意一步迭代中的全局的上界和下界可以用下面的公式计算得出：

这个很好理解：

$LB = \mathbf{c^Ty_{k+1}^{*}} + \eta_{k+1}^{*}$ ，是只考虑了一部分随机场景，松弛了一部分约束，因此对其求，一定是更小(或者等于最小)，因此它提供了全局最优解的一个下界；直观理解就是少考虑了约束，又是min问题，因此是个下界；
$UB=\mathbf{c^Ty_k^* + \mathcal{Q}(y_k^*)}$ ，是相当于将原问题拆分成两个相互独立的问题进行分别求解，然后将各自的目标函数相加，分开优化的可行解，一定也是联合优化的可行解；而分开优化的解，对于联合优化而言，只是一个可行解而已，并不一定是最优解，联合优化的解，一定好于或者等于分开优化的解。因此为Two-stage RO原模型的最优解提供了一个上界。
当算法迭代过程中，达到 $UB=LB$ ，则算法终止，得到了最优解。

注意到，上述算法描述中， $\pi_k^*$ 和 $u_k^*$ 都是在他们各自的可行域中的极点(或者离散的点)。这样一来，Benders-dual算法的复杂度如下：

图4-Benders-dual算法的复杂度

The column-and-constraint generation algorithm

上一节我们介绍了Benders-dual method，本节我们来介绍column-and-constraint generation algorithm。
文献[1]提出的Column-and-constraint generation algorithm可以看做是一种Benders Decomposition的拓展，但是与Benders Decomposition又有显著的不同。具体的不同点，下文有详细介绍。
我们首先来描述一下Column-and-constraint generation algorithm的基本思想，然后再详细阐述该算法的详细原理。
Column-and-constraint generation algorithm的基本思想如下：

首先，只考虑第一阶段的决策变量和只包含第一阶段决策变量的约束，这个版本可以看做是整个问题的一个松弛版本。然后我们想办法将松弛了的部分的影响逐渐添加回来。我们来看第二阶段，第二阶段的目标旨在找到worst case。假设我们直接可以找到worst case，我们只需要直接将worst case下的参数，及其对应的约束带入到主问题中即可等价地解决Two-stage RO问题。但是直接识别worst case一般是不可能的。不过，我们可以通过固定第一阶段的决策变量，然后求解其对应的第二阶段的子问题，并且找到一个非常可能是worst case的 ，之后将这个非常可能是worst case的 对应的场景下的决策变量和约束全部加回到主问题中。循环这个操作，主问题中考虑的可能是最坏场景的场景就会越来越多，解就会被逐渐改进，最终达到最优解。

由于主问题变量和约束数量增多，下界会上升，而主问题提供给子问题的也在变化，会使得全局上界得到改善(也就是下降)。持续迭代这个过程，上界和下界不断得到改进，直到算法收敛。

这就是Column-and-constraint generation algorithm的基本思想。总结一句：

首先将不确定变量松弛了，然后通过求解子问题识别最有可能是worst case的，并将其对应的约束和决策变量加回到主问题中，将松弛掉的部分的影响再以秋后算账的方式(也就是recourse)加回来。循环该过程，直到算法收敛为止。

原文中的描述如下：

In this section, we present another cutting plane procedure to solve two-stage RO problems. Because the generated cutting planes are defined by a set of created recourse decision variables in the forms of constraints of the recourse problem, the whole procedure is a column-and-constraint generation (C&CG) procedure. To make our exposition simple, we first mention an observation when is a finite discrete set. Let $\ \ \mathcal{U} = {u_1, \cdots , u_r }$ and ${\mathbf{x^1, \cdots , x^r} }$ be the corresponding recourse decision variables. Then, the two-stage RO in (1) in Section 2 can be reformulated as the following:(原文)

注意，这里，每一个被识别的(也就是有可能是worst case的)不确定的场景，都会对应该场景下的一系列的决策变量，是一个向量，也就是 $\mathbf{x^r}$ ；类似于随机规划中一个场景下的一套决策。
举个例子，假设只有4个变量，分别为 $x_1, x_2, x_3, x_4$ ：

比如是第1种场景，假设求解完子问题后，得到相应的 $\mathbf{x^1} = [2, 4, 1, 4]$ ，也就是 $x_1^1 = 2, x_2^1 = 4, x_3^1 = 1, x_4^1 = 4$ ；
比如是第2种场景，假设求解完子问题后，得到相应的 $\mathbf{x^2} = [1, 1, 2, 2]$ ，也就是 $x_1^2 = 1, x_2^2 = 1, x_3^2= 2, x_4^2 = 2$ ；
$\cdots$

基于上面的描述，则原来的two-stage RO模型（为了查看方便，我们再写一遍），

$\\ \begin{aligned} \min_{\mathbf{y}} \quad & \mathbf{c^Ty} +\max_{u \in \mathcal{U}} \min_{ \mathbf{x}\in {\mathbf{F}(\mathbf{y}, u )}} \mathbf{b}^T \mathbf{x} \\ s.t.\quad & \mathbf{Ay \geqslant d}, \qquad \mathbf{y \in S_y} \end{aligned}$
where
$\\ \mathbf{F}(\mathbf{y}, u ) = \{ \mathbf{x} \in \mathbf{S}_{\mathbf{x}} : \mathbf{Gx \geqslant h - Ey - M}u \}$
with
$\\ \mathbf{S_y} \subseteq \mathbb{R}{+}^{n}$ and $\mathbf{S_x} \subseteq \mathbb{R}{+}^{m}$ .

可以被reformulation为下面的形式

$\\ \begin{aligned} \min_{\mathbf{y}} \quad &\mathbf{c^T y} + \eta && && (1) \\ s.t. \quad & \mathbf{Ay \geqslant d},&& && (2) \\ & \eta \geqslant \mathbf{b^T x}^l, \quad && l=1, \cdots, r, && (3) \\ & \mathbf{Ey + G} \mathbf{x}^l \geqslant \mathbf{h - M}u_l, \quad && l = 1, \cdots, r , && (4) \\ & \mathbf{y \in S_y}, && && (5) \\ & \mathbf{x}^l \in \mathbf{S_x}, && l=1, \cdots, r. && (6) \end{aligned}$

其中， $l=1, \cdots, r$ 表示 $r$ 个被子问题识别的场景。上述reformulated的模型中， 实际上已经是一个已知的值了，对应一种被识别的场景。
上面的reformulation的形式中，约束(3), (4),(6)就是reformulation引入的约束，也就是把两阶段鲁棒模型，转化为一阶段的优化问题了。并且 ,  其实就是第二阶段的目标函数值。约束(3)就是说，  是第 $l$ 个极点对应的第二阶段的目标函数，而  是全局最优解中，第二阶段的目标值，因此  一定 $\geqslant\mathbf{b^T x}^l$ 。同时由于目标函数优化方向为  ，即： $\min_{\mathbf{y}} \mathbf{c^T y} + \eta$ ，因此  就正好会取到所有  中的最大值，因此就是等价的转化。
同样的，第二阶段问题识别的场景 $u^l$ 对应的约束也需要引入主问题。也就是约束  ，也必须要被加入到主问题中。这是保证在给定的场景(  )下，决策变量需要满足相应的约束条件，需要可行。在这里我们再次强调，reformulation的形式中  已经是给定的参数了，因为我们通过算法，已经将其识别了。
除了不断加割平面 $\eta\geqslant \mathbf{b^T x}^l$ 和约束  外，我们还需要相应地向主问题添加新的决策变量：  。因此，随着算法的迭代，主问题中的决策变量  、以及约束都在不断增多。也就是 $\mathbf{x}^l \in \mathbf{S_x}, \quad l=1, \cdots, r$ 是随着迭代次数的增加，一直在增加的。正是因为主问题中，行和列都在随着迭代的进行而增加，因此该算法才叫column-and-constraint generation (C&CG)。

说明：
① $\mathbf{x}^l \in \mathbf{S_x}, \quad l=1, \cdots, r.$ 是随着算法的迭代不断增加的。这就是column generation；
② 并且，随着迭代的进行，约束也在不断增加，也就是在被不断引入，这就是约束生成，constraint generation。

这部分原文是这么描述的

As a result, solving a two-stage RO problem reduces to solve an equivalent (probably large-scale) mixed integer program. When the uncertainty set is very large or is a polyhedron, developing the equivalent formulation by enumerating all the possible uncertain scenarios in U and deriving its optimal value is not practically feasible. Nevertheless, based on constraints in (7) (约束 (7)也就是 ), it is straightforward that a formulation based on a partial enumeration, i.e., a formulation defined over a subset of , provides a valid relaxation (and, consequently, a lower bound) to the original two-stage RO (or its equivalent formulation). Hence, by expanding a partial enumeration by adding non-trivial scenarios gradually, stronger lower bounds can be expected. With this observation in mind, we were motivated to design a column-andconstraint generation procedure that expands a subset of byidentifying and including significant scenarios, i.e., generating the corresponding recourse decision variables and (7)–(8) on the fly.

上面提到，在C&CG的每一步的迭代中，我们都需要向主问题加入下面的割平面和约束，即约束(3)–(4)：

和

$\mathbf{Ey + G} \mathbf{x}^l \geqslant \mathbf{h - M}u_l, \quad l = 1, \cdots, r$

因此，求解一个两阶段的RO问题可以简化为解决一个等价的（可能是大规模的）一阶段的优化模型。当不确定性集合非常大或者是一个多面体时，通过列举中所有可能的不确定性场景，并得出等价的模型实际上是不可行的。然而，基于（3）中的约束条件（也就是 )，很明显，这是一个部分枚举的公式，即（3）是定义在的一个子集上的公式，它为原始的两阶段RO(或其等价模型)提供了一个有效的松弛(因此，也是一个下界)。所以，通过逐步增加有意义 (或者有可能为worst case)的情景来扩大部分枚举，可以预期会有更紧的下界。基于这些启示，作者设计了列与约束生成的算法，通过识别和包括重要的场景来扩展一个的一个子集(也就是被识别出来的重要的场景)，即在on-the-fly的过程中生成相应的追索决策变量和割平面（3），以及场景对用的约束（4）。

这里需要注意，为什么只需要考虑一部分重要的场景（即有可能为worst case的场景）即可呢？

举一个很简单的例子，假设有 $x_1, x_2$ 两个决策变量以及有2个不确定参数 $u_1, u_2$ ，他两的取值分别在 $u_1 \in {0, 1 }, u_2 \in {0, 2 }$ 。这样一来，就是有4种场景，因为也就是 $[u_1, u_2] = [0, 0]$ 或者 $[0, 2]$ 或者 $[1, 0]$ 或者 $[1, 2]$ 。但是在RO中，实际上只考虑最坏情况，只可能是两个不确定参数都取到了最坏值，因此只需要考虑的就是 $[u_1, u_2] = [1, 2]$ ，其余的场景都可以不考虑。那么在最终的等价模型中，只需要识别出重要的场景即可，其余的场景是不会在最优解中binding的，因为他们不是最坏的情形，所以不需要考虑。

这部分原文是这么描述的

Similar to the Benders-dual method, this column-andconstraint generation procedure is implemented in a mastersubproblem framework. We assume that an oracle can solve the following subproblem in the procedure. It can either derive an optimal solution $(u^?, \mathbf{x}^?)$ with a finite optimal value or identify some $u^? ∈ \mathcal{U}$ for which the second-stage decision problem is infeasible. in the latter case is set to $+\infty$ by convention（按照惯例）.

我们将上面描述的内容，总结为C&CG算法的伪代码。在给出伪代码之前，我们先把C&CG算法中全局的上界和下界的计算再搬过来复习一下：

C&CG下的全局上下界的计算与Benders-dual method中的相同

C&CG算法的伪代码如下：

图5-C&CG算法的伪代码1

图6-C&CG算法的伪代码2

这里就需要注意了，在算法迭代的起始，我们可以只求解(11)中带有第一个约束的部分，并且设置 $\eta=0$ ，得到，将带入到中进行求解，得到一个解。
但是这里需要注意，文中是假设是可以首先通过对偶，再用一个oracle来顺利求解的。这个地方其实我是有一些疑惑的。但是不影响理解算法。也就是说，通过这个oracle，求解第二阶段的模型SP，然后得到cut，加回到MP中去。这里的oracle需要说明一下，这个概念来源于理论计算机，中文名称为谕示机（Oracle Machine），是带有“黑箱”的图灵机，由图灵本人提出，“黑箱”就是一个谕示，经过一个谕示就可以得到一个问题的判断结果。
这部分原文是这么描述的

Note that constraints (12)–(13) generated in Step 4(a) serve as optimality cuts and constraints (14) generated in Step 4(b) serve as feasibility cuts. In fact, because constraint (12) with $\mathbf{x}^{k+1}$ for an infeasible scenario is also valid, we can simply generate both (12) and (13) for any identified scenario. Therefore, it yields a unified approach to deal with optimality and feasibility. Next, if the second-stage problem is LP and the relatively complete recourse assumption holds, this algorithm terminates in a finite number of iterations (see A2 in the Electronic Companion [11] for the proof).(原文)

图7-最优切

在步骤4(a)中添加的约束，其中包括了约束(7)( $\eta \geqslant \mathbf{b^T x}^l$ ),是最优的割平面，可以看做是benders中的optimality cut。原文中给出了Benders-Dual 与 CCG 的一些不同点，归纳如下：

①主问题中的决策变量。CCG通过在每次迭代中引入一组新变量 $x^l$ 来增加解空间的维度，而Bders-Dual的决策变量不变；
②适用性：CCG提供了一种处理可行性问题的通用方法，对第二阶段的变量类型没有要求；而目前的Benders-Dual是特定于问题的，要求第二阶段的决策问题为线性规划；
③CCG复杂度更低：Benders-Dual的复杂度为 $O(pq)$ ，而CCG的复杂度为 $O(p)$ ；
④CCG添加的cut相较于Benders-Dual相对来说更好的实现recourse(补救) 7-

图7-Benders-Dual 与 CCG 的不同点

本文的参考论文假设不确定变量属于一个不确定的场景集合，即是一个离散的集合，我们可以比较容易地通过枚举比较每种情景下的解的质量；接下来就是考虑很复杂的情况，即多面体不确定集polyhedral uncertainty set。处理多面体不确定集的求解方法主要有outer approximation algorithm和mixed integer linear reformulations，前者基于启发式，后者利用不确定集的特殊结构将双线性规划转化为等价的MIP。然而，准确地求解具有一般多面体不确定性集的两阶段反问题仍然是一个具有挑战性的问题。
这里我们把前面提到的求解方法再次回顾一遍：

$\\ \begin{aligned} \mathbf{SP_2}: \mathcal{Q}(\mathbf{y}) = \{{ \max_{ u \in \mathcal{U}} \min_{x}\ }\mathbf{b^T x}:\mathbf{G} \mathbf{x} \geqslant \mathbf{h - Ey-M}u,\mathbf{x}\in \mathbf{S_x}\} \end{aligned}$

求解的可选方法：
1. 启发式算法
2. 通过将的inner level的部分取对偶，将变成一阶段的bilinear的二次规划模型，然后使用Gurobi，COPT等求解器来求解；
3. 使用Karush-Kuhn-Tucker(KKT)条件将inner level的模型做等价转化，然后直接变为single leve的MIP进行求解。
本文着重介绍基于KKT条件的解法。这里一定注意：KKT条件应用的前提是强对偶成立。

1.问题是凸的，且强对偶成立 $\Leftrightarrow$ KKT条件成立。(KKT条件成立的充要条件)
2.问题是非凸的，且强对偶成立 $\Rightarrow$ KKT条件成立。(KKT条件成立的必要条件)

以为例子，我们将inner level的模型（inner level是凸的，且强对偶显然成立，因此KKT条件成立是充要条件）：

$\\ \begin{aligned} \min_{ \mathbf{x}} \quad & \mathbf{b}^T \mathbf{x} \\ s.t.\quad& \mathbf{Gx \geqslant h - E\bar{y} - M}u , \qquad \\ & \mathbf{S_x} \subseteq \mathbb{R}_{+}^{m} \end{aligned}$
通过拉格朗日对偶结合KKT条件，等价转化成下面的KKT条件方程组： $\\ \begin{aligned} &\mathbf{G} \mathbf{x} \geqslant \mathbf{h - Ey-M}u&& \\ &\mathbf{G^T} \mathbf{\pi} \leqslant \mathbf{b} && \\ & (\mathbf{G} \mathbf{x}-\mathbf{h + Ey + M}u)_i\mathbf{\pi}_i=0, \quad \forall i && \\ & \mathbf{(b-G^T\pi)}_jx_j=0, \quad \forall j && \\ & \mathbf{u} \in \mathcal{U}, \quad \mathbf{x} \in \mathbf{S_x}, \quad \mathcal{\pi}\ge0. && \end{aligned}$
求解这个方程组，就等价于求解了inner level的模型。

这里需要强调一句：

KKT条件，实际上就是把一个有目标的约束规划问题，转化成了一个无目标的方程组。求解无目标的方程组，得到的解就是有目标的约束规划问题的最优解。将解带入目标函数，得到目标函数值，就是有目标的约束规划问题的最优值。
相信很多小伙伴一定很疑惑，上面的KKT条件方程组是怎么来的呢？别着急，我们在后文中有完整的推导过程。

然后我们将outer level的部分考虑进来，则整个bi-level的就等价为下面single level的模型：

的等价形式:
$\\ \begin{aligned} \max_{\mathbf{x}, u, \pi} \quad &\mathbf{b^T x} && &&(7) \\ s.t. \quad &\mathbf{G} \mathbf{x} \geqslant \mathbf{h - Ey-M}u&& &&(8) \\ &\mathbf{G^T} \mathbf{\pi} \leqslant \mathbf{b} && &&(9) \\ & (\mathbf{G} \mathbf{x}-\mathbf{h + Ey + M}u)_i\mathbf{\pi}_i=0, \quad \forall i && &&(10) \\ & \mathbf{(b-G^T\pi)}_jx_j=0, \quad \forall j && &&(11) \\ & \mathbf{u} \in \mathcal{U}, \quad \mathbf{x} \in \mathbf{S_x}, \quad \mathcal{\pi}\ge0. && &&(12) \end{aligned}$

在构造的新模型中含有两个等式约束，即(10)和(11)，有没有熟悉的感觉，没错，(10)和(11)就是互补松弛性条件。
上述模型中，约束(10)和(11)是非线性的，但是由于是乘积等于0的形式，因此可以进行等价线性化，从而将其进一步等价转化成MIP，从而直接可以求解。
具体转化方法，见下面截自论文中的图片。

图8-线性化转化方法-1

图9-线性化转化方法-2

补充：为什么用KKT

由于后面需要用到KKT条件求解，但是一些小伙伴可能对KKT条件并不熟悉，因此我们先在这里补充讲解一下KKT条件。
KKT(Karush-Kuhn-Tucker)条件，是非线性规划领域里最重要的理论成果之一，是确定某点为极值点的必要条件。对于凸优化，KKT就是优化极值点的充分必要条件。但实际上这话是不严谨的，还缺少一个regularity条件，在使用KKT条件之前需要验证regularity条件，否则就无法保证KKT条件给出的结论一定成立。（证明请看王源师兄的推文——【学界】关于KKT条件的深入探讨）
首先，我们考虑下面的约束优化问题
$\\ \begin{aligned} \min \quad z=&f(x_1,x_2,···,x_n) && \\ s.t. \quad &g_1(x_1,x_2,···,x_n)\le b_1&& \\ & g_2(x_1,x_2,···,x_n)\le b_2 \\ & \quad \quad \quad \quad \quad ·& \\ & \quad \quad \quad \quad \quad ·& \\ & \quad \quad \quad \quad \quad ·& \\ & g_m(x_1,x_2,···,x_n)\le b_m&& \\ \end{aligned}$
其中 $f,g_1,···,g_m$ 可微且一阶导数连续。

如果有是一个最优解，那么肯定是满足这个约束条件的，其次还一定存在一组非负实数( $\bar{\lambda}_1,\bar{\lambda}_2,···,\bar{\lambda}_m$ )使得下式成立
$\\ \begin{aligned} & \frac{\partial f(\bar{x})}{\partial{x_j}}-\sum^m_{i=1}\bar{\lambda}_i\frac{\partial g_i(\bar{x})}{\partial{x_j}}=0, &&j=1,2,···,n \\ &\bar{\lambda}_i[b_i-g_i(\bar{x})]=0, &&i=1,2,···,m; \\ & \bar{\lambda}_i\ge 0 &&i=1,2,···,m \end{aligned}$
其中 $\lambda_i$ 是引入的广义拉格朗日乘子，就是对应约束的影子价格，上式也就是这个优化问题的KKT条件。

下面我们用一个具体的例子来加深对KKT条件的理解。
除此之外，前文介绍到，KKT条件成立的一个必要条件是，强对偶成立。那么什么情况下强对偶成立呢？这就不得不提到Slater&#39;s condition。
Slater&#39;s条件

Slater&#39;s condition:如果满足原问题是凸优化问题，并且至少存在一个相对内点（什么叫相对内点？
答：就是一个可以让所有不等式约束都取严格 $>$ 或者的可行点），那么该问题就具有强对偶性。

例如，考虑以下凸优化问题 $\\ \begin{aligned} \min \quad &f_0(x)&& \\ s.t. \quad &f_i(x)\le 0,& \forall i=1,2,···,m& \\ & Ax=b \\ \end{aligned}$
该问题中如果至少存在一个可行解 $x'$ 满足：使不等式约束 $f_i(x')$ 取严格，即 $f_i(x')<0, \forall i=1,2,···,m$ ，那么该问题的强对偶性成立。
综上，KKT条件成立的前提：

① 目标函数和约束函数均可微且一阶导数连续；
② 问题强对偶成立(满足Slater&#39;s condition)；

KKT条件求解非线性规划的小例子

考虑下面的二次规划问题：

$\\ \begin{aligned} \max \quad z=&x_1^2-2x_2^2+27x_1+45x_2-10x_3 && \\ s.t. \quad &x_1+x_2-x_3\le 0&& \\ & x_3-17.25\le 0 \\ \end{aligned}$
首先，我们来验证该问题是否为convex或者concave。我们得到目标函数 $\\ \max z=x_1^2-2x_2^2+27x_1+45x_2-10x_3$
的海塞矩阵为
$\\ H=\left[ \begin{matrix} -2& 0& 0\\ 0& -4& 0\\ 0& 0& 0\\ \end{matrix} \right]$
所以该问题的 $一阶顺序主子式=-2<0$ ， $二阶顺序主子式=8>0$ ， $三阶主子式=0$ ，因此海塞矩阵是半负定的，所以该三元函数是concave的，并且concave函数的极小值就是全局最小值。

我们首先写出上述约束规划的拉格朗日函数 $\begin{aligned} \mathcal{L}(\mathbf{x}, \lambda_1, \lambda_2) = x_1^2-2x_2^2+27x_1+45x_2-10x_3 + \lambda_1(x_1+x_2-x_3) +\lambda_2 (x_3-17.25) \end{aligned}$
其对应的KKT条件为 $\\ \begin{aligned} &\frac{\partial \mathcal{L}}{\partial x_1} = 0 \rightarrow \quad -2x_1+27-\lambda_1=0&& (1) \\ &\frac{\partial \mathcal{L}}{\partial x_2} = 0 \rightarrow \quad -4x_2+45-\lambda_1=0&& (2) \\ &\frac{\partial \mathcal{L}}{\partial x_3} = 0 \rightarrow \quad -10+\lambda_1-\lambda_2=0&& (3) \\ &\frac{\partial \mathcal{L}}{\partial \lambda_1} = 0 \rightarrow \quad -x_1+x_2-x_3=0&& (4) \\ &\frac{\partial \mathcal{L}}{\partial \lambda_2} = 0 \rightarrow \quad x_3-17.25=0&& (5) \\ &\text{互补松弛条件1}\rightarrow \quad \lambda_1(-x_1-x_2+x_3)=0&& (6) \\ &\text{互补松弛条件2}\rightarrow \quad \lambda_2(x_3-17.25)=0&& (7) \\ &\lambda_1,\lambda_2\ge 0&& (8) \\ \end{aligned}$
下面我们来求解KKT条件方程组，以得到最优解。我们分下面几种情况讨论：

① $\lambda_1=0,\lambda_2=0$ ，很容易发现条件(3)是不满足的，所以无解；
②   $\lambda_1=0, \lambda_2>0$ ,条件(3)也是不能够满足的，所以无解；
③   $\lambda_1>0, \lambda_2=0$ ,根据条件(3)我们可以得到，再根据条件(1)-(2)可以得到 $x_1=\frac{17}{2},x_2=\frac{35}{4}$ ，然后由条件(4)得到 $x_3=\frac{69}{4}=17.25$ ，那么我们就找到了一组最优解 $(x_1=\frac{17}{2},x_2=\frac{35}{4},x_3=\frac{69}{4})$ ，目标函数值为 $z=225.375$
④   $\lambda_1>0, \lambda_2>0$ ，根据条件(5)可以得到 $x_3=17.25$ ，再结合上约束(1)-(2)与(4)，可以得到 $x_1= \frac{17}{2},x_2=\frac{35}{4}$ ，进而得到了  ，但是与条件(3)矛盾，因此无解。
所以最后得到的最优解是 $(x_1=\frac{17}{2},x_2=\frac{35}{4},x_3=\frac{69}{4},\lambda_1=10,\lambda_2=0)$ ，最优目标值为 $\begin{aligned} z=&x_1^2-2x_2^2+27x_1+45x_2-10x_3 \\ =&\frac{17}{2}\times \frac{17}{2}-2\times \frac{35}{4}\times \frac{35}{4}+27\times \frac{17}{2}+45\times \frac{35}{4}-10\times \frac{69}{4} \\ =&369.875 \end{aligned}$

与之前描述的相同，求解KKT条件方程组，就等价于求解了原来的约束最优化问题。
推导的inner level模型的KKT条件

我们来看second-stage的部分，也就是

$\\ \begin{aligned} \max_{u \in \mathcal{U}} \min_{ \mathbf{x}\in {\mathbf{F}(\mathbf{y}, u )}} \quad & \mathbf{b}^T \mathbf{x} \\ s.t. \quad & \mathbf{Gx \geqslant h - Ey - M}u , \qquad u \in \mathcal{U} \hspace{1cm} \rightarrow \pi && \\ & \mathbf{x} \geqslant \mathbf{0} \hspace{1cm} \rightarrow \theta && \end{aligned}$
其中和 $\theta$ 分别是两个约束的对偶变量。这里需要注意， $\mathbf{x} \geqslant \mathbf{0}$ 也要视作约束。

我们将约束变成 $\leqslant$ ，则为

$\\ \mathbf{h - Ey - M}u -\mathbf{Gx} \leqslant 0, \qquad u \in \mathcal{U} \\ \mathbf{-x} \leqslant \mathbf{0}$

只考虑inner level的部分。我们将约束松弛松弛到目标函数中去，构造拉格朗日函数：

$\\ \begin{aligned} \mathcal{L}(\mathbf{x}, \pi, \theta) = \mathbf{b}^T \mathbf{x} + (\mathbf{h - Ey - M}u -\mathbf{Gx}) \pi - \mathbf{x}^T \theta \end{aligned}$

根据KKT条件，得到下面的式子：

$\begin{aligned} \frac{\partial \mathcal{L}}{\partial x} = 0 \rightarrow &\qquad \mathbf{b}^T - \mathbf{G} \pi - \theta= 0 && (1) \\ \frac{\partial \mathcal{L}}{\partial \pi} = 0 \rightarrow &\qquad \mathbf{h - Ey - M}u -\mathbf{Gx} = 0 && (2) \\ \frac{\partial \mathcal{L}}{\partial \theta} = 0 \rightarrow & \qquad \mathbf{x} = 0 && (3) \\ \text{complementary slackness conditions} \rightarrow &\qquad (\mathbf{h - Ey - M}u -\mathbf{Gx}) \pi_i= 0 && (4) \\ \text{complementary slackness conditions} \rightarrow &\qquad \mathbf{x}^T \theta= \mathbf{0} && (5) \\ & \pi_i \geqslant 0 && (6) \\ & \theta_i \geqslant 0 && (7) \end{aligned}$
(根据上面的文本，(2),(3)其实是不需要的，但是根据一个博客以及一些参考文献，确实需要(2)和(3))，博客地址: https://www.cnblogs.com/90zeng/p/Lagrange_duality.html#top

根据(1)和(5)，可得 $\begin{aligned} (\mathbf{b}^T - \mathbf{G} \pi) \mathbf{x^T} = \mathbf{0} \end{aligned}$
根据(1)和(7)，可得

$\\ \begin{aligned} (\mathbf{b}^T - \mathbf{G} \pi) = \theta \geqslant 0 \end{aligned}$
也就是
$\\ \begin{aligned} \mathbf{b}^T - \mathbf{G} \pi \geqslant 0 \qquad \rightarrow \qquad \mathbf{G^T} \pi \leqslant \mathbf{b} \end{aligned}$

此外，原问题约束为

$\\ \begin{aligned} & \mathbf{Gx \geqslant h - Ey - M}u , \qquad u \in \mathcal{U} \hspace{1cm} \rightarrow \pi \\ & \mathbf{x} \geqslant \mathbf{0} \hspace{1cm} \rightarrow \theta \end{aligned}$

因此，使用KKT条件处理inner level之后的的等价形式为

$\\ \begin{aligned} \max_{u \in \mathcal{U}} \quad & \mathbf{b}^T \mathbf{x} && && (1) \\ s.t.\quad & \mathbf{Gx \geqslant h - Ey - M}u , \qquad && u \in \mathcal{U} && (2) \\ & \mathbf{G^T} \pi \leqslant \mathbf{b} && && (3) \\ & (\mathbf{Gx - h - Ey - M}u)_i \pi_i = 0 , \qquad && u \in \mathcal{U}, \forall i && (4) \\ & (\mathbf{b - G^T} \pi)_j x_j = 0 && \forall j && (5) \\ & \mathbf{x} \geqslant \mathbf{0} && && (6) \\ & \pi \geqslant 0 \end{aligned}$

约束(4)和(5)是非线性约束。但是由于乘积为0，因此可以进行等价线性化。引入0-1辅助变量

$v_i$ : 如果 $v_i = 0$ ，则 $\pi_i=0$ ；如果 $v_i = 1$ ，则 $\pi_i >0$ ，即相当于 $\\(\mathbf{Gx - h - Ey - M}u)_i=0$
$w_i$ : 如果 $w_i = 0$ ，则 $x_i=0$ ；如果 $w_i = 1$ ，则 $x_i >0$ ，即相当于 $\\(\mathbf{b - G^T} \pi)_i=0$

上述非线性模型可以线性化为下面的等价形式：

该模型为一个MIP，可以直接很方便的进行求解。
经过了如上一系列的操作，我们来总结一下完整的使用C&CG算法求解Two-stage RO的步骤。
C&CG算法求解Two-stage RO的完整简洁步骤

Two-stage RO原问题为

$\\ \begin{aligned} \min_{\mathbf{y}} \quad & \mathbf{c^Ty} +\max_{u \in \mathcal{U}} \min_{ \mathbf{x}\in {\mathbf{F}(\mathbf{y}, u )}} \mathbf{b}^T \mathbf{x} \\ s.t.\quad & \mathbf{Ay \geqslant d}, \\ & \mathbf{Gx \geqslant h - Ey - M}u , \qquad u \in \mathcal{U} \\ & \mathbf{S_y} \subseteq \mathbb{R}_{+}^{n} \\ & \mathbf{S_x} \subseteq \mathbb{R}_{+}^{m} \end{aligned}$

我们使用C&CG求解该Two-stage RO模型，需要将其分解为如下的主问题Master Problem：

$\\ \begin{aligned} \min_{\mathbf{y}} \quad &\mathbf{c^T y} + \eta && \\ s.t. \quad & \mathbf{Ay \geqslant d}&& \\ & \eta \geqslant \mathbf{b^T x}^l, \quad && l=1, \cdots, r \\ & \mathbf{Ey + G} \mathbf{x}^l \geqslant \mathbf{h - M}u_l, \quad && l = 1, \cdots, r \\ & \mathbf{y \in S_y} \\ & \mathbf{x}^l \in \mathbf{S_x}, && l=1, \cdots, r. \end{aligned}$

以及如下的子问题Subproblem：

$\\ \begin{aligned} \max_{u \in \mathcal{U}} \min_{ \mathbf{x}\in {\mathbf{F}(\mathbf{y}, u )}} \quad & \mathbf{b}^T \mathbf{x} \\ s.t. \quad & \mathbf{Gx \geqslant h - Ey - M}u , \qquad u \in \mathcal{U} \\ & \mathbf{S_x} \subseteq \mathbb{R}_{+}^{m} \end{aligned}$

通过对inner level使用KKT条件进行reformulation，则Subproblem可以变化为

$\\ \begin{aligned} \max_{u \in \mathcal{U}} \quad & \mathbf{b}^T \mathbf{x} && \\ s.t. \quad & \mathbf{Gx \geqslant h - Ey - M}u , \qquad && u \in \mathcal{U} \\ & \mathbf{G^T} \pi \leqslant \mathbf{b} && \\ & (\mathbf{Gx - h - Ey - M}u)_i \pi_i = 0 , \qquad && u \in \mathcal{U}, \forall i \\ & (\mathbf{b - G^T} \pi)_j x_j = 0 && \forall j \\ & \mathbf{x} \geqslant \mathbf{0} \\ & \pi \geqslant 0 \end{aligned}$

将subproblem等价线性化，则变成

C&CG算法实际上就是不断迭代求解主问题Master Problem，以及重构后的子问题subproblem,从而实现快速求解Two-stage RO模型的。
Case study: robust location-transportation problem

Two-stage robust location-transportation problem

这部分原文是这么描述的

Consider the following location-transportation problem. To supply a commodity to customers, it will be first stored at potential facilities and then be transported to $n$ customers. The fixed cost of the building facilities at site is and the unit capacity cost is for $i = 1, \cdots,m$ . The demand is $d_j$ for $j = 1, . . . , n$ , and the unit transportation cost between and $j$ is for $i − j$ pair. The maximal allowable capacity of the facility at site is $K_i$ and $\sum_{i}{K_i} \geqslant \sum_{j}{d_j}$ ensures feasibility. Let be the facility location variable, $z_i \in \mathbb{R}+$ be the capacity variable, and be the transportation variable. Then, the nominal formulation of this location-transportation problem is as follows:

总结一下：设施选址和存储上限在第一阶段确定，运输安排将在第二阶段确定，以满足客户需求。参数

:: 运输成本
：fixed cost
: unit capacity cost

决策变量

: facility location variable
$z_i \in \mathbb{R}_+$ ：the capacity variable
: transportation variable

插一句题外话：这里为什么要引入 $z_i$ 呢，其实就是为了让第二阶段的问题变成一个线性规划。因为第二阶段的问题，是不涉及到 $y$ 的，只涉及到，而是连续变量。

模型如下：

$\\ \begin{aligned} \min_{\mathbf{y,z,x}} \quad & \sum_{i}{f_iy_i} + \sum_{i}{a_iz_i} + \sum_{i}\sum_{j}{c_{ij}x_{ij}} && && (22) \\ s.t. \quad& z_i \leqslant K_i y_i && \forall i && (23) \\ & \sum_{j}{x_{ij}} \leqslant z_i && \forall i && (24) \\ & \sum_{i}{x_{ij}} \geqslant d_j && \forall j && (25) \\ & y_{i} \in \{0,1 \}, z_i \geqslant 0, && \forall i && (26) \\ & x_{ij} \geqslant 0, && \forall i, \forall j && (26-2) \end{aligned}$
这里一定要注意，需要保证第二阶段的模型一定是可行的，所以是有一些额外约束的 $\sum_{i} \left(\sum_j {x_{ij}} \right) \leqslant \sum_i {z_i} \\ \sum_j {d_j} \leqslant \sum_{j} \left(\sum_i {x_{ij}}\right)$
因此，我们有
$\\ \sum_j {d_j} \leqslant \sum_i {z_i}$
i.e., $\\ \sum_i {z_i} \geqslant 206 + 274 + 220 + 40 \times 1.8 = 772$

另外在实际问题中，因为需求是未知的，还需要处理这种不确定性。这部分原文是这么描述的

The objective function in (22) is to minimize the overall cost, including the fixed cost, capacity cost, and transportation cost. Constraints in (23) and (24) require that capacity can be installed only at a site with a built facility, and the supply cannot exceed the capacity. Constraints in (25) guarantee that the demand is satisfied. In practice, the demand is unknown before any facility is built and capacity is installed. A popular way to capture that uncertainty is as follows [7,15,1]:

$\\ \begin{aligned} \mathbf{D} = \{ \mathbf{d}: d_j = \underline{d}_j + g_j \tilde{d}_j , g_j \in [0,1], \sum_j {g_j \leqslant \Gamma}, j=1, \cdots, n\} \end{aligned}$
其实就是

图10-需求不确定性-1

这里，注意 $\mathbf{S_y} = { \mathbf{(y,z)}\in {0,1}^m \times \mathbb{R}_{+}^m: (23) }$ 里面的符号。

图11-需求不确定性-2

Experimental results and discussion

这部分原文是这么描述的

Next, we employed both C&CG and Benders-dual methods to study this two-stage robust problem. The detailed formulations of master and sub-problems are omitted here but provided in A4 in the Electronic Companion [11]. In all of our experiments, CPLEX 12.4 was used as the solver to the master problem and the oracle to the linearized subproblem. For both the master problem and subproblems, the optimality tolerance was set to 104. Both the C&CG and Benders-dual algorithms were implemented in C++ on a desktop Dell OPTIPLEX 760 (Intel Core 2 Duo CPU, 3.0 GHz, 3.25 GB of RAM) in a Windows 7 environment. We first study dynamic behaviors of the C&CG and Bendersdual methods on a small scale. An illustrative problem is given with three potential facilities, three customers, and a general polyhedral uncertainty set. The deterministic formulation is presented as follows:

下面是确定性的模型：

$\\ \begin{aligned} \min \quad &400y_0 + 414y_1 + 326y_2 \\ &+18z_0 + 25z_1 + 20z_2\\ &+22x_{00} + 33x_{01} + 24x_{02} + 33x_{10} + 23x_{11} \\ &+ 30x_{12}+20x_{20} + 25x_{21} + 27x_{22} \\ s.t. \quad & z_i ≤ 800y_i, && \forall i= 0, 1, 2; \\ & \sum_{j}{x_{ij}} \leqslant z_i, && \forall i= 0, 1, 2; \\ & \sum_{i}{x_{ij}} \geqslant d_j, && \forall j= 0, 1, 2; \\ & y_{i} \in \{0,1 %C}, z_i \geqslant 0, && \forall i= 0, 1, 2; \\ & x_{ij} \geqslant 0, && \forall i= 0, 1, 2;, \forall j= 0, 1, 2; \end{aligned}$

但是如果需求是不确定的，则需要使用鲁棒优化。文献中将其建模为一个两阶段鲁棒优化模型。
其中，不确定集合的定义如下：

图11-不确定集合的定义

其实就是

$\\ \begin{aligned} \mathbf{D} = & \{ \mathbf{d}: \} \\ s.t.\quad & d_0 = 206 + 40g_0, \\ & d_1 = 274 + 40g_1, \\ & d_2 = 220 + 40g_2, \\ & 0 ≤ g_0 ≤ 1, 0 ≤ g_1 ≤ 1, 0 ≤ g_2 ≤ 1, \\ & g_0 + g_1 + g_2 ≤ 1.8, \\ &g_0 + g_1 ≤ 1.2 \end{aligned}$

两阶段鲁棒优化的模型如下

接下来我们使用前面介绍的C&CG算法来求解上述两阶段鲁棒优化模型。
主问题为：

$\\ \begin{aligned} \mathbf{MP_2}: \quad &\min_{\mathbf{y,z \in S_y}, \eta}\quad \sum_{i}{f_iy_i} + \sum_{i}{a_iz_i} + \eta \\ s.t. \quad & z_i ≤ 800y_i, && \forall i= 0, 1, 2; \\ & \sum_i {z_i} \geqslant 772 && \\ & \eta \geqslant \sum_{i}\sum_{j}{c_{ij}x_{ij}^{l}}, \quad && l=1, \cdots, r \\ & \sum_{j}{x_{ij}^l} \leqslant z_i, && \forall i= 0, 1, 2; l=1, \cdots, r \\ & \sum_{i}{x_{ij}^l} \geqslant d_j^{l*}, && \forall j= 0, 1, 2; l=1, \cdots, r \\ & \eta \in \mathbb{R}. \\ & y_{i} \in \{0,1 \}, z_i \geqslant 0, && \forall i= 0, 1, 2; \\ & x_{ij} \geqslant 0, && \forall i= 0, 1, 2;, \forall j= 0, 1, 2; \end{aligned}$

子问题为 (注意，其中可以看做输入参数，由上一阶段的模型给出)

$\\ \begin{aligned} \mathbf{SP_2}: \mathcal{Q}(\mathbf{y}) =& \{ \max_{\mathbf{d \in D}} \min_{\mathbf{x}}\sum_{i}\sum_{j}{c_{ij}x_{ij}}\} \\ s.t.\quad & \sum_{j}{x_{ij}} \leqslant z_i^*, && \forall i= 0, 1, 2; \qquad \rightarrow \pi \\ & \sum_{i}{-x_{ij}} \leqslant -d_j, && \forall j= 0, 1, 2; \qquad \rightarrow \theta \\ & x_{ij} \geqslant 0, && \forall i= 0, 1, 2;, \forall j= 0, 1, 2; \end{aligned}$

我们对其进行对偶，得到

$\\ \begin{aligned} \mathbf{SP_2}: \mathcal{Q}(\mathbf{y}) =& \max_{\mathbf{d \in D}, \pi, \theta} \,\,\sum_{i}{z_i^* \pi_i} - \sum_{j}{d_j \theta_j} \\ s.t.\quad & \pi_i - \theta_j \leqslant c_{ij}, && \forall i= 0, 1, 2; \forall j= 0, 1, 2; \qquad \rightarrow x_{ij} \\ & \pi_i \leqslant 0, && \forall i= 0, 1, 2; \\ & \theta_j \leqslant 0, && \forall j= 0, 1, 2; \end{aligned}$
注意，这里，还需要加上

因此就是 $\\ \begin{aligned} \mathbf{SP_2}: \mathcal{Q}(\mathbf{y}) =& \max_{\mathbf{d \in D}, \pi, \theta} \,\,\sum_{i}{z_i^* \pi_i} - \sum_{j}{d_j \theta_j} \\ s.t.\quad & \pi_i - \theta_j \leqslant c_{ij}, && \forall i= 0, 1, 2; \forall j= 0, 1, 2; \qquad \rightarrow x_{ij} \\ & d_j = \underline{d}_j + g_j \tilde{d}_j && \forall j= 0, 1, 2; \\ & \sum_j {g_j \leqslant \Gamma}, \qquad && j=1, \cdots, n \\ & g_j \in [0,1], \qquad && j=1, \cdots, n \\ & \pi_i \leqslant 0, && \forall i= 0, 1, 2; \\ & \theta_j \leqslant 0, && \forall j= 0, 1, 2; \\ & d_j \geqslant 0, && \forall j= 0, 1, 2; \end{aligned}$

所以在初始情况下，我们发现，该目标函数是一个非线性的，也就是bi-linear的。直接使用Gurobi求解该模型也是可以的，但是效率应该会比求解MIP版本的低一些。
我们这里还是使用前文介绍的KKT条件的办法来求解。
我们写出KKT条件(注意，互补松弛的部分一定要是写成 $\geqslant 0$ 的形式的左端项)

$\begin{aligned} \mathbf{SP_2}: \mathcal{Q}(\mathbf{y}) =& \max \sum_{i}\sum_{j}{c_{ij}x_{ij}} \\ s.t.\quad & \sum_{j}{x_{ij}} \leqslant z_i^*, && \forall i= 0, 1, 2; \qquad \rightarrow \pi \\ & \sum_{i}{-x_{ij}} \leqslant -d_j, && \forall j= 0, 1, 2; \qquad \rightarrow \theta \\ &\pi_i - \theta_j \leqslant c_{ij}, && \forall i= 0, 1, 2; \forall j= 0, 1, 2; \qquad \rightarrow x_{ij} \\ & \left(z_i^* - \sum_{j}{x_{ij}} \right) \pi_i = 0, && \forall i= 0, 1, 2; \qquad \qquad (a) \\ & \left(\sum_{i}{x_{ij}} -d_j \right) \theta_j = 0, && \forall j= 0, 1, 2; \qquad \qquad (b) \\ & (c_{ij} - \pi_i+ \theta_j)x_{ij} =0, && \forall i= 0, 1, 2; \forall j= 0, 1, 2; \qquad \qquad (c) \\ & d_j = \underline{d}_j + g_j \tilde{d}_j && \forall j= 0, 1, 2; \\ & \sum_j {g_j \leqslant \Gamma}, \qquad && j=1, \cdots, n \\ & g_j \in [0,1], \qquad && j=1, \cdots, n \\ & \pi_i \leqslant 0, && \forall i= 0, 1, 2; \\ & \theta_j \leqslant 0, && \forall j= 0, 1, 2; \\ & d_j \geqslant 0, && \forall j= 0, 1, 2; \end{aligned}$

其中，约束(a)和(b)和( c )是非线性的，我们可以将其进行等价线性化处理。下面的三个约束

$\\ \begin{aligned} & \left(z_i^* - \sum_{j}{x_{ij}} \right) \pi_i = 0, && \forall i= 0, 1, 2; \qquad \qquad (a) \\ & \left(\sum_{i}{x_{ij}} -d_j \right) \theta_j = 0, && \forall j= 0, 1, 2; \qquad \qquad (b) \\ & (c_{ij} - \pi_i+ \theta_j)x_{ij} =0, && \forall i= 0, 1, 2; \forall j= 0, 1, 2; \qquad \rightarrow x_{ij} \end{aligned}$

可以等价线性化为:

$\\ \begin{aligned} &-\pi_i \leqslant Mv_i, \\ & z_i^* - \sum_{j}{x_{ij}} \leqslant M(1-v_i) \\ &-\theta_j \leqslant Mw_j, \\ & \sum_{i}{x_{ij}} -d_j \leqslant M(1-w_j) \\ &x_{ij} \leqslant Mh_{ij}, \\ & c_{ij} - \pi_i+ \theta_j \leqslant M(1-h_{ij}) \\ & v_i, w_j, h_{i,j} \in \{0, 1\}. \end{aligned}$
其中引入的3个辅助变量的含义与之前介绍的类似。

综上，最终等价形式为

$\begin{aligned} \mathbf{MP_2}: \quad &\min_{\mathbf{y,z \in S_y}, \eta}\quad \sum_{i}{f_iy_i} + \sum_{i}{a_iz_i} + \eta \\ s.t. \quad & z_i ≤ 800y_i, && \forall i= 0, 1, 2; \\ & \sum_i {z_i} \geqslant 772 && \\ & \eta \geqslant \sum_{i}\sum_{j}{c_{ij}x_{ij}^{l}}, \quad && l=1, \cdots, r \\ & \sum_{j}{x_{ij}^l} \leqslant z_i^{l*}, && \forall i= 0, 1, 2; l=1, \cdots, r \\ & \sum_{i}{x_{ij}^l} \geqslant d_j^{l*}, && \forall j= 0, 1, 2; l=1, \cdots, r \\ & \eta \in \mathbb{R}. \\ & y_{i} \in \{0,1 \}, z_i \geqslant 0, && \forall i= 0, 1, 2; \\ & x_{ij} \geqslant 0, && \forall i= 0, 1, 2;, \forall j= 0, 1, 2; \end{aligned}$

子问题经过KKT条件转化之后的等价形式为

$\begin{aligned} \mathbf{SP_2}: \mathcal{Q}(\mathbf{y}) =& \max \sum_{i}\sum_{j}{c_{ij}x_{ij}} \\ s.t.\quad & \sum_{j}{x_{ij}} \leqslant z_i^*, && \forall i= 0, 1, 2; \qquad \rightarrow \pi \\ & \sum_{i}{-x_{ij}} \leqslant -d_j, && \forall j= 0, 1, 2; \qquad \rightarrow \theta \\ &\pi_i - \theta_j \leqslant c_{ij}, && \forall i= 0, 1, 2; \forall j= 0, 1, 2; \qquad \rightarrow x_{ij} \\ &-\pi_i \leqslant Mv_i, && \forall i= 0, 1, 2; \\ & z_i^* - \sum_{j}{x_{ij}} \leqslant M(1-v_i) && \forall i= 0, 1, 2; \\ &-\theta_j \leqslant Mw_j, && \forall j= 0, 1, 2; \\ & \sum_{i}{x_{ij}} -d_j \leqslant M(1-w_j), && \forall j= 0, 1, 2; \\ &x_{ij} \leqslant Mh_{ij}, && \forall i= 0, 1, 2; \forall j= 0, 1, 2; \\ & c_{ij} - \pi_i+ \theta_j \leqslant M(1-h_{ij}) && \forall i= 0, 1, 2; \forall j= 0, 1, 2; \\ & v_i, w_j, h_{i,j} \in \{0, 1\}. && \forall i= 0, 1, 2; \forall j= 0, 1, 2; \\ & d_j = \underline{d}_j + g_j \tilde{d}_j && \forall j= 0, 1, 2; \\ & \sum_j {g_j \leqslant \Gamma}, \qquad && j=1, \cdots, n \\ & g_j \in [0,1], \qquad && j=1, \cdots, n \\ & \pi_i \leqslant 0, && \forall i= 0, 1, 2; \\ & \theta_j \leqslant 0, && \forall j= 0, 1, 2; \\ & d_j \geqslant 0, && \forall j= 0, 1, 2; \end{aligned}$

论文中给出了这个算例的数值实验结果，如下：

The upper and lower bounds of the two algorithms are presented inTable 1, which clearly shows the superiority of the C&CG method.（CCG与Benders Dual的对比）

图12-CCG与Benders Dual的对比

接下来我们使用Python调用Gurobi来实现C&CG算法，并且求解上述算例。
Python调用Gurobi求解确定性的模型

确定性的location-transportation problem的模型

$\\ \begin{aligned} \min \quad &400y_0 + 414y_1 + 326y_2 \\ &+18z_0 + 25z_1 + 20z_2\\ &+22x_{00} + 33x_{01} + 24x_{02} + 33x_{10} + 23x_{11} \\ &+ 30x_{12}+20x_{20} + 25x_{21} + 27x_{22} \\ s.t. \quad & z_i ≤ 800y_i, && \forall i= 0, 1, 2; \\ & \sum_{j}{x_{ij}} \leqslant z_i, && \forall i= 0, 1, 2; \\ & \sum_{i}{x_{ij}} \geqslant d_j, && \forall j= 0, 1, 2; \\ & y_{i} \in \{0,1 \}, z_i \geqslant 0, && \forall i= 0, 1, 2; \\ & x_{ij} \geqslant 0, && \forall i= 0, 1, 2;, \forall j= 0, 1, 2; \end{aligned}$

Python调用Gurobi完整代码

首先定义好参数
from gurobipy import *
import numpy as np

&#34;&#34;&#34; The input parameter &#34;&#34;&#34;
facility_num = 3
customer_num = 3
fixed_cost = [400, 414, 326]
unit_capacity_cost = [18, 25, 20]
trans_cost = [[22, 33, 24],
            [33, 23, 30],
            [20, 25, 27]]
max_capacity = 800

demand_nominal = [206, 274, 220]
demand_var = [40, 40, 40]然后建模求解
&#34;&#34;&#34; The deterministic model  &#34;&#34;&#34;

&#34;&#34;&#34; Create variables  &#34;&#34;&#34;
model = Model(&#39;Location-transportation model&#39;)
x = {}
z = {}
y = {}

for i in range(facility_num):
y = model.addVar(vtype = GRB.BINARY, name = &#39;y_%s&#39;%(i))
z = model.addVar(lb = 0, vtype = GRB.CONTINUOUS, name = &#39;y_%s&#39;%(i))
for j in range(customer_num):
      x[i, j] = model.addVar(lb = 0, vtype = GRB.CONTINUOUS, name = &#39;x_%s_%s&#39;%(i, j))

&#34;&#34;&#34; Set objective &#34;&#34;&#34;
obj = LinExpr()
for i in range(facility_num):
obj.addTerms(fixed_cost, y)
obj.addTerms(unit_capacity_cost, z)
for j in range(customer_num):
      obj.addTerms(trans_cost, x[i, j])
model.setObjective(obj, GRB.MINIMIZE)

&#34;&#34;&#34; Add Constraints  &#34;&#34;&#34;
# cons 1
for i in range(facility_num):
model.addConstr(z <= max_capacity * y)

# cons 2
for i in range(facility_num):
expr = LinExpr()
for j in range(customer_num):
      expr.addTerms(1, x[i, j])
model.addConstr(expr <= z)

# cons 3
for j in range(facility_num):
expr = LinExpr()
for i in range(customer_num):
      expr.addTerms(1, x[i, j])
model.addConstr(expr >= demand_nominal[j])

&#34;&#34;&#34; solve the model and output &#34;&#34;&#34;
model.optimize()
print(&#39;Obj = {}&#39;.format(model.ObjVal))
print(&#39;-----  location ----- &#39;)
for key in z.keys():
print(&#39;facility : {}, location: {}, capacity: {}&#39;.format(key, y[key].x, z[key].x))运行结果如下
Root relaxation: objective 3.312739e+04, 5 iterations, 0.00 seconds

Nodes | Current Node |    Objective Bounds    |    Work
Expl Unexpl |  Obj  Depth IntInf | Incumbent BestBd Gap | It/Node Time

   0    0 33127.3886 0 2       - 33127.3886    -    - 0s
H 0    0                   33504.000000 33127.3886  1.12%    - 0s

Explored 1 nodes (5 simplex iterations) in 0.01 seconds
Thread count was 16 (of 16 available processors)

Solution count 1: 33504

Optimal solution found (tolerance 1.00e-04)
Best objective 3.350400000000e+04, best bound 3.350400000000e+04, gap 0.0000%
Obj = 33504.0
-----  location -----
facility : 0, location: 1.0, capacity: 490.0
facility : 1, location: -0.0, capacity: 0.0
facility : 2, location: 1.0, capacity: 282.0Python调用Gurobi求解Two-stage RO model： C&CG algorithm
C&CG算法的伪代码

图13-CCG伪代码-1

图14-CCG伪代码-2

C&CG算法的主问题和子问题

C&CG算法求解Two-stage RO的主问题的最终形式为

$\\ \begin{aligned} \mathbf{MP_2}: \quad &\min_{\mathbf{y,z \in S_y}, \eta}\quad \sum_{i}{f_iy_i} + \sum_{i}{a_iz_i} + \eta \\ s.t. \quad & z_i ≤ 800y_i, && \forall i= 0, 1, 2; \\ & \sum_i {z_i } \geqslant 772 \\ & \eta \geqslant \sum_{i}\sum_{j}{c_{ij}x_{ij}^{l}}, \quad && l=1, \cdots, r \\ & \sum_{j}{x_{ij}^l} \leqslant z_i^{l*}, && \forall i= 0, 1, 2; l=1, \cdots, r \\ & \sum_{i}{x_{ij}^l} \geqslant d_j^{l*}, && \forall j= 0, 1, 2; l=1, \cdots, r \\ & \eta \in \mathbb{R}. \\ & y_{i} \in \{0,1 \}, z_i \geqslant 0, && \forall i= 0, 1, 2; \\ & x_{ij} \geqslant 0, && \forall i= 0, 1, 2;, \forall j= 0, 1, 2; \end{aligned}$
C&CG算法求解Two-stage RO的子问题的最终形式为

$\\ \begin{aligned} \mathbf{SP_2}: \mathcal{Q}(\mathbf{y}) =& \max \sum_{i}\sum_{j}{c_{ij}x_{ij}} \\ s.t.\quad & \sum_{j}{x_{ij}} \leqslant z_i^*, && \forall i= 0, 1, 2; \qquad \rightarrow \pi \\ & \sum_{i}{-x_{ij}} \leqslant -d_j, && \forall j= 0, 1, 2; \qquad \rightarrow \theta \\ &\pi_i - \theta_j \leqslant c_{ij}, && \forall i= 0, 1, 2; \forall j= 0, 1, 2; \qquad \rightarrow x_{ij} \\ &-\pi_i \leqslant Mv_i, && \forall i= 0, 1, 2; \\ & z_i^* - \sum_{j}{x_{ij}} \leqslant M(1-v_i) && \forall i= 0, 1, 2; \\ &-\theta_j \leqslant Mw_j, && \forall j= 0, 1, 2; \\ & \sum_{i}{x_{ij}} -d_j \leqslant M(1-w_j), && \forall j= 0, 1, 2; \\ &x_{ij} \leqslant Mh_{ij}, && \forall i= 0, 1, 2; \forall j= 0, 1, 2; \\ & c_{ij} - \pi_i+ \theta_j \leqslant M(1-h_{ij}) && \forall i= 0, 1, 2; \forall j= 0, 1, 2; \\ & v_i, w_j, h_{i,j} \in \{0, 1\}. && \forall i= 0, 1, 2; \forall j= 0, 1, 2; \\ & d_j = \underline{d}_j + g_j \tilde{d}_j && \forall j= 0, 1, 2; \\ & \sum_j {g_j \leqslant \Gamma}, \qquad && j=1, \cdots, n \\ & g_j \in [0,1], \qquad && j=1, \cdots, n \\ & \pi_i \leqslant 0, && \forall i= 0, 1, 2; \\ & \theta_j \leqslant 0, && \forall j= 0, 1, 2; \\ & d_j \geqslant 0, && \forall j= 0, 1, 2; \end{aligned}$
Python调用Gurobi实现C&CG算法

由于篇幅原因，本推文仅展示部分代码，完整代码见文末获取方式。
from gurobipy import *
import numpy as np

&#34;&#34;&#34; The input parameter &#34;&#34;&#34;
facility_num = 3
customer_num = 3
fixed_cost = [400, 414, 326]
unit_capacity_cost = [18, 25, 20]
trans_cost = [[22, 33, 24],
            [33, 23, 30],
            [20, 25, 27]]
max_capacity = 800

demand_nominal = [206, 274, 220]
demand_var = [40, 40, 40]

&#34;&#34;&#34; build initial master problem &#34;&#34;&#34;
&#34;&#34;&#34; Create variables &#34;&#34;&#34;
master = Model(&#39;master problem&#39;)
x_master = {}
z = {}
y = {}
eta = master.addVar(lb=0, vtype=GRB.CONTINUOUS, name=&#39;eta&#39;)

for i in range(facility_num):
y = master.addVar(vtype=GRB.BINARY, name=&#39;y_%s&#39; % (i))
z = master.addVar(lb=0, vtype=GRB.CONTINUOUS, name=&#39;z_%s&#39; % (i))

&#34;&#34;&#34; Set objective &#34;&#34;&#34;
obj = LinExpr()
for i in range(facility_num):
obj.addTerms(fixed_cost, y)
obj.addTerms(unit_capacity_cost, z)
obj.addTerms(1, eta)

master.setObjective(obj, GRB.MINIMIZE)

&#34;&#34;&#34; Add Constraints  &#34;&#34;&#34;
# cons 1
for i in range(facility_num):
master.addConstr(z <= max_capacity * y)

&#34;&#34;&#34; Add initial value Constraints  &#34;&#34;&#34;
# create new variables x
iter_cnt = 0
for i in range(facility_num):
for j in range(customer_num):
      x_master[iter_cnt, i, j] = master.addVar(lb=0
                                             , ub=GRB.INFINITY
                                             , vtype=GRB.CONTINUOUS
                                             , name=&#39;x_%s_%s_%s&#39; % (iter_cnt, i, j))
# create new constraints
expr = LinExpr()
for i in range(facility_num):
for j in range(customer_num):
      expr.addTerms(trans_cost[j], x_master[iter_cnt, i, j])
master.addConstr(eta >= expr)

expr = LinExpr()
for i in range(facility_num):
expr.addTerms(1, z)
master.addConstr(expr >= 772)  # 206 + 274 + 220 + 40 * 1.8

&#34;&#34;&#34; solve the model and output &#34;&#34;&#34;
master.optimize()
print(&#39;Obj = {}&#39;.format(master.ObjVal))
print(&#39;-----  location ----- &#39;)
for key in z.keys():
print(&#39;facility : {}, location: {}, capacity: {}&#39;.format(key, y[key].x, z[key].x))

def print_sub_sol(model, d, g, x):
d_sol = {}
if(model.status != 2):
      print(&#39;The problem is infeasible or unbounded!&#39;)
      print(&#39;Status: {}&#39;.format(model.status))
      d_sol[0] = 0
      d_sol[1] = 0
      d_sol[2] = 0
else:
      print(&#39;Obj(sub) : {}&#39;.format(model.ObjVal), end=&#39;\t | &#39;)
      for key in d.keys():
         # print(&#39;demand: {}, perturbation = {}&#39;.format(d[key].x, g[key].x))
         d_sol[key] = d[key].x
return d_sol

&#34;&#34;&#34; Column-and-constraint generation &#34;&#34;&#34;

LB = -np.inf
UB = np.inf
iter_cnt = 0
max_iter = 30
cut_pool = {}
eps = 0.001
Gap = np.inf

z_sol = {}
for key in z.keys():
z_sol[key] = z[key].x
# print(z_sol)

&#34;&#34;&#34; solve the master problem and update bound &#34;&#34;&#34;
master.optimize()

&#34;&#34;&#34;
Update the Lower bound
&#34;&#34;&#34;
LB = master.ObjVal
print(&#39;LB: {}&#39;.format(LB))

&#39;&#39;&#39; create the subproblem &#39;&#39;&#39;
subProblem = Model(&#39;sub problem&#39;)
x = {}  # transportation decision variables in subproblem
d = {}  # true demand
g = {}  # uncertainty part: var part
pi = {}  # dual variable
theta = {}  # dual variable
v = {}  # aux var
w = {}  # aux var
h = {}  # aux var
big_M = 100000
for i in range(facility_num):
pi = subProblem.addVar(lb=-GRB.INFINITY, ub=0, vtype=GRB.CONTINUOUS, name=&#39;pi_%s&#39; % i)
v = subProblem.addVar(vtype=GRB.BINARY, name=&#39;v_%s&#39; % i)
for j in range(customer_num):
w[j] = subProblem.addVar(vtype=GRB.BINARY, name=&#39;w_%s&#39; % j)
g[j] = subProblem.addVar(lb=0, ub=1, vtype=GRB.CONTINUOUS, name=&#39;g_%s&#39; % j)
theta[j] = subProblem.addVar(lb=-GRB.INFINITY, ub=0, vtype=GRB.CONTINUOUS, name=&#39;theta_%s&#39; % j)
d[j] = subProblem.addVar(lb=0, ub=GRB.INFINITY, vtype=GRB.CONTINUOUS, name=&#39;d_%s&#39; % j)
for i in range(facility_num):
for j in range(customer_num):
      h[i, j] = subProblem.addVar(vtype=GRB.BINARY, name=&#39;h_%s_%s&#39; % (i, j))
      x[i, j] = subProblem.addVar(lb=0, ub=GRB.INFINITY, vtype=GRB.CONTINUOUS, name=&#39;x_%s_%s&#39; % (i, j))

&#34;&#34;&#34; set objective &#34;&#34;&#34;
sub_obj = LinExpr()
for i in range(facility_num):
for j in range(customer_num):
      sub_obj.addTerms(trans_cost[j], x[i, j])
subProblem.setObjective(sub_obj, GRB.MAXIMIZE)

&#34;&#34;&#34; add constraints to subproblem &#34;&#34;&#34;
# cons 1
for i in range(facility_num):
expr = LinExpr()
for j in range(customer_num):
      expr.addTerms(1, x[i, j])
subProblem.addConstr(expr <= z_sol, name=&#39;sub_capacity_1_z_%s&#39; % i)

# cons 2
for j in range(facility_num):
expr = LinExpr()
for i in range(customer_num):
      expr.addTerms(1, x[i, j])
subProblem.addConstr(expr >= d[j])

# cons 3
for i in range(facility_num):
for j in range(customer_num):
      subProblem.addConstr(pi - theta[j] <= trans_cost[j])

&#34;&#34;&#34; demand constraints &#34;&#34;&#34;
for j in range(customer_num):
subProblem.addConstr(d[j] == demand_nominal[j] + g[j] * demand_var[j])

subProblem.addConstr(g[0] + g[1] + g[2] <= 1.8)
subProblem.addConstr(g[0] + g[1] <= 1.2)

&#34;&#34;&#34; logic constraints &#34;&#34;&#34;
# logic 1
....

# logic 2
....

# logic 3
....

subProblem.write(&#39;SP.lp&#39;)
subProblem.optimize()
d_sol = {}

print(&#39;\n\n\n *******          C&CG starts       *******  &#39;)
print(&#39;\n **             Initial Solution          ** &#39;)

d_sol = print_sub_sol(subProblem, d, g, x)

&#34;&#34;&#34;
Update the initial Upper bound
&#34;&#34;&#34;
# UB = min(UB, subProblem.ObjVal + master.ObjVal - eta.x)
# print(&#39;UB (iter {}): {}&#39;.format(iter_cnt, UB))

# close the outputflag
master.setParam(&#39;Outputflag&#39;, 0)
subProblem.setParam(&#39;Outputflag&#39;, 0)
&#34;&#34;&#34;
Main loop of CCG algorithm
&#34;&#34;&#34;
while (UB - LB > eps and iter_cnt <= max_iter):
iter_cnt += 1
# print(&#39;\n\n --- iter : {} --- \n&#39;.format(iter_cnt))
print(&#39;\n iter : {} &#39;.format(iter_cnt), end=&#39;\t | &#39;)

# create new variables x
for i in range(facility_num):
      for j in range(customer_num):
         x_master[iter_cnt, i, j] = master.addVar(lb=0
                                                   , ub=GRB.INFINITY
                                                   , vtype=GRB.CONTINUOUS
                                                   , name=&#39;x_%s_%s_%s&#39; % (iter_cnt, i, j))

# if subproblem is frasible and bound, create variables xk+1 and add the new constraints
if (subProblem.status == 2 and subProblem.ObjVal < 1000000000):
      # create new constraints
      expr = LinExpr()
      for i in range(facility_num):
         for j in range(customer_num):
            expr.addTerms(trans_cost[j], x_master[iter_cnt, i, j])
      master.addConstr(eta >= expr)

      # create worst case related constraints
      # cons 2
      ....

      # cons 3
      ....

      # solve the resulted master problem
      master.optimize()
      print(&#39;Obj(master): {}&#39;.format(master.ObjVal), end=&#39;\t | &#39;)
      &#34;&#34;&#34; Update the LB &#34;&#34;&#34;
      LB = master.ObjVal
      print(&#39;LB (iter {}): {}&#39;.format(iter_cnt, LB), end=&#39;\t | &#39;)

      &#34;&#34;&#34; Update the subproblem &#34;&#34;&#34;
      # first, get z_sol from updated master problem
      for key in z.keys():
         z_sol[key] = z[key].x

      # change the coefficient of subproblem
      ....

      # add new constraints
      # cons 1
      ....

      # logic 1
      ....

      &#34;&#34;&#34; Update the lower bound &#34;&#34;&#34;
      subProblem.optimize()
      d_sol = print_sub_sol(subProblem, d, g, x)

      &#34;&#34;&#34;
      Update the Upper bound
      &#34;&#34;&#34;
      if (subProblem.status == 2):
         UB = min(UB, subProblem.ObjVal + master.ObjVal - eta.x)
         # print(&#39;eta = {}&#39;.format(eta.x))
      print(&#39;UB (iter {}): {}&#39;.format(iter_cnt, UB), end=&#39;\t | &#39;)
      Gap = round(100 * (UB - LB) / UB, 2)
      print(&#39;eta = {}&#39;.format(eta.x), end=&#39;\t | &#39;)
      print(&#39; Gap: {} %  &#39;.format(Gap), end=&#39;\t&#39;)

      # If the subproblem is unbounded
if (subProblem.status == 4):
      # create worst case related constraints
      # cons 2
      ....

      # cons 3
      ....

      # solve the resulted master problem
      master.optimize()
      print(&#39;Obj(master): {}&#39;.format(master.ObjVal))

      &#34;&#34;&#34; Update the LB &#34;&#34;&#34;
      LB = master.ObjVal
      print(&#39;LB (iter {}): {}&#39;.format(iter_cnt, LB))

      &#34;&#34;&#34; Update the subproblem &#34;&#34;&#34;
      # first, get z_sol from updated master problem
      for key in z.keys():
         z_sol[key] = z[key].x

      # change the coefficient of subproblem
      ....

      # add new constraints
      # cons 1
      ....

      # logic 1
      ....

      &#34;&#34;&#34; Update the lower bound &#34;&#34;&#34;
      subProblem.optimize()
      d_sol = print_sub_sol(subProblem, d, g, x)

      &#34;&#34;&#34;
      Update the Upper bound
      &#34;&#34;&#34;
      if (subProblem.status == 2):
         UB = min(UB, subProblem.ObjVal + master.ObjVal - eta.x)
         print(&#39;eta = {}&#39;.format(eta.x))
      print(&#39;UB (iter {}): {}&#39;.format(iter_cnt, UB))
      Gap = round(100 * (UB - LB) / UB, 2)
      print(&#39;---- Gap: {} % ---- &#39;.format(Gap))

master.write(&#39;finalMP.lp&#39;)
print(&#39;\n\nOptimal solution found !&#39;)
print(&#39;Opt_Obj : {}&#39;.format(master.ObjVal))
print(&#39; **  Final Gap: {} %  **  &#39;.format(Gap))
print(&#39;\n  ** Solution  **  &#39;)
for i in range(facility_num):
print(&#39; {}: {},\t  {}: {} &#39;.format(y.varName, y.x, z.varName, z.x), end=&#39;&#39;)
print(&#39;\t actual demand: {}: {}&#39;.format(d.varName, d.x), end=&#39;&#39;)
print(&#39;\t perturbation in worst case: {}: {}&#39;.format(g.varName, g.x))
print(&#39;\n  ** Transportation solution  **  &#39;)
for i in range(facility_num):
for j in range(customer_num):
      if (x[i, j].x > 0):
         print(&#39;trans: {}: {}, cost : {} \t &#39;.format(x[i, j].varName, x[i, j].x, trans_cost[j] * x[i, j].x),
               end=&#39;&#39;)
print()
*******          C&CG starts       *******

**             Initial Solution          **
Obj(sub) : 20942.0 |
iter : 1 | Obj(master): 33680.0  | LB (iter 1): 33680.0  | Obj(sub) : 18034.0 | UB (iter 1): 33696.0  | eta = 18018.0    |  Gap: 0.05 %
iter : 2 | Obj(master): 33680.0  | LB (iter 2): 33680.0  | Obj(sub) : 18430.0 | UB (iter 2): 33680.0  | eta = 18430.0    |  Gap: 0.0 %

Optimal solution found !
Opt_Obj : 33680.0
**  Final Gap: 0.0 %  **

  ** Solution  **
y_0: 1.0, z_0: 458.0    actual demand: d_0: 206.0 perturbation in worst case: g_0: 0.0
y_1: 0.0, z_1: 0.0 actual demand: d_1: 314.0 perturbation in worst case: g_1: 1.0
y_2: 1.0, z_2: 314.0    actual demand: d_2: 252.0 perturbation in worst case: g_2: 0.8000000000000007

  ** Transportation solution  **
trans: x_0_0: 206.0, cost : 4532.0 trans: x_0_2: 252.0, cost : 6048.0

trans: x_2_1: 314.0, cost : 7850.0

Process finished with exit code 0参考文献

1.Zeng B, Zhao L. Solving two-stage robust optimization problems using a column-and-constraint generation method[J]. Operations Research Letters, 2013, 41(5): 457-461.
2.https://www.cnblogs.com/90zeng/p/Lagrange_duality.html#top
3.【学界】关于KKT条件的深入探讨，https://zhuanlan.zhihu.com/p/33229011

完整集成版本代码获取方式

需要完整代码的小伙伴，可以通过下面的方式获取完整代码：
代码获取方式：

转发该推文至朋友圈，集齐30赞，并将集赞截图发送至公众号后台，即可免费获得。

转发该推文至朋友圈，集齐15赞，并打赏10元，之后将集赞截图和打赏截图发送至公众号后台，即可获得。

不想集赞的小伙伴，可以直接打赏20元，并将打赏截图发送至公众号后台，即可获得。

注意，小编将在2天之后，统一将代码发送给各位小伙伴。2天之内不发送，望各位读者周知。
立。
<hr/>欢迎关注我们的微信公众号 运小筹
<hr/>往期推文合集
第94篇：优化 | 手把手教你用C++调用Gurobi求解CVRP：环境配置、数学模型及完整代码
第93篇：OR Talk Pool：近期研讨会、暑期学校和学术报告推荐
第92篇：优化求解器 | Gurobi的MVar类：矩阵建模利器、求解对偶问题的备选方案 (附详细案例+代码)
第91篇：【求解器】Gurobi中非线性约束的对偶值的获取：QCP、QCQP、SOCP | 以及求解器最优性理论证明
第90篇：优化 | 二部图最大匹配问题的精确算法详解(HK算法和匈牙利算法)：一份让您满意的【理论介绍+代码实现】学习笔记
第89篇：优化算法 | Benders Decomposition: 一份让你满意的【入门-编程实战-深入理解】的学习笔记
第88篇：优化 | 史上最【皮】数学题求解的新尝试：一种求解器的视角 (Python调用Gurobi实现)
第87篇：优化 | 寻找新的建模视角——从直观解释对偶问题入手：以Cutting Stock Problem的对偶问题为例
第86篇：ORers‘ Bling Chat【03】 | 【高光聊天记录集锦】：运小筹读者群里那些热烈的讨论
第85篇：非线性优化 | 非线性问题线性化以及求解的详细案例及Python+Gurobi求解
第84篇：ORers&#39; Bling Chat | 【高光聊天记录集锦-02】：运小筹读者群里那些热烈的讨论
第83篇：Machine-Learning–Based Column Selection for Column Generation
第82篇：最新！205页运小筹优化理论学习笔记发布(2021年9月--2022年3月)！
第81篇：【鲁棒优化】| 补充证明:为什么最优解一定满足取值等于绝对值(论文笔记：The Price of Robustness)
第80篇：【鲁棒优化】| 论文笔记：The Price of Robustness - 列不确定性模型部分的推导和一些思考
第79篇：ORers&#39; Bling Chat | 【高光聊天记录集锦-01】：运小筹读者群里那些热烈的讨论
第78篇：优化| 手把手教你学会杉数求解器(COPT)的安装、配置与测试
第77篇：【教学视频】优化 | 线性化(2)：连续变量 * 0-1变量的线性化
第76篇：【教学视频】优化 | 线性化：两个0-1变量相乘的线性化
第75篇：强化学习实战 | DQN和Double DQN保姆级教程：以Cart-Pole为例
第74篇：强化学习| 概念梳理：强化学习、马尔科夫决策过程与动态规划
第73篇：强化学习实战 | Q-learning求解最短路(SPP)问题
第72篇：鲁棒优化 | 以Coding入门鲁棒优化：以一个例子引入（二）
第71篇：鲁棒优化|基于ROME编程入门鲁棒优化：以一个例子引入(一)
第70篇：优化|含分式的非线性规划求解: 基于Outer Approximation的Branch-and-cut 算法及其Java实现
第69篇：科研工具 | 手把手教你玩转文献管理神器：Endnote
第68篇：相约深圳 | 2022 INFORMS 服务科学国际会议·征稿通知
第67篇：鲁棒优化｜利用rome求解鲁棒优化简单案例：以CVaR不确定性集为例
第66篇：机器学习 | 交通流特征工程小技巧和思考
第65篇：最新！145页运小筹优化理论学习笔记发布(2021年4月--9月)！
第64篇：优化 | 手把手教你用Python调用SCIP求解最优化模型
第63篇：优化 | 随机规划案例：The value of the stochastic solution
第62篇：工具 | draw.io: 科研流程示意图必备大杀器
第61篇：优化 | 开源求解器SCIP的安装与入门教程（Windows+Python）
第60篇：优化|Gurobi处理非线性目标函数和约束的详细案例
第59篇：让出租车更“懂”您的出行
第58篇：测试算例下载：《运筹优化常用模型、算法及案例实战：代码手册》
第57篇：鲁棒优化|分布式鲁棒优化转化为SOCP案例及代码实现(Python+RSOME)
第56篇：鲁棒优化 | 分布式鲁棒优化简单案例及实战(RSOME+Gurobi)
第55篇：【尾款+发货】|【代码手册】《运筹优化常用模型、算法及案例实战：Python+Java
第54篇：深度强化学习之：PPO训练红白机1942
第53篇：简单装配线平衡问题
第52篇：【视频讲解】CPLEX的OPL语言：可视化的优化模型求解解决方案
第51篇：算法 | 基于英雄联盟寻路背景的A星算法及python实现
第50篇：【转发】清华大学深圳国际研究生院2021年物流工程与管理项目优秀大学生夏令营报名通知
第49篇：优化 | 精确算法之分支定界介绍和实现（附Python代码）
第48篇：【顶刊论文速递】综述：服务科学视角下的金融科技
第47篇：【重新发布】|《运筹优化常用模型、算法及案例实战：Python+Java实现》【代码手册】开始预购啦！！！
第46篇：智慧交通可视化：手把手教你用Python做出行数据可视化-osmnx 包入门教程
第45篇：优化 | Pick and delivery problem的介绍与建模实现（二）
第44篇：推一个知乎学弱猹的公众号
第43篇：元启发式算法 | 遗传算法(GA)解决TSP问题(Python实现)
第42篇：优化|视频详解Python调用Gurobi的应用案例：TSP
第41篇：最新！213页运小筹优化理论系列笔记发布！
第40篇：运小筹第四期月刊发布！
第39篇：开源交通仿真平台SimMobility的安装教程
第38篇：浅谈 | P问题、NP问题、NPC问题、NP-Hard问题
第37篇：一份掏心掏肺的信息素养笔记分享
第36篇：强化学习｜Q-learning (王者荣耀视角)
第35篇：优化｜高级建模方法（Gurobi）：线性化表达小技巧
第34篇：深度强化学习介绍 | Deep Q Network——理论与代码实现
第33篇：优化 | 列生成算法及Java调用cplex实现
第32篇：优化 | Pick and delivery problem的简介与建模实现（一）
第31篇：最新！运小筹月刊-1月份发布！
第30篇：“Learn to Improve”（L2I）：ICLR文章分享 | 运用RL解决VRP问题
第29篇：线性规划求解小程序 | 基于Python 的Cplex 求解器图形化界面
第28篇：运筹学与管理科学TOP期刊揭秘 —TR系列
第27篇：Julia安装与配置Jupyter Notebook
第26篇：期刊征文| IEEE TRANSACTIONS应对COVID-19的特刊征文消息速递
第25篇：两阶段随机规划(Two-stage Stochastic Programming)：一个详细的例子
第24篇：最新！运小筹月刊-12月份发布！
第23篇：Python调用Gurobi：Assignment Problem简单案例
第22篇：初识随机规划：一个小小例子
第21篇：机器学习运用到VRP的若干小知识
第20篇：运筹学与管理科学TOP期刊揭秘 —Service Science
第19篇：手把手教你用Python实现Dijkstra算法（伪代码+Python实现）
第18篇：运筹学与管理科学大揭秘—TOP期刊主编及研究方向一览
第17篇：优化 | 手把手教你用Python实现动态规划Labeling算法求解SPPRC问题
第16篇：代码 | 运小筹GitHub项目上线啦，快来标星收藏吧！！
第15篇：最新！运小筹月刊首次发布！
第14篇：优化| 手把手教你用Python实现Dijkstra算法（附Python代码和可视化）
第13篇：优化|手把手教你用Python调用Gurobi求解最短路问题
第12篇：优化| 手把手教你用Java实现单纯形法
第11篇：优化| 手把手教你用Python调用Gurobi求解VRPTW
第10篇：优化 | 两阶段启发式算法求解带时间窗车辆路径问题（附Java代码）
第9篇：Java调用cplex求解运输问题
第8篇：优化 | 如何优雅地写出大规模线性规划的对偶
第7篇：优化 | TSP中两种不同消除子环路的方法及Callback实现（Python调用Gurobi求解）
第6篇：元启发式算法 | 禁忌搜索(Tabu Search)解决TSP问题(Python代码实现)
第5篇：论文代码复现 | 无人机与卡车联合配送(Python+Gurobi)
第4篇：优化|Shortest Path Problem及其对偶问题的一些探讨(附Python调用Gurobi实现)
第3篇：可视化|Python+OpenStreetMap的交通数据可视化(二)：OpenStreetMap+Python画城市交通路网图
第2篇：优化|最速下降法：详细案例+Python实现
第1篇：Python+networkX+OpenStreetMap实现交通数据可视化(一)：用OpenStreetMap下载地图数据

<hr/>作者：刘兴禄，清华大学，清华-伯克利深圳学院（博士在读）
审校：张瑞三，四川大学，硕士在读，邮箱：zrssum3@stu.scu.edu.cn 知乎ID：MunSum3

DomDomm · 发表于 2022-7-20 09:54

读者有问题可以直接私聊我

		自动登录	找回密码
密码			立即注册

鲁棒优化| C&CG算法求解两阶段鲁棒优化：全网最完整、最 ...

本帖子中包含更多资源