马尔可夫排队网络——Python分析

马尔科夫排队网络（Markovian Queueing Networks）是一类特殊的排队网络，假设系统中的到达过程和服务时间均遵循指数分布，系统状态之间的转移遵循马尔可夫性质。这些假设使得马尔科夫排队网络可以通过解析方法进行分析，从而为实际系统的设计和性能优化提供理论依据。通过理论推导和模型构建，马尔科夫排队网络可以定量分析系统性能，识别瓶颈，并提出改进方案。尽管实际系统可能具有更复杂的特性（如非指数分布的到达和服务时间、优先级队列等），但马尔科夫模型提供了一个重要的基础和起点，能够帮助我们理解和解决复杂的排队问题。在实际应用中，马尔科夫排队网络被广泛用于电信网络、计算机系统、制造系统、交通工程和医疗服务等领域。例如，在呼叫中心中，客户呼叫到达和服务时间通常被建模为泊松过程和指数分布。通过使用马尔科夫排队模型，可以确定平均等待时间、座席利用率和放弃率，从而优化客户服务和资源配置。

一、马尔科夫排队网络模型

如果一个排队系统只是一个节点，那么多个节点就可以组成一个排队网络。每个节点都含有一个服务机构和排队机构，是一个简单的排队系统，当顾客离开某个排队系统节点就进入下一个节点。例如数据包在通信网中传输。马尔可夫排队网络指的是各节点外部到达顾客流是泊松流，各节点的服务时间是负指数分布的网络系统。

1.1 马尔科夫排队系统分类特点

根据系统的结构和客户流动的方式，马尔科夫排队系统可分为两类：开马尔科夫排队网络和闭马尔科夫排队网络。

开马尔科夫排队网络	闭马尔科夫排队网络
外部到达率：客户可以从外部进入系统，通常到达过程服从泊松分布	固定客户数量：系统中客户的总数是固定的，这些客户只能在网络内部的各个节点之间转移
内部转移：客户在不同服务节点之间转移，转移过程符合马尔科夫性质（记忆无关性）	无外部到达：没有客户从外部进入系统
离开系统：客户在完成服务后可以离开系统，即有明确的离开机制	无离开机制：客户不能离开系统，只能在内部循环
不固定的客户数量：系统中的客户数量可以随时间变化，没有固定的总数量限制	内部转移：客户在节点之间的转移同样符合马尔科夫性质
应用场景：呼叫中心：客户电话呼入呼出；计算机网络：数据包的到达和转发	应用场景：计算机系统中的进程调度：固定数量的任务在不同处理器之间调度。

1.2 开马尔科夫排队网络模型

• 假设
  外部到达过程：客户到达系统的过程是泊松过程，且到达率为$\lambda_i$
  服务时间分布：服务时间服从指数分布，服务速率为 $\mu_i$
  转移过程：客户在不同节点之间的转移概率矩阵为 $P = [p_{ij}]$ ，其中$p_{ij}$ 表示客户从节点$i$ 转移到节点$j$ 的概率
  状态转移：系统状态的转移服从马尔可夫性质（记忆无关性）。

• 系统绩效公式

流平衡方程：对于每个节点$i$ ，总到达率$\lambda_i$ 包括外部到达率和其他节点转发到该节点的客户，即：

$$\lambda_i = \lambda_i^{\text{ext}} + \sum_{j \neq i} \lambda_j p_{ji}$$

平均队列长度$L_i$：使用 M/M/1 排队模型公式：

$$L_i = \frac{\lambda_i}{\mu_i - \lambda_i}$$

平均停留时间 ( W )：

$$W = \sum_{i} \frac{\lambda_i L_i}{\sum_{i} \lambda_i}$$

平均吞吐量：系统在稳态下的平均吞吐量等于系统的总到达率：

$$\lambda_{\text{total}} = \sum_{i} \lambda_i$$

二、开马尔科夫排队网路实例

一个包含三个节点的马尔科夫排队网络，其中：

• 到达率：$\lambda_1 = \lambda_2 = 15$分组/秒，$\lambda_3 = 20$分组/秒。
• 各节点之间的转发矩阵 $[r_{ij}]$ 如最上面图所示。
• 每个节点的服务速率 $U = 50$分组/秒。

计算：网络中的平均分组数；每一分组在网络中的总平均停留时间；平均吞吐量。

• 理论分析

首先，根据转发矩阵$[r_{ij}]$，我们可以求出每个节点的总输入率$\lambda_i$

• 转发矩阵$[r_{ij}]$：

$$r_{ij} = \begin{pmatrix} 0 & \frac{1}{3} & \frac{1}{3} \\ \frac{1}{6} & 0 & \frac{1}{3} \\ \frac{1}{4} & \frac{1}{8} & 0 \end{pmatrix}$$

• 流平衡方程

对于一个包含多个节点的马尔科夫排队网络，每个节点的流平衡方程可以表示为进入节点的流量等于离开节点的流量。

• 节点1的流平衡方程：
  进入节点1的流量包括外部到达率 $\lambda_1$ 和其他节点转移到节点1的流量：

$$\lambda_1 + \pi_2 \cdot r_{21} + \pi_3 \cdot r_{31} = \pi_1 \cdot \left(1 - r_{11}\right)$$

因为 $r_{11}$= 0

$$\lambda_1 + \pi_2 \cdot r_{21} + \pi_3 \cdot r_{31} = \pi_1$$

• 节点2的流平衡方程：

进入节点2的流量包括外部到达率 λ2\lambda_2λ2​ 和其他节点转移到节点2的流量：

$$\lambda_2 + \pi_1 \cdot r_{12} + \pi_3 \cdot r_{32} = \pi_2 \cdot \left(1 - r_{22}\right)$$

因为$r_{22} = 0$

$$\lambda_2 + \pi_1 \cdot r_{12} + \pi_3 \cdot r_{32} = \pi_2$$

• 节点3的流平衡方程：

进入节点3的流量包括外部到达率 $\lambda_3$和其他节点转移到节点3的流量：

$$\lambda_3 + \pi_1 \cdot r_{13} + \pi_2 \cdot r_{23} = \pi_3 \cdot \left(1 - r_{33}\right)$$

因为$r_{33} = 0$

$$\lambda_3 + \pi_1 \cdot r_{13} + \pi_2 \cdot r_{23} = \pi_3$$

• 构建线性方程组

将上面的平衡方程表示为矩阵形式：

$$\begin{cases} \lambda_1 + \pi_2 \cdot r_{21} + \pi_3 \cdot r_{31} = \pi_1 \\ \lambda_2 + \pi\_1 \cdot r_{12} + \pi_3 \cdot r_{32} = \pi_2 \\ \lambda_3 + \pi_1 \cdot r_{13} + \pi_2 \cdot r_{23} = \pi_3 \end{cases}$$

这可以写成矩阵形式 $A \cdot \pi=b$

$$\begin{pmatrix} 1 & -r_{21} & -r_{31} \\ -r_{12} & 1 & -r_{32} \\ -r_{13} & -r_{23} & 1 \end{pmatrix} \begin{pmatrix} \pi_1 \\ \pi_2 \\ \pi_3 \end{pmatrix}= \begin{pmatrix} \lambda_1 \\ \lambda_2 \\ \lambda_3 \end{pmatrix}$$

• 性能指标

平均分组数$L_i$：利用 M/M/1 排队模型公式

$$L_i = \frac{\lambda_i}{\mu_i - \lambda_i}$$

其中 $\mu_i = U$。

总平均停留时间$W$：

$$W = \sum_{i=1}^{3} \frac{\lambda_i L_i}{\sum_{i=1}^{3} \lambda_i}$$

平均吞吐量：由于系统达到稳定状态，平均吞吐量等于到达系统的总分组率，即

$$\lambda = \sum_{i=1}^{3} \lambda_i$$

import numpy as np

# 输入已知参数
lambda_vals = np.array([15, 15, 20])
U = 50
r = np.array([
    [0, 1/3, 1/3],
    [1/6, 0, 1/3],
    [1/4, 1/8, 0]
])

# 节点数
n = len(lambda_vals)

# 构建系数矩阵A
A = np.zeros((n, n))
for i in range(n):
    for j in range(n):
        if i != j:
            A[i, j] = -r[j, i]  # 系数为转移概率的负值
            A[i, i] += r[j, i]  # 对角线元素为到达该节点的概率和

# 输出系数矩阵A进行检查
print("系数矩阵A:")
print(A)

# 增加对角线元素以包括服务速率U
for i in range(n):
    A[i, i] += 1  # 对角线增加1以包括自身服务速率影响

# 输出调整后的系数矩阵A进行检查
print("调整后的系数矩阵A:")
print(A)

# 构建向量b
b = lambda_vals

# 输出向量b进行检查
print("向量b:")
print(b)

# 求解线性方程组
pi = np.linalg.solve(A, b)

# 计算网络中的平均分组数
L = np.sum(pi)

# 每一分组在网络中的总平均停留时间
W = L / np.sum(lambda_vals)

# 平均吞吐量
X = np.sum(lambda_vals)

# 输出结果
pi_rounded = np.round(pi, 2)
L_rounded = round(L, 2)
W_rounded = round(W, 2)
X_rounded = round(X, 2)

print("节点平均分组数 pi:")
print(pi_rounded)
print("网络中的平均分组数 L:")
print(L_rounded)
print("每一分组在网络中的总平均停留时间 W:")
print(W_rounded)
print("平均吞吐量 X:")
print(X_rounded)

节点平均分组数 pi:
[15.61 15.42 18.21]
网络中的平均分组数 L:
49.24
每一分组在网络中的总平均停留时间 W:
0.98
平均吞吐量 X:
50

三、数据中心可靠性分析

数据中心	马尔科夫链

问题：

自动诊断系统使用“ping”来检测故障，并在检测到故障时通知工作人员。
机器有两种状态：运行中或故障。
故障的机器：等待$N$个维修人员中的一个进行维修或正在维修过程中。
如果下一次故障的时间和维修时间都呈指数分布，我们可以使用马尔可夫链来解决这个排队问题！

性能问题:

在任何给定时间，恰好有 $\mathrm{j}(1 \leq \mathrm{j} \leq \mathrm{M})$ 台机器在运行的概率是多少?
在任何给定时间，至少有 $\mathrm{j}(1 \leq \mathrm{j} \leq \mathrm{M})$ 台机器在运行的概率是多少?
为了保证至少有 $j(1 \leq \mathrm{j} \leq \mathrm{M})$ 台机器以给定的概率在运行，需要多少维修人员?
维修团队的规模对修复机器的平均时间 (MTTR) 有什么影响? 维修团队的规模对可以预期在任何给定时间运行的机器的百分比有什么影响?
维修团队成员修复机器所需的平均时间（即，他们的技能水平) 对总MTTR（即，包括等待维修的时间）有什么影响? 维修人员的技能水平如何影响运行中的机器的百分比?

3.1 模型构建

假设：机器独立失败，并且以相同的速率失败，所有维修人员具有相同的技能水平，因此维修速率也相同。用$\lambda$表示机器的故障率。因此，$\frac{1}{\lambda}$是平均故障时间。用µ表示维修速率。
状态：令状态$k$表示故障的机器数量（故障的机器需要维修）。状态转移：当一台机器失败时，从状态$k$转换到$k + 1$;当一台机器被修复时，从状态$k$转换到$k-1$；在状态$k$时，有$M-k$台机器在运行，每台的故障率为$\lambda$

• 在状态$k$的总故障率$λ_k$，是：

  $$\lambda_k = (M - k) \lambda \quad 其中k = 0, ..., M - 1$$

• 在状态$k$的总维修率$μ_k$，取决于所有$N$个维修人员是否都忙：

$$\mu_k = \begin{cases} k\mu & \text{if } k = 1, ..., N \\ N \mu & \text{if } k = N + 1, ..., M \end{cases}$$

• 求解模型（马尔可夫链）以找到$P_k$，即处于状态$k$的稳态概率$（k=0，...，M）$，也就是说$k$台机器故障。

$$p_k= \begin{cases}p_0\left(\frac{\lambda}{\mu}\right)^k\binom{M}{k} & k=1, \ldots, N \\ p_0\left(\frac{\lambda}{\mu}\right)^k\binom{M}{k} \frac{N^{N-k} k!}{N!} & k=(N+1), \ldots, M\end{cases}$$

其中 $p_0$ 可以通过满足 $\sum_{k=0}^M p_k=1$ 获得。因此，

$$p_0=\left[\sum_{k=0}^N\left(\frac{\lambda}{\mu}\right)^k\binom{M}{k}+\sum_{k=N+1}^M\left(\frac{\lambda}{\mu}\right)^k\binom{M}{k} \frac{N^{N-k} k!}{N!}\right]^{-1}$$

• 平均机器故障率（平均故障率）：

$$\bar{X}_f=\sum_{k=0}^{M-1} \lambda_k \times p_k=\sum_{k=0}^{M-1}(M-k) \lambda p_k$$

• 维修平均时间：

$$\mathrm{MTTR}=\frac{M}{\bar{X}_f}-\mathrm{MTTF}=\frac{M}{\bar{X}_f}-1 / \lambda$$

• 平均故障机器数:

$$N_f=\bar{X}_f \times \mathrm{MTTR}=M-\bar{X}_f \times \mathrm{MTTF}=M-\bar{X}_f / \lambda$$

• 平均运行机器数:

$$N_o=M-N_f=\bar{X}_f \times \mathrm{MTTF}=\bar{X}_f / \lambda$$

3.2 $N = 5, 15, 25$时故障状态$j$的概率演化

import math
from scipy.special import comb
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

def calculate_p0(M, N, lambda_, mu):
    sum_1 = sum((lambda_/mu)**k * comb(M, k) for k in range(N+1))
    sum_2 = sum((lambda_/mu)**k * comb(M, k) * (N**(N-k) * math.factorial(k) / math.factorial(N)) for k in range(N+1, M+1))
    p0 = 1 / (sum_1 + sum_2)
    return p0

def calculate_pk(M, N, lambda_, mu, p0):
    pks = []
    for k in range(M+1):
        if k <= N:
            pk = p0 * (lambda_/mu)**k * comb(M, k)
        else:
            pk = p0 * (lambda_/mu)**k * comb(M, k) * (N**(N-k) * math.factorial(k) / math.factorial(N))
        pks.append(pk)
    return pks

# Parameters
M = 120
lambda_ = 1  # Assuming λ = 1 for calculation
mu = 2  # Assuming μ = 1 for calculation
Ns = [5, 15, 25]

# Calculations for each N
probabilities = {}
for N in Ns:
    p0 = calculate_p0(M, N, lambda_, mu)
    pks = calculate_pk(M, N, lambda_, mu, p0)
    probabilities[N] = pks

# Setting font sizes
plt.rcParams.update({'font.size': 24})  # Globally setting the font size

# Plotting the probabilities for each state
plt.figure(figsize=(12, 8))
for N in Ns:
    plt.plot(range(M+1), probabilities[N], label=f'N={N}', linewidth=3)

plt.xlabel('State (j)', fontsize=24)
plt.ylabel('Probability', fontsize=24)
plt.title('State Probabilities for Different Values of N', fontsize=24)
plt.legend(fontsize=18)
plt.grid(True)
plt.show()

3.3 系统稳态分析

这里选择$\lambda=1$和$\mu=2$，系统稳态结果见下表：

$N$	5	10	15	20	25	30	35	40	45	50
平均机器故障率 (X_f)	10.00	20.00	30.00	40.00	50.00	60.00	69.83	77.03	79.51	79.95
维修平均时间(MTTR)	11.00	5.00	3.00	2.00	1.40	1.00	0.72	0.56	0.51	0.50
平均故障机器数 (N_f)	110.00	100.00	90.00	80.00	70.00	60.00	50.17	42.97	40.49	40.05
平均运行机器数 (N_o)	10.00	20.00	30.00	40.00	50.00	60.00	69.83	77.03	79.51	79.95
运行机器率 (N_o/M * 100)	8.33	16.67	25.00	33.33	41.67	50.00	58.19	64.19	66.25	66.63

import math
from scipy.special import comb
import pandas as pd

def calculate_p0(M, N, lambda_, mu):
    sum_1 = sum((lambda_/mu)**k * comb(M, k) for k in range(N+1))
    sum_2 = sum((lambda_/mu)**k * comb(M, k) * (N**(N-k) * math.factorial(k) / math.factorial(N)) for k in range(N+1, M+1))
    p0 = 1 / (sum_1 + sum_2)
    return p0

def calculate_pk(M, N, lambda_, mu, p0):
    pks = []
    for k in range(M+1):
        if k <= N:
            pk = p0 * (lambda_/mu)**k * comb(M, k)
        else:
            pk = p0 * (lambda_/mu)**k * comb(M, k) * (N**(N-k) * math.factorial(k) / math.factorial(N))
        pks.append(pk)
    return pks

def calculate_metrics_with_ratio(M, N, lambda_, mu):
    p0 = calculate_p0(M, N, lambda_, mu)
    pks = calculate_pk(M, N, lambda_, mu, p0)
    X_f = sum((M-k) * lambda_ * pks[k] for k in range(M))
    MTTR = M / X_f - 1 / lambda_
    N_f = M - X_f / lambda_
    N_o = X_f / lambda_
    N_o_ratio = (N_o / M) * 100
    return X_f, MTTR, N_f, N_o, N_o_ratio

# Parameters
M = 120
lambda_ = 1  # Assuming λ = 1 for calculation
# Increasing mu to improve repair efficiency
mu = 2  # New μ value
Ns = range(5, 51, 5)  # Same range for N

# Calculations with the adjusted Ns range and the additional metric
results_with_ratio_adjusted = []
for N in Ns:
    X_f, MTTR, N_f, N_o, N_o_ratio = calculate_metrics_with_ratio(M, N, lambda_, mu)
    results_with_ratio_adjusted.append((N, X_f, MTTR, N_f, N_o, N_o_ratio))

# Output results with the adjusted Ns range and the additional metric
df_with_ratio_adjusted = pd.DataFrame(results_with_ratio_adjusted, columns=['N', 'Average Failure Rate (X_f)', 'MTTR', 'Average Failed Machines (N_f)', 'Average Operating Machines (N_o)', 'Operating Machines Ratio (N_o/M * 100)'])
df_transposed_with_ratio_adjusted = df_with_ratio_adjusted.set_index('N').transpose()
print(df_transposed_with_ratio_adjusted)

为了增加$N_o / M * 100$，可以考虑以下几种方法：
增加维修人员的数量 (N)：增加维修人员的数量会减少平均修理时间 (MTTR)，从而增加机器的平均运行数量 (N_o)。随着维修人员数量的增加，更多的机器能够得到及时维修并恢复运行。
提高维修效率 (增加 μ)：提高每个维修人员的维修效率，即增加维修率𝜇，也会减少 MTTR，从而增加 N_o。
减少故障率 (减少 λ)：降低机器的故障率 λ 会减少故障机器的数量，从而增加平均运行机器的数量 (N_o)。

总结

马尔科夫排队网络模型用于描述和分析系统中客户排队和服务行为，分为开马尔科夫排队网络和闭马尔科夫排队网络。开网络允许客户从外部进入和离开，客户数量不固定，适用于呼叫中心和计算机网络等场景。闭网络中客户数量固定，无外部到达和离开，适用于制造系统和进程调度等场景。通过平衡方程和转移概率矩阵分析，马尔科夫排队网络能够量化系统性能，识别瓶颈，提供优化方案。尽管假设简单，它们为解决复杂排队问题提供了重要基础。

参考文献

1. 马尔可夫排队模型
2. 马尔可夫排队网络
3. 第四章 马尔科夫排队系统