{'fraction_non_alphanumeric': 0.05403647597626081, 'fraction_numerical': 0.011950895639684562, 'mean_word_length': 3.7906010645203168, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 13, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Efficiency of simulation in monotone hyper-stable queueing networks

We consider Jackson queueing networks with finite buffer constraints (JQN) and analyze the efficiency of sampling from their stationary distribution. In the context of exact sampling, the monotonicity structure of JQNs ensures that such efficiency is of the order of the coupling time (or meeting time) of two extremal sample paths. In the context of approximate sampling, it is given by the mixing time. Under a condition on the drift of the stochastic process underlying a JQN, which we call hyper-stability, in our main result we show that the coupling time is polynomial in both the number of queues and buffer sizes. Then, we use this result to show that the mixing time of JQNs behaves similarly up to a given precision threshold. Our proof relies on a recursive formula relating the coupling times of trajectories that start from network states having “distance one”, and it can be used to analyze the coupling and mixing times of other Markovian networks, provided that they are monotone. An illustrative example is shown in the context of JQNs with blocking mechanisms.


Introduction
The stationary behavior of Markovian queueing networks (QN) can be computed quite efficiently only under specific assumptions that yield the so called productform property [5,11]. This property means that the stationary probability distribution of network states can be written, up to a normalizing constant, as the product of a number of simple terms, where each term is associated with a different network node (or queue), and provides a way to compute the stationary behavior of a QN that is much faster than the solution of the global-balance equations of the underlying Markov chain [6]. In several cases, product-form QNs have a restricted modeling power because they often assume that nodes have infinite buffer sizes or that the behavior of a network node does not depend on the state of other nodes; e.g., [10]. Examples of phenomena that do not yield, in general, the product-form property are loss or blocking mechanisms due to finite-buffer constraints or state-dependent routing. On the other hand, the stationary behavior of non-product-form QNs is extremely difficult to obtain. While it is possible to obtain it through the solution of a set of linear equations, i.e., the global-balance equations of the underlying Markov chain, the huge size of their state space makes this approach prohibitively expensive from a computational standpoint. For instance, for a QN with 10 nodes where each node has a buffer of 10 units, the number of such equations is not smaller than 11 10 . In this setting, simulation is a useful approach to obtain robust measures and insights on the stationary behavior.
Existing research in the simulation of Markov chains relies on two approaches: simulation into the future (or Monte Carlo simulation) and simulation from the past (or coupling from the past-CFTP). These are briefly summarized in the following.
Simulation into the future generates a trajectory from one state of the chain according to its transition matrix until it is believed that the proportion of visits on any network state is sufficiently close to its corresponding stationary probability. The mixing time [15], i.e., the point in time where the Markov chain is close to its stationary behavior up to a given precision threshold, is a desirable quantity of interest for deciding when to stop the simulation, and it is used to measure the efficiency of simulation into the future. In practice, this type of simulation has the drawbacks of producing excessively-long simulations and approximate estimates of stationary measures.
On the other hand, simulation from the past (or CFTP) generates a trajectory for each state of the chain according to its transition matrix and iteratively goes backward in time until all trajectories have collided in a singleton [19]. The remarkable property of this type of simulation is that such a singleton is perfectly distributed according to the stationary distribution of the chain. Then, to infer the stationary behavior of the chain, one can get several independent samples by running several such simulations and use the law of large numbers. It is important to note that running a trajectory from each state of the chain sounds impractical. However, in monotone chains, which is the case of the QNs studied in this article, it suffices to generate a trajectory only for each extremal state, and there will be only two of such states in our case. This holds true because monotonicity ensures that the trajectories from non-extremal states are sandwiched between the ones from extremal states.
The major issue behind CFTP is the understanding of the time it takes to produce one sample, which gives a quantitative estimate of the amount of resources which is needed to simulate the Markov chain of interest. In the context of monotone chains, this time is of the order of the coupling time (also known as meeting time), i.e., the point in time where the trajectories starting from all possible states collapse into a single one [19].
The mixing time and the coupling time of a Markov chain are related (see Sect. 4). In this article, we are interested in bounding these two quantities for a class of important QNs.

Related work
Since the structure of QNs is arbitrary, exact expressions for the mixing and coupling times are typically difficult to obtain and, therefore, one seeks bounds. In some cases, the mixing time of Markovian QNs with finite state-space can be bounded through general results that involve the second largest eigenvalue modulus of the transition matrix of the chain (see [7,Theorem 3.2]). However, this numerical approach usually hides the understanding of its qualitative dependence with respect to general input parameters.
Qualitative bounds on both the coupling and mixing time of closed Jackson QNs with infinite buffer sizes [6] have been proposed in [8,13,14], where the product-form of such networks is the essential tool used for their derivation; see also the references therein for other results concerning product-form QNs. However, the problem becomes much more difficult when buffer sizes are finite because in general the product-form property is lost. In this context, open Jackson QNs with finite buffers (or JQN in the following) have been studied during the last decade. Upper bounds on the mean coupling times of such JQNs are derived in [2], where conditions are given to prove that (i) the coupling time is linear in the total number of queues and exponential in the size of the largest buffer size, and (ii) the coupling time is linear in the sum of the buffer sizes of all queues, provided that the total number of queues is constant.
Linearity in the buffer sizes holds true if the JQN is acyclic, as has also been shown in [9] in the context of state-dependent routing, or if some other strong conditions on the JQN parameters are satisfied. These conditions, which include restrictions on network topology, will be discussed in Sect. 3.2. It is conjectured in [2,20] that linearity in the buffer sizes holds even when the network contains cycles. As numerical experiments reveal, this indicates that the mean coupling time of Markovian QNs can be much smaller than the size of the state space, which is given by the product of the buffer size of each queue, implying that CFTP is efficient for simulating the behavior of JQNs.

Our contribution
In this article, we consider JQNs where jobs that try to join a saturated queue are lost. These QNs are intractable in general, as the only exact solution method that is known in the literature relies on the numerical solution of the global-balance equations of the underlying Markov chain. This amounts to solving M i=1 (C i + 1) linear equations (one for each state), where M is the number of queues and C i is the buffer size of queue i, for i = 1, . . . , M. As a matter of fact, the stationary distribution of this class of QNs does not have a product-form; e.g., [2,9]. To keep the product-form property in networks of queues with finite buffers, one needs special conditions on the topology, the value of the parameters and on the blocking policies used to deal with saturated queues (see for example [4] for a rather exhaustive treatment of such questions). In the general case treated here, no product-form solution exists and the stationary behavior is prohibitively expensive to compute. This motivates our investigation of the efficiency of simulation.
We study the mean coupling time of JQNs in a regime where both the number of queues and buffer sizes vary. In our main result, we give a sufficient condition to prove that the coupling time of JQNs is polynomial in M and the C i 's. This condition lets us deal with networks having arbitrary topology, for which the best bound available in the literature is exponential in the C i 's, provided that both M and the C i 's vary [2]. Exploiting classical arguments, we then use this result to prove that the mixing time of JQNs behaves similarly up to a precision threshold.
The starting idea of our proof relies on a new recursive formula on the coupling times of special trajectories having "distance one" at most and follows by stochastic comparison arguments. This allows for a tractable decomposition of the problem that yields our bound. Our approach also provides a new framework for obtaining bounds on the coupling and mixing times of other monotone queueing systems with arbitrary topology, and an example is given in the context of JQNs with blocking [4].
The article is organized as follows. In Sect. 2, we give a probabilistic description of JQNs in terms of job movements and an equivalent definition in terms of discrete events. We also briefly review the CFTP method to provide the necessary background. In Sect. 3, we derive our polynomial bound on the coupling time, and, in Sect. 4, we exploit this bound to derive a bound on the mixing time. Section 5 shows how our approach applies to other queueing network models, and Sect. 6 draws the conclusions of our work. A preliminary version of this study appeared in [3].

Queueing network model
We consider JQNs with M queues. The vector C = (C 1 , . . . , C M ) denotes the buffer size of each queue. If not otherwise specified, indices i, j, k will implicitly range from 1 to M.
An infinite stream of jobs that follow a Poissonian process with rate λ joins the JQN from an external source. The probability that a job joins queue i, upon arrival to the network, is p 0i . Thus, i p 0i > 0. In queue i, each job requires some processing for an exponentially distributed amount of time with mean service rate μ i . The service discipline of each queue i is work-conserving. Upon completion of service at queue i, a job is sent to queue j with probability p i j , and it is accepted if queue j has an available slot (i.e., if it is non-saturated), otherwise it is lost. The probability that a job leaves the network after service at i is p i0 , which can be also interpreted as the probability that a job is sent to a queue with buffer of size zero. Since each job eventually leaves the network, i p i0 > 0. The JQN with M queues described above can be seen as a discrete-event system with a single type of events, namely a i j , (i, j ∈ {0, 1, . . . , M}) corresponding to the service of one job in queue i that then joins queue j. The dummy queue 0 corresponds to the outside world. An event of type a 0 j is an exogenous arrival in queue j and an event of type a i0 corresponds to the departure of a job from queue i. If queue i is empty or if queue j is full, then event a i j is disabled. The set of all events is denoted by A .
The rate of event a i j is γ i j and, in view of Remark 1, it is independent of M and C, for any i, j. Using the previous description of a JQN, if i, j = 0, then γ i j = μ i p i j , γ 0 j = λp 0 j , and γ i0 = μ i p i0 . The total event rate Γ def = i≥0, j≥0 γ i j is finite (we fix γ 00 = 0).
We also denote by α i the expected number of events occurring between two consecutive events involving queue i. By definition, we have The continuous-time Markov chain described above can be transformed into a discrete-time Markov chain x n with the same stationary distribution by uniformization by constant Γ . In the following, this discrete chain is assumed to be irreducible and aperiodic. The evolution of the Markov chain x n can be written under the form x n+1 = φ(x n , u n ) where u n is a random variable over the event space that takes value a i j with probability γ i j /Γ , and the transition function φ : S × A → S is defined as follows: Let φ (n) : S × A n → S denote the function whose output is the state of the chain after n iterations starting in state x ∈ S . That is: This notation can be extended to sets of states: for E ⊆ S , In the following, |E| denotes the cardinality of set E.
By definition, the function φ is monotone for all event a i j ∈ A . This implies that the trajectories of the Markov chain starting from ordered initial states (using the component-wise ordering) stay ordered: if x 0 ≤ x (1) 0 then is distributed according to its stationary distribution as soon as it is reduced to a singleton, and the mean backward coupling time, defined Theorem 1 has an algorithmic counterpart, namely a coupling from the past (CFTP) algorithm, given as Algorithm 1. In this algorithm, it is implicit that the events u −i are stored and reused along the iterations of the repeat cycle.
Provided that n is initialized properly, the average time efficiency of Algorithm 1 is O(ET ) where T is the number of iterations in the last for loop (called couplingtime in the following); see [19].
The notation used in the remainder of the article is summarized in Table 1 for quick reference.

Bound on coupling time
In this section, we present our main result (Theorem 2) on the efficiency of CFTP applied to the queueing network model described above. Let be the coupling time of x and y, i.e., the random variable of the time where the trajectories starting in states x ∈ S and y ∈ S of a JQN with M queues and buffer sizes C meet for the first time.
Probability that a job is forwarded to j upon completion at i, Probability that a job joins queue i upon arrival to the network, Probability that a job leaves the network after completion at i, μ i Mean job service rate at queue i, Coupling time of the trajectories starting in states x and y (see (2)) The mean coupling time Eτ (0, C) is related to the efficiency of simulation from the past. In particular, it is of the order of the average time efficiency of Algorithm 1, as discussed above. The aim of this article is to bound Eτ (0, C) from above with respect to JQNs with arbitrary topology.
In the following, we denote by ≤ st the usual stochastic order [17]: for two random variables X 1 and X 2 , X 1 ≤ st X 2 if and only if Pr(X 1 ≥ x) ≤ Pr(X 2 ≥ x). The following proposition, known as Strassen's theorem, gives an equivalent characterization of the ≤ st -order and will be the basic tool used in our proofs.

Proposition 1
For two random variables X 1 and X 2 , X 1 ≤ st X 2 if and only if there exists a pair of random variables (X 1 ,X 2 ) defined on a common probability space such that X 1 = distX1 , X 2 = distX2 , and Pr(X 1 ≤X 2 ) = 1.
The following proposition is the starting point of our proof technique.
Proposition 2 For all i and y ∈ S : y i > 0 where Proof (of Proposition 2). First, we introduce the following lemma.
Lemma 1 For any possible event a ∈ A , k = 1, . . . , M, and any state y such that y, y − e k ∈ S , Proof In fact, we have the following cases for every y ∈ S and k such that y−e k ∈ S : under the same sequence of events), for some k (see Fig. 1 for the latter case where k = i and y = C). Our upper bound on τ (0, y) follows by excluding the possibility of the former case. Taking into account only the latter case, the coupling of the trajectories from y and 0 occurs before τ (0, y − e i ) plus the additional coupling time of the trajectories starting in x * and x * + e k , where x * , x * + e k ∈ S . Taking the worst-case k, this establishes the stochastic comparison in inequality (3).
Proposition 2 provides a recursive framework to analyze τ (0, C) in terms of the coupling times of trajectories having "distance one" at most, i.e., τ (x * , x * + e k ). The fact that the trajectories from x * and x * + e k remain at distance one at each time step follows by Lemma 1, used in the proof of Proposition 2. As we show in the following, this decomposition of the problem is very useful because these trajectories appear easier to analyze than directly τ (0, C).
When a single queue is considered, i.e., M = 1, the next proposition is known for the mean hitting time from state x to state zero, defined as E min{n : φ (n) (x, u 1→n ) = 0} [9].
For an M/M/1/C queue, it is clear that the mean hitting time from C to 0 bounds from above its mean coupling time of the trajectories from states C and 0. This observation, together with previous proposition, will be exploited in Proposition 4.
Proof Relation (6) follows immediately by coupling the sample paths of both processes (φ (n) (x, u 1→n )) n∈N and (φ (n) ∞ (x, u 1→n )) n∈N under the same sequence of events. Then, using (6) and the list of transitions in the proof of Lemma 1, one can easily verify that if φ (n) u 1→n ). This proves (7).
Let Γ i be an upper bound on the maximum mean rate in which jobs can arrive in queue i. For any network state, the minimal upper bound that one can choose is which corresponds to the situation where all queues are non-empty. Let also The following proposition provides a bound on the mean coupling time of the trajectories that start from states x and x − e k .

Proposition 4
Let a JQN be given such that ρ i < 1, ∀i. Then, for any x ∈ S and k 1 such that x − e k 1 ∈ S , and (V k 1 (1), . . . , V k 1 (M)) is the unique solution of the linear system Therefore, for n > 0, we define and , 1} for all n, and therefore I (n) is the index i of the unique queue such that φ (n) ∞ (x − e k , u 1→n ) i = 1 if such index exists, and zero otherwise. The process (I (n)) n∈N starts in k 1 and remains in k 1 until when one of cases ii).c and iii).b occur (these cases are defined in the proof of Lemma 1). For instance, case iii).b corresponds to a routing event a k 1 k 2 , for some k 2 > 0 such that ∞ (x, a 1→n ) k 2 < C k 2 = ∞. If case iii).b occurs at time n, then I (n + 1) = k 2 . If ii).c occurs at time n, then I (n + 1) = 0, and it will remain to zero for all n > n + 1 because both trajectories from x and x − e k 1 met. Therefore, if H is the random variable of the total number of jumps of the process (I (n)) n∈N , then the ergodicity of the Markov chain φ

E(T (h)|H = h)Pr(H = h)
We now derive a bound for

Lemma 3 For t ≤ h,
Proof Between times T (t − 1) and T (t) − 1, I (n) is constant, which means that there exists i such that φ (n) Assume t = 1. Since T (0) = 0 and using that ρ i < 1 for all i, it is clear that E(T (1)|H = h) can be upper bounded, up to the multiplicative constant α k i (see (1)), by the mean hitting time from x k 1 to 0 of an M/M/1 queue with arrival rate Γ k 1 and service rate μ k 1 . This holds true by simply coupling such M/M/1 queue with the one-dimensional process (φ (n) ∞ (x, u 1→n )) n∈N moves along dimension k 1 . The stochastic dominance of such M/M/1 queue holds because Γ k 1 is the maximum rate in which jobs can join queue k 1 . The hitting time to zero corresponds to the occurrence of one of cases ii).c and iii).b, which trigger a jump of I (n). Therefore, using also Proposition 3, we obtain Assume t > 1. Using the same coupling argument above, E(T (t)−T (t −1)|H = h) can be upper bounded by the mean hitting time from some initial state explicitely expressed in the following to 0 of an M/M/1 queue with arrival rate Γ k t and service rate μ k t , provided that k t = I (n) when T (t − 1) ≤ n < T (t). The initial state of such M/M/1 queue must be greater than or equal to the random variable of the value of φ (n) ∞ (x, u 1→n ) k t at the time n = T (t − 1) of the (t − 1)-th jump of (I (n)) n∈N . We denote such random variable by J (t − 1). Conditioning on J (t − 1) and using again Proposition 3, we have The following lemma provides a bound on EJ (t − 1) and concludes the proof of (21).

Lemma 4 EJ
Proof For all t, we recall that we defined k t as the value of I (n) when T (t − 1) ≤ n < T (t). First, we observe that φ  [16] and follows again by coupling the sample paths of both vector processes under the same sequence of events. Therefore, for any time n (random or not), It is well known, e.g., [1], that for the stable M/M/1 queue Y i (n), for any non-random time n. Using (25) in the former, we get for any non-random time n. The proof is thus concluded if we show that (27) holds true even when n is the random time T (t − 1), which is what we prove now. We recall that In (31), we have used that T (t − 1) and Y k t (n) are independent random variables: in fact, by construction, T (t − 1) is a random variable that only depends on some queue k t−1 = k t that hits zero (because either case ii).c or iii).b occurs at time T (t − 1) − 1) and all the M queues Y i (n)'s are independent.
Using (21) in (20), we obtain Therefore, it suffices to understand the magnitude of EH = h h Pr(H = h). Now, let I (t) be the value of (I (n)) n∈N at the time where it makes the t-th jump. Since the jumps of (I (n)) n∈N are only due to occurrences of cases ii).c and iii).b, ifĨ (t) = j then the probability thatĨ (t +1) = i is , for all j, i. Given these jump probabilities, we can interpret the stochastic process (Ĩ (t)) t∈N as a discrete-time Markov chain where state zero is absorbing. Given this interpretation, we let V k 1 (i) denote the mean number of visits that (Ĩ (t)) t∈N makes to state i when it is started in k 1 . This means EH = M i=0 V k 1 (i). From the theory of absorbing Markov chains [12,Chap. 3], EH is called "mean time to absorption" and it is known that the V k 1 (i)'s are given by the solution of (12). It is clear that V k 1 (0) = 1 because (Ĩ (t)) t∈N remains in state zero after having visited (absorbing) state zero. Substituting this in (33) proves (10).
We are now in a position to state and prove our main result.
Theorem 2 Let a JQN be given such that ρ i < 1, ∀i. Then, where c is given by (11).
Proof Recall that u −n+1→0 is a sequence of n independent events from time −n + 1 to time 0. For simplicity of notation, let x * (y) where n * def = n * (y) is the backward coupling time between 0 and y, i.e., min{n ∈ N | φ (n) (0, u −n+1→0 ) = φ (n) (y, u −n+1→0 )}. By definition, n * (y) has the same distribution as τ (y, 0), and therefore n * < ∞ almost surely by Theorem 1. We have The first inequality comes from the fact that φ(·) is monotone; the following equality comes from Theorem 1; the last inequality follows from (6).
Using Propositions 2 and 4, if ρ i < 1 for all i, then In (47), we have used that the mean stationary number of jobs in queue i of a Jackson network with infinite buffers is geometric. Now, for y = C − e i and following the recursion along dimension i we have Solving the recursion along each dimension, we get (34).
Theorem 2 provides sufficient conditions to characterize the qualitative behavior of the mean coupling time of JQNs when the network size, i.e., the number of queues and the buffer sizes, varies. In the following remark, we derive such qualitative behavior.

Remark 1
We recall that the V k 1 (i)'s, as defined in (12), are interpreted as the mean number of visits performed to each queue of our JQN by a job that starts visiting the network from queue k 1 ; see, e.g., [18,Chap. 7]. Suppose we are given a sequence of JQNs indexed by r where M (r ) and C (r ) are increasing in r ∈ N. Furthermore, in each JQN of the sequence, assume that i (solution of (12)) for all i, j ≥ 0, are uniformly bounded by a constant independent of the sequences (M (r ) ) r ≥0 and (C (r ) ) r ≥0 . As a consequence, c = O((M (r ) ) 2 ), and Theorem 2 shows that the mean coupling time of the r -th JQN is Formula (50) shows that the coupling time of JQNs is polynomial in both M and C. This improves the bounds developed in [2] (a comparison with those bounds follows below).
The following remark comments on the assumption of Theorem 2.
Remark 2 For a network with arbitrary topology and service rates, the assumption ρ i = Γ i /μ i < 1 corresponds to what can be called an "hyper-stable" queue. Indeed, in such a case, the arrival rate in queue i is always smaller than its service rate, even when all input queues in i are not empty and send jobs to i at full rate. This does not mean that the network is in a light-load regime. In fact, by varying the external arrival rates and the routing probabilities of some downstream queues, each queue can have any load between 0 and 1 within this assumption. From a technical standpoint, this assumption implies that the drift along each dimension i of the state space of the underlying Markov chain transitions can be bounded uniformly (from above) by a positive constant, i.e., Γ i . This structure is essentially exploited in the proof of Proposition 4. This bounded drift assumption suggests that our main result can be proved using some Lyapunov function on the expected time to empty the system. However, this approach cannot work as we discuss in Sect. 3.1.

Time to empty the system
It is easy to see that τ (0, C) can be stochastically bounded from above by the hitting time from C to 0, i.e., min{n : φ (n) (C, u 1→n ) = 0}. In turn, using a simple coupling argument, we have and, using Kac's lemma, E min{n > 0 : φ (n) (0, u 1→n ) = 0} = 1/π(0), where π(0) denotes the stationary probability of being in state 0. However, it is known that 1/π(0) grows exponentially in M: in fact, if the C i 's are sufficiently large and the overall mean arrival rate to i, say λ i , is such that Remark 3 The argument above shows that we cannot use an argument based on a Lyapunov function and on the expected time to empty the system to prove polynomial bounds on Eτ (0, C) when M and the C i 's increase, even when the drift to 0 of the underlying Markov chain is uniformly bounded as in the assumption of Theorem 2.

Comparison with [2]
The main difference of our bound (34) with respect to the ones in [2] is that (34) is of the order of a polynomial in both M and C (see (50)). In other words, it is not exponential in either M or the C i 's. For the case of JQN with cycles, Proposition 4.3 in [2] gives some conditions to prove that the coupling time of JQN is linear in the buffer size of a single queue. In particular, the constant of proportionality is not specified and, most importantly, the number of queues M is considered as a constant (in contrast with our approach). If all the buffer sizes are allowed to grow, these conditions are stronger than the ones in Theorem 2. In fact, there it is required hyperstability (as we do) plus the condition on the network topology that each queue can either receive jobs from outside or send jobs to outside, i.e., γ i0 > 0 or γ 0i > 0, for all i. a precision threshold. This notion is very useful for getting probabilistic guarantees for Monte Carlo simulation. More precisely, let P be the transition matrix of the chain and let d(n) def = max x∈S ||P n (x, ·) − π(·)|| T V , where || · || T V is the total variation distance. The mixing time with uncertainty , m( ), is the time it takes for the total variation distance between the distribution of the chain at time n and the stationary distribution π , to be less than or equal to , uniformly on the initial state. Namely, The mixing time is usually defined with uncertainty 1/4: T mi x def = m(1/4) (see [15]). If we consider the Markov chain defined by a monotone QN with finite buffer, then its mixing time is bounded from above by its coupling time. This is stated in the following proposition, whose proof is based on rather classical arguments that relate the mixing and the coupling time.

Proposition 5 Let a monotone QN with buffer vector C be given. Then, for all
Proof Theorem 5.2 in [15] says that ||P n (x, ·) − P n (y, ·)|| T V ≤ Pr(τ (x, y) > n). From this point, let s(n) be the total variation between to transient measures of the chain. We have where the last equality comes from the monotonicity of the chain. Now, the mixing can also be bounded since by definition d(n) ≤ s(n), where the penultimate inequality is the Markov inequality. By replacing by 1/4, one gets T mi x ≤ 4Eτ (0, C). Again, for an arbitrary integer k, d(km( )) ≤ s(km( )).
Thus, by means of Proposition 5, it is possible to exploit the results on the coupling time presented in Sect. 3 to bound the mixing time.

Queueing networks with blocking
The proof technique presented in Sect. 3 can be applied to other types of QNs with finite buffers, provided that they remain monotone and hyper-stable. For example, the same approach can be used for other types of QNs -with other types of mechanisms to deal with full buffers (e.g,. with blocking), -with state-dependent routing (e.g., join-the-shortest-queue), -with stations having service rates depending on the number of jobs in their queues (e.g., −/M/k/C queues), and -cases where jobs belong to multiple classes.
In all these cases, the coupling time can be shown to be polynomial both in the number of queues and buffer sizes. Bounds on the mean mixing time trivially follow by using Proposition 5, which only requires the monotonicity of the chain.
Here, we provide a detailed treatment of the blocking case. There exist many blocking mechanisms, and the most classical ones include blocking after service, blocking before service, repetitive service or re-circulate blocking. They are all similar in terms of coupling and mixing times. We consider the case with repetitive service, where upon completion of service at queue i, a job is sent to queue j with probability p i j and if j has no available slots, queue i repeats the service for that job. Afterward, a new target queue is selected according to the routing probabilities of i in a i.i.d. manner. As mentioned before, this type of blocking is called Repetitive-Service (RS) blocking, and even for this class of JQNs no product-form is known for the stationary distribution; see [4] for a background and further details.

Remark 4
The evolution of a JQN with RS blocking is thus identical to the model introduced in Sect. 2 except for routing events to queues that are full, for which the state remains unchanged (a state is still represented by the number of jobs in each queue).
Furthermore, quantity Γ i , defined in (8), still represents the minimal upper bound on the mean arrival rate of jobs in queue i with respect to any network state. We also observe that JQNs with RS blocking are monotone (as shown in [2]). Proof For any possible event a ∈ A , k = 1, . . . , M, and any state y such that y, y−e k ∈ S , (4) still holds true. In fact, for every y ∈ S and k such that y−e k ∈ S , cases i.a)-ii.c) of the proof of Proposition 2 holds for JQNs with RS blocking because their dynamics are equivalent to the ones of Sect. 2 (within these cases). For each routing events a = a i j , we have the following cases: iii.a) If y : 1 < y i ≤ C i and k = i, then φ(y, a) = φ(y − e k , a) + e k iii.b) If y : y i = 1 and k = i and y j < C j , then φ(y, a) = φ(y − e k , a) + e j iii.c') If y : y i = 1 and k = i and y j = C j , then φ(y, a) = φ(y − e k , a) + e k iii.d) If y : y j < C j and k = j, then φ(y, a) = φ(y − e k , a) + e k iii.e') If y : y j = C j and k = j and y i = 0, then φ(y, a) = φ(y − e k , a) + e j iii.f) If y : y j = C j and k = j and y i = 0, then φ(y, a) = φ(y − e k , a) + e k iii.g) If k = i and k = j, then φ(y, a) = φ(y − e k , a) + e k , i.e., only cases iii.c) and iii.e) change. Therefore, (4) holds true and, at time τ (0, y−e i ), both trajectories from 0 and y will be at distance one.
The proofs of Proposition 4 and Theorem 2 are based on a stochastic dominance property of JQNs with infinite buffers with respect to their finite counterparts (see Lemma 2). In JQNs with RS blocking, this property is lost and thus, those proofs cannot be applied directly to obtain the same bound. However, a quadratic bound in the C i 's can be derived as we show in the following.
The proof of the following proposition proceeds on the same lines of the proof of Proposition 4.

Proposition 7 Let a JQN with RS blocking be given such that
Then, for any x ∈ S and k 1 such that x − e k 1 ∈ S , where c is and (V k 1 (1), . . . , V k 1 (M)) is the unique solution of the linear system Proof Similarly to (14), let for n > 0 and I (0) 1}. By construction, it is clear that τ (x − e k 1 , x) = min{n : I (n) = 0}. Therefore, in the following, we study E min{n : I (n) = 0}.
For h ≥ 1, let T (h) be as in the proof of Proposition 7 be the random variable "time where the process (I (n)) n∈N makes the h-th jump if it exists, and ∞ otherwise"," and let T (0) = 0.
Using the same argument in (20), we have and therefore we derive a bound for E T (h ) − T (h − 1)|H = h and on Pr (H = h), where we recall that H is the random variable of the total number of jumps of the process (I (n)) n∈N .
As in Lemma 3, E(T (t) − T (t − 1)|H = h) can be upper bounded by the mean hitting time from some state to 0 of an M/M/1/C k t with arrival rate Γ k t and service rate μ k t , provided that k t = I (n) when T (t − 1) ≤ n < T (t). Using Proposition 3, in the worst case we trivially have and substituting in (59) we obtain Now, letĨ (t) be the value of (I (n)) n∈N at the time of its t-th jump. The jumps of (I (n)) n∈N to state k > 0 are only due to occurrences of cases iii).b and iii).c' (defined in the proof of Proposition 6), which means that at the time n + 1 of a jump of I we have either φ (n) (x, u 1→n ) I (n) = 1 or φ (n) (x, u 1→n ) I (n) = C I (n) . Therefore, ifĨ (t) = j and n t + 1 is the time of the t-th jump, then Pr(Ĩ (t + 1) = i|Ĩ (t) = j and φ (n t ) (x, u 1→n t ) I (n t ) = 1) = γ ji M k=0,k = j γ jk Pr(Ĩ (t + 1) = i|Ĩ (t) = j and φ (n t ) (x, u 1→n t ) I (n t ) = C I (n t ) ) = γ i j M k=0,k = j γ k j for all i, j. Let V k 1 (i) denote the mean number of visits that (Ĩ (t)) t∈N makes to state i when it is started in k 1 . The above conditional jump probabilities imply and we can upper bound EH easily as done in the proof of Proposition 4 by summing the V k 1 (i)'s. It is clear that V k 1 (0) = 1, because (I (n)) n∈N remains in state zero after having visited zero, and upper bounds on the V k 1 (i)'s are given by (57), which uses (64). Substituting this in (59) proves (55).
We are now in a position to prove the following result.

Theorem 3 Let a JQN with RS blocking be given such that
Following the recursion along dimension i, we have Solving the recursion in each dimension, we find Eτ (0, C) ≤ Mc max j C j i C i .
Based on the same conditions of Remark 1 and the ones of Theorem 3, we conclude that the qualitative behavior of the mean coupling time when both M and C increase is

Conclusions
We have proposed an approach to obtain bounds on the coupling and mixing times of a class of Markovian queueing networks, which are related to the efficiency of sampling from their stationary distribution in an exact or approximate manner, respectively. Our results provide conditions to show that the coupling time of Jackson queueing networks with finite buffers and arbitrary topology grows slowly (polynomially) when both the number of queues and the size of all buffers increase. Our bounds significantly improve the best bounds known in the literature and extend the ones found for acyclic Jackson QNs in [2,20]. We have also shown that the mixing time of finite Jackson QNs behaves in a similar way up to a factor depending on the accuracy of the desired samples. With minor variations in the proof-technique used, we extended our approach to a class of queueing networks with blocking obtaining similar (polynomial in M and the C i 's) bounds. Other possible extensions include monotone networks where (i) routing can be state-dependent, e.g., join-the-shortest-queue; (ii) stations have service rates depending on the number of jobs in their queues, e.g., −/M/k/C queues; and (iii) cases where jobs belong to multiple classes. Again, these follow by adapting (with minor variations) the proofs of Propositions 2 and 4 and Theorem 2 or 3. In conclusion, our bounds promote both CFTP and Monte Carlo methods as efficient tools for evaluating the stationary performance of monotone queueing networks. While the former produces exact samples, the latter is more efficient, and one can choose between these two techniques depending on which of these properties suits one's needs better.