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:

• Exact solutions: qsys_mm1(), qsys_mmk(), qsys_mg1(), qsys_gm1()

• G/G/1 approximations: qsys_gig1_approx_kingman(), qsys_gig1_approx_kimura(), qsys_gig1_approx_marchal()

• Multi-server approximations: qsys_gigk_approx(), qsys_gigk_approx_kingman()

• Loss systems: qsys_mm1k_loss(), qsys_mg1k_loss(), qsys_mg1k_loss_mgs()

• Performance bounds: qsys_gig1_ubnd_kingman()

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)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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')

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)

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 (numpy.ndarray) – Initial probability vector (1 x n)

• T (numpy.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)

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)

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)

Analyze a MAP/M/1 queue.

Parameters:

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

• D1 (numpy.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)

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)

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)

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)

Bases: Enum

Queueing system topology types.

STANDARD = 'standard'

RETRIAL = 'retrial'

RENEGING = 'reneging'

RETRIAL_RENEGING = 'retrial_reneging'

class BmapMatrix(D0, D_batch)

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__()

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)

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__()

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)

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__()

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)

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)

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)

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)

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)

Initialize retrial queue analyzer.

Parameters:

• sn (Any) – NetworkStruct with queue configuration

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

analyze()

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)

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

Detect the type of queueing system from topology.

Returns:

QueueType enum indicating retrial, reneging, or combination

Return type:

QueueType

extract_bmap()

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

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

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)

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)

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)

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