算法模拟

电池问题-随机模拟的方法论

字数统计: 2.1k字   |   阅读时长: 7分
, , ,

摘要: 随机模拟入门，参考《随机模拟方法与应用》

【对算法，数学，计算机感兴趣的同学，欢迎关注我哈，阅读更多原创文章】
我的网站：潮汐朝夕的生活实验室
我的公众号：算法题刷刷
我的知乎：潮汐朝夕
我的github：FennelDumplings
我的leetcode：FennelDumplings

随机模拟方法

人们面临的大量实际问题都存在不确定性。例如交通情况，商品库存，金融市场，服务排队，彩票中奖等等。这些问题的共同特点是需要用概率论来对这些系统进行建模，然后分析系统的行为特征。但是建立的随机模型的求解非常困难，甚至无法获得解析解。随机模拟方法在此时就很重要了。

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

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

这里推荐一本国内的书《随机模拟方法与应用》，非常适合入门，大致上会涉及以下一些内容。

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

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

问题: 电池问题

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

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

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

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

解析解

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

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

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

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

在图中，我们定义状态：

S0: 新电池开始供电，备用电池已就绪，也就是初始状态
S1: 新电池开始供电，备用电池开始充电，也就是过程中会不断重复的状态
S2: 无新电池可用，也就是供电停止的状态

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

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

解得 x = 14

随机模拟

step1: 描述系统

输入：新开始供电的电池的耗尽时间 r
状态：S0(2块备用电池就绪，也就是初始状态), S1(1块备用电池就绪，也就是正常换电池的状态), S2(无备用电池就绪，也就是停止供电的状态)
随机事件：换电池，ti 表示第 i 次随机事件(换电池)的时间点。
输出：第 m 次随机事件时，状态处于 S2，则 tm 为停止时间

step2: 设置变量

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

step3: 运行规则

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

step4: 模拟系统

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

代码如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
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: 抽样与统计

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

1
2
3
4
5
6
7
8
9
10
11
12
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]))

模拟结果

1
平均停止供电时间: 13.991030

step6: 解释结果

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

1
2
3
4
5
6
7
8
9
10
11
12
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()
  • 平均停止时间随着试验次数的变化图

  • 停止时间的分布图

随机模拟的方法论

在上面的例题中，我们已经实践了随机模拟的方法论。这里总结如下

  • 描述系统
    • 系统的输入、状态、输出。
    • 随机事件是什么
  • 设置变量
    • 为系统的输入、状态、输出设置变量
    • 为随机事件设定随机数
  • 运行规则
    • 系统状态如何变化
    • 随机数如何产生
  • 模拟系统
    • 给定系统初始情况
    • 给出系统的输出
  • 抽样与统计
    • 大量重复模拟试验
    • 对结果进行统计
  • 解释结果
    • 解释模拟结果
    • 必要时改变某些设定重新模拟
Share