Product-Form Queueing Networks

pfqn_qlen_joint_moments(L, N, Z=None, pairs=None, route='auto', lg_source=None, method='ca', options=None)[source]

Joint moments of the queue-length vector of a closed product-form network.

The coordinates are (station, class) pairs. Two pairs sharing a class give the cross-station covariance of that class; two pairs sharing a station give the cross-class covariance at that station, which is what a class-oriented method of moments (pfqn_comomrm and its relatives) is positioned to deliver.

The result is EXACT: the survival array covers the support, since a queue length is bounded by the population of its class.

Parameters:

L – Service demand matrix of the QUEUEING stations (M x R). Delay (infinite-server) stations belong in Z, their marginals following a different law. Load-dependent or multiserver stations are out of scope for both routes and must not be passed here.
N – Population vector (R,).
Z – Think time vector (R,), zeros if omitted.
pairs – Sequence of (station, class) 0-based pairs, one per dimension of the returned arrays. Defaults to every class of every station.
route – ‘tail’ uses the single-class survival identity and only touches the original network; ‘pmf’ uses the complementary-network joint distribution and works for any number of classes; ‘auto’ picks ‘tail’ when R = 1 and ‘pmf’ otherwise.
lg_source – Where log G comes from. None calls pfqn_nc. A callable is invoked ONCE per network with (Lsub, pops), pops being a (P x R) integer array, and must return P values of log G with NaN where it cannot serve; those are filled in by pfqn_nc. An ndarray is read as a precomputed table indexed by population, which is what a convolution sweep produces for free. Note that the ‘pmf’ route queries the COMPLEMENTARY network, so a table must be its table, not the original network’s.
method – Method passed to pfqn_nc for unserved populations.
options – Options passed to pfqn_nc.

Returns:

‘tail’ (survival), ‘binomial’, ‘factorial’, ‘raw’, ‘central’, ‘cumulant’, plus ‘mean’, the covariance matrix ‘cov’, and ‘info’ holding the route, the number of populations requested, how many the source served, how many pfqn_nc evaluations were needed, and the pairs used.

Return type:

Dictionary with the joint arrays over the selected coordinates

Raises:

ValueError – If the arguments are inconsistent, if a pair is out of range, or if the ‘tail’ route is requested with several classes.

Example

out = pfqn_qlen_joint_moments(L, N, Z, pairs=[(0, 0), (0, 1)]) cov01 = out[‘cov’][0, 1]

pfqn_nc_resolved_method(method)[source]

The algorithm pfqn_nc actually runs for METHOD.

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

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

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

PANACEA asymptotic expansion for load-independent closed networks.

McKenna-Mitra normal-usage expansion whose coefficients are linear combinations of pseudonetwork partition functions evaluated by convolution.

Parameters:

L (ndarray) – Service demand matrix (M x R)
N (ndarray) – Population vector (R,)
Z (ndarray) – Think time vector (R,)
terms (int) – Number of terms in the normal-usage asymptotic series (1, 2, or 3; default 3), as selectable in the original PANACEA package (Ramakrishnan-Mitra, BSTJ 61(10):2849-2872, 1982)

Returns:

Tuple (G, lG) - normalizing constant and its log, both NaN when the model is not in normal usage

Return type:

pfqn_panaceald(L, N, Z=None, mu=None, terms=3)[source]

PANACEA asymptotic expansion for load-dependent closed networks.

Mitra-McKenna (JACM 33(3):568-592, 1986) load-dependent PANACEA: the expansion coefficients A_n are linear combinations of partition functions of a pseudonetwork whose load dependence is the phi(n) transform of the original {f(n)}. See _kb/03-api-layer.md (pfqn/ family, pfqn_panaceald).

Parameters:

L (ndarray) – Service demand matrix (M x R)
N (ndarray) – Population vector (R,)
Z (ndarray) – Think time vector (R,) or matrix (D x R), summed over rows
mu (ndarray) – Load-dependent rate matrix (M x sum(N))
terms (int) – Number of terms in the normal-usage asymptotic series (1, 2 or 3)

Returns:

Tuple (G, lG) - normalizing constant and its log, both NaN when the model is not in normal usage

Return type:

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 (ndarray) – Service demand matrix (M x R) where M is stations, R is classes
N (ndarray) – Population vector (1 x R or R,) - number of jobs per class
Z (ndarray) – Think time vector (1 x R or R,) - think time per class (default 0)

Returns:

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

Return type:

Tuple (G, lG, X) where

References

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

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

pfqn_ls(L, N, Z=None, I=100000)[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 (ndarray) – Service demand matrix (M x R)
N (ndarray) – Population vector (R,)
Z (ndarray) – Think time vector (R,) (default: zeros)
I (int) – Number of samples for Monte Carlo integration (default: 100000)

Returns:

G: Estimated normalizing constant lG: log(G)

Return type:

Tuple (G, lG) where

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

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

Markov chain Monte Carlo estimate of the throughputs and queue lengths.

Estimates the class throughputs X(r) = G(N-e_r)/G(N) and the mean queue lengths Q(i,r) of a CLOSED multiclass product-form (BCMP, no type changes) network by the REGULARIZATION algorithm of Chen and O’Cinneide (ACM TOMACS 8(3), 1998).

The three steps of the paper are:

I. CONSTRUCT THE REGULARIZED NETWORK. Write rho(i,r) for the surrogate traffic intensity of class r at station i – here the service demand, since rho = lambda/mu is a visit ratio over a service rate – and rho(r) = sum_i rho(i,r). The regularized network has the same stations, classes and populations, UNIT service rates at every station, the processor-sharing discipline, and a routing matrix that depends on the destination only, P*(i->m | class r) = rho(m,r)/rho(r). By Theorem 2.1 it is a REVERSIBLE chain with the SAME steady-state distribution as the original network, and its throughputs satisfy Theta*(r) = rho(r)*Theta(r).

II. SIMULATE IT at service-completion epochs. With Y(i,r) the number of class-r jobs at station i, Y(i) their total and Psi_i(k) = min(s_i,k) the number of busy servers:

r(i,r) = Y(i,r)/Y(i) * Psi_i(Y(i)),   r(r) = sum_i r(i,r),
r      = sum_i Psi_i(Y(i)),

the next completion is of class r at station i with probability r(i,r)/r, and the conditional expected time to it is 1/r. Equation (10) of the paper is the holding-time weighted ratio estimator Theta*(r) = sum_t r(r,t)/r(t) / sum_t 1/r(t), and the same weights give the time-average queue lengths, which need no transformation at all because the two networks share their steady state.

TRANSFORM BACK: X(r) = Theta*(r)/rho(r).

Because P* forgets the station of origin and every station serves at unit rate, the regularized chain has neither the slowly mixing routing chain nor the customer-trapping slow station that make the original chain converge slowly. The paper proves O(N^2*M^3) mixing in two special cases (Section 4) and reports the general behaviour experimentally (Section 5).

Delay (infinite-server) demand enters as ONE extra station with s = inf and demand Z. Aggregating infinite-server stations that way is exact in the product form, since their joint term is multinomial in the per-class totals.

Confidence: the run is split into non-overlapping batches (Schmeiser 1982, 30 by default, the count used in the tables of the paper), the batch means of the ratio estimator give a standard error, and the intervals are the paper’s two-sigma ones. The estimator is a ratio of correlated averages, so it carries an O(1/samples) bias on top of the initialization bias; the paper ignores both, this implementation additionally discards a warm-up fraction (10% by default).

Parameters:

L ((M, R) array) – Per-class service demands at the M queueing stations.
N ((R,) array) – Closed population vector; finite and integer.
Z ((R,) or (K, R) array, optional) – Aggregated think times, summed over rows; None or zeros if the model has no delay.
s ((M,) array, optional) – Number of servers at each queueing station, inf for an infinite server; None means all stations single-server.
options (dict or options object, optional) – Fields samples (default 1e5), seed, and inside config the batch count mcmc_batches (30) and the warm-up fraction mcmc_burnin (0.1).

Returns:

Throughputs, queue lengths and their two-sigma intervals.

Return type:

PfqnMcmcResult

Examples

>>> mu = np.array([0.2, 0.5, 0.8]); sets = [[0, 1, 2], [0, 1], [0, 2], [1, 2]]
>>> L = np.zeros((3, 4)); Z = np.zeros(4)
>>> for c, st in enumerate(sets):
...     L[st, c] = 1.0 / mu[st]
...     if c > 0:
...         Z[c] = 1.0 / 0.5
>>> res = pfqn_mcmc(L, 3 * np.ones(4), Z, options={'samples': 100000, 'seed': 23000})

class PfqnMcmcResult(X, Q, Xse, Xlo, Xhi, Qse, Qlo, Qhi, batches, samples, burnin)[source]

Bases: NamedTuple

Estimates of pfqn_mcmc() with their batch-means confidence intervals.

There is deliberately no normalizing constant here: the estimator is a ratio of holding-time weighted averages that yields G(N-e_r)/G(N) directly, and G itself never enters the algorithm.

Create new instance of PfqnMcmcResult(X, Q, Xse, Xlo, Xhi, Qse, Qlo, Qhi, batches, samples, burnin)

X: ndarray: (R,) throughput estimates G(N-e_r)/G(N)

Q: ndarray: (M, R) mean queue lengths at the queueing stations

Xse: ndarray: (R,) batch-means standard error of X

Xlo: ndarray: (R,) lower end of the two-sigma interval for X

Xhi: ndarray: (R,) upper end of the two-sigma interval for X

Qse: ndarray: (M, R) batch-means standard error of Q

Qlo: ndarray: (M, R) lower end of the two-sigma interval for Q

Qhi: ndarray: (M, R) upper end of the two-sigma interval for Q

batches: int: batches the run was split into

samples: int: service completions simulated after warm-up

burnin: int: completions discarded as warm-up

pfqn_clw(L, N, Z=None, m=None, l=None, gamma=None, euler=True, euler_n=11, euler_m=20, euler_tol=1e-10, euler_maxm=160, beta=None, dimred=True, dimred_maxd=4)[source]

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

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

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

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

Both of the paper’s accelerations are applied. Dimension reduction by decomposition (Sec. 3, Sec. 5.4) removes from the interdependence graph of the factors the subset D minimizing |D| + max_i |S_i(D)| (eq. 3.3): with the D variables fixed on their contours the remaining factors share no variable, so each connected component is inverted separately and the results multiplied. Euler summation (Sec. 2.4, eq. 2.22) replaces the nearly alternating inner sum of (2.3) by the Euler sum of its first n+m+1 terms, applied once for k >= 0 and once for k < 0, so prod_j K_j becomes prod_j min(n+m+1, K_j) in the cost (eq. 2.26); the order m is doubled until the paper’s own estimate |E(m,n) - E(m,n+1)| falls under euler_tol.

Parameters:

L (ndarray) – (q’ x p) single-server relative traffic intensities, L[i,j]=rho_{ji}.
N (ndarray) – (p,) closed-chain population vector K.
Z (ndarray) – (p,) aggregate infinite-server relative intensities rho_{j0}. Default 0.
m (ndarray) – (q’,) queue multiplicities m_i. Default ones.
l (ndarray) – (p,) inner lattice parameters l_j indexed by chain. Default by inversion depth: 1 at depth 1, 2 at depths 2-3, 3 deeper.
gamma (ndarray) – (p,) aliasing parameters gamma_j. Default by depth: 11, 13, 13, 15.
euler (bool) – apply Euler summation where K_j > euler_n + euler_m.
euler_n (int) – terms summed exactly before averaging (n in eq. 2.22).
euler_m (int) – starting order of the Euler averaging (m in eq. 2.22).
euler_tol (float) – relative tolerance on |E(m,n) - E(m,n+1)|.
euler_maxm (int) – largest Euler order reached by doubling.
beta (ndarray) – (p,) multipliers on the scale parameters alpha_j, the manual tuning of page 956 (the paper uses 0.8 <= beta <= 1.2 on its largest examples). Default ones.
dimred (bool) – apply dimension reduction by decomposition.
dimred_maxd (int) – largest |D| examined when minimizing (3.3).

Returns:

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

Return type:

Tuple (G, lG)

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

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

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

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

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

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

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

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

Parameters:

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

Returns:

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

Return type:

Tuple (G, lG)

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

pfqn_perm(A, m=None)[source]

Permanent of a demand matrix, with optional column multiplicities.

The pfqn_ entry point of the permanent library. It exists because the product-form joint queue-length probability of the per-station TOTAL populations is a permanent of the demand matrix replicated once per job, which is a normalizing-constant quantity rather than a general-purpose linear algebra one; see pfqn_jointmarg.

Orientation is chosen before repeated lines are grouped. That is a correctness concern, not an optimisation: perm(A) is transpose-invariant but the Ryser sum is not, and exploiting repeated rows silently expands the transpose.

Parameters:

A (ndarray) – Square matrix, or the matrix of distinct columns when m is given
m (ndarray) – Multiplicity of each column of A; sum(m) must be the order of the expanded matrix

Returns:

The permanent value; 1.0 for the empty matrix

Return type:

References

H. J. Ryser, “Combinatorial Mathematics”, Carus Mathematical Monographs 14, Mathematical Association of America, 1963.

pfqn_qdamva(L, N, Z=None, mu=None, Q0=None, tol=1e-6, maxiter=10000)[source]

QD-AMVA on a closed multiclass product-form network.

Parameters:

L (ndarray) – (M x R) service demand matrix.
N (ndarray) – (R,) population vector, finite.
Z (ndarray | None) – (R,) think time vector; None means no think time.
mu (ndarray | None) – (M x smax) queue-dependent rate multipliers; None means none.
Q0 (ndarray | None) – (M x R) initial guess; None means the reference’s demand split.
tol (float) – convergence tolerance on the queue lengths.
maxiter (int) – maximum number of iterations.

Returns:

(M x R) mean queue lengths. X: (R,) per-class throughputs. U: (M x R) per-class utilizations, carrying the g scaling. iter: number of iterations performed. R: (M x R) per-class residence times, Q = X * R.

Return type:

Q

pfqn_qsa(L, N, Z=None, type=None, tol=1e-10, maxiter=100, levels=3, QN0=None)[source]

Queue-Shift Approximation for a closed product-form network.

Parameters:

L – Service demand matrix (M x R)
N – Population vector (R,)
Z – Think time vector (default 0)
type – Scheduling strategy per station; SchedStrategy.INF marks a delay centre, whose demand enters the cycle time without a queueing term (the paper’s DC set)
tol (float) – Residual tolerance of the Newton iteration (default 1e-10)
maxiter (int) – Maximum Newton iterations (default 100)
levels (int) – 2 for the two-level QSA of eq. (14), 3 for eq. (16)
QN0 – Warm start for the Bard-Schweitzer initialization (M x R)

Returns:

Tuple of (Q, U, W, C, X, totiter); W holds residence times and C cycle times, matching pfqn_linearizer.

Return type:

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

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

Bard Large Customer Population (LCP) approximate MVA.

Y. Bard, “Some extensions to multiclass queueing network analysis”, in Performance of Computer Systems, North-Holland, 1979. The arrival-instant queue length is estimated by the time-averaged one WITHOUT removing the arriving customer,

A_k^(c)(N) = Q_k(N - 1_c) ~= Q_k(N) = sum_s Q_ks(N),

since with a large population one customer less cannot change the mean queue lengths appreciably. Setting the Bard-Schweitzer proportional term Q_kc(N)/N_c to zero recovers this algorithm, so LCP is uniformly more pessimistic than pfqn_bs and is inaccurate at small populations.

Returns (XN, QN, UN, RN, it), matching pfqn_bs.

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

Chow Second Approximation (SA) approximate MVA.

W.-M. Chow, “Approximations for large scale closed queueing networks”, Perform. Eval. 3(1), 1983. The arrival-instant queue length is written exactly as

A_k^(c)(N) = Q_k(N - 1_c) = Q_k(N) (1 + theta_ck), theta_ck = [Q_k(N - 1_c) - Q_k(N)] / Q_k(N),

and the theta-terms are estimated ONCE, off the Bard LCP solution, before the fixed point is run. variant='forward' uses Qhat(N + 1_c) and 'backward' uses Qhat(N - 1_c); Chow reports the forward form to be the more accurate, so it is the default. Setting every theta to zero recovers pfqn_lcp.

Returns (XN, QN, UN, RN, it).

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

Eager Looping approximate MVA bounds.

D. L. Eager, “Bounding Algorithms for Queueing Network Models of Computer Systems”, Ph.D. thesis, Tech. Rept. CSRG-156, University of Toronto, 1984. Looping supplies the initial pessimistic and optimistic estimates that the multiple-class performance bound hierarchy starts from, so it carries a pair of bounds rather than a single fixed point. A HEAP H_j is the class-j congestion the current queue-length LOWER BOUNDS have not accounted for; it is charged back at the pessimistic inflation factor V_c = max_k D_ck or the optimistic one L_c = min_k D_ck. The whole bracket rests on Q_jk(N - 1_c) being a lower bound, which is why the queue lengths are seeded from Little’s law at the station,

Q_jk(N - 1_c) = X_j(N - 1_c) R_jk(N - 1_c) >= n_j D_jk / (Z_j + U_j),

with R_jk >= D_jk and U_j any UPPER bound on R_j(N - 1_c): the level-0 PBH bound B_j, or the pessimistic R_j(N) of the current iterate, whichever is smaller, since response time is nondecreasing in the population. The refinement is then monotone and every iterate is a bound. It previously used the convolution identity of Zahorjan (1980) with one class-level ratio at every station, which is not a per-station under-estimate; the queue lengths stopped being lower bounds, the heaps clamped to zero, and R^(opt) collapsed onto R^(pess). The level-0 multiple-class PBH bounds on the mean response time are

J_j(n) = sum_k D_jk, B_j(n) = sum_k D_jk + (sum(n) - 1) max_k D_jk,

i.e. an arriving customer queues behind nobody, respectively behind every other customer in the network at its own worst centre.

Returns (Xlo, Xup, QN, RN, it) with Xlo the pessimistic and Xup the optimistic throughput bound.

pfqn_pam(L, N, Z=None, variant='pamb')[source]

Hsieh-Lam Proportional Approximation Methods (PAMB/PAMI/PAMT).

C. T. Hsieh, S. S. Lam, “PAM - A noniterative approximate solution method for closed multichain queueing networks”, ACM SIGMETRICS Perform. Eval. Rev. 16(1), 1988. The three variants are NONITERATIVE: the queue lengths are seeded by the proportion of a class demand that falls at each centre,

E_ck = D_ck / sum_i D_ci, Q_ck(N) = E_ck N_c,

and the MVA equations are then unrolled a fixed number of times. ‘pamb’ applies the last MVA step; ‘pami’ additionally scales a class down wherever it would drive a centre past full utilization; ‘pamt’ seeds at N - 1_i - 1_j and applies the last TWO MVA steps before that capping. The seed spreads the whole class population over the queueing centres and ignores Z, exactly as published: PAM buys speed, not accuracy.

Returns (XN, QN, UN, RN).

pfqn_clust(L, N, Z=None, subnets=None, localclasses=None, inner='lin', tol=1e-6, maxiter=1000)[source]

de Souza e Silva-Lavenberg-Muntz Clustering Approximation (CA).

E. de Souza e Silva, S. S. Lavenberg, R. R. Muntz, “A clustering approximation technique for queueing network models with a large number of chains”, IEEE Trans. Computers C-35(5), 1986. The network is covered by subnetworks whose union is the whole network but which need not be disjoint. Every class visiting a subnetwork S is either LOCAL to S, and is then solved inside it, or FOREIGN, and is then seen only through the utilization it leaves behind. Each subnetwork is solved by an ordinary approximate MVA algorithm with two replacements: the complement of S is collapsed into a per-class delay P_c and the foreign classes into a per-centre utilization U_k,

X_c(N) = N_c / (sum_{k in S} R_ck(N) + Z_c + P_c), Q_k(N) = [sum_{c in LC(S)} R_ck(N) X_c(N) + U_k] / (1 - U_k).

Choosing the PE algorithm for every subnetwork reproduces global PE exactly, so the useful setting is Linearizer inside, PE outside.

When no decomposition is supplied the criterion of the paper is applied automatically: the cheap PAMB estimate of the centre utilizations is taken, every class is attached to the centre where it loads the most, classes sharing that centre form one cluster, and the subnetwork of a cluster is the set of centres its classes visit.

Returns (XN, QN, UN, RN, it).

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

de Souza e Silva-Muntz Improved Linearizer (IL).

E. de Souza e Silva, R. R. Muntz, “A note on the computational cost of the Linearizer algorithm for queueing networks”, IEEE Trans. Computers 39(6), 1990. Linearizer evaluates the arrival-instant queue length as

A_k^(c)(n) = sum_i (n_i - delta_c^(i)) [Q_ik(n)/n_i + Delta^(i)_ck],

re-summing the C Delta-terms at every Core iteration, at every one of the C+1 populations: O(K C^3) per refresh pass. IL splits that sum into the part that moves with the Core iterate and the part that does not,

A_k^(c)(n) = sum_i (n_i - delta_c^(i)) Q_ik(n)/n_i + xi_ck(n), xi_ck(N) = sum_i (N_i - delta_c^(i)) Delta^(i)_ck, xi_ck(N - 1_j) = xi_ck(N) - Delta^(j)_ck,

so the C K aggregates xi are computed ONCE per refresh pass and each Core iteration then costs O(K C) instead of O(K C^2). Because the split is an identity and not an approximation, the fixed point is the one Linearizer reaches: pfqn_dmlin and pfqn_linearizer agree to round-off.

type_sched is accepted for signature parity and unused: the Linearizer family in LINE treats every station as single-server PS.

Returns (Q, U, W, T, C, X, totiter), matching pfqn_linearizer.

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

Linearizer approximate MVA algorithm.

Parameters:

L (ndarray) – Demand matrix (M x R) - rows are stations, columns are classes
N (ndarray) – Population vector (R,)
Z (ndarray) – Think time vector (R,)
sched_type (List[str]) – List of scheduling strategies per station (‘FCFS’, ‘PS’, etc.)
tol (float) – Convergence tolerance; ‘cn’ or NaN selects the published Linearizer termination test of Chandy and Neuse, Commun. ACM 25(2), 1982, p.129: each Core call stops when max_{i,r}|dQ(i,r)|/N_r falls below pfqn_cntol evaluated at the population Core is running at, rather than on the Frobenius norm of dQ. See pfqn_cntol.
maxiter (int) – Maximum iterations

Returns:

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

Return type:

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

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

General-form linearizer approximate MVA.

Parameters:

L (ndarray) – Demand matrix (M x R)
N (ndarray) – Population vector (R,)
Z (ndarray) – Think time vector (R,)
sched_type (List[str]) – List of scheduling strategies per station
tol (float) – Convergence tolerance; ‘cn’ or NaN selects the published Linearizer termination test of Chandy and Neuse, Commun. ACM 25(2), 1982, p.129: each Core call stops when max_{i,r}|dQ(i,r)|/N_r falls below pfqn_cntol evaluated at the population Core is running at, rather than on the Frobenius norm of dQ. See pfqn_cntol.
maxiter (int) – Maximum iterations
alpha (float) – Linearization parameter (scalar)

Returns:

Same as pfqn_linearizer

Return type:

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

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

Extended general-form linearizer with class-specific parameters.

Parameters:

L (ndarray) – Demand matrix (M x R)
N (ndarray) – Population vector (R,)
Z (ndarray) – Think time vector (R,)
sched_type (List[str]) – List of scheduling strategies per station
tol (float) – Convergence tolerance; ‘cn’ or NaN selects the published Linearizer termination test of Chandy and Neuse, Commun. ACM 25(2), 1982, p.129: each Core call stops when max_{i,r}|dQ(i,r)|/N_r falls below pfqn_cntol evaluated at the population Core is running at, rather than on the Frobenius norm of dQ. See pfqn_cntol.
maxiter (int) – Maximum iterations
alpha (ndarray) – Class-specific linearization parameters (R,)
npasses (int) – Number of Delta refresh passes (3 is the Chandy-Neuse fixed rule; pfqn_scat passes 1)

Returns:

Tuple of (Q, U, W, T, C, X, iterations), where T is the throughput PER REFERENCE VISIT (X broadcast over the stations the class visits), not the per-station throughput: visits are folded into L here and cannot be recovered from it. MATLAB returns no T at all and has the caller build T = V .* X; the slot exists because three Python callers unpack seven values. See the comment at the return statement.

Return type:

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

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

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

Parameters:

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

Return type:

PfqnMomlin

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

Bases: object

Result container for pfqn_momlin().

Variables:

Q (np.ndarray (M x R)) – Mean queue length.
X (np.ndarray (1 x R)) – Throughput per class.
R (U,) – Utilization and residence time.
QVar (np.ndarray (M x R)) – Queue-length variance.
QCov (np.ndarray (M x R x M x R)) – Queue-length covariance, QCov[i, r, j, s] = Cov[n_ir, n_js].
dQ (np.ndarray (M x R x M x R)) – Demand derivatives, dQ[i, r, j, s] = dQ_ir/dD_js.

class SchedStrategy(*values)[source]

Bases: Enum

Scheduling strategies for queueing stations.

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

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

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

Parameters:

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

Return type:

PfqnSens

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

Bases: object

Result container for pfqn_sens().

Variables:

R (X, Q, U,) – Base MVA measures. X is (1 x R) system throughput per class, Q (M x R) mean queue length, U (M x R) utilization, R (M x R) residence time per visit.
params (list of dict) – One entry per differentiation parameter p with keys type (‘L’ or ‘Z’), station (i, or -1 for Z) and jobclass (r).
dX (np.ndarray (R x P))
dR (dQ, dU,) – Derivative of each base measure w.r.t. parameter p.
QCov (np.ndarray (M x R x M x R)) – Queue-length covariance, QCov[i, r, j, s] = Cov[n_ir, n_js] = D_js dQ_ir/dD_js. The same-station blocks (i == j) come from pfqn_sens_mva(); the cross-station ones are read off the Jacobian.
QVar (np.ndarray (M x R)) – Queue-length variance, QVar[i, r] = QCov[i, r, i, r].
QTotVar (np.ndarray (M,)) – Variance of the total queue length per station, QTotVar[i] = Var[sum_r n_ir].
QCovAsym (float) – Roundoff-level residual of the moment recursion, see pfqn_sens_mva().

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

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

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

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

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

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

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

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

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

Parameters:

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

Return type:

PfqnSensMva

Notes

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

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

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

Bases: object

Result container for pfqn_sens_mva().

Variables:

R (X, Q, U,) – Base MVA measures, identical entry by entry to pfqn_mva(L, N, Z, mi). X is (1 x R) system throughput per class, Q (M x R) mean queue length, U (M x R) utilization, R (M x R) residence time per visit.
QCov (np.ndarray (M x R x R)) – QCov[i, r, s] = Cov[n(i,r), n(i,s)], the queue-length covariance of classes r and s at station i. Symmetric in (r, s).
QVar (np.ndarray (M x R)) – QVar[i, r] = QCov[i, r, r] = Var[n(i,r)].
QTotVar (np.ndarray (M,)) – QTotVar[i] = Var[sum_r n(i,r)], i.e. sum_{r,s} QCov[i, r, s]. This is Theorem 3 of de Souza e Silva and Muntz, obtained here without a capacity derivative.
QCovAsym (float) – max |W(r,s) - W(s,r)| over the covariance entries before symmetrization. The two triangles come from differentiating two different classes’ MVA equations, so this is an independent residual of the recursion and should sit at roundoff; a large value signals a bug.

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

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

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

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

Parameters:

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

Returns:

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

Return type:

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

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

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

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

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

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

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

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

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

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

Parameters:

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

Return type:

PfqnSensMvaldmx

Notes

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

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

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

Bases: object

Result container for pfqn_sens_mvaldmx().

Variables:

R (X, Q, U,) – Base measures, identical to pfqn_mvaldmx(lambda, D, N, Z, mu, S). X is (R,) throughput per class, Q (M x R) mean queue length, U (M x R) utilization, R (M x R) residence time.
QCov (np.ndarray (M x R x R)) – QCov[i, r, s] = Cov[n(i,r), n(i,s)], the same-station block.
QCovFull (np.ndarray (M x R x M x R)) – QCovFull[i, r, j, s] = Cov[n(i,r), n(j,s)].
QVar (np.ndarray (M x R)) – QVar[i, r] = Var[n(i,r)].
QTotVar (np.ndarray (M,)) – QTotVar[i] = Var[sum_r n(i,r)].
QCovAsym (float) – max |QCovFull[i, r, j, s] - QCovFull[j, s, i, r]| before symmetrization. The two entries are produced by differentiating two different classes’ equations, so this residual is an independent check of the recursion and should sit at roundoff.

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

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

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

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

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

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

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

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

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

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

Parameters:

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

Return type:

PfqnSensMom

Notes

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

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

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

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

Bases: object

Result container for pfqn_sens_mom().

Variables:

R (X, Q, U,) – Base MVA measures, identical entry by entry to pfqn_mva(L, N, Z, mi). X is (1 x R), Q, U and R are (M x R).
m (np.ndarray (M, G)) – m[i, g] = E[Q_(i,g)], the mean queue length of group g at station i. COLLAPSED to (M,) in the default single-group case, where it is the per-station total.
dm (np.ndarray (M, G, M, G)) – dm[i, g, j, g2], the scaled first derivative of m_(i,g) with respect to the parameter of (j, g2). Collapsed to (M, M) when G == 1.
d2m (np.ndarray (M, G)) – d2m[i, g], the scaled pure second derivative with respect to the parameter of (i, g). Only that entry is needed by (3.2); the mixed second derivatives are not required for moments of a single Q_(i,g) and would cost an extra factor M*G to carry. Collapsed to (M,) when G == 1.
Var (np.ndarray (M, G)) – Var[Q_(i,g)]. Collapsed to (M,) when G == 1.
Cov (np.ndarray (M, G, M, G)) – Cov[i, g, j, g2] = Cov[Q_(i,g), Q_(j,g2)]. Collapsed to (M, M) when G == 1, where the group index carries no information.
M2 (np.ndarray (M, G)) – E[Q_(i,g)^2]. Collapsed to (M,) when G == 1.
M3 (np.ndarray (M, G)) – E[Q_(i,g)^3]. Collapsed to (M,) when G == 1.
Skew (np.ndarray (M, G)) – Skewness of Q_(i,g), i.e. the third central moment divided by Var^(3/2). NaN where the variance is zero (a deterministic queue length). Collapsed to (M,) when G == 1.
CovAsym (float) – max |x_j dm_i/dx_j - x_i dm_j/dx_i| before symmetrization of Cov. The two are distinct expressions that must agree, so this is a live residual of the recursion; expect roundoff.

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

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

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

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

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

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

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

Parameters:

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

Return type:

PfqnSensRespt

Notes

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

pfqn_busyp(alpha, mu, P, N, subnet, n, gamma=None, tol=1e-12)[source]

Mean busy period of order n for the subnetwork subnet.

The busy period of order n is the time from the instant a job entering the subnetwork finds n-1 jobs in it up to the next instant when fewer than n jobs remain in it.

Parameters:

alpha (array (J,)) – Relative arrival rates, the solution of x*P = x for a closed network and of x = gamma + x*P for an open one.
mu (array (J, K) or callable) – Load-dependent service rates, mu[j, k-1] with k jobs at node j, or a callable mu(j, kvec) when the rates do not saturate.
P (array (J, J)) – Routing matrix.
N (int or float) – Population; numpy.inf for an open network.
subnet (sequence of int) – Zero-based indexes of the nodes forming the subnetwork.
n (int or sequence of int) – Busy period order(s), 1 <= n <= N.
gamma (array (J,), optional) – External arrival rates; required for an open network.
tol (float) – Relative tolerance of the open-network tail truncation.

Returns:

b (np.ndarray) – Mean busy period duration(s), same shape as n.
lG (np.ndarray) – Log normalizing constants of the subnetwork.
lH (np.ndarray) – Log normalizing constants of the complement, empty for an open network.

pfqn_busyp_multiclass(alpha, mu, P, N, subnet, n, gamma=None, phi=None, tol=1e-12, jobclass=-1)[source]

Mean busy period of order n for the subnetwork, multichain.

Parameters:

alpha (array (J, R)) – Relative arrival rates, one column per chain.
mu (array (J, R)) – Service rates, the chain-r rate at node j.
P (array (J, J) or sequence of R such arrays) – Routing, shared by every chain or one matrix per chain.
N (array (R,)) – Population per chain; numpy.inf entries for an open chain.
subnet (sequence of int) – Zero-based node indexes forming the subnetwork.
n (int or sequence of int) – Busy period order(s), counting the jobs of every chain.
gamma (array (J, R), optional) – External arrival rates; required for an open network.
phi (array (J, K), optional) – Dimensionless load-dependent scaling; None means a single server.
tol (float) – Relative tolerance of the open-network tail truncation.
jobclass (int) – Zero-based chain whose own jobs are counted, or -1 to count every chain. A per-class order is bounded by that chain’s population, not by the total.

Returns:

Mean duration(s), and the log normalizing constants of the subnetwork and of its complement over the lattice (lH empty when open).

Return type:

b, lG, lH

pfqn_busyp_clw(alpha, mu, P, N, subnet, n, gamma=None, isdelay=None, method='clw')[source]

Mean busy period of order n for the subnetwork, via NC point evaluations.

Parameters:

alpha (array (J, R)) – Relative arrival rates, one column per chain.
mu (array (J, R)) – Service rates, the chain-r rate at node j.
P (array (J, J) or sequence of R such arrays) – Routing, shared by every chain or one matrix per chain.
N (array (R,)) – Population per chain; numpy.inf entries for an open chain.
subnet (sequence of int) – Zero-based node indexes forming the subnetwork.
n (int or sequence of int) – Busy period order(s), counting the jobs of every chain.
gamma (array (J, R), optional) – External arrival rates; required for an open network.
isdelay (array (J,) of bool, optional) – Infinite-server nodes; the rest are single servers.
method (str) – Method name of the normalizing-constant method used for the point evaluations.

Returns:

b – Mean busy period duration(s).

Return type:

float or np.ndarray

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

Bases: object

Result container for pfqn_sens_respt().

Variables:

U (X, Q,) – Base measures at population N. X is (1 x R), Q and U are (M x R).
W (np.ndarray (M, R)) – W[i, l] = E[W_(i,l)], the mean sojourn time per visit of a class-l job at station i. Zero where class l does not visit i.
WM (np.ndarray (M, R, tmax)) – WM[i, l, t-1] = E[W_(i,l)^t].
WVar (np.ndarray (M, R)) – Var[W_(i,l)] = E[W^2] - E[W]^2. Requires tmax >= 2.
WSkew (np.ndarray (M, R)) – Skewness of W_(i,l). Requires tmax >= 3; NaN if the variance is zero.
m (np.ndarray (M,)) – E[Q_i] at population N, the total queue length.
Var (np.ndarray (M,)) – Var[Q_i] at population N.
p (np.ndarray (M, max(b))) – p[i, j] = P[Q_i = j] at population N, for j = 0..b_i-1. These are the only marginal probabilities the b-server recursion needs, so the matrix is RAGGED: row i is meaningful only up to column b_i and is zero-padded out to max(b). A padded entry is not P[Q_i = j]; it is simply not computed. Read row i as p[i, :b[i]].
Wresid (np.ndarray (M, R)) – The residence time w_i(l) of the MVA recursion. The identity W[i, l] = Wresid[i, l] / V[i, l] is an independent check of the t = 1 case of (4.5).

pfqn_respt_ps_moments(S, N, Z, method='auto')[source]

Sojourn-time moments at the PS station of a closed terminal-driven system.

The system is a bank of terminals in series with a single processor-sharing CPU, with class-dependent exponential think times (mean Z[r]) and class-dependent exponential service times (mean S[r]), and N[r] jobs of class r cycling between the two.

Two routes to the moments are implemented, both from Mitra and Morrison (1983):

'exact': solves the linear system c'[A - q_J I] = -pi'B of Proposition 3 on the state space {n : 0 <= n <= K}, K being the population vector with the tagged class decremented by one. The moments are then E[W_J] = sum_n c(n) and (q_J/2) E[W_J^2] = sum_n (n'1+1) c(n). Exact to solver precision, at the cost of a linear solve of dimension prod_r (K[r]+1).
'asymptotic': evaluates the two leading terms of the asymptotic expansion in inverse powers of the large parameter Nexp = max_r Z[r]/S[r], E[W_J^2] ~ c0 + c1/Nexp, of Proposition 6. The cost is a linear system of dimension R, the number of classes, and is therefore independent of the populations. Note that the expansion parameter is the think-to-service ratio and NOT the population, so a model with short think times is expanded in a small parameter no matter how many jobs it holds.
'auto': (default) takes the exact route when the state space has at most AUTO_MAX states and the asymptotic route otherwise.

The asymptotic route requires the normal-usage condition alpha > 0. Where it fails and the exact route is not affordable, the entry of W and W2 is NaN and the result records ‘unavailable’; asking for ‘asymptotic’ explicitly in that regime raises rather than returning a blank.

Parameters:

S (array_like (R,)) – Per-class mean service times at the PS station, positive.
N (array_like (R,)) – Per-class populations, non-negative integers.
Z (array_like (R,)) – Per-class mean think times, positive where N > 0.
method (str) – ‘auto’ (default), ‘exact’ or ‘asymptotic’.

Returns:

W (np.ndarray (R,)) – Per-class mean sojourn times at the PS station.
W2 (np.ndarray (R,)) – Per-class second moments of the sojourn time.
out (PfqnResptPsMoments) – Route taken and expansion diagnostics.
A class with ``N[r] = 0` has no sojourn time` and its entries are NaN.

See also

qsys_mm1_ps: the open counterpart, exact in closed form.

class PfqnResptPsMoments(R)[source]

Bases: object

Result container for pfqn_respt_ps_moments().

Variables:

method (list of str) – Per-class route taken: ‘exact’, ‘asymptotic’, ‘unavailable’ or ‘none’.
c1 (c0,) – Coefficients of the expansion E[W^2] ~ c0 + c1/expansionParam, NaN on the exact route.
alpha (np.ndarray (R,)) – Per-class unutilized fraction 1 - sum_r lambda_r/q_r of the CPU in the corresponding open system.
nstates (np.ndarray (R,)) – Size of the exact state space that the tagged class would need.
expansionParam (float) – The large parameter Nexp = max_r Z(r)/S(r).

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

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

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

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

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

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

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

Parameters:

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

Return type:

PfqnSensLinearizer

Notes

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

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

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

Bases: object

Result container for pfqn_sens_linearizer().

Variables:

W (X, Q, U,) – Approximate base measures. X is (1 x R); Q, U and W are (M x R).
m (np.ndarray (M,)) – Approximate E[Q_i], the total queue length at station i.
dm (np.ndarray (M, M)) – dm[i, h] = x_h dm_i/dx_h, the scaled first derivative.
d2m (np.ndarray (M,)) – d2m[i] = x_i^2 d^2m_i/dx_i^2.
Skew (Var, Cov, M2, M3,) – The moments of (3.2), formed exactly as in pfqn_sens_mom() but from the approximate derivatives. Var, M2, M3 and Skew are (M,); Cov is (M x M).
CovAsym (float) – Raw asymmetry of Cov before symmetrization. Unlike the exact routines, this is NOT expected to sit at roundoff: the Linearizer fixed point does not enforce the symmetry that the product form guarantees, so this is a useful measure of the approximation error.
iter (int) – Total CORE iterations performed.

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

Returns:

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

XN: Class throughputs (1 x R).
QN: Mean queue lengths (M x R).
UN: Utilizations (M,), per STATION, i.e. 1-P_j(0). Matches MATLAB
    pfqn_mvald, which likewise reports utilization per station and
    not per class.

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

Return type:

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

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. Standard arrival theorem throughout; for the interlocked-flow correction of Franks (1999), Ch. 4, Eq. (4.7) call pfqn_mvams_ilock instead.

Parameters:

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

Returns:

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

XN: Class throughputs (1 x R).
QN: Mean queue lengths (M x R).
UN: Utilizations. Per STATION-CLASS (M x R) on every branch except
    the closed multiserver one, which delegates to pfqn_mvald and so
    reports per STATION (M,). This shape inconsistency is inherited
    from MATLAB pfqn_mvams, which behaves identically; the only
    caller (api/solvers/mva/handler.py) discards UN and recomputes
    utilization analytically, as solver_mva.m does.

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

lG: Log normalizing constant.

Return type:

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

References

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

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

MVA entry point for models carrying the interlocked-flow correction.

The interlock of Franks (1999), Ch. 4, Eq. (4.7) is defined only for closed single-server models, so that is the one shape accepted here; anything else is refused rather than served without the correction. Models with no interlock go to pfqn_mvams.

Parameters:

lambda_arr (ndarray) – Arrival rate vector (R,). Must be all zero (closed model).
L (ndarray) – Service demand matrix (M x R).
N (ndarray) – Population vector (R,). Must be finite (closed model).
Z (ndarray) – Think time vector (R,).
mi (ndarray | None) – Queue replication factors (M,) (default: ones).
S (ndarray | None) – Number of servers per station (M,) (default: ones). Must be all one.
IL (ndarray | None) – Interlock matrix (R x R), see pfqn_mva_ilock. Required.

Returns:

Tuple of (XN, QN, UN, CN, lG); lG is always NaN.

Return type:

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

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

Returns:

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

Return type:

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

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

Lower ABA (asymptotic bound analysis) bound on throughput.

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

Parameters:

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

Returns:

Lower bound on throughput.

Return type:

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

Upper asymptotic bound on throughput (Zahorjan-Balanced).

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

Parameters:

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

Returns:

Upper bound on throughput.

Return type:

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

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

Parameters:

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

Returns:

Lower bound on mean queue length at station i.

Return type:

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

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

Parameters:

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

Returns:

Upper bound on mean queue length at station i.

Return type:

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

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

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

Parameters:

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

Returns:

Lower bound on throughput.

Return type:

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

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

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

Parameters:

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

Returns:

Upper bound on throughput.

Return type:

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

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

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

Parameters:

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

Returns:

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

Return type:

Xlo

pfqn_harel_bounds(rho, N, Z=0.0, maxUB=0)[source]

Harel-Namn-Sturm throughput bounds of a single-class closed network.

These are the SHARP bounds of Harel, Namn and Sturm, “Simple bounds for closed queueing networks” (Queueing Systems 31, 1999), distinct from the ‘sb’ family in the BA solver: ‘sb’ uses only the first three power sums in closed form, whereas this family evaluates the normalizing constant exactly at small populations and extrapolates from it. Both cite the same paper; they are different results in it and neither subsumes the other.

With the power sums A_i = sum_j rho_j^i,

G(n) = h_n(rho), the complete homogeneous symmetric polynomial, TH(n) = G(n-1)/G(n), the exact throughput at population n, LB = N / (A_1 + (N-1) (A_N/A_1)^{1/(N-1)}), UB(n) = N / (A_1 + ((N-1)/(n-1)) (n/TH(n) - A_1)), 2 <= n <= N.

G(n) IS the normalizing constant of the closed load-independent network at population n, so it must equal pfqn_ca on the same demands and TH(n) must equal the exact pfqn_mva throughput at population n. G is evaluated by the Newton-Girard recurrence n G(n) = sum_{i=1..n} A_i G(n-i); the n <= 7 ceiling on the extrapolation point is kept from the reference.

Parameters:

rho (ndarray) – (k,) relative utilizations, all strictly positive.
N (int) – population, at least 1.
Z (float) – think time; must be zero.
maxUB (int) – largest extrapolation point; defaults to min(N, 7) when <= 0.

Returns:

Tuple (LB, UB, TH) with UB[n-1] the upper bound extrapolated from population n (UB[0] unset) and TH[n-1] the exact throughput at population n, n = 1..maxUB.

Return type:

Tuple[float, ndarray, ndarray]

pfqn_harel_lb(rho, N, Z=0.0)[source]

Harel-Namn-Sturm throughput lower bound of a single-class closed network.

LB = N / (A_1 + (N-1) (A_N/A_1)^{1/(N-1)}) with A_i = sum_j rho_j^i, from Harel, Namn and Sturm, “Simple bounds for closed queueing networks” (Queueing Systems 31, 1999). A nonzero think time is refused.

Parameters:

rho (ndarray) – (k,) relative utilizations, all strictly positive.
N (int) – population, at least 1.
Z (float) – think time; must be zero.

Returns:

The throughput lower bound at population N.

Return type:

pfqn_harel_ub(rho, N, n, Z=0.0)[source]

Harel-Namn-Sturm throughput upper bound of a single-class closed network.

Extrapolated from the EXACT throughput TH(n) = G(n-1)/G(n) at the small population n, UB(n) = N / (A_1 + ((N-1)/(n-1)) (n/TH(n) - A_1)). G is evaluated by the Newton-Girard recurrence; the n <= 7 ceiling is kept from the reference implementation. A nonzero think time is refused.

Parameters:

rho (ndarray) – (k,) relative utilizations, all strictly positive.
N (int) – population, at least 1.
n (int) – extrapolation point, 2 <= n <= min(N, 7).
Z (float) – think time; must be zero.

Returns:

The throughput upper bound at population N.

Return type:

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

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

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

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

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

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

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

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

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

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

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

pfqn_ssd(L, N, Z=0.0, nservers=None)[source]

Server-Station Disaggregation bounds (Suri-Dallery 1986, Thm 5), single-class multiserver. Returns (Xlo, Xhi).

With Z>0 the queueing terms carry the terminal-workload correction of Lazowska et al. 1984, Table 5.2; adding Z without it is not a bound.

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

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

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

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

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

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

pfqn_scb(L, N)[source]

Single-class bounds of multiclass networks (Dowdy et al. 1992, JACM 39(1)).

Returns (Xlo, Xhi, Ulo, Uhi), a bracket on the total throughput and on the per-device utilizations of the UNKNOWN multiclass system whose single-class counterpart has demand vector L at population N.

SEMANTICS DIFFER FROM EVERY OTHER pfqn_* BOUND. aba/bjb/gb/… bracket the exact solution OF THE GIVEN MODEL; this brackets the multiclass system that the given single-class model aggregates. The lower side is therefore the EXACT single-class solution, not an approximation of it.

Theorem 2 / Corollary 2: aggregating an R-class model into its single-class counterpart can only understate performance, U_k,1 <= U_k,R and X_1 <= X_R, and Corollary 1 makes the utilization ratio uniform, U_k,R/U_k,1 = X_R/X_1 for every k. Theorem 3 (their Expression 3) caps the relative throughput error at (m-1)/(N+m-1), m = min(N,K), independently of the demands. The single-server capacity U_k,R <= 1 caps the same ratio at 1/(X_1*max(L)), tight on the paper’s own worst case, so both are applied. L holds queueing stations only: Theorem 3 rests on the delay-free balanced-network throughput, so a delay station is not admitted.

pfqn_scbgap(N, K, r=None, undominated=False)[source]

Maximum relative throughput error of merging r of N classes (Dowdy 1992).

Demand-free bound on the relative throughput error incurred when r of the N single-customer classes of a closed product-form network are merged into one class. With r = N (the default) this is the full single-class aggregation error of their Theorem 3, at most 50%; with r < N it is the partial-aggregation error of their Theorem 4. The bound never reads the demands, so it can be attached as a certified error bar to any result computed on merged chains.

General case, dominating classes allowed (Expression 4, and with r = N Expression 3): e = (min(r,K)-1)/(r+min(r,K)-1). Undominated case, every customer placing the same total demand (Theorem 5 and its comment (3), which lifts the N = R restriction): e = r(r-1)/(min(N,K)(2r-1)), valid for r <= K only, smaller than the general case by the factor r/min(N,K) and equal to it at r = K. THE DOMAIN IS NOT COSMETIC: Theorem 5 gives each of its R classes a dedicated device, so r never exceeds K there, and comment (3) states the generalization for r < K. Evaluated at r > K the expression climbs past the general bound and past the 50% cap of Theorem 3, i.e. it stops being a bound, so r > K is refused rather than returned.

pfqn_usumbound(R, K, N)[source]

Upper bound on sum_k U_k in a closed R-class network (Dowdy 1992, Thm 6).

sum_k U_k,R <= (H-1) + (K-H+1)(N-H+1)/(K+N-2H+1), H = min(R,K). Demand-free and nondecreasing in R, which is what makes it invertible into a lower bound on the number of necessary classes; see pfqn_minclasses. At R >= min(N,K) it reaches min(N,K), the trivial one-busy-server-per-device cap. The paper’s worked case is K = 2, N = 3, R = 1, giving 2N/(N+1) = 1.5.

pfqn_minclasses(Usum, K, N)[source]

Lower bound on the class count from a measured utilization sum (Dowdy 1992).

Smallest number of customer classes R consistent with an observed sum of device utilizations, obtained by inverting the demand-free Expression (6) bound of pfqn_usumbound, which is nondecreasing in R. Only measured quantities are needed – the utilizations, the device count and the population – so the answer is available BEFORE any class-specific demand has been characterized. An upper bound on R is meaningless (extra classes can always be introduced by splitting) and none is returned.

Returns NaN when Usum exceeds min(N,K) and so is unattainable by ANY class structure, which signals a measurement error rather than a workload needing more classes. The paper’s example: K = 2, N = 3, Usum = 1.6 -> 2, since a single class admits at most 2N/(N+1) = 1.5.

pfqn_explicit(L, N, tol=None, method='auto', maxloss=np.inf)[source]

Explicit closed-form normalizing constant of a multiclass closed network.

Evaluates the two explicit expressions of Casale, “Accelerating Performance Inference over Closed Systems by Asymptotic Methods”, ACM SIGMETRICS 2017, Eqs. (15) and (16). Both instantiate the divided-difference form of Corollary 3.2,

G(N) = sum_{0<=t<=N} (-1)^(|N|-|t|)/(N_1!…N_R!) prod_r C(N_r,t_r) g_t(|N|)

by substituting a closed form for the single-class constant g_t(|N|) at the induced demands theta_k(t) = sum_r t_r L(k,r). Eq. (15) is Gordon’s partial fraction and needs the induced demands PAIRWISE DISTINCT; Eq. (16) is the general partial-fraction expansion over the distinct values and their multiplicities, and reduces to Eq. (15) when every multiplicity is one. The choice is automatic: Eq. (16) is used as soon as two induced demands are closer than tol relative to the largest one at that t.

SINGLE CLASS. At R=1 the multiclass constant IS the single-class constant at demands L, so the outer sum is skipped: g_t(N) = t^N g_1(N) and sum_t (-1)^(N-t) t^N/(t!(N-t)!) = S(N,N) = 1. Running the difference anyway would add N alternating terms, and their cancellation, to a closed form that carries none of them. What is left is O(K^2) work at any population.

Only single-server load-independent queues are admissible: infinite servers need the integral form of Corollary 3.4 and load-dependent rates need the load-dependent generalization of the outer sum.

NUMERICS. Both expressions alternate in sign with terms far larger than the result, so they are evaluated as signed log-sum-exps: this removes the floating-point RANGE problem but not the cancellation, which is what makes multiprecision arithmetic necessary on all but small models.

Parameters:

L – Service demand matrix (KxR) of single-server load-independent queues.
N – Population vector (1xR).
tol (float | None) – Relative tolerance declaring two induced demands redundant (default: machine epsilon).
method (str) – ‘auto’ (default), ‘distinct’ to force Eq. (15), ‘repeated’ to force Eq. (16).
maxloss (float) – Cancellation budget in decimal digits. Finite values turn the warnings into a silent REFUSAL (lG=nan) once the budget is exceeded, for callers that hold a fallback; default inf keeps the warnings.

Returns:

the logarithm of the normalizing constant, the constant, the expression actually used (‘distinct’ or ‘repeated’), and the decimal digits lost to cancellation.

Return type:

(lG, G, method, lossDigits)

pfqn_explicit_ld(L, N, mu=None, tol=None, method='auto', maxloss=np.inf)[source]

Explicit closed-form normalizing constant of a multiclass LLD network.

Load-dependent counterpart of pfqn_explicit. It evaluates the same divided-difference form of Casale (SIGMETRICS 2017), Corollary 3.2,

G(N) = sum_{0<=t<=N} (-1)^(|N|-|t|)/(N_1!…N_R!) prod_r C(N_r,t_r) h_t(|N|)

but substitutes for the single-class constant h_t(|N|) the LIMITED LOAD-DEPENDENT closed form of Casale, Harrison and Ong (Perform. Eval. 2021), Theorem 1, Eq. (8),

h_theta(N) = sum_{0<=v<s} g_sigma(N-|v|) prod_k phi_k(v_k) phi_k(v_k) = theta_k^v_k / prod_{t=1..v_k} alpha_k(t) * (1 - alpha_k(v_k)/alpha_k(s_k))

at the induced demands theta_k(t) = sum_r t_r L(k,r). Here alpha_k(.) = mu(k,.) is the load-dependent scaling of station k, s_k the population past which it stays constant, sigma_k = theta_k/alpha_k(s_k) the SCALED demands, and g_sigma the FIXED-RATE single-class constant at those scaled demands, which is exactly what pfqn_explicit evaluates in closed form (Eqs. 15 and 16). The result is therefore explicit throughout, with no recursion over population.

Two conventions of Theorem 1 are not those of the equilibrium distribution and are easy to get wrong. alpha_k(0) is taken as ZERO inside the bracket of phi_k, so that phi_k(0) = 1, even though the state probabilities use alpha_k(0) = 1; and g_sigma(n) = 0 for n < 0, which caps the outer sum at |v| <= |N|. With alpha_k(n) = min(n,s_k) the expression collapses to Gordon’s multi-server formula, Oper. Res. 38(5), 1990, Eq. (29), but unlike that one it needs neither a multi-server shape nor distinct scaled demands.

LIMITED LOAD DEPENDENCE. Theorem 1 holds for any s_k with alpha_k(n) = alpha_k(s_k) for all n >= s_k, and a LARGER s_k is always admissible, so s_k is detected here as the smallest index whose value the tail of mu(k,:) repeats to within tol. A station whose rates never settle (an infinite server, mu(k,n) = n) gets s_k = |N|, which is still exact: populations above |N| do not occur, so redefining alpha_k there changes nothing. It is merely expensive, since the inner sum costs prod_k s_k terms, capped by |v| <= |N|. Think time is not admissible: a delay would have to enter g_sigma, whose closed form covers queues only.

NUMERICS. Both sums alternate in sign with terms far larger than the result, so they are evaluated as signed log-sum-exps. phi_k is sign-definite when alpha_k increases, as a multi-server station does, and changes sign where alpha_k decreases, so a decreasing rate function costs digits in the inner sum too.

SINGLE CLASS. At R=1 the divided difference is the identity, since h_theta(N) is homogeneous of degree N in theta exactly as in the fixed-rate case, so the outer sum is skipped and Theorem 1 is evaluated once at theta = L.

Parameters:

L – Service demand matrix (MxR).
N – Population vector (1xR).
mu – Load-dependent rate matrix (Mx sum(N)), alpha_i(j) = mu[i,j-1]; default all ones.
tol (float | None) – Relative tolerance declaring two scaled demands redundant, and the rate tail constant (default: machine epsilon).
method (str) – ‘auto’ (default), ‘distinct’ to force Eq. (15), ‘repeated’ to force Eq. (16).
maxloss (float) – Cancellation budget in decimal digits. Finite values turn the warnings into a silent REFUSAL (lG=nan) once the budget is exceeded, for callers that hold a fallback; default inf keeps the warnings.

Returns:

the logarithm of the normalizing constant, the constant, the expression used for g_sigma (‘distinct’ or ‘repeated’), and the decimal digits lost to cancellation.

Return type:

(lG, G, method, lossDigits)

pfqn_rgf(L, N, Z=0.0)[source]

Exact normalizing constant of a single-class closed product-form network by convolving per-node generating-function sequences.

A GROUP of m stations sharing the same demand p collapses into the single negative-binomial sequence r(k) = C(k+m-1,k) p^k, so the whole group costs one sequence rather than m convolution passes; the delay contributes the Poisson sequence Z^k/k!. Cost O(G N^2) against Buzen’s O(M N), so RGF is the cheaper route on heavily replicated models with moderate populations (G N < M). The recursion runs entirely in the log domain.

Parameters:

L – Service demand vector (M,) of the queueing stations.
N – Population (nonnegative integer scalar).
Z (float) – Think time (scalar, default 0).

Returns:

Tuple (G, lG, lg) with lg the vector of log g(0), …, log g(N).

Return type:

Tuple[float, float, ndarray]

pfqn_rgfmc(L, N, Z=None, tol=1e-12, maxterms=1000000, maxcancel=15.0)[source]

Exact multiclass normalizing constant by recursion on generating functions.

Eliminates one class at a time by residues (Harrison-Coury 2002 Thm 1, algorithmised as Harrison-Lee 2004 eqs. 4-5) until a single class is left, then finishes with the multiplicity-aware convolution of Coury-Harrison 1997 Property 1, memoised by load vector as Harrison-Lee sec. 3.4.

THINK TIMES ARE NOT IN EITHER PAPER. Their generating function is the rational prod_i (1 - rho_i z)^-m_i, and an infinite server multiplies it by the entire exp(sum_r Z_r z_r), which destroys the residues-sum-to-zero identity (Bertozzi-McKenna 1993 fact IV, p. 246) the recursion rests on. The delay is therefore carried by their truncation, eqs. (3.19)-(3.21): only the first k_r+1 Taylor coefficients of exp(Z_r z_r) can reach [z_r^k_r], so replacing it by that polynomial is EXACT, not an approximation, and restores a rational integrand. The cost is that the eliminated class’s population re-enters the term count, which is precisely the population-insensitivity Harrison-Lee sec. 4 advertises; the class kept for the base case pays nothing, so the base is chosen to be a populous one only when conditioning allows.

The recursion is exact in exact arithmetic but is an ALTERNATING sum over residues, so near-coincident loads over the eliminated class destroy significance. maxcancel bounds the nats of cancellation tolerated and the routine REFUSES beyond it rather than returning a confidently wrong lG.

Parameters:

L – Service demand matrix (M x R).
N – Population vector (R,), nonnegative integers.
Z – Think time vector (R,), default zeros.
tol (float) – relative tolerance for calling two affine forms proportional.
maxterms (int) – cap on residue terms carried between eliminations.
maxcancel (float) – nats of cancellation tolerated before refusing.

Returns:

Tuple (G, lG).

Return type:

pfqn_gerasimov(L, N, Z=None, tol=1e-12, maxterms=200000)[source]

Exact normalizing constant of a closed multiclass product-form network by ITERATED RESIDUES of its rational generating function, one class at a time.

Gerasimov (1995) evaluates

G(N_1,…,N_R) = (2 pi i)^-R int_G1 … int_GR
prod_s z_s^(N_s-1) prod_i (1 - sum_s x_is/z_s)^-1

by residues, and gives the resulting CLOSED FORM only for R = 1 (Thm 1-2) and R = 2 (Thm 3 for simple poles, Thm 4 for multiple ones), stating that “for three or more classes of customers, the normalizing constants can be found by numerical methods”. This routine implements the residue elimination itself, so the closed form is produced for ANY R; at R = 2 it reproduces Thm 3/4 term by term.

Written as a coefficient of the u_s = 1/z_s series,

G(N) = [prod_s u_s^(N_s)] exp(sum_s Z_s u_s) prod_i (1 - sum_s x_is u_s)^-1,

every factor is AFFINE in u, so singling out u_r gives f = A - B u_r with A affine in the surviving variables. Partial fractions in u_r map a sum of products of affine powers into another one with one variable fewer, and R-1 such steps leave a univariate coefficient extraction. At R = 2 the single step returns one term per station i, with outer factor x_i2^(N_2+M-1) / prod_{k!=i}(x_i2-x_k2), a pole of order N_2+1 at x_i1 and simple poles at the paper’s z_1ik = (x_k1 x_i2 - x_i1 x_k2)/(x_i2 - x_k2): exactly Thm 3, with the multiple poles of Thm 4 (his xi_i < M) handled by the same step. Tied x_i2, vanishing x_i2 and identical station rows, all outside the paper’s hypotheses, are ordinary cases here.

Cost. Let M be the number of stations and order the populations N_(1) <= … <= N_(R). The first elimination turns the single input term into M, and every later one multiplies the count by C(S+M-1,M-1) + M-1, where S is the total population already eliminated: a pole of order S+1 has to be differentiated against the M-1 remaining ones. The innermost extraction then convolves M series of length N_(1). Hence R = 1 costs O(M N), Buzen’s own cost; R = 2 costs O(M^2 N_(1)^2), INDEPENDENT OF N_(2); and R >= 3 costs the same times prod_{r=3}^{R} C(N_(r)+M-1, M-1). The R = 2 line is the reason to reach for this method: a population removed by residues enters only as a pole ORDER, i.e. through binomial coefficients, so it costs nothing at all. On a 4-station two-class model at N = [6, 20000] this returns lG in 0.4 ms where pfqn_ca needs 1.8 s, to the same 1.3e-16. For R >= 3 the term count is polynomial in the populations of degree (M-1)(R-2) and exponential in R, which is why the paper stops at two classes and why maxterms exists.

Conditioning. The sum is alternating, exactly as the paper writes it, and two decisions keep it usable: near-coincident poles are merged under a RELATIVE tolerance, so they are one multiple pole rather than two nearly cancelling simple ones, and the class left for the innermost extraction is the one with the SMALLEST population, because that population is the degree the final, sign-indefinite series is carried to. Measured on 372 random models against pfqn_ca: median 2.0e-16, p90 4.6e-15, p99 1.1e-11, worst 1.5e-09. On an ill-conditioned demand matrix pfqn_ca or pfqn_nc are still the safer routes to the same number.

Parameters:

L – service demand matrix (M x R), L[i,r] = demand of class r at station i.
N – population vector (R,), nonnegative integers.
Z – think time vector (R,), default zeros. A delay contributes the entire factor exp(sum_s Z_s u_s), handled exactly by convolving its Poisson coefficients into each elimination.
tol (float) – relative tolerance for declaring two affine forms proportional, hence one pole rather than two.
maxterms (int) – cap on the number of residue terms carried between eliminations. Exceeding it is an error, not a truncation: a truncated residue sum is not a bound or an approximation of G, it is a wrong number.

Returns:

(G, lG) the normalizing constant and its logarithm.

Return type:

pfqn_dnc(L, N)[source]

Normalizing constant and throughput at a REAL-VALUED population, by partial-fraction inversion of the network generating function.

With distinct loads x_1..x_G of multiplicities m_1..m_G the generating function prod_g (1-x_g u)^{-m_g} expands as:

G(n) = sum_g sum_{j=1..m_g} A_gj C(n+j-1,j-1) x_g^n,

every term of which is analytic in n, so evaluating at a real n interpolates the integral normalizing constants exactly and gives a smooth throughput curve X(N) = G(N-1)/G(N) through the integral points. For all-distinct loads A_g = prod_{l!=g} x_g/(x_g - x_l) is used directly; with repeated loads the coefficients are recovered from G(0..M-1).

Only the queueing part admits this continuation: the delay sequence Z^n/n! is entire and has no partial-fraction expansion, so a think time is not accepted here. Use pfqn_nintmva for nonintegral populations with a delay.

Parameters:

L – Service demand vector (M,) of the queueing stations.
N – Population (real nonnegative scalar; may be fractional).

Returns:

Tuple (X, G, lG).

Return type:

Tuple[float, float, float]

pfqn_nintmva(L, N, Z=0.0)[source]

Exact MVA recursion started from the FRACTIONAL base n0 = N - floor(N), giving mean performance measures at a real-valued population (“aMVA”).

The recursion is the standard Reiser-Lavenberg one, stepped in unit increments from n = n0 (where the arrival-theorem term is taken as 0, the network below the base being empty) up to n = N. At integer N the base is 0 and the recursion is bit-identical to exact MVA; at fractional N it interpolates smoothly through the integral points.

Unlike pfqn_dnc this accepts a think time. It is single-class: for fractional multiclass populations use pfqn_bs, which accepts them directly.

Parameters:

L – Service demand vector (M,) of the queueing stations.
N – Population (real nonnegative scalar; may be fractional).
Z (float) – Think time (scalar, default 0).

Returns:

Tuple (X, Q, U, R).

pfqn_mva_interval(L, N, Z=0.0)[source]

Exact hull of single-class MVA over an input box.

Single-class MVA is monotone in every input: the throughput decreases in each demand and in the think time and increases in the population, the per-station queue length and residence time increase in the own demand and in the population and decrease in the other demands and in the think time, and the totals increase in every demand and in the population and decrease in the think time (Luthi and Haring 1998, Theorems 2-5, Table 1). By their Theorem 1 the exact range of a function monotone in each argument is attained at the endpoints of the input box, so each bound below is one ordinary MVA call at the corner that the sign pattern selects. This is the algorithm of their Fig. 2 and it costs 2*(m+2) MVA calls, m being the number of thick demand intervals; evaluating the MVA recursion in interval arithmetic instead would be a valid but far wider enclosure, since every input recurs at each step (the dependency problem, 14x too wide on the paper’s own example).

The returned interval is the exact hull of MVA over the input box, not a bound on the true network: it holds conditionally on the demands lying in the box, and says nothing about the accuracy of MVA itself. It must therefore not be composed with the brackets of SolverBA, which bracket the exact solution of a model whose demands are known.

Delay stations are folded into Z, exactly as in pfqn_mva: a delay demand interval enters as a term of the think-time interval, and the hull of the sum is the sum of the hulls when the delays vary independently. Load-independent single-server queueing stations only, one class only; the monotonicity theorems cover no other case.

Parameters:

L – Service demand intervals (M, 2), column 0 lower, column 1 upper. An (M,) vector is read as a thin box.
N – Population interval [nlo, nup], or a scalar for a thin population.
Z – Think time interval [zlo, zup], or a scalar (default 0).

Returns:

PfqnMvaIntervalResult with fields X, Q, U, R, Rtot, Qtot.

class PfqnMvaIntervalResult(X, Q, U, R, Rtot, Qtot)[source]

Bases: NamedTuple

Interval-valued mean performance measures. Every field is [lower, upper].

Create new instance of PfqnMvaIntervalResult(X, Q, U, R, Rtot, Qtot)

X: ndarray: Alias for field number 0

Q: ndarray: Alias for field number 1

U: ndarray: Alias for field number 2

R: ndarray: Alias for field number 3

Rtot: ndarray: Alias for field number 4

Qtot: ndarray: Alias for field number 5

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

Approximate MVA whose arrival-instant queue lengths come from the THROUGHPUT ELASTICITIES rather than from a population-shift heuristic.

Let E_mkc = (D_mk/X_c) dX_c/dD_mk be the elasticity of the class-c throughput with respect to the class-k demand at station m. Tay shows that the elasticities satisfy the R linear equations

E_mkj sum_t B_tj Q_jt (1+Q_jt) =

-[(delta_jk + Q_jm) B_mk Q_km

sum_{c!=j} E_mkc sum_t B_tc Q_jt Q_ct]

with B_ir = 1/(1 + D_ir X_r/N_r), and that the arrival-instant queue length is then simply Q_km^(r) = Q_km + E_mkr, which closes the MVA recursion R_rm = D_rm (1 + sum_k Q_km^(r)). One R x R solve per (station, class) pair per iteration.

Delay stations enter through Z only. They are “AS” servers in the survey’s notation (d_t = 0), so they contribute Z_j X_j to the denominator of the elasticity equations but nothing to its numerator.

Parameters:

L – Service demand matrix (M x R).
N – Population vector (R,).
Z – Think time vector (R,). Default: zeros.
tol (float) – Convergence tolerance on the queue lengths.
maxiter (int) – Maximum number of iterations.
QN0 (ndarray | None) – Initial guess for the queue lengths (M x R).

Returns:

Tuple (XN, QN, UN, RN, it, QNarr). QNarr[m,k,r] is the class-k queue length at station m as seen by an arriving class-r job: the auxiliary quantity the method is tabulated on, NOT the queue length of the model re-solved at N - e_r (the same object only for an exact solution).

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

SCAT approximate MVA.

Parameters:

L (ndarray) – Demand matrix (M x R)
N (ndarray) – Population vector (R,)
Z (ndarray) – Think time vector (R,)
sched_type (List[str] | None) – Scheduling strategy per station; accepted for interface parity with pfqn_linearizer, but the residence-time recursion is discipline-independent
tol (float) – Convergence tolerance
maxiter (int) – Maximum inner iterations
QN0 (ndarray) – (M x R) warm start for the Bard-Schweitzer initialization

Returns:

(Q, U, W, T, C, X, iterations), where T is the throughput PER REFERENCE VISIT.

Return type:

Same tuple as pfqn_linearizer

pfqn_hst(L, N, Z=0.0, ist=None)[source]

Robustness certificate for a single-class closed product-form solution: how far the predicted throughput can move when the homogeneous-service-time (HST) assumption fails at one station.

HST states that the mean service time at station i does not depend on the queue length there. Suri perturbs it to S_i(n) = S_i (1 + a_n), one relative deviation per queue-length level, and shows (eq. 3.11) that to first order:

[(1/X0) dX0/da_n] = c_n = P(n_i >= n+1)/u_i - P(n_i >= n),

with u_i = L_i X0 and P(n_i >= n) = L_i^n G(N-n)/G(N). The naive certificate |dX0/X0| <= (sum_n |c_n|) d follows from |a_n| <= d alone, and by Lemma 3.1 that total equals Q_i(N) - Q_i(N-1).

That bound is loose because an operationally consistent perturbation must leave the observed mean service time unchanged, sum_n p_n a_n = 0. The constrained problem (P1):

max |sum_n c_n a_n|  s.t.  |a_n| <= d,  sum_n p_n a_n = 0

is a one-constraint linear program, solved here exactly: its optimum sets a_n = +/-d according to whether c_n/p_n exceeds a threshold, with at most one fractional coordinate.

Parameters:

L – Service demand vector (M,) of the queueing stations.
N – Population (nonnegative integer scalar).
Z (float) – Think time (scalar, default 0).
ist (int | None) – 0-based station index the perturbation applies to (default: the bottleneck, argmax L).

Returns:

Dict with keys station, X, U, Q, Pgeq, p, c, total, worst, astar.

Return type:

dict

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

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

Multiclass (L with >1 column) raises NotImplementedError.

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

Logistic Expansion (LE) asymptotic approximation for normalizing constant.

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

Parameters:

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

Returns:

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

Return type:

Tuple of (Gn, lGn)

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

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

The common corrected asymptotic expansion (LE-KT), computed on the cheaper side.

The corrected logistic expansion (pfqn_ble) and the corrected Knessl-Tier expansion (pfqn_bkt) are ONE estimator, evaluated in M-1 and in R dimensions. With a think time their stationary points are one point in dual coordinates, xi_r = N_r/(Z_r + v u’L_r) being the class throughputs of the LE fixed point and v u_k = 1/(1-U_k) the M/M/1 factor of the KT saddle, and Sylvester’s identity exchanges the R x R Hessian determinant for the M x M one, after which every 2 pi cancels; they agree to the accuracy of the two saddle-point solvers (~1e-7 nats, 1e-14 with polished saddles). Without a think time the LE branch integrates the radius exactly as Gamma(N+M) while KT Laplaces it, so the two differ by the constant (1-log(2 pi)/2) - r(N+M), r the Stirling remainder of a Gamma direction; the common estimator is defined as the KT value, and the LE side here carries M(1-log(2 pi)/2) - r(N+M) rather than pfqn_ble’s (M-1)(1-log(2 pi)/2). The route is pfqn_lekt_route. See _kb/03-api-layer.md.

Parameters:

L (ndarray) – Service demand matrix (MxR).
N (ndarray) – Population vector (1xR).
Z (ndarray) – Think time vector (1xR), optional.

Returns:

(Gn, lGn), the normalizing constant and its logarithm.

Return type:

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

The side pfqn_lekt computes on: ‘kt’ when R <= M or a class self-loops (one nonzero demand and no think time, which pfqn_kt extracts exactly), ‘le’ otherwise. The KT side is an R-dimensional convex solve and an R x R determinant, the LE side an M-dimensional fixed point and an (M-1) x (M-1) one.

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

Logistic expansion with the eps->0 bias correction (BLE).

Cas17 Theorem 4.1 holds for eps >= eps_N > 0, where the K(1+eps*N) self-looping populations make the integrand concentrate. Evaluated at eps->0, as pfqn_le does, the curvature at the saddle tends to 1 rather than growing with N, so Laplace’s method has no asymptotic regime there and carries an O(1) relative bias of e/sqrt(2*pi) PER LAPLACED DIRECTION. The count is the exponent on sqrt(2*pi) in the branch taken: M-1 with Z=0, where the radial integral is exact as Gamma(N+M), and M with Z>0, where the radius is Laplaced too. Measured over the 1562 models of the Cas17 dataset (Zenodo 546873, sec5.3.1, sigma=100) the Z>0 deficit is M to within 0.01 units. The published expansion is NOT in error; the correction is EMPIRICAL and is not part of Cas17. See _kb/03-api-layer.md.

Parameters:

L (ndarray) – Service demand matrix (MxR).
N (ndarray) – Population vector (1xR).
Z (ndarray) – Think time vector (1xR), optional.

Returns:

(Gn, lGn), the normalizing constant and its logarithm.

Return type:

pfqn_aghq(L, N, Z=None, q=3)[source]

Adaptive Gauss-Hermite quadrature of the McKenna-Mitra integral.

Rescaling the simplex integral by the LE mode and curvature, w = w* + A^-1/2 z, and applying the q-node probabilists’ Gauss-Hermite rule in each of the M-1 directions gives a convergent rule whose q=1 member is pfqn_le itself (single node at the mode, weight sqrt(2*pi)), to the tolerance of the shared fixed point. Cost is q^(M-1) evaluations, which confines the method to small M.

A tensor rule is not invariant to the choice of A^-1/2: any B with B B’ = inv(A) is admissible and they place the nodes differently. The principal-axis frame from the eigendecomposition is used, as in the reference results; where two curvatures are close to equal the frame is close to arbitrary and two valid rules can part company well above their own error, converging back together as q grows. Do not compare across codebases node by node.

With Z>0 the radius is integrated numerically and the rule is applied to the M-1 simplex directions, so every node costs one radial quadrature. q=1 there is LE with an exact radius, NOT pfqn_le’s own Z>0 branch.

Parameters:

L (ndarray) – Service demand matrix (M x R).
N (ndarray) – Population vector (R,).
Z (ndarray | None) – Think time vector (R,). Optional.
q (int) – Nodes per simplex direction (default 3).

Returns:

Tuple of (Gn, lGn).

Return type:

Reference:: J. McKenna, D. Mitra. “Integral representations and asymptotic expansions for closed Markovian queueing networks: normal usage.” BSTJ 61(5), 1982. 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 (ndarray) – Service demand matrix (M x R).
N (ndarray) – Population vector (R,).
Z (ndarray | None) – Think time vector (R,). Optional.
order (int | None) – Degree of cubature rule (default: ceil((sum(N)-1)/2)).
atol (float) – Absolute tolerance (default: 1e-8).

Returns:

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

Return type:

Tuple of (Gn, lGn)

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

pfqn_mci(D, N, Z=None, I=100000, variant='imci')[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 (ndarray) – Service demand matrix (M x R).
N (ndarray) – Population vector (R,).
Z (ndarray | None) – Think time vector (R,). Optional, defaults to zeros.
I (int) – Number of samples (default: 100000).
variant (str) – MCI variant - ‘mci’, ‘imci’ (improved), ‘amci’, ‘lhsmci’ or ‘rm’ (repairman). Default: ‘imci’. ‘amci’ and ‘lhsmci’ use the ‘imci’ tilt and differ only in how the uniforms are drawn. ‘amci’ draws ANTITHETIC pairs (u, 1-u). This does NOT reliably reduce variance here: the tilted integrand is not monotone in the exponential draws (the tilt term -(1-gamma)V decreases while the N log(VD+Z) term increases), so the pair correlation is not systematically negative; measured variance ratios against ‘imci’ range from 0.54 to 1.6 across models. It is kept because it is the Ross-Wang construction, not because it is the better default. ‘lhsmci’ stratifies each coordinate by Latin hypercube sampling, which IS reliably variance-reducing on the same models (ratios 0.0 to 0.48, exact quadrature in the limit of one station) at O(I log I) extra cost.

Returns:

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

Return type:

Tuple of (G, lG, lZ)

Reference:: Implementation based on MonteQueue methodology.

pfqn_grnmol(L, N)[source]

Normalizing constant using Grundmann-Moeller quadrature.

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

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

Parameters:

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

Returns:

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

Return type:

Tuple of (G, lG)

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

pfqn_le_fpi(L, N)[source]

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

Public wrapper for the internal _pfqn_le_fpi function.

Parameters:

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

Returns:

Mode location vector u (M,).

Return type:

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

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

Public wrapper for the internal _pfqn_le_fpiZ function.

Parameters:

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

Returns:

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

Return type:

Tuple of (u, v)

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

Compute Hessian matrix (no think time case).

Public wrapper for the internal _pfqn_le_hessian function.

Parameters:

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

Returns:

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

Return type:

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

Compute Hessian matrix (with think time case).

Public wrapper for the internal _pfqn_le_hessianZ function.

Parameters:

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

Returns:

Hessian matrix (M x M).

Return type:

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

Returns:

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

Return type:

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

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

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

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

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

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

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

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

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

Parameters:

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

Returns:

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

Return type:

Examples

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

See also

pfqn_is, pfqn_ncld, pfqn_nc

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

Normalizing constant and mean measures for mixed open/closed networks with limited load dependence.

Parameters:

lam (ndarray) – Arrival rate vector (R,) - 0 on closed classes
D (ndarray) – Service demand matrix (M x R)
N (ndarray) – Population vector (R,) - inf for open classes
Z (ndarray | None) – Think time vector (R,), optional
mu (ndarray | None) – Load-dependent rate matrix (M x >= sum(N_closed)), optional
S (ndarray | None) – Number of servers per station (M,), kept for signature parity
options (Dict[str, Any] | None) – Solver options forwarded to pfqn_ncld

Returns:

PfqnNcldmxResult with the closed-conditional constant (G, lG), the open-class normalizing prefactor lGopen, the method used for the closed-conditional solve, and the mean throughputs XN (1 x R) and queue lengths QN (M x R).

Return type:

PfqnNcldmxResult

class PfqnNcldmxResult(G, lG, lGopen, method='default', XN=None, QN=None)[source]

Bases: object

Result of a mixed limited load-dependent normalizing constant computation.

G: float

lG: float

lGopen: float

method: str = 'default'

XN: ndarray | None = None

QN: ndarray | None = None

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:

L (ndarray) – Service demands at all stations (M x R)
N (ndarray) – Number of jobs for each class (1 x R)
mu (ndarray) – Load-dependent scalings (M x Ntot)
options (Dict[str, Any] | None) – Solver options

Returns:

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

Return type:

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.

The recursion is

g(m,n,t) = g(m-1,n,1) + L_m/mu(m,t) * g(m,n-1,t+1)

whose right-hand side, at a fixed t, lives entirely at t+1. Sweeping t downward therefore advances every population at once through one shifted logaddexp, which is why the population axis is a vector here rather than the innermost loop of a scalar triple loop. Only the t=1 slice is carried from one station to the next. The work is still O(M*Ntot^2), but it is O(M*Ntot) array operations instead of O(M*Ntot^2) interpreted ones, and the values are bit-identical to the scalar order of accumulation.

Parameters:

L (ndarray) – Service demands at all stations (M x 1)
N (ndarray) – Number of jobs (scalar or 1x1 array)
mu (ndarray) – Load-dependent scaling factors (M x Ntot)
options (Dict[str, Any] | None) – Solver options (unused, for API compatibility)

Returns:

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

Raises:

RuntimeError – If multiclass model is detected

Return type:

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

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

Same recursion, same arithmetic and bit-identical results to pfqn_gldsingle, but with the rate-offset axis truncated at the limited load-dependence threshold instead of at the population. Unrolling

g(m,n,t) = g(m-1,n,1) + L_m/mu(m,t) * g(m,n-1,t+1)

gives

g(m,n,t) = sum_{j=0..n} prod_{i=0..j-1} L_m/alpha_m(t+i) * g(m-1,n-j,1)

so once t >= s_m, where s_m is the population past which alpha_m stays constant, every factor is alpha_m(s_m), the product collapses to (L_m/alpha_m(s_m))^j and

g(m,n,t) = g(m,n,s_m) for all t >= s_m

The N-s_m upper slices are therefore duplicates of one another. Sweeping t from s_m down instead of from N down keeps every value the answer reads.

The t = s_m slice is the only one that cannot be a shifted logaddexp of the slice above it, since it reads ITSELF at n-1; it is accumulated in a scalar loop, which is O(Nval) rather than the O(Nval) vector operations the original spends on the collapsed slices. Cost drops from O(M*Nval^2) to O(Nval*sum_k s_k), linear in the population on a multiserver model, and the log-domain branch remains a sum of nonnegative terms so no digits are lost to cancellation.

A station whose rates never settle, an infinite server alpha(n)=n being the usual case, gets s_k = Nval and costs what it costs in pfqn_gldsingle; the saving is over the other stations. An arbitrary rate matrix is accepted and simply yields s_k = Nval throughout, at which point this is pfqn_gldsingle.

Parameters:

L (ndarray) – Service demands at all stations (M x 1)
N (ndarray) – Number of jobs (scalar or 1x1 array)
mu (ndarray) – Load-dependent scaling factors (M x Ntot)
options (Dict[str, Any] | None) – Solver options (unused, for API compatibility)

Returns:

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

Raises:

RuntimeError – If multiclass model is detected

Return type:

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

Normalizing constant of a multiclass LIMITED load-dependent closed model.

Same recursion, same arithmetic and the same result as pfqn_gld, but with the rate shift saturated at the limited load-dependence threshold, which makes the recursion’s state space finite and lets it be memoised. This is what pfqn_lldsingle does to pfqn_gldsingle, one level up: there the rate offset is an index into a table, here it is the shift pfqn_mushift applies.

pfqn_gld peels the last station and advances its rate lattice one job at a time,

g(m,n,j) = g(m-1,n,0) + sum_r L[m,r]/alpha_m(j+1) * g(m,n-e_r,j+1)

with j the number of shifts row m has taken, so that pfqn_mushift’s leading element is alpha_m(j+1). Once j >= s_m-1, where s_m is the population past which alpha_m stays constant, every remaining entry of the row is alpha_m(s_m) and a further shift LEAVES THE ROW UNCHANGED over the columns the recursion can still read. Saturating j at s_m-1 therefore returns the same value and makes the state (m, n, j) repeat, at which point one memo answers what pfqn_gld recomputes down an exponential tree.

COST. The state space is M * prod_r(N_r+1) * max_k s_k, against pfqn_gld’s unmemoised recursion, which revisits the same states exponentially often. Without the saturation a memo would still be bounded, but by M * prod_r(N_r+1) * (Ntot+1): the threshold is what replaces the population by the server count, exactly as in pfqn_lldsingle.

Every terminal case of pfqn_gld is delegated back to it on the materialised block, so the two agree to the last bit rather than to a tolerance. A SYMBOLIC rate matrix is passed straight through to pfqn_gld: locating the threshold means comparing rates, which has no truth value on a symbol.

Parameters:

L (ndarray) – Service demands at all stations (M x R)
N (ndarray) – Number of jobs for each class (1 x R)
mu (ndarray) – Load-dependent scalings (M x Ntot)
options (Dict[str, Any] | None) – Solver options

Returns:

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

Return type:

pfqn_xia(L, N, s)[source]

Xia’s asymptotic approximation of the load-dependent normalizing constant.

The demands are first rescaled so that the largest per-server utilization rho_i = L_i/s_i is one. The stations that attain it are the bottleneck set B; they saturate and contribute the M/M/s saturated term, while every other station contributes its finite-capacity Erlang-like partial sum F(u,k) = sum_{j<k} u^j/j! + (u^k/k!)/(1 - u/k), the closed form of the geometric tail beyond the k-th server. The result is

log G ~ -log((|B|-1)!) - N log(c) + sum_{b in B} [s_b log L_b - log(s_b!)]

sum_{k not in B} log F(L_k, s_k),

with c the rescaling factor. The leading behaviour in N enters ONLY through -N log(c): this is the large-population limit, so the approximation does not resolve the O(1) corrections a finite population carries.

A non-bottleneck station with u > k gives a NEGATIVE F, whose logarithm is not real. Only an infinite F (u == k exactly) is dropped, matching the reference: suppressing a negative term would quietly return a plausible number for a model the expansion does not cover. The condition cannot arise when every station has one server.

Parameters:

L (ndarray) – Service demand vector (M,).
N (int) – Closed population (scalar).
s (ndarray) – Server counts (M,).

Returns:

Logarithm of the approximate normalizing constant.

Return type:

pfqn_mushift(mu, k)[source]

Shift a load-dependent scaling vector by one position.

Used in recursive normalizing constant computations.

Parameters:

mu (ndarray) – Load-dependent scalings matrix (M x N)
k (int) – Row index to shift

Returns:

Shifted mu matrix (M x N-1)

Return type:

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

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

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

Parameters:

L (ndarray) – Service demands at all stations (M x R)
N (ndarray) – Number of jobs for each class (1 x R)
Z (ndarray) – Think times for each class (1 x R)
mu (ndarray) – Load-dependent scalings (M x Ntot)
options (Dict[str, Any] | None) – Solver options

Returns:

PfqnComomrmLdResult with G, lG, and marginal probabilities

Return type:

PfqnComomrmLdResult

pfqn_fnc(alpha, c=None)[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 (ndarray) – Load-dependent scalings (M x N)
c (ndarray | None) – Scaling constants (1 x M), optional. If None, auto-selected.

Returns:

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

Return type:

PfqnFncResult

pfqn_ncoi(Z, N, mu=None, visits=None, options=None)[source]

Normalizing constant of a closed OI + single-delay product-form network.

The OI stations are analyzed by the balanced-fairness recursion of Bonald and Proutiere (2003) combined with the multichain convolution over stations. For a single OI station with rank rate mu(supp(n)) the balance function is Phi(0) = 1, Phi(n) = (1/mu(n)) sum_{r: n_r>0} Phi(n - e_r), and G(N) is the convolution of the per-station balance functions with the multinomial delay factor F_Z(n) = prod_r Z_r^{n_r}/n_r!, g_0 = F_Z, g_m(n) = sum_{0<=x<=n} Phi_m(x) g_{m-1}(n-x), G(N) = g_M(N).

This is a MACROSTATE routine: everything is tabulated over the count lattice 0 <= n <= N, never over orderings, which is legitimate because an OI rate is permutation-invariant so Phi closes on the count vector. With a non-empty swap graph that closure fails; use pfqn_pas_nc() instead.

Cost: O(M R L) for the balance functions and O(M prod_r (N_r+1)(N_r+2)/2) for the convolutions, with L = prod_r (N_r+1).

Parameters:

Z ((R,) think-time demand vector of the aggregated delay node.)
N ((R,) closed population vector, finite.)
mu (list of callables, one per OI station. Each ``mu[m](n)` returns the`) – total service rate of station m given the per-class occupancy (count) vector n. A station state with non-positive rate is unreachable and is assigned a zero balance value. May be empty/None for a pure delay network.
options (accepted for signature parity; unused.)

Returns:

(G, lG, Gtab) – lattice of normalizing constants, Gtab[dot(n, strides)] = G(n) for every 0 <= n <= N with strides = cumprod([1, N[:-1]+1]). The convolution produces this table anyway, so a caller needing G at more than one population must index this output, NOT re-call the routine per population – the latter costs a needless factor prod_r (N_r+1).

Return type:

normalizing constant, its natural log, and the WHOLE

pfqn_ncjd(Z, N, mu=None, visits=None, options=None)[source]

Joint-dependent name of pfqn_ncoi().

The two names denote the SAME routine because the balanced-fairness recursion mu_i(n) Phi_i(n) = sum_{r: n_r>0} v_{i,r} Phi_i(n - e_r) never inspects the structure of mu_i: it evaluates the handle at the full count vector n. Order independence (mu_i constant on each support) is a modelling restriction that buys insensitivity and a physical reading of Phi, not something the convolution uses, so any joint-dependent scaling eta_i(n) is admissible here.

Use pfqn_ncoi() when the model is genuinely order independent and the name should say so; use this one when the rate is a general joint dependence. See pfqn_clwjd() for the transform route, which needs the rate to saturate at a finite cutoff and is not merely a renaming.

Arguments and returns are exactly those of pfqn_ncoi().

pfqn_clwoi(Z, N, mu=None, visits=None, options=None)[source]

Normalizing constant of a closed delay + order-independent network.

Inverts the multichain generating function with the lattice-Poisson algorithm of Choudhury, Leung and Whitt (J. ACM 42(5):935-970, 1995). This is the transform counterpart of the convolution routine pfqn_ncoi(); both return the same G(N) and differ in cost.

An OI station factor is rational and available in closed form: splitting the count lattice by support S, on which mu_i(n) = mu_{i,S} is constant,

(mu_{i,S} - sum_{r in S} v_{i,r} z_r) F_{i,S}(z)

= sum_{r in S} v_{i,r} z_r F_{i,S-r}(z), F_{i,{}} = 1,

F_i(z) = sum_S F_{i,S}(z).

The singularities are the |S| hyperplanes sum_{r in S} v_{i,r} z_r = mu_{i,S}, one per support, and the restrictive static scaling of CLW eqs. 5.41-5.46 runs on the expanded constraint matrix that lists one row per (station, nonempty support) pair.

Cost: prod_r 2 l_r N_r contour points, each O(M R 2^R), against O(M prod_r (N_r+1)(N_r+2)/2) for the convolution of pfqn_ncoi(). The inversion is linear rather than quadratic in each population and returns G at the single population N.

Parameters:

Z ((R,) think-time demand vector of the aggregated delay node.)
N ((R,) closed population vector, finite.)
mu (list of callables, one per OI station. Each ``mu[m](n)` returns the`) – total service rate of station m at the per-class occupancy vector n and must depend on n only through its support. May be None for a pure delay network.
visits ((M x R) array or list of (R,) vectors of class visit ratios) – weighting the balance recursion. Default: unit visits.
options (dict with optional keys ``l` (inner lattice parameters) and`) – gamma (aliasing parameters). Defaults follow CLW.

Returns:

(G, lG)

Return type:

normalizing constant G(N), inf on overflow, and its natural log.

Raises:

ValueError – If a rate varies inside a support. Every handle is verified exhaustively on the count lattice before the inversion, because a violation would otherwise return a plausible but wrong G(N).

pfqn_clwjd(Z, N, mu=None, visits=None, lcut=None, options=None)[source]

Normalizing constant of a closed delay + limited joint-dependent network.

Joint-dependent generalization of pfqn_clwoi(), which is the case lcut = 1. Station i has a rate that reads the whole per-class occupancy but saturates coordinatewise: with a cutoff vector l_i,

mu_i(n) = c_{i,t}, t = (min(n_1, l_{i,1}), …, min(n_R, l_{i,R})),

so the clipped vector t ranges over a finite box and the station factor is rational, with

(mu_{i,t} - sum_{r: t_r = l_{i,r}} v_{i,r} z_r) F_{i,t}(z)
= sum_{r: t_r >= 1} v_{i,r} z_r F_{i,t-e_r}(z), F_{i,0} = 1.

The singular hyperplanes are indexed by the SATURATED sets, at most 2^R per station however large the cutoffs are. G(N) is then recovered by R nested lattice-Poisson inversions exactly as in pfqn_clwoi().

Cost: prod_r 2 l_r N_r contour points, each O(M R prod_r (lcut_{i,r}+1)). With lcut = N the region box is the whole lattice and the convolution of pfqn_ncjd() wins outright; the inversion pays off when the joint dependence saturates early.

Parameters:

Z ((R,) think-time demand vector of the aggregated delay node.)
N ((R,) closed population vector, finite.)
mu (list of callables, one per LJD station, each taking the per-class) – occupancy vector. May be None for a pure delay network.
visits ((M x R) array or list of (R,) vectors of class visit ratios.)
lcut ((M x R) matrix of per-station per-class saturation cutoffs >= 1, or a) – scalar/row broadcast to every station. Entries are clipped to N, which is exact. Default: N (no truncation).
options (dict with optional keys ``l`` and gamma.)

Returns:

(G, lG)

Return type:

normalizing constant G(N), inf on overflow, and its natural log.

pfqn_oi_fnc(Phi, N=None, f=None, options=None)[source]

OI generalization of the load-dependent functional server.

Builds an auxiliary OI station whose balance function Psi satisfies the convolution identity (Psi * Phi)(n) = (1 + f(n)) Phi(n), then inverts Psi to the FNC rate mu_f(n) = (sum_{r: n_r>0} Psi(n-e_r)) / Psi(n).

Parameters:

Phi (balance function of the existing OI station over the lattice, an) – R-dimensional array of shape (N_1+1, …, N_R+1), or a flat column-major vector.
N ((R,) closed population vector. Optional when Phi is a full) – R-dimensional array (then N = shape(Phi) - 1).
f (target queue-dependent function f(n), f(0)=0 (default f = sum(n)).)
options (accepted for signature parity; unused.)

Returns:

(muf, Psi, mu) – balance array (lattice shape), and tabulated FNC rate array.

Return type:

callable rate handle (Inf outside the lattice), FNC

pfqn_oi_insvc(oirate, N, options=None)[source]

Conditional mean number of in-service jobs per class at an OI station.

This is the quantity underlying the LINE utilization convention at order-independent stations, U_r = E[sir_r] / c with c the number of servers and sir_r the number of class-r jobs receiving a strictly positive service rate.

In an OI station the state is the ordered list c = (c_1,...,c_n) of job classes (position 1 = head) and the job in position p is served at the rank rate increment Delta_p(c) = mu(c_1..c_p) - mu(c_1..c_{p-1}), so the total rate telescopes to mu(c). Position p is in service when Delta_p(c) > 0, and sir_r(c) = #{p : c_p = r, Delta_p(c) > 0}. Note that sir_r counts JOBS, not servers: a single job served concurrently by several compatible servers counts once. This matches the definition used by the exact CTMC solver (State.to_marginal, PAS branch) and by LDES.

Because mu is permutation-invariant, the unnormalized weight of an ordering c of the multiset n factorizes over its prefixes as w(c) = prod_p 1/mu(n(c_1..c_p)), and Phi(n) = sum_c w(c) obeys the balanced-fairness recursion (condition on the tail element):

Phi(0) = 1,   Phi(n) = (1/mu(n)) sum_{r: n_r>0} Phi(n - e_r)

Conditioning the same way and using sir_r(c) = sir_r(c_1..c_{|n|-1}) + [c_{|n|} = r] * 1{mu(n) > mu(n - e_r)} gives the companion recursion for Xi_r(n) = sum_c w(c) sir_r(c):

Xi_r(0) = 0
Xi_r(n) = (1/mu(n)) [ sum_{s: n_s>0} Xi_r(n - e_s)
                      + 1{n_r > 0} 1{mu(n) > mu(n - e_r)} Phi(n - e_r) ]

Given n every ordering carries the same class-weight factor, so the conditional law of the ordering is w(c)/Phi(n) and E[sir_r | n] = Xi_r(n)/Phi(n) =: g_r(n), a function of the count vector alone. The station mean then follows from the count marginal pM as E[sir_r] = sum_n pM(n) g_r(n), or in normalizing-constant form from the functional-server identity of pfqn_oi_fnc() applied to f(n) = g_r(n) (note g_r(0) = 0, as required).

Parameters:

oirate (function ``mu(n)` returning the OI total service rate for the`) – per-class count vector n (length R). mu(0) is taken as 0.
N ((R,) closed population vector, finite.)
options (accepted for signature parity, currently unused.)

Returns:

g ((prod(N+1), R) table, column-major over the lattice 0 <= n <= N, with) – g[1 + sum(n * stride), r] = E[sir_r | n].
Xi ((prod(N+1), R) table with the sir-weighted balance ``Xi_r(n)`.`)
Phi ((prod(N+1),) table with the OI balance function ``Phi(n)`.`)

Return type:

Tuple[ndarray, ndarray, ndarray]

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

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

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

Parameters:

N ((R,) closed population vector (macrostate), finite.)
mu (list of exactly two callables. mu[m](n) is the total OI rank rate of) – station m given the per-class occupancy (count) vector n (permutation-invariant, but NOT a function of supp(n) alone); this is the count-based svcRateFun of an OI/PAS node.
H ((R, R) swap-graph adjacency (H[b, a] != 0 forbids a before b). Empty ->) – pure OI (all orderings feasible).
options (dict or options object; fields ``samples` (default 1e4),`) – seed (optional), verbose (default False), qlen (default True). qlen=False estimates ONLY the normalizing constant: the per-class prefix counts and their coefficients are neither accumulated nor allocated, and Q comes back as zeros. The ordering is drawn from the same stream either way, so G is unchanged to the last bit – this is for the callers that want G(N - e_r) and discard the rest.

Returns:

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

Return type:

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

pfqn_pas_nc(Z, N, mu=None, prec=None, options=None)[source]

Normalizing constant G_C of one communicating class of a closed P&S network.

With a non-empty swap graph the ordered-state chain is reducible (Comte and Dorsman, 2021, arXiv:2009.12299): the recurrent communicating classes are the placement-order-adhering sets and the product form pi(c) = prod_m Phi_m(c_m)/G_C holds per class. This routine returns G_C.

It is a MICROSTATE routine: it walks the ordered chains position by position, because with a placement order the reachable set is a set of ORDERINGS that does not collapse onto the count lattice. For the plain OI case (empty prec) the lattice does suffice and pfqn_ncoi() returns the same G at far lower cost.

Method: build station M’s chain head-first; appending class r at position k = sum(occ)+1 is admissible iff no class already placed at that station must come after r, and contributes 1/mu_M(occ+e_r); the chain may be finalized (recursing to station M-1) only when occ is a placement-order ideal at full multiplicity. Once every P&S station is peeled the residual population sits at the delay node with weight prod_r Z_r^{N_r}/N_r!.

Cost: one node per feasible ordered prefix; with an empty prec that is sum_{b<=N} C(|b|+M-1, M-1) |b|!/prod_r b_r!, factorial in sum(N). A placement order prunes the orderings, which is what makes the microstate walk affordable in the P&S case.

Parameters:

Z ((R,) think-time demand vector of the aggregated delay node; None/empty) – for a network with no delay station.
N ((R,) closed population vector, finite.)
mu (list of callables, one per P&S station. ``mu[m](n)` returns the total`) – service rate of station m for the per-class occupancy vector n.
prec (placement order. Either a list of M (R, R) precedence matrices, one) – per station, or a single (R, R) matrix broadcast to every station, with prec[m][i, j] != 0 iff class i must be placed before class j at station m (the closure returned by pas_placement(), fed by the global DAG of pas_swap2order()). None or all-zero means no order: every ordering is feasible and G is the plain OI constant. NOTE the orientation: around a cycle each downstream station traverses its chain in the opposite direction, so downstream stations take the TRANSPOSE of the upstream order ([P, P.T] for a two-station cycle). Passing the same P to both stations of a cycle silently returns a smaller, wrong G.
options (accepted for signature parity; unused.)

Returns:

(G, lG)

Return type:

normalizing constant G_C of the communicating class and its log.

pas_placement(H)[source]

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

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

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

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

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

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

Parameters:

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

Returns:

H

Return type:

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

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

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

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

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

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

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

Parameters:

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

Returns:

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

Return type:

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

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

Joint-dependent name of pfqn_mvaoi().

The two names denote the SAME routine because the recursion evaluates the rate handle at a full occupancy vector, mu_i(s_i + e_r) with s_i the shift already committed at the bottom of station i, and never inspects the structure of mu_i. This is the “third form” of the Conditional MVA of Casale, “A Note on Stable Flow-Equivalent Aggregation in Closed Networks” (QUESTA 2009): a rate depending on the full per-class occupancy vector.

Unlike the AMVA joint-dependence route, which evaluates eta at the MEAN arrival-instant vector 1 + E[Q] and therefore collapses a support indicator to 1, this routine evaluates the rate at exact integer occupancies and is exact for the balanced-fair station.

Arguments and returns are exactly those of pfqn_mvaoi().

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

Exact marginal load-dependent MVA for OI networks.

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

Parameters:

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

Returns:

(XN, QN)

Return type:

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

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

Bases: object

Result of normalizing constant computation.

G: float

lG: float

method: str = 'default'

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

Bases: object

Result of COMOM load-dependent computation.

G: float

lG: float

prob: ndarray

class PfqnFncResult(mu, c)[source]

Bases: object

Result of functional server scaling computation.

mu: ndarray

c: 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 (ndarray) – Service demand matrix (M x R)
mu (ndarray | None) – Load-dependent rate matrix (M x Ntot), optional - pass None if not used
gamma (ndarray | None) – Class-dependent service rate matrix (M x R), optional - pass None if not used
tol (float) – Tolerance for floating point comparison (default 1e-14)

Returns:

PfqnUniqueResult containing reduced matrices and mapping information

Return type:

PfqnUniqueResult

pfqn_expand(QN, UN, CN, mapping)[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 (ndarray) – Queue lengths from reduced model (M’ x R)
UN (ndarray) – Utilizations from reduced model (M’ x R)
CN (ndarray) – Cycle times from reduced model (M’ x R)
mapping (ndarray) – Mapping vector from pfqn_unique (length M), mapping[i] = unique station index

Returns:

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

Return type:

Tuple[ndarray, ndarray, ndarray]

pfqn_combine_mi(mi, mapping, M_unique)[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 (ndarray) – User-provided multiplicity vector (1 x M_original or M_original,)
mapping (ndarray) – Mapping vector from pfqn_unique (length M_original)
M_unique (int) – Number of unique stations

Returns:

Combined multiplicity vector (1 x M_unique)

Return type:

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.

Variables:

L_unique (numpy.ndarray) – Reduced demand matrix (M’ x R) with M’ <= M unique stations
mu_unique (numpy.ndarray | None) – Reduced load-dependent rates (M’ x Ntot), None if mu was empty
gamma_unique (numpy.ndarray | None) – Reduced class-dependent rates (M’ x R), None if gamma was empty
mi (numpy.ndarray) – Multiplicity vector (1 x M’), mi[j] = count of stations mapping to unique station j
mapping (numpy.ndarray) – Mapping vector (1 x M), mapping[i] = unique station index for original station i

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

L_unique: ndarray: Alias for field number 0

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

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

mi: ndarray: Alias for field number 3

mapping: ndarray: Alias for field number 4

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 (ndarray) – Queue population vector (M,)
lldscaling (ndarray | None) – Load-dependent scaling matrix (M x Nmax), optional
nservers (ndarray | None) – Number of servers per station (M,), optional

Returns:

Scaling factor vector (M,)

Return type:

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:

References

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

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

Sanitize and preprocess network parameters for NC solvers.

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

Parameters:

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

Returns:

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

Return type:

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

References

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

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

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

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

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

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

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

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

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

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

Parameters:

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

Returns:

Scaling factor vector (M,)

Return type:

References

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

pfqn_sdrcoeff(sdr)[source]

Validate an SDR structure and return its derived coefficients.

The structure carries, with branch index 1 reserved for the complement M-V:

entry, departure center indices of e and d of Q(V,V) branch[b] center indices of branch b, b >= 2; branch[0] unused entryOf[b] center index of the branch entry e(b) departureOf[b] center index of the branch departure d(b) level[b] the unique t with B_b in V_t - V_{t+1} C length T, the coefficients C_t of eq. (11) d T x B, d[t][b] read for 2 <= b <= B, 1 <= t <= level[b]

The population bounds mmax and vmax are consequences of the coefficients, not independent inputs: with C_t < 0 the routing enforces them by itself.

pfqn_sdrprob(sdr, n)[source]

SDR routing probabilities of eq. (10).

Returns (P, Ped) where P[b] is the probability of proceeding from the entry center e of Q(V,V) to the entry center of branch b (P[0] is zero, branch index 1 being the complement M-V) and Ped = 1 - sum(P) is the probability of proceeding directly to the departure center d, that is of being denied entry into Q(V,V) and returned to e.

These probabilities are chain independent: they read the total branch and subnetwork populations, not the per-chain ones. The chain-dependent form of eq. (1) has no published product form and is not implemented.

A branch population with delta_tb(m_b) < 0 lies beyond the bound that SDR enforces itself and is unreachable, so the probability there is zero.

pfqn_sdr(S, xi, N, sdr, alpha=None)[source]

Exact product form of eq. (16).

S and xi are M x J: the mean service times 1/mu_ij and the coefficients xi_ij of Section 3.2. They are required separately rather than as their product because under SDR the xi are not visit ratios, so the per-center throughputs cannot be recovered from the demands alone.

alpha is an optional M x sum(N) matrix of load-dependent rate scalings, alpha[i, k-1] = alpha_i(k). Defaults to a fixed-rate center; use k for an infinite server and min(k, c) for a c-server center.

Returns (Q, X, U, R, G, lG, prob, states), all M x J except the last three. X holds the per-center chain throughputs, U = X * S the mean number in service, and R = Q / X the response time at the center.

pfqn_sdrvisits(sdr, P)[source]

Coefficients xi of Section 3.2.

P is M x M x J: P[x, y, j] is the state-independent probability that a chain j customer leaving center x proceeds to center y. The state-dependent arcs out of the entry center are not part of P and are ignored if present.

Three rules fix the coefficients: the complement M-V obeys the ordinary traffic equations with the whole SDR subnetwork collapsed into a single e -> d arc of probability one; every branch obeys its own traffic equations driven by an injection of xi_e at its entry center; and xi_e = 1.

The paper states xi_ij = xi_ej for the branch entry and departure centers and works out only single-center branches. The traffic equations above are the reading that extends it: they return xi_{d(b)} = xi_e because a customer leaves a branch only through d(b), and xi_{e(b)} = xi_e whenever e(b) takes no internal feedback. They have been checked against a brute-force CTMC on a branch that does take such feedback, where the paper’s literal rule fails.

These xi are not relative visit counts: the rate at which customers enter a branch is state dependent, so a ratio of two xi carries no flow meaning.

pfqn_sdrmva(S, xi, N, sdr, alpha=None)[source]

Section 4 MVA and convolution of a network with state-dependent routing.

Krzesinski (1987), Performance Evaluation 7:125-143, Section 4. Same signature and outputs as pfqn_sdr, which evaluates eq. (16) exactly by state enumeration, so the two are directly comparable. This routine costs O(J T M (V_1…V_J)^2) rather than the size of the state space.

Restrictions, both from the paper: every SDR branch must be a SINGLE centre, and every C_t must be negative. A C_t other than -1 is rescaled internally, which leaves eqs. (10) and (16) unchanged because the factors telescope.

Returns (Q, X, U, R, lG).

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 (ndarray) – Array of non-negative numbers
Returns:: Array of log(n!) values
Return type:: 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:

oner(n, s)[source]

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

Parameters:

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

Returns:

Copy of n with n[s] -= 1

Return type:

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:

multichoosecon(n, S)[source]

Pick vectors of S elements from the available units in vector n.

Twin of MATLAB multichoosecon.m. Unlike multichoose, the count drawn from category i is capped by n[i], so the enumeration never proposes a job of a class the station does not hold.

Parameters:

n (ndarray) – Per-category availability
S (int) – Number of units to draw

Returns:

Matrix (C x len(n)) of draws, one per row; empty when S exceeds sum(n)

Return type:

matchrow(matrix, row)[source]

Find the index of a row in a matrix.

Parameters:

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

Returns:

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

Return type:

int

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

CoMoM algorithm for computing the normalizing constant.

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

Parameters:

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

Returns:

Logarithm of the normalizing constant

Return type:

lG

References

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

pfqn_comomrm(L, N, Z, m=1, atol=1e-8)[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 (ndarray) – Service demand matrix (1 x R) - single station
N (ndarray) – Population vector (R,)
Z (ndarray) – Think time vector (R,)
m (int) – Replication factor (default: 1)
atol (float) – Absolute tolerance

Returns:

ComomResult with lG (log normalizing constant) and lGbasis

Return type:

ComomResult

References

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

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

Original CoMoM implementation for repairman model.

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

Parameters:

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

Returns:

ComomResult with lG and lGbasis

Return type:

ComomResult

References

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

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

CoMoM for multi-server repairman model.

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

Parameters:

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

Returns:

ComomResult with lG and lGbasis

Return type:

ComomResult

References

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

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

ProCoMoM algorithm for computing marginal queue-length probabilities.

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

Parameters:

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

Returns:

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

Q: Mean queue length vector (M,).

Return type:

Pr

References

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

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

Projected CoMoM method for normalizing constant.

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

Parameters:

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

Returns:

Logarithm of the normalizing constant

Return type:

lG

References

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

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

Bases: NamedTuple

Result of CoMoM normalizing constant computation.

Create new instance of ComomResult(lG, lGbasis)

lG: float: Alias for field number 0

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

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

McKenna-Mitra integral form using scipy.integrate.

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

Parameters:

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

Returns:

G: Normalizing constant
lG: Logarithm of normalizing constant

Return type:

Tuple of (G, lG)

References

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

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

McKenna-Mitra integral with Gauss-Legendre quadrature.

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

Parameters:

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

Returns:

G: Normalizing constant
lG: Logarithm of normalizing constant

Return type:

Tuple of (G, lG)

References

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

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

McKenna-Mitra integral with Gauss-Laguerre quadrature.

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

Parameters:

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

Returns:

G: Normalizing constant
lG: Logarithm of normalizing constant

Return type:

Tuple of (G, lG)

References

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

pfqn_mmsample2(L, N, Z, m=1, n_samples=10000)[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 (ndarray) – Service demand vector (R,)
N (ndarray) – Population vector (R,)
Z (ndarray) – Think time vector (R,)
m (int) – Replication factor (default: 1)
n_samples (int) – Number of Monte Carlo samples (default: 10000)

Returns:

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

Return type:

Tuple of (G, lG)

References

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

logsumexp(x)[source]

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

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

pfqn_schmidt(D, N, S, sched, v=None)[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 (ndarray) – Service demand matrix (M x R)
N (ndarray) – Population vector (R,)
S (ndarray) – Number of servers per station (M,) or (M x R)
sched (ndarray) – Scheduling discipline per station (M,) - SchedStrategy values
v (ndarray | None) – Visit ratio matrix (M x R), optional (default: ones)

Returns:

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

Return type:

Tuple of (XN, QN, UN, CN)

References

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

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

Extended Schmidt MVA algorithm with queue-aware alpha corrections.

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

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

Parameters:

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

Returns:

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

Return type:

Tuple of (XN, QN, UN, CN)

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

Bases: object

Result from Schmidt’s exact MVA.

XN: ndarray

QN: ndarray

UN: ndarray

CN: 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 (ndarray) – Current population vector or max population
Nmax (ndarray | None) – Maximum population per class (if incrementing)

Returns:

Next population vector, or array of -1 if done

Return type:

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

Hash population vector to linear index.

Parameters:

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

Returns:

Linear index (1-based for MATLAB compatibility)

Return type:

int

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

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

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

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

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

Parameters:

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

Returns:

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

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

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

References

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

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

RECAL (REcursive CALculation) method for normalizing constant.

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

Parameters:

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

Returns:

G: Normalizing constant lG: Logarithm of normalizing constant

Return type:

Tuple of (G, lG)

References

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

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

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

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

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

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

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

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

Parameters:

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

Returns:

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

Return type:

Tuple of (XN, QN, UN, CN)

References

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

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

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

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

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

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

tau_k(v,i)  = T_ik^-1 sum_{n=0}^{k-1} P^{k-1}_i(n,v+1_i)
                         * mu_i(n+v_i+1)/(n+v_i+1)                   (21)
tau_k(v,i)  = T_ik^-1,                                   i IS        (22)
L^k_{jk}(v) = theta_jk tau_k(v,j)^-1 / sum_m theta_mk tau_k(v,m)^-1  (23)
lambda^k_k(v) = tau_k(v,j(k)) L^k_{j(k)k}(v)                         (24)
P^k_j(n,v)  = L^k_{jk}(v) P^{k-1}_j(n-1,v+1_j)
              + sum_{m != j} L^k_{mk}(v) P^{k-1}_j(n,v+1_m)          (25)

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

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

c_i(k,v)      = sum_{n=0}^{k-1} P^{k-1}_i(n,v+1_i)
                    * mu_i(n+v_i+1)/(n+v_i+1)
L^k_{jk}(v)   = (a_jk/c_j) / sum_m (a_mk/c_m)
lambda^k_k(v) = 1 / sum_m (a_mk/c_m)

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

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

Parameters:

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

Returns:

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

XN: Per-class throughput at the reference station (R,)
QN: Per-class mean queue-length at the QLD centers (M x R)
UN: Utilization of each center (M,), 1 - P_j(0). PER-STATION, not
    per-class, as in pfqn_mvald and pfqn_dac: for a load-dependent
    center the per-class product XN[r]*L[j,r] of pfqn_mvac is NOT
    the utilization

CN: Per-class cycle time exclusive of think time (R,),
    N[r]/XN[r]-Z[r], as in pfqn_mvald and pfqn_dac. NOT the (M x R)
    per-station residence time of pfqn_mvac: the whole LD family
    reports a cycle time here

pij: Marginal queue-length probabilities (M x (sum(N)+1)),
     pij[j,n] = P(n jobs at center j)

Return type:

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

References

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

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

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

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

Parameters:

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

Returns:

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

Return type:

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

References

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

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

Compute effective capacity terms for MVALDMX solver.

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

Parameters:

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

Returns:

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

Return type:

Tuple of (EC, E, Eprime, Lo)

References

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

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

Load-dependent MVA for multiserver mixed networks.

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

Parameters:

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

Returns:

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

Return type:

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

References

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

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

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

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

Parameters:

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

Returns:

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

Return type:

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

References

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

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

Linearizer for mixed open/closed queueing networks.

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

Parameters:

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

Returns:

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

Return type:

QN

References

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

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

Conway (1989) multiserver Linearizer approximation for FCFS queues.

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

Parameters:

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

Returns:

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

Return type:

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

References

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

ljd_linearize(nvec, cutoffs)[source]

Convert per-class population vector to linearized index.

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

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

Parameters:

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

Returns:

1-based linearized index

Return type:

int

References

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

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

Helper function for infinite radius computation with logistic transformation.

Used in normalizing constant computation via integration methods.

Parameters:

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

Returns:

Evaluated function value for integration

Return type:

References

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

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

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

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

Parameters:

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

Returns:

Evaluated function value for integration

Return type:

References

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

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

Knessl-Tier asymptotic expansion for normalizing constant.

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

Parameters:

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

Returns:

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

Return type:

Tuple of (G, lG, X, Q)

References

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

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

Knessl-Tier expansion corrected for the Stirling remainder (BKT).

pfqn_kt extracts N from the generating function of G by steepest descent. On the demand-free integral the exact coefficient is [u^N] exp(Z u) = Z^N/N!, whereas the expansion returns N log Z - (N log N - N + log(2 pi N)/2), Stirling’s approximation of log(N!) in place of log(N!). So KT lies ABOVE the exact value by the remainder s(N) per Laplaced class direction, and BKT subtracts sum_r s(N_r). The remainder is evaluated exactly from gammaln: truncating it at 1/(12 N) loses an order of magnitude (on the 1562 models of Cas17 sec5.3.1 the median |error| is 0.083 nats for KT, 1.9e-4 for the truncation and 1.4e-5 for the exact remainder).

Only the classes pfqn_kt actually Laplaces are corrected: a class with no jobs is dropped by its recursion and a self-looping class (one nonzero demand and no think time) has its coefficient extracted exactly, so neither carries a remainder. The predicate here is pfqn_kt’s own. See _kb/03-api-layer.md.

Parameters:

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

Returns:

(Gn, lGn), the normalizing constant and its logarithm.

Return type:

References

Original MATLAB: matlab/src/api/pfqn/pfqn_bkt.m Knessl and Tier, IEEE Trans. Computers 41(4):480-488, 1992.

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

Birman-Kogan saddle point normalizing constant with bottleneck detection.

Stations that serve a single chain and appear only once (the paper’s dedicated single servers) stay outside the exponent as O(1) algebraic factors, so their poles may be crossed by the saddle point; Algorithm 1 detects those chains and pins their coordinate on the pole. The remaining stations are the paper’s large groups of identical stations.

Parameters:

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

Returns:

G, lG: normalizing constant and its logarithm X: chain throughputs (the saddle point coordinates) U: utilizations (M x R) A: chains whose dedicated station is not saturated (eq. 29) B: chains whose dedicated station is a bottleneck (eq. 30)

Return type:

Tuple of (G, lG, X, U, A, B)

pfqn_bkue(L, N, Z=0.0)[source]

Birman-Kogan uniform (van der Waerden) expansion for a single chain.

The plain saddle point loses accuracy once the saddle approaches the dominant pole of the integrand, which is the regime where the station holding that pole saturates. The uniform expansion keeps the pole and the saddle in one formula through the complementary error function.

Parameters:

L (ndarray) – Service demand vector (M,), single class
N (float) – Population (scalar)
Z (float) – Think time (scalar)

Returns:

Tuple of (G, lG)

Return type:

pfqn_bklc(L, N, Z=None, method='mva', tol=1e-10, maxiter=1000)[source]

Birman-Kogan load concealment algorithm (Algorithm 2).

Chain l is solved on its own with every station slowed by the residual capacity the other chains leave it, A_i = 1 - sum_{k != l} L(i,k) X_k. Sweeping the chains in Gauss-Seidel order and iterating to a fixed point is the load concealment algorithm.

Parameters:

L (ndarray) – Service demand matrix (M x R)
N (ndarray) – Population vector (R,)
Z (ndarray | None) – Think time vector (R,), optional
method (str) – single chain solver, ‘mva’ (default) or ‘ue’
tol (float) – convergence tolerance on the throughputs
maxiter (int) – maximum number of sweeps

Returns:

Tuple of (X, Q, U, it)

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

Returns:

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

Return type:

AbAmvaResult

pfqn_ab_core(K, M, population, nservers, sched_type, v, s, maxiter, D_frac, l_in, fcfs_schmidt=False, marginal_prob_method='ab')[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 (ndarray) – Population vector (K,)
nservers (ndarray) – Number of servers at each station (M,)
sched_type (ndarray) – Scheduling strategies for each station (M,)
v (ndarray) – Visit ratio matrix (M x K)
s (ndarray) – Service time matrix (M x K)
maxiter (int) – Maximum iterations
D_frac (ndarray) – Fractional changes matrix (M x K x K)
l_in (ndarray) – Initial queue length matrix (M x K)
fcfs_schmidt (bool) – Whether to use Schmidt formula for FCFS stations
marginal_prob_method (str) – Method for marginal probability (‘ab’ or ‘scat’)

Returns:

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

Return type:

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

References

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

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

Bases: object

Result from Akyildiz-Bolch AMVA algorithm.

QN: ndarray

UN: ndarray

RN: ndarray

CN: ndarray

XN: ndarray

totiter: int

pfqn_rd(L, N, Z, mu=None, options=None)[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 (ndarray) – Service demand matrix (M x R)
N (ndarray) – Population vector (1 x R)
Z (ndarray) – Think time vector (1 x R)
mu (ndarray | None) – Load-dependent rate matrix (M x sum(N)). If None, assumes load-independent rates (all 1.0).
options (RdOptions | None) – Solver options

Returns:

lGN: Logarithm of normalizing constant Cgamma: Gamma correction factor

Return type:

Tuple of (lGN, Cgamma) where

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

Bases: object

Options for RD algorithm.

tol: float = 1e-06

method: str = 'default'

class RdResult(lGN, Cgamma)[source]

Bases: object

Result from Reduction Heuristic algorithm.

lGN: float

Cgamma: float

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

Logarithm of the normalizing constant of a limited load-dependent model.

Parameters:

L (ndarray) – Service demand matrix (M x R)
N (ndarray) – Population vector (R,)
Z (ndarray | None) – Think time vector (R,) or matrix (D x R) - optional
alpha (ndarray | None) – Load-dependent rate matrix (M x sum(N))

Returns:

Logarithm of the normalizing constant

Return type:

lG

pfqn_nre_full(L, N, Z=None, alpha=None, options=None, vfix=None)[source]

The full form of the reference’s four outputs, named alike in the JAR and the C++ port.

Parameters:

L (ndarray) – Service demand matrix (M x R)
N (ndarray) – Population vector (R,)
Z (ndarray | None) – Think time vector (R,) or matrix (D x R) - optional
alpha (ndarray | None) – Load-dependent rate matrix (M x sum(N))
options (Dict[str, Any] | None) – Solver options
vfix (ndarray | None) – Tilt to use instead of solving the saddle-point equation, None for the standard estimator. Supplying the tilt obtained at a nearby population makes numerator and denominator of a ratio share one expansion point, which is the Tierney-Kadane arrangement.

Returns:

The constant, the saddlepoint term alone and the tilt actually used

Return type:

PfqnNreResult

class PfqnNreResult(lG, G, lGs, vsad=None)[source]

Bases: object

The reference’s [lG,G,lGs,vsad].

lG - lGs is the Edgeworth correction, so a caller wanting the plain saddlepoint estimate reads lGs rather than re-deriving it. vsad is None on the shortcut arms that never solve a saddle point, matching the reference’s empty [].

lG: float

G: float

lGs: float

vsad: ndarray | None = None

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

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

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

Parameters:

L (ndarray) – Service demand matrix (M x R)
N (ndarray) – Population vector (R,)
Z (ndarray) – Think time vector (R,) - optional
alpha (ndarray) – Load-dependent rate matrix (M x Ntot); omitted means the load-independent mu_i(n) = 1 at every station

Returns:

Logarithm of normalizing constant

Return type:

lG

References

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

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

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

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

Parameters:

L (ndarray) – Service demand matrix (M x R)
N (ndarray) – Population vector (R,)
Z (ndarray) – Think time vector (R,) - optional
alpha (ndarray) – Load-dependent rate matrix (M x Ntot); omitted means the load-independent mu_i(n) = 1 at every station

Returns:

Logarithm of normalizing constant

Return type:

lG

References

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

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

Laplace approximation for normalizing constant.

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

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

Parameters:

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

Returns:

Logarithm of normalizing constant approximation

Return type:

lG

References

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

laplaceapprox(h, x0, tol=1e-5)[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 (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 (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: