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

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

2、数&=(2：器罷，利用MATLAB 建立传递函数模型，微分方程，并转换为状态空间模型。(2)已知某系统的状态方程的系数矩阵为：_00_d=1001利用 MATLAB 建立状态空间模型，并将其转换为传递函数模型和零极点模型。MATLAB 转换为传递函数模型和状态空间模型。2.系统的离散、连接、降阶给定连续系统的传递函数 a=(2：；8(3：罷，将该连续系统的传递函数用零阶重构器和一阶重构器转换为离散型传递函数，抽样时 001 0a=b=c= 11-2-31 102(3)已知系统的零极点传递函数为G(s)=2(s+1)(s+2)(s+3)(s+4)，利用2间 T=1 秒。(2)该系统与系统

3、 H(s)=1分别串联并联负反馈连s2+6s+5接，求出组成的新系统的传递函数模型。(3) 将串联组成的新系统进行降阶处理，求出降阶后系统的模型，并用 plot 图形比较降阶前后系统的阶跃响应。要求：将以上过程用 MATLAB 编程(M 文件)实现，运行输出结果。三、实验说明关于系统建模的主要 MATLAB 函数 1建立传递函数模型：tf 函数：格式：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

4、,T)%产生离散时间状态空间模型 3建立零极点模型的函数：zpk格式：sys=zpk(z,p,k)4模型转换函数：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模型的连接串联：sys=series(sysl,sys2)并联：sys=parallel(sys1,sys2)反

5、馈连接： sys=feedback(sys1,sys2,sign)%负反馈时 sign 可忽略； 正反馈时为 1。6系统扩展：把若干个子系统组成系统组。格式：sys=append(sysl,sys2,.)7模型降阶(1)基于平衡的状态空间实现-balreal格式：sysb=balreal(sys)sysh,g,T,Ti=balreal(sys)sys为原系统，sysb(sysh)为平衡实现系统，g 为平衡对角线矩阵，T 为状态变换矩阵，Ti是前者的逆矩阵。两种格式的区别：前者只给出原系统的一个平衡的状态空间实现，而后者还给出平衡实现的对角线矩阵 g， 从中可以看出哪个状态变量该保留，哪个状态变

6、量该删去，从而实现降阶。(2)降阶的实现一 modred格式：rsys=modred(sys,elim)rsys=modred(sys,elim,mde)rsys=modred(sys,elim,del)强调：这里的 sys应是函数 balreal()变换的模型，elim为待消去的状态，mde指降阶中保持增益匹配，del指降阶中不保持增益匹配。48连续系统模型离散化函数：C2DMConversionofcontinuousLTIsystemstoiscrete-time.格式：Ad,Bd,Cd,Dd二 C2DM(A,B,C,D,Ts,method)将连续系统状态空间离散系统状态空间method

7、:zoh零阶重构器 zeroorderholdfoh一阶重构器 firstorderholdNUMd,DENd=C2DM(NUM,DEN,Ts,method)将连续系统传递函数离散系统传递函数G(s)=NUM(s)/DEN(s)toG(z)=NUMd(z)/DENd(z).四、实验报告要求四、实验报告要求1整理好经过运行并证明是正确的程序，必要的地方加上注释2给出实验的结果。5实验二连续系统的数字仿真实验二连续系统的数字仿真一、计算机仿真一、计算机仿真在计算机支持下进行的现代仿真技术称为计算机仿真。仿真不单纯是对模型的实验，它包括建立模型、仿真运行和分析研究仿真结果，即建模实验分析的全过程。M

8、ATLAB 提供各种用于系统仿真的函数，用户可以通过 m 文件调用指令和函数进行系统仿真，也可以通过 Simulink 工具箱，进行面向系统结构方框图的系统仿真。这两种方式可解决任意复杂系统的动态仿真问题，前者编辑灵活，而后者直观性强，实现可视化编辑。内容：连续系统仿真：数值积分法、离散相似法离散事件系统仿真SIMULINK 动态仿真二、基于数值积分法的连续系统仿真二、基于数值积分法的连续系统仿真1.数值积分法的 MATLAB 函数MATLAB 的工具箱提供了各种数值积分方法函数：格式：T,Y=solver(F,TSPAN,Yo,OPTIONS)solver为微分方程的求解函数名。F为系统模型

9、文件名，模型为yQ 二f(t,y)TSPAN=ToTfinal为积分区间，初值一终值，Yo 为系统输出初始值，即 To 时刻的初值列向量；OPTIONS 设置积分相对允误RelTol和绝对允误AbsTol， 缺省时，RelTol=1e-3,AbsTol=1e-6.输出参数 T 和 Y 为列向量，T 为时刻向量，Y 表示不同时刻的函6数值。系统模型函数的编写格式是固定的，如果其格式没有按照要求去编写则将得出错误的求解结果，系统模型函数的引导语句为：functionxdot=模型函数名（t,x，附加参数）其中 t为时间变量，x 为状态变量，xdot 为状态变量的导数。如果有附加参数需要传递，则可以

10、列出，中间用逗号分开。solver:ode23Runge-Kutta 法三阶积分算法、二阶误差估计、变积分步长的低阶算法ode45Runge-Kutta 法，变步长的中等阶次积分算法odell3 变阶的 Adams-Bashforth-Moulton，多步长ode15s改进的 Gear法，用于刚性方程的求解。例：求微分方程 x=px+5，0t10，x0=1先建立一个系统模型文件(m 文件函数)dfun.mfunctiony=dfun(t,x)y=sqrt(x)+5;然后建立 m文件 mp2-1%mp2-1t,x=ode23(dfun,010,1)plot(t,x)结果：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),则可以选择一组状态变量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)三、基于离散相似法的连续系统仿真三、基于离散相似法的连续系统仿真所谓离散相似法是首先将连续系统模型离散化，得到等价的或相似的离散化的模型，然后对相似的离散模型进行仿真计算。根据这一原理，首先应将连续时间系统模型转换为等价的离散时间系统模型。连续系统离散化处理是通过转移矩阵法；采样和信号保持器；变换法(如双线性变换)来实现的。1.转移矩阵法的实现：如果连续系

13、统的状态空间模型为：8|x=Ax+buy=Cx+Du则其离散状态空间模型为：x(k+1)二(T)x(k)+O(T T)u(k)my(k)二 Cx(k)+Du(k)其中(T)二eAT状态转移矩阵(矩阵指数)(T)=fe eA(T-t)Bdtm0由此可知，利用状态方程离散化时的主要问题是如何计算(T)、(T)。对于一阶、二阶环节，(T)、(T)可以用解析方法求出mm来，而对于高阶及多输入多输出系统，就要采用数值解法。MATLAB提供了计算矩阵指数的函数 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)%将符号运算转换为数值9结果：phi=1.00000.095200.9048phim1=-9/10+exp(-1/10)1-exp(-1/10)phim=0.00480.09522采样和信号保持器以及双线性变换法的实现：采样和信号保持器以及双线性变换法的实现：MATLAB 还提供了通过采样和信号保持器以及双线性变化法将连续系统模型转换为离散时间系统模型的函数 C2D，调用格式为sysd=c2d(sys,Ts,method)其中，sys为线性连续时间系统；Ts为采样时间；sysd 为等价的离散时间系统。method 为离散化方法，可以选用：1zoh为零阶保持器2foh为一阶保持器4

16、prewarp 为改进的双线性变换法5matched使连续和离散系统具有匹配的 DC 增益器，采样周期为TS二0.1s，求其离散化系统模型，并比较离散前后系统阶跃响应。tustion为双线性变换法,Z=eT1+Ts21-Ts2s例：连续系统传递函数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);运行结果：Transferfunction:0.0115z人3+0.0456z人2-0.0562z-0

17、.009104zA6-1.629zA5+0.6703zA4Samplingtime:0.1离散前后系统阶跃响应比较11四、实验内容四、实验内容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 存在极显著的直线回归关系。pplot（0,1,0,m+（m0）上面的命令画一条（0,0）到（0,m）的直线，当 m0 直线平移 1 以保证所画直线在模块框区域内。此外，还要把 drawingcoordinate 设置为norma

18、lized。（5）按 OK 键封装关闭封装编辑器。五实验内容：五实验内容：1.将摄氏温度转换为华氏温度的公式模型（学号mp1.mdl）9将摄氏温度转换为华氏温度的公式：T=9T+32F5cT 为华氏温度，T 为摄氏温度。FC252. 创建简单的连续系统模型(学号mp2.mdl)给定微分方程建模：x(t)=-2x(t)+u(t)u(t)是幅度为 1 频率为 1 弧度/分的方波。3. 给定如下系统(2)建立子系统，创建后的系统模型如下：4对 VanderPol 方程x=x,X二一卩(x2-1)x-x1221211将给定的微分方程模型建立 Simulink 模型。 (学号 mp4.mdl)提示：第一个方程可以认为是将x(t)信号作为一个积分器的输入端，IniOut2+SineWaveMuxScopeTrjnsftrFenTriniftrFcn1Sum(1)创建此系统模型(学子系统创建后的系统模型Out1262这样积分器的输出端则将成为x(t)信号。类似的，x(t)信号本身也12可以认为是一个积分器的输出，在积分器的输入端信号应该为叽Xj1)X2-X!，在构造该信号时还应该使用信号乘积的处理