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

Halcom

基于细菌觅食算法的PID整定：
百度网盘链接：https://pan.baidu.com/s/1khndoz4gj7EI2gabUldaeQ  提取码：v8l0
具体链接在halcom.cn论坛，联系人QQ：3283892722
该论坛是一个学习交流平台，我会逐一的和大家分享学习。
欢迎大家录制视频，并提交给我，我来设置视频，你可在论坛进行打赏分享。
视频专用播放器：http://halcom.cn/forum.php?mod=viewthread&tid=258&extra=page%3D1

基于细菌觅食算法的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.                 % 比较 个体间比较
    
53.                 if fitness(j)<fitnesszbest
    
54.                     fitnesszbest = fitness(j);
    
55.                     zbest =  pop(j,:);
    
56.                 end
    
57.             end   % sizepop  种群数量
    
58.             
    
59.             % 记录最优适应度值
    
60.             NcSizepop = NcSizepop+1;
    
61.             fitness_iter(NcSizepop) = fitnesszbest;
    
62.             
    
63.             K_p(1,NcSizepop) = zbest(1);
    
64.             K_i(1,NcSizepop) = zbest(2);
    
65.             K_d(1,NcSizepop) = zbest(3);
    
66.             
    
67.         end       % Nc       趋化次数
    
68.         
    
69.         % 复制操作
    
70.         [maxF,index] = sort(fitness,'descend');  % 降序排列
    
71.         for Nre2 = 1:Sr   % 将最大适应度值的Sr个种群，进行更新
    
72.             pop(index(Nre2),:) = popmin + (popmax-popmin).*rand(1,nvar);
    
73.             fitness(index(Nre2)) = fun(pop(index(Nre2),:));
    
74.             C(index(Nre2),:) = Cmin + (Cmax-Cmin).*rand(1,nvar);  % 步长
    
75.             % 比较 个体间比较
    
76.             if fitness(index(Nre2))<fitnesszbest
    
77.                 fitnesszbest = fitness(index(Nre2));
    
78.                 zbest =  pop(index(Nre2),:);
    
79.             end
    
80.         end
    
81.     end   % Nre  复制操作
    
82.    
    
83.     for j=1:sizepop     % 种群
    
84.         if Ped>rand
    
85.             pop(j,:) = popmin + (popmax-popmin).*rand(1,nvar);
    
86.             fitness(j) = fun(pop(j,:));
    
87.             % 比较 个体间比较
    
88.             if fitness(j)<fitnesszbest
    
89.                 fitnesszbest = fitness(j);
    
90.                 zbest =  pop(j,:);
    
91.             end
    
92.         end
    
93.     end
    
94.    
    
95. end       % Ned   驱散(迁移)次数
    
96. %% 显示最优解
    
97. disp('最优解')
    
98. disp(zbest)
    
99. fprintf('\n')
    
100. disp(['最优解Kp=', num2str(zbest(1))])
     
101. disp(['最优解Ki=', num2str(zbest(2))])
     
102. disp(['最优解Kd=', num2str(zbest(3))])
     
103. 
     
104. figure('color',[1,1,1])
     
105. plot(fitness_iter,'ro-','linewidth',2)
     
106. % loglog(fitness_iter,'ro-','linewidth',2)
     
107. axis tight
     
108. grid on
     
109. 
     
110. figure('color',[1,1,1])     % 绘制PID控制器参数变化曲线
     
111. plot(K_p,'b-','LineWidth',2); hold on
     
112. plot(K_i,'k-.','LineWidth',2)
     
113. plot(K_d,'r--','LineWidth',2)
     
114. title('Kp、Ki、Kd 优化曲线','fontsize',10);
     
115. xlabel('迭代次数','fontsize',10);   ylabel('参数值','fontsize',10);
     
116. set(gca,'Fontsize',10);
     
117. legend('Kp','Ki','Kd',1);
     
118. grid on;hold off;
     
119. 
     
120. %% 响应曲线
     
121. % [timef, r, yout]=PID_Fun_1_response(zbest);
     
122. [timef, r, yout]=PID_Fun_2_response(zbest);
     
123. figure('color',[1,1,1])
     
124. plot(timef, r,'b-.','linewidth',2)
     
125. hold on;
     
126. plot(timef, yout,'r-o','linewidth',2)
     
127. hold off;
     
128. legend('输入信号','阶跃响应曲线')
     
129. xlabel('time(s)');
     
130. 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));

复制代码

Halcom

4个对应的整定传递函数：

1. >> ts=0.001;
   
2. sys=tf([1.0],[1 1])
   
3. 
   
4. sys =
   
5. 
   
6.     1
   
7.   -----
   
8.   s + 1
   
9. 
   
10. Continuous-time transfer function.
    
11. 
    
12. >> sys=tf([1.0],[1 1],'ioDelay',0.2)
    
13. 
    
14. sys =
    
15. 
    
16.                   1
    
17.   exp(-0.2*s) * -----
    
18.                 s + 1
    
19. 
    
20. Continuous-time transfer function.
    
21. 
    
22. >> sys=tf([1.6],[1 2.584 1.6])
    
23. 
    
24. sys =
    
25. 
    
26.           1.6
    
27.   -------------------
    
28.   s^2 + 2.584 s + 1.6
    
29. 
    
30. Continuous-time transfer function.
    
31. 
    
32. >> sys=tf([1.6],[1 2.584 1.6],'ioDelay',0.1)
    
33. 
    
34. sys =
    
35. 
    
36.                         1.6
    
37.   exp(-0.1*s) * -------------------
    
38.                 s^2 + 2.584 s + 1.6
    
39. 
    
40. Continuous-time transfer function.

复制代码

wzshtuhao

不错，很有用！

pineapple

PID_Fun_2_response(zbest)这个函数代码里没有，不是完整版吗

Halcom

pineapple 发表于 2023-5-26 09:52
PID_Fun_2_response(zbest)这个函数代码里没有，不是完整版吗

function BsJ=PID_Fun_2(Kpidi)
这个函数不是提供了吗？
你修改以下输出即可