内点惩罚函数法子程序

风花雪月

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

本程序由heizi友情提供

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

风花雪月

dfpopt.c

1. #include "hjfgf.c"
   
2. #include "stdlib.h"
   
3. void gradient(double x[],double g[],int n)
   
4. {int i;
   
5.   double af,f1,f2,dltx=0.000001;
   
6.   for(i=0;i<n;i++)
   
7.   g[i]=0;
   
8.   f1=objf(x);
   
9.   for(i=0;i<n;i++)
   
10.   {af=*(x+i);
    
11.    *(x+i)=af+dltx;
    
12.    f2=objf(x);
    
13.    g[i]=(f2-f1)/dltx;
    
14.    *(x+i)=af;
    
15.   }
    
16. }
    
17. double oneoptim(double x0[],double s[],double h0,double epsg,int n,double x[])
    
18. {double *a,*b,ff;
    
19.   a=(double *)malloc(n*sizeof (double ));
    
20.   b=(double *)malloc(n*sizeof (double ));
    
21.   jtf(x0,h0,s,n,a,b);
    
22.   ff=gold(a,b,epsg,n,x);
    
23.   free(a);
    
24.   free(b);
    
25.   return(ff);
    
26. }
    
27. double dfpopt(double xx[],double h0,double eps,double epsg,int n)
    
28. {int i,j,k;
    
29.   double ae,zcc;
    
30.   double *s,*x,*ay[2],*df[2],*zd[2],*zc[2],*zh[2];
    
31.   s=(double *)malloc(n*sizeof (double ));
    
32.   x=(double *)malloc(n*sizeof (double ));
    
33.   for(i=0;i<2;i++)
    
34.   {ay[i]=(double*)malloc(n*sizeof(double));
    
35.    df[i]=(double*)malloc(n*sizeof(double));
    
36.    zd[i]=(double*)malloc(n*sizeof(double));
    
37.    zc[i]=(double*)malloc(n*sizeof(double));
    
38.    zh[i]=(double*)malloc(n*n*sizeof(double));
    
39.   }
    
40.   for(i=0;i<n;i++)
    
41.   *(ay[0]+i)=xx[i];
    
42.   L1:
    
43.    k=0;
    
44.    for (i=0;i<n;i++)
    
45.    for (j=0;j<n;j++)
    
46.     {*(zh[0]+i*n+j)=0;
    
47.      if(i==j)
    
48.       *(zh[0]+i*n+j)=1;
    
49.     }
    
50.   for(i=0;i<n;i++)
    
51.   {*(x+i)=*(ay[0]+i);
    
52.    *(ay[1]+i)=*(ay[0]+i);
    
53.   }
    
54.   gradient(x,df[0],n);
    
55.   for(i=0;i<n;i++)
    
56.   {*(s+i)=0;
    
57.    *(df[1]+i)=*(df[0]+i);
    
58.    for(j=0;j<n;j++)
    
59.    *(s+i)=*(s+i)-(*(zh[0]+i*n+j))*(*(df[0]+j));
    
60.   }
    
61.   L2:
    
62.   oneoptim(x,s,h0,epsg,n,ay[0]);
    
63.   for(i=0;i<n;i++)
    
64.    *(x+i)=*(ay[0]+i);
    
65.   gradient(x,df[0],n);
    
66.   ae=0;
    
67.   for(i=0;i<n;i++)
    
68.    ae=ae+(*(df[0]+i))*(*(df[0]+i));
    
69.   if(ae<=eps)
    
70.    {for(i=0;i<n;i++)
    
71.     *(xx+i)=*(x+i);
    
72.     free(s);
    
73.     free(x);
    
74.     for(i=0;i<2;i++)
    
75.     {free(ay[i]);
    
76.      free(df[i]);
    
77.      free(zd[i]);
    
78.      free(zc[i]);
    
79.      free(zh[i]);
    
80.     }
    
81.     return(objf(xx));
    
82.    }
    
83.    if(k==n) goto L1;
    
84.    zcc=0;
    
85.    for(i=0;i<n;i++)
    
86.    {*(zd[0]+i)=*(ay[0]+i)-(*(ay[1]+i));
    
87.     *(zd[1]+i)=*(df[0]+i)-(*(df[1]+i));
    
88.     *(df[1]+i)=*(df[0]+i);
    
89.     zcc=zcc+(*(zd[0]+i))*(*(zd[1]+i));
    
90.    }
    
91.    for(i=0;i<n;i++)
    
92.     for(j=0;j<n;j++)
    
93.     *(zh[1]+i*n+j)=*(zd[0]+i)*(*(zd[0]+j))/zcc;
    
94.    for(i=0;i<n;i++)
    
95.     { *(zc[0]+i)=0;
    
96.       for(j=0;j<n;j++)
    
97.        *(zc[0]+i)=*(zc[0]+i)+(*(zd[1]+j))*(*(zh[0]+j*n+i));
    
98.     }
    
99.      zcc=0;
    
100.     for(i=0;i<n;i++)
     
101.      zcc=zcc+(*(zc[0]+i))*(*(zd[1]+i));
     
102.      for(i=0;i<n;i++)
     
103.      {*(zc[0]+i)=0;
     
104.       *(zc[1]+i)=0;
     
105.       for (j=0;j<n;j++)
     
106.        {*(zc[0]+i)=*(zc[0]+i)+(*(zh[0]+i*n+j))*(*(zd[1]+j));
     
107.         *(zc[1]+i)=*(zc[1]+i)+(*(zh[1]+j))*(*(zh[0]+j*n+i));
     
108.        }
     
109.      }
     
110.     for(i=0;i<n;i++)
     
111.     for(j=0;j<n;j++)
     
112.     *(zh[0]+i*n+j)=*(zh[0]+i*n+j)+(*(zh[1]+i*n+j))-(*(zc[0]+i))*(*(zc[1]+j))/zcc;
     
113.     for(i=0;i<n;i++)
     
114.      {*(s+i)=0;
     
115.       for(j=0;j<n;j++)
     
116.        *(s+i)=*(s+i)-(*(zh[0]+i*n+j))*(*(df[0]+j));
     
117.      }
     
118.      k=k+1;
     
119.      goto L2;
     
120. }

复制代码

风花雪月

funct.c

1. #include "stdio.h"
   
2. #include "stdlib.h"
   
3. #include "math.h"
   
4. const int kkg=3;
   
5. double r0;
   
6. double f(double x[])
   
7. {double ff;
   
8. ff=pow((x[0]-8),2)+pow((x[1])-8,2);
   
9. return(ff);
   
10. }
    
11. void strain(double  x[],double g[])
    
12. {   g[0]=x[0]-1;
    
13.     g[1]=x[1]-1;
    
14.     g[2]=11-x[0]-x[1];
    
15. }
    
16. double objf(double p[])
    
17. {int i;
    
18. double ff,sg,*g;
    
19. g=(double *)malloc(kkg *sizeof(double));
    
20. sg=0;
    
21. strain(p,g);
    
22. for(i=0;i<kkg;i++)
    
23.   {if(*(g+i)>0) sg=sg+r0/(*(g+i));
    
24.    else sg=sg+r0*(1e+10);
    
25.   }
    
26.   free(g);
    
27.   ff=f(p)+sg;
    
28.   return(ff);
    
29. }

复制代码

风花雪月

hjfgf.c

1. #include "jtf.c"
   
2. double gold(double a[],double b[],double eps,int n,double xx[])
   
3. {int i;
   
4. double f1,f2,*x[2],ff,q,w;
   
5. for(i=0;i<2;i++)
   
6.   x[i]=(double *)malloc(n*sizeof(double));
   
7. for(i=0;i<n;i++)
   
8.   {*(x[0]+i)=a[i]+0.618*(b[i]-a[i]);
   
9.    *(x[1]+i)=a[i]+0.382*(b[i]-a[i]);
   
10.   }
    
11.   f1=objf(x[0]);
    
12.   f2=objf(x[1]);
    
13.   do
    
14.    {if(f1>f2)
    
15.      {for(i=0;i<n;i++)
    
16.       {b[i]=*(x[0]+i);
    
17.        *(x[0]+i)=*(x[1]+i);
    
18.        }
    
19.      f1=f2;
    
20.      for(i=0;i<n;i++)
    
21.       *(x[1]+i)=a[i]+0.382*(b[i]-a[i]);
    
22.      f2=objf(x[1]);
    
23.      }
    
24.     else
    
25.      { for(i=0;i<n;i++)
    
26.        {a[i]=*(x[1]+i);
    
27.        *(x[1]+i)=*(x[0]+i);}
    
28.      f2=f1;
    
29.     for(i=0;i<n;i++)
    
30.      *(x[0]+i)=a[i]+0.618*(b[i]-a[i]);
    
31.     f1=objf(x[0]);
    
32.      }
    
33.   q=0;
    
34.   for(i=0;i<n;i++)
    
35.    q=q+(b[i]-a[i])*(b[i]-a[i]);
    
36.   w=sqrt(q);
    
37.   }while(w>eps);
    
38.   for(i=0;i<n;i++)
    
39.    xx[i]=0.5*(a[i]+b[i]);
    
40.   ff=objf(xx);
    
41.   for(i=0;i<2;i++)
    
42.   free(x[i]);
    
43.   return(ff);
    
44. }

复制代码

风花雪月

jtf.c

1. #include "funct.c"
   
2. void jtf(double x0[],double h0,double s[],int n,double a[],double b[])
   
3. {int i;
   
4. double *x[3],h,f1,f2,f3;
   
5. for(i=0;i<3;i++)
   
6.   x[i]=(double *)malloc(n*sizeof(double));
   
7.   h=h0;
   
8. for(i=0;i<n;i++)
   
9.   *(x[0]+i)=x0[i];
   
10. f1=objf(x[0]);
    
11. for(i=0;i<n;i++)
    
12.   *(x[1]+i)=*(x[0]+i)+h*s[i];
    
13.   f2=objf(x[1]);
    
14. if(f2>=f1)
    
15.   {h=-h0;
    
16.     for(i=0;i<n;i++)
    
17.     *(x[2]+i)=*(x[0]+i);
    
18.    f3=f1;
    
19.     for(i=0;i<n;i++)
    
20.     {*(x[0]+i)=*(x[1]+i);
    
21.      *(x[1]+i)=*(x[2]+i);
    
22.     }
    
23.    f1=f2;
    
24.    f2=f3;
    
25.    }
    
26.    for(;;)
    
27.    {h=2*h;
    
28.      for(i=0;i<n;i++)
    
29.      *(x[2]+i)=*(x[1]+i)+h*s[i];
    
30.    f3=objf(x[2]);
    
31.    if(f2<f3) break;
    
32.    else
    
33.     { for(i=0;i<n;i++)
    
34.        {*(x[0]+i)=*(x[1]+i);
    
35.         *(x[1]+i)=*(x[2]+i);
    
36.        }
    
37.       f1=f2;
    
38.       f2=f3;
    
39.     }
    
40.    }
    
41.    if(h<0)
    
42.     for(i=0;i<n;i++)
    
43.     {a[i]=*(x[2]+i);
    
44.      b[i]=*(x[0]+i);
    
45.     }
    
46.    else
    
47.     for(i=0;i<n;i++)
    
48.     {a[i]=*(x[0]+i);
    
49.      b[i]=*(x[2]+i);
    
50.      }
    
51.    for(i=0;i<3;i++)
    
52.    free(x[i]);
    
53. }

复制代码

风花雪月

mldfhs1.c

1. #include "dfpopt.c"
   
2. main()
   
3. {int i;
   
4. double p[]={3,4};
   
5. double fom,fxo,x[2];
   
6. double c=0.1;
   
7. double  r0=120;
   
8. fom=100;
   
9. do
   
10.   { fxo=dfpopt(p,0.2,0.0001,0.000001,2);
    
11.     if(fabs(fom-fxo)>0.001)
    
12.     {fom=fxo;
    
13.      r0=c*r0;
    
14.      for(i=0;i<2;i++)
    
15.      *(p+i)=x[i];
    
16.     }
    
17.    else
    
18.     {printf("\nx[0]=%f,x[1]=%f,ff=%f",x[0],x[1],fxo);
    
19.      return; }
    
20.     }while(1);
    
21.     getch();
    
22. }

复制代码

风花雪月

powell.c

1. #include "hjfgf.c"
   
2. double oneoptim(double x0[],double s[],double h0,double epsg,int n,double x[])
   
3. {double *a,*b,ff;
   
4. a=(double *)malloc(n*sizeof(double));
   
5. b=(double *)malloc(n*sizeof(double));
   
6. jtf(x0,h0,s,n,a,b);
   
7. ff=gold(a,b,epsg,n,x);
   
8. free(a);
   
9. free(b);
   
10. return (ff);
    
11. }
    
12. double powell(double p[],double h0,double eps,double epsg,int n,double x[])
    
13. {int i,j,m;
    
14. double *xx[4],*ss,*s;
    
15. double f,f0,f1,f2,f3,fx,dlt,df,sdx,q,d;
    
16. ss=(double *)malloc(n*(n+1)*sizeof(double));
    
17. s=(double *)malloc(n*sizeof(double));
    
18. for(i=0;i<n;i++)
    
19.   {for(j=0;j<=n;j++)
    
20.    *(ss+i*(n+1)+j)=0;
    
21.    *(ss+i*(n+1)+i)=1;
    
22.   }
    
23.   for(i=0;i<4;i++)
    
24.   xx[i]=(double *)malloc(n*sizeof(double));
    
25.   for(i=0;i<n;i++)
    
26.   *(xx[0]+i)=p[i];
    
27.   for(;;)
    
28.   {for(i=0;i<n;i++)
    
29.     {*(xx[1]+i)=*(xx[0]+i);
    
30.      x[i]=*(xx[1]+i);
    
31.     }
    
32.    f0=f1=objf(x);
    
33.    dlt=-1;
    
34.    for(j=0;j<n;j++)
    
35.     {for(i=0;i<n;i++)
    
36.       {*(xx[0]+i)=x[i];
    
37.        *(s+i)=*(ss+i*(n+1)+j);
    
38.       }
    
39.      f=oneoptim(xx[0],s,h0,epsg,n,x);
    
40.      df=f0-f;
    
41.      if(df>dlt)
    
42.       {dlt=df;
    
43.        m=j;
    
44.       }
    
45.     }
    
46.     sdx=0;
    
47.     for(i=0;i<n;i++)
    
48.      sdx=sdx+fabs(x[i]-(*(xx[1]+i)));
    
49.     if(sdx<eps)
    
50.      {free(ss);
    
51.       free(s);
    
52.       for(i=0;i<4;i++)
    
53.       free(xx[i]);
    
54.       return(f);
    
55.      }
    
56.     for(i=0;i<n;i++)
    
57.      *(xx[2]+i)=x[i];
    
58.     f2=f;
    
59.     for(i=0;i<n;i++)
    
60.     {*(xx[3]+i)=2*(*(xx[2]+i)-(*(xx[1]+i)));
    
61.      x[i]=*(xx[3]+i);
    
62.     }
    
63.     fx=objf(x);
    
64.     f3=fx;
    
65.     q=(f1-2*f2+f3)*(f1-f2-dlt)*(f1-f2-dlt);
    
66.     d=0.5*dlt*(f1-f3)*(f1-f3);
    
67.     if((f3<f1)||(q<d))
    
68.      {if(f2<=f3)
    
69.        for(i=0;i<n;i++)
    
70.         *(xx[0]+i)=*(xx[2]+i);
    
71.       else
    
72.        for(i=0;i<n;i++)
    
73.         *(xx[0]+i)=*(xx[3]+i);
    
74.      }
    
75.     else
    
76.     {for(i=0;i<n;i++)
    
77.      {*(ss+(i+1)*(n+1))=x[i]-(*(xx[1]+i));
    
78.       *(s+i)=*(ss+(i+1)*(n+1));
    
79.      }
    
80.    f=oneoptim(xx[0],s,h0,epsg,n,x);
    
81.    for(i=0;i<n;i++)
    
82.     *(xx[0]+i)=x[i];
    
83.    for(j=m+1;j<=n;j++)
    
84.    for(i=0;i<n;i++)
    
85.     *(ss+i*(n+1)+j-1)=*(ss+i*(n+1)+j);
    
86.     }
    
87.   }
    
88. }

复制代码

donglei06

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

风花雪月

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

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

donglei06

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

donglei06

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

donglei06

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

donglei06

目标函数 objf可能有问题吧

xwc6134

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

博闻

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