机械优化设计黄金分割法 外推法.docx
《机械优化设计黄金分割法 外推法.docx》由会员分享,可在线阅读,更多相关《机械优化设计黄金分割法 外推法.docx(23页珍藏版)》请在冰豆网上搜索。
机械优化设计黄金分割法外推法
大学
机械优化设计部分程序
1.外推法
2.黄金分割法
3.二次插值法
4.坐标轮换法
5.随机方向法
6.四杆机构优化设计
1.外推法
源程序:
#include
#include
#defineR0.01
doublefun(doublex)
{doublem;
m=x*x-10*x+36;
returnm;
}
voidmain()
{
doubleh0=R,y1,y2,y3,x1,x2,x3,h;
x1=0;h=h0;x2=h;
y1=fun(x1);y2=fun(x2);
if(y2>y1)
{h=-h;
x3=x1;
y3=y1;
x1=x2;
y1=y2;
x2=x3;
y2=y3;
}
x3=x2+h;y3=fun(x3);
while(y3{h*=2.0;
x1=x2;
y1=y2;
x2=x3;
y2=y3;
x3=x2+h;
y3=fun(x3);
}
printf("fun(%f)=%f,fun(%f)=%f,fun(%f)=%f\n",x1,y1,x2,y2,x3,y3);
}
运行过程及结果:
fun(2.560000)=16.953600,
fun(5.120000)=11.014400,
fun(10.240000)=38.457600
2.黄金分割法
源程序:
#include
#include
#definef(x)x*x*x*x-5*x*x*x+4*x*x-6*x+60
doublehj(double*a,double*b,doublee,int*n)
{
doublex1,x2,s;
if(fabs((*b-*a)/(*b))<=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=hj(a,b,e,n);
}
returns;
}
voidmain()
{
doubles,a,b,e,m;
intn=0;
printf("输入a,b值和精度e值\n");
scanf("%lf%lf%lf",&a,&b,&e);
s=hj(&a,&b,e,&n);
m=(a+b)/2;
printf("a=%lf,b=%lf,s=%lf,m=%lf,n=%d\n",a,b,s,m,n);
}
运行过程及结果:
输入a,b值和精度e值
-3
5
0.0001
a=3.279466,b=3.279793,s=22.659008,m=3.279629,n=21
3.二次插值法
源程序:
#include
#include
intmain(void)
{
doublea1,a2,a3,ap,y1,y2,y3,yp,c1,c2,m;
doublej[3];
inti,h=1;
voidfinding(doublea[3]);
finding(j);
a1=j[0];
a2=j[1];
a3=j[2];
m=0.001;
doublef(doublex);
y1=f(a1);
y2=f(a2);
y3=f(a3);
for(i=1;1>=1;i++)
{
c1=(y3-y1)/(a3-a1);
c2=((y2-y1)/(a2-a1)-c1)/(a2-a3);
ap=0.5*(a1+a3-c1/c2);
yp=f(ap);
if(fabs((y2-yp)/y2)break;
elseif((ap-a2)*h>0)
{
if(y2>=yp){
a1=a2;y1=y2;
a2=ap;y2=yp;}
else{
a3=ap;y3=yp;}
}
elseif(y2>=yp){
a3=a2;y3=y2;
a2=ap;y2=yp;}
else{a1=ap;y1=yp;}
}
doublex,y;
if(y2<=yp){
x=a2;y=y2;}
else{
x=ap;y=yp;}
printf("a*=%f\n",x);
printf("y*=%f\n",y);
return0;
}
doublef(doublex)
{
doubley;
y=x*x-10*x+36;
returny;
}
voidfinding(doublea[3])
{
inth,i;
doubley[3];
a[0]=0;
h=1;
a[1]=h;
y[0]=f(a[0]);y[1]=f(a[1]);
if(y[1]>y[0])
{
h=-h;
a[2]=a[0];y[2]=y[0];
do{
a[0]=a[1];a[1]=a[2];
y[0]=y[1];y[1]=y[2];
a[2]=a[1]+h;y[2]=f(a[2]);
h=2*h;
}while(y[2]}
else{
for(i=1;i>=1;i++){
a[2]=a[1]+h;y[2]=f(a[2]);
if(y[2]>=y[1])
break;
h=2*h;
a[0]=a[1];y[0]=y[1];
a[1]=a[2];y[1]=y[2];}
}
return;
}
运行过程及结果:
a*=5.000000
y*=11.000000
4.坐标轮换法
源程序:
#include
#include
#include
floatfun1(floatx,floata,floatb)
{floaty;
y=x+a*b;
returny;
}
floatfun2(floatx,floaty)
{floatz;
z=4*(x-5)*(x-5)+(y-6)*(y-6);
returnz;
}
main()
{floatd[100][3],x[100][3],xx[3],ax[100][3];
floata1,a2,a3,h,t,y1,y2,y3,e,a,b,l,fi;
inti,k;
printf("输入初始点坐标\n");
scanf("%f%f",&x[0][1],&x[0][2]);
e=0.000001;
l=0.618;
x[2][1]=x[0][1];
x[2][2]=x[0][2];
k=0;
k--;
do
{x[0][1]=x[2][1];
x[0][2]=x[2][2];
k++;
for(i=1;i<=2;i++)
{
if(i==1)
{d[i][1]=1;
d[i][2]=0;
}
else
{d[i][1]=0;
d[i][2]=1;
}
h=0.1;
a1=0;
a2=h;
x[i][1]=fun1(x[i-1][1],d[i][1],a1);
x[i][2]=fun1(x[i-1][2],d[i][2],a1);
y1=fun2(x[i][1],x[i][2]);
x[i][1]=fun1(x[i-1][1],d[i][1],a2);
x[i][2]=fun1(x[i-1][2],d[i][2],a2);
y2=fun2(x[i][1],x[i][2]);
if(y2>y1)
{h=-h;
a3=a1;
y3=y1;
a1=a2;
a2=a3;
y1=y2;
y2=y3;
}
a3=a2+h;
x[i][1]=fun1(x[i-1][1],d[i][1],a3);
x[i][2]=fun1(x[i-1][2],d[i][2],a3);
y3=fun2(x[i][1],x[i][2]);
do
{a1=a2;
y1=y2;
a2=a3;
y2=y3;
a3=a2+h;
x[i][1]=fun1(x[i-1][1],d[i][1],a3);
x[i][2]=fun1(x[i-1][2],d[i][2],a3);
y3=fun2(x[i][1],x[i][2]);
}
while(y3for(;a1>a3;)
{t=a3;
a3=a1;
a1=t;
t=y1;
y3=y1;
y1=t;
}
a=a1;
b=a3;
a1=b-l*(b-a);
a2=a+l*(b-a);
x[i][1]=fun1(x[i-1][1],d[i][1],a1);
x[i][2]=fun1(x[i-1][2],d[i][2],a1);
y1=fun2(x[i][1],x[i][2]);
x[i][1]=fun1(x[i-1][1],d[i][1],a2);
x[i][2]=fun1(x[i-1][2],d[i][2],a2);
y2=fun2(x[i][1],x[i][2]);
if(b<1e-3)
{
for(;fabs(b-a)>e;)
{
if(y1>=y2)
{a=a1;
a1=a2;
y1=y2;
a2=a+l*(b-a);
x[i][1]=fun1(x[i-1][1],d[i][1],a2);
x[i][2]=fun1(x[i-1][2],d[i][2],a2);
y2=fun2(x[i][1],x[i][2]);
}
else
{b=a2;
a2=a1;
y2=y1;
a1=b-l*(b-a);
x[i][1]=fun1(x[i-1][1],d[i][1],a1);
x[i][2]=fun1(x[i-1][2],d[i][2],a1);
y1=fun2(x[i][1],x[i][2]);
}
}
}
else
{
for(;fabs((b-a)/b)>=e||fabs((y2-y1)/y2)>=e;)
{
if(y1>=y2)
{a=a1;
a1=a2;
y1=y2;
a2=a+l*(b-a);
x[i][1]=fun1(x[i-1][1],d[i][1],a2);
x[i][2]=fun1(x[i-1][2],d[i][2],a2);
y2=fun2(x[i][1],x[i][2]);
}
else
{b=a2;
a2=a1;
y2=y1;
a1=b-l*(b-a);
x[i][1]=fun1(x[i-1][1],d[i][1],a1);
x[i][2]=fun1(x[i-1][2],d[i][2],a1);
y1=fun2(x[i][1],x[i][2]);
}
}
}
ax[k][i]=0.5*(a+b);
x[i][1]=fun1(x[i-1][1],d[i][1],ax[k][i]);
x[i][2]=fun1(x[i-1][2],d[i][2],ax[k][i]);
}
}
while(sqrt(pow((x[2][1]-x[0][1]),2)+pow((x[2][2]-x[0][2]),2))>=1e-6);
xx[1]=x[2][1];
xx[2]=x[2][2];
fi=fun2(xx[1],xx[2]);
printf("最优解为\nx1*=%f\nx2*=%f\n
f*=%f\nk=%d\n",xx[1],xx[2],fi,k);
}
运行过程及结果:
输入初始点坐标
8
9
最优解为
x1*=5.000000
x2*=6.000000
f*=0.000000
k=2
5.随机方向法
源程序:
#include
#include
#include
floatf(floatx,floaty)
{
floatz;
z=(x-2)*(x-2)+(y-1)*(y-1);
returnz;
}
floatg1(floatx,floaty)
{
floatz;
z=x*x-y;
returnz;
}
floatg2(floatx,floaty)
{
floatz;
z=x+y-2;
returnz;
}
voidmain()
{
inti,j;
floatk=8,c=0.000001,a0=-3,b0=3,a1=-3,b1=3;
floatx[10],x0[10],xl[10],e[10],r[10],d[10],h,fl,f0,fx;
while(g1(x0[0],x0[1])>0||g2(x0[0],x0[1])>0)
{
x0[0]=a0+(rand()/32767.00)*(b0-a0);
x0[1]=a1+(rand()/32767.00)*(b1-a1);
}
fl=f(x0[0],x0[1]);
f0=f(x0[0],x0[1]);
while
(1)
{
h=0.01;
j=1;
r[0]=-1+(rand()/32767.00)*(1-(-1));
r[1]=-1+(rand()/32767.00)*(1-(-1));
e[0]=r[0]/sqrt(r[0]*r[0]+r[1]*r[1]);
e[1]=r[1]/sqrt(r[0]*r[0]+r[1]*r[1]);
x[0]=x0[0]+h*e[0];
x[1]=x0[1]+h*e[1];
if(g1(x[0],x[1])<=0&&g2(x[0],x[1])<=0)
{
fx=f(x[0],x[1]);
if(fx{
fl=fx;
for(i=0;i<2;i++)
{d[i]=e[i];xl[i]=x[i];}
}
}
while(j<=k)
{
j++;
r[0]=-1+(rand()/32767.00)*(1-(-1));
r[1]=-1+(rand()/32767.00)*(1-(-1));
e[0]=r[0]/sqrt(r[0]*r[0]+r[1]*r[1]);
e[1]=r[1]/sqrt(r[0]*r[0]+r[1]*r[1]);
x[0]=x0[0]+h*e[0];
x[1]=x0[1]+h*e[1];
if(g1(x[0],x[1])<=0&&g2(x[0],x[1])<=0)
{
fx=f(x[0],x[1]);
if(fx{
fl=fx;
for(i=0;i<2;i++)
{d[i]=e[i];xl[i]=x[i];}
}
}
}
x[0]=xl[0];
x[1]=xl[1];
while
(1)
{
h=1.3*h;
x[0]=x[0]+h*d[0];
x[1]=x[1]+h*d[1];
if(g1(x[0],x[1])>0||g2(x[0],x[1])>0)
break;
fx=f(x[0],x[1]);
if(fxelsebreak;
}
do
{
x[0]=x[0]-h*d[0];
x[1]=x[1]-h*d[1];
h=0.7*h;
if(hbreak;
x[0]=x[0]+h*d[0];
x[1]=x[1]+h*d[1];
if(g1(x[0],x[1])>0||g2(x[0],x[1])>0)
continue;
fx=f(x[0],x[1]);
}
while(fx>=fl);
if(fabs((f0-fx)/f0)>=c)
{
x0[0]=x[0];
x0[1]=x[1];
fl=fx;
f0=fx;
}
else
break;
}
printf("输出最优解为\nx1*=%f,x2*=%f,
y*=%f\n",x[0],x[1],fx);
}
运行过程及结果:
输出最优解为
x1*=0.995421,x2*=1.004521,y*=1.009200
6.四杆机构优化设计
源程序:
#include
#include
#include
#definePai3.1415926
intg(floatl1,floatl2)
{
if((-l1<=0)
&&(-l2<=0)
&&(6-l1-l2<=0)
&&(1-l2-4<=0)
&&(l2-l1-4<=0)
&&(l1*l1+l2*l2-1.414*l1*l2-16<=0)
&&(36-l1*l1-l2*l2-1.414*l1*l2<=0))
return
(1);
else
return(0);
}
floatfun(floatx0[2])
{
floatf,a[31],b[31],r[31],p[31],q[31],w[31],x1[2];
inti;
p[0]=acos(((1+x0[0])*(1+x0[0])-x0[1]*x0[1]+25)/(10+10*x0[0]));
q[0]=acos(((1+x0[0])*(1+x0[0])-x0[1]*x0[1]-25)/(10*x0[1]));
f=0;
for(i=1;i<=30;i++)
{
p[i]=p[0]+(Pai/60)*i;
r[i]=sqrt(26-10*cos(p[i]));
a[i]=acos((r[i]*r[i]+x0[1]*x0[1]-x0[0]*x0[0])/(2*r[i]*x0[1]));
b[i]=acos((r[i]*r[i]+24)/(10*r[i]));
q[i]=Pai-a[i]-b[i];
w[i]=q[0]+(2*(p[i]-p[0])*(p[i]-p[0]))/(3*Pai);
f=f+(Pai/60)*(q[i]-w[i])*(q[i]-w[i])*(p[i]-p[i-1]);
}
returnf;
}
voidmain()
{
floata,q,f,fl,f0,l[2],z[2],d0[100],d1[100],x[2],xi[2],fx,m0,m1,e;
inti,j,n,k;
printf("输入精度");
scanf("%f",&e);
do{z[0]=0+5*(rand()/32767.00);
z[1]=0+5*(rand()/32767.00);
}
while(g(z[0],z[1])==0);
for(i=0;i<=99;i++)
{
d0[i]=-1+2*(rand()/32767.00);
}
for(j=0;j<=99;j++)
{
d1[j]=-1+2*(rand()/32767.00);
}
f0=fun(z);
fl=fun(z);
ss:
a=0.01;
for(i=0,j=0;i<=99&&j<=99;i++,j++)
{
n=1/sqrt((d0[i])*(d0[i])+d1[j]*d1[j]);
d0[i]=n*d0[i];
d1[j]=n*d1[j];
x[0]=z[0]+a*d0[i];
x[1]=z[1]+a*d1[j];
if(g(x[0],x[1])==1)
{
f=fun(x);
if(f{
fl=f;
m0=d0[i];
m1=d1[j];
l[0]=x[0];
l[1]=x[1];
}
}
}
x[0]=l[0];
x[1]=l[1];
do
{a=1.3*a;
x[0]=x[0]+a*m0;
x[1]=x[1]+a*m1;
if(g(z[0],z[1])==0)
break;
f=fun(x);
if(ffl=f;
elsebreak;}
while(g(z[0],z[1])==1);
do{
x[0]=x[0]-a*m0;
x[1]=x[1]-a*m1;
a=0.7*a;
if(a<0.00001)break;
x[0]=x[0]+a*m0;
x[1]=x[1]+a*m1;
if(g(z[0],z[1])==1)
f=fun(x);
}
while(f>=fl);
if(fabs((f0-f)/f0){xi[0]=x[0];xi[1]=x[1];
fx=f;
printf("最优解为\nx1*=%f\nx2*=%f\n
fx=%f\n",xi[0],xi[1],fx);
}
else
{f0=f;
fl=f;
z[0]=x[0];
z[1]=x[1];
gotoss;
}
}
运行过程及结果:
输入精度0.001
最优解为
x1*=4.161386
x2*=2.311257
fx=0.000021