外点罚函数优化实例

外点罚函数优化实例

一.优化问题

如图1所示，为某一桁架的一部分，杆2距。点30cm处有一支点C。为了固定

桁架，现欲在杆1和2上设置支点A和B,用来连接杆3（可拆卸1已知当桁架

固定时，杆1和2成直角；而且，杆1右边有一段长为20cm的重要部位，不能设

置支点。卸去杆3.收起桁架时，支点A的位置不能高于BC段中点D。求取支点

A.B的位置，使得杆3的长度尽量小，以节省材料。

C■-------

杆2

D,■—A，

/

/B,•-T软

/杆:V,

宴

杆1I/

AL20’0

一XI一

图1桁架结构示意图

二.数学模型

设A.B两点距离。点的长度分别为口和口，而桁架固定时杆1和2成直角。所

以杆3的长度为口。

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

设口取口。因此数学模型为：

极小化目标函数

niin/(X)=X12+x}

约束条件

4(X)=F+20V0

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

三.求解数学模型

(1)外点罚函数法求解

构造外点法罚函数，如下：

22

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

程序流程图如图2所示：

图2外点罚函数法程序流程图

程序步骤：

①选择适当的初始罚因子口、初始点口、收敛精度□和罚因子系数C。在本程序

中分别取口,□,□,c=8e令迭代步数k=0o

②采用牛顿法求无约束问题minX,")的极值点X'(M⑹)。

目标函数等值残

25

图4目标函数等值线图

(2)Matlab优化工具fmincon求解

利用Matlab文件编辑器为目标函数编写M文件(goalfun.m):

functionf=goalfun(x)

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

编写约束条件的M文件(confun.m):

function[c,ceq]=confun(x)

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

20-x(l)];

ceq=[];

编写主函数的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

运行结果：

x=[20,10]

fval=500

同样地，利用Matlab优化工具解得，支点A.B与O的距离分别为20cm、

10cmz杆3的最小长度为Clem。

四.结论分析

Q)罚因子系数c对外点罚函数法的影响

本次程序中,c取值为8,运行步数k=llo若取c=4,则运行步数k=16;取c=16,

则运行步数k=9;取c=32,则运行步数k=8;取c=64r则运行步数k=70由此

可知，罚因子系数c的大小会影响程序的迭代次数k°c的值取得越大，运行步数k

越小，程序收敛速度越快，效率越高。但对于C的其他一些取值，如5.7、9等，会

导致罚函数形态变坏，使迭代出现问题，导致程序运行失败。因此，需选取合适的

罚因子系数c。

(2)外点罚函数法与Matlab优化工具fmincon的比较

通过。pt.m的运行结果可知,fmincon的运行步数k=3,这一运行效率明显比外点

法函数法的效率高。对于相同的收敛精度口和初始点，虽然优化结果相同，但是M

文件WaiDianNiuDun.m中外点罚函数法的运行效率仍有待提高。

附一外点罚函数matlab程序

closeall

clearall

cic

symsxlx2M;%M为罚因子。

m(l)=l;

c=8;%c为递增系数。赋初值。

a(l)=20;

b(l)=20;

f=xl人2+x2人2+M*((20-xl)A2+(2*xl-x2-30)八2);%夕卜点罚函数

f0(l)=500;

%求偏导、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')；

%外点法M迭代循环

fork=l:100

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

%牛顿法求最优值

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)%最优值收敛条

件

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)%罚因子迭代收敛条件

%输出最优点坐标，罚因子迭代次数，最优值

a(k+l)

b(k+l)

k

fO(k+l)

break;

else

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

end

end

%绘制目标函数曲线图

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('目标函数曲线图');

xlabel('xl');

ylabel('x2');

%绘制目标函数等值线图，并画出优化路径

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('目标函数等值线')；

xlabel('xl')；

ylabel('x2');

holdon

plot(a,b,'r+-');%绘制优化路径