非线性规划中的序列内点反射信赖域混合算法基础题
题目描述
考虑一个带不等式约束的非线性规划问题：
min f(x)
s.t. c_i(x) ≤ 0, i = 1, 2, ..., m
x ∈ R^n
其中，目标函数 f(x) 和约束函数 c_i(x) 都是二阶连续可微的，并且可能非凸。序列内点反射信赖域混合算法是一种结合了内点法（维持解的可行性）、反射技术（处理边界附近的震荡）和信赖域法（控制步长和方向质量）的迭代求解框架。本基础题旨在理解该混合算法的核心迭代步骤，特别是如何构造子问题并确保收敛。
给定一个具体例子：
min f(x) = (x_1 - 3)^2 + (x_2 - 2)^2
s.t. c_1(x) = x_1^2 + x_2^2 - 4 ≤ 0
c_2(x) = -x_1 ≤ 0
c_3(x) = -x_2 ≤ 0
初始点 x⁰ = (0.5, 0.5)，该点严格可行（即所有约束 c_i(x⁰) < 0）。内点参数 μ⁰ = 10。信赖域初始半径 Δ⁰ = 1.0。请你详细解释序列内点反射信赖域混合算法从初始点 x⁰ 出发，计算第一个迭代点 x¹ 的完整过程。
解题过程
我们的核心思路是：在当前迭代点 x^k，我们构造一个结合了目标函数和约束的“代理模型”（通常是二次模型），并在一个信赖域内求解这个子问题。子问题的解（称为试验步）可能位于边界或导致违反约束，因此引入反射操作将其“拉回”可行域。最终，我们根据实际目标函数下降量与模型预测下降量的比值，来决定是否接受该步长，并更新信赖域半径和内点参数。
步骤 1: 构造障碍函数与二次近似模型
-
内点法思想引入障碍函数：内点法通过在目标函数中添加一个障碍项来阻止迭代点穿越约束边界。在迭代 k，我们使用对数障碍函数：
B(x; μ^k) = f(x) - μ^k * Σ_{i=1}^{m} log(-c_i(x))
其中 μ^k > 0 是障碍参数。当 x 趋近于边界时，-c_i(x) → 0⁺，log(-c_i(x)) → -∞，障碍项 μ^k * Σ log(-c_i(x)) → +∞，从而在函数值上惩罚接近边界的点，使其保持在可行域内部。 -
计算障碍函数在 x^k 处的梯度与 Hessian：我们需要在当前点 x^k 构造障碍函数 B(x; μ^k) 的二次近似（泰勒展开）。为此，需计算梯度 ∇B(x^k; μ^k) 和 Hessian ∇²B(x^k; μ^k)。计算过程如下：
- 目标函数部分：∇f(x) = [2(x_1-3)， 2(x_2-2)]^T。在 x⁰=(0.5,0.5) 处，∇f(x⁰) = [2*(0.5-3)， 2*(0.5-2)]^T = [-5, -3]^T。
∇²f(x) 是常数矩阵：[[2, 0], [0, 2]]。 - 约束函数 c_1(x) = x_1^2 + x_2^2 - 4。梯度 ∇c_1(x) = [2x_1, 2x_2]^T，在 x⁰ 处为 [1, 1]^T。Hessian ∇²c_1(x) = [[2, 0], [0, 2]]。
约束函数 c_2(x) = -x_1。梯度 ∇c_2(x) = [-1, 0]^T， Hessian 为零矩阵。
约束函数 c_3(x) = -x_2。梯度 ∇c_3(x) = [0, -1]^T， Hessian 为零矩阵。 - 在 x⁰ 处，约束函数值：c_1(x⁰) = 0.5^2+0.5^2-4 = -3.5， c_2(x⁰) = -0.5， c_3(x⁰) = -0.5。
- 障碍项梯度：∇[-μ^k Σ log(-c_i(x))] = μ^k Σ [∇c_i(x) / (-c_i(x))]。
对于 i=1: ∇c_1 / (-c_1) = [1, 1]^T / 3.5 ≈ [0.2857, 0.2857]^T
i=2: ∇c_2 / (-c_2) = [-1, 0]^T / 0.5 = [-2, 0]^T
i=3: ∇c_3 / (-c_3) = [0, -1]^T / 0.5 = [0, -2]^T
求和：Σ ∇c_i / (-c_i) ≈ [0.2857-2+0, 0.2857+0-2]^T ≈ [-1.7143, -1.7143]^T
乘以 μ⁰=10：障碍项梯度 ≈ 10 * [-1.7143, -1.7143]^T = [-17.143, -17.143]^T - 障碍函数梯度 ∇B(x⁰; μ⁰) = ∇f(x⁰) + 障碍项梯度 = [-5, -3]^T + [-17.143, -17.143]^T = [-22.143, -20.143]^T。
- 障碍项 Hessian：∇²[-μ^k Σ log(-c_i(x))] = μ^k Σ [ (∇c_i ∇c_i^T) / (c_i)^2 + ∇²c_i / (-c_i) ]。我们通常忽略二阶项 ∇²c_i / (-c_i) 或使用拟牛顿法（如BFGS）近似整个 Hessian，以简化计算并保持正定性。在基础题中，我们常使用拟牛顿近似（例如初始化为单位矩阵）或只计算一阶信息并忽略二阶项，构造一个简单的正定模型。 这里为了清晰，我们采用一个简化策略：使用目标函数的精确 Hessian 加上一个由梯度信息构造的对角占正定矩阵来近似障碍项的 Hessian。更常用的做法是使用拟牛顿法（BFGS）更新一个矩阵 H^k 来近似 ∇²B。我们假设在第一步，我们初始化 H⁰ = I（单位矩阵）。所以，我们的二次模型为：
m^k(s) = B(x^k; μ^k) + (∇B(x^k; μ^k))^T s + (1/2) s^T H^k s
其中 s = x - x^k 是试探步，H⁰ = I。
- 目标函数部分：∇f(x) = [2(x_1-3)， 2(x_2-2)]^T。在 x⁰=(0.5,0.5) 处，∇f(x⁰) = [2*(0.5-3)， 2*(0.5-2)]^T = [-5, -3]^T。
步骤 2: 求解信赖域子问题
- 构造信赖域子问题：在当前点 x^k，我们希望找到一个步长 s，使得模型 m^k(s) 在信赖域 ‖s‖ ≤ Δ^k 内尽可能减小。子问题为：
min_s m^k(s) = B(x^k; μ^k) + (∇B(x^k; μ^k))^T s + (1/2) s^T H^k s
s.t. ‖s‖ ≤ Δ^k
这里 ‖·‖ 通常指2-范数。由于 H⁰ = I 是正定矩阵，该子问题是严格凸的，其全局解 s* 可以通过求解 (H^k + λI) s* = -∇B(x^k; μ^k) 对某个 λ ≥ 0 来获得（即信赖域问题的精确或近似解）。在基础题中，我们可以采用折线法（Cauchy点与牛顿步的折中）或直接求解这个线性系统并调整 λ 以满足约束。为了简化，我们这里计算最速下降方向（负梯度方向）的截断步作为近似解：
s_c = - (Δ^k / ‖∇B‖) * (∇B / ‖∇B‖) = - (Δ^k / ‖∇B‖) ∇B
其中 ∇B = ∇B(x⁰; μ⁰) = [-22.143, -20.143]^T， ‖∇B‖ ≈ sqrt((-22.143)^2 + (-20.143)^2) ≈ sqrt(490.3 + 405.7) ≈ sqrt(896) ≈ 29.93。
Δ⁰ = 1.0，所以 s_c = - (1.0 / 29.93) * [-22.143, -20.143]^T ≈ [0.740, 0.673]^T。
这个步长 s_c 满足 ‖s_c‖ = 1.0。于是，试验点 x_trial = x⁰ + s_c = (0.5+0.740, 0.5+0.673) = (1.240, 1.173)。
步骤 3: 反射操作（Reflection）
- 检查试验点可行性：计算在 x_trial 处的约束函数值：c_1(x_trial) = 1.240^2 + 1.173^2 - 4 ≈ 1.5376 + 1.3759 - 4 = -1.0865 (可行，<0)。c_2 = -1.240 < 0, c_3 = -1.173 < 0。所以 x_trial 是严格可行的。如果试验点违反了某些约束（c_i(x_trial) > 0），则需要将其“反射”回可行域。反射操作通常沿着约束边界的法向方向将点推回可行域内部。一个简单的策略是：如果 c_i(x_trial) > 0，则沿着负梯度方向 -∇c_i(x_trial) 回退，直到 c_i ≤ 0。更系统的方法可能涉及求解一个小的线性或二次规划来找到最近的可行点。由于本例中 x_trial 可行，反射步骤可跳过。但算法框架中通常预设此步。
步骤 4: 接受性检验与参数更新
-
计算实际下降与预测下降：
- 实际下降：Ared^k = B(x^k; μ^k) - B(x_trial; μ^k)
首先计算 B(x⁰; μ⁰) = f(x⁰) - μ⁰ Σ log(-c_i(x⁰))
f(x⁰) = (0.5-3)^2 + (0.5-2)^2 = 6.25 + 2.25 = 8.5
-Σ log(-c_i) = -[log(3.5) + log(0.5) + log(0.5)] = -[1.2528 + (-0.6931) + (-0.6931)] = -(-0.1334) = 0.1334? 我们仔细计算：
log(-c_1) = log(3.5) ≈ 1.2528
log(-c_2) = log(0.5) ≈ -0.6931
log(-c_3) = log(0.5) ≈ -0.6931
和 = 1.2528 -0.6931 -0.6931 = -0.1334
所以 -μ⁰ Σ log(-c_i) = -10 * (-0.1334) = 1.334
B(x⁰; μ⁰) = 8.5 + 1.334 = 9.834
然后计算 B(x_trial; μ⁰)：
f(x_trial) = (1.240-3)^2 + (1.173-2)^2 = (-1.76)^2 + (-0.827)^2 = 3.0976 + 0.6839 = 3.7815
c_1 = -1.0865, c_2 = -1.240, c_3 = -1.173
log(-c_1) = log(1.0865) ≈ 0.0831, log(-c_2) = log(1.240) ≈ 0.2151, log(-c_3) = log(1.173) ≈ 0.1595
和 = 0.0831+0.2151+0.1595 = 0.4577
-μ⁰ Σ log(-c_i) = -10 * 0.4577 = -4.577
B(x_trial; μ⁰) = 3.7815 - 4.577 = -0.7955
实际下降 Ared^k = 9.834 - (-0.7955) = 10.6295 - 预测下降：Pred^k = m^k(0) - m^k(s_c) = [B(x⁰) + 0 + 0] - [B(x⁰) + ∇B^T s_c + 1/2 s_c^T H⁰ s_c] = -[∇B^T s_c + 1/2 s_c^T s_c]
∇B^T s_c = [-22.143, -20.143] · [0.740, 0.673] = (-22.1430.740) + (-20.1430.673) = -16.386 -13.556 = -29.942
s_c^T s_c = 0.740^2+0.673^2 = 0.5476+0.4529 = 1.0005 ≈ 1.0
所以 Pred^k = -[-29.942 + 0.5*1.0] = -[-29.942 + 0.5] = -[-29.442] = 29.442 - 比值 ρ^k = Ared^k / Pred^k = 10.6295 / 29.442 ≈ 0.361
- 实际下降：Ared^k = B(x^k; μ^k) - B(x_trial; μ^k)
-
决定是否接受试验点并更新信赖域半径：
- 通常设定阈值 η₁, η₂，例如 η₁=0.01, η₂=0.9。
- 如果 ρ^k ≥ η₁，则认为模型是足够精确的，接受该步：x^{k+1} = x_trial。本例中 ρ^k ≈ 0.361 > 0.01，所以接受，即 x¹ = (1.240, 1.173)。
- 更新信赖域半径 Δ^{k+1}：
如果 ρ^k ≥ η₂（本例 0.361 < 0.9），则保持半径不变或略增，例如 Δ^{k+1} = Δ^k = 1.0。
如果 ρ^k < η₁，则拒绝步长，缩小信赖域，例如 Δ^{k+1} = 0.5 * Δ^k。
这里我们保持 Δ¹ = 1.0。
-
更新障碍参数 μ：通常，如果步长被接受，并且迭代点趋近于边界，我们可以减小 μ 以降低障碍项的权重，使解更接近原问题的最优解。一个简单的更新规则是：μ^{k+1} = σ * μ^k，其中 σ ∈ (0,1)，例如 σ=0.1。但通常在算法初期不会每步都减少。在基础题中，我们可以在完成一定次数的成功迭代后更新。为简化，假设我们保持 μ¹ = μ⁰ = 10 进行下一步。
步骤 5: 迭代与终止
- 以上完成了从 x⁰ 到 x¹ 的一次完整迭代。算法会重复步骤1-4，直到满足某个收敛准则，例如 ‖∇B(x^k; μ^k)‖ 足够小，或者 μ^k 足够小且步长 ‖s^k‖ 足够小。
总结：序列内点反射信赖域混合算法的主要步骤包括：1) 构造包含障碍函数的二次近似模型；2) 在信赖域内求解该模型得到试验步；3) 若试验点不可行，则进行反射操作将其映射回可行域；4) 计算实际下降与预测下降的比值，据此决定是否接受试验点并更新信赖域半径；5) 适时减小障碍参数 μ。该框架平衡了可行性保持（内点法）、全局收敛性（信赖域法）和边界处理（反射技术）。