哪位能告诉我正确的复合形算法啊?谢谢!

sun1993

我倒. 手头上三本关于优化的书, 介绍的复合形法都有点不同,
试着编了三个程序验证一下哪个比较好, 优化出的结果都离正确值偏差很多.
搞不懂怎么会这样!
请大家帮帮我!谢谢!

风花雪月

复合形法怎么会不一样的？不能理解

1. /* 以下为复合形法法子程序 */
   
2. 
   
3. #include "stdlib.h"
   
4. 
   
5. #include "math.h"
   
6. 
   
7. #include "stdio.h"
   
8. 
   
9. float objfx(float x[]);
   
10. 
    
11. void constraint(float x[],float g[]);
    
12. int gau(float x[],float g[],int kg)
    
13. 
    
14. {
    
15. 
    
16.        int i;
    
17. 
    
18.        constraint(x,g);
    
19. 
    
20.        for(i=0;i<kg;i++)
    
21. 
    
22.        {
    
23. 
    
24.          if(g[i]<0)
    
25. 
    
26.          goto s333;
    
27. 
    
28.        }
    
29. 
    
30.        return 1;
    
31. 
    
32. s333:return 0;
    
33. 
    
34. }
    
35. void xcent(int n,int k,int ll,int lh,float x0[],float xcom[][100])
    
36. 
    
37. {
    
38. 
    
39.        int i,l;
    
40. 
    
41.        float xs;
    
42. 
    
43.        for(i=0;i<n;i++)
    
44. 
    
45.        {
    
46. 
    
47.          xs=0;
    
48. 
    
49.          for(l=0;l<ll;l++)
    
50. 
    
51.          {
    
52. 
    
53.     if(l!=lh)
    
54. 
    
55.     xs=xs+xcom[i][l];
    
56. 
    
57.          }
    
58. 
    
59.          if(lh>-1)
    
60. 
    
61.     x0[i]=xs/(ll-1);
    
62. 
    
63.          else
    
64. 
    
65.     x0[i]=xs/ll;
    
66. 
    
67.        }
    
68. 
    
69. }
    
70. void fxse(int n,int k,float x[],float xcom[][100],float fxk[])
    
71. 
    
72. {
    
73. 
    
74.        int l,lp,lp1,i;
    
75. 
    
76.        float w;
    
77. 
    
78.        for(l=0;l<k-1;l++)
    
79. 
    
80.        for(lp=0;lp<k-l;lp++)
    
81. 
    
82.        {
    
83. 
    
84.          lp1=lp+1;
    
85. 
    
86.          if(fxk[lp]<=fxk[lp1])
    
87. 
    
88.          {
    
89. 
    
90.     w=fxk[lp];
    
91. 
    
92.     fxk[lp]=fxk[lp1];
    
93. 
    
94.     fxk[lp1]=w;
    
95. 
    
96.     for(i=0;i<n;i++)
    
97. 
    
98.     {
    
99. 
    
100.       x[i]=xcom[i][lp];
     
101. 
     
102.       xcom[i][lp]=xcom[i][lp1];
     
103. 
     
104.       xcom[i][lp1]=x[i];
     
105. 
     
106.     }
     
107. 
     
108.          }
     
109. 
     
110.        }
     
111. 
     
112. }
     
113. 
     
114. 
     
115. void complex(int n,int k,int kg,float ep,float x[],float bl[],float bu[],
     
116. 
     
117.         float xcom[][100],float *f,int *nf,int *ng)
     
118. 
     
119. {
     
120. 
     
121.        int i,iw,l,ll,lh,it;
     
122. 
     
123.        float fx,fx0,sdx,fxh,fxr,alp;
     
124. 
     
125.        float *x0=(float*)calloc(n,sizeof(float));
     
126. 
     
127.        float *xh=(float*)calloc(n,sizeof(float));
     
128. 
     
129.        float *xr=(float*)calloc(n,sizeof(float));
     
130. 
     
131.        float *fxk=(float*)calloc(k,sizeof(float));
     
132. 
     
133.        float *g=(float*)calloc(kg,sizeof(float));
     
134. 
     
135. s5:    for(i=0;i<n;i++)
     
136. 
     
137.        x[i]=bl[i]+rand()/40000.0*(bu[i]-bl[i]);
     
138. 
     
139.        iw=gau(x,g,kg);
     
140. 
     
141.        *ng=*ng+1;
     
142. 
     
143.        if(iw==0)goto s5;
     
144. 
     
145.        for(i=0;i<n;i++)
     
146. 
     
147.        xcom[i][0]=x[i];
     
148. 
     
149.        for(l=1;l<k;l++)
     
150. 
     
151.        for(i=0;i<n;i++)
     
152. 
     
153.        xcom[i][l]=bl[i]+rand()/50000.0*(bu[i]-bl[i]);
     
154. 
     
155.        lh=-1;
     
156. 
     
157.        for(ll=1;ll<k;ll++)
     
158. 
     
159.        {
     
160. 
     
161.          xcent(n,k,ll,lh,x0,xcom);
     
162. 
     
163.          iw=gau(x0,g,kg);
     
164. 
     
165.          *ng=*ng+1;
     
166. 
     
167.          if(iw==0)goto s5;
     
168. 
     
169.          for(i=0;i<n;i++)
     
170. 
     
171.          x[i]=xcom[i][ll+1];
     
172. 
     
173. s24:     iw=gau(x,g,kg);
     
174. 
     
175.          *ng=*ng+1;
     
176. 
     
177.          if(iw==0)
     
178. 
     
179.          {
     
180. 
     
181.     for(i=0;i<n;i++)
     
182. 
     
183.     x[i]=x0[i]+0.5*(x[i]-x0[i]);
     
184. 
     
185.     goto s24;
     
186. 
     
187.          }
     
188. 
     
189.          else
     
190. 
     
191.          {
     
192. 
     
193.     for(i=0;i<n;i++)
     
194. 
     
195.     xcom[i][ll+1]=x[i];
     
196. 
     
197.          }
     
198. 
     
199.        }
     
200. 
     
201.        for(l=0;l<k;l++)
     
202. 
     
203.        {
     
204. 
     
205.          for(i=0;i<n;i++)
     
206. 
     
207.          x[i]=xcom[i][l];
     
208. 
     
209.          fx=objfx(x);
     
210. 
     
211.          *nf=*nf+1;
     
212. 
     
213.          fxk[l]=fx;
     
214. 
     
215.        }
     
216. 
     
217.        it=0;
     
218. 
     
219. s14: it=it+1;
     
220. 
     
221.        printf("\n\n    +++ ITERATION COMPUTE +++");
     
222. 
     
223.        printf("\n       ITER    = %d",it);
     
224. 
     
225.        lh=-1;
     
226. 
     
227.        xcent(n,k,k,lh,x0,xcom);
     
228. 
     
229.        fx0=objfx(x0);
     
230. 
     
231.        *nf=*nf+1;
     
232. 
     
233.        iw=gau(x0,g,kg);
     
234. 
     
235.        *ng=*ng+1;
     
236. 
     
237.        printf("\n       Fmid    = %f",fx0);
     
238. 
     
239.        for(i=0;i<n;i++)
     
240. 
     
241.          printf("\n       X(%d)mid=%f",i,x0[i]);
     
242. 
     
243.        for(i=0;i<kg;i++)
     
244. 
     
245.          printf("\n       G(%d)mid=%f",i,g[i]);
     
246. 
     
247.        sdx=0;
     
248. 
     
249.        for(l=0;l<k;l++)
     
250. 
     
251.          sdx=sdx+(fx0-fxk[l])*(fx0-fxk[l]);
     
252. 
     
253.        sdx=sqrt(sdx/(float)k);
     
254. 
     
255.        if(sdx<ep)goto s38;
     
256. 
     
257.        fxse(n,k,x,xcom,fxk);
     
258. 
     
259.        lh=0;
     
260. 
     
261. s22: fxh=fxk[lh];
     
262. 
     
263.        for(i=0;i<n;i++)
     
264. 
     
265.          xh[i]=xcom[i][lh];
     
266. 
     
267.        xcent(n,k,k,lh,x0,xcom);
     
268. 
     
269.        iw=gau(x0,g,kg);
     
270. 
     
271.        *ng=*ng+1;
     
272. 
     
273.        if(iw==0)goto s36;
     
274. 
     
275.        alp=1.3;
     
276. 
     
277. s12: for(i=0;i<n;i++)
     
278. 
     
279.        xr[i]=x0[i]+alp*(x0[i]-xh[i]);
     
280. 
     
281.        iw=gau(xr,g,kg);
     
282. 
     
283.        *ng=*ng+1;
     
284. 
     
285.        if(iw==0)
     
286. 
     
287.        {
     
288. 
     
289.          alp=alp*0.5;
     
290. 
     
291.          goto s12;
     
292. 
     
293.        }
     
294. 
     
295.        fxr=objfx(xr);
     
296. 
     
297.        *nf=*nf+1;
     
298. 
     
299.        if(fxr>=fxh)
     
300. 
     
301.        {
     
302. 
     
303.          if(alp>1.0e-4)
     
304. 
     
305.          {
     
306. 
     
307.     alp=alp*0.5;
     
308. 
     
309.     goto s12;
     
310. 
     
311.          }
     
312. 
     
313.          lh=lh+1;
     
314. 
     
315.          if(lh<3)goto s22;
     
316. 
     
317.        }
     
318. 
     
319.        for(i=0;i<n;i++)
     
320. 
     
321.        xcom[i][lh]=xr[i];
     
322. 
     
323.        fxk[lh]=fxr;
     
324. 
     
325.        goto s14;
     
326. 
     
327. s36: for(i=0;i<n;i++)
     
328. 
     
329.        {
     
330. 
     
331.          bl[i]=xcom[i][k];
     
332. 
     
333.          bu[i]=x0[i];
     
334. 
     
335.        }
     
336. 
     
337.        goto s5;
     
338. 
     
339. s38: for(i=0;i<n;i++)
     
340. 
     
341.        x[i]=x0[i];
     
342. 
     
343.        *f=objfx(x);
     
344. 
     
345.        *nf=*nf+1;
     
346. 
     
347.        free(x0);
     
348. 
     
349.        free(xh);
     
350. 
     
351.        free(xr);
     
352. 
     
353.        free(g);
     
354. 
     
355.        free(fxk);
     
356. 
     
357. }

复制代码

来自:http://hi.baidu.com/mec_auto/blo ... 999f305d600850.html