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

Matlab仿真PID控制(带M文件、simulink截图和参数分析)在直接数字域设计中,我们常常需要用到PID算法,而PID算法投入单片机使用时,往往需要硬件的支持,在调试时非常麻烦。本文通过Matlab仿真的手段实现PID,方便了开发者对系统的设计和实时调试。

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函数建立系统结构体如下:

sys=tf(0.998,[0.021,1]);   %建立被控对象传递函数,即式4.1
只听到从架构师办公室传来架构君的声音:
自明及晦,所行几里?夜光何德,死则又育?厥利维何,而顾菟在腹?女岐无合,夫焉取九子?伯强何处?惠气安在?何阖而晦?何开而明?角宿未旦,曜灵安藏?不任汩鸿,师何以尚之?佥曰“何忧,何不课而行之?”鸱龟曳衔,鲧何听焉?顺欲成功,帝何刑焉?永遏在羽山,夫何三年不施?伯禹愎鲧,夫何以变化?纂就前绪,遂成考功。有谁来对上联或下联?

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

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

此代码由Java架构师必看网-架构君整理
ts=0.005; %采样时间=0.005s dsys=c2d(sys,ts,'z'); %离散化

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

[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)=z0.78810.2114
在PID仿真的过程中我们需要求解出时域表达式 ,因此需要借助差分方程解决,对于以下的Z变换:

\begin{equation}
Y(z)=dsys\cdot 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)=Kp\cdot e(k)+Kd\cdot (e(k)-e_1)+Ki\cdot Ee u(k)=Kpe(k)+Kd(e(k)e1)+KiEe
PID控制器构造完毕,我们需要通过r(k)和y(k)得到e(k),再通过e(k)得出u(k),进而再求解出y(k),再结合r(k)求解出e(k),…以此循环,求解出离散的响应点。
详细的代码如下:

此代码由Java架构师必看网-架构君整理
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控制曲线如下:
PID控制

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

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

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

输出的子图矩阵如下:
PID子图矩阵
可以发现,修改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设为一个子过程调用):

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

5.simulink仿真

需要的模块名称(不区分大小写)如下:

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

以表中第一组PI参数为例:
在这里插入图片描述
得到的示波器曲线如下:
在这里插入图片描述

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

本文来源非线性光学元件,由架构君转载发布,观点不代表Java架构师必看的立场,转载请标明来源出处:https://javajgs.com/archives/209481
0

发表评论