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()
Discrete time (slotted): qsys_geogeo1(), qsys_geoxgeo1(), qsys_geoxgeo1_moments()

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 Discrete time (slotted): qsys_geogeo1, qsys_geoxgeo1

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_mdc_crommelin(lambda_arr, s, c, truncation=-1)[source]

Solve M/D/c via Crommelin’s embedded DTMC.

Parameters:

lambda_arr (float) – Poisson arrival rate.
s (float) – Deterministic service time (>0).
c (int) – Number of servers (>=1).
truncation (int) – Optional state-space truncation level. If <=0, chosen automatically to keep dense LU manageable.

Returns:

mean_queue_length: E[N] (number in system) mean_waiting_queue: Lq = E[(N-c)+] mean_waiting_time: Wq = Lq / lambda mean_sojourn_time: W = Wq + s utilization: rho = lambda * s / c

Return type:

dict with keys

qsys_dmc(lambda_arr, mu, c, truncation=-1, quad_steps=200)[source]

Solve D/M/c via embedded DTMC at arrival epochs.

Parameters:

lambda_arr (float) – Arrival rate (deterministic interarrivals of mean 1/lambda).
mu (float) – Service rate per server.
c (int) – Number of servers.
truncation (int) – Optional state-space truncation. Auto-chosen if <=0.
quad_steps (int) – Number of trapezoidal-rule steps for cycle integration.

Returns:

dict with mean_queue_length, mean_waiting_queue, mean_waiting_time, mean_sojourn_time, utilization.

Return type:

Dict

qsys_phm1(alpha, T, mu)[source]

Solve PH/M/1 via the GI/M/1 sigma-root.

Parameters:

alpha – PH entry probability vector (length k).
T – PH sub-generator matrix (k x k), row sums <= 0.
mu (float) – Exponential service rate (>0).

Returns:

dict with mean_queue_length, mean_waiting_queue, mean_waiting_time, mean_sojourn_time, utilization, sigma.

Return type:

Dict

qsys_phmc(alpha, T, mu, c, max_iter=50000, tol=1e-14)[source]

Solve PH/M/c via matrix-geometric.

Parameters:

alpha – PH entry probability vector (length k).
T – PH sub-generator matrix (k x k), row sums <= 0.
mu (float) – Exponential service rate per server (>0).
c (int) – Number of servers (>=1).
max_iter (int) – Maximum iterations for R fixed-point.
tol (float) – Convergence tolerance for R.

Returns:

dict with mean_queue_length, mean_waiting_queue, mean_waiting_time, mean_sojourn_time, utilization.

Return type:

Dict

qsys_geogeo1(a, s, convention=LAS_DA)[source]

Analyze a discrete-time Geo/Geo/1 queue.

In each slot an arrival occurs with probability a and, if the server is engaged, a service completion occurs with probability s, independently of everything else. The system content at slot boundaries is a discrete birth-death chain with a geometric stationary distribution.

The two conventions are not two systems. They are one system, Daduna’s LA-rule (events at the end of their slot) with the D/A-rule (departure resolved before arrival), observed at two instants. With X(t+1) = X(t) - D(t) + A(t), LAS_DA is the law of X, taken after both events, and EAS is the law of Y(t) = X(t) - D(t), taken after the departure and before the arrival. They are one departure apart, so the mean contents differ by exactly a and the sojourn times by one slot; the queueing delay is the same under both.

Parameters:

a (float) – Per-slot arrival probability, 0 < a < s
s (float) – Per-slot service completion probability, 0 < s <= 1
convention (str) – Observation epoch, ‘LAS_DA’ (default) or ‘EAS’

Returns:

convention: epoch used
arrivalProb, serviceProb: the inputs
utilization: a/s
throughput: a
emptyProb: probability the system is empty at the epoch
ratio: geometric decay ratio r = a(1-s)/(s(1-a))
meanQueueLength, meanWaitingQueue
meanSojournTime, meanWaitingTime, meanServiceTime
pmf: callable n -> stationary probability of n jobs
analyzer: ‘qsys_geogeo1’

Return type:

dict with keys

Example

>>> r = qsys_geogeo1(0.2, 0.5)
>>> print(f"{r['meanSojournTime']:.4f}")
2.6667

qsys_geoxgeo1(a, beta, s, convention=LAS_DA)[source]

Analyze a discrete-time Geo^X/Geo/1 queue with geometric batch sizes.

In each slot a batch arrives with probability a; the batch size is geometric on {1,2,…} with parameter beta, so E[X] = 1/beta and E[X(X-1)] = 2(1-beta)/beta**2. Stability requires lambda = a*E[X] < s. At beta == 1 the batch is always a single job and the result equals qsys_geogeo1().

Parameters:

a (float) – Per-slot probability that a batch arrives, 0 < a <= 1
beta (float) – Batch-size geometric parameter, 0 < beta <= 1
s (float) – Per-slot service completion probability, 0 < s <= 1
convention (str) – Observation epoch, ‘LAS_DA’ (default) or ‘EAS’

Returns:

dict, see qsys_geoxgeo1_moments()

Return type:

Dict[str, object]

Example

>>> r = qsys_geoxgeo1(0.1, 0.5, 0.9)
>>> print(f"{r['arrivalRate']:.4f}")
0.2000

qsys_geoxgeo1_moments(a, batch_mean, batch_second_factorial, s, convention=LAS_DA)[source]

Analyze a discrete-time Geo^X/Geo/1 queue for an arbitrary batch law.

The batch enters the solution only through its first two factorial moments, so specifying those is fully general. With A(z) = 1-a+a*X(z) the pgf of the number of jobs arriving in one slot, the slot-boundary content obeys X(t+1) = X(t) - D(t) + A(t) with the departure resolved first, giving:

P(z) = p0 s (z-1) A(z) / ( z - A(z)(s+(1-s)z) ),  p0 = 1 - lambda/s
E[N] = lambda + ( a E[X(X-1)]/2 + lambda(1-s) ) / (s - lambda)

Parameters:

a (float) – Per-slot probability that a batch arrives, 0 < a <= 1
batch_mean (float) – E[X], at least 1 since an arriving batch carries a job
batch_second_factorial (float) – E[X(X-1)], non-negative and at least batch_mean**2 - batch_mean
s (float) – Per-slot service completion probability, 0 < s <= 1
convention (str) – Observation epoch, ‘LAS_DA’ (default) or ‘EAS’

Returns:

convention, batchArrivalProb, batchMean, batchSecondFactorialMoment, serviceProb
arrivalRate, throughput: lambda = a*E[X]
utilization: lambda/s
boundaryEmptyProb: 1 - lambda/s, the empty probability AT THE SLOT BOUNDARY under both conventions
meanQueueLength, meanWaitingQueue
meanSojournTime, meanWaitingTime, meanServiceTime
pgf: callable (z, A_of_z) -> P(z), defined for 0 < z <= 1
analyzer: ‘qsys_geoxgeo1’

Return type:

dict with keys

No pmf is returned: for a general batch law the stationary distribution has no elementary closed form, so only the generating function is exact.

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_mmck(lambda_val, mu, c, K)[source]

Exact closed-form analysis of an M/M/c/K queue (finite capacity K, c servers).

Port of MATLAB qsys_mmck.m. Stationary distribution (truncated Erlang form):: a = lambda/mu, rho = a/c p_n = a^n/n! * p0 for 0 <= n <= c p_n = a^c/c! * rho^(n-c) * p0 for c <= n <= K

with p0 normalizing the (K+1)-point distribution.

Parameters:

lambda_val (float) – Poisson arrival rate (> 0)
mu (float) – Per-server exponential service rate (> 0)
c (int) – Number of servers (>= 1)
K (int) – System capacity, total jobs allowed (K >= c)

Returns:

dict with L, Lq, W, Wq, rho, P0 plus MATLAB-style aliases (meanQueueLength, meanQueueLengthQ, meanWaitingTime, meanSojournTime, utilization, throughput, lossProbability, queueLengthDist).

Return type:

Dict[str, float]

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(sigma, mu)[source]

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

Matches MATLAB qsys_gm1(sigma, mu) and JAR Qsys_gm1: the number of customers found by an arrival is geometric with parameter sigma, so the mean response time (time in system) is W = 1/(mu*(1-sigma)).

Parameters:

sigma (float) – Root in (0,1) of sigma = A*(mu*(1-sigma)), where A* is the Laplace-Stieltjes transform of the interarrival-time distribution.
mu (float) – Service rate

Returns:

{‘W’: mean response time}

Return type:

dict

Note

To obtain sigma from the first two moments of the interarrival time, use qsys_gg1(lambda_val, mu, ca2, 1.0), which fits a two-moment renewal process and solves the fixed point.

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_rq(rho, mu, cs2, IaFun)[source]

Robust Queueing (RQ) approximation for a single G/GI/1 queue partially characterized by its arrival rate, index of dispersion for counts (IDC) and the first two moments of the service time. Implements the mean steady-state workload

Z* = sup_{x>=0} { -(1-rho) x + sqrt( 2 rho x (I_a(x) + c2_s) / mu ) }

and the derived steady-state performance measures.

Reference:: W. Whitt and W. You (2018), “A Robust Queueing Network Analyzer Based on Indices of Dispersion”, eqs. (13),(16)-(18).

Parameters:

rho (float) – Traffic intensity lambda/mu (0<rho<1)
mu (float) – Service rate
cs2 (float) – Service SCV c2_s
IaFun – callable, IaFun(x) -> arrival IDC I_a(x) at time argument x>0

Returns:

Tuple (Z, W, Q, X) with mean workload E[Z], waiting time E[W], queue length E[Q] (waiting + in service), and number in system E[X].

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: alias of qsys_gig1_ubnd_kingman (‘gig1.kingman’ method).

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 (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 diffusion approximation for G/G/1 with instantaneous-return boundary:

p(0) = 1-rho, p(n) = rho*(1-rhat)*rhat^(n-1), n>=1 rhat = exp(-2*(1-rho)/(rho*ca^2+cs^2))

hence E[N] = rho/(1-rhat) and the mean response time (time in system) is W = E[N]/lambda = 1/(mu*(1-rhat)).

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

References

Gelenbe, E. (1975). “On approximate computer system models”. Journal of the ACM 22(2), 261-269.

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

Kimura’s diffusion-interpolation approximation for G/G/1:: Wq = rho*(ca^2+cs^2)/(mu*(1-rho)*(1+ca^2))

exact for M/M/1 and M/G/1. The returned W adds the mean service time (response time, time in system).

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

References

Kimura, T. (1986). “A two-moment approximation for the mean waiting time in the GI/G/s queue”. Management Science 32(6), 751-763.

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

Myskja’s third-moment approximation for G/G/1:: Wq = rho/(2*mu*(1-rho))*((1+cs^2)+(q0/qa)^(1/rho-rho)*(1/rho)*(ca^2-1))

exact for M/G/1 (ca=1). The returned W adds the mean service time (response time, time in system).

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:

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_myskja2(lambda_val, mu, ca, cs, q0, qa)[source]

Modified Myskja (Myskja2) approximation for G/G/1, returning the mean response time (time in system). For ca=1 the interpolation parameter theta is a 0/0 form, so the exact M/G/1 result is returned instead (also the interpolation anchor of the method).

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:

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_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. In the G/M/1 case, the interarrival-time distribution is fitted from (lambda, ca2) by a two-moment renewal process (H2 with balanced means for ca2>1, mixed Erlang for ca2<1) and sigma is the root of sigma = A*(mu*(1-sigma)), with A* the interarrival-time LST.

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 MATLAB: matlab/src/api/qsys/qsys_gg1.m

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]

GI/G/k approximation by interpolation of the M/M/k, M/D/k and D/M/k queues (Cosmetatos 1982; Page 1982):

Wq = [ca^2*cs^2 + ca^2*(1-cs^2)*phi1/2

(1-ca^2)*cs^2*phi3/2] * Wq(M/M/k)

where phi1 and phi3 are the Cosmetatos (1975) correction factors for M/D/k and D/M/k, with the safeguards of Whitt (1993). The D/D/k corner has Wq=0. The interpolation requires ca^2<=1 and cs^2<=1; outside this region the Lee-Longton scaling Wq = ((ca^2+cs^2)/2)*Wq(M/M/k) is used.

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 (time in system) rhohat: Effective utilization

Return type:

Tuple of (W, rhohat) where

References

Cosmetatos, G.P. (1975). “Approximate explicit formulae for the average queueing time in the processes (M/D/r) and (D/M/r)”. INFOR 13, 328-331. Page, E. (1982). “Tables of waiting times for M/M/n, M/D/n and D/M/n and their use to give approximate waiting times in more general queues”. J. Opl. Res. Soc. 33, 453-473.

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

GI/G/k approximation of Whitt (1993), eqs. (2.16)-(2.25):: Wq = phi(rho,ca^2,cs^2,k) * ((ca^2+cs^2)/2) * Wq(M/M/k)

where phi interpolates the Cosmetatos M/D/k (phi1) and D/M/k (phi3) correction factors. Exact for M/M/k; reduces to the Cosmetatos M/D/k approximation for cs=0. Implements eq. (2.25) as printed, which was validated against the paper’s Tables 5-7 (New column).

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 (time in system) rhohat: Effective utilization

Return type:

Tuple of (W, rhohat) where

References

Whitt, W. (1993). “Approximations for the GI/G/m queue”. Production and Operations Management 2(2), 114-161.

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

Kingman’s upper bound on the mean waiting time of a G/G/1 queue:: Wq <= lambda*(sa^2+ss^2)/(2*(1-rho)),

with sa^2=ca^2/lambda^2 and ss^2=cs^2/mu^2. The returned W adds the mean service time, so it upper-bounds the mean response time (time in system).

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

Kingman, J.F.C. (1962). “Some inequalities for the queue GI/G/1”. Biometrika 49(3/4), 315-324. 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_mm1_dps(lambda_vec, mu_vec, w_vec, tol=1e-10, max_cutoff=2048)[source]

Numerically exact M/M/1 Discriminatory Processor Sharing (DPS) queue.

Solves the multiclass DPS continuous-time Markov chain on the per-class population vector (n_1..n_K): arrivals lambda_k, class-k service completion rate mu_k * n_k * w_k / sum_j n_j * w_j. The state space is truncated at a total population level chosen from the geometric tail bound (the total-count process is stochastically dominated by the M/M/1 with rate min_k mu_k), and the truncation level is doubled until the mean queue lengths are stable to the requested tolerance – so the result is exact to solver precision and conserves the M/M/1 total for equal service rates by construction.

Parameters:

lambda_vec (numpy.ndarray) – Per-class Poisson arrival rates (K,)
mu_vec (numpy.ndarray) – Per-class exponential service rates (K,)
w_vec (numpy.ndarray) – Per-class DPS weights (K,), positive
tol (float) – Convergence tolerance on the per-class mean counts
max_cutoff (int) – Hard bound on the total-population truncation level

Returns:

per-class mean response times (K,) via Little’s law, and the total utilization sum_k lambda_k/mu_k.

Return type:

Tuple of (T, rho)

qsys_mg1_srpt(lambda_vec, mu_vec, cs_vec)[source]

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

SRPT is a size-based policy: it always serves the job with the smallest remaining processing time, preempting whenever a shorter job arrives. The class-conditional mean response time follows the Schrage-Miller formula (Bansal-Harchol-Balter, SIGMETRICS 2001, Sec. 4, Eqs (1)-(3), after Schrage-Miller 1966). For a job of size x:

E[T(x)] = E[W(x)] + E[R(x)] E[W(x)] = lambda*(m2(x) + x^2*(1-F(x))) / (2*(1-rho(x))^2) E[R(x)] = integral_0^x dt/(1-rho(t))

with f the mixture job-size density, F its CDF, rho(x)=lambda*int_0^x t f(t)dt and m2(x)=int_0^x t^2 f(t)dt. The per-class mean is E[T_r]=int_0^inf E[T(x)] f_r(x) dx; since E[T(x)] depends only on the job size (SRPT is size-based) this is exact. Integrals use cumulative trapezoidal quadrature on a common grid. Each class is matched to its (mean=1/mu, scv=cs^2): exponential for cs=1, a two-phase balanced hyperexponential for cs>1, and a Tijms Erlang-(k-1)/Erlang-k mixture for cs<1. The fully exponential case reproduces the exact M/M/1/SRPT result.

Matches MATLAB qsys_mg1_srpt.m.

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 (original class order) rho: System load measure Q/(1+Q) with Q = sum(lambda.*W)

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=None)[source]

Exact M/G/1/K loss probability via the Markov chain embedded at service-start epochs (transform-free analysis in the spirit of Niu-Cooper).

State: number of customers waiting in the queue immediately after a service start, q in {0,…,K-2} (capacity K includes the job in service; just after a departure at most K-1 jobs remain, one of which enters service). With a_j = P(j Poisson arrivals during a service time):

q=0if no arrival occurs during the service the system empties and
the next service starts with the next arrival (q’=0), so both a_0 and a_1 lead to q’=0 and j>=2 arrivals lead to q’=j-1;

q>=1: q’ = q-1+j, with arrivals beyond the free capacity lost
(aggregated in the last column).

The loss probability follows from the renewal-reward argument: E[cycle] = E[S] + sigma_0*a_0/lambda, lambda_eff = 1/E[cycle], P_loss = 1 - lambda_eff/lambda = 1 - 1/(rho + sigma_0*a_0)

where sigma is the stationary distribution at service-start epochs.

Parameters:

lambda_val (float) – Arrival rate
service_pdf (Callable[[float], float]) – Probability density function of service time f(t)
K (int) – Buffer capacity (including the customer in service)
max_t (float | None) – Maximum integration time (default: smallest horizon covering the service-time distribution mass to within 1e-10)

Returns:

sigma0: Stationary probability of an empty queue at service-start: epochs

rho: Offered load lossprob: Probability of loss

Return type:

Tuple of (sigma0, 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:

Moment-based: qsys_mxm1(lambda_batch, mu, E_X, E_X2)
PMF-based: qsys_mxm1(lambda_batch, mu, batch_sizes, pmf)
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:

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)[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:

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

qsys_ldps_workload(lambd, B, alpha, N, t=None, ngrid=None)[source]

Stationary distribution of the quantity of work in a single-stage load-dependent processor sharing station with Poisson arrivals and blocking.

Assumed model (Cohen 1979, Sect. 9; the model of Sect. 7 with one stage): a single service stage fed by a Poisson arrival stream of rate lambd; a blocking capacity N, so that a request arriving when N requests are already present is lost and leaves no trace on the state; generalized processor sharing, so that when x requests are present each accrues service at rate f(x) and the stage completes work at total rate x*f(x); and required service times i.i.d. with absolutely continuous distribution B of finite mean beta. The station is parametrized by the LINE load-dependent total rate scaling alpha(x)=x*f(x), the argument of setLoadDependence at a PS station.

With psi the total amount of service still to be given to the requests present, Cohen eqs. (9.1)-(9.3) give

Pr{psi < y} = sum_{h=0}^{N} p_h Psi^{h*}(y) p_h = (rho^h/h!) phi(h) / sum_k (rho^k/k!) phi(k), rho = lambd*beta phi(h) = 1/prod_{k=1}^{h} f(k), phi(0)=1 Psi(y) = int_0^y (1-B(v))/beta dv

with Psi^{h*} the h-fold convolution of Psi and Psi^{0*} degenerate at zero. Substituting f(k)=alpha(k)/k the factorial cancels, leaving p_h proportional to rho^h/prod_{k=1}^{h} alpha(k), the familiar load-dependent birth-death form. Psi is the equilibrium (residual life) distribution of B, so psi is a mixture of h-fold convolutions of residual service times with an atom p_0 at zero.

This is the model of Cohen (1979) Sect. 9 only. It is not the weighted GPS/DPS discipline of SchedStrategy.GPS, whose per-class weights this formula does not represent.

Parameters:

lambd – rate of the Poisson arrival stream (finite, positive)
B – required service time Distribution (continuous, finite positive mean)
alpha – rate scaling alpha(n)=n*f(n) for n=1..N (finite, positive)
N – blocking capacity (finite positive integer)
t – optional grid at which the CDF is returned. Default: an automatically sized grid covering the bulk of the distribution.
ngrid – optional number of points of the internal uniform quadrature grid on which the convolutions are formed. Default 2001. Accuracy is second order in the step for an absolutely continuous B, the case Cohen assumes, and falls back to first order when B has an atom so that 1-B is discontinuous; raise ngrid for those.

Returns:

(F, t, p) where F[j] = Pr{psi <= t[j]}, t is the grid F is reported on, and p[h] = Pr{x = h} for h = 0..N. Note F[0] = p[0] = Pr{psi = 0} when t[0] = 0, since the workload has an atom at zero.