基于线性规划的随机规划两阶段问题的L形算法求解示例

一、 问题描述

想象你是一家电力公司的决策者，面临一个经典的两阶段决策问题。

1. 第一阶段（此时决策）： 在不确定的未来风电出力（即能发多少电）具体数值已知之前，你需要提前决定从火电厂购买多少电量。假设每单位电量的价格为 \(c\) 元。设购买量为 \(x\)， 成本为 \(c x\)。这是一个**“此时”决策**。

2. 第二阶段（事后决策）： 第二天，风电的实际出力 \(omega\) 被观测到。如果实际风电出力高，你可能买多了火电，需要支付惩罚 \(q_1\) 来处理多余电量（如储存或弃用）。如果风电出力低，你买少了，需要启动昂贵的紧急备用电源（如燃气轮机）来弥补缺口，单位成本为 \(q_2\)。这是一个**“事后”决策**，取决于第一阶段决策 \(x\) 和随机实现 \(omega\)。

你的目标是：找到一个第一阶段购买量 \(x\)，使得第一阶段成本与第二阶段期望成本之和最小。

形式化模型：

• 设风电出力 \(omega\) 是一个随机变量，服从已知的概率分布。
• 设第二阶段，实际风电出力为 \(omega\)， 则电力缺口/盈余为 \(omega - x\)。
  • 如果 \(omega - x > 0\)， 表示风电多，有盈余，惩罚为 \(q_1 (omega - x)\)。
  • 如果 \(omega - x < 0\)， 表示风电少，有缺口，惩罚为 \(q_2 (x - omega)\)。
• 目标：\(min_{x ge 0} left{ c x + mathbb{E}_{omega}[Q(x, omega)] right}\)
  其中，\(Q(x, omega) = min_{y^+, y^- ge 0} left{ q_1 y^+ + q_2 y^- : y^+ - y^- = omega - x right}\)。

问题的挑战： 随机变量 \(omega\) 的取值空间可能非常大（甚至是连续的），直接计算期望 \(mathbb{E}_{omega}[Q(x, omega)]\) 非常复杂，因为 \(Q(x, omega)\) 本身是一个优化问题（第二阶段子问题）的最优值。L形算法就是为了高效求解此类两阶段随机线性规划而设计的经典方法。

二、 核心思想：Benders分解 / L形算法

L形算法的本质是一种特殊的Benders分解算法，应用于两阶段随机规划。

1. 分解：

   • 将原问题分解为一个主问题（只包含第一阶段变量 \(x\) 和一个代表“期望未来成本”的代理变量 \(theta\)）。
   • 和多个子问题（对应每个可能的随机场景 \(omega\)， 给定一个试探解 \(hat{x}\) 后，求解第二阶段的成本 \(Q(hat{x}, omega)\)）。

2. 迭代思想：

   • 首先，忽略未来成本的复杂性，在主问题中假设 \(theta\) 可以取任意值（通常先设其下界为 \(-infty\)），得到一个初始的第一阶段解 \(hat{x}\)。
   • 然后，用这个 \(hat{x}\) 去“测试”各个子问题。子问题的解会返回两种信息：
     • 最优性割： 如果子问题有最优解，它会提供一个关于 \(Q(x, omega)\) 的线性下界估计，这个下界是 \(x\) 的函数。将其添加到主问题，可以“教导”主问题更准确地估计未来成本。
     • 可行性割： 如果给定的 \(hat{x}\) 导致某个子问题不可行（例如，某些约束无法满足），子问题会返回一个线性不等式，禁止主问题再产生这类不可行的 \(x\)。
   • 主问题在添加了这些割平面后重新求解，得到新的 \(x\) 和 \(theta\)。不断迭代，直到主问题中的 \(theta\) 与通过子问题计算出的真实期望成本足够接近。

三、 逐步求解过程

我们用一个简化例子说明。假设风电出力 \(omega\) 只有两种可能场景，概率各0.5：

• 场景1：高出力，\(omega_1 = 8\) 单位
• 场景2：低出力，\(omega_2 = 2\) 单位
  参数：\(c=1\), \(q_1=1\)（盈余惩罚）, \(q_2=3\)（缺口惩罚）。
  设 \(x\) 为第一阶段购买量，\(y_1^+, y_1^-\) 为场景1的盈余和缺口变量，\(y_2^+, y_2^-\) 为场景2的变量。

两阶段问题的确定等价形式为：

\[ min quad 1 cdot x + 0.5(1 cdot y_1^+ + 3 cdot y_1^-) + 0.5(1 cdot y_2^+ + 3 cdot y_2^-) \]

\[ text{subject to:} quad x + y_1^+ - y_1^- = 8 \]

\[ qquad qquad qquad quad x + y_2^+ - y_2^- = 2 \]

\[ qquad qquad qquad quad x, y_1^+, y_1^-, y_2^+, y_2^- ge 0 \]

现在我们用L形算法来求解。

步骤1：构造主问题（MP）

主问题只包含第一阶段变量 \(x\) 和一个表示“期望第二阶段成本”的连续变量 \(theta\)。

\[ text{(MP) } min quad c x + theta \]

\[ text{s.t.} quad x ge 0, quad theta in mathbb{R} quad text{(初始无其他约束)} \]

初始时，因为没有关于 \(theta\) 的任何信息，我们可以简单地设其下界为 \(-infty\)（或一个很大的负数）。在第一次迭代中，主问题的目标会驱使 \(theta\) 趋向负无穷，这在优化中通常通过添加一个“无限制”的代理变量或直接先忽略其下界来处理。更常见的初始化是设 \(theta = 0\) 或求解一个松弛问题。为了直观，我们从一个简单的试探解开始。

实际上，算法可以从一个试探的 \(x\) 值开始。我们这里从求解一个极度松弛的主问题开始：如果我们暂时不考虑 \(theta\) 的任何约束，目标函数 \(min x + theta\) 在 \(x=0\) 和 \(theta to -infty\) 时趋于负无穷。这显然不合理。为了启动算法，我们需要一个初始试探解 \(x^k\)。一个常见且安全的选择是求解期望值问题（用 \(omega\) 的期望值代替随机变量），得到一个初始 \(x\)。

1. \(bar{omega} = 0.5*8 + 0.5*2 = 5\)
2. 求解 \(min { c x + Q(x, 5) }\)， 这等价于解一个确定性问题：
   \(min x + 1*y^+ + 3*y^-\)， s.t. \(x + y^+ - y^- = 5\)， 所有变量 \(ge 0\)。
   容易分析，最优解是让 \(x = 5\)， 使得 \(y^+=y^-=0\)， 总成本为 \(5\)。
   我们取 \(x^1 = 5\) 作为第一次迭代的试探解。

步骤2：求解子问题（给定 \(x^1=5\)）

对于每个场景 \(s\)， 求解子问题 \(Q(x^1, omega_s)\)：

• 子问题模型 (对给定 \(x\) 和 \(omega\)):

\[ Q(x, omega) = min_{y^+, y^- ge 0} { q_1 y^+ + q_2 y^- } \]

\[ text{s.t.} quad y^+ - y^- = omega - x \]

场景1 (\(omega_1=8\))：
约束：\(y_1^+ - y_1^- = 8 - 5 = 3\)
目标：\(min 1*y_1^+ + 3*y_1^-\)
显然，最优解是让 \(y_1^+ = 3, y_1^- = 0\)， 成本 \(Q_1 = 1*3 = 3\)。
此子问题是可行且有界的。其对偶变量 \(pi_1\) 是约束 \(y^+ - y^- = 3\) 的影子价格。由于原问题变量系数 \(c=(1,3)\)， 约束矩阵 \(A=(1, -1)\)， 根据线性规划对偶理论，最优对偶解 \(pi_1\) 需满足互补松弛条件。简单计算可得，\(pi_1\) 取值在 \(q_1\) 和 \(-q_2\) 之间。更直接地，我们可以通过观察得到：该子问题的最优基是 \(y^+\) 在基中。对偶解 \(pi_1 = q_1 = 1\)。

场景2 (\(omega_2=2\))：
约束：\(y_2^+ - y_2^- = 2 - 5 = -3\)
目标：\(min 1*y_2^+ + 3*y_2^-\)
显然，最优解是让 \(y_2^+ = 0, y_2^- = 3\)， 成本 \(Q_2 = 3*3 = 9\)。
最优基是 \(y^-\) 在基中。对偶解 \(pi_2 = -q_2 = -3\)。

步骤3：生成最优性割并更新主问题

计算当前迭代的真实期望成本：
\(mathbb{E}[Q] = 0.5*3 + 0.5*9 = 6\)

计算主问题目标值下界（基于当前主问题解 \(x^1=5, theta^1\)）：
当前主问题中 \(theta^1\) 是无约束的，在上次求解时，为了得到 \(x^1=5\)， 我们是从一个期望值模型得到的，并不是从MP求解的。现在我们正式建立MP的第一次迭代后形式。我们需要从子问题的最优解生成一个Benders最优性割。

对于两阶段线性问题，最优性割的一般形式为：

\[ theta ge sum_{s} p_s (pi_s (omega_s - x)) \]

其中 \(p_s\) 是场景概率，\(pi_s\) 是该场景子问题的最优对偶解。

代入我们的数据：

• 场景1: \(p_1=0.5, pi_1=1, omega_1=8\)
• 场景2: \(p_2=0.5, pi_2=-3, omega_2=2\)

计算割平面：
\(theta ge 0.5 * [1*(8 - x)] + 0.5 * [(-3)*(2 - x)]\)
\(theta ge 0.5*(8 - x) + 0.5*(-6 + 3x)\)
\(theta ge 0.5*(8 - x -6 + 3x)\)
\(theta ge 0.5*(2 + 2x)\)
\(theta ge 1 + x\)

因此，我们得到第一个最优性割： \(theta ge 1 + x\)

更新主问题 (MP)：

\[ min quad x + theta \]

\[ text{s.t.} quad theta ge 1 + x \]

\[ quad x ge 0 \]

步骤4：求解新的主问题

将约束 \(theta ge 1 + x\) 代入目标：\(min x + (1 + x) = min 2x + 1\)， 在 \(x ge 0\) 下，最优解在 \(x=0\) 处取得。
此时，\(theta = 1 + 0 = 1\)。
得到新的试探解：\(x^2 = 0, quad theta^2 = 1\)。
主问题目标值下界 \(LB = 0 + 1 = 1\)。

步骤5：求解新的子问题（给定 \(x^2=0\)）

场景1 (\(omega_1=8\))：
约束：\(y_1^+ - y_1^- = 8 - 0 = 8\)
最优解：\(y_1^+ = 8, y_1^- = 0\)， 成本 \(Q_1 = 8\)， 对偶解 \(pi_1 = 1\)。

场景2 (\(omega_2=2\))：
约束：\(y_2^+ - y_2^- = 2 - 0 = 2\)
最优解：\(y_2^+ = 2, y_2^- = 0\)， 成本 \(Q_2 = 2\)， 对偶解 \(pi_2 = 1\)。（注意：这里缺口为0，因为 \(omega - x = 2 > 0\)， 是盈余而非缺口，所以用便宜的 \(y^+\) 处理）。

计算真实期望成本： \(mathbb{E}[Q] = 0.5*8 + 0.5*2 = 5\)

计算主问题上界 (UB)： 这是原问题一个可行解 \((x=0)\) 对应的总成本：\(c*x + mathbb{E}[Q] = 1*0 + 5 = 5\)。
当前最优上界 \(UB = 5\)。

当前最优下界 (LB)： 来自主问题最优值，\(LB = 1\)。
\(UB > LB\)， 未收敛，继续。

步骤6：生成第二个最优性割

用新的对偶解 \(pi_1=1, pi_2=1\) 生成割：
\(theta ge 0.5*[1*(8-x)] + 0.5*[1*(2-x)]\)
\(theta ge 0.5*(8 - x + 2 - x)\)
\(theta ge 0.5*(10 - 2x)\)
\(theta ge 5 - x\)

更新主问题 (MP)：

\[ min quad x + theta \]

\[ text{s.t.} quad theta ge 1 + x quad text{(割1)} \]

\[ quad theta ge 5 - x quad text{(割2)} \]

\[ quad x ge 0 \]

步骤7：求解新的主问题

目标为 \(min x + theta\)， 约束为 \(theta ge max(1+x, 5-x)\)。
我们需要最小化 \(x + max(1+x, 5-x)\)。
令 \(f(x) = x + max(1+x, 5-x)\)。
考虑两条直线：

• 情况A: 当 \(1+x ge 5-x\) 即 \(2x ge 4\)， \(x ge 2\) 时，\(f(x) = x + (1+x) = 2x+1\)， 在 \(x=2\) 处取最小值 \(5\)。
• 情况B: 当 \(1+x le 5-x\) 即 \(x le 2\) 时，\(f(x) = x + (5-x) = 5\)， 为常数 \(5\)。

因此，全局最小值在 \(x in [0, 2]\) 时取得，值为 \(5\)。为简便，取 \(x^3 = 2\)， 则 \(theta^3 = max(1+2, 5-2) = max(3, 3) = 3\)。
新的试探解：\(x^3 = 2, quad theta^3 = 3\)。
主问题目标值下界 \(LB = 2 + 3 = 5\)。

步骤8：求解新的子问题（给定 \(x^3=2\)）

场景1 (\(omega_1=8\))：
约束：\(y_1^+ - y_1^- = 8-2=6\)
最优：\(y_1^+=6, y_1^-=0\)， 成本 \(Q_1=6\)， \(pi_1=1\)。

场景2 (\(omega_2=2\))：
约束：\(y_2^+ - y_2^- = 2-2=0\)
最优：\(y_2^+=0, y_2^-=0\)， 成本 \(Q_2=0\)， \(pi_2\) 可以是 \([-3, 1]\) 间的任意值（因为原问题变量都为0，对偶可行域是一个区间）。为了生成最强的割（最紧的下界），我们通常取能使割平面尽可能“支撑”期望成本函数的对偶解，即最优多面体的极值点。这里我们取 \(pi_2 = 1\)（另一个极值是 \(-3\)，试试哪个生成更紧的割）。

真实期望成本： \(mathbb{E}[Q] = 0.5*6 + 0.5*0 = 3\)
计算上界： 对应总成本 = \(c*x + mathbb{E}[Q] = 1*2 + 3 = 5\)。与当前下界 \(LB=5\) 相等。

收敛！ 因为 \(UB = LB = 5\)。

四、 结果与总结

最优解为：第一阶段购买电量 \(x^* = 2\) 单位。
最优总期望成本为 \(5\)。
此时：

• 如果遇到高风电场景 (\(omega=8\))， 盈余 \(6\) 单位，支付惩罚 \(6\)。
• 如果遇到低风电场景 (\(omega=2\))， 供需平衡，无额外成本。
  第一阶段成本为 \(2\)， 期望总成本为 \(2 + 0.5*6 + 0.5*0 = 5\)。

L形算法流程总结：

1. 初始化 主问题（只有第一阶段变量和 \(theta\)），设定上下界 \(UB=+infty, LB=-infty\)。
2. 求解主问题， 得到试探解 \((hat{x}, hat{theta})\)， 更新下界 \(LB = c hat{x} + hat{theta}\)。
3. 对每个随机场景， 用 \(hat{x}\) 求解第二阶段的子问题。
   • 如果所有子问题可行，计算总期望成本 \(TC = c hat{x} + sum p_s Q(hat{x}, omega_s)\)， 更新上界 \(UB = min(UB, TC)\)。从子问题的最优对偶解生成最优性割，加入主问题。
   • 如果有子问题不可行，则生成可行性割加入主问题。
4. 检查收敛：如果 \(UB - LB\) 小于预设容差，停止。否则，返回步骤2。

这个算法的强大之处在于它将一个大规模随机规划问题（场景多时变量爆炸）分解为一系列较小的问题迭代求解，主问题提供决策，子问题提供反馈（割平面），逐步逼近最优解，是求解两阶段随机线性规划的主流方法之一。