蒙特卡罗法模拟

wyy下午发给了我一道思考题蛮有意思的。赚了一杯奶茶的米。

问题:某个城市急救站,平均一小时会接到4个呼救电话(平均间隔15分钟)。一次救护行动平均需要35分钟(救护车开出到送完病人回来)。用蒙特卡罗法模拟其工作过程。急救站要配备多少辆救护车,才能保证有人呼救时,有90%的概率站上有车可出?(等待时间在5分钟以内可以接受,不算等)

q

蒙特卡罗法

蒙特卡罗方法又称统计模拟法、随机抽样技术,是一种随机模拟方法,以概率和统计理论方法为基础的一种计算方法,是使用随机数(或更常见的伪随机数)来解决很多计算问题的方法。将所求解的问题同一定的概率模型相联系,用电子计算机实现统计模拟或抽样,以获得问题的近似解。为象征性地表明这一方法的概率统计特征,故借用赌城蒙特卡罗命名。


将题目抽象化简

  1. 一个色子的1,2点视作有求助电话;3,4,5,6视作无求助电话。每5min判断一次。一小时内的期望为4次。

  2. 两个色子随机出每次出车的服务时间10min~60min

和数 2 3 4 5 6 7 8 9 10 11 12
概率 1/36 2/36 3/36 4/36 5/36 6/36 5/36 4/36 3/36 2/36 1/36

Python实现概率分布

1.

1
2
3
4
5
6
7
8
9
10
11
12
13
import random

# 服务时间概率相关变量
numbers = [10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60]
probabilities = [1/36, 2/36, 3/36, 4/36, 5/36, 6/36, 5/36, 4/36, 3/36, 2/36, 1/36]

# 依据probabilities输出每次的服务时间
def service_time():
random_number = random.random()

cumulative_probabilities = [sum(probabilities[:i+1]) for i in range(len(probabilities))]
index = next(i for i, p in enumerate(cumulative_probabilities) if p > random_number)
return numbers[index]

2.

1
2
3
4
5
6
7
# 1/3的概率有电话打来
def sendcar():
random_number = random.random()
if random_number < 1/3:
return True
else:
return False

具体算法

1
2
# listTime为每分钟有多少救护车的列表
listTime=[carsNumber]*timeLength

利用列表存储每个时刻可以使用车的数量。

当收到求助电话的时刻,等待,直到有车辆可供使用,记录下等待的时长与车辆即将服务的时长。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
# 在时间线上找
for waitTime in range(currentTime+5, timeLength):
if listTime[waitTime]>0 and waitTime-currentTime<10:# 打完电话(currentTime+5)后 等待时间小于5分钟
validResponse+=1
service_time_value = str(service_time())
beginTime = waitTime
finishTime = beginTime + eval(service_time_value)
print(f"有效发车 电话到达时间{currentTime} | 实际发车时间{beginTime} | 服务持续时间{service_time_value} | 服务完成时间{finishTime}", end='\n')

# 根据题目情景,打电话开始时刻5分钟后派遣出可用车辆
for j in range(waitTime , waitTime + eval(service_time_value)):
if j>=timeLength:
break
listTime[j] -= 1`#重点

currentTime+=1
break

完整代码

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
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
import random

# 服务时间概率相关变量
numbers = [10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60]
probabilities = [1/36, 2/36, 3/36, 4/36, 5/36, 6/36, 5/36, 4/36, 3/36, 2/36, 1/36]

# 依据probabilities输出每次的服务时间
def service_time():
random_number = random.random()

cumulative_probabilities = [sum(probabilities[:i+1]) for i in range(len(probabilities))]
index = next(i for i, p in enumerate(cumulative_probabilities) if p > random_number)
return numbers[index]

# 1/3的概率有电话打来
def sendcar():
random_number = random.random()
if random_number < 1/3:
return True
else:
return False
# 救护车数量和时间长度(分钟)
carsNumber=4
timeLength=24

# listTime为每分钟有多少救护车的列表
listTime=[carsNumber]*timeLength

currentTime = 0 #当前时间
validResponse=0 #有效发车次数
invalidResponse=0 #无效发车次数

# 主函数部分
while currentTime < timeLength:
# 在5t时刻
if(currentTime%5==0):
# 判断是否有电话打来
if sendcar():
# 在时间线上找
for waitTime in range(currentTime+5, timeLength):
if listTime[waitTime]>0 and waitTime-currentTime<10:# 打完电话(currentTime+5)后 等待时间小于5分钟
validResponse+=1
service_time_value = str(service_time())
beginTime = waitTime
finishTime = beginTime + eval(service_time_value)
print(f"有效发车 电话到达时间{currentTime} | 实际发车时间{beginTime} | 服务持续时间{service_time_value} | 服务完成时间{finishTime}", end='\n')

# 根据题目情景,打电话开始时刻5分钟后派遣出可用车辆
for j in range(waitTime , waitTime + eval(service_time_value)):
if j>=timeLength:
break
listTime[j] -= 1

currentTime+=1
break

elif listTime[waitTime]>0 and waitTime-currentTime>=10:# 打完电话(currentTime+5)后 等待时间大于5分钟
invalidResponse+=1
service_time_value = str(service_time())
beginTime = waitTime
finishTime = beginTime + eval(service_time_value)
print(f"无效发车 电话到达时间{currentTime} | 实际发车时间{beginTime} | 服务持续时间{service_time_value} | 服务完成时间{finishTime}", end='\n')

# 根据题目情景,打电话开始时刻5分钟后派遣出可用车辆
for j in range(waitTime , waitTime + eval(service_time_value)):
if j>=timeLength:
break
listTime[j] -= 1

currentTime+=1
break
# 在限定时间内找不到可用车辆
else:
print(f"找不到可用无效发车 电话到达时间{currentTime} ", end='\n')
invalidResponse+=1
currentTime+=1

# 在5t该时刻无电话打来
else:
currentTime+=1

# 不在5t时刻
else:
currentTime += 1

print("有效发车:",validResponse)
print("无效发车:",invalidResponse)

# 有效发车/无效发车+有效发车
rate = validResponse/(validResponse+invalidResponse)

print(f"当救护车辆为{carsNumber}时运行{timeLength}分钟,有效发车率为{rate:.3f}", end='\n')

运行结果

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
32
33
34
35
有效发车: 30288
无效发车: 4913
当救护车辆为4时运行525600分钟,有效发车率为0.860

有效发车: 30122
无效发车: 4900
当救护车辆为4时运行525600分钟,有效发车率为0.860

有效发车: 30278
无效发车: 4884
当救护车辆为4时运行525600分钟,有效发车率为0.861

有效发车: 30335
无效发车: 4591
当救护车辆为4时运行525600分钟,有效发车率为0.869

有效发车: 30173
无效发车: 4620
当救护车辆为4时运行525600分钟,有效发车率为0.867

有效发车: 34040
无效发车: 976
当救护车辆为5时运行525600分钟,有效发车率为0.972

有效发车: 34044
无效发车: 1002
当救护车辆为5时运行525600分钟,有效发车率为0.971

有效发车: 33988
无效发车: 939
当救护车辆为5时运行525600分钟,有效发车率为0.973

有效发车: 34129
无效发车: 1035
当救护车辆为5时运行525600分钟,有效发车率为0.971