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

实用文档/35第一章模拟化设计基础数字控制系统的设计有两条道路，一是模拟化设计，一是直接数字设计。如果已经有成熟的模拟控制器，可以节省很多时间和部分试验费用，只要将模拟控制器离散化即可投入应用。如果模拟控制器还不存在，可以利用已有的模拟系统的设计经验，先设计出模拟控制器，再进行离散化。将模拟控制器离散化，如果用手工进行，计算量比较大。借助数学软件MATLAB控制工具箱，可以轻松地完成所需要的全部计算步骤。如果需要的话，还可以使用MATLAB的SIMULINK工具箱，进行模拟仿真。第一节步骤步骤1模拟控制器的处理在数字控制系统中，总是有传输特性为零阶保持器的数模转换器（DAC），因此，如果模拟控制器尚未设计，则应以下图的方式设计模拟控制器，即在对象前面加上一个零阶保持器，形成一个新对象，然后针对这个新对象求模拟控制器D(s)。事实上，模拟控制器一般是已经设计好的，无法或不方便更改了，离散化后的系统只好作为近似设计了。然而，按照上述思路，可否将已有的控制器除以一个零阶保持器再离散化呢？还没有这方面的实际经验。以下假设选定的G(s)，D(s)如下图，而且不对G(s)作添加保持器的预处理。步骤2离散化模拟控制器离散化模拟控制器之前，先要确定离散化算法和采样时间。离散化算法有好几种，第二章中有详细的论述，现假定采用双线性变换法。确定采样时间，需要考虑被控对象的特性，计算机的性能，以及干扰信号的影响等，初步可按采样时间T<0.1Tp，Tp为被控对象时间常数，或T=(0.125～0.25)τ，为被控对象的纯滞后，初步确定后再综合平衡其它因素，当然这需要一定的经验，现在假定取0.05秒。假设模拟控制器为，在MATLAB中，用c2d函数进行离散化，过程为：ds=zpk(-2,-15,8)ds=zpk(-2,-15,8)%建立模拟控制器的s传递函数dz=c2d(ds,0.05,'tustin')%将模拟控制器按tustin方法转换为z传递函数的数字控制器dz=c2d(ds,0.05,'tustin')%将模拟控制器按tustin方法转换为z传递函数的数字控制器转换结果为：步骤3检验数字控制器的性能数字控制器的性能项目比较多，我们仅以直流增益，频率特性，零极点分布说明。直流增益dcgain(dz) 返回直流增益1.0667频率特性bode(ds,'r',dz,'g') 伯德图，见下页左图零极点分布pzmap(dz) 零极点分布图，见下页右图步骤4离散化控制对象为了进行模拟仿真，需要对控制对象进行离散化，由于步骤1所说的原因，应把被控对象视为零阶保持器与原对象的串连，即应对进行离散化，这时可在c2d函数中使用零阶保持器(zoh)方法，如果认为不需要添加零阶保持器，即直接对G(s)离散化，则应在c2d函数中使用冲击响应不变法（imp）。借用零阶保持器(zoh)方法，将对象带一阶保持器离散化的过程如下：%模拟控制器D(s)转换为D(z)的过程见前gs%模拟控制器D(s)转换为D(z)的过程见前gs=zpk([],[0,-2],20)%建立对象的s传递函数g1z=c2d(gs,0.05,'zoh')%借用c2d函数进行带零阶保持器的对象的离散化转换结果为：步骤5模拟仿真求离散系统的闭环传递函数和连续系统的闭环传递函数。离散系统的闭环传递函数为：连续系统的闭环传递函数为：用MATLAB算TRCZ与TRCS：trcz=dztrcz=dz*g1z/(1+dz*g1z)trcs=ds*gs/(1+ds*gs)结果为：用MATLAB函数STEP画阶跃响应图形：holdonholdon%图形保持step(trcs,'r',2)%画连续系统的阶跃响应图,红色,终止时间为2秒step(trcz,'b',2)%画离散系统的阶跃响应图,兰色,终止时间为2秒响应图形为：步骤6求数字控制器的时域表达式上面已经求出，连续传递函数的tustin离散式为，或。对上式取z反变换，得时域表达式，根据此式就可以写出计算的程序代码来了。除上述步骤之外，在编写程序代码时，还需要考虑几个问题：①ADC位数。②数据类型在控制算法中有3种数据类型可以采用：无符号整数，定点数，浮点数。我们知道，无符号整数运算既是最简单和速度最快的运算，又是定点数运算和浮点数运算的基础。在数据动态范围比较小的情况下，应尽可能用无符号整数运算代替定点数运算和浮点数运算。③数值计算误差数值计算引入的误差有3个方面，一是定点数字长不够或者浮点数有效数字过少，一是两个相近的数相减，一是加减乘除次数过多。在程序设计中，应优化算法以避免计算误差的引入。

第二节在MATLAB中离散化1.建立s降幂传递函数①建立多项式型s降幂传递函数方法1.sys=tf(num,den)方法2.s=tf('s')，再令sys=f(s)若G0以零极点形式给出，即已知，仍可用tf函数建立多项式型s降幂传递函数，但需用多项式乘法函数conv配合，G0=tf(conv([3],[0.51]),conv(conv([10],[11]),[0.251]))②建立零极点型s传递函数方法1.sys=zpk(z,p,k)方法2.s=zpk('s')，再令sys=f(s)。2.传递函数的转换①将任意形式的s传递函数转换为s降幂传递函数 sys=tf(sys)②将任意形式的s传递函数转换为s零极点传递函数 sys=zpk(sys)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)当已有多项式形式的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的分项分式，[anadak]=residue(num,den)。分项结果可能出现共轭复数，在这种情况下应将含有共轭复数的分式合并成二次有理质分式。④根据分项结果手工写出z降幂多项式形式的并联表达式。⑤“手工”将z降幂多项式形式的并联表达式改写成z-1降幂多项式形式的并联表达式。⑥“手工”对z-1降幂多项式形式的并联表达式中的每一个分式项降阶，即将每一个分式项变形，使得分式项的分子的阶次比分母的阶次低1阶，变形完毕再将全部常数项合并。举例设有z-1降幂传递函数改写成z降幂传递函数，取分项矢量[anadak]=residue(nz,dz)得an=[0.5101-0.1971]ad=[1.00000.2700]ak=0.1手工写分项分式，令，令又令故验证实际上，大多数数字控制器的传递函数都是一阶或者二阶的，所以需要分解的并不是很多。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-，可按如下方法处理：1.在tf函数中使用属性’inputdely’或者’iodely’，例如： >>tf([1-1],[135],'inputdelay',0.35) 将返回以下形式的传递函数s–1 exp(-0.35*s)*s^2+3s+5 使用这个方法不能建立形如的传递函数，因为带延时的传递函数不能与不带延时的传递函数相加，但可以使用c2d进行离散化，但要求延时时间t必须是采样时间t的整数倍，若不是整数倍，则在转换时不理会延时环节，例如： 用tf函数建立2个传递函数，主体部分相同，但一个无输入延时，一个有输入延时0.35s，>>a=tf([1-1],[145])s-1s^2+4s+5>>a1=tf([1-1],[145],'iodelay',0.35)s-1exp(-0.35*s)*s^2+4s+5若采样时间为0.05，因为延时时间是采样时间的整数倍，转换结果的主体部分完全一样：>>c2d(a,0.05,'imp')z^2-1.039z+9.146e-018samplingtime：0.05z^2-1.807z+0.8187>>c2d(a1,0.05,'imp')z^2-1.039z+9.146e-018z^(-7)*samplingtime：0.05z^2-1.807z+0.8187若采样时间为0.1，因为延时时间不是采样时间的整数倍，结果的主体部分不一样：>>c2d(a,0.1,'imp')z^2-1.06z+4.349e-018samplingtime：0.1z^2-1.629z+0.6703>>c2d(a1,0.1,'imp')0.768z-0.851z^(-3)*samplingtime：0.1z^2-1.629z+0.67032.将e-Ts有理化设，因e-Ts的一阶有理表达式是，故。为了对ts进行离散化，首先使用tf函数建立lti模型的ts。若，则。3.在离散化时使用恒等式设，因，采样时间为Δt，若T=，则，因为离散化时总是认为，故取。根据以上假设，用MATLAB的c2d函数对g(s)进行离散化，则，进而。则。第四节控制函数分类以下函数在controltoolbox中，这里所述仅限于siso模型1创建多项式形式的传递函数sys=tf(num,den,t)创建一个z降幂离散传递函数sys，t是采样时间，行矢量num和den同上。sys=tf(sys)把一个任意的lti模型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创建任意形式的传递函数s=tf('s')指定多项式传递函数变量为s变量z=tf('z',t)指定多项式传递函数变量为z变量，t为采样时间形如的传递函数，既不能直接用tf函数建立，也不能直接用zpk函数建立。定义了s=tf('s')后，写出赋值式f=10/(s*(0.25*s+1)*(0.05*s+1))，回车后即可得到多项式型传递函数，再令f1=zpk(f)，可得。s=zpk('s')指定零极点传递函数变量为s变量z=zpk('z',t)指定零极点传递函数变量为z变量，t为采样时间定义了s=zpk('s')，写出赋值式f=10/(s*(0.25*s+1)*(0.05*s+1))，回车后即可得到零极点型传递函数。4建立z-1降幂离散传递函数tz-1tsys=filt(m)返回增益离散传递函数。5连续函数离散化sysd=c2d(sysc,t,method)把连续传递函数sysc转换成采样时间为t的离散传递函数，字符串method为离散化方法：'zoh'(零阶保持，即阶跃响应不变)，'foh'(一阶保持)，'imp'(冲击响应不变,v6以上版本)，'tustin'(双线性近似)，'prewarp'(带预畸变的双线性近似)，'matched'(零极点匹配)。注：①缺省的方法是'zoh'②'foh'的算法是而不是，这一点可以通过验证证实，验证方法是，令，但此结果中，的分子和分母将含有公因式，所以应进一步用zpk函数把表示成零极点形式，然后用“手工”的方法写出不含公因式的来，可以看到，最后的结果与用foh方法得到的结果完全一致。t。6取多项式模型传递函数的分子和分母的系数矢量[num,den]=tfdata(sys,'v')对于siso模型sys返回作为分子和分母系数的单行矩阵num和den。[num,den,t]=tfdata(sys,'v')同上，同时返回采样时间t。7取零极点模型传递函数的零点和极点的单行矩阵[z,p,k]=zpkdata(sys,'v')返回lti模型sys的零极点矢量z和p，增益k。8取e-近似式[num,den]=pade(t,n)返回e-的n阶pade近似式，行矢量num和den是s的降幂多项式系数。e-Ts的一阶有理表达式是。e-的高阶有理表达式是9画阶跃响应图step(sys)画出由tf，zpk，orss等函数创建的lti模型sys的阶跃响应图。step(sys,tfinal)画出lti模型sys从t=0到t=tfinal的阶跃响应图。对于未指定采样时间的离散模型，tfinal被解释为采样的数目。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就是一个尺寸为[ltnynu]的阵列，而y(:,:,j)给出第j个输入通道的阶跃响应。10画脉冲响应图impulse脉冲响应函数，用法与step相同11画频率响应图--伯德图(连续或离散)bode(sys)画伯德图bode(sys,{wmin,wmax})在频率wmin，wmax(inradians/second)之间画伯德图bode(sys,w)按指定的频率矢量w(inradians/second)画伯德图bode(sys1,sys2,...)画多个lti模型sys1，sys2，...的伯德图bode(sys1,sys2,...,w)按指定的频率矢量w(inradians/second)画多个lti模型sys1，sys2，...的伯德图以下函数在符号工具箱symbolick中，需注意，在使用这些函数前，要对所使用的变量进行符号说明，例如：symsat%是用空格分隔而不是用逗号分隔a=sin(t)symsat%是用空格分隔而不是用逗号分隔a=sin(t)L=laplace(a)12福里哀变换 fourier反福里哀变换 ifourier13拉普拉斯变换 laplace反拉普拉斯变换 ilaplace14z变换 ztrans反z变换 iztrans注：以上3种变换必须是符号表达式，例如：symstlaplace(sin(t))15改善公式的可读性 pretty16多项式转换为符号表达式 poly2sym17化简符号表达式 simplify18取符号表达式的分子和分母 numden以下函数在符号工具箱polyfun中，在公式变换中可能会用到19部分分式展开 [anadak]=residue(n,d)n和d分别为原分式的分子和分母矢量，an和ad分别为分项式的分子和分母矢量，ak为整式部分。这是一个数学公式，在数字控制器程序设计中，利用部分分式展开的方法，把高于2阶的分式变换为不高于2阶的分式之和，从而把高阶传递函数算法变为低阶传递函数并联的算法。在使用这个方法时，不论在分子矢量中还是在分母矢量中，如果有共轭复数出现，则应将其整合为2阶质因式。20多项式乘法 c=conv(a,b)a×b=c21多项式除法 [q,r]=deconv(b,a)b÷a=q...r，即b=conv(a,q)+r

第二章离散化算法连续传递函数离散化的核心环节，就是将控制器的s传递函数转换为z传递函数。离散化后，系统应该仍有好的稳定性，好的控制精度，而不是要求转换前后，两个数学公式等值。因此，离散化方法有多种。本章对这些方法的转换原理和由来进行了演绎。离散化后得到的离散传递函数的稳定性，没有进行讨论，仅列了一张表进行比较。离散的结果，虽然可能会控制精度降低，应该认为，这不是主要问题。真正影响控制精度的因素，主要还是采样周期的长短。一般情况下，由连续到离散的设计最好多实验几种方法（通过仿真，得出满意的结果）。因为匹配零、极点映射法、双线性变换法都能得出比较满意的结果，初步设计时，可以试用这些方法。而其实，后向差分近似法也是合理的选择。但MATLAB的c2d函数中没有这一方法，在该方法的介绍之后，给出了一个可由MATLAB引用的m文件函数。摘要3种保持器法 保持器法即将s函数串联上一个保持器后取z变换，也可以从“对输入信号的响应不变”的角度导出。无保持器法 （阶跃响应不变） imp零阶保持器法 （阶跃响应不变） zoh一阶保持器法 （斜坡输入不变） foh （一阶保持器） MATLAB中无此方法4种近似法 近似法将s与z的无理关系近似地化为有理关系，主要应用于传递函数从连续到离散的变换。后向差分近似法 MATLAB的C2D函数中没有这个方法前向差分近似法 MATLAB的C2D函数中没有这个方法双线性近似法 tustin预畸双线性近似法 ，1为进行预畸变的频率 prevarp增益匹配：或1种匹配法 匹配法完全从控制学的角度看问题，也是应用于传递函数从连续到离散的变换。零极点匹配法若 ， matched则 ，或者 ，确定增益kz：a令，b若D(s)分子有s因子，例如，可以令，也可以令比较离散化方法变换公式映射关系特点冲击响应不变法Z变换法imp脉冲响应采样值相同；容易产生频率混迭现象，阶跃响应不变法零阶保持器法zoh阶跃响应采样值相同；稳定增益不变。斜坡响应不变法一阶保持器法foh向前差分法D(s)稳定，D(z)可能不稳定；等效精度差。向后差分法变换计算简单；②如果D(s)稳定，D(z)稳定；③离散滤波器的过程特性及频率特性有一定的失真，需要较小的采样周期T。双线性变换法tustinD(s)稳定，D(z)也稳定；低频特失真，但无频率混迭现象。稳定增益不变；具有串联特性。预修正双线性变换法prevarp有前一种变换的特点;还能保证在关键频率ω1处，幅频特性不变。零极点匹配法matchedZ域与s域零极点位置一一对应；当没有零点时，补充z=1的零点可避免频率混迭现象。

第一节冲击响应不变法(imp，无保持器直接z变换法)公式推导1—无保持器，直接转换：公式推导2—按冲击响应不变的原则：（冲击响应不变，中的1表示）第二节阶跃响应不变法(zoh，零阶保持器z变换法)，式中T为采样周期，所以公式推导1—按串联零阶保持器的原则：因为零阶保持器的s传递函数为，故公式推导2—按阶跃响应不变的原则：设阶跃信号，则，，按照下一节对斜坡响应不变法的推导，立即可以写出第三节斜坡响应不变法(foh，一阶保持器z变换法)，，从一阶保持器的思路推出以上2个公式中T为采样周期，两式差异请参看“附录两种一阶离散化方法的结果的比较”。公式推导1—按斜坡响应不变的原则：设连续滤波器为D(s)，输入为斜坡函数e(t)=t，则，采样输出为us(kT)。按z变换的定义，us(kT)的变换就是的z变换，故。又设离散滤波器的传递函数为D(z)，它的输出为uz(KT)，按斜坡响应不变的要求，应有uz(kT)=us(kT)，故，但离散滤波器的输出为U(z)=E(z)D(z)，而斜坡函数e(t)=t的z变换为，求Uz(z)与E(z)的比值，得斜坡响应不变法离散公式。公式推导2—按串联一阶保持器的原则：一阶保持器的传递函数为，则一阶保持器法离散公式为，故。根据台劳近似式，用取代中的Ts+1，则得到。显然，反过来也可以将近似为。由此可以认为，“斜坡响应不变公式”和“串联一阶保持器公式”互为近似式。经验证，MATLAB的C2D函数中’foh’方法使用的就是“斜坡响应不变法公式”。使用台劳近似式的条件是，T很小，时间单位也很小。第四节后向差分代换法(backwarddifference)公式推导1—近似微分：设连续传递函数，则，取反变换得，将中的微分用后向差分代替，得，对上式取z变换，得，整理后得离散传递函数，比较和，知。公式推导2—近似积分：设连续传递函数，则，反变换得，进行近似积分，而且把输出看作由两部分组成，一部分是前一个时刻面积的累加值，一部分是本次的面积值。如果认为本次的面积是本次采样值与采样周期的乘积，则式(B)取z变换，得合并同类项，得离散传递函数比较和，知。公式推导3—有理化z与s的关系：，由此式求得。%差分近似法functiony=c2d1(ds,t0,str)%函数c2d1补充MATLAB控制类函数c2d之不足,可以将连续传递函数以向前差分法和向后差分法变换为离散传递函数%因为双线性变换法与向前差分法和向后差分法同属近似法,所以c2d1函数也包括了双线性变换法%输出参数y：LTI模型,离散传递函数%输入参数ds：LTI模型,连续传递函数%t0：数值,采样时间%str：字符串,变换方法backdiff=向后差分fordiff=向前差分tustin=双线性symsz;%检查采样时间Ts是否大于0ift0<=0error('Ts<0orTs=0unallowed');end%根据输入的"方法",确定代入公式ifstrcmp(str,'backdiff')s=(z-1)/(t0*z);string='coefficientofresultofbackwarddifference:';string1='zero-poleofresultofbackwarddifference:';elseifstrcmp(str,'fordiff')s=(z-1)/t0;string='coefficientofresultofforwarddifference:';string1='zero-poleofresultofforwarddifference:';elseifstrcmp(str,'tustin')s=(2/t0)*((z-1)/(z+1));string='coefficientofresultofdoublelinear:';string1='zero-poleofresultofdoublelinear:';elsedisp('');error('pleasechekmethodparameter');end%%以下将模型的传递函数转换为符号表达式,以便s=f(z)代入到ds中[num,den]=tfdata(ds,'v'); %取模型连续传递函数的分子和分母(多项式)nums=poly2sym(num,'s'); %将多项式转换为符号dens=poly2sym(den,'s');tran=nums/dens; %生成一个符号分式m=compose(tran,s); %代入m=simplify(m); %化简[nd]=numden(m); %取符号表达式的分子和分母n=sym2poly(n); %取符号表达式的分子多项式系数，按降幂排列d=sym2poly(d); %取符号表达式的分母多项式系数，按降幂排列%%以下显示分子和分母的系数,是辅助性的disp(string)n=n/d(1); %使分母的最高项的系数为1d=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);[zpk]=zpkdata(y,'v')%%以下画出连续传递函数和离散传递函数的伯德图bode(ds,'r',y,'g');zpk(z,p,k,t0)%程序结束第五节前向差分代换法(forwarddifference)公式推导1—近似微分：设连续传递函数，则，反变换得，将中的微分用前向差分代替，得。将所有的采样值提前一个周期，得，。取z变换，得，故。整理后得离散传递函数。比较和，知。也可以不进行时间位移，直接从进行z变换，同样可得，。。公式推导2—近似积分：连续传递函数则，(式B)，取z变换，得，合并同类项，得离散传递函数。比较和，得所以，前向差分法就是积分法，并且以前一次采样值为计算依据，所以可以称为前次采样积分，式(B)的图形是公式推导3—有理化z与s的关系：，由此式求得。第六节双线性近似法(tustin)公式推导1—近似微分：连续传递函数，则，它的微分表达式是。对上式取z变换，即，故，比较和，得。公式推导2—近似积分：连续传递函数，它的积分表达式是，，进行近似积分，而且把输出看作由两部分组成，一部分是前一个时刻面积的累加值，一部分是本次的面积值。如果认为本次面积是前次采样值与本次采样值，两次采样值的平均值平均值与采样周期的乘积，则(式B),取z变换，得,合并同类项，得。比较和，得。所以，双线性近似法也是积分法，与差分法不同的是它是梯形积分，式(B)的图形是：公式推导3—有理化z与s的关系：，由此式求得。公式推导4—使用的台劳级数：由，得故。公式推导5—使用ths压缩s平面：脉冲响应不变法的主要缺点是频谱交叠产生的混淆，这是从S平面到Z平面的标准变换的多值对应关系导致的。为了克服这一缺点，设想通过两个步骤建立S平面与Z平面一一对应的单值关系。双线性变换的映射关系第一步将整个S平面压缩到S1平面的一条横带里，当由时，由，。将这一关系推广到整个s平面，则得到s平面到平面的映射关系。第二步通过标准变换关系将此横带变换到整个Z平面上去令，得S平面与Z平面的单值映射关系，，或者验证：⑴变换的单值性：S平面Z平面S1平面⑵变换的稳定性：代入z表达式，得和。当时，|z|<1即s左半平面映射在单位圆内，s右半平面映射在单位圆外，因此，稳定的模拟滤波器通过双线性变换后，所得到的数字滤波器也是稳定的。第七节预畸双线性法(prevarp)第一步：求预修正频率，1为预防畸变的频率，为修正后的频率。第二步：预畸变。第三步：变换。综合以上3步：，1为预防畸变的频率。变换后，稳态增益将变化，需以或的原则进行增益匹配。举例：若，则若，则第八节零极点匹配法(matched)D(S)D(Z)若，则。分子上因子(Z+1)n-m的作用是：把处的零点投射到来单位园上，使D(Z)的分母和分子的阶数就相同。（如果不加这个因子，D(Z)脉冲响应会产生n-m个采样时间的延迟，对系统造成不利影响?）。或者。后一个公式比前一个公式少一个零点，MATLAB函数c2d的matched方法是按后一个公式进行变换的（但脉冲响应具有一个采样时间的延迟?）。确定D(z)的增益kz的方法：a令与的稳态增益相等求得，即令，b若D(s)分子有s因子，例如，可依高频段增益相等原则确定增益,即，也可选择某关键频率处的幅频相等，即。对零极点匹配法做的描述，可以使用符号数学编写m文件程序，“手工式”地进行验证。下面的程序以连续传递函数为例，运行此程序后，可以看到，文中描述的方法与c2d函数的mtched方法是一致的，也说明文本中没有排版错误，这样应该更有于助理解和使用零极点匹配法。%程序开始symsszkzt%预置采样周期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('Sampletime:',num2str(tt),'seconds'))%对dz的说明display(['Sampletime:',num2str(tt),'seconds'])%对dz的说明,用方括号代替strcat,以保留空格display('Discrete-timetransferfunction.')%续对dz的说明ds0=tf(1,[1,1.8,1.8,1])%比较用c2d展开的结果dz0=c2d(ds0,tt,'matched')%程序结束

第三章时域化算法控制系统中控制算法一般是指控制器的传递函数，但要用计算机实现，还应该将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—双中间变量向后递推设离散传递函数为，则输出的Z变换表达式为，。反变换得输出u的时间表达式：。根据U(k)表达式，画出直接算法1的信流图：特点：2n个中间变量，取得输入量e(k)后，要进行2n+1次乘法，n次加法，n次减法才能得出输出量u(k)。for(i=1;i<n+1;i++){e(i)=0;u(i)=0;} /*初始化*/u=b(0)*e; /*u(k)为本次输出值，e(k)为本次采样值*/for(i=1;i<n+1;i++)u=u+b(i)*e(i)+(-a(i))*u(i);for(i=n;i>1;i--){e(i)=e(i-1);u(i)=u(i-1);} /*数据延时*/e(1)=e;u(1)=u;第二节直接算法2—双中间变量向前递推设离散传递函数为，中间变量，输出量。第一步以和积递推的形式展开M0(Z)中的和号∑：上式中，令，，...，。则中间变量，反变换得。M1(Z)到Mn(Z)的反变换为：，，...，。第二步以和积递推的形式展开输出量U(Z)中的和号∑：，上式中，令，，...，，则输出量，反变换为。N1(Z)到Nn(Z)反变换为：，，...，。用和画出直接算法2的信流图：每一个mi(k)和ni(k)，都是在被更新后的下一个时刻被引用，因此，信流图中的节点mi(k-1)和ni(k-1)分别直接由mi(k)和ni(k)兼任，并用z-1因子表明信息的延迟。这样，直接算法2的信流图简化为下图：for(i=1;i<n+1;i++){m(i)=0;n(i)=0;} /*初始化*/m(0)=e+m(1); /*求本次u(k)*/u=b(0)*m(0)+n(1)for(i=1,i<n,i++)m(i)=(-a(i)*m(0))+m(i+1); /*求下一次的m(n)*/m(n)=-a(n)*m(0);for(i=1;i<n;i++)n(i)=(b(i)*m(0))+n(i+1); /*求下一次的m(n)*/n(n)=-b(n)*m(0);第三节直接算法3—单中间变量向后递推设离散传递函数为，由得中间变量，直接写出中间变量的差分方程。由得输出变量,直接写出输出变量的差分方程。画出直接算法3的信流图：for(i=1;i<n+1;i++)m(i)=0; /*初始化*/for(i=1,i<n+1,i++)m(0)=m(0)+(-a(i)*m(i)); /*求本次m(0)*/m(0)=e+m(0);for(i=1;i<n+1;i++)u(0)=u(0)+b(i)*m(i) /*求本次u(k)*/u=b(0)*m(0)+u(0);for(i=n;i>0;i--)m(i)=m(i-1); /*中间变量m(i)延时*/第四节直接算法4—单中间变量向前递推(简约快速算法)设离散传递函数为，则以和积递推的形式展开和号∑：。令，，...，，则，反变换得。M1(Z)到Mn(Z)的反变换为：，，...，，。用u(k)和mi(k)的表达式，画出直接算法4的信流图：和直接算法2一样，每一个mi(k)都是在被更新后的下一个时刻被引用，因此，信流图中的节点mi(k-1)可以直接由mi(k)兼任，并用z-1因子表明信息的延迟。这样，直接算法4的信流图简化为下图：特点：n个中间变量，在求出输出信号u(k)并从输出口输出后，还需要为下一次输出进行2n+2次乘法和2n+1次加法，显然，这一点并不影响输出的实时性。for(i=1;i<n+1;i++)m(i)=0; /*初始化*/u=b(0)*e+m(1); /*求本次u*/for(i=1;i<n;i++)m(i)=b(i)*e+(-a(i))*u+m(i+1); /*对于i<n，求下一次的m(i)*/m(n)=b(n)*e+(-a(n))*u; /*对于i=n，求下一次的m(i)*/第五节串联算法设离散传递函数为，将传递函数改写成零极点模型，按照系数相近的原则对分子和分母配对，再使分子分母中的常数项为1，得传递函数的串联算法表达式如下：在串联算法表达式中，所有的因式只有一阶分式和二阶分式两种形式。求一阶分式的时域表达式：因为，所以。从画出串联算法一阶分式的信流图,图中：m(1)=0; /*初始化*/u(k)=e(k)+m(1); /*求本次的u(k)*/m(1)=b(1)*e(k)-a(1)*u(k); /*求下一次的m(k-1)*/求二阶分式的时域表达式：令则故且从画出串联算法二阶分式的信流图：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)=b(2)*e(k)-a(2)*u(k); /*求下一次的m(2)*/第六节并联算法设离散传递函数为，通过部分分式的分解，将传递函数改写若干个部分分式的和，即，在并联算法表达式中，所有的项只有一阶分式和二阶分式两种形式。求并联算法一阶分式的时域表达式：画出并联算法一阶分式的信流图：m(1)=0; /*初始化*/u(k)=b(0)*e(k)+m(1); /*求本次的u(k)*/m(1)=-a(1)*u(k); /*求下一次的m(1)*/求并联算法二阶分式的时域表达式：，令，，则，故，由M1(Z),M2(Z)反变换得，。画出并联算法二阶分式的信流图：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

第四章数字PID控制算法PID控制算法是模拟手动调节过程而构造出来的，它有3个特点：①控制算式中比例项，积分项，微分项物理意义明确。②控制算式中比例带，积分时间，微分时间等3个参数互相独立，因而可以在闭环中用凑试法确定这3个参数，并且不必预先知道控制对象的特性。③将一阶对象或二阶对象对象与控制器形成闭环系统，假定闭环系统的传递函数为，则控制器的传递函数总是PID传递函数。由于有这3个特点，PID控制算法具有很强的普适性，因而被广泛地应用在“任意”场合。不过，若受控对象的纯滞后比较大，PID控制算法可能会出现积分饱和等问题，因此，PID控制算法出现了一些变形。在工业生产过程中，例如对温度，压力，流量等参数进行控制，使用PID控制算法有比较成熟的经验。第一节微分方程和差分方程1微分方程：。2积分分离微分方程：。3位置式差分方程：。4增量形式差分方程：令得。输出变量： 纯量，对于实际的物理量，设阀门全开时为1。输入变量： 纯量，为误差值，对于实际的物理量，设的最大值为1。参变量： 比例系数，纯量，忽略积分项和微分项，它使输出为1时输入量的变化，它是比例带的倒数。 积分常数，时间量，在忽略比例项和微分项且时，它是使输出为1时需要的积分时间。 微分常数，时间量，为微分输入的作用时间。第二节不完全微分由微分方程，得s传递函数，将微分通道串联一个一阶环节，变为，其中同，即为的若干分之一。所以不完全微分的传递函数为。经过离散化得。以下分别计算各个通道：比例通道：令，则，积分通道：令，则，，微分通道：令，则，变形后得化为时域表达上式中令和，则以上比例，积分，微分，各个通道都是位置式的，同样可以令变换为增量式。第三节参数选择采样时间受控物理量采样时间(s)备注温度15--20取纯滞后时间常数压力3--10优先选用6--8s流量1--5优先选用1--2s液位6--8优先选用7s成份15--20优先选用18s当系统中仅惯性时间常数起作用时，。当系统中纯滞后时间占有一定位置时，应选择。当系统中纯滞后时间占主导地位时，应选择。控制参数受控物理量特点(%)TI(分钟)TD(分钟)温度对象为多容，滞后较大，须用微分20-603.0-100.5-3压力对象为多容，滞后一般较大，不用微分30-700.4-3--流量对象时间常数小并有噪声，应较大，TI应较短，不用微分40-1000.1-1--液位在允许有余差时，不必用积分，不用微分20-80--成份第四节c51框架这仅仅是一个PID控制程序的基本架构，在实际应用时，还有很多细节问题需要处理。#include<string.h>#include<stdio.h>/*============================================================================================*/typedefstructPID{doubleSetPoint; //设定目标DesiredValuedoubleProportion; //比例常数ProportionalConstdoubleIntegral; //积分常数IntegralConstdoubleDerivative; //微分常数DerivativeConstdoubleLastError; //Error[-1]doublePrevError; //Error[-2]doubleSumError; //SumsofErrors}PID;/*============================================================================================*/doublePIDCalc(PID*pp,doubleNextPoint) //PID计算{doubledError,Error;Error=pp->SetPoint-NextPoint; //偏差pp->SumError+=Error; //积分dError=pp->LastError-pp->PrevError; //当前微分pp->PrevError=pp->LastError;pp->LastError=Error;return(pp->Proportion*Error //比例项+pp->Integral*pp->SumError //积分项+pp->Derivative*dError //微分项);}/*===========================================================================================*/voidPIDInit(PID*pp) //初始化PID数据结构{memset(pp,0,sizeof(PID));}/*============================================================================================*/doublesensor(void) //读AD转换器，这里给出的是一个空函数{return100.0; //这个返回值是虚设的}/*============================================================================================*/vo