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

第一题1.1 题目

求函数f(X)=x

14-2x

1

2x

2

-2x

1

x

2

+3x

1

2+4x

2

2+4.5x

1

-4x

2

+5的极小值，初始点为

X(0)=[-2，2]T，误差ε不大于0.001。

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

1.2 建立数学模型

Find x

1 , x

2

min f(X)=x

14-2x

1

2x

2

-2x

1

x

2

+3x

1

2+4x

2

2+4.5x

1

-4x

2

+5

初始点 X(0)=[-2，2]T, ε≤0.001 1.3 运行结果

迭

代值

迭代次

数

x1 x2 f

1 -0.500 -0.58333333333333 8.33333333333333

2 -0.566 0.372 3.293

3 -0.65287432092310 0.39898283419758 2.97857158805445

4 -0.653 0.43 2.97849714375511

5 -0.659 0.43 2.97849714338108

通过牛顿法迭代5次可得出结果，当x1=-0.659，x2=0.43时，目标函数最优值fmin=2.97849714338108，且满足0.01

ξ≤。

1.4 迭代曲线

1.5 检验结果

用Matlab自带优化程序检验

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

[x,fval]=fminsearch('x(1)^4-5*x(1)^2*x(2)-2*x(1)*x(2)+4*x(1)^2+6*x(2) ^2+4.5*x(1)-4*x(2)+5',x0)

x =-0.656 0.45

fval =2.97849714628600

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

0.001。

1.6 讨论

（1）由以上迭代曲线可知，牛顿法迭代收敛速度很快，本优化经过迭代3次后目标函数值趋于平稳。也可采用黄金分割法，变尺度法等其他方法优化。由于本题比较简单，不必采用变尺度法来优化。

（2）采用Matlab编程解决了求导和计算海森阵比较复杂的难题，编程简单方便。

1.7 Matlab源程序

function ZY32

format long

syms x1 x2 %定义符号变量x1,x2

f=x1^4-2*x1^2*x2-5*x1*x2+6*x1^2+7*x2^2+4.5*x1-4*x2+5 %定义函数f df=[diff(f,x1);diff(f,x2)]

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

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

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

%diff(diff(f,x1),x2)用于对函数f中变量下x1,x2求偏导 f3=diff(diff(f,x2),x1);

f4=diff(f,x2,2); %diff(f,x2,2)用于对函数f中变量x2求而二次偏导 ddf=[f1,f2;f3,f4]; %求函数f的海森阵

x1=-2;x2=2;

td=eval(df); %计算梯度初值

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

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

i=0;

eps=0.1

while eps>0.01;

i=i+1;

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

x=[x1;x2]-inv(hs)*td; %记录自变量迭代值

x1=x(1);

x2=x(2);

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

td=eval(df);

hs=eval(ddf);

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

end

disp('迭代次数');

i

xx=xx

fmin=last

figure(1)

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

grid on;

hold on;

xlabel('迭代次数');

ylabel('f(x)');

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

figure(2)

plot(xx(:,1));

grid on;

xlabel('迭代次数');

ylabel('x1迭代值');

title('x1变化图');

figure(3)

plot(xx(:,2));

grid on;

xlabel('迭代次数');

ylabel('x2迭代值');

title('x2变化图');

第二题2.1 数学模型

Find x

1 , x

2

min f (X)=2x

12+2x

2

2-4x

1

-6x

2

s.t. g

1(X)=x

1

+5x

2

-5≤0

g 2(X)=2x

1

2-x

2

≤0

g 3(X)=-x

1

≤0

g 4(X)=-x

2

≤0

初始点X(0)=[0,0.75]T, ε=0.001

2.2 优化方法

此问题可规结为强约束问题，此题用罚函数法对在求解强约束非线性优化问题求解。用功能函数constr可以解决此类问题，本题中目标函数是非线性的，约束是线性的。

2.3 优化结果

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

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

验证优化结果正确。

2.4 Matlab 源程序

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

Function questions2.m

funf='f=x(1)^2+4*x(2)^2-7*x(1)-5*x(2) ; ％定义目标函数

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

funi=[funf fung]; ％组合矩阵

x0=[0,0]; %初值；

options=1; ％显示过程 vlb=[0 0]; ％下边界x； vub=[ ]; %无上边界；

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

2.5 计算结果

18 -5.503

9

0.6589 0.8682

4.7634e

-005

-9.526

9e-006

-0.6589 -0.8682

19 -5.503

9

0.6589 0.8682

-8.8818

e-016

-1.110

2e-016

-0.6589 -0.8682

20 -5.490

5

0.6565 0.8682 -0.0024

-0.006

2

-0.6565 -0.8682

21 -5.504

4

0.6589 0.8680 -0.0013

2.6130

e-004

-0.6589 -0.8680

22 -5.503

9

0.6589 0.8682

-1.0000

e-004

0.0000 -0.6589 -0.8682

2.6 迭代曲线

2.7 讨论

（1）利用Matlab优化工具箱里边的功能函数condtr来优化求解目标函数f，能很快实现优化，达到目标，方便快捷。

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

（3）优化结果满足精度要求ε＜0.001。

（4）可见Matlab很容易对目标实现优化，对编程人员要求不高。只要知道如何使用命令。

第三题

3.1 题目

平面铰链四连杆机构，各杆长度为1l ，2l ，3l ，4l 。已知145l l =；主动杆1l 输入角为?，从动杆3l 的输出角为ψ；摇杆3l 在右极位时，杆1、杆2伸直，主动杆1l 在初始位置0?，从动杆3l 在初始位置0ψ。试设计四连杆机构的各杆长度，使其输出角ψ在曲柄从0?转到?+900?时，最佳再现给定函数：

2002*()/(3)ψψ??π=+-，且要求最小传动角不小于45ｏ，即≥45ｏ。

3.2 数学模型

3.2.1 确定设计变量

机构位置决定于4个杆长与主动杆转角?，再现角位移的机构与杆件的绝对长度无关，只决定于相对长度；转角连续变化，故非设计变量。因此计算时取曲柄长度为单位长度，即1l =1；而其他的杆长按比例取为1l 的倍数，机架长45l =，而初始角0?可以按下列几何关系求得：

???

???+-++=42123242210)(2)()(arccos l l l l l l l ?

??

?

???+-+=4

3242322102)()(arccos l l l l l l ψ

故仅有2l 、3l 为独立变量，是二维最优化设计问题，设其设计变量为：

??

????=??????=3221l l x x X

3.2.2 确定目标函数

此再现运动的机械运动优化中，目标函数为机构实际运动轨迹与预定运动轨迹均方根最小，位置取若干离散点，即

[]min )()(2

→-=∑i i i f X f w X f

本题则取目标函数为：

[]min )()(2

→-=∑i e i X X f ψψ )~1(P i =

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

π??ψψ5/)(4200-+=i ei

60/2/00π?π??i P i i +=+= )30,,1,0(Λ=i

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

i i i i ei X X βαπβαπψ+---=)()({

下半圈

上半圈

)2()0(π?ππ?≤<≤

)]2/()arccos[()(221222x r x x r X i i i -+=α )]2/()arccos[()(421242l r l l r X i i i -+=β

)]10/()24arccos[(2

i i r r +=

5.05.0412

421)cos 1026()cos 2(i i i l l l l r ??-=-+=

3.2.3 确定约束条件

根据对传动角的约束要求及曲柄与机架处于共线位置时

m ax

γ和min γ机构的尺

寸的关系：[]ο

1352)(arccos max 3

22142322max

=≤??

????+-+=γγl l l l l l

[]ο

45

2)(arccos max 3

22142322min

=≥??

????--+=γγl l l l l l 以及

]

cos[cos max max γγ≥；

]cos[cos min min γγ≤，得约束条件：

0)(2)(214212

22

11≤++---=l l x x x x X g

0)(2)(214212

2212≤---+=l l x x x x X g

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

0)()(14213≤++--=l l x x X g

0)()(14214≤---=l l x x X g

0)()(14215≤--+-=l l x x X g

01)(16≤-=x X g 01)(27≤-=x X g

取5141==l l 时，，进一步分析可行域可知和)(1X g )(2X g （两椭圆方程）起约束作用，其他约束均在其外不起作用。 3.2.4 数学模型

find []221,R x x X T

∈=

min ∑=-=30

2])([)(i ei i X X f ψψ

..t s 0362)(212

2211≤+---=x x x x X g 0162)(212

22

12≤--+=x x x x X g

3.3 优化方法

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

3.4 优化结果

最小目标函数值f=0.0133

自变量X * =[4.4846 ， 1.9224]T 约束条件值g=[0.0000， -4.3848]T 优化结果达到要求.

3.5 Matlab 程序

由于某些符号在Matlab 中无法表示，符号有如下变动，其中in0设为0?，in 设为i ?， out0设为i ψ，out 设为ei ψ，a 设为)(X i α，b 设为)(X i β，r 设为i r ，f 设为所求函数)(X f 。 3.5.1 目标函数)(X f

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

function f=cpp3(x)

L1=1;L4=5; %定义杆1，4长度 L2=x(1);L3=x(2); %定义杆2，3为变量

format long;

%初始输入角

inp=acos(((L1+L2).^2+(L4.^2-L3.^2))/(2*(L1+L2)*L4));

%初始输出角

outp=acos(((L1+L2).^2-(L4.^2+L3.^2))/(2*L3*L4));

pianc=0; %给定初值。

for i=1: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});

if 0

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

elseif pi

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

else

disp('超出范围')

end

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

pinc=piancha{i}+pinc;

end

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

%定义约束条件g1,g2;

3.5.2 约束函数值

function [g,z]=yueshu3(x)

g=[36-x(1).^2-x(2).^2-1.4142*x(1)*x(2);

x(1).^2+x(2).^2-16-1.4142*x(1)*x(2)]

z=[];

x %输出自变量值；

3.5.3 求解目标函数

命名文件CPP33.m.

x=[5 3] ; % x为初始值

options = optimset ('LargeScale' , 'off');

[x ,f] =fmincon ( @ cpp3, x ,[], [], [], [], [], [], @yueshu3, options) 在command window中输入cpp33.m, 回车即可。

3.6 计算结果

数

1 0.313

2 5

3 -19.2130 -3.2130

2 0.1630 5.0972 2.6247 -15.7902 -2.0495

3 0.738

4 5.0472 1.4130 -1.5563 1.3849

4 0.0211 5.0722 2.0188 -8.2842 -0.6785

5 0.0204 5.017

6 2.0210 -7.6010 1.0804

6 0.0173 4.7873 1.9538 -3.9638 -2.4920

7 0.0147 4.5284 1.8818 -0.0986 -4.0034

8 0.0140 4.5109 1.8922 -0.0002 -4.1423

9 0.0140 4.5109 1.8922 -0.0001 -4.1423

10 0.0140 4.5109 1.8922 -0.0002 -4.1423

11 0.0133 4.4868 1.9199 -0.0004 -4.3646

12 0.0133 4.4846 1.9224 -0.0000 -4.3848

13 0.0133 4.4846 1.9224 0.0001 -4.3848

14 0.0133 4.4846 1.9224 -0.0000 -4.3848

3.7 迭代曲线

3.8 讨论

（1）利用Matlab的优化工具箱能够很方便的实现目标优化。。

（2）由图可看出迭代4次后，目标函数基本上达到最小值，并且满足精度要求。

（3）因为用到了for循环连续重复调用函数所以计算速度较慢。这是Matlab 的缺陷。但仍是解决问题的好方法。