求解二次规划问题的拉格朗日及有效集方法.docx

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

2、录】摘要等式约束凸二次规划的解法1.1问题描述1.2拉格朗日方法求解等式约束二次规划问题1.2.1拉格朗日方法的推导1.2.2拉格朗日方法的应用一般凸二次规划问题的解法2.1问题描述2.2有效集法求解一般凸二次规划问题2.2.1有效集方法的理论推导2.2.2有效集方法的算法步骤-11 -11 -4.1拉格朗日方法的matlab程序4.2有效集方法的Matlab程序1等式约束凸二次规划的解法1.1问题描述我们考虑如下的二次规划问题 . 1 T 口 Tmin -x Hx + c x,z 八 2 (1.1) i s.t. Ax = b其中对称正定， R行满秩，c,x亡Rn，b亡Rm1.2拉格朗日方法

3、求解等式约束二次规划问题1.2.1拉格朗日方法的推导首先写出拉格朗日函数:1L(x, A) = xtHx+cTx A (Ax-b)，(1.2)2jL(x,k)=0, J丄(x,Q =0，得到方程组Hx - A 冬=-c,-Ax = -b.将上述方程组写成分块矩阵形式:H -A我们称伤处方程组的系数矩阵-at-0足二阶充分条件，即dTHd aOM 5,d 工 O,Ad =0,则线性方程组(1.4)的系数矩阵非奇异，即方程组(1.4 )有唯一解。其中，方程 组(1.4)为(1.1 )对应的齐次方程组：下面，我们来推导方程 是非奇异的，故可设其逆为由恒等式可得LHa(1.3)叮/0 (14)-的求解

4、公式。根据定理 1,拉格朗日矩阵必然-At0Tg -bt/i-B CH-atG-Bt r 1 n 0nxm-A0BC|_0mxi 1 m ”HG + ATB = ln AG =0m 河HBt -AtC =0恤ABt =lm .于是由上述四个等式得到矩阵 G, B, C的表达式1 1 T 1 T 1 1G =H-H A (AH A ) AH4 T 4 A.B =(AH A ) AH ,C = -(AH AT)因此，由(1.3)可得解得表达式ELF -叮0B -Gc + BTb C j-bBc-Cb其中，G, B,C 分别由(1.5),(1.6),(1.7) 给出。(1.5)(16)(1.7)(1

5、.8)F面给出X和几的另一种等价表达式。设xk是问题(1.1)的任一可行点，即Xk满足AXk =b。而在此点处目标函数的梯度为gk =Nf(Xk) = HXk 中c，利用 Xk和gk，可将(1.8)改写为x_ xk - Ggk Lji Bgk .(1.9)1.2.2拉格朗日方法的应用(1)拉格朗日方法的Matlab程序见附录。(2)利用拉格朗日方法求解下列问题：mins.t.2 2 2Xr + 2x2 +X3-2x1X2 +X3,xi +X2 + X3 = 4,2X1 X2 + X3=2.解容易写出f2-20-f0-H =-240c =00021aJ1-1b邛L2J利用Matlab程序求解该问

6、题可以结果如下:k 9090909090909091.9545454545454550.135*3637Ian =2. 630303636363636-L 363636363636363fval =3, 97；272727272：2S2 一般凸二次规划问题的解法2.1问题描述考虑一般二次规划min 1xT HcTx,2ajxbi =0,i- E=1,，l,s.t.! ajx-bi O,iJ =l +1,m(2.1)其中H是n阶对称阵。记l(x*)=i laTx*-bi =0,i迂1，F面的定理给出了问题(2.1)的一个最优性充要条件。定理2 x是二次规划问题(2.1)的局部极小点当且仅当(1)

7、存在Z e Rm，使得Hx +c2；入iNj 送 Zia0,任 igT *ai X bi = 0,i 亡 E,T *ai X -bi 0,i I,入* 0,i M* =0,飞 I I(x*)* , , *= 0,i |(x )且 M 0.记S =d Rn 0 |dTa0, E;dTa0, I(x*);dTai 则对于任意的d亡S，均有dTHd 0.容易发现，问题(2.1)是凸二次规划的充要条件是H半正定。此时，定理 2 的第二部分自然满足。注意到凸优化问题的局部极小点也是全局极小点的性质， 我们有下面的定理：定理3 X*是凸二次规划的全局极小点的充要条件是 X*满足KT条件，即存在力迂Rm，使

8、得Hx + C无几 idi 2 入 i ai = 0,任 iSaT X -bj =0,i 忘 E,aTx -bj 0,i 引，人 0,i 引；k* =0,i 引 I(x*).2.2有效集法求解一般凸二次规划问题2.2.1有效集方法的理论推导首先引入下面的定理，它是有效集方法理论基础。定理4设X是一般凸二次规划问题的全局极小点，且在 X处的有效集为S*) = eU|(x*)，则x*也是下列等式约束凸二次规划I . 1 T . TI min 2X HC X (2.2)st. a：X b = 0,i 亡 S(x ).的全局极小点。从上述定理可以发现，有效集方法的最大难点是事先一般不知道有效集S(x*

9、)，因此只有想办法构造一个集合序列去逼近它，即从初始点 X0出发，计算-6 -有效集S(X0)，解对应的等式约束子问题。重复这一做法，得到有效集序列S(Xk), k二。,1，使之S(Xk)T S(X )，以获得原问题的最优解。基于上述定理，我们分4步来介绍有效集方法的算法原理和实施步骤。第一步，形成子问题并求出搜索方向dk.设Xk是问题(2.1)的一个可行点，据 此确定相应的有效集 Sk = E U |(xk)，其中I (Xk) =i |aiTXk - bi = 0,i忘I.求解相 应的子问题 . 1 T + T(2.3)(2.4)min - X Hx + c X, T 2s.t. ai X

10、- bi = 0, i 忘 Sk.上述问题等价于1T Tmin qk(d) =-d Hd +gkd,2s.t. ajd = 0, i 忘 Sk.其中X =Xk +4 -GXk +c.设求出问题(2.4)的全局极小点为dk,Xk是对应的拉格 朗日乘子。 (1)若Xk +dk是冋题(2.1)的可行点,则令叭 Hxkdt =Xk +dk. 若Xk +dk不是问题(2.1)的可行点，则通过线性搜索求出下降最好的可行点。注意到目标函数是凸二次函数，那么这一点应该在可行域的边界上达到。 因 此只要求出满足可行条件的最大步长otk即可。当 i - Sk 时，对于任意的 叭 0，都有 ajdk =0和 aiT

11、(Xk Ekdk) =aiTXk = bi，此时，叭0不受限制。当浮Sk时，即第i个约束是严格的不等式约束，此时要 求 a k 满足 aT (Xk +ctkdk)3bi，即GkaTdk b aTXk,i 忘注意到上式右端非正，故当aTdk二0时，上式恒成立。而当aTdO时，由上式-7 -可解得ai dk故有s=7k=min 卩；兀 laTdk o.I ai dk合并第(1)(2)可得o-k = min1,ot k (2.5).第三步，修正Sk.当ak =1时，有效集不变，即Sy = Sk.而当otkCl时,. T-bik -aikXk 叫k aikdk故aT(Xakdkbi ，因此在Xk十处增

12、加了一个有效约束，即k kSk 十：=Sk Ui k.第四步，考虑dk =O的情形。此时Xk是问题(2.3)的全局极小点。若这是对应的不等式约束的拉格朗日乘子均为非负， 则Xk也是问题(2.1 )的全局极小点，迭代终止。否则，如果对应的不等式约束的拉格朗日乘子有负的分量， 那么需要重新寻找一个下降可行方向。设几jk v。 jk引(Xk)，现在要求一个下降可行方向 dk，满足gidk 0,巧-I(Xk)，为简便计算，按下述方式选取 dk :aTk(Xk +dk)Abjk, aT(Xk+dk)=bjNj-Sk,jH jf adk AO,bdk5Sk,jHjk, 3另一方面，注意到xk是子问题(2.

13、3)的全局极小点，故有Hxk +c-送几:ai =0,idgk =Ak,其中Ak = ( ai )i Sk，几 k =(人)i 旳.从而，gldkAldk.由(2.6)知Aldk = (aTdk)ej =(aTkdk)ejk，jSk于是有gidk=沐(ajk dk )ej = Xjk(ajk dk) 0.上式表明，由(2.6)确定的dk疋曰一个下降可行方向。因此，令Sk=Sk jk，贝M修正后的子问题4min qk(d) = dTHd +gTd,2s.t. aTd = 0,i 壬 Sk的全局极小点必然是原问题的一个下降可行方向。222有效集方法的算法步骤经过上面的分析和推导，我们现在可写出有效

14、集方法的算法步骤:(1 )选取初值。给定初始可行点xo - Rn，令k =0.(2)解子问题。确定相应的有效集 Sk =eU |(xk).求解子问题min qk(dH1dTHgTd,s.t. ajd =0,i 忘 Sk,得极小点dk和拉格朗日乘子向量 几k右dk K0转步骤(4);否则，转步骤(3)检验终止准则。计算拉格朗日乘子“k = Bkgk,其中gk =HXk+c, Bk =(人HAT)AkHt 宀=)赳.k)t =m(in)0k)i.若仏k)t X0，则Xk是全局极小点，停算。否则，若(S)t ，则令Sk=Skt,转步骤(2)。(4)确定步长ak .令ak = min1,ot k，其中

15、a k = min Jb Xk | aT d k 0，.0 I ai dk i令 Xk+：=Xk +akdk.(5)若 Gk =1，则令 SkH1=Sk ；否则，若 a k 0, X2 二 0.解首先确定有关数据:利用Matlab计算可得结果如下:ex it flag =O0.50002.25003总结与体会通过本次实验，笔者对求解等式约束凸二次规划问题的拉格朗日方法及一般 约束条件下凸二次规划问题的有效集方法有了较深的认识和了解。感谢阮老师的悉心教诲和指导，使得笔者对最优化课程中的理论推导、 算法 步骤及编程都比较熟悉，对以后的科研工作有很好的指导和借鉴意义。4附录4.1拉格朗日方法的mat

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

17、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)有效集方法函数%本程序主要适用于求解一般约束条件下的凸二次规划问题fun cti on x,lamk,exitflag,out put=qp act(H,c,Ae,be,Ai,bi,xO)%功能：用有效集方法解一般约束二次规划问题%min f(x)=0.5*x*H*x+c*x,% s.t. a_i*x- b_i=0,

18、(i=1, ,l),% 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是对应的乘子向量；out put是结构变量 % 输出极小值f(x),迭代次数k等信息，exitflag是算法终止类型 %初始化ep silo n=1.0e-9;err=1.0e-6;k=0;x=x0 ;n=len gth(x);kmax=1.0e3;n e=le ngth(be); ni=le ng

19、th(bi);lamk=zeros (ne+n i,1);in dex=ones(n i,1);for i=1: niif Ai(i,:)*xbi(i)+e psilo nin dex(i)=0;endend%算法主程序while k0Aee=Ae;endfor j=1: niif in dex(j)0Aee=Aee;Ai(j,:);endendgk=H*x+c;m1, n1=size(Aee);dk,lamk=qsub p( H,gk,Aee,zeros(m1,1);if n orm(dk)v=eiTy=0.0;if len gth(lamk) ney,jk=mi n(lamk( ne+1:

20、le ngth(lamk);endif y=0exitflag=0;elseexitflag=1;for i=1: niif in dex(i )&(n e+sum(i ndex(1:i)=jk in dex(i)=0;break;endendendk=k+1;elseexitflag=1;%求步长al pha=1.0;tm=1.0;for i=1: niif (in dex(i)=0)&(Ai(i,:)*dk0) tm1=(bi(i)-Ai(i,:)*x)/(Ai(i,:)*dk);if tm1tmtm=tm1;ti=i;endendendalp ha=min(alp ha,tm);x=x+a Ip ha*dk;%修正有效集if tm0rb=Ae*gi nvH*c+be;lambda=pin v(Ae*gi nvH*Ae)*rb; x=gi nvH*(Ae*lambda-c);elsex=-gi nvH*c;lambda=0;end(3)有效集方法求解一般约束凸二次规划问题主程序cic;clear;H=2 -1;-1 4;c=-1 -10;Ae=;be=;Ai=-3 -2;1 0;0 1;bi=-6 0 0;xO=O 0;x,lambda,exitflag,out pu t=qp act(H,c,Ae,be,Ai,bi,xO)