基于粒子群优化算法的PID控制器参数整定

创作时间:

作者:

@小白创作中心

基于粒子群优化算法的PID控制器参数整定

引用

CSDN

https://blog.csdn.net/qq_42249050/article/details/111239679

本文介绍了一种基于粒子群优化算法（PSO）整定PID控制器参数的方法。文章详细介绍了PSO算法的基本原理，并通过具体代码示例展示了如何使用PSO算法优化PID控制器参数。

0. 背景

PID控制器是工业控制中应用最广泛的控制器之一，其性能很大程度上取决于参数KP、KI和KD的设置。传统的PID参数整定方法往往依赖于经验或试错法，而基于优化算法的整定方法则能够更系统地搜索参数空间，找到最优解。本文将介绍如何使用粒子群优化算法（PSO）来整定PID控制器参数。

1. 算法介绍

1.1 标准的粒子群算法PSO

粒子群优化算法（PSO）是一种群体智能优化算法，最初由Kennedy和Eberhart于1995年提出。PSO算法模拟鸟群觅食行为，通过粒子之间的信息交流来寻找最优解。算法中，每个粒子代表一个潜在的解，粒子在解空间中不断更新位置和速度，以逼近全局最优解。

PSO算法的关键参数包括惯性权重w、认知学习因子c1和社交学习因子c2。粒子的位置更新公式如下：

其中，v表示粒子速度，x表示粒子位置，pbest表示粒子历史最优位置，gbest表示全局最优位置，r1和r2是[0,1]之间的随机数。

1.2 算法举例

为了验证PSO算法的有效性，我们使用Rastrigin函数作为测试函数。Rastrigin函数是一个多峰函数，常用于测试优化算法的性能。函数定义如下：

使用PSO算法优化Rastrigin函数的MATLAB代码如下：

%% 清空环境及变量
close all;
clear;
clc;
tic;
%% 目标函数
f=@(x)20+(x(:,1)).^2+(x(:,2)).^2+-10*(cos(2*pi*x(:,1))+cos(2*pi*x(:,2)));
[Pg,fmin]=PSO(f,2,100,200,40,-40,40,-40);
Pg,
fmin,
toc;

PSO函数的实现如下：

function [Pg,fmin]=PSO(f,dimension,n,m,xmax,xmin,vmax,vmin)
%全局粒子群算法,f为目标函数,dimension为维度,n为代数,m为种群规模
%位置限幅为[xmin,xmax]，速度限幅为[vmin,vmax]
    Savef=zeros(n+1,1);
    SaveData=zeros(m,dimension,10);%记录11代种群的位置
%% 固定惯性因子 m=100
    w=1;   
%     w=1.0-[1:1:m]*(1.0)/m; 
    c1=2;
    c2=2;   %速度惯性系数
    dt=0.3; %位移仿真间隔
    
    v=zeros(m,dimension);%初始速度为0
    
    X=(xmax-xmin)*rand(m,dimension)+xmin;%初始位置满足(-xmax,xmax)内均匀分布
    
    P=X;%P为每个粒子每代的最优位置
    
    last_f=f(X);
    
    [fmin,min_i]=min(last_f);%Pg为所有代中的最优位置 
    
    Pg=X(min_i,:);
    Savef(1)=fmin;
    N=0;
    
    for i=1:n
        
        v=w*v+c1*rand*(P-X)+c2*rand*(ones(m,1)*Pg-X);
        
        v=(v>vmax).*vmax+(v>=vmin & v<=vmax).*v+(v<vmin).*vmin;
        X=X+v*dt;
        X=(X>xmax).*xmax+(X>=xmin & X<=xmax).*X+(X<xmin).*xmin;
        
        new_f=f(X);%新的目标函数值
        
        update_j=find(new_f<last_f);
        P(update_j,:)=X(update_j,:);%修正每个粒子的历史最优值
        [new_fmin,min_i]=min(new_f);
        new_Pg=X(min_i,:);
        Pg=(new_fmin<fmin)*new_Pg+(new_fmin>=fmin)*Pg;
        last_f=new_f;%保存当前的函数值
        fmin=min(new_fmin,fmin);%更新函数最小值
         Savef(i)=fmin;
         
         if mod(i,floor(n/10))==0%每10代记录一次种群参数
             N=N+1;
            SaveData(:,:,N)=X;
         end
%         w=w-i/n*0.7*w;
    end
    
    for j=1:10
        figure(j)
        plot(SaveData(:,1,j),SaveData(:,2,j),'o');
        xlim([-xmax,xmax])
        ylim([-xmax,xmax])
        axis tight
    end
    
    figure
    plot(Savef','b-')
    disp(Pg)
    disp(fmin)
    
end

2. PID参数整定

2.1 M文件编写传递函数的PID参数整定

在PID参数整定中，我们通常关注响应曲线的超调量、上升时间、调节时间等性能指标。为了简化问题，我们选择调节时间和超调量作为评价函数的依据。评价函数定义如下：

使用PSO算法优化PID参数的MATLAB代码如下：

%---------------粒子群优化寻PID参数----------------%
close all;
clear;
clc;
tic;
dimension=3;
n=100;
m=50;%PID为3位参数，n是迭代次数，m为种群规模
% [Pg,fmin] = PSO(dimension,n,m);
%% PSO 优化函数---PSO算法优化函数
%全局粒子群算法,f为目标函数,dimension为维度,n为代数,m为种群规模
%     w=1;
    w=1.0-[1:1:n]*(1.0)/n; 
    
    c1=2;
    c2=2; %速度惯性系数
    
    sigma_data=zeros(1,n);
    ts_data=zeros(1,n);
    dt=0.3;%位移仿真间隔
    
    vmax=1;%速度限幅
    
    xmin=[0.2,0,0];
    xmax=[15,50,2];%位置限幅
    
    
    v=zeros(m,dimension);%初始速度为0
    X=(xmax-xmin).*rand(m,dimension)+xmin;%初始位置满足(xmin,xmax)内均匀分布
    P=X;%P为每个粒子每代的最优位置
    last_f=PID_access(X); %%目标适应度函数优化
    [fmin,min_i]=min(last_f);%Pg为所有代中的最优位置 
    Pg=X(min_i,:);
    N=0;
    
    for i=1:n
        v=w(i)*v+c1*rand*(P-X)+c2*rand*(ones(m,1)*Pg-X);
%         v=w(i)*v+c1*rand*(P-X)+c2*rand*(ones(m,1)*Pg-X);
        v=(v>vmax).*vmax+(v>=-vmax & v<=vmax).*v+(v<-vmax).*(-vmax);
        X=X+v*dt;
        X=(X>xmax).*xmax+(X>=xmin & X<=xmax).*X+(X<xmin).*(xmin);
        
        new_f=PID_access(X);%新的目标函数值
        update_j=find(new_f<last_f);
        P(update_j,:)=X(update_j,:);%修正每个粒子的历史最优值
        
        [new_fmin,min_i]=min(new_f);
        new_Pg=X(min_i,:);
        Pg=(new_fmin<fmin)*new_Pg+(new_fmin>=fmin)*Pg; %%三个PID优化参数的优化区间设计
        
        last_f=new_f;%保存当前的函数值
        fmin=min(new_fmin,fmin);%更新函数最小值
        Savef(i)=fmin;
        
        [~,~,ts,sigma]= PID_sim(Pg(1),Pg(2),Pg(3),true); %%传递函数设计
        
        sigma_data(1,i)=sigma;
        ts_data(1,i)=ts;
        hold on
    end
    legend('1','2','3','4','5','6','7','8','9','10','11','12','13','14','15')
    title('响应曲线随迭代次数的变化')
    set(gca,'XLim',[0 1.0]);%X轴的数据显示范围
    set(gca,'YLim',[0 1.2]);%Y轴的数据显示范围
    
    figure
    plot(1:1:n,Savef');
    title('算法优化最小适应度函数变化');
    
    figure
    plot(ts_data)
    title('调节时间随迭代次数的变化')
    
    figure
    plot(sigma_data)
    title('超调量随迭代次数的变化')
    
toc;
%% 根据参数得出评价度的函数，程序中为评价度越低越好
function y=PID_access(para)
%PID调节性能的指标参数
kp=para(:,1);
ki=para(:,2);
kd=para(:,3);
[~,~,ts,sigma]=PID_sim(kp,ki,kd,false);
y=log(ts/5e-2+1)+log(sigma/1e-2+1);
end
%% PID并行仿真引擎
function [f_infty,tp,ts,sigma]=PID_sim(kp,ki,kd,debug)
    %kp,ki,kd为PID参数,T0为采样时间,total_t为仿真时间
    
    Tt=5e-4;%仿真采样周期
%     Tt=1e-3;
    T0=1e-2;%控制采样周期
    
    Tf=1e-3;%微分环节的滤波器系数
    alp=Tf/(Tt+Tf);
    total_t=1;
    N=floor(total_t/Tt);%样本总数
    M=numel(kp);% 返回仿真序列数
    
    kp=reshape(kp,M,1);
    ki=reshape(ki,M,1);
    kd=reshape(kd,M,1);
    
%     sys=tf(400,[1,1,0]);  
    sys=tf(0.998,[0.021,1]);   %建立被控对象传递函数
    dsys=c2d(sys,Tt,'z');      %离散化
    [num,den]=tfdata(dsys,'v');   %
    
    e_1=zeros(M,1);      %前一时刻的偏差      
    Ee=zeros(M,1);       %累积偏差
    u_1=zeros(M,1);    %前一时刻的控制量
    y_1=zeros(M,1);       %前一时刻的输出
    ud_1=zeros(M,1);     %前一时刻的微分输出
    
    %预分配内存
    r=1;
    y=zeros(M,N);
    u=zeros(M,N);
    
    for k=1:1:N
        y(:,k)=-1*den(2).*y_1+num(2)*u_1+num(1).*u(:,k);%系统响应输出序列
        e0=r-y(:,k);   %误差信号
        ud=alp.*ud_1+(1-alp)/T0.*kd.*(e0-e_1);
        u(:,k)=kp.*e0+T0*ki.*Ee+ud; %系统PID控制器输出序列
        Ee=Ee+e0;    %误差的累加和
        u_1=u(:,k);    	%前一个的控制器输出值
        y_1=y(:,k);    	%前一时刻的系统响应输出值
        e_1=e0;		%前一个误差信号的值
        ud_1=ud;
    end
    
    if debug %非调试模式下不显示也不打印图像
        plot(linspace(0,total_t,N),y)
    end
    [f_infty,tp,ts,sigma]=parameter_cal(y,Tt,2e-2,debug);%取得阶跃响应指标
    
end
%% 计算仿真曲线参数
function [f_infty,tp,ts,sigma]=parameter_cal(y,T0,delt_err,debug)
%y是系统输出序列
%T0是采样时间，可以结合时间序列点序号算出实际时间
%delt_err是调节时间的误差精度
    M=size(y,1);
    N=size(y,2);  %M为运算维度，N为时间序列长度
    f_infty=y(:,N);%稳态值序列
    err=y-f_infty*ones(1,N);%通过稳态值计算误差序列
    ferr=fliplr(abs(err));%倒序并取绝对值
    [~,ts_i]=max(ferr>delt_err*f_infty,[],2);
    ts_i=N*ones(M,1)-ts_i;
    ts=ts_i*T0;%调节时间
    [fp,tp]=max(y,[],2);%峰值和函数峰值
    tp=tp*T0;
    tp(abs(fp-f_infty)<1e-5)=NaN; %过阻尼无超调，没有峰值时间
    sigma=(fp-f_infty)./f_infty;
    
    if debug && M==1 %非调试模式下不显示
        disp(['系统稳态值为',num2str(f_infty)])
        disp(['系统超调量为',num2str(sigma*100),'%'])
            if isnan(tp)
                disp('系统不存在峰值时间')
                else
                disp(['系统峰值时间为',num2str(tp),'s'])
            end
        disp(['系统的调节时间为',num2str(ts),'s'])
    end
    
end

2.2 总结

可以对某个系统运用辨识算法得到精确的传递函数，基于此运用上述算法进行优化PID控制器参数，同理可推广至其他控制器的参数优化。
对于传递函数的编写，除运用m文件进行编写外，还可以通过Simulink仿真模型进行搭建。如此，也可以通过算法与Simulink模型进行在线的参数整定。后面将重点讲述下数据交互的方式与问题。