《建模与仿真及其医学应用》.docx

1、建模与仿真及其医学应用建模与仿真及其医学应用实 验 讲 义天津医科大学生物医学工程系2004年实验一 系统建模的MATLAB实现一、实验目的：1学习MATLAB基本知识。2掌握数学模型的MATLAB实现：时域模型、状态空间模型和零极点模型。3学习用MATLAB实现系统外部模型到内部模型的转换。4学习用MATLAB实现系统模型的连接：串联、并联、反馈连接。5了解模型降阶的MATLAB实现。二、实验内容1系统的实现、外部模型到内部模型的转换(1) 给定连续系统的传递函数，利用MATLAB建立传递函数模型，微分方程，并转换为状态空间模型。(2)已知某系统的状态方程的系数矩阵为： 利用MATLAB建立

2、状态空间模型，并将其转换为传递函数模型和零极点模型。(3)已知系统的零极点传递函数为，利用MATLAB转换为传递函数模型和状态空间模型。2系统的离散、连接、降阶（1）给定连续系统的传递函数，将该连续系统的传递函数用零阶重构器和一阶重构器转换为离散型传递函数，抽样时间T=1秒。（2）该系统与系统分别串联并联负反馈连接，求出组成的新系统的传递函数模型。（3）将串联组成的新系统进行降阶处理，求出降阶后系统的模型，并用plot图形比较降阶前后系统的阶跃响应。要求：将以上过程用MATLAB编程（M文件）实现，运行输出结果。三、实验说明关于系统建模的主要MATLAB函数1建立传递函数模型：tf函数 ：格式

3、：sys=tf(num,den)num=bm,bm-1,b0 分子多项式系数den=an,an-1,a0 分母多项式系数2建立状态空间模型：ss函数 ：格式：sys=ss(a,b,c,d) %a,b,c,d为状态方程系数矩阵sys=ss(a,b,c,d,T) %产生离散时间状态空间模型3建立零极点模型的函数：zpk格式：sys=zpk(z,p,k)4模型转换函数：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) %i

4、u指定是哪个输入z,p,k=ss2zp(a,b,c,d,iu)num,den=zp2tf(z,p,k)a,b,c,d=zp2ss(z,p,k)5模型的连接串联：sys=series(sys1,sys2)并联：sys=parallel(sys1,sys2)反馈连接：sys=feedback(sys1,sys2,sign)%负反馈时sign可忽略；正反馈时为1。6系统扩展：把若干个子系统组成系统组。格式：sys=append(sys1,sys2,)7模型降阶（1）基于平衡的状态空间实现-balreal格式：sysb=balreal(sys)sysh,g,T,Ti=balreal(sys)sys为原

5、系统，sysb(sysh)为平衡实现系统，g为平衡对角线矩阵，T为状态变换矩阵，Ti是前者的逆矩阵。两种格式的区别：前者只给出原系统的一个平衡的状态空间实现，而后者还给出平衡实现的对角线矩阵g，从中可以看出哪个状态变量该保留，哪个状态变量该删去，从而实现降阶。（2）降阶的实现modred格式：rsys=modred(sys,elim)rsys=modred(sys,elim,mde)rsys=modred(sys,elim,del)强调：这里的sys应是函数balreal()变换的模型，elim为待消去的状态，mde指降阶中保持增益匹配，del 指降阶中不保持增益匹配。8连续系统模型离散化函数

6、：C2DM Conversion of continuous LTI systems to iscrete-time.格式：Ad,Bd,Cd,Dd=C2DM(A,B,C,D,Ts,method)将连续系统状态空间离散系统状态空间method: zoh 零阶重构器 zero order hold foh 一阶重构器 first order holdNUMd,DENd = C2DM(NUM,DEN,Ts,method) 将连续系统传递函数离散系统传递函数G(s) = NUM(s)/DEN(s) to G(z) = NUMd(z)/DENd(z).四、实验报告要求1 整理好经过运行并证明是正确的程序

7、，必要的地方加上注释。2 给出实验的结果。实验二 连续系统的数字仿真一、计算机仿真在计算机支持下进行的现代仿真技术称为计算机仿真。仿真不单纯是对模型的实验，它包括建立模型、仿真运行和分析研究仿真结果，即建模实验分析的全过程。MATLAB提供各种用于系统仿真的函数，用户可以通过m 文件调用指令和函数进行系统仿真，也可以通过Simulink工具箱，进行面向系统结构方框图的系统仿真。这两种方式可解决任意复杂系统的动态仿真问题，前者编辑灵活，而后者直观性强，实现可视化编辑。内容：连续系统仿真：数值积分法、离散相似法离散事件系统仿真SIMULINK动态仿真二、基于数值积分法的连续系统仿真1数值积分法的M

8、ATLAB函数MATLAB的工具箱提供了各种数值积分方法函数：格式：T,Y=solver(F,TSPAN,Yo,OPTIONS)solver为微分方程的求解函数名。F为系统模型文件名，模型为TSPAN=To Tfinal为积分区间，初值终值，Yo为系统输出初始值，即To时刻的初值列向量；OPTIONS设置积分相对允误RelTol和绝对允误AbsTol，缺省时，RelTol=1e-3, AbsTol=1e-6.输出参数T和Y为列向量，T为时刻向量，Y表示不同时刻的函数值。系统模型函数的编写格式是固定的，如果其格式没有按照要求去编写则将得出错误的求解结果，系统模型函数的引导语句为：function

9、 xdot=模型函数名(t,x,附加参数)其中t为时间变量，x为状态变量，xdot为状态变量的导数。如果有附加参数需要传递，则可以列出，中间用逗号分开。solver:ode23 Runge-Kutta法 三阶积分算法、二阶误差估计、变积分步长的低阶算法ode45 Runge-Kutta法，变步长的中等阶次积分算法ode113 变阶的Adams-Bashforth-Moulton，多步长ode15s 改进的Gear法，用于刚性方程的求解。例：求微分方程，先建立一个系统模型文件（m文件函数）dfun.mfunction y=dfun(t,x)y=sqrt(x)+5;然后建立m文件mp2-1%mp2

10、-1t,x=ode23(dfun, 0 10 , 1)plot(t,x)结果： 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对于高阶常微分方程，,则可以选择一组状态变量，将

11、原高阶微分方程模型变换成以下的一阶微分方程组形式：例：可变换成functiom y=vdp_eq(t,x,mu)y=x(2);-mu*(x(1).2-1).*x(2)-x(1)三、基于离散相似法的连续系统仿真所谓离散相似法是首先将连续系统模型离散化，得到等价的或相似的离散化的模型，然后对相似的离散模型进行仿真计算。根据这一原理，首先应将连续时间系统模型转换为等价的离散时间系统模型。连续系统离散化处理是通过转移矩阵法；采样和信号保持器；变换法（如双线性变换）来实现的。1转移矩阵法的实现：如果连续系统的状态空间模型为：则其离散状态空间模型为：其中 状态转移矩阵（矩阵指数）由此可知，利用状态方程离散

12、化时的主要问题是如何计算、。对于一阶、二阶环节，、可以用解析方法求出来，而对于高阶及多输入多输出系统，就要采用数值解法。MATLAB提供了计算矩阵指数的函数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

13、 -1; % Define 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)%将符号运算转换为数值结果：phi =1.0000 0.0952 0 0.9048phim1 = -9/10+exp(-1/10) 1-exp(-1/10) phim = 0.0048 0.09522采样和信号保持器以及双线性变换法的

14、实现：MATLAB还提供了通过采样和信号保持器以及双线性变化法将连续系统模型转换为离散时间系统模型的函数C2D，调用格式为sysd = c2d (sys, Ts, method)其中，sys为线性连续时间系统；Ts为采样时间；sysd为等价的离散时间系统。method为离散化方法，可以选用： zoh 为零阶保持器 foh为一阶保持器 tustion为双线性变换法， prewarp为改进的双线性变换法 matched使连续和离散系统具有匹配的DC增益例：连续系统传递函数，采用一阶采样保持器，采样周期为，求其离散化系统模型，并比较离散前后系统阶跃响应。用MATLAB编写程序： %mp2-3sysc

15、 = tf ( l -1 , 14 5 , td , 0.35 );%time delaysysd = c2d ( sysc, 0.l, foh )step ( sysc, sysd );运行结果：Transfer function:Sampling time: 0.1 离散前后系统阶跃响应比较四、实验内容1 求解方程在不同值的解，=1，；初值, ；=2，；初值, ；=1000，；初值, 。2. 给定系统k/s1/(s+a)u x1 x2 y用转移矩阵法仿真，其中k=2;a=1;T=0.01;x1(0)=0.1;x2(0)=0.7，仿真时间为0.2 秒，求系统的阶跃响应。采用零阶和一阶采样保持

16、器，求其离散化系统模型，给出系统阶跃响应。采用双线性变换法，求其离散化系统模型，给出系统阶跃响应。实验三 最小二乘法及数据拟合建模的回归分析一、实验目的:1掌握用最小二乘建立回归数学模型。2学习通过几个数据拟合的回归分析来判断曲线（直线）拟合的精度，通过回归分析来判断模型建立是否正确。3应用建立的模型进行预测。二、基本原理和方法1建立回归数学模型在进行建模和仿真分析时，人们经常面临用已知系统实测数据应用数学模型描述对应系统，即对数据进行拟合。拟合的目的是寻找给定的曲线（直线），它在某种准则下最佳地拟合数据。最佳拟合要在什么准则下的最佳？以及用什么样的曲线模型去拟合。常用的拟合方法之一是多项式的

17、最小二乘拟合，其准则是最小误差平方和准则，所用的拟合曲线为多项式。本实验在Matlab平台上，以多项式最小二乘拟合为例，掌握回归模型的建立（包括参数估计和模型建立）和用模型进行预测的方法，并学习回归分析的基本方法。2在MATLAB里，用于求解最小二乘多项式拟合问题的函数如下： polyfit 最小二乘拟合 p=polyfit(x,y,n) 对输入数据y的n阶最小二乘拟合多项式p(x)的系数 Y=polyval(p,x) 求多项式的函数值Y 以下是一个多项式拟合的例子。已知 x=0,1,2,3,4,5,6,7,8,9,10 共11个点（自变量），实测数据y=-0.447, 1.978, 3.28

18、, 6.16, 7.08, 7.34, 7.66, 9.56, 9.48, 9.30, 11.2求：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阶的预测方程 的系数 p= b2 b1 b0 z=polyval(p,x); %求预测的y值 （z表示） p2=polyfit(x,y,8) %求8阶的预测方程 z1=polyval(p2,x); plot(x,y,om,x,z,:*rx,z1, :+b)

19、图中：”0” 代表散点图 “+”代表8阶预测方程 “*”代表2阶预测方程 图1 散点图与2阶预测方程3回归模型的检验回归模型的检验是判断数据拟合的好坏即模型建立的正确与否，为建立模型和应用模型提供支持。在MATLAB平台，用于回归检验的语句如下： b,bint,r,rint,stats=regress(y，X，) 其中， 为回归系数 e 随机误差（均值为0，方差） b：回归系数的估计值 bint：回归系数的置信区间 r：残差 rint：残差的置信区间 stats：统计量 R2 F P 用此语句，可以得到回归系数，复相关系数R，方差检验F值，P。rcoplot(r,rint) %打印残差分布图在

20、图中，若残差的置信区间不包含零点，则视为异常点，将其剔除后重新计算。b=leastsq(函数名,b0) % 非线性最小二乘法拟合， b0为初始值s=sqrt(sum(函数名.* 函数名)./（n-2） %计算剩余标准差s=sqrt(sum(y-Y).2)./(n-2) %计算剩余标准差4预测值： 用经过检验的数学模型即可预测数据。即把x代入回归方程对y进行估计，该估计值为。以下用一个例子说明回归模型的检验与预测： 有人研究了黏虫孵化历期平均温度（x,oC）与历期天数（y,天）之间关系，试验资料列入下表，求直线回归方程，并进行检验。x,历期平均温度oC 11.8 14.7 15.6 16.8 1

21、7.1 18.8 19.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）与历期天数（y,天）的散点图从图中可见y与x基本，因此用一元线性回归。X=ones(8,1),x;b,bint,r,rint,stats=regress(y,X)；b,bint,stats,rcoplot(r,rint)运行

22、结果：b = 57.0393 -2.5317bint = 45.9035 68.1750 -3.1851 -1.8782stats = 0.9374 89.8675 0.0001得到以下结果：复相关系数的平方 R2=0.9374, 则R=0.96820.8,F=89.8675 存在极显著的直线回归关系。 p10-4,因此此回归模型y接受。预测历期平均温度为21.05oC时，孵化历期天数为多少天？（天） 图三 残差分布图 （全部观测点）观察图中的残差分布，除第一点外，其余残差的置信区间均包含零点，第一个点应视为异常点，将其剔除后重新计算，可得：b = 45.2899 -1.8863bint =

23、38.5660 52.0138 -2.2670 -1.5057stats = 0.9701 162.2546 0.0001图三 残差分布图 （除掉第一个异常点观测点） 观察图中的残差分布，残差的置信区间均包含零点，置信区间明显缩小，R与F明显增大。重新预测预测历期平均温度为21.05oC时，孵化历期天数 显然，此估计值应更加符合真实数据。三、实验内容：1 用给定的多项式 产生一组数据，再在yi上添加随机干扰，然后用和添加了随机干扰的yi作3次多项式拟合，与原系数比较。如果作2次或4次多项式拟合，结果如何？2给定数据（xi,yi）见下表画出散点图观察二者的关系，试建立合适的回归模型：直线回归方程

24、；二次曲线；对数曲线；线性化的对数曲线等.并作回归分析,得到最佳。得到的最佳回归方程做出散点图和曲线。xi23457810yi106.42109.20109.58109.50110.00109.93110.49xi111415151819yi110.59110.60110.90110.76111.00111.20实验四 Simulink动态仿真SIMULINK是MATLAB重要软件包，用于对动态系统建模和仿真，它适用于连续系统和离散系统，也适用线性系统和非线性系统。它采用系统模块直观地描述系统典型环节，因此可十分方便地建立系统模型而不需要花较多时间编程。正由于这些特点，SIMULINK广泛流行

25、，被认为是最受欢迎的仿真软件。一、SIMULINK的特点启动：在命令窗口输入Simulink或单击“simulink library Browser”按钮或使用菜单命令filenewmodel，均可打开“simulink library Browser”窗口，里面有许多Simulink基本模块，用户可以直接调用这些模块。1支持图形用户界面。在命令窗口键入“thermo”(这是一个房间热力学仿真演示程序)。里面包含许多模块，模块之间用直线相连，组成模型，模型造好后，可以仿真运行，等待结果。（单击start按钮或菜单命令simulinkstart,然后单击thermo plots模块，查看仿真结果

26、室内温度、室外温度和热量曲线，这种和实际示波器输出相似的图形化显示，给用户一个很友好的界面）。2层次性 双击“house”模块，出现“thermo/house”窗口，表示house模块是由窗口右边所示的一些模块连接而成。像house这样由由几个相互关联的模块组合而成的模块在simulink中称为子系统（后面讲）。而这种一个模块又由许多模块组成的特性，这是simulink的层次性。为了和子系统相区别，这里把模型本身称为模型的顶层系统。层次性的好处是建立的模型在结构上非常清晰整齐，一目了然，更重要的是，这一特性使得用户可以选择是采用从上而下建模，还是从下而上建模。因此在建模时一定要有层次。 模型浏

27、览器可以用来观看模型的层次结构。（菜单view/model browser）。2 封装子系统功能用户自定义该子系统的图表和设置参数对话框。House房屋形状的图标就是封装后的结果。还可以做其他尝试，体会simulink的特点。如：把温度设为80（预期的室内温度），观察仿真结果的变化。标签为daily temperature variation的正弦模块是设置日温度变化，试改变幅度值参数，观察仿真结果的变化。模块库（simulink library）是存放simulink模块的地方，可以在浏览器窗口中选择所需的模块。Simulink这一标题必定会有，是基本模块库，其他例如通信模块库（commun

28、ication blockset），数字信号处理（DSP blockset）等，要安装了相应的工具箱才会有。可以单击每个标题前的+号，查看库里的内容。Simulink可分为continuous, discrete, function&table(函数和平台), math, nonlinear, signal&system, sinks（接收器）, sources（源）等子库。Sinks里面是一些信号的接收器，如scope(示波器),display(显示器),XY Graph(XY图形).二创建一个简单的模型示范模型的功能是对一个正弦信号进行积分，并显示积分的结果。Mp1.mdl图4-1 mp1

29、模型在模型结构已设计好的前提下，用SIMULINK建立模型的过程可以概括为：在SIMULINK的模块库中找到所需的模块，并把它们拖曳到模型窗口，将这些模块排列好，然后用直线将各个模块连接起来。步骤如下：（1） 启动SIMULINK模块库浏览器窗口；（2） 新建一个空白模型：单击库浏览器工具栏上的空白按钮。在SIMULINK里，模型是保存在模型文件里的，后缀名为.mdl.（3） 在模块库浏览器窗口找到所需的模块，并将模块拖曳到空白的模型窗口中的适当位置（均在SIMULINK库）。正弦发生器（source子库的sine wave）、积分器（continuous子库的integrator）、复用器（signals & systems子库的MUX）、示波器（sinks子库的scope）。（4） 用直线将模块连接起来，注意一个模块的输入端只能和另一个模块的输出端连接。（5） 保存，*.mdl。（6） 运行(simulation/start)、查看结果（双击scope）。修改参数，查看运行的结果。三基本操作1对模块操作：任务操作选择一个模块单击鼠标选择多个模块Shift+单击鼠标或者用方框包含选择对象移动模块拖动模块