机械优化设计研究生大作业

第一题

1.1题目

,222

求函数f(X)=x,-2xlx2-2x,x2+3x1+4x2+4.5X,-4X2+5的极小值，初始点为

X@=[-2,2]T,误差£不大于0.001。

注：此问题为无约束非线性规划问题的求解。

1.2建立数学模型

FindX],x2

42+22+

minf(X)=X1-2XIX2-2XIX23XI+4X24.5XI-4X2+5

初始点X(0)=[-2,2]T,£WO.001

1-3运行结果

X1x2f

迭

数\

1-0.50000000000000-0.583333333333338.33333333333333

2-0.566062176165800.370682210708123.00457834092293

3-0.652874320923100.398982834197582.97857158805445

4-0.650848939158930.401358609903932.97849714375511

5-0.650839107315690.401356625245232.97849714338108

通过牛顿法迭代5次可得出结果，当xl=-0.65083910731569,

x2=0.40135662524523时，目标函数最优值fmin=2.97849714338108,且满足

JM0.01。

1.4迭代曲线

他)函数迭代值

X1变化图

-0.48

■0.5

Q

。

姻Q52

W54

枳Q56

I58

0.6

£62

64

Q

。

■0.66

11.522.533.544.5

迭代次数

X2变化图

0.6

0.4

a2

出0

卫

根Q2

xe

■0.4

-0.6

-0.8

34

迭代次数

1.5检驷结果

用Matlab自带优化程序检验

程序为：为x0=[-2,2];

[x,fval]=fminsearch(，x(l)*4-5*x(l)*2*x(2)-2*x(1)*x(2)+4*x(l)*2+6*x(2)

"2+4.5*x(l)-4*x(2)+5),x0)

x=-0.650866586874660.40137142333985

fval=2.97849714628600

经检验用牛顿法进行迭代优化结果是正确的，优化结果达到精度要求，eW

0.OOE

1.6讨论

(1)由以上迭代曲线可知，牛顿法迭代收敛速度很快，本优化经过迭代3

次后目标函数值趋于平稳。也可采用黄金分割法，变尺度法等其他方法优化。由

于本题比较简单，不必采用变尺度法来优化。

(2)采用Matlab编程解决了求导和计算海森阵比较复杂的难题，编程简单

方便。

1.7Matlab源程序

functionZY32

formatlong

symsxlx2%定义符号变量xl,x2

f=xl"4-2*xl〜2*x2-5*xl*x2+6*xl-2+7*x2~2+4.5*xl-4*x2+5%定义函数f

df=[diff(f,xl);diff(f,x2)]

%diff(f,xl)用于对函数f中变量xl求偏导

%diff(f,x2)用于对函数f中变量x2求偏导

fl=diff(f,xl,2);%diff(f,xl,2)用于对函数f中变量xl求而二次偏导

f2=diff(diff(f,xl),x2);

%diff(diff(f,xl),x2)用于对函数f中变量下xl,x2求偏导

f3=diff(diff(f,x2),xl);

f4=diff(f,x2,2);%diff(f,x2,2)用于对函数f中变量x2求而二次偏导

ddf=[fl,f2;f3,f4];%求函数f的海森阵

xl=-2;x2=2;

td=eval(df);肮十算梯度初值

hs=eval(ddf);%计算海森阵初值

%eval命令用于将符号变量转化为数值变量

i=0;

eps=0.1

whileeps>0.01;

i=i+l;

eps=sqrt(td(l).2+td(2)2);%求梯度的模

x=[xl;x2]-inv(hs)*td;肮己录自变量迭代值

xl=x⑴；

x2=x⑵；

xx(i,:)=x';%记录迭代过程中间值

td=eval(df);

hs=eval(ddf);

last(i)=eval(f);断己录迭代过程中间值

end

disp('迭代次数')；

xx=xx

fmin=last

figure(l)

plot(fmin);%画目标函数变化图

gridon;

holdon;

xlabelC迭代次数');

ylabel('f(x)');

title，f(x)函数迭代值');

figure(2)

plot(xx(:,1));

gridon;

xlabel('迭代次数')；

ylabelCxl迭代值');

titleCxl变化图');

figure(3)

plot(xx(:,2));

gridon;

xlabel('迭代次数')；

ylabelCx2迭代值');

titleCx2变化图')；

第二题

2.1数学模型

FindXi,x2

22

minf(X)=2XI+2X2-4X-6X2

s.t.gi(X)=XI+5X2-5^0

2

g2（X）=2X,-X2^0

g3（X）=-xW0

g，1（X）=-x2<0

初始点内）=[0,0.75『,e=0.001

2.2优化方法

此问题可规结为强约束问题，此题用罚函数法对在求解强约束非线性优化问

题求解。用功能函数constr可以解决此类问题，本题中目标函数是非线性的，

约束是线性的。

2.3优化结果

最小目标函数值f=-5.5039,

自变量取值：X=[0.6589,0.8682]T

验证优化结果正确。

2.4Matlab源程序

在commandwindow输入question2,然后回车即行;

Functionquestions2.m

funf=，f=x(l)-2+4*x(2)-2-7*x(l)-5*x(2);%定义目标函数

fung=，g=[x⑴+5*x(2)-5;2*x(1)^2-x(2)];%定义约束条件

funi=[funffung];%组合矩阵

xO=[O,0];%初瞰

options=l;%显示过程

vlb=[00];%下边界x；

vub=[];%无上边界；

[x,options]=constr(funi,xO,options,vlb,vub)

2.5计算结果

迭

代目标值自变量自变量

约束&约束g2约束g3约束g.,

次fXiX1

数

11010-10-1

-1.000-0.000-1.0000-1.0001.0000e-00

21.0000-1

00e-00808

-1.000-5.0000-1.000

301.00000-1.0000

0e-0080

4-1050050-50

-12.75

52.50000.5000012-2.5000-0.5000

00

-8.687

61.25000.750002.3750-1.2500-0.7500

5

-8.687-1.0000

71.25000.75002.3750-1.2500-0.7500

5e-008

-8.687-5.0000

81.25000.75002.3750-1.2500-0.7500

5e-008

-5.936

90.71520.4511-2.02930.5720-0.7152-0.4511

5

-5.516

100.66090.867800.0059-0.6609-0.8678

4

-5.5161.OOOOe

110.66090.86780.0059-0.6609-0.8678

4-008

-5.5165.OOOOe

120.66090.86780.0059-0.6609-0.8678

4-008

-5.5038.5830

130.65890.86820-0.6589-0.8682

9e-006

-5.5031.OOOOe8.6094

140.65890.8682-0.6589-0.8682

9-008e-006

-5.503-5.00008.5930

150.65890.8682-0.6589-0.8682

9e-008e-006

-5.5031.8429

160.65890.86820-0.6589-0.8682

9e-011

-5.5038.2958e2.1864

170.65890.8682-0.6589-0.8682

9-007e-006

-5.5034.7634e-9.526

180.65890.8682-0.6589-0.8682

9-0059e-006

-5.503-8.8818-1.110

190.65890.8682-0.6589-0.8682

9e-0162e-016

-5.490-0.006

200.65650.8682-0.0024-0.6565-0.8682

52

-5.5042.6130

210.65890.8680-0.0013-0.6589-0.8680

4e-004

-5.503-1.0000

220.65890.86820.0000-0.6589-0.8682

9e-004

2.6迭代曲线

2.7讨论

(1)利用Matlab优化工具箱里边的功能函数condtr来优化求解目标函数

f，能很快实现优化，达到目标，方便快捷。

(2)由图可见迭代五次达到最低端，10次后趋于平稳。

(3)优化结果满足精度要求£<0.001o

(4)可见Matlab很容易对目标实现优化，对编程人员要求不高。

只要知道如何使用命令。

第三题

3.1题目

平面钱链四连杆机构，各杆长度为/2,h，/4o已知，4=5/1；主动杆输

入角为夕，从动杆4的输出角为〃；摇杆&在右极位时，杆1、杆2伸直，主动

杆乙在初始位置生，从动杆。在初始位置“°。试设计四连杆机构的各杆长度，

使其输出角”在曲柄从0。转到外+9"时，最佳再现给定函数：

”=〃0+2*(0-%)2/(3万)，且要求最小传动角不小于45°,即245°。

3.2数学模型

3.2.1确定设计变量

机构位置决定于4个杆长与主动杆转角(p,再现角位移的机构与杆件的绝对

长度无关，只决定于相对长度；转角连续变化，故非设计变量。因此计算时取曲

柄长度为单位长度，即/尸1;而其他的杆长按比例取为乙的倍数，机架长乙=5,

而初始角州可以按下列几何关系求得：

+/)2+(/2-/?)

(p＜、=arccos~2——4----

L2(/1+/2)Z4

12/3/4

故仅有4、4为独立变量，是二维最优化设计问题,设其设计变量为:

x=X'=，］

3.2.2确定目标函数

此再现运动的机械运动优化中，目标函数为机构实际运动轨迹与预定运动轨

迹均方根最小，位置取若干离散点，即

/(X)=ZT〃X)-42fmin

本题则取目标函数为：

/(X)=Z"(X)-忆Jfmin(i=1~P)

。取30,则，,=31,式中期望输出角的离散值为：

匕产%+4(。,一。0>/5万

9=Wo兀12P=(Po+E/60i-(0,1,…,30)

实际输出角的计算公式为：

产一名(X)-Bi(0<(p.t<万)上半圈

I//={,„

“万一%(X)+Bi(乃<(p.t<2%)下半圈

由三角关系式得：

222

%(X)=arccos[(f；+x2—^,)/(2/;x2)]

月(X)=arccos比2+〃_/：)/(2而]

=arccos[(^2+24)/(IOr;)]

205

r.t=(『+/4-2l}l4cos0产=(26-lOcos^,)

3.2.3确定约束条件

根据对传动角的约束要求及曲柄与机架处于共线位置时4-x和7min机构的尺

寸的关系：小=arcc十鸟—卜鼠」心

/2/2_“_/\2

341

/min=arccos-—^[/max]=45°以及COS/max>cos[/max]；

乙I)Io

COS/min<COS[/min],得约束条件：

_

gi(X)=-X|—x2~~5/2^x2+(Z4+1\)~<0

2

g2(X)=X]~+x2~—42.x{x2—(/4—/])W0

满足曲柄存在条件，得约束条件：

g3(X)=—X]+&+/|)W0

^4(X)=x1-x2-(/4-/,)<0

gs(X)=—X]+x2—(/4—l])<0

g6(X)=lf«0

^7(X)=l-x2<0

取4=1时，乙=5,进一步分析可行域可知.(X)和g2(X)(两椭圆方程)

起约束作用，其他约束均在其外不起作用。

3.2.4数学模型

2

findX=[^)，%27G7?

30

min/(X)=ZW(X)-仁12

/=0

22

s.t.^,(X)=—x,—x2—V2X,A:2+36<0

2

g2(X)=X：+x2—V2X]X2—16<0

3.3优化方法

用Matlab自带函数fmincon来优化求解，

3.4优化结果

最小目标函数值f=0.0133

自变量X*=[4.4846,1.9224]T

约束条件值g=[0.0000，-4.3848]T

优化结果达到要求.

3.5Matlab程序

由于某些符号在Matlab中无法表示，符号有如下变动，其中inO设为牝，

in设为9，outO设为心,out设为匕,，a设为%(X),b设为4(X),r设为公,

f设为所求函数/(X)。

3.5.1目标函数/(X)

将/(X)程序命名为cpp3.m

functionf=cpp3(x)

L1=1;L4=5;/定义杆1,4长度

L2=x(l);L3=x(2);%定义杆2,3为变量

formatlong;

%初始输入角

inp=acos(((Ll+L2)."2+(L4.'2-L3.-2))/(2*(Ll+L2)*L4));

%初始输出角

outp=acos(((Ll+L2).*2-(L4.-2+L3."2))/(2*L3*L4));

pianc=0;%给定初值。

fori=l:30

in{i}=sym(inp+i*pi/60);%在不同位置的输入角

out{i}=sym(outp+2*(in{i}-inp)."2/(3*pi));%输出角

r{i}=sym(sqrt(26-10*cos(eval(in{i}))));%符号值

a[i}=sym(acos(((eval(r{i}))."2+x(2)."2-x(1).2)/(2*eval(r{i})*x(2

))))；

b{i}=sym(acos(((eval(r{i}))."2+24)/(10*(eval(r{i})))));

A{i}=eval(a{i});%将符号转换成为数值

B{i}=eval(b{i});

NI{i}=eval(in{i});

TUO{i}=eval(out{i});

if0<NI{i}<=pi

real{i}=pi-A{i}-B{i};

elseifpi<NI{i}<=2*pi

rl{i}=pi-A{i}+B{i};

else

dispC超出范围')

end

piancha{i}=(rl{i}-TU0{i})."2;%计算偏差

pinc=piancha{i}+pinc;

end

pine%用于显示优化过程中的目标函数值

%定义约束条件gl,g2;

3.5.2约束函数值

function[g,z]=yueshu3(x)

g=[36-x(l)."2-x(2)."2-1.4142*x(l)*x(2);

x(1).~2+x(2)「2T6T.4142*x(1)*x(2)]

z=[