Queueing and Stochastic Networks

Motivation

Queueing control problems (admission control, routing, service rate selection) are naturally modelled as Markov decision chains with a countably infinite state space and unbounded costs, for instance a holding cost that grows with the queue length. For such models the long-run average cost criterion is often the relevant one, and the central question is whether an optimal stationary policy exists and can be computed from an average cost optimality equation (ACOE). Chapter 7 of Linn I. Sennott, Stochastic Dynamic Programming and the Control of Queueing Systems (Wiley, 1999, doi:10.1002/9780470317037) develops a verifiable set of conditions, the (SEN) assumptions, under which an average cost optimality inequality (ACOI) holds and yields an optimal stationary policy. The inequality may be strict (Example 7.3.1), and an optimal policy may induce a Markov chain without positive recurrent states.

Sections 7.4 and 7.5 answer two practical questions: when is the ACOI in fact an equation, and how can (SEN) be checked in a concrete model? The answer culminates in the (BOR) assumptions, which require only one well-behaved stationary policy and the finiteness of a set of low-cost states.

According to the book's bibliographic notes (p. 163): the (BOR) assumptions modify a line of development due to Borkar (SIAM J. Control Optim. 22, 1984, and 27, 1989; monograph 1991) and are weaker than his original conditions; the proof that (BOR) implies (SEN) is from Cavazos-Cadena and Sennott (Oper. Res. Letters 11, 1992), and the version of (BOR) used here is from Sennott (Prob. Eng. Inform. Sci. 7, 1993). Proposition 7.5.5 and the (CAV*) assumptions go back to Cavazos-Cadena (Kybernetika 25, 1989); Proposition 7.5.3 and Corollary 7.5.4 to Sennott (Oper. Res. 37, 1989).

Setting

A Markov decision chain consists of a countable state space $S$, finite nonempty action sets $A_i$, nonnegative finite costs $C(i,a)$ and transition probabilities $P_{ij}(a)$. A policy $\theta$ may use the whole history and randomize. For $\alpha\in(0,1)$ the discount value function is $V_\alpha(i)=\inf_\theta V_{\theta,\alpha}(i)$, the infimum of $\sum_t\alpha^tE_\theta[C(X_t,A_t)\mid X_0=i]$; the average cost of $\theta$ is $J_\theta(i)=\limsup_n\frac1nE_\theta[\sum_{t<n}C(X_t,A_t)\mid X_0=i]$ and the minimum average cost is $J(i)=\inf_\theta J_\theta(i)$. All of these lie in $[0,\infty]$.

For a distinguished state $z$ the relative value is $h_\alpha(i)=V_\alpha(i)-V_\alpha(z)$. The (SEN) assumptions are: (SEN1) $(1-\alpha)V_\alpha(z)$ is bounded on $(0,1)$; (SEN2) $h_\alpha\le M$ for a finite function $M\ge0$; (SEN3) $h_\alpha\ge-L$ for a finite constant $L\ge0$. Under (SEN), $J=\lim_{\alpha\to1^-}(1-\alpha)V_\alpha(i)$ is a finite constant, and a limit function $h$ is a pointwise limit of $h_{\beta_n}$ along some $\beta_n\to1^-$. The ACOI and ACOE read

$$J+h(i)\ \ge\ (\text{resp. }=)\ \min_{a\in A_i}\Big\{C(i,a)+\sum_jP_{ij}(a)h(j)\Big\},\qquad i\in S.$$

For a nonempty set $G$ the first passage time is $T=\min\{n\ge1:X_n\in G\}$. The class $\Re(i,G)$ consists of the policies that, from $i$, enter $G$ with probability one in finite expected time $m_{iG}(\theta)$; $\Re^*(i,G)$ adds a finite expected first passage cost $c_{iG}(\theta)=E_\theta[\sum_{t<T}C(X_t,A_t)]$. A (randomized) stationary policy $d$ is $z$ standard if the Markov chain it induces has $m_{iz}<\infty$ and $c_{iz}<\infty$ for every $i$; it then has a single positive recurrent class $R_d\ni z$ and a finite constant average cost $J_d$.

Formalization targets

Goal: Theorem 7.5.6

Assume (BOR): (BOR1) a $z$ standard policy $d$ exists; (BOR2) for some $\varepsilon>0$ the set $D=\{i: C(i,a)\le J_d+\varepsilon\text{ for some }a\}$ is finite; (BOR3) every $i\in D-R_d$ can be reached from $z$ by some $\theta_i\in\Re^*(z,i)$. Then (SEN) holds and every limit function satisfies the ACOE; every average cost optimal stationary policy $e$ has a positive recurrent state in

$$D(e)=\{i: C(i,e)\le J+\varepsilon\},$$

at most $|D(e)|$ positive recurrent classes and no null recurrent class; and a policy realizing the minimum in the ACOE satisfies $e\in\Re^*(i,D(e)\cap R(e))$ for every $i$.

Milestones

• Lemma 7.4.1: $h_\alpha(i)\le c_{iz}(\theta_i)$ for $\theta_i\in\Re^*(i,z)$, hence (SEN2).
• Lemma 7.4.2: $h(i)\le c_{iG}(\theta)-Jm_{iG}(\theta)+E_\theta[h(X_T)]$ for $\theta\in\Re(i,G)$ under an integrability condition.
• Theorem 7.4.3: four sufficient conditions for equality in the ACOI at a state.
• Lemma 7.5.2: $J_d=(1-\alpha)\sum_{i\in R}\pi_i(d)V_{d,\alpha}(i)$ for a $z$ standard $d$.
• Proposition 7.5.3: a $z$ standard policy gives (SEN1–2).
• Corollary 7.5.4: on $S=\{0,1,\dots\}$, increasing $V_\alpha$ plus a $0$ standard policy gives (SEN), with nonnegative increasing limit functions.
• Proposition 7.5.5: an optimal stationary policy has a positive recurrent state of cost at most $J+\varepsilon$, reachable from $i$, when (7.33) holds.
• Corollaries 7.5.9 and 7.5.10: the (CAV) and (CAV*) conditions imply (BOR).

Significance

Theorem 7.5.6 reduces the verification of the ACOE for a queueing model to three checks that do not involve the discount value function: exhibit one stationary policy with finite mean return times and costs to a fixed state (typically a stable "serve at maximal rate" policy), check that low costs occur on a finite set (automatic when the holding cost grows without bound, Corollaries 7.5.9–7.5.10), and check reachability of finitely many states. Its conclusions go beyond existence: optimal stationary policies induce chains with positive recurrent classes located in a known finite set, and ACOE-realizing policies reach them in finite expected time and cost. This is what makes value iteration and approximating-sequence methods in later chapters of the book applicable to these models.

The results are proved in the book. The present mission produces machine-checked statements of the first passage calculus for general (history-dependent, randomized) policies, of (SEN) and limit functions, and of the chain of implications from (CAV*) to the ACOE. No machine-checked version of these statements is known.

Difficulty

The obvious approach to the ACOE is to pass to the limit $\alpha\to1^-$ in the discount optimality equation. Exchanging this limit with $\sum_jP_{ij}(a)h_\alpha(j)$ requires a dominating function, and (SEN2) only gives a pointwise bound $M$ whose expectation may be infinite; Fatou's lemma then yields only the inequality. Obtaining equality requires tracking first passages to sets and showing that the discrepancy $\Phi$ vanishes along them, which in turn needs finiteness of $c_{iG}$ that is not assumed but has to be derived. On the recurrence side, the average cost criterion is a limit superior of Cesàro averages over a countable state space, and mass can escape to infinity; the finiteness of the set $D$ is what prevents an optimal policy from spending its time in transient or null recurrent states, and turning that into positive recurrence requires the renewal-type identities of Appendix C.

Formalization scope

States form a countable type $S$; action sets are nonempty Finsets; costs are in ℝ≥0; transition probabilities are ℝ≥0∞-valued with row sums one on admissible actions. Policies are general: a history is a state sequence and an action sequence, and all probabilities and expectations (hitting probabilities, $m_{iG}$, $c_{iG}$, $P_\theta(X_T=j)$, $Q^{(n)}_{ij}$) are computed from the history probabilities of the process. $V_\alpha$, $J_\theta$, $m_{iG}$ and $c_{iG}$ take values in $[0,\infty]$; $m_{iG}=\infty$ when $G$ is missed with positive probability; the first passage time satisfies $T\ge1$. $h_\alpha$ and $\sum_jP_{ij}(a)h(j)$ are in the extended reals, with the book's convention that a function bounded below has an expectation in $(-\infty,+\infty]$. Limit functions are real valued. Positive recurrence, communicating classes and steady state probabilities $\pi_j=(m_{jj})^{-1}$ are the notions for the chain induced by a (randomized) stationary policy. $J_d$ is the average cost of $d$ from $z$.

A formalization in which the ACOE is asserted for some convenient function instead of every limit function, or in which $|D(e)|$ is a natural-number cardinality that vanishes on infinite sets, would trivialize part of the goal; the statements quantify over all limit functions and use Set.encard.

A complete development needs: history-dependent policies and their path laws on countable spaces; first passage decompositions (strong Markov property at $T$); Abelian limits of $\sum_t\alpha^tP(T=t)$; Fatou and dominated convergence for series; and the renewal reward theorem for positive recurrent classes (Appendix C of the book). The first passage and Markov chain layer is reusable beyond this mission. Proofs of individual milestones, and sharper statements of the Appendix C facts they use, are welcome.

Selected references

• L. I. Sennott, Stochastic Dynamic Programming and the Control of Queueing Systems, Wiley Series in Probability and Statistics, John Wiley & Sons, 1999. doi:10.1002/9780470317037
• V. S. Borkar, "On minimum cost per unit time control of Markov chains", SIAM J. Control Optim. 22 (1984), 965–978.
• V. S. Borkar, "Control of Markov chains with long-run average cost criterion: the dynamic programming equations", SIAM J. Control Optim. 27 (1989), 642–657.
• V. S. Borkar, Topics in Controlled Markov Chains, Pitman Research Notes in Mathematics 240, Longman, 1991.
• R. Cavazos-Cadena, "Weak conditions for the existence of optimal stationary policies in average Markov decision chains with unbounded costs", Kybernetika 25 (1989), 145–156.
• R. Cavazos-Cadena and L. I. Sennott, "Comparing recent assumptions for the existence of average optimal stationary policies", Oper. Res. Letters 11 (1992), 33–37.
• L. I. Sennott, "The average cost optimality equation and critical number policies", Prob. Eng. Inform. Sci. 7 (1993).
• L. I. Sennott, "Average cost optimal stationary policies in infinite state Markov decision processes with unbounded costs", Operations Research 37 (1989), 626–633. doi:10.1287/opre.37.4.626
• K. L. Chung, Markov Chains with Stationary Transition Probabilities, 2nd ed., Springer, 1967.

Motivation

Almost every model in queueing theory is analysed through a Markov chain. The number of customers in an M/M/c queue is a continuous-time birth–death chain; the number left behind by departing customers of an M/G/1 queue is a discrete-parameter chain on $\{0,1,2,\dots\}$ (the imbedded Markov chain); networks of queues are chains on vectors of queue lengths. Before any steady-state formula (Erlang's formulas, the Pollaczek–Khintchine formula, product forms) can be used, one has to know that the chain has a steady state at all: that it is positive recurrent, so that a stationary distribution exists and equals the limiting distribution.

Chapter 1 of Gross, Shortle, Thompson and Harris, Fundamentals of Queueing Theory (4th ed., Wiley 2008, DOI 10.1002/9781118625651), collects the two ingredients the rest of the book stands on: the Poisson process with its exponential interarrival times (§§1.7–1.8), and the classification theory of discrete-parameter Markov chains (§1.9), ending with Foster's criterion (Theorem 1.2), a sufficient condition for positive recurrence in terms of a drift inequality. The criterion goes back to F. G. Foster, On the stochastic matrices associated with certain queuing processes, Ann. Math. Statist. 24 (1953) (DOI 10.1214/aoms/1177728976), and is the ancestor of the Foster–Lyapunov method used for stability of queueing networks and stochastic systems.

This mission is the first of a series formalizing the book chapter by chapter.

Setting

A homogeneous discrete-parameter Markov chain on $\{0,1,2,\dots\}$ is given by a transition matrix $P=\{p_{ij}\}$ with $p_{ij}\ge0$ and $\sum_j p_{ij}=1$ for every $i$. The $m$-step transition probabilities $p_{ij}^{(m)}$ are the entries of $P^m$.

The first-passage probability $f_{ij}^{(n)}$ is the probability that the chain started in $i$ enters $j$ for the first time at step $n\ge1$; for $i=j$ it is the probability of first return at step $n$. The return probability is $f_{jj}=\sum_{n\ge1}f_{jj}^{(n)}$ and the mean recurrence time is $m_{jj}=\sum_{n\ge1}n f_{jj}^{(n)}\in[0,\infty]$. A state is positive recurrent if $f_{jj}=1$ and $m_{jj}<\infty$; the chain is positive recurrent if every state is.

The chain is irreducible if for every pair of states $(i,j)$ some $p_{ij}^{(n)}$ is positive, and aperiodic if for every state $k$ the greatest common divisor of $\{n\ge1:p_{kk}^{(n)}>0\}$ is $1$. A stationary distribution is a probability vector $\pi$ with $\pi=\pi P$, i.e. $\pi_j=\sum_i\pi_i p_{ij}$ for every $j$.

For the Poisson part, $T_0,T_1,\dots$ are independent interarrival times, each exponentially distributed with rate $\lambda>0$; the arrival epochs are $S_n=T_0+\dots+T_{n-1}$, and $N(t)=\#\{n\ge1:S_n\le t\}$ counts the arrivals in $[0,t]$.

Formalization targets

Goal: Theorem 1.2 (Foster's criterion)

An irreducible, aperiodic chain is positive recurrent if there exist $x_j\ge0$ with

$$\sum_{j=0}^\infty p_{ij}x_j\le x_i-1\quad(i\ne0),\qquad\sum_{j=0}^\infty p_{0j}x_j<\infty .$$

Milestones: the Markov chain theorems

• Theorem 1.1(a). In an irreducible, positive recurrent chain, $\pi_j=1/m_{jj}$ is a stationary distribution, and it is the only one.
• Theorem 1.1(c). If moreover the chain is aperiodic and all moments of $\pi$ are finite, then $\lim_{m\to\infty}p_{ij}^{(m)}=\pi_j$ for all $i,j$.

Milestones: the Poisson process and the exponential distribution

• Eqs. (1.11)–(1.14). The unique solution of $p_0'=-\lambda p_0$, $p_n'=-\lambda p_n+\lambda p_{n-1}$ with $p_0(0)=1$, $p_n(0)=0$ is $p_n(t)=(\lambda t)^n e^{-\lambda t}/n!$.
• Eq. (1.15). With exponential interarrival times,

$$\Pr\{N(t)\le n\}=\int_t^\infty\frac{\lambda(\lambda x)^n}{n!}e^{-\lambda x}dx=\sum_{i=0}^n\frac{(\lambda t)^ie^{-\lambda t}}{i!}.$$

• Eq. (1.16). Given $N(L)=k$, the arrival epochs have density $k!/L^k$ on $\{0<t_1<\dots<t_k<L\}$.
• Eq. (1.17) and its converse (p.21). The exponential law satisfies $\Pr\{T\le t_1\mid T\ge t_0\}=\Pr\{0\le T\le t_1-t_0\}$, and it is the only continuous distribution on $[0,\infty)$ that does.
• Nonhomogeneous Poisson law (p.22). With a continuous rate $\lambda(t)$ the forward equations have the unique solution $p_n(t)=e^{-m(t)}m(t)^n/n!$, $m(t)=\int_0^t\lambda(s)\,ds$.

Significance

Foster's criterion reduces positive recurrence, a statement about return times, to exhibiting one test function $x$ with negative drift outside a single state. In the book it is the tool that establishes the existence of steady state for imbedded chains of the M/G/1 and G/M/1 queues (Chapter 5); its generalizations are the standard stability proofs for queueing networks. Theorem 1.1 then supplies what positive recurrence buys: the stationary distribution exists, is unique, equals $1/m_{jj}$, and is the limit of the transition probabilities. The Poisson results justify the "Markovian" arrivals and services of Chapters 2–4.

All of these results are classical and proved in the literature; the book states Theorems 1.1 and 1.2 without proof. The Prove2Me platform already holds machine-checked versions of related Markov chain theorems in other missions (Levin–Peres–Wilmer's and Durrett's countable-chain convergence theorems), stated with different definitions and hypotheses. What this mission adds is a formal development in the book's own terms — first-passage probabilities $f_{jj}^{(n)}$, mean recurrence times $m_{jj}$, gcd periodicity — on which the later missions of the series (imbedded chains, birth–death processes) can build, together with a formal proof of Foster's criterion, which is not on the platform.

Difficulty

For Foster's criterion the natural first step, taking expectations of the drift inequality along the chain, only shows that the expected value of $x$ decreases while the chain stays away from $0$. Turning that into a bound on the expected return time to $0$ requires an optional-stopping or telescoping argument over a random time, with the value $x$ possibly unbounded, and a separate argument that positive recurrence of state $0$ propagates to all states of an irreducible chain. The book's hypotheses include aperiodicity, which the argument does not use.

For Theorem 1.1, identifying the stationary distribution with $1/m_{jj}$ requires relating the matrix powers $P^n$ to the first-passage probabilities (a renewal decomposition), and uniqueness over countably many states needs care with infinite sums. For the Poisson results, the difficulty is measure-theoretic: the distribution of the sum of $n+1$ exponential variables, and conditioning on the event $\{N(L)=k\}$ for the order-statistics property.

Formalization scope

States are natural numbers; the transition matrix is a real function $p:\mathbb N\times\mathbb N\to\mathbb R$ with nonnegative entries and rows summing to one (as a convergent series). The return probability and the mean recurrence time are valued in $[0,\infty]$, so null recurrence ($m_{jj}=\infty$) is representable. Irreducibility is the per-pair notion. Stationary equations are stated componentwise with convergent series.

In Foster's criterion the series $\sum_j p_{ij}x_j$ are required to converge for every $i$, which is the book's condition $\sum_j p_{0j}x_j<\infty$ together with the finiteness implicit in the inequalities for $i\ne0$; $x$ is real-valued and nonnegative. Dropping the convergence requirement would let a divergent row series (whose Lean sum is $0$) satisfy the inequality vacuously; allowing $x_j=\infty$ would make the hypothesis trivially satisfiable. Neither is permitted.

The closed forms stated explicitly are: $\pi_j=1/m_{jj}$ (Theorem 1.1(a), with both existence and uniqueness), the Poisson probabilities $(\lambda t)^ne^{-\lambda t}/n!$ (1.14), the Erlang tail integral and the Poisson CDF (1.15), the density $k!/L^k$ (1.16), and $e^{-m(t)}m(t)^n/n!$ for the nonhomogeneous law. Equations (1.14) and the nonhomogeneous law are stated as "solves the equations with the initial conditions if and only if equals the closed form", so both existence and uniqueness are asserted.

The Poisson results use random variables on a probability space, with Mathlib's expMeasure for the exponential law and cond for conditional probability. The derivation of the forward equations from the $o(\Delta t)$ axioms of §1.7 is not formalized; the Poisson law is reached from the equations and, separately, from exponential interarrival times.

Not formalized: Theorem 1.1(b) and the word "ergodic" in 1.1(c), which rest on the book's informal notion of ergodicity; Theorem 1.3, whose phrase "for Theorem 1.1 to be valid" for a continuous-time chain is not pinned down.

The Markov chain definitions are reusable by every later mission that studies an imbedded chain. Contributions welcome: proofs of the milestones, and supporting lemmas (Chapman–Kolmogorov, renewal decomposition of $p_{jj}^{(n)}$, class properties of recurrence).

Selected references

• D. Gross, J. F. Shortle, J. M. Thompson, C. M. Harris, Fundamentals of Queueing Theory, 4th ed., Wiley, 2008. https://doi.org/10.1002/9781118625651
• F. G. Foster, On the stochastic matrices associated with certain queuing processes, Annals of Mathematical Statistics 24 (1953), 355–360. https://doi.org/10.1214/aoms/1177728976

Why birth–death queues and Erlang's formulas

Every call center, hospital ward, cloud server pool and telephone exchange that is sized by formula is sized by one of a handful of explicit expressions from Markovian queueing theory. The two oldest are A. K. Erlang's: the Erlang-B (loss) formula of 1917, which gives the fraction of calls lost when $c$ trunks carry an offered load of $r$ erlangs, and the Erlang-C formula, which gives the probability that a customer of a $c$-server queue must wait. Both are still the default dimensioning rules of telecommunications and call-center workforce management (Gans, Koole & Mandelbaum 2003).

This mission formalizes Chapter 2, §§2.1–2.10, of Gross, Shortle, Thompson and Harris, Fundamentals of Queueing Theory, 4th ed. (Wiley 2008), which derives these formulas from a single result about birth–death processes and closes with the modern answer to the staffing question.

Timeline. Erlang (1917) obtained the loss formula; Vaulot (1927), Pollaczek (1932), Palm (1938) and Kosten (1948) completed its proof for general service times. Halfin and Whitt (1981) showed that in the $M/M/n$ queue the delay probability converges to a limit strictly between $0$ and $1$ exactly when the number of servers exceeds the offered load by an amount of order $\sqrt n$. This is the quality-and-efficiency-driven (QED) regime on which square-root staffing rests.

Setting

A birth–death process is a continuous-time Markov chain on the states $n \in \{0, 1, 2, \dots\}$ that moves from $n$ to $n+1$ at rate $\lambda_n \ge 0$ (a birth, or arrival) and, for $n \ge 1$, from $n$ to $n-1$ at rate $\mu_n > 0$ (a death, or departure). A steady-state solution is a probability sequence $\{p_n\}$ ($p_n \ge 0$, $\sum_n p_n = 1$) solving the global balance equations (2.1):

$$(\lambda_n + \mu_n)p_n = \lambda_{n-1}p_{n-1} + \mu_{n+1}p_{n+1}\ (n \ge 1), \qquad \lambda_0 p_0 = \mu_1 p_1.$$

The queues of the chapter are birth–death processes with particular rates. The $M/M/1$ queue has $\lambda_n = \lambda$, $\mu_n = \mu$ and traffic intensity $\rho = \lambda/\mu$. The $M/M/c$ queue has $\lambda_n = \lambda$, $\mu_n = \min(n, c)\mu$ (2.30), offered load $r = \lambda/\mu$ and $\rho = r/c$. The $M/M/c/c$ loss system is the same with $\lambda_n = 0$ for $n \ge c$. The $M/M/\infty$ queue has $\mu_n = n\mu$.

The explicit functions are the Erlang-B formula

$$B(c, r) = \frac{r^c/c!}{\sum_{i=0}^{c} r^i/i!},$$

the Erlang-C formula, defined for $\rho = r/c < 1$,

$$C(c, r) = \frac{r^c/(c!(1-\rho))}{r^c/(c!(1-\rho)) + \sum_{n=0}^{c-1} r^n/n!},$$

and, with $\phi$, $\Phi$ the standard normal density and distribution function,

$$\alpha(\beta) = \frac{\phi(\beta)}{\phi(\beta) + \beta\Phi(\beta)}.$$

Formalization targets

Goal: the Halfin–Whitt theorem (§2.4, p.75)

For offered loads $0 < r_n < n$,

$$\lim_{n\to\infty} C(n, r_n) = \alpha \in (0,1) \iff \lim_{n\to\infty} \frac{n - r_n}{\sqrt n} = \beta > 0, \qquad \alpha = \alpha(\beta).$$

It is stated as three facts: $\alpha$ maps $(0, \infty)$ into $(0, 1)$; each $\alpha \in (0, 1)$ has exactly one preimage $\beta > 0$; and for every $\beta > 0$ the two limits are equivalent.

Milestones

1. (2.3)–(2.4): the steady-state solution of a general birth–death process, $p_n = p_0\prod_{i=1}^n \lambda_{i-1}/\mu_i$, and its existence if and only if $1 + \sum_{n\ge1}\prod_{i=1}^n \lambda_{i-1}/\mu_i < \infty$.
2. (2.9): $M/M/1$, $p_n = (1-\rho)\rho^n$, existing iff $\rho < 1$.
3. (2.31)–(2.32): the $M/M/c$ law, existing iff $\lambda/(c\mu) < 1$.
4. (2.33): $L_q = r^c\rho\,p_0/(c!(1-\rho)^2)$.
5. (2.37)–(2.38): $1 - \sum_{n<c} p_n = C(c, r)$.
6. (2.52)–(2.53): the $M/M/c/c$ law and $p_c = B(c, r)$.
7. (2.54): $B(c, r) = rB(c-1, r)/(c + rB(c-1, r))$, $B(0, r) = 1$.
8. (2.55): $C(c, r) = cB(c, r)/(c - r + rB(c, r))$.
9. (2.57): $M/M/\infty$, $p_n = r^n e^{-r}/n!$.

Significance

The results. Items 1–9 are the working formulas of Markovian capacity planning: a stationary law for each basic model and the measures read off from it. (2.54) and (2.55) are how $B$ and $C$ are computed in practice, since the factorials of the closed forms overflow for $c > 170$. The Halfin–Whitt theorem is the reason the rule $c \approx r + \beta\sqrt r$ holds a fixed service level, and it is the entry point to the QED heavy-traffic literature (diffusion limits of many-server queues, Garnett–Mandelbaum–Reiman, Gamarnik–Momčilović).

Formalizing them. All results are classical and proved in the literature. The book states the Halfin–Whitt theorem without proof, and (2.54)–(2.55) are left to exercises. The formalization would supply machine-checked versions of the Erlang identities and of the Halfin–Whitt limit theorem. No Lean development of either was found on the platform when this mission was drafted. A related Erlang-B statement from Kelly and Yudovina is on the platform, stated with detailed balance on a finite state space.

Difficulty

The stationary laws are induction plus geometric and exponential series, and the Erlang identities are finite algebra. The difficulty is concentrated in the goal. $C(n, r_n)$ is a ratio of a Poisson-type tail to a truncated exponential sum in which both $n$ and $r_n$ grow. The naive route, substituting Stirling's formula term by term, fails: the sums have $\Theta(\sqrt n)$ significant terms, each of relative size $\exp(-k^2/2n)$, and the error has to be controlled uniformly over them. The converse direction also requires showing that $\alpha(\cdot)$ is strictly monotone. Without that, convergence of $C(n, r_n)$ does not force convergence of $(n - r_n)/\sqrt n$.

Formalization scope

Rates are real sequences indexed by $\mathbb N$, and a steady-state solution is a real sequence with HasSum p 1, nonnegative entries, and the balance equations (2.1) exactly as printed (global balance, not detailed balance). Every "the steady-state solution is X" is stated in both halves: X is a steady-state solution, and every steady-state solution equals X; the book's existence conditions ($\rho < 1$, $\lambda/(c\mu) < 1$, convergence of the series) are part of the statements. The $M/M/c/c$ system is the $\mathbb N$-indexed process with $\lambda_n = 0$ for $n \ge c$, as §2.5 sets it up; the statement records that states above $c$ carry no mass.

The closed forms that are fixed in Lean: $\prod_{i=1}^n \lambda_{i-1}/\mu_i$ over Finset.Icc 1 n; $B(c, r)$ and $C(c, r)$ exactly as displayed above; $\phi$ = gaussianPDFReal 0 1, $\Phi$ = the CDF of gaussianReal 0 1; $W_q(0) = \sum_{n=0}^{c-1} p_n$, as evaluated on p.69; $L_q = \sum_{n > c}(n - c)p_n$ as a convergent series.

$C(c, r)$ is a total function in Lean, but its value for $r \ge c$ carries no meaning. The goal assumes $0 < r_n < n$ for $n \ge 1$, the book's standing condition $\rho < 1$. A statement about some other function with the same limiting behaviour, or with $B$ and $C$ left abstract, would not be this mission. Neither would one-directional or existence-only versions of the stationary laws.

Not included: the waiting-time distributions (2.28) and (2.39), which need an FCFS waiting-time model with arrival-point probabilities; the $M/M/c/K$ measures (2.45)–(2.48); finite-source and state-dependent models (§§2.8–2.10). Useful contributions beyond the milestones are Poisson tail estimates at the $\sqrt n$ scale and monotonicity of $\alpha(\beta)$. Both are reusable in other many-server heavy-traffic statements.

Selected references

• D. Gross, J. F. Shortle, J. M. Thompson, C. M. Harris, Fundamentals of Queueing Theory, 4th ed., Wiley, 2008. https://doi.org/10.1002/9781118625651
• S. Halfin, W. Whitt, Heavy-traffic limits for queues with many exponential servers, Operations Research 29(3), 567–588, 1981. https://doi.org/10.1287/opre.29.3.567
• N. Gans, G. Koole, A. Mandelbaum, Telephone call centers: tutorial, review, and research prospects, Manufacturing & Service Operations Management 5(2), 79–141, 2003. https://doi.org/10.1287/msom.5.2.79.16071
• A. K. Erlang, Solution of some problems in the theory of probabilities of significance in automatic telephone exchanges, Elektroteknikeren 13, 1917 (English translation in The Life and Works of A. K. Erlang, 1948).
• F. P. Kelly, E. Yudovina, Stochastic Networks, Cambridge University Press, 2014. https://doi.org/10.1017/CBO9781139565363

Motivation

Steady-state formulas describe a queue that has been running forever. Many practical questions are about a queue that has not: a call centre just after opening, a server just after a reset, a system under a burst of load. For these, the relevant quantity is the transient distribution $p_n(t) = \Pr\{N(t) = n\}$ of the number $N(t)$ in the system at a finite time $t$. It is also what determines how fast the steady state is approached, and it is needed for the busy period: the length of time a server stays busy once a customer arrives at an idle server.

For the single-server Markovian queue M/M/1 the transient distribution has an explicit closed form in modified Bessel functions. Its history is short and well documented. Ledermann and Reuter (1954) obtained it by spectral analysis of the birth–death process. Bailey (1954) found it by generating functions and Laplace transforms, and Champernowne (1956) by combinatorial methods. Bailey's route is the standard textbook derivation, and it is the one Gross, Shortle, Thompson and Harris outline in §2.11 of Fundamentals of Queueing Theory (4th ed., 2008). Abate and Whitt (1989) showed that computing with the resulting series is numerically delicate, which is one reason for having the formula pinned down exactly.

This mission formalizes §§2.11–2.12 of that book: the transient laws of M/M/1/1, M/M/1 and M/M/∞, and the M/M/1 busy period.

Setting

Customers arrive in a Poisson stream of rate $\lambda > 0$. Each service takes an exponential time of rate $\mu > 0$, and $\rho = \lambda/\mu$. The number in the system is a continuous-time Markov chain on $\{0, 1, 2, \dots\}$, and its state probabilities $p_n(t)$ satisfy the forward (differential–difference) equations. For M/M/1 started with $N(0) = i$ they are, for $t \ge 0$,

$$p_n'(t) = -(\lambda+\mu)p_n(t) + \lambda p_{n-1}(t) + \mu p_{n+1}(t)\ (n > 0), \qquad p_0'(t) = -\lambda p_0(t) + \mu p_1(t), \tag{2.72}$$

with $p_n(0) = 1$ if $n = i$ and $0$ otherwise. The other systems are variants:

• M/M/1/1, no waiting room: two states and equations (2.70).
• M/M/∞, ample service: the death rate in state $n$ is $n\mu$, giving (2.76).
• The busy-period system: (2.72) with $0$ made absorbing ($\lambda_0 = 0$) and $N(0) = 1$. Its $p_0(t)$ is the distribution function of the busy period $T_{bp}$.

A family $(p_n)$ solves a system on $[0,\infty)$ when each $p_n$ has, at every $t \ge 0$, the prescribed derivative (a right derivative at $t = 0$). It is a probability solution when $p_n(t) \ge 0$ and $\sum_n p_n(t) = 1$ for every $t \ge 0$. The modified Bessel function of the first kind is

$$I_n(y) = \sum_{k=0}^{\infty} \frac{(y/2)^{n+2k}}{k!\,(n+k)!}, \qquad I_{-n} = I_n,$$

and the Laplace transform of $f$ is $\bar f(s) = \int_0^\infty e^{-st} f(t)\,dt$ for $\operatorname{Re} s > 0$.

Formalization targets

Goal: the transient M/M/1 law, (2.75)

With $y = 2t\sqrt{\lambda\mu}$,

$$p_n(t) = e^{-(\lambda+\mu)t}\Big[\rho^{(n-i)/2} I_{n-i}(y) + \rho^{(n-i-1)/2} I_{n+i+1}(y) + (1-\rho)\rho^n \sum_{j=n+i+2}^{\infty} \rho^{-j/2} I_j(y)\Big].$$

The goal asserts five things for every $\lambda, \mu > 0$ and every $i$, with no restriction on $\rho$:

1. the series converges;
2. these functions solve (2.72);
3. they meet the initial condition;
4. they form a probability distribution at every $t$;
5. they are the only probability solution.

Milestones

1. (2.71): the M/M/1/1 solution $p_1(t) = \frac{\lambda}{\lambda+\mu}(1-e^{-(\lambda+\mu)t}) + p_1(0)e^{-(\lambda+\mu)t}$, and the matching formula for $p_0$.
2. (2.74) and Rouché's theorem: for $\operatorname{Re} s > 0$, the quadratic $(\lambda+\mu+s)z - \mu - \lambda z^2$ has exactly one zero in $|z| < 1$, namely $z_1 = (\lambda+\mu+s-\sqrt{(\lambda+\mu+s)^2-4\lambda\mu})/(2\lambda)$.
3. The transform of $p_0$: $\bar p_0(s) = z_1^{i+1}/(\mu(1-z_1))$.
4. The limit of (2.75): $p_n(t) \to (1-\rho)\rho^n$ if $\rho < 1$, and $p_n(t) \to 0$ if $\rho \ge 1$.
5. (2.77), M/M/∞: started empty, $p_n(t) = a(t)^n e^{-a(t)}/n!$ with $a(t) = (1-e^{-\mu t})\lambda/\mu$. The statement says that this family solves (2.76), is the unique probability solution, and has generating function $\exp((z-1)a(t))$.
6. The busy-period transform: $\bar p_0(s) = 2\mu/(s[\lambda+\mu+s+\sqrt{(\lambda+\mu+s)^2-4\lambda\mu}])$.
7. The busy-period density: $p_0'(t) = \sqrt{\mu/\lambda}\,e^{-(\lambda+\mu)t} I_1(2\sqrt{\lambda\mu}\,t)/t$.
8. (2.79): for $\lambda < \mu$, $E[T_{bp}] = 1/(\mu-\lambda)$ and $E[T_{bc}] = 1/\lambda + 1/(\mu-\lambda)$.

Significance

The formula (2.75) is the exact finite-time law of the most basic queue. It gives the rate at which M/M/1 approaches equilibrium, and it gives the distribution of the queue under overload ($\rho \ge 1$), where no steady state exists. It is the reference against which numerical transient methods, such as the uniformization of Chapter 8 of the same book, are checked. The busy-period density and its mean (2.79) enter server-utilisation and vacation models, and the Laplace-transform method used here recurs in the M/G/1 analysis of Chapter 5.

All of these results are classical and proved in the literature. None of them is machine-checked, as far as the platform's catalogue and Mathlib show. The chain from a countable system of linear ODEs, through generating functions and a root-location argument, to a Bessel series is a standard pattern in applied probability, and a formal version of it is what this mission adds. The formal statements also make explicit what the book leaves implicit: the sense in which the equations hold at $t = 0$, and the class in which the solution is unique.

Difficulty

The forward equations (2.72) form an infinite linear system. The obvious approach is to treat it like a finite system of ODEs, whose solution is a matrix exponential, and read off (2.75). That fails for two reasons. The generator is an infinite matrix, so its exponential needs a functional-analytic setting. And uniqueness is not automatic for infinite systems: it needs a class, such as probability solutions, and an argument that works in that class.

The Bessel form is a second, independent difficulty. The transform $\bar p_0(s)$ is fixed by a root-location argument in the complex plane. Inverting the transform, or verifying (2.75) directly, requires manipulating the three-term Bessel recurrence and exchanging infinite sums. The tail sum $\sum_{j} \rho^{-j/2} I_j$ has to be controlled uniformly enough to be differentiated term by term. For $\rho \ge 1$ the factor $(1-\rho)$ is non-positive, so the nonnegativity of $p_n(t)$ is not visible from the formula.

Formalization scope

Conventions committed to:

• Parameters. Rates are real with $\lambda, \mu > 0$, and $\rho = \lambda/\mu$. States are ℕ (Fin 2 for M/M/1/1).
• Solutions. "Solves on $[0,\infty)$" is HasDerivWithinAt on Set.Ici 0 at every $t \ge 0$. Uniqueness is asserted among solutions that are probability distributions at every time.
• Special functions. Half-integer powers of $\rho$ are real powers, and $I_{-m} = I_m$ is part of the definition. Laplace transforms are complex Bochner integrals over $(0,\infty)$, and each statement also asserts the integrability it needs. Square roots with positive real part are hypotheses $r^2 = (\lambda+\mu+s)^2 - 4\lambda\mu$, $\operatorname{Re} r > 0$.

The closed forms stated exactly as in the book are:

• (2.71);
• $z_1$ and $z_2$ of (2.74);
• $\bar p_0(s) = z_1^{i+1}/(\mu(1-z_1))$;
• (2.75), with the Bessel series of p.101;
• the M/M/∞ law and (2.77);
• the busy-period transform and density of p.102;
• (2.79).

The book derives (2.79) by a steady-state ratio argument valid for M/G/1. Here it is stated for M/M/1, as the mean of the explicit density.

A statement of (2.75) that only asserts the right-hand side is well defined, or checks only $n = 0$, is ruled out: the goal requires the ODE system, the initial condition, the probability property and uniqueness. For the same reason, the M/M/∞ law is tied to the system (2.76) and does not reduce to a Taylor expansion.

Needed infrastructure that Mathlib lacks:

• modified Bessel functions of integer order;
• Laplace transforms;
• a Rouché-type zero count or a direct root-location lemma;
• uniqueness for countable linear ODE systems with bounded or linearly growing rates.

The Bessel and Laplace definitions, and the uniqueness lemma for birth–death forward equations, are reusable beyond this mission. Contributions of those as separate lemmas are welcome.

Selected references

• D. Gross, J. F. Shortle, J. M. Thompson, C. M. Harris, Fundamentals of Queueing Theory, 4th ed., Wiley, 2008, §§2.11–2.12, pp.97–103. https://doi.org/10.1002/9781118625651
• N. T. J. Bailey, "A continuous time treatment of a simple queue using generating functions", J. Royal Statistical Society B 16 (1954) 288–291. https://doi.org/10.1111/j.2517-6161.1954.tb00172.x
• W. Ledermann, G. E. H. Reuter, "Spectral theory for the differential equations of simple birth and death processes", Phil. Trans. Royal Society A 246 (1954) 321–369. https://doi.org/10.1098/rsta.1954.0001
• D. G. Champernowne, "An elementary method of solution of the queueing problem with a single server and constant parameters", J. Royal Statistical Society B 18 (1956) 125–128. https://doi.org/10.1111/j.2517-6161.1956.tb00217.x
• J. Abate, W. Whitt, "Calculating time-dependent performance measures for the M/M/1 queue", IEEE Trans. Communications 37 (1989) 1102–1104. https://doi.org/10.1109/26.41165

Motivation

In many service systems a customer who finds every server busy does not join a queue. A caller who hears a busy signal hangs up and redials later; a request rejected by a saturated server is resent after a timeout; an aircraft that cannot land circles and tries again. These retrial queues are the subject of a substantial literature in telephone traffic engineering, computer networks and call-centre design, surveyed in the monograph of Falin and Templeton (1997) and the bibliography of Artalejo (1999). Their analysis is harder than that of ordinary queues: the blocked customers form an orbit whose size is part of the state, so even the simplest model is a two-dimensional Markov chain, and explicit stationary distributions are rare.

This mission is the fourth of a series formalizing Gross, Shortle, Thompson and Harris, Fundamentals of Queueing Theory (4th ed., Wiley 2008). Its goal is the explicit stationary distribution of the single-server retrial queue, Eq. (3.57) of §3.5.1, one of the few retrial models solvable in closed form. Chapter 3 of the book treats Markovian queues that are not birth–death processes: bulk arrivals, bulk service, Erlang phases, priority disciplines and retrials. The milestones also collect three capstone formulas from the chapter's other sections: the bulk-input queue, the partial-batch bulk-service queue, and Cobham's formula for nonpreemptive priorities (Cobham, 1954).

Setting

In the $M/M/1$ retrial queue customers arrive according to a Poisson process with rate $\lambda$ and are served one at a time by a single server, with exponential service times of mean $1/\mu$. An arrival that finds the server busy enters the orbit and stays there for an exponential time with mean $1/\gamma$, after which it tries again; each customer in orbit retries independently. No customer leaves because of impatience. With $N_s(t) \in \{0,1\}$ the number in service and $N_o(t)$ the number in orbit, the pair is a continuous-time Markov chain on states $\{i, n\}$, $i \in \{0,1\}$, $n \in \{0,1,2,\dots\}$. Writing $p_{i,n}$ for the steady-state probability of $\{i,n\}$, the rate-balance equations are

$$\begin{aligned} (\lambda + n\gamma)p_{0,n} &= \mu p_{1,n}, && n \ge 0, && (3.47)\\ (\lambda+\mu)p_{1,n} &= \lambda p_{0,n} + (n+1)\gamma p_{0,n+1} + \lambda p_{1,n-1}, && n \ge 1, && (3.48)\\ (\lambda+\mu)p_{1,0} &= \lambda p_{0,0} + \gamma p_{0,1}. && && (3.49) \end{aligned}$$

Following the book's convention (§1.9, and the footnote on p.118), a steady-state solution is a nonnegative solution of these equations whose total mass $\sum_n (p_{0,n} + p_{1,n})$ equals $1$. The traffic intensity is $\rho = \lambda/\mu$, and the partial generating functions are $P_0(z) = \sum_n z^n p_{0,n}$ and $P_1(z) = \sum_n z^n p_{1,n}$.

The other models of the mission use the same convention. In the bulk-input queue $M^{[X]}/M/1$, batches arrive at rate $\lambda$ with batch-size probabilities $c_n = \Pr\{X = n\}$, $n \ge 1$, and batch-size generating function $C(z) = \sum_n c_n z^n$. In the partial-batch bulk-service queue $M/M^{[K]}/1$, single arrivals come at rate $\lambda$ and the server serves up to $K$ customers together in an exponential time of mean $1/\mu$. In the nonpreemptive priority queue there are $r$ classes with rates $\lambda_k$ and $\mu_k$, loads $\rho_k = \lambda_k/\mu_k$ and cumulative loads $\sigma_k = \rho_1 + \cdots + \rho_k$.

Formalization targets

Goal: the stationary distribution (3.57)

For $\lambda, \mu, \gamma > 0$ and $\rho < 1$, the numbers

$$p_{0,n} = (1-\rho)^{(\lambda/\gamma)+1}\frac{\rho^n}{n!\,\gamma^n}\prod_{i=0}^{n-1}(\lambda+i\gamma), \qquad p_{1,n} = (1-\rho)^{(\lambda/\gamma)+1}\frac{\rho^{n+1}}{n!\,\gamma^n}\prod_{i=1}^{n}(\lambda+i\gamma)$$

form a steady-state solution of (3.47)–(3.49), and every steady-state solution equals them.

Milestones on the retrial queue

The generating functions satisfy (3.50)–(3.52) on $(-1,1)$, including the separable equation

$$P_0'(z) = \frac{\lambda\rho}{\gamma(1-\rho z)}P_0(z),$$

their closed form is (3.55),

$$P_0(z) = (1-\rho z)\left(\frac{1-\rho}{1-\rho z}\right)^{(\lambda/\gamma)+1}, \qquad P_1(z) = \rho\left(\frac{1-\rho}{1-\rho z}\right)^{(\lambda/\gamma)+1},$$

and the mean orbit size is (3.58), $L_o = \frac{\rho^2}{1-\rho}\cdot\frac{\mu+\gamma}{\gamma}$.

Milestones from the rest of Chapter 3

The bulk-input generating function (3.3), $p_0 = 1 - \rho$ with $\rho = \lambda\mathrm E[X]/\mu$, and the mean (3.4); the unique root $r_0 \in (0,1)$ of $\mu r^{K+1} - (\lambda+\mu)r + \lambda = 0$ and the geometric law $p_n = (1-r_0)r_0^n$ (3.9); and Cobham's formula (3.41)/(3.43), the unique solution of the linear system (3.40).

Significance

The closed form (3.57) makes every performance measure of the $M/M/1$ retrial queue explicit. The server is busy a fraction $\rho$ of the time, exactly as without retrials. The mean orbit size (3.58) is the $M/M/1$ mean queue length multiplied by $(\mu+\gamma)/\gamma$, and the mean time in orbit (3.59) follows from Little's law. These formulas quantify the cost of retrials against an ordinary queue and are the reference case against which approximations for multi-server retrial systems are checked.

The results are classical and proved in the book, partly through exercises (Problems 3.39–3.41). None of them is formalized in any proof assistant, as far as the platform's catalogue shows: there is no retrial, bulk or priority queue on Prove2Me. The mission produces machine-checked statements and, once solved, proofs of the chapter's main closed forms. It also produces a small reusable layer: generating functions of probability sequences on the closed unit disc, and the "probability solution of the balance equations" pattern for chains with countable state spaces.

Difficulty

The derivation in the book is formal. It differentiates power series term by term, divides by $1 - z$, integrates $\ln P_0$, and fixes the constant by setting $z = 1$, without justifying any of these steps. A formal proof has to show that the series converge and are differentiable on $(-1,1)$, that the differential equation determines $P_0$ up to a constant, and that the values at $z = 1$ are the limits of the values inside the disc (Abel's theorem). The uniqueness half of the goal is the hardest part. The book never proves it; it follows from the ODE argument only once every step is shown to hold for an arbitrary probability solution. Verifying that (3.57) solves (3.47)–(3.49) is only the easy half. The same pattern recurs in the bulk-input queue, where $z = 1$ is a removable singularity of (3.3). In the bulk-service queue the root $r_0$ is only characterized as the unique root in $(0,1)$, so existence and uniqueness of the root are part of the claim.

Formalization scope

A steady-state solution is a pair p0 p1 : ℕ → ℝ (resp. one sequence p : ℕ → ℝ) that is pointwise nonnegative, has total mass $1$ as a HasSum, and solves the book's balance equations exactly as printed, global balance and not detailed balance. Every "the steady-state solution is X" is stated with both halves: X is a steady-state solution, and every steady-state solution equals X. Stating only that (3.57) solves (3.47)–(3.49), without normalization or uniqueness, would be a trivializing formalization. So would taking $r_0$ as a given root in (3.9), or taking the $W_q^{(i)}$ in (3.41) as numbers assumed to satisfy it. None of these is used. The closed forms instantiated are (3.52), (3.55), (3.57), (3.58), (3.3), (3.4), (3.9), (3.41) and (3.43), each written out in full, with the real power $(1-\rho)^{(\lambda/\gamma)+1}$ as Real.rpow.

The conventions are as follows. The retrial generating functions take real arguments, on $(-1,1)$ for the differential equations and on $[-1,1]$ for the closed form. The bulk-input generating function takes complex arguments with $|z| \le 1$, $z \ne 1$, because (3.3) is $0/0$ at $z = 1$. The condition $\rho < 1$ is a hypothesis of every retrial statement. For the bulk-service queue the book's unnamed condition is stated as $\lambda < K\mu$. For bulk input, $\mathrm E[X] < \infty$ is assumed throughout, and the mean (3.4) is asserted under the further condition $\mathrm E[X^2] < \infty$, which it requires. For Cobham's formula only the algebraic content is formalized; the mean-value argument that yields (3.40) and (3.42) is not.

Needed infrastructure: power series of summable nonnegative sequences on the closed unit disc (convergence, term-by-term differentiation, Abel continuity), the binomial series $(1 - x)^{-a} = \sum_n \frac{a(a+1)\cdots(a+n-1)}{n!}x^n$ for real $a$, and uniqueness of invariant probability vectors for irreducible chains. All of this is reusable beyond the mission. Proofs of any milestone, of the easy half of the goal, or of the needed series facts are welcome contributions.

Selected references

• D. Gross, J. F. Shortle, J. M. Thompson, C. M. Harris, Fundamentals of Queueing Theory, 4th ed., Wiley, 2008, §§3.1, 3.2.0.1, 3.4.2, 3.5.1. https://doi.org/10.1002/9781118625651
• G. I. Falin, J. G. C. Templeton, Retrial Queues, Chapman & Hall, 1997. https://doi.org/10.1007/978-1-4899-2977-8
• J. R. Artalejo, Accessible bibliography on retrial queues, Mathematical and Computer Modelling 30 (1999) 1–6. https://doi.org/10.1016/S0895-7177(99)00128-4
• A. Cobham, Priority assignment in waiting line problems, Journal of the Operations Research Society of America 2 (1954) 70–76. https://doi.org/10.1287/opre.2.1.70

Motivation

Networks of queues model systems in which a job visits several service stations in turn: jobs in a computer system alternating between CPU and disks, machines cycling between operation and repair, parts routed through a job shop. In a closed network no job enters or leaves; a fixed population of $N$ customers circulates among $k$ nodes. Closed networks are the standard model of multiprogrammed computer systems and of machine-repair and finite-source systems, and they are the setting of chapter 4 of Gross, Shortle, Thompson and Harris, Fundamentals of Queueing Theory (4th ed., Wiley 2008, doi:10.1002/9781118625651).

The chapter's results form a short line of computational ideas:

• Jackson (1957, 1963) showed that open networks of exponential servers with Markovian routing have a product-form steady state; Gordon and Newell (1967) gave the closed-network version, (4.15)–(4.18) of the book.
• Buzen (1973) gave a convolution recursion for the normalizing constant $G(N)$ and for marginal distributions, (4.19)–(4.22).
• Reiser and Lavenberg (1980) introduced mean-value analysis (MVA), which computes mean queue lengths, waiting times and throughputs population by population without ever forming $G(N)$, (4.23)–(4.25); the book presents it following Bruell and Balbo (1980).
• The book closes the section with a recursion for the full marginal distributions, (4.26), which it proves from the product form (pp.207–209).

This mission formalizes that line, ending at (4.26).

Setting

A closed Jackson network has nodes $i = 1, \dots, k$, each with a single server whose service times are exponential with rate $\mu_i > 0$. A customer finishing service at node $i$ moves to node $j$ with probability $r_{ij}$; the routing matrix $R = (r_{ij})$ has nonnegative entries and rows summing to one, and it is irreducible: every node can be reached from every other. The state is $\bar n = (n_1, \dots, n_k)$, the number of customers at each node, with $n_1 + \cdots + n_k = N$; this state space is finite.

The steady-state distribution $p_{\bar n}$ is the probability vector on the state space that solves the flow-balance equations (4.14),

$$\sum_{j=1}^{k}\sum_{\substack{i=1\\ i\ne j}}^{k} \mu_i r_{ij}\, p_{\bar n;i^+j^-} = \sum_{i=1}^{k}\mu_i(1-r_{ii})\,p_{\bar n},$$

where $\bar n;i^+j^-$ has one more customer at $i$ and one fewer at $j$, and terms with a negative subscript or with $\mu_i$ at an empty node vanish. The traffic equations (4.16) are $\mu_i\rho_i = \sum_j \mu_j r_{ji}\rho_j$; they determine $\rho = (\rho_1, \dots, \rho_k)$ up to a positive factor. The normalizing constant is

$$G(N) = \sum_{n_1+\cdots+n_k=N}\rho_1^{n_1}\cdots\rho_k^{n_k},$$

and more generally, with $f_i(n) = \rho_i^{\,n}/a_i(n)$ for $c_i$-server nodes ((4.13)), $G(N) = \sum \prod_i f_i(n_i)$ and Buzen's function $g_m(n) = \sum_{n_1+\cdots+n_m=n}\prod_{i\le m} f_i(n_i)$.

For each population $N$ write $p_i(n, N) = \Pr\{N_i = n\}$ for the marginal distribution at node $i$, $\bar P_i(n; N) = \Pr\{N_i \ge n\}$, $L_i(N)$ for the mean number at node $i$, and

$$\lambda_i(N) = \Pr\{\text{server busy at node } i\}\cdot\mu_i$$

for the throughput of node $i$.

Formalization targets

Goal: the marginal recursion (4.26)

For every node $i$,

$$p_i(0,0) = 1, \qquad p_i(n, N) = \frac{\lambda_i(N)}{\mu_i}\,p_i(n-1, N-1) \quad (n, N \ge 1).$$

It involves only the steady-state distributions and quantities computed from them; it holds for every irreducible routing matrix and every choice of rates.

Milestones

1. Product form (4.14)–(4.16). For any positive solution $\rho$ of (4.16), a probability distribution solves (4.14) if and only if $p_{\bar n} = G(N)^{-1}\rho_1^{n_1}\cdots\rho_k^{n_k}$.
2. Buzen's algorithm (4.19)–(4.21). $G(N) = g_k(N)$, $g_m(n) = \sum_{i=0}^{n} f_m(i)\,g_{m-1}(n-i)$, $g_1 = f_1$, $g_m(0) = 1$.
3. Marginal at the last node (4.22). $p_k(n) = f_k(n)\,g_{k-1}(N-n)/G(N)$ for $0 \le n \le N$.
4. Complementary marginal (p.208). $\bar P_i(n_i; N) = \rho_i^{\,n_i}G(N-n_i)/G(N)$.
5. Mean-value analysis (4.23)–(4.25). $L_i(0) = 0$; $L_i(N) = \lambda_i(N)W_i(N)$ with $W_i(N) = (1 + L_i(N-1))/\mu_i$; and for $v$ solving $v_i = \sum_j v_j r_{ji}$ with $v_l = 1$, $\lambda_l(N) = N/\sum_i v_iW_i(N)$ and $\lambda_i(N) = \lambda_l(N)v_i$.

Significance

The product form reduces a $\binom{N+k-1}{N}$-state Markov chain to the constants $G(0), \dots, G(N)$, and Buzen's recursion computes them in $O(kN^2)$ operations. Mean-value analysis goes further and avoids $G(N)$, whose magnitude can overflow or underflow for large populations; it is the method used in capacity planning of computer systems. The recursion (4.26) extends MVA from means to full marginal distributions, so a single pass over $N$ yields every nodal distribution.

All of these results are classical and proved in the literature; the book proves (4.26) itself. What the mission adds is a machine-checked development of them from the global balance equations: the product form with its uniqueness, the convolution identities, the marginal formulas, and the correctness of the MVA iteration as stated by the book, all over one shared definition layer. A search of the platform on 2026-09-28 found no formal statement of Buzen's algorithm or of MVA. The platform has Kelly's closed migration process theorem (KellyStochasticNetworks.closed_migration_equilibrium), which shows that the unnormalized product form satisfies the equilibrium equations under Kelly's conventions; the normalization, uniqueness and everything downstream of the product form are new here.

Difficulty

The combinatorial identities (Buzen's recursion, the tail marginal) are reindexings of finite sums over compositions of $N$; in Lean the work is in bijections between the state spaces $\{n_1+\cdots+n_k = N\}$ for different $k$ and $N$. The substantive step is uniqueness in the product-form theorem: the global balance equations have a one-dimensional solution space only because the chain on the $N$-customer states is irreducible on the population level, which is a property of the network chain and not of the routing matrix alone. The goal and MVA also need a positive solution of the traffic equations, which is not among the hypotheses and has to come from irreducibility of $R$. The book's own intuitive derivation of MVA via the arrival theorem is not the route the statements require; they are stated in terms of the steady-state distributions alone.

Formalization scope

Nodes are Fin k (book node $i$ is index $i-1$); states are n : Fin k → ℕ with ∑ i, n i = N, collected in a Finset, and all sums are finite. A distribution is a real function on $\mathbb N^k$ that is nonnegative, vanishes off the $N$-customer states and sums to one there. The balance equations are (4.14) verbatim with the book's boundary convention (p.188), not detailed balance. All results except Buzen's algorithm and (4.22) are for single-server nodes, as in the book; (4.13)'s multiserver factor $a_i(n)$ enters only (4.19)–(4.22).

Closed forms instantiated in the statements: the product form $G(N)^{-1}\prod_i\rho_i^{n_i}$ ((4.15)); $G(N)$ as the explicit sum (4.18)/(4.19); $a_i(n)$ from (4.13); $g_m$ from (4.20); $p_k(n) = f_k(n)g_{k-1}(N-n)/G(N)$ ((4.22)); $\bar P_i(n;N) = \rho_i^nG(N-n)/G(N)$ (p.208); $W_i(N) = (1+L_i(N-1))/\mu_i$ ((4.23)); $\lambda_l(N) = N/\sum_i v_iW_i(N)$ (MVA step (iii)(b)).

Two trivializing formalizations are ruled out: $\lambda_i(N)$ in (4.26) and (4.24) is the throughput computed from the steady-state distribution, not a free constant (which would make (4.26) a definition); and $g_m$ is defined by the sum (4.20), so the recursion (4.21) is a theorem rather than rfl. The product-form statement is an equivalence, so it asserts both that the product form is a steady state and that it is the only one.

Needed infrastructure: bijections between compositions of $N$ into $k$ and $k-1$ parts, uniqueness of stationary distributions of irreducible finite continuous-time chains (stated directly via the balance equations), and existence of positive solutions of $v = vR$ for irreducible stochastic $R$. The last two are reusable beyond this mission. Contributions welcome: proofs of the milestones in any order, and helper lemmas on these three points.

Not formalized: open Jackson networks (4.11) and Burke's theorem (4.5)–(4.6), multiclass networks (§4.2.1), the multiserver recursion (4.27) and cyclic queues (§4.4).

Selected references

• D. Gross, J. F. Shortle, J. M. Thompson, C. M. Harris, Fundamentals of Queueing Theory, 4th ed., Wiley, 2008, §4.3, pp.195–209. https://doi.org/10.1002/9781118625651
• J. R. Jackson, "Jobshop-like queueing systems", Management Science 10(1), 1963. https://doi.org/10.1287/mnsc.10.1.131
• W. J. Gordon, G. F. Newell, "Closed queuing systems with exponential servers", Operations Research 15(2), 1967. https://doi.org/10.1287/opre.15.2.254
• J. P. Buzen, "Computational algorithms for closed queueing networks with exponential servers", Communications of the ACM 16(9), 1973. https://doi.org/10.1145/362342.362345
• M. Reiser, S. S. Lavenberg, "Mean-value analysis of closed multichain queuing networks", Journal of the ACM 27(2), 1980. https://doi.org/10.1145/322186.322195
• S. C. Bruell, G. Balbo, Computational Algorithms for Closed Queueing Networks, North-Holland, 1980.