控制系统数字仿真与cad第7章控制系统的计算机辅助分析

1本章内容1利用MATLAB求取系统的各种特性函数；2利用MATLAB分析系统的稳定性；3利用MATLAB求取系统在典型和任意输入信号作用下的时域响应；4利用MATLAB绘制系统的根轨迹，在根轨迹上可确定任意点的根轨迹增益K值，从而得到系统稳定的根轨迹增益K值范围；5利用MATLAB绘制系统的BODE图、NICHOLS图和NYQUIST图等，并求取系统的幅值裕量和相位裕量；6利用MATLAB分析具有时间延迟系统的频率特性；7求取频率响应数据，且根据频率响应数据辩识系统的模型参数8利用MATLAB分析系统的能控性和能观测性，并能对不完全能控或不完全能观测的系统进行结构分解。第7章控制系统的计算机辅助分析2711矩阵函数1矩阵行列式矩阵AAIJ的行列式定义为ADETA1KA1K1A2K2ANKNMATLAB求矩阵行列式函数的调用格式为DETA计算矩阵的行列式有多种算法,在MATLAB中采用的方法为LU分解法。71系统的特性函数32矩阵求逆对于一个已知的NN维非奇异方阵A来说，如果有一个同样大小的C矩阵满足ACCAI式中I为单位阵，则称C矩阵为A矩阵的逆矩阵，并记作CA1。MATLAB提供了一个求取逆矩阵的函数INV，其调用格式为INVA43矩阵的迹假设一个方阵为AAIJ,I,J1,2,N则矩阵A的迹定义为亦即矩阵的迹为该矩阵对角线上各个元素之和。由代数理论可知矩阵的迹和该矩阵的特征值之和是相同的。在MATLAB中提供了求取矩阵迹的函数TRACE，其调用方法为TRACEANIIIAATR154矩阵的秩对于NM维的矩阵A，若矩阵所有的列向量中共有RC个线性无关，则称矩阵的列秩为RC，如果RCM，则称A为列满秩矩阵，相应地，若矩阵A的行向量中有RR个是线性无关的，则称矩阵A的行秩为RR,如果RRN，则称A为行满秩矩阵。MATLAB提供了一个内部函数RANK来用数值方法求取一个已知矩阵的秩，其调用格式为KRANKA65矩阵的三角分解矩阵的三角分解又称为LU分解，它的目的是将一个矩阵A分解成一个下三角矩阵L和一个上三角矩阵U的乘积，亦即可以写成ALU。在MATLAB下也给出了矩阵的LU分解函数LU，该函数的调用格式为L,ULUA76矩阵的奇异值分解CONDA7矩阵的范数NNORMA,选项选项意义无矩阵的最大奇异值，即A22与默认方式相同，亦为A21矩阵的1范数,即A1INF或INF矩阵的无穷范数,即AFRO矩阵的F范数,即AFSQRTATAIIINF只可用于向量,AMINAI数值P对向量可取任何整数,而对矩阵只可取1，2，INF或FRO表72矩阵范数函数的选项定义88矩阵的特征值与特征向量V，DEIGA其中A为要处理的矩阵，D为一个对角矩阵，其对角线上的元素为矩阵A的特征值，而每个特征值对应的V矩阵的列为该特征值的特征向量，该矩阵是一个满秩矩阵，它满足AVVD，且每个特征向量各元素的平方和（即2范数）均为1。9矩阵的特征多项式、特征方程和特征根MATLAB提供了求取矩阵特征多项式系数的函数POLY），其调用格式为PPOLYA9其中A为给定的矩阵，返回值P为一个行向量，其各个分量为矩阵A的降幂排列的特征多项式系数。即PA0A1ANMATLAB中根据矩阵特征多项式求特征根的函数为ROOTS，其调用格式为VROOTSP其中P为特征多项式的系数向量，而V为特征多项式的解，即原始矩阵的特征根。10712数值分解1数据处理如果给定一组数据XI，I1,2,N,则可利用MATLAB将这些数据用一个向量表示出来，即XX1,X2,XN利用MATLAB的函数MAX和MIN便可求出这组数据的最大和最小值，命令格式如下XM,IMAXXXM,IMINX112按实部或幅值对特征值进行排序MATLAB的函数ESORT和DSORT的调用格式为S,NDXESORTP或S,NDXDSORTP其中ESORTP对连续系统，根据实部按递减顺序对矢量P中的复特征值进行排序，NDX为索引矢量123多项式及多项式矩阵求值如果FX函数为下面的一个多项式FXA0XNA1XN1AN1XAN则可以求出该函数的导数函数为FXNA0XN1N1A1XN2AN1在MATLAB中提供了多项式求值函数POLYVAL和多项式求导的函数POLYDER它们的调用格式分别为FX0POLYVALP,X0及P1POLYDERP134数值积分对于函数FX的定积分可利用MATLABQUAD函数来求取定积分，调用格式为YQUAD函数名,A,B,TOLBAXXFYD145非线性方程求解对于非线性方程FX,T0MATLAB提供了函数FZERO和FSOLVE可以方便地求得解，它们的调用格式分别为XFZERO函数名，初值和XFSOLVE函数名，初值156函数的极值MATLAB中提供了两个基于单纯形算法求解函数最小值的函数FMIN和FMINS，它们分别对应于单变量和多变量的最优化问题的求解，其调用格式分别为XFMIN函数名，初值，选项XFMINS函数名，初值，选项167函数图形的绘制MATLAB还允许调用FPLOT函数来直接绘制出函数的图形，该函数的调用格式为FPLOT函数名,T0,TF,N其中T0,TF为绘制区间,N为点数，N默认时将取N25，函数名既可为自定义的任意M函数，也可为基本数学函数。17713求系统的阻尼系数和自然频率MATLAB函数DAMP和DDAMP的调用格式为N,DAMPAN,PDDAMPAN,PDDAMPA,TS18714求控制系统的增益MATLAB的函数DCGAIN的调用格式为KDCGAINNUM,DEN或KDCGAINA,B,C,D其中返回值K为系统增益，即放大系数。对于单变量系统，K为标量；对于多变量系统，K为向量或矩阵。193显示打印线性系统MATLAB的PRINTSYS函数可按特殊格式打印出状态空间和传递函数表示的系统。PRINTSYSNUM,DEN,S显示打印连续系统的传递函数，默认方式；PRINTSYSNUM,DEN,Z显示打印离散系统的脉冲传递函数；PRINTSYSA,B,C,D显示打印状态空间形式的系数矩阵；PRINTSYSA,B,C,D,ULABELS,YLABELS,XLABDS用ULABELS,YLABELS,XLABDS中指定的符号标记出系统矩阵A,B,C,D。2072控制系统的稳定性分析1利用极点判断系统的稳定性判断一个线性系统稳定性的一种最有效的方法是直接求出系统所有的极点，然后根据极点的分布情况来确定系统的稳定性，对于极点的求取我们在上节中已作过介绍，下面举例说明其判断方法。2利用特征值判断系统的稳定性系统的特征方程|SIA|SNA1SN1AN1SAN0的根称为系统的特征值，即系统的闭环极点。当然判断系统的稳定性同样可利用特征值来判断。213用李雅普诺夫第二法来判断系统的稳定性线性定常连续系统72在平衡状态XE0处，渐近稳定的充要条件是对任给的一个正定对称矩阵Q，存在一个正定的对称矩阵P，且满足矩阵方程ATPPAQ而标量函数VXXTPX是这个系统的一个二次型形式的李雅普诺夫函数。MATLAB提供了李雅普诺夫方程的求解函数LYAP，其调用格式为PLYAPA,QAXX2273控制系统的时域分析1任意信号函数生成任意信号函数GENSIG的调用格式为U,TGENSIGTYPE,TA或U,TGENSIGTYPE,TA,TF,T其中第一式产生一个类型为TYPE的信号序列UT，周期为TA，TYPE为以下标识字符串之一SIN正弦波；SQUARE方波；PULSE脉冲序列；第二式同时定义信号序列UT的持续时间TF和采样时间T。23例711生成一个周期为5秒，持续时间为30秒，采样时间为01秒的方波。解MATLAB窗口中执行以下命令可得图73所示结果。U,TGENSIGSQUARE,5,30,01PLOTT,UAXIS0,30051524图73252连续系统的单位阶跃响应单位阶跃响应函数STEP的调用格式为Y,X,TSTEPNUM,DEN,T或Y,X,TSTEPA,B,C,D,IU,T式中T为选定的仿真时间向量，函数返回值Y为系统在各个仿真时刻的输出所组成的矩阵；而X为自动选择的状态变量的时间响应数据。如只想绘制出系统的阶跃响应曲线，则可以由如下的格式调用此函数STEPNUM,DEN,TSTEPA,B,C,D,T26例713对于典型二阶系统试绘制出无阻尼自然振荡频率N6，阻尼比分别为02,04,10,20时系统的单位阶跃响应曲线。2222NNNSSSG27解MATLAB程序为EXAMPLE7_13MWN6ZETA020210,20FIGURE1HOLDONFORIZETANUMWN2DEN1,2IWN,WN2STEPNUM,DENENDTITLESTEPRESPONSEHOLDOFF执行后可得如图75所示的单位阶跃响应曲线。图75283离散系统的单位阶跃响应离散系统的单位阶跃响应函数DSTEP的调用格式为Y,XDSTEPNUM,DEN,N或Y,XDSTEPG,H,C,D,IU,N式中N为选定的取样点个数，当N省略时，取样点数由函数自动选取，其余参数定义同前。29例715已知二阶离散系统试求其单位阶跃响应。解MATLAB程序为EXAMPLE7_15MNUM23415DEN11608DSTEPNUM,DENTITLEPISCRETESTEPRESPONSE执行后得如图77所示的单位阶跃响应曲线。80615143222ZZZZZG图77304单位脉冲响应单位脉冲响应函数IMPULSE和DIMPULSE与单位阶跃函数STEP和DSTEP的调用格式完全一致，这里就不一一列写了。315系统的零输入响应对于连续系统由初始状态所引引起的响应，即零输入响应，可由函数INITIAL来求得，其调用格式为Y,X,TINITIALA,B,C,D,X0或Y,X,TINITIALA,B,C,D,X0,T其中X0为初始状态，其余参数定义同前。326任意输入函数的响应连续系统对任意输入函数的响应可利用MATLAB的函数LSIM求取，其调用格式为Y,XLSIMNUM,DEN,U,T或Y,XLSIMA,B,C,D,IU,U,T其中U为给定输入序列构成的矩阵，它的每列对应一个输入，每行对应一个新的时间点，其行数与时间T的长度相等。其他用法同STEP函数。3374根轨迹法所谓根轨迹是指，当开环系统的某一参数从零变化到无穷大时，闭环系统的特征方程根在S平面上所形成的轨迹。一般地，将这一参数选作开环系统的增益K，而在无零极点对消时，闭环系统特征方程的根就是闭环传递函数的极点。341绘制系统的零极点图PZMAP函数可绘制系统的零极点图，其调用格式为P,ZPZMAPA,B,C,D或P,ZPZMAPP,Z，P,ZPZMAPNUM,DEN其中列向量P为系统极点位置，列向量Z为系统的零点位置。可通过PZMAPP,Z绘制出零极点图，图中的极点用“”表示，零点用“O”表示。352绘制系统的根轨迹对于图712所示的负反馈系统，其特征方程可表示为或利用RLOCUS函数可绘制出当开环增益K由0至变化时，闭环系统的特征根在S平面变化的轨迹，该函数的调用格式为R,KRLOCUSNUM,DENR,KRLOCUSNUM,DEN,K或R,KRLOCUSA,B,C,DR,KRLOCUSA,B,C,D,K其中返回值R为系统的闭环极点，K为相应的增益。01SHSG01SDENSNUMKGSHS图71236在系统分析过程中，常常希望确定根轨迹上某一点处的增益值K，这时可利用MATLAB中的RLOCFIND函数，在使用此函数前要首先得到系统的根轨迹，然后再执行如下命令K,POLESRLOCFINDNUM,DEN或K,POLESRLOCFINDA,B,C,D执行上述命令后，将在屏幕上的图形中生成一个十字光标，使用鼠标移动它至所希望的位置，然后敲击鼠标左键即可得到该极点的位置坐标值POLES以及它所对应的增益K值。37例720已知某负反馈系统的开环传递函数为试绘制系统根轨迹，并分析系统稳定的K值范围。解MATLAB的程序为EXAMPLE7_20MNUM1DENCONV1,0,CONV1,1,1,2RLOCUSNUM,DEN,K,POLESRLOCFINDNUM,DEN21SSSKSHSG38执行以上程序，并移动鼠标到根轨迹与虚轴的交点处单击鼠标左键后可得如图713所示的根轨迹和如下结果SELECTED_POINT0000014142IK60000POLES300000000014142I0000014142I图71339由此可见根轨迹与虚轴交点处的增益K6，这说明当K6时，系统不稳定；利用RLOCFIND函数也可找出根轨迹从实轴上的分离点处的增益K038,这说明当0SUBPLOT2,1,1SEMILOGXW,20LOG10MAGSUBPLOT2,1,2SEMILOGXW,PHASE在同一个窗口上同时绘制出系统的BODE图了，其中前一条命令中对幅值向量MAG求取分贝DB值。44如果只想绘制出系统的BODE图，而对获得幅值和相位的具体数值并不感兴趣，则可以采用如下简单的调用格式BODENUM,DEN,BODEA,B,C,D,IU,）或更简单地BODENUM,DENBODEA,B,C,D,IU45例722已知二阶系统的开环传递函数为绘制出取不同值时,N在01到10之间的频率范围上的BODE图。解当N6，取02,04,06,08,10时二阶系统的BODE图可直接采用BODE函数得到。MATLAB程序为EXAMPLE7_22MWN6ZETA020210WLOGSPACE1,1FIGURE1NUMWN22NN22N2SSSG46FORKZETADEN12KWNWN2MAG,PHASE,W1BODENUM,DEN,WSUBPLOT2,1,1HOLDONSEMILOGXW1,MAGSUBPLOT2,1,2HOLDONSEMILOGXW1,PHASEENDSUBPLOT2,1,1GRIDONGRIDGRIDTITLEBODEPLOTXLABELFREQUENCYRAD/SESYLABELGAINDBSUBPLOT2,1,2GRIDONGRIDXLABELFRQUENCYRAD/SELYLABELPHASEDEGHOLDOFF执行后得如图715所示BODE图。从图中可以看出，当0时，相角趋于0，当时，相角趋于1800，当N时，相角等于900，此时的幅值也最大。图715473幅值裕量和相位裕量在判断系统稳定性时，常常需要求出系统的幅值裕量和相位裕量。利用MATLAB控制系统工具箱提供的MARGIN函数可以求出系统的幅值裕量与相位裕量，该函数的调用格式为GM,PM,WCG,WCPMARGINNUM,DEN或GM,PM,WCG,WCPMARGINA,B,C,D式中GM和PM分别为系统的幅值裕量和相位裕量，而WCG和WCP分别为幅值裕量和相位裕量处相应的频率值。48例723给定系统的开环状态空间表达式为求系统的幅值裕量和相位裕量，并画出BODE图。001875156210005422048213562100001000010TXTYTUTXTX49EXAMPLE7_23MWLOGSPACE1,1A0100001000016252138204254B0001C1562187500D0GM,PM,WCG,WCPMARGINA,B,C,DBODEA,B,C,D,1,WDISP幅值裕量,NUM2STR20LOG10GM,DB,相位裕量,NUM2STRPM,度。TITLEGAINMARGIN,NUM2STRGM,PHASEMARGINNUM2STRPM50执行后得如下数据及图716所示BODE图。图71651GM44922PM230705WCG126466WCP58275幅值裕量KG44922130492DB，相位裕量R230706度。52除了根据系统模型直接求取幅值和相位裕量之外，MATLAB的控制系统工具箱中还提供了由幅值和相位相应数据来求取裕量的方法，这时函数的调用格式为GM,PM,WCG,WCPMARGINMAG,PHASE,式中频率响应可以是由BODE函数获得的幅值和相位向量，也可以是系统的实测幅值与相位向量，为相应的频率点向量。534频率响应值MATLAB也提供了直接的求取频率响应数据的函数FREQRESP,其调用格式为FFREQRESPNUM,DEN,SQRT1或FFREQRESPA,B,C,D,IU,SQRT1式中F为频率响应，为给定的频率范围向量。545系统的奈奎斯特图（NYQUIST图）NYQUIST函数的调用格式为RE,IM,NYQUISTNUM,DENRE,IM,NYQUISTNUM,DEN,RE,IM,NYQUISTA,B,C,DRE,IM,NYQUISTA,B,C,D,IURE,IM,NYQUISTA,B,C,D,IU,其中返回值RE,IM和分别为频率特性的实部向量、虚部向量和对应的频率向量，有了这些值就可利用命令PLOTRE,IM来直接绘出系统的奈奎斯特图。55当然也可使用下面的简单命令来直接绘出系统的奈魁斯特图。NYQUISTNUM,DEN,或NYQUISTA,B,C,D更简单地NYQUISTNUM,DEN或NYQUISTA,B,C,D,IU它的使用方法基本同BODE函数的用法。56例726已知系统的开环传递函数为绘制NYQUIST图，并判断系统的稳定性。解MATLAB命令为NUM05DEN12105NYQUISTNUM,DEN执行后可得如图717所示的曲线，由于NYQUIST曲线没有包围1,J0点，且P0，所以由GSHS构成的单位负反馈闭环系统稳定。5025023SSSSHSG57图717586尼柯尔斯图NICHOLS图函数调用格式为MAG,PHASE,NICHOLSNUM,DEN,或MAG,PHASE,NICHOLSA,B,C,D,IU,可见该函数的调用格式以及返回的值与BODE函数完全一致，事实上虽然它们使用的算法不同,但这两个函数得出的结果还是基本一致的。但NICHOLS图的绘制方式和BODE图是不同的，它可由以下命令绘制PLOTPHASE,20LOG10MAG当然，NICHOLS图也可采用与BODE图类似的简单命令来直接绘制。59例729已知单位负反馈的开环传递函数为试绘制NICHOLS图。解MATLAB程序为EXAMPLE7_29MNUM1DENCONV1,0,CONV1,1,02,1WLOGSPACE1,1,400MAG,PHASENICHOLSNUM,DEN,WPLOTPHASE,20LOG10MAGNGRID执行后可得如图721所示的NICHOLS图。12011SSSSHSG60图721617离散系统的频率分析离散时间系统的频率分析也可以调用相应的MATLAB控制系统工具箱函数来完成，这些函数是以连续系统的函数名前加一字母D来命名的，例如离散时间系统的BODE图可以由DBODE函数求出，离散时间系统的NYQUIST图可以由DNYQUIST函数来求出，以及离散时间系统NICHOLS图可以由DNICHOLS函数来求出，其实在MATLAB的控制系统工具箱中这样的函数命名方式是相当普遍的，它们的调用格式与连续系统类似。628时间延迟系统的频率分析函数的调用格式为NUM,DENPADE,N9利用频率响应数据辩识系统模型该函数的调用格式为NUM,DENINVFREQSF,M,N或NUM,DENINVFREQSMAGEXPSQRT1PHASE,M,N6376系统的能控性和能观测性分析761系统的能控性和能观测性