请选择 进入手机版 | 继续访问电脑版

Hello Mat

 找回密码
 立即注册
查看: 4438|回复: 0

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

[复制链接]

1277

主题

1500

帖子

85

金钱

管理员

Rank: 9Rank: 9Rank: 9

积分
22528
发表于 2019-8-11 14:40:23 | 显示全部楼层 |阅读模式
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. [bestf, bestindex]=min(fitness);
  55. zbest=pop(bestindex,:);   % 全局最佳
  56. gbest=pop;                % 个体最佳
  57. fitnessgbest=fitness;              % 个体最佳适应值
  58. fitnesszbest=bestf;               % 全局最佳适应值

  59. %% 迭代寻优
  60. iter = 0;
  61. y_fitness = zeros(1,maxiter );   % 预先产生4个空矩阵
  62. K_p = zeros(1,maxiter );
  63. K_i = zeros(1,maxiter );
  64. K_d = zeros(1,maxiter );
  65. while( (iter < maxiter ) && (fitnesszbest > MinFit) )
  66.     disp(['当前迭代次数: ', num2str(iter) ])
  67.     for j=1:sizepop
  68.         % 速度更新
  69.         % 引入模拟退火算法思想
  70.         [minF,minindex] = min(fitness);
  71.         Pg = pop(minindex,:);      % 从所有粒子中选取一个全局最优的某个粒子
  72.         deltaE = fitness(j) - minF;
  73.         if exp(-deltaE/T)>rand     % 以exp(-dC/T)概率接受新个体
  74.             V(j,:) = V(j,:) + c1*rand*(gbest(j,:) - pop(j,:)) + c2*rand*(zbest - pop(j,:));
  75.         else
  76.             V(j,:) = V(j,:) + c1*rand*(gbest(j,:) - pop(j,:)) + c2*rand*(Pg - pop(j,:));
  77.         end
  78.         % 速度越界限制
  79.         if V(j,:)>Vmax,
  80.             V(j,:)=Vmax;
  81.         end
  82.         if V(j,:)<Vmin,
  83.             V(j,:)=Vmin;
  84.         end
  85.         
  86.         % 位置更新
  87.         pop(j,:)=pop(j,:)+w*V(j,:);
  88.         pop(j,:) = lb_ub(pop(j,:),Lb,Ub);
  89.         % 适应值
  90. %         [fitness(j), fitB(j) ]=fun(pop(j,:));
  91.         fitness(j) = fun(pop(j,:));
  92.         
  93.         % 压缩因子
  94.         C = Cs - (Cs-Ce)*(iter/maxiter)^2;
  95.         xgima = 2/abs( 2-C - sqrt(C^2-4*C+eps) );
  96.         newpop = pop(j,:)+xgima*V(j,:);
  97.         newpop = lb_ub(newpop,Lb,Ub);
  98. %         [newfitness, newfitB ]= fun(newpop);
  99.         newfitness = fun(newpop);
  100.         if(newfitness<fitness(j))
  101.             pop(j,:) = newpop;
  102.             fitness(j) = newfitness;
  103. %             fitB(j) = newfitB;
  104.         end

  105.         % 个体最优更新
  106.         if fitness(j) < fitnessgbest(j)
  107.             gbest(j,:) = pop(j,:);
  108.             fitnessgbest(j) = fitness(j);
  109.         end
  110.         % 群体最优更新
  111.         if fitness(j) < fitnesszbest
  112.             zbest = pop(j,:);
  113.             fitnesszbest = fitness(j);
  114. %             fitBzbest = fitB(j);
  115.         end
  116.         
  117.         % 引入跳跃算子思想
  118.         Pm = exp( -(fitnessgbest(j)-fitnesszbest)/T );
  119.         if(Pm>rand)
  120.             newpop = pop(j,:) + 0.5*rand*(Ub-Lb);
  121.         end
  122.         newpop = lb_ub(newpop,Lb,Ub);
  123. %         [newfitness, newfitB ]= fun(newpop);
  124.         newfitness = fun(newpop);
  125.         if(newfitness<fitnesszbest)
  126.             pop(j,:) = newpop;
  127.             fitness(j) = newfitness;
  128.             zbest = newpop;
  129.             fitnesszbest = newfitness;
  130. %             fitBzbest = newfitB;
  131.         end
  132.         
  133.     end
  134.     iter = iter+1;                      % 迭代次数更新
  135.     y_fitness(1,iter) = fitnesszbest;         % 为绘图做准备
  136.     K_p(1,iter) = zbest(1);
  137.     K_i(1,iter) = zbest(2);
  138.     K_d(1,iter) = zbest(3);
  139.     T = q*T;   % 更新当前温度
  140. end
  141. %% 绘图
  142. disp(['IPSO最优解', num2str(zbest)])
  143. disp(['IPSO最优解对应的目标值', num2str(fitnesszbest)])
  144. fprintf('\n')

  145. figure(1)      % 绘制性能指标ITAE的变化曲线
  146. plot(y_fitness,'LineWidth',3)
  147. title('最优个体适应值','fontsize',10);
  148. xlabel('迭代次数','fontsize',10);ylabel('适应值','fontsize',10);
  149. set(gca,'Fontsize',10);
  150. grid on

  151. figure(2)      % 绘制PID控制器参数变化曲线
  152. plot(K_p,'b-','LineWidth',2); hold on
  153. plot(K_i,'k-.','LineWidth',2)
  154. plot(K_d,'r--','LineWidth',2)
  155. title('Kp、Ki、Kd 优化曲线','fontsize',10);
  156. xlabel('迭代次数','fontsize',10);   ylabel('参数值','fontsize',10);
  157. set(gca,'Fontsize',10);
  158. legend('Kp','Ki','Kd',1);
  159. grid on;hold off;

  160. % PID响应
  161. if choosed==1
  162.     [rint, yout]=PID_Fun_1_response(zbest);
  163. elseif(choosed==2)
  164.     [rint, yout]=PID_Fun_2_response(zbest);
  165. elseif(choosed==3)
  166.     [rint, yout]=PID_Fun_3_response(zbest);
  167. elseif(choosed==4)
  168.     [rint, yout]=PID_Fun_4_response(zbest);
  169. elseif(choosed==5)
  170.     [timef, rint, yout]=PID_Fun_5_response(zbest);
  171. end
  172. figure('color',[1,1,1])
  173. plot(timef, rint,'b-','linewidth',2);hold on;
  174. plot(timef, yout,'r.-','linewidth',2);
  175. axis tight;grid on;
  176. hold off;

  177. save IPSO.mat
  178. rmpath(genpath('./PID_funtion'))
  179. % 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. load('PSO.mat')
  5. figure(1)      % 绘制性能指标ITAE的变化曲线
  6. plot(y_fitness,'b-', 'LineWidth',3); hold on;
  7. clear
  8. load('IPSO.mat')
  9. plot(y_fitness,'r-', 'LineWidth',3)
  10. legend('PSO', '改进的PSO')
  11. title('最优个体适应值','fontsize',10);
  12. xlabel('迭代次数','fontsize',10);ylabel('适应值','fontsize',10);
  13. grid on

  14. clear
  15. load('PSO.mat')
  16. figure('color',[1,1,1])
  17. plot(timef, yout,'b.-','linewidth',2);hold on;
  18. axis tight;grid on;
  19. clear
  20. load('IPSO.mat')
  21. plot(timef, yout,'r.-','linewidth',2);
  22. plot(timef, rint,'g-','linewidth',2);
  23. legend('PSO', '改进的PSO')
  24. title('阶跃响应曲线')
  25. xlabel('(s)'); ylabel('响应值')
  26. hold off;
复制代码







本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
算法QQ  3283892722
群智能算法链接http://halcom.cn/forum.php?mod=forumdisplay&fid=73
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

Python|Opencv|MATLAB|Halcom.cn ( 蜀ICP备16027072号 )

GMT+8, 2024-3-5 01:45 , Processed in 0.201994 second(s), 26 queries .

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表