外点罚函数优化实例

1、外点罚函数优化实例一、优化问题如图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、外点法罚函数，如下： 程序流程图如图2所示：给定、c、k=0i=0求与Hessian矩阵输出和YNi=i+1k=k+1YN结束牛顿法求的极值点图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、数曲线图与目标函数等值线图分别如图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 allclear allclcx0=20,20;lb=;ub=;options=optimset(

4、LargeScale,off,display,iter,tolx,1e-6);x,fval,exitflag,output=fmincon(goalfun,x0,confun,options);xfval运行结果：x=20，10fval=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的大小会影

5、响程序的迭代次数k。c的值取得越大，运行步数k越小，程序收敛速度越快，效率越高。但对于c的其他一些取值，如5、7、9等，会导致罚函数形态变坏，使迭代出现问题，导致程序运行失败。因此，需选取合适的罚因子系数c。 (2)外点罚函数法与Matlab优化工具fmincon的比较通过opt.m的运行结果可知，fmincon的运行步数k=3，这一运行效率明显比外点法函数法的效率高。对于相同的收敛精度和初始点，虽然优化结果相同，但是M文件WaiDianNiuDun.m中外点罚函数法的运行效率仍有待提高。附一 外点罚函数matlab程序close allclear allclcsyms x1 x2 M; %M

6、为罚因子。m(1)=1; c=8; %c为递增系数。赋初值。a(1)=20;b(1)=20; f=x12+x22+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); %

7、求解梯度值和Hessian矩阵 f2=subs(fx2); f11=subs(fx1x1); f12=subs(fx1x2); f21=subs(fx2x1); f22=subs(fx2x2); if(double(sqrt(f12+f22)=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 endif(double(sqrt(a(k+1)-a(k)2

8、+(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);endend% 绘制目标函数曲线图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 endendfigure(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=mes