Queueing Systems 

qsys_mapphc(D0, D1, alpha, S, c, max_num_comp=500, num_w_moms=3, w_points=None)[source]

Analyze a MAP/PH/c FCFS queue exactly.

Parameters:

D0 – arrival MAP of order ma
D1 – arrival MAP of order ma
alpha – PH service initial vector (ms,)
S – PH service sub-generator (ms, ms)
c (int) – number of servers, identical so the service law is shared
max_num_comp (int) – cap on the queue length probabilities returned
num_w_moms (int) – how many waiting-time moments to return
w_points (Sequence[float] | None) – times at which to evaluate P(Wq > t)

Returns:

MapPhcResult

Return type:

MapPhcResult

qsys_mmapgk1(MMAP, svc, w_points=None, num_w_moms=3, tol=1e-12, iter_max=10000)[source]

Analyze an MMAP[K]/G[K]/1 FCFS queue.

Parameters:

MMAP (Sequence) – LINE convention [D0, D1, D^(1), …, D^(K)] with D1 = sum_k D^(k)
svc (Sequence) – K service laws, each a LINE Distribution, a [D0, D1] phase-type pair, or a dict {‘lst’: handle, ‘moments’: [E[S], E[S^2], …]} with at least num_w_moms+1 moments
w_points (Sequence[float] | None) – times at which to evaluate the per-type waiting time CDF
num_w_moms (int) – how many per-type waiting time moments to return

Returns:

MmapGk1Result

Return type:

MmapGk1Result

class MmapGk1Result(lambdas, arrivalRate, utilization, idleVector, waitLST, waitMoments, meanWaitingTime, meanSojournTime, meanQueueLength, waitCDF, waitPoints, analyzer)[source]

Bases: object

Result structure for the MMAP[K]/G[K]/1 analysis.

lambdas: ndarray

arrivalRate: float

utilization: float

idleVector: ndarray

waitLST: Callable[[complex], ndarray]

waitMoments: ndarray

meanWaitingTime: ndarray

meanSojournTime: ndarray

meanQueueLength: float

waitCDF: ndarray | None

waitPoints: ndarray | None

analyzer: str

class MapPhcResult(meanQueueLength, meanWaitingTime, meanSojournTime, utilization, queueLengthDist, waitingTimeMoments, waitingTimeCCDF, waitingTimePoints, probWait, phaseCount, analyzer)[source]

Bases: object

Result structure for the exact MAP/PH/c analysis.

meanQueueLength: float

meanWaitingTime: float

meanSojournTime: float

utilization: float

queueLengthDist: ndarray

waitingTimeMoments: ndarray

waitingTimeCCDF: ndarray | None

waitingTimePoints: ndarray | None

probWait: float

phaseCount: int

analyzer: str

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

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

Parameters:

D0 (ndarray) – MAP hidden transition matrix.
D1 (ndarray) – MAP arrival transition matrix.
s (float) – Deterministic service time.
max_num_comp (int) – Maximum number of queue length components.
num_steps (int) – Number of waiting time points per interval.

Returns:

Performance metrics (see qsys_mapdc).

Return type:

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:

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:

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:

qsys_bmapm1(D, mu, uniformization=None, max_iter=10000, tolerance=1e-12, max_level=None, tail_tolerance=1e-10)[source]

Analyze a single-server queue fed by a batch Markovian arrival process and with exponential service of rate mu.

Parameters:

D (Sequence[ndarray]) – list of BMAP matrices [D0, D1, …, DK]. D0 carries the hidden transitions, Dk (k >= 1) the transitions that release a batch of k customers.
mu (float) – exponential service rate.
uniformization (float | None) – uniformization constant q used to randomize the generator into a discrete-time M/G/1-type chain. It must dominate every total outflow rate; by default it is chosen as max_i(-D0[i,i]) + mu.
max_iter (int) – maximum functional iterations for G.
tolerance (float) – convergence tolerance for G.
max_level (int | None) – level truncation used for the queue-length distribution (default: adaptive).
tail_tolerance (float) – relative truncation target for the level distribution.

Beyond the usual performance measures the result exposes the intermediate matrix-analytic quantities themselves, so that the algorithm can be inspected and taught rather than only its output:

theta        - stationary vector of the BMAP phase process, sum_k D_k
lambda       - mean arrival rate, theta * sum_k k*D_k * e
rho          - offered load lambda/mu
q            - uniformization constant actually used
A0, A1, Bk   - randomized blocks: A0 = (mu/q)I is a service completion
               (level down by one), A1 = (1/q)(D0 - mu*I) + I keeps the
               level, Bk[k] = (1/q)D_(k+1) raises the level by k+1
B0           - boundary local block (1/q)D0 + I, used at level 0 where
               no service can complete
A            - A0 + A1 + sum_k Bk[k], the phase process of the chain
alpha        - stationary vector of A
G            - minimal non-negative solution of
               G = A0 + A1*G + sum_k Bk[k]*G^(k+1)
drift        - alpha*(sum_k k*Bk[k])*e - alpha*A0*e. The queue is
               stable iff this is strictly negative
decayRate    - geometric decay rate of the level probabilities,
               measured as the limiting ratio pi_(n+1)/pi_n. Reported
               rather than derived from a spectral convention so that
               it is unambiguous
levelProb    - level probabilities pi_n as rows (level 0 first)
pi0          - probability the system is empty (equals 1-rho exactly)

Example:

# Example 6.4 of Bolch et al.
D0 = np.array([[-2, 0.5], [1/3, -3]])
D1 = np.array([[0.25, 0.5], [1/3, 1.0]])
D2 = np.array([[0.25, 0.5], [1.0, 1/3]])
result = qsys_bmapm1([D0, D1, D2], 11)

See also qsys_mapm1, qsys_mapph1, qsys_bmapphnn_retrial.

qsys_mm1_ps(lam, mu)[source]

Sojourn-time moments of the multiclass M/M/1-PS queue.

Class j arrives in a Poisson stream of rate lam[j] and requires an exponential amount of service with rate mu[j]. The processor is shared equally by all jobs in service, so the class of a job affects its sojourn time both through its own service rate and through the mix of rates of the jobs it shares the processor with. With alpha = 1 - sum_j lam[j]/mu[j] the unutilized fraction of the processor, the moments of the sojourn time W_r of a tagged class-r job are:

E[W_r]   = 1/(alpha*mu[r])
E[W_r^2] = 2/(alpha*mu[r])**2
           * (1 - sum_j lam_j (mu_j-mu_r)/(mu_j(mu_j+mu_r)))
           / (1 - sum_j lam_j/(mu_j+mu_r))

which is equation (7) of Mitra and Morrison (1983). Both are exact, not asymptotic: the open system is the N -> infinity limit of the closed terminal-driven system whose moments that paper expands in 1/N, and the leading term of the expansion is exact in the limit. For a single class the second moment reduces to the classical 4/(mu^2 (1-rho)^2 (2-rho)) of Coffman, Muntz and Trotter (1970).

Parameters:

lam (array_like (R,)) – Per-class Poisson arrival rates, non-negative.
mu (array_like (R,)) – Per-class exponential service rates, positive.

Returns:

W (np.ndarray (R,)) – Per-class mean sojourn times.
W2 (np.ndarray (R,)) – Per-class second moments of the sojourn time.
alpha (float) – Unutilized fraction of the processor, 1 - sum_j lam_j/mu_j.

qsys_mg1_ps(lam, svc, svcparam, x=None, s=None, t=None, nterms=41, pdf=None)[source]

Sojourn time distribution of the M/G/1-PS queue.

Parameters:

lam (float) – Poisson arrival rate, positive.
svc (array_like (n,) or callable) – Phase-type initial probability vector alpha, or a callable bhat(tau) returning the service LST, which must accept complex arguments.
svcparam (array_like (n,n) or float) – Phase-type subgenerator T when svc is a vector, or the mean service time m1 when svc is a callable.
x (array_like, optional) – Service requirements to condition on.
s (array_like, optional) – Transform arguments at which to tabulate the LST.
t (array_like, optional) – Times at which to evaluate the sojourn time distribution.
nterms (int) – Function evaluations per numerical Laplace inversion, odd.
pdf (callable, optional) – Service density, needed to remove the conditioning when svc is a callable. Filled in automatically on the phase-type path.

Returns:

With keys rho, m1, m2, lstCond, lstExcess, lstUncond, dominantRoot, x, s, t, lstCondVal, lstUncondVal, atomCond, atomUncond, meanCond, m2Cond, varCond, meanUncond, m2Uncond, varUncond, pdfCond, cdfCond, pdfUncond, cdfUncond.

Return type:

dict

qsys_mm1(lambda_val, mu)[source]

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

Parameters:

lambda_val (float) – Arrival rate (lambda)
mu (float) – Service rate

Returns:

Performance measures including:

L: Mean number in system
Lq: Mean number in queue
W: Mean response time (time in system)
Wq: Mean waiting time (time in queue)
rho: Utilization (lambda/mu)

Return type:

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:

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:

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:

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:

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:

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:

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_gigk_rqt(lambda_val, mu, Gamma_a, Gamma_s, k=1, alpha_a=2.0, alpha_s=2.0)[source]

Robust Queueing Theory worst-case system time of a G/G/k FCFS queue.

The uncertainty sets are

U^a = {T : (sum_{i=k+1}^n T_i - (n-k)/lambda)/(n-k)^(1/alpha_a) >= -Gamma_a} U^s = {X : (sum_{i=k}^n X_i - (n-k+1)/mu)/(n-k+1)^(1/alpha_s) <= Gamma_s}

with alpha=2 the finite-variance regime and alpha in (1,2) the heavy-tailed one. The returned W is the closed-form bound of Theorem 3 (Theorem 8 when the two tail coefficients differ, with alphabar = min(alpha_a,alpha_s)),

W <= (ab-1)/ab^(ab/(ab-1)) lambda^(1/(ab-1))
(Gamma_a+Gamma_s/k^(1/ab))^(ab/(ab-1)) / (1-rho)^(1/(ab-1)) + k/lambda,

which for k=1 reduces to Theorem 2 and, at alphabar=2, to the Kingman-like form (lambda/4)(Gamma_a+Gamma_s)^2/(1-rho) + 1/lambda. Sworst is the exact worst case over the uncertainty sets, eq. (45), the supremum over the integer x >= 1 of

x/mu + Gamma_s x^(1/alpha_s) - k(x-1)/lambda + Gamma_a (k(x-1))^(1/alpha_a).

The arrival deviation ADDS to the worst case, since the adversary shortens the interarrival times: the sign printed in eq. (12) is easily misread as a subtraction of the whole arrival bracket.

W is a SYSTEM time (waiting plus service), and its additive term is k/lambda rather than the mean service time 1/mu.

Parameters:

lambda_val (float) – Arrival rate
mu (float) – Service rate of each server
Gamma_a (float) – Variability parameter of the arrival uncertainty set
Gamma_s (float) – Variability parameter of the service uncertainty set
k (int) – Number of servers
alpha_a (float) – Arrival tail coefficient in (1,2]
alpha_s (float) – Service tail coefficient in (1,2]

Returns:

Tuple of (W, rhohat, Sworst)

Return type:

Tuple[float, float, float]

qsys_gig1_rqt(lambda_val, mu, Gamma_a, Gamma_s, alpha_a=2.0, alpha_s=2.0)[source]

Robust Queueing Theory worst-case system time of a G/G/1 FCFS queue, the single-server case of qsys_gigk_rqt (Theorem 2 and eq. 12).

Parameters:

lambda_val (float) – Arrival rate
mu (float) – Service rate
Gamma_a (float) – Variability parameter of the arrival uncertainty set
Gamma_s (float) – Variability parameter of the service uncertainty set
alpha_a (float) – Arrival tail coefficient in (1,2]
alpha_s (float) – Service tail coefficient in (1,2]

Returns:

Tuple of (W, rhohat, Sworst)

Return type:

Tuple[float, float, float]

qsys_gigk_rqt_gamma(rho, mu, Gamma_a, sigma_s, k=1, alpha_a=2.0, regime='independent')[source]

Service variability parameter of the RQT framework, from the first two moments, by the adaptation of Section 7.1:

Gamma_s = (2 (theta0 + theta1 sigma_s^2/k + theta2 Gamma_a^2 rho^2 k))^((a-1)/a)
          - Gamma_a k^((a-1)/a)

where (theta0,theta1,theta2) are regressed so that the worst-case system time of Theorem 3 approximates the MEAN system time of the corresponding stochastic queue. The arrival side needs no adaptation: Gamma_a = sigma_a for an external renewal stream. Since the last term cancels Gamma_a at alpha=2, the adaptation acts on the sum Gamma_a + Gamma_s/k^(1/alpha) that Theorem 3 reads.

THE FACTOR 2 IS NOT IN THE PRINTED FORMULA and is restored here. Section 7.1 states that the form is motivated by Kingman’s bound, which the alpha=2 bound of Theorem 3 reproduces when (Gamma_a+Gamma_s)^2 = 2(sigma_a^2+sigma_s^2); the published thetas are all near unity, i.e. corrections to that bound rather than a substitute for its factor 2. Dropping the factor puts M/M/1 about 40% BELOW its exact mean system time at rho=0.9, contradicting the errors of at most 9.5% that Tables 2-3 report; restoring it gives +4.7%.

CAUTION: the form is not dimensionally homogeneous, since theta0 is an additive constant on a scale of variances, so it is only valid in the time unit the regression was run in. It is evaluated here in units of the mean service time, 1/mu = 1, and converted back.

Parameters:

rho (float) – Traffic intensity lambda/(k*mu)
mu (float) – Service rate of each server, which sets the time unit
Gamma_a (float) – Variability parameter of the arrival uncertainty set
sigma_s (float) – Standard deviation of the service time
k (int) – Number of servers
alpha_a (float) – Effective arrival tail coefficient in (1,2]
regime (str) – Adaptation regime of Table 1, ‘independent’, ‘normal’ or ‘pareto’

Returns:

The service variability parameter Gamma_s

Return type:

float

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/java/jline/api/qsys/Qsys_gig1_lbnd.java

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 (ndarray) – Vector of arrival rates per priority class (class 1 = highest)
mu_vec (ndarray) – Vector of service rates per priority class
cs_vec (ndarray) – Vector of coefficients of variation per priority class

Returns:

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

Return type:

Tuple of (W, rho)

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 (ndarray) – Per-class Poisson arrival rates (K,)
mu_vec (ndarray) – Per-class exponential service rates (K,)
w_vec (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 (ndarray) – Vector of arrival rates per class
mu_vec (ndarray) – Vector of service rates per class
cs_vec (ndarray) – Vector of coefficients of variation per class

Returns:

W: Vector of mean response times per class (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 (ndarray) – Vector of arrival rates per class
mu_vec (ndarray) – Vector of service rates per class
cs_vec (ndarray) – Vector of coefficients of variation per class

Returns:

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

Return type:

Tuple of (W, rho)

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 (ndarray) – Vector of arrival rates per class
mu_vec (ndarray) – Vector of service rates per class
cs_vec (ndarray) – Vector of coefficients of variation per class

Returns:

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

Return type:

Tuple of (W, rho)

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 (ndarray) – Vector of arrival rates per class
mu_vec (ndarray) – Vector of service rates per class
cs_vec (ndarray) – Vector of coefficients of variation per class

Returns:

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

Return type:

Tuple of (W, rho)

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 (ndarray) – Vector of arrival rates per class
mu_vec (ndarray) – Vector of service rates per class
cs_vec (ndarray) – Vector of coefficients of variation per class

Returns:

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

Return type:

Tuple of (W, rho)

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=0 : if 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

qsys_mapg1k(D0, D1, svc, K, tol=1e-12, nmax=200000)[source]

Exact analysis of a MAP/G/1/K queue with tail drop.

Markovian arrivals, arbitrary service time distribution F, and a finite buffer of K packets (the position held by the packet in transmission included). Unlike qsys_mapg1 the service time is NOT fitted to a phase-type distribution: F enters exactly, through the functionals A_m and Q_m evaluated by uniformization of the arrival MAP.

Parameters:

D0 – MAP parameter matrices (M x M), D0 + D1 an irreducible generator
D1 – MAP parameter matrices (M x M), D0 + D1 an irreducible generator
svc (Dict[str, Any]) – service time descriptor dict with key ‘type’: ‘gamma’ : keys alpha (shape), theta (scale) ‘det’ : key d (constant service time) ‘ph’ : keys alpha (1 x p), T (p x p subgenerator) ‘density’: key pdf (callable), optional key tmax
K (int) – buffer size in packets, K >= 1
tol (float) – uniformization truncation tolerance
nmax (int) – cap on the uniformization order

Returns:

Dict with p0, pK, lossProbability, throughput, lambda, meanServiceTime, utilization, rho, nmax, sigma, pKvec, p0vec, plevel, meanQueueLength.

Return type:

Method:

The chain embedded at departure epochs is used, in the state (n,j): n = 0..K-1 packets left behind by a departure, j = MAP phase. With A_m the matrix of “m arrivals during a service, phase i -> j”,

n >= 1: n’ = n-1+min(m, K-n), overflow sum_{m>=K-n} A_m n == 0: the phase first jumps by (-D0)^{-1}*D1 (the idle period ends at an arrival), the service then proceeds as from n=1.

Its stationary law sigma gives, by Markov renewal reward, the cycle mean, p0 and pK, where Q_m is the expected time within a service with exactly m arrivals so far. Time-stationary p0 and pK follow, so no PASTA assumption is needed on the MAP side.

References

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

qsys_mmapg1k(D0, D1c, svc, K, tol=1e-12, nmax=200000)[source]

Exact per-class throughput and loss ratio of an MMAP[K]/G/1/K queue.

Two classes of equal arrival rate but different interarrival variability or autocorrelation receive different loss ratios. Aggregate-only finite-buffer analyses cannot express it: they return a single blocking probability p and set T_k = lambda_k*(1-p), making the loss ratio identical by construction.

Parameters:

D0 – M x M hidden transition matrix of the arrival MMAP
D1c (Sequence) – sequence of R matrices, D1c[k] = M x M arrival matrix of class k
svc (Dict[str, Any]) – service time descriptor, see qsys_mapg1k
K (int) – buffer size in packets, K >= 1

Returns:

Dict with throughput, lossRatio, lambda (all per class), the aggregate quantities, and the level/phase quantities of the driving MAP model.

Return type:

Method:

The aggregate MAP {D0, sum_k D1c[k]} drives qsys_mapg1k, whose embedded chain returns the joint law of buffer level and MAP phase. A class-k arrival leaves phase i at rate (D1c[k]*e)_i, so

lambda_k = pi*D1c[k]*e, L_k = (pKvec*D1c[k]*e)/lambda_k.

This is exact: no independence between classes is assumed and no PASTA argument is used, the phase resolution of pKvec doing the work.

Assumes a single server and a service law that is iid and independent of class.

References

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

qsys_mapg1k_perflow(MAPS, svc, K, tol=1e-12, nmax=200000)[source]

Per-flow throughput and loss ratio of a FIFO buffer fed by N flows.

Flow n is described by its own MAP, so two flows may share an arrival rate and still differ in the shape and autocorrelation of their interarrival times. The buffer holds K packets including the one in transmission.

Parameters:

MAPS (Sequence) – sequence of N pairs (D0n, D1n); the orders M_n may differ
svc (Dict[str, Any]) – service time descriptor, see qsys_mapg1k
K (int) – buffer size in packets, K >= 1

Returns:

Dict with throughput, lossRatio, lambda, p0 and pK per flow, plus the aggregate quantities and rho.

Return type:

Method:

The exact model of N flows would need prod_n M_n * (K+1) states. Instead one model per flow is solved: flow n is kept exactly as MAP_n while the other N-1 flows are replaced by a single Poisson stream of rate lambda - lambda_n, justified by the Palm-Khinchin limiting theorem on the superposition of many point processes. The superposition yields

D0 = D0n - lambdaBar_n*I, D1 = D1n + lambdaBar_n*I,

which is passed to qsys_mapg1k. The sweep is O(N*(K*M)^3) against the O(M^(3N)*K^3) of the exact joint model.

References

Original MATLAB: matlab/src/api/qsys/qsys_mapg1k_perflow.m Chydzinski, A. Applied System Innovation 2026, 9, 112, Theorem 1.

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

Bases: object

Result structure for queue analysis.

meanQueueLength: float

meanWaitingTime: float

meanSojournTime: float

utilization: float

queueLengthDist: ndarray | None = None

queueLengthMoments: ndarray | None = None

sojournTimeMoments: ndarray | None = None

analyzer: str = 'native'

ph_to_map(alpha, T)[source]

Convert a PH distribution to its equivalent MAP representation.

For a PH renewal process, the MAP has:: D0 = T (transitions within the PH, no arrival) D1 = t * alpha where t = -T*e (exit rates times restart distribution)

Parameters:

alpha (ndarray) – Initial probability vector (1 x n)
T (ndarray) – Sub-generator matrix (n x n)

Returns:

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

Return type:

Tuple of (D0, D1)

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

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

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

Parameters:

alpha (ndarray) – Arrival PH initial probability vector (1 x n)
T (ndarray) – Arrival PH generator matrix (n x n)
beta (ndarray) – Service PH initial probability vector (1 x m)
S (ndarray) – Service PH generator matrix (m x m)
numQLMoms (int) – Number of queue length moments to compute (default: 3)
numQLProbs (int) – Number of queue length probabilities (default: 100)
numSTMoms (int) – Number of sojourn time moments (default: 3)

Returns:

QueueResult with queue performance metrics

Return type:

References

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

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

Analyze a MAP/PH/1 queue.

Parameters:

D0 (ndarray) – MAP hidden transition matrix
D1 (ndarray) – MAP observable transition matrix
beta (ndarray) – Service PH initial probability vector
S (ndarray) – Service PH generator matrix
numQLMoms (int) – Number of queue length moments
numQLProbs (int) – Number of queue length probabilities
numSTMoms (int) – Number of sojourn time moments

Returns:

QueueResult with queue performance metrics

Return type:

References

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

qsys_mapm1(D0, D1, mu)[source]

Analyze a MAP/M/1 queue.

Parameters:

D0 (ndarray) – MAP hidden transition matrix
D1 (ndarray) – MAP observable transition matrix
mu (float) – Service rate

Returns:

QueueResult with queue performance metrics

Return type:

References

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

qsys_mapmc(D0, D1, mu, c, max_num_comp=1000)[source]

Analyze a MAP/M/c queue.

The chain is a level-dependent QBD in the number in system: above level c the c servers are all busy, so the blocks repeat as A0 = c*mu*I (down), A1 = D0 - c*mu*I (local), A2 = D1 (up) and the tail is matrix-geometric in R. Below level c the departure rate is level dependent and the boundary vector comes from the Gaver, Jacobs and Latouche backward recursion. The waiting time is phase type; its generator is the fixed point of a Sylvester equation, and the representation is the time reversal with respect to the arrival-epoch vector.

Parameters:

D0 (ndarray) – MAP hidden transition matrix
D1 (ndarray) – MAP observable transition matrix
mu (float) – Service rate per server
c (int) – Number of servers
max_num_comp (int) – Cap on the number of queue length probabilities

Returns:

QueueResult with queue performance metrics

Return type:

References

Perez, Van Velthoven, Van Houdt, Q-MAM, ValueTools 2008 (Q_CT_MAP_M_C) Gaver, Jacobs, Latouche, Adv. Appl. Probab. 16:715-731, 1984 Asmussen, Moller, Queueing Systems 37(1):9-29, 2001 Original MATLAB: matlab/src/api/qsys/qsys_mapmc.m

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

Analyze a MAP/MAP/1 queue.

Both arrival and service processes are Markovian Arrival Processes.

Parameters:

D0_arr (ndarray) – Arrival MAP hidden transition matrix
D1_arr (ndarray) – Arrival MAP observable transition matrix
D0_srv (ndarray) – Service MAP hidden transition matrix
D1_srv (ndarray) – Service MAP observable transition matrix

Returns:

QueueResult with queue performance metrics

Return type:

References

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

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

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

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

Parameters:

D0 (ndarray) – MAP hidden transition matrix (n x n)
D1 (ndarray) – MAP arrival transition matrix (n x n)
service_moments (ndarray) – First k raw moments of service time [E[S], E[S^2], …] (k = 2 or 3 for best accuracy)
num_ql_moms (int) – Number of queue length moments to compute (default: 3)
num_ql_probs (int) – Number of queue length probabilities (default: 100)
num_st_moms (int) – Number of sojourn time moments to compute (default: 3)

Returns:

QueueResult with queue performance metrics

Return type:

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

D0: ndarray

D_batch: List[ndarray]

__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.

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

beta: ndarray

S: ndarray

__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

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

A0: ndarray

A1: ndarray

A2: ndarray

B0: ndarray

max_level: int

__post_init__()[source]: Validate QBD state space.

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

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]

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.

class RetrialQueueResult(queue_type, L_orbit, N_server, utilization, throughput, P_idle, P_empty_orbit, stationary_dist=None, truncation_level=0, trunc_error=inf, converged=False, iterations=0, error=inf, analyzer='LINE:qsys_bmapphnn_retrial')[source]

Bases: object

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

Variables:

queue_type (line_solver.api.qsys.retrial.QueueType) – Type of queue (retrial, reneging, etc.)
L_orbit (float) – Expected number of customers in orbit
N_server (float) – Expected number of customers being served
utilization (float) – Server utilization
throughput (float) – System throughput
P_idle (float) – Probability that server is idle
P_empty_orbit (float) – Probability that orbit is empty
stationary_dist (numpy.ndarray | None) – Stationary distribution vector
truncation_level (int) – QBD truncation level used
trunc_error (float) – Relative orbit-truncation error estimate at truncation_level
converged (bool) – Whether numerical solution converged
iterations (int) – Number of iterations to convergence
error (float) – Final error estimate
analyzer (str) – Name of the engine that produced the result

queue_type: QueueType

L_orbit: float

N_server: float

utilization: float

throughput: float

P_idle: float

P_empty_orbit: float

stationary_dist: ndarray | None = None

truncation_level: int = 0

trunc_error: float = inf

converged: bool = False

iterations: int = 0

error: float = inf

analyzer: str = 'LINE:qsys_bmapphnn_retrial'

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

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

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

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

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, 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, 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': fixed orbit truncation level. When None or
    non-positive (default) the level is chosen adaptively: it is
    doubled until the mass retained at the top level contributes
    less than 'TailTolerance' of the mean orbit length. A fixed
    level disables the adaptive refinement.

'Tolerance': convergence tolerance (default: 1e-10).
'TailTolerance': relative orbit-truncation error target
    (default: 1e-6).

'MaxDim': cap on the total generator dimension explored by the
    adaptive refinement (default: 2e5).

'MaxBlockSize': cap on the per-level block size V*d
    (default: 5000). Exceeding it is an error: the phase-type
    service order and the server count make the level block
    intractable.

'RetrialPolicy': RetrialPolicy.LINEAR (default), where the
    aggregate retrial rate is (orbit size)*alpha, or
    RetrialPolicy.CONSTANT, where the orbit retries as a whole at
    rate alpha whenever it is non-empty.

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

is_retrial: bool

station_idx: int | None = None

node_idx: int | None = None

source_idx: int | None = None

class_idx: int | None = None

error_msg: str = ''

N: int | None = None

alpha: float = 0.1

gamma: float = 0.0

p: float = 0.0

R: int | None = None

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.

is_reneging: bool

source_idx: int | None = None

queue_idx: int | None = None

class_idx: int | None = None

n_servers: int | None = None

service_rate: float | None = None

error_msg: str = ''

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[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[ndarray | None, 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:

'tol': Convergence tolerance (default 1e-10)
'verbose': Print progress messages (default False)
'config': Dict with 'mapmsg_quantization' (default 11),
    'orbit_maxlevel' (fixed orbit truncation level; default None,
    meaning the engine chooses it adaptively) and 'orbit_tailtol'
    (relative orbit-truncation error target, default 1e-6)

Returns:

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

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
perf: the engine result object, exposing the matrix-analytic
    internals (see SolverMAM.getMAMResult)

Return type:

Tuple[ndarray, ndarray, ndarray, ndarray, ndarray, ndarray, int, Any]

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.

qsys_lindley_moment(lambda_, mu, Wn, m)[source]

One conditional Lindley moment for exponential primitives.

Evaluates E[max(Wn + S - A, 0)**m] with A ~ Exp(lambda_) and S ~ Exp(mu). This is the algorithm shared by qsys_mm1_lindley(), which calls it once per moment order, and qsys_hh1_lindley(), which mixes it over the arrival and service phases.

See qsys_mm1_lindley() for the derivation and for why the upper incomplete gamma function reduces to a finite sum here.

Parameters:

lambda – Arrival rate, positive
mu (float) – Service rate, positive
Wn – Current waiting times
m (int) – Moment order, at least 1

Returns:

The conditional moment at each entry of Wn.

Return type:

ndarray

qsys_mm1_lindley(lambda_, mu, Wn, mmax=2)[source]

Conditional waiting-time moments of the M/M/1 Lindley recursion.

One step of W_{n+1} = max(W_n + S_n - A_n, 0) with A_n ~ Exp(lambda_) and S_n ~ Exp(mu): given the waiting time of customer n, the exact conditional moments of the waiting time of customer n+1.

The m-th conditional moment is

E[W_{n+1}^m | W_n] = lambda mu/(lambda+mu) [ S + T ], S = sum_{k=0}^{m} C(m,k) W_n^k (m-k)! / mu^(m-k+1), T = (-1)^m e^{-lambda W_n} (Gamma(m+1,-lambda W_n) - m!) / lambda^(m+1),

where the density of S_n - A_n is the asymmetric Laplace density lambda mu/(lambda+mu) times e^{-mu x} for x > 0 and e^{lambda x} for x < 0. Because m+1 is a positive integer, the upper incomplete gamma function admits the finite form Gamma(m+1,x) = m! e^{-x} sum_{k=0}^{m} x^k/k!, valid at the negative argument -lambda W_n needed here. Substituting it cancels the growing exponential and leaves the numerically stable

T = (-1)^m m! ( sum_{k=0}^{m} (-lambda W_n)^k/k! - e^{-lambda W_n} ) / lambda^(m+1),

which is what this function evaluates. No incomplete gamma routine is needed.

The mean is returned from the equivalent explicit form W_n + (lambda-mu)/(lambda mu) + mu e^{-lambda W_n}/(lambda(lambda+mu)), and the variance as the second moment less the squared mean.

Verified against 4e6 Monte Carlo replications to 5e-4 relative error for m = 1, 2, 3.

Parameters:

lambda – Arrival rate, positive
mu (float) – Service rate, positive
Wn – Current waiting times, finite and nonnegative
mmax (int) – Highest moment order, raised to 2 so the variance is available

Returns:

Dict with keys mean, var, moments (shape (len(Wn), mmax)), mmax and analyzer.

Raises:

ValueError – If a rate is not positive or a waiting time is negative.

Return type:

qsys_hh1_lindley(lambda_, pa, mu, ps, Wn, mmax=2)[source]

Conditional waiting-time moments of the Hl/Hn/1 Lindley recursion.

Hyperexponential primitives are mixtures of exponentials, so conditioning on the arrival phase i and the service phase j reduces one Lindley step to the M/M/1 step of qsys_mm1_lindley() at rates lambda_[i] and mu[j], and the conditional moment is the corresponding mixture

E[W_{n+1}^m | W_n] = sum_i sum_j pa[i] ps[j] E_{ij}[W_{n+1}^m | W_n].

Phases are drawn independently for each customer, which is what makes the mixture exact rather than an approximation; a Markov-modulated arrival stream would not decompose this way.

Note that the variance is not the corresponding mixture of the per-phase variances, because the phase is itself random: it is recovered here from the first two mixed raw moments, which adds the between-phase spread of the means.

Verified against 4e6 Monte Carlo replications to 2e-3 relative error for m = 1, 2.

Parameters:

lambda – Arrival phase rates, positive
pa (Sequence[float]) – Arrival phase probabilities, nonnegative and summing to 1
mu (Sequence[float]) – Service phase rates, positive
ps (Sequence[float]) – Service phase probabilities, nonnegative and summing to 1
Wn – Current waiting times, finite and nonnegative
mmax (int) – Highest moment order, raised to 2 so the variance is available

Returns:

The same dict as qsys_mm1_lindley().

Raises:

ValueError – If the phase vectors are inconsistent or unnormalized.

Return type:

qsys_mm1_tandem_lindley(lambda_, mu1, mu2, Wk, Wk1)[source]

Conditional downstream waiting time in an M/M/1 tandem.

The conditional mean waiting time of customer n+1 at the downstream station of a two-station single-server tandem queue, given that customer n waited Wk upstream and Wk1 downstream.

The point of the tandem recursion is that the interarrival time at the downstream station is the interdeparture time upstream, not an independent draw. With A ~ Exp(lambda_) the interarrival time upstream and S1, S1' the service times upstream of customers n and n+1, that interdeparture time is

D = max(A - Wk - S1, 0) + S1’,

an idle period followed by the next service, and the downstream Lindley step is W2_{n+1} = (Wk1 + S2 - D)^+ with S2 ~ Exp(mu2) independent of D.

Because A is exponential, max(A - Wk - S1, 0) is zero with probability 1-q and Exp(lambda_) with probability q = e^{-lambda Wk} mu1/(lambda+mu1), the probability the upstream server goes idle, so D is either Exp(mu1) or the sum of Exp(mu1) and Exp(lambda_). Averaging the downstream step over both cases needs only two elementary transforms of:

g(d) = E[(Wk1 + S2 - d)^+] = Wk1 - d + 1/mu2      for d <= Wk1,
                           = e^{-mu2 (d-Wk1)}/mu2  for d > Wk1,

namely J(c) = int_0^inf e^{-cu} g(u) du and Jw(c) = int_0^inf u e^{-cu} g(u) du, both closed form, giving:

E[W2_{n+1} | Wk, Wk1] = (1-q) mu1 J(mu1) + q C,
C = lambda mu1 (J(mu1) - J(lambda))/(lambda-mu1)  if lambda != mu1,
  = mu1^2 Jw(mu1)                                 if lambda == mu1.

As Wk grows the upstream server never idles, q vanishes, and the mean tends to mu1 J(mu1) = E[g(S1')], as it must.

Two caveats, both inherited from the reference and both quantified here.

First, this is exact for the step taken in isolation, that is when the conditioning pair is independent of the four primitives that drive the step. In a running tandem it is not: the downstream wait Wk1 was itself determined by an interdeparture time containing S1, so conditioning on (Wk, Wk1) is not conditioning on a Markov state of the tandem. Measured against a 4e6-customer simulation of the real tandem at lambda = 0.8, mu1 = mu2 = 1, the formula is within 0.4% to 1.3% away from the empty state and 7% at Wk = Wk1 = 0, where the entanglement is strongest. Treat it as exact for one isolated step and as a good approximation in a running tandem.

Second, this closed form was derived rather than transcribed from the reference’s theorem 4, because that theorem rests on its proposition 1, which omits a service-time difference and so does not describe a tandem queue; see qsys_tandem_lindley(). The two differ: at lambda = 0.8, mu1 = mu2 = 1 and Wk = Wk1 = 0 the published route gives 0.3016 against 0.3457 here, the latter matching simulation of the step to 6e-4 relative error.

As in the reference, the upstream interarrival time is taken to be Exp(lambda_), which by Burke’s theorem is also the stationary interdeparture law, so the same formula is applied at any pair of consecutive stations of a longer M/M/1 tandem, with the caveat above compounding.

Parameters:

lambda – External arrival rate upstream, positive
mu1 (float) – Upstream service rate, positive
mu2 (float) – Downstream service rate, positive
Wk – Current upstream waiting times, finite and nonnegative
Wk1 – Current downstream waiting times, same shape as Wk or scalar

Returns:

Dict with keys mean, interdepMean, idleProb and analyzer.

Raises:

ValueError – If a rate is not positive or the two wait arrays disagree.

Return type:

qsys_tandem_lindley(A, S, W0=None)[source]

Tandem network Lindley recursion on a sample path.

Propagates the waiting times of a series of K single-server FCFS stations in tandem, driven by the primitives of the sample path. A is the length-N vector of interarrival times at the first station, A[n] separating customers n and n+1, and S is the (N, K) matrix of service times, S[n, k] being the service time of customer n at station k.

At the first station this is Lindley’s recursion, W[n+1, 0] = max(W[n, 0] + S[n, 0] - A[n], 0). Downstream the interarrival time is not a primitive: the arrival epoch of customer n at station k is its departure epoch from station k-1, so the interarrival time at station k is the interdeparture time upstream. Writing G[n, k] for the interarrival time at station k between customers n and n+1, with G[n, 0] = A[n], the exact interdeparture identity is

G[n, k+1] = G[n, k] + W[n+1, k] - W[n, k] + S[n+1, k] - S[n, k],

equivalently and more transparently

G[n, k+1] = max(G[n, k] - W[n, k] - S[n, k], 0) + S[n+1, k],

an idle period at station k followed by the next customer’s service there. The recursion at station k is then W[n+1, k] = max(W[n, k] + S[n, k] - G[n, k], 0).

Nothing here is distributional, so the recursion is exact for arbitrary interarrival and service times, dependent or not, and is the reference a simulated tandem sample path can be checked against directly. It reproduces a direct event-driven tandem simulation to 1e-12 over four stations.

Note that proposition 1 of the reference states this identity without the S[n+1, k] - S[n, k] term, which makes it wrong as a sample-path identity: the omitted difference has mean zero, so the mean interdeparture time survives, but individual waiting times do not. Implementing it as published gives station-1 waiting times that are correct and downstream ones that are not, by up to several mean service times. The form above is used instead.

Parameters:

A – Interarrival times at the first station
S – Service times, shape (N, K)
W0 (Sequence[float] | None) – Waiting times of customer 1 at each station, None for an empty network

Returns:

Dict with keys W, G, T, departure and analyzer, every matrix indexed [customer, station]. The last row of G is NaN, there being no customer N+1 to separate from.

Raises:

ValueError – If the shapes disagree or a primitive is negative.

Return type:

qsys_tandem_ub_ciucu(x, lst, p, mu, dlst=None)[source]

Tail bounds for a GI/Hn/1 -> ./Hn/1 tandem of two FCFS single servers.

Both stations serve the same hyperexponential law Y, Z ~ sum_i p_i Exp(mu_i), a scalar p = 1 giving exponential service, and the arrivals are renewal with a light-tailed interarrival time supplied through its Laplace-Stieltjes transform E[e^{-s X}].

With theta the positive root of E[e^{theta (Y-X)}] = 1 and alpha = E[X e^{-theta X}], the test function

gamma(u,v) = 1{0<=u<=v} [1 - A e^{-theta u} - (B + C u + D v) e^{-theta v}]

satisfies the integral inequality of Theorem 1(b) of the reference once the five sufficient conditions of its Lemma 4 fix A, B, C and D,

A = 1, C = theta sum_i p_i/(mu_i-theta) / sum_i p_i mu_i/(mu_i-theta)^2, D = max(-C E[U e^{theta V}]/E[V e^{theta V}], 0), U = Y-X, V = Z-X, B = C (1/mu_1 - alpha E[e^{theta Z}]) if D = 0,

= (C+D)/(mu_1-theta) - theta/mu_1 if D > 0,

with mu_1 the smallest service rate. Corollary 2 then turns gamma into

P(S > x) <= sum_i p_i { e^{-mu_i x}

mu_i/(mu_i-theta) (A+B) (e^{-theta x} - e^{-mu_i x})

mu_i/(mu_i-theta)^2 (C+D) (((mu_i-theta)x-1) e^{-theta x}

e^{-mu_i x}) }

and the corresponding closed form for W when the service is exponential. E[V e^{theta V}] is positive at any stable load, so D is always well defined: h(s) = E[e^{s(Z-X)}] is convex with h(0) = h(theta) = 1, hence h'(theta) > 0.

In the M/M/1 -> ./M/1 case the five inequalities hold as equalities, so gamma is the exact joint distribution and both bounds are exact, P(S > x) = (1 + theta x) e^{-theta x}. Away from it the bound stays sharp: against an exact CTMC reference for the Erlang(2)/M/1 -> ./M/1 tandem it is within 2% at P(S>x) = 1e-2 and within 0.6% at 5e-10, with the correct asymptotic slope theta^2/(mu(1-alpha mu)). Accuracy degrades with service variability, to about a factor of two at CV(Y) = 2.

Parameters:

x – Thresholds at which the tails are bounded, nonnegative
lst (Callable[[float], float]) – Interarrival transform, s -> E[e^{-s X}] for s >= 0
p – Service phase probabilities, nonnegative and summing to one
mu – Service phase rates, positive
dlst (Callable[[float], float] | None) – s -> E[X e^{-s X}], minus the derivative of lst; None to obtain it by a Richardson-extrapolated central difference, which costs four extra transform evaluations and loses roughly four digits

Returns:

Dict with keys ‘S’ (bound on P(S>x), capped at one), ‘W’ (bound on P(W>x), NaN unless the service is exponential), ‘theta’, ‘alpha’, and the coefficients ‘A’, ‘B’, ‘C’, ‘D’.

Return type:

Example

>>> from math import exp
>>> r = qsys_tandem_ub_ciucu([5, 10], lambda s: exp(-s * 4 / 3), 1.0, 1.0,
...                          lambda s: (4 / 3) * exp(-s * 4 / 3))
>>> float(round(r['S'][0], 6))
0.493699

qsys_mgisrgi_whitt(lambda_val, mu, s, r, patience, wPoints=None, maxQueue=100000, tol=1e-14, invMethod='euler', invN=41)[source]

Engineering solution of the M/GI/s/r+GI queue.

Poisson arrivals at rate lambda_val, iid general service times of mean 1/mu, s servers, r extra waiting spaces and iid patience times with a general distribution.

The general patience law is replaced by state-dependent Markovian abandonment, a customer jth from the end of the queue abandoning at rate delta_j = h(j/lambda) for the patience hazard h (eq. 3.3), because such a customer has been waiting for about j/lambda; the general service law is replaced by an exponential of the same mean (Section 5). What is left is a birth-and-death process, solved exactly. Only the hazard NEAR THE ORIGIN matters, not the mean or the tail of the patience law.

Parameters:

lambda_val (float) – arrival rate
mu (float) – service rate of one server, the reciprocal of the mean service time
s (int) – number of servers
r (float) – extra waiting spaces, float('inf') for an unbounded queue
patience (float | Callable[[float], float] | Dict[str, Callable[[float], float]]) – scalar rate (exponential patience, then the answer is exact and the model is Erlang A), a callable hazard h(t), or a dict {'ccdf': G} using the integrated form of eq. (3.6)
wPoints (Sequence[float] | None) – times at which to return the waiting-time cdfs
maxQueue (int) – truncation level used when r is infinite
tol (float) – relative tail tolerance for that truncation
invMethod (str) – Laplace inversion method for the cdfs
invN (int) – number of inversion nodes

Returns:

Dict with the steady-state distribution queueLengthDist, the probabilities probLoss/probNoWait/probServed/probAbandon, the moments meanNumber/varNumber/meanQueueLength/ varQueueLength/meanWaitServed/varWaitServed/ meanWaitAbandon/varWaitAbandon/meanWait/secondMomentWait, the rates utilization/throughput/abandonRate, the abandonment rates abandonRates/totalAbandonRates, and, when wPoints is given, cdfWaitServed/cdfWaitAbandon/cdfWait.

Return type:

References

W. Whitt (2005). Engineering solution of a basic call-center model. Management Science 51(2), 221-235.

qsys_erlanga(lambda_val, mu, theta, s, r=float('inf'), **kwargs)[source]

Exact analysis of the Erlang A model M/M/s/r+M.

Poisson arrivals at rate lambda_val, exponential service of rate mu at each of s servers and exponential patience of rate theta. The number in system is the birth-and-death process with death rate min(k,s)*mu + (k-s)^+ * theta, so every measure is exact: this is the case in which the state-dependent Markovian approximation of qsys_mgisrgi_whitt() reproduces the model rather than approximating it (eq. 7.12 of the reference). theta = 0 recovers M/M/s/r, and then a finite r is required whenever lambda_val >= s*mu.

Parameters:

lambda_val (float) – arrival rate
mu (float) – service rate of one server
theta (float) – abandonment rate of a waiting customer
s (int) – number of servers
r (float) – extra waiting spaces, infinite by default
**kwargs – passed through to qsys_mgisrgi_whitt()

Returns:

The dict returned by qsys_mgisrgi_whitt().

Return type:

References

W. Whitt (2005). Engineering solution of a basic call-center model. Management Science 51(2), 221-235, Section 7 and eq. (7.12). The model itself is due to C. Palm (1937, 1957).

qsys_ggisgi_fluid(lambda_val, mu, s, patienceCcdf, servingCcdf=None, agePoints=None, tol=1e-12, maxTime=None)[source]

Steady state of the G/GI/s+GI fluid model.

Scale the content by s and let s grow. Customers become quanta of fluid but their sojourns do not shrink, so the ages survive the limit: the state is the density b(x) of fluid in service of age x and the density q(x) of fluid waiting of age x. With rho = lambda/(s*mu),

rho <= 1: b(x) = rho G^c(x), q = 0, no wait, no abandonment;
rho > 1: b(x) = G^c(x), q(x) = rho F^c(x) on [0,w],

the queue boundary w solving F^c(w) = 1/rho (eq. 3.6): fluid that survives its patience for w enters service, so the surviving fraction must equal the fraction 1/rho the servers can absorb.

Parameters:

lambda_val (float) – arrival rate
mu (float) – service rate of one server
s (int) – number of servers
patienceCcdf (Callable[[float], float]) – F^c(t) = P(patience > t)
servingCcdf (Callable[[float], float] | None) – G^c(x) = P(service > x), needed only for the in-service age density; defaults to the exponential of rate mu
agePoints (Sequence[float] | None) – ages at which to return the two densities
tol (float) – bisection tolerance for w
maxTime (float | None) – largest age searched for w; the search grows automatically when this is None

Returns:

Dict with regime, trafficIntensity, offeredWait, meanWait, meanWaitServed, meanWaitAbandon, probAbandon, meanQueueLength, meanNumberInService, meanNumber, utilization, throughput, abandonRate, and, when agePoints is given, serviceAgeDensity and queueAgeDensity.

Return type:

References

W. Whitt (2006). Fluid models for multiserver queues with abandonments. Operations Research 54(1), 37-54.

qsys_mmk_qed(lambda_val, mu, s)[source]

Halfin-Whitt QED approximation for the M/M/s queue.

Let s grow with the offered load a = lambda/mu so that the server slack beta = (1-rho)sqrt(s) = (s-a)/sqrt(s) stays fixed. The delay probability then has the non-degenerate limit alpha(beta): servers are busy a fraction 1 - beta/sqrt(s) of the time, so efficiency tends to 1, and yet the delay probability tends to a constant strictly between 0 and 1, so quality does not collapse.

Useful even though M/M/s is exactly solvable, because Erlang C needs a sum of s terms a^j/j! that overflows in double precision well before the thousands of servers a large contact centre or thread pool has.

Parameters:

lambda_val (float) – arrival rate
mu (float) – service rate of one server
s (int) – number of servers

Returns:

Dict with offeredLoad, trafficIntensity, beta, probDelay, meanWaitDelayed, meanWait, meanQueueLength, meanNumber and utilization. An overloaded model (beta <= 0) has no QED limit: probDelay is 1 and the waiting-time fields are infinite.

Return type:

References

S. Halfin, W. Whitt (1981). Heavy-traffic limits for queues with many exponential servers. Operations Research 29(3), 567-588.

qsys_mmk_qed_alpha(beta)[source]

The Halfin-Whitt delay-probability function alpha(beta) = [1 + beta*Phi(beta)/phi(beta)]^-1 for beta > 0, with phi and Phi the standard normal density and cdf.

It is the limit of the Erlang C delay probability of the M/M/s queue as s -> inf with beta = (1-rho)sqrt(s) held fixed, decreasing strictly from 1 at beta = 0 to 0 as beta -> inf, which is what makes it invertible for staffing. Non-positive beta returns 1: with no server slack every arrival is delayed.

Evaluated as phi/(phi + beta*Phi) rather than as the reciprocal of 1 + beta*Phi/phi: the two are the same function, but the quotient Phi/phi overflows once phi underflows (beta beyond about 38), whereas this form degrades to 0/(0+beta) = 0, the correct limit.

Parameters:: beta (float | ndarray) – the QED server-slack parameter, scalar or array
Returns:: alpha(beta), of the same shape as the input.
Return type:: float | ndarray

References

S. Halfin, W. Whitt (1981). Heavy-traffic limits for queues with many exponential servers. Operations Research 29(3), 567-588.

qsys_mmk_qed_staffing(lambda_val, mu, target, criterion='delay', exact=False, maxServers=10**7)[source]

Square-root staffing of the M/M/s queue.

Invert alpha(beta) = target for the server slack and staff s = ceil(a + beta*sqrt(a)) with a = lambda/mu: the base a erlangs of work plus a cushion that grows only as the square root of the load. Doubling the load needs only sqrt(2) times the cushion, which is why large service systems can be both highly utilized and responsive.

Parameters:

lambda_val (float) – arrival rate
mu (float) – service rate of one server
target (Any) – the target, read according to criterion: a probability for 'delay', a time for 'meanwait', or a dict with keys deadline and level for 'servicelevel'
criterion (str) – 'delay' (P(W>0) <= target), 'meanwait' (E[W] <= target) or 'servicelevel' (P(W <= deadline) >= level)
exact (bool) – walk s until the EXACT Erlang C measure meets the target, starting from the square-root answer
maxServers (int) – cap on that walk

Returns:

Dict with numServers, beta, betaTarget, offeredLoad, probDelay, meanWait, exactUsed and, for the service-level criterion, serviceLevel.

Return type:

References

S. Halfin, W. Whitt (1981). Heavy-traffic limits for queues with many exponential servers. Operations Research 29(3), 567-588. The staffing form is the standard reading of that limit; see also W. Whitt (2007), Naval Research Logistics 54(5), 476-484.

qsys_mtginf(lambdaFun, serviceCcdf, ES, tvals, startTime=-np.inf, ES2=None, servicePdf=None, tol=1e-12, panels=4000, maxAge=1e12)[source]

Exact time-varying analysis of the Mt/G/infinity queue.

With a non-homogeneous Poisson arrival rate lambda(t) and iid service times S, the number in system at time t is POISSON with mean

\[m(t) = E\left[\int_{t-S}^{t}\lambda(u)du\right] = E[S]\,E[\lambda(t-S_e)] = \int_0^\infty \lambda(t-x)P(S>x)dx\]

where S_e is the stationary-excess (equilibrium) law of S, with density P(S>x)/E[S]. This is exact, not an approximation: infinitely many servers mean customers never interact, so the model is a Poisson random measure and the whole distribution is known.

THE PHYSICS. Writing the mean as E[S] E[lambda(t - S_e)] says the time-varying load is the stationary load E[S]lambda(t) subjected to a TIME LAG and a SPACE SHIFT: to first order m(t) ~ E[S] lambda(t - E[S_e]) with E[S_e] = E[S^2]/(2E[S]), so peak congestion lags peak arrival rate, and by more than the mean service time when the service law is variable. The pointwise stationary approximation E[S]lambda(t) is the zeroth-order term, which is why it misses the lag.

Parameters:

lambdaFun (Callable[[Any], Any]) – the arrival rate, ideally array-aware; must accept arguments in the past when startTime is infinite
serviceCcdf (Callable[[Any], Any]) – G^c(x) = P(S > x)
ES (float) – the mean service time
tvals (Sequence[float]) – the times at which to evaluate
startTime (float) – time the system started empty; the default -inf assumes the arrival rate has been running forever
ES2 (float | None) – the second moment of the service time, for the lag approximation
servicePdf (Callable[[Any], Any] | None) – the service density, used for the exact departure rate; when absent the departure rate comes from the flow balance m'(t) = lambda(t) - delta(t) by a central difference
tol (float) – service-tail cut for the age integral
panels (int) – Simpson panels for that integral
maxAge (float) – cap on the age integrated over

Returns:

Dict with times, meanNumber (the Poisson mean m(t)), varNumber (equal to it), departureRate, arrivalRate, offeredLoadPSA (the pointwise stationary approximation E[S]lambda(t)) and, when ES2 is given, meanLag (E[S_e]) and lagApproximation (E[S]lambda(t-E[S_e])).

Return type:

References

S. G. Eick, W. A. Massey, W. Whitt (1993). The physics of the Mt/G/inf queue. Operations Research 41(4), 731-742.

qsys_gig1_bnds_extremal(lambda_val, mu, ca, cs, K=4000, N=2000, skipTight=False)[source]

Extremal two-moment bounds for the GI/GI/1 queue.

Two moments do not determine E[W]; they determine a SET of possible values, and the width of that set is the honest uncertainty in any two-moment approximation. The extremal laws attain its ends: the lower end with deterministic interarrival times and a three-point service law on multiples of that interval (closed form, eq. 2.12), the upper end asymptotically with TWO-POINT laws, an interarrival law with an atom at 0 and a service law whose upper atom runs to infinity as its probability vanishes. Making an interarrival time larger only empties the queue once; making a service time larger delays everyone behind it, which is why the two ends look so different.

The upper end is reduced to a D(1/p)/RS(D(rho),p)/1 model with p = 1/(1+ca^2) and evaluated by Spitzer’s identity with the negative binomial pmf, so it is a truncated numerical limit rather than a formula; the closed-form companion (eq. 3.4) is within about 1% of it.

Parameters:

lambda_val (float) – arrival rate
mu (float) – service rate
ca (float) – coefficient of variation of the interarrival time
cs (float) – coefficient of variation of the service time
K (int) – truncation of the negative binomial value
N (int) – truncation of the random-walk length
skipTight (bool) – skip the O(K*N) tight bound and return the closed forms only

Returns:

trafficIntensity, lowerBound, upperBound, upperBoundClosed, upperBoundDaley, upperBoundKingman, heavyTraffic, delta, relativeWidth, tightComputed.

Return type:

Dict of TIMES IN QUEUE (add 1/mu for response times)

References

Y. Chen, W. Whitt (2020). Algorithms for the upper bound mean waiting time in the GI/GI/1 queue. Queueing Systems 94, 327-356.

qsys_gtmtst_fluid(lambdaFun, sFun, muFun, patienceCcdf, T, dt=None, B0=0.0, w0=0.0, sPrimeFun=None, patiencePdf=None, lambdaPast=None)[source]

The Gt/Mt/st+GI many-server fluid queue.

Time-varying arrival rate lambda(t), time-varying staffing s(t), exponential service at the time-varying rate mu(t), general patience with complementary cdf F^c, and unlimited waiting room.

THE MODEL ALTERNATES BETWEEN TWO REGIMES and the whole algorithm is the bookkeeping of that alternation:

UNDERLOADED: the queue is empty and every arrival enters service at once, so the system is the infinite-server fluid model and B obeys B'(t) = lambda(t) - mu(t)B(t) (eq. 18 of the reference, in its Mt form). It ends when B reaches s while lambda exceeds the rate Gamma(t) = s'(t) + s(t)mu(t) at which capacity frees up (eq. 15).
OVERLOADED: every server is busy, B(t) = s(t), fluid enters service at exactly Gamma(t), and the queue is described by its BOUNDARY WAITING TIME w(t), the age of the oldest fluid still waiting. Content of age x is what arrived x ago and has not yet abandoned, q(t,x) = lambda(t-x)F^c(x), and the boundary moves by the delay differential equation (eq. 21)

w’(t) = 1 - Gamma(t) / [lambda(t-w(t)) F^c(w(t))].

It ends when w returns to 0 with lambda no longer above Gamma (eq. 14).

WHY w AND NOT Q. The queue content is a functional of w, but not the other way round: two systems with the same Q and different age profiles abandon at different rates. Tracking the boundary keeps the age profile exact, which is what makes a general patience law admissible at all.

Parameters:

lambdaFun – arrival rate lambda(t)
sFun – staffing s(t), a positive function or a constant
muFun – service rate mu(t), a function or a constant
patienceCcdf (Callable[[Any], Any]) – F^c(x) = P(patience > x)
T (float) – horizon; the model is solved on [0,T]
dt (float) – grid step, default T/2000
B0 (float) – fluid in service at time 0
w0 (float) – boundary waiting time at time 0, 0 for an empty queue
sPrimeFun – s’(t); differentiated numerically from sFun when absent
patiencePdf (Callable[[Any], Any] | None) – the patience density, for the abandonment rate; differenced from the ccdf when absent
lambdaPast – the arrival rate before time 0, needed only when the queue starts non-empty; defaults to lambdaFun evaluated at negative times

Returns:

times, regime (1 overloaded, 0 underloaded), B (fluid in service), Q (fluid in queue), X = B+Q, w (boundary waiting time), v (potential waiting time), sigma (service completion rate), alpha (abandonment rate), utilization (B/s), arrivalRate, staffing, capacityRate (Gamma).

Return type:

Dict on the grid

References

Y. Liu, W. Whitt (2012). The Gt/GI/st+GI many-server fluid queue. Queueing Systems 71, 405-444; Y. Liu, W. Whitt (2014). Algorithms for time-varying networks of many-server fluid queues. INFORMS Journal on Computing 26(1), 59-73.

qsys_ggnm_diffusion(lambda_val, mu, n, m, ca, cs, serviceCcdf=None, tol=1e-12, panels=4000)[source]

Diffusion approximation for the G/GI/n/m queue.

A general arrival process characterized by its rate and its variability parameter ca^2, iid general service times of mean 1/mu and SCV cs^2, n servers and m extra waiting spaces.

THE APPROXIMATION IS ONE DIFFUSION WITH TWO REGIONS. Below the staffing level the queue behaves like an infinite-server system, whose limit is NORMAL with variance-to-mean ratio the ASYMPTOTIC PEAKEDNESS

z = 1 + (ca^2 - 1) omega_G, omega_G = int G^c(x)^2 dx / int G^c(x) dx

(eqs. 1.6-1.7); above it the queue behaves like a single-server queue, whose limit is EXPONENTIAL with variability v = (ca^2 + cs^2)/2 (eq. 3.7). The steady-state law is therefore a normal piece spliced to an exponential piece, and every measure below is an integral of that density (eq. 3.14).

WHAT z SAYS. The service-time distribution enters the delay probability ONLY through omega_G, which is 1 for deterministic service, 1/2 for exponential, and falls toward 0 as service gets more variable. So when ca^2 = 1 the delay probability does not depend on the service law at all (z = 1), which is the long-standing M/GI/n approximation by M/M/n; away from ca^2 = 1 it does, and this quantifies how much.

The delay probability is alpha(beta/sqrt(z)) with the Halfin-Whitt function alpha when m is infinite (eq. 3.10), so this generalizes qsys_mmk_qed().

Parameters:

lambda_val (float) – arrival rate
mu (float) – service rate of one server
n (int) – number of servers
m (float) – extra waiting spaces; float('inf') for an unbounded queue
ca (float) – coefficient of variation of the interarrival time
cs (float) – coefficient of variation of the service time
serviceCcdf (Callable[[float], float] | None) – G^c(x) = P(S > x); the exponential of rate mu by default
tol (float) – service-tail cut for the peakedness integral
panels (int) – Simpson panels for it

Returns:

Dict with beta (the QED server slack), gamma (the scaled waiting room), peakedness (z), variability (v), probDelay, probBlock, meanQueueLength (customers waiting), meanNumber (in system), meanWait, utilization and throughput.

Return type:

References

W. Whitt (2004). A diffusion approximation for the G/GI/n/m queue. Operations Research 52(6), 922-941.

qsys_ggingi_tga(lambda_val, mu, n, ca, cs, patienceCcdf, patiencePdf=None, serviceCcdf=None)[source]

Truncated Gaussian approximation (TGA-G) for the G/GI/n+GI queue.

A general stationary arrival process of rate lambda_val and variability ca^2, iid general service of mean 1/mu and variability cs^2, n servers, unlimited waiting room, and iid general patience.

THE APPROXIMATION IS A FLUID CENTRE PLUS A GAUSSIAN FLUCTUATION, TRUNCATED. In the efficiency-driven regime (rho > 1 held fixed as n grows) the fluid limit gives the centre – all servers busy, waiting time w = F^-1(1-1/rho), queue Q = lambda int_0^w F^c – and the many-server central limit theorem gives a NORMAL fluctuation of order sqrt(n) around it. Adding the two directly can produce negative queues and negative waits, so both are TRUNCATED at zero, which is what makes the formulas usable down to moderate overload; the paper reports good accuracy for rho > 1.02 and abandonment rates below 2.

Three independent sources of variability enter separately, which is what lets the exponential-service formula be generalized: the service law appears only as the factor (cs+1)rho in sigma_W^2 (eq. 24), reducing to the exponential case at cs = 1.

An UNDERLOADED model (rho <= 1) has no queue in the limit; the number in system is then normal with the infinite-server variance, whose variance-to-mean ratio is the asymptotic peakedness of qsys_ggnm_diffusion().

Parameters:

lambda_val (float) – arrival rate
mu (float) – service rate of one server
n (int) – number of servers
ca (float) – coefficient of variation of the interarrival time
cs (float) – coefficient of variation of the service time
patienceCcdf (Callable[[float], float]) – F^c(x) = P(patience > x)
patiencePdf (Callable[[float], float] | None) – the patience density; differenced from the ccdf when absent
serviceCcdf (Callable[[float], float] | None) – G^c(x) = P(S > x), used only in the underloaded branch

Returns:

Dict with regime, trafficIntensity, fluidWait, fluidQueueLength, meanWait, varWait, meanQueueLength, varQueueLength, meanNumberInService, meanNumber, probDelay, probAbandon, sigmaW, sigmaX.

Return type:

References

Y. Liu, W. Whitt, Y. Yu (2016). Approximations for heavily-loaded G/GI/n+GI queues. Naval Research Logistics 63(3), 187-217.

qsys_mtgs0_mol(lambdaFun, serviceCcdf, ES, s, tvals, startTime=-np.inf, delay=False, ES2=None, **kwargs)[source]

Modified-offered-load (MOL) and pointwise-stationary (PSA) approximations for a time-varying multiserver system.

THE ONE IDEA. A stationary loss system with offered load a blocks with probability B(s,a). In a time-varying system the question is WHICH LOAD to put in that formula. PSA uses the instantaneous one, lambda(t)E[S]. MOL uses the offered load of the corresponding INFINITE-SERVER system,

m(t) = E[S] E[lambda(t - S_e)] = int_0^inf lambda(t-x)P(S>x)dx,

which is exact for that system and therefore carries the TIME LAG and the smoothing that the finite-server system also has. MOL is then B(s, m(t)). The difference between the two is precisely the lag: PSA peaks when the arrival rate peaks, MOL peaks later, and the real system peaks later too.

WHY IT WORKS. The blocking system differs from the infinite-server one only in what happens at the ceiling, and the ceiling does not change the AGE structure of the load much when blocking is not extreme. That is why the approximation is asymptotically correct in the many-server regime and degrades when blocking is heavy.

WHAT TO EXPECT. Measured against the exact time-varying birth-death chain on a sinusoidal rate, MOL cuts the mean RELATIVE error roughly threefold (0.13 against 0.44 at s = 100), because it gets the phase right. It does not always win on ABSOLUTE error: that is dominated by the peak of the cycle, where both approximations are weakest. Under constant input MOL is exact, reducing to the stationary Erlang formula.

Parameters:

lambdaFun (Callable[[Any], Any]) – the arrival rate; must accept arguments in the past when startTime is infinite
serviceCcdf (Callable[[Any], Any]) – G^c(x) = P(S > x)
ES (float) – the mean service time
s (int) – number of servers
tvals (Sequence[float]) – times at which to evaluate
startTime (float) – time the system started empty; -inf assumes an infinite past
delay (bool) – use Erlang C rather than Erlang B, i.e. approximate the DELAY probability of an Mt/M/s queue rather than the blocking probability of an Mt/G/s/0 loss system
ES2 (float | None) – second moment of the service time, passed through for the time lag
**kwargs – passed to qsys_mtginf()

Returns:

Dict with times, offeredLoad (m(t)), instantLoad (lambda(t)E[S]), probBlockMOL, probBlockPSA, meanBusyMOL (the carried load m(t)(1-B) for the loss model), and, when ES2 is given, meanLag.

Return type:

References

W. A. Massey, W. Whitt (1994). An analysis of the modified offered load approximation for the nonstationary Erlang loss model. Annals of Applied Probability 4(4), 1145-1160; W. Whitt (1991). The pointwise stationary approximation for Mt/Mt/s queues is asymptotically correct as the rates increase. Management Science 37(3), 307-314.

erlang_b(s, a)[source]

Erlang B blocking probability with s servers and offered load a, by the recursion B_j = a B_{j-1}/(j + a B_{j-1}), which never forms a^s/s! and so never overflows.

Parameters:

s (int) – number of servers
a (float) – offered load in erlangs

Returns:

The probability that all servers are busy.

Return type:

float

erlang_c(s, a)[source]

Erlang C delay probability with s servers and offered load a, from the same recursion; 1 when the load saturates the servers.

Parameters:

s (int) – number of servers
a (float) – offered load in erlangs

Returns:

The probability that an arrival waits.

Return type:

float

qsys_maxima_twomoment(n, mean, cs2, q=None, exactFitted=True)[source]

Approximate the maximum of n iid non-negative variables from two moments.

THE SHAPE OF THE ANSWER. For a law with an exponential-like tail the maximum of n samples grows like c~^2 (log n + ...): doubling n adds a constant, it does not scale the answer. What the two moments buy is the SLOPE c~^2 of that logarithm and an offset eta:

x_n(q) = c~^2 [ log(n eta) - log log(1/q) ], E[M_n] = c~^2 [log(n eta) + gamma]

with, for cs2 >= 1, c~^2 = cs2 and eta = (cs2+1)/(2 cs2^2) from the H2 representative, and for cs2 < 1 the shifted-exponential representative c~^2 = sqrt(cs2), eta = exp((1-sqrt(cs2))/sqrt(cs2)).

WHEN NOT TO USE IT. The extreme-value form needs n past a threshold n* ~ cs2/q, because with a highly variable law most of the n samples come from the short component and only about n p of them can contend for the maximum. Measured against exact maxima, the closed form is within a few percent for n >= 100 at cs2 = 4 and 16, and useless at n = 10 for cs2 = 16 – which is exactly what n* predicts.

AND WHEN TWO MOMENTS ARE NOT ENOUGH. Below cs2 = 1 the maximum is genuinely family-dependent: an Erlang and a shifted exponential with the same two moments have maxima that differ by tens of percent and diverge as n grows, because their tails decay at different rates. The paper’s own caution. exactFitted therefore also returns the maximum computed exactly from the fitted representative, which is the reliable route it recommends.

Parameters:

n (int) – the number of samples
mean (float) – the mean of the underlying law
cs2 (float) – its squared coefficient of variation
q (float | None) – a quantile level in (0,1); the mean is returned when absent
exactFitted (bool) – also compute the maximum exactly from the fitted representative distribution, by integrating 1-F^n

Returns:

Dict with value (the closed-form mean or quantile), slope (c~^2 mean), eta, threshold (n*), reliable (whether n >= n*), family, and, when requested, exactFittedValue.

Return type: