《建模与仿真及其医学应用》(精)

《《建模与仿真及其医学应用》(精)》由会员分享，可在线阅读，更多相关《《建模与仿真及其医学应用》(精)（27页珍藏版）》请在装配图网上搜索。

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函数 ：格式：sys=tf(num

3、,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) %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模型的连接串联：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为原系统，sysb(sys

5、h)为平衡实现系统，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连续系统模型离散化函数：C2DM Conve

6、rsion 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数值积分法的MATLAB函数MATL

8、AB的工具箱提供了各种数值积分方法函数：格式：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 xdot=模型函数名

9、(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-1t,x=ode23

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

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

15、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. 给定系统1/(s+a)k/su 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, 6.16, 7.0

18、8, 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)图中：”0” 代表散点

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

22、393 -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 存在极显著的直线回归关系。 pplot(0,1,0,m+(m0)上面的命令画一条(0,0)到(0,m)的直线，当m0直线平移1以保证所画直线在模块框区域内。此外，还要把drawing coordinate设置为normalized。（5）按OK键封装关闭封装编辑器。五实验内容：1将摄氏温度转换为华氏温度的公式模型（学号_mp1.mdl）将摄氏

23、温度转换为华氏温度的公式：TF为华氏温度，TC为摄氏温度。2创建简单的连续系统模型(学号_mp2.mdl)给定微分方程建模：u(t)是幅度为1频率为1弧度/分的方波。3给定如下系统(1) 创建此系统模型(学号_mp3.mdl)(2) 建立子系统，创建后的系统模型如下：子系统创建后的系统模型4对Van der Pol方程 将给定的微分方程模型建立Simulink模型。（学号_mp4.mdl）提示：第一个方程可以认为是将信号作为一个积分器的输入端，这样积分器的输出端则将成为信号。类似的，信号本身也可以认为是一个积分器的输出，在积分器的输入端信号应该为，在构造该信号时还应该使用信号乘积的处理。此外，结果输出既可以使用输出端口模块（simulink/signals &systems/out1），也可以使用示波器输出。输入适当参数和初始值：mu=1 (mu即)x01=1,x02=-2然后启动仿真命令，仿真结果将赋给MATLAB工作空间中的保留变量tout和yout。给出仿真结果的时间响应曲线。注：积分模块simulink/continuous/integrator乘积(加法)模块simulink/math/product(sum)常数模块simulink/sources/constant增益模块simulink/math/gain