控制系統(tǒng)仿真及MATLAB語言-連續(xù)系統(tǒng)的離散化方法.ppt

第四章 連續(xù)系統(tǒng)的離散化方法,4.1 常微分方程的數(shù)值解法,一. 數(shù)值求解的基本概念,設(shè)微分方程為,則求解方程中函數(shù)x(t)問題的常微分方程初值問題,所謂數(shù)值求解就是要在時間區(qū)間a, b中取若干離散點,求出微分方程在這些時刻的近似值,取前兩項近似：,這種方法的幾何意義就是把f(t,x)在區(qū)間tk,tk+1內(nèi)的曲邊面積用矩形面積近似代替。計算簡單，計算量小，而且可以自啟動。當h很小時，造成的誤差是允許的。該算法具有一階精度。,取k=0,1,2,N,從t0開始，逐點遞推求解t1時的y1， t2時的y2,直至tn時的yn，稱之為歐拉遞推公式。,矩形面積,1. 歐拉法,歐拉法的特點：導(dǎo)出簡單，幾何意義明顯，便于理解，能說明構(gòu)造數(shù)值解法一般計算公式的基本思想。通常用它來說明有關(guān)的基本概念。,例 設(shè)系統(tǒng)方程為,用Euler法求其數(shù)值解（取步長 , ),遞推公式為,則,已知方程的解析解為 精確解和解析解作比較：,誤差在 數(shù)量級，精度較差。,2. 龍格庫塔法,基本思想：取Taylor級數(shù)展開式前三項近似求解，并利用線性組合代替導(dǎo)數(shù)的求解。 既可避免計算高階導(dǎo)數(shù)，又可提高數(shù)值積分的精度，這就是Runge-Kutta法的基本思想。,2. 龍格庫塔法,r為精度階次，ai為待定系數(shù)，由精度確定；ki用下式表示,線性組合,等各階導(dǎo)數(shù)不易計算，用下式中ki的線性組合代替,1)當r=1時：,與Taloy展開式相比較，可得a1=1,則上式成為,歐拉遞推公式,2)當r=2時：,將 在點 展成Taylor級數(shù),與臺勞公式的二階展開近似公式相比，可得以下關(guān)系：,三個方程，四個未知數(shù)，解不唯一,各個系數(shù)的幾種取法見書上。,3) r=4時，四階龍格庫塔公式-最常用：,仿真中遇到的大多數(shù)工程實際問題，四階龍格庫塔法以能滿足精度要求，其截斷誤差o(h5) 與h5同數(shù)量級。該法可以自啟動。,4)、狀態(tài)空間四階龍格-庫塔遞推式 若單輸入單輸出系統(tǒng)的狀態(tài)空間表達式為：,在仿真中，對于n階系統(tǒng)，狀態(tài)方程可以寫成一階微分方程,根據(jù)四階龍格-庫塔公式,有,T=tk時刻的xi值,T=tk+h時刻的xi值,另,狀態(tài)方程的四階龍格-庫塔公式如下：,RK法的特點：,1 需要存儲的數(shù)據(jù)少，占用的存儲空間少； 2 只需知道初值，即可啟動遞推公式進行計 算，可自啟動； 3 容易實現(xiàn)變步長運算。 4 每積分一步需要計算多次右函數(shù)，計算量 大。,基于龍格庫塔法，MATLAB提供了求常微分方程數(shù)值解 的函數(shù)，一般調(diào)用格式為：,t, x=ode23(xfun, t0, tf ,x0),t, x=ode45(xfun, t0, tf ,x0),常微分方 程函數(shù)名,起始 時間,終止 時間,初始狀 態(tài)向量,輸入,輸出,4/5階龍格-庫塔算法,2/3階龍格-庫塔算法,3.常微分方程Matlab求解,解: 令 y1=x，y2=x,1、建立M-文件vdp.m如下： function dy=vdp(t,y) dy=zeros(2,1); dy(1)=y(2); dy(2)=2*(1-y(1)2)*y(2)-y(1);,2、取t0=0，tf=20，輸入命令： T,Y=ode45(vdp,0 10,1;1); plot(T,Y(:,1),-, T,Y(:,2),3、結(jié)果,解 1、建立m-文件rigid.m如下： function dy=rigid(t,y) dy=zeros(3,1); dy(1)=y(2)*y(3); dy(2)=-y(1)*y(3); dy(3)=-0.51*y(1)*y(2);,2、取t0=0，tf=12，輸入命令： T,Y=ode45(rigid,0 12,0 1 1); plot(T,Y(:,1),-,T,Y(:,2),*,T,Y(:,3),+),3、結(jié)果如圖,圖中，y1的圖形為實線，y2的圖形為“*”線，y3的圖形為“+”線.,4.2 數(shù)值算法的穩(wěn)定性及求解原則,1.數(shù)值算法的穩(wěn)定性,特征根在s平面的左半平面，系統(tǒng)穩(wěn)定。 (1)歐拉法： 穩(wěn)定： (2)梯形法： 恒穩(wěn),2.數(shù)值算法的選擇原則,Matlab提供了微分方程數(shù)值求解的一般方法,作為仿真算法的使用者,可不必考慮算法具體實現(xiàn),而應(yīng)關(guān)心各種方法在使用中會出現(xiàn)的問題,以及如何在仿真中恰當?shù)倪x用這些方法. 一般,選用數(shù)值算法從以下幾個方面考慮: (1)精度 受算法和h影響 截斷誤差+舍入誤差=累計誤差 (2)計算速度 受算法和h影響 算法簡單，速度就快些。 (3)穩(wěn)定性 受h影響，一般 h（2-3）,系統(tǒng)最小時間,4.3 數(shù)值算法中的“病態(tài)”問題,1 “病態(tài)”常微分方程,例：,其中,采用四階龍格庫塔法,h=0.01時，計算時間長,h=0.04時，誤差很大,當h0.05后，曲線發(fā)散振蕩，數(shù)值不穩(wěn)定，完全失去意義,系統(tǒng)矩陣的特征值差異較大,一般線性常微分方程組：,的系數(shù)矩陣A的特征值具有如下特征：,則稱為“病態(tài)”方程。,定義：,2 控制系統(tǒng)仿真中的“病態(tài)”問題,1 病態(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取得很小。,2 病態(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取得很大。,3 對于病態(tài)問題的仿真需要尋求更加合理的算法，以解決病態(tài)系統(tǒng)帶來的選取計算步長與計算精度，計算時間之間的矛盾。,3 “病態(tài)”系統(tǒng)的仿真方法,采用穩(wěn)定性好，計算精度高的數(shù)值算法，并且允許計算步長能根據(jù)系統(tǒng)性能動態(tài)變化的情況在一定范圍內(nèi)作相應(yīng)的變化，采用隱式吉爾法,該法已經(jīng)證明對病態(tài)方程求解過程是數(shù)值穩(wěn)定的。,隱式吉爾法從理論上十分適應(yīng)于病態(tài)系統(tǒng) ，但需要解決好以下問題,(1) 自啟動 r階多步算式無法自啟動，需要用單步法求出前r步值,(2) 預(yù)估迭代 迭代方法要求收斂性良好，否則在大步長時會造成數(shù) 值發(fā)散。,(3) 變步長 初始階段采用小步長，隨后可逐步放大步長。,對不同精度要求的系統(tǒng)仿真，要考慮變階次問題，即為減小每一步計算的截斷誤差，以提高精度，應(yīng)選用較高的階次，而當精度較低時，為減少工作量，則應(yīng)選取較低的階次。仿真時應(yīng)根據(jù)估計誤差 與給定的誤差精度相比較改變步長或階次來重新計算。,4.4 連續(xù)系統(tǒng)狀態(tài)方程的離散化,上章所述的連續(xù)系統(tǒng)數(shù)學(xué)模型的離散化，是通過數(shù)值積分法實現(xiàn)的，盡管面向結(jié)構(gòu)圖的仿真方法是按環(huán)節(jié)給定參數(shù)，但是在計算時還是按整個系統(tǒng)進行離散化，這就不便于引進非線性環(huán)節(jié)以進行非線性系統(tǒng)的仿真。在本節(jié)，將介紹連續(xù)系統(tǒng)離散模型的建立和仿真。,數(shù)值積分法疊代求解,h改變時，疊代過程重復(fù)求解，費時繁瑣 不能對非線性環(huán)節(jié)單獨考慮。,連續(xù)系統(tǒng)離散化,思想：用差分方程描述連續(xù)系統(tǒng)的狀態(tài)方程模型 （因為差分方程的主要特點就是方程中各變量由各相鄰時刻的變 化量制約，這相當于遞推方程）,1、連續(xù)系統(tǒng)的離散化,設(shè)連續(xù)系統(tǒng)狀態(tài)方程為,則由現(xiàn)代控制理論基礎(chǔ)知，狀態(tài)變量X(t)的解為,而 t = ( k + 1 )T 時，可表示為,當系統(tǒng)輸入u(t)給定時，可求出系統(tǒng)離散化狀態(tài)方程的解。一般， u(t)未知，通常采用兩種方法近似處理： (1)令u(kT+t)u(kT) (0tT) 相當在系統(tǒng)輸入端加一個采樣開關(guān)和零階保持器 X(k+1)T) = GX(kT) + Hu(kT) G = e A T，為t = T 時的狀態(tài)轉(zhuǎn)移矩陣,離散后的狀態(tài)空間表達式為：,、Matlab表示,已知連續(xù)系統(tǒng)狀態(tài)方程為,在采樣周期T下離散后的狀態(tài)空間表達可表示為：,在Matlab中，若已知連續(xù)系統(tǒng)狀態(tài)方程各陣模型參數(shù) (A、B、C、D) 以及采樣周期T，則語句： G，H = c2d (A,B,T) 返回的矩陣G 、H 就是所要求的( T ) 、m ( T ) 。 此外， Matlab還提供了功能更強的求取連續(xù)系統(tǒng)離散化矩陣函數(shù)c2dm(),他容許調(diào)用時選用離散化變換方式，并且得到的是標準的離散化狀態(tài)方程。 G,H,C,D=c2dm (A,B,C,D,T,選項),表 離散化變換方式選項,2.離散函數(shù)的連續(xù)化 在MATLAB中也提供了從離散化系統(tǒng)轉(zhuǎn)換為連續(xù)系統(tǒng)各系數(shù)矩陣求取的功能函數(shù)，其調(diào)用格式分別如下 A ,B=d2c(G ,H ,T) 或 A ,B ,C ,D=d2cm (G ,H ,C,D ,T ,選項) 其中選項同上。,例 對連續(xù)系統(tǒng)在采樣周期T=0.1時進行離散化。,解：利用以下Matlab命令，對系統(tǒng)進行離散化 k=6; Z=-3; P=-1；-2；-5； T=0.1; A,B,C,D=zp2ss(Z,P,k); G1,H1=c2d(A,B,T); G2,H2,C2,D2=c2dm(A,B,C,D,T,zoh); G3,H3,C3,D3=c2dm(A,B,C,D,T,foh);,3、脈沖傳遞函數(shù)