最新控制系統(tǒng)仿真及MATLAB語言-連續(xù)系統(tǒng)的離散化方法精品ppt課件講課稿

1、控制系統(tǒng)仿真及MATLAB語言-連續(xù)系統(tǒng)的離散化方法4.1 常微分方程的數(shù)值解法例例 設(shè)系統(tǒng)方程為 用Euler法求其數(shù)值解（取步長 , )1 . 0h 2,01xxx 10 t遞推公式為1,10.1nnnnnnxxhftxxx則4628. 01 . 01, 0 . 17519. 0819. 01 . 01819. 01 . 01, 3 . 0819. 091. 09 . 01 . 01, 2 . 09 . 01 . 01, 1 . 01, 099101022331122001100 yyytyyytyyytyyytyt已知方程的解析解為 精確解和解析解作比較： 誤差在 數(shù)量級，精度較差。 2

2、10ty11t00.10.20.30.40.51.0精確解10.909 0.833 0.769 0.666 0.625 0.5數(shù)值解10.90.819 0.752 0.659 0.6470.4632. 龍格庫塔法龍格庫塔法 基本思想：基本思想：取Taylor級數(shù)展開式前三項近似求解，并利用線性組合代替導(dǎo)數(shù)的求解。 既可避免計算高階導(dǎo)數(shù)，又可提高數(shù)值積分的精度，這就是Runge-Kutta法的基本思想。2. 龍格庫塔法龍格庫塔法kktxffr為精度階次，為精度階次，ai為待定系數(shù)，由精度確定；為待定系數(shù)，由精度確定；ki用下用下式表示式表示線性組合線性組合1121(,) ,2,3,iikkjjk

3、f tb h xhbkir21()2!kkkkktxkhxxh ffff11rkkiiixxha k2(2)0000()( )( )( )2hx thx th x txt等各階導(dǎo)數(shù)不易計算，用下式中等各階導(dǎo)數(shù)不易計算，用下式中kiki的線性組合代替的線性組合代替1)當(dāng)r=1時：11 1,kkxxha k1( ,),kkkf tx與Taloy展開式相比較，可得a1=1,則上式成為11( ,),kkkkkxxhkxhf tx歐拉遞推公式2)當(dāng)r=2時：12121kkkkkft ,xkftb h,xhb k2121kkktxkfbhfb hk f12 1kkf tbh,xhbkkkt ,x將 在點

4、展成Taylor級數(shù)與臺勞公式的二階展開近似公式相比，可得以下關(guān)系：與臺勞公式的二階展開近似公式相比，可得以下關(guān)系：12222 111 21 2aab aa b三個方程，四個未知數(shù)，解不唯一三個方程，四個未知數(shù)，解不唯一各個系數(shù)的幾種取法各個系數(shù)的幾種取法見書上。見書上。3) r=4時，四階龍格庫塔公式時，四階龍格庫塔公式-最常用：最常用：仿真中遇到的大多數(shù)工程實際問題，四階龍格庫塔法以能滿足精度要求，其截斷誤差o(h5) 與h5同數(shù)量級。該法可以自啟動。1123412132432262222kkkkkkkkkkhxx(KKKK )Kf t ,xhhKft,xKhhKft,xKKf th,xh

5、K 4)、狀態(tài)空間四階龍格、狀態(tài)空間四階龍格-庫塔遞推式庫塔遞推式 若單輸入單輸出系統(tǒng)的狀態(tài)空間表達(dá)式為：若單輸入單輸出系統(tǒng)的狀態(tài)空間表達(dá)式為：XAXBUYCXDU 在仿真中，對于在仿真中，對于n階系統(tǒng)，狀態(tài)方程可以寫成一階微分方程階系統(tǒng)，狀態(tài)方程可以寫成一階微分方程1 122123( ,)(1,2, )iiiinniinxa xa xa xbuf t x x xxin 根據(jù)四階龍格根據(jù)四階龍格-庫塔公式庫塔公式,有有1123411 122211113112241133(22)(1,2, )6( )()()()222()()()222()(kkiiiiiikkkiiiinnikkkiiiinn

6、iikkkiiiinniikkkiiiinnihxxkkkkinka xa xa xbu thhhkaxkaxkbu thhhkaxkaxkbu tkaxhkaxhk)()ikbu thT=tT=tk k時刻的時刻的x xi i值值T=tT=tk k+h+h時刻的時刻的x xi i值值 另另111121211211212222313233414244 , , ,kkkTkkk TnnTTnnTTnnxxxxxxkkkkkkkkkkkkk+1k1xxKKKK 狀態(tài)方程的四階龍格狀態(tài)方程的四階龍格-庫塔公式如下：庫塔公式如下：23423243(22)6( )()()22()()22()()kkkk

7、hu thhu thhu thu thk+1k11kk1kkk+1k+1xxKKKKKAxBKA xKBKA xKBKA xKByCxRK法的特點：法的特點： 1 需要存儲的數(shù)據(jù)少，占用的存儲空間少；2 只需知道初值，即可啟動遞推公式進(jìn)行計 算，可自啟動；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ù)名程函數(shù)名起始起始時間時間終止終止時間時間初始狀初始狀態(tài)向量態(tài)向

8、量輸入輸入輸出輸出4/54/5階龍格階龍格- -庫塔算法庫塔算法2/32/3階龍格階龍格- -庫塔算法庫塔算法3.3.常微分方程常微分方程MatlabMatlab求解求解解解: 令令 y1=x，y2=x1、建立、建立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,Y(:,1),-, T,YT,Y(:,2)3、結(jié)果、結(jié)果1)0(, 1)0

9、(, 0)0(51. 0321213312321yyyyyyyyyyyy解解 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、結(jié)果如圖、結(jié)果如圖圖中，y1的圖形為實線，y2的圖形為“*”線，y3的圖形為“+”線.

10、4.2 數(shù)值算法的穩(wěn)定性及求解原則1.1.數(shù)值算法的穩(wěn)定性數(shù)值算法的穩(wěn)定性 特征根在特征根在s s平面的左半平面，系統(tǒng)穩(wěn)定。平面的左半平面，系統(tǒng)穩(wěn)定。(1)歐拉法： 穩(wěn)定：(2)梯形法： 恒穩(wěn)11zh( )(1)( )zX zhX z1z ( )(1)(1)( )22hhzX zX z2.2.數(shù)值算法的選擇原則數(shù)值算法的選擇原則 Matlab提供了微分方程數(shù)值求解的一般方法,作為仿真算法的使用者,可不必考慮算法具體實現(xiàn),而應(yīng)關(guān)心各種方法在使用中會出現(xiàn)的問題,以及如何在仿真中恰當(dāng)?shù)倪x用這些方法. 一般,選用數(shù)值算法從以下幾個方面考慮:(1)精度受算法和h影響截斷誤差+舍入誤差=累計誤差(2)計算

11、速度受算法和h影響算法簡單，速度就快些。(3)穩(wěn)定性 受h影響，一般 h（2-3）系統(tǒng)最小時間系統(tǒng)最小時間4.3 數(shù)值算法中的數(shù)值算法中的“病態(tài)病態(tài)”問題問題 1 “病態(tài)”常微分方程例：(0)(1,0, 1)TXAXX其中123( ,)TTXx x x-211920A= 192120404040采用四階龍格庫塔法h=0.01時，計算時間長時，計算時間長h=0.04時，誤差很大時，誤差很大當(dāng)當(dāng)h0.05后，曲線發(fā)散振蕩，數(shù)值不穩(wěn)定，完全失去意義后，曲線發(fā)散振蕩，數(shù)值不穩(wěn)定，完全失去意義一般線性常微分方程組：一般線性常微分方程組：00( )( )( ),( )X tAX tBU tX tX的系數(shù)矩

12、陣的系數(shù)矩陣A的特征值具有如下特征：的特征值具有如下特征：11Re()0max Re()min Re()iiii ni n 則稱為則稱為“病態(tài)病態(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)分量衰減最慢的部分，它決定了整個系統(tǒng)的動態(tài)過渡過程時間的長短。一般與系統(tǒng)中具有最小時間常數(shù)Tmax的環(huán)節(jié)有關(guān)，要求計算步長h取得很大。3 對于病態(tài)

13、問題的仿真需要尋求更加合理的算法，以解決病態(tài)系統(tǒng)帶來的選取計算步長與計算精度，計算時間之間的矛盾。3 “病態(tài)”系統(tǒng)的仿真方法 采用穩(wěn)定性好，計算精度高的數(shù)值算法，并且允許計算步長能根據(jù)系統(tǒng)性能動態(tài)變化的情況在一定范圍內(nèi)作相應(yīng)的變化，采用隱式吉爾法101111kkkrk rkxa xa xa xhf該法已經(jīng)證明對病態(tài)方程求解過程是數(shù)值穩(wěn)定的。隱式吉爾法從理論上十分適應(yīng)于病態(tài)系統(tǒng) ，但需要解決好以下問題(1) 自啟動 r階多步算式無法自啟動，需要用單步法求出前r步值(2) 預(yù)估迭代 迭代方法要求收斂性良好，否則在大步長時會造成數(shù) 值發(fā)散。(3) 變步長 初始階段采用小步長，隨后可逐步放大步長。 對

14、不同精度要求的系統(tǒng)仿真，要考慮變階次問題，即為減小每一步計算的截斷誤差，以提高精度，應(yīng)選用較高的階次，而當(dāng)精度較低時，為減少工作量，則應(yīng)選取較低的階次。仿真時應(yīng)根據(jù)估計誤差 與給定的誤差精度相比較改變步長或階次來重新計算。 4.4 連續(xù)系統(tǒng)狀態(tài)方程的離散化上章所述的連續(xù)系統(tǒng)數(shù)學(xué)模型的離散化，上章所述的連續(xù)系統(tǒng)數(shù)學(xué)模型的離散化，是通過數(shù)值積分法實現(xiàn)的，盡管面向結(jié)是通過數(shù)值積分法實現(xiàn)的，盡管面向結(jié)構(gòu)圖的仿真方法是按環(huán)節(jié)給定參數(shù)，但構(gòu)圖的仿真方法是按環(huán)節(jié)給定參數(shù)，但是在計算時還是是在計算時還是按整個系統(tǒng)進(jìn)行離散化按整個系統(tǒng)進(jìn)行離散化，這就不便于引進(jìn)非線性環(huán)節(jié)以進(jìn)行非線這就不便于引進(jìn)非線性環(huán)節(jié)以進(jìn)行非

15、線性系統(tǒng)的仿真。在本節(jié)，將介紹連續(xù)系性系統(tǒng)的仿真。在本節(jié)，將介紹連續(xù)系統(tǒng)離散模型的建立和仿真。統(tǒng)離散模型的建立和仿真。開 始輸入開環(huán)傳函分母、分子系數(shù)jica ,輸入反饋函數(shù) V、狀態(tài)初始值 X0求開環(huán)狀態(tài)方程系數(shù)陣 A,B,C求閉環(huán)狀態(tài)方程系數(shù)陣BvCAAb輸入初始時間0t、終止時間ft、計算步長 h、輸入幅值 r求龍格庫塔法各次斜率4 , 3 , 2 , 1, jKj求11,kkyX輸出數(shù)據(jù)、曲線1kkXXft到否？NY結(jié) 束 h改變時，疊代過程重復(fù)求解，費時繁瑣不能對非線性環(huán)節(jié)單獨考慮。連續(xù)系統(tǒng)離散化連續(xù)系統(tǒng)離散化思想：用差分方程描述連續(xù)系統(tǒng)的狀態(tài)方程模型思想：用差分方程描述連續(xù)系統(tǒng)的狀

16、態(tài)方程模型（因為差分方程的主要特點就是方程中各變量由各相鄰時刻的變（因為差分方程的主要特點就是方程中各變量由各相鄰時刻的變 化量制約，這相當(dāng)于遞推方程）化量制約，這相當(dāng)于遞推方程）1 1、連續(xù)系統(tǒng)的離散化、連續(xù)系統(tǒng)的離散化DuCXyBuAXX 設(shè)連續(xù)系統(tǒng)狀態(tài)方程為設(shè)連續(xù)系統(tǒng)狀態(tài)方程為00)(XtX其中其中為狀態(tài)初始值為狀態(tài)初始值.則由現(xiàn)代控制理論基礎(chǔ)知，狀態(tài)變量則由現(xiàn)代控制理論基礎(chǔ)知，狀態(tài)變量X(t)的解為的解為()00()( )kTAkTA kTkTeedXXButtAtAdee0)(0)()(uBXtX而而 t = ( k + 1 )T 時，可表示為時，可表示為()0(1) )() ()T

17、TTkTekTedkTAAXXBu當(dāng)系統(tǒng)輸入當(dāng)系統(tǒng)輸入u(t)u(t)給定時，可求出系統(tǒng)離散化狀態(tài)方程給定時，可求出系統(tǒng)離散化狀態(tài)方程的解。一般，的解。一般， u(t)u(t)未知，通常采用兩種方法近似處未知，通常采用兩種方法近似處理：理：(1)(1)令令u(kT+t)u(kT+t)u(kT) (0u(kT) (0t tT)T)相當(dāng)在系統(tǒng)輸入端加一個采樣開關(guān)和零階保持器相當(dāng)在系統(tǒng)輸入端加一個采樣開關(guān)和零階保持器 X X(k k+1)+1)T T) = ) = GXGX( (kTkT) + ) + HuHu( (kTkT) )G G = = e e A TA T，為，為t t = = T T 時

18、的狀態(tài)轉(zhuǎn)移矩陣時的狀態(tài)轉(zhuǎn)移矩陣()0TTHedAB離散后的狀態(tài)空間表達(dá)式為：離散后的狀態(tài)空間表達(dá)式為：(1) ()()()()()x kTGx kTHu kTy kTCx kTDu kT、MatlabMatlab表示表示DuCXyBuAXX已知連續(xù)系統(tǒng)狀態(tài)方程為已知連續(xù)系統(tǒng)狀態(tài)方程為在采樣周期在采樣周期T下離散后的狀態(tài)空間表達(dá)可表示為：下離散后的狀態(tài)空間表達(dá)可表示為：(1) ()()()()()x kTGx kTHu kTy kTCx kTDu kT()0G,TTTeHedAAB其中，在在MatlabMatlab中，若已知連續(xù)系統(tǒng)狀態(tài)方程各陣模型參數(shù)中，若已知連續(xù)系統(tǒng)狀態(tài)方程各陣模型參數(shù) (

19、(A A、B B、C C、D D) ) 以及采樣周期以及采樣周期T T，則語句：，則語句： G G，H H = c2d ( = c2d (A A, ,B B, ,T)返回的矩陣返回的矩陣G G 、H H 就是所要求的就是所要求的( ( T T ) ) 、m m ( ( T T ) ) 。此外，此外， MatlabMatlab還提供了功能更強的求取連續(xù)系統(tǒng)離還提供了功能更強的求取連續(xù)系統(tǒng)離散化矩陣函數(shù)散化矩陣函數(shù)c2dm(),c2dm(),他容許調(diào)用時選用離散化變換他容許調(diào)用時選用離散化變換方式，并且得到的是標(biāo)準(zhǔn)的離散化狀態(tài)方程。方式，并且得到的是標(biāo)準(zhǔn)的離散化狀態(tài)方程。 G G, ,H H, ,C C, ,D D=c2dm (=c2dm (A A, ,B B, ,C C, ,D D, ,T T,選項選項) ) 表表 離散化變換方式選項離散化變換方式選項選選 項項說說 明明Zoh假設(shè)輸入端加一個采樣開關(guān)和零階保持器假設(shè)輸入端加一個采樣開關(guān)和零階保持器Foh假設(shè)輸入端加一個采樣開關(guān)和一階保持器。假設(shè)輸入端加一個采樣開關(guān)和一階保持器。Tustin采用雙線性變換（采用雙線性變換（Tustin算法）方法算法）方法Prewarp