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

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

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

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

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

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

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

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

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

這就是歐拉公式。以x(t0)=x0作為初始值，應用歐拉公式，就可以一步步求出每一時刻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)就求出來了。歐拉法的計算雖然比較簡單，但精度較低。圖3-1為歐拉法的幾何解釋。3.1.2梯形法

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

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

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

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

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

,tk-2

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

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

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

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

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

定義式中的W

，W0和Wc陣為連接矩陣，W反映了各典型環(huán)節(jié)輸入輸出間的連接關系，W0反映了系統(tǒng)的參考輸入與各環(huán)節(jié)輸入間的連接關系，Wc反映了系統(tǒng)的輸出與各環(huán)節(jié)輸出間的關系。一般也將系統(tǒng)中各典型環(huán)節(jié)的系數(shù)寫成如下矩陣的形式（假設系統(tǒng)由n個典型環(huán)節(jié)組成）(3-21)3.3.3確定系統(tǒng)的狀態(tài)方程典型環(huán)節(jié)和連接矩陣確定后，便可求得系統(tǒng)的狀態(tài)空間表達式，推導過程如下。假設系統(tǒng)由n個典型環(huán)節(jié)組成，則根據(jù)典型環(huán)節(jié)的傳遞函數(shù)有

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

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

（3-22）

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

Tf和計算步長h開始給定輸入信號求H，H-1和Q陣值求A，B陣根據(jù)龍格-庫塔法求狀態(tài)方程的解計算系統(tǒng)輸出yT=TfyN輸出結(jié)果圖3-10面向結(jié)構圖的數(shù)字仿真程序框圖例3-2假設某一系統(tǒng)由四個典型環(huán)節(jié)組成，如圖3-11所示。求輸出量y的動態(tài)響應。解由圖可得各環(huán)節(jié)的輸入與輸出以及系統(tǒng)的輸出與環(huán)節(jié)的輸出間關系為根據(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('仿真時間Tf=');h=input('計算步長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