matlab（四连杆优化设计）

机械优化设计在matlab中的应用

东南大学机械工程学院**

一 优化设计目的：

在生活和工作中，人们对于同一个问题往往会提出多个解决方案，并通过各方面的论证从中提取最佳方案。最优化方法就是专门研究如何从多个方案中科学合理地提取出最佳方案的科学。由于优化问题无所不在，目前最优化方法的应用和研究已经深入到了生产和科研的各个领域，如土木工程、机械工程、化学工程、运输调度、生产控制、经济规划、经济管理等，并取得了显著的经济效益和社会效益。

二 优化设计步骤：

1.机械优化设计的全过程一般可以分为如下几个步骤：

1)建立优化设计的数学模型；

2)选择适当的优化方法；

3)编写计算机程序；

4)准备必要的初始数据并伤及计算；

5)对计算机求得的结果进行必要的分析。

其中建立优化设计数学模型是首要的和关键的一步，它是取得正确结果的前提。优化方法的选取取决于数学模型的特点，例如优化问题规模的大小，目标函数和约束函数的性态以及计算精度等。在比较各种可供选用的优化方法时，需要考虑的一个重要因素是计算机执行这些程序所花费的时间和费用，也即计算效率。

2.建立数学模型的基本原则与步骤

word/media/image1.gif 设计变量的确定；

设计变量是指在优化设计的过程中，不断进行修改，调整，一直处于变化的参数称为设计变量。设计变量的全体实际上是一组变量，可用一个列向量表示：

x=742d716efc83f24c8b40823861cf9442.png。

word/media/image2.gif 目标函数的建立；

选择目标函数是整个优化设计过程中最重要的决策之一。当对某以设计性能有特定的要求，而这个要求有很难满足时，则针对这一性能进行优化会得到满意的效果。目标函数是设计变量的函数，是一项设计所追求的指标的数学反映，因此它能够用来评价设计的优劣。

目标函数的一般表达式为：

f(x)=15e98e44a0804018a4f3c1199773c223.png，要根据实际的设计要求来设计目标函数。

word/media/image3.gif 约束条件的确定。

一个可行性设计必须满足某些设计限制条件，这些限制条件称为约束条件，简称约束。

由若干个约束条件构成目标函数的可行域，而可行域内的所有设计点都是满足设计要求的，一般情况下，其设计可行域可表示为

在可行域中，任意设计点满足全部约束条件，称为可行解，但不是最优解，而优化设 计就是要求出目标函数在可行域的最优解。

三 实例分析 （机械优化设计P241页例8-5）

设计一曲柄摇杆机构如图，要求：

曲柄dc488a398c287323d6a37d25ba87d147.png从44a4aafcdaf3a992ba992eee0b96ff8b.png

80c29e0435e64860d40f03f3e70cf4ce.png且已知429b9d8c1937ae0b57568ec6471a6de4.png范围内变化。

分析：

1） 设计变量的确定

决定机构尺寸的各杆长度，以及当摇杆按已知运动规律开始运行时，曲柄所载的位置角bf2bc469c46fe7bbab66107f42f74731.png应列为设计变量，即：

X=c1db2d1529bfc9d4c2f82691fc8364ff.png=a93a9fa2989c15f61b9a0fbb719d0596.png

考虑到机构的杆长按比例变化时，不会改变其运动规律，因此在计算时常取f7d04059b94df4461ca5b36815343b23.png

，而其他杆长则按比例取为dc488a398c287323d6a37d25ba87d147.png的倍数。若取曲柄的初始位置角为极位角，则bf2bc469c46fe7bbab66107f42f74731.png及相应摇杆392b6d73c62722c51eff4b8f8b38072a.png位置角e6cbb8dd4a432f373957dab75c941e07.png均为杆长的函数，几何图形关系如右图,其关系式为：

bf2bc469c46fe7bbab66107f42f74731.png=arcos[780cceebd88e64b9e9c821b3c530819e.png] （1）

e6cbb8dd4a432f373957dab75c941e07.png=arcos[a05e2b1b7f5e7767d616187778af40b0.png] （2）

将43a5e7187717bebef99995527c24caff.png的长度代入上式（1），（2）得到：

bf2bc469c46fe7bbab66107f42f74731.png=arcos[f7904ea923f8e6050c5da182a9deb35f.png]

e6cbb8dd4a432f373957dab75c941e07.png=arcos[1979d28dc650cf46b34b9d016b02a377.png]

因此，只有91b936f1d3197385d97675e8a00d8cf8.png为独立变量，设计变量减少，故最后的设计变量为：

X=1fde3d3b7214cb4e281783b5005ef0e3.png=bf9a80fb25bbd69834b3c00d19ac5476.png

2) 目标函数的建立

目标函数可根据已知的运动规律与机构实际运动规律之间的偏差最小为指标来建立，即

f(x)=de9119fb96a4aac34ced465cb1ebb890.png → min

式中 83a8e63da82ec68a9712b38d8d507b9c.png—期望输出角，6318bffb479c44c353dac44fad86c70a.png；

m — 输入角等分数；

1915e87837adb2afcde76f9bffd7c821.png—实际输出角，由下图得：

a) 0≤9f673ae7eda2453cbf9896bad9f4fa92.png<db8a3cb40ed3ef2387211e8c2fa2c3ec.png b)ec173a73a4b77bfd6cb1bff32fd4ed1f.png

3050d7342a96ff8da40f18d1a1ccc4e7.png （3）

式中

810db9a979ffe0c03a833309f8fbe4ce.png

edf3153eff750ecbffb16bf06b108ad9.png

6176b70f13c871da510ca83e1b11f501.png

将输入角分成30等分，并用近似公式计算，可得目标函数的表达式：

f(x)=2db4b56ea15358f57928172480bd3bad.png

由题意知，传动角的变化范围是5dca31e4fe609e2f38bbd1b9dce573f5.png，则上式中变量的最后形式可以写成：

d29ac7cdaab7066e7c5d8d9af43640fe.png

ce9cc0702ee2d372ec6f03d5e4239889.png

1c57f3a37c709cecd3e51079c2d7ea9c.png5a701b17cf01d3abc6124e0551e0c51e.png

9ae52ffae4bcac38335b30a4fa1f1c69.png

34dcd7b4c7c4baa1529f92e029d435d1.png

将8ce2d72956a70b4d818622a223a1dbdb.png带入（4） （5）得：

03425042ea8a2ecc27fc9377dcb42e8f.png

83a8e63da82ec68a9712b38d8d507b9c.png为当a66d02b02867a74f0896c93eeaebc27e.png时的理想输出角，其值在题目中已经给出：

d4bbcc6cc0957c0e99e3e77de97f5910.png

3）约束条件的确定

word/media/image1.gif曲柄摇杆机构应满足曲柄存在条件，可得

69109eccdbf08310c60491dc4af3fe26.png

298b0e4cbf3aad627877f853301e0a2c.png

83307248d663f9c9c4610e5208e96f6d.png

8f932f03b07bb4fb616bcc82772f5cca.png

34bf8924defc01dcedad49d0064d57f8.png

ac665d4c101405e2c31b1e8577a4d4d1.png

word/media/image2.gif曲柄摇杆机构的传动角6cb2abde9ccca7a7a5eb292698839476.png可得

a9b4897c3eeffcf8f580cd1d983f0dd0.png

1aba39f7ab20cea607a3ced55e8ebc7f.png

把约束条件简化(6bb2a547d1b26cd5df693c768298554e.png eecaf73e13c3f1b4eabfc2159329bd17.png=5 6ba9e36846ff158dd8680ababe644f73.png=a8c75315c6f595f8d1edb68e6d2f1c7a.png=06fd3218fdb55c0209cda7c084ae27e0.png ，eb69de4537bf711182fe429101048c7c.png)

72244fc81446dc4ad28082eecf1ffa71.png

44f1577f5a638c8df987627efa72df11.png

d373b36392ad6dff34afad7794994c24.png

8fdb9b849b883593ef13960d6fceac8a.png

1c02954a1f92a4fdb92ec0cc7c7f2725.png

2ab58dd70e3564111f1d451b931eb63c.png

6a039c5561617bf858447f28ae486aad.png

771d2cf539d88e5dbc4468aebfb8acef.png

其中9ba52365524736fa5e370523d63f084c.png满足条件，故最后一共有两个设计变量（92c46774ab2052d03d930b4715ecc2d4.png），7个约束条件。

4）优化计算

word/media/image1.gif此问题的图解见上图，有7个约束条件构成了改优化模型的可行域，而最优解在可行域内。

word/media/image2.gif优化方法选择：

该问题属于一般的约束非线性最优化类型，可以使用matlab优化工具箱里面的‘fmincon’函数进行求解。

word/media/image9.giffmincon里面算法的选择：fmincon里面一共提供了‘largescale’，‘ 'medium-scale’两种算法，由于此问题只有两个设计变量，维数较低，故采用‘medium-scale’算法。‘medium-scale’算法是采用SQP，算法中Hessian阵可以通过BFGS迭代，初始Hessian阵任给。注意BFGS公式中q项是需要计算目标函数梯度得到的。所以Hessian矩阵的近似计算是需要用到有限差分法。在采用‘'medium-scale’算法时，需提供其设计变量的初始点 b57899cdb003a4427061fb8108ba83e8.png的信息，而初始点的选择也将影响计算得收敛性和收敛速度，如果初始点选择得不恰当，可能最后函数不能收敛，得不到计算结果。

word/media/image10.gif精度的控制：为了得到更加精确地解，需要设置优化函数的控制精度，函数本身默认精度为1e-4，精度比较低，通过options的设置将精度提高到 1e-9,这样得到的结果更精确。

以上两点通过设置options参数即可：

options=optimset('largescale','off','display','off','Algorithm','active-set','TolFun',1e-9);

所有的程序编好以后，在命令窗口输入：youhua

得到的matlab的运行结果如下：

>> youhua

x =

4.1574

2.2909 %最优解

fval =

5.1899e-004 %目标函数最优点的值

exitflag =

5 %标志值，’5’表示重要方向导数小于规定的容许范围并且约束违背小于options.ToLCon

output =

iterations: 12 %迭代次数

funcCount: 40 %函数的评价次数

lssteplength: 1

stepsize: 7.6955e-005

algorithm: 'medium-scale: SQP, Quasi-Newton, line-search' %采用的中型算法

firstorderopt: 1.0832e-006 %一阶最优性条件

constrviolation: -1.0852e-006

message: [1x780 char] %跳出信息

lambda =

lower: [2x1 double]

upper: [2x1 double]

eqlin: [0x1 double]

eqnonlin: [0x1 double]

ineqlin: [5x1 double]

ineqnonlin: [2x1 double]

grad = %函数在最优点处梯度信息

1.0e-003 *

0.4888

0.4445

hessian = %函数在最优点处海塞矩阵

0.0016 0.0075

0.0075 0.0468

>>

5）结果分析

word/media/image1.gif采用fmincon求解的最优值：

ba3a4b7a23ddc0fd33af1465c33f9c2b.png=[4.1574; 2.2909];

f98b621effd5c32c8cebfbf808704a51.png；

采用算法：中型算法（mediun-scale）。

这与课本给出的最优解：ba3a4b7a23ddc0fd33af1465c33f9c2b.png=[4.1286;2.3325]，a48266edfc6acff7d1f0fda8564afba9.png=0.0156相比，计算精度更高，最优解的数值更精确，故计算准确度高。

word/media/image2.gif用matlab绘制输入——输出曲线关系图

上图中（单位为“度”）蓝色的线代表曲柄摇杆机构的实际输出角与输入角的关系，红色的线代表理想输出角与输入角的关系。

可以看出：

实际输出和理论输出曲线之间存在线性误差，其最大线性误差为d91e6fd56c47d780fc33f62a2b32b666.png，误差在允许的范围之内，故结果的可信度也较大，运用matlab优化工具箱计算所得结果正确。

小结

通过结合实际问题的分析，计算，求解，更加深入地了解和掌握机械优化设计的过程和步骤，比较重要的步骤是数学模型的建立，以及设计变量的选取，以及数学模型的尺度变换，根据机构实际工作需要，建立目标函数的约束条件等等，当数学模型建好以后，剩下的工作可以再matlab里面完成，而matlab里面的优化工具箱，给用户提供了多种优化函数，使用者只需要将数学模型按要求编写成子程序嵌入已有的优化程序即可。

在设计过程中也遇到一些困难，比如说在在用matlab计算时，计算机已知处于busy状态，得不到函数的最优解，最后反复的检查，终于找的了其原因，是由于初始点选择不恰当引起的，如果初始点选择得好，可以节省计算时间和计算空间，故初始点的选取比较重要。

附录

1. 编写目标函数M文件myfun.m：

function f=myfun(x)

f=0; %函数f赋初值

a0=acos(((1+x(1))^2-x(2)^2+25)/(10*(1+x(1)))); %初始计算点曲柄和摇杆的角度

b0=acos(((1+x(1))^2-x(2)^2-25)/(10*x(2)));

i=2;

while(i<=31) %设置迭代次数为30次

a(i)=a0+(pi/2)*(i/30); % 计算曲柄各分度的角度值

b(i)=b0+2*(a(i)-a0)^2/(3*pi); % 计算摇杆各分度的角度值

r=sqrt(26-10*cos(a(i)));

c(i)=acos((r^2+x(2)^2-x(1)^2)/(2*x(2)*r));

d(i)=acos((r^2+24)/(10*r));

if a(i)<=pi

e(i)=pi-c(i)-d(i); %计算摇杆输出的实际值

else if a(i)<=2*pi

e(i)=pi-c(i)+d(i);

end

end

a(1)=a0;

f=f+((b(i)-e(i))^2)*(a(i)-a(i-1)); %目标函数的计算

i=i+1;

end

2. 编写非线性不等式约束M文件constrain.m：

function[c ceq]=constrain(x)

c=[36-x(1)^2-x(2)^2-1.414*x(1)*x(2); x(1)^2+x(2)^2-1.414*x(1)*x(2)-16]; %非线性不等式约束

ceq=[];

3. 调用fmincon优化函数，建立youhua.m文件：

lb=[1;1]; %设计变量的下界

x0=[4;2]; %迭代初始点

A=[-1,0;0,-1;-1,-1;1,-1;-1,1]; %线性不等式约束

b=[-1;-1;-6;4;4];

options=optimset('largescale','off','display','off','Algohm','active-set','TolFun',1e-9); %采用中型算法，设计精度为1e-9

[x,fval,exitflag,output]=fmincon(@myfun,x0,A,b,[],[],lb,[],@constrain,options) %调用fmincon函数

本文来源：https://www.2haoxitong.net/k/doc/e3828390bf1e650e52ea551810a6f524cdbfcb00.html