Estimating Tail Probability in MMPP/D/1 Queue with Importance Sampling by Service Rate Adjustments

Abstract

The Asynchronous Transfer Mode (ATM) is an efficient technology for call relays, and it transmits information from multiple services including data, video, or voice. This information is conveyed at ATM multiplexers in small fixed-size packets called cells. The acceptable cell loss probability at ATM multiplexers is about ${10}^{-12}$. Important Sampling (IS) is an efficient method for estimating tiny probabilities that cannot be achieved by traditional Monte Carlo (MC) methods. This research presents a novel approach for evaluating the tail probability in the MMPP/D/1 queue system utilizing importance sampling simulation in the ATM network. To generate more rare events, a virtual queue is implemented in the dequeue process by decreasing the processing rate in the queue. In this way, the tail probability can be estimated on a real-time network.

1. Introduction

Asynchronous Transfer Mode (ATM) networks are connection-oriented networks for cell relay that support some services such as video, voice streaming, and data communications. At a multiplexer in an ATM network, the admissible cell loss probability needs to be less than ${10}^{-12}$ for real-time services such as voice and video. Evaluating this extremely small probability by simulation is a challenging task. Monte Carlo (MC) is an effective technique for performing simulations of problems that cannot be solved theoretically. However, the disadvantage of the Monte Carlo method is that it consumes a large amount of simulation time. The Importance Sampling technique is useful for simulations in the fields of communication and engineering [1,2]. As an improvement of the Monte Carlo method, IS simulation generated the target events more frequently by using another distribution density function.

It is crucial to ascertain the most advantageous settings for the distribution in the simulation. The distribution that produces the estimate with the lowest variance is called the optimal simulation distribution. The determination of optimal distribution in IS simulation has been studied in many papers, including heuristic and theoretical methods. In [3], the authors proposed a heuristic method to search parameters for optimal simulation distribution. Another method was proposed in [1], using an annealing algorithm to obtain optimal simulation distribution. These methods search the large domain of parameters, so they require computation time. The author in [4] introduced a theoretical approach that utilizes big deviation theory to determine the most efficient simulation distribution for the tail probability in the MMPP/D/1 queue system. Then [5] solved the general case of n states of the MMPP/D/1 queue. These theoretical methods have less computing time compared with the heuristic methods.

The MMPP, also known as the Markov-Modulated Poisson Process, is a stochastic process that combines both a Poisson process and a Markov chain. In this process, the intensity of the Poisson process is determined by the current state of the Markov chain. The MMPP is capable of modeling the sudden increase in video, voice, and data traffic on the Internet. Let us examine the two-state Markov-Modulated Poisson Process (MMPP). The symbol ${\lambda }_i$ represents the mean arrival rate of the Poisson process in state i, where i can take on the values of either 0 or 1. Additionally, $r_{ij}$ represents the probability of transitioning from state i to state j, where both i and j can be either 0 or 1. Figure 1 depicts the representation of a two-state Markov-Modulated Poisson Process (MMPP).

The MMPP/D/1 queue model is characterized by an arrival process following the MMPP model, deterministic service times, and an unlimited buffer capacity. Let $\mu$ be the service rate of the MMPP/D/1 queue. Therefore, $1/\mu$ represents both the duration of serving a single packet and the size of a time slot.

The MMPP/D/1 queue is a suitable model for simulating the arrival of packets in an ATM network multiplexer. Due to the minuscule likelihood of packet loss in an ATM network, determining it using the traditional Monte Carlo method is a time-consuming process. We employed the IS (Importance Sampling) technique to compute the probability of packet loss in the MMPP/D/1 queue through an ATM network simulation. By decreasing the processing speed at the queue, we may carry out estimations in real time.

Kobayashi [6] proposed a method to estimate the Tail Probability of FIFO queue length by using Importance Sampling. A virtual queue, ISQL counter, was implemented to decrease the service rate. Consequently, the queue length increases quickly and tail events happen more frequently. In this way, the estimation of the Tail Probability can be calculated faster and more exactly.

In our earlier research, we developed an IS approach for evaluating the possibility of packet loss in the MMPP/D/1 queue, with the intention of reducing the processing rate as proposed by Kobayashi (2014) and cited in Hung (2018). The most efficient distribution of IS is achieved by reducing the service rate. The MMPP/D/1 queue system is more intricate compared with the FIFO queue. The traditional approach [4,5] modified the pace at which events occur in order to produce a greater number of rare events. This necessitates altering user traffic, which is currently impossible within the existing network infrastructure.

In this paper, we have more extensive research than our, previous study [7]. The first section, Section 4.1, provides mathematical proof explaining the optimal processing rate formula. The second section, Section 4.3.2, explains the relationship of the MMPP’s states in the IS estimation; as a result, these states will be removed to gain more simulation quickly and can be applied in an online (real-time) network. These are two main contributions in this paper.

The next sections of this work are structured in the following manner: In Section 3, a review of the fundamental concepts of Monte Carlo simulation and Importance Sampling simulation is provided. Section 4 provides a comprehensive explanation of the modeling of the two-state MMPP/D/1 queue. In Section 5, we assess the effectiveness of our suggested method using simulations. Section 6 presents the conclusions and further research.

2. Variance Reduction in Importance Sampling

As we know, the Monte Carlo method takes a long time to simulate, especially since the accuracy of this method is not reliable in cases where the probability is extremely small. Consider the general case of estimating the d-dimensional integral as follows:

$I_{\varphi }={\int }_{{\mathbb{R}}^d}\varphi (x)f(x)dx$

where $f:{\mathbb{R}}^d\to \mathbb{R}$ is the probability density function of a random variable X and $\varphi :{\mathbb{R}}^d\to \mathbb{R}$ is a positive integral function. Let us consider the simple case when $\varphi ={\mathbf{1}}_A$ with A is a subset of ${\mathbb{R}}^d$, as follows:

${\mathbf{1}}_A=\left\{\begin{array}{cc}\hfill 1,\hfill & \hfill \mathrm{if}X\in A,\hfill \\ \hfill 0,\hfill & \hfill \mathrm{if}X\notin A\hfill \end{array}\right.$

Thus, $I_{\varphi }$ is the probability that $X\in A$ with X is distributed according to the probability density function f. The Monte Carlo estimation of the integral $I_{\varphi }$ is given by

$I_{\varphi }^{MC}={\displaystyle \frac{1}{N}}{\displaystyle \sum _{i=1}^N}\varphi (X_i)$

where $X_1,X_2,\cdots ,X_N$ are random samples generated by function f. The variance of $I_{\varphi }^{MC}$ is given by

$Var(I_{\varphi }^{MC})={\displaystyle \frac{1}{N}}\left({\int }_{{\mathbb{R}}^d}\varphi {(x)}^{2}f(x)dx-I_{\varphi }^2\right)$

Since $\varphi (x)={\mathbf{1}}_A$, we have $\varphi {(x)}^2=\varphi (x),\forall x$, and

${\int }_{{\mathbb{R}}^d}\varphi {(x)}^{2}f(x)dx={\int }_{{\mathbb{R}}^d}\varphi (x)f(x)dx$

We have the variance of $I_{\varphi }^{MC}$ only depending on N and $I_{\varphi }$

$Var(I_{\varphi }^{MC})={\displaystyle \frac{1}{N}}\left(I_{\varphi }-I_{\varphi }^2\right)$

The standard deviation of $I_{\varphi }^{MC}$, ${\sigma }_{I_{\varphi }^{MC}}$ is given by

${\sigma }_{I_{\varphi }^{MC}}={\displaystyle \frac{1}{\sqrt{N}}}\sqrt{\left(I_{\varphi }-I_{\varphi }^2\right)}$

When the event $\{X\in A\}$ is a rare event, the probability $I_{\varphi }$ becomes very small, and the standard deviation of $I_{\varphi }^{MC}$ tends towards $\sqrt{{\displaystyle \frac{I_{\varphi }}{N}}}$. In this case, unless the sample size N takes very large values, $\sqrt{{\displaystyle \frac{I_{\varphi }}{N}}}$ is much larger than $I_{\varphi }$. Therefore, the Monte Carlo estimation is not adapted to estimate such low probabilities.

The IS technique is the alternative method to reduce the variance of Monte Carlo techniques without increasing the sample size N. In IS simulation, an auxiliary PDF h is used to generate the samples $X_1,X_2,\cdots ,X_N$, and then estimate $I_{\varphi }$ as follows:

$I_{\varphi }^{IS}={\displaystyle \frac{1}{N}}{\displaystyle \sum _{i=1}^N}\varphi (X_i){\displaystyle \frac{f(X_i)}{h(X_i)}}$

Then, the expected value of $I_{\varphi }^{IS}$ is given by

$E(I_{\varphi }^{IS})={\displaystyle \frac{1}{N}}{\displaystyle \sum _{i=1}^N}{\int }_{{\mathbb{R}}^d}\varphi (x_i){\displaystyle \frac{f(x_i)}{h(x_i)}}h(x_i)d{x}_i$

The term $I_{\varphi }^{IS}$ is an unbiased estimator of $I_{\varphi }$ since

$E(I_{\varphi }^{IS})={\int }_{{\mathbb{R}}^d}\varphi (x)f(x)dx=I_{\varphi }$

Denote $w(X)={\displaystyle \frac{f(X)}{h(X)}}$; the variance of $I_{\varphi }^{IS}$ given by

$\begin{array}{ccc}\hfill Var(I_{\varphi }^{IS})& =& {\displaystyle \frac{1}{N}}\left({\int }_{{\mathbb{R}}^d}\varphi {(x)}^{2}w{(x)}^{2}h(x)dx-I_{\varphi }^2\right)\hfill \\ & =& \frac{1}{N}(E(\varphi {(X)}^{2}w{(X)}^2)-I_{\varphi }^2)\hfill \end{array}$

Thus, the variance of the IS estimate notably depends on the choice of h function. If h is well chosen, then the variance of the IS estimate can be very low. However, conversely, if h is chosen and not adapted, the variance of the IS estimate can be higher than the Monte Carlo estimate. Therefore, by choosing the adapted h function, IS simulation can reduce the variance of Monte Carlo simulation.

3. MMPP/D/1 Simulation for Rare Event

This section examines the use of Monte Carlo simulation and Importance Sampling simulation to assess the the potential of the queue length exceeding a certain threshold. The variable Q denotes the length of the stationary queue, whereas $P(Q>q)$ reflects the likelihood that Q surpasses the value of q. As the value of q tends towards infinity, the probability $P(Q>q)$ declines dramatically, and the event $\{Q>q\}$ is considered unusual. Monte Carlo simulation involves rare occurrences that have a low probability of occurring, which requires a significant amount of time for the simulation to run. To tackle this problem, the conventional IS technique boosts the packet arrival rate to more frequently generate the event $\{Q>q\}$ in comparison with MC simulation. The main advantage of IS simulation is its capacity to reduce simulation time while simultaneously improving accuracy. The acronyms and their meanings are explained in Table 1.

3.1. Monte Carlo Simulation for Queue

Examine the MMPP (Markov-Modulated Poisson Process) as the input traffic, which consists of two states, each representing a Poisson process. Analyze the MMPP/D/1 queue system. Define $Q_t$ as the steady-state queue length and $S_t$ as the state of a Markov-Modulated Poisson Process (MMPP) at the t-th time slot. A time slot refers to the specific duration of time allocated for the processing of a packet. The Monte Carlo estimation for events $\{Q>q\}$ is performed by

(1)$P(Q>q)={\displaystyle \frac{1}{N}}{\displaystyle \sum _{t=1}^N}{\delta }_q(Q_t)$

where N is the number of time slots and ${\delta }_q(Q_t)$ is the indicator function given by

(2)${\delta }_q(Q_t)=\left\{\begin{array}{cc}\hfill 1,\hfill & \hfill \mathrm{if}{Q}_t>q,\hfill \\ \hfill 0,\hfill & \hfill \mathrm{if}{Q}_t\le q\hfill \end{array}\right.$

If the event $\{Q>q\}$ in Equation (1) is considered to be a rare occurrence, it is necessary to have a significantly large number of samples N. Therefore, to estimate the probability of a rare event, MC simulation takes much time to obtain the target value [8]. This is the limitation of MC simulation for estimating the probability of rare events.

3.2. IS Simulation for Queue

The goal of the IS method is to generate as many rare events as possible. In queue simulation, to increase the number of packets in the queue, one increases the arrival rate. This way, if the processing speed stays the same, the queue will quickly fill up and rare events will become more numerous. In [1], the author proposed the DIS (Dynamic Importance Sampling) method in which time slots are divided into Regenerative Cycle (RC) which described in Figure 2. Each RC is the time interval from when the queue length is greater than 0 to when its length returns to 0 (Figure 3). DIS is an important concept in IS simulation.

Now, we will discuss the process of estimating the tail probability using importance sampling (IS) simulation. Let $Q_{t_1}^{\prime },Q_{t_2}^{\prime },Q_{t_3}^{\prime },\dots$ represent the queue length associated with the time slot sequence $t_1,t_2,t_3,\dots$ according to the IS simulation distribution. Next, we will present the IS estimate of the tail probability as follows:

(3)$P_{IS}(Q>q)=\frac{\frac{1}{N}{\sum }_{k=1}^N{\sum }_{t_k=1}^{T_k}{\delta }_q(Q_{t_k}^{\prime })W_{t_k}}{\frac{1}{M}{\sum }_{k=1}^M{\sum }_{t_k=1}^{T_k^{\prime }}{W}_{t_k}}$

(4)$\begin{array}{c}\hfill W_{t_k}={\displaystyle \prod _{m_k=1}^{t_k}}\frac{P(Q_{m_k-1}^{\prime },Q_{m_k}^{\prime })}{P^{\prime }(Q_{m_k-1}^{\prime },Q_{m_k}^{\prime })}\end{array}$

where

• M and N represent the quantities of RCs;

• $t_k$ represents a specific time interval within the kth RC;

• $T_k$ and $T_k^{\prime }$ represent the length of the kth RC;

• $W_{t_k}$ represents the weighting function;

• The parts $P(Q_{m_k-1}^{\prime },Q_{m_k}^{\prime })$ and $P^{\prime }(Q_{m_k-1}^{\prime },Q_{m_k}^{\prime })$ are the state transition probability from $(m_k-1)$ to $m_k$, respectively.

In the next section, we will present how to determine the optimal IS distribution.

3.3. The Optimal IS Distribution

In IS simulation, the optimal distribution determination is most important. If we choose the good optimal distribution, the variance of the estimator will be smaller and the estimate result will be more exact. There are some approaches to determine the optimal simulation distribution in [1,3]. However, these methods are heuristic, exhaustive search in parameter space to determine the optimal distribution; therefore, they take much computation time. In this paper, based on the research in [4], we determine the optimal simulation distribution for IS simulation. This approach is theoretical; hence, it takes less computation time.

The state of the MMPP and the queue length are denoted by $S_t$ and $Q_t$, where t is the time slot, $t=0,1,2,\dots$. Denote $X_t=(Q_t,S_t)$; then $X_t$ satisfies the formula

$\begin{array}{c}{X}_0=X^0(\mathrm{initial}\mathrm{state})\\ X_{t+1}=X_t+\Delta (X_t,W)\end{array}$

where W is an independent and identically distributed random vector; therefore, a random process $X_t$ forms a Markov chain [5]. In this paper, for simplicity, $X_t$ is represented by X at an arbitrary t time slot and $M_X(\theta )$ is the moment-generating function of $\Delta (X,W)$. $M_X(\theta )$ is defined by

(5)$M_X(\theta )=E[e^{\theta \Delta (X,W)}]$

where E is the expectation with respect to $\Delta (X,W)$ and the parameter $\theta$.

$P(X_t,X_{t+1})$ is defined as the probability of moving from state $X_t$ to state $X_{t+1}$, according to the original distribution. The Markov chain is governed by the state transition probability, determined by the parameter $\theta$, denoting the probability of transitioning between states.

(6)$P_{\theta }(X_t,X_{t+1})=\frac{e^{\theta \Delta (X,W)}}{M_X(\theta )}P(X_t,X_{t+1})$

As stated in reference [9], the value $\theta ={\theta }^{*}$ that solves the equation $M_X(\theta )=1$ is considered the optimal distribution. The simulation distribution that is most advantageous, abbreviated as $P^{*}\equiv P_{{\theta }^{*}}$, is given by

(7)$P^{*}(X_t,X_{t+1})=e^{{\theta }^{*}\Delta (X,W)}P(X_t,X_{t+1})$

In [4], the authors used the large deviation theory to find the optimal IS distribution. In the next section, we will determine the optimal distribution for our IS simulation based on [4].

4. Proposed Method for Estimating the Tail Probability of the MMPP/D/1 Queue

4.1. Optimal Parameter for Estimation

Now, we are identifying the optimal parameter for estimating the Importance Sampling (IS) of the MMPP/D/1 queue. The queue length and the state of the Markov-Modulated Poisson Process (MMPP) at time slot t are denoted by $Q_t$ and $S_t$, respectively. The process $X_t=(Q_t,S_t)$ denotes a series of states at various time intervals, with t having a range of values of 0, 1, 2, and so forth. The probability of transitioning from state $X_t$ to state $X_{t+1}$ is represented by $P(X_t,X_{t+1})$. Due to the deterministic service time of the MMPP/D/1 queue system, the packet length remains constant. The disparity in the length of the queue between the tth and $(t+1)$st time slots is represented as y and can be calculated as the difference between $Q_{t+1}$ and $Q_t$.

We will provide a comprehensive elucidation on the process of determining the different queue lengths, represented as y, in the MMPP/D/1 queue. At time slot $(t+1)$, the state of the Markov-Modulated Poisson Process (MMPP) is represented as $S_{t+1}=j$, indicating that $A_j$ packets have been received. In this context, $A_j$ denotes a Poisson random variable that corresponds to state j and has an arrival rate of ${\lambda }_j$. If the value of $Q_t$ is greater than zero, then the value of y is equal to $A_j$ minus one. The equation $y=A_j-1$ represents the relationship between the number of packets served, denoted by y, and the number of consecutive slots, denoted by $A_j$. This equation accounts for the fact that there is exactly one packet served between two consecutive slots. If $Q_t=0$, then y is equal to $A_j$.

Denoting by $h=j-i,z=(y,h)$, we have the state transition probability $P(X_t,X_t+z)$ as follows:

(8)$P(X_t,X_t+z)=\left\{\begin{array}{cc}\hfill {\displaystyle \frac{{({\lambda }_j/\mu )}^{y+1}}{(y+1)!}}{e}^{-{\lambda }_j/\mu }{r}_{ij},\hfill & \hfill \mathrm{if}{Q}_t>0,\hfill \\ \hfill {\displaystyle \frac{{({\lambda }_j/\mu )}^y}{y!}}{e}^{-{\lambda }_j/\mu }{r}_{ij},\hfill & \hfill \mathrm{if}{Q}_t=0\hfill \end{array}\right.$

To obtain the optimal distribution simulation, we solve the equation

(9)$M_{X_t}({\theta }_0,{\theta }_1)=1$

According [4], we have ${\theta }_0={\theta }_0^{*}$, and ${\theta }_1={\theta }_1^{*}$ are the solutions of Equation (9). Hence, the optimal distribution simulation $P^{*}(X_t,X_t+z)$ is given by

(10)$P^{*}(X_t,X_t+z)=e^{{\theta }_0^{*}y}{e}^{{\theta }_1^{*}(j-i)}P(X_t,X_t+z)$

In case of $Q>0$, from (10), we have

${\displaystyle \frac{{({\lambda }_j^{*}/{\mu }^{*})}^{y+1}}{(y+1)!}}{e}^{-{\lambda }_j^{*}/{\mu }^{*}}{r}_{ij}^{*}={\displaystyle \frac{{({\lambda }_j/\mu )}^{y+1}}{(y+1)!}}{e}^{-{\lambda }_j/\mu }{r}_{ij}$

(11)$⟺{({\lambda }_j^{*}/{\mu }^{*})}^{y+1}{e}^{-{\lambda }_j^{*}/{\mu }^{*}}{r}_{ij}^{*}={({\lambda }_j/\mu )}^{y+1}{e}^{-{\lambda }_j/\mu }{r}_{ij}$

To compare both sides of (11), we have the optimal parameter ${\lambda }_i^{*}/{\mu }^{*}$ and $r_{ij}^{*}$ as follows:

(12)$\left\{\begin{array}{c}{\lambda }^{*}/{\mu }^{*}=e^{{\theta }_0}\lambda /\mu \hfill \\ r_{ij}^{*}=e^{-{\theta }_0-{\lambda }_j/\mu -{\lambda }_j^{*}/{\mu }^{*}+{\theta }_1(j-i)}{r}_{ij}.\hfill \end{array}\right.$

Similarly, we have the optimal parameter in case of $Q=0$ as follows:

(13)$\left\{\begin{array}{c}{\lambda }^{*}/{\mu }^{*}=\lambda /\mu \hfill \\ r_{ij}^{*}=r_{ij}.\hfill \end{array}\right.$

In (12), we can change the optimal service rate ${\lambda }^{*}$ and optimal service rate ${\mu }^{*}$ to obtain the optimal distribution simulation. In this case, the value of the arrival rate is kept unchanged; Equation (12) becomes

(14)$\left\{\begin{array}{c}{\lambda }_i^{*}={\lambda }_i\hfill \\ {\mu }^{*}=e^{-{\theta }_0^{*}}\mu \hfill \end{array}\right.$

Thus, the optimal processing rate ${\mu }^{*}=e^{-{\theta }_0^{*}}\mu$ was determined for the IS simulation.

4.2. The ISDC (Importance Sampling Dequeue Counter)

In the previous section, we determined the most favorable parameters for simulating the distribution of IS simulation. As a result, the service rate decreased by a factor of $e^{{\theta }_0^{*}}$. The updated service rate ${\mu }^{*}$ of the proposed method is as follows:

(15)${\mu }^{*}=e^{-{\theta }_0^{*}}\mu$

where $\mu$ is the conventional service rate of IS and MC simulations.

An important concept is called Importance Sampling Dequeue Counter (ISDC). The Algorithm 1 describes the operation of ISDC. The $ISDC$ counts the number of dequeue events. Denote $Q^{*}$ as the queue length in IS by using the optimal distribution. Next, the value of $Q^{*}$ will be computed using ISDC. Figure 4 illustrates the functioning of ISDC. The ISDC algorithm is defined as follows. Initially, the virtual IS queue length, denoted as $Q^{*}$, and the variable $counter$ are both set to 0. Additionally, the value of $step$ is determined as $e^{-{\theta }_0^{*}}$. The variable $counter$ will tally the occurrences of dequeue events in the IS virtual queue.

Algorithm 1 Algorithm of ISDC

1:. Set the initial values: 2:.       $Q^{*}=0;$ 3:.       $counter=0;$ 4:.       $step=e^{-{\theta }_0^{*}};$ 5:. For each packet that comes to the queue: 6:.       //increase the virtual queue length by 1 7:.       $Q^{*}=Q^{*}+1;$ 8:. For each packet that exits the queue: 9:.       $counter+=step;$ 10:.       If counter >= 1 then 11:.            $Q^{*}=Q^{*}-1;$ 12:.            $counter-=1;$ 13:.       End if

During the $enqueue$ event, the virtual queue length is incremented by 1, resulting in $Q^{*}=Q^{*}+1$, which corresponds to the addition of one packet to the virtual IS queue length $Q^{*}$. In this work, it is specified that the queue length is measured in packets. In this case, we add a single packet to the queue $Q^{*}$ as there is only one packet that has arrived.

During the $dequeue$ event, the algorithm operates in the following manner. Initially, the value of $counter$ is incremented by $step$. Next, evaluate the value of $counter$ in relation to 1. If $counter$ is more than or equal to 1, then one packet will be processed, and both $Q*$ and $counter$ will be reduced by 1.

By using $counter$ and $step$, the "virtual queue" reduces the processing rate $e^{{\theta }_0^{*}}$ times.

4.3. IS Estimation for Tail Probability

4.3.1. Conventional IS Estimation

In traditional importance sampling (IS) simulation, once the optimal parameters are identified, the tail probability estimate can be obtained using Equations (3) and (4), with ${\lambda }^{*}=e^{{\theta }_0^{*}}\lambda$.

(16)$P_{IS}^{*}=\frac{\frac{1}{N}{\sum }_{k=1}^N{\sum }_{t_k=1}^{T_k}{\delta }_q(Q_{t_k}^{*})W_{t_k}^{*}}{\frac{1}{M}{\sum }_{k=1}^M{\sum }_{t_k=1}^{T_k^{\prime }}{W}_{t_k}^{*}}$

(17)$\begin{array}{cc}\hfill W_{t_k}^{*}& \hfill ={\displaystyle \prod _{m_k=1}^{t_k}}\frac{P(Q_{m_k-1}^{*},Q_{m_k}^{*})}{P^{*}(Q_{m_k-1}^{*},Q_{m_k}^{*})}\hfill \end{array}$

From (10) and (17), the weighting function [10] becomes

(18)$W_{t_k}^{*}={\displaystyle \prod _{m_k=1}^{t_k}}{e}^{-{\theta }_0^{*}y}{e}^{-{\theta }_1^{*}(j-i)}$

where

• y represents the discrepancy in the lengths of the queues between the $(m_k-1)$th and $m_k$th time slots, given by the equation $y=Q_{m_k}^{*}-Q_{m_k-1}^{*}$;

• ${\theta }_0^{*}$ and ${\theta }_1^{*}$ are the ideal parameters;

• i and j represent the states of the Markov-Modulated Poisson Process (MMPP) at the $(m_k-1)$-th and $m_k$-th time slots, respectively;

• $T_k$ and $T_k^{\prime }$ represent the length of the RC kth;

• $t_k$ represents a specific time slot within a resource constraint;

• M and N represent the quantities of RCs in the simulation.

A time slot refers to the duration required to process a single packet. Therefore, y represents the disparity in the length of the IS queue between two successive time intervals.

4.3.2. IS Estimation

Let us consider Equation (18). In this equation, the parameters ${\theta }_0^{*}$ and ${\theta }_1^{*}$ are determined and calculated in the above sections based on the input data ${\lambda }_0,{\lambda }_1,r_{01},r_{10}$. However, the parameters $y,j,i$ must be calculated in real time. At the dequeue event, we can determine the queue length and compute the difference queue length y between successive dequeue events. However, determining the state of the MMPP in real time is very difficult. Nevertheless, to facilitate the calculation of the weighting function in real time, we removed the part $e^{-{\theta }_1^{*}(j-i)}$, and then (18) becomes

(19)$W_{t_k}^{*}={\displaystyle \prod _{m_k=1}^{t_k}}{e}^{-{\theta }_0^{*}y}$

Why we can remove it? To explain this, we investigated the events $\{j>i\}$ and $\{j<i\}$ during the simulation (in case of $\{j=i\}$, the value of $e^{-{\theta }_1^{*}(j-i)}=1$; hence, for each of the time slots that the state of the MMPP does not change, the value of the weighting function $W_{t_k}^{*}$ is not also changed). In fact, in the two-state MMPP/D/1 queue, if the state of the MMPP changes from state 0 to state 1, then it will be changed again from state 1 to state 0. This means that the number of transitions from state 0 to state 1 and the number of transitions from state 1 to state 0 are equivalent. On the other hand, the number of events $\{j>i\}$ and the number of events $\{j<i\}$ are equivalent. Therefore, the part $e^{-{\theta }_1^{*}(j-i)}$ can be removed from the weighting function without much impact on the final result. In the next section, we will examine and contrast the effectiveness of the proposed method with those of the traditional MC and IS methods.

5. Experiments

5.1. Performance Evaluation

We execute the proposed IS algorithm and evaluate the precision and duration of the tail probability calculation of the MMPP/D/1 queue using NS-2. NS-2 is a discrete event simulator utilized to simulate real-time network traffic and topology in order to conduct analysis. Because NS-2 is an open-source project and supports several algorithms in routing and queuing, we can modify the source code for our simulation easily.

Next, we compare the outcomes of our suggested approach with those of the MC method and the standard IS method. We conduct two scenarios to simulate the MMPP/D/1 queue. In Scenario-1, we consider the special case of the MMPP with two Poisson processes having the same arrival rate ${\lambda }_1={\lambda }_2$. Since the two-state MMPP has the same Poisson process, the MMPP/D/1 queue becomes the M/D/1 queue. Then, in Scenario-2, we consider the normal case of the MMPP with two Poisson processes with different arrival rates ${\lambda }_1\ne {\lambda }_2$.

For the two-state MMPP, the average arrival rate is calculated by

(20)${\lambda }_{mean}=\frac{{\lambda }_0{r}_{10}+{\lambda }_1{r}_{01}}{r_{01}+r_{10}}$

with the variables ${\lambda }_0$ and ${\lambda }_1$ representing the arrival rates of two Poisson processes. The variables $r_{01}$ and $r_{10}$ denote the probabilities of transitioning from state 0 to state 1 and vice versa.

5.1.1. Scenario-1

In this scenario, we set the arrival rate of the two-state MMPP to be the same, ${\lambda }_0={\lambda }_1$; then the MMPP/D/1 queue becomes the M/D/1 queue. The parameters for MC, conventional IS, and the proposed IS simulation are shown in Table 2, Table 3, Table 4, respectively.

These parameters ${\lambda }_0^{*},{\lambda }_1^{*},r_{01}^{*},r_{10}^{*}$ are calculated by Equation (12), and ${\lambda }_{mean}^{*}$ is calculated by

(21)${\lambda }_{mean}^{*}=\frac{{\lambda }_0^{*}{r}_{10}^{*}+{\lambda }_1^{*}{r}_{01}^{*}}{r_{01}^{*}+r_{10}^{*}}$

In this scenario, we run 3 cases of simulation (Case 1, Case 2, Case 3) with means of arrival rate for MC simulation of $0.7$, $0.5$, and $0.3,$ respectively. These parameters are chosen to spread evenly in the range [0, 1].

5.1.2. Scenario-2

For this scenario, we perform three simulations with different parameters, namely, Case 4, Case 5, and Case 6. Table 5 presents the parameters utilized in the Monte Carlo simulation, which is based on the original distribution.

Table 6 displays the parameters used in the optimal distribution simulation of the Importance Sampling (IS) scenario.

The parameters for the proposed technique are displayed in Table 7. In contrast with the parameters used in Monte Carlo simulation, the suggested technique preserves the values of ${\lambda }_0,{\lambda }_1,r_{01},r_{10}$ and only modifies the service rate parameter $\mu$. The service rate ${\mu }^{*}$ of the proposed IS simulation is determined using Equation (15).

5.2. Evaluate Performance

5.2.1. The Accuracy of Estimation

The simulation results for Scenario-1 and Scenario-2 are displayed in Figure 5 and Figure 6, respectively. In MC, conventional IS, and the proposed IS simulation, there are ${10}^6$ [cycles] of regenerative cycles. Figure 5 and Figure 6 demonstrate that the performance of the proposed method is comparable to that of the MC and traditional IS. Otherwise, the estimated value of the IS approach can be produced at a value smaller than ${10}^{-25}$, whereas the MC technique only estimates the probability $P(Q>q)$ with a minimum value of roughly ${10}^{-6}$. When q is sufficiently large, the event $\{Q>q\}$ does not happen frequently, making it impossible to estimate the value using the MC approach, or requiring a significant amount of time to acquire the estimated value. Consequently, in these situations, the IS approach is helpful.

Simulation time: We now want to evaluate the simulation times of the MC method, conventional IS, and proposed IS simulation. In Case 4, we will simulate to estimate the tail probability $P(Q>q)$ for $q=20[packets]$ for all three methods. The number of RCs in this experiment is ${10}^7$ for MC simulation and ${10}^2$ for both conventional IS and the proposed IS simulation. To calculate the sample mean and sample variance, we ran the simulation 30 times to obtain the estimated value of the tail probability. The outcome is presented in Table 8, where the estimation variances are nearly the same. The simulation time using Monte Carlo (MC) is 1136.81 s, whereas it is 0.39 s using conventional Importance Sampling (IS) and 0.35 s using the proposed IS simulation. Therefore, the proposed IS simulation has a simulation time that is approximately 3200 times faster than the MC simulation and is slightly shorter than traditional IS. The reason can be explained as follows: First, it is the difference between the calculation of the weight function of conventional IS and the proposed IS. Conventional IS simulation uses Equation (18) to calculate the weighting function. This equation is more complex than Equation (19) used in the proposed IS. In Equation (18), it takes time to determine the states {$i,j$} of the MMPP. Moreover, the second part $e^{-{\theta }_1^{*}(j-i)}$ of this equation also takes more time to calculate. Hence, the simulation time of conventional IS is longer than that of the proposed IS simulation.

5.2.2. Speed of Convergence Comparison

Next, we want to show the convergence of the estimation by MC simulation and the proposed IS simulation. To perform this, we simulate the parameters in Case 4 by changing a variety of numbers of cycles. For MC simulation, the number of cycles is ${10}^6$ and ${10}^7$ cycles, and for the proposed IS simulation, the number of cycles is ${10}^2$, ${10}^3$, and ${10}^4$ cycles. In both of the two methods, we set the queue threshold as $q=20[packets]$. Figure 7 shows the convenience of the simulations. We can see that the speed of convergence of the proposed IS is faster than the speed of convergence of MC. The results demonstrate that the proposed IS significantly reduces simulation time compared with MC simulation.

5.3. Discussion on the States of the MMPP in Conventional IS Simulation

We now discuss the removing state information of the MMPP in conventional IS simulation. As mentioned in Section 4.3, we can remove state $\{i,j\}$ from Equation (18), but the proposed IS estimated value is not changed much. In this section, we want to make an experiment to show the estimated value of the rare event probability. We run IS simulation without the state with the parameters in Scenario-2. The result shows that the estimated values of conventional IS before and after the removal of the state information are equivalent. The reason is explained in Section 4.3. On the other hand, we also compare the simulation time of conventional IS without the state information with that of conventional IS. In Case 4, we will simulate to estimate the probability $P(Q>q)$ with $q=20[packets]$, and the number of RCs is ${10}^2$ cycles. The result in Table 9 shows that the simulation time of conventional IS without the state information is slightly faster than the simulation time of conventional IS (about $5\%$ faster). The reason is that the formula of the weighting function of conventional IS without states is simpler than the weighting function formula of conventional IS (the part $e^{-{\theta }_1^{*}(j-i)},(i,j=0,1)$ is removed).

6. Conclusions

This work presents a proposed simulation of an IS (Importance Sampling) technique to speed up the Monte Carlo simulation for the MMPP/D/1 queue. In traditional Importance Sampling (IS), the arrival rate is augmented in order to generate a significant number of rare occurrences. Our idea does not include increasing the arrival rate, but rather decreasing the service rate. Using this method, we may calculate the chance of the queue length exceeding a certain value in an MMPP/D/1 queue.

The authors in [6] also introduced a method to estimate the tail probability of FIFO queue length in real time. Nevertheless, the authors employed a Poisson process for the arrival process, but in this research, we utilized an MMPP for the arrival process. The determination of the service rate in the MMPP arrival process is more intricate compared with a Poison procedure in a FIFO queue.

Our technique provides the benefit of not altering user traffic, making it suitable for use in the actual network. The results of the aforementioned simulation demonstrate that our suggested technique achieves accuracy that is comparable to traditional IS and MC methods. However, our proposed method’s simulation time is approximately 3200 times faster than that of MC.

In addition, in conventional IS, four parameters need to be changed: the arrival rates ${\lambda }_0^{*},{\lambda }_1^{*}$ and the state transition probabilities $r_{01}^{*},r_{10}^{*}$. However, we only change the service rate; therefore, our proposed method can be performed more easily than conventional IS.

We also investigated the effect of the state information in conventional IS simulation. We show that the formula for the weighting function (18) does not depend on states i and j. In this way, we can perform an IS simulation for the two-state MMPP/D/1 queue. Otherwise, if we remove the state information from Equation (18), then the simulation time will be slightly reduced.

In our future work, we will study n-state MMPP/D/1 in cases of finite buffer models and more complicated models. In addition, we will study some effective IS algorithms, such as Adaptive Importance Sampling (AIS) for queues.

Author Contributions

Conceptualization, N.N.H.; methodology, N.N.H. and P.T.G.; writing, N.N.H.; review and editing, N.H.A. and P.T.G.; funding acquisition, N.H.A. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Data is contained within the article.

Acknowledgments

Special thanks go to the editor and anonymous reviewers of this paper for their constructive comments.

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Two-state MMPP.

Figure 2. The Regenerative Cycle.

Figure 3. The Dynamic Importance Sampling method.

Figure 4. The virtual queue.

Figure 5. Scenario-1: performance comparison of estimates by MC, conventional IS, and the proposed method ([Forumla omitted. See PDF.]).

Figure 6. Scenario-2: performance comparison of estimates by MC, conventional IS, and the proposed ([Forumla omitted. See PDF.]).

Figure 7. Convergence comparison of estimates by MC and the proposed IS (in Case 4 with q = 20 [packets]).

Table 1

The acronyms and meaning.

Acronyms	Meaning
ATM	Asynchronous Transfer Mode
MC	Monte Carlo
IS	Importance Sampling
MMPP	Markov-Modulated Poisson Process
ISDC	Importance Sampling Dequeue Counter
RC	Regenerative Cycle

Table 2

Scenario-1: parameters for MC simulation.

Case	${\mathit{\lambda }}_{\mathbf{0}}$	${\mathit{\lambda }}_{\mathbf{1}}$	${\mathit{r}}_{\mathbf{01}}$	${\mathit{r}}_{\mathbf{10}}$	$\mathit{\mu }$	${\mathit{\lambda }}_{\mathit{mean}}$
1	0.7	0.7	0.1	0.4	1.0	0.7
2	0.5	0.5	0.7	0.4	1.0	0.5
3	0.3	0.3	0.1	0.4	1.0	0.3

Table 3

Scenario-1: IS simulation’s parameters.

Case	${\mathit{\lambda }}_0^{*}$	${\mathit{\lambda }}_1^{*}$	${\mathit{r}}_{01}^{*}$	${\mathit{r}}_{10}^{*}$	$\mathit{\mu }$	${\mathit{\lambda }}_{\mathit{mean}}^{*}$
1	1.38	1.38	0.1	0.4	1.0	1.38
2	1.76	1.76	0.7	0.4	1.0	1.76
3	2.36	2.36	0.1	0.4	1.0	2.36

Table 4

Scenario-1: Parameters for the proposed IS simulation.

Case	${\mathit{\lambda }}_0$	${\mathit{\lambda }}_1$	${\mathit{r}}_{01}$	${\mathit{r}}_{10}$	${\mathit{\theta }}_0^{*}$	${\mathit{\mu }}^{*}$
1	0.7	0.7	0.1	0.4	0.68	0.51
2	0.5	0.5	0.1	0.4	1.26	0.28
3	0.3	0.3	0.1	0.4	2.06	0.13

Table 5

Scenario-2: parameters for MC simulation.

Case	${\mathit{\lambda }}_0$	${\mathit{\lambda }}_1$	${\mathit{r}}_{01}$	${\mathit{r}}_{10}$	$\mathit{\mu }$	${\mathit{\lambda }}_{\mathit{mean}}$
4	0.8	0.5	0.1	0.4	1.0	0.74
5	0.6	0.3	0.1	0.4	1.0	0.54
6	0.4	0.3	0.7	0.3	1.0	0.33

Table 6

Scenario-2: IS simulation’s parameters.

Case	${\mathit{\lambda }}_0^{*}$	${\mathit{\lambda }}_1^{*}$	${\mathit{r}}_{01}^{*}$	${\mathit{r}}_{10}^{*}$	$\mathit{\mu }$	${\mathit{\lambda }}_{\mathit{mean}}^{*}$
4	1.38	0.86	0.07	0.50	1.0	1.31
5	1.71	0.85	0.04	0.63	1.0	1.66
6	2.65	1.99	0.57	0.43	1.0	2.27

Table 7

Scenario-2: parameters for the proposed IS simulation.

Case	${\mathit{\lambda }}_0$	${\mathit{\lambda }}_1$	${\mathit{r}}_{01}$	${\mathit{r}}_{10}$	${\mathit{\theta }}_0^{*}$	${\mathit{\mu }}^{*}$
4	0.8	0.5	0.1	0.4	0.54	0.58
5	0.6	0.3	0.1	0.4	1.05	0.35
6	0.4	0.3	0.7	0.3	1.89	0.15

Table 8

The simulation time’s comparison.

	MC Method	Conventional IS Method	Proposed IS Method
RCs	${10}^7$	${10}^2$	${10}^2$
Mean of estimate	$5.58\times {10}^{-6}$	$6.19\times {10}^{-6}$	$5.80\times {10}^{-6}$
Sample variance	$1.14\times {10}^{-12}$	$3.09\times {10}^{-13}$	$4.12\times {10}^{-13}$
Simulation time [s]	1136.81	0.39	0.35

Table 9

Simulation time of conventional IS without state information.

	Conventional IS	Conventional IS without States
Number of RCs	${10}^2$	${10}^2$
Mean of estimate	$6.19\times {10}^{-6}$	$6.53\times {10}^{-6}$
Sample variance	$3.09\times {10}^{-13}$	$7.17\times {10}^{-13}$
Simulation time [sec]	0.39	0.37

References

1. Devetsikiotis, M.; Townsend, J. Statistical optimization of dynamic importance sampling parameters for efficient simulation of communication networks. IEEE/ACM Trans. Netw.; 1993; 1, pp. 293-305. [DOI: https://dx.doi.org/10.1109/90.234852]

2. Bucklew, J. Large Deviation Techniques in Decision, Simulation, and Estimation; Wiley Series in Probability and Statistics; Wiley: Hoboken, NJ, USA, 1990.

3. Wang, Q.; Frost, V. Efficient estimation of cell blocking probability for ATM systems. IEEE/ACM Trans. Netw.; 1993; 1, pp. 230-235. [DOI: https://dx.doi.org/10.1109/90.222929]

4. Ogawa, K.; Nakagawa, K. The optimal IS Simulation Distribution of MMPP/D/1 Queueing. IEICE Trans. Commun.; 1997; J80-B-I, pp. 64-73.

5. Kenji, N. The Importance Sampling Simulation of MMPP/D/1 Queueing. IEICE Trans. Fundam. Electron. Commun. Comput. Sci.; 1997; E80-A, pp. 2238-2244.

6. Kobayashi, N. On-line Estimation by Importance Sampling for the Tail Probability of FIFO Queue Length. Proceedings of the 10-th International Workshop on Rare Event Simulation (RESIM ’10); Amsterdam, The Netherlands, 27–29 August 2014.

7. Ngoc, H.N.; Nakagawa, K. Online Estimation for Packet Loss Probability of MMPP/D/1 Queuing by Importance Sampling. Proceedings of the 9th International Symposium on Information and Communication Technology (SoICT ’18); Danang City, Vietnam, 6–7 December 2018; Association for Computing Machinery: New York, NY, USA, 2018; pp. 145-149. [DOI: https://dx.doi.org/10.1145/3287921.3287928]

8. Morio, J. Importance sampling: How to approach the optimal density?. Eur. J. Phys.; 2010; 31, L41. [DOI: https://dx.doi.org/10.1088/0143-0807/31/2/L01]

9. Cottrell, M.; Fort, J.C.; Malgouyres, G. Large deviations and rare events in the study of stochastic algorithms. IEEE Trans. Autom. Control; 1983; 28, pp. 907-920. [DOI: https://dx.doi.org/10.1109/TAC.1983.1103345]

10. Srinivasan, R. Importance Sampling-Applications in Communications and Detection; Springer Science & Business Media: Berlin/Heidelberg, Germany, 2013; Volume 15, XIV, 242. [DOI: https://dx.doi.org/10.1007/978-3-662-05052-1]