1 理论基础
PID控制器应用广泛,其一般形式为
可见,PID控制器的性能取决于Kp、Ki、Kd这3个参数是否合理,因此,优化PID控制器参数具有重要意义。目前,PID控制器参数主要是人工调整,这种方法不仅费时,而且不能保证获得最佳的性能。PSO已经广泛应用于函数优化、神经网络训练、模式分类、模糊系统控制以及其他应用领域,使用PSO进行PID控制器参数的优化设计。
2 背景
2.1 问题描述
PID控制器的系统结构图如图14-1所示。
在Simulink环境下建立的模型如图14-2所示。
图14-2中,微分环节由一个一阶环节近似,输出端口1即为式(14-2)所示的ITAE指标,通过将时间及误差绝对值的乘积进行积分后得到。
2.2 解题思路及步骤
1.优化设计过程
利用粒子群算法对PID控制器的参数进行优化设计,其过程如图14-3所示。
图14-3中,粒子群算法与Simulink模型之间连接的桥梁是粒子(即PID控制器参数)和该粒子对应的适应值(即控制系统的性能指标)。优化过程如下:PSO产生粒子群(可以是初始化粒子群,也可以是更新后的粒子群),将该粒子群中的粒子依次赋值给PID控制器的参数Kp、Ki、Kd,然后运行控制系统的Simulink模型,得到该组参数对应的性能指标,该性能指标传递到PSO中作为该粒子的适应值,最后判断是否可以退出算法。
2.粒子群算法实现
粒子群算法的基本原理在前面已经述及,进一步地,粒子在搜索空间中的速度和位置根据以下公式确定:
PSO的流程如下:
(1)初始化粒子群,随机产生所有粒子的位置和速度,并确定粒子的Pt和Gt。
(2)对每个粒子,将其适应值与该粒子所经历过的最优位置P,的适应值进行比较,如较好,则将其作为当前的Pt。
(3)对每个粒子,将其适应值与整个粒子群所经历过的最优位置Gt的适应值进行比较,如较好,则将其作为当前的Gt。
(4)按式(14-4)和式(14-5)更新粒子的速度和位置。
(5)如果没有满足终止条件(通常为预设的最大迭代次数和适应值下限值),则返回步骤(2);否则,退出算法,得到最优解。
3 MATLAB程序实现
3.1 Simulink部分的程序实现
图14-3所示的PSO优化PID过程示意图,其右侧的Simulink部分的程序实现如下:
function z=PSO_PID(x) assignin('base','Kp',x(1)); assignin('base','Ki',x(2)); assignin('base','Kd',x(3)); [t_time,x_state,y_out]=sim('PID_Model',[0,20]); z=y_out(end,1);
其中,x为PSO中传递过来的粒子,首先,调用函数assignin将x(1)、x(2)、x(3)的值赋值给Workspace中的Kp、Ki、Kd,该语句实现了图14-3中从PSO部分到Simulink部分的参数传递;然后,调用函数sim对图14-2所示的模型进行仿真,其中,PID_Model为Simulink模型的文件名,[0,20]为仿真时间,返回的y_out即为输出端子1的值;最后,将性能指标ITAE赋值给z,以实现图14-3中从Simulink部分到PSO部分的参数传递。
3.2 PSO部分的程序实现
设置PSO的参数为:惯性因子w=0.6,加速常数c?=c?=2,维数为3(有3个待优化参数),粒子群规模为100,待优化函数为14.3.1节中的函数PSO_PID,最大迭代次数为100,最小适应值为0.1,速度范围为[-1,1],3个待优化参数范围均为[0,300]。代码如下:
%% 清空环境 clear clc %% 参数设置 w = 0.6; % 惯性因子 c1 = 2; % 加速常数 c2 = 2; % 加速常数 Dim = 3; % 维数 SwarmSize = 100; % 粒子群规模 ObjFun = @PSO_PID; % 待优化函数句柄 MaxIter = 100; % 最大迭代次数 MinFit = 0.1; % 最小适应值 Vmax = 1; Vmin = -1; Ub = [300 300 300]; Lb = [0 0 0]; %% 粒子群初始化 Range = ones(SwarmSize,1)*(Ub-Lb); Swarm = rand(SwarmSize,Dim).*Range + ones(SwarmSize,1)*Lb % 初始化粒子群 VStep = rand(SwarmSize,Dim)*(Vmax-Vmin) + Vmin % 初始化速度 fSwarm = zeros(SwarmSize,1); for i=1:SwarmSize fSwarm(i,:) = feval(ObjFun,Swarm(i,:)); % 粒子群的适应值 end %% 个体极值和群体极值 [bestf bestindex]=min(fSwarm); zbest=Swarm(bestindex,:); % 全局最佳 gbest=Swarm; % 个体最佳 fgbest=fSwarm; % 个体最佳适应值 fzbest=bestf; % 全局最佳适应值 %% 迭代寻优 iter = 0; y_fitness = zeros(1,MaxIter); % 预先产生4个空矩阵 K_p = zeros(1,MaxIter); K_i = zeros(1,MaxIter); K_d = zeros(1,MaxIter); while( (iter < MaxIter) && (fzbest > MinFit) ) for j=1:SwarmSize % 速度更新 VStep(j,:) = w*VStep(j,:) + c1*rand*(gbest(j,:) - Swarm(j,:)) + c2*rand*(zbest - Swarm(j,:)); if VStep(j,:)>Vmax, VStep(j,:)=Vmax; end if VStep(j,:)<Vmin, VStep(j,:)=Vmin; end % 位置更新 Swarm(j,:)=Swarm(j,:)+VStep(j,:); for k=1:Dim if Swarm(j,k)>Ub(k), Swarm(j,k)=Ub(k); end if Swarm(j,k)<Lb(k), Swarm(j,k)=Lb(k); end end % 适应值 fSwarm(j,:) = feval(ObjFun,Swarm(j,:)); % 个体最优更新 if fSwarm(j) < fgbest(j) gbest(j,:) = Swarm(j,:); fgbest(j) = fSwarm(j); end % 群体最优更新 if fSwarm(j) < fzbest zbest = Swarm(j,:); fzbest = fSwarm(j); end end iter = iter+1; % 迭代次数更新 y_fitness(1,iter) = fzbest; % 为绘图做准备 K_p(1,iter) = zbest(1); K_i(1,iter) = zbest(2); K_d(1,iter) = zbest(3); end %% 绘图输出 figure(1) % 绘制性能指标ITAE的变化曲线 plot(y_fitness,'LineWidth',2) title('最优个体适应值','fontsize',18); xlabel('迭代次数','fontsize',18);ylabel('适应值','fontsize',18); set(gca,'Fontsize',18); figure(2) % 绘制PID控制器参数变化曲线 plot(K_p) hold on plot(K_i,'k','LineWidth',3) plot(K_d,'--r') title('Kp、Ki、Kd 优化曲线','fontsize',18); xlabel('迭代次数','fontsize',18);ylabel('参数值','fontsize',18); set(gca,'Fontsize',18); legend('Kp','Ki','Kd',1);
其中,MaxIter和MinFit即终止条件;Vmax和Vmin分别为速度的上限和下限;Ub(i)和Lb(i)分别为第i个待优化参数的上限和下限。粒子群的初始化采用与遗传算法相似的方法(用函数rand且保证粒子在上下限范围内),迭代过程采用while进行大循环,速度更新和位置更新按照式 (14-4)和式(14-5)进行(且在迭代 更新过程中,若超出了限值,则将其 设为限制)。
3.3 结果分析
运行以上代码,得到优化过程如图14-4和图14-5所示,前者为PID控制器3个参数Ki、Kp、Kd的变化曲线,后者为性能指标ITAE的变化曲线。得到的最优控制器参数及性能指标为
Kp= 33.6469,Ki=0.1662, Kd=38.7990, ITAE =1.0580
将以上参数代回图14-2所示的模型,得到的单位阶跃响应曲线如图14-6所示。由图14-5可知,算法优化过程中,性能指标ITAE不断减小,PSO不断寻找更优的参数。
关注公众号:“电击小子程高兴的MATLAB小屋”获取完整资源