Matlab仿真PID控制（带M文件、simulink截图和参数分析）

2022-09-06 11:16:12 浏览数 (3)

大家好，又见面了，我是你们的朋友全栈君。

文章目录

  • 0.符号说明
  • 1.如何根据连续系统建立差分方程
    • 1.1.获取连续系统的传递函数
    • 1.2.获取离散系统的传递函数
    • 1.3.转换为差分方程
  • 2.基本PID控制原理
  • 3.比较PID输出，分析参数产生的影响
  • 4.改进PID算法（遇限削弱积分法）
  • 5.simulink仿真

0.符号说明

  1. y(k)——系统响应输出的离散值
  2. u(k)——数字PID控制输出的离散值
  3. r(k)——期望输出的离散值（事先已知），在本例中为常数（即阶跃输入）
  4. e(k)——e(k)=r(k)-y(k)，为期望值-实际值，是单位负反馈的误差比较信号

注：图片来源于百度百科

1.如何根据连续系统建立差分方程

1.1.获取连续系统的传递函数

线性定常系统的控制中，PID是个非常常见的控制方式，如果可以通过Matlab仿真出PID的控制效果图，那么对系统设计时的实时调试将会容易得多。在这里我们将会以一个利用系统辨识参数的PID设计为为例展示Matlab仿真PID的过程。 首先需要对一个未知的系统的参数进行辨识，以延迟环节可以忽略不计的电机调速系统为例。将时间戳导入xdata向量，对应的时刻转速导入ydata向量，进行系统辨识

链接：Matlab的系统辨识

我们就以上文链接中辨识的系统传递函数为例： G ( s ) = 0.998 0.021 s 1 G(s)=frac{0.998}{0.021s 1} G(s)=0.021s 10.998​因此通过tf函数建立系统结构体如下：

代码语言：javascript复制
sys=tf(0.998,[0.021,1]);   %建立被控对象传递函数，即式4.1

1.2.获取离散系统的传递函数

由于是数字PID仿真，我们需要选取一个采样时间，本案例选用的是0.005s（注意，采样周期应该小于系统纯滞后时间的0.1倍）。在对其进行数字PID控制前，我们需要将这个系统离散化：

代码语言：javascript复制
ts=0.005;  %采样时间=0.005s
dsys=c2d(sys,ts,'z');      %离散化

dsys即我们根据采样周期离散化的Z变换系统。首先我们需要提取这个Z变化d那系统的参数方便后面的计算：

代码语言：javascript复制
[num,den]=tfdata(dsys,'v');%'v'代表强制以向量的格式（默认为元胞数组）输出num和den

1.3.转换为差分方程

求解出的Z变换表达式为 d s y s = n u m ( 1 ) ⋅ z n u m ( 2 ) d e n ( 1 ) ⋅ z d e n ( 2 ) = 0.2114 z − 0.7881 dsys=frac{num(1)cdot z num(2)}{den(1)cdot z den(2)}=frac{0.2114}{z-0.7881} dsys=den(1)⋅z den(2)num(1)⋅z num(2)​=z−0.78810.2114​ 在PID仿真的过程中我们需要求解出时域表达式 ，因此需要借助差分方程解决，对于以下的Z变换：

begin{equation} Y(z)=dsyscdot U(z)=frac{num(2)}{den(1)cdot z den(2)}cdot U(z) label{eq:Sample1} end{equation}

begin{equation} zY(z) den(2)Y(z)=num(1)zU(z) num(2)U(z) label{eq:Sample2} end{equation} 对上式进行反Z变换，可以得到以下的差分方程：

begin{equation} y(k 1) den(2)y(k)=num(1)u(k 1) num(2)u(k) label{eq:Sample3} end{equation}

begin{equation} y(k 1)=-den(2)y(k) num(1)u(k 1) num(2)u(k) label{eq:Sample4} end{equation} 位置型PID仿真时实际上可以不需要保存前一个数据(u(k)和y(k))，增量型PID必须要保存前一个数据。这里我们使用了位置型PID，但仍然利用 u 1 u_1 u1​和 y 1 y_1 y1​保存了上一个数据，仅仅是为了演示这一过程。begin{equation} y(k 1)=-den(2)y(k) num(1)u(k 1) num(2)u(k) end{equation} 可以转换为下面的式子： begin{equation} y(k)=-den(2)y_1 num(1)u(k) num(2)u_1 label{eq:Sample5} end{equation} 我们的差分方程就这样建立完毕。注意，此差分方程仅仅是描述系统模型的运算规律的，和我们的控制无关。因此是y(k)和u(k)的映射关系。我们下面的控制则是利用负反馈信号e(k)导出u(k)的输出，求解的是控制器u(k)的序列值。

2.基本PID控制原理

以位置型PID控制为例。将连续的PID控制转换为数字式时，微分环节被用差分代替，积分环节被累加和代替，比例环节则保持不变。差分的实现非常简单，只需要用 e ( k 1 ) − e ( k ) e(k 1)-e(k) e(k 1)−e(k)即 e ( k ) − e 1 e(k)-e_1 e(k)−e1​等效即可。积分的实现在每一次运算的后面都累加原来的误差，即Ee=Ee e_1;即可。PID的控制器输出 u ( k ) = K p ⋅ e ( k ) K d ⋅ ( e ( k ) − e 1 ) K i ⋅ E e u(k)=Kpcdot e(k) Kdcdot (e(k)-e_1) Kicdot Ee u(k)=Kp⋅e(k) Kd⋅(e(k)−e1​) Ki⋅Ee PID控制器构造完毕，我们需要通过r(k)和y(k)得到e(k)，再通过e(k)得出u(k)，进而再求解出y(k)，再结合r(k)求解出e(k),…以此循环，求解出离散的响应点。 详细的代码如下：

代码语言：javascript复制
ts=0.005;  %采样时间=0.005s
sys=tf(0.998,[0.021,1]);   %建立被控对象传递函数，即式4.1
dsys=c2d(sys,ts,'z');      %离散化
[num,den]=tfdata(dsys,'v');   %
e_1=0;      %前一时刻的偏差      
Ee=0;       %累积偏差
u_1=0.0;    %前一时刻的控制量
y_1=0;       %前一时刻的输出
%PID参数
kp=0.22;    
ki=0.13;
kd=0;
u=zeros(1,1000);%预先分配内存
time=zeros(1,1000);%时刻点（设定1000个）
for k=1:1:1000
    time(k)=k*ts;   %时间参数
    r(k)=1500;      %期望值
    y(k)=-1*den(2)*y_1 num(2)*u_1 num(1)*u(k);%系统响应输出序列
    e(k)=r(k)-y(k);   %误差信号
    u(k)=kp*e(k) ki*Ee kd*(e(k)-e_1); %系统PID控制器输出序列
    Ee=Ee e(k);    %误差的累加和
    u_1=u(k);    	%前一个的控制器输出值
    y_1=y(k);    	%前一个的系统响应输出值
    e_1=e(k);		%前一个误差信号的值
end
%（仅绘制过渡过程的曲线，x坐标限制为[0,1]）
p1=plot(time,r,'-.');xlim([0,1]);hold on;%指令信号的曲线（即期望输入）
p2=plot(time,y,'--');xlim([0,1]);%不含积分分离的PID曲线
hold on;

输出的PID控制曲线如下：

3.比较PID输出，分析参数产生的影响

一个基本的PID就完成了。下面如果我们想要知道修改PID的三个参数kp,ki,kd会带来什么效果，只需要在程序中修改即可。为了方便起见，我们建立一个PID的数组，kp,ki,kd每次都取数组的一个值，然后设定一个大循环开始循环仿真。再利用subplot输出子图的方式将所有的PID效果都输出到一个图进行对比。该代码根据上述代码修改已经很容易，PID比较图的代码如下：

代码语言：javascript复制
close all
PID=[0.22,0.13,0;
    0.4,0.13,0;
    0.4,0.25,0;
    0.8,0.23,0.4;
    0.8,0.2,1;
    0.7,0.2,0.9];%初始化PID参数
for pid=1:1:6
ts=0.005;  %采样时间=0.005s
sys=tf(0.998,[0.021,1]);   %建立被控对象传递函数，即式4.1
dsys=c2d(sys,ts,'z');      %离散化
[num,den]=tfdata(dsys,'v');   %
e_1=0;      %前一时刻的偏差      
Ee=0;       %累积偏差
u_1=0.0;    %前一时刻的控制量
y_1=0;       %前一时刻的输出
%PID参数
kp=PID(pid,1);    
ki=PID(pid,2);
kd=PID(pid,3);
u=zeros(1,1000);
time=zeros(1,1000);
for k=1:1:1000
    time(k)=k*ts;   %时间参数
    r(k)=1500;      %给定量
    y(k)=-1*den(2)*y_1 num(2)*u_1 num(1)*u(k);
    e(k)=r(k)-y(k);   %偏差
    u(k)=kp*e(k) ki*Ee kd*(e(k)-e_1);   
    Ee=Ee e(k);    
    u_1=u(k);    
    y_1=y(k);    
    e_1=e(k);
end
subplot(2,3,pid);
p1=plot(time,r,'-.');xlim([0,1]);hold on;
p2=plot(time,y,'--');xlim([0,1]);
title(['Kp=',num2str(kp),' Ki=',num2str(ki),' Kd= ',num2str(kd)]);
hold on;
end

输出的子图矩阵如下：

可以发现，修改Kp会造成上升时间的缩短，但是有可能也会带来较大的超调。积分的增加是一个严重的滞后环节，会减小相位裕度，也会带来超调（超调量并不是绝对的，相对于较小的Kp可能会产生较大的超调，而Kp较大时超调会减小（例如第一行的1图和2图的对比））。然而积分的引入也是必要的，否则将会很长时间无法削弱误差e(k)（例如第二行第二个图）。微分的引入相当于一个超前校正，会减少超调，但是过渡的微分很可能会造成尾部振荡，系统逐渐变得不稳定。因此微分和积分之间需要一个平衡，当满足这个平衡的时候，系统几乎没有振荡，同时响应速度也较快。（第一行的图3是积分过多，产生超调，第二行的图1和图3就比较理想） 综合上述，PID的调节经验可以归结为以下几点：

  • Kp较小时，系统对微分和积分环节的引入较为敏感，积分会引起超调，微分可能会引起振荡，而振荡剧烈的时候超铁也会增加。
  • Kp增大时，积分环节由于滞后产生的超调逐渐减小，此时如果想要继续减少超调可以适当引入微分环节。继续增大Kp系统可能会不太稳定，因此在增加Kp的同时引入Kd减小超调，可以保证在Kp不是很大的情况下也能取得较好的稳态特性和动态性能。
  • Kp较小时，积分环节不宜过大，Kp较大时积分环节也不宜过小（否则调节时间会非常地长），在下面这个例子中我们还会介绍到，当使用分段PID，在恰当的条件下分离积分，可以取得更好的控制效果。原因在于在稳态误差即将满足要求时，消除了系统的滞后。因此系统超调会明显减少。本例中采样的抗积分饱和的方法是遇限削弱积分法。

4.改进PID算法（遇限削弱积分法）

遇限削弱积分法的原理是 当 u ( k ) > u m a x u(k)>u_{max} u(k)>umax​时，若e(k)>0即输出值还未到达指定值，则认为积分会带来滞后，不再积分。 当 u ( k ) < 0 u(k)<0 u(k)<0时，若e(k)<0即输出值超过了指定值，则认为积分会带来滞后，不再积分。 在本案例中认为 u m a x = r ( k ) u_{max}=r(k) umax​=r(k) 改进PID算法如下（需要些两个循环，当然也可以用一个循环，将其中的PID设为一个子过程调用）：

代码语言：javascript复制
close all
ts=0.005;  %采样时间=0.005s
sys=tf(0.998,[0.021,1]);   %建立被控对象传递函数，即式4.1
dsys=c2d(sys,ts,'z');      %离散化
[num,den]=tfdata(dsys,'v');   %
e_1=0;      %前一时刻的偏差      
Ee=0;       %累积偏差
u_1=0.0;    %前一时刻的控制量
y_1=0;       %前一时刻的输出
%PID参数
kp=0.22;    
ki=0.13;
kd=0;
u=zeros(1,1000);
time=zeros(1,1000);
for k=1:1:1000
    time(k)=k*ts;   %时间参数
    r(k)=1500;      %给定量
    y(k)=-1*den(2)*y_1 num(2)*u_1 num(1)*u(k);
    e(k)=r(k)-y(k);   %偏差
    u(k)=kp*e(k) ki*Ee kd*(e(k)-e_1);   
    Ee=Ee e(k);    
    u_1=u(k);    
    y_1=y(k);    
    e_1=e(k);
end
p1=plot(time,r,'-.');xlim([0,1]);hold on;
p2=plot(time,y,'--');xlim([0,1]);
hold on;
a=1;%控制积分分离的二值数
e_1=0;Ee=0;u_1=0.0;y_1=0;%重新初始化       
for k=1:1:1000
    time(k)=k*ts;   %时间参数
    r(k)=1500;      %给定量
    y(k)=-1*den(2)*y_1 num(2)*u_1;
    e(k)=r(k)-y(k);   %偏差
    u(k)=kp*e(k) ki*Ee kd*(e(k)-e_1);   
     if ((u(k)>r(k)) && (e(k)>0))||((u(k)<0) && (e(k)<0))
         a=0;
     else 
         a=1;
     end     
    Ee=Ee a*e(k);    
    u_1=u(k);    
    y_1=y(k);    
    e_1=e(k);
end
p3=plot(time,y,'-');xlim([0,1]);
title('含积分分离与不含积分分离的对比');
legend([p1,p2,p3],'指令信号','不含积分分离','含积分分离');

输出的曲线对比图如下：

可以发现，系统的超调量明显减少了，调节时间也减少了一点。原因在于我们采用了分段PID的手段，既消除了稳态误差还削弱了积分环节带来的滞后影响。

5.simulink仿真

需要的模块名称（不区分大小写）如下：

  • gain（参数分别为0.22和0.13/0.005）
  • sum（参数分别为”|±”和”| ”）
  • integrator
  • scope 注意：本文使用的是离散PID仿真，而simulink使用的是连续系统仿真，转换PID参数时P参数不变，I参数应该除以仿真间隔Ts=0.005，D参数应该乘Ts。

以表中第一组PI参数为例：

得到的示波器曲线如下：

希望本文对您有帮助，谢谢阅读。

发布者：全栈程序员栈长，转载请注明出处：https://javaforall.cn/134823.html原文链接：https://javaforall.cn

matlab equation label pid simulink

1 人点赞