第七章 系統(tǒng)時間響應(yīng)及其仿真

第七章系統(tǒng)時間響應(yīng)及其仿真第1頁，課件共31頁，創(chuàng)作于2023年2月第七章系統(tǒng)時間響應(yīng)及其仿真7.1仿真算法對系統(tǒng)的時間響應(yīng)進行動態(tài)仿真，采用什么樣的仿真算法是一個至關(guān)重要的問題。對連續(xù)時間系統(tǒng)進行數(shù)字動態(tài)仿真，主要是兩種方法：基于數(shù)值積分的仿真方法;基于離散相似法的仿真方法。由于后者涉及到離散控制系統(tǒng)理論，因此本節(jié)重點介紹基于數(shù)字積分的連續(xù)系統(tǒng)仿真方法。第2頁，課件共31頁，創(chuàng)作于2023年2月7.1仿真算法7.1.1數(shù)值求解的基本概念設(shè)微分方程為則求解方程中函數(shù)y(t)問題，就是已知初值的常微分方程求解問題。所謂數(shù)值求解就是要在時間區(qū)間[a,b]中取若干離散點

求出微分方程在這些時刻的近似值常微分方程數(shù)值求解的基本方法是數(shù)值積分法。第3頁，課件共31頁，創(chuàng)作于2023年2月7.1仿真算法7.1.2數(shù)值積分的基本原理積分區(qū)間的劃分將區(qū)間[a,b]分成N個小區(qū)間,時間間隔h（)也稱為積分步長,在第k個間隔t=[tk,tk+1]內(nèi)積分：則可用yk（k=0,1,…N）作為解y(t)的近似值，如圖所示。abtky0ykyt數(shù)值積分圖解tk+1第4頁，課件共31頁，創(chuàng)作于2023年2月7.1仿真算法7.1.2數(shù)值積分的基本原理數(shù)值積分的展開式為避免(2)式中的積分項，將y在tk

，以h為增量展開成Taylor級數(shù)：式(3)是一個遞推公式。積分值與實際微分方程解的誤差取決于步長h和計算所用的階數(shù)，它是數(shù)值積分的基礎(chǔ)。第5頁，課件共31頁，創(chuàng)作于2023年2月7.1仿真算法7.1.2數(shù)值積分的基本原理有關(guān)概念單步法和多步法單步法指計算yk+1值只需利用tk時刻的信息，也稱為自啟動算法；多步法在計算yk+1值時，則需利用tk，tk-1，…時刻的信息。顯示法和隱式法顯示法在計算yk+1時所需數(shù)據(jù)均已算出；隱式法在計算yk+1時需用到tk+1時刻的數(shù)據(jù)，該算法必須借助予估公式。定步長和變步長定步長為積分步長在仿真運行過程中始終不變；變步長指在仿真運行過程中自動修改步長。第6頁，課件共31頁，創(chuàng)作于2023年2月7.1仿真算法7.1.3數(shù)值積分的幾個算法歐拉算法在(3)式中取前兩項：可得歐拉算法：t0t1t2t3hy(t)y0y1yt歐拉近似解歐拉法【說明】歐拉法是用一條過各點的切線取代曲線來逼近精確解。該算法簡單，計算量小，但精度較低。第7頁，課件共31頁，創(chuàng)作于2023年2月7.1仿真算法7.1.3數(shù)值積分的幾個算法梯度法梯度法是歐拉法的改進。與歐拉法相比，梯度法是用兩個點(tm，ym)、(tm+1,ym+1)的斜率的平均值來確定下一點的y值。由于上式計算時需要用到y(tǒng)m+1的值，而ym+1不能預(yù)先知道，故梯度法需要和歐拉法結(jié)合使用,即用歐拉法對ym+1

進行予估，再由梯度法計算ym+1第8頁，課件共31頁，創(chuàng)作于2023年2月7.1仿真算法7.1.3數(shù)值積分的幾個算法龍格-庫塔法龍格-庫塔法的基本思想歐拉算法的精度較低，主要是其微分方程解y

的Taylor展開式所取的項數(shù)太少。顯然為了提高計算精度，應(yīng)當(dāng)取泰勒公式（3）更高階項。雖然增加高階項可提高計算精度，但也同時帶來了需要計算高階導(dǎo)數(shù)的困難。龍格-庫塔法的關(guān)鍵是利用低階導(dǎo)數(shù)構(gòu)成的曲線去擬合含有高階導(dǎo)數(shù)的曲線，從而避免了計算高階導(dǎo)數(shù)的問題。第9頁，課件共31頁，創(chuàng)作于2023年2月7.1仿真算法7.1.3數(shù)值積分的幾個算法龍格-庫塔法二階龍格-庫塔（RK）法取(3)式的前三項，則有設(shè)原微分方程(1)式解具有以下形式:式中，a1，a2，b1，b2為待定系數(shù)。第10頁，課件共31頁，創(chuàng)作于2023年2月將(8)式中K2按二元函數(shù)展開成泰勒級數(shù)，并取前三項將K1，K2代入(8)式：比較(6-10)、(6-7)式：第11頁，課件共31頁，創(chuàng)作于2023年2月顯然由(11)式并不能唯一確定a1,a2,b1,b2，因為只有三個方程。因此對于同一種算法可以有不同的表現(xiàn)形式?！菊f明】由于該算法只取到泰勒展開式的二階導(dǎo)數(shù)項，所以稱為二階龍格-庫塔法。但由(8)~(12)式可知，算法并沒有用y的二階導(dǎo)數(shù)。若設(shè)a1=a2，則即二階RK法公式為第12頁，課件共31頁，創(chuàng)作于2023年2月7.1仿真算法7.1.3數(shù)值積分的幾個算法龍格-庫塔法龍格-庫塔（RK）法的一般形式

式中，

i為待定權(quán)系數(shù)，ai,bij為待定系數(shù)，r為使用Ki的個數(shù)(即級數(shù))，Ki為所取各點導(dǎo)數(shù)f的值。

Ki的個數(shù)與yk+1泰勒展開式所取的項數(shù)有關(guān)(即RK算法的階數(shù))，同時還與計算區(qū)間內(nèi)所取導(dǎo)數(shù)值的點數(shù)有關(guān)。第13頁，課件共31頁，創(chuàng)作于2023年2月7.1仿真算法7.1.3數(shù)值積分的幾個算法龍格-庫塔法四階RK公式四階RK公式用到了y的泰勒展開式的四階導(dǎo)數(shù)。在RK算法的一般公式（13）中，取r=4可得：由于(14)式在同級的RK算法中，計算精度較高，計算量較少，而在系統(tǒng)仿真的數(shù)值積分中應(yīng)用十分廣泛。稱之為四階四級RK公式。第14頁，課件共31頁，創(chuàng)作于2023年2月7.1仿真算法7.1.3數(shù)值積分的幾個算法Gear算法“病態(tài)”常微分方程(剛性方程)的系數(shù)矩陣A的特征值具有如下特征：則稱為“病態(tài)”方程。第15頁，課件共31頁，創(chuàng)作于2023年2月7.1仿真算法7.1.3數(shù)值積分的幾個算法Gear算法控制系統(tǒng)仿真中的“病態(tài)”問題病態(tài)系統(tǒng)中絕對值最大的特征值對應(yīng)于系統(tǒng)動態(tài)性能解中瞬態(tài)分量衰減最快的部分，它反映了系統(tǒng)的動態(tài)響應(yīng)和系統(tǒng)的反應(yīng)靈敏度。一般與系統(tǒng)中具有最小時間常數(shù)Tmin的環(huán)節(jié)有關(guān)，要求計算步長h取得很小。病態(tài)系統(tǒng)中絕對值最小的特征值對應(yīng)于系統(tǒng)動態(tài)性能解中瞬態(tài)分量衰減最慢的部分,它決定了整個系統(tǒng)的動態(tài)過渡過程時間的長短。一般與系統(tǒng)中具有最大時間常數(shù)Tmax的環(huán)節(jié)有關(guān)，要求計算步長h取得很大。對于病態(tài)問題的仿真需要尋求更加合理的算法，以解決病態(tài)系統(tǒng)帶來的選取計算步長與計算精度、計算時間之間的矛盾。第16頁，課件共31頁，創(chuàng)作于2023年2月7.1仿真算法7.1.3數(shù)值積分的幾個算法Gear算法Gear算法Gear算法適用于病態(tài)系統(tǒng)的仿真，該算法類似于四階RK算法第17頁，課件共31頁，創(chuàng)作于2023年2月7.1仿真算法7.1.4數(shù)值積分方法的選擇在選擇積分方法時應(yīng)考慮以下幾個問題。計算精度數(shù)值積分方法所得到的離散數(shù)值解只是精確解的近似，其誤差來自兩個方面，即舍入誤差和局部截斷誤差。舍入誤差：由計算機字長有限而造成的計算時的舍入誤差，它隨計算次數(shù)的增加而增加。因此舍入誤差與計算步長

h

成反比。局部截斷誤差：由積分方法和階次的限制而引起的誤差。這種誤差與h成正比。截斷誤差舍入誤差總誤差eh

誤差與積分步長

顯然選擇一個合適的積分步長可使總誤差達到最小。第18頁，課件共31頁，創(chuàng)作于2023年2月7.1.4數(shù)值積分方法的選擇積分步長的選擇和控制積分步長的選擇原則在保證數(shù)值積分穩(wěn)定性和精度的前提下，盡可能選則較大的積分步長，以減少仿真計算次數(shù)和仿真時間。固定步長與變步長固定步長：在整個仿真計算過程中，積分步長h始終不變。其算法簡單，但很難保證步長最優(yōu)。此外，h還應(yīng)與模型的信號響應(yīng)情況有關(guān)，例如在穩(wěn)態(tài)時，可取較大的步長，見上圖。變步長：在仿真計算過程中根據(jù)計算誤差的大小來改變步長。其目的是在保證一定計算精度的前提下，盡可能選擇較大步長。第19頁，課件共31頁，創(chuàng)作于2023年2月7.2系統(tǒng)仿真的MATLAB函數(shù)7.2.1數(shù)值積分方法的MATLAB函數(shù)對于用數(shù)值方法求解常系數(shù)微分方程(OrdinaryDifferentialEquation,簡寫為ODE)或微分方程組，MATLAB提供了七種解函數(shù)，最常用的是ODE45（四階RK算法，單步、變步長，用五階RK算法估算局部截斷誤差），其調(diào)用格式為：[T,Y]=ode45(‘f’,tspan,y0)【說明1】‘f’為常微分方程（組）或系統(tǒng)模型的文件名；

tspan=[t0,tfinal]即積分時間初值和終值；

y0是積分初值；

T為計算時間點的時間向量；

Y為相應(yīng)的微分方程解數(shù)據(jù)向量或矩陣。第20頁，課件共31頁，創(chuàng)作于2023年2月7.2系統(tǒng)仿真的MATLAB函數(shù)7.2.1數(shù)值積分方法的MATLAB函數(shù)【說明2】對于剛性微分方程（特征值數(shù)值相差較大），可用ode15s，其調(diào)用格式與ode45相同。ode函數(shù)只能用于求解一階微分方程或一階微分方程組。若系統(tǒng)的數(shù)學(xué)模型為高階微分方程，則應(yīng)將高階微分方程轉(zhuǎn)化成一階微分方程組。因此在用MATLAB的ode函數(shù)求解微分方程時，應(yīng)首先建立描述系統(tǒng)模型的一階微分方程（組）函數(shù)‘f’。第21頁，課件共31頁，創(chuàng)作于2023年2月【例2】已知二階微分方程求時間區(qū)間t=[0，20]微分方程的解。解：（1）將微分方程表示為一階微分方程組【說明】這種描述系統(tǒng)微分方程的函數(shù)與ODE函數(shù)配套使用，其格式是固定的。dy為2*1數(shù)組，其維數(shù)等于微分方程的階數(shù)。（2）建立描述系統(tǒng)微分方程的m-函數(shù)文件vdp.mfunctiondy=vdp(t,y)dy=zeros(2,1);%生成2行1列的零陣dy(1)=y(2);%dy(2)=(1-y(1)^2)*y(2)-y(1);%第22頁，課件共31頁，創(chuàng)作于2023年2月（3）編寫MATLAB主程序[T,Y]=ode45('VDPd',[020],[0,1]);%調(diào)用ode45產(chǎn)生離散點時間向量和解向量plot(T,Y(:,1),'r-',T,Y(:,2),'b:')title('Solution')xlabel('times'),ylabel('PositionY')legend('y1','y2')運行結(jié)果如右圖所示。其中y1(紅線)為微分方程的解。第23頁，課件共31頁，創(chuàng)作于2023年2月7.2系統(tǒng)仿真的MATLAB函數(shù)7.2.2時間響應(yīng)仿真的MATLAB函數(shù)對于線性時不變系統(tǒng)，MATLAB直接提供了在各種輸入作用下的時間響應(yīng)函數(shù)。階躍響應(yīng)仿真函數(shù)（STEP）基本調(diào)用格式對于LTI連續(xù)(或離散)時間系統(tǒng)，以下調(diào)用格式可用于繪制系統(tǒng)單位階躍響應(yīng)曲線。step(sys)step(sys,Tfinal)其中，sys為系統(tǒng)模型（傳遞函數(shù)）；Tfinal為仿真終止時間，若省略則由系統(tǒng)默認(rèn)。第24頁，課件共31頁，創(chuàng)作于2023年2月【例4】已知系統(tǒng)模型，求其單位階躍響應(yīng)。sys=tf([1,-1],[1,1,5])subplot(1,2,1),step(sys,20)subplot(1,2,2),step(sys)建立系統(tǒng)模型指定階躍響應(yīng)時間不指定階躍響應(yīng)時間第25頁，課件共31頁，創(chuàng)作于2023年2月7.2.2時間響應(yīng)仿真的MATLAB函數(shù)階躍響應(yīng)仿真函數(shù)（STEP）多系統(tǒng)階躍響應(yīng)調(diào)用格式在同一幅圖中繪制多個系統(tǒng)的單位階躍響應(yīng)曲線，可用以下調(diào)用格式：這種調(diào)用格式，還可定義每個系統(tǒng)響應(yīng)曲線的顏色、線型和標(biāo)志，例如返回仿真輸出的調(diào)用格式其中，Y為輸出響應(yīng)，T為仿真時間向量。這種調(diào)用格式不繪制仿真曲線圖。step(sys1,sys2,…)step(sys1,'r',sys2,'y--',sys3,'gx')

[Y,T]=step(sys)第26頁，課件共31頁，創(chuàng)作于2023年2月7.2系統(tǒng)仿真的MATLAB函數(shù)7.2.2時間響應(yīng)仿真的MATLAB函數(shù)脈沖響應(yīng)仿真函數(shù)（IMPULSE）IMPULSE函數(shù)用來計算LTI系統(tǒng)的單位脈沖響應(yīng)。其調(diào)用格式與STEP函數(shù)相同。impulse(sys)impulse(sys,Tfinal)impulse(sys1,sys2,…)[Y,T]=impulse(sys)第27頁，課件共31頁，創(chuàng)作于2023年2月7.2.2時間響應(yīng)仿真的MATLAB函數(shù)信號發(fā)生器和任意輸入響應(yīng)函數(shù)MATLAB也可計算LTI系統(tǒng)在任意輸入作用下的時間響應(yīng)。信號發(fā)生器函數(shù)GENSIGGENSIG可為系統(tǒng)時間響應(yīng)產(chǎn)生周期輸入信號，其調(diào)用格式為：其中,Type為信號類型：‘sin’—正弦波

‘square’—方波

'pulse'—周期脈沖波Tau為信號周期；U為信號值向量;T為與U對應(yīng)的時間向量;Tf為信號的時間區(qū)間;Ts為采樣周期。[U,T]=gensig(Type,Tau)[U,T]=gensig(Type,Tau,Tf,Ts)第28頁，課件共31頁，創(chuàng)作于202