运筹学练习Python精解——排队论

练习1

自动取款机问题：银行计划安置取款机，A 机价格和平均服务率都是 B 机的 2 倍，应购置 1 台 A 机还是 2 台 B 机？顾客平均每分钟到达 1 位，A 型机的平均服务时间为 0.9，B 型机为 1.8 分钟，顾客到达间隔和服务时间都服从指数分布。

模型参数

• $\lambda$: 顾客到达率（顾客每分钟到达 1 位）
• $\mu_A$: A 型机的服务率（每分钟服务顾客数），$\mu_A = \frac{1}{0.9} = 1.111$
• $\mu_B$: B 型机的服务率（每分钟服务顾客数），$\mu_B = \frac{1}{1.8} = 0.556$

对于 1 台 A 型机，使用 M/M/1 模型：

$$\rho_A = \frac{\lambda}{\mu_A}$$

对于 2 台 B 型机，使用 M/M/2 模型：

$$\rho_B = \frac{\lambda}{2 \mu_B}$$

需要计算的指标
系统中平均顾客数 $L$；系统中平均等待时间$W$；队列中平均顾客数$L_q$；队列中平均等待时间 $W_q$

def mm1_metrics(arrival_rate, service_rate):分钟
    rho = arrival_rate / service_rate
    L = round(rho / (1 - rho), 2)
    W = round(L / arrival_rate, 2)
    L_q = round(L - rho, 2)
    W_q = round(L_q / arrival_rate, 2)
    return L, W, L_q, W_q

def mm2_metrics(arrival_rate, service_rate):
    rho = arrival_rate / (2 * service_rate)
    if rho >= 1:
        raise ValueError("System is unstable (rho >= 1)")
    L_q = round((arrival_rate**2 * (1 + (arrival_rate / (2 * service_rate)))) / 
                (2 * service_rate * (2 * service_rate - arrival_rate)), 2)
    L = round(L_q + (arrival_rate / service_rate), 2)
    W = round(L / arrival_rate, 2)
    W_q = round(L_q / arrival_rate, 2)
    return L, W, L_q, W_q

# 参数设置
arrival_rate = 1  # 顾客到达率（每分钟 1 位）
service_rate_A = 1 / 0.9  # A 型机服务率（每分钟服务 1.111 位顾客）
service_rate_B = 1 / 1.8  # B 型机服务率（每分钟服务 0.556 位顾客）

# 计算 1 台 A 型机的指标
metrics_A = mm1_metrics(arrival_rate, service_rate_A)
print("1 台 A 型机的指标 (L, W, L_q, W_q):", metrics_A)

# 计算 2 台 B 型机的指标
metrics_B = mm2_metrics(arrival_rate, service_rate_B)
print("2 台 B 型机的指标 (L, W, L_q, W_q):", metrics_B)
#1 台 A 型机的指标 (L, W, L_q, W_q): (9.0, 9.0, 8.1, 8.1)
#2 台 B 型机的指标 (L, W, L_q, W_q): (17.19, 17.19, 15.39, 15.39)

练习2

银行经理正试图通过提供更好的服务来提高顾客满意度。管理层期待顾客平均等待时间小于2分钟，平均队列长度（等待队列的长度）是2人或者更少。银行估算，每天大约为150名顾客提供服务。以下是现有的到达和服务时间的数据，试对该系统进行仿真分析。

到达时间表

到达间隔（min.）	0	1	2	3	4	5
概率	0.1	0.15	0.1	0.35	0.25	0.05

服务时间表

服务时间（min.）	1	2	3	4
概率	0.25	0.2	0.4	0.15

import numpy as np
import matplotlib.pyplot as plt

# 参数设置
n = 150
ta = np.array([5, 4, 3, 2, 1, 0])
pa = np.array([0.05, 0.25, 0.35, 0.10, 0.15, 0.10])
ts = np.array([4, 3, 2, 1])
ps = np.array([0.15, 0.40, 0.20, 0.25])

# 累积概率
pacum = np.cumsum(pa)
pscum = np.cumsum(ps)

# 计算到达时刻
Tarrival = np.random.rand(n)
for i in range(len(pa)):
    Tarrival[Tarrival < pacum[i]] = ta[i]
Tarrival = np.cumsum(Tarrival)

# 计算服务时间
Tservice = np.random.rand(n)
for i in range(len(ps)):
    Tservice[Tservice < pscum[i]] = ts[i]

# 初始化变量
Tstart = np.zeros(n)
Tleave = np.zeros(n)
Twait = np.zeros(n)
line = np.zeros(n)

# 初始化第一位顾客
Tstart[0] = Tarrival[0]
Tleave[0] = Tstart[0] + Tservice[0]
Twait[0] = Tleave[0] - Tarrival[0] - Tservice[0]
line[0] = 0

# 计算每位顾客的开始服务时刻、离开时刻、等待时间和队长
for i in range(1, n):
    Tstart[i] = max(Tleave[i-1], Tarrival[i])
    Tleave[i] = Tstart[i] + Tservice[i]
    Twait[i] = Tleave[i] - Tarrival[i] - Tservice[i]
    
    k = i - 1
    while k >= 0 and Tarrival[i] < Tleave[k]:
        line[i] += 1
        k -= 1

# 计算排队系统指标
average_waiting_time = np.mean(Twait)
average_queue_length = np.mean(line)
average_service_time = np.mean(Tservice)
average_arrival_rate = 1 / np.mean(np.diff(np.insert(Tarrival, 0, 0)))

# 输出排队系统指标
print(f"平均等待时间 (Wq): {average_waiting_time:.2f} 分钟")
print(f"平均队列长度 (Lq): {average_queue_length:.2f} 人")
print(f"平均服务时间 (Ts): {average_service_time:.2f} 分钟")
print(f"平均到达率 (λ): {average_arrival_rate:.2f} 每分钟")

# 绘制直方图
plt.figure(figsize=(12, 6))

plt.subplot(1, 2, 1)
plt.hist(Twait, bins=20, edgecolor='black')
plt.title('Histogram of Waiting Time')
plt.xlabel('Waiting Time')
plt.ylabel('Frequency')

plt.subplot(1, 2, 2)
plt.hist(line, bins=20, edgecolor='black')
plt.title('Histogram of Queue Length')
plt.xlabel('Queue Length')
plt.ylabel('Frequency')

plt.tight_layout()
plt.show()

平均等待时间 (Wq): 2.83 分钟
平均队列长度 (Lq): 1.33 人
平均服务时间 (Ts): 2.53 分钟
平均到达率 (λ): 0.36 每分钟

根据上面指标判断该系统没有达到要求，可适当增加服务台或者其他资源。

练习3

某修理厂为设备检修服务。已知检验的设备（顾客）到达服从Poisson分布，每天到达率$\lambda=42$台，当需要等待时，每台造成的损失为400元。服务(检修)时间服从负指数分布,平均每天服务率$\mu=20$台。每设置一个检修人员每天的服务成本为160元。问设立几个检修人员才能使平均总费用最小？

这个问题是一个经典的排队论优化问题。我们需要确定最优的修理人员数量，以使得总成本最小。总成本包括等待的损失和服务的成本。这里采用 M/M/c 模型 来求解。

M/M/c 模型的参数和公式

• 到达率$\lambda = 42$ 台/天
• 每个检修人员的服务率 $\mu = 20$ 台/天
• 每个检修人员每天的服务成本 = 160 元
• 每台设备等待的损失 = 400 元

系统绩效公式

• 系统利用率：$\rho = \frac{\lambda}{c\mu}$
• 计算概率 $P_0$： $$P_0 = \left[ \sum_{n=0}^{c-1} \frac{(\lambda/\mu)^n}{n!} + \frac{(\lambda/\mu)^c}{c!} \cdot \frac{1}{1 - \rho} \right]^{-1}$$
• 平均队列长度 $L_q$： $$L_q = \frac{P_0 (\lambda/\mu)^c \rho}{c! (1 - \rho)^2}$$
• 平均等待时间 $W_q$： $$W_q = \frac{L_q}{\lambda}$$
• 平均总费用： $$\text{总费用} = c \times 160 + \lambda \times W_q \times 400$$

import numpy as np
from scipy.special import factorial

# 参数
arrival_rate = 42  # λ，每天到达率
service_rate = 20  # μ，每天服务率
cost_per_wait = 400  # 每台设备等待的损失
cost_per_server = 160  # 每个检修人员每天的服务成本

def calculate_total_cost(c, arrival_rate, service_rate, cost_per_wait, cost_per_server):
    rho = arrival_rate / (c * service_rate)
    if rho >= 1:
        return float('inf')  # 系统不稳定，返回无穷大成本
    P0_inv_sum = sum((arrival_rate/service_rate)**n / factorial(n) for n in range(c))
    P0_inv_sum += (arrival_rate/service_rate)**c / (factorial(c) * (1 - rho))
    P0 = 1 / P0_inv_sum
    Lq = (P0 * (arrival_rate/service_rate)**c * rho) / (factorial(c) * (1 - rho)**2)
    Wq = Lq / arrival_rate
    total_cost = c * cost_per_server + arrival_rate * Wq * cost_per_wait
    return total_cost

# 找到最优的检修人员数量
min_cost = float('inf')
optimal_c = 1
for c in range(1, 101):  # 假设检修人员数量不会超过100
    total_cost = calculate_total_cost(c, arrival_rate, service_rate, cost_per_wait, cost_per_server)
    if total_cost < min_cost:
        min_cost = total_cost
        optimal_c = c

print(f"最优的检修人员数量: {optimal_c}")
print(f"最小平均总费用: {min_cost:.2f} 元")

最优的检修人员数量: 4
最小平均总费用: 728.17 元

练习4

考虑一个收款台的排队系统。某商店只有一个收款台，顾客到达收款台的时间间隔服从平均时间为10秒钟的负指数分布，负指数分布为：

$$f\left( x \right) =\begin{cases} \frac{1}{\lambda}e^{-\frac{x}{\lambda}}\,\,x>0\\ 0 x\leqslant 0\\ \end{cases}$$

每个顾客的服务时间服从均值为6.5秒，标准差为1.2秒的正态分布。利用计算机仿真计算顾客在收款台的平均逗留时间，系统的服务强度。

import numpy as np
import matplotlib.pyplot as plt

# 设置中文字体（以微软雅黑为例，根据系统中已安装的中文字体进行调整）
plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']
plt.rcParams['axes.unicode_minus'] = False  # 处理负号显示问题

# 模拟参数
lambda_arrival = 1 / 10  # 平均到达间隔时间为10秒，即到达率 λ = 1/10
mean_service = 6.5  # 平均服务时间为6.5秒
std_service = 1.2  # 服务时间的标准差为1.2秒

# 模拟时长（秒）
sim_time = 3600  # 模拟1小时，即3600秒

# 初始化仿真时钟和事件列表
clock = 0  # 初始化时钟
arrival_times = []  # 记录顾客到达时间
service_times = []  # 记录顾客开始服务时间
departure_times = []  # 记录顾客离开时间

while clock < sim_time:
    # 模拟顾客到达时间间隔
    inter_arrival_time = np.random.exponential(1 / lambda_arrival)
    clock += inter_arrival_time
    
    # 记录顾客到达时间
    arrival_times.append(clock)
    
    # 模拟顾客服务时间
    service_time = np.random.normal(mean_service, std_service)
    
    # 如果顾客到达时收款台空闲，则该顾客开始服务的时间即为到达时间
    if len(departure_times) == 0 or arrival_times[-1] > departure_times[-1]:
        service_start_time = arrival_times[-1]
    else:
        service_start_time = departure_times[-1]
    
    # 记录顾客开始服务的时间
    service_times.append(service_start_time)
    
    # 计算顾客离开的时间
    departure_time = service_start_time + service_time
    departure_times.append(departure_time)

# 计算顾客在系统中的平均逗留时间
waiting_times = np.array(service_times) - np.array(arrival_times)
mean_waiting_time = np.mean(waiting_times)

# 计算系统的服务强度（ρ）
arrival_rate = len(arrival_times) / sim_time
service_rate = 1 / mean_service
rho = arrival_rate / service_rate

print(f"顾客在收款台的平均逗留时间: {mean_waiting_time:.2f} 秒")
print(f"系统的服务强度 (ρ): {rho:.2f}")

# 绘制等待时间的直方图
plt.hist(waiting_times, bins=30, density=True, alpha=0.75)
plt.xlabel('等待时间（秒）')
plt.ylabel('概率密度')
plt.title('顾客在收款台的等待时间分布')
plt.grid(True)
plt.show()

练习5

某单位电话交换台有一台200门内线的总机，已知在上班8h的时间内，有20％的内线分机平均每40min要一次外线电话，80％的分机平均隔120min要一次外线。又知外线打入内线的电话平均每分钟1次。假设与外线通话的时间平均为3min，并且上述时间均服从负指数分布，如果要求电话的通话率为95％，问该交换台应设置多少条外线？

• 内线打外线，其服务强度为：

$$\lambda_{1} = \left( \frac{60}{40} \times 0.2 + \frac{60}{120} \times 0.8 \right) \times 200 = 140$$

• 外线打内线，其服务强度为：

$$\lambda_{2} = 1 \times 60 = 60$$

• 总服务强度

$$\lambda = \lambda_{1} + \lambda_{2} = 140 + 60 = 200$$

• 系统损失率
  这是一个损失制服务系统，按题目要求，系统损失率不能超过5%，即：

$$P_{\text{lost }} \leq 0.05$$

• 外线数量
  外线是整数，在满足条件下，条数越小越好。

import math

def calculate_rho(lamda, mu, s):
    return lamda / (mu * s)

def calculate_p_lost(rho, s):
    if rho >= 1:
        return 1  # 如果 rho >= 1，丢失概率为1
    else:
        sum_part = sum((rho ** k) / math.factorial(k) for k in range(s))
        p_lost = (rho ** s / math.factorial(s)) / (sum_part + (rho ** s / math.factorial(s)))
        return p_lost

def find_min_s(lamda, mu, target_p_lost):
    s = 1
    while True:
        rho = calculate_rho(lamda, mu, s)
        p_lost = calculate_p_lost(rho, s)
        if p_lost < target_p_lost:
            return s
        s += 1

# 模型参数
lamda = 200  # 到达率
mu = 60 / 3  # 服务率，这里假设单位是1/分钟

# 目标丢失概率
target_p_lost = 0.05

# 查找满足丢失概率要求的最小服务器数量 s
min_s = find_min_s(lamda, mu, target_p_lost)

print(f"满足丢失概率小于 {target_p_lost} 的最小服务器数量为: {min_s}")
#满足丢失概率小于 0.05 的最小服务器数量为: 11

练习6

某公司医务室为职工检查身体，职工的到达服从泊松分布，每小时平均到达50人，若职工不能按时体检，造成的损失为每小时每人平均60元。体检所花时间服从指数分布负指数分布，平均每小时服务率为$\mu$，每人的体检费为30元，试确定使公司总支出最少的参数$\mu$。

成本模型设计：
- 到达率 $\lambda$：每小时平均到达50人，因此 $\lambda = 50$
- 服务率 $\mu$：需要找到使总成本最小化的最佳服务率。
- 损失成本：如果职工不能按时体检，造成的损失为每小时每人平均60元。
- 体检费用：每人的体检费为30元。
- 总支出：由等待损失成本和体检费用组成。
计算总支出： 总支出可以表示为等待损失成本和体检费用的总和。等待损失成本是因为职工不能按时体检而导致的损失，计算方式为未服务人数乘以等待时间乘以损失率，而体检费用为服务人数乘以体检费用。
目标： 找到使总支出最少的服务率 $\mu$

import numpy as np

def calculate_rho(lamda, mu):
    return lamda / mu

def calculate_p0(lamda, mu):
    rho = calculate_rho(lamda, mu)
    return 1 / (1 + rho)

def calculate_average_customers(lamda, mu):
    rho = calculate_rho(lamda, mu)
    if rho >= 1:
        return float('inf')  # 如果 rho >= 1，系统不稳定，返回无穷大
    p0 = calculate_p0(lamda, mu)
    L = rho * p0 / (1 - rho)
    return L

def total_cost(lamda, loss_cost, checkup_cost, mu):
    L = calculate_average_customers(lamda, mu)
    Cs = checkup_cost
    Cw = loss_cost
    total_cost = Cs * mu + Cw * L
    return total_cost

# 模型参数
arrival_rate = 50  # 每小时到达率
loss_cost = 60  # 每小时每人的等待损失成本
checkup_cost = 30  # 每人的体检费用

# 寻找使总成本最小的服务率
best_cost = float('inf')
best_mu = None
for mu in np.linspace(0.1, 100, 1000):  # 从0.1到100的范围内尝试不同的服务率
    cost = total_cost(arrival_rate, loss_cost, checkup_cost, mu)
    if cost < best_cost:
        best_cost = cost
        best_mu = mu

print(f"使公司总支出最少的服务率 μ 为: {best_mu:.2f}")
print(f"最小总支出为: {best_cost:.2f}")

#使公司总支出最少的服务率 μ 为: 57.10，最小总支出为: 1938.27