合肥工业大学《机械优化设计》课程实践研究报告班级:11级机设学号:2011姓名:授课老师:王卫荣日期:2014年4月日一、研究报告内容:1、λ=0.618的证明、一维搜索程序作业;黄金分割法要求插入点α1、α2的位置相对于区间{a,b}两端点具有对称性,α1=b−λ(b−a)α2=a+λ(b−a)除对称要求外,黄金分割法还要求在保留下来的区间内再插入一点所形成的区间新三段,与原来区间的三段具有相同的比例分布。设原区间(a,b)长度为1,保留下来的区间(a,α2)长度为λ,区间缩短率为λ。为了保持相同的比例分布,新插入点α3应在λ(1-λ)位置上,α1在原区间的1-λ位置应相当于在保留区间的λ2位置。故有1-λ=λ2即λ2+λ-1=0取方程正数解,得λ=√5−12≈0.6180.618法C语言程序:#includestdio.h#includemath.hfloatm=0.618;floatfun(floatt){floaty;y=cos(t);returny;}main(){floata,b,eps;printf(\min=);scanf(%f,&a);%输入函数下限%printf(\max=);scanf(%f,&b);%输入函数上限%floatt1,t2,t,f1,f2,min;printf(eps=);scanf(%f,&eps);%输入精度%while((b-a)/b=eps){t1=a+(1-m)*(b-a);t2=a+m*(b-a);f1=fun(t1);f2=fun(t2);if(f1=f2){a=t1;t1=t2;f1=f2;t2=a+m*(b-a);f2=fun(t2);}else{b=t2;t2=t1;f2=f1;t1=a+(1-m)*(b-a);f1=fun(t1);}}t=(a+b)/2;min=fun(t);printf(最优点t=%f\n,t);%输出最优点t%printf(最优值f=%f\n,min);}%输出最优值f%1.Y=cos(t)2.y=(t-2)*(t-2)+32.单位矩阵程序作业程序如下所示:#includestdio.hvoidmain(void){inti,j;intn;inta=1;intb=0;printf(阶数n=);scanf(%d,&n);for(i=1;i=n;i++){for(j=1;j=n;j++)if(i==j)printf(%2d,a);elseprintf(%2d,b);printf(\n);}}3连杆机构问题+自行选择小型机械设计问题或其他工程优化问题;一、连杆机构问题:问题描述:图1现优化一曲柄连杆机构,如图1所示,已知曲柄长度L1为44mm,机架长度L4为220mm,,要求当曲柄的转角在[φ0,φ0+π/2]时,对应的摇杆的输出角为Ψi,且两者满足对应函数关系Ψi=Ψ0+(φ0-φi)2,φ0和Ψ0分别对应于四连杆在初始位置时曲柄和摇杆的位置角。要求机构传动角的范围是[π/4,3π/4],优化该问题使得从动件的一系列实际输出角与期望实现函数Ψ=f(φ)的值的平方偏差之和最小。模型建立1、设计变量曲柄摇杆机构按照原动件和从动件的对应关系可知其有5个独立参数,对于图1分别为曲柄长度L1,连杆长度L2,摇杆长度L3,机架长度L4,曲柄初始角φ0和摇杆的初始角Ψ0,由于L1和L4已知,且由图1的几何关系知:所以φ0和Ψ0已不再是独立参数,而是杆长的函数。经上分析独立变量只有L2和L3。因此,选择连杆长度L2和摇杆长度L3作为设计变量。即:X=[L2L3]T=[X1X2]T2、目标函数图2图3由上面图2和图3中机构的几何关系可得如下的运动规律:S为角度区间的分段数;Ψsi为机构的实际输出角,计算式为:Ψsi=根据图中的角度关系求得:所以根据本机构设计问题,以机构实际输出角Ψsi与理论输出角Ψi的平方偏差最小原则来建立目标函数。优化目标函数表达式:3、约束条件根据已知条件,该机构的约束条件有两方面:一是传动运动过程中的传动角应大于45度且小于135度;二是保证四杆机构满足曲柄存在条件。(1)保证传动角约束图4图5根据图4和图5中机构处于最大传动角和最小传动角时的连杆几何关系,由余弦定理知将L2=X1,L3=X2,L1=44,L4=220代入上面两式得(2)、曲柄存在条件由机械原理的知识可知,曲柄存在的条件为:将已知的杆长和设计变量代入上述条件得:经分析上述杆长条件不起约束作用,实际起作用的约束只g1(X)和g2(X),所以最终的数学模型为:优化设计1、程序运行结果CommandWindow:Workspace:根据程序运行结果可知,exitflag的值为1,说明目标函数收敛到局部最优解,优化效果较为理想,此时目标函数f的值为0.0121rad2,连杆的长度为181.5602mm,摇杆的长度为102.4337mm。2、结果分析当曲柄在[φ0,φ0+]范围内转动时,摇杆输出角与期望实现函数Ψ=f(φ)的平方偏差值之和最小为0.0121rad2,最优点位于约束条件g2(X)=0上。二、桁架优化实例1、问题如图所示,为某一特殊桁架的一部分,杆2距O点40cm处有一支点C。为了固定桁架,现欲在杆1和2上设置支点A和B,用来连接杆3(可拆卸)。已知当桁架固定时,杆1和2成60度角;而且,杆1右边有一段长为25cm的重要部位,不能设置支点。卸去杆3、收起桁架时,支点A的位置不能高于BC段中点D。求取支点A、B的位置,使得杆3的长度尽量小,以节省材料。图1桁架结构示意图2、数学模型设A、B两点距离O点的长度分别为1x和2x,而桁架固定时杆1和2成60度角。所以杆3的长度为2122213xxxxl。由图1可知,251x且)40(2121xx,即0251x且040221xx。设TxxX],[21,取222123)(xxlXf。因此,数学模型为:极小化目标函数:212221)(minxxxxXf约束条件:025)(11xXg0402)(212xxXg3、求解数学模型外点罚函数法求解构造外点法罚函数,如下:})]0,402[max()]0,25{[(max(),(22121)(212221)(xxxMxxxxMXkk程序流程图如图2所示:图2外点罚函数法程序流程图程序步骤:①选择适当的初始罚因子)0(M、初始点)0(X、收敛精度和罚因子系数c。在本程序中分别取1)0(M,]20,25[)0(X,610,c=8。令迭代步数k=0。②采用牛顿法求无约束问题),(min)(kMX的极值点)()(*kMX。③检验迭代终止准则,若满足||)()(||)1()(*kkMXMX及|)]([)]([)]([|)1(*)1(*)(*kkkMXfMXfMXfk=0i=0求)()(iX与Hessian矩阵)()(iXH||)(||)(iX)()(*)(ikXMX||)()(||)1()(*kkMXMX|)]([)]([)]([|)1(*)1(*)(*kkkMXfMXfMXf输出*X和)(*XfYN)()]([)(1)()()1(iiiiXXHXXi=i+1)()(*)0(kMXX)()1(kkcMMk=k+1YN结束牛顿法求),(min)(kMX的极值点)()(*kMX给定)0(X、)0(M、c、则停止迭代计算,输出最优点)()(**kMXX;否则,转入步骤④。④取)()1(kkcMM,)()(*)0(kMXX,k=k+1,转入步骤②继续迭代。具体程序请看附一。运行结果:X*=[25.0000,10.0000]f(X*)=475.0000因此,A、B两支点与O的距离分别为25cm、10cm,杆3的最小长度为475cm。目标函数曲线图与目标函数等值线图分别如图3和图4所示。优化路径如图4所示。图3目标函数曲线图图4目标函数等值线图4、结论分析罚因子系数c对外点罚函数法的影响本次程序中,c取值为8,运行步数k=11。若取c=4,则运行步数k=16;取c=16,则运行步数k=9;取c=32,则运行步数k=8;取c=64,则运行步数k=7。由此可知,罚因子系数c的大小会影响程序的迭代次数k。c的值取得越大,运行步数k越小,程序收敛速度越快,效率越高。但对于c的其他一些取值,如5、7、9等,会导致罚函数形态变坏,使迭代出现问题,导致程序运行失败。因此,需选取合适的罚因子系数c。附程序:closeallclearallclcsymsx1x2M;%M为罚因子。m(1)=1;c=8;%c为递增系数。赋初值。a(1)=25;b(1)=20;f=x1^2+x2^2-x1*x2+M*((25-x1)^2+(2*x1-x2-40)^2);%外点罚函数f0(1)=200;%求偏导、Hessian元素fx1=diff(f,'x1');fx2=diff(f,'x2');fx1x1=diff(fx1,'x1');fx1x2=diff(fx1,'x2');fx2x1=diff(fx2,'x1');fx2x2=diff(fx2,'x2');%外点法M迭代循环fork=1:100x1=a(k);x2=b(k);M=m(k);%牛顿法求最优值forn=1:100f1=subs(fx1);%求解梯度值和Hessian矩阵f2=subs(fx2);f11=subs(fx1x1);f12=subs(fx1x2);f21=subs(fx2x1);f22=subs(fx2x2);if(double(sqrt(f1^2+f2^2-f1*f2))=1e-6)%最优值收敛条件a(k+1)=double(x1);b(k+1)=double(x2);f0(k+1)=double(subs(f));break;elseX=[x1x2]'-inv([f11f12;f21f22])*[f1f2]';x1=X(1,1);x2=X(2,1);endendif(double(sqrt((a(k+1)-a(k))^2+(b(k+1)-b(k))^2-(a(k+1)-a(k))*(b(k+1)-b(k))))=1e-6)&&(double(abs((f0(k+1)-f0(k))/f0(k)))=1e-6)%罚因子迭代收敛条件%输出最优点坐标,罚因子迭代次数,最优值a(k+1)b(k+1)kf0(k+1)break;elsem(k+1)=c*m(k);endend%绘制目标函数曲线图xx1=0:0.5:50;xx2=0:0.5:50;fori=1:length(xx1)forj=1:length(xx2)if((2*xx1(i)-xx2(j)-40=0)&&(25-xx1(i)=0))Z(i,j)=xx1(i)^2+xx2(j)^2-xx1(i)*xx2(j);elseZ(i,j)=0;endendendfigure(1);surf(xx1,xx2,Z);axis([05005002500])title('目标函数曲线图');xlabel('x1');ylabel('x2');%绘制目标函数等值线图,并画出优化路径figure(2);x11=-5:0.5:35;x12=-5:0.5:35;[xx11,xx12]=meshgrid(x11,x12);F=xx11^2+xx12^2-xx11*xx12;axis([-550-550]);contour(xx11,xx12,F);title('目标函数等值线');xlabel('x1');ylabel('x2');holdonplot(a,b,'r+-');%绘制优化路径自行选择小型机械设计问题或其他工程优化问题;某航空公司每天有三个航班服务于A,B,C,H四个城市,其中城市H是可供转机使用的,三个航班的出发地-目的地分别为AH,HB,HC,可搭乘旅客的最大数量分别为120人,100人,110人,机票的价格分头等舱和经济舱两类.经过市场