Abstract
In queueing models, phase-type servers are very useful since they can be used to approximate any service distributions of real-life telecommunication networks. However, the corresponding Markov chains are in most cases very large and hard to solve without the help of an efficient simulation model. A simulation framework, based on perfect sampling, had already been developed to address the evaluation of large queueing models. Perfect sampling enables us to compute unbiased samples of the stationary distribution based on a backward coupling of the Markov chain trajectories, starting from all possible states. We use an optimized algorithm which consists in computing a set of extremal envelopes which bound all trajectories to obtain the coupling. When the chain is monotone, the envelopes are directly obtained by the trajectories coming from the extremal states. Otherwise, envelopes are upper and lower bounds over the whole set of trajectories.
In this article we prove that phase-type systems are structurally non monotone even with the simplest phase-type distributions like Erlang, hyper-exponential or Coxian-type distributions. Then we provide a generic method for the computation of the upper and lower envelopes for general finite phase-type services in queueing networks. This allows the perfect sampling of such systems and its efficiency is illustrated through several examples.
Keywords: Markov chains, perfect simulation, phase-type distributions.
Introduction
Using queueing models of real-life networks for dimensioning and performance evaluation can be dated back to the pioneering work of Erlang. The famed Erlang formulas have been used for assessing the performance of telecommunication network since them. The simplest formulas are based on the assumption that the incoming traffic is Poisson and the service distribution is exponential.
While the Poisson assumption for the incoming traffic can be justified by the fact that incoming traffic is often the superposition of a large number of independent processes, the exponential distribution of service only comes from the fact that it leads to simple mathematical models.
Phase-type distribution for services is a class of distributions that is both easy to manipulate from a mathematical point of view because it is directly amenable to a Markov chain model and of great generality since this family is dense among all distribution over R+. Of course, when a network designer wants to analyse a given distribution of service (say 7), far from an exponential, the number of phases required for a phase-type distribution to approach 7 increases. The number of states in the corresponding Markov model also increases, making it harder and harder to solve.
In this paper, we investigate how perfect simulation can be used to generate samples from the stationary distribution of a Markov chain that models a network of queues whose service distributions are phase-type. This problem contains two difficulties. The first one comes from the fact that the service events are not monotone in this case. This prevents us from using the classical extreme point simulation algorithm. Instead, we show how envelope simulations, introduced in [1], can be used in this framework. The second difficulty comes from the fact that envelopes may take an considerable amount of time to couple under a certain range of parameters. We have modified the Markov chain (while keeping the same stationary distribution) to make envelopes couple fast in all cases.
Perfect sampling with envelopes for non-monotone chains is presented in Section 2. The model of phase-type server, as well as the construction of the new chain is defined in Section 3. Here we show that a queueing network containing phase-type services is not monotone. We then give the construction of the envelopes that enables us to apply the algorithm previously presented. In section 4, we prove that queueing networks containing phase-type servers do couple under the presented model. Section 5 displays some numerical evidence to study the link between the topology of the graph of phases and the coupling time of the perfect simulation. The influence of the arrival rates as well as of the service rates is also discussed and tested numerically.
Simulation Method
Perfect Sampling
Let
be an irreducible and aperiodic discrete time Markov chain with a finite state space X and a transition matrix
Let
denote the steady state distribution of the chain:
The evolution of the Markov chain can always be described by a stochastic recurrence sequence
an independent and identically distributed sequence of events belonging to a set of events E. The transition function
verifies the property that
for every pair of states
X x X and a random event
Given an initial state
and a sequence of events
the sequence of states
resulting from the stochastic recurrence sequence is called a trajectory of the Markov chain.
Let
denote the function whose output is the state of the chain after n iterations and starting in state
That is:
This notation can be extended to sets of states: for
In the following, |A| denotes the cardinality of set A.
Theorem 1 ([5]). There exists 
such that
The system couples if t = 1. In that case, the value of 
is steady state distributed.
From theorem 1, perfect sampling consists in simulating the MC starting with all states in X. While the end state of trajectories at time 0 is not unique (i.e. while trajectories are not coupled), simulation is started from a further time in the past. Coupling is obtain with any time period greater than the coupling time t*,
The main drawback is that computing the trajectories with all possible initial states is often too costly in practice, because of state-space explosion.
Several approaches have been used to overcome this problem. The main one has been presented in [5]. When the state space X is partially ordered by an order < and when the function &(-, e) is monotone for all e, that is
then the Markov chain is monotone. Then, one only needs to compute trajectories issued from extremal states (maximal and minimal states). When coupling of the extrema occurs, all the trajectories couple by a sandwiching of the extremal trajectories, because of monotonicity.
Monotone perfect sampling is an efficient method to estimate precisely the stationary distribution of a Markov chain. However, as soon as there is one non-monotone case, it becomes unusable. In the next subsection we introduced a new algorithm, initially proposed in [1], which enables the steady-state estimation of a wide class of non-monotone systems.
Bounding Envelopes
In the following, the state space X is endowed with a partial order, denoted With no loss of generality, this partial ordered set can be included into a larger set X, where the order < becomes a complete lattice (i.e. where sup and inf exist). This is typically the case by setting X to be the Dedekind-MacNeille completion of X (i.e. the smallest lattice containing X), see [2].
We consider a new transition function r, called the envelope, defined on the Cartesian product of the completed state space:
for all to, M in X and e in £,
Let us call
the top (resp. bottom) element of X.
The process
is a Markov chain over the state space X x X. However, the upper envelope alone is not a Markov chain (it depends on the lower envelope), neither is the lower one.
Theorem 2 ([1]). Assume that (Mn, mn) hits the diagonal D (i.e. states of the form (x, x)) in finite time:
then Te is a backward coupling time of the initial Markov chain Xn. The state defined by
has the steady-state distribution of the Markov chain.
Algorithm 1. Envelope Perfect Sampling Algorithm
The trajectories of m, M as constructed in Algorithm 1 are the envelopes of the chain. All states of the chain stay in
Therefore, as soon as m = M (after Te steps), the chain has coupled at state m = M that is an unbiased sample of n.
The gain, compared to the classical perfect sampling algorithm, is that the complexity does not depend on the size of the state space. However, the coupling time of envelopes Te might be larger than the coupling time t* for the initial chain. Therefore, the efficiency of envelopes depends on the comparison of Te with t*.
In the framework of queuing networks, this comparison is often in favor of envelopes (see [1]). For example, batch arrivals, negative customers, fork and join nodes all produce events that are non-monotone but for which envelopes are easy to compute and have a short coupling time.
In the next section, we design the envelopes of generalized phase-type servers.
Phase-Type Server
A phase-type distribution, as defined in [4], is described by a Continuous Time Markov Chain H with an absorbing state (called F), where a state of the MC represents the current phase of the server: A client begins its service in an initial phase randomly chosen according to an initial distribution a. Then it goes through the states of H (called phases) until it reaches the absorbing state, finishing its service. The transition parameters of H are given as follows: The time spent in each phase i is exponentially distributed with parameter Tj. The transitions of the MC are given by a transition matrix P where Pj is the probability to move from phase i to phase j so that the total time of an execution of the MC (the time to reach the absorbing state) follows the defined phase-type distribution. An example is displayed in Figure 1(a).
Fig. 1. Phase-type server example
From this definition, we construct another continuous time Markov chain that models the succession of services. We first define an intermediate chain, where the absorbing state F is connected to all initial states using the initial distribution a: the server moves from phase F to an initial phase to treat the next client. This Markov chain has a drawback: the server will spend some time (exponentially distributed) in state F before starting a new service. However, in most queueing models, service of the next customer starts immediately after the end of the previous one.
Fig. 2. Phase-type tandem example
Therefore, to obtain the CTMC describing the succession of services, the state F is removed and all the predecessors of F inherit its transitions. This CTMC is called the phase-type server in the following. The construction of the phase-type server is given in Figure 1(b) for the running example.
A phase-type tandem queue is composed of an input queue q of capacity C and an output queue q’ of capacity C’ where service times of the server in the input queue (server of q) follows a phase-type distribution. We represent the evolution of the system by a continuous time Markov chain {(Xq,Xp,Xq/)n}ngN, where Xq GXq = {0,…, C} is the number of clients in the input queue, Xp GXp = {1,…, F} is the phase into which the current client is treated and Xq> G Xq> = {0,…, C’} is the number of clients in the output queue. The state space X is the Cartesian product of all components;X = Xq x Xp x Xq/.
When a service occurs in the first server, the client which is served moves from the input queue to the output queue, so the state of the three components is modified simultaneously, so that the three components are synchronized.
However, to construct the global Markov chain for the whole system, one problem remains. Indeed, the state of the first server when the first queue is empty (i.e. when Xq = 0) is not defined.
The example displayed in Figure 2 (a phase-type server of 3 phases with an input queue of capacity 2) will be used in the following to illustrate how this synchronization problem can be solved.
A first solution is to aggregate all phases when Xq = 0 into one unique rest state for the server, this model will be called Single Rest State (SRS, Figure 2(b)) in the following. When a customer arrives while the server is in its rest state, the server chooses a new phase according to the initial distribution a.
Another alternative is to keep the same phase-type server dynamics when Xq = 0. This will be called the Multiple Rest State (MRS, Figure 2(c)) model. Here, the idle server changes its phase as if it were serving a customer. To keep the same probabilistic behavior as with the SRS solution, as soon as a customer arrives, the server chooses a new initial phase according to the initial distribution a, independently of its current phase.
Although the SRS solution may look more natural, we will use the MRS model because the runtime of EPSA is exponentially smaller than with the first solution when the arrival rate is smaller than the service rate (see Section 4.1 for further discussions on this).
A complete description of the phase-type tandem is obtained by synchronizing the services of the phase-type server with the arrivals in the output queue. If the output queue is full, we assume that the client is rejected and lost by the system (overflow policy).
To describe the server dynamics, we associate to each transition of the Phase-type server an event. In Figure 2, dashed arrows represent the end of service of a client, as well as a phase transition from a final phase to an initial phase.
If the transition from phase i to phase j is an end of service, the corresponding event is denoted sij, otherwise it is denoted e ij. The rates of events e ij and s ij are given by
We give the transition function of these two classes of events:
As for arrivals (with rate A ), when the input queue is empty, the server must be reset to an initial phase according to a and that must be included into the transition function of arrival events. Thus, there is one arrival event ai for each initial phase i (ai > 0) with rate
The transition function of an arrival event ai is:
Finally, to model the services of clients in the output queue, we define a service event d with rate ^ and transition function
Non Monotonicity of a Phase-Type Server
We propose to show that phase-type tandems are not monotone on an example: a 2 phase Coxian service. The system consists of an M/Cox2/1 queue sequentially composed with an M/M/1 queue (Figure 3). The state space of the system is
Fig. 3. Tandem queueing network with Coxian service
With the construction defined in the previous section, the events of the phase-type tandem are: arrival in phase 1 denoted by ai, phase 1 end of service denoted by e12, phase 1 end of service with skip of the second phase denoted by s11, phase 2 end of service denoted by s21, departure from the output queue denoted by d.
Proposition 1. Let
be a partial order defined on X by:
a partial order on
Then there is no order
such that the tandem network with Cox-2 service be monotone.
Proof. In the following, the notation
means that the states x and y are not comparable.
For every choice of <p, there exists a counterexample that shows that the system is not monotone.
Envelopes for Phase-Type Service Events
The phase transition, service and arrival events we have defined are non-monotone even for a very simple phase-type distribution. The aim of this section is to define the envelope functions of these events.
We define an order < on X which is the product order of totally ordered components (input queue, phase, output queue). The order on Xq and Xq/ is the natural order on integer. On the other hand, we define an order <p on Xp that will ensure the coupling of the envelopes (proved in section 4) as follows. We choose an arbitrary vertex of the phase-type server graph, that we call imax and we construct the order from a covering in-tree of the graph where imax is the root of the in-tree. For example, on the graph of Figure 1(b), we choose phase 6 to be the root and the in-tree represented by bold arrows on the Figure. This tree describes the Hass diagram of a partial order. The order is an arbitrary linear extension of this previous partial order.
For example, from the graph of Figure 1(b), we can choose the order : 1 -<,„
As mentioned before, the product order
Then, we can define the envelope of events eij (moving from phase i to phase j) and sij (end of service of a client in phase i and moving to phase j).Note that we always suppose that
Otherwise, all states
are unchanged by
so the envelope is also unchanged. We define the envelope function of an event eij by: If
then
If
the marginal envelope functions on the phase component
and 
are
Since an event eij is just a phase transition, the marginal envelope function on q and q’ components always keep the value of the upper envelope M and the lower envelope m.
For an event Sj, the envelopes on the phase component are the same as for an event ej. The marginal envelope functions on q and q’ components of an event Sij are:
If
then
In the same way, the envelope of an event ai, an arrival in the input queue with beginning in the initial phase i if xq = 0, is given by:
Notice that the cost of computing envelopes is in 0(1) here. However, this is not sufficient to ensure the efficiency of the method because envelopes may never couple or have a large coupling time. In the next section, we show that envelopes always couple under the proposed model. In section 5, the mean coupling time is studied for several phase-type examples by numerical experiments.
