建模与仿真及其医学应用精

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

实验讲义

天津医科大学生物医学工程系

2004年

实验一 系统建模的MATLAB 实现

一、实验目的：

1．学习MATLAB 基本知识。

2．掌握数学模型的MATLAB 实现：时域模型、状态空间模型和零极点模型。

3．学习用MATLAB 实现系统外部模型到内部模型的转换。

4．学习用MATLAB 实现系统模型的连接：串联、并联、反馈连接。

5．了解模型降阶的MATLAB 实现。

二、实验内容

1．系统的实现、外部模型到内部模型的转换

(1)给定连续系统的传递函数)

1343)(32()52)(8()(22++++++=s s s s s s s G ，利用MATLAB 建立传递函数模型，微分方程，并转换为状态空间模型。

(2)已知某系统的状态方程的系数矩阵为：

??????--=3210a ??????=1101b ??????????=210011c ????

??????=100010d 利用MATLAB 建立状态空间模型，并将其转换为传递函数模型和零极点模型。

(3)已知系统的零极点传递函数为)4)(3)(2()1(2)(++++=s s s s s G ，利用MATLAB 转换为传递函数模型和状态空间模型。

2．系统的离散、连接、降阶

（1）给定连续系统的传递函数)

1343)(32()52)(8()(22++++++=s s s s s s s G ，将该连

续系统的传递函数用零阶重构器和一阶重构器转换为离散型传递函数，抽样时间T=1秒。

（2）该系统与系统5

61)(2++=s s s H 分别①串联②并联③负反馈连接，求出组成的新系统的传递函数模型。

（3）将串联组成的新系统进行降阶处理，求出降阶后系统的模型，并用plot 图形比较降阶前后系统的阶跃响应。

要求：将以上过程用MATLAB 编程（M 文件）实现，运行输出结果。

三、实验说明—关于系统建模的主要MATLAB 函数

1．建立传递函数模型：tf 函数 ：

格式：sys=tf(num,den)

num=[b m ,b m-1,……,b 0] 分子多项式系数

den=[a n ,a n-1,……,a 0] 分母多项式系数

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,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．连续系统模型离散化

函数：

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 hold

②[NUMd,DENd] = C2DM(NUM,DEN,Ts,'method')

将连续系统传递函数—离散系统传递函数

G(s) = NUM(s)/DEN(s) to G(z) = NUMd(z)/DENd(z).

四、实验报告要求

1．整理好经过运行并证明是正确的程序，必要的地方加上注释。2．给出实验的结果。

实验二 连续系统的数字仿真

一、计算机仿真

在计算机支持下进行的现代仿真技术称为计算机仿真。

仿真不单纯是对模型的实验，它包括建立模型、仿真运行和分析研究仿真结果，即建模——实验——分析的全过程。

MATLAB 提供各种用于系统仿真的函数，用户可以通过m 文件调用指令和函数进行系统仿真，也可以通过Simulink 工具箱，进行面向系统结构方框图的系统仿真。这两种方式可解决任意复杂系统的动态仿真问题，前者编辑灵活，而后者直观性强，实现可视化编辑。 内容：连续系统仿真：数值积分法、离散相似法

离散事件系统仿真

SIMULINK 动态仿真

二、基于数值积分法的连续系统仿真

1．数值积分法的MATLAB 函数

MATLAB 的工具箱提供了各种数值积分方法函数：

格式：[T,Y]=solver(‘F’,TSPAN,Yo,OPTIONS)

solver 为微分方程的求解函数名。

F 为系统模型文件名，模型为()()y t f t y ,'=

TSPAN=[To Tfinal]为积分区间，初值—终值，

Yo 为系统输出初始值，即To 时刻的初值列向量；

OPTIONS 设置积分相对允误’RelTol’和绝对允误’AbsTol’，

缺省时，RelTol=1e-3, AbsTol=1e-6.

输出参数T 和Y 为列向量，T 为时刻向量，Y 表示不同时刻的函数值。

系统模型函数的编写格式是固定的，如果其格式没有按照要求去编写则将得出错误的求解结果，系统模型函数的引导语句为：

function xdot=模型函数名(t,x,附加参数)

其中t 为时间变量，x 为状态变量，xdot 为状态变量的导数。如果有附加参数需要传递，则可以列出，中间用逗号分开。

solver:

ode23 Runge-Kutta 法 三阶积分算法、二阶误差估计、变积分步长

的低阶算法

ode45 Runge-Kutta 法，变步长的中等阶次积分算法

ode113 变阶的Adams-Bashforth-Moulton ，多步长

ode15s 改进的Gear 法，用于刚性方程的求解。 例：求微分方程5+=?

x x ，100≤≤t ，10=x

先建立一个系统模型文件（m 文件函数）dfun.m

function y=dfun(t,x)

y=sqrt(x)+5; 然后建立m 文件mp2-1

%mp2-1

[t,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 2

3.6778

3.7589 3

4.0341

4.7589 4

5.3214

5.7589 57.4815

6.7589 70.4721

7.7589 84.2612

8.7589 98.8230

9.7589 114.1365

10.0000 117.9384

2．对于高阶常微分方程，),...,,()1()(-=n n y y

y t f y ,则可以选择一组状态变量)1(21,....,,-===n x y x y

x y x ，将原高阶微分方程模型变换成以下的一阶微分方程组形式：

()???????===n n x x x t f x

x x x x ,...,,213221 例：0)1(2=+'-+''y y y y μ

可变换成1221221)1(,x x x x x x

---==μ functiom y=vdp_eq(t,x,mu)

y=[x(2);-mu*(x(1).^2-1).*x(2)-x(1)]

三、基于离散相似法的连续系统仿真

所谓离散相似法是首先将连续系统模型离散化，得到等价的或相似的离散化的模型，然后对相似的离散模型进行仿真计算。

根据这一原理，首先应将连续时间系统模型转换为等价的离散时间系统模型。连续系统离散化处理是通过①转移矩阵法；②采样和信号保持器；③变换法（如双线性变换）来实现的。

1．转移矩阵法的实现：

如果连续系统的状态空间模型为：

???+=+=Du

Cx y bu Ax x 则其离散状态空间模型为：

?

??+=Φ+Φ=+)()()()()()()()1(k Du k Cx k y k u T k x T k x m 其中AT e T =Φ)( 状态转移矩阵（矩阵指数）

?-=ΦT

t T A m Bdt e T 0

)()(

由此可知，利用状态方程离散化时的主要问题是如何计算)(T Φ、)(T m Φ。对于一阶、二阶环节，)(T Φ、)(T m Φ可以用解析方法求出来，而对于高阶及多输入多输出系统，就要采用数值解法。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.

例：Bu Ax x

+= ??????-=1010A ，??

????=10b 求)(T Φ、)(T m Φ。

%mp2-2

A = [0 1 ; 0 -1]; % Define system matrices

B = [0 ; 1];

t=0.1

syms tau % Define tau to be symbolic

phi = expm(A*t) % Symbolically calculate e^(A*t)

phim1= int(expm(A*(t-tau)),tau,0,t)*B

phim=sym2poly(phim1)%将符号运算转换为数值

结果：phi =1.0000 0.0952

0 0.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'为双线性变换法，

2s T 12

s T 1e Z s s s T s -+==

④ 'prewarp'为改进的双线性变换法 ⑤ 'matched'使连续和离散系统具有匹配的DC 增益 例：连续系统传递函数

s

e s s s s H 35.02541)(-++-=，采用一阶采样保持器，采样周期为

s T s 1.0=，求其离散化系统模型，并比较离散前后系

统阶跃响应。

用MATLAB 编写程序：

%mp2-3

sysc = tf ( [l -1] , [14 5] , 'td' , 0.35 );%time delay

sysd = c2d ( sysc, 0.l, 'foh' )

step ( sysc, sysd );

运行结果： Transfer function:

4^6703.05^629.16^009104

.00562.02^0456.03^0115.0z z z z z z +---+ Sampling time: 0.1

离散前后系统阶跃响应比较

四、实验内容

1．求解方程0)1(2=+'-+''y y y y μ在不同μ值的解，

①μ=1，200≤≤t ；初值2.0)0(-=y , 7.0)0(-='y ；

②μ=2，200≤≤t ；初值2.0)0(-=y , 7.0)0(-='y ；

③μ=1000，30000≤≤t ；初值2)0(=y , 0)0(='y 。

2. 给定系统

y ①用转移矩阵法仿真，其中k=2;a=1;T=0.01;x 1(0)=0.1;x 2(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 )1n (p x )n (p x )2(p x )1(p Y 1n n +++++=-

以下是一个多项式拟合的例子。

已知 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=[-.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阶的预测方程 2210x b x b b y ++= 的系数 p= b 2 b 1 b 0

z=polyval(p,x); %求预测的y 值 （z 表示y ）

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 ，α)

其中，01122m m y x x x e =β+β+β++β+

)m ,2,1,0i (i =β为回归系数

e 随机误差（均值为0，方差2σ）

b ：回归系数的估计值 ?b =β

bint ：回归系数的置信区间

r ：残差 y y r -=

rint ：残差的置信区间

stats ：统计量 R 2 F P

用此语句，可以得到回归系数，复相关系数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 进行估计，该估计值为?y。

以下用一个例子说明回归模型的检验与预测：

有人研究了黏虫孵化历期平均温度（x,o C）与历期天数（y,天）之间关系，试验资料列入下表，求直线回归方程，并进行检验。

先作出(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.7]'

plot(x,y,'r+')

图2 历期平均温度（x,o C）与历期天数（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.9035 68.1750

-3.1851 -1.8782

stats = 0.9374 89.8675 0.0001

得到以下结果：

复相关系数的平方 R 2=0.9374, 则R=0.9682>0.8,

F=89.8675 存在极显著的直线回归关系。 p<10-4,因此此回归模型y 接受。

预测历期平均温度为21.05o C 时，孵化历期天数为多少天？ 7470.305.21)5317.2(0393.57bx a y =?-+=+= （天）

图三 残差分布图 （全部观测点）

观察图中的残差分布，除第一点外，其余残差的置信区间均包含零点，第一个点应视为异常点，将其剔除后重新计算，可得：

b = 45.2899

-1.8863

bint = 38.5660 52.0138

-2.2670 -1.5057

stats = 0.9701 162.2546 0.0001

图三残差分布图（除掉第一个异常点观测点）

观察图中的残差分布，残差的置信区间均包含零点，置信区间明显缩小，R与F明显增大。重新预测预测历期平均温度为21.05o C时，孵化历期天数

?

y=

a

bx

-

45

=

+

=

+

)

.

05

.5

5832

.

21

2899

8863

(

.1

显然，此估计值应更加符合真实数据。

三、实验内容：

1．用给定的多项式 3x 5x 6x y 23-+-=产生一组数据

)n ,,2,1i ,y ,x (i i =，再在y i 上添加随机干扰，然后用i x 和添加了随机干扰的y i 作3次多项式拟合，与原系数比较。如果作2次或4次多项式拟合，结果如何？

2．给定数据（x i ,y i ）见下表①画出散点图观察二者的关系，②试建

立合适的回归模型：直线回归方程；二次曲线；对数曲线；线性化的对数曲线等.并作回归分析,得到最佳。③得到的最佳回归方程做出散点图和曲线。

实验四 Simulink动态仿真

SIMULINK是MATLAB重要软件包，用于对动态系统建模和仿

．．．．．．．．真．，它适用于连续系统和离散系统，也适用线性系统和非线性系统。它采用系统模块直观地描述系统典型环节，因此可十分方便地建立系统模型而不需要花较多时间编程。正由于这些特点，SIMULINK广泛流行，被认为是最受欢迎的仿真软件。

一、SIMULINK的特点

启动：在命令窗口输入Simulink或单击“simulink library Browser”按钮或使用菜单命令file\new\model，均可打开“simulink library Browser”窗口，里面有许多Simulink基本模块，用户可以直接调用这些模块。

1．支持图形用户界面。

在命令窗口键入“thermo”(这是一个房间热力学仿真演示程序)。

里面包含许多模块，模块之间用直线相连，组成模型，模型造好后，可以仿真运行，等待结果。（单击start按钮或菜单命令simulink\start,然后单击thermo plots模块，查看仿真结果—室内温度、室外温度和热量曲线，这种和实际示波器输出相似的图形化显示，给用户一个很友好的界面）。

2．层次性

双击“house”模块，出现“thermo/house”窗口，表示house 模块是由窗口右边所示的一些模块连接而成。像house这样由由几个相互关联的模块组合而成的模块在simulink中称为子系统（后面讲）。而这种一个模块又由许多模块组成的特性，这是simulink的层次性。为了和子系统相区别，这里把模型本身称为模型的顶层系统。