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

Hello Mat

 找回密码
 立即注册

QQ登录

只需一步,快速开始

查看: 294|回复: 0

基于细菌觅食算法的PID整定

[复制链接]

691

主题

813

帖子

2万

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
29493
发表于 2019-5-23 22:40:45 | 显示全部楼层 |阅读模式
基于细菌觅食算法的PID整定
  1. clc,clear,close all
  2. warning off
  3. format longG
  4. % BFO优化算法 参数
  5. Nc = 50;                   % 趋化次数
  6. Ns = 4;                    % 游动次数
  7. Nre = 4;                   % 复制次数
  8. Ned = 2;                   % 驱散(迁移)次数
  9. sizepop = 20;              % 种群数量
  10. Sr = ceil( sizepop/2 );    % 复制(分裂)次数
  11. Ped = 0.25;                             % 细菌驱散(迁移)概率
  12. nvar = 3;              % 3个未知量
  13. popmin = [0,0,0];      % x min
  14. popmax = [30,30,30];   % x max
  15. Cmin = [-2,-2,-2];     % 最小步长
  16. Cmax = [2,2,2];        % 最大步长
  17. C(:,:) = 0.001*ones(sizepop,nvar);            % 翻转选定方向后,单个细菌前进的步长
  18. %% 适应度函数,也就是我们的目标函数
  19. % fun = @(x)PID_Fun_1(x);
  20. fun = @(x)PID_Fun_2(x);
  21. %% 初始化种群
  22. for i=1:sizepop
  23.     pop(i,:) = popmin + (popmax-popmin).*rand(1,nvar);   % 初始化个体
  24.     fitness(i) = fun(  pop(i,:) );              % 适应度值
  25.     C(i,:) = Cmin + (Cmax-Cmin).*rand(1,nvar);  % 步长
  26. end
  27. %% 记录一组最优值
  28. [bestfitness,bestindex]=min(fitness);  % 目标值越小越好
  29. zbest=pop(bestindex,:);   % 全局最佳
  30. fitnesszbest=bestfitness; % 全局最佳适应度值
  31. NcSizepop = 0;            % 记录最优适应度值(函数值)
  32. %% 迭代寻优
  33. for i = 1:Ned                   % 驱散(迁移)次数
  34.     disp(['当前驱散(迁移)次数 i=', num2str(i), ' , 驱散(迁移)次数 Ned=', num2str(Ned)])
  35.     for k = 1:Nre               % 复制次数
  36.         disp(['复制次数 k=', num2str(k), ' , 复制次数 Nre=', num2str(Nre)])
  37.         for m = 1:Nc            % 趋化次数
  38.             disp(['趋化次数 m=', num2str(m), ' , 趋化次数 Nc=', num2str(Nc)])
  39.             for j=1:sizepop     % 种群
  40.       
  41.                 % 翻转
  42.                 delta = 2*rand(1,nvar)-0.5;
  43.                 pop(j,:) = pop(j,:) + C(j,:).*delta./(sqrt( delta*delta' ));
  44.                
  45.                 % 取值范围约束
  46.                 pop(j,:) = lb_ub(pop(j,:), popmin, popmax);
  47.                
  48.                 % 更新当前适应度值
  49.                 fitness(j) = fun(pop(j,:));

  50.                 % 适应度更新
  51.                 % 比较 个体间比较
  52.                 if fitness(j)<fitnesszbest
  53.                     fitnesszbest = fitness(j);
  54.                     zbest =  pop(j,:);
  55.                 end
  56.             end   % sizepop  种群数量
  57.             
  58.             % 记录最优适应度值
  59.             NcSizepop = NcSizepop+1;
  60.             fitness_iter(NcSizepop) = fitnesszbest;
  61.             
  62.             K_p(1,NcSizepop) = zbest(1);
  63.             K_i(1,NcSizepop) = zbest(2);
  64.             K_d(1,NcSizepop) = zbest(3);
  65.             
  66.         end       % Nc       趋化次数
  67.         
  68.         % 复制操作
  69.         [maxF,index] = sort(fitness,'descend');  % 降序排列
  70.         for Nre2 = 1:Sr   % 将最大适应度值的Sr个种群,进行更新
  71.             pop(index(Nre2),:) = popmin + (popmax-popmin).*rand(1,nvar);
  72.             fitness(index(Nre2)) = fun(pop(index(Nre2),:));
  73.             C(index(Nre2),:) = Cmin + (Cmax-Cmin).*rand(1,nvar);  % 步长
  74.             % 比较 个体间比较
  75.             if fitness(index(Nre2))<fitnesszbest
  76.                 fitnesszbest = fitness(index(Nre2));
  77.                 zbest =  pop(index(Nre2),:);
  78.             end
  79.         end
  80.     end   % Nre  复制操作
  81.    
  82.     for j=1:sizepop     % 种群
  83.         if Ped>rand
  84.             pop(j,:) = popmin + (popmax-popmin).*rand(1,nvar);
  85.             fitness(j) = fun(pop(j,:));
  86.             % 比较 个体间比较
  87.             if fitness(j)<fitnesszbest
  88.                 fitnesszbest = fitness(j);
  89.                 zbest =  pop(j,:);
  90.             end
  91.         end
  92.     end
  93.    
  94. end       % Ned   驱散(迁移)次数
  95. %% 显示最优解
  96. disp('最优解')
  97. disp(zbest)
  98. fprintf('\n')
  99. disp(['最优解Kp=', num2str(zbest(1))])
  100. disp(['最优解Ki=', num2str(zbest(2))])
  101. disp(['最优解Kd=', num2str(zbest(3))])

  102. figure('color',[1,1,1])
  103. plot(fitness_iter,'ro-','linewidth',2)
  104. % loglog(fitness_iter,'ro-','linewidth',2)
  105. axis tight
  106. grid on

  107. figure('color',[1,1,1])     % 绘制PID控制器参数变化曲线
  108. plot(K_p,'b-','LineWidth',2); hold on
  109. plot(K_i,'k-.','LineWidth',2)
  110. plot(K_d,'r--','LineWidth',2)
  111. title('Kp、Ki、Kd 优化曲线','fontsize',10);
  112. xlabel('迭代次数','fontsize',10);   ylabel('参数值','fontsize',10);
  113. set(gca,'Fontsize',10);
  114. legend('Kp','Ki','Kd',1);
  115. grid on;hold off;

  116. %% 响应曲线
  117. % [timef, r, yout]=PID_Fun_1_response(zbest);
  118. [timef, r, yout]=PID_Fun_2_response(zbest);
  119. figure('color',[1,1,1])
  120. plot(timef, r,'b-.','linewidth',2)
  121. hold on;
  122. plot(timef, yout,'r-o','linewidth',2)
  123. hold off;
  124. legend('输入信号','阶跃响应曲线')
  125. xlabel('time(s)');
  126. ylabel('Amp');
复制代码
PID整定对象:
  1. function BsJ=PID_Fun_2(Kpidi)
  2. ts=0.001;
  3. sys=tf([1.0],[1 1],'ioDelay',0.2);
  4. dsys=c2d(sys,ts,'z');
  5. [num,den]=tfdata(dsys,'v');
  6. u_1=0.0;
  7. y_1=0.0;
  8. x=[0,0,0]';
  9. B=0;
  10. error_1=0;
  11. tu=1;
  12. s=0;
  13. P=1000;
  14. for k=1:1:P
  15.     timef(k)=k*ts;
  16.     r(k)=1;
  17.     u(k)=Kpidi(1)*x(1)+Kpidi(2)*x(3)+Kpidi(3)*x(2);
  18. %     if u(k)>=10
  19. %         u(k)=10;
  20. %     end
  21. %     if u(k)<=-10
  22. %         u(k)=-10;
  23. %     end
  24.     yout(k)=-den(2)*y_1+num(2)*u_1;
  25.     error(k)=r(k)-yout(k);
  26.     %-------------------Return of PID parameters----------------
  27.     u_1=u(k);
  28.     y_1=yout(k);
  29.     x(1)=error(k);                % Calculating  P
  30.     x(2)=(error(k)-error_1)/ts;   % Calculating  D
  31.     x(3)=x(3)+error(k)*ts;        % Calculating  I
  32.     error_2=error_1;
  33.     error_1=error(k);
  34.     if s==0
  35.         if yout(k)>0.95 && yout(k)<1.05
  36.             tu=timef(k);
  37.             s=1;
  38.         end
  39.     end
  40. end
  41. % for i=1:1:P
  42. %     Ji(i)=0.999*abs(error(i))+0.01*u(i)^2*0.1;
  43. %     B=B+Ji(i);
  44. %     if i>1
  45. %         erry(i)=yout(i)-yout(i-1);
  46. %         if erry(i)<0
  47. %             B=B+100*abs(erry(i));
  48. %         end
  49. %     end
  50. % end
  51. % BsJ=B+0.2*tu*10;
  52. % BsJ = abs(error(P));
  53. BsJ = sum(abs(error));
复制代码








程序:算法QQ  3283892722
群智能算法视频:X元,链接http://halcom.cn/forum.php?mod=forumdisplay&fid=73
回复

使用道具 举报

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

本版积分规则


Python|Opencv|MATLAB|Halcom.cn  

GMT+8, 2019-6-17 08:45 , Processed in 0.129907 second(s), 25 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.

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