matlab线性控制系统的计算机辅助分析

1、1第第4 4章章 线性控制系统的线性控制系统的计算机辅助设计计算机辅助设计2 早期的控制系统分析过程复杂而耗时，如想得到一个系统的冲激响应曲线，首先需要编写一个求解微分方程的子程序，然后将已经获得的系统模型输入计算机，通过计算机的运算获得冲激响应的响应数据，然后再编写一个绘图程序，将数据绘制成可供工程分析的响应曲线。 MATLAB控制系统工具箱和SIMULINK辅助环境的出现，给控制系统分析带来了福音。 控制系统的分析包括系统的稳定性分析、时域分析、频域分析及根轨迹分析。34.1 线性系统的定性分析对于连续时间系统，4.1.1 线性系统稳定性分析要使x(t)有界，则要eAt有界，即A矩阵的所有

2、特征根均有负实部。故如果闭环极点全部在S平面左半平面，则系统是稳定。tttAttAdBuetxetx00)()()()(0)(在有界信号u(t)的激励下，其状态变量的解析解为4 对于离散时间系统，)()()()()() 1(kTDukTCxkTykTGukTFxTkx其状态变量的解析解为)()0()(101kTGuFxFkTxkiikk要使x(kT)有界，则要Fk有界，即F矩阵的所有特征根的模均小于1。故如果系统全部极点都位于Z平面的单位圆内，则系统是稳定的。 若连续时间系统的全部零极点都位于S左半平面；或若离散时间系统的全部零极点都位于Z平面单位圆内，则系统是最小相位系统。MATLAB中，可

3、用eig( )函数直接求取系统的特征根，也可用pzmap( )函数直接绘制系统的零极点。5例1：系统传函为4032010958411812467284224494536546364032018576022208812266436380598251418)(2345678234567ssssssssssssssssG试判断其稳定性。num=18 514 5982 36380 122664 222088 185760 40320;den=1 36 546 4536 22449 67284 118124 109584 40320;G=tf(num,den);eig(G)ans = -8.0000 -

4、7.0000 -6.0000 -5.0000 -4.0000 -3.0000 -2.0000 -1.0000可见，系统是稳定的。6例2：离散系统受控对象的传函为)951229425. 0)(535261429. 0)(1()71390672. 04040929. 3(00147635. 0)(2zzzzzzH控制器模型为 试分析单位负反馈下的闭环系统稳定性。z=tf(z,0.1);G=0.00147635*(z2+3.4040929*z+ 0.71390672)/(z-1)*(z-0.535261429)*(z-0.951229425);Gc=1.5*(z-0.5)/(z+0.8);GG=fe

5、edback(G*Gc,1);eig(GG) abs(eig(GG)ans = -0.7991 0.7991 0.9745 + 0.0782i 0.9776 0.9745 - 0.0782i 0.9776 0.5344 0.5344 )8 . 0/()5 . 0(5 . 1)(zzzGc7pzmap(GG)84.1.2 线性系统的线性相似变换 相似变换是对方阵进行的一种运算，它把A变成P-1AP，而可逆矩阵P称为进行这一变换的相似变换矩阵。 这种变换的重要意义在于简化对矩阵的各种运算，其方法是先通过相似变换，将矩阵变成与之等价的对角矩阵，再对对角矩阵进行运算，从而将比较复杂的矩阵的运算转化为比

6、较简单的对角矩阵的运算。 MATLAB控制系统工具箱提供了ss2ss( )函数完成状态方程模型的相似变换： G1=ss2ss(G, T) 其中，G为原始的状态方程模型，T为变换矩阵。9., 的特征向量的对应于特征值就是的列向量而的特征值是若iiiApPA163053064A设例3：110101102P令.200010001 1APP则有10例4：设系统的状态方程为实际应用中，若不要求将A变换为对角阵，则P也可用任意非奇异矩阵。 )(01724)()(1000)(10355024100001000010)(txtytutxtx变换矩阵P为反对角矩阵，反对角线上的元素为1，其余元素为0。另：首先介

7、绍fliplr( )函数，其变换矩阵行元素的左右顺序。如A = 1 2 3 4 5 6 7 8 9 10 11 12fliplr(A)= 4 3 2 1 8 7 6 5 12 11 10 911A=0 1 0 0;0 0 1 0;0 0 0 1;-24 -50 -35 -10;G1=ss(A,0;0;0;1,24 7 1 0,0);P=fliplr(eye(4);G2=ss2ss(G1,P)a = x1 x2 x3 x4 x1 -10 -35 -50 -24 x2 1 0 0 0 x3 0 1 0 0 x4 0 0 1 0b = u1 x1 1 x2 0 x3 0 x4 0c = x1 x2

8、x3 x4 y1 0 1 7 24d = u1 y1 0124.1.3 线性系统的可控性分析 现代控制理论用状态方程和输出方程描述系统，输出和输入构成系统的外部变量，而状态为系统的内部变量，系统就好比是一块集成电路芯片，内部结构可能十分复杂，物理量很多，而外部只有少数几个引脚，对电路内部物理量的控制和观测都只能通过这为数不多的几个引脚进行。这就存在着系统内的所有状态是否都受输入控制和所有状态是否都可以从输出反映出来的问题。如果系统所有状态变量的运动都可以通过有限的控制点的输入来使其由任意的初态达到任意设定的终态，则称系统是可控的，更确切的说是状态可控的；如果系统中所有的状态都是可控的，则称该系

9、统为完全可控的系统。否则，就称系统不可控。131. 线性系统的可控性判定 可通过构造可控性判定矩阵 若Tc为满秩矩阵，则系统为完全可控的。如果该矩阵不是满秩矩阵，则它的秩为系统的可控状态的个数。 可控性判定矩阵由Tc=ctrb(A, B)函数构造。 rank( )函数可求出矩阵的秩。例5：,12BABAABBTnc)(48446496)(5 . 311 . 04 . 15 . 029 . 06 . 05 . 163 . 62 . 015 . 17 . 02 . 2) 1(kTukTxTkx试判断系统的可控性。14若系统的Gram矩阵是非奇异矩阵，则该系统是完全可控的。 Gram矩阵为Gram矩

10、阵是以下Lyapunov方程的解求解该Lyapunov方程可用，lyap(A, B*B) 若调用函数不能求出方程的解，则该系统不完全可控。 控制系统的可控Gram矩阵还可以由 Gc=gram(G,c)直接求出。teBBeLtATtAGramTd0TTGramGramBBALAL15例5：已知采样周期为0.1s，试判断系统稳定性。 num=0.1324 -0.5743 0.3879 -0.0889; den=1 -3.233 3.9869 -2.2209 0.4723; G=tf(num,den,Ts,0.1); Lc=gram(ss(G),c);4723. 02209. 29869. 3233

11、. 30889. 03879. 05743. 03124. 0)(23423zzzzzzzzG162. 可控性阶梯分解 对于不完全可控的系统，可对其进行可控性阶梯分解，即构造变换矩阵T，将状态方程变换为这样即将系统的可控子空间，与不可控子空间分离出来。 Ac, Bc, Cc, Tc=ctrbf(A, B, C)函数，可将系统变换为可控性阶梯模型，其中，Tc为相似变换矩阵。例6：请将例5中的系统进行可控性阶梯分解。A=; B=; C =; Ac, Bc, Cc, Tc=ctrbf(A, B, C)173. 可控标准型及其MATLAB实现 若系统完全可控，则可利用矩阵Tc将其变换为第一可控规范型

12、，其系数阵之间满足关系),(111cccCBA111cccccccCCTATATBTB1211111aaaaAnnnc00011cB11BCACABCBCnc其中，18例6：已知系统 的系数阵为),(CBA011,101,400140002CBA试判断其可控性。若完全可控，则求其第一可控规范型。 A=2 0 0;0 4 1;0 0 4; B=1;0;1; C=1 1 0; Tc=ctrb(A,B); rank(Tc)Ac=inv(Tc)*A*Tc; Bc=inv(Tc)*B; Cc=C*Tc;Ac = 0 0 32 1 0 -32 0 1 10Bc = 1 0 0Cc = 1 3 12结果19

13、4.1.4 线性系统的可观性分析 如果系统所有的状态变量任意形式的运动均可由有限测量点的输出完全确定出来，则称系统可观测；如果系统中所有的状态都是可观测的，则称该系统为完全可观测的系统。反之，则称系统不可观测。1. 线性系统的可观测性判定 可通过构造可观测性判定矩阵 若To为满秩矩阵，则系统为完全可观测的。 可观测性判定矩阵由To=obsv(A, C)函数构造。 rank( )函数可求出矩阵的秩。120nCACACACT20A=4 4 4;-11 -12 -12;13 14 13;B=1;-1;0;C=1 1 1; To=obsv(A,C) rank(To)111,011,1314131212

14、11444CBA例7：已知系统的系数阵为试判断它的可观性。To = 1 1 1 6 6 5 23 22 17rank(To)=3故系统可观。结果212. 可观标准型及其MATLAB实现 若系统完全可观，则可利用矩阵P将其变换为第一可观规范型，其系数阵之间满足关系1211111aaaaAnnno其中，111oooCCPPBBPAAPBCACABCBBno11001 1oCoTP122例：上例中，求其第一可控规范型。 P=inv(To); Ao=inv(P)*A*P; Bo=inv(P)*B; Co=C*P;Ao = -0.0000 1.0000 0.0000 0.0000 -0.0000 1.0

15、000 4.0000 -8.0000 5.0000Bo = 0.0000 0.0000 1.0000Co = 1.0000 0.0000 -0.0000结果233. 系统的可观性分解 对于状态不完全可观的系统，同样可对其进行可观性分解。在MATLAB中可调用obsvf( )函数。 AO, BO, CO, T, K=obsvf(A, B, C)其中，T为相似变换阵，K为可观子阵的阶次向量。对系统进行可观性分解后得到相应可观子系统其中，(Ao, Co)为可观子对。24例8：试确定系统的可观性并进行可观性分解。11,01,2112CBAA=-2 1;1 -2; B=1;0; C=1 -1;n=len

16、gth(A); Q=C;C*A; r=rank(Q);if r= =n disp(system is observable.) Qelse disp(system is unobservable.) disp(rank) rend system is unobservable.rankr = 1结果25AO,BO,CO,T,K=obsvf(A,B,C)disp(observable submatrix)AO=AO(2,2)CO=CO(1,2)AO = -1.0000 0 -0.0000 -3.0000BO = 0.7071 -0.7071CO = 0 -1.4142T = 0.7071 0.7

17、071 -0.7071 0.7071K = 1 0observable submatrixAO = CO = -3.0000 -1.4142结果：264.2 控制系统的时域分析 系统的过渡过程性能（如上升时间、调节时间、超调量及稳态误差）常用典型输入作用下的时间响应来描述。响应是指零初始值条件下某种典型的输入函数作用下对象的响应，其响应是时间t的函数，故称为时域响应。控制系统常用的输入函数为单位阶跃函数和脉冲激励函数（即冲激函数）。1. 单位阶跃函数1(t) 0100)( 1tttsesdtettLstst11)( 1)( 1 00X(t)1t272. 单位冲激函数 0t0t0t0st1dte

18、 ) t ()t (Lt)(t)(t0Bttx)(X(t)t 斜坡函数0,0, 0)(tBtttxB=1时称为单位斜坡函数。 其拉氏变换后的像函数为：2)(sBtxL284.2.1 线性定常系统状态方程的解x(0)=x0是系统的初始状态。问题：对给定的控制输入和初始状态，如何确定任意时刻的系统状态和输出；状态的变化行为。首先考虑两边同时左乘e-At,得根据矩阵微积分知识，上式进一步有：BUeXedtdAtAt)(4.1)(4.2)(4.3)(4.4)29两边同时在t0,t区间积分，得两边同时左乘eAt,并整理得即：当初始时刻t0=0时，上式可变为ttAttAtdBUetXe00)()(ttAA

19、tAtdBUetXetXe00)()()(0ttAAtttAdBUeetXetX00)()()(0)(ttdBUttXtttX0)()()()()(00tdBUtXttX0)()()0()()(由此可知，非齐次状态方程的解由两部分组成，第一部分是在初始状态X(t0)作用下的自由运动，第二部分为在系统输入U(t)的作用下的强制运动。(4.5)(4.6)(4.8)(4.9)(4.10)30当U(t)为几种典型的控制输入时，则有如下形式。1.冲激函数输入，即时)()(tktUBkeXeBkdIeXeBkdeeXedBkeeXetXAtAtAtAttAAtAttAAtAt00000000 )( )(

20、)()()()(0BkXetXAt即：BkeXeBkdIeXeBkdeeXedBkeeXetXAtAtAtAttAAtAttAAtAt00000000 )( )( )()(4.11)31BkAIeXeBkAeIeXeBkdeeXeBkdeeXetXAtAtAtAtAttAAtAttAAtAt10100000)( )(2. 阶跃信号输入，即)( 1)(tktU3. 斜坡信号输入，即U(t)=kt，可以求得：BktAIeAXetXAtAt120)()(4.12)(4.13)32例9：求以下系统在单位阶跃函数作用下的状态响应由kkAttAktAAtIe!1! 2122或得ttttttttAteeee

21、eeeee22222222)(11AsILeAt则tttttttttttteeeeeeeeeeeetX2222222121 110012123122212)(4.14)(4.15)33grid; t=0:0.1:10; x1=0.5-exp(-t)+0.5*exp(-2*t); x2=exp(-t)-exp(-2*t);plot(t,x1,x,t,x2,*); xlabel(时间轴)ylabel(x代表x1，-*代表x2)系统状态轨迹图34另： 当输入信号为kitdiitddtddetcu(t)04342)sin()cos(1可将系统原始的状态方程转化为其中，00000000001000000

22、144132ddddBBdBdAA(4.16)(4.17)(4.18)350032DDdDdCC )()()()()()(34321txtxtxtxtxtxxknnnnnT334241,),sin(),cos(11knntdntdnxxtdextdex!01)0()0(10kcccxxkT则，系统状态方程的解析解为以下为MATLAB编写的ss_augment( )，用来求取系统的增广状态方程模型。其中，cc=c0,c1,ck，dd=d1,d2,d3,d4(4.19)(4.20)(4.21)36function Ga,Xa=ss_augment(G,cc,dd,X)G=ss(G); Aa=G.a

23、; Ca=G.c; Xa=X; Ba=G.b; D=G.d;if (length(dd)0 & sum(abs(dd)1e-5), if (abs(dd(4)1e-5), Aa=Aa dd(2)*Ba,dd(3)*Ba;. zeros(2,length(Aa),dd(1), -dd(4);dd(4),dd(1); Ca=Ca dd(2)*D dd(3)*D; Xa=Xa;1;0;Ba=Ba;0;0; else Aa=Aa dd(2)*B;zeros(1,length(Aa) dd(1); Ca=Ca dd(2)*D;Xa=Xa;1;Ba=B;0; endend37if (length(

24、cc)0 & sum(abs(cc)1e-5), M=length(cc); Aa=Aa Ba zeros(length(Aa),M-1); zeros(M-1,length(Aa)+1) eye(M-1); zeros(1,length(Aa)+M); Ca=Ca D zeros(1,M-1);Xa=Xa;cc(1);ii=1; for i=2:M ii=ii*i; %阶乘运算 Xa(length(Aa)+i)=cc(i)*ii; end;endGa=ss(Aa,zeros(size(Ca),Ca,D);38例10: 系统的状态方程模型为cc=2;dd=-4 0 2 3;A=-5 2

25、 0 0;0 -4 0 0;-3 2 -4 -1;-3 2 0 -4;x0=1;2;0;1;B=1;2;3;4;C=1 1 1 1;D=0;G=ss(A,B,C,D);Ga,xx0=ss_augment(G,cc,dd,x0)其中状态变量初值xT(0)=1,2,0,1。假设系统的输入信号u(t)=2+2e-4tsin(3t)，则由ss_augment( )得系统的增广状态方程则由下面语句可直接获得生成信号的解析解syms t; y=Ga.c*expm(Ga.a*t)*xx0; latex(y)3911211121.)(nnnnmnmmasasasabsbsbsbsG4.2.2 基于部分分式展开

26、方法求解 控制系统常用到并联系统，这时就要对系统函数进行分解，使其表现为一些基本控制单元的和的形式。如)()()(sUsGsY若系统输出的Laplace式子的分母多项式的根pi都是不重复的，则mmpsrpsrpsrsY2211)(tpmtptpmererersYLty21211)()(若其中的j个根为m重根pj，则该部分的Laplace反变换为tpmmjjjjetrmtrr)!1(1! 11(111(4.22)(4.23)40部分分式展开：num=2,0,9,1;den=1,1,4,4; r,p,k=residue(num,den)44192)(233ssssssG12225. 0225. 0

27、2)(sisiisisGp= 0.0000+2.0000i 0.0000-2.0000i -1.0000k= 2r= 0.0000-0.2500i 0.0000+0.2500i -2.0000结果表达式：41 MATLAB环境中，函数r,p,k=residue(num,den)对两个多项式进行部分展开，以及把传函分解为微分单元的形式。 向量num和den是按s的降幂排列的多项式系数。部分分式展开后，k为展开余项。例11: 考虑系统的传函为2450351024247)(23423ssssssssG系统的输入信号为单位阶跃信号，则系统的输出信号为)24503510(24247)(23423ssss

28、sssssY将输出信号用部分分式展开,则有42 num=1 7 24 24;den=1 10 35 50 24 0;r,p,k=residue(num,den);r,pans = -1.0000 -4.0000 2.0000 -3.0000 -1.0000 -2.0000 -1.0000 -1.0000 1.0000 0故,由上面结果写出系统输入表达式为12)(234tttteeeety此例为系统传函只含有实数极点的情况，下面的例子中系统输出还含有复数极点。)24503510(24247)(23423sssssssssY4318181123)(234ssssssG例12: 考虑系统的传函为 n

29、um=1 3;den=1 2 11 18 18 0;r,p,k=residue(num,den);r,pans = 0.0020 + 0.0255i 0.0000 + 3.0000i 0.0020 - 0.0255i 0.0000 - 3.0000i -0.0853 + 0.0088i -1.0000 + 1.0000i -0.0853 - 0.0088i -1.0000 - 1.0000i 0.1667 0 故1667. 0)0088. 00853. 0()0088. 00853. 0()0255. 0002. 0()0255. 0002. 0()()1()1(33tjtjjtjtejeje

30、jejty44若想获得输出的物理表达形式，)/(tan) 1(,2)(122abiRbaiR则可借助函数pfrac( )function R,P,K=pfrac(num,den)R,P,K=residue(num,den);for i=1:length(R), if abs(imag(P(i)eps a=real(R(i);b=imag(R(i); R(i)=-2*sqrt(a2+b2);R(i+1)=-atan2(a,b); elseif abs(imag(P(i)eps R(i)=real(R(i); endend)1(sin()()()()()(iRdteiRebjaebjacttdjc

31、tdjc其中，45 离散系统的解析解离散系统的输出信号为mmpzrpzrpzrzY1212111)(nmmmnnpprpprpprzYLny)1()1()1()()(2221111则若考虑采样周期T的因素，则系统输入信号的解析解为kTmmmkTkTpprpprpprzYLkTy)1()1()1()()(2221111(4.24)(4.25)(4.26)46D=conv(1 -1/3,conv(1 -1/4,conv(1 1/5,1 -1); %分母N=0 0 conv(1 -1/2,1 0); %分子前面要补足零N=N(end:-1:1);D=D(end:-1:1); %将z变换式逆序排列R,

32、P,K=residue(N,D);R,P,-R./P)5/1)(4/1)(3/1()2/1()(zzzzzG例13：离散系统的传函为假设系统的输入为阶跃信号，其z变换为z/(z-1)，则其输入可有下列语句计算ans = 12.1528 -5.0000 2.4306 35.5556 4.0000 -8.8889 -16.8750 3.0000 5.6250 -0.8333 1.0000 0.83338333. 0)31(6250. 5)41(8889. 8)51(4306. 2)(nnnny可写出47下例为系统的输出含有重根的情况。例14：离散系统的传函为)3/1()2/1()25()(3zzz

33、zG则其阶跃响应的解析解为D=conv(1 -1/2,conv(1 -1/2,conv(1 -1/2,.conv(1,-1/3,1 -1);N=0 0 0 5 -2 0;R,P,K=residue(N(end:-1:1),D(end:-1:1);R,Pans = 324.0000 3.0000 -240.0000 2.0000 -96.0000 2.0000 192.0000 2.0000 -36.0000 1.0000于是有136)2(192)2(9622403324)(1312111zzzzzzY48 时间延迟系统的解析解 对于带有时间延迟的系统G(s)e-Ls和H(z)z-k，若将其直接

34、进行部分分式展开比较困难。可先不考虑时间延迟，得出系统的解析解y(t)和y(n)，然后分别用t-L或n-k代替得出的解析解中的t或n，即可得到时间延迟系统的解析解。例15：带有时间延迟的系统535)3/1()2/1()25()(zzzzzzG36)726012()21()31(108 36)2)(1()21()2(2/192 .) 1()21()2(96)21(2240)31(3324)(232nnnnnnynnnnnn将上例的解的n替换为n-5，即可494.2.3 二阶系统的阶跃响应及阶跃响应指标 假设系统的开环模型为 ，由单位负反馈构造出的闭环系统模型为)2()(2nnosssG2222)

35、(nnnsssG于是，二阶系统的阶跃响应y(t)为)(21)()()(2dntdntdndndneety其中，nd21根据的取值不同，系统输入也不同。50s2s1 当=0，y(t)=1-cos(nt)为无阻尼振荡。 当01，tnnetty)1 (1)()11(121)(2)1(2)1(222ttneety称为过阻尼振荡。1S1,2S1S2152wn=1;yy=;t=0:.1:12;zet=0:0.1:1,2,3,5;for z=zet if z= =0,y=1-cos(wn*t); elseif (z0 & z1, dd=sqrt(z2-1);lam1=-z-dd;lam2=-z+dd

36、; y=1-0.5*wn*(exp(lam1*t)/lam1-exp(lam2*t)/lam2)/dd; end yy=yy;y;endplot(t,yy)当n=1rad/sec，选择不同的，由下列命令可以得到系统在不同阻尼比下的阶跃响应曲线。53不同阻尼比下的阶跃响应54dt)(h)(9 .0h)(5 .0h)(1 .0hrtptst)(h稳态值稳态值 线性系统典型的阶跃响应曲线，如下图所示误差带：误差带：5，2（1）延迟时间）延迟时间dt（2）上升时间）上升时间rt（3）峰值时间）峰值时间pt（4）st其中，系统的稳态值可由dcgain( )函数获得。（5）%100)()()(%hhthp

37、554.3 线性系统的数字仿真分析4.3.1 线性系统的时域响应1、求取系统单位阶跃响应 step( ) 格式1： Y, T=step(G)格式2： Y,T=step(G,t)格式3： Y, T=step(G,iu)格式4： Y,T=step(G,iu,t)说明：G为tf( ),zpk( ),ss( )中任一种模型。如果用户在调用step( )函数时不返回任何向量，则将自动地绘出阶跃响应输出曲线。对于带有返回参数的将不绘制曲线,其中Y是输出向量, T是时间向量 。t为用户设定的时间向量,一般可以由t=0:step:end等步长地产生出来。对于MIMO系统,iu表示第iu个输入到所有输出的冲激响

38、应曲线. 56num=7 7;den=conv(conv(1 0,1 3),1 4 5);g=tf(num,den);gg=feedback(g,1,-1);y,t,x=step(gg);plot(t,y)例15 已知系统框图)54)(3() 1(7)(2ssssssG其中，输入以下MATLAB命令，可得系统的单位阶跃响应曲线如下图所示。57图4.1 阶跃响应指标显示58例16：设连续系统的数学模型为sessG5 . 03) 1(1)(选择采用周期T=0.1s，试比较各种离散方法的时域响应的差异。G=tf(1,1 3 3 1,iodelay,.5); G1=c2d(G,0.1,zoh);G2=

39、c2d(G,0.1,foh); G3=c2d(G,0.1,tustin);step(G,-r,G1,-g,G2,:b,G3,-.y)59602、 impulse 求连续系统的单位脉冲响应。 格式1： Y, T=impulse(G)格式2： Y,T=impulse(G, t)格式3： Y, T=impulse(G, iu)格式4： Y,T=impulse(G, iu, t)说 明： impulse ( )函数与step( )函数调用格式完全一致。例17:系统sejsjssssssssG2)1)(1)(5)(4)(5 . 3()3)(2)(1(8)(则其脉冲响应曲线由以下语句获得G=zpk(-1;

40、-2;-3,-1+i;-1-i;-3.5;-4;-5,8,iodelay,2);impulse(G,8)61012345678-Impulse ResponseTime (sec)Amplitude624.3.2 任意输入下系统的响应193. 2318. 0109. 1361. 03378. 0107. 2924. 0148. 478. 11134. 0)(29. 123 . 0272. 0sessesssesGssslsim( )函数求任意输入信号时系统的响应 格式1： Y,T=lsim(sys1,u,t) 格式2： Y,T=lsim(

41、sys2,u,t,x0) 说明： u为输入信号.t为等间隔时间向量. sys1为tf( )或zpk( )模型。 sys2为ss( )模型,x0为初始条件。例18：系统其两路输入分别为)2cos()sin()(),13sin(1)(21tttutetut63MATLAB语句为： g11=tf(0.1134,1.78,4.48,1);g12=tf(0.924,2.07,1);g21=tf(0.3378,0.361,1.09,1);g22=tf(-0.318,2.93,1);G=g11,g12;g21,g22;G.ioDelay=0.72,0;0.3,1.29;t=0:.1:15;u=1-exp(-

42、t).*sin(3*t+1),sin(t).*cos(t+2)lsim(G,u,t);64-1-0.500.511.5To: Out(1)051015-1-0.500.511.5To: Out(2)Linear Simulation ResultsTime (sec)Amplitude在两路输入下，系统的时域响应曲线为：6538440010420200)(234sssssG系统的开环传函为综合应用 （1）求单位负反馈时，其闭环传递函数，并绘制输出量阶跃响应曲线和脉冲响应曲线。（2）选择函数的状态变量将其化为状态方程模型，并绘制状态变量的阶跃响应曲线和脉冲响应曲线。close all; clea

43、r;numo=0 0 0 0 200; deno=1 20 140 400 384; %系统开环传函numc=numo; n=length(deno);denc=zeros(1,n);denc=numo+deno; %系统闭环传函disp(System Closed Loop Transfer Function is:)numcdenc66t=0:.05:3;y=step(numc,denc,t); %系统输出量的阶跃响应yy=impulse(numc,denc,t); %系统输出量的脉冲响应plot(t,y);title(System Step Response);xlabel(Time-s

44、ec);ylabel(Response-value);grid;pauseplot(t,yy);title(system impulse respone);xlabel(time-sec);ylabel(respone-value);grid;pauseA,B,C,D=tf2ss(numc,denc);disp(system state-space model is:);67A， B， C， Dys,x=step(A,B,C,D,1,t); %系统状态变量的阶跃响应yys,xx=impulse(A,B,C,D,1,t); %系统状态变量的阶跃响应plot(t,x(:,1),t,x(:,2),+

45、,t,x(:,3),.-,t,x(:,4),.)title(system state-variables step response);xlabel(time-sec);ylabel(response-value);grid;pauseplot(t,xx(:,1),t,xx(:,2),+,t,xx(:,3),.-,t,xx(:,4),.)title(system state-variables impulse response);xlabel(time-sec);ylabel(response-value);grid;68Notice:仿真时间t的选择：对于典型二阶系统根据其响应时间的估算公式

46、 可以确定。对于高阶系统往往其响应时间很难估计，一般采用试探的方法，把t选大一些，看看响应曲线的结果，最后再确定其合适的仿真时间。一般来说，先不指定仿真时间，由MATLAB自己确定，然后根据结果，最后确定合适的仿真时间。在指定仿真时间时，步长的不同会影响到输出曲线的光滑程度，一般不易取太大。nswt43694.4 控制系统的根轨迹分析方法 所谓根轨迹(Root Locus)是指，当系统的某个参数(如开环增益K)由零连续变化到无穷大时, 闭环特征根在复平面上运动的轨迹。一般来说，在无零极点对消时，闭环系统特征根就是闭环传递函数的极点。 根轨迹法(Root Locus Technique)是分析和

47、设计线性定常控制系统的图解方法，使用十分简便。利用它可以对系统进行各种性能分析。一、根轨迹分析方法的概念70R(s) 150(s.sKC(s)图 4-1 控制系统框图 将图4-1所示系统的开环传递函数转化为 )2() 15 . 0()(sskssKsG其中，k=2K，式(4.27)便是绘制根轨迹所用的开环传递函数的标准形式。(4.27) 下面结合图4-1所示的二阶系统的例子, 介绍有关根轨迹的基本概念。71 由式(4.27)可得两开环极点分别为 p1=0, p2= -2，无开环零点。将这两个开环极点绘于图4-2上，并用“”表示。0)(1sG即 02)(2ksssD 由式(4.27)可得闭环系统

48、的特征方程为 得闭环系统的特征根(闭环极点)为 ksks11,1121 从上式看出，闭环系统极点s1, s2与标准化参数k存在一定的关系，参数k变化会引起闭环系统极点s1, s2发生变化。他们之间的关系可由图4-2表示。从图4-2可以看出：(1) 当k=0时，s1， s2与p1， p2重合，即开环极点和闭环极点重合;72(2) 当0k1时，s1， s2均为(-2, 0)区间内的负实数; (3) 当k=1时，s1 = s2= -1， 即两闭环极点重合; (4) 当1k时， ，即两闭环极点互为共轭； (5) 当k时，s1， s2将沿着直线= -1趋于无穷远处。jksjks11,1121k k k

49、1k 0k 0120k 2k 3k 2k 3js二阶系统根轨迹图4-273 从图4-2可以看出: 无论k取何值, 由图4-2表示的控制系统的闭环极点均位于复平面的左半平面, 因此系统是闭环稳定的; 而k=1(K=0.5)时，此二阶系统由过阻尼状态过渡到欠阻尼状态的分界点。并且从图中可以看出, 根轨迹是连续且对称于实轴的, 这也是根轨迹的一个特性。 需要指出的是, 绘制根轨迹时选择的可变参数可以是系统的任何参量, 但实际中最常用的是系统的开环增益。以系统的开环增益为可变参数绘制的根轨迹常规根轨迹. 根据根轨迹，可得以下系统性能： 稳定性 当开环增益K从零到无穷大变化时，图中的根轨迹不会越过虚轴进

50、入s右半平面，因此这个系统对所有的K值都是稳定的。如果根轨迹越过虚轴进入右半s平面，则其交点的K值就是临界稳定开环增益。74 稳态性能开环系统在坐标原点有一个极点，因此根轨迹上的K值就是静态速度误差系数，如果给定系统的稳态误差要求，则可由根轨迹确定闭环极点容许的范围。 动态性能当0K0.5时，闭环极点为复数极点，系统为欠阻尼系统，单位阶跃响应为阻尼振荡过程，且超调量与K成正比。75二、根轨迹分析函数 通常来说，绘制系统的根轨迹是很繁琐的事情，因此在MATLAB中，专门提供了绘制根轨迹的有关函数。rlocus：求系统根轨迹。rlocfind：计算给定一组根的根轨迹增益。sgrid：在连续系统根轨

51、迹图和零极点图中绘制出阻尼系数和自然频率栅格。761、根轨迹图绘制rlocus( )函数 rlocus(a,b,c,d)或者rlocus(num,den)：根据SISO开环系统的状态空间描述模型和传递函数模型，直接在屏幕上绘制出系统的根轨迹图。开环增益的值从零到无穷大变化。 rlocus(a,b,c,d,k)或rlocus(num,den,k)： 通过指定开环增益k的变化范围来绘制系统的根轨迹图。 r=rlocus(num,den,k) 或者r,k=rlocus(num,den) ：不在屏幕上直接绘出系统的根轨迹图，而根据开环增益变化矢量k ，返回闭环系统特征方程1k*num(s)/den(s

52、)=0的根r，它有length(k)行，length(den)-1列，每行对应某个k值时的所有闭环极点。或者同时返回k与r。772、rlocfind()函数 MATLAB提供了函数rlocfind( )来找出给定的一组根（闭环极点）对应的根轨迹增益。其用法如下：k,p=rlocfind(a,b,c,d)或者k,p=rlocfind(num,den)它要求在屏幕上已经绘制好有关的根轨迹图。然后，此命令将产生一个光标以用来选择希望的闭环极点。命令执行结果：k为对应选择点处根轨迹开环增益；p为对应K的系统闭环特征根。 不带输出参数项k,p时，同样可以执行，只是此时只将k的值返回到缺省变量ans中。7

53、83、sgrid( )函数 sgrid：在现存的屏幕根轨迹或零极点图上绘制出自然振荡频率wn和阻尼比矢量z对应的格线。 sgrid(new)：是先清屏，再画格线。 sgrid(z,wn)：绘制由用户指定的阻尼比矢量z、自然振荡频率wn的格线。79三、三、根轨迹分析应用实例 例 19 设一单位负反馈系统的开环传递函数如下：) 3)(2() 1()(ssssKsG试绘制该系统的根轨迹。解解 使用MATLAB绘制此根轨迹的程序如下：num=1 1;den=conv(1 0,conv(1 2,1 3);G=tf(num,den);rlocus(G);title(); xlabel(Re); ylabe

54、l(Im);根轨迹曲线如图4-3所示。80图 4-3 例 19 的MATLAB仿真结果 81例例 20 设离散系统的开环传递函数为 4723. 02209. 29869. 3233. 30889. 03879. 05743. 03124. 0)(23423zzzzzzzzG系统的采样周期为T=0.1s，试画出系统的根轨迹图。 num=0.3124 -0.5743 0.3879 -0.0889;den=1 -3.233 3.9869 -2.2209 0.4723;G=tf(num,den,Ts,0.1);rlocus(G)sgrid另：如图4-4所示，单击左侧与单位圆相交的点，可以得到系统的临界

55、增益值为8。则当k8，闭环系统的全部极点均在单位圆内。82图 4-4 例 20 根轨迹图 83若系统的传函不变，但增加了6步纯滞后，则G.ioDelay=6; rlocus(G),grid,84图4-5由图4-5可见，增加滞后环节后，系统的稳定范围将缩小。85)65)(9 . 0s8 . 1(045. 005. 0)(22sssssG例21：绘制如下系统的根轨迹图请在根轨迹图中用rlocfind函数在图中选择极点位置，得到反馈增益的值。G=tf(0.05,0.045,conv(1,-1.8,0.9,1,5,6)rlocus(G)K=rlocfind(G) %执行后出现下面信息，等待选择极点位置

56、Select a point in the graphics windowselected_point= -3.0166 - 0.0186iK = 3.6981864.5 线性系统的频域分析 频率响应是指系统对正弦输入信号的稳态响应，从频率响应中可以得出带宽、增益、转折频率、闭环稳定性等系统特征。 频率特性是指系统在正弦信号作用下，稳态输出与输入之比对频率的关系特性。频率特性函数与传递函数有直接的关系，记为：为相频特性为幅频特性其中)()()()()()()()()()()(wwwwXwXwAewAjwXjwXjwGioiowjio若以为横轴，A()为纵轴，则可构造幅频特性图；若以为横轴， 为

57、纵轴，则可构造相频特性图。)(87q求取系统对数频率特性图（bode图）：bode( )q求取系统nyquist图（幅相曲线图或极坐标图）：nyquist( ) 频域分析法是应用频率特性研究控制系统的一种典型方法。采用这种方法可直观地表达出系统的频率特性（相频和幅频），分析方法比较简单，物理概念比较明确，对于诸如防止结构谐振、抑制噪声、改善系统稳定性和暂态性能等问题，都可以从系统的频率特性上明确地看出其物理实质和解决途经。通常将频率特性用曲线的形式进行表示，包括对数频率特性曲线和幅相频率特性曲线（简称幅相曲线），MATLAB提供了绘制这两种曲线的函数。884.5.1 单变量系统的频域分析 若将

58、s=j代入线性定常系统的传递函数G(s) ，则有)()()(jQPjG用P()为横轴，Q()为纵轴，将G(j)在复平面表示出来即为Nyquist图。以下为MATLAB中nyquist( )函数格式： nyquist(a,b,c,d)：绘制出系统的一组Nyquist曲线，每条曲线相应于连续状态空间系统a,b,c,d的输入/输出组合对。其中频率范围由函数自动选取，而且在响应快速变化的位置会自动采用更多取样点。 nyquist(a,b,c,d,iu)：可得到从系统第iu个输入到所有输出的极坐标图。一、Nyquist图89 nyquist(num,den)：可绘制出以连续时间多项式传递函数表示的系统的

59、极坐标图。 nyquist(a,b,c,d,iu,w)或nyquist(num,den,w)：可利用指定的角频率矢量w绘制出系统的极坐标图。 当不带返回参数时，直接在屏幕上绘制出系统的极坐标图（图上用箭头表示w的变化方向，负无穷到正无穷） 。当带输出变量re,im,w引用函数时，可得到系统频率特性函数的实部re和虚部im及角频率点w矢量（为正的部分）。可以用plot(re,im)绘制出对应w从负无穷到零变化的部分。90 通过系统的开环频率特性函数Nyquist图，可以判断系统的稳定性，这就是著名的Nyquist判据： 1. 若系统开环稳定（开环传递函数没有极点在s右半平面），则闭环系统稳定（闭

60、环传递函数没有极点在s右半平面）的条件是开环频率特性函数G(j )轨迹不包围（-1, j0）点。 2. 若系统开环不稳定（开环传递函数有n个极点在s右半平面），则闭环系统稳定的条件是，当从-连续变化到+时，开环频率特性函数G(j )的Nyquist图逆时针方向包围复平面上的（-1, j0)点n圈。91例21：有二阶系统 32152)(22sssssG现要得到系统的Nyquist曲线，可输入 -1-0.500.511.522.53-2-1.5-1-0.500.511.52Nyquist PlotReal AxisImaginary Axisnum=2 5 1;den=1 2 3;nyquist(num,den);ti