1、1,本章内容(1) 利用MATLAB求取系统的各种特性函数；(2) 利用MATLAB分析系统的稳定性；(3) 利用MATLAB求取系统在典型和任意输入信号作用下的时域响应；(4) 利用MATLAB绘制系统的根轨迹，在根轨迹上可确定任意点的根轨迹增益K值，从而得到系统稳定的根轨迹增益K值范围；(5) 利用MATLAB绘制系统的Bode图、Nichols图和Nyquist图等，并求取系统的幅值裕量和相位裕量；(6) 利用MATLAB分析具有时间延迟系统的频率特性；(7) 求取频率响应数据，且根据频率响应数据辩识系统的模型参数(8) 利用MATLAB分析系统的能控性和能观测性，并能对不完全能控或不完
2、全能观测的系统进行结构分解。,第7章 控制系统的计算机辅助分析,2,7.1.1 矩阵函数,1 矩阵行列式矩阵Aaij的行列式定义为 Adet(A)=(-1)ka1k1 a2k2anknMATLAB求矩阵行列式函数的调用格式为det(A) 计算矩阵的行列式有多种算法,在MATLAB 中采用的方法为LU分解法。,7.1 系统的特性函数,3,2 矩阵求逆对于一个已知的nn维非奇异方阵A来说，如果有一个同样大小的C矩阵满足ACCAI 式中 I为单位阵，则称C矩阵为A矩阵的逆矩阵，并记作CA1 。MATLAB提供了一个求取逆矩阵的函数inv( )，其调用格式为inv(A),4,3 矩阵的迹假设一个方阵为
3、 Aaij,i,j=1,2,n;则矩阵A的迹定义为亦即矩阵的迹为该矩阵对角线上各个元素之和。由代数理论可知矩阵的迹和该矩阵的特征值之和是相同的。在MATLAB中提供了求取矩阵迹的函数trace( )，其调用方法为trace(A),5,4 矩阵的秩对于nm维的矩阵A，若矩阵所有的列向量中共有 rc个线性无关，则称矩阵的列秩为rc，如果rc=m， 则称A为列满秩矩阵，相应地，若矩阵A的行向量中有rr个是线性无关的，则称矩阵A的行秩为rr, 如果rrn，则称A为行满秩矩阵。 MATLAB提供了一个内部函数rank( )来用数值方法求取一个已知矩阵的秩，其调用格式为k=rank(A),6,5 矩阵的三
4、角分解矩阵的三角分解又称为LU分解，它的目的是将一个矩阵A分解成一个下三角矩阵L和一个上三角矩阵U的乘积，亦即可以写成ALU。在MATLAB下也给出了矩阵的LU分解函数lu( )，该函数的调用格式为L,U=lu(A),7,6 矩阵的奇异值分解cond(A) 7 矩阵的范数N=norm(A,选项),表7-2 矩阵范数函数的选项定义,8,8 矩阵的特征值与特征向量V，D=eig(A) 其中 A为要处理的矩阵，D为一个对角矩阵，其对角线上的元素为矩阵A的特征值，而每个特征值对应的V矩阵的列为该特征值的特征向量，该矩阵是一个满秩矩阵，它满足AVVD，且每个特征向量各元素的平方和（即2范数）均为1。 9
5、 矩阵的特征多项式、特征方程和特征根MATLAB提供了求取矩阵特征多项式系数的函数poly( ），其调用格式为P=poly(A),9,其中 A为给定的矩阵，返回值P为一个行向量，其各个分量为矩阵A的降幂排列的特征多项式系数。即 P= a 0 a 1 a nMATLAB中根据矩阵特征多项式求特征根的函数为roots( )，其调用格式为V=roots(P) 其中 P为特征多项式的系数向量，而V为特征多项式的解，即原始矩阵的特征根。,10,7.1.2 数值分解,1. 数据处理 如果给定一组数据xi，i=1,2,n, 则可利用MATLAB将这些数据用一个向量表示出来，即 x=x1,x2,xn 利用MA
6、TLAB的函数max( )和min( )便可求出这组数据的最大和最小值，命令格式如下 xM,i=max(x)xm,i=min(x),11,2. 按实部或幅值对特征值进行排序MATLAB的函数esort( )和dsort( )的调用格式为s,ndx=esort(P) 或 s,ndx=dsort(P) 其中 esort(P)对连续系统，根据实部按递减顺序对矢量P中的复特征值进行排序，ndx为索引矢量,12,3. 多项式及多项式矩阵求值 如果f(x)函数为下面的一个多项式f(x)=a0xn+ a1x(n-1)+an-1x+an 则可以求出该函数的导数函数为f (x)=n a0 xn-1+(n-1)a
7、1xn-2+an-1在MATLAB中提供了多项式求值函数polyval( ) 和多项式求导的函数polyder( )它们的调用格式分别为f(x0)=polyval(p,x0) 及 p1=polyder(p),13,4. 数值积分 对于函数f(x)的定积分可利用MATLABquad( )函数来求取定积分，调用格式为y=quad(函数名,a,b,tol),14,5. 非线性方程求解 对于非线性方程 f(x,t)=0 MATLAB提供了函数fzero( )和fsolve( )可以方便地求得解，它们的调用格式分别为x=fzero(函数名，初值) 和 x=fsolve(函数名，初值),15,6. 函数的
8、极值MATLAB中提供了两个基于单纯形算法求解函数最小值的函数fmin( )和fmins( )，它们分别对应于单变量和多变量的最优化问题的求解，其调用格式分别为x=fmin(函数名，初值，选项)x=fmins(函数名，初值，选项),16,7. 函数图形的绘制 MATLAB还允许调用fplot( )函数来直接绘制出函数的图形，该函数的调用格式为 fplot(函数名,t0,tf,N) 其中 t0,tf为绘制区间,N为点数，N默认时将取N25，函数名既可为自定义的任意m函数，也可为基本数学函数。,17,7.1.3 求系统的阻尼系数和自然频率,MATLAB函数damp( )和ddamp( ) 的调用格
9、式为n,=damp(A)n, p=ddamp(A)n, p=ddamp(A,Ts),18,7.1.4 求控制系统的增益,MATLAB的函数dcgain ( )的调用格式为K=dcgain(num,den) 或 K=dcgain(A,B,C,D) 其中 返回值K为系统增益，即放大系数。对于单变量系统，K为标量；对于多变量系统，K为向量或矩阵。,19,3.显示打印线性系统MATLAB的printsys( )函数可按特殊格式打印出状态空间和传递函数表示的系统。 printsys(num,den,s) %显示打印连续系统的传递函数，默认方式； printsys(num,den,z) %显示打印离散系统
10、的脉冲传递函数； printsys(A,B,C,D) %显示打印状态空间形式的系数矩阵； printsys(A,B,C,D,ulabels,ylabels,xlabds) %用ulabels,ylabels,xlabds中指定的符号标记出系统矩阵A,B,C,D。,20,7.2 控制系统的稳定性分析,1.利用极点判断系统的稳定性判断一个线性系统稳定性的一种最有效的方法是直接求出系统所有的极点，然后根据极点的分布情况来确定系统的稳定性，对于极点的求取我们在上节中已作过介绍，下面举例说明其判断方法。 2.利用特征值判断系统的稳定性系统的特征方程 |sI-A|=sn+a1sn-1+an-1s+an0
11、的根称为系统的特征值，即系统的闭环极点。当然判断系统的稳定性同样可利用特征值来判断。,21,3.用李雅普诺夫第二法来判断系统的稳定性,线性定常连续系统(7-2)在平衡状态xe=0处，渐近稳定的充要条件是：对任给的一个正定对称矩阵Q，存在一个正定的对称矩阵P，且满足矩阵方程ATPPAQ 而标量函数V(x)=xTPx是这个系统的一个二次型形式的李雅普诺夫函数。MATLAB提供了李雅普诺夫方程的求解函数lyap( )，其调用格式为P=lyap(A ,Q),22,7.3 控制系统的时域分析,1.任意信号函数 生成任意信号函数gensig( )的调用格式为u,t=gensig(type,Ta) 或 u,
12、t=gensig(type,Ta,Tf,T) 其中 第一式产生一个类型为type的信号序列u(t)，周期为Ta，type为以下标识字符串之一：sin正弦波；square方波；pulse脉冲序列；第二式同时定义信号序列u(t)的持续时间Tf和采样时间T。,23,例7-11 生成一个周期为5秒，持续时间为30秒，采样时间为0.1秒的方波。 解 Matlab窗口中执行以下命令可得图7-3所示结果。 u,t=gensig(square,5,30,0.1); plot(t,u) axis(0,30 0.5 1.5),24,图7-3,25,2.连续系统的单位阶跃响应单位阶跃响应函数step( )的调用格式
13、为y,x,t=step(num,den,t) 或 y,x,t=step(A,B,C,D,iu,t) 式中 t为选定的仿真时间向量，函数返回值y为系统在各个仿真时刻的输出所组成的矩阵；而x为自动选择的状态变量的时间响应数据。如只想绘制出系统的阶跃响应曲线，则可以由如下的格式调用此函数step(num,den,t)step(A,B,C,D,t),26,例7-13 对于典型二阶系统试绘制出无阻尼自然振荡频率n=6，阻尼比分别为0.2,0.4,1.0,2.0时系统的单位阶跃响应曲线。,27,解 MATLAB程序为 %Example7_13.m wn=6;zeta=0.2:0.2:1.0,2.0; fi
14、gure(1);hold on for I=zeta num=wn.2; den=1,2*I*wn,wn.2; step(num,den);end title(Step Response);hold off 执行后可得如图7-5所示的单位阶跃响应曲线。,图7-5,28,3.离散系统的单位阶跃响应离散系统的单位阶跃响应函数dstep( )的调用格式为y,x=dstep(num,den,n) 或 y,x=dstep(G,H,C,D,iu,n) 式中 n为选定的取样点个数，当n省略时，取样点数由函数自动选取，其余参数定义同前。,29,例7-15 已知二阶离散系统试求其单位阶跃响应。 解 MATLAB
15、程序为 %Example7_15.m num=2 -3.4 1.5; den=1 -1.6 0.8; dstep(num,den); title(Piscrete Step Response) 执行后得如图7-7所示的单位阶跃响应曲线。,图7-7,30,4.单位脉冲响应 单位脉冲响应函数impulse( ) 和dimpulse( )与单位阶跃函数step( )和dstep( )的调用格式完全一致，这里就不一一列写了。,31,5.系统的零输入响应对于连续系统由初始状态所引引起的响应，即零输入响应，可由函数initial( )来求得，其调用格式为y,x,t=initial(A,B,C,D,x0)
16、或 y,x,t=initial(A,B,C,D,x0,t) 其中 x0为初始状态，其余参数定义同前。,32,6. 任意输入函数的响应连续系统对任意输入函数的响应可利用MATLAB的函数lsim( )求取，其调用格式为y,x=lsim(num,den, u,t) 或 y,x=lsim(A,B,C,D,iu,u,t) 其中 u为给定输入序列构成的矩阵，它的每列对应一个输入，每行对应一个新的时间点，其行数与时间t的长度相等。其他用法同step( )函数。,33,7.4 根轨迹法,所谓根轨迹是指，当开环系统的某一参数从零变化到无穷大时，闭环系统的特征方程根在s平面上所形成的轨迹。一般地，将这一参数选作
17、开环系统的增益K，而在无零极点对消时，闭环系统特征方程的根就是闭环传递函数的极点。,34,1.绘制系统的零极点图pzmap( )函数可绘制系统的零极点图，其调用格式为p,z=pzmap(A,B,C,D) 或 p,z=pzmap(p,z)，p,z=pzmap(num, den) 其中 列向量P为系统极点位置，列向量z为系统的零点位置。可通过pzmap(p,z)绘制出零极点图，图中的极点用 “”表示，零点用 “o”表示。,35,2.绘制系统的根轨迹,对于图7-12所示的负反馈系统，其特征方程可表示为或 利用rlocus( )函数可绘制出当开环增益K由0至变化时，闭环系统的特征根在s平面变化的轨迹，
18、该函数的调用格式为r,K=rlocus(num,den) r,K=rlocus(num,den, K) 或 r,K=rlocus(A,B,C,D)r,K=rlocus(A,B,C,D, K) 其中 返回值r为系统的闭环极点， K为相应的增益。,36,在系统分析过程中，常常希望确定根轨迹上某一点处的增益值K，这时可利用MATLAB中的rlocfind( )函数，在使用此函数前要首先得到系统的根轨迹，然后再执行如下命令K,poles=rlocfind(num,den) 或 K,poles=rlocfind(A,B,C,D)执行上述命令后，将在屏幕上的图形中生成一个十字光标，使用鼠标移动它至所希望的
19、位置，然后敲击鼠标左键即可得到该极点的位置坐标值poles以及它所对应的增益K值。,37,例7-20 已知某负反馈系统的开环传递函数为试绘制系统根轨迹，并分析系统稳定的K值范围。 解 MATLAB的程序为 %Example7_20.m num=1;den=conv(1,0,conv(1,1,1,2); rlocus(num,den), K,poles=rlocfind(num,den);,38,执行以上程序，并移动鼠标到根轨迹与虚轴的交点处单击鼠标左键后可得如图7-13所示的根轨迹和如下结果,selected_point = 0.0000 - 1.4142i K =6.0000 poles =
20、-3.0000 0.0000 +1.4142i 0.0000 - 1.4142i,图7-13,39,由此可见根轨迹与虚轴交点处的增益K=6，这说明当K 6时，系统不稳定；利用rlocfind( )函数也可找出根轨迹从实轴上的分离点处的增益K =0.38, 这说明当0 K 0.38时，系统为单调衰减稳定，当0.38 K 6时系统为振荡衰减稳定的。,40,7.5 控制系统的频域分析,频率响应研究系统的频率行为，从频率响应中可得带宽、增益、转折频率和闭环系统稳定性等系统特征。MATLAB的控制系统工具箱提供了多种求取线性系统频率响应曲线的函数，如表7-5所示。,41,1. 产生频率向量频率向量可由l
21、ogspace( )函数来构成。此函数的调用格式为 =logspace(m,n,npts) 此命令可生成一个以10为底的指数向量(10m10n )，点数由npts任意选定。,42,2.系统的伯德图（Bode图）,连续系统的伯德图可利用bode( )函数来绘制，其调用格式为mag,phase,=bode(num,den)mag,phase,=bode(num,den,)mag,phase,=bode(A,B,C,D)mag,phase,=bode(A,B,C,D,iu)mag,phase,=bode(A,B,C,D,iu,) 式中 num,den和A,B,C,D分别为系统的开环传递函数和状态方程
22、的参数，而为频率点构成的向量。,43,有了这些数据就可以利用下面的MATLAB命令subplot(2,1,1);semilogx(w,20*log10(mag)subplot(2,1,2);semilogx(w,phase)在同一个窗口上同时绘制出系统的Bode图了，其中前一条命令中对幅值向量mag求取分贝(dB)值。,44,如果只想绘制出系统的Bode图，而对获得幅值和相位的具体数值并不感兴趣，则可以采用如下简单的调用格式bode(num,den,)bode(A,B,C,D,iu,） 或更简单地bode(num,den) bode(A,B,C,D,iu),45,例7-22 已知二阶系统的开环
23、传递函数为绘制出取不同值时, n在0.1到10之间的频率范围上的Bode图。 解 当n=6，取0.2,0.4,0.6,0.8,1.0时二阶系统的Bode图可直接采用bode( )函数得到。MATLAB程序为 %Example7_22.m wn=6;zeta=0.2:0.2:1.0; w=logspace(-1,1);figure(1);num=wn.2;,46,for k=zeta den=1 2*k*wn wn.2; mag,phase, w1=bode(num,den, w); subplot(2,1,1);hold on semilogx(w1, mag); subplot(2,1,2)
24、;hold on semilogx(w1,phase); end subplot(2,1,1);grid on;grid;grid; title(Bode plot) xlabel(Frequency (rad/ses);ylabel(Gain dB); subplot(2,1,2);grid on;grid; xlabel(Frquency (rad/sel);ylabel(Phase deg);hold off,执行后得如图7-15所示Bode图。 从图中可以看出，当0时，相角趋于0，当时，相角趋于 -1800，当=n时，相角等于 -900，此时的幅值也最大。,图7-15,47,3.幅值裕
25、量和相位裕量在判断系统稳定性时，常常需要求出系统的幅值裕量和相位裕量。利用MATLAB控制系统工具箱提供的margin( )函数可以求出系统的幅值裕量与相位裕量，该函数的调用格式为Gm,Pm,Wcg,Wcp=margin(num,den) 或 Gm,Pm,Wcg,Wcp=margin(A,B,C,D) 式中 Gm和Pm分别为系统的幅值裕量和相位裕量，而Wcg 和Wcp分别为幅值裕量和相位裕量处相应的频率值。,48,例7-23 给定系统的开环状态空间表达式为 求系统的幅值裕量和相位裕量，并画出Bode图。,49,%Example7_23.m w=logspace(-1,1); A=0 1 0 0
26、;0 0 1 0;0 0 0 1;-62.5 -213.8 -204.2 -54; B=0;0;0;1; C=1562 1875 0 0;D=0; Gm,Pm,Wcg,Wcp=margin(A,B,C,D) bode(A,B,C,D,1,w); disp(幅值裕量=,num2str(20*log10(Gm),dB,相位裕量=,num2str(Pm),度。); title(Gain margin=,num2str(Gm), Phase margin=num2str(Pm),50,执行后得如下数据及图7-16所示Bode图。,图7-16,51,Gm =4.4922 Pm =23.0705 Wcg
27、=12.6466 Wcp =5.8275 幅值裕量kg=4.4922=13.0492dB，相位裕量r =23.0706度。,52,除了根据系统模型直接求取幅值和相位裕量之外，MATLAB的控制系统工具箱中还提供了由幅值和相位相应数据来求取裕量的方法，这时函数的调用格式为 Gm,Pm,Wcg,Wcp=margin(mag,phase,) 式中 频率响应可以是由bode( )函数获得的幅值和相位向量，也可以是系统的实测幅值与相位向量，为相应的频率点向量。,53,4.频率响应值MATLAB也提供了直接的求取频率响应数据的函数freqresp( ),其调用格式为F=freqresp(num,den,s
28、qrt(-1)*) 或 F=freqresp(A,B,C,D,iu,sqrt(-1)*) 式中 F为频率响应，为给定的频率范围向量。,54,5.系统的奈奎斯特图（Nyquist图）,nyquist( )函数的调用格式为Re,Im,=nyquist(num,den)Re,Im,=nyquist(num,den,)Re,Im,=nyquist(A,B,C,D)Re,Im,=nyquist(A,B,C,D,iu)Re,Im,=nyquist(A,B,C,D,iu,) 其中 返回值Re,Im和分别为频率特性的实部向量、虚部向量和对应的频率向量，有了这些值就可利用命令plot(Re,Im)来直接绘出系统
29、的奈奎斯特图。,55,当然也可使用下面的简单命令来直接绘出系统的奈魁斯特图。nyquist(num,den,) 或 nyquist(A,B,C,D) 更简单地 nyquist(num,den) 或 nyquist(A,B,C,D,iu) 它的使用方法基本同bode( )函数的用法。,56,例7-26 已知系统的开环传递函数为绘制Nyquist图，并判断系统的稳定性。 解 MATLAB命令为 num=0.5;den=1 2 1 0.5; nyquist(num,den) 执行后可得如图7-17所示的曲线，由于Nyquist曲线没有包围(-1,j0)点，且P0，所以由G(s)H(s)构成的单位负反
30、馈闭环系统稳定。,57,图7-17,58,6.尼柯尔斯图(Nichols图)函数调用格式为mag,phase,=nichols(num,den,) 或 mag,phase,=nichols(A,B,C,D,iu,) 可见该函数的调用格式以及返回的值与bode( )函数完全一致，事实上虽然它们使用的算法不同, 但这两个函数得出的结果还是基本一致的。但Nichols图的绘制方式和Bode图是不同的，它可由以下命令绘制 plot(phase,20*log10(mag)当然，Nichols图也可采用与Bode图类似的简单命令来直接绘制。,59,例7-29 已知单位负反馈的开环传递函数为试绘制Nicho
31、ls图。 解 MATLAB程序为 %Example7_29.m num=1;den=conv(1,0,conv(1,1,0.2,1); w=logspace(-1,1,400); mag,phase=nichols(num,den,w); plot(phase,20*log10(mag) ngrid 执行后可得如图7-21所示的Nichols图。,60,图7-21,61,7.离散系统的频率分析 离散时间系统的频率分析也可以调用相应的MATLAB控制系统工具箱函数来完成，这些函数是以连续系统的函数名前加一字母d来命名的，例如离散时间系统的Bode图可以由dbode( )函数求出，离散时间系统的N
32、yquist图可以由dnyquist( )函数来求出，以及离散时间系统Nichols图可以由dnichols( )函数来求出，其实在MATLAB的控制系统工具箱中这样的函数命名方式是相当普遍的，它们的调用格式与连续系统类似。,62,8.时间延迟系统的频率分析 函数的调用格式为num,den=pade(,n) 9.利用频率响应数据辩识系统模型该函数的调用格式为num,den=invfreqs(F,m,n) 或 num,den=invfreqs(mag.*exp(sqrt(-1)*phase),m,n),63,7.6 系统的能控性和能观测性分析,7.6.1 系统的能控性和能观测性对于n阶线性定常系
33、统,64,能控性矩阵为 Uc=B AB A2B An-1B当rank Uc=n时，系统的状态完全能控，否则系统不能控。 能观测性矩阵为当rankVon时，系统的状态完全能观测，否则系统状态不能观测。,65,在MATLAB中，可利用ctrb( )和obsv( )函数直接求出能控性和能观测性矩阵，从而确定系统的状态能控性和能观测性。它们的调用格式分别为 Uc=ctrb(A,B)Vo=obsv(A,C) 其中 A,B,C为系统的各矩阵，Uc和Vo分别为能控性矩阵和能观测性矩阵,66,另外，利用系统的能控性和能观测性Gram矩阵和 是否满秩，也可判别系统的能控性和能观测性。在MATLAB中，求系统的能
34、控和能观测性Gram矩阵的函数gram( )的调用格式为Wc=gram(A,B) 和 Wo=gram(A,C),67,7.6.2 将系统按能控和不能控性进行分解,MATLAB的控制系统工具箱中提供了这种分解的函数ctrbf( )，其调用格式为=ctrbf(A,B,C)其中 T为相似变换阵，K是长度为n的一个矢量，其元素为各个块的秩，sum(K)可求出A中能控部分的秩， 对应与转换后系统的A,B,C。,68,7.6.3 将系统按能观测性和不能观测性进行分解,MATLAB控制系统工具箱中，能观测性分解函数obsvf( )的调用格式为 = obsvf(A,B,C )其中 T为相似变换阵，k是长度为n的一个矢量，其元素为各个块的秩，sum(K)可求出A中能观测部分的秩， 对应与转换后系统的A,B,C。,69,7.7 系统模型的降阶,7.7.1 平衡实现该函数的调用格式为 balreal 这一函数可由给定的状态空间表达式(A,B,C)，直接求出平衡系统的Gram矩阵G和变换矩阵T及平衡实现的模型。,70,7.7.2 模型降阶函数的调用格式为modred( elim) 式中 为平衡实现系统的状态空间表达式模型参数，elim为要消去的状态变量序号，当然为获得较好的近似，消去的状态变量应该选为Gram矩阵元素的数值较小的变量，返回的A1,B1,C1,D1为所获得的降阶模型。,
