求解二次規(guī)劃問題的拉格朗日及有效集方法

1、求解二次規(guī)劃問題的拉格朗日及有效集方法最優(yōu)化方法課程實驗報告學 院:數(shù)學與統(tǒng)計學院班 級：碩2041班學 號：指導教師：同組人：錢東東求解二次規(guī)劃問題的拉格朗日及有效集方法摘要二次規(guī)劃師非線性優(yōu)化中的一種特殊情形，它的目標函數(shù)是二次實函數(shù)，約 束函數(shù)都是線性函數(shù)。由于二次規(guī)劃比較簡單，便于求解（僅次于線性規(guī)劃）， 并且一些非線性優(yōu)化問題可以轉(zhuǎn)化為求解一些列的二次規(guī)劃問題，因此二次規(guī)劃 的求解方法較早引起人們的重視，稱為求解非線性優(yōu)化的一個重要途徑。二次規(guī) 劃的算法較多，本文僅介紹求解等式約束凸二尺規(guī)劃的拉格朗日方法以及求解一 般約束凸二次規(guī)劃的有效集方法。關(guān)鍵字：二次規(guī)劃，拉格朗日方法，有效集

2、方法?！灸夸洝?TOC o 1-5 h z 摘要-1- HYPERLINK l bookmark16 o Current Document 1等式約束凸二次規(guī)劃的解法-3- HYPERLINK l bookmark19 o Current Document 1.1問題描述-3- HYPERLINK l bookmark22 o Current Document 1.2拉格朗日方法求解等式約束二次規(guī)劃問題-3- HYPERLINK l bookmark25 o Current Document 1.2.1拉格朗日方法的推導-3- HYPERLINK l bookmark31 o Current

3、Document 1.2.2拉格朗日方法的應用-4- HYPERLINK l bookmark36 o Current Document 2一般凸二次規(guī)劃問題的解法-5- HYPERLINK l bookmark39 o Current Document 2.1問題描述-5- HYPERLINK l bookmark43 o Current Document 2.2有效集法求解一般凸二次規(guī)劃問題-6- HYPERLINK l bookmark46 o Current Document 2.2.1有效集方法的理論推導-6- HYPERLINK l bookmark49 o Current Doc

4、ument 2.2.2有效集方法的算法步驟-9-2.2.3有效集方法的應用-10-3總結(jié)與體會-11-4附錄-11- HYPERLINK l bookmark61 o Current Document 4.1拉格朗日方法的matlab程序-11- HYPERLINK l bookmark73 o Current Document 4.2有效集方法的Matlab程序-11-1等式約束凸二次規(guī)劃的解法1.1問題描述我們考慮如下的二次規(guī)劃問題f 1 ，-min 2xtHx + ctx, ( s.t. Ax = b其中H e Rnxn對稱正定，A e Rmxn行滿秩，c,x G Rn，b e Rm。1

5、.2拉格朗日方法求解等式約束二次規(guī)劃問題1.2.1拉格朗日方法的推導首先寫出拉格朗日函數(shù)：L(x,人)=2 xtHx + ctx -Xt (Ax - b)， (1.2)令V L(x, X) = 0,V/(x, X) = 0，得到方程組Hx - ATX = -c,-Ax=-b.將上述方程組寫成分塊矩陣形式：；At x X=- c -b.(1.3)我們稱傷處方程組的系數(shù)矩陣H -AT-A 0為拉格朗日矩陣。下面的定理給出了線性方程組(1.1)有唯一解的充分條件。定理1設H e Rmxn對稱正定，A e Rm、n彳丁滿秩。若在問題(1.1)的解x*處滿 足二階充分條件，即dTHd 0, Vd e R

6、n, d。0, Ad = 0,則線性方程組(1.4)的系數(shù)矩陣非奇異，即方程組(1.4 )有唯一解。其中，方程 組(1.4)為(1.1)對應的齊次方程組：H-At3-a0V=0 (1.4).故可設其逆為下面，我們來推導方程（1.3）的求解公式。根據(jù)定理1，拉格朗日矩陣必然 是非奇異的一 H-At=一 G-Bt -A0-Bc由恒等式H-At G-Bt -=_ In0 一nxm-a0 JL- BC _0mxnIm可得HG + At B = I-HBt - AtC = 0nxm-AG = 0mxnABt = Im于是由上述四個等式得到矩陣G,B,C的表達式G = H-1 - H-iAt (AH-i

7、At )-i AH-1,(1.5)B = (AH -i At ) -i AH -i,(1.6)(1.7)C = -( AH -i At ) -i.因此，由（1.3）可得解得表達式X-Bt 一r- c - Gc + BTb=-Bc _-b=Bc - Cb(1.8)其中，G,B,C分別由（1.5）,（1.6）,（1.7）給出。下面給出X和廠的另一種等價表達式。設x是問題（1.1）的任一可行點，即 kXk滿足Axk = n 而在此點處目標函數(shù)的梯度為gk=W（七）=Hx +c，利用Xk(1.9)和gk，可將（1.8）改寫為廠1.2.2拉格朗日方法的應用（1）拉格朗日方法的Matlab程序見附錄。（2

8、）利用拉格朗日方法求解下列問題：min X2 + 2X2 + 尤2 2x x + x , s.t. x + x + x = 4, 2x - x + x = 2.解容易寫出一 2-20O114H =-240,c =0,A =,b =2-112002111*利用Matlab程序求解該問題可以結(jié)果如下：K 二1.9090909090909091.9545454545454550. 136363636363637Lam =2.636363636363636-1.363636363036363val 二3.9772727272727282 一般凸二次規(guī)劃問題的解法2.1問題描述考慮一般二次規(guī)劃I 1 ，

9、min xtHx + ctx,(2.1)s.t. aTx-b = 0,i e E = 1,l,aTx-b 0,i e I = l +1,mii其中H是n階對稱陣。記I(x*) = ilaTx*-bi = 0,i e I，下面的定理給出了問題 (2.1)的一個最優(yōu)性充要條件。定理2 *是二次規(guī)劃問題(2.1)的局部極小點當且僅當(1)存在X* e Rm，使得Hx+c&a=0,ieEieIaTx* - b = 0, i e E,iiaTx* - b 0,i e I,X* 0, i e I; X* = 0, i e I I (x *) ii記S = d e Rn 0ldTa. = 0,i e E;d

10、Ta. 0,i e I(x*);dTa. = 0,i e I(x*)且X* 0.則對于任意的d e S，均有dTHd 0.容易發(fā)現(xiàn)，問題(2.1)是凸二次規(guī)劃的充要條件是H半正定。此時，定理2 的第二部分自然滿足。注意到凸優(yōu)化問題的局部極小點也是全局極小點的性質(zhì)， 我們有下面的定理：定理3 x*是凸二次規(guī)劃的全局極小點的充要條件是x*滿足KT條件，即存在X* e Rm，使得Hx* + c-ZX*ai -ZX*ai = 0,ieEieIaTx* - b = 0,i e E, iiaTx* - b 0,i e I,iiX* 0,i e I;X* = 0,i e I I(x*).ii2.2有效集法求

11、解一般凸二次規(guī)劃問題2.2.1有效集方法的理論推導首先引入下面的定理，它是有效集方法理論基礎。定理4設x*是一般凸二次規(guī)劃問題的全局極小點，且在x*處的有效集為S(x*)= E I(x*)，則x*也是下列等式約束凸二次規(guī)劃min xtHx + ctx,c、 0, i e I.則令a = 1, x = x + d .(2)若xk + dk不是問題(2.1)的可行點，則通過線性搜索求出下降最好的可行 點。注意到目標函數(shù)是凸二次函數(shù)，那么這一點應該在可行域的邊界上達到。因 此只要求出滿足可行條件的最大步長a即可。k當i e S時，對于任意的a 0，都有aTd = 0和aT (x +a d ) = a

12、Tx = b， kki ki k k k i k i此時，ak 0不受限制。當i冬S/寸，即第i個約束是嚴格的不等式約束，此時要 求a 滿足aT (x +a d ) b，艮口a aTd b 一 aTx ,i e S .k i k i i kk注意到上式右端非正，故當aTdk 0時，上式恒成立。而當aTdk 0時，由上式-7 -可解得a b ”.kaTd故有=a = min 一氣 X | aTd .k aTd l k J、 l k/合并第(1)(2)可得第三步，修正S .當a k故 aT 3 +a d ) = blklkaTdj k=0, Vj e E; aTd jk 0, Vj e I(xk)

13、，為簡便計算，按下述方式選取dk :a = min1,a k(2.5).=1時，有效集不變，即七：=S而當a k 1時，b - aTxak ak ，kaTdk k， l k k因此在xk+(處增加了一個有效約束，即 第四步，考慮dk = 0的情形。此時x是問題(2.3)的全局極小點。若這是對應的不等式約束的拉格朗日乘子均為非負，則七也是問題(2.1)的全局極小點，迭代 終止。否則，如果對應的不等式約束的拉格朗日乘子有負的分量，那么需要重新 尋找一個下降可行方向。kk設七廣0, jk e I (x )，現(xiàn)在要求一個下降可行方向 d，滿足gTd b ,jk k kjkaT (x + d ) = b

14、 , Vj e S , j。j ,j k k jkk|aTdk 0,ak = 0：Vj e Sk, j 牛 j,a*)另一方面，注意到xk是子問題(2.3)的全局極小點，故有Hx + c Z 人k a = 0,leSk其中氣=(氣)心,廣。:)心kk從而，gTd = XtAtd .由(2.6)知 k k k k kATd =Z (aTd )e = (aT d )e于是有jeSkJJJkJkgTd = Xt (qt d )e, =Xk (aT d ) 0，則x是全局極小點，停算。否則，若(X ) 0，則令S := S t， ktkk tk k轉(zhuǎn)步驟(2)。確定步長a .令a = min1,ak，

15、其中=min心k否則，若a V1，k則令 Sk+1 = Sk/J，其b 一件 k | aTd v 0.aTd， k J中匕滿足 b - aT xk aT djk k（6）令k := k +1，轉(zhuǎn)步驟（1）。2.2.3有效集方法的應用有效集法的Matlab程序見附錄。用有效集方法求解下列二次規(guī)劃問題:min f (x) = x 2 xx + 2x 2 x 一 10 x , 11 2212 0, x2 0.解首先確定有關(guān)數(shù)據(jù)：- 3-2-62-1-1H =,c =,Ae = , be = , Ai =10,bi =0-14-1011 1010利用Matlab計算可得結(jié)果如下:X =0. 50002

16、. 250Cesitflag =0lambda =output -0.750Cfva.1: -13.7500iter: 83總結(jié)與體會通過本次實驗，筆者對求解等式約束凸二次規(guī)劃問題的拉格朗日方法及一般 約束條件下凸二次規(guī)劃問題的有效集方法有了較深的認識和了解。感謝阮老師的悉心教誨和指導，使得筆者對最優(yōu)化課程中的理論推導、算法 步驟及編程都比較熟悉，對以后的科研工作有很好的指導和借鑒意義。4附錄4.1拉格朗日方法的matlab程序拉格朗日算法函數(shù)%本程序用拉格朗日方法求解燈飾約束條件的二次規(guī)劃問題。function x,lam,fval=qlag(H,A,b,c)%功能：用拉格朗日方法求解等式約

17、束二次規(guī)劃：% min f(x)=0.5*xHx+cx，s.t. Ax=b%輸入：H,c分別是目標函數(shù)的矩陣和向量，A,b分別是約束條件中的矩陣和向量%輸出：(x,lam)是KT點，fval是最優(yōu)值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=B*c-C*b;fval=0.5*x*H*x+c*x;拉格朗日算法求解等式約束凸二次規(guī)劃問題主程序：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

18、,fval=qlag(H,A,b,c)4.2有效集方法的Matlab程序(1)有效集方法函數(shù)%本程序主要適用于求解一般約束條件下的凸二次規(guī)劃問題function x,lamk,exitflag,output=qpact(H,c,Ae,be,Ai,bi,x0)%功能：用有效集方法解一般約束二次規(guī)劃問題%min f(x)=0.5*x*H*x+c*x,% s.t. a_i*x-b_i=0,(i=1,.,l),%a_i*x-b_i=0,(i=l+1,.,m)%輸入：x0是初始點，H,c分別是目標函數(shù)二次矩陣和向量；%Ae=(a_1,.,a_l),be=(b_1,.,b_l);%Ai=(a_l+1,.,

19、a_m),bi=(b_l+1,.,b_m).%輸出：x是最優(yōu)解，lambda是對應的乘子向量；output是結(jié)構(gòu)變量%輸出極小值f(x)，迭代次數(shù)k等信息，exitflag是算法終止類型%初始化epsilon=1.0e-9;err=1.0e-6;k=0;x=x0;n=length(x);kmax=1.0e3;ne=length(be);ni=length(bi);lamk=zeros(ne+ni,1);index=ones(ni,1);for i=1:niif Ai(i,:)*xbi(i)+epsilonindex(i)=0;endend%算法主程序while k0Aee=Ae;endfor j=1:niif index(j)0Aee=Aee;Ai(j,:);endendg