实验二利用MATLAB进行时域分析

1、实验二 利用 MATLAB 进行时域分析本实验内容包含以下三个部分：基于MATLAB的线性系统稳定性分析、基于MATLAB的线性系统动态性能分析、和MATALB进行控制系统时域分析的一些其它实例。一、 基于MATLAB的线性系统稳定性分析线性系统稳定的充要条件是系统的特征根均位于S平面的左半部分。系统的零极点模型可以直接被用来判断系统的稳定性。另外，MATLAB语言中提供了有关多项式的操作函数，也可以用于系统的分析和计算。(1) 直接求特征多项式的根设p为特征多项式的系数向量，则MATLAB函数roots()可以直接求出方程 p=0在复数范围内的解v,该函数的调用格式为：v=roots(p)例

2、 已知系统的特征多项式为：5 3 2x 3x 2x x 1特征方程的解可由下面的MATLAB命令得出。>> p=1,0,3,2,1,1;v=roots(p)结果显示：v =+利用多项式求根函数roots(), 可以很方便的求出系统的零点和极点，然后根据零极点分析系统稳定性和其它性能。(2) 由根创建多项式如果已知多项式的因式分解式或特征根，可由MATLAB函数poly()直接得出特征多项式系数向量，其调用格式为：p=poly(v)如上例中：v=+; ;+;>> p=poly(v)结果显示p =由此可见，函数roots()与函数poly()是互为逆运算的。(3) 多项式求

3、值在MATLAB中通过函数polyval()可以求得多项式在给定点的值，该函数的调用格式为:polyval(p,v)对于上例中的p值，求取多项式在x点的值，可输入如下命令：>> p=1,0,3,2,1,1;x=1polyval(p,x)结果显示x =1ans =(4) 部分分式展开考虑下列传递函数：M (s)num b°sn biSn 1bnN(s)den a0sn a1 sn 1an式中a。0，但是ai和bj中某些量可能为零。MATLAB函数可将器展开成部分分式，直接求出展开式中的留数'极点和余项。该函数的调用格式为:r,p,k residue( nu m,de

4、 n)则总的部分分式展开由下式给出:N(s)k(s) s P(n)M (s)r(1)r(2)応 s p(1) s p(2)式中 p(1)P1,p(2)P2,，P5)Pn，为极点,r(1)r1 ,r(2) Q,r(n)rn为各极点的留数，k(s)为余项。例设传递函数为:G(s)2s3 5s2 3s 6s3 6s211s 6该传递函数的部分分式展开由以下命令获得:>> nu m=2,5,3,6;den=1,6,11,6;r,p,k=residue( nu m,de n)命令窗口中显示如下结果中留数为列向量r，极点为列向量p,余项为行向量 k。由此可得出部分分式展开式:64G(s)32该

5、函数也可以逆向调用,把部分分式展开转变回多项式器之比的形式，命令格式为:nu m,de n=residue(r,p,k)对上例有：>> nu m,de n=residue(r,p,k)结果显示num =den =应当指出，如果p(j)=p(j+1)= -=p(j+m-1),则极点p(j)是一个m重极点。在这种情况下，部分分式展开式将包括下列诸项:r(j)s p( j)r(j 1)s p(j) 2r(j m 1)p(j)例设传递函数为:G(s)s2 2s 3(s 1)3s2 2s 3s3 3s2 3s 1则部分分式展开由以下命令获得:>> v=-1,-1,-1num=0,

6、1,2,3;den=poly(v);r,p,k=residue( nu m,de n)结果显示v =-1-1 -1其中由poly()命令将分母化为标准降幕排列多项式系数向量由上可得展开式为：den, k=为空矩阵。G(s)1 02s 1 (s 1)2(s 1)3(5) 由传递函数求零点和极点。在MATLAB控制系统工具箱中，给出了由传递函数对象G求出系统零点和极点的函数，其调用格式分别为：Z=tzero(G)P=1注意：式19中要求的G必须是零极点模型对象，且出现了矩阵的点运算“”和大括号訂表示的矩阵元素，详细内容参阅后面章节。例已知传递函数为：2、6.8s61.2s95.2G(s) 3s47

7、.5s3 22s2 19.5s输入如下命令：num=,;den=1,22,0;G=tf( nu m,de n);G1=zpk(G);Z=tzero(G)P=1结果显示Z =-7-2(6) 零极点分布图在MATLAB中，可利用pzmap()函数绘制连续系统的零、极点图，从而分析系统的稳定性，该函数调用格式为：pzmap( nu m,de n)例给定传递函数：3s4 2s3 5s2 4s 6s5 3s4 4s3 2s2 7s 2利用下列命令可自动打开一个图形窗口，显示该系统的零、极点分布图，如图3-1所示。>> nu m=3,2,5,4,6;den=1,3,4,2,7,2;pzmap

8、(nu m,de n)title( 1 Pole-Zero Map1) % 图形标题。1.5cfxa wanraay-1.5-2-ILX-一-rr-1.5-1-0.5Real Axis0Pole-Zero Map5 oo5 o-0.5系统动态特性分析。控制系统软件包提供了控制系统工程需要的基本的时域与频域分析工具函数。连续时间系统分析函数Impulse脉冲响应Step阶跃响应Lsim任意输入的仿真Bode波特图Nyquist奈奎斯特图Lyap李雅普诺夫方程Gram可控性与可观性(1)时域响应解析算法部分分式展开法用拉氏变换法求系统的单位阶跃响应，可直接得出输出c(t)随时间t变化的规律，对于高

9、阶系统，输出的拉氏变换象函数为：1num 1C(s) G(s)sden s对函数c(s)进行部分分式展开，我们可以用num,den,0来表示c(s)的分子和分母。例给定系统的传递函数：、s3 7s2 24s 24G(s) 3s410s335s250s 24用以下命令对9®进行部分分式展开。s>> nu m,de n,0num=1,7,24,24den=1,10,35,50,24r,p,k=residue( nu m,de n, 0)输出结果为num =24241103550r =den =240k =输出函数c(s)为：C( s) =num/ ( den*s )num =

10、1 72424ans =24 01103550拉氏变换得：L=laplace(sA3+7*s+24*s+24)/(sA4+10*sA3+35*sA2+50*s+24),s,t)(2 )单位阶跃响应的求法:控制系统工具箱中给出了一个函数step()来直接求取线性系统的阶跃响应，如果已知传递函数为:G(s)numden则该函数可有以下几种调用格式:step (nu m,de n)step (nu m,de n,t)或(b)step(G)(c)step(G,t)(d)该函数将绘制出系统在单位阶跃输入条件下的动态响应图，同时给出稳态值。对于式(b)和(d)，t为图像显示的时间长度，是用户指定的时间向量

11、。式(a)和(c)的显示时间由系统根据输出曲线的形状自行设定。如果需要将输出结果返回到MATLAB工作空间中，则采用以下调用格式:c=step(G)此时，屏上不会显示响应曲线，必须利用plot()命令去查看响应曲线。plot可以根据两个或多个给定的向量绘制二维图形，详细介绍可以查阅后面的章节。例已知传递函数为:G(s)25s2 4s 251 .40.3050.4UllbSt即住讯皿梵 GT G(S>23CSss-25)UU.51.522.5I me sec0图3-2 MATLAB绘制的响应曲线>>num=0,0,25;den=1,4,25;step( nu m,de n)gr

12、id %绘制网格线title( 1Unit-Step Response of G(s)=25(s2+4s+25) 1) % 图像标题我们还可以用下面的语句来得出阶跃响应曲线>> G=tf(0,0,25,1,4,25);t=0:5;%从0到5每隔取一个值。c=step(G,t); %动态响应的幅值赋给变量cplot(t,c) %绘二维图形，横坐标取t,纵坐标取C。Css=dcgai n(G) % 求取稳态值。系统显示的图形类似于上一个例子，在命令窗口中显示了如下结果Css=1(3 )求阶跃响应的性能指标MATLAB提供了强大的绘图计算功能，可以用多种方法求取系统的动态响应指标。我们首

13、先介绍一种最简单的方法一一游动鼠标法。对于例16,在程序运行完毕后，用鼠标左键 点击时域响应图线任意一点， 系统会自动跳出一个小方框， 小方框显示了这一点的横坐标 (时 间)和纵坐标(幅值)。按住鼠标左键在曲线上移动，可以找到曲线幅值最大的一点一一即曲线最大峰值， 此时小方框中显示的时间就是此二阶系统的峰值时间， 根据观察到的稳态值 和峰值可以计算出系统的超调量。 系统的上升时间和稳态响应时间可以依此类推。 这种方法 简单易用，但同时应注意它不适用于用plot()命令画出的图形。另一种比较常用的方法就是用编程方式求取时域响应的各项性能指标。与上一段介绍的游动鼠标法相比，编程方法稍微复杂，但通过

14、下面的学习，读者可以掌握一定的编程技巧， 能够将控制原理知识和编程方法相结合， 自己编写一些程序， 获取一些较为复杂的性能指标。通过前面的学习， 我们已经可以用阶跃响应函数step( )获得系统输出量，若将输出量返回到变量 y 中，可以调用如下格式y,t=step(G)该函数还同时返回了自动生成的时间变量 t ，对返回的这一对变量 y 和 t 的值进行计算， 可以得到时域性能指标。 峰值时间 (timetopeak) 可由以下命令获得：Y,k=max(y);timetopeak=t(k)应用取最大值函数 max()求出y的峰值及相应的时间，并存于变量Y和k中。然后在变量 t 中取出峰值时间，并

15、将它赋给变量 timetopeak 。 最大 (百分比 )超调量 (percentovershoot) 可由以下命令得到：C=dcgain(G);Y,k=max(y);percentovershoot=100*(Y-C)/Cdcgain()函数用于求取系统的终值，将终值赋给变量C,然后依据超调量的定义，由Y 和C计算出百分比超调量。 上升时间(risetime)可利用MATLAB中控制语句编制 M文件来获得。首先简单介绍 一下循环语句 while 的使用。while 循环语句的一般格式为：while< 循环判断语句 >循环体end其中 ,循环判断语句为某种形式的逻辑判断表达式。当表

16、达式的逻辑值为真时， 就执行循环体内的语句； 当表达式的逻辑值为假时， 就退出 当前的循环体。 如果循环判断语句为矩阵时， 当且仅当所有的矩阵元素非零时， 逻辑表达式 的值为真。为避免循环语句陷入死循环， 在语句内必须有可以自动修改循环控制变量的命令。要求出上升时间，可以用 while 语句编写以下程序得到：C=dcgain(G);n=1;while y(n)<Cn=n+1;endrisetime=t(n)在阶跃输入条件下， y 的值由零逐渐增大，当以上循环满足 y=C 时，退出循环，此时对 应的时刻，即为上升时间。对于输出无超调的系统响应，上升时间定义为输出从稳态值的10%上升到 90

17、%所需时间，则计算程序如下：C=dcgain(G);n=1;while y(n)<*Cn=n+1;endm=1;while y(n)<*Cm=m+1;endrisetime=t(m)-t (n) 调节时间(setllingtime)可由while语句编程得到：C=dcgai n( G);i=le ngth(t);while(y(i)>*C)&(y(i)<*C)i=i-1;endsetlli ngtime=t(i)i的上限值。用向量长度函数length()可求得t序列的长度，将其设定为变量例已知二阶系统传递函数为：G(s)3(s 13i)(s 1 3i)利用下面的

18、程序可得到阶跃响应如图3-3及性能指标数据。G(S)=3SA2+2S+10)>> nu m=0,0,3;den=1,2,10;step( nu m,de n)gridtitle( 1Un it-Step Respo nse of G(s)=3/(sA2+2s+10) 1)或者G=tf(0,0,3,1,2,10);t=0:5;c=step(G,t);plot(t,c)>> Css=dcga in(G)Css =G=tf(0,0,3,1,2,10);>> y,t=step(G)y =0>> Y,k=max(y);>> timetopeak=t(k)timetopeak =>> C=dcgain(G);>> Y,k=max(