控制系統(tǒng)數(shù)字仿真課件

第4章 控制系統(tǒng)數(shù)字仿真數(shù)字仿真就是采用數(shù)學(xué)模型，在數(shù)字計算機(jī)上借助數(shù)值解法所進(jìn)行的仿真實(shí)驗。所謂數(shù)值解法，就是尋求y(t)在a,b區(qū)間內(nèi)的一系列離散節(jié)點(diǎn)t1t2tmtm+1上的近似值y1, y2, ，ym, ym+1, ，即求取ym+1 y(tm+1)。相鄰節(jié)點(diǎn)的間距h=tm+1-tm稱為步長，這里假定h為定值，即tm=t0+mh，m=0, 1, 2,。本章主要講述數(shù)字仿真的基本理論和方法。4.1 連續(xù)系統(tǒng)數(shù)值積分方法系統(tǒng)的動態(tài)特性通常用一階微分方程組來描述，也即狀態(tài)空間表達(dá)式。一般來說，只有極少數(shù)的微分方程能用到初等方法得其解析解（或用解析的方法得到精確解），多數(shù)只能用近似數(shù)值求解。利用計算機(jī)求微分方程主要使用數(shù)值積分法，它是系統(tǒng)仿真的一個重要方法。在這里，我們主要研究一階微分方程的形式，如： （1）， 求y(t)解： 當(dāng)t=tm+1，t0=tm時 （2）數(shù)值積分法時在已知初值的情況下，對f(t, y)進(jìn)行近似積分，從而對y(t)進(jìn)行數(shù)值求解的方法。下面介紹幾種在數(shù)字仿真常用的數(shù)字積分法。4.1.1 歐拉法歐拉法又稱為折線法，是最簡單，也是最早的一種數(shù)值計算方法。對于式（2），如果積分間隔h=tm+1 - tm取得足夠小，使得在tm與tm+1之間的f(t, y)可近似看做常數(shù)f(tm, ym)。這樣式（2）可化為：即 （3）（3）式即為歐拉公式。歐拉公式的幾何解釋：對于微分方程（1）的解y(t)看作是一條曲線，在任一步長內(nèi)，用一段直線代替函數(shù)y(t)的曲線，此直線段得斜率等于該函數(shù)在步長起點(diǎn)的斜率?；谏鲜龅膸缀谓忉專覀儚某跏键c(diǎn) (t0, y0)出發(fā)向前推進(jìn)(t1, y1)點(diǎn)，(t2, y2)點(diǎn)，t0t1t2tm+1t0y0y1y2ym+1y圖中陰影部分即為誤差。歐拉法的特點(diǎn)是：計算簡單，但精度較低。例：用歐拉法求解初值問題（），h=0.1。解：因為則，該題解為：，將準(zhǔn)確解y(tm)與近似解ym一起放入下表，可得：tm ym y(tm)tm ym y(tm) 0.1 1.1000 1.09540.2 1.1918 1.18320.3 1.2774 1.26490.4 1.3582 1.34160.5 1.4351 1.4142 0.6 1.5090 1.48320.7 1.5803 1.59420.8 1.6498 1.61250.9 1.7178 1.67331.0 1.7848 1.7321由此表可以看出歐拉公式的精度很差。4.1.2 后退的歐拉法若用f(tm+1,ym+1)來代替f(tm,ym)，則（3）式可變?yōu)椋?（4）則（4）式稱為后退的歐拉公式。后退的歐拉公式是隱式的（因為（4）式右邊的ym+1是未知的），此時通常需要用迭代法求解，即： （5）t0t1t2tm+1t0y0y1y2ym+1y后退的歐拉公式的幾何解釋：在任一步長內(nèi)，用一段直線代替函數(shù)y(t)的曲線，此直線段得斜率等于該函數(shù)在步長終點(diǎn)的斜率。例：用后退的歐拉法求解初值問題（），h=0.1。解：后退的歐拉公式和歐拉公式的精度相同，都是一階精度。4.1.3 梯形法比較歐拉公式和后退的歐拉公式可知，如果對這兩種方法進(jìn)行算術(shù)平均，即可大大消除主要誤差，從而獲得更大的精度，這種方法通常稱為梯形法，其計算公式為： （6）同后退的歐拉公式一樣，梯形公式也是隱式的。此時通常采用歐拉公式先預(yù)報一個，再將預(yù)報的代入（6）式進(jìn)行校正，求出ym+1。梯形法的迭代公式為： （7）（7）式又被稱為預(yù)估-校正公式。顯然，梯形法要比歐拉法和后退的歐拉法精度更高，但計算量比歐拉法大。例：用梯形法求解初值問題（），h=0.1。解：tm ym y(tm)tm ym y(tm) 0.1 1.0959 1.09540.2 1.1841 1.18320.3 1.2662 1.26490.4 1.3434 1.34160.5 1.4164 1.4142 0.6 1.4860 1.48320.7 1.5525 1.59420.8 1.6153 1.61250.9 1.6782 1.67331.0 1.7379 1.73214.1.4 龍格庫塔法則：我們把稱作tm,tm+1上的平均斜率，由此可見，只要給平均斜率提供一種算法，便可導(dǎo)出一種計算公式。當(dāng)=0時，平均斜率為f(tm,ym)，為歐拉公式；當(dāng)=1時，平均斜率為f(tm+1,ym+1)，為后退的歐拉公式；若采用tm和tm+1兩點(diǎn)斜率值的算術(shù)平均作為平均斜率，則平均斜率為，此時為梯形公式（精度更高）。這給我們啟示，如果設(shè)法在tm,tm+1區(qū)間內(nèi)多預(yù)測幾個點(diǎn)的斜率值，然后將它們進(jìn)行加權(quán)平均作為平均斜率，則可構(gòu)造出精度更高的計算公式，這就是龍格庫塔的基本思想。（1）二階龍格庫塔法隨意考查區(qū)間tm,tm+1區(qū)間內(nèi)一點(diǎn)：，用tm和tm+p兩個點(diǎn)的斜率值K1，K2線性組合得到平均斜率，則 （8）因為：所以：（8）又因為： （9）（根據(jù)泰勒公式）比較（8）式跟（9）式，可得： （10）滿足（10）式的一族公式（8）統(tǒng)稱為二階龍格庫塔公式。當(dāng)p=1時，此時，此時：即梯形公式是二階龍格庫塔的一個特例。另外給出一個常用的二階龍格庫塔公式：（）（2）三階龍格庫塔法為了進(jìn)一步提高精度，取tm,tm+1區(qū)間內(nèi)三個點(diǎn)，即tm，tm+p= tm +ph，tm+q= tm +qh，其中，0pq1，用這三點(diǎn)的斜率值K1，K2，K3的線性組合得到平均斜率，此時： （11）經(jīng)過推導(dǎo)可得： （12）滿足（12）式的一族公式（11）統(tǒng)稱為三階龍格庫塔公式。另外給出一個常用的三階龍格庫塔公式：（3）四階龍格庫塔法同理，取 tm,tm+1區(qū)間內(nèi)四個點(diǎn)的斜率值進(jìn)行線性組合來得到平均斜率，即為四階龍格庫塔法。因其計算很復(fù)雜，直接給出最常用的四階龍格庫塔公式： （13）例：用四階龍格庫塔法求解初值問題（），h=0.2。解：tm ym(1) ym(2)y(tm)0.10.20.30.40.50.60.70.80.91.01.09591.18411.26621.34341.41641.48601.55251.61531.67821.73791.18321.34171.48331.61251.73211.18321.34161.48321.61251.7321四階龍格庫塔法的計算量比梯形法大將近一倍，但由于放大了步長，使得計算量幾乎相同。由此可以看出選擇算法的重要性。關(guān)于數(shù)值積分法的幾點(diǎn)討論：（1）單步法和多步法以上介紹的數(shù)值積分法都有一個特點(diǎn)，每次計算只需用到前一次的計算結(jié)果，而不需要更前面的計算結(jié)果，我們稱這類方法為單步法。優(yōu)點(diǎn)：占用內(nèi)存少能自啟動（只要知道初值就可自行進(jìn)行計算）容易實(shí)現(xiàn)變步長計算 多步法比單步法更精確（因為多步法利用的信息量更大）（2）顯式和隱式 顯式算法可直接利用遞推公式求解下一步結(jié)果，而隱式算法需要采用迭代法。隱式算法計算過程復(fù)雜、計算速度慢，但通常精度較高。（3）數(shù)值積分的穩(wěn)定性例：， 0t1.5其精確解為：取步長h=0.1，則歐拉法：y1.5=-1.0922510-4四階龍格庫塔法：y1.5=3.95730101精確解為：y(1.5)=9.5417310-21可見此時數(shù)值積分的結(jié)果是不正確的，或是不穩(wěn)定的。原因：數(shù)值積分是一種近似積分，它在每一步的遞推中將引入誤差，所以若誤差的積累越積越大，計算將出現(xiàn)不穩(wěn)定?，F(xiàn)在我們以下面的微分方程為例，來分析上面所述的各種數(shù)值積分法的穩(wěn)定性， 0表示有正的連接關(guān)系，wij0表示有負(fù)的連接關(guān)系；wii0，表示各環(huán)節(jié)i有自反饋。采用連接矩陣表示各環(huán)節(jié)的連接關(guān)系，使得對復(fù)雜系統(tǒng)結(jié)構(gòu)的仿真變得簡明方便。4.3 面向系統(tǒng)結(jié)構(gòu)圖的數(shù)字仿真4.3.1典型閉環(huán)系統(tǒng)的數(shù)字仿真1.典型閉環(huán)系統(tǒng)的結(jié)構(gòu)形式典型閉環(huán)系統(tǒng)結(jié)構(gòu)如下圖所示：y - e r G(s) v 圖中，V為系統(tǒng)的反饋系數(shù)，設(shè)其為常數(shù)。設(shè)圖中開環(huán)傳函為： （14）2.系統(tǒng)仿真模型與求解思路由（14）式可得能控標(biāo)準(zhǔn)型為：式中，其中，A、C陣中的，為（14）式首一化后的分子分母多項式系數(shù)，即， 。令 u=e則 u=e=r-v y令，則 求解思路：（1）歐拉公式：（2）梯形公式：（3）四階龍格庫塔公式K1，K2，K3，K4為狀態(tài)變量X在t=tk時的四組斜率，由此可求出tk+1時刻的狀態(tài)變量Xk+1，從而求得tk+1時刻的輸出：yk+1=CXk+1。所以只要給出了狀態(tài)變量的初值，就可以由上面的遞推公式，求出各時刻的狀態(tài)變量Xk和輸出yk。系統(tǒng)結(jié)構(gòu)圖開環(huán)傳函開環(huán)狀態(tài)空間表達(dá)式閉環(huán)狀態(tài)空間表達(dá)式采用數(shù)值積分法（如四階龍格庫塔法）進(jìn)行仿真3.仿真程序框圖與實(shí)現(xiàn)（1）程序框圖（參見書P133頁圖4.9）（2）程序 輸入數(shù)據(jù)a=a0,a1,.,an;b=b0,b1,.,bm; X0=x0,x1,.,xn;V=V0;n=n0;T0=t0;Tf=tf;h=h0;r=r0;形成開閉環(huán)系數(shù)矩陣b=b/a(1);a=a/a(1);A=a(2:n+1);A=rot90(rot90(eye(n-1,n);-fliplr(A);B=zeros(1,n-1),1;m1=length(b);C=fliplr(b),zeros(1,n-m1);A1=A-B*V*C;X=X0;y=0;t=T0;運(yùn)算求解N=round(Tf -T0)/h);for i=1:N K1=A1*X+B*r; K2=A1*(X+h*K1/2)+B*r; K3=A1*(X+h*K2/2)+B*r; K4=A1*(X+h*K3)+B*r; X=X+h*(K1+2*K2+2*K3+K4)/6; y=y;C*X; t=t;t(i)+h;end輸出結(jié)果t,yplot(t,y)例題（書P135）程序如下：程序1：k=1;a=conv(1 0 0,conv(0.25 1,0.25 1);b=2*k, k;X0=0 0 0 0;V=1;n=4;T0=0;Tf=300;h=0.3;r=1;b=b/a(1);a=a/a(1);A=a(2:n+1);A=rot90(rot90(eye(n-1,n);-fliplr(A);B=zeros(1,n-1),1;m1=length(b);C=fliplr(b),zeros(1,n-m1)A1=A-B*V*CX=X0;y=0;t=T0;N=round(Tf-T0)/h);for i=1:N K1=A1*X+B*r; K2=A1*(X+h*K1/2)+B*r; K3=A1*(X+h*K2/2)+B*r; K4=A1*(X+h*K3)+B*r; X=X+h*(K1+2*K2+2*K3+K4)/6; y=y,C*X; t=t,t(i)+h;endplot(t,y)程序2：r=1;num1=2 1;den1=0.0625 0.5 1 0 0;num2=1;den2=1;num,den=feedback(num1,den1,num2,den2);A,B,C,D=tf2ss(num,den)Tf=100; % input(仿真時間Tf=);h=0.25; % input(計算步長h=);X=zeros(length(A),1);y=0;t=0;for i=1:Tf/h K1=A*X+B*r; K2=A*(X+h*K1/2)+B*r; K3=A*(X+h*K2/2)+B*r; K4=A*(X+h*K3)+B*r; X=X+h*(K1+2*K2+2*K3+K4)/6; y=y;C*X; t=t;t(i)+h;endplot(t,y)4.3.2復(fù)雜連接閉環(huán)系統(tǒng)的數(shù)字仿真基本思路：與實(shí)際系統(tǒng)的結(jié)構(gòu)圖相對應(yīng)，只需將各環(huán)節(jié)的參數(shù)及各環(huán)節(jié)間的連接方式輸入計算機(jī)，仿真程序就能自動求出閉環(huán)系統(tǒng)的狀態(tài)空間表達(dá)式，從而運(yùn)用數(shù)值積分法進(jìn)行求解。（1）典型環(huán)節(jié)比例環(huán)節(jié)： 積分環(huán)節(jié)： 比例積分環(huán)節(jié)： 一階慣性環(huán)節(jié)： 一階超前滯后環(huán)節(jié)：二階振蕩環(huán)節(jié)： 對于二階振蕩環(huán)節(jié)， y u y u- y u 可見，我們通常遇到的都是一階環(huán)節(jié)，即使是二階振蕩環(huán)節(jié)，也可用兩個一階環(huán)節(jié)等效連接得到，因此在本章中，我們認(rèn)為所遇到的環(huán)節(jié)均為一階環(huán)節(jié)。典型的一階環(huán)節(jié)（通用）yi ui 用此格式可描述所有環(huán)節(jié)的輸入輸出關(guān)系。設(shè)輸入向量為U=u1,u2,unT輸出向量為Y =y1,y2,ynT則系統(tǒng)中所有環(huán)節(jié)的輸入輸出關(guān)系可表示為：(A+BS)Y=(C+DS)U （15）其中，（2）仿真求解把系統(tǒng)中各環(huán)節(jié)描述出來后，還需要進(jìn)行數(shù)值積分求解，此時需要把各環(huán)節(jié)的相互關(guān)系表示出來，即用上節(jié)所講的連接絕陣。舉例說明：y r - - 1 2 3 3 2 即 U=WY+WOr （16）, 連接矩陣及模型參數(shù)確定之后，求取閉環(huán)系統(tǒng)的狀態(tài)空間表達(dá)式（仿真模型）。將（16）式代入（15）式，得：設(shè)Q=B - DW，且逆矩陣存在，H=CW A，則為不增加求解的復(fù)雜性，應(yīng)使DW0=0，則設(shè)，則求解一階微分方程組的矩陣形式，就可采用數(shù)值積分法進(jìn)行求解。應(yīng)注意兩點(diǎn)： 為保證Q-1的存在，應(yīng)按照bi0的原則確定每個典型環(huán)節(jié)，即避免以純比例、純微分環(huán)節(jié)作為典型環(huán)節(jié)。 為保證DWO=0，應(yīng)使含有微分項系數(shù)（即 di0）的環(huán)節(jié)不直接與外加輸入相連。（3）仿真程序框圖與實(shí)現(xiàn) 系統(tǒng)模型參數(shù)矩陣輸入處理一般將系統(tǒng)中各典型環(huán)節(jié)的系數(shù)整理成如下矩陣的形式（假設(shè)系統(tǒng)由n個典型環(huán)節(jié)組成）：，A=diag(P(:,1)或 a=a1,a2,anA=diag(a)然后由程序自動轉(zhuǎn)換生成A、B、C、D陣（diag函數(shù)）。 連接矩陣輸入方法系統(tǒng)連接矩陣W和輸入連接矩陣WO中大部分元素為零，非零元素較少，為加快輸入速度，通常采用只輸入非零元素的方式來構(gòu)成WIJ陣（由W和WO陣中的非零元素及對應(yīng)的環(huán)節(jié)號所構(gòu)成的矩陣），然后由程序自動轉(zhuǎn)換成W和WO。具體作法是建立非零元素矩陣WIJ陣（m3）i j wij 其中m為W和WO中非零元素的個數(shù)，i為被作用環(huán)節(jié)，j為作用環(huán)節(jié)（外加輸入用0來表示），wij為對應(yīng)的連接系數(shù)i j wij 程序框圖（參見書P141頁圖4.16） 程序?qū)崿F(xiàn)P=a1,b1,c1,d1; a2,b2,c2,d2; an,bn,cn,dnWIJ=i, i, wij ; ; Y0=y10, y20,., yn0 ;L=L0 ;n=n0 ;T0=t0 ;Tf =tf ;h=h0 ;r=r0 ;nout=nput;A=diag(P(:,1);B=diag(P(:,2);C=diag(P(:,3);D=diag(P(:,4);m=length(WIJ(:,1);W0=zeros(n,1);W=zeros(n,n);for k=1:m if WIJ(k,2)=0 W0(WIJ(k,1)=WIJ(k,3); else W(WIJ(k,1),WIJ(k,2)=WIJ(k,3); endendQ=B-D*W;Q1=inv(Q);H=C*W-A;A1=Q1*H;B1=Q1*C*W0;Y=Y0;y=Y(nout);t=T0;e=0;N=round(Tf-T0)/h*L);for i=1:Nfor j=1:L K1=A1*Y+B1*r; K2=A1*(Y+h*K1/2)+B1*r; K3=A1*(Y+h*K2/2)+B1*r; K4=A1*(Y+h*K3)+B1*r; Y=Y+h*(K1+2*K2+2*K3+K4)/6;en y=y;Y(nout); t=t;t(i)+h*L; e=e;r-Y(nout);endt, yplot(t,y,t,e)4.4 連續(xù)系統(tǒng)的離散相似法4.4.1連續(xù)系統(tǒng)狀態(tài)方程的離散化該狀態(tài)方程的解為：當(dāng)t=kT時， 當(dāng)t=KT時，設(shè)t=-kT，即= t+kT，則d = dt 若采用零階保持器，則 u(t+kT)= u(kT)，0t0? 返回調(diào)用|u|s? MATLAB程序：function y=saturation(u,s) if(abs(u)=s) if(u0) y=s; else y=-s; end else y=u; end0 -s s u y 2.死區(qū)非線性數(shù)學(xué)描述為： 其中，s為死區(qū)值。N N Y Y y=u-s y=u+s y=0 u0? 返回調(diào)用|u|s? MATLAB程序：function y=dead(u,s) if(abs(u)=s) if(u0) y=u-s; else y=u+s; end else y=0; end0 -s s u y 3.滯環(huán)非線性數(shù)學(xué)描述為： 其中，s為滯環(huán)寬度；yb表示非線性環(huán)節(jié)的前一次輸出值，其值需要判斷