6-PSO_PID--粒子群算法--适应度函数有更改（采样时间）

Halcom

6-PSO_PID--粒子群算法--适应度函数有更改（采样时间）
百度网盘链接：https://pan.baidu.com/s/134vHa0M2gPMENML-t-AvnA 提取码：5hlg
具体链接在halcom.cn论坛，联系人QQ：3283892722
该论坛是一个学习交流平台，我会逐一的和大家分享学习。
欢迎大家录制视频，并提交给我，我来设置视频，你可在论坛进行打赏分享。
视频专用播放器：http://halcom.cn/forum.php?mod=viewthread&tid=258&extra=page%3D1
程序如下：

1. %% 清空环境
   
2. clc,clear,close all
   
3. warning off
   
4. format longG
   
5. addpath(genpath('./PID_funtion'))
   
6. %% PSO参数设置
   
7. w = 0.6;      % 惯性因子
   
8. c1 = 2;       % 加速常数
   
9. c2 = 2;       % 加速常数
   
10. Dim = 3;            % 维数，未知量个数nvar, N_VAR, N_PAR
    
11. SwarmSize = 30;     % 粒子群规模
    
12. MaxIter = 50;       % 最大迭代次数
    
13. Vmax = 1;
    
14. Vmin = -1;
    
15. Ub = [15 10 5];  % 上限
    
16. Lb = [0.5 0 0];  % 下限
    
17. 
    
18. %% 适应度函数
    
19. fun = @(x)PID_controller(x);
    
20. 
    
21. %% 粒子群初始化
    
22. Range = ones(SwarmSize,1)*(Ub-Lb);
    
23. pop = rand(SwarmSize,Dim).*Range + ones(SwarmSize,1)*Lb;    % 初始化粒子群
    
24. V = rand(SwarmSize,Dim)*(Vmax-Vmin) + Vmin;               % 初始化速度
    
25. fitness = zeros(SwarmSize,1);
    
26. for i=1:SwarmSize
    
27.     fitness(i,:) = fun(pop(i,:));                         % 粒子群的适应值
    
28. end
    
29. 
    
30. %% 个体极值和群体极值
    
31. [bestf, bestindex]=min(fitness);
    
32. zbest=pop(bestindex,:);   % 全局最佳
    
33. gbest=pop;                % 个体最佳
    
34. fitnessgbest=fitness;              % 个体最佳适应值
    
35. fitnesszbest=bestf;               % 全局最佳适应值
    
36. 
    
37. %% 迭代寻优
    
38. iter = 0;
    
39. y_fitness = zeros(1,MaxIter);   % 预先产生4个空矩阵
    
40. K_p = zeros(1,MaxIter);
    
41. K_i = zeros(1,MaxIter);
    
42. K_d = zeros(1,MaxIter);
    
43. while( (iter < MaxIter) )
    
44.     disp(['当前迭代次数：',num2str(iter)]) % 迭代次数
    
45.     for j=1:SwarmSize
    
46.         % 速度更新
    
47.         V(j,:) = w*V(j,:) + c1*rand*(gbest(j,:) - pop(j,:)) + c2*rand*(zbest - pop(j,:));
    
48.         for k=1:Dim
    
49.             if V(j,k)>Vmax,
    
50.                 V(j,k)=Vmax;
    
51.             end
    
52.             if V(j,k)<Vmin,
    
53.                 V(j,k)=Vmin;
    
54.             end
    
55.         end
    
56.         % 位置更新
    
57.         pop(j,:)=pop(j,:)+0.5*V(j,:);
    
58.         for k=1:Dim
    
59.             if pop(j,k)>Ub(k),
    
60.                 pop(j,k)=Ub(k);
    
61.             end
    
62.             if pop(j,k)<Lb(k),
    
63.                 pop(j,k)=Lb(k);
    
64.             end
    
65.         end
    
66.         % 适应值
    
67.         fitness(j,:) =fun(pop(j,:));
    
68.         % 个体最优更新
    
69.         if fitness(j) < fitnessgbest(j)
    
70.             gbest(j,:) = pop(j,:);
    
71.             fitnessgbest(j) = fitness(j);
    
72.         end
    
73.         % 群体最优更新
    
74.         if fitness(j) < fitnesszbest
    
75.             zbest = pop(j,:);
    
76.             fitnesszbest = fitness(j);
    
77.         end
    
78.     end
    
79.     iter = iter+1;                      % 迭代次数更新
    
80.     y_fitness(1,iter) = fitnesszbest;         % 为绘图做准备
    
81.     K_p(1,iter) = zbest(1);
    
82.     K_i(1,iter) = zbest(2);
    
83.     K_d(1,iter) = zbest(3);
    
84. end
    
85. %% 绘图
    
86. disp(['PSO最优解', num2str(zbest)])
    
87. disp(['PSO最优解对应的目标值', num2str(fitnesszbest)])
    
88. fprintf('\n')
    
89. 
    
90. figure(1)      % 绘制性能指标ITAE的变化曲线
    
91. plot(y_fitness,'LineWidth',3)
    
92. title('最优个体适应值','fontsize',10);
    
93. xlabel('迭代次数','fontsize',10);ylabel('适应值','fontsize',10);
    
94. set(gca,'Fontsize',10);
    
95. grid on
    
96. 
    
97. figure(2)      % 绘制PID控制器参数变化曲线
    
98. plot(K_p,'b-','LineWidth',2); hold on
    
99. plot(K_i,'k-.','LineWidth',2)
    
100. plot(K_d,'r--','LineWidth',2)
     
101. title('Kp、Ki、Kd 优化曲线','fontsize',10);
     
102. xlabel('迭代次数','fontsize',10);   ylabel('参数值','fontsize',10);
     
103. set(gca,'Fontsize',10);
     
104. legend('Kp','Ki','Kd',1);
     
105. grid on;hold off;
     
106. 
     
107. % PID响应
     
108. [BSJ, rint, yout]=PID_controller_response(zbest);
     
109. 
     
110. figure('color',[1,1,1])
     
111. plot(rint,'b-','linewidth',2);hold on;
     
112. plot(yout,'r.-','linewidth',2);
     
113. axis tight;grid on;
     
114. hold off;
     
115. title('阶跃响应曲线')
     
116. 
     
117. save PSO_result.mat
     
118. rmpath(genpath('./PID_funtion'))

复制代码

PID控制系统如下（细节、细节、细节）：

1. function BsJ=PID_controller(Kpidi)
   
2. % PID控制误差--适应度函数
   
3. % 输入：
   
4. %       Kpidi---包含Kp、Ki、Kd三个控制参数
   
5. % 输出：
   
6. %       BsJ---控制累计误差
   
7. %% 二阶控制系统
   
8. ts=1;                   % 采样周期
   
9. % sys=tf([-2,1],[0.2 1.2 1.0]);
   
10. % dsys=c2d(sys,ts,'z');
    
11. % [num,den]=tfdata(dsys,'v'); % z变换后的分子分母系数
    
12. % 减小计算时间
    
13. num = [0,-0.361018550145424,0.988879914151562];
    
14. den = [1,-0.374617388170528,0.00247875217666640];
    
15. 
    
16. %常规PID
    
17. u_11=0;u_21=0;u_31=0;u_41=0;u_51=0;
    
18. y_11=0;y_21=0;y_31=0;
    
19. error_21=0;
    
20. error_11=0;
    
21. ei1=0;
    
22. ei=0;
    
23. 
    
24. u_1=0;u_2=0;u_3=0;u_4=0;u_5=0;
    
25. y_1=0;y_2=0;y_3=0;
    
26. error_2=0;
    
27. error_1=0;
    
28. 
    
29.     kp1=Kpidi(1);
    
30.     ki1=Kpidi(2);
    
31.     kd1=Kpidi(3);
    
32.    
    
33. for k=1:1:200
    
34.     time(k)=k*ts;
    
35.    
    
36.     %  rin(k)=bb(k);
    
37.     rin(k)=1;
    
38.     yout1(k)=-den(2)*y_11-den(3)*y_21+num(2)*u_11+num(3)*u_21;
    
39.     error1(k)=(rin(k)-yout1(k))*0.05;
    
40.     ei1=ei1+error1(k)*ts;
    
41. 
    
42.     u1(k)= u_11+kp1*error1(k)+ kd1*(error1(k)-error_11)/ts+ki1*ei1;
    
43.     if u1(k)>=10       % Restricting the output of controller
    
44.         u1(k)=10;
    
45.     end
    
46.     if u1(k)<=-10
    
47.         u1(k)=-10;
    
48.     end
    
49.     u_51=u_41;u_41=u_31;u_31=u_21;
    
50.     u_21=u_11;u_11=u1(k);
    
51.     y_31=y_21;
    
52.     y_21=y_11;y_11=yout1(k);
    
53.     error_21=error_11;
    
54.     error_11=error1(k);
    
55. 
    
56. end
    
57. 
    
58. BsJ = sum(abs(error1));

复制代码