论坛 › 算法及编程语言 › 内点惩罚函数法子程序

风花雪月 发表于 2007-5-10 01:28

内点惩罚函数法子程序

本程序包含6个C文件
dfpopt.c
funct.c
hjfgf.c
jtf.c
mldfhs1.c
powell.c

本程序由heizi友情提供

[ 本帖最后由 风花雪月 于 2008-1-4 08:55 编辑 ]

风花雪月 发表于 2007-5-10 01:29

dfpopt.c
#include "hjfgf.c"
#include "stdlib.h"
void gradient(double x[],double g[],int n)
{int i;
double af,f1,f2,dltx=0.000001;
for(i=0;i<n;i++)
g=0;
f1=objf(x);
for(i=0;i<n;i++)
{af=*(x+i);
   *(x+i)=af+dltx;
   f2=objf(x);
   g=(f2-f1)/dltx;
   *(x+i)=af;
}
}
double oneoptim(double x0[],double s[],double h0,double epsg,int n,double x[])
{double *a,*b,ff;
a=(double *)malloc(n*sizeof (double ));
b=(double *)malloc(n*sizeof (double ));
jtf(x0,h0,s,n,a,b);
ff=gold(a,b,epsg,n,x);
free(a);
free(b);
return(ff);
}
double dfpopt(double xx[],double h0,double eps,double epsg,int n)
{int i,j,k;
double ae,zcc;
double *s,*x,*ay,*df,*zd,*zc,*zh;
s=(double *)malloc(n*sizeof (double ));
x=(double *)malloc(n*sizeof (double ));
for(i=0;i<2;i++)
{ay=(double*)malloc(n*sizeof(double));
   df=(double*)malloc(n*sizeof(double));
   zd=(double*)malloc(n*sizeof(double));
   zc=(double*)malloc(n*sizeof(double));
   zh=(double*)malloc(n*n*sizeof(double));
}
for(i=0;i<n;i++)
*(ay+i)=xx;
L1:
   k=0;
   for (i=0;i<n;i++)
   for (j=0;j<n;j++)
    {*(zh+i*n+j)=0;
   if(i==j)
      *(zh+i*n+j)=1;
    }
for(i=0;i<n;i++)
{*(x+i)=*(ay+i);
   *(ay+i)=*(ay+i);
}
gradient(x,df,n);
for(i=0;i<n;i++)
{*(s+i)=0;
   *(df+i)=*(df+i);
   for(j=0;j<n;j++)
   *(s+i)=*(s+i)-(*(zh+i*n+j))*(*(df+j));
}
L2:
oneoptim(x,s,h0,epsg,n,ay);
for(i=0;i<n;i++)
   *(x+i)=*(ay+i);
gradient(x,df,n);
ae=0;
for(i=0;i<n;i++)
   ae=ae+(*(df+i))*(*(df+i));
if(ae<=eps)
   {for(i=0;i<n;i++)
    *(xx+i)=*(x+i);
    free(s);
    free(x);
    for(i=0;i<2;i++)
    {free(ay);
   free(df);
   free(zd);
   free(zc);
   free(zh);
    }
    return(objf(xx));
   }
   if(k==n) goto L1;
   zcc=0;
   for(i=0;i<n;i++)
   {*(zd+i)=*(ay+i)-(*(ay+i));
    *(zd+i)=*(df+i)-(*(df+i));
    *(df+i)=*(df+i);
    zcc=zcc+(*(zd+i))*(*(zd+i));
   }
   for(i=0;i<n;i++)
    for(j=0;j<n;j++)
    *(zh+i*n+j)=*(zd+i)*(*(zd+j))/zcc;
   for(i=0;i<n;i++)
    { *(zc+i)=0;
      for(j=0;j<n;j++)
       *(zc+i)=*(zc+i)+(*(zd+j))*(*(zh+j*n+i));
    }
   zcc=0;
    for(i=0;i<n;i++)
   zcc=zcc+(*(zc+i))*(*(zd+i));
   for(i=0;i<n;i++)
   {*(zc+i)=0;
      *(zc+i)=0;
      for (j=0;j<n;j++)
       {*(zc+i)=*(zc+i)+(*(zh+i*n+j))*(*(zd+j));
      *(zc+i)=*(zc+i)+(*(zh+j))*(*(zh+j*n+i));
       }
   }
    for(i=0;i<n;i++)
    for(j=0;j<n;j++)
    *(zh+i*n+j)=*(zh+i*n+j)+(*(zh+i*n+j))-(*(zc+i))*(*(zc+j))/zcc;
    for(i=0;i<n;i++)
   {*(s+i)=0;
      for(j=0;j<n;j++)
       *(s+i)=*(s+i)-(*(zh+i*n+j))*(*(df+j));
   }
   k=k+1;
   goto L2;
}

风花雪月 发表于 2007-5-10 01:29

funct.c
#include "stdio.h"
#include "stdlib.h"
#include "math.h"
const int kkg=3;
double r0;
double f(double x[])
{double ff;
ff=pow((x-8),2)+pow((x)-8,2);
return(ff);
}
void strain(doublex[],double g[])
{   g=x-1;
    g=x-1;
    g=11-x-x;
}
double objf(double p[])
{int i;
double ff,sg,*g;
g=(double *)malloc(kkg *sizeof(double));
sg=0;
strain(p,g);
for(i=0;i<kkg;i++)
{if(*(g+i)>0) sg=sg+r0/(*(g+i));
   else sg=sg+r0*(1e+10);
}
free(g);
ff=f(p)+sg;
return(ff);
}

风花雪月 发表于 2007-5-10 01:30

hjfgf.c
#include "jtf.c"
double gold(double a[],double b[],double eps,int n,double xx[])
{int i;
double f1,f2,*x,ff,q,w;
for(i=0;i<2;i++)
x=(double *)malloc(n*sizeof(double));
for(i=0;i<n;i++)
{*(x+i)=a+0.618*(b-a);
   *(x+i)=a+0.382*(b-a);
}
f1=objf(x);
f2=objf(x);
do
   {if(f1>f2)
   {for(i=0;i<n;i++)
      {b=*(x+i);
       *(x+i)=*(x+i);
       }
   f1=f2;
   for(i=0;i<n;i++)
      *(x+i)=a+0.382*(b-a);
   f2=objf(x);
   }
    else
   { for(i=0;i<n;i++)
       {a=*(x+i);
       *(x+i)=*(x+i);}
   f2=f1;
    for(i=0;i<n;i++)
   *(x+i)=a+0.618*(b-a);
    f1=objf(x);
   }
q=0;
for(i=0;i<n;i++)
   q=q+(b-a)*(b-a);
w=sqrt(q);
}while(w>eps);
for(i=0;i<n;i++)
   xx=0.5*(a+b);
ff=objf(xx);
for(i=0;i<2;i++)
free(x);
return(ff);
}

风花雪月 发表于 2007-5-10 01:30

jtf.c
#include "funct.c"
void jtf(double x0[],double h0,double s[],int n,double a[],double b[])
{int i;
double *x,h,f1,f2,f3;
for(i=0;i<3;i++)
x=(double *)malloc(n*sizeof(double));
h=h0;
for(i=0;i<n;i++)
*(x+i)=x0;
f1=objf(x);
for(i=0;i<n;i++)
*(x+i)=*(x+i)+h*s;
f2=objf(x);
if(f2>=f1)
{h=-h0;
    for(i=0;i<n;i++)
    *(x+i)=*(x+i);
   f3=f1;
    for(i=0;i<n;i++)
    {*(x+i)=*(x+i);
   *(x+i)=*(x+i);
    }
   f1=f2;
   f2=f3;
   }
   for(;;)
   {h=2*h;
   for(i=0;i<n;i++)
   *(x+i)=*(x+i)+h*s;
   f3=objf(x);
   if(f2<f3) break;
   else
    { for(i=0;i<n;i++)
       {*(x+i)=*(x+i);
      *(x+i)=*(x+i);
       }
      f1=f2;
      f2=f3;
    }
   }
   if(h<0)
    for(i=0;i<n;i++)
    {a=*(x+i);
   b=*(x+i);
    }
   else
    for(i=0;i<n;i++)
    {a=*(x+i);
   b=*(x+i);
   }
   for(i=0;i<3;i++)
   free(x);
}

风花雪月 发表于 2007-5-10 01:31

mldfhs1.c
#include "dfpopt.c"
main()
{int i;
double p[]={3,4};
double fom,fxo,x;
double c=0.1;
doubler0=120;
fom=100;
do
{ fxo=dfpopt(p,0.2,0.0001,0.000001,2);
    if(fabs(fom-fxo)>0.001)
    {fom=fxo;
   r0=c*r0;
   for(i=0;i<2;i++)
   *(p+i)=x;
    }
   else
    {printf("\nx=%f,x=%f,ff=%f",x,x,fxo);
   return; }
    }while(1);
    getch();
}

风花雪月 发表于 2007-5-10 01:31

powell.c
#include "hjfgf.c"
double oneoptim(double x0[],double s[],double h0,double epsg,int n,double x[])
{double *a,*b,ff;
a=(double *)malloc(n*sizeof(double));
b=(double *)malloc(n*sizeof(double));
jtf(x0,h0,s,n,a,b);
ff=gold(a,b,epsg,n,x);
free(a);
free(b);
return (ff);
}
double powell(double p[],double h0,double eps,double epsg,int n,double x[])
{int i,j,m;
double *xx,*ss,*s;
double f,f0,f1,f2,f3,fx,dlt,df,sdx,q,d;
ss=(double *)malloc(n*(n+1)*sizeof(double));
s=(double *)malloc(n*sizeof(double));
for(i=0;i<n;i++)
{for(j=0;j<=n;j++)
   *(ss+i*(n+1)+j)=0;
   *(ss+i*(n+1)+i)=1;
}
for(i=0;i<4;i++)
xx=(double *)malloc(n*sizeof(double));
for(i=0;i<n;i++)
*(xx+i)=p;
for(;;)
{for(i=0;i<n;i++)
    {*(xx+i)=*(xx+i);
   x=*(xx+i);
    }
   f0=f1=objf(x);
   dlt=-1;
   for(j=0;j<n;j++)
    {for(i=0;i<n;i++)
      {*(xx+i)=x;
       *(s+i)=*(ss+i*(n+1)+j);
      }
   f=oneoptim(xx,s,h0,epsg,n,x);
   df=f0-f;
   if(df>dlt)
      {dlt=df;
       m=j;
      }
    }
    sdx=0;
    for(i=0;i<n;i++)
   sdx=sdx+fabs(x-(*(xx+i)));
    if(sdx<eps)
   {free(ss);
      free(s);
      for(i=0;i<4;i++)
      free(xx);
      return(f);
   }
    for(i=0;i<n;i++)
   *(xx+i)=x;
    f2=f;
    for(i=0;i<n;i++)
    {*(xx+i)=2*(*(xx+i)-(*(xx+i)));
   x=*(xx+i);
    }
    fx=objf(x);
    f3=fx;
    q=(f1-2*f2+f3)*(f1-f2-dlt)*(f1-f2-dlt);
    d=0.5*dlt*(f1-f3)*(f1-f3);
    if((f3<f1)||(q<d))
   {if(f2<=f3)
       for(i=0;i<n;i++)
      *(xx+i)=*(xx+i);
      else
       for(i=0;i<n;i++)
      *(xx+i)=*(xx+i);
   }
    else
    {for(i=0;i<n;i++)
   {*(ss+(i+1)*(n+1))=x-(*(xx+i));
      *(s+i)=*(ss+(i+1)*(n+1));
   }
   f=oneoptim(xx,s,h0,epsg,n,x);
   for(i=0;i<n;i++)
    *(xx+i)=x;
   for(j=m+1;j<=n;j++)
   for(i=0;i<n;i++)
    *(ss+i*(n+1)+j-1)=*(ss+i*(n+1)+j);
    }
}
}

donglei06 发表于 2007-5-16 21:32

疑问

照这些,怎么运行出不来结果啊

风花雪月 发表于 2007-5-17 15:16

程序应该是没有问题的，但是是否能用于你的具体问题不好说

由于我现在没有装C编译器，无法帮你调试

donglei06 发表于 2007-5-17 22:40

请教你个问题！

我的 是个有约束条件 的 多元函数，应该直接套用吧，可老提示溢出，我已经用了double了 ，真是 没有办法了 ，希望那位高手能给我 发个能出结果的 程序，小弟将不胜感激！我的 邮箱是donglei06@163.com

donglei06 发表于 2007-5-18 16:12

请教!

那几个约束条件怎么没有用上啊,结果不受约束啊,所求值是极小值,参数值也不杂约束条件内啊!

donglei06 发表于 2007-5-21 08:00

疑问

照你的 程序,运行结果上,参变量不受约束条件限制啊!

donglei06 发表于 2007-5-21 08:01

疑问

目标函数 objf可能有问题吧

xwc6134 发表于 2007-5-21 17:11

是啊，运行结果不受条件约束，怎么办

博闻 发表于 2007-6-29 00:25

谢谢楼主了
我会去\试试看的

查看完整版本: 内点惩罚函数法子程序