數(shù)值積分法仿真

1、 本章主要教學(xué)內(nèi)容 本章主要介紹控制系統(tǒng)數(shù)學(xué)模型的相關(guān)知識(shí)，通過本章的學(xué)習(xí)，讀者應(yīng)掌握以下內(nèi)容：求解常微分方程數(shù)值解的一般方法數(shù)值積分法的基本概念及其常用方法以系統(tǒng)微分方程或傳遞函數(shù)作為數(shù)學(xué)模型的仿真過程及程序設(shè)計(jì)方法以系統(tǒng)動(dòng)態(tài)結(jié)構(gòu)圖作為數(shù)學(xué)模型的仿真過程及程序設(shè)計(jì)方法仿真步長(zhǎng)的選擇與系統(tǒng)仿真精度和穩(wěn)定性的對(duì)應(yīng)關(guān)系快速仿真算法的概念、特點(diǎn)及其應(yīng)用 第6章 數(shù)值積分法仿真 16.1 數(shù)值積分法 6.1.1 概述 數(shù)字仿真模型、算法及仿真工具 控制系統(tǒng)的數(shù)字仿真是利用數(shù)字計(jì)算機(jī)作為仿真工具，采用數(shù)學(xué)上的各種數(shù)值算法求解控制系統(tǒng)運(yùn)動(dòng)的微分方程，得到被控物理量的運(yùn)動(dòng)規(guī)律。通常，計(jì)算機(jī)模擬被控對(duì)象是用一

2、定的仿真算法來實(shí)現(xiàn)被控對(duì)象的運(yùn)動(dòng)規(guī)律，這是基于被控對(duì)象的數(shù)學(xué)模型來完成的。 控制系統(tǒng)的數(shù)學(xué)模型經(jīng)過合理的近似及簡(jiǎn)化，大多數(shù)建立為常微分方程的表達(dá)形式。由于數(shù)學(xué)計(jì)算的難度和實(shí)際系統(tǒng)的復(fù)雜程度，在實(shí)際中遇到的大部分微分方程難以得到其解析解，通常都是通過數(shù)字計(jì)算機(jī)采用數(shù)值計(jì)算的方法來求取其數(shù)值解。在高級(jí)仿真軟件（例如MATLAB）環(huán)境下，已提供了功能十分強(qiáng)大、且能保證相應(yīng)精度的數(shù)值求解的功能函數(shù)或程序段，使用者僅需要按規(guī)定的語言規(guī)格調(diào)用即可，而無需從數(shù)值算法的底層考慮其編程實(shí)現(xiàn)過程。26.1.2 離散化原理 在數(shù)字計(jì)算機(jī)上對(duì)連續(xù)系統(tǒng)進(jìn)行仿真時(shí)，首先遇到的問題是，數(shù)字計(jì)算機(jī)的數(shù)值及時(shí)間都是離散的（計(jì)算

3、精度，指令執(zhí)行時(shí)間），而被仿真系統(tǒng)的數(shù)值和時(shí)間是連續(xù)的，后者如何用前者來實(shí)現(xiàn)？ 設(shè)系統(tǒng)模型為： ,其中u(t)為輸入變量，y(t)為系統(tǒng)狀態(tài)變量。令仿真時(shí)間間隔為h,離散化后的輸入變量為u(tk), 其中tk表示t=kh。如果 u(tk) u(tk), y(tk) y(tk)，則認(rèn)為兩模型等價(jià)，稱為相似原理。 對(duì)仿真建模方法有三個(gè)基本要求： 1、穩(wěn)定性，若原系統(tǒng)是穩(wěn)定的，則離散化后的仿真模型也得是穩(wěn)定的2、準(zhǔn)確性，絕對(duì)或相對(duì)誤差小于規(guī)定誤差3、快速性，數(shù)字仿真是一步步推進(jìn)的，由某個(gè)初始值y0出發(fā)，依次計(jì)算出y1、y2yk，每一步計(jì)算所需時(shí)間決定了仿真速度。3 6.1.3 數(shù)值積分法 一般情況下

4、，在控制系統(tǒng)仿真中最常用、最基本的求解常微分方程數(shù)值解的方法主要是數(shù)值積分法。 對(duì)于形如 的系統(tǒng)，已知系統(tǒng)狀態(tài)變量y的初值y0,現(xiàn)要計(jì)算y隨時(shí)間變化的過程y(t)，對(duì)微分方程的積分可以寫作：右圖所示曲線下的面積就是y(t),由于難以得到積分的數(shù)值表達(dá)式，所以采用近似的方法，常用有三種形式：歐拉法梯形法龍格一庫塔法 4歐拉公式,采用矩形面積近似積分結(jié)果，即當(dāng)t= tk+1時(shí)hk=tk+1-tk,若步距不變，則hk=h.為了提高精度，人們提出了“梯形法”，其中最簡(jiǎn)單的是亞當(dāng)姆斯預(yù)報(bào)校正公式 ，先用歐拉法估計(jì)近似值，然后用梯形法進(jìn)行校正：5龍格-庫塔法的基本原理在連續(xù)系統(tǒng)仿真中，主要的數(shù)值計(jì)算工作是

5、求解一階微分方程：若已知y的初值y0,再按離散化原理，對(duì)上式我們可以寫成：再對(duì)上式的右端函數(shù)f(t,y)(為任意非線性函數(shù))在tk附近展開成泰勒級(jí)數(shù)，依照展開的階次不同我們就構(gòu)成了不同的龍格-庫塔公式。6二階龍格庫塔公式,記在tk時(shí)刻的狀態(tài)變量為yk,并定義兩個(gè)附加向量型變量 :7四階龍格庫塔公式 ：不論幾階RK法，它們的計(jì)算公式都是由兩部分組成，即上一步的結(jié)果yk和步長(zhǎng)h乖以tk至tk+1時(shí)間間隔間各外推點(diǎn)的導(dǎo)數(shù)ki的加權(quán)平均和8有了上面的數(shù)學(xué)算法，就可以用MATLA編寫出該算法的函數(shù)： functiontout,yout = rk4(ode) t0= tspan(1);th= tspan(

6、2); if length(tspan) 0,本步失敗，但仍采用上式，重新進(jìn)行積分。196.3 面向微分方程的仿真程序設(shè)計(jì)6.3.1 通用仿真程序的一般結(jié)構(gòu)及工作原理1. 通用仿真程序的基本結(jié)構(gòu)以數(shù)字計(jì)算機(jī)作為仿真工具，使用適當(dāng)?shù)乃惴ㄕZ言來編制通用的仿真程序，可以針對(duì)不同的系統(tǒng)進(jìn)行相應(yīng)的仿真處理。按常規(guī)組成結(jié)構(gòu)，通用仿真程序可分為3個(gè)層次，即主程序塊、功能程序塊、基本子程序塊。各模塊功分析能如下：（1）主程序：完成仿真邏輯控制，實(shí)現(xiàn)各功能模塊的調(diào)用、模型的轉(zhuǎn)換、系統(tǒng)運(yùn)行、輸入輸出的控制等。20（2）初始化程序：完成各類初始數(shù)據(jù)的準(zhǔn)備工作，如設(shè)置工作單元、給定變量初值和系統(tǒng)仿真參數(shù)等。（3）運(yùn)行

7、程序：實(shí)現(xiàn)系統(tǒng)運(yùn)行的控制，調(diào)用數(shù)值積分法完成仿真算法處理，得出系統(tǒng)的響應(yīng)結(jié)果。 （4）輸出程序：按用戶指定的輸出形式，可以在顯示器、打印機(jī)、繪圖儀等設(shè)備上將仿真的結(jié)果以數(shù)據(jù)、動(dòng)態(tài)曲線、圖形等方式輸出。 212. 仿真程序流程框圖 22MATLAB提供了兩個(gè)常微分方程求解的函數(shù)ode23()、ode45(),分別采用了二階三級(jí)的RKF方法和四階五級(jí)的RKF法，并采用自適應(yīng)變步長(zhǎng)的求解方法，即當(dāng)解的變化較慢時(shí)采用較大的計(jì)算步長(zhǎng)，從而使得計(jì)算速度很快；當(dāng)解的變化較快時(shí)，步長(zhǎng)會(huì)自動(dòng)變小，從而使計(jì)算精度提高。SyntaxT,Y = solver(odefun,tspan,y0)T,Y = solver(

8、odefun,tspan,y0,options)where solver is one of ode45, ode23, ode113, ode15s, ode23s, ode23t, or ode23tb. ArgumentsodefunA function that evaluates the right-hand side of the differential equations. All solvers solve systems of equations in the form or problems that involve a mass matrix, . The ode23s

9、 solver can solve only equations with constant mass matrices. ode15s and ode23t can solve problems with a mass matrix that is singular, i.e., differential-algebraic equations (DAEs). tspanA vector specifying the interval of integration, t0,tf. To obtain solutions at specific times (all increasing or

10、 all decreasing), use tspan = t0,t1,.,tf.y0A vector of initial conditions.optionsOptional integration argument created using the odeset function. See odeset for details.p1,p2.Optional parameters that the solver passes to odefun and all the functions specified in options23Description T,Y = solver(ode

11、fun,tspan,y0) with tspan = t0 tf integrates the system of differential equations from time t0 to tf with initial conditions y0. Function f = odefun(t,y), for a scalar t and a column vector y, must return a column vector f corresponding to . Each row in the solution array Y corresponds to a time retu

12、rned in column vector T. To obtain solutions at the specific times t0, t1,.,tf (all increasing or all decreasing), use tspan = t0,t1,.,tf. T,Y = solver(odefun,tspan,y0,options) solves as above with default integration parameters replaced by property values specified in options, an argument created w

13、ith the odeset function. Commonly used properties include a scalar relative error tolerance RelTol (1e-3 by default) and a vector of absolute error tolerances AbsTol (all components are 1e-6 by default). See odeset for details.Example 1. An example of a nonstiff system is the system of equations des

14、cribing the motion of a rigid body without external forces. 24To simulate this system, create a function rigid containing the equations function dy = rigid(t,y)dy = zeros(3,1); % a column vectordy(1) = y(2) * y(3);dy(2) = -y(1) * y(3);dy(3) = -0.51 * y(1) * y(2);In this example we change the error t

15、olerances using the odeset command and solve on a time interval 0 12 with an initial condition vector 0 1 1 at time 0. options = odeset(RelTol,1e-6,AbsTol,1e-6 1e-6 1e-5);T,Y = ode45(rigid,0 12,0 1 1,options);Plotting the columns of the returned array Y versus T shows the solution plot(T,Y(:,1),-,T,

16、Y(:,2),-.,T,Y(:,3),.)2526Example 2. An example of a stiff system is provided by the van der Pol equations in relaxation oscillation. The limit cycle has portions where the solution components change slowly and the problem is quite stiff, alternating with regions of very sharp change where it is not

17、stiff. To simulate this system, create a function vdp1000 containing the equations function dy = vdp1000(t,y)dy = zeros(2,1); % a column vectordy(1) = y(2);dy(2) = 1000*(1 - y(1)2)*y(2) - y(1);For this problem, we will use the default relative and absolute tolerances (1e-3 and 1e-6, respectively) an

18、d solve on a time interval of 0 3000 with initial condition vector 2 0 at time 0. T,Y = ode15s(vdp1000,0 3000,2 0);plot(T,Y(:,1),-o)2728Plot(t,y(:,2)296.6 面向結(jié)構(gòu)圖的仿真程序設(shè)計(jì)面向結(jié)構(gòu)圖的線性系統(tǒng)仿真基本思想為：（1） 把一個(gè)復(fù)雜的高階線性系統(tǒng)化成由若干典型環(huán)節(jié)組成的模擬結(jié)構(gòu)圖表示。（2） 將各典型環(huán)節(jié)參數(shù)以及系統(tǒng)各環(huán)節(jié)的連接關(guān)系輸入計(jì)算機(jī)。（3） 仿真程序?qū)⑤斎氲南到y(tǒng)模型自動(dòng)轉(zhuǎn)化為狀態(tài)空間描述，即狀態(tài)方程形式。（4） 調(diào)用數(shù)值積分法求解

19、，并輸出仿真結(jié)果。 30典型環(huán)節(jié)的確定及算法描述典型環(huán)節(jié)的選擇是重要的一個(gè)步驟，它應(yīng)具備下述兩個(gè)原則：（1） 典型性由它可方便地組成其它任何形式的動(dòng)態(tài)環(huán)節(jié)。（2） 簡(jiǎn)易性由它組成的系統(tǒng)簡(jiǎn)便，計(jì)算機(jī)編程容易實(shí)現(xiàn)。常見的動(dòng)態(tài)環(huán)節(jié)根據(jù)控制理論可知，在實(shí)際控制系統(tǒng)中比較常見的動(dòng)態(tài)環(huán)節(jié)主要有以下五種：（1） 積分環(huán)節(jié) （2） 比例積分環(huán)節(jié) （3） 慣性環(huán)節(jié) （4） 一階超前（或滯后）環(huán)節(jié) （5） 二階振蕩環(huán)節(jié) 316.5 快速仿真算法 在對(duì)系統(tǒng)進(jìn)行仿真時(shí)，會(huì)碰到較高階次的控制系統(tǒng)，由于采用的計(jì)算機(jī)檔次不高會(huì)影響到仿真計(jì)算速度，占用較長(zhǎng)的機(jī)時(shí)；在參數(shù)尋優(yōu)時(shí)往往需要對(duì)控制系統(tǒng)進(jìn)行反復(fù)的仿真計(jì)算，也將使計(jì)算過

20、程加長(zhǎng)；此外，系統(tǒng)的實(shí)時(shí)仿真也會(huì)對(duì)仿真的快速性提出較高的要求。對(duì)于前面所討論的數(shù)值積分法由于有相應(yīng)的計(jì)算工作量，單純加大仿真步長(zhǎng)會(huì)影響到系統(tǒng)仿真的精度和穩(wěn)定性問題。 本節(jié)介紹幾種常用的快速仿真方法，采用這些方法來編制仿真計(jì)算子程序，可以彌補(bǔ)數(shù)值積分法仿真在速度上的缺陷，便于在實(shí)際工程中系統(tǒng)仿真時(shí)合理地加以選擇，達(dá)到提高系統(tǒng)仿真速度的最終目的。326.5.1 時(shí)域矩陣法 時(shí)域矩陣法是一種在時(shí)域內(nèi)采用無窮矩陣進(jìn)行系統(tǒng)仿真的算法，它每一步的計(jì)算量較小，而且與系統(tǒng)階次無關(guān)，適合于系統(tǒng)的快速仿真。 采用時(shí)域矩陣法來分析和討論系統(tǒng)的動(dòng)態(tài)性能具備下述特點(diǎn)：（1）時(shí)域矩陣法多用于采樣控制系統(tǒng)，由于采用脈沖過程

21、函數(shù)g（t）來計(jì)算系統(tǒng)的閉環(huán)響應(yīng)，不會(huì)因系統(tǒng)階次的增加而加大計(jì)算工作量，從而提高了仿真速度；但有時(shí)求解高階系統(tǒng)的脈沖過程函數(shù)g（t）會(huì)有一定的難度。（2）由于每個(gè)采樣時(shí)刻的g（k）是準(zhǔn)確計(jì)算出來的，所以采用時(shí)域矩陣法仿真時(shí)系統(tǒng)的采樣周期（或仿真步距）可以選得大些。（3）時(shí)域矩陣法可推廣到非線性系統(tǒng)的快速仿真。 336.5.2 增廣矩陣法 增廣矩陣法是將系統(tǒng)的控制量增廣到狀態(tài)變量中，使原來的非奇次常微分方程變?yōu)橐粋€(gè)齊次方程。6.5.3 替換法 快速仿真的系統(tǒng)通常比較關(guān)注系統(tǒng)仿真的速度應(yīng)該達(dá)到規(guī)定的要求，而對(duì)精度一般不做太高的要求。對(duì)于一個(gè)高階系統(tǒng)，如果能從它的傳遞函數(shù)G（s）直接推導(dǎo)出與之相匹配的并且允許較大采樣周期T的脈沖傳遞函數(shù)G（z），然后由G（z）獲得仿真