机械优化设计上机实践报告

1、机械优化设计上机实践报告班 级 ： 机械（茅以升）101姓 名 : 学 号 : 1004010510 成 绩 : 指导教师 : 张 迎 辉 日 期 : 2013.11.20 1 一维搜索方法上机实践报告1、写出所选择的一维搜索算法的基本过程、原理（可附流程图说明）。 （一）进退法1. 算法原理进退法是用来确定搜索区间（包含极小值点的区间）的算法，其理论依据是：为单谷函数（只有一个极值点），且为其极小值点的一个搜索区间，对于任意，如果，则为极小值的搜索区间，如果，则为极小值的搜索区间。因此，在给定初始点，及初始搜索步长的情况下，首先以初始步长向前搜索一步，计算。（1） 如果则可知搜索区间为，其中

2、待求，为确定，后退一步计算，为缩小系数，且，直接找到合适的，使得，从而确定搜索区间。（2） 如果则可知搜索区间为，其中待求，为确定，前进一步计算，为放大系数，且，知道找到合适的，使得，从而确定搜索区间。2. 算法步骤用进退法求一维无约束问题的搜索区间（包含极小值点的区间）的基本算法步骤如下：（1） 给定初始点，初始步长，令，；（2） 令，置；（3） 若，则转步骤（4），否则转步骤（5）；（4） 令，令，转步骤（2）；（5） 若，则转步骤（6）否则转步骤（7）；（6） 令，转步骤（2）；（7） 令，停止计算，极小值点包含于区间 （二）黄金分割法1、黄金分割法基本思路：黄金分割法适用于a，b区间上

3、的任何单股函数求极小值问题，对函数除要求“单谷”外不做其他要求，甚至可以不连续。因此，这种方法的适应面非常广。黄金分割法也是建立在区间消去法原理基础上的试探方法，即在搜索区间a，b内适当插入两点a1，a2，并计算其函数值。a1，a2将区间分成三段，应用函数的单谷性质，通过函数值大小的比较，删去其中一段，是搜索区间得以缩小。然后再在保留下来的区间上作同样的处理，如此迭代下去，是搜索区间无限缩小，从而得到极小点的数值近似解。2 黄金分割法的基本原理一维搜索是解函数极小值的方法之一，其解法思想为沿某一已知方向求目标函数的极小值点。一维搜索的解法很多，这里主要采用黄金分割法（0.618法）。该方法用不

4、变的区间缩短率0.618代替斐波那契法每次不同的缩短率，从而可以看成是斐波那契法的近似，实现起来比较容易，也易于人们所接受。图1 黄金分割法是用于一元函数f(x)在给定初始区间a,b内搜索极小点*的一种方法。它是优化计算中的经典算法，以算法简单、收敛速度均匀、效果较好而著称，是许多优化算法的基础，但它只适用于一维区间上的凸函数6，即只在单峰区间内才能进行一维寻优，其收敛效率较低。其基本原理是：依照“去劣存优”原则、对称原则、以及等比收缩原则来逐步缩小搜索区间7。具体步骤是：在区间a,b内取点：a1 ，a2 把a,b分为三段。如果f(a1)>f(a2)，令a=a1,a1=a2,a2=a+r

5、*(b-a)；如果f(a1)<f(a2) ，令b=a2，a2=a1,a1=b-r*(b-a),如果(b-a)/b和(y1-y2)/y2都大于收敛精度重新开始。因为a,b为单峰区间，这样每次可将搜索区间缩小0.618倍或0.382倍，处理后的区间都将包含极小点的区间缩小，然后在保留下来的区间上作同样的处理，如此迭代下去，将使搜索区a,b逐步缩小，直到满足预先给定的精度时，即获得一维优化问题的近似最优解。黄金分割法原理如图所示， 3 程序流程如下：开始4 实验所编程序框图r=0.618给定a=-3,b=5,收敛精度=0.001结束a*=(a+b)/2a1=b-r*(b-a)y1=f(a1)b

6、=a2a2=a1 y2=y1a2=a+r*(b-a)y2=f(a2)a=a1a1=a2 y1=y2y1>=y2a1=b-r*(b-a) y1=f(a1)a2=a+r*(b-a) y2=f(a2)否是(b-a)/b<和 （y2-y1）/y2<?否是算例1：minf(x)= x*x+2*x(1)C+程序如下：#include <math.h>#include <stdio.h> #define f(x) x*x+2*xdouble calc(double *a,double *b,double e,int *n) double x1,x2,s; if(fa

7、bs(*b-*a)<=e) s=f(*b+*a)/2); else x1=*b-0.618*(*b-*a); x2=*a+0.618*(*b-*a); if(f(x1)>f(x2) *a=x1; else *b=x2; *n=*n+1; s=calc(a,b,e,n); return s; main() double s,a,b,e; int n=0; scanf("%lf %lf %lf",&a,&b,&e); s=calc(&a,&b,e,&n); printf("a=%lf,b=%lf,s=%lf,

8、n=%dn",a,b,s,n); 2、程序运行结果：算例2：minf=x2-10*x+36 理论最优解：x*=5.0，f（x*）=11.0（1）MATLAB程序清单：function f=myfun_yi(x)f=x2-10*x+36>> fminbnd(myfun_yi,1,12)（2）运行结果：>> fminbnd(myfun_yi,1,12)f = 11.0407f = 18.8309f = 12.9691f = 11f = 11.0000f = 11.0000ans = 5（3）结果分析：由迭代程序f=11.0，ans=5，与理论结果相等算例3：min

9、f=x4-5*x3+4*x2-6*x+60 理论最优解：x*=3.2796，f（x*）=22.6590（1）MATLAB程序清单：function f=myfun_yi(x)f=x4-5*x3+4*x2-6*x+60>> fminbnd(myfun_yi,1,12)（2）运行结果：>> fminbnd(myfun_yi,1,12)f = 165.3948f = 1.5836e+03f = 24.8730f = 35.9194f = 23.9089f = 22.7621f = 31.7507f = 22.6673f = 22.6594f = 22.6590f = 22.6

10、590f = 22.6590f = 22.6590ans = 3.2796（3）结果分析：由迭代程序得f =22.659，ans =3.2796，与理论最优解相等2 无约束优化搜索方法上机实践报告1、写出所选择的无约束优化搜索算法的基本过程、原理（可附流程图说明）。 鲍威尔改进方法鲍威尔（Powell）法是直接利用函数值来构造共轭方向的一种方法 在鲍威尔基本算法中，每一轮迭代都用连结始点和终点所产生出的搜索方向去替换原向量组中的第一个向量，而不管它的“好坏”，这是产生向量组线性相关的原因所在。 在改进的算法中首先判断原向量组是否需要替换。如果需要替换，还要进一步判断原向量组中哪个向量最坏，然后

11、再用新产生的向量替换这个最坏的向量，以保证逐次生成共轭方向。2、程序计算结果分析：中间各步骤的结果分析及与理论计算结果分析对比。算例1：minf=4*(x(1)-5)2+(x(2)-6)2 初始点：x0=8;9,f(x0)=45 最优解：x*=5;6,f(x*)=0(1)MATLAB程序清单：function f=myfun_wuyueshu(x)f=4*(x(1)-5)2+(x(2)-6)2>> x,fval=fminunc(myfun_wuyueshu,x0)(2)运行结果：f = 45Warning: Gradient must be provided for trust-r

12、egion algorithm; using line-search algorithm instead. > In fminunc at 367 f = 45.0000f = 45.0000f = 23.5625f = 23.5625f = 23.5625f = 2.6958f = 2.6958f = 2.6958f = 1.3788f = 1.3788f = 1.3788f = 0.0054f = 0.0054f = 0.0054f = 6.4975e-05f = 6.4973e-05f = 6.4975e-05f = 6.1579e-09f = 6.1522e-09f = 6.14

13、43e-09f = 1.7876e-12f = 1.8627e-12f = 1.5586e-12Local minimum found.Optimization completed because the size of the gradient is less thanthe default value of the function tolerance.<stopping criteria details>x = 5.0000 6.0000fval = 1.7876e-12(3)结果分析：由迭代程序得x = 5.0000； 6.0000,fval =1.7876e-12，与理论

14、最优解相等。算例2：minf=(x(1)2+x(2)-11)2+(x(1)+x(2)2-7)2 初始点：x0=1;1,f(x0)=106 最优解：x*=3;2,f(x*)=0(1)MATLAB程序清单：function f=myfun_wuyueshu(x)f=(x(1)2+x(2)-11)2+(x(1)+x(2)2-7)2>> x,fval=fminunc(myfun_wuyueshu,x0)(2)运行结果：>> x0=1;1x0 = 1 1>> x,fval=fminunc(myfun_wuyueshu,x0)f = 106Warning: Gradie

15、nt must be provided for trust-region algorithm; using line-search algorithm instead. > In fminunc at 367 f = 106.0000f = 106.0000f = 29.5430f = 29.5430f = 29.5430f = 1.7450e+04f = 1.7450e+04f = 1.7450e+04f = 90.3661f = 90.3661f = 90.3661f = 0.3575f = 0.3575f = 0.3575f = 0.0179f = 0.0179f = 0.0179

16、f = 0.0064f = 0.0064f = 0.0064f = 1.0048e-06f = 1.0044e-06f = 1.0049e-06f = 4.8639e-09f = 4.8567e-09f = 4.8781e-09f = 5.2125e-12f = 5.8703e-12f = 5.7870e-12Local minimum found.Optimization completed because the size of the gradient is less thanthe default value of the function tolerance.<stopping

17、 criteria details>x = 3.0000 2.0000fval = 5.2125e-12(3)结果分析：由迭代程序得x=3;2,fval = 5.2125e-12,与理论最优解相等算例3:ff=x0*x0+2*x1*x1-4*x0-2*x0*x1;（1）鲍威尔改进算法C+程序清单：#include "stdio.h"#include "stdlib.h"#include "math.h"double objf(double x)double ff;ff=x0*x0+2*x1*x1-4*x0-2*x0*x1;ret

18、urn(ff);void jtf(double x0 ,double h0,double s ,int n,double a ,double b )int i; double *x3,h,f1,f2,f3; for (i=0;i<3;i+) xi=(double *)malloc (n*sizeof(double); h=h0; for(i=0;i<n;i+)*(x0+i)=x0i; f1=objf(x0); for(i=0;i<n;i+)*(x1+i)=*(x0+i)+h*si; f2=objf(x1); if(f2>=f1) h= -h0; for (i=0;i&l

19、t;n;i+)*(x2+i)=*(x0+i); f3=f1; for(i=0;i<n;i+) *(x0+i)= *(x1+i);*(x1+i)= *(x2+i); f1=f2; f2=f3;for(;)h=2. *h;for(i=0;i<n;i+) *(x2+i)=* (x1+i) +h*si;f3= objf(x2);if(f2<f3) break;else for(i=0;i<n;i+) *(x0+i)= *(x1+i);*(x1+i)= *(x2+i); f1=f2; f2=f3; if(h<0. )for(i=0;i<n;i+) ai=*(x2+i)

20、; bi=*(x0+i); else for(i=0;i<n;i+) ai=*(x0+i); bi=*(x2+i); for(i=0;i<3;i+) free(xi);double gold(double a,double b,double eps,int n,double xx) int i; double f1,f2,*x2,ff,q,w; for(i=0;i<2;i+) xi=(double*)malloc (n*sizeof(double); for(i=0;i<n;i+) *(x0+i)=ai+0.618*(bi-ai); *(x1+i)=ai+0.382*(

21、bi-ai); f1=objf(x0); f2=objf(x1); do if(f1>f2) for(i=0;i<n;i+) bi=*(x0+i); *(x0+i)=*(x1+i); f1=f2; for(i=0;i<n;i+) *(x1+i)=ai+0.382*(bi-ai); f2=objf(x1); else for(i=0;i<n;i+) ai=*(x1+i); *(x1+i)=*(x0+i); f2=f1; for(i=0;i<n;i+) *(x0+i)=ai+0.618*(bi-ai); f1=objf(x0); q=0; for(i=0;i<n

22、;i+) q=q+(bi-ai)*(bi-ai); w=sqrt(q); while(w>eps); for(i=0;i<n;i+) xxi=0.5*(ai+bi); ff=objf(xx); for(i=0;i<2;i+) free(xi); return(ff);double oneoptim(double x0,double s,double h0,double epsg,int n,double x)double *a,*b,ff;a=(double *)malloc(n*sizeof(double);b=(double *)malloc(n*sizeof(doubl

23、e);jtf(x0,h0,s,n,a,b);ff=gold(a,b,epsg,n,x);free(a);free(b);return(ff);double powell(double p,double h0,double eps,double epsg,int n,double x)int i,j,m;double *xx4,*ss,*s;double f,f0,f1,f2,f3,fx,dlt,df,sdx,q,d;ss=(double *)malloc(n*(n+1)*sizeof(double);s=(double *)malloc(n*sizeof(double);for (i=0;i&

24、lt;n;i+)for (j=0;j<=n;j+)*(ss+i*(n+1)+j)=0;*(ss+i*(n+1)+i)=1;for (i=0;i<4;i+)xxi=(double *)malloc(n*sizeof(double);for (i=0;i<n;i+)*(xx0+i)=pi;for(;)for (i=0;i<n;i+)*(xx1+i)=*(xx0+i);xi=*(xx1+i);f0=f1=objf(x);dlt=-1;for (j=0;j<n;j+)for (i=0;i<n;i+)*(xx0+i)=xi;*(s+i)=*(ss+i*(n+1)+j)

25、;f=oneoptim(xx0,s,h0,epsg,n,x);df=f0-f;if(df>dlt)dlt=df;m=j;sdx=0.;for (i=0;i<n;i+)sdx=sdx+fabs(xi-(*(xx1+i);if(sdx<eps)free(ss);free(s);for (i=0;i<4;i+)free(xxi);return(f);for (i=0;i<n;i+)*(xx2+i)=xi;f2=f;for (i=0;i<n;i+)*(xx3+i)=2.*(*(xx2+i)-(*(xx1+i);xi=*(xx3+i);fx=objf(x);f3=fx

26、;q=(f1-2*f2+f3)*(f1-f2-dlt)*(f1-f2-dlt);d=0.5*dlt*(f1-f3)*(f1-f3);if(f3<f1)|(q<d)if(f2<=f3)for (i=0;i<n;i+)*(xx0+i)=*(xx2+i);elsefor (i=0;i<n;i+)*(xx0+i)=*(xx3+i);elsefor (i=0;i<n;i+)*(ss+(i+1)*(n+1)=xi-(*(xx1+i);*(s+i)=*(ss+(i+1)*(n+1);f=oneoptim(xx0,s,h0,epsg,n,x);for(i=0;i<n;

27、i+)*(xx0+i)=xi;for (j=m+1;j<=n;j+)for (i=0;i<n;i+)*(ss+i*(n+1)+j-1)=*(ss+i*(n+1)+j);void main()double p=1,1;double ff,x2,x1,x2,f;ff=powell(p,0.3,0.001,0.0001,2,x);printf("shuchuzuiyoujie：n");x1=x1;x2=x2;f=ff;printf("x1=%f,x2=%f,f=%fn",x1,x2,f);getchar();（2）运行结果为：3约束优化搜索方法上机

28、实践报告1、写出所选择的约束优化搜索算法的基本过程、原理（可附流程图说明）。2、程序计算结果分析：中间各步骤的结果分析及与理论计算结果分析对比。算例1:minf=（x(1)-2)2+(x(2)-1)2; s.t g1(x)=x(1)2-x(2)<=0 g2(x)=x(1)+x(2)-2 <=0 初始点：x0=3;3,f(x0)=5 最优解：x*=1;1 f(x*)=1 (1)MATLAB程序清单：function f=myfun_constrain(x)f=（x(1)-2)2+(x(2)-1)2;function c,ceq=mycon(x)c=x(1)2-x(2);x(1)+x(

29、2)-2 ceq=>> x,fval=fmincon(myfun_constrain,x0,A,b,mycon)（2） 运行结果：>> A=-1,0;0,-1b=0;0x0=3;3A = -1 0 0 -1b = 0 0x0 = 3 3>> x,fval=fmincon(myfun_constrain,x0,A,b,mycon)Warning: The default trust-region-reflective algorithm does not solve problems with the constraints you havespecified

30、. FMINCON will use the active-set algorithm instead. For information on applicable algorithms, seeChoosing the Algorithm in the documentation. > In fmincon at 486f = 5c = 6 4ceq = f = 5.0000c = 6.0000 4.0000ceq = f = 5.0000c = 6.0000 4.0000ceq = f = 2.0000c = 1.0000 -1.0000ceq = f = 2.0000c = 1.0

31、000 -1.0000ceq = f = 2.0000c = 1.0000 -1.0000ceq = f = 1.0000c = 1.0e-15 * 0.9992 0.4441ceq = f = 1.0000c = 1.0e-07 * 0.2980 0.1490ceq = f = 1.0000c = 1.0e-07 * -0.1490 0.1490ceq = Local minimum found that satisfies the constraints.Optimization completed because the objective function is non-decreas

32、ing in feasible directions, to within the default value of the function tolerance,and constraints are satisfied to within the default value of the constraint tolerance.<stopping criteria details>Active inequalities (to within options.TolCon = 1e-06): lower upper ineqlin ineqnonlin 1 2x = 1.000

33、0 1.0000fval = 1.0000 (3)结果分析：由迭代程序得x =1.0000; 1.0000,fval =1.0000,与理论最优解相等 算例2.minf=1000-x(1)2-2*x(2)2-x(3)2-x(1)*x(2)-x(1)*x(3); S.t g1(x)=x(1)2+x(2)2+x(3)2-25<=0 g2(x)8*x(1)+14*x(2)+7*x(3)-56<=0 g3(x)=-x(1)<=0 g4(x)=-x(2)<=0 g5(x)=-x(3)<=0 初始点：x0=2;2;2,f(x0)=976 最优解：x*=3.512;0.217;

34、3.552,f(x*)=961.715（1）MATLAB程序清单：function f=myfun_constrain(x)f=1000-x(1)2-2*x(2)2-x(3)2-x(1)*x(2)-x(1)*x(3);function c,ceq=mycon(x)c=x(1)2+x(2)2+x(3)2-25;8*x(1)+14*x(2)+7*x(3)-56ceq=>> x,fval=fmincon(myfun_constrain,x0,A,b,mycon)（2）运行结果>> A=-1,0,0;0,-1,0;0,0,-1b=0;0;0x0=2;2;2A = -1 0 0

35、0 -1 0 0 0 -1b = 0 0 0x0 = 2 2 2>> x,fval=fmincon(myfun_constrain,x0,A,b,mycon)Warning: The default trust-region-reflective algorithm does not solve problems with the constraints you havespecified. FMINCON will use the active-set algorithm instead. For information on applicable algorithms, see

36、Choosing the Algorithm in the documentation. > In fmincon at 486c = -13 2ceq = c = -13.0000 2.0000ceq = c = -13.0000 2.0000ceq = c = -13.0000 2.0000ceq = c = -5.9320 0ceq = c = -5.9320 0.0000ceq = c = -5.9320 0.0000ceq = c = -5.9320 0.0000ceq = c = 1.1562 0.0000ceq = c = 1.1562 0.0000ceq = c = 1.

37、1562 0.0000ceq = c = 1.1562 0.0000ceq = c = 31.5713 0.0000ceq = c = 8.1851 0.0000ceq = c = 8.1851 0.0000ceq = c = 8.1851 0.0000ceq = c = 8.1851 0.0000ceq = c = 3.2553 0.0000ceq = c = 3.2553 0.0000ceq = c = 3.2553 0.0000ceq = c = 3.2553 0.0000ceq = c = 1.0789 -0.0000ceq = c = 1.0789 0.0000ceq = c = 1

38、.0789 0.0000ceq = c = 1.0789 0.0000ceq = c = 0.4793 0.0000ceq = c = 0.4793 0.0000ceq = c = 0.4793 0.0000ceq = c = 0.4793 0.0000ceq = c = 0.0039 -0.0000ceq = c = 0.0039 0.0000ceq = c = 0.0039 0.0000ceq = c = 0.0039 0.0000ceq = c = 1.0e-06 * 0.4884 0.0000ceq = c = 1.0e-06 * 0.8559 0.4187ceq = c = 1.0e

39、-06 * 0.4948 0.2086ceq = c = 1.0e-06 * 0.8644 0.3705ceq = Local minimum possible. Constraints satisfied.fmincon stopped because the predicted change in the objective functionis less than the default value of the function tolerance and constraints are satisfied to within the default value of the cons

40、traint tolerance.<stopping criteria details>Active inequalities (to within options.TolCon = 1e-06): lower upper ineqlin ineqnonlin 1 2x = 3.5120 0.2170 3.5523fval = 961.7152>>（3）结果分析：由迭代程序得x= 3.5120; 0.2170; 3.5523,fval =961.7152,与理论结果相等。 4平面连杆机构中再现已知运动规律的优化设计上机实践报告1、写出所针对此实际的优化问题，所确定的优化

41、目标函数和约束条件。例8.5 设计一曲柄摇杆机构（如图9所示），要求曲柄从转到时，摇杆的转角最佳再现已知的运动规律： ， 且已知,为极位角，其传动角允许在范围内变化。2、MATLAB程序清单：function f=objfun(x)n=30;L1=1;L2=5;fx=0;fa0=acos(L1+x(1)2-x(2)2+L22)/(2*(L1+x(1)*L2); %Çú±ú³õʼ½Çpu0=acos(L1+x(1)2-x(2)2-L22)/(2*(x(2)*L2); %Ò¡¸Ë³õʼ½Çfor i=1:n fai=fa0+0.5*pi*i/n; pui=pu0+2*(fai-fa0)2/(3*pi); ri=sqrt(L12+L22-2*L1*L2*cos(fai); alfai=acos(ri2+x(2)2-x(1)2)/(2*ri*x(2); betai=acos(ri2+L22-L12)/(2*ri*L2); if fai>0&am