《建模與仿真及其醫(yī)學(xué)應(yīng)用》(精)

1、建模與仿真及其醫(yī)學(xué)應(yīng)用建模與仿真及其醫(yī)學(xué)應(yīng)用實(shí)驗(yàn)講義實(shí)驗(yàn)講義天津醫(yī)科大學(xué)生物醫(yī)學(xué)工程系天津醫(yī)科大學(xué)生物醫(yī)學(xué)工程系20042004 年年1實(shí)驗(yàn)一系統(tǒng)建模的實(shí)驗(yàn)一系統(tǒng)建模的 MATLABMATLAB 實(shí)現(xiàn)實(shí)現(xiàn)一、實(shí)驗(yàn)?zāi)康模阂?、?shí)驗(yàn)?zāi)康模?學(xué)習(xí) MATLAB 基本知識。2. 掌握數(shù)學(xué)模型的 MATLAB 實(shí)現(xiàn)：時域模型、狀態(tài)空間模型和零極點(diǎn)模型。3. 學(xué)習(xí)用 MATLAB 實(shí)現(xiàn)系統(tǒng)外部模型到內(nèi)部模型的轉(zhuǎn)換。4. 學(xué)習(xí)用 MATLAB 實(shí)現(xiàn)系統(tǒng)模型的連接：串聯(lián)、并聯(lián)、反饋連接。5. 了解模型降階的 MATLAB 實(shí)現(xiàn)。二、實(shí)驗(yàn)內(nèi)容二、實(shí)驗(yàn)內(nèi)容1.系統(tǒng)的實(shí)現(xiàn)、外部模型到內(nèi)部模型的轉(zhuǎn)換給定連續(xù)系統(tǒng)的傳遞函

2、數(shù)&=(2：器罷，利用MATLAB 建立傳遞函數(shù)模型，微分方程，并轉(zhuǎn)換為狀態(tài)空間模型。(2)已知某系統(tǒng)的狀態(tài)方程的系數(shù)矩陣為：_00_d=1001利用 MATLAB 建立狀態(tài)空間模型，并將其轉(zhuǎn)換為傳遞函數(shù)模型和零極點(diǎn)模型。MATLAB 轉(zhuǎn)換為傳遞函數(shù)模型和狀態(tài)空間模型。2.系統(tǒng)的離散、連接、降階給定連續(xù)系統(tǒng)的傳遞函數(shù) a=(2：；8(3：罷，將該連續(xù)系統(tǒng)的傳遞函數(shù)用零階重構(gòu)器和一階重構(gòu)器轉(zhuǎn)換為離散型傳遞函數(shù)，抽樣時 001 0a=b=c= 11-2-31 102(3)已知系統(tǒng)的零極點(diǎn)傳遞函數(shù)為G(s)=2(s+1)(s+2)(s+3)(s+4)，利用2間 T=1 秒。(2)該系統(tǒng)與系統(tǒng)

3、 H(s)=1分別串聯(lián)并聯(lián)負(fù)反饋連s2+6s+5接，求出組成的新系統(tǒng)的傳遞函數(shù)模型。(3) 將串聯(lián)組成的新系統(tǒng)進(jìn)行降階處理，求出降階后系統(tǒng)的模型，并用 plot 圖形比較降階前后系統(tǒng)的階躍響應(yīng)。要求：將以上過程用 MATLAB 編程(M 文件)實(shí)現(xiàn)，運(yùn)行輸出結(jié)果。三、實(shí)驗(yàn)說明關(guān)于系統(tǒng)建模的主要 MATLAB 函數(shù) 1建立傳遞函數(shù)模型：tf 函數(shù)：格式：sys=tf(num,den)num=bm,bm-1,b0分子多項(xiàng)式系數(shù)den=an,an-1,a0分母多項(xiàng)式系數(shù)2建立狀態(tài)空間模型：ss 函數(shù)：格式：sys=ss(a,b,c,d)%a,b,c,d為狀態(tài)方程系數(shù)矩陣 sys=ss(a,b,c,d

4、,T)%產(chǎn)生離散時間狀態(tài)空間模型 3建立零極點(diǎn)模型的函數(shù)：zpk格式：sys=zpk(z,p,k)4模型轉(zhuǎn)換函數(shù)：tf2sstf2zpss2tfss2zpzp2tfzp2ss%2 為 to 的意思格式：a,b,c,d=tf2ss(num,den)z,p,k=tf2zp(num,den)num,den=ss2tf(a,b,c,d,iu)%iu 指定是哪個輸入z,p,k=ss2zp(a,b,c,d,iu)3num,den=zp2tf(z,p,k)a,b,c,d=zp2ss(z,p,k)5模型的連接串聯(lián)：sys=series(sysl,sys2)并聯(lián)：sys=parallel(sys1,sys2)反

5、饋連接： sys=feedback(sys1,sys2,sign)%負(fù)反饋時 sign 可忽略； 正反饋時為 1。6系統(tǒng)擴(kuò)展：把若干個子系統(tǒng)組成系統(tǒng)組。格式：sys=append(sysl,sys2,.)7模型降階(1)基于平衡的狀態(tài)空間實(shí)現(xiàn)-balreal格式：sysb=balreal(sys)sysh,g,T,Ti=balreal(sys)sys為原系統(tǒng)，sysb(sysh)為平衡實(shí)現(xiàn)系統(tǒng)，g 為平衡對角線矩陣，T 為狀態(tài)變換矩陣，Ti是前者的逆矩陣。兩種格式的區(qū)別：前者只給出原系統(tǒng)的一個平衡的狀態(tài)空間實(shí)現(xiàn)，而后者還給出平衡實(shí)現(xiàn)的對角線矩陣 g， 從中可以看出哪個狀態(tài)變量該保留，哪個狀態(tài)變

6、量該刪去，從而實(shí)現(xiàn)降階。(2)降階的實(shí)現(xiàn)一 modred格式：rsys=modred(sys,elim)rsys=modred(sys,elim,mde)rsys=modred(sys,elim,del)強(qiáng)調(diào)：這里的 sys應(yīng)是函數(shù) balreal()變換的模型，elim為待消去的狀態(tài)，mde指降階中保持增益匹配，del指降階中不保持增益匹配。48連續(xù)系統(tǒng)模型離散化函數(shù)：C2DMConversionofcontinuousLTIsystemstoiscrete-time.格式：Ad,Bd,Cd,Dd二 C2DM(A,B,C,D,Ts,method)將連續(xù)系統(tǒng)狀態(tài)空間離散系統(tǒng)狀態(tài)空間method

7、:zoh零階重構(gòu)器 zeroorderholdfoh一階重構(gòu)器 firstorderholdNUMd,DENd=C2DM(NUM,DEN,Ts,method)將連續(xù)系統(tǒng)傳遞函數(shù)離散系統(tǒng)傳遞函數(shù)G(s)=NUM(s)/DEN(s)toG(z)=NUMd(z)/DENd(z).四、實(shí)驗(yàn)報告要求四、實(shí)驗(yàn)報告要求1整理好經(jīng)過運(yùn)行并證明是正確的程序，必要的地方加上注釋2給出實(shí)驗(yàn)的結(jié)果。5實(shí)驗(yàn)二連續(xù)系統(tǒng)的數(shù)字仿真實(shí)驗(yàn)二連續(xù)系統(tǒng)的數(shù)字仿真一、計(jì)算機(jī)仿真一、計(jì)算機(jī)仿真在計(jì)算機(jī)支持下進(jìn)行的現(xiàn)代仿真技術(shù)稱為計(jì)算機(jī)仿真。仿真不單純是對模型的實(shí)驗(yàn)，它包括建立模型、仿真運(yùn)行和分析研究仿真結(jié)果，即建模實(shí)驗(yàn)分析的全過程。M

8、ATLAB 提供各種用于系統(tǒng)仿真的函數(shù)，用戶可以通過 m 文件調(diào)用指令和函數(shù)進(jìn)行系統(tǒng)仿真，也可以通過 Simulink 工具箱，進(jìn)行面向系統(tǒng)結(jié)構(gòu)方框圖的系統(tǒng)仿真。這兩種方式可解決任意復(fù)雜系統(tǒng)的動態(tài)仿真問題，前者編輯靈活，而后者直觀性強(qiáng)，實(shí)現(xiàn)可視化編輯。內(nèi)容：連續(xù)系統(tǒng)仿真：數(shù)值積分法、離散相似法離散事件系統(tǒng)仿真SIMULINK 動態(tài)仿真二、基于數(shù)值積分法的連續(xù)系統(tǒng)仿真二、基于數(shù)值積分法的連續(xù)系統(tǒng)仿真1.數(shù)值積分法的 MATLAB 函數(shù)MATLAB 的工具箱提供了各種數(shù)值積分方法函數(shù)：格式：T,Y=solver(F,TSPAN,Yo,OPTIONS)solver為微分方程的求解函數(shù)名。F為系統(tǒng)模型

9、文件名，模型為yQ 二f(t,y)TSPAN=ToTfinal為積分區(qū)間，初值一終值，Yo 為系統(tǒng)輸出初始值，即 To 時刻的初值列向量；OPTIONS 設(shè)置積分相對允誤RelTol和絕對允誤AbsTol， 缺省時，RelTol=1e-3,AbsTol=1e-6.輸出參數(shù) T 和 Y 為列向量，T 為時刻向量，Y 表示不同時刻的函6數(shù)值。系統(tǒng)模型函數(shù)的編寫格式是固定的，如果其格式?jīng)]有按照要求去編寫則將得出錯誤的求解結(jié)果，系統(tǒng)模型函數(shù)的引導(dǎo)語句為：functionxdot=模型函數(shù)名（t,x，附加參數(shù)）其中 t為時間變量，x 為狀態(tài)變量，xdot 為狀態(tài)變量的導(dǎo)數(shù)。如果有附加參數(shù)需要傳遞，則可以

10、列出，中間用逗號分開。solver:ode23Runge-Kutta 法三階積分算法、二階誤差估計(jì)、變積分步長的低階算法ode45Runge-Kutta 法，變步長的中等階次積分算法odell3 變階的 Adams-Bashforth-Moulton，多步長ode15s改進(jìn)的 Gear法，用于剛性方程的求解。例：求微分方程 x=px+5，0t10，x0=1先建立一個系統(tǒng)模型文件(m 文件函數(shù))dfun.mfunctiony=dfun(t,x)y=sqrt(x)+5;然后建立 m文件 mp2-1%mp2-1t,x=ode23(dfun,010,1)plot(t,x)結(jié)果：tx01.00000.0

11、1331.08030.08001.48900.27202.72630.56854.78001.03568.30351.758914.34052.758923.67783.758934.03414.758945.32145.758957.48156.758970.472177.758984.26128.758998.82309.7589114.136510.0000117.93842.對于高階常微分方程，y(n)=f(t,y,y，y(”一1),則可以選擇一組狀態(tài)變量x=y,x=y,.,x=y(n-1)，將原高階微分方程模型變換成以12x下的一階微分方程組形式：x x=x12x x=x x2.3x

12、 x=f f(/,x,x,.x)n12n例：y+嘰y2-1)yr+y=0可變換成 x=x,x=一卩(x2-1)x-x122121functiomy=vdp_eq(t,x,mu)y=x(2);-mu*(x(1).人2-1).*x(2)-x(1)三、基于離散相似法的連續(xù)系統(tǒng)仿真三、基于離散相似法的連續(xù)系統(tǒng)仿真所謂離散相似法是首先將連續(xù)系統(tǒng)模型離散化，得到等價的或相似的離散化的模型，然后對相似的離散模型進(jìn)行仿真計(jì)算。根據(jù)這一原理，首先應(yīng)將連續(xù)時間系統(tǒng)模型轉(zhuǎn)換為等價的離散時間系統(tǒng)模型。連續(xù)系統(tǒng)離散化處理是通過轉(zhuǎn)移矩陣法；采樣和信號保持器；變換法(如雙線性變換)來實(shí)現(xiàn)的。1.轉(zhuǎn)移矩陣法的實(shí)現(xiàn)：如果連續(xù)系

13、統(tǒng)的狀態(tài)空間模型為：8|x=Ax+buy=Cx+Du則其離散狀態(tài)空間模型為：x(k+1)二(T)x(k)+O(T T)u(k)my(k)二 Cx(k)+Du(k)其中(T)二eAT狀態(tài)轉(zhuǎn)移矩陣(矩陣指數(shù))(T)=fe eA(T-t)Bdtm0由此可知，利用狀態(tài)方程離散化時的主要問題是如何計(jì)算(T)、(T)。對于一階、二階環(huán)節(jié)，(T)、(T)可以用解析方法求出mm來，而對于高階及多輸入多輸出系統(tǒng)，就要采用數(shù)值解法。MATLAB提供了計(jì)算矩陣指數(shù)的函數(shù) expm，EXPMMatrixexponential.EXPM(X)isthematrixexponentialofX.EXPMiscompute

14、dusingascalingandsquaringalgorithmwithaPadeapproximation.EXPM1,EXPM2andEXPM3arealternativemethods.例：x 二AxAx+BuBu1-A， b 二0-11求(T)、(T)。m%mp2-2A=01;0-1;%DefinesystemmatricesB=0;1;t=0.1symstau%Definetautobesymbolicphi=expm(A*t)%SymbolicallycalculateeA(A*t)phim1=int(expm(A*(t-tau),tau,0,t)*Bphim=sym2poly

15、(phim1)%將符號運(yùn)算轉(zhuǎn)換為數(shù)值9結(jié)果：phi=1.00000.095200.9048phim1=-9/10+exp(-1/10)1-exp(-1/10)phim=0.00480.09522采樣和信號保持器以及雙線性變換法的實(shí)現(xiàn)：采樣和信號保持器以及雙線性變換法的實(shí)現(xiàn)：MATLAB 還提供了通過采樣和信號保持器以及雙線性變化法將連續(xù)系統(tǒng)模型轉(zhuǎn)換為離散時間系統(tǒng)模型的函數(shù) C2D，調(diào)用格式為sysd=c2d(sys,Ts,method)其中，sys為線性連續(xù)時間系統(tǒng)；Ts為采樣時間；sysd 為等價的離散時間系統(tǒng)。method 為離散化方法，可以選用：1zoh為零階保持器2foh為一階保持器4

16、prewarp 為改進(jìn)的雙線性變換法5matched使連續(xù)和離散系統(tǒng)具有匹配的 DC 增益器，采樣周期為TS二0.1s，求其離散化系統(tǒng)模型，并比較離散前后系統(tǒng)階躍響應(yīng)。tustion為雙線性變換法,Z=eT1+Ts21-Ts2s例：連續(xù)系統(tǒng)傳遞函數(shù)H(s)二二1e-0.35ss2+4s+5，采用一階采樣保持10用 MATLAB 編寫程序：%mp2-3sysc=tf(l-1,145,td,0.35);%timedelaysysd=c2d(sysc,0.l,foh)step(sysc,sysd);運(yùn)行結(jié)果：Transferfunction:0.0115z人3+0.0456z人2-0.0562z-0

17、.009104zA6-1.629zA5+0.6703zA4Samplingtime:0.1離散前后系統(tǒng)階躍響應(yīng)比較11四、實(shí)驗(yàn)內(nèi)容四、實(shí)驗(yàn)內(nèi)容1求解方程 y+y(y2-l)y+y=0在不同卩值的解，1卩=1,0t20；初值 y(0)=0.2,yr(0)=0.7；2卩=2,0t20；初值 y(0)=0.2,y(0)=0.7；3卩=1000,0t0.8,F=89.8675 存在極顯著的直線回歸關(guān)系。pplot（0,1,0,m+（m0）上面的命令畫一條（0,0）到（0,m）的直線，當(dāng) m0 直線平移 1 以保證所畫直線在模塊框區(qū)域內(nèi)。此外，還要把 drawingcoordinate 設(shè)置為norma

18、lized。（5）按 OK 鍵封裝關(guān)閉封裝編輯器。五實(shí)驗(yàn)內(nèi)容：五實(shí)驗(yàn)內(nèi)容：1.將攝氏溫度轉(zhuǎn)換為華氏溫度的公式模型（學(xué)號mp1.mdl）9將攝氏溫度轉(zhuǎn)換為華氏溫度的公式：T=9T+32F5cT 為華氏溫度，T 為攝氏溫度。FC252. 創(chuàng)建簡單的連續(xù)系統(tǒng)模型(學(xué)號mp2.mdl)給定微分方程建模：x(t)=-2x(t)+u(t)u(t)是幅度為 1 頻率為 1 弧度/分的方波。3. 給定如下系統(tǒng)(2)建立子系統(tǒng)，創(chuàng)建后的系統(tǒng)模型如下：4對 VanderPol 方程x=x,X二一卩(x2-1)x-x1221211將給定的微分方程模型建立 Simulink 模型。 (學(xué)號 mp4.mdl)提示：第一個方程可以認(rèn)為是將x(t)信號作為一個積分器的輸入端，IniOut2+SineWaveMuxScopeTrjnsftrFenTriniftrFcn1Sum(1)創(chuàng)建此系統(tǒng)模型(學(xué)子系統(tǒng)創(chuàng)建后的系統(tǒng)模型Out1262這樣積分器的輸出端則將成為x(t)信號。類似的，x(t)信號本身也12可以認(rèn)為是一個積分器的輸出，在積分器的輸入端信號應(yīng)該為嘰Xj1)X2-X!，在構(gòu)造該信號時還應(yīng)該使用信號乘積的處理