(完整word版)求解二次规划问题的拉格朗日及有效集方法

1、求解二次规划问题的拉格朗日及有效集方法最优化方法课程实验报告学 院：数学与统计学院班 级：硕2041班姓 名：王选学 号：3112054028指导教师：阮小娥同组人：钱东东求解二次规划问题的拉格朗日及有效集方法求解二次规划问题的拉格朗日 及有效集方法摘要二次规划师非线性优化中的一种特殊情形， 它的目标函数是二次实函数，约 束函数都是线性函数。由于二次规划比较简单，便于求解（仅次于线性规划），并且一些非线性优化问题可以转化为求解一些列的二次规划问题，因此二次规划的求解方法较早引起人们的重视，称为求解非线性优化的一个重要途径。 二次规 划的算法较多，本文仅介绍求解等式约束凸二尺规划的拉格朗日方法以

2、及求解一 股约束凸二次规划的有效集方法。关键字：二次规划，拉格朗日方法，有效集方法。-13 -【目录】摘要-1 -1等式约束凸二次规划的解法 -3 -1.1 问题描述-3 -1.2 拉格朗日方法求解等式约束二次规划问题 -3 -1.2.1 拉格朗日方法的推导 -3 -1.2.2 拉格朗日方法的应用 -4 -2 一般凸二次规划问题的解法 -5 -2.1 问题描述-5 -2.2 有效集法求解一般凸二次规划问题 -6 -2.2.1 有效集方法的理论推导 -6 -2.2.2 有效集方法的算法步骤 -9 -2.2.3 有效集方法的应用 -10 -3总结与体会-11 -4附录-11 -4.1 拉格朗日方法

3、的 matlab程序-11 -4.2 有效集方法的 Matlab程序 -11 -1等式约束凸二次规划的解法1.1问题描述我们考虑如下的二次规划问题E C 1 TTmin x Hx c x,/2(1.1)s.t. Ax=b其中 HwRn>n对称正定，AwRm"!1满秩，c,xWRn, bRm。1.2拉格朗日方法求解等式约束二次规划问题1.2.1 拉格朗日方法的推导首先写出拉格朗日函数：1 TT TL(x, £) = 一 x Hx +c x 九(Ax -b) , (1.2)2 xL(x2) = 0,九L(x,九)=0 ,得到方程组Hx -AT - -c, -Ax - -b

4、.-AT0=c .(1.3)IL-b将上述方程组写成分块矩阵形式:H一-乂我们称伤处方程组的系数矩阵-ATI0为拉格朗日矩阵。下面的定理给出了线性方程组(1.1)有唯一解的充分条件。定理1设HWRmM对称正定, 足二阶充分条件，即AWRmM行满秩。若在问题（1.1）的解x*处满dTHd 0, -dRn,d = 0,Ad - 0,则线性方程组（1.4）的系数矩阵非奇异，即方程组（1.4）有唯一解。其中，方程 组（1.4）为（1.1）对应的齐次方程组：-H -AT Vd八、1 = 0 (1.4).、A 0下面，我们来推导方程(1.3)的求解公式。根据定理 1,拉格朗日矩阵必然 是非奇异的，故可设其

5、逆为H -ATG -BT=I.jA 0 _ B C _由恒等式:H -ATG -BT- In Onxml- A 0 B C10mMIm .可得HG ATB =In , -HBT - ATC =0nm ,- AG =0mn , ABT = Im .于是由上述四个等式得到矩阵 G, B, C的表达式- -1a T T G = H -H A (AH A ) AH ,(1.5)1TliB =(AH A ) AH ,(1.6)_1 T 1C = -(AH A ) .(1.7)-BT -c 二-Gc BTbC lb b. |L Bc-Cb '因此，由(1.3)可得解得表达式(1.8)其中，G, B

6、,C 分别由(1.5),(1.6),(1.7) 给出下面给出X和工的另一种等价表达式。设Xk是问题(1.1)的任一可行点，即xk满足Axk=b。而在此点处目标函数的梯度为gk = $f (xk) = Hxk +c ,利用xk和gk,可将(1.8)改写为(1.9)1.2.2拉格朗日方法的应用(1)拉格朗日方法的Matlab程序见附录(2)利用拉格朗日方法求解下列问题：mins.t.2 c 22 cx12x2 x3 - 2x1x2 x3,解容易写出2x1 - x2 x3=2.一2-2-20102PI0L11A =2 -1利用Matlab程序求解该问题可以结果如下:1.909090909090909

7、1.9545454545454550. 1353636363636371am =2. 636363636363636-1.363636363636363fval =3,9772727272727282 一般凸二次规划问题的解法2.1问题描述考虑一般二次规划，S.t.min 1xT Hx cTx, 2aTx-bi =0,i E =1, ,l,(2.1)aiTx - bi _0,i I =l 1, ,m其中H是n阶对称阵。记I(x )=i |a：x -> =0” I,卜面的定理给出了问题(2.1)的一个最优性充要条件。定理2 x是二次规划问题(2.1)的局部极小点当且仅当x1 x2 x3 =

8、 4,(1)存在九* W Rm ,使得,一 *一一* 一一*_Hx +c-£ 九a -£ >a ai =0,i舌ieT * IC . L/ai x bi =0,i = E,. T *,ai x 一也 A0,i = I,* * *九i 之 0,i w I;九i = 0,i w I I (x )记S=d w Rn 0 |dTai =0,i w E;dTai 之 0,iw I(x*);dTai = 0,i w I(x *)且九* > 0.则对于任意的dwS,均有dTHd之0.容易发现，问题(2.1)是凸二次规划的充要条件是H半正定。此时，定理 2 的第二部分自然满足。

9、注意到凸优化问题的局部极小点也是全局极小点的性质， 我们有下面的定理：定理3 x是凸二次规划的全局极小点的充要条件是x满足KT条件，即存在九W Rm ,使得*Hx +c-Z iai -£ iai 0, i佳归T *一 ai x -8=0” E, T *一ai x - bi 之 Qi = I, _ *内至 0,i W I; % =0,i W I I(x ).2.2有效集法求解一般凸二次规划问题2.2.1 有效集方法的理论推导首先引入下面的定理，它是有效集方法理论基础。. *定理4设x是一般凸二次规划问题的全局极小点，且在 x处的有效集为S(x*) = E U I (x*),则x*也是下

10、列等式约束凸二次规划(2.2)1 Tu . T min - x Hx 十c x,2T.*st. ai x - b = 0,i S(x ).的全局极小点。从上述定理可以发现，有效集方法的最大难点是事先一般不知道有效集*、S(x),因此只有想办法构造一个集合序列去逼近它，即从初始点X0出发，计算有效集S(X0),解对应的等式约束子问题。重复这一做法，得到有效集序列*.S(Xk)，k=01，使之S(Xk)T S(X '以获得原问题的最优解。基于上述定理，我们分4步来介绍有效集方法的算法原理和实施步骤。第一步，形成子问题并求出搜索方向dk.设Xk是问题(2.1)的一个可行点，据此确定相应的有效

11、集 Sk =EUl(Xk),其中I(Xk)=i%-bi =0,iw I.求解相 应的子问题1 TT(2.3)min x Hx c x,T2一一s.t. ai x - bi = 0, i Sk.上述问题等价于1 TT(2.4)min qk(d) =gd Hd gkd,st. aiTd =0,i Sk.其中x =xk +d,gk =GXk +c.设求出问题(2.4)的全局极小点为dk,%是对应的拉格 朗日乘子。第二步，进行线性搜索确定步长因子 外.假设dk #0 ,分两种情形讨论。(1)若Xk +dk是问题(2.1)的可行点,即aT (Xk +dk) -bi =0,i E 及 aT(Xk +dk)

12、-bi 0,i I.则令'' k = 1, xk 1 = Xk dk .(2)若Xk +dk不是问题(2.1)的可行点，则通过线性搜索求出下降最好的可行 点。注意到目标函数是凸二次函数，那么这一点应该在可行域的边界上达到。因此只要求出满足可行条件的最大步长外即可。当i w Sk时，对于任意的二卜至0 ,者B有aTdk = 0和aT (Xk +c(kdk) =aTXk = bi , 此时，«k之0不受限制。当TSk时，即第i个约束是严格的不等式约束，此时要求外满足aT(Xk +Adk)之bi ,即aTdk -bi -aTXk,i Sk.注意到上式右端非正，故当aTdk之

13、0时，上式包成立。而当aTdk <0时，由上式可解得h -aTXka%故有bi - a xkaTdk| aTdk:二 0合并第(1)(2)可得:k = min1, ： k(2.5).第三步，修正Sk.当外=1时，有效集不变，即5小=5而当外<:1时,bik - aiTXka：dk故aT(xk +akdk) =bik ,因此在Xk+处增加了一个有效约束，即Sk i ： = Sk i k.第四步，考虑dk=0的情形。止匕时Xk是问题(2.3)的全局极小点。若这是对应 的不等式约束的拉格朗日乘子均为非负， 则Xk也是问题(2.1)的全局极小点，迭代 终止。否则，如果对应的不等式约束的拉格

14、朗日乘子有负的分量， 那么需要重新 寻找一个下降可行方向。设九八<0, jk -I(Xk),现在要求一个下降可行方向dk,满足gTdk<0且 aTdk =QW亡E;aTdk >0,Vj I(Xk),为简便计算，按下述方式选取 dk ：耳国九)bjk，aT (Xk dk) = bj, j Sk, j 二 jk,T ajkdk0,T(2.6)ajdk =0jSk, j - jk,另一方面，注意到Xk是子问题(2.3)的全局极小点，故有HXk c ? /"k ai =0,iSk即gk 二 Ak k,其中Ak - (ai )i Sk， ' k - ( ' i

15、 )i Sk .从而，g：dk =7；ATdk.由(2.6)知ATdk 八(aTdk)ej =(a；kdk)ejk,j Sk于是有gTdk = T(aTkdk)ej = 'kk (aLdk):二 0.上式表明，由(2.6)确定的dk是一个下降可行方向。因此，令SJ=Sk jk,则修 正后的子问题min qk(d) =1dTHd gTd,T 2°s.t. aTd =0,i Sk'的全局极小点必然是原问题的一个下降可行方向。2.2.2有效集方法的算法步骤经过上面的分析和推导，我们现在可写出有效集方法的算法步骤：(1)选取初值。给定初始可行点xo w Rn,令k:=0.(2

16、)解子问题。确定相应的有效集Sk = E U I(xk).求解子问题min qk(d) = 2dTHd +gTd,T _s.t. ai d =0,i Sk,得极小点dk和拉格朗日乘子向量 九k.若dk#0转步骤(4)；否则，转步骤(3)。(3)检验终止准则。计算拉格朗日乘子,k = Bkgk ,其中gk =也卜 c,Bk =(AkH,AT)AkH 入=(a)Sk.令( k)t = min( J.i, I(Xk)若(?，k)t >0,则Xk是全局极小点，停算。否则，若(%)t <0,则令Sk=Sk t, 转步骤(2)。(4)确定步长 ctk .令01k = min1,ctk,其中a

17、k = min Jbi _ ai xk i T . nT | ai d k< 0 >.ai dk令xk 1 ： = xk >.1kdk.(5)若% =1 ,则令 Sk-=Sk；否则，若外<1,则令S =SkUjk,其中jk满足bjk - aTkXkaLdk(6)令 k:=k+1,转步骤(1)2.2.3有效集方法的应用(1)有效集法的Matlab程序见附录。(2)用有效集方法求解下列二次规划问题: 22min f(x) =x1 -x1x2 +2x2 - x1 -10x2,«st.-3x1 -2x2 >-6x1 之 0, x2 之 0.解首先确定有关数据：H

18、l"c二 一1卜10Ae =, be=, Ai =-31-2 10 , bi =1一-610J .利用Matlab计算可得结果如下:ex it fl ag =O.5OOO2.2500lambda 二0.7500output -fval: -13. 7500iter: 33总结与体会通过本次实验，笔者对求解等式约束凸二次规划问题的拉格朗日方法及一般 约束条件下凸二次规划问题的有效集方法有了较深的认识和了解。感谢阮老师的悉心教诲和指导，使得笔者对最优化课程中的理论推导、算法步骤及编程都比较熟悉，对以后的科研工作有很好的指导和借鉴意义。4附录4.1拉格朗日方法的matlab程序(1)拉格朗

19、日算法函数%本程序用拉格朗日方法求解灯饰约束条件的二次规划问题。function x,lam,fval=qlag(H,A,b,c)%功能：用拉格朗日方法求解等式约束二次规划：% min f(x)=0.5*x'Hx+c'x , s.t. Ax=b%输入：H,c分别是目标函数的矩阵和向量，A,b分别是约束条件中的矩阵和向量%输出：(x,lam)是KT点，fval是最优值IH=inv(H);AHA=A*IH*A'IAHA=inv(AHA);AIH=A*IH;G=IH-AIH'*IAHA*AIH;B=IAHA*AIH;C=-IAHA;x=B'*b-G*c;lam

20、=B*c-C*b;fval=0.5*x'*H*x+c'*x;(2)拉格朗日算法求解等式约束凸二次规划问题主程序：H=2 -2 0;-2 4 0;0 0 2;c=0 0 1'A=1 1 1;2 -1 1;b=4 2'x,lam,fval=qlag(H,A,b,c)4.2有效集方法的Matlab程序(1)有效集方法函数%本程序主要适用于求解一般约束条件下的凸二次规划问题function x,lamk,exitflag,output=qpact(H,c,Ae,be,Ai,bi,x0)%功能：用有效集方法解一般约束二次规划问题%min f(x)=0.5*x'*H

21、*x+c'*x,% s.t. a'_i*x-b_i=0,(i=1,1)，%a'_i*x-b_i>=0,(i=l+1,，m)%输入：x0是初始点，H,c分别是目标函数二次矩阵和向量；%Ae= (a_1,.,a_l) ,be=(b_1,.,b_l)；%Ai=(a_l+1,a_m),bi=(b_l+1,b_m)'.%输出：x是最优解，lambda是对应的乘子向量；output是结构变量% 输出极小值f(x)，迭代次数k等信息，exitflag是算法终止类型%初始化epsilon=1.0e-9;err=1.0e-6;k=0;x=x0;n=length(x);kma

22、x=1.0e3;ne=length(be);ni=length(bi);lamk=zeros(ne+ni,1);index=ones(ni,1);for i=1:niif Ai(i,:)*x>bi(i)+epsilonindex(i)=0;endend%算法主程序while k<=kmax%求解子问题Aee二口；if ne>0Aee=Ae;endfor j=1:niif index(j)>0Aee=Aee;Ai(j,:);endendgk=H*x+c;m1,n1=size(Aee);dk,lamk尸qsubp(H,gk,Aee,zeros(m1,1);if norm(dk)<=erry=0.0;if length(lamk)>ney,jk=min(lamk(ne+1:length(lamk);endif y>=0exitflag=0;elseexitfl