2022年机械优化设计大作业2

高等流体力学班级：机设15学硕班学号：13姓名：张湘楠授课老师：毕新胜日期：2016年7月1日

一、研究报告内容：1、λ＝0.618的证明、一维搜索程序作业；2、单位矩阵程序作业；3、连杆机构问题＋自行选择小型机械设计问题或其他工程优化问题；（1）分析优化对象，根据设计问题的要求，选择设计变量，确立约束条件，建立目标函数，建立优化设计的数学模型并编制问题程序；（2）选择适当的优化方法，简述方法原理，进行优化计算；（3）进行结果分析，并加以说明。4、写出课程实践心得体会，附列程序文本。5、为响应学校2014年度教学工作会议的改革要求，探索新的课程考核评价方法，特探索性设立一开放式考核项目，占总成绩的5%。试用您自己认为合适的方式（书面）表达您在本门课程学习方面的努力、进步与收获。（考评将重点关注您的独创性、简洁性与可验证性）。二、研究报告要求1、报告命名规则：学号－姓名－《机械优化设计》课程实践报告.doc2、报告提交邮址：（收到回复，可视为提交成功）。追求：问题的工程性，格式的完美性，报告的完整性。不追求：问题的复杂性，方法的惟一性。评判准则：独一是好，先交为好；切勿拷贝。目录：λ＝0.618的证明、一维搜索程序作业关于的证明……………………4一维搜索的作业采用matlab进行编程……………5采用C语言进行编程………………7单位矩阵程序作业采用matlab的编程…………………9采用c语言进行编程…………………9机械优化工程实例连杆机构…………………11自选机构…………………16课程实践心得……………………20附列程序文本……………………21进步，努力，建议………………25一、λ＝0.618的证明、一维搜索程序作业①关于的证明黄金分割法要求插入点，的位置相对于区间两端具有对称性，即其中为待定常数。此外，黄金分割法还要求在保留下来的区间内再插入一点所形成的区间新三段，与原来的区间三段具有相同的比例分布。黄金分割法还要求在保留下来的区间内再插一点所形成的区间新三段，与原来的区间三段有相同的比例分布。bα1α211-λλ图一设原区间的长度为1，如图一所示，保留下来的区间长度为，区间缩短率为。为了保持相同的分别比例。插入新点应在位置上，在原区间的1位置应相当于在保留区间的位置。bα1α211-λλ图一故有：aα2aα22α3α1λλ(1-λ)λ2图二

编写0.618的程序，并计算下列问题程序框图程序框图（1）采用MATLAB进行编程%%%fun.m:黄金分割法求极值点%输入数据%a–搜索区间下限%b–搜索区间上限%e–精度%输出数据：%x–极小值点%其他常量：%c1,c2,c3–区间%r–黄金分割比例0.618%%%定义函数functionx=fun(a,b,e)r=0.618;c1=b-r*(b-a);c2=a+r*(b-a);y1=f(c1);y2=f(c2);while(abs((b-a)/b)>e)&&(abs((y2-y1)/y2)>e)ify1>=y2a=c1;c1=c2;y1=y2;c2=a+r*(b-a);y2=f(c2);elseb=c2;c2=c1;y2=y1;c1=b-r*(b-a);y1=f(c1);endendx=0.5*(a+b);end对f函数的确立functiony=f(x)y=(x-2)^2+3;endfunctiony=f(x)y=cosxend如果要计算y=(x-2)^2+3;的黄金分割法，则需要将图所示的f脚本中的函数写成如①所示如果要计算y=cosx；黄金分割法，则需要将图所示的f脚本中的函数写成如②所示fun函数表示对matlab的主程序语言。函数NO.2运算结果：函数NO.1运算结果：（2）采用C语言进行编程#include<stdio.h>#defineK0.618doublef(double);/*****函数值计算函数声明*****/voidmain(void){ doublea,b,size; doublea1,a2; intI; printf(“请输入区间两端点（端点值应大于0）:”); scanf(“%lf,%lf”,&a,&b);/*****输入端点值*****/ printf(“请输入精度:”); scanf(“%lf”,&size);/*****输入精度*****/ printf(“区间为(%lf,%lf),精度为%lf\n”,a,b,size); printf(“序号\ta1\t\ta2\t\tf(a1)\t\tf(a2)\n”); for(i=0;i<64;i++) printf(“-“); printf(“\n”); i=0;while((b-a)>size)/*****用精度控制循环次数*****/ { i++; a1=b-K*(b-a);/*****按0.618法插入两点*****/ a2=a+K*(b-a);printf(“%2d:\t%f\ta%f\t%f\t%f\n”,I,a1,a2,f(a1),f(a2)); /*****输出每次计算后a1,a2,f(a1),f(a2)的值*****/ if(f(a1)>=f(a2)) a=a1; else b=a2; } printf(“所求极小值点为:x=%lf\t极小值f(x)=%f\n”,a,f(b));}doublef(doublex)/*****函数值计算函数*****/{ doublef; f=(x-2)*(x-2)+3; return(f);}对于y=cosx，须在程序中加一个#include（math.h）头程序，以示我要调用函数语句。再将倒数第三行函数改为y=cos（x），再次编译运算即可。运行如图所示我们可以看出，大约在x=2处取到极值。与理论相符。运行如图所示我们可以看出，大约在π处取到极值。与理论相符。二、用简单的语句写一个单位矩阵①采用matlab构造：Matlab构造的矩阵最是简便，因为本身就有一个单位矩阵的函数。由于matlab是基于C语言而设立的一个数学运用软件，所以他的集成度非常高，而且具有很好的开放性，于是我们得知eye函数，将其构造。Eg：若要构造一个3维的单位矩阵，则输入eye（3）若要构造一个n为的单位矩阵，这输入eye（n）②用C语言构造程序如图所示：#include<stdio.h>intmain(){intI,j,n;while(scanf(“%d”,&n)==1){for(I=0;I<n;++i){for(j=0;j<n;++j)printf(“%d“,I==j);printf(“\n”);}}return0;}三、机械优化设计工程实例①连杆机构问题（1）连杆机构问题描述图1机构简图设计一曲柄连杆摇杆机构，要求曲柄从时，摇杆的转角最佳再现已知的运动规律：且=1，=5，为极位角，其传动角允许在范围内变化。（2）数学模型的建立设计变量：这里有两个独立参数和。因此设计变量为目标函数：将输入角分成30等分，并用近似公式计算，可得目标函数的表达式约束条件：GX(1)=-X(1)0 GX(2)=-X(2)0 GX(3)=-(X(1)+X(2))+6.00 GX(4)=-(X(2)+4.0)+X(1)0 GX(5)=-(4.0+X(1))+X(2)0 GX(6)=-(1.4142*X(1)*X(2)-X(1)**2-X(2)**2)-16.00 GX(7)=-(X(1)**2+X(2)**2+1.4142*X(1)*X(2))+36.00（3）程序编制C ====================== SUBROUTINEFFX(N,X,FX)C ====================== DIMENSIONX(N) COMMON/ONE/I1,I2,I3,I4,NFX,I6 NFX=NFX+1 P0=ACOS(((1.0+X(1))**2-X(2)**2+25.0)/(10.0*(1.0+X(1)))) Q0=ACOS(((1.0+X(1))**2-X(2)**2-25.0)/(10.0*X(2))) T=90.0*3.1415926/(180.0*30.0) FX=0.0 DO10K=0,30 PI=P0+K*T QE=Q0+2.0*(PI-P0)**2/(3.0*3.1415926) D=SQRT(26.0-10.0*COS(PI)) AL=ACOS((D*D+X(2)*X(2)-X(1)*X(1))/(2.0*D*X(2))) BT=ACOS((D*D+24.0)/(10.0*D))IF(PI.GE.0.0.AND.PI.LT.3.1415926)THENQI=3.1415926-AL-BTELSEQI=3.1415926-AL+BTENDIFIF(K.NE.0.OR.K.NE.30)THENFX=FX+(QI-QE)**2*TELSEFX=FX+(QI-QE)**2*T/2.0ENDIFCONTINUE RETURN ENDC ========================= SUBROUTINEGGX(N,KG,X,GX)C ========================= DIMENSIONX(N),GX(KG) GX(1)=-X(1) GX(2)=-X(2) GX(3)=-(X(1)+X(2))+6.0 GX(4)=-(X(2)+4.0)+X(1) GX(5)=-(4.0+X(1))+X(2) GX(6)=-(1.4142*X(1)*X(2)-X(1)**2-X(2)**2)-16.0 GX(7)=-(X(1)**2+X(2)**2+1.4142*X(1)*X(2))+36.0 RETURN ENDC ========================= SUBROUTINEHHX(N,KH,X,HX)C ========================= DIMENSIONX(N),HX(KH) X(1)=X(1) RETURN END输入数据如下：2,7,04.3,3.22.,0.2,0.01,0.00001,0.000010,1,03.5,2.8,5.5,4.5（4）可执行程序的生成将程序文本保存成文件t001.txt，存放在OPT_EXAM\FORTRAN\0-fortran目录下。使用DOS操作界面，OPT_EXAM\FORTRAN\0-fortran目录下，输入命令asumt.for+t001.txtt001.for回车。生成t001.for文件。for1t001;回车；④for2回车；⑤linkt001;回车，即生成可执行程序t001.exe。⑥从t001.dat文件中加载数据，运行t001.exe。t001<t001.dat>t001.rtf回车得到结果。操作过程如图：操作步骤（5）结果与分析结果显示如下：==============PRIMARYDATA==============N=2KG=7KH=0X:.4300000E+01.3200000E+01FX:.7847605E-02GX:-.4300000E+01-.3200000E+01-.1500000E+01-.2900000E+01-.5100000E+01-.6729390E+01-.1218939E+02X:.4300000E+01.3200000E+01FX:.7847605E-02GX:-.4300000E+01-.3200000E+01-.1500000E+01-.2900000E+01-.5100000E+01-.6729390E+01-.1218939E+02HX:PEN=.3974390E+01R=.2000000E+01C=.2000000E+00T0=.1000000E-01EPS1=.1000000E-04EPS2=.1000000E-04==============ITERATIONCOMPUTE==============*****IRC=1R=.2000000E+01PEN=.3974390E+01*****IRC=2R=.4000000E+00PEN=.3119361E+01*****IRC=3R=.8000001E-01PEN=.6407908E+00*****IRC=4R=.1600000E-01PEN=.1439420E+00*****IRC=5R=.3200000E-02PEN=.3733451E-01*****IRC=6R=.6400000E-03PEN=.8986510E-02*****IRC=7R=.1280000E-03PEN=.2346609E-02*****IRC=8R=.2560000E-04PEN=.9042178E-03*****IRC=9R=.5120001E-05PEN=.5396702E-03*****IRC=10R=.1024000E-05PEN=.4410887E-03*****IRC=11R=.2048000E-06PEN=.4128061E-03*****IRC=12R=.4096001E-07PEN=.4071001E-03*****IRC=13R=.8192003E-08PEN=.4014154E-03*****IRC=14R=.1638401E-08PEN=.3992769E-03*****IRC=15R=.3276801E-09PEN=.3987240E-03*****IRC=16R=.6553602E-10PEN=.3984505E-03*****IRC=17R=.1310720E-10PEN=.3983411E-03*****IRC=18R=.2621441E-11PEN=.3982889E-03===============OPTIMUMSOLUTION==============IRC=18ITE=35ILI=94NPE=447NFX=482NGR=0R=.2621441E-11PEN=.3982696E-03X:.4135127E+01.2315381E+01FX:.3982623E-03GX:-.4135127E+01-.2315381E+01-.4505072E+00-.2180254E+01-.5819746E+01-.7079845E+01-.3633157E-03HX:Stop-Programterminated.从结果中得知当X1=4.135127，X2=2.315381时，目标函数取得最小值,摇杆的转角最佳再现已知的运动规律：。②自选结构模型设计某带式输送机减速器的高速级齿轮传动。已知高速级输入功率P1=10kW,小齿轮转速n1=960r/min,传动比i=3.2。齿轮材料和热处理:大齿轮45号钢(调质)硬度为217～255HBS,小齿轮40Cr(调质)硬度为241～286HBS,工作寿命15年,假设每年工作300天,两班制,带式输送机工作平稳,转向不变。常规设计方案采用直齿圆柱齿轮:m=2.5,z1=30,Φd=1。解：（1）设计变量，单级直齿圆柱齿轮传动的中心距:f齿宽：f将m,z1,X=m，2）目标函数根据多目标优化的线性加权法建立体积最小的目标函数:f(x)=ω1·f1(x)+ω2·f2(x)=ω1·1+ⅈx1x其中:ω1,ω2是加权系数,且ω1+ω2=1,分别根据设计时径向和轴向安装位置的要求设定;取ω1=1表示要求中心距最小,取ω2=1则表示要求齿宽最小。3）约束条件模数的限制:对于传递动力的齿轮,通常要求模数不少于1.5-2,得约束条件:G1G（2）小齿轮齿数的限制:小齿轮齿数z1应不大于产生根切的最小齿数17,得约束条件G（3）齿宽系数的限制:由于ϕdmin≤ϕd≤ϕdG（4）齿面接触强度的限制，根据公式并查表得约束条件：GG（5）齿根弯曲强度的限制，根据公式查表得约束条件：GG4建立数学优化模型高速级齿轮传动多目标优化设计的数学模型为:（ω1取0.6，ω2取0.4）Fun(x)=min[ω1⋅1+3.2x1x22+X=m，GGGGGGGGG5编写程序并运行结果目标函数M文件：functionf=zhwm(x)f=0.6*2.1*x(1)*x(2)+0.4*x(1)*x(2)*x(3);约束函数M文件：function[cceq]=zhwy(x)c(1)=1.04*10^7-2.916*10^5*(x(1)*x(2))^3*x(3);c(2)=1.04*10^7-8.95*10^6*(x(1)*x(2))^3*x(3);c(3)=1.51*10^6-303.57*x(1)^3*x(2)^2*x(3);c(4)=1.42*10^6-2445.92*x(1)^3*x(2)^2*x(3);ceq=[];优化函数M文件：x0=[2321];lb=[1.5170.7];ub=[2inf1.15];u=[];[x,fval]=fmincon(@zhwm,x0,[],[],[],[],lb,ub,@zhwy)约束函数目标函数经过Matlab优化并圆整后的齿轮参数如下：设计参数模数m/mm齿数z齿宽系数ϕ中心距/mm普通设计2.5301157.5MATLAB优化1.791127.43771.149984.5254圆整后参数1.75301117.6经过计算，最小体积为87.15。四、课程学习心得从懒懒散散的寒假到匆匆忙忙的考试周，仿佛只是在弹指一挥间。然而我们就在这段时间内，又学了几门课程，而这些课程中，最有趣，最好玩的当然就是《机械优化设计了》《机械优化设计》是一门理论性非常强的一门课，刚开始的时候不得要领，但是在王卫荣老师的带领下，逐渐登堂入室，渐渐有所了解。首先，这门课所依托的最强武器便是数学。数学即是方法，思维，和逻辑。其次所依托的便是算法程序语言。归根结底的还是数学，因为数学赋予了算法以逻辑的力量，分析的过程，以及所能解决的一切问题。还好我的数学基础还算可以，逻辑思维还算通顺，再加上老师所给予的方法，学习这门课程，一路上虽跌跌撞撞，但也有所收获。所谓收获，从方法论的角度上讲，即时通过一个通用模型，解决一系列的问题，这是数学建模给我们的要求吗，同时也是《优化》这门课所赋予的核心思想。大二的时候，特别想搞建模，可是由于种种原因，其中最重要的原因就是懒惰，与之失之交臂。《优化》这门课则弥补了我的部分遗憾。因为没有建模的要求高，于是上手比较容易，同时又与我们自己的机械方面的课程紧密地结合在一起，过程轻车熟路。于是《优化》真是一门寓学于乐的科目，只恨课程太多时间太少，不然真要好好地研究一番。说道大作业，那可是真的难，第一问和第二问还好一点，第三问就不行了。说来惭愧，连杆机构问题是借鉴了上届学长的做法才学会的。因为太难了，所以就没有用fortran语言做这些问题，全部换成了matlab。真要感谢《优化设计》这门课程，如果不是这门课程的话，即使MATLAB已经学过了，我估计也不会回过头来复习，是《优化设计》让我重新拾起matlab，并且我可以说，凡是在优化大作业中使用过的算法，程序，我已经忘不掉了。因为fortran语言是一种相对陌生的语言，所以使用fortran变成的话有些力不从心，而我又是一个刨根究底的人，对于0.618的算法和单位矩阵的生成，我都是用的两种算法，这样保证思维的发散性。Matlab真的是一款很强大的软体，他讲机会用到的所有常见数学模型全部都涵括了，所以，我们，亲切地称它为，万能数学工具。当学习变成一种乐趣的时候，每当我们想做这种游戏的时候，时间，精力，甚至思路，仿佛全部有了，这就是积极所带来的意义。学习《优化设计》，传达出这样的信号：重要的不是考试的结果，而是解决问题的过程。所以，我觉得《优化设计》这门课的评价方式，在我进入大学以来，是第一次遇到，大作业，小作业的比值，竟然超过了考试成绩的比值。换句话说，老师认为，即使考试能考的很好，但是大作业，小作业不行，说明这只是一个会考试的工具，对于解决实际问题的能力还是有所欠缺。对！这样的评价很接地气。《优化设计》这门课已经结课了，考试已经结束了。把大作业交上去，仿佛这门课与我们就没有什么瓜葛了一样。不，绝对不是。解决问题的能力才刚刚开始。五、附列程序文本黄金分割法求极值点的matlab程序%%%fun.m:黄金分割法求极值点%输入数据%a–搜索区间下限%b–搜索区间上限%e–精度%输出数据：%x–极小值点%其他常量：%c1,c2,c3–区间%r–黄金分割比例0.618%%%定义函数functionx=fun(a,b,e)r=0.618;c1=b-r*(b-a);c2=a+r*(b-a);y1=f(c1);y2=f(c2);while(abs((b-a)/b)>e)&&(abs((y2-y1)/y2)>e)ify1>=y2a=c1;c1=c2;y1=y2;c2=a+r*(b-a);y2=f(c2);elseb=c2;c2=c1;y2=y1;c1=b-r*(b-a);y1=f(c1);endendx=0.5*(a+b);end对f函数的确立①functiony=f(x)y=(x-2)^2+3;end②functiony=f(x)y=cosxend黄金分割法求极值点的C语言程序#include<stdio.h>#defineK0.618doublef(double);/*****函数值计算函数声明*****/voidmain(void){ doublea,b,size; doublea1,a2; inti; printf("请输入区间两端点（端点值应大于0）:"); scanf("%lf,%lf",&a,&b);/*****输入端点值*****/ printf("请输入精度:"); scanf("%lf",&size);/*****输入精度*****/ printf("区间为(%lf,%lf),精度为%lf\n",a,b,size); printf("序号\ta1\t\ta2\t\tf(a1)\t\tf(a2)\n"); for(i=0;i<64;i++) printf("-"); printf("\n"); i=0;while((b-a)>size)/*****用精度控制循环次数*****/ { i++; a1=b-K*(b-a);/*****按0.618法插入两点*****/ a2=a+K*(b-a);printf("%2d:\t%f\ta%f\t%f\t%f\n",i,a1,a2,f(a1),f(a2)); /*****输出每次计算后a1,a2,f(a1),f(a2)的值*****/ if(f(a1)>=f(a2)) a=a1; else b=a2; } printf("所求极小值点为:x=%lf\t极小值f(x)=%f\n",a,f(b));}doublef(doublex)/*****函数值计算函数*****/{ doublef; f=(x-2)*(x-2)+3; return(f);}对于y=cosx，须在程序中加一个#include（math.h）头程序，以示我要调用函数语句。再将倒数第三行函数改为y=cos（x），再次编译运算即可。单位矩阵的matlab程序eye（x）单位矩阵的C语言程序#include<stdio.h>intmain(){intI,j,n;while(scanf(“%d”,&n)==1){for(I=0;I<n;++i){for(j=0;j<n;++j)printf(“%d“,I==j);printf(“\n”);}}return0;连杆机构的fortran语言C ====================== SUBROUTINEFFX(N,X,FX)C ====================== DIMENSIONX(N) COMMON/ONE/I1,I2,I3,I4,NFX,I6 NFX=NFX+1 P0=ACOS(((1.0+X(1))**2-X(2)**2+25.0)/(10.0*(1.0+X(1)))) Q0=ACOS(((1.0+X(1))**2-X(2)**2-25.0)/(10.0*X(2))) T=90.0*3.1415926/(180.0*30.0) FX=0.0 DO10K=0,30 PI=P0+K*T QE=Q0+2.0*(PI-P0)**2/(3.0*3.1415926) D=SQRT(26.0-10.0*COS(PI)) AL=ACOS((D*D+X(2)*X(2)-X(1)*X(1))/(2.0*D*X(2))) BT=ACOS((D*D+24.0)/(10.0*D))IF(PI.GE.0.0.AND.PI.LT.3.1415926)THENQI=3.1415926-AL-BTELSEQI=3.1415926-AL+BTENDIFIF(K.NE.0.OR.K.NE.30)THENFX=FX+(QI-QE)**2*TELSEFX=FX+(QI-QE)**2*T/2.0ENDIF10CONTINUE RETURN ENDC ========================= SUBROUTINEGGX(N,KG,X,GX)C ========================= DIMENSIONX(N),GX(KG) GX(1)=-X(1) GX(2)=-X(2) GX(3)=-(X(1)+X(2))+6.0 GX(4)=-(X(2)+4.0)+X(1) GX(5)=-(4.0+X(1))+X(2) GX(6)=-(1.4142*X(1)*X(2)-X(1