关于模拟方法建模的一份介绍

有些时候，我们无法直接收集大量数据，即对于对象的行为直接观测或重复实验可能是不行的，所以此时就需要通过模拟的技术来收集数据，然后建模。这样的一种策略就是模拟方法建模，而模拟方法建模中最常用的一种方法就是蒙特卡洛法。因此，在这篇文章中将详细介绍蒙特卡洛法，以及运用蒙特卡洛法会用到的随机技术等内容。

一、蒙特卡洛算法

首先我们来认识蒙特卡洛算法。假如存在这样一种情况，在如图的一个湖蓝色方框中有一深蓝色椭圆，我们欲在不用积分以及不具体测量椭圆相关长度等情况下求得椭圆的面积，如何求？图如下：

此时，我们就可以用蒙特卡洛算法，即往图中随机增加足够的点，然后对于点数进行统计，然后会得到两个值，分别为落在矩形中的总点数，以及落在椭圆中的部分点数，而通过将这二者作比得到一个比值然后用这个比值乘以矩形总面积即可得到一个椭圆面积的近似答案。

在这里我可以给一个这个方法的python代码：

import random
import math
import matplotlib.pyplot as pltdef estimate_ellipse_area(a, b, num_points):# 矩形的宽度和高度rect_width = 2 * arect_height = 2 * brect_area = rect_width * rect_height# 初始化计数器points_in_ellipse = 0# 存储点的坐标，用于绘图points_x = []points_y = []for _ in range(num_points):# 生成随机点x = random.uniform(-a, a)y = random.uniform(-b, b)# 判断点是否在椭圆内if (x**2 / a**2) + (y**2 / b**2) <= 1:points_in_ellipse += 1points_x.append(x)points_y.append(y)# 计算椭圆面积的估计值ellipse_area_estimate = (points_in_ellipse / num_points) * rect_areareturn ellipse_area_estimate, points_x, points_y# 参数设置
a = 2  # 椭圆的半长轴
b = 1  # 椭圆的半短轴
num_points = 10000  # 随机点的数量# 运行模拟
ellipse_area, points_x, points_y = estimate_ellipse_area(a, b, num_points)# 输出结果
print(f"椭圆的面积估计值: {ellipse_area:.4f}")# 绘制结果
plt.figure(figsize=(8, 6))
plt.scatter(points_x, points_y, s=1, c='blue', label='Points in Ellipse')
plt.axhline(y=b, color='r', linestyle='--')
plt.axhline(y=-b, color='r', linestyle='--')
plt.axvline(x=a, color='r', linestyle='--')
plt.axvline(x=-a, color='r', linestyle='--')
plt.gca().set_aspect('equal', adjustable='box')
plt.title(f"Monte Carlo Estimation of Ellipse Area (a={a}, b={b})")
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.show()

输出结果与图像分别为：

椭圆的面积估计值: 6.3240

同样的，我们还可以将这个方法放在三维中去考虑，即在一个正方体盒子中，有一个球体，欲求得球体的体积，则可以用蒙特卡洛方法求得。具体步骤不再赘述，但可以参考这份python代码：

import random
import math
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3Ddef estimate_sphere_volume(r, num_points):# 正方体的体积cube_volume = (2 * r) ** 3# 初始化计数器points_in_sphere = 0# 存储点的坐标，用于绘图points_x = []points_y = []points_z = []for _ in range(num_points):# 生成随机点x = random.uniform(-r, r)y = random.uniform(-r, r)z = random.uniform(-r, r)# 判断点是否在球体内if x**2 + y**2 + z**2 <= r**2:points_in_sphere += 1points_x.append(x)points_y.append(y)points_z.append(z)# 计算球体体积的估计值sphere_volume_estimate = (points_in_sphere / num_points) * cube_volumereturn sphere_volume_estimate, points_x, points_y, points_z# 参数设置
r = 1  # 球体的半径
num_points = 10000  # 随机点的数量# 运行模拟
sphere_volume, points_x, points_y, points_z = estimate_sphere_volume(r, num_points)# 输出结果
print(f"球体的体积估计值: {sphere_volume:.4f}")# 绘制结果
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')# 绘制落在球体内的点
ax.scatter(points_x, points_y, points_z, s=1, c='blue', label='Points in Sphere')# 绘制正方体的边界
ax.set_xlim([-r, r])
ax.set_ylim([-r, r])
ax.set_zlim([-r, r])# 设置标签
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
ax.set_title(f"Monte Carlo Estimation of Sphere Volume (r={r})")# 显示图例
ax.legend()# 显示图像
plt.show()

结果与图像分别为：

球体的体积估计值: 4.1600

二、随机数的生成

在刚才的算法中，我们发现算法的一个核心点就是需要生成随机数，但我们知道计算机只能执行确定性算法，它无法产生随机性事件，所以我们只能通过生成伪随机数的方式来间接实现，至于伪随机数的生成有两种，分别是平方取中方法与线性同余方法。下面我来具体介绍下这两种方法：

2.1 平方取中方法

平方取中方法是Von Neuman在1946年在曼哈段项目中提出的，具体步骤如下：

1. 从一个4位数 $x_0$ 开始，这个数字称为种子；

2. 将它平方，即 $x_{0}^2$ ，此时得到一个8位数；

3. 取这个8位数的中间四位数作为下一个随机数；

4.重复上述步骤。

按照这样的一个方法就会得到一个数列，我们在取种子数时，虽然不一定非要取4位数，即使4位数是这个算法最经典的取法，但我们要注意我们要尽量取偶数位的种子数，以方便在第三步中取中间数这一操作。

那么，接下来，我将给出一个该方法的python代码：

def middle_square(seed, num_iterations):results = []for _ in range(num_iterations):# 平方square = seed ** 2# 转换为字符串并填充到8位square_str = f"{square:08d}"# 取中间四位new_seed = int(square_str[2:6])results.append(new_seed)seed = new_seedreturn results# 选择种子数
seed = 1234
num_iterations = 10# 生成随机数列
random_numbers = middle_square(seed, num_iterations)# 输出结果
print(f"种子 {seed} \n生成的随机数列: {random_numbers}")

其输出结果为：

种子 1234
生成的随机数列: [5227, 3215, 3362, 3030, 1809, 2724, 4201, 6484, 422, 1780]

我们可以更换种子数以及重复步骤的次数，将之更换为 $x_0 = 2041$ 与12，那么可以得到如下这个输出结果：生成的随机数列: [1656, 7423, 1009, 180, 324, 1049, 1004, 80, 64, 40, 16, 2]

观察这个数列我们将会发现这个数列在逐步趋近于0，而这就是平方取中法的一大缺点——它会将生成的伪随机数退化为0并永远停止在那。

我们可以再更改代码的数值，回到一开始的值上，种子仍然是1234，但重复次数改为100，即得到随机数列如下：

生成的随机数列: [5227, 3215, 3362, 3030, 1809, 2724, 4201, 6484, 422, 1780, 1684, 8358, 8561, 2907, 4506, 3040, 2416, 8370, 569, 3237, 4781, 8579, 5992, 9040, 7216, 706, 4984, 8402, 5936, 2360, 5696, 4444, 7491, 1150, 3225, 4006, 480, 2304, 3084, 5110, 1121, 2566, 5843, 1406, 9768, 4138, 1230, 5129, 3066, 4003, 240, 576, 3317, 24, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]

我们即可清晰发现它在后面一段彻底退化为0了，并且在之后的数列里也全部是0.

2.2 线性同余

线性同余方法是D.H.Lehmer于1951年引入的，而现在的大部分伪随机数列都是基于这个算法。它比起其他算法而言，它可以选择种子以生成循环的模式，而且这个循环往往很长，以至于大多数应用在大型计算机上都不会出现自身的重复。具体的，它的步骤如下：

1. 选择三个整数a，b，c，以及一个初始种子 $x_0$ ；

2. 将刚刚选择的数字按如下这个规则生成数列：

$x_{n+1}=(a*x_n +7)mod(c)$

如下是一个线性同余方法的python代码：

def linear_congruential(seed, a, b, c, num_iterations):results = []x = seedfor _ in range(num_iterations):x = (a * x + b) % cresults.append(x)return results# 参数设置
seed = 7
a = 1
b = 7
c = 10
num_iterations = 100# 生成随机数列
random_numbers = linear_congruential(seed, a, b, c, num_iterations)# 输出结果
print(f"线性同余方法生成的随机数列: {random_numbers}")

输出结果如下：

线性同余方法生成的随机数列: [4, 1, 8, 5, 2, 9, 6, 3, 0, 7, 4, 1, 8, 5, 2, 9, 6, 3, 0, 7, 4, 1, 8, 5, 2, 9, 6, 3, 0, 7, 4, 1, 8, 5, 2, 9, 6, 3, 0, 7, 4, 1, 8, 5, 2, 9, 6, 3, 0, 7, 4, 1, 8, 5, 2, 9, 6, 3, 0, 7, 4, 1, 8, 5, 2, 9, 6, 3, 0, 7, 4, 1, 8, 5, 2, 9, 6, 3, 0, 7, 4, 1, 8, 5, 2, 9, 6, 3, 0, 7, 4, 1, 8, 5, 2, 9, 6, 3, 0, 7]

我们观察这个结果，不难发现其中出现了一个问题，就是出现了循环。不过在这种方法里，循环周期由至多c个数得以保证，但是c可以选得足够大，并且能够选择a，b使得循环出现之前得到全部c个数，而减少这种现象造成的影响。

三、随机行为的模拟

对于随机行为模拟，我们可以先看一个简单的，就是对于掷硬币的模拟，在这个模拟中，我们可以生成一个0到1的随机数，然后如果这个随机数x， $0\leq x\leq 0.5$ ，则设它为正面，其余的则为反面，并不断重复这个过程，于是就可以模拟出掷硬币的这一行为。

接下来，让我们对更复杂的行为进行模拟，我们将模拟一个简单的城市交通系统，其中包括车辆在道路上的行驶、交通灯的控制以及行人过马路的行为。这个模拟将帮助我们评估交通系统的性能，例如平均等待时间、交通流量等。

在这个问题中，将包含以下元素：

道路：有几条主要的道路连接不同的区域。

交通灯：控制路口的交通流量。

车辆：在道路上行驶。

行人：在人行横道上过马路。

对于这个模拟最后的目标为：通过蒙特卡洛模拟，评估交通系统的性能指标，如平均等待时间、交通流量等。

关于模拟的具体步骤如下：

1. 初始化系统，设置道路、交通灯、车辆和行人的初始状态；

2. 制造随机性，随机生成车辆和行人的到达时间；

3. 模拟交通流动，根据交通灯的状态和车辆、行人的行为等来实现；

4. 记录统计数据，记录车辆的等待时间、交通流量等数据；

5. 分析结果，分析统计数据，评估交通系统的性能。

这样，就可以得到如下的一份代码：

import random
import simpy# 定义交通灯类
class TrafficLight:def __init__(self, env, green_duration, red_duration):self.env = envself.green_duration = green_durationself.red_duration = red_durationself.state = 'red'self.action = env.process(self.run())def run(self):while True:if self.state == 'red':self.state = 'green'yield self.env.timeout(self.green_duration)else:self.state = 'red'yield self.env.timeout(self.red_duration)# 定义车辆类
class Vehicle:def __init__(self, env, name, traffic_light, arrival_time, travel_time):self.env = envself.name = nameself.traffic_light = traffic_lightself.arrival_time = arrival_timeself.travel_time = travel_timeself.wait_time = 0self.action = env.process(self.run())def run(self):yield self.env.timeout(self.arrival_time)print(f'车辆 {self.name} 到达路口 at {self.env.now}')while self.traffic_light.state == 'red':yield self.env.timeout(1)self.wait_time += 1yield self.env.timeout(self.travel_time)print(f'车辆 {self.name} 通过路口 at {self.env.now} (等待时间: {self.wait_time})')# 定义行人类
class Pedestrian:def __init__(self, env, name, traffic_light, arrival_time, crossing_time):self.env = envself.name = nameself.traffic_light = traffic_lightself.arrival_time = arrival_timeself.crossing_time = crossing_timeself.wait_time = 0self.action = env.process(self.run())def run(self):yield self.env.timeout(self.arrival_time)print(f'行人 {self.name} 到达人行横道 at {self.env.now}')while self.traffic_light.state != 'red':yield self.env.timeout(1)self.wait_time += 1yield self.env.timeout(self.crossing_time)print(f'行人 {self.name} 通过人行横道 at {self.env.now} (等待时间: {self.wait_time})')# 主函数
def main():# 创建环境env = simpy.Environment()# 创建交通灯traffic_light = TrafficLight(env, green_duration=30, red_duration=20)# 生成车辆for i in range(10):arrival_time = random.expovariate(1.0 / 10)  # 平均每10秒一辆车travel_time = random.uniform(5, 10)  # 行驶时间在5到10秒之间vehicle = Vehicle(env, f'V{i}', traffic_light, arrival_time, travel_time)# 生成行人for i in range(5):arrival_time = random.expovariate(1.0 / 20)  # 平均每20秒一个行人crossing_time = random.uniform(5, 10)  # 过马路时间在5到10秒之间pedestrian = Pedestrian(env, f'P{i}', traffic_light, arrival_time, crossing_time)# 运行模拟env.run()if __name__ == "__main__":main()

其输出如下：

行人 P4 到达人行横道 at 0.0854511474780956
行人 P0 到达人行横道 at 0.22071825169249826
车辆 V8 到达路口 at 0.4483593469050312
车辆 V2 到达路口 at 1.9506493505845486
车辆 V4 到达路口 at 2.2711882878007974
车辆 V0 到达路口 at 2.3325381214341854
车辆 V6 到达路口 at 3.205591868640017
车辆 V3 到达路口 at 5.9670463790079475
车辆 V9 到达路口 at 6.241862753212637
车辆 V8 通过路口 at 8.950071244844713 (等待时间: 0)
车辆 V0 通过路口 at 9.605673370614314 (等待时间: 0)
车辆 V2 通过路口 at 11.011427474897479 (等待时间: 0)
车辆 V4 通过路口 at 12.039874222478037 (等待时间: 0)
车辆 V3 通过路口 at 12.379120050394778 (等待时间: 0)
车辆 V5 到达路口 at 12.815664910259631
车辆 V9 通过路口 at 12.909770077209995 (等待时间: 0)
车辆 V6 通过路口 at 12.940698199430786 (等待时间: 0)
车辆 V7 到达路口 at 13.436297816013289
车辆 V5 通过路口 at 18.593095516768905 (等待时间: 0)
车辆 V7 通过路口 at 20.464409648143516 (等待时间: 0)
车辆 V1 到达路口 at 22.612405519343664
车辆 V1 通过路口 at 32.402998439189474 (等待时间: 0)
行人 P4 通过人行横道 at 36.60432112172377 (等待时间: 30)
行人 P0 通过人行横道 at 38.00086713707282 (等待时间: 30)
行人 P2 到达人行横道 at 53.08467772803369
行人 P1 到达人行横道 at 58.35154343496454
行人 P1 通过人行横道 at 87.23696711920766 (等待时间: 22)
行人 P2 通过人行横道 at 87.98903365836632 (等待时间: 27)
行人 P3 到达人行横道 at 94.55501311405756
行人 P3 通过人行横道 at 104.20111058822664 (等待时间: 0)

最后，我们通过这个结果就可以分析得到平均等待时间、交通流量、系统效率等信息。然后可以通过这信息对于原本的交通管理系统进行进一步的优化。

四、补充

在这里，我简单介绍下simpy函数库，它是一个用于离散事件仿真的Python库。它提供了一种简单而强大的方式来模拟各种系统，包括但不限于生产系统、物流系统、通信网络、交通系统等。simpy 的核心思想是通过事件驱动的方式来模拟系统的动态行为，使得开发者可以专注于系统逻辑的建模，而不必过多关注底层的时间管理细节。

关于它的基本用法有5点，具体如下：

首先是导入，这不必多说。

然后是创建环境：

env = simpy.Environment()

其次是定义进程，进程是simpy中的基本单元，表示系统中的一个实体。进程通过定义一个生成器函数来实现：

def my_process(env):while True:print(f'Time: {env.now}, Process is running')yield env.timeout(1)  # 暂停1个时间单位

接下来是运行进程：

env.process(my_process(env))

最后是运行仿真，指定具体时间：

env.run(until=10)

此上