1、毕业设计 系统仿真的 MATLAB 实现（仅供参考）由于计算机技术的高速发展,我们可以借助计算机完成系统的数字仿真.综前所述,数字仿真实质上是根据被研究的真实系统的模型,利用计算机进行实验研究的一种方法.仿真的主要过程是:建立模型、仿真运行和分析研究仿真结果 .仿真运行就是借助一定的算法,获得系统的有关信息.MATLAB 是一种面向科学与工程计算的高级语言,它集科学计算、自动控制、信号处理、神经网络和图像处理等学科的处理功能于一体,具有极高的编程效率.MATLAB 是一个高度集成的系统,MATLAB 提供的 Simulink 是一个用来对动态系统进行建模、仿真和分析的软件包,它支持线性和非线性
2、系统,能够在连续时间域、离散时间域或者两者的混合时间域里进行建模,它同样支持具有多种采样速率的系统.在过去几年里,Simulink 已经成为数学和工业应用中对动态系统进行建模时使用得最为广泛的软件包.MATLAB 仿真有两种途径:(1)MATLAB 可以在 SIMULINK 窗口上进行面向系统结构方框图的系统仿真;(2) 用户可以在MATLAB 的 COMMAND 窗口下,用运行 m 文件,调用指令和各种用于系统仿真的函数,进行系统仿真.这两种方式可解决任意复杂系统的动态仿真问题,前者编辑灵活,而后者直观性强,实现可视化编辑.下面介绍在 MATLAB 上实现几类基本仿真.7.1 计算机仿真的步
3、骤在学习计算机仿真以前,让我们先总结一下计算机仿真的步骤.计算机仿真,概括地说是一个“建模实验 分析“ 的过程,即仿真不单纯是对模型的实验 ,还包括从建模到实验再到分析的全过程.因此进行一次完整的计算机仿真应包括以下步骤:(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文件,还有一个函数是 ode1.C,是直接用 C 语言编写的.函数 ode23( )是用 Runge-Kutta 法求解微分方程.它是一种采用三阶积分算法、二阶误差估计、变积分步长的低阶积分算法,调用格式为T, Y = ode23
6、 ( F, TSPAN, YO, OPTIONS )其中,F 为系统模型文件名,模型为 y = f( t, y )形式;TSPAN = To TFINAL 为积分计算时间,初值为 To,终值为 TFINAL;YO 为系统输出初始值;OPTIONS 选项积分计算相对允差 RelTol 和绝对允差 AbsTol,当缺省时,Reltol=1e-3, AbsTol=1e-6T 为计算点时间向量,Y 为微分方程的解.函数 ode45( )也是用 Runge-Kutta 法求解微分方程,它是变步长的一种中等阶次积分算法,调用格式为T, Y = ode45 ( F , TSPAN, YO, OPTIONS
7、)各项含义同上.函数 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 中,系统仿真有变步距和固定步距两种工作方式 .在变步距仿真中,可选用 od
8、e45、ode23、ode113、ose15s和 ode23s.在固定步距仿真中,可选用材 ode5、ode4 、ode3、ode2 和 ode1.ode5 是 5 阶 Runge-Kutta 法,ode4 是 4 阶 Runge-Kutta 法,odel 是 Euler 法.【例 7.1】求微分方程.用 MATLAB 编写程序:建立一个 m 函数文件 dfun.mfuncion y=dhfn ( t,x )y=sqrt ( 3*x )+5 ;在 MATLAB COMMAND 窗口下运行% MATLAB PROGRAM 7-1 t, x =ode23 ( dfun,0 l0,1) ;Plot
9、( t,x ) ;7.2.2 基于数值积分法的连续系统的数字仿真对于一个 n 阶微分方程表示的连续系统,也可以用具有 n 个状态变量的状态空间表达式来描述:状态方程实际上是 n 个一阶微分方程组成的方程组.如果应用四阶 Runge-Kutta 法进行仿真计算.则该式可以改写为式中,K1,K2,K3,K4 均为 n 维列向量.系统输出为yk=CXk+DU( tk )类似地,可将各种数值积分方法用于连续系统的仿真中.【例 7.2】用 MATLAB 编写图 7.2 所示液压控制系统的仿真程序.图 7.2 系统方框图用 MATLAB 编写仿真程序,采用四阶 Runge-Kutta 法:% MATLAB
10、 PROGRAM 7-2Input system data 调入数据文件hynat;% Input system function;调入系统模型ypfun1 = valve ;ypfun = hysys ;% Initialization %初始化yref = 5000;Referent value of system outputx0 = 0 0Initial value of servo valvey0 = 0 0 0Initial value of hydraulic cylinderu0 = 0;t0 = 0;Start time of simulationtfinal = 1;End
11、 time of simulationtsamp = 0.01;Sample periodh = 0.001;Simulation step sizeY1imit = 512;Ilimit = 40;max_epoch = fix ( tfinal / h )-1;t = t0;u1 = u0;x = x0;y = y0;tout = zeros ( max_epoch, 1 );uout = zeros ( max_epoch, 1);yd = zeros ( max_epoch, 1 );yout = zeros ( max_epoch, legth ( y ) );i = 1;tout
12、( i ) = t;uout ( i ) = 1;yd ( i ) = Kf * y ( l );yout ( i, : ) = y;% The main loopfor 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 ); %模块 valve 仿真sv3 = feval (ypfun1, t+h/2, u, x+h*sv2/2, av, hv );sv4 = feval (ypfu
13、n1, t+h, u, x+h*sv3, av, hv);x = x+h*(svl+2*sv2+2*sv3+sv4)/6;vo = cv*x; % Compute output of cylinderxl = feval (ypfun, t, vo, y, a, b);s2 = feval (ypfun, t+h/2, vo, y+h*s1/2, a, h)模块 hysya 仿真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
14、*s3+s4)/6;i=i+lt=t+h;tout(i) =t;uout(i) =u;yd(i)=Kf*y(l);yout(i, :) = y;% Discrete control process %离散控制量计算if abs (round (t/tsamp)-t/tsamp)Ylimitxl=ylimit ;else if ulIlimitx4=Ilimit;else if u4 rand (3, 4)ans =0.9528 0.5982 0.8368 0.37590.7041 0.8407 0.5187 0.89860.9539 0.4428 0.0222 0.4290a = zero (
15、2, 2)zeros 函数可以起到声明变量的作用a =0 00rand (size(a)ans =0.1996 0.53850.3031 0.9201(3-1)*rand(2, 3)+l %产生1,3区间上的均匀分布随机数ans =2.0506 1.0689 2.53741.6137 2.4307 1.11902.randn其作用是产生标准正态分布的随机数或随机矢量,其调用形式和 rand 类似.【例 7.5】randn(1,2)ans =-0.4326 -1.6656randn(1, 2) / 2+1 %产生 N(1, 0.25)的正态分布.ans =l.0627 1.14383.randp
16、erm其作用是产生一个随机置换,调用形式p = randperm(n) 产生 1:n 的一个随机置换.【例 7.6】p = randperm(4)p =2 3 4 1以 rand 和 randn 为基础可以产生给定分布的随机数,例如泊松分布和二项分布.M 文件函数 posong 演示了如何产生泊松分布的随机数.function p = 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 =
17、 r th;p = p+dv;end在这里我们对它的 MATLAB 的实现,简要地说明一下.参数 a 指泊松分布的参数,n, m 则指返回的随机矢量的维数.这个程序中有几个地方需要注意:(1)MATLAB 支持矢量的逻辑运算,运算结果是一个元素为 0 或 1 的矢量.如 “r th“这条语句.(2)函数 nnz 的作用是统计一个矢量非 0 元素的个数.在这里的作用是判断是不是 r 中所有的元素小于 exp(-a).(3)运算符“.*“ 表示数组元素对元素的乘法 ,如a1, a2, ,an.*b1, b2, , bn = al*bl, a2*b2, , an*bn总之这个程序比较充分的利用了 M
18、ATLAB 基于向量的特点,产生二项分布随机数 M 文件函数,请大家自己尝试.其实,像上面的这些工作可以由 MATLAB 完成,MATLAB 的统计工具箱提供了一个各种分布随机数的发生器函数.安装MATLAB 的统计工具箱,就可以使用这些函数.表 7.1 列出了这些函数.你可以通过 help 命令在 MATLAB 界面上来获得这些函数的详尽使用方法.下面作为示例来看看 random 的用法.4.Random它的作用是产生一个由 name 参数指定的分布的随机数.其调用形式为:R = random (name, A, M, N ), name 指定分布的形式,A 说明该分布的参数,M,N 规定产
19、生的随机数的矩阵维数,M 为行数,N 为列数,它们的缺省值都是1(MATLAB 所有函数涉及维数的参数 ,顺序都是行在前,列在后).由于 name 指定的分布可能不止由一个参数来表示,所以调用形式要做相应的变化.R = random(name, A, B, M, N),指定的分布含有两个参数 .R = random(name, A, B, C, M, N),指定的分布含有三个参数 .例如,泊松分布只含有一个参数,就要选用第一种形式来调用.random ( poiss, 4, 2, 2)ans =41表 7.1 MATLAB 统计工具箱里的随机数发生器7.4.2 离散事件系统仿真的一个实例报童问
20、题仿真报童问题是一个古典的概率统计分析问题,虽然问题本身并不复杂,但作为一个演示实例,可以反映出离散事件系统计算机仿真的很多特征.1.报童问题一报童从报刊发行中心定报后零售,每卖一份报纸可赚钱 a 元,若定报后卖不出去,则可再退回发行处,此时每退一份报要赔钱 b 元.虽然每天卖出报的份数是随机的,但报童可根据以往卖报情况的统计来获得每天卖 k 份的概率 p(k),试求报童每天期望受益达到最大的定报量 z.2.数学摸型设报童每天订报 z 份 ,而报纸每天卖出 y 份,我们假设 y 的分布为考虑到报童每天的损失有如下两种情形.(1)供过于求.因退货造成的平均损失为:(2)供不应求.因缺货造成的平均
21、损失为所以,每天的期望损失费(也可以从总收益的角度来考虑)为现在我们的目标是求出使得每天期望损失最小的定报量,换言之,就是使报童的每天期望总收益达到最大.写成一个目标函数的形式约束条件如 z 的取值范围 ,要受到报童的资本多少的影响.只有在特殊的概率分布情况下,我们才可以推导出 C(z)的解析形式,并通过求极值的方法来求解.但在实际的应用中,这样的思路往往是行不通的.但是可以通过枚举所有可能的订报量,求出对应的平均损失,进行比较求出满足条件的 z,这里搜索域通常是有限的.以上就是一个比较简单的计算机仿真方法.3.报童问题的计算机仿真对于给定每一订报量 Z 值,利用离散随机变量采样算法产生给定分
22、布的随机数 R,用来表示报童当天卖出的报纸数,从而可以计算出一天的损失以及一个阶段的平均损失.这里比较关键的一点是如何产生服从给定分布的随机变量,这个内容在本教材中不再详尽介绍.而且在实际的应用中,分布并非总是给定了的,需要我们收集数据,并从中辨识分布,进行参数估计.图7.5 是根据上述思路设计出的仿真流程框图.其中各变量含义如下:Tm一轮试验的预定模拟天数 T一轮实验的仿真天数累计值Z订报量 Z最优订报量G定报量 Z 之上界 S1损失值之累计值S最小损失值这里,a 和 b 是这个问题的两个参数 .图 7.5 报童问题计算机仿真流程框图4.计算机仿真根据图 7.5 所示框图,我们不难写出报童问
23、题的仿真程序,其 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
24、,b=0.4 的条件下,在 MATLAB 命令窗口运行仿真程序就可以得出仿真的结果.例如 z, s = baotong (5, 10, 0.2, 0.4)z =4s = 0.2400此外,我们还可以改变参数 a、b 的值,观察相应最优值的变化,用 MATLAB 可以很方便的画出 C-a,C-b 的变化曲线.6.报童问题模拟系统的推广和应用报童仿真问题的改进推广和应用可以有以下几个方面:(1)求每天的卖出报数服从任意分布的情况下,使报童收益最大意义下的最优订报量 Z;(2)对报纸总发行量进行测算;(3)适当修改仿真系统,可将其用于其它系统,例如企业的订货和库存策略研究.7.5 SIMULINK
25、动态仿真SIMULINK 是 MATLAB 重要软件包,用于对动态系统仿真,它适用于连续系统和离散系统,也适用线性系统和非线性系统.它采用系统模块直观地描述系统典型环节,因此可十分方便地建立系统模型而不需要花较多时间编程.正由于这些特点,SIMULINK 广泛流行,被认为是最受欢迎的仿真软件 .本节对 SIMULINK 和它的 MATLAB 实现,作一个简单介绍,为大家今后运用 SIMULINK 打下基础.SIMULINK 实际上是面向结构的系统仿真软件 .利用它进行系统仿真的步骤 如下:启动 SIMULINK,进入 SIMULINK 窗口;在 SIMULINK 窗口下,借助 SIMULINK
26、 模块库,创建系统框图模型并调整模块参数 ;设置仿真参数后,启动仿真;输出仿真结果.7.5.1 SIMULINK 窗口及模型库用户首先进入 MATLAB COMMAND 窗口,键入 Simulink,立即弹出 Simulink 模块库窗口,如图 7.6 所示.Library:Simulink 窗口是 Simulink 模块库窗口 ,共有 7 个标准模块子库.一个 MATLAB 演示程序图标.这七个标准模块子库是 Source(信号源)、Sink(显示输出)、Discrete(离散)、Linear(线性)、Nonlinear( 非线性)、Connections(连接)、 Blocksets &
27、Toolboxes.每个模块子库又包含许多常用的模块 ,供用户建模时使用,SIMULINK 提供的库模块清单见附录.图 7.6 Simulink 模块库窗口7.5.2 系统框图模型的建立系统框图模型建立的过程如下:1.建立模型窗口建立新的模型窗口常有四种方法:(1)在 MATLAB COMMAND 窗口下,键入 Simulink,弹出 Simulink 模块库窗口同时,也弹出一个 Untitled 窗口,该窗口为未取名的模型窗口.用户可在该窗口下建立新的系统框图模型.(2)在 Simulink 窗口下,用鼠标选取菜单 File 中 New 子菜单的 Mode 后,会弹出一个 untitled
28、窗口,如图 7.7 所示.该窗口供用户建立系统框图模型.图 7.7 模型窗口(3)若模型文件已存在,Simulink 窗口下,选择菜单 File 中 Open 命令,输入文件名,即打开一个已存在的方框图模型.(4)若模型文件已存在,在 MATLAB COMMAND 窗口下,选择 File中Run Script命令,输入文件名,也可打开一个已存在的方框图模型.(5)若模型文件已存在,在 MATLAB COMMAND 窗口下,键入文件名,也可打开一个已存在的模型.2.选取模块或模块组在建立框图模型过程中,需进行如拷贝、删除模块等操作.必须首先选择模块或模块组,具体操作如下:(1)在模型或模块库的窗
29、口内,找出所需模块图标,用鼠标左键单击.图标四角出现黑圆点,表示该模块已被选中.(2)在模型或模块库窗口内,用鼠标左键在窗口矩形边界两个对角单击一下,即生成一个边界框将所需几个框块图标包围,松开鼠标,则边界框内模型和连接线出现黑圆点,表示这些模型(包括在连接线)均被选中.用同样方法可以选取一个系统框图模型的全部模块.3.模块拷贝及删除用户在建立自己模型的时候,常常需要从 Simulink 模块库、其他模块库或其他模型窗口中复制所需的模块并移动至自己的模型窗口内.有两种操作方法:鼠标拖动方法和菜单命令法 .鼠标拖动法如下:打开模块库窗口或模型窗口.将鼠标移至要拷贝的模块图标上,按下鼠标左键并保持
30、.移动鼠标将模块图标拖至目标模型窗口一定位置.松开鼠标左键,模块图标保留在目标模型窗口内,模块拷贝完成.菜单命令复制方法如下:打开模块库窗口或模型窗口,按上述方法选取要复制的模块.从Edit菜单中选取Copy 命令.用鼠标单击目标模型适当位置.(4)从Edit菜单中选取Paste命令,要拷贝的模块出现在目标模型窗口内,模块复制完成.模块的删除有两种方法:(1)选取要删除的模型,从Edit 菜单中选取Clear 或Cut命令,用Cut 命令删除的模块允许使用Paste粘贴在另一个地方.(2)选取要删除的模块,并按Del键.SIMULINK 允许模块更名,图标大小改变、模块图标移动、模块图标旋转等
31、操作.4.模块参数设置用鼠标双击待设置参数的模块图标.打开模块对话框,按对话框栏目中提供的信息,输入或改变模块参数.按Close模块参数设置或修改完成.各模块对话框需设置参数不同.大家打开对话框后,一目了然.5.模块连接线模块之间的连接线是信号线,每根连接线都表示标量或向量信号的传输,连接线的箭头表示信号流向.连接线把一个模块的输出端口和另一个模块的输入端口连接起来,也可以利用分支线把一个模块的输出端口和几个模块的输入端口连接 起来.模块间连接线生成方法如下:将鼠标置于某模块的输出端口上,立即呈现出一个十字形光标.拖动十字光标至另一模块的输入端口,在两模块之间生成一条带箭头的连线,模块连线完成
32、.分支线的生成方法如下:用鼠标键单击被选取好信号线,按鼠标右键对准该信号线后显示十字形光标,拖动它会生成信号线的一条分支线.按此法,可生成多条分支线.6.模型文件取名及保存一旦把模型窗口上各模型连接起来,一个系统方框图模型建立工作就已完成.选择模型窗口File菜单中Save as命令,弹出对话框,填入模型文件名.如 hexam,按保存(S),模型以文件名 hexam.mdl 保存在用户的子目录内.按上面方法建立一个简单系统框图模型,如图 7.8 所示,文件名为 hexam.图 7.8 简单系统框图模型7.5.3 系统仿真运行系统仿真运行常有两种方法进行,一种是利用 SIMULINK 菜单命令,
33、二是在 MABLAB COMMAND 窗口下输入有关命令.1.SIMULINK 模型窗口下的仿真运行在 SIMULINK 模型窗口下进行仿真操作简单、直观 ,不必记忆命令的语法规则,人机交互方式选择或修改仿真参数,模型参数等.具体操作如下:打开系统模型窗口.(2)从菜单Simulation中选取Parameters,弹出 Simulation Parameters 对话框,如图 7.9 所示.图 7.9 仿真参数对话框Simulation Parameters 对话框用三页管理仿真参数 . Solver 页,设置或修改参数有:仿真时间:开始时间,终止时间仿真方法类型:变步距(积分算法一):discrete ode45, ode23, odel13, odel5s, ode23s定步距(积分算法二):ode5,ode4, ode3, ode2, ode1, discrete
