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

实验二利用MATLAB进行时域分析本实验内容包含以下三个部分：基于MATLAB的线性系统稳定性分析、基于MATLAB的线性系统动态性能分析、和MATALB进行控制系统时域分析的一些其它实例。一、基于MATLAB的线性系统稳定性分析线性系统稳定的充要条件是系统的特征根均位于S平面的左半部分。系统的零极点模型可以直接被用来判断系统的稳定性。另外，MATLAB语言中提供了有关多项式的操作函数，也可以用于系统的分析和计算。（1）直接求特征多项式的根设p为特征多项式的系数向量，则MATLAB函数roots()可以直接求出方程p=0在复数范围内的解v,该函数的调用格式为：v=roots(p)例3.1已知系统的特征多项式为：123235xxxx特征方程的解可由下面的MATLAB命令得出。p=[1,0,3,2,1,1];v=roots(p)结果显示：v=0.3202+1.7042i0.3202-1.7042i-0.72090.0402+0.6780i0.0402-0.6780i利用多项式求根函数roots(),可以很方便的求出系统的零点和极点，然后根据零极点分析系统稳定性和其它性能。（2）由根创建多项式如果已知多项式的因式分解式或特征根，可由MATLAB函数poly()直接得出特征多项式系数向量，其调用格式为：p=poly(v)如上例中：v=[0.3202+1.7042i;0.3202-1.7042i;-0.7209;0.0402+0.6780i;0.0402-0.6780i];p=poly(v)结果显示p=1.00000.00013.00002.00010.99980.9999由此可见，函数roots()与函数poly()是互为逆运算的。（3）多项式求值在MATLAB中通过函数polyval()可以求得多项式在给定点的值，该函数的调用格式为：polyval(p,v)对于上例中的p值，求取多项式在x点的值，可输入如下命令：p=[1,0,3,2,1,1];x=1polyval(p,x)结果显示x=1ans=8(4)部分分式展开考虑下列传递函数：nnnnnnasasabsbsbdennumsNsM110110)()(式中0a0，但是ia和jb中某些量可能为零。MATLAB函数可将)()(sNsM展开成部分分式，直接求出展开式中的留数、极点和余项。该函数的调用格式为：),(,,dennumresiduekpr则)()(sNsM的部分分式展开由下式给出：)()()()2()2()1()1()()(sknpsnrpsrpsrsNsM式中1)1(pp,2)2(pp,…,npnp)(,为极点，1)1(rr,2)2(rr,…,nrnr)(为各极点的留数，)(sk为余项。例3.2设传递函数为：61166352)(2323sssssssG该传递函数的部分分式展开由以下命令获得：num=[2,5,3,6];den=[1,6,11,6];[r,p,k]=residue(num,den)命令窗口中显示如下结果r=-6.0000-4.00003.0000p=-3.0000-2.0000-1.0000k=2中留数为列向量r，极点为列向量p，余项为行向量k。由此可得出部分分式展开式：2132436)(ssssG该函数也可以逆向调用，把部分分式展开转变回多项式)()(sNsM之比的形式，命令格式为：[num,den]=residue(r,p,k)对上例有：[num,den]=residue(r,p,k)结果显示num=2.00005.00003.00006.0000den=1.00006.000011.00006.0000应当指出，如果p(j)=p(j+1)=…=p(j+m-1),则极点p(j)是一个m重极点。在这种情况下，部分分式展开式将包括下列诸项：mjpsmjrjpsjrjpsjr)()1()()1()()(2例3.3设传递函数为：13332)1(32)(23232sssssssssG则部分分式展开由以下命令获得：v=[-1,-1,-1]num=[0,1,2,3];den=poly(v);[r,p,k]=residue(num,den)结果显示v=-1-1-1r=1.00000.00002.0000p=-1.0000-1.0000-1.0000k=[]其中由poly()命令将分母化为标准降幂排列多项式系数向量den,k=[]为空矩阵。由上可得展开式为：0)1(2)1(011)(32ssssG（5）由传递函数求零点和极点。在MATLAB控制系统工具箱中，给出了由传递函数对象G求出系统零点和极点的函数，其调用格式分别为：Z=tzero(G)P=G.P{1}注意：式19中要求的G必须是零极点模型对象，且出现了矩阵的点运算“.”和大括号｛｝表示的矩阵元素，详细内容参阅后面章节。例3.4已知传递函数为：sssssssG5.19225.72.952.618.6)(2342输入如下命令：num=[6.8,61.2,95.2];den=[1,7.5,22,19.5,0];G=tf(num,den);G1=zpk(G);Z=tzero(G)P=G1.P{1}结果显示Z=-7-2P=0-3.0000+2.0000i-3.0000-2.0000i-1.5000（6）零极点分布图在MATLAB中，可利用pzmap()函数绘制连续系统的零、极点图，从而分析系统的稳定性，该函数调用格式为：pzmap(num,den)例3.5给定传递函数:2724364523)(2345234ssssssssssG利用下列命令可自动打开一个图形窗口，显示该系统的零、极点分布图，如图3-1所示。num=[3,2,5,4,6];den=[1,3,4,2,7,2];pzmap(num,den)title(¹Pole-ZeroMap¹)%图形标题。-2-1.5-1-0.500.5-1.5-1-0.500.511.5Pole-ZeroMapRealAxisImaginaryAxis二.系统动态特性分析。控制系统软件包提供了控制系统工程需要的基本的时域与频域分析工具函数。连续时间系统分析函数Impulse脉冲响应Step阶跃响应Lsim任意输入的仿真Bode波特图Nyquist奈奎斯特图Lyap李雅普诺夫方程Gram可控性与可观性（1）时域响应解析算法――部分分式展开法用拉氏变换法求系统的单位阶跃响应，可直接得出输出c(t)随时间t变化的规律，对于高阶系统，输出的拉氏变换象函数为：sdennumssGsC11)()(对函数c(s)进行部分分式展开，我们可以用num,[den,0]来表示c(s)的分子和分母。例3.6给定系统的传递函数：2450351024247)(23423ssssssssG用以下命令对ssG)(进行部分分式展开。num,[den,0]num=[1,7,24,24]den=[1,10,35,50,24][r,p,k]=residue(num,[den,0])输出结果为num=172424den=110355024r=-1.00002.0000-1.0000-1.00001.0000p=-4.0000-3.0000-2.0000-1.00000k=[]输出函数c(s)为：C（s）=num/（den*s）num=172424ans=1103550240拉氏变换得：L=laplace((s^3+7*s+24*s+24)/(s^4+10*s^3+35*s^2+50*s+24),s,t)（2）单位阶跃响应的求法：控制系统工具箱中给出了一个函数step()来直接求取线性系统的阶跃响应，如果已知传递函数为：dennumsG)(则该函数可有以下几种调用格式：step(num,den)(a)step(num,den,t)(b)或step(G)(c)step(G,t)(d)该函数将绘制出系统在单位阶跃输入条件下的动态响应图，同时给出稳态值。对于式(b)和(d)，t为图像显示的时间长度，是用户指定的时间向量。式(a)和(c)的显示时间由系统根据输出曲线的形状自行设定。如果需要将输出结果返回到MATLAB工作空间中，则采用以下调用格式：c=step(G)此时，屏上不会显示响应曲线，必须利用plot()命令去查看响应曲线。plot可以根据两个或多个给定的向量绘制二维图形，详细介绍可以查阅后面的章节。例3.7已知传递函数为：25425)(2sssG利用以下MATLAB命令可得阶跃响应曲线如图3-2所示。num=[0,0,25];den=[1,4,25];step(num,den)grid%绘制网格线title(¹Unit-StepResponseofG(s)=25/(s^2+4s+25)¹)%图像标题我们还可以用下面的语句来得出阶跃响应曲线G=tf([0,0,25],[1,4,25]);t=0:0.1:5;%从0到5每隔0.1取一个值。c=step(G,t);%动态响应的幅值赋给变量cplot(t,c)%绘二维图形，横坐标取t,纵坐标取c。Css=dcgain(G)%求取稳态值。系统显示的图形类似于上一个例子，在命令窗口中显示了如下结果Css=1（3）求阶跃响应的性能指标MATLAB提供了强大的绘图计算功能，可以用多种方法求取系统的动态响应指标。我们图3-2MATLAB绘制的响应曲线首先介绍一种最简单的方法――游动鼠标法。对于例16，在程序运行完毕后，用鼠标左键点击时域响应图线任意一点，系统会自动跳出一个小方框，小方框显示了这一点的横坐标（时间）和纵坐标（幅值）。按住鼠标左键在曲线上移动，可以找到曲线幅值最大的一点――即曲线最大峰值，此时小方框中显示的时间就是此二阶系统的峰值时间，根据观察到的稳态值和峰值可以计算出系统的超调量。系统的上升时间和稳态响应时间可以依此类推。这种方法简单易用，但同时应注意它不适用于用plot()命令画出的图形。另一种比较常用的方法就是用编程方式求取时域响应的各项性能指标。与上一段介绍的游动鼠标法相比，编程方法稍微复杂，但通过下面的学习，读者可以掌握一定的编程技巧，能够将控制原理知识和编程方法相结合，自己编写一些程序，获取一些较为复杂的性能指标。通过前面的学习，我们已经可以用阶跃响应函数step()获得系统输出量，若将输出量返回到变量y中，可以调用如下格式[y,t]=step(G)该函数还同时返回了自动生成的时间变量t，对返回的这一对变量y和t的值进行计算，可以得到时域性能指标。①峰值时间(timetopeak)可由以下命令获得：[Y,k]=max(y);timetopeak=t(k)应用取最大值函数max()求出y的峰值及相应的时间，并存于变量Y和k中。然后在变量t中取出峰值时间，并将它赋给变量timetopeak。②最大(百分比)超调量(percentovershoot)可由以下命令得到：C=dcgain(G);[Y,k]=max(y);percentovershoot=100*(Y-C)/Cdcgain()函数用于求取系统的终值，将终值赋给变量C，然后依据超调量的定义，由Y和C计算出百分比超调量。③上升时间(risetime)可利用MATLAB中控制语句编制M文件来获得。首先简单介绍一下循环语句while的使用。while循环语句的一般格式为：while循环判断语句循环体end其中,循环判断语句为某种形式的逻辑判断表达式。当表达式的逻辑值为真时，就执行循环体内的语句；当表达式的逻辑值为假时，就退出当前的循环体。如果循环判断语句为矩阵时，当且仅当所有的矩阵元素非零时，逻辑表达式的值为真。为避免循环语句陷入死循环，在语句内必须有可以自动修改循环控制变量的命令。要求出上升时间，可以用while语句编写