序列仿射缩放内点法基础题

题目描述：
考虑非线性规划问题：

\[\begin{aligned} \min_{x \in \mathbb{R}^2} \quad & f(x) = (x_1 - 3)^2 + (x_2 - 2)^2 \\ \text{s.t.} \quad & x_1 + 2x_2 \leq 6 \\ & x_1 \geq 0, \quad x_2 \geq 0 \end{aligned} \]

我们要求用序列仿射缩放内点法求解。初始点取可行内点 \(x^{(0)} = (1, 2)^T\)，初始对偶变量（拉格朗日乘子）取 \(y^{(0)} = 1\)（对应不等式约束），初始障碍参数 \(\mu^{(0)} = 1\)，缩放参数 \(\sigma = 0.5\)，迭代至多2步，并给出每步的迭代点、对偶变量、障碍参数和近似目标函数值。

解题过程循序渐进讲解

1. 理解问题与算法思想
   序列仿射缩放内点法是一种求解线性或非线性凸规划的内点法。其核心思想是：

   • 在目标函数中加入对数障碍函数，将不等式约束 \(x \geq 0\) 和 \(g(x) \leq 0\) 转化为可微的惩罚项，形成障碍问题。
   • 通过迭代减小障碍参数 \(\mu\)，使得障碍问题的解逐渐逼近原问题的最优解。
   • 每一步中，利用当前点处的仿射缩放变换（即用决策变量的倒数构成缩放矩阵），将问题转化为一个缩放空间中的线性或二次近似子问题，求解该子问题得到搜索方向。
   • 沿搜索方向进行线搜索得到新迭代点，并更新障碍参数和对偶变量。

2. 写出障碍问题形式
   将不等式约束 \(x_1 + 2x_2 \leq 6\) 改写为 \(g(x) = x_1 + 2x_2 - 6 \leq 0\)，并引入松弛变量 \(s = -g(x) \geq 0\)。由于 \(x_1, x_2 \geq 0\)，直接使用对数障碍。障碍问题为：

\[ \min_{x} \quad B(x; \mu) = f(x) - \mu \left[ \ln x_1 + \ln x_2 + \ln(6 - x_1 - 2x_2) \right] \]

其中 \(\mu > 0\) 是障碍参数。当 \(\mu \to 0\) 时，\(B(x; \mu)\) 的最优解趋近于原问题最优解。

3. 构造迭代步的二次近似子问题
   在第 \(k\) 步迭代，当前点为 \(x^{(k)}\)，当前障碍参数为 \(\mu^{(k)}\)，当前对偶变量估计为 \(y^{(k)}\)（对应约束 \(g(x) \leq 0\)）。定义缩放矩阵：

\[ D^{(k)} = \text{diag}\left( \frac{1}{x_1^{(k)}}, \frac{1}{x_2^{(k)}}, \frac{1}{s^{(k)}} \right) \]

其中 \(s^{(k)} = 6 - x_1^{(k)} - 2x_2^{(k)} > 0\) 是松弛变量。但注意我们只有三个非负变量：\(x_1, x_2, s\)。在决策空间 \(\mathbb{R}^2\) 中，我们只需处理 \(x\) 的缩放。

实际上，更常用的形式是构造一个二次近似子问题，其搜索方向 \(d\) 由以下线性系统给出（通过一阶最优性条件线性化得到）：

\[ \begin{bmatrix} \nabla^2 f(x^{(k)}) + \mu^{(k)} (X^{(k)})^{-2} & \nabla g(x^{(k)}) \\ \nabla g(x^{(k)})^T & -S^{(k)} / \mu^{(k)} \end{bmatrix} \begin{bmatrix} d \\ \Delta y \end{bmatrix} = - \begin{bmatrix} \nabla f(x^{(k)}) - \mu^{(k)} (X^{(k)})^{-1} e - y^{(k)} \nabla g(x^{(k)}) \\ g(x^{(k)}) + s^{(k)} \end{bmatrix} \]

其中：

• \(X^{(k)} = \text{diag}(x_1^{(k)}, x_2^{(k)})\)，
• \(e = (1,1)^T\)，
• \(S^{(k)} = s^{(k)}\) 是标量（对应一个不等式约束），
• \(y^{(k)}\) 是当前对偶变量，
• \(\Delta y\) 是对偶变量的增量。

但对于我们的简单问题，我们可以直接推导障碍函数梯度并求解牛顿方向。

4. 第一步迭代（k=0）
   • 已知：\(x^{(0)} = (1, 2)^T\)，\(y^{(0)} = 1\)，\(\mu^{(0)} = 1\)，\(\sigma = 0.5\)。
   • 计算相关量：

\[ f(x) = (x_1-3)^2 + (x_2-2)^2, \quad \nabla f(x) = \begin{pmatrix} 2(x_1-3) \\ 2(x_2-2) \end{pmatrix} \]

\[ g(x) = x_1 + 2x_2 - 6, \quad \nabla g(x) = \begin{pmatrix} 1 \\ 2 \end{pmatrix} \]

 在 $x^{(0)} = (1,2)^T$：

\[ f(x^{(0)}) = (1-3)^2 + (2-2)^2 = 4, \quad \nabla f(x^{(0)}) = \begin{pmatrix} 2(1-3) \\ 2(2-2) \end{pmatrix} = \begin{pmatrix} -4 \\ 0 \end{pmatrix} \]

\[ g(x^{(0)}) = 1 + 4 - 6 = -1, \quad s^{(0)} = -g(x^{(0)}) = 1 > 0 \]

• 障碍函数梯度：

\[ \nabla B(x; \mu) = \nabla f(x) - \mu \begin{pmatrix} 1/x_1 \\ 1/x_2 \end{pmatrix} + \mu \frac{1}{6 - x_1 - 2x_2} \begin{pmatrix} 1 \\ 2 \end{pmatrix} \]

 注意最后一项符号：$\frac{\partial}{\partial x} [-\mu \ln(6-x_1-2x_2)] = \mu \frac{1}{6-x_1-2x_2} \nabla g(x)$，因为 $g(x) \leq 0$，$6-x_1-2x_2 = -g(x) = s > 0$。
 代入数值：

\[ \nabla B(x^{(0)}; \mu^{(0)}) = \begin{pmatrix} -4 \\ 0 \end{pmatrix} - 1 \cdot \begin{pmatrix} 1/1 \\ 1/2 \end{pmatrix} + 1 \cdot \frac{1}{1} \begin{pmatrix} 1 \\ 2 \end{pmatrix} = \begin{pmatrix} -4 \\ 0 \end{pmatrix} - \begin{pmatrix} 1 \\ 0.5 \end{pmatrix} + \begin{pmatrix} 1 \\ 2 \end{pmatrix} = \begin{pmatrix} -4 \\ 1.5 \end{pmatrix} \]

• 障碍函数海森矩阵：

\[ \nabla^2 f(x) = \begin{pmatrix} 2 & 0 \\ 0 & 2 \end{pmatrix} \]

\[ H_B = \nabla^2 f(x) + \mu \begin{pmatrix} 1/x_1^2 & 0 \\ 0 & 1/x_2^2 \end{pmatrix} + \mu \frac{1}{(6-x_1-2x_2)^2} \begin{pmatrix} 1 & 2 \\ 2 & 4 \end{pmatrix} \]

 代入数值：

\[ H_B^{(0)} = \begin{pmatrix} 2 & 0 \\ 0 & 2 \end{pmatrix} + 1 \cdot \begin{pmatrix} 1/1^2 & 0 \\ 0 & 1/2^2 \end{pmatrix} + 1 \cdot \frac{1}{1^2} \begin{pmatrix} 1 & 2 \\ 2 & 4 \end{pmatrix} = \begin{pmatrix} 2 & 0 \\ 0 & 2 \end{pmatrix} + \begin{pmatrix} 1 & 0 \\ 0 & 0.25 \end{pmatrix} + \begin{pmatrix} 1 & 2 \\ 2 & 4 \end{pmatrix} = \begin{pmatrix} 4 & 2 \\ 2 & 6.25 \end{pmatrix} \]

• 求解牛顿方向 \(d^{(0)}\)：

\[ H_B^{(0)} d^{(0)} = -\nabla B(x^{(0)}; \mu^{(0)}) \]

\[ \begin{pmatrix} 4 & 2 \\ 2 & 6.25 \end{pmatrix} \begin{pmatrix} d_1 \\ d_2 \end{pmatrix} = - \begin{pmatrix} -4 \\ 1.5 \end{pmatrix} = \begin{pmatrix} 4 \\ -1.5 \end{pmatrix} \]

 解线性方程组：
 从第一式：$4 d_1 + 2 d_2 = 4$ → $d_2 = 2 - 2 d_1$
 代入第二式：$2 d_1 + 6.25 (2 - 2 d_1) = -1.5$
 $2 d_1 + 12.5 - 12.5 d_1 = -1.5$
 $-10.5 d_1 = -14$ → $d_1 = 4/3 \approx 1.3333$
 于是 $d_2 = 2 - 2 \times 4/3 = 2 - 8/3 = -2/3 \approx -0.6667$

• 线搜索：为保证迭代点严格内点（\(x>0, s>0\)），需选择步长 \(\alpha^{(0)}\) 使 \(x^{(0)} + \alpha d^{(0)} > 0\) 且 \(s = 6 - x_1 - 2x_2 > 0\)。
  计算最大步长：

\[ \alpha_{\max} = \min\left( \min_{i: d_i < 0} \frac{-x_i^{(0)}}{d_i}, \min_{d_s < 0} \frac{-s^{(0)}}{d_s} \right) \]

 其中 $d_s = -\nabla g(x^{(0)})^T d = -(1, 2) \cdot (d_1, d_2) = - (d_1 + 2 d_2) = - (4/3 - 4/3) = 0$
 因为只有 $d_2 = -2/3 < 0$，所以：

\[ \alpha_{\max} = \frac{-x_2^{(0)}}{d_2} = \frac{-2}{-2/3} = 3 \]

 由于 $d_s = 0$，s 的变化率为零，不影响可行性。通常取步长 $\alpha = \min(1, 0.99 \alpha_{\max}) = \min(1, 2.97) = 1$。

• 更新迭代点：

\[ x^{(1)} = x^{(0)} + \alpha d^{(0)} = \begin{pmatrix} 1 \\ 2 \end{pmatrix} + 1 \cdot \begin{pmatrix} 4/3 \\ -2/3 \end{pmatrix} = \begin{pmatrix} 7/3 \\ 4/3 \end{pmatrix} \approx (2.3333, 1.3333) \]

• 更新对偶变量：在序列仿射缩放内点法中，对偶变量通常由一阶最优性条件近似更新。对偶变量对应不等式约束 \(g(x) \leq 0\) 的拉格朗日乘子，其近似值为：

\[ y^{(1)} = \frac{\mu^{(0)}}{s^{(0)}} = \frac{1}{1} = 1 \]

 更精确的更新可用牛顿步中解出的对偶增量，但这里简单采用此近似。

• 更新障碍参数：

\[ \mu^{(1)} = \sigma \mu^{(0)} = 0.5 \times 1 = 0.5 \]

• 计算近似目标函数值（原目标函数）：

\[ f(x^{(1)}) = (7/3 - 3)^2 + (4/3 - 2)^2 = (-2/3)^2 + (-2/3)^2 = 4/9 + 4/9 = 8/9 \approx 0.8889 \]

5. 第二步迭代（k=1）
   • 当前：\(x^{(1)} = (7/3, 4/3)^T\)，\(y^{(1)} = 1\)，\(\mu^{(1)} = 0.5\)。
   • 计算相关量：

\[ g(x^{(1)}) = 7/3 + 2 \times 4/3 - 6 = 7/3 + 8/3 - 6 = 15/3 - 6 = 5 - 6 = -1, \quad s^{(1)} = 1 \]

\[ \nabla f(x^{(1)}) = \begin{pmatrix} 2(7/3-3) \\ 2(4/3-2) \end{pmatrix} = \begin{pmatrix} 2(-2/3) \\ 2(-2/3) \end{pmatrix} = \begin{pmatrix} -4/3 \\ -4/3 \end{pmatrix} \]

• 障碍函数梯度：

\[ \nabla B(x^{(1)}; \mu^{(1)}) = \begin{pmatrix} -4/3 \\ -4/3 \end{pmatrix} - 0.5 \begin{pmatrix} 3/7 \\ 3/4 \end{pmatrix} + 0.5 \cdot \frac{1}{1} \begin{pmatrix} 1 \\ 2 \end{pmatrix} = \begin{pmatrix} -4/3 \\ -4/3 \end{pmatrix} - \begin{pmatrix} 1.5/7 \\ 1.5/4 \end{pmatrix} + \begin{pmatrix} 0.5 \\ 1 \end{pmatrix} \]

 计算数值：
 第一分量：$-1.3333 - 0.2143 + 0.5 = -1.0476$
 第二分量：$-1.3333 - 0.375 + 1 = -0.7083$
 即 $\nabla B \approx (-1.0476, -0.7083)^T$

• 障碍函数海森矩阵：

\[ H_B^{(1)} = \begin{pmatrix} 2 & 0 \\ 0 & 2 \end{pmatrix} + 0.5 \begin{pmatrix} (3/7)^2 & 0 \\ 0 & (3/4)^2 \end{pmatrix} + 0.5 \cdot \frac{1}{1^2} \begin{pmatrix} 1 & 2 \\ 2 & 4 \end{pmatrix} = \begin{pmatrix} 2 & 0 \\ 0 & 2 \end{pmatrix} + 0.5 \begin{pmatrix} 9/49 & 0 \\ 0 & 9/16 \end{pmatrix} + 0.5 \begin{pmatrix} 1 & 2 \\ 2 & 4 \end{pmatrix} \]

 计算：
 第一项+第三项：$\begin{pmatrix} 2 & 0 \\ 0 & 2 \end{pmatrix} + \begin{pmatrix} 0.5 & 1 \\ 1 & 2 \end{pmatrix} = \begin{pmatrix} 2.5 & 1 \\ 1 & 4 \end{pmatrix}$
 加第二项：$\begin{pmatrix} 2.5 & 1 \\ 1 & 4 \end{pmatrix} + \begin{pmatrix} 0.5 \times 9/49 \approx 0.0918 & 0 \\ 0 & 0.5 \times 9/16 = 0.28125 \end{pmatrix} = \begin{pmatrix} 2.5918 & 1 \\ 1 & 4.2813 \end{pmatrix}$

• 解牛顿方向：

\[ \begin{pmatrix} 2.5918 & 1 \\ 1 & 4.2813 \end{pmatrix} \begin{pmatrix} d_1 \\ d_2 \end{pmatrix} = - \begin{pmatrix} -1.0476 \\ -0.7083 \end{pmatrix} = \begin{pmatrix} 1.0476 \\ 0.7083 \end{pmatrix} \]

 由第一式：$2.5918 d_1 + d_2 = 1.0476$ → $d_2 = 1.0476 - 2.5918 d_1$
 代入第二式：$d_1 + 4.2813 (1.0476 - 2.5918 d_1) = 0.7083$
 $d_1 + 4.2813 \times 1.0476 - 4.2813 \times 2.5918 d_1 = 0.7083$
 计算：$4.2813 \times 1.0476 \approx 4.485$，$4.2813 \times 2.5918 \approx 11.097$
 方程：$d_1 + 4.485 - 11.097 d_1 = 0.7083$
 $-10.097 d_1 = 0.7083 - 4.485 = -3.7767$
 $d_1 \approx 0.3741$，则 $d_2 = 1.0476 - 2.5918 \times 0.3741 \approx 1.0476 - 0.9695 = 0.0781$

• 线搜索：计算 \(d_s = -(d_1 + 2 d_2) = -(0.3741 + 0.1562) = -0.5303\)
  由于 \(d_1, d_2 > 0\)，对 \(x>0\) 无限制；但 \(d_s < 0\)，需满足 \(s + \alpha d_s > 0\)：
  \(\alpha_{\max} = -s^{(1)} / d_s = -1 / (-0.5303) \approx 1.886\)
  取步长 \(\alpha = \min(1, 0.99 \times 1.886) = 1\)
• 更新迭代点：

\[ x^{(2)} = x^{(1)} + \alpha d^{(1)} = \begin{pmatrix} 7/3 \\ 4/3 \end{pmatrix} + \begin{pmatrix} 0.3741 \\ 0.0781 \end{pmatrix} \approx (2.6667, 1.4114) \]

• 更新对偶变量：

\[ y^{(2)} = \frac{\mu^{(1)}}{s^{(1)}} = \frac{0.5}{1} = 0.5 \]

• 更新障碍参数：

\[ \mu^{(2)} = \sigma \mu^{(1)} = 0.5 \times 0.5 = 0.25 \]

• 近似目标函数值：

\[ f(x^{(2)}) = (2.6667-3)^2 + (1.4114-2)^2 = (-0.3333)^2 + (-0.5886)^2 \approx 0.1111 + 0.3465 = 0.4576 \]

6. 结果与解释
   经过两步迭代，点从 \((1,2)\) 移动到约 \((2.6667, 1.4114)\)，目标函数值从 4 降到约 0.4576，逐渐接近理论最优点（此问题理论最优解在 \(x^* = (3, 1.5)\)，满足约束 \(x_1+2x_2=6\)，目标值为 0）。障碍参数从 1 降到 0.25，对偶变量从 1 降到 0.5，反映了内点法逐步逼近最优解的过程。