非线性规划中的内点反射梯度-信赖域混合算法基础题
题目描述：
考虑如下带有非线性不等式约束的优化问题：
最小化 \(f(x) = (x_1 - 3)^2 + (x_2 - 4)^2\)
约束条件为：
\(c_1(x) = x_1^2 + x_2^2 - 9 \leq 0\) （可行域在半径为3的圆内）
\(c_2(x) = -x_1 - x_2 + 1 \leq 0\) （可行域在直线右上方）
要求使用内点反射梯度-信赖域混合算法求解。该算法结合了内点法的障碍函数思想、反射梯度法的可行性保持机制，以及信赖域法的局部模型管理策略。目标是找到满足约束的局部最优点，并理解混合策略如何平衡可行性、最优性和收敛速度。
解题过程循序渐进讲解：
1. 算法核心思想
本算法是三种方法的混合：
- 内点法：通过障碍函数（如对数障碍）将约束问题转化为一系列无约束子问题，确保迭代点始终在可行域内部。
- 反射梯度法：当迭代点靠近边界时，梯度方向可能指向可行域外，反射梯度通过“反射”操作调整搜索方向，使其沿边界切向移动，保持可行性。
- 信赖域法：在每一步构建目标函数的局部二次模型，并在一个信赖域内求解该模型，通过比较实际下降量与预测下降量动态调整信赖域半径，保证全局收敛性。
混合策略的优势：内点法处理约束，反射梯度改善边界附近的搜索方向，信赖域控制步长可靠性。
2. 问题具体化与初始化
给定目标函数 \(f(x) = (x_1 - 3)^2 + (x_2 - 4)^2\)，这是一个以 (3,4) 为圆心的二次函数。约束 \(c_1(x)\) 表示圆内区域，\(c_2(x)\) 表示半平面 \(x_1 + x_2 \geq 1\)。可行域是圆与半平面的交集。
初始点选择：必须在可行域内部，例如 \(x^{(0)} = (0.5, 0.5)\)，满足 \(c_1 = -8.5 < 0\)，\(c_2 = 0 > 0\)? 注意：检查 \(c_2(0.5,0.5) = -1+1=0\)，恰好满足等式约束，但内点法要求严格内点（约束值 < 0），因此调整为 \(x^{(0)} = (1.0, 1.0)\)，此时 \(c_1 = -7 < 0\)，\(c_2 = -1 < 0\)，是严格内点。
设置初始信赖域半径 \(\Delta_0 = 0.5\)，障碍参数 \(\mu = 1\)（后续会减小），收敛容差 \(\epsilon = 10^{-6}\)。
3. 构造障碍函数
内点法使用对数障碍函数将约束问题转化为：
\(B(x; \mu) = f(x) - \mu \sum_{i=1}^2 \ln(-c_i(x))\)
其中 \(\mu > 0\) 是障碍参数，\(-\ln(-c_i(x))\) 在 \(c_i(x) \to 0^-\) 时趋于无穷大，从而阻止迭代点越界。
在当前例子中：
\(B(x; \mu) = (x_1-3)^2 + (x_2-4)^2 - \mu [\ln(9 - x_1^2 - x_2^2) + \ln(x_1 + x_2 - 1)]\)
注意：\(c_1 \leq 0\) 即 \(9 - x_1^2 - x_2^2 \geq 0\)，所以 \(\ln(9 - x_1^2 - x_2^2)\) 定义合理；\(c_2 \leq 0\) 即 \(x_1 + x_2 - 1 \geq 0\)，所以 \(\ln(x_1 + x_2 - 1)\) 定义合理。
4. 计算反射梯度方向
标准梯度方向为 \(\nabla B(x; \mu)\)，但靠近边界时该方向可能指向可行域外。反射梯度法修正如下：
- 计算约束的梯度：\(\nabla c_1 = (2x_1, 2x_2)^T\)，\(\nabla c_2 = (-1, -1)^T\)。
- 识别“积极约束”：当 \(-c_i(x) < \delta\)（δ 是一个小正数，如 0.1）时，认为该约束接近活跃。假设当前点 \(x^{(k)}\) 满足 \(-c_1(x)\) 很小（即靠近圆边界），而 \(c_2\) 不活跃。
- 定义反射算子：对于接近活跃的约束，将其梯度单位化得到法向量 \(n_i = \nabla c_i / \|\nabla c_i\|\)。然后计算反射方向：
\(d_{\text{refl}} = d - 2 (d^T n_i) n_i\)，其中 \(d = -\nabla B\) 是负梯度方向。这相当于将 \(d\) 关于边界切平面做镜面反射。 - 最终搜索方向 \(d_k\) 取为反射后的方向（如果约束活跃）或原始负梯度方向。
5. 信赖域子问题构建
在点 \(x^{(k)}\) 处，用二次模型近似障碍函数：
\(m_k(s) = B(x^{(k)}; \mu) + \nabla B(x^{(k)}; \mu)^T s + \frac{1}{2} s^T H_k s\)
其中 \(H_k\) 是 \(B\) 的Hessian近似（可用BFGS更新或精确计算）。约束为 \(\|s\| \leq \Delta_k\)（信赖域半径）。
但这里方向使用了反射梯度 \(d_k\)，因此子问题可简化为：沿 \(d_k\) 方向在信赖域内求最优步长 \(\alpha\)：
最小化 \(m_k(\alpha d_k) = B(x^{(k)}) + \alpha \nabla B^T d_k + \frac{1}{2} \alpha^2 d_k^T H_k d_k\)
满足 \(|\alpha| \|d_k\| \leq \Delta_k\)。
这是一个一维问题，可直接求解：最优 \(\alpha^* = \text{argmin}_{\alpha} m_k(\alpha)\)，且 \(|\alpha| \leq \Delta_k / \|d_k\|\)。如果 \(d_k^T H_k d_k > 0\)，则 \(\alpha^* = -\frac{\nabla B^T d_k}{d_k^T H_k d_k}\)，并截断到信赖域边界。
6. 接受步长与更新信赖域半径
计算候选点 \(x^+ = x^{(k)} + \alpha^* d_k\)。检查可行性：必须满足 \(c_i(x^+) < 0\)。如果不可行，则缩小 \(\alpha^*\) 或减少 \(\Delta_k\)。
计算实际下降量：\(\text{ared} = B(x^{(k)}) - B(x^+)\)
预测下降量：\(\text{pred} = m_k(0) - m_k(\alpha^* d_k)\)
比率 \(\rho_k = \frac{\text{ared}}{\text{pred}}\)
更新规则：
- 若 \(\rho_k > 0.75\)，说明模型拟合好，扩大信赖域：\(\Delta_{k+1} = \min(2\Delta_k, \Delta_{\max})\)
- 若 \(0.25 \leq \rho_k \leq 0.75\)，保持 \(\Delta_{k+1} = \Delta_k\)
- 若 \(\rho_k < 0.25\)，缩小信赖域：\(\Delta_{k+1} = 0.5 \Delta_k\)
如果 \(\rho_k > 0.1\)，接受步长：\(x^{(k+1)} = x^+\)；否则拒绝步长，保持 \(x^{(k)}\) 不变，仅缩小信赖域。
7. 障碍参数更新与收敛判断
每完成一定迭代次数（如5次）或达到子问题收敛后，减小障碍参数：\(\mu_{\text{new}} = 0.1 \mu\)。这逐渐降低障碍项的影响，使迭代点逼近原问题的最优解。
收敛条件：当 \(\mu < \epsilon\) 且 \(\|\nabla f(x^{(k)}) + \sum \lambda_i \nabla c_i(x^{(k)})\| < \epsilon\)，其中 \(\lambda_i = \mu / (-c_i(x))\) 是拉格朗日乘子估计。此时满足 KKT 条件，算法停止。
8. 针对本问题的数值演示（概要）
从 \(x^{(0)} = (1,1)\) 开始，\(\mu=1\)，\(\Delta_0=0.5\)。
第一步：计算 \(B\) 的梯度，由于点靠近直线边界 \(c_2=0\)（因为 \(x_1+x_2-1=1\)，值较小），反射梯度方向会偏向沿边界切向。在信赖域内求 \(\alpha^*\)，得到新点。
迭代过程会沿边界移动，最终趋近最优点。本例理论最优点是约束 \(c_1=0\) 和 \(c_2=0\) 的交点中使 \(f\) 最小的点，即解方程：
\(x_1^2+x_2^2=9\)，\(x_1+x_2=1\)，代入得 \(2x_1^2 - 2x_1 - 8=0\)，解出 \(x_1=2.5\)，\(x_2=-1.5\) 或 \(x_1=-1.5\)，\(x_2=2.5\)。计算 \(f\) 可知 \((2.5, -1.5)\) 更接近 (3,4)，但需检查可行性：\(c_2(2.5,-1.5) = -2.5+1.5+1=0\)，可行。进一步比较无约束最优点 (3,4)（不可行），因此约束最优点就是 \((2.5,-1.5)\)。
算法会从内部逼近该点，最终 \(\mu \to 0\)，梯度条件满足。
9. 算法特点总结
- 内点性：迭代点始终严格可行，适合要求可行解的应用。
- 反射机制：避免边界上的锯齿现象，加速收敛。
- 信赖域：保证每一步的可靠性，特别适用于非凸或曲率变化大的问题。
- 混合策略：兼具全局收敛性与边界处理能力，比单一方法更鲁棒。
通过以上步骤，你可以理解该混合算法的运作机制，并能应用到其他非线性约束问题中。