機械優(yōu)化設計上機實踐報告

1、機械優(yōu)化設計上機實踐報告班 級 ： 機械（茅以升）101姓 名 : 學 號 : 1004010510 成 績 : 指導教師 : 張 迎 輝 日 期 : 2013.11.20 1 一維搜索方法上機實踐報告1、寫出所選擇的一維搜索算法的基本過程、原理（可附流程圖說明）。 （一）進退法1. 算法原理進退法是用來確定搜索區(qū)間（包含極小值點的區(qū)間）的算法，其理論依據是：為單谷函數（只有一個極值點），且為其極小值點的一個搜索區(qū)間，對于任意，如果，則為極小值的搜索區(qū)間，如果，則為極小值的搜索區(qū)間。因此，在給定初始點，及初始搜索步長的情況下，首先以初始步長向前搜索一步，計算。（1） 如果則可知搜索區(qū)間為，其中

2、待求，為確定，后退一步計算，為縮小系數，且，直接找到合適的，使得，從而確定搜索區(qū)間。（2） 如果則可知搜索區(qū)間為，其中待求，為確定，前進一步計算，為放大系數，且，知道找到合適的，使得，從而確定搜索區(qū)間。2. 算法步驟用進退法求一維無約束問題的搜索區(qū)間（包含極小值點的區(qū)間）的基本算法步驟如下：（1） 給定初始點，初始步長，令，；（2） 令，置；（3） 若，則轉步驟（4），否則轉步驟（5）；（4） 令，令，轉步驟（2）；（5） 若，則轉步驟（6）否則轉步驟（7）；（6） 令，轉步驟（2）；（7） 令，停止計算，極小值點包含于區(qū)間 （二）黃金分割法1、黃金分割法基本思路：黃金分割法適用于a，b區(qū)間上

3、的任何單股函數求極小值問題，對函數除要求“單谷”外不做其他要求，甚至可以不連續(xù)。因此，這種方法的適應面非常廣。黃金分割法也是建立在區(qū)間消去法原理基礎上的試探方法，即在搜索區(qū)間a，b內適當插入兩點a1，a2，并計算其函數值。a1，a2將區(qū)間分成三段，應用函數的單谷性質，通過函數值大小的比較，刪去其中一段，是搜索區(qū)間得以縮小。然后再在保留下來的區(qū)間上作同樣的處理，如此迭代下去，是搜索區(qū)間無限縮小，從而得到極小點的數值近似解。2 黃金分割法的基本原理一維搜索是解函數極小值的方法之一，其解法思想為沿某一已知方向求目標函數的極小值點。一維搜索的解法很多，這里主要采用黃金分割法（0.618法）。該方法用不

4、變的區(qū)間縮短率0.618代替斐波那契法每次不同的縮短率，從而可以看成是斐波那契法的近似，實現起來比較容易，也易于人們所接受。圖1 黃金分割法是用于一元函數f(x)在給定初始區(qū)間a,b內搜索極小點*的一種方法。它是優(yōu)化計算中的經典算法，以算法簡單、收斂速度均勻、效果較好而著稱，是許多優(yōu)化算法的基礎，但它只適用于一維區(qū)間上的凸函數6，即只在單峰區(qū)間內才能進行一維尋優(yōu)，其收斂效率較低。其基本原理是：依照“去劣存優(yōu)”原則、對稱原則、以及等比收縮原則來逐步縮小搜索區(qū)間7。具體步驟是：在區(qū)間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為單峰區(qū)間，這樣每次可將搜索區(qū)間縮小0.618倍或0.382倍，處理后的區(qū)間都將包含極小點的區(qū)間縮小，然后在保留下來的區(qū)間上作同樣的處理，如此迭代下去，將使搜索區(qū)a,b逐步縮小，直到滿足預先給定的精度時，即獲得一維優(yōu)化問題的近似最優(yōu)解。黃金分割法原理如圖所示， 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 理論最優(yōu)解：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 理論最優(yōu)解：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，與理論最優(yōu)解相等2 無約束優(yōu)化搜索方法上機實踐報告1、寫出所選擇的無約束優(yōu)化搜索算法的基本過程、原理（可附流程圖說明）。 鮑威爾改進方法鮑威爾（Powell）法是直接利用函數值來構造共軛方向的一種方法 在鮑威爾基本算法中，每一輪迭代都用連結始點和終點所產生出的搜索方向去替換原向量組中的第一個向量，而不管它的“好壞”，這是產生向量組線性相關的原因所在。 在改進的算法中首先判斷原向量組是否需要替換。如果需要替換，還要進一步判斷原向量組中哪個向量最壞，然后

11、再用新產生的向量替換這個最壞的向量，以保證逐次生成共軛方向。2、程序計算結果分析：中間各步驟的結果分析及與理論計算結果分析對比。算例1：minf=4*(x(1)-5)2+(x(2)-6)2 初始點：x0=8;9,f(x0)=45 最優(yōu)解：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、最優(yōu)解相等。算例2：minf=(x(1)2+x(2)-11)2+(x(1)+x(2)2-7)2 初始點：x0=1;1,f(x0)=106 最優(yōu)解：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,與理論最優(yōu)解相等算例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約束優(yōu)化搜索方法上機

28、實踐報告1、寫出所選擇的約束優(yōu)化搜索算法的基本過程、原理（可附流程圖說明）。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 最優(yōu)解：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,與理論最優(yōu)解相等 算例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 最優(yōu)解：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平面連桿機構中再現已知運動規(guī)律的優(yōu)化設計上機實踐報告1、寫出所針對此實際的優(yōu)化問題，所確定的優(yōu)化

41、目標函數和約束條件。例8.5 設計一曲柄搖桿機構（如圖9所示），要求曲柄從轉到時，搖桿的轉角最佳再現已知的運動規(guī)律： ， 且已知,為極位角，其傳動角允許在范圍內變化。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