外点惩罚函数法优化子程序

风花雪月

本程序包含5个C文件
mwdfhs.c
powell.c
funct.c
jtf.c
hjfgf.c


本程序由heizi友情提供

风花雪月

mwdfhs.c代码如下：

1. #include "powell.c"
   
2. void main()
   
3. {int i;
   
4. double p[]={7,8};
   
5. double fom,fxo,c,x[2];
   
6. c=10;
   
7. r0=0.1;
   
8. fom=100;
   
9. do
   
10. {fxo=powell(p,0.1,0.0001,0.0001,2,x);
    
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("输出最优点及其目标函数值：\n");
    
19.     printf("x[0]=%f,x[1]=%f,ff=%f",x[0],x[1],fxo);
    
20.     return;
    
21.    }
    
22. }while(1);
    
23. }

复制代码

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. }

复制代码

风花雪月

funct.c代码如下：

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

复制代码

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. }

复制代码

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. }

复制代码

mas3593

请教一下楼主，为什么我把这几个程序连起来了调试过了，但运行时就会弹出对话框，忽略以后输出不正确

风花雪月

原帖由 mas3593 于 2006-12-26 10:34 发表
请教一下楼主，为什么我把这几个程序连起来了调试过了，但运行时就会弹出对话框，忽略以后输出不正确

说明详细的错误信息

dxawicso

请教楼主,同样是编译通过,但运行出错.
可否请楼主自己编译运行一下,了解具体情况

heijian20

在Turbo C环境下可以调试通过的

qxkice2004

lz能不能解释一下这几个子程序的含义和用法啊，感激！

wanghuaipan

我的问题同上

程瑞岩

有些东西现在都看不懂了。

nari

编译通过,但运行出错.
可否请楼主自己编译运行一下,了解具体情况