非线性规划中的改进内点反射梯度法基础题

题目描述
考虑如下带不等式约束的非线性规划问题：

\[\begin{aligned} \min_{x \in \mathbb{R}^n} \quad & f(x) \\ \text{s.t.} \quad & c_i(x) \le 0, \quad i = 1, 2, \dots, m \end{aligned} \]

其中，目标函数 \(f: \mathbb{R}^n \to \mathbb{R}\) 连续可微，约束函数 \(c_i: \mathbb{R}^n \to \mathbb{R}\) 也连续可微，且可行域记为 \(\Omega = \{x \mid c_i(x) \le 0, i=1,\dots,m\}\)。
内点反射梯度法是一种结合内点法与梯度投影思想的算法，其基本思路是从一个可行内点出发，沿负梯度方向移动，若步进后越出可行域边界，则通过“反射”操作将迭代点拉回可行域内部。本题目中的“改进”主要体现在对反射策略的调整和步长的自适应选择上。
现在给定一个具体算例：

\[\begin{aligned} \min \quad & f(x) = (x_1 - 3)^2 + (x_2 - 2)^2 \\ \text{s.t.} \quad & c_1(x) = x_1 + x_2 - 4 \le 0 \\ & c_2(x) = -x_1 \le 0 \\ & c_3(x) = -x_2 \le 0 \end{aligned} \]

该问题的最优解为 \(x^* = (2, 2)^T\)，最优值 \(f(x^*) = 1\)。
请逐步阐述如何使用改进的内点反射梯度法求解此问题，包括算法的完整步骤、关键参数的设置、以及如何确保在迭代过程中保持可行性并收敛到最优解。

解题过程循序渐进讲解

第一步：理解问题与算法框架

1. 问题分析

   • 目标函数 \(f(x)\) 是凸二次函数，约束是线性不等式，因此该问题是凸优化问题，局部最优即全局最优。
   • 可行域是一个三角形区域：\(x_1 \ge 0, x_2 \ge 0, x_1 + x_2 \le 4\)。最优解位于约束 \(c_1(x)=0\) 的边界上。

2. 算法核心思想

   • 内点法思想：迭代点始终保持在严格可行域内部（即 \(c_i(x) < 0\)），避免过早触及边界导致数值困难。
   • 反射操作：当沿负梯度方向移动时，若步进后点 \(x_{\text{new}}\) 越出边界（即某个 \(c_i(x_{\text{new}}) \ge 0\)），不直接截断到边界，而是进行“反射”，即模拟光线在边界上的反射，将点弹回可行域内部，同时尽量保持下降方向的有效性。
   • 改进点：
     (1) 反射规则更精细，根据违反约束的程度和梯度方向调整反射向量；
     (2) 步长采用自适应机制，结合 Armijo 条件保证充分下降，并利用可行域内点的距离信息动态调整。

第二步：算法详细步骤

1. 初始化

   • 选取严格内点 \(x^0\)。对于本例，可取 \(x^0 = (1, 1)^T\)，满足所有约束严格不等式（\(c_1= -2<0, c_2=-1<0, c_3=-1<0\)）。
   • 设置参数：初始步长 \(\alpha_0 = 1\)，收缩因子 \(\beta = 0.5\)，Armijo 条件参数 \(\sigma = 0.01\)，容许误差 \(\epsilon = 10^{-6}\)。
   • 迭代计数器 \(k := 0\)。

2. 计算梯度
   计算当前点 \(x^k\) 的目标函数梯度：

\[ \nabla f(x^k) = \begin{bmatrix} 2(x_1^k - 3) \\ 2(x_2^k - 2) \end{bmatrix} \]

在本例中，梯度方向指向无约束最优点 (3,2)，但受约束限制，实际最优解是 (2,2)。

3. 确定搜索方向

   • 基本搜索方向为负梯度方向：\(d^k = -\nabla f(x^k)\)。
   • 但需注意：在靠近边界时，负梯度方向可能指向可行域外。本算法采用“反射梯度方向”作为实际搜索方向，具体在下一步中确定。

4. 反射操作与可行性保持

   • 尝试步进：计算试探点 \(x_{\text{try}} = x^k + \alpha_k d^k\)。
   • 检查可行性：计算所有约束函数值 \(c_i(x_{\text{try}})\)。
     • 若所有 \(c_i(x_{\text{try}}) < 0\)，则 \(x_{\text{try}}\) 是可行内点，直接进入步长搜索。
     • 若存在某个约束 \(c_j(x_{\text{try}}) \ge 0\)，则说明沿 \(d^k\) 方向移动会越出第 \(j\) 个约束边界。此时进行反射：
       (1) 计算边界交点：通过线性插值或求解 \(c_j(x^k + t d^k)=0\) 得到步长 \(t_b \in (0, \alpha_k)\)，使得 \(x_b = x^k + t_b d^k\) 正好在边界上。
       (2) 反射方向计算：在边界点 \(x_b\) 处，计算约束梯度 \(\nabla c_j(x_b)\)。反射方向 \(d_r\) 由入射方向 \(d^k\) 关于边界切平面的反射得到，公式为：

\[ d_r = d^k - 2 \frac{d^k \cdot \nabla c_j(x_b)}{\|\nabla c_j(x_b)\|^2} \nabla c_j(x_b) \]

   注意：这里使用了向量反射公式，反射后方向与边界“镜像”对称，有助于探索可行域内部。  
   (3) 新的试探点：从边界点沿反射方向移动剩余步长：$x_{\text{try}} = x_b + (\alpha_k - t_b) d_r$。  
   (4) 若反射后仍不可行（即又触及其他边界），可重复反射或缩减步长。本算法中，若反射后仍不可行，则直接进入步长缩减。

5. 自适应步长搜索
   • 用上述步骤（可能包含反射）得到试探点 \(x_{\text{try}}\) 后，检查 Armijo 条件：

\[ f(x_{\text{try}}) \le f(x^k) + \sigma \nabla f(x^k)^T (x_{\text{try}} - x^k) \]

 若条件满足，则接受该点，令 $x^{k+1} = x_{\text{try}}$，并适当增大步长（如 $\alpha_{k+1} = 1.2 \alpha_k$），以加速收敛。  

• 若条件不满足，则按 \(\alpha_k := \beta \alpha_k\) 缩小步长，重新从第4步开始尝试（即重新计算试探点，可能需要新的反射）。
• 步长自适应机制：初始步长设为1，若连续几次迭代都接受，则适当增大步长；若频繁发生反射或 Armijo 条件不满足，则减小步长。

6. 收敛判断
   • 计算当前点的 KKT 条件误差：对于内点法，可检查梯度在可行域内的投影。简单起见，可用梯度范数结合约束违反量：

\[ \eta_k = \|\nabla f(x^k) + \sum_{i=1}^m \lambda_i^k \nabla c_i(x^k)\| + \sum_{i=1}^m |\min(-c_i(x^k), \lambda_i^k)| \]

 其中 $\lambda_i^k$ 是约束的拉格朗日乘子估计，可用边界附近的约束违反程度近似，例如 $\lambda_i^k = \max(0, -\mu / c_i(x^k))$，$\mu>0$ 是小参数。  

• 若 \(\eta_k < \epsilon\) 或迭代次数超过上限，则停止，输出当前点作为近似最优解。否则，令 \(k:=k+1\)，返回第2步。

第三步：应用到给定算例

1. 迭代示例（第1轮）
   • 初始点 \(x^0=(1,1)\), \(f=5\), 梯度 \(\nabla f = (-4, -2)^T\)，负梯度方向 \(d^0 = (4, 2)^T\)。
   • 尝试步长 \(\alpha_0=1\): \(x_{\text{try}} = (5, 3)\)，违反约束 \(c_1\) 和 \(c_2, c_3\)。
   • 检测到最先碰到的边界是 \(c_1=0\)（因为沿 \(d^0\) 方向运动时，\(c_1\) 从负变零最快）。计算交点：
     解 \(c_1(x^0 + t d^0) = (1+4t) + (1+2t) - 4 = 6t - 2 = 0\)，得 \(t_b = 1/3\)。
     边界点 \(x_b = (1+4/3, 1+2/3) = (7/3, 5/3)\)。
   • 计算反射：约束梯度 \(\nabla c_1 = (1,1)^T\)，归一化因子为 \(\|\nabla c_1\|^2=2\)。
     入射向量在边界点处仍为 \(d^0=(4,2)^T\)，计算反射：

\[ d_r = \begin{pmatrix}4\\2\end{pmatrix} - 2 \frac{4+2}{2} \begin{pmatrix}1\\1\end{pmatrix} = \begin{pmatrix}4\\2\end{pmatrix} - 6 \begin{pmatrix}1\\1\end{pmatrix} = \begin{pmatrix}-2\\-4\end{pmatrix} \]

• 剩余步长为 \(1 - 1/3 = 2/3\)，反射点 \(x_{\text{try}} = (7/3, 5/3) + (2/3)(-2, -4) = (1, -1)\)。
  但此点违反 \(c_3\)（即 \(x_2 \ge 0\)），说明反射后碰到另一边界。此时算法可进行二次反射，但为简化，我们选择缩减步长。

2. 步长缩减
   • 令 \(\alpha_0 := 0.5\)，重新尝试。此时 \(t_b\) 会变化，但过程类似。通过几次调整，最终可找到一个步长，使得反射后点仍在可行域内且满足 Armijo 条件。实际上，对于此凸问题，适当小的步长可保证可行性。
   • 经过若干次迭代，点列将收敛到最优解 (2,2)。在最优解处，梯度 \(\nabla f = (-2, 0)^T\)，约束 \(c_1\) 为积极约束，梯度可表示为 \(\nabla f = -\lambda_1 \nabla c_1\)，其中 \(\lambda_1=2\)，满足 KKT 条件。

第四步：收敛性与改进说明

• 收敛性：由于目标函数凸、约束线性，且在反射操作中保证了迭代点始终是可行内点，算法产生的序列会收敛到全局最优解。改进的反射策略减少了在边界处的振荡，自适应步长提高了收敛速度。
• 与经典梯度投影法的区别：梯度投影法通常将点投影到可行域，计算投影较费时；而反射法通过几何反射快速将点拉回内部，适合简单边界形状（如线性约束）。
• 关键改进点总结：
  1. 反射方向计算使用了精确的向量反射公式，使得反射后方向更合理地指向可行域内部。
  2. 步长选择结合了可行性检测与 Armijo 条件，动态调整以平衡效率与稳定性。
  3. 可处理多个约束被违反的情形，通过顺序反射或步长缩减确保可行性。

第五步：算法伪代码

输入：初始内点 x0，参数 α0=1，β=0.5，σ=0.01，ε=1e-6，最大迭代次数 K。
设置 k=0，xk=x0，α=α0。
while k < K 且 收敛误差 > ε：
   计算梯度 gk = ∇f(xk)。
   设搜索方向 dk = -gk。
   设置步长 αc = α。
   while True：
      尝试点 xtry = xk + αc dk。
      检查可行性：若存在 i 使 ci(xtry) ≥ 0，则找到第一个这样的约束 j，计算交点步长 tb，边界点 xb，反射方向 dr，并令 xtry = xb + (αc - tb) dr。
      若 xtry 仍不可行，则令 αc = β αc 并继续循环。
      若 f(xtry) ≤ f(xk) + σ gk^T (xtry - xk)：
         接受点：x_{k+1} = xtry。
         适当增加步长：α = min(1.2α, 1)。
         break
      否则：
         缩减步长：αc = β αc。
   结束内层循环。
   更新收敛误差（如梯度投影范数）。
   k = k+1。
输出：xk 作为近似最优解。

通过以上步骤，改进的内点反射梯度法能够系统地求解给定的约束优化问题，并在保持可行性的前提下逐步逼近最优解。