Computing the Stationary Distribution of Queueing Systems with Random Resource Requirements via Fast Fourier Transform

Abstract

Queueing systems with random resource requirements, in which an arriving customer, in addition to a server, demands a random amount of resources from a shared resource pool, have proved useful to analyze wireless communication networks. The stationary distributions of such queuing systems are expressed in terms of truncated convolution powers of the cumulative distribution function of the resource requirements. Discretization of the cumulative distribution function and the application of the fast Fourier transform are a traditional way of calculating convolutions. We suggest finding truncated convolution powers of the cumulative distribution functions by calculating the convolution powers of the truncated cumulative distribution functions via fast Fourier transform. This radically decreases computational complexity. We introduce the concept of resource load and investigate the accuracy of the proposed method at low and high resource loads. It is shown that the proposed method makes it possible to quickly and accurately calculate truncated convolution powers required for the analysis of queuing systems with random resource requirements.

1. Introduction

When queuing theory is applied to modeling modern information technology systems, one should take into account various system’s features, such as the reliability of a radio channel in wireless communication networks. The latest models from references [1,2] are widely used for the analysis of such networks. In these models, the analytical solution is often cumbersome or difficult to solve and there is a need for approximate [3] or special [4] computational algorithms.

Resource queuing systems, in which an arriving customer, in addition to a server, demands a random amount of resources from a shared resource pool, have proved useful to analyze wireless communication networks. The steady state distributions of such systems are expressed in terms of truncated convolution powers of the cumulative distribution function (CDF) of the resource requirement. The problem of calculating convolution powers of a CDF arises in solving many problems of queuing theory. Its solution is facilitated if the Laplace–Stieltjes transform (LST) of the CDF is known. In this case, the convolution can be found using tables [5,6], or by numerically reversing the LST [7,8]. The method of numerical transform inversion can be successfully applied for the analysis of some resource queueing systems with discrete CDF of resource requirements, such as product-form loss networks [9,10,11].

In practice, the need often arises for calculating convolutions of an exotic continuous CDF with unknown LST [12,13,14]. An alternative method in such cases is the numerical calculation of convolutions by discretization of the initial distribution functions and then using the discrete Fourier transform (DFT) [15]. When computing convolutions numerically, it is convenient to choose some discretization interval $\tau$ and to approximate the CDF of a continuous nonnegative random variable using its values at 0, $\tau$, $2\tau$, $3\tau ,$ .... Ackroyd et al. [16,17] used two very simple approximations $F_1(x)=G\left(\left[\frac{x}{\tau }\right]\tau \right)$ and $F_2(x)=G\left(\left[\frac{x}{\tau }\right]\tau +\tau \right)$, where $[x]$ denotes the integral part of $x$, thus putting upper and lower bounds on CDF $G(x)$. In reference [17] it was shown that the n-fold convolutions of these approximations always satisfy the inequalities $F_1^{(n)}(x)\le G^{(n)}(x)\le F_2^{(n)}(x)$. A more accurate approximation, $F_3(x)=G\left(\left[\frac{x}{\tau }\right]\tau +0.5\tau \right)$, was used in reference [18] to compute the CDF of stationary waiting time in a single server queue with general inter-arrival and service time distributions by means of the fast Fourier transform (FFT). An algorithm for computing the high-order convolution of CDFs representable as the mixture of continuous and discrete components was proposed in reference [19]. Schaller and Temnov [20] analyze an error in the quintile function of the n-fold convolution and applied FFT for the convolution of heavy-tailed distributions. A general-purpose efficient and precise algorithm based on FFT for convolution of two CDFs is considered in reference [21]. In insurance mathematics, recursive schemes to compute compound distributions have been in use for a long time [22]. A comparison of recursion schemes and algorithms based on FFT can be found in references [20,23,24]. Errors of compound distributions computations are investigated in references [25,26].

For a wide range of Markovian resource queueing systems, the stationary distribution of the total amount of occupied resources is given by a truncated compound distribution [27]. In order to obtain the stationary distribution for a resource queueing system of capacity $L$ it is necessary to compute $L$ truncated convolution powers of the resource requirement CDF. The existing computational methods [28,29] calculate some performance metrics for the resource systems but not the entire stationary distribution. In this paper we develop an algorithm for the computation of sequences of truncated convolution powers of continuous CDF and truncated compound distributions via FFT. Instead of truncation of convolution powers we calculate convolution powers of truncated CDFs, which radically decreases computational complexity. We also propose a new approximation of continuous CDF by arithmetic CDF and compare it with $F_3(x)$. In the next section we describe a class of Markovian resource queuing systems to which the results of this paper are applicable. The computation method is detailed in Section 3. Section 4 presents several numerical experiments, in which we compute series of truncated convolutions and analyze the accuracy of the proposed method. The results are discussed in Section 5.

In the remainder of the paper we adopt the following notation: $R_{+}$ denotes the set of nonnegative real numbers, $G^{\ast n}(x)$ denotes the n-fold convolution of CDF $G(x)$, $a\ast b$ denotes the convolution of sequences $a=(a_i)$ and $b=(b_i)$, and ${\delta }_{i,j}=1$ if $i=j$ and ${\delta }_{i,j}=0$ otherwise.

2. Markovian Resource Queueing Systems

Consider a queuing system that can be described by a process $X(t)$ with a finite state space $X$, whose paths are continuous on the right and have limits on the left. Customers arrive and leave the system one by one. We use $0<a_1<a_2<\dots$ to denote the arrival times and $0<d_1<d_2<\dots$ to denote the departure times. We assume that the process $X(t)$ describes the system so that for any of its states $i\in X$ it is possible to determine the number $\nu (i)$ of customers in the system. We denote the maximum number of customers in the system by L and split the state space $X$ into disjoint subsets $X_k=\{i\in X|\nu (i)=k\}$, $k=0,1,\dots ,L.$

In a resource queuing system, an arriving customer requires a certain amount of resources. $V$ denotes the maximum amount of the resource and $r_n$ the resource quantities requested by $n$-th customer. Information about resources occupied at time t is stored as the vector $s(t)=(s_1(t),s_2(t),\dots ,s_k(t))$ of the length $k=\nu (X(t))$, where $s_i(t)$ is the amount of resource occupied by the customer with the number $i$. The customer arriving at time $a_j$ is admitted to the system only if $\nu \left(a_j-0\right)<L$ and the requested resource quantity can be provided, i.e., if $\sigma (a_j-0)+r_j\le V,$ where $\sigma (t)=s_1(t)+s_2(t)+\dots +s_k(t)$ is the total amount of resources occupied at time $t$.

All customers in the system are numbered. If the system accepts the customer arriving at $a_n$, it will be assigned a number ${\phi }_n$ from the interval $1\le {\phi }_n\le \nu (X(a_n-0))+1$. The numbers ${\phi }_n,{\phi }_n+1,\dots ,k,$ assigned to the customers before $a_n$, will be increased by one. If, at the service completion time $d_n$, the system leaves the customer with a number ${\psi }_n$, the numbers ${\psi }_n+1,{\psi }_n+2,\dots ,\nu (X(d_n-0)),$ assigned to the customers before $d_n,$ will be reduced by one. For example, in the “First In First Out” service discipline, the numbers of incoming and serviced customers are determined by ${\phi }_n=\nu (X(a_n-0))+1$ and ${\psi }_n=1$, and for the “Last In First Out” service discipline, we have ${\phi }_n=\nu (X(a_n-0))+1$ and ${\psi }_n=\nu (X(d_n-0))$.

We assume that the process $Y(t)=(X(t),s(t))$ is a homogeneous Markov jump process with a state space $Y=\{(i,x)|i\in X,x\in S_{v(i)}\}$, with $s_k=\{(x_1,x_2,\mathrm{...},x_k)\in R_{+}^k|x_1+x_2+\mathrm{...}+x_k\le V\}$. The sojourn time in each state $(i,x)$ has an exponential distribution with a parameter depending only on the state i of the system. Additionally, the vector of occupied resources $s(t)$ can change only upon customer arrivals and departures. With $p_i$, $i\in X$ we denote the stationary probabilities of the system states in the case of unlimited resources. For the Markovian resource queueing systems considered here, probabilities $p_i$ yield a solution to the steady-state equations with a block tri-diagonal generator [27].

3. A Method for Computing the Stationary Distribution

Consider the stationary probabilities of the process $Y(t)$,

$$P_i=\underset{t\to \infty }{\lim }P\{X(t)=i\},i\in X_0,P_i(x)=\underset{t\to \infty }{\lim }P\{X(t)=i,s(t)\le x\}, i\in X_k, X\in S_k,0<k\le L.$$

(1)

For a wide range of resource queueing systems described in Section 1, the stationary distribution of $Y(t)$ has the product form:

$$P_i=c{p}_i, i\in X_0, P_i\left(x_1,\dots ,x_k\right)=c{p}_{i}G\left(x_1\right)\cdots G\left(x_k\right), i\in X_k,\left(x_1,...,x_k\right)\in S_k,0<k\le L.$$

(2)

where $G(x)$ is the CDF of the customer resource requirements and $c$ is the normalizing constant:

$$c={\left({\displaystyle \sum _{k=0}^L{\displaystyle \sum _{i\in X_k}{p}_i{G}^{\ast k}(V)}}\right)}^{-1}$$

(3)

Necessary and sufficient conditions for the product-form (Equation (2)) of the stationary distribution (Equation (1)) can be found in reference [27].

It follows from Equation (2) that the stationary distribution

$$Q(x)=\underset{t\to \infty }{\lim }P\{\sigma (t)\le x\},0\le x\le V$$

of the total amount of occupied resources is given by the truncated compound distribution:

$$Q(x)=c{\displaystyle \sum _{k=0}^L{q}_k{G}^{\ast k}(x)},0\le x\le V$$

where

$$q_k={\displaystyle \sum _{i\in X_k}{p}_i},0\le k\le L$$

is the stationary distribution of the number of customers in the system.

In order to find the stationary distribution (Equations (2) and (3)) of resource queueing systems, one should compute a series of convolutions $G^{\ast n}(x)$, $n=2,3,\dots ,L,$ in the interval $0\le x\le V$. These functions are related by the following recursion formula:

$$G^{\ast 0}(x)=1,G^{\ast n}(x)={\displaystyle \underset{0}{\overset{x}{\int }}{G}^{\ast (n-1)}(x-t)d}G(t),n=1,2,\dots ,L,0\le x\le V.$$

(4)

For computing $G^{\ast n}(x)$ we approximate $G(x)$ by some step function $F(x)$ with jumps of size $f_k$ at $\tau k$, $k=0,1,2,\dots ,K$, where $\tau K=V$. The convolutions $F^{\ast n}(x)$ are also step functions and their jumps $f_k^{\ast n}$ at $\tau k$ are related by the recursion formula:

$$f_k^{\ast 0}={\delta }_{k,0},f_k^{\ast n}={\displaystyle \sum _{i=0}^k{f}_{k-i}^{\ast (n-1)}}{f}_i,n=1,2,\dots ,L,0\le k\le K.$$

(5)

Thus, the task of computing approximately the integrals in Equation (4) depends on computing the convolution powers $f^{\ast n}=(f_k^{\ast n})$ of the sequence $f=(f_k)$.

3.1. Computing Convolutions via the Discrete Fourier Transform

The method for computing convolutions by means of the forward and inverse discrete Fourier transforms (DFT) is well known [15]. Let $a=(a_i,i=0,1,\dots ,I),$ and $b=(b_j,j=0,1,\dots ,J)$ be two finite sequences. Choose a whole number $M>I+J$ and let $A=(A_m)$ and $B=(B_m)$ denote the DFTs of the sequences $a$ and $b$ over the range from 0 to $M-1$:

$$A_m={\displaystyle \sum _{r=0}^I{a}_r}{e}^{2\pi jmr/M}, B_m={\displaystyle \sum _{r=0}^J{b}_r}{e}^{2\pi jmr/M}, m=0,1,\dots ,M-1,$$

where $j=\sqrt{-1}$. Now the elements $c_k$ of the convolution $c=a\ast b$ are given by the inverse DFT:

$$c_k=\frac{1}{M}{\displaystyle \sum _{m=0}^{M-1}{A}_m{B}_m}{e}^{-2\pi jmk/M}, k=0,1,\dots ,I+J.$$

3.2. Computing the Truncated Convolutions

The procedure described above allows computing consecutively convolutions (Equation (5)). This requires the calculation of the DFTs ${\mathsf{\Phi }}^{(n)}=({\mathsf{\Phi }}_m^{(n)})$ of the convolution $f^{\ast n}$,

$${\mathsf{\Phi }}_m^{(n)}={\displaystyle \sum _{r=0}^{nK}{f}_r^{\ast n}}{e}^{2\pi jmr/M_n},m=0,1,\dots ,M_n-1,$$

over the range of length $M_n>nK,$ which increases with $n$.

Equation (5) can be quickened further in the following way. Note that there is no need to calculate the values of $f_k^{\ast n}$ for $k>K$ since they are not present in Equations (4) and (5). For any given sequence $a=(a_0,a_1,\dots ,a_w)$ of length $w\ge K$, let $Tr(a)=(a_0,a_1,\dots ,a_K)$ denote the sequence of its first K + 1 members. Obviously, probabilities of Equation (2) will remain unchanged if we replace $f^{\ast n}$ with their truncations $g^{(n)}=Tr(f^{\ast n})$ in Equations (2) and (3). It follows from Equation (5) that the truncated sequences $g^{(n)}$ satisfy the recursion relation:

$$g^{(1)}=Tr(f), g^{(n)}=Tr(g^{(n-1)}\ast g^{(1)}), n=2,3,\dots ,L.$$

(6)

There are two options for calculating truncated convolution powers $f^{\ast n}$, $n=2,3,\dots ,L$: calculation of convolution powers using the recursion of Equation (5) with their subsequent truncation, and the calculation of the truncated convolution powers using the recursion of Equation (6). In the first case the forward and inverse DFTs deal with sequences with length $M_n>nK$, $n=2,3,\dots ,L$, while in the second case DFTs deal with sequences of equal size $M_2$. For large L it radically decreases computational time because the total length (L − 1)M₂ of sequences $g^{\left(n\right)}$, $n=2,3,\dots ,L$, is much smaller than the total length $M_2+M_3+\dots +M_L$ of sequences $f^{\ast n}$, $n=2,3,\dots ,L$. Computation time of sequences $g^{(n)}$ can be substantially reduced further using FFT with $M_2$ set equal to a power of 2 [15], for example:

$$M_2=2^m,m=1+[{\log }_2(1+K)].$$

(7)

3.3. Computing the Stationary Distribution

From the above it follows that one can compute the stationary distribution of a resource queueing system via FFT in the following steps:

• Choose a whole number $K$ and a discretization step size $\tau =V/K$.
• Choose a step function $F(x)$ with jumps of size $f_k$ at $\tau k$, $k=0,1,\dots ,K$ that would approximate CDF $G(x)$.
• Apply FFT and compute the sequences $g^{(n)}=(f_0^{\ast n},f_1^{\ast n},\dots ,f_K^{\ast n})$, $n=2,3,\dots ,L.$ using formulae (6).
• Define the functions

  $$F^{\ast n}(x)={\displaystyle \sum _{k=0}^{\left[\frac{x}{\tau }\right]}{f}_k^{\ast n}},0\le x\le V, n=2,3,\dots ,L.$$
• Obtain the stationary state probabilities $p_i$, $i\in X$, of the system with unlimited resources.
• Use distribution $p_i$, $i\in X$, and the approximations $G^{\ast n}(x)\approx F^{\ast n}(x)$ given by Equations (2) and (4) to compute the stationary distribution of the system with limited resources.

4. Numerical Examples

To approximate the resource requirement CDF $G(x)$, it is convenient to use one of the previously mentioned functions (see Figure 1):

$$F_1(x)=G\left(\left[\frac{x}{\tau }\right]\tau \right),F_2(x)=G\left(\left[\frac{x}{\tau }\right]\tau +\tau \right),F_3(x)=G\left(\left[\frac{x}{\tau }\right]\tau +0.5\tau \right),0\le x<V.$$

We use $F_3(x)$ because it approximates $G(x)$ closer than $F_1(x)$ and $F_2(x)$. We also consider a new approximation:

$$F_4(x)=\frac{1}{4\tau }\left(G(k\tau )+2G((k+0.5)\tau )+G((k+1)\tau )\right),k\tau \le x<(k+1)\tau ,k=0,1,\dots ,K-1,F_4(V)=G(V),$$

which, in the interval $k\tau \le x<(k+1)\tau$, equals the average value of $G(x)$ obtained approximately using the trapezoidal rule [30]:

$$\frac{1}{\tau }{\displaystyle \underset{k\tau }{\overset{(k+1)\tau }{\int }}G(t)dt}\approx \frac{1}{4\tau }\left(G(k\tau )+2G((k+0.5)\tau )+G((k+1)\tau )\right).$$

The plot of $F_4(x)$ is similar to $F_3(x)$ and, therefore, not shown in Figure 1.

To verify the accuracy of the proposed method, it is necessary to compare the obtained approximate values of truncated convolution powers with their exact values $G^{\ast n}(x)$, $0\le x<V.$ The gamma distribution

$${\mathsf{\Gamma }}_{\alpha ,\beta }(x)=\frac{1}{\mathsf{\Gamma }(\alpha )}{\displaystyle \underset{0}{\overset{\beta x}{\int }}{e}^{-t}}{t}^{\alpha -1}dt,$$

is well suited for these purposes because it has a useful property facilitating computation of exact values of convolution powers:

$${\mathsf{\Gamma }}_{\alpha ,\beta }^{\ast n}(x)={\mathsf{\Gamma }}_{n\alpha ,\beta }(x).$$

We compare the exact values of $G^{\ast n}(x)={\mathsf{\Gamma }}_{\alpha ,\beta }^{\ast n}(x)$ with their approximations $F_3^{\ast n}(x)$ and $F_4^{\ast n}(x)$ obtained using the method proposed in Section 3.2. The gamma distribution was calculated via the incog procedure from reference [31], and the convolutions were found via FFT using the ffako procedure from reference [32].

Now we consider the case where the total amount of system resources $V=100$ and the interval $[0,V]$ is divided into $K={10}^5$ subintervals. Figure 2 shows the approximation errors

$${\epsilon }_i^{(n)}=\underset{1\le k\le K}{\max }|G^{\ast n}(k\tau )-F_i^{*n}(k\tau )|,i=3,4,$$

as functions of $n$ for various values of the mean $m=\alpha /\beta$ and the square of the coefficient of variation $s=1/\alpha$ of the resource requirements. The time needed to compute the series of truncated convolutions of length 250 on a PC with an Intel i5-7200 2.5 GHz processor did not exceed 20 min. For s < 1, the approximation errors ${\epsilon }_3^{(n)}$ and ${\epsilon }_4^{(n)}$ were under 0.0001 and the difference between them was negligible. Therefore, only results for s > 1 are shown in Figure 2. Computational errors grew with the coefficient of variation and have maximum value mainly for small $n.$

Table 1. Maximum computational errors when calculating the convolutions for $s=3$.

	m	n = 2	n = 3	n = 4	n = 5
${\epsilon }_3^{(n)}$	0.5	0.006027	0.000468	0.000196	0.000154
	1	0.003797	0.000234	0.000097	0.000076
	2	0.002392	0.000117	0.000048	0.000037
${\epsilon }_4^{(n)}$	0.5	0.004003	0.000275	0.000165	0.000126
	1	0.002522	0.000142	0.000084	0.000065
	2	0.001589	0.000073	0.000043	0.000033

shows the maximum computational errors for the cases depicted in Figure 2.

5. Discussion

The levels of light and heavy load in resource queueing systems can vary significantly since they depend on specific for these systems performance metrics. When speaking of resource load, we refer to:

$$R_n=\frac{mn}{V},$$

which allows estimation of the extent to which the system’s resources suffice when it contains $n$ customers. The resource load grows with $n$. For any given resource load $R$ the corresponding maximum system capacity under which the resource load does not exceed $R$ is given by $L=R\frac{V}{m}$.

In the plots of Figure 2, three ranges of $n$ can be clearly distinguished: those of light, moderate, and heavy resource load. Light and high load areas are the areas at the beginning and the end of the n axis respectively, where the approximation error rapidly decreases. The light load area is the area with poorest computational accuracy, especially if the variance of the resource requirements is large (see Table 1). Apparently, this is because for $R_n>1$ the values of $G^{\ast n}(x)$ and $F_i^{\ast n}(x)$ in the interval $[0,V]$ approach 0 as $n$ grows large. In the range of moderate resource load, which is located between light and high load areas, the computational accuracy is mostly high, although it may decrease when approaching the heavy load area. The proposed discretization technique via function $F_4(x)$ yields better or similar results as via function $F_3(x)$. All in all, FFT permits computing long series of truncated convolution powers with sufficient accuracy.