外点罚函数优化实例文档格式.docx

外点罚函数优化实例文档格式.docx

《外点罚函数优化实例文档格式.docx》由会员分享，可在线阅读，更多相关《外点罚函数优化实例文档格式.docx（9页珍藏版）》请在冰豆网上搜索。

因此，数学模型为：

极小化目标函数

约束条件

三、求解数学模型

（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的最小长度为

四、结论分析

（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+-'

%绘制优化路径