基于MATLAB的可行方向法求極值問題講解

1、基于MATLAB可行方向法求極值的實現(xiàn)姓名：xxx學(xué)號：xxx（北京理工大學(xué)機(jī)械與車輛學(xué)院車輛工程，北京 100081）摘要：在工程實際的優(yōu)化設(shè)計中，隨著設(shè)計變量數(shù)和約束條件數(shù)的增加，隨機(jī)方向搜索法和復(fù)合形法等直接優(yōu)化解法的求解效率會偏低?？尚蟹较蚍?，顧名思義，一種始終在可行域內(nèi)尋找下降方向的搜索法，以其收斂速度快、效果好的優(yōu)點(diǎn)已成為求解約束非線性問題的一種有代表性的直接解法，同時也是求解大型約束優(yōu)化問題的主要方法之一。本文將簡單介紹可行方向法的數(shù)學(xué)思想，采用線性規(guī)劃法和約束最優(yōu)步長法編寫MATLAB程序，最后通過算例完成對優(yōu)化問題的求解。關(guān)鍵字：可行方向法；MATLAB；優(yōu)化方法。1. 可行

2、方向法的基本數(shù)學(xué)思想1.1可行方向法的搜索策略可行方向法迭代的第一步都是從可行域的某一初始點(diǎn)出發(fā)，沿負(fù)梯度方向移至某一個或J個起作用約束面的交集上。以后的搜索路線和迭代計算可根據(jù)約束函數(shù)和目標(biāo)函數(shù)的不同性狀，分別采用以下三種不同策略繼續(xù)搜索。1） 由點(diǎn)出發(fā)，沿可行方向作一維最優(yōu)化搜索，若所得新點(diǎn)在可行域內(nèi)，則再沿方向作一維最優(yōu)化搜索；若所得的新點(diǎn)不在可行域內(nèi)，則將它移至約束面上再反復(fù)重復(fù)上述步驟，若，則停止迭代，如圖1.1所示。2） 由點(diǎn)出發(fā)，沿可行方向作一維最優(yōu)化搜索，若所得新點(diǎn)在可行域外，則沿可行方向以最大步長到達(dá)另一個約束面上一點(diǎn)，將該點(diǎn)作為迭代點(diǎn)進(jìn)行反復(fù)搜索，直至滿足給出的K-T條件，

3、如圖1.2所示。3） 沿著約束面進(jìn)行搜索。對于只具有線性約束的非線性規(guī)劃問題，如圖1.3所示，從點(diǎn)出發(fā)，沿約束面移動，在有限的幾步內(nèi)即可搜索到約束最優(yōu)點(diǎn)；對于非線性約束函數(shù)，如圖1.4所示情況，就是沿著約束面的切線移動，但這樣將會進(jìn)入非可行域，使問題變得復(fù)雜。此時，應(yīng)將進(jìn)入非可行域的新點(diǎn)X設(shè)法調(diào)整回約束邊界。調(diào)整方法為：先規(guī)定約束面容差，建立新的約束邊界，然后將已經(jīng)離開約束面的X點(diǎn)沿起作用約束的負(fù)梯度方向返回約束面上，計算公式為式中，為調(diào)整步長，可以用試探法確定，也可以用下式估計，使其回到約束邊界上。 (1-1)圖1.1圖1.2圖1.3圖1.4不管采用哪種搜索方法都要進(jìn)行兩條決策：一是產(chǎn)生一個

4、適用的可行方向；二是沿方向確定一個不會越出可行域的適當(dāng)?shù)牟介L因子。1.2產(chǎn)生可行方向和步長的數(shù)學(xué)原理可行方向要保證沿該方向作微小移動后所得的新點(diǎn)是可行點(diǎn)，且目標(biāo)函數(shù)值有所下降。顯然可行方向應(yīng)滿足可行和下降兩個條件。這里從數(shù)學(xué)原理的角度分析可行性和下降性。1.2.1可行性設(shè)為可行域D中的一個點(diǎn)，即。對于某一方向來說，若存在數(shù)，使其對于任意的A，式1-2均成立，則稱方向?qū)c(diǎn)滿足可行性條件。 （1-2）若點(diǎn)在可行域內(nèi)，即為內(nèi)點(diǎn)，則任何方向都滿足可行性條件。當(dāng)點(diǎn)在某一約束邊界上時，如圖1.5所示，該約束面就是起作用約束。當(dāng)點(diǎn)位于一個約束面上時，如圖1.6所示，作出初起作用約束的梯度和切線，可見只要與起

5、作用約束的梯度成直角或鈍角，則可指向可行域內(nèi)，滿足可行性條件。當(dāng)同時處于J個約束面時，就要求與這J個約束面的梯度均垂直或交成鈍角，如圖1.6所示，寫成數(shù)學(xué)表達(dá)式如下：式中，J為起作用約束的個數(shù)。上式稱為可行性條件。圖1.5圖1.61.2.2下降性沿一方向搜索時，要求下降愈快愈好。對某一點(diǎn)來說，負(fù)梯度方向為最速下降方向。如果負(fù)梯度方向是可行方向，那么負(fù)梯度方向是最有利的方向；否則為了保證目標(biāo)函數(shù)值有所下降，至少要使搜索方向與目標(biāo)函數(shù)的負(fù)梯度方向成銳角或與梯度方向成鈍角，寫成數(shù)學(xué)表達(dá)式為該式稱為下降性條件。如圖1.7所示，可行下降方向顯然位于點(diǎn)約束面的切線與目標(biāo)函數(shù)等值線的切線所圍成的扇形區(qū)域內(nèi)，

6、推廣到一般的情況就是可行方向在目標(biāo)函數(shù)超等值面的超切面和J個起作用約束的超切面所圍成的超錐體內(nèi)。該區(qū)域稱為可行下降方向區(qū)。圖1.7綜上所述，當(dāng)點(diǎn)位于J個起作用的約束面上時，滿足 （1-3）式中，即為可行方向。1.3可行方向的產(chǎn)生本文采用線性規(guī)劃法來產(chǎn)生可行方向，用約束最優(yōu)步長確定最優(yōu)步長，為后面編制MATLAB程序做理論基礎(chǔ)。本節(jié)先闡述線性規(guī)劃法，下節(jié)再對約束最優(yōu)步長法進(jìn)行描述。線性規(guī)劃法對包含線性和非線性的不等式約束的最優(yōu)化問題都適用，但不允許有等式約束。其基本原理是將具有一階連續(xù)偏導(dǎo)數(shù)的目標(biāo)函數(shù)和約束條件在點(diǎn)用Taylor展開式展成線性近似函數(shù)（一次項），并用這些線性近似函數(shù)代替目標(biāo)函數(shù)和

7、它的約束條件，使得問題線性化。這樣成為約束條件變?yōu)橛么嫔鲜街械腦，得到s.t. 式中，、為常數(shù)，只有搜索方向和是未知量。當(dāng)?shù)c(diǎn)位于J個起作用約束邊界的交點(diǎn)上時，問題就變成求解線性規(guī)劃問題：s.t. 式中，只起到方向作用，故規(guī)定其向量的模是有界的，也就是規(guī)定每個分量的絕對值不大于1，即由于可行方向還得滿足式1-3的可行性和下降性條件，故確定可行方向的線性規(guī)劃數(shù)學(xué)模型如下：s.t. (1-4)求解后，可得到迭代方向，作為下一步的搜索方向。1.4最優(yōu)步長的確定由線性規(guī)劃法求的可行方向后，由出發(fā)，找到該方向上的可行新點(diǎn)。新點(diǎn)必須滿足兩個條件，一是新點(diǎn)必須還是可行點(diǎn)，二是目標(biāo)函數(shù)值具有最大的下降量。

8、約束最優(yōu)步長即可滿足這兩點(diǎn)，可用以下一維搜索法求解：s.t. 式中，。2. 可行方向法的迭代過程和算法流程圖2.1迭代過程1)輸入收斂精度；2)形成可行初始點(diǎn)，令k=1;3)確定起作用約束下標(biāo)集合，；4)判斷E是否為空集，若是空集則說明在可行域內(nèi)，轉(zhuǎn)步驟5）；否則說明在約束邊界上，轉(zhuǎn)步驟6)；5)判斷是否滿足終止條件，若滿足，轉(zhuǎn)步驟10)；若不滿足，則取搜索方向為，轉(zhuǎn)步驟8)；6)采用線性規(guī)劃法，用式1-4求解后，得到可行方向；7)判斷是否成立，若成立，則轉(zhuǎn)步驟10)輸出，停機(jī)結(jié)束；若不成立，則繼續(xù)下一步；8)采用約束最優(yōu)步長法進(jìn)行一維搜索，得到新點(diǎn);9)令k=k+1，轉(zhuǎn)步驟3)；10)輸出，

9、停機(jī)結(jié)束。2.2算法流程圖算法流程圖如圖2-1。圖 2-1 可行方向法的算法流程圖3. 可行方向法的MATLAB程序function kxfxf(x0)eps=1.0e-6;%剛開始給的x0為行向量x0=transpose(x0);funcsz=length(x0);G1 = G;for k=1:1:50 %'f(x)梯度：' sf=diff_val(x0); sf=eval('sf'); %'f(x)梯度范數(shù)：' norm_s=norm(diff_val(x0); GL=length(G); G1 = G; G_copy=G; for i=1:

10、1:GL g=subs(G1(i,:),x,x0); G1(i,:)=g; end G_zero=eval('G1'); for i=GL:-1:1 if abs(G_zero(i,:)>0.1 G_zero(i,:)=; G_copy(i,:)=; end end I=size(G_zero,1); add=-ones(I,1); %根據(jù)E是否為空集(I=0)確定不同情況下的搜索方向 if I=0 if norm_s<=3 '最優(yōu)解：' '自變量值:' x0 '目標(biāo)函數(shù)值：' f=fval(x0) '迭代次數(shù)

11、：' k'=' return else d0=-diff_val(x0); %線性搜索 d0=transpose(d0); vmax = gmin(x0,d0); h=fmin(x0,d0,vmax); x0=x0+h*d0; end else %線性規(guī)劃 grad=jacobian(G_copy,x); G_zero=subs(grad,x,x0); Ac=sf;G_zero; lb=-1*ones(sz,1); ub=ones(sz,1); p=zeros(I+1,1); dz = size(1,sz); c=sf*dz' dz,mn,m1,m2,m3=li

12、nprog(c,Ac,p,lb,ub); z0 = sf*dz; z0=abs(z0); if z0<0.01 '最優(yōu)解：' '自變量值:' x0 '目標(biāo)函數(shù)值：' f=fval(x0) '迭代次數(shù)：' k'=' return else d0=dz; %線性搜索 vmax = gmin(x0,d0); h=fmin(x0,d0,vmax); x0=x0+h*d0; end endend程序中用到一些回調(diào)函數(shù)，簡單列舉如下，.m文件詳細(xì)內(nèi)容見附錄一。func記錄函數(shù)及其自變量信息fval(x0)計算函數(shù)在x0處

13、得函數(shù)值gmin(x0,d0) 求函數(shù)在滿足約束條件下的最小值fmin(x0,d0,vmax)求函數(shù)在0,vmax上的最小值diff_val(x0)計算函數(shù)在x0處得導(dǎo)數(shù)值。4.利用MATLAB求解實例4.1實例本文以課本例題5.3為例來驗證可行方向法在MATLAB里的實現(xiàn)。已知目標(biāo)函數(shù)，約束條件為s.t. 初始點(diǎn)為，求的極小值，計算精度。首先在func.m文件中輸入目標(biāo)函數(shù)和約束條件，如下：f=4*x1+x22-12;g1=x12+x22-25;g2=-x2;g3=-x1;G = g1;g2;g3;在MATLAB命令窗口里輸入如下命令：>> x0=3 2>> kxfx

14、f(x0)計算結(jié)果為：最優(yōu)解：自變量值:x0 = 1.0e-03 *ans =目標(biāo)函數(shù)值：f = -11.999735470616054ans =迭代次數(shù)：k = 3與書中提供的精確解基本一致，最優(yōu)解的相對誤差為， 說明此方法有效可行。4.2實例結(jié)果分析此處主要從以下三個方面來分析此算法及計算結(jié)果。1) 可靠性。通過對其他多元多次函數(shù)（教科書例5.1）進(jìn)行求極值計算，結(jié)果都比較精確，說明此算法可靠性、通用性好。通過對教課書例5.1進(jìn)行求解，可得以下結(jié)果：x0 = 6.999853314348210目標(biāo)函數(shù)值：f = 2.000586785240415迭代次數(shù)：k = 3與課本上精確解相比，幾乎

15、完全吻合。2) 有效性。對于本例，目標(biāo)函數(shù)不是很復(fù)雜，其求解速率在1s內(nèi)可完成，迭代次數(shù)為k=3次。3) 簡便性。此程序邏輯清晰，簡明易懂，對計算機(jī)的要求也不是很高，其簡便性不言而明。參考文獻(xiàn)1 李志鋒。機(jī)械優(yōu)化設(shè)計。高等教育出版社。2 孫靖民。機(jī)械優(yōu)化設(shè)計。機(jī)械工業(yè)出版社。1999.5。3 蒲俊，吉家鋒，伊良忠。MATLAB 6.0數(shù)學(xué)手冊。浦東電子出版社。2002.1。附錄一 回調(diào)函數(shù)內(nèi)容1. func 記錄函數(shù)及其自變量信息根據(jù)所求解的目標(biāo)函數(shù)和相應(yīng)的約束條件不同而不同，由用戶進(jìn)行手動輸入。示例見4.1。2. fval(x0) 計算函數(shù)在x0處得函數(shù)值function f_val=fva

16、l(x0)x0=transpose(x0);func;f_val=subs(f,x,x0);3. gmin(x0,d0) 求函數(shù)在滿足約束條件下的最小值function h=gmin(x0,d0)syms h;funcn = length(G);alpha = zeros(1,n);a=x0+h*d0;g_val=subs(G,x,a);k = 1;for i=1:1:n alpha_temp = solve(g_val(i); if length(alpha_temp)>1 if imag(alpha_temp) = 0 for j=1:1:length(alpha_temp) if

17、alpha_temp(j,1)>0 | alpha_temp(j,1) = 0 alpha(1,k) = alpha_temp(j,1); k = k+1; end end else end else if alpha_temp > 0 | alpha_temp = 0 alpha(1,k) = alpha_temp; k = k+1; else alpha(:,k) = ; end endendh_temp = alpha(1,1);len = length(alpha);for i=1:1:len if h_temp > alpha(1,i) h_temp = alpha(1,i); endendh=h_temp;4. fmin(x0,d0