建模与仿真及其医学应用.docx
《建模与仿真及其医学应用.docx》由会员分享，可在线阅读，更多相关《建模与仿真及其医学应用.docx（21页珍藏版）》请在冰豆网上搜索。
建模与仿真及其医学应用
《建模与仿真及其医学应用》
实验讲义
天津医科大学生物医学工程系
2004年
实验一系统建模的MATLAB实现
一、实验目的：
1．学习MATLAB基本知识。
2．掌握数学模型的MATLAB实现：
时域模型、状态空间模型和零极点模型。
3．学习用MATLAB实现系统外部模型到内部模型的转换。
4．学习用MATLAB实现系统模型的连接：
串联、并联、反馈连接。
5．了解模型降阶的MATLAB实现。
二、实验内容
1．系统的实现、外部模型到内部模型的转换
（1）给定连续系统的传递函数
，利用MATLAB建立传递函数模型，微分方程，并转换为状态空间模型。
（2）已知某系统的状态方程的系数矩阵为：
利用MATLAB建立状态空间模型，并将其转换为传递函数模型和零极点模型。
（3）已知系统的零极点传递函数为
，利用MATLAB转换为传递函数模型和状态空间模型。
2．系统的离散、连接、降阶
（1）给定连续系统的传递函数
，将该连续系统的传递函数用零阶重构器和一阶重构器转换为离散型传递函数，抽样时间T=1秒。
（2）该系统与系统
分别①串联②并联③负反馈连接，求出组成的新系统的传递函数模型。
（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,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）]
[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（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．连续系统模型离散化
函数：
C2DMConversionofcontinuousLTIsystemstoiscrete-time.
格式：
①[Ad,Bd,Cd,Dd]=C2DM（A,B,C,D,Ts,'method'）
将连续系统状态空间—离散系统状态空间
'method':
'zoh'零阶重构器zeroorderhold
'foh'一阶重构器firstorderhold
②[NUMd,DENd]=C2DM（NUM,DEN,Ts,'method'）
将连续系统传递函数—离散系统传递函数
G（s）=NUM（s）/DEN（s）toG（z）=NUMd（z）/DENd（z）.
四、实验报告要求
1．整理好经过运行并证明是正确的程序，必要的地方加上注释。
2．给出实验的结果。
实验二连续系统的数字仿真
一、计算机仿真
在计算机支持下进行的现代仿真技术称为计算机仿真。
仿真不单纯是对模型的实验，它包括建立模型、仿真运行和分析研究仿真结果，即建模——实验——分析的全过程。
MATLAB提供各种用于系统仿真的函数，用户可以通过m文件调用指令和函数进行系统仿真，也可以通过Simulink工具箱，进行面向系统结构方框图的系统仿真。
这两种方式可解决任意复杂系统的动态仿真问题，前者编辑灵活，而后者直观性强，实现可视化编辑。
内容：
连续系统仿真：
数值积分法、离散相似法
离散事件系统仿真
SIMULINK动态仿真
二、基于数值积分法的连续系统仿真
1．数值积分法的MATLAB函数
MATLAB的工具箱提供了各种数值积分方法函数：
格式：
[T,Y]=solver（‘F’,TSPAN,Yo,OPTIONS）
solver为微分方程的求解函数名。
F为系统模型文件名，模型为
TSPAN=[ToTfinal]为积分区间，初值—终值，
Yo为系统输出初始值，即To时刻的初值列向量；
OPTIONS设置积分相对允误’RelTol’和绝对允误’AbsTol’，
缺省时，RelTol=1e-3,AbsTol=1e-6.
输出参数T和Y为列向量，T为时刻向量，Y表示不同时刻的函数值。
系统模型函数的编写格式是固定的，如果其格式没有按照要求去编写则将得出错误的求解结果，系统模型函数的引导语句为：
functionxdot=模型函数名（t,x,附加参数）
其中t为时间变量，x为状态变量，xdot为状态变量的导数。
如果有附加参数需要传递，则可以列出，中间用逗号分开。
solver:
ode23Runge-Kutta法三阶积分算法、二阶误差估计、变积分步长的低阶算法
ode45Runge-Kutta法，变步长的中等阶次积分算法
ode113变阶的Adams-Bashforth-Moulton，多步长
ode15s改进的Gear法，用于刚性方程的求解。
例：
求微分方程
，
，
先建立一个系统模型文件（m文件函数）dfun.m
functiony=dfun（t,x）
y=sqrt（x）+5;
然后建立m文件mp2-1
%mp2-1
[t,x]=ode23（'dfun',[010],1）
plot（t,x）
结果：
tx
01.0000
0.01331.0803
0.08001.4890
0.27202.7263
0.56854.7800
1.03568.3035
1.758914.3405
2.758923.6778
3.758934.0341
4.758945.3214
5.758957.4815
6.758970.4721
7.758984.2612
8.758998.8230
9.7589114.1365
10.0000117.9384
2．对于高阶常微分方程，
则可以选择一组状态变量
，将原高阶微分方程模型变换成以下的一阶微分方程组形式：
例：
可变换成
functiomy=vdp_eq（t,x,mu）
y=[x
（2）;-mu*（x
（1）.^2-1）.*x
（2）-x
（1）]
三、基于离散相似法的连续系统仿真
所谓离散相似法是首先将连续系统模型离散化，得到等价的或相似的离散化的模型，然后对相似的离散模型进行仿真计算。
根据这一原理，首先应将连续时间系统模型转换为等价的离散时间系统模型。
连续系统离散化处理是通过①转移矩阵法；②采样和信号保持器；③变换法（如双线性变换）来实现的。
1．转移矩阵法的实现：
如果连续系统的状态空间模型为：
则其离散状态空间模型为：
其中
状态转移矩阵（矩阵指数）
由此可知，利用状态方程离散化时的主要问题是如何计算
、
。
对于一阶、二阶环节，
、
可以用解析方法求出来，而对于高阶及多输入多输出系统，就要采用数值解法。
MATLAB提供了计算矩阵指数的函数——expm，
EXPMMatrixexponential.
EXPM（X）isthematrixexponentialofX.EXPMiscomputedusingascalingandsquaringalgorithmwithaPadeapproximation.
EXPM1,EXPM2andEXPM3arealternativemethods.
例：
，
求
、
。
%mp2-2
A=[01;0-1];%Definesystemmatrices
B=[0;1];
t=0.1
symstau%Definetautobesymbolic
phi=expm（A*t）%Symbolicallycalculatee^（A*t）
phim1=int（expm（A*（t-tau））,tau,0,t）*B
phim=sym2poly（phim1）%将符号运算转换为数值
结果：
phi=1.00000.0952
00.9048
phim1=[-9/10+exp（-1/10）]
[1-exp（-1/10）]
phim=0.0048
0.0952
2．采样和信号保持器以及双线性变换法的实现：
MATLAB还提供了通过采样和信号保持器以及双线性变化法将连续系统模型转换为离散时间系统模型的函数C2D，调用格式为
sysd=c2d（sys,Ts,method）
其中，sys为线性连续时间系统；Ts为采样时间；sysd为等价的离散时间系统。
method为离散化方法，可以选用：
①'zoh'为零阶保持器
②'foh'为一阶保持器
③'tustion'为双线性变换法，
④'prewarp'为改进的双线性变换法
⑤'matched'使连续和离散系统具有匹配的DC增益
例：
连续系统传递函数
，采用一阶采样保持器，采样周期为
，求其离散化系统模型，并比较离散前后系统阶跃响应。
用MATLAB编写程序：
%mp2-3
sysc=tf（[l-1],[145],'td',0.35）;%timedelay
sysd=c2d（sysc,0.l,'foh'）
step（sysc,sysd）;
运行结果：
Transferfunction:
Samplingtime:
0.1
离散前后系统阶跃响应比较
四、实验内容
1．求解方程
在不同
值的解，
①
=1，
；初值
；
②
=2，
；初值
；
③
=1000，
；初值
。
2.给定系统
ux1x2y
①用转移矩阵法仿真，其中k=2;a=1;T=0.01;x1（0）=0.1;x2（0）=0.7，
仿真时间为0.2秒，求系统的阶跃响应。
②采用零阶和一阶采样保持器，求其离散化系统模型，给出系统阶跃响应。
③采用双线性变换法，求其离散化系统模型，给出系统阶跃响应。
实验三最小二乘法及数据拟合建模的回归分析
一、实验目的:
1．掌握用最小二乘建立回归数学模型。
2．学习通过几个数据拟合的回归分析来判断曲线（直线）拟合的精度，通过回归分析来判断模型建立是否正确。
3．应用建立的模型进行预测。
二、基本原理和方法
1．建立回归数学模型
在进行建模和仿真分析时，人们经常面临用已知系统实测数据应用数学模型描述对应系统，即对数据进行拟合。
拟合的目的是寻找给定的曲线（直线），它在某种准则下最佳地拟合数据。
最佳拟合要在什么准则下的最佳？
以及用什么样的曲线模型去拟合。
常用的拟合方法之一是多项式的最小二乘拟合，其准则是最小误差平方和准则，所用的拟合曲线为多项式。
本实验在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.08,7.34,7.66,9.56,9.48,9.30,11.2
求：
2阶的预测方程，并用8阶的预测方程与之比较。
x=linspace（0,1,11）;
y=[-.4471.9783.286.167.087.347.669.569.489.3011.2];
p=polyfit（x,y,2）
%求2阶的预测方程
的系数p=b2b1b0
z=polyval（p,x）;%求预测的y值（z表示
）
p2=polyfit（x,y,8）%求8阶的预测方程
z1=polyval（p2,x）;
plot（x,y,'om',x,z,':
*r'x,z1,':
+b'）
图中：
”0”代表散点图“+”代表8阶预测方程
“*”代表2阶预测方程
图1散点图与2阶预测方程
3．回归模型的检验
回归模型的检验是判断数据拟合的好坏即模型建立的正确与否，为建立模型和应用模型提供支持。
在MATLAB平台，用于回归检验的语句如下：
[b,bint,r,rint,stats]=regress（y，X，α）
其中，
为回归系数
e随机误差（均值为0，方差
）
b：
回归系数的估计值
bint：
回归系数的置信区间
r：
残差
rint：
残差的置信区间
stats：
统计量R2FP
用此语句，可以得到回归系数，复相关系数R，方差检验F值，P。
rcoplot（r,rint）%打印残差分布图
在图中，若残差的置信区间不包含零点，则视为异常点，将其剔除后重新计算。
b=leastsq（‘函数名’,b0）%非线性最小二乘法拟合，b0为初始值
s=sqrt（sum（‘函数名’.*‘函数名’））./（n-2）%计算剩余标准差
s=sqrt（sum（（y-Y）.^2）./（n-2））%计算剩余标准差
4．预测值：
用经过检验的数学模型即可预测数据。
即把x代入回归方程对y进行估计，该估计值为
。
以下用一个例子说明回归模型的检验与预测：
有人研究了黏虫孵化历期平均温度（x,oC）与历期天数（y,天）之间关系，试验资料列入下表，求直线回归方程，并进行检验。
x,历期平均温度oC11.814.715.616.817.118.819.520.4
Y，孵化历期（天）30.117.316.713.611.910.78.36.7
先作出（xi,yi）的散点图
x=[11.814.715.616.817.118.819.520.4]';
y=[30.117.316.713.611.910.78.36.7]'
plot（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.0393
-2.5317
bint=45.903568.1750
-3.1851-1.8782
stats=0.937489.86750.0001
得到以下结果：
复相关系数的平方R2=0.9374,则R=0.9682>0.8,
F=89.8675存在极显著的直线回归关系。
p<10-4,因此此回归模型y接受。
预测历期平均温度为21.05oC时，孵化历期天数为多少天？
（天）
图三残差分布图（全部观测点）
观察图中的残差分布，除第一点外，其余残差的置信区间均包含零点，第一个点应视为异常点，将其剔除后重新计算，可得：
b=45.2899
-1.8863
bint=38.566052.0138
-2.2670-1.5057
stats=0.9701162.25460.0001
图三残差分布图（除掉第一个异常点观测点）
观察图中的残差分布，残差的置信区间均包含零点，置信区间明显缩小，R与F明显增大。
重新预测预测历期平均温度为21.05oC时，孵化历期天数
显然，此估计值应更加符合真实数据。
三、实验内容：
1．用给定的多项式
产生一组数据
，再在yi上添加随机干扰，然后用
和添加了随机干扰的yi作3次多项式拟合，与原系数比较。
如果作2次或4次多项式拟合，结果如何？
2．给定数据（xi,yi）见下表①画出散点图观察二者的关系，②试建立合适的回归模型：
直线回归方程；二次曲线；对数曲线；线性化的对数曲线等.并作回归分析,得到最佳。
③得到的最佳回归方程做出散点图和曲线。
xi
2
3
4
5
7
8
10
yi
106.42
109.20
109.58
109.50
110.00
109.93
110.49
xi
11
14
15
15
18
19
yi
110.59
110.60
110.90
110.76
111.00
111.20
实验四Simulink动态仿真
SIMULINK是MATLAB重要软件包，用于对动态系统建模和仿真，它适用于连续系统和离散系统，也适用线性系统和非线性系统。
它采用系统模块直观地描述系统典型环节，因此可十分方便地建立系统模型而不需要花较多时间编程。
正由于这些特点，SIMULINK广泛流行，被认为是最受欢迎的仿真软件。
一、SIMULINK的特点
启动：
在命令窗口输入Simulink或单击“simulinklibraryBrowser”按钮或使用菜单命令\model，均可打开“simulinklibraryBrowser”窗口，里面有许多Simulink基本模块，用户可以直接调用这些模块。
1．支持图形用户界面。
在命令窗口键入“thermo”（这是一个房间热力学仿真演示程序）。
里面包含许多模块，模块之间用直线相连，组成模型，模型造好后，可以仿真运行，等待结果。
（单击start按钮或菜单命令simulink\start,然后单击thermoplots模块，查看仿真结果—室内温度、室外温度和热量曲线，这种和实际示波器输出相似的图形化显示，给用户一个很友好的界面）。
2．层次性
双击“house”模块，出现“thermo/house”窗口，表示house模块是由窗口右边所示的一些模块连接而成。
像house这样由由几个相互关联的模块组合而成的模块在simulink中称为子系统（后面讲）。
而这种一个模块又由许多模块组成的特性，这是simulink的层次性。
为了和子系统相区别，这里把模型本身称为模型的顶层系统。
层次性的好处是建立的模型在结构上非常清晰整齐，一目了然，更重要的是，这一特性使得用户可以选择是采用从上而下建模，还是从下而上建模。
因此在建模时一定要有层次。
模型浏览器可以用来观看模型的层次结构。
（菜单view/modelbrowser）。
2．封装子系统功能
用户自定义该子系统的图表和设置参数对话框。
House房屋形状的图标就是封装后的结果。
还可以做其他尝试，体会simulink的特点。
如：
把温度设为80（预期的室内温度），观察仿真结果的变化。
标签为dailytemperaturevariation的正弦模块是设置日温度变化，试改变幅度值参数，观察仿真结果的变化。
模块库（simulinklibrary）是存放simulink模块的地方，可以在浏览器窗口中选择所需的模块。
Simulink这一标题必定会有，是基本模块库，其他例如通信模块库（communicationblockset），数字信号处理（DSPblockset）等，要安装了相应的工具箱才会有。
可以单击每个标题前的+号，查看库里的内容。
Simulink可分为continuous,discrete,function&table（函数和平台）,math,nonlinear,signal&system,sinks（接收器）,sources（源）等子库。
Sinks里面是一些信号的接收器，如scope（示波器）,display（显示器）,XYGraph（XY图形）.
二．创建一个简单的模型
示范模型的功能是对一个正弦信号进行积分，并显示积分的结果。
Mp1.mdl
图4-1mp1模型
在模型结构已设计好的前提下，用SIMULINK建立模型的过程可以概括为：
在SIMULINK的模块库中找到所需的模块，并把它们拖曳到模型窗口，将这些模块排列好，然后用直线将各个模块连接起来。
步骤如下：
（1）启动SIMULINK模块库浏览器窗口；
（2）新建一个空白模型：
单击库浏览器工具栏上的空白按钮。
在SIMULINK里，模型是保存在模型文件里的，后缀名为.mdl.
（3）在模块库浏览器窗口找到所需的模块，并将模块拖曳到空白的模型窗口中的适当位置（均在SIMULINK库）。
正弦发生器（source子库的sinewave）、积分器（continuous子库的integrator）、复用器（signals&systems子库的MUX）、示波器（sinks子库的scope）。
（4）用直线将模块连接起来，注意一个模块的输入端只能和另一个模块的输出端连接。
（5）保存，*.mdl。
（6）运行（simulation/start）、查看结果（双击scope）。
修改参数，查看运行的结果。
三．基本操作
1．对模块操作：
任务
操作
选择一个模块
单击鼠标
选择多个模块
Shift+单击鼠标或者用方框包含选择对象
移动模块
拖动模块
复制模块
Ctrl+鼠标左键，然后拖动鼠标
在模块间连接
鼠标
升级会员 