连续传递函数离散化的方法与原理

目录

第一章 模拟化设计基础 第一节 步骤

第二节 在MATLAB中离散化 第三节 延时e-Ts环节的处理 第四节 控制函数分类 第二章 离散化算法 摘要 比较

第一节 冲击响应不变法(imp,无保持器直接z变换法) 第二节 阶跃响应不变法(zoh,零阶保持器z变换法) 第三节 斜坡响应不变法(foh,一阶保持器z变换法) 第四节 后向差分近似法 第五节 前向差分近似法 第六节 双线性近似法(tustin) 第七节 预畸双线性法(prevarp) 第八节 零极点匹配法(matched) 第三章 时域化算法

第一节 直接算法1—双中间变量向后递推 第二节 直接算法2—双中间变量向前递推 第三节 直接算法3—单中间变量向后递推

第四节 直接算法4—单中间变量向前递推(简约快速算法) 第五节 串联算法 第六节 并联算法 第四章 数字PID控制算法 第一节 微分方程和差分方程 第二节 不完全微分 第三节 参数选择 第四节 c51框架 第五章 保持器 第一节 零阶保持器 第二节 一阶保持器

附录 两种一阶离散化方法的结果的比较

1 1 3 5 6 10 10 11 11 11 11 12 14 15 17 18 19 19 20 21 21 22 23 24 25 25 26 27 33 33 30 31

第一章 模拟化设计基础

数字控制系统的设计有两条道路，一是模拟化设计，一是直接数字设计。如果已经有成熟的模拟控制器，可以节省很多时间和部分试验费用，只要将模拟控制器离散化即可投入应用。如果模拟控制器还不存在，可以利用已有的模拟系统的设计经验，先设计出模拟控制器，再进行离散化。

将模拟控制器离散化，如果用手工进行，计算量比较大。借助数学软件MATLAB控制工具箱，可以轻松地完成所需要的全部计算步骤。如果需要的话，还可以使用MATLAB的SIMULINK工具箱，进行模拟仿真。

第一节 步骤

步骤1 模拟控制器的处理

在数字控制系统中，总是有传输特性为零阶保持器的数模转换器（DAC），因此，如果模拟控制器尚未设计，则应以下

1eTs图的方式设计模拟控制器，即在对象前面加上一个零阶保持器，形成一个新对象G(s)，然后针对这个新对象求模拟

s控制器D(s)。事实上，模拟控制器一般是已经设计好的，无法或不方便更改了，离散化后的系统只好作为近似设计了。 然而，按照上述思路，可否将已有的控制器除以一个零阶保持器再离散化呢？还没有这方面的实际经验。

x-模拟控制器e1-es-TsD(s)G(s)u对象

以下假设选定的G(s)，D(s)如下图，而且不对G(s)作添加保持器的预处理。

x-es+2D(s)=8.s+1520G(s)=s(s+2)u

步骤2 离散化模拟控制器

离散化模拟控制器之前，先要确定离散化算法和采样时间。离散化算法有好几种，第二章中有详细的论述，现假定采用双线性变换法。确定采样时间，需要考虑被控对象的特性，计算机的性能，以及干扰信号的影响等，初步可按采样时间T<0.1Tp，Tp为被控对象时间常数，或T=(0.125～0.25)τ，为被控对象的纯滞后，初步确定后再综合平衡其它因素，当然这需要一定的经验，现在假定取0.05秒。 假设模拟控制器为D(s)8

转换结果为：

x-es+2D(s)=8.s+15D(z)=6.1091(z-0.9048)z-0.4545G(s)=20s(s+2)us2，在MATLAB中，用c2d函数进行离散化，过程为：

s15ds=zpk(-2,-15,8) %建立模拟控制器的s传递函数 dz=c2d(ds,0.05,'tustin') %将模拟控制器按tustin方法转换为z传递函数的数字控制

步骤3 检验数字控制器的性能

数字控制器的性能项目比较多，我们仅以直流增益，频率特性，零极点分布说明。 直流增益 dcgain(dz)

返回直流增益1.0667

以连续传递函数为基础的数字控制设计1/35

频率特性 bode(ds,'r',dz,'g') 伯德图，见下页左图 零极点分布 pzmap(dz) 步骤4 离散化控制对象

为了进行模拟仿真，需要对控制对象进行离散化，由于步骤1所说的原因，应把被控对象视为零阶保持器与原对象的串

零极点分布图，见下页右图

1eTs连，即应对G(s)进行离散化，这时可在c2d函数中使用零阶保持器(zoh)方法，如果认为不需要添加零阶保持器，即

s直接对G(s)离散化，则应在c2d函数中使用冲击响应不变法（imp）。 借用零阶保持器(zoh)方法，将对象G(s)

转换结果为：

e-D(z)=6.1091(z-0.9048)z-0.45451-eG1(z)=Z(s-Ts20带一阶保持器离散化的过程如下：

s(s2)...... %模拟控制器D(s)转换为D(z)的过程见前 gs=zpk([ ],[0,-2],20) %建立对象的s传递函数 g1z=c2d(gs,0.05,'zoh') %借用c2d函数进行带零阶保持器的对象的离散化 xD(z)G(s))uG1(z)=0.024187(z-0.9672)(z-1)(z-0.9048)

步骤5 模拟仿真

求离散系统的闭环传递函数和连续系统的闭环传递函数。 离散系统的闭环传递函数为：TRCZ 连续系统的闭环传递函数为：TRCS 用MATLAB算TRCZ与TRCS：

trcz=dz*g1z/(1+dz*g1z) trcs=ds*gs/(1+ds*gs) D(z)G1(z)

1D(z)G1(z)D(s)G(s)

1D(s)G(s) 结果为： TRCZ 0.14776(z-0.9048)(z-0.9048)(z-1)(z-0.4545)(z+0.9672)(z-0.4545)(z-0.9047)(z-0.9048)(z-1)(z2-1.307z+0.5975)

以连续传递函数为基础的数字控制设计2/35

TRCS160s(s+2)2(s+15)s(s+15)(s+2)2(s2+15s+160)

用MATLAB函数STEP画阶跃响应图形：

响应图形为：

hold on %图形保持 step(trcs,'r',2) %画连续系统的阶跃响应图,红色,终止时间为2秒 step(trcz,'b',2) %画离散系统的阶跃响应图,兰色,终止时间为2秒

步骤6 求数字控制器的时域表达式 上面已经求出， 连续传递函数D(s)8s2的tustin离散式为 s15

.091e(k)5.527e(k1)0.4545u(k1)，根据此式就可以写出计算u(k)的 对上式取z反变换，得时域表达式u(k)61程序代码来了。

除上述步骤之外，在编写程序代码时，还需要考虑几个问题： ① ADC位数

ADC位数是一个硬件问题，在系统设计时，就应该结合控制算法，仔细分析ADC位数对控制精度的影响。 ② 数据类型

在控制算法中有3种数据类型可以采用：无符号整数，定点数，浮点数。我们知道，无符号整数运算既是最简单和速度最快的运算，又是定点数运算和浮点数运算的基础。在数据动态范围比较小的情况下，应尽可能用无符号整数运算代替定点数运算和浮点数运算。

浮点数运算是一整套运算，包括加，减，乘，除，对阶，规格化，溢出处理等一系列子程序，使用浮点数的程序，将占用比较多的代码空间和比较长的运行时间。不论使用汇编语言还是c语言，对于是否使用浮点数运算，都应进行比较仔细的酌斟。 ③ 数值计算误差

数值计算引入的误差有3个方面，一是定点数字长不够或者浮点数有效数字过少，一是两个相近的数相减，一是加减乘除次数过多。在程序设计中，应优化算法以避免计算误差的引入。

11U(z)61091.(z0.9048)6.1091-5.527z1U(z)61091E.(z)5.527E(z)z0.4545U(z)z，或 。 D(z)E(z)z0.454510.4545z1以连续传递函数为基础的数字控制设计3/35

第二节 在MATLAB中离散化

1. 建立s降幂传递函数 ① 建立多项式型s降幂传递函数 方法1. sys = tf(num,den) 方法2. s = tf('s')，再令sys = f(s) 例： 已知G015s.30.25s31.25s2s，用tf函数建立多项式型s降幂传递函数，G0tf([.153],[0.25125.10])，

若G0以零极点形式给出，即已知G03(0.5s1)，仍可用tf函数建立多项式型s降幂传递函数，但需用多项式乘

s(s1)(0.25s1)法函数conv配合，G0=tf(conv([3],[0.5 1]),conv(conv([1 0],[1 1]),[0.25 1])) ② 建立零极点型s传递函数 方法1. sys = zpk(z,p,k)

方法2. s = zpk('s')，再令sys = f(s)。 2. 传递函数的转换

① 将任意形式的s传递函数转换为s降幂传递函数 sys = tf(sys) ② 将任意形式的s传递函数转换为s零极点传递函数 3. 建立z传递函数

① 将连续传递函数转换为离散传递函数 sysd = c2d(sysc,t,method) ② 建立多项式形式的z传递函数 方法1. sys = tf(num,den,t) 方法2. z = tf('z',t)，再令sys = f(z) ③ 建立零极点z传递函数 方法1. sys = zpk(z,p,k,t)

方法2. z = zpk('z',t)，再令sys = f(z) 4. 建立z-1格式的传递函数

直接根据分子和分母建立z-1格式的传递函数

sys_z = filt(num,den,t)

sys = zpk(sys)

当已有多项式形式的z降幂传递函数时，按以下步骤：

① 取z降幂传递函数a的分子多项式系数 [num,den] = tfdata(sys_s,'v') ② 建立z-1格式的传递函数

sys_z = filt(num,den,t)

5. 将多项式形式的高阶z-1降幂传递函数转换为并联传递函数，按以下步骤： ① “手工”将z-1降幂传递函数a改写成多项式形式的z降幂传递函数b。 ② 取z降幂传递函数b的分子多项式系数num和分母多项式系数den。

③ 利用residue函数取z降幂传递函数b的分项分式，[an ad ak]=residue(num,den)。

分项结果可能出现共轭复数，在这种情况下应将含有共轭复数的分式合并成二次有理质分式。 ④ 根据分项结果手工写出z降幂多项式形式的并联表达式。

⑤ “手工”将z降幂多项式形式的并联表达式改写成z-1降幂多项式形式的并联表达式。

⑥ “手工”对z-1降幂多项式形式的并联表达式中的每一个分式项降阶，即将每一个分式项变形，使得分式项的分子的阶次比分母的阶次低1阶，变形完毕再将全部常数项合并。 举例

设有z-1降幂传递函数 tz_0101.0186z.0.0864z211.27z10.27z2

以连续传递函数为基础的数字控制设计4/35

改写成z降幂传递函数 tz_010.1z20.186z0.0864z21.27z0.27，

取分项矢量 [an ad ak]=residue(nz,dz) 得 an=[0.5101 -0.1971] ad=[1.0000 0.2700] ak=0.1

0.51010.51010.5101z10.5101 手工写分项分式 ， 令 0.5101tz_p1z11z11z11z1

0.730.19710.1971z10.73， 令 tz_p20.73z0.2710.27z110.27z110.27z1 又令 tz_p30.10.51010.730.3199

1tz_p20.3199 故 tz_0tz_p3tz_p 验证 0.31990.51011z10.731z10.51011z10.731z1

0.10.186z10.08637z211.27z10.27z2 实际上，大多数数字控制器的传递函数都是一阶或者二阶的，所以需要分解的并不是很多。 6. 将高阶z-1降幂传递函数生成串联传递函数

使用zpk函数 sys_zpk= zpk(sys_pl)

零极点的概念是相对于z而不是相对于z-1说的，但对于以z-1为变量的降幂传递函数sys_pl来说，仍然可以用zpk(sys_pl)

生成以z-1为变量的因式积形式传递函数，权且也称为零极点形式。zpk函数对于z降幂传递函和z-1降幂传都能得到合理的结果，原因是zpk函数的作用就是把分子多项式和分母多项式分别进行因式分解。

第三节 延时环节e-Ts的处理

在建立s传递函数的LTI模型时，对于延时环节e-T s，可按如下方法处理： 1. 在tf函数中使用属性’inputdely’或者’iodely’ ，例如：

>> tf([1 -1],[1 3 5],'inputdelay',0.35) 将返回以下形式的传递函数 exp(-0.35*s) * --------------

s – 1 s^2 + 3 s + 5

1eTs使用这个方法不能建立形如的传递函数，因为带延时的传递函数不能与不带延时的传递函数相加，但可以使

s用tf函数建立2个传递函数，主体部分相同，但一个无输入延时，一个有输入延时0.35s，

用c2d进行离散化，但要求延时时间t必须是采样时间t的整数倍，若不是整数倍，则在转换时不理会延时环节，例如：

>> a=tf([1 -1],[1 4 5]) s - 1 ---------------- s^2 + 4 s + 5 >> a1=tf([1 -1],[1 4 5],'iodelay',0.35) s - 1

exp(-0.35*s) * --------------- s^2 + 4 s + 5

若采样时间为0.05，因为延时时间是采样时间的整数倍，转换结果的主体部分完全一样： >> c2d(a,0.05,'imp') z^2 - 1.039 z + 9.146e-018

------------------------------ sampling time： 0.05

以连续传递函数为基础的数字控制设计5/35

z^2 - 1.807 z + 0.8187 >> c2d(a1,0.05,'imp') z^2 - 1.039 z + 9.146e-018

z^(-7) * ----------------------------- sampling time： 0.05 z^2 - 1.807 z + 0.8187

若采样时间为0.1，因为延时时间不是采样时间的整数倍，结果的主体部分不一样： >> c2d(a,0.1,'imp') z^2 - 1.06 z + 4.349e-018

----------------------------- sampling time： 0.1 z^2 - 1.629 z + 0.6703 >> c2d(a1,0.1,'imp') 0.768 z - 0.851

z^(-3) * -------------------- sampling time： 0.1 z^2 - 1.629 z + 0.6703 2. 将e-Ts有理化

设tsG(s)eTs，因e-Ts的一阶有理表达式是 eTs1221Tsss2T，故tsG(s)eTsG(s)Tn(s)。 122d(s)1Tsss2TT

为了对ts进行离散化，首先使用tf函数建立lti模型的ts。若tstf(n,d)，则tzc2d(ts,t,'imp')。 3. 在离散化时使用恒等式zneTstTst

me 设tsG(s)1Ts，因zne，采样时间为Δt，若T=mΔt，则tsG(s)z，因为离散化时总是认为T1t，

故取tsG(s)z。

0z 根据以上假设，用MATLAB的c2d函数对g(s)进行离散化，则Gz_0c2d(Gs,t,'imp')，进而GzGz_则sys_zsys_zazn1。

。

第四节 控制函数分类

以下函数在control toolbox中，这里所述仅限于siso模型 1 创建多项式形式的传递函数

sys = tf(num,den) 创建一个s降幂多项式连续传递函数sys，分子多项式系数和分母多项式系数分别为num和den。 sys = tf(num,den,t) 创建一个z降幂离散传递函数sys， sys = tf 创建一个空的tf对象。 sys = tf(m) 指定静态增益m。 2 创建零极点形式的传递函数

sys = zpk(z,p,k) 创建一个零极点模型的连续传递函数sys，零极点矢量分别是z和p，增益是k。

sys = zpk(z,p,k,t) 创建一个零极点模型的离散传递函数sys，零极点矢量分别是z和p，增益是k，采样时间是t。在零极点对象中，如果没有零点，则z=[]。 sys = zpk 创建一个空的零极点对象。 sys = zpk(d) 指定静态增益d。 3 创建任意形式的传递函数

以连续传递函数为基础的数字控制设计6/35

t是采样时间，行矢量num和den同上。

sys = tf(sys) 把一个任意的lti模型sys转换成多项式传递函数，例如把零极点模型转换成多项式传递函数。

s = tf('s') 指定多项式传递函数变量为s变量

z = tf('z',t) 指定多项式传递函数变量为z变量，t为采样时间 形如

10的传递函数，既不能直接用tf函数建立，也不能直接用zpk函数建立。

s(0.25s1)(0.05s1) 定义了s=tf('s')后，写出赋值式f=10/(s*(0.25*s+1)*(0.05*s+1))，回车后即可得到多项式型传递函数

f10800，再令f1=zpk(f)，可得f。 10.0125s^30.3s^2ss(s20)(s4) s = zpk('s') 指定零极点传递函数变量为s变量

z = zpk('z',t) 指定零极点传递函数变量为z变量，t为采样时间

定义了s = zpk('s')，写出赋值式f=10/(s*(0.25*s+1)*(0.05*s+1))，回车后即可得到零极点型传递函数f14 建立z-1降幂离散传递函数

sys = filt(num,den,t) 返回z-1降幂离散传递函数，num和den分别是分子和分母多项式系数，t是采样时间。 sys = filt(m) 返回增益离散传递函数。 5 连续函数离散化

sysd = c2d(sysc,t,method)把连续传递函数sysc转换成采样时间为

t的离散传递函数，字符串method为离散化方法：

800。

s(s20)(s4)'zoh'(零阶保持，即阶跃响应不变)，'foh'(一阶保持)，'imp'(冲击响应不变,v6以上版本)，'tustin' (双线性近似)，'prewarp'(带预畸变的双线性近似)，'matched'(零极点匹配)。

注： ① 缺省的方法是'zoh'

22D(s)(1z1)Ts1(1z1) ② 'foh'的算法是D(z)Z[2]而不是D(z)Z[D(s)，这一点可以通过验证证]TsTz1s2D(s)z22z1实，验证方法是，令D(z)c2d(2,，但此结果中，D(z)的分子和分母将含有公因式，所以应进一T,'imp')Tzs步用zpk函数把D(z)表示成零极点形式，然后用“手工”的方法写出不含公因式的D(z)来，可以看到，最后的结果与用foh方法得到的结果完全一致。

③ 当使用'prewarp'方法时，临界频率wc(in rad/sec)作为第四个输入来指定，如sysd = c2d(sysc,t,'prewarp',wc) 。另有1种形式是，opt = c2dOptions('Method','tustin','PrewarpFrequency',.5), c2d(ds,.05,opt)。 6 取多项式模型传递函数的分子和分母的系数矢量

[num,den] = tfdata(sys,'v') 对于siso模型sys返回作为分子和分母系数的单行矩阵num和den。 [num,den,t] = tfdata(sys,'v') 同上，同时返回采样时间7 取零极点模型传递函数的零点和极点的单行矩阵

[z,p,k] = zpkdata(sys,'v') 返回lti模型sys的零极点矢量z和p，增益k。 8 取e-T s近似式

s的n阶pade近似式，行矢量num和den是s的降幂多项式系数。 [num,den] = pade(t,n) 返回e-T·

t。

e-Ts的一阶有理表达式是 eTs121Tss2T。 121Tss2T231(Ts)(Ts)1Ts...2848 231(Ts)(Ts)1Ts...2848s的高阶有理表达式是 e e-T·

Ts9 画阶跃响应图

step(sys) 画出由tf，zpk，or ss等函数创建的lti模型sys的阶跃响应图。

step(sys,tfinal) 画出lti模型sys从t=0到t=tfinal的阶跃响应图。对于未指定采样时间的离散模型，tfinal被解释为采样

以连续传递函数为基础的数字控制设计7/35

的数目。

step(sys,t) 使用用户提供的矢量t画阶跃响应图。对于离散时间模型，t的形式应该是ti:t:tf，在这里，t是采样时间。对于连续时间模型，t的形式应该是ti:dt:tf，在这里，dt变成对于连续系统的近似离散化的采样时间。Ti是开始时间，tf是终止时间。因为阶跃输入总是假定在t=0开始，所以通常不考虑ti和tf，即只使用一个终止时间t。

step(sys1,sys2,...,t) 在一个单个的图上画出多个lti模型sys1，sys2，... 的阶跃响应图，时间矢量t是可选择，还可以以step(sys1,r,sys2,y,sys3,gx)的方式对每一个系统指定颜色，线型和标记。

[y,t] = step(sys) 返回用于仿真的时间t的输出响应y，但并没有图形画在屏幕上，如果sys有ny输出和nu输入和lt = length(t)，y就是一个尺寸为[lt ny nu]的阵列，而y(:,:,j)给出第j个输入通道的阶跃响应。 10 画脉冲响应图

impulse 脉冲响应函数，用法与step相同 11 画频率响应图--伯德图(连续或离散) bode(sys) 画伯德图

bode(sys,{wmin,wmax}) 在频率wmin，wmax(in radians/second)之间画伯德图 bode(sys,w) 按指定的频率矢量w(in radians/second)画伯德图 bode(sys1,sys2,...) 画多个lti模型sys1，sys2，...的伯德图

bode(sys1,sys2,...,w) 按指定的频率矢量w(in radians/second)画多个lti模型sys1，sys2，...的伯德图 以下函数在符号工具箱symbolick中，需注意，在使用这些函数前，要对所使用的变量进行符号说明，例如：

12 福里哀变换

fourier

syms a t %是用空格分隔而不是用逗号分隔 a=sin(t) L=laplace(a) 反福里哀变换 ifourier 13 拉普拉斯变换 laplace 反拉普拉斯变换 14 z变换

ztrans

ilaplace

反z变换 iztrans

注：以上3 种变换必须是符号表达式，例如： syms t laplace(sin(t))

15 改善公式的可读性 pretty 16 多项式转换为符号表达式 poly2sym 17 化简符号表达式

simplify

numden

18 取符号表达式的分子和分母

以下函数在符号工具箱polyfun中，在公式变换中可能会用到 19 部分分式展开 [an ad ak] = residue(n,d)

n和d分别为原分式的分子和分母矢量，an和ad分别为分项式的分子和分母矢量，ak为整式部分。这是一个数学公式，在数字控制器程序设计中，利用部分分式展开的方法，把高于2阶的分式变换为不高于2阶的分式之和，从而把高阶传递函数算法变为低阶传递函数并联的算法。在使用这个方法时，不论在分子矢量中还是在分母矢量中，如果有共轭复数出现，则应将其整合为2阶质因式。 20 多项式乘法 21 多项式除法

c = conv(a,b) a×b=c

[q,r] = deconv(b,a) b÷a=q...r，即b = conv(a,q) + r

以连续传递函数为基础的数字控制设计8/35

第二章 离散化算法

连续传递函数离散化的核心环节，就是将控制器的s传递函数转换为z传递函数。离散化后，系统应该仍有好的稳定性，好的控制精度，而不是要求转换前后，两个数学公式等值。因此，离散化方法有多种。本章对这些方法的转换原理和由来进行了演绎。

离散化后得到的离散传递函数的稳定性，没有进行讨论，仅列了一张表进行比较。离散的结果，虽然可能会控制精度降低，应该认为，这不是主要问题。真正影响控制精度的因素，主要还是采样周期的长短。

一般情况下，由连续到离散的设计最好多实验几种方法（通过仿真，得出满意的结果）。因为匹配零、极点映射法、双线性变换法都能得出比较满意的结果，初步设计时，可以试用这些方法。

而其实，后向差分近似法也是合理的选择。但MATLAB的c2d函数中没有这一方法，在该方法的介绍之后，给出了一个可由MATLAB引用的m文件函数。

摘要

3种保持器法 保持器法即将s函数串联上一个保持器后取z变换，也可以从“对输入信号的响应不变”的角度导出。 无保持器法

D(z)Z[D(s)] （阶跃响应不变）

D(z)D(z)imp

zoh

foh

零阶保持器法 一阶保持器法

D(z)2(z1)z11Z[D(s)] zs（阶跃响应不变）

2(z1)1Z[2D(s)] Tzs（斜坡输入不变）

Tz2Ts1Z[2D(s)]（一阶保持器）

sMATLAB中无此方法

4种近似法 近似法将s与z的无理关系s 后向差分近似法 D(z)D(s)|z1 Tzlnz近似地化为有理关系，主要应用于传递函数从连续到离散的变换。 TsMATLAB的C2D函数中没有这个方法

前向差分近似法 D(z)D(s)| 双线性近似法

z1sT MATLAB的C2D函数中没有这个方法

D(z)D(s)|s2z1 tustin Tz1z1 Tz1tan121 预畸双线性近似法 D(z)D(s)|s，1为进行预畸变的频率 prevarp

|D(z)|e 增益匹配： |D(s)|s0z1或|D(j1)||D(j1)|

1种匹配法 匹配法完全从控制学的角度看问题，也是应用于传递函数从连续到离散的变换。 零极点匹配法 若 则

D(s)D(z)Ks(szsz2)...(szm)1)((spsp2)...(spn)1)(， matched

nmKz(zez1t)(zez2t)...(zezmt)(z1)(zep1t)(zep2t)...(zepnt)(zep1t)(zep2t)...(zepnt)， ，

或者 D(z)nm1Kz(zez1t)(zez2t)...(zezmt)(z1)|D(z)| 确定增益kz： a 令|D(s)|s0z1，

b 若D(s)分子有s因子，例如D(s)s， s1j1|D(z)| 可以令 |D(s)|esz1，也可以令 |D(j1)||D(以连续传递函数为基础的数字控制设计9/35

)|

比较

离散化方法 变换公式 映射关系 特点 冲击响应不变法 Z变换法 imp 脉冲响应采样值相同； D(z)Z[D(s)] 容易产生频率混迭现象， Ws2，为采样角频率。 T阶跃响应不变法 零阶保持器法 zoh z1D(s)D(z)Z[] zsWs21Ws2阶跃响应采样值相同； 稳定增益不变。 斜坡响应不变法 一阶保持器法 foh D(z)2(z1)Tz22Ts1Z[2D(s)]s D(z)(z1)1Z[2D(s)] Tzs向前差分法 D(z)D(s)|sz1T 1T1D(s)稳定，D(z)可能不稳定；等效精度差。 变换计算简单； 向后差分法 D(z)D(s)|z1 sTz②如果D(s)稳定，D(z)稳定； 1③离散滤波器的过程特性及频率特性有一定的失真，需要 较小的采样周期T。 D(s)稳定，D(z)也稳定； 双线性变换法 tustin D(z)D(s)|s2z1 Tz1低频特失真，但无频率混迭现象。 稳定增益不变； 具有串联特性。 1 预修正双线性变换法 D(z)D(s)|s prevarp z1Tz11tan()21 有前一种变换的特点; 还能保证在关键频率ω1处，幅频特性不变。 Z域与s域零极点位置一一对零极点匹配法 matched 应； 1当没有零点时，补充z=1的零 点可避免频率混迭现象。

以连续传递函数为基础的数字控制设计10/35

第一节 冲击响应不变法(imp，无保持器直接z变换法)

D(z)Z[D(s)] 公式推导1— 无保持器，直接转换：

D(z)Z[D(s)]

公式推导2— 按冲击响应不变的原则：

D(z)Z[1D(s)]Z[D(s)] （冲击响应不变，Z[1D(s)]中的1表示L((0))） 1第二节 阶跃响应不变法(zoh，零阶保持器z变换法)

1eTsz1D(s)Ts1 D(z)Z[D(s)]Z[]，式中T为采样周期，所以ez szs公式推导1— 按串联零阶保持器的原则：

1eTs1eTsz11 因为零阶保持器的s传递函数为，故D(z)Z[D(s)]Z[D(s)]

sszs公式推导2— 按阶跃响应不变的原则： 设阶跃信号e(t)1，则E(s)1z，E(z)，按照下一节对斜坡响应不变法的推导，立即可以写出

z1s1Z[D(s)]Z[E(s)D(s)]z11D(z)sZ[D(s)]

zE(z)zsz1第三节 斜坡响应不变法(foh，一阶保持器z变换法)

2(z1)1D(z)Z[2D(s)]，按斜坡响应不变原则导出，也可按一阶保持器原则配合e-Ts的台劳近似式得出，foh方法使用此式 Tzs 221Ts1eTs2(1z1)Ts1(z1)Ts1D(z)Z[()D(s)]Z[D(s)]Z[D(s)]，从一阶保持器的思路推出 TsTs2Tz2s2 以上2个公式中T为采样周期，两式差异请参看“附录 两种一阶离散化方法的结果的比较”。 公式推导1— 按斜坡响应不变的原则：

设连续滤波器为D(s)，输入为斜坡函数e(t)=t，则E(s)1s2，采样输出为us(kT)。按z变换的定义，us(kT)的z变换就是

u(t)的z变换，故US(z)Z[1s2D(s)]。

又设离散滤波器的传递函数为D(z)，它的输出为uz(KT)，按斜坡响应不变的要求，应有uz(kT)=us(kT)，故

Tz1UZ(z)US(z)Z[2D(s)]，但离散滤波器的输出为U(z)=E(z)D(z)，而斜坡函数e(t)=t的z变换为E(z)，求Uz(z) 与2s(z1)1Z[2D(s)]2UZ(z)(z1)1sZ[2D(s)]。 E(z)的比值，得斜坡响应不变法离散公式 D(z)TzE(z)Tzs2(z1)以连续传递函数为基础的数字控制设计11/35

公式推导2— 按串联一阶保持器的原则：

1Ts1eTs2 一阶保持器的传递函数为h(s)()，则一阶保持器法离散公式为D(z)Z[h(s)D(s)]，故

Ts1Ts1eTs2(1z1)2Ts1(z1)2Ts1()D(s)]Z[2D(s)]Z[2D(s)]。 D(z)Z[TsTsTz2s 根据台劳近似式eTs21(z1)2Ts11(z1)1Z[2D(s)]中的Ts+1，则得到D(z)，用Ts取代D(z) Z[2D(s)]。21TsTzTzses2(z1)1(z1)2Ts1 显然，反过来也可以将D(z)Z[2D(s)]。 Z[2D(s)]近似为D(z)TzTz2ss 由此可以认为，“斜坡响应不变公式”和“串联一阶保持器公式”互为近似式。经验证，MATLAB的C2D函数中’foh’方法使用的就是“斜坡响应不变法公式”。

使用台劳近似式的条件是，T很小，时间单位也很小。

第四节 后向差分代换法(backward difference)

D(z)D(s)|sz1Tz 公式推导1— 近似微分： 设连续传递函数D(s) 将

U(s)1du(t)，则U(s)sE(s)，取反变换得e(t)， E(s)sdtu(kT)u((k1)T)du(t)e(kT)， e(t)中的微分用后向差分代替，得

dtTU(z)U(z)z1 对上式取z变换，得 E(z)，

T 整理后得离散传递函数 D(z)U(z)T， E(z)1z11z1z1 比较D(s)和D(z)，知 s。 TTz公式推导2— 近似积分：

tE(s)U(s)1，则U(s)t)e(t)dt，进行近似积分，而且把输出看作由两 设连续传递函数D(s)，反变换得 u(E(s)ss0部分组成，一部分是前一个时刻面积的累加值，一部分是本次的面积值。如果认为本次的面积是本次采样值与采样周期的乘积，则

u(kT)u((k1)T)Te(kT) 式(B) 取z变换，得 U(z)U(z)zTE(z) 合并同类项，得离散传递函数

1D(z)U(z)T E(z)1z11z1z1 比较D(s)和D(z)，知 s。 TTz 所以，后向差分法就是积分法，并且以当前采样值为计算依据，所以可以称为本次采样积分，式(B)的图形见下页。 公式推导3—有理化z与s的关系： zeTs|由台劳级数1Ts1z1，由此式求得s。 1TsTz以连续传递函数为基础的数字控制设计12/35

u((k-1)T)T*e(kT)后向差分u(kT)=u((k-1)T)+Te(kT)01T2T(k-1)TkT

% 差分近似法 function y=c2d1(ds,t0,str)

% 函数 c2d1 补充MATLAB控制类函数 c2d 之不足,可以将连续传递函数以向前差分法和向后差分法变换为离散传递函数 % 因为双线性变换法与向前差分法和向后差分法同属近似法,所以 c2d1 函数也包括了双线性变换法 % 输出参数 y ： LTI模型,离散传递函数 % 输入参数 ds ： LTI模型,连续传递函数 % t0 ： 数值,采样时间

% str ： 字符串,变换方法 backdiff=向后差分 fordiff=向前差分 tustin=双线性 syms z;

% ---------检查采样时间Ts是否大于0 if t0<=0

error('Ts<0 or Ts=0 unallowed'); end

% ---------根据输入的\"方法\确定代入公式 if strcmp(str,'backdiff') s=(z-1)/(t0*z);

string='coefficient of result of backward difference:'; string1='zero-pole of result of backward difference:'; elseif strcmp(str,'fordiff') s=(z-1)/t0;

string='coefficient of result of forward difference:'; string1='zero-pole of result of forward difference:'; elseif strcmp(str,'tustin') s=(2/t0)*((z-1)/(z+1));

string='coefficient of result of double linear:'; string1='zero-pole of result of double linear:'; else disp(' ');

error('please chek method parameter'); end

% ---------%以下将模型的传递函数转换为符号表达式,以便s=f(z)代入到ds中 [num,den]=tfdata(ds,'v'); %取模型连续传递函数的分子和分母(多项式) nums=poly2sym(num,'s'); dens=poly2sym(den,'s'); tran=nums/dens;

%生成一个符号分式

以连续传递函数为基础的数字控制设计13/35

%将多项式转换为符号

m=compose(tran,s); %代入 m=simplify(m); %化简

[n d]=numden(m);

%取符号表达式的分子和分母

n=sym2poly(n); %取符号表达式的分子多项式系数，按降幂排列 d=sym2poly(d); %取符号表达式的分母多项式系数，按降幂排列 % ---------%以下显示分子和分母的系数,是辅助性的 disp(string)

n=n/d(1); %使分母的最高项的系数为1 d=d/d(1); n_str=num2str(n); d_str=num2str(d);

n_disp=strcat('fn = ',n_str); d_disp=strcat('fd = ',d_str); disp(n_disp); disp(d_disp);

% ---------%以下生成离散传递函数 y=tf(n,d,t0);

% ---------%以下求出离散传递函数的零点和极点 disp(string1); [z p k]=zpkdata(y,'v')

% ---------%以下画出连续传递函数和离散传递函数的伯德图 bode(ds,'r',y,'g'); zpk(z,p,k,t0) % 程序结束

第五节 前向差分代换法(forward difference)

D(z)D(s)|z1sT 公式推导1— 近似微分： 设连续传递函数D(s) 将

du(t)U(s)1e(t)， ，则U(s)sE(s)，反变换得

E(s)sdtdu(t)u((k1)T)u(kT)e(t)中的微分用前向差分代替，得 e(kT)。 dtTu(kT)u((k1)T)e((k1)T)，

T 将所有的采样值提前一个周期，得

u(kT)u((k1)T)Te((k1T))。 取z变换，得 U(z)U(z)zTE(z)z1z)TE(z)z 故 U(z)(111，

1。

U(z)Tz1T 整理后得离散传递函数 D(z)。 E(z)1z1z1以连续传递函数为基础的数字控制设计14/35

比较D(s)和D(z)，知 sz1。 T 也可以不进行时间位移，直接从 u((k1)T)u(kT)Te(kT) 进行z变换，同样可得 U(z)zU(z)TE(z)， D(z) s公式推导2— 近似积分：

tE(s)U(s)1，则U(s) 设连续传递函数D(s)，它的积分表达式是u(t)e(t)dt，进行近似积分，而且把输出看作E(s)ss0U(z)T。 E(z)z1z1。 T由两部分组成，一部分是前一个时刻面积的累加值，一部分是本次的面积值。如果认为本次的面积是前次采样值与与采样周期的乘积，则 u(kT)u((k1)T)Te((k1)T) (式B)， 取z变换，得 U(z)U(z)zTE(z)z11，

U(z)Tz1 合并同类项，得离散传递函数 D(z)。 E(z)1z1 比较D(s)和D(z)，得 s1z1Tz1z1 T 所以，前向差分法就是积分法，并且以前一次采样值为计算依据，所以可以称为前次采样积分，式(B)的图形是

u((k-1)T)T*e((k-1)T)前向差分u(kT)=u((k-1)T)+Te((k-1)T)01T2T(k-1)TkT

公式推导3— 有理化z与s的关系：

s ze|由台劳级数1Ts，由此式求得

Tsz1。 T第六节 双线性近似法(tustin)

D(z)D(s)|2z1 sTz1 公式推导1— 近似微分： 设连续传递函数D(s) 将

U(s)1du(t)，则U(s)e(t)。 sE(s)，它的微分表达式是E(s)sdtdu(t)u(kT)u((k1)T)e((k1)T)e(kT)(t)用前采样值与本次采样值的平均值代替，得用后向差分代替，e，

dtT2U(z)U(z)z1E(z)z1E(z) 对上式取z变换 ， T21z) 即 U(z)(1T1E(z)(z1)， 2以连续传递函数为基础的数字控制设计15/35

故 D(z)U(z)T1z1， E(z)21z121z12z1 比较D(s)和D(z)，得 s。 T1z1Tz1公式推导2— 近似积分： 设连续传递函数D(s)tU(s)1，它的积分表达式是，u(t)e(t)dt，进行近似积分，而且把输出看作由两部分组E(s)s0成，一部分是前一个时刻面积的累加值，一部分是本次的面积值。如果认为本次面积是前次采样值与本次采样值，两次采样值的平均值平均值与采样周期的乘积，则u(kT)u((k1)T)T1e((k1)T)e(kT) (式B),

2 取z变换，得 U(z)U(z)zT 合并同类项，得

E(z)z1E(z),

2U(z)T1z1 D(z)。 E(z)21z121z12z1 比较D(s)和D(z)，得 s。 T1z1Tz1 所以，双线性近似法也是积分法，与差分法不同的是它是梯形积分，式(B)的图形是：

u((k-1)T)T*e((k-1)T)+e(kT)2双线性近似u(kT)=u((k-1)T)+T*0T1T2T(k-1)TkTe((k-1)T)+e(kT)2

公式推导3— 有理化z与s的关系： zeTsTse2Tse2|由台劳级数Ts2，由此式求得s2z1。 TsTz1121z12z1 故s。 z1Tz1公式推导4— 使用lnz的台劳级数： 由ze公式推导5— 使用th s压缩s平面：

Ts，得Tslnz2 脉冲响应不变法的主要缺点是频谱交叠产生的混淆，这是从S平面到Z平面的标准变换ze为了克服这一缺点，设想通过两个步骤建立S平面与Z平面一一对应的单值关系。 JJ1sT的多值对应关系导致的。jyT-xTs平面s1平面z平面 双线性变换的映射关系 第一步 将整个S平面压缩到S1平面的一条横带里 以连续传递函数为基础的数字控制设计16/35

考察双曲正切变换j平面1轴上的2jTth（1）0时，由0，也就是说将s平面的轴压缩到s1，当1由T2TT一段上。将这一关系推广到整个s平面，则得到s平面到s1平面的映射关系

TTs21T2j1T2j1T221es1T1T0jthj1，可得1。 ，系数的来源是，当时，sth（）2T2T2TT2T1es1T 第二步 通过标准变换关系将此横带变换到整个Z平面上去

Ts21zs21T 令ze，得S平面与Z平面的单值映射关系，s， 或者 z

TT1z11s211 验证：

⑴ 变换的单值性：当zejT1(式中1T)时，代入s21z121ej1T2j1Tth()j，即S的

T1z1T1ej1TT2虚轴映射到Z平面正好是单位圆。 S平面 Z平面 S1平面 T2T2TT(1)())j2222 ⑵ 变换的稳定性：sj 代入z表达式，得 z 和|z|。当0时，

TTTT22(1)j(1)()2222(1|z|<1即s左半平面映射在单位圆内，s右半平面映射在单位圆外，因此，稳定的模拟滤波器通过双线性变换后，所得到的数字滤波器也是稳定的。

第七节 预畸双线性法(prevarp)

第一步：求预修正频率 1*2t*tan（1） ，1为预防畸变的频率， 1为修正后的频率。 T2s*1 第二步：预畸变 D(s)D()。

s 第三步：变换 D(z)D(*)|。

1s2z1,*2tan(1T)Tz11T2 综合以上3步： D(z)D(s)|s1tan(1T2)z1 z1，1为预防畸变的频率 。 j1|D(z)||D(j)||D(e)|的原则进行增益匹配。 变换后，稳态增益将变化，需以|D(s)|或s0z11 举例： 若 ts1 若 ts210.3533z0.3533ts11,,'prewarp',)1， 则c2d(

z0.2934s11s20.2s1， 则c2d(ts21,,'prewarp',)10.212z20.424z0.212z20.9967z0.8448

第八节 零极点匹配法(matched)

D(S)  D(Z) 若 D(s)Ks(sz1)(sz2)...(szm)，

(sp1)(sp2)...(spn)nmKz(zez1t)(zez2t)...(zezmt)(z1) 则 D(z)。

(zep1t)(zep2t)...(zepnt)以连续传递函数为基础的数字控制设计17/35

分子上因子(Z+1)n-m的作用是：把s处的零点投射到来单位园上，使D(Z)的分母和分子的阶数就相同。（如果不加这个因子，D(Z)脉冲响应会产生n-m个采样时间的延迟，对系统造成不利影响?）。

nm1Kz(zez1t)(zez2t)...(zezmt)(z1) 或者 D(z)。

(zep1t)(zep2t)...(zepnt) 后一个公式比前一个公式少一个零点，MATLAB函数c2d的matched方法是按后一个公式进行变换的（但脉冲响应具有一个采样时间的延迟?）。 确定D(z)的增益kz的方法：

a 令D(s)与D(z)的稳态增益相等求得，即令|D(s)|s0|D(z)|z1， b 若D(s)分子有s因子，例如D(s)s，可依高频段增益相等原则确定增益,即|D(s)|s|D(z)|z1，也可选s1j1择某关键频率ω1处的幅频相等，即|D(j1)||D(e)|。

对零极点匹配法做的描述，可以使用符号数学编写m文件程序，“手工式”地进行验证。下面的程序以连续传递函数

1s31.8s21.8s1为例，运行此程序后，可以看到，文中描述的方法与c2d函数的mtched方法是一致的，也说明文本中

没有排版错误，这样应该更有于助理解和使用零极点匹配法。 %程序开始

syms s z kz t %预置采样周期 tt=0.7 %预置采样周期

s1=1 %建立ds的分子 s2=vpa(s^3+1.8*s^2+1.8*s+1) %建立ds的分母 s2r=solve(s2) %求ds的极点

z1=(z+1)^2 %建立dz分子的零点表达式,不包含增益

z1=vpa(expand(z1),5) %展开dz分子的零点表达式为多项式表达式,不包含增益 z2=(z-exp(t*s2r(1)))*(z-exp(t*s2r(2)))*(z-exp(t*s2r(3))) %建立dz分母的极点表达式

z2=vpa(subs(z2,'t',tt),5) %将dz分母中周期t用预置的值tt代替，不改变极点表达式形式 z2=vpa(expand(z2),5) %展开dz分母的极点表达式为多项式表达式

kz=limit(s1/s2,s,0)/limit(z1/z2,z,1) %计算dz的静态增益 dz=vpa(eval(kz*z1/z2),4) %获得dz的表达式

%display(strcat('Sample time: ',num2str(tt),' seconds')) %对dz的说明

display(['Sample time: ',num2str(tt),' seconds']) %对dz的说明,用方括号代替strcat,以保留空格 display('Discrete-time transfer function.') %续对dz的说明

ds0=tf(1,[1,1.8,1.8,1]) %比较用c2d展开的结果 dz0=c2d(ds0,tt,'matched') %程序结束

以连续传递函数为基础的数字控制设计18/35

第三章 时域化算法

控制系统中控制算法一般是指控制器的传递函数，但要用计算机实现，还应该将z域的离散传递函数，转换为时域的差分方程。由于差分方程已经清晰地显示出输入输出的时序关系，甚至可以在一定程度上把它等同于计算机程序流程图。 传递函数由z域向t域转换，也有几算法。如果把控制器的传递函数看作一体，叫做直接算法；如果把控制器的传递函数分解为几个传递函数的积，叫做串联算法；如果把控制器的传递函数分解为几个传递函数的和，则叫做并联算法；但串联算法和并联算法的z域分解式的转换，仍然要归结为直接算法。

通过不同形式的从z域到t域的数学转换，可以得到，同一个z传递函数的不同形式的输出差分表达式。对于直接算法，共可获得4种形式的输出差分表达式，而串联算法和并联算法形成的一阶分式和二阶分式的转换，同样归结为直接算法。 信流图对时域的输入输出关系的描述，比差方程形象，因此，对每一种算法都根据数学公式画出了它的信流图。参照信流图所表达的信号之间的关系，将给程序设计带来很大的方便。在信流图之后，给出了c语言的参考语句。

离散传递函数有z格式和z-1格式两种，在将离散传递函数变换为差分方程前，要将z格式的传递函数改写为z-1格式的传递函数。

远离当前时刻的过去时刻叫“后”： ......  m(-2) m(-1) m(0) 当前时刻和靠近当前时刻的过去时刻叫“前”： ......  m(-2) m(-1) m(0) 公式和信流图中，E(Z)，U(Z)，e(k)，u(k)表示输入量和输出量，Mi(Z)，Ni(Z)，mi(k)，ni(k)表示中间变量，ai，bi表示系 数。 程序中，输入量和输出量仍用e(k)，u(k)表示，输入输出的延时形式e(k-i)，u(k-i)和中间变量mi(k)，ni(k)，分别改为用 数组形式e(i)，u(i)，m(i)，n(i)表示，系数ai，bi也改用数组形式a(i)，b(i)表示。 信流图路径上的z-1表示，流经该路径的信息，在当前时刻（k时刻）到达指向的节点，将被延迟一个周期，在下一个时刻（k-1时刻）被引用。z-1本身没有数值意义。 给出的c语言程序参考语句，仅为说明算法，很不严谨，如：初始化语句未与主体语句隔离，主体语句未作无限循 环，系数没有列出数据表，当写成函数的形式时，可能需要设置静态变量，等等。 第一节 直接算法1 — 双中间变量向后递推

设离散传递函数为

U(Z)E(Z)1aiZi1i0nibiZn，则输出的Z变换表达式为，

iU(Z)biE(Z)ZiaiU(Z)Zi

i0i1nn b0E(Z)b1E(Z1)b2E(Z2)...bnE(Zn)a1U(Z1)a2U(Z2)...anU(Zn)。 反变换得输出u的时间表达式：

u(k)b0e(k)b1e(k1)b2e(k2)...bne(kn)a1u(k1)a2u(k2)...anu(kn)。 根据U(k)表达式，画出直接算法1的信流图：

e(k)e(k)z-1b0b1b2bn-1bn1-a1-a2-an-1-anu(k)u(k-1)u(k-2)z-1z-1e(k-1)z-1e(k-2)z-1e(k-n+1)z-1e(k-n)-1u(k-n+1)z-1u(k-n)z

特点：2n个中间变量，取得输入量e(k)后，要进行2n+1次乘法，n次加法，n次减法才能得出输出量u(k)。 for(i=1;i以连续传递函数为基础的数字控制设计19/35

u=b(0)*e; /* u(k)为本次输出值，e(k)为本次采样值 */ for(i=1;i1;i--) { e(i)=e(i-1); u(i)=u(i-1); } /* 数据延时*/ e(1)=e; u(1)=u;

第二节 直接算法2 — 双中间变量向前递推

设离散传递函数为

U(Z)M0(Z)U(Z)E(Z)E(Z)M0(Z)1aiZii1i0nibiZn11aiZii1ni0biZi，

n 中间变量 M0(Z)E(Z)(aiM0(Z)Zi。 Zi)，输出量 U(Z)biM(Z)i1nni0 第一步 以和积递推的形式展开M0(Z)中的和号∑：

iM0(Z)E(Z)(aiM0(Z)Z)E(Z)((...((anM0(Z)Zi1n1

an1M0(Z))ZMn1(Z)1...)Z1a2M0(Z))ZM2(Z)1a1M0(Z))ZM1(Z)1。

Mn(Z)11 上式中，令 Mn(Z)anM0(Z)，Mn1(Z)Mn(Z)Zan1M0(Z)，...，M1(Z)M2(Z)ZaM(Z)。 101 则中间变量 M0(Z)E(Z)M1(Z)Z，反变换得 m0(k)e(k)m1(k1)。

M1(Z)到Mn(Z)的反变换为：mn(k)anm0(k)，mn1(k)mn(k1)an1m0(k)，...，m1(k)mn(k1)a1m0(k)。 第二步 以和积递推的形式展开输出量U(Z)中的和号∑：

(...((bnM0(Z)Zbn1M0(Z))Z...)Z U(Z)b0M(Z)(Nn(Z)111b2M0(Z))ZbM(Z))Z10N2(Z)11，

Nn1(Z)1 上式中，令 Nn(Z)bnM0(Z)， Nn1(Z)Nn(Z)Zbn1M0(Z)，...，N1(Z)N2(Z)ZbM(Z)， 101 则输出量 U(Z)b0M0(Z)N1(Z)Z，反变换为 u(k)b0m0(k)n1(k1)。

1N1(Z)N1(Z)到Nn(Z)反变换为：nn(k)bnm0(k)，nn1(k)nn(k1)bn1m0(k)，...，n1(Z)n2(k1)b1m0(k)。 用m0(k)e(k)m1(k1)和u(k)b0m0(k)n1(k1)画出直接算法2的信流图：

e(k)m1(k-1)m0(k)z-1m1(k)z-1m2(k)z-1mn-1(k)z-11-a1-a2-an-1b0b1b2bn-1bnu(k)n1(k)n2(k)nn-1(k)nn(k)m2(k-1)mn-1(k-1)mn(k-1)mn(k)-an

每一个mi(k)和ni(k)，都是在被更新后的下一个时刻被引用，因此，信流图中的节点mi(k-1)和ni(k-1) 分别直接由mi(k)和ni(k)兼任，并用z-1因子表明信息的延迟。这样，直接算法2的信流图简化为下图：

e(k)m0(k)z-1m1(k)z-1m2(k)z-1mn-1(k)z-1mn(k)1-a1z-1-a2z-1-an-1z-1-anz-1b0b1z-1b2z-1bn-1z-1bnz-1u(k)n1(k)n2(k)z-1z-1nn-1(k)z-1nn(k)z-1

以连续传递函数为基础的数字控制设计20/35

for(i=1;i/* 求本次u(k) */

u=b(0)*m(0)+n(1)

for(i=1,ifor(i=1;i/* 求下一次的m(n) */

第三节 直接算法3 — 单中间变量向后递推

U(Z)M(Z)U(Z) 设离散传递函数为E(Z)E(Z)M(Z)i0nibiZn1aiZii111aiZii1ni0biZi，

n 由

M(Z)E(Z)11aiZi1ni 得中间变量 M(Z)E(Z)aiM(Z)Zi1ni，

直接写出中间变量的差分方程 m(k)e(k)aim(ki)e(k)a1m(k1)a2m(k2)...anm(kn)。

i1n 由

nU(Z)ni 得输出变量 U(Z)biM(Z)Zi, biZM(Z)i0i0n 直接写出输出变量的差分方程 u(k)bim(ki)b0m(k)b1m(k1)...bnm(kn)。

i0 画出直接算法3的信流图：

e(k)1z-1z-1m(k)m(k-1)m(k-2)z-1z-1b0b1b2bn-1bnu(k)-a1-a2-an-1-anz-1z-1-1m(k-n+1)zz-1m(k-n)

for(i=1;ifor(i=1,ifor(i=1;ifor(i=n;i>0;i--) m(i)=m(i-1); /* 中间变量m(i)延时 */

/* 求本次u(k) */

第四节 直接算法4 — 单中间变量向前递推(简约快速算法)

U(Z)E(Z)i0nibiZn 设离散传递函数为

1aiZii1， 则 U(Z)b0E(Z)(biE(Z)aiU(Z))Zi

i1n 以和积递推的形式展开和号∑：

((bnE(Z)anU(Z))Z U(Z)b0E(Z)((...1bn1E(Z)an1U(Z))Z1...)Z1b1E(Z)a1U(Z))Z1。

Mn(Z)以连续传递函数为基础的数字控制设计21/35

Mn1(Z)M1(Z)

1 令 Mn(Z)bnE(Z)anU(Z)，Mn1(Z)Mn(Z) Zbn1E(Z)an1U(Z)，...，M1(Z)M2(Z)ZbE(Z)aU(Z)，1111 则 U(Z)b0E(Z)M1(Z)Z，反变换得 u(k)b0e(k)m1(k1)。

M1(Z)到Mn(Z)的反变换为： m1(k)b1e(k)a1u(k)m2(k1)，m2(k)b2e(k)a2u(k)m3(k1)，... ， mn1(k)bn1e(k)an1u(k)mn(k1)，mn(k)bne(k)anu(k)。 用u(k)和mi(k)的表达式，画出直接算法4的信流图：

e(k)b0m1(k-1)b1b2m1(k)m2(k)zz-11-a1m2(k-1)-1u(k)-a2mn-1(k-1)-1bn-1bnmn-1(k)mn(k)zz-an-1mn(k-1)-1-an

和直接算法2一样，每一个mi(k)都是在被更新后的下一个时刻被引用，因此，信流图中的节点mi(k-1)可以直接由mi(k)兼任，并用z-1因子表明信息的延迟。这样，直接算法4的信流图简化为下图：

e(k)b0b1z-1b2z-1bn-1z-1bnz-1m1(k)m2(k)z-11-a1z-1-a2z-1-an-1z-1-anz-1z-1u(k)-1mn-1(k)z-1mn(k)z

特点：n个中间变量，在求出输出信号u(k)并从输出口输出后，还需要为下一次输出进行2n+2次乘法和2n+1次加法，显然，这一点并不影响输出的实时性。 for(i=1;ifor(i=1;i/* 初始化 */

第五节 串联算法

设离散传递函数为

U(Z)E(Z)1(aiZ)i1i0ni(biZ)n，

i 将传递函数改写成零极点模型，按照系数相近的原则对分子和分母配对，再使分子分母中的常数项为1，得传递函数的串联算法表达式如下：

11b21Z2U(Z)1b01Z11b0mZ11b1bb2nZ211Z1nZc...... 11E(Z)1a01Z11a0mZ11aa21Z21aa2nZ211Z1nZ 在串联算法表达式中，所有的因式只有一阶分式和二阶分式两种形式。

1U(Z)1b1Z 求一阶分式 的时域表达式： 1E(Z)1aZ1 因为 U(Z)E(Z)bE(Z)ZaU(Z)Z1111，

以连续传递函数为基础的数字控制设计22/35

所以 u(k)e(k)(b1e(k1)a1u(k1))。

从u(k)画出串联算法一阶分式的信流图,图中m1(k1)b1e(k1)a1u(k1)：

e(k)1b1z-1-1m1(k)z1-a1z-1u(k) m(1)=0; /* 初始化 */

/* 求下一次的 m(k-1) */

u(k)=e(k)+m(1); /* 求本次的u(k) */ m(1)=b(1)*e(k)-a(1)*u(k);

12U(Z)1b1Zb2Z 求二阶分式 的时域表达式： 12E(Z)1aZaZ12 U(Z)E(Z)bE(Z)Zb2E(Z)Z112aU(Z)Za2U(Z)Z1112

E(Z)((b2E(Z)a2U(Z))Z1b1E(Z)a1U(Z))Z 令 M2(Z)b2E(Z)a2U(Z)

M1(Z)M2(Z)ZbE(Z)aU(Z) 11 则 U(Z)E(Z)M1(Z)Z11

故 u(k)e(k)m1(k1)

且 m1(k)m2(k1)b1e(k)a1u(k) m2(k)b2e(k)a2u(k)

从u(k)画出串联算法二阶分式的信流图：

e(k)1b1z-1b2z-1-1m1(k)z-1m2(k)z1-a1z-1-a2z-1u(k) m(1)=0;m(2)=0; /* 初始化 */ u(k)=e(k)+m(1); /* 求本次的u(k) */ m(1)=b(1)*e(k)-a(1)*u(k)+m(2); m(2)=b(2)*e(k)-a(2)*u(k);

/* 求下一次的 m(1) */ /* 求下一次的 m(2) */

第六节 并联算法

设离散传递函数为

U(Z)E(Z)1(aiZ)i1i0ni(biZ)n，

i 通过部分分式的分解，将传递函数改写若干个部分分式的和，即

1U(Z)b01b0mbbn0bn1Z110b11Zc...... ， 1212E(Z)1a01Z11a0mZ11aZaZ1aZaZ1112n1n2 在并联算法表达式中，所有的项只有一阶分式和二阶分式两种形式。

以连续传递函数为基础的数字控制设计23/35

求并联算法一阶分式

U(Z)b0的时域表达式： 1E(Z)1aZ11 U(Z)b0E(Z)aU(Z)Z1

u(k)b0e(k)a1u(k1) 画出并联算法一阶分式的信流图：

e(k)b0m1(k)z-11-a1z-1u(k) m(1)=0; /* 初始化 */

u(k)=b(0)*e(k)+m(1); /* 求本次的u(k) */ m(1)=-a(1)*u(k); /* 求下一次的 m(1) */

1U(Z)b0b1Z 求并联算法二阶分式 的时域表达式： 12E(Z)1aZaZ12 U(Z)b0E(Z)bE(Z)ZaU(Z)Za2U(Z)Z11112 ，

b0E(Z)((a2U(Z))Z1bZ)aZ))Z1E(1U(11 令 M2(Z)a2U(Z)，M1(Z)M2(Z)ZbE(Z)aU(Z)， 111 则 U(Z)b0E(Z)M1(Z)Z，

故 u(k)b0e(k)m1(k1)，

由M1(Z),M2(Z)反变换得 m1(k)m2(k1)a1u(k)，m2(k)a2u(k)。 画出并联算法二阶分式的信流图：

e(k)b0b1z-1m1(k)z-1m2(k)z-11-a1z-1-a2z-1u(k) m(1)=0;m(2)=0; /* 初始化 */ u(k)=e(k)+m(1); /* 求本次的u(k) */ m(1)=b(1)*e(k)-a(1)*u(k)+m(2);

/* 求下一次的 m(1) */

m(2)=-a(2)*u(k); /* 求下一次的 m(2) */2

以连续传递函数为基础的数字控制设计24/35

第四章 数字PID控制算法

PID控制算法是模拟手动调节过程而构造出来的，它有3个特点： ① 控制算式中比例项，积分项，微分项物理意义明确。

② 控制算式中比例带，积分时间，微分时间等3个参数互相独立，因而可以在闭环中用凑试法确定这3个参数，并且不必预先知道控制对象的特性。

③ 将一阶对象或二阶对象对象与控制器形成闭环系统，假定闭环系统的传递函数为e传递函数。

由于有这3个特点，PID控制算法具有很强的普适性，因而被广泛地应用在“任意”场合。不过，若受控对象的纯滞后比较大，PID控制算法可能会出现积分饱和等问题，因此，PID控制算法出现了一些变形。

在工业生产过程中，例如对温度，压力，流量等参数进行控制，使用PID控制算法有比较成熟的经验。

，则控制器的传递函数总是PID

第一节 微分方程和差分方程

1 微分方程：

1tde(t) u(t)kP[e(t)t)dtTD]。 e(TdtI02 积分分离微分方程： u(t)kP[e(t)1tde(t)t)dtTD]。 e(TdtI03 位置式差分方程： u(kT)KP[e(kT)TkTe(jT)D(e(kT)e((k1)T)]。 TTj1I4 增量形式差分方程：

令 u(kT)u(kT)u((k1T))

T) 得 u(kT)KP[(e(kT)e((k1)输出变量： 输入变量： 参变量：

TTe(kT)D(e(kT)2e((k1)T)e((k2)T)]。 TITu(t) 纯量，对于实际的物理量，设阀门全开时为1。

e(t) 纯量，为误差值，对于实际的物理量，设e(t)的最大值为1。

kP 比例系数，纯量，忽略积分项和微分项，它使输出为1时输入量的变化，它是比例带的倒数。

TI 积分常数，时间量，在忽略比例项和微分项且kP1时，它是使输出为1时需要的积分时间。 TD 微分常数，时间量，为微分输入的作用时间。

第二节 不完全微分

由微分方程 u(t)kP[e(t)1tde(t)t)dtTD] ， e(TdtI0E(s)TDsE(s)] ， TIs1Ts，变为DE(s)，

Tfs1Tfs1得s传递函数 U(s)kP[E(s)将微分通道 TDsE(s)串联一个一阶环节 其中同TfTD，即Td的若干分之一。 f为TTk以连续传递函数为基础的数字控制设计25/35

所以不完全微分的传递函数为 U(s)kP[E(s) kPE(s)E(s)TsDE(s)] TIsTfs1kPkTsE(s)PDE(s) TIsTfs1 UP(s)UI(s)UD(s)。 经过离散化得 u(k)uP(k)uI(k)uD(k)。 以下分别计算各个通道：

比例通道： 令 UP(s)kPE(s)， 则 UP(z)kPE(z)， uP(k)kPe(k) 积分通道： 令 UI(s) 则 UI(z)1E(s)， TIskP1z1TITE(z)，

uI(kT)微分通道： 令 UD(s)Te(kT)uI((k1)T)， TITDsE(s)，

Tfs11z1kPTDTE(z)， 则 UD(z)11zT1fT 变形后得 UD(z)kPTDTfTE(z)kPTDTfTE(z)z1T1fUD(z)z TfT1uD((k1)T) TTf 化为时域表达 uD(kT) 上式中令 KDkPTDTfTe(kT)kPTDTfTe((k1)T)kPTDTT和f，1 TTTTTff1)[e(kT)e((k1)T)]uD((k1)T) 则 uD(kT)KD(k)u(k)u(k1)变换为增量式。 以上比例，积分，微分，各个通道都是位置式的，同样可以令u(第三节 参数选择

采样时间

受控物理量 温度 压力 流量 液位 15--20 3--10 1-- 5 6-- 8 采样时间(s) 取纯滞后时间常数 优先选用6--8s 优先选用1--2s 优先选用7s 以连续传递函数为基础的数字控制设计26/35

备 注 成份 15--20 优先选用18s 当系统中仅惯性时间常数起作用时，s10m。 当系统中纯滞后时间占有一定位置时，应选择T01.。 当系统中纯滞后时间占主导地位时，应选择T01.。 控制参数 受控物理量 温度 压力 流量 液位 成份 特点 对象为多容，滞后较大，须用微分 对象为多容，滞后一般较大，不用微分 对象时间常数小并有噪声，应较大，TI应较短，不用微分 在允许有余差时，不必用积分，不用微分 (%) 20-60 30-70 40-100 20-80 TI(分钟) 3.0-10 0.4-3 0.1-1 TD(分钟) 0.5-3 -- -- -- 第四节 c51框架

这仅仅是一个PID控制程序的基本架构，在实际应用时，还有很多细节问题需要处理。 #include #include

/*============================================================================================*/ typedef struct PID {

double SetPoint; // 设定目标 Desired Value

double Proportion;

// 比例常数 Proportional Const

double Integral; // 积分常数 Integral Const double Derivative;

double LastError; double PrevError; double SumError; } PID;

/*============================================================================================*/ double PIDCalc( PID *pp, double NextPoint ) //PID计算 {

double dError,Error;

Error = pp->SetPoint - NextPoint; pp->SumError += Error; // 积分 dError = pp->LastError - pp->PrevError; pp->PrevError = pp->LastError; pp->LastError = Error; return (pp->Proportion * Error

// 比例项

// 积分项 // 当前微分 // 偏差

// Error[-1] // Error[-2] // Sums of Errors

// 微分常数 Derivative Const

+ pp->Integral * pp->SumError

+ pp->Derivative * dError // 微分项 );

以连续传递函数为基础的数字控制设计27/35

}

/*===========================================================================================*/ void PIDInit (PID *pp) // 初始化PID 数据结构 {

memset ( pp,0,sizeof(PID)); }

/*============================================================================================*/ double sensor (void) // 读AD转换器，这里给出的是一个空函数 {

return 100.0; }

/*============================================================================================*/ void actuator(double rDelta) // 向DA转换器输出，也是给出一个空函数 { }

/*============================================================================================*/ void main(void) // 主函数 {

PID sPID; double rOut;

// 变量定义，PID 数据结构 // 变量定义，PID 响应值 (输出) // 这个返回值是虚设的

double rIn; // 变量定义，PID 反馈 (输入)

PIDInit ( &sPID ); // 初始化PID 数据结构 sPID.Proportion = 0.5; // 初始设计PID参数 sPID.Integral = 0.5; sPID.Derivative = 0.0; sPID.SetPoint = 100.0;

for (;;) {

// 进行PID 处理

// 设定“设定点”

//在这里写检查是否改变PID参数的语句和改变PID参数的语句 rIn = sensor (); // 读反馈值 rOut = PIDCalc ( &sPID,rIn ); actuator ( rOut ); } }

// PID 迭代

// 系统输出

以连续传递函数为基础的数字控制设计28/35

第五章 保持器

第一节 零阶保持器 特

点 零阶保持器是一个连续环节，它对输入的脉冲在一个周期内按恒值外推后输出。 e(t) t零阶保持器 脉冲响应： 设在某一k=0时刻e(0T)k0，则零阶保持器单位脉冲输入的响应为 t)k0 h0(0tT 其余时刻 h1(t)0t)k0[u(t)u(tT)]，u(t)传递函数： 用延时函数将零保持器的脉冲响应函数表示为两个阶跃函数的代数和 h(式中，

tT)是纯滞后为T的单位阶跃函数。 是单位阶跃函数，u(1eTs h(t)的拉氏变换为H(s)k0，这是零阶保持器的输出值传递函数。因为脉冲输入e0=e(t)(kT)，

sH(s)1eTst)k0，并且，在时刻kT时e(所以脉冲输入的拉氏变换为E0(s)k0，所以零阶保持器的传递函数Ha(s)。

E0(s)s第二节 一阶保持器

特 点： 一阶保持器是一个连续环节，它对输入的脉冲在一个周期内按线性外推后输出，图中，采样周期的斜率是k000k0，如图中虚线（-T,k0）所示，退出采样周期后的斜率是，如图中（k0,T）所示。 TTe(t) t一阶保持器 脉冲响应式： 设在某一k=0时刻e(0T)k0，则零阶保持器单位脉冲输入的响应为 t)k0(1 h0((t)k0(0 h1t)0 h2(10tt)k0(1)TT01(tT))k0(1t)TT0tT Tt2T 其余时刻

传递函数： 为了求得一阶保持器的传递函数，将脉冲响应图形表达为6个波形的代数和：

t)k0(u(t) h(ttTt2Tu(t)2u(tT)2u(tT)u(t2T)u(t2T))。 TTT以连续传递函数为基础的数字控制设计29/35

式中，u(t)是单位阶跃函数，u(tT)，u(t2T)分别是纯滞后为T和2T的单位阶跃函数，此式的拉氏变112eTs2eTse2Tse2Ts换为 k0()。

sTs2ssTs2Ts2Ts11eTs2()，这是一阶保持器输出值的输出表达式，因为脉冲输入 将上式整理后得H(s)k0Tse0=e(t)(kT)，而脉冲输入的拉氏变换为E0(s)=k0，所以一阶保持器的传递函数为Ha(s) --------------------------------------------------- 附录

H(s)Ts11eTs2()。 E0(s)Ts两种一阶离散化方法的结果的比较

(在MATLAB v6.5上运行)

%第一套\"MATLAB一阶离散化\"方法用到的公式

1(1eTs)s2TeTs1(1z1)s2Tz1

22Ts1(1eTs)1(1z1)%第二套\"标准一阶离散化\"方法用到的两部分 2TTs2s%两种\"一阶离散化\"方法的对比 Clc %清除窗口 Clear %清除工作区

s0=tf(1,[1,1,1]); %原传递函数 t=0.05; z=tf('z',t);

s11=tf(1,[1,0,0]); %第一套\"一阶离散化\"方法的两部分 z11=tf([1,-2,1],[t,0],t);

s12=tf([t,1],[1,0,0]);

%第二套\"一阶离散化\"方法的两部分

%设定

z12=tf([1,-2,1],[t,0,0],t);

z0=c2d(s0,t,'foh');

%与第一套\"一阶离散化\"对应的MATLAB的\"foh\"方法

%进行第一套\"一阶离散化\" %进行第二套\"一阶离散化\"

z1=(c2d(s0*s11,t,'imp'))*z11; z2=(c2d(s0*s12,t,'imp'))*z12;

zpk0=zpk(z0); zpk1=zpk(z1); zpk2=zpk(z2);

%以下提取零极点和增益 z_x(1)=z1; z_x(2)=z2; for r=1:2

[zx,px,kx]=zpkdata(z_x(r),'v'); x=length(zx);

%观察分子分母上的公因子

以连续传递函数为基础的数字控制设计30/35

y=length(px);

%以下z约去分子和分母中的公因数 for i=x:-1:1 for j=y:-1:1

if abs(zx(i)-px(j))<0.0001 %这一个数值不能太小 zx(i)=[]; px(j)=[]; y=y-1; break; end end end

%得到最简的零极点形式的传递函数 zpkn(r)=zpk(zx,px,kx,t); end

disp('************************************ zpk0') zpk(zpk0)

disp('************************************ zpkn(1)') zpk(z_x(1)) zpkn(1)

disp('************************************ zpkn(2)') zpk(z_x(2)) zpkn(2)

%以下绘图部分,仅供参考------------------------------- figure figure figure figure(1) step(z0) figure(2) step(zpkn(1)) figure(3) step(zpkn(2)) ◆ 下面的包对象中包含，本文使用的，用autoCAD制作的图形：离散化图，信流图，比较图等。用“好压”压缩软件压缩。 离散化图和信流图.zip 以连续传递函数为基础的数字控制设计31/35