机械优化设计程序

ttksd

　　今天终于交了机械优化设计作业，程序贴出来，那位要用就拿去吧，资源共享了，O(∩_∩)O
　　下面是利用外点惩罚函数(Powell)法求解三维目标函数最优解与最优值的程序，用C++编写。

1. 　　#include
   
2. 
   
3. 　　#include
   
4. 
   
5. 　　double lamta[10]={0, 1.0 ,0 ,0 ,0 ,1 ,0 ,0 ,0 ,1};//鲍威尔方法初始化方向，线性无关
   
6. 
   
7. 　　double lamta1[3]={0, 0 , 0};//暂存新的搜索方向
   
8. 
   
9. 　　double x1[4]={0, 0 ,0, 0 };//x1到x3用于存储各共轭方向的点
   
10. 
    
11. 　　double x2[4]={0, 0 ,0, 0 };
    
12. 
    
13. 　　double x3[4]={0, 0 ,0, 0 };
    
14. 
    
15. 　　double x4[4]={0, 0 ,0, 0 };//x4用于中间判断
    
16. 
    
17. 　　double x5[4]={0, 0 ,0, 0 };//x5用存放于更换方向后产生的新点
    
18. 
    
19. 　　int m=0;//标志
    
20. 
    
21. 　　double x_[4]={0, 0, 0, 0};//暂存鲍威尔最优解
    
22. 
    
23. 　　double x0[4]={0, 2, 2 , 2};//初值
    
24. 
    
25. 　　double c=10;//递减系数
    
26. 
    
27. 　　double e=0.00001;//精度控制
    
28. 
    
29. 　　double r0=1;//初始惩罚因子
    
30. 
    
31. 　　double r=1;
    
32. 
    
33. 　　//函数声明部分
    
34. 
    
35. 　　void Powell(double r); //鲍威尔方法函数
    
36. 
    
37. 　　double fxy(double x1,double x2,double x3,double r); //待求函数
    
38. 
    
39. 　　double ysearch(double x); //一维搜索的目标函数
    
40. 
    
41. 　　void search(double &a,double &b,double h); //区间搜索
    
42. 
    
43. 　　double yellowcut(double &a,double &b); //黄金分割
    
44. 
    
45. 　　void sort(double *p,int size);//选择法排序
    
46. 
    
47. 　　void main() //约束优化方法主函数入口
    
48. 
    
49. 　　{
    
50. 
    
51. 　　cout<<"请输入精度"<
    
52. 
    
53. 　　cin>>e;
    
54. 
    
55. 　　changyan:Powell(r);
    
56. 
    
57. 　　double cmpare[4];
    
58. 
    
59. 　　int flag1=0;
    
60. 
    
61. 　　for (int i=1;i<=3;i++)
    
62. 
    
63. 　　{
    
64. 
    
65. 　　cmpare[i]=x_[i]-x0[i];
    
66. 
    
67. 　　if (fabs(cmpare[i])
    
68. 
    
69. 　　{flag1++;}
    
70. 
    
71. 　　}
    
72. 
    
73. 　　if (flag1==3)
    
74. 
    
75. 　　{
    
76. 
    
77. 　　cout<<"最优解为："<<"x1="<
    
78. 
    
79. 　　}
    
80. 
    
81. 　　else
    
82. 
    
83. 　　{
    
84. 
    
85. 　　for (int j=1;j<=3;j++)
    
86. 
    
87. 　　{
    
88. 
    
89. 　　x0[j]=x_[j];
    
90. 
    
91. 　　}
    
92. 
    
93. 　　r=c*r;
    
94. 
    
95. 　　goto changyan;
    
96. 
    
97. 　　}
    
98. 
    
99. 　　}
    
100. 
     
101. 　　//子函数定义部分
     
102. 
     
103. 　　double fxy(double x1,double x2,double x3,double r)//待求函数
     
104. 
     
105. 　　{
     
106. 
     
107. 　　double m,n,p;
     
108. 
     
109. 　　m=(-x1>0)?(-x1):0;
     
110. 
     
111. 　　n=(-x2>0)?(-x2):0;
     
112. 
     
113. 　　p=(-x3>0)?(-x3):0;
     
114. 
     
115. 　　return //惩罚函数
     
116. 
     
117. 　　1000-x1*x1-2*x2*x2-x3*x3-x1*x2-x1*x3+r*(m*m+n*n+p*p)+r*((x1*x1+x2*x2+x3*x3-25)*(x1*x1+x2*x2+x3*x3-25)+(8*x1+14*x2+7*x3-56)*(8*x1+14*x2+7*x3-56));
     
118. 
     
119. 　　}
     
120. 
     
121. 　　void Powell(double r) //鲍威尔方法函数定义
     
122. 
     
123. 　　{
     
124. 
     
125. 　　double det=0.0001; //迭代精度
     
126. 
     
127. 　　int k;
     
128. 
     
129. 　　my1: for (k=1;k<=3;k++)
     
130. 
     
131. 　　{
     
132. 
     
133. 　　m=3*k-2;
     
134. 
     
135. 　　double a=0,b=0,xo=0;
     
136. 
     
137. 　　search(a,b,1); //完成区间搜索
     
138. 
     
139. 　　double temp;
     
140. 
     
141. 　　temp=yellowcut(a,b);//黄金分割法
     
142. 
     
143. 　　int n=3*k-2;
     
144. 
     
145. 　　for (int i=1;i<=3;i++)
     
146. 
     
147. 　　{
     
148. 
     
149. 　　switch (k)
     
150. 
     
151. 　　{
     
152. 
     
153. 　　case 1:x1[i]=x0[i]+temp*lamta[n++];break;
     
154. 
     
155. 　　case 2:x2[i]=x1[i]+temp*lamta[n++];break;
     
156. 
     
157. 　　case 3:x3[i]=x2[i]+temp*lamta[n++];break;
     
158. 
     
159. 　　default :break;
     
160. 
     
161. 　　}
     
162. 
     
163. 　　}
     
164. 
     
165. 　　}
     
166. 
     
167. 　　double cmp[4];
     
168. 
     
169. 　　int flag=0;
     
170. 
     
171. 　　for (int i=1;i<=3;i++)
     
172. 
     
173. 　　{
     
174. 
     
175. 　　cmp[i]=x3[i]-x0[i];
     
176. 
     
177. 　　if (fabs(cmp[i])
     
178. 
     
179. 　　{flag++;}}
     
180. 
     
181. 　　if (flag==3) //找到最优解
     
182. 
     
183. 　　{
     
184. 
     
185. 　　x_[1]=x3[1];
     
186. 
     
187. 　　x_[2]=x3[2];
     
188. 
     
189. 　　x_[3]=x3[3];
     
190. 
     
191. 　　}
     
192. 
     
193. 　　else
     
194. 
     
195. 　　{
     
196. 
     
197. 　　double fy[4];
     
198. 
     
199. 　　fy[0]=fxy(x0[1],x0[2],x0[3],r);
     
200. 
     
201. 　　fy[1]=fxy(x1[1],x1[2],x1[3],r);
     
202. 
     
203. 　　fy[2]=fxy(x2[1],x2[2],x2[3],r);
     
204. 
     
205. 　　fy[3]=fxy(x3[1],x3[2],x3[3],r); double fyy[3];
     
206. 
     
207. 　　for (int ii=0;ii<3;ii++)
     
208. 
     
209. 　　{fyy[ii]=fy[ii]-fy[ii+1];}
     
210. 
     
211. 　　sort(fyy,3);
     
212. 
     
213. 　　for (ii=1;ii<=3;ii++)
     
214. 
     
215. 　　{x4[ii]=2*x3[ii]-x0[ii];}
     
216. 
     
217. 　　double f0,f3,f4;
     
218. 
     
219. 　　f0=fy[0];
     
220. 
     
221. 　　f3=fy[3];
     
222. 
     
223. 　　f4=fxy(x4[1],x4[2],x4[3],r);
     
224. 
     
225. 　　if ((f0+f4-2*f3)/2>=fyy[2])
     
226. 
     
227. 　　{
     
228. 
     
229. 　　if (f3
     
230. 
     
231. 　　{
     
232. 
     
233. 　　for (int t=1;t<=3;t++)
     
234. 
     
235. 　　{x0[t]=x3[t];}
     
236. 
     
237. 　　}
     
238. 
     
239. 　　else
     
240. 
     
241. 　　{
     
242. 
     
243. 　　for (int t=1;t<=3;t++)
     
244. 
     
245. 　　{x0[t]=x4[t];
     
246. 
     
247. 　　}}
     
248. 
     
249. 　　goto my1;
     
250. 
     
251. 　　}
     
252. 
     
253. 　　else{
     
254. 
     
255. 　　for (int t=0;t<3;t++)
     
256. 
     
257. 　　{lamta1[t]=x3[t+1]-x0[t+1];}
     
258. 
     
259. 　　m=0; //switch 标志!
     
260. 
     
261. 　　double aa=0,bb=0;
     
262. 
     
263. 　　search(aa,bb,1);
     
264. 
     
265. 　　double temp1;
     
266. 
     
267. 　　temp1=yellowcut(aa,bb);
     
268. 
     
269. 　　for (int i=1;i<=3;i++)
     
270. 
     
271. 　　{x5[i]=x3[i]+temp1*lamta1[i-1];}
     
272. 
     
273. 　　for (i=1;i<=3;i++)
     
274. 
     
275. 　　{x0[i]=x5[i];}
     
276. 
     
277. 　　for (i=1;i<=6;i++)
     
278. 
     
279. 　　{lamta[i]=lamta[i+3];}
     
280. 
     
281. 　　for (i=1;i<=3;i++)
     
282. 
     
283. 　　{
     
284. 
     
285. 　　lamta[6+i]=lamta1[i-1];}
     
286. 
     
287. 　　goto my1;
     
288. 
     
289. 　　}}}
     
290. 
     
291. 　　double ysearch(double x) //一维搜索的目标函数
     
292. 
     
293. 　　{
     
294. 
     
295. 　　switch (m)
     
296. 
     
297. 　　{
     
298. 
     
299. 　　case 1: return fxy(x0[1]+x*lamta[m],x0[2]+x*lamta[m+1],x0[3]+x*lamta[m+2],r);break;
     
300. 
     
301. 　　case 4: return fxy(x1[1]+x*lamta[m],x1[2]+x*lamta[m+1],x1[3]+x*lamta[m+2],r);break;
     
302. 
     
303. 　　case 7: return fxy(x2[1]+x*lamta[m],x2[2]+x*lamta[m+1],x2[3]+x*lamta[m+2],r);break;
     
304. 
     
305. 　　case 0: return fxy(x3[1]+x*lamta1[0],x3[2]+x*lamta1[1],x3[3]+x*lamta1[2],r);break;//更改方向后的一维搜索
     
306. 
     
307. 　　default:return 0; break;
     
308. 
     
309. 　　}
     
310. 
     
311. 　　}
     
312. 
     
313. 　　void search(double &a,double &b,double h) //区间搜索
     
314. 
     
315. 　　{double a1,a2,a3,y1,y2,y3;
     
316. 
     
317. 　　h=1;
     
318. 
     
319. 　　a1=a,y1=ysearch(a1);
     
320. 
     
321. 　　a2=a+h,y2=ysearch(a2);
     
322. 
     
323. 　　if(y2>=y1){
     
324. 
     
325. 　　h=-h,a3=a1,y3=y1;
     
326. 
     
327. 　　a1=a2,y1=y2,a2=a3,y2=y3;}
     
328. 
     
329. 　　a3=a2+h,y3=ysearch(a3);
     
330. 
     
331. 　　while(y3<=y2){
     
332. 
     
333. 　　h=2*h;
     
334. 
     
335. 　　a1=a2,y1=y2,a2=a3,y2=y3;
     
336. 
     
337. 　　a3=a2+h,y3=ysearch(a3);
     
338. 
     
339. 　　}
     
340. 
     
341. 　　if(h<0)a=a3,b=a1;
     
342. 
     
343. 　　else a=a1,b=a3;}
     
344. 
     
345. 　　double yellowcut(double &a,double &b){
     
346. 
     
347. 　　double e; //黄金分割法求解
     
348. 
     
349. 　　e=0.001;
     
350. 
     
351. 　　double c,fc;
     
352. 
     
353. 　　c=a+0.382*(b-a);
     
354. 
     
355. 　　fc=ysearch(c);
     
356. 
     
357. 　　double d,fd;
     
358. 
     
359. 　　double xo;
     
360. 
     
361. 　　d=a+0.618*(b-a);
     
362. 
     
363. 　　fd=ysearch(d);
     
364. 
     
365. 　　label2: if (fc<=fd)
     
366. 
     
367. 　　{b=d;
     
368. 
     
369. 　　d=c;
     
370. 
     
371. 　　fd=fc;
     
372. 
     
373. 　　c=a+0.382*(b-a);
     
374. 
     
375. 　　fc=ysearch(c);}
     
376. 
     
377. 　　else
     
378. 
     
379. 　　{a=c;
     
380. 
     
381. 　　c=d;
     
382. 
     
383. 　　fc=fd;
     
384. 
     
385. 　　d=a+0.618*(b-a);
     
386. 
     
387. 　　fd=ysearch(d);}
     
388. 
     
389. 　　if ((b-a)<=e)
     
390. 
     
391. 　　{xo=(a+b)/2;}
     
392. 
     
393. 　　else
     
394. 
     
395. 　　goto label2;
     
396. 
     
397. 　　return xo;
     
398. 
     
399. 　　}
     
400. 
     
401. 　　void sort(double *p,int size){//选择法排序
     
402. 
     
403. 　　int i,j;
     
404. 
     
405. 　　double k;
     
406. 
     
407. 　　for(i=0;i
     
408. 
     
409. 　　for(j=i+1;j
     
410. 
     
411. 　　if(*(p+i)>*(p+j)){k=*(p+i);*(p+i)=*(p+j);*(p+j)=k;}
     
412. 
     
413. 　　}
     

复制代码

转自：http://blog.sina.com.cn/s/blog_49d083e50100cojl.html