机械优化设计实验报告(23页).doc

-机械优化设计实验报告-第 23 页机械优化设计实验报告目录1.进退法确定初始区间31.1 进退法基本思路31.2 进退法程序框图31.3 题目41.4 源程序代码及运行结果42.黄金分割法52.2黄金分割法流程图52.3 题目52.4 源程序代码及结果53.牛顿型法63.1牛顿型法基本思路63.2 阻尼牛顿法的流程图63.3 题目73.4 源程序代码及结果74.鲍威尔法84.1 鲍威尔法基本思路84.2 鲍威尔法流程图843 题目94.4 源程序代码及结果95. 复合形法165.1 复合行法基本思想165.3 源程序代码及结果166. 外点惩罚函数法246.1解题思路：246.2 流程框图246.3 题目246.4 源程序代码及结果247.机械设计实际问题分析307.2计算过程如下307.3 源程序编写328.报告总结331.进退法确定初始区间1.1 进退法基本思路：按照一定的规则试算若干个点，比较其函数值的大小，直至找到函数值按“高-低-高”变化的单峰区间。1.2 进退法程序框图1.3 题目：用进退法求解函数的搜索区间1.4 源程序代码及运行结果#include <stdio.h>#include <math.h>main()float h,h0,y1,y2,y3,a1=0,a2,a3,fa2,fa3;scanf("h0=%f,y1=%f",&h0,&y1);h=h0;a2=h;y2=a2*a2-7*a2+10;if (y2>y1)h=-h;a3=a1;y3=y1;loop:a1=a2;y1=y2;a2=a3;y2=y3;a3=a2+2*h;y3=a3*a3-7*a3+10;if (y3<y2)goto loop;elseprintf("a1=%f,a2=%f,a3=%f,y1=%f,y2=%f,y3=%fn",a1,a2,a3,y1,y2,y3);搜索区间为0 62.黄金分割法2.1黄金分割法基本思路：通过不断的缩短单峰区间的长度来搜索极小点的一种有效方法。按（） 缩小 比较大小 确定取舍区间。 2.2黄金分割法流程图2.3 题目：对函数，给定搜索区间时，试用黄金分割法求极小点2.4 源程序代码及结果：f=inline('x2-7*x+9')a=0;b=8;eps=0.001;a1=b-0.618*(b-a);y1=f(a1);a2=a+0.618*(b-a);y2=f(a2);while(abs(b-a)>eps) if(y1>=y2) a=a1; a1=a2; y1=y2; a2=a+0.618*(b-a); y2=f(a2); else b=a2;a2=a1;y2=y1; a1=b-0.618*(b-a); y1=f(a1); endendxxx=0.5*(a+b)f = Inline function: f(x) = x2-7*x+9xxx = 3.49973.牛顿型法 3.1牛顿型法基本思路：在邻域内用一个二次函数 来近似代替原目标函数，并将 的极小点作为对目标函数求优的下一个迭代点。经多次迭代，使之逼近目标函数的极小点。3.2 阻尼牛顿法的流程图：3.3 题目：用牛顿阻尼法求函数的极小点3.4 源程序代码及结果： k=0; ptol=1.0e-5;xk=input('input x0:') itcl=1;1; while norm(itcl)>=ptol f1=4*xk(1,1)3-24*xk(1,1)2+50*xk(1,1)-4*xk(2,1)-32;-4*xk(1,1)+8*xk(2,1); G=12*xk(1,1)2-48*xk(1,1)+50,-4;-4,8; dk=-inv(G)*f1; a=-(dk'*f1)/(dk'*G*dk); xk=xk+a*dk; itcl=a*dk; k=k+1; endf=(xk(1,1)-2)4+(xk(1,1)-2*xk(2,1)2; fprintf('n ÓÃ×èÄáÅ£¶Ù·¨µü´ú %d ´ÎºóµÃµ½ ¼«Ð¡µã x*¼°¼«Ð¡Öµ f Ϊ:n',k); disp(xk); disp(f); 结果显示：input x0:1;1 用阻尼牛顿法迭代 27 次后得到 极小点 x*及极小值 f 为: 2.0000 1.00001.3270e-0194.鲍威尔法 4.1 鲍威尔法基本思路：在不用导数的前提下，在迭代中逐次构造G的共轭方向。 4.2 鲍威尔法流程图：43 题目：求函数f（x）= x0*x0+x1*x1-x0*x1-10*x0-4*x1+60的最优点，收敛精度=0.0014.4 源程序代码及结果： #include "stdio.h"#include "stdlib.h"#include "math.h"double objf(double x)double ff;ff=x0*x0+x1*x1-x0*x1-10*x0-4*x1+60;return(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<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);bi=*(x0+i);elsefor(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*(bi-ai);f1=objf(x0);f2=objf(x1);doif(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;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(double);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<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);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;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;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,2;double ff,x2;ff=powell(p,0.3,0.001,0.0001,2,x);printf("x0=%f,x1=%f,ff=%fn",x0,x1,ff);getchar();5. 复合形法 5.1 复合行法基本思想：在可行域中选取K个设计点 （n+1K2n）作为初始复合形的顶点。比较各顶点目标函数值的大小，去掉目标函数值最大的顶点(称最坏点)，以坏点以外其余各点的中心为映射中心，用坏点的映射点替换该点，构成新的复合形顶点。 反复迭代计算，使复合形不断向最优点移动和收缩，直至收缩到复合形的顶点与形心非常接近，且满足迭代精度要求为止。 5.2 题目：求函数f(x)=(x1-5)*(x1-5)+4*(x2-6)*(x2-6)的最优点，约束条件为g1(x)=64-x1*x1-x2*x20；g2(x)=x2-x1-100；g3(x)=x1-100；收敛精度自定义；5.3 源程序代码及结果：#include <stdio.h> #include <stdlib.h> #include <time.h> #include <math.h> #define E0 1e-5 /*复合形法收敛控制精度*/ double *apply(int,int); /*申请矩阵空间*/ double f(double *); /*目标函数*/ double *g(double *); /*约束函数*/ bool judge(double *); /*可行点的判断*/ int main() int n,k; int i,j,k1; int l; double temporary; double restrain; /*收敛条件*/ double reflect; /*反射系数*/ srand(unsigned)time(NULL); printf("请输入目标函数的维数 n:"); /*输入已知数据*/ scanf("%d",&n); printf("请输入复合形的顶点数 k:"); scanf("%d",&k); double *x=apply(k,n); /*存放复合形顶点*/ double *y=(double *)calloc(k,sizeof(double); /*存放目标函数值*/ double *p=(double *)calloc(3,sizeof(double); /*存放约束函数值*/ double *a=(double *)calloc(n,sizeof(double); /*存放设计变量的下限*/ double *b=(double *)calloc(n,sizeof(double); /*存放设计变量的上限*/ double *x_c=(double *)calloc(n,sizeof(double); /*存放可行点中心*/ double *x_r=(double *)calloc(n,sizeof(double); /*存放最坏点的反射点*/ printf("请输入选定的第一个可行点 x1(包含%d 个数):",n); for(i=0;i<n;i+) scanf("%lf",*x+i); printf("请输入初选变量的下限 a(包含%d 个数):",n); for(i=0;i<n;i+) scanf("%lf",a+i); printf("请输入初选变量的上限 b(包含%d 个数):",n); for(i=0;i<n;i+) scanf("%lf",b+i); printf("输出输入结果为:nn=%d,k=%d,x1=(",n,k); /*输出已知数据*/ for(i=0;i<n-1;i+) printf("%.5lf ",*(*x+i); printf("%.5lf)na=(",*(*x+n-1);for(i=0;i<n-1;i+) printf("%f ",*(a+i); printf("%.5lf),b=(",*(a+n-1); for(i=0;i<n-1;i+) printf("%f ",*(b+i); printf("%.5lf)n",*(b+n-1); L1: for(i=1;i<k;i+) /*随机得到其余(k-1)个可行点*/for(j=0;j<n;j+) *(*(x+i)+j)=*(a+j)+(double)(rand()%10000)/10000*(*(b+j)-*(a+j); l=1; for(i=1;i<k;i+) /*找出可行点的个数 l，并把可行点放在前 l 个位置上*/ if(judge(*(x+i)for(j=1;j<k;j+) if(!judge(*(x+j)for(k1=0;k1<n;k1+) temporary=*(*(x+i)+k1); *(*(x+i)+k1)=*(*(x+j)+k1);*(*(x+j)+k1)=temporary; break;l+; for(i=0;i<l-1;i+) /*把前 l 个可行点按目标函数值从大到小排序*/for(j=i+1;j<l;j+)if(f(*(x+i)<f(*(x+j)for(k1=0;k1<n;k1+)temporary=*(*(x+i)+k1);*(*(x+i)+k1)=*(*(x+j)+k1);*(*(x+j)+k1)=temporary;for(i=0;i<n;i+) /*求可行点中心*/ *(x_c+i)=0;for(i=0;i<l;i+) for(j=0;j<n;j+) *(x_c+j)+=*(*(x+i)+j);for(i=0;i<n;i+) *(x_c+i)/=l;if(!judge(x_c) /*判断可行点中心是否可行*/ for(i=0;i<n;i+)*(a+i)=*(*(x+l-1)+i);*(b+i)=*(x_c+i);goto L1;else for(i=l;i<k;i+) /*将不可行点可行化*/ dofor(j=0;j<n;j+) *(*(x+i)+j)=*(x_c+j)+0.5*(*(*(x+i)+j)-*(x_c+j);while(!judge(*(x+i);L2:for(i=0;i<k-1;i+) /*将可行点按目标函数值从大到小排序*/ for(j=i+1;j<k;j+) if(f(*(x+i)<f(*(x+j) for(k1=0;k1<n;k1+) temporary=*(*(x+i)+k1);*(*(x+i)+k1)=*(*(x+j)+k1); *(*(x+j)+k1)=temporary;restrain=0; /*求收敛条件*/ for(i=0;i<k;i+)restrain+=(f(*(x+i)-f(*(x+k-1)*(f(*(x+i)-f(*(x+k-1);restrain=sqrt(1.0/(k-1)*restrain); if(restrain<E0) /*判断收敛条件*/ printf("n 求得约束最优点为:( ");for(i=0;i<n;i+) printf("%.5f ",*(*(x+k-1)+i);printf(")n 目标函数的最优解为:%.5fn",f(*(x+k-1); return 0;elseL3:for(i=0;i<n;i+) /*计算除去最坏点*x 外的(k-1)个顶点的中心*/ *(x_c+i)=0; for(i=1;i<k;i+)for(j=0;j<n;j+)*(x_c+j)+=*(*(x+i)+j); for(i=0;i<n;i+) *(x_c+i)/=k-1;reflect=1.3; L4:for(i=0;i<n;i+) /*求反射点*/ *(x_r+i)=*(x_c+i)+reflect*(*(x_c+i)-*(*x+i);if(!judge(x_r) reflect*=0.5;goto L4;else if(f(x_r)<f(*x) for(i=0;i<n;i+) *(*x+i)=*(x_r+i); goto L2; else if(reflect<=1e-10) for(i=0;i<n;i+) *(*x+i)=*(*(x+1)+i);goto L3; else reflect*=0.5;goto L4; double *apply(int row,int col) /*申请矩阵空间*/ int i; double *x=(double*)calloc(row*col,sizeof(double); double *y=(double *)calloc(row,sizeof(double *); if(!x | !y) printf("内存分配失败!"); exit(1); for(i=0;i<row;i+) *(y+i)=x+i*col; return y; double f(double *x) /*目标函数*/ return (*x-5)*(*x-5)+4*(*(x+1)-6)*(*(x+1)-6); double *g(double *x) /*约束函数*/ double *p=(double *)calloc(3,sizeof(double); if(!p) printf("内存分配失败!"); exit(1); *p=64-(*x)*(*x)-(*(x+1)*(*(x+1); *(p+1)=*(x+1)-*x-10; *(p+2)=*x-10; return p; bool judge(double *x) /*可行点的判断*/ int i; double *p=(double *)calloc(3,sizeof(double); p=g(x); for(i=0;i<3;i+) if (*(p+i)>0) break; if(i=3) return true; else return false; 6. 外点惩罚函数法6.1解题思路：外点法是从可行域的外部构造一个点序列去逼近原约束问题的最优解。外点法可以用来求解含不等式和等式约束的优化问题。外点惩罚函数的形式为： 6.2 流程框图： 6.3 题目：求函数f(x)=(x1-5)*(x1-5)+4*(x2-6)*(x2-6)的最优点，约束条件：g1(x)=64-x1*x1-x2*x20；g2(x)=x2-x1-100；g3(x)=x1-100；收敛精度=0.00001；6.4 源程序代码及结果： #include <stdio.h>#include<iostream.h>#include<math.h>double lamta10=0, 1.0 ,0 ,0 ,0 ,1 ,0 ,0 ,0 ,1;/鲍威尔方法初始化方向，线性无关double lamta13=0, 0 , 0;/暂存新的搜索方向double x14=0, 0 ,0, 0 ;/x1到x3用于存储各共轭方向的点double x24=0, 0 ,0, 0 ;double x34=0, 0 ,0, 0 ;double x44=0, 0 ,0, 0 ;/x4用于中间判断double x54=0, 0 ,0, 0 ;/x5用存放于更换方向后产生的新点int m=0;/标志double x_4=0, 0, 0, 0;/暂存鲍威尔最优解double x04=0, 2, 2 , 2;/初值 double c=10;/递减系数double e=0.00001;/精度控制double r0=1;/初始惩罚因子double r=1;/函数声明部分void Powell(double r); /鲍威尔方法函数double fxy(double x1,double x2,double x3,double r); /待求函数double ysearch(double x); /一维搜索的目标函数void search(double &a,double &b,double h); /区间搜索double yellowcut(double &a,double &b); /黄金分割void sort(double *p,int size);/选择法排序void main() /约束优化方法主函数入口 cout<<"请输入精度"<<endl; cin>>e;changyan:Powell(r); double cmpare4; int flag1=0; for (int i=1;i<=3;i+) cmparei=x_i-x0i; if (fabs(cmparei)<e) flag1+; if (flag1=3) printf("x1=%lf x2=%lfn",x_1,x_2);/ cout<<"最优解为："<<"x1="<<x_1<<" "<<"x2="<<x_2<<" "<<"x3="<<x_3<<endl ; cout<<"最小值为"<<fxy(x_1,x_2,x_3,r)<<endl; else for (int j=1;j<=3;j+) x0j=x_j; r=c*r; goto changyan;/子函数定义部分double fxy(double x1,double x2,double x3,double r)/待求函数 double m,n,p; m=(64-x1*x1-x2*x2)>0)?(64-x1*x1-x2*x2):0; n=(x2-x1-10)>0)?(x2-x1-10):0; p=(x1-10)>0)?(x1-10):0; return /惩罚函数(x1-5)*(x1-5)+4*(x2-6)*(x2-6)+r*(m*m+n*n+p*p)+r*(x3*x3);void Powell(double r) /鲍威尔方法函数定义 double det=0.0001; /迭代精度 int k;my1: for (k=1;k<=3;k+) m=3*k-2; double a=0,b=0,xo=0; search(a,b,1); /完成区间搜索 double temp; temp=yellowcut(a,b);/黄金分割法 int n=3*k-2; for (int i=1;i<=3;i+) switch (k) case 1:x1i=x0i+temp*lamtan+;break; case 2:x2i=x1i+temp*lamtan+;break; case 3:x3i=x2i+temp*lamtan+;break; default :break; double cmp4; int flag=0; for (int i=1;i<=3;i+) cmpi=x3i-x0i; if (fabs(cmpi)<det) flag+; if (flag=3) /找到最优解 x_1=x31; x_2=x32; x_3=x33; else double fy4; fy0=fxy(x01,x02,x03,r); fy1=fxy(x11,x12,x13,r); fy2=fxy(x21,x22,x23,r); fy3=fxy(x31,x32,x33,r); double fyy3; for (int ii=0;ii<3;ii+) fyyii=fyii-fyii+1; sort(fyy,3); for (ii=1;ii<=3;ii+) x4ii=2*x3ii-x0ii; double f0,f3,f4; f0=fy0; f3=fy3; f4=fxy(x41,x42,x43,r); if (f0+f4-2*f3)/2>=fyy2) if (f3<f4) for (int t=1;t<=3;t+) x0t=x3t; else for (int t=1;t<=3;t+) x0t=x4t; goto my1; else for (int t=0;t<3;t+) lamta1t=x3t+1-x0t+1; m=0; /switch 标志！ double aa=0,bb=0; search(aa,bb,1); double temp1; temp1=yellowcut(aa,bb); for (int i=1;i<=3;i+) x5i=x3i+temp1*lamta1i-1; for (i=1;i<=3;i+) x0i=x5i; for (i=1;i<=6;i+) lamtai=lamtai+3; for (i=1;i<=3;i+) lamta6+i=lamta1i-1; goto my1;double ysearch(double x) /一维搜索的目标函数switch (m)case 1: return fxy(x01+x*lamtam,x02+x*lamtam+1,x03+x*lamtam+2,r);break;case 4: return fxy(x11+x*lamtam,x12+x*lamtam+1,x13+x*lamtam+2,r);break;case 7: return fxy(x21+x*lamtam,x22+x*lamtam+1,x23+x*lamtam+2,r);break;case 0: return fxy(x31+x*lamta10,x32+x*lamta11,x33+x*lamta12,r);break;/