matlab在建模仿真中的应用
合集下载
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
2013/5/13 19
u (u1 , u2 ur )T y ( y1 , y2 ym )T
2013/5/13
20
动态模型
结构图
可控标准型
线性系统
2013/5/13
21
2013/5/13
22
动态模型
例：线性微分方程求解
1898年德国生理学家 Frank提出了弹性腔模型（ Windkessel circulatory model），在心舒张期血 管的压力可被描述为：
d C E k1C B ( k 2 k3 )C E k 4 C M dt d C M k 3C E k 4 C M dt
CB是系统输入 k1,k2,k3,k4是参数，CE，CM是输出
2013/5/13
27
2013/5/13
28
二、求解过程—求解线性微分方程组
y1=subs(c.cm) y2=subs(c.ce)
2013/5/13 15
动态模型的表达形式
生理系统 动态模型
连续时间 模型
离散时间 模型
微分方程
传递函数
状态方程
结构图
差分方程
系统函数
状态方程
结构图
2013/5/13
16
动态模型
微分方程
动态模型
传递函数
d y dy du 5 6 y 2 8u dt 2 dt dt
2
Y (s) 2s 8 U ( s ) s 2 5s 6
r=dsolve('C*Dp+1/R*p=0','p(t0)=P0','t')
r =P0/exp(-1/R/C*t0)*exp(-1/R/C*t) P0*exp(-(-t0+t)/R/C)
simple(r)
2013/5/13
25
2013/5/13
26
二、求解过程—求解线性微分方程组
FDG in blood ( V1，CB ) k1 k2 FDG in tissue (V2， CE ) k3 k4 FDG-6-P in tissue (V2， CM )
2013/5/13
17
2013/5/13
18
wk.baidu.com
3
动态模型
状态方程
动态模型
将一个高阶微分方程表示为一组一阶方程，利用矩阵求解
Ax Bu x y Cx Du
状态方程 输出方程
x ( x1 , x2 xn )T
状态变量 输入向量 输出向量
Ax Bu x y Cx Du
x1 x x 2 xn
b1 b b 2 bn
2013/5/13
5
2013/5/13
6
1
例：线性稳态模型
2.稳态线性模型
求解 线性代数方程组数值求解理论 非迭代法——高斯法， 迭代法——雅可比法 Matlab程序实现 [R，j]=rref([A,b]) X=A\b 在生物技术和代谢工程中需要对微生物细胞的培养和生长进行 设计，也就是要设计生物反应器。生物反应器利用生化底物为 微生物生长提供基本化学元素。以下是一个表示微生物需氧生 长过程的典型化学反应（反应模型）
其中： b是摩擦系数 Re是雷诺数
2013/5/13 13
2013/5/13
14
二、动态模型
真实的生理系统是动态的，它们的状态一直随时间或 空间的不同而发生变化，当系统的暂态过程结束后， 模型即变为稳态模型了。 一个自变量常微分方程 多个自变量偏微分方程 方法： 欧拉法 龙格-库塔法
常微分方程求解的符号函数 r = dsolve('eq1,eq2,...', 'cond1,cond2,...', 'v') eq是带求解的微分方程 cond是求解的初始条件 v是自变量 eq的表示 D表示对自变量的微分求导，Dy，D2y，D5y
C
dP(t ) P (t ) 0, t0 t T dt R
C dP(t ) P(t ) 0, t0 t T dt R
( t t0 ) RC
可利常微分方程求解
P P0 e
非线性系统
可利用matlab函数
2013/5/13 23 2013/5/13 24
4
二、求解过程—求解线性微分方程
dsolve函数
呼吸熵RQ ——每消耗单位摩尔氧气产 生的二氧化碳摩尔数
二、求解过程—求解稳态方程
1 0.8 0 0 3 1.7 2 0 0.4 0 1 0.15 0 a 2 2 b 6 1 c 0.6 0 d 0
2013/5/13
二、求解过程—求解线性微分方程组
已知输入
已知参数
(y(10)-y(13))/(c(7)*rp(3))];
33
2013/5/13
34
二、求解过程—求解线性微分方程组
global c c=[0.65,1.4,0.003,0.15,0.003,0.0004,0.0002]; global l l=[6e-4,0.0075,0.048,0.014,0.09,0.12,0.4]; global r r=[0.02,0.1,0.36,0.2,0.5,0.6,0.9]; global rp rp=[300,24,27,0.8]; t_final=20; h_opt=odeset; y0=[0;0;0;0;0;0;0;0;0;0;0;0;0]; [t,y]=ode15s(@arteryeqs,[0,t_final],y0,h_opt); plot(t,y(:,4));axis([19.2,20,0,600]);title('Q1') plot(t,y(:,6));axis([19.2,20,-10,20]);title('Q3')
目的：描述生理系统的稳态行为 模型形式：稳定状态下参量的方程 线性：线性代数方程组 非线性：非线性代数方程组
矩阵形式
方法：建立守恒定律的数学公式
Ax b
a11 a A 21 an1
a12 a1n a22 a2 n an 2 ann
9
2013/5/13
A=[0.8,0,1,0;0,-3,1.7,2;-2,0,0.4,1;0,1,0.15,0]; b=[2;6;-0.6;0]; x=inv(A)*b % x=A\b % r=rref([A,b])
2.稳态非线性模型
非线性代数方程 包含变量的非线性函数或变量的高次项
方程组：dy/dt=f(y,x) X用连续函数描述样条插值 定义全局变量
集中参数模型的控制方程31
2013/5/13
32
function dy=arteryeqs(t,y) global c; global l; global r; global rp; T=0.8; q=[0,100,240,450…,0]; ts=[0,0.0093,0.0186, …0.42,0.5,0.6,0.8]; x=spline(ts,q,mod(t,T)); dy=[(x-y(4)-y(6)-y(7))/c(1); (y(6)-y(8)-y(9))/c(3); (y(9)-y(10))/c(6); (y(1)-y(5)*rp(4)-y(4)*r(2))/l(2); (y(4)-y(5))/(rp(4)*c(2)); (y(1)-y(2)-y(6)*r(3))/l(3); (y(1)-y(11)*rp(1)-y(7)*r(4))/l(4); (y(2)-y(12)*rp(2)-y(8)*r(5))/l(5); (y(2)-y(3)-y(9)*r(6))/l(6); (y(3)-y(13)*rp(3)-y(10)*r(7))/l(7); (y(7)-y(11))/(c(4)*rp(1)); (y(8)-y(12))/(c(5)*rp(2));
1 4.07 ln(Re b ) 0.6 b
sym x solve('1/sqrt(x)-4.07*log(2e5*sqrt(x))+0.6') ans = 0.00083 roots（p） fzero（@func,x0） 编写牛顿迭代法函数
1 4.07 ln(Re b ) 0.6 b
C2 H 5OH aO2 bNH 3 cCH1.7 N 0.15O0.4 dH 2O eCO2
呼吸熵RQ——每消耗单位摩尔氧气产生的二氧化碳摩尔数
C : 2 c ( RQ )a H : 6 3b 1.7c 2d O :1 2a 0.4c d 2( RQ)a N : b 0.15c
线性模型 非线性模型
如何估计参数
动态模型
空间无关的 空间相关的 线性的 非线性
2013/5/13
3
2013/5/13
4
一、稳态模型
1.稳态问题
2、稳态线性模型
线性代数方程组
a11 x1 a12 x2 b1 a21 x1 a22 x2 b2 an1 x1 an 2 x2 bn
RQ通过实验可测为已知量=0.8 用matlab求解 高斯消元法 Guass-Jordan消元法 迭代法
2013/5/13 10
1 RQ 0 0 3 1.7 2 0 0.4 0 1 0.15
0 a 2 2 b 6 1 c 1 2 RQ 0 d 0
例：线性微分方程组求解
颈内动脉
主动脉弓其他分支
左颈总动脉
颈外动脉
颈外动脉后分支
升主动脉
降主动脉
2013/5/13
29
2013/5/13
30
5
定义一个函数来描述这一方程组 输入变量：Q0x 输出矩阵 y=[P1,P3,P6,Q1,P2,Q3,Q4,Q5,Q6,Q7,P4,P5,P7]
带颈动脉分岔血管分支的体动脉系统集中参数模型
2013/5/13
7
2013/5/13
8
二、求解过程—求解稳态方程
生物反应器模型
C2 H 5OH aO2 bNH 3 cCH1.7 N 0.15O0.4 dH 2O eCO2
C : 2 c ( RQ )a H : 6 3b 1.7c 2d O :1 2a 0.4c d 2( RQ)a N : b 0.15c
二、求解过程—求解线性微分方程组
c=dsolve('Dce=k1*cb(k2+k3)*ce+k4*cm','Dcm=k3*cek4*cm','ce(0)=0','cm(0)=0','t'); ？已知输入求输出
t=0:120; cb=1*exp(-0.14*t); k1=1.96;k2=1.022;k3=0.149;k4=0.01;
2013/5/13
35
2013/5/13
36
6
例：非线性微分方程组的求解
酶促反应求解 假设有一种酶E，通过与底物S形成中间产物ES来催化底 物转化成产物P，其反应式为： 初始条件
常数参量
描述该酶促反应动力学过程的微分方程组如下：
生理系统建模与仿真中的问题
如何建模
生物系统建模与仿真
matlab在建模仿真中的应用
——用数学公式描述生理系统
推导法 集总参数法 系统辨识法
2013/5/13
1 2013/5/13 2
生理系统建模与仿真中的问题
建模与仿真问题数学问题 如何求解模型
模型类型
稳态模型
——已知输入求输出 ——已知输入和输出求方程参数
ln( x) x 2 x sin x 4 x 5x 4 0
3
非线性问题的处理方法： 换元或代换法：非线性线性 迭代法：线性插值法，牛顿法
2013/5/13
11
2013/5/13
12
2
例：稳态非线性模型
生物医学中的流体传输动力学模型，通常血管中的血流是 层流，但当血管某处流速过大时有可能发生湍流，湍流时 要用摩擦系数确定流体的速度变化，管道的摩擦系数可以 由下面的模型描述： Matlab求解非线性模型
u (u1 , u2 ur )T y ( y1 , y2 ym )T
2013/5/13
20
动态模型
结构图
可控标准型
线性系统
2013/5/13
21
2013/5/13
22
动态模型
例：线性微分方程求解
1898年德国生理学家 Frank提出了弹性腔模型（ Windkessel circulatory model），在心舒张期血 管的压力可被描述为：
d C E k1C B ( k 2 k3 )C E k 4 C M dt d C M k 3C E k 4 C M dt
CB是系统输入 k1,k2,k3,k4是参数，CE，CM是输出
2013/5/13
27
2013/5/13
28
二、求解过程—求解线性微分方程组
y1=subs(c.cm) y2=subs(c.ce)
2013/5/13 15
动态模型的表达形式
生理系统 动态模型
连续时间 模型
离散时间 模型
微分方程
传递函数
状态方程
结构图
差分方程
系统函数
状态方程
结构图
2013/5/13
16
动态模型
微分方程
动态模型
传递函数
d y dy du 5 6 y 2 8u dt 2 dt dt
2
Y (s) 2s 8 U ( s ) s 2 5s 6
r=dsolve('C*Dp+1/R*p=0','p(t0)=P0','t')
r =P0/exp(-1/R/C*t0)*exp(-1/R/C*t) P0*exp(-(-t0+t)/R/C)
simple(r)
2013/5/13
25
2013/5/13
26
二、求解过程—求解线性微分方程组
FDG in blood ( V1，CB ) k1 k2 FDG in tissue (V2， CE ) k3 k4 FDG-6-P in tissue (V2， CM )
2013/5/13
17
2013/5/13
18
wk.baidu.com
3
动态模型
状态方程
动态模型
将一个高阶微分方程表示为一组一阶方程，利用矩阵求解
Ax Bu x y Cx Du
状态方程 输出方程
x ( x1 , x2 xn )T
状态变量 输入向量 输出向量
Ax Bu x y Cx Du
x1 x x 2 xn
b1 b b 2 bn
2013/5/13
5
2013/5/13
6
1
例：线性稳态模型
2.稳态线性模型
求解 线性代数方程组数值求解理论 非迭代法——高斯法， 迭代法——雅可比法 Matlab程序实现 [R，j]=rref([A,b]) X=A\b 在生物技术和代谢工程中需要对微生物细胞的培养和生长进行 设计，也就是要设计生物反应器。生物反应器利用生化底物为 微生物生长提供基本化学元素。以下是一个表示微生物需氧生 长过程的典型化学反应（反应模型）
其中： b是摩擦系数 Re是雷诺数
2013/5/13 13
2013/5/13
14
二、动态模型
真实的生理系统是动态的，它们的状态一直随时间或 空间的不同而发生变化，当系统的暂态过程结束后， 模型即变为稳态模型了。 一个自变量常微分方程 多个自变量偏微分方程 方法： 欧拉法 龙格-库塔法
常微分方程求解的符号函数 r = dsolve('eq1,eq2,...', 'cond1,cond2,...', 'v') eq是带求解的微分方程 cond是求解的初始条件 v是自变量 eq的表示 D表示对自变量的微分求导，Dy，D2y，D5y
C
dP(t ) P (t ) 0, t0 t T dt R
C dP(t ) P(t ) 0, t0 t T dt R
( t t0 ) RC
可利常微分方程求解
P P0 e
非线性系统
可利用matlab函数
2013/5/13 23 2013/5/13 24
4
二、求解过程—求解线性微分方程
dsolve函数
呼吸熵RQ ——每消耗单位摩尔氧气产 生的二氧化碳摩尔数
二、求解过程—求解稳态方程
1 0.8 0 0 3 1.7 2 0 0.4 0 1 0.15 0 a 2 2 b 6 1 c 0.6 0 d 0
2013/5/13
二、求解过程—求解线性微分方程组
已知输入
已知参数
(y(10)-y(13))/(c(7)*rp(3))];
33
2013/5/13
34
二、求解过程—求解线性微分方程组
global c c=[0.65,1.4,0.003,0.15,0.003,0.0004,0.0002]; global l l=[6e-4,0.0075,0.048,0.014,0.09,0.12,0.4]; global r r=[0.02,0.1,0.36,0.2,0.5,0.6,0.9]; global rp rp=[300,24,27,0.8]; t_final=20; h_opt=odeset; y0=[0;0;0;0;0;0;0;0;0;0;0;0;0]; [t,y]=ode15s(@arteryeqs,[0,t_final],y0,h_opt); plot(t,y(:,4));axis([19.2,20,0,600]);title('Q1') plot(t,y(:,6));axis([19.2,20,-10,20]);title('Q3')
目的：描述生理系统的稳态行为 模型形式：稳定状态下参量的方程 线性：线性代数方程组 非线性：非线性代数方程组
矩阵形式
方法：建立守恒定律的数学公式
Ax b
a11 a A 21 an1
a12 a1n a22 a2 n an 2 ann
9
2013/5/13
A=[0.8,0,1,0;0,-3,1.7,2;-2,0,0.4,1;0,1,0.15,0]; b=[2;6;-0.6;0]; x=inv(A)*b % x=A\b % r=rref([A,b])
2.稳态非线性模型
非线性代数方程 包含变量的非线性函数或变量的高次项
方程组：dy/dt=f(y,x) X用连续函数描述样条插值 定义全局变量
集中参数模型的控制方程31
2013/5/13
32
function dy=arteryeqs(t,y) global c; global l; global r; global rp; T=0.8; q=[0,100,240,450…,0]; ts=[0,0.0093,0.0186, …0.42,0.5,0.6,0.8]; x=spline(ts,q,mod(t,T)); dy=[(x-y(4)-y(6)-y(7))/c(1); (y(6)-y(8)-y(9))/c(3); (y(9)-y(10))/c(6); (y(1)-y(5)*rp(4)-y(4)*r(2))/l(2); (y(4)-y(5))/(rp(4)*c(2)); (y(1)-y(2)-y(6)*r(3))/l(3); (y(1)-y(11)*rp(1)-y(7)*r(4))/l(4); (y(2)-y(12)*rp(2)-y(8)*r(5))/l(5); (y(2)-y(3)-y(9)*r(6))/l(6); (y(3)-y(13)*rp(3)-y(10)*r(7))/l(7); (y(7)-y(11))/(c(4)*rp(1)); (y(8)-y(12))/(c(5)*rp(2));
1 4.07 ln(Re b ) 0.6 b
sym x solve('1/sqrt(x)-4.07*log(2e5*sqrt(x))+0.6') ans = 0.00083 roots（p） fzero（@func,x0） 编写牛顿迭代法函数
1 4.07 ln(Re b ) 0.6 b
C2 H 5OH aO2 bNH 3 cCH1.7 N 0.15O0.4 dH 2O eCO2
呼吸熵RQ——每消耗单位摩尔氧气产生的二氧化碳摩尔数
C : 2 c ( RQ )a H : 6 3b 1.7c 2d O :1 2a 0.4c d 2( RQ)a N : b 0.15c
线性模型 非线性模型
如何估计参数
动态模型
空间无关的 空间相关的 线性的 非线性
2013/5/13
3
2013/5/13
4
一、稳态模型
1.稳态问题
2、稳态线性模型
线性代数方程组
a11 x1 a12 x2 b1 a21 x1 a22 x2 b2 an1 x1 an 2 x2 bn
RQ通过实验可测为已知量=0.8 用matlab求解 高斯消元法 Guass-Jordan消元法 迭代法
2013/5/13 10
1 RQ 0 0 3 1.7 2 0 0.4 0 1 0.15
0 a 2 2 b 6 1 c 1 2 RQ 0 d 0
例：线性微分方程组求解
颈内动脉
主动脉弓其他分支
左颈总动脉
颈外动脉
颈外动脉后分支
升主动脉
降主动脉
2013/5/13
29
2013/5/13
30
5
定义一个函数来描述这一方程组 输入变量：Q0x 输出矩阵 y=[P1,P3,P6,Q1,P2,Q3,Q4,Q5,Q6,Q7,P4,P5,P7]
带颈动脉分岔血管分支的体动脉系统集中参数模型
2013/5/13
7
2013/5/13
8
二、求解过程—求解稳态方程
生物反应器模型
C2 H 5OH aO2 bNH 3 cCH1.7 N 0.15O0.4 dH 2O eCO2
C : 2 c ( RQ )a H : 6 3b 1.7c 2d O :1 2a 0.4c d 2( RQ)a N : b 0.15c
二、求解过程—求解线性微分方程组
c=dsolve('Dce=k1*cb(k2+k3)*ce+k4*cm','Dcm=k3*cek4*cm','ce(0)=0','cm(0)=0','t'); ？已知输入求输出
t=0:120; cb=1*exp(-0.14*t); k1=1.96;k2=1.022;k3=0.149;k4=0.01;
2013/5/13
35
2013/5/13
36
6
例：非线性微分方程组的求解
酶促反应求解 假设有一种酶E，通过与底物S形成中间产物ES来催化底 物转化成产物P，其反应式为： 初始条件
常数参量
描述该酶促反应动力学过程的微分方程组如下：
生理系统建模与仿真中的问题
如何建模
生物系统建模与仿真
matlab在建模仿真中的应用
——用数学公式描述生理系统
推导法 集总参数法 系统辨识法
2013/5/13
1 2013/5/13 2
生理系统建模与仿真中的问题
建模与仿真问题数学问题 如何求解模型
模型类型
稳态模型
——已知输入求输出 ——已知输入和输出求方程参数
ln( x) x 2 x sin x 4 x 5x 4 0
3
非线性问题的处理方法： 换元或代换法：非线性线性 迭代法：线性插值法，牛顿法
2013/5/13
11
2013/5/13
12
2
例：稳态非线性模型
生物医学中的流体传输动力学模型，通常血管中的血流是 层流，但当血管某处流速过大时有可能发生湍流，湍流时 要用摩擦系数确定流体的速度变化，管道的摩擦系数可以 由下面的模型描述： Matlab求解非线性模型