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

课程设计报告学生姓名:学号：学院:班级:题目: 专业方向课程设计仿真组-利用双线性变换求其离散传递函数指导教师：孟杰姜文娟职称: 副教授讲师2017年 1月 13日目录1.题目背景与意义 (1)2.设计题目介绍 (1)2.1 设计内容和要求 (1)2.2设计工作任务及工作量的要求 (1)3 双线性变换法 (1)3.1 双线性变换的定义 (1)3.2 双线性变换法的优缺点 (2)3.2.1 双线性变换法的优点 (2)3.2.2 双线性变换法的缺点 (2)3.2 双线性变换的原理 (2)3.3双线性变换的主要特性 (2)4 设计步骤 (3)5 理论计算 (5)6 结果分析 (6)参考文献 (10)附录1 程序清单 (11)附录2 计算机实现程序框图 (17)1.题目背景与意义本课程设计以自动控制理论、现代控制理论、MATLAB 及应用等知识为基础，利用双线性变换求连续系统对应的离散化的系统，目的是使学生在现有的控制理论的基础上，学会用MATLAB 语言编写控制系统的离散化的程序，通过上机实习加深对课堂所学知识的理解，掌握一种能方便地对系统进行离散化的设计工具。

2.设计题目介绍2.1 设计内容和要求1 在理论上对连续系统采用双线性变换求离散化推导出算法和计算公式2 画出计算机实现算法的框图3 编写程序并调试和运行4 以下面的系统为例，进行计算 已知系统闭环传递函数)2)(1(4)(++=s s s s G ，利用双线性变换求其离散传递函数。

5 分析运算结果（离散化步长对系统性能的影响）6 程序应具有一定的通用性，对不同参数能有兼容性。

2.2设计工作任务及工作量的要求1 本次课程设计要求每周学生至少见指导教师2次，其中集中辅导答疑部不于3次。

2 设计说明书的格式按设计说明书格式要求，采用word 软件排版，计算机打印。

（具体包括：封皮、目录、正文、参考文献等）3 程序清单用A4纸打印后，作为附录订装在说明书后面。

4 框图和其他图表放在正文中。

3 双线性变换法3.1 双线性变换的定义双线性变换法又称突斯汀（Tustin ）法，是一种基于梯形积分规则的数字积分变换方法。

3.2 双线性变换法的优缺点3.2.1 双线性变换法的优点双线性变换的主要优点:靠频率的严重非线性关系得到S 平面与Z 平面的单值一一对应关系，整个j Ω轴单值对应于单位圆一周，其中ω和Ω为非线性关系。

在零频率附近，Ω~ω接近于线性关系，Ω进一步增加时，ω增长变得缓慢，(ω终止于折叠频率处)，所以双线性变换不会出现由于高频部分超过折叠频率而混淆到低频部分去的现象。

3.2.2 双线性变换法的缺点双线性变换法的缺点:Ω与ω的非线性关系，导致数字滤波器的幅频响应相对于模拟滤波器的幅频响应有畸变，(使数字滤波器与模拟滤波器在响应与频率的对应关系上发生畸变)。

例如，一个模拟微分器，它的幅度与频率是线性关系，但通过双线性变换后，就不可能得到数字微分器。

3.2 双线性变换的原理双线性变换公式可以从梯形积分公式中直接推导出来。

按这种替换公式进行替换，可以保证)(z G 的稳定性，而且，具有一定的仿真精度。

已知梯形积分公式为：)(211++++=k k k k x x Tx x 即： xz Tx z )1(2)1+=-（ 则有： 1121-+=-z z T s即： 112+-=z z T s （3.1）（3.1）式称为双线性变换公式，也可写成为：2121sTsT z -+=（3.2）3.3双线性变换的主要特性1、若)(s D 稳定，则)(z D 一定稳定；2、变换前后，稳态增益不变；3、双线性变换后的阶次不变，且分子、分母具有相同的阶次。

4 设计步骤双线性替换不仅可以精度极高的仿真模型，而且能利用计算机程序来实现这种替换。

下面来介绍一种程序替换法。

折线形系统的传递函数为nn n n nn n n b s b s b s b a s a s a s a s E s U s G ++++++++==----111011100......)()()( （4.1） 在双线性替换下得到的Z 传递函数为：；nn n n nn n n e z e z e z e d z d z d z d z E z U z G ++++++++==----111011100......)()()( （4.2） 现需要由i i b a ，（i=0，1，......，n ）确定i i e d ，（i=0，1，......，n ），若直接将双线性替换公式代入)(0s G ，可得：nn n n n n nn n n n n b z z T b z z T b z z T b a z z T a z z T a z z T a z G ++-+++-++-++-+++-++-=------)11)(2(...)11()2()11()2()11)(2(...)11()2()11()2()(11110111100 （4.3） 将其分子、分母同时乘以nz ）（1+，可得：)()()1(...)1()1(2)1(2)1(...)1()1(2)1(2)(111011100z B z A z b z z T b z T b z a z z T a z T a z G n n n n n n n n n n nn =++++-+-++++-+-=----）（）（）（）（ （4.4） 将)(z A ，)(z B 写成向量形式，()⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡--++⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=---n n nn n n nz z z z T TTT a a a a z A )1(...)1()1()1()2(0)2(...)2(0)2(...)(1110011由于nz ）（1+，，......，nz ）（1-均为n 阶多项式，可得到：⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡--+-++---1........................)1(...)1(1...)1(111101011110001001z z z x x x x x x x x x x x x z z z z z z n n nn n n in i i n n ii i n n n）（）（）（ 矩阵ij x 为(n+1)×（n+1）阶；其中第一行诸元素为nz ）（1+的展开式的各系数，第一列诸元素为1，阶次n 确定后，这些元素均为已知，并可以证明，其余n ×n个元素可由下式求得：1,11,,1,------=j i j i j i j i x x x x i,j=1，2，......，n从而可得：()[][]Tn nj i n n nz z z x T Tdiag a a a a z A 1...)2( (2)1...)(1,011--⎥⎦⎤⎢⎣⎡=这样就得到了)(0z G 分子分母的各系数的表达式：()()A H a a a a d d d n n n **=-01110......（4.5）其中：⎥⎦⎤⎢⎣⎡=n T T diag H )2( (2)1[]j i x A ,=5 理论计算)2)(1(4)(++=s s s s Gss s s G 234)(23++=（5.1）将112+-=z z T s 代入（5.1）式得上下同时乘以31）（+z将A(z)，B(z)写成向量形式：由于，，2)1)(1(-+z z ，均为3阶多项式，可得到：这样就得到了G(z)的分子、分母的各系数的表达式：()()⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡------⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎣⎡=1331111111111331)2()2(210004323210T TTd d d d⎥⎦⎤⎢⎣⎡-+-+----++=32323232812424124241248124T T T T T T T T T T T T 采样周期取T=1时，()[]0832243210-=e e e e6 结果分析分析离散化步长对系统性能的影响：由于对于连续系统离散化的问题，采样周期和步长一样。

所以在这里仅分析T 对系统性能的影响。

以为例，分别，通过MATLAB 程序求其性能指标p t ，r t ，s t 和超调量pos 。

表7.1 不同的T 对系统性能的影响性能指标 s t r t posp t 连续系统 37.5515 1.7458 70.0161 2.9975 T=0.2 32.8000 2 69.9109 3.2000 T=0.5 39 2.5000 66.2127 3.5000 T=1 48 3 62.5573 4 T=5 255 10 29.7149 10 T=10 700 20 14.7137 40 T=1002001.55572600(a) 连续系统(b) T=0.2(c) T=0.5(e) T=5(g) T=100图7.1 连续系统及T 取不同值时单位阶跃响应曲线由以上结果可知，当T 的取值越大时，系统的性能指标，s t ，r t 都随之增大，超调量随之减小，超过某一值时，系统不在稳定；T 的取值越小，与原系统的性能指标越接近。

参考文献[1] 孙增晰.计算机辅助设计.北京：清华大学出版社，1995[2] 郑大钟.线性系统理论.北京：清华大学出版社，2004[3] 姚全珠.软件技术基础.北京：电子工业出版社，2006[4] 王小丹.基于MATLAB 的系统分析与设计.西安：西安电子科技大学出版社，2007附录1 程序清单1 利用双线性变换求其离散传递函数的计算机实现num=input('请输入分子多项式num=');den=input('请输入分母多项式den=');T=input('请输入步长T=');sys=tf(num,den)n=length(den)-1 %计算传递系统阶数for i=1:(n+1)/2d=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)disp(den)a(1)=1;b=[1,1];c=[1,1];for j=2:n+1;a(j)=a(j-1)*(2/T);end%构造内部向量以1为首项2/T为倍数的等比序列H=diag(a) %显示对角矩阵%求出z+1的n次方的序列的系数矩阵for i=1:n-1;b=conv(b,c);endfor i=1:n+1;A(i,1)=1;endfor j=1:n+1;A(1,j)=b(j);endfor 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);endenddisp(A) %显示矩阵Adnum=num*H*A %得到离散后的分子系数矩阵dden=den*H*A %得到离散后的分母系数矩阵sys=tf(dnum,dden,T)运行程序，显示如下：> tustin请输入分子多项式num=[0 0 0 4]请输入分母多项式den=[1 3 2 0]请输入步长T=[1]sys =4-----------------s^3 + 3 s^2 + 2 sContinuous-time transfer function.n =34 0 0 00 2 3 1H =1 0 0 00 2 0 00 0 4 00 0 0 81 3 3 11 1 -1 -11 -1 -1 11 -3 3 -1 dnum =4 12 12 4 dden =24 -32 8 0 sys =4 z^3 + 12 z^2 + 12 z + 4-------------------------24 z^3 - 32 z^2 + 8 zSample time: 1 seconds Discrete-time transfer function.2 分析离散化步长对系统性能的影响的计算机实现num=input('请输入分子多项式num=');den=input('请输入分母多项式den=');G0=tf(num,den);[bnum,bden]=cloop(num,den,-1); %计算闭环系统的分子分母多项式sys=tf(bnum,bden);figure(1) %在图1上画图step(sys,150) %画时间为150s的单位阶跃响应曲线[y,t]=step(sys,150); %求出150s内系统的y值和t值[mp,f]=max(y); %求y的最大值及其所在位置tp=t(f) %求峰值时间l=length(t); %计算t的长度ys=y(l);pos=100*(mp-ys)/ys %求超调量posfor i=1:lif y(i)>=1tr=t(i)breakendend%求上升时间trfor j=l:-1:1if y(j)<=0.98*ys,y(j)>=1.02*ysts=t(j)breakendend%求调节时间tsn=length(den)-1;for i=1:(n+1)/2d=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.2,0.5,1,5,10,100]; %T分别取0.2，0.5，1，5，10 ，100for w=1:length(Ts) %T在取不同值是分别执行下面的求值过程T=Ts(w)a(1)=1;b=[1,1];c=[1,1];for j=2:n+1;a(j)=a(j-1)*(2/T);end%构造内部向量以1为首项2/T为倍数的等比序列H=diag(a) %显示对角矩阵for i=1:n-1;b=conv(b,c);endfor j=1:n+1;A(j,1)=1;endfor i=1:n+1;A(1,i)=b(i);endfor j=2:n+1;for i=2:n+1;A(j,i)=A(j-1,i)-A(j,i-1)-A(j-1,i-1);endenddisp(A) %显示矩阵Adnum=num*H*A %得到离散后的分子系数矩阵dden=den*H*A %得到离散后的分母系数矩阵sys=tf(dnum,dden,T)[dbnum,dbden]=cloop(dnum,dden,-1);figure(2)dstep(dbnum,dbden,150);[y,t]=dstep(dbnum,dbden);[mp,f]=max(y); %求y的最大值及其所在位置tp=f*T %求峰值时间l=length(y); %计算t的长度ys=y(l);pos=100*(mp-ys)/ys %求超调量posfor i=1:lif y(i)>=1tr=i*T %求上升时间break;endendfor i=l:-1:1if y(i)<=0.98*ys,y(i)>=1.02*ys;ts=i*T %求调节时间breakendendend附录2 计算机实现程序框图。