matlab(四连杆优化设计)

机械优化设计在matlab中的应用
东南大学机械工程学院**
一优化设计目的：
在生活和工作中，人们对于同一个问题往往会提出多个解决方案，并通过各方面的论证从中提取最佳方案。

最优化方法就是专门研究如何从多个方案中科学合理地提取出最佳方案的科学。

由于优化问题无所不在，目前最优化方法的应用和研究已经深入到了生产和科研的各个领域，如土木工程、机械工程、化学工程、运输调度、生产控制、经济规划、经济管理等，并取得了显著的经济效益和社会效益。

二优化设计步骤：
1.机械优化设计的全过程一般可以分为如下几个步骤：
1)建立优化设计的数学模型；
2)选择适当的优化方法；
3)编写计算机程序；
4)准备必要的初始数据并伤及计算；
5)对计算机求得的结果进行必要的分析。

其中建立优化设计数学模型是首要的和关键的一步，它是取得正确结果的前提。

优化方法的选取取决于数学模型的特点，例如优化问题规模的大小，目标函数和约束函数的性态以及计算精度等。

在比较各种可供选用的优化方法时，需要考虑的一个重要因素是计算机执行这些程序所花费的时间和费用，也即计算效率。

2.建立数学模型的基本原则与步骤
设计变量的确定；
设计变量是指在优化设计的过程中，不断进行修改，调整，一直处于变化的参数称为设计变量。

设计变量的全体实际上是一组变量，可用一个列向量表示：
x=。

目标函数的建立；
选择目标函数是整个优化设计过程中最重要的决策之一。

当对某以设计性能有特定的要求，而这个要求有很难满足时，则针对这一性能进行优化会得到满意的效果。

目标函数是设计变量的函数，是一项设计所追求的指标的数学反映，因此它能够用来评价设计的优劣。

目标函数的一般表达式为：
f(x)=
，要根据实际的设计要求来设计目标函数。

约束条件的确定。

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

由若干个约束条件构成目标函数的可行域，而可行域内的所有设计点都是满足设计要求的，一般情况下，其设计可行域可表示为
在可行域中，任意设计点满足全部约束条件，称为可行解，但不是最优解，而优化设计就是要求出目标函数在可行域的最优解。

三实例分析（机械优化设计P241页例8-5）
设计一曲柄摇杆机构如图，要求：
曲柄
从
且已知
范围内变化。

分析：
1）设计变量的确定
决定机构尺寸的各杆长度，以及当摇杆按已知运动规律开始运行时，曲柄所载的位置角
应列为设计变量，即：
X=
=
考虑到机构的杆长按比例变化时，不会改变其运动规律，因此在计算时常取
，而其他杆长则按比例取为
的倍数。

若取曲柄的初始位置角为极位角，则
及相应摇杆
位置角
均为杆长的函数，几何图形关系如右图,其关系式为：
=arcos[
] （1）
=arcos[
] （2）
将
的长度代入上式（1），（2）得到：
=arcos[
]
=arcos[
]
因此，只有
为独立变量，设计变量减少，故最后的设计变量为：
X=
=
2) 目标函数的建立
目标函数可根据已知的运动规律与机构实际运动规律之间的偏差最小为指标来建立，即
f(x)=
→ min
式中
—期望输出角，
；
m —输入角等分数；
—实际输出角，由下图得：
a) 0≤
<
b)
（3）
式中
将输入角分成30等分，并用近似公式计算，可得目标函数的表达式：f(x)=
由题意知，传动角的变化范围是
，则上式中变量的最后形式可以写成：
将
带入（4）（5）得：
为当
时的理想输出角，其值在题目中已经给出：
3）约束条件的确定
曲柄摇杆机构应满足曲柄存在条件，可得
曲柄摇杆机构的传动角
可得
把约束条件简化(
=5
=
=
，
)
其中
满足条件，故最后一共有两个设计变量（
），7个约束条件。

4）优化计算
此问题的图解见上图，有7个约束条件构成了改优化模型的可行域，而最优解在可行域内。

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

fmincon里面算法的选择：fmincon里面一共提供了‘largescale’，‘ 'medium-scale’两种算法，由于此问题只有两个设计变量，维数较低，故采用‘medium-scale’算法。

‘medium-scale’算法是采用SQP，算法中Hessian阵可以通过BFGS迭代，初始Hessian阵任给。

注意BFGS公式中q项是需要计算目标函数梯度得到的。

所以Hessian矩阵的近似计算是需要用到有限差分法。

在采用‘'medium-scale’算法时，需提供其设计变量的初始点
的信息，而初始点的选择也将影响计算得收敛性和收敛速度，如果初始点选择得不恰当，可能最后函数不能收敛，得不到计算结果。

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

以上两点通过设置options参数即可：
options=optimset('largescale','off','display','off','Algorithm','ac tive-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）结果分析
采用fmincon求解的最优值：
=[4.1574; 2.2909];
；
采用算法：中型算法（mediun-scale）。

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

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

可以看出：
实际输出和理论输出曲线之间存在线性误差，其最大线性误差为
，误差在允许的范围之内，故结果的可信度也较大，运用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','activ e-set','TolFun',1e-9); %采用中型算法，设计精度为1e-9
[x,fval,exitflag,output]=fmincon(@myfun,x0,A,b,[],[],lb,[],@constra in,options) %调用fmincon函
数。