2.随机模拟的方法论-灵析社区

由于随机模拟往往可验证概率问题的结果，因此它在后面解题中基本都会用到。

本章首先介绍随机模拟方法的背景，然后以一个例题来阐述实际场景中随机模拟怎么用、具体流程有哪些，最后进行总结。

2.1随机模拟方法

大量实际问题都存在不确定性，如：交通路况、商品库存变动、金融市场交易量变化、服务排队情况等。这些问题的共同点是要求我们用概率方法对系统进行建模，然后分析系统的行为特征。但要求解这类随机模型非常困难，我们甚至无法获得解析解。而随机模拟对于解决上述问题就显得尤为重要了。

随机模拟方法是通过仿真随机系统的运行来获得系统的状态变化与输出结果的大量数据，进而对所得数据进行统计分析，估算系统行为的某些特征，并将估计的误差控制在一定范围内。这里的仿真其实是另一个单独的话题，里面也有很多分类，在这里我们只关注对随机系统的仿真。

随机模拟方法除了在解析解很困难或没有解析解的情况下是很有力的方法，在有解析解的时候，还可以用于验证解析解的正确性。

随机模拟有一个基本的方法论，它是一种用随机模拟解决问题的一个思路框架。有了这个方法论，拿到实际问题的时候就可以快速形成思路。这其中会分为描述系统、设置变量、运行规则、模拟系统、抽样与统计、解释结果这几步。
正确地描述系统的输入输出，以及运行规则，需要对系统进行数学建模。数学建模常见的模型有图论模型（例如二分图匹配）、概率统计模型（例如逻辑回归模型）、动态模型（例如常微分方程）、优化模型（例如线性规划）等。
对系统建模是需要准确找到随机事件是什么。这需要一些数学基础，主要是概率论和随机过程。
对系统中的随机事件，需要正确地为该随机事件配置随机数，这需要随机数采样的技术，常见的连续型和离散型分布可以直接生成，而一些特殊分布需要一些算法，例如逆变换法，接受拒绝法，抽样多维联合分布法。
蒙特卡洛法是对系统进行模拟的重要方法，主要包括马尔可夫链蒙特卡洛法，蒙特卡洛优化，蒙特卡洛积分。
随机服务系统，随机游走，以及元胞自动机是随机模拟的比较大的应用。

下面我们以一个实际问题来看一下随机模拟的方法论。

2.2随机模拟方法论实践：电池问题

考虑一个由充电电池构成的供电系统，一共有两个电池和一个充电器。其中一个电池给设备供电，另一个电池备用。

电池的耗尽时间为 1, 2, 3, 4, 5, 6 小时的其中一种情况，并且随机。耗尽的电池充满电需要 2.5 小时。

初始状态两个电池都是充满电的。

问：设备可以持续工作多长时间。

解析解

刚开始的时候，设备上的电池无论随机到的耗尽时间是多少，都会继续供电，因为备用电池在初始的时候是就绪的。

第一次更换电池后，新电池开始供电，备用电池开始充电。此时新电池随机到的耗尽时间就很关键了。

由于充电时间是 2.5，因此如果随机到的耗尽时间为 1, 2，那么耗尽后供电就会终止。如果随机到的耗尽时间是 3, 4, 5, 6，那么耗尽时备用电池已就绪，可以回到前面刚刚更换电池的状态，也就是【新电池开始供电，备用电池开始充电】的状态。

这个过程可以用有向图建模

在图中，我们定义状态：

S0: 新电池开始供电，备用电池已就绪，也就是初始状态

S1: 新电池开始供电，备用电池开始充电，也就是过程中会不断重复的状态

S2: 无新电池可用，也就是供电停止的状态

状态的转移概率，以及状态转移时对应的期望时间间隔标注在图中。

定义 x 为从初始到停止供电的总时间的期望，y 为从第一次换电池到停止供电的时间的期望。根据图中的状态转移关系，我们可以写出

随机模拟

随机模拟大致分为六个步骤。

step1: 描述系统

输入：新开始供电的电池的耗尽时间 r。

状态：S0（2 块备用电池就绪，也就是初始状态），S1（1 块备用电池就绪，也就是正常换电池的状态），S2（无备用电池就绪，也就是停止供电的状态）。

随机事件：换电池，ti 表示第 i 次随机事件（换电池）的时间点。

输出：第 m 次随机事件时，状态处于 S2，则 tm 为停止时间

step2: 设置变量

r 是 1 ~ 6 的离散均匀分布。用随机数发生器产生一系列随机数作为输入。

step3: 运行规则

产生一系列随机数后，随机事件的发生也就取决于这些随机数。

step4: 模拟系统

首先导入必要的包。其中 font 是字体对象，在 step6 中画图的时候需要用到。

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.font_manager import FontProperties
font = FontProperties(fname="/home/ppp/anaconda3/envs/python-3.6/lib/python3.6/site-packages/matplotlib/mpl-data/fonts/ttf/simhei.ttf")
import seaborn as sns

初始 t0 = 0, s = 0，然后模拟并产生输出即可。

代码如下

class PowerSupplySystem:
    def __init__(self, s0=0, t0=0.0, N=int(1e5)):
        self.N = N
        self.charging_time = 2.5
        self.s0 = s0
        self.t0 = t0
        self.reset()

    def reset(self):
        self.t = self.t0
        self.s = self.s0
        self.rs = np.random.randint(1, 7, self.N)
        self.i = 0

    def exhaust_time(self):
        nxt_t = self.rs[self.i]
        self.i += 1
        return nxt_t

    def __call__(self):
        np.random.seed()
        self.t += self.exhaust_time()
        self.s = 1
        while self.s != 2 and self.i < self.N:
            t = self.exhaust_time()
            self.t += t
            if t < self.charging_time:
                self.s = 2
        tmp = self.t
        self.reset()
        return tmp

system = PowerSupplySystem() 创建一个供电系统实例后，调用该实例后，会进行一次试验。并返回停止供电时间。

step5: 抽样与统计

重复模拟试验，对结果进行统计。

system = PowerSupplySystem()
T = int(1e5)
ts = np.zeros(T)
for i in range(T):
    ts[i] = system()

avg_ts = np.zeros(T)
avg_ts[0] = ts[0]
for i in range(1, ts.shape[0]):
    avg_ts[i] = (avg_ts[i - 1] * i + ts[i]) / (i + 1)

print("平均停止供电时间: {:.6f}".format(avg_ts[-1]))

模拟结果

平均停止供电时间: 13.991030

step6: 解释结果

下面我们画一下平均停止时间随着试验次数的变化图，以及停止时间的分布图。

fig = plt.figure()
ax1 = fig.add_subplot(1, 2, 1)
ax2 = fig.add_subplot(1, 2, 2)
ax1.plot(avg_ts)
ax1.set_title("供电系统", fontproperties=font)
ax1.set_xlabel("试验次数", fontproperties=font)
ax1.set_ylabel("平均断电时间", fontproperties=font)
ax2.hist(ts, bins=int(ts.max()))
ax2.set_title("供电系统", fontproperties=font)
ax2.set_xlabel("断电时间", fontproperties=font)
ax2.set_ylabel("频率", fontproperties=font)
plt.show()

平均停止时间随着试验次数的变化图

停止时间的分布图

2.3随机模拟的方法论总结

通过电池问题，我们简单实现了随机模拟过程，现进一步总结如下：

描述系统

系统的输入、状态、输出。

随机事件是什么

设置变量

为系统的输入、状态、输出设置变量

为随机事件设定随机数

运行规则

系统状态如何变化

随机数如何产生

模拟系统

给定系统初始情况

给出系统的输出

抽样与统计

大量重复模拟试验

对结果进行统计

解释结果

解释模拟结果

必要时改变某些设定重新模拟