《控制系統(tǒng)數(shù)字仿真與cad》第3章連續(xù)系統(tǒng)的數(shù)字仿真

1、 本章內(nèi)容(1) 熟悉在數(shù)字計算機(jī)仿真技術(shù)中常用的幾種數(shù)值積分法， 特別是四階龍格-庫塔法；(2) 典型環(huán)節(jié)及其系數(shù)矩陣的確定；(3) 各連接矩陣的確定；(4) 利用MATLAB在四階龍格-庫塔法的基礎(chǔ)上，對以狀態(tài) 空間表達(dá)式和方框圖描述的連續(xù)系統(tǒng)進(jìn)行仿真；l 用數(shù)字計算機(jī)來仿真或模擬一個連續(xù)控制系統(tǒng)的目的就是求解系統(tǒng)的數(shù)學(xué)模型。由控制理論知，一個n階連續(xù)系統(tǒng)可以被描述成由n個積分器組成的模擬結(jié)構(gòu)圖。因此利用數(shù)字計算機(jī)來進(jìn)行連續(xù)系統(tǒng)的仿真，從本質(zhì)上講就是要在數(shù)字計算機(jī)上構(gòu)造出n個數(shù)字積分器，也就是讓數(shù)字計算機(jī)進(jìn)行n次數(shù)值積分運(yùn)算?？梢?，連續(xù)系統(tǒng)數(shù)字仿真中的最基本的算法是數(shù)值積分算法。l 連續(xù)系

2、統(tǒng)通常把數(shù)學(xué)模型化為狀態(tài)空間表達(dá)式，為了對n階連續(xù)系統(tǒng)在數(shù)字計算機(jī)上仿真及求解，就要采用數(shù)值積分法來求解系統(tǒng)數(shù)學(xué)模型中的n個一階微分方程。l 設(shè)n階連續(xù)系統(tǒng)所包含的n個一階微分方程中的第i個一階微分方程為l （3-1）l所謂數(shù)值積分法，就是要逐個求出區(qū)間a ,b內(nèi)若干個離散點(diǎn)a t0 t1 t0時，x(t)是未知的，因此式（3-2）右端的積分是求不出的。為了解決這個問題，把積分間隔取得足夠小，使得在tk與tk+1之間的f(t,x(t)可以近似看作常數(shù)f(tk,x(tk)，這樣便得到用矩形公式積分的近似公式l或簡化為l l這就是歐拉公式。htxtftxtxkkkk)(,()()(1hxtfxxk

3、kkk),(1l即 k =0,x1x0+f(t0,x0)hl k=1, x2x1+f(t1,x1)hl l k=n-1, xnxn-1+f(tn-1,xn-1)hl 這樣式(3-1)的解x(t)就求出來了。歐拉法的計算雖然比較簡單，但精度較低。圖3-1為歐拉法的幾何解釋。l由圖3-2知，用矩形面積tkabtk+1代替積分，其誤差就是圖中陰影部分。為了提高精度，現(xiàn)用梯形面積tkactk+1來代替積分，即l于是可得梯形法的計算公式為 111(, () ( , )(,)2kktkkkkthf t xt dtf t xf tx111 ( ,)(,)2kkkkkkhxxf txf txl 由于上式右邊包

4、含未知量xk+1,所以每一步都必須通過迭代求解，每一步迭代的初值xk+1(0)通常采用歐拉公式來計算，因此梯形法每一步迭代公式為l l （3-3）l式中 迭代次數(shù)R=0,1,2,(0)1(1)()111( ,) ( ,)(,)2kkkkRRkkkkkkxxhf txhxxf txf txl3.1.3 3.1.3 預(yù)估校正法預(yù)估校正法l雖然梯形法比歐拉法精確，但是由于每一步都要進(jìn)行多次疊代，計算量大，為了簡化計算，有時只對式（3-3）進(jìn)行一次疊代就可以了，因此可得l通常稱這類方法為預(yù)估校正方法。它首先根據(jù)歐拉公式計算出xk+1的預(yù)估值xk+1(0)，然后再對它進(jìn)行校正，以得到更準(zhǔn)確的近似值xk+

5、1(1)。),(),(2),()0(11)1(1k)0(1kkkkkkkkkxtfxtfhxxxthfxxl3.1.4 3.1.4 龍格庫塔法龍格庫塔法l 根據(jù)泰勒級數(shù)將式（3-1）在tk+1=tk+h時刻的解xk+1=x(tk+h) 在tk附近展開，有l(wèi) l （3-5）l 可以看出，提高截斷誤差的階次，便可提高其精度，但是由于計算各階導(dǎo)數(shù)相當(dāng)麻煩，所以直接采用泰勒級數(shù)公式是不適用的，為了解決提高精度問題，龍格和庫塔兩人先后提出了間接使用泰勒級數(shù)公式的方法，即用函數(shù)值f (t,x)的線性組合來代替f (t,x)的導(dǎo)數(shù)，然后按泰勒公式確定其中的系數(shù), 這樣既能避免計算f (t,x)的導(dǎo)數(shù)，又可以

6、提高數(shù)值計算精度，其方法如下。)( 0! 21 1)(21 ppkpkkkkhxphxhxhxxl因l故式（3-5）可寫成l （3-6） l為了避免計算式（3-6）中的各階導(dǎo)數(shù)項(xiàng)，可令xk+1由以下多項(xiàng)式表示。l （3-7）kxkxxttxxttkkkkkfffdtdxxftfdtxtdfxfxtfxkkkkk )(),(),()(0)(! 2121pkxkkkkhfffhhfxxkvmmmkkkahxx11l式中 am為待定因子，v為使用f函數(shù)值的個數(shù)，km滿足下列方程l （3-8）l即：l將式（3-7）展開成h的冪級數(shù)并與微分方程式（3-1）精確解式（3-6）逐項(xiàng)比較,便可求得式（3-7）

7、和式（3-8）中的系數(shù)am ,bmj和cm等。vmchkbxhctfkmjjmjkmkm, 2 , 1, 0),(111),(),(),(23213133121221hkbhkbxhctfkhkbxhctfkxtfkkkkkkkl現(xiàn)以v=2為例，來說明這些參數(shù)的確定方法。設(shè)v=2，則有l(wèi)（3-9）l l將k1和k2在同一點(diǎn)（tk ,xk）上用二元函數(shù)展開為),(),()(12122122111hkbxhctfkxtfkkakahxxkkkkkkkkkfxtfk),(1)( 0)( 0),(3212312122hff hbf hcfhxfhkbtfhcxtfkkxkkxxttxxttkkkkkk

8、kl將k1和k2代入式（3-9）整理后可得l l（3-10）l將上式與式（3-6）逐項(xiàng)進(jìn)行比較，可得以下關(guān)系式l若取 l則)( 0)()(3212222211hffbafcahfaahxxkxkkkkk21,21, 12122221bacaaa12c1,212121baal于是由式3-9可得l (3-11)l 由于式（3-11）只取到泰勒級數(shù)展開式的h2項(xiàng)，故稱這種方法為兩階龍格庫塔法，其截斷誤差為0(h3)。)(2211kkhxxkk),(),(121hkxhtfkxtfkkkkkl 同理當(dāng)v=4時，仿照上述方法可得如下四階龍格-庫塔公式),()2,2()2,2(),()22(6342312

9、143211hkxhtfkkhxhtfkkhxhtfkxtfkkkkkhxxkkkkkkkkkk(3-12)l 通過上述龍格-庫塔法的介紹，可以把以上介紹的幾種數(shù)值積分法統(tǒng)一起來，它們都是基于在初值附近展開成泰勒級數(shù)的原理，所不同的是取泰勒級數(shù)多少項(xiàng)。歐拉公式僅取到h項(xiàng)，梯形法與二階龍格庫塔法相同，均取到h2項(xiàng)，四階龍格庫塔法取到h4項(xiàng)。從理論上講，取得的項(xiàng)數(shù)愈多，計算精度愈高，但計算量愈大，愈復(fù)雜，計算誤差也將增加，因此要適當(dāng)?shù)倪x擇。目前在數(shù)字仿真中，最常用的是四階龍格庫塔法，其截斷誤差為(h5), 已能滿足仿真精度的要求。l3.1.5 3.1.5 關(guān)于仿真數(shù)值積分法的幾點(diǎn)討論關(guān)于仿真數(shù)值積

10、分法的幾點(diǎn)討論l1單步法和多步法l 解初值問題的數(shù)值解法的共同特點(diǎn)是步進(jìn)式，即從最初一點(diǎn)或幾點(diǎn)出發(fā)，每一步根據(jù)xk一點(diǎn)或前面幾點(diǎn)xk-1 , xk-2 ,來計算新的xk+1的值，這樣逐步推進(jìn)。l當(dāng)從tk推進(jìn)到tk+1只需用tk時刻的數(shù)據(jù)時，稱為單步法，例如歐拉法和龍格庫塔法。l 相反，需要用到tk以及過去時刻tk-1 ,tk-2 ,的數(shù)據(jù)時稱為多步法。線性多步法的一般形式是l （3-13）l 多步法不能從t=0自啟動，通常需要選用相同階次精度的單步法來啟動，獲得所需前k步數(shù)據(jù)后，方可轉(zhuǎn)入相應(yīng)多步法，因多步法利用信息量大，因而比單步法更精確。)(0111101nknkknknkkkfffhxxx

11、xl2顯式和隱式l 在計算xk+1時公式右端所用到的數(shù)據(jù)均已知時，稱為顯式算法。例如歐拉法、龍格庫塔法和式（3-13）中-1=0的情況。相反，在算式右端中隱含有未知量xk+1時，稱為隱式算法。例如梯形法、預(yù)估校正法和式（3-13）中-10的情況。l 顯式算法利用前幾步計算結(jié)果即可進(jìn)行遞推求解下步結(jié)果，因而易于計算。而隱式計算需要迭代法，先用另一同階次顯式公式估計出一個初值xk+1(0),并求得fk+1,然后再用隱式求得校正值xk+1(1), 若未達(dá)到所需精度要求，則再次迭代求解，直到兩次迭代值xk+1 (i)和xk+1(i+1) 之間的誤差在要求的范圍內(nèi)為止，故隱式算法精度高，對誤差有較強(qiáng)的抑

12、制作用。盡管隱式算法計算過程復(fù)雜，計算速度慢，但有時基于對精度、數(shù)值穩(wěn)定性等考慮，仍然經(jīng)常被使用，如求解病態(tài)方程等問題。l3.數(shù)值穩(wěn)定性與仿真誤差 數(shù)值積分法求解微分方程，實(shí)質(zhì)上是通過差分方程作為遞推公式進(jìn)行的，因此，在將微分方程差分化的變化過程中，應(yīng)保持原系統(tǒng)穩(wěn)定的特征，即要求用于計算的差分方程是穩(wěn)定的。但是，在計算機(jī)逐次計算時，初始數(shù)據(jù)的誤差及計算過程的舍入誤差等都會使誤差不斷積累。如果這種仿真誤差積累能夠抑制，不會隨計算時間增加而無限增大，則可以認(rèn)為相應(yīng)的計算方法是數(shù)值穩(wěn)定的。反之，則是數(shù)值不穩(wěn)定的。l（1）截斷誤差由于仿真模型僅是原系統(tǒng)模型的一種逼近，以及各種數(shù)值積分法的計算都是近似的

13、算法。通常計算步長愈小，截斷誤差也愈小。l（2）舍入誤差由于計算機(jī)的精度有限（有限位數(shù)）所產(chǎn)生。通常計算步長愈小，計算次數(shù)愈多, 舍入誤差愈大。l 對截斷誤差而言，計算步長愈小愈好，但太小不但會增加計算時間，而且由于舍入誤差的增加，不一定能達(dá)到提高精度的目的，甚至可能出現(xiàn)數(shù)值不穩(wěn)情況。顯然計算步長太大，不但精度不能滿足要求，而且計算步長超過該算法的判穩(wěn)條件時，也會出現(xiàn)不穩(wěn)定情況。由此可見，計算步長只能在某一范圍內(nèi)選擇，圖中的h0為最佳計算步長。l 一般控制系統(tǒng)的輸出動態(tài)響應(yīng)在開始段變化較快，到最后變化將會很緩慢。這時，計算可以采用變步長的方法，即在開始階段步長取得小一些，在最后階段取得大一些，

14、這樣即可以保證計算的精度，也可以加快計算的速度。l 對于一般工程計算，計算精度要求并不太高，故常用定步長的方法。作為經(jīng)驗(yàn)數(shù)據(jù)，當(dāng)采用四階龍格庫塔法作數(shù)值積分計算時，取計算步長l h=tr/10或ts/40l式中 tr系統(tǒng)在階躍函數(shù)作用下的上升時間；ts系統(tǒng)在階躍函數(shù)作用下的過渡過程時間。若系統(tǒng)有多個回路，則應(yīng)按反應(yīng)最快的回路考慮。自動控制系統(tǒng)常常是由許多環(huán)節(jié)組成的。應(yīng)用上節(jié)介紹的數(shù)字仿真方法對系統(tǒng)分析和研究 時，首先需要求出總的傳遞函數(shù)，再轉(zhuǎn)化為狀態(tài)空間表達(dá)式的形式，然后對其求解。當(dāng)改變系統(tǒng)某一環(huán)節(jié)的參數(shù)時，尤其是要改變小閉環(huán)中某一環(huán)節(jié)的參數(shù)時，以上整個過程又需重新計算，這對研究對象參數(shù)變化對

15、整個控制系統(tǒng)的影響是十分不便的。為了克服這些缺點(diǎn)，同時大多數(shù)從事自動化工作的科技人員更習(xí)慣于用結(jié)構(gòu)圖的形式來分析和研究控制系統(tǒng)，為此產(chǎn)生了面向結(jié)構(gòu)圖的仿真方法。該方法只需將各環(huán)節(jié)的參數(shù)及各環(huán)節(jié)間的連接方式輸入計算機(jī)，仿真程序就能自動求出閉環(huán)系統(tǒng)的狀態(tài)空間表達(dá)式。 本節(jié)主要介紹由典型環(huán)節(jié)參數(shù)和連接關(guān)系構(gòu)成閉環(huán)系統(tǒng)的狀態(tài)方程的方法，而動態(tài)響應(yīng)的計算，仍采用四階龍格-庫塔法。l 這種方法與上節(jié)介紹的方法相比，有以下幾個主要優(yōu)點(diǎn)：l）便于研究各環(huán)節(jié)參數(shù)對系統(tǒng)的影響；l）可以得到每個環(huán)節(jié)的動態(tài)響應(yīng)；l）可對多入多出系統(tǒng)進(jìn)行仿真。l下面具體介紹面向結(jié)構(gòu)圖的仿真方法。l3.3.1 3.3.1 典型環(huán)節(jié)的確定

16、典型環(huán)節(jié)的確定l 一個控制系統(tǒng)可能由各種各樣的環(huán)節(jié)所組成，但比較常見環(huán)節(jié)有：（1）比例環(huán)節(jié)：G(s)=k （2）積分環(huán)節(jié)： （3）比例積分環(huán)節(jié)：（4）慣性環(huán)節(jié)：（5）超前滯后環(huán)節(jié)：（6）二階振蕩環(huán)節(jié)： ( )kG ss2121( )kk skG skss( )1kG sTs211( )1T sG skTs22( )21kG sT sTsl 為了編制比較簡單而且通用的仿真程序，必須恰當(dāng)選擇仿真環(huán)節(jié)。在這里選用圖3-7所示的典型環(huán)節(jié)作為仿真環(huán)節(jié)，即l式中 u為典型環(huán)節(jié)的輸入，x為典型環(huán)節(jié)的輸出。l 利用這個典型環(huán)節(jié)，只要改變a ,b ,c 和d參數(shù)的值，便可分別表示以上所述的各一階環(huán)節(jié)，至于二階振

17、蕩環(huán)節(jié)，則可用l兩個一階環(huán)節(jié)等效連l接得到，如圖3-8所示。bsadscsUsXsG)()()(ux圖3-7 典型環(huán)節(jié)kux圖3-8 二階振蕩環(huán)節(jié)的等效結(jié)構(gòu)圖l3.3.2 3.3.2 連接矩陣連接矩陣l一個控制系統(tǒng)用典型環(huán)節(jié)來描述時，必須用連接矩陣把各個典型環(huán)節(jié)連接起來。所謂連接矩陣，就是用矩陣的形式表示各個典型環(huán)節(jié)之間的關(guān)系。下面介紹連接矩陣的建立方法, 假設(shè)多輸入多輸出系統(tǒng)的結(jié)構(gòu)圖如圖3-9所示。同理，三階及三階以上的環(huán)節(jié)也完全可以用若干個一階環(huán)節(jié)等效連接得到。由此可見，任何一個復(fù)雜的控制系統(tǒng)都可以用若干個典型環(huán)節(jié)來組成。l 圖中帶數(shù)字的方框表示典型環(huán)節(jié)， 表示比例系數(shù)。235, IIII

18、IIIV235Vl 由圖可得各環(huán)節(jié)的輸入與各環(huán)節(jié)的輸出間的關(guān)系，以及系統(tǒng)的輸出與環(huán)節(jié)的輸出間關(guān)系分別為l和 453433323222121551xuxurxxurxxurxu3241xyxyl寫成矩陣形式l和32154321325543210000001000100010100000100001000010000rrrxxxxxuuuuu54321210010001000 xxxxxyyl或?qū)懗蒷 (3-20)l l定義式中的W ，W0 和 Wc陣為連接矩陣，W反映了各典型環(huán)節(jié)輸入輸出間的連接關(guān)系，W0反映了系統(tǒng)的參考輸入與各環(huán)節(jié)輸入間的連接關(guān)系， Wc反映了系統(tǒng)的輸出與各環(huán)節(jié)輸出間的關(guān)系。x

19、WyrWWxuc0l 一般也將系統(tǒng)中各典型環(huán)節(jié)的系數(shù)寫成如下矩陣的形式（假設(shè)系統(tǒng)由n個典型環(huán)節(jié)組成）l (3-21)nnnndcbadcbadcbaP22221111l3.3.3 3.3.3 確定系統(tǒng)的狀態(tài)方程確定系統(tǒng)的狀態(tài)方程l 典型環(huán)節(jié)和連接矩陣確定后，便可求得系統(tǒng)的狀態(tài)空間表達(dá)式，推導(dǎo)過程如下。l 假設(shè)系統(tǒng)由n個典型環(huán)節(jié)組成，則根據(jù)典型環(huán)節(jié)的傳遞函數(shù)有l(wèi) l （ i=1,2,）l即l l（ i=1,2,）sbasdcsUsXsGiiiiiii)()()()()()()(sUsdcsXsbaiiiiii112200000000,0000nnabababAB112200000000,0000

20、nncdcdcdCD() ( ) () ( )ssssA B XC+D Ul 將式（3-20）中的上式進(jìn)行拉氏變換后代入式（3-22）中可得l對上式兩邊取拉氏反變換得l（3-23）l 若參考輸入向量r=r1 r2 rmT中的r1,r2,rm均為階躍函數(shù)，則上式可簡化為l （3-24）000() ( ) ()( )( )()( ) () ( )( )( )sssssssssssA+B XC+D WXWRB-DW XCWA XCWRDW R00(B-DW)x(CW - A)x CW rDW r0()()B - DW xCW - A xCW rl令l則式（3-24）可寫成l若H的逆存在，則有l(wèi)再令l

21、可得 （3-25）l 上式即為閉環(huán)系統(tǒng)的狀態(tài)方程，它是一個典型的狀態(tài)方程，利用前面介紹的求解方法可方便地求出各典型環(huán)節(jié)的輸出響應(yīng)，最后根據(jù)式（3-20）中的第二式便可求出系統(tǒng)的輸出響應(yīng)。H = B - DW ,Q = CW - A0Hx = Qx+ CW r0-1-1x = H Qx+ H CW r0-1-1A = H Q,B = HCWBrAxx在建立系統(tǒng)的各典型環(huán)節(jié)時應(yīng)注意以下兩點(diǎn)：（1）為保證H的逆H-1存在，應(yīng)嚴(yán)格按照的原則，確定每個典型環(huán)節(jié)。即避免以純比例、純微分環(huán)節(jié)作為典型環(huán)節(jié)。（2）在輸入向量不全為階躍函數(shù)的情況下，只要在確定典型環(huán)節(jié)時，注意使含有微分項(xiàng)系數(shù)（即 ）的環(huán)節(jié)不直接與

22、參考輸入連接，也可避免式3-23中出現(xiàn)r的導(dǎo)數(shù)。0ib 0id l 面向結(jié)構(gòu)圖的數(shù)字仿真程序框圖如圖-10所示，其程序清單通過下例給出。給定典型環(huán)節(jié)參數(shù)和連接矩陣輸入仿真時間Tf Tf和計算步長h開始給定輸入信號求H， H-1和Q陣值求A，B陣根據(jù)龍格-庫塔法求狀態(tài)方程的解計算系統(tǒng)輸出yT=TfyN輸出結(jié)果l例例3-23-2 假設(shè)某一系統(tǒng)由四個典型環(huán)節(jié)組成，如圖3-11所示。求輸出量y的動態(tài)響應(yīng)。r=10圖3-111 . 05 . 0sss122s1010sy+-x1 u2x2 u3x3 u4u1x4l解解 由圖可得各環(huán)節(jié)的輸入與輸出以及系統(tǒng)的輸出與環(huán)節(jié)的輸出間關(guān)系為rxxxxuuuu000101000010000110004321432143211000 xxxxyl根據(jù)以上兩式和各典型環(huán)節(jié)的系數(shù)值，可得如下連接矩陣和系數(shù)矩陣1000,000101000010000110000cWWW0101100212011015 . 011 . 04444333322221111dcbadcbadcbadcbaPl仿真程序如下l%