Non-Product-Form Networks

Approximations for intractable queueing networks.

The npfqn module provides approximation methods for non-product-form queueing networks that do not satisfy BCMP conditions.

Key function categories:

Traffic merging: npfqn_traffic_merge(), npfqn_traffic_merge_cs()
Traffic splitting: npfqn_traffic_split_cs()
Service approximations: npfqn_nonexp_approx()

Non-Product-Form Queueing Network (NPFQN) algorithms.

Native Python implementations for approximating performance of non-product-form queueing networks.

Key algorithms:: npfqn_nonexp_approx: Non-exponential distribution approximation npfqn_traffic_merge: Merge multiple MMAP traffic flows npfqn_traffic_merge_cs: Merge traffic flows with class switching npfqn_traffic_split_cs: Split traffic flows with class switching

npfqn_nonexp_approx(method, sn, ST, V, SCV, Tin, Uin, gamma, nservers)[source]

Approximates non-product-form queueing networks using the specified method.

This function adjusts service times and other parameters to account for non-exponential service time distributions in product-form analysis.

Parameters:

method (str) – Approximation method (“default”, “none”, “hmva”, “interp”)
sn (NetworkStruct) – Network structure
ST (ndarray) – Service time matrix (M x K)
V (ndarray | None) – Visit ratios matrix (M x K), optional
SCV (ndarray) – Squared coefficient of variation matrix (M x K)
Tin (ndarray) – Initial throughput matrix (M x K)
Uin (ndarray) – Initial utilization matrix (M x K)
gamma (ndarray) – Gamma correction matrix (M x 1)
nservers (ndarray) – Number of servers matrix (M x 1)

Returns:

NpfqnNonexpApproxResult with updated matrices

Raises:

ValueError – If unknown approximation method is specified

Return type:

NpfqnNonexpApproxResult

class NpfqnNonexpApproxResult(ST, gamma, nservers, rho, scva, scvs, eta)[source]

Bases: object

Result of non-exponential approximation.

ST: ndarray

gamma: ndarray

nservers: ndarray

rho: ndarray

scva: ndarray

scvs: ndarray

eta: ndarray

npfqn_traffic_merge(MMAPa, config_merge='default', config_compress=None)[source]

Merge multiple MMAP traffic flows.

Combines multiple MMAPs using specified aggregation strategies. Supports various merge configurations for different network topologies.

Parameters:

MMAPa (Dict[int, List[ndarray] | None]) – Dictionary of MMAP traffic flows to be merged. Keys are integer indices, values are MMAP lists [D0, D1, …].
config_merge (str) – Merge configuration. Options: - “default”, “super”: Use MMAP superposition - “mixture”: Apply mixture fitting after superposition - “interpos”: Interposition method (falls back to super)
config_compress (str | None) – Compression configuration. Options: - None, “none”: No compression - “default”: Apply default compression

Returns:

Merged and normalized MMAP, or None if input is empty.

Raises:

RuntimeError – If unsupported merge configuration is provided.

Return type:

List[ndarray] | None

Example

>>> mmap1 = [np.array([[-1.0]]), np.array([[1.0]])]
>>> mmap2 = [np.array([[-2.0]]), np.array([[2.0]])]
>>> result = npfqn_traffic_merge({0: mmap1, 1: mmap2})

npfqn_traffic_merge_cs(MMAPs, prob, config='default')[source]

Merge MMAP traffic flows with class switching.

Combines multiple MMAPs while applying class switching transformations based on the probability matrix.

Parameters:

MMAPs (Dict[int, List[ndarray]]) – Dictionary of MMAP traffic flows indexed by source.
prob (ndarray) –
Class switching probability matrix ((n*R) x R) where:
- n is the number of sources
- R is the number of classes
- prob[(i-1)*R + r, s] = probability that class r from source i becomes class s in the merged stream
config (str) – Merge configuration (“default” or “super”).

Returns:

Merged MMAP with class switching applied, or None if empty.

Return type:

List[ndarray] | None

Algorithm:

Apply mmap_mark to each source MMAP to encode class switching
Superpose all marked MMAPs
Return result

Example

>>> mmap1 = [np.array([[-1.0]]), np.array([[0.5]]), np.array([[0.5]])]
>>> mmap2 = [np.array([[-2.0]]), np.array([[1.0]]), np.array([[1.0]])]
>>> prob = np.array([[0.8, 0.2], [0.3, 0.7], [0.6, 0.4], [0.5, 0.5]])
>>> result = npfqn_traffic_merge_cs({0: mmap1, 1: mmap2}, prob)

npfqn_traffic_split_cs(MMAP_input, P)[source]

Split MMAP traffic flows with class switching.

Decomposes a single MMAP into multiple MMAPs based on routing probabilities and class switching matrix.

Parameters:

MMAP_input (List[ndarray]) – Input MMAP as list [D0, D1, D2, …].
P (ndarray) –
Class switching probability matrix (R x J) where:
- R is the number of arrival classes
- J = M * R where M is the number of destinations
- P[r, (jst-1)*R + s] = probability that class r arrival becomes class s at destination jst

Returns:

Dictionary mapping destination index to split MMAP. Keys are 0-indexed destination indices.

Return type:

Dict[int, List[ndarray] | None]

Algorithm:

Parse dimensions from P matrix
For each destination, create MMAP with weighted marking matrices
Normalize each result

Example

>>> mmap = [np.array([[-3.0]]), np.array([[1.0]]), np.array([[2.0]])]
>>> P = np.array([[0.8, 0.2, 0.5, 0.5], [0.3, 0.7, 0.6, 0.4]])
>>> result = npfqn_traffic_split_cs(mmap, P)

npfqn_sqd(sn, N=None, calibration_mode=0, server_blocking_time=True, neighbor_mode='downstream', v1_policy='compound', initial_v1=None)[source]

Solve a BAS closed network from its NetworkStruct.

Parameters:

sn (NetworkStruct) – chain-aggregated network structure
N (int | None) – total closed-class population (defaults to sn.nclosedjobs)
calibration_mode (int) – 0=base, 1=fixed heuristic, 2=blocking-aware
server_blocking_time (bool) – add a manufacturing-blocking term to server time
neighbor_mode (str) – ‘downstream’ (routed) or ‘ownserver’ blocking aggregation
v1_policy (str) – ‘compound’ or ‘fresh’ load-dependent rate-scale update
initial_v1 (ndarray | None) – optional per-station initial V1 (None = _INITIAL_V1)

Returns:

NpfqnSqdResult with per-station X, Q, U, R.

Return type:

NpfqnSqdResult

class NpfqnSqdResult(X, Q, U, R)[source]

Bases: object

Per-station results of the BAS approximation (single chain).

Variables:

X (numpy.ndarray) – (M,) per-station throughput
Q (numpy.ndarray) – (M,) per-station queue length
U (numpy.ndarray) – (M,) per-station utilization
R (numpy.ndarray) – (M,) per-station residence time

X: ndarray

Q: ndarray

U: ndarray

R: ndarray

npfqn_rqna_weight(t)[source]

Canonical RBM correlation weight function w*(t) used by the RQNA.

w*(t) = 1 - (1 - c*(t))/(2 t), where c*(t) is the correlation function of the stationary version of canonical reflected Brownian motion (drift -1, diffusion coefficient 1),

c*(t) = 2(1 - 2t - t^2) Phi^c(sqrt(t)) + 2 sqrt(t) phi(sqrt(t)) (1 + t),

with Phi^c the standard-normal complementary cdf and phi its density. The weight is monotonically increasing with w*(0)=0 and w*(Inf)=1.

Reference: Whitt and You (2018), eqs. (24)-(25).

Parameters:: t – scalar or array of nonnegative time arguments
Returns:: Weight(s) w*(t) in [0,1], same shape as t (scalar in -> float out).

npfqn_traffic_rqt(lambda0, Gamma0, alpha0, F)[source]

Effective arrival process perceived at each node of a single-class open queueing network under the Robust Queueing Theory calculus.

The characterization composes three operators: passage through a queue with adversarial servers leaves the uncertainty set unchanged (robust Burke, Theorem 4), superposition merges sets by Theorem 5, and thinning by a fraction f scales the rate by f and the variability by f^(-1/alpha) (Theorem 6). The resulting equations are:

lambda_j = lambda0_j + sum_i lambda_i f_ij
Gamma_j  = (1/lambda_j) [ 1{a0_j=ab_j} (lambda0_j Gamma0_j)^(p_j)
            + sum_i 1{ab_i=ab_j} (lambda_i Gamma_i)^(p_i) f_ij ]^(1/p_j)

with p_j = ab_j/(ab_j-1) and ab_j the minimum tail coefficient among the streams feeding j: the heaviest tail upstream dominates. Both are solved exactly rather than iteratively. The rate equations are the usual traffic equations, and in the variables z_j = (lambda_j Gamma_j)^(p_j) the variability equations are linear as well, so each is one linear system; ab is obtained by propagating the minimum to a fixed point.

Parameters:

lambda0 – (J,) external arrival rate at each node, 0 where there is none
Gamma0 – (J,) variability parameter of each external arrival process, which for a renewal stream is the interarrival standard deviation
alpha0 – (J,) tail coefficient in (1,2] of each external arrival process
F – (J,J) routing probability matrix, F[i,j] = fraction of the jobs leaving node i that go to node j (row sums <= 1)

Returns:

Tuple of (lambda, Gamma, alpha), each (J,)

Return type:

Tuple[ndarray, ndarray, ndarray]

npfqn_bnd_bpt(lambda0, mu, P, station_of, c=None)[source]

Lower bound on sum_r c[r]*x[r] over the achievable region.

x[r] is the mean sojourn time of class r; the bound is valid for EVERY non-idling scheduling policy.

A “class” here is a buffer with its own exponential service rate and its own Markovian routing, so a station serving several customer types owns one class per type. The network is open: class r receives external Poisson arrivals at rate lambda0[r] and, on completing service, becomes class r' with probability P[r, r'] or leaves with the row deficit.

METHOD. Uniformize the chain and let R(t) = sum_r f(r) n_r(t) for an arbitrary vector f. The steady-state balance of E[R^2] is an identity quadratic in f; since it holds for every f, the two sides’ coefficient matrices agree entrywise. Diagonal entries give one equation per class, off-diagonal entries one per unordered pair, in the variables x_r, I(r,l) = E[1{sigma(r) busy with r} n_l] and N(i,l) = E[1{station i idle} n_l]. A third block states that the events “station i serves class r” and “station i idle” are mutually exclusive and exhaustive, so their terms sum to E[n_l] = lambda_l x_l. Minimizing over this polyhedron is a relaxation, hence a lower bound.

EXACT ON M/M/1: the LP reduces to mu*I11 - lambda^2*x = lambda and I11 + N11 = lambda*x with N11 >= 0, whence x >= 1/(mu-lambda) with equality.

NOT INCLUDED, DELIBERATELY. The valid inequality I(r,r) >= rho_r would tighten the relaxation but is not part of the reference’s characterization, and reproducing the reference’s published bounds is the acceptance test.

Parameters:

lambda0 (Sequence[float]) – external Poisson arrival rate into each class (0 if none).
mu (Sequence[float]) – exponential service rate of each class.
P – K x K routing, P[r, r'] = P(class r becomes r’ after service); row sums must not exceed 1.
station_of (Sequence[int]) – zero-based station index of each class.
c (Sequence[float] | None) – objective weights; None means all ones.

Returns:

NpfqnBndBptResult.

Return type:

NpfqnBndBptResult

Reference:: D. Bertsimas, I. Paschalidis, J. Tsitsiklis (1994). Optimization of multiclass queueing networks: polyhedral and nonlinear characterizations of achievable performance. Annals of Applied Probability 4(1), 43-75. See also D. Bertsimas (1995), Queueing Systems 21, 337-389, Theorem 9, which restates the same characterization.

class NpfqnBndBptResult(zlb, x, lambda_, rho, rho_station, nvars, nrows)[source]

Bases: object

Outcome of one achievable-region LP solve.

Variables:

zlb – lower bound on sum_r c[r]*x[r].
x – the x block of the LP optimizer, length K. Only the objective value is a bound; an individual x[r] is a vertex coordinate, not a bound on class r, unless c is the r-th unit vector.
lambda – effective arrival rate of each class.
rho – per-class utilization lambda_r/mu_r.
rho_station – per-station utilization.
nvars – number of LP variables.
nrows – number of LP rows.

zlb

x

lambda_

rho

rho_station

nvars

nrows

npfqn_bnd_bgt(lambda_, mu, sigma, J=None)[source]

Upper bound the steady-state queue lengths of a multitype network.

The bound is valid for EVERY work-conserving Markovian policy.

MODEL. J single-server stations; I customer types; type i arrives as a Poisson stream of rate lambda_[i] and passes through stages k = 0..len(mu[i])-1, stage k being served at station sigma[i][k] at exponential rate mu[i][k]. Class (i,k) is the buffer of type i at stage k; N = sum_i len(mu[i]) is the number of classes.

METHOD. Solve the Down-Meyn global-stability linear program GLP[dm], eq. (25)-(28) of the reference, in the piecewise-linear Lyapunov function Phi(x) = max_j L^j'x:

L^j(i,1) lambda_i + mu(i,k) (L^j(i,k+1) - L^j(i,k)) + V_j <= -gamma
                                             for (i,k) in station j
mu(i,k) (L^j(i,k+1) - L^j(i,k)) <= V_j       for (i,k) not in j
(1/(J-1)) sum_{j' != j} L^j'(i,k) >= L^j(i,k)  for (i,k) not in j
L, V, gamma >= 0

with L^j(i,Ji+1) = 0. A feasible solution with gamma > 0 certifies that EVERY work-conserving policy is stable, and a smoothed Phi is then a Lyapunov function with drift gamma/4 and an explicit exception parameter, giving the reference’s Theorem 4 bound:

E[L^j'Q] <= 16 N J^2 (J-1) (Lmax+gamma)^3 / gamma^2
            + 8 (Lmax + gamma/2)^2 / gamma  =: U

for every j, whence E[Q(i,k)] <= U / max_j L^j(i,k).

THE RATES ARE RESCALED so that sum_i lambda_i + sum_{i,k} mu(i,k) = 1, the uniformization the reference imposes before Theorem 4. Queue lengths are counts and are unaffected by the time scale.

NORMALIZATION, WHICH THE REFERENCE LEAVES OPEN. GLP[dm] is homogeneous and so is the bound, so this routine fixes L^j(i,k) <= 1 and MAXIMIZES gamma, then breaks ties among gamma-optimal solutions by maximizing sum L: a degenerate optimum can otherwise zero some L^j(i,k) and report an infinite bound for a class for no reason.

THE BOUND IS LOOSE, and knowingly so: the exception parameter carries (Lmax+gamma)^3/gamma^2 and dominates as soon as J > 1. What is sharp is the STABILITY CERTIFICATE gamma > 0 and the geometric tail RATE.

Reference:: D. Bertsimas, D. Gamarnik, J. N. Tsitsiklis (2001). Performance of multiclass Markovian queueing networks via piecewise linear Lyapunov functions. Annals of Applied Probability 11(4), 1384-1428, Section 5.1 (GLP[dm] of Down and Meyn 1997, and Theorem 4).

class NpfqnBndBgtResult(**kw)[source]

Bases: object

Outcome of one GLP[dm] solve and of the Theorem 4 bound.

Variables:

Qub – list of arrays, Qub[i][k] upper bounds E[Q(i,k)]; inf where the LP optimum leaves max_j L^j(i,k) = 0.
gamma – the drift certificate; strictly positive on success.
Lmax – max over j and (i,k) of L.
L – the Lyapunov coefficients, shape (J, N).
V – the per-station slack V_j.
B – the exception parameter of the smoothed Lyapunov function.
U – the Theorem 4 bound on E[L^j'Q], the same for every j.
tail_ratio – geometric decay ratio of the tail bound.
tail_step – step of the tail bound, 2(Lmax+gamma/2).
rho – per-class nominal load.
rho_station – per-station nominal load.
scale – the uniformization divisor applied to lambda and mu.
class_station (class_type, class_stage,) – per-class index maps.

Qub

gamma

Lmax

L

V

B

U

tail_ratio

tail_step

rho

rho_station

scale

class_type

class_stage

class_station

npfqn_traffic_idc(lambda0, P, c2a0, a0IdcFun, mu, cs2, sIdcFun, corrections=None)[source]

Traffic variability equations for the RQNA. Assembles and solves:

the limiting variability equations (eq. 42/44) for the asymptotic total-arrival variability parameters c2_{a,i} = I_{a,i}(Inf);
a solver (ctx.IaFun) for the time-dependent IDC equations (eq. 40/43), returning I_{a,i}(t) for all internal arrival flows, using the default correction terms alpha_{i,j} (eq. 34) and beta_i (eqs. 38-39) and tuning function h(rho)=rho^2.

Models a single-class open network of K single-server FCFS queues with Markovian routing P (P[i,j]=p_{i,j}).

Parameters:

lambda0 – (K,) external arrival rate into each queue
P – (K,K) routing matrix among queues
c2a0 – (K,) asymptotic IDC (SCV) of each external arrival process
a0IdcFun – callable a0IdcFun(t) -> (K,) external arrival IDC I_{a,0,i}(t)
mu – (K,) service rate at each queue
cs2 – (K,) service SCV c2_{s,i}
sIdcFun – callable sIdcFun(t) -> (K,) service IDC I_{s,i}(t)
corrections – optional dict with ‘alpha’/’beta’ bool toggles

Returns:

_TrafficIdcContext with lambda, rho, Xi, c2a, c2d, c2aij, c2x fields and method IaFun(t).

npfqn_traffic_split_rr(sn)[source]

Deterministic (round-robin) split degree of the departure stream of each station-class.

kRR[i, r] = k > 1 means that the class-r departures of station i are dispatched one-in-k by a round-robin node, so that a downstream flow carrying a fraction p of them is the k-fold convolution thinned with probability q = k*p and has SCV 1+p*(d2-k), against the Markovian 1+p*(d2-1). kRR[i, r] = 1 marks an ordinary probabilistic split.

Only two topologies admit the deterministic rule: the station dispatches round-robin itself, or it feeds with probability one a router that does and whose pointer no other flow advances. Anything else falls back to k=1.

Parameters:: sn – NetworkStruct
Returns:: (M, K) array of split degrees
Return type:: ndarray

npfqn_gtmtst_fluid(lambdaFuns, sFuns, muFuns, patienceCcdfs, P, T, dt=None, B0=None, w0=None, tol=1e-6, maxIter=100)[source]

A time-varying open network of many-server fluid queues with abandonment.

Each queue is the Gt/Mt/st+GI fluid queue of qsys_gtmtst_fluid(); the departure flow of queue i is routed to queue j with the (possibly time-varying) proportion P[i][j], whatever is left leaving the network.

THE NETWORK IS A FIXED POINT. The total arrival rate of queue j is lambda_j(t) = lambda_j^0(t) + sum_i sigma_i(t) P_ij(t) with sigma_i = mu_i B_i the service completion rate (eqs. 23-24), and sigma_i itself depends on lambda_i. The iteration starts from the external rates alone and adds one more traversal of the network per round, so the nth iterate is the fluid that has made n transitions; the map is a monotone contraction, so the rates increase to the fixed point.

Parameters:

lambdaFuns (Sequence) – external arrival rate of each queue
sFuns (Sequence) – staffing of each queue
muFuns (Sequence) – service rate of each queue
patienceCcdfs (Sequence) – patience ccdf of each queue
P – routing proportions, either an m x m array or a callable P(t) returning one
T (float) – horizon
dt (float) – grid step
B0 (Sequence[float] | None) – initial fluid in service at each queue
w0 (Sequence[float] | None) – initial boundary waiting time at each queue
tol (float) – sup-norm tolerance on the arrival-rate iteration
maxIter (int) – cap on the iterations

Returns:

Dict with times, queues (the per-queue dicts of qsys_gtmtst_fluid()), arrivalRates (the converged total rates, one row per queue), iterations and residual.

Return type:

Dict[str, Any]

References

Y. Liu, W. Whitt (2014). Algorithms for time-varying networks of many-server fluid queues. INFORMS Journal on Computing 26(1), 59-73.

npfqn_feedback_elim(P, rho, cs2=None, lambda_=None, immediateOnly=False)[source]

Eliminate near-immediate feedback from an open queueing network.

WHY FEEDBACK BREAKS DECOMPOSITION. A parametric decomposition treats the arrival stream at each station as if it were renewal. Feedback destroys that badly: a customer that leaves a busy station and comes straight back arrives exactly when the station is busy, so the flow is strongly correlated with the queue it feeds. The fix is not to model the correlation but to REMOVE the feedback, by folding the repeated visits into the service time.

THE TRANSFORMATION. With feedback probability p at a station, a customer is served a geometric number of times, so the effective service is S_p = sum_{i=1}^{N} S_i with N geometric of mean 1/(1-p). Hence

effective mean service E[S]/(1-p),
effective service SCV p + (1-p)cs^2 (eq. 37 and the line after it),
fresh arrival rate lambda(1-p),
per-visit waiting time = (1-p) times the wait in the modified system.

The modified system has the same heavy-traffic limits for queue length, workload, waiting time and external departures, so the elimination is asymptotically exact rather than merely plausible.

NEAR-IMMEDIATE, NOT JUST IMMEDIATE. Feedback rarely returns a customer in one hop. What matters is whether it returns WITHOUT PASSING A BUSIER STATION: a detour through a station of lower traffic intensity is fast on the time scale of the busy station, so it behaves like immediate feedback. The probability computed here is therefore the probability of returning to station i through stations of strictly smaller rho only.

Parameters:

P – routing matrix, substochastic, P[i][j] from station i to j
rho (Sequence[float]) – traffic intensity of each station, which fixes what counts as “near-immediate”
cs2 (Sequence[float] | None) – service SCV of each station; the modified SCVs are returned when given
lambda – arrival rate of each station; the modified rates are returned when given
immediateOnly (bool) – keep only the self-loops P[i][i], i.e. immediate feedback in the strict sense of Section 4.1

Returns:

Dict with feedbackProb (p-hat per station), modifiedScv, modifiedRates, modifiedRouting (the routing with the eliminated feedback removed and the remaining rows renormalized), and visitInflation (1/(1-p), the mean visits per customer).

Return type:

Dict[str, Any]

References

W. Whitt, W. You (2022). A robust queueing network analyzer based on indices of dispersion. Naval Research Logistics 69(1), 36-56, Section 4.

npfqn_dps_morrison(N, Z, S, w)[source]

Evaluate Morrison’s two-term approximation.

Parameters:

N (array_like) – Per-class populations, finite and positive.
Z (array_like) – Per-class mean think times, finite and positive.
S (array_like) – Per-class mean DPS service times, finite and positive.
w (array_like) – Per-class DPS weights, finite and positive.

Returns:

Mean queue lengths, sojourn times and throughputs, with the intermediate constants.

Return type:

NpfqnDpsMorrisonResult

class NpfqnDpsMorrisonResult(Q=None, R=None, X=None, Qlead=None, Rlead=None, sigma=None, W=None, rho=0.0, a=0.0, cB=0.0, cC=0.0, cD=0.0, cH=0.0, cI=0.0, cJ=0.0, cK=0.0, cL=0.0, cM=0.0, cQ=0.0, delta=0.0, cR=0.0, cS=0.0, cU=0.0, cA=0.0, cV=0.0)[source]