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

第四章連續(xù)系統(tǒng)的離散化方法 4 1常微分方程的數(shù)值解法 一 數(shù)值求解的基本概念 設(shè)微分方程為 則求解方程中函數(shù)x t 問(wèn)題的常微分方程初值問(wèn)題 所謂數(shù)值求解就是要在時(shí)間區(qū)間 a b 中取若干離散點(diǎn) 求出微分方程在這些時(shí)刻的近似值 取前兩項(xiàng)近似 這種方法的幾何意義就是把f t x 在區(qū)間 tk tk 1 內(nèi)的曲邊面積用矩形面積近似代替 計(jì)算簡(jiǎn)單 計(jì)算量小 而且可以自啟動(dòng) 當(dāng)h很小時(shí) 造成的誤差是允許的 該算法具有一階精度 取k 0 1 2 N 從t0開(kāi)始 逐點(diǎn)遞推求解t1時(shí)的y1 t2時(shí)的y2 直至tn時(shí)的yn 稱(chēng)之為歐拉遞推公式 矩形面積 1 歐拉法 歐拉法的特點(diǎn) 導(dǎo)出簡(jiǎn)單 幾何意義明顯 便于理解 能說(shuō)明構(gòu)造數(shù)值解法一般計(jì)算公式的基本思想 通常用它來(lái)說(shuō)明有關(guān)的基本概念 例設(shè)系統(tǒng)方程為 用Euler法求其數(shù)值解 取步長(zhǎng) 遞推公式為 則 已知方程的解析解為精確解和解析解作比較 誤差在數(shù)量級(jí) 精度較差 2 龍格庫(kù)塔法 基本思想 取Taylor級(jí)數(shù)展開(kāi)式前三項(xiàng)近似求解 并利用線性組合代替導(dǎo)數(shù)的求解 既可避免計(jì)算高階導(dǎo)數(shù) 又可提高數(shù)值積分的精度 這就是Runge Kutta法的基本思想 2 龍格庫(kù)塔法 r為精度階次 ai為待定系數(shù) 由精度確定 ki用下式表示 線性組合 等各階導(dǎo)數(shù)不易計(jì)算 用下式中ki的線性組合代替 1 當(dāng)r 1時(shí) 與Taloy展開(kāi)式相比較 可得a1 1 則上式成為 歐拉遞推公式 2 當(dāng)r 2時(shí) 將在點(diǎn)展成Taylor級(jí)數(shù) 與臺(tái)勞公式的二階展開(kāi)近似公式相比 可得以下關(guān)系 三個(gè)方程 四個(gè)未知數(shù) 解不唯一 各個(gè)系數(shù)的幾種取法 見(jiàn)書(shū)上 3 r 4時(shí) 四階龍格庫(kù)塔公式 最常用 仿真中遇到的大多數(shù)工程實(shí)際問(wèn)題 四階龍格庫(kù)塔法以能滿足精度要求 其截?cái)嗾`差o h5 與h5同數(shù)量級(jí) 該法可以自啟動(dòng) 4 狀態(tài)空間四階龍格 庫(kù)塔遞推式若單輸入單輸出系統(tǒng)的狀態(tài)空間表達(dá)式為 在仿真中 對(duì)于n階系統(tǒng) 狀態(tài)方程可以寫(xiě)成一階微分方程 根據(jù)四階龍格 庫(kù)塔公式 有 T tk時(shí)刻的xi值 T tk h時(shí)刻的xi值 另 狀態(tài)方程的四階龍格 庫(kù)塔公式如下 RK法的特點(diǎn) 1需要存儲(chǔ)的數(shù)據(jù)少 占用的存儲(chǔ)空間少 2只需知道初值 即可啟動(dòng)遞推公式進(jìn)行計(jì)算 可自啟動(dòng) 3容易實(shí)現(xiàn)變步長(zhǎng)運(yùn)算 4每積分一步需要計(jì)算多次右函數(shù) 計(jì)算量大 基于龍格 庫(kù)塔法 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階龍格 庫(kù)塔算法 2 3階龍格 庫(kù)塔算法 3 常微分方程Matlab求解 解 令y1 x y2 x 1 建立M 文件vdp m如下 functiondy 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 010 1 1 plot T Y 1 T Y 2 3 結(jié)果 解1 建立m 文件rigid m如下 functiondy 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 012 011 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)的問(wèn)題 以及如何在仿真中恰當(dāng)?shù)倪x用這些方法 一般 選用數(shù)值算法從以下幾個(gè)方面考慮 1 精度受算法和h影響截?cái)嗾`差 舍入誤差 累計(jì)誤差 2 計(jì)算速度受算法和h影響算法簡(jiǎn)單 速度就快些 3 穩(wěn)定性受h影響 一般h 2 3 系統(tǒng)最小時(shí)間 4 3數(shù)值算法中的 病態(tài) 問(wèn)題 1 病態(tài) 常微分方程 例 其中 采用四階龍格庫(kù)塔法 h 0 01時(shí) 計(jì)算時(shí)間長(zhǎng) h 0 04時(shí) 誤差很大 當(dāng)h 0 05后 曲線發(fā)散振蕩 數(shù)值不穩(wěn)定 完全失去意義 系統(tǒng)矩陣的特征值差異較大 一般線性常微分方程組 的系數(shù)矩陣A的特征值具有如下特征 則稱(chēng)為 病態(tài) 方程 定義 2控制系統(tǒng)仿真中的 病態(tài) 問(wèn)題 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ì)算步長(zhǎng)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)過(guò)渡過(guò)程時(shí)間的長(zhǎng)短 一般與系統(tǒng)中具有最小時(shí)間常數(shù)Tmax的環(huán)節(jié)有關(guān) 要求計(jì)算步長(zhǎng)h取得很大 3對(duì)于病態(tài)問(wèn)題的仿真需要尋求更加合理的算法 以解決病態(tài)系統(tǒng)帶來(lái)的選取計(jì)算步長(zhǎng)與計(jì)算精度 計(jì)算時(shí)間之間的矛盾 3 病態(tài) 系統(tǒng)的仿真方法 采用穩(wěn)定性好 計(jì)算精度高的數(shù)值算法 并且允許計(jì)算步長(zhǎng)能根據(jù)系統(tǒng)性能動(dòng)態(tài)變化的情況在一定范圍內(nèi)作相應(yīng)的變化 采用隱式吉爾法 該法已經(jīng)證明對(duì)病態(tài)方程求解過(guò)程是數(shù)值穩(wěn)定的 隱式吉爾法從理論上十分適應(yīng)于病態(tài)系統(tǒng) 但需要解決好以下問(wèn)題 1 自啟動(dòng)r階多步算式無(wú)法自啟動(dòng) 需要用單步法求出前r步值 2 預(yù)估迭代迭代方法要求收斂性良好 否則在大步長(zhǎng)時(shí)會(huì)造成數(shù)值發(fā)散 3 變步長(zhǎng)初始階段采用小步長(zhǎng) 隨后可逐步放大步長(zhǎng) 對(duì)不同精度要求的系統(tǒng)仿真 要考慮變階次問(wèn)題 即為減小每一步計(jì)算的截?cái)嗾`差 以提高精度 應(yīng)選用較高的階次 而當(dāng)精度較低時(shí) 為減少工作量 則應(yīng)選取較低的階次 仿真時(shí)應(yīng)根據(jù)估計(jì)誤差與給定的誤差精度相比較改變步長(zhǎng)或階次來(lái)重新計(jì)算 4 4連續(xù)系統(tǒng)狀態(tài)方程的離散化 上章所述的連續(xù)系統(tǒng)數(shù)學(xué)模型的離散化 是通過(guò)數(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í) 疊代過(guò)程重復(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 0 t T 相當(dāng)在系統(tǒng)輸入端加一個(gè)采樣開(kāi)關(guān)和零階保持器 X k 1 T GX kT Hu kT G eAT 為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 則語(yǔ)句 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ù)求解 連續(xù)離散function nd dd tfc2d nc dc T A B C D tf2ss nc dc G H c2d A B T nd dd ss2tf G H C D 離散連續(xù)function nc dc tfd2c nd dd T A B C D tf2ss nd dd G H d2c A B T nc dc ss2tf G H C D