排队论时延分析

实验 排队论时延分析

一，预习内容

实验概述

本次实验利用ns-3的离散事件仿真机制构建M/M/1和M/G/1排队系统，通过实验和数据分析，掌握这两个经典排队模型的特性。同时，将学习ns-3中随机变量和随机变量流的使用方法。本次实验仅涉及离散事件仿真，不涉及具体网络协议。
实验目标

• 掌握Little 定律，M/M/1 和 M/G/1 排队系统的定义及其时延特性；
• 了解ns-3 的离散事件仿真原理；
• 熟悉ns-3 中随机变量（流）的使用，了解蒙特卡洛仿真方法；
• 进一步提高文档和文献查阅能力，以及数据处理和可视化能力。

预备知识

通信网络中的时延是指数据分组从一端发出到另一端正确接收所花费的时间。端到端时延是通信网络中两个最重要的性能指标之一（另一个是吞吐量），是影响用户体验的重要因素。网络时延通常包括数据处理时延、传输时延、传播时延、排队时延等。处理时延是指分组发送前处理相关信息所需的时间，包括数据校验、控制信息处理、路由计算等，取决于计算芯片的处理速率；传输时延是指将分组中比特全部传输至链路上的时间，取决于网络端口的物理传输速率；传播时延是指数据在链路上传播的时间，取决于链路距离和信号的传播速度；排队时延是指分组到达一个节点的传输队列到该分组实际被传输的时间，它受到该节点分组到达速率、排队规则、发送速率等多重影响。

处理时延、传输时延、传播时延是相对确定的，主要取决于网络硬件的物理条件，而排队时延则取决于网络实时条件和主观设计。因此，分析排队时延是评估和优化网络设计的关键。排队论（Queueing Theory）是分析排队时延的主要工具，它是概率论下的一个数学分支，由丹麦电信工程师厄尔朗（Erlang）首先提出。排队论广泛应用于电信、计算机网络、交通工程、生产运筹等用户资源共享的随机服务系统的性能分析和设计。

排队论用A/S/C 表示法来建模排队系统，其中

• A表示到达过程
• S表示服务过程
• C 表示服务者的个数
  M/M/1和M/G/1是最常用的排队系统模型。M/M/1表示到达间隔和服务时间均服从指数分布，服务者数为1。指数分布是一个无记忆（Memoryless）分布，因此用M表示。到达间隔服从指数分布的随机过程又称为泊松（Poisson）过程（因为此时到达用户数服从泊松分布）（无记忆特性）。M/G/1则表示到达过程是泊松过程，但服务时间服从一般性（Generic）分布。

几个重要的数量指标

1. 平均队长
   $$Ls=\frac{\rho }{1-\rho }=\frac{\lambda }{\mu -\lambda }$$
2. 平均排队长
   $$Lq=L-\rho =\frac{{\lambda }^2}{\mu (\mu -\lambda )}$$
3. 平均逗留时间
   $$Ws=\frac{1}{\mu -\lambda }$$
4. 平均等待时间
   $$Wq=Ws-\frac{1}{\mu }=\frac{\lambda }{\mu (\mu -\lambda )}$$
5. 重要关系
   $$Ls=\lambda WsLq=\lambda Wq$$
6. 忙期平均时长
   $$B=\frac{\rho }{(\mu \times (1-\rho ))}$$
7. 系统利用率
   $$\rho =\frac{\lambda }{\mu }$$

蒙特卡洛仿真方法

蒙特卡罗方法也称统计模拟方法，是指使用随机数（或者更常见的伪随机数）来解决很多计算问题的方法。他的工作原理就是两件事：不断抽样、逐渐逼近。

实验内容

本次实验基于提供的queue-system.cc文件，将其置于scratch/目录下。运行命令./ns3 run "queue-system --numPackets=100"，确认仿真脚本可以运行。
上述仿真脚本构建一个先到先服务（FCFS）排队系统，分别设置分组到达间隔和服务时间的分布，并依赖随机变量的取值递归地调度到达和离开事件，实现排队系统的仿真。

问题 3.1

阅读脚本，指出服务事件（QueueSystem::Service函数）是何时初次被调度？
Service函数首次被调度的时机是在Arrival函数中，当队列从空状态变为非空状态时。具体地，在Arrival函数的以下代码段中：

if (m_queue->GetNPackets () == 1)   
{  
  // 当队列中第一个数据包到达时，调度Service函数  
  Simulator::Schedule (Seconds (m_var_deq->GetValue ()), &QueueSystem::Service, this);  
}

问题 3.2

由main()函数入手阅读仿真脚本，理解代码并简要文字描述enqueueVar和dequeueVar变量的作用

enqueueVar 和 dequeueVar 是用于模拟队列系统中数据包到达时间间隔和数据包服务时间间隔的随机变量对象。

enqueueVar 用于设置数据包到达的时间间隔的概率分布，它通过指数随机变量来模拟数据包的到达率。在仿真过程中，根据指定的到达率 lambda，enqueueVar 生成数据包到达的时间间隔，以便模拟实际情况下数据包的到达情况。

dequeueVar 则用于设置数据包服务时间的概率分布，它同样使用指数随机变量来模拟数据包的服务时间间隔。根据指定的服务率 mu，dequeueVar 生成数据包服务的时间间隔，以模拟数据包在队列中的服务情况。

这两个随机变量对象的设置使得队列系统能够更好地模拟实际情况下数据包到达和服务的随机性，从而更准确地评估队列系统的性能和行为。

问题 3.3

试使用Python 或 Matlab 工具处理得到的 queuelog.dat，绘出该日志中到达时间间隔的经验（empirical）分布函数。试说明经验分布函数与累积分布函数在概念上的区别和联系。

import numpy as np
import matplotlib.pyplot as plt
from statsmodels.distributions.empirical_distribution import ECDF

# 从文件中读取数据
data = np.genfromtxt('queuelog.txt', delimiter=' ', dtype=str)

# 提取到达时间间隔
arrivals = data[data[:, 1] == '+'][:, 0].astype(float)
intervals = np.diff(arrivals)

# 计算经验分布函数
ecdf = ECDF(intervals)

# 绘制经验分布函数
plt.plot(ecdf.x, ecdf.y, marker='o', linestyle='-')
plt.xlabel('Arrival Time Interval')
plt.ylabel('Empirical Cumulative Probability')
plt.title('Empirical Cumulative Distribution Function of Arrival Time Interval')
plt.grid(True)
plt.show()

这段代码首先从文件中读取数据，提取到达时间间隔，然后使用 ECDF 类计算经验分布函数，并绘制图形。

经验分布函数与累积分布函数的区别在于，经验分布函数是通过观测数据计算出的样本累积分布函数，用于描述随机变量的经验性质；而累积分布函数则是对随机变量的理论分布进行描述，给出了随机变量小于或等于某个值的概率。它们的联系在于，当样本数量充分大时，经验分布函数会趋近于理论的累积分布函数。

到达时间点

import numpy as np
import matplotlib.pyplot as plt
from statsmodels.distributions.empirical_distribution import ECDF

# 从文件中读取数据
data = np.genfromtxt('queuelog.txt', delimiter=' ', dtype=str)

# 提取到达时间间隔
arrivals = data[data[:, 1] == '+'][:, 0].astype(float)
#intervals = np.diff(arrivals)

# 计算经验分布函数
ecdf = ECDF(arrivals)

# 绘制经验分布函数
plt.plot(ecdf.x, ecdf.y, marker='.', linestyle='-')
plt.xlabel('Arrival Time')
plt.ylabel('Empirical Cumulative Probability')
plt.title('Empirical Cumulative Distribution Function of Arrival Time')
plt.grid(True)
plt.show()

问题 3.4

处理queuelog-i.dat，分别计算10000个分组的停留时间（即到达时间和离开时间之差），进而计算平均停留时间。同时，计算离开（-）事件队列长度的平均值。

#计算平均停留时间 python code
#计算离开（-）事件队列长度的平均值

import numpy as np
sum_all = 0
queue_lenth_all = 0
#从文件中读取数据

for i in range(16):
    data = np.genfromtxt('queuelog{}.txt'.format(i), delimiter=' ', dtype=str)
    data = np.genfromtxt('queuelog{}.txt'.format(i), delimiter=' ', dtype=str)
    #提取到达时间间隔和离开时间间隔
    arrivals = data[data[:,1] == '+'][:,0].astype(float)
    leaves = data[data[:,1] == '-'][:,0].astype(float)
    queue_lenth = data[data[:, 1] == '-'][:,2].astype(float)

    sum_stop = 0
    for i in range(len(leaves)):
        sum_stop += leaves[i]-arrivals[i]
    sum_all += sum_stop / (len(leaves))

    sum1=0
    for i in range(len(queue_lenth)):
        sum1  += queue_lenth[i]
    queue_lenth_all += (sum1/len(queue_lenth))    

print('平均停留时间：{}'.format(sum_all/16))
print('离开事件队列平均长度：{}'.format(queue_lenth_all/16))

平均停留时间：1.0143895302725001
离开事件队列平均长度：1.0190312499999998

问题3.5

计算每次仿真的分组平均停留时间和离开事件时队列长度的平均值，进而求得该100 次仿真的平均停留时间和平均队列长度的平均值，并与M/M/1排队系统的平均停留时间和平均用户数的理论值对比。

平均停留时间理论值
$$Ws=\frac{1}{\mu -\lambda }=1$$
平均用户数理论值
$$Ls=\frac{\rho }{1-\rho }=\frac{\lambda }{\mu -\lambda }=1$$

有模拟和理论值对比，数值大致接近。

问题 3.6

计算到达（+）事件时队列长度的平均值，与离开事件时的队列长度平均值比较，说明差异？两者哪一个才是正确的系统平均用户数？为什么？

#计算离开（-）事件队列长度的平均值

import numpy as np
queue_lenth_all = 0
#从文件中读取数据

for i in range(16):
    data = np.genfromtxt('queuelog{}.txt'.format(i), delimiter=' ', dtype=str)
    data = np.genfromtxt('queuelog{}.txt'.format(i), delimiter=' ', dtype=str)
    
    queue_lenth = data[data[:, 1] == '+'][:,2].astype(float)

    sum1=0
    for i in range(len(queue_lenth)):
        sum1  += queue_lenth[i]
    queue_lenth_all += (sum1/len(queue_lenth)) 

print('到达事件队列平均长度：{}'.format(queue_lenth_all/16))

到达事件队列平均长度：2.0190312500000003

到达事件队列平均长度：2.0190312500000003，离开事件队列平均长度：1.0190312499999998。系统平均用户数理论值对比可以看出离开事件队列长度与系统平均用户数接近，由此离开事件队列长度才是正确的系统平均用户数。到达与离开事件队列长度之间相差大约1倍左右。

原因：首先，系统平均用户数是指在长时间内，系统中顾客数的平均值，这包括了正在接受服务的顾客和等待服务的顾客。它反映了系统在稳定状态下的负载情况。

其次，离开事件队列长度是指当顾客完成服务并离开系统时，系统中剩余的顾客数。在稳态条件下，系统的到达率和离开率大致相等，这意味着系统中的顾客数在长时间内会趋于一个稳定值。因此，离开事件队列长度实际上也反映了系统在稳定状态下的顾客数情况。

问题 3.7

将问题3.5前仿真命令中的mu由2改为10，分别取lambada为1,3,5,7,9运行仿真。绘出以lambda/mu为横轴，分别以对应的平均停留时间和平均队列长度为纵轴的曲线。将结果与M/M/1排队系统的理论结果做比较，验证Little定理

Little定理
利特尔法则，基于等候理论，可用于一个稳定的，非抢占式的系统中。其内容为：
在一个稳定（稳态）的系统（L）中，长期的平均顾客人数，等于长期的有效抵达率lambda，乘以顾客在这个系统中平均的停留时间。可以用一个代数式来表达：
$$Ls=\lambda Ws$$
$$N=\lambda T$$

• 非常一般性的场景
• 稳态状态下

./ns3 run "queue-system --numPackets=10000 --lambda=3 --mu=10 --RngRun=$RANDOM

#计算平均停留时间 python code
#计算离开（-）事件队列长度的平均值

import numpy as np
#从文件中读取数据

for j in range(1,10,2):
    data = np.genfromtxt('sec3.7_{}.txt'.format(j), delimiter=' ', dtype=str)
#    data = np.genfromtxt('sec3.7_{}.txt'.format(j), delimiter=' ', dtype=str)
    #提取到达时间间隔和离开时间间隔
    arrivals = data[data[:,1] == '+'][:,0].astype(float)
    leaves = data[data[:,1] == '-'][:,0].astype(float)
    queue_lenth = data[data[:, 1] == '-'][:,2].astype(float)
    
    sum_all = 0
    queue_lenth_all = 0
    sum_stop = 0
    
    for i in range(len(leaves)):
        sum_stop += leaves[i]-arrivals[i]
    sum_all += sum_stop / (len(leaves))

    sum1=0
    for i in range(len(queue_lenth)):
        sum1  += queue_lenth[i]
    queue_lenth_all += (sum1/len(queue_lenth))    

    print('lambda={}平均停留时间：{}'.format(j,sum_all))
    print('lambda={}平均队列长度：{}'.format(j,queue_lenth_all))

lambda=1平均停留时间：0.11039593699999757
lambda=1平均队列长度：0.1165
lambda=3平均停留时间：0.1416565460799991
lambda=3平均队列长度：0.4202
lambda=5平均停留时间：0.19700832930000037
lambda=5平均队列长度：0.9708
lambda=7平均停留时间：0.2990537547000008
lambda=7平均队列长度：2.0326
lambda=9平均停留时间：0.9792642193000034
lambda=9平均队列长度：8.786

平均停留时间理论值
$$Ws=\frac{1}{\mu -\lambda }$$
平均队列长度理论值
$$Ls=\frac{\rho }{1-\rho }=\frac{\lambda }{\mu -\lambda }$$

import numpy as np  
import matplotlib.pyplot as plt  
  
# 假设服务率μ为1  
mu = 10  
  
# 已知的λ值和对应的平均停留时间和平均队列长度  
lambdas = [1, 3, 5, 7, 9]  
average_waiting_times = [0.11039593699999757, 0.1416565460799991, 0.19700832930000037, 0.2990537547000008, 0.9792642193000034]  
average_queue_lengths = [0.1165, 0.4202, 0.9708, 2.0326, 8.786]  
  
# 计算理论上的平均停留时间和平均队列长度  
theoretical_waiting_times = [1 / (mu - l) for l in lambdas if l < mu]  
theoretical_queue_lengths = [l / (mu - l) for l in lambdas if l < mu]  

# 绘制平均停留时间曲线  
plt.figure(figsize=(10, 5))  
plt.subplot(1, 2, 1)  
plt.plot(lambdas, average_waiting_times, label='average_waiting_times')  
plt.plot(lambdas[:len(theoretical_waiting_times)], theoretical_waiting_times, label='theoretical_waiting_times', linestyle='--')  
plt.xlabel('λ ')  
plt.ylabel('average_waiting_times')  
plt.title('average_waiting_times about lambada')  
plt.legend()  
  
# 绘制平均队列长度曲线  
plt.subplot(1, 2, 2)  
plt.plot(lambdas, average_queue_lengths, label='average_queue_lengths')  
plt.plot(lambdas[:len(theoretical_queue_lengths)], theoretical_queue_lengths, label='theoretical_queue_lengths', linestyle='--')  
plt.xlabel('λ')  
plt.ylabel('average_queue_lengths')  
plt.title('average_queue_lengths about lambada')  
plt.legend()  
  
# 显示图形  
plt.tight_layout()  
plt.show()

问题 3.8

试通过处理某个到达速率为λ=1，服务速率为µ=2时的 queuelog.dat，计算出对应的M/M/1排队系统的忙期平均时长。试查阅文献，找出M/M/1排队系统的平均忙期时长的理论公式，与仿真结果对比。

队列长度由 0 变为 1 的时刻至队列长度重新变为 0 的时刻之间的阶段被称为排队系统的忙期（busyperiod）。忙期时长是排队论中一个重要分析指标。M/M/1 队列的忙期时长是一个典型的更新（renewal）过程，即每个忙期时长服从独立同分布。

对于M/M/1模型，忙期平均时长B可以通过一下公式计算：
$$B=\frac{\rho }{(\mu \times (1-\rho ))}=0.5$$

import numpy as np


sum_busy = 0
# 从文件中读取数据
data = np.genfromtxt('queuelog.txt', delimiter=' ', dtype=str)


# 提取到达时间间隔
busy_time = data[data[:, 2] != '0'][:, 0].astype(float)

busy_time_intervals = np.diff(busy_time)

for i in range(len(busy_time_intervals)):
    sum_busy += busy_time_intervals[i]

aver_busytime = sum_busy/len(busy_time)
print('平均时长为：{}'.format(aver_busytime))

平均时长为：0.6641316519535938

问题 3.9

查阅ns-3 在线文档，用 SetAttribute函数设置该均匀分布随机变量的上下界（即Min 和 Max属性）分别为0和2/µ，使得其均值与原指数分布一致（1/µ）。取λ=1，µ =2，重复问题3.5。计算每次仿真的分组平均停留时间和离开事件时队列长度的平均值，进而求得该100 次仿真的平均停留时间和平均队列长度的平均值，并与M/M/1排队系统的平均停留时间和平均用户数的理论值对比。

将 queue-system.cc 中的 dequeueVar 变量的分布由指数分布改为其他分布，我们可以仿真M/G/1排队系统。例如，可将 dequeueVar变量由指数分布改为均匀分布。Ns-3中内置的均匀分布随机变量流名为 UniformRandomVariable，即修改相应行为：

// 设置服务时间的均匀分布，使其均值等于1/mu  
Ptr<UniformRandomVariable> dequeueVar = CreateObject<UniformRandomVariable> ();  
dequeueVar->SetStream (2);  
double minUniform = 0.0; // 均匀分布的最小值  
double maxUniform = 1.0 / mu; // 均匀分布的最大值，使得均值等于1/mu  
dequeueVar->SetAttribute ("Min", DoubleValue (minUniform));  
dequeueVar->SetAttribute ("Max", DoubleValue (maxUniform));  
qsys->SetServiceDist (dequeueVar);

./ns3 run "queue-system --numPackets=10000 --lambda=1 --mu=2 --RngRun=$RANDOM"

#计算平均停留时间 python code
#计算离开（-）事件队列长度的平均值

import numpy as np
sum_all = 0
queue_lenth_all = 0
#从文件中读取数据

for i in range(1):
    data = np.genfromtxt('sec3.9.txt', delimiter=' ', dtype=str)

    #提取到达时间间隔和离开时间间隔
    arrivals = data[data[:,1] == '+'][:,0].astype(float)
    leaves = data[data[:,1] == '-'][:,0].astype(float)
    queue_lenth = data[data[:, 1] == '-'][:,2].astype(float)

    sum_stop = 0
    for i in range(len(leaves)):
        sum_stop += leaves[i]-arrivals[i]
    sum_all += sum_stop / (len(leaves))

    sum1=0
    for i in range(len(queue_lenth)):
        sum1  += queue_lenth[i]
    queue_lenth_all += (sum1/len(queue_lenth))    

print('平均停留时间：{}'.format(sum_all))
print('离开事件队列平均长度：{}'.format(queue_lenth_all))

平均停留时间：0.30212277208000055
离开事件队列平均长度：0.3005

M/M/1排队系统的平均停留时间和平均用户数的理论值
平均停留时间理论值
$$Ws=\frac{1}{\mu -\lambda }=1$$
平均用户数理论值
$$Ls=\frac{\rho }{1-\rho }=\frac{\lambda }{\mu -\lambda }=1$$

问题3.10

将 dequeueVar 变量改为常值“随机”变量（ConstantRandomVariable），设置其Constant 属性为 1/µ。取 λ=1，µ=2，重复问题3.5。

  // Set service time probability distribution
  // Exponential service time
  Ptr<ExponentialRandomVariable> dequeueVar = CreateObject<ExponentialRandomVariable> ();
  dequeueVar->SetStream (2);
  dequeueVar->SetAttribute ("Mean", DoubleValue (1/mu));
  qsys->SetServiceDist (dequeueVar);

  Ptr<ConstantRandomVariable> dequeueVar = CreateObject<ConstantRandomVariable> ();
  dequeueVar->SetStream (2);
  double constantValue = 1.0 / mu;  
  dequeueVar->SetAttribute ("Constant", DoubleValue(constantValue));
  qsys->SetServiceDist (dequeueVar);
  
./ns3 run "queue-system --numPackets=10000 --lambda=1 --mu=2 --RngRun=$RANDOM"

由运行结果可知，仿真结果与理论结果相差较大

#计算平均停留时间 python code
#计算离开（-）事件队列长度的平均值

import numpy as np
sum_all = 0
queue_lenth_all = 0
#从文件中读取数据

for i in range(1):
    data = np.genfromtxt('sec3.10.txt', delimiter=' ', dtype=str)

    #提取到达时间间隔和离开时间间隔
    arrivals = data[data[:,1] == '+'][:,0].astype(float)
    leaves = data[data[:,1] == '-'][:,0].astype(float)
    queue_lenth = data[data[:, 1] == '-'][:,2].astype(float)

    sum_stop = 0
    for i in range(len(leaves)):
        sum_stop += leaves[i]-arrivals[i]
    sum_all += sum_stop / (len(leaves))

    sum1=0
    for i in range(len(queue_lenth)):
        sum1  += queue_lenth[i]
    queue_lenth_all += (sum1/len(queue_lenth))    

print('平均停留时间：{}'.format(sum_all))
print('离开事件队列平均长度：{}'.format(queue_lenth_all))

平均停留时间：0.7372172213000007
离开事件队列平均长度：0.7194

问题3.11

理论计算上述两种随机变量的二阶矩，代入M/G/1 队列的 P-K 公式计算平均停留时间和平均队列长度的理论值，验证仿真结果。

M/G/1排队系统与M/M/1系统的主要差别是服务时间为一般性的独立同分布的服务时间。

M/G/1排队系统
假定第$i$个用户的服务时间为$X_i$,$X_i$是独立同分布的，并且与到达间隔相互独立。令$X=X_1,X_2,\cdots$，则服务时间的均值和二阶矩为：
$$平均服务时间=\overline{X}=EX=\frac{1}{\mu }$$
$$服务时间的二阶矩=\overline{X^2}=E{X}^2$$
M/G/1排队系统的平均等待时间为
$$W=\frac{\lambda \overline{X^2}}{2(1-\rho )}$$
式中，$\rho =\frac{\lambda }{\mu }=\lambda \overline{X}$。该式称为$P-K$(Pollaczek-Khinchin)公式。
根据$P-K$公式，可以得到该系统的平均时延为
$$T=\overline{X}+W=\overline{X}+\frac{\lambda \overline{X^2}}{2(1-\rho )}$$
应用Little定理，有平均队列长度为
$$N_Q=\lambda W=\frac{{\lambda }^2\overline{X^2}}{2(1-\rho )}$$

均匀分布随机变量均值$\overline{X}=\frac{1}{\mu }=\frac{1}{2}=\frac{a+b}{2}=0.5$，二阶矩$\overline{X^2}=E{X}^2-E(X{)}^2=0.33$

常值“随机”变量均值$\overline{X}=\frac{1}{\mu }=0.5$，二阶矩$\overline{X^2}=\frac{1}{{\mu }^2}=0.25$

代入M/G/1 队列的 P-K 公式计算平均停留时间和平均队列长度的理论值

均匀分布随机变量

平均停留时间为
$$W=\frac{\lambda \overline{X^2}}{2(1-\rho )}=\frac{1\times 0.25}{2\times (1-0.5)}=0.25$$

平均队列长度为
$$N_Q=\lambda W=\frac{{\lambda }^2\overline{X^2}}{2(1-\rho )}=\frac{1\times 0.25}{2\times (2\times 0.5)}=0.25$$

常值“随机”变量
有$\rho =\frac{\lambda }{\mu }=0.5$
$$W=\frac{\rho }{2\mu (1-\rho )}$$

平均停留时间为
$$W=\frac{\rho }{2\mu (1-\rho )}=\frac{0.5}{2\times 2\times 0.5}=0.25$$
平均队列长度为
$$N_Q=\lambda W=1\times 0.25=0.25$$