机械优化设计黄金分割法 外推法

.

郑州大学

机械优化设计部分程序

1.外推法

2.黄金分割法

3.二次插值法

4.坐标轮换法

5.随机方向法

6.四杆机构优化设计

.

1.外推法

源程序：

#include

#include

#define R 0.01

double fun(double x)

{ double m;

m=x*x-10*x+36;

return m;

}

void main()

{

double h0=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

#define f(x) x*x*x*x-5*x*x*x+4*x*x-6*x+60 double hj(double *a,double *b,double e,int *n) {

double x1,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);

}

return s;

}

void main()

{

double s,a,b,e,m;

int n=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.279 629,n=21

3.二次插值法

源程序：

#include

#include

int main(void)

{

double a1,a2,a3,ap,y1,y2,y3,yp,c1,c2,m; double j[3];

int i,h=1;

void finding(double a[3]);

finding(j);

a1=j[0];

a2=j[1];

a3=j[2];

m=0.001;

double f(double x);

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;

else if((ap-a2)*h>0)

{

if(y2>=yp){

a1=a2;y1=y2;

a2=ap;y2=yp;}

else{

a3=ap;y3=yp;}

}

else if(y2>=yp){

a3=a2;y3=y2;

a2=ap;y2=yp;}

else{a1=ap;y1=yp;}

}

double x,y;

if(y2<=yp){

x=a2;y=y2;}

else{

x=ap;y=yp;}

printf("a*=%f\n",x);

printf("y*=%f\n",y); return 0;

}

double f(double x)

{

double y;

y=x*x-10*x+36;

return y;

}

void finding(double a[3]) {

int h,i;

double y[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

float fun1(float x,float a,float b)

{float y;

y=x+a*b;

return y;

}

float fun2(float x,float y)

{float z;

z=4*(x-5)*(x-5)+(y-6)*(y-6);

return z;

}

main()

{float d[100][3],x[100][3],xx[3],ax[100][3]; float a1,a2,a3,h,t,y1,y2,y3,e,a,b,l,fi;

int i,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(y3

for(;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

float f(float x,float y)

{

float z;

z=(x-2)*(x-2)+(y-1)*(y-1);

return z;

}

float g1(float x,float y)

{

float z;

z=x*x-y;

return z;

}

float g2(float x,float y)

{

float z;

z=x+y-2;

return z;

}

void main()

{

int i,j;

float

k=8,c=0.000001,a0=-3,b0=3,a1=-3,b1=3;

float

x[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(fx

else break;

}

do

.

{

x[0]=x[0]-h*d[0];

x[1]=x[1]-h*d[1];

h=0.7*h;

if(h

break;

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

#define Pai 3.1415926

int g(float l1,float l2)

{

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);

}

float fun(float x0[2])

{

floatf,

a[31],b[31],r[31],p[31],q[31],w[31],x1[2];

int i;

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*P ai);

f=f+(Pai/60)*(q[i]-w[i])*(q[i]-w[i])*(p[i]-p [i-1]);

}

return f;

}

void main()

{

float

a,q,f,fl,f0,l[2],z[2],d0[100],d1[100],x[2],xi[2],f x,m0,m1,e;

int i,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(f

fl=f;

else break;}

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];

goto ss;

}

}

运行过程及结果：

输入精度0.001

最优解为

x1*=4.161386

x2*=2.311257

fx=0.000021