系统仿真的MATLAB实现.docx

1、系统仿真的MATLAB实现第七章 系统仿真的 MATLAB 实现由于计算机技术的高速发展， 我们可以借助计算机完成系统的数字仿真。 综前 所述，数字仿真实质上是根据被研究的真实系统的模型， 利用计算机进行实验研究 的一种方法。仿真的主要过程是：建立模型、仿真运行和分析研究仿真结果。仿真 运行就是借助一定的算法，获得系统的有关信息。MATLAB 是一种面向科学与工程计算的高级语言，它集科学计算、自动控制、 信号处理、神经网络和图像处理等学科的处理功能于一体，具有极高的编程效率。 MATLAB 是一个高度集成的系统， MATLAB 提供的 Simulink 是一个用来对动态系 统进行建模、 仿真和

2、分析的软件包， 它支持线性和非线性系统， 能够在连续时间域、 离散时间域或者两者的混合时间域里进行建模， 它同样支持具有多种采样速率的系 统。在过去几年里， Simulink 已经成为数学和工业应用中对动态系统进行建模时使 用得最为广泛的软件包。MATLAB仿真有两种途径：（1） MATLAB可以在SIMULINK 窗口上进行面向 系统结构方框图的系统仿真； （ 2）用户可以在 MATLAB 的 COMMAND 窗口下， 用运行 m 文件，调用指令和各种用于系统仿真的函数，进行系统仿真。这两种方 式可解决任意复杂系统的动态仿真问题，前者编辑灵活，而后者直观性强，实现可 视化编辑。下面介绍在 M

3、ATLAB 上实现几类基本仿真。7.1 计算机仿真的步骤在学习计算机仿真以前，让我们先总结一下计算机仿真的步骤。 计算机仿真，概括地说是一个“建模实验分析”的过程，即仿真不单纯是 对模型的实验， 还包括从建模到实验再到分析的全过程。 因此进行一次完整的计算 机仿真应包括以下步骤：（1）列举并列项目 每一项研究都应从说明问题开始，问题由决策者提供或由熟悉问题的分析者 提供。（2）设置目标及完整的项目计划 目标表示仿真要回答的问题、系统方案的说明。项目计划包括人数、研究费用 以及每一阶段工作所需时间。（3）建立模型和收集数据 模型和实际系统没有必要一一对应， 模型只需描述实际系统的本质或者描述系

4、统中所研究部分的本质。因此，最好从简单的模型开始，然后进一步建立更复杂的 模型。（ 4）编制程序和验证 利用数学公式、逻辑公式和算法等来表示实际系统的内部状态和输入 /输出的关系。建模者必须决定是采用通用语言如 MATLAB 、FORTRAN 、C 还是专用仿真 语言来编制程序。在本教材中，我们选择的是 MATLAB 和其动态仿真工具 Simulink 。（5）确认确认指确定模型是否精确地代表实际系统。 它不是一次完成， 而是比较模型和 实际系统特性的差异，不断对模型进行校正的迭代过程。（6）实验设计 确定仿真的方案、初始化周期的长度、仿真运行的长度以及每次运行的重复 次数。（7）生产性运行和

5、分析通常用于估计被仿真系统设计的性能量度。 利用理论定性分析、 经验定性分析 或系统历史数据定量分析来检验模型的正确性， 利用灵敏度分析等手段来检验模型 的稳定性。8）文件清单和报表结果9）实现图7.1是计算机仿真的程序图图7.1计算机仿真程序流图7.2基于数值积分法的连续系统仿真7.2.1数值积分法的MATLAB实现MATLAB的工具箱提供了各种数值积分方法函数，这些函数是 ODE23、ODE45、ODE113和ODE15S。这些函数均是 m文件，还有一个函数是 odel.C，是 直接用C语言编写的。函数ode23( 是用Runge-Kutta法求解微分方程。它是一种采用三阶积分算法、 二阶

6、误差估计、变积分步长的低阶积分算法，调用格式为T, Y = ode23 ( F, TSPAN, YO, OPTIONS )其中，F为系统模型文件名，模型为y = f( t, y )形式;TSPAN = To TFINAL为积分计算时间，初值为 To,终值为TFINAL ;YO为系统输出初始值；OPTIONS选项积分计算相对允差RelTol和绝对允差AbsTol,当缺省时，Reltol = 1e-3, AbsTol= 1e-6T为计算点时间向量，丫为微分方程的解。函数ode45()也是用Runge-Kutta法求解微分方程，它是变步长的一种中等阶 次积分算法，调用格式为T, Y = ode45

7、( F , TSPAN, YO, OPTIONS )各项含义同上。函数ode113()是变阶的Adams-Bashforth-Moulton，用变阶方法解微分方程，采用多步法，调用格式为T, Y = ode113 ( F, TSPAN, YO, OPTIONS )各项含义同上。函数odel5s()采用改进的Gear法解微分方程，调用格式为T, Y = odel5s ( F, TSPAN, YO, OPTIONS )各项含义同上。MATLAB还提供了解微分方程函数 ode23s。函数ode1.C是用Euler法求解微 分方程。在SIMULINK中，系统仿真有变步距和固定步距两种工作方式。在变步距

8、仿真中，可选用ode45 ode23 ode113 ose15s和 ode23s在固定步距仿真中，可选 用材 ode5、ode4 ode3 ode2和 ode1。ode5 是 5 阶 Runge-Kutta 法，ode4 是 4 阶 Runge-Kutta法, odel 是 Euler法。【例7.1】求微分方程x = = x 5, 0 _ t _ 10, x0 =1。用MATLAB编写程序：建立一个 m函数文件dfun.mfuncion y= dhfn ( t, x )y = sqrt ( 3*x ) + 5 ;在MATLAB COMMAND 窗口下运行% MATLAB PROGRAM 7-1

9、t, x = ode23 ( dfun, 0 10, 1);Plot ( t, x );7.2.2基于数值积分法的连续系统的数字仿真对于一个n阶微分方程表示的连续系统，也可以用具有n个状态变量的状态空间表达式来描述：X 二AX BU丫二 CX DU状态方程X =AX BU实际上是n个一阶微分方程组成的方程组。如果应用四阶Runge-Kutta法进行仿真计算。则该式可以改写为XkXk -(K! 2K2 2K3 K4)6K AX k BU k心二 A(Xk -KJ BU(tk 7)2 2K3 二A(X k 7K2 BU(tk -)2 2K4 =A(Xk ha) BU(t k h)式中，Ki，K2，

10、K3，K4均为n维列向量。系统输出为yk = CX k+ DU( tk )类似地，可将各种数值积分方法用于连续系统的仿真中。【例7.2】用MATLAB编写图7.2所示液压控制系统的仿真程序。图7.2系统方框图用MATLAB编写仿真程序，采用四阶Run ge-Kutta 法:% MATLAB PROGRAM 7-2%调入数据文件% * Input system data * hynat;% In put system fun cti on; ypfu n1 = valve ;ypfu n 二hysys ;% In itializatio n yref = 5000;x0 = 0 0;y0 = 0

11、 0 0;u0 = 0;t0 = 0;tfinal = 1;tsamp = 0.01;h = 0.001;Y1imit = 512;Ilimit = 40;%调入系统模型%初始化% Refere nt value of system output% In itial value of servo valve% In itial value of hydraulic cyli nder% Start time of simulatio n% End time of simulation% Sample period% Simulatio n step sizemax_epoch = fix (

12、tfinal / h ) 1; t = t0;u1 = u0;x = x0;y = y0;%模块 valve 仿真%模块 hysya 仿真%离散控制量计算tout = zeros ( max_epoch, 1 ); uout = zeros ( max_epoch, 1); yd = zeros ( max_epoch, 1 );yout = zeros ( max_epoch, legth ( y ) );i = 1; tout ( i ) = t; uout ( i ) = 1; yd ( i ) = Kf * y ( l ); yout ( i, : ) = y;% The main l

13、oop for i = 1 : max_epoch % Compute output of valvesvl = feval ( ypfun1, t, u, x, av, bv );sv2 = fcval (ypfun1, t+h/2, u, x+h*sv1/2 , av, hv ); sv3 = feval (ypfun1, t+h/2, u, x+h*sv2/2, av, hv ); sv4 = feval (ypfun1, t+h, u, x+h*sv3, av, hv);x = x+h*(svl+2*sv2+2*sv3+sv4)/6; vo = cv*x;% Compute outpu

14、t of cylinderxl = feval (ypfun, t, vo, y, a, b);s2 = feval (ypfun, t+h/2, vo, y+h*s1/2, a, h);s3 = feval (ypfun, t+h/2, vo, y+h*s2/2, a, b); s4 = feval (ypfun, t+h u, vo, y+h*s3, a, b);y = y+h* (sl+2*s2+2*s3+s4)/6; i=i+l t=t+h; tout(i) =t; uout(i) =u; yd(i)=Kf*y(l); yout(i, :) = y;% Discrete control

15、 process if abs (round (t/tsamp)-t/tsamp)Ylimitxl=ylimit ; else if ulIlimitx4=Ilimit;else if u4Ilimit x4= Ilimit;else x4 =u4; endendu = Kq*x4;endend % save data to file hout = tout uout yout; save hout.dat hout -ascii; plot (tout, yd, y); gird;% for discrete section% for main loop%存储仿真数据m 函数文件valve

16、.mfunction xd = valve (t, u, x, av, bv) xd = av*x+hv*u;m函数文件hysys. Mfunction yd=hysys(t, u, y, a, b) yd = a* y+h*u上例的仿真结果如图7.3所示图7.3例7.2系统仿真结果7.3基于离散相似法的连续系统仿真前面讨论的连续系统仿真采用的方法是数值积分法。对于连续系统仿真还有一 种方法是离散相似法。所谓离散相似法是首先将连续系统模型离散化，得到等价的 或相似的离散化的模型，然后对相似的离散模型进行仿真计算。根据这一原理，首 先应将连续时间系统模型转换为等价的离散时间系统模型。连续系统离散

17、化处理是 通过采样保持器来实现的。从一个连续时间系统转换为相应的离散时间模型，两者 的等价特性取决于采样保持器。MATLAB提供了将连续系统模型转换为离散时间系统模型的函数 C2D，调用格式为sysd = c2d (sys, Ts)其中，sys为线性连续时间系统；Ts为采样时间；sysd为等价的离散时间系统当采样保持器为零阶保持器时；sysd = c2d (sys, Ts, method)其中，method为离散化方法，MATLAB提供以下几种离散化方法，可以选用:1zoh 为零阶保持器2foh为一阶保持器3讪。n为双线性变换法，齐=需4prewarp为改进的双线性变换法5matched使连续

18、和离散系统具有匹配的 DC增益【例7.3】连续系统传递函数H (s)二十 1 e35s，采样一阶采样保持器，采样s2 +4s+5周期为Ts =0.1s，求其离散化系统模型，并比较离散前后系统阶跃响应用MATLAB编写程序：% MATLAB PROGRAM 7-3sys = tf ( I -1 , 14 5 , td , 0.35 );sysd = c2d ( sysc, 0.l, for)step ( sysc, sysd );运行结果：Tran sfer fun cti on:0.0 1 13 0.0 4 52 -0.0 5 622-0.0 0 9 1 0 4 zW -1.6 2 9八5+0

19、.6 7 0ZMSampling time: 0.1015启 rand (3, 4) ans =0.95280.59820.83680.37590.70410.84070.51870.89860.95390.44280.02220.4290a = zero (2, 2)%zeros函数可以起到声明变量的作用a =0 000 rand (size(a)ans =0.1996 0.5385 0.3031 0.9201%产生1,3区间上的均匀分布随机数(3-1)*rand(2, 3)+lans =2.0506 1.0689 2.53741.6137 2.4307 1.11902 randn其作用是产

20、生标准正态分布的随机数或随机矢量，其调用形式和 rand 类似例 7.5】randn(1,2) ans =%产生 N(1, 0.25)的正态分布-0.4326 -1.6656 randn(1, 2) / 2+1 ans =l.0627 1.14383 randperm其作用是产生一个随机置换，调用形式p = randperm(n) 产生 1: n 的一个随机置换。【例 7.6】p = randperm(4)p =2 3 4 1以 rand 和 randn 为基础可以产生给定分布的随机数， 例如泊松分布和二项分布 M 文件函数 posong 演示了如何产生泊松分布的随机数。function p

21、= posong(a, n, m) w = ones(n, m)*(-a);th = exp(w); r = ones(n, m) dv = r th;p = zeros(n, m);while nnz(dv ) 0r = r.*(rand(n, m).*dv);dv = r th; p = p+dv;end在这里我们对它的 MATLAB 的实现，简要地说明一下。参数 a 指泊松分布的 参数，n, m则指返回的随机矢量的维数。这个程序中有几个地方需要注意：(1)MATLAB 支持矢量的逻辑运算，运算结果是一个元素为 0 或 1 的矢量。 如“ r th”这条语句。(2)函数nnz的作用是统计一

22、个矢量非 0元素的个数。在这里的作用是判断 是不是 r 中所有的元素小于 exp(-a) 。(3)运算符“.*”表示数组元素对元素的乘法，女口 a1, a2,，an*b1, b2,bn = al*bl, a2*b2,an*bn总之这个程序比较充分的利用了 MATLAB 基于向量的特点，产生二项分布随 机数 M 文件函数，请大家自己尝试。其实，像上面的这些工作可以由 MATLAB 完成， MATLAB 的统计工具箱提供 了一个各种分布随机数的发生器函数。安装 MATLAB 的统计工具箱，就可以使用 这些函数。表 7.1 列出了这些函数。你可以通过help命令在MATLAB界面上来获得这些函数的详

23、尽使用方法。 下 面作为示例来看看random的用法。4Random它的作用是产生一个由n ame参数指定的分布的随机数。其调用形式为：R = random (name, A, M, N ),name指定分布的形式，A说明该分布的参数，M , N规定产生的随机数的矩阵维 数， M 为行数， N 为列数，它们的缺省值都是 1(MATLAB 所有函数涉及维数的 参数，顺序都是行在前，列在后)。由于name指定的分布可能不止由一个参数来表 示，所以调用形式要做相应的变化。R = random(name. A, B, M, N)，指定的分布含有两个参数。R = random(name, A, B, C

24、, M, N)，指定的分布含有三个参数。例如，泊松分布只含有一个参数，就要选用第一种形式来调用random ( poiss, 4, 2, 2) ans =6 46 1表7.1 MATLAB统计工具箱里的随机数发生器功能说團betfimdBeta分布binomd:进制分带Expmd捋数分布FmdF分布gaiTun側马井帝gcomd几何分布hygcmd却凡何分布nhinnid负二冷廿布ncfmd非申心F舟布nctmd非申心1号布iwrmrnd泊松分布randjum一揩定的分布mylmd瑞剧廿布tn idT分布unidmdunifmd均匀分布WEibrnd I 唏怕并布7.4.2离散事件系统仿真的一

25、个实例 一一报童问题仿真报童问题是一个古典的概率统计分析问题， 虽然问题本身并不复杂，但作为一个演示实例，可以反映出离散事件系统计算机仿真的很多特征。1 报童问题一报童从报刊发行中心定报后零售， 每卖一份报纸可赚钱a元，若定报后卖不 出去，则可再退回发行处，此时每退一份报要赔钱 b元。虽然每天卖出报的份数是 随机的，但报童可根据以往卖报情况的统计来获得每天卖 k份的概率p(k)，试求报童每天期望受益达到最大的定报量 z。2 数学摸型设报童每天订报z份，而报纸每天卖出y份，我们假设y的分布为p(y =k) = Pk k =0,1,2,考虑到报童每天的损失有如下两种情形。(1)供过于求。因退货造成

26、的平均损失为：QOG (z- k)pkk=0(2)供不应求。因缺货造成的平均损失为QOC2 二 a (k -z)pkk=z 1所以，每天的期望损失费(也可以从总收益的角度来考虑)为C(z) = Ci C2现在我们的目标是求出使得每天期望损失最小的定报量， 换言之，就是使报童的每天期望总收益达到最大。写成一个目标函数的形式z =Cm?n(z)约束条件如z的取值范围，要受到报童的资本多少的影响。只有在特殊的概率分布情况下，我们才可以推导出C(z)的解析形式，并通过求 极值的方法来求解。但在实际的应用中，这样的思路往往是行不通的。但是可以通 过枚举所有可能的订报量，求出对应的平均损失，进行比较求出满

27、足条件的 z,这里搜索域通常是有限的。以上就是一个比较简单的计算机仿真方法。3.报童问题的计算机仿真对于给定每一订报量Z值，利用离散随机变量采样算法产生给定分布的随机数 R,用来表示报童当天卖出的报纸数，从而可以计算出一天的损失以及一个阶段的 平均损失。这里比较关键的一点是如何产生服从给定分布的随机变量， 这个内容在本教材中不再详尽介绍。而且在实际的应用中，分布并非总是给定了的，需要我们收集数据，并从中辨识分布，进行参数估计。图 7.5是根据上述思路设计出的仿真流程框图。其中各变量含义如下:S 最小损失值这里，a和b是这个问题的两个参数4计算机仿真根据图 7.5 所示框图，我们不难写出报童问题

28、的仿真程序，其 MATLAB 源码 如下。function superz, supers = baotong ( tm, g, a, b)z = 1;supers =1000;while zr;s = sum(z-r)*b).*dv);s = s+sum(r-z)*a).*(l-dv); aver_s =s/tm;if supers=aver_ssupers = aver_s;superz = z;end;z = z +1;end;上面的代码，是在随机变量均匀分布的假设下编写的，对于其他的分布，大家 只要用相应的随机数发生器进行替换即可5仿真结果和输出数据分析在分布为均匀分布，参数 a=0.2， b=0.4 的条件下，在 MATLAB 命令窗口运 行仿真程序就可以得出