外点罚函数优化实例Word文档下载推荐.docx
《外点罚函数优化实例Word文档下载推荐.docx》由会员分享，可在线阅读，更多相关《外点罚函数优化实例Word文档下载推荐.docx（8页珍藏版）》请在冰豆网上搜索。
三、求解数学模型
（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
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,'
fx1x2=diff（fx1,'
fx2x1=diff（fx2,'
fx2x2=diff（fx2,'
%外点法M迭代循环
fork=1:
100
x1=a（k）;
x2=b（k）;
M=m（k）;
%牛顿法求最优值
forn=1:
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
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）
else
m（k+1）=c*m（k）;
end
%绘制目标函数曲线图
xx1=0:
0.5:
50;
xx2=0:
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;
Z（i,j）=0;
figure
（1）;
surf（xx1,xx2,Z）;
axis（[05005004500]）
title（'
目标函数曲线图'
xlabel（'
ylabel（'
%绘制目标函数等值线图，并画出优化路径
figure
（2）;
x11=-5:
25;
x12=-5:
[xx11,xx12]=meshgrid（x11,x12）;
F=xx11.^2+xx12.^2;
axis（[-525-525]）;
contour（xx11,xx12,F）;
目标函数等值线'
holdon
plot（a,b,'
r+-'
%绘制优化路径
THANKS!
!
致力为企业和个人提供合同协议，策划案计划书，学习课件等等
打造全网一站式需求
欢迎您的下载，资料仅供参考
升级会员 