控制系统的设计与仿真MATLAB控制系统设计与仿真教学课件PPT课件

《控制系统的设计与仿真MATLAB控制系统设计与仿真教学课件PPT课件》由会员分享，可在线阅读，更多相关《控制系统的设计与仿真MATLAB控制系统设计与仿真教学课件PPT课件（238页珍藏版）》请在装配图网上搜索。

1、4.1 系统建模与仿真框图的创建 问题描述 众所周知,现实世界中存在着各种不同的控制系统。对于线性时不变（LTI）系统,一般可以分为连续和离散系统。MATLAB中为用户提供了丰富的针对各种系统的建模手段。图4.1显示了MATLAB中各种线性时不变（LTI）系统之间的转换关系。 第1页/共238页图4.1 连续与离散系统的关系示意图 第2页/共238页 图4.1中显示了MATLAB可以完成离散和连续系统的建模,并且同一系统可以表示成连续系统,也可以表示成离散系统,它们之间可以以状态方程形式进行转化。这一节将结合一个具体实例来演示MATLAB中各种模型创建和相互之间进行转化的方法,以及如何用Sim

2、ulink进行连续系统的仿真。首先给出实例的源程序MODLDEMO.M,然后根据不同的主题结合实例进行讲述。 第3页/共238页 例4.1 对于Mass Spring Dashpot机械系统： .( )( )( )m y tc yky tu t(4.1) 试建立该系统的连续和离散模型,并进行时域和频域仿真。 解：程序源代码如下： %MODLDEMO.M 演示各种建模与仿真（时域和频域）技巧 clearall,closeall %程序开始,清空工作空间,关闭所有窗口 deletemodldemo.out,diarymodldemo.out %打开二进制文件 第4页/共238页 disp(*MOD

3、LDEMO.OUT* DiaryFileforMODLDEMO.M),disp() m=1 %各种系统常数定义 k=1 %单位kg/s2 c= %单位kg/s 第5页/共238页 %创建系统的状态空间模型 disp(StateSpaceModels) km=k/m; A1=01;-km-c(1)/m,A2=01;-km-c(2)/m A3=01;-km-c(3)/m,A4=01;-km-c(4)/m B=01/m,C=10,D=0 sys1s=ss(A1,B,C,D);sys2s=ss(A2,B,C,D); sys3s=ss(A3,B,C,D);sys4s=ss(A4,B,C,D); 第6页/

4、共238页 %仿真系统的脉冲和阶跃响应（时域） t=0:.2:15; y1=impulse(sys1s,t);y2=impulse(sys2s,t); y3=impulse(sys3s,t);y4=impulse(sys4s,t); figure(1) subplot(221),plot(t,y1,r),title(CriticalDampingImpulse),grid xlabel(Time),ylabel(SystemResponse) subplot(222),plot(t,y2,r),title(OverDampingImpulse),grid xlabel(Time),ylabel

5、(SystemResponse) 第7页/共238页 subplot(223),plot(t,y3,r), title(UnderDampingImpulse),grid xlabel(Time),ylabel(SystemResponse) subplot(224),plot(t,y4,r), title(NoDampingImpulse),grid xlabel(Time),ylabel(SystemResponse) 第8页/共238页 % y1=step(sys1s,t);y2=step(sys2s,t);y3=step(sys3s,t); y4=step(sys4s,t); figu

6、re(2) subplot(221),plot(t,y1,r), title(CriticalDampingStep),grid xlabel(Time),ylabel(SystemResponse) subplot(222),plot(t,y2,r), title(OverDampingStep),grid第9页/共238页 xlabel(Time),ylabel(SystemResponse) subplot(223),plot(t,y3,r), title(UnderDampingStep),grid xlabel(Time),ylabel(SystemResponse) subplot

7、(224),plot(t,y4,r), title(NoDampingStep),grid xlabel(Time),ylabel(SystemResponse) disp(hitanykeytocontinue),pause第10页/共238页 %在Matlab中进行模型转化。 对于m=k=1,系统简化的传递函数为 %G(s)=1/s2+cs+1=z(s)/p(s) disp(TransferFunctionForm) sys1t=tf(sys1s),sys2t=tf(sys2s) sys3t=tf(sys3s),sys4t=tf(sys4s) disp(hitanykeytocontinu

8、e),pause第11页/共238页 %将传递函数形式转化成零极点形式 disp(ZeroPoleGainForm) sys1z=zpk(sys1t),sys2z=zpk(sys2t) sys3z=zpk(sys3t),sys4z=zpk(sys4t) disp(hitanykeytocontinue),pause %演示留数定理的使用 disp(ResidueForm) n1,d1=tfdata(sys1t);n2,d2=tfdata(sys2t); n3,d3=tfdata(sys3t);n4,d4=tfdata(sys4t); 第12页/共238页 %注意tfdata命令定义元胞数组 %

9、将元胞数组转化成向量形式 n1=n11,d1=d11,n2=n21,d2=d21, n3=n31,d3=d31,n4=n41,d4=d41, %使用留数定理 r1,pr1,kr=residue(n1,d1),r2,pr2,kr=residue(n2,d2) r3,pr3,kr=residue(n3,d3),r4,pr4,kr=residue(n4,d4) disp(hitanykeytocontinue),pause第13页/共238页 %进行频域仿真,这里使用nyquist函数（参考bode,freqs等函数） w=logspace(-2,2,100); re,im=nyquist(sys1

10、s,w);%SS形式 re1(:,1)=re(1,1,:);im1(:,1)=im(1,1,:);g1=re1+i*im1; mag1=20*log10(abs(g1);phase1=angle(g1)*180/pi; re,im=nyquist(sys2s,w);%SS形式 re2(:,1)=re(1,1,:);im2(:,1)=im(1,1,:);g2=re2+i*im2; mag2=20*log10(abs(g2);phase2=angle(g2)*180/pi; re,im=nyquist(sys3t,w);%TF形式 re3(:,1)=re(1,1,:);im3(:,1)=im(1,

11、1,:);g3=re3+i*im3; 第14页/共238页 mag3=20*log10(abs(g3);phase3=angle(g3)*180/pi; re,im=nyquist(sys4t,w);%TF形式 re4(:,1)=re(1,1,:);im4(:,1)=im(1,1,:);g4=re4+i*im4; mag4=20*log10(abs(g4);phase4=angle(g4)*180/pi; %幅值曲线 figure(3) subplot(221),semilogx(w,mag1,r),title(CriticalDamping), xlabel(frequency),ylabe

12、l(|G( jw)|indb) subplot(222),semilogx(w,mag2,r),title(OverDamping), xlabel(frequency),ylabel(|G( jw)|indb)第15页/共238页 subplot(223),semilogx(w,mag3,r),title(UnderDamping), xlabel(frequency),ylabel(|G( jw)|indb) subplot(224),semilogx(w,mag4,r),title(NoDamping), xlabel(frequency),ylabel(|G( jw)|indb) di

13、sp(hitanykeytocontinue),pause %频率曲线 figure(4) subplot(221),semilogx(w,phase1,r), title(CriticalDamping), xlabel(frequency),ylabel(angle)第16页/共238页 subplot(222),semilogx(w,phase2,r), title(OverDamping), xlabel(frequency),ylabel(angle) subplot(223),semilogx(w,phase3,r), title(UnderDamping), xlabel(fre

14、quency),ylabel(angle) subplot(224),semilogx(w,phase4,r), title(NoDamping), xlabel(frequency),ylabel(angle) disp(hitanykeytocontinue),pause第17页/共238页 %绘制系统Nichols图（对数坐标形式） figure(5) subplot(221),plot(phase1,mag1,r+), title(CriticalDamping), xlabel(angle),ylabel(|G( jw)|indb), gtext(w=0) subplot(222),

15、plot(phase2,mag2,r+), title(OverDamping), xlabel(angle),ylabel(|G( jw)|indb),gtext(w=0) subplot(223),plot(phase3,mag3,r+),title(UnderDamping), 第18页/共238页 xlabel(angle),ylabel(|G( jw)|indb),gtext(w=0) subplot(224),plot(phase4,mag4,r+),title(NoDamping), xlabel(angle),ylabel(|G( jw)|indb),gtext(w=0) di

16、sp(hitanykeytocontinue),pause第19页/共238页 %绘制Nyquist图(实部-虚部形式) figure(6) subplot(221),plot(re1,im1,r+), title(CriticalDamping), xlabel(Real),ylabel(Imag), gtext(w=0) subplot(222),plot(re2,im2,r+), title(OverDamping), xlabel(Real),ylabel(Imag), gtext(w=0)第20页/共238页 subplot(223),plot(re3,im3,r+),title(U

17、nderDamping), xlabel(Real),ylabel(Imag),gtext(w=0) subplot(224),plot(re4,im4,r+),title(NoDamping), xlabel(Real),ylabel(Imag),gtext(w=0) disp(hitanykeytocontinue),pause 第21页/共238页 %直接使用Nyquist函数绘制Nyquist图 figure(7) subplot(2,2,1),nyquist(sys1s,w),title(CriticalDamping), subplot(2,2,2),nyquist(sys2s,w

18、),title(OverDamping), subplot(2,2,3),nyquist(sys3s,w),title(UnderDamping), 第22页/共238页 subplot(2,2,4),nyquist(sys4s,w),title(NoDamping), disp(hitanykeytocontinue), disp(hitcntrl-ctogetoutofthefileforinteractiveanalysis) ,pause %使用Simulink框图（文件名为mdemosl.mdl,如图4.2所示）进行仿真 clearall%清空所有内存第23页/共238页 m=1,k

19、=1,c=1.2%定义常数 disp(DatafromgraphicalSimulinkmodel) %计算线性化模型 A,B,C,D=linmod(mdemosl),syss=ss(A,B,C,D); sysz1=zpk(syss),syssm=minreal(syss),sysz2=zpk(syssm) diaryoff%关闭二进制文件第24页/共238页图4.2 Mass Spring Dashpot系统仿真框图 第25页/共238页 物理模型 例4.1 中研究的对象是一个简单的物质交换机械系统,可以用微分方程表示成22( )( )( )( )ddmy tcy tky tu tdtdt(

20、4.2) 其中,y(t)是系统的瞬时交换的质量,k和c分别为比例常数。如果取状态变量x1=y和x2=dy/dt,则可以得到系统的状态方程形式：第26页/共238页 可以看出,矩阵A、B、C和D可以看成LTI系统的标准状态方程形式,于是式（4.3）可以写成112212010( )/1/( )10 xxdu txk mc mxmdtxy tx(4.3) dXAXBUdtYCXDU第27页/共238页 时域仿真 一般对控制系统进行时域仿真可以采用impulse、step和lsim函数,这些函数用来处理系统的状态空间描述。例如,对于输入信号u(t),系统的动态仿真可表示成 sys=ss(A,B,C,D

21、) Y,T,X=lsim(sys,U,t,xo) 第28页/共238页 当然,这些函数也可以处理系统的频域表达形式,这时,LTI系统可以通过下式表示 sys=tf(num,den) 作为一个典型的例子,例4.1使用了impulse和step函数来仿真系统的时域特性。图4.3和图4.4分别为例4.1执行的仿真结果。第29页/共238页图4.3 典型二阶系统的脉冲响应曲线第30页/共238页图4.4 典型二阶系统的阶跃响应曲线 第31页/共238页 模型转换 LTI系统的频域描述可以用传递函数表示为 Y(s)=G(s)U(s) (4.5) 系统的传递函数矩阵为1( )G sC sIABD(4.6)

22、 对于例4.1研究的简单SISO机械系统,其传递函数可以写成简单的标量形式 21( )mG sckssmm(4.7) 第32页/共238页 但是对于多输入多输出（MIMO）系统而言,其传递函数描述就有些复杂了,这时MATLAB中的模型转换函数可以发挥作用,它可以完成系统在状态方程形式与传递函数形式之间的互换,同时也可以将传递函数形式转换成零极点-增益形式。相关的函数包括 sys1=ss(A,B,C,D) sys2=tf(sys1) sys3=zpk(sys2)第33页/共238页 我们也可以采用ssdata、tfdata和zpkdata等命令将存储在与一个指定LTI对象相联的数据结构中的信息抽

23、取出来。例如 num,den=tfdata(sys2) 返回LTI对象sys2的分子和分母多项式系数,num与den为相应的元胞数组,其行数为输出的维数,列数等于输入的维数。其中第i行第j列元素表示从第j个输入到第i个输出的传递函数。第34页/共238页 另一种从数据结构中得到元胞数组的方法是使用MATLAB的celldisp命令。简单的显示数据的信息,可以使用如下的命令: fieldnames(sys2) num1=sys2.num,den1=sys2.den celldisp(num1),celldisp(den1) 也可以对零极点-增益形式完成显示的操作。例如 Z,P,K=zpkdata

24、(sys3) 第35页/共238页 将获取LTI系统sys3每一个IO通道的零极点和增益大小。元胞数组Z、P和矩阵K 的行与列分别与输出和输入的维数相同。其中第i行第j列元素表示从第j个输入到第i个输出传递函数的零极点和增益。 对于单输入单输出（SISO）系统,其传递函数与零极点-增益形式可以简化成普通的分数形式，即12123( )()()( )( )()()()num sszszG sKden sspspsp(4.8) 第36页/共238页 留数定理 我们也可以将系统写成几个分数相加的形式,例如对于SISO的机械系统,G(s)可以写成1121221,21( )( )( )( )42B srr

25、mG sk sckA sspspssmmccmkpm (4.9)(4.10) 第37页/共238页 r1,r2可以通过各种不同的方法计算得到。对于上述问题,MATLAB可用residue函数来完成这一运算，即 r,p,k=residue(B,A) 其中B和A 为包含多项式系数的行向量,而r和p 是包含留数和极点的列向量。如果B(s)比A(s)的维数大,则k(s)不为零。第38页/共238页 频域仿真 系统的频域仿真在概念上是非常直观的,但是计算起来常常比较复杂。作为练习,读者可以针对不同的值,计算下面的表达式 G(j)=CjI-A-1 B+D (4.11) 然后通常采用下列三种方法来绘制频域曲

26、线：Bode图、Nichols图和Nyquist曲线。第39页/共238页 MATLAB中的一些函数用来获取系统的频域信号。首先必须产生一个频率的向量。采样点的坐标通常采用对数形式,即从10d1到10d2共n个点,可以通过下面的命令来完成: w=logspace(d1,d2,n) bode和nyquist函数可以用来计算每一个频率所对应的G( j)，即 MAG,PHASE=bode(sys,w) RE,IM=nyquist(sys,w)第40页/共238页图4.5典型二阶系统的Bode幅值曲线 第41页/共238页图4.6 典型二阶系统的Bode频率曲线 第42页/共238页图 4.7第43页

27、/共238页图 4.8第44页/共238页 仿真框图及系统建立 许多设计系统都可以由一些基本的组件和框图中的反馈回路组成。在有些情况下,寻找系统的等价描述和相应的状态空间矩阵是比较困难的。幸运的是,我们可以借助MATLAB从系统的Simulink仿真框图直接建立它的状态空间描述。这些工作可以通过MATLAB中的控制工具箱或Simulink的图形仿真界面来完成。 第45页/共238页 为了演示这一过程,同样考虑一下例4.1所述的简单机械系统。首先建立该系统的Simulink仿真框图,然后自动创建原系统的状态空间和整个系统的传递函数形式。下面将状态方程展开,并且进行Laplace变换1212212

28、211( )( )11( )( )( )dxxX sXsdtsdkcxxxudtmmmkmmXsX sU sccssmm (4.12) 第46页/共238页 上述拉普拉斯变换可以对应于如图4.10所示的基本模块。 现在将这些模块连接起来,定义输出为Y(s)=X1(s),最后得到如图4.2所示的仿真框图。下面我们可以使用linmod函数来计算LTI系统的状态矩阵： A,B,C,D=linmod(mdemosl) 计算的结果与前面得到的结果一致。 第47页/共238页图4.9 MATLAB计算产生的典型二阶系统的Nyquist曲线 第48页/共238页图4.10 Laplace变换下的基本模块第4

29、9页/共238页4.2 控制系统设计 在经典控制系统设计中通常以线性系统模型为研究对象。对于一个线性时不变（LTI）系统,其状态方程可以描述为dXAXBdtyCX(4.13) (4.14) 第50页/共238页 这里已经假定系统的输出没有显式地包含输入变量u（即D=0）。系统（4.14）也可以表示为传递函数形式： Y(s)=G(s)U(s)G(s)=C(sI-A)-1B (4.15) 一个LTI系统的控制系统方框图如图4.11所示。 第51页/共238页图4.11 系统的线性状态方程模型 第52页/共238页 经典比例控制器 在经典控制系统的例子中,首先来看一下图4.12所示的一个简单的闭环系

30、统。对于一个SISO系统而言,系统传递函数G(s)仅仅是式（4.15）所示的标量函数,该传递函数嵌入在图4.12所示的方框图中。 反馈回路包含传感器传递函数H(s),而控制器部分只有简单的增益环节Kc组成,rd是闭环系统期望的响应或参考点。 第53页/共238页图4.12 SISO系统的经典比例控制器框图第54页/共238页 图4.12所示闭环系统的传递函数可以写成( )( )( )( )( )( ) ( )( )( )( )( )( )1( )( )ccdcdcdcY sG s K E sG s KRsH s Y sK G sY sRsG s RsK G s H s(4.16) 其中,Gc(

31、s)为闭环传递函数,Kc是经典比例增益。对于单位反馈情况有H(s)=1,Gc(s)可以简化为( )( )1( )cccK G sG sK G s(4.17) 下面是标量输入函数的时域表示 u(t)=Kc(rd(t)-y(t)=Kc(rd(t)-CTX(t) (4.18)第55页/共238页 从而式（4.14）可以写成()TTcdcccddXAXK BrK BC XAK BCXK Brdt(4.19) 这里的参考点rd成为系统的一个独立输入变量。第56页/共238页 既然控制器只有唯一的参数Kc需要确定,因此该系统的控制器设计比较简单。闭环系统的暂态响应由状态方程系数矩阵的特征值或者整个系统的根

32、极点确定。我们可以在时域中通过选择合适的控制参数Kc,使得(A-KcBCT)的特征值产生期望的暂态响应（上升时间、最大超调量等）。与此类似,也可以在传递函数中通过选择合适的控制参数Kc来设计式（4.17）系统的根极点位置。这两种设计方法是等价的。我们知道Gc(s)的极点是1+KcG(s)的根,因此可以将极点配置方程看作控制增益Kc的根。运用根轨迹方法可以确定满足设计要求的控制参数。 第57页/共238页 状态反馈控制器 上述经典控制器的主要不足是系统仅有唯一的控制参数Kc可供调整,而对于N维控制系统,系统开环矩阵具有N个特征值或者开环传递函数具有N个极点,即 det(A-I)=0 或 det(

33、sI-A)=0 (4.20) 要想将所有这些系统根极点调整到需要的位置,控制器至少需要N个独立变量,因此仅仅将系统输出信号进行反馈将不能满足控制器设计的要求。一个自然的想法就是将系统的所有状态变量X都进行反馈,这就产生了状态反馈控制器。 第58页/共238页 对于SISO系统,状态反馈后的系统输入变成 u(t)=rd(t)-KTsX(t) (4.21) Ks称为系统的反馈系数。 这样,闭环系统的状态方程可以写成()TTdssddXAXBrBK XABKXBrdt(4.22) 闭环系统的框图如图4.13所示。同时,图4.11所示的状态反馈系统变成图4.14所示的仿真框图。 第59页/共238页

34、图4.13 SISO系统的状态反馈控制器第60页/共238页图 4.14第61页/共238页 完全可控性 为了设计具有状态观测器的状态反馈控制器,让我们首先熟悉有关系统可控性的定义。 假设一个SISOLTI系统由式（4.23）描述TdXAXBuyC Xdt(4.23) 第62页/共238页 如果该系统能够构造一个无约束的输入信号u(t),使得系统能够在有限的时间间隔内（t0ttf）由初始状态运动到任何其它的状态,则可以说系统在t0时刻是可控的。如果系统的每个状态都是可控的,则称该系统是完全可控的。 不失一般性,假设X(tf)=0,t0=0,则()0( )(0)( )tAtA tX te XeB

35、ud(4.24) 根据完全可控性的定义,有 ()0( )0(0)( )ffftAtA tX teXeBud(4.25) 第63页/共238页或者 01011000(0)( )( )(0)( ) ( )fftANAkkkNNtkkkkkkXe BudeaAXA BaudA B 根据Sylvester积分公式 有 (4.26) (4.27) (4.28) 第64页/共238页 当如下矩阵非奇异时,系统满足完全可控的条件： M=B AB A2BAN-1 B (4.30)或者012121(0)NNXBABA BAB(4.29) 第65页/共238页 极点配置 设计状态反馈控制器的最简单方法是采用极点配

36、置。其基本思想是首先确定闭环系统N个根极点的期望位置,然后设计适当的反馈增益,从而将系统的极点调整到期望的位置。 如果系统是完全可控的,则这一过程完全可以表示成包含N个未知参数的N个方程组的求解。所需要设计的反馈控制增益就是该方程组的解。第66页/共238页 如果系统比较简单,则完全可以通过手工计算完成系统的极点配置,但无论是手工计算,还是通过MATLAB函数自动计算,其基本步骤都是相同的，如下所示： （1）检查系统的可控矩阵是否满秩。 （2）确定闭环系统的期望极点,1,2,N。 （3）确定希望配置的极点位置后,可以建立期望的特征方程。1121()()()0NNNNsssss第67页/共238

37、页 （4）最后建立闭环系统的特征方程,即(sI-(A-BKTs)=0,将（3）、（4）步建立的方程联立,由于其多项式的系数相等,由此可以建立N个位置参数的N个方程组,从而可以唯一地确定系统的反馈增益矩阵KTs。第68页/共238页 例4.2 假定SISOLTI系统的状态方程为01020.61dXAXBuABsdt 闭环系统的期望极点为1,2=-1.82.4j,试设计确定系统状态反馈的增益矩阵。 解：首先观察开环系统的极点 2120.6020.6ssIAss第69页/共238页 可以看出,系统开环极点为s1,2=4.539,系统是不稳定的。 闭环系统的期望极点是由期望的系统暂态响应特性（上升时间

38、、读者可以验证,1 , 2=-1.82.4j的闭环极点将产生较好的动态特性（大约10%的最大超调量和大约0.6s的上升时间）。因此,期望的闭环系统极点是不唯一的。 下面在已经确定期望闭环系统极点的情况下来设计系统的反馈增益矩阵。 Step1：验证系统的可控性。第70页/共238页 矩阵M的秩等于N,因此系统满足完全可控性条件。 Step23：计算期望的特征方程 (s-1)(s-2)=(s+1.8-j2.4)(s+1.8+j2.4)=s2+1s+2=0 Step4：计算闭环系统的特征方程010001111020.60MBAB det()(0TTsssIABKsIABK第71页/共238页1212

39、1200010120.6TsTsBKkkkkABKkk 因此 221121()20.6020.6Tsss IABKsk skksk第72页/共238页 带全观测器的状态反馈控制 设计状态反馈控制器的主要问题是要求系统的所有状态变量都是可测的。然而对于一个实际系统而言,有些状态的信号值很难测量甚至不可能直接通过传感器进行测量,或者虽然可以进行直接测量,但在经济上却要增加相应的成本。这样,如果不能得到系统的全状态向量,前面讲述的状态反馈控制就不可能实现。第73页/共238页 解决以上问题的方法是利用系统某种数学形式的仿真来估计不能测量的状态值,这种方法称之为系统的状态观测器设计。 下面假定以SIS

40、OLTI系统为研究对象,这意味着系统有唯一的可控变量和唯一的可测量。同时,假定系统输出y(t)是唯一能够测量的量,它将被引入到状态观测器中来提高状态值的估计过程。这里采用 来表示状态向量X(t)的在t时刻的估计值。( )X t第74页/共238页 状态观测器的框图如图4.15所示（注意变量xc表示 。该观测器使用u(t)和y(t)作为输入量,并且输出系统状态关于时间的估计值。从框图中可以看出( )X t( )()()TdX tAXBuL yC XALCXBuLydt(4.31) 这里的L为未知的增益,它是根据该子系统期望的暂态响应特性确定的,称为状态观测器的增益矩阵。对于SISO系统,L是长度

41、为N的列向量。 第75页/共238页图4.15 SISO系统的状态观测器模型 第76页/共238页 观测器的设计过程与前面讲述的标准状态反馈控制器类似。这里的观测器的增益选择应使状态观测器的特征值是稳定的,同时使得观测器的动态变化快于整个闭环系统的动态属性。观测器的特征值由下式确定： det(sI-(A-LCT)=0 (4.32) 在状态反馈控制系统中加入状态观测器,可得到图4.16所示的系统框图。对于该系统,系统输入为 ( )( )( )Tdsu tr tKX t(4.33) 第77页/共238页 如果系统模型与状态观测器模型都采用相同的状态空间矩阵A、B、CT进行描述,则对于被研究对象有T

42、dXAXBuyC Xdt 将式（4.33）的输入代入上面的方程,则得到系统的完整模型(4.34) TdsdXAXBrBKXdt(4.35) 第78页/共238页图4.16 具有全状态观测器的状态反馈 控制系统框图（SISO系统） 第79页/共238页 对于状态观测器,将式（4.33）代入式（4.31）可以得到观测器的完整模型()()TTTdsTTTsddXALCXBrBKXLC XdtdXABKLCXBrLC Xdt或者 (4.36) (4.37) 定义误差向量 ()TEXXdE ALCEdt将式（4.37）代入得到误差向量的动态模型 (4.38) (4.39) 第80页/共238页 完全可观

43、性 如果系统的每个状态X(t0)都可以通过y(t)一段时间的观测值来确定,则该系统被称为是完全可观的。考虑式（4.23）定义的SISOLTI系统,其时域解为()0( )(0)( )tTAtTA ty tC e XCeBud(4.40) 假设u()=0,为方便计算,上式可以写成( )(0)TAty tC e X(4.41)第81页/共238页 其中,CTeAt已知,y(t)可测。因此状态向量X(0)可以通过y(t)的观测值间接计算得到。 对于SISO系统,方程（4.41）仅有一个方程,但包含N个未知参数。然而,由于该方程独立于时间变量,因此,在多个时刻对y(t)进行测量,可以获得多个类似的方程,

44、将它们联立,就可以唯一确定系统的初始状态X(0)。 下面给方程（4.41）两边同时乘以已知的系数矩阵,得到)()(0)()()(0)TTTTTAtTTAtTTAtTAtTAtTTA tA tA tTAtC eyC eC e XC eeCeCeCyeCC e X(4.42) (4.43) (4.44)变换方程有 第82页/共238页 将整个观测时间内的所有方程进行联立,得到(0)QWX(4.45) 其中 00( )TTffttAATAQeCydWeCC e d(4.46) 最后求解方程（4.45）,得到 1(0)XW Q 如果W是非奇异的,X(0)可以由y(t)的观测值唯一确定,从而系统是完全可

45、观的。 再次使用Sylvester的积分公式,得到第83页/共238页21201101011011( )( )TTNTTTAtTkTkNkTNNA tkTTTkkNCC AC eNt C A aaaC AC AeCt A CCA CA CAC定义 21NTTTTHCA CA CAC(4.48) (4.49) (4.50) 第84页/共238页 观测器增益的确定 确定观测增益矩阵同样采用极点配置的方法。然而在这里,我们确定的是状态观测器误差方程的极点位置。误差极点位置的选择比较随意,但误差动态变化应该比被控系统的动态变化快一些。如果系统完全可观,则(A-LCT)的N个特征值的位置应该唯一确定观测

46、器增益矩阵的N个元素。设计的过程如下： （1）检查系统可观矩阵是否奇异。第85页/共238页 （2）为误差方程 指定期望的极点位置（1,2,N）。这些极点位置与系统的主导极点相比较应该更靠近复平面的左手边。 （3）根据期望的极点位置创建期望的特征方程。 （4）最后创建误差方程的特征方程,从而得到含有N个位置参数的方程组。 同样，以例4.2为例来说明SISO系统状态观测器的设计过程。 EXX第86页/共238页 例4.4用MATLAB的place函数重新设计例4.2的状态反馈控制器与例4.3中的全状态观测器。 解：程序代码如下： %SFSOTEST.MSISOLTI系统的状态反馈控制器与全状态观

47、测器的设计 clearall,closeall,nfig=0; %打开二进制文件保存结果 deletesfsotest.out diarysfsotest.out disp()第87页/共238页 disp() %PartI.创建线性系统模型,显示其开环系统是不稳定的Setupbasedataforthelinear A=01;20.60;B=01;%建立系统的状态空间模型 C=10;D=0; disp(StateSpaceMatricesforthePlant) A,B,C,D第88页/共238页 %compute eigenvalues of state matrix for open l

48、oop plant disp(Eigenvalues of theOpen Loop Plant);%计算开环系统状态矩阵特征值 ev=eig(A) %PartII.加入状态反馈控制器以稳定系统,对状态1的输出量进行仿真 %检查系统的可控性 disp(Controllability Matrix for thissystem),M=ctrb(A,B) disp(RankofControllabilityMatrix),rank(M) clp=-1.8+2.4j-1.8-2.4j;%计算状态反馈控制增益 Ks=place(A,B,clp);第89页/共238页 disp(Desiredclose

49、dlooppolesforstatefeedbackcontroller);clpdisp(Statefeedbackgainsneededtogivedesiredpoles);Ks disp(Calculatedeigenvaluesofsystemwithstatefeedback); eig(A-B*Ks) Nv=-1.0/(C*inv(A-B*Ks)*B);%计算Nv disp(SetpointgainforzeroSSerror);Nv to=0;tf=5;nfig=0;%对被控系统+控制器进行仿真 t=linspace(to,tf,101); syscl1=ss(A-B*Ks,B

50、*Nv,C,D); y1,t,x1=step(syscl1,t);第90页/共238页 nfig=nfig+1;figure(nfig)%绘制相关结果曲线 subplot(2,1,1),plot(t,x1(:,1),r-,t,x1(:,2),g-),grid, title(StatesforStateFeedbackTestCase) xlabel(Time),ylabel(StateVariables) legend(x1(t),x2(t) %PartIII.加入状态反馈控制器与全状态观测器,仿真状态1的阶跃响应特性 %检查系统的完全可观性 disp(ObservabilityMatrixf

51、orthissystem),H=obsv(A,C) disp(RankofObservabilityMatrix),rank(H)第91页/共238页 %针对指定的观测器极点计算观测器增益 op=3*clp;%观测器的速度是闭环系统的3倍 L=place(A,C,op);L=L; disp(Desiredobserverpolesforstatefeedbackcontroller);op disp(Estimatorgainsneededtogivedesiredpoles);L disp(Calculatedeigenvaluesofestimatorsystem);eig(A-L*C)

52、A11=A;A12=-B*Ks;B1=B*Nv; A21=L*C;A22=A-L*C-B*Ks;B2=B*Nv; zz=0; AB=A11A12;A21A22;BB=B1;B2; CB=Czz*C;第92页/共238页 %对控制系统进行仿真 syscl2=ss(AB,BB,CB,D);y2,t,x2=step(syscl2,t); %分离系统状态与观测器的估计状态 nn=max(size(A); xp2=x2(:,1:nn);xe2=x2(:,nn+1:2*nn); %绘制全观测器情况下的计算结果 subplot(2,1,2),plot(t,xp2(:,1),r-,t,xp2(:,2),g-)

53、,grid, title(StatesforStateFeedbackwithFullObserverTestCase) xlabel(Time),ylabel(StateVariables) legend(x1(t),x2(t)第93页/共238页 %绘制状态估计的误差 nfig=nfig+1;figure(nfig) plot(t,xp2(:,1)-xe2(:,1),r-,t,xp2(:,2)-xe2(:,2),g-),gri d, title(DifferenceBetweenPlantandObserverStates) xlabel(Time),ylabel(ErrorinState

54、Variables) legend(e1(t),e2(t) Diaryoff %关闭二进制文件 第94页/共238页图4.17 系统模型与观测器动态性能演示 第95页/共238页图4.17 系统模型与观测器动态性能演示 第96页/共238页图4.18 状态误差的动态曲线 第97页/共238页 对偶系统 回忆前面论述的系统可控与可观性的定义，可知系统的可控性要求其可控矩阵满秩,反映的是状态矩阵A与输入矩阵B之间的关系；而系统的可观性要求其可观矩阵满秩,反映的是状态矩阵A与输出矩阵CT之间的关系。这两个概念实际上反映了控制系统的对偶原则。 对于一般的MIMO系统 *dXAXBUdtdZA ZC V

55、dtYCXWB Z(4.52) (4.54) (4.53) 第98页/共238页 对于系统1,其完全状态可控的充要条件是M矩阵满秩，即 M1=B AB A2B AN-1B (4.55) 系统完全状态可观性的充要条件是H*矩阵满秩，即 H*1=C* A* C* A* 2 C*A*N-1 C* (4.56) 对于系统2,其完全状态可控的充要条件是M矩阵满秩，即 M2=C* A* C* A*2 C* A*N-1C* (4.57)第99页/共238页 系统完全状态可观性的充要条件是H*矩阵满秩，即 H*2=B AB A2BAN-1B (4.58) 综上所述,给定系统的可观性可以通过其对偶系统的可控性来

56、检验,而研究系统的可控性则可以通过其对偶系统的可观性来研究,这些性质称为系统的对偶原则。第100页/共238页4.3 控制系统的时域仿真 前一节以经典控制器和状态反馈控制器为例讲述了控制系统控制参数的设计过程。设计过程主要依据的是系统的时域特性,即闭环系统状态矩阵的特征值决定了闭环系统的暂态响应特性。控制器设计的目标是选择控制器的增益,使得闭环系统状态矩阵的特征值位于期望的极点位置。 第101页/共238页 控制器参数确定以后,下一步需要对闭环系统进行仿真。在设计过程中研究的对象一般是系统的线性化模型。而在仿真过程中,应该尽可能准确地再现实际的系统模型,这就常常要求以实际的时变或非线性系统为仿

57、真的对象。这一节主要论述控制器参数确定以后的系统仿真过程。为简单起见,仍然以SISO系统为研究对象。其线性与非线性模型分别为 TdXAXBuyC Xdt(4.59) (, , )TdXF X u tyC Xdt(4.60) 第102页/共238页 经典比例控制器设计 单位反馈回路的简单比例控制系统框图如图4.19所示。系统控制输入为()cduK ry(4.61) 使用线性化模型的闭环仿真方程为()TTcdcccddXAXK BrK BC XAK BCXK Brdt(4.62) 写成标准的状态方程形式 dTTccdXAXBryCXdtAAK BCBK BCC(4.63) 第103页/共238页

58、如果采用非线性模型进行仿真,必须使用MATLAB中标准的ODE求解器来完成。为此,用户必须定义一个函数文件,MATLAB中的ODE求解器调用该函数来完成非线性系统的仿真。该函数文件包含下面的内容： （1）指定t时刻的参考输入rd(t)。 （2）计算系统t时刻的输出y(t)=CTX(t)。 （3）确定t时刻的输入u(t)=Kc(rd(t)-y(t)。 （4）计算t时刻状态变量的导数 ( )( ), ( ), )dX tF X t u t tdt第104页/共238页图4.19 SISO系统的经典比例单位反馈控制框图 第105页/共238页 状态反馈控制器 带全状态观测器的状态反馈控制框图如图4.

59、20所示。该框图类似于图4.16,不同之处在于这里的模型包括一个附加的稳态状态增益模块,其中包含一个归一化的增益变量Nr。 我们再来看看该系统线性化模型与非线性模型的仿真方程,并且重新调整控制器的增益（KTs,L和Nr等）。控制规律写成Tr dsuN rKX(4.64) 第106页/共238页 该系统的闭环模型包括实际的系统状态X(t)以及估计的系统状态 。这样,该系统具有2N个未知参数。对于线性系统模型,该系统的完整模型为被控系统模型( )X tTr dsdXAXBN rBKXdt(4.65) 观测器模型 ()TTTsr ddXABKLCXBN rLC Xdt(4.66) 写成标准的状态方程

60、形式为 ddZAZBryC Zdt第107页/共238页图4.20 具有全状态观测器的状态反馈控制系统 第108页/共238页 混合控制器 下面将研究一种将状态反馈控制运用到经典控制器中的例子。在这个例子中,状态反馈制器为经典控制器提供修正的参考信号。例如,对于图4.18的经典控制系统,如果选择的增益Kc不满足系统的设计要求,我们可能会手动调节参考点rd,以使系统具有更好的暂态响应特性。可以设想一种控制方案来自动为参考信号rd(t)提供实时的修正信号rdm(t),控制框图如图4.20所示。 第109页/共238页 图中左边的sum模块的输出为修正后的参考点信号,在此之后,该系统看上去就像是具有

61、单位反馈回路和比例增益Kc的经典控制器。在这个控制方案中,状态观测器和状态反馈回路的目标是提供一个修正的参考信号,来改善经典控制系统的控制效果。这种方法称为状态反馈辅助的经典控制方法（SFACC）,也称为混合的控制器。 第110页/共238页 图4.21显示的控制系统实际上可以看作是前面讲述的状态反馈控制系统,其控制框图可以认为是从图4.20的标准形式通过某些模块变换发展而来。图4.224.24显示了其中的变换,主要步骤包括: （1）在系统输入前加上一个比例增益模块,并且在相应的引入回路中抵消它的影响（如图4.22所示）。第111页/共238页图4.21 具有嵌入经典控制器的状态反馈控制框图（

62、SFACC） 第112页/共238页 （2）对系统输出y构造一个负反馈,它可以对估计输出 的正反馈加以抵消（如图4.23所示）。 （3）分离出经典控制器部分（图4.24中的阴影部分）。 （4）最后定义修正的状态反馈增益Ksm和修正的稳态状态归一因子Nnn。yTTssmcrnncKKCKNNK(4.68) (4.69) 第113页/共238页图4.22 在SFC中加入经典控制增益 第114页/共238页图4.23 加入经典负反馈回路 第115页/共238页图4.24 分离出经典比例控制器部分 第116页/共238页 图4.21所示的SFACC控制系统与传统的状态反馈控制系统SFC（如图4.20所

63、示）在功能上是一致的。仅有的近似存在于第(2)步,这一步假设估计的输出与实际的系统输出完全抵消。 如果状态观测器设计成相对控制系统有足够快的反应时间,这种近似是可以满足要求的（对于线性模型）。 第117页/共238页 从以上的讨论中可以看出,SFC和SFACC具有相同的暂态响应。然而,SFACC算法比前者具有更好的控制效果。并且SFACC可以方便地在原先经典控制系统的基础上进行扩展,而不需要破坏原来已经存在的控制系统。这让用户可以从经典控制器直接过渡到现代控制系统。 为了对SFACC控制系统进行仿真,我们假定控制规律为()cdmTdmnn dsmTcnn dcsmcuK ryrN rKXuK

64、N rK KXK y修正后的参考信号为 将式（4.71）代入式（4.70），有 （4.70） (4.71) (4.72) 第118页/共238页 与SAC控制系统类似,SFACC闭环系统既包含实际状态也包含估计状态,因此,可以获得2N个未知参数。对于线性系统模型,有 系统模型为()()()()()Tcnn dcsmcTTccsmcnn dTTTTcnn dcsmcdXAXBuAXB K N rK KXK ydtdXAK BCXK BKXK BN rdtdXALCXBuLydtALCXB K N rK KXK yLC X(4.73) (4.74) (4.75) 第119页/共238页或者 ()(

65、)TTTTcsmccnn ddXALCK BKXLCK BCXK BN rdt(4.76) 写成标准的状态方程形式为 ddZAZBryCZdt其中 11122122111221220cnnTcnnTTccsmTTTTccsmXK BNAAZABC CK BNAAXAAK BCAK BKALCK BCAAK BKLC (4.78) 第120页/共238页 参考输入增益的调整 在前面的框图中,为了减小比例控制器带来的稳态误差,引入了一个输入归一化模块。该模块的增益值可以根据下面的算式进行确定（图4.25显示了系统的简化模型）。 线性模型TTr dsdXAXBuyC XdtuN rK X(4.79)

66、 (4.80) 控制律 第121页/共238页 系统稳态时,我们希望yss与参考点rdss相同。简化系统稳态时的状态方程为0()Tssssssr dsssssTsssr dssAXBuAXBN rBK XABKXBN r(4.81) 求解该方程,得稳态时的状态向量为1()Tsssr dsXABKBN r (4.82) 第122页/共238页图4.25 具有SS增益模块的状态反馈控制系统第123页/共238页 稳态输出为1()TTsssr dsyCABKBN r (4.83) 最后将yss=rdss代入,得到Nr的计算公式为11()rTTsNCABKB(4.84) 第124页/共238页4.4 实例:倒摆系统的建模与仿真 引言 这一节我们先看一个实际的例子。图4.26是某个倒摆系统的示意图,倒摆通过转动关节安装在驱动小车上,通过对小车施加一定的驱动,使倒摆保持一定的姿态。这是姿态控制问题的典型代表。第125页/共238页图4.26 驱动小车上的倒摆示意图 第126页/共238页 基本方程 假定倒摆由无质量的轻杆和质量为m的小球组成,小车的质量为M,系统所受的外力包括小球受到的重力和对小车