非线性规划中的序列二次约束二次规划(Sequential Quadratically Constrained Quadratic Programming, SQCQP)算法基础题
题目描述
考虑一个具有非线性目标函数和非线性约束的一般形式优化问题:
\[\begin{aligned} \min_{x \in \mathbb{R}^n} \quad & f(x) \\ \text{s.t.} \quad & c_i(x) \le 0, \quad i = 1, 2, \ldots, m \\ & h_j(x) = 0, \quad j = 1, 2, \ldots, p \end{aligned} \]
其中,\(f(x)\) 是二次连续可微的非线性目标函数,\(c_i(x)\) 是二次连续可微的非线性不等式约束,\(h_j(x)\) 是二次连续可微的非线性等式约束。我们的目标是找到满足约束条件并使目标函数最小化的变量 \(x^*\)。
我们要求使用序列二次约束二次规划(SQCQP) 算法来求解此问题。SQCQP 算法是序列二次规划(SQP)方法的一种扩展,其特点是在每一步迭代中,不仅用二次模型近似目标函数,还将非线性约束近似为二次约束,从而在每一步求解一个二次约束二次规划(QCQP)子问题。本题将引导你理解SQCQP的基本框架、推导子问题构造、并完成算法实现步骤。
解题过程
第一步:理解SQCQP的基本思想
SQP方法通过在当前迭代点 \(x_k\) 处,用二次模型(泰勒展开到二阶)近似目标函数,用线性模型(一阶泰勒展开)近似约束函数,构造一个二次规划(QP)子问题。SQCQP在此基础上更进一步:它意识到对于高度非线性的约束,线性近似可能不够准确,导致子问题的可行域与原问题相差较大,进而使算法收敛缓慢甚至失败。因此,SQCQP将非线性约束用二次模型(泰勒展开到二阶)近似,从而每一步求解的是一个QCQP子问题,其可行域能更好地逼近原问题的可行域,尤其适用于约束曲率较大的问题。
核心迭代格式:在点 \(x_k\) 处,构造如下QCQP子问题:
\[\begin{aligned} \min_{d \in \mathbb{R}^n} \quad & f(x_k) + \nabla f(x_k)^T d + \frac{1}{2} d^T H_k d \\ \text{s.t.} \quad & c_i(x_k) + \nabla c_i(x_k)^T d + \frac{1}{2} d^T \nabla^2 c_i(x_k) d \le 0, \quad i=1,\ldots,m \\ & h_j(x_k) + \nabla h_j(x_k)^T d + \frac{1}{2} d^T \nabla^2 h_j(x_k) d = 0, \quad j=1,\ldots,p \end{aligned} \]
其中:
- \(d = x - x_k\) 是搜索方向。
- \(H_k\) 是拉格朗日函数的Hessian或其近似(例如BFGS更新)。
- 约束中的二阶项 \(\frac{1}{2} d^T \nabla^2 c_i(x_k) d\) 和 \(\frac{1}{2} d^T \nabla^2 h_j(x_k) d\) 使得约束近似为二次(可能非凸)。
第二步:构造拉格朗日函数与Hessian近似
为了构造目标函数的二次模型,我们通常使用拉格朗日函数的Hessian,因为它能同时考虑目标函数和约束的曲率信息。定义拉格朗日函数:
\[\mathcal{L}(x, \lambda, \mu) = f(x) + \sum_{i=1}^m \lambda_i c_i(x) + \sum_{j=1}^p \mu_j h_j(x) \]
其中 \(\lambda_i \ge 0\) 和 \(\mu_j\) 是对应约束的拉格朗日乘子。
在SQCQP中,我们取 \(H_k\) 为 \(\nabla_{xx}^2 \mathcal{L}(x_k, \lambda_k, \mu_k)\) 或其拟牛顿近似(例如BFGS更新)。使用精确Hessian能提供二阶收敛速度,但计算成本高;拟牛顿近似能降低计算量,具有超线性收敛。
Hessian更新(BFGS):若使用拟牛顿法,则 \(H_k\) 通过BFGS公式更新:
\[H_{k+1} = H_k - \frac{H_k s_k s_k^T H_k}{s_k^T H_k s_k} + \frac{y_k y_k^T}{s_k^T y_k} \]
其中 \(s_k = x_{k+1} - x_k\),\(y_k = \nabla_x \mathcal{L}(x_{k+1}, \lambda_{k+1}, \mu_{k+1}) - \nabla_x \mathcal{L}(x_k, \lambda_{k+1}, \mu_{k+1})\)。
第三步:求解QCQP子问题
子问题是一个二次约束二次规划(QCQP)。由于约束是二次的,子问题可能非凸,因此求解比QP更复杂。常用求解方法包括:
- 内点法:适用于中小规模问题,能处理非凸QCQP,通过引入屏障函数将问题转化为一系列无约束或仅含等式约束的问题。
- 序列凸近似(SCA):在QCQP内部进一步迭代,将二次约束线性化,转化为一系列QP子问题求解。
- 增广拉格朗日法:将QCQP转化为一系列无约束问题求解。
简化处理:在基础题中,为简化,我们假设在 \(x_k\) 附近,约束的二次项较小,子问题仍是凸的,或使用信赖域技巧控制步长 \(d\) 的大小,保证子问题可解。实际中,可添加信赖域约束 \(\|d\| \le \Delta_k\) 将子问题变为凸的。
第四步:线搜索与信赖域策略
得到搜索方向 \(d_k\) 后,需确定步长 \(\alpha_k\) 以获得新迭代点 \(x_{k+1} = x_k + \alpha_k d_k\)。
线搜索:使用Armijo条件或 Wolfe 条件,确保目标函数充分下降并满足约束可行性。例如,使用 \(l_1\) 精确罚函数:
\[\Phi_1(x; \rho) = f(x) + \rho \left( \sum_{i=1}^m \max(0, c_i(x)) + \sum_{j=1}^p |h_j(x)| \right) \]
其中 \(\rho > 0\) 是罚参数。线搜索要求 \(\Phi_1(x_{k+1}; \rho) \le \Phi_1(x_k; \rho) + \eta \alpha_k D \Phi_1(x_k, d_k; \rho)\),其中 \(D \Phi_1\) 是方向导数。
信赖域:在子问题中加入约束 \(\|d\| \le \Delta_k\),控制近似模型的精度。根据实际下降量与预测下降量的比值 \(\rho_k\) 调整信赖域半径 \(\Delta_k\):
\[\rho_k = \frac{\text{实际下降量}}{\text{预测下降量}} = \frac{f(x_k) - f(x_k + d_k)}{-[\nabla f(x_k)^T d_k + \frac{1}{2} d_k^T H_k d_k]} \]
若 \(\rho_k\) 接近1,模型拟合好,增大 \(\Delta_k\);若 \(\rho_k\) 很小,减小 \(\Delta_k\)。
第五步:更新拉格朗日乘子
在得到新迭代点 \(x_{k+1}\) 后,需更新乘子 \(\lambda\) 和 \(\mu\)。可通过求解子问题的KKT条件获得乘子估计,或使用如下更新公式:
\[\lambda_i^{k+1} = \max(0, \lambda_i^k + \alpha_k c_i(x_k)), \quad \mu_j^{k+1} = \mu_j^k + \alpha_k h_j(x_k) \]
这实质上是梯度投影法对乘子的更新。
第六步:算法步骤总结
- 初始化:给定初始点 \(x_0\),初始乘子 \(\lambda_0, \mu_0\),初始Hessian近似 \(H_0 = I\)(单位阵),初始信赖域半径 \(\Delta_0 > 0\),设定容忍误差 \(\epsilon > 0\),令 \(k = 0\)。
- 构造QCQP子问题:在当前点 \(x_k\),计算梯度 \(\nabla f(x_k)\)、\(\nabla c_i(x_k)\)、\(\nabla h_j(x_k)\) 和Hessian \(\nabla^2 c_i(x_k)\)、\(\nabla^2 h_j(x_k)\)。用BFGS更新 \(H_k\)(若使用拟牛顿法)。构造子问题:
\[ \begin{aligned} \min_{d} \quad & \nabla f(x_k)^T d + \frac{1}{2} d^T H_k d \\ \text{s.t.} \quad & c_i(x_k) + \nabla c_i(x_k)^T d + \frac{1}{2} d^T \nabla^2 c_i(x_k) d \le 0, \ \forall i \\ & h_j(x_k) + \nabla h_j(x_k)^T d + \frac{1}{2} d^T \nabla^2 h_j(x_k) d = 0, \ \forall j \\ & \|d\| \le \Delta_k \quad \text{(信赖域约束)} \end{aligned} \]
- 求解子问题:使用内点法或序列凸近似求解该QCQP,得到搜索方向 \(d_k\) 及对应的乘子估计 \(\lambda_k^q, \mu_k^q\)。
- 计算步长:通过线搜索确定步长 \(\alpha_k \in (0,1]\),满足 \(\Phi_1(x_k + \alpha_k d_k; \rho) \le \Phi_1(x_k; \rho) + \eta \alpha_k D \Phi_1(x_k, d_k; \rho)\)。
- 更新迭代点:\(x_{k+1} = x_k + \alpha_k d_k\)。
- 更新乘子:\(\lambda_{k+1} = \lambda_k^q\),\(\mu_{k+1} = \mu_k^q\)(或使用梯度投影更新)。
- 更新信赖域半径:计算比值 \(\rho_k\),按如下规则调整 \(\Delta_{k+1}\):
\[ \Delta_{k+1} = \begin{cases} \Delta_k / 2, & \text{if } \rho_k < 0.25 \\ \min(2\Delta_k, \bar{\Delta}), & \text{if } \rho_k > 0.75 \\ \Delta_k, & \text{otherwise} \end{cases} \]
其中 \(\bar{\Delta}\) 是最大半径。
8. 检查收敛:若 \(\|d_k\| < \epsilon\) 且约束违反度小于 \(\epsilon\),则停止;否则令 \(k = k+1\),返回步骤2。
关键点说明
- SQCQP 通过二阶近似约束,能更准确地捕捉可行域边界曲率,适合处理高非线性约束问题。
- 子问题为QCQP,求解比QP复杂,但现代内点法可有效求解中小规模问题。
- 结合信赖域可保证全局收敛性,避免因模型不准确导致迭代发散。
- 拟牛顿更新Hessian可避免计算二阶导数,适用于大规模问题。
此算法是SQP的自然扩展,在工程优化、经济模型等有广泛应用。通过上述步骤,你可逐步实现SQCQP算法求解一般非线性规划问题。