信号与系统教程.docx

9.1线性系统的状态方程9.1.1 状态变量与状态方程为了说明状态变量与和状态方程的概念，首先来研究图9-1所示的二阶电网络。图中us为电压源。由于电容电流和电感电压分别为ict=Cducdtult=Ldildt则由KCLhe KVL可列出如下方程ducdt=il-ucR2Ldildt=us-RliL-Uc整理可得ducdt=-1CR2Uc+1CiLdildt=-1LUc-R1LiL+1LUs若指定个电感电压uL为输出，则由方程uL=-uc-RliL+us式（9-3）形式的一阶微分方程称为状态方程，其中uc(t), iL(t)称为状态变量；式（9-4）是以输入信号和状态变量表示的代数方程，它称为输出方程。方程式（9-3）和（9-4）表明，如果电路在t=t0时刻的状态uc(t)和 iL(t)为已知，那么根据tt0时给定的输入us(t)就可以唯一的确定方程组式（9-3）的解uc(t)和 iL(t)。再由所得的状态变量和tt0时的输入，就可以确定tt0时的输出。 一般而言，连续动态系统在某一时刻t0的状态，是描述该系统所必须的最少的一组数x1t0,x2t0,xn(t0),根据这组数和tt0时给定的输入就可以唯一的确定在tt0的任意时刻的状态。状态变量x1t0,x2t0,xn(t0)，是描述状态随时间变化的一组变量，它们在t0时的值就组成了系统在该时刻的状态。状态变量方程简称状态方程，它是用状态变量和激励（有时为零）表示的一组独立的一阶微分方程；而输出方程是用状态变量和激励（优势还可能有激励的某些函数）表示的袋鼠方程。通常将系统的状态方程和输出方程总称为动态方程。在电路系统中，一般取独立电容上电压和电感中电流为状态变量。下面以图9-2所示的二姐电路为例，进一步说明直观列写状态方程的方法。选择uc(t)和 iL(t)为状态变量，对 包含电感的回路，有KVL，有diLdt-uct+ust=0对于节点a，由KCL，有Cducdt=-1RCuc-1CiL+1Cis整理可得Cducdt=-1RCuc-1CiL+1Cisdildt=1Luc(t)-1Lus(t)写成矩阵形式为 ducdtdildt=-1RC-1C1L0uc(t)iL(t)+1C00-1Lisus状态变量xt=x1(t)x2(t)=uc(t)iL(t), A=-1RC-1C1L0, B=1C00-1Lft=is usT则上式可写为X=Axt+Bf(t)式中，X表示状态变量的一阶导数，x(t)称为状态变量，A称为状态变量的系数矩阵，对线性时不变系统，A和B为常数矩阵。式（9-5）为有外加输入系统的状态方程的标准形式。 由上可知，状态方程式一组一阶微分方程，只要知道其实状态，就可以求取uc(t)和 iL(t)，随之该电路的其他量也可以确定。因此uc(t)和 iL(t)是该电路中最少的一组状态变量。对于一般的电路，直观列写状态方程的步骤如下：第一步：选择独立的电容上电压和电感中电流为状态变量。第二步：对于电容相连的节点列写KCL方程，对包含电感的回路列写KVL方程。第三部：消去非状态变量，整理成标准形式的状态方程x=Axt+Bf(t) 式中，Bf(t)是与外加信号有关的项，B为常数矩阵。例9-1对于图9-3（a）所示电路，试写出其状态方程，并以uR为输出写出输出方程。解 首先选择电容上电压和电感上电流为状态变量，即xt=uc1 uc2 iLT画出与电路相对于的图，如图9-3（b）所示。对于C1和C2支路相连的节点分别列出KCL方程，即C1duc1dt+uc1-usR+iL=0C2duc2dt+iL+iL=0对电感所在的回路列写KVL方程，即LdiLdt-uc2-uc1+us=0整理可得状态方程duc1dt=1C1-uc1R-iL+usRduc2dt=1C2-iL-isdiLdt=1Luc1+uc2-us写出矩阵形式为duc1dtduc2dtdiLdt=-1C1R0-1C00-1C21L1L0uc1uc2iL+-1C1R00-1C2-1L0usis即有标准形式x=Axt+Bft输出方程urt=ust-uc1t对于一般系统，如果已知其模拟框图，也可以写出它们的状态方程。例如图9-4所示的三种简单例子，只要把每个积分器的输出变量设为状态变量，即可容易地写出状态方程和输出方程。对于图9-4（a）所示的反馈积分放大环节，设积分器输出为状态变量x，则有状态方程x=ax+bft 对于图9-4（b）所示的两级反馈环节级联情况，有状态方程x1=-a1x1+x2x2=-a2x2+ftx1x2=-a110-a2x1x2+01ft输出方程为yt=x1对于图9-4（c）所示的两级反馈环节并联情况，有状态方程x1=-3x1+2f(t)x2=-2x2+2ft即x1x2=-310-2x1x2+22ft输出方程yt=x1-x29.1.2 动态方程的一般形式 利用状态，状态变量，状态方程和输出方程的概念，我们研究系统一般情况。设线性系统的n个状态变量为x1t,x2t,xn(t),系统的p个输入为f1t,f2t,fn(t),则系统的状态方程可一般的写为dx1dt=a11x1+a1nxn+b11f1+b12f2+b1pfpdx2dt=a21x1+a2nxn+b21f1+b22f2+b2pfpdxndt=an1x1+annxn+bn1f1+bn2f2+bnpfp式中，各系数aij和bij是由系统参数决定的系数，对于线性时不变系统，它们都是常数；对于线性时变系统，它们是时间的函数。 把上述方程组些微矩阵形式为x1x2.xn= a11 a12 a1n a21 a22 a2n .an1 an2 ann x1x2.xn+ b11 b12 b1n b21 b22 b2n .bn1 bn2 bnn f1f2.fn可简记为矩阵方程x=Axt+Bft式中xt=x1 x2 xnTx(t)=x1 x2 xnTf(t)=f1 f2 fnT分别是状态变量的倒数，状态变量和激励组成的列矩阵。A= a11 a12 a1n a21 a22 a2n .an1 an2 ann B= b11 b12 b1n b21 b22 b2n .bn1 bn2 bnn 分别是系数矩阵，其中A是n*n阶方阵。 如果网络有q个输出y1t,y2t,yq(t),那么，它们中的每一个都是由状态变量和激励表示的代数方程。其矩阵形式为y1y2.yn= c11 c12 c1n c21 c22 c2n .cq1 cq2 cqn x1x2.xn+ d11 d12 d1p d21 d22 d2p .dq1 dq2 dqp f1f2.fp上式可简记为yt=Cxt+Dft式中，系数矩阵C= c11 c12 c1n c21 c22 c2n .cq1 cq2 cqn D= d11 d12 d1p d21 d22 d2p .dq1 dq2 dqp 上述系数矩阵，对于线性时不变系统，它们都是常数矩阵。式（9-7）和式（9-8）分别是状态方程和输出方程的一般形式。由于x(t)可看做是n维空间的向量称为状态向量，故状态变量分析也常称为状态空间分析。应用状态空方程和输出方程的概念，可以研究许多复杂的工程问题。如倒立振子的运动，机器人双脚的行走，飞机和火箭的升空等。图9-5是火箭升空控制的示意图。其中x1,x2,x3和x4为状态变量。 9.2 状态方程的解9.2.1 状态方程的时域求解 这里研究状态方程在时间域的求解方法。有状态方程的一半形式x=Axt+Bf(t)可改写为x-Axt=Bf(t)它与一阶电路的微分方程yt-ayt=bft相似。将a换为A，则状态方程的解可写为x(t)=eAtx0+0teAt-Bfd或者表示为xt=eAtx0+eAt*Bft式中，x(0)为其实状态，以后设状态连续，不分0_和0+。式中，eAt是一矩阵指数函数，通常称为状态转移矩阵函数。观察上式可知，为求x(t)，关键是设法求出eAt。 在时间域中求解eAt并非易事。为了突出应用，我们将在后面介绍较简便的s域求解法。一般情况下，设动态系统为xt=Axt+Bft可以证明yt=CeAt*Bft+Dft系统的冲击响应矩阵为ht=CeAtB+Dt式中，t为对角阵。 对于给定的状态方程，一种比较方便的方法是借助计算机用逐次逼近法或者龙格-库塔法等求解。下面署名逐次逼近法。设状态方程xt=axt+bft将x(t)展开泰勒级数并取一次项近似，得 xt=x0+dx(0)dt.t. x2t=xt+dx(t)dt.t. x(k+1)t=xkt+dx(kt)dt.t又因为在kt时状态方程变为 xkt=axkt+bf(kt)将式（9-14）带入式（9-13），得xk+1t=axkt+axktt+Bf(kt) t.把上式推广于矩阵形式的状态方程，得xk+1t=Axkt+Axktt+Bf(kt) t.上式为逐次逼近法的一半公式，当步长t走足够小时，求解结果非常接近真实解。9.2.2 状态方程的s域求解 拉普拉斯变换是求解常系数线性微分方程的有力工具，现用于求解状态方程。 由函数矩阵积分的定义，状态向量x(t)的拉氏变换为简记为Xs=Lxt同样地，输入和输出向量的拉式变化简记做Fs=LxtYs=Lyt根据拉氏变换的积分性质LLdxtdt=sXs-x0应用于状态方程xt=Axt+Bft则有变换sXs-x0=AXs+BFs即sXs-AXs=x0+BFs利用单位矩阵I，又可以写为sI-AXs=x0+BFs将上式两端前乘以sI-A-1，得Xs=sI-A-1x0+BFs=sI-A-1x0+sI-A-1BFs上式是状态向量x(t)的拉氏变换。由此可知，其第一项的反变换是状态向量的零输入解，第二项的反变换是状态向量的零状态解。 对照式（9-9，且考到x(0)是常数矩阵，应有eAtx0=-1sI-A-1x0=-1sI-A-1x0于是的状态转移矩阵eAt=-1sI-A-1=-1adjsI-AdetsI-A上式给我们提供了另一种求eAt的方法。为了方便，我们定义分解矩阵s=eAt=sI-A-1于是式（9-17）可写为Xs=sx0+sBFs对于输出方程yt=Cxt+Dft取拉式变换，得Ys=CXs+DFs将式（9-20）带入上式，可得Ys=Csx0+CsB+DFs有上式可知，零状态响应的拉氏变换 Yzss=CsB+DFs=H(s)F(s)式中，系统函数矩阵Hs=CsB+D它是一个q*p阶矩阵（q是输出书，p是输入数）。由此可得Lht=H(s)例9-2 如有矩阵A=-100-2试用s域方法求状态转移矩阵eAt。解 由已知矩阵A可得sI-A=s+100s+2(sI-A)-1=1s+1(s+2)s+200s+2=1s+1001s+2因此eAt=-1sI-A-1=e-t00e-2t例9-3 设有系统状态方程和输出方程分别为x1(t)x2t=120-1x1(t)x2t+0110f1(t)f2ty1(t)y2t=110-1x1(t)x2t+1010f1(t)f2t其起始状态x1(0)x20=1-1输入向量f1(t)ft=(t)t试求状态变量x(t)和响应y(t).解 由状态方程和输出方程知系数矩阵A=120-1，B=0110 C=110-1，D=1010所以sI-A=s-1-20s+1detsI-A=s-1s+1adjsI-A=s+120s-1故分解矩阵s=sI-A-1=adjsI-AdetsI-A=1s-12s+1(s-1)01s+1又输入向量的变换Fs=(t)t=1s1故得状态向量的象函数Xs=sx0+sBFs=1s+1-1s+1+s2+s+2ss-1(s+1)1ss+1取拉式反变换，得xt=x1(t)x2t=e-t-e-t+2et+e-t-21-e-t而又可以得输出向量的像函数矩阵Ys=Csx0+CsB+DFs=01s+1+2s-11s+1取反变换，的输出 y1(t)y2(t)=0e-t+2ete-t (t0) 例9-4 设有二姐电路系统的状态方程和输出方程分别为ducdtdiLdt=-31-21uciL+0120usisucuL=10-10uciL+0010usis试求其系统函数矩阵H(s).解 由方程知A=-31-20故sI-A=s+3-12sadjsI-A=s1-2s+3detsI-A=s+3-12s=s2+3s+2=s+1s+2故分解矩阵s=sI-A-1=adjsI-AdetsI-A=ss+1(s+2)1s+1(s+2)-2s+1(s+2)s+3s+1(s+2)把s和系数矩阵C，D，B代入式（9-28），得系统函数矩阵Hs=CsB+D=10-10ss+1(s+2)1s+1(s+2)-2s+1(s+2)s+3s+1(s+2)0120+0010=2s+1(s+2)ss+1(s+2)1-2s+1(s+2)-ss+1(s+2)取上式的反变换可得冲击响矩阵。 以上讨论重点在求解的概念和分析过程。实际中已有简单的计算机分析方法，这里不在叙述。下面介绍状态轨迹的概念。选取两个状态变量x1和x2构成一个状态平面，从初始状态起到t期间，x1t1,x2(t1), x1t2,x2(t2),对应该平面的各点，连接这些点的轨迹称为状态轨迹。例如，设有状态方程ddtuciL=01C-1L-RLuciL当R=1,L=1,C=1F, uc=2 V, iL0=2A,并取t=0.01s，在t=03s中求解，其状态轨迹如图9-6所示。 9.4 离散系统的状态变量分析 线性离散时间系统可以用输入-输出发分析，也可以用状态变量发分析。用输入-输出法分析时，描述系统的是n阶差分方程；用状态变量分析时，则是一阶线性差分方程组。9.4.1 离散系统的状态方程 对于线性时不变离散系统，与连续系统相似，也可以表示为规范的状态方程形式。其状态方程和输出方程都是离散的状态变量和输入序列的线性组合。状态方程的一半形式为x1n+1=a11x1(n)+a1mxm(n)+b11f1(n)+b12f2(n)+b1pfp(n)x2n+1=a21x1(n)+a2mxm(n)+b21f1(n)+b22f2(n)+b2pfp(n)xmn+1=am1x1(n)+ammxm(n)+bm1f1(n)+bm2f2(n)+bmmfp(n)若有q个输出，则输出方程的一般形式为yn+1=c11x1(n)+c1mxm(n)+d11f1(n)+d12f2(n)+d1pfp(n)y2n+1=c21x1(n)+c2mxm(n)+d21f1(n)+d22f2(n)+d2pfp(n)ymn+1=cq1x1(n)+cqmxm(n)+dq1f1(n)+dq2f2(n)+dqpfp(n)写成矩阵形式，分别为xn+1=Axn+Bfnyn=Cxn+Dfn式中xn=x1n x2n xm(n)Tfn=f1n f2n fp(n)Tyn=y1n y2n yq(n)T各系数矩阵A= a11 a12 a1n a21 a22 a2n .an1 an2 ann ，B= b11 b12 b1n b21 b22 b2n .bn1 bn2 bnn C= c11 c12 c1n c21 c22 c2n .cq1 cq2 cqn ，D= d11 d12 d1p d21 d22 d2p .dq1 dq2 dqp 9.4.2 离散状态方程的时域求解 求解离散状态方程的简单方法之一是所谓的迭代法（递推法）。例如，对于状态方程x1nx2n=12 014 14x1nx2n+10fn设起始状态x1(0)x2(0)=00输入序列为fn=n令n=0时，有x1(1)x2(1)=12 014 1400+101=10令n=1时，有x1(2)x2(2)=12 014 1410+100=1214令n=2时，有x1(3)x2(3)=12 014 141214+100=14316依次类推，如此不断进行，就可以求的状态变量的解x11,x21，x3(1)，.; x11,x21，x3(1),.。但这种方法很难得到闭式解。 对于一般状态方程，可以证明xn=Anx0+k=0n-1An-1-kBfn输出方程的解yn=CAnx0+k=0n-1CAn-1-kBfk+Dfn值得注意的是。式（9-32）和式（9-33）中的零输入解的n0；零状态响应中k=0n-1（.）的n1；输入序列f(.)为n0。9.4.3 状态昂藏的z域求解法对于LTI离散系统的状态方程xn+1=Axn+Bfn设状态变量矩阵x(n)和输入序列矩阵f(n)的z变换分别为Xz=zx(n)Fz=zf(n)对式（9-34）取z变换，得zXz-zx0=AXz+BFz也可以写为zI-AXz=zx0+BFz上式两端前乘以zI-A-1，得Xz=zI-A-1zx0+zI-A-1BFz取上式的z反变换，可得x(n)其中第一项的反变换为零输入解，第二项的反变换为零状态解。比较可知，状态转移矩阵An=z-1zI-A-1z这又提供了一种求An的方法。 对于输出输出方程yn=Cxn+Dfn对其z变换，得Yz=CXz+DFz把X(z)代入，有Yz=CzI-A-1zx0+CzI-A-1BFz+DFz= CzI-A-1zx0+CzI-A-1B+DFz =CzI-A-1zx0+HzFz式中，第一项为零输入响应的象函数，第二项为零状态响应的象函数。其中Hz=CzI-A-1B+D为系统函数矩阵。例9-5 设矩阵A=5 -61 0试用z域方法求An。解 特征矩阵zI-A=z1 00 1-5 -61 0=z-5 -6-1 z它的逆矩阵(zI-A)-1=adj(zI-A)det(zI-A)=1z2-5z+6z -61 z-5=zz2-5z+6 -6z2-5z+61z2-5z+6 z-5z2-5z+6An=z-1sI-A-1z=z-1zz2-5z+6 -6z2-5z+61z2-5z+6 z-5z2-5z+6故有An=3n+1-2n+1 -6(3)n+6(2)n3n-2n -2(3)n+3(2)n9.5 MATLAB方法用于状态变量分析9.5.1 状态方程的解例 9-6 已知RLC并联电路的状态方程为diLdt=1Lucducdt=-1ciL-1RCuc设初始值iL0=1A, uc0=1V,在下列情况下求解iL和uc，并画出iL-uc平面的状态轨道。(a) R=0.4,L=0.1 H,C=0.1F;(b) R=5/6,L0.1 H,C=0.1F。解ex.m:R=input(电阻R=)；%以交互方式输入电阻R的值L=input(电感L=);%以交互方式输入电感L的值C=input(电容C=);%以交互方式输入电容C的值A=0 1/L;-1/C -1/(R*C);%状态矩阵ASyms t;%符号变量tF=expm(A*t);%计算exp(A*t) X0=1;1;%电流和电业的初始值i(0)和u(0)X=F*x0%求解电流i和电压uT=0:0.02:2I=subs(X(1,:),t, t)Subplot(1,3,1),plot(t,I);%绘制i(t)的曲线Title(i(t)的曲线)；Grid on;axis aquare;U=subs(X(2,:), t, t);Subplot(1,3,2),plot(t,U); 绘制u(t)的曲线Title(u(t)的曲线)；Grid on;axis square;Subplot(1,3,3),plot(I,U);%绘制i(t)-u(t)的状态轨道Title(i(t)-u(t)的状态轨道)；Xlabel(i(t);ylabel(u(t);Grid on;axis square;(a) 电阻R=0.4,电感L=0.1，电容C=0.1X=-exp(-20*t)+2*exp(-5*t)-exp(-5*t)+2*exp(-20*t)即iLt=-e-20t+2e-5tt,uct=(-e-5t+2e-20t)t(b) 电阻R=5/6电感L=0.1电容C=0.1X=exp(-6*t)*cos(8*t)+2*exp(-6*t)*sin(8*t)-2*exp(-6*t)*sin(8*t)+exp(-6*t)*cos(8*t)即iLt=e-6tcos8t+2sin8tt,uc(t)=e-6t(-2sin8t+cos8t)t例9-7 设有状态方程x1x2=120-1x1x2+0110f1(t)f2(t)y1(t)y2(t)=110-1x1x2+1000f1(t)f2(t)x1(t)x2(t)=1-1f1t=t, f2t=t试求x1t和x2(t)，y1t和y2(t)，打印个曲线。解ex.mA=1 2;0 -1;B=0 1;1 ;0;C=1 1;0 -1;D=1 0;1 0;X0=1 -1;T=0:dt:2;F(:,1)=oned(length(t),1);%f1(t)=u(t)F(:,2)=1,zeros(length(t)-1,1);%f2(t)=dirac(t)Sys=ss(A,B,C,D);%状态方程模型y t0 x=lsim(sys,f,t,x0);%求解状态方程Subplot(2,2,1),plot(t,x(:,2);%绘制x1(t)Xlable(t),ylabel(x1(t);Subplot(2,2,2),plot(t,x(:,2);%绘制x2(t)Xlable(t),ylabel(x2(t);Subplot(2,2,3),plot(t,x(:,2);%绘制y1(t)Xlable(t),ylabel(y1(t);Subplot(2,2,4),plot(t,x(:,2);%绘制y2(t)Xlable(t),ylabel(y1(t);运行结果如图9-12所示。例9-8设有状态方程（洛伦兹LORenz方程x=-ax+ayy=x-y-xzz=xy-z）当a=0, =8/3, =28时出现混沌。设x(0)=3,y(0)=3.5,z(0)=4.2.求x(t),y(t),z(t)波形