M/M/c Queue Model: Analysis & Extensions

• M/M/c Queue Model is a continuous-time, Markovian multi-server queue characterized by Poisson arrivals, exponential service times, and Erlang’s C formula.
• The model enables closed-form computation of key metrics such as server utilization, waiting probability, and mean sojourn time, supporting robust system design.
• Extensions include priority schemes, state-dependent service, and extreme-value analysis, broadening its application in telecommunications and operational research.

The $M/M/c$ queue is a classical continuous-time, Markovian multi-server queueing model central to queueing theory, performance analysis, and operations research. Defined by Poisson arrivals, $c$ parallel servers, and independent exponential service times, it underlies a wide range of performance models in e.g. telecommunications, call centers, and service systems. Recent research continues to extend, compute, and apply its probabilistic properties: from joint queue-length distributions and extreme-value statistics to busy-period durations, multi-priority systems, and state-dependent service mechanisms.

1. Formal Model Description

Let arrivals form a Poisson process of rate $\lambda$. There are $c$ identical servers, each independently serving one customer at exponential rate $\mu$. The system allows infinite waiting room; arrivals that find fewer than $c$ customers in service immediately begin service, while those finding all $c$ servers occupied join a first-come-first-served backlog queue. Stability requires $\rho := \lambda/(c\mu) < 1$. The system state is $X(t) = n$ ($0 \le n < \infty$), the total number of customers present. The process is birth–death, with rates: $c$0

$c$1

Balance equations yield the steady-state probabilities $c$2. The normalization is: $c$3 The explicit forms are: $c$4 The probability that an arrival must wait (all servers busy) is: $c$5 This is Erlang’s $c$6 formula and is foundational across contexts requiring multi-server delay analysis.

2. Key Performance Metrics and First-Order Properties

The principal performance metrics and their closed forms are as follows:

Measure	Expression	Comments
Utilization per server	$c$7	$c$8 for stability
Mean number in queue ($c$9)	$\lambda$0	Excludes those in service
Mean queueing time ($\lambda$1)	$\lambda$2	Little’s law
Mean system size ($\lambda$3)	$\lambda$4	$\lambda$5 plus those in service
Mean sojourn time ($\lambda$6)	$\lambda$7	Includes service time
Probability all servers idle ($\lambda$8)	$\lambda$9	Explicit sum

These explicit steady-state forms enable closed-form, numerically stable computation for arbitrary $c$0, rendering the $c$1 queue tractable for analytic, design, and simulation contexts (Finch, 2019).

3. Distributional and Extreme-Value Properties

Beyond stationary means, recent work has addressed the high-level behavior of the queue, particularly the distribution of the maximum queue length over a finite interval. Let $c$2 (excluding those in service). The Poisson clumping heuristic provides an accurate asymptotic approximation for this extreme-value distribution.

Let $c$3 and define $c$4. The probability that the maximum does not exceed $c$5 is, for large $c$6,

$c$7

where

$c$8

Higher-order (finite-$c$9) corrections introduce an additional bracketed term inside the exponential, improving the fit for moderate $\mu$0 ($\mu$1–$\mu$2). The limit law is a discrete Gumbel and governs extremes for general $\mu$3 (Finch, 2019). Mean and variance estimates are given by: $\mu$4

$\mu$5

These expressions match simulation results to a high degree of accuracy.

4. Extensions: Priority Classes and State-Dependent Service

Priority Extensions

Multi-class and multi-priority $\mu$6 models have been extended to both non-preemptive and preemptive regimes.

Non-preemptive priority with $\mu$7 classes: Poisson arrivals at rates $\mu$8, exponential services, $\mu$9 servers. Zuk and Kirszenblat (Zuk et al., 2023) derive explicit forms for the joint queue-length generating function $c$0 as a product of root-ratio terms, and invert it via multidimensional FFT to recover the joint stationary distribution. The fixed-point iteration for the vectorized stationary balance equations is robust and probability-preserving, with convergence slowed only as $c$1 or $c$2 increases. The methodology yields all joint and marginal performance measures and has been numerically validated to high precision.

Preemptive two-class priority: The continuous-time Markov structure yields Laplace transforms for all finite-time transition probabilities via spectral methods, including explicit closed forms for the vertical boundary (states with $c$3 high-priority customers). The general case is recursively tractable with Ramaswami’s formula adapted from $c$4 Markov chains (Selen et al., 2016).

State-Dependent and Time-Dependent Service

Service rates affected by experienced waiting time are considered in (D'Auria et al., 2021): each arrival is "tagged" by whether its virtual waiting time at entry exceeds a threshold $c$5, and is served at either rate $c$6 (if $c$7) or $c$8 (if $c$9). The model is formulated via a Markov process tracking the virtual queueing time and server occupancy, matched to an embedded process with piecewise ODE structure. The stationary distribution is characterized recursively and in closed form as a mixture of exponentials for the waiting time, with rates determined as zeros of characteristic polynomials tied to the underlying Markov dynamics. When $c$0, performance can change drastically, with slowdown ($c$1) causing heavy-tailed delays.

5. Partial Busy Periods and Time-Dependent Behavior

The partial busy period ($c$2): the time from system entry by an arrival (when queue was empty) until the system empties again, is analyzed explicitly for general $c$3 in (Zuk et al., 2023). For $c$4 the result connects to classic Karlin–McGregor integral representations. For arbitrary $c$5, two main techniques are provided:

• Spectral method: The Laplace transform $c$6 is written in terms of resolvent sums over eigenvalues of the finite-dimensional generator, leading to a survival function as a sum of exponentials (from simple poles) plus a branch cut integral.
• Algebraic method: The closed-form representation involves a sum over Marcum Q-functions, with parameters and coefficients arising from roots of a uniquely defined "cut polynomial" of degree $c$7.

Large-$c$8 asymptotics distinguish two regimes: fixed-traffic intensity ($c$9) gives a mixture of empty and exponentially-distributed periods; the $\rho := \lambda/(c\mu) < 1$0 scaling yields mixtures tied to zeros of Kummer functions.

Efficient computation leverages eigenvalue routines (e.g., LAPACK’s {\tt dstevd}), Chebyshev quadrature, and FFT-based polynomial root solvers. These approaches allow full distributional computation for $\rho := \lambda/(c\mu) < 1$1 up to $\rho := \lambda/(c\mu) < 1$2.

6. Analytical and Computational Methods

Explicit forms for steady-state and transient distributions are available in the classical case; closed-form, efficiently computable representations for time-dependent, multi-class, and state-dependent extensions require:

• Reduction to birth–death or block-tridiagonal Markov chains;
• Laplace transform (spectral) analysis, inversion via residue and cut contours;
• Explicit, finite-sum representations (Marcum-Q, Bessel series, or mixture of exponentials) using parameter polynomials derived from model structure;
• Fixed-point or matrix-analytic recursions (Ramaswami-like algorithms) for balance equations in high-dimensional spaces;
• Accelerated numerical inversion (FFT-based or Gauss–Chebyshev quadrature) for generating functions and transforms, maintaining probability normalization and exploiting structure for tractability;
• Direct simulation for empirical validation, as in the extreme-value study of (Finch, 2019).

Accuracy, robustness, and computational efficiency are achieved up to $\rho := \lambda/(c\mu) < 1$3 and, for priority models, $\rho := \lambda/(c\mu) < 1$4–$\rho := \lambda/(c\mu) < 1$5.

7. Applications and Implications

The $\rho := \lambda/(c\mu) < 1$6 model and its extensions directly address:

• Performance guarantees in service systems (delay probabilities, maximal queue length, mean waiting times);
• Regenerative simulation interval estimation and confidence procedures (busy period distributions) (Zuk et al., 2023);
• Service system design; e.g., setting $\rho := \lambda/(c\mu) < 1$7 to meet SLA thresholds (via Erlang $\rho := \lambda/(c\mu) < 1$8 or joint queue-length measures);
• Priority queueing in telecommunications, health care, or call centers (preemptive vs non-preemptive delay, fairness, overflow probabilities);
• Queueing-time dependent dynamical control policies and their impact.

The collective body of work reaffirms the tractability and flexibility of the $\rho := \lambda/(c\mu) < 1$9 queue for both theoretical and applied queueing analysis, while delineating the frontiers of what is computable: joint queue-lengths for many classes require FFT inversion; busy periods demand spectral tools; maxima require clumping heuristics. Generalizations for non-Markovian service and arrival, reneging, or blocking remain challenging.

8. Context and Future Research

Recent advances yield the most explicit known formulas for joint queue-lengths (Zuk et al., 2023), partial busy periods (Zuk et al., 2023), and time-dependent multi-class probabilities (Selen et al., 2016) for $X(t) = n$0 systems. These results fill longstanding gaps, especially for $X(t) = n$1 and general $X(t) = n$2.

A plausible implication is that similar algebraic–spectral hybrid techniques could be adapted to $X(t) = n$3 or $X(t) = n$4 queues, though the loss of Markov structure presents substantial obstacles. As system parameters scale (large $X(t) = n$5, near-critical $X(t) = n$6), asymptotic analysis and numerically stable algorithms remain an ongoing research frontier.

The $X(t) = n$7 queue retains a central role in queueing theory both as a reference model and as a basis for extendable, analytically tractable specializations in modern stochastic network analysis.