API Documentation

This section provides detailed API documentation for all LINE Solver Python modules.

API Submodules

Cache API (line_solver.api.cache)

Cache Analysis Algorithms.

Native Python implementations for analyzing cache systems, including exact recursive methods, singular perturbation methods, importance sampling, TTL-based caches, and various approximation techniques.

Key algorithms:

cache_erec: Exact recursive normalizing constant cache_prob_erec: Exact recursive state probabilities cache_spm: Singular perturbation method cache_miss: Miss rate computation cache_is: Importance sampling cache_ttl_*: TTL-based cache analysis cache_rrm_*: Random replacement model

cache_erec(gamma, m)[source]

Compute the cache normalizing constant using exact recursive method.

This method serves as a wrapper that calls the auxiliary function to perform the actual computation.

Parameters:

  • gamma (numpy.ndarray) – Cache access factors matrix (n x h), where n is number of items and h is number of cache levels.

  • m (numpy.ndarray) – Cache capacity vector (1 x h or h,).

Returns:

Normalizing constant E.

Return type:

float

cache_erec_aux(gamma, m, k)[source]

Auxiliary method for computing cache normalizing constant using exact recursive method.

This method performs the core computation recursively, adjusting the size of the input matrix.

Parameters:

  • gamma (numpy.ndarray) – Cache access factors matrix (n x h).

  • m (numpy.ndarray) – Cache capacity vector (h,).

  • k (int) – Current number of rows in the recursive step.

Returns:

Normalizing constant for the given configuration.

Return type:

float

cache_prob_erec(gamma, m)[source]

Compute cache state probabilities using exact recursive method.

This method calculates the probabilities of the cache being in different states based on the cache access factors and capacity.

Parameters:

  • gamma (numpy.ndarray) – Cache access factors matrix (n x h), where n is number of items and h is number of cache levels.

  • m (numpy.ndarray) – Cache capacity vector (1 x h or h,).

Returns:

Matrix (n x h+1) containing cache state probabilities. Column 0 is miss probability, columns 1..h are hit probabilities at each cache level.

Return type:

numpy.ndarray

cache_mva(gamma, m)[source]

Mean Value Analysis for cache systems.

Computes cache performance metrics using MVA approach.

Parameters:
Returns:

  • pi: Steady-state probabilities

  • pi0: Miss probabilities per item

  • pij: Hit probabilities per item per level (n x h)

  • x: Throughput vector

  • u: Utilization vector

  • E: Normalizing constant

Return type:

Tuple of (pi, pi0, pij, x, u, E) containing

cache_xi_iter(gamma, m, tol=1e-12, max_iter=1000)[source]

Compute cache xi terms using iterative method from Gast-van Houdt.

This method calculates the xi values which are important for understanding the distribution of items in the cache. The algorithm assumes that the access factors are monotone with the list index.

Parameters:

  • gamma (numpy.ndarray) – Cache access factors matrix (n x h).

  • m (numpy.ndarray) – Cache capacity vector (h,).

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

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

Returns:

Vector of xi values (h,).

Return type:

numpy.ndarray

cache_spm(gamma, m)[source]

Approximate the normalizing constant using SPM method.

Computes the normalizing constant of the cache steady state distribution using the singular perturbation method (also known as ray integration).

Parameters:
Returns:

  • Z: Approximated normalizing constant

  • lZ: Log of normalizing constant

  • xi: Xi terms vector

Return type:

Tuple of (Z, lZ, xi) where

cache_prob_spm(gamma, m)[source]

Compute cache miss probabilities using singular perturbation method.

Uses the SPM (ray integration) method to compute miss probabilities for cache systems.

Parameters:
Returns:

Vector of miss probabilities for each item.

Return type:

numpy.ndarray

cache_prob_fpi(gamma, m, tol=1e-8, max_iter=1000)[source]

Compute cache hit probabilities using fixed point iteration (FPI method).

Matches MATLAB cache_prob_fpi: Uses fixed point iteration to compute cache hit probability distribution.

Parameters:

  • gamma (numpy.ndarray) – Cache access factors matrix (n x h).

  • m (numpy.ndarray) – Cache capacity vector (h,).

  • tol (float) – Convergence tolerance (unused, kept for API compatibility).

  • max_iter (int) – Maximum iterations (unused, kept for API compatibility).

Returns:

  • Column 0: miss probabilities

  • Columns 1:h: hit probabilities at each level

Return type:

Matrix (n x h+1) where

cache_miss(gamma, m, lambd=None)[source]

Compute miss rates for a cache system using exact recursive method.

Parameters:

  • gamma (numpy.ndarray) – Item popularity probabilities (n x h matrix)

  • m (numpy.ndarray) – Cache capacity vector (h,)

  • lambd (numpy.ndarray | None) – Optional arrival rates per user per item (u x n x h+1)

Returns:

  • M: Global miss rate

  • MU: Per-user miss rate (u,) or None

  • MI: Per-item miss rate (n,) or None

  • pi0: Per-item miss probability (n,) or None

Return type:

Tuple of (M, MU, MI, pi0) where

References

Original MATLAB: matlab/src/api/cache/cache_miss.m

cache_xi_fp(gamma, m, xi_init=None, tol=1e-14, max_iter=10000)[source]

Fixed-point iteration for computing cache performance metrics.

Computes cache performance metrics including Lagrange multipliers, miss probabilities, and hit probabilities.

Parameters:

  • gamma (numpy.ndarray) – Item popularity probabilities (n x h)

  • m (numpy.ndarray) – Cache capacity vector (h,)

  • xi_init (numpy.ndarray | None) – Optional initial guess for Lagrange multipliers

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

  • max_iter (int) – Maximum iterations (default: 10000)

Returns:

  • xi: Converged Lagrange multipliers (h,)

  • pi0: Miss probability per item (n,)

  • pij: Hit probability per item per list (n x h)

  • it: Number of iterations

Return type:

Tuple of (xi, pi0, pij, it) where

References

Original MATLAB: matlab/src/api/cache/cache_xi_fp.m

cache_miss_fpi(gamma, m, lambd=None)[source]

Compute cache miss rates using fixed-point iteration method.

Parameters:

  • gamma (numpy.ndarray) – Item popularity probabilities (n x h)

  • m (numpy.ndarray) – Cache capacity vector (h,)

  • lambd (numpy.ndarray | None) – Optional arrival rates per user per item (u x n x h+1)

Returns:

  • M: Global miss rate

  • MU: Per-user miss rate or None

  • MI: Per-item miss rate or None

  • pi0: Per-item miss probability or None

Return type:

Tuple of (M, MU, MI, pi0) where

References

Original MATLAB: matlab/src/api/cache/cache_miss_fpi.m

cache_miss_spm(gamma, m, lambd=None)[source]

Compute cache miss rates using singular perturbation method.

Parameters:

  • gamma (numpy.ndarray) – Item popularity probabilities (n x h)

  • m (numpy.ndarray) – Cache capacity vector (h,)

  • lambd (numpy.ndarray | None) – Optional arrival rates per user per item (u x n x h+1)

Returns:

  • M: Global miss rate

  • MU: Per-user miss rate or None

  • MI: Per-item miss rate or None

  • pi0: Per-item miss probability or None

  • lE: Log of normalizing constant

Return type:

Tuple of (M, MU, MI, pi0, lE) where

References

Original MATLAB: matlab/src/api/cache/cache_miss_spm.m

cache_mva_miss(p, m, R)[source]

Compute cache miss rates using Mean Value Analysis.

Parameters:
Returns:

  • M: Global miss rate

  • Mk: Per-item miss rate (n,)

Return type:

Tuple of (M, Mk) where

References

Original MATLAB: matlab/src/api/cache/cache_mva_miss.m

cache_is(gamma, m, samples=100000)[source]

Importance sampling estimation of cache normalizing constant.

Estimates the normalizing constant for cache models using Monte Carlo importance sampling.

Parameters:

  • gamma (numpy.ndarray) – Item popularity probabilities (n x h matrix)

  • m (numpy.ndarray) – Cache capacity vector (h,)

  • samples (int) – Number of Monte Carlo samples (default: 100000)

Returns:

  • E: Normalizing constant estimate

  • lE: Log of normalizing constant

Return type:

Tuple of (E, lE) where

References

Original MATLAB: matlab/src/api/cache/cache_is.m

cache_prob_is(gamma, m, samples=100000)[source]

Importance sampling estimation of cache hit probabilities.

Estimates cache hit probability distribution using Monte Carlo importance sampling.

Parameters:

  • gamma (numpy.ndarray) – Item popularity probabilities (n x h matrix)

  • m (numpy.ndarray) – Cache capacity vector (h,)

  • samples (int) – Number of Monte Carlo samples (default: 100000)

Returns:

prob[i, 0] = miss probability for item i prob[i, 1+j] = hit probability for item i at level j

Return type:

Cache hit probability matrix (n x h+1)

References

Original MATLAB: matlab/src/api/cache/cache_prob_is.m

cache_miss_is(gamma, m, lambd=None, samples=100000)[source]

Importance sampling estimation of cache miss rates.

Computes global, per-user, and per-item miss rates using Monte Carlo importance sampling.

Parameters:

  • gamma (numpy.ndarray) – Item popularity probabilities (n x h matrix)

  • m (numpy.ndarray) – Cache capacity vector (h,)

  • lambd (numpy.ndarray | None) – Optional arrival rates per user per item (u x n x h+1)

  • samples (int) – Number of Monte Carlo samples (default: 100000)

Returns:

  • M: Global miss rate

  • MU: Per-user miss rate or None

  • MI: Per-item miss rate or None

  • pi0: Per-item miss probability or None

  • lE: Log of normalizing constant

Return type:

Tuple of (M, MU, MI, pi0, lE) where

References

Original MATLAB: matlab/src/api/cache/cache_miss_is.m

logmeanexp(x)[source]

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

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

cache_t_lrum(gamma, m)[source]

Compute characteristic times for LRU(m) cache levels.

Uses fixed-point iteration (fsolve) to compute the characteristic time for each level of an LRU(m) cache.

Parameters:

  • gamma (numpy.ndarray) – Item popularity probabilities or arrival rates (n x h)

  • m (numpy.ndarray) – Cache capacity vector (h,)

Returns:

Characteristic time for each cache level (h,)

Return type:

numpy.ndarray

References

Original MATLAB: matlab/src/api/cache/cache_t_lrum.m

cache_t_hlru(gamma, m)[source]

Compute characteristic times for hierarchical LRU cache levels.

Uses fixed-point iteration (fsolve) to compute the characteristic time for each level of a hierarchical LRU cache.

Parameters:

  • gamma (numpy.ndarray) – Item popularity probabilities or arrival rates (n x h)

  • m (numpy.ndarray) – Cache capacity vector (h,)

Returns:

Characteristic time for each cache level (h,)

Return type:

numpy.ndarray

References

Original MATLAB: matlab/src/api/cache/cache_t_hlru.m

cache_ttl_lrum(gamma, m)[source]

Compute steady-state probabilities for TTL-based LRU(m) cache.

Parameters:

  • gamma (numpy.ndarray) – Arrival rates per item per list. Can be (n x h) or (1 x n x h+1). If 3D, uses gamma[0, :, 1:] as the (n x h) input.

  • m (numpy.ndarray) – Cache capacity vector (h,)

Returns:

prob[:, 0] = miss probability prob[:, 1:] = hit probability at each level

Return type:

Steady-state probability distribution (n x h+1)

References

Original MATLAB: matlab/src/api/cache/cache_ttl_lrum.m

cache_ttl_hlru(gamma, m)[source]

Compute steady-state probabilities for TTL-based hierarchical LRU cache.

Parameters:

  • gamma (numpy.ndarray) – Arrival rates per item per list. Can be (n x h) or (1 x n x h+1). If 3D, uses gamma[0, :, 1:] as the (n x h) input.

  • m (numpy.ndarray) – Cache capacity vector (h,)

Returns:

prob[:, 0] = miss probability prob[:, 1] = hit probability at level h

Return type:

Steady-state probability distribution (n x 2)

References

Original MATLAB: matlab/src/api/cache/cache_ttl_hlru.m

cache_ttl_lrua(lambd, R, m, seed=23000)[source]

Compute steady-state probabilities for TTL-LRU cache with arrival routing.

Uses fixed-point iteration with DTMC solving for cache systems with multiple users, items, and levels with routing.

Parameters:

  • lambd (numpy.ndarray) – Arrival rates per user per item per list (u x n x h+1)

  • R (list) – Routing probability structure. Can be either: - 1D list: R[i] is the (h+1 x h+1) routing matrix for item i - 2D list: R[v][i] is the routing matrix for user v, item i

  • m (numpy.ndarray) – Cache capacity vector (h,)

  • seed (int) – Random seed for initialization (default: 23000)

Returns:

Steady-state probability distribution (n x h+1)

Return type:

numpy.ndarray

References

Original MATLAB: matlab/src/api/cache/cache_ttl_lrua.m

cache_rrm_meanfield_ode(t, x, lambd, m, n, h)[source]

ODE function for RRM mean-field cache dynamics.

Defines the differential equations for the Random Replacement Model mean-field cache dynamics.

Parameters:

  • t (float) – Time variable (unused, for ODE solver compatibility)

  • x (numpy.ndarray) – State vector of length n*(h+1), representing probabilities

  • lambd (numpy.ndarray) – Arrival rates per item (n,)

  • m (numpy.ndarray) – Cache capacity vector (h,)

  • n (int) – Number of items

  • h (int) – Number of cache levels

Returns:

Time derivative of state vector

Return type:

numpy.ndarray

References

Original MATLAB: matlab/src/api/cache/cache_rrm_meanfield_ode.m

cache_rrm_meanfield(lambd, m, t_end=10000.0, seed=23000)[source]

Solve RRM mean-field steady state using ODE integration.

Computes the steady-state probability distribution for a cache with random replacement policy using mean-field ODE dynamics.

Parameters:

  • lambd (numpy.ndarray) – Arrival rates per item (n,)

  • m (numpy.ndarray) – Cache capacity vector (h,)

  • t_end (float) – End time for ODE integration (default: 10000.0)

  • seed (int) – Random seed for initial conditions (default: 23000)

Returns:

  • prob: Steady-state probability matrix (n x h+1)

  • missrate: Global miss rate (lambda * miss_prob)

  • missratio: Miss ratio (missrate / sum(lambda))

Return type:

Tuple of (prob, missrate, missratio) where

References

Original MATLAB: matlab/src/api/cache/cache_rrm_meanfield.m

cache_gamma_lp(lambd, R)[source]

Compute gamma parameters for cache models using linear programming approach.

Computes item popularity probabilities at each cache level based on arrival rates and routing probabilities.

Parameters:

  • lambd (numpy.ndarray) – Arrival rates per user per item per list (u x n x h+1)

  • R (list) – Routing probability structure (list of lists, R[v][i] is matrix for user v, item i)

Returns:

  • gamma: Item popularity probabilities at each level (n x h)

  • u: Number of users

  • n: Number of items

  • h: Number of cache levels

Return type:

Tuple of (gamma, u, n, h) where

References

Original MATLAB: matlab/src/api/cache/cache_gamma_lp.m

Product-Form Queueing Networks (line_solver.api.pfqn)

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

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 (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_nc(L, N, Z=None, method='ca')[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’: Reduced decomposition - ‘ls’: Linearizer

Returns:

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

Return type:

Tuple[float, float]

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

PANACEA algorithm (hybrid convolution/MVA).

Currently implemented as wrapper around convolution algorithm.

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

class SchedStrategy(*values)[source]

Bases: Enum

Scheduling strategies for queueing stations.

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 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)[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 (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)[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 asymptotic bound on throughput (Zahorjan-Balanced).

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

Parameters:
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:
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:
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:
Returns:

Upper bound on throughput.

Return type:

float

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

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:
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:
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:
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:
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:
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:
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_gld(L, N, mu, options=None)[source]

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

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

Return type:

PfqnNcResult

pfqn_gldsingle(L, N, mu, options=None)[source]

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

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

  • mu (numpy.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)[source]

Shift a load-dependent scaling vector by one position.

Used in recursive normalizing constant computations.

Parameters:

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

  • k (int) – Row index to shift

Returns:

Shifted mu matrix (M x N-1)

Return type:

numpy.ndarray

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

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

PfqnComomrmLdResult with G, lG, and marginal probabilities

Return type:

PfqnComomrmLdResult

pfqn_fnc(alpha, c=None)[source]

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 (numpy.ndarray) – Load-dependent scalings (M x N)

  • c (numpy.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')[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:
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)[source]

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 (numpy.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:

numpy.ndarray

References

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

pfqn_ljdfun(n, ljdscaling=None)[source]

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

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

Returns:

Scaling factor vector (M,)

Return type:

numpy.ndarray

References

Original MATLAB: matlab/src/api/pfqn/pfqn_ljdfun.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:
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:
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:
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 Convolution 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:
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:
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:
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:
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) - 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)[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:
Returns:

Linear index (1-based for MATLAB compatibility)

Return type:

int

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_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_mvaldmx_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:
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)[source]

Load-dependent MVA for multiserver mixed networks.

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

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

Linearizer for mixed open/closed queueing networks.

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

Parameters:

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

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

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

  • type_sched (numpy.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)[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 table lookups in LJD 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

ljd_delinearize(idx, cutoffs)[source]

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 (numpy.ndarray) – Per-class cutoffs [N1, N2, …, NK]

Returns:

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

Return type:

numpy.ndarray

References

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

ljcd_interpolate(nvec, cutoffs, table, K)[source]

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 (numpy.ndarray) – Continuous population vector [n1, n2, …, nK] (clamped to cutoffs)

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

  • table (numpy.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)[source]

Helper function for infinite radius computation with logistic transformation.

Used in normalizing constant computation via integration methods.

Parameters:
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:
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:
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 Reduced Decomposition algorithm.

lGN: float
Cgamma: float
pfqn_nrl(L, N, Z=None, alpha=None)[source]

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

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

Parameters:
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]

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

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

Parameters:
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:
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.

Matrix-Analytic Methods (line_solver.api.mam)

Matrix-Analytic Methods (MAM) for MAP/PH distributions.

Native Python implementations for analyzing Markovian Arrival Processes (MAPs), Phase-Type (PH) distributions, and related matrix-analytic methods.

Key algorithms:

map_piq: CTMC steady-state of MAP map_pie: Embedded DTMC steady-state map_lambda: Arrival rate computation map_mean, map_var, map_scv: Moment computations solver_mam_map_bmap_1: MAP/BMAP/1 queue solver using GI/M/1-type ETAQA solver_mam_bmap_map_1: BMAP/MAP/1 queue solver using M/G/1-type ETAQA

map_infgen(D0, D1)[source]

Compute the infinitesimal generator of a MAP.

The generator Q = D0 + D1 represents the underlying CTMC.

Parameters:

  • D0 (numpy.ndarray) – Hidden transition matrix (non-arrival transitions)

  • D1 (numpy.ndarray) – Visible transition matrix (arrival transitions)

Returns:

Infinitesimal generator matrix Q = D0 + D1

Return type:

numpy.ndarray

map_piq(D0, D1=None)[source]

Compute steady-state distribution of the underlying CTMC of a MAP.

Solves πQ = 0 where Q = D0 + D1 is the generator.

Parameters:

  • D0 (numpy.ndarray) – Hidden transition matrix, or stacked [D0, D1] if D1 is None

  • D1 (numpy.ndarray) – Visible transition matrix (optional)

Returns:

Steady-state probability vector π

Return type:

numpy.ndarray

map_pie(D0, D1=None)[source]

Compute equilibrium distribution of embedded DTMC.

The embedded DTMC has transition matrix P = (-D0)^{-1} * D1. Its steady-state is π_e = π * D1 / (π * D1 * e).

Parameters:

  • D0 (numpy.ndarray) – Hidden transition matrix, or stacked [D0, D1] if D1 is None

  • D1 (numpy.ndarray) – Visible transition matrix (optional)

Returns:

Equilibrium distribution of embedded DTMC

Return type:

numpy.ndarray

map_lambda(D0, D1=None)[source]

Compute the arrival rate (λ) of a MAP.

The arrival rate is λ = π * D1 * e where π is the steady-state and e is the column vector of ones.

Parameters:

  • D0 (numpy.ndarray) – Hidden transition matrix, or stacked [D0, D1] if D1 is None

  • D1 (numpy.ndarray) – Visible transition matrix (optional)

Returns:

Arrival rate λ

Return type:

float

map_mean(D0, D1=None)[source]

Compute mean inter-arrival time of a MAP.

The mean is 1/λ where λ is the arrival rate.

Parameters:

  • D0 (numpy.ndarray) – Hidden transition matrix, or stacked [D0, D1] if D1 is None

  • D1 (numpy.ndarray) – Visible transition matrix (optional)

Returns:

Mean inter-arrival time

Return type:

float

map_var(D0, D1=None)[source]

Compute variance of inter-arrival times of a MAP.

Var[X] = E[X²] - E[X]²

Uses map_moment for consistency with Kotlin implementation.

Parameters:

  • D0 (numpy.ndarray) – Hidden transition matrix, or stacked [D0, D1] if D1 is None

  • D1 (numpy.ndarray) – Visible transition matrix (optional)

Returns:

Variance of inter-arrival times

Return type:

float

map_scv(D0, D1=None)[source]

Compute squared coefficient of variation (SCV) of a MAP.

SCV = Var[X] / E[X]² = (E[X²] - E[X]²) / E[X]²

Parameters:

  • D0 (numpy.ndarray) – Hidden transition matrix, or stacked [D0, D1] if D1 is None

  • D1 (numpy.ndarray) – Visible transition matrix (optional)

Returns:

Squared coefficient of variation

Return type:

float

map_moment(D0, D1, k)[source]

Compute the k-th moment of inter-arrival time distribution.

E[X^k] = k! * π_e * (-D0)^{-k} * e

where π_e is the embedded DTMC steady-state.

Parameters:
Returns:

k-th raw moment

Return type:

float

map_scale(D0, D1, factor)[source]

Scale a MAP by a given factor.

Scaling changes the time scale: D0’ = factor * D0, D1’ = factor * D1.

Parameters:
Returns:

Tuple of scaled (D0’, D1’)

Return type:

Tuple[numpy.ndarray, numpy.ndarray]

map_normalize(D0, D1)[source]

Normalize a MAP to have unit mean inter-arrival time.

Parameters:
Returns:

Tuple of normalized (D0’, D1’)

Return type:

Tuple[numpy.ndarray, numpy.ndarray]

map_isfeasible(D0, D1, tolerance=1e-10)[source]

Check if (D0, D1) form a valid MAP.

A valid MAP requires: - D0 has non-positive diagonal and non-negative off-diagonal - D1 has non-negative elements - D0 + D1 is a valid generator (row sums = 0)

Parameters:
Returns:

True if valid MAP

Return type:

bool

exp_map(lambda_rate)[source]

Create a MAP representation of an exponential distribution.

Parameters:

lambda_rate (float) – Rate parameter (λ > 0)

Returns:

Tuple of (D0, D1) matrices representing Exp(λ)

Return type:

Tuple[numpy.ndarray, numpy.ndarray]

erlang_map(k, lambda_rate)[source]

Create a MAP representation of an Erlang-k distribution.

Parameters:

  • k (int) – Number of phases (k >= 1)

  • lambda_rate (float) – Overall rate parameter

Returns:

Tuple of (D0, D1) matrices representing Erlang(k, k*λ)

Return type:

Tuple[numpy.ndarray, numpy.ndarray]

hyperexp_map(probs, rates)[source]

Create a MAP representation of a hyperexponential distribution.

Parameters:
Returns:

Tuple of (D0, D1) matrices representing hyperexponential

Return type:

Tuple[numpy.ndarray, numpy.ndarray]

map_exponential(mean)[source]

Create a MAP representation of an exponential distribution (MATLAB-style).

This is a MATLAB-compatible wrapper that takes mean instead of rate.

Parameters:

mean (float) – Mean inter-arrival time (= 1/λ)

Returns:

Tuple of (D0, D1) matrices representing Exp(1/mean)

Return type:

Tuple[numpy.ndarray, numpy.ndarray]

Examples

>>> D0, D1 = map_exponential(2)  # Poisson process with rate λ=0.5
map_erlang(mean, k)[source]

Create a MAP representation of an Erlang-k distribution (MATLAB-style).

This is a MATLAB-compatible wrapper that takes mean as first argument.

Parameters:

  • mean (float) – Mean inter-arrival time

  • k (int) – Number of phases

Returns:

Tuple of (D0, D1) matrices representing Erlang-k with given mean

Return type:

Tuple[numpy.ndarray, numpy.ndarray]

Examples

>>> D0, D1 = map_erlang(2, 3)  # Erlang-3 with mean 2
map_hyperexp(probs, means)[source]

Create a MAP representation of a hyperexponential distribution (MATLAB-style).

This is a MATLAB-compatible wrapper that takes means instead of rates.

Parameters:

  • probs (numpy.ndarray) – Probability vector for choosing each phase

  • means (numpy.ndarray) – Mean values for each phase (1/rates)

Returns:

Tuple of (D0, D1) matrices representing hyperexponential

Return type:

Tuple[numpy.ndarray, numpy.ndarray]

map_gamma(mean, n)[source]

Create a MAP approximation of a Gamma distribution.

Uses Erlang-n with matching mean as an approximation.

Parameters:

  • mean (float) – Mean of the distribution

  • n (int) – Number of phases (shape parameter)

Returns:

Tuple of (D0, D1) matrices

Return type:

Tuple[numpy.ndarray, numpy.ndarray]

map_sumind(maps)[source]

Compute the sum of independent MAPs.

Creates a MAP representing the sum (concatenation) of independent random variables represented by the input MAPs.

Parameters:

maps (list) – List of MAPs, each as (D0, D1) tuple

Returns:

Tuple of (D0, D1) matrices representing the sum

Return type:

Tuple[numpy.ndarray, numpy.ndarray]

Examples

>>> # Sum of exponential and Erlang-2
>>> MAP1 = map_exponential(1.0)
>>> MAP2 = map_erlang(1.0, 2)
>>> D0, D1 = map_sumind([MAP1, MAP2])
map_cdf(D0, D1, points)[source]

Compute cumulative distribution function of inter-arrival times.

F(t) = 1 - π_e * exp(D0*t) * e

Parameters:
Returns:

CDF values at specified points

Return type:

numpy.ndarray

Examples

>>> map_cdf(D0, D1, 1.0)  # Returns P(T <= 1)
>>> map_cdf(D0, D1, [1.0, 5.0])  # Returns [P(T<=1), P(T<=5)]
map_pdf(D0, D1, points)[source]

Compute probability density function of inter-arrival times.

f(t) = π_e * exp(D0*t) * (-D0) * e

Parameters:
Returns:

PDF values at specified points

Return type:

numpy.ndarray

Examples

>>> map_pdf(D0, D1, [0.5, 1.0, 2.0])
map_sample(D0, D1, n_samples, rng=None)[source]

Generate random samples from a MAP distribution.

Parameters:
Returns:

Array of inter-arrival times

Return type:

numpy.ndarray

map_skew(D0, D1)[source]

Compute skewness of inter-arrival times.

Parameters:
Returns:

Skewness of inter-arrival times

Return type:

float

map_kurt(D0, D1)[source]

Compute kurtosis of inter-arrival times.

Parameters:
Returns:

Kurtosis of inter-arrival times

Return type:

float

map_acf(D0, D1, lags=1)[source]

Compute autocorrelation coefficients of inter-arrival times.

Parameters:
Returns:

Array of autocorrelation coefficients at specified lags

Return type:

numpy.ndarray

Examples

>>> map_acf(D0, D1)  # lag-1 autocorrelation
>>> map_acf(D0, D1, np.arange(1, 11))  # first 10 autocorrelations
map_acfc(D0, D1, kset, u)[source]

Compute autocorrelation of counting process at given lags.

Parameters:
Returns:

Autocorrelation coefficients at specified lags

Return type:

numpy.ndarray

map_idc(D0, D1)[source]

Compute the asymptotic index of dispersion.

I = SCV * (1 + 2 * sum_{k=1}^{inf} rho_k)

where SCV is the squared coefficient of variation and rho_k is the lag-k autocorrelation coefficient.

Parameters:
Returns:

Asymptotic index of dispersion

Return type:

float

map_count_mean(D0, D1, t)[source]

Compute mean of counting process at resolution t.

Parameters:
Returns:

Mean arrivals in (0, t]

Return type:

numpy.ndarray

map_count_var(D0, D1, t)[source]

Compute variance of counting process at resolution t.

Parameters:
Returns:

Variance of arrivals in (0, t]

Return type:

numpy.ndarray

Reference:

He and Neuts, “Markov chains with marked transitions”, 1998

map_varcount(D0, D1, tset)[source]

Compute variance of counting process (alternative implementation).

Parameters:
Returns:

Variance at each time point

Return type:

numpy.ndarray

map_count_moment(D0, D1, t, orders)[source]

Compute power moments of counts at resolution t.

Uses numerical differentiation of the moment generating function.

Parameters:
Returns:

Power moments of counts

Return type:

numpy.ndarray

map_mmpp2(mean, scv, skew=-1, acf1=-1)[source]

Fit an MMPP(2) as a MAP.

Parameters:

  • mean (float) – Mean inter-arrival time

  • scv (float) – Squared coefficient of variation

  • skew (float) – Skewness (-1 for automatic minimization)

  • acf1 (float) – Lag-1 autocorrelation (-1 for maximum feasible)

Returns:

Tuple of (D0, D1) matrices

Return type:

Tuple[numpy.ndarray, numpy.ndarray]

Examples

>>> D0, D1 = map_mmpp2(1, 2, -1, 0.2)  # Minimal skewness, ACF=0.2
>>> D0, D1 = map_mmpp2(1, 2, -1, -1)   # Minimal skewness, max ACF
map_gamma2(D0, D1)[source]

Compute the second largest eigenvalue of embedded DTMC.

This is the autocorrelation decay rate.

Parameters:
Returns:

Second largest eigenvalue (gamma_2)

Return type:

float

map_rand(k=2)[source]

Generate a random MAP of order k.

Parameters:

k (int) – Order of the MAP (default: 2)

Returns:

Tuple of (D0, D1) matrices

Return type:

Tuple[numpy.ndarray, numpy.ndarray]

map_randn(k, mu=(1.0, 1.0), sigma=(0.5, 0.5))[source]

Generate a random MAP with normally distributed elements.

Parameters:
Returns:

Tuple of (D0, D1) matrices

Return type:

Tuple[numpy.ndarray, numpy.ndarray]

map_renewal(D0, D1)[source]

Remove all correlations from a MAP, creating a renewal process.

Parameters:
Returns:

Renewal MAP with same CDF but no correlations

Return type:

Tuple[numpy.ndarray, numpy.ndarray]

map_embedded(D0, D1)[source]

Compute embedded discrete-time transition matrix.

P = (-D0)^{-1} * D1

Parameters:
Returns:

Transition matrix of embedded DTMC

Return type:

numpy.ndarray

map_sum(D0, D1, n)[source]

Create MAP for sum of n IID random variables.

Parameters:
Returns:

Tuple of (D0_new, D1_new) for the sum

Return type:

Tuple[numpy.ndarray, numpy.ndarray]

map_super(D0_a, D1_a, D0_b, D1_b)[source]

Create superposition of two MAPs.

Parameters:
Returns:

Superposed MAP (D0, D1)

Return type:

Tuple[numpy.ndarray, numpy.ndarray]

map_mixture(alpha, maps)[source]

Create probabilistic mixture of MAPs.

Parameters:

  • alpha (numpy.ndarray) – Probability vector for choosing each MAP

  • maps (list) – List of MAPs as (D0, D1) tuples

Returns:

Mixture MAP (D0, D1)

Return type:

Tuple[numpy.ndarray, numpy.ndarray]

map_max(D0_a, D1_a, D0_b, D1_b)[source]

Create MAP for max(X, Y) where X ~ MAP_A and Y ~ MAP_B.

Parameters:
Returns:

MAP for the maximum (D0, D1)

Return type:

Tuple[numpy.ndarray, numpy.ndarray]

map_timereverse(D0, D1)[source]

Compute time-reversed MAP.

Parameters:
Returns:

Time-reversed MAP (D0_r, D1_r)

Return type:

Tuple[numpy.ndarray, numpy.ndarray]

map_mark(D0, D1, prob)[source]

Mark arrivals from a MAP according to given probabilities.

Parameters:
Returns:

MMAP as list [D0, D1_class1, D1_class2, …]

Return type:

list

map_stochcomp(D0, D1, retain_idx)[source]

Stochastic complementation to reduce MAP order.

Parameters:
Returns:

Reduced MAP (D0_new, D1_new)

Return type:

Tuple[numpy.ndarray, numpy.ndarray]

map_kpc(maps)[source]

Kronecker product composition of MAPs.

Parameters:

maps (list) – List of MAPs as (D0, D1) tuples

Returns:

Composed MAP (D0, D1)

Return type:

Tuple[numpy.ndarray, numpy.ndarray]

map_bernstein(f, n=20)[source]

Convert distribution to MAP via Bernstein approximation.

Parameters:

  • f – PDF function handle

  • n (int) – Number of phases (default: 20)

Returns:

MAP representation (D0, D1)

Return type:

Tuple[numpy.ndarray, numpy.ndarray]

map_pntiter(D0, D1, na, t, M=None)[source]

Compute probability of na arrivals in interval [0, t].

Uses iterative bisection method (Neuts and Li).

Parameters:

  • D0 (numpy.ndarray) – Hidden transition matrix

  • D1 (numpy.ndarray) – Visible transition matrix

  • na (int) – Number of arrivals

  • t (float) – Time interval length

  • M (int | None) – Bisection depth (auto-computed if None)

Returns:

Matrix of probabilities

Return type:

numpy.ndarray

map_pntquad(D0, D1, na, t)[source]

Compute probability of na arrivals using ODE solver.

Parameters:
Returns:

Matrix P_n(t)

Return type:

numpy.ndarray

map2_fit(e1, e2, e3=-1.0, g2=0.0)[source]

Fit a MAP(2) distribution to moments and autocorrelation.

Based on: A. Heindl, G. Horvath, K. Gross “Explicit inverse characterization of acyclic MAPs of second order”

Parameters:

  • e1 (float) – First moment E[X]

  • e2 (float) – Second moment E[X^2]

  • e3 (float) – Third moment E[X^3], or -1 for automatic selection, -2 for minimum, -3 for maximum, or negative fraction for interpolation

  • g2 (float) – Autocorrelation decay rate (gamma_2)

Returns:

  • MAP: Tuple (D0, D1) if successful, None if failed

  • error_code: 0 for success, >0 for various errors

Return type:

Tuple of (MAP, error_code) where

Error codes:

0: Success 10: Mean out of bounds 20: Correlated exponential 30: h2 out of bounds 40: h3 out of bounds 51-54: g2 out of bounds

map_joint(D0, D1, a, i)[source]

Compute joint moments of a MAP.

E[(X_{a1})^{i1} * (X_{a1+a2})^{i2} * …]

Parameters:
Returns:

Joint moment

Return type:

float

map_issym(D0, D1=None)[source]

Check if MAP contains symbolic elements.

In Python, this always returns False as we use numpy arrays.

Parameters:
Returns:

False (Python arrays are always numeric)

Return type:

bool

map_feastol()[source]

Get tolerance for feasibility checks.

Returns:

Tolerance exponent (10^(-feastol))

Return type:

int

map_feasblock(E1, E2, E3, G2, opt='')[source]

Fit the most similar feasible MAP(2).

Parameters:

  • E1 (float) – First moment (mean)

  • E2 (float) – Second moment (or SCV if opt=’scv’)

  • E3 (float) – Third moment

  • G2 (float) – Autocorrelation decay rate

  • opt (str) – ‘scv’ if E2 is SCV instead of second moment

Returns:

Feasible MAP (D0, D1)

Return type:

Tuple[numpy.ndarray, numpy.ndarray]

map_block(E1, E2, E3, G2, opt='')[source]

Construct a MAP(2) from moments and autocorrelation.

Parameters:

  • E1 (float) – First moment

  • E2 (float) – Second moment (or SCV if opt=’scv’)

  • E3 (float) – Third moment

  • G2 (float) – Autocorrelation decay rate

  • opt (str) – ‘scv’ if E2 is SCV

Returns:

MAP (D0, D1)

Return type:

Tuple[numpy.ndarray, numpy.ndarray]

map_largemap()[source]

Get threshold for “large” MAP where exact computation is expensive.

Returns:

Return type:

Order threshold (default

mmap_infgen(D0, D_list)[source]

Compute the infinitesimal generator of an MMAP.

Parameters:

  • D0 (numpy.ndarray) – Hidden transition matrix (non-arrival transitions)

  • D_list (List[numpy.ndarray]) – List of arrival marking matrices [D1, D2, …, Dm]

Returns:

Infinitesimal generator Q = D0 + D1 + D2 + … + Dm

Return type:

numpy.ndarray

mmap_normalize(D0, D_list)[source]

Normalize an MMAP to ensure -diag(D0) = rowsums(D1 + D2 + … + Dm).

This ensures D0 represents true hidden transitions (negative diagonal).

Parameters:
Returns:

Tuple of (normalized_D0, normalized_D_list)

Return type:

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

mmap_super_safe(mmap_list, maxorder=1000, method='default')[source]

Safely superpose multiple MMAPs into a single MMAP.

Combines multiple MMAPs while ensuring the order of the resulting MMAP does not exceed maxorder. If combining would exceed maxorder, alternative methods are used (fitting to exponential or simplified MAP).

Parameters:

  • mmap_list (List[Tuple[numpy.ndarray, List[numpy.ndarray]]]) – List of (D0, D_list) tuples, one per MMAP

  • maxorder (int) – Maximum allowed order (state space size) for result

  • method (str) – Combination method (“default” or “match”)

Returns:

Tuple of (D0_super, D_list_super) for the superposed MMAP

Return type:

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

Algorithm:

  1. Sort MMAPs by squared coefficient of variation (SCV) of unmarked process

  2. Iteratively combine: start with simplest, add others in order

  3. If Kronecker product would exceed maxorder, use simpler approximations

mmap_mark(D0, D_list, prob)[source]

Reclassify arrivals in an MMAP according to a probability matrix.

Converts an MMAP with K classes to one with R classes based on a KxR probability matrix that describes how arrivals are reclassified.

Parameters:

  • D0 (numpy.ndarray) – Hidden transition matrix

  • D_list (List[numpy.ndarray]) – List of K arrival marking matrices [D1, D2, …, DK]

  • prob (numpy.ndarray) – K x R probability matrix where prob[k,r] is the probability that a class-k arrival in the original MMAP becomes a class-r arrival in the new MMAP

Returns:

Tuple of (D0_new, D_list_new) where D_list_new has R matrices

Return type:

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

Example

If prob = [[0.8, 0.2], [0.3, 0.7]], a 2-class MMAP becomes a 2-class MMAP where 80% of class-1 arrivals are now class-1, 20% become class-2, etc.

mmap_scale(D0, D_list, M, max_iter=30)[source]

Scale the mean inter-arrival times of an MMAP.

Adjusts all matrices to achieve specified mean inter-arrival times.

Parameters:

  • D0 (numpy.ndarray) – Hidden transition matrix

  • D_list (List[numpy.ndarray]) – List of arrival marking matrices

  • M (float | numpy.ndarray) – Desired mean inter-arrival times. Can be: - Single float: uniform scaling of all classes - Array of K values: individual scaling per class

  • max_iter (int) – Maximum iterations for optimization (only for vector M)

Returns:

Tuple of (D0_scaled, D_list_scaled)

Return type:

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

Algorithm:

For single M: Scale all matrices uniformly by ratio = M_old / M For vector M: Use iterative coordinate descent to find scaling factors

mmap_hide(D0, D_list, types)[source]

Hide (remove) specified arrival classes from an MMAP.

The hidden classes are set to zero matrices, effectively removing them from observation while maintaining the underlying stochastic process.

Parameters:
Returns:

Tuple of (D0_hidden, D_list_hidden) where specified classes are zero

Return type:

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

Example

mmap_hide(D0, [D1, D2, D3], types=[1]) hides class 2 (index 1)

mmap_compress(D0, D_list, method='default', target_order=None)[source]

Compress an MMAP using various approximation methods.

Reduces the state space of an MMAP while preserving key statistical properties (moments, inter-arrival time distribution, etc.).

Parameters:

  • D0 (numpy.ndarray) – Hidden transition matrix

  • D_list (List[numpy.ndarray]) – List of arrival marking matrices

  • method (str) – Compression method. Options: - “default”, “mixture”: Mixture fitting (order 1) - “exponential”: Fit to single-state Poisson

  • target_order (int | None) – Target state space size (overrides method if specified)

Returns:

Tuple of (D0_compressed, D_list_compressed)

Return type:

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

Notes

The “default” mixture method fits each class independently to MMPP(2) if possible, then combines them. This is fast and often effective.

More advanced methods (m3pp, mamap2) can be added as extensions.

mmap_exponential(lambda_rate, nclasses=None)[source]

Create an exponential MMAP (Poisson arrival process).

Parameters:

  • lambda_rate (float | numpy.ndarray) – Arrival rate (scalar) or rate per class (vector)

  • nclasses (int | None) – Number of classes (if lambda_rate is scalar)

Returns:

Tuple of (D0, D_list) representing exponential MMAP

Return type:

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

mmap_issym(mmap)[source]

Check if an MMAP contains symbolic elements.

In Python/NumPy, we don’t have built-in symbolic support like MATLAB, so this always returns False for numpy arrays.

Parameters:

mmap (List[numpy.ndarray]) – List of matrices [D0, D1, D2, …, Dm]

Returns:

False (symbolic computation not supported in numpy)

Return type:

bool

mmap_isfeasible(mmap, tol=None)[source]

Check whether an MMAP is feasible up to the given tolerance.

Checks: - Elements are real (no imaginary parts) - D0 + D1 rows sum to zero (generator property) - Diagonal of D0 is negative - Off-diagonal of D0 is non-negative - All D1c matrices are non-negative - D1 = D11 + D12 + … + D1C

Parameters:

  • mmap (List[numpy.ndarray]) – List of matrices [D0, D1, D2, …, Dm]

  • tol (float) – Tolerance for feasibility checks (default: 1e-6)

Returns:

True if MMAP is feasible, False otherwise

Return type:

bool

mmap_lambda(mmap)[source]

Compute the arrival rate of each class in an MMAP.

Parameters:

mmap (List[numpy.ndarray]) – List of matrices [D0, D1, D2, …, Dm]

Returns:

Array of arrival rates, one per class

Return type:

numpy.ndarray

mmap_count_lambda(mmap)[source]

Compute the arrival rate of the counting process for each class.

The rate for class k is λ_k = π * D_k * e where π is the steady-state of the underlying CTMC.

Parameters:

mmap (List[numpy.ndarray]) – List of matrices [D0, D1, D2, …, Dm]

Returns:

Array of arrival rates, one per class

Return type:

numpy.ndarray

mmap_pie(mmap)[source]

Compute the stationary probability of the DTMC embedded at restart instants after an arrival of each class.

Parameters:

mmap (List[numpy.ndarray]) – List of matrices [D0, D1, D2, …, Dm]

Returns:

Matrix of shape (m, n) where row c is the steady-state for class c

Return type:

numpy.ndarray

mmap_pc(mmap)[source]

Compute the arrival probabilities of each class for the given MMAP.

Parameters:

mmap (List[numpy.ndarray]) – List of matrices [D0, D1, D2, …, Dm]

Returns:

Array of arrival probabilities, one per class

Return type:

numpy.ndarray

mmap_embedded(mmap)[source]

Compute the embedded DTMC transition matrices for each class.

For class k, the embedded matrix is Pc_k = (-D0)^{-1} * D_{2+k}

Parameters:

mmap (List[numpy.ndarray]) – List of matrices [D0, D1, D2, …, Dm]

Returns:

List of embedded DTMC transition matrices, one per class

Return type:

List[numpy.ndarray]

mmap_sample(mmap, n_samples, pi=None, seed=None)[source]

Generate samples from a Marked MAP.

Parameters:

  • mmap (List[numpy.ndarray]) – List of matrices [D0, D1, D2, …, Dm]

  • n_samples (int) – Number of random samples to collect

  • pi (numpy.ndarray) – Initial probability (if None, uses steady-state)

  • seed (int) – Random seed for reproducibility

Returns:

  • T: Array of inter-arrival times

  • A: Array of class labels (1-indexed for MATLAB compatibility)

Return type:

Tuple of (T, A) where

mmap_rand(order, classes)[source]

Generate a random MMAP with given order and number of classes.

Parameters:

  • order (int) – Number of states

  • classes (int) – Number of arrival classes

Returns:

List of matrices [D0, D1, D21, D22, …, D2K]

Return type:

List[numpy.ndarray]

mmap_timereverse(mmap)[source]

Compute the time-reversed MMAP.

Parameters:

mmap (List[numpy.ndarray]) – List of matrices [D0, D1, D2, …, Dm]

Returns:

Time-reversed MMAP

Return type:

List[numpy.ndarray]

mmap_sum(mmap, n)[source]

Create an MMAP representing the sum of n identical MMAPs (n-fold sum).

This creates an MMAP whose inter-arrival times are distributed as the sum of n inter-arrival times from the original MMAP.

Parameters:
Returns:

Summed MMAP

Return type:

List[numpy.ndarray]

mmap_super(mmap_a, mmap_b=None, opt='default')[source]

Superpose two or more MMAPs.

Parameters:

  • mmap_a (List[numpy.ndarray]) – First MMAP (or list of MMAPs if mmap_b is None)

  • mmap_b (List[numpy.ndarray]) – Second MMAP (optional)

  • opt (str) – “default” - each class in MMAPa and MMAPb is distinct in result “match” - class c in both maps is mapped to class c in result

Returns:

Superposed MMAP

Return type:

List[numpy.ndarray]

mmap_mixture(alpha, maps)[source]

Create a probabilistic mixture of MAPs.

Each MAP in the list is selected with probability alpha[i].

Parameters:

  • alpha (numpy.ndarray) – Array of mixing probabilities

  • maps (List) – List of MAPs (each MAP is [D0, D1])

Returns:

Mixed MMAP with len(maps) classes

Return type:

List[numpy.ndarray]

mmap_max(mmap_a, mmap_b, k)[source]

Create MMAP for the maximum of arrivals from two synchronized MMAPs.

This models a synchronization queue where arrivals from both sources must be paired before release.

Parameters:
Returns:

MMAP for the synchronized (maximum) process

Return type:

List[numpy.ndarray]

mmap_maps(mmap)[source]

Extract K MAPs, one for each class of the MMAP[K] process.

For class k, the MAP is {D0 + D1 - D_{2+k}, D_{2+k}}

Parameters:

mmap (List[numpy.ndarray]) – List of matrices [D0, D1, D21, D22, …, D2K]

Returns:

List of MAPs, one per class

Return type:

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

mmap_count_mean(mmap, t)[source]

Compute the mean of the counting process at resolution t.

Parameters:

  • mmap (List[numpy.ndarray]) – List of matrices [D0, D1, D2, …, Dm]

  • t (float) – Time period for counting

Returns:

Array of mean counts, one per class

Return type:

numpy.ndarray

mmap_count_var(mmap, t)[source]

Compute the variance of the counting process at resolution t.

Parameters:

  • mmap (List[numpy.ndarray]) – List of matrices [D0, D1, D2, …, Dm]

  • t (float) – Time period for counting

Returns:

Array of variances, one per class

Return type:

numpy.ndarray

mmap_count_idc(mmap, t)[source]

Compute the per-class Index of Dispersion of Counts at resolution t.

IDC = Var[N(t)] / E[N(t)]

Parameters:
Returns:

Array of IDC values, one per class

Return type:

numpy.ndarray

mmap_count_mcov(mmap, t)[source]

Compute the count covariance between each pair of classes at time scale t.

Parameters:
Returns:

m x m covariance matrix

Return type:

numpy.ndarray

mmap_idc(mmap)[source]

Compute the asymptotic Index of Dispersion of Counts for each class.

This is the limit of IDC(t) as t -> infinity.

Parameters:

mmap (List[numpy.ndarray]) – List of matrices [D0, D1, D2, …, Dm]

Returns:

Array of asymptotic IDC values, one per class

Return type:

numpy.ndarray

mmap_sigma(mmap)[source]

Compute one-step class transition probabilities.

p_{i,j} = P(C_k = j | C_{k-1} = i)

Parameters:

mmap (List[numpy.ndarray]) – List of matrices [D0, D1, D2, …, Dm]

Returns:

C x C matrix of transition probabilities

Return type:

numpy.ndarray

mmap_sigma2(mmap)[source]

Compute two-step class transition probabilities.

p_{i,j,h} = P(C_k = h | C_{k-1} = j, C_{k-2} = i)

Parameters:

mmap (List[numpy.ndarray]) – List of matrices [D0, D1, D2, …, Dm]

Returns:

C x C x C 3D array of transition probabilities

Return type:

numpy.ndarray

mmap_forward_moment(mmap, orders, norm=True)[source]

Compute the theoretical forward moments of an MMAP.

Forward moments are E[X^k | previous arrival was class c].

Parameters:

  • mmap (List[numpy.ndarray]) – List of matrices [D0, D1, D2, …, Dm]

  • orders (numpy.ndarray) – Array of moment orders to compute

  • norm (bool) – If True, normalize by class probability

Returns:

Matrix of shape (C, len(orders)) where element (c, k) is the order-k forward moment for class c

Return type:

numpy.ndarray

mmap_backward_moment(mmap, orders, norm=True)[source]

Compute the theoretical backward moments of an MMAP.

Backward moments are E[X^k | next arrival will be class c].

Parameters:

  • mmap (List[numpy.ndarray]) – List of matrices [D0, D1, D2, …, Dm]

  • orders (numpy.ndarray) – Array of moment orders to compute

  • norm (bool) – If True, normalize by class probability

Returns:

Matrix of shape (C, len(orders)) where element (c, k) is the order-k backward moment for class c

Return type:

numpy.ndarray

mmap_cross_moment(mmap, k)[source]

Compute the k-th order moment of inter-arrival times between class pairs.

Computes E[X^k | previous = class i, next = class j] for all i, j.

Parameters:

  • mmap (List[numpy.ndarray]) – List of matrices [D0, D1, D2, …, Dm]

  • k (int) – Order of the moment

Returns:

C x C matrix where element (i, j) is the k-th moment for transitions i->j

Return type:

numpy.ndarray

mmap_modulate(P, HT, MMAP)[source]

Modulate MMAPs in continuous time according to holding time distributions.

Creates an MMAP that switches between different MMAPs according to a continuous-time Markov modulated process.

Parameters:

  • P (numpy.ndarray) – J x J transition probability matrix between states

  • HT (List) – List of phase-type distributions for holding time in each state

  • MMAP (List) – List of MMAPs, one for each state

Returns:

Modulated MMAP

Return type:

List[numpy.ndarray]

ldqbd(Q0, Q1, Q2, options=None)[source]

Solve a level-dependent QBD process.

Parameters:

  • Q0 (List[numpy.ndarray]) – List of upward transition matrices [Q0^(0), Q0^(1), …, Q0^(N-1)] Q0^(n) has shape (states_n, states_{n+1})

  • Q1 (List[numpy.ndarray]) – List of local transition matrices [Q1^(0), Q1^(1), …, Q1^(N)] Q1^(n) has shape (states_n, states_n)

  • Q2 (List[numpy.ndarray]) – List of downward transition matrices [Q2^(1), Q2^(2), …, Q2^(N)] Q2^(n) has shape (states_n, states_{n-1})

  • options (LdqbdOptions | None) – LdqbdOptions instance (uses defaults if None)

Returns:

LdqbdResult containing rate matrices R and stationary distribution pi

Return type:

LdqbdResult

Algorithm:

  1. Compute rate matrices backward using continued fraction recursion

  2. Compute stationary distribution forward using R matrices

class LdqbdResult(R, pi)[source]

Bases: object

Result of LDQBD solver.

R

List of rate matrices R^(1), R^(2), …, R^(N)

Type:

List[numpy.ndarray]

pi

Stationary distribution vector [pi_0, pi_1, …, pi_N]

Type:

numpy.ndarray

R: List[numpy.ndarray]
pi: numpy.ndarray
class LdqbdOptions(epsilon=1e-10, max_iter=1000, verbose=False)[source]

Bases: object

Options for LDQBD solver.

epsilon

Convergence tolerance (default: 1e-10)

Type:

float

max_iter

Maximum number of iterations (default: 1000)

Type:

int

verbose

Print debug information (default: False)

Type:

bool

epsilon: float = 1e-10
max_iter: int = 1000
verbose: bool = False
class MAPBMAP1Result(mean_queue_length, utilization, mean_response_time, throughput, pi, R, mean_batch_size)[source]

Bases: object

Result of MAP/BMAP/1 queue analysis.

mean_queue_length: float

Mean queue length E[N]

utilization: float

Server utilization rho

mean_response_time: float

Mean response time E[R]

throughput: float

Throughput (arrival rate)

pi: numpy.ndarray

Aggregated stationary probabilities [pi0, pi1, piStar]

R: numpy.ndarray

R matrix

mean_batch_size: float

Mean batch size of service

solver_mam_map_bmap_1(C0, C1, D)[source]

Solve a MAP/BMAP/1 queue using GI/M/1-type matrix-analytic methods.

The MAP is specified by matrices (C0, C1) where: - C0: transitions without arrivals - C1: transitions triggering arrivals

The BMAP for service is specified by matrices {D0, D1, D2, …, DK} where: - D0: transitions without service completions - Dk: transitions triggering batch service of k customers (k >= 1)

The GI/M/1-type structure for MAP/BMAP/1 is: ```

B1 A0 0 0 … B2 A1 A0 0 …

Q = B3 A2 A1 A0 …

…

```

Where:

A0 = C1 otimes I_ms (MAP arrival, level +1) A1 = C0 otimes I_ms + I_ma otimes D0 (phase changes, level 0) A_{k+1} = I_ma otimes D_k (batch size k service, level -k)

Parameters:
Returns:

MAPBMAP1Result with performance metrics

Return type:

MAPBMAP1Result

class BMAPMAP1Result(mean_queue_length, utilization, mean_response_time, throughput, pi, G, mean_batch_size)[source]

Bases: object

Result of BMAP/MAP/1 queue analysis.

mean_queue_length: float

Mean queue length E[N]

utilization: float

Server utilization rho

mean_response_time: float

Mean response time E[R]

throughput: float

Throughput (total customer arrival rate)

pi: numpy.ndarray

Aggregated stationary probabilities [pi0, pi1, piStar]

G: numpy.ndarray

G matrix

mean_batch_size: float

Mean batch size of arrivals

solver_mam_bmap_map_1(D, S0, S1)[source]

Solve a BMAP/MAP/1 queue using M/G/1-type matrix-analytic methods.

The BMAP is specified by matrices {D0, D1, D2, …, DK} where: - D0: transitions without arrivals - Dk: transitions triggering batch size k arrivals (k >= 1)

The MAP for service is specified by matrices (S0, S1) where: - S0: transitions without service completions - S1: transitions triggering service completions

The M/G/1-type structure for BMAP/MAP/1 is: ```

A0 = I_ma otimes S1 (service completion, level -1) A1 = D0 otimes I_ms + I_ma otimes S0 (phase changes, level 0) A_{k+1} = D_k otimes I_ms (batch arrival size k, level +k)

```

Parameters:

  • D (List[numpy.ndarray]) – BMAP matrices as list [D0, D1, D2, …, DK]

  • S0 (numpy.ndarray) – MAP service matrix for transitions without service completions

  • S1 (numpy.ndarray) – MAP service matrix for service completions

Returns:

BMAPMAP1Result with performance metrics

Return type:

BMAPMAP1Result

class QBDResult(R, G=None, U=None, eta=None)[source]

Bases: object

Result of QBD analysis.

G: numpy.ndarray | None = None
U: numpy.ndarray | None = None
eta: float | None = None
R: numpy.ndarray
qbd_R(B, L, F, iter_max=100000, tol=1e-12)[source]

Compute QBD rate matrix R using successive substitutions.

Solves the matrix quadratic equation:

R^2 * A_{-1} + R * A_0 + A_1 = 0

where A_{-1} = B, A_0 = L, A_1 = F.

Parameters:

  • B (numpy.ndarray) – Backward transition block A_{-1}

  • L (numpy.ndarray) – Local transition block A_0

  • F (numpy.ndarray) – Forward transition block A_1

  • iter_max (int) – Maximum iterations (default: 100000)

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

Returns:

Rate matrix R

Return type:

numpy.ndarray

References

Original MATLAB: matlab/src/api/mam/qbd_R.m

qbd_R_logred(B, L, F, iter_max=1000, tol=1e-14)[source]

Compute QBD rate matrix R using logarithmic reduction.

Uses the logarithmic reduction algorithm which has quadratic convergence compared to linear convergence of successive substitutions.

Parameters:

  • B (numpy.ndarray) – Backward transition block A_{-1}

  • L (numpy.ndarray) – Local transition block A_0

  • F (numpy.ndarray) – Forward transition block A_1

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

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

Returns:

Rate matrix R

Return type:

numpy.ndarray

References

Original MATLAB: matlab/src/api/mam/qbd_R_logred.m Latouche & Ramaswami, Ch. 8

qbd_rg(B, L, F, method='logred', iter_max=1000, tol=1e-14)[source]

Compute both R and G matrices for a QBD process.

G is the minimal non-negative solution to:

A_1 * G^2 + A_0 * G + A_{-1} = 0

R is the minimal non-negative solution to:

R^2 * A_{-1} + R * A_0 + A_1 = 0

Parameters:

  • B (numpy.ndarray) – Backward transition block A_{-1}

  • L (numpy.ndarray) – Local transition block A_0

  • F (numpy.ndarray) – Forward transition block A_1

  • method (str) – ‘logred’ or ‘successive’ (default: ‘logred’)

  • iter_max (int) – Maximum iterations

  • tol (float) – Convergence tolerance

Returns:

QBDResult with R, G, U, and eta (caudal characteristic)

Return type:

QBDResult

References

Original MATLAB: matlab/src/api/mam/qbd_rg.m

qbd_blocks_mapmap1(D0_arr, D1_arr, D0_srv, D1_srv)[source]

Construct QBD blocks for a MAP/MAP/1 queue.

Builds the backward (B), local (L), and forward (F) transition blocks for the QBD representation of a MAP/MAP/1 queue.

Parameters:
Returns:

Tuple of (B, L, F) QBD blocks

Return type:

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

References

Original MATLAB: matlab/src/api/mam/qbd_mapmap1.m

qbd_bmapbmap1(MAPa, pbatcha, MAPs)[source]

Compute QBD blocks for a BMAP/BMAP/1 queue.

Constructs the QBD (Quasi-Birth-Death) transition blocks for a BMAP/BMAP/1 queue with batch arrivals.

Parameters:
Returns:

A0: Local transition block A_1: Downward transition block A1_list: List of upward transition blocks for each batch size B0: Initial boundary local block B1_list: List of boundary upward blocks for each batch size

Return type:

Tuple of (A0, A_1, A1_list, B0, B1_list) where

References

Original MATLAB: matlab/src/api/mam/qbd_bmapbmap1.m

qbd_mapmap1(MAPa, MAPs, util=None)[source]

Analyze a MAP/MAP/1 queue using QBD methods.

Solves a MAP/MAP/1 queue using Quasi-Birth-Death process methods, computing throughput, queue length, utilization, and other metrics.

Parameters:
Returns:

Tuple of (XN, QN, UN, pqueue, R, eta, G, A_1, A0, A1, U, MAPs_scaled) where:

XN: System throughput QN: Mean queue length UN: Utilization pqueue: Queue length distribution R: Rate matrix R eta: Caudal characteristic (spectral radius of R) G: Rate matrix G A_1: Downward transition block A0: Local transition block A1: Upward transition block U: Matrix U MAPs_scaled: Scaled service process

Return type:

Tuple[float, float, float, numpy.ndarray, numpy.ndarray, float | None, numpy.ndarray, numpy.ndarray, numpy.ndarray, numpy.ndarray, numpy.ndarray, Tuple[numpy.ndarray, numpy.ndarray]]

References

Original MATLAB: matlab/src/api/mam/qbd_mapmap1.m

qbd_raprap1(RAPa, RAPs, util=None)[source]

Analyze a RAP/RAP/1 queue using QBD methods.

Solves a RAP/RAP/1 queue (Rational Arrival Process) using QBD methods, computing throughput, queue length, utilization, and other metrics.

Note: This is a simplified implementation. For full RAP analysis, specialized algorithms like those in SMCSolver are needed.

Parameters:
Returns:

XN: System throughput QN: Mean queue length UN: Utilization pqueue: Queue length distribution R: Rate matrix R eta: Caudal characteristic G: Rate matrix G B: Backward transition block L: Local transition block F: Forward transition block

Return type:

Tuple of (XN, QN, UN, pqueue, R, eta, G, B, L, F) where

References

Original MATLAB: matlab/src/api/mam/qbd_raprap1.m

qbd_setupdelayoff(lambda_val, mu, alpharate, alphascv, betarate, betascv)[source]

Analyze queue with setup delay and turn-off phases.

Performs queue-length analysis for a queueing system with setup delay (warm-up) and turn-off periods using QBD methods.

The system operates as follows: 1. When empty and job arrives, server enters setup phase 2. After setup, server becomes active and serves jobs 3. When queue empties, server enters turn-off phase 4. After turn-off, server becomes idle

Parameters:

  • lambda_val (float) – Arrival rate

  • mu (float) – Service rate

  • alpharate (float) – Rate of setup delay phase

  • alphascv (float) – Squared coefficient of variation for setup delay

  • betarate (float) – Rate of turn-off phase

  • betascv (float) – Squared coefficient of variation for turn-off period

Returns:

Average queue length QN

Return type:

float

References

Original MATLAB: matlab/src/api/mam/qbd_setupdelayoff.m

map_ccdf_derivative(MAP, i)[source]

Compute derivative at zero of a MAP’s complementary CDF.

Calculates the i-th derivative at t=0 of the complementary cumulative distribution function (CCDF) of the inter-arrival time distribution.

Formula: ν_i = π_e * D_0^i * e

where π_e is the embedded stationary vector and e is the column vector of ones.

Parameters:

  • MAP (List[numpy.ndarray]) – Markovian Arrival Process as [D0, D1]

  • i (int) – Order of the derivative

Returns:

Value of the i-th derivative at zero

Return type:

float

References

Original MATLAB: matlab/src/api/mam/map_ccdf_derivative.m A. Horvath et al., “A Joint Moments Based Analysis of Networks of MAP/MAP/1 Queues”

map_jointpdf_derivative(MAP, iset)[source]

Compute partial derivative at zero of a MAP’s joint PDF.

Calculates the partial derivative at t=0 of the joint probability density function of consecutive inter-arrival times.

For index set {i_1, i_2, …, i_k}: γ = π_e * D_0^{i_1} * D_1 * D_0^{i_2} * D_1 * … * D_0^{i_k} * D_1 * e

Parameters:
Returns:

Value of the partial derivative at zero

Return type:

float

References

Original MATLAB: matlab/src/api/mam/map_jointpdf_derivative.m A. Horvath et al., “A Joint Moments Based Analysis of Networks of MAP/MAP/1 Queues”

map_factorial_moment(MAP, k)[source]

Compute the k-th factorial moment of a MAP.

The k-th factorial moment is computed using derivatives of the inter-arrival time distribution.

Parameters:

  • MAP (List[numpy.ndarray]) – Markovian Arrival Process as [D0, D1]

  • k (int) – Order of the factorial moment

Returns:

k-th factorial moment

Return type:

float

References

Based on MAP moment formulas from matrix-analytic methods

map_joint_moment(MAP, k, l)[source]

Compute the (k,l)-th joint moment of consecutive inter-arrival times.

E[X_n^k * X_{n+1}^l] for a MAP with inter-arrival times X_n.

Parameters:

  • MAP (List[numpy.ndarray]) – Markovian Arrival Process as [D0, D1]

  • k (int) – Power for first inter-arrival time

  • l (int) – Power for second inter-arrival time

Returns:

Joint moment E[X_n^k * X_{n+1}^l]

Return type:

float

References

Based on joint moment formulas for MAPs

map_m1ps_cdf_respt(C, D, mu, x, epsilon=1e-11, epsilon_prime=1e-10)[source]

Compute complementary sojourn time CDF for MAP/M/1-PS queue.

Based on:

Masuyama, H., & Takine, T. (2003). Sojourn time distribution in a MAP/M/1 processor-sharing queue.

Parameters:

  • C (numpy.ndarray | list) – MAP C matrix (transitions without arrivals).

  • D (numpy.ndarray | list) – MAP D matrix (transitions with arrivals).

  • mu (float) – Service rate (must be > 0).

  • x (numpy.ndarray | list) – Time points for CDF evaluation.

  • epsilon (float) – Queue length truncation parameter.

  • epsilon_prime (float) – Uniformization truncation parameter.

Returns:

Complementary CDF values W_bar(x) = Pr[W > x] at each point in x.

Return type:

numpy.ndarray

Example

>>> # Exponential arrivals (Poisson process with rate lambda=0.5)
>>> C = np.array([[-0.5]])
>>> D = np.array([[0.5]])
>>> mu = 1.0
>>> x = np.array([0.0, 0.5, 1.0, 2.0, 5.0])
>>> cdf = map_m1ps_cdf_respt(C, D, mu, x)
map_compute_R(C, D, mu)[source]

Compute rate matrix R for MAP/M/1 queue.

Computes the minimal nonnegative solution of the matrix equation:

D + R(C - mu*I) + mu*R^2 = 0

This matrix is used in the analysis of MAP/M/1 queues based on quasi-birth-death processes.

Parameters:

  • C (numpy.ndarray | list) – M x M matrix governing MAP transitions without arrivals

  • D (numpy.ndarray | list) – M x M matrix governing MAP transitions with arrivals

  • mu (float) – Service rate (scalar)

Returns:

M x M rate matrix (minimal nonnegative solution)

Return type:

numpy.ndarray

References

Original MATLAB: matlab/src/api/map/map_compute_R.m Masuyama, H., & Takine, T. (2003). Sojourn time distribution in a MAP/M/1 processor-sharing queue.

map_m1ps_h_recursive(C, D, mu, N, K)[source]

Recursive computation of h_{n,k} coefficients for MAP/M/1-PS.

The h_{n,k} vectors satisfy the recursion (Theorem 1 in paper):

h_{n,0} = e (vector of ones), for n = 0, 1, … h_{n,k+1} = 1/(theta+mu) * [n*mu/(n+1) * h_{n-1,k} + (theta*I + C) * h_{n,k}

  • D * h_{n+1,k}]

where h_{-1,k} = 0 for all k

Parameters:

  • C (numpy.ndarray | list) – M x M matrix governing MAP transitions without arrivals

  • D (numpy.ndarray | list) – M x M matrix governing MAP transitions with arrivals

  • mu (float) – Service rate (scalar)

  • N (int) – Maximum value of n to compute (determines rows)

  • K (int) – Maximum value of k to compute (determines columns)

Returns:

2D list h[n][k] of (M,1) arrays containing h_{n,k}

Return type:

list

References

Original MATLAB: matlab/src/api/map/map_m1ps_h_recursive.m Masuyama, H., & Takine, T. (2003).

map_m1ps_sojourn(C, D, mu, x, epsilon=1e-11, epsilon_prime=1e-10)[source]

Compute sojourn time distribution in MAP/M/1-PS queue.

Alias for map_m1ps_cdf_respt.

Parameters:

  • C (numpy.ndarray | list) – MAP C matrix (transitions without arrivals)

  • D (numpy.ndarray | list) – MAP D matrix (transitions with arrivals)

  • mu (float) – Service rate (must be > 0)

  • x (numpy.ndarray | list) – Time points for CDF evaluation

  • epsilon (float) – Queue length truncation parameter

  • epsilon_prime (float) – Uniformization truncation parameter

Returns:

Complementary CDF values W_bar(x) = Pr[W > x]

Return type:

numpy.ndarray

References

Original MATLAB: matlab/src/api/map/map_m1ps_sojourn.m

mmdp_isfeasible(Q, R, tol=1e-10)[source]

Check if (Q, R) defines a valid MMDP.

Requirements: - Q must be a valid generator (square, row sums = 0, proper signs) - R must be diagonal with non-negative entries - Q and R must have compatible dimensions

Parameters:

  • Q (numpy.ndarray) – Generator matrix (n x n)

  • R (numpy.ndarray) – Rate matrix (n x n diagonal) or rate vector (n,)

  • tol (float) – Numerical tolerance for validation

Returns:

True if (Q, R) defines a valid MMDP, False otherwise

Return type:

bool

Examples

>>> Q = np.array([[-0.5, 0.5], [0.3, -0.3]])
>>> R = np.array([[2.0, 0], [0, 5.0]])
>>> mmdp_isfeasible(Q, R)
True

Markov Chain API (line_solver.api.mc)

Markov Chain analysis algorithms.

Native Python implementations for continuous-time and discrete-time Markov chain analysis.

Key algorithms:

ctmc_solve: CTMC steady-state distribution ctmc_transient: CTMC transient analysis ctmc_uniformization: Uniformization method for transient analysis ctmc_stochcomp: Stochastic complementation dtmc_solve: DTMC steady-state distribution

ctmc_solve(Q)[source]

Solve for steady-state probabilities of a CTMC.

Computes the stationary distribution π by solving πQ = 0 with normalization constraint Σπ = 1.

Handles reducible CTMCs by decomposing into strongly connected components and solving each separately.

Parameters:

Q (numpy.ndarray) – Infinitesimal generator matrix (row sums should be zero)

Returns:

Steady-state probability distribution (1D array)

Return type:

numpy.ndarray

ctmc_solve_reducible(Q)[source]

Solve reducible CTMCs by converting to DTMC via uniformization.

Parameters:

Q (numpy.ndarray) – Infinitesimal generator matrix (possibly reducible)

Returns:

Steady-state probability vector

Return type:

numpy.ndarray

ctmc_makeinfgen(Q)[source]

Convert a matrix into a valid infinitesimal generator for a CTMC.

An infinitesimal generator has: - Row sums equal to zero - Non-positive diagonal elements - Non-negative off-diagonal elements

Parameters:

Q (numpy.ndarray) – Candidate infinitesimal generator matrix

Returns:

Valid infinitesimal generator matrix with corrected diagonal

Return type:

numpy.ndarray

ctmc_transient(Q, initial_dist, time_points, method='expm')[source]

Compute transient probabilities of a CTMC.

Calculates time-dependent state probabilities π(t) for specified time points using matrix exponential methods.

Parameters:

  • Q (numpy.ndarray) – Infinitesimal generator matrix

  • initial_dist (numpy.ndarray) – Initial probability distribution π(0)

  • time_points (float | numpy.ndarray) – Array of time points to evaluate, or single time value

  • method (str) – ‘expm’ for matrix exponential, ‘ode’ for ODE solver

Returns:

Transient probabilities at each time point. Shape: (len(time_points), n) if multiple times, (n,) if single time

Return type:

numpy.ndarray

ctmc_uniformization(Q, lambda_rate=None)[source]

Uniformize CTMC generator matrix.

Converts CTMC to an equivalent uniformized discrete-time chain for numerical analysis and simulation purposes.

The uniformized DTMC has transition matrix P = I + Q/λ where λ is the uniformization rate (max exit rate).

Parameters:

  • Q (numpy.ndarray) – Infinitesimal generator matrix

  • lambda_rate (float | None) – Uniformization rate (optional, auto-computed if None)

Returns:

  • ‘P’: Uniformized transition matrix

  • ’lambda’: Uniformization rate

Return type:

dict containing

ctmc_randomization(Q, initial_dist, time_points, precision=1e-10)[source]

Compute CTMC transient probabilities using randomization.

Uses Jensen’s randomization method (uniformization) to compute transient probabilities by converting the CTMC to a uniformized DTMC.

This method is numerically stable and avoids matrix exponentials.

Parameters:

  • Q (numpy.ndarray) – Infinitesimal generator matrix

  • initial_dist (numpy.ndarray) – Initial probability distribution

  • time_points (numpy.ndarray) – Array of time points to evaluate

  • precision (float) – Numerical precision for truncation (Poisson tail)

Returns:

Transient probabilities at each time point

Return type:

numpy.ndarray

ctmc_stochcomp(Q, keep_states, eliminate_states=None)[source]

Compute stochastic complement of CTMC.

Reduces the CTMC by eliminating specified states while preserving the steady-state distribution restricted to the kept states.

Parameters:

  • Q (numpy.ndarray) – Infinitesimal generator matrix

  • keep_states (numpy.ndarray) – States to retain in reduced model (array of indices)

  • eliminate_states (numpy.ndarray | None) – States to eliminate (optional, inferred if None)

Returns:

  • ‘S’: Stochastic complement (reduced generator)

  • ’Q11’: Submatrix for kept states

  • ’Q12’: Transitions from kept to eliminated

  • ’Q21’: Transitions from eliminated to kept

  • ’Q22’: Submatrix for eliminated states

  • ’T’: Transient contribution matrix

Return type:

dict containing

ctmc_timereverse(Q, pi=None)[source]

Compute time-reversed CTMC generator.

The time-reversed generator Q* has elements: Q*_{ij} = π_j * Q_{ji} / π_i

Parameters:

  • Q (numpy.ndarray) – Original infinitesimal generator matrix

  • pi (numpy.ndarray | None) – Steady-state distribution (optional, computed if None)

Returns:

Time-reversed generator matrix

Return type:

numpy.ndarray

ctmc_rand(n, density=0.3, max_rate=10.0)[source]

Generate random CTMC generator matrix.

Parameters:

  • n (int) – Number of states

  • density (float) – Sparsity density (0 to 1, default 0.3)

  • max_rate (float) – Maximum transition rate (default 10.0)

Returns:

Random infinitesimal generator matrix

Return type:

numpy.ndarray

ctmc_simulate(Q, initial_state, max_time, max_events=10000, seed=None)[source]

Simulate CTMC sample path using Gillespie algorithm.

Generates a realization of the continuous-time Markov chain using the next-reaction method.

Parameters:

  • Q (numpy.ndarray) – Infinitesimal generator matrix

  • initial_state (int) – Starting state (integer index)

  • max_time (float) – Maximum simulation time

  • max_events (int) – Maximum number of transitions (default: 10000)

  • seed (int | None) – Random seed for reproducibility (optional)

Returns:

  • ‘states’: Array of visited states

  • ’times’: Array of transition times

  • ’sojourn_times’: Time spent in each state

Return type:

dict with

ctmc_isfeasible(Q, tolerance=1e-10)[source]

Check if matrix is a valid CTMC infinitesimal generator.

Validates: - Off-diagonal elements are non-negative - Row sums are zero - Diagonal elements are non-positive

Parameters:

  • Q (numpy.ndarray) – Candidate generator matrix

  • tolerance (float) – Numerical tolerance (default: 1e-10)

Returns:

True if matrix is valid CTMC generator

Return type:

bool

ctmc_ssg(sn, options=None)[source]

Generate complete CTMC state space for a queueing network.

Creates all possible network states including those not reachable from the initial state. For open classes, a cutoff parameter limits the maximum population to keep state space finite.

The state space is aggregated to show per-station-class job counts.

Parameters:

  • sn (Any) – NetworkStruct object (from getStruct())

  • options (Dict | None) – Solver options dict with fields: - cutoff: Population cutoff for open classes (required if open) - config.hide_immediate: Hide immediate transitions (default True)

Returns:

  • state_space: Complete state space matrix (rows=states, cols=state components)

  • state_space_aggr: Aggregated state space (rows=states, cols=stations*classes)

  • state_space_hashed: Hashed state indices for lookup

  • node_state_space: Dictionary of per-node state spaces

  • sn: Updated network structure with space field populated

Return type:

CtmcSsgResult containing

References

MATLAB: matlab/src/api/mc/ctmc_ssg.m

ctmc_ssg_reachability(sn, options=None)[source]

Generate reachable CTMC state space for a queueing network.

Creates only the states reachable from the initial state through valid transitions. This is more efficient than ctmc_ssg for networks with constrained reachability.

Parameters:

  • sn (Any) – NetworkStruct object (from getStruct())

  • options (Dict | None) – Solver options dict with fields: - config.hide_immediate: Hide immediate transitions (default True)

Returns:

  • state_space: Reachable state space matrix

  • state_space_aggr: Aggregated state space (per station-class)

  • state_space_hashed: Hashed state indices

  • node_state_space: Dictionary of per-node state spaces

  • sn: Updated network structure

Return type:

CtmcSsgResult containing

References

MATLAB: matlab/src/api/mc/ctmc_ssg_reachability.m

class CtmcSsgResult(state_space, state_space_aggr, state_space_hashed, node_state_space, sn)[source]

Bases: object

Result from CTMC state space generation.

state_space: numpy.ndarray
state_space_aggr: numpy.ndarray
state_space_hashed: numpy.ndarray
node_state_space: Dict[int, numpy.ndarray]
sn: Any
dtmc_solve(P)[source]

Solve for steady-state probabilities of a DTMC.

Computes the stationary distribution π by solving π(P - I) = 0 with normalization constraint Σπ = 1.

This leverages the CTMC solver by treating (P - I) as an infinitesimal generator.

Parameters:

P (numpy.ndarray) – Transition probability matrix (row stochastic)

Returns:

Steady-state probability distribution (1D array)

Return type:

numpy.ndarray

dtmc_solve_reducible(P, pin=None)[source]

Solve reducible DTMCs with transient states.

Handles DTMCs with multiple recurrent classes and transient states by: 1. Decomposing into strongly connected components (SCCs) 2. Identifying recurrent vs transient SCCs 3. Computing limiting distribution considering absorption from transient states

For a reducible DTMC with a single transient SCC, this computes the limiting distribution when starting from the transient states (e.g., class switching networks where jobs start in a transient class).

Parameters:

  • P (numpy.ndarray) – Transition probability matrix (possibly reducible)

  • pin (numpy.ndarray) – Initial probability vector (optional)

Returns:

Steady-state probability vector

Return type:

numpy.ndarray

dtmc_makestochastic(A)[source]

Convert matrix to row-stochastic transition matrix.

Normalizes each row to sum to 1. Rows with zero sum are replaced with uniform distribution.

Parameters:

A (numpy.ndarray) – Input matrix to normalize

Returns:

Row-stochastic matrix

Return type:

numpy.ndarray

dtmc_isfeasible(P, tolerance=1e-10)[source]

Check if matrix is a valid DTMC transition matrix.

Validates: - All elements are non-negative - All row sums equal 1

Parameters:

  • P (numpy.ndarray) – Candidate transition matrix

  • tolerance (float) – Numerical tolerance (default: 1e-10)

Returns:

True if matrix is valid DTMC transition matrix

Return type:

bool

dtmc_simulate(P, initial_state, num_steps, seed=None)[source]

Simulate DTMC sample path.

Generates a realization of the discrete-time Markov chain for a specified number of steps.

Parameters:

  • P (numpy.ndarray) – Transition probability matrix

  • initial_state (int) – Starting state index

  • num_steps (int) – Number of simulation steps

Returns:

Array of visited states (length num_steps + 1)

Return type:

numpy.ndarray

dtmc_rand(n, density=0.5)[source]

Generate random DTMC transition matrix.

Parameters:

  • n (int) – Number of states

  • density (float) – Sparsity density (0 to 1, default 0.5)

Returns:

Random transition probability matrix

Return type:

numpy.ndarray

dtmc_timereverse(P, pi=None)[source]

Compute time-reversed DTMC transition matrix.

The time-reversed chain has transition probabilities: P*_{ij} = π_j * P_{ji} / π_i

Parameters:

  • P (numpy.ndarray) – Original transition matrix

  • pi (numpy.ndarray | None) – Steady-state distribution (optional, computed if None)

Returns:

Time-reversed transition matrix

Return type:

numpy.ndarray

dtmc_stochcomp(P, keep_states, eliminate_states=None)[source]

Compute stochastic complement of DTMC.

Reduces the DTMC by eliminating specified states while preserving the steady-state distribution restricted to the kept states.

Parameters:

  • P (numpy.ndarray) – Transition probability matrix

  • keep_states (numpy.ndarray) – States to retain in reduced model

  • eliminate_states (numpy.ndarray | None) – States to eliminate (optional, inferred if None)

Returns:

Reduced transition matrix (stochastic complement)

Return type:

numpy.ndarray

dtmc_transient(P, initial_dist, steps)[source]

Compute transient probabilities of a DTMC.

Calculates π(n) = π(0) * P^n for each step from 0 to steps.

Parameters:

  • P (numpy.ndarray) – Transition probability matrix

  • initial_dist (numpy.ndarray) – Initial probability distribution π(0)

  • steps (int) – Number of time steps

Returns:

Array of shape (steps+1, n) with transient probabilities

Return type:

numpy.ndarray

dtmc_hitting_time(P, target_states)[source]

Compute mean hitting times to target states.

Calculates the expected number of steps to reach any target state from each starting state.

Parameters:
Returns:

Array of mean hitting times from each state

Return type:

numpy.ndarray

class CourtoisResult(p, Qperm, Qdec, eps, epsMAX, P, B, q)[source]

Bases: object

Result of Courtois decomposition.

p: numpy.ndarray
Qperm: numpy.ndarray
Qdec: numpy.ndarray
eps: float
epsMAX: float
P: numpy.ndarray
B: numpy.ndarray
q: float
class KMSResult(p, p_1, Qperm, eps, epsMAX, pcourt)[source]

Bases: object

Result of KMS aggregation-disaggregation.

p: numpy.ndarray
p_1: numpy.ndarray
Qperm: numpy.ndarray
eps: float
epsMAX: float
pcourt: numpy.ndarray
class TakahashiResult(p, p_1, pcourt, Qperm, eps, epsMAX)[source]

Bases: object

Result of Takahashi aggregation-disaggregation.

p: numpy.ndarray
p_1: numpy.ndarray
pcourt: numpy.ndarray
Qperm: numpy.ndarray
eps: float
epsMAX: float
ctmc_courtois(Q, MS, q=None)[source]

Courtois decomposition for near-completely decomposable CTMCs.

Decomposes a large CTMC into macrostates and computes approximate steady-state probabilities using hierarchical aggregation.

Parameters:

  • Q (numpy.ndarray) – Infinitesimal generator matrix

  • MS (List[List[int]]) – List where MS[i] is the list of state indices in macrostate i

  • q (float | None) – Randomization coefficient (optional)

Returns:

CourtoisResult with approximate solution and diagnostics

Return type:

CourtoisResult

References

Original MATLAB: matlab/src/api/mc/ctmc_courtois.m Courtois, “Decomposability: Queueing and Computer System Applications”, 1977

ctmc_kms(Q, MS, numSteps=10)[source]

Koury-McAllister-Stewart aggregation-disaggregation method.

Iteratively refines the Courtois decomposition solution using aggregation and disaggregation steps.

Parameters:

  • Q (numpy.ndarray) – Infinitesimal generator matrix

  • MS (List[List[int]]) – List where MS[i] is the list of state indices in macrostate i

  • numSteps (int) – Number of iterative steps (default: 10)

Returns:

KMSResult with refined solution

Return type:

KMSResult

References

Original MATLAB: matlab/src/api/mc/ctmc_kms.m Koury, McAllister, Stewart, “Iterative Methods for Computing Stationary Distributions of Nearly Completely Decomposable Markov Chains”, 1984

ctmc_takahashi(Q, MS, numSteps=10)[source]

Takahashi’s aggregation-disaggregation method.

Iteratively refines the Courtois decomposition solution using a different aggregation-disaggregation scheme.

Parameters:

  • Q (numpy.ndarray) – Infinitesimal generator matrix

  • MS (List[List[int]]) – List where MS[i] is the list of state indices in macrostate i

  • numSteps (int) – Number of iterative steps (default: 10)

Returns:

TakahashiResult with refined solution

Return type:

TakahashiResult

References

Original MATLAB: matlab/src/api/mc/ctmc_takahashi.m Takahashi, “A Lumping Method for Numerical Calculations of Stationary Distributions of Markov Chains”, 1975

ctmc_multi(Q, MS, MSS)[source]

Multigrid aggregation-disaggregation method.

Two-level hierarchical decomposition using nested macrostates.

Parameters:

  • Q (numpy.ndarray) – Infinitesimal generator matrix

  • MS (List[List[int]]) – List where MS[i] is the list of state indices in macrostate i

  • MSS (List[List[int]]) – List where MSS[i] is the list of macrostate indices in macro-macrostate i

Returns:

TakahashiResult with multigrid solution

Return type:

TakahashiResult

References

Original MATLAB: matlab/src/api/mc/ctmc_multi.m

Queueing Systems (line_solver.api.qsys)

Native Python implementations for queueing system analysis.

This module provides pure Python/NumPy implementations for analyzing single queueing systems, including basic queues (M/M/1, M/M/k, M/G/1), G/G/1 approximations, MAP-based queues, and scheduling disciplines.

Key algorithms:

Basic queues: qsys_mm1, qsys_mmk, qsys_mg1, qsys_gm1, qsys_mminf, qsys_mginf G/G/1 approximations: Allen-Cunneen, Kingman, Marchal, Whitt, Heyman, etc. G/G/k approximations: qsys_gigk_approx MAP/D queues: qsys_mapdc, qsys_mapd1 MAP/PH queues: qsys_phph1, qsys_mapph1, qsys_mapm1, qsys_mapmc, qsys_mapmap1 Scheduling: qsys_mg1_prio, qsys_mg1_srpt, qsys_mg1_fb, etc. Loss systems: qsys_mm1k_loss, qsys_mg1k_loss, qsys_mxm1

qsys_mapdc(D0, D1, s, c, max_num_comp=1000, num_steps=1, verbose=0)[source]

Analyze MAP/D/c queue (MAP arrivals, deterministic service, c servers).

Uses Non-Skip-Free (NSF) Markov chain analysis embedding at deterministic service intervals. Multiple arrivals can occur per interval.

Parameters:

  • D0 (numpy.ndarray) – MAP hidden transition matrix (n x n).

  • D1 (numpy.ndarray) – MAP arrival transition matrix (n x n).

  • s (float) – Deterministic service time (positive scalar).

  • c (int) – Number of servers.

  • max_num_comp (int) – Maximum number of queue length components (default 1000).

  • num_steps (int) – Number of waiting time distribution points per interval (default 1).

  • verbose (int) – Verbosity level (default 0).

Returns:

Performance metrics including:

  • mean_queue_length: Mean number of customers in system

  • mean_waiting_time: Mean waiting time in queue

  • mean_sojourn_time: Mean sojourn time (waiting + service)

  • utilization: Server utilization (per server)

  • queue_length_dist: Queue length distribution P(Q=n)

  • waiting_time_dist: Waiting time CDF at discrete points

  • analyzer: Analyzer identifier

Return type:

dict

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

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

Parameters:

  • D0 (numpy.ndarray) – MAP hidden transition matrix.

  • D1 (numpy.ndarray) – MAP arrival transition matrix.

  • s (float) – Deterministic service time.

  • max_num_comp (int) – Maximum number of queue length components.

  • num_steps (int) – Number of waiting time points per interval.

Returns:

Performance metrics (see qsys_mapdc).

Return type:

dict

qsys_mm1(lambda_val, mu)[source]

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

Parameters:

  • lambda_val (float) – Arrival rate (lambda)

  • mu (float) – Service rate

Returns:

Performance measures including:

  • L: Mean number in system

  • Lq: Mean number in queue

  • W: Mean response time (time in system)

  • Wq: Mean waiting time (time in queue)

  • rho: Utilization (lambda/mu)

Return type:

dict

Example

>>> result = qsys_mm1(0.5, 1.0)
>>> print(f"Utilization: {result['rho']:.2f}")
Utilization: 0.50
qsys_mmk(lambda_val, mu, k)[source]

Analyze M/M/k queue (Poisson arrivals, k exponential servers).

Parameters:

  • lambda_val (float) – Arrival rate (lambda)

  • mu (float) – Service rate per server

  • k (int) – Number of parallel servers

Returns:

Performance measures including:

  • L: Mean number in system

  • Lq: Mean number in queue

  • W: Mean response time

  • Wq: Mean waiting time

  • rho: Utilization per server (lambda/(k*mu))

  • P0: Probability of empty system

Return type:

dict

Example

>>> result = qsys_mmk(2.0, 1.0, 3)
>>> print(f"Utilization: {result['rho']:.2f}")
Utilization: 0.67
qsys_mg1(lambda_val, mu, cs)[source]

Analyze M/G/1 queue using Pollaczek-Khinchine formula.

Parameters:

  • lambda_val (float) – Arrival rate (lambda)

  • mu (float) – Service rate (mean service time = 1/mu)

  • cs (float) – Coefficient of variation of service time (std/mean)

Returns:

Performance measures including:

  • L: Mean number in system

  • Lq: Mean number in queue

  • W: Mean response time

  • Wq: Mean waiting time

  • rho: Utilization (lambda/mu)

Return type:

dict

Example

>>> result = qsys_mg1(0.5, 1.0, 1.0)  # cs=1 is exponential (M/M/1)
qsys_gm1(lambda_val, mu, ca)[source]

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

Parameters:

  • lambda_val (float) – Arrival rate (lambda)

  • mu (float) – Service rate

  • ca (float) – Coefficient of variation of inter-arrival time

Returns:

Performance measures including:

  • L: Mean number in system

  • Lq: Mean number in queue

  • W: Mean response time

  • Wq: Mean waiting time

  • rho: Utilization

Return type:

dict

Note

Uses approximation based on the coefficient of variation.

qsys_mminf(lambda_val, mu)[source]

Analyze M/M/inf queue (infinite servers / delay station).

Parameters:

  • lambda_val (float) – Arrival rate (lambda)

  • mu (float) – Service rate

Returns:

Performance measures including:

  • L: Mean number in system (= lambda/mu)

  • Lq: Mean number in queue (= 0)

  • W: Mean time in system (= 1/mu)

  • Wq: Mean waiting time (= 0)

  • P0: Probability of empty system

Return type:

dict

qsys_mginf(lambda_val, mu, k=None)[source]

Analyze M/G/inf queue (infinite servers, general service).

Performance is independent of service time distribution shape. Number of customers follows Poisson distribution.

Parameters:

  • lambda_val (float) – Arrival rate (lambda)

  • mu (float) – Service rate (mean service time = 1/mu)

  • k (int | None) – Optional state for probability computation

Returns:

Performance measures including:

  • L: Mean number in system

  • Lq: Mean number in queue (= 0)

  • W: Mean time in system (= 1/mu)

  • Wq: Mean waiting time (= 0)

  • P0: Probability of empty system

  • Pk: Probability of k customers (if k provided)

Return type:

dict

qsys_gig1_approx_allencunneen(lambda_val, mu, ca, cs)[source]

Allen-Cunneen approximation for G/G/1 queue.

Matches MATLAB qsys_gig1_approx_allencunneen.m exactly.

Parameters:

  • lambda_val (float) – Arrival rate

  • mu (float) – Service rate

  • ca (float) – Coefficient of variation of inter-arrival time

  • cs (float) – Coefficient of variation of service time

Returns:

W: Mean response time (time in system) rhohat: Effective utilization (so that M/M/1 formulas still hold)

Return type:

Tuple of (W, rhohat) where

qsys_gig1_approx_kingman(lambda_val, mu, ca, cs)[source]

Kingman’s upper bound approximation for G/G/1 queue.

Note: In MATLAB, ‘gig1’ and ‘gig1.kingman’ map to qsys_gig1_ubnd_kingman. This function provides the same formula.

Parameters:

  • lambda_val (float) – Arrival rate

  • mu (float) – Service rate

  • ca (float) – Coefficient of variation of inter-arrival time

  • cs (float) – Coefficient of variation of service time

Returns:

W: Mean response time (time in system) rhohat: Effective utilization (so that M/M/1 formulas still hold)

Return type:

Tuple of (W, rhohat) where

qsys_gig1_approx_marchal(lambda_val, mu, ca, cs)[source]

Marchal’s approximation for G/G/1 queue.

Matches MATLAB qsys_gig1_approx_marchal.m exactly. Note: MATLAB formula uses ca (not ca^2) in the numerator factor.

Parameters:

  • lambda_val (float) – Arrival rate

  • mu (float) – Service rate

  • ca (float) – Coefficient of variation of inter-arrival time

  • cs (float) – Coefficient of variation of service time

Returns:

W: Mean response time (time in system) rhohat: Effective utilization (so that M/M/1 formulas still hold)

Return type:

Tuple of (W, rhohat) where

qsys_gig1_approx_whitt(lambda_val, mu, ca, cs)[source]

Whitt’s approximation for G/G/1 queue.

Uses QNA (Queueing Network Analyzer) approximation. Note: No direct MATLAB counterpart (qsys_gig1_approx_whitt.m does not exist).

Parameters:

  • lambda_val (float) – Arrival rate

  • mu (float) – Service rate

  • ca (float) – Coefficient of variation of inter-arrival time

  • cs (float) – Coefficient of variation of service time

Returns:

W: Mean response time (time in system) rhohat: Effective utilization (so that M/M/1 formulas still hold)

Return type:

Tuple of (W, rhohat) where

qsys_gig1_approx_heyman(lambda_val, mu, ca, cs)[source]

Heyman’s approximation for G/G/1 queue.

Matches MATLAB qsys_gig1_approx_heyman.m exactly.

Parameters:

  • lambda_val (float) – Arrival rate

  • mu (float) – Service rate

  • ca (float) – Coefficient of variation of inter-arrival time

  • cs (float) – Coefficient of variation of service time

Returns:

W: Mean response time (time in system) rhohat: Effective utilization (so that M/M/1 formulas still hold)

Return type:

Tuple of (W, rhohat) where

qsys_gig1_approx_kobayashi(lambda_val, mu, ca, cs)[source]

Kobayashi’s approximation for G/G/1 queue.

Matches MATLAB qsys_gig1_approx_kobayashi.m exactly.

Parameters:

  • lambda_val (float) – Arrival rate

  • mu (float) – Service rate

  • ca (float) – Coefficient of variation of inter-arrival time

  • cs (float) – Coefficient of variation of service time

Returns:

W: Mean response time (time in system) rhohat: Effective utilization

Return type:

Tuple of (W, rhohat) where

qsys_gig1_approx_klb(lambda_val, mu, ca, cs)[source]

Kraemer-Langenbach-Belz (KLB) approximation for G/G/1 queue.

Matches MATLAB qsys_gig1_approx_klb.m exactly.

Parameters:

  • lambda_val (float) – Arrival rate

  • mu (float) – Service rate

  • ca (float) – Coefficient of variation of inter-arrival time

  • cs (float) – Coefficient of variation of service time

Returns:

W: Mean response time (time in system) rhohat: Effective utilization (so that M/M/1 formulas still hold)

Return type:

Tuple of (W, rhohat) where

qsys_gig1_approx_gelenbe(lambda_val, mu, ca, cs)[source]

Gelenbe’s approximation for G/G/1 queue.

Matches MATLAB qsys_gig1_approx_gelenbe.m exactly. Note: MATLAB returns only W (no rhohat).

Parameters:

  • lambda_val (float) – Arrival rate

  • mu (float) – Service rate

  • ca (float) – Coefficient of variation of inter-arrival time

  • cs (float) – Coefficient of variation of service time

Returns:

Mean response time W

Return type:

float

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

Kimura’s approximation for G/G/1 queue.

Matches MATLAB qsys_gig1_approx_kimura.m exactly. Note: MATLAB uses sigma (=rho) as first parameter, returns only W.

Parameters:

  • lambda_val (float) – Arrival rate (used as sigma = lambda/mu = rho)

  • mu (float) – Service rate

  • ca (float) – Coefficient of variation of inter-arrival time

  • cs (float) – Coefficient of variation of service time

Returns:

Mean response time W

Return type:

float

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

Myskja’s approximation for G/G/1 queue.

Uses third moments of arrival and service distributions for improved accuracy.

Reference: Myskja, A. (1991). “An Experimental Study of a H₂/H₂/1 Queue”. Stochastic Models, 7(4), 571-595.

Parameters:

  • lambda_val (float) – Arrival rate

  • mu (float) – Service rate

  • ca (float) – Coefficient of variation of inter-arrival time

  • cs (float) – Coefficient of variation of service time

  • q0 (float) – Lowest relative third moment for given mean and SCV

  • qa (float) – Third relative moment E[X^3]/6/E[X]^3 of inter-arrival time

Returns:

Waiting time W

Return type:

float

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

Modified Myskja (Myskja2) approximation for G/G/1 queue.

An improved version of Myskja’s approximation with better accuracy for a wider range of parameter combinations.

Reference: Myskja, A. (1991). “An Experimental Study of a H₂/H₂/1 Queue”. Stochastic Models, 7(4), 571-595.

Parameters:

  • lambda_val (float) – Arrival rate

  • mu (float) – Service rate

  • ca (float) – Coefficient of variation of inter-arrival time

  • cs (float) – Coefficient of variation of service time

  • q0 (float) – Lowest relative third moment for given mean and SCV

  • qa (float) – Third relative moment E[X^3]/6/E[X]^3 of inter-arrival time

Returns:

Waiting time W

Return type:

float

qsys_gg1(lambda_val, mu, ca2, cs2)[source]

G/G/1 queue analysis using exact methods for special cases and Allen-Cunneen approximation for the general case.

Parameters:

  • lambda_val (float) – Arrival rate

  • mu (float) – Service rate

  • ca2 (float) – Squared coefficient of variation of inter-arrival time

  • cs2 (float) – Squared coefficient of variation of service time

Returns:

W: Mean response time (time in system) rhohat: Effective utilization

Return type:

Tuple of (W, rhohat) where

References

Original JAR: jar/src/main/kotlin/jline/api/qsys/Qsys_gg1.kt

qsys_gig1_lbnd(lambda_val, mu, ca, cs)[source]

Fundamental theoretical lower bounds for G/G/1 queues.

These are the minimum possible values that performance measures cannot fall below for any realization of the arrival and service processes.

Parameters:

  • lambda_val (float) – Arrival rate

  • mu (float) – Service rate

  • ca (float) – Coefficient of variation of inter-arrival time

  • cs (float) – Coefficient of variation of service time

Returns:

W: Lower bound on mean response time (= 1/mu) rhohat: Effective utilization

Return type:

Tuple of (W, rhohat) where

References

Original JAR: jar/src/main/kotlin/jline/api/qsys/Qsys_gig1_lbnd.kt

qsys_gigk_approx(lambda_val, mu, ca, cs, k)[source]

Approximation for G/G/k queue.

Matches MATLAB qsys_gigk_approx.m formula using alpha-factor correction.

Parameters:

  • lambda_val (float) – Arrival rate

  • mu (float) – Service rate per server

  • ca (float) – Coefficient of variation of inter-arrival time

  • cs (float) – Coefficient of variation of service time

  • k (int) – Number of servers

Returns:

W: Approximate mean response time rhohat: Effective utilization

Return type:

Tuple of (W, rhohat) where

qsys_gigk_approx_cosmetatos(lambda_val, mu, ca, cs, k)[source]

G/G/k queue approximation using the Cosmetatos method.

Adjusts M/M/k results based on the variability of arrival and service processes.

Parameters:

  • lambda_val (float) – Arrival rate

  • mu (float) – Service rate per server

  • ca (float) – Coefficient of variation of inter-arrival time

  • cs (float) – Coefficient of variation of service time

  • k (int) – Number of servers

Returns:

W: Approximate mean response time rhohat: Effective utilization

Return type:

Tuple of (W, rhohat) where

References

Original JAR: jar/src/main/kotlin/jline/api/qsys/Qsys_gigk_approx_cosmetatos.kt

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

G/G/k queue approximation using Whitt’s QNA method.

Uses the Queue Network Analyzer (QNA) methodology to provide accurate estimates for multi-server queues with general distributions.

Parameters:

  • lambda_val (float) – Arrival rate

  • mu (float) – Service rate per server

  • ca (float) – Coefficient of variation of inter-arrival time

  • cs (float) – Coefficient of variation of service time

  • k (int) – Number of servers

Returns:

W: Approximate mean response time rhohat: Effective utilization

Return type:

Tuple of (W, rhohat) where

References

Original JAR: jar/src/main/kotlin/jline/api/qsys/Qsys_gigk_approx_whitt.kt

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

Kingman’s upper bound for G/G/1 queue waiting time.

This is an upper bound approximation that provides conservative estimates for the mean waiting time.

Parameters:

  • lambda_val (float) – Arrival rate

  • mu (float) – Service rate

  • ca (float) – Coefficient of variation of inter-arrival time

  • cs (float) – Coefficient of variation of service time

Returns:

W: Upper bound on mean response time rhohat: Effective utilization (so M/M/1 formulas still hold)

Return type:

Tuple of (W, rhohat) where

References

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

qsys_gigk_approx_kingman(lambda_val, mu, ca, cs, k)[source]

Kingman’s approximation for G/G/k queue waiting time.

Extends Kingman’s approximation to multi-server queues using M/M/k waiting time as a base.

Parameters:

  • lambda_val (float) – Arrival rate

  • mu (float) – Service rate per server

  • k (int) – Number of servers

  • ca (float) – Coefficient of variation of inter-arrival time

  • cs (float) – Coefficient of variation of service time

Returns:

W: Approximate mean response time rhohat: Effective utilization

Return type:

Tuple of (W, rhohat) where

References

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

qsys_mg1_prio(lambda_vec, mu_vec, cs_vec)[source]

Analyze M/G/1 queue with non-preemptive (Head-of-Line) priorities.

Matches MATLAB qsys_mg1_prio.m exactly.

Parameters:

  • lambda_vec (numpy.ndarray) – Vector of arrival rates per priority class (class 1 = highest)

  • mu_vec (numpy.ndarray) – Vector of service rates per priority class

  • cs_vec (numpy.ndarray) – Vector of coefficients of variation per priority class

Returns:

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

Return type:

Tuple of (W, rho)

qsys_mg1_srpt(lambda_vec, mu_vec, cs_vec)[source]

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

Matches MATLAB qsys_mg1_srpt.m exactly: - Exponential case: preemptive priority with cumulative residual E_R_k - General case: Schrage-Miller formula with numerical integration

Parameters:

  • lambda_vec (numpy.ndarray) – Vector of arrival rates per class

  • mu_vec (numpy.ndarray) – Vector of service rates per class

  • cs_vec (numpy.ndarray) – Vector of coefficients of variation per class

Returns:

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

Return type:

Tuple of (W, rho)

qsys_mg1_fb(lambda_vec, mu_vec, cs_vec)[source]

Analyze M/G/1 queue with Foreground-Background (FB/LAS) scheduling.

Matches MATLAB qsys_mg1_fb.m exactly: - Exponential case: numerical integration of E[T(x)] * f_k(x) - General case: class-based approximation

Parameters:

  • lambda_vec (numpy.ndarray) – Vector of arrival rates per class

  • mu_vec (numpy.ndarray) – Vector of service rates per class

  • cs_vec (numpy.ndarray) – Vector of coefficients of variation per class

Returns:

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

Return type:

Tuple of (W, rho)

qsys_mg1_lrpt(lambda_vec, mu_vec, cs_vec)[source]

Analyze M/G/1 queue with Longest Remaining Processing Time (LRPT).

Matches MATLAB qsys_mg1_lrpt.m exactly: - Exponential case: numerical integration of E[T(x)] * f_k(x) - General case: preemptive priority with descending service time ordering

Parameters:

  • lambda_vec (numpy.ndarray) – Vector of arrival rates per class

  • mu_vec (numpy.ndarray) – Vector of service rates per class

  • cs_vec (numpy.ndarray) – Vector of coefficients of variation per class

Returns:

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

Return type:

Tuple of (W, rho)

qsys_mg1_psjf(lambda_vec, mu_vec, cs_vec)[source]

Analyze M/G/1 queue with Preemptive Shortest Job First (PSJF).

Matches MATLAB qsys_mg1_psjf.m exactly: - Exponential case: numerical integration of E[T(x)] * f_k(x) - General case: class-based truncated moment formula

Parameters:

  • lambda_vec (numpy.ndarray) – Vector of arrival rates per class

  • mu_vec (numpy.ndarray) – Vector of service rates per class

  • cs_vec (numpy.ndarray) – Vector of coefficients of variation per class

Returns:

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

Return type:

Tuple of (W, rho)

qsys_mg1_setf(lambda_vec, mu_vec, cs_vec)[source]

Analyze M/G/1 queue with Shortest Expected Time First (SETF).

Matches MATLAB qsys_mg1_setf.m exactly: SETF = FB/LAS + residual service time penalty (non-preemptive).

Parameters:

  • lambda_vec (numpy.ndarray) – Vector of arrival rates per class

  • mu_vec (numpy.ndarray) – Vector of service rates per class

  • cs_vec (numpy.ndarray) – Vector of coefficients of variation per class

Returns:

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

Return type:

Tuple of (W, rho)

qsys_mm1k_loss(lambda_val, mu, K)[source]

Compute loss probability for M/M/1/K queue.

Uses the closed-form formula for the M/M/1/K loss system where customers are rejected when the buffer (capacity K) is full.

Parameters:

  • lambda_val (float) – Arrival rate

  • mu (float) – Service rate

  • K (int) – Buffer capacity (including customer in service)

Returns:

lossprob: Probability that an arriving customer is rejected rho: Offered load (lambda/mu)

Return type:

Tuple of (lossprob, rho)

References

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

qsys_mg1k_loss(lambda_val, service_pdf, K, max_t=100.0)[source]

Compute loss probability for M/G/1/K queue.

Uses the Niu-Cooper transform-free analysis method for computing the stationary distribution and loss probability of M/G/1/K queues.

Parameters:

  • lambda_val (float) – Arrival rate

  • service_pdf (Callable[[float], float]) – Probability density function of service time f(t)

  • K (int) – Buffer capacity

  • max_t (float) – Maximum integration time (default: 100)

Returns:

sigma: Normalizing constant term rho: Offered load lossprob: Probability of loss

Return type:

Tuple of (sigma, rho, lossprob)

References

Original MATLAB: matlab/src/api/qsys/qsys_mg1k_loss.m Niu-Cooper, “Transform-Free Analysis of M/G/1/K”, 1993

qsys_mg1k_loss_mgs(lambda_val, mu, mu_scv, K)[source]

Compute loss probability for M/G/1/K using MacGregor Smith approximation.

Matches MATLAB qsys_mg1k_loss_mgs.m exactly.

Parameters:

  • lambda_val (float) – Arrival rate

  • mu (float) – Service rate

  • mu_scv (float) – Squared coefficient of variation of service time

  • K (int) – Buffer capacity

Returns:

lossprob: Probability of loss rho: Offered load

Return type:

Tuple of (lossprob, rho)

References

Original MATLAB: matlab/src/api/qsys/qsys_mg1k_loss_mgs.m J. MacGregor Smith, “Optimal Design and Performance Modelling of M/G/1/K Queueing Systems”

qsys_mxm1(lambda_batch, mu, E_X_or_batch_sizes, E_X2_or_pmf, mode=None)[source]

Analyze MX/M/1 queue with batch arrivals.

Matches MATLAB qsys_mxm1.m exactly.

Three input formats:

  1. Moment-based: qsys_mxm1(lambda_batch, mu, E_X, E_X2)

  2. PMF-based: qsys_mxm1(lambda_batch, mu, batch_sizes, pmf)

  3. Variance: qsys_mxm1(lambda_batch, mu, E_X, Var_X, ‘variance’)

Parameters:

  • lambda_batch (float) – Batch arrival rate

  • mu (float) – Service rate

  • E_X_or_batch_sizes – Mean batch size (scalar) or array of batch sizes

  • E_X2_or_pmf – Second moment of batch size, PMF, or variance

  • mode (str | None) – Optional ‘variance’ flag for variance-based input

Returns:

W: Mean time in system Wq: Mean waiting time in queue U: Server utilization Q: Mean queue length (including service)

Return type:

Tuple of (W, Wq, U, Q)

References

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

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

Bases: object

Result structure for queue analysis.

analyzer: str = 'native'
queueLengthDist: numpy.ndarray | None = None
queueLengthMoments: numpy.ndarray | None = None
sojournTimeMoments: numpy.ndarray | None = None
meanQueueLength: float
meanWaitingTime: float
meanSojournTime: float
utilization: float
ph_to_map(alpha, T)[source]

Convert a PH distribution to its equivalent MAP representation.

For a PH renewal process, the MAP has:

D0 = T (transitions within the PH, no arrival) D1 = t * alpha where t = -T*e (exit rates times restart distribution)

Parameters:
Returns:

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

Return type:

Tuple of (D0, D1)

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

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

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

Parameters:

  • alpha (numpy.ndarray) – Arrival PH initial probability vector (1 x n)

  • T (numpy.ndarray) – Arrival PH generator matrix (n x n)

  • beta (numpy.ndarray) – Service PH initial probability vector (1 x m)

  • S (numpy.ndarray) – Service PH generator matrix (m x m)

  • numQLMoms (int) – Number of queue length moments to compute (default: 3)

  • numQLProbs (int) – Number of queue length probabilities (default: 100)

  • numSTMoms (int) – Number of sojourn time moments (default: 3)

Returns:

QueueResult with queue performance metrics

Return type:

QueueResult

References

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

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

Analyze a MAP/PH/1 queue.

Parameters:

  • D0 (numpy.ndarray) – MAP hidden transition matrix

  • D1 (numpy.ndarray) – MAP observable transition matrix

  • beta (numpy.ndarray) – Service PH initial probability vector

  • S (numpy.ndarray) – Service PH generator matrix

  • numQLMoms (int) – Number of queue length moments

  • numQLProbs (int) – Number of queue length probabilities

  • numSTMoms (int) – Number of sojourn time moments

Returns:

QueueResult with queue performance metrics

Return type:

QueueResult

References

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

qsys_mapm1(D0, D1, mu)[source]

Analyze a MAP/M/1 queue.

Parameters:
Returns:

QueueResult with queue performance metrics

Return type:

QueueResult

References

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

qsys_mapmc(D0, D1, mu, c)[source]

Analyze a MAP/M/c queue.

Parameters:

  • D0 (numpy.ndarray) – MAP hidden transition matrix

  • D1 (numpy.ndarray) – MAP observable transition matrix

  • mu (float) – Service rate per server

  • c (int) – Number of servers

Returns:

QueueResult with queue performance metrics

Return type:

QueueResult

References

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

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

Analyze a MAP/MAP/1 queue.

Both arrival and service processes are Markovian Arrival Processes.

Parameters:

  • D0_arr (numpy.ndarray) – Arrival MAP hidden transition matrix

  • D1_arr (numpy.ndarray) – Arrival MAP observable transition matrix

  • D0_srv (numpy.ndarray) – Service MAP hidden transition matrix

  • D1_srv (numpy.ndarray) – Service MAP observable transition matrix

Returns:

QueueResult with queue performance metrics

Return type:

QueueResult

References

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

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

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

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

Parameters:

  • D0 (numpy.ndarray) – MAP hidden transition matrix (n x n)

  • D1 (numpy.ndarray) – MAP arrival transition matrix (n x n)

  • service_moments (numpy.ndarray) – First k raw moments of service time [E[S], E[S^2], …] (k = 2 or 3 for best accuracy)

  • num_ql_moms (int) – Number of queue length moments to compute (default: 3)

  • num_ql_probs (int) – Number of queue length probabilities (default: 100)

  • num_st_moms (int) – Number of sojourn time moments to compute (default: 3)

Returns:

QueueResult with queue performance metrics

Return type:

QueueResult

References

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

Note

Uses the MMAPPH1FCFS solver from BuTools after fitting the general service distribution to a PH distribution.

class QueueType(*values)[source]

Bases: Enum

Queueing system topology types.

STANDARD = 'standard'
RETRIAL = 'retrial'
RENEGING = 'reneging'
RETRIAL_RENEGING = 'retrial_reneging'
class BmapMatrix(D0, D_batch)[source]

Bases: object

Batch Markovian Arrival Process matrix representation.

D₀ = drift matrix (no arrivals) D₁, D₂, …, Dₖ = batch arrival matrices for batch sizes 1, 2, …, K

Properties:

order: Dimension of the BMAP num_batches: Maximum batch size arrival_rate: Overall arrival rate

__post_init__()[source]

Validate BMAP structure.

property arrival_rate: float

Compute overall arrival rate of the BMAP.

lambda = theta * (-D0) * e where theta is stationary distribution of BMAP Markov chain and e is unit vector.

property fundamental_arrival_rate: float

Arrival rate from fundamental matrix.

D0: numpy.ndarray
D_batch: List[numpy.ndarray]
class PhDistribution(beta, S)[source]

Bases: object

Phase-Type (PH) distribution representation.

Beta = initial probability vector (shape: m,) S = transient generator matrix (shape: m x m)

where m is the number of phases.

Properties:

mean: Mean of PH distribution (1/mu in queue notation) scv: Squared coefficient of variation num_phases: Number of phases

__post_init__()[source]

Validate PH distribution structure.

property mean: float

Compute mean of PH distribution.

E[X] = -beta * S^{-1} * e

property mean_squared: float

Compute second moment of PH distribution.

E[X²] = 2 * beta * S^{-2} * e

property scv: float

Squared coefficient of variation of PH distribution.

SCV = (E[X²] / E[X]²) - 1

beta: numpy.ndarray
S: numpy.ndarray
class QbdStatespace(A0, A1, A2, B0, max_level)[source]

Bases: object

Quasi-Birth-Death Markov chain state space representation.

For a QBD process, states are of the form (n, i) where: - n = level (number of retrying customers in orbit) - i = phase (service phase or phase-type stage)

Generator matrix has block structure: Q = | B_0 A_0 0 0 … |

A_2 A_1 A_0 0 … |
0 A_2 A_1 A_0 … |
… |
Properties:

max_level: Maximum retrial orbit size (truncation level) phase_dim: Number of phases at each level total_states: Total state space dimension

__post_init__()[source]

Validate QBD state space.

property binomial_state_dimension: int

Compute state space dimension using binomial coefficient formula.

For a retrial queue with N total customers and m servers, dimension = C(N + m - 1, m - 1)

This is a rough upper bound for QBD truncation.

get_level_phase(idx)[source]

Map linear state index to (level, phase).

Parameters:

idx (int) – Linear state index

Returns:

Tuple (level, phase)

Return type:

Tuple[int, int]

get_state_index(level, phase)[source]

Map (level, phase) to linear state index.

Parameters:

  • level (int) – Retrial orbit size (0 ≤ level ≤ max_level)

  • phase (int) – Phase within level (0 ≤ phase < phase_dim)

Returns:

Linear state index

Return type:

int

A0: numpy.ndarray
A1: numpy.ndarray
A2: numpy.ndarray
B0: numpy.ndarray
max_level: int
class RetrialQueueResult(queue_type, L_orbit, N_server, utilization, throughput, P_idle, P_empty_orbit, stationary_dist=None, truncation_level=0, converged=False, iterations=0, error=numpy.inf)[source]

Bases: object

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

queue_type

Type of queue (retrial, reneging, etc.)

Type:

line_solver.api.qsys.retrial.QueueType

L_orbit

Expected number of customers in orbit

Type:

float

N_server

Expected number of customers being served

Type:

float

utilization

Server utilization

Type:

float

throughput

System throughput

Type:

float

P_idle

Probability that server is idle

Type:

float

P_empty_orbit

Probability that orbit is empty

Type:

float

stationary_dist

Stationary distribution vector

Type:

numpy.ndarray | None

truncation_level

QBD truncation level used

Type:

int

converged

Whether numerical solution converged

Type:

bool

iterations

Number of iterations to convergence

Type:

int

error

Final error estimate

Type:

float

converged: bool = False
iterations: int = 0
stationary_dist: numpy.ndarray | None = None
truncation_level: int = 0
queue_type: QueueType
L_orbit: float
N_server: float
utilization: float
throughput: float
P_idle: float
P_empty_orbit: float
class RetrialQueueAnalyzer(sn, options=None)[source]

Bases: object

Framework for analyzing BMAP/PH/N/N retrial queues.

This class provides the foundation for future full solver implementation, including topology detection, parameter extraction, and QBD setup.

Example

analyzer = RetrialQueueAnalyzer(model) queue_type = analyzer.detect_queue_type() if queue_type == QueueType.RETRIAL:

result = analyzer.analyze()

Initialize retrial queue analyzer.

Parameters:

  • sn (Any) – NetworkStruct with queue configuration

  • options (Dict | None) – Analysis options (tolerance, max iterations, etc.)

__init__(sn, options=None)[source]

Initialize retrial queue analyzer.

Parameters:

  • sn (Any) – NetworkStruct with queue configuration

  • options (Dict | None) – Analysis options (tolerance, max iterations, etc.)

analyze()[source]

Analyze the retrial queue.

This is the main entry point for analysis: 1. Detect queue type 2. Extract arrival and service parameters 3. Build QBD state space 4. Solve for stationary distribution using matrix-analytic methods 5. Compute performance metrics

Returns:

RetrialQueueResult with performance metrics

Return type:

RetrialQueueResult

build_qbd_statespace(bmap, ph_service, retrial_params)[source]

Build QBD state space for the retrial queue.

This constructs the generator matrix blocks for the QBD process following the BMAP/PH/N/N retrial queue formulation.

Parameters:

  • bmap (BmapMatrix) – BMAP arrival process

  • ph_service (PhDistribution) – PH service distribution

  • retrial_params (Dict[str, float]) – Retrial parameters (alpha, gamma, p, R, N)

Returns:

QbdStatespace instance or None if construction fails

Return type:

QbdStatespace | None

detect_queue_type()[source]

Detect the type of queueing system from topology.

Returns:

QueueType enum indicating retrial, reneging, or combination

Return type:

QueueType

extract_bmap()[source]

Extract BMAP parameters from arrival process.

LINE stores arrival processes in MAP format: {D0, D1, D2, …} where D0 is the “hidden” generator and D1, D2, … are arrival matrices.

Returns:

BmapMatrix instance or None if arrival is not BMAP/MAP

Return type:

BmapMatrix | None

extract_ph_service()[source]

Extract Phase-Type service parameters.

LINE stores PH in MAP format: {D0, D1} D0 = T (subgenerator matrix) D1 = S0 * alpha (exit rate times initial prob)

Returns:

PhDistribution instance or None if service is not PH

Return type:

PhDistribution | None

extract_retrial_parameters()[source]

Extract retrial-specific parameters.

Returns:

  • alpha: Retrial rate (rate at which customers retry)

  • gamma: Orbit impatience rate (reneging rate)

  • p: Batch rejection probability

  • R: Threshold for admission control

  • N: Number of servers

Return type:

Dict with keys

qsys_bmapphnn_retrial(arrival_matrix, service_params, N, retrial_params=None, options=None)[source]

Analyze BMAP/PH/N/N bufferless retrial queue.

Implements the algorithm from Dudin et al., “Analysis of BMAP/PH/N-Type Queueing System with Flexible Retrials Admission Control”, Mathematics 2025, 13(9), 1434.

Parameters:

  • arrival_matrix (Dict[str, numpy.ndarray]) – Dict with ‘D0’, ‘D1’, … for BMAP matrices. D0: hidden transition matrix (V x V). D1, …, DK: arrival matrices for batch sizes 1, …, K.

  • service_params (Dict[str, numpy.ndarray]) – Dict with ‘beta’ (initial prob vector, 1xM) and ‘S’ (PH subgenerator matrix, MxM).

  • N (int) – Number of servers (also capacity, hence bufferless).

  • retrial_params (Dict[str, float] | None) – Dict with ‘alpha’ (retrial rate per customer), ‘gamma’ (impatience/abandonment rate), ‘p’ (batch rejection probability), ‘R’ (admission threshold, scalar or 1xV).

  • options (Dict | None) – Dict with optional keys: ‘MaxLevel’: max orbit level for truncation (default: auto). ‘Tolerance’: convergence tolerance (default: 1e-10). ‘Verbose’: print progress (default: False).

Returns:

RetrialQueueResult with performance metrics.

Return type:

RetrialQueueResult

References

Dudin, A., Klimenok, V., & Vishnevsky, V. (2020). Port from: matlab/src/api/qsys/qsys_bmapphnn_retrial.m

qsys_is_retrial(sn)[source]

Check if network is a valid BMAP/PH/N/N bufferless retrial queue.

Validates that the network structure matches the requirements for the BMAP/PH/N/N retrial queue solver: - Single bufferless queue (capacity == number of servers) - Retrial drop strategy configured - BMAP/MAP arrival process at source - PH/Exp service at queue - Open class model

Based on: Dudin et al., “Analysis of BMAP/PH/N-Type Queueing System with Flexible Retrials Admission Control”, Mathematics 2025, 13(9), 1434.

Parameters:

sn (Any) – NetworkStruct object

Returns:

is_retrial: True if network is valid BMAP/PH/N/N retrial topology retrial_info: RetrialInfo with parameters for the retrial solver

Return type:

Tuple of (is_retrial, retrial_info) where

References

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

class RetrialInfo(is_retrial, station_idx=None, node_idx=None, source_idx=None, class_idx=None, error_msg='', N=None, alpha=0.1, gamma=0.0, p=0.0, R=None)[source]

Bases: object

Information about a valid retrial queue topology.

N: int | None = None
R: int | None = None
alpha: float = 0.1
class_idx: int | None = None
error_msg: str = ''
gamma: float = 0.0
node_idx: int | None = None
p: float = 0.0
source_idx: int | None = None
station_idx: int | None = None
is_retrial: bool

Stochastic Networks (line_solver.api.sn)

Service Network (SN) utilities.

Native Python implementations for stochastic network structure analysis, validation, and parameter extraction.

Key classes:

NetworkStruct: Data structure summarizing network characteristics SnGetDemandsResult: Result of sn_get_demands_chain calculation

Key functions:

sn_get_demands_chain: Aggregate class-level parameters into chain-level sn_has_*: Network property predicates sn_is_*: Model type checks sn_get_*: Parameter extraction sn_validate: Network validation

class MatrixArray(input_array)[source]

Bases: ndarray

Numpy array subclass with .get() and .set() methods for API compatibility.

This class provides compatibility with the wrapper mode that uses JLine’s Matrix class which has get(i, j) and set(i, j, value) methods.

Create MatrixArray from existing array.

__array_finalize__(obj)[source]

Handle view casting and new-from-template.

__getitem__(key)[source]

Override indexing to handle 2D indexing on 1D arrays.

This provides compatibility with MATLAB-style row/column vectors where a 1D array can be indexed as (0, j) or (i, 0).

static __new__(cls, input_array)[source]

Create MatrixArray from existing array.

__setitem__(key, value)[source]

Override item setting to handle 2D indexing on 1D arrays.

get(i, j=None)[source]

Get element at index (i, j) or just i if 1D.

Parameters:

  • i – Row index (or element index for 1D)

  • j – Column index (optional, for 2D arrays)

Returns:

Element value at the specified index

set(i, j, value=None)[source]

Set element at index (i, j) or just i if 1D.

Parameters:

  • i – Row index (or element index for 1D)

  • j – Column index or value (for 1D arrays)

  • value – Value to set (optional, for 2D arrays)

class NetworkStruct(nstations=0, nstateful=0, nnodes=0, nclasses=0, nchains=0, nclosedjobs=0, njobs=<factory>, nservers=<factory>, cap=None, classcap=None, rates=<factory>, scv=<factory>, phases=None, phasessz=None, phaseshift=None, visits=<factory>, nodevisits=<factory>, inchain=<factory>, chains=<factory>, refstat=<factory>, refclass=<factory>, sched=<factory>, schedparam=None, routing=<factory>, rt=None, rtnodes=None, nodetype=<factory>, isstation=<factory>, isstateful=<factory>, isstatedep=None, nodeToStation=<factory>, nodeToStateful=<factory>, stationToNode=<factory>, stationToStateful=<factory>, statefulToNode=<factory>, statefulToStation=<factory>, state=<factory>, stateprior=<factory>, space=<factory>, lldscaling=None, cdscaling=None, ljdscaling=None, ljdcutoffs=None, classprio=None, classdeadline=None, isslc=None, issignal=None, signaltype=None, syncreply=None, connmatrix=None, nodenames=<factory>, classnames=<factory>, mu=None, phi=None, proc=None, pie=None, procid=None, lst=None, fj=None, droprule=None, nregions=0, region=None, regionrule=None, regionweight=None, regionsz=None, sync=None, gsync=None, nodeparam=None, routing_weights=None, reward=None, rtorig=None, csmask=None, nvars=None)[source]

Bases: object

Data structure summarizing network characteristics.

This class is the Python equivalent in native Python. It contains all parameters needed by solvers to analyze a queueing network.

nstations

Number of stations (queues, delays, sources, joins, places)

Type:

int

nstateful

Number of stateful nodes

Type:

int

nnodes

Total number of nodes

Type:

int

nclasses

Number of job classes

Type:

int

nchains

Number of chains (routing chains)

Type:

int

nclosedjobs

Total number of jobs in closed classes

Type:

int

njobs

(1, K) Population per class (inf for open classes)

Type:

numpy.ndarray

nservers

(M, 1) Number of servers per station

Type:

numpy.ndarray

rates

(M, K) Service rates

Type:

numpy.ndarray

scv

(M, K) Squared coefficient of variation

Type:

numpy.ndarray

visits

Dict[int, ndarray] - Chain ID -> (M, K) visit ratios

Type:

Dict[int, numpy.ndarray]

inchain

Dict[int, ndarray] - Chain ID -> class indices in chain

Type:

Dict[int, numpy.ndarray]

chains

(K, 1) Chain membership per class

Type:

numpy.ndarray

refstat

(K, 1) Reference station per class

Type:

numpy.ndarray

refclass

(1, C) Reference class per chain

Type:

numpy.ndarray

sched

Dict[int, SchedStrategy] - Station ID -> scheduling strategy

Type:

Dict[int, int]

routing

(N, K) Routing strategy matrix

Type:

numpy.ndarray

rt

Routing probability matrix

Type:

numpy.ndarray | None

nodetype

List[NodeType] - Node types

Type:

List[int]

isstation

(N, 1) Boolean mask for stations

Type:

numpy.ndarray

isstateful

(N, 1) Boolean mask for stateful nodes

Type:

numpy.ndarray

nodeToStation

(N, 1) Node index -> station index mapping

Type:

numpy.ndarray

nodeToStateful

(N, 1) Node index -> stateful index mapping

Type:

numpy.ndarray

stationToNode

(M, 1) Station index -> node index mapping

Type:

numpy.ndarray

stationToStateful

(M, 1) Station index -> stateful index mapping

Type:

numpy.ndarray

statefulToNode

(S, 1) Stateful index -> node index mapping

Type:

numpy.ndarray

statefulToStation

(S, 1) Stateful index -> station index mapping

Type:

numpy.ndarray

state

Dict State per stateful node

Type:

Dict[int, numpy.ndarray]

lldscaling

(M, Nmax) Load-dependent scaling matrix

Type:

numpy.ndarray | None

cdscaling

Class-dependent scaling functions

Type:

Dict | None

cap

(M, 1) Station capacities

Type:

numpy.ndarray | None

classcap

(M, K) Per-class capacities

Type:

numpy.ndarray | None

connmatrix

(N, N) Connection matrix

Type:

numpy.ndarray | None

nodenames

List[str] - Node names

Type:

List[str]

classnames

List[str] - Class names

Type:

List[str]

__post_init__()[source]

Ensure arrays are MatrixArray (numpy arrays with .get()/.set() methods).

__repr__()[source]

String representation.

__setattr__(name, value)[source]

Override to convert numpy arrays to MatrixArray for API compatibility.

cap: numpy.ndarray | None = None
cdscaling: Dict | None = None
classcap: numpy.ndarray | None = None
classdeadline: numpy.ndarray | None = None
classprio: numpy.ndarray | None = None
connmatrix: numpy.ndarray | None = None
copy()[source]

Create a deep copy of this NetworkStruct.

csmask: numpy.ndarray | None = None
droprule: Dict | None = None
fj: numpy.ndarray | None = None
get_chain_population(chain_id)[source]

Get total population in a chain.

Parameters:

chain_id (int) – Chain index (0-based)

Returns:

Total number of jobs in the chain

Return type:

float

get_closed_class_indices()[source]

Get indices of closed classes.

get_open_class_indices()[source]

Get indices of open classes.

get_scheduling_at_station(station_id)[source]

Get scheduling strategy at a station.

Parameters:

station_id (int) – Station index (0-based)

Returns:

SchedStrategy value

Return type:

int

get_stateful_indices()[source]

Get indices of stateful nodes.

Returns:

Array of node indices that are stateful

Return type:

numpy.ndarray

get_station_indices()[source]

Get indices of station nodes.

Returns:

Array of node indices that are stations

Return type:

numpy.ndarray

get_total_population()[source]

Get total population across all closed classes.

gsync: Dict | None = None
has_class_dependence()[source]

Check if model has class-dependent scaling.

has_closed_classes()[source]

Check if model has closed (finite population) classes.

has_load_dependence()[source]

Check if model has load-dependent service rates.

has_multi_server()[source]

Check if any station has multiple servers.

has_open_classes()[source]

Check if model has open (infinite population) classes.

is_closed_chain(chain_id)[source]

Check if a chain is closed (finite population).

Parameters:

chain_id (int) – Chain index (0-based)

Returns:

True if chain is closed, False if open

Return type:

bool

is_open_chain(chain_id)[source]

Check if a chain is open (infinite population).

Parameters:

chain_id (int) – Chain index (0-based)

Returns:

True if chain is open, False if closed

Return type:

bool

is_valid()[source]

Check if structure is valid.

Returns:

True if structure passes validation, False otherwise

Return type:

bool

issignal: numpy.ndarray | None = None
isslc: numpy.ndarray | None = None
isstatedep: numpy.ndarray | None = None
ljdcutoffs: numpy.ndarray | None = None
ljdscaling: list | None = None
lldscaling: numpy.ndarray | None = None
lst: Dict | None = None
mu: Dict | None = None
nchains: int = 0
nclasses: int = 0
nclosedjobs: int = 0
nnodes: int = 0
nodeparam: Dict | None = None
nregions: int = 0
nstateful: int = 0
nstations: int = 0
nvars: numpy.ndarray | None = None
property obj

Return self for compatibility with wrapper code that accesses .obj

phases: numpy.ndarray | None = None
phaseshift: numpy.ndarray | None = None
phasessz: numpy.ndarray | None = None
phi: Dict | None = None
pie: Dict | None = None
proc: Dict | None = None
procid: Dict | None = None
region: List | None = None
regionrule: numpy.ndarray | None = None
regionsz: numpy.ndarray | None = None
regionweight: numpy.ndarray | None = None
reward: Dict | None = None
routing_weights: Dict | None = None
rt: numpy.ndarray | None = None
rtnodes: numpy.ndarray | None = None
rtorig: Dict | None = None
schedparam: numpy.ndarray | None = None
signaltype: List | None = None
sync: Dict | None = None
syncreply: numpy.ndarray | None = None
validate()[source]

Validate structural consistency.

Raises:

ValueError – If structural consistency is violated

njobs: numpy.ndarray
nservers: numpy.ndarray
rates: numpy.ndarray
scv: numpy.ndarray
visits: Dict[int, numpy.ndarray]
nodevisits: Dict[int, numpy.ndarray]
inchain: Dict[int, numpy.ndarray]
chains: numpy.ndarray
refstat: numpy.ndarray
refclass: numpy.ndarray
sched: Dict[int, int]
routing: numpy.ndarray
nodetype: List[int]
isstation: numpy.ndarray
isstateful: numpy.ndarray
nodeToStation: numpy.ndarray
nodeToStateful: numpy.ndarray
stationToNode: numpy.ndarray
stationToStateful: numpy.ndarray
statefulToNode: numpy.ndarray
statefulToStation: numpy.ndarray
state: Dict[int, numpy.ndarray]
stateprior: Dict[int, numpy.ndarray]
space: Dict[int, numpy.ndarray]
nodenames: List[str]
classnames: List[str]
class NodeType(*values)[source]

Bases: IntEnum

Node types in a queueing network.

NOTE: Values must match lang/base.py NodeType enum.

static toText(node_type)[source]

Convert node type to text representation.

SOURCE = 0
SINK = 1
QUEUE = 2
DELAY = 3
JOIN = 5
CACHE = 6
ROUTER = 7
CLASSSWITCH = 8
PLACE = 9
TRANSITION = 10
LOGGER = 11
FINITE_CAPACITY_REGION = 12
class SchedStrategy(*values)[source]

Bases: IntEnum

Scheduling strategies.

LCFSPI = 3
HOL = 9
LPS = 17
SETF = 18
class RoutingStrategy(*values)[source]

Bases: IntEnum

Routing strategies.

Values must match MATLAB’s RoutingStrategy constants for JMT compatibility.

RL = 7
class DropStrategy(*values)[source]

Bases: IntEnum

Drop strategies for finite capacity.

WAITQ = 0
DROP = 1
BAS = 2
sn_get_demands_chain(sn)[source]

Calculate new queueing network parameters after aggregating classes into chains.

This function computes chain-level demands, service times, visit ratios, and other parameters by aggregating class-level data based on chain membership.

Parameters:

sn (NetworkStruct) – NetworkStruct object for the queueing network model

Returns:

  • Lchain: (M, C) chain-level demand matrix

  • STchain: (M, C) chain-level service time matrix

  • Vchain: (M, C) chain-level visit ratio matrix

  • alpha: (M, K) class-to-chain weighting matrix

  • Nchain: (1, C) population per chain

  • SCVchain: (M, C) chain-level squared coefficient of variation

  • refstatchain: (C, 1) reference station per chain

Return type:

SnGetDemandsResult containing chain parameters

class SnGetDemandsResult(Lchain, STchain, Vchain, alpha, Nchain, SCVchain, refstatchain)[source]

Bases: object

Result of sn_get_demands_chain calculation.

Lchain

(M, C) Chain-level demand matrix

Type:

numpy.ndarray

STchain

(M, C) Chain-level service time matrix

Type:

numpy.ndarray

Vchain

(M, C) Chain-level visit ratio matrix

Type:

numpy.ndarray

alpha

(M, K) Class-to-chain weighting matrix

Type:

numpy.ndarray

Nchain

(1, C) Population per chain

Type:

numpy.ndarray

SCVchain

(M, C) Chain-level squared coefficient of variation

Type:

numpy.ndarray

refstatchain

(C, 1) Reference station per chain

Type:

numpy.ndarray

Lchain: numpy.ndarray
STchain: numpy.ndarray
Vchain: numpy.ndarray
alpha: numpy.ndarray
Nchain: numpy.ndarray
SCVchain: numpy.ndarray
refstatchain: numpy.ndarray
sn_deaggregate_chain_results(sn, Lchain, ST, STchain, Vchain, alpha, Qchain, Uchain, Rchain, Tchain, Cchain, Xchain)[source]

Calculate class-based performance metrics from chain-level performance measures.

This function disaggregates chain-level performance metrics (queue lengths, utilizations, response times, throughputs) to class-level metrics using the aggregation factors (alpha).

Parameters:

  • sn (NetworkStruct) – NetworkStruct object for the queueing network model

  • Lchain (numpy.ndarray) – (M, C) Service demands per chain

  • ST (numpy.ndarray | None) – (M, K) Mean service times per class (optional, computed from rates if None)

  • STchain (numpy.ndarray) – (M, C) Mean service times per chain

  • Vchain (numpy.ndarray) – (M, C) Mean visits per chain

  • alpha (numpy.ndarray) – (M, K) Class aggregation coefficients

  • Qchain (numpy.ndarray | None) – (M, C) Mean queue-lengths per chain (optional)

  • Uchain (numpy.ndarray | None) – (M, C) Mean utilization per chain (optional)

  • Rchain (numpy.ndarray) – (M, C) Mean response time per chain

  • Tchain (numpy.ndarray) – (M, C) Mean throughput per chain

  • Cchain (numpy.ndarray | None) – (1, C) Mean system response time per chain (optional, not yet supported)

  • Xchain (numpy.ndarray) – (1, C) Mean system throughput per chain

Returns:

  • Q: (M, K) queue lengths

  • U: (M, K) utilizations

  • R: (M, K) response times

  • T: (M, K) throughputs

  • C: (1, K) system response times

  • X: (1, K) system throughputs

Return type:

SnDeaggregateResult containing class-level performance metrics

Raises:

RuntimeError – If Cchain is provided (not yet supported)

class SnDeaggregateResult(Q, U, R, T, C, X)[source]

Bases: object

Result of sn_deaggregate_chain_results calculation.

Q

(M, K) Class-level queue lengths

Type:

numpy.ndarray

U

(M, K) Class-level utilizations

Type:

numpy.ndarray

R

(M, K) Class-level response times

Type:

numpy.ndarray

T

(M, K) Class-level throughputs

Type:

numpy.ndarray

C

(1, K) Class-level system response times

Type:

numpy.ndarray

X

(1, K) Class-level system throughputs

Type:

numpy.ndarray

Q: numpy.ndarray
U: numpy.ndarray
R: numpy.ndarray
T: numpy.ndarray
C: numpy.ndarray
X: numpy.ndarray
class ProductFormParams(lam, D, N, Z, mu, S, V)[source]

Bases: NamedTuple

Result of sn_get_product_form_params calculation.

Create new instance of ProductFormParams(lam, D, N, Z, mu, S, V)

D: numpy.ndarray

Alias for field number 1

N: numpy.ndarray

Alias for field number 2

S: numpy.ndarray

Alias for field number 5

V: numpy.ndarray

Alias for field number 6

Z: numpy.ndarray

Alias for field number 3

lam: numpy.ndarray

Alias for field number 0

mu: numpy.ndarray

Alias for field number 4

sn_get_product_form_params(sn)[source]

Extract standard product-form parameters from the network structure.

This function extracts class-level parameters from a network structure for use in product-form queueing network analysis.

Parameters:

sn (NetworkStruct) – NetworkStruct object

Returns:

  • lam: Arrival rates for open classes

  • D: Service demands at queueing stations

  • N: Population vector

  • Z: Think times (service demands at delay stations)

  • mu: Load-dependent service capacity scaling factors

  • S: Number of servers at queueing stations

  • V: Visit ratios

Return type:

ProductFormParams containing

References

MATLAB: matlab/src/api/sn/sn_get_product_form_params.m

sn_get_residt_from_respt(sn, RN, WH=None)[source]

Compute residence times from response times.

This function converts response times to residence times by accounting for visit ratios at each station.

Parameters:

  • sn (NetworkStruct) – NetworkStruct object

  • RN (numpy.ndarray) – Average response times (M, K)

  • WH (Dict | None) – Residence time handles (optional)

Returns:

Average residence times (M, K)

Return type:

WN

References

MATLAB: matlab/src/api/sn/sn_get_residt_from_respt.m

sn_get_state_aggr(sn)[source]

Get aggregated state representation.

Parameters:

sn (NetworkStruct) – NetworkStruct object

Returns:

Dictionary mapping stateful node index to aggregated state

Return type:

Dict[int, numpy.ndarray]

References

MATLAB: matlab/src/api/sn/sn_get_state_aggr.m

sn_set_arrival(sn, station_idx, class_idx, rate)[source]

Set arrival rate for a class at a station.

Parameters:

  • sn (NetworkStruct) – NetworkStruct object (modified in place)

  • station_idx (int) – Station index (0-based)

  • class_idx (int) – Class index (0-based)

  • rate (float) – Arrival rate

References

MATLAB: matlab/src/api/sn/sn_set_arrival.m

sn_set_service(sn, station_idx, class_idx, rate, scv=1.0)[source]

Set service rate for a class at a station.

Parameters:

  • sn (NetworkStruct) – NetworkStruct object (modified in place)

  • station_idx (int) – Station index (0-based)

  • class_idx (int) – Class index (0-based)

  • rate (float) – Service rate

  • scv (float) – Squared coefficient of variation (default 1.0 for exponential)

References

MATLAB: matlab/src/api/sn/sn_set_service.m

sn_set_servers(sn, station_idx, nservers)[source]

Set number of servers at a station.

Parameters:

  • sn (NetworkStruct) – NetworkStruct object (modified in place)

  • station_idx (int) – Station index (0-based)

  • nservers (int) – Number of servers

References

MATLAB: matlab/src/api/sn/sn_set_servers.m

sn_set_population(sn, class_idx, njobs)[source]

Set population for a class.

Parameters:

  • sn (NetworkStruct) – NetworkStruct object (modified in place)

  • class_idx (int) – Class index (0-based)

  • njobs (float) – Number of jobs (inf for open class)

References

MATLAB: matlab/src/api/sn/sn_set_population.m

sn_set_priority(sn, class_idx, priority)[source]

Set priority for a class.

Parameters:

  • sn (NetworkStruct) – NetworkStruct object (modified in place)

  • class_idx (int) – Class index (0-based)

  • priority (int) – Priority level (higher = more priority)

References

MATLAB: matlab/src/api/sn/sn_set_priority.m

sn_set_routing(sn, source_node, dest_node, source_class, dest_class, prob)[source]

Set routing probability between nodes and classes.

Parameters:

  • sn (NetworkStruct) – NetworkStruct object (modified in place)

  • source_node (int) – Source node index (0-based)

  • dest_node (int) – Destination node index (0-based)

  • source_class (int) – Source class index (0-based)

  • dest_class (int) – Destination class index (0-based)

  • prob (float) – Routing probability

References

MATLAB: matlab/src/api/sn/sn_set_routing.m

sn_refresh_visits(sn)[source]

Refresh visit ratios from routing matrix.

This function solves traffic equations to compute visit ratios at each station from the routing probability matrix.

Parameters:

sn (NetworkStruct) – NetworkStruct object (modified in place)

References

MATLAB: matlab/src/api/sn/sn_refresh_visits.m

sn_set_fork_fanout(sn, fork_node_idx, fan_out)[source]

Set fork fanout (tasksPerLink) for a Fork node.

Updates the fanOut field in nodeparam for a Fork node.

Parameters:

  • sn (NetworkStruct) – NetworkStruct object

  • fork_node_idx (int) – Node index of the Fork node (0-based)

  • fan_out (int) – Number of tasks per output link (>= 1)

Returns:

Modified NetworkStruct

Raises:

ValueError – If the specified node is not a Fork node

Return type:

NetworkStruct

References

MATLAB: matlab/src/api/sn/sn_set_fork_fanout.m

sn_set_service_batch(sn, rates, scvs=None, auto_refresh=False)[source]

Set service rates for multiple station-class pairs.

Batch update of service rates. NaN values are skipped (not updated). More efficient than calling sn_set_service multiple times.

Parameters:

  • sn (NetworkStruct) – NetworkStruct object

  • rates (numpy.ndarray) – Matrix of new rates (nstations x nclasses), NaN = skip

  • scvs (numpy.ndarray | None) – Matrix of new SCVs (optional)

  • auto_refresh (bool) – If True, refresh process fields (default False)

Returns:

Modified NetworkStruct

Return type:

NetworkStruct

References

MATLAB: matlab/src/api/sn/sn_set_service_batch.m

sn_nonmarkov_toph(sn, options=None)[source]

Convert non-Markovian distributions to Phase-Type using approximation.

This function scans all service and arrival processes in the network structure and converts non-Markovian distributions to Markovian Arrival Processes (MAPs) using the specified approximation method.

Supported non-Markovian distributions: - GAMMA: Gamma distribution - WEIBULL: Weibull distribution - LOGNORMAL: Lognormal distribution - PARETO: Pareto distribution - UNIFORM: Uniform distribution - DET: Deterministic (converted to Erlang)

Parameters:

  • sn (NetworkStruct) – NetworkStruct object (from getStruct())

  • options (Dict[str, Any] | None) – Solver options dict with fields: - config.nonmkv: Method for conversion (‘none’, ‘bernstein’) - config.nonmkvorder: Number of phases for approximation (default 20) - config.preserveDet: Keep deterministic distributions (for MAP/D/c)

Returns:

Modified NetworkStruct with converted processes

Return type:

NetworkStruct

References

MATLAB: matlab/src/api/sn/sn_nonmarkov_toph.m

class ChainParams(lambda_vec, D, N, Z, mu, S, V)[source]

Bases: object

Chain-aggregated product-form parameters.

lambda_vec: numpy.ndarray
D: numpy.ndarray
N: numpy.ndarray
Z: numpy.ndarray
mu: numpy.ndarray
S: numpy.ndarray
V: numpy.ndarray
sn_get_arvr_from_tput(sn, TN, TH=None)[source]

Compute average arrival rates at stations from throughputs.

Calculates the average arrival rate at each station in steady-state from the station throughputs and routing matrix.

Parameters:
Returns:

Average arrival rates at stations (M x R)

Return type:

AN

References

Original MATLAB: matlab/src/api/sn/sn_get_arvr_from_tput.m

sn_get_node_arvr_from_tput(sn, TN, TH=None, AN=None)[source]

Compute node arrival rates from station throughputs.

This function handles: - Station nodes: Uses station arrival rates directly - Cache nodes: Only requesting classes arrive (not hit/miss classes) - Non-station nodes (ClassSwitch, Sink): Uses nodevisits-based computation

Parameters:
Returns:

Node arrival rates (I x R)

Return type:

ANn

References

Original MATLAB: matlab/src/api/sn/sn_get_node_arvr_from_tput.m

sn_get_node_tput_from_tput(sn, TN, TH=None, ANn=None)[source]

Compute node throughputs from station throughputs.

This function handles: - Station nodes: Uses station throughputs directly - Cache nodes: Uses actual hit/miss probabilities if available - Non-station nodes: Uses routing matrix (rtnodes) for computation

Parameters:
Returns:

Node throughputs (I x R)

Return type:

TNn

References

Original MATLAB: matlab/src/api/sn/sn_get_node_tput_from_tput.m

sn_get_product_form_chain_params(sn)[source]

Extract product-form parameters aggregated by chain.

Extracts parameters from a network structure and aggregates them by chain for use in product-form analysis methods.

Parameters:

sn (NetworkStruct) – Network structure

Returns:

ChainParams with lambda_vec, D, N, Z, mu, S, V

Return type:

ChainParams

References

Original MATLAB: matlab/src/api/sn/sn_get_product_form_chain_params.m

sn_set_routing_prob(sn, from_stateful, from_class, to_stateful, to_class, prob, auto_refresh=False)[source]

Set a routing probability between two stateful node-class pairs.

Updates a single entry in the rt matrix.

Parameters:

  • sn (NetworkStruct) – Network structure

  • from_stateful (int) – Source stateful node index (0-based)

  • from_class (int) – Source class index (0-based)

  • to_stateful (int) – Destination stateful node index (0-based)

  • to_class (int) – Destination class index (0-based)

  • prob (float) – Routing probability [0, 1]

  • auto_refresh (bool) – If True, refresh visit ratios (default False)

Returns:

Modified network structure

Return type:

NetworkStruct

References

Original MATLAB: matlab/src/api/sn/sn_set_routing_prob.m

sn_is_closed_model(sn)[source]

Check if the network model is closed (all finite populations).

A closed model has all finite job populations.

Parameters:

sn (NetworkStruct) – NetworkStruct object

Returns:

True if the network is a closed model

Return type:

bool

sn_is_open_model(sn)[source]

Check if the network model is open (all infinite populations).

An open model has only infinite (open) job classes.

Parameters:

sn (NetworkStruct) – NetworkStruct object

Returns:

True if the network is an open model

Return type:

bool

sn_is_mixed_model(sn)[source]

Check if the network model is mixed (both open and closed classes).

Parameters:

sn (NetworkStruct) – NetworkStruct object

Returns:

True if the network has both open and closed classes

Return type:

bool

sn_is_population_model(sn)[source]

Check if the network model is a population model.

A population model uses only delay-like scheduling strategies (INF, PS, PSPRIO, DPS, GPS, GPSPRIO, DPSPRIO, EXT), has no priorities, and no fork-join topology.

Parameters:

sn (NetworkStruct) – NetworkStruct object

Returns:

True if model is population-based

Return type:

bool

sn_has_closed_classes(sn)[source]

Check if the network has closed (finite population) classes.

Parameters:

sn (NetworkStruct) – NetworkStruct object

Returns:

True if network has at least one closed class

Return type:

bool

sn_has_open_classes(sn)[source]

Check if the network has open (infinite population) classes.

Parameters:

sn (NetworkStruct) – NetworkStruct object

Returns:

True if network has at least one open class

Return type:

bool

sn_has_mixed_classes(sn)[source]

Check if the network has both open and closed classes.

Parameters:

sn (NetworkStruct) – NetworkStruct object

Returns:

True if network has both open and closed classes

Return type:

bool

sn_has_single_class(sn)[source]

Check if the network has exactly one class.

Parameters:

sn (NetworkStruct) – NetworkStruct object

Returns:

True if network has exactly one class

Return type:

bool

sn_has_multi_class(sn)[source]

Check if the network has multiple classes.

Parameters:

sn (NetworkStruct) – NetworkStruct object

Returns:

True if network has more than one class

Return type:

bool

sn_has_multiple_closed_classes(sn)[source]

Check if the network has multiple closed classes.

Parameters:

sn (NetworkStruct) – NetworkStruct object

Returns:

True if network has more than one closed class

Return type:

bool

sn_has_single_chain(sn)[source]

Check if the network has exactly one chain.

Parameters:

sn (NetworkStruct) – NetworkStruct object

Returns:

True if network has exactly one chain

Return type:

bool

sn_has_multi_chain(sn)[source]

Check if the network has multiple chains.

Parameters:

sn (NetworkStruct) – NetworkStruct object

Returns:

True if network has more than one chain

Return type:

bool

sn_has_fcfs(sn)[source]

Check if the network has any FCFS (First-Come First-Served) stations.

Parameters:

sn (NetworkStruct) – NetworkStruct object

Returns:

True if network has at least one FCFS station

Return type:

bool

sn_has_ps(sn)[source]

Check if the network has any PS (Processor Sharing) stations.

Parameters:

sn (NetworkStruct) – NetworkStruct object

Returns:

True if network has at least one PS station

Return type:

bool

sn_has_inf(sn)[source]

Check if the network has any INF (Infinite Server/Delay) stations.

Parameters:

sn (NetworkStruct) – NetworkStruct object

Returns:

True if network has at least one INF station

Return type:

bool

sn_has_lcfs(sn)[source]

Check if the network has any LCFS (Last-Come First-Served) stations.

Parameters:

sn (NetworkStruct) – NetworkStruct object

Returns:

True if network has at least one LCFS station

Return type:

bool

sn_has_lcfspr(sn)[source]

Check if the network has any LCFS-PR (LCFS Preemptive Resume) stations.

Parameters:

sn (NetworkStruct) – NetworkStruct object

Returns:

True if network has at least one LCFS-PR station

Return type:

bool

sn_has_lcfs_pr(sn)[source]

Check if the network has any LCFS-PR (LCFS Preemptive Resume) stations.

This is an alias for sn_has_lcfspr, matching the MATLAB function name.

Parameters:

sn (NetworkStruct) – NetworkStruct object

Returns:

True if network has at least one LCFS-PR station

Return type:

bool

sn_has_lcfs_pi(sn)[source]

Check if the network has any LCFS-PI (LCFS Preemptive Identical) stations.

Parameters:

sn (NetworkStruct) – NetworkStruct object

Returns:

True if network has at least one LCFS-PI station

Return type:

bool

sn_has_siro(sn)[source]

Check if the network has any SIRO (Service In Random Order) stations.

Parameters:

sn (NetworkStruct) – NetworkStruct object

Returns:

True if network has at least one SIRO station

Return type:

bool

sn_has_dps(sn)[source]

Check if the network has any DPS (Discriminatory Processor Sharing) stations.

Parameters:

sn (NetworkStruct) – NetworkStruct object

Returns:

True if network has at least one DPS station

Return type:

bool

sn_has_dps_prio(sn)[source]

Check if the network has any DPS with priority stations.

Parameters:

sn (NetworkStruct) – NetworkStruct object

Returns:

True if network has at least one DPS-PRIO station

Return type:

bool

sn_has_gps(sn)[source]

Check if the network has any GPS (Generalized Processor Sharing) stations.

Parameters:

sn (NetworkStruct) – NetworkStruct object

Returns:

True if network has at least one GPS station

Return type:

bool

sn_has_gps_prio(sn)[source]

Check if the network has any GPS with priority stations.

Parameters:

sn (NetworkStruct) – NetworkStruct object

Returns:

True if network has at least one GPS-PRIO station

Return type:

bool

sn_has_ps_prio(sn)[source]

Check if the network has any PS with priority stations.

Parameters:

sn (NetworkStruct) – NetworkStruct object

Returns:

True if network has at least one PS-PRIO station

Return type:

bool

sn_has_hol(sn)[source]

Check if the network has any HOL (Head of Line) priority stations.

Parameters:

sn (NetworkStruct) – NetworkStruct object

Returns:

True if network has at least one HOL station

Return type:

bool

sn_has_lps(sn)[source]

Check if the network has any LPS (Least Progress Scheduling) stations.

Parameters:

sn (NetworkStruct) – NetworkStruct object

Returns:

True if network has at least one LPS station

Return type:

bool

sn_has_setf(sn)[source]

Check if the network has any SETF (Shortest Elapsed Time First) stations.

Parameters:

sn (NetworkStruct) – NetworkStruct object

Returns:

True if network has at least one SETF station

Return type:

bool

sn_has_sept(sn)[source]

Check if the network has any SEPT (Shortest Expected Processing Time) stations.

Parameters:

sn (NetworkStruct) – NetworkStruct object

Returns:

True if network has at least one SEPT station

Return type:

bool

sn_has_lept(sn)[source]

Check if the network has any LEPT (Longest Expected Processing Time) stations.

Parameters:

sn (NetworkStruct) – NetworkStruct object

Returns:

True if network has at least one LEPT station

Return type:

bool

sn_has_sjf(sn)[source]

Check if the network has any SJF (Shortest Job First) stations.

Parameters:

sn (NetworkStruct) – NetworkStruct object

Returns:

True if network has at least one SJF station

Return type:

bool

sn_has_ljf(sn)[source]

Check if the network has any LJF (Longest Job First) stations.

Parameters:

sn (NetworkStruct) – NetworkStruct object

Returns:

True if network has at least one LJF station

Return type:

bool

sn_has_polling(sn)[source]

Check if the network has any polling stations.

Parameters:

sn (NetworkStruct) – NetworkStruct object

Returns:

True if network has at least one polling station

Return type:

bool

sn_has_homogeneous_scheduling(sn, strategy)[source]

Check if the network uses an identical scheduling strategy at every station.

Parameters:

  • sn (NetworkStruct) – NetworkStruct object

  • strategy (int) – SchedStrategy value to check for

Returns:

True if all stations use the specified strategy

Return type:

bool

sn_has_multi_class_fcfs(sn)[source]

Check if the network has an FCFS station that serves multiple classes.

Parameters:

sn (NetworkStruct) – NetworkStruct object

Returns:

True if any FCFS station serves more than one class

Return type:

bool

sn_has_multi_class_heter_fcfs(sn)[source]

Check if network has multiclass heterogeneous FCFS stations.

A heterogeneous FCFS station has different service rates for different classes. Uses MATLAB’s range() check: max(rates) - min(rates) > 0 across all classes at each FCFS station.

Parameters:

sn (NetworkStruct) – NetworkStruct object

Returns:

True if network has FCFS stations with heterogeneous class rates

Return type:

bool

sn_has_multi_class_heter_exp_fcfs(sn)[source]

Check if network has multiclass heterogeneous exponential FCFS stations.

Returns true if any FCFS station has heterogeneous rates AND all service time SCVs at that station are approximately 1.0 (exponential).

Parameters:

sn (NetworkStruct) – NetworkStruct object

Returns:

True if network has FCFS stations with heterogeneous exponential service

Return type:

bool

sn_has_multi_server(sn)[source]

Check if the network has any multi-server stations.

Parameters:

sn (NetworkStruct) – NetworkStruct object

Returns:

True if any station has more than one server

Return type:

bool

sn_has_load_dependence(sn)[source]

Check if the network has load-dependent service.

Parameters:

sn (NetworkStruct) – NetworkStruct object

Returns:

True if network has load-dependent scaling

Return type:

bool

sn_has_joint_dependence(sn)[source]

Check if the network has joint-dependent service rates.

This includes both LJD (Limited Joint Dependence) and LJCD (Limited Joint Class Dependence) scaling.

Parameters:

sn (NetworkStruct) – NetworkStruct object

Returns:

True if network has joint-dependent scaling

Return type:

bool

sn_has_fork_join(sn)[source]

Check if the network uses fork and/or join nodes.

Parameters:

sn (NetworkStruct) – NetworkStruct object

Returns:

True if network has fork-join topology

Return type:

bool

sn_has_priorities(sn)[source]

Check if the network uses class priorities.

In LINE, priority 0 is default (no priority). Values > 0 indicate priority classes are in use.

Parameters:

sn (NetworkStruct) – NetworkStruct object

Returns:

True if any class has priority > 0

Return type:

bool

sn_has_class_switching(sn)[source]

Check if the network has class switching.

Class switching is indicated by the number of classes differing from the number of chains.

Parameters:

sn (NetworkStruct) – NetworkStruct object

Returns:

True if number of classes differs from number of chains

Return type:

bool

sn_has_fractional_populations(sn)[source]

Check if the network has fractional (non-integer) populations.

Parameters:

sn (NetworkStruct) – NetworkStruct object

Returns:

True if any class has fractional population

Return type:

bool

sn_has_sd_routing(sn)[source]

Check if the network has state-dependent routing strategies.

State-dependent routing strategies violate the product-form assumption. These include Round-Robin, Weighted Round-Robin, Join Shortest Queue, Power of K Choices, and Reinforcement Learning.

Product-form requires state-independent (Markovian) routing. PROB and RAND are product-form compatible.

Parameters:

sn (NetworkStruct) – NetworkStruct object

Returns:

True if network has state-dependent routing

Return type:

bool

sn_has_product_form(sn)[source]

Check if the network has a known product-form solution.

A network has product form if: - All stations use INF, PS, FCFS, LCFS-PR, or EXT scheduling - No multiclass heterogeneous FCFS - No priorities - No fork-join

Parameters:

sn (NetworkStruct) – NetworkStruct object

Returns:

True if network has product-form solution

Return type:

bool

sn_has_product_form_not_het_fcfs(sn)[source]

Check if network has product form except for heterogeneous FCFS.

This checks: - All stations use INF, PS, FCFS, LCFSPR, or EXT scheduling - No priorities, no fork-join, no state-dependent routing - At FCFS stations, all active class SCVs are approximately 1 (exponential)

Parameters:

sn (NetworkStruct) – NetworkStruct object

Returns:

True if network would have product form without heterogeneous FCFS

Return type:

bool

sn_has_product_form_except_multi_class_heter_exp_fcfs(sn)[source]

Check if network has product form except for multiclass heterogeneous exponential FCFS.

Parameters:

sn (NetworkStruct) – NetworkStruct object

Returns:

True if network would have product form without multiclass heter exp FCFS

Return type:

bool

sn_is_state_valid(sn)[source]

Check if the network state is valid.

Parameters:

sn (NetworkStruct) – NetworkStruct object

Returns:

True if state is valid

Return type:

bool

sn_print(sn, file=None)[source]

Print comprehensive information about a NetworkStruct object.

This function displays all fields, matrices, lists, and maps in a formatted manner useful for debugging and inspection of network structures.

Parameters:

  • sn (NetworkStruct) – Network structure to inspect

  • file – Output file (default: sys.stdout)

References

MATLAB: matlab/src/api/sn/sn_print.m

sn_print_routing_matrix(sn, onlyclass=None, file=None)[source]

Print the routing matrix of the network.

This function displays the routing probabilities between nodes and classes in a human-readable format.

Parameters:

  • sn (NetworkStruct) – Network structure

  • onlyclass (Any | None) – Optional filter for a specific class (object with ‘name’ attribute)

  • file – Output file (default: sys.stdout)

References

MATLAB: matlab/src/api/sn/sn_print_routing_matrix.m

sn_refresh_process_fields(sn, station_idx, class_idx)[source]

Refresh process fields based on rate and SCV values.

Updates mu, phi, proc, pie, phases based on current rate and SCV values. - SCV = 1.0: Exponential (1 phase) - SCV < 1.0: Erlang approximation - SCV > 1.0: Hyperexponential(2) approximation

Parameters:

  • sn (NetworkStruct) – Network structure (modified in place)

  • station_idx (int) – Station index (0-based)

  • class_idx (int) – Class index (0-based)

Returns:

Modified network structure

Return type:

NetworkStruct

References

MATLAB: matlab/src/api/sn/sn_refresh_process_fields.m

sn_rtnodes_to_rtorig(sn)[source]

Convert node routing matrix to the original routing matrix format.

This function converts the node-level routing matrix to the original routing matrix format, excluding class-switching nodes.

Parameters:

sn (NetworkStruct) – Network structure

Returns:

rtorigcell: Dictionary representation {(r,s): ndarray} rtorig: Sparse/dense matrix representation

Return type:

Tuple of (rtorigcell, rtorig) where

References

MATLAB: matlab/src/api/sn/sn_rtnodes_to_rtorig.m

Non-Product-Form Networks (line_solver.api.npfqn)

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

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

Key algorithms:

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

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

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

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

Parameters:

  • method (str) – Approximation method (“default”, “none”, “hmva”, “interp”)

  • sn (NetworkStruct) – Network structure

  • ST (numpy.ndarray) – Service time matrix (M x K)

  • V (numpy.ndarray | None) – Visit ratios matrix (M x K), optional

  • SCV (numpy.ndarray) – Squared coefficient of variation matrix (M x K)

  • Tin (numpy.ndarray) – Initial throughput matrix (M x K)

  • Uin (numpy.ndarray) – Initial utilization matrix (M x K)

  • gamma (numpy.ndarray) – Gamma correction matrix (M x 1)

  • nservers (numpy.ndarray) – Number of servers matrix (M x 1)

Returns:

NpfqnNonexpApproxResult with updated matrices

Raises:

ValueError – If unknown approximation method is specified

Return type:

NpfqnNonexpApproxResult

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

Bases: object

Result of non-exponential approximation.

ST: numpy.ndarray
gamma: numpy.ndarray
nservers: numpy.ndarray
rho: numpy.ndarray
scva: numpy.ndarray
scvs: numpy.ndarray
eta: numpy.ndarray
npfqn_traffic_merge(MMAPa, config_merge='default', config_compress=None)[source]

Merge multiple MMAP traffic flows.

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

Parameters:

  • MMAPa (Dict[int, List[numpy.ndarray] | None]) – Dictionary of MMAP traffic flows to be merged. Keys are integer indices, values are MMAP lists [D0, D1, …].

  • config_merge (str) – Merge configuration. Options: - “default”, “super”: Use MMAP superposition - “mixture”: Apply mixture fitting after superposition - “interpos”: Interposition method (falls back to super)

  • config_compress (str | None) – Compression configuration. Options: - None, “none”: No compression - “default”: Apply default compression

Returns:

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

Raises:

RuntimeError – If unsupported merge configuration is provided.

Return type:

List[numpy.ndarray] | None

Example

>>> mmap1 = [np.array([[-1.0]]), np.array([[1.0]])]
>>> mmap2 = [np.array([[-2.0]]), np.array([[2.0]])]
>>> result = npfqn_traffic_merge({0: mmap1, 1: mmap2})
npfqn_traffic_merge_cs(MMAPs, prob, config='default')[source]

Merge MMAP traffic flows with class switching.

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

Parameters:

  • MMAPs (Dict[int, List[numpy.ndarray]]) – Dictionary of MMAP traffic flows indexed by source.

  • prob (numpy.ndarray) –

    Class switching probability matrix ((n*R) x R) where: - n is the number of sources - R is the number of classes - prob[(i-1)*R + r, s] = probability that class r from source i

    becomes class s in the merged stream

  • config (str) – Merge configuration (“default” or “super”).

Returns:

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

Return type:

List[numpy.ndarray] | None

Algorithm:

  1. Apply mmap_mark to each source MMAP to encode class switching

  2. Superpose all marked MMAPs

  3. Return result

Example

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

Split MMAP traffic flows with class switching.

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

Parameters:

  • MMAP_input (List[numpy.ndarray]) – Input MMAP as list [D0, D1, D2, …].

  • P (numpy.ndarray) –

    Class switching probability matrix (R x J) where: - R is the number of arrival classes - J = M * R where M is the number of destinations - P[r, (jst-1)*R + s] = probability that class r arrival

    becomes class s at destination jst

Returns:

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

Return type:

Dict[int, List[numpy.ndarray] | None]

Algorithm:

  1. Parse dimensions from P matrix

  2. For each destination, create MMAP with weighted marking matrices

  3. Normalize each result

Example

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

Polling Systems (line_solver.api.polling)

Polling System Analysis Algorithms.

Native Python implementations for analyzing polling/vacation queue systems with various disciplines.

Key algorithms:

polling_qsys_exhaustive: Exhaustive polling discipline polling_qsys_gated: Gated polling discipline polling_qsys_1limited: 1-Limited polling discipline

polling_qsys_exhaustive(arvMAPs, svcMAPs, switchMAPs)[source]

Compute exact mean waiting times for exhaustive polling system.

In exhaustive polling, the server continues to serve a queue until it becomes empty before moving to the next queue.

Based on Takagi, ACM Computing Surveys, Vol. 20, No. 1, 1988, eq (15).

Parameters:
Returns:

Array of mean waiting times for each queue.

Return type:

numpy.ndarray

polling_qsys_gated(arvMAPs, svcMAPs, switchMAPs)[source]

Compute exact mean waiting times for gated polling system.

In gated polling, the server serves all customers present at the beginning of a visit period.

Based on Takagi, ACM Computing Surveys, Vol. 20, No. 1, 1988, eq (20).

Parameters:
Returns:

Array of mean waiting times for each queue.

Return type:

numpy.ndarray

polling_qsys_1limited(arvMAPs, svcMAPs, switchMAPs)[source]

Compute exact mean waiting times for 1-limited polling system.

In 1-limited polling, the server serves at most one customer from each queue before moving to the next queue.

Based on Takagi, ACM Computing Surveys, Vol. 20, No. 1, 1988, eq (20).

Parameters:
Returns:

Array of mean waiting times for each queue.

Return type:

numpy.ndarray

Loss Networks (line_solver.api.lossn)

Loss Network Analysis Algorithms.

Native Python implementations for analyzing loss networks using Erlang formulas and related methods.

Key algorithms:

lossn_erlangfp: Erlang fixed-point algorithm for loss networks erlang_b: Erlang B blocking probability erlang_c: Erlang C delay probability

lossn_erlangfp(nu, A, c, tol=1e-8, max_iter=1000)[source]

Erlang fixed point approximation for loss networks.

Calls (jobs) on route (class) r arrive according to Poisson rate nu_r. Call service times on route r have unit mean.

The link capacity requirements are:

sum_r A[j,r] * n[j,r] < c[j]

for all links j, where n[j,r] counts calls on route r on link j.

Parameters:

  • nu (numpy.ndarray) – Arrival rates vector (R,) for each route.

  • A (numpy.ndarray) – Capacity requirement matrix (J, R) - A[j,r] is capacity required on link j by route r.

  • c (numpy.ndarray) – Capacity vector (J,) - c[j] is capacity of link j.

  • tol (float) – Convergence tolerance.

  • max_iter (int) – Maximum iterations.

Returns:

  • qlen: Mean queue-length for each route (R,)

  • loss: Loss probability for each route (R,)

  • eblock: Blocking probability for each link (J,)

  • niter: Number of iterations

Return type:

Tuple of (qlen, loss, eblock, niter) where

Example

>>> nu = np.array([0.3, 0.1])
>>> A = np.array([[1, 1], [1, 4]])  # 2 links, 2 routes
>>> c = np.array([1, 3])
>>> qlen, loss, eblock, niter = lossn_erlangfp(nu, A, c)
erlang_b(offered_load, servers)[source]

Compute Erlang B blocking probability.

The Erlang B formula gives the probability that an arriving call is blocked in an M/M/c/c loss system.

Parameters:

  • offered_load (float) – Traffic intensity (arrival rate * service time).

  • servers (int) – Number of servers.

Returns:

Blocking probability.

Return type:

float

Example

>>> erlang_b(10.0, 12)  # Offered load 10 Erlang, 12 channels
0.1054...
erlang_c(offered_load, servers)[source]

Compute Erlang C delay probability.

The Erlang C formula gives the probability that an arriving call must wait in an M/M/c queue.

Parameters:

  • offered_load (float) – Traffic intensity (arrival rate * service time).

  • servers (int) – Number of servers.

Returns:

Delay probability (probability of waiting).

Return type:

float

Example

>>> erlang_c(10.0, 12)  # Offered load 10 Erlang, 12 agents

Layered Stochastic Networks (line_solver.api.lsn)

Layered Stochastic Network (LSN) utilities.

Native Python implementations for layered queueing network analysis.

Key functions:

lsn_max_multiplicity: Compute maximum multiplicity for tasks in a layered network.

Key classes:

LayeredNetworkStruct: Structure representing a layered network. LayeredNetworkElement: Enumeration of layered network element types.

class LayeredNetworkElement(*values)[source]

Bases: IntEnum

Types of elements in a layered network.

TASK = 1
ENTRY = 2
ACTIVITY = 3
PROCESSOR = 4
HOST = 5
class LayeredNetworkStruct(dag, mult, type, isref)[source]

Bases: object

Structure representing a layered network for LSN analysis.

dag

Directed acyclic graph adjacency matrix (n x n). dag[i,j] > 0 indicates an edge from node i to node j.

Type:

numpy.ndarray

mult

Multiplicity (max concurrent instances) for each node (n,).

Type:

numpy.ndarray

type

Node type for each node (n,), using LayeredNetworkElement values.

Type:

numpy.ndarray

isref

Reference task flags (n,). Non-zero indicates a reference task.

Type:

numpy.ndarray

dag: numpy.ndarray
mult: numpy.ndarray
type: numpy.ndarray
isref: numpy.ndarray
kahn_topological_sort(adjacency)[source]

Perform Kahn’s algorithm for topological sorting.

Parameters:

adjacency (numpy.ndarray) – Adjacency matrix where adjacency[i,j] > 0 means edge i -> j.

Returns:

List of node indices in topological order.

Raises:

ValueError – If the graph contains a cycle.

Return type:

List[int]

lsn_max_multiplicity(lsn)[source]

Compute the maximum multiplicity for each task in a layered network.

This function uses flow analysis based on Kahn’s topological sorting algorithm to determine the maximum sustainable throughput for each task, considering both the incoming flow and the multiplicity constraints.

Parameters:

lsn (LayeredNetworkStruct) – The layered network structure containing task dependencies and constraints.

Returns:

Matrix of maximum multiplicities for each task in the network (n x 1).

Return type:

numpy.ndarray

Algorithm:

  1. Build binary adjacency graph from DAG

  2. Apply Kahn’s topological sort to determine processing order

  3. Initialize inflow from reference tasks

  4. For each node in topological order: - outflow = min(inflow, multiplicity constraint) - Propagate outflow to downstream nodes

  5. Handle unreachable tasks (infinite multiplicity)

Example

>>> lsn = LayeredNetworkStruct(
...     dag=np.array([[0, 1, 0], [0, 0, 1], [0, 0, 0]]),
...     mult=np.array([2, 3, 5]),
...     type=np.array([1, 1, 1]),
...     isref=np.array([1, 0, 0])
... )
>>> max_mult = lsn_max_multiplicity(lsn)

Trace Analysis (line_solver.api.trace)

Trace Analysis Functions.

Native Python implementations for statistical analysis of empirical trace data including means, variances, correlations, and index of dispersion.

Key functions:

trace_mean: Mean of trace data trace_var: Variance of trace data trace_scv: Squared coefficient of variation trace_acf: Autocorrelation function trace_summary: Comprehensive summary statistics

trace_mean(trace)[source]

Compute the arithmetic mean of trace data.

Parameters:

trace (numpy.ndarray | list) – Array of trace values.

Returns:

Mean value of the trace.

Return type:

float

Example

>>> trace_mean([1.0, 2.0, 3.0, 4.0, 5.0])
3.0
trace_var(trace)[source]

Compute the variance of trace data.

Uses population variance (ddof=0) for consistency with Kotlin.

Parameters:

trace (numpy.ndarray | list) – Array of trace values.

Returns:

Variance of the trace.

Return type:

float

Example

>>> trace_var([1.0, 2.0, 3.0, 4.0, 5.0])
2.0
trace_scv(trace)[source]

Compute the squared coefficient of variation (SCV).

SCV = Var(X) / E[X]^2

Parameters:

trace (numpy.ndarray | list) – Array of trace values.

Returns:

Squared coefficient of variation.

Return type:

float

Example

>>> trace_scv([1.0, 2.0, 3.0])  # Var=0.667, Mean=2, SCV=0.167
trace_acf(trace, lags=None)[source]

Compute the autocorrelation function at specified lags.

Parameters:
Returns:

Array of autocorrelation values at each lag.

Return type:

numpy.ndarray

Example

>>> trace = np.random.randn(100)
>>> acf = trace_acf(trace, [1, 2, 3])
trace_gamma(trace, limit=1000)[source]

Estimate the autocorrelation decay rate of a trace.

Parameters:
Returns:

Array containing [GAMMA, RHO0, RESIDUALS].

Return type:

numpy.ndarray

Example

>>> gamma, rho0, residuals = trace_gamma(trace_data)
trace_iat2counts(trace, scale)[source]

Compute the counting process from inter-arrival times.

Parameters:
Returns:

Array of counts after scale units of time from each arrival.

Return type:

numpy.ndarray

Example

>>> iat = [0.5, 0.3, 0.8, 0.2, 0.4]
>>> counts = trace_iat2counts(iat, 1.0)
trace_idi(trace, kset, option=None, n=1)[source]

Compute the Index of Dispersion for Intervals.

Parameters:

  • trace (numpy.ndarray | list) – Array of trace values.

  • kset (numpy.ndarray | list) – Set of k values to compute IDI for.

  • option (str) – Aggregation option (None, ‘aggregate’, ‘aggregate-mix’).

  • n (int) – Aggregation parameter.

Returns:

Tuple of (IDI values, support values).

Return type:

Tuple[numpy.ndarray, numpy.ndarray]

Example

>>> idi, support = trace_idi(trace_data, [10, 20, 50])
trace_idc(trace)[source]

Compute the Index of Dispersion for Counts.

Asymptotically equal to IDI.

Parameters:

trace (numpy.ndarray | list) – Array of trace values.

Returns:

IDC value.

Return type:

float

Example

>>> idc = trace_idc(inter_arrival_times)
trace_pmf(X)[source]

Compute the probability mass function of discrete data.

Parameters:

X (numpy.ndarray | list) – Array of discrete values.

Returns:

Tuple of (PMF values, unique values).

Return type:

Tuple[numpy.ndarray, numpy.ndarray]

Example

>>> pmf, values = trace_pmf([1, 2, 2, 3, 3, 3])
trace_shuffle(trace)[source]

Shuffle trace data randomly.

Parameters:

trace (numpy.ndarray | list) – Array of trace values.

Returns:

Shuffled trace array.

Return type:

numpy.ndarray

Example

>>> shuffled = trace_shuffle([1, 2, 3, 4, 5])
trace_joint(trace, lag, order)[source]

Compute joint moments E[X^{k_1}_{i} * X^{k_2}_{i+j} * …].

Parameters:
Returns:

Joint moment value.

Return type:

float

Example

>>> jm = trace_joint(trace, [0, 1], [1, 1])  # E[X_i * X_{i+1}]
trace_iat2bins(trace, scale)[source]

Compute counts in bins with specified timescale.

Parameters:
Returns:

Tuple of (counts per bin, bin membership for each element).

Return type:

Tuple[numpy.ndarray, numpy.ndarray]

Example

>>> counts, bins = trace_iat2bins(iat_data, 1.0)
trace_summary(trace)[source]

Compute comprehensive summary statistics for a trace.

Parameters:

trace (numpy.ndarray | list) – Array of trace values.

Returns:

Array containing [MEAN, SCV, MAD, SKEW, KURT, Q25, Q50, Q75, P95, MIN, MAX, IQR, ACF1, ACF2, ACF3, ACF4, IDC_SCV_RATIO].

Return type:

numpy.ndarray

Example

>>> summary = trace_summary(trace_data)
>>> print(f"Mean: {summary[0]}, SCV: {summary[1]}")
trace_bicov(trace, grid)[source]

Compute bicovariance of a trace.

Bicovariance measures third-order correlation structure at multiple lag combinations.

Parameters:
Returns:

  • bicov: Array of bicovariance values

  • bicov_lags: 2D array of lag combinations (each row is [1, i, j])

Return type:

Tuple of (bicov, bicov_lags) where

Example

>>> trace = np.random.randn(1000)
>>> bicov, lags = trace_bicov(trace, [1, 2, 3, 4, 5])
mtrace_mean(trace, ntypes, types)[source]

Compute the mean of a trace, divided by types.

Parameters:
Returns:

Array containing the mean values for each type.

Return type:

numpy.ndarray

Example

>>> trace = [1.0, 2.0, 3.0, 4.0]
>>> types = [0, 1, 0, 1]
>>> means = mtrace_mean(trace, 2, types)
mtrace_var(trace, ntypes, types)[source]

Compute the variance of a trace, divided by types.

Parameters:
Returns:

Array containing the variance values for each type.

Return type:

numpy.ndarray

Example

>>> var_by_type = mtrace_var(trace, 2, types)
mtrace_count(trace, ntypes, types)[source]

Count elements per type in a trace.

Parameters:
Returns:

Array containing counts for each type.

Return type:

numpy.ndarray

Example

>>> counts = mtrace_count(trace, 2, types)
mtrace_sigma(T, L)[source]

Compute one-step class transition probabilities from a marked trace.

Computes P(C_k = j | C_{k-1} = i) empirically from the trace.

Parameters:
Returns:

C x C matrix where element (i,j) is the probability of observing class j after class i.

Return type:

numpy.ndarray

mtrace_sigma2(T, L)[source]

Compute two-step class transition probabilities from a marked trace.

Computes P(C_k = h | C_{k-1} = j, C_{k-2} = i) empirically.

Parameters:
Returns:

C x C x C 3D array of transition probabilities.

Return type:

numpy.ndarray

mtrace_cross_moment(T, L, k)[source]

Compute the k-th order moment of inter-arrival times between class pairs.

Parameters:
Returns:

C x C matrix where element (i,j) is E[T^k | C_{t-1} = i, C_t = j]

Return type:

numpy.ndarray

mtrace_forward_moment(T, A, orders, norm=True)[source]

Compute forward moments of a marked trace.

Forward moment for class c is E[T_{k+1}^order | C_k = c].

Parameters:
Returns:

Matrix of shape (C, len(orders)) containing forward moments.

Return type:

numpy.ndarray

mtrace_backward_moment(T, A, orders, norm=True)[source]

Compute backward moments of a marked trace.

Backward moment for class c is E[T_k^order | C_k = c].

Parameters:
Returns:

Matrix of shape (C, len(orders)) containing backward moments.

Return type:

numpy.ndarray

mtrace_cov(T, A)[source]

Compute lag-1 covariance between classes in a marked trace.

Parameters:
Returns:

C x C matrix of 2x2 covariance matrices.

Return type:

numpy.ndarray

mtrace_pc(T, L)[source]

Compute the probability of arrival for each class.

Parameters:
Returns:

Array of class probabilities.

Return type:

numpy.ndarray

mtrace_summary(T, L)[source]

Compute comprehensive summary statistics for a marked trace.

Parameters:
Returns:

  • M: First 5 aggregate moments

  • ACF: Autocorrelation for lags 1-100

  • F1, F2: Forward moments of order 1 and 2

  • B1, B2: Backward moments of order 1 and 2

  • C1, C2: Cross moments of order 1 and 2

  • Pc: Class probabilities

  • Pab: One-step transition probabilities

Return type:

Dictionary containing

mtrace_split(T, L)[source]

Split a multi-class trace into per-class traces.

For each class, computes inter-arrival times between consecutive events of that class.

Parameters:
Returns:

List of arrays, one per class, containing inter-arrival times.

Return type:

list

mtrace_merge(t1, t2)[source]

Merge two traces into a single marked (multi-class) trace.

Parameters:
Returns:

  • T: Merged inter-arrival times

  • L: Class labels (1 for t1, 2 for t2)

Return type:

Tuple of (T, L) where

mtrace_joint(T, A, i)[source]

Compute class-dependent joint moments.

Computes E[(X^(a)_j)^i[0] * (X^(a)_{j+1})^i[1]] for all classes a.

Parameters:
Returns:

Array of joint moments, one per class.

Return type:

numpy.ndarray

mtrace_moment(T, A, orders, after=False, norm=False)[source]

Compute class-dependent moments of a multi-class trace.

Parameters:

  • T (numpy.ndarray | list) – Array of inter-arrival times

  • A (numpy.ndarray | list) – Array of class labels

  • orders (numpy.ndarray | list) – Array of moment orders to compute

  • after (bool) – If True, compute moments of Bucholz variables (forward) If False, compute moments of Horvath variables (backward)

  • norm (bool) – If True, normalize by class probability

Returns:

Matrix of shape (C, len(orders)) containing moments per class.

Return type:

numpy.ndarray

mtrace_moment_simple(T, A, k)[source]

Simple interface to compute k-th order moments per class.

Parameters:
Returns:

Array of k-th order moments, one per class.

Return type:

numpy.ndarray

mtrace_bootstrap(T, A, n_samples=100, seed=None)[source]

Generate bootstrap samples from a marked trace.

Parameters:
Returns:

Tuple of (T_boot, A_boot) containing bootstrap samples.

Return type:

Tuple[numpy.ndarray, numpy.ndarray]

mtrace_iat2counts(T, L, scale)[source]

Compute counting process from marked inter-arrival times.

For each class, counts arrivals in windows of specified scale.

Parameters:
Returns:

Array of counts per class per window.

Return type:

numpy.ndarray

Reinforcement Learning (line_solver.api.rl)

Reinforcement Learning (RL) agents for queueing control.

Native Python implementations of RL agents for queueing network control and optimization.

Key classes:

RLTDAgent: Temporal Difference learning agent RLEnvironment: Queueing network environment

Port from:

matlab/src/api/rl/rl_td_agent.m matlab/src/api/rl/rl_env.m

class RLEnvironment(model, queue_indices, source_indices, state_size, gamma=0.99)[source]

Bases: object

Reinforcement Learning environment for queueing networks.

This class wraps a queueing network model and provides an interface for RL agents to interact with it through sampling and state updates.

model

The queueing network model

gamma

Discount factor for future rewards

queue_indices

Indices of queue nodes in model.nodes

source_indices

Indices of source nodes in model.nodes

state_size

Maximum state size to consider

action_size

Number of possible actions (number of queues)

Initialize the RL environment.

Parameters:

  • model (Any) – Queueing network model

  • queue_indices (List[int]) – Indices of queue nodes in model.nodes

  • source_indices (List[int]) – Indices of source nodes in model.nodes

  • state_size (int) – Maximum number of jobs per queue to consider

  • gamma (float) – Discount factor (0 < gamma <= 1)

__init__(model, queue_indices, source_indices, state_size, gamma=0.99)[source]

Initialize the RL environment.

Parameters:

  • model (Any) – Queueing network model

  • queue_indices (List[int]) – Indices of queue nodes in model.nodes

  • source_indices (List[int]) – Indices of source nodes in model.nodes

  • state_size (int) – Maximum number of jobs per queue to consider

  • gamma (float) – Discount factor (0 < gamma <= 1)

is_in_state_space(state)[source]

Check if current state is within the defined state space.

Parameters:

state (ndarray) – Array of queue lengths

Returns:

True if all queue lengths are within state_size bounds

Return type:

bool

is_in_action_space(state)[source]

Check if actions are valid from current state.

Parameters:

state (ndarray) – Array of queue lengths

Returns:

True if all queues can accept new jobs (not at capacity)

Return type:

bool

sample()[source]

Sample the next event from the environment.

Uses the SSA solver to sample a single system event.

Returns:

Tuple of (time_delta, departure_node_index)

Return type:

Tuple[float, int]

update(new_state)[source]

Update the model state after an event.

Parameters:

new_state (ndarray) – New queue lengths for each queue

reset()[source]

Reset the environment to initial state.

Returns:

Initial state (zeros)

Return type:

ndarray

class RLTDAgent(learning_rate=0.05, epsilon=1.0, epsilon_decay=0.99)[source]

Bases: object

Temporal Difference (TD) learning agent for queueing control.

Implements average-reward TD learning for optimal routing decisions in queueing networks. The agent learns a value function V(s) that estimates the long-run average cost from each state.

learning_rate

Step size for value function updates

epsilon

Exploration rate for epsilon-greedy policy

epsilon_decay

Decay factor for exploration rate

V

Value function array

Q

Q-function array (state-action values)

Initialize the TD agent.

Parameters:

  • learning_rate (float) – Learning rate (step size) for updates

  • epsilon (float) – Initial exploration rate (0 to 1)

  • epsilon_decay (float) – Decay factor applied to epsilon each episode

__init__(learning_rate=0.05, epsilon=1.0, epsilon_decay=0.99)[source]

Initialize the TD agent.

Parameters:

  • learning_rate (float) – Learning rate (step size) for updates

  • epsilon (float) – Initial exploration rate (0 to 1)

  • epsilon_decay (float) – Decay factor applied to epsilon each episode

reset(env)[source]

Reset agent and environment.

Parameters:

env (RLEnvironment) – The RL environment

get_value_function()[source]

Get the learned value function.

get_q_function()[source]

Get the learned Q-function.

solve(env, num_episodes=10000, verbose=True)[source]

Train the agent using TD learning.

Runs TD(0) learning for the specified number of episodes, learning an optimal routing policy for the queueing network.

Parameters:

  • env (RLEnvironment) – The RL environment

  • num_episodes (int) – Number of training episodes

  • verbose (bool) – Whether to print progress

Returns:

  • ‘V’: Learned value function

  • ’Q’: Learned Q-function

  • ’mean_cost_rate’: Final estimated average cost rate

Return type:

Dictionary with training results

class RLEnvironmentGeneral(model, queue_indices, action_node_indices, state_size, gamma=0.99)[source]

Bases: object

General RL environment for queueing networks with flexible action spaces.

Unlike RLEnvironment which assumes sources dispatch to queues, this class supports arbitrary action nodes with configurable routing destinations. Actions are dispatching decisions at specific nodes that route jobs to connected downstream nodes.

model

The queueing network model

gamma

Discount factor for future rewards

queue_indices

Indices of queue nodes in model.nodes

nqueues

Number of queues

action_node_indices

Indices of nodes where routing actions are needed

state_size

Maximum state size to consider

action_space

Dict mapping action_node -> list of possible destination nodes

MATLAB: matlab/src/api/rl/rl_env_general.m

Initialize the general RL environment.

Parameters:

  • model (Any) – Queueing network model

  • queue_indices (List[int]) – Indices of queue nodes in model.nodes

  • action_node_indices (List[int]) – Indices of nodes where dispatch actions are taken

  • state_size (int) – Maximum number of jobs per queue to consider

  • gamma (float) – Discount factor (0 < gamma <= 1)

__init__(model, queue_indices, action_node_indices, state_size, gamma=0.99)[source]

Initialize the general RL environment.

Parameters:

  • model (Any) – Queueing network model

  • queue_indices (List[int]) – Indices of queue nodes in model.nodes

  • action_node_indices (List[int]) – Indices of nodes where dispatch actions are taken

  • state_size (int) – Maximum number of jobs per queue to consider

  • gamma (float) – Discount factor (0 < gamma <= 1)

is_in_state_space(state)[source]

Check if the given state is within the defined state space.

Parameters:

state (ndarray) – Array of queue lengths (one per queue)

Returns:

True if all queue lengths are within state_size bounds

Return type:

bool

is_in_action_space(state)[source]

Check if actions are valid from the given state.

Parameters:

state (ndarray) – Array of queue lengths

Returns:

True if at least one queue can accept a new job

Return type:

bool

sample()[source]

Sample the next event from the environment.

Uses the SSA solver to sample a single system event.

Returns:

Tuple of (time_delta, departure_node, arrival_node, sample_data)

Return type:

Tuple[float, int, int, Any]

update(sample_data)[source]

Update the model state using the sample event data.

Applies the state transitions from the sampled events to the model.

Parameters:

sample_data (Any) – Sample data from the SSA solver

reset()[source]

Reset the environment to initial state.

Returns:

Initial state (zeros)

Return type:

ndarray

class RLTDAgentGeneral(learning_rate=0.1, epsilon=1.0, epsilon_decay=0.9999)[source]

Bases: object

General TD learning agent for queueing control with flexible policies.

Supports multiple solve methods: - solve(): TD control with tabular value function - solve_for_fixed_policy(): TD learning (evaluation) with fixed heuristic - solve_by_hashmap(): TD control with hash-map value function - solve_by_linear(): TD control with linear function approximation - solve_by_quad(): TD control with quadratic function approximation

Works with RLEnvironmentGeneral for arbitrary action spaces.

MATLAB: matlab/src/api/rl/rl_td_agent_general.m

Initialize the general TD agent.

Parameters:

  • learning_rate (float) – Learning rate for updates

  • epsilon (float) – Initial exploration rate (0 to 1)

  • epsilon_decay (float) – Decay factor applied to epsilon each episode

__init__(learning_rate=0.1, epsilon=1.0, epsilon_decay=0.9999)[source]

Initialize the general TD agent.

Parameters:

  • learning_rate (float) – Learning rate for updates

  • epsilon (float) – Initial exploration rate (0 to 1)

  • epsilon_decay (float) – Decay factor applied to epsilon each episode

reset(env)[source]

Reset agent and environment.

Parameters:

env (RLEnvironmentGeneral) – The RL environment

get_value_function()[source]

Get the learned value function.

solve_for_fixed_policy(env, num_episodes=10000, verbose=True)[source]

TD learning for value function evaluation with heuristic routing.

Evaluates the value function under the existing (model-defined) routing policy without modifying routing decisions.

Parameters:

  • env (RLEnvironmentGeneral) – The general RL environment

  • num_episodes (int) – Number of training episodes

  • verbose (bool) – Whether to print progress

Returns:

Learned value function array

Return type:

ndarray

MATLAB: rl_td_agent_general.solve_for_fixed_policy

solve(env, num_episodes=10000, verbose=True)[source]

TD control with tabular value function.

Learns an optimal routing policy using epsilon-greedy exploration. At each action node departure, the agent selects the best downstream destination based on the current value function.

Parameters:

  • env (RLEnvironmentGeneral) – The general RL environment

  • num_episodes (int) – Number of training episodes

  • verbose (bool) – Whether to print progress

Returns:

Learned value function array

Return type:

ndarray

MATLAB: rl_td_agent_general.solve

solve_by_hashmap(env, num_episodes=10000, verbose=True)[source]

TD control using hash-map value function.

Uses a dictionary to store value estimates only for visited states, with an ‘external’ default for unvisited states.

Parameters:

  • env (RLEnvironmentGeneral) – The general RL environment

  • num_episodes (int) – Number of training episodes

  • verbose (bool) – Whether to print progress

Returns:

X: State features matrix (n_states x (1 + nqueues)) Y: Value estimates (n_states x 1)

Return type:

Tuple of (X, Y) where

MATLAB: rl_td_agent_general.solve_by_hashmap

solve_by_linear(env, num_episodes=10000, verbose=True)[source]

TD control with linear function approximation.

Learns a linear value function: v(q1,…,qn) = w0 + w1*q1 + … + wn*qn

Parameters:

  • env (RLEnvironmentGeneral) – The general RL environment

  • num_episodes (int) – Number of training episodes

  • verbose (bool) – Whether to print progress

Returns:

X: State features matrix Y: Value estimates coefficients: Linear regression coefficients

Return type:

Tuple of (X, Y, coefficients) where

MATLAB: rl_td_agent_general.solve_by_linear

solve_by_quad(env, num_episodes=10000, verbose=True)[source]

TD control with quadratic function approximation.

Learns a quadratic value function: v(q1,…,qn) = sum_{i,j} w_{ij} * q_i * q_j

Parameters:

  • env (RLEnvironmentGeneral) – The general RL environment

  • num_episodes (int) – Number of training episodes

  • verbose (bool) – Whether to print progress

Returns:

X_quad: Augmented features (linear + quadratic terms) Y: Value estimates coefficients: Quadratic regression coefficients

Return type:

Tuple of (X_quad, Y, coefficients) where

MATLAB: rl_td_agent_general.solve_by_quad

Utilities Module (line_solver.utils)

Mock module that can handle any attribute access

Constants Module (line_solver.constants)

Constants and enumerations for LINE queueing network models.

This module defines the various constants, enumerations, and strategies used throughout LINE for specifying model behavior, including:

  • Scheduling strategies (FCFS, LCFS, PS, etc.)

  • Routing strategies (PROB, RAND, etc.)

  • Node types (SOURCE, QUEUE, SINK, etc.)

  • Job class types (OPEN, CLOSED)

  • Solver types and options

  • Activity precedence types for layered networks

  • Call types and drop strategies

These constants ensure type safety and consistency across the API.

class ActivityPrecedenceType(*values)[source]

Bases: Enum

Types of activity precedence relationships in layered networks.

These specify how activities are ordered and synchronized: - PRE_SEQ: Sequential prerequisite (must complete before) - PRE_AND: AND prerequisite (all must complete before) - PRE_OR: OR prerequisite (any must complete before) - POST_SEQ: Sequential post-condition - POST_AND: AND post-condition - POST_OR: OR post-condition - POST_LOOP: Loop post-condition - POST_CACHE: Cache post-condition

class CallType(*values)[source]

Bases: Enum

Types of calls between tasks in layered networks.

  • SYNC: Synchronous call (caller waits for response)

  • ASYNC: Asynchronous call (caller continues immediately)

  • FWD: Forward call (caller terminates, response goes to caller’s caller)

class DropStrategy(*values)[source]

Bases: Enum

Strategies for handling queue overflow and capacity limits.

  • WaitingQueue: Jobs wait in a waiting queue when capacity is exceeded

  • Drop: Jobs are dropped (lost) when capacity is exceeded

  • BlockingAfterService: Jobs are blocked after service completion

class SignalType(*values)[source]

Bases: Enum

Types of signals for signal classes in G-networks and related models.

NEGATIVE

Removes a job from the destination queue (G-network negative customer)

REPLY

Triggers a reply action

CATASTROPHE

Removes ALL jobs from the destination queue

NEGATIVE = 1
REPLY = 2
CATASTROPHE = 3
class RemovalPolicy(*values)[source]

Bases: Enum

Removal policies for negative signals in G-networks.

RANDOM

Select job uniformly at random from all jobs at the station

FCFS

Remove the oldest job (first arrived)

LCFS

Remove the newest job (last arrived)

RANDOM = 1
class EventType(*values)[source]

Bases: Enum

Types of events in discrete-event simulation.

  • INIT: Initialization event

  • LOCAL: Local processing event

  • ARV: Job arrival event

  • DEP: Job departure event

  • PHASE: Phase transition event in multi-phase processes

  • READ: Cache read event

  • STAGE: Staging area event

READ = 6
class JobClassType(*values)[source]

Bases: Enum

Types of job classes in queueing networks.

  • OPEN: Open class (jobs arrive from outside the system)

  • CLOSED: Closed class (fixed population circulating in the system)

  • DISABLED: Disabled class (not currently active)

class JoinStrategy(*values)[source]

Bases: Enum

Strategies for join node synchronization in fork-join networks.

  • STD: Standard join (wait for all parallel branches)

  • PARTIAL: Partial join (proceed when some branches complete)

  • Quorum: Quorum-based join (wait for minimum number of branches)

  • Guard: Guard condition join (custom completion criteria)

class MetricType(*values)[source]

Bases: Enum

Types of performance metrics that can be computed.

SysDropR = 11
SysPower = 13
Tard = 24
SysTard = 25
class Metric(metric_type, job_class, station=None)[source]

Bases: object

An output metric of a Solver, such as a performance index.

class TranResult(t, metric)[source]

Bases: object

Container for transient result time series with attribute access.

class NodeType(*values)[source]

Bases: Enum

Types of nodes in queueing network models.

static from_line(obj)[source]
class ProcessType(*values)[source]

Bases: Enum

Types of stochastic processes for arrivals and service times.

static fromString(obj)[source]
class RoutingStrategy(*values)[source]

Bases: Enum

Strategies for routing jobs between network nodes. Values must match MATLAB’s RoutingStrategy constants for JMT compatibility.

RL = 7
class SchedStrategy(*values)[source]

Bases: Enum

Scheduling strategies for service stations.

Values match lang/base.py SchedStrategy for consistency.

LCFSPI = 3
HOL = 9
LPS = 17
SETF = 18
FCFSPR = 22
EDF = 23
JOIN = 25
EDD = 27
SRPT = 28
SRPTPRIO = 29
LCFSPRIO = 30
LCFSPRPRIO = 31
LCFSPIPRIO = 32
FCFSPRPRIO = 33
FCFSPIPRIO = 34
PSJF = 36
FB = 37
LAS = 38
LRPT = 39
static fromString(obj)[source]
static fromLINEString(sched)[source]
static toID(sched)[source]
class SchedStrategyType(*values)[source]

Bases: Enum

Categories of scheduling strategies by preemption behavior.

class ServiceStrategy(*values)[source]

Bases: Enum

Service strategies defining service time dependence.

class SolverType(*values)[source]

Bases: Enum

Types of solvers available in LINE.

DES = 3
class TimingStrategy(*values)[source]

Bases: Enum

Timing strategies for transitions in Petri nets.

class VerboseLevel(*values)[source]

Bases: Enum

Verbosity levels for LINE solver output.

class PollingType(*values)[source]

Bases: Enum

Polling strategies for polling systems.

static fromString(obj)[source]
class HeteroSchedPolicy(*values)[source]

Bases: Enum

Scheduling policies for heterogeneous multiserver queues.

ORDER = 1
ALIS = 2
ALFS = 3
FAIRNESS = 4
FSF = 5
RAIS = 6
static fromString(obj)[source]
class GlobalConstants[source]

Bases: object

Global constants and configuration for the LINE solver.

Immediate = 100000000.0
classmethod getInstance()[source]

Get the singleton instance of GlobalConstants.

classmethod get_instance()

Get the singleton instance of GlobalConstants.

classmethod getVerbose()[source]

Get the current verbosity level.

classmethod get_verbose()

Get the current verbosity level.

classmethod setVerbose(verbosity)[source]

Set the verbosity level for solver output.

classmethod set_verbose(verbosity)

Set the verbosity level for solver output.

classmethod getConstants()[source]

Get a dictionary of all global constants.

classmethod get_constants()

Get a dictionary of all global constants.