Queueing Systems

Single-station queueing system analysis.

The qsys module provides exact and approximate formulas for single queueing systems such as M/M/c, M/G/1, G/G/1, and their variants.

Key function categories:

Native Python implementations for queueing system analysis.

This module provides pure Python/NumPy implementations for analyzing single queueing systems, including basic queues (M/M/1, M/M/k, M/G/1), G/G/1 approximations, MAP-based queues, and scheduling disciplines.

Key algorithms:

Basic queues: qsys_mm1, qsys_mmk, qsys_mg1, qsys_gm1, qsys_mminf, qsys_mginf G/G/1 approximations: Allen-Cunneen, Kingman, Marchal, Whitt, Heyman, etc. G/G/k approximations: qsys_gigk_approx MAP/D queues: qsys_mapdc, qsys_mapd1 MAP/PH queues: qsys_phph1, qsys_mapph1, qsys_mapm1, qsys_mapmc, qsys_mapmap1 Scheduling: qsys_mg1_prio, qsys_mg1_srpt, qsys_mg1_fb, etc. Loss systems: qsys_mm1k_loss, qsys_mg1k_loss, qsys_mxm1

qsys_mapdc(D0, D1, s, c, max_num_comp=1000, num_steps=1, verbose=0)[source]

Analyze MAP/D/c queue (MAP arrivals, deterministic service, c servers).

Uses Non-Skip-Free (NSF) Markov chain analysis embedding at deterministic service intervals. Multiple arrivals can occur per interval.

Parameters:

  • D0 (numpy.ndarray) – MAP hidden transition matrix (n x n).

  • D1 (numpy.ndarray) – MAP arrival transition matrix (n x n).

  • s (float) – Deterministic service time (positive scalar).

  • c (int) – Number of servers.

  • max_num_comp (int) – Maximum number of queue length components (default 1000).

  • num_steps (int) – Number of waiting time distribution points per interval (default 1).

  • verbose (int) – Verbosity level (default 0).

Returns:

Performance metrics including:

  • mean_queue_length: Mean number of customers in system

  • mean_waiting_time: Mean waiting time in queue

  • mean_sojourn_time: Mean sojourn time (waiting + service)

  • utilization: Server utilization (per server)

  • queue_length_dist: Queue length distribution P(Q=n)

  • waiting_time_dist: Waiting time CDF at discrete points

  • analyzer: Analyzer identifier

Return type:

dict

qsys_mapd1(D0, D1, s, max_num_comp=1000, num_steps=1)[source]

Analyze MAP/D/1 queue (single-server convenience function).

Parameters:

  • D0 (numpy.ndarray) – MAP hidden transition matrix.

  • D1 (numpy.ndarray) – MAP arrival transition matrix.

  • s (float) – Deterministic service time.

  • max_num_comp (int) – Maximum number of queue length components.

  • num_steps (int) – Number of waiting time points per interval.

Returns:

Performance metrics (see qsys_mapdc).

Return type:

dict

qsys_mm1(lambda_val, mu)[source]

Analyze M/M/1 queue (Poisson arrivals, exponential service).

Parameters:

  • lambda_val (float) – Arrival rate (lambda)

  • mu (float) – Service rate

Returns:

Performance measures including:

  • L: Mean number in system

  • Lq: Mean number in queue

  • W: Mean response time (time in system)

  • Wq: Mean waiting time (time in queue)

  • rho: Utilization (lambda/mu)

Return type:

dict

Example

>>> result = qsys_mm1(0.5, 1.0)
>>> print(f"Utilization: {result['rho']:.2f}")
Utilization: 0.50
qsys_mmk(lambda_val, mu, k)[source]

Analyze M/M/k queue (Poisson arrivals, k exponential servers).

Parameters:

  • lambda_val (float) – Arrival rate (lambda)

  • mu (float) – Service rate per server

  • k (int) – Number of parallel servers

Returns:

Performance measures including:

  • L: Mean number in system

  • Lq: Mean number in queue

  • W: Mean response time

  • Wq: Mean waiting time

  • rho: Utilization per server (lambda/(k*mu))

  • P0: Probability of empty system

Return type:

dict

Example

>>> result = qsys_mmk(2.0, 1.0, 3)
>>> print(f"Utilization: {result['rho']:.2f}")
Utilization: 0.67
qsys_mg1(lambda_val, mu, cs)[source]

Analyze M/G/1 queue using Pollaczek-Khinchine formula.

Parameters:

  • lambda_val (float) – Arrival rate (lambda)

  • mu (float) – Service rate (mean service time = 1/mu)

  • cs (float) – Coefficient of variation of service time (std/mean)

Returns:

Performance measures including:

  • L: Mean number in system

  • Lq: Mean number in queue

  • W: Mean response time

  • Wq: Mean waiting time

  • rho: Utilization (lambda/mu)

Return type:

dict

Example

>>> result = qsys_mg1(0.5, 1.0, 1.0)  # cs=1 is exponential (M/M/1)
qsys_gm1(lambda_val, mu, ca)[source]

Analyze G/M/1 queue (general arrivals, exponential service).

Parameters:

  • lambda_val (float) – Arrival rate (lambda)

  • mu (float) – Service rate

  • ca (float) – Coefficient of variation of inter-arrival time

Returns:

Performance measures including:

  • L: Mean number in system

  • Lq: Mean number in queue

  • W: Mean response time

  • Wq: Mean waiting time

  • rho: Utilization

Return type:

dict

Note

Uses approximation based on the coefficient of variation.

qsys_mminf(lambda_val, mu)[source]

Analyze M/M/inf queue (infinite servers / delay station).

Parameters:

  • lambda_val (float) – Arrival rate (lambda)

  • mu (float) – Service rate

Returns:

Performance measures including:

  • L: Mean number in system (= lambda/mu)

  • Lq: Mean number in queue (= 0)

  • W: Mean time in system (= 1/mu)

  • Wq: Mean waiting time (= 0)

  • P0: Probability of empty system

Return type:

dict

qsys_mginf(lambda_val, mu, k=None)[source]

Analyze M/G/inf queue (infinite servers, general service).

Performance is independent of service time distribution shape. Number of customers follows Poisson distribution.

Parameters:

  • lambda_val (float) – Arrival rate (lambda)

  • mu (float) – Service rate (mean service time = 1/mu)

  • k (int | None) – Optional state for probability computation

Returns:

Performance measures including:

  • L: Mean number in system

  • Lq: Mean number in queue (= 0)

  • W: Mean time in system (= 1/mu)

  • Wq: Mean waiting time (= 0)

  • P0: Probability of empty system

  • Pk: Probability of k customers (if k provided)

Return type:

dict

qsys_mmcc_retrial_fp(lambda_val, mu, c, tol=1e-10, maxiter=10000)[source]

Fixed-point approximation for M/M/c/c retrial queue.

Customers arrive at rate lambda to a system with c servers (no waiting room), each with service rate mu. Blocked customers join an orbit and retry. Under the assumption that the retrial rate is small relative to the service rate, the total arrival flow (fresh + retrial) is approximated by a Poisson process with rate lambda + r, where r satisfies the fixed-point equation:

r = (lambda + r) * B((lambda + r) / mu, c)

and B(a, c) is the Erlang-B blocking probability for offered load a and c servers.

Parameters:

  • lambda_val (float) – Arrival rate

  • mu (float) – Service rate per server

  • c (int) – Number of servers (= capacity, no waiting room)

  • tol (float) – Convergence tolerance (default: 1e-10)

  • maxiter (int) – Maximum iterations (default: 10000)

Returns:

Performance measures including:

  • blocProb: Blocking probability

  • r: Additional arrival rate due to retrials

  • niter: Number of iterations to converge

  • rho: Offered load (lambda / (c * mu))

  • L: Mean number of busy servers

Return type:

dict

References

Cohen (1957), fixed-point approximation for M/M/c/c retrial queues. Phung-Duc, “Retrial Queueing Models: A Survey on Theory and Applications”, 2019, Eq. (1).

Example

>>> result = qsys_mmcc_retrial_fp(2.0, 1.0, 3)
>>> print(f"Blocking: {result['blocProb']:.4f}")
qsys_gig1_approx_allencunneen(lambda_val, mu, ca, cs)[source]

Allen-Cunneen approximation for G/G/1 queue.

Matches MATLAB qsys_gig1_approx_allencunneen.m exactly.

Parameters:

  • lambda_val (float) – Arrival rate

  • mu (float) – Service rate

  • ca (float) – Coefficient of variation of inter-arrival time

  • cs (float) – Coefficient of variation of service time

Returns:

W: Mean response time (time in system) rhohat: Effective utilization (so that M/M/1 formulas still hold)

Return type:

Tuple of (W, rhohat) where

qsys_gig1_approx_kingman(lambda_val, mu, ca, cs)[source]

Kingman’s upper bound approximation for G/G/1 queue.

Note: In MATLAB, ‘gig1’ and ‘gig1.kingman’ map to qsys_gig1_ubnd_kingman. This function provides the same formula.

Parameters:

  • lambda_val (float) – Arrival rate

  • mu (float) – Service rate

  • ca (float) – Coefficient of variation of inter-arrival time

  • cs (float) – Coefficient of variation of service time

Returns:

W: Mean response time (time in system) rhohat: Effective utilization (so that M/M/1 formulas still hold)

Return type:

Tuple of (W, rhohat) where

qsys_gig1_approx_marchal(lambda_val, mu, ca, cs)[source]

Marchal’s approximation for G/G/1 queue.

Matches MATLAB qsys_gig1_approx_marchal.m exactly. Note: MATLAB formula uses ca (not ca^2) in the numerator factor.

Parameters:

  • lambda_val (float) – Arrival rate

  • mu (float) – Service rate

  • ca (float) – Coefficient of variation of inter-arrival time

  • cs (float) – Coefficient of variation of service time

Returns:

W: Mean response time (time in system) rhohat: Effective utilization (so that M/M/1 formulas still hold)

Return type:

Tuple of (W, rhohat) where

qsys_gig1_approx_whitt(lambda_val, mu, ca, cs)[source]

Whitt’s approximation for G/G/1 queue.

Uses QNA (Queueing Network Analyzer) approximation. Note: No direct MATLAB counterpart (qsys_gig1_approx_whitt.m does not exist).

Parameters:

  • lambda_val (float) – Arrival rate

  • mu (float) – Service rate

  • ca (float) – Coefficient of variation of inter-arrival time

  • cs (float) – Coefficient of variation of service time

Returns:

W: Mean response time (time in system) rhohat: Effective utilization (so that M/M/1 formulas still hold)

Return type:

Tuple of (W, rhohat) where

qsys_gig1_approx_heyman(lambda_val, mu, ca, cs)[source]

Heyman’s approximation for G/G/1 queue.

Matches MATLAB qsys_gig1_approx_heyman.m exactly.

Parameters:

  • lambda_val (float) – Arrival rate

  • mu (float) – Service rate

  • ca (float) – Coefficient of variation of inter-arrival time

  • cs (float) – Coefficient of variation of service time

Returns:

W: Mean response time (time in system) rhohat: Effective utilization (so that M/M/1 formulas still hold)

Return type:

Tuple of (W, rhohat) where

qsys_gig1_approx_kobayashi(lambda_val, mu, ca, cs)[source]

Kobayashi’s approximation for G/G/1 queue.

Matches MATLAB qsys_gig1_approx_kobayashi.m exactly.

Parameters:

  • lambda_val (float) – Arrival rate

  • mu (float) – Service rate

  • ca (float) – Coefficient of variation of inter-arrival time

  • cs (float) – Coefficient of variation of service time

Returns:

W: Mean response time (time in system) rhohat: Effective utilization

Return type:

Tuple of (W, rhohat) where

qsys_gig1_approx_klb(lambda_val, mu, ca, cs)[source]

Kraemer-Langenbach-Belz (KLB) approximation for G/G/1 queue.

Matches MATLAB qsys_gig1_approx_klb.m exactly.

Parameters:

  • lambda_val (float) – Arrival rate

  • mu (float) – Service rate

  • ca (float) – Coefficient of variation of inter-arrival time

  • cs (float) – Coefficient of variation of service time

Returns:

W: Mean response time (time in system) rhohat: Effective utilization (so that M/M/1 formulas still hold)

Return type:

Tuple of (W, rhohat) where

qsys_gig1_approx_gelenbe(lambda_val, mu, ca, cs)[source]

Gelenbe’s approximation for G/G/1 queue.

Matches MATLAB qsys_gig1_approx_gelenbe.m exactly. Note: MATLAB returns only W (no rhohat).

Parameters:

  • lambda_val (float) – Arrival rate

  • mu (float) – Service rate

  • ca (float) – Coefficient of variation of inter-arrival time

  • cs (float) – Coefficient of variation of service time

Returns:

Mean response time W

Return type:

float

qsys_gig1_approx_kimura(lambda_val, mu, ca, cs)[source]

Kimura’s approximation for G/G/1 queue.

Matches MATLAB qsys_gig1_approx_kimura.m exactly. Note: MATLAB uses sigma (=rho) as first parameter, returns only W.

Parameters:

  • lambda_val (float) – Arrival rate (used as sigma = lambda/mu = rho)

  • mu (float) – Service rate

  • ca (float) – Coefficient of variation of inter-arrival time

  • cs (float) – Coefficient of variation of service time

Returns:

Mean response time W

Return type:

float

qsys_gig1_approx_myskja(lambda_val, mu, ca, cs, q0, qa)[source]

Myskja’s approximation for G/G/1 queue.

Uses third moments of arrival and service distributions for improved accuracy.

Reference: Myskja, A. (1991). “An Experimental Study of a H₂/H₂/1 Queue”. Stochastic Models, 7(4), 571-595.

Parameters:

  • lambda_val (float) – Arrival rate

  • mu (float) – Service rate

  • ca (float) – Coefficient of variation of inter-arrival time

  • cs (float) – Coefficient of variation of service time

  • q0 (float) – Lowest relative third moment for given mean and SCV

  • qa (float) – Third relative moment E[X^3]/6/E[X]^3 of inter-arrival time

Returns:

Waiting time W

Return type:

float

qsys_gig1_approx_myskja2(lambda_val, mu, ca, cs, q0, qa)[source]

Modified Myskja (Myskja2) approximation for G/G/1 queue.

An improved version of Myskja’s approximation with better accuracy for a wider range of parameter combinations.

Reference: Myskja, A. (1991). “An Experimental Study of a H₂/H₂/1 Queue”. Stochastic Models, 7(4), 571-595.

Parameters:

  • lambda_val (float) – Arrival rate

  • mu (float) – Service rate

  • ca (float) – Coefficient of variation of inter-arrival time

  • cs (float) – Coefficient of variation of service time

  • q0 (float) – Lowest relative third moment for given mean and SCV

  • qa (float) – Third relative moment E[X^3]/6/E[X]^3 of inter-arrival time

Returns:

Waiting time W

Return type:

float

qsys_gg1(lambda_val, mu, ca2, cs2)[source]

G/G/1 queue analysis using exact methods for special cases and Allen-Cunneen approximation for the general case.

Parameters:

  • lambda_val (float) – Arrival rate

  • mu (float) – Service rate

  • ca2 (float) – Squared coefficient of variation of inter-arrival time

  • cs2 (float) – Squared coefficient of variation of service time

Returns:

W: Mean response time (time in system) rhohat: Effective utilization

Return type:

Tuple of (W, rhohat) where

References

Original JAR: jar/src/main/kotlin/jline/api/qsys/Qsys_gg1.kt

qsys_gig1_lbnd(lambda_val, mu, ca, cs)[source]

Fundamental theoretical lower bounds for G/G/1 queues.

These are the minimum possible values that performance measures cannot fall below for any realization of the arrival and service processes.

Parameters:

  • lambda_val (float) – Arrival rate

  • mu (float) – Service rate

  • ca (float) – Coefficient of variation of inter-arrival time

  • cs (float) – Coefficient of variation of service time

Returns:

W: Lower bound on mean response time (= 1/mu) rhohat: Effective utilization

Return type:

Tuple of (W, rhohat) where

References

Original JAR: jar/src/main/kotlin/jline/api/qsys/Qsys_gig1_lbnd.kt

qsys_gigk_approx(lambda_val, mu, ca, cs, k)[source]

Approximation for G/G/k queue.

Matches MATLAB qsys_gigk_approx.m formula using alpha-factor correction.

Parameters:

  • lambda_val (float) – Arrival rate

  • mu (float) – Service rate per server

  • ca (float) – Coefficient of variation of inter-arrival time

  • cs (float) – Coefficient of variation of service time

  • k (int) – Number of servers

Returns:

W: Approximate mean response time rhohat: Effective utilization

Return type:

Tuple of (W, rhohat) where

qsys_gigk_approx_cosmetatos(lambda_val, mu, ca, cs, k)[source]

G/G/k queue approximation using the Cosmetatos method.

Adjusts M/M/k results based on the variability of arrival and service processes.

Parameters:

  • lambda_val (float) – Arrival rate

  • mu (float) – Service rate per server

  • ca (float) – Coefficient of variation of inter-arrival time

  • cs (float) – Coefficient of variation of service time

  • k (int) – Number of servers

Returns:

W: Approximate mean response time rhohat: Effective utilization

Return type:

Tuple of (W, rhohat) where

References

Original JAR: jar/src/main/kotlin/jline/api/qsys/Qsys_gigk_approx_cosmetatos.kt

qsys_gigk_approx_whitt(lambda_val, mu, ca, cs, k)[source]

G/G/k queue approximation using Whitt’s QNA method.

Uses the Queue Network Analyzer (QNA) methodology to provide accurate estimates for multi-server queues with general distributions.

Parameters:

  • lambda_val (float) – Arrival rate

  • mu (float) – Service rate per server

  • ca (float) – Coefficient of variation of inter-arrival time

  • cs (float) – Coefficient of variation of service time

  • k (int) – Number of servers

Returns:

W: Approximate mean response time rhohat: Effective utilization

Return type:

Tuple of (W, rhohat) where

References

Original JAR: jar/src/main/kotlin/jline/api/qsys/Qsys_gigk_approx_whitt.kt

qsys_gig1_ubnd_kingman(lambda_val, mu, ca, cs)[source]

Kingman’s upper bound for G/G/1 queue waiting time.

This is an upper bound approximation that provides conservative estimates for the mean waiting time.

Parameters:

  • lambda_val (float) – Arrival rate

  • mu (float) – Service rate

  • ca (float) – Coefficient of variation of inter-arrival time

  • cs (float) – Coefficient of variation of service time

Returns:

W: Upper bound on mean response time rhohat: Effective utilization (so M/M/1 formulas still hold)

Return type:

Tuple of (W, rhohat) where

References

Original MATLAB: matlab/src/api/qsys/qsys_gig1_ubnd_kingman.m

qsys_gigk_approx_kingman(lambda_val, mu, ca, cs, k)[source]

Kingman’s approximation for G/G/k queue waiting time.

Extends Kingman’s approximation to multi-server queues using M/M/k waiting time as a base.

Parameters:

  • lambda_val (float) – Arrival rate

  • mu (float) – Service rate per server

  • k (int) – Number of servers

  • ca (float) – Coefficient of variation of inter-arrival time

  • cs (float) – Coefficient of variation of service time

Returns:

W: Approximate mean response time rhohat: Effective utilization

Return type:

Tuple of (W, rhohat) where

References

Original MATLAB: matlab/src/api/qsys/qsys_gigk_approx_kingman.m

qsys_mg1_prio(lambda_vec, mu_vec, cs_vec)[source]

Analyze M/G/1 queue with non-preemptive (Head-of-Line) priorities.

Matches MATLAB qsys_mg1_prio.m exactly.

Parameters:

  • lambda_vec (numpy.ndarray) – Vector of arrival rates per priority class (class 1 = highest)

  • mu_vec (numpy.ndarray) – Vector of service rates per priority class

  • cs_vec (numpy.ndarray) – Vector of coefficients of variation per priority class

Returns:

W: Vector of mean response times per priority class rho: System utilization (rhohat = Q/(1+Q) format)

Return type:

Tuple of (W, rho)

qsys_mg1_srpt(lambda_vec, mu_vec, cs_vec)[source]

Analyze M/G/1 queue with Shortest Remaining Processing Time (SRPT).

Matches MATLAB qsys_mg1_srpt.m exactly: - Exponential case: preemptive priority with cumulative residual E_R_k - General case: Schrage-Miller formula with numerical integration

Parameters:

  • lambda_vec (numpy.ndarray) – Vector of arrival rates per class

  • mu_vec (numpy.ndarray) – Vector of service rates per class

  • cs_vec (numpy.ndarray) – Vector of coefficients of variation per class

Returns:

W: Vector of mean response times per class rho: System utilization (rhohat format)

Return type:

Tuple of (W, rho)

qsys_mg1_fb(lambda_vec, mu_vec, cs_vec)[source]

Analyze M/G/1 queue with Foreground-Background (FB/LAS) scheduling.

Matches MATLAB qsys_mg1_fb.m exactly: - Exponential case: numerical integration of E[T(x)] * f_k(x) - General case: class-based approximation

Parameters:

  • lambda_vec (numpy.ndarray) – Vector of arrival rates per class

  • mu_vec (numpy.ndarray) – Vector of service rates per class

  • cs_vec (numpy.ndarray) – Vector of coefficients of variation per class

Returns:

W: Vector of mean response times per class rho: System utilization (rhohat format)

Return type:

Tuple of (W, rho)

qsys_mg1_lrpt(lambda_vec, mu_vec, cs_vec)[source]

Analyze M/G/1 queue with Longest Remaining Processing Time (LRPT).

Matches MATLAB qsys_mg1_lrpt.m exactly: - Exponential case: numerical integration of E[T(x)] * f_k(x) - General case: preemptive priority with descending service time ordering

Parameters:

  • lambda_vec (numpy.ndarray) – Vector of arrival rates per class

  • mu_vec (numpy.ndarray) – Vector of service rates per class

  • cs_vec (numpy.ndarray) – Vector of coefficients of variation per class

Returns:

W: Vector of mean response times per class rho: System utilization (rhohat format)

Return type:

Tuple of (W, rho)

qsys_mg1_psjf(lambda_vec, mu_vec, cs_vec)[source]

Analyze M/G/1 queue with Preemptive Shortest Job First (PSJF).

Matches MATLAB qsys_mg1_psjf.m exactly: - Exponential case: numerical integration of E[T(x)] * f_k(x) - General case: class-based truncated moment formula

Parameters:

  • lambda_vec (numpy.ndarray) – Vector of arrival rates per class

  • mu_vec (numpy.ndarray) – Vector of service rates per class

  • cs_vec (numpy.ndarray) – Vector of coefficients of variation per class

Returns:

W: Vector of mean response times per class rho: System utilization (rhohat format)

Return type:

Tuple of (W, rho)

qsys_mg1_setf(lambda_vec, mu_vec, cs_vec)[source]

Analyze M/G/1 queue with Shortest Expected Time First (SETF).

Matches MATLAB qsys_mg1_setf.m exactly: SETF = FB/LAS + residual service time penalty (non-preemptive).

Parameters:

  • lambda_vec (numpy.ndarray) – Vector of arrival rates per class

  • mu_vec (numpy.ndarray) – Vector of service rates per class

  • cs_vec (numpy.ndarray) – Vector of coefficients of variation per class

Returns:

W: Vector of mean response times per class rho: System utilization (rhohat format)

Return type:

Tuple of (W, rho)

qsys_mm1k_loss(lambda_val, mu, K)[source]

Compute loss probability for M/M/1/K queue.

Uses the closed-form formula for the M/M/1/K loss system where customers are rejected when the buffer (capacity K) is full.

Parameters:

  • lambda_val (float) – Arrival rate

  • mu (float) – Service rate

  • K (int) – Buffer capacity (including customer in service)

Returns:

lossprob: Probability that an arriving customer is rejected rho: Offered load (lambda/mu)

Return type:

Tuple of (lossprob, rho)

References

Original MATLAB: matlab/src/api/qsys/qsys_mm1k_loss.m

qsys_mg1k_loss(lambda_val, service_pdf, K, max_t=100.0)[source]

Compute loss probability for M/G/1/K queue.

Uses the Niu-Cooper transform-free analysis method for computing the stationary distribution and loss probability of M/G/1/K queues.

Parameters:

  • lambda_val (float) – Arrival rate

  • service_pdf (Callable[[float], float]) – Probability density function of service time f(t)

  • K (int) – Buffer capacity

  • max_t (float) – Maximum integration time (default: 100)

Returns:

sigma: Normalizing constant term rho: Offered load lossprob: Probability of loss

Return type:

Tuple of (sigma, rho, lossprob)

References

Original MATLAB: matlab/src/api/qsys/qsys_mg1k_loss.m Niu-Cooper, “Transform-Free Analysis of M/G/1/K”, 1993

qsys_mg1k_loss_mgs(lambda_val, mu, mu_scv, K)[source]

Compute loss probability for M/G/1/K using MacGregor Smith approximation.

Matches MATLAB qsys_mg1k_loss_mgs.m exactly.

Parameters:

  • lambda_val (float) – Arrival rate

  • mu (float) – Service rate

  • mu_scv (float) – Squared coefficient of variation of service time

  • K (int) – Buffer capacity

Returns:

lossprob: Probability of loss rho: Offered load

Return type:

Tuple of (lossprob, rho)

References

Original MATLAB: matlab/src/api/qsys/qsys_mg1k_loss_mgs.m J. MacGregor Smith, “Optimal Design and Performance Modelling of M/G/1/K Queueing Systems”

qsys_mxm1(lambda_batch, mu, E_X_or_batch_sizes, E_X2_or_pmf, mode=None)[source]

Analyze MX/M/1 queue with batch arrivals.

Matches MATLAB qsys_mxm1.m exactly.

Three input formats:

  1. Moment-based: qsys_mxm1(lambda_batch, mu, E_X, E_X2)

  2. PMF-based: qsys_mxm1(lambda_batch, mu, batch_sizes, pmf)

  3. Variance: qsys_mxm1(lambda_batch, mu, E_X, Var_X, ‘variance’)

Parameters:

  • lambda_batch (float) – Batch arrival rate

  • mu (float) – Service rate

  • E_X_or_batch_sizes – Mean batch size (scalar) or array of batch sizes

  • E_X2_or_pmf – Second moment of batch size, PMF, or variance

  • mode (str | None) – Optional ‘variance’ flag for variance-based input

Returns:

W: Mean time in system Wq: Mean waiting time in queue U: Server utilization Q: Mean queue length (including service)

Return type:

Tuple of (W, Wq, U, Q)

References

Original MATLAB: matlab/src/api/qsys/qsys_mxm1.m

class QueueResult(meanQueueLength, meanWaitingTime, meanSojournTime, utilization, queueLengthDist=None, queueLengthMoments=None, sojournTimeMoments=None, analyzer='native')[source]

Bases: object

Result structure for queue analysis.

analyzer: str = 'native'
queueLengthDist: numpy.ndarray | None = None
queueLengthMoments: numpy.ndarray | None = None
sojournTimeMoments: numpy.ndarray | None = None
meanQueueLength: float
meanWaitingTime: float
meanSojournTime: float
utilization: float
ph_to_map(alpha, T)[source]

Convert a PH distribution to its equivalent MAP representation.

For a PH renewal process, the MAP has:

D0 = T (transitions within the PH, no arrival) D1 = t * alpha where t = -T*e (exit rates times restart distribution)

Parameters:
Returns:

D0: MAP hidden transition matrix D1: MAP observable transition matrix

Return type:

Tuple of (D0, D1)

qsys_phph1(alpha, T, beta, S, numQLMoms=3, numQLProbs=100, numSTMoms=3)[source]

Analyze a PH/PH/1 queue using matrix-analytic methods.

Converts the arrival PH to MAP representation and uses the MMAPPH1FCFS solver from BuTools.

Parameters:

  • alpha (numpy.ndarray) – Arrival PH initial probability vector (1 x n)

  • T (numpy.ndarray) – Arrival PH generator matrix (n x n)

  • beta (numpy.ndarray) – Service PH initial probability vector (1 x m)

  • S (numpy.ndarray) – Service PH generator matrix (m x m)

  • numQLMoms (int) – Number of queue length moments to compute (default: 3)

  • numQLProbs (int) – Number of queue length probabilities (default: 100)

  • numSTMoms (int) – Number of sojourn time moments (default: 3)

Returns:

QueueResult with queue performance metrics

Return type:

QueueResult

References

Original MATLAB: matlab/src/api/qsys/qsys_phph1.m

qsys_mapph1(D0, D1, beta, S, numQLMoms=3, numQLProbs=100, numSTMoms=3)[source]

Analyze a MAP/PH/1 queue.

Parameters:

  • D0 (numpy.ndarray) – MAP hidden transition matrix

  • D1 (numpy.ndarray) – MAP observable transition matrix

  • beta (numpy.ndarray) – Service PH initial probability vector

  • S (numpy.ndarray) – Service PH generator matrix

  • numQLMoms (int) – Number of queue length moments

  • numQLProbs (int) – Number of queue length probabilities

  • numSTMoms (int) – Number of sojourn time moments

Returns:

QueueResult with queue performance metrics

Return type:

QueueResult

References

Original MATLAB: matlab/src/api/qsys/qsys_mapph1.m

qsys_mapm1(D0, D1, mu)[source]

Analyze a MAP/M/1 queue.

Parameters:
Returns:

QueueResult with queue performance metrics

Return type:

QueueResult

References

Original MATLAB: matlab/src/api/qsys/qsys_mapm1.m

qsys_mapmc(D0, D1, mu, c)[source]

Analyze a MAP/M/c queue.

Parameters:

  • D0 (numpy.ndarray) – MAP hidden transition matrix

  • D1 (numpy.ndarray) – MAP observable transition matrix

  • mu (float) – Service rate per server

  • c (int) – Number of servers

Returns:

QueueResult with queue performance metrics

Return type:

QueueResult

References

Original MATLAB: matlab/src/api/qsys/qsys_mapmc.m

qsys_mapmap1(D0_arr, D1_arr, D0_srv, D1_srv)[source]

Analyze a MAP/MAP/1 queue.

Both arrival and service processes are Markovian Arrival Processes.

Parameters:

  • D0_arr (numpy.ndarray) – Arrival MAP hidden transition matrix

  • D1_arr (numpy.ndarray) – Arrival MAP observable transition matrix

  • D0_srv (numpy.ndarray) – Service MAP hidden transition matrix

  • D1_srv (numpy.ndarray) – Service MAP observable transition matrix

Returns:

QueueResult with queue performance metrics

Return type:

QueueResult

References

Original MATLAB: matlab/src/api/qsys/qsys_mapmap1.m

qsys_mapg1(D0, D1, service_moments, num_ql_moms=3, num_ql_probs=100, num_st_moms=3)[source]

Analyze a MAP/G/1 queue using BuTools MMAPPH1FCFS.

The general service time distribution is fitted to a Phase-Type (PH) distribution using moment matching before analysis.

Parameters:

  • D0 (numpy.ndarray) – MAP hidden transition matrix (n x n)

  • D1 (numpy.ndarray) – MAP arrival transition matrix (n x n)

  • service_moments (numpy.ndarray) – First k raw moments of service time [E[S], E[S^2], …] (k = 2 or 3 for best accuracy)

  • num_ql_moms (int) – Number of queue length moments to compute (default: 3)

  • num_ql_probs (int) – Number of queue length probabilities (default: 100)

  • num_st_moms (int) – Number of sojourn time moments to compute (default: 3)

Returns:

QueueResult with queue performance metrics

Return type:

QueueResult

References

Original MATLAB: matlab/src/api/qsys/qsys_mapg1.m

Note

Uses the MMAPPH1FCFS solver from BuTools after fitting the general service distribution to a PH distribution.

class QueueType(*values)[source]

Bases: Enum

Queueing system topology types.

STANDARD = 'standard'
RETRIAL = 'retrial'
RENEGING = 'reneging'
RETRIAL_RENEGING = 'retrial_reneging'
class BmapMatrix(D0, D_batch)[source]

Bases: object

Batch Markovian Arrival Process matrix representation.

D₀ = drift matrix (no arrivals) D₁, D₂, …, Dₖ = batch arrival matrices for batch sizes 1, 2, …, K

Properties:

order: Dimension of the BMAP num_batches: Maximum batch size arrival_rate: Overall arrival rate

__post_init__()[source]

Validate BMAP structure.

property arrival_rate: float

Compute overall arrival rate of the BMAP.

lambda = theta * (-D0) * e where theta is stationary distribution of BMAP Markov chain and e is unit vector.

property fundamental_arrival_rate: float

Arrival rate from fundamental matrix.

D0: numpy.ndarray
D_batch: List[numpy.ndarray]
class PhDistribution(beta, S)[source]

Bases: object

Phase-Type (PH) distribution representation.

Beta = initial probability vector (shape: m,) S = transient generator matrix (shape: m x m)

where m is the number of phases.

Properties:

mean: Mean of PH distribution (1/mu in queue notation) scv: Squared coefficient of variation num_phases: Number of phases

__post_init__()[source]

Validate PH distribution structure.

property mean: float

Compute mean of PH distribution.

E[X] = -beta * S^{-1} * e

property mean_squared: float

Compute second moment of PH distribution.

E[X²] = 2 * beta * S^{-2} * e

property scv: float

Squared coefficient of variation of PH distribution.

SCV = (E[X²] / E[X]²) - 1

beta: numpy.ndarray
S: numpy.ndarray
class QbdStatespace(A0, A1, A2, B0, max_level)[source]

Bases: object

Quasi-Birth-Death Markov chain state space representation.

For a QBD process, states are of the form (n, i) where: - n = level (number of retrying customers in orbit) - i = phase (service phase or phase-type stage)

Generator matrix has block structure: Q = | B_0 A_0 0 0 … |

A_2 A_1 A_0 0 … |
0 A_2 A_1 A_0 … |
… |
Properties:

max_level: Maximum retrial orbit size (truncation level) phase_dim: Number of phases at each level total_states: Total state space dimension

__post_init__()[source]

Validate QBD state space.

property binomial_state_dimension: int

Compute state space dimension using binomial coefficient formula.

For a retrial queue with N total customers and m servers, dimension = C(N + m - 1, m - 1)

This is a rough upper bound for QBD truncation.

get_level_phase(idx)[source]

Map linear state index to (level, phase).

Parameters:

idx (int) – Linear state index

Returns:

Tuple (level, phase)

Return type:

Tuple[int, int]

get_state_index(level, phase)[source]

Map (level, phase) to linear state index.

Parameters:

  • level (int) – Retrial orbit size (0 ≤ level ≤ max_level)

  • phase (int) – Phase within level (0 ≤ phase < phase_dim)

Returns:

Linear state index

Return type:

int

A0: numpy.ndarray
A1: numpy.ndarray
A2: numpy.ndarray
B0: numpy.ndarray
max_level: int
class RetrialQueueResult(queue_type, L_orbit, N_server, utilization, throughput, P_idle, P_empty_orbit, stationary_dist=None, truncation_level=0, converged=False, iterations=0, error=numpy.inf)[source]

Bases: object

Analysis results for BMAP/PH/N/N retrial queue.

queue_type

Type of queue (retrial, reneging, etc.)

Type:

line_solver.api.qsys.retrial.QueueType

L_orbit

Expected number of customers in orbit

Type:

float

N_server

Expected number of customers being served

Type:

float

utilization

Server utilization

Type:

float

throughput

System throughput

Type:

float

P_idle

Probability that server is idle

Type:

float

P_empty_orbit

Probability that orbit is empty

Type:

float

stationary_dist

Stationary distribution vector

Type:

numpy.ndarray | None

truncation_level

QBD truncation level used

Type:

int

converged

Whether numerical solution converged

Type:

bool

iterations

Number of iterations to convergence

Type:

int

error

Final error estimate

Type:

float

converged: bool = False
iterations: int = 0
stationary_dist: numpy.ndarray | None = None
truncation_level: int = 0
queue_type: QueueType
L_orbit: float
N_server: float
utilization: float
throughput: float
P_idle: float
P_empty_orbit: float
class RetrialQueueAnalyzer(sn, options=None)[source]

Bases: object

Framework for analyzing BMAP/PH/N/N retrial queues.

This class provides the foundation for future full solver implementation, including topology detection, parameter extraction, and QBD setup.

Example

analyzer = RetrialQueueAnalyzer(model) queue_type = analyzer.detect_queue_type() if queue_type == QueueType.RETRIAL:

result = analyzer.analyze()

Initialize retrial queue analyzer.

Parameters:

  • sn (Any) – NetworkStruct with queue configuration

  • options (Dict | None) – Analysis options (tolerance, max iterations, etc.)

__init__(sn, options=None)[source]

Initialize retrial queue analyzer.

Parameters:

  • sn (Any) – NetworkStruct with queue configuration

  • options (Dict | None) – Analysis options (tolerance, max iterations, etc.)

analyze()[source]

Analyze the retrial queue.

This is the main entry point for analysis: 1. Detect queue type 2. Extract arrival and service parameters 3. Build QBD state space 4. Solve for stationary distribution using matrix-analytic methods 5. Compute performance metrics

Returns:

RetrialQueueResult with performance metrics

Return type:

RetrialQueueResult

build_qbd_statespace(bmap, ph_service, retrial_params)[source]

Build QBD state space for the retrial queue.

This constructs the generator matrix blocks for the QBD process following the BMAP/PH/N/N retrial queue formulation.

Parameters:

  • bmap (BmapMatrix) – BMAP arrival process

  • ph_service (PhDistribution) – PH service distribution

  • retrial_params (Dict[str, float]) – Retrial parameters (alpha, gamma, p, R, N)

Returns:

QbdStatespace instance or None if construction fails

Return type:

QbdStatespace | None

detect_queue_type()[source]

Detect the type of queueing system from topology.

Returns:

QueueType enum indicating retrial, reneging, or combination

Return type:

QueueType

extract_bmap()[source]

Extract BMAP parameters from arrival process.

LINE stores arrival processes in MAP format: {D0, D1, D2, …} where D0 is the “hidden” generator and D1, D2, … are arrival matrices.

Returns:

BmapMatrix instance or None if arrival is not BMAP/MAP

Return type:

BmapMatrix | None

extract_ph_service()[source]

Extract Phase-Type service parameters.

LINE stores PH in MAP format: {D0, D1} D0 = T (subgenerator matrix) D1 = S0 * alpha (exit rate times initial prob)

Returns:

PhDistribution instance or None if service is not PH

Return type:

PhDistribution | None

extract_retrial_parameters()[source]

Extract retrial-specific parameters.

Returns:

  • alpha: Retrial rate (rate at which customers retry)

  • gamma: Orbit impatience rate (reneging rate)

  • p: Batch rejection probability

  • R: Threshold for admission control

  • N: Number of servers

Return type:

Dict with keys

qsys_bmapphnn_retrial(arrival_matrix, service_params, N, retrial_params=None, options=None)[source]

Analyze BMAP/PH/N/N bufferless retrial queue.

Implements the algorithm from Dudin et al., “Analysis of BMAP/PH/N-Type Queueing System with Flexible Retrials Admission Control”, Mathematics 2025, 13(9), 1434.

Parameters:

  • arrival_matrix (Dict[str, numpy.ndarray]) – Dict with ‘D0’, ‘D1’, … for BMAP matrices. D0: hidden transition matrix (V x V). D1, …, DK: arrival matrices for batch sizes 1, …, K.

  • service_params (Dict[str, numpy.ndarray]) – Dict with ‘beta’ (initial prob vector, 1xM) and ‘S’ (PH subgenerator matrix, MxM).

  • N (int) – Number of servers (also capacity, hence bufferless).

  • retrial_params (Dict[str, float] | None) – Dict with ‘alpha’ (retrial rate per customer), ‘gamma’ (impatience/abandonment rate), ‘p’ (batch rejection probability), ‘R’ (admission threshold, scalar or 1xV).

  • options (Dict | None) – Dict with optional keys: ‘MaxLevel’: max orbit level for truncation (default: auto). ‘Tolerance’: convergence tolerance (default: 1e-10). ‘Verbose’: print progress (default: False).

Returns:

RetrialQueueResult with performance metrics.

Return type:

RetrialQueueResult

References

Dudin, A., Klimenok, V., & Vishnevsky, V. (2020). Port from: matlab/src/api/qsys/qsys_bmapphnn_retrial.m

qsys_is_retrial(sn)[source]

Check if network is a valid BMAP/PH/N/N bufferless retrial queue.

Validates that the network structure matches the requirements for the BMAP/PH/N/N retrial queue solver: - Single bufferless queue (capacity == number of servers) - Retrial drop strategy configured - BMAP/MAP arrival process at source - PH/Exp service at queue - Open class model

Based on: Dudin et al., “Analysis of BMAP/PH/N-Type Queueing System with Flexible Retrials Admission Control”, Mathematics 2025, 13(9), 1434.

Parameters:

sn (Any) – NetworkStruct object

Returns:

is_retrial: True if network is valid BMAP/PH/N/N retrial topology retrial_info: RetrialInfo with parameters for the retrial solver

Return type:

Tuple of (is_retrial, retrial_info) where

References

Original MATLAB: matlab/src/api/qsys/qsys_is_retrial.m

class RetrialInfo(is_retrial, station_idx=None, node_idx=None, source_idx=None, class_idx=None, error_msg='', N=None, alpha=0.1, gamma=0.0, p=0.0, R=None)[source]

Bases: object

Information about a valid retrial queue topology.

N: int | None = None
R: int | None = None
alpha: float = 0.1
class_idx: int | None = None
error_msg: str = ''
gamma: float = 0.0
node_idx: int | None = None
p: float = 0.0
source_idx: int | None = None
station_idx: int | None = None
is_retrial: bool
class RenegingInfo(is_reneging, source_idx=None, queue_idx=None, class_idx=None, n_servers=None, service_rate=None, error_msg='')[source]

Bases: object

Information about a valid reneging queue topology.

class_idx: int | None = None
error_msg: str = ''
n_servers: int | None = None
queue_idx: int | None = None
service_rate: float | None = None
source_idx: int | None = None
is_reneging: bool
detect_reneging_topology(sn)[source]

Detect if model is suitable for MAP/M/s+G (MAPMsG) reneging solver.

Requirements: - Open model, single class - Single queue station with reneging/patience configured - MAP/BMAP arrival at source - Exponential service at queue (single-phase PH) - FCFS scheduling

Parameters:

sn (Any) – NetworkStruct object

Returns:

Tuple of (is_reneging, reneging_info)

Return type:

Tuple[bool, RenegingInfo]

References

Original MATLAB: matlab/src/solvers/MAM/solver_mam_retrial.m (detectRenegingTopology)

has_reneging_patience(sn)[source]

Check if model has reneging/patience configured on any queue station.

Returns True if any queue station has ImpatienceType.RENEGING configured with a patience distribution.

Parameters:

sn (Any) – NetworkStruct object

Returns:

True if reneging patience is configured

Return type:

bool

References

Original MATLAB: matlab/src/solvers/MAM/solver_mam_analyzer.m (hasRenegingPatience)

extract_bmap_matrices(proc)[source]

Extract BMAP matrices {D0, D1, …} from LINE process representation.

LINE stores arrival processes in MAP format: {D0, D1, D2, …} where D0 is the “hidden” generator and D1, D2, … are arrival matrices. May also be in PH format {alpha, T} which is converted to MAP.

Parameters:

proc (Any) – Process representation from sn.proc[station][class]

Returns:

List of numpy arrays [D0, D1, …] or None if extraction fails

Return type:

List[numpy.ndarray] | None

References

Original MATLAB: matlab/src/solvers/MAM/solver_mam_retrial.m (extractBMAPMatrices)

extract_ph_params(proc)[source]

Extract PH parameters (beta, S) from LINE process representation.

LINE stores PH in MAP format: {D0, D1} where D0=T (subgenerator), D1=S0*alpha (exit rate times initial prob).

Parameters:

proc (Any) – Process representation from sn.proc[station][class]

Returns:

Tuple of (beta, S) where beta is initial probability vector and S is subgenerator matrix. Returns (None, None) if extraction fails.

Return type:

Tuple[numpy.ndarray | None, numpy.ndarray | None]

References

Original MATLAB: matlab/src/solvers/MAM/solver_mam_retrial.m (extractPHParams)

convert_patience_to_regimes(patience_proc, options=None)[source]

Convert patience distribution to piecewise-constant abandonment regimes for MAPMsG.

Converts a patience distribution (in MAP/PH format) to boundary levels and abandonment function values for the MRMFQ solver.

Parameters:

  • patience_proc (Any) – Patience distribution from sn.patienceProc[station, class]

  • options (Dict | None) – Dict with optional ‘mapmsg_quantization’ key (default 11)

Returns:

boundary_levels: Array of regime boundary time points ga: Array of abandonment probabilities at each regime quantization: Number of regimes

Return type:

Tuple of (boundary_levels, ga, quantization) where

References

Original MATLAB: matlab/src/solvers/MAM/solver_mam_retrial.m (convertPatienceToRegimes)

solver_mam_retrial(sn, options=None)[source]

Solve queueing models with customer impatience (retrial or reneging).

Dispatches to either: 1. RETRIAL: BMAP/PH/N/N bufferless retrial solver 2. RENEGING: MAP/M/s+G solver (MAPMsG)

Parameters:

  • sn (Any) – NetworkStruct object

  • options (Dict | None) – Solver options dict with optional keys: ‘iter_max’: Maximum truncation level (default 150) ‘tol’: Convergence tolerance (default 1e-10) ‘verbose’: Print progress messages (default False) ‘config’: Dict with ‘mapmsg_quantization’ (default 11)

Returns:

QN: (M, K) queue lengths UN: (M, K) server utilizations RN: (M, K) response times TN: (M, K) throughputs CN: (1, K) cycle times XN: (1, K) system throughputs totiter: iteration/truncation level count

Return type:

Tuple of (QN, UN, RN, TN, CN, XN, totiter) where

References

Original MATLAB: matlab/src/solvers/MAM/solver_mam_retrial.m