matlab系统仿真.doc

第七章 系统仿真的MATLAB实现由于计算机技术的高速发展，我们可以借助计算机完成系统的数字仿真。综前所述，数字仿真实质上是根据被研究的真实系统的模型，利用计算机进行实验研究的一种方法。仿真的主要过程是：建立模型、仿真运行和分析研究仿真结果。仿真运行就是借助一定的算法，获得系统的有关信息。MATLAB是一种面向科学与工程计算的高级语言，它集科学计算、自动控制、信号处理、神经网络和图像处理等学科的处理功能于一体，具有极高的编程效率。MATLAB是一个高度集成的系统，MATLAB提供的Simulink是一个用来对动态系统进行建模、仿真和分析的软件包，它支持线性和非线性系统，能够在连续时间域、离散时间域或者两者的混合时间域里进行建模，它同样支持具有多种采样速率的系统。在过去几年里，Simulink已经成为数学和工业应用中对动态系统进行建模时使用得最为广泛的软件包。 MATLAB仿真有两种途径：（1）MATLAB可以在SIMULINK窗口上进行面向系统结构方框图的系统仿真；（2）用户可以在MATLAB的COMMAND窗口下，用运行m文件，调用指令和各种用于系统仿真的函数，进行系统仿真。这两种方式可解决任意复杂系统的动态仿真问题，前者编辑灵活，而后者直观性强，实现可视化编辑。 下面介绍在MATLAB上实现几类基本仿真。7.1 计算机仿真的步骤在学习计算机仿真以前，让我们先总结一下计算机仿真的步骤。计算机仿真，概括地说是一个“建模实验分析”的过程，即仿真不单纯是对模型的实验，还包括从建模到实验再到分析的全过程。因此进行一次完整的计算机仿真应包括以下步骤：（1）列举并列项目每一项研究都应从说明问题开始，问题由决策者提供或由熟悉问题的分析者 提供。（2）设置目标及完整的项目计划目标表示仿真要回答的问题、系统方案的说明。项目计划包括人数、研究费用以及每一阶段工作所需时间。（3）建立模型和收集数据模型和实际系统没有必要一一对应，模型只需描述实际系统的本质或者描述系统中所研究部分的本质。因此，最好从简单的模型开始，然后进一步建立更复杂的模型。 （4）编制程序和验证利用数学公式、逻辑公式和算法等来表示实际系统的内部状态和输入/输出的关系。建模者必须决定是采用通用语言如MATLAB、FORTRAN、C还是专用仿真语言来编制程序。在本教材中，我们选择的是MATLAB和其动态仿真工具Simulink。（5）确认确认指确定模型是否精确地代表实际系统。它不是一次完成，而是比较模型和实际系统特性的差异，不断对模型进行校正的迭代过程。（6）实验设计 确定仿真的方案、初始化周期的长度、仿真运行的长度以及每次运行的重复 次数。 （7）生产性运行和分析 通常用于估计被仿真系统设计的性能量度。利用理论定性分析、经验定性分析或系统历史数据定量分析来检验模型的正确性，利用灵敏度分析等手段来检验模型的稳定性。 （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 ( F, TSPAN, YO, OPTIONS )其中，F为系统模型文件名，模型为y = f( t, y )形式；TSPAN = To TFINAL 为积分计算时间，初值为To，终值为TFINAL；YO为系统输出初始值；OPTIONS选项积分计算相对允差 RelTol 和绝对允差 AbsTol，当缺省时， Reltol1e-3, AbsTol1e-6T为计算点时间向量，Y为微分方程的解。函数ode45( )也是用Runge-Kutta法求解微分方程，它是变步长的一种中等阶次积分算法，调用格式为T, Y = ode45 ( 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中，系统仿真有变步距和固定步距两种工作方式。在变步距仿真中，可选用ode45、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 ydhfn ( t，x ) ysqrt ( 3*x )5 ; 在MATLAB COMMAND窗口下运行 % MATLAB PROGRAM 7-1 t, x ode23 ( dfun，0 l0，1) ; Plot ( t，x ) ;7.2.2 基于数值积分法的连续系统的数字仿真对于一个n阶微分方程表示的连续系统，也可以用具有n个状态变量的状态空间表达式来描述：状态方程实际上是n个一阶微分方程组成的方程组。如果应用四阶Runge-Kutta法进行仿真计算。则该式可以改写为式中，K1，K2，K3，K4均为n维列向量。系统输出为ykCXkDU( tk )类似地，可将各种数值积分方法用于连续系统的仿真中。【例7.2】用MATLAB编写图7.2所示液压控制系统的仿真程序。图7.2 系统方框图用MATLAB编写仿真程序，采用四阶Runge-Kutta法：% MATLAB PROGRAM 7-2% * Input system data * %调入数据文件hynat;% Input system function; %调入系统模型ypfun1 = valve ;ypfun = hysys ;% Initialization %初始化yref = 5000; % Referent value of system outputx0 = 0 0; % Initial value of servo valvey0 = 0 0 0; % Initial value of hydraulic cylinderu0 = 0; t0 = 0; % Start time of simulation tfinal = 1; % End 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 ( i ) = t;uout ( i ) = 1;yd ( i ) = Kf * y ( l );yout ( i, : ) = y;% The main loop 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 ); %模块valve仿真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 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*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.3759 0.7041 0.8407 0.5187 0.8986 0.9539 0.4428 0.0222 0.4290a = zero (2, 2) %zeros函数可以起到声明变量的作用a = 0 00 0rand (size(a)ans = 0.1996 0.5385 0.3031 0.9201(3-1)*rand(2, 3)+l %产生1,3区间上的均匀分布随机数ans = 2.0506 1.0689 2.5374 1.6137 2.4307 1.11902randn其作用是产生标准正态分布的随机数或随机矢量，其调用形式和rand类似。【例7.5】randn(1,2)ans =-0.4326 -1.6656randn(1, 2) / 2+1 %产生N(1, 0.25)的正态分布。ans = l.0627 1.14383randperm其作用是产生一个随机置换，调用形式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 = 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总之这个程序比较充分的利用了MATLAB基于向量的特点，产生二项分布随机数M文件函数，请大家自己尝试。其实，像上面的这些工作可以由MATLAB完成，MATLAB的统计工具箱提供了一个各种分布随机数的发生器函数。安装MATLAB的统计工具箱，就可以使用这些函数。表7.1列出了这些函数。你可以通过help命令在MATLAB界面上来获得这些函数的详尽使用方法。下面作为示例来看看random的用法。4Random它的作用是产生一个由name参数指定的分布的随机数。其调用形式为：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, M, N)，指定的分布含有三个参数。例如，泊松分布只含有一个参数，就要选用第一种形式来调用。random ( poiss, 4, 2, 2)ans =6 46 1表7.1 MATLAB统计工具箱里的随机数发生器7.4.2 离散事件系统仿真的一个实例报童问题仿真报童问题是一个古典的概率统计分析问题，虽然问题本身并不复杂，但作为一个演示实例，可以反映出离散事件系统计算机仿真的很多特征。 1报童问题 一报童从报刊发行中心定报后零售，每卖一份报纸可赚钱a元，若定报后卖不出去，则可再退回发行处，此时每退一份报要赔钱b元。虽然每天卖出报的份数是随机的，但报童可根据以往卖报情况的统计来获得每天卖k份的概率p(k)，试求报童每天期望受益达到最大的定报量z。 2数学摸型 设报童每天订报z份，而报纸每天卖出y份，我们假设y的分布为考虑到报童每天的损失有如下两种情形。（1）供过于求。因退货造成的平均损失为：（2）供不应求。因缺货造成的平均损失为所以，每天的期望损失费（也可以从总收益的角度来考虑）为现在我们的目标是求出使得每天期望损失最小的定报量，换言之，就是使报童的每天期望总收益达到最大。写成一个目标函数的形式约束条件如z的取值范围，要受到报童的资本多少的影响。只有在特殊的概率分布情况下，我们才可以推导出C(z)的解析形式，并通过求极值的方法来求解。但在实际的应用中，这样的思路往往是行不通的。但是可以通过枚举所有可能的订报量，求出对应的平均损失，进行比较求出满足条件的z，这里搜索域通常是有限的。以上就是一个比较简单的计算机仿真方法。3报童问题的计算机仿真对于给定每一订报量Z值，利用离散随机变量采样算法产生给定分布的随机数R，用来表示报童当天卖出的报纸数，从而可以计算出一天的损失以及一个阶段的平均损失。这里比较关键的一点是如何产生服从给定分布的随机变量，这个内容在本教材中不再详尽介绍。而且在实际的应用中，分布并非总是给定了的，需要我们收集数据，并从中辨识分布，进行参数估计。图7.5是根据上述思路设计出的仿真流程框图。其中各变量含义如下：Tm一轮试验的预定模拟天数 T一轮实验的仿真天数累计值Z订报量 Z最优订报量G定报量Z之上界 S1损失值之累计值S最小损失值这里，a和b是这个问题的两个参数。图7.5 报童问题计算机仿真流程框图4计算机仿真根据图7.5所示框图，我们不难写出报童问题的仿真程序，其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命令窗口运行仿真程序就可以得出仿真的结果。例如 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动态仿真 SIMULINK是MATLAB重要软件包，用于对动态系统仿真，它适用于连续系统和离散系统，也适用线性系统和非线性系统。它采用系统模块直观地描述系统典型环节，因此可十分方便地建立系统模型而不需要花较多时间编程。正由于这些特点，SIMULINK广泛流行，被认为是最受欢迎的仿真软件。本节对SIMULINK和它的MATLAB实现，作一个简单介绍，为大家今后运用SIMULINK打下基础。SIMULINK实际上是面向结构的系统仿真软件。利用它进行系统仿真的步骤 如下：（1） 启动SIMULINK，进入SIMULINK窗口；（2） 在SIMULINK窗口下，借助SIMULINK模块库，创建系统框图模型并调整模块参数；（3） 设置仿真参数后，启动仿真；（4） 输出仿真结果。7.5.1 SIMULINK窗口及模型库用户首先进入MATLAB COMMAND窗口，键入Simulink，立即弹出Simulink模块库窗口，如图7.6所示。Library：Simulink窗口是Simulink模块库窗口，共有7个标准模块子库。一个MATLAB演示程序图标。这七个标准模块子库是Source（信号源）、Sink（显示输出）、Discrete（离散）、Linear（线性）、Nonlinear（非线性）、Connections（连接）、Blocksets & Toolboxes。每个模块子库又包含许多常用的模块，供用户建模时使用，SIMULINK提供的库模块清单见附录。图7.6 Simulink模块库窗口7.5.2 系统框图模型的建立 系统框图模型建立的过程如下： 1建立模型窗口 建立新的模型窗口常有四种方法： （1）在MATLAB COMMAND窗口下，键入Simulink，弹出Simulink模块库窗口同时，也弹出一个Untitled窗口，该窗口为未取名的模型窗口。用户可在该窗口下建立新的系统框图模型。（2）在Simulink窗口下，用鼠标选取菜单 File 中 New 子菜单的 Mode 后，会弹出一个untitled窗口，如图7.7所示。该窗口供用户建立系统框图模型。图7.7 模型窗口 （3）若模型文件已存在，Simulink窗口下，选择菜单 File 中 Open 命令，输入文件名，即打开一个已存在的方框图模型。（4）若模型文件已存在，在MATLAB COMMAND窗口下，选择 File中Run Script命令，输入文件名，也可打开一个已存在的方框图模型。（5）若模型文件已存在，在MATLAB COMMAND窗口下，键入文件名，也可打开一个已存在的模型。2选取模块或模块组 在建立框图模型过程中，需进行如拷贝、删除模块等操作。必须首先选择模块或模块组，具体操作如下： （1）在模型或模块库的窗口内，找出所需模块图标，用鼠标左键单击。图标四角出现黑圆点，表示该模块已被选中。 （2）在模型或模块库窗口内，用鼠标左键在窗口矩形边界两个对角单击一下，即生成一个边界框将所需几个框块图标包围，松开鼠标，则边界框内模型和连接线出现黑圆点，表示这些模型（包括在连接线）均被选中。 用同样方法可以选取一个系统框图模型的全部模块。 3模块拷贝及删除用户在建立自己模型的时候，常常需要从Simulink模块库、其他模块库或其他模型窗口中复制所需的模块并移动至自己的模型窗口内。有两种操作方法：鼠标拖动方法和菜单命令法。 鼠标拖动法如下：（1） 打开模块库窗口或模型窗口。（2） 将鼠标移至要拷贝的模块图标上，按下鼠标左键并保持。（3） 移动鼠标将模块图标拖至目标模型窗口一定位置。（4） 松开鼠标左键，模块图标保留在目标模型窗口内，模块拷贝完成。 菜单命令复制方法如下：（1） 打开模块库窗口或模型窗口，按上述方法选取要复制的模块。（2） 从Edit菜单中选取Copy命令。（3） 用鼠标单击目标模型适当位置。（4）从Edit菜单中选取Paste命令，要拷贝的模块出现在目标模型窗口内，模块复制完成。 模块的删除有两种方法：（1）选取要删除的模型，从Edit菜单中选取Clear或Cut命令，用Cut命令删除的模块允许使用Paste粘贴在另一个地方。（2）选取要删除的模块，并按Del键。 SIMULINK允许模块更名，图标大小改变、模块图标移动、模块图标旋转等操作。 4模块参数设置 用鼠标双击待设置参数的模块图标。打开模块对话框，按对话框栏目中提供的信息，输入或改变模块参数。 按Close模块参数设置或修改完成。 各模块对话框需设置参数不同。大家打开对话框后，一目了然。 5模块连接线 模块之间的连接线是信号线，每根连接线都表示标量或向量信号的传输，连接线的箭头表示信号流向。连接线把一个模块的输出端口和另一个模块的输入端口连接起来，也可以利用分支线把一个模块的输出端口和几个模块的输入端口连接 起来。 模块间连接线生成方法如下：（1） 将鼠标置于某模块的输出端口上，立即呈现出一个十字形光标。（2） 拖动十字光标至另一模块的输入端口，在两模块之间生成一条带箭头的连线，模块连线完成。 分支线的生成方法如下： 用鼠标键单击被选取好信号线，按鼠标右键对准该信号线后显示十字形光标，拖动它会生成信号线的一条分支线。按此法，可生成多条分支线。 6模型文件取名及保存 一旦把模型窗口上各模型连接起来，一个系统方框图模型建立工作就已完成。 选择模型窗口File菜单中Save as命令，弹出对话框，填入模型文件名。如hexam，按保存(S)，模型以文件名hexam.mdl保存在用户的子目录内。按上面方法建立一个简单系统框图模型，如图7.8所示，文件名为hexam。图7.8 简单系统框图模型7.5.3 系统仿真运行 系统仿真运行常有两种方法进行，一种是利用SIMULINK菜单命令，二是在MABLAB COMMAND窗口下输入有关命令。1SIMULINK模型窗口下的仿真运行 在SIMULINK模型窗口下进行仿真操作简单、直观，不必记忆命令的语法规则，人机交互方式选择或修改仿真参数，模型参数等。具体操作如下：（1） 打开系统模型窗口。（2）从菜单Simulation中选取Parameters，弹出Simulation Parameters对话框，如图7.9所示。图7.9 仿真参数对话框Simulation Parameters对话框用三页管理仿真参数。 Solver页，设置或修改参数有：仿真时间：开始时间，终止时间仿真方法类型： 变步距（积分算法一）：discrete ode45, ode23, odel13, odel5s, ode23s 定步距（积分算法二）：ode5，ode4, ode3, ode2, ode1, discrete 初始步距 最大步距 Workspace I/O页设置或修改参数有： 从工作空间调用输入 把数据存入工作空间（时间、状态、输出） 状态值（初始值调入，初始值存放） 存储选择 Diagnostics页设置或选择项目：仿真过程警告信息（3）应用仿真参数在通过Simulation Parameters对话框设置或修改好仿真参数后，准备应用到用户模型仿真中，按对话中Apply键，再按对话框中Close键，关闭对话框。（4）仿真开始和停止从Simulation菜单中选取Start命令。仿真立即开始运行，这时Start变为Stop。若要停止仿真，从Simulation菜单中选择Stop。仿真运行立即停止。若要使仿真运行暂停，可从Simulation菜单中选择Pause，这时Pause变为Continue。若要使仿真继续运行，选择Continue。2MATLAB的COMMAND窗口下的仿真运行 在MATLAB COMMAND窗口下运行SIMULINK模型仿真也可进行，方法是在COMMAND窗口下调用函数SIM。调用格式为t, x, y = sim ( model, timespan, option, ut, p1, pn)其中，t为返回的仿真时间向量；y为返回的输入矩阵，从输出模块输出的数据；model为框图模型文件名或MEX-文件或M-文件S函数；timespan = tStart，tFinal；Options为仿真参数选项，由函数Simset( )产生；ut为加在第一级输入模块的外部输入，如一个函数u(t)等；p1, pn为S-函数参数，如Sys = Sfun(t，x，u，flay，p1,p2, pn)仅用于M-文件和MEX-文件系统。例如，图7.8所示系统，在MATLAB COMMAND窗口下仿真t, y = Sim ( hexm );plot ( t, y )比较在SIMILINK模型窗口小的仿真结果，可见两种方法仿真结果相同。由于在MATLAB的COMMAND窗口下运行SIMULINK模型不如直接在SIMULINK模型窗口下仿真方便、直观，因此不推荐使用。利用Mum模块和Scope模块可显示多个变量的运动轨迹，如图7.10所示。图7.10 显示多个变量轨迹7.5.4 仿真结果输出SIMULINK可通过下列三种方法之一绘制输出轨迹曲线。1利用Scope模块在仿真期间，Scope模块动态地显示输出轨迹的一个简单例子如图7.10所示。Scope模块还具有放大、缩小功能以便更清晰显示感兴趣的区域，Scope模块还可将数据存储在工作空间（Workspace）内。利用XY图形模型，可显示两个变量之间变化轨迹。2利用变量返回值利用Connection模块子库的Out模块，建立如图7.11的模型。仿真运行后，输出量yout保存在工作空间内，可供用户调用，绘制输出轨迹。图7.11 Out模块应用也可从Simulation菜单调用Simulation Parameters对话框的Workspace I/O页面，定义时间变量为tout，输出变量为yout。仿真结束后，在MATLAB COMMAND窗口下键入plot ( tout, yout )观察输出轨迹。3利用To Workspace模块建立如图7.12所示的SIMULINK模型，模型中采用两个To Workspace模块，一个用于记录时间变量t，一个用来记录变量y。变量名t和y通过打开To Workspace模块对话框填入。To Workspace只能接受一个输入变量。仿真结束后，在MATLAB COMMAND窗口下，键入plot ( t, y )可观察输出y的轨迹。仿真结果可以以数据文件的形式存起来，有两种方法：（1）利用MATLAB COMMAND窗口，键入 ty=