0 引言

稳定序列二次规划(sSQP)方法作为序列二次规划(SQP)方法的重要扩展,致力于考虑子问题模型、假设条件、技术构造、收敛性质(全局收敛或收敛速度)与数值效果等方面的研究.尤其是对相对弱条件下的收敛速度与数值效果.值得注意的是,对于退化约束优化问题,当对应原始解的拉格朗日乘子不唯一时,快速收敛性的实现难度极大.自1998年,Wright^[1]首次提出求解退化不等式约束优化的sSQP方法以来,以 Wright,Hager,Izmailov,Fernández,Solodov,Gill及Robinson等为代表的sSQP方法及理论研究发展迅速,取得系列成果.文献^[1]的方法在每次迭代时需求解一个目标函数二次的极大极小问题,可等价转化为原始对偶空间上的稳定二次规划(QP)问题。该作者在二阶充分最优性条件(SOSC)、Mangasarian-Fromovitz约束规格(MFCQ)和严格互补等条件下,证明了sSQP方法的局部二次收敛性.文献^[2]进一步改进文献^[1]的sSQP方法,在SOSC和MFCQ等较弱的条件下,证明算法具备局部二次收敛性.Hager^[3]对收敛性假设条件进行研究,仅在SOSC条件下获得文献^[1]中算法的局部收敛性,并提出用一对不等式约束来表示等式约束,从而可将方法推广到一般约束优化问题.相比较,传统SQP方法对乘子唯一性有着较苛刻的要求,只能依靠严格MFCQ条件来保证.特别地,在等式约束情形下,迫使线性无关约束规格(LICQ)成立^[4].文献^[5]将sSQP方法推广到求解含等式与不等式约束的变分问题,当初始点充分靠近原始对偶稳定点时,仅需SOSC假设条件,实现了算法的超线性收敛性.对于等式约束优化问题,文献^[6]在非临界乘子的条件下建立了超线性收敛的sSQP方法.文献^[7]基于线性方程组,在无需LICQ和严格互补等较强的假设条件下,提出一个二次收敛的牛顿型sSQP方法,并给出了该方法超线性收敛的重要定理.文献^[8]将拟牛顿型sSQP方法拓广到变分不等式问题,通过对Bregman距离极小化来更新矩阵信息,这样仅需要SOSC便可证明算法的超线性收敛性.

上述文献主要聚焦在最优解的局部范围内构造有效算法,并研究不同假设条件对其收敛速度的影响,建立了一批局部超线性收敛或二次收敛的sSQP方法.这些文献都侧重于局部收敛速度的分析,对算法全局优化策略(全局收敛性质)未做深入研究.而此问题正是实际应用和优化学者们希冀解决的问题,也是衡量最优化算法有效性的重要指标之一.近年来，尽管国内外不少学者对全局化sSQP方法潜心研究,但成果有限.本文主要介绍如下几类的全局sSQP算法:罚函数型sSQP方法^[9-12]、滤子型sSQP方法^[13]和IR型sSQP方法^[14],并对约束优化问题sSQP方法的进一步研究进行探索性思考.

1 函数型sSQP方法

当前,对sSQP方法全局化策略的研究仍是一项极其具有挑战性的工作,其中约束优化sSQP方法全局化时罚函数的选取极其困难,这对于最优性或下降性影响甚大.Gill等^[9]引入原始对偶广义增广拉格朗日(AL)函数,提出了一个求解等式约束加简单界约束优化问题的全局收敛sSQP方法.Izmailov等^[11]充分利用AL方法的鲁棒性,在没有任何积极约束识别策略下构造了一个sSQP方法,并证明了算法的全局收敛性和局部收敛速度.随后,他们又在文献^[12]中以原始对偶精确罚函数^[15]为效益函数,提出了一个求解等式约束优化问题的sSQP方法,并分析了算法的全局收敛性和超线性收敛速度.下面参考文献^[12],介绍一个具体的罚函数型sSQP算法.

考虑等式约束优化问题:

$\begin{align} & min~f\left( x \right) \\ & s.t.~h\left( x \right)=0 \\ \end{align}$

(1)

其中f(x):Rⁿ→R,h(x)=(h₁(x),…,h_l(x))^T:Rⁿ→R^l至少是二阶可微函数.令问题(1)的Lagrange函数为

$L\left( x,\lambda \right)=f\left( x \right)+\left\langle \lambda ,h\left( x \right) \right\rangle ,$

(2)

其中,符号<·,·> 表示向量内积运算.则问题(1)对应的一阶最优性条件为

$\frac{\partial L}{\partial x}\left( x,\lambda \right)=0,h\left( x \right)=0$

(3)

进而,定义函数Φ:Rⁿ×R^l→ Rⁿ×R^l,

$\Phi \left( x,\lambda \right)=\left( \frac{\partial L}{\partial x}\left( x,\lambda \right),h\left( x \right) \right)$

(4)

对一个给定的原始对偶迭代点(x,λ)∈Rⁿ×R^l及稳定参数σ>0,考虑如下原始对偶空间上的sSQP子问题,并产生sSQP方向(ξ,η)，

$\begin{array}{*{35}{l}} \underset{\left( \xi ,\eta \right)}{\mathop{\min }}\,\left\langle f\prime \left( x \right),\xi \right\rangle +\frac{1}{2}\left\langle \frac{{{\partial }^{2}}L}{\partial {{x}^{2}}}\left( x,\lambda \right)\xi ,\xi \right\rangle +\frac{\sigma }{2}{{\left\| \lambda +\eta \right\|}^{2}} \\ s.t.h\left( x \right)+h\prime \left( x \right)\xi -\sigma \eta =0, \\ \end{array}$

(5)

其中h′(x)=(▽h₁(x),…,▽h_l(x))^T.

基于上述sSQP方向(搜索方向),文献^[12]采用如下原始对偶效益 函数^[15]φ_c1,c2:Rⁿ×R^l→ R,

$\begin{align} & {{\varphi }_{{{c}_{1}},{{c}_{2}}}}\left( x,\lambda \right)=\text{ }L\left( x,\lambda \right)+\frac{{{c}_{1}}}{2}{{\left\| h\left( x \right) \right\|}^{2}}+ \\ & \frac{{{c}_{2}}}{2}{{\left\| \frac{\partial L}{\partial x}\left( x,\lambda \right) \right\|}^{2}} \\ \end{align}$

(6)

此处选取适当罚参数c₁>0,c₂>0,可以保证由子问题(5)产生的sSQP方向是φ_c1,c2(x,λ)的下降方向.

由文献^[16]和文献^[17]知,原始对偶效益函数(6)是一个精确罚函数.即如果c₂>0充分小,c₁>0充分大,则φ_c1,c2(x,λ)任意稳定点都是问题(1)的稳定点;反之,如果最优性系统(3)的原始对偶解(x,λ)满足LICQ和SOSC等条件,则它必是罚函数φ_c1,c2(x,λ)的严格局部极小解.

为便于分析和陈述算法的收敛性条件,下面复述二阶充分条件(SOSC)(文献^[12])以及临界Lagrangian乘子和非临界Lagrangian乘子^{[6, 18-19]}等定义.记M(x)是由问题(1)在稳定点x处对应的Lagrangian乘子构成的集合.记号ker A,im A分别表示矩阵A的零空间和列空间,即ker A={d:Ad=0},im A={d:d=Ay}.

定义1 λ称∈M(x)满足二阶充分条件(SOSC),如果

$\left\langle \frac{{{\partial }^{2}}L}{\partial {{x}^{2}}}\left( \bar{x},\bar{\lambda } \right)\xi ,\xi \right\rangle >0,\forall \xi \in ker~h\prime (\bar{x})\backslash \left\{ 0 \right\}.$

(7)

定义2 λ称 ∈M(x)是临界乘子,如存在ξ∈ker h′(x)\{0},使得 $\frac{{{\partial }^{2}}L}{\partial {{x}^{2}}}\left( \bar{x},\bar{\lambda } \right)\xi \in im{{\left( h\prime (\bar{x}) \right)}^{T}}$ ,否则,称是非临界乘子.

显然,如λ∈M(x)满足SOSC,则其是非临界乘子,故非临界乘子条件弱于SOSC.

算法1 (文献^[12]中sSQP算法)

初始步 选取参数σ>0,c₁>0,c₂>0,C₁>0,C₂>0,δ>0,ρ>0,q>1,θ∈[0,1].固定连续函数(除了0点外函数值几乎处处为正)ψ₁,ψ₂,ψ₃:R₊→ R₊.选取原始对偶初始点(x⁰,λ⁰)∈Rⁿ×R^l,置k:=0.

步骤1 由(4)式计算Φ_k$\triangleq $Φ(x^k,λ^k),若Φ_k=0,终止.

步骤2 令σ_k=min{σ,‖Φ_k‖},σ=σ_k,(x,λ)=(x^k,λ^k),求解子问题(5)的KKT系统得到稳定点d^k=(ξ^k,η^k).若该系统无解,则进入步骤4.

步骤3 (i)如果

$\left\langle {{\varphi }_{{{c}_{1}},{{c}_{2}}}}\prime ({{x}^{k}},{{\lambda }^{k}}),{{d}^{k}} \right\rangle \le -\rho {{\left\| {{d}^{k}} \right\|}^{q}},$

(8)

则进入步骤5.

(ii)如果

$\left\| h({{x}^{k}}) \right\|\ge {{\psi }_{1}}({{\sigma }_{k}})$

(9)

且

$\left\langle h({{x}^{k}}),h\prime ({{x}^{k}}){{\xi }^{k}} \right\rangle \le -{{\psi }_{2}}(\|h({{x}^{k}})\|),$

(10)

则令c₁=c_1,k+δ,其中c_1,k=c₁(ρ‖d^k‖^q;x^k,λ^k;ξ^k,η^k).如果c₁≤C₁,进入步骤5;否则,令c₁=C₁,进入步骤4.

(iii)如果

$\left\| \frac{\partial L}{\partial x}({{x}^{k}},{{\lambda }^{k}}) \right\|\ge {{\psi }_{3}}({{\sigma }_{k}}),$

(11)

则令c₂=c_2,k+δ,其中c _2,k=c₂(ρ‖d^k‖^q;x^k,λ^k;ξ^k,η^k,c₁). 如果c₂≤C₂,进入步骤5;否则,令c₂=C₂.

步骤4 选取(n+l)阶对称正定矩阵Q_k,计算方向d^k=-Q_kφ′_c1,c2(x^k,λ^k).

步骤5 计算α_k=θ^j,此处j为满足不等式

$\begin{align} & {{\varphi }_{{{c}_{1}},{{c}_{2}}}}\left( \left( {{x}^{k}},{{\lambda }^{k}} \right)+{{\theta }^{j}}{{d}^{k}} \right)\le {{\varphi }_{{{c}_{1}},{{c}_{2}}}}\left( {{x}^{k}},{{\lambda }^{k}} \right) \\ & +{{\theta }^{j}}\left\langle \varphi {{\prime }_{{{c}_{1}},{{c}_{2}}}}({{x}^{k}},{{\lambda }^{k}}),{{d}^{k}} \right\rangle \\ \end{align}$

(12)

的最小非负整数.

步骤6 令(x^k+1,λ^k+1)=(x^k,λ^k)+α_k d^k,k:=k+1,返回步骤1.

在算法1的迭代过程中,常量值C₁,C₂的设置对该算法的执行有着重要的防御作用.一般地,这两个数值应该取足够大,方能保证罚参数c₁,c₂的更新法则不受影响.为了增强算法设计的有效性,当子问题(5)KKT系统无解或步骤3的任一情形都不被执行时,算法1执行步骤4的拟牛顿步防御措施,以修正sSQP方向使之具有下降性,从而算法具有良好的适定性和全局收敛性质.当初始点充分靠近稳定点时,仅需要非临界Lagrangian乘子假设,即可证明算法1的超线性收敛.对退化测试问题,算法1表现出了良好的数值效果.

2 滤子型sSQP方法

滤子法是近20年提出的求解非线性约束优化的一种有效方法.该方法的提出是为了避免在实际问题中设置罚参数.早期滤子法的研究归功于Fletcher和Leyffer^[20],随后该方法因良好的数值效果而备受关注,近期代表成果见文献^[21-23]等.文献^[13]对稳定QP子问题进行修正,并结合双滤子技术^[21]提出了求解一般约束优化问题的滤子型sSQP方法.下面对文献^[13]中的算法进行分析.

考虑如下一般非线性规划问题:

$\begin{array}{*{35}{l}} min~f\left( x \right) \\ s.t.~{{h}_{\varepsilon }}\left( x \right)=0, \\ {{h}_{I}}\left( x \right)\le 0 \\ \end{array}$

(13)

其中,h_ε(x)=(h_i(x),i∈ε),h_I(x)=(h_i(x),i∈I),ε= {1,2,·s,l},I= {l+1,l+2,·s,m},f:Rⁿ→ R,h_i:Rⁿ→ R,i∈ε∪I为连续可微函数.定义问题(13)的Lagrangian函数L(x,μ)为

$L\left( x,\mu \right)=f\left( x \right)+\sum\limits_{i\in \varepsilon \cup I}^{{}}{{}}{{\mu }_{i}}~{{h}_{i}}\left( x \right).$

(14)

这里μ=(μ₁,μ₂,·s,μ_m)^T是Lagrangian乘子向量.

对于当前给定的原始对偶估计迭代点对(x^k,μ^L,k),求解如下修正的稳定QP子问题:

$\begin{array}{*{35}{l}} \begin{align} & \underset{\left( d,\Delta \lambda \right)}{\mathop{min}}\,\text{ }f\prime {{\left( {{x}^{k}} \right)}^{T}}d+\frac{1}{2}{{d}^{T}}{{B}^{k}}d+\frac{{{\sigma }^{L,k}}}{2}\left\| {{\mu }^{L,k}}+ \right. \\ & {{\left. \Delta \lambda \right\|}^{2}} \\ \end{align} \\ s.t.{{h}_{\varepsilon }}\left( {{x}^{k}} \right)+h{{\prime }_{\varepsilon }}{{\left( {{x}^{k}} \right)}^{T}}d-{{\sigma }^{L,k}}\Delta {{\lambda }_{\varepsilon }}=0, \\ {{h}_{I}}\left( {{x}^{k}} \right)+h{{\prime }_{I}}{{\left( {{x}^{k}} \right)}^{T}}d-{{\sigma }^{L,k}}\Delta {{\lambda }_{I}}\le 0. \\ \end{array}$

(15)

其中,σ^L,k是稳定参数,B^k是▽_xx²L(x^k,μ^L,k)的近似对称正定矩阵.当B^k=▽_xx²L(x^k,μ^k),μ^L,k=μ^k,σ^L,k=σ(x^k,μ^k)时,其中

$\sigma \left( {{x}^{k}},{{\mu }^{k}} \right)=\left\| \left( \begin{align} & {{\nabla }_{x}}L\left( {{x}^{k}},{{\mu }^{k}} \right) \\ & {{h}_{\varepsilon }}\left( {{x}^{k}} \right) \\ & \min \left\{ -{{h}_{I}}\left( {{x}^{k}} \right),\mu _{I}^{k} \right\} \\ \end{align} \right) \right\|,$

(16)

问题(15)是一个传统的稳定QP子问题模型. 源于算法良好的收敛性质,设计了两种乘子更新方式,其中μ^L,k保证全局收敛,μ^k则满足局部收敛的需求.类似地,考虑B^k和σ^L,k的更新准则.

基于问题(15)的解(d^k,Δλ^k),定义原始对偶搜索方向(d^k,Δμ^k)(此处Δμ^k=Δλ^k+μ^L,k-μ^k),此方向对算法收敛性理论分析有着重要的作用,其等价于求解如下相容的稳定QP子问题:

$\begin{array}{*{20}{l}} \begin{array}{l} \mathop {min}\limits_{\left( {d,\Delta \lambda } \right)} {\mkern 1mu} {\rm{ }}f\prime {\left( {{x^k}} \right)^T}d + \frac{1}{2}{d^T}{B^k}d + \frac{{{\sigma ^{L,k}}}}{2}\left\| {{\mu ^{L,k}} + } \right.\\ {\left. {\Delta \mu } \right\|^2} \end{array}\\ \begin{array}{l} s.t.{h_\varepsilon }\left( {{x^k}} \right) + h{\prime _\varepsilon }{\left( {{x^k}} \right)^T}d - {\sigma ^{L,k}}\left( {\mu _\varepsilon ^k + \Delta {\mu _\varepsilon } - } \right.\\ \left. {\mu _\varepsilon ^{L,k}} \right) = 0, \end{array}\\ {{h_{\rm I}}\left( {{x^k}} \right) + h{\prime _{\rm I}}{{\left( {{x^k}} \right)}^T}d - {\sigma ^{L,k}}\left( {\mu _\varepsilon ^k + \Delta {\mu _\varepsilon } - \mu _\varepsilon ^{L,k}} \right) \le 0.} \end{array}$

(17)

在设计非线性规划问题(13)的有效算法时,除需要对目标函数进行极小化,同时要不断降低约束违反度函数:

$\tilde{h}\left( x \right)=\sum\limits_{i=1}^{l}{{}}\left| {{h}_{i}}\left( x \right) \right|+\sum\limits_{i=l+1}^{m}{{}}max\{{{h}_{i}}\left( x \right),0\}.$

(18)

然而,搜索方向(d^k,Δμ^k)难以直接达到这样的要求.因此,引入如下辅助目标函数Φ(x,μ)和松弛约束违反度函数 p(x,μ):

$\Phi \left( x,\mu \right)=f\left( x \right)+\frac{{{\sigma }^{L,k}}}{2}\|\mu {{\|}^{2}},$

(19)

$\begin{align} & p\left( x,\mu \right)=\sum\limits_{i=1}^{l}{{}}\left| {{h}_{i}}\left( x \right)-{{\sigma }^{L,k}}\left( {{\mu }_{i}}-{{\mu }_{i}}^{L,k} \right) \right|+ \\ & \sum\limits_{i=l+1}^{m}{{}}max\{{{h}_{i}}\left( x \right)-{{\sigma }^{L,k}}({{\mu }_{i}}-{{\mu }_{i}}^{L,k}),0\}. \\ \end{align}$

(20)

进而引进双滤子技术(基于文献^[21],但有所改进),包括“全局滤子”和“局部滤子”.前者致力于算法收敛于KKT点或稳定点,后者以实现其局部收敛速度.对于第k次迭代,“全局滤子”和“局部滤子”分别记为F_g^k,F_l^k.

定义3^[13] 称点对( $\hat{x},\hat{\mu }$ )被“全局滤子”元(x^l,μ^l)或(p(x^l,μ^l),Φ(x^l,μ^l))接受,如果它满足如下两个不等式之一:

$p({{x}^{l}},{{\mu }^{l}})-p\left( \hat{x},\hat{\mu } \right)\ge {{\gamma }_{1}}~{{\left[ p\left( \hat{x},\hat{\mu } \right) \right]}^{2}},$

(21)

$\Phi ({{x}^{l}},{{\mu }^{l}})-\Phi \left( \hat{x},\hat{\mu } \right)\ge {{\gamma }_{2}}{{\left[ p\left( \hat{x},\hat{\mu } \right) \right]}^{2}},$

(22)

此处γ₁,γ₂∈(0,1).

定义4^[13] 称点对( $\hat{x},\hat{\mu }$ )为“全局滤子”元,如果它被“全局滤子”所有元素(p(x^l,μ^l),Φ(x^l,μ^l))接受.

定义5^[13] 称点对( $\hat{x},\hat{\mu }$ )被“局部滤子”元(x^l,μ^l)或(h(x^l),f(x^l))接受,如果它满足下面两个不等式之一:

$\begin{align} & \tilde{h}({{x}^{l}})-\tilde{h}(\hat{x})\ge max\{{{\gamma }_{1}}\tilde{h}{{(\hat{x})}^{2}},min\{{{\gamma }_{3}},\sigma ({{x}^{l}}, \\ & {{\mu }^{l}})\}\}, \\ \end{align}$

(23)

$\begin{align} & f({{x}^{l}})-f(\hat{x})\ge max\{{{\gamma }_{2}}{{(\hat{x})}^{2}},min\{{{\gamma }_{3}},\sigma ({{x}^{l}}, \\ & {{\mu }^{l}})\}\}, \\ \end{align}$

(24)

此处γ₃>0是一个常数.

定义6^[13] 称点对( $\hat{x},\hat{\mu }$ )为“局部滤子”元,如果它被“局部滤子”所有元素(h(x^l),f(x^l))接受.

对于问题(13)及当前迭代点对(x^k,μ^k),求解稳定QP子问题(17)得到搜素方向(d^k,Δμ^k),并沿着此方向和步长α,分别定义函数Φ的实际下降量和线性预测下降量:

$\begin{align} & \Delta {{\Phi }^{k}}\left( \alpha \right)=\Phi \left( {{x}^{k}},{{\mu }^{k}} \right)-\Phi \left( {{x}^{k}}+\alpha {{d}^{k}},{{\mu }^{k}}+\alpha \Delta {{\mu }^{k}} \right), \\ & \Delta l_{\Phi }^{k}\left( \alpha \right)=-\alpha \left[ {{f}^{'}}{{\left( {{x}^{k}} \right)}^{\text{T}}}{{d}^{k}}+{{\sigma }^{L,k}}{{\left( {{\mu }^{k}} \right)}^{\text{T}}}{{\mu }^{k}} \right]. \\ \end{align}$

令 ${\hat{x}}$ =x^k+αd^k, $\hat{\mu }~$ =μ^k+αΔμ^k,当试探步长α=1时,定义KKT误差下降条件:

$v\sigma \left( \hat{x},\hat{\mu } \right)\le \tau {{\sigma }^{L,k}},\tau \in \left( 0,\frac{1}{2} \right),v>0.~$

(25)

此外,文献^[13]结合“全局滤子”接受条件,进一步利用回溯技术选择了恰当的步长α^k,判断条件和步长下界.

开关条件:

$\Delta l_{\Phi }^{k}({{\alpha }^{k}})\ge \delta {{[p({{x}^{k}},{{\mu }^{k}})]}^{2}},\delta \in \left( 0,\infty \right).~$

(26)

充分下降条件:

$\Delta {{\Phi }^{k}}({{\alpha }^{k}})\text{ }\ge \tau \Delta l_{\Phi }^{k}({{\alpha }^{k}}),~$

(27)

$\left\{ \begin{array}{l} min\{ {\gamma _1},\frac{{ - min\{ {\gamma _2}p({x^k},{\mu ^k}),\delta {{[p({x^k},{\mu ^k})]}^2}\} }}{{f\prime {{({x^k})}^T}{d^k} + {\sigma ^{L,k}}{{({\mu ^k})}^T}\Delta {\mu ^k}\} }},\\ 若f\prime {({x^k})^T}{d^k} + {\sigma ^{L,k}}{({\mu ^k})^T}\Delta {\mu ^k} < 0;\\ \alpha ,{\rm{否则}}. \end{array} \right.$

(28)

其中,常数α∈(0,1).

对于原始对偶迭代点对( $\hat{x},\hat{\mu }$ ),如果它被局部滤子接受,或者满足KKT误差下降条件(25),或者满足如下近似松弛一阶条件:

$\psi \left( \hat{x},\hat{\mu } \right)\le \in _{\sigma }^{k}$

(29)

其中,

$\begin{array}{l} \psi \left( {\hat x,\hat \mu } \right): = \\ \left\| {\left( \begin{array}{l} {\nabla _x}L\left( {\hat x,\hat \mu } \right)\\ {h_\varepsilon }\left( {\hat x} \right) - {\sigma ^{L,k}}({{\hat \mu }_\varepsilon } - \mu _\varepsilon ^{L,k})\\ \min \left( { - {h_I}\left( {\hat x} \right) + {\sigma ^{L,k}}\left( {{{\hat \mu }_\varepsilon } - \mu _I^{L,k},{{\hat \mu }_I}} \right)} \right) \end{array} \right)} \right\|, \end{array}$

(30)

则序列{∈_σ^k}下降趋于0.构造Lagrange乘子μ^L,k和稳定参数σ^L,k的更新准则:

$\left\{ \begin{align} & {{\mu }^{L,k+1}}:=\hat{\mu },若{{\left\| {\hat{\mu }} \right\|}_{\infty }}~\le {{\mu }_{max}}, \\ & {{\mu }^{L,k}},否则. \\ \end{align} \right.$

(31)

${{\sigma }^{L,k+1}}=min\left\{ \beta {{\sigma }^{L,k}},\upsilon \sigma \left( \hat{x},\hat{\mu } \right) \right\},$

(32)

其中,μ_max>0,β∈(0,1)是两个常数.

下面给出求解问题(13)的滤子型sSQP算法.

算法2 (文献^[13]中sSQP算法之外循环)

初始步 选取初始点x⁰∈Rⁿ,参数∈_σ⁰>0,u_l>0,u_g>0,δ>0,γ₃>0,μ_max>0, α∈(0,1),γ₁∈(0,1),γ₂∈(0,1),r∈(0,1),τ∈ (0, $\frac{1}{2}$ ),β∈(0,1),μ^L,0=0,μ⁰=0,υ>0,σ^L,0>0,F_l⁰={(u_l,-∞)},F_g⁰={(u_g,-∞)},Flag=GLOBAL,置k:=0.

步骤1  如果(x^k,μ^k)是问题(13)的KKT点或不可行稳定点,则终止.

步骤2 求解稳定QP子问题(17),产生原始对偶解(d^k,Δμ^k),令α=1, ${\hat{x}}$ =x^k+αd^k, ${\hat{\mu }}$ =μ^k+αΔμ^k.

步骤3 如果( $\hat{x},\hat{\mu }$ )被局部滤子接受或满足(25)式,则由准则(31)和(32)更新乘子μ^L,k+1和稳定参数σ^L,k+1,置Flag=LOCAL,F_g^k+1={(u_g,-∞)}.否则,进入步骤4.

步骤4 启动算法3(算法2之内循环),获得点对( $\hat{x},\hat{\mu }$ )和步长α.如果α<α_min^k,则令μ^k+1=μ^L,k,进入可行性恢复阶段并找到新的迭代点x^k+1,使得(x^k+1,μ^k+1)被全局滤子F_g^k接受.μ^L,k+1,σ^L,k+1,∈_σ^k+1,B^k+1保持不变.如果p(x^k,μ^k)>0,则令k:=k+1,将(p(x^k,μ^k),Φ(x^k,μ^k))加入F_g^k.

如果点对( $\hat{x},\hat{\mu }$ )满足(29)式,则由准则(31)和(32)更新μ^L,k+1和σ^L,k+1,令F_g^k+1={(u_g,-∞)},∈_σ^k+1=∈_σ^k/10.否则,令μ^L,k+1=μ^L,k,σ^L,k+1=σ^L,k.

步骤5 如果($\hat x,\hat \mu $)被局部滤子接受,则将($\tilde h\left( {\hat x} \right),f\left( {\hat x} \right)$)加入F_l^k,更新局部滤子.否则,进入步骤6.

步骤6 令x^k+1=,μ^k+1=,α^k=α,k:=k+1,且Flag=GLOBAL,利用拟牛顿公式更新矩阵B^k+1,返回步骤1.

算法3(算法2之内循环)

步骤1 令=x^k+α d^k,=μ^k+αΔμ^k.如果α≥α_min^k,则进入步骤2;否则,终止.

步骤2 如果($\hat x,\hat \mu $)被全局滤子接受,则进入步骤3.否则令α=r α,返回步骤1.

步骤3 如果开关条件(26)不成立或充分下降条件(27)成立,则终止.否则令α=r α,返回步骤1.

利用双滤子技术,由算法2产生的所有迭代点均能被滤子接受(全局滤子或局部滤子),这对整个算法的收敛性质分析起到至关重要的作用.在理论上,不需要任何约束规格,即可证明:滤子型sSQP算法2不仅具有全局收敛性(存在收敛子列或收敛于KKT点,或收敛于不可行稳定点),而且在SOSC下,算法达到超线性收敛速度.

3 IR型sSQP方法

sSQP方法的研究最早可以追溯到20世纪末,虽然在适当的假设条件或技术下早期的研究成果实现了算法快速的收敛速度,但在理论上能实现全局收敛的成果较少.滤子型sSQP方法需储存滤子信息,会造成存储量大,而且当迭代点远离可行域,或者在迭代过程中当试探步长比既定下界要小时,为寻找一个能被全局滤子接受的新迭代点,需要执行可行性恢复阶段,这无疑增加算法计算成本,从而影响数值效果.然而,IR技术^[24-27]可以减少可行性恢复阶段的复杂性,对sSQP算法的理论分析有着重要的作用.文献^[14]提出了求解含等式且变量有界约束优化的IR型sSQP方法.下面详细介绍此算法.

考虑问题(1),且满足约束条件x∈Ω={x∈Rⁿ|a≤x≤b},此处a,b∈Rⁿ.定义问题(1)的自然残差σ:Rⁿ×R^l→R，

$\sigma \left( x,\lambda \right)=\text{ }\left\| \left[ \begin{align} & {{P}_{\Omega }}\left( \text{ }x-\frac{\partial L}{\partial x}\left( x,\lambda \right) \right)-x \\ & h\left( x \right) \\ \end{align} \right] \right\|,$

(33)

其中,P_Ω表示在Ω上的正交投影,L:Rⁿ×R^l→R是Lagrangian函数(2).

对于原始对偶迭代点对(x^k,λ^k),罚参数ρ_k>0及Q_k，其中Q_k是问题(1)Lagrangian Hesse的近似,考虑如下拟牛顿型sSQP子问题:

$\begin{align} & \underset{\left( \xi ,\lambda \right)}{\mathop{min}}\,\left\langle f\prime \left( {{x}^{k}} \right),\xi -{{x}^{k}} \right\rangle +\frac{1}{2}{{Q}_{k}}\left( \xi -{{x}^{k}} \right),\xi \\ & -{{x}^{k}}\text{ }+\frac{1}{2{{\rho }_{k}}}{{\left\| \lambda \right\|}^{2}} \\ & s.t.h\left( {{x}^{k}} \right)+h\prime {{\left( {{x}^{k}} \right)}^{T}}\left( \xi -{{x}^{k}} \right)-\frac{1}{{{\rho }_{k}}}\left( \lambda -{{\lambda }^{k}} \right)=0, \\ & \xi \in \Omega . \\ \end{align}$

(34)

为构造全局收敛算法,引入下面的辅助函数F_k(x,λ),H_k(x,λ):Rⁿ × R^l→ R:

${{F}_{k}}\left( x,\lambda \right)=f\left( x \right)+\frac{1}{2{{\rho }_{k}}}{{\left\| \lambda \right\|}^{2}},$

(35)

${{H}_{k}}\left( x,\lambda \right)=h\left( x \right)-\frac{1}{{{\rho }_{k}}}(\lambda -{{\lambda }^{k}}).$

(36)

选取点Y^k(x)=(x,λ^k+ρ_kh(x)),易知H_k(Y^k(x))=0,∀x∈Rⁿ,且对于Y^k:=Y^k(x^k),Y:=(x,λ),sSQP子问题(34)等价于如下QP子问题:

$\begin{align} & \underset{Y}{\mathop{\min }}\,\left\langle {{F}^{\prime }}_{k}\left( {{Y}^{k}} \right),Y-{{Y}^{k}} \right\rangle +\frac{1}{2}\left\langle \left[ \begin{matrix} {{Q}_{k}} & 0 \\ 0 & \frac{1}{{{\rho }_{k}}}I \\ \end{matrix} \right](Y \right. \\ & \left. -{{Y}^{k}}),Y-{{Y}^{k}} \right\rangle \\ & s.t.{{H}^{\prime }}_{k}{{\left( {{Y}^{k}} \right)}^{T}}\left( Y-{{Y}^{k}} \right)=0,Y\in \Omega \times {{R}^{l}}. \\ \end{align}$

(37)

注意到H_k(Y^k)=0.不难发现上述QP子问题(37)对应于优化问题:

$\begin{array}{*{35}{l}} \underset{Y}{\mathop{\min }}\,{{F}_{k}}\left( Y \right) \\ s.t.{{H}_{k}}\left( Y \right)=0,Y\in \Omega \times {{R}^{l}} \\ \end{array}$

(38)

的QP近似子问题.通过近似求解问题(37)可获得sSQP方法的局部收敛性质^{[3, 5]}.

对每次迭代,IR方法包括恢复阶段和极小化阶段.在恢复阶段,给定迭代点X^k,计算恢复点Y^k,避免了目标函数值以及可行性恶化.在极小化阶段,基于可行性和最优性构造了一个含罚参数的效益函数,并沿着一阶可行方向进行线搜索.

下面给出IR型sSQP算法的具体步骤.

算法4(文献^[14]中sSQP算法)

初始步 选取参数γ∈(0,1),r∈(0,1),ε>0,序列{∈_k}满足∈_k $\searrow $ 0,α_L,α_U>0.选取任意近似初始点X⁰=(x⁰,λ⁰)∈Ω×R^l,ρ₀>0,ψ_-1=σ(X⁰),置k:=0.对当前参数ρ_k,定义

$\begin{align} & {{\Pi }_{k}}\left( x,\lambda \right)=(x,max\{-{{\alpha }_{L}}\sqrt{{{\rho }_{k}}}e,min\{\lambda , \\ & {{\alpha }_{U}}\sqrt{{{\rho }_{k}}}e\}\}). \\ \end{align}$

(39)

步骤1 计算σ(X^k),若σ(X^k)≤,则终止.

步骤2 (i)令X^k,0=(x^k,0,λ^k,0)=(x^k,λ^k),θ₀∈(0,1),Q_k,0为对称正定矩阵,j:=0.

(ii)计算Y^k,j=(x^k,j,λ^k,j+ρ_k h(x^k,j)),求解如下QP子问题,得到最优解D^k,j∈ Rⁿ×R^l.

$\begin{align} & \underset{D}{\mathop{min}}\,\left\langle {{F}^{\prime }}_{k}\left( {{Y}^{k,j}} \right),D \right\rangle +\frac{1}{2}\left[ \begin{matrix} {{Q}_{k}} & 0 \\ 0 & \frac{1}{{{\rho }_{k}}}I \\ \end{matrix} \right]D,D \\ & s.t.{{H}^{\prime }}_{k}{{\left( {{Y}^{k,j}} \right)}^{T}}D=0,{{Y}^{k,j}}=D\in \Omega \times \text{ }{{R}^{l}}. \\ \end{align}$

(40)

如果‖D^k,j‖∈_k,则令X^k+1=Π_k(Y^k,j+D^k,j),进入步骤3.

(iii) 令θ^j+1∈{2^-i:i∈N∪{0}}为满足下面不等式的最大值:

$\begin{align} & {{\Phi }_{k}}\left( {{Y}^{k,j}},{{\theta }^{j+1}} \right)-{{\Phi }_{k}}\left( {{X}^{k,j}},{{\theta }^{j+1}} \right)\le \\ & \frac{\bar{r}-1}{2}\|{{H}_{k}}({{X}^{k,j}})\|, \\ \end{align}$

(41)

此处Φ_k(X,θ)=θF_k(X)+(1-θ)‖H_k(X)‖是效益函数.

(iv) 令t_j∈{2^-i:i∈N∪{0}}为满足下面不等式的最大值:

$\begin{align} & {{\Phi }_{k}}\left( {{Y}^{k,j}}+{{t}_{j}}{{D}^{k,j}},{{\theta }^{j+1}} \right)-{{\Phi }_{k}}\left( {{X}^{k,j}},{{\theta }^{j+1}} \right)\le \\ & \frac{\bar{r}-1}{2}\|{{H}_{k}}({{X}^{k,j}})\|. \\ \end{align}$

(42)

(v) 令X^k,j+1=Y^k,j+t_jD^k,j,选取对称正定矩阵Q_k,j+1,令j:=j+1,返回步骤2(ii).

步骤3 更新罚参数,令ψ_k=min{ψ_k-1,σ(X^k)},如果‖h(x^k+1‖>γ‖h(x^k‖)且σ(X^k+1)>γψ_k,则令ρ_k+1=10ρ_k.否则令ρ_k+1=ρ_k.令k:=k+1,返回步骤1.

事实上，基于IR方法的思想,算法4无需实质性地修正sSQP子问题(34)的结构，只需求解子问题(37)得到非精确解，然后巧妙利用类似增广拉格朗日(AL)罚参数更新策略(步骤3),即可获得算法的全局收敛性质.在严格MFCO和SOSC等条件下,可保证罚参数序列是有界的,进而算法继承了良好的局部收敛速度(线性收敛).

4 展望

本文对求解约束优化问题的若干sSQP方法的思想及算法作了一个较详细的概述,主要介绍了早期sSQP方法的局部收敛成果;近年来实现了全局收敛的罚函数型sSQP方法、滤子型sSQP方法和非精确恢复(IR)型sSQP方法.对sSQP方法还有如下一些问题值得思考和研究.

① 如何修正原始对偶罚函数型sSQP算法,在适当的假设条件下,将方法拓展到一般约束优化问题,并制成软件包使之得到广泛应用;

② 现有的滤子型sSQP算法考虑内点法求解一个与传统QP子问题相当的修正稳定QP子问题,无疑会在严格内点的选择和计算上遇到困难;加之滤子技术需要进入可行性恢复阶段,也会增加计算成本.因此,如何改进算法以减少计算量,有待进一步探讨;

③ 目前基于IR技术构造sSQP方法时,需要严格MFCQ和SOSC等条件以实现算法线性收敛,没有继承传统sSQP方法优点(具有快速收敛速度的鲁棒性),如何修正罚参数以更新技术等,建立超线性收敛的IR型sSQP方法值得进一步深入研究;

④ 求解大规模稀梳优化与工程优化 (如庞杂的电气工程中的潮流问题、机组合问题)的高效的sSQP方法有必要进一步探索.