1、第4章 控制系統的設計與仿真 4.1 系統建模與仿真框圖的創(chuàng)建4.2 控制系統設計4.3 控制系統的時域仿真4.4 實例:倒擺系統的建模與仿真 4.1 系統建模與仿真框圖的創(chuàng)建 4.1.1 問題描述 眾所周知,現實世界中存在著各種不同的控制系統。對于線性時不變（LTI）系統,一般可以分為連續(xù)和離散系統。MATLAB中為用戶提供了豐富的針對各種系統的建模手段。圖4.1顯示了MATLAB中各種線性時不變（LTI）系統之間的轉換關系。 圖4.1 連續(xù)與離散系統的關系示意圖 圖4.1中顯示了MATLAB可以完成離散和連續(xù)系統的建模,并且同一系統可以表示成連續(xù)系統,也可以表示成離散系統,它們之間可以以狀

2、態(tài)方程形式進行轉化。這一節(jié)將結合一個具體實例來演示MATLAB中各種模型創(chuàng)建和相互之間進行轉化的方法,以及如何用Simulink進行連續(xù)系統的仿真。首先給出實例的源程序MODLDEMO.M,然后根據不同的主題結合實例進行講述。 例4.1 對于Mass Spring Dashpot機械系統： (4.1) 試建立該系統的連續(xù)和離散模型,并進行時域和頻域仿真。 解：程序源代碼如下： %MODLDEMO.M 演示各種建模與仿真（時域和頻域）技巧 clearall,closeall %程序開始,清空工作空間,關閉所有窗口 deletemodldemo.out,diarymodldemo.out %打開二

3、進制文件 disp(*MODLDEMO.OUT*DiaryFileforMODLDEMO.M),disp()m=1 %各種系統常數定義k=1 %單位kg/s2c=2.02.51.20.0 %單位kg/s %創(chuàng)建系統的狀態(tài)空間模型disp(StateSpaceModels)km=k/m;A1=01;-km-c(1)/m,A2=01;-km-c(2)/mA3=01;-km-c(3)/m,A4=01;-km-c(4)/mB=01/m,C=10,D=0sys1s=ss(A1,B,C,D);sys2s=ss(A2,B,C,D);sys3s=ss(A3,B,C,D);sys4s=ss(A4,B,C,D);

4、 %仿真系統的脈沖和階躍響應（時域）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),gridxlabel(Time),ylabel(SystemResponse)subplot(222),plot(t,y2,r),title(OverDampingImpulse),gridxlabel(Time),ylabel(SystemRespon

5、se) subplot(223),plot(t,y3,r),title(UnderDampingImpulse),gridxlabel(Time),ylabel(SystemResponse)subplot(224),plot(t,y4,r),title(NoDampingImpulse),gridxlabel(Time),ylabel(SystemResponse) %y1=step(sys1s,t);y2=step(sys2s,t);y3=step(sys3s,t);y4=step(sys4s,t);figure(2)subplot(221),plot(t,y1,r),title(Crit

6、icalDampingStep),gridxlabel(Time),ylabel(SystemResponse)subplot(222),plot(t,y2,r),title(OverDampingStep),gridxlabel(Time),ylabel(SystemResponse)subplot(223),plot(t,y3,r),title(UnderDampingStep),gridxlabel(Time),ylabel(SystemResponse)subplot(224),plot(t,y4,r),title(NoDampingStep),gridxlabel(Time),yla

7、bel(SystemResponse)disp(hitanykeytocontinue),pause%在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(hitanykeytocontinue),pause%將傳遞函數形式轉化成零極點形式disp(ZeroPoleGainForm)sys1z=zpk(sys1t),sys2z=zpk(sys2t)s

8、ys3z=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); %注意tfdata命令定義元胞數組%將元胞數組轉化成向量形式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=r

9、esidue(n2,d2)r3,pr3,kr=residue(n3,d3),r4,pr4,kr=residue(n4,d4)disp(hitanykeytocontinue),pause %進行頻域仿真,這里使用nyquist函數（參考bode,freqs等函數）w=logspace(-2,2,100);re,im=nyquist(sys1s,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)

10、;%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,1,:);g3=re3+i*im3; 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=re

11、4+i*im4;mag4=20*log10(abs(g4);phase4=angle(g4)*180/pi;%幅值曲線figure(3)subplot(221),semilogx(w,mag1,r),title(CriticalDamping),xlabel(frequency),ylabel(|G(jw)|indb)subplot(222),semilogx(w,mag2,r),title(OverDamping),xlabel(frequency),ylabel(|G(jw)|indb)subplot(223),semilogx(w,mag3,r),title(UnderDamping),

12、xlabel(frequency),ylabel(|G(jw)|indb)subplot(224),semilogx(w,mag4,r),title(NoDamping),xlabel(frequency),ylabel(|G(jw)|indb)disp(hitanykeytocontinue),pause%頻率曲線figure(4)subplot(221),semilogx(w,phase1,r),title(CriticalDamping),xlabel(frequency),ylabel(angle)subplot(222),semilogx(w,phase2,r),title(Over

13、Damping),xlabel(frequency),ylabel(angle)subplot(223),semilogx(w,phase3,r),title(UnderDamping),xlabel(frequency),ylabel(angle)subplot(224),semilogx(w,phase4,r),title(NoDamping),xlabel(frequency),ylabel(angle)disp(hitanykeytocontinue),pause%繪制系統Nichols圖（對數坐標形式）figure(5)subplot(221),plot(phase1,mag1,r+

14、),title(CriticalDamping),xlabel(angle),ylabel(|G(jw)|indb),gtext(w=0)subplot(222),plot(phase2,mag2,r+),title(OverDamping),xlabel(angle),ylabel(|G(jw)|indb),gtext(w=0)subplot(223),plot(phase3,mag3,r+),title(UnderDamping), xlabel(angle),ylabel(|G(jw)|indb),gtext(w=0)subplot(224),plot(phase4,mag4,r+),t

15、itle(NoDamping),xlabel(angle),ylabel(|G(jw)|indb),gtext(w=0)disp(hitanykeytocontinue),pause%繪制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)subplo

16、t(223),plot(re3,im3,r+),title(UnderDamping),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 %直接使用Nyquist函數繪制Nyquist圖figure(7)subplot(2,2,1),nyquist(sys1s,w),title(CriticalDamping),subplot(2,2,2),ny

17、quist(sys2s,w),title(OverDamping),subplot(2,2,3),nyquist(sys3s,w),title(UnderDamping), subplot(2,2,4),nyquist(sys4s,w),title(NoDamping),disp(hitanykeytocontinue),disp(hitcntrl-ctogetoutofthefileforinteractiveanalysis),pause%使用Simulink框圖（文件名為mdemosl.mdl,如圖4.2所示）進行仿真clearall%清空所有內存m=1,k=1,c=1.2%定義常數di

18、sp(DatafromgraphicalSimulinkmodel)%計算線性化模型A,B,C,D=linmod(mdemosl),syss=ss(A,B,C,D);sysz1=zpk(syss),syssm=minreal(syss),sysz2=zpk(syssm)diaryoff%關閉二進制文件圖4.2 Mass Spring Dashpot系統仿真框圖 4.1.2 物理模型 例4.1 中研究的對象是一個簡單的物質交換機械系統,可以用微分方程表示成(4.2) 其中,y(t)是系統的瞬時交換的質量,k和c分別為比例常數。如果取狀態(tài)變量x1=y和x2=dy/dt,則可以得到系統的狀態(tài)方程形式

19、： 可以看出,矩陣A、B、C和D可以看成LTI系統的標準狀態(tài)方程形式,于是式（4.3）可以寫成(4.3) 4.1.3 時域仿真 一般對控制系統進行時域仿真可以采用impulse、step和lsim函數,這些函數用來處理系統的狀態(tài)空間描述。例如,對于輸入信號u(t),系統的動態(tài)仿真可表示成 sys=ss(A,B,C,D) Y,T,X=lsim(sys,U,t,xo) 當然,這些函數也可以處理系統的頻域表達形式,這時,LTI系統可以通過下式表示 sys=tf(num,den) 作為一個典型的例子,例4.1使用了impulse和step函數來仿真系統的時域特性。圖4.3和圖4.4分別為例4.1執(zhí)行的

20、仿真結果。圖4.3 典型二階系統的脈沖響應曲線圖4.4 典型二階系統的階躍響應曲線 4.1.4 模型轉換 LTI系統的頻域描述可以用傳遞函數表示為 Y(s)=G(s)U(s) (4.5)系統的傳遞函數矩陣為(4.6) 對于例4.1研究的簡單SISO機械系統,其傳遞函數可以寫成簡單的標量形式 (4.7) 但是對于多輸入多輸出（MIMO）系統而言,其傳遞函數描述就有些復雜了,這時MATLAB中的模型轉換函數可以發(fā)揮作用,它可以完成系統在狀態(tài)方程形式與傳遞函數形式之間的互換,同時也可以將傳遞函數形式轉換成零極點-增益形式。相關的函數包括 sys1=ss(A,B,C,D) sys2=tf(sys1)

21、sys3=zpk(sys2) 我們也可以采用ssdata、tfdata和zpkdata等命令將存儲在與一個指定LTI對象相聯的數據結構中的信息抽取出來。例如 num,den=tfdata(sys2) 返回LTI對象sys2的分子和分母多項式系數,num與den為相應的元胞數組,其行數為輸出的維數,列數等于輸入的維數。其中第i行第j列元素表示從第j個輸入到第i個輸出的傳遞函數。 另一種從數據結構中得到元胞數組的方法是使用MATLAB的celldisp命令。簡單的顯示數據的信息,可以使用如下的命令: fieldnames(sys2) num1=sys2.num,den1=sys2.den cell

22、disp(num1),celldisp(den1) 也可以對零極點-增益形式完成顯示的操作。例如 Z,P,K=zpkdata(sys3) 將獲取LTI系統sys3每一個IO通道的零極點和增益大小。元胞數組Z、P和矩陣K 的行與列分別與輸出和輸入的維數相同。其中第i行第j列元素表示從第j個輸入到第i個輸出傳遞函數的零極點和增益。 對于單輸入單輸出（SISO）系統,其傳遞函數與零極點-增益形式可以簡化成普通的分數形式，即(4.8) 4.1.5 留數定理 我們也可以將系統寫成幾個分數相加的形式,例如對于SISO的機械系統,G(s)可以寫成(4.9)(4.10) r1,r2可以通過各種不同的方法計算得

23、到。對于上述問題,MATLAB可用residue函數來完成這一運算，即 r,p,k=residue(B,A) 其中B和A 為包含多項式系數的行向量,而r和p 是包含留數和極點的列向量。如果B(s)比A(s)的維數大,則k(s)不為零。 4.1.6 頻域仿真 系統的頻域仿真在概念上是非常直觀的,但是計算起來常常比較復雜。作為練習,讀者可以針對不同的值,計算下面的表達式 G(j)=CjI-A-1 B+D (4.11) 然后通常采用下列三種方法來繪制頻域曲線：Bode圖、Nichols圖和Nyquist曲線。 MATLAB中的一些函數用來獲取系統的頻域信號。首先必須產生一個頻率的向量。采樣點的坐標通

24、常采用對數形式,即從10d1到10d2共n個點,可以通過下面的命令來完成: w=logspace(d1,d2,n) bode和nyquist函數可以用來計算每一個頻率所對應的G(j)，即 MAG,PHASE=bode(sys,w) RE,IM=nyquist(sys,w)圖4.5典型二階系統的Bode幅值曲線 圖4.6 典型二階系統的Bode頻率曲線 圖 4.7圖 4.8 4.1.7 仿真框圖及系統建立 許多設計系統都可以由一些基本的組件和框圖中的反饋回路組成。在有些情況下,尋找系統的等價描述和相應的狀態(tài)空間矩陣是比較困難的。幸運的是,我們可以借助MATLAB從系統的Simulink仿真框圖直

25、接建立它的狀態(tài)空間描述。這些工作可以通過MATLAB中的控制工具箱或Simulink的圖形仿真界面來完成。 為了演示這一過程,同樣考慮一下例4.1所述的簡單機械系統。首先建立該系統的Simulink仿真框圖,然后自動創(chuàng)建原系統的狀態(tài)空間和整個系統的傳遞函數形式。下面將狀態(tài)方程展開,并且進行Laplace變換(4.12) 上述拉普拉斯變換可以對應于如圖4.10所示的基本模塊。 現在將這些模塊連接起來,定義輸出為Y(s)=X1(s),最后得到如圖4.2所示的仿真框圖。下面我們可以使用linmod函數來計算LTI系統的狀態(tài)矩陣： A,B,C,D=linmod(mdemosl) 計算的結果與前面得到的

26、結果一致。 圖4.9 MATLAB計算產生的典型二階系統的Nyquist曲線 圖4.10 Laplace變換下的基本模塊4.2 控制系統設計 在經典控制系統設計中通常以線性系統模型為研究對象。對于一個線性時不變（LTI）系統,其狀態(tài)方程可以描述為(4.13) (4.14) 這里已經假定系統的輸出沒有顯式地包含輸入變量u（即D=0）。系統（4.14）也可以表示為傳遞函數形式： Y(s)=G(s)U(s)G(s)=C(sI-A)-1B (4.15) 一個LTI系統的控制系統方框圖如圖4.11所示。 圖4.11 系統的線性狀態(tài)方程模型 4.2.1 經典比例控制器 在經典控制系統的例子中,首先來看一下

27、圖4.12所示的一個簡單的閉環(huán)系統。對于一個SISO系統而言,系統傳遞函數G(s)僅僅是式（4.15）所示的標量函數,該傳遞函數嵌入在圖4.12所示的方框圖中。 反饋回路包含傳感器傳遞函數H(s),而控制器部分只有簡單的增益環(huán)節(jié)Kc組成,rd是閉環(huán)系統期望的響應或參考點。 圖4.12 SISO系統的經典比例控制器框圖 圖4.12所示閉環(huán)系統的傳遞函數可以寫成(4.16) 其中,Gc(s)為閉環(huán)傳遞函數,Kc是經典比例增益。對于單位反饋情況有H(s)=1,Gc(s)可以簡化為(4.17) 下面是標量輸入函數的時域表示 u(t)=Kc(rd(t)-y(t)=Kc(rd(t)-CTX(t) (4.1

28、8) 從而式（4.14）可以寫成(4.19) 這里的參考點rd成為系統的一個獨立輸入變量。 既然控制器只有唯一的參數Kc需要確定,因此該系統的控制器設計比較簡單。閉環(huán)系統的暫態(tài)響應由狀態(tài)方程系數矩陣的特征值或者整個系統的根極點確定。我們可以在時域中通過選擇合適的控制參數Kc,使得(A-KcBCT)的特征值產生期望的暫態(tài)響應（上升時間、最大超調量等）。與此類似,也可以在傳遞函數中通過選擇合適的控制參數Kc來設計式（4.17）系統的根極點位置。這兩種設計方法是等價的。我們知道Gc(s)的極點是1+KcG(s)的根,因此可以將極點配置方程看作控制增益Kc的根。運用根軌跡方法可以確定滿足設計要求的控制

29、參數。 4.2.2 狀態(tài)反饋控制器 上述經典控制器的主要不足是系統僅有唯一的控制參數Kc可供調整,而對于N維控制系統,系統開環(huán)矩陣具有N個特征值或者開環(huán)傳遞函數具有N個極點,即 det(A-I)=0 或 det(sI-A)=0 (4.20) 要想將所有這些系統根極點調整到需要的位置,控制器至少需要N個獨立變量,因此僅僅將系統輸出信號進行反饋將不能滿足控制器設計的要求。一個自然的想法就是將系統的所有狀態(tài)變量X都進行反饋,這就產生了狀態(tài)反饋控制器。 對于SISO系統,狀態(tài)反饋后的系統輸入變成 u(t)=rd(t)-KTsX(t) (4.21) Ks稱為系統的反饋系數。 這樣,閉環(huán)系統的狀態(tài)方程可以

30、寫成(4.22) 閉環(huán)系統的框圖如圖4.13所示。同時,圖4.11所示的狀態(tài)反饋系統變成圖4.14所示的仿真框圖。 圖4.13 SISO系統的狀態(tài)反饋控制器圖 4.14 4.2.3 完全可控性 為了設計具有狀態(tài)觀測器的狀態(tài)反饋控制器,讓我們首先熟悉有關系統可控性的定義。 假設一個SISOLTI系統由式（4.23）描述(4.23) 如果該系統能夠構造一個無約束的輸入信號u(t),使得系統能夠在有限的時間間隔內（t0ttf）由初始狀態(tài)運動到任何其它的狀態(tài),則可以說系統在t0時刻是可控的。如果系統的每個狀態(tài)都是可控的,則稱該系統是完全可控的。 不失一般性,假設X(tf)=0,t0=0,則(4.24)

31、 根據完全可控性的定義,有 (4.25) 或者 根據Sylvester積分公式 有 (4.26) (4.27) (4.28) 當如下矩陣非奇異時,系統滿足完全可控的條件： M=B AB A2BAN-1 B (4.30)或者(4.29) 4.2.4 極點配置 設計狀態(tài)反饋控制器的最簡單方法是采用極點配置。其基本思想是首先確定閉環(huán)系統N個根極點的期望位置,然后設計適當的反饋增益,從而將系統的極點調整到期望的位置。 如果系統是完全可控的,則這一過程完全可以表示成包含N個未知參數的N個方程組的求解。所需要設計的反饋控制增益就是該方程組的解。 如果系統比較簡單,則完全可以通過手工計算完成系統的極點配置,

32、但無論是手工計算,還是通過MATLAB函數自動計算,其基本步驟都是相同的，如下所示： （1）檢查系統的可控矩陣是否滿秩。 （2）確定閉環(huán)系統的期望極點,1,2,N。 （3）確定希望配置的極點位置后,可以建立期望的特征方程。 （4）最后建立閉環(huán)系統的特征方程,即(sI-(A-BKTs)=0,將（3）、（4）步建立的方程聯立,由于其多項式的系數相等,由此可以建立N個位置參數的N個方程組,從而可以唯一地確定系統的反饋增益矩陣KTs。 例4.2 假定SISOLTI系統的狀態(tài)方程為 閉環(huán)系統的期望極點為1,2=-1.82.4j,試設計確定系統狀態(tài)反饋的增益矩陣。 解：首先觀察開環(huán)系統的極點 可以看出,系

33、統開環(huán)極點為s1,2=4.539,系統是不穩(wěn)定的。 閉環(huán)系統的期望極點是由期望的系統暫態(tài)響應特性（上升時間、讀者可以驗證,1,2=-1.82.4j的閉環(huán)極點將產生較好的動態(tài)特性（大約10%的最大超調量和大約0.6s的上升時間）。因此,期望的閉環(huán)系統極點是不唯一的。 下面在已經確定期望閉環(huán)系統極點的情況下來設計系統的反饋增益矩陣。 Step1：驗證系統的可控性。 矩陣M的秩等于N,因此系統滿足完全可控性條件。 Step23：計算期望的特征方程 (s-1)(s-2)=(s+1.8-j2.4)(s+1.8+j2.4)=s2+1s+2=0 Step4：計算閉環(huán)系統的特征方程因此 4.2.5 帶全觀測器

34、的狀態(tài)反饋控制 設計狀態(tài)反饋控制器的主要問題是要求系統的所有狀態(tài)變量都是可測的。然而對于一個實際系統而言,有些狀態(tài)的信號值很難測量甚至不可能直接通過傳感器進行測量,或者雖然可以進行直接測量,但在經濟上卻要增加相應的成本。這樣,如果不能得到系統的全狀態(tài)向量,前面講述的狀態(tài)反饋控制就不可能實現。 解決以上問題的方法是利用系統某種數學形式的仿真來估計不能測量的狀態(tài)值,這種方法稱之為系統的狀態(tài)觀測器設計。 下面假定以SISOLTI系統為研究對象,這意味著系統有唯一的可控變量和唯一的可測量。同時,假定系統輸出y(t)是唯一能夠測量的量,它將被引入到狀態(tài)觀測器中來提高狀態(tài)值的估計過程。這里采用 來表示狀態(tài)

35、向量X(t)的在t時刻的估計值。 狀態(tài)觀測器的框圖如圖4.15所示（注意變量xc表示 。該觀測器使用u(t)和y(t)作為輸入量,并且輸出系統狀態(tài)關于時間的估計值。從框圖中可以看出(4.31) 這里的L為未知的增益,它是根據該子系統期望的暫態(tài)響應特性確定的,稱為狀態(tài)觀測器的增益矩陣。對于SISO系統,L是長度為N的列向量。 圖4.15 SISO系統的狀態(tài)觀測器模型 觀測器的設計過程與前面講述的標準狀態(tài)反饋控制器類似。這里的觀測器的增益選擇應使狀態(tài)觀測器的特征值是穩(wěn)定的,同時使得觀測器的動態(tài)變化快于整個閉環(huán)系統的動態(tài)屬性。觀測器的特征值由下式確定： det(sI-(A-LCT)=0 (4.32)

36、 在狀態(tài)反饋控制系統中加入狀態(tài)觀測器,可得到圖4.16所示的系統框圖。對于該系統,系統輸入為 (4.33) 如果系統模型與狀態(tài)觀測器模型都采用相同的狀態(tài)空間矩陣A、B、CT進行描述,則對于被研究對象有 將式（4.33）的輸入代入上面的方程,則得到系統的完整模型(4.34) (4.35) 圖4.16 具有全狀態(tài)觀測器的狀態(tài)反饋 控制系統框圖（SISO系統） 對于狀態(tài)觀測器,將式（4.33）代入式（4.31）可以得到觀測器的完整模型或者 (4.36) (4.37) 定義誤差向量 將式（4.37）代入得到誤差向量的動態(tài)模型 (4.38) (4.39) 4.2.6 完全可觀性 如果系統的每個狀態(tài)X(t

37、0)都可以通過y(t)一段時間的觀測值來確定,則該系統被稱為是完全可觀的。考慮式（4.23）定義的SISOLTI系統,其時域解為(4.40) 假設u()=0,為方便計算,上式可以寫成(4.41) 其中,CTeAt已知,y(t)可測。因此狀態(tài)向量X(0)可以通過y(t)的觀測值間接計算得到。 對于SISO系統,方程（4.41）僅有一個方程,但包含N個未知參數。然而,由于該方程獨立于時間變量,因此,在多個時刻對y(t)進行測量,可以獲得多個類似的方程,將它們聯立,就可以唯一確定系統的初始狀態(tài)X(0)。 下面給方程（4.41）兩邊同時乘以已知的系數矩陣,得到(4.42) (4.43) (4.44)變

38、換方程有 將整個觀測時間內的所有方程進行聯立,得到(4.45) 其中 (4.46) 最后求解方程（4.45）,得到 如果W是非奇異的,X(0)可以由y(t)的觀測值唯一確定,從而系統是完全可觀的。 再次使用Sylvester的積分公式,得到定義 (4.48) (4.49) (4.50) 4.2.7 觀測器增益的確定 確定觀測增益矩陣同樣采用極點配置的方法。然而在這里,我們確定的是狀態(tài)觀測器誤差方程的極點位置。誤差極點位置的選擇比較隨意,但誤差動態(tài)變化應該比被控系統的動態(tài)變化快一些。如果系統完全可觀,則(A-LCT)的N個特征值的位置應該唯一確定觀測器增益矩陣的N個元素。設計的過程如下： （1）

39、檢查系統可觀矩陣是否奇異。 （2）為誤差方程 指定期望的極點位置（1,2,N）。這些極點位置與系統的主導極點相比較應該更靠近復平面的左手邊。 （3）根據期望的極點位置創(chuàng)建期望的特征方程。 （4）最后創(chuàng)建誤差方程的特征方程,從而得到含有N個位置參數的方程組。 同樣，以例4.2為例來說明SISO系統狀態(tài)觀測器的設計過程。 例4.4用MATLAB的place函數重新設計例4.2的狀態(tài)反饋控制器與例4.3中的全狀態(tài)觀測器。 解：程序代碼如下： %SFSOTEST.MSISOLTI系統的狀態(tài)反饋控制器與全狀態(tài)觀測器的設計 clearall,closeall,nfig=0; %打開二進制文件保存結果 de

40、letesfsotest.out diarysfsotest.out disp()disp(*SFSOTEST.OUT*DiaryFileforSFSOTEST.M)disp()%PartI.創(chuàng)建線性系統模型,顯示其開環(huán)系統是不穩(wěn)定的SetupbasedataforthelinearA=01;20.60;B=01;%建立系統的狀態(tài)空間模型C=10;D=0;disp(StateSpaceMatricesforthePlant)A,B,C,D%compute eigenvalues of state matrix for open loop plantdisp(Eigenvalues of the

41、Open Loop Plant);%計算開環(huán)系統狀態(tài)矩陣特征值ev=eig(A)%PartII.加入狀態(tài)反饋控制器以穩(wěn)定系統,對狀態(tài)1的輸出量進行仿真%檢查系統的可控性disp(Controllability Matrix for thissystem),M=ctrb(A,B)disp(RankofControllabilityMatrix),rank(M)clp=-1.8+2.4j-1.8-2.4j;%計算狀態(tài)反饋控制增益Ks=place(A,B,clp);disp(Desiredclosedlooppolesforstatefeedbackcontroller);clpdisp(State

42、feedbackgainsneededtogivedesiredpoles);Ksdisp(Calculatedeigenvaluesofsystemwithstatefeedback);eig(A-B*Ks)Nv=-1.0/(C*inv(A-B*Ks)*B);%計算Nvdisp(SetpointgainforzeroSSerror);Nvto=0;tf=5;nfig=0;%對被控系統+控制器進行仿真t=linspace(to,tf,101);syscl1=ss(A-B*Ks,B*Nv,C,D);y1,t,x1=step(syscl1,t);nfig=nfig+1;figure(nfig)%繪

43、制相關結果曲線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.加入狀態(tài)反饋控制器與全狀態(tài)觀測器,仿真狀態(tài)1的階躍響應特性%檢查系統的完全可觀性disp(ObservabilityMatrixforthissystem),H=obsv(A,C)disp(RankofObservabilityMatrix),rank(H)%針對指定的觀測器極點計算

44、觀測器增益op=3*clp;%觀測器的速度是閉環(huán)系統的3倍L=place(A,C,op);L=L;disp(Desiredobserverpolesforstatefeedbackcontroller);opdisp(Estimatorgainsneededtogivedesiredpoles);Ldisp(Calculatedeigenvaluesofestimatorsystem);eig(A-L*C)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;

45、%對控制系統進行仿真syscl2=ss(AB,BB,CB,D);y2,t,x2=step(syscl2,t);%分離系統狀態(tài)與觀測器的估計狀態(tài)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-),grid,title(StatesforStateFeedbackwithFullObserverTestCase)xlabel(Time),ylabel(StateVariables)legend(x1(t),x2(t)%繪制狀

46、態(tài)估計的誤差nfig=nfig+1;figure(nfig)plot(t,xp2(:,1)-xe2(:,1),r-,t,xp2(:,2)-xe2(:,2),g-),grid,title(DifferenceBetweenPlantandObserverStates)xlabel(Time),ylabel(ErrorinStateVariables)legend(e1(t),e2(t)Diaryoff %關閉二進制文件 圖4.17 系統模型與觀測器動態(tài)性能演示 圖4.17 系統模型與觀測器動態(tài)性能演示 圖4.18 狀態(tài)誤差的動態(tài)曲線 4.2.8 對偶系統 回憶前面論述的系統可控與可觀性的定義，可

47、知系統的可控性要求其可控矩陣滿秩,反映的是狀態(tài)矩陣A與輸入矩陣B之間的關系；而系統的可觀性要求其可觀矩陣滿秩,反映的是狀態(tài)矩陣A與輸出矩陣CT之間的關系。這兩個概念實際上反映了控制系統的對偶原則。 對于一般的MIMO系統 (4.52) (4.54) (4.53) 對于系統1,其完全狀態(tài)可控的充要條件是M矩陣滿秩，即 M1=B AB A2B AN-1B (4.55)系統完全狀態(tài)可觀性的充要條件是H*矩陣滿秩，即 H*1=C* A* C* A* 2 C*A*N-1 C* (4.56) 對于系統2,其完全狀態(tài)可控的充要條件是M矩陣滿秩，即 M2=C* A* C* A*2 C* A*N-1C* (4.

48、57) 系統完全狀態(tài)可觀性的充要條件是H*矩陣滿秩，即 H*2=B AB A2BAN-1B (4.58) 綜上所述,給定系統的可觀性可以通過其對偶系統的可控性來檢驗,而研究系統的可控性則可以通過其對偶系統的可觀性來研究,這些性質稱為系統的對偶原則。4.3 控制系統的時域仿真 前一節(jié)以經典控制器和狀態(tài)反饋控制器為例講述了控制系統控制參數的設計過程。設計過程主要依據的是系統的時域特性,即閉環(huán)系統狀態(tài)矩陣的特征值決定了閉環(huán)系統的暫態(tài)響應特性。控制器設計的目標是選擇控制器的增益,使得閉環(huán)系統狀態(tài)矩陣的特征值位于期望的極點位置。 控制器參數確定以后,下一步需要對閉環(huán)系統進行仿真。在設計過程中研究的對象一

49、般是系統的線性化模型。而在仿真過程中,應該盡可能準確地再現實際的系統模型,這就常常要求以實際的時變或非線性系統為仿真的對象。這一節(jié)主要論述控制器參數確定以后的系統仿真過程。為簡單起見,仍然以SISO系統為研究對象。其線性與非線性模型分別為 (4.59) (4.60) 4.3.1 經典比例控制器設計 單位反饋回路的簡單比例控制系統框圖如圖4.19所示。系統控制輸入為(4.61) 使用線性化模型的閉環(huán)仿真方程為(4.62) 寫成標準的狀態(tài)方程形式 (4.63) 如果采用非線性模型進行仿真,必須使用MATLAB中標準的ODE求解器來完成。為此,用戶必須定義一個函數文件,MATLAB中的ODE求解器調

50、用該函數來完成非線性系統的仿真。該函數文件包含下面的內容： （1）指定t時刻的參考輸入rd(t)。 （2）計算系統t時刻的輸出y(t)=CTX(t)。 （3）確定t時刻的輸入u(t)=Kc(rd(t)-y(t)。 （4）計算t時刻狀態(tài)變量的導數 圖4.19 SISO系統的經典比例單位反饋控制框圖 4.3.2 狀態(tài)反饋控制器 帶全狀態(tài)觀測器的狀態(tài)反饋控制框圖如圖4.20所示。該框圖類似于圖4.16,不同之處在于這里的模型包括一個附加的穩(wěn)態(tài)狀態(tài)增益模塊,其中包含一個歸一化的增益變量Nr。 我們再來看看該系統線性化模型與非線性模型的仿真方程,并且重新調整控制器的增益（KTs,L和Nr等）。控制規(guī)律寫

51、成(4.64) 該系統的閉環(huán)模型包括實際的系統狀態(tài)X(t)以及估計的系統狀態(tài) 。這樣,該系統具有2N個未知參數。對于線性系統模型,該系統的完整模型為被控系統模型(4.65) 觀測器模型 (4.66) 寫成標準的狀態(tài)方程形式為 圖4.20 具有全狀態(tài)觀測器的狀態(tài)反饋控制系統 4.3.3 混合控制器 下面將研究一種將狀態(tài)反饋控制運用到經典控制器中的例子。在這個例子中,狀態(tài)反饋制器為經典控制器提供修正的參考信號。例如,對于圖4.18的經典控制系統,如果選擇的增益Kc不滿足系統的設計要求,我們可能會手動調節(jié)參考點rd,以使系統具有更好的暫態(tài)響應特性。可以設想一種控制方案來自動為參考信號rd(t)提供實

52、時的修正信號rdm(t),控制框圖如圖4.20所示。 圖中左邊的sum模塊的輸出為修正后的參考點信號,在此之后,該系統看上去就像是具有單位反饋回路和比例增益Kc的經典控制器。在這個控制方案中,狀態(tài)觀測器和狀態(tài)反饋回路的目標是提供一個修正的參考信號,來改善經典控制系統的控制效果。這種方法稱為狀態(tài)反饋輔助的經典控制方法（SFACC）,也稱為混合的控制器。 圖4.21顯示的控制系統實際上可以看作是前面講述的狀態(tài)反饋控制系統,其控制框圖可以認為是從圖4.20的標準形式通過某些模塊變換發(fā)展而來。圖4.224.24顯示了其中的變換,主要步驟包括: （1）在系統輸入前加上一個比例增益模塊,并且在相應的引入回

53、路中抵消它的影響（如圖4.22所示）。圖4.21 具有嵌入經典控制器的狀態(tài)反饋控制框圖（SFACC） （2）對系統輸出y構造一個負反饋,它可以對估計輸出 的正反饋加以抵消（如圖4.23所示）。 （3）分離出經典控制器部分（圖4.24中的陰影部分）。 （4）最后定義修正的狀態(tài)反饋增益Ksm和修正的穩(wěn)態(tài)狀態(tài)歸一因子Nnn。(4.68) (4.69) 圖4.22 在SFC中加入經典控制增益 圖4.23 加入經典負反饋回路 圖4.24 分離出經典比例控制器部分 圖4.21所示的SFACC控制系統與傳統的狀態(tài)反饋控制系統SFC（如圖4.20所示）在功能上是一致的。僅有的近似存在于第(2)步,這一步假設估

54、計的輸出與實際的系統輸出完全抵消。 如果狀態(tài)觀測器設計成相對控制系統有足夠快的反應時間,這種近似是可以滿足要求的（對于線性模型）。 從以上的討論中可以看出,SFC和SFACC具有相同的暫態(tài)響應。然而,SFACC算法比前者具有更好的控制效果。并且SFACC可以方便地在原先經典控制系統的基礎上進行擴展,而不需要破壞原來已經存在的控制系統。這讓用戶可以從經典控制器直接過渡到現代控制系統。 為了對SFACC控制系統進行仿真,我們假定控制規(guī)律為修正后的參考信號為 將式（4.71）代入式（4.70），有 （4.70） (4.71) (4.72) 與SAC控制系統類似,SFACC閉環(huán)系統既包含實際狀態(tài)也包含

55、估計狀態(tài),因此,可以獲得2N個未知參數。對于線性系統模型,有 系統模型為(4.73) (4.74) (4.75) 或者 (4.76) 寫成標準的狀態(tài)方程形式為 其中 (4.78) 4.3.4 參考輸入增益的調整 在前面的框圖中,為了減小比例控制器帶來的穩(wěn)態(tài)誤差,引入了一個輸入歸一化模塊。該模塊的增益值可以根據下面的算式進行確定（圖4.25顯示了系統的簡化模型）。 線性模型(4.79) (4.80) 控制律 系統穩(wěn)態(tài)時,我們希望yss與參考點rdss相同。簡化系統穩(wěn)態(tài)時的狀態(tài)方程為(4.81) 求解該方程,得穩(wěn)態(tài)時的狀態(tài)向量為(4.82) 圖4.25 具有SS增益模塊的狀態(tài)反饋控制系統 穩(wěn)態(tài)輸出

56、為(4.83) 最后將yss=rdss代入,得到Nr的計算公式為(4.84) 4.4 實例:倒擺系統的建模與仿真 4.4.1 引言 這一節(jié)我們先看一個實際的例子。圖4.26是某個倒擺系統的示意圖,倒擺通過轉動關節(jié)安裝在驅動小車上,通過對小車施加一定的驅動,使倒擺保持一定的姿態(tài)。這是姿態(tài)控制問題的典型代表。圖4.26 驅動小車上的倒擺示意圖 4.4.2 基本方程 假定倒擺由無質量的輕桿和質量為m的小球組成,小車的質量為M,系統所受的外力包括小球受到的重力和對小車水平方向的驅動力u(t)。x(t)和(t)分別表示小車的水平坐標和倒擺偏離垂直方向的角度,如圖4.27所示。 根據牛頓運動學第二定律,在

57、水平x軸方向上系統滿足下面的方程：(4.85) 圖4.27 小球受到的力矢量圖 其中,小球受到的重力為常值,其中心的坐標表示為(xG,yG),滿足 xG=x+lsin yG=lcos (4.86) l為倒擺的桿長。將式（4.86）代入式（4.85）,得到(4.87) 由下面的基本關系(4.88) (4.89) 得到 同樣,通過小球的力矩平衡關系,可以得到小球的平衡方程：(4.90) (4.91) 由式（4.86）（4.89）,可以計算得到Fx和Fy的表達式:(4.92)(4.93) 將式（4.92）、（4.93）代入式（4.91）中,兩邊消去l：最后得到 (4.94) 方程（4.90）和（4.

58、94）構成了倒擺系統的動力學模型。顯然,從數學的角度上看,該系統為明顯的非線性系統。但是對小車施加驅動力的目的是保持倒擺在垂直方向上的姿態(tài),因此,我們關注的是倒擺在垂直方向附近的動態(tài)變化,為此,可以將系統在該參考位置（=0）進行線性化,這樣可以對簡化得到的線性模型進行控制器的設計。在后面的內容中,將分別對線性化前后的模型進行比較和研究。 4.4.3 非線性狀態(tài)方程描述 為了對倒擺系統的非線性模型進行數值仿真,需要首先將模型表示成標準的狀態(tài)方程形式：由式（4.94）可以得到 (4.95)(4.96) 將其代入式（4.90）,有 (4.97) (4.98) 將其代入式（4.90）中 同樣,可以由式

59、（4.94）得到(4.99) 上式滿足式（4.95）的標準狀態(tài)方程的形式。假定系統輸出為倒擺的角度和小車的x軸坐標,則系統的輸出方程為(4.102) 4.4.4 線性狀態(tài)方程描述 為了得到倒擺在垂直位置附近的線性化模型,我們只要對式（4.97）描述的非線性狀態(tài)方程進行線性化即可。根據前幾節(jié)論述的方法,該系統的線性化模型為(4.103) 這里的參考狀態(tài)為倒擺靜止在垂直方向,并且小車的驅動力為零,即Z0=0,u0=0。下面來逐項計算Jacobian矩陣。 Jz(Z0,u0)的第一列由 決定,其中第一和第三個函數與z1不直接相關,其偏導數等于零。對于第二項,有 同樣,對第四項有 Jz(Z0,u0)的

60、第二列由fi/z2|z0,u0決定,式(4.101)顯示f1/z2=1,f3/z2=0,對于第二和第四項有Jz(Z0,u0)的第三和第四列中僅有的非零項為 4.4.5 倒擺系統的MATLAB仿真 1.倒擺系統的開環(huán)與閉環(huán)仿真 這一部分的仿真主要研究下面三方面的內容： (1)對倒擺非線性模型與線性化模型不穩(wěn)定動態(tài)特性的比較。 (2)線性化模型可以由MATLAB函數linmod直接計算得到，以及如何用S函數表示非線性模型。 (3)列舉一些不夠理想的經典反饋控制方案。 仿真計算的M文件包括INVPN1.M、INVPNNL1.M和INVPNNL2.M等等。其中,INVPN1.M為主文件,主要功能是控制