matlab模拟技术精品资料.docx

《matlab模拟技术精品资料.docx》由会员分享，可在线阅读，更多相关《matlab模拟技术精品资料.docx（40页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、1.1.1 模拟技术在许多数学方法中，一般都要用到解析论证和数值计算的技巧。但是，许多现实的系统是很复杂的，其中的随机性因素往往难以用数学公式表示出来，也就很难使用数学推导或数值计算机的手段来分析系统、预测系统的性能。因此，产生了对系统进行模拟的技术。在使用计算机对系统进行模拟之前，一般要清楚模拟的一般步骤和方法。后面将从较小的模拟实例，对模拟技术进行简要的介绍。模拟过程的一般过程为：（1） 分析问题，收集资料。需要搞清楚问题要达到的目标，根据问题的性质收集有关随机性因素的资料。这里用得较多的知识为概率统计方面。在这个阶段，还应当估计一下待建立的模拟系统的规模和条件，说明哪些是可以控制的变量，

2、哪些是不可控制的变量。（2） 建立模拟模型，编制模拟程序。按照一般的建模方法，对问题进行适当的假设。也就是说，模拟模型未必要将被模拟系统的每个细节全部考虑。模拟模型的优劣将通过与实际系统有关资料的比较来评价。如果一个“粗糙”的模拟模型已经比较符合实际系统的情况，也就没有必要建立费时、复杂的模型。当然，如果开始建立的模型比较简单，与实际系统相差较大，那么可以在建立了简单模型后，逐步加入一些原先没有考虑的因素，直到模型达到预定的要求为止。编写模拟程序之前，要现画出程序框图或写出算法步骤。然后选择合适的计算机语言，编写模拟程序。模拟实现的工具较多，如数学软件类：Matlab、Mathematica、

3、Maple、MathCAD等，还有其它的高级语言工具，如：Visual C+、C+ Builder、Delphi、Borland C+3.1等。（3） 运行模拟程序，计算结果。为了减小模拟结果的随机性偏差，一般要多次运行模拟程序，还有就是增加模拟模型的时段次数。（4） 分析模拟结果，并检验。模拟结果一般说来反映的是统计特性，结果的合理性、有效性，都需要结合实际的系统来分析，检验。以便提出合理的对策、方案。以上步骤是一个反复的过程，在时间和步骤上是彼此交错的。比如模型的修改和改进，都需要重新编写和改动模拟程序。模拟结果的不合理，则要求检查模型，并修改模拟程序。1.1.2 模拟时间利用计算机进行模

4、拟时，有两种控制模拟时间的方法。一种是固定时间增量法，另一种是可变时间增量法，又叫面向事件法。固定时间增量法，是选用一段合适的时间作单位，然后每隔一个单位时间就计算一次有关参数的值，到达预定的模拟时间后，模拟程序结束。在编写这种程序时，一般可以建立一个“模拟时钟”变量。程序的主体框架一般时个大的循环，循环变量，则为模拟时间；在每个循环体内，就是对每个时段作处理。例如，有些排队论模型，可能就是以每隔一段时间（一天或者一个月）进行处理。采用可变时间增量法编写的模拟程序也有一个“模拟时钟”变量，但它是在一个事件发生时，“模拟时钟”才向前推进。需要注意的是，该模拟方法每一步经过的时间是可变的，而且会自

5、动寻找下一个最早使系统状态发生变化的事件。整个模拟直到“模拟时钟”到达指定的时间长度为止。可以参考有关离散系统仿真的内容。1.1.3 模拟语言运用计算机进行模拟，还要选用适当的程序设计的语言。当然，象C/C+、Pascal、Fortran这样的高级语言是可以用于模拟的实现，特别是在面向对象这样的程序设计思想的引入，使得采用这样的语言实现模拟要更方便些。但是，用这样的语言编写的程序一般都很长，而且编写复杂，调试费时。因此人们研究初了许多专门用于模拟的语言，如GPSS、SIMULA等。这里本书不讨论其他的模拟语言，而主要讲解如何使用Matlab语言编写模拟程序。Matlab不仅数值计算功能强大，而

6、且由于其语言的简洁和高级，编写的代码少，而且容易调试，实现模拟模型很快。1.1.4 随机数的模拟计算机模拟主要用于模拟复杂过程或现象的一些方法。采用计算机模拟一个实验或一个过程，那么用不同的数据重复计算机模拟就能得出统计学结论。使用这种研究方法得到的结论可能在数学上不很精确，但其精确性对于我们了解所模拟的过程已经足够了。考虑落在单位区间（0，1）中的一个实数序列。简单的说，如果这些数是杂乱地分布于整个区间中，且其排列次序似乎也无章法可循，那么这一序列就称为是随机地。下列序列的数就不是随机的：（1） 该序列的数是单调增加的或单调减小的；（2） 后一个数是关于前一个数的连续函数，如；1.1.5 随

7、机数的产生大多数计算机系统都有随机数生成器，MATLAB也有自己的随机数生成器，但是这些系统产生的随机数一般称为伪随机数，它们不是真正随机的。有许多产生随机数的算法，下面介绍一种算法还是比较令人满意的。产生均匀分布于开区间(0,1)中的随机数。算法如下：取一个整数，使得，对，计算：这里的是范围在中的整数。初始整数称为这个序列的种子。1和这个梅森(Mersenne)质数之间的任何一个整数都可取为种子。下面就用Matlab编写实现该算法的程序，函数名为mrand：function r=mrandglobal LL=mod(16807*L,2147483647);r=L*4.656612875245

8、9e-10;调用说明：如果要调用该函数mrand，需要在程序中添加如下2条命令：global LL=3第一句：global L 声明变量L为全局变量第二句：给L赋予初值下面集中随机数也是常用的，都可以借助上面的函数mrand实现。fix:取整函数，向靠近0的方向取整（1） 要产生在(a,b)上均匀分布的随机数x，则x=(b-a)*mrand+a（2） 产生集合0,1,2,n中的随机数I，则I=fix(n+1)*mrand)（3） 产生从j到k(jk)的随机整数X，则X=fix(k-j+1)*mrand)+j注意：一般说来直接采用Matlab中的rand就可以产生（0，1）之间均匀分布的随机数。

9、下面通过自定义的Matlab函数用于模拟随机变量，当然也可以改写为非函数实现。1.1.6 模拟均匀分布随机变量的函数function r=rnd_u(a,b)%产生在a,b间均匀分布的随机数r=a+(b-a)*rand;return1.1.7 模拟指数分布随机变量的函数function r=rnd_beta(lamada)%模拟指数分布%lamad表示指数分布的参数r = -log(rand)/lamada;return1.1.8 模拟正态分布随机变量的函数function r=rnd_normal(arg_mean,arg_segema)%arg_mean 均值%arg_segema 标准差

10、r = arg_mean + arg_segema*randn;%不是randreturn1.2 蒙特卡罗模拟法蒙特卡罗（Monte Carlo）方法，也称为随机模拟（random simulation）。基本思想：为了解决数学、物理、工程技术等方面的问题，首先建立一个概率模型或随机过程，使它的参数等于问题的解；然后通过对模型或过程的观察或抽样试验来计算所求参数的统计特征，最后给出所求解的近似值。1.2.1 模拟寻求近似圆周率平行四边形ABCD的面积为圆的面积现在模拟产生在正方形ABCD中均匀分布的点n个，如果这n个点中有m个点在该圆内，则圆的面积与正方形ABCD的面积之比可近似为m/n；即通

11、过模拟求圆周率的程序如下：n=yesinput(请输入产生点的个数：,500);m=0;for i=1:n if (-1+2*rand)2+(-1+2*rand)2=1 m=m+1; endendmypi=4*m/n下面是一组运行结果：n得近似100003.1028500003.13951000003.14182000003.1403注意：以上运行结果如果n得输入都相同，但结果都一般都会完全相同，这是由于产生的点是随机的，自然结果也就不同了。1.2.2 用蒙特卡罗法估算定积分例：求定积分。是该定积分的精确解。思想：对于，如果，则可以通过模拟的方法计算其定积分。构造一个矩形包含了曲边梯形，产生n

12、（足够大）个在举行区域内的点，如果落在由函数构成的曲边梯形内的点为m个，则所求定积分为求解程序：n=106;a=0;b=1;d=max(a,b)+1;m=0;for i=1:n, x=a+rand*(b-a); y=d*rand; if y=x2, m=m+1; end ends=m/n*d*(b-a)运行上面的程序，求得结果为 0.33321600000000这个结果非常接近于真实值。下面编写一个通用的程序，前提是被积函数通用函数：function s=cmmmonitixing(funname,a,b,d,n)if nargin=5, n=10000;%缺省endif n10000, n=

13、10000;%至少产生1万个点endif nargin4 d=max(a,b)+1;%缺省endm=0;%计数器for i=1:n, x=a+rand*(b-a); y=d*rand; if y=feval(funname,x), m=m+1; end ends=m/n*d*(b-a)%计算近似定积分为了计算上面的定积分，按如下步骤：第一步：编写被积函数的M文件function y=myfunx2(x)y=x2;第二步：调用上面的通用程序输入：cmmmonitixing(myfunx2,0,1)ans = 0.33680000000000cmmmonitixing(myfunx2,0,1,3,

14、100000)ans = 0.33429000000000通过上面的调用可以看出，当n取得越大，所计算得到的定积分更接近于真实值。1.2.3 用蒙特卡罗法估计体积有这样一个问题：位于锥面S1:上方和球面S2:内部的区域的体积。首先分析，这块区域包含于以和为界的盒子内，这个盒子的体积为8。理论上，所求区域的体积为球面S2面积的一般加上圆锥的体积，设为V.因此，需要生成在这个盒子内生成随机点，并将落在这块区域内的随机点个数于生成的随机点的总数之比乘以8（盒子的体积），即为所求区域的面积。实现程序如下：numinner=0;times=yesinput(输入模拟次数:,1000,100 10000)

15、for i=1:times, x=2.0*rand - 1.0; y=2.0*rand - 1.0; z=2.0*rand; if (x*x+y*y=z*z) & x*x+y*y=z*(2-z) numinner=numinner+1; end endsprintf(模拟到%d次后，近似体积%.5f,times,numinner/times*8)模拟结果：模拟到5000次后，近似体积3.064001.3 案例：渡口模型1.3.1 问题描述一个渡口的渡船营运者拥有一只甲板长32米，可以并排停放两列车辆的渡船。他在考虑怎样在甲板上安排过河车辆的位置，才能安全地运过最多数量的车辆。分析：怎样安排过河

16、车辆，关心一次可以运多少辆各类车。准备工作： 观察数日，发现每次情况不尽相同，得到下列数据和情况： (1) 车辆随机到达，形成一个等待上船的车列；(2) 来到车辆，轿车约占40，卡车约占55,摩托车约占5；(3) 轿车车身长为3.55.5米，卡车车身长为810米。1.3.2 问题分析这是一个机理较复杂的随机问题，是遵循“先到先服务”的随机排队问题。解决方法：采用模拟模型方法。因此需考虑以下问题：(1) 应该怎样安排摩托车？ (2) 下一辆到达的车是什么类型？(3) 怎样描述一辆车的车身长度？ (4) 如何安排到达车辆加入甲板上两列车队中的哪一列中去？本实验主要模拟装载车辆的情况，暂时不考虑渡船

17、的安全。1.3.3 模型建立设到达的卡车、轿车长度分别为随机变量。结合实际，这里不妨假设卡车、轿车的车身长度均服从正态分布。由于卡车车身长为810m，所以卡车车长的均值为m，由概率知识中的“”原则，其标准差为，所以得到。同理可得。1.3.4 模拟程序设计由以上的分析，程序设计时的应划分的主要模块（函数）如下：（1） 确定下一辆到达车辆的类型；（2） 根据车的类型确定到达车辆的长度；（3） 根据一定的停放规则，确定放在哪一列。1.3.5 模型求解结果及分析（一）运行结果程序名为sim_dukou，运行程序，输出结果如下：sim_dukou输入模拟次数：1000平均每次渡船上的车数mean_n =

18、 5.4840 3.9180 0.5160（二）结果分析上面为运行一次模拟程序，模拟次数为1000次的模拟结果。从模拟结果，你能得出什么结论？发现摩托车的平均数量不到1辆，因此从另外一方面看，忽略摩托车的长度是合理的。统计结果显示平均每次渡口时船上卡车、轿车、摩托车数量分别为5.484、3.918、0.516辆。1.3.6 模拟程序function sim_dukou%渡口模型的模拟%程序设计: 张勇%2005-2-1n=input(输入模拟次数：);if isempty(n) | (n500) n=500;endN=zeros(1,3);%依次为卡车数量、轿车数量、摩托车数量for i=1:

19、n, isfull=0; L=0 ,0;%第一列长度,第二列长度 while isfull, id=makeid; N(id) = N(id) + 1; newlen=getlength(id); isfull,pos=getiffull(L,newlen); if isfull, L(pos)=L(pos)+newlen; end%if end%whileend%fordisp(平均每次渡船上的车数)mean_n=N/nfunction r=makeid%模拟下一辆到达车的类型t=rand;if t=0.55, r=1;%到达卡车elseif tL(2), if L(2)+newlen32,

20、 pos=2; else full=1; endelse if L(1)+newlen32, pos=1; else full=1; end end1.3.7 思考题（1） 装载车辆时，如何做到“安全”？应该如何制定相应的策略？（2） 请将这些策略加入到模拟程序中，再观察结果。如果甲板长度及停放的列数不一定是2列，请如何修改程序使得程序更通用些？1.4 案例：核反应堆屏蔽层设计问题1.4.1 问题描述与分析核反应堆屏蔽层是用一定厚度的铅包围反应堆，用以阻挡或减弱反应堆发出的各种射线。在各种射线中，中子对人体伤害极大，因此，在屏蔽层的设计中，了解中子穿透屏蔽层的概率，对反应堆的安全运行至关重要。

21、问题背景 假定屏蔽层是理想的均匀平板。一个中子进入屏蔽层后运动的物理过程：中子以初速度和方向角射入屏蔽层，运动一段距离后与铅核发生碰撞，中子获得新的速度及方向。再游动一段距离后，与铅核发生第二次碰撞，并获得新的状态，如此等等，经过若干次碰撞后，出现下述情况之一时中子终止运动过程：（1）中子被弹回反应堆；（2）中子穿透屏蔽层；（3）第n次碰撞后，中子被屏蔽层吸收。为使屏蔽层的厚度达到安全设计要求，在计算机上对中子在屏蔽层的运动过程进行模拟。1.4.2 模型假设：1 假定屏蔽层平行板厚度为D=3d，其中d为两次碰撞之间中子的平均游动距离；2 假设在第10 次碰撞以后，中子速度下降到为某一很小数值而

22、终止运动（被引收）.因为每次碰撞后，中子因损失一部分能量而速度下降。3 假定中子在屏蔽层内相继两次碰撞之间游动的距离服从指数分布；4 中子经碰撞后的弹射角 U(0, 2).1.4.3 中子运动的数学描述弹射角i 第i 次碰撞后中子的运动方向与x 轴正向的夹角.xi 第i 次碰撞后中子所处位置与屏蔽层内壁的距离。Ri 中子在第 i 次碰撞前后的游动距离.中子在屏蔽层里随机游动，第 i 次碰撞以后，按照它的位置坐标 xi，可能有以下三种情况发生：（1）xi0，中子返回反应堆；（2）xiD，中子穿透屏蔽层；（3）0xiD，若iD); c(2) = c(2) + length(t); if lengt

23、h(t)0, xx=xx,x(t); yy=yy,y(t); x(t)=; y(t)=; end %c(1)返回反应堆数量 t=find(x0, xx=xx,x(t); yy=yy,y(t); x(t)=; y(t)=; end if j=N, c(3) = c(3) + length(x); end end xx=xx,x; yy=yy,y; check = sum(c) - n %right-0 check2=length(xx)-n%right-0 bili=c/n %return%将此行注释或删除后才能执行下列代码 %重绘所有中子，当n增大时，会花较长时间 plot(xx,yy,r.)

24、 hold on line(0,0,-H,H) hold on line(D,D,-H,H) text(-4*d,H-5*d,返回) text(-4*d,H-9*d,sprintf(%6.2f%),c(1)/n*100) text(D/2-d,H-5*d,吸收) text(D/2-2*d,H-9*d,sprintf(%6.2f%),c(3)/n*100) text(D+2*d,H-5*d,穿透) text(D+d,H-9*d,sprintf(%6.2f%),c(2)/n*100) hold off1.4.7 思考题请写出该模拟程序的算法描述或者流程图。1.5 案例：理发店系统研究一个理发店有两

25、位服务员A和B，顾客们随机到达店内，其中60的顾客仅需剪发，每位花5分钟时间，另外40顾客既要剪发又要洗发，每位用时8分钟。理发店是个含有多种随机因素的系统，请对该系统进行模拟，并对其进行评判。可供参考内容：“排队论”，“系统模拟”，“离散系统模拟”，“事件调度法”1.5.1 问题分析理发店系统包含诸多随机因素，为了对其进行评判就是要研究其运行效率，从理发店自身利益来说，要看服务员工作负荷是否合理，是否需要增加员工等考虑。从顾客角度讲，还要看顾客的等待时间，顾客的等待队长，如等待时间过长或者等待的人过多，则顾客会离开。理发店系统是一个典型的排队系统，可以用排队论有关知识来研究。1.5.2 模型

26、假设：1 60的顾客只需剪发，40的顾客既要剪发，又要洗发；2 每个服务员剪发需要的时间均为5分钟，既剪发又洗发则花8分钟；3 顾客的到达间隔时间服从指数分布；4 服务中服务员不休息。1.5.3 变量说明：u：剪发时间（单位：分钟），u=5m；v: 既剪发又理发花的时间（单位：分钟），v=8m；T： 顾客到达的间隔时间，是随机变量，服从参数为的指数分布，（单位：分钟）T0：顾客到达的平均间隔时间（单位：秒），T0；1.5.4 模型建立由于该系统包含诸多随机因素，很难给出解析的结果，因此可以借助计算机模拟对该系统进行模拟。考虑一般理发店的工作模式，一般是上午9：00开始营业，晚上10：00左右结

27、束，且一般是连续工作的，因此一般营业时间为13小时左右。这里以每天运行12小时为例，进行模拟。这里假定顾客到达的平均间隔时间T0服从均值3分钟的指数分布，则有3小时到达人数约为人，6小时到达人数约为人，10小时到达人数约为人，这里模拟顾客到达数为60人的情况。（如何选择模拟的总人数或模拟总时间）1.5.5 系统模拟：根据系统模拟的一般方法，需要考虑系统的如下数据、参数。1. 状态（变量）（1） 等待服务的顾客数；（2） A是否正在服务；（3） B是否正在服务；2. 实体：两名服务员、顾客们3. 事件： （1） 一名新顾客的到达；（2） A开始服务；（3） A结束服务；（4） B开始服务；（1）

28、 B结束服务；4. 活动：（1） 顾客排队时间（2） 顾客们到达的间隔时间（3） A的服务时间（4） B的服务时间；在系统模拟时，为了研究系统的整体情况，这里考虑顾客到达后不离开，且等待队长不限。要考虑如果服务员均空闲时，顾客先选择谁服务？要考虑模拟的时间设置还有顾客数目。模拟终止条件是根据顾客数目还是根据营业时间终止？1.5.6 系统模拟算法设计自行设计算法finished=0;初始化运行时钟while finished=0if 产生的顾客数不到规定数目时 then，产生该顾客的有关数据；将顾客加入等待队列；else运行时钟继续；endif处理服务员的状态（包括工作状态，空闲时间）；获得服务

29、员的服务优先顺序；根据服务员优先顺序从等待队列中安排服务；endwhile注：其它实现方法参考，离散系统仿真算法：事件调度法。1.5.7 系统模拟程序顾客到达的间隔时间T的计算机产生方法，利用T=，%理发店系统的模拟（案例分析之一）%关键词：面向事件的计算机模拟技术clear allcurclock=0;%当前时刻，动态变化totalcustomer=0;%总共服务的顾客数numsrv=2;srvstatus=zeros(numsrv,5);%服务员有关数据%srvstatus 第1列：服务状态（0空闲，1正在服务）；第2列：当前服务顾客编号；% 第3列：当前服务结束时刻；第4列：服务员空闲时

30、间；第5列：服务的顾客总数endtime =0;%结束时间waiting=;%等待队列数据%waiting 第1列：顾客编号；第2列：顾客到达时刻；第3列：顾客开始接受服务时刻；% 第4列：接受服务时间；第5列：顾客结束服务时刻；第6列：间隔时间cur=zeros(1,6);%当前产生顾客的数据，对应关系同waitingavgwaitlen=;%平均等待队长avgwaittime =;% 平均等待时间ujiange=5;%平均间隔时间finished=0;numsimucustumer=yesinput(输入等待模拟的顾客数：,10,10 1000);while finished=0, if

31、totalcustomer numsimucustumer %产生一个顾客的到达及其有关性质的数据 totalcustomer = totalcustomer+1; jiange= -log(rand)*ujiange;%与上一个顾客的到达的间隔时间 curclock = curclock + jiange; cur(1)= totalcustomer ;% 第1列：顾客编号 cur(2) = curclock;%第2列：顾客到达时刻cur(6) = jiange; 第6列：间隔时间%下面产生接受服务时间（可改进模型） if rand0.6, %产生顾客有关性质：这里是产生接受服务时间 cur

32、(4) = 5; else cur(4) = 8; end %放入等待队列 if isempty(waiting), waiting= cur; else m,n=size(waiting); waiting(m+1,:)= cur; end else curclock = curclock + (-log(rand)*ujiange); end%if totalcustomer %分配等待队列（看是否有服务员空闲，如果有则分配；否则继续执行）%处理服务员的服务状态 for i=1:numsrv, if srvstatus(i,1)=1 & srvstatus(i,3) curclock, s

33、rvstatus(i,4)= 0;%没有休息（正在忙） else srvstatus(i,4)= curclock-srvstatus(i,3);%目前已经空闲的时间 end end %处理服务员服务的先后顺序（依据空闲时间）（精细处理） tmp=srvstatus(:,4); for i=1:numsrv, value,id=max(tmp); b(i)=id; tmp(id)=0;%已经排序了 end %此时等待队列必然不为空 for j=1:numsrv, i=b(j);%确定服务员的序号 if(srvstatus(i,1)=0) %找一个顾客开始服务，同时计算该顾客什么时候接受服务，结

34、束服务； m,n=size(waiting); if m=0, break; end if waiting(1,5)=0,%还没有开始接受服务 waiting(1,3)= curclock; waiting(1,5)= waiting(1,3)+waiting(1,4);%结束时刻 srvstatus(i,1)=1;%设置为忙状态 srvstatus(i,2)=waiting(1,1);%顾客编号 srvstatus(i,3)= waiting(1,5);%结束时刻 srvstatus(i,5)=srvstatus(i,5)+1;%又服务了一个顾客 %计算等待时间 avgwaittime(en

35、d+1) = waiting(1,3)-waiting(1,2); disp(sprintf(间隔时间(%8.2f) 顾客编号：%5d 接受服务员（%4d）服务(到达时刻%10.2f),waiting(1,6),waiting(1,1),i,waiting(1,2) endtime=max(endtime,waiting(1,5) waiting(1,:)=;%从等待队列中离开 end end%if end%for m,n=size(waiting); %计算队长（这里的计算式子可以参考排队论有关术语进行确定） if totalcustomer =numsimucustumer,%队列为空，结

36、束 finished=1; end end%whiledisp(服务顾客数：)disp(srvstatus(:,5)disp(平均队长);disp(mean(avgwaitlen);disp(运行时间（分钟，小时）);disp(sprintf(%8.f%8.f,curclock,curclock/60);disp(平均等待时间（分钟）);disp(mean(avgwaittime );disp(结束时间（分钟）);disp(endtime );figure title(平均队长)bar(avgwaitlen) figure title(平均等待时间);bar(avgwaittime)1.6 实验题目1.6.1 实验：赶上火车的概率1.6.1.1 问题描述如图，一列火车从A站开往B站，某人每天赶往B站上这趟火车。 他已了解到：（1） 火车从A站到B站的运行时间是均值为30分钟，标准差为2分钟的随机变量；（2） 火车在下午大约1点离开A站，离开时刻的频率分布如下： 人火车BA出发时刻午后1:00午后1:05午后1:10频率0.70.20.1此人到达B站的时刻频率分布为：时刻午后1:28午后1:30午后1:32午后1:34频率0.30.40.20.1问他能赶上火车的概率是多少？1.6.1.2 变量