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

第三章連續(xù)系統(tǒng)的數(shù)字仿真本章內(nèi)容(1)

熟悉在數(shù)字計(jì)算機(jī)仿真技術(shù)中常用的幾種數(shù)值積分法，特別是四階龍格-庫塔法；(2)

典型環(huán)節(jié)及其系數(shù)矩陣的確定；(3)

各連接矩陣的確定；(4)

利用MATLAB在四階龍格-庫塔法的基礎(chǔ)上，對(duì)以狀態(tài)空間表達(dá)式和方框圖描述的連續(xù)系統(tǒng)進(jìn)行仿真；用數(shù)字計(jì)算機(jī)來仿真或模擬一個(gè)連續(xù)控制系統(tǒng)的目的就是求解系統(tǒng)的數(shù)學(xué)模型。由控制理論知，一個(gè)n階連續(xù)系統(tǒng)可以被描述成由n個(gè)積分器組成的模擬結(jié)構(gòu)圖。因此利用數(shù)字計(jì)算機(jī)來進(jìn)行連續(xù)系統(tǒng)的仿真，從本質(zhì)上講就是要在數(shù)字計(jì)算機(jī)上構(gòu)造出n個(gè)數(shù)字積分器，也就是讓數(shù)字計(jì)算機(jī)進(jìn)行n次數(shù)值積分運(yùn)算?？梢?，連續(xù)系統(tǒng)數(shù)字仿真中的最基本的算法是數(shù)值積分算法。3.1數(shù)值積分法連續(xù)系統(tǒng)通常把數(shù)學(xué)模型化為狀態(tài)空間表達(dá)式，為了對(duì)n階連續(xù)系統(tǒng)在數(shù)字計(jì)算機(jī)上仿真及求解，就要采用數(shù)值積分法來求解系統(tǒng)數(shù)學(xué)模型中的n個(gè)一階微分方程。設(shè)n階連續(xù)系統(tǒng)所包含的n個(gè)一階微分方程中的第i個(gè)一階微分方程為

（3-1）所謂數(shù)值積分法，就是要逐個(gè)求出區(qū)間［a,b］內(nèi)若干個(gè)離散點(diǎn)a≤t0

<t1<…<tn≤b處的近似值x(t1),x(t2),…,x(tn)。

當(dāng)t>t0時(shí)，x(t)是未知的，因此式（3-2）右端的積分是求不出的。為了解決這個(gè)問題，把積分間隔取得足夠小，使得在tk與tk+1之間的f(t,x(t))可以近似看作常數(shù)f(tk,x(tk))，這樣便得到用矩形公式積分的近似公式或簡化為

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

由上可知?dú)W拉公式中的積分是用矩形面積f(tk,xk)h來近似的。圖3-2為梯形法的幾何意義。由圖3-2知，用矩形面積tkabtk+1代替積分，其誤差就是圖中陰影部分。為了提高精度，現(xiàn)用梯形面積tkactk+1來代替積分，即于是可得梯形法的計(jì)算公式為由于上式右邊包含未知量xk+1,所以每一步都必須通過迭代求解，每一步迭代的初值xk+1(0)通常采用歐拉公式來計(jì)算，因此梯形法每一步迭代公式為

（3-3）式中迭代次數(shù)R=0,1,2,…3.1.3預(yù)估－校正法雖然梯形法比歐拉法精確，但是由于每一步都要進(jìn)行多次疊代，計(jì)算量大，為了簡化計(jì)算，有時(shí)只對(duì)式（3-3）進(jìn)行一次疊代就可以了，因此可得通常稱這類方法為預(yù)估－校正方法。它首先根據(jù)歐拉公式計(jì)算出xk+1的預(yù)估值xk+1(0)，然后再對(duì)它進(jìn)行校正，以得到更準(zhǔn)確的近似值xk+1(1)。3.1.4龍格－庫塔法根據(jù)泰勒級(jí)數(shù)將式（3-1）在tk+1=tk+h時(shí)刻的解xk+1=x(tk+h)在tk附近展開，有

（3-5）可以看出，提高截?cái)嗾`差的階次，便可提高其精度，但是由于計(jì)算各階導(dǎo)數(shù)相當(dāng)麻煩，所以直接采用泰勒級(jí)數(shù)公式是不適用的，為了解決提高精度問題，龍格和庫塔兩人先后提出了間接使用泰勒級(jí)數(shù)公式的方法，即用函數(shù)值f(t,x)的線性組合來代替f(t,x)的導(dǎo)數(shù)，然后按泰勒公式確定其中的系數(shù),這樣既能避免計(jì)算f(t,x)的導(dǎo)數(shù)，又可以提高數(shù)值計(jì)算精度，其方法如下。因故式（3-5）可寫成（3-6）為了避免計(jì)算式（3-6）中的各階導(dǎo)數(shù)項(xiàng)，可令xk+1由以下多項(xiàng)式表示。（3-7）式中am為待定因子，v為使用f函數(shù)值的個(gè)數(shù)，km滿足下列方程（3-8）即：將式（3-7）展開成h的冪級(jí)數(shù)并與微分方程式（3-1）精確解式（3-6）逐項(xiàng)比較,便可求得式（3-7）和式（3-8）中的系數(shù)am,bmj和cm等?，F(xiàn)以v=2為例，來說明這些參數(shù)的確定方法。設(shè)v=2，則有（3-9）

將k1和k2在同一點(diǎn)（tk,xk）上用二元函數(shù)展開為將k1和k2代入式（3-9）整理后可得

（3-10）將上式與式（3-6）逐項(xiàng)進(jìn)行比較，可得以下關(guān)系式若取則于是由式3-9可得(3-11)由于式（3-11）只取到泰勒級(jí)數(shù)展開式的h2項(xiàng)，故稱這種方法為兩階龍格－庫塔法，其截?cái)嗾`差為0(h3)。同理當(dāng)v=4時(shí)，仿照上述方法可得如下四階龍格-庫塔公式(3-12)通過上述龍格-庫塔法的介紹，可以把以上介紹的幾種數(shù)值積分法統(tǒng)一起來，它們都是基于在初值附近展開成泰勒級(jí)數(shù)的原理，所不同的是取泰勒級(jí)數(shù)多少項(xiàng)。歐拉公式僅取到h項(xiàng)，梯形法與二階龍格庫塔法相同，均取到h2項(xiàng)，四階龍格庫塔法取到h4項(xiàng)。從理論上講，取得的項(xiàng)數(shù)愈多，計(jì)算精度愈高，但計(jì)算量愈大，愈復(fù)雜，計(jì)算誤差也將增加，因此要適當(dāng)?shù)倪x擇。目前在數(shù)字仿真中，最常用的是四階龍格庫塔法，其截?cái)嗾`差為０(h5),已能滿足仿真精度的要求。3.1.5關(guān)于仿真數(shù)值積分法的幾點(diǎn)討論1．單步法和多步法解初值問題的數(shù)值解法的共同特點(diǎn)是步進(jìn)式，即從最初一點(diǎn)或幾點(diǎn)出發(fā)，每一步根據(jù)xk一點(diǎn)或前面幾點(diǎn)xk-1,xk-2,…來計(jì)算新的xk+1的值，這樣逐步推進(jìn)。當(dāng)從tk推進(jìn)到tk+1只需用tk時(shí)刻的數(shù)據(jù)時(shí)，稱為單步法，例如歐拉法和龍格庫塔法。相反，需要用到tk以及過去時(shí)刻tk-1

,tk-2

,…的數(shù)據(jù)時(shí)稱為多步法。線性多步法的一般形式是（3-13）多步法不能從t=0自啟動(dòng)，通常需要選用相同階次精度的單步法來啟動(dòng)，獲得所需前k步數(shù)據(jù)后，方可轉(zhuǎn)入相應(yīng)多步法，因多步法利用信息量大，因而比單步法更精確。2．顯式和隱式在計(jì)算xk+1時(shí)公式右端所用到的數(shù)據(jù)均已知時(shí)，稱為顯式算法。例如歐拉法、龍格庫－塔法和式（3-13）中β-1=0的情況。相反，在算式右端中隱含有未知量xk+1時(shí)，稱為隱式算法。例如梯形法、預(yù)估校正法和式（3-13）中β-1≠0的情況。顯式算法利用前幾步計(jì)算結(jié)果即可進(jìn)行遞推求解下步結(jié)果，因而易于計(jì)算。而隱式計(jì)算需要迭代法，先用另一同階次顯式公式估計(jì)出一個(gè)初值xk+1(0),并求得fk+1,然后再用隱式求得校正值xk+1(1),若未達(dá)到所需精度要求，則再次迭代求解，直到兩次迭代值xk+1

(i)和xk+1(i+1)

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

仿真誤差與數(shù)值計(jì)算方法、計(jì)算機(jī)的精度以及計(jì)算步長的選擇有關(guān)。當(dāng)計(jì)算方程和計(jì)算機(jī)確定以后，則僅與計(jì)算步長有關(guān)，所以在仿真中計(jì)算步長是一個(gè)重要的參數(shù)。仿真誤差一般有如下兩種：（1）截?cái)嗾`差――由于仿真模型僅是原系統(tǒng)模型的一種逼近，以及各種數(shù)值積分法的計(jì)算都是近似的算法。通常計(jì)算步長愈小，截?cái)嗾`差也愈小。（2）舍入誤差――由于計(jì)算機(jī)的精度有限（有限位數(shù)）所產(chǎn)生。通常計(jì)算步長愈小，計(jì)算次數(shù)愈多,舍入誤差愈大。對(duì)截?cái)嗾`差而言，計(jì)算步長愈小愈好，但太小不但會(huì)增加計(jì)算時(shí)間，而且由于舍入誤差的增加，不一定能達(dá)到提高精度的目的，甚至可能出現(xiàn)數(shù)值不穩(wěn)情況。顯然計(jì)算步長太大，不但精度不能滿足要求，而且計(jì)算步長超過該算法的判穩(wěn)條件時(shí)，也會(huì)出現(xiàn)不穩(wěn)定情況。由此可見，計(jì)算步長只能在某一范圍內(nèi)選擇，圖中的h0為最佳計(jì)算步長。一般控制系統(tǒng)的輸出動(dòng)態(tài)響應(yīng)在開始段變化較快，到最后變化將會(huì)很緩慢。這時(shí)，計(jì)算可以采用變步長的方法，即在開始階段步長取得小一些，在最后階段取得大一些，這樣即可以保證計(jì)算的精度，也可以加快計(jì)算的速度。對(duì)于一般工程計(jì)算，計(jì)算精度要求并不太高，故常用定步長的方法。作為經(jīng)驗(yàn)數(shù)據(jù)，當(dāng)采用四階龍格－庫塔法作數(shù)值積分計(jì)算時(shí)，取計(jì)算步長h=tr/10或ts/40式中tr－系統(tǒng)在階躍函數(shù)作用下的上升時(shí)間；ts－系統(tǒng)在階躍函數(shù)作用下的過渡過程時(shí)間。若系統(tǒng)有多個(gè)回路，則應(yīng)按反應(yīng)最快的回路考慮。3.3面向系統(tǒng)結(jié)構(gòu)圖的仿真自動(dòng)控制系統(tǒng)常常是由許多環(huán)節(jié)組成的。應(yīng)用上節(jié)介紹的數(shù)字仿真方法對(duì)系統(tǒng)分析和研究時(shí)，首先需要求出總的傳遞函數(shù)，再轉(zhuǎn)化為狀態(tài)空間表達(dá)式的形式，然后對(duì)其求解。當(dāng)改變系統(tǒng)某一環(huán)節(jié)的參數(shù)時(shí)，尤其是要改變小閉環(huán)中某一環(huán)節(jié)的參數(shù)時(shí)，以上整個(gè)過程又需重新計(jì)算，這對(duì)研究對(duì)象參數(shù)變化對(duì)整個(gè)控制系統(tǒng)的影響是十分不便的。為了克服這些缺點(diǎn)，同時(shí)大多數(shù)從事自動(dòng)化工作的科技人員更習(xí)慣于用結(jié)構(gòu)圖的形式來分析和研究控制系統(tǒng)，為此產(chǎn)生了面向結(jié)構(gòu)圖的仿真方法。該方法只需將各環(huán)節(jié)的參數(shù)及各環(huán)節(jié)間的連接方式輸入計(jì)算機(jī)，仿真程序就能自動(dòng)求出閉環(huán)系統(tǒng)的狀態(tài)空間表達(dá)式。本節(jié)主要介紹由典型環(huán)節(jié)參數(shù)和連接關(guān)系構(gòu)成閉環(huán)系統(tǒng)的狀態(tài)方程的方法，而動(dòng)態(tài)響應(yīng)的計(jì)算，仍采用四階龍格-庫塔法。這種方法與上節(jié)介紹的方法相比，有以下幾個(gè)主要優(yōu)點(diǎn)：）便于研究各環(huán)節(jié)參數(shù)對(duì)系統(tǒng)的影響；）可以得到每個(gè)環(huán)節(jié)的動(dòng)態(tài)響應(yīng)；）可對(duì)多入多出系統(tǒng)進(jìn)行仿真。下面具體介紹面向結(jié)構(gòu)圖的仿真方法。3.3.1典型環(huán)節(jié)的確定一個(gè)控制系統(tǒng)可能由各種各樣的環(huán)節(jié)所組成，但比較常見環(huán)節(jié)有：1）比例環(huán)節(jié)：G(s)=k2）積分環(huán)節(jié)：3）比例積分環(huán)節(jié)：4）慣性環(huán)節(jié)：5）超前滯后環(huán)節(jié)：6）二階振蕩環(huán)節(jié)：

為了編制比較簡單而且通用的仿真程序，必須恰當(dāng)選擇仿真環(huán)節(jié)。在這里選用圖3-7所示的典型環(huán)節(jié)作為仿真環(huán)節(jié)，即式中u為典型環(huán)節(jié)的輸入，x為典型環(huán)節(jié)的輸出。利用這個(gè)典型環(huán)節(jié)，只要改變a,b,c和d參數(shù)的值，便可分別表示以上所述的各一階環(huán)節(jié)，至于二階振蕩環(huán)節(jié)，則可用兩個(gè)一階環(huán)節(jié)等效連接得到，如圖3-8所示。圖3-7典型環(huán)節(jié)圖3-8二階振蕩環(huán)節(jié)的等效結(jié)構(gòu)圖3.3.2連接矩陣一個(gè)控制系統(tǒng)用典型環(huán)節(jié)來描述時(shí)，必須用連接矩陣把各個(gè)典型環(huán)節(jié)連接起來。所謂連接矩陣，就是用矩陣的形式表示各個(gè)典型環(huán)節(jié)之間的關(guān)系。下面介紹連接矩陣的建立方法,假設(shè)多輸入多輸出系統(tǒng)的結(jié)構(gòu)圖如圖3-9所示。同理，三階及三階以上的環(huán)節(jié)也完全可以用若干個(gè)一階環(huán)節(jié)等效連接得到。由此可見，任何一個(gè)復(fù)雜的控制系統(tǒng)都可以用若干個(gè)典型環(huán)節(jié)來組成。圖中帶數(shù)字的方框表示典型環(huán)節(jié)，表示比例系數(shù)。IIIIIIIVV由圖可得各環(huán)節(jié)的輸入與各環(huán)節(jié)的輸出間的關(guān)系，以及系統(tǒng)的輸出與環(huán)節(jié)的輸出間關(guān)系分別為和寫成矩陣形式和或?qū)懗?3-20)

定義式中的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)系。一般也將系統(tǒng)中各典型環(huán)節(jié)的系數(shù)寫成如下矩陣的形式（假設(shè)系統(tǒng)由n個(gè)典型環(huán)節(jié)組成）(3-21)3.3.3確定系統(tǒng)的狀態(tài)方程典型環(huán)節(jié)和連接矩陣確定后，便可求得系統(tǒng)的狀態(tài)空間表達(dá)式，推導(dǎo)過程如下。假設(shè)系統(tǒng)由n個(gè)典型環(huán)節(jié)組成，則根據(jù)典型環(huán)節(jié)的傳遞函數(shù)有

（i=1,2,…,ｎ）即

（i=1,2,…,ｎ）寫成矩陣形式：

（3-22）

式中：將式（3-20）中的上式進(jìn)行拉氏變換后代入式（3-22）中可得對(duì)上式兩邊取拉氏反變換得（3-23）若參考輸入向量r=[r1r2…rm]T中的r1,r2,…,rm均為階躍函數(shù)，則上式可簡化為（3-24）令則式（3-24）可寫成若H的逆存在，則有再令可得（3-25）上式即為閉環(huán)系統(tǒng)的狀態(tài)方程，它是一個(gè)典型的狀態(tài)方程，利用前面介紹的求解方法可方便地求出各典型環(huán)節(jié)的輸出響應(yīng)，最后根據(jù)式（3-20）中的第二式便可求出系統(tǒng)的輸出響應(yīng)。在建立系統(tǒng)的各典型環(huán)節(jié)時(shí)應(yīng)注意以下兩點(diǎn)：（1）為保證H的逆H-1存在，應(yīng)嚴(yán)格按照的原則，確定每個(gè)典型環(huán)節(jié)。即避免以純比例、純微分環(huán)節(jié)作為典型環(huán)節(jié)。（2）在輸入向量不全為階躍函數(shù)的情況下，只要在確定典型環(huán)節(jié)時(shí)，注意使含有微分項(xiàng)系數(shù)（即）的環(huán)節(jié)不直接與參考輸入連接，也可避免式3-23中出現(xiàn)r的導(dǎo)數(shù)。3.3.4面向結(jié)構(gòu)圖的數(shù)字仿真程序面向結(jié)構(gòu)圖的數(shù)字仿真程序框圖如圖３-10所示，其程序清單通過下例給出。給定典型環(huán)節(jié)參數(shù)和連接矩陣輸入仿真時(shí)間Tf

Tf和計(jì)算步長h開始給定輸入信號(hào)求H，H-1和Q陣值求A，B陣根據(jù)龍格-庫塔法求狀態(tài)方程的解計(jì)算系統(tǒng)輸出yT=TfyN輸出結(jié)果圖3-10面向結(jié)構(gòu)圖的數(shù)字仿真程序框圖例3-2假設(shè)某一系統(tǒng)由四個(gè)典型環(huán)節(jié)組成，如圖3-11所示。求輸出量y的動(dòng)態(tài)響應(yīng)。解由圖可得各環(huán)節(jié)的輸入與輸出以及系統(tǒng)的輸出與環(huán)節(jié)的輸出間關(guān)系為根據(jù)以上兩式和各典型環(huán)節(jié)的系數(shù)值，可得如下連接矩陣和系數(shù)矩陣仿真程序如下%Example3.2r=10;P=[0.110.51;0110;2120;101100];W=[000-1;1000;0100;0010];W0=[1;0;0;0];Wc=[0001];Tf=input('仿真時(shí)間Tf=');h=input('計(jì)算步長h=');A1=diag(P(:,1));B1=diag(P(:,2));C1=diag(P(:,3));D1=diag(P(:,4));H=B1-D1*W;

Q=C1*W-A1;A=inv(H)*Q;B=inv(H)*C1*W0;x=[zeros(length(A),1)];y=[z