利用双线性变换求其离散传递函数.doc

1 设计背景11.1 设计目的11.2 设计内容和要求11.3 设计工作任务及工作量的要求12 双线性变换及其原理22.1 双线性变换的定义22.2 双线性变换的原理22.2.1 公式的推导22.2.2 公式的验证22.2.2 设计步骤42.3 双线性变换的主要特性63 计算机实现程序框图74 理论计算85 程序验证106 结果分析11参考文献13附表 程序清单141 设计背景1.1 设计目的本课程设计以自动控制理论、现代控制理论、MATLAB及应用等知识为基础，利用双线性变换求连续系统对应的离散化的系统，目的是使学生在现有的控制理论的基础上，学会用MATLAB语言编写控制系统的离散化的程序，通过上机实习加深对课堂所学知识的理解，掌握一种能方便地对系统进行离散化的设计工具。1.2 设计内容和要求1 在理论上对连续系统采用双线性变换求离散化推导出算法和计算公式。2 画出计算机实现算法的框图。3 编写程序并调试和运行。4 以下面的系统为例，进行计算。已知系统闭环传递函数，利用双线性变换求其离散传递函数。5 分析运算结果（离散化步长对系统性能的影响）。6 程序应具有一定的通用性，对不同参数能有兼容性。1.3 设计工作任务及工作量的要求1 本次课程设计要求每周学生至少见指导教师4次，其中集中辅导答疑部不于3次。2 设计说明书的格式按设计说明书格式要求，采用word软件排版，计算机打印。（具体包括：封皮、目录、正文、参考文献等）3 程序清单用A4纸打印后，作为附录订装在说明书后面。4 框图和其他图表放在正文中。2 双线性变换及其原理2.1 双线性变换的定义双线性变换法又称突斯汀（Tustin）法，是一种基于梯形积分规则的数字积分变换方法。2.2 双线性变换的原理2.2.1 公式的推导双线性替换公式（或称“突斯汀（Tustin）”公式），它可以从梯形积分公式中直接推导出来。按这种替换公式进行替换，可以保证G(z)的稳定性，而且，具有一定的仿真精度。已知梯形积分公式为： 即： 则有： 即： （2.1）（2.1）式称为双线性替换公式，也可写成为 （2.2）2.2.2 公式的验证由（2.2）式可得, （2.3）由（2.3）式可知，若0，则0，则1。这就是说，若采用（2.1）式的替换公式，Z平面上的单位圆，映射到S平面上将是整个左半平面，其逆也真，如图2.1所示。因此，利用（2.1）式的替换公式，如果原来G(s)稳定，那么G(z)也是稳定的。例如，假定传递函数图2.1 双线性变换的映射关系 （2.4）则根据（2.1）式进行替换，得： （2.5）将（2.5）式写成差分方程式，得 （2.6）为了说明利用双线性替换公式（2.1）式所得的仿真模型具有一定的精度，我们用频率分析法来分析一下G(s)与G(z)两者之间的误差。根据（2.4）式的G(s)，用代入，可得： （2.7）它的幅度与相频的具体数据如表2.1所示。对于（2.5）式所示的G(z)，同样可以用代入，可得： （2.8）用代入，同样可以求出它的幅频特性与相频特性。当T1s时，它们的具体数值也列于表2.1中。比较表2.1中的两组数据可知：当T1s时，由G(s)与G(z)所获得的频率特性在 rad/s之内二者是十分接近的 (幅频的误差仅1.5%，相频的误差在5之内)。这说明，用双线性替换公式所获得的仿真模型不仅稳定，而且有一定的仿真精度。表2.1 用双线性替换法所获得的仿真模型与实际连续系统的比较(rad/s)0.10.30.60.81.01.21.5幅频特性连续系统0.099010.27520.44120.48780.500.49180.4615线性化替换法0.09910.2770.4470.4930.4980.4760.417相频特性连续系统748.5856.6028.0712.680-10.39-22.62线性化替换法78.5756.3626.519.57-5.06-17.67-33.562.2.2 设计步骤双线性替换不仅可以获得精度极爱搞的仿真模型，而且能利用计算机程序来实现这种替换。下面来介绍一种程序替换法。设线性系统的传递函数为 （2.9）在双线性替换下得到的Z传递函数为： （2.10）现需要由ai，bi(i=0，1，n)确定di，ei(i=0，1，n)，若直接将双线性替换公式代入，可得： 将其分子、分母同乘以，可得： =将，写成向量形式 由于，均为n阶多项式，可得到： 矩阵为（n+1）(n+1)阶；其中第一行诸元素为的展开式的各系数，第一列诸元素为1，阶次n确定后，这些元素均为已知，并可以证明，其余nn个元素可由下式求得： i，j=1，2,，n （2.11）从而可得： diag diag这样就得到了分子分母的各系数的表达式： （2.12） （2.13）其中 ， （2.14）2.3 双线性变换的主要特性 若D(s)稳定，则D(z)一定稳定变换前后，稳态增益不变。双线性变换后D(z)的阶次不变，且分子、分母具有相同的阶次。3 计算机实现程序框图开始输入分子分母多项式num，den和步长T判断该系统是几阶系统显示S传递函数将原系统的num,den的向量顺序颠倒，得到新的num,den生成以（1,2/T,(2/T)n）为对角线的矩阵H生成Xi,j矩阵A离散系统分子多项式dnum=num*H*A;离散系统分母多项式dden=den*H*A;显示Z传递函数结束图3.1 仿真程序流程图4 理论计算 (4.1)将代入（4.1）式得上下同时乘以=将A(z),B(z)写成向量形式由于，,均为3阶多项式，可得到： 这样就得到了分子分母的各系数的表达式： 取步长T=1时 5 程序验证运行程序后，显示如下：请输入分子多项式num=0,0,0,4请输入分母多项式den=1,3,2,0请输入步长T=1 Transfer function: 4-s3 + 3 s2 + 2 s n = 3 4 0 0 0 0 2 3 1H = 1 0 0 0 0 2 0 0 0 0 4 0 0 0 0 8 1 3 3 1 1 1 -1 -1 1 -1 -1 1 1 -3 3 -1dnum = 4 12 12 4dden = 24 -32 8 0 Transfer function:4 z3 + 12 z2 + 12 z + 4- 24 z3 - 32 z2 + 8 z Sampling time: 1程序出来的结果与理论计算的结果相同，所以程序没有错误，并且能够实现双线性变换的仿真。6 结果分析分析离散化步长对系统性能的影响：由于对于连续系统离散化的问题，这里的采样时间和步长一般是一样的。所以之后在这里分析T对系统性能的影响。在这里我以上面例题为例，分别取T=0.5,1,2,3,4，通过matlab程序求其性能指标tp,tr,ts和超调量pos，程序详见附录程序清单。结果如表6.1，图像如图6.1：表6.1 不同的T对系统性能的影响 性能指标 tptrtspos连续系统3.01261.84137.489569.9872T=0.53.52.541.567.4572T=1434862.5573T=2648037.3825T=36612638.2260T=48818434.6113(a) （b） (c) (d) (e) (f)图6.1 连续系统及T取不同值时的单位阶跃响应曲线由以上结果可知，当T的取值越大，系统的性能指标tp,tr,ts都随之增大，超调量随之减小。T的取值越小，与原系统的性能指标越相近。参考文献1孙增晰.计算机辅助设计.北京：清华大学出版社，19952胡寿松.自动控制原理.北京：科学出版社，20113MATLAB及应用.吉林：东北电力大学出版社，20114王小丹.基于MATLAB的系统分析与设计.西安：西安电子科技大学出版社，2007附表 程序清单1.实现双线性变换的matlab程序如下：num=input(请输入分子多项式num=); %键盘输入分子多项式den=input(请输入分母多项式den=); %键盘输入分母多项式T=input(请输入采样周期T=); %键盘输入采样周期Tsys=tf(num,den) %构成S传递函数n=length(den)-1 %计算为n阶系统for i=1:(n+1)/2 d=num(i);num(i)=num(n-i+2);num(n-i+2)=d; e=den(i);den(i)=den(n-i+2);den(n-i+2)=e;end %颠倒原系统的分子分母多项式内向量的顺序disp(num) %显示颠倒后的numdisp(den) %显示颠倒后的dena(1)=1;b=1,1;c=1,1;for i = 2:n+1, a(i) =a(i-1)*(2/T);end %构造内部向量以1为首项2/T为倍数的等比序列aH=diag(a) %显示对角线为序列a的矩阵Hfor i=1:n-1 b=conv(b,c);end %求出（z+1）的n次方的系数序列bfor i=1:n+1, A(i,1)=1;end %矩阵A的第一列为1for j=1:n+1, A(1,j)=b(j);end %矩阵A的第一行为b的系数for i=2:n+1, for j=2:n+1, A(i,j)=A(i-1,j)-A(i,j-1)-A(i-1,j-1); endend %矩阵A的其他行，列的值disp(A) %显示矩阵Adnum=num*H*A %得到离散系统的分子多项式dden=den*H*A %得到离散系统的分母多项式sys=tf(dnum,dden,T) %构成Z传递函数2.求离散闭环系统的性能指标的matlab程序如下：num=0,0,0,4;den=1,3,2,0;G0=tf(num,den);numb,denb=cloop(num,den,-1); %计算闭环系统的分子分母多项式sys=tf(numb,denb);figure(1) %在图1上画图step(sys,120) %画时间为120秒的单位响应曲线y,t=step(sys,120); %求出120秒内系统的y值和t值mp,f=max(y); %求y的最大值及其所在位置tp=t(f) %求峰值时间tpg=length(t); %计算t的长度yss=y(g);pos=100*(mp-yss)/yss %求超调量posfor i=1:g if y(i)=1 tr=t(i) break; endend %求上升时间trfor i=g:-1:1 if y(i)=1.02*yss ts=t(i) break; endend %求调节时间tsn=length(den)-1;for i=1:(n+1)/2 d=num(i);num(i)=num(n-i+2);num(n-i+2)=d; e=den(i);den(i)=den(n-i+2);den(n-i+2)=e;endTs=0.5,1,2,3,4; %T分别取0.5,1,2,3,4for w=1:length(Ts) %T在取不同的值时分别执行下面的求值过程 T=Ts(w) a(1)=1; b=1,1;c=1,1; for i = 2:n+1, a(i) =a(i-1)*(2/T); end for i=1:n-1 b=conv(b,c); end H=diag(a); for i=1:n+1, A(i,1)=1; end for j=1:n+1, A(1,j)=b(j); end for i=2:n+1, for j=2:n+1, A(i,j)=A(i-1,j)-A(i,j-1)-A(i-1,j-1); end end dnum=num*H*A; dden=den*H*A; sys=tf(dnum,dden,T); dnumb,ddenb=cloop(dnum,dden,-1); %求其闭环系统的分子分母 figure(w+1) %在图w+1上画图 dstep(dnumb,ddenb,120) %画闭环系统的单位阶跃响应 y,t=dstep(dnum