第4章线性控制系统的计算机辅助设计

1、本章主要内容线性系统定性分析线性系统时域响应解析解法线性系统的数字仿真分析根轨迹分析线性系统频域分析4.1 4.1 线性系统性质分析-主要内容 线性系统稳定性分析线性反馈系统内部稳定性分析线性系统的相似变换线性系统可控性分析线性系统可观测性分析Kalman分解系统状态方程的标准型系统的范数测度及求解4.1 4.1 线性系统性质分析-稳定性分析(1)(1) ()()()()()x kTFx kTGu kTy kTCx kTDu kT00()()0( )( )( )tA t tA ttx tex teBud( )( )( )( )( )( )x tAx tBu ty tCx tDu t110()(

2、0)()kkk iix kTF xFGu iT 4.1 4.1 线性系统性质分析-稳定性分析(2)基于MATLAB的稳定性直接判定状态方程模型 由eig(A)可求出所有特征根 离散系统：abs(eig(A)传递函数模型：完全同样方法基于MATLAB的图解判定法连续系统：pzmap(G)离散系统：pzmap(G)，同时画出单位圆4.1 4.1 线性系统性质分析-稳定性分析(3)例4-1 单位负反馈构成的闭环系统稳定性判定直接分析方法num=10 50 100 100 40; den=1 21 184 870 2384 3664 2496 0; G=tf(num,den);GG=feedback(

3、G,1); pzmap(GG) eig(GG)用zpk(GG)命令得出零极点模型，从而判断稳定性432765432105010010040( )21184870238436642496ssssG ssssssss4.1 4.1 线性系统性质分析-稳定性分析(4)例4-2 高阶离散单位负反馈系统模型den=1 -1 0.25 0.25 -0.125; num=6 -0.6 -0.12; H=tf(num,den, Ts,0.1); z=tf(z, Ts,0.1); Gc=0.3*(z-0.6)/(z+0.8); GG=feedback(H*Gc,1); pzmap(GG), abs(eig(GG

4、)2432c60.60.12( ),250.6( )0.3,0.10.8zzH zzzzzzG zTz4.1 4.1 线性系统性质分析-线性反馈系统内部稳定性(1)输入、输出稳定是不够的，因为若内部信号可能过大，对系统作硬件破坏应该引入内部稳定性概念，保证内部信号也是稳定的。带有扰动的线性反馈控制系统如下图所示c( )G sud1xr2x( )G sy( )H svn3x4.1 4.1 线性系统性质分析-线性反馈系统内部稳定性(2)从输入信号(r,d,n)到内部输出信号(x1,x2,x3)都稳定的系统称为内部稳定系统传递函数矩阵 其中逐一判断每个子传递函数的稳定性很繁琐1

5、2cc3c1( )( )( )1( )1( )( )( )( )( )( )1xG s H sH srxG sG s H sdM sxG s G sG sn c( )1( )( )( )M sG s G s H s 4.1 4.1 线性系统性质分析-线性反馈系统内部稳定性(3)闭环系统内部稳定的充要条件为 1+H(s)G(s)Gc(s)没有不稳定的零点 H(s)G(s)Gc(s)没有不稳定零极点对消第一个条件等效于输入输出稳定性判断第2个条件即可可以编写MATLAB函数判断内部稳定性 key=intstable(G,Gc,H)4.1 4.1 线性系统性质分析-线性反馈系统内部稳定性(4)判定的

6、MATLAB函数function key=intstable(G,Gc,H) GG=minreal(feedback(G*Gc,H);G0=H*G*Gc;G01=minreal(G0);p=eig(GG);z0=eig(G0);z1=eig(G01);zz=setdiff(z0,z1);if(G.Ts1), %离散系统判定 key=any(abs(p)1); if key=0, key=2*any(abs(zz)1);endelse, %连续系统判断 key=any(real(p)0); if key=0, key=2*any(real(zz)0);endend1:输入输出不稳定；2：输入输出

7、稳定、但内部不稳定；0：内部稳定4.1 4.1 线性系统性质分析-线性相似变换(1)系统的状态方程表示称为系统实现不同状态选择下，状态方程不唯一相似变换非奇异矩阵T状态变换z=T-1x新状态方程模型1( )( )( ),(0)(0)( )( )( )ttttz tAz tBu tzTxy tC z tDu t4.1 4.1 线性系统性质分析-线性相似变换(2)状态变换公式MATLAB求解方法G1=ss2ss(G,T)11,ttttATATBT BCCTDD输入输出都要是状态方程对象4.1 4.1 线性系统性质分析-线性相似变换(3)例4-3 已知系统和转换矩阵A=0 1 0 0;0 0 1 0

8、;0 0 0 1;-24 -50 -35 -10; G1=ss(A,0;0;0;1,24 24 7 1,0); T=fliplr(eye(4); G2=ss2ss(G1,T)0100010010( )010001( )( )012450351011( )242471( )x tx tu tTy tx t 4.1 4.1 线性系统性质分析-线性相似变换(4)变换结果：相似变换能改变系统的结构引入相似变换矩阵可以将已知系统转换成其他的形式10355024110000( )( )( )0100000100( )172424( )z tz tu ty tz t 4.1 4.1 线性系统性质分析-可控性

9、分析(1)可控性定义假设系统由状态方程(A,B,C,D)给出，对任意的初始时刻t0，如果状态空间中任一状态xi(t)可以从初始状态xi(t0)处，由有界的输入信号u(t)的驱动下，在有限时间tf内能够到达任意预先指定的状态xi(tf)，则称此状态是可控的；如果系统中所有的状态都是可控的，则称该系统为完全可控的系统；系统的可控性就是指系统内部的状态是不是可以由外部输出信号控制的性质；4.1 4.1 线性系统性质分析-可控性分析(2)可控性判定矩阵方法 若矩阵Tc满秩，则系统完全可控 基于MATLAB的判定方法： rank(T) 构造可控性判定矩阵： Tc=ctrb(A,B)21 ,ncTB AB

10、 A BAB4.1 4.1 线性系统性质分析-可控性分析(3)例4-4 离散状态方程的可控性A=-2.2 -0.7 1.5 -1;0.2 -6.3 6 -1.5;0.6 -0.9 -2 -0.5;1.4 -0.1 -1 -3.5; B=6 9;4 6;4 4;8 4; Tc=ctrb(A,B); rank(Tc)或者Tc=B,A*B,A2*B,A3*B; %用直接方法建立可控性矩阵 rank(Tc)1690.26.361.546(1) ()()0.60.920.5484x kTx kTu kT4.1 4.1 线性系统性质分析-可控性分析(4)由Gram矩

11、阵的非奇异判定可控性 引入可控Gram矩阵： 该矩阵满足Lyapunov方程 MATLAB求解：Lc=lyap(A,B*B) 可控Gram矩阵还可以由Gc=gram(G,c)直接求出0TAtTA tcLeBB edtTTccALL ABB 4.1 4.1 线性系统性质分析-可控性分析(5)例4-5 求Gram矩阵num=6 -0.6 -0.12;den=1 -1 0.25 0.25 -0.125; H=tf(num,den, Ts,0.1); Lc=gram(ss(H), c)243260.60.12( )25zzH zzzzz4.1 4.1 线性系统性质分析-可控性分析

12、(6)可控性阶梯分解 对于不完全可控的系统进行阶梯分解 阶梯标准型 MATLAB函数调用 Ac,Bc,Cc,Tc=ctrbf(A,B,C) 若原系统状态方程完全可控，则不必分解2100,cccccccAABCCCBAA4.1 4.1 线性系统性质分析-可控性分析(7)例4-6 不完全可控系统A=-2.2 -0.7 1.5 -1;0.2 -6.3 6 -1.5; . 0.6 -0.9 -2 -0.5; 1.4 -0.1 -1 -3.5; B=6 9;4 6;4 4;8 4;C=1 2 3 4; Ac,Bc,Cc,Tc=ctrbf(A,B,C);1690.26.361.546(1

13、) ()()0.60.920.5484x kTx kTu kT4.1 4.1 线性系统性质分析-可观测性分析(1)可观测性定义假设系统由状态方程(A,B,C,D)给出,对任意的初始时刻t0,如果状态空间中任一状态xi(t)在任意有限时刻tf的状态xi(tf)可以由输出信号在这一时间区间内tt0, tf的值精确地确定出来，则称此状态是可观测的。如果系统中所有的状态都是可观测的，则称该系统为完全可观测的系统。系统的可观测性就是指系统内部的状态是不是可以由系统输出信号重建起来的性质。4.1 4.1 线性系统性质分析-可观测性分析(2)可观测性与可控性是对偶问题可观测性判定矩阵方

14、法 系统(A,C)的可观测性问题等同于(AT ,CT) 的可控性问题 MATLAB求解可观测矩阵obsv(), obsvf()21onCCATCACA4.1 4.1 线性系统性质分析-可观测性分析(3)Gram矩阵判定可观测性 MATLAB求解：gram(G, o) Gram矩阵满足Lyapunov方程 0TA tTAtoLeC CedtTTooA LL AC C 4.1 4.1 线性系统性质分析-kalman规范分解(1)Kalman规范分解,1,2,3,13,2,3,4,4,2,0000000( )( )( )00( )00( )c oc oc oc oc oc oc oc oAAAz t

15、z tu tBAAAABAAy tCCz t4.1 4.1 线性系统性质分析-kalman规范分解(2)子空间 为既不可控又不可观测的子空间， 为不可控但可观测的子空间，和 分别为可控但不可观测的子空间和既可控又可观测的子空间。既可控又可观的子空间就是前面提及的最小实现模型,(,0,0)c oA,(,0,)c oc oAC,(,0)c oc oAB,(,)c oc oc oABC4.1 4.1 线性系统性质分析-系统状态方程标准型的MATLAB求解(1)常用标准型单变量系统的标准型可控标准型实现可观测标准型实现Jordan标准型实现多变量系统Leunberge标准型侧重点：如何用MATLAB直

16、接获取标准型4.1 4.1 线性系统性质分析-系统状态方程标准型的MATLAB求解(2)单变量系统的标准型 可控标准型 可观测标准型ccccxA xB uyC xDu12120100000000101nnxxuaaaybbb x ooooxA xB uyC xD u1122330001000100010,0,1nnababxaxbuabyx4.1 4.1 线性系统性质分析-系统状态方程标准型的MATLAB求解(3)可控标准型和可观测标准型互为对偶形式 可控可观测标准型MATLAB实现 function Gs=sscanform(G,type) switch type case ctrl G=t

17、f(G);Gs=; G.num1=G.num1/G.den1(1); G.den1=G.den1/G.den1(1); d=G.num1(1); G1=G; G1.ioDelay=0; G1=G1-d; num=G1.num1;den=G1.den1;n=length(G.den1)-1;,TTTTcocococoAABCCBDD4.1 4.1 线性系统性质分析-系统状态方程标准型的MATLAB求解(4) A=zeros(n-1,1) eye(n-1);-den(end:-1:2); B=zeros(n-1,1);1; C=num(end:-1:2); D=d; Gs=ss(A,B,C,D,T

18、s,G.Ts,ioDelay,G.ioDelay); case obsv Gc=sscanform(G,ctrl); Gs=ss(Gc.a,Gc.c,Gc.b,Gc.d,Ts,G.Ts,ioDelay,G.ioDelay); otherwise error(Only options ctrl and obsv are applicable.) end 4.1 4.1 线性系统性质分析-系统状态方程标准型的MATLAB求解(5)Jordan标准型 假设系统矩阵A的特征根为1 ,2 , , n ，第i个特征根i 对应特征向量为v i ，则 矩阵A对应的模态矩阵定义为 MATLAB变换： G1,T=

19、canon(G, modal) ,1,2,iiiAvvin 12kJJJz=Tx4.1 4.1 线性系统性质分析-系统状态方程标准型的MATLAB求解(6)多变量系统的可控标准型实现-Luenberger标准型 可控性判定矩阵 构造矩阵1211111122,ppSb AbAb bAbAb11211TTTllLl提取此行提取此行 : 保证前面各列线性无关的最大指数值，即最大可控性指数is4.1 4.1 线性系统性质分析-系统状态方程标准型的MATLAB求解(7)得出Luenberger变换矩阵编写Luenberger 实现的MATLAB函数: T=luenberger(A,B)111212111

20、TTTll ATlA11ss,ATAT BT B4.1 4.1 线性系统性质分析-系统状态方程标准型的MATLAB求解(8)MATLAB函数清单function T=luenberger (A,B)n=size(A,1); p=size(B,2); S=; sigmas=; k=1;for i=1:p for j=0:n-1 S=S, Aj*B(:,i) ; if rank(S)=k, k=k+1; else, sigmas(i)=j-1; S=S(:,1:end-1); break; end, end if kn, break; endendk=k-1; 4.1 4.1 线性系统性质分析-系

21、统状态方程标准型的MATLAB求解(9)if kn % 如果不是完全可控，则用随机数补足 while rank(S)=n, S(:,k+1:n)=rank(n,n-k); endendL=inv(S); iT=;for i=1:p for j=0:sigmas(i) iT=iT; L(i+sum(sigmas(1:i),:)*Aj; end,endif knum=6 0 2 8 10;den=2 0 6 4 8; G=tf(num,den); Gs=sscanform(G,obsv)例4-8A=15,6,-12,9;4,14,8,-4;2,4,10,-2;9,6,-12,15; B=3 3;2

22、 2;-2 -2;3 9;T=luenberger(A,B) A1=inv(T)*A*T, B1=inv(T)*B424262810( )648sssG ssss156129334148422( )( )( )241022296121539x tx tu t4.1 4.1 线性系统性质分析-系统状态方程标准型的MATLAB求解(12)系统的范数测度及求解 连续系统的范数定义 H2范数： H范数： H范数是系统频域响应幅值的峰值2211( ) ()2pjijiG sG jdj ( )sup ()G sG j4.1 4.1 线性系统性质分析-系统状态方程标准型的MATLAB求解(13)离散系统的范

23、数定义范数的MATLAB求解 norm(G) 或 norm(G,2) norm(G,inf)221( ) ()pjiiG zG ed( )sup ()jG zG e4.1 4.1 线性系统性质分析-系统状态方程标准型的MATLAB求解(14)例4-9 已知离散系统模型A=-2.2 -0.7 1.5 -1;0.2 -6.3 6 -1.5;0.6 -0.9 -2 -0.5;1.4 -0.1 -1 -3.5; B=6 9;4 6;4 4;8 4; C=1 2 3 4; G=ss(A,B,C,0 0); norm(G,2), norm(G,inf)1690.26.361.546(1)

24、 ()()0.60.920.5484()1234 ()x kTx kTu kTy kTx kT 4.2 4.2 线性系统时域响应解析解法给线性系统一个激励信号，输出是什么？有两大类方法解析解方法 求解微分方程、差分方程解析解数值解方法主要内容基于状态方程的解析解方法基于传递函数部分方式展开的解析解方法二阶系统的解析解方法4.2 4.2 线性系统时域响应解析解法-基于状态方程的解析解方法(1)状态方程模型解析解求解难点( )( )( )( )( )( )x tAx tBu ty tCx tDu t00()()0( )( )( )tA t tA ttx tex teBud0(

25、)( )tA tteBud4.2 4.2 线性系统时域响应解析解法-基于状态方程的解析解方法(2)状态增广方法消除矩阵B，变成自治系统单位阶跃信号 ，若假设有另外一个状态变量 ，则其导数为增广状态方程自治系统可以直接求解析解 ( )1( )u tt1( )( )nxtu t1( )0nxt11( )( )( )( )00nnttABxtxtxx( )(0)Att = exx4.2 4.2 线性系统时域响应解析解法-基于状态方程的解析解方法(3)一般输入信号的系统增广一般输入信号模型引入增广状态变量11224340( )( )( )cos()sin()md tiiiu tu tu tctedd

26、tdd t111424(1)3131cos()sin()( ),( )d tnd tnmnn mxed txed txu txut4.2 4.2 线性系统时域响应解析解法-基于状态方程的解析解方法(4)增广状态方程模型其中，解析解( )(0)Att = exx231441000001000000,000Ad Bd BBddddA 12343( )( )( )( )( )( )( )nnnnn mtxtxttxtxtxt xx( )(0)Att =exx01(0)10(0)!mccc mxx4.2 4.2 线性系统时域响应解析解法-基于状态方程的解析解方法(5)MATLAB实现函数 functi

27、on Ga,Xa=ss_augment(G,cc,dd,X) G=ss(G); Aa=G.a; 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)

28、; Ca= Ca dd(2)*D; Xa=Xa; 1; Ba=B;0; end end4.2 4.2 线性系统时域响应解析解法-基于状态方程的解析解方法(6) if(length(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)+

29、i)=cc(i)*ii; end, end Ga=ss(Aa,zeros(size(Ca),Ca,D);调用格式信号描述00,ss_augment( , , ,)aGG xc d x010124,mc ccd d ddc =d4.2 4.2 线性系统时域响应解析解法-基于状态方程的解析解方法(7)例4-10 连续系统模型初值输入信号求解析解191616191211617190( )( )( )201716201201616192( )2,1,0,0 ( )ttu ty tt xxxT(0)0,1,1,2x3( )22sin(2 )tu tet4.2 4.2 线性系统时域响应解析解法-基于状态方

30、程的解析解方法(8)系统增广cc=2; dd=-3,0,2,2; x0=0;1;1;2; A=-19,-16,-16,-19; 21,16,17,19;20,17,16,20;-20,-16,-16,-19; B=1;0;1;2; C=2 1 0 0; D=0; G=ss(A,B,C,D); Ga, xx0=ss_augment(G,cc,dd,x0); Ga.a, xx0增广模型 -19-16-16-1902102116171900012017162002011( ) =( ),(0) =-20-16-16-1904220000-3-20100002-30000000002tt xxx4.2

31、 4.2 线性系统时域响应解析解法-基于状态方程的解析解方法(9)解析解求解syms t; y=Ga.c*expm(Ga.a*t)*xx0; latex(y);解析解的数学形式323312711913577( )54574cos(2 )sin(2 )4884tttttty tteeet eetet4.2 4.2 线性系统时域响应解析解法-基于部分分式展开方法求解(1)连续系统的解析解法 输入信号的Laplace变换U(s) 输出信号的Laplace变换Y(s)=G(s) U(s) 无重根时部分分式展开 由Laplace反变换求解析解112112121( )mmmmnnnnnbsb sb sbG

32、 ssa sa sasa1212( )mmrrrY sspspsp12-112( ) ( )mp tp tp tmy tY srerer eL4.2 4.2 线性系统时域响应解析解法-基于部分分式展开方法求解(2)有重根时相应项的解析解为部分分式的MATLAB求解：r,p,K=residue(num,den)112()()jjj mmjjjrrrspspsp11111111!(1)!111!(1)!jjjjmp tp tp tjmjjp tmjjjmrtr erteemrrtrtem-+-+-+-+-轾犏=+犏-臌LL4.2 4.2 线性系统时域响应解析解法-基于部分分式展开方法求解(3)例4

33、-11 系统的传递函数为 ，输入信号为阶跃信号R(s)=1/s，则输出信号的Laplace变换为MATLAB求解：num=1 7 3 4; den=1 7 17 17 6; R,P,K=residue(num, den, 0); R,P解析解：32432734( )717176sssG sssss+=+325432734( )717176sssY ssssss+=+32( )2.583395.753.50.667tttty teeete-=-+-+4.2 4.2 线性系统时域响应解析解法-基于部分分式展开方法求解(4)例4-12 带有复数极点的系统阶跃响应解析解num=1 3; den=1 2

34、 11 18 18; r,p,k=residue(num,den,0);r,p解析解4323( )2111818sG sssss+=+33( 1)( 1)( )(0.0020.0255 )(0.0020.0255 )( 0.08530.0088 )( 0.08530.0088 )0.1667jtjtj tj ty tj ej ej ej e- +- -=+-+ -+ -+4.2 4.2 线性系统时域响应解析解法-基于部分分式展开方法求解(5)解析解的进一步化简 基于Euler公式的化简 其中， 新MATLAB函数： r,p,K=pfrac(num,den) 若p(i)为实数，则(r(i),p(

35、i)同residue() 若p(i)为复数( )，则r(i),r(i+1)对返回A和()()()()sin()jtjttabj eabj eAetswswswF+-+-=+222,arctan(/ )Aaba bF= -+=-( )p ijsw=4.2 4.2 线性系统时域响应解析解法-基于部分分式展开方法求解(6)MATLAB函数function R,P,K=pfrac(num,den)R,P,K=residue(num,den);for i=1:length(R), if imag(P(i)eps a=real(R(i); b=imag(R(i); R(i)=-2*sqrt(a2+b2);

36、 R(i+1)=-atan2(a,b); elseif abs(imag(P(i)num=1 3; den=1 2 11 18 18; r,p,k=pfrac(num,den,0);r,p解析解4323( )2111818sG sssss+=+( )0.0511sin(30.0768)0.1715sin(1.4677)0.1677ty ttet-= -+4.2 4.2 线性系统时域响应解析解法-基于部分分式展开方法求解(7)基于Laplace变换的求解参见附录A步骤： 定义符号变量 描述原函数表达式 调用ilaplace()函数求解 结果化简，如simple()函数求解举例4.2 4.2 线性

37、系统时域响应解析解法-基于部分分式展开方法求解(8)例1 求系统 的阶跃响应解析解 MATLAB求解syms s; G=(s3+7*s2+3*s+4)/(s4+7*s3+17*s2+17*s+6); y=ilaplace(G/s) latex(y)解析解32432734( )717176sssG sssss+=+3223123( )7/293124ttty tetee-骣=+ -+-桫4.2 4.2 线性系统时域响应解析解法-基于部分分式展开方法求解(9)例2 求系统 的阶跃响应解析解 MATLAB求解：syms s; G=(s+3)/(s4+2*s3+11*s2+18*s+18); ilap

38、lace(G/s), latex(ans)解析解4323( )2111818sG sssss+=+1131293( )cos(3 )sin(3 )cossin2552556170170tty tttetet-=-+-4.2 4.2 线性系统时域响应解析解法-基于部分分式展开方法求解(10)离散系统的解析解法无重根时 部分分式展开 解析解-111nqqzpp p-轾骣犏= -犏-桫臌L1211112( )mmrrrY zzpzpzp-=+-L-1121122111( )( )nnnmmmrrry nY zpppppp骣骣骣鼢珑鼢= -珑鼢珑鼢珑桫桫桫LL对z-1进行展开4.2 4.2 线性系统时

39、域响应解析解法-基于部分分式展开方法求解(11)考虑采样周期例4-14 系统的离散传递函数为 ，输入为阶跃信号 D=conv(1 -1/2,conv(1 -1/4, conv(1 1/5, 1 -1); %分母 N=0 0 conv(1 -1/3, 1 0); %分子，因为要逆序排列，前面补足零 N=N(end:-1:1); D=D(end:-1:1); %因为需要对z-1进行展开，所以逆序排列 R,P,K=residue(N,D); R,P,-R./P-1121122111()( )kTkTkTmmmrrry kTY zpppppp骣骣骣鼢珑鼢= -珑鼢珑鼢珑桫桫桫LL(1/3)( )(1/

40、2)(1/4)(1/5)zG zzzz-=-+4.2 4.2 线性系统时域响应解析解法-基于部分分式展开方法求解(12)输出信号解析解Z变换求解步骤定义符号变量调用iztrans()函数求解化简11117.05473.95063.80951.4815( )5421Y zzzzz-=+-111( )1.41090.98771.90481.4815542nnny n骣骣骣鼢珑=-+鼢珑鼢珑桫桫桫-4.2 4.2 线性系统时域响应解析解法-基于部分分式展开方法求解(13)利用符号运算工具箱求解syms z G=(z-1/3)/(z-1/2)*(z-1/4)*(z+1/5) iztrans(G*z/(

41、z-1)求解结果408080040( )(1/2)(1/4)( 1/5)218156727nnny n = -+-+4.2 4.2 线性系统时域响应解析解法-基于部分分式展开方法求解(14)有重根时，部分分式表达式的Z反变换例4-15 系统的离散传函为 ，求阶跃响应的解析解部分分式展开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=residue(N(end:-1:1), D(end:-1:1); R P -11( 1)(1)(2)(1)()(1)!()mmnmqqnnnmzpmp-+

42、轾-犏=+-犏-臌LL352( )(1/2) (1/3)zG zzz-=-4.2 4.2 线性系统时域响应解析解法-基于部分分式展开方法求解(15)得到部分分式展开后Z变换表达式解析解11121313242409619236( )32(2)(2)1Y zzzzzz-=+-232324 1240 1961( )(1)3 322( 2)2192/2 1(1)(2)36( 2)211108( 126072)3632nnnnnny nnnnnn骣骣骣-鼢珑=+鼢珑鼢珑桫桫桫-骣+桫-骣骣鼢珑= -+-+鼢珑鼢珑桫桫4.2 4.2 线性系统时域响应解析解法-基于部分分式展开方法求解(16)符号运算求解s

43、yms z; G=(5*z-2)/(z-1/2)3/(z-1/3);iztrans(G*z/(z-1)解析解更直观，不建议用前者求解，而直接采用Z反变换的符号运算方法求解21111108726036123222nnnnnn骣骣骣骣鼢鼢珑珑-+-+-鼢鼢珑珑鼢鼢珑珑桫桫桫桫4.2 4.2 线性系统时域响应解析解法-基于部分分式展开方法求解(17)时间延迟系统的解析解法连续系统模型 : 求解G(s)的解析解，用t-L代替t即可离散系统传递函数 : 求解H(z)的解析解，用n-k替代n即可例4-16无延迟解析解：有延迟解析解： ( )LsG s e-( )kH z z-55352( )(1/2) (

44、1/3)zG z zzzz-=-211( )108( 126072)3632nny nnn骣骣鼢珑= -+-+鼢珑鼢珑桫桫55255211( )108 12(5)60(5)72)36 1(5)3211108( 126072)36 1(5)32nnnny nnnnnnn-骣骣鼢珑= -+-+-鼢珑鼢珑桫桫骣骣鼢珑= -+-+-鼢珑鼢珑桫桫4.2 4.2 线性系统时域响应解析解法-二阶系统的阶跃响应及指标(1)二阶系统开环模型单位负反馈的闭环控制系统模型阶跃响应的解析解其中2nn( )( +2)oG ss swzw=2n22nn( )+2oG ssswzww=+ndnd()()2ndndnd( )

45、12tteey tzwwzwwwwzwwzww-+-骣=+-+-桫2dn1wz w=-4.2 4.2 线性系统时域响应解析解法-二阶系统的阶跃响应及指标(2)根据 的不同取值，分为如下几种情况：u无阻尼振荡 ，u欠阻尼振荡 ，u临界阻尼振荡u过阻尼振荡z0z=n( )1cos()y ttw=-01z22(1)(1)n222( )12111tteey tzzzzwzzzzz-+-骣=-+-桫4.2 4.2 线性系统时域响应解析解法-二阶系统的阶跃响应及指标(3)二阶系统阶跃响应曲线wn=1; yy=;t=0:0.1:12;zet=0:0.1:1,2,3,5;for z=zet if z=0, y

46、=1-cos(wn*t); elseif (z0 & z1, dd=sqrt(z2-1); lam1=-z-dd; lam2=-z+dd; y=1-0.5*wn*(exp(lam1*t)/lam1-exp(lam2*t)/lam2)/dd; end yy=yy;y;endplot(t,yy)4.2 4.2 线性系统时域响应解析解法-二阶系统的阶跃响应及指标(4)三维曲面绘制：i=find(zetG=tf(10 20,10 23 26 23 10,ioDelay,1); step(G,30);利用MATLAB 提供的功能，可以从曲线上得到更多的信息，如超调量等4321020( )1023

47、262310ssG sessss-+=+4.3 4.3 线性系统的数字仿真分析-线性系统的阶跃响应与脉冲响应(3)MATLAB求解解析解num=G.num1; den=G.den1; r,p=pfrac(num,den,0); r,p得到解析解数值解精度比较：y,t=step(tf(num,den); y0=r(1)*exp(p(1)*t)+r(2)*t.exp(p(2)*t)+r(5)*exp(p(5)*t)+ r(3)*exp(real(p(3)*t).*sin(imag(p(3)*t+r(4); norm(y-y0) 0.15( )1.76470.58821.248sin(0.98870

48、.1897)2ttty teteet-= -+4.3 4.3 线性系统的数字仿真分析-线性系统的阶跃响应与脉冲响应(4)例4-18 离散化采样周期 T=0.01，0.1，0.5，1.2 G=tf(1,1 0.2 1,ioDelay,1); G1=c2d(G,0.01,zoh);G2=c2d(G,0.1); G3=c2d(G,0.5);G4=c2d(G,1.2); step(G,-,G2,-,G3,:,G4,-.,10)对得出的曲线进行比较，发现若采样周期选择过大，离散化后可能丢失原来系统的信息21( )0.21sG sess-=+4.3 4.3 线性系统的数字仿真分析-线性系统的阶跃响应与脉冲

49、响应(5)例4-19 求多变量系统阶跃响应MATLAB实现g11=tf(0.1134,1.78 4.48 1,ioDelay,0.72);g12=tf(0.924,2.07 1);g21=tf(0.3378,0.361 1.09 1,ioDelay,0.3);g22=tf(-0.318,2.93 1,ioDelay,1.29);G=g11, g12;g21, g22; step(G)0.7220.31.2920.11340.9241.784.4812.071( )0.33780.3180.3611.0912.931sssesssG seesss-轾犏犏+犏=犏-犏犏+臌4.3 4.3 线性系统

50、的数字仿真分析-线性系统的阶跃响应与脉冲响应(6)阶跃响应曲线是在两路输入单独作用下得出的存在系统耦合4.3 4.3 线性系统的数字仿真分析-线性系统的阶跃响应与脉冲响应(7)系统解耦引入前置补偿矩阵 对延迟环节引入Pad近似 n1,d1=paderm(0.72,0,2); g11.ioDelay=0; g11=tf(n1,d1)*g11; n1,d1=paderm(0.3,0,2); g21.ioDelay=0; g21=tf(n1,d1)*g21; n1,d1=paderm(1.29,0,2); g22.ioDelay=0; g22=tf(n1,d1)*g22; G=g11,g12;g21

51、,g22;补偿后系统的模型：Kp=0.1134,0.924;0.3378,-0.318; step(G*Kp);p0.11340.9240.33780.318K轾犏=犏-臌4.3 4.3 线性系统的数字仿真分析-线性系统的阶跃响应与脉冲响应(8)解耦效果较好！若要使得多变量系统能直接设计，在设计前必须解耦系统的脉冲响应曲线impluse()函数与step()函数调用结构一致G=tf(10 20,10 23 26 23 10,ioDelay,1); impulse(G,30);4.3 4.3 线性系统的数字仿真分析-任意输入下的系统响应(1)可以利用step()和impulse()函数求解输出信

52、号计算：Y(s)=G(s)R(s)如R(s)已知，则可以直接求解例4-20 斜坡响应斜坡信号的Laplace变换为1/s2，故系统响应可由两种方法求出：1. G(s)/s系统的阶跃响应2. G(s)/s2系统的脉冲响应4321020( )1023262310ssG sessss-+=+4.3 4.3 线性系统的数字仿真分析-任意输入下的系统响应(2)MATLAB求解：G=tf(10 20,10 23 26 23 10,ioDelay,1); s=tf(s);step(G/s); %或者 impulse(G/s2)其他输入的响应可以由lsim()函数求取：lsim(G,u,t)例4-21 多变量

53、系统输入信号0.7220.31.2920.11340.9241.784.4812.071( )0.33780.3180.3611.0912.931sssesssG seesss-轾犏犏+犏=犏-犏犏+臌12( )1sin(31),( )sin( )cos(2)tu tetu ttt-=-+=+4.3 4.3 线性系统的数字仿真分析-任意输入下的系统响应(3)MATLAB求解g11=tf(0.1134,1.78 4.48 1, ioDelay, 0.72);g12=tf(0.924,2.07 1);g21=tf(0.3378,0.361 1.09 1, ioDelay, 0.3);g22=tf(

54、-0.318,2.93 1, ioDelay, 1.29);G=g11,g12;g21,g22; t=0:0.1:15;u=1-exp(-t).*sin(3*t+1),sin(t).*cos(t+2);lsim(G,u,t);小节：多变量系统的时域响应可以这样求解比较容易理解曲线含义两个信号共同作用下系统的时域响应，与多变量系统阶跃响应的概念不同4.3 4.3 线性系统的数字仿真分析-降阶模型的时域分析及比较(1)前面介绍了降阶方法，但未比较效果例4-22MATLAB求解G=tf(1,7,11,5,1,7,21,37,30); Gr=pademod(G,1,2); Gr1=opt_app(G,

55、1,2,0) step(G,Gr,Gr1)通过阶跃响应曲线判断降阶模型的效果，可发现Pad近似算法的结果不是很理想，而次最优降阶算法效果较好。324327115( )7213730sssG sssss+=+4.3 4.3 线性系统的数字仿真分析-降阶模型的时域分析及比较(2)例4-23 非最小相位系统MATLAB求解G=tf(10 -60 110 60,1 17 82 130 100); Gr=opt_app(G,1,2,1); Gr1=opt_app(G,1,2,0); %获得带有时间延迟和不带时间延迟的最优降阶模型 step(G,Gr,Gr1)通过阶跃响应曲线判断降阶模型的效果324321

56、06011060( )1782130100sssG sssss-+=+4.4 根轨迹分析(1)单变量开环传递函数G(S)控制器增益K单位负反馈闭环系统特征方程1+kG(s)=0对k的不同取值，则可能绘制出每个特征根变化的曲线，这样的曲线称为系统的根轨迹。根轨迹用开环信息研究闭环特性4.4 根轨迹分析(2)MATLAB求解rlocus(G) rlocus(G,K) R,K=rlocus(G) rlocus(G1,-,G2,-.b,G3,:r)该函数可以用于单变量不含有时间延迟的连续、离散系统的根轨迹绘制，也可以用于带有时间延迟的单变量离散系统的根轨迹绘制4.4 根轨迹分析(3)例4-24 开环系

57、统MATLAB求解num=1 4 8;den=1,18,120.3,357.5,478.5,306; G=tf(num,den); rlocus(G)两个问题：1. 如何求解临界增益？2. 闭环系统稳定性如何变换2543248( )18120.3357.5478.5306ssG ssssss+=+4.4 根轨迹分析(4)例4-25根轨迹求解：s=tf(s); G=10/(s*(s+3)*(s2+3*s+4); rlocus(G), grid点击根轨迹图，求出阻尼比在=0.707处的增益求出临界增益下的阶跃响应曲线：k=0.524;step(feedback(G*k,1)210( )(3)(34

58、)G ss sss=+浙江大学电气工程学院系统系 包哲静4.4 根轨迹分析(5)4.4 根轨迹分析(6)例4-26 绘制离散系统(T=0.1秒)根轨迹z=tf(z,Ts,0.1);G=-0.95*(z+0.51)*(z+0.68)*.(z+1.3)* (z2-0.84*z+0.196)/.(z+0.66)*(z+0.96)*.(z2-0.52*z+0.1117)*.(z2+1.36*z+0.7328);rlocus(G),grid2220.95(0.51)(0.68)(1.3)(0.840.196)( )(0.66)(0.96)(0.520.1117)(1.360.7328)zzzzzG zz

59、zzzzz-+-+=+-+4.4 根轨迹分析(7)例4-27 离散系统模型T=0.1秒z=tf(z,Ts,0.1);G=0.52*(z-0.49)*(z2+.1.28*z+0.4385)/.(z-0.78)*(z+0.29)*.(z2+0.7*z+0.1586);rlocus(G)根据根轨迹曲线求取系统的临界增益220.52(0.49)(1.280.4385)( )(0.78)(0.29)(0.70.1586)zzzG zzzzz-+=-+4.4 根轨迹分析(8)系统假设有6步的纯滞后，则z=tf(z,Ts,0.1); G=0.52*(z-0.49)*(z2+1.28*z+0.4385)/(z

60、-0.78)*(z+0.29)*(z2+0.7*z+0.1586); G.ioDelay=6; rlocus(G)根据根轨迹曲线求取系统的临界增益，发现延迟后系统的临界增益减小。4.4 根轨迹分析(9)例4-28 延迟状态方程考虑无延迟系统根轨迹：A=-0.99,1.16,1.76,-0.16;-2.03, .-2.3,2.9,-2.45;-0.48,-3.96, .-2.05,-0.91;-0.43,1.23,2.26,-1.2; B=-1.3;-0.73;-0.57;0.62; C=-1.34,0.13,-1.11,0; G=ss(A,B,C,0); rlocus(G)0.991.161.732.