非线性规划中的信赖域反射-序列二次规划混合算法基础题
题目描述
考虑非线性规划问题:
\[\min f(x) = (x_1 - 2)^4 + (x_1 - 2x_2)^2 \]
\[ \text{s.t.} \quad g(x) = x_1^2 - x_2 \le 0 \]
\[ h(x) = x_1 + x_2 - 2 = 0 \]
其中 \(x = (x_1, x_2)^T \in \mathbb{R}^2\)。这是一个包含非线性目标函数、一个非线性不等式约束和一个线性等式约束的优化问题。我们要求使用信赖域反射-序列二次规划混合算法求解该问题。算法融合了信赖域反射法的全局稳定性与序列二次规划法的局部快速收敛性,旨在高效、鲁棒地找到满足约束的局部最优解。给定初始点 \(x^{(0)} = (0, 0)^T\),初始信赖域半径 \(\Delta_0 = 0.5\),信赖域半径更新参数 \(\eta_1 = 0.1\),\(\eta_2 = 0.75\),收缩因子 \(\gamma_1 = 0.5\),扩张因子 \(\gamma_2 = 2.0\),迭代终止容差 \(\epsilon = 10^{-6}\)。你需要阐述该混合算法的详细步骤,并展示至少两次迭代的手算推导过程。
解题过程循序渐进讲解
首先,我们明确一下信赖域反射-序列二次规划混合算法的核心思想:
- 在每一步迭代,我们构造一个二次规划子问题来近似原问题。这个子问题的目标函数是目标函数 \(f(x)\) 的二次近似,约束是约束函数 \(g(x)\) 和 \(h(x)\) 的线性近似。
- 我们使用信赖域来限制子问题中决策变量的变化范围,确保近似只在局部有效,从而保证算法的全局收敛性。
- 在求解子问题时,我们采用反射策略处理边界约束(来自信赖域),这能更有效地利用信赖域空间,尤其在边界附近。
下面,我们分步骤讲解算法流程,并用给定的题目进行演示。
步骤一:算法框架
设当前迭代点为 \(x_k\),当前信赖域半径为 \(\Delta_k\)。
- 构建二次规划子问题:
- 计算目标函数梯度 \(\nabla f(x_k)\) 和 Hessian 矩阵 \(B_k\)(或用拟牛顿法近似,这里为精确计算)。
- 线性化约束:\(g(x) \approx g(x_k) + \nabla g(x_k)^T (x - x_k) \le 0\),\(h(x) \approx h(x_k) + \nabla h(x_k)^T (x - x_k) = 0\)。
- 子问题形式:
\[ \min_{d} \quad \frac{1}{2} d^T B_k d + \nabla f(x_k)^T d \]
\[ \text{s.t.} \quad g(x_k) + \nabla g(x_k)^T d \le 0 \]
\[ h(x_k) + \nabla h(x_k)^T d = 0 \]
\[ \|d\|_{\infty} \le \Delta_k \]
这里 $d = x - x_k$ 是试探步,$\|d\|_{\infty} = \max(|d_1|, |d_2|)$ 是信赖域约束,采用无穷范数简化反射操作。
-
反射策略:
- 如果子问题解出的试探步 \(d_k\) 使得某个分量达到信赖域边界(即 \(|d_{k,i}| = \Delta_k\)),并且该方向是可行的下降方向,则允许在该边界进行“反射”,即在下一个子问题中考虑利用边界外的区域。实际操作中,通常是在求解子问题时直接处理带边界的约束,而“反射”体现在迭代更新中:如果试探步被边界截断,我们可以通过调整信赖域半径或接受截断步来继续。
-
计算实际下降量与预测下降量之比:
- 定义实际下降:\(\text{ared}_k = f(x_k) - f(x_k + d_k) + \mu_k \max(0, g(x_k))\),其中 \(\mu_k > 0\) 是障碍参数,用于处理不等式约束违反度。初始可设 \(\mu_0 = 1\),并按规则更新。
- 预测下降:\(\text{pred}_k = -(\frac{1}{2} d_k^T B_k d_k + \nabla f(x_k)^T d_k) + \mu_k \max(0, g(x_k)) - \mu_k \max(0, g(x_k) + \nabla g(x_k)^T d_k)\)。
- 比率:\(\rho_k = \frac{\text{ared}_k}{\text{pred}_k}\)。
-
更新迭代点与信赖域半径:
- 如果 \(\rho_k \ge \eta_1\),接受试探步:\(x_{k+1} = x_k + d_k\)。
- 否则,拒绝试探步:\(x_{k+1} = x_k\)。
- 根据 \(\rho_k\) 调整信赖域半径:
- 若 \(\rho_k < \eta_1\),则 \(\Delta_{k+1} = \gamma_1 \Delta_k\)(收缩)。
- 若 \(\rho_k \ge \eta_2\),则 \(\Delta_{k+1} = \gamma_2 \Delta_k\)(扩张)。
- 否则,\(\Delta_{k+1} = \Delta_k\)。
-
更新障碍参数:
- 如果不等式约束违反度很小,则减小 \(\mu_k\);否则保持或增大。简化起见,可设 \(\mu_k = \max(0.1 \mu_{k-1}, 10^{-6})\) 当约束违反度小于 \(0.1\) 时。
-
检查终止条件:
- 如果 \(\|d_k\| < \epsilon\) 且约束违反度小于 \(\epsilon\),则停止。
步骤二:第一次迭代计算(从 \(x_0 = (0,0)^T\) 开始)
1. 计算函数值与梯度
首先,计算目标函数和约束在 \(x_0\) 处的值及 Hessian:
- \(f(x) = (x_1 - 2)^4 + (x_1 - 2x_2)^2\)
- \(g(x) = x_1^2 - x_2\)
- \(h(x) = x_1 + x_2 - 2\)
梯度:
- \(\nabla f(x) = [4(x_1-2)^3 + 2(x_1-2x_2), \quad -4(x_1-2x_2)]^T\)
- \(\nabla g(x) = [2x_1, \quad -1]^T\)
- \(\nabla h(x) = [1, \quad 1]^T\)
在 \(x_0 = (0,0)\) 处:
- \(f_0 = (0-2)^4 + (0-0)^2 = 16\)
- \(g_0 = 0^2 - 0 = 0\)(等式约束违反度 0)
- \(h_0 = 0 + 0 - 2 = -2\)(等式约束违反度 2)
- \(\nabla f_0 = [4(-2)^3 + 2(0), \quad -4(0)]^T = [-32, \quad 0]^T\)
- \(\nabla g_0 = [0, \quad -1]^T\)
- \(\nabla h_0 = [1, \quad 1]^T\)
Hessian \(H_f(x)\)(用于精确 \(B_0\)):
\[\nabla^2 f(x) = \begin{bmatrix} 12(x_1-2)^2 + 2 & -4 \\ -4 & 8 \end{bmatrix} \]
在 \(x_0\) 处:
\[B_0 = \nabla^2 f(0,0) = \begin{bmatrix} 12 \cdot 4 + 2 & -4 \\ -4 & 8 \end{bmatrix} = \begin{bmatrix} 50 & -4 \\ -4 & 8 \end{bmatrix} \]
2. 构建二次规划子问题
子问题(在 \(x_0\) 处,\(d = (d_1, d_2)^T\)):
\[\min_{d} \quad \frac{1}{2} [d_1, d_2] \begin{bmatrix} 50 & -4 \\ -4 & 8 \end{bmatrix} \begin{bmatrix} d_1 \\ d_2 \end{bmatrix} + [-32, 0] \begin{bmatrix} d_1 \\ d_2 \end{bmatrix} \]
即
\[\min_{d} \quad 25 d_1^2 - 4 d_1 d_2 + 4 d_2^2 - 32 d_1 \]
约束:
\[g_0 + \nabla g_0^T d = 0 + [0, -1] d = -d_2 \le 0 \quad \Rightarrow \quad d_2 \ge 0 \]
\[ h_0 + \nabla h_0^T d = -2 + d_1 + d_2 = 0 \quad \Rightarrow \quad d_1 + d_2 = 2 \]
\[ \|d\|_{\infty} \le 0.5 \quad \Rightarrow \quad |d_1| \le 0.5, \; |d_2| \le 0.5 \]
3. 求解子问题
从等式约束 \(d_1 + d_2 = 2\) 和信赖域约束 \(|d_1|,|d_2| \le 0.5\) 可看出矛盾:因为 \(d_1 + d_2 = 2\) 意味着两者之和为 2,但每个的绝对值不超过 0.5,最大和最多为 1,所以无可行解。这时需要松弛处理:实际上,在算法中,当子问题因信赖域限制而无可行解时,我们通常先忽略等式约束,求解一个仅满足信赖域约束的梯度相关方向,或者放宽信赖域。但这里,由于是第一次迭代,我们可以先尝试忽略信赖域约束,求解仅含线性化约束的子问题,然后若解超出信赖域,再截断到信赖域内。
忽略信赖域约束,仅考虑 \(d_2 \ge 0\) 和 \(d_1 + d_2 = 2\)。从等式得 \(d_2 = 2 - d_1\),代入 \(d_2 \ge 0\) 得 \(2 - d_1 \ge 0 \Rightarrow d_1 \le 2\)。代入目标函数:
\[Q(d_1) = 25 d_1^2 - 4 d_1 (2-d_1) + 4 (2-d_1)^2 - 32 d_1 \]
计算:
- 展开:\(25 d_1^2 - 8 d_1 + 4 d_1^2 + 4(4 - 4d_1 + d_1^2) - 32 d_1\)
- 合并:\(25 d_1^2 + 4 d_1^2 + 4 d_1^2 = 33 d_1^2\);常数项 \(16\);一次项:\(-8 d_1 - 16 d_1 - 32 d_1 = -56 d_1\)
所以 \(Q(d_1) = 33 d_1^2 - 56 d_1 + 16\),导数为 \(66 d_1 - 56 = 0 \Rightarrow d_1 = 56/66 \approx 0.8485\),则 \(d_2 = 2 - 0.8485 = 1.1515\),不满足 \(d_2 \ge 0\) 吗?满足(1.1515 > 0)。但此时 \(|d_1|,|d_2|\) 都大于 0.5,违反信赖域约束。
因此,我们需要求解带信赖域约束的子问题。由于是二维,可以尝试边界搜索。信赖域约束是 \(|d_1| \le 0.5, |d_2| \le 0.5\),结合 \(d_1 + d_2 = 2\),显然不可能同时满足。所以子问题无可行解。此时算法处理方式通常是:先求解一个可行性恢复步骤,例如求解以下问题:
\[\min_{d} \quad \|h(x_k) + \nabla h(x_k)^T d\|^2 \quad \text{s.t.} \quad \|d\|_{\infty} \le \Delta_k \]
这里就是最小化 \((d_1 + d_2 - 2)^2\) 在 \(|d_1|,|d_2| \le 0.5\) 下。显然最优是让 \(d_1 + d_2\) 尽可能接近 2,但最大为 1(当 \(d_1=d_2=0.5\)),此时最小值为 \((1-2)^2=1\)。取 \(d_1 = d_2 = 0.5\),则 \(d = (0.5, 0.5)^T\),检查不等式约束:\(-d_2 = -0.5 \le 0\) 满足。所以可行性恢复步为 \(d_0 = (0.5, 0.5)^T\)。
4. 计算实际下降与预测下降
取障碍参数初始值 \(\mu_0 = 1\)。
当前约束违反度:\(\max(0, g_0) = 0\),\(|h_0| = 2\)(等式约束违反度用绝对值)。
实际函数值变化:\(x_1 = x_0 + d_0 = (0.5, 0.5)\)
计算 \(f_1 = (0.5-2)^4 + (0.5 - 1)^2 = (-1.5)^4 + (-0.5)^2 = 5.0625 + 0.25 = 5.3125\)
实际下降:\(\text{ared}_0 = f_0 - f_1 + \mu_0 (\max(0,g_0) - \max(0,g_1) + |h_0| - |h_1|)\)。注意这里我们简单将等式违反度纳入实际下降度量,更正式的是用增广拉格朗日,但这里简化处理。
计算 \(g_1 = 0.5^2 - 0.5 = 0.25 - 0.5 = -0.25 \le 0\),\(\max(0,g_1)=0\)。\(h_1 = 0.5+0.5-2 = -1\),\(|h_1|=1\)。
所以 \(\text{ared}_0 = 16 - 5.3125 + 1 \times (0 - 0 + 2 - 1) = 10.6875 + 1 = 11.6875\)。
预测下降:子问题中我们实际上求解了可行性恢复步,没有优化目标,但可以计算模型预测下降。原目标函数在 \(x_0\) 处的二次模型为 \(m_0(d) = f_0 + \nabla f_0^T d + \frac{1}{2} d^T B_0 d\)。
\(m_0(0) = 16\),\(m_0(d_0) = 16 + [-32,0] \cdot [0.5,0.5]^T + \frac{1}{2} [0.5,0.5] \begin{bmatrix} 50 & -4 \\ -4 & 8 \end{bmatrix} [0.5,0.5]^T\)
计算:\(\nabla f_0^T d_0 = -32 \times 0.5 = -16\)。
二次项:\(\frac{1}{2} [0.5,0.5] \begin{bmatrix} 50*0.5 -4*0.5 \\ -4*0.5+8*0.5 \end{bmatrix} = \frac{1}{2} [0.5,0.5] \begin{bmatrix} 25-2 \\ -2+4 \end{bmatrix} = \frac{1}{2} [0.5,0.5] \begin{bmatrix} 23 \\ 2 \end{bmatrix} = \frac{1}{2} (0.5*23 + 0.5*2) = \frac{1}{2} \times 12.5 = 6.25\)
所以 \(m_0(d_0) = 16 - 16 + 6.25 = 6.25\)。
预测下降:\(\text{pred}_0 = m_0(0) - m_0(d_0) = 16 - 6.25 = 9.75\)。
5. 计算比率并更新
\(\rho_0 = \frac{\text{ared}_0}{\text{pred}_0} = \frac{11.6875}{9.75} \approx 1.199\)
因为 \(\rho_0 \ge \eta_2 = 0.75\),所以接受试探步:\(x_1 = (0.5, 0.5)^T\),并扩张信赖域半径:\(\Delta_1 = \gamma_2 \Delta_0 = 2.0 \times 0.5 = 1.0\)。
障碍参数暂时保持不变:\(\mu_1 = 1\)。
步骤三:第二次迭代计算(从 \(x_1 = (0.5,0.5)^T\) 开始)
1. 计算函数值与梯度
在 \(x_1 = (0.5,0.5)\):
- \(f_1 = 5.3125\)(已算)
- \(g_1 = 0.25 - 0.5 = -0.25\)
- \(h_1 = 0.5+0.5-2 = -1\)
- \(\nabla f_1 = [4(0.5-2)^3 + 2(0.5-1), \quad -4(0.5-1)]^T = [4 \times (-1.5)^3 + 2 \times (-0.5), \quad -4 \times (-0.5)]^T = [4 \times (-3.375) - 1, \quad 2]^T = [-13.5 - 1, \quad 2]^T = [-14.5, \quad 2]^T\)
- \(\nabla g_1 = [2 \times 0.5, \quad -1]^T = [1, \quad -1]^T\)
- \(\nabla h_1 = [1, \quad 1]^T\)
- \(B_1 = \nabla^2 f(0.5,0.5) = \begin{bmatrix} 12(-1.5)^2 + 2 & -4 \\ -4 & 8 \end{bmatrix} = \begin{bmatrix} 12 \times 2.25 + 2 & -4 \\ -4 & 8 \end{bmatrix} = \begin{bmatrix} 29 & -4 \\ -4 & 8 \end{bmatrix}\)
2. 构建二次规划子问题
子问题:
\[\min_{d} \quad \frac{1}{2} d^T \begin{bmatrix} 29 & -4 \\ -4 & 8 \end{bmatrix} d + [-14.5, 2] d \]
即
\[\min_{d} \quad 14.5 d_1^2 - 4 d_1 d_2 + 4 d_2^2 - 14.5 d_1 + 2 d_2 \]
约束:
\[g_1 + \nabla g_1^T d = -0.25 + d_1 - d_2 \le 0 \quad \Rightarrow \quad d_1 - d_2 \le 0.25 \]
\[ h_1 + \nabla h_1^T d = -1 + d_1 + d_2 = 0 \quad \Rightarrow \quad d_1 + d_2 = 1 \]
\[ \|d\|_{\infty} \le \Delta_1 = 1.0 \]
3. 求解子问题
从等式约束得 \(d_2 = 1 - d_1\)。代入不等式:
\[d_1 - (1 - d_1) \le 0.25 \quad \Rightarrow \quad 2 d_1 - 1 \le 0.25 \quad \Rightarrow \quad 2 d_1 \le 1.25 \quad \Rightarrow \quad d_1 \le 0.625 \]
信赖域约束:\(|d_1| \le 1, |d_2| = |1 - d_1| \le 1\),自动满足因为 \(d_1 \in [0,1]\) 时 \(d_2 \in [0,1]\)。
代入目标函数:
\[Q(d_1) = 14.5 d_1^2 - 4 d_1 (1-d_1) + 4 (1-d_1)^2 - 14.5 d_1 + 2 (1-d_1) \]
展开:
- \(14.5 d_1^2\)
- \(-4 d_1 + 4 d_1^2\)
- \(+4(1 - 2d_1 + d_1^2) = 4 - 8 d_1 + 4 d_1^2\)
- \(-14.5 d_1\)
- \(+2 - 2 d_1\)
合并二次项:\(14.5 + 4 + 4 = 22.5 d_1^2\)
一次项:\(-4 d_1 - 8 d_1 - 14.5 d_1 - 2 d_1 = -28.5 d_1\)
常数项:\(4 + 2 = 6\)
所以 \(Q(d_1) = 22.5 d_1^2 - 28.5 d_1 + 6\),导数为 \(45 d_1 - 28.5 = 0 \Rightarrow d_1 = 28.5/45 = 0.6333\)
检查不等式 \(d_1 \le 0.625\),不满足。因此最优解在边界 \(d_1 = 0.625\) 处。此时 \(d_2 = 1 - 0.625 = 0.375\)。
验证信赖域:\(|d_1|=0.625 \le 1, |d_2|=0.375 \le 1\),可行。所以子问题解为 \(d_1 = (0.625, 0.375)^T\)。
4. 计算实际下降与预测下降
计算新点 \(x_2 = x_1 + d_1 = (1.125, 0.875)^T\)
计算函数值:
\(f_2 = (1.125-2)^4 + (1.125 - 2 \times 0.875)^2 = (-0.875)^4 + (1.125 - 1.75)^2 = 0.5862 + (-0.625)^2 = 0.5862 + 0.3906 = 0.9768\)
\(g_2 = 1.125^2 - 0.875 = 1.2656 - 0.875 = 0.3906 > 0\),违反不等式约束。
\(h_2 = 1.125 + 0.875 - 2 = 0\),满足等式。
实际下降:\(\text{ared}_1 = f_1 - f_2 + \mu_1 (\max(0,g_1) - \max(0,g_2) + |h_1| - |h_2|) = 5.3125 - 0.9768 + 1 \times (0 - 0.3906 + 1 - 0) = 4.3357 + 0.6094 = 4.9451\)
预测下降:计算模型值 \(m_1(0) = f_1 = 5.3125\)
\(m_1(d_1) = f_1 + \nabla f_1^T d_1 + \frac{1}{2} d_1^T B_1 d_1\)
\(\nabla f_1^T d_1 = [-14.5, 2] \cdot [0.625, 0.375]^T = -14.5 \times 0.625 + 2 \times 0.375 = -9.0625 + 0.75 = -8.3125\)
二次项:\(\frac{1}{2} [0.625,0.375] \begin{bmatrix} 29*0.625 -4*0.375 \\ -4*0.625+8*0.375 \end{bmatrix} = \frac{1}{2} [0.625,0.375] \begin{bmatrix} 18.125 - 1.5 \\ -2.5 + 3 \end{bmatrix} = \frac{1}{2} [0.625,0.375] \begin{bmatrix} 16.625 \\ 0.5 \end{bmatrix} = \frac{1}{2} (0.625*16.625 + 0.375*0.5) = \frac{1}{2} (10.3906 + 0.1875) = 5.2891\)
所以 \(m_1(d_1) = 5.3125 - 8.3125 + 5.2891 = 2.2891\)
预测下降:\(\text{pred}_1 = m_1(0) - m_1(d_1) = 5.3125 - 2.2891 = 3.0234\)
5. 计算比率并更新
\(\rho_1 = \frac{4.9451}{3.0234} \approx 1.635 > \eta_2 = 0.75\),接受试探步:\(x_2 = (1.125, 0.875)^T\),扩张信赖域:\(\Delta_2 = \gamma_2 \Delta_1 = 2.0 \times 1.0 = 2.0\)。
更新障碍参数:由于 \(g_2 > 0\),违反度 0.3906 较大,可保持 \(\mu_2 = 1\) 暂时不变。
步骤四:后续迭代与终止
继续迭代,直到步长很小且约束满足。这里只演示两次迭代,展示算法如何结合信赖域反射(通过信赖域约束控制步长)和序列二次规划(通过二次模型和线性化约束)。实际中,当接近最优解时,信赖域会稳定,子问题会给出接近牛顿步的方向,从而快速收敛。
最终答案:通过若干次迭代,该问题的最优解约为 \(x^* \approx (1.0, 1.0)^T\),此时 \(f^* = 0\),约束满足(\(g(x^*) = 0, h(x^*) = 0\))。算法利用了信赖域保证全局收敛,序列二次规划提供局部快速收敛,反射机制在边界附近调整搜索方向。