5式中， ○
ri2?l32?l22ri2?x22?x126 ?i?arccos()?arccos() ○
2rl2rix2i3ri2?l42?l12ri2?247 ?i?arccos()?arccos() ○
2rl10ri4iri?l12?l42?2l1l4cos?i?26?10cos?i ○8
上式中r?BD
（四）、约束条件
1、曲柄存在条件
l1?l2;l1?l3,l1?l4?l2?l3l2?(l4?l1)?l3,l3?(l4?l1)?l2
2、曲柄与机架共线位置时的传动角（连杆BC和摇杆CD之间的夹角）:
最小传动角rmin?min?BCD?45? 最大转动角rmax?max?BCD?135?
由上面的分析可以算出：
rmin
?l22?l32??l42?l12???x2?x2?16??2???1?arccos???452l2l32x1x2?????? ○9
rmax?l22?l32??l42?l12???x2?x2?36??2???1?arccos???1352l2l32x1x2?????? ○10
四、利用MATLAB进行优化方案的解决
通过上面分析，可以将输入角分成 30 等分（m=30），并且转化为标准形式，
得到曲柄摇杆机构优化设计的标准数学模型为:
11 f?x?????Ei??i??min ○
i?1m2x??l2l3???x1x2? ○12
g1?x??1?x1?0??g2?x??1?x2?0??g3?x??6?x1?x2?0? ○13 g4?x??x1?x2?4?0??g5?x??x2?x1?4?0??g6?x??x12?x22?1.414x1x2?36?0?gx?36?x2?x2?1.414xx?01212?7??机械优化设计中的问题，大多数属于约束优化问题，此为非线性约束优化
问题，所以运用 MATLAB 来处理有约束的非线性多元函数的优化问题。 1、目标函数程序的建立
%目标函数的建立
function f=constraint(x); s=30;qb=1;jj=5;fx=0;
fa0=acos(((qb+x(1))^2+jj^2-x(2)^2)/(2*(qb+x(1))*jj));%曲柄初始角 pu0=acos(((qb+x(1))^2-x(2)^2-jj^2)/(2*x(2)*jj));%摇杆初始角 for i=1:s
fai=fa0+0.5*pi*(i/s);
pui=pu0+((2*(fai-fa0)^2)/(3*pi)); ri=sqrt(qb^2+jj^2-2*qb*jj*cos(fai));
alfi=acos((ri^2+x(2)^2-x(1)^2)/(2*ri*x(2))); bati=acos((ri^2+jj^2-qb^2)/(2*ri*jj)); if fai>0 && fai<=pi psi=pi-alfi-bati; else
psi=pi-alfi+bati; end
fx=fx+(pui-psi)^2; end f=fx;
2、约束函数程序的建立
%约束函数的建立
function [c,ceq]=fun1(x); qb=1;jj=5;
c(1)=45-rad2deg(acos((x(1)^2+x(2)^2-(jj-qb)^2)/(2*x(1)*x(2)))); %最小传动角约束
c(2)=45-180+rad2deg(acos((x(1)^2)+x(2)^2-(jj+qb)^2)/(2*x(1)*x(2))); %最大传动角约束 ceq=[];
约束函数建立
function [c,ceq]=fun1(x); qb=1;jj=5;
c(1)=45-rad2deg(acos((x(1)^2+x(2)^2-(jj-qb)^2)/(2*x(1)*x(2))));
%最小传动角约束
c(2)=45-180+rad2deg(acos((x(1)^2)+x(2)^2-(jj+qb)^2)/(2*x(1)*x(2))); %最大传动角约束 ceq=[];
3、运行的主程序
%运行的主程序 clc;clear; x0=[5,5];
A=[-1,0;0,-1;1,-1;-1,1;-1,-1];b=[-1;-1;4;4;-6]; Aeq=[];beq=[];
vlb=[1,1];vub=[5,5];
options=optimset('LargeScale','off','TolFun',1e-12);
[x,fval]=fmincon(@constraint,x0,A,b,Aeq,beq,vlb,vub,@fun1,options) x,fval
4、运行结果
Local minimum possible. Constraints satisfied.
fmincon stopped because the size of the current search direction is less than twice the default value of the step size tolerance and constraints are satisfied to within the default value of the constraint tolerance.
Active inequalities (to within options.TolCon = 1e-06): lower upper ineqlin ineqnonlin 5 x =
3.4554 2.5446 fval = 0.0035 x =
3.4554 2.5446 fval = 0.0035
即l1?1,l2?3.4554,l3?2.5446,l4?5，偏差fval为0.0035
5、因为MATLAB没有内置的rad2deg m文件，所以在运行时出现了一连串的错误提示，上面的截屏是在添加rad2deg m文件后的运行结果
function jdeg=rad2deg(rad) rad=rad*180/pi; deg=fix(rad);
xsbf=(rad-deg)*60; degf=fix(xsbf);
degm=(xsbf-degf)*60;
jdeg=deg+degf/100.0+degm/10000.0;
五、对MATLAB软件的运行结果的分析
1、把MATLAB运行的实际各杆所组合成的曲柄摇杆机构的输出角与期望的输出角运用MATLAB进行对比，编程的曲柄摇杆机构的输出角与输入角的关系线图如图3
图3曲柄摇杆机构的输出角与输入角的关系线图
根据图3可以知道，两条曲线在两点处相交，两条线的距离也是很接近，这个就可以说明MATLAB产生的结果是正确的、可靠的，我的设计方法是正确的。 把输入角变为50等分后，编程后的运行结果为： x =
3.4495 2.5505 fval =
0.0055 即l1?1,l2?3.4495,l3?2.55.5,l4?5，偏差fval为0.0055
再次运用此结果进行实际的和期望的曲柄摇杆机构的输出角与期望的输出角线图的编程，出现的结果更图3相同，这就再次验证了设计方法和过程的正确性。
2、利用MATLAB对曲柄摇杆机构在给定区间内的传动角变化进行编程，运行的图线如图4所示
图4曲柄摇杆机构在给定区间内的传动角变化曲线图
六、小结
为期8周时间的工程优化与MATLAB实现结束了，这8周的时间对我来说是短暂的，但是蛮有收获的，由以前的只知道设计但是不知道优化开始了解优化设计的重要性了，由以前的只是知道MATLAB这个软件到我开始能独自编写一些数学模型的程序了，由以前遇到什么事都是得过且过变得要自己独立去思考一个问题了，由以前的很少看一些科技文献到现在可以坐下来看有两篇了。这段时间的成长是有的，但更多的是对讨论和探讨一些问题的缺憾，对于曲柄摇杆机构结合MATLAB进行优化的一些结果自己仍然感到不满意，主要体现在对结果的综合性分析，对软件运行的结果进行正确性和主导性的深思，比如：设计思路中对设计变量的个数上自己的思维还是有些被限制，通过其他同学的有关设计变量的问题了解，对于不同的设计方法，出现的结果也是不同的；还有就是约束函数的到底主要是那几个起主要约束还是搞的不明白，这就暴露出我在机械原理课程上的不足。
在曲柄摇杆机构的设计中暴露出的一些问题真好反应了我对机械原理以及机械设计的掌握不够透彻的缺点，这也让我知道了自己的不足。工程优化与MATLAB实现这门课让我更懂得去独立思考，去发现一些以前没有注意到的不足，鼓励去思考某一个问题，并对这一个问题有自己的独到的见解。
参考文献
1、基于 M A T L A B 的平面四杆机构运动学分析 雷经发 1， 扈静 2， 魏涛 2 [1.安徽建筑工业学院 机械与电气工程系，合肥 2 3 0 6 0 1； 2.合肥工业大学 机械与汽车工程学院，合肥 2 3 0 0 0 9]
2、基于MATLAB的四杆机构的优化设计及仿真 刘菊蓉，吴金强，王旭飞，穆塔里夫[新