机械优化设计黄金分割法 外推法.docx
《机械优化设计黄金分割法 外推法.docx》由会员分享，可在线阅读，更多相关《机械优化设计黄金分割法 外推法.docx（23页珍藏版）》请在冰豆网上搜索。
机械优化设计黄金分割法外推法
郑州大学
机械优化设计部分程序
1.外推法
2.黄金分割法
3.二次插值法
4.坐标轮换法
5.随机方向法
6.四杆机构优化设计
1.外推法
源程序：
#include<>
#include<>
#defineR
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*=;
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=,
fun=,
fun=
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-*a）;
x2=*a+*（*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
a=,b=,s=,m=,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=;
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=*（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*=
y*=
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=;
l=;
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=;
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]=*（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*=
x2*=
f*=
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=,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（）/*（b0-a0）;
x0[1]=a1+（rand（）/*（b1-a1）;
}
fl=f（x0[0],x0[1]）;
f0=f（x0[0],x0[1]）;
while
（1）
{
h=;
j=1;
r[0]=-1+（rand（）/*（1-（-1））;
r[1]=-1+（rand（）/*（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（）/*（1-（-1））;
r[1]=-1+（rand（）/*（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=*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=*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*=,x2*=,y*=
6.四杆机构优化设计
源程序：
#include<>
#include<>
#include<>
#definePai
intg（floatl1,floatl2）
{
if（（-l1<=0）
&&（-l2<=0）
&&（6-l1-l2<=0）
&&（1-l2-4<=0）
&&（l2-l1-4<=0）
&&（l1*l1+l2**l1*l2-16<=0）
&&（36-l1*l1-l2**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（）/;
z[1]=0+5*（rand（）/;
}
while（g（z[0],z[1]）==0）;
for（i=0;i<=99;i++）
{
d0[i]=-1+2*（rand（）/;
}
for（j=0;j<=99;j++）
{
d1[j]=-1+2*（rand（）/;
}
f0=fun（z）;
fl=fun（z）;
ss:
a=;
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=*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=*a;
if（ax[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;
}
}
运行过程及结果：
输入精度
最优解为
x1*=
x2*=
fx=
升级会员 