序列二次规划算法

序列二次规划法求解一般线性优化问题:12min(x)h(x)0,iE{1,...,m}s.t.(x)0,i{1,...,m}iifgI(1.1)基本思想：在每次迭代中通过求解一个二次规划子问题来确定一个下降方向，通过减少价值函数来获取当前迭代点的移动步长，重复这些步骤直到得到原问题的解。1.1等式约束优化问题的Lagrange-Newton法考虑等式约束优化问题min(x)s.t.h(x)0,E{1,...,m}jfj(1.2)其中:,nfRR:()nihRRiE都为二阶连续可微的实函数.记1()((),...,())Tmhxhxhx.则(1.3)的Lagrange函数为：1(,)()*()()*()mTiiiLxufxuhxfxuhx（1.3）其中12(,,...,)Tmuuuu为拉格朗日乘子向量。约束函数()hx的Jacobi矩阵为：1()()((),...,())TTmAxhxhxhx.对（1.3）求导数，可以得到下列方程组：(,)()A()*(,)0(,)()TxuLxufxxuLxuLxuhx（1.4）现在考虑用牛顿法求解非线性方程（1.4）.(,)Lxu的Jacobi矩阵为：(,)()(,)()0TWxuAxNxuAx（1.5）其中221(,)L(,)()*()mxxiiiWxuxufxuhx是拉格朗日函数L(,)xu关于x的Hessen矩阵.(,)Nxu也称为K-T矩阵。对于给定的点(,)kkkzxu，牛顿法的迭代格式为：1kkkzzz.其中kk(d,v)kz是线性方程组kkkk(,)()(x)A(x)u*()0(x)kkkkTTkkdWxuAxfAxvh（1.6）的解。注意：只要kA(x)行满秩且(,)kkWxu是正定的，那么（1.6）的系数矩阵非奇异，且方程组有唯一解。引理1：已知矩阵,nnnmURSR，则对任意满足*0TSx的非零向量x都有0TxUx的充要条件是存在常数*0，使得对任意的*都有*(U*S*S)0,0TTnxxxR.证明略。鉴于方程组（1.6）的求解数值不稳定，故考虑将它转化成一个严格凸二次规划问题.转化的条件是(1.4)的解点*x处的最优性二阶充分条件成立，即对满足*()*0TAxd的任一向量0d，成立***(,)*0TdWxud。再由引理1知：当0充分小时，1(*,*)(*)(*)2TWxuAxAx正定。考虑（1.6）中的(,)kkWxu用一个正定矩阵来代替，记1(,)(,)()()2kTkkkkkBxuWxuAxAx则当**(,)(,)kkxuxu时，矩阵**B(,)xu正定。（1.6）的第一个展开式为k(,)*d(x)*(x)(x)*TTkkkkkkkWxuAvfAu将上式变形为：k11[(,)()()]*d(x)*[()](x)22kkTTkkkkkkkkWxuAxAxAvuAxdf令~1:()2kTkkkkuvuAxd后得：~k(,)*d(x)*(x)TkkkkkBxuAuf.因此，（1.6）等价于k~k(x)(,)()*()0(x)kkkkTkkdfBxuAxAxhu（1.7）进一步，可以把方程（1.7）转换成如下严格凸二次规划：kkkkk1min(d)(x,u)df(x)2..(x)A(x)d0TTkqdBdsth（1.8）方程（1.7）和（1.8）具有同解的。1.2一般形式的约束优化问题将1.1节中构造二次规划子问题求解等式约束优化问题的思想推广到一般形式的约束优化问题(1.5)。在给定点kk(x,u,)kkz后，将约束函数线性化，并对拉格朗日函数进行二次多项式近似，得到下列二次规划子问题：kkkkkkk1min(x,u)df(x)2(x)(x)d0,iE..(x)(x)d0,iTTTiiTiidWdhhstggI（1.9）其中212kkkkkk{1,...,m},I{1,...m},WW(x,u,)(x,u,)kxxEL，拉格朗日函数为1211(,)()*()*g()mmiiiiiiLxufxuhxx.于是，迭代点kx的校正步kd以及新的拉格朗日乘子估计量11,kku可以分别定义为问题的一个K-T点*x和相应的拉格朗日乘子向量*,*u。定理1：给定约束优化问题（1.1）的最优解*d和相应的拉格朗日乘子*,*0u.假定在*x处，下面的条件成立：（1）有效约束的Jacobi矩阵(x*)SJ行满秩，其中(*)E(*)SxIx;（2）严格互补松弛条件成立，即******(x)0,0,(x)0,(x)0;iiiiiiggg（3）二阶最优性充分条件成立，即对满足*()*0TAxd的任一向量0d，成立****(,,)*0TdWxud.那么若kkk(x,u,)充分靠近***(x,u,)，则二次规划问题(1.9)存在一个局部极小点*d，使得其对应的有效约束指标集*()Sd与原问题在*x处的有效指标集*()Sx是相同的。注意：在构造二次规划子问题时，需要计算拉格朗日函数在迭代点kx处的Hessen矩阵，计算量过大。为了克服这个缺陷，韩世平基于牛顿-拉格朗日法提出了一种利用对称正定矩阵kB来代替拉格朗日矩阵kW的序列二次规划法。对于一般约束优化问题（1.1），在迭代点kkkkz(x,u,)，构造下列形式的二次规划子问题：kkkkkkk1min(x,u)df(x)2(x)(x)d0,iE..(x)(x)d0,iTTTiiTiidBdhhstggI（1.10）并且用（1.10）的解kd作为原问题变量x在第k次迭代过程中的搜索方向。其中kd有一个好的性质是它许多罚函数（价值函数）的下降方向。例如，对于L1精确罚函数：1(x,)f(x)[|h(x)||[g(x)]_|]iiiEiIP其中0为罚参数，g(x)]_max{0,g(x)}ii。为了保证SQR方法的全局收敛性，通常借助价值函数来确定搜索步长。用来衡量一维搜索的好坏。算法（一般约束优化问题的SQP方法）Step0:给定初始点12000(,u,)RRR,mmnx对称正定矩阵0nBR.计算00(x)ETAh，00(x)ITAg，000EIAAA.选择参数1(0,),(0,1),2容许误差120,1,令:0.kStep1:求解子问题（1.10）得最优解kd.Step2:若k11||d||且k1k12||||||()_||hg,stop,得到(1.1)的一个近似KT点,,kkkxu（）.Step3:对于某种价值函数(,)x，选择罚参数k，使得kd是该函数在kx处的下降方向。Step4:Armijo搜索.令km是使下列不等式成立的最小非负整数m:mm'kkkkkkkk(d,)(,)(,;d),xxx令km1:,:.kkkkkaxxadStep5:计算1111111(),(),EkETITkkkkkIkAAhxAhxAA以及最小二乘乘子1111111kTkkkkkuAAAfStep6:校正矩阵kB为1kB.令11111,(,,)(,,)kkkkxkkkxkkksadyLxuLxu1TTkkkkkkkkTTkkkkkBssBzzBBsBssz其中(1)Bskkkkkkzy参数k定义为kkkk..........1...........,sy0.2s0.8s,sy0.2ssTkkkTkTkkkkkkTTkkkkkBsBsBsBssyStep7:令:1,kk转1.注意：（step1）利用K-T条件，问题(1.10)等价于1k2kkk(,,)B(A)(A)()0,(,,)()A0,0,g()A0,[g()A]0.ETITkkkEkIIkkHdudufxHduhxdxdxd（1.11）第三式是2m维互补问题，定义光滑函数222(,,)2ababab其中0为光滑参数.令21(,,)((,,),...,(,,))Tmababab，其中222kk(,,)[g()(A)][g()(A)]2IIiiikiiikiabxdxd其中(A)Iki表示AIk的第i行.记12(,,,)RRRRmmnzdu，那么（1.11）问题等价于12(,,)(z):(,,,)0(,,)(,,)HduHHduHdud则(z)H的Jacobi矩阵为'2110000()()(z)000(z)0(z)ETITkkkEkIkBAAHAvDAD其中21(,,)(,...,)Tmvdvv,iv由下式确定：222k2[g()(A)]2iIiikivxd而221121(z)diag((z),...,(z)),(z)diag((z),...,(z))mmDaaDbb，其中(z),(z)iiab由下式确定：222kk222k1,[g()(A)]2g()(A)1.[g()(A)]2iiIiikiIikiiIiikiaxdxdbxd给定参数(0,1)，定义非负函数(z)||(z)||min{1,||(z)||}.HH（step3）中选择价值函数111(,)()[||()||||()_||]xfxhxgx可令max{||u||,||||}kk，g(x)_max{0,g(x)}ii任意选择一个0，定义罚参数的修正规则为111111,(2),kkkk（step6）因为BFGS修正公式需要满足曲率条件，即0Tkksy,而ky可能不满足这一条件，为此有必要对ky进行修正。参考文献[1]马昌凤.最优化方法及其matlab程序设计.科学出版社.2010.