8-改进PSO的三阶系统PID整定

Halcom

8-改进PSO的三阶系统PID整定
百度网盘视频链接：https://pan.baidu.com/s/1GE72YSMENrHLcAJLx8rirg 提取码：nn5n
录制的视频是算法底层原理讲解，底层代码实现，方便大家真正掌握算法实质，开发出更加出色的算法。录制视频的初衷是：避免读者朋友利用大把时间学习已有常见算法，本系列视频旨在让读者朋友快速深入了解这些常见算法原理，有更多的时间去研究更加高大上算法（价值）。

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

1.约束因子的改进
带压缩因子的PSO优化算法，Clerc和Kennedy提出的带压缩因子的PSO优化算法通过选取合适参数，可确保PSO算法的收敛性，并可取消对速度的边界限制。速度和位置更新公式如下：
v_(i,j) (k+1)=χ[v_(i,j) (k)+c_1 r_1 (p_(i,j) (k)-x_(i,j) (k))+c_2 r_2 (p_(g,j) (k)-x_(i,j) (k))]
x_(i,j) (k+1)=x_(i,j) (k)+v_(i,j) (k+1),j=1,⋯,n
χ=2/|2-C-√(C^2-4C)|
C=c_1+c_2,C>4

2.学习因子的更新
粒子在迭代初期应该有较强的学习能力，所以学习因子应该随迭代次数的增加而逐渐减小，更新公式如下：
c=c_s-(c_s-c_e )*〖(iter/Maxiter)〗^2
在上述公式中，cs和ce分别是学习因子的上限和下限，学习因子的大小随迭代次数的增加呈抛物线式下降。

3.引入模拟退火算法思想
通过将模拟退火算法的机制融合，采用2轮盘赌输策略，从所有粒子中选取一个全局最优的某个粒子，来替代 Pg，记作 P'g． 计算当前温度T下，当前粒子 Pi和 P'g目标函数值之差 ΔE，exp(－ΔE /T) 即跳突概率．
引入模拟退火算法的概率判断 Metropolis 准则:
P=exp⁡(-∆E/T)
P 为接受差解的概率; ∆E为前后2代目标函数的差值；T为模拟退火算法中的温度参数，其随迭代次数的增加而不断下降。
速度更新公式如下:
v_(i,j) (k+1)=χ[v_(i,j) (k)+c_1 r_1 (p_(i,j) (k)-x_(i,j) (k))+c_2 r_2 (〖p'〗_(g,j) (k)-x_(i,j) (k))]
x_(i,j) (k+1)=x_(i,j) (k)+v_(i,j) (k+1),j=1,⋯,n

4.引入跳跃算子思想
提出一种自适应跳跃算子，通过比较pbest与gbest的相似程度，赋予粒子不同的跳跃概率。在t时刻，第i个粒子跳出当前位置的概率公式为：
P_m=exp⁡(f(〖gbest〗^t )-f(〖pbest_i〗^t))
跳跃公式为：
x_i^t=x_i^t+rand∙(Ub-Lb)
其中：D为搜索空间的维数；rand为[0,1]之间的随机数；Ub和Lb为问题的上下限。从公式中可以看出当粒子的pbest与种群gbest越接近时，粒子产生跳跃的概率较大，也就是说可以使停止更新的粒子，以较大的概率跳出极值点，对于距离种群 gbest 较远的粒子，跳跃算子对其的作用较小。

主程序如下：

1. %% 清空环境
   
2. clc,clear,close all
   
3. warning off
   
4. format longG
   
5. addpath(genpath('./PID_funtion'))
   
6. %% 参数设置
   
7. w = 0.6;      % 惯性因子
   
8. % PSO 参数
   
9. c1 = 2.4995;  
   
10. c2 = 2.4995;
    
11. C = c1+c2;
    
12. if C<=4
    
13.     msgbox('c1+c2<=4！！！！')
    
14.     return;
    
15. end
    
16. xgima = 2/abs( 2-C - sqrt(C^2-4*C+eps) );
    
17. Cs = 8;   % C的上限
    
18. Ce = 4;   % C的下限
    
19. Dim = 3;            % 维数
    
20. sizepop = 20;       % 粒子群规模
    
21. maxiter  = 30;      % 最大迭代次数
    
22. MinFit = 0.1;       % 最小适应值
    
23. Vmax = 1;           % 速度最大值
    
24. Vmin = -1;          % 速度最小值
    
25. %       P  I  D
    
26. Ub = [100 100 100];     % PID取值 最大值
    
27. Lb = [0.00 0.00 0.00];   % PID取值 最小值
    
28. %% 模拟退火算法
    
29. T = 500;          % initiation temperature
    
30. q = 0.5;          % 降温ratio
    
31. %% 适应度函数
    
32. choosed = 5;
    
33. if choosed==1
    
34.     fun = @(x)PID_Fun_1(x);
    
35. elseif(choosed==2)
    
36.     fun = @(x)PID_Fun_2(x);
    
37. elseif(choosed==3)
    
38.     fun = @(x)PID_Fun_3(x);
    
39. elseif(choosed==4)
    
40.     fun = @(x)PID_Fun_4(x);
    
41. elseif(choosed==5)
    
42.     fun = @(x)PID_Fun_5(x);
    
43. end
    
44. %% 粒子群初始化
    
45. Range = ones(sizepop,1)*(Ub-Lb);
    
46. pop = rand(sizepop,Dim).*Range + ones(sizepop,1)*Lb;    % 初始化粒子群
    
47. V = rand(sizepop,Dim)*(Vmax-Vmin) + Vmin;                 % 初始化速度
    
48. fitness = zeros(sizepop,1);
    
49. for i=1:sizepop
    
50. %     [fitness(i,:), fitB(i)] = fun(pop(i,:));                         % 粒子群的适应值
    
51.     fitness(i,:) = fun(pop(i,:));                         % 粒子群的适应值
    
52. end
    
53. 
    
54. %% 个体极值和群体极值
    
55. [bestf, bestindex]=min(fitness);
    
56. zbest=pop(bestindex,:);   % 全局最佳
    
57. gbest=pop;                % 个体最佳
    
58. fitnessgbest=fitness;              % 个体最佳适应值
    
59. fitnesszbest=bestf;               % 全局最佳适应值
    
60. 
    
61. %% 迭代寻优
    
62. iter = 0;
    
63. y_fitness = zeros(1,maxiter );   % 预先产生4个空矩阵
    
64. K_p = zeros(1,maxiter );
    
65. K_i = zeros(1,maxiter );
    
66. K_d = zeros(1,maxiter );
    
67. while( (iter < maxiter ) && (fitnesszbest > MinFit) )
    
68.     disp(['当前迭代次数: ', num2str(iter) ])
    
69.     for j=1:sizepop
    
70.         % 速度更新
    
71.         % 引入模拟退火算法思想
    
72.         [minF,minindex] = min(fitness);
    
73.         Pg = pop(minindex,:);      % 从所有粒子中选取一个全局最优的某个粒子
    
74.         deltaE = fitness(j) - minF;
    
75.         if exp(-deltaE/T)>rand     % 以exp(-dC/T)概率接受新个体
    
76.             V(j,:) = V(j,:) + c1*rand*(gbest(j,:) - pop(j,:)) + c2*rand*(zbest - pop(j,:));
    
77.         else
    
78.             V(j,:) = V(j,:) + c1*rand*(gbest(j,:) - pop(j,:)) + c2*rand*(Pg - pop(j,:));
    
79.         end
    
80.         % 速度越界限制
    
81.         if V(j,:)>Vmax,
    
82.             V(j,:)=Vmax;
    
83.         end
    
84.         if V(j,:)<Vmin,
    
85.             V(j,:)=Vmin;
    
86.         end
    
87.         
    
88.         % 位置更新
    
89.         pop(j,:)=pop(j,:)+w*V(j,:);
    
90.         pop(j,:) = lb_ub(pop(j,:),Lb,Ub);
    
91.         % 适应值
    
92. %         [fitness(j), fitB(j) ]=fun(pop(j,:));
    
93.         fitness(j) = fun(pop(j,:));
    
94.         
    
95.         % 压缩因子
    
96.         C = Cs - (Cs-Ce)*(iter/maxiter)^2;
    
97.         xgima = 2/abs( 2-C - sqrt(C^2-4*C+eps) );
    
98.         newpop = pop(j,:)+xgima*V(j,:);
    
99.         newpop = lb_ub(newpop,Lb,Ub);
    
100. %         [newfitness, newfitB ]= fun(newpop);
     
101.         newfitness = fun(newpop);
     
102.         if(newfitness<fitness(j))
     
103.             pop(j,:) = newpop;
     
104.             fitness(j) = newfitness;
     
105. %             fitB(j) = newfitB;
     
106.         end
     
107. 
     
108.         % 个体最优更新
     
109.         if fitness(j) < fitnessgbest(j)
     
110.             gbest(j,:) = pop(j,:);
     
111.             fitnessgbest(j) = fitness(j);
     
112.         end
     
113.         % 群体最优更新
     
114.         if fitness(j) < fitnesszbest
     
115.             zbest = pop(j,:);
     
116.             fitnesszbest = fitness(j);
     
117. %             fitBzbest = fitB(j);
     
118.         end
     
119.         
     
120.         % 引入跳跃算子思想
     
121.         Pm = exp( -(fitnessgbest(j)-fitnesszbest)/T );
     
122.         if(Pm>rand)
     
123.             newpop = pop(j,:) + 0.5*rand*(Ub-Lb);
     
124.         end
     
125.         newpop = lb_ub(newpop,Lb,Ub);
     
126. %         [newfitness, newfitB ]= fun(newpop);
     
127.         newfitness = fun(newpop);
     
128.         if(newfitness<fitnesszbest)
     
129.             pop(j,:) = newpop;
     
130.             fitness(j) = newfitness;
     
131.             zbest = newpop;
     
132.             fitnesszbest = newfitness;
     
133. %             fitBzbest = newfitB;
     
134.         end
     
135.         
     
136.     end
     
137.     iter = iter+1;                      % 迭代次数更新
     
138.     y_fitness(1,iter) = fitnesszbest;         % 为绘图做准备
     
139.     K_p(1,iter) = zbest(1);
     
140.     K_i(1,iter) = zbest(2);
     
141.     K_d(1,iter) = zbest(3);
     
142.     T = q*T;   % 更新当前温度
     
143. end
     
144. %% 绘图
     
145. disp(['IPSO最优解', num2str(zbest)])
     
146. disp(['IPSO最优解对应的目标值', num2str(fitnesszbest)])
     
147. fprintf('\n')
     
148. 
     
149. figure(1)      % 绘制性能指标ITAE的变化曲线
     
150. plot(y_fitness,'LineWidth',3)
     
151. title('最优个体适应值','fontsize',10);
     
152. xlabel('迭代次数','fontsize',10);ylabel('适应值','fontsize',10);
     
153. set(gca,'Fontsize',10);
     
154. grid on
     
155. 
     
156. figure(2)      % 绘制PID控制器参数变化曲线
     
157. plot(K_p,'b-','LineWidth',2); hold on
     
158. plot(K_i,'k-.','LineWidth',2)
     
159. plot(K_d,'r--','LineWidth',2)
     
160. title('Kp、Ki、Kd 优化曲线','fontsize',10);
     
161. xlabel('迭代次数','fontsize',10);   ylabel('参数值','fontsize',10);
     
162. set(gca,'Fontsize',10);
     
163. legend('Kp','Ki','Kd',1);
     
164. grid on;hold off;
     
165. 
     
166. % PID响应
     
167. if choosed==1
     
168.     [rint, yout]=PID_Fun_1_response(zbest);
     
169. elseif(choosed==2)
     
170.     [rint, yout]=PID_Fun_2_response(zbest);
     
171. elseif(choosed==3)
     
172.     [rint, yout]=PID_Fun_3_response(zbest);
     
173. elseif(choosed==4)
     
174.     [rint, yout]=PID_Fun_4_response(zbest);
     
175. elseif(choosed==5)
     
176.     [timef, rint, yout]=PID_Fun_5_response(zbest);
     
177. end
     
178. figure('color',[1,1,1])
     
179. plot(timef, rint,'b-','linewidth',2);hold on;
     
180. plot(timef, yout,'r.-','linewidth',2);
     
181. axis tight;grid on;
     
182. hold off;
     
183. 
     
184. save IPSO.mat
     
185. rmpath(genpath('./PID_funtion'))
     
186. % rmpath(genpath('./test_function'))

复制代码

取值范围约束

1. function pop = lb_ub(pop,popmin,popmax)
   
2. j=1;
   
3. for i=1:length(popmin)
   
4.     if pop(j,i)>popmax(i)
   
5.         pop(j,i)=popmax(i);
   
6.     end
   
7.     if pop(j,i)<popmin(i)
   
8.         pop(j,i)=popmin(i);
   
9.     end
   
10. end
    
11. end

复制代码

对比画图：

1. clc,clear,close all
   
2. warning off
   
3. format longG
   
4. 
   
5. load('PSO.mat')
   
6. figure(1)      % 绘制性能指标ITAE的变化曲线
   
7. plot(y_fitness,'b-', 'LineWidth',3); hold on;
   
8. clear
   
9. load('IPSO.mat')
   
10. plot(y_fitness,'r-', 'LineWidth',3)
    
11. legend('PSO', '改进的PSO')
    
12. title('最优个体适应值','fontsize',10);
    
13. xlabel('迭代次数','fontsize',10);ylabel('适应值','fontsize',10);
    
14. grid on
    
15. 
    
16. clear
    
17. load('PSO.mat')
    
18. figure('color',[1,1,1])
    
19. plot(timef, yout,'b.-','linewidth',2);hold on;
    
20. axis tight;grid on;
    
21. clear
    
22. load('IPSO.mat')
    
23. plot(timef, yout,'r.-','linewidth',2);
    
24. plot(timef, rint,'g-','linewidth',2);
    
25. legend('PSO', '改进的PSO')
    
26. title('阶跃响应曲线')
    
27. xlabel('(s)'); ylabel('响应值')
    
28. hold off;

复制代码