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)

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 (ndarray) – Service demand matrix (M x R) where M is stations, R is classes

• N (ndarray) – Population vector (1 x R or R,) - number of jobs per class

• Z (ndarray) – Think time vector (1 x R or R,) - think time per class (default 0)

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

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 (ndarray) – Service demands at each station (1D array of length M)

• Z (float) – Think time (default 0)

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

Balanced System (Asymptotic) analysis for product-form networks.

Provides asymptotic approximations based on bottleneck analysis. Fast but less accurate than exact MVA for small populations.

Parameters:

• L (ndarray) – Service demand matrix (M x R)

• N (ndarray) – Population vector

• Z (ndarray) – Think time vector (default 0)

Returns:

Same format as pfqn_mva

Return type:

Tuple[ndarray, ndarray, ndarray, ndarray, ndarray, ndarray, ndarray]

pfqn_aql(L, N, Z=None, max_iter=1000, tol=1e-6)

Approximate Queue Length (AQL) algorithm.

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

Parameters:

• L (ndarray) – Service demand matrix (M x R)

• N (ndarray) – Population vector

• Z (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[ndarray, ndarray, ndarray, ndarray, ndarray, ndarray, ndarray]

pfqn_sqni(L, N, Z=None)

Single Queue Network Interpolation (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 (ndarray) – Service demand vector (1 x R or R,) - demands at the bottleneck queue

• N (ndarray) – Population vector (1 x R or R,) - number of jobs per class

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

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 (ndarray) – Service demand matrix (M x R) - rows are stations, columns are classes

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

QD-Linearizer (QDLIN) Approximate MVA.

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

Parameters:

• L (ndarray) – Service demand matrix (M x R)

• N (ndarray) – Population vector (R,)

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

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 (ndarray) – Service demand matrix (M x R)

• N (ndarray) – Population vector (R,)

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

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 (ndarray) – Service demand matrix (M x R)

• N (ndarray) – Population vector (R,)

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

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 (ndarray) – Service demand matrix (M x R) where M is stations, R is classes

• N (ndarray) – Population vector (1 x R) - number of jobs per class

• Z (ndarray) – Think time vector (1 x R) - default: zeros

• tol (float) – Convergence tolerance (default 1e-6)

• max_iter (int) – Maximum iterations (default 1000)

• QN (ndarray) – Initial queue length matrix (M x R) - default: uniform distribution

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

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 (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 (ndarray) – Service demand matrix (M x R)

• N (ndarray) – Population vector (1 x R)

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

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 (ndarray) – Service demand matrix (M x R) where M is stations, R is classes

• N (ndarray) – Population vector (1 x R or R,) - number of jobs per class

• Z (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_nc(L, N, Z=None, method='ca')

Normalizing constant computation dispatcher.

Selects appropriate algorithm based on method parameter.

Parameters:

• L (ndarray) – Service demand matrix (M x R)

• N (ndarray) – Population vector (R,)

• Z (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’: Reduced decomposition - ‘ls’: Linearizer

Returns:

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

Return type:

Tuple[float, float]

pfqn_panacea(L, N, Z=None)

PANACEA algorithm (hybrid convolution/MVA).

Currently implemented as wrapper around convolution algorithm.

Parameters:

• L (ndarray) – Service demand matrix

• N (ndarray) – Population vector

• Z (ndarray) – Think time vector

Returns:

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

Return type:

Tuple[float, float]

pfqn_propfair(L, N, Z=None)

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 (ndarray) – Service demand matrix (M x R) where M is stations, R is classes

• N (ndarray) – Population vector (1 x R or R,) - number of jobs per class

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

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 (ndarray) – Service demand matrix (M x R)

• N (ndarray) – Population vector (R,)

• Z (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_linearizer(L, N, Z, sched_type, tol=1e-8, maxiter=1000)

Linearizer approximate MVA algorithm.

Parameters:

• L (ndarray) – Demand matrix (M x R) - rows are stations, columns are classes

• N (ndarray) – Population vector (R,)

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

General-form linearizer approximate MVA.

Parameters:

• L (ndarray) – Demand matrix (M x R)

• N (ndarray) – Population vector (R,)

• Z (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[ndarray, ndarray, ndarray, ndarray, ndarray, ndarray, int]

pfqn_egflinearizer(L, N, Z, sched_type, tol=1e-8, maxiter=1000, alpha=None)

Extended general-form linearizer with class-specific parameters.

Parameters:

• L (ndarray) – Demand matrix (M x R)

• N (ndarray) – Population vector (R,)

• Z (ndarray) – Think time vector (R,)

• sched_type (List[str]) – List of scheduling strategies per station

• tol (float) – Convergence tolerance

• maxiter (int) – Maximum iterations

• alpha (ndarray) – Class-specific linearization parameters (R,)

Returns:

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

Return type:

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

class SchedStrategy(*values)

Bases: Enum

Scheduling strategies for queueing stations.

pfqn_mvald(L, N, Z, mu, stabilize=True)

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 (ndarray) – Service demand matrix (M x R) - rows are stations, columns are classes.

• N (ndarray) – Population vector (R,).

• Z (ndarray) – Think time vector (R,).

• mu (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 x R). 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)

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 (ndarray) – Arrival rate vector (R,). Use 0 for closed classes.

• L (ndarray) – Service demand matrix (M x R).

• N (ndarray) – Population vector (R,). Use np.inf for open classes.

• Z (ndarray) – Think time vector (R,).

• mi (ndarray | None) – Queue replication factors (M,) (default: ones).

• S (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 (M x R). CN: Cycle times (1 x R). 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)

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 (ndarray) – Arrival rate vector (R,). Use 0 for closed classes.

• L (ndarray) – Service demand matrix (M x R).

• N (ndarray) – Population vector (R,). Use np.inf for open classes.

• Z (ndarray) – Think time vector (R,).

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

Lower asymptotic bound on throughput (Zahorjan-Balanced).

Provides a simple lower bound on system throughput for single-class closed queueing networks.

Parameters:

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

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

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

Parameters:

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

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

Parameters:

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

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

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 (ndarray) – Service demand vector (M,).

• N (int | float) – Population (scalar).

• Z (float) – Think time.

Returns:

Upper bound on throughput.

Return type:

float

pfqn_le(L, N, Z=None)

Laguerre 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 (ndarray) – Service demand matrix (M x R).

• N (ndarray) – Population vector (R,).

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

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 (ndarray) – Service demand matrix (M x R).

• N (ndarray) – Population vector (R,).

• Z (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')

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 (ndarray) – Service demand matrix (M x R).

• N (ndarray) – Population vector (R,).

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

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 (ndarray) – Service demand matrix (M x R).

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

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

Public wrapper for the internal _pfqn_le_fpi function.

Parameters:

• L (ndarray) – Service demand matrix (M x R).

• N (ndarray) – Population vector (R,).

Returns:

Mode location vector u (M,).

Return type:

ndarray

pfqn_le_fpiZ(L, N, Z)

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

Public wrapper for the internal _pfqn_le_fpiZ function.

Parameters:

• L (ndarray) – Service demand matrix (M x R).

• N (ndarray) – Population vector (R,).

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

Compute Hessian matrix (no think time case).

Public wrapper for the internal _pfqn_le_hessian function.

Parameters:

• L (ndarray) – Service demand matrix (M x R).

• N (ndarray) – Population vector (R,).

• u (ndarray) – Mode location vector (M,).

Returns:

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

Return type:

ndarray

pfqn_le_hessianZ(L, N, Z, u, v)

Compute Hessian matrix (with think time case).

Public wrapper for the internal _pfqn_le_hessianZ function.

Parameters:

• L (ndarray) – Service demand matrix (M x R).

• N (ndarray) – Population vector (R,).

• Z (ndarray) – Think time vector (R,).

• u (ndarray) – Mode location vector (M,).

• v (float) – Scale factor.

Returns:

Hessian matrix (M x M).

Return type:

ndarray

pfqn_ncld(L, N, Z, mu, options=None)

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 (ndarray) – Service demands at all stations (M x R)

• N (ndarray) – Number of jobs for each class (1 x R)

• Z (ndarray) – Think times for each class (1 x R)

• mu (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_gld(L, N, mu, options=None)

Compute normalizing constant of a load-dependent closed queueing network.

Uses the generalized convolution algorithm for computing normalizing constants in load-dependent closed queueing networks.

Parameters:

• L (ndarray) – Service demands at all stations (M x R)

• N (ndarray) – Number of jobs for each class (1 x R)

• mu (ndarray) – Load-dependent scalings (M x Ntot)

• options (Dict[str, Any] | None) – Solver options

Returns:

PfqnNcResult with G (normalizing constant) and lG (log)

Return type:

PfqnNcResult

pfqn_gldsingle(L, N, mu, options=None)

Compute normalizing constant for single-class load-dependent model.

Auxiliary function used by pfqn_gld to compute the normalizing constant in a single-class load-dependent model using dynamic programming.

Parameters:

• L (ndarray) – Service demands at all stations (M x 1)

• N (ndarray) – Number of jobs (scalar or 1x1 array)

• mu (ndarray) – Load-dependent scaling factors (M x Ntot)

• options (Dict[str, Any] | None) – Solver options (unused, for API compatibility)

Returns:

PfqnNcResult with G (normalizing constant) and lG (log)

Raises:

RuntimeError – If multiclass model is detected

Return type:

PfqnNcResult

pfqn_mushift(mu, k)

Shift a load-dependent scaling vector by one position.

Used in recursive normalizing constant computations.

Parameters:

• mu (ndarray) – Load-dependent scalings matrix (M x N)

• k (int) – Row index to shift

Returns:

Shifted mu matrix (M x N-1)

Return type:

ndarray

pfqn_comomrm_ld(L, N, Z, mu, options=None)

Run the COMOM normalizing constant method on a load-dependent repairman model.

Implements the Convolution Method of Moments (COMOM) for computing normalizing constants in load-dependent repairman queueing models.

Parameters:

• L (ndarray) – Service demands at all stations (M x R)

• N (ndarray) – Number of jobs for each class (1 x R)

• Z (ndarray) – Think times for each class (1 x R)

• mu (ndarray) – Load-dependent scalings (M x Ntot)

• options (Dict[str, Any] | None) – Solver options

Returns:

PfqnComomrmLdResult with G, lG, and marginal probabilities

Return type:

PfqnComomrmLdResult

pfqn_fnc(alpha, c=None)

Compute scaling factor of a load-dependent functional server.

Used to calculate the mean queue length in load-dependent systems by computing functional scaling factors from load-dependent service rate parameters.

Parameters:

• alpha (ndarray) – Load-dependent scalings (M x N)

• c (ndarray | None) – Scaling constants (1 x M), optional. If None, auto-selected.

Returns:

PfqnFncResult with mu (functional server scalings) and c (scaling constants)

Return type:

PfqnFncResult

class PfqnNcResult(G, lG, method='default')

Bases: object

Result of normalizing constant computation.

method: str = 'default'

G: float

lG: float

class PfqnComomrmLdResult(G, lG, prob)

Bases: object

Result of COMOM load-dependent computation.

G: float

lG: float

prob: ndarray

class PfqnFncResult(mu, c)

Bases: object

Result of functional server scaling computation.

mu: ndarray

c: ndarray

pfqn_unique(L, mu=None, gamma=None, tol=1e-14)

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 (ndarray) – Service demand matrix (M x R)

• mu (ndarray | None) – Load-dependent rate matrix (M x Ntot), optional - pass None if not used

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

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 (ndarray) – Queue lengths from reduced model (M’ x R)

• UN (ndarray) – Utilizations from reduced model (M’ x R)

• CN (ndarray) – Cycle times from reduced model (M’ x R)

• mapping (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[ndarray, ndarray, ndarray]

pfqn_combine_mi(mi, mapping, M_unique)

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 (ndarray) – User-provided multiplicity vector (1 x M_original or M_original,)

• mapping (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:

ndarray

class PfqnUniqueResult(L_unique, mu_unique, gamma_unique, mi, mapping)

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: ndarray

Alias for field number 0

gamma_unique: ndarray | None

Alias for field number 2

mapping: ndarray

Alias for field number 4

mi: ndarray

Alias for field number 3

mu_unique: ndarray | None

Alias for field number 1

pfqn_lldfun(n, lldscaling=None, nservers=None)

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 (ndarray) – Queue population vector (M,)

• lldscaling (ndarray | None) – Load-dependent scaling matrix (M x Nmax), optional

• nservers (ndarray | None) – Number of servers per station (M,), optional

Returns:

Scaling factor vector (M,)

Return type:

ndarray

References

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

pfqn_mu_ms(N, m, c)

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:

ndarray

References

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

pfqn_nc_sanitize(lam, L, N, Z, atol=1e-8)

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 (ndarray) – Arrival rate vector

• L (ndarray) – Service demand matrix (M x R)

• N (ndarray) – Population vector

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

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

Computes the scaling factor at each station based on the class-dependent population state.

Parameters:

• nvec (ndarray) – Population state matrix (M x R) or vector (M,)

• cdscaling (List | None) – List of class-dependent scaling functions, one per station. Each function takes a class population vector and returns a scalar.

Returns:

Scaling factor vector (M,)

Return type:

ndarray

References

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

pfqn_ljdfun(n, ljdscaling=None)

AMVA-QD limited joint-dependence function.

Computes scaling factors for stations with limited joint-dependence, where the service rate depends on the total population across stations.

Parameters:

• n (ndarray) – Queue population vector (M,)

• ljdscaling (ndarray | None) – Limited joint-dependence scaling matrix (M x Nmax), optional

Returns:

Scaling factor vector (M,)

Return type:

ndarray

References

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

factln(n)

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)

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

Parameters:

arr (ndarray) – Array of non-negative numbers

Returns:

Array of log(n!) values

Return type:

ndarray

softmin(a, b, alpha=20.0)

Compute a smooth approximation to min(a, b).

Uses the log-sum-exp trick for numerical stability.

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)

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

Parameters:

• n (ndarray) – Population vector

• s (int) – Index to decrement (0-based)

Returns:

Copy of n with n[s] -= 1

Return type:

ndarray

multichoose(r, n)

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:

ndarray

matchrow(matrix, row)

Find the index of a row in a matrix.

Parameters:

• matrix (ndarray) – 2D array to search in

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

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 (ndarray) – Service demand matrix (M x R)

• N (ndarray) – Population vector (R,)

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

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 (ndarray) – Service demand matrix (1 x R) - single station

• N (ndarray) – Population vector (R,)

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

Original CoMoM implementation for repairman model.

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

Parameters:

• L (ndarray) – Service demand matrix (1 x R)

• N (ndarray) – Population vector (R,)

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

CoMoM for multi-server repairman model.

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

Parameters:

• L (ndarray) – Service demand matrix (1 x R)

• N (ndarray) – Population vector (R,)

• Z (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_procomom2(L, N, Z=None, atol=1e-8)

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 (ndarray) – Service demand matrix (M x R)

• N (ndarray) – Population vector (R,)

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

Bases: NamedTuple

Result of CoMoM normalizing constant computation.

Create new instance of ComomResult(lG, lGbasis)

lG: float

Alias for field number 0

lGbasis: ndarray | None

Alias for field number 1

pfqn_mmint2(L, N, Z, m=1)

McKenna-Mitra integral form using scipy.integrate.

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

Parameters:

• L (ndarray) – Service demand vector (R,)

• N (ndarray) – Population vector (R,)

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

McKenna-Mitra integral with Gauss-Legendre quadrature.

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

Parameters:

• L (ndarray) – Service demand vector (R,)

• N (ndarray) – Population vector (R,)

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

McKenna-Mitra integral with Gauss-Laguerre quadrature.

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

Parameters:

• L (ndarray) – Service demand vector (R,)

• N (ndarray) – Population vector (R,)

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

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 (ndarray) – Service demand vector (R,)

• N (ndarray) – Population vector (R,)

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

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

Parameters:

x (ndarray) – Array of values

Returns:

log(sum(exp(x)))

Return type:

float

pfqn_schmidt(D, N, S, sched, v=None)

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 (ndarray) – Service demand matrix (M x R)

• N (ndarray) – Population vector (R,)

• S (ndarray) – Number of servers per station (M,) or (M x R)

• sched (ndarray) – Scheduling discipline per station (M,) - SchedStrategy values

• v (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) - returns zeros (not computed) 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)

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 (ndarray) – Service demand matrix (M x R)

• N (ndarray) – Population vector (R,)

• S (ndarray) – Number of servers per station (M,) or (M x R)

• sched (ndarray) – Scheduling discipline per station (M,) - SchedStrategy values

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

Bases: object

Result from Schmidt’s exact MVA.

XN: ndarray

QN: ndarray

UN: ndarray

CN: ndarray

pprod(N, Nmax=None)

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 (ndarray) – Current population vector or max population

• Nmax (ndarray | None) – Maximum population per class (if incrementing)

Returns:

Next population vector, or array of -1 if done

Return type:

ndarray

hashpop(nvec, Nc, C, prods)

Hash population vector to linear index.

Parameters:

• nvec (ndarray) – Population vector

• Nc (ndarray) – Maximum population per class

• C (int) – Number of classes

• prods (ndarray) – Precomputed products for hashing

Returns:

Linear index (1-based for MATLAB compatibility)

Return type:

int

pfqn_recal(L, N, Z=None, m0=None)

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 (ndarray) – Service demand matrix (M x R)

• N (ndarray) – Population vector (R,)

• Z (ndarray | None) – Think time vector (R,), optional (default: zeros)

• m0 (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_mvaldmx(lam, D, N, Z, mu=None, S=None)

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 (ndarray) – Arrival rate vector (R,) - non-zero for open classes

• D (ndarray) – Service demand matrix (M x R)

• N (ndarray) – Population vector (R,) - inf for open classes

• Z (ndarray) – Think time vector (R,)

• mu (ndarray | None) – Load-dependent rate matrix (M x Nt), optional

• S (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_mvaldmx_ec(lam, D, mu)

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 (ndarray) – Arrival rate vector (R,)

• D (ndarray) – Service demand matrix (M x R)

• mu (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_mvaldmx_ec.m

pfqn_mvaldms(lam, D, N, Z, S)

Load-dependent MVA for multiserver mixed networks.

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

Parameters:

• lam (ndarray) – Arrival rate vector (R,)

• D (ndarray) – Service demand matrix (M x R)

• N (ndarray) – Population vector (R,)

• Z (ndarray) – Think time vector (R,)

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

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

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

Parameters:

• L (ndarray) – Service demand matrix (M x R)

• N (ndarray) – Population vector (R,)

• Z (ndarray) – Think time vector (R,)

• nservers (ndarray) – Number of servers per station (M,)

• type_sched (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(lam, L, N, Z, nservers, type_sched=None, tol=1e-8, maxiter=1000, method='egflin')

Linearizer for mixed open/closed queueing networks.

Handles networks with both open and closed classes using a decomposition approach.

Parameters:

• lam (ndarray) – Arrival rate vector (R,) - inf for closed classes

• L (ndarray) – Service demand matrix (M x R)

• N (ndarray) – Population vector (R,) - inf for open classes

• Z (ndarray) – Think time vector (R,)

• nservers (ndarray) – Number of servers per station (M,)

• type_sched (ndarray | None) – Scheduling strategy per station (M,), optional

• tol (float) – Convergence tolerance (default: 1e-8)

• maxiter (int) – Maximum iterations (default: 1000)

• method (str) – Linearizer variant (‘lin’, ‘gflin’, ‘egflin’, default: ‘egflin’)

Returns:

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

Return type:

Tuple of (QN, UN, WN, TN, CN, XN, totiter)

References

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

pfqn_conwayms(L, N, Z, nservers, type_sched=None, tol=1e-8, maxiter=1000)

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 (ndarray) – Service demand matrix (M x R)

• N (ndarray) – Population vector (R,)

• Z (ndarray) – Think time vector (R,)

• nservers (ndarray) – Number of servers per station (M,)

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

Convert per-class population vector to linearized index.

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

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

Parameters:

• nvec (ndarray) – Per-class populations [n1, n2, …, nK]

• cutoffs (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

ljd_delinearize(idx, cutoffs)

Convert linearized index back to per-class population vector.

Inverse of ljd_linearize: recovers the multi-dimensional population vector from a linear index.

Parameters:

• idx (int) – 1-based linearized index

• cutoffs (ndarray) – Per-class cutoffs [N1, N2, …, NK]

Returns:

Per-class populations [n1, n2, …, nK]

Return type:

ndarray

References

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

ljcd_interpolate(nvec, cutoffs, table, K)

Multi-linear interpolation for LJCD scaling tables.

Performs multi-linear interpolation of a throughput value from an LJCD scaling table for non-integer population vectors.

For K classes, interpolates between 2^K corner points of the hypercube containing the population vector using multi-linear interpolation.

Parameters:

• nvec (ndarray) – Continuous population vector [n1, n2, …, nK] (clamped to cutoffs)

• cutoffs (ndarray) – Per-class cutoffs [N1, N2, …, NK]

• table (ndarray) – Linearized throughput table (1-D array indexed by ljd_linearize)

• K (int) – Number of classes

Returns:

Interpolated throughput value

Return type:

float

References

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

infradius_h(x, L, N, alpha)

Helper function for infinite radius computation with logistic transformation.

Used in normalizing constant computation via integration methods.

Parameters:

• x (ndarray) – Logistic transformation parameters

• L (ndarray) – Service demand matrix (M x R)

• N (ndarray) – Population vector (R,)

• alpha (ndarray) – Load-dependent rate matrix

Returns:

Evaluated function value for integration

Return type:

ndarray

References

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

infradius_hnorm(x, L, N, alpha, logNormConstScale)

Normalized helper function for infinite radius computation.

Like infradius_h but with normalization for numerical stability.

Parameters:

• x (ndarray) – Logistic transformation parameters

• L (ndarray) – Service demand matrix (M x R)

• N (ndarray) – Population vector (R,)

• alpha (ndarray) – Load-dependent rate matrix

• logNormConstScale (float) – Logarithm of normalizing scale factor

Returns:

Normalized evaluated function value

Return type:

ndarray

References

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

pfqn_kt(L, N, Z=None)

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 (ndarray) – Service demand matrix (M x R)

• N (ndarray) – Population vector (R,)

• Z (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')

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

Parameters:

• D (ndarray) – Service time matrix (M x K)

• N (ndarray) – Population vector (1 x K)

• V (ndarray) – Visit ratio matrix (M x K)

• nservers (ndarray) – Number of servers at each station (M x 1)

• sched (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')

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 (ndarray) – Population vector (K,)

• nservers (ndarray) – Number of servers at each station (M,)

• sched_type (ndarray) – Scheduling strategies for each station (M,)

• v (ndarray) – Visit ratio matrix (M x K)

• s (ndarray) – Service time matrix (M x K)

• maxiter (int) – Maximum iterations

• D_frac (ndarray) – Fractional changes matrix (M x K x K)

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

Bases: object

Result from Akyildiz-Bolch AMVA algorithm.

QN: ndarray

UN: ndarray

RN: ndarray

CN: ndarray

XN: ndarray

totiter: int

pfqn_rd(L, N, Z, mu=None, options=None)

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 (ndarray) – Service demand matrix (M x R)

• N (ndarray) – Population vector (1 x R)

• Z (ndarray) – Think time vector (1 x R)

• mu (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')

Bases: object

Options for RD algorithm.

method: str = 'default'

tol: float = 1e-06

class RdResult(lGN, Cgamma)

Bases: object

Result from Reduced Decomposition algorithm.

lGN: float

Cgamma: float

pfqn_nrl(L, N, Z=None, alpha=None)

Normal Radius-Logistic (NRL) approximation for normalizing constant.

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

Parameters:

• L (ndarray) – Service demand matrix (M x R)

• N (ndarray) – Population vector (R,)

• Z (ndarray) – Think time vector (R,) - optional

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

Normal Radius-Probit (NRP) approximation for normalizing constant.

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

Parameters:

• L (ndarray) – Service demand matrix (M x R)

• N (ndarray) – Population vector (R,)

• Z (ndarray) – Think time vector (R,) - optional

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

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 (ndarray) – Service demand vector (R,) - single station

• N (ndarray) – Population vector (R,)

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

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

Compute numerical Hessian matrix using finite differences.

Uses central differences for improved accuracy.

Parameters:

• f (Callable) – Function to differentiate (maps x -> scalar)

• x0 (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:

ndarray

pfqn_stdf(L, N, Z, S, fcfs_nodes, rates, tset)

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 (ndarray) – Load matrix (M x R) - service demands

• N (ndarray) – Population vector (R,) - number of jobs per class

• Z (ndarray) – Think time vector (R,) or matrix (1 x R)

• S (ndarray) – Number of servers at each station (M,)

• fcfs_nodes (ndarray) – Array of FCFS station indices (1-indexed as in MATLAB)

• rates (ndarray) – Service rates matrix (M x R)

• tset (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], 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)

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 (ndarray) – Load matrix (M x R) - service demands

• N (ndarray) – Population vector (R,) - number of jobs per class

• Z (ndarray) – Think time vector (R,) or matrix (1 x R)

• S (ndarray) – Number of servers at each station (M,)

• fcfs_nodes (ndarray) – Array of FCFS station indices (0-indexed)

• rates (ndarray) – Service rates matrix (M x R)

• tset (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], ndarray]

References

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