外點(diǎn)罰函數(shù)優(yōu)化實(shí)例

。-可編輯修改-外點(diǎn)罰函數(shù)優(yōu)化實(shí)例一、優(yōu)化問(wèn)題如圖1所示，為某一桁架的一部分，桿2距O點(diǎn)30cm處有一支點(diǎn)C。為了固定桁架，現(xiàn)欲在桿1和2上設(shè)置支點(diǎn)A和B，用來(lái)連接桿3（可拆卸）。已知當(dāng)桁架固定時(shí)，桿1和2成直角；而且，桿1右邊有一段長(zhǎng)為20cm的重要部位，不能設(shè)置支點(diǎn)。卸去桿3、收起桁架時(shí)，支點(diǎn)A的位置不能高于BC段中點(diǎn)D。求取支點(diǎn)A、B的位置，使得桿3的長(zhǎng)度盡量小，以節(jié)省材料。圖1桁架結(jié)構(gòu)示意圖二、數(shù)學(xué)模型設(shè)A、B兩點(diǎn)距離O點(diǎn)的長(zhǎng)度分別為和，而桁架固定時(shí)桿1和2成直角。所以桿3的長(zhǎng)度為。由圖1可知，且，即且。設(shè)，取。因此，數(shù)學(xué)模型為：極小化目標(biāo)函數(shù)約束條件三、求解數(shù)學(xué)模型(1)外點(diǎn)罰函數(shù)法求解構(gòu)造外點(diǎn)法罰函數(shù)，如下：程序流程圖如圖2所示：給定、給定、、c、k=0i=0求與Hessian矩陣輸出和YNi=i+1k=k+1YN結(jié)束牛頓法求的極值點(diǎn)圖2外點(diǎn)罰函數(shù)法程序流程圖程序步驟：①選擇適當(dāng)?shù)某跏剂P因子、初始點(diǎn)、收斂精度和罰因子系數(shù)c。在本程序中分別取，，，c=8。令迭代步數(shù)k=0。②采用牛頓法求無(wú)約束問(wèn)題的極值點(diǎn)。③檢驗(yàn)迭代終止準(zhǔn)則，若滿足及則停止迭代計(jì)算，輸出最優(yōu)點(diǎn)；否則，轉(zhuǎn)入步驟④。④取，，k=k+1，轉(zhuǎn)入步驟②繼續(xù)迭代。具體程序請(qǐng)看附一。運(yùn)行結(jié)果：X*=[20.0000，10.0000]f(X*)=500.0000因此，A、B兩支點(diǎn)與O的距離分別為20cm、10cm，桿3的最小長(zhǎng)度為cm。目標(biāo)函數(shù)曲線圖與目標(biāo)函數(shù)等值線圖分別如圖3和圖4所示。優(yōu)化路徑如圖4所示。圖3目標(biāo)函數(shù)曲線圖圖4目標(biāo)函數(shù)等值線圖(2)Matlab優(yōu)化工具fmincon求解利用Matlab文件編輯器為目標(biāo)函數(shù)編寫M文件(goalfun.m)：functionf=goalfun(x)f=x(1)^2+x(2)^2;編寫約束條件的M文件(confun.m)：function[c,ceq]=confun(x)c=[2*x(1)-x(2)-30;20-x(1)];ceq=[];編寫主函數(shù)的M文件(opt.m)：closeallclearallclcx0=[20,20];lb=[];ub=[];options=optimset('LargeScale','off','display','iter','tolx',1e-6);[x,fval,exitflag,output]=fmincon('goalfun',x0,[],[],[],[],[],[],'confun',options);xfval運(yùn)行結(jié)果：x=[20，10]fval=500同樣地，利用Matlab優(yōu)化工具解得，支點(diǎn)A、B與O的距離分別為20cm、10cm，桿3的最小長(zhǎng)度為cm。四、結(jié)論分析(1)罰因子系數(shù)c對(duì)外點(diǎn)罰函數(shù)法的影響本次程序中，c取值為8，運(yùn)行步數(shù)k=11。若取c=4，則運(yùn)行步數(shù)k=16；取c=16，則運(yùn)行步數(shù)k=9；取c=32，則運(yùn)行步數(shù)k=8；取c=64，則運(yùn)行步數(shù)k=7。由此可知，罰因子系數(shù)c的大小會(huì)影響程序的迭代次數(shù)k。c的值取得越大，運(yùn)行步數(shù)k越小，程序收斂速度越快，效率越高。但對(duì)于c的其他一些取值，如5、7、9等，會(huì)導(dǎo)致罰函數(shù)形態(tài)變壞，使迭代出現(xiàn)問(wèn)題，導(dǎo)致程序運(yùn)行失敗。因此，需選取合適的罰因子系數(shù)c。(2)外點(diǎn)罰函數(shù)法與Matlab優(yōu)化工具fmincon的比較通過(guò)opt.m的運(yùn)行結(jié)果可知，fmincon的運(yùn)行步數(shù)k=3，這一運(yùn)行效率明顯比外點(diǎn)法函數(shù)法的效率高。對(duì)于相同的收斂精度和初始點(diǎn)，雖然優(yōu)化結(jié)果相同，但是M文件WaiDianNiuDun.m中外點(diǎn)罰函數(shù)法的運(yùn)行效率仍有待提高。附一外點(diǎn)罰函數(shù)matlab程序closeallclearallclcsymsx1x2M;%M為罰因子。m(1)=1;c=8;%c為遞增系數(shù)。賦初值。a(1)=20;b(1)=20;f=x1^2+x2^2+M*((20-x1)^2+(2*x1-x2-30)^2);%外點(diǎn)罰函數(shù)f0(1)=500;%求偏導(dǎo)、Hessian元素fx1=diff(f,'x1');fx2=diff(f,'x2');fx1x1=diff(fx1,'x1');fx1x2=diff(fx1,'x2');fx2x1=diff(fx2,'x1');fx2x2=diff(fx2,'x2');%外點(diǎn)法M迭代循環(huán)fork=1:100x1=a(k);x2=b(k);M=m(k);%牛頓法求最優(yōu)值forn=1:100f1=subs(fx1);%求解梯度值和Hessian矩陣f2=subs(fx2);f11=subs(fx1x1);f12=subs(fx1x2);f21=subs(fx2x1);f22=subs(fx2x2);if(double(sqrt(f1^2+f2^2))<=1e-6)%最優(yōu)值收斂條件a(k+1)=double(x1);b(k+1)=double(x2);f0(k+1)=double(subs(f));break;elseX=[x1x2]'-inv([f11f12;f21f22])*[f1f2]';x1=X(1,1);x2=X(2,1);endendif(double(sqrt((a(k+1)-a(k))^2+(b(k+1)-b(k))^2))<=1e-6)&&(double(abs((f0(k+1)-f0(k))/f0(k)))<=1e-6)%罰因子迭代收斂條件%輸出最優(yōu)點(diǎn)坐標(biāo)，罰因子迭代次數(shù)，最優(yōu)值a(k+1)b(k+1)kf0(k+1)break;elsem(k+1)=c*m(k);endend%繪制目標(biāo)函數(shù)曲線圖xx1=0:0.5:50;xx2=0:0.5:50;fori=1:length(xx1)forj=1:length(xx2)if((2*xx1(i)-xx2(j)-30<=0)&&(20-xx1(i)<=0))Z(i,j)=xx1(i)^2+xx2(j)^2;elseZ(i,j)=0;endendendfigure(1);surf(xx1,xx2,Z);axis([05005004500])title('目標(biāo)函數(shù)曲線圖');xlabel('x1');ylabel('x2');%繪制目標(biāo)函數(shù)等值線圖，并畫出優(yōu)化路徑figure(2);x11=-5:0.5:25;x12=-5:0.5:25;[xx11,xx12]=meshgrid(x11,x12);F=xx11.^2+xx12.^2;axis