合肥工业大学机械优化设计大作业

给定X(0)、M(0)、c、? k=0 i=0 k=k+1 求??(X)与Hessian矩阵H(X(i)) (i)i=i+1 X(0)?X*(M(k)) M(k?1)?cM(k) X(i?1)?X(i)?[H(X(i))]?1??(X(i))||??(X(i))||?? Y N 牛顿法求min?(X,M(k))的极 X(M*(k))?X(i)值点X(M*(k)) ||X*(M(k))?X(M(k?1))||??N | f[X*(M(k))]?f[X*(M(k?1))]|??*(k?1)f[X(M)]Y 输出X和*f(X*) 结束

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

程序步骤：

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

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

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

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

X*=[25.0000，10.0000] f(X*)=475.0000

因此，A、B两支点与O的距离分别为25cm、10cm，杆3的最小长度为475cm。

目标函数曲线图与目标函数等值线图分别如图3和图4所示。优化路径如图4所示。

图3 目标函数曲线图

图4 目标函数等值线图

4、结论分析

罚因子系数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。 附程序：

close all clear all clc

syms x1 x2 M; %M为罚因子。 m(1)=1;

c=8; %c为递增系数。赋初值。 a(1)=25;

b(1)=20;

f=x1^2+x2^2-x1*x2+M*((25-x1)^2+(2*x1-x2-40)^2); %外点罚函数 f0(1)=200;

%求偏导、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-f1*f2))<=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-(a(k+1)-a(k))*

(b(k+1)-b(k))))<=1e-6)&&(double(abs((f0(k+1)-f0(k))/f0(k)))<=1e-6) %罚因子迭代收敛条件