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

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

2、換為傳遞函數(shù)模型和零極點模型。(3)已知系統(tǒng)的零極點傳遞函數(shù)為，利用MATLAB轉換為傳遞函數(shù)模型和狀態(tài)空間模型。2系統(tǒng)的離散、連接、降階（1）給定連續(xù)系統(tǒng)的傳遞函數(shù)，將該連續(xù)系統(tǒng)的傳遞函數(shù)用零階重構器和一階重構器轉換為離散型傳遞函數(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)的階躍響應。要求：將以上過程用MATLAB編程（M文件）實現(xià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模型轉換函數(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應是函數(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' 零階重構器 zero order hold 'foh' 一階重構器 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) = NUM

7、d(z)/DENd(z).四、實驗報告要求1 整理好經(jīng)過運行并證明是正確的程序，必要的地方加上注釋。2 給出實驗的結果。實驗二 連續(xù)系統(tǒng)的數(shù)字仿真一、計算機仿真在計算機支持下進行的現(xiàn)代仿真技術稱為計算機仿真。仿真不單純是對模型的實驗，它包括建立模型、仿真運行和分析研究仿真結果，即建模實驗分析的全過程。MATLAB提供各種用于系統(tǒng)仿真的函數(shù)，用戶可以通過m 文件調(diào)用指令和函數(shù)進行系統(tǒng)仿真，也可以通過Simulink工具箱，進行面向系統(tǒng)結構方框圖的系統(tǒng)仿真。這兩種方式可解決任意復雜系統(tǒng)的動態(tài)仿真問題，前者編輯靈活，而后者直觀性強，實現(xiàn)可視化編輯。內(nèi)容：連續(xù)系統(tǒng)仿真：數(shù)值積分法、離散相似法離散事件系

8、統(tǒng)仿真SIMULINK動態(tài)仿真二、基于數(shù)值積分法的連續(xù)系統(tǒng)仿真1數(shù)值積分法的MATLAB函數(shù)MATLAB的工具箱提供了各種數(shù)值積分方法函數(shù)：格式：T,Y=solver(F,TSPAN,Yo,OPTIONS)solver為微分方程的求解函數(shù)名。F為系統(tǒng)模型文件名，模型為TSPAN=To Tfinal為積分區(qū)間，初值終值，Yo為系統(tǒng)輸出初始值，即To時刻的初值列向量；OPTIONS設置積分相對允誤RelTol和絕對允誤AbsTol，缺省時，RelTol=1e-3, AbsTol=1e-6.輸出參數(shù)T和Y為列向量，T為時刻向量，Y表示不同時刻的函數(shù)值。系統(tǒng)模型函數(shù)的編寫格式是固定的，如果其格式?jīng)]有按

9、照要求去編寫則將得出錯誤的求解結果，系統(tǒng)模型函數(shù)的引導語句為：function xdot=模型函數(shù)名(t,x,附加參數(shù))其中t為時間變量，x為狀態(tài)變量，xdot為狀態(tài)變量的導數(shù)。如果有附加參數(shù)需要傳遞，則可以列出，中間用逗號分開。solver:ode23 Runge-Kutta法 三階積分算法、二階誤差估計、變積分步長的低階算法ode45 Runge-Kutta法，變步長的中等階次積分算法ode113 變階的Adams-Bashforth-Moulton，多步長ode15s 改進的Gear法，用于剛性方程的求解。例：求微分方程，先建立一個系統(tǒng)模型文件（m文件函數(shù)）function y=dfun

10、(t,x)y=sqrt(x)+5;然后建立m文件mp2-1%mp2-1t,x=ode23('dfun', 0 10 , 1)plot(t,x)結果： t x 2對于高階常微分方程，,則可以選擇一組狀態(tài)變量，將原高階微分方程模型變換成以下的一階微分方程組形式：例：可變換成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ù)這一原理，首先應將連續(xù)時間系統(tǒng)模型轉換為等價的離散時間系

11、統(tǒng)模型。連續(xù)系統(tǒng)離散化處理是通過轉移矩陣法；采樣和信號保持器；變換法（如雙線性變換）來實現(xiàn)的。1轉移矩陣法的實現(xiàn)：如果連續(xù)系統(tǒng)的狀態(tài)空間模型為：則其離散狀態(tài)空間模型為：其中 狀態(tài)轉移矩陣（矩陣指數(shù)）由此可知，利用狀態(tài)方程離散化時的主要問題是如何計算、。對于一階、二階環(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 squa

12、ring algorithm with a Pade approximation.EXPM1, EXPM2 and EXPM3 are alternative methods.例：， 求、。%mp2-2A = 0 1 ; 0 -1; % Define system matricesB = 0 ; 1;syms 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)%將符號運算轉換為數(shù)

13、值phim1 = -9/10+exp(-1/10) 1-exp(-1/10) 2采樣和信號保持器以及雙線性變換法的實現(xiàn)：MATLAB還提供了通過采樣和信號保持器以及雙線性變化法將連續(xù)系統(tǒng)模型轉換為離散時間系統(tǒng)模型的函數(shù)C2D，調(diào)用格式為sysd = c2d (sys, Ts, method)其中，sys為線性連續(xù)時間系統(tǒng)；Ts為采樣時間；sysd為等價的離散時間系統(tǒng)。method為離散化方法，可以選用： 'zoh '為零階保持器 'foh'為一階保持器 'tustion'為雙線性變換法， 'prewarp'為改進的雙線性變換法 &

14、#39;matched'使連續(xù)和離散系統(tǒng)具有匹配的DC增益例：連續(xù)系統(tǒng)傳遞函數(shù)，采用一階采樣保持器，采樣周期為，求其離散化系統(tǒng)模型，并比較離散前后系統(tǒng)階躍響應。用MATLAB編寫程序： %mp2-3sysc = tf ( l -1 , 14 5 , 'td' , 0.35 );%time delaysysd = c2d ( sysc, 0.l, 'foh' )step ( sysc, sysd );運行結果：Transfer function: 離散前后系統(tǒng)階躍響應比較四、實驗內(nèi)容1 求解方程在不同值的解，=1，；初值, ；=2，；初值, ；=1000，

15、；初值, 。2. 給定系統(tǒng)1/(s+a)k/su x1 x2 y用轉移矩陣法仿真，其中k=2;a=1;T=0.01;x1(0)=0.1;x2，仿真時間為0.2 秒，求系統(tǒng)的階躍響應。采用零階和一階采樣保持器，求其離散化系統(tǒng)模型，給出系統(tǒng)階躍響應。采用雙線性變換法，求其離散化系統(tǒng)模型，給出系統(tǒng)階躍響應。實驗三 最小二乘法及數(shù)據(jù)擬合建模的回歸分析一、實驗目的:1掌握用最小二乘建立回歸數(shù)學模型。2學習通過幾個數(shù)據(jù)擬合的回歸分析來判斷曲線（直線）擬合的精度，通過回歸分析來判斷模型建立是否正確。3應用建立的模型進行預測。二、基本原理和方法1建立回歸數(shù)學模型在進行建模和仿真分析時，人們經(jīng)常面臨用已知系統(tǒng)實

16、測數(shù)據(jù)應用數(shù)學模型描述對應系統(tǒng)，即對數(shù)據(jù)進行擬合。擬合的目的是尋找給定的曲線（直線），它在某種準則下最佳地擬合數(shù)據(jù)。最佳擬合要在什么準則下的最佳？以及用什么樣的曲線模型去擬合。常用的擬合方法之一是多項式的最小二乘擬合，其準則是最小誤差平方和準則，所用的擬合曲線為多項式。本實驗在Matlab平臺上，以多項式最小二乘擬合為例，掌握回歸模型的建立（包括參數(shù)估計和模型建立）和用模型進行預測的方法，并學習回歸分析的基本方法。2在MATLAB里，用于求解最小二乘多項式擬合問題的函數(shù)如下： polyfit 最小二乘擬合 p=polyfit(x,y,n) 對輸入數(shù)據(jù)y的n階最小二乘擬合多項式p(x)的系數(shù) Y

17、=polyval(p,x) 求多項式的函數(shù)值Y 以下是一個多項式擬合的例子。已知 x=0,1,2,3,4,5,6,7,8,9,10 共11個點（自變量），實測數(shù)據(jù)y=-0.447,求：2階的預測方程，并用8階的預測方程與之比較。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階的預測方程 的系數(shù) p= b2 b1 b0 z=polyval(p,x); %求預測的y值 （z表示） p2=polyfit(x,y,8) %求8階的預測方程 z1=po

18、lyval(p2,x); plot(x,y,'om',x,z,':*r'x,z1, ':+b')圖中：”0” 代表散點圖 “+”代表8階預測方程 “*”代表2階預測方程 圖1 散點圖與2階預測方程3回歸模型的檢驗回歸模型的檢驗是判斷數(shù)據(jù)擬合的好壞即模型建立的正確與否，為建立模型和應用模型提供支持。在MATLAB平臺，用于回歸檢驗的語句如下： b,bint,r,rint,stats=regress(y，X，) 其中， 為回歸系數(shù) e 隨機誤差（均值為0，方差） b：回歸系數(shù)的估計值 bint：回歸系數(shù)的置信區(qū)間 r：殘差 rint：殘差的置信區(qū)間

19、stats：統(tǒng)計量 R2 F P 用此語句，可以得到回歸系數(shù)，復相關系數(shù)R，方差檢驗F值，P。rcoplot(r,rint) %打印殘差分布圖在圖中，若殘差的置信區(qū)間不包含零點，則視為異常點，將其剔除后重新計算。b=leastsq(函數(shù)名,b0) % 非線性最小二乘法擬合， b0為初始值s=sqrt(sum(函數(shù)名.* 函數(shù)名)./（n-2） %計算剩余標準差s=sqrt(sum(y-Y).2)./(n-2) %計算剩余標準差4預測值： 用經(jīng)過檢驗的數(shù)學模型即可預測數(shù)據(jù)。即把x代入回歸方程對y進行估計，該估計值為。以下用一個例子說明回歸模型的檢驗與預測： 有人研究了黏蟲孵化歷期平均溫度（x,o

20、C）與歷期天數(shù)（y,天）之間關系，試驗資料列入下表，求直線回歸方程，并進行檢驗。x,歷期平均溫度oC 先作出(xi,yi)的散點圖 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.7'plot(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)運行結

21、果： 得到以下結果：復相關系數(shù)的平方 R2374, 則R=>0.8,F=89.8675 存在極顯著的直線回歸關系。 p<10-4,因此此回歸模型y接受。oC時，孵化歷期天數(shù)為多少天？（天） 圖三 殘差分布圖 （全部觀測點）觀察圖中的殘差分布，除第一點外，其余殘差的置信區(qū)間均包含零點，第一個點應視為異常點，將其剔除后重新計算，可得：圖三 殘差分布圖 （除掉第一個異常點觀測點） oC時，孵化歷期天數(shù) 顯然，此估計值應更加符合真實數(shù)據(jù)。三、實驗內(nèi)容：1 用給定的多項式 產(chǎn)生一組數(shù)據(jù)，再在yi上添加隨機干擾，然后用和添加了隨機干擾的yi作3次多項式擬合，與原系數(shù)比較。如果作2次或4次多項式

22、擬合，結果如何？2給定數(shù)據(jù)（xi,yi）見下表畫出散點圖觀察二者的關系，試建立合適的回歸模型：直線回歸方程；二次曲線；對數(shù)曲線；線性化的對數(shù)曲線等.并作回歸分析,得到最佳。得到的最佳回歸方程做出散點圖和曲線。xi23457810yixi111415151819yi實驗四 Simulink動態(tài)仿真SIMULINK是MATLAB重要軟件包，用于對動態(tài)系統(tǒng)建模和仿真，它適用于連續(xù)系統(tǒng)和離散系統(tǒng)，也適用線性系統(tǒng)和非線性系統(tǒng)。它采用系統(tǒng)模塊直觀地描述系統(tǒng)典型環(huán)節(jié)，因此可十分方便地建立系統(tǒng)模型而不需要花較多時間編程。正由于這些特點，SIMULINK廣泛流行，被認為是最受歡迎的仿真軟件。一、SIMULINK

23、的特點啟動：在命令窗口輸入Simulink或單擊“simulink library Browser”按鈕或使用菜單命令filenewmodel，均可打開“simulink library Browser”窗口，里面有許多Simulink基本模塊，用戶可以直接調(diào)用這些模塊。1支持圖形用戶界面。在命令窗口鍵入“thermo”(這是一個房間熱力學仿真演示程序)。里面包含許多模塊，模塊之間用直線相連，組成模型，模型造好后，可以仿真運行，等待結果。（單擊start按鈕或菜單命令simulinkstart,然后單擊thermo plots模塊，查看仿真結果室內(nèi)溫度、室外溫度和熱量曲線，這種和實際示波器輸出

24、相似的圖形化顯示，給用戶一個很友好的界面）。2層次性 雙擊“house”模塊，出現(xiàn)“thermo/house”窗口，表示house模塊是由窗口右邊所示的一些模塊連接而成。像house這樣由由幾個相互關聯(lián)的模塊組合而成的模塊在simulink中稱為子系統(tǒng)（后面講）。而這種一個模塊又由許多模塊組成的特性，這是simulink的層次性。為了和子系統(tǒng)相區(qū)別，這里把模型本身稱為模型的頂層系統(tǒng)。層次性的好處是建立的模型在結構上非常清晰整齊，一目了然，更重要的是，這一特性使得用戶可以選擇是采用從上而下建模，還是從下而上建模。因此在建模時一定要有層次。 模型瀏覽器可以用來觀看模型的層次結構。（菜單view/m

25、odel browser）。2 封裝子系統(tǒng)功能用戶自定義該子系統(tǒng)的圖表和設置參數(shù)對話框。House房屋形狀的圖標就是封裝后的結果。還可以做其他嘗試，體會simulink的特點。如：把溫度設為80（預期的室內(nèi)溫度），觀察仿真結果的變化。標簽為daily temperature variation的正弦模塊是設置日溫度變化，試改變幅度值參數(shù)，觀察仿真結果的變化。模塊庫（simulink library）是存放simulink模塊的地方，可以在瀏覽器窗口中選擇所需的模塊。Simulink這一標題必定會有，是基本模塊庫，其他例如通信模塊庫（communication blockset），數(shù)字信號處理（

26、DSP blockset）等，要安裝了相應的工具箱才會有?？梢詥螕裘總€標題前的+號，查看庫里的內(nèi)容。Simulink可分為continuous, discrete, function&table(函數(shù)和平臺), math, nonlinear, signal&system, sinks（接收器）, sources（源）等子庫。Sinks里面是一些信號的接收器，如scope(示波器),display(顯示器),XY Graph(XY圖形).二創(chuàng)建一個簡單的模型示范模型的功能是對一個正弦信號進行積分，并顯示積分的結果。圖4-1 mp1模型在模型結構已設計好的前提下，用SIMULIN

27、K建立模型的過程可以概括為：在SIMULINK的模塊庫中找到所需的模塊，并把它們拖曳到模型窗口，將這些模塊排列好，然后用直線將各個模塊連接起來。步驟如下：（1） 啟動SIMULINK模塊庫瀏覽器窗口；（2） 新建一個空白模型：單擊庫瀏覽器工具欄上的空白按鈕。在SIMULINK里，模型是保存在模型文件里的，后綴名為.mdl.（3） 在模塊庫瀏覽器窗口找到所需的模塊，并將模塊拖曳到空白的模型窗口中的適當位置（均在SIMULINK庫）。正弦發(fā)生器（source子庫的sine wave）、積分器（continuous子庫的integrator）、復用器（signals & systems子庫的

28、MUX）、示波器（sinks子庫的scope）。（4） 用直線將模塊連接起來，注意一個模塊的輸入端只能和另一個模塊的輸出端連接。（5） 保存，*.mdl。（6） 運行(simulation/start)、查看結果（雙擊scope）。修改參數(shù)，查看運行的結果。三基本操作1對模塊操作：任務操作選擇一個模塊單擊鼠標選擇多個模塊Shift+單擊鼠標或者用方框包含選擇對象移動模塊拖動模塊復制模塊Ctrl+鼠標左鍵，然后拖動鼠標在模塊間連接鼠標左鍵斷開模塊間連接Shift+拖動2建立模型注釋：在模型空白處雙擊鼠標，然后輸入注釋文字。 刪除注釋：shift+選中注釋，然后按刪除鍵。四子系統(tǒng)創(chuàng)建及封裝對于一個

29、復雜系統(tǒng)，可以創(chuàng)建一些子系統(tǒng)，每個子系統(tǒng)完成系統(tǒng)的部分特定功能，使復雜系統(tǒng)模型框圖更加清晰，一目了然。創(chuàng)建子系統(tǒng)有兩種途徑（1）增加一個子系統(tǒng)模塊到你的模型；（2）組合已存在的模塊建立子系統(tǒng)。1通過子系統(tǒng)模塊來建立子系統(tǒng)下面以最簡單的子系統(tǒng)為例說明子系統(tǒng)模塊創(chuàng)建一般方法和步驟。 一個子系統(tǒng)方程表示為 y = mx + b式中，x為輸入，y為輸出，m, b為常數(shù)。要求用戶能通過對話框方便修改m, b值。 步驟：（1）將Subsystem模塊復制到模型窗口，它的位置在simulink/signals & systems；（2）雙擊Subsystem模型，打開子系統(tǒng)編輯窗口；（3）在Subs

30、ystem窗口下，建立子系統(tǒng)模型。這里參數(shù)m取名為Slope，參數(shù)b取名為Intercept。子系統(tǒng)輸入和輸出要用輸入模塊In和輸出模塊Out。（4）關閉子系統(tǒng)窗口，保存模型。圖4-2 子系統(tǒng)模型2組合已存在的模塊建立子系統(tǒng)（1） 選取要組合成子系統(tǒng)的模塊：正確的方法是用方框包含待選擇對象。（2） 用edit/creat subsystem命令產(chǎn)生子系統(tǒng)。3子系統(tǒng)封裝利用SIMULINK的封裝（MASK）功能，用戶可以創(chuàng)建一個SIMULINK模塊庫中沒有的子系統(tǒng)（Subsystem）模塊，并建立模塊對話框的圖標。這種功能使SIMULINK功能不斷擴展并能滿足各學科領域的計算機仿真問題。同時，也可使復雜系統(tǒng)的框圖模型得到簡化。步驟：（1）建立子系統(tǒng)。如果不封裝，雙擊子系統(tǒng)會看到子系統(tǒng)的內(nèi)部結構。為了改變m和b兩個參數(shù)，必須分別打開Gain和Constant的參數(shù)對話框逐個設置，如果子系統(tǒng)內(nèi)的模塊很多。這項工作是十分煩瑣的。封裝可以簡化這