蒙特卡洛模型应用

普大帝

原发表于数学中国微信公众号，关注数学中国微信公众号可查看更多。请用电脑查看，手机版图片可能无法正常显示。
& k, O) k* U. i- L

很多实际问题要么没有具体求解方法，要么求解方法非常复杂，这时往往可以采取蒙特卡洛方法（即随机模拟方法）解决问题.例如，想求解一个不规则图形的面积，用矩形将不规则图形框住，将点随机地放入矩阵中，矩阵的面积乘以点落入不规则图形内的频率就是不规则图形面积的估计值.
0 }# H) B+ o' Z  U  P+ m1 }

蒙特卡洛是世界闻名的城市，蒙特卡洛方法借用了城市的名称，是以概率和统计理论方法为基础的一种计算机模拟方法，利用随机数来解决复杂的计算问题，将所求解的问题同一定的概率模型相联系，用计算机实现统计模拟或抽样，以获得问题的近似解.随机试验次数越多，获得的近似解的精度也越高.

由于要得到高精度的结果需要进行海量实验，所以传统蒙特卡洛方法一直没有得到广泛的应用.随着计算机技术的飞速发展，现在已经可以利用计算机来做大量实验，因此蒙特卡洛方法在近代得到了长足发展，在金融工程学、宏观经济学、计算物理学（如粒子输运计算、量子热力学计算、空气动力学计算）等领域都有广泛应用.

下面我们就以蒙特卡洛方法的几个经典例子来说明其思想，体会其用法.

2.2.1 投针算圆周率问题

问题描述 早在1777年，蒲丰（Georges Louis Leclere de Buffon,1707—1788）提出求解圆周率的一个另类思路，其过程是首先在纸上画一组间距为a的平行横线，然后随机往里投针，针长b（b＜a），再计算针与线相交的频率，从而估计出圆周率，示意图如图2.2所示.

图2.2 蒲丰投针计算示意图

思路分析 针是随机投到纸上的，其与线的夹角φ在0～π之间是等可能取到的.另外，针的中点到最近一条线的垂直距离x是在0~a/2上等可能取到的，从而可知，任一针的位置可用针与线的夹角φ和针的中心到最近一条线的垂直距离x来代表，记为（φ,x），显然（φ,x）在区域S上是均匀分布的.

模型建立 图2.2（a）中，φ为针与线相交的角度，x为针中点到与其最近的平行线的距离，显然，（x,φ）代表了针的角度以及与线的位置关系.注意到（x,φ）在区域S上是均匀分布的，而针与线相交的数学关系是.图2.2（b）所示是（x,φ）的变化区域.区域G上的曲线就是.而区域G的面积为

，而整个落针区域的面积为απ/2所以，点（x,φ）如果落在区域G，则针线相交，这样，利用几何概型，计算得到针线相交的概率P如下:p=2b/π×α

求出π的值为

其中N为总投针次数，n为针与平行线相交次数.以针线相交的频数估计针线相交的概率，从而得到要估计的量π.

模型求解 当初蒲丰的实验就是真实地投3 408次，并数出有多少针与平行线相交，得到相交的频率后代入计算式，从而求得π的估计值为355/113，很接近π了.读者可以亲自动手试验一下.

当然，我们完全可以利用计算机进行模拟实验.利用R程序，下载animation包，然后加载该包，利用其中的buffon.needle（）函数就可以随机模拟蒲丰投针试验以及相应的π估计值的变化.

buffon.needle（namx=50,interval=0）函数可以模拟500次蒲丰投针，结果如图2.3所示.

图2.3 蒲丰投针演示图

可以发现随着投针次数的增加，频率逐渐趋于概率，π的估计值也逐渐向π靠近，由于投针次数只有500次，结果与π尚有差距.随着投针次数逐渐增加，π的估计值会越来越接近真实值.理论上已经证明，超过80万次的投针试验，可以保证π的3位有效数字准确，也就是得到π为3.14.

2.2.2 交通路口堵车问题

城市交通是城市管理的重要部分，堵车问题是一个难治的“城市病”.要想制定治理对策，首先要了解堵车情况.

问题描述 如何通过蒙特卡洛方法来模拟一个车辆通过火车道口的交通情况?

思路分析 因为道路车辆的运行情况很随机，这时蒙特卡洛模型就是一个很有力的工具.一般道路的交叉路口以十字路口比较常见.十字路口有直行和左右转弯的车流.为了帮助读者了解这个问题，我们先从简单的情形开始思考，例如考虑只有直行的车流的情况.我们可以模拟一下，在路口亮红灯期间，会停下多少车?在下一次亮绿灯期间，这些车能不能及时通过道口?

假设以下情况：

（1）路口有来回两路车流，它们的到达数服从泊松分布，其泊松强度是一样的，如果不一样，取强度大的，记为a；

（2）有两个速度，车流在绿灯时通过路口的平均速度为v，红灯转为绿灯后，车辆由停转开，通过道口的平均速度为u，显然v＞u.假定u=50m/min，道口宽度为50m；

（3）红灯亮的持续时间为c，绿灯亮的持续时间为d（在模拟例子中，假设c为1min,d为5min）；

（4）黄灯亮的持续时间忽略，假定驾驶者看到红灯亮后可以立即停车；

（5）模拟开始时没有候车；

（6）模拟开始时刚由绿灯转为红灯；

（7）暂不考虑出现交通事故的可能.

模型建立 在亮绿灯期间，车子以速度v正常驶过路口，不会拥堵.亮红灯期间，车子停在路口，车子到达数按泊松分布，即在等候红灯的时间c里，到达路口并停下等候的车为n辆的概率为[插图],n=0,1,2,…，这里a是到达强度，需要根据路口情况事先确定.假如a=2,c=1，那么在亮红灯时间段，到达路口等红灯的车辆数的概率如表2.1所示.

表2.1 亮红灯时路口车量数目出现的概率分布

模型求解 此处仅给出求解流程，如下所示.

（1）抽取[0,1]之间均匀分布的随机数，确定这次模拟路口停红灯的车数，例如，抽到0.732，则这个数落在区间（0.671,0.857）的范围里，所以这次模拟停车数为3.

（2）计算红灯转为绿灯后，在绿灯延续期间d（如题设5min）内，这部车以速度u（如题设50m/min）通过道口（如题设宽度50m）共需时间t=（50/50）×3（min），如果t＞d，那么道口发生堵塞，在本次模拟中t=3min，没有发生堵塞.

（3）抽取随机数很多次，如10 000次，记下其中发生多少次堵塞，从而估算出路口发生堵塞的概率.

同样的思想可以模拟更复杂的情况，例如：

（1）一个多向交通交叉路口的车流情况（如十字路口），从而确定红绿灯的最佳交换时间；

（2）车辆到达强度是时间的函数，如早晚高峰时期，强度很大，夜间强度很低的情况；

（3）如果知道事故发生的概率，也可以模拟事故发生后，疏通交通需要的平均时间.

2.2.3 电梯问题

利用蒙特卡洛方法，可以对一些很难写出模型的表达式，或者表达式过于复杂，计算不可行的问题进行计算机随机模拟，从而得到在模型条件下的近似结果，以帮助我们了解实际问题.

问题描述 高层商务楼中一般配备了多部电梯，如何安排各部电梯的运行方式，使得既能保证大楼内各公司员工的正常工作和出行，又能降低能耗，节约成本?在一般高层商务楼中，经常采用的是分层或单双层的运行方式，或者某部电梯直达某层以上的方法，试建立一个适合的电梯运行方案（高峰时），并具体评价这些方案的优劣.

思路分析 评价电梯运行方案往往以电梯高峰期运行时间为依据.一般来说，可以预估电梯可能停靠楼层数、电梯运载次数、电梯停靠时间等参数来计算电梯高峰期运行总时间.但这种估计的方法十分粗略，可能与实际结果相差巨大.我们的目的是模拟电梯一次循环所需的平均时间，并设计电梯停靠方案以使这个时间最短.这里的主要随机量是各楼层乘客的到达数.因此可以考虑采用蒙特卡洛方法对电梯上下楼的方案进行随机模拟.

模型建立 在做出符合实际的模型假设后，对电梯上下楼的运行情况进行模拟.若参数设置合理，得到的结果将与实际情况十分吻合.对此，做出以下模型假设（文中的英文字符为参数）：

（1）高层商务楼一共有fl层，每层有people人上班；一共有n部电梯，每部电梯最多承载elvn人；

（2）假设1～4层为商场，顾客上下楼不通过商务电梯，即认为1～4楼没有人等电梯；

（3）上下班时，一段时间内等电梯的员工数服从泊松分布，设平均每分钟到达poisson人；

（4）因为刚上班时选择下楼的人数很少，可以忽略下楼人数对电梯运营的影响.上班时认为不会有人下楼，同理，下班时认为不会有人上楼；

（5）电梯从启动到匀速运行需要时间，电梯加速比匀速多用的时间记为elva，而电梯匀速经过一层楼需要elvrun秒（s），电梯加速需要elvrun-elva秒（s）；

（6）电梯开关门需要时间.假设电梯开关门时间一样，为elvclose秒（s）.员工进入电梯需要时间，假设每个人进入电梯需要elvout秒（s）；

（7）不会有两部以上的电梯同时抵达同一层接员工；

（8）上班时，当电梯运载超过90%的员工上楼后，则认为高峰期结束，下班同理.

首先对下班高峰时期进行分析.选择3种电梯运行方案：方案一为n部电梯每层楼都能抵达；方案二为单双层运行，即部分电梯只停单层，部分电梯只停双层；方案三为高低层电梯，部分电梯只达某高层以上，部分电梯只能到达某高层以下.

方案一的程序框图如图2.4所示，其中“时间”代表电梯由停靠到运行或者由运行到停靠所需要的时间.

若假设高层商务楼一共有25层，每层有100人上班，一共有6部电梯，每部电梯最多承载20人.电梯从启动到匀速需要时间，假设电梯能在一层楼的距离内由静止加速到匀速状态再减速到零，所需要的时间为1.5s，而电梯匀速经过一层楼需要1s，电梯开关门需要时间，开门需1s，关门需1s.员工进入电梯也需要时间，假设每个人进入电梯需要0.4s.

模型求解 采用方案一，编程输入相关参数即可对方案一进行模拟计算.

方案一基本编程思路如下.

由于电梯在上下行的过程中存在两种状态：停靠或运行，所以以电梯改变状态为循环依据.如果某电梯的状态即将发生改变，则一次循环结束.

图2.4 方案一程序框图

以此为思想，首先计算每个电梯到其状态改变所需要的时间，找出最短的时间p以及对应的电梯k；其次计算电梯i到新的改变状态需要的时间，其余电梯到其状态改变的时间要减去时间p.按以上方法进行循环，对电梯运载进行模拟，此外还需要考虑电梯载满客后去一楼等特殊情况，停止准则为电梯运载超过90%的员工进入电梯.

类似地，对于方案二和方案三也可以编写相应的MATLAB程序，只需要对电梯k将要抵达的楼层进行一定的限制即可.

此外，对于方案三，还需要找出最优临界层数使得电梯高峰期运行时间最短.对当前问题的设定数据，编程可得最优临界层数为16楼，即3部电梯负责5～15楼，3部电梯负责16～25楼.

结果解读 用MATLAB进行3次模拟，每次模拟对3种方案各模拟10 000次后取平均，单位为分钟，结果如表2.2所示。

表2.2 3种方案模拟结果

显然，对题设数据，方案三运行时间最少，电梯分高低层的方法值得采纳.方案一的优势在于每部电梯都能自由抵达各楼层，劣势在于低楼层员工不容易疏散，会出现高楼层的员工都下楼了而低楼层还存在很多等待的员工；方案二时间最多的原因是1～4层没有人，商务楼高25层，所以负责单数楼层的电梯要比负责双数楼层的电梯多运载一层的人数；方案三可以缓解方案一的劣势，让15楼以下的员工能先行离开.不过，还是会出现较多员工聚集在5楼或16楼的情况.

如果模型假设更贴合实际，利用蒙特卡洛方法可以得到更接近实际的结果.上面提供的MATLAB程序可以通过修改各参数，得到不同情况下方案一的所需时间.

2 a/ e: h% ?5 l9 Q
& j/ y9 `$ r( Q( ]; w: z
2 f" C3 Q" Y1 d7 B  J" J

, Y" O$ t3 Y: L4 H) O! w& T

' m* m  @9 r  ~