合肥工业大学《机械优化设计》课程实践研究报告班级：学号：姓名：授课教师：日期：2016年 11 月 7 日目录作业要求 (2)一、λ＝0.618的证明、一维搜索程序作业 (3)1、0.618法的基本思想 (3)2、关于0.618法中参数λ＝0.618的证明 (4)3、一维搜索程序作业 (5)二、单位矩阵程序作业 (8)三、注释最佳再现给定运动规律连杆机构优化设计 (10)问题模型子程序 (10)四、连杆机构问题＋其他工程优化问题 (12)1、连杆机构问题 (12)2、其他工程问题： (15)五、课程实践心得体会 (18)作业要求1、λ＝0.618的证明、一维搜索程序作业；2、单位矩阵程序作业；3、注释最佳再现给定运动规律连杆机构优化设计问题模型子程序；4、连杆机构问题 ＋ 自行选择小型机械设计问题或其他工程优化问题；（1）分析优化对象，根据设计问题的要求，选择设计变量，确立约束条件，建立目标函数，建立优化设计的数学模型并编制问题程序；（2）选择适当的优化方法，简述方法原理，进行优化计算； （3）进行结果分析，并加以说明。

5、写出课程实践心得体会，附列程序文本。

一、λ＝0.618的证明、一维搜索程序作业1、0.618法的基本思想“0.618法”，又称为黄金分割法，是常用的一种一维搜索试探方法，适用于[,]a b 区间上的任何单调函数求极小值问题。

0.618法是建立在区间消去法原理基础上的试探方法，即在搜索区间[,]a b 内适当插入两点1a 、1b ，且11a b ，如下图所示。

通过比较函数值1()f a 与1()f b 的大小，应用函数的单调性，可得出以下两种情况：1） 若11()()f a f b <，则取1[,]a b 为缩短后的区间。

2） 若11()()f a f b >，则取1[,]a b 为缩短后的区间。

然后在保留下来的区间上进行同样的处置，如此迭代下去，使搜索区间无限缩小，从而得到极小点的数值近似解。

2、关于0.618法中参数λ＝0.618的证明0.618法要求插入点1α，2α的位置相对于区间[,]a b 两端点具有对称性，即12()()b b a a b a αλαλ=--=+-假设[,][0,1]a b =，根据以上公式，得出分割后的区间如下图所示：进行再次分割时，0.618法要求在保留下来的区间内再插入一点，所形成的区间新三段与原来区间的三段具有相同的比例分布。

假设保留下来的区间为2[,]a α，区间长度为λ。

为了保持相同的比例分布，根据以上公式计算，新插入点3α应在(1)λλ-位置上，1α在原区间的1λ-位置相当于在保留区间的2λ位置。

所谓0.618法，就是使整段长与较长段的长度比值等于较长段与较短段长度的比值，即：11λλλ=- 通过计算解得0.618λ≈。

若保留下来的区间为1[,]b α，根据插入点的对称性，也能推得同样的λ值。

3、一维搜索程序作业0.618法的搜索过程如下：1) 给出初始搜索区间及收敛精度，将λ代入0.618。

2) 按坐标点计算公式计算1α和2α，并计算其对应的函数值。

3) 根据区间消去法原理缩短搜索区间。

4) 检查区间是否缩短到足够小和函数值是否收敛到足够近，如果条件不满足，则返回第二步。

5) 如果条件满足，取最后两实验点的平均值作为极小点的数值近似解。

程序框图如下：根据以上思路，下面借助C++,运用0.618法求解正弦函数的极小值。

初始区间为程序代码如下：// 0.618.cpp : 定义控制台应用程序的入口点。

#include "stdafx.h"#include<iostream>#include<iomanip>#include <cmath>#define PI 3.1415926using namespace std;float main(){float a=0,b=2*PI,t; //t为计算精度float a1,a2,y1,y2,A,min; //A为极小点，min为所对应的极小值float r=0.618;int p=0; //p为迭代次数cout<<"请输入计算精度，如0.1：";cin>>t;a1=b-r*(b-a);a2=a+r*(b-a);y1=sin(a1);y2=sin(a2);while(abs((b-a)/b)>=t || abs((y2-y1)/y2)>=t){if(y1>=y2){a=a1;a1=a2;y1=y2;a2=a+r*(b-a);y2=sin(a2);}else{b=a2;a2=a1;y2=y1;a1=b-r*(b-a);y1=sin(a1);}p++;}A=(a+b)/2;min=sin(A);cout<<"\n迭代次数为："<< p<<endl;cout<<"在【0，2π】区间内正弦函数的最小值为："<<setprecision(9) <<min<<endl;system("pause");return 0;}程序运行结果如下二、单位矩阵程序作业➢作业要求：编写生成n n 阶单位矩阵的程序，要求维数n在程序运行时从键盘输入。

➢C++程序代码如下：#include "stdafx.h"#include<iostream>using namespace std;int main(){int i,j,n;int a[200][200];cout<<"请输入n*n单位矩阵的阶数(n小于200)：";cin>>n;for(i=0;i<n;i++){for(j=0;j<n;j++)if(i==j)a[i][j]=1;elsea[i][j]=0;}cout<<" N ";for(i=0;i<n;i++){cout.width(3);cout<<i+1;}cout<<"\n";for(i=0;i<n+1;i++)cout<<"---";cout<<"\n";for(i=0;i<n;i++){cout.width(2);cout<<i+1<<":";for(j=0;j<n;j++){cout.width(3);cout<<a[i][j];}cout<<"\n";}system("pause");return 0;}➢程序运行结果如下：三、注释最佳再现给定运动规律连杆机构优化设计问题模型子程序➢作业要求对给定的最佳再现给定运动规律连杆机构优化设计问题模型子程序（FORTRAN语言）进行必要的注释。

➢注释结果连杆机构问题函数子程序目标函数==============SUBROUTINE FFX(N,X,FX) /定义目标函数（子函数）FFX/======================DIMENSION X(N) /定义一维数组X/COMMON /ONE/ I1,I2,I3,I4,NFX,I6 /开辟公共区域ONE，其中有整型变量I1、I2、I3、I4、I6，实型变量NFX/NFX=NFX+1P0=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) 将输入角30等份FX=0.0DO 10 K=0,30 /DO语句实现循环功能，循环初值为0，终值为30/PI=P0+K*TQE=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) THEN /IF条件语句：若 PI>0 ‘且’ PI<π，则执行件/QI=3.1415926-AL-BT /THEN块，条件满足选择此句/ELSEQI=3.1415926-AL+BT /ELSE块，条件不满足选择此句/ENDIF /结束IF选择结构/IF(K.NE.0 .OR. k.NE.30) THENFX=FX+(QI-QE)**2*TELSEFX=FX+(QI-QE)**2*T/2.0ENDIF10 CONTINUEEND /程序停止运行/不等约束=================SUBROUTINE GGX(N,KG,X,GX) /计算不等式约束函数子程序GGX/=========================DIMENSION X(N),GX(KG) /一维数组X、GX/GX(1)=1.0-X(1)GX(2)=1.0-X(2)GX(3)=1.0-5.0GX(4)=(1.0+X(1))-(X(2)+5.0)GX(5)=(1.0+X(2))-(X(1)+5.0)GX(6)=(1.0+5.0)-(X(1)+X(2))GX(7)=-(1.4142*X(1)*X(2)-X(1)**2-X(2)**2)-16.0GX(8)=-(X(1)**2+X(2)**2+1.4142*X(1)*X(2))+36.0 /不等式约束条件GX(1)、GX(2)、GX(3)、GX(4)、GX(5)、GX(6)、GX(7)、GX(8)/END /使程序停止运行/等式约束=================SUBROUTINE HHX(N,KH,X,HX) /计算等式约束函数子程序HHX/=========================DIMENSION X(N),HX(KH) /定义一维数组X、HX/X(1)=X(1)END /使程序停止运行/四、连杆机构问题 ＋其他工程优化问题1、连杆机构问题设计一曲柄连杆摇杆机构，要求曲柄从转到时，摇杆的转角最佳再现已知的运动规律：且，，为极位角，其传动角允许在范围内变化。

解：1、 设计变量的确定已知141,5l l ==，根据余弦定理有：2222212432301242()(1)25arccos arccos 210l l l l l l l l l l ϕ++-++-==⨯+⨯+极位角（）（1）222221243230343()(1)25arccos arccos 210l l l l l l l l l ψ+--+--==位置角因此，只有23,l l 为独立变量，则设计变量为1223()()T T x x x l l ==2、目标函数的建立目标函数可根据已知的运动规律与机构实际运动规律之间的偏差最小为指标来建立，先假设在和间进行30等分，则可得目标函数的表达式：3021()()min i Ei i f x ψψ==-→∑式中：02i i i i i i i παβϕπψπαβπϕπ--≤≤⎧=⎨-+≤≤⎩（）（）2223232222414arccos224arccos arccos210i i i i i i i i i r r l l rl r l l r rl r αβ==+-=+-+==2002()3Ei i ψψϕϕπ=+- 3、约束条件的确定曲柄存在的条件：1122131214()0()0()0g x l l g x l l l g x l l =-≤⎫⎪=-≤⎬⎪=-≤⎭为最短杆 414235123461324()0()0()0g x l l l l g x l l l l g x l l l l =+--≤⎫⎪=+--≤⎬⎪=+--≤⎭最短与最长杆之和小于或等于另两杆之和传动角要求：22223147232222341823()()arccos 13502()()45arccos 02o ol l l l g x l l l l l l g x l l +-+=-≤+--=-≤4、优化计算结果 运用MATLAB 编程求解：定义目标函数M 文件function f=myobj(x) f=0;% f 为目标函数Ao=acos(((1+x(1))^2-x(2)^2+25)/(10*(1+x(1)))); %Ao--初始极位角，x （1）--L2，x （2）--L3Bo=acos(((1+x(1))^2-x(2)^2-25)/(10*x(2))); %Bo--初始位置角for Ai=Ao:(pi/2)/30:(Ao+pi/2)%极位角Ai从Ao开始，分100等分，一直转到Ao+90°的角度Bei=Bo+2*(Ai-Ao)^2/(3*pi);%Bei为理想情况下曲柄每转过一个角度时，摇杆转过的角度r=sqrt(26-10*cos(Ai));m=acos((r^2+x(2)^2-x(1)^2)/(2*r*x(2)));n=acos((r^2+24)/(10*r));if(Ai>=0 && Ai<=pi)Bi=pi-m-n;elseBi=pi-m+n;endf=f+(Bi-Bei)^2; %Bi为曲柄每转过一个角度时，摇杆实际转过的角度end●定义非线性约束条件M文件function[c,ceq]=mycon(x)c=[acos((x(1)^2+x(2)^2-36)/(2*x(1)*x(2)))-135*(pi/180);45*(pi/180)-acos((x(1)^2+x(2)^2-16)/(2*x(1)*x(2)))];%c为非线性不等式约束（传动角要求）ceq=[];%ceq为非线性等式约束（无）●调用主程序x0=[6,3];A=[-1 0;0 -1;-1 -1;1 -1;-1 1];b=[-1;-1;-6;4;4];[x,fval]=fmincon(@myobj,x0,A,b,[],[],[],[],@nonlcon)程序运行结果如下：通过MATLAB 编程求解得到该问题的最优解为*(4.1287 2.3225)0.0076T x f ==2、其他工程问题：运装计划：某航空公司的运输机分前、后舱两部分装运客货，前舱容积为3160m ，最大装载量为10t ，后舱容积为3320m ，最大装载量为15t 。