外点罚函数优化实例.doc
《外点罚函数优化实例.doc》由会员分享,可在线阅读,更多相关《外点罚函数优化实例.doc(8页珍藏版)》请在冰豆网上搜索。
外点罚函数优化实例
一、优化问题
如图1所示,为某一桁架的一部分,杆2距O点30cm处有一支点C。
为了固定桁架,现欲在杆1和2上设置支点A和B,用来连接杆3(可拆卸)。
已知当桁架固定时,杆1和2成直角;而且,杆1右边有一段长为20cm的重要部位,不能设置支点。
卸去杆3、收起桁架时,支点A的位置不能高于BC段中点D。
求取支点A、B的位置,使得杆3的长度尽量小,以节省材料。
图1桁架结构示意图
二、数学模型
设A、B两点距离O点的长度分别为和,而桁架固定时杆1和2成直角。
所以杆3的长度为。
由图1可知,且,即且。
设,取。
因此,数学模型为:
极小化目标函数
约束条件
三、求解数学模型
(1)外点罚函数法求解
构造外点法罚函数,如下:
程序流程图如图2所示:
给定、、c、
k=0
i=0
求与Hessian矩阵
输出和
Y
N
i=i+1
k=k+1
Y
N
结束
牛顿法求的极
值点
图2外点罚函数法程序流程图
程序步骤:
①选择适当的初始罚因子、初始点、收敛精度和罚因子系数c。
在本程序中分别取,,,c=8。
令迭代步数k=0。
②采用牛顿法求无约束问题的极值点。
③检验迭代终止准则,若满足
及
则停止迭代计算,输出最优点;否则,转入步骤④。
④取,,k=k+1,转入步骤②继续迭代。
具体程序请看附一。
运行结果:
X*=[20.0000,10.0000]
f(X*)=500.0000
因此,A、B两支点与O的距离分别为20cm、10cm,杆3的最小长度为cm。
目标函数曲线图与目标函数等值线图分别如图3和图4所示。
优化路径如图4所示。
图3目标函数曲线图
图4目标函数等值线图
(2)Matlab优化工具fmincon求解
利用Matlab文件编辑器为目标函数编写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=[];
编写主函数的M文件(opt.m):
closeall
clearall
clc
x0=[20,20];
lb=[];
ub=[];
options=optimset('LargeScale','off','display','iter','tolx',1e-6);
[x,fval,exitflag,output]=fmincon('goalfun',x0,[],[],[],[],[],[],'confun',options);
x
fval
运行结果:
x=[20,10]
fval=500
同样地,利用Matlab优化工具解得,支点A、B与O的距离分别为20cm、10cm,杆3的最小长度为cm。
四、结论分析
(1)罚因子系数c对外点罚函数法的影响
本次程序中,c取值为8,运行步数k=11。
若取c=4,则运行步数k=16;取c=16,则运行步数k=9;取c=32,则运行步数k=8;取c=64,则运行步数k=7。
由此可知,罚因子系数c的大小会影响程序的迭代次数k。
c的值取得越大,运行步数k越小,程序收敛速度越快,效率越高。
但对于c的其他一些取值,如5、7、9等,会导致罚函数形态变坏,使迭代出现问题,导致程序运行失败。
因此,需选取合适的罚因子系数c。
(2)外点罚函数法与Matlab优化工具fmincon的比较
通过opt.m的运行结果可知,fmincon的运行步数k=3,这一运行效率明显比外点法函数法的效率高。
对于相同的收敛精度和初始点,虽然优化结果相同,但是M文件WaiDianNiuDun.m中外点罚函数法的运行效率仍有待提高。
附一外点罚函数matlab程序
closeall
clearall
clc
symsx1x2M;%M为罚因子。
m
(1)=1;
c=8;%c为递增系数。
赋初值。
a
(1)=20;
b
(1)=20;
f=x1^2+x2^2+M*((20-x1)^2+(2*x1-x2-30)^2);%外点罚函数
f0
(1)=500;
%求偏导、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');
%外点法M迭代循环
fork=1:
100
x1=a(k);x2=b(k);M=m(k);
%牛顿法求最优值
forn=1:
100
f1=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)%最优值收敛条件
a(k+1)=double(x1);b(k+1)=double(x2);f0(k+1)=double(subs(f));
break;
else
X=[x1x2]'-inv([f11f12;f21f22])*[f1f2]';
x1=X(1,1);x2=X(2,1);
end
end
if(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)%罚因子迭代收敛条件
%输出最优点坐标,罚因子迭代次数,最优值
a(k+1)
b(k+1)
k
f0(k+1)
break;
else
m(k+1)=c*m(k);
end
end
%绘制目标函数曲线图
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;
else
Z(i,j)=0;
end
end
end
figure
(1);
surf(xx1,xx2,Z);
axis([05005004500])
title('目标函数曲线图');
xlabel('x1');
ylabel('x2');
%绘制目标函数等值线图,并画出优化路径
figure
(2);
x11=-5:
0.5:
25;
x12=-5:
0.5:
25;
[xx11,xx12]=meshgrid(x11,x12);
F=xx11.^2+xx12.^2;
axis([-525-525]);
contour(xx11,xx12,F);
title('目标函数等值线');
xlabel('x1');
ylabel('x2');
holdon
plot(a,b,'r+-');%绘制优化路径