线性规划的原始-对偶内点法求解随机线性规划问题示例

字数 13695 2025-12-07 12:22:48

线性规划的原始-对偶内点法求解随机线性规划问题示例

题目描述

考虑以下随机线性规划问题：

目标：最小化 \(c^T x\)
约束：\(Ax = b, \quad x \ge 0\)
其中：
- 部分系数（例如，成本向量 \(c\)、约束矩阵 \(A\) 的某些元素、或右端向量 \(b\) 的某些分量）是随机变量，其概率分布已知。
- 我们希望在满足约束的情况下，最小化“期望总成本”，即最小化 \(\mathbb{E}[c]^T x\)，同时要求决策变量 \(x\) 是“此时必须确定、不依赖随机变量实现”的决策（在随机规划中，这被称为“第一阶段”决策或“此时此地”决策）。

为了具体化，我们设定一个简单但典型的例子。假设一个生产商需要提前决定两种产品（产品1和产品2）的生产量 \(x_1, x_2\)，以应对未来的随机需求。已知每件产品的生产成本是确定的，分别为5和4。但是，市场需求是随机的，且如果产量超过实际需求，会有单位库存持有成本（随机，与需求相关）；如果产量不足，则会有单位缺货惩罚成本（随机）。生产受资源约束。我们想最小化“期望总成本”（生产成本 + 期望库存持有与缺货成本之和）。

经过建模，该问题可以表述为一个两阶段随机线性规划，其中第一阶段决策是 \(x_1, x_2\)，第二阶段针对每种随机场景优化库存/缺货决策。在特定假设下（例如，随机参数只出现在目标函数的某些成本系数中，且其期望可计算），该问题可简化为一个确定性线性规划，其目标系数为对应随机参数的期望值。

为了演示算法，我们使用如下具体的确定性等价线性规划模型：

最小化：

\[2x_1 + 3x_2 \]

约束于：

\[\begin{aligned} x_1 + 2x_2 &\ge 4 \quad (\text{需求约束，系数为期望需求}) \\ 3x_1 + x_2 &\ge 5 \quad (\text{资源约束，期望资源消耗要求}) \\ x_1, x_2 &\ge 0 \]

注意：这里的系数 2 和 3 已经是考虑了生产、库存、缺货等随机成本后的期望单位成本。约束右端 4 和 5 是考虑了随机性后的期望要求值。

我们将使用原始-对偶内点法来求解这个线性规划。原始-对偶内点法是求解线性规划的一类高效内点法，它通过求解由原始可行性、对偶可行性和互补松弛条件组成的非线性方程组，并沿中心路径逼近最优解。

解题过程

步骤1：将问题转化为标准形式

原始-对偶内点法通常处理标准形式的线性规划：

\[\begin{aligned} \min \quad & c^T x \\ \text{s.t.} \quad & Ax = b, \\ & x \ge 0 \end{aligned} \]

其中，\(A \in \mathbb{R}^{m \times n}\)，通常假设 \(A\) 行满秩。

我们的例子是：

\[\begin{aligned} \min \quad & 2x_1 + 3x_2 \\ \text{s.t.} \quad & x_1 + 2x_2 \ge 4, \\ & 3x_1 + x_2 \ge 5, \\ & x_1, x_2 \ge 0. \end{aligned} \]

引入松弛变量 \(x_3, x_4 \ge 0\)，将不等式化为等式：

\[\begin{aligned} x_1 + 2x_2 - x_3 &= 4, \\ 3x_1 + x_2 - x_4 &= 5. \end{aligned} \]

现在，令

\[x = \begin{pmatrix} x_1 \\ x_2 \\ x_3 \\ x_4 \end{pmatrix}, \quad c = \begin{pmatrix} 2 \\ 3 \\ 0 \\ 0 \end{pmatrix}, \quad A = \begin{pmatrix} 1 & 2 & -1 & 0 \\ 3 & 1 & 0 & -1 \end{pmatrix}, \quad b = \begin{pmatrix} 4 \\ 5 \end{pmatrix}. \]

问题变为标准形式：

\[\begin{aligned} \min \quad & c^T x \\ \text{s.t.} \quad & Ax = b, \\ & x \ge 0. \end{aligned} \]

步骤2：写出对偶问题及最优性条件

对偶问题为：

\[\begin{aligned} \max \quad & b^T y \\ \text{s.t.} \quad & A^T y + s = c, \\ & s \ge 0. \end{aligned} \]

其中 \(y \in \mathbb{R}^m\) 是对偶变量，\(s \in \mathbb{R}^n\) 是对偶松弛变量。

Karush-Kuhn-Tucker (KKT) 最优性条件（对于线性规划，也是充要条件）为：

原始可行性：\(Ax = b, \quad x \ge 0\).
对偶可行性：\(A^T y + s = c, \quad s \ge 0\).
互补松弛：\(x_i s_i = 0, \quad i = 1, \dots, n\).

步骤3：原始-对偶内点法的基本思路

内点法通过在可行域内部迭代，逼近边界上的最优解。原始-对偶内点法引入一个“中心路径”参数 \(\mu > 0\)，将互补松弛条件松弛为：

\[x_i s_i = \mu, \quad i = 1, \dots, n. \]

这称为“中心路径条件”。当 \(\mu \to 0\) 时，满足该条件的点序列趋近于最优解。

因此，对于给定的 \(\mu > 0\)，我们求解非线性方程组：

\[F_\mu(x, y, s) = \begin{pmatrix} Ax - b \\ A^T y + s - c \\ XSe - \mu e \end{pmatrix} = 0, \]

其中：

\(X = \operatorname{diag}(x_1, \dots, x_n)\) 是以 \(x\) 为对角元素的对角矩阵。
\(S = \operatorname{diag}(s_1, \dots, s_n)\) 是以 \(s\) 为对角元素的对角矩阵。
\(e = (1,1,\dots,1)^T \in \mathbb{R}^n\).
\(XSe\) 表示向量 \((x_1 s_1, x_2 s_2, \dots, x_n s_n)^T\)。

原始-对偶内点法通过牛顿法求解上述方程组，并逐渐减小 \(\mu\)。

步骤4：牛顿步的推导

在当前点 \((x^k, y^k, s^k)\)，我们希望对 \(F_\mu = 0\) 进行线性化，以得到搜索方向 \((\Delta x, \Delta y, \Delta s)\)。

将 \(F_\mu = 0\) 在 \((x^k, y^k, s^k)\) 处线性化：

\[\begin{pmatrix} A & 0 & 0 \\ 0 & A^T & I \\ S^k & 0 & X^k \end{pmatrix} \begin{pmatrix} \Delta x \\ \Delta y \\ \Delta s \end{pmatrix} = - \begin{pmatrix} Ax^k - b \\ A^T y^k + s^k - c \\ X^k S^k e - \mu e \end{pmatrix}. \]

这个系统称为“牛顿方程”。

通常，我们将目标参数 \(\mu\) 取为 \(\mu^k = \sigma \cdot \frac{(x^k)^T s^k}{n}\)，其中 \(\sigma \in [0,1]\) 称为“中心参数”，\((x^k)^T s^k\) 是当前对偶间隙的估计。当 \(\sigma = 0\) 时，我们直接向最优解方向推进；当 \(\sigma = 1\) 时，我们向中心路径方向推进，以保持在可行域内部。

步骤5：具体求解牛顿方向

从牛顿方程中，我们可以解出方向。由第三行方程：

\[S^k \Delta x + X^k \Delta s = -X^k S^k e + \mu e. \]

令 \(r_b = Ax^k - b\)（原始残差），\(r_c = A^T y^k + s^k - c\)（对偶残差），\(r_{xs} = -X^k S^k e + \mu e\)。

我们可以先消去 \(\Delta s\)。从第三行方程解出：

\[\Delta s = X^{-1} (r_{xs} - S \Delta x) = -S \Delta x + X^{-1} r_{xs}. \]

但更常用的方式是直接从原方程推导。实际上，标准做法是：从第一、二行方程和第三行方程联立，得到关于 \(\Delta x\) 和 \(\Delta y\) 的方程。

具体推导（可跳过，记住结果）：
由第二行方程：\(A^T \Delta y + \Delta s = -r_c\)。
将第三行方程写为：\(S \Delta x + X \Delta s = r_{xs}\)。
从第三行解出：\(\Delta s = X^{-1} (r_{xs} - S \Delta x)\)。
代入第二行方程：

\[A^T \Delta y + X^{-1} (r_{xs} - S \Delta x) = -r_c. \]

整理：

\[A^T \Delta y - X^{-1} S \Delta x = -r_c - X^{-1} r_{xs}. \]

注意到 \(X^{-1} S = D^{-2}\)，其中 \(D = X^{1/2} S^{-1/2}\)，但更直接的方法是乘以 \(X S^{-1}\) 来消去。实际上，通常构造正规方程：

从第三行：\(\Delta s = -S \Delta x + X^{-1} r_{xs}\)（不精确，需验证维度，实际上应解为 \(\Delta s = -S X^{-1} X \Delta x + X^{-1} r_{xs} = -S X^{-1} (X \Delta x) + X^{-1} r_{xs}\)）。更标准的方法是引入缩放变量。但为简单起见，我们采用常用结论：求解以下方程得到 \(\Delta y\)：

\[A D^2 A^T \Delta y = p, \]

其中 \(D^2 = X S^{-1}\)，\(p\) 是某个向量。然后再回代求 \(\Delta x\) 和 \(\Delta s\)。

在实际的原始-对偶内点法实现中，我们通常直接解下面的对称正定系统（通过Cholesky分解等）：

\[(A D^2 A^T) \Delta y = A D^2 (r_c - X^{-1} r_{xs}) + r_b, \]

其中 \(D^2 = X S^{-1}\)。

然后计算：

\[\Delta x = D^2 (A^T \Delta y - r_c + X^{-1} r_{xs}), \]

\[ \Delta s = -r_c - A^T \Delta y. \]

步骤6：迭代步骤与算法框架

初始化：选择初始点 \((x^0, y^0, s^0)\) 满足 \(x^0 > 0, s^0 > 0\)（严格内点）。通常可以采用一些启发式方法找到初始内点，或者使用“初始点策略”如“人工变量法”。
迭代：对于 \(k = 0, 1, 2, \dots\)，执行：
a. 计算对偶间隙 \(\mu^k = \frac{(x^k)^T s^k}{n}\)。
b. 如果 \(\mu^k < \epsilon\)（预设精度，例如 \(10^{-8}\)），则停止，当前点近似最优。
c. 选择中心参数 \(\sigma_k \in [0,1]\)，例如 \(\sigma_k = 0.1\)。
d. 计算牛顿方向 \((\Delta x^k, \Delta y^k, \Delta s^k)\) 对应 \(\mu = \sigma_k \mu^k\)。
e. 计算步长 \(\alpha^k\)：为了保持 \(x > 0, s > 0\)，我们取：

\[ \alpha_{\text{pri}} = \min\left(1, \, \gamma \cdot \min_{\Delta x_i^k < 0} \frac{x_i^k}{-\Delta x_i^k} \right), \]

\[ \alpha_{\text{dual}} = \min\left(1, \, \gamma \cdot \min_{\Delta s_i^k < 0} \frac{s_i^k}{-\Delta s_i^k} \right), \]

  其中 $ \gamma \in (0,1) $，例如 0.995，是为了避免直接走到边界。

f. 更新：

\[ x^{k+1} = x^k + \alpha_{\text{pri}} \Delta x^k, \quad y^{k+1} = y^k + \alpha_{\text{dual}} \Delta y^k, \quad s^{k+1} = s^k + \alpha_{\text{dual}} \Delta s^k. \]

输出近似最优解 \((x^*, y^*, s^*)\)。

步骤7：应用于我们的例子

我们手工模拟几步，以展示过程。由于手工计算较繁琐，我们只展示第一轮迭代的概念，具体数值会适当简化。

给定问题：

\[A = \begin{pmatrix} 1 & 2 & -1 & 0 \\ 3 & 1 & 0 & -1 \end{pmatrix}, \quad b = \begin{pmatrix} 4 \\ 5 \end{pmatrix}, \quad c = (2, 3, 0, 0)^T. \]

步骤7.1：初始化

我们需要一个严格内点 \((x^0, y^0, s^0)\) 满足 \(x^0 > 0, s^0 > 0, Ax^0 = b, A^T y^0 + s^0 = c\)。通常可以通过解一个辅助问题得到，但这里我们“猜测”一个内点。注意，我们需要 \(Ax = b\) 和 \(A^T y + s = c\) 同时成立，并且 \(x, s > 0\)。

我们尝试找 \(x^0\)：从约束 \(x_1 + 2x_2 - x_3 = 4, \, 3x_1 + x_2 - x_4 = 5\)。取 \(x_1=2, x_2=2\)，则 \(x_3 = 2+4-4=2\)，\(x_4 = 6+2-5=3\)。所以 \(x^0 = (2,2,2,3)^T > 0\)，且满足 \(Ax^0 = b\)。

再找对偶可行内点：需要 \(A^T y + s = c\) 且 \(s > 0\)。取 \(y = (0,0)^T\)，则 \(s = c - A^T y = c = (2,3,0,0)^T\)。但此时 \(s_3 = 0, s_4 = 0\)，不满足严格大于零。所以我们调整 \(y\)。我们希望 \(s > 0\)，即

\[\begin{pmatrix} 1 & 3 \\ 2 & 1 \\ -1 & 0 \\ 0 & -1 \end{pmatrix} \begin{pmatrix} y_1 \\ y_2 \end{pmatrix} + s = \begin{pmatrix} 2 \\ 3 \\ 0 \\ 0 \end{pmatrix}. \]

取 \(y_1 = -1, y_2 = -1\)，则

\[A^T y = \begin{pmatrix} 1\cdot(-1)+3\cdot(-1) \\ 2\cdot(-1)+1\cdot(-1) \\ -1\cdot(-1)+0\cdot(-1) \\ 0\cdot(-1)+(-1)\cdot(-1) \end{pmatrix} = \begin{pmatrix} -4 \\ -3 \\ 1 \\ 1 \end{pmatrix}. \]

于是 \(s = c - A^T y = (2,3,0,0)^T - (-4,-3,1,1)^T = (6, 6, -1, -1)^T\)，不满足 \(s > 0\)（后两个分量为负）。看来手工构造严格对偶可行内点不太容易。实际上，原始-对偶内点法在实现时通常采用“不可行启动”版本，即允许初始点不满足等式约束，但要求 \(x>0, s>0\)。我们这里为简化，假设我们有一个初始内点（可能通过求解一个扩大问题的技巧得到），具体数值略过，我们直接假设：

\[x^0 = (2, 2, 2, 3)^T, \quad y^0 = (0.5, 0.5)^T, \quad s^0 = (1, 1, 1, 1)^T. \]

验证：\(Ax^0 = (1*2+2*2-2, 3*2+1*2-3) = (4,5) = b\)，满足原始可行。对偶可行性：\(A^T y^0 = (1*0.5+3*0.5, 2*0.5+1*0.5, -1*0.5+0*0.5, 0*0.5-1*0.5) = (2, 1.5, -0.5, -0.5)^T\)。则 \(s^0 = (1,1,1,1)^T\)，于是 \(A^T y^0 + s^0 = (3, 2.5, 0.5, 0.5)^T\)，而 \(c = (2,3,0,0)^T\)，不相等。所以这个点对偶不可行。但这不影响演示算法步骤，因为不可行原始-对偶内点法可以处理。我们仅演示一步计算。

步骤7.2：计算残差和对偶间隙

当前点：

\[x^0 = (2,2,2,3)^T, \quad y^0 = (0.5,0.5)^T, \quad s^0 = (1,1,1,1)^T. \]

计算：
原始残差 \(r_b = Ax^0 - b = 0\)（因为原始可行）。
对偶残差 \(r_c = A^T y^0 + s^0 - c = (3, 2.5, 0.5, 0.5)^T - (2,3,0,0)^T = (1, -0.5, 0.5, 0.5)^T\)。
对偶间隙 \(\mu^0 = \frac{(x^0)^T s^0}{4} = \frac{2+2+2+3}{4} = \frac{9}{4} = 2.25\)。

步骤7.3：选择中心参数并计算 \(\mu\)

取中心参数 \(\sigma = 0.1\)，则 \(\mu = \sigma \mu^0 = 0.225\)。

步骤7.4：计算 \(r_{xs}\) 和 \(D^2\)

\[X^0 = \operatorname{diag}(2,2,2,3), \quad S^0 = \operatorname{diag}(1,1,1,1). \]

\[ X^0 S^0 e = (2,2,2,3)^T. \]

\[ r_{xs} = -X^0 S^0 e + \mu e = (-2+0.225, -2+0.225, -2+0.225, -3+0.225)^T = (-1.775, -1.775, -1.775, -2.775)^T. \]

\[ D^2 = X^0 (S^0)^{-1} = \operatorname{diag}(2,2,2,3). \]

步骤7.5：计算牛顿方程右端项并求解 \(\Delta y\)

我们需要解：

\[(A D^2 A^T) \Delta y = A D^2 (r_c - (X^0)^{-1} r_{xs}) + r_b. \]

计算 \((X^0)^{-1} r_{xs} = (-1.775/2, -1.775/2, -1.775/2, -2.775/3)^T = (-0.8875, -0.8875, -0.8875, -0.925)^T\)。

然后 \(r_c - (X^0)^{-1} r_{xs} = (1 - (-0.8875), -0.5 - (-0.8875), 0.5 - (-0.8875), 0.5 - (-0.925))^T = (1.8875, 0.3875, 1.3875, 1.425)^T\)。

计算 \(A D^2\)：

\[A D^2 = \begin{pmatrix} 1&2&-1&0 \\ 3&1&0&-1 \end{pmatrix} \operatorname{diag}(2,2,2,3) = \begin{pmatrix} 2&4&-2&0 \\ 6&2&0&-3 \end{pmatrix}. \]

然后 \(A D^2 A^T\)：

\[A D^2 A^T = \begin{pmatrix} 2&4&-2&0 \\ 6&2&0&-3 \end{pmatrix} \begin{pmatrix} 1&3 \\ 2&1 \\ -1&0 \\ 0&-1 \end{pmatrix} = \begin{pmatrix} 2*1+4*2+(-2)*(-1)+0*0 & 2*3+4*1+(-2)*0+0*(-1) \\ 6*1+2*2+0*(-1)+(-3)*0 & 6*3+2*1+0*0+(-3)*(-1) \end{pmatrix} = \begin{pmatrix} 2+8+2 & 6+4 \\ 6+4 & 18+2+3 \end{pmatrix} = \begin{pmatrix} 12 & 10 \\ 10 & 23 \end{pmatrix}. \]

右端项向量：

\[A D^2 (r_c - (X^0)^{-1} r_{xs}) = \begin{pmatrix} 2&4&-2&0 \\ 6&2&0&-3 \end{pmatrix} \begin{pmatrix} 1.8875 \\ 0.3875 \\ 1.3875 \\ 1.425 \end{pmatrix} = \begin{pmatrix} 2*1.8875+4*0.3875+(-2)*1.3875+0*1.425 \\ 6*1.8875+2*0.3875+0*1.3875+(-3)*1.425 \end{pmatrix} = \begin{pmatrix} 3.775+1.55-2.775 \\ 11.325+0.775-4.275 \end{pmatrix} = \begin{pmatrix} 2.55 \\ 7.825 \end{pmatrix}. \]

加上 \(r_b = 0\)，所以右端项为 \((2.55, 7.825)^T\)。

解方程：

\[\begin{pmatrix} 12 & 10 \\ 10 & 23 \end{pmatrix} \begin{pmatrix} \Delta y_1 \\ \Delta y_2 \end{pmatrix} = \begin{pmatrix} 2.55 \\ 7.825 \end{pmatrix}. \]

计算行列式：\(12*23 - 10*10 = 276-100=176\)。逆矩阵：

\[\frac{1}{176} \begin{pmatrix} 23 & -10 \\ -10 & 12 \end{pmatrix}. \]

所以：

\[\Delta y = \frac{1}{176} \begin{pmatrix} 23*2.55 -10*7.825 \\ -10*2.55+12*7.825 \end{pmatrix} = \frac{1}{176} \begin{pmatrix} 58.65-78.25 \\ -25.5+93.9 \end{pmatrix} = \frac{1}{176} \begin{pmatrix} -19.6 \\ 68.4 \end{pmatrix} = \begin{pmatrix} -0.11136 \\ 0.38864 \end{pmatrix} \approx \begin{pmatrix} -0.1114 \\ 0.3886 \end{pmatrix}. \]

步骤7.6：计算 \(\Delta x\) 和 \(\Delta s\)

公式：

\[\Delta x = D^2 (A^T \Delta y - r_c + (X^0)^{-1} r_{xs}). \]

计算 \(A^T \Delta y\)：

\[A^T \Delta y = \begin{pmatrix} 1&3 \\ 2&1 \\ -1&0 \\ 0&-1 \end{pmatrix} \begin{pmatrix} -0.1114 \\ 0.3886 \end{pmatrix} = \begin{pmatrix} -0.1114+1.1658 \\ -0.2228+0.3886 \\ 0.1114+0 \\ 0-0.3886 \end{pmatrix} = \begin{pmatrix} 1.0544 \\ 0.1658 \\ 0.1114 \\ -0.3886 \end{pmatrix}. \]

然后 \(A^T \Delta y - r_c + (X^0)^{-1} r_{xs} = (1.0544, 0.1658, 0.1114, -0.3886)^T - (1, -0.5, 0.5, 0.5)^T + (-0.8875, -0.8875, -0.8875, -0.925)^T = (1.0544-1-0.8875, 0.1658+0.5-0.8875, 0.1114-0.5-0.8875, -0.3886-0.5-0.925) = (-0.8331, -0.2217, -1.2761, -1.8136)^T\)。

乘以 \(D^2 = \operatorname{diag}(2,2,2,3)\)：

\[\Delta x = (2*(-0.8331), 2*(-0.2217), 2*(-1.2761), 3*(-1.8136))^T = (-1.6662, -0.4434, -2.5522, -5.4408)^T. \]

然后计算 \(\Delta s = -r_c - A^T \Delta y\)：

\[\Delta s = -(1, -0.5, 0.5, 0.5)^T - (1.0544, 0.1658, 0.1114, -0.3886)^T = (-2.0544, 0.3342, -0.6114, -0.1114)^T. \]

步骤7.7：计算步长并更新

对于原始步长，检查 \(\Delta x\) 的负分量：\(\Delta x_1 = -1.6662, \Delta x_2 = -0.4434, \Delta x_3 = -2.5522, \Delta x_4 = -5.4408\)，全为负。所以

\[\min_{\Delta x_i < 0} \frac{x_i^0}{-\Delta x_i} = \min\left( \frac{2}{1.6662}, \frac{2}{0.4434}, \frac{2}{2.5522}, \frac{3}{5.4408} \right) = \min(1.200, 4.511, 0.784, 0.551) = 0.551. \]

取 \(\gamma = 0.995\)，则 \(\alpha_{\text{pri}} = \min(1, 0.995*0.551) = 0.548\)。

对偶步长类似，检查 \(\Delta s\) 的负分量：\(\Delta s_1 = -2.0544, \Delta s_3 = -0.6114, \Delta s_4 = -0.1114\)。计算：

\[\min\left( \frac{1}{2.0544}, \frac{1}{0.6114}, \frac{1}{0.1114} \right) = \min(0.487, 1.636, 8.977) = 0.487. \]

则 \(\alpha_{\text{dual}} = \min(1, 0.995*0.487) = 0.484\)。

更新：

\[x^1 = x^0 + \alpha_{\text{pri}} \Delta x = (2,2,2,3)^T + 0.548 * (-1.6662, -0.4434, -2.5522, -5.4408)^T = (2-0.912, 2-0.243, 2-1.398, 3-2.979) = (1.088, 1.757, 0.602, 0.021)^T. \]

\[ y^1 = y^0 + \alpha_{\text{dual}} \Delta y = (0.5,0.5) + 0.484*(-0.1114, 0.3886) = (0.5-0.0539, 0.5+0.1881) = (0.4461, 0.6881). \]

\[ s^1 = s^0 + \alpha_{\text{dual}} \Delta s = (1,1,1,1) + 0.484*(-2.0544, 0.3342, -0.6114, -0.1114) = (1-0.994, 1+0.1617, 1-0.2959, 1-0.0539) = (0.006, 1.1617, 0.7041, 0.9461). \]

步骤8：继续迭代

重复步骤2-7，直到对偶间隙 \(\mu^k\) 足够小。最终，我们将逼近最优解。

通过求解原问题（简单线性规划，可图解或单纯形法），我们可以验证最优解。原问题为：

\[\min 2x_1+3x_2 \quad \text{s.t.} \quad x_1+2x_2 \ge 4, \; 3x_1+x_2 \ge 5, \; x_1,x_2 \ge 0. \]

解两个等式：\(x_1+2x_2=4\) 和 \(3x_1+x_2=5\)。解得交点：从第二个方程 \(x_2=5-3x_1\)，代入第一个：\(x_1+2(5-3x_1)=4 \Rightarrow x_1+10-6x_1=4 \Rightarrow -5x_1=-6 \Rightarrow x_1=1.2, x_2=5-3*1.2=1.4\)。目标值 \(2*1.2+3*1.4=2.4+4.2=6.6\)。这是最优解。对应到标准形式，\(x_3 = x_1+2x_2-4=1.2+2.8-4=0\)，\(x_4=3*1.2+1.4-5=3.6+1.4-5=0\)。所以最优 \(x^* = (1.2, 1.4, 0, 0)^T\)。

对偶问题：最大化 \(4y_1+5y_2\)，约束 \(y_1+3y_2 \le 2, 2y_1+y_2 \le 3, y_1 \ge 0, y_2 \ge 0\)。在最优解处，前两个约束紧（因为 \(x_1, x_2 > 0\)），解 \(y_1+3y_2=2, 2y_1+y_2=3\)，解得 \(y_1=1.4, y_2=0.2\)。对偶松弛变量 \(s = c - A^T y = (2,3,0,0)^T - (1.4+0.6, 2.8+0.2, -1.4, -0.2)^T = (0,0,1.4,0.2)^T\)。

所以，原始-对偶内点法迭代最终会逼近 \(x^* = (1.2,1.4,0,0)^T\)，\(y^* = (1.4,0.2)^T\)，\(s^* = (0,0,1.4,0.2)^T\)，满足互补松弛。

总结

原始-对偶内点法是求解线性规划（包括来自随机规划的确定性等价形式）的有效算法。它通过牛顿法求解松弛的KKT条件，并沿中心路径逼近最优解。本示例展示了该方法在解决一个简单随机线性规划（期望值模型）中的应用，包括问题转化、算法步骤和一轮迭代的详细计算。

线性规划的原始-对偶内点法求解随机线性规划问题示例题目描述考虑以下随机线性规划问题：目标：最小化 \( c^T x \) 约束：\( Ax = b, \quad x \ge 0 \) 其中：部分系数（例如，成本向量 \( c \)、约束矩阵 \( A \) 的某些元素、或右端向量 \( b \) 的某些分量）是随机变量，其概率分布已知。我们希望在满足约束的情况下，最小化“期望总成本”，即最小化 \( \mathbb{E}[ c ]^T x \)，同时要求决策变量 \( x \) 是“此时必须确定、不依赖随机变量实现”的决策（在随机规划中，这被称为“第一阶段”决策或“此时此地”决策）。为了具体化，我们设定一个简单但典型的例子。假设一个生产商需要提前决定两种产品（产品1和产品2）的生产量 \( x_ 1, x_ 2 \)，以应对未来的随机需求。已知每件产品的生产成本是确定的，分别为5和4。但是，市场需求是随机的，且如果产量超过实际需求，会有单位库存持有成本（随机，与需求相关）；如果产量不足，则会有单位缺货惩罚成本（随机）。生产受资源约束。我们想最小化“期望总成本”（生产成本 + 期望库存持有与缺货成本之和）。经过建模，该问题可以表述为一个两阶段随机线性规划，其中第一阶段决策是 \( x_ 1, x_ 2 \)，第二阶段针对每种随机场景优化库存/缺货决策。在特定假设下（例如，随机参数只出现在目标函数的某些成本系数中，且其期望可计算），该问题可简化为一个确定性线性规划，其目标系数为对应随机参数的期望值。为了演示算法，我们使用如下具体的确定性等价线性规划模型：最小化： \[ 2x_ 1 + 3x_ 2 \] 约束于： \[ \begin{aligned} x_ 1 + 2x_ 2 &\ge 4 \quad (\text{需求约束，系数为期望需求}) \\ 3x_ 1 + x_ 2 &\ge 5 \quad (\text{资源约束，期望资源消耗要求}) \\ x_ 1, x_ 2 &\ge 0 \] 注意：这里的系数 2 和 3 已经是考虑了生产、库存、缺货等随机成本后的期望单位成本。约束右端 4 和 5 是考虑了随机性后的期望要求值。我们将使用原始-对偶内点法来求解这个线性规划。原始-对偶内点法是求解线性规划的一类高效内点法，它通过求解由原始可行性、对偶可行性和互补松弛条件组成的非线性方程组，并沿中心路径逼近最优解。解题过程步骤1：将问题转化为标准形式原始-对偶内点法通常处理标准形式的线性规划： \[ \begin{aligned} \min \quad & c^T x \\ \text{s.t.} \quad & Ax = b, \\ & x \ge 0 \end{aligned} \] 其中，\( A \in \mathbb{R}^{m \times n} \)，通常假设 \( A \) 行满秩。我们的例子是： \[ \begin{aligned} \min \quad & 2x_ 1 + 3x_ 2 \\ \text{s.t.} \quad & x_ 1 + 2x_ 2 \ge 4, \\ & 3x_ 1 + x_ 2 \ge 5, \\ & x_ 1, x_ 2 \ge 0. \end{aligned} \] 引入松弛变量 \( x_ 3, x_ 4 \ge 0 \)，将不等式化为等式： \[ \begin{aligned} x_ 1 + 2x_ 2 - x_ 3 &= 4, \\ 3x_ 1 + x_ 2 - x_ 4 &= 5. \end{aligned} \] 现在，令 \[ x = \begin{pmatrix} x_ 1 \\ x_ 2 \\ x_ 3 \\ x_ 4 \end{pmatrix}, \quad c = \begin{pmatrix} 2 \\ 3 \\ 0 \\ 0 \end{pmatrix}, \quad A = \begin{pmatrix} 1 & 2 & -1 & 0 \\ 3 & 1 & 0 & -1 \end{pmatrix}, \quad b = \begin{pmatrix} 4 \\ 5 \end{pmatrix}. \] 问题变为标准形式： \[ \begin{aligned} \min \quad & c^T x \\ \text{s.t.} \quad & Ax = b, \\ & x \ge 0. \end{aligned} \] 步骤2：写出对偶问题及最优性条件对偶问题为： \[ \begin{aligned} \max \quad & b^T y \\ \text{s.t.} \quad & A^T y + s = c, \\ & s \ge 0. \end{aligned} \] 其中 \( y \in \mathbb{R}^m \) 是对偶变量，\( s \in \mathbb{R}^n \) 是对偶松弛变量。 Karush-Kuhn-Tucker (KKT) 最优性条件（对于线性规划，也是充要条件）为：原始可行性：\( Ax = b, \quad x \ge 0 \). 对偶可行性：\( A^T y + s = c, \quad s \ge 0 \). 互补松弛：\( x_ i s_ i = 0, \quad i = 1, \dots, n \). 步骤3：原始-对偶内点法的基本思路内点法通过在可行域内部迭代，逼近边界上的最优解。原始-对偶内点法引入一个“中心路径”参数 \( \mu > 0 \)，将互补松弛条件松弛为： \[ x_ i s_ i = \mu, \quad i = 1, \dots, n. \] 这称为“中心路径条件”。当 \( \mu \to 0 \) 时，满足该条件的点序列趋近于最优解。因此，对于给定的 \( \mu > 0 \)，我们求解非线性方程组： \[ F_ \mu(x, y, s) = \begin{pmatrix} Ax - b \\ A^T y + s - c \\ XSe - \mu e \end{pmatrix} = 0, \] 其中： \( X = \operatorname{diag}(x_ 1, \dots, x_ n) \) 是以 \( x \) 为对角元素的对角矩阵。 \( S = \operatorname{diag}(s_ 1, \dots, s_ n) \) 是以 \( s \) 为对角元素的对角矩阵。 \( e = (1,1,\dots,1)^T \in \mathbb{R}^n \). \( XSe \) 表示向量 \( (x_ 1 s_ 1, x_ 2 s_ 2, \dots, x_ n s_ n)^T \)。原始-对偶内点法通过牛顿法求解上述方程组，并逐渐减小 \( \mu \)。步骤4：牛顿步的推导在当前点 \( (x^k, y^k, s^k) \)，我们希望对 \( F_ \mu = 0 \) 进行线性化，以得到搜索方向 \( (\Delta x, \Delta y, \Delta s) \)。将 \( F_ \mu = 0 \) 在 \( (x^k, y^k, s^k) \) 处线性化： \[ \begin{pmatrix} A & 0 & 0 \\ 0 & A^T & I \\ S^k & 0 & X^k \end{pmatrix} \begin{pmatrix} \Delta x \\ \Delta y \\ \Delta s \end{pmatrix} \begin{pmatrix} Ax^k - b \\ A^T y^k + s^k - c \\ X^k S^k e - \mu e \end{pmatrix}. \] 这个系统称为“牛顿方程”。通常，我们将目标参数 \( \mu \) 取为 \( \mu^k = \sigma \cdot \frac{(x^k)^T s^k}{n} \)，其中 \( \sigma \in [ 0,1 ] \) 称为“中心参数”，\( (x^k)^T s^k \) 是当前对偶间隙的估计。当 \( \sigma = 0 \) 时，我们直接向最优解方向推进；当 \( \sigma = 1 \) 时，我们向中心路径方向推进，以保持在可行域内部。步骤5：具体求解牛顿方向从牛顿方程中，我们可以解出方向。由第三行方程： \[ S^k \Delta x + X^k \Delta s = -X^k S^k e + \mu e. \] 令 \( r_ b = Ax^k - b \)（原始残差），\( r_ c = A^T y^k + s^k - c \)（对偶残差），\( r_ {xs} = -X^k S^k e + \mu e \)。我们可以先消去 \( \Delta s \)。从第三行方程解出： \[ \Delta s = X^{-1} (r_ {xs} - S \Delta x) = -S \Delta x + X^{-1} r_ {xs}. \] 但更常用的方式是直接从原方程推导。实际上，标准做法是：从第一、二行方程和第三行方程联立，得到关于 \( \Delta x \) 和 \( \Delta y \) 的方程。具体推导（可跳过，记住结果）：由第二行方程：\( A^T \Delta y + \Delta s = -r_ c \)。将第三行方程写为：\( S \Delta x + X \Delta s = r_ {xs} \)。从第三行解出：\( \Delta s = X^{-1} (r_ {xs} - S \Delta x) \)。代入第二行方程： \[ A^T \Delta y + X^{-1} (r_ {xs} - S \Delta x) = -r_ c. \] 整理： \[ A^T \Delta y - X^{-1} S \Delta x = -r_ c - X^{-1} r_ {xs}. \] 注意到 \( X^{-1} S = D^{-2} \)，其中 \( D = X^{1/2} S^{-1/2} \)，但更直接的方法是乘以 \( X S^{-1} \) 来消去。实际上，通常构造正规方程：从第三行：\( \Delta s = -S \Delta x + X^{-1} r_ {xs} \)（不精确，需验证维度，实际上应解为 \( \Delta s = -S X^{-1} X \Delta x + X^{-1} r_ {xs} = -S X^{-1} (X \Delta x) + X^{-1} r_ {xs} \)）。更标准的方法是引入缩放变量。但为简单起见，我们采用常用结论：求解以下方程得到 \( \Delta y \)： \[ A D^2 A^T \Delta y = p, \] 其中 \( D^2 = X S^{-1} \)，\( p \) 是某个向量。然后再回代求 \( \Delta x \) 和 \( \Delta s \)。在实际的原始-对偶内点法实现中，我们通常直接解下面的对称正定系统（通过Cholesky分解等）： \[ (A D^2 A^T) \Delta y = A D^2 (r_ c - X^{-1} r_ {xs}) + r_ b, \] 其中 \( D^2 = X S^{-1} \)。然后计算： \[ \Delta x = D^2 (A^T \Delta y - r_ c + X^{-1} r_ {xs}), \] \[ \Delta s = -r_ c - A^T \Delta y. \] 步骤6：迭代步骤与算法框架初始化：选择初始点 \( (x^0, y^0, s^0) \) 满足 \( x^0 > 0, s^0 > 0 \)（严格内点）。通常可以采用一些启发式方法找到初始内点，或者使用“初始点策略”如“人工变量法”。迭代：对于 \( k = 0, 1, 2, \dots \)，执行： a. 计算对偶间隙 \( \mu^k = \frac{(x^k)^T s^k}{n} \)。 b. 如果 \( \mu^k < \epsilon \)（预设精度，例如 \( 10^{-8} \)），则停止，当前点近似最优。 c. 选择中心参数 \( \sigma_ k \in [ 0,1] \)，例如 \( \sigma_ k = 0.1 \)。 d. 计算牛顿方向 \( (\Delta x^k, \Delta y^k, \Delta s^k) \) 对应 \( \mu = \sigma_ k \mu^k \)。 e. 计算步长 \( \alpha^k \)：为了保持 \( x > 0, s > 0 \)，我们取： \[ \alpha_ {\text{pri}} = \min\left(1, \, \gamma \cdot \min_ {\Delta x_ i^k < 0} \frac{x_ i^k}{-\Delta x_ i^k} \right), \] \[ \alpha_ {\text{dual}} = \min\left(1, \, \gamma \cdot \min_ {\Delta s_ i^k < 0} \frac{s_ i^k}{-\Delta s_ i^k} \right), \] 其中 \( \gamma \in (0,1) \)，例如 0.995，是为了避免直接走到边界。 f. 更新： \[ x^{k+1} = x^k + \alpha_ {\text{pri}} \Delta x^k, \quad y^{k+1} = y^k + \alpha_ {\text{dual}} \Delta y^k, \quad s^{k+1} = s^k + \alpha_ {\text{dual}} \Delta s^k. \] 输出近似最优解 \( (x^ , y^ , s^* ) \)。步骤7：应用于我们的例子我们手工模拟几步，以展示过程。由于手工计算较繁琐，我们只展示第一轮迭代的概念，具体数值会适当简化。给定问题： \[ A = \begin{pmatrix} 1 & 2 & -1 & 0 \\ 3 & 1 & 0 & -1 \end{pmatrix}, \quad b = \begin{pmatrix} 4 \\ 5 \end{pmatrix}, \quad c = (2, 3, 0, 0)^T. \] 步骤7.1：初始化我们需要一个严格内点 \( (x^0, y^0, s^0) \) 满足 \( x^0 > 0, s^0 > 0, Ax^0 = b, A^T y^0 + s^0 = c \)。通常可以通过解一个辅助问题得到，但这里我们“猜测”一个内点。注意，我们需要 \( Ax = b \) 和 \( A^T y + s = c \) 同时成立，并且 \( x, s > 0 \)。我们尝试找 \( x^0 \)：从约束 \( x_ 1 + 2x_ 2 - x_ 3 = 4, \, 3x_ 1 + x_ 2 - x_ 4 = 5 \)。取 \( x_ 1=2, x_ 2=2 \)，则 \( x_ 3 = 2+4-4=2 \)，\( x_ 4 = 6+2-5=3 \)。所以 \( x^0 = (2,2,2,3)^T > 0 \)，且满足 \( Ax^0 = b \)。再找对偶可行内点：需要 \( A^T y + s = c \) 且 \( s > 0 \)。取 \( y = (0,0)^T \)，则 \( s = c - A^T y = c = (2,3,0,0)^T \)。但此时 \( s_ 3 = 0, s_ 4 = 0 \)，不满足严格大于零。所以我们调整 \( y \)。我们希望 \( s > 0 \)，即 \[ \begin{pmatrix} 1 & 3 \\ 2 & 1 \\ -1 & 0 \\ 0 & -1 \end{pmatrix} \begin{pmatrix} y_ 1 \\ y_ 2 \end{pmatrix} + s = \begin{pmatrix} 2 \\ 3 \\ 0 \\ 0 \end{pmatrix}. \] 取 \( y_ 1 = -1, y_ 2 = -1 \)，则 \[ A^T y = \begin{pmatrix} 1\cdot(-1)+3\cdot(-1) \\ 2\cdot(-1)+1\cdot(-1) \\ -1\cdot(-1)+0\cdot(-1) \\ 0\cdot(-1)+(-1)\cdot(-1) \end{pmatrix} = \begin{pmatrix} -4 \\ -3 \\ 1 \\ 1 \end{pmatrix}. \] 于是 \( s = c - A^T y = (2,3,0,0)^T - (-4,-3,1,1)^T = (6, 6, -1, -1)^T \)，不满足 \( s > 0 \)（后两个分量为负）。看来手工构造严格对偶可行内点不太容易。实际上，原始-对偶内点法在实现时通常采用“不可行启动”版本，即允许初始点不满足等式约束，但要求 \( x>0, s>0 \)。我们这里为简化，假设我们有一个初始内点（可能通过求解一个扩大问题的技巧得到），具体数值略过，我们直接假设： \[ x^0 = (2, 2, 2, 3)^T, \quad y^0 = (0.5, 0.5)^T, \quad s^0 = (1, 1, 1, 1)^T. \] 验证：\( Ax^0 = (1 2+2 2-2, 3 2+1 2-3) = (4,5) = b \)，满足原始可行。对偶可行性：\( A^T y^0 = (1 0.5+3 0.5, 2 0.5+1 0.5, -1 0.5+0 0.5, 0 0.5-1 0.5) = (2, 1.5, -0.5, -0.5)^T \)。则 \( s^0 = (1,1,1,1)^T \)，于是 \( A^T y^0 + s^0 = (3, 2.5, 0.5, 0.5)^T \)，而 \( c = (2,3,0,0)^T \)，不相等。所以这个点对偶不可行。但这不影响演示算法步骤，因为不可行原始-对偶内点法可以处理。我们仅演示一步计算。步骤7.2：计算残差和对偶间隙当前点： \[ x^0 = (2,2,2,3)^T, \quad y^0 = (0.5,0.5)^T, \quad s^0 = (1,1,1,1)^T. \] 计算：原始残差 \( r_ b = Ax^0 - b = 0 \)（因为原始可行）。对偶残差 \( r_ c = A^T y^0 + s^0 - c = (3, 2.5, 0.5, 0.5)^T - (2,3,0,0)^T = (1, -0.5, 0.5, 0.5)^T \)。对偶间隙 \( \mu^0 = \frac{(x^0)^T s^0}{4} = \frac{2+2+2+3}{4} = \frac{9}{4} = 2.25 \)。步骤7.3：选择中心参数并计算 \( \mu \) 取中心参数 \( \sigma = 0.1 \)，则 \( \mu = \sigma \mu^0 = 0.225 \)。步骤7.4：计算 \( r_ {xs} \) 和 \( D^2 \) \[ X^0 = \operatorname{diag}(2,2,2,3), \quad S^0 = \operatorname{diag}(1,1,1,1). \] \[ X^0 S^0 e = (2,2,2,3)^T. \] \[ r_ {xs} = -X^0 S^0 e + \mu e = (-2+0.225, -2+0.225, -2+0.225, -3+0.225)^T = (-1.775, -1.775, -1.775, -2.775)^T. \] \[ D^2 = X^0 (S^0)^{-1} = \operatorname{diag}(2,2,2,3). \] 步骤7.5：计算牛顿方程右端项并求解 \( \Delta y \) 我们需要解： \[ (A D^2 A^T) \Delta y = A D^2 (r_ c - (X^0)^{-1} r_ {xs}) + r_ b. \] 计算 \( (X^0)^{-1} r_ {xs} = (-1.775/2, -1.775/2, -1.775/2, -2.775/3)^T = (-0.8875, -0.8875, -0.8875, -0.925)^T \)。然后 \( r_ c - (X^0)^{-1} r_ {xs} = (1 - (-0.8875), -0.5 - (-0.8875), 0.5 - (-0.8875), 0.5 - (-0.925))^T = (1.8875, 0.3875, 1.3875, 1.425)^T \)。计算 \( A D^2 \)： \[ A D^2 = \begin{pmatrix} 1&2&-1&0 \\ 3&1&0&-1 \end{pmatrix} \operatorname{diag}(2,2,2,3) = \begin{pmatrix} 2&4&-2&0 \\ 6&2&0&-3 \end{pmatrix}. \] 然后 \( A D^2 A^T \)： \[ A D^2 A^T = \begin{pmatrix} 2&4&-2&0 \\ 6&2&0&-3 \end{pmatrix} \begin{pmatrix} 1&3 \\ 2&1 \\ -1&0 \\ 0&-1 \end{pmatrix} = \begin{pmatrix} 2 1+4 2+(-2) (-1)+0 0 & 2 3+4 1+(-2) 0+0 (-1) \\ 6 1+2 2+0* (-1)+(-3) 0 & 6 3+2 1+0 0+(-3) (-1) \end{pmatrix} = \begin{pmatrix} 2+8+2 & 6+4 \\ 6+4 & 18+2+3 \end{pmatrix} = \begin{pmatrix} 12 & 10 \\ 10 & 23 \end{pmatrix}. \] 右端项向量： \[ A D^2 (r_ c - (X^0)^{-1} r_ {xs}) = \begin{pmatrix} 2&4&-2&0 \\ 6&2&0&-3 \end{pmatrix} \begin{pmatrix} 1.8875 \\ 0.3875 \\ 1.3875 \\ 1.425 \end{pmatrix} = \begin{pmatrix} 2 1.8875+4 0.3875+(-2) 1.3875+0 1.425 \\ 6 1.8875+2 0.3875+0 1.3875+(-3)* 1.425 \end{pmatrix} = \begin{pmatrix} 3.775+1.55-2.775 \\ 11.325+0.775-4.275 \end{pmatrix} = \begin{pmatrix} 2.55 \\ 7.825 \end{pmatrix}. \] 加上 \( r_ b = 0 \)，所以右端项为 \( (2.55, 7.825)^T \)。解方程： \[ \begin{pmatrix} 12 & 10 \\ 10 & 23 \end{pmatrix} \begin{pmatrix} \Delta y_ 1 \\ \Delta y_ 2 \end{pmatrix} = \begin{pmatrix} 2.55 \\ 7.825 \end{pmatrix}. \] 计算行列式：\( 12 23 - 10 10 = 276-100=176 \)。逆矩阵： \[ \frac{1}{176} \begin{pmatrix} 23 & -10 \\ -10 & 12 \end{pmatrix}. \] 所以： \[ \Delta y = \frac{1}{176} \begin{pmatrix} 23 2.55 -10 7.825 \\ -10 2.55+12 7.825 \end{pmatrix} = \frac{1}{176} \begin{pmatrix} 58.65-78.25 \\ -25.5+93.9 \end{pmatrix} = \frac{1}{176} \begin{pmatrix} -19.6 \\ 68.4 \end{pmatrix} = \begin{pmatrix} -0.11136 \\ 0.38864 \end{pmatrix} \approx \begin{pmatrix} -0.1114 \\ 0.3886 \end{pmatrix}. \] 步骤7.6：计算 \( \Delta x \) 和 \( \Delta s \) 公式： \[ \Delta x = D^2 (A^T \Delta y - r_ c + (X^0)^{-1} r_ {xs}). \] 计算 \( A^T \Delta y \)： \[ A^T \Delta y = \begin{pmatrix} 1&3 \\ 2&1 \\ -1&0 \\ 0&-1 \end{pmatrix} \begin{pmatrix} -0.1114 \\ 0.3886 \end{pmatrix} = \begin{pmatrix} -0.1114+1.1658 \\ -0.2228+0.3886 \\ 0.1114+0 \\ 0-0.3886 \end{pmatrix} = \begin{pmatrix} 1.0544 \\ 0.1658 \\ 0.1114 \\ -0.3886 \end{pmatrix}. \] 然后 \( A^T \Delta y - r_ c + (X^0)^{-1} r_ {xs} = (1.0544, 0.1658, 0.1114, -0.3886)^T - (1, -0.5, 0.5, 0.5)^T + (-0.8875, -0.8875, -0.8875, -0.925)^T = (1.0544-1-0.8875, 0.1658+0.5-0.8875, 0.1114-0.5-0.8875, -0.3886-0.5-0.925) = (-0.8331, -0.2217, -1.2761, -1.8136)^T \)。乘以 \( D^2 = \operatorname{diag}(2,2,2,3) \)： \[ \Delta x = (2* (-0.8331), 2* (-0.2217), 2* (-1.2761), 3* (-1.8136))^T = (-1.6662, -0.4434, -2.5522, -5.4408)^T. \] 然后计算 \( \Delta s = -r_ c - A^T \Delta y \)： \[ \Delta s = -(1, -0.5, 0.5, 0.5)^T - (1.0544, 0.1658, 0.1114, -0.3886)^T = (-2.0544, 0.3342, -0.6114, -0.1114)^T. \] 步骤7.7：计算步长并更新对于原始步长，检查 \( \Delta x \) 的负分量：\( \Delta x_ 1 = -1.6662, \Delta x_ 2 = -0.4434, \Delta x_ 3 = -2.5522, \Delta x_ 4 = -5.4408 \)，全为负。所以 \[ \min_ {\Delta x_ i < 0} \frac{x_ i^0}{-\Delta x_ i} = \min\left( \frac{2}{1.6662}, \frac{2}{0.4434}, \frac{2}{2.5522}, \frac{3}{5.4408} \right) = \min(1.200, 4.511, 0.784, 0.551) = 0.551. \] 取 \( \gamma = 0.995 \)，则 \( \alpha_ {\text{pri}} = \min(1, 0.995* 0.551) = 0.548 \)。对偶步长类似，检查 \( \Delta s \) 的负分量：\( \Delta s_ 1 = -2.0544, \Delta s_ 3 = -0.6114, \Delta s_ 4 = -0.1114 \)。计算： \[ \min\left( \frac{1}{2.0544}, \frac{1}{0.6114}, \frac{1}{0.1114} \right) = \min(0.487, 1.636, 8.977) = 0.487. \] 则 \( \alpha_ {\text{dual}} = \min(1, 0.995* 0.487) = 0.484 \)。更新： \[ x^1 = x^0 + \alpha_ {\text{pri}} \Delta x = (2,2,2,3)^T + 0.548 * (-1.6662, -0.4434, -2.5522, -5.4408)^T = (2-0.912, 2-0.243, 2-1.398, 3-2.979) = (1.088, 1.757, 0.602, 0.021)^T. \] \[ y^1 = y^0 + \alpha_ {\text{dual}} \Delta y = (0.5,0.5) + 0.484* (-0.1114, 0.3886) = (0.5-0.0539, 0.5+0.1881) = (0.4461, 0.6881). \] \[ s^1 = s^0 + \alpha_ {\text{dual}} \Delta s = (1,1,1,1) + 0.484* (-2.0544, 0.3342, -0.6114, -0.1114) = (1-0.994, 1+0.1617, 1-0.2959, 1-0.0539) = (0.006, 1.1617, 0.7041, 0.9461). \] 步骤8：继续迭代重复步骤2-7，直到对偶间隙 \( \mu^k \) 足够小。最终，我们将逼近最优解。通过求解原问题（简单线性规划，可图解或单纯形法），我们可以验证最优解。原问题为： \[ \min 2x_ 1+3x_ 2 \quad \text{s.t.} \quad x_ 1+2x_ 2 \ge 4, \; 3x_ 1+x_ 2 \ge 5, \; x_ 1,x_ 2 \ge 0. \] 解两个等式：\( x_ 1+2x_ 2=4 \) 和 \( 3x_ 1+x_ 2=5 \)。解得交点：从第二个方程 \( x_ 2=5-3x_ 1 \)，代入第一个：\( x_ 1+2(5-3x_ 1)=4 \Rightarrow x_ 1+10-6x_ 1=4 \Rightarrow -5x_ 1=-6 \Rightarrow x_ 1=1.2, x_ 2=5-3 1.2=1.4 \)。目标值 \( 2 1.2+3 1.4=2.4+4.2=6.6 \)。这是最优解。对应到标准形式，\( x_ 3 = x_ 1+2x_ 2-4=1.2+2.8-4=0 \)，\( x_ 4=3 1.2+1.4-5=3.6+1.4-5=0 \)。所以最优 \( x^* = (1.2, 1.4, 0, 0)^T \)。对偶问题：最大化 \( 4y_ 1+5y_ 2 \)，约束 \( y_ 1+3y_ 2 \le 2, 2y_ 1+y_ 2 \le 3, y_ 1 \ge 0, y_ 2 \ge 0 \)。在最优解处，前两个约束紧（因为 \( x_ 1, x_ 2 > 0 \)），解 \( y_ 1+3y_ 2=2, 2y_ 1+y_ 2=3 \)，解得 \( y_ 1=1.4, y_ 2=0.2 \)。对偶松弛变量 \( s = c - A^T y = (2,3,0,0)^T - (1.4+0.6, 2.8+0.2, -1.4, -0.2)^T = (0,0,1.4,0.2)^T \)。所以，原始-对偶内点法迭代最终会逼近 \( x^* = (1.2,1.4,0,0)^T \)，\( y^* = (1.4,0.2)^T \)，\( s^* = (0,0,1.4,0.2)^T \)，满足互补松弛。总结原始-对偶内点法是求解线性规划（包括来自随机规划的确定性等价形式）的有效算法。它通过牛顿法求解松弛的KKT条件，并沿中心路径逼近最优解。本示例展示了该方法在解决一个简单随机线性规划（期望值模型）中的应用，包括问题转化、算法步骤和一轮迭代的详细计算。