《控制系統(tǒng)仿真及MATLAB語言-連續(xù)系統(tǒng)的離散化方法》由會員分享，可在線閱讀，更多相關(guān)《控制系統(tǒng)仿真及MATLAB語言-連續(xù)系統(tǒng)的離散化方法（44頁珍藏版）》請在裝配圖網(wǎng)上搜索。

1、第四章 連續(xù)系統(tǒng)的離散化方法 4.1 常微分方程的數(shù)值解法 00 ( , ) () dx f t x dt x t x 一 . 數(shù)值求解的基本概念 設(shè)微分方程為 則求解方程中函數(shù) x(t)問題的常微分方程初值問題 所謂數(shù)值求解就是要在時間區(qū)間 a, b中取若干離散點 求出微分方程在這些時刻的近似值 0 1 2 Nx x x x 01 Na t t t b 2 ( 2 ) ( ) 0 0 0 0 0( ) ( ) ( ) ( ) ( )2 ! ! k khhx t h x t h x t x t x t k nnn tth 1 取前兩項近似： 1 ( , )k k k kx h f t x x

2、這種方法的幾何意義就是把 f(t,x)在 區(qū)間 tk,tk+1內(nèi)的曲邊面積用矩形面 積近似代替。計算簡單，計算量小， 而且可以自啟動。當 h很小時，造成 的誤差是允許的。該算法具有一階 精度。 取 k=0,1,2,N, 從 t0開始，逐點遞推求解 t1時的 y1， t2時的 y2, 直至 tn時的 yn，稱之為歐拉遞推公式。 矩形面積 1. 歐拉法 歐拉法的特點： 導(dǎo)出簡單，幾何意 義明顯，便于理解，能說明構(gòu)造數(shù) 值解法一般計算公式的基本思想 。 通常用它來說明有關(guān)的基本概念。 例 設(shè)系統(tǒng)方程為 用 Euler法求其數(shù)值解（取步長 , ) 1.0h 2 , 0 1x x x 10 t 遞推公式

3、為 1 , 1 0 . 1n n n n n nx x h f t x x x 則 4628.01.01,0.1 7519.0819.01.01819.01.01,3.0 819.091.09.01.01,2.0 9.01.01,1.0 1,0 991010 2233 1122 0011 00 yyyt yyyt yyyt yyyt yt 已知方程的解析解為 精確解和解析解作比較： 誤差在 數(shù)量級， 精度較差。 210 ty 1 1 t 0 0.1 0.2 0.3 0.4 0.5 1.0 精確解 1 0.909 0.833 0.769 0.666 0.625 0.5 數(shù)值解 1 0.9 0.8

4、19 0.752 0.659 0.647 0.463 2. 龍格庫塔法 基本思想： 取 Taylor級數(shù)展開式前三項 近似求解，并利用線性組合代替導(dǎo)數(shù)的 求解。 既可避免計算高階導(dǎo)數(shù)，又可提高 數(shù)值積分的精度，這就是 Runge-Kutta 法的基本思想。 2. 龍格庫塔法 kktxff r為精度階次， ai為待定系數(shù)，由精度確定； ki用下 式表示 線性組合 1 12 1 ( , ) , 2 , 3 , i i k k j j k f t b h x h b k i r 2 1 ()2! kkk k k t x k hx x h f f f f 1 1 r k k i i i x x h a

5、 k 2 ( 2 ) 0 0 0 0( ) ( ) ( ) ( )2 hx t h x t h x t x t 等各階導(dǎo)數(shù)不易計算，用下式中 ki的線性組合代替 1)當 r=1時： 1 1 1 ,kkx x h a k 1 ( , ) ,kkk f t x 與 Taloy展開式相比較，可得 a1=1,則上式成為 11 ( , ) ,k k k k kx x h k x h f t x 歐拉遞推公式 2)當 r=2時： 1 2 1 2 1 kk kk k f t , x k f t b h , x h b k 2 1 2 1kkk t xk f b h f b h k f 1 2 1kkf t

6、b h , x h b k kkt ,x將 在點 展成 Taylor級數(shù) 與臺勞公式的二階展開近似公式相比，可得以下關(guān)系： 12 22 21 1 12 12 aa ba ab 三個方程，四個未知數(shù)，解不唯一 各個系數(shù)的幾種取法 見書上。 3) r=4時，四階龍格庫塔公式 -最常用： 仿真中遇到的大多數(shù)工程實際問題，四階龍格庫塔法以能 滿足精度要求，其截斷誤差 o(h5) 與 h5同數(shù)量級。該法可以 自啟動。 1 1 2 3 4 1 21 32 43 22 6 22 22 kk kk kk kk kk h x x ( K K K K ) K f t , x hh K f t , x K hh K

7、 f t , x K K f t h , x hK 4)、狀態(tài)空間四階龍格 -庫塔遞推式 若單輸入單輸出系統(tǒng)的狀態(tài)空間表達式為： X A X B U Y C X D U 在仿真中，對于 n階系統(tǒng)，狀態(tài)方程可以寫成一階微分方程 1 1 2 2 1 2 3( , , , , , ) ( 1 , 2 , , ) i i i in n i in x a x a x a x b u f t x x x x i n 根據(jù)四階龍格 -庫塔公式 ,有 1 1 2 3 4 1 1 1 2 2 2 1 1 1 1 3 1 1 2 2 4 1 1 3 3 ( 2 2 ) ( 1 , 2 , , ) 6 () ( )

8、 ( ) ( ) 2 2 2 ( ) ( ) ( ) 2 2 2 ( ) ( kk i i i i i i k k k i i i in n i k kk i i i in n i i k kk i i i in n i i k kk i i i in n i h x x k k k k i n k a x a x a x b u t h h h k a x k a x k b u t h h h k a x k a x k b u t k a x hk a x hk ) ( ) ik b u t h T=tk時刻的 xi值 T=tk+h時刻的 xi值 另 1 1 1 1 2 1 2 1 1

9、2 1 1 2 1 2 2 2 2 3 1 3 2 3 3 4 1 4 2 4 4 , , , k k k T k k k T nn TT nn x x x x x x k k k k k k k k k k k k k + 1 k 1 xx KK 狀態(tài)方程的四階龍格 -庫塔公式如下： 2 3 4 2 32 43 ( 2 2 ) 6 () ( ) ( ) 22 ( ) ( ) 22 ( ) ( ) k k k k h ut hh ut hh ut h u t h k + 1 k 1 1k k1 k k k + 1 k + 1 x x K K K K K A x B K A x K B K A

10、x K B K A x K B y C x RK法的特點： 1 需要存儲的數(shù)據(jù)少，占用的存儲空間少； 2 只需知道初值，即可啟動遞推公式進行計 算，可自啟動； 3 容易實現(xiàn)變步長運算。 4 每積分一步需要計算多次右函數(shù)，計算量 大。 基于龍格庫塔法， MATLAB提供了求常微分方程數(shù)值解 的函數(shù)，一般調(diào)用格式為： t, x=ode23(xfun, t0, tf ,x0) t, x=ode45(xfun, t0, tf ,x0) 常微分方 程函數(shù)名 起始 時間 終止 時間 初始狀 態(tài)向量 輸入 輸出 4/5階龍格 -庫塔算法 2/3階龍格 -庫塔算法 3.常微分方程 Matlab求解 例 1 經(jīng)

11、典的微分方程 2 2 2 2 1 0 01 01 d x dx ( x ) x dt dt x ( ) ; x ( ) 解 : 令 y1=x， y2=x 則微分方程變?yōu)橐浑A微分方程組： 12 2 2 1 2 1 12 21 0 2 0 0 y y y ( y ) y y y ( ) , y ( ) 1、建立 M-文件 vdp.m如下： function dy=vdp(t,y) dy=zeros(2,1); dy(1)=y(2); dy(2)=2*(1-y(1)2)*y(2)-y(1); 2、取 t0=0， tf=20，輸入命令： T,Y=ode45(vdp,0 10,1;1); plot(T,

12、Y(:,1),-, T,Y(:,2) 3、結(jié)果 練習(xí) 解微分方程組 . 1)0(,1)0(,0)0( 51.0 321 213 312 321 yyy yyy yyy yyy 解 1、建立 m-文件 rigid.m如下： function dy=rigid(t,y) dy=zeros(3,1); dy(1)=y(2)*y(3); dy(2)=-y(1)*y(3); dy(3)=-0.51*y(1)*y(2); 2、 取 t0=0， tf=12，輸入命令： T,Y=ode45(rigid,0 12,0 1 1); plot(T,Y(:,1),-,T,Y(:,2),*,T,Y(:,3),+) 3、

13、結(jié)果如圖 圖中， y1的圖形為實線， y2的圖形為“ *”線， y3的圖形為“ +”線 . 4.2 數(shù)值算法的穩(wěn)定性及求解原則 1.數(shù)值算法的穩(wěn)定性 特征根在 s平面的左半平面，系統(tǒng)穩(wěn)定。 (1)歐拉法： 穩(wěn)定： (2)梯形法： 恒穩(wěn) 11zh ( ) ( 1 ) ( )zX z h X z 1z ( ) ( 1 ) ( 1 ) ( ) 22 hhz X z X z 2.數(shù)值算法的選擇原則 Matlab提供了微分方程數(shù)值求解的一般方法 ,作為仿真算 法的使用者 ,可不必考慮算法具體實現(xiàn) ,而應(yīng)關(guān)心各種方法在 使用中會出現(xiàn)的問題 ,以及如何在仿真中恰當?shù)倪x用這些方法 . 一般 ,選用數(shù)值算法從以

14、下幾個方面考慮 : (1)精度 受算法和 h影響 截斷誤差 +舍入誤差 =累計誤差 (2)計算速度 受算法和 h影響 算法簡單，速度就快些。 (3)穩(wěn)定性 受 h影響，一般 h （ 2-3） 系統(tǒng)最小時間 4.3 數(shù)值算法中的 “ 病態(tài) ” 問題 1 “病態(tài)”常微分方程 例： ( 0) (1 , 0 , 1 ) T X AX X 其中 1 2 3( , , )TTX x x x - 2 1 1 9 2 0 A = 1 9 2 1 2 0 4 0 4 0 4 0 采用四階龍格庫塔法 h=0.01時，計算時間長 h=0.04時，誤差很大 當 h0.05后，曲線發(fā)散振蕩，數(shù)值不穩(wěn)定，完全失去意義 系

15、統(tǒng)矩陣的特征值差異較大 一般線性常微分方程組： 00( ) ( ) ( ) , ( )X t A X t B U t X t X 的系數(shù)矩陣 A的特征值具有如下特征： 11 Re ( ) 0 m a x Re ( ) m in Re ( ) i iiinin 則稱為“病態(tài)”方程。 定義： 2 控制系統(tǒng)仿真中的“病態(tài)”問題 1 病態(tài)系統(tǒng)中絕對值最大的特征值對應(yīng)于系統(tǒng)動態(tài)性能 解中瞬態(tài)分量衰減最快的部分，它反映了系統(tǒng)的動態(tài)響 應(yīng)和系統(tǒng)的反應(yīng)靈敏度。一般與系統(tǒng)中具有最小時間常 數(shù) Tmin的環(huán)節(jié)有關(guān)，要求計算步長 h取得很小。 2 病態(tài)系統(tǒng)中絕對值最小的特征值對應(yīng)于系統(tǒng)動態(tài)性能 解中瞬態(tài)分量衰減最慢

16、的部分，它決定了整個系統(tǒng)的動 態(tài)過渡過程時間的長短。一般與系統(tǒng)中具有最小時間常 數(shù) Tmax的環(huán)節(jié)有關(guān)，要求計算步長 h取得很大。 3 對于病態(tài)問題的仿真需要尋求更加合理的算法，以解 決病態(tài)系統(tǒng)帶來的選取計算步長與計算精度，計算時間 之間的矛盾。 3 “病態(tài)”系統(tǒng)的仿真方法 采用穩(wěn)定性好，計算精度高的數(shù)值算法，并且允許計 算步長能根據(jù)系統(tǒng)性能動態(tài)變化的情況在一定范圍內(nèi)作相 應(yīng)的變化，采用 隱式吉爾法 1 0 1 1 1 1k k k r k r kx a x a x a x h f 該法已經(jīng)證明對病態(tài)方程求解過程是數(shù)值穩(wěn)定的。 隱式吉爾法從理論上十分適應(yīng)于病態(tài)系統(tǒng) ，但需要解決好 以下問題 (

17、1) 自啟動 r階多步算式無法自啟動，需要用單步法求出前 r步值 (2) 預(yù)估迭代 迭代方法要求收斂性良好，否則在大步長時會造成數(shù) 值發(fā)散。 (3) 變步長 初始階段采用小步長，隨后可逐步放大步長。 對不同精度要求的系統(tǒng)仿真，要考慮變階次問題， 即為減小每一步計算的截斷誤差，以提高精度，應(yīng)選用 較高的階次，而當精度較低時，為減少工作量，則應(yīng)選 取較低的階次。仿真時應(yīng)根據(jù)估計誤差 與給定的誤差精 度相比較改變步長或階次來重新計算。 4.4 連續(xù)系統(tǒng)狀態(tài)方程的離散化 上章所述的連續(xù)系統(tǒng)數(shù)學(xué)模型的離散化， 是通過數(shù)值積分法實現(xiàn)的，盡管面向結(jié) 構(gòu)圖的仿真方法是按環(huán)節(jié)給定參數(shù)，但 是在計算時還是 按整個

18、系統(tǒng)進行離散化 ， 這就不便于引進非線性環(huán)節(jié)以進行非線 性系統(tǒng)的仿真。在本節(jié)，將介紹連續(xù)系 統(tǒng)離散模型的建立和仿真。 開 始 輸入開環(huán)傳函分母、分子系數(shù) ji ca , 輸入反饋函數(shù) V 、狀態(tài)初始值 X 0 求開環(huán)狀態(tài)方程系 數(shù)陣 A , B , C 求閉環(huán)狀態(tài)方程系數(shù)陣 B v CAA b 輸入初始時間 0 t 、終止時間 f t 、 計算步長 h 、輸入幅值 r 求龍格 庫塔法各次斜率 4,3,2,1, jK j 求 11 , kk yX 輸出數(shù)據(jù)、曲線 1 kk XX f t 到否 ？ N Y 結(jié) 束 數(shù)值積分法疊代求解 h改變時，疊代過程重復(fù)求解，費時繁瑣 不能對非線性環(huán)節(jié)單獨考慮。

19、 連續(xù)系統(tǒng)離散化 思想：用差分方程描述連續(xù)系統(tǒng)的狀態(tài)方程模型 （因為差分方程的主要特點就是方程中各變量由各相鄰時刻的變 化量制約，這相當于遞推方程） 1、連續(xù)系統(tǒng)的離散化 DuCXy BuAXX 設(shè)連續(xù)系統(tǒng)狀態(tài)方程為 00 )( XtX 其中 為狀態(tài)初始值 . 則由現(xiàn)代控制理論基礎(chǔ)知，狀態(tài)變量 X(t)的解為 () 0 0 ( ) ( ) kT A k T A k Tk T e e d X X B u t tAtA dee 0 )( 0 )()( uBXtX 而 t = ( k + 1 )T 時，可表示為 ( ) 0 ( ( 1 ) ) ( ) ( ) T TTk T e k T e d k

20、T AAX X B u 當系統(tǒng)輸入 u(t)給定時，可求出系統(tǒng)離散化狀態(tài)方程 的解。一般， u(t)未知，通常采用兩種方法近似處 理： (1)令 u(kT+t) u(kT) (0 t T) 相當在系統(tǒng)輸入端加一個采樣開關(guān)和零階保持器 X(k+1)T) = GX(kT) + Hu(kT) G = e A T，為 t = T 時的狀態(tài)轉(zhuǎn)移矩陣 () 0 T TH e d A B 離散后的狀態(tài)空間表達式為： ( 1 ) ( ) ( ) ( ) ( ) ( ) x k T G x k T H u k T y k T C x k T D u k T 、 Matlab表示 DuCXy BuAXX 已知連續(xù)

21、系統(tǒng)狀態(tài)方程為 在采樣周期 T下離散后的狀態(tài)空間表達可表示為： ( 1 ) ( ) ( ) ( ) ( ) ( ) x k T G x k T H u k T y k T C x k T D u k T () 0 G, T TTe H e d AA B其 中 ， 在 Matlab中，若已知連續(xù)系統(tǒng)狀態(tài)方程各陣模型參數(shù) (A、 B、 C、 D) 以及采樣周期 T，則語句： G， H = c2d (A,B,T) 返回的矩陣 G 、 H 就是所要求的 ( T ) 、 m ( T ) 。 此外， Matlab還提供了功能更強的求取連續(xù)系統(tǒng)離 散化矩陣函數(shù) c2dm(),他容許調(diào)用時選用離散化變換 方式

22、，并且得到的是標準的離散化狀態(tài)方程。 G,H,C,D=c2dm (A,B,C,D,T, 選項 ) 表 離散化變換方式選項 選 項 說 明 Zoh 假設(shè)輸入端加一個采樣開關(guān)和零階保持器 Foh 假設(shè)輸入端加一個采樣開關(guān)和一階保持器 。 Tustin 采用雙線性變換 （ Tustin算法 ） 方法 Prewarp 采用改進的 Tustin變換方法 Matched 采用 SISO系統(tǒng)的零極點匹配法 2.離散函數(shù)的連續(xù)化 在 MATLAB中也提供了從離散化系統(tǒng)轉(zhuǎn)換為 連續(xù)系統(tǒng)各系數(shù)矩陣求取的功能函數(shù) ， 其調(diào) 用格式分別如下 A ,B=d2c(G ,H ,T) 或 A ,B ,C ,D=d2cm (G

23、 ,H ,C,D ,T ,選項 ) 其中選項同上。 例 對連續(xù)系統(tǒng)在采樣周期 T=0.1時進行離散化。 63 1 2 5 ( s )G ( s ) ( s ) ( s ) ( s ) 解：利用以下 Matlab命令，對系統(tǒng)進行離散化 k=6; Z=-3; P=-1； -2； -5； T=0.1; A,B,C,D=zp2ss(Z,P,k); G1,H1=c2d(A,B,T); G2,H2,C2,D2=c2dm(A,B,C,D,T,zoh); G3,H3,C3,D3=c2dm(A,B,C,D,T,foh); 3、脈沖傳遞函數(shù)求解 連續(xù) 離散 function nd,dd=tfc2d(nc,dc,T) A,B,C,D=tf2ss(nc,dc); G,H=c2d(A,B,T); nd,dd=ss2tf(G,H,C,D); 離散 連續(xù) function nc,dc=tfd2c(nd,dd,T) A,B,C,D=tf2ss(nd,dd); G,H=d2c(A,B,T); nc,dc=ss2tf(G,H,C,D);