Abstract
This paper investigates the performance of a queueing model with multiple finite queues and a single server. Departures from the queues are synchronised or coupled which means that a service completion leads to a departure in every queue and that service is temporarily interrupted whenever any of the queues is empty. We focus on the numerical analysis of this queueing model in a Markovian setting: the arrivals in the different queues constitute Poisson processes and the service times are exponentially distributed. Taking into account the state space explosion problem associated with multidimensional Markov processes, we calculate the terms in the series expansion in the service rate of the stationary distribution of the Markov chain as well as various performance measures when the system is (i) overloaded and (ii) under intermediate load. Our numerical results reveal that, by calculating the series expansions of performance measures around a few service rates, we get accurate estimates of various performance measures once the load is above 40% to 50%.
1. Introduction
Numerical methods for queueing systems involving multiple queues like queueing networks [1], polling systems [2], priority queues [3], and fork-join queues [4] often suffer from the state space explosion problem. State space explosion refers to the problem that multidimensionality of Markov processes leads to processes with a very large state space. Indeed, the size of the state space of a multidimensional Markov process is the product of the number of states in each of its dimensions. Once a few dimensions are involved, the state space becomes very large and direct solution techniques for Markov processes fail. For some particular types of Markov processes, a solution can be readily found, but this depends on structural properties of the Markov chain at hand. We mention Markov chains with product form solutions (like Jackson networks) [5] and -type and -type Markov processes [6] as particular examples. However many queueing problems do not possess these structural properties, thereby requiring nonstandard solution techniques.
This is the case for the queueing system investigated in this paper. We consider a queueing system with queues in parallel as depicted in Figure 1. Customers in all queues receive service simultaneously and there is a departure from every queue upon service completion. Moreover, whenever one of the queues is empty, the server remains idle. That is, an empty queue completely blocks service for all other queues. This queueing system is a natural abstraction for an assembly operation with in-house production. The queues represent inventories for semifinished products which are replenished by in-house production facilities. The final assembly requires all semifinished products and therefore the assembly operation is halted once any of the inventories is completely depleted. Finally, the service time of the coupled queueing system represents the assembly time.
We study the service-coupled queueing system under Markovian assumptions. That is, we assume independent Poisson arrivals to all queues with arrival rates , respectively, and independent exponentially distributed service times with rate . Even for these simplified assumptions, the analysis of the coupled queueing system is challenging. First, one cannot impose the often simplifying assumption that queues have infinite capacity as the resulting Markov process is either null recurrent if all arrival rates are equal or transient if this is not the case; see [7] for the coupled queueing system with only two queues. Secondly, the state space of the Markov process for the system with queues of capacity is such that a direct solution of the Markov chain is not numerically feasible for moderate and . Finally, matrix-analytic methods for neither -type nor -type queueing systems apply, and there is no product form solution.
To overcome these challenges, literature proposes two alternative approaches, both focusing on approximations for various performance measures of the coupled queueing system. The first approach aims at decomposing the queueing system into a number of independent queueing systems which can be analysed in isolation [8]. Such an analysis approximates the interaction between the different queues by a simpler process which in turn facilitates the analysis. The interaction process is parametrised such that the simplified interaction process corresponds to the expected interaction by the queue in isolation. Alternatively, the system can be studied approximately by means of series expansion techniques if one limits the study to a subset of the parameter space. This is the case in [9, 10] where the coupled queueing system was studied in overload. In these papers it was shown that the terms of the Maclaurin series expansion of the steady-state distribution in the service rate can be obtained at low computational cost. The series expansion of the performance measures can then be easily obtained from the calculated steady-state distribution. However, the numerical approach advocated there only leads to good results when the service rate is close to or, equivalently, when the system is considerably overloaded.
Series expansion techniques for Markov chains go by different names in literature, including perturbation techniques, the power series algorithm, and light-traffic approximations. While the naming is not absolute, perturbation methods are mainly motivated by sensitivity analysis of the results with respect to some system parameter. In particular singular perturbations where the perturbation does not preserve the class-structure of the nonperturbed chain have received considerable attention in literature [11–13]. The power series algorithm transforms a Markov chain of interest in a set of Markov chains parametrised by a variable . For , not only is the chain easily solved, but one can also obtain the series expansion in . For one gets the original Markov chain such that the series expansion can be used to approximate the solution of the original Markov chain, provided the convergence region of the series expansion includes [14–17]. Finally, light-traffic approximation often corresponds to a series expansion in the arrival rate at a queue. For an overview on the technique of series expansions in stochastic systems, we further refer the reader to the surveys in [18, 19].
The present contribution builds on the results of [9, 10] but considers the service-coupled queueing system when the load of the system is lower. In the context of assembly systems, the overload situation is only natural if assembly is the bottleneck in the production/assembly system. In case production is the actual bottleneck, the assembly queues are not overloaded and the results of [9, 10] do not apply. However, it is still worth investigating the assembly system in this case as assembly will be interrupted more often due to a lack of semifinished products.
Balancing computational cost and accuracy, we investigate the use of Taylor series expansions to calculate the performance measures for a wider range of the service rate. In contrast to the Maclaurin series expansions in [9, 10], the terms in the Taylor series expansion around some service rate cannot be obtained directly. Therefore we rely on iterative solution methods to solve for the terms in the Taylor series expansion. So, in contrast to the power series algorithm, our approach does not primarily aim for simplifying the solution of the Markov chain but aims for obtaining the solution in a wide subset of the parameter space at once and relies on iterative procedures to do so.
For any iterative method, a good initial guess of the solution can reduce the number of required iterations considerably. In the present setting, such an initial guess is available if one considers a sequence of Taylor series expansions around increasing values of the service rate starting at . As shown in [9, 10], the initial series expansion around can be calculated efficiently. For higher , the expansion around the preceding -value can be used to get an initial guess.
The remainder of this paper is organised as follows. The model at hand and the numerical evaluation method are described in the next section. We then illustrate our approach by numerical examples in Section 3, prior to drawing conclusions in Section 4.
2. Performance Analysis
We consider a queueing system with finite capacity queues as depicted in Figure 1. We denote the capacity of the th queue by . The arrival process to the th queue is assumed to be a Poisson process with a fixed rate , the arrival processes to the different queues being mutually independent. As mentioned above, service is coupled. This means that there are simultaneous departures from all queues with rate as long as all queues are nonempty, while there are no departures when any of the queues is empty.
In view of the Markovian assumptions on both arrival and service processes, the state (in the Markovian sense) of the queueing system is completely described by the numbers of customers in the different queues. That is, the state of the system is described by a vector , where denotes the number of customers in the th queue and where is the state space. We have the following state transitions from state :(i)Arrival in queue (for ): when , the arrival rate in queue is , the new state being . Here is a vector of zeroes, apart from its th element which is one. There are no arrivals in queue when .(ii)Departure: when all queues are nonempty () there is a departure from all queues with rate . The new state is , where is a vector of ones.
Given the summary of the possible transitions above, the balance equations of the Markov process are readily retrieved. For , let be the steady-state probability vector of the queueing system. Equating the total probability flow out of and into state , we then have the following set of balance equations, for , where denotes the indicator function of the event and where we have assumed for to simplify notation. Since already for a moderate number of queues the state space is prohibitively large to compute the stationary distribution directly, we rely on a series expansion approach in the remainder.
As the system of (1) is finite, we find by Cramer’s rule that the stationary probabilities can be expressed as rational functions of with at most distinct poles and no other singularities. Here is the size of the state space . Denoting the set of singularities by , this observation implies that, for any , the Taylor series expansion in of around converges to the correct value in a neighbourhood of . For further reference, let be the th term in the Taylor series expansion in of around . Hence, in a neighbourhood of , we have
First, when is close to , we approximate the stationary probabilities by their Maclaurin series expansion in as investigated in [9]. Plugging the expansion (2) for in the balance equations (1) and comparing terms in equal powers of , we obtain for and . For , we find by the normalisation condition for . For and , we further find which shows that for (by evaluation of the expression in lexicographical order). The normalisation condition then further yields , such that for . The th-order terms are trivial and the higher order terms can be calculated one by one in lexicographical order of by expressions (3) and (4) above. The numerical complexity of finding the terms of a single order for all is at most. However, one easily verifies that for all lexicographically smaller than , which further reduces the computational complexity of finding the th order terms to . Note that, for large , is considerably smaller than .
While the terms in the Maclaurin series expansion can be calculated efficiently, the resulting expansion only converges to the exact solution in a neighbourhood of as, in general, the region of convergence of the series expansion will be finite. Therefore, we now consider Taylor series expansions around to get results for a wider range of the service rate.
Plugging the series expansion (2) in the balance equations (1) and isolating terms in , we get, for ,for , whereas the normalisation condition yields
In contrast to the Maclaurin expansion above, the system of equations (7)-(8) cannot be solved easily. Therefore, we rely on iterative solution methods to find the solution of this system of equations. More specifically, we use weighted Jacobi iteration which calculates the terms in the series expansion by iteratively evaluating Here denotes the weight of the weighted Jacobi iteration. For each term and , we evaluate for and approximate by for sufficiently large. In practice, we stop iterating when the corresponding terms in the series expansion of the mean and second-order moment of the queue content (cf. infra) converge (up to 6 to 8 significant digits).
This iterative approach is computationally feasible as the number of possible transitions from a state is far less than the number of states (the generator matrix is sparse). More precisely, the number of transitions is related to the number of queues such that the numerical complexity of a single iteration for finding the th order terms for all is .
If is within the radius of convergence of the preceding expansion, say around , we use the preceding expansion to get a first approximation for as to reduce the number of iterations till convergence. That is, we choose If is not within the radius of convergence of the preceding expansion, we set with and with That is, we approximate the coupled queueing system, by a queueing system with independent queues with service rate , where is a crude approximation for the probability that at least one queue is empty. In addition, we set for .
Once the terms in the series expansion are found, we can find approximations for various performance measures. For instance, the th order expansion of the th moment of the queue content is calculated as where denotes the queue content of the th queue and with . In particular, the mean and variance of can be approximated as Note that the above approximation for the variance is not the th order series expansion of the variance as the approximation of the square of the mean also contains terms in for . By numerical experimentation, we found that these higher order terms hardly influence the results.
Analogously, let the system content be defined as the total number of customers in all queues; then we can approximate the mean and variance of the system content as Again the same remark applies to the approximation of the variance.
The effective load is defined as the fraction of time that the server is serving. As the server is serving whenever all queues are nonempty, we find the following th order expansion of the effective load ,
Finally, let the blocking probability be the fraction of customers that cannot be accepted upon arrival in the queueing system. The effective load allows for calculating the blocking probability in the th queue. Indeed, noting that all accepted customers must be served, we have or, equivalently, Notice that only depends on the queue capacity through . The latter is influenced by the capacities of all the different queues, which particularly implies that the capacity of one queue influences the blocking probabilities of the other queues.
3. Numerical Results
We now evaluate our numerical approximation approach by some numerical examples. We focus on the mean and standard deviation of the queue content as well as the blocking probability. Noting that, in a coupled queueing system with nonequal arrival loads, the performance is mainly determined by the queues with the lowest loads (the queues with higher load can be neglected when studying the overall performance), we first focus on a coupled queueing system with an equal arrival rate in all queues. Without loss of generality, we set (as we can scale to investigate a different ). We consider queues, each having capacity .
Figures 2, 3, and 4 depict the mean queue content versus the service rate , the blocking probability versus the service rate , and the standard deviation of the queue content versus , respectively. Note that we have the same blocking probability and the same mean and variance of the queue content for every queue due to symmetry and that we approximate the standard deviation of the queue content by with given in (14). Each figure shows the th-, th-, and th-order approximation on a separate subfigure, and we combine the Maclaurin expansion around with the approximation around for all performance measures. For visual reference, the point is marked on all the figures with a cross. The order of the expansion refers to both the order of the expansion around and the order of the expansion around . In addition, we show simulation results for the performance measures at hand, which allows for evaluating the accuracy of the approximations. We used uniformization to simulate the queueing system (based on the balance equations) and generated samples, for each simulation point. We calculated the confidence interval by means of the batch means method but omitted the confidence intervals from the plots as the obtained upper and lower bounds are visually indiscernible.
For the coupled queueing system under study with queues of capacity , the Markov chain has states. The figures show that the approximations of the mean queue content and the blocking probability are already fairly accurate for the th-order expansion (), whereas the standard deviation of the queue content requires some more terms (). As the order of the expansions further increases, the approximations even more closely approximate the performance measures at hand. The figures further reveal that the match is very good in a limited region (of or of ), while the approximations quickly grow to very large values outside this region. This is not unexpected as the region of convergence is finite for sure ( is a rational function of , cf. supra). While the sharp deterioration of the approximation prevents one to extend the results outside the region of convergence of the series expansion, it does give a clear indication where the approximation is accurate. Overall, we find that the th-order approximations for the mean queue content, the blocking probability, and the standard deviation are accurate for loads above ( below ).
The effect of increasing on the mean queue content and on the blocking probability confirms intuition. If the service speed increases, the mean content decreases and as it is less likely that the queues are full, the blocking probability decreases as well. The decrease is fast for low and slower for larger , the change of the decay rate being around (or a load of ) for the blocking probability and just above for the mean queue content. For the standard deviation, we observe that it increases with .
Next, we study an example with nonequal arrival rates at the different queues. In particular, we consider a system with queues, each having capacity , which results in a Markov chain with states. In order to investigate the impact of nonequal arrival loads, we consider a system with two arrival rates: arrival rate for half of the queues and arrival rate for the remaining queues.
Figures 5, 6, and 7 depict the mean queue content, the blocking probability, and the standard deviation of the queue content versus the service rate , respectively, for queues with arrival rate as well as queues with arrival rate . We again depict approximations of order , , and on different subfigures. For every order , we consider the expansion around and the expansion around , the point being marked with a cross on all plots. The plots again reveal that the approximations are quite accurate, especially the 30th-order approximation which is again accurate for up to . An increase of leads to a decrease of the mean queue content and of the blocking probability as for the symmetric case, while it leads to an increase of the standard deviation of the queue content. Also, the queues with the highest arrival rate () have higher mean queue content and blocking probability as there are more arrivals, which also leads to a reduction of the standard deviation of the queue content, as the more heavily loaded queue is close to full most of the time.
As a final example, we assess the impact of the number of queues involved. To this end, we compare the performance of the queueing system with , , and queues. All queues have capacity and equal arrival rate . Figure 8 shows the 30th-order approximations (in and ) for the mean queue content and the blocking probability as a function of the service rate . We can readily observe that adding queues lead to performance degradation (higher mean queue content and higher blocking probability), especially when the system is not in overload. This is not unexpected as it is more likely that one of the queues is empty in systems with more queues. For coupled queueing systems in overload, the number of queues involved has hardly any impact on performance though. In overload, it is unlikely that queues are empty, so the number of queues does not matter.
4. Conclusions
In this paper we presented a numerical approach for the performance evaluation of coupled queueing systems. The study was motivated by an assembly-like system, where inventory replenishments can be modelled by Poisson processes. The presented method focuses on coupled queueing systems working under intermediate load and builds on a previously designed method for such systems in overload. We showed that the region where an accurate estimation is obtained can be extended to lower loads by iteratively calculating the terms of the Taylor series expansion of the steady-state probability vector.
An important contribution of the study is that the problem is tackled numerically, while existing analysis methods for large-scale queueing systems mainly rely on simulation. We showed that our analysis method allows for performance evaluation under intermediate load, although the specific region of accuracy may vary depending on the system size and structure.
Conflicts of Interest
The authors declare that they have no conflicts of interest.