PID传递函数模块及Matlab调试和传递函数离散化

1 simulink中微分模块、传递函数模块使用

大多数物理系统可以用微分方程来描述，因此可以用连续系统模拟。最简单的模型为线性模型和定常模型。

在simulink中，用来模拟连续系统的模块有四种：增益模块、求和模块、微分模块、积分模块。另外，传递函数模块也常常用来模拟物理系统和控制器。

1.1 积分、微分模块

积分模块：

  • 定义：计算输入信号从起始时间到当前时间对时间的积分，即：对输入信号积分。
  • 需要初始化条件。
  • 连续状态。

微分模块：

  • 定义：计算输入对时间的变化率。
  • 根据输出的返回的差值来拟合输入变化的速率。

实例：在simulink中选择正弦信号作为输入信号，并选择微分模块作为微分程序，利用示波器输出微分后的信号波形。

输出波形：

输入波形为黄色，输出波形为蓝色。输入波形为正弦信号，由于微分模块作用，因此输出模块为余弦信号，在本例中由于0时刻，初始值为0，因此输出信号波形在0时刻有突变。


1.2 传递函数模块使用及举例应用

传递函数模块表示法频繁用于控制系统设计和系统的动态模拟，传递函数的定义为系统在零初始状态下输出的Laplace变换与输入的Laplace变换之比。因此，传递函数是一种描述系统动力学输入输出关系的简便方法。

在本例中，以二阶低通滤波器为例，在simulink中应用传递函数模块，阶跃信号作为输入信号，并用示波器展示输出信号波形。

其中，传递函数为：(-1) / (s^2 + 3s + 1)，并在simulink中的传递函数模块设置分子、分母参数。

输出波形：随时间发生衰减，并最终稳定于-1。


1.3 Matlab中波特图使用方法

首先，在Matlab命令窗口中输入help bode学习下如何使用波特图。得到help文档：

在该文档中可以看出，在使用波特图之前，需要先定义传递函数，上图红箭头指的是help文档给我们举的例子，以下是我们利用上面的传递函数得到的波特图，传递函数：(-1) / (s^2 + 3s + 1)

首先，在命令行窗口，如上图箭头指示那样，我先令g等于传递函数表达式，然后利用bode(g)指令绘制波特图，并得到如下图：

从这个图中可以看到：随着输入信号频率的增加，对应的通过的幅值信号衰减比较严重，而在频率较低的情况下，信号比较容易通过，因此为低通滤波器。


2 PID控制参数Matlab调试

2.1 序

首先最重要的是了解每个参数调节了系统响应的哪些属性，通过观察响应从而调节参数改变属性。

PID作用概述：

  1. P产生响应速度和力度，过小响应慢，过大会产生振荡，是I和D的基础。
  2. I在有系统误差和外力作用时消除偏差、提高进度，同时也会增加响应速度，产生过冲，过大会产生振荡。
  3. D抑制过冲和振荡，过小系统会过冲，过大会减慢响应速度。D的另外一个作用时抵抗外界的突发干扰，阻止系统的突变。

同时调节顺序是：P > I > D

PID调节目标：

  • 衰减比在4-10之间最佳，也就是响应曲线的前两个幅值B：B1的比值在4-10之间。
  • 稳态误差趋近于0。
  • 系统响应越快越好。


2.2 P参数选取

tip：在第一步牢记P产生响应速度和力度，过小响应慢，过大会产生振荡，是I和D的基础。

如果想自己调试尝试可以打开matlab，运行simulink，照着下面的图进行连接，如果想直接应用可以往后看。

图中的系统为一个PID控制二阶系统。

拿上面的系统进行举例，首先设定P=0.1，I=0，D=0观察响应。可以看到图像没有超调，说明P产生的响应速度和力度太小了。

P=1，I=0，D=0观察系统响应，超调量出现但是只有一个波形，同时也就意味着调节时间太慢了，继续加大P。

P=10，I=0，D=0，此时调节时间显著下降，可以看到此时的数量级已经调整完成，也就是P参数只需要微调。

P=100，I=0，D=0，系统开始变得振荡。

如果继续加大P，系统会达到一个临界值，产生等幅振荡，最后开始发散，如下图。


2.3 I的调节

tip：I在有系统误差和外力作用时消除、提高精度，同时也会增加响应速度，产生过冲，过大会产生振荡。

I主要调节稳态输出，消除扰动。由于系统没有扰动输入因此看不到I对于消除扰动的效果。P=10，I=10，D=0，此时I过大导致系统振荡加剧。

P=10，I=1，D=0，此时响应波形基本符合预期。观察稳态输出约为0.963左右。

P=10，I=0.1，D=0，可以看到几乎响应波形没有变化。说明在没有扰动的情况下I只要不过大影响不大。但是稳态输出变化为0.916。

P=10，I=0，D=0，稳态输出变为0.91左右。

最终我们可以通过I少量调节稳态输出的值，最终将稳态误差消除。关于I对波形影响的作用总结如下图：


2.4 D的调节

tip：D抑制过冲和振荡，过小系统会过冲，过大会减慢响应速度。D的另一个作用是抵抗外界的突发干扰，阻止系统的突变。

P=10，I=0.1，D=10，可以看到将所有的冲击都消除掉了。

P=10，I=0.1，D=1，消除冲击减弱，此时显然衰减比不符合要求。

P=10，I=0.1，D=1，此时基本符合要求。


2.5 总结

首先调节P的数量级达到一个只有2个左右明显峰值的波形，再调节I找到不会波形振荡也不会没有超调的的区间，在区间内找到一个I将稳态误差尽可能消除。最终使用D来控制衰减比和波形的峰值、超调量。最后根据要求的稳态值、调节时间、超调量、上升时间、峰值时间等指标进行微调达到目标。

参数整定找最佳，从小到大顺序查，
先是比例后积分，最后再把微分加，
曲线振荡很频繁，比例度盘要放大，
曲线漂浮绕大湾，比例度盘往小扳，
曲线偏离回复慢，积分时间往下降，
曲线波动周期长，积分时间再加长，
曲线振荡频率快，先把微分降下来，
动差大来波动慢，微分时间应加长，
理想曲线两个波，前高后低4比1，
一看二调多分析，调节质量不会低 。


3 传递函数如何离散化转化为差分方程

3.1 前言

我们在Matlab/Simulink做实验仿真的对象一般习惯使用s域传递函数，但是真正转化到单片机中时，又无法识别传递函数，因此需要进行传递函数的离散化，转化为差分方程来处理。

离散化的目的：

  • 将s域下的传递函数转换为离散的z域函数。
  • 将离散域下的函数转化为差分方程，然后在单片机中实现。

3.2 传递函数形式

一阶惯性环节形如以下形式：

如何将上式改写成差分方程，需要按照下述流程。S域传递函数——Z域离散函数——差分方程。第一步先将其离散化。


3.3 离散化方法

3.3.1 一阶向前差分

将传递函数中的s使用上式进行替换即可，接下来是推导过程。（T为离散化采样时间）


3.3.2 一阶向后差分


3.3.3 双线性变换

一阶向前、向后差分法的依据是积分的矩形法则，有时效果并不好，故出现了依据积分的梯形法则的双线性变换法（又称塔斯汀法）。


3.3.4 零阶保持器法

3.3.5 一阶保持器法

3.3.6 脉冲响应不变法


至此，将传递函数中的s全部替换，得到传递函数的离散域表示（即差分方程）。我们也可以使用matlab中的c2d函数来对传递函数离散化，调用格式如下。

1
[numZ denZ]=c2d(num,den,‘tustin’) %bilinear method

c2d()中第一个参数为s域传递函数，第二个参数为采样间隔时间，第三个参数表示离散的方法。

离散方法包括以下几种方法：

  • zoh零阶保持器
  • foh一阶保持器
  • tustin双线性变换法
  • matched零极点匹配法
  • impulse脉冲响应不变法

注意：得到的差分方程的表示如下，n必须大于等于m(因为未来时刻的信息无法获知)

与s域传递函数类似，z呈降幂排列，但是首相幂为0，无正幂项。

式中：

表示前n个采样时刻的信号值。

当从连续域到离散域的离散效果不好时，有如下改进的方法：

  • 选取更适合的离散化方法。
  • 提高采样频率。
  • 修正连续域设计，如增加稳定裕度指标等。

3.4 实例说明

3.4.1 实例一

假设当前我们获得的传递函数模型如下：

采用方法一：将模型中的s用方法一替代。

化简得：

接下来根据以下转化关系：

上式可写成：

到此差分转换结束。其余变换可参考同样方法。


3.4.2 实例二


3.5 Matlab辅助完成

s域的传递函数在matlab中可如此表示：

通过如下方式获得Z变换形式：

分子分母同除以z，获得以下形式：

进一步：

其他情况可参考。暂时到这里结束，提醒一下，单片机中采用差分方程切记设置初值，初始化初值。

1
2
3
4
5
% 例2
ts=0.001;%采样时间=0.001s 
sys=tf(400,[1,50,0]);%建立被控对象传递函数 
dsys=c2d(sys,ts,'z');%把传递函数离散化（问题1） 
[num,den]=tfdata(dsys,'v');%离散化后提取分子、分母

3.6 使用不同离散化方法的影响

下面我们首先使用不同离散化方法对传函G(s) = 1/(0.2s+1)进行离散，离散化采样时间Ts=0.05：

然后对比不同离散化方法后的bode图和未离散化的系统bode图：

从bode图中可以看出，前向差分法离散化以后偏离实际的系统更远，不能保证系统的稳定性。而后向差分法和双线性变换法，更接近未离散化的bode图，所以后向差分法和双线性变换法在工程中较为常用。


3.6.1 传递函数到单片机实现的例子

假设一个低通滤波器的传函为：

我们先用matlab看一下该传函的特性：

1
2
3
4
5
6
>> num=1;
den=[6.687e-4,1];
sys=tf(num,den);
p=bodeoptions;
p.FreqUnits= 'Hz';
bode(sys,p)

看一下它的bode图：

从图中可以看到该低通滤波器截止频率为237Hz，截止频率处相移位45度。

我们也用simulink搭建一个仿真，验证下当该滤波器输入一个237Hz的正弦波时，输出是不是一个为-3dB衰减，且相移45度的正弦波。

波形如下：

图中sine Wave和sine Wave1两个波形相差90度，而滤波出来的Transfer Fcn波形可以正好相差45度，衰减-3dB。

我们现在用matlab对他进行离散化：

1
2
3
4
5
6
7
8
9
10
11
num=1;
den=[6.687e-4,1];
sys=tf(num,den);
ts=64e-6;
>>  dsys=c2d(sys,ts,'z') % 注意连续转离散（c2d），和离散转连续（d2c）的用法

dsys =
 
   0.09127
  ----------
  z - 0.9087

然后我们用simulink再搭建他对应的离散传递函数下的系统特性有没有变化：

采样时间为ts=64e^-6，波形如下：

从图中可以看出滤波出来的Discrete Transfer Fcn波形正好相差45度，衰减-3dB，所以它在离散下特性没变。

下面将离散域下的函数转化为差分方程：

假设输入为x，输出为y。

转化为上面的差分方程后，就可以在单片机中实现了，我们也可以在simulink中用sfun仿真，看下转化为数字滤波后的特性是否发生变化。同样输入幅值为1，频率为237Hz的正弦波，仿真如下：

sfun如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
function y = fcn(u)
persistent yl;  %定义一个静态变量
if isempty(yl)  %给静态变量赋初值
   yl=0;
end
persistent ul;
if isempty(ul)
    ul=0;
end

y=0.9086*yl+0.09136*ul;
yl=y;
ul=u;

仿真波形如下：

从图中可以看出滤波出来的Matlab Function波形正好相差45度，衰减-3dB。

注意：可能有的人会问，一阶低通滤波器的差分方程的形式不应该是：

是的，这个是用一阶后项差分法离散后得出的差分方程。

因为没有找到，matlab中传递函数怎么进行一阶后项差分离散（也就是用 s = (1-z^-1) / Ts 替代），其实我们也可以定义符号变量进行替代。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
>> syms s z Ts;%定义符号变量 
s=(1-z^-1)/Ts;%一阶后项差分
sys=1/(s+1);
subs(sys,[s],[z])%传递函数中用z替代s
 
ans =
 
-1/((1/z - 1)/Ts - 1) %这里我们可以对符号表达式进行各种化简
                      %参见matlab符号表达式的化简常用命令函数
 simplify(ans) %对结果进行化简
                     
 
ans =
 
(Ts*z)/(z + Ts*z - 1)
 
>> pretty(ans)%美化下结果
    Ts z
------------
z + Ts z - 1

化成差分方程的形式：(化简后的公式，第二个+号改为×号)

以上整个就是传递函数的离散化的单片机实现过程了。


3.7 二阶以上传函的离散化（使用状态方程的方式）

传递函数要编程实现，就必须先进行离散化。前面我们介绍了一阶传递函数的离散化，其过程比较简单。那么对于二阶及二阶以上的传递函数怎么离散化呢？我们可以把传递函数转化为状态方程，这样就可以变成多个一阶的传递函数，然后再进行离散化。

传递函数离散化步骤：

  1. 传递函数转化为状态方程。
  2. 状态方程离散化。

3.7.1 传递函数转化为状态方程

这里我们举一个简单的例子，把下面的传递函数转化为状态方程的形式：

这里我们借助matlab工具求解，代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
>> num=[2 1];
>> den=[1 7 9];
>> sys=tf(num,den)

sys =
 
     2 s + 1
  -------------
  s^2 + 7 s + 9
 
Continuous-time transfer function.

>> [A B C D]=tf2ss(num,den)

A =

    -7    -9
     1     0

B =

     1
     0

C =

     2     1

D =

     0

>>


3.7.2 状态方程离散化

我们把上面的状态方程离散化，写出matlab function模块仿真如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
function y = fcn(u,Ts)

persistent x1 x2;
if isempty(x1)
    x1=0;
end
if isempty(x2)
    x2=0;
end

x1 = x1 + Ts*x2;
x2 = x2 + Ts*(-9*x1-7*x2+u);
y = x1+2*x2;

end

3.7.3 传递函数离散化的仿真验证

我们把上面例子中的传递函数，状态方程和离散化的代码进行仿真验证下查看其阶跃响应：

三种方式的阶跃响应结果完全重合如下：