机械优化设计程序
今天终于交了机械优化设计作业，程序贴出来，那位要用就拿去吧，资源共享了，O(∩_∩)O下面是利用外点惩罚函数(Powell)法求解三维目标函数最优解与最优值的程序，用C++编写。
#include
#include
double lamta={0, 1.0 ,0 ,0 ,0 ,1 ,0 ,0 ,0 ,1};//鲍威尔方法初始化方向，线性无关
double lamta1={0, 0 , 0};//暂存新的搜索方向
double x1={0, 0 ,0, 0 };//x1到x3用于存储各共轭方向的点
double x2={0, 0 ,0, 0 };
double x3={0, 0 ,0, 0 };
double x4={0, 0 ,0, 0 };//x4用于中间判断
double x5={0, 0 ,0, 0 };//x5用存放于更换方向后产生的新点
int m=0;//标志
double x_={0, 0, 0, 0};//暂存鲍威尔最优解
double x0={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;
int flag1=0;
for (int i=1;i<=3;i++)
{
cmpare=x_-x0;
if (fabs(cmpare)
{flag1++;}
}
if (flag1==3)
{
cout<<"最优解为："<<"x1="<
}
else
{
for (int j=1;j<=3;j++)
{
x0=x_;
}
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=x0+temp*lamta;break;
case 2:x2=x1+temp*lamta;break;
case 3:x3=x2+temp*lamta;break;
default :break;
}
}
}
double cmp;
int flag=0;
for (int i=1;i<=3;i++)
{
cmp=x3-x0;
if (fabs(cmp)
{flag++;}}
if (flag==3) //找到最优解
{
x_=x3;
x_=x3;
x_=x3;
}
else
{
double fy;
fy=fxy(x0,x0,x0,r);
fy=fxy(x1,x1,x1,r);
fy=fxy(x2,x2,x2,r);
fy=fxy(x3,x3,x3,r); double fyy;
for (int ii=0;ii<3;ii++)
{fyy=fy-fy;}
sort(fyy,3);
for (ii=1;ii<=3;ii++)
{x4=2*x3-x0;}
double f0,f3,f4;
f0=fy;
f3=fy;
f4=fxy(x4,x4,x4,r);
if ((f0+f4-2*f3)/2>=fyy)
{
if (f3
{
for (int t=1;t<=3;t++)
{x0=x3;}
}
else
{
for (int t=1;t<=3;t++)
{x0=x4;
}}
goto my1;
}
else{
for (int t=0;t<3;t++)
{lamta1=x3-x0;}
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=x3+temp1*lamta1;}
for (i=1;i<=3;i++)
{x0=x5;}
for (i=1;i<=6;i++)
{lamta=lamta;}
for (i=1;i<=3;i++)
{
lamta=lamta1;}
goto my1;
}}}
double ysearch(double x) //一维搜索的目标函数
{
switch (m)
{
case 1: return fxy(x0+x*lamta,x0+x*lamta,x0+x*lamta,r);break;
case 4: return fxy(x1+x*lamta,x1+x*lamta,x1+x*lamta,r);break;
case 7: return fxy(x2+x*lamta,x2+x*lamta,x2+x*lamta,r);break;
case 0: return fxy(x3+x*lamta1,x3+x*lamta1,x3+x*lamta1,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;
double xo;
d=a+0.618*(b-a);
fd=ysearch(d);
label2: if (fc<=fd)
{b=d;
d=c;
fd=fc;
c=a+0.382*(b-a);
fc=ysearch(c);}
else
{a=c;
c=d;
fc=fd;
d=a+0.618*(b-a);
fd=ysearch(d);}
if ((b-a)<=e)
{xo=(a+b)/2;}
else
goto label2;
return xo;
}
void sort(double *p,int size){//选择法排序
int i,j;
double k;
for(i=0;i
for(j=i+1;j
if(*(p+i)>*(p+j)){k=*(p+i);*(p+i)=*(p+j);*(p+j)=k;}
}
转自：http://blog.sina.com.cn/s/blog_49d083e50100cojl.html
页:
[1]