《建模与仿真及其医学应用》实验讲义天津医科大学生物医学工程系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）该系统与系统561)(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.mfunction y=dfun(t,x)y=sqrt(x)+5; 然后建立m 文件mp2-1%mp2-1[t,x]=ode23('dfun', [0 10] , 1)plot(t,x)结果： t x0 1.00000.0133 1.08030.0800 1.48900.2720 2.72630.5685 4.78001.0356 8.30351.7589 14.34052.7589 23.67783.7589 34.03414.7589 45.32145.7589 57.48156.7589 70.47217.7589 84.26128.7589 98.82309.7589 114.136510.0000 117.93842．对于高阶常微分方程，),...,,()1()(-=n n y yy t f y ,则可以选择一组状态变量)1(21,....,,-===n x y x yx y x ，将原高阶微分方程模型变换成以下的一阶微分方程组形式：()⎪⎪⎩⎪⎪⎨⎧===n n x x x t f xx 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．转移矩阵法的实现：如果连续系统的状态空间模型为：⎩⎨⎧+=+=DuCx y bu Ax x 则其离散状态空间模型为：⎩⎨⎧+=Φ+Φ=+)()()()()()()()1(k Du k Cx k y k u T k x T k x m 其中AT e T =Φ)( 状态转移矩阵（矩阵指数）⎰-=ΦTt 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-2A = [0 1 ; 0 -1]; % Define 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.09520 0.9048phim1 =[ -9/10+exp(-1/10) ][ 1-exp(-1/10) ]phim = 0.00480.09522．采样和信号保持器以及双线性变换法的实现：MATLAB 还提供了通过采样和信号保持器以及双线性变化法将连续系统模型转换为离散时间系统模型的函数C2D ，调用格式为sysd = c2d (sys, Ts, method)其中，sys 为线性连续时间系统；Ts 为采样时间；sysd 为等价的离散时间系统。