matlab随机模拟

MATLAB 中的计算机模拟 第一节 导言 计算机科学技术的迅猛发展，给许多学科带来了巨大的影响．计算机不但使问题 的求解变得更加方便、快捷和精确，而且使得解决实际问题的领域更加广泛．计 算机适合于解决那些规模大、难以解析化以及不确定的数学模型．例如对于一些 带随机因素的复杂系统，用分析方法建模常常需要作许多简化假设，与面临的实 际问题可能相差甚远，以致解答根本无法应用，这时模拟几乎成为人们的唯一选 择．在历届的美国和中国大学生的数学建模竞赛（MCM）中，学生们经常用到 计算机模拟方法去求解、检验等．计算机模拟(computer simulation)是建模过程 中较为重要的一类方法． 本章将讨论如何利用计算机技术对连续系统和离散系统进行模拟．由于计算机以 重复性运算见长，故它为研究模拟方法提供了极为合适的手段．计算机模拟是一 种广义的数值计算方法．通过本章的学习，你将会了解蒙特卡洛方法的思想；初 步掌握对连续系统或离散系统进行模拟的方法；掌握由实际问题怎样去建立计算 机模拟模型以及应用MATLAB编程语言进行计算． 第二节 引例 第三节 随机变量的抽样 第四节 连续系统的模拟 第五节 离散系统的模拟 第六节 范例 第七节 实验 第二节 引例：葡丰投针问题 在用传统方法难以解决的问题中，有很大一部分可以用概率模型进行描述．由于这类模 型含有不确定的随机因素，分析起来通常比确定性的模型困难．有的模型难以作定量分析， 得不到解析的结果，或者是虽有解析结果，但计算代价太大以至不能使用．在这种情况下， 可以考虑采用 Monte Carlo 方法。下面通过例子简单介绍 Monte Carlo 方法的基本思想． Monte Carlo方法是计算机模拟的基础，它的名字来源于世界著名的赌城——摩纳哥的 蒙特卡洛，其历史起源于 1777 年法国科学家蒲丰提出的一种计算圆周 π 的方法——随机投 针法，即著名的蒲丰投针问题。 1) Monte Carlo 方法的基本思想是首先建立一个概率模型，使所求问题的解正好是该模型 的参数或其他有关的特征量．然后通过模拟一统计试验，即多次随机抽样试验（确定 m 和 n）， 统计出某事件发生的百分比．只要试验次数很大，该百分比便近似于事件发生的概率．这实 际上就是概率的统计定义．利用建立的概率模型，求出要估计的参数．蒙特卡洛方法属于试 验数学的一个分支． MATLAB语言编程实现 l=1; n=1000; d=2； m=0； for k＝l：n x＝unifrnd（0，d／2）； p＝unifrnd（0，pi）； if )sin(15.0 yx ××< m=m＋1 elsc end end p=m/n pi_m=1/p 运行，即得结果． 蒙特卡洛方法适用范围很广泛，它既能求解确定性的问题，也能求解随机性的问题以及 科学研究中的理论问题．例如利用蒙特卡洛方法可以近似地计算定积分，即产生数值积分问 题． 任意曲边梯形面积的近似计算 一个古老的问题：用一堆石头测量一个水塘的面积．应该怎样做呢？测量方法如下：假 定水塘位于一块面积已知的矩形农田之中．如图 8．2 所示．随机地向这块农田扔石头使得 它们都落在农田内．被扔到农田中的石头可能溅上了水，也可能没有溅上水，估计被“溅上 水的”石头量占总的石头量的百分比．试想如何利用这估计的百分比去近似计算该水塘面 积？ 结合图 8．2 中的图形(1)分析，只要已知各种参数及函数（a，b，H，f(x)），有以下两种 方法可近似计算水塘面积． 1．随机投点法 1）赋初值：试验次数 n=0，成功次数 m=0；规定投点试验的总次数 N； 2）随机选择m个数对 ,1,, miyx ii << ，其中 Hybxa ii <<<< 0, ，置 n＝n＋l； 3）判断 ，若是，转 4，否则停止计算； Nn ≤ 4）判断条件 (示一块溅水的石头)是否成立，若成立则置m=m+1，转 2， 否则转 2； )( ii xfy < 5）计算水塘面积的近似值 NmabHS /)( ×−×= ． 2．平均值估计法 1）产生[a,b]区间的均匀随机数 ;,,2,1, Nixi "= 2) 计算 ;,,2,1),( Nixf i "= 3）计算 ∑ = −= N i ixfN abS 1 )()( 。 该方法的特点是估计函数 f(x)在[a,b]上的平均值，面积近似等于该平均值乘以(b-a)． 第三节 随机变量的抽样 根据随机变量遵循的分布规律使用一种方法获取它的具体值，叫做抽样(sampling).随机 变量的抽样方法很多，不同的分布采用的方法不尽相同。区间[0,1]上均匀分布的随机变量的 抽样是其中最基本的方法．在计算机上，其他各种分布的随机抽样方法都是在它的基础上产 生的． 计算机上产生的随机数是按照确定的算法产生的，它遵循一定的规律，显然不是真正随 机的，因此这种随机数我们叫做伪随机数(random number).只要伪随机数能通过一系列的统 计检验，就可以把它们当作真正的随机数放心地使用，而不会引起太大的误差． 产生均匀分布的伪随机数的常用方法有平方取中法、线性同余法和广义同余法等．由于 目前计算机上常用的高级语言(如C,Pascal,Fortran等)都有产生均匀分布随机数的系统函数， 我们可以直接使用而不必关心其实现原理，在此不作专门介绍．首先介绍MATLAB软件中 产生区间[0，1]上的均匀分布随机数的系统函数R=rand(n)，它产生 nn× 阶均匀随机矩阵．一 般R=rand(m,n)，产生 阶均匀随机矩阵． nm× 下面简单介绍由均匀随机数以及某种算法产生其他分布随机数的方法．用 表示 独立同分布于U[0，1]的随机数列． ",, 21 rr 直接抽样法 设连续的随机变量 X 有分布函数 ∫ ∞−= x dttfxF )()( 其中f（t）是X的密度函数。 因为 0＜F(x)＜1，令 r=F(x)，若函数 F(x)的反函数存在，对么对随机变量 X 的抽样可由 公式 产生出． )(1 rFx −= 例子 近似抽样法 近似抽样方法的步骤是设一组独立同分布的随机变量 利用中心极限定 理 "" ,,,, 21 nrrr 时当 ∞→ −∑ = nN nDr Err n i n i ii ),1,0(~ / 1 1 设 为[0,1]区间上的均匀随机数列，则",, 21 rr 12/1,2/1 == ii DrEr ，从而有，当 n 充 分大时， )1,0(~ 12 1 2 11 1 N n r n n i i∑ = − 常用的是 n＝12 的情形， ，其中 是正态 分布的 随机数，即由 12 个区间[0,1]上均匀随机数产生一个标准正态随机数．如果还要产生一般正 态分布的随机数，则只需用变换 ∑∑ = − = −=−= 6 1 122 12 1 12 6 i ii i i rrru 12u ασ += uz .另外，还有一种利用二维正态变换产生的标准 正态分布的随机数，其变换如下： 时当 0 2sinln2 2cosln2 1 21 21 ≠ ⎪⎩ ⎪⎨ ⎧ −= −= r rry rrx π π MATLAB 中各种分布下产生随机数的命令 常见的分布函数 MATLAB 语句 均匀分布 U[0,1] R=rand(m,n) 均匀分布 U[a,b] R=unifrnd(a,b,m,n) 指数分布 )(λE n)m,,exprnd(R λ= 正态分布 N(mu,sigma) R=normrnd(mu,sigma,m,n) 二项分布 B(n,p) R=binornd(n,p,m,n1) 了松分布 )(λP n)m,,poissrnd(R λ= 以上语句均产生 的矩阵． nm× 第四节 连续系统的模拟 系统模拟是在整个运行过程中对系统的仿真，是非常有效的和广泛使用的分析、研究复 杂系统的技术．在一定假设条件下，利用数学运算模拟系统的运行，称为数学模拟．现代的 数学模拟都是在计算机上进行的，因此也称为计算机模拟，简称模拟(simulation) 模拟分为静态模拟(static simulation)和动态模拟(dynamic simulation).数值积分中的蒙特 卡洛方法是典型的静态模拟．动态模拟又分为连续系统模拟和离散系统模拟． 状态随着时间连续变化的系统，称为连续系统(continuous system).对连续系统的计算机 模拟是近似地获取系统状态在一些离散时刻点上的数值．在一定假设条件下，利用数学运算 模拟系统的运行过程．连续系统模型一般是微分方程，它在数值模拟中最基本的算法是数值 积分算法．例如有一系统可用微分方程来描述 ),( ytf dt dy = 已知输出量 y 的初始条件 ，现在要求出输出量 y 随时间变化的过程 。最直 观的想法是：首先将时间离散化，令 00 )( yty = )(ty kkk tth −= +1 ，称为第 k 步的计算步距（一般是等间距 的），然后按以下算法计算状态变量 在各时刻 上的近似值 )(ty 1+kt ),)(,()( 111 kkkkkkk ttytfyyty −+=≈ +++ 其中初始点 按照这种作法即可求出整个 的曲线．这种最简单的数值 积分算法称为欧拉法．除此之外，还有其他一些算法． ",2,1),,( 00 =kyt )(ty 因此，连续系统模拟方法是：首先确定系统的连续状态变量，然后将它在时间上进行离 散化处理，并由此模拟系统的运行状态． 第五节 离散系统的模拟 离散系统(discrete system)是指系统状态只在有限的时间点或可数的时间点上有随机事 件驱动的系统．例如排队系统（queue system），显然状态量的变化只是在离散的随机时间点 上发生．假设离散系统状态的变化是在一个时间点上瞬间完成的． 为了模拟离散系统，必须设置一个模拟时钟(simulate clock)，它能将时间从一个时刻向 另一个时刻进行推进，并且能随时反映系统时间的当前值．其中，模拟时间推进方式有两种， 下次事件推进法和均匀间隔时间推进法．常用的是下次事件推进法．其过程是：置模拟时钟 的初值为 0，跳到第一个事件发生的时刻，计算系统的状态，产生未来事件并加入到队列中 去，跳到下一事件，计算系统状态，……，重复这一过程直到满足某个终止条件为止．为了 学习离散系统的模拟方法，举一个最简单的例子，以便帮助理解． 单服务排队系统 通常顾客到达时刻、顾客服务完毕并离去时刻等均视为随机事件（瞬间完成）。系统关 心的指标通常是：顾客的平均等候时间、服务效率等．定义程序事件为模拟运行到 150 个时 间单位（分钟）结束．该系统的模型一般用框图来描述，然后编制程序，模拟在一定时 间范围内系统运行的活动过程． 请看下列简图所表示的时间推进方式： 引入以下符号： 一第 i 个顾客到达的时刻； ix it 一相邻两个顾客到达的时间间隔 )( 1 iii xxt −= + 一第 i 个顾客接受服务的时间； is 一第 i 个顾客的排队等待时间； iD 一第 i 个顾客接受服务后离开的时刻ic )( iiii Dsxc ++= ； 在任意时刻 t，系统的状态可以用排队等候的顾客数目和服务员是否在工作来描述．排 队等候的顾客数目称为队长，记作 L(t)，为非负整数．服务员的状态用 S（t）表示，当服务 员工作时，令 S(t)=1；服务员空闲时，令 S(t)=0 系统的性能指标通常用排队长度、等待时间和服务利用率等来衡量．由于它们随时间改 变，一般用一段时间内的平均值作为数量指标．有以下三个指标： 1）平均队长 指队长 L(t)在[0,T]内的平均值，计算公式为 ∫= T dttLTL 0 )(1 2）顾客的平均等待时间 指每个顾客平均等待的时间长度，记作W ． 3) 服务利用率 指服务员工作时间在 T 中的比例， ∫= T dttSTU 0 )(1 为了简化问题，假设在上述模型下，系统的性能指标只有一个，即顾客的平均等待时 间．考虑用模拟方法来求W ，若系统能模拟出每位顾客的等待时间序列 ， 则 },,,{ 21 nDDD " ∑ = = n i iDn W 1 1 具体模拟步骤如下； 第 1 步 调查并收集和处理数据，记录顾客到达时刻、等待时间和服务时间．假定顾 客到达的间隔时间服从指数分布（均值为 10 分钟）；每个顾客的服务时间服从均匀分布 U[10,15] 第 2 步 构造模拟模型．输人因素：顾客的到达间隔时间和服务时间；排队规则：先 到先服务；一个服务机构． 第 3 步 模拟实验．设置模拟时钟及总的运行时间 T，如 8 小时等．推进原则按下次 事件推进或均匀间隔推进． 用MATLAB编制的程序 x=0;D=zeros(1,n)； t=exprnd(10,1,n); ％产生指数分布(0,1)的随机数，数组大小 n×1 s=unifrmd（10,15,1,n）； ％产生均匀分布U[10，15]的辩，数组大小 n×1 x=x+t(1) ％第一位顾客的到达时刻 for i＝2：n y=x+t(i) ％第 2－n 位顾客的到达时刻 j=i 一 l； c=x＋s(j) ＋D（j） %计算顾客离开时刻 if c＜y ％比较相邻两顾客的离开、到达时刻的大小 h＝0； else h＝c-y； end D(i)＝h； ％输出第 i 个顾客的等待时间 x=y； end E(D)=mean(D) ％计算数值的平均值 计算的一组结果如下表 8．1 近似的平均等待时间：E（D）=5.5l 67（分钟） 另一段更简洁的程序 function meantime=serve(n) arrive=zeros(1,0) for i=2:n arrive(i)=arrive(i-1)+exprnd(0.1); end wait=zeros(1,n) for i=1:n if (i= =1) wait(i)=0; else servetime=unifrnd(10,15); if(arrive(i-1)+servetime+wait(i-1)>arrive(i)) wait(i)=arrive(i-1)+servetime+wait(i-1)-arrive(i); else wait(i)=0; end end end meantime=mean(wait); 第六节 范例 追逐问题 1．问题提出 在图 8．4 中，假设正方形 ABCD 的四个顶点处各站一人．在某一时刻，四人同时以匀 速 v 沿顺时针方向追逐下一个人，并且在任意时刻他们始终保持追逐的方向是对准追逐目 标，例如，A 追逐 B，任意时刻 A 始终向着 B 追．可以证明四人的运动轨迹将按螺旋曲线 状汇合于中心 O． 怎样证明呢？有两种证明方法．一是分别求出四人的运动轨迹曲线解析式，求证四条曲 线在某时刻相交于一点．另一方法则是用计算机模拟将四人的运动轨迹直观地表示在图形 上． 2．建立模型及模拟方法 模拟步骤： 1）建立平面直角坐标系． 2）以时间间隔 t∆ 进行采样，在每一时t计算每个人在下一时t+ t∆ 时的坐标． 3）不妨设甲的追逐对象是乙，在时间t时，甲的坐标为 ，乙的坐 ),( ii yx 标 为 ),( 22 yx ． 甲 在 t+ t∆ 时 的 坐 标 为 ),sin,cos( 11 θθ tvytvx ∆+∆+ 其中 2122121212 )()(,sin,cos yyxxdd yy d xx −+−=−=−= θθ 同理，乙在 t+ t∆ 时的坐标为 )sin,cos( 22 θθ tvytvx ∆+∆+ ． 4）选取足够小的 ，模拟到t∆ tvd ∆< 时为止． 5）连接四人在各时刻的位置，就得到所求的轨迹． 连续系统模拟的特点是首先选定一个时间步长（通常是等间距的）；其次按时间顺序推 进，每推进一个时间步长，就对系统的活动和状态按预定的规则和目的进行考察。分析、计 算、记录，直到预定模拟结束条件（通常是时间条件）为止． 3．MATLAB 实现 根据以上模拟步骤，可编出 MATLAB 程序(simu2.m)如下： ％取 v=1,t=12,A,B,C,D 点的坐标分另为（0，10），（10，10），（10，0），(0, 0) v=1; dt=0.05; d=20； x=[0 0 0 10 10 10 10 0]； x（9）=x（l）； x（10）=x（2）； hold axis(‘equal’) axis([0 10 0 10]); for k=1：2：7 plot(x(k)，x(k+1)，’．’ ) end while(d＞0．1) for i=1：2：7 d=sqrt（（x（i）-x(i+1)）^2+(x(i+1)-x(i+3))^2）; x(i)=x(i)+v*dt*(x(i+2)-x(i))/d; x(i+1)=x(i+1)+v*dt*(x(i+3)-x(i+1))/d; plot(x(i),x(i+1),’.’) end x（9）= x（l）；x（10）= x（2）； end hold 运行上述程序(simu2, 回车)可得到图 8．4． 企业生产的库存系统 1．问题提出 在销售部门、工厂等领域中都存在库存问题，库存太多造成浪费以及资金积压，库存少 了不能满足需求也造成损失，部门的工作人员需要决定何时进货，进多少，使得所花费的平 均费用最少，而收益最大，这就是库存问题． 某企业生产易变质的产品．当天生产的产 品必须当天售出，否则就会变质．该产品单位成本为 2．5 元，单位产品售价为 5 元．企业 为避免存货过多而造成损失，拟从以下两种库存中选出一个较优的方案． 方案甲：按前一天的销售量作为当天的货存量； 方案乙：按前二天的平均销售量作为当天的货存量； 假定市场对该产品的每天需求量是一个随机变量，但从以往的统计分析得知它服从正态 分布 )4.22,135( 2N 2．建立数学模型 计算机模拟的基本思路： 第一，需要获得市场对该产品的需要量的数据； 第二，计算出按照两种不同方案经 T 天后企业的利润值； 第三，比较大小，从中选出一个较优的方案． 为了增加可读性和方便地使用 MATLAB 语言编程，各变量的符号解释如下： D 一一一每天的需求量； Q1 一方案甲当天的货存量； Q2 一方案乙当天的货存量; 一方案甲前一天的销售量； 1S 一方案乙前一天的销售量； 21S 一方案乙前二天的销售量； 22S 一方案甲当天实际销售量； 3S 一方案乙当天实际销售量； 4S 一方案甲当天的利润； 1L 一方案乙当天的利润； 2L 一方案甲累计总利润； 1TL 一方案乙累计总利润； 2TL T 一一一预定模拟天数． 建立数学模拟框图，如图 8．5． 3．求解 编制 MATLAB 程序(kucum.m) function [TL1，TL2]=kucum（T，S1，S21，S22） TL1=0；TL2=0；k=l； while k＜T Q1=S1； Q2=（S21＋S22）／2； D=normrnd（135，22．4）； If D＜Q1 S3=Q1; else S3=D； end if D＜Q2 S4=Q2 else S4=D； end L1=5*S3—2．5*Q1； L2=5*S4 一 2．5*Q2； TL1=TL1+L1； TL2=TL2+L2; K=k＋1； end S1=S3； S22=S21； S21=S4； 运行程序(kutest．m) T=10；S1=136；S21=136；S22=148； 「TLI，TL2」=kucun（T，S1，S21，S22） 计算结果： TL1=3．5289e＋003 TL2=3．5823 e＋003 反复运行程序，因为多数情形是 TL2＜TL2，所以认为方案乙较优． 第七节 实验 实验一：导弹跟踪问题 某军一导弹基地发现正北方向 120 千米处海面上有一艘敌艇以 90 千米/小时的速度向正 东方向行驶．该基地立即发射导弹跟踪追击敌艇，导弹速度为 450 千米/小时，自动导航系 统使导弹在任一时刻都能对准敌艇． l）试问导弹在何时何处击中敌艇？ 2）如果当基地发射导弹的同时，敌艇立即由仪器发觉．假定敌艇为高速快艇，它即刻 以 135 千米/小时的速度向与导弹方向垂直的方向逃逸，问导弹何时何地击中敌艇？ 3）敌艇与导弹方向成何夹角逃逸才好？从结论中你能得到些什么启示？ 思考：敌艇的逃跑策略是什么？在什么样的情形下可能逃脱？ 实验二：盐水浓度问题 某水池有 2000 立方米的水，其中含盐 2 千克．以每分钟 6 千克的速度向池中注人含盐 量为 0．5 千克／立方米的盐水，同时又以每分钟 4 立方米的速度从水池流出搅拌均匀的盐 水．每隔 10 分钟计算水池中水的体积、含盐量和含盐率，列出一表．从表中查出含盐量达 到 0．2 千克／立方米时用去多少时间． 实验三：定货策略问题 在物资的供应过程中，由于到货与销售不可能做到同步、同量，故总要保持一定的库存 储备．如果库存过多，就会造成积压浪费以及保管费用的上升；如果库存过少，会造成缺货．如 何选择库存和订货策略，就是一个需要研究的问题．现要研究以下问题： 某自行车商店的仓库管理人员采取一种简单的订货策略，当库存降低到 P 辆自行车时 就向厂家订货 Q 辆．如果某一天的需求量超过了库存量，商店就有销售损失和信誉损失， 但如果库存量过多，将会导致资金积压和保管费增加．若现在已有如表 8．2 中的五种库存 策略，试比较选择一种策略以使花费最少．已知该问题的条件为 l）从发出订货到收到货物需隔 3 天； 2）每辆自行车保管费为 0．75 元／天，每辆自行车的缺货损失为 1．80 元／天， 每次的订货费为 75 元； 3）每天自行车的需求量服从 0 到 99 之间的均匀分布； 4）原始库存为 115 辆，并假设第一天没有发出订货． 实验四：机器看管系统 一个机器看管系统有 m 台机器，并由 c 个工人共同负责看管与修理．并假设 1）各台机器的质量相同，机器的连续运转时间相互独立且服从同一负指数分布，平均 寿命为 ; )0(/1 >vv 2）每个工人技术相同，且修理时间相互独立并服从同一负指数分布，平均修理时间为 )0(/1 >uu 3）工人对故障机器的修理与其他机器连续运转是否正常无关，修复后的机器寿命分布 与新的一样； 4）机器停止运转每单位时间的损失费为 元，工人单位时间的产值为 元． 1C 2C 若机器的等待时间为 E ，工人总的空闲时间为 F ，则系统总的损失费为 ．试求当机器数 m 固定时，为使系统的总损失费最小，应配 FCECS ** 21 += 备多少工人为最优？假设已知 m=86，1/v=500 小时，1/u=34 小时， 46.31 =C 元／小时， 元／小时 2.32 =C