【Simulink】PSO优化算法整定PID控制器参数（一）一一一高阶不稳定系统

最新推荐文章于 2024-08-16 14:31:50 发布

置顶昔时扬尘处

最新推荐文章于 2024-08-16 14:31:50 发布

阅读量1.2w

点赞数 23

分类专栏：基于智能算法的PID参数整定研究

本文链接：https://blog.csdn.net/qq_42249050/article/details/111239679

版权

粒子群优化 PID控制器参数整定 MATLAB 传递函数

关键词由CSDN通过智能技术生成

基于智能算法的PID参数整定研究专栏收录该内容

23 篇文章 150 订阅

订阅专栏

0背景

写在前面：
1.本代码基于MATLAB2019a版本，如果低版本或者不同版本可能会报错，mdl文件或者slx文件打开可能会失败；
2.附上代码并详细介绍；
3.slx文件下载链接：见评论区

现在给大家介绍几种算法整定PID控制器参数的方法，讲到底还是基于数据驱动的，运用算法在有限解集空间内按一定的评价函数去搜索所想要的最优解。优化算法常见的有模拟退火、遗传算法、粒子群算法、蚁群算法和多目标优化算法等。本博客主要以粒子群算法PSO为主要讲解目标，所参考文献及资料见文末，代码可运行成功，复现论文中的程序及代码。

我们都知道PID控制器中的偏差和相应的系数KP、KI、KD组成了它的3个重要的环节。这3个系数的不同组合会影响控制系统响应的性能指标的值，所以KP、KI、KD这3个系数的优化是PID控制的关键。如图1所示，为PID控制器原理图，其中被控对象可大可小，也可以是精确的被控传递函数，也可以是无法表示的Simulink混杂系统(离散和连续的混合)，当然也可以是整个半实物仿真在线运行系统等。本博客主要是以传递函数为主的被控对象，运用PSO优化算法来进行整定与优化PID的三个控制参数。

在这里插入图片描述
图1 PID控制器原理图

其中，由图可推导出相应的控制器表达式，如下所示。
在这里插入图片描述

由上可看出系统响应的稳态误差ess、上升时间tr、超调量等性能指标的值，不同的KP、KI、KD组合会得到不同的性能指标。稳态误差越小、上升时间越短、超调量越小的系统响应越好。PID控制器的参数优化就是找到那组使得这些性能指标更好的KP、KI、KD系数的值。从而让系统的响应结果满足一些实际的需要。在这里如果不太懂PID控制器的同学，可以参考我之前写的博客基于遗传算法的PID参数整定研究（一）。其中就讲到了PID控制参数与控制性能之间的影响关系，以及所常见的一些参数整定方法，在此不再赘述。

1算法介绍

1.1标准的粒子群算法PSO

1995年提出的粒子群算法（PSO）是一种群体迭代算法，**所需要考虑的参数较少，惯性因子和学习因子，相比于遗传算法少了交叉和变异，仅存在选择。**作一个形象的比喻那就是，鸟群中的每一只鸟通过与群体中的其它同类相互交流从而交换信息来更新自己的飞行速度和所在位置，逐渐逼近食物的附近从而找到食物所在的位置，研究者们将鸟群中的个体与粒子相对应，将鸟群群体与粒子群相对应，将鸟群寻找食物过程中的规律对应于迭代进化式（2）、式（3），将鸟群寻找食物过程中的每一个粒子对应于待优化问题的可行解，初始化的粒子群是在规定可行解所在的范围内随机产生的，通过式（2），式（3）可以看出在其迭代更新机制是由以下三个部分所组成的：
在这里插入图片描述

式（2）中第一部分的系数w叫作惯性权重值，第二部分的系数由两个值组成第一个代表的是粒子个体所经历的最好位置对新一代种群位置影响的程度，我们把这个表示程度的系数叫作学习因子用c1表示，另一个是随机均匀分布在[0 1]范围内的随机变量用n表示，第三部分的系数组成同样由两个部分组成，一个代表的是粒子群所经历的最好位置对新一代种群位置影响的程度，我们把它也叫作学习因子用c2表示，另一个是随机均匀分布在[0 1]范围内的随机变量用n表示。

在PSO中，学习因子c1、c2反映了粒子群之间的信息交流；与此同时，c1、c2也决定了粒子自身经验信息及其它粒子经验信息对粒子运动轨迹的影响；c1取值较大，则粒子会在局部过多徘徊；c2值较大会使粒子早熟，陷入局部最优。

1.2算法举例

在本算法举例中，以确定的目标函数中去验证标准的算法运行的可行性。其中经常使用一个函数来测试遗传算法，这个函数就是Rastrigin函数。函数如下所示：
在这里插入图片描述

如何利用遗传算法来寻找Rastrigin函数的最小值？

以遗传算法GA工具箱来优化所得结果如下：

命令行窗口输入optimtool(‘ga’)打开优化工具的图形界面

optimtool('ga')

在这里插入图片描述
在 Fitness function 函数区域输入函数句柄@rastriginsfcn
在 Number of variables输入独立变量的个数2

在这里插入图片描述
当算法运行时，在Current iteration 区域会显示当前迭代的次数，算法运行的最后的结果是：

Objective function value: 0.3226433098693917

算法停止的原因是：average change in the fitness value less than options.TolFun.（适应度函数值的平均变化小于TolFun这个参数）。

在命令行中，输入下面的命令：

[x fval exitflag] = ga(@rastriginsfcn, 2)

算法运行的结果如下所示：
在这里插入图片描述

以标准粒子群PSO来优化所得结果如下：
在这里插入图片描述

主函数

%% 清空环境及变量
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

在这里可以进行惯性因子w的改进设计，并进行比较来进行验证上述的结论。

在这里插入图片描述

在这里插入图片描述
由上可知，固定惯性因子w=1.0下，各个粒子的收敛并没有最终收敛到最后一个点，在附近还存在着零散的点，因此将惯性因子w进行改进，假定按线性递减的方式去进行优化上述算法，如下：

 w=1.0-[1:1:m]*(1.0)/m;

以改进的粒子群PSO来优化所得结果如下：

在这里插入图片描述

优化的结果并没有上述的好，Rastrigin函数极小值是位于原点附近。

在这里插入图片描述

2PID参数整定

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

在自适应度函数的作用下，PSO算法会选择到我们满意的参数。在调整PID参数时，我们通常会关注响应曲线的超调、上升时间、调节时间、峰值时间等等。利用这些参数组合为评价函数，让评价函数越小，仿真得到的参数越好，以此作为解的好坏的对比依据（相当于鸟离食物的距离）为了简便起见，我们只选取了调节时间和超调量σ作为我们的比较依据，我们将评价函数设定为：
在这里插入图片描述
代码如下

%---------------粒子群优化寻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

以标准的粒子群PSO来优化PID控制参数所得结果如下：
在这里插入图片描述

在这里插入图片描述

2.2总结

1.可以对某个系统运用辨识算法得到精确的传递函数，基于此运用上述算法进行优化PID控制器参数，同理可推广至其他控制器的参数优化。
2.对于传递函数的编写，除运用m文件进行编写外，还可以通过Simulink仿真模型进行搭建。如此，也可以通过算法与Simulink模型进行在线的参数整定。后面将重点讲述下数据交互的方式与问题。
3.时间过得很快，学习不应止步。临近毕业，有一小段时间没更新博客在此博主也很遗憾，将继续更新下去。

在这里插入图片描述

3参考文献

[1]苏攀,张伟,李强,李世港.基于改进粒子群算法的PID参数优化研究[J].软件导刊,2020,19(10):94-97.
[2]张继荣,张天.基于改进粒子群算法的PID控制参数优化[J].计算机工程与设计,2020,41(04):1035-1040.
[3]https://blog.csdn.net/weixin_44044411/article/details/91353491
[4]https://yarpiz.com/440/ytea101-particle-swarm-optimization-pso-in-matlab-video-tutorial
[5]https://blog.csdn.net/weixin_43145941/article/details/104758286