外点罚函数优化实例

外点罚函数优化实例

一、优化问题

如图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 点的长度分别为1x 和2x ，而桁架固定时杆1和2成直角。

所以杆3的长度为22213x x l +=。

由图1可知，201≥x 且)30(2

121+≤x x ，即0201≤+-x 且030221≤--x x 。 设T x x X ],[21=，取222123)(x x l X f +==。因此，数学模型为：

极小化目标函数

2221)(min x x X f +=

约束条件

020)(11≤+-=x X g

0302)(212≤--=x x X g

三、求解数学模型

(1)外点罚函数法求解

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

})]0,302[m ax ()]0,20{[(m ax (),(22121)(2221)(--+-++=x x x M x x M X k k φ

程序流程图如图2所示：

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

程序步骤：

①选择适当的初始罚因子)0(M 、初始点)0(X 、收敛精度ε和罚因子系数c 。在本程序中分别取1)0(=M ，]20,20[)0(=X ，610-=ε，c=8。令迭代步数k=0。 ②采用牛顿法求无约束问题),(m in )(k M X φ的极值点)()(*k M X 。

③检验迭代终止准则，若满足

ε≤--||)()(||)1()(*k k M X M X 及 ε≤---|)]

([)]([)]([|)1(*)1(*)(*k k k M X f M X f M X f 则停止迭代计算，输出最优点)()(**k M X X =；否则，转入步骤④。

④取)()1(k k cM M =+，)()(*)0(k M X X =，k=k+1，转入步骤②继续迭代。 具体程序请看附一。

运行结果：

X *=[20.0000，10.0000]

f(X *)=500.0000

因此，A 、B 两支点与O 的距离分别为20cm 、10cm ，杆3的最小长度为510cm 。 目标函数曲线图与目标函数等值线图分别如图3和图4所示。优化路径如图4所示。

图3 目标函数曲线图

图4 目标函数等值线图

(2)Matlab优化工具fmincon求解

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

function f=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)：

close all

clear all

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，10cm。

杆3的最小长度为5

四、结论分析

(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程序

close all

clear all

clc

syms x1 x2 M; %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迭代循环

for k=1:100

x1=a(k);x2=b(k);M=m(k);

%牛顿法求最优值

for n=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=[x1 x2]'-inv([f11 f12;f21 f22])*[f1 f2]';

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)))<=1 e-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;

for i=1:length(xx1)

for j=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([0 50 0 50 0 4500])

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([-5 25 -5 25]);

contour(xx11,xx12,F);

title('目标函数等值线'); xlabel('x1');

ylabel('x2');

hold on

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