内点法优化程序

donglei06

内点法优化程序

1. #include<math.h>
   
2. #include<stdlib.h>
   
3. #include<stdio.h>
   
4. const int kkg=8;
   
5. double r0;
   
6. void main()
   
7. {double f(double x[]);
   
8. double objf(double p[]);
   
9. void jtf(double x0[],double h0, double s[],int n, double a[],double b[]);
   
10. double gold(double a[],double b[],double eps,int n, double xx[]);
    
11. double oneoptim(double x0[],double s[],double h0, double epsg,int n, double x[]);
    
12. double powell(double p[],double h0, double eps, double epsg,int n, double x[]);
    
13. int i;
    
14. double p[]={21,35,108,16};
    
15. double fom,fxo,c,x[4];
    
16. c=0.5;
    
17. r0=45800;
    
18. fom=10000;
    
19. do
    
20. {
    
21. fxo=powell(p,0.1,1500,650,4,x);
    
22. if(fabs(fom-fxo)>1800)
    
23. {fom=fxo;
    
24. r0=c*r0;
    
25. for(i=0;i<4;i++)
    
26. *(p+i)=x[i];
    
27. }
    
28. else
    
29. {printf("x[0]=%f,x[1]=%f,x[2]=%f,x[3]=%f\n",x[0],x[1],x[2],x[3]);
    
30. return;
    
31. }
    
32. }while(1);
    
33. }
    
34. double f(double x[])
    
35. { double ff;
    
36. ff=0.8322*x[2]*(3.7889*x[0]*x[0]+5.8174*x[0]+3.9422)*pow((0.9715*x[3]*x[0]-0.5*x[2]),2)
    
37. /(cos(0.0087*x[1])*x[0]*x[0]);
    
38. return(ff);
    
39. }
    
40. void strain(double x[], double g[])
    
41. {g[0]=4.842e-6*x[0]*x[2]*x[3]*x[3]-1.0;
    
42. g[1]=4.86e-4*x[0]*x[3]*sqrt(x[2])-1.0;
    
43. g[2]=x[0]-12;
    
44. g[3]=x[1]-30;
    
45. g[4]=45-x[1];
    
46. g[5]=x[3]-1.5;
    
47. g[6]=x[2]-0.2429*x[3]*x[0];
    
48. g[7]=0.3238*x[3]*x[0]-x[2];
    
49. }
    
50. double objf(double p[])
    
51. {int i;
    
52. double ff,sg,*g;
    
53. g=(double*)malloc(kkg*sizeof(double));
    
54. sg=0;
    
55. for(i=0;i<kkg;i++)
    
56. {if(*(g+i)>0)
    
57. sg=sg+r0/(*(g+i));
    
58. else
    
59. sg=sg+r0*(1e+10);
    
60. }
    
61. free(g);
    
62. ff=f(p)+sg;
    
63. return(ff);
    
64. }
    
65. void jtf(double x0[],double h0, double s[],int n, double a[],double b[])
    
66. {int i;
    
67. double *x[3],h,f1,f2,f3;
    
68. for(i=0;i<3;i++)
    
69. x[i]=(double*)malloc(n*sizeof(double));
    
70. h=h0;
    
71.   for(i=0;i<n;i++)
    
72.    *(x[0]+i)=x0[i];
    
73.    f1=objf(x[0]);
    
74.    for(i=0;i<n;i++)
    
75.    *(x[1]+i)=*(x[0]+i)+h*s[i];
    
76.    f2=objf(x[1]);
    
77. if(f2>=f1)
    
78.    {h=-h0;
    
79.      for(i=0;i<n;i++)
    
80.      *(x[2]+i)=*(x[0]+i);
    
81.      f3=f1;
    
82.      for(i=0;i<n;i++)
    
83.       {*(x[0]+i)=*(x[1]+i);
    
84.       *(x[1]+i)=*(x[2]+i);
    
85.       }
    
86.      f1=f2;
    
87.      f2=f3;
    
88.    }
    
89. for(;;)
    
90.   {h=2*h;
    
91.    for(i=0;i<n;i++)
    
92.    *(x[2]+i)=*(x[1]+i)+h*s[i];
    
93.    f3=objf(x[2]);
    
94.    if(f2<f3)
    
95.      break;
    
96. else
    
97. {for(i=0;i<n;i++)
    
98.   {*(x[0]+i)=*(x[1]+i);
    
99.    *(x[1]+i)=*(x[2]+i);
    
100.   }
     
101. f1=f2;
     
102. f2=f3;
     
103. }
     
104. }
     
105. if(h<0)
     
106. for(i=0;i<n;i++)
     
107.   {a[i]=*(x[2]+i);
     
108.    b[i]=*(x[0]+i);
     
109.   }
     
110. else
     
111. for(i=0;i<n;i++)
     
112. {a[i]=*(x[0]+i);
     
113.   b[i]=*(x[2]+i);
     
114. }
     
115. for(i=0;i<3;i++)
     
116.   free(x[i]);
     
117. }
     
118. double gold(double a[],double b[],double eps,int n, double xx[])
     
119. {int i;
     
120. double f1,f2,*x[2],ff,q,w;
     
121.   for(i=0;i<2;i++)
     
122.   x[i]=(double*)malloc(n*sizeof(double));
     
123. for(i=0;i<n;i++)
     
124.   {*x[0]=a[i]+0.618*(b[i]-a[i]);
     
125.    *(x[1]+i)=a[i]+0.382*(b[i]-a[i]);
     
126.   }
     
127. f1=objf(x[0]);
     
128. f2=objf(x[1]);
     
129. do
     
130. {if(f1>f2)
     
131.    { for(i=0;i<n;i++)
     
132.      {b[i]=*(x[0]+i);
     
133.       *(x[0]+i)=*(x[1]+i);
     
134.   }
     
135. f1=f2;
     
136. for(i=0;i<n;i++)
     
137.    *(x[1]+i)=a[i]+0.382*(b[i]-a[i]);
     
138.    f2=objf(x[1]);
     
139. }
     
140. else
     
141. { for(i=0;i<n;i++)
     
142.    {a[i]=*(x[1]+i);
     
143.      *(x[1]+i)= *(x[0]+i);
     
144.     }
     
145.   f2=f1;
     
146. for(i=0;i<n;i++)
     
147. *x[0]=a[i]+0.618*(b[i]-a[i]);
     
148. f1=objf(x[0]);
     
149. }
     
150. q=0;
     
151. for(i=0;i<n;i++)
     
152.   q=q+(b[i]-a[i])*(b[i]-a[i]);
     
153.   w=sqrt(q);
     
154. }while(w>eps);
     
155.   for(i=0;i<n;i++)
     
156.   xx[i]=0.5*(a[i]+b[i]);
     
157.   ff=objf(xx);
     
158.   for(i=0;i<2;i++)
     
159.   free(x[i]);
     
160.   return(ff);
     
161. }
     
162. double oneoptim(double x0[],double s[],double h0,double epsg,int n,double x[])
     
163. {double *a,*b,ff;
     
164. a=(double*)malloc(n*sizeof(double));
     
165. b=(double*)malloc(n*sizeof(double));
     
166. jtf(x0,h0,s,n,a,b);
     
167. ff=gold(a,b,epsg,n,x);
     
168. free(a);
     
169. free(b);
     
170. return(ff);
     
171. }
     
172. double powell(double p[],double h0,double eps,double epsg,int n,double x[])
     
173. {int i,j,m;
     
174. double *xx[4],*ss,*s;
     
175. double f,f0,f1,f2,f3,fx,dlt,df,sdx,q,d;
     
176. ss=(double*)malloc(n*(n+1)*sizeof(double));
     
177. s=(double*)malloc(n*sizeof(double));
     
178. for(i=0;i<n;i++)
     
179. {for(j=0;j<=n;j++)
     
180.   *(ss+i*(n+1)+j)=0;
     
181.   *(ss+i*(n+1)+i)=1;
     
182. }
     
183. for(i=0;i<4;i++)
     
184.   xx[i]=(double*)malloc(n*sizeof(double));
     
185.   for(i=0;i<n;i++)
     
186.   *(xx[0]+i)=p[i];
     
187. for(;;)
     
188. {for(i=0;i<n;i++)
     
189. {*(xx[1]+i)=*(xx[0]+i);
     
190.   x[i]=*(xx[1]+i);
     
191.   }
     
192. f0=f1=objf(x);
     
193. dlt=-1;
     
194. for(j=0;j<n;j++)
     
195.   {for(i=0;i<n;i++)
     
196.    {*(xx[0]+i)=x[i];
     
197.     *(s+i)=*(ss+i*(n+1)+j);
     
198.    }
     
199. f=oneoptim(xx[0],s,h0,epsg,n,x);
     
200. df=f0-f;
     
201. if(df>dlt)
     
202. {dlt=df;
     
203.   m=j;
     
204. }
     
205. }
     
206. sdx=0;
     
207. for(i=0;i<n;i++)
     
208. sdx=sdx+fabs(x[i]-(*(xx[1]+i)));
     
209. if(sdx<eps)
     
210.   {free(ss);
     
211.    free(s);
     
212.    for(i=0;i<4;i++)
     
213.     free(xx[i]);
     
214.    return(f);
     
215. }
     
216. for(i=0;i<n;i++)
     
217.   *(xx[2]+i)=x[i];
     
218. f2=f;
     
219. for(i=0;i<n;i++)
     
220.   {*(xx[3]+i)=2*(*(xx[2]+i)-(*(xx[1]+i)));
     
221.    x[i]=*(xx[3]+i);
     
222.   }
     
223. fx=objf(x);
     
224. f3=fx;
     
225. q=(f1-2*f2+f3)*(f1-f2-dlt)*(f1-f2-dlt);
     
226. d=0.5*dlt*(f1-f3)*(f1-f3);
     
227. if((f3<f1)||(q<d))
     
228.    {if(f2<=f3)
     
229.       for(i=0;i<n;i++)
     
230.        *(xx[0]+i)=*(xx[2]+i);
     
231.       else
     
232.         for(i=0;i<n;i++)
     
233.          *(xx[0]+i)=*(xx[3]+i);
     
234.    }
     
235. else
     
236.   {for(i=0;i<n;i++)
     
237.    {*(ss+(i+1)*(n+1))=x[i]-(*(xx[1]+i));
     
238.     *(s+i)=*(ss+(i+1)*(n+1));
     
239.    }
     
240.   f=oneoptim(xx[0],s,h0,epsg,n,x);
     
241.   for(i=0;i<n;i++)
     
242.    *(xx[0]+i)=x[i];
     
243.   for(j=m+1;j<=n;j++)
     
244.   for(i=0;i<n;i++)
     
245.    *(ss+i*(n+1)+j-i)=*(ss+i*(n+1)+j);
     
246. }
     
247. }
     
248. }

复制代码

运行结果是:Floating point error:Over flow

[ 本帖最后由 风花雪月 于 2007-6-15 11:15 编辑 ]

hnsqzrj

能运行么

总有一天

对不起，我觉得还挺长的，建议楼主以后加上注释

风花雪月

自己好好找找是什么地方溢出，我运行了一下，好像需要运行很长时间

估计这样的问题，别人是没有那么多时间帮你调试的