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

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