58-带罚函数的自适应粒子群算法

Halcom

百度网盘链接：
视频讲解链接：https://pan.baidu.com/s/1Xs--i-MwTU-_YUchq-iPPQ
https://pan.baidu.com/s/1xVCDEpAM2QXjC2Mvy5LdgA

具体链接在halcom.cn论坛，联系人QQ：3283892722
该论坛是一个学习交流平台，我会逐一的和大家分享学习。
欢迎大家录制视频，并提交给我，我来设置视频，你可在论坛进行打赏分享。
视频专用播放器：http://halcom.cn/forum.php?mod=viewthread&tid=258&extra=page%3D1


带罚函数的自适应粒子群算法：简单示例，

1. % 改进的快速粒子群优化算法 (APSO):
   
2. function apso
   
3. 
   
4. Lb=[0.1 0.1  0.1  0.1]; %下边界
   
5. Ub=[2.0 10.0 10.0 2.0]; %上边界
   
6. % 默认参数
   
7. para=[25 150 0.95]; %[粒子数，迭代次数，gama参数]
   
8. 
   
9. % APSO 优化求解函数
   
10. [gbest,fmin]=pso_mincon(@cost,@constraint,Lb,Ub,para);
    
11. 
    
12. % 输出结果
    
13. Bestsolution=gbest % 全局最优个体
    
14. fmin
    
15. 
    
16. %% 目标函数
    
17. function f=cost(x)
    
18. f=1.10471*x(1)^2*x(2)+0.04811*x(3)*x(4)*(14.0+x(2));
    
19. 
    
20. % 非线性约束
    
21. function [g,geq]=constraint(x)
    
22. % 不等式限制条件
    
23. Q=6000*(14+x(2)/2);
    
24. D=sqrt(x(2)^2/4+(x(1)+x(3))^2/4);
    
25. J=2*(x(1)*x(2)*sqrt(2)*(x(2)^2/12+(x(1)+x(3))^2/4));
    
26. alpha=6000/(sqrt(2)*x(1)*x(2));
    
27. beta=Q*D/J;
    
28. tau=sqrt(alpha^2+2*alpha*beta*x(2)/(2*D)+beta^2);
    
29. sigma=504000/(x(4)*x(3)^2);
    
30. delta=65856000/(30*10^6*x(4)*x(3)^3);
    
31. F=4.013*(30*10^6)/196*sqrt(x(3)^2*x(4)^6/36)*(1-x(3)*sqrt(30/48)/28);
    
32. 
    
33. g(1)=tau-13600;
    
34. g(2)=sigma-30000;
    
35. g(3)=x(1)-x(4);
    
36. g(4)=0.10471*x(1)^2+0.04811*x(3)*x(4)*(14+x(2))-5.0;
    
37. g(5)=0.125-x(1);
    
38. g(6)=delta-0.25;
    
39. g(7)=6000-F;
    
40. % 如果没有等式约束，则置geq=[]；
    
41. geq=[];
    
42. 
    
43. %%  APSO Solver
    
44. function [gbest,fbest]=pso_mincon(fhandle,fnonlin,Lb,Ub,para)
    
45. if nargin<=4,
    
46.     para=[20 150 0.95];
    
47. end
    
48. n=para(1);% 粒子种群大小
    
49. time=para(2); %时间步长，迭代次数
    
50. gamma=para(3); %gama参数
    
51. scale=abs(Ub-Lb); %取值区间
    
52. % 验证约束条件是否合乎条件
    
53. if abs(length(Lb)-length(Ub))>0,
    
54.     disp('Constraints must have equal size');
    
55.     return
    
56. end
    
57. 
    
58.   alpha=0.2; % alpha=[0,1]粒子随机衰减因子
    
59.   beta=0.5;  % 收敛速度(0->1)=(slow->fast);
    
60. 
    
61. % 初始化粒子群
    
62. best=init_pso(n,Lb,Ub);
    
63. 
    
64. fbest=1.0e+100;
    
65. % 迭代开始
    
66. for t=1:time,     
    
67.    
    
68. %寻找全局最优个体
    
69.   for i=1:n,   
    
70.     fval=Fun(fhandle,fnonlin,best(i,:));
    
71.     % 更新最有个体
    
72.     if fval<=fbest,
    
73.         gbest=best(i,:);
    
74.         fbest=fval;
    
75.     end
    
76.         
    
77.   end
    
78. 
    
79. % 随机性衰减因子
    
80. alpha=newPara(alpha,gamma);
    
81. 
    
82. % 更新粒子位置
    
83.   best=pso_move(best,gbest,alpha,beta,Lb,Ub);  
    
84. 
    
85. % 结果显示
    
86.         str=strcat('Best estimates: gbest=',num2str(gbest));
    
87.         str=strcat(str,'  iteration='); str=strcat(str,num2str(t));
    
88.         disp(str);
    
89. 
    
90.     fitness1(t)=fbest;
    
91.     plot(fitness1,'r','Linewidth',2)
    
92.     grid on
    
93.     hold on
    
94.     title('适应度')
    
95. end
    
96. 
    
97. % 初始化粒子函数
    
98. function [guess]=init_pso(n,Lb,Ub)
    
99. ndim=length(Lb);
    
100. for i=1:n,
     
101.     guess(i,1:ndim)=Lb+rand(1,ndim).*(Ub-Lb);
     
102. end
     
103. 
     
104. %更新所有的粒子 toward (xo,yo)
     
105. function ns=pso_move(best,gbest,alpha,beta,Lb,Ub)
     
106. % 增加粒子在上下边界区间内的随机性
     
107. n=size(best,1); ndim=size(best,2);
     
108. scale=(Ub-Lb);
     
109. for i=1:n,
     
110.     ns(i,:)=best(i,:)+beta*(gbest-best(i,:))+alpha.*randn(1,ndim).*scale;
     
111. end
     
112. ns=findrange(ns,Lb,Ub);
     
113. 
     
114. % 边界函数
     
115. function ns=findrange(ns,Lb,Ub)
     
116. n=length(ns);
     
117. for i=1:n,
     
118.   % 下边界约束
     
119.   ns_tmp=ns(i,:);
     
120.   I=ns_tmp<Lb;
     
121.   ns_tmp(I)=Lb(I);
     
122.   
     
123.   % 上边界约束
     
124.   J=ns_tmp>Ub;
     
125.   ns_tmp(J)=Ub(J);
     
126.   
     
127.   %更新粒子
     
128.   ns(i,:)=ns_tmp;
     
129. end
     
130. 
     
131. % 随机性衰减因子
     
132. function alpha=newPara(alpha,gamma);
     
133. alpha=alpha*gamma;
     
134. 
     
135. % 带约束的d维目标函数的求解
     
136. function z=Fun(fhandle,fnonlin,u)
     
137. % 目标
     
138. z=fhandle(u);
     
139. 
     
140. z=z+getconstraints(fnonlin,u); % 非线性约束
     
141. 
     
142. function Z=getconstraints(fnonlin,u)
     
143. % 罚常数 >> 1
     
144. PEN=10^15;
     
145. lam=PEN; lameq=PEN;
     
146. 
     
147. Z=0;
     
148. % 非线性约束
     
149. [g,geq]=fnonlin(u);
     
150. 
     
151. %通过不等式约束建立罚函数
     
152. for k=1:length(g),
     
153.     Z=Z+ lam*g(k)^2*getH(g(k));
     
154. end
     
155. % 等式条件约束
     
156. for k=1:length(geq),
     
157.    Z=Z+lameq*geq(k)^2*geteqH(geq(k));
     
158. end
     
159. 
     
160. % Test if inequalities
     
161. function H=getH(g)
     
162. if g<=0,
     
163.     H=0;
     
164. else
     
165.     H=1;
     
166. end
     
167. 
     
168. % Test if equalities hold
     
169. function H=geteqH(g)
     
170. if g==0,
     
171.     H=0;
     
172. else
     
173.     H=1;
     
174. end
     

复制代码

参考：
《MATLAB优化算法案例分析与应用》第13章（有约束函数极值APSO寻优）
《MATLAB优化算法案例分析与应用（进阶篇）》第18章（汽车动力传动参数优化设计）

德帝

想问下，这段代码有些地方看不懂，有讲解资料或者视频吗，谢谢

Halcom

德帝 发表于 2017-12-1 19:34
想问下，这段代码有些地方看不懂，有讲解资料或者视频吗，谢谢

代码先单步调试看看

季天翼

机器码：FF61C512CD36ED7C-0-EB3A9010F521F17C

季天翼

机器码：FF61C512CD36ED7C-0-EB3A9010F521F17C

季天翼

机器码：FF61C512CD36ED7C-0-EB3A9010F521F17C

jackian

非常好的学习资源，感谢分享

jackian

机器码：FF48E712CD53F57C-0-FF48C202D90BF17C

hou666666

非常好的学习资源，感谢分享

sunlei0417

观摩一下大佬写的代码