蚁群算法ACO的PID参数整定

Halcom

蚁群算法ACO的PID参数整定：
代码参考链接：链接：https://pan.baidu.com/s/15nqkVj9-WuwI-2HQ1gB7pg 提取码：howq

1. clc,clear,close all
   
2. warning off
   
3. feature jit off
   
4. load('U.mat')
   
5. load('Y.mat')
   
6. % 对Y进行抽样
   
7. Y1 = Y(1:1000:end);
   
8. % ACO 参数
   
9. maxiter = 100;  % 迭代次数
   
10. sizepop = 10;  % 种群数量
    
11. popmin = [-1,-1,-1, 0.1, 0.1, 0.1 ]; % x_min
    
12. popmax = [50,50,50, 15.0, 15.0, 15.0 ]; % x_max
    
13. Rou = 0.8;     % 信息素增量强度
    
14. P0 = 0.1;      % 转移概率   
    
15. nvar = 6;      % nvar个自变量
    
16. % 初始化种群
    
17. for i=1:sizepop
    
18.     pop(i,:) = popmin + (popmax-popmin).*rand(1,nvar)*0.5;
    
19.     fitness(i) = BJ_Fun2( pop(i,:), Y1 );
    
20. end
    
21. % 记录一组最优值
    
22. [bestfitness,bestindex]=min(fitness);
    
23. zbest=pop(bestindex,:);   % 全局最佳
    
24. gbest=pop;                % 个体最佳
    
25. fitnessgbest=fitness;     % 个体最佳适应度值
    
26. fitnesszbest=bestfitness; % 全局最佳适应度值
    
27. %% 迭代寻优
    
28. for i=1:maxiter
    
29.     lamda = 1/i;       % 信息素挥发因子
    
30.     [bestfit,flag] =min(fitness);
    
31.     for j=1:sizepop
    
32.         Pt(j) = (fitness(j)-bestfit)./bestfit;
    
33.     end
    
34.     for j=1:sizepop
    
35.         % 转移概率值
    
36.         if Pt(j)<P0
    
37.             pop(j,:) = pop(j,:) + (2*rand(1, nvar)-1)*lamda/2;
    
38.         else
    
39.             pop(j,:) = pop(j,:) + (popmax-popmin).*(rand(1, nvar)-0.5);
    
40.         end   
    
41.         
    
42.         % 越界处理
    
43.         for k=1:nvar
    
44.             if pop(j,k)>popmax(k)
    
45.                 pop(j,k)=popmax(k);
    
46.             end
    
47.             if pop(j,k)<popmin(k)
    
48.                 pop(j,k)=popmin(k);
    
49.             end
    
50.         end
    
51.         
    
52.         % 适应度值
    
53.         fitness1 = BJ_Fun2(pop(j,:), Y1);
    
54.         fitness2 = BJ_Fun2(gbest(j,:), Y1);
    
55.         if fitness1<fitness2
    
56.             gbest(j,:) = pop(j,:);
    
57.         else
    
58.             pop(j,:) = gbest(j,:);
    
59.         end
    
60.         
    
61.         fitness(j) = (1-Rou)*fitness(j) + BJ_Fun2(pop(j,:), Y1);
    
62.         
    
63.         if fitness(j) < fitnesszbest
    
64.             fitnesszbest = fitness(j);
    
65.             fitnessgbest(j) = BJ_Fun2(pop(j,:), Y1);
    
66.         end
    
67.     end
    
68.     [fitness_iter(i),index]= min(fitnessgbest);
    
69. end
    
70. zbest = pop(index,:);  % 最优个体
    
71. save zbest.mat zbest
    
72. %% 显示
    
73. disp('最优解')
    
74. disp(['最优解： ', num2str( zbest )])
    
75. fprintf('\n')
    
76. disp(['最优解对应的目标函数值： ', num2str( bestfitness )])
    
77. disp('传递函数');
    
78. sys = tf( zbest(4), [1, zbest(5), zbest(4)] )
    
79. %% 画图显示
    
80. figure('color',[1,1,1])
    
81. plot(fitness_iter,'ro-','linewidth',2)
    
82. grid on
    
83. xlabel('迭代次数');
    
84. ylabel('适应度函数值')
    
85. %% 响应曲线
    
86. Ygap = 100;
    
87. [rint_BFO, yout_BFO] = BJ_Fun_response2(zbest, Y(1:Ygap:end));
    
88. figure;hold on;
    
89. plot(rint_BFO.*max(Y(1:Ygap:end)),'r-','linewidth',2);
    
90. plot(yout_BFO.*max(Y(1:Ygap:end)),'b.-','linewidth',2);
    
91. legend('阶跃控制信号','ACO算法优化的信号逼近')
    
92. axis tight;
    
93. hold off

复制代码

适应度函数如下：

1. function fitness = BJ_Fun2(xx, Y1)
   
2. ts=0.001;
   
3. KP=xx(1); KI=xx(2); KD=xx(3);
   
4. aa=xx(4); bb=xx(5); cc=xx(6);
   
5. % sys = tf( [1.6], [1,1.5,1.6] );
   
6. sys = tf( aa, [1, bb, aa] );
   
7. dsys = c2d(sys, ts, 'z');
   
8. [num,den] = tfdata( dsys, 'v' );
   
9. u_1 = 0; u_2=0;
   
10. y_1 = 0; y_2=0;
    
11. error_1 = 0;
    
12. x = [0,0,0]';
    
13. for k=1:length(Y1)
    
14.     timef(k) = k*ts;
    
15.     r(k) = Y1(k)./(max(Y1)+eps);
    
16.     u(k) = KP*x(1)+KI*x(2)+KD*x(3);
    
17.     yout(k) = -den(2)*y_1 -den(3)*y_2 + num(2)*u_1 + num(3)*u_2;
    
18.     error(k) = r(k) - yout(k);
    
19.     % PID参数返回值
    
20.     u_2 = u_1;
    
21.     u_1 = u(k);
    
22.     y_2 = y_1;
    
23.     y_1 = yout(k);
    
24.    
    
25.     x(1) = error(k);              % Kp
    
26.     x(2) = x(2) + error(k)*ts;    % KI
    
27.     x(3) = (error(k)-error_1)/ts; % KD
    
28.    
    
29.     error_1 = error(k);
    
30. end
    
31. fitness = sum( abs(error) ) + 300*min( abs( error ) );

复制代码

响应度曲线：

1. function [r, yout] = BJ_Fun_response2(xx, Y1)
   
2. KP=xx(1); KI=xx(2); KD=xx(3);
   
3. aa=xx(4); bb=xx(5); cc=xx(6);
   
4. ts=0.001;
   
5. % sys = tf( [1.6], [1,1.5,1.6] );
   
6. sys = tf( aa, [1, bb, aa] );
   
7. dsys = c2d(sys, ts, 'z');
   
8. [num,den] = tfdata( dsys, 'v' );
   
9. u_1 = 0; u_2=0;
   
10. y_1 = 0; y_2=0;
    
11. error_1 = 0;
    
12. x = [0,0,0]';
    
13. for k=1:length(Y1)
    
14.     timef(k) = k*ts;
    
15.     r(k) = Y1(k)./(max(Y1)+eps);
    
16.     u(k) = KP*x(1)+KI*x(2)+KD*x(3);
    
17.     yout(k) = -den(2)*y_1 -den(3)*y_2 + num(2)*u_1 + num(3)*u_2;
    
18.     error(k) = r(k) - yout(k);
    
19.     % PID参数返回值
    
20.     u_2 = u_1;
    
21.     u_1 = u(k);
    
22.     y_2 = y_1;
    
23.     y_1 = yout(k);
    
24.    
    
25.     x(1) = error(k);              % Kp
    
26.     x(2) = x(2) + error(k)*ts;    % KI
    
27.     x(3) = (error(k)-error_1)/ts; % KD
    
28.    
    
29.     error_1 = error(k);
    
30. end
    
31. % fitness = sum( abs(error) ) + 300*min( abs( error ) );

复制代码