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

一.優(yōu)化問(wèn)題

如圖1所示，為某一桁架的一部分，桿2距。點(diǎn)30cm處有一支點(diǎn)C。為了固定

桁架，現(xiàn)欲在桿1和2上設(shè)置支點(diǎn)A和B,用來(lái)連接桿3（可拆卸1已知當(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é)省材料。

C■-------

桿2

D,■—A，

/

/B,?-T軟

/桿:V,

宴

桿1I/

AL20’0

一XI一

圖1桁架結(jié)構(gòu)示意圖

二.數(shù)學(xué)模型

設(shè)A.B兩點(diǎn)距離。點(diǎn)的長(zhǎng)度分別為口和口，而桁架固定時(shí)桿1和2成直角。所

以桿3的長(zhǎng)度為口。

由圖1可知,□且口即口且口。

設(shè)口取口。因此數(shù)學(xué)模型為：

極小化目標(biāo)函數(shù)

niin/(X)=X12+x}

約束條件

4(X)=F+20V0

^2(X)=2Xj-x2-30<0

三.求解數(shù)學(xué)模型

(1)外點(diǎn)罰函數(shù)法求解

構(gòu)造外點(diǎn)法罰函數(shù)，如下：

22

0(X,M⑹)=%；+君+例"){[(max(20-x,,0)J+[皿雙2XI-X2-30,0)J)

程序流程圖如圖2所示：

圖2外點(diǎn)罰函數(shù)法程序流程圖

程序步驟：

①選擇適當(dāng)?shù)某跏剂P因子口、初始點(diǎn)口、收斂精度□和罰因子系數(shù)C。在本程序

中分別取口,□,□,c=8e令迭代步數(shù)k=0o

②采用牛頓法求無(wú)約束問(wèn)題minX,")的極值點(diǎn)X'(M⑹)。

目標(biāo)函數(shù)等值殘

25

圖4目標(biāo)函數(shù)等值線圖

(2)Matlab優(yōu)化工具fmincon求解

利用Matlab文件編輯器為目標(biāo)函數(shù)編寫(xiě)M文件(goalfun.m):

functionf=goalfun(x)

f=xCL)八2+x⑵八2;

編寫(xiě)約束條件的M文件(confun.m):

function[c,ceq]=confun(x)

c=[2*x(l)-x⑵-30;

20-x(l)];

ceq=[];

編寫(xiě)主函數(shù)的M文件(opt.m):

closeall

clearall

cic

xO=[20,20];

lb=[]；

ub=[];

,,,,,,

options=optimset(LargeScale/'off',display',iter',tolx,le-6);

,

[x,fval,exitflag,output]=fmincon(goalfun'/xO,[],[],[],[],[],[],'confun',optio

ns);

x

fval

運(yùn)行結(jié)果：

x=[20,10]

fval=500

同樣地，利用Matlab優(yōu)化工具解得，支點(diǎn)A.B與O的距離分別為20cm、

10cmz桿3的最小長(zhǎng)度為Clem。

四.結(jié)論分析

Q)罰因子系數(shù)c對(duì)外點(diǎn)罰函數(shù)法的影響

本次程序中,c取值為8,運(yùn)行步數(shù)k=llo若取c=4,則運(yùn)行步數(shù)k=16;取c=16,

則運(yùn)行步數(shù)k=9;取c=32,則運(yùn)行步數(shù)k=8;取c=64r則運(yùn)行步數(shù)k=70由此

可知，罰因子系數(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ò)。pt.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程序

closeall

clearall

cic

symsxlx2M;%M為罰因子。

m(l)=l;

c=8;%c為遞增系數(shù)。賦初值。

a(l)=20;

b(l)=20;

f=xl人2+x2人2+M*((20-xl)A2+(2*xl-x2-30)八2);%夕卜點(diǎn)罰函數(shù)

f0(l)=500;

%求偏導(dǎo)、Hessian元素

fxl=diff(f；xl');

fx2=diff(f,'x2')；

fxlxl=diff(fxl,'xl');

fxlx2=diff(fxl,'x2');

fx2xl=diff(fx2,'xl');

fx2x2=diff(fx2；x2')；

%外點(diǎn)法M迭代循環(huán)

fork=l:100

xl=a(k);x2=b(k);M=n(k);

%牛頓法求最優(yōu)值

forn=l:100

fl=subs(fxl);%求解梯度值和Hessian矩陣

f2=subs(fx2);

fll=subs(fxlxl);

fl2=subs(fxlx2);

f21=subs(fx2xl);

f22=subs(fx2x2);

if(double(sqrt(flA2+f2A2))<=le-6)%最優(yōu)值收斂條

件

a(k+l)=double(xl);b(k+l)=double(x2);f0(k+l)=double(subs(f));

break;

else

X=[xlx2],-inv([fllfl2;f21f22])*[flf2]';

xl=X(l,l);x2=X(2,l);

end

end

if(double(sqrt((a(k+l)-a(k))A2+(b(k+l)-b(k))A2))<=le-6)&&(double(abs((fO(k+l)-fO(

k))/f0(k)))<=le-6)%罰因子迭代收斂條件

%輸出最優(yōu)點(diǎn)坐標(biāo)，罰因子迭代次數(shù)，最優(yōu)值

a(k+l)

b(k+l)

k

fO(k+l)

break;

else

m(k+l)=c*m(k);

end

end

%繪制目標(biāo)函數(shù)曲線圖

xxl=0:0.5:50;

xx2=0:0.5:50;

fori=l:length(xxl)

forj=l:length(xx2)

if((2*xxl(i)-xx2(j)-30<=0)&&(20-xxl(i)<=0))

Z(iJ)=xxl(i)A2+xx2(j)A2;

else

Z(iJ)=0;

end

end

end

figure(l);

surf(xxl,xx2zZ);

axis([05005004500])

title('目標(biāo)函數(shù)曲線圖');

xlabel('xl');

ylabel('x2');

%繪制目標(biāo)函數(shù)等值線圖，并畫(huà)出優(yōu)化路徑

figure(2);

xll=-5:0.5:25;

xl2=-5:0.5:25;

[xxll,xxl2]=meshgrid(xll,xl2);

F=xxll.A2+xxl2.A2;

axis([-525-525]);

contour(xxll,xxl2,F);

title('目標(biāo)函數(shù)等值線')；

xlabel('xl')；

ylabel('x2');

holdon

plot(a,b,'r+-');%繪制優(yōu)化路徑