非线性规划中的序列二次规划-积极集-乘子法混合算法进阶题：带非凸不等式与线性等式约束的大规模稀疏优化问题
字数 4585 2025-12-20 08:25:24

非线性规划中的序列二次规划-积极集-乘子法混合算法进阶题：带非凸不等式与线性等式约束的大规模稀疏优化问题

题目描述

考虑一个具有大规模稀疏结构的非线性规划问题，其一般形式如下：

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

其中：

  • 目标函数 \(f: \mathbb{R}^n \to \mathbb{R}\) 是二阶连续可微的，可能是非凸的。
  • 不等式约束 \(c_i: \mathbb{R}^n \to \mathbb{R}\) 也是二阶连续可微的，可能是非凸的，且数量 \(m\) 很大（例如 \(m \sim 10^4-10^6\)）。
  • 线性等式约束 \(A x = b\) 中，矩阵 \(A \in \mathbb{R}^{p \times n}\) 是稀疏的，且 \(p \ll n\)。
  • 问题的总体结构是大规模且稀疏的，这意味着目标函数和约束的Hessian矩阵和Jacobian矩阵具有大量零元素，适合利用稀疏线性代数技术。

要求：设计一个结合序列二次规划（SQP）、积极集法（Active Set）和乘子法（Method of Multipliers）的混合算法，以高效求解该问题。算法需要在每一步求解一个中等规模的二次规划子问题，并自适应地调整积极集和乘子估计，同时利用问题的稀疏性以降低计算成本。

解题过程

我们将循序渐进地设计并讲解这个混合算法的核心步骤。

步骤1：问题重构与增广拉格朗日函数

为了处理等式和不等式约束，我们首先引入松弛变量 \(s \geq 0\) 将不等式转化为等式：\(c_i(x) + s_i = 0\)。但为保持算法效率，我们更常用增广拉格朗日方法。定义增广拉格朗日函数（仅对不等式约束进行增广，等式约束用乘子直接处理）：

\[\mathcal{L}(x, \lambda, \mu; \rho) = f(x) + \lambda^\top (Ax - b) + \sum_{i=1}^m \left[ \mu_i c_i(x) + \frac{\rho}{2} \max(0, c_i(x) + \frac{\mu_i}{\rho})^2 \right], \]

其中：

  • \(\lambda \in \mathbb{R}^p\) 是对应线性等式约束的拉格朗日乘子。
  • \(\mu \in \mathbb{R}^m\) 是对应不等式约束的乘子，且要求 \(\mu \geq 0\)。
  • \(\rho > 0\) 是罚参数。

注：增广项是针对不等式约束的二次罚函数形式，它使得在迭代中约束违反较大时惩罚增强，帮助引导迭代点走向可行域。

步骤2：序列二次规划（SQP）子问题构造

在每次迭代 \(k\)，给定当前点 \(x_k\)、乘子估计 \(\lambda_k\) 和 \(\mu_k\)，以及罚参数 \(\rho_k\)，我们构造一个二次规划（QP）子问题来近似原问题。这个QP基于增广拉格朗日函数的二阶近似。

  1. 梯度计算：

    • 计算目标梯度：\(g_k = \nabla f(x_k)\)。
    • 计算等式约束Jacobian：\(A\)（常数，无需重算）。
    • 计算不等式约束梯度：\(a_i = \nabla c_i(x_k)\)，\(i=1,\dots,m\)。

  2. Hessian近似：
    由于原问题可能非凸，我们使用拟牛顿法（如BFGS）来维护一个正定近似 \(B_k \approx \nabla^2_{xx} \mathcal{L}(x_k, \lambda_k, \mu_k; \rho_k)\)。利用稀疏结构，我们只存储和更新Hessian的非零元素。

  3. QP子问题：
    令搜索方向为 \(d = x - x_k\)。子问题为：

\[ \begin{aligned} \min_{d \in \mathbb{R}^n} \quad & \frac{1}{2} d^\top B_k d + g_k^\top d \\ \text{s.t.} \quad & A d = 0, \quad \text{(来自线性等式约束的线性化)} \\ & c_i(x_k) + a_i^\top d \leq 0, \quad i=1,\dots,m. \end{aligned} \]

这是一个带线性等式和不等式约束的QP。由于 \(m\) 很大，直接求解这个QP代价高昂。因此，我们需要结合积极集法来减少每次迭代中活跃的约束数量。

步骤3：积极集策略

积极集法的核心思想是，在最优解处，只有一部分不等式约束是起作用的（即 \(c_i(x^*) = 0\)），其余的可以忽略。在迭代中，我们预测一个“工作集” \(\mathcal{W}_k \subseteq \{1,\dots,m\}\)，它包含可能活跃的约束。

  1. 初始工作集：可以根据当前点 \(x_k\) 的约束违反程度来选择。例如，对于某个阈值 \(\epsilon > 0\)，定义

\[ \mathcal{W}_k = \{ i \mid c_i(x_k) \geq -\epsilon \}. \]

由于问题稀疏，许多约束函数 \(c_i\) 只依赖于 \(x\) 的少部分分量，因此工作集大小远小于 \(m\)。

  1. 简化QP：
    只考虑工作集中的约束，得到简化QP：

\[ \begin{aligned} \min_{d} \quad & \frac{1}{2} d^\top B_k d + g_k^\top d \\ \text{s.t.} \quad & A d = 0, \\ & c_i(x_k) + a_i^\top d \leq 0, \quad i \in \mathcal{W}_k. \end{aligned} \]

这个QP的规模大大减小，可以利用稀疏QP求解器（如基于内点法或积极集法的稀疏求解器）高效求解。

  1. 工作集更新：
    求解简化QP得到方向 \(d_k\)。然后检查是否有一些不在工作集中的约束被违反（即 \(c_i(x_k + d_k) > 0\)）。如果是，将这些约束加入工作集，重新求解QP；否则，当前工作集可能合适。同时，如果工作集中某个约束对应的乘子估计 \(\mu_i\) 为负（或接近零），且约束不活跃，则可以从工作集中移除。

步骤4：乘子更新与步长接受

得到QP的解 \(d_k\) 后，我们更新迭代点：

\[x_{k+1} = x_k + \alpha_k d_k, \]

其中步长 \(\alpha_k \in (0,1]\) 通过线搜索确定，以确保足够的下降和可行性改善。我们使用增广拉格朗日函数作为价值函数（merit function）：

\[\phi(x; \rho) = f(x) + \lambda^\top (Ax - b) + \frac{\rho}{2} \sum_{i=1}^m \max\left(0, c_i(x) + \frac{\mu_i}{\rho}\right)^2. \]

线搜索条件（Armijo条件）为：

\[\phi(x_{k+1}; \rho_k) \leq \phi(x_k; \rho_k) + \eta \alpha_k D\phi(d_k), \]

其中 \(D\phi(d_k)\) 是价值函数在方向 \(d_k\) 的方向导数，\(\eta \in (0,1)\) 是常数。

然后，更新乘子：

  • 等式约束乘子：\(\lambda_{k+1} = \lambda_k + \rho_k (A x_{k+1} - b)\)。
  • 不等式约束乘子：\(\mu_{i,k+1} = \max\left(0, \mu_{i,k} + \rho_k c_i(x_{k+1})\right)\)，\(i=1,\dots,m\)。

步骤5：自适应调整罚参数

为了平衡约束满足和最优性，我们自适应调整罚参数 \(\rho\)。一个常见策略是基于约束违反的变化：令约束违反量为

\[v_k = \|Ax_k - b\| + \left\| \max(0, c(x_k)) \right\|. \]

如果 \(v_k\) 相比上一步没有足够减少（例如，减少因子小于 \(\tau \in (0,1)\)），则增加罚参数：\(\rho_{k+1} = \beta \rho_k\)，其中 \(\beta > 1\)；否则保持 \(\rho_{k+1} = \rho_k\)。这迫使迭代点更快地趋向可行域。

步骤6：算法流程与收敛条件

将以上步骤整合为完整算法：

  1. 初始化：给定初始点 \(x_0\)，乘子 \(\lambda_0, \mu_0 \geq 0\)，罚参数 \(\rho_0 > 0\)，容忍度 \(\epsilon > 0\)，以及Hessian近似 \(B_0 = I\)。
  2. 主循环（对于 \(k=0,1,2,\dots\)）：
    a. 计算梯度 \(g_k\) 和约束值 \(c_i(x_k)\)。
    b. 根据阈值选择工作集 \(\mathcal{W}_k\)。
    c. 求解简化QP得到方向 \(d_k\)。
    d. 执行线搜索确定步长 \(\alpha_k\)，更新 \(x_{k+1} = x_k + \alpha_k d_k\)。
    e. 更新乘子 \(\lambda_{k+1}, \mu_{k+1}\) 和罚参数 \(\rho_{k+1}\)。
    f. 更新Hessian近似 \(B_{k+1}\) 使用拟牛顿公式（例如稀疏BFGS）。
    g. 检查收敛条件：如果 \(\|d_k\| < \epsilon\) 且约束违反 \(v_k < \epsilon\)，则停止。
  3. 输出：近似最优解 \(x_{k+1}\) 和乘子 \(\lambda_{k+1}, \mu_{k+1}\)。

步骤7：利用稀疏性

由于问题是大规模稀疏的，算法实现中需注意：

  • 只计算和存储梯度、Jacobian和Hessian的非零元素。
  • 使用稀疏矩阵格式（如CSR）来存储 \(A\) 和约束梯度。
  • 在求解QP时，利用其稀疏结构，使用稀疏QP求解器（如OSQP、qpOASES的稀疏版本）。
  • 拟牛顿更新也采用稀疏更新策略，只更新Hessian中与非零梯度分量对应的部分。

总结

本算法结合了SQP的局部快速收敛、积极集法对大规模约束的降维处理、乘子法的全局收敛保证，并通过自适应调整罚参数和稀疏计算，高效求解大规模稀疏非凸约束优化问题。其核心优势在于每一步只处理一个中等规模的QP，且利用稀疏性显著降低了计算和存储成本。

非线性规划中的序列二次规划-积极集-乘子法混合算法进阶题：带非凸不等式与线性等式约束的大规模稀疏优化问题 题目描述 考虑一个具有大规模稀疏结构的非线性规划问题，其一般形式如下： \[\begin{aligned} \min_ {x \in \mathbb{R}^n} \quad & f(x) \\ \text{s.t.} \quad & c_ i(x) \leq 0, \quad i=1,\dots,m, \\ & Ax = b, \end{aligned}\] 其中： 目标函数 \(f: \mathbb{R}^n \to \mathbb{R}\) 是二阶连续可微的，可能是非凸的。 不等式约束 \(c_ i: \mathbb{R}^n \to \mathbb{R}\) 也是二阶连续可微的，可能是非凸的，且数量 \(m\) 很大（例如 \(m \sim 10^4-10^6\)）。 线性等式约束 \(A x = b\) 中，矩阵 \(A \in \mathbb{R}^{p \times n}\) 是稀疏的，且 \(p \ll n\)。 问题的总体结构是大规模且稀疏的，这意味着目标函数和约束的Hessian矩阵和Jacobian矩阵具有大量零元素，适合利用稀疏线性代数技术。 要求 ：设计一个结合序列二次规划（SQP）、积极集法（Active Set）和乘子法（Method of Multipliers）的混合算法，以高效求解该问题。算法需要在每一步求解一个中等规模的二次规划子问题，并自适应地调整积极集和乘子估计，同时利用问题的稀疏性以降低计算成本。 解题过程 我们将循序渐进地设计并讲解这个混合算法的核心步骤。 步骤1：问题重构与增广拉格朗日函数 为了处理等式和不等式约束，我们首先引入松弛变量 \(s \geq 0\) 将不等式转化为等式：\(c_ i(x) + s_ i = 0\)。但为保持算法效率，我们更常用增广拉格朗日方法。定义增广拉格朗日函数（仅对不等式约束进行增广，等式约束用乘子直接处理）： \[ \mathcal{L}(x, \lambda, \mu; \rho) = f(x) + \lambda^\top (Ax - b) + \sum_ {i=1}^m \left[ \mu_ i c_ i(x) + \frac{\rho}{2} \max(0, c_ i(x) + \frac{\mu_ i}{\rho})^2 \right ], \] 其中： \(\lambda \in \mathbb{R}^p\) 是对应线性等式约束的拉格朗日乘子。 \(\mu \in \mathbb{R}^m\) 是对应不等式约束的乘子，且要求 \(\mu \geq 0\)。 \(\rho > 0\) 是罚参数。 注 ：增广项是针对不等式约束的二次罚函数形式，它使得在迭代中约束违反较大时惩罚增强，帮助引导迭代点走向可行域。 步骤2：序列二次规划（SQP）子问题构造 在每次迭代 \(k\)，给定当前点 \(x_ k\)、乘子估计 \(\lambda_ k\) 和 \(\mu_ k\)，以及罚参数 \(\rho_ k\)，我们构造一个二次规划（QP）子问题来近似原问题。这个QP基于增广拉格朗日函数的二阶近似。 梯度计算 ： 计算目标梯度：\(g_ k = \nabla f(x_ k)\)。 计算等式约束Jacobian：\(A\)（常数，无需重算）。 计算不等式约束梯度：\(a_ i = \nabla c_ i(x_ k)\)，\(i=1,\dots,m\)。 Hessian近似 ： 由于原问题可能非凸，我们使用拟牛顿法（如BFGS）来维护一个正定近似 \(B_ k \approx \nabla^2_ {xx} \mathcal{L}(x_ k, \lambda_ k, \mu_ k; \rho_ k)\)。利用稀疏结构，我们只存储和更新Hessian的非零元素。 QP子问题 ： 令搜索方向为 \(d = x - x_ k\)。子问题为： \[ \begin{aligned} \min_ {d \in \mathbb{R}^n} \quad & \frac{1}{2} d^\top B_ k d + g_ k^\top d \\ \text{s.t.} \quad & A d = 0, \quad \text{(来自线性等式约束的线性化)} \\ & c_ i(x_ k) + a_ i^\top d \leq 0, \quad i=1,\dots,m. \end{aligned} \] 这是一个带线性等式和不等式约束的QP。由于 \(m\) 很大，直接求解这个QP代价高昂。因此，我们需要结合积极集法来减少每次迭代中活跃的约束数量。 步骤3：积极集策略 积极集法的核心思想是，在最优解处，只有一部分不等式约束是起作用的（即 \(c_ i(x^* ) = 0\)），其余的可以忽略。在迭代中，我们预测一个“工作集” \(\mathcal{W}_ k \subseteq \{1,\dots,m\}\)，它包含可能活跃的约束。 初始工作集 ：可以根据当前点 \(x_ k\) 的约束违反程度来选择。例如，对于某个阈值 \(\epsilon > 0\)，定义 \[ \mathcal{W}_ k = \{ i \mid c_ i(x_ k) \geq -\epsilon \}. \] 由于问题稀疏，许多约束函数 \(c_ i\) 只依赖于 \(x\) 的少部分分量，因此工作集大小远小于 \(m\)。 简化QP ： 只考虑工作集中的约束，得到简化QP： \[ \begin{aligned} \min_ {d} \quad & \frac{1}{2} d^\top B_ k d + g_ k^\top d \\ \text{s.t.} \quad & A d = 0, \\ & c_ i(x_ k) + a_ i^\top d \leq 0, \quad i \in \mathcal{W}_ k. \end{aligned} \] 这个QP的规模大大减小，可以利用稀疏QP求解器（如基于内点法或积极集法的稀疏求解器）高效求解。 工作集更新 ： 求解简化QP得到方向 \(d_ k\)。然后检查是否有一些不在工作集中的约束被违反（即 \(c_ i(x_ k + d_ k) > 0\)）。如果是，将这些约束加入工作集，重新求解QP；否则，当前工作集可能合适。同时，如果工作集中某个约束对应的乘子估计 \(\mu_ i\) 为负（或接近零），且约束不活跃，则可以从工作集中移除。 步骤4：乘子更新与步长接受 得到QP的解 \(d_ k\) 后，我们更新迭代点： \[ x_ {k+1} = x_ k + \alpha_ k d_ k, \] 其中步长 \(\alpha_ k \in (0,1 ]\) 通过线搜索确定，以确保足够的下降和可行性改善。我们使用增广拉格朗日函数作为价值函数（merit function）： \[ \phi(x; \rho) = f(x) + \lambda^\top (Ax - b) + \frac{\rho}{2} \sum_ {i=1}^m \max\left(0, c_ i(x) + \frac{\mu_ i}{\rho}\right)^2. \] 线搜索条件（Armijo条件）为： \[ \phi(x_ {k+1}; \rho_ k) \leq \phi(x_ k; \rho_ k) + \eta \alpha_ k D\phi(d_ k), \] 其中 \(D\phi(d_ k)\) 是价值函数在方向 \(d_ k\) 的方向导数，\(\eta \in (0,1)\) 是常数。 然后，更新乘子： 等式约束乘子：\(\lambda_ {k+1} = \lambda_ k + \rho_ k (A x_ {k+1} - b)\)。 不等式约束乘子：\(\mu_ {i,k+1} = \max\left(0, \mu_ {i,k} + \rho_ k c_ i(x_ {k+1})\right)\)，\(i=1,\dots,m\)。 步骤5：自适应调整罚参数 为了平衡约束满足和最优性，我们自适应调整罚参数 \(\rho\)。一个常见策略是基于约束违反的变化：令约束违反量为 \[ v_ k = \|Ax_ k - b\| + \left\| \max(0, c(x_ k)) \right\|. \] 如果 \(v_ k\) 相比上一步没有足够减少（例如，减少因子小于 \(\tau \in (0,1)\)），则增加罚参数：\(\rho_ {k+1} = \beta \rho_ k\)，其中 \(\beta > 1\)；否则保持 \(\rho_ {k+1} = \rho_ k\)。这迫使迭代点更快地趋向可行域。 步骤6：算法流程与收敛条件 将以上步骤整合为完整算法： 初始化 ：给定初始点 \(x_ 0\)，乘子 \(\lambda_ 0, \mu_ 0 \geq 0\)，罚参数 \(\rho_ 0 > 0\)，容忍度 \(\epsilon > 0\)，以及Hessian近似 \(B_ 0 = I\)。 主循环 （对于 \(k=0,1,2,\dots\)）： a. 计算梯度 \(g_ k\) 和约束值 \(c_ i(x_ k)\)。 b. 根据阈值选择工作集 \(\mathcal{W} k\)。 c. 求解简化QP得到方向 \(d_ k\)。 d. 执行线搜索确定步长 \(\alpha_ k\)，更新 \(x {k+1} = x_ k + \alpha_ k d_ k\)。 e. 更新乘子 \(\lambda_ {k+1}, \mu_ {k+1}\) 和罚参数 \(\rho_ {k+1}\)。 f. 更新Hessian近似 \(B_ {k+1}\) 使用拟牛顿公式（例如稀疏BFGS）。 g. 检查收敛条件：如果 \(\|d_ k\| < \epsilon\) 且约束违反 \(v_ k < \epsilon\)，则停止。 输出 ：近似最优解 \(x_ {k+1}\) 和乘子 \(\lambda_ {k+1}, \mu_ {k+1}\)。 步骤7：利用稀疏性 由于问题是大规模稀疏的，算法实现中需注意： 只计算和存储梯度、Jacobian和Hessian的非零元素。 使用稀疏矩阵格式（如CSR）来存储 \(A\) 和约束梯度。 在求解QP时，利用其稀疏结构，使用稀疏QP求解器（如OSQP、qpOASES的稀疏版本）。 拟牛顿更新也采用稀疏更新策略，只更新Hessian中与非零梯度分量对应的部分。 总结 本算法结合了SQP的局部快速收敛、积极集法对大规模约束的降维处理、乘子法的全局收敛保证，并通过自适应调整罚参数和稀疏计算，高效求解大规模稀疏非凸约束优化问题。其核心优势在于每一步只处理一个中等规模的QP，且利用稀疏性显著降低了计算和存储成本。