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

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)=\frac{0.998}{0.021s+1}$$因此通过tf函数建立系统结构体如下：

sys=tf(0.998,[0.021,1]);   %建立被控对象传递函数，即式4.1

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

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

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变换表达式为$dsys=\frac{num(1)\cdot z+num(2)}{den(1)\cdot z+den(2)}=\frac{0.2114}{z-0.7881}$
在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$和$y_1$保存了上一个数据，仅仅是为了演示这一过程。\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)-e_1$等效即可。积分的实现在每一次运算的后面都累加原来的误差，即Ee=Ee+e_1;即可。PID的控制器输出$$u(k)=Kp\cdot e(k)+Kd\cdot (e(k)-e_1)+Ki\cdot Ee$$
PID控制器构造完毕，我们需要通过r(k)和y(k)得到e(k)，再通过e(k)得出u(k)，进而再求解出y(k)，再结合r(k)求解出e(k),…以此循环，求解出离散的响应点。
详细的代码如下：

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比较图的代码如下：

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_{max}$时，若e(k)>0即输出值还未到达指定值，则认为积分会带来滞后，不再积分。
当$u(k)<0$时，若e(k)<0即输出值超过了指定值，则认为积分会带来滞后，不再积分。
在本案例中认为$u_{max}=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的手段，既消除了稳态误差还削弱了积分环节带来的滞后影响。

5.simulink仿真

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

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

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

得到的示波器曲线如下：

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