机械优化设计程序

机械优化设计程序 (2009-04-24 19:30:52)

转载 标签： 机械优化设计 程序 powell 外点惩罚函数法 c 杂谈

今天终于交了机械优化设计作业，程序贴出来，那位要用就拿去吧，资源共享了，O(∩_∩)O 下面是利用外点惩罚函数（Powell）法求解三维目标函数最优解与最优值的程序，用C++编写。

#include #include double lamta[10]={0, 1.0 ,0 ,0 ,0 ,1 ,0 ,0 ,0 ,1};//鲍威尔方法初始化方向，线性无关 double lamta1[3]={0, 0 , 0};//暂存新的搜索方向 double x1[4]={0, 0 ,0, 0 };//x1到x3用于存储各共轭方向的点 double x2[4]={0, 0 ,0, 0 }; double x3[4]={0, 0 ,0, 0 }; double x4[4]={0, 0 ,0, 0 };//x4用于中间判断 double x5[4]={0, 0 ,0, 0 };//x5用存放于更换方向后产生的新点 int m=0;//标志 double x_[4]={0, 0, 0, 0};//暂存鲍威尔最优解 double x0[4]={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<<"请输入精度"< cin>>e; changyan:Powell(r); double cmpare[4]; int flag1=0; for (int i=1;i<=3;i++) { cmpare[i]=x_[i]-x0[i]; if (fabs(cmpare[i]) {flag1++;} } if (flag1==3) { cout<<"最优解为："<<"x1="< } else { for (int j=1;j<=3;j++) { x0[j]=x_[j]; } r=c*r; goto changyan; } } //子函数定义部分 double fxy(double x1,double x2,double x3,double r)//待求函数 { double m,n,p; m=(-x1>0)?(-x1):0; n=(-x2>0)?(-x2):0; p=(-x3>0)?(-x3):0; return //惩罚函数 1000-x1*x1-2*x2*x2-x3*x3-x1*x2-x1*x3+r*(m*m+n*n+p*p)+r*((x1*x1+x2*x2+x3*x3-25)*(x1*x1+x2*x2+x3*x3-25)+(8*x1+14*x2+7*x3-56)*(8*x1+14*x2+7*x3-56)); } 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:x1[i]=x0[i]+temp*lamta[n++];break; case 2:x2[i]=x1[i]+temp*lamta[n++];break; case 3:x3[i]=x2[i]+temp*lamta[n++];break; default :break; } } } double cmp[4]; int flag=0; for (int i=1;i<=3;i++) { cmp[i]=x3[i]-x0[i]; if (fabs(cmp[i]) {flag++;}} if (flag==3) //找到最优解 { x_[1]=x3[1]; x_[2]=x3[2]; x_[3]=x3[3];

} else { double fy[4]; fy[0]=fxy(x0[1],x0[2],x0[3],r); fy[1]=fxy(x1[1],x1[2],x1[3],r); fy[2]=fxy(x2[1],x2[2],x2[3],r); fy[3]=fxy(x3[1],x3[2],x3[3],r); double fyy[3]; for (int ii=0;ii<3;ii++) {fyy[ii]=fy[ii]-fy[ii+1];} sort(fyy,3); for (ii=1;ii<=3;ii++) {x4[ii]=2*x3[ii]-x0[ii];} double f0,f3,f4; f0=fy[0]; f3=fy[3]; f4=fxy(x4[1],x4[2],x4[3],r); if ((f0+f4-2*f3)/2>=fyy[2]) { if (f3 { for (int t=1;t<=3;t++) {x0[t]=x3[t];} } else { for (int t=1;t<=3;t++) {x0[t]=x4[t]; }} goto my1; } else{ for (int t=0;t<3;t++) {lamta1[t]=x3[t+1]-x0[t+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++) {x5[i]=x3[i]+temp1*lamta1[i-1];} for (i=1;i<=3;i++) {x0[i]=x5[i];} for (i=1;i<=6;i++) {lamta[i]=lamta[i+3];} for (i=1;i<=3;i++) { lamta[6+i]=lamta1[i-1];} goto my1; }}}

double ysearch(double x) //一维搜索的目标函数 { switch (m) { case 1: return fxy(x0[1]+x*lamta[m],x0[2]+x*lamta[m+1],x0[3]+x*lamta[m+2],r);break; case 4: return fxy(x1[1]+x*lamta[m],x1[2]+x*lamta[m+1],x1[3]+x*lamta[m+2],r);break; case 7: return fxy(x2[1]+x*lamta[m],x2[2]+x*lamta[m+1],x2[3]+x*lamta[m+2],r);break; case 0: return fxy(x3[1]+x*lamta1[0],x3[2]+x*lamta1[1],x3[3]+x*lamta1[2],r);break;//更改方向后的一维搜索 default:return 0; break; } } void search(double &a,double &b,double h) //区间搜索 {double a1,a2,a3,y1,y2,y3; h=1; a1=a,y1=ysearch(a1); a2=a+h,y2=ysearch(a2); if(y2>=y1){ h=-h,a3=a1,y3=y1; a1=a2,y1=y2,a2=a3,y2=y3;} a3=a2+h,y3=ysearch(a3); while(y3<=y2){ h=2*h; a1=a2,y1=y2,a2=a3,y2=y3; a3=a2+h,y3=ysearch(a3); } if(h<0)a=a3,b=a1; else a=a1,b=a3;} double yellowcut(double &a,double &b){ double e; //黄金分割法求解 e=0.001; double c,fc; c=a+0.382*(b-a); fc=ysearch(c); double d,fd;