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

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

2、換為傳遞函數(shù)模型和零極點模型。(3)已知系統(tǒng)的零極點傳遞函數(shù)為，利用MATLAB轉(zhuǎn)換為傳遞函數(shù)模型和狀態(tài)空間模型。2系統(tǒng)的離散、連接、降階（1）給定連續(xù)系統(tǒng)的傳遞函數(shù)，將該連續(xù)系統(tǒng)的傳遞函數(shù)用零階重構(gòu)器和一階重構(gòu)器轉(zhuǎn)換為離散型傳遞函數(shù)，抽樣時間T=1秒。（2）該系統(tǒng)與系統(tǒng)分別串聯(lián)并聯(lián)負反饋連接，求出組成的新系統(tǒng)的傳遞函數(shù)模型。（3）將串聯(lián)組成的新系統(tǒng)進行降階處理，求出降階后系統(tǒng)的模型，并用plot圖形比較降階前后系統(tǒng)的階躍響應(yīng)。要求：將以上過程用MATLAB編程（M文件）實現(xiàn)，運行輸出結(jié)果。三、實驗說明關(guān)于系統(tǒng)建模的主要MATLAB函數(shù)1建立傳遞函數(shù)模型：tf函數(shù) ：格式：sys=tf(num

3、,den)num=bm,bm-1,b0 分子多項式系數(shù)den=an,an-1,a0 分母多項式系數(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,T) %產(chǎn)生離散時間狀態(tài)空間模型3建立零極點模型的函數(shù)：zpk格式：sys=zpk(z,p,k)4模型轉(zhuǎn)換函數(shù)：tf2ss tf2zp ss2tf ss2zp zp2tf zp2ss %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

4、,k=ss2zp(a,b,c,d,iu)num,den=zp2tf(z,p,k)a,b,c,d=zp2ss(z,p,k)5模型的連接串聯(lián)：sys=series(sys1,sys2)并聯(lián)：sys=parallel(sys1,sys2)反饋連接：sys=feedback(sys1,sys2,sign)%負反饋時sign可忽略；正反饋時為1。6系統(tǒng)擴展：把若干個子系統(tǒng)組成系統(tǒng)組。格式：sys=append(sys1,sys2,)7模型降階（1）基于平衡的狀態(tài)空間實現(xiàn)-balreal格式：sysb=balreal(sys)sysh,g,T,Ti=balreal(sys)sys為原系統(tǒng)，sysb(sys

5、h)為平衡實現(xiàn)系統(tǒng)，g為平衡對角線矩陣，T為狀態(tài)變換矩陣，Ti是前者的逆矩陣。兩種格式的區(qū)別：前者只給出原系統(tǒng)的一個平衡的狀態(tài)空間實現(xiàn)，而后者還給出平衡實現(xiàn)的對角線矩陣g，從中可以看出哪個狀態(tài)變量該保留，哪個狀態(tài)變量該刪去，從而實現(xiàn)降階。（2）降階的實現(xiàn)modred格式：rsys=modred(sys,elim)rsys=modred(sys,elim,mde)rsys=modred(sys,elim,del)強調(diào)：這里的sys應(yīng)是函數(shù)balreal()變換的模型，elim為待消去的狀態(tài)，mde指降階中保持增益匹配，del 指降階中不保持增益匹配。8連續(xù)系統(tǒng)模型離散化函數(shù)：C2DM Conve

6、rsion of continuous LTI systems to iscrete-time.格式：Ad,Bd,Cd,Dd=C2DM(A,B,C,D,Ts,method)將連續(xù)系統(tǒng)狀態(tài)空間離散系統(tǒng)狀態(tài)空間method: zoh 零階重構(gòu)器 zero order hold foh 一階重構(gòu)器 first order holdNUMd,DENd = C2DM(NUM,DEN,Ts,method) 將連續(xù)系統(tǒng)傳遞函數(shù)離散系統(tǒng)傳遞函數(shù)G(s) = NUM(s)/DEN(s) to G(z) = NUMd(z)/DENd(z).四、實驗報告要求1 整理好經(jīng)過運行并證明是正確的程序，必要的地方加上注釋。

7、2 給出實驗的結(jié)果。實驗二 連續(xù)系統(tǒng)的數(shù)字仿真一、計算機仿真在計算機支持下進行的現(xiàn)代仿真技術(shù)稱為計算機仿真。仿真不單純是對模型的實驗，它包括建立模型、仿真運行和分析研究仿真結(jié)果，即建模實驗分析的全過程。MATLAB提供各種用于系統(tǒng)仿真的函數(shù)，用戶可以通過m 文件調(diào)用指令和函數(shù)進行系統(tǒng)仿真，也可以通過Simulink工具箱，進行面向系統(tǒng)結(jié)構(gòu)方框圖的系統(tǒng)仿真。這兩種方式可解決任意復(fù)雜系統(tǒng)的動態(tài)仿真問題，前者編輯靈活，而后者直觀性強，實現(xiàn)可視化編輯。內(nèi)容：連續(xù)系統(tǒng)仿真：數(shù)值積分法、離散相似法離散事件系統(tǒng)仿真SIMULINK動態(tài)仿真二、基于數(shù)值積分法的連續(xù)系統(tǒng)仿真1數(shù)值積分法的MATLAB函數(shù)MATL

8、AB的工具箱提供了各種數(shù)值積分方法函數(shù)：格式：T,Y=solver(F,TSPAN,Yo,OPTIONS)solver為微分方程的求解函數(shù)名。F為系統(tǒng)模型文件名，模型為TSPAN=To Tfinal為積分區(qū)間，初值終值，Yo為系統(tǒng)輸出初始值，即To時刻的初值列向量；OPTIONS設(shè)置積分相對允誤RelTol和絕對允誤AbsTol，缺省時，RelTol=1e-3, AbsTol=1e-6.輸出參數(shù)T和Y為列向量，T為時刻向量，Y表示不同時刻的函數(shù)值。系統(tǒng)模型函數(shù)的編寫格式是固定的，如果其格式?jīng)]有按照要求去編寫則將得出錯誤的求解結(jié)果，系統(tǒng)模型函數(shù)的引導(dǎo)語句為：function xdot=模型函數(shù)名

9、(t,x,附加參數(shù))其中t為時間變量，x為狀態(tài)變量，xdot為狀態(tài)變量的導(dǎo)數(shù)。如果有附加參數(shù)需要傳遞，則可以列出，中間用逗號分開。solver:ode23 Runge-Kutta法 三階積分算法、二階誤差估計、變積分步長的低階算法ode45 Runge-Kutta法，變步長的中等階次積分算法ode113 變階的Adams-Bashforth-Moulton，多步長ode15s 改進的Gear法，用于剛性方程的求解。例：求微分方程，先建立一個系統(tǒng)模型文件（m文件函數(shù)）dfun.mfunction y=dfun(t,x)y=sqrt(x)+5;然后建立m文件mp2-1%mp2-1t,x=ode23

10、(dfun, 0 10 , 1)plot(t,x)結(jié)果： t x 0 1.0000 0.0133 1.0803 0.0800 1.4890 0.2720 2.7263 0.5685 4.7800 1.0356 8.3035 1.7589 14.3405 2.7589 23.6778 3.7589 34.0341 4.7589 45.3214 5.7589 57.4815 6.7589 70.4721 7.7589 84.2612 8.7589 98.8230 9.7589 114.1365 10.0000 117.93842對于高階常微分方程，,則可以選擇一組狀態(tài)變量，將原高階微分方程模型變換

11、成以下的一階微分方程組形式：例：可變換成functiom y=vdp_eq(t,x,mu)y=x(2);-mu*(x(1).2-1).*x(2)-x(1)三、基于離散相似法的連續(xù)系統(tǒng)仿真所謂離散相似法是首先將連續(xù)系統(tǒng)模型離散化，得到等價的或相似的離散化的模型，然后對相似的離散模型進行仿真計算。根據(jù)這一原理，首先應(yīng)將連續(xù)時間系統(tǒng)模型轉(zhuǎn)換為等價的離散時間系統(tǒng)模型。連續(xù)系統(tǒng)離散化處理是通過轉(zhuǎn)移矩陣法；采樣和信號保持器；變換法（如雙線性變換）來實現(xiàn)的。1轉(zhuǎn)移矩陣法的實現(xiàn)：如果連續(xù)系統(tǒng)的狀態(tài)空間模型為：則其離散狀態(tài)空間模型為：其中 狀態(tài)轉(zhuǎn)移矩陣（矩陣指數(shù)）由此可知，利用狀態(tài)方程離散化時的主要問題是如何計

12、算、。對于一階、二階環(huán)節(jié)，、可以用解析方法求出來，而對于高階及多輸入多輸出系統(tǒng)，就要采用數(shù)值解法。MATLAB提供了計算矩陣指數(shù)的函數(shù)expm，EXPM Matrix exponential.EXPM(X) is the matrix exponential of X. EXPM is computed using a scaling and squaring algorithm with a Pade approximation.EXPM1, EXPM2 and EXPM3 are alternative methods.例：， 求、。%mp2-2A = 0 1 ; 0 -1; % Defi

13、ne system matricesB = 0 ; 1;t=0.1syms tau % Define tau to be symbolicphi = expm(A*t) % Symbolically calculate e(A*t)phim1= int(expm(A*(t-tau),tau,0,t)*Bphim=sym2poly(phim1)%將符號運算轉(zhuǎn)換為數(shù)值結(jié)果：phi =1.0000 0.0952 0 0.9048phim1 = -9/10+exp(-1/10) 1-exp(-1/10) phim = 0.0048 0.09522采樣和信號保持器以及雙線性變換法的實現(xiàn)：MATLAB還提

14、供了通過采樣和信號保持器以及雙線性變化法將連續(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為離散化方法，可以選用： zoh 為零階保持器 foh為一階保持器 tustion為雙線性變換法， prewarp為改進的雙線性變換法 matched使連續(xù)和離散系統(tǒng)具有匹配的DC增益例：連續(xù)系統(tǒng)傳遞函數(shù)，采用一階采樣保持器，采樣周期為，求其離散化系統(tǒng)模型，并比較離散前后系統(tǒng)階躍響應(yīng)。用MATLAB編寫程序： %mp2-3sysc = tf ( l -

15、1 , 14 5 , td , 0.35 );%time delaysysd = c2d ( sysc, 0.l, foh )step ( sysc, sysd );運行結(jié)果：Transfer function:Sampling time: 0.1 離散前后系統(tǒng)階躍響應(yīng)比較四、實驗內(nèi)容1 求解方程在不同值的解，=1，；初值, ；=2，；初值, ；=1000，；初值, 。2. 給定系統(tǒng)1/(s+a)k/su x1 x2 y用轉(zhuǎn)移矩陣法仿真，其中k=2;a=1;T=0.01;x1(0)=0.1;x2(0)=0.7，仿真時間為0.2 秒，求系統(tǒng)的階躍響應(yīng)。采用零階和一階采樣保持器，求其離散化系統(tǒng)模型

16、，給出系統(tǒng)階躍響應(yīng)。采用雙線性變換法，求其離散化系統(tǒng)模型，給出系統(tǒng)階躍響應(yīng)。實驗三 最小二乘法及數(shù)據(jù)擬合建模的回歸分析一、實驗?zāi)康?1掌握用最小二乘建立回歸數(shù)學(xué)模型。2學(xué)習(xí)通過幾個數(shù)據(jù)擬合的回歸分析來判斷曲線（直線）擬合的精度，通過回歸分析來判斷模型建立是否正確。3應(yīng)用建立的模型進行預(yù)測。二、基本原理和方法1建立回歸數(shù)學(xué)模型在進行建模和仿真分析時，人們經(jīng)常面臨用已知系統(tǒng)實測數(shù)據(jù)應(yīng)用數(shù)學(xué)模型描述對應(yīng)系統(tǒng)，即對數(shù)據(jù)進行擬合。擬合的目的是尋找給定的曲線（直線），它在某種準則下最佳地擬合數(shù)據(jù)。最佳擬合要在什么準則下的最佳？以及用什么樣的曲線模型去擬合。常用的擬合方法之一是多項式的最小二乘擬合，其準則是

17、最小誤差平方和準則，所用的擬合曲線為多項式。本實驗在Matlab平臺上，以多項式最小二乘擬合為例，掌握回歸模型的建立（包括參數(shù)估計和模型建立）和用模型進行預(yù)測的方法，并學(xué)習(xí)回歸分析的基本方法。2在MATLAB里，用于求解最小二乘多項式擬合問題的函數(shù)如下： polyfit 最小二乘擬合 p=polyfit(x,y,n) 對輸入數(shù)據(jù)y的n階最小二乘擬合多項式p(x)的系數(shù) Y=polyval(p,x) 求多項式的函數(shù)值Y 以下是一個多項式擬合的例子。已知 x=0,1,2,3,4,5,6,7,8,9,10 共11個點（自變量），實測數(shù)據(jù)y=-0.447, 1.978, 3.28, 6.16, 7.0

18、8, 7.34, 7.66, 9.56, 9.48, 9.30, 11.2求：2階的預(yù)測方程，并用8階的預(yù)測方程與之比較。x=linspace(0,1,11); y=-.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.30 11.2; p=polyfit(x,y,2) %求2階的預(yù)測方程 的系數(shù) p= b2 b1 b0 z=polyval(p,x); %求預(yù)測的y值 （z表示） p2=polyfit(x,y,8) %求8階的預(yù)測方程 z1=polyval(p2,x); plot(x,y,om,x,z,:*rx,z1, :+b)圖中：”0” 代表散點

19、圖 “+”代表8階預(yù)測方程 “*”代表2階預(yù)測方程 圖1 散點圖與2階預(yù)測方程3回歸模型的檢驗回歸模型的檢驗是判斷數(shù)據(jù)擬合的好壞即模型建立的正確與否，為建立模型和應(yīng)用模型提供支持。在MATLAB平臺，用于回歸檢驗的語句如下： b,bint,r,rint,stats=regress(y，X，) 其中， 為回歸系數(shù) e 隨機誤差（均值為0，方差） b：回歸系數(shù)的估計值 bint：回歸系數(shù)的置信區(qū)間 r：殘差 rint：殘差的置信區(qū)間 stats：統(tǒng)計量 R2 F P 用此語句，可以得到回歸系數(shù)，復(fù)相關(guān)系數(shù)R，方差檢驗F值，P。rcoplot(r,rint) %打印殘差分布圖在圖中，若殘差的置信區(qū)間

20、不包含零點，則視為異常點，將其剔除后重新計算。b=leastsq(函數(shù)名,b0) % 非線性最小二乘法擬合， b0為初始值s=sqrt(sum(函數(shù)名.* 函數(shù)名)./（n-2） %計算剩余標準差s=sqrt(sum(y-Y).2)./(n-2) %計算剩余標準差4預(yù)測值： 用經(jīng)過檢驗的數(shù)學(xué)模型即可預(yù)測數(shù)據(jù)。即把x代入回歸方程對y進行估計，該估計值為。以下用一個例子說明回歸模型的檢驗與預(yù)測： 有人研究了黏蟲孵化歷期平均溫度（x,oC）與歷期天數(shù)（y,天）之間關(guān)系，試驗資料列入下表，求直線回歸方程，并進行檢驗。x,歷期平均溫度oC 11.8 14.7 15.6 16.8 17.1 18.8 19

21、.5 20.4Y，孵化歷期（天） 30.1 17.3 16.7 13.6 11.9 10.7 8.3 6.7先作出(xi,yi)的散點圖x=11.8 14.7 15.6 16.8 17.1 18.8 19.5 20.4;y=30.1 17.3 16.7 13.6 11.9 10.7 8.3 6.7plot(x,y,r+) 圖2 歷期平均溫度（x,oC）與歷期天數(shù)（y,天）的散點圖從圖中可見y與x基本，因此用一元線性回歸。X=ones(8,1),x;b,bint,r,rint,stats=regress(y,X)；b,bint,stats,rcoplot(r,rint)運行結(jié)果：b = 57.0393 -2.5317bint = 45.9035 68.1750 -3.1851 -1.8782stats = 0.9374 89.8675 0.0001得到以下結(jié)果：復(fù)相關(guān)系數(shù)的平方 R2=0.9374, 則R=0.96820.8,F=89.8675 存在極顯著的直線回歸關(guān)系。 pplot(0,1,0