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 (ndarray) – MAP hidden transition matrix (n x n).

  • D1 (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 (ndarray) – MAP hidden transition matrix.

  • D1 (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_gig1_approx_allencunneen(lambda_val, mu, ca, cs)[source]

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

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:

Performance measures including L, Lq, W, Wq, rho

Return type:

dict

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

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

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:

Performance measures including L, Lq, W, Wq, rho

Return type:

dict

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

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

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:

Performance measures including L, Lq, W, Wq, rho

Return type:

dict

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

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

Uses QNA (Queueing Network Analyzer) approximation.

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:

Performance measures including L, Lq, W, Wq, rho

Return type:

dict

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

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

A simple two-moment approximation.

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:

Performance measures including W, Wq, rho, rhohat

Return type:

dict

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

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

Uses an exponential decay formulation for the modified utilization.

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:

Performance measures including W, Wq, rho, rhohat

Return type:

dict

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

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

A simple two-moment approximation.

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:

Performance measures including W, Wq, rho

Return type:

dict

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

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

Modified utilization-based approximation.

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:

Performance measures including W, Wq, rho

Return type:

dict

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

Approximation for G/G/k queue.

Uses Allen-Cunneen approximation extended to multi-server.

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:

Performance measures including L, Lq, W, Wq, rho

Return type:

dict

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.

Computes per-class mean waiting times for multiple priority classes using the Pollaczek-Khinchine formula extended for priorities.

For K priority classes (class 1 = highest priority), the waiting time for class k is:

W_k = (B_0 / (1 - sum_{i=1}^{k-1} rho_i)) * (1 / (1 - sum_{i=1}^{k} rho_i))

Parameters:

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

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

  • cs_vec (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)

Raises:

  • ValueError – If inputs have different lengths or are not positive

  • ValueError – If system is unstable (utilization >= 1)

References

Original MATLAB: matlab/src/api/qsys/qsys_mg1_prio.m Kleinrock, “Queueing Systems, Volume I: Theory”, Section 3.5

qsys_mg1_srpt(lambda_vec, mu_vec, cs_vec)[source]

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

For SRPT, smaller jobs have preemptive priority over larger jobs. When all service time distributions are exponential (cs=1), SRPT reduces to preemptive priority based on service rates.

Parameters:

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

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

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

Returns:

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

Return type:

Tuple of (W, rho)

References

Original MATLAB: matlab/src/api/qsys/qsys_mg1_srpt.m Schrage and Miller, “The queue M/G/1 with the shortest remaining processing time discipline”, Operations Research, 1966

qsys_mg1_fb(lambda_val, mu_val, cs_val)[source]

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

In FB scheduling, the job with the least attained service gets priority. This is also known as Least Attained Service (LAS) scheduling.

For FB, the conditional response time for a job of size x is:

E[T|x] = x / (1 - rho) for exponential service E[T|x] = integral formula for general service

Parameters:

  • lambda_val (float) – Arrival rate

  • mu_val (float) – Service rate

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

Returns:

W: Mean response time rho: System utilization (rhohat format)

Return type:

Tuple of (W, rho)

References

Original MATLAB: matlab/src/api/qsys/qsys_mg1_fb.m Kleinrock, “Queueing Systems, Volume II”

qsys_mg1_lrpt(lambda_vec, mu_vec, cs_vec)[source]

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

LRPT gives preemptive priority to jobs with the longest remaining processing time. This is the opposite of SRPT.

Parameters:

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

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

  • cs_vec (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)

References

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

qsys_mg1_psjf(lambda_vec, mu_vec, cs_vec)[source]

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

PSJF gives preemptive priority to the job with the smallest expected total service time (not remaining time).

For exponential service, PSJF equals SRPT.

Parameters:

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

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

  • cs_vec (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)

References

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

qsys_mg1_setf(lambda_vec, mu_vec, cs_vec)[source]

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

SETF is non-preemptive and prioritizes based on expected service time.

Parameters:

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

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

  • cs_vec (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)

References

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

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, cs, K)[source]

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

Uses the MacGregor Smith approximation which provides good accuracy for moderate coefficient of variation.

Parameters:

  • lambda_val (float) – Arrival rate

  • mu (float) – Service rate

  • cs (float) – 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 MacGregor Smith, “Queueing with capacity and performance”

qsys_mxm1(lambda_vec, mu, batch_sizes=None)[source]

Analyze M^X/M/1 queue (batch arrivals).

Computes performance metrics for a queue with batch Poisson arrivals and exponential service.

Parameters:

  • lambda_vec (ndarray) – Vector of arrival rates for each batch size (lambda_vec[k] = rate of batches of size k+1)

  • mu (float) – Service rate

  • batch_sizes (ndarray | None) – Optional array of batch sizes (default: 1, 2, 3, …)

Returns:

L: Mean number in system W: Mean response time rho: Utilization

Return type:

Tuple of (L, W, rho)

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: ndarray | None = None
queueLengthMoments: ndarray | None = None
sojournTimeMoments: 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:

  • alpha (ndarray) – Initial probability vector (1 x n)

  • T (ndarray) – Sub-generator matrix (n x n)

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 (ndarray) – Arrival PH initial probability vector (1 x n)

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

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

  • S (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 (ndarray) – MAP hidden transition matrix

  • D1 (ndarray) – MAP observable transition matrix

  • beta (ndarray) – Service PH initial probability vector

  • S (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:

  • D0 (ndarray) – MAP hidden transition matrix

  • D1 (ndarray) – MAP observable transition matrix

  • mu (float) – Service rate

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 (ndarray) – MAP hidden transition matrix

  • D1 (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 (ndarray) – Arrival MAP hidden transition matrix

  • D1_arr (ndarray) – Arrival MAP observable transition matrix

  • D0_srv (ndarray) – Service MAP hidden transition matrix

  • D1_srv (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 (ndarray) – MAP hidden transition matrix (n x n)

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

  • service_moments (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: ndarray
D_batch: List[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: ndarray
S: 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: ndarray
A1: ndarray
A2: ndarray
B0: 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=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
error: float = inf
iterations: int = 0
stationary_dist: 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. Currently returns a placeholder result. Full implementation would: 1. Detect queue type 2. Extract arrival and service parameters 3. Build QBD state space 4. Solve for stationary distribution 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.

Parameters:
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.

Returns:

BmapMatrix instance or None if arrival is not BMAP/MAP

Return type:

BmapMatrix | None

extract_ph_service()[source]

Extract Phase-Type service parameters.

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

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.

Parameters:

  • arrival_matrix (Dict[str, ndarray]) – Dict with ‘D0’, ‘D1’, … for BMAP

  • service_params (Dict[str, float]) – Dict with ‘beta’ (initial prob) and ‘S’ (gen matrix) for PH

  • N (int) – Number of customers (for closed queue)

  • retrial_params (Dict[str, float] | None) – Dict with ‘alpha’ (retrial rate), ‘gamma’ (impatience), ‘p’ (rejection prob), ‘R’ (threshold)

  • options (Dict | None) – Solver options

Returns:

RetrialQueueResult with performance metrics

Raises:

NotImplementedError – Full solver not yet implemented

Return type:

RetrialQueueResult

References

Dudin, A., Klimenok, V., & Vishnevsky, V. (2020). “The Theory of Quasi-Birth-Death Processes and Its Applications”. Springer.

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