連續(xù)系統(tǒng)數(shù)字仿真1_08

1、1第第3 3章章 連續(xù)系統(tǒng)的數(shù)字仿真連續(xù)系統(tǒng)的數(shù)字仿真 本章內(nèi)容本章內(nèi)容(1) 熟悉在數(shù)字計(jì)算機(jī)仿真技術(shù)中常用的幾種熟悉在數(shù)字計(jì)算機(jī)仿真技術(shù)中常用的幾種數(shù)值積分法數(shù)值積分法， 特別是四階龍格特別是四階龍格-庫(kù)塔法；庫(kù)塔法；(2) 典型環(huán)節(jié)典型環(huán)節(jié)及其系數(shù)矩陣的及其系數(shù)矩陣的確定確定；(3) 各各連接矩陣連接矩陣的確定；的確定；(4) 利用利用MATLAB在四階龍格在四階龍格-庫(kù)塔法的基礎(chǔ)上，對(duì)庫(kù)塔法的基礎(chǔ)上，對(duì)以狀態(tài)以狀態(tài) 空間表達(dá)式和方框圖描述的連續(xù)系統(tǒng)空間表達(dá)式和方框圖描述的連續(xù)系統(tǒng)進(jìn)行進(jìn)行仿真仿真；(5) 了解以了解以增廣矩陣法增廣矩陣法為基礎(chǔ)的連續(xù)系統(tǒng)的為基礎(chǔ)的連續(xù)系統(tǒng)的快速仿真快速仿

2、真方法方法。 2 用數(shù)字計(jì)算機(jī)來(lái)用數(shù)字計(jì)算機(jī)來(lái)仿真仿真或模擬一個(gè)連續(xù)控或模擬一個(gè)連續(xù)控制系統(tǒng)的目的就是制系統(tǒng)的目的就是求解系統(tǒng)的數(shù)學(xué)模型求解系統(tǒng)的數(shù)學(xué)模型。由控制理論知，一個(gè)由控制理論知，一個(gè)n n階連續(xù)系統(tǒng)可以被描階連續(xù)系統(tǒng)可以被描述成由述成由n n個(gè)積分器組成的模擬結(jié)構(gòu)圖。因此個(gè)積分器組成的模擬結(jié)構(gòu)圖。因此利用數(shù)字計(jì)算機(jī)來(lái)進(jìn)行連續(xù)系統(tǒng)的仿真，利用數(shù)字計(jì)算機(jī)來(lái)進(jìn)行連續(xù)系統(tǒng)的仿真，從本質(zhì)上講就是要在數(shù)字計(jì)算機(jī)上構(gòu)造出從本質(zhì)上講就是要在數(shù)字計(jì)算機(jī)上構(gòu)造出n n個(gè)數(shù)字積分器，也就是個(gè)數(shù)字積分器，也就是讓數(shù)字計(jì)算機(jī)進(jìn)行讓數(shù)字計(jì)算機(jī)進(jìn)行n n次數(shù)值積分運(yùn)算次數(shù)值積分運(yùn)算?？梢?jiàn)，連續(xù)系統(tǒng)數(shù)字仿?？梢?jiàn)，連

3、續(xù)系統(tǒng)數(shù)字仿真中的最基本的算法是數(shù)值積分算法。真中的最基本的算法是數(shù)值積分算法。343.1 3.1 數(shù)值積分法數(shù)值積分法 連續(xù)系統(tǒng)通常連續(xù)系統(tǒng)通常把數(shù)學(xué)模型化為狀態(tài)空間表達(dá)式把數(shù)學(xué)模型化為狀態(tài)空間表達(dá)式，為，為了對(duì)了對(duì)n階連續(xù)系統(tǒng)在數(shù)字計(jì)算機(jī)上仿真及求解，就要階連續(xù)系統(tǒng)在數(shù)字計(jì)算機(jī)上仿真及求解，就要采用數(shù)值積分法來(lái)求解系統(tǒng)數(shù)學(xué)模型中的采用數(shù)值積分法來(lái)求解系統(tǒng)數(shù)學(xué)模型中的n個(gè)一階微個(gè)一階微分方程。分方程。 設(shè)設(shè)n階連續(xù)系統(tǒng)由以下階連續(xù)系統(tǒng)由以下n個(gè)一階微分方程組成個(gè)一階微分方程組成 （3-13-1） 所謂數(shù)值積分法，就是要逐個(gè)求出區(qū)間所謂數(shù)值積分法，就是要逐個(gè)求出區(qū)間a ,b內(nèi)若內(nèi)若干個(gè)離散點(diǎn)干

4、個(gè)離散點(diǎn)a t0 t1 t0 0時(shí)，時(shí)，x( (t) )是未知的，因此式（是未知的，因此式（3-23-2）右端的積分是求不出的。為了解決這個(gè)問(wèn)題，右端的積分是求不出的。為了解決這個(gè)問(wèn)題，把積分間隔取得足夠小，使得在把積分間隔取得足夠小，使得在tk與與tk+1+1之間之間的的f( (t, ,x( (t)可以近似看作常數(shù)可以近似看作常數(shù)f( (tk,x( (tk)，這樣便得到用矩形公式積分的近似公式這樣便得到用矩形公式積分的近似公式或簡(jiǎn)化為或簡(jiǎn)化為 這就是這就是歐拉公式歐拉公式。htxtftxtxkkkk)(,()()(1hxtfxxkkkk),(17 以以x(t0)=x0作為初始值，應(yīng)用歐拉公式

5、，作為初始值，應(yīng)用歐拉公式，就可以一步步求出每一時(shí)刻就可以一步步求出每一時(shí)刻tk的的xk值，值， 即即 k =0,=0,x1 1x0 0+ +f( (t0 0, ,x0 0) )h k=1, =1, x2 2x1 1+ +f(t1 1, ,x1 1) )h k= =n-1-1, xn nxn-1n-1+ +f( (tn-1n-1, ,xn-1n-1) )h 這樣式這樣式(3-1)(3-1)的解的解x( (t) )就求出來(lái)了。歐拉法的計(jì)就求出來(lái)了。歐拉法的計(jì)算雖然比較簡(jiǎn)單，但精度算雖然比較簡(jiǎn)單，但精度較低。圖較低。圖3-13-1為歐拉法的為歐拉法的幾何解釋。幾何解釋。xx0 x2x1t0t2t1

6、(t1 ,x1)(t2 ,x2)t圖3 -1 歐 拉法的幾何解釋83.1.2 3.1.2 梯形法梯形法 由上可知由上可知?dú)W拉公式歐拉公式中的積分是用中的積分是用矩形面積矩形面積f( (tk, ,xk) )h 來(lái)來(lái)近似近似的。的。 由圖由圖3-23-2知，用矩形面積知，用矩形面積tkabtk+1+1代替積分，其誤差就代替積分，其誤差就是圖中陰影部分。為了提高精是圖中陰影部分。為了提高精度現(xiàn)度現(xiàn)用梯形面積用梯形面積tkactk+1+1來(lái)代來(lái)代替積分替積分，即，即 于是可得梯形法的計(jì)算公式為于是可得梯形法的計(jì)算公式為 ),(),(2)(, (111kkkkttxtfxtfhdttxtfkk),(),

7、(2111kkkkkkxtfxtfhxxtdtdxf kt1ktabc圖3 -2 梯 形法的幾何意義9 由于上式由于上式右邊包含未知量右邊包含未知量xk+1+1, ,所以每一步都所以每一步都必須通過(guò)必須通過(guò)迭代求解迭代求解，每一步迭代的初值每一步迭代的初值xk+1+1(0)(0)通常采用歐拉公式來(lái)計(jì)算通常采用歐拉公式來(lái)計(jì)算，因此梯形法每一步，因此梯形法每一步迭代公式為迭代公式為 (3-3(3-3）式中式中 迭代次數(shù)迭代次數(shù)R=0,1,2,=0,1,2,),(),(2),()(11)1(1)0(1RkkkkkRkkkkkxtfxtfhxxxthfxx103.1.3 3.1.3 預(yù)估校正法預(yù)估校正

8、法雖然雖然梯形法梯形法比歐拉法精確，但是由于每比歐拉法精確，但是由于每一步都要進(jìn)行多次疊代，計(jì)算量大，為了一步都要進(jìn)行多次疊代，計(jì)算量大，為了簡(jiǎn)化計(jì)算，有時(shí)只對(duì)式（簡(jiǎn)化計(jì)算，有時(shí)只對(duì)式（3-33-3）進(jìn)行）進(jìn)行一次疊一次疊代代就可以了，因此可得就可以了，因此可得 通常稱(chēng)這類(lèi)方法為預(yù)估校正方法。它首通常稱(chēng)這類(lèi)方法為預(yù)估校正方法。它首先根據(jù)歐拉公式計(jì)算出先根據(jù)歐拉公式計(jì)算出xk+1的預(yù)估值的預(yù)估值xk+1(0)(0)，然后再對(duì)它進(jìn)行校正，以得到更準(zhǔn)確的，然后再對(duì)它進(jìn)行校正，以得到更準(zhǔn)確的近似值近似值xk+1(1)(1)。),(),(2),()0(11)1(1k)0(1kkkkkkkkkxtfxtf

9、hxxxthfxx113.1.4 3.1.4 龍格庫(kù)塔法龍格庫(kù)塔法 根據(jù)泰勒級(jí)數(shù)將式（根據(jù)泰勒級(jí)數(shù)將式（3-13-1）在）在tk+1=tk+h時(shí)刻的解時(shí)刻的解xk+1=x(tk+h) 在在tk附近展開(kāi)，有附近展開(kāi)，有 3-53-5） 可以看出，提高截?cái)嗾`差的階次，便可提高其精可以看出，提高截?cái)嗾`差的階次，便可提高其精度，但是由于計(jì)算各階導(dǎo)數(shù)相當(dāng)麻煩，所以直接采度，但是由于計(jì)算各階導(dǎo)數(shù)相當(dāng)麻煩，所以直接采用泰勒級(jí)數(shù)公式是不適用的，為了解決提高精度問(wèn)用泰勒級(jí)數(shù)公式是不適用的，為了解決提高精度問(wèn)題，龍格和庫(kù)塔兩人先后提出了間接使用泰勒級(jí)數(shù)題，龍格和庫(kù)塔兩人先后提出了間接使用泰勒級(jí)數(shù)公式的方法，即用函

10、數(shù)值公式的方法，即用函數(shù)值f (t,x)的線性組合來(lái)代替的線性組合來(lái)代替f (t,x)的導(dǎo)數(shù)，然后按泰勒公式確定其中的系數(shù)的導(dǎo)數(shù)，然后按泰勒公式確定其中的系數(shù), , 這這樣既能避免計(jì)算樣既能避免計(jì)算f (t,x)的導(dǎo)數(shù)，又可以提高數(shù)值計(jì)算的導(dǎo)數(shù)，又可以提高數(shù)值計(jì)算精度，其方法如下。精度，其方法如下。)( 0! 21 1)(21 ppkpkkkkhxphxhxhxx12因因故式（故式（3-53-5）可寫(xiě)成）可寫(xiě)成 （3-63-6） 為了避免計(jì)算式（為了避免計(jì)算式（3-63-6）中的各階導(dǎo)數(shù)項(xiàng)，可令）中的各階導(dǎo)數(shù)項(xiàng)，可令xk+1由由以下多項(xiàng)式表示。以下多項(xiàng)式表示。 （3-73-7）kxkxxttx

11、xttkkkkkfffdtdxxftfdtxtdfxfxtfxkkkkk )(),(),()(0)(! 2121pkxkkkkhfffhhfxxkvmmmkkkahxx1113 式中式中 am為待定因子，為待定因子，v為使用為使用f函數(shù)值的個(gè)數(shù)，函數(shù)值的個(gè)數(shù)，km滿足下列方程滿足下列方程 （3-83-8）即：即： 將式（將式（3-73-7）展開(kāi)成）展開(kāi)成h的冪級(jí)數(shù)并與微分方程式的冪級(jí)數(shù)并與微分方程式（3-13-1）精確解式（）精確解式（3-63-6）逐項(xiàng)比較）逐項(xiàng)比較, ,便可求得式便可求得式（3-73-7）和式（）和式（3-83-8）中的系數(shù)）中的系數(shù)am ,bmj和和cm等。等。vmchk

12、bxhctfkmjjmjkmkm, 2 , 1, 0),(111),(),(),(23213133121221hkbhkbxhctfkhkbxhctfkxtfkkkkkkk14現(xiàn)以現(xiàn)以v=2=2為例，來(lái)說(shuō)明這些參數(shù)的確定方法。為例，來(lái)說(shuō)明這些參數(shù)的確定方法。設(shè)設(shè)v=2=2，則有，則有 （3-93-9） 將將k1 1和和k k2 2在同一點(diǎn)（在同一點(diǎn)（tk ,xk）上用二元函數(shù)展開(kāi)為）上用二元函數(shù)展開(kāi)為),(),()(12122122111hkbxhctfkxtfkkakahxxkkkkkkkkkfxtfk),(1)( 0)( 0),(3212312122hff hbf hcfhxfhkbtfh

13、cxtfkkxkkxxttxxttkkkkkkk15將將k1 1和和k2 2代入式（代入式（3-93-9）整理后可得）整理后可得 （3-103-10）將上式與式（將上式與式（3-63-6）逐項(xiàng)進(jìn)行比較，可得以下關(guān)系式）逐項(xiàng)進(jìn)行比較，可得以下關(guān)系式若取若取 則則)( 0)()(3212222211hffbafcahfaahxxkxkkkkk21,21, 12122221bacaaa12c1,212121baa16于是可得于是可得 (3-11)(3-11) 由于式（由于式（3-113-11）只取到泰勒級(jí)數(shù)展開(kāi)式的）只取到泰勒級(jí)數(shù)展開(kāi)式的h2 2項(xiàng)，故稱(chēng)這種方法為項(xiàng)，故稱(chēng)這種方法為兩階龍格庫(kù)塔兩階龍

14、格庫(kù)塔法法，其截?cái)嗾`差為，其截?cái)嗾`差為0(0(h3 3) )。)(2211kkhxxkk),(),(121hkxhtfkxtfkkkkk17 同理當(dāng)同理當(dāng)v=4=4時(shí)，仿照上述方法可得如下時(shí)，仿照上述方法可得如下四階四階龍格龍格- -庫(kù)塔公式庫(kù)塔公式),()2,2()2,2(),()22(6342312143211hkxhtfkkhxhtfkkhxhtfkxtfkkkkkhxxkkkkkkkkkk18 通過(guò)上述龍格通過(guò)上述龍格- -庫(kù)塔法的介紹，可以把以上庫(kù)塔法的介紹，可以把以上介紹的幾種數(shù)值積分法統(tǒng)一起來(lái)，它們都是基介紹的幾種數(shù)值積分法統(tǒng)一起來(lái)，它們都是基于在初值附近展開(kāi)成泰勒級(jí)數(shù)的原理，所

15、不同于在初值附近展開(kāi)成泰勒級(jí)數(shù)的原理，所不同的是取泰勒級(jí)數(shù)多少項(xiàng)。的是取泰勒級(jí)數(shù)多少項(xiàng)。歐拉公式僅取到歐拉公式僅取到h項(xiàng)項(xiàng)，梯形法與二階龍格庫(kù)塔法相同均取到梯形法與二階龍格庫(kù)塔法相同均取到h2 2項(xiàng)，項(xiàng)，四四階龍格庫(kù)塔法取到階龍格庫(kù)塔法取到h4 4項(xiàng)項(xiàng)。從理論上講，取得的。從理論上講，取得的項(xiàng)數(shù)愈多，計(jì)算精度愈高，但計(jì)算量愈大，愈項(xiàng)數(shù)愈多，計(jì)算精度愈高，但計(jì)算量愈大，愈復(fù)雜，計(jì)算誤差也將增加，因此要適當(dāng)?shù)倪x擇。復(fù)雜，計(jì)算誤差也將增加，因此要適當(dāng)?shù)倪x擇。目前在數(shù)字仿真中，目前在數(shù)字仿真中，最常用的是四階龍格庫(kù)塔最常用的是四階龍格庫(kù)塔法，其截?cái)嗾`差為法，其截?cái)嗾`差為( (h5 5),), 已能滿足

16、仿真精度已能滿足仿真精度的要求。的要求。193.1.5 3.1.5 關(guān)于仿真數(shù)值積分法的幾點(diǎn)討關(guān)于仿真數(shù)值積分法的幾點(diǎn)討1 1單步法和多步法單步法和多步法 解解初值問(wèn)題初值問(wèn)題的數(shù)值解法的的數(shù)值解法的共同特點(diǎn)共同特點(diǎn)是是步進(jìn)式步進(jìn)式，即，即從最初一點(diǎn)或幾點(diǎn)出發(fā)，每一步根據(jù)從最初一點(diǎn)或幾點(diǎn)出發(fā)，每一步根據(jù)xk一點(diǎn)或前面一點(diǎn)或前面幾點(diǎn)幾點(diǎn)xk-1 , xk-2 ,來(lái)計(jì)算新的來(lái)計(jì)算新的xk+1的值，這樣逐步推的值，這樣逐步推進(jìn)。進(jìn)。 當(dāng)從當(dāng)從tk推進(jìn)到推進(jìn)到tk+1只需用只需用tk時(shí)刻的數(shù)據(jù)時(shí)，稱(chēng)為時(shí)刻的數(shù)據(jù)時(shí)，稱(chēng)為單步單步法，例如歐拉法和龍格庫(kù)塔法。法，例如歐拉法和龍格庫(kù)塔法。20 相反，需要用到

17、相反，需要用到tk以及以及過(guò)去時(shí)刻過(guò)去時(shí)刻tk-1 ,tk-2 ,的數(shù)據(jù)的數(shù)據(jù)時(shí)稱(chēng)為時(shí)稱(chēng)為多步法多步法。線性多步法的一般形式是。線性多步法的一般形式是 （3-133-13） 多步法不能從多步法不能從t=0=0自啟動(dòng)，通常需要選用相同階次自啟動(dòng)，通常需要選用相同階次精度的單步法來(lái)啟動(dòng)，獲得所需前精度的單步法來(lái)啟動(dòng)，獲得所需前k步數(shù)據(jù)后，方可步數(shù)據(jù)后，方可轉(zhuǎn)入相應(yīng)多步法，因多步法利用信息量大，因而比單轉(zhuǎn)入相應(yīng)多步法，因多步法利用信息量大，因而比單步法更精確。步法更精確。)(0111101nknkknknkkkfffhxxxx212 2顯式和隱式顯式和隱式 在計(jì)算在計(jì)算xk+1時(shí)公式時(shí)公式右端所用到

18、的數(shù)據(jù)均已知右端所用到的數(shù)據(jù)均已知時(shí)，稱(chēng)時(shí)，稱(chēng)為為顯式顯式算法。例如算法。例如歐拉法、龍格庫(kù)塔法歐拉法、龍格庫(kù)塔法和式（和式（3-133-13）中中-1-1=0=0的情況。相反，在算式右端中隱含有未知量的情況。相反，在算式右端中隱含有未知量xk+1時(shí)，稱(chēng)為時(shí)，稱(chēng)為隱式隱式算法。例如算法。例如梯形法、預(yù)估校正法梯形法、預(yù)估校正法和和式（式（3-133-13）中）中-1-100的情況。的情況。 顯式算法利用前幾步計(jì)算結(jié)果即可進(jìn)行遞推求解下顯式算法利用前幾步計(jì)算結(jié)果即可進(jìn)行遞推求解下步結(jié)果，因而易于計(jì)算。而隱式計(jì)算需要迭代法，先步結(jié)果，因而易于計(jì)算。而隱式計(jì)算需要迭代法，先用另一同階次顯式公式估計(jì)出

19、一個(gè)初值用另一同階次顯式公式估計(jì)出一個(gè)初值xk+1(0)(0), ,并求得并求得fk+1, ,然后再用隱式求得校正值然后再用隱式求得校正值xk+1(1)(1), , 若未達(dá)到所需若未達(dá)到所需精度要求，則再次迭代求解，直到兩次迭代值精度要求，則再次迭代求解，直到兩次迭代值xk+1 ( (i) )和和xk+1(i+1) 之間的誤差在要求的范圍內(nèi)為止，故之間的誤差在要求的范圍內(nèi)為止，故隱式隱式算法算法精度高精度高，對(duì)誤差有較強(qiáng)的抑制作用對(duì)誤差有較強(qiáng)的抑制作用。223.3.數(shù)值穩(wěn)定性與仿真誤差數(shù)值穩(wěn)定性與仿真誤差 仿真誤差與數(shù)值計(jì)算方法、計(jì)算機(jī)的精度以及仿真誤差與數(shù)值計(jì)算方法、計(jì)算機(jī)的精度以及計(jì)算步長(zhǎng)

20、的選擇有關(guān)。當(dāng)計(jì)算方程和計(jì)算機(jī)確定以計(jì)算步長(zhǎng)的選擇有關(guān)。當(dāng)計(jì)算方程和計(jì)算機(jī)確定以后，則僅與計(jì)算步長(zhǎng)有關(guān)，所以在仿真中計(jì)算步長(zhǎng)后，則僅與計(jì)算步長(zhǎng)有關(guān)，所以在仿真中計(jì)算步長(zhǎng)是一個(gè)重要的參數(shù)。仿真誤差一般有如下兩種：是一個(gè)重要的參數(shù)。仿真誤差一般有如下兩種：（1 1）截?cái)嗾`差截?cái)嗾`差由于仿真模由于仿真模型僅是原系統(tǒng)型僅是原系統(tǒng)模型模型的一種的一種逼近逼近，以及各種數(shù)值積分法的計(jì)算都以及各種數(shù)值積分法的計(jì)算都是近似的算法。通常計(jì)算步長(zhǎng)是近似的算法。通常計(jì)算步長(zhǎng)愈小，截?cái)嗾`差也愈小。愈小，截?cái)嗾`差也愈小。（2 2）舍入誤差舍入誤差由于由于計(jì)算機(jī)計(jì)算機(jī)的精度有限的精度有限（有限位數(shù)）所產(chǎn)（有限位數(shù)）所產(chǎn)生

21、。通常計(jì)算步長(zhǎng)愈小，計(jì)算生。通常計(jì)算步長(zhǎng)愈小，計(jì)算次數(shù)愈多次數(shù)愈多, , 舍入誤差愈大。舍入誤差愈大。圖3 -3 仿 真誤差曲線h 0h誤差總誤差截?cái)嗾`差舍入誤差023 對(duì)截?cái)嗾`差而言，計(jì)算步長(zhǎng)愈小愈好，但太小不但對(duì)截?cái)嗾`差而言，計(jì)算步長(zhǎng)愈小愈好，但太小不但會(huì)增加計(jì)算時(shí)間，而且由于舍入誤差的增加，不一定會(huì)增加計(jì)算時(shí)間，而且由于舍入誤差的增加，不一定能達(dá)到提高精度的目的，甚至可能出現(xiàn)數(shù)值不穩(wěn)情況。能達(dá)到提高精度的目的，甚至可能出現(xiàn)數(shù)值不穩(wěn)情況。顯然計(jì)算步長(zhǎng)太大，不但精度不能滿足要求，而且計(jì)顯然計(jì)算步長(zhǎng)太大，不但精度不能滿足要求，而且計(jì)算步長(zhǎng)超過(guò)該算法的判穩(wěn)條件時(shí)，也會(huì)出現(xiàn)不穩(wěn)定情算步長(zhǎng)超過(guò)該算法

22、的判穩(wěn)條件時(shí)，也會(huì)出現(xiàn)不穩(wěn)定情況。由此可見(jiàn)，計(jì)算步長(zhǎng)只能在某一范圍內(nèi)選擇，圖況。由此可見(jiàn)，計(jì)算步長(zhǎng)只能在某一范圍內(nèi)選擇，圖中的中的h0 0為最佳計(jì)算步長(zhǎng)。為最佳計(jì)算步長(zhǎng)。 一般一般控制系統(tǒng)控制系統(tǒng)的輸出的輸出動(dòng)態(tài)響應(yīng)動(dòng)態(tài)響應(yīng)在在開(kāi)始段變化較快開(kāi)始段變化較快，到到最后變化將會(huì)很緩慢最后變化將會(huì)很緩慢。這時(shí)，計(jì)算可以。這時(shí)，計(jì)算可以采用變步長(zhǎng)采用變步長(zhǎng)的方法的方法，即在開(kāi)始階段步長(zhǎng)取得小一些，在最后階段，即在開(kāi)始階段步長(zhǎng)取得小一些，在最后階段取得大一些，這樣即可以保證計(jì)算的取得大一些，這樣即可以保證計(jì)算的精度精度，也可以加，也可以加快計(jì)算的快計(jì)算的速度速度。24 對(duì)于一般工程計(jì)算，計(jì)算精度要求并不

23、對(duì)于一般工程計(jì)算，計(jì)算精度要求并不太高，故常用太高，故常用定步長(zhǎng)定步長(zhǎng)的方法。作為的方法。作為經(jīng)驗(yàn)經(jīng)驗(yàn)數(shù)據(jù)，當(dāng)采用四階龍格庫(kù)塔法作數(shù)值數(shù)據(jù)，當(dāng)采用四階龍格庫(kù)塔法作數(shù)值積分計(jì)算時(shí)，取計(jì)算步長(zhǎng)積分計(jì)算時(shí)，取計(jì)算步長(zhǎng) h=tr r/10/10或或ts/40/40 式中式中 tr系統(tǒng)在階躍函數(shù)作用下的上升系統(tǒng)在階躍函數(shù)作用下的上升時(shí)間；時(shí)間；ts系統(tǒng)在階躍函數(shù)作用下的過(guò)渡系統(tǒng)在階躍函數(shù)作用下的過(guò)渡過(guò)程時(shí)間。若系統(tǒng)有多個(gè)回路，則應(yīng)按過(guò)程時(shí)間。若系統(tǒng)有多個(gè)回路，則應(yīng)按反應(yīng)最快的回路考慮。反應(yīng)最快的回路考慮。253.2 3.2 連續(xù)系統(tǒng)的數(shù)字仿真程序連續(xù)系統(tǒng)的數(shù)字仿真程序若系統(tǒng)的狀態(tài)空間表達(dá)式為若系統(tǒng)的狀態(tài)

24、空間表達(dá)式為 （3-143-14） （3-153-15）其中其中 A: nn; b: n1; C:1nCxybuAxx26 假設(shè)在仿真中，數(shù)值積分法采用四階龍格假設(shè)在仿真中，數(shù)值積分法采用四階龍格庫(kù)塔方法，因?qū)τ趲?kù)塔方法，因?qū)τ趎 n階系統(tǒng)，其狀態(tài)方程式階系統(tǒng)，其狀態(tài)方程式（3-143-14）可寫(xiě)成以下）可寫(xiě)成以下n個(gè)一階微分方程個(gè)一階微分方程 (3-16)(3-16)ubxaxaxaxininiii.2211), 2 , 1(), (321nixxxxtfni27 故根據(jù)式（故根據(jù)式（3-123-12）可得）可得求解以上一階微求解以上一階微分方程組的四階龍格庫(kù)塔公式分方程組的四階龍格庫(kù)塔公式

25、如下如下), 2 , 1()22(64321)()1(nikkkkhxxiiiikiki)()()()2()2()2()2()2()2()(33114221131111222111htubhkxahkxakhtubkhxakhxakhtubkhxakhxaktubxaxaxakkiikninikiikiikninikiikiikninikiikikninkikii式中式中 xik為為t=tk時(shí)刻的時(shí)刻的xi值，值，xik+1為為t=tk+h時(shí)時(shí)刻的刻的xi值。值。28令令TnTnTnTnTknkkkTknkkkkkkKkkkKkkkKkkkKxxxxxxxx424144323133222122

26、1211112111211129則式（則式（3-173-17）可寫(xiě)成如下矩陣的形式）可寫(xiě)成如下矩陣的形式根據(jù)式（根據(jù)式（3-153-15）可得）可得t=tk+1時(shí)刻的輸出時(shí)刻的輸出)()()2()2()2()2()()22(6342312143211htbuhKxAKhtbuKhxAKhtbuKhxAKtbuAxKKKKKhxxkkkkkkkkkk11kkCxy30給定系統(tǒng)參數(shù)并根據(jù)需要轉(zhuǎn)化為狀態(tài)空間表達(dá)式輸入仿真時(shí)間和計(jì)算步長(zhǎng)根據(jù)龍格庫(kù)塔法求 K1,K2,K3,K4求 和 值1kx1kyt=Tf輸出結(jié)果Ny圖3 -4 連 續(xù)系統(tǒng)的數(shù)字仿真程序框圖開(kāi)始給定輸入信號(hào)31例例3-13-1 假設(shè)單變

27、量系統(tǒng)如下圖所示。假設(shè)單變量系統(tǒng)如下圖所示。 試根據(jù)四階龍格庫(kù)塔法，試根據(jù)四階龍格庫(kù)塔法，求輸出量求輸出量y的動(dòng)的動(dòng)態(tài)響應(yīng)。態(tài)響應(yīng)。yr=20.40.002+ 圖3 -5-5 .638 .2132 .2045410562. 110875. 123466sssss32解解 仿真程序如下仿真程序如下ex3_1ex3_1 取仿真時(shí)間：取仿真時(shí)間：Tf f=5=5； 計(jì)算步長(zhǎng)：計(jì)算步長(zhǎng)：h=0.02 =0.02 在在MATLABMATLAB環(huán)境下執(zhí)行以上程序可得如圖環(huán)境下執(zhí)行以上程序可得如圖3-63-6所示仿真所示仿真曲線。曲線。33圖3-634 圖中帶數(shù)字的方框表示典型環(huán)節(jié)，圖中帶數(shù)字的方框表示典型

28、環(huán)節(jié)，, , ,表表示比例系數(shù)。示比例系數(shù)。V352r1r2r3+-y1y2x3 u4u3x2x1u2u1x4x5u5圖3 -9 多 變量系統(tǒng)結(jié)構(gòu)圖 +353.3 3.3 面向系統(tǒng)結(jié)構(gòu)圖的仿真面向系統(tǒng)結(jié)構(gòu)圖的仿真 這種方法與上節(jié)介紹的方法相比，有以下這種方法與上節(jié)介紹的方法相比，有以下幾個(gè)主要優(yōu)點(diǎn)：幾個(gè)主要優(yōu)點(diǎn)：1）便于研究各環(huán)節(jié)參數(shù)對(duì)系統(tǒng)的影響；便于研究各環(huán)節(jié)參數(shù)對(duì)系統(tǒng)的影響；2）可以得到每個(gè)環(huán)節(jié)的動(dòng)態(tài)響應(yīng)；可以得到每個(gè)環(huán)節(jié)的動(dòng)態(tài)響應(yīng)；3）可對(duì)多變量系統(tǒng)進(jìn)行仿真?？蓪?duì)多變量系統(tǒng)進(jìn)行仿真。下面具體介紹面向結(jié)構(gòu)圖的仿真方法。下面具體介紹面向結(jié)構(gòu)圖的仿真方法。363.3.1 3.3.1 典型環(huán)節(jié)的

29、確定典型環(huán)節(jié)的確定 一個(gè)控制系統(tǒng)可能由各種各樣的環(huán)節(jié)所一個(gè)控制系統(tǒng)可能由各種各樣的環(huán)節(jié)所組成，但比較常見(jiàn)環(huán)節(jié)有：組成，但比較常見(jiàn)環(huán)節(jié)有：（1）比例環(huán)節(jié)：比例環(huán)節(jié)： （2）積分環(huán)節(jié)：積分環(huán)節(jié)： （3）比例積分環(huán)節(jié)：比例積分環(huán)節(jié)：（4 4）慣性環(huán)節(jié)：）慣性環(huán)節(jié)：（5）超前滯后環(huán)節(jié)：超前滯后環(huán)節(jié)：（6 6）二階振蕩環(huán)節(jié)：）二階振蕩環(huán)節(jié)： 37 為了編制比較簡(jiǎn)單而且通用的仿真程序必須恰當(dāng)選為了編制比較簡(jiǎn)單而且通用的仿真程序必須恰當(dāng)選擇仿真環(huán)節(jié)。在這里選用圖擇仿真環(huán)節(jié)。在這里選用圖3-73-7所示的典型環(huán)節(jié)作為所示的典型環(huán)節(jié)作為仿真環(huán)節(jié)，即仿真環(huán)節(jié)，即式中式中 u為典型環(huán)節(jié)的輸入，為典型環(huán)節(jié)的輸入，x為

30、典型環(huán)節(jié)的輸出。為典型環(huán)節(jié)的輸出。 利用這個(gè)典型環(huán)節(jié)，只要利用這個(gè)典型環(huán)節(jié)，只要改變改變a ,b ,c 和和d 參數(shù)的參數(shù)的值值，便可分別表示以上所述的各一階環(huán)節(jié)，至于二階，便可分別表示以上所述的各一階環(huán)節(jié)，至于二階振蕩環(huán)節(jié)，則可用振蕩環(huán)節(jié)，則可用兩個(gè)一階環(huán)節(jié)等效連兩個(gè)一階環(huán)節(jié)等效連接得到，如圖接得到，如圖3-83-8所示。所示。bsadscsUsXsG)()()(xu圖3-7 典型環(huán)節(jié)bsadscx21TsTs1ku+-圖3-8 二階振蕩環(huán)節(jié)的等效結(jié)構(gòu)圖383.3.2 3.3.2 連接矩陣連接矩陣一個(gè)控制系統(tǒng)用典型環(huán)節(jié)來(lái)描述時(shí)，必一個(gè)控制系統(tǒng)用典型環(huán)節(jié)來(lái)描述時(shí)，必須用連接矩陣把各個(gè)典型環(huán)節(jié)連

31、接起來(lái)。須用連接矩陣把各個(gè)典型環(huán)節(jié)連接起來(lái)。所謂連接矩陣，就是用矩陣的形式表示各所謂連接矩陣，就是用矩陣的形式表示各個(gè)典型環(huán)節(jié)之間的關(guān)系。下面介紹連接矩個(gè)典型環(huán)節(jié)之間的關(guān)系。下面介紹連接矩陣的建立方法陣的建立方法, , 假設(shè)多輸入多輸出系統(tǒng)的假設(shè)多輸入多輸出系統(tǒng)的結(jié)構(gòu)圖如圖結(jié)構(gòu)圖如圖3-93-9所示。所示。39 圖中帶數(shù)字的方框表示典型環(huán)節(jié)，圖中帶數(shù)字的方框表示典型環(huán)節(jié)，, , ,表表示比例系數(shù)。示比例系數(shù)。V352r1r2r3+-y1y2x3 u4u3x2x1u2u1x4x5u5圖3 -9 多 變量系統(tǒng)結(jié)構(gòu)圖 +40 由圖可得由圖可得各環(huán)節(jié)的輸入各環(huán)節(jié)的輸入與與各環(huán)節(jié)的輸出各環(huán)節(jié)的輸出間的間

32、的關(guān)系關(guān)系以及系統(tǒng)的以及系統(tǒng)的輸出與環(huán)節(jié)的輸出輸出與環(huán)節(jié)的輸出間間關(guān)系關(guān)系分別為分別為和和 453433323222121551xuxurxxurxxurxu3241xyxy41寫(xiě)成寫(xiě)成矩陣矩陣形式形式和和32154321325543210000001000100010100000100001000010000rrrxxxxxuuuuu54321210010001000 xxxxxyy42 或?qū)懗苫驅(qū)懗?(3-20)(3-20) 定義式中的定義式中的W ，W0 和和 Wc陣為連接矩陣陣為連接矩陣，W反映了各典型環(huán)節(jié)輸入輸出間的連接關(guān)系，反映了各典型環(huán)節(jié)輸入輸出間的連接關(guān)系，W0反映了系統(tǒng)的參考

33、輸入與各環(huán)節(jié)輸入間的反映了系統(tǒng)的參考輸入與各環(huán)節(jié)輸入間的連接關(guān)系，連接關(guān)系， Wc c反映了系統(tǒng)的輸出與各環(huán)節(jié)輸反映了系統(tǒng)的輸出與各環(huán)節(jié)輸出間的關(guān)系。出間的關(guān)系。xWyrWWxuc043 一般也將系統(tǒng)中一般也將系統(tǒng)中各典型環(huán)節(jié)各典型環(huán)節(jié)的系數(shù)寫(xiě)成如下的系數(shù)寫(xiě)成如下矩陣的形式矩陣的形式（假設(shè)系統(tǒng)由（假設(shè)系統(tǒng)由n n個(gè)典型環(huán)節(jié)組成）個(gè)典型環(huán)節(jié)組成） (3-21)(3-21)nnnndcbadcbadcbaP22221111443.3.3 3.3.3 確定系統(tǒng)的狀態(tài)方程確定系統(tǒng)的狀態(tài)方程 典型環(huán)節(jié)和連接矩陣確定后典型環(huán)節(jié)和連接矩陣確定后，便可，便可求求得得系系統(tǒng)的統(tǒng)的狀態(tài)空間表達(dá)式狀態(tài)空間表達(dá)式，推

34、導(dǎo)過(guò)程推導(dǎo)過(guò)程如下。如下。 假設(shè)系統(tǒng)由假設(shè)系統(tǒng)由n n個(gè)典型環(huán)節(jié)組成，則根據(jù)典個(gè)典型環(huán)節(jié)組成，則根據(jù)典型環(huán)節(jié)的傳遞函數(shù)有型環(huán)節(jié)的傳遞函數(shù)有 （ i=1,2,）即即 sbasdcsUsXsGiiiiiii)()()()()()()(sUsdcsXsbaiiiiii45寫(xiě)成矩陣形寫(xiě)成矩陣形 （3-223-22）式中：式中：nnbbbBaaaA000000,0000002121nndddDcccC000000,0000002121)()()()(sUsDCsXsBA46 將式（將式（3-203-20）中的上式進(jìn)行拉氏變換后代）中的上式進(jìn)行拉氏變換后代入式（入式（3-223-22）中可得）中可得 對(duì)上

35、式兩邊取拉氏反變換得對(duì)上式兩邊取拉氏反變換得（3-233-23） 若參考輸入向量若參考輸入向量r=r1 1 r2 2 rm T T中中的的, ,r2 2,rm均為階躍函數(shù)，則上式可簡(jiǎn)化為均為階躍函數(shù)，則上式可簡(jiǎn)化為 （3-243-24）)()()()()()()()()()()(000ssRWDsRWCsXAWCssXWDBsRWsWXsDCsXsBArWDrWCxAWCxWDB00)()(rWCxAWCxWDB0)()(47令令則式（則式（3-243-24）可寫(xiě)成）可寫(xiě)成若若H的逆存在，則有的逆存在，則有再令再令可得可得 （3-253-25） 上式即為閉環(huán)系統(tǒng)的狀態(tài)方程，它是一個(gè)典型的上式即

36、為閉環(huán)系統(tǒng)的狀態(tài)方程，它是一個(gè)典型的狀態(tài)方程，利用前面介紹的求解方法可方便地求出各狀態(tài)方程，利用前面介紹的求解方法可方便地求出各典型環(huán)節(jié)的輸出響應(yīng)，最后根據(jù)式（典型環(huán)節(jié)的輸出響應(yīng)，最后根據(jù)式（3-203-20）中的第二）中的第二式便可求出系統(tǒng)的輸出響應(yīng)。式便可求出系統(tǒng)的輸出響應(yīng)。AWCQWDBH,rWCQxxH0rWCHQxHx011011,WCHBQHABrAxx483.3.4 3.3.4 面向結(jié)構(gòu)圖的數(shù)字仿真程序面向結(jié)構(gòu)圖的數(shù)字仿真程序 面向結(jié)構(gòu)圖面向結(jié)構(gòu)圖的數(shù)字仿真的數(shù)字仿真程序框圖如程序框圖如圖圖-10-10所所示，其程序示，其程序清單通過(guò)下清單通過(guò)下例給出。例給出。給定環(huán)節(jié)參數(shù)及連接

37、矩陣輸入仿真時(shí)間和計(jì)算步長(zhǎng)求 和 Q 陣值1HH ，求A ,B 陣t =Tf輸出結(jié)果根據(jù)龍格庫(kù)塔法求狀態(tài)方程的解Ny圖3 -10 面 向結(jié)構(gòu)圖的數(shù)字仿真程序框圖開(kāi)始給定輸入信號(hào)49 例例3-23-2 假設(shè)某一系統(tǒng)由四個(gè)典型環(huán)節(jié)組成，假設(shè)某一系統(tǒng)由四個(gè)典型環(huán)節(jié)組成，如圖如圖3-113-11所示。求輸出量所示。求輸出量y的動(dòng)態(tài)響應(yīng)。的動(dòng)態(tài)響應(yīng)。r=10圖3-111 . 05 . 0sss122s1010sy+-x1 u2x2 u3x3 u4u1x450 解解 由圖可得各由圖可得各環(huán)節(jié)環(huán)節(jié)的的輸入與輸出輸入與輸出以及系統(tǒng)的以及系統(tǒng)的輸出與環(huán)節(jié)的輸出間關(guān)系為輸出與環(huán)節(jié)的輸出間關(guān)系為rxxxxuuuu0

38、00101000010000110004321432143211000 xxxxy51 根據(jù)以上兩式和各典型環(huán)節(jié)的系數(shù)值，根據(jù)以上兩式和各典型環(huán)節(jié)的系數(shù)值，可得如下可得如下連接矩陣和系數(shù)矩陣連接矩陣和系數(shù)矩陣1000,000101000010000110000cWWW0101100212011015 . 011 . 04444333322221111dcbadcbadcbadcbaP52取仿真時(shí)間取仿真時(shí)間: : Tf=10；計(jì)算步長(zhǎng)；計(jì)算步長(zhǎng): : h=0.05 在在MATLABMATLAB環(huán)境下執(zhí)行以上程序可得如圖環(huán)境下執(zhí)行以上程序可得如圖3-123-12所示仿真曲線。所示仿真曲線。仿真程

39、序如下仿真程序如下ex3_2.mex3_2.m53圖3-12543.4 3.4 連續(xù)系統(tǒng)的快速仿真連續(xù)系統(tǒng)的快速仿真前面介紹過(guò)的兩種連續(xù)系統(tǒng)的數(shù)字仿真前面介紹過(guò)的兩種連續(xù)系統(tǒng)的數(shù)字仿真方法，當(dāng)系統(tǒng)比較復(fù)雜并要求滿足較高的方法，當(dāng)系統(tǒng)比較復(fù)雜并要求滿足較高的計(jì)算精度時(shí)，計(jì)算工作量較大，計(jì)算速度計(jì)算精度時(shí)，計(jì)算工作量較大，計(jì)算速度較慢，有時(shí)不能滿足實(shí)時(shí)仿真的要求，為較慢，有時(shí)不能滿足實(shí)時(shí)仿真的要求，為了解決這個(gè)問(wèn)題，下面介紹一種連續(xù)系統(tǒng)了解決這個(gè)問(wèn)題，下面介紹一種連續(xù)系統(tǒng)的快速數(shù)字仿真方法的快速數(shù)字仿真方法-增廣矩陣法增廣矩陣法。553.4.1 3.4.1 增廣矩陣法的基本原理增廣矩陣法的基本原理設(shè)連續(xù)系統(tǒng)的狀態(tài)方程為設(shè)連續(xù)系統(tǒng)的狀態(tài)方程為( (齊次齊次) ) （3-263-26）則它的解為則它的解為將將eAt展開(kāi)成泰勒級(jí)數(shù)，即展開(kāi)成泰勒級(jí)數(shù)，即則有則有)0()()()(0 xtxtAxtxt)0()(xetxAt3322! 31! 21tAtAAtIeAt)0(! 31! 21)(3322xtAtAAtItx56 可以證明，如果取可以證明，如果取eAt的泰勒級(jí)數(shù)的前五項(xiàng)，則式的泰勒級(jí)數(shù)的前五項(xiàng)，則式（3-273-27）的計(jì)算精度與四階龍格庫(kù)塔法相同，設(shè)計(jì)）的計(jì)算精度與四階龍格庫(kù)塔法相同，設(shè)計(jì)算步長(zhǎng)為算步長(zhǎng)為T(mén), , 則式（則式（3-273-27）可寫(xiě)成）可寫(xiě)成 上式括號(hào)中只