外点罚函数优化实例.docx
《外点罚函数优化实例.docx》由会员分享，可在线阅读，更多相关《外点罚函数优化实例.docx（9页珍藏版）》请在冰豆网上搜索。
外点罚函数优化实例
外点罚函数优化实例
一、优化问题
如图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、
图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+-'）;%绘制优化路径
升级会员 