序列二次规划算法

序列二次规划法求解一般线性优化问题:12min (x)h (x)0,i E {1,...,m }s.t.(x)0,i {1,...,m }i i f g I =∈=⎧⎨≥∈=⎩ (1.1) 基本思想：在每次迭代中通过求解一个二次规划子问题来确定一个下降方向，通过减少价值函数来获取当前迭代点的移动步长，重复这些步骤直到得到原问题的解。

1.1等式约束优化问题的Lagrange-Newton 法考虑等式约束优化问题min (x)s.t.h (x)0,E {1,...,m}j f j =∈=(1.2)其中:,n f R R →:()n i h R R i E →∈都为二阶连续可微的实函数. 记1()((),...,())T m h x h x h x =. 则(1.3)的Lagrange 函数为： 1(,)()*()()*()mT i i i L x u f x u h x f x u h x ==-=-∑（1.3）其中12(,,...,)T m u u u u =为拉格朗日乘子向量。

约束函数()h x 的Jacobi 矩阵为：1()()((),...,())T T m A x h x h x h x =∇=∇∇.对（1.3）求导数，可以得到下列方程组：(,)()A()*(,)0(,)()T x u L x u f x x u L x u L x u h x ∇⎡⎤⎡⎤∇-∇===⎢⎥⎢⎥∇-⎣⎦⎣⎦（1.4）现在考虑用牛顿法求解非线性方程（1.4）.(,)L x u ∇的Jacobi 矩阵为：(,)()(,)()0T W x u A x N x u A x ⎛⎫-= ⎪-⎝⎭（1.5）其中221(,)L(,)()*()mxx iii W x u x u f x u h x ==∇=∇-∇∑是拉格朗日函数L(,)x u 关于x 的Hessen 矩阵.(,)N x u 也称为K-T 矩阵。

对于给定的点(,)k k k z x u =，牛顿法的迭代格式为：1k k k z z z +=+∆. 其中k k (d ,v )k z ∆=是线性方程组k k k k (,)()(x )A(x )u *()0(x )k k k k T T k k d W x u A x f A x v h ⎛⎫-⎛⎫-∇+⎛⎫= ⎪ ⎪ ⎪ ⎪-⎝⎭⎝⎭⎝⎭（1.6）的解。

注意：只要k A(x )行满秩且(,)k k W x u 是正定的，那么（1.6）的系数矩阵非奇异，且方程组有唯一解。

引理1：已知矩阵,n nn m U RS R ⨯⨯∈∈，则对任意满足*0T S x =的非零向量x 都有0Tx Ux >的充要条件是存在常数*0σ>，使得对任意的*σσ≥都有*(U *S*S )0,0T T n x x x R σ+>∀≠∈.证明略。

鉴于方程组（1.6）的求解数值不稳定，故考虑将它转化成一个严格凸二次规划问题.转化的条件是(1.4)的解点*x 处的最优性二阶充分条件成立，即对满足*()*0TA x d =的任一向量0d ≠，成立***(,)*0Td W x u d >。

再由引理1知：当0τ>充分小时，1(*,*)(*)(*)2T W x u A x A x τ+正定。

考虑（1.6）中的(,)k k W x u 用一个正定矩阵来代替，记1(,)(,)()()2k T k k k k k B x u W x u A x A x τ=+则当**(,)(,)k k x u x u →时，矩阵**B(,)x u 正定。

（1.6）的第一个展开式为k (,)*d (x )*(x )(x )*T T k k k k k k k W x u A v f A u -=-∇+将上式变形为：k 11[(,)()()]*d (x )*[()](x )22k k T T k k k k k k k k W x u A x A x A v u A x d f ττ+-++=-∇ 令~1:()2k Tk k k k u v u A x d τ=++后得：~k (,)*d (x )*(x )T k k k k k B x u A u f -=-∇. 因此，（1.6）等价于k ~k (x )(,)()*()0(x )k k k k T k k d f B x u A x A x h u ⎛⎫⎛⎫∇-⎛⎫⎪=- ⎪ ⎪ ⎪ ⎪⎝⎭⎝⎭⎝⎭（1.7）进一步，可以把方程（1.7）转换成如下严格凸二次规划：k k k k k 1min (d)(x ,u )d f(x )2..(x )A(x )d 0TT k q d B ds t h =+∇+=（1.8）方程（1.7）和（1.8）具有同解的。

1.2一般形式的约束优化问题将1.1节中构造二次规划子问题求解等式约束优化问题的思想推广到一般形式的约束优化问题(1.5)。

在给定点k k (x ,u ,)k k z λ=后，将约束函数线性化，并对拉格朗日函数进行二次多项式近似，得到下列二次规划子问题：k k k k k k k 1min (x ,u )d f(x )2(x )(x )d 0,i E ..(x )(x )d 0,i T T Ti i Ti i d W dh h s t g g I+∇⎧+∇=∈⎪⎨+∇=∈⎪⎩ （1.9）其中212k k k k k k {1,...,m },I {1,...m },W W(x ,u ,)(x ,u ,)k xx E L λλ====∇，拉格朗日函数为1211(,)()*()*g ()m m i i i i i i L x u f x u h x x λ===--∑∑.于是，迭代点k x 的校正步k d 以及新的拉格朗日乘子估计量11,k k u λ++可以分别定义为问题的一个K-T 点*x 和相应的拉格朗日乘子向量*,*u λ。

定理1：给定约束优化问题（1.1）的最优解*d 和相应的拉格朗日乘子*,*0u λ≥.假定在*x 处，下面的条件成立：（1） 有效约束的Jacobi 矩阵(x*)S J 行满秩，其中(*)E(*)S x I x =;（2） 严格互补松弛条件成立，即******(x )0,0,(x )0,(x )0;i i i i i i g g g λλλ≥≥=+>（3） 二阶最优性充分条件成立，即对满足*()*0TA x d =的任一向量0d ≠，成立****(,,)*0T d W x u d λ>.那么若k k k (x ,u ,)λ充分靠近***(x ,u ,)λ，则二次规划问题(1.9)存在一个局部极小点*d ，使得其对应的有效约束指标集*()S d 与原问题在*x 处的有效指标集*()S x 是相同的。

注意：在构造二次规划子问题时，需要计算拉格朗日函数在迭代点k x 处的Hessen 矩阵，计算量过大。

为了克服这个缺陷，韩世平基于牛顿-拉格朗日法提出了一种利用对称正定矩阵k B 来代替拉格朗日矩阵k W 的序列二次规划法。

对于一般约束优化问题（1.1），在迭代点k k k k z (x ,u ,)λ=，构造下列形式的二次规划子问题：k k k k k k k 1min (x ,u )d f(x )2(x )(x )d 0,i E ..(x )(x )d 0,i T T Ti i Ti i d B dh h s t g g I+∇⎧+∇=∈⎪⎨+∇=∈⎪⎩ （1.10）并且用（1.10）的解k d 作为原问题变量x 在第k 次迭代过程中的搜索方向。

其中k d 有一个好的性质是它许多罚函数（价值函数）的下降方向。

例如，对于L1精确罚函数：1(x,)f(x)[|h (x)||[g (x)]_|]i i i Ei IP σσ∈∈=++∑∑其中0σ>为罚参数，g (x)]_max{0,g (x)}i i =-。

为了保证SQR 方法的全局收敛性，通常借助价值函数来确定搜索步长。

用来衡量一维搜索的好坏。

算法（一般约束优化问题的SQP 方法）Step 0: 给定初始点12000(,u ,)R R R ,mmn x λ∈⨯⨯对称正定矩阵0nB R ∈.计算0(x )ETA h =∇，00(x )I TA g =∇，000E I A A A ⎡⎤=⎢⎥⎢⎥⎣⎦.选择参数1(0,),(0,1),2ηρ∈∈容许误差120,1,εε令:0.k =Step 1: 求解子问题（1.10）得最优解k d .Step 2: 若k 11||d ||ε≤且k 1k 12||||||()_||h g ε+≤,stop,得到(1.1)的一个近似KT 点,,k k k x u λ（）.Step 3: 对于某种价值函数(,)x σΦ，选择罚参数k σ，使得k d 是该函数在k x 处的下降方向。

Step 4: Armijo 搜索. 令k m 是使下列不等式成立的最小非负整数m :m m 'k k k k k k k k (d ,)(,)(,;d ),x x x ρσσηρσΦ+-Φ≤Φ令km 1:,:.k k k k k a x x a d ρ+==+Step 5: 计算1111111(),(),E k ETI Tk k k k k I k A Ah x Ah x A A +++++++⎡⎤=∇=∇=⎢⎥⎢⎥⎣⎦以及最小二乘乘子1111111k T k k k k k u A A A f λ-++++++⎡⎤⎡⎤=∇⎢⎥⎣⎦⎣⎦Step 6: 校正矩阵k B 为1k B +.令11111,(,,)(,,)k k k k x k k k x k k k s a d y L x u L x u λλ+++++==∇-∇1T Tk k k k k kk k T T k k k k kB s s B z z B B s B s s z +=-+其中(1)B s k k k k k k z y θθ=+-参数k θ定义为k k k k ..........1...........,s y 0.2s 0.8s ,s y 0.2s s Tk k kT k T k k kk k k T T k k k k kB s B s B s B s s y θ⎧≥⎪=⎨<⎪-⎩ Step 7: 令:1,k k =+转1.注意：（step 1）利用K-T 条件，问题(1.10)等价于1k 2k k k (,,)B (A )(A )()0,(,,)()A 0,0,g()A 0,[g()A ]0.E T I Tk k k E k I I k k H d u d u f x H d u h x d x d x d λλλλλ=--+∇==+=≥+≥+= （1.11）第三式是2m 维互补问题，定义光滑函数(,,)a b a b εΦ=+ 其中0ε>为光滑参数.令21(,,)((,,),...,(,,))T m a b a b a b εεεΦ=ΦΦ，其中k (,,)[g ()(A )]Ii i i k i a b x d ελΦ=++其中(A )I k i 表示A Ik 的第i 行.记12(,,,)R R R Rmm n z d u ελ+=∈⨯⨯⨯，那么（1.11）问题等价于12(,,)(z):(,,,)0(,,)(,,)H d u H H d u H d u d ελελλελ⎡⎤⎢⎥⎢⎥===⎢⎥⎢⎥Φ⎣⎦则(z)H 的Jacobi 矩阵为'2110000()()(z)000(z)0(z)E TI T kk k E k I kB A A H A v D A D ⎡⎤⎢⎥--⎢⎥=⎢⎥⎢⎥⎣⎦其中21(,,)(,...,)T m v d v v εελ=∇Φ=,i v 由下式确定：i v =而221121(z)diag((z),...,(z)),(z)diag((z),...,(z))m m D a a D b b ==，其中(z),(z)i i a b 由下式确定：11i Ii a b ==-给定参数(0,1)γ∈，定义非负函数(z)||(z)||min{1,||(z)||}.H H βγ=（step 3）中选择价值函数111(,)()[||()||||()_||]x f x h x g x σσΦ=++可令max{||u ||,||||}k k τλ=，g (x)_max{0,g (x)}i i =- 任意选择一个0δ>，定义罚参数的修正规则为111111,(2),k k k k σστδστδστδ------⎧≥+⎪=⎨+<+⎪⎩ （step 6）因为BFGS 修正公式需要满足曲率条件，即0Tk k s y >,而k y 可能不满足这一条件，为此有必要对k y 进行修正。