非线性规划中的逐步二次逼近信赖域反射混合算法进阶题
题目描述
考虑一个非线性规划问题：
\[\begin{aligned} \min_{x \in \mathbb{R}^n} \quad & f(x) \\ \text{s.t.} \quad & c_i(x) = 0, \quad i \in \mathcal{E} \\ & c_i(x) \leq 0, \quad i \in \mathcal{I} \end{aligned} \]
其中，目标函数 \(f: \mathbb{R}^n \to \mathbb{R}\) 是连续可微的，约束函数 \(c_i: \mathbb{R}^n \to \mathbb{R}\) 也是连续可微的。\(\mathcal{E}\) 和 \(\mathcal{I}\) 分别表示等式和不等式约束的指标集。该问题可能具有高度的非凸性（目标或约束函数非凸），传统方法如序列二次规划（SQP）可能在迭代过程中遇到Hessian矩阵非正定、子问题不可行或收敛缓慢等问题。本题要求运用一种结合“逐步二次逼近”（也称为序列二次近似）与“信赖域反射”思想的混合算法来求解该问题。这个算法的核心在于，在每一步迭代中，它不仅用二阶模型（二次函数）来近似原始问题，而且利用信赖域来限制步长范围，确保逼近的准确性。同时，它引入“反射”策略来处理约束边界，避免过早陷入不可行域。此方法特别适合处理具有非线性、非凸约束的中等规模优化问题。
解题过程
我们的目标是设计一个迭代算法，在每一步 \(x_k\) 处，构造一个近似的二次规划子问题，在信赖域内求解该子问题，并根据反射策略调整试探步，以决定是否接受该步长并更新迭代点、信赖域半径等参数。
步骤1：构建局部二次逼近模型
在当前迭代点 \(x_k\) 处，对目标函数 \(f\) 和约束函数 \(c_i\) 进行二阶泰勒展开，得到近似二次函数：
\[\begin{aligned} \hat{f}_k(x_k + d) &= f(x_k) + \nabla f(x_k)^\top d + \frac{1}{2} d^\top B_k d \\ \hat{c}_{i,k}(x_k + d) &= c_i(x_k) + \nabla c_i(x_k)^\top d + \frac{1}{2} d^\top H_{i,k} d, \quad i \in \mathcal{E} \cup \mathcal{I} \end{aligned} \]
其中：
- \(d\) 是待求的步长向量。
- \(B_k\) 是目标函数 Hessian 矩阵 \(\nabla^2 f(x_k)\) 或其近似（例如 BFGS 更新得到的拟牛顿矩阵），以确保正定性。
- \(H_{i,k}\) 是约束函数 Hessian 矩阵 \(\nabla^2 c_i(x_k)\) 或其近似。在实际应用中，为了简化计算，常常只保留一阶项，即令 \(H_{i,k} = 0\)，这就退化为标准的 SQP 子问题。但在此进阶题中，我们要求保留二阶项以获得更好的局部逼近精度，这被称为“逐步二次逼近”（Sequential Quadratic Approximation, SQA）。
步骤2：定义信赖域子问题
为了控制近似模型的有效范围，我们引入信赖域约束 \(\|d\|_{\infty} \leq \Delta_k\)，其中 \(\Delta_k > 0\) 是当前信赖域半径。于是，在 \(x_k\) 处的局部逼近子问题为：
\[\begin{aligned} \min_{d \in \mathbb{R}^n} \quad & \hat{f}_k(x_k + d) \\ \text{s.t.} \quad & \hat{c}_{i,k}(x_k + d) = 0, \quad i \in \mathcal{E} \\ & \hat{c}_{i,k}(x_k + d) \leq 0, \quad i \in \mathcal{I} \\ & \|d\|_{\infty} \leq \Delta_k \end{aligned} \]
这是一个带二次约束的二次规划（QCQP）。由于约束是凸的（如果 \(H_{i,k}\) 正定），但目标函数中的 \(B_k\) 可能非正定，该子问题可能非凸，求解难度较大。一种常见的简化是忽略约束的二次项（即令 \(H_{i,k}=0\)），使其变为一个带线性约束的二次规划（QP），这就是标准的信赖域 SQP 子问题。但在本算法中，我们为了更精确的逼近，可以选择性地对某些关键约束保留二次项，或者通过正则化确保 \(B_k\) 正定，从而使子问题是一个凸的 QCQP，可用内点法有效求解。
步骤3：引入反射步处理边界
信赖域反射（Trust Region Reflection）是一种处理边界约束的策略，但其思想可推广到一般约束。基本想法是：当试探步 \(d_k\) 导致新的试验点 \(x_k + d_k\) 违反某个约束时，不直接拒绝该点，而是尝试对其进行“反射”，即沿着约束边界的法向进行反射，得到一个修正点 \(x_k^+\)，使其回到可行域或至少减少约束违反度。
具体操作如下：
- 求解上述信赖域子问题，得到试探步 \(d_k\)。
- 计算试验点 \(x_k^+ = x_k + d_k\)。
- 如果 \(x_k^+\) 满足所有约束，则直接进入下一步。否则，对违反的约束（例如 \(c_i(x_k^+) > 0\) 对于某个 \(i \in \mathcal{I}\)），计算一个反射步。以简单边界约束 \(l \leq x \leq u\) 为例，若某个分量 \(x_j^+ < l_j\)，则将其反射为 \(x_j^{ref} = 2l_j - x_j^+\)，并确保其不超过上界 \(u_j\)。对于一般非线性约束，反射操作更复杂，常用的是沿约束梯度方向进行投影或反射，使其回到可行域边界附近。一种简化方法是，对违反的约束引入一个惩罚项，将原问题转化为一个无约束或简单约束问题，但这里我们采用几何反射概念：令违反约束的指标集为 \(I_v\)，计算反射点：
\[x_k^{ref} = x_k + d_k - 2 \sum_{i \in I_v} \lambda_i \nabla c_i(x_k) \]
其中，\(\lambda_i\) 是拉格朗日乘子估计，可通过求解子问题的 KKT 条件得到。反射步的目的是“反弹”回可行域内部。
- 将反射后的点 \(x_k^{ref}\) 作为新的试验点，重新评估约束违反度。如果仍不可行，可重复反射或缩小信赖域半径。
步骤4：计算实际下降量与预测下降量之比，更新信赖域半径
定义实际下降量（Actual Reduction, Ared）和预测下降量（Predicted Reduction, Pred）：
\[\begin{aligned} \text{Ared}_k &= f(x_k) - f(x_k^+) + \rho \left( \|c(x_k)\|_1 - \|c(x_k^+)\|_1 \right) \\ \text{Pred}_k &= \hat{f}_k(x_k) - \hat{f}_k(x_k^+) + \rho \left( \|c(x_k)\|_1 - \|\hat{c}(x_k^+)\|_1 \right) \end{aligned} \]
其中，\(\rho > 0\) 是惩罚参数，用于平衡目标函数下降和约束违反度的减少。\(\|c(x)\|_1 = \sum_{i \in \mathcal{E}} |c_i(x)| + \sum_{i \in \mathcal{I}} \max(0, c_i(x))\) 是约束违反的1-范数度量。这里 \(\hat{c}(x_k^+)\) 是二次逼近模型在试验点的约束值，而 \(c(x_k^+)\) 是原始约束的实际值。
然后，计算比值：
\[r_k = \frac{\text{Ared}_k}{\text{Pred}_k} \]
这个比值衡量了二次逼近模型的准确度：
- 如果 \(r_k\) 接近1（例如 \(r_k > 0.75\)），说明模型拟合很好，可以接受该步长，并可能增大信赖域半径：\(\Delta_{k+1} = \min(2\Delta_k, \Delta_{\max})\)。
- 如果 \(r_k\) 在中等范围（例如 \(0.1 \leq r_k \leq 0.75\)），说明模型一般，接受步长但保持信赖域半径不变：\(\Delta_{k+1} = \Delta_k\)。
- 如果 \(r_k\) 很小（例如 \(r_k < 0.1\)），说明模型拟合很差，拒绝该步长，缩小信赖域半径：\(\Delta_{k+1} = 0.5 \Delta_k\)，并重新求解子问题（或反射后重试）。
步骤5：更新迭代点与Hessian近似
如果步长被接受（即 \(r_k \geq 0.1\)），则更新迭代点：\(x_{k+1} = x_k^+\)。否则，保持 \(x_{k+1} = x_k\)。
接下来，更新目标函数的Hessian近似矩阵 \(B_k\)。通常使用拟牛顿法（如BFGS）更新公式：
\[B_{k+1} = B_k + \frac{y_k y_k^\top}{y_k^\top s_k} - \frac{B_k s_k s_k^\top B_k}{s_k^\top B_k s_k} \]
其中，\(s_k = x_{k+1} - x_k\)，\(y_k = \nabla_x L(x_{k+1}, \lambda_{k+1}) - \nabla_x L(x_k, \lambda_{k+1})\)，这里 \(L(x, \lambda) = f(x) + \sum_{i \in \mathcal{E} \cup \mathcal{I}} \lambda_i c_i(x)\) 是拉格朗日函数，\(\lambda_{k+1}\) 是当前乘子估计（来自于子问题的解）。这样可以确保 \(B_{k+1}\) 满足拟牛顿方程，并保持正定性。
步骤6：判断收敛
检查是否满足收敛条件。常见的收敛准则包括：
- 梯度条件：原始可行性（约束违反度）和对偶可行性（梯度条件）的范数小于给定容差。
\[ \|c(x_k)\|_{\infty} \leq \epsilon_c, \quad \|\nabla_x L(x_k, \lambda_k)\|_{\infty} \leq \epsilon_g \]
- 步长条件：步长大小足够小，\(\|x_{k+1} - x_k\|_{\infty} \leq \epsilon_x\)。
- 目标变化：\(|f(x_{k+1}) - f(x_k)| \leq \epsilon_f\)。
如果满足收敛条件，则算法终止，输出 \(x_{k+1}\) 作为近似最优解。否则，令 \(k := k+1\)，返回步骤1继续迭代。
步骤7：整体算法流程总结
- 初始化：给定初始点 \(x_0\)，初始信赖域半径 \(\Delta_0 > 0\)，初始 Hessian 近似 \(B_0 = I\)（单位矩阵），惩罚参数 \(\rho\)，收敛容差 \(\epsilon_c, \epsilon_g, \epsilon_x, \epsilon_f\)，最大迭代次数 \(K\)。设 \(k=0\)。
- 构造子问题：在当前点 \(x_k\)，构建二次逼近模型 \(\hat{f}_k, \hat{c}_{i,k}\)，形成带信赖域约束的 QCQP 子问题。
- 求解子问题：使用凸优化求解器（如内点法）求解该子问题，得到步长 \(d_k\) 和拉格朗日乘子估计 \(\lambda_k\)。
- 反射处理：计算试验点 \(x_k^+ = x_k + d_k\)，检查约束违反。若有违反，进行反射操作得到修正点 \(x_k^{ref}\)，令 \(x_k^+ := x_k^{ref}\)。重复此步骤直到满足可行性或达到最大反射次数。
- 计算比值：计算实际下降量与预测下降量之比 \(r_k\)。
- 更新信赖域半径：根据 \(r_k\) 的值，按步骤4的规则调整 \(\Delta_{k+1}\)。
- 接受或拒绝步长：若 \(r_k \geq 0.1\)，则令 \(x_{k+1} = x_k^+\)；否则，令 \(x_{k+1} = x_k\)。
- 更新Hessian近似：利用 BFGS 公式更新 \(B_{k+1}\)。
- 检查收敛：若满足收敛条件或 \(k \geq K\)，则停止；否则，令 \(k := k+1\)，返回步骤2。
总结
这个“逐步二次逼近信赖域反射混合算法”的核心优势在于：它通过二阶模型提高了逼近精度，利用信赖域控制步长可靠性，再通过反射机制处理约束边界，增强了算法的鲁棒性和收敛速度。它特别适合处理非凸、非线性约束问题，其中反射步骤能有效避免算法在约束边界附近振荡或陷入不可行域。在实现时，需要仔细调整参数（如初始信赖域半径、惩罚参数 \(\rho\)、反射策略的系数等），并对子问题的求解选择合适的凸优化算法，以确保整体效率。