机械优化设计惩罚函数内点法.doc
《机械优化设计惩罚函数内点法.doc》由会员分享,可在线阅读,更多相关《机械优化设计惩罚函数内点法.doc(6页珍藏版)》请在冰豆网上搜索。
#include
#include
#definem12
doublef(doublex[],doubler);
voidjintuifa(doubleab[m][m],intn,doublex0[],doubleh,intij,doublea[],doubleb[],doubler0);
voidhongjinfa(intn,doublea[],doubleb[],doubleflag,doublex[],doubler0);
voidbaoweifa(intn,doublex0[],doubleh,doubleflag,doublea[],doubleb[],doublex[],doubler0);
doublefahansu(doublett)
{
doublety;
if(tt<0)ty=-tt;elsety=0;
returnty;
}
doubleyuanhansu(doublex[])
{
doubles;
//s=x[0]*x[0]+x[1]*x[1];
s=x[0]*x[0]+x[1]*x[1]+x[2]*x[2]+x[3]*x[3];
returns;
}
doublef(doublex[],doubler)
{
doubles,t,t2;
//t=1-x[0];
t=1-x[0];t2=2-x[1];
//s=yuanhansu(x)-r*log(fahansu(t));
s=yuanhansu(x)-r*log(fahansu(t))-r*log(fahansu(t2));
returns;
}
voidjintuifa(doubleab[m][m],intn,doublex0[],doubleh,intij,doublea[],doubleb[],doubler0)
{
inti,j,z;
doublex1[m],x2[m],x3[m],f1,f2,f3;
doubles[m];
for(i=0;i{
s[i]=ab[i][ij];
}
for(i=0;i{
x1[i]=x0[i];
x2[i]=x0[i]+(double)h*s[i];
}
f1=f(x1,r0);
f2=f(x2,r0);
if(f2>=f1)
{
h=(-1)*h;
for(i=0;ix3[i]=x1[i];
f3=f1;
for(i=0;ix1[i]=x2[i];
f1=f2;
for(i=0;ix2[i]=x3[i];
f2=f3;
}
for(i=0;ix3[i]=x2[i]+(double)h*s[i];
f3=f(x3,r0);
while(f3{
h=2*h;
for(i=0;ix1[i]=x2[i];
f1=f2;
for(i=0;ix2[i]=x3[i];
f2=f3;
for(i=0;ix3[i]=x2[i]+(double)h*s[i];
f3=f(x3,r0);
}
for(i=0;i{
if(x1[i]{
a[i]=x1[i];
b[i]=x3[i];
}
else
{
a[i]=x3[i];
b[i]=x1[i];
}
}
}
voidhongjinfa(intn,doublea[],doubleb[],doubleflag,doublex[],doubler0)
{
inti;
doublex1[m],x2[m],f1,f2;
while
(1)
{
for(i=0;ix1[i]=b[i];
f1=f(x1,r0);
for(i=0;ix2[i]=a[i];
f2=f(x2,r0);
if(fabs((f2-f1)/f2)<=flag)
{
for(i=0;ix[i]=(double)(b[i]+a[i])/2.0;
break;
}
else
{
for(i=0;ix1[i]=b[i]-(double)0.618*(b[i]-a[i]);
for(i=0;ix2[i]=a[i]+(double)0.618*(b[i]-a[i]);
if(f(x1,r0)>f(x2,r0))
for(i=0;ia[i]=x1[i];
else
for(i=0;ib[i]=x2[i];
}
}
}
voidbaoweifa(intn,doublex0[],doubleh,doubleflag,doublea[],doubleb[],doublex[],doubler0)
{
inti,j,k,r,sum,p,e;
doublex1[m],x2[m],abc[m],kr[m],f0,f1,f2,fn[m],c[m],q,wo;
k=0;
doubleab[m][m];
for(i=0;i{
for(j=0;j{
if(i==j)
ab[i][j]=1;
else
ab[i][j]=0;
}
}
while
(1)
{
wo=f(x0,r0);
for(p=0;p{
q=f(x0,r0);
jintuifa(ab,n,x0,h,p,a,b,r0);
hongjinfa(n,a,b,flag,x1,r0);
for(e=0;e{
x0[e]=x1[e];
}
fn[p]=(double)(q-f(x1,r0));
}
k=k+1;
for(i=0;ix2[i]=(double)2*x1[i]-x0[i];
for(i=0;ic[i]=(double)x1[i]-x0[i];
r=0;
for(i=0;i{
if(fn[0]<=fn[i])
{
fn[0]=fn[i];
r=i;
}
}
f0=f(x0,r0);
f1=f(x1,r0);
f2=f(x2,r0);
if(f2>=f0||((double)f0-(double)2*f1+(double)f2)*((double)f0-(double)f1-(double)fn[0])*((double)f0-(double)f1-(double)fn[0])>=0.5*fn[0]*(double)(f0-f2)*(double)(f0-f2))
{
if(f1for(i=0;ix0[i]=x1[i];
else
for(i=0;ix0[i]=x2[i];
}
else
{
intij;
ij=n-1;
for(i=0;ikr[i]=x1[i];
for(i=0;i{
for(j=0;j{
if(j>=r)
ab[i][j]=ab[i][j+1];
}
ab[i][n-1]=c[i];
}
jintuifa(ab,n,kr,h,ij,a,b,r0);
hongjinfa(n,a,b,flag,x1,r0);
for(i=0;ix0[i]=x1[i];
}
if((fabs(f(x0,r0)-wo))/f(x0,r0)<=flag)
break;
}
for(i=0;ix[i]=x0[i];
}
intmain()
{
intv,i,n;
doubleh,flag,x0[m],a[m],b[m],x[m];
doublefom,fxo,c,r0;
c=0.5;
r0=100.0;
fom=100;
printf("请输入维数:
\n");
scanf("%d",&n);
printf("请输入初始点:
");
for(i=0;i{
printf("\nx0[%d]=",i);
scanf("%lf"