Product-Form Queueing Networks

MVA, convolution, and normalizing constant methods.

The pfqn module contains algorithms for product-form queueing networks, including Mean Value Analysis (MVA), convolution, and normalizing constant methods.

Key function categories:

Mean Value Analysis: pfqn_mva(), pfqn_mvald(), pfqn_mvams(), pfqn_mvamx()
Normalizing constants: pfqn_nc(), pfqn_ncld(), pfqn_fnc()
Convolution algorithms: pfqn_ca(), pfqn_bs(), pfqn_comom()
Linearizer methods: pfqn_linearizer(), pfqn_linearizerms(), pfqn_gflinearizer()
Bounds and approximations: pfqn_qzgblow(), pfqn_qzgbup(), pfqn_xzabalow(), pfqn_xzabaup()
LCFS queueing networks: pfqn_lcfsqn_mva(), pfqn_lcfsqn_nc(), pfqn_lcfsqn_ca()
Maximum entropy: me_oqn() - Maximum Entropy for open queueing networks

Product-form queueing network (PFQN) algorithms.

Native Python implementations of analytical algorithms for product-form queueing networks, including Mean Value Analysis (MVA), normalizing constant methods, and various approximation techniques.

Key algorithms:: pfqn_mva: Standard Mean Value Analysis pfqn_ca: Convolution Algorithm pfqn_nc: Normalizing Constant methods pfqn_bs: Balanced System analysis pfqn_aql: Approximate queue lengths

pfqn_mva(L, N, Z=None, mi=None)[source]

Mean Value Analysis for multi-class closed product-form network.

Implements the exact MVA algorithm using population recursion. Computes exact performance measures for closed product-form networks with load-independent stations.

Parameters:

L (numpy.ndarray) – Service demand matrix (M x R) where M is stations, R is classes
N (numpy.ndarray) – Population vector (1 x R or R,) - number of jobs per class
Z (numpy.ndarray) – Think time vector (1 x R or R,) - think time per class (default 0)
mi (numpy.ndarray) – Multiplicity vector (1 x M or M,) - servers per station (default 1)

Returns:

XN: Throughputs per class (1 x R)
CN: Response times per class (1 x R) - total cycle time
QN: Queue lengths (M x R)
UN: Utilizations (M x R)
RN: Residence times (M x R)
TN: Node throughputs (M x R)
AN: Arrival rates (M x R)

Return type:

Tuple of (XN, CN, QN, UN, RN, TN, AN) where

pfqn_mva_single_class(N, L, Z=0.0, mi=None)[source]

Mean Value Analysis for single-class closed network.

Simplified MVA for single customer class, with optional multi-server stations specified via mi (multiplicity).

Parameters:

N (int) – Number of customers
L (numpy.ndarray) – Service demands at each station (1D array of length M)
Z (float) – Think time (default 0)
mi (numpy.ndarray | None) – Number of servers at each station (default all 1)

Returns:

‘X’: Throughput
’Q’: Queue lengths (array of length M)
’R’: Residence times (array of length M)
’U’: Utilizations (array of length M)
’lG’: Log of normalizing constant

Return type:

dict with keys

pfqn_bs(L, N, Z=None, tol=1e-6, maxiter=1000, QN0=None, type_sched=None)[source]

Bard-Schweitzer Approximate Mean Value Analysis (MVA).

Iterative approximate MVA algorithm that uses the (N-1)/N correction for the arrival theorem, providing good accuracy for most networks.

Parameters:

L (numpy.ndarray) – Service demand matrix (M x R)
N (numpy.ndarray) – Population vector
Z (numpy.ndarray) – Think time vector (default 0)
tol (float) – Convergence tolerance (default 1e-6)
maxiter (int) – Maximum iterations (default 1000)
QN0 (numpy.ndarray) – Initial queue lengths (default: uniform distribution)
type_sched (numpy.ndarray) – Scheduling strategy per station (default: PS)

Returns:

XN: System throughputs (1 x R) QN: Mean queue lengths (M x R) UN: Utilizations (M x R) RN: Residence times (M x R) it: Number of iterations performed

Return type:

Tuple (XN, QN, UN, RN, it) matching MATLAB’s pfqn_bs

pfqn_aql(L, N, Z=None, max_iter=1000, tol=1e-6, QN0=None)[source]

Approximate Queue Length (AQL) algorithm.

Uses iterative approximation to compute queue lengths for large populations where exact MVA would be computationally expensive.

Parameters:

L (numpy.ndarray) – Service demand matrix (M x R)
N (numpy.ndarray) – Population vector
Z (numpy.ndarray) – Think time vector (default 0)
max_iter (int) – Maximum iterations (default 1000)
tol (float) – Convergence tolerance (default 1e-6)

Returns:

Same format as pfqn_mva

Return type:

Tuple[numpy.ndarray, numpy.ndarray, numpy.ndarray, numpy.ndarray, numpy.ndarray, numpy.ndarray, numpy.ndarray]

pfqn_sqni(L, N, Z=None)[source]

Square-root Non-iterative (SQNI) approximate MVA.

Implements a fast approximation for multi-class closed queueing networks that reduces the system to single-queue representations with interpolation-based corrections.

This method is particularly efficient for networks where one station dominates (bottleneck analysis), providing a good trade-off between accuracy and computational speed.

Parameters:

L (numpy.ndarray) – Service demand vector (1 x R or R,) - demands at the bottleneck queue
N (numpy.ndarray) – Population vector (1 x R or R,) - number of jobs per class
Z (numpy.ndarray) – Think time vector (1 x R or R,) - think time per class (default 0)

Returns:

Q: Queue lengths (2 x R) - first row for queue, second placeholder
U: Utilizations (2 x R) - first row for queue, second placeholder
X: Throughputs (1 x R)

Return type:

Tuple of (Q, U, X) where

Reference:: Based on the SQNI method for approximate MVA analysis.

pfqn_qd(L, N, ga=None, be=None, Q0=None, tol=1e-6, max_iter=1000)[source]

Queue-Dependent (QD) Approximate MVA.

Implements the QD-AMVA algorithm that uses queue-dependent correction factors to improve accuracy of approximate MVA for closed networks.

The algorithm iteratively computes queue lengths using: - A correction factor delta = (N_tot - 1) / N_tot - Per-class correction factor delta_r = (N_r - 1) / N_r - Optional scaling functions ga(A) and be(A) for advanced corrections

Parameters:

L (numpy.ndarray) – Service demand matrix (M x R) - rows are stations, columns are classes
N (numpy.ndarray) – Population vector (R,) - number of jobs per class
ga (callable) – Gamma scaling function ga(A) -> array(M,) (default: ones) A is the arrival queue seen by class, A[k] = 1 + delta * sum(Q[k,:])
be (callable) – Beta scaling function be(A) -> array(M, R) (default: ones) A is the arrival queue per class, A[k,r] = 1 + delta_r * Q[k,r]
Q0 (numpy.ndarray) – Initial queue length estimate (M x R) (default: proportional)
tol (float) – Convergence tolerance (default 1e-6)
max_iter (int) – Maximum iterations (default 1000)

Returns:

Q: Mean queue lengths (M x R) X: Class throughputs (R,) U: Utilizations (M x R) iter: Number of iterations performed

Return type:

Tuple of (Q, X, U, iter) where

Reference:: Schweitzer, P.J. “Approximate analysis of multiclass closed networks of queues.” Proceedings of the International Conference on Stochastic Control and Optimization (1979).

pfqn_qdlin(L, N, Z=None, tol=1e-6, max_iter=1000)[source]

QD-Linearizer (QDLIN) Approximate MVA.

Combines Queue-Dependent (QD) correction with Linearizer iteration for improved accuracy in multi-class closed networks.

Parameters:

L (numpy.ndarray) – Service demand matrix (M x R)
N (numpy.ndarray) – Population vector (R,)
Z (numpy.ndarray) – Think time vector (R,) (default: zeros)
tol (float) – Convergence tolerance (default 1e-6)
max_iter (int) – Maximum iterations (default 1000)

Returns:

Q: Mean queue lengths (M x R) U: Utilizations (M x R) R: Residence times (M x R) X: Class throughputs (1 x R) C: Cycle times (1 x R) iter: Number of iterations performed

Return type:

Tuple of (Q, U, R, X, C, iter) where

pfqn_qli(L, N, Z=None, tol=1e-6, max_iter=1000)[source]

Queue-Line (QLI) Approximate MVA (Wang-Sevcik).

Implements the Wang-Sevcik Queue-Line approximation which provides improved accuracy for multi-class networks by better estimating the queue length seen by arriving customers.

Parameters:

L (numpy.ndarray) – Service demand matrix (M x R)
N (numpy.ndarray) – Population vector (R,)
Z (numpy.ndarray) – Think time vector (R,) (default: zeros)
tol (float) – Convergence tolerance (default 1e-6)
max_iter (int) – Maximum iterations (default 1000)

Returns:

Q: Mean queue lengths (M x R) U: Utilizations (M x R) R: Residence times (M x R) X: Class throughputs (1 x R) C: Cycle times (1 x R) iter: Number of iterations performed

Return type:

Tuple of (Q, U, R, X, C, iter) where

Reference:: Wang, W. and Sevcik, K.C. “Performance Models for Multiprogrammed Systems.” IBM Research Report RC 5925 (1976).

pfqn_fli(L, N, Z=None, tol=1e-6, max_iter=1000)[source]

Fraction-Line (FLI) Approximate MVA (Wang-Sevcik).

Implements the Wang-Sevcik Fraction-Line approximation, an alternative to Queue-Line that uses a different formula for estimating the queue length seen by arriving customers.

Parameters:

L (numpy.ndarray) – Service demand matrix (M x R)
N (numpy.ndarray) – Population vector (R,)
Z (numpy.ndarray) – Think time vector (R,) (default: zeros)
tol (float) – Convergence tolerance (default 1e-6)
max_iter (int) – Maximum iterations (default 1000)

Returns:

Q: Mean queue lengths (M x R) U: Utilizations (M x R) R: Residence times (M x R) X: Class throughputs (1 x R) C: Cycle times (1 x R) iter: Number of iterations performed

Return type:

Tuple of (Q, U, R, X, C, iter) where

Reference:: Wang, W. and Sevcik, K.C. “Performance Models for Multiprogrammed Systems.” IBM Research Report RC 5925 (1976).

pfqn_bsfcfs(L, N, Z=None, tol=1e-6, max_iter=1000, QN=None, weight=None)[source]

Bard-Schweitzer approximate MVA for FCFS scheduling with weighted priorities.

Implements AMVA with FCFS approximation where classes can have relative priority weights affecting the expected waiting times.

Parameters:

L (numpy.ndarray) – Service demand matrix (M x R) where M is stations, R is classes
N (numpy.ndarray) – Population vector (1 x R) - number of jobs per class
Z (numpy.ndarray) – Think time vector (1 x R) - default: zeros
tol (float) – Convergence tolerance (default 1e-6)
max_iter (int) – Maximum iterations (default 1000)
QN (numpy.ndarray) – Initial queue length matrix (M x R) - default: uniform distribution
weight (numpy.ndarray) – Weight matrix (M x R) for relative priorities - default: ones

Returns:

XN: System throughput per class (1 x R) QN: Mean queue lengths (M x R) UN: Utilizations (M x R) RN: Residence times (M x R) it: Number of iterations performed

Return type:

Tuple of (XN, QN, UN, RN, it) where

Reference:: Bard, Y. and Schweitzer, P.J. “Analyzing Closed Queueing Networks with Multiple Job Classes and Multiserver Stations.” Performance Evaluation Review 7.1-2 (1978).

pfqn_joint(n, L, N, Z=None, lGN=None)[source]

Compute joint queue-length probability distribution.

Computes the joint probability for a given queue-length state vector in a closed product-form queueing network.

Parameters:

n (numpy.ndarray) – Queue-length state vector (M,) for total or (M x R) for per-class - If 1D (M,): n[i] is the total number of jobs at station i - If 2D (M x R): n[i,r] is the number of class-r jobs at station i
L (numpy.ndarray) – Service demand matrix (M x R)
N (numpy.ndarray) – Population vector (1 x R)
Z (numpy.ndarray) – Think time vector (1 x R) - default: zeros
lGN (float) – Log normalizing constant (optional, computed if not provided)

Returns:

Joint probability of state n

Return type:

pjoint

Examples

# Total queue-lengths (Z > 0) >>> p = pfqn_joint([2, 1], [[10, 2], [5, 4]], [2, 2], [91, 92])

# Per-class queue-lengths >>> p = pfqn_joint([[1, 0], [0, 1]], [[10, 2], [5, 4]], [2, 2], [91, 92])

pfqn_ca(L, N, Z=None)[source]

Convolution Algorithm for normalizing constant computation.

Computes the normalizing constant G(N) for a closed product-form queueing network using Buzen’s convolution algorithm.

Parameters:

L (numpy.ndarray) – Service demand matrix (M x R) where M is stations, R is classes
N (numpy.ndarray) – Population vector (1 x R or R,) - number of jobs per class
Z (numpy.ndarray) – Think time vector (1 x R or R,) - think time per class (default 0)

Returns:

G: Normalizing constant
lG: log(G)

Return type:

Tuple (G, lG) where

pfqn_is(L, N, Z=None, options=None)[source]

Importance-sampling (IS) estimate of the normalizing constant of a closed LOAD-INDEPENDENT product-form queueing network with M single-server queues of per-class demand L and an aggregated delay of think time Z.

This is the load-independent case of pfqn_ld_is() (capacities mu_i(k)=1), and the ordinary-network counterpart of the order-independent pfqn_oi_is and the pass-and-swap pfqn_pas_is: all four are the same sample-an-ordering estimator, differing only in the per-position factor of each station’s balance function. For a single-server queue that factor is the demand of the class at that position, L(i,q_p); for the delay it is Z(q_p)/p; for an OI/P&S station it is the reciprocal rank rate 1/mu_i(supp(q_1..q_p)).

With ell = sum(N), an ordering c of all ell jobs is drawn by placing a uniformly random present class at each step (probability p(c) = product of the reciprocal branching factors), and the sum over ALL ways of cutting c into contiguous per-station segments is computed exactly by dynamic programming:

G(N) = E_{C~p}[ S(C)/p(C) ],
S(c) = sum_{cuts} prod_m prod_p L(m, seg_m(p))

which is unbiased for the exact constant of pfqn_nc().

Parameters:

L ((M, R) array) – Per-class service demands at the M single-server queues.
N ((R,) array) – Closed population vector, finite.
Z ((R,) array, optional) – Aggregated think time (delay) demand; None or zeros if none.
options (dict or options object, optional) – Fields samples (default 1e4) and seed (optional).

Returns:

(G, lG) – IS estimate of the normalizing constant and its logarithm.

Return type:

tuple of float

Examples

>>> L = np.array([[0.5, 0.3], [0.2, 0.4]]); N = np.array([3, 2]); Z = np.array([1.0, 1.0])
>>> G, lG = pfqn_is(L, N, Z, {'samples': 100000})

See also

pfqn_ld_is, pfqn_nc

pfqn_nc(L, N, Z=None, method='ca', options=None)[source]

Normalizing constant computation dispatcher.

Selects appropriate algorithm based on method parameter.

Parameters:

L (numpy.ndarray) – Service demand matrix (M x R)
N (numpy.ndarray) – Population vector (R,)
Z (numpy.ndarray) – Think time vector (R,) (default: zeros)
method (str) – Algorithm to use: - ‘ca’, ‘exact’: Convolution algorithm - ‘default’: Auto-select based on problem size - ‘le’: Leading eigenvalue asymptotic - ‘cub’: Controllable upper bound - ‘imci’: Importance sampling Monte Carlo integration - ‘panacea’: Hybrid convolution/MVA - ‘propfair’: Proportionally fair allocation - ‘mmint2’: Gauss-Legendre quadrature - ‘gleint’: Gauss-Legendre integration - ‘sampling’: Monte Carlo sampling - ‘kt’: Knessl-Tier expansion - ‘comom’: Conditional moments - ‘rd’: Reduction heuristic - ‘ls’: Linearizer

Returns:

Tuple (G, lG) - normalizing constant and its log

Return type:

Tuple[float, float]

pfqn_nc_resolved_method(method)[source]

The algorithm pfqn_nc actually runs for METHOD.

pfqn_nc substitutes a different algorithm for some method names, so the requested name is not always the one that ran. Callers that report the method to the user must resolve it through here, otherwise the banner names an algorithm that never executed.

Currently only ‘comom’ substitutes: the native pfqn_comomrm port is not numerically robust for R>1, so convolution is used instead, which is exact for the single-station product-form models CoMoM-RM targets. The substitution is unconditional, hence resolvable without solving. Note this is BROADER than MATLAB, whose pfqn_nc reports ‘ca’ only for the R==1 case and genuinely runs CoMoM-RM for R>1 with a single queue.

pfqn_panacea(L, N, Z=None)[source]

PANACEA algorithm (hybrid convolution/MVA).

Currently implemented as wrapper around convolution algorithm.

Parameters:

L (numpy.ndarray) – Service demand matrix
N (numpy.ndarray) – Population vector
Z (numpy.ndarray) – Think time vector

Returns:

Tuple (G, lG) - normalizing constant and its log

Return type:

Tuple[float, float]

pfqn_propfair(L, N, Z=None)[source]

Proportionally Fair allocation approximation for normalizing constant.

Estimates the normalizing constant using a convex optimization program that is asymptotically exact in models with single-server PS queues only.

This method is based on Schweitzer’s approach and Walton’s proportional fairness theory for multi-class networks.

Parameters:

L (numpy.ndarray) – Service demand matrix (M x R) where M is stations, R is classes
N (numpy.ndarray) – Population vector (1 x R or R,) - number of jobs per class
Z (numpy.ndarray) – Think time vector (1 x R or R,) - think time per class (default 0)

Returns:

G: Estimated normalizing constant
lG: log(G)
X: Asymptotic throughputs per class (1 x R)

Return type:

Tuple (G, lG, X) where

References

Schweitzer, P. J. (1979). Approximate analysis of multiclass closed networks of queues. In Proceedings of the International Conference on Stochastic Control and Optimization.

Walton, N. (2009). Proportional fairness and its relationship with multi-class queueing networks.

pfqn_ls(L, N, Z=None, I=100000)[source]

Logistic sampling approximation for normalizing constant.

Approximates the normalizing constant using importance sampling from a multivariate normal distribution fitted at the leading eigenvalue mode.

This method is particularly effective for large networks where convolution becomes computationally expensive.

Parameters:

L (numpy.ndarray) – Service demand matrix (M x R)
N (numpy.ndarray) – Population vector (R,)
Z (numpy.ndarray) – Think time vector (R,) (default: zeros)
I (int) – Number of samples for Monte Carlo integration (default: 100000)

Returns:

G: Estimated normalizing constant lG: log(G)

Return type:

Tuple (G, lG) where

Reference:: G. Casale. “Accelerating performance inference over closed systems by asymptotic methods.” ACM SIGMETRICS 2017.

pfqn_clw(L, N, Z=None, m=None, l=None, gamma=None)[source]

Choudhury-Leung-Whitt normalization constant by numerical inversion of the generating function (JACM 42(5):935-970, 1995).

Computes g(K) of a multichain closed product-form network with single-server and (optionally) infinite-server queues by numerically inverting its p-dimensional generating function (eq. 4.5)

G(z) = exp(sum_j rho_{j0} z_j) / prod_i (1 - sum_j rho_{ji} z_j)^{m_i}

where j=1..p indexes chains, i=1..q’ the distinct single-server queues with multiplicity m_i. g(K) is recovered by p nested one-dimensional lattice-Poisson inversions (eq. 2.3) with restrictive static scaling (eqs. 5.41-5.46) and log-domain recovery (eq. 7.1).

Parameters:

L (numpy.ndarray) – (q’ x p) single-server relative traffic intensities, L[i,j]=rho_{ji}.
N (numpy.ndarray) – (p,) closed-chain population vector K.
Z (numpy.ndarray) – (p,) aggregate infinite-server relative intensities rho_{j0}. Default 0.
m (numpy.ndarray) – (q’,) queue multiplicities m_i. Default ones.
l (numpy.ndarray) – (p,) inner lattice parameters l_j. Default 1,2,2,3,3,…
gamma (numpy.ndarray) – (p,) aliasing parameters gamma_j. Default 11,13,13,15,15,…

Returns:

normalization constant (inf if it overflows double) and its natural logarithm (always finite).

Return type:

Tuple (G, lG)

Note: exact nested inversion of cost prod_j 2 l_j K_j; practical for moderate populations and few chains. The paper’s Euler summation and dimension reduction speed-ups are not applied here.

pfqn_clw_lld(L, N, Z=None, mu=None, l=None, gamma=None)[source]

Choudhury-Leung-Whitt normalization constant by numerical inversion of the generating function (JACM 42(5):935-970, 1995), extended to limited load-dependent (LLD) stations via the per-center transforms of Bertozzi and McKenna (SIAM Review 35(2):239-268, 1993).

The generating function is (Bertozzi-McKenna eqs. 2.17/2.23)

G(z) = exp(sum_j rho_{j0} z_j) prod_i F_i(sum_j rho_{ji} z_j)

where F_i is the transform of the station factor of queue i (eq. 2.16) with load-dependent rate scalings S_i(k) = mu[i,k]. For an LLD queue, S_i(k) = c_i constant for k >= l_i, and F_i is the rational function (eq. 2.19)

F_i(x) = [c_i + sum_{n=1}^{l_i-1} (c_i - S_i(n))
/ prod_{k=1}^n S_i(k) * x^n] / (c_i - x),

analytic except for a simple pole at x = c_i. Multiserver and load-independent queues are special cases. Since g(K) depends on S_i(k) only for k <= sum(K), general load-dependent input is truncated to LLD at sum(K) without loss of exactness.

g(K) is recovered by p nested one-dimensional lattice-Poisson inversions (CLW eq. 2.3) with restrictive static scaling adapted from CLW eqs. 5.41-5.46 (each queue normalized by its pole c_i, simple pole) and log-domain recovery (eq. 7.1).

Parameters:

L (numpy.ndarray) – (q’ x p) single-server relative traffic intensities, L[i,j]=rho_{ji}.
N (numpy.ndarray) – (p,) closed-chain population vector K.
Z (numpy.ndarray) – (p,) aggregate infinite-server relative intensities rho_{j0}. Default 0.
mu (numpy.ndarray) – (q’ x n) load-dependent rate scalings mu[i,k] = S_i(k+1); if fewer than sum(N) columns are given the last column is extended (LLD assumption). Default ones (all queues load-independent).
l (numpy.ndarray) – (p,) inner lattice parameters l_j. Default 1,2,2,3,3,…
gamma (numpy.ndarray) – (p,) aliasing parameters gamma_j. Default 11,13,13,15,15,…

Returns:

normalization constant (inf if it overflows double) and its natural logarithm (always finite).

Return type:

Tuple (G, lG)

Note: cost is prod_j 2 l_j K_j contour points, each of cost O(sum_i l_i); practical for moderate populations and few chains.

pfqn_linearizer(L, N, Z, sched_type, tol=1e-8, maxiter=1000, QN0=None)[source]

Linearizer approximate MVA algorithm.

Parameters:

L (numpy.ndarray) – Demand matrix (M x R) - rows are stations, columns are classes
N (numpy.ndarray) – Population vector (R,)
Z (numpy.ndarray) – Think time vector (R,)
sched_type (List[str]) – List of scheduling strategies per station (‘FCFS’, ‘PS’, etc.)
tol (float) – Convergence tolerance
maxiter (int) – Maximum iterations

Returns:

Q: Queue lengths (M x R) U: Utilizations (M x R) W: Waiting times (M x R) T: Station throughputs (M x R) C: Response times (1 x R) X: Class throughputs (1 x R) iterations: Number of iterations

Return type:

Tuple of (Q, U, W, T, C, X, iterations)

pfqn_gflinearizer(L, N, Z, sched_type, tol=1e-8, maxiter=1000, alpha=1.0, QN0=None)[source]

General-form linearizer approximate MVA.

Parameters:

L (numpy.ndarray) – Demand matrix (M x R)
N (numpy.ndarray) – Population vector (R,)
Z (numpy.ndarray) – Think time vector (R,)
sched_type (List[str]) – List of scheduling strategies per station
tol (float) – Convergence tolerance
maxiter (int) – Maximum iterations
alpha (float) – Linearization parameter (scalar)

Returns:

Same as pfqn_linearizer

Return type:

Tuple[numpy.ndarray, numpy.ndarray, numpy.ndarray, numpy.ndarray, numpy.ndarray, numpy.ndarray, int]

pfqn_egflinearizer(L, N, Z, sched_type, tol=1e-8, maxiter=1000, alpha=None, QN0=None)[source]

Extended general-form linearizer with class-specific parameters.

Parameters:

L (numpy.ndarray) – Demand matrix (M x R)
N (numpy.ndarray) – Population vector (R,)
Z (numpy.ndarray) – Think time vector (R,)
sched_type (List[str]) – List of scheduling strategies per station
tol (float) – Convergence tolerance
maxiter (int) – Maximum iterations
alpha (numpy.ndarray) – Class-specific linearization parameters (R,)

Returns:

Tuple of (Q, U, W, T, C, X, iterations)

Return type:

Tuple[numpy.ndarray, numpy.ndarray, numpy.ndarray, numpy.ndarray, numpy.ndarray, numpy.ndarray, int]

pfqn_momlin(L, N, Z=None, tol=1e-8, maxiter=1000)[source]

Approximate first and second queue-length moments of a closed product-form network, scalable to large populations and many classes.

Parameters:

L ((M, R) array) – Service demand matrix.
N ((R,) array) – Closed population per class.
Z ((R,) array, optional) – Think time per class (default zeros).
tol (float) – Convergence tolerance on the queue-length fixed point.
maxiter (int) – Maximum iterations.

Return type:

PfqnMomlin

class PfqnMomlin(Q, X, U, R, QVar, QCov, dQ)[source]

Bases: object

Result container for pfqn_momlin().

Q

Mean queue length.

Type:: np.ndarray (M x R)

X

Throughput per class.

Type:: np.ndarray (1 x R)

U, R

Utilization and residence time.

Type:: np.ndarray (M x R)

QVar

Queue-length variance.

Type:: np.ndarray (M x R)

QCov

Queue-length covariance, QCov[i, r, j, s] = Cov[n_ir, n_js].

Type:: np.ndarray (M x R x M x R)

dQ

Demand derivatives, dQ[i, r, j, s] = dQ_ir/dD_js.

Type:: np.ndarray (M x R x M x R)

class SchedStrategy(*values)[source]

Bases: Enum

Scheduling strategies for queueing stations.

pfqn_sens(L, N, Z=None, mi=None)[source]

Exact derivatives of {X, Q, U, R} w.r.t. demands L and think times Z.

Derivatives are analytic (exact to machine precision), not finite differences. The CoMoM-backed kernel is selected for the repairman model (a single single-server queue, M=1, plus a think-time delay where every populated class has positive think time and demand); it is polynomial in the number of classes R. Every other model uses forward-mode differentiation of the exact MVA recursion.

Parameters:

L ((M, R) array) – Service demand matrix, L[i, r] = visits_ir / rate_ir.
N ((R,) array) – Population per class.
Z ((R,) array, optional) – Think time per class (default zeros).
mi ((M,) array, optional) – Station residence multiplicity (default ones).

Return type:

PfqnSens

class PfqnSens(X, Q, U, R, params, dX, dQ, dU, dR, QCov=None, QVar=None, QTotVar=None, QCovAsym=None)[source]

Bases: object

Result container for pfqn_sens().

X, Q, U, R

Base MVA measures. X is (1 x R) system throughput per class, Q (M x R) mean queue length, U (M x R) utilization, R (M x R) residence time per visit.

Type:: np.ndarray

params

One entry per differentiation parameter p with keys type (‘L’ or ‘Z’), station (i, or -1 for Z) and jobclass (r).

Type:: list of dict

dX

Type:: np.ndarray (R x P)

dQ, dU, dR

Derivative of each base measure w.r.t. parameter p.

Type:: np.ndarray (M x R x P)

QCov

Queue-length covariance, QCov[i, r, j, s] = Cov[n_ir, n_js] = D_js dQ_ir/dD_js. The same-station blocks (i == j) come from pfqn_sens_mva(); the cross-station ones are read off the Jacobian.

Type:: np.ndarray (M x R x M x R)

QVar

Queue-length variance, QVar[i, r] = QCov[i, r, i, r].

Type:: np.ndarray (M x R)

QTotVar

Variance of the total queue length per station, QTotVar[i] = Var[sum_r n_ir].

Type:: np.ndarray (M,)

QCovAsym

Roundoff-level residual of the moment recursion, see pfqn_sens_mva().

Type:: float

pfqn_sens_mva(L, N, Z=None, mi=None)[source]

Exact second moments of the queue lengths of a closed product-form network.

The moments are obtained by an MVA-type recursion evaluated on the same population lattice as pfqn_mva(), so no derivative of the model is ever formed and the cost is O(M*R^2) per lattice point rather than the O(M^2*R^2) of the differentiated-MVA kernel used by pfqn_sens().

The recursion is obtained by differentiating the Reiser-Lavenberg MVA equation Q(j,v|N) = X(v|N) L(j,v) (mi(j) + Qtot(j|N-e_v)) with respect to the visit ratio theta(i,k) of class k at station i and rescaling. Writing W(k,i;v,j|N) = Cov[n(i,k), n(j,v)] at population N,

W(k,i;v,j|N) = Q(j,v|N) (Q(i,k|N-e_v) - Q(i,k|N))

[i==j & k==v] Q(j,v|N)

X(v|N) L(j,v) sum_t W(k,i;t,j|N-e_v)

with W(.|0) = 0. This routine evaluates the same-station case i==j, which is self-contained: the inner sum then only involves same-station terms, so a single scalar Ssum(j,k|N) = sum_t W(k,j;t,j|N) carried along the lattice closes the recursion. The cross-station case i != j is not self-contained (it couples every station pair) and costs as much as the full Jacobian, so it is left to pfqn_sens().

Setting i == j and mi == 1 reproduces Corollary 1 of de Souza e Silva and Muntz, i.e. its equations (2.9a) for the variance and (2.10) for the covariance; the station multiplicity mi cancels identically because X(v|N) L(j,v) (mi(j) + Qtot(j|N-e_v)) = Q(j,v|N) is the MVA equation for any mi. The equivalent statement for an infinite-server station, equation (2.9b), is recovered automatically because LINE folds the delay into the think time Z, which enters only through X(v|N) and carries no queue-length moment of its own.

Parameters:

L ((M, R) array) – Service demand matrix, L[i, r] = visits_ir / rate_ir.
N ((R,) array) – Population per class. Must be finite (closed populations only).
Z ((R,) array, optional) – Think time per class (default zeros).
mi ((M,) array, optional) – Server multiplicity per station (default ones).

Return type:

PfqnSensMva

Notes

Restricted to closed populations. Mixed and load-dependent models are handled by pfqn_sens_mvaldmx().

For a station of multiplicity mi[i] > 1, which LINE treats as mi[i] identical replicas sharing the demand row L[i, :], the moments returned are those of the aggregate queue length over the replicas.

class PfqnSensMva(X, Q, U, R, QCov, QVar, QTotVar, QCovAsym)[source]

Bases: object

Result container for pfqn_sens_mva().

X, Q, U, R

Base MVA measures, identical entry by entry to pfqn_mva(L, N, Z, mi). X is (1 x R) system throughput per class, Q (M x R) mean queue length, U (M x R) utilization, R (M x R) residence time per visit.

Type:: np.ndarray

QCov

QCov[i, r, s] = Cov[n(i,r), n(i,s)], the queue-length covariance of classes r and s at station i. Symmetric in (r, s).

Type:: np.ndarray (M x R x R)

QVar

QVar[i, r] = QCov[i, r, r] = Var[n(i,r)].

Type:: np.ndarray (M x R)

QTotVar

QTotVar[i] = Var[sum_r n(i,r)], i.e. sum_{r,s} QCov[i, r, s]. This is Theorem 3 of de Souza e Silva and Muntz, obtained here without a capacity derivative.

Type:: np.ndarray (M,)

QCovAsym

max |W(r,s) - W(s,r)| over the covariance entries before symmetrization. The two triangles come from differentiating two different classes’ MVA equations, so this is an independent residual of the recursion and should sit at roundoff; a large value signals a bug.

Type:: float

pfqn_sens_ldmx_ec(lam, D, mu)[source]

Effective capacity terms of the mixed load-dependent MVA and their exact derivatives with respect to the open-class load.

Computes the terms EC, E and Eprime of the mixed load-dependent MVA of Bruell-Balbo-Afshari, exactly as pfqn_ldmx_ec() does, and additionally their exact analytic derivatives with respect to the open-class load Lo[i] of each station.

Lo[i] = sum_r lambda[r] D[i, r] is the only channel through which a service demand enters E, Eprime and EC: the load-dependent rates mu are independent of the demands. Station i’s terms depend on Lo[i] alone, so a single derivative per station is enough, and the chain rule then yields the derivative with respect to any demand-scaling parameter. This is the factorization behind equations (19), (21) and (24)-(31) of Akyildiz and Strelen, which are reproduced here term by term.

Parameters:

lam ((R,) array) – Arrival rate vector. Zero for closed classes.
D ((M, R) array) – Service demand matrix.
mu ((M, Nt) array) – Load-dependent rate matrix, limited load dependence.

Returns:

EC ((M, Nt) array) – Effective capacity matrix.
E ((M, 1 + Nt) array) – E-function values, E[i, n] holding the MATLAB E(i, 1+n).
Eprime ((M, 1 + Nt) array) – E-prime function values, indexed as E.
Lo ((M,) array) – Open class load vector.
dEC ((M, Nt) array) – dEC[i, n] = dEC(i, n)/dLo(i).
dE ((M, 1 + Nt) array) – dE[i, n] = dE(i, n)/dLo(i).
dEprime ((M, 1 + Nt) array) – dEprime[i, n] = dEprime(i, n)/dLo(i).

Return type:

Tuple[numpy.ndarray, numpy.ndarray, numpy.ndarray, numpy.ndarray, numpy.ndarray, numpy.ndarray, numpy.ndarray]

pfqn_sens_mvaldmx(lam, D, N, Z, mu=None, S=None)[source]

Exact second moments of the queue lengths of a mixed open/closed product-form network with limited load-dependent service rates.

This is the load-dependent and mixed counterpart of pfqn_sens_mva(), which is restricted to closed load-independent models.

The method is the moment analysis of Akyildiz and Strelen. Their Theorem 1, equation (11), states that multiplying a queue-length moment by one further factor Q_jT costs one derivative with respect to a parameter y_j that scales the service demands s_ir of the classes r in T at station j:

E[Q_jT^k …] = d/dy_j E[Q_jT^(k-1) …]|_{y_j=1}

nbar_jT E[Q_jT^(k-1) …]

Taking k=2 and T={s} gives the second moment, hence

Cov[n(i,r), n(j,s)] = d nbar(i,r) / dy_(j,s) |_{y=1}

which is evaluated here by forward-mode differentiation of the mixed load-dependent MVA of Bruell-Balbo-Afshari, i.e. of exactly the recursion implemented by pfqn_mvaldmx(). The differentiated equations are (13) for the residence times, (15)-(17) for the conditional marginal probabilities, (18) for the throughputs, (19) and (24)-(31) for the effective capacities (delegated to pfqn_sens_ldmx_ec()), (32) for the closed-class queue lengths and (33) for the open-class ones.

Because a demand-scaling parameter perturbs the whole network through the closed-class throughputs, the derivatives must be propagated for every parameter, so the cross-station covariances come out at no extra cost and are returned in QCovFull. This is unlike pfqn_sens_mva(), whose cheaper same-station recursion cannot reach them.

Parameters:

lam ((R,) array) – Arrival rate vector. Must be zero on closed classes.
D ((M, R) array) – Service demand matrix.
N ((R,) array) – Population vector. inf entries denote open classes.
Z ((R,) array) – Think time vector.
mu ((M, sum(N)) array, optional) – Load-dependent rate matrix, limited load dependence (default ones).
S ((M,) array, optional) – Number of servers per station. Accepted for signature compatibility with pfqn_mvaldmx(), which likewise does not read it: the multiserver behaviour is carried entirely by the rates mu.

Return type:

PfqnSensMvaldmx

Notes

Open classes are supported: an open class contributes to the load Lo[i] that drives the effective capacities, and equation (21) supplies the corresponding dLo/dy.

The moments of an open class are those of its queue length at a station, which is finite even though its population is infinite.

class PfqnSensMvaldmx(X, Q, U, R, QCov, QCovFull, QVar, QTotVar, QCovAsym)[source]

Bases: object

Result container for pfqn_sens_mvaldmx().

X, Q, U, R

Base measures, identical to pfqn_mvaldmx(lambda, D, N, Z, mu, S). X is (R,) throughput per class, Q (M x R) mean queue length, U (M x R) utilization, R (M x R) residence time.

Type:: np.ndarray

QCov

QCov[i, r, s] = Cov[n(i,r), n(i,s)], the same-station block.

Type:: np.ndarray (M x R x R)

QCovFull

QCovFull[i, r, j, s] = Cov[n(i,r), n(j,s)].

Type:: np.ndarray (M x R x M x R)

QVar

QVar[i, r] = Var[n(i,r)].

Type:: np.ndarray (M x R)

QTotVar

QTotVar[i] = Var[sum_r n(i,r)].

Type:: np.ndarray (M,)

QCovAsym

max |QCovFull[i, r, j, s] - QCovFull[j, s, i, r]| before symmetrization. The two entries are produced by differentiating two different classes’ equations, so this residual is an independent check of the recursion and should sit at roundoff.

Type:: float

pfqn_sens_mom(L, N, Z=None, mi=None, groups=None)[source]

Exact moments E[Q_i], E[Q_i^2], E[Q_i^3] and the covariances Cov[Q_i,Q_j] of the TOTAL queue lengths Q_i = sum_r n(i,r) of a closed product-form (BCMP) queueing network.

The method is the moment analysis of Strelen. Its Theorem 3.1 states that one further factor Q_i in a moment costs one differentiation with respect to x_i, the reciprocal of the capacity of station i, i.e. a parameter that scales the service times of ALL classes at station i:

E[Q_i^j] = m_i E[Q_i^(j-1)] + x_i d/dx_i E[Q_i^(j-1)],   Q_i^0 = 1

Iterating from E[Q_i^0] = 1 gives, with m_i = E[Q_i], equation (3.2):

Var[Q_i]     = x_i dm_i/dx_i
Cov[Q_i,Q_j] = x_j dm_i/dx_j = x_i dm_j/dx_i
E[Q_i^2]     = x_i dm_i/dx_i + m_i^2
E[Q_i^3]     = x_i^2 d^2m_i/dx_i^2 + (x_i + 3 x_i m_i) dm_i/dx_i + m_i^3

so the third moment requires the SECOND derivative of the MVA recursion, which is what this routine adds over pfqn_sens_mva() and pfqn_sens() (both first order only). The derivatives are obtained by second-order forward-mode differentiation of the Reiser-Lavenberg recursion, i.e. by carrying, for each parameter, the value together with its first and second derivative along the population lattice (Theorem 3.2 for one class, Theorem 3.5 for several).

The parameter need not scale a whole column. Theorem 1 of Akyildiz and Strelen states the same recursion for a parameter that scales the service times of an arbitrary class subset T at station i, and the moments it generates are then those of Q_(i,T) = sum_(r in T) n(i,r). groups supplies that subset structure: it partitions the classes, and the routine reports the moments of each group’s queue length at each station. The three useful settings are:

groups = ones(R)   the whole column: per-station TOTALS (default, this
                   is Strelen's x_i)
groups = 1..R      one class per group: PER-CLASS moments, so that even
                   the third moment is per class
groups = chain(r)  one group per chain: PER-CHAIN moments

Strelen states only the first; the generalization is Akyildiz and Strelen’s T. The per-class setting reproduces pfqn_sens_mva()’s second moments exactly.

Parameters:

L ((M, R) array) – Service demand matrix, L[i, r] = visits_ir / rate_ir.
N ((R,) array) – Population per class. Must be finite (closed populations only).
Z ((R,) array, optional) – Think time per class (default zeros).
mi ((M,) array, optional) – Server multiplicity per station (default ones).
groups ((R,) array, optional) – Class-to-group map, a partition of the classes into G = max(groups) groups labelled consecutively 1..G with no empty group. Default ones(R), i.e. one group holding every class, the per-station total. Labels are 1-based, as in MATLAB.

Return type:

PfqnSensMom

Notes

Restricted to closed populations, as is the moment analysis of the reference. Mixed and load-dependent second moments are in pfqn_sens_mvaldmx().

Moments of the sojourn times at FCFS centers are built on top of these queue-length moments by pfqn_sens_respt(), following Theorem 4.1.

The exact recursion costs O(prod(N+1)) lattice points; pfqn_sens_linearizer() approximates the same quantities in polynomial time.

class PfqnSensMom(X, Q, U, R, m, dm, d2m, Var, Cov, M2, M3, Skew, CovAsym)[source]

Bases: object

Result container for pfqn_sens_mom().

X, Q, U, R

Base MVA measures, identical entry by entry to pfqn_mva(L, N, Z, mi). X is (1 x R), Q, U and R are (M x R).

Type:: np.ndarray

m

m[i, g] = E[Q_(i,g)], the mean queue length of group g at station i. COLLAPSED to (M,) in the default single-group case, where it is the per-station total.

Type:: np.ndarray (M, G)

dm

dm[i, g, j, g2], the scaled first derivative of m_(i,g) with respect to the parameter of (j, g2). Collapsed to (M, M) when G == 1.

Type:: np.ndarray (M, G, M, G)

d2m

d2m[i, g], the scaled pure second derivative with respect to the parameter of (i, g). Only that entry is needed by (3.2); the mixed second derivatives are not required for moments of a single Q_(i,g) and would cost an extra factor M*G to carry. Collapsed to (M,) when G == 1.

Type:: np.ndarray (M, G)

Var

Var[Q_(i,g)]. Collapsed to (M,) when G == 1.

Type:: np.ndarray (M, G)

Cov

Cov[i, g, j, g2] = Cov[Q_(i,g), Q_(j,g2)]. Collapsed to (M, M) when G == 1, where the group index carries no information.

Type:: np.ndarray (M, G, M, G)

M2

E[Q_(i,g)^2]. Collapsed to (M,) when G == 1.

Type:: np.ndarray (M, G)

M3

E[Q_(i,g)^3]. Collapsed to (M,) when G == 1.

Type:: np.ndarray (M, G)

Skew

Skewness of Q_(i,g), i.e. the third central moment divided by Var^(3/2). NaN where the variance is zero (a deterministic queue length). Collapsed to (M,) when G == 1.

Type:: np.ndarray (M, G)

CovAsym

max |x_j dm_i/dx_j - x_i dm_j/dx_i| before symmetrization of Cov. The two are distinct expressions that must agree, so this is a live residual of the recursion; expect roundoff.

Type:: float

pfqn_sens_respt(S, V, N, Z=None, b=None, tmax=3)[source]

Exact raw moments E[W_(i,l)^t], t = 1..tmax, of the sojourn time of a class-l job at an FCFS b-server center i of a closed product-form queueing network, together with the variance of that sojourn time.

This is Theorem 4.1 of the reference. Its mechanism is the arrival theorem of Lavenberg-Reiser and Sevcik-Mitrani: a class-l job arriving at center i finds j jobs already there with probability p_i(j, N - 1_l). Conditioning the sojourn time on j and inverting the Laplace transform of the conditional density gives:

E[W_(i,l)^t] = t!/mu^t + sum_{tau=0..t} a_(t,tau)(0) E[Qt_i^tau]
               - sum_{j=0..b-1} p_i(j,N-1_l) sum_{tau=0..t} a_(t,tau)(0) j^tau

where mu = 1/S(i) is the rate of each of the b servers, Qt_i is the total queue length at center i at population N - 1_l (so its moments are those of pfqn_sens_mom() evaluated one job down in class l), and the coefficients a_(t,tau)(0) depend only on b and mu, not on the network (Remark 4.3 of the reference). The moments E[Qt_i^tau] up to tau = 3 need the second derivative of the MVA recursion, so this routine carries a second-order forward-mode pass exactly as pfqn_sens_mom() does, but over the b-server recursion (4.1)-(4.2) rather than the single-server one.

For b = 1 the coefficients a_(t,0)(0) vanish identically and the double-sum correction disappears, so no marginal probabilities are needed (Remark 4.2 of the reference); the routine still evaluates the general expression, which reduces to that case on its own.

Only FCFS centers are covered. The reference is explicit that the sojourn-time distribution at PS and LCFS centers is in general not known, so no analogue exists there. FCFS in a BCMP network further requires the service time to be exponential and class-independent, which is why this routine takes a per-station service time S(i) and a separate visit-ratio matrix V rather than a demand matrix: the sojourn time is per visit, so the per-visit rate mu = 1/S(i) must be known and cannot be recovered from the demand L(i,l) = S(i)*V(i,l) alone.

Parameters:

S ((M,) array) – Service time at each station, common to all classes.
V ((M, R) array) – Visit ratio matrix. The demand is L[i, r] = S[i] * V[i, r].
N ((R,) array) – Population per class. Must be finite (closed populations only).
Z ((R,) array, optional) – Think time per class (default zeros).
b ((M,) array, optional) – Number of servers at each station (default ones).
tmax (int, optional) – Highest sojourn-time moment to return, 1..3 (default 3). The coefficients a_(t,tau)(0) are tabulated in the reference up to t = 3.

Return type:

PfqnSensRespt

Notes

Restricted to closed populations. Load-dependent rates are not covered here; the b-server dependence is the only state dependence, and it is carried exactly by (4.1)-(4.2).

class PfqnSensRespt(X, Q, U, m, Var, p, W, WM, Wresid, WVar, WSkew)[source]

Bases: object

Result container for pfqn_sens_respt().

X, Q, U

Base measures at population N. X is (1 x R), Q and U are (M x R).

Type:: np.ndarray

W

W[i, l] = E[W_(i,l)], the mean sojourn time per visit of a class-l job at station i. Zero where class l does not visit i.

Type:: np.ndarray (M, R)

WM

WM[i, l, t-1] = E[W_(i,l)^t].

Type:: np.ndarray (M, R, tmax)

WVar

Var[W_(i,l)] = E[W^2] - E[W]^2. Requires tmax >= 2.

Type:: np.ndarray (M, R)

WSkew

Skewness of W_(i,l). Requires tmax >= 3; NaN if the variance is zero.

Type:: np.ndarray (M, R)

m

E[Q_i] at population N, the total queue length.

Type:: np.ndarray (M,)

Var

Var[Q_i] at population N.

Type:: np.ndarray (M,)

p

p[i, j] = P[Q_i = j] at population N, for j = 0..b_i-1. These are the only marginal probabilities the b-server recursion needs, so the matrix is RAGGED: row i is meaningful only up to column b_i and is zero-padded out to max(b). A padded entry is not P[Q_i = j]; it is simply not computed. Read row i as p[i, :b[i]].

Type:: np.ndarray (M, max(b))

Wresid

The residence time w_i(l) of the MVA recursion. The identity W[i, l] = Wresid[i, l] / V[i, l] is an independent check of the t = 1 case of (4.5).

Type:: np.ndarray (M, R)

pfqn_sens_linearizer(L, N, Z=None, tol=None, maxiter=200)[source]

Approximate moments E[Q_i], Var[Q_i], Cov[Q_i,Q_j], E[Q_i^2] and E[Q_i^3] of the per-station total queue lengths of a closed product-form queueing network, by the LINEARIZER-2 / LINEARIZER-3 algorithms of the reference (Section 5).

Motivation. The exact moment analysis of pfqn_sens_mom() evaluates the MVA recursion on the whole population lattice, so it costs O(prod(N+1)) and is unusable once the populations are large. The Linearizer replaces that lattice by a fixed point over a handful of populations, and the reference observes that the same trick applies to the derivatives: differentiate the Linearizer equations, append the differentiated equations to the originals, and iterate all of them together. This routine does that, carrying both the first and the second derivative, so it returns everything (3.2) needs, including the third moment. Carrying only the first derivative is the reference’s LINEARIZER-2; carrying the second as well is its LINEARIZER-3.

The approximation. CORE (equations (5.1)-(5.2)) estimates the queue lengths at population n - 1_l from those at n by:

v_i(l)          = m_i^(n)(l) / n(l)
m_i^(n-1_l')(l) = (n - 1_l')_l * ( v_i(l) + delta_i(l',l) )

and substitutes them into the exact MVA equations. Setting the delta terms to zero gives Bard-Schweitzer; Linearizer instead estimates them from (5.3), delta_i^(N)(l',l) = v_i^(N-1_l')(l) - v_i^(N)(l), by running CORE at each of the N - 1_l populations, and holds them fixed across populations (the heuristic (5.4)). Differentiating (5.1)-(5.3) gives (5.5)-(5.8), which are carried alongside.

Accuracy. The reference reports, over 51 networks including 34 stress cases, relative errors below 2.1% on E[Q], 4.1% on E[Q^2] and 6.2% on E[Q^3]. The tests measure the error against the exact pfqn_sens_mom() on models small enough for both, and assert bands of that order rather than machine precision: this routine is an approximation and is expected to disagree with the exact answer.

Parameters:

L ((M, R) array) – Service demand matrix, L[i, r] = visits_ir / rate_ir.
N ((R,) array) – Population per class. Must be finite (closed populations only).
Z ((R,) array, optional) – Think time per class (default zeros).
tol (float, optional) – Convergence tolerance of the CORE fixed point on the mean queue lengths. Default: the test of the reference, 1/(4000+16*sum(n)), which is also applied to the variances at 1e-3.
maxiter (int, optional) – Maximum CORE iterations (default 200).

Return type:

PfqnSensLinearizer

Notes

Single-server stations plus an optional delay Z, matching pfqn_linearizer().

The exact counterpart is pfqn_sens_mom(); the per-class exact second moments are in pfqn_sens_mva().

class PfqnSensLinearizer(X, Q, U, W, m, dm, d2m, Var, Cov, M2, M3, Skew, CovAsym, iter)[source]

Bases: object

Result container for pfqn_sens_linearizer().

X, Q, U, W

Approximate base measures. X is (1 x R); Q, U and W are (M x R).

Type:: np.ndarray

m

Approximate E[Q_i], the total queue length at station i.

Type:: np.ndarray (M,)

dm

dm[i, h] = x_h dm_i/dx_h, the scaled first derivative.

Type:: np.ndarray (M, M)

d2m

d2m[i] = x_i^2 d^2m_i/dx_i^2.

Type:: np.ndarray (M,)

Var, Cov, M2, M3, Skew

The moments of (3.2), formed exactly as in pfqn_sens_mom() but from the approximate derivatives. Var, M2, M3 and Skew are (M,); Cov is (M x M).

Type:: np.ndarray

CovAsym

Raw asymmetry of Cov before symmetrization. Unlike the exact routines, this is NOT expected to sit at roundoff: the Linearizer fixed point does not enforce the symmetry that the product form guarantees, so this is a useful measure of the approximation error.

Type:: float

iter

Total CORE iterations performed.

Type:: int

pfqn_mvald(L, N, Z, mu, stabilize=True)[source]

Exact MVA for load-dependent closed queueing networks.

This algorithm extends standard MVA to handle stations where the service rate depends on the number of jobs present.

Parameters:

L (numpy.ndarray) – Service demand matrix (M x R) - rows are stations, columns are classes.
N (numpy.ndarray) – Population vector (R,).
Z (numpy.ndarray) – Think time vector (R,).
mu (numpy.ndarray) – Load-dependent rate matrix (M x Ntot) where mu[i,k] is the service rate at station i when k jobs are present.
stabilize (bool) – Force non-negative probabilities (default: True).

Returns:

XN: Class throughputs (1 x R). QN: Mean queue lengths (M x R). UN: Utilizations (M,), per STATION, i.e. 1-P_j(0). Matches MATLAB

pfqn_mvald, which likewise reports utilization per station and not per class.

CN: Cycle times (1 x R). lGN: Log normalizing constant evolution. isNumStable: True if numerically stable. pi: Marginal queue-length probabilities (M x Ntot+1).

Return type:

Tuple of (XN, QN, UN, CN, lGN, isNumStable, pi)

pfqn_mvams(lambda_arr, L, N, Z, mi=None, S=None)[source]

General-purpose MVA for mixed networks with multiserver nodes.

This function handles networks with open/closed classes and multi-server stations, routing to the appropriate specialized algorithm.

Parameters:

lambda_arr (numpy.ndarray) – Arrival rate vector (R,). Use 0 for closed classes.
L (numpy.ndarray) – Service demand matrix (M x R).
N (numpy.ndarray) – Population vector (R,). Use np.inf for open classes.
Z (numpy.ndarray) – Think time vector (R,).
mi (numpy.ndarray | None) – Queue replication factors (M,) (default: ones).
S (numpy.ndarray | None) – Number of servers per station (M,) (default: ones).

Returns:

XN: Class throughputs (1 x R). QN: Mean queue lengths (M x R). UN: Utilizations. Per STATION-CLASS (M x R) on every branch except

the closed multiserver one, which delegates to pfqn_mvald and so reports per STATION (M,). This shape inconsistency is inherited from MATLAB pfqn_mvams, which behaves identically; the only caller (api/solvers/mva/handler.py) discards UN and recomputes utilization analytically, as solver_mva.m does.

CN: Residence times per STATION-CLASS (M x R), as in MATLAB: pfqn_mvams and pfqn_mva. Note this is NOT what pfqn_mva returns as its own CN in Python (that is the (1 x R) cycle time); the residence time is pfqn_mva’s RN.

lG: Log normalizing constant.

Return type:

Tuple of (XN, QN, UN, CN, lG)

References

Original MATLAB: matlab/src/api/pfqn/pfqn_mvams.m

pfqn_mvamx(lambda_arr, L, N, Z, mi=None)[source]

Exact MVA for mixed open/closed single-server networks.

Handles networks with both open classes (infinite population, external arrivals) and closed classes (fixed population, no external arrivals).

Parameters:

lambda_arr (numpy.ndarray) – Arrival rate vector (R,). Use 0 for closed classes.
L (numpy.ndarray) – Service demand matrix (M x R).
N (numpy.ndarray) – Population vector (R,). Use np.inf for open classes.
Z (numpy.ndarray) – Think time vector (R,).
mi (numpy.ndarray | None) – Queue replication factors (M,) (default: ones).

Returns:

XN: Class throughputs (1 x R). QN: Mean queue lengths (M x R). UN: Utilizations (M x R). CN: Cycle times (M x R). lGN: Log normalizing constant.

Return type:

Tuple of (XN, QN, UN, CN, lGN)

pfqn_xzabalow(L, N, Z)[source]

Lower ABA (asymptotic bound analysis) bound on throughput.

Returns N / (Z + sum(L)*N), the ABA lower throughput bound for single-class closed queueing networks. This is NOT the classical Zahorjan-Balanced (balanced job bounds) lower bound; that one is pfqn_xzgsblow, which is tighter.

Parameters:

L (numpy.ndarray) – Service demand vector (M,).
N (int | float) – Population (scalar).
Z (float) – Think time.

Returns:

Lower bound on throughput.

Return type:

float

pfqn_xzabaup(L, N, Z)[source]

Upper asymptotic bound on throughput (Zahorjan-Balanced).

Provides a simple upper bound on system throughput for single-class closed queueing networks based on bottleneck analysis.

Parameters:

L (numpy.ndarray) – Service demand vector (M,).
N (int | float) – Population (scalar).
Z (float) – Think time.

Returns:

Upper bound on throughput.

Return type:

float

pfqn_qzgblow(L, N, Z, i)[source]

Lower asymptotic bound on queue length (Zahorjan-Gittelsohn-Bryant).

Parameters:

L (numpy.ndarray) – Service demand vector (M,).
N (int | float) – Population (scalar).
Z (float) – Think time.
i (int) – Station index (0-based).

Returns:

Lower bound on mean queue length at station i.

Return type:

float

pfqn_qzgbup(L, N, Z, i)[source]

Upper asymptotic bound on queue length (Zahorjan-Gittelsohn-Bryant).

Parameters:

L (numpy.ndarray) – Service demand vector (M,).
N (int | float) – Population (scalar).
Z (float) – Think time.
i (int) – Station index (0-based).

Returns:

Upper bound on mean queue length at station i.

Return type:

float

pfqn_xzgsblow(L, N, Z)[source]

Lower asymptotic bound on throughput (Zahorjan-Gittelsohn-Schweitzer-Bryant).

Provides a tighter lower bound than pfqn_xzabalow by accounting for queue length bounds.

Parameters:

L (numpy.ndarray) – Service demand vector (M,).
N (int | float) – Population (scalar).
Z (float) – Think time.

Returns:

Lower bound on throughput.

Return type:

float

pfqn_xzgsbup(L, N, Z)[source]

Upper asymptotic bound on throughput (Zahorjan-Gittelsohn-Schweitzer-Bryant).

Provides a tighter upper bound than pfqn_xzabaup by accounting for queue length bounds.

Parameters:

L (numpy.ndarray) – Service demand vector (M,).
N (int | float) – Population (scalar).
Z (float) – Think time.

Returns:

Upper bound on throughput.

Return type:

float

pfqn_mwrbb(V, S, N, Z=None, sched=None, prio=None)[source]

Majumdar-Woodside robust box bounds on throughput for closed multiclass queueing networks with mixed scheduling disciplines.

Computes distribution-insensitive (NBUE) upper and lower bounds on the per-class system throughput of a closed multiclass queueing network, per S. Majumdar and C.M. Woodside, “Robust bounds and throughput guarantees for closed multiclass queueing networks”, Performance Evaluation 32 (1998) 101-136. The upper bound intersects the no-contention bound (eq. 2) with the utilization-based bound (eq. 3) and is discipline-independent. The lower bound is the multiclass throughput guarantee of Theorem 2 (eq. 15): X_c >= N_c / (Z_c + sum_k V_kc (S_kc + d_kc+)), where d_kc+ depends on the discipline at station k – FIFO (Theorem 1 / Lemma 1), processor sharing (Lemma 2), preemptive priority (Lemma 3), non-preemptive priority (Lemmas 4-5). The coupled inequalities are resolved by the interval- narrowing fixed point reproducing the BNR-Prolog robust box bounds; for a single FIFO class it reduces to the Muntz-Wong bounds. Only queueing stations are passed; Z aggregates the pure-delay stations.

Parameters:

V (numpy.ndarray) – (K, C) mean visits of class c at queueing station k.
S (numpy.ndarray) – (K, C) mean service demand per visit of class c at station k.
N (numpy.ndarray) – (C,) population of class c.
Z (numpy.ndarray) – (C,) think time of class c (default zeros).
sched (numpy.ndarray) – (K,) discipline code per station (0=FIFO, 1=PS, 2=non-preemptive priority, 3=preemptive priority, 4=ABA full-contention discipline-independent); default all FIFO.
prio (numpy.ndarray) – (C,) class priority, lower value = higher priority; default equal.

Returns:

(C,) lower bound on class throughput (Theorem 2). Xup: (C,) upper bound on class throughput (eqs. 2-3). Wlo: (K, C) per-visit residence time consistent with the lower bound.

Return type:

Xlo

pfqn_pbh(L, N, Z=0.0, level=1)[source]

Performance Bound Hierarchy (Eager-Sevcik 1983), single-class.

Returns (Xlo, Xhi, Qlo, Qhi). Level-level throughput/queue bounds; level 1 (Z=0) equals the BJB optimistic bound, and the bracket tightens to exact MVA as level -> N.

pfqn_pbk(L, N, Z=0.0, k=1)[source]: Iterative PB(k) proportional bounds (Eager-Sevcik / CMS08). Backed by the PBH recursion at level k. Returns (Xlo, Xhi).

pfqn_bjbk(L, N, Z=0.0, k=1)[source]: Iterative BJB(k) balanced job bounds (CMS08). BJB(1) recovers the noniterative balanced job bound. Returns (Xlo, Xhi).

pfqn_cbh(L, N, Z=0.0, level=2)[source]

Convolutional Bound Hierarchy (Dowdy et al. 1984), single-class.

level exactly-convolved servers (1..M); the bracket tightens monotonically and equals exact at level M. Returns (Xlo, Xhi).

pfqn_mcub(L, N, Z=None)[source]

Multiclass Composite Upper Bound (Kerola 1986). L is M x R.

Returns (Xub, Xlb): Xub the per-class composite UPPER bound (eqs 13-16), Xlb the per-class multiclass Balanced Job Bounds LOWER bound (eq 10) that seeds it. Both are 1 x R arrays.

pfqn_ssd(L, N, Z=0.0, nservers=None)[source]: Server-Station Disaggregation bounds (Suri-Dallery 1986, Thm 5), single-class multiserver. Returns (Xlo, Xhi).

pfqn_sib(L, N, Z=0.0, level=3)[source]

Successively Improving Bounds (Srinivasan 1985), single-class, Z=0 only.

Returns (Xlo, Xhi, Wlo, Whi). Raises ValueError for Z>0 (delay needs the Section-3.2 demand substitution, not yet implemented).

pfqn_ldbcmp(L, N, Z=0.0, c=None, tol=1e-10)[source]

Anselmi-Cremonesi (2008) lower throughput bound for closed single-class BCMP networks with load-dependent stations. Returns (Xlo, Rhi, Qhat).

c[i]=0 marks a fixed-rate (LI) station; c[i]>0 a Heffes LD station with open queue (c[i]+1)*rho/(1-rho). Bottleneck assumed fixed-rate. NaN if N < Qhat.

pfqn_marie(L, N, Z=None, scv=None, tol=1e-8, maxiter=1000, nservers=None)[source]

Marie’s method (single-class). L is M x 1 demands, N scalar population, Z scalar think time, scv per-station SCV (M,). Returns (X, Q, U, C, it, mu): X scalar chain throughput, Q/U/C per-station (M,), it iterations, mu the converged LD multiplier matrix (M x N).

Multiclass (L with >1 column) raises NotImplementedError.

pfqn_le(L, N, Z=None)[source]

Logistic Expansion (LE) asymptotic approximation for normalizing constant.

Provides an asymptotic estimate of the normalizing constant for closed product-form queueing networks. Useful for large populations where exact methods become computationally expensive.

Parameters:

L (numpy.ndarray) – Service demand matrix (M x R).
N (numpy.ndarray) – Population vector (R,).
Z (numpy.ndarray | None) – Think time vector (R,). Optional.

Returns:

Gn: Estimated normalizing constant. lGn: Logarithm of normalizing constant.

Return type:

Tuple of (Gn, lGn)

Reference:: G. Casale. “Accelerating performance inference over closed systems by asymptotic methods.” ACM SIGMETRICS 2017.

pfqn_cub(L, N, Z=None, order=None, atol=1e-8)[source]

Cubature method for normalizing constant using Grundmann-Moeller rules.

Uses numerical integration over simplices to compute the normalizing constant exactly (for sufficient order) or approximately.

Parameters:

L (numpy.ndarray) – Service demand matrix (M x R).
N (numpy.ndarray) – Population vector (R,).
Z (numpy.ndarray | None) – Think time vector (R,). Optional.
order (int | None) – Degree of cubature rule (default: ceil((sum(N)-1)/2)).
atol (float) – Absolute tolerance (default: 1e-8).

Returns:

Gn: Estimated normalizing constant. lGn: Logarithm of normalizing constant.

Return type:

Tuple of (Gn, lGn)

Reference:: G. Casale. “Accelerating performance inference over closed systems by asymptotic methods.” ACM SIGMETRICS 2017.

pfqn_mci(D, N, Z=None, I=100000, variant='imci')[source]

Monte Carlo Integration (MCI) for normalizing constant estimation.

Provides a Monte Carlo estimate of the normalizing constant for closed product-form queueing networks.

Parameters:

D (numpy.ndarray) – Service demand matrix (M x R).
N (numpy.ndarray) – Population vector (R,).
Z (numpy.ndarray | None) – Think time vector (R,). Optional, defaults to zeros.
I (int) – Number of samples (default: 100000).
variant (str) – MCI variant - ‘mci’, ‘imci’ (improved), or ‘rm’ (repairman). Default: ‘imci’.

Returns:

G: Estimated normalizing constant. lG: Logarithm of normalizing constant. lZ: Individual random sample log values.

Return type:

Tuple of (G, lG, lZ)

Reference:: Implementation based on MonteQueue methodology.

pfqn_grnmol(L, N)[source]

Normalizing constant using Grundmann-Moeller quadrature.

Computes the normalizing constant for closed product-form queueing networks using Grundmann-Moeller cubature rules on simplices.

This is an exact method that uses polynomial quadrature to compute the normalizing constant integral representation.

Parameters:

L (numpy.ndarray) – Service demand matrix (M x R).
N (numpy.ndarray) – Population vector (R,).

Returns:

G: Normalizing constant. lG: Logarithm of normalizing constant.

Return type:

Tuple of (G, lG)

Reference:: Grundmann, A. and Moller, H.M. “Invariant Integration Formulas for the N-Simplex by Combinatorial Methods”, SIAM J Numer. Anal. 15 (1978), pp. 282-290.

pfqn_le_fpi(L, N)[source]

Fixed-point iteration to find mode location (no think time).

Public wrapper for the internal _pfqn_le_fpi function.

Parameters:

L (numpy.ndarray) – Service demand matrix (M x R).
N (numpy.ndarray) – Population vector (R,).

Returns:

Mode location vector u (M,).

Return type:

numpy.ndarray

pfqn_le_fpiZ(L, N, Z)[source]

Fixed-point iteration to find mode location (with think time).

Public wrapper for the internal _pfqn_le_fpiZ function.

Parameters:

L (numpy.ndarray) – Service demand matrix (M x R).
N (numpy.ndarray) – Population vector (R,).
Z (numpy.ndarray) – Think time vector (R,).

Returns:

u: Mode location vector (M,). v: Scale factor.

Return type:

Tuple of (u, v)

pfqn_le_hessian(L, N, u)[source]

Compute Hessian matrix (no think time case).

Public wrapper for the internal _pfqn_le_hessian function.

Parameters:

L (numpy.ndarray) – Service demand matrix (M x R).
N (numpy.ndarray) – Population vector (R,).
u (numpy.ndarray) – Mode location vector (M,).

Returns:

Hessian matrix (M-1 x M-1).

Return type:

numpy.ndarray

pfqn_le_hessianZ(L, N, Z, u, v)[source]

Compute Hessian matrix (with think time case).

Public wrapper for the internal _pfqn_le_hessianZ function.

Parameters:

L (numpy.ndarray) – Service demand matrix (M x R).
N (numpy.ndarray) – Population vector (R,).
Z (numpy.ndarray) – Think time vector (R,).
u (numpy.ndarray) – Mode location vector (M,).
v (float) – Scale factor.

Returns:

Hessian matrix (M x M).

Return type:

numpy.ndarray

pfqn_ncld(L, N, Z, mu, options=None)[source]

Main method to compute normalizing constant of a load-dependent model.

Provides the main entry point for computing normalizing constants in load-dependent queueing networks with automatic method selection and preprocessing.

Parameters:

L (numpy.ndarray) – Service demands at all stations (M x R)
N (numpy.ndarray) – Number of jobs for each class (1 x R)
Z (numpy.ndarray) – Think times for each class (1 x R)
mu (numpy.ndarray) – Load-dependent scalings (M x Ntot)
options (Dict[str, Any] | None) – Solver options with keys: - method: ‘default’, ‘exact’, ‘rd’, ‘comomld’, etc. - tol: Numerical tolerance

Returns:

PfqnNcResult with G (normalizing constant), lG (log), and method used

Return type:

PfqnNcResult

pfqn_ld_is(L, N, Z=None, mu=None, options=None)[source]

Importance-sampling (IS) estimate of the normalizing constant of a closed LOAD-DEPENDENT product-form queueing network. Load-dependent counterpart of pfqn_pas_is / pfqn_oi_is: the same sample-an-ordering estimator, with the order-independent rank rate replaced by the load-dependent capacity.

Identity. Every product-form station’s balance function is the sum, over the orderings q of a given per-class count vector n, of an ordered product of a per-position factor:

F_i(n) = |n|!/prod_r(n_r!) * prod_r L(i,r)^{n_r} / prod_{k=1}^{|n|} mu_i(k)
       = sum_{q: |q|=n} prod_{p=1}^{|n|} L(i,q_p) / mu_i(p)

since the multiset has |n|!/prod_r(n_r!) orderings, each contributing the same ordered product. The delay (infinite-server) node is the special case mu_Z(k)=k, giving F_Z(n)=prod_r Z_r^{n_r}/n_r!; a single-server queue is mu_i(k)=1; a c-server queue is mu_i(k)=min(k,c).

Consequently, with ell = sum(N) and a “cut vector” splitting an ordering c of all ell jobs into S contiguous segments (one per station):

G(N) = sum_{c} sum_{cuts} prod_{m=1}^{S} w_m(seg_m),
w_m(q) = prod_{p=1}^{|q|} L(m,q_p) / mu_m(p)

because summing over the orderings of each segment independently reproduces prod_m F_m(n_m), and each count split is realized exactly once.

Estimator. An ordering c is drawn by placing, at each step, a uniformly random present class; p(c) is the product of the reciprocal branching factors. For the sampled c the inner sum over ALL cut vectors is computed exactly by the dynamic program A_0(0)=1, A_m(k) = sum_{j<=k} A_{m-1}(j) * w_m(c_{j+1..k}), so S(c)=A_S(ell) in O(S*ell^2) time (no cut enumeration). Then G = E_{C~p}[S(C)/p(C)] is unbiased, estimated by the sample mean.

Parameters:

L ((M, R) array) – Per-class service demands at the M queueing stations.
N ((R,) array) – Closed population vector, finite.
Z ((R,) array, optional) – Aggregated think time (delay) demand; None or zeros if none.
mu ((M, ell) array or sequence of callables, optional) – Load-dependent capacities; mu[i][k-1] is the capacity of station i holding k jobs. None for the load-independent case mu(i,k)=1 (see pfqn_is()).
options (dict or options object, optional) – Fields samples (default 1e4) and seed (optional).

Returns:

PfqnNcResult with ``G` the IS estimate` of the normalizing constant and
lG = log(G).

Return type:

PfqnNcResult

Examples

>>> L = np.array([[0.5, 0.3], [0.2, 0.4]]); N = np.array([3, 2]); Z = np.array([1.0, 1.0])
>>> mu = np.array([[1, 2, 2, 2, 2], [1, 1, 1, 1, 1]], dtype=float)
>>> res = pfqn_ld_is(L, N, Z, mu, {'samples': 100000, 'seed': 7})

pfqn_pas_is(N, mu, H=None, options=None)[source]

Importance-sampling estimate of G_C and mean queue lengths of a closed two-station P&S tandem with swap graph H.

G_C = sum_{c in D} sum_{k=0}^{ell} Phi_1(c[:k]) Phi_2(reverse(c[k:])), where D is the set of orderings non-decreasing w.r.t. H and Phi_m(q) = prod_p 1/mu_m(supp(q[:p])). Orderings are drawn from D by placing, at each step, a uniformly random placement-order-minimal present class (auto-normalized IS, notebook generator IS_3); E[xi] = G_C[xi]/G_C[1] reuses the same samples. Taking xi = class-r count in the prefix gives the station-1 mean queue length of class r.

Parameters:

N ((R,) closed population vector (macrostate), finite.)
mu (list of exactly two callables. mu[m](n) is the total OI rank rate of) – station m given the per-class occupancy (count) vector n (depends only on supp(n)); this is the count-based svcRateFun of an OI/PAS node.
H ((R, R) swap-graph adjacency (H[b, a] != 0 forbids a before b). Empty ->) – pure OI (all orderings feasible).
options (dict or options object; fields ``samples` (default 1e4),`) – seed (optional), verbose (default False).

Returns:

(G, lG, Q) – station-1 mean queue lengths and Q[1] = N - Q[0].

Return type:

G is the IS estimate of G_C, lG = log(G), Q is (2, R) with the

pas_placement(H)[source]

Placement-order logic of a P&S / OI network with swap graph H.

An ordering c is feasible iff it is non-decreasing w.r.t. H, i.e. class a never appears before class b whenever H[b, a] != 0 (Comte and Dorsman, 2021); equivalently H[b, a] != 0 means b must precede a. The transitive closure P[i, j] = 1 iff i must precede j makes the full order explicit.

Parameters:: H ((R, R) swap-graph adjacency. Empty/all-zero -> no constraint.)
Returns:: (P, placeable) – placeable(x) returns the array of class indices drawable next given the remaining per-class count vector x (present classes with no remaining predecessor still to be placed).
Return type:: P is the (R, R) precedence closure (or None if H is empty);

pas_swap2order(swap, listRate, N0=None)[source]

Global placement-order DAG H of a closed two-station P&S tandem 1->2->1.

With a non-empty swap graph the ordered chain is reducible; the recurrent communicating class is the set of splits of the orderings that are the linear extensions of a single placement partial order (Comte and Dorsman, 2021). pfqn_pas_is samples those orderings from H, so it needs this GLOBAL order. The order is class-level (multiplicity-independent): enumerate the reachable class from the all-in-queue-1 single-job-per-class state; each reachable state (l1; l2) exposes the full ordering c = l1 + reverse(l2); then H[i, j] = 1 iff i precedes j in every such c (forced precedence).

Parameters:

swap ((R, R) swap graph, or list [G1, G2] of the two per-queue graphs) – (G[a, b] != 0 means class a chases class b). Empty/all-zero -> H = 0.
listRate (list of two callables; listRate[m](c) is the total service rate) – of queue m on the ordered prefix c (0-based). Prunes zero-rate (non- head) completions.
N0 ((R,) minimal probing population; defaults to ones(R).)

Returns:

H

Return type:

(R, R) global placement-order DAG; H[i, j] = 1 iff i must precede j.

pfqn_mvaoi(Z, N, mu, Dli=None, options=None)[source]

Mean-value analysis of a closed product-form OI network.

Mean-value counterpart of pfqn_oi_nc() and the marginal form pfqn_mvaoi_marg(): for a closed product-form network of an aggregated infinite-server (delay) node, any number of load-independent (LI) single-server product-form queues, and any number of order-independent (OI) stations, it returns the same exact per-class throughput and queue-lengths WITHOUT computing any normalizing constant or joint marginal, using only mean quantities. It is the composition-dependent generalization of the Conditional MVA (CMVA) of Casale, “A Note on Stable Flow-Equivalent Aggregation in Closed Networks” (QUESTA 2009), extended to MULTIPLE OI stations by carrying one rate-shift vector s_i per OI station i (row i of the shift matrix S).

Throughout, r and s index job classes; i indexes OI stations; j indexes LI queues. State (S, Nn) is processed by increasing sum(Nn); each OI station keeps its own D^i, rho^i and Q^i recursions driven by the common throughput X^{(S)}(Nn), and the population conservation aggregates every station’s contribution:

Nn_r = X_r Z_r + sum_j Q^{(j)}_r + sum_i Q^{(i)}_r,

with the LI queue term Q^{(j)}_r = X_r D_{j,r}(1 + sum_s Q^{(j)}_s(Nn - e_r)).

Parameters:

Z ((R,) think-time demand vector of the aggregated delay node.)
N ((R,) closed population vector, finite.)
mu (callable or list of callables ``mu_i(n)` returning the OI total service`) – rate of station i for the per-class occupancy (count) vector n. A bare callable is accepted as the single-station shorthand.
Dli ((J, R) per-class demand matrix of the LI single-server queues; None or) – empty when J = 0.
options (accepted for signature parity; unused.)

Returns:

(X, Qoi, Qli, Qdelay, Soi) – (K, R), LI queue-lengths (J, R), delay queue-length (R,) = X * Z, and the per-class mean number of IN-SERVICE jobs at each OI station (K, R). Soi[i, r] = E[sir_r] counts the class-r jobs receiving a strictly positive rank rate (see pfqn_oi_insvc()); the utilization of OI station i is Soi[i, r] / c_i. Unlike X/Qoi/Qli, which are pure mean-value quantities, Soi is a distributional statistic and is obtained from the OI count marginal assembled from the zero-shift throughputs X^{(0)}(k) already cached by the mean-value recursion above (no normalizing constant is formed).

Return type:

per-class throughput (R,), OI queue-lengths

pfqn_mvaoi_marg(D, N, isDelay, mu)[source]

Exact marginal load-dependent MVA for OI networks.

Marginal-distribution counterpart of pfqn_mvaoi(). Carries, for each OI station, its joint count-vector marginal pM_i(n | k) and closes the per-class throughput by population conservation. Handles delay + LI product-form queues + any number of OI stations.

Parameters:

D ((M, R) per-class demand at every station (OI rows ignored).)
N ((R,) closed population vector, finite.)
isDelay ((M,) True for infinite-server (delay) stations.)
mu (length-M list; ``mu[i]` is the OI rate callable` of the count vector n, or) – None for non-OI stations.

Returns:

(XN, QN)

Return type:

per-class throughput (R,) and per-station queue-lengths (M, R).

class PfqnNcResult(G, lG, method='default')[source]

Bases: object

Result of normalizing constant computation.

method: str = 'default'

G: float

lG: float

class PfqnComomrmLdResult(G, lG, prob)[source]

Bases: object

Result of COMOM load-dependent computation.

G: float

lG: float

prob: numpy.ndarray

class PfqnFncResult(mu, c)[source]

Bases: object

Result of functional server scaling computation.

mu: numpy.ndarray

c: numpy.ndarray

pfqn_unique(L, mu=None, gamma=None, tol=1e-14)[source]

Consolidate replicated stations into unique stations with multiplicity.

Identifies stations with identical demand rows L[i,:] and (if present) identical load-dependent rates mu[i,:] or class-dependent rates gamma[i,:]. Returns reduced matrices with only unique stations plus a multiplicity vector.

Parameters:

L (numpy.ndarray) – Service demand matrix (M x R)
mu (numpy.ndarray | None) – Load-dependent rate matrix (M x Ntot), optional - pass None if not used
gamma (numpy.ndarray | None) – Class-dependent service rate matrix (M x R), optional - pass None if not used
tol (float) – Tolerance for floating point comparison (default 1e-14)

Returns:

PfqnUniqueResult containing reduced matrices and mapping information

Return type:

PfqnUniqueResult

pfqn_expand(QN, UN, CN, mapping)[source]

Expand per-station metrics from reduced model to original dimensions.

Expands performance metrics computed on a reduced model (with unique stations) back to the original model dimensions by replicating values according to mapping.

Parameters:

QN (numpy.ndarray) – Queue lengths from reduced model (M’ x R)
UN (numpy.ndarray) – Utilizations from reduced model (M’ x R)
CN (numpy.ndarray) – Cycle times from reduced model (M’ x R)
mapping (numpy.ndarray) – Mapping vector from pfqn_unique (length M), mapping[i] = unique station index

Returns:

Tuple of (QN_full, UN_full, CN_full) in original dimensions (M x R)

Return type:

Tuple[numpy.ndarray, numpy.ndarray, numpy.ndarray]

pfqn_combine_mi(mi, mapping, M_unique)[source]

Combine user-provided multiplicity vector with detected replica multiplicity.

For each unique station j, sums the mi values of all original stations mapping to it.

Parameters:

mi (numpy.ndarray) – User-provided multiplicity vector (1 x M_original or M_original,)
mapping (numpy.ndarray) – Mapping vector from pfqn_unique (length M_original)
M_unique (int) – Number of unique stations

Returns:

Combined multiplicity vector (1 x M_unique)

Return type:

numpy.ndarray

class PfqnUniqueResult(L_unique, mu_unique, gamma_unique, mi, mapping)[source]

Bases: NamedTuple

Result class for pfqn_unique containing all output matrices and mapping information.

L_unique

Reduced demand matrix (M’ x R) with M’ <= M unique stations

Type:: numpy.ndarray

mu_unique

Reduced load-dependent rates (M’ x Ntot), None if mu was empty

Type:: numpy.ndarray | None

gamma_unique

Reduced class-dependent rates (M’ x R), None if gamma was empty

Type:: numpy.ndarray | None

mi

Multiplicity vector (1 x M’), mi[j] = count of stations mapping to unique station j

Type:: numpy.ndarray

mapping

Mapping vector (1 x M), mapping[i] = unique station index for original station i

Type:: numpy.ndarray

Create new instance of PfqnUniqueResult(L_unique, mu_unique, gamma_unique, mi, mapping)

L_unique: numpy.ndarray: Alias for field number 0

gamma_unique: numpy.ndarray | None: Alias for field number 2

mapping: numpy.ndarray: Alias for field number 4

mi: numpy.ndarray: Alias for field number 3

mu_unique: numpy.ndarray | None: Alias for field number 1

pfqn_lldfun(n, lldscaling=None, nservers=None)[source]

AMVA-QD load and queue-dependent scaling function.

Computes the scaling factor for load-dependent queueing stations, accounting for multi-server stations and general load-dependent service rate scaling.

Parameters:

n (numpy.ndarray) – Queue population vector (M,)
lldscaling (numpy.ndarray | None) – Load-dependent scaling matrix (M x Nmax), optional
nservers (numpy.ndarray | None) – Number of servers per station (M,), optional

Returns:

Scaling factor vector (M,)

Return type:

numpy.ndarray

References

Original MATLAB: matlab/src/api/pfqn/pfqn_lldfun.m

pfqn_mu_ms(N, m, c)[source]

Compute load-dependent rates for m identical c-server FCFS stations.

Calculates the effective service rate as a function of the number of jobs in the system for a network of m identical stations, each with c parallel servers.

Parameters:

N (int) – Maximum population
m (int) – Number of identical stations
c (int) – Number of servers per station

Returns:

Load-dependent service rate vector (1 x N)

Return type:

numpy.ndarray

References

Original MATLAB: matlab/src/api/pfqn/pfqn_mu_ms.m

pfqn_nc_sanitize(lam, L, N, Z, atol=1e-8)[source]

Sanitize and preprocess network parameters for NC solvers.

Removes empty/ill-defined classes, rescales demands for numerical stability, and reorders classes by think time.

Parameters:

lam (numpy.ndarray) – Arrival rate vector
L (numpy.ndarray) – Service demand matrix (M x R)
N (numpy.ndarray) – Population vector
Z (numpy.ndarray) – Think time vector
atol (float) – Absolute tolerance for numerical comparisons

Returns:

lambda: Sanitized arrival rates
L: Sanitized service demands (rescaled)
N: Sanitized populations
Z: Sanitized think times (rescaled)
lGremaind: Log normalization factor from removed classes

Return type:

Tuple of (lambda, L, N, Z, lGremaind) where

References

Original MATLAB: matlab/src/api/pfqn/pfqn_nc_sanitize.m

pfqn_cdfun(nvec, cdscaling=None, class_idx=0)[source]

AMVA-QD class-dependence function for queue-dependent scaling.

Returns, for every station i, the reciprocal of the class-dependent scaling

beta_{i,r}(n_i1, …, n_iR)

evaluated at the per-class population vector nvec[i, :], for class r = class_idx.

cdscaling[i] is a callable of the per-class population vector at station i. It may return either

a scalar, i.e. a chain-independent scaling beta_i(n) shared by every class (the common case, and the historical contract), or
an array of length R, i.e. the per-class scalings [beta_{i,1}(n), …, beta_{i,R}(n)], of which element class_idx is taken.

The per-class form expresses Sauer’s chain-dependent service rates mu_{r,i}(n) (Sauer 1983, “Computational Algorithms for State-Dependent Queueing Networks”, eq. (40)), so a single class-dependence mechanism covers both the chain-independent and the chain-specific cases.

An empty (None) entry means station i declares no class dependence and is left at the neutral scaling 1.

Parameters:

nvec (numpy.ndarray) – Population state matrix (M x R) or vector (M,)
cdscaling (List | None) – List of class-dependence callables, one per station
class_idx (int) – 0-based class index selecting beta_{i,r} (default: 0)

Returns:

Scaling factor vector (M,)

Return type:

numpy.ndarray

References

Original MATLAB: matlab/src/api/pfqn/pfqn_cdfun.m

factln(n)[source]

Compute log(n!) using log-gamma function.

Parameters:: n (float) – Non-negative number
Returns:: log(n!) = log(Gamma(n+1))
Return type:: float

factln_vec(arr)[source]

Compute log(n!) element-wise for an array.

Parameters:: arr (numpy.ndarray) – Array of non-negative numbers
Returns:: Array of log(n!) values
Return type:: numpy.ndarray

softmin(a, b, alpha=20.0)[source]

Compute a smooth approximation to min(a, b) using weighted average.

Matches MATLAB formula: (x*exp(-alpha*x) + y*exp(-alpha*y)) / (exp(-alpha*x) + exp(-alpha*y))

Parameters:

a (float) – First value
b (float) – Second value
alpha (float) – Smoothing parameter (larger = sharper approximation)

Returns:

Smooth approximation of min(a, b)

Return type:

float

oner(n, s)[source]

Return a copy of n with position s reduced by 1.

Parameters:

n (numpy.ndarray) – Population vector
s (int) – Index to decrement (0-based)

Returns:

Copy of n with n[s] -= 1

Return type:

numpy.ndarray

multichoose(r, n)[source]

Generate all combinations with repetition.

Returns all ways to choose n items from r categories with repetition, where the result is a matrix with each row being a combination.

Parameters:

r (int) – Number of categories
n (int) – Number of items to choose

Returns:

Matrix (C x r) where C = C(n+r-1, r-1) is the number of combinations

Return type:

numpy.ndarray

matchrow(matrix, row)[source]

Find the index of a row in a matrix.

Parameters:

matrix (numpy.ndarray) – 2D array to search in
row (numpy.ndarray) – 1D array to find

Returns:

1-based index of the matching row, or 0 if not found

Return type:

int

pfqn_comom(L, N, Z=None, atol=1e-8)[source]

CoMoM algorithm for computing the normalizing constant.

Implements the Composite Method of Moments algorithm for computing normalizing constants in closed product-form queueing networks.

Parameters:

L (numpy.ndarray) – Service demand matrix (M x R)
N (numpy.ndarray) – Population vector (R,)
Z (numpy.ndarray) – Think time vector (R,), default zeros
atol (float) – Absolute tolerance

Returns:

Logarithm of the normalizing constant

Return type:

lG

References

Original MATLAB: matlab/src/api/pfqn/pfqn_comom.m

pfqn_comomrm(L, N, Z, m=1, atol=1e-8)[source]

CoMoM for finite repairman model.

Computes the normalizing constant for a closed network with a single queueing station and delay stations (repairman model).

Parameters:

L (numpy.ndarray) – Service demand matrix (1 x R) - single station
N (numpy.ndarray) – Population vector (R,)
Z (numpy.ndarray) – Think time vector (R,)
m (int) – Replication factor (default: 1)
atol (float) – Absolute tolerance

Returns:

ComomResult with lG (log normalizing constant) and lGbasis

Return type:

ComomResult

References

Original MATLAB: matlab/src/api/pfqn/pfqn_comomrm.m

pfqn_comomrm_orig(L, N, Z, m=1, atol=1e-8)[source]

Original CoMoM implementation for repairman model.

This is the original implementation of CoMoM without optimizations. Kept for reference and validation.

Parameters:

L (numpy.ndarray) – Service demand matrix (1 x R)
N (numpy.ndarray) – Population vector (R,)
Z (numpy.ndarray) – Think time vector (R,)
m (int) – Replication factor (default: 1)
atol (float) – Absolute tolerance

Returns:

ComomResult with lG and lGbasis

Return type:

ComomResult

References

Original MATLAB: matlab/src/api/pfqn/pfqn_comomrm_orig.m

pfqn_comomrm_ms(L, N, Z, m, c, atol=1e-8)[source]

CoMoM for multi-server repairman model.

Computes the normalizing constant for a repairman model where the queueing station has multiple servers.

Parameters:

L (numpy.ndarray) – Service demand matrix (1 x R)
N (numpy.ndarray) – Population vector (R,)
Z (numpy.ndarray) – Think time vector (R,)
m (int) – Number of identical stations
c (int) – Number of servers per station
atol (float) – Absolute tolerance

Returns:

ComomResult with lG and lGbasis

Return type:

ComomResult

References

Original MATLAB: matlab/src/api/pfqn/pfqn_comomrm_ms.m

pfqn_procomom(L, N, Z=None, atol=1e-14)[source]

ProCoMoM algorithm for computing marginal queue-length probabilities.

Computes the marginal queue-length probability distribution at each station using the Probabilistic Class-Oriented Method of Moments. Uses matrix recursion with SVD/QR decomposition for numerical stability, with automatic perturbation on rank deficiency.

Parameters:

L (numpy.ndarray) – Service demand matrix (M x R).
N (numpy.ndarray) – Population vector (R,).
Z (numpy.ndarray) – Think time vector (R,), default zeros.
atol (float) – Absolute numerical tolerance (default 1e-14).

Returns:

Marginal probability matrix (M x sumN+1).: Pr[k, j] = P(n_k = j) for station k, queue length j.

Q: Mean queue length vector (M,).

Return type:

Pr

References

Original MATLAB: matlab/src/api/pfqn/pfqn_procomom.m

pfqn_procomom2(L, N, Z=None, atol=1e-8)[source]

Projected CoMoM method for normalizing constant.

Uses a projection-based approach to compute the normalizing constant, which can be more efficient for certain model structures.

Parameters:

L (numpy.ndarray) – Service demand matrix (M x R)
N (numpy.ndarray) – Population vector (R,)
Z (numpy.ndarray) – Think time vector (R,), default zeros
atol (float) – Absolute tolerance

Returns:

Logarithm of the normalizing constant

Return type:

lG

References

Original MATLAB: matlab/src/api/pfqn/pfqn_procomom2.m

class ComomResult(lG, lGbasis=None)[source]

Bases: NamedTuple

Result of CoMoM normalizing constant computation.

Create new instance of ComomResult(lG, lGbasis)

lG: float: Alias for field number 0

lGbasis: numpy.ndarray | None: Alias for field number 1

pfqn_mmint2(L, N, Z, m=1)[source]

McKenna-Mitra integral form using scipy.integrate.

Computes the normalizing constant using numerical integration of the McKenna-Mitra integral representation.

Parameters:

L (numpy.ndarray) – Service demand vector (R,)
N (numpy.ndarray) – Population vector (R,)
Z (numpy.ndarray) – Think time vector (R,)
m (int) – Replication factor (default: 1)

Returns:

G: Normalizing constant
lG: Logarithm of normalizing constant

Return type:

Tuple of (G, lG)

References

Original MATLAB: matlab/src/api/pfqn/pfqn_mmint2.m

pfqn_mmint2_gausslegendre(L, N, Z, m=1)[source]

McKenna-Mitra integral with Gauss-Legendre quadrature.

Uses Gauss-Legendre quadrature for improved accuracy in computing the McKenna-Mitra integral.

Parameters:

L (numpy.ndarray) – Service demand vector (R,)
N (numpy.ndarray) – Population vector (R,)
Z (numpy.ndarray) – Think time vector (R,)
m (int) – Replication factor (default: 1)

Returns:

G: Normalizing constant
lG: Logarithm of normalizing constant

Return type:

Tuple of (G, lG)

References

Original MATLAB: matlab/src/api/pfqn/pfqn_mmint2_gausslegendre.m

pfqn_mmint2_gausslaguerre(L, N, Z, m=1)[source]

McKenna-Mitra integral with Gauss-Laguerre quadrature.

Uses Gauss-Laguerre quadrature which is naturally suited for integrals with exponential decay.

Parameters:

L (numpy.ndarray) – Service demand vector (R,)
N (numpy.ndarray) – Population vector (R,)
Z (numpy.ndarray) – Think time vector (R,)
m (int) – Replication factor (default: 1)

Returns:

G: Normalizing constant
lG: Logarithm of normalizing constant

Return type:

Tuple of (G, lG)

References

Original MATLAB: matlab/src/api/pfqn/pfqn_mmint2_gausslaguerre.m

pfqn_mmsample2(L, N, Z, m=1, n_samples=10000)[source]

Monte Carlo sampling approximation for normalizing constant.

Uses Monte Carlo sampling to approximate the McKenna-Mitra integral. Useful for very large populations where quadrature becomes expensive.

Parameters:

L (numpy.ndarray) – Service demand vector (R,)
N (numpy.ndarray) – Population vector (R,)
Z (numpy.ndarray) – Think time vector (R,)
m (int) – Replication factor (default: 1)
n_samples (int) – Number of Monte Carlo samples (default: 10000)

Returns:

G: Normalizing constant (approximate)
lG: Logarithm of normalizing constant

Return type:

Tuple of (G, lG)

References

Original MATLAB: matlab/src/api/pfqn/pfqn_mmsample2.m

logsumexp(x)[source]

Compute log(sum(exp(x))) in a numerically stable way.

Parameters:: x (numpy.ndarray) – Array of values
Returns:: log(sum(exp(x)))
Return type:: float

pfqn_schmidt(D, N, S, sched, v=None)[source]

Schmidt’s exact MVA for networks with general scheduling disciplines.

Implements Schmidt’s exact Mean Value Analysis algorithm for product-form queueing networks with PS, FCFS, or INF scheduling disciplines, including support for multi-server stations.

Parameters:

D (numpy.ndarray) – Service demand matrix (M x R)
N (numpy.ndarray) – Population vector (R,)
S (numpy.ndarray) – Number of servers per station (M,) or (M x R)
sched (numpy.ndarray) – Scheduling discipline per station (M,) - SchedStrategy values
v (numpy.ndarray | None) – Visit ratio matrix (M x R), optional (default: ones)

Returns:

XN: System throughput (M x R) - same per station for closed networks QN: Mean queue lengths (M, R) UN: Utilization (M, R), per station-class, D*X/nservers CN: Cycle times / response times (M, R)

Return type:

Tuple of (XN, QN, UN, CN)

References

Original MATLAB: matlab/src/api/pfqn/pfqn_schmidt.m

pfqn_schmidt_ext(D, N, S, sched, v=None)[source]

Extended Schmidt MVA algorithm with queue-aware alpha corrections.

A queue-aware version of the Schmidt algorithm that precomputes alpha values for improved accuracy in networks with class-dependent FCFS scheduling.

Reference:: R. Schmidt, “An approximate MVA algorithm for exponential, class-dependent multiple server stations,” Performance Evaluation, vol. 29, no. 4, pp. 245-254, 1997.

Parameters:

D (numpy.ndarray) – Service demand matrix (M x R)
N (numpy.ndarray) – Population vector (R,)
S (numpy.ndarray) – Number of servers per station (M,) or (M x R)
sched (numpy.ndarray) – Scheduling discipline per station (M,) - SchedStrategy values
v (numpy.ndarray | None) – Visit ratio matrix (M x R), optional (default: ones)

Returns:

XN: System throughput (M x R) QN: Mean queue lengths (M, R) UN: Utilization (M, R) CN: Cycle times / response times (M, R)

Return type:

Tuple of (XN, QN, UN, CN)

class SchmidtResult(XN, QN, UN, CN)[source]

Bases: object

Result from Schmidt’s exact MVA.

XN: numpy.ndarray

QN: numpy.ndarray

UN: numpy.ndarray

CN: numpy.ndarray

pprod(N, Nmax=None)[source]

Generate population vectors for MVA recursion.

When called with one argument, initializes to zero vector. When called with two arguments, increments to next population vector. Returns -1 when iteration is complete.

Parameters:

N (numpy.ndarray) – Current population vector or max population
Nmax (numpy.ndarray | None) – Maximum population per class (if incrementing)

Returns:

Next population vector, or array of -1 if done

Return type:

numpy.ndarray

hashpop(nvec, Nc, C, prods)[source]

Hash population vector to linear index.

Parameters:

nvec (numpy.ndarray) – Population vector
Nc (numpy.ndarray) – Maximum population per class
C (int) – Number of classes
prods (numpy.ndarray) – Precomputed products for hashing

Returns:

Linear index (1-based for MATLAB compatibility)

Return type:

int

pfqn_dac(L, N, Z=None, mu=None)[source]

DAC (Distribution Analysis by Chain) method for joint queue-length distributions.

Computes the joint queue-length distribution of a closed product-form network by a chain-by-chain recursion over a related network in which every chain holds a single customer, a transformation that leaves the aggregate queue-length distribution unchanged. Given the distribution of a network with k-1 such chains, adding one customer of a chain with demands r gives

c_j = sum_{n=1..k} (n/mu_j(n)) * P_j^{k-1}(n-1) lambda_k = 1 / sum_j r_j c_j P^k(n) = lambda_k * sum_j r_j (n_j/mu_j(n_j)) * P^{k-1}(n-e_j)

where lambda_k is the throughput of the customer being added and 1/c_j is the throughput of a chain visiting center j only. The recursion conserves probability mass by construction, hence it is numerically stable.

Parameters:

L (numpy.ndarray) – Service demand matrix (M x R)
N (numpy.ndarray) – Population vector (R,)
Z (numpy.ndarray | None) – Think time vector (R,), optional (default: zeros). If sum(Z)>0 an extra infinite-server station is appended, so that states has M+1 columns and its last column holds the think-station population.
mu (numpy.ndarray | None) – Load-dependent rate matrix (M x Nt), Nt=sum(N), optional (default: ones, i.e. single-server fixed rate). Use mu[j,:]=1..Nt for infinite server and mu[j,n]=min(n+1,c) for a c-server station.

Returns:

Pjoint: Probability of the aggregate state in the corresponding row: of states, of length nchoosek(Nt+J-1, J-1)

states: Aggregate states (S x J), states[s,j] = jobs at center j XN: Throughput of chain r (R,) QN: Mean number of chain-r customers at station j (M x R) UN: Utilization of station j (M,), i.e. 1-P_j(0) CN: Cycle time of chain r, exclusive of think time (R,) pi: Marginal probabilities, pi[j,n] = P(n jobs at station j),

shape (M, Nt+1)

Return type:

Tuple of (Pjoint, states, XN, QN, UN, CN, pi)

References

E. de Souza e Silva, “Distribution Analysis of Product Form Queueing Networks”, UCLA Computer Science Department, CSD-870023, April 1987.

pfqn_recal(L, N, Z=None, m0=None)[source]

RECAL (REcursive CALculation) method for normalizing constant.

Computes the normalizing constant G(N) using the RECAL recursive method, which is efficient for networks with moderate population sizes.

Parameters:

L (numpy.ndarray) – Service demand matrix (M x R)
N (numpy.ndarray) – Population vector (R,)
Z (numpy.ndarray | None) – Think time vector (R,), optional (default: zeros)
m0 (numpy.ndarray | None) – Initial multiplicity vector (M,), optional (default: ones)

Returns:

G: Normalizing constant lG: Logarithm of normalizing constant

Return type:

Tuple of (G, lG)

References

Original MATLAB: matlab/src/api/pfqn/pfqn_recal.m Conway, A.E. and Georganas, N.D. “RECAL - A New Efficient Algorithm for the Exact Analysis of Multiple-Chain Closed Queueing Networks”, JACM, 1986.

pfqn_mvac(L, N, Z=None)[source]

MVAC (Mean Value Analysis by Chain) for closed product-form networks.

Exact mean value analysis of a closed multichain product-form queueing network composed of single-server fixed-rate (SSFR) queues and infinite-server (IS) centers. Unlike the classic MVA recursion of pfqn_mva, which recurs on the population vector and costs O(prod(N+1)), MVAC recurs on the chains: each class is reduced to single-customer chains and the removed chains are replaced by self-looping single-customer (SCSL) chains pinned at a service center, so the multiplicity vector v = (v_1,…,v_J), with v_j the number of SCSL chains at center j, indexes the recursion in place of the population vector. MVAC is thus attractive for networks with few centers and many chains, and is the mean-value counterpart of the RECAL normalizing-constant recursion of pfqn_recal. Since no normalizing constant is formed, MVAC does not suffer the floating-point underflow/overflow that complicates RECAL and convolution.

With j, i indexing centers (j = 1,…,J1 SSFR and j = J1+1,…,J IS), k, l indexing single-customer chains, a_jk the relative utilization (service demand) of chain k at center j, a_k = sum_j a_jk, K the total number of single-customer chains, and I_k = {v : sum_j v_j = K - k}, the recursion of Section II of the paper reads

lambda^k_k(v) = 1 / (a_k + sum_{j=1}^{J1} a_jk (L^{k-1}_j(v) + v_j)) (10) L^k_{jk}(v) = lambda^k_k(v) a_jk (1 + L^{k-1}_j(v) + v_j), j SSFR (9a) L^k_{jk}(v) = lambda^k_k(v) a_jk, j IS (9b) L^k_i(v) = sum_j L^k_{jk}(v) L^{k-1}_i(v + 1_j) + L^k_{ik}(v) (7) L^k_{il}(v) = sum_j L^k_{jk}(v) L^{k-1}_{il}(v + 1_j), l = 1..k-1 (6)

with L^0_j(v) = 0, where L^k_j(v) is the mean number of customers at center j (SCSL customers excluded), L^k_{jl}(v) the mean number of chain-l customers at center j, and lambda^k_k(v) the throughput of chain k, all for the network with normalizing constant G_k(v). Equation (10) is the arrival-theorem closure obtained by summing (9a)-(9b) over all centers, since chain k holds a single customer. The measures of the original network are read off at k = K and v = 0.

Part 1 of the basic step evaluates (10), (9) and (7) and yields the measures of chain K; part 2 evaluates (6) and yields those of the chains that visit at least one IS center, whose throughput then follows from Little’s law at that center. Chains visiting only SSFR centers require a re-execution of part 1 with their label interchanged with K, which is cheap because the levels below the interchanged label are unaffected and are reused. Classes with N_r > 1, and classes with identical demand columns, collapse into a single subset of identical single-customer chains: only one representative per subset is analyzed and its per-chain measures are scaled by the class population, so the cost depends on the number D of distinct chains, not K.

Parameters:

L (numpy.ndarray) – Service demand matrix of the SSFR queues (M x R)
N (numpy.ndarray) – Population vector (R,), finite and nonnegative
Z (numpy.ndarray | None) – Service demand matrix of the IS centers, (R,) or (Mz x R), one row per IS center, optional (default: zeros)

Returns:

XN: Per-class throughput at the reference station (R,) QN: Per-class mean queue-length at the SSFR queues (M x R) UN: Per-class utilization, XN[r] * L[i,r] (M x R) CN: Per-class residence time, QN[i,r] / XN[r] (M x R)

Return type:

Tuple of (XN, QN, UN, CN)

References

A. E. Conway, E. de Souza e Silva and S. S. Lavenberg, “Mean Value Analysis by Chain of Product Form Queueing Networks”, IEEE Trans. Computers, 38(3):432-442, 1989. Original MATLAB: matlab/src/api/pfqn/pfqn_mvac.m

pfqn_mvacld(L, N, Z=None, mu=None)[source]

MVAC for closed product-form networks with queue-length dependent centers.

Exact mean value analysis by chain of a closed multichain product-form queueing network that may contain queue-length dependent (QLD) service centers. This is the Section V extension of Conway, de Souza e Silva and Lavenberg (1989); pfqn_mvac implements Sections II-IV, which cover single-server fixed-rate (SSFR) and infinite-server (IS) centers only.

Where pfqn_mvac propagates the MEAN queue-lengths L^k_j(v) through eq. (7) and closes the recursion with the arrival-theorem identity (10), the QLD extension propagates the MARGINAL queue-length DISTRIBUTIONS P^k_j(n,v) instead. That is forced by load dependence – the rate seen by a job depends on the whole occupancy, so a mean no longer suffices – but it also SIMPLIFIES the recursion: eq. (21)-(25) read level k-1 only at the shifted vectors v + 1_i, so the basic step sweeps v in I_k alone, where pfqn_mvac must sweep the larger I_k u … u I_K. The marginals come almost for free and are returned as a first-class output.

Notation follows pfqn_mvac: j, i index centers (j = 1,…,J1 the QLD centers of L, j = J1+1,…,J the IS centers of Z), k, l index the single-customer chains, a_jk = theta_jk T_jk is the demand of chain k at center j, K = sum(N) and I_k = {v : sum_j v_j = K - k} with v_j the number of self-looping single-customer (SCSL) chains pinned at center j. P^k_j(n,v) is the probability of n customers at center j – EXCLUDING the v_j SCSL customers there – in the network with normalizing constant G_k(v). Writing tau_k(v,i) for the throughput of an SCSL chain that replaces chain k at center i:

tau_k(v,i) = T_ik^-1 sum_{n=0}^{k-1} P^{k-1}_i(n,v+1_i)

mu_i(n+v_i+1)/(n+v_i+1) (21)

tau_k(v,i) = T_ik^-1, i IS (22) L^k_{jk}(v) = theta_jk tau_k(v,j)^-1 / sum_m theta_mk tau_k(v,m)^-1 (23) lambda^k_k(v) = tau_k(v,j(k)) L^k_{j(k)k}(v) (24) P^k_j(n,v) = L^k_{jk}(v) P^{k-1}_j(n-1,v+1_j)

sum_{m != j} L^k_{mk}(v) P^{k-1}_j(n,v+1_m) (25)

with P^0_j(0,v) = 1 and P^{k-1}_j(n,.) = 0 for n < 0 or n > k-1. Eq. (21) is just “the mean rate at which an SCSL chain is served”: given n other customers the processor-sharing rate share is mu_i(n+v_i+1)/(n+v_i+1), averaged over the distribution of those others. The queueing discipline may be assumed PS with no loss of generality, since product-form measures do not depend on it.

This implementation writes (23)-(24) in the reference-station-free form

c_i(k,v) = sum_{n=0}^{k-1} P^{k-1}_i(n,v+1_i)

mu_i(n+v_i+1)/(n+v_i+1)

L^k_{jk}(v) = (a_jk/c_j) / sum_m (a_mk/c_m) lambda^k_k(v) = 1 / sum_m (a_mk/c_m)

which follows from theta_jk tau_k(v,j)^-1 = a_jk/c_j and theta_{j(k)k} = 1, so only the demands a_jk are needed and the visit ratios never appear separately. For an IS center c_i = 1 identically, which is exactly (22). Eq. (25) is self-normalizing, sum_n P^k_j(n,v) = sum_m L^k_{mk}(v) = 1, so no normalizing constant is formed and the recursion involves only positive quantities: unlike the classic load-dependent MVA of pfqn_mvald it cannot produce negative probabilities and needs no stabilization.

Parts 2 and 3 are unchanged from pfqn_mvac, since eq. (6) holds verbatim in the presence of QLD centers: part 2 resolves the chains that visit at least one IS center, and the chains that visit no IS center are resolved by re-executing part 1 with their label interchanged with K.

Parameters:

L (numpy.ndarray) – Service demand matrix of the QLD centers (M x R)
N (numpy.ndarray) – Population vector (R,), finite and nonnegative
Z (numpy.ndarray | None) – Demand matrix of the IS centers, (R,) or (Mz x R), one row per center, optional (default: zeros)
mu (numpy.ndarray | None) – Load-dependent rates (M x Nt) with Nt >= sum(N); mu[j,n-1] is the total service rate of center j with n jobs present. mu[j,:] = 1 is a single-server fixed-rate queue, mu[j,n-1] = min(n,c) a c-server queue, mu[j,n-1] = n an infinite server. Optional (default: ones, i.e. all centers SSFR, in which case results agree with pfqn_mvac)

Returns:

XN: Per-class throughput at the reference station (R,) QN: Per-class mean queue-length at the QLD centers (M x R) UN: Utilization of each center (M,), 1 - P_j(0). PER-STATION, not

per-class, as in pfqn_mvald and pfqn_dac: for a load-dependent center the per-class product XN[r]*L[j,r] of pfqn_mvac is NOT the utilization

CN: Per-class cycle time exclusive of think time (R,),: N[r]/XN[r]-Z[r], as in pfqn_mvald and pfqn_dac. NOT the (M x R) per-station residence time of pfqn_mvac: the whole LD family reports a cycle time here
pij: Marginal queue-length probabilities (M x (sum(N)+1)),: pij[j,n] = P(n jobs at center j)

Return type:

Tuple of (XN, QN, UN, CN, pij)

References

A. E. Conway, E. de Souza e Silva and S. S. Lavenberg, “Mean Value Analysis by Chain of Product Form Queueing Networks”, IEEE Trans. Computers, 38(3):432-442, 1989, Section V. Original MATLAB: matlab/src/api/pfqn/pfqn_mvacld.m

pfqn_mvaldmx(lam, D, N, Z, mu=None, S=None)[source]

Load-dependent MVA for mixed open/closed networks with limited load dependence.

Implements the MVALDMX algorithm for analyzing mixed queueing networks with load-dependent service rates using limited load dependence.

Parameters:

lam (numpy.ndarray) – Arrival rate vector (R,) - non-zero for open classes
D (numpy.ndarray) – Service demand matrix (M x R)
N (numpy.ndarray) – Population vector (R,) - inf for open classes
Z (numpy.ndarray) – Think time vector (R,)
mu (numpy.ndarray | None) – Load-dependent rate matrix (M x Nt), optional
S (numpy.ndarray | None) – Number of servers per station (M,), optional

Returns:

XN: System throughput (R,) QN: Mean queue lengths (M, R) UN: Utilization (M, R) CN: Cycle times (M, R) lGN: Logarithm of normalizing constant Pc: Marginal queue-length probabilities (M, 1+Ntot, prod(1+Nc))

Return type:

Tuple of (XN, QN, UN, CN, lGN, Pc)

References

Original MATLAB: matlab/src/api/pfqn/pfqn_mvaldmx.m

pfqn_ldmx_ec(lam, D, mu)[source]

Compute effective capacity terms for MVALDMX solver.

Calculates the effective capacity E, E’, and EC terms needed for load-dependent MVA with limited load dependence.

Parameters:

lam (numpy.ndarray) – Arrival rate vector (R,)
D (numpy.ndarray) – Service demand matrix (M x R)
mu (numpy.ndarray) – Load-dependent rate matrix (M x Nt)

Returns:

EC: Effective capacity matrix (M x Nt) E: E-function values (M x (1+Nt)) Eprime: E-prime function values (M x (1+Nt)) Lo: Open class load vector (M,)

Return type:

Tuple of (EC, E, Eprime, Lo)

References

Original MATLAB: matlab/src/api/pfqn/pfqn_ldmx_ec.m

pfqn_mvaldms(lam, D, N, Z, S)[source]

Load-dependent MVA for multiserver mixed networks.

Wrapper for pfqn_mvaldmx that adjusts utilizations to account for multi-server stations.

Parameters:

lam (numpy.ndarray) – Arrival rate vector (R,)
D (numpy.ndarray) – Service demand matrix (M x R)
N (numpy.ndarray) – Population vector (R,)
Z (numpy.ndarray) – Think time vector (R,)
S (numpy.ndarray) – Number of servers per station (M,)

Returns:

XN: System throughput (R,) QN: Mean queue lengths (M, R) UN: Utilization (M, R) - adjusted for multiservers CN: Cycle times (M, R) lGN: Logarithm of normalizing constant

Return type:

Tuple of (XN, QN, UN, CN, lGN)

References

Original MATLAB: matlab/src/api/pfqn/pfqn_mvaldms.m

pfqn_linearizerms(L, N, Z, nservers, type_sched=None, tol=1e-8, maxiter=1000, QN0=None)[source]

Multiserver Linearizer (Krzesinski/Conway/De Souza-Muntz).

Extends the Linearizer algorithm to handle multi-server stations in product-form queueing networks.

Parameters:

L (numpy.ndarray) – Service demand matrix (M x R)
N (numpy.ndarray) – Population vector (R,)
Z (numpy.ndarray) – Think time vector (R,)
nservers (numpy.ndarray) – Number of servers per station (M,)
type_sched (numpy.ndarray | None) – Scheduling strategy per station (M,), optional (default: PS)
tol (float) – Convergence tolerance (default: 1e-8)
maxiter (int) – Maximum iterations (default: 1000)

Returns:

Q: Mean queue lengths (M, R) U: Utilization (M, R) R: Residence times (M, R) C: Cycle times (R,) X: System throughput (R,) totiter: Total iterations performed

Return type:

Tuple of (Q, U, R, C, X, totiter)

References

Original MATLAB: matlab/src/api/pfqn/pfqn_linearizerms.m

pfqn_linearizermx(lambda_arr, L, N, Z, nservers, sched_type, tol=1e-8, maxiter=1000, method='egflin', QN0=None)[source]

Linearizer for mixed open/closed queueing networks.

This function extends the linearizer algorithm to handle networks with both open classes (with external arrivals) and closed classes (with fixed populations).

Parameters:

lambda_arr (numpy.ndarray) – Arrival rate vector (R,). For closed classes, should be 0 or inf.
L (numpy.ndarray) – Service demand matrix (M x R)
N (numpy.ndarray) – Population vector (R,). Inf for open classes, finite for closed.
Z (numpy.ndarray) – Think time vector (R,) or matrix
nservers (numpy.ndarray) – Number of servers per station (M,)
sched_type (List[str]) – Scheduling strategy per station (list of strings)
tol (float) – Convergence tolerance
maxiter (int) – Maximum iterations
method (str) – Linearizer variant (‘lin’, ‘gflin’, ‘egflin’)

Returns:

Mean queue lengths (M x R) UN: Utilization (M x R) WN: Waiting times (M x R) TN: Throughputs (M x R) CN: Cycle times (1 x R) XN: System throughput (R,) totiter: Total iterations

Return type:

QN

References

MATLAB: matlab/src/api/pfqn/pfqn_linearizermx.m

pfqn_conwayms(L, N, Z, nservers, type_sched=None, tol=1e-8, maxiter=1000, QN0=None)[source]

Conway (1989) multiserver Linearizer approximation for FCFS queues.

Implements the algorithm from Conway (1989), “Fast Approximate Solution of Queueing Networks with Multi-Server Chain-Dependent FCFS Queues”.

Parameters:

L (numpy.ndarray) – Service demand matrix (M x R)
N (numpy.ndarray) – Population vector (R,)
Z (numpy.ndarray) – Think time vector (R,)
nservers (numpy.ndarray) – Number of servers per station (M,)
type_sched (numpy.ndarray | None) – Scheduling strategy per station (M,), optional (default: FCFS)
tol (float) – Convergence tolerance (default: 1e-8)
maxiter (int) – Maximum iterations (default: 1000)

Returns:

Q: Mean queue lengths (M, R) U: Utilization (M, R) R: Residence times (M, R) C: Cycle times (R,) X: System throughput (R,) totiter: Total iterations performed

Return type:

Tuple of (Q, U, R, C, X, totiter)

References

Conway, A. E., “Fast Approximate Solution of Queueing Networks with Multi-Server Chain-Dependent FCFS Queues”, Performance Evaluation, Vol. 8, 1989, pp. 141-159.

ljd_linearize(nvec, cutoffs)[source]

Convert per-class population vector to linearized index.

Maps a multi-dimensional population vector to a single linear index for efficient lookups in tabulated scaling tables.

Index formula: idx = 1 + n1 + n2*(N1+1) + n3*(N1+1)*(N2+1) + …

Parameters:

nvec (numpy.ndarray) – Per-class populations [n1, n2, …, nK]
cutoffs (numpy.ndarray) – Per-class cutoffs [N1, N2, …, NK]

Returns:

1-based linearized index

Return type:

int

References

Original MATLAB: matlab/src/api/pfqn/ljd_linearize.m

infradius_h(x, L, N, alpha)[source]

Helper function for infinite radius computation with logistic transformation.

Used in normalizing constant computation via integration methods.

Parameters:

x (numpy.ndarray) – Logistic transformation parameters
L (numpy.ndarray) – Service demand matrix (M x R)
N (numpy.ndarray) – Population vector (R,)
alpha (numpy.ndarray) – Load-dependent rate matrix

Returns:

Evaluated function value for integration

Return type:

numpy.ndarray

References

Original MATLAB: matlab/src/api/pfqn/infradius_h.m

infradius_hnorm(x, L, N, alpha)[source]

Helper function for infinite radius computation with normal CDF (probit) transformation.

Uses normcdf/normpdf transformation instead of logistic (used in infradius_h).

Parameters:

x (numpy.ndarray) – Normal CDF transformation parameters
L (numpy.ndarray) – Service demand matrix (M x R)
N (numpy.ndarray) – Population vector (R,)
alpha (numpy.ndarray) – Load-dependent rate matrix

Returns:

Evaluated function value for integration

Return type:

numpy.ndarray

References

Original MATLAB: matlab/src/api/pfqn/infradius_hnorm.m

pfqn_kt(L, N, Z=None)[source]

Knessl-Tier asymptotic expansion for normalizing constant.

Computes the normalizing constant using Knessl-Tier’s asymptotic expansion, which is particularly accurate for large populations.

Parameters:

L (numpy.ndarray) – Service demand matrix (M x R)
N (numpy.ndarray) – Population vector (R,)
Z (numpy.ndarray | None) – Think time vector (R,), optional (default: zeros)

Returns:

G: Normalizing constant lG: Logarithm of normalizing constant X: System throughput (R,) Q: Mean queue lengths (M, R)

Return type:

Tuple of (G, lG, X, Q)

References

Original MATLAB: matlab/src/api/pfqn/pfqn_kt.m

pfqn_ab_amva(D, N, V, nservers, sched, fcfs_schmidt=False, marginal_prob_method='ab')[source]

Akyildiz-Bolch AMVA method for multi-server BCMP networks.

Parameters:

D (numpy.ndarray) – Service time matrix (M x K)
N (numpy.ndarray) – Population vector (1 x K)
V (numpy.ndarray) – Visit ratio matrix (M x K)
nservers (numpy.ndarray) – Number of servers at each station (M x 1)
sched (numpy.ndarray) – Scheduling strategies for each station (M x 1)
fcfs_schmidt (bool) – Whether to use Schmidt formula for FCFS stations
marginal_prob_method (str) – Method for marginal probability (‘ab’ or ‘scat’)

Returns:

AbAmvaResult containing queue lengths, utilization, residence times, cycle times, throughput, and iteration count.

Return type:

AbAmvaResult

pfqn_ab_core(K, M, population, nservers, sched_type, v, s, maxiter, D_frac, l_in, fcfs_schmidt=False, marginal_prob_method='ab')[source]

Akyildiz-Bolch core method for multi-server BCMP networks.

Public wrapper for the internal core algorithm of the Akyildiz-Bolch linearizer method.

Parameters:

K (int) – Number of classes
M (int) – Number of stations
population (numpy.ndarray) – Population vector (K,)
nservers (numpy.ndarray) – Number of servers at each station (M,)
sched_type (numpy.ndarray) – Scheduling strategies for each station (M,)
v (numpy.ndarray) – Visit ratio matrix (M x K)
s (numpy.ndarray) – Service time matrix (M x K)
maxiter (int) – Maximum iterations
D_frac (numpy.ndarray) – Fractional changes matrix (M x K x K)
l_in (numpy.ndarray) – Initial queue length matrix (M x K)
fcfs_schmidt (bool) – Whether to use Schmidt formula for FCFS stations
marginal_prob_method (str) – Method for marginal probability (‘ab’ or ‘scat’)

Returns:

QN: Queue lengths (M x K) UN: Utilization (M x K) RN: Residence times (M x K) CN: Cycle times (K,) XN: Throughput (K,) totiter: Total iterations

Return type:

Tuple of (QN, UN, RN, CN, XN, totiter)

References

Akyildiz, I.F. and Bolch, G., “Mean Value Analysis Approximation for Multiple Server Queueing Networks”, Performance Evaluation, 1988.

class AbAmvaResult(QN, UN, RN, CN, XN, totiter)[source]

Bases: object

Result from Akyildiz-Bolch AMVA algorithm.

QN: numpy.ndarray

UN: numpy.ndarray

RN: numpy.ndarray

CN: numpy.ndarray

XN: numpy.ndarray

totiter: int

pfqn_rd(L, N, Z, mu=None, options=None)[source]

Reduction Heuristic (RD) method for load-dependent networks.

Computes the logarithm of the normalizing constant using the reduction heuristic method, which handles load-dependent service rates.

Parameters:

L (numpy.ndarray) – Service demand matrix (M x R)
N (numpy.ndarray) – Population vector (1 x R)
Z (numpy.ndarray) – Think time vector (1 x R)
mu (numpy.ndarray | None) – Load-dependent rate matrix (M x sum(N)). If None, assumes load-independent rates (all 1.0).
options (RdOptions | None) – Solver options

Returns:

lGN: Logarithm of normalizing constant Cgamma: Gamma correction factor

Return type:

Tuple of (lGN, Cgamma) where

class RdOptions(tol=1e-06, method='default')[source]

Bases: object

Options for RD algorithm.

method: str = 'default'

tol: float = 1e-06

class RdResult(lGN, Cgamma)[source]

Bases: object

Result from Reduction Heuristic algorithm.

lGN: float

Cgamma: float

pfqn_nrl(L, N, Z=None, alpha=None)[source]

Norlund-Rice Logit (NRL) approximation for normalizing constant.

Computes the logarithm of the normalizing constant using Laplace approximation with logistic transformation.

Parameters:

L (numpy.ndarray) – Service demand matrix (M x R)
N (numpy.ndarray) – Population vector (R,)
Z (numpy.ndarray) – Think time vector (R,) - optional
alpha (numpy.ndarray) – Load-dependent rate matrix (M x Ntot) - optional

Returns:

Logarithm of normalizing constant

Return type:

lG

References

Original MATLAB: matlab/src/api/pfqn/pfqn_nrl.m

pfqn_nrp(L, N, Z=None, alpha=None)[source]

Norlund-Rice Probit (NRP) approximation for normalizing constant.

Computes the logarithm of the normalizing constant using Laplace approximation with probit (normalized) transformation.

Parameters:

L (numpy.ndarray) – Service demand matrix (M x R)
N (numpy.ndarray) – Population vector (R,)
Z (numpy.ndarray) – Think time vector (R,) - optional
alpha (numpy.ndarray) – Load-dependent rate matrix (M x Ntot) - optional

Returns:

Logarithm of normalizing constant

Return type:

lG

References

Original MATLAB: matlab/src/api/pfqn/pfqn_nrp.m

pfqn_lap(L, N, Z)[source]

Laplace approximation for normalizing constant.

Computes the logarithm of the normalizing constant using the classical Laplace approximation with root finding.

This method uses a saddle-point approximation to estimate the normalizing constant integral representation.

Parameters:

L (numpy.ndarray) – Service demand vector (R,) - single station
N (numpy.ndarray) – Population vector (R,)
Z (numpy.ndarray) – Think time vector (R,)

Returns:

Logarithm of normalizing constant approximation

Return type:

lG

References

Original MATLAB: matlab/src/api/pfqn/pfqn_lap.m

laplaceapprox(h, x0, tol=1e-5)[source]

Laplace approximation for multidimensional integrals.

Approximates I = ∫ h(x) dx using the Laplace method:: I ≈ h(x0) * sqrt((2π)^d / det(-H))

where H is the Hessian of log(h) at x0.

Parameters:

h (Callable) – Function to integrate (must be positive)
x0 (numpy.ndarray) – Point for Laplace approximation (typically the mode)
tol (float) – Tolerance for numerical Hessian computation

Returns:

I: Approximate integral value H: Hessian matrix at x0 logI: Logarithm of integral value (more stable)

Return type:

Tuple (I, H, logI) where

num_hess(f, x0, tol=1e-5)[source]

Compute numerical Hessian matrix using finite differences.

Uses central differences for improved accuracy.

Parameters:

f (Callable) – Function to differentiate (maps x -> scalar)
x0 (numpy.ndarray) – Point at which to compute Hessian
tol (float) – Step size for finite differences

Returns:

Hessian matrix (d x d) where d = len(x0)

Return type:

numpy.ndarray

pfqn_stdf(L, N, Z, S, fcfs_nodes, rates, tset)[source]

Compute sojourn time distribution for multiserver FCFS nodes.

Implements McKenna’s method for computing the response time distribution at multiserver FCFS stations in closed queueing networks.

Parameters:

L (numpy.ndarray) – Load matrix (M x R) - service demands
N (numpy.ndarray) – Population vector (R,) - number of jobs per class
Z (numpy.ndarray) – Think time vector (R,) or matrix (1 x R)
S (numpy.ndarray) – Number of servers at each station (M,)
fcfs_nodes (numpy.ndarray) – Array of FCFS station indices (1-indexed as in MATLAB)
rates (numpy.ndarray) – Service rates matrix (M x R)
tset (numpy.ndarray) – Time points at which to evaluate the distribution

Returns:

Dictionary with (station, class) tuples as keys and response time distribution arrays as values. Each array has shape (len(tset), 2) with:

column 0: CDF values

column 1: time points

Return type:

Dict[Tuple[int, int], numpy.ndarray]

Examples

>>> L = np.array([[0.5, 0.3], [0.2, 0.4]])
>>> N = np.array([2, 3])
>>> Z = np.array([1.0, 1.0])
>>> S = np.array([2, 1])
>>> fcfs_nodes = np.array([0])  # 0-indexed
>>> rates = np.array([[1.0, 1.0], [2.0, 2.0]])
>>> tset = np.linspace(0.1, 5.0, 50)
>>> RD = pfqn_stdf(L, N, Z, S, fcfs_nodes, rates, tset)

pfqn_stdf_heur(L, N, Z, S, fcfs_nodes, rates, tset)[source]

Heuristic sojourn time distribution for multiserver FCFS nodes.

Implements a variant of McKenna’s 1987 method that uses per-class service rate weighting for improved accuracy in multiclass networks.

Unlike pfqn_stdf which uses a simpler Erlang model for the waiting component, this heuristic accounts for class-dependent waiting times based on the expected queue composition.

Parameters:

L (numpy.ndarray) – Load matrix (M x R) - service demands
N (numpy.ndarray) – Population vector (R,) - number of jobs per class
Z (numpy.ndarray) – Think time vector (R,) or matrix (1 x R)
S (numpy.ndarray) – Number of servers at each station (M,)
fcfs_nodes (numpy.ndarray) – Array of FCFS station indices (0-indexed)
rates (numpy.ndarray) – Service rates matrix (M x R)
tset (numpy.ndarray) – Time points at which to evaluate the distribution

Returns:

Dictionary with (station, class) tuples as keys and response time distribution arrays as values. Each array has shape (len(tset), 2) with:

column 0: CDF values

column 1: time points

Return type:

Dict[Tuple[int, int], numpy.ndarray]

References

McKenna, J. “Mean Value Analysis for networks with service-time-dependent queue disciplines.” Performance Evaluation, 1987.

pfqn_cftp(L, N, S=None, nsamples=1, method='cftp')[source]

Exact (perfect) stationary state sampling for closed single-class multiserver product-form networks via monotone Coupling From The Past.

Parameters:

L (numpy.ndarray) – Service demands (stations,), L[i] = theta_i / mu_i
N (int) – Total closed population (scalar)
S (numpy.ndarray) – Number of servers per station (default 1; use np.inf for infinite server). If None, defaults to ones.
nsamples (int) – Number of independent samples to draw (default 1)
method (str) – ‘cftp’ for exact/perfect sampling (default) or ‘approx’ for rapidly-mixing approximate sampler M_A

Returns:

Empirical mean queue length per station (1 x stations) X: Sampled states, one per row (nsamples x stations), each row sums to N T: Per-sample coalescence horizon (‘cftp’) or mixing steps used

(‘approx’), as (nsamples,)

Return type:

Q

Raises:

ValueError – If less than 2 stations or non-positive demands