api.pfqn

pfqn_explicit_ld(L, N, mu, tol, method, maxloss)

Explicit closed-form normalizing constant of a multiclass limited load-dependent network.

Load-dependent counterpart of pfqn_explicit. It evaluates the same divided-difference form of Casale, “Accelerating Performance Inference over Closed Systems by Asymptotic Methods”, ACM SIGMETRICS 2017, Corollary 3.2,

G(N) = sum_{0<=t<=N} (-1)^(|N|-|t|)/(N_1!…N_R!) prod_r nchoosek(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, “Facilitating Load-Dependent Queueing Analysis Through Factorization”, 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: pfqn_divdiff_ld is the same outer sum carried over pfqn_gldsingle’s O(M|N|^2) recursion instead.

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: this removes the floating-point RANGE problem but not the cancellation. 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. The fourth output reports the decimal digits lost and a warning is raised once the loss exceeds what double precision carries.

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); default all ones.
tol – Relative tolerance declaring two scaled demands redundant, and the rate tail constant (default: eps).
method – ‘auto’ (default), ‘distinct’ to force Eq. (15), ‘repeated’ to force Eq. (16).
maxloss – 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:

lG – Logarithm of the normalizing constant. G: Normalizing constant. method: Expression actually used for g_sigma, ‘distinct’ (Eq. 15) or ‘repeated’ (Eq. 16). lossDigits: Decimal digits lost to cancellation.

pfqn_explicit(L, N, tol, method, maxloss)

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 nchoosek(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), induced demands PAIRWISE DISTINCT. Gordon’s partial-fraction formula g_t(|N|) = sum_k theta_k^(|N|+K-1) / prod_{i~=k}(theta_k-theta_i), with the convention 0/0 = 0. It is O(K) per term of the outer sum.

Eq. (16), induced demands REDUNDANT (repeated). With K’ distinct induced demands theta_j of multiplicity m_j, the general partial-fraction expansion is used instead,

g_t(|N|) = sum_j (-1)^(m_j-1) theta_j^(|N|+K-m_j) * sum_{r>=0, |r|=m_j-1} (-1)^r_j nchoosek(|N|+r_j,r_j) prod_{k~=j} nchoosek(m_k+r_k-1,r_k) theta_k^r_k / (theta_j-theta_k)^(m_k+r_k)

which reduces to Eq. (15) when every m_j is one. It costs sum_j nchoosek(K’+m_j-2,m_j-1) evaluations per term of the outer sum.

The choice between the two is automatic: the induced demands are scanned over every t of the outer sum and Eq. (16) is used as soon as two of them are closer than tol relative to the largest induced demand at that t, Eq. (15) otherwise. tol defaults to machine precision (eps).

Only single-server load-independent queues are admissible: infinite servers need the integral form of Corollary 3.4 and load-dependent rates need pfqn_explicit_ld, which keeps this closed form as its inner kernel, or pfqn_divdiff_ld, which generalizes the outer sum over a recursion instead.

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. The fourth output reports the decimal digits lost, and a warning is raised once the loss exceeds what double precision carries.

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, which is why the single-class route is the cheap one on large N.

Parameters:

L – Service demand matrix (KxR) of single-server load-independent queues.
N – Population vector (1xR).
tol – Relative tolerance declaring two induced demands redundant (default: eps).
method – ‘auto’ (default), ‘distinct’ to force Eq. (15), ‘repeated’ to force Eq. (16).
maxloss – 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:

lG – Logarithm of the normalizing constant. G: Normalizing constant. method: Expression actually used, ‘distinct’ (Eq. 15) or ‘repeated’ (Eq. 16). lossDigits: Decimal digits lost to cancellation.

pfqn_divdiff_ld(L, N, Z, mu, options)

Exact load-dependent normalizing constant by divided differences of a single-class constant.

Evaluates the multiclass load-dependent normalizing constant as the N-th order finite difference, in the class populations, of a SINGLE-CLASS load-dependent constant:

G(N) = sum_{0<=n<=N} (-1)^(|N|-|n|) / (N_1! … N_R!) * prod_r nchoosek(N_r,n_r) * Gld_n(|N|)

where |N| = sum_r N_r, |n| = sum_r n_r and Gld_n(|N|) is the normalizing constant of the single-class model over the same M stations, with the same load-dependent rates alpha_i(.) = mu(i,.), total population |N|, and aggregated demands rho_i(n) = sum_r n_r * L(i,r). Since nchoosek(N_r,n_r)/N_r! = 1/(n_r! (N_r-n_r)!), the coefficient is evaluated here in the equivalent factorial form, which needs no binomial.

The identity holds for arbitrary rate functions alpha_i(.), hence it covers load-independent (alpha=1), infinite-server (alpha(j)=j), multiserver (alpha(j)=min(j,s_i)) and limited load-dependent stations alike. Proof: expand rho_i(n)^k_i multinomially in Gld_n(|N|); the R-fold difference operator annihilates every monomial whose degree in n_r is below N_r, and since the total degree is |N| the only surviving monomial is prod_r n_r^N_r, whose coefficient is the multiclass constant times prod_r N_r!.

Cost is prod_r (N_r+1) evaluations of pfqn_gldsingle, i.e. O(M |N|^2) time each. The O(1) space of the theoretical statement is attained only where the single-class constant itself has a closed form (e.g. Gordon’s formula in the multiserver case); this implementation uses the standard recursion instead, so the space is that of pfqn_gldsingle.

NUMERICS. The sum alternates in sign and its terms are much larger than the result, so it is evaluated as a signed log-sum-exp: this removes the floating-point RANGE problem but not the cancellation. The third output reports the decimal digits lost to cancellation and a warning is raised once the loss exceeds what double precision carries; multiprecision arithmetic is needed beyond that point.

Parameters:

L – Service demand matrix (MxR).
N – Population vector (1xR).
Z – Think time vector (1xR or DxR), folded in as an infinite server.
mu – Load-dependent rate matrix (Mx sum(N)), alpha_i(j) = mu(i,j).
options – Solver options.

Returns:

lG – Logarithm of the normalizing constant. G: Normalizing constant. lossDigits: Decimal digits lost to cancellation in the alternating sum.

pfqn_qdamva(L, N, Z, mu, Q0, tol, maxiter)

QD-AMVA, the queue-dependent AMVA of Casale-Perez-Wang (IFIP PERFORMANCE 2015), on a closed multiclass product-form network. Schweitzer/Bard core in which the class-r demand at station k is scaled by the queue-dependence term g_k evaluated at the arrival-instant total queue length, g = pfqn_lldfun(1 + delta*sum(Q(k,:)), mu).

SETTING MU TO A CONSTANT ROW RECOVERS PLAIN SCHWEITZER AMVA ONLY FOR A SINGLE CLASS. pfqn_lldfun does skip a constant row, so g == 1 there, but the residence time that remains is 1 + delta*sum(Q(k,:)) with ONE aggregate delta = (sum(N)-1)/sum(N) applied to the whole arrival-instant queue, where Bard- Schweitzer shrinks the TAGGED class alone: 1 + sum_{s~=r} Q(k,s) + (N(r)-1)/N(r) * Q(k,r). The two coincide iff K == 1. Measured over 40 random three-class instances, pfqn_qdamva(L,N,Z,ones) departs from pfqn_bs by up to 0.217 in absolute queue length, and is the LESS accurate of the two on single-server multiclass models (mean relative error on Q 0.069 against 0.056 at R = 3), the aggregate delta buying nothing once g == 1. This is the QD-AMVA closure, not a defect of the implementation, but do not use the function as a Schweitzer oracle for K > 1.

MU IS A DIMENSIONLESS RATE MULTIPLIER, NOT A RATE. mu(k,n) is the factor by which station k serves faster when it holds n jobs. Two traps follow from pfqn_lldfun: - it SKIPS a station whose mu row is constant (its range(…)>0 gate), so a single-server station must be ones(1,smax) and a c-server station min(1:smax, c). Passing a c-server station a constant row silently returns g=1, i.e. a single server. - smax = size(mu,2) must be at least ceil(sum(N)) or its interp1 clamps the population and the top of the rate curve is never reached.

Delay stations are carried in Z, not as rows of L. Closed classes only: an infinite N(r) is not supported.

Parameters:

L – (M x R) service demand matrix.
N – (1 x R) population vector, finite.
Z – (1 x R) think time vector.
mu – (M x smax) queue-dependent rate multipliers.
Q0 – (M x R) initial guess for the queue lengths.
tol – Convergence tolerance on the queue lengths (default 1e-6).
maxiter – Maximum number of iterations (default 1e4).

Returns:

Q – (M x R) mean queue lengths. X: (1 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.

[Q,X,U,ITER,R] = PFQN_QDAMVA(L,N,Z,MU,Q0,TOL,MAXITER)

pfqn_sdrvisits(sdr, P)

Section 3.2 coefficients xi of a network with state-dependent routing.

Parameters:

sdr – State-dependent routing structure.
P – State-independent routing probabilities per chain.

Returns:

xi – Coefficients xi_ij of the product form.

XI = PFQN_SDRVISITS(SDR, P)

Coefficients xi_ij of Krzesinski (1987), Section 3.2, for a network with state-dependent routing. P is an MxMxJ array of the state-independent (SIR) routing probabilities: P(x,y,j) is the probability that a chain j customer leaving center x proceeds to center y. The state-dependent arcs out of the entry center e of Q(V,V) are not part of P and are ignored if present.

The coefficients are fixed by three rules:

the entry center e and the departure center d of Q(V,V) satisfy xi_ej = xi_dj, and the centers of the complement M-V obey the ordinary traffic equations in which the whole SDR subnetwork acts as a single arc from e to d carrying probability one;
the entry center e(b) and the departure center d(b) of every branch, and every center inside a branch, obey the branch’s own traffic equations driven by an injection of xi_ej at e(b). Because a customer leaves a branch only through d(b), this returns xi_{d(b)j} = xi_ej, and it returns xi_{e(b)j} = xi_ej whenever e(b) receives no internal feedback, which covers every single-center branch. The paper states the identity 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 to a branch holding several centers;
the normalization xi_ej = 1.

These xi are not relative visit counts. Under SDR the rate at which customers enter a branch depends on the network state, so the ratio of two xi carries no flow interpretation; they are the solution of the transformed balance equations of Appendix A.

pfqn_sdrprob(sdr, n)

Krzesinski state-dependent routing probabilities.

Parameters:

sdr – State-dependent routing structure.
n – Per-station total population vector.

Returns:

P – Routing probabilities from the entry center to each branch entry. Ped: Probability of routing from the entry center to the departure center.

[P, PED] = PFQN_SDRPROB(SDR, N)

State-dependent routing (SDR) probabilities of Krzesinski (1987), “Multiclass Queueing Networks with State-Dependent Routing”, Performance Evaluation 7:125-143, eq. (10):

P_{e,e(b)}(N) = delta_tb(m_b) * prod_{s=1}^{t} omega_{s-1,s}(v_s)/omega_ss(v_s)

for the branch b at level t = level(b), and zero whenever omega_ss(v_s) = 0 for any s <= t. Here m_b is the total population of branch b, v_s the total population of the subnetwork Q(V_s,V_s), and

delta_tb(m) = C_t m + d_tb, omega_tt(v) = C_t v + D_tt, omega_{t-1,t}(v) = C_{t-1} v + D_{t-1,t}, omega_{0,1}(v) = 1.

N is the 1xM vector of total station populations. P is 1xB with P(b) the probability of proceeding from the entry center e of Q(V,V) to the entry center e(b) of branch b; P(1) is zero because branch index 1 denotes the complement M-V, which is not reached from e. 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 (Sec. 2.5).

These probabilities are chain independent: they are functions of the total branch and subnetwork populations, not of the per-chain populations. The chain-dependent form of eq. (1) has no published product form (the paper defers it to an unpublished IBM report) and is refused elsewhere.

A branch population m_b with delta_tb(m_b) < 0 lies beyond the bound that SDR itself enforces and so is unreachable; the probability returned there is zero, consistent with eq. (10) never routing a customer into such a state.

pfqn_sdrmva(S, xi, N, sdr, alpha)

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

Parameters:

S – Mean service times per station and chain.
xi – Relative visit counts per station and chain.
N – Chain population vector.
sdr – State-dependent routing structure.
alpha – Optional load-dependent rate scalings.

Returns:

Q – Mean queue lengths per station and chain. X: Mean throughputs per station and chain. U: Mean utilizations per station and chain. R: Mean response times per station and chain. lG: Logarithm of the normalizing constant.

[Q, X, U, R, LG] = PFQN_SDRMVA(S, XI, N, SDR, ALPHA)

Mean value analysis and convolution of a closed multiclass network with the state-dependent routing of Krzesinski (1987), “Multiclass Queueing Networks with State-Dependent Routing”, Performance Evaluation 7:125-143, Section 4.

The solution proceeds level by level, outermost subnetwork last:

Sec. 4.2.1 a modified MVA of the centers of V_t - V_{t+1} in isolation,: whose arrival theorem carries the SDR admission coefficient delta_ti(n) = d_ti - n and the center’s own queue-length distribution;
Sec. 4.2.2 convolution of that level with the already solved inner: subnetwork Q(V,V_{t+1}), weighted by Omega_{t-1,t}/Omega_tt;

Sec. 4.2.3 the same convolution re-normalizes the inner queue lengths; Sec. 4.3 a final convolution against the complement M-V.

Same signature and same outputs as PFQN_SDR, which evaluates the product form (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, and is the one to use when the populations make enumeration impractical.

RESTRICTIONS, both from the paper itself:

every SDR branch must be a SINGLE center. The paper’s Section 4 is stated that way and defers the general case to an unpublished technical report; PFQN_SDR has no such restriction.
every C_t must be negative. Section 2.5 assumes it, and Section 4 is written throughout for C_t = -1. A structure with C_t < 0 but not -1 is rescaled internally, which leaves both the routing probabilities of eq. (10) and the weights of eq. (16) unchanged: replacing (C_t, d_tb) by (C_t/k_t, d_tb/k_t) scales delta_tb and omega_tt by 1/k_t and omega_{t-1,t} by 1/k_{t-1}, and the factors telescope away.

pfqn_gerasimov(L, N, Z, options)

Gerasimov’s residue (closed-form) normalizing constant of a multiclass closed product-form network, generalized to R classes.

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

1. Gerasimov, “On Normalizing Constants in Multiclass Queueing Networks”, Operations Research 43(4):704-711, 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 (see the correspondence below).

Parameters:

L – Service demand matrix (M x R), L(i,r) = demand of class r at queueing station i.
N – Population vector (1 x R), nonnegative integers.
Z – Think time vector (1 x R). Default: zeros. A delay contributes the entire factor exp(sum_s Z_s u_s), which is handled exactly by convolving its Poisson coefficients into each elimination.
options – Struct with optional fields: .tol relative tolerance for declaring two affine forms

proportional, hence one pole rather than two. Default: 1e-12.

.maxterms cap on the number of residue terms carried between
eliminations. Default: 2e5. 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 – Normalizing constant. lG: Logarithm of the normalizing constant.

The recursion:

Write the same object 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,

and eliminate one class at a time. Every factor is AFFINE in u, f(u) = c_0 + sum_s c_s u_s, so when u_r is singled out it reads f = A - B u_r with A affine in the surviving variables. Partial fractions in u_r give, for a pole of order m_j at u_r = A_j/B_j,

[u_r^n] prod_j (A_j - B_j u_r)^-m_j

= sum_j sum_{k=0}^{m_j-1} (-B_j)^-k C(n+m_j-k-1,n) B_j^n
A_j^-(n+m_j-k) [t^k] prod_{l~=j} (C_jl - B_l t)^-m_l,

C_jl = (A_l B_j - B_l A_j)/B_j,

and C_jl, A_j are again AFFINE. The class-r elimination therefore maps a sum of products of affine powers into another one with one variable fewer: the structure is closed under the residue step, and R-1 steps leave a univariate coefficient extraction.

Why this is Gerasimov’s two-class formula:

At R = 2 the first (and only) step has all m_j = 1, so it returns one term per station i,

x_i2^(N_2+M-1) / prod_{k~=i}(x_i2-x_k2)

(1-x_i1 u_1)^-(N_2+1) prod_{k~=i} (1 - z_1ik u_1)^-1,

z_1ik = (x_k1 x_i2 - x_i1 x_k2)/(x_i2 - x_k2),

which is exactly the paper’s outer factor and its set Omega_i = {x_i1, z_1ij}. The pole of order N_2+1 at x_i1 is the paper’s tau_ik = taubar_ik + N_2, and coincidences inside Omega_i (the xi_i < M case of Thm 4) are the multiple poles the step already handles. The paper’s own hypotheses are relaxed in three places: x_i2 = x_k2 (his z_1ik is then undefined) leaves an affine form with no constant term, x_i2 = 0 leaves a factor with no pole in u_2, and identical station rows merge into one factor of doubled multiplicity. None of the three is a special case 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 O(M N), Buzen’s own cost; R = 2 O(M^2 N_(1)^2), INDEPENDENT OF N_(2); R >= 3 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 options.maxterms exists.

Conditioning:

The sum is an alternating one over residues, 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, whereas the eliminated ones cancel nothing; basing on N = 100 instead of N = 6 in one four-station model cost 39 nats of lG. What remains, 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.

pfqn_mva_interval(L, N, Z)

Exact interval-valued MVA for single-class closed product-form networks.

Exact output intervals of single-class MVA when the service demands, the think time and the population are known only up to intervals.

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

Reference: J. Luthi, G. Haring, “Mean value analysis for queueing network models with intervals as input parameters”, Performance Evaluation 32(3):185-215, 1998.

Parameters:

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

Returns:

X – Throughput interval (1 x 2). Q: Mean queue-length intervals per station (M x 2). U: Utilization enclosures per station (M x 2), exact where the demand

interval is thin, capped at 1 elsewhere.

R: Residence-time intervals per station (M x 2). Rtot: Total response-time interval (1 x 2). Qtot: Interval of the total number of jobs at the stations (1 x 2).

pfqn_ldbcmp(L, N, Z, c, varargin)

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

Lower bound on system throughput (and upper bound on response time) for a closed, single-class BCMP network with load-dependent stations, via the asymptotic closed-open equivalence of Anselmi and Cremonesi, “Bounding the Performance of BCMP Networks with Load-Dependent Stations” (2008). The bound (their eq. 15) exploits the monotonicity of system throughput and the fact that a closed BCMP network is, in the limit N -> inf, equivalent to the open network obtained by removing the bottleneck and injecting arrivals at rate 1/D_max. It is applicable when N >= Qhat and is asymptotically exact; Algorithm 1 refines it to a monotone fixed point.

Parameters:

L – Fixed-rate (limiting) service demand vector (M x 1). For a Heffes LD station, L(i) is the limiting demand D_i = lim_n D_i(n).
N – Total population (scalar).
Z – Think time (scalar; modeled as a non-bottleneck delay station).
c – Per-station Heffes load-dependence coefficient (M x 1, default 0). c(i)=0 marks a fixed-rate (LI) station with open queue rho_i/(1-rho_i); c(i)>0 a Heffes LD station with open queue (c(i)+1)*rho_i/(1-rho_i) (their eq. 22-23). The bottleneck is assumed fixed-rate (population transform (7) reduces to N’=N).
varargin – Optional trailing arguments, in order: tol, the fixed-point tolerance for Algorithm 1 (default 1e-10).

Returns:

Xlo – Lower bound on system throughput X(N); NaN if N < Qhat. Rhi: Upper bound on system response+think time, N/Xlo (Little). Qhat: Sum of non-bottleneck limiting queue lengths (eq. 11).

pfqn_harel_lb(rho, N, Z)

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

Parameters:

rho – Relative utilizations (k x 1), all strictly positive.
N – Closed population, at least 1.
Z – Think time; must be zero.

Returns:

LB – Throughput lower bound at population N.

LB = PFQN_HAREL_LB(RHO, N, Z)

Lower bound alone of PFQN_HAREL_BOUNDS,

LB = N / (A_1 + (N-1) (A_N/A_1)^{1/(N-1)}), 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.

See also PFQN_HAREL_BOUNDS, PFQN_HAREL_UB.

pfqn_cub_evals(M, order, Z)

NEVALS = PFQN_CUB_EVALS(M, ORDER, Z)

Number of integrand evaluations pfqn_cub performs on an M-station model at the given cubature order. The Grundmann-Moeller rule of degree ORDER on the (M-1)-simplex evaluates sum_{d=0..ORDER} nchoosek(M-1+2d, M-1) points, and a non-zero think time makes pfqn_cub repeat the whole rule at each of its v-quadrature steps. pfqn_nc prices CUB against this before selecting it.

pfqn_mvaoi_marg(D, N, isDelay, mu)

[XN, QN] = PFQN_MVAOI_MARG(D, N, ISDELAY, MU)

Exact marginal load-dependent MVA for a closed product-form network of infinite-server (delay) and load-independent (single-server, product-form) stations plus ANY number of order-independent (OI) stations. An OI station is a class-dependent load-dependent server whose total service rate mu_i(n) is a permutation-invariant function of the per-class count vector n. The network is product-form and is solved exactly by the load-dependent MVA recursion that carries, for EACH OI station i, its joint count-vector marginal distribution pM_i(n | k):

pM_i(n | k) = (1/mu_i(n)) * sum_r X_r(k) * pM_i(n - e_r | k - e_r), n ~= 0 pM_i(0 | k) = 1 - sum_{n ~= 0} pM_i(n | k)

All OI-station marginals share the common per-class throughput X_r(k); the recursion is exact per station because in product form pM_i(n|k) = Phi_i(n) G_{-i}(k-n)/G(k) with X_r(k) = G(k-e_r)/G(k) and the balanced-fairness identity Phi_i(n) = (1/mu_i(n)) sum_r Phi_i(n-e_r).

Because the OI rate is class-dependent, the mean-value response-time formula is not exact; instead the per-class throughput X_r(k) is closed at each population level by population conservation

X_r(k) * A_r(k) + sum_i QM_ir(k; X) = k_r, A_r(k) = sum_{i not OI} R_ir(k)

where QM_ir(k) = sum_n n_r pM_i(n | k) is read off each OI station’s exact marginal. This yields Q, X matching the exact CTMC / normalizing-constant (NC) results for any number of OI stations. This is the marginal-distribution counterpart of PFQN_MVAOI (the mean-value CMVA form).

Parameters:

D - (M x R) – Rows of OI stations are ignored (rate comes from MU).
N - (1 x R)
isDelay - (1 x M) logical, true for infinite-server (delay)
mu - (1 x M) cell; mu{i} is a function handle mu_i(n) – total service rate for the per-class occupancy vector n (1 x R) at OI station i, and [] for non-OI stations.

Returns:

XN - (1 x R) per-class throughput X_r = G(N-e_r)/G(N). QN - (M x R) per-class mean queue-length at every station.

Reference:: Reiser, Lavenberg (1980). Mean-Value Analysis of Closed Multichain Queuing Networks. JACM 27(2). Load-dependent extension: Bruell, Balbo, Afshari (1984). OI stations: Casale, Comte, Dorsman (2026).

pfqn_dmlin(L, N, Z, type, tol, maxiter, QN0, npasses)

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

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). Time drops to O(K C^2) with the space unchanged at O(K C^2), and, 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. That is what makes IL dominate AQL (pfqn_aql), which buys the same cost by aggregating the queue lengths themselves and does change the answer.

Parameters:

L – Service demand matrix (stations x classes).
N – Population vector.
Z – Think time vector.
type – Scheduling strategy type per station (accepted, unused: the Linearizer family in LINE treats every station as single-server PS).
tol – Convergence tolerance (default: 1e-8).
maxiter – Maximum number of iterations (default: 1000).
QN0 – (M x R) queue lengths that warm-start the Bard-Schweitzer initialization; empty for the default cold start.
npasses – Number of xi refresh passes (default 3, the Chandy-Neuse rule).

Returns:

Q – Mean queue lengths. U: Utilization. W: Residence times. C: Cycle times. X: System throughput. totiter: Total iterations performed.

pfqn_sens_respt(S, V, N, Z, b, tmax)

Exact moments of the sojourn time of a job at FCFS multiserver centers of a closed product-form queueing network.

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.

Reference: J. C. Strelen, “Moment Analysis for Closed Queuing Networks and its Linearizer”, Performance Evaluation 11:127-142, 1990, Theorem 4.1 with equations (4.1)-(4.5) and Remarks 4.2-4.3.

Parameters:

S – Service time at each station (M x 1), common to all classes.
V – Visit ratio matrix (M x R). The demand is L(i,r) = S(i)*V(i,r).
N – Population vector (1 x R).
Z – Think time vector (1 x R). Default: zeros.
b – Number of servers at each station (M x 1). Default: ones.
tmax – 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.

Returns:

res –

A struct with:

.X (1 x R), .Q (M x R), .U (M x R) base measures at population N. .W (M x 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 (M x R x tmax) WM(i,l,t) = E[W_(i,l)^t]. .WVar (M x R) Var[W_(i,l)] = E[W^2] - E[W]^2. Requires tmax >= 2. .WSkew (M x R) skewness of W_(i,l). Requires tmax >= 3; NaN if the

variance is zero.

.m (M x 1) E[Q_i] at population N, the total queue length. .Var (M x 1) Var[Q_i] at population N. .p (M x max(b)) p(i,1+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,1:b(i)).

.Wresid (M x R) the residence time w_i(l) of the MVA recursion. The
identity W(i,l) = w_i(l)/V(i,l) is an independent check of the t = 1 case of (4.5) and is asserted by pfqn_sens_respt_validate.

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

Validation harness for pfqn_sens_mva. Checks the MVA-like moment recursion of de Souza e Silva and Muntz (1988), Corollary 1, against three independent references:

brute-force enumeration of the closed product-form equilibrium distribution (ground truth, mi==1); B. the differentiated-MVA Jacobian of pfqn_sens, via the identity Cov[n(i,r),n(i,s)] = L(i,s) * dQ(i,r)/dL(i,s) (covers mi>1); C. pfqn_mva for the base measures X, Q, U, R.

It also reports the raw asymmetry of QCov before symmetrization: the recursion computes W(k,j;t,j) and W(t,j;k,j) by numerically distinct expressions, so their agreement is a nontrivial check of the formula.

pfqn_sens_mom_validate()

Validation harness for pfqn_sens_mom, the higher-moment analysis of Strelen (1990).

Five references: A. brute-force enumeration of the closed product-form distribution, which is ground truth for m, Var, Cov, E[Q^2] and E[Q^3]. B. pfqn_sens_mva. Summing its per-class covariance matrix at station i over all class pairs must give Var[Q_i], since Var[sum_r n(i,r)] = sum_{r,s} Cov[n(i,r),n(i,s)]. This ties the per-station-total moments of Strelen to the finer per-class moments of de Souza e Silva and Muntz. C. pfqn_mva for the base measures. D. Cov symmetry: x_j dm_i/dx_j and x_i dm_j/dx_i are computed by different derivative tracks and must agree. E. the published table of Example 3.4 of the reference (the Kobayashi central-server model), which pins the second derivative against numbers the author printed rather than against our own code.

Reference: J. C. Strelen, “Moment Analysis for Closed Queuing Networks and its Linearizer”, Performance Evaluation 11:127-142, 1990.

pfqn_sens_linearizer_validate()

Validation harness for pfqn_sens_linearizer, the LINEARIZER-2 / LINEARIZER-3 moment approximation of Strelen (1990), Section 5.

This routine is an APPROXIMATION, so it must not be held to machine precision. The checks are therefore of three kinds:

accuracy bands against the exact pfqn_sens_mom, on models small enough for the exact lattice. The reference reports relative errors below 2.1% on E[Q], 4.1% on E[Q^2] and 6.2% on E[Q^3] over its own 51 networks; the bands asserted here are of that order. This is the only meaningful statement of correctness for an approximation: it must track the exact answer, not equal it. B. exactness where the approximation degenerates. At a population of one job the CORE estimate of the queue lengths one job down is identically zero whatever the delta terms are, so the Linearizer equations coincide with the exact MVA and every moment must match pfqn_sens_mom to roundoff. This pins the derivative algebra independently of the heuristic. C. structural invariants that hold for any population: the mean queue lengths conserve the population, and the mean queue lengths agree with LINE’s own pfqn_linearizer, which runs the same heuristic without derivatives.

Reference: J. C. Strelen, “Moment Analysis for Closed Queuing Networks and its Linearizer”, Performance Evaluation 11:127-142, 1990.

pfqn_sens_ldmx_ec(lambda, D, mu)

Computes the effective capacity 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 the reference, which are reproduced here term by term.

Reference: I. F. Akyildiz and J. C. Strelen, “Moment Analysis for Load-Dependent Mixed Product Form Queueing Networks”, IEEE Trans. Communications 39(6):828-832, 1991.

Parameters:

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

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 x 1). dEC: dEC(i,n)/dLo(i), same shape as EC. dE: dE(i,1+n)/dLo(i), same shape as E. dEprime: dEprime(i,1+n)/dLo(i), same shape as Eprime.

[EC,E,EPRIME,LO,DEC,DE,DEPRIME] = PFQN_SENS_LDMX_EC(LAMBDA,D,MU)

pfqn_schmidt(D, N, S, sched, v)

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

Parameters:

D – Service demand matrix.
N – Population vector.
S – Number of servers per station (matrix or vector).
sched – Scheduling discipline per station.
v – Visit ratio matrix (optional, defaults to ones).

Returns:

XN – System throughput. QN: Mean queue lengths. UN: Utilization (M x R), per station-class, D*X/nservers. CN: Cycle times. T: Results table.

[XN,QN,UN,CN] = PFQN_SCHMIDT(D,N,S,SCHED,V)

pfqn_qlen_joint_moments(L, N, Z, pairs, route, lGsrc, options)

out = pfqn_qlen_joint_moments(L, N, Z, pairs, route, lGsrc, options)

Joint moments of the queue-length vector of a closed product-form network, obtained from normalizing constants.

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. Two exact routes reach the joint survival array, and both end in the same conversion, the tail edge of the house of moments (api/moment) followed by the joint central-moment and cumulant conversions:

SINGLE CLASS (R = 1), route ‘tail’. The survival probabilities are ratios of normalizing constants of the network itself,

P(n_i >= k_i for all i) = (prod_i L_i^k_i) * G(N - sum_i k_i) / G(N)

which holds because a load-independent single-class station has the geometric occupancy L_i^n. Only N+1 constants of the ORIGINAL model are needed, which is why any normalizing-constant algorithm serves it.

MULTICLASS, route ‘pmf’. The geometric factorization fails, since a multiclass load-independent station carries the multinomial occupancy f_i(n_i) = (|n_i|)! prod_r L_ir^n_ir / n_ir!. What holds instead is the joint law of the selected stations in terms of the COMPLEMENTARY network, the model with those stations deleted and the think times kept,

P(n_i = m_i, i in S) = prod_i f_i(m_i) * G_(S^c)(N - sum_i m_i) / G(N)

The survival array is the reverse cumulative sum of that array, exactly, since the box covers the support.

Neither the factorial nor the raw moments have a one-constant closed form; the survival array is the queue-length functional that does. The normalizing-constant algorithm is INJECTED rather than called at a fixed site: the whole set of populations is known before any evaluation, so it is emitted in one batch and an algorithm that produces several constants in one pass serves it without recomputation.

Input:

L: service demand matrix of the QUEUEING stations (MxR). Delay stations: belong in Z, their marginals following a different law. Load-dependent and multiserver stations are out of scope for both routes

N: population vector (1xR) Z: think time vector (1xR), zeros if empty pairs: Px2 matrix of 1-based (station,class) pairs, one per dimension of

the returned arrays. Defaults to every class of every station

route: ‘auto’ (default), ‘tail’ or ‘pmf’ lGsrc: where log G comes from. Empty calls pfqn_nc. A function handle is

invoked ONCE per network as lGsrc(Lsub, pops), pops being a PxR matrix of populations, and must return P values of log G with NaN where it cannot serve; those are filled in by pfqn_nc. A numeric array is read as a table indexed by population, which is what a convolution sweep produces for free. The ‘pmf’ route queries the COMPLEMENTARY network, so a table must be its table

options: options struct passed to pfqn_nc. Defaults to: SolverNC.defaultOptions with method ‘exact’, since an approximate normalizing constant would silently make the moments approximate

Output:

out: struct with the joint arrays over the selected coordinates (tail,: binomial, factorial, raw, central, cumulant), the mean vector, the covariance matrix cov, and info holding route, points, served, evals and the pairs used

Example

out = pfqn_qlen_joint_moments([2 1], [4 3], [0.5 0.8], [1 1; 1 2]); cov12 = out.cov(1,2);

Reference: M. Reiser and S. S. Lavenberg. Mean-value analysis of closed multichain queuing networks. Journal of the ACM, 27(2):313-322, 1980.

pfqn_procomom2(L, N, Z, mu, m)

Product-form CoMoM for 2-station repairman model (queue + delay).

Parameters:

L – Service demand vector.
N – Population vector.
Z – Think time vector.
mu – Load-dependent rates (optional).
m – Replication factor (default: 1).

Returns:

pk – Marginal state probabilities. lG: Logarithm of normalizing constant. G: Normalizing constant. T: Transfer matrices. F: Product transfer matrix. B: Combined transfer matrix.

Marginal state probabilities for the queue in a model consisting of a queueing station and a delay station only.

pfqn_pas_is(N, mu, H, options)

[G, LG, Q] = PFQN_PAS_IS(N, MU, H, OPTIONS)

Importance-sampling (IS) estimate of the normalizing constant of a SINGLE communicating class of a cyclic two-station pass-and-swap (P&S) queueing network with swap graph H. This is the Monte-Carlo counterpart of the exact microstate convolution PFQN_PAS_NC: it estimates the same per-communicating-class constant G_C but scales to populations where the exact enumeration of the feasible orderings becomes expensive.

Model. Two OI/P&S stations (1 = upstream, 2 = downstream of the cycle) hold all N jobs (no delay). With a non-empty swap graph the ordered-state chain is reducible; the recurrent communicating class is the set of splits of the orderings that are non-decreasing w.r.t. the placement partial order induced by H (Comte & Dorsman, 2021, arXiv:2009.12299). Writing D for that set of orderings and Phi_m for the balanced-fairness balance function of station m,

G_C = sum_{c in D} sum_{k=0}^{ell} Phi_1(c_{1..k}) Phi_2(c_{ell..k+1}),

where c_{1..k} is the length-k prefix placed at station 1 and c_{ell..k+1} the reversed suffix placed at station 2. Along a fixed ordering q the balanced- fairness value is the ordered product of reciprocal rank rates,

Phi_m(q) = prod_{p=1}^{|q|} 1 / mu_m(n(q_{1..p})), n(.) = prefix counts,

evaluated at the per-class COUNT vector of each prefix (OI property P1 makes mu permutation-invariant, i.e. a function of the counts – not of the support alone, which differs as soon as any class holds two or more jobs).

Auto-normalized IS (notebook generator IS_3). Orderings c are drawn from D by placing, at each step, a uniformly random placement-order-minimal present class; the draw probability p(c) is the product of the reciprocal branching factors. Then, for any coefficient xi,

G_C[xi] = E_{C~p}[ (sum_k xi(C,k) Phi_1(C_{1..k}) Phi_2(C_{ell..k+1})) / p(C) ],

and E[xi] = G_C[xi]/G_C[1] reuses the SAME samples for numerator and denominator (auto-normalized IS; A. Owen, MCM notes; convergence per Agapiou et al. 2017). Taking xi = number of class-r jobs in the prefix yields the mean queue length of class r at station 1.

Parameters:

N - (1 x R) closed population vector (the macrostate)
mu - cell {1 x 2} of function handles. mu{m} (n) – rank rate of station m given the per-class occupancy (count) vector n (1 x R); OI, so it depends only on supp(n), i.e. the sum of the capacities of the servers compatible with the present classes. This is exactly the svcRateFun stored on an OI/PAS node.
H - (R x R) – non-decreasing w.r.t. H: class a may not precede class b whenever H(b,a) ~= 0. The empty/all-zero graph reduces D to all orderings (pure OI); the estimate then targets the OI constant of PFQN_NCOI.
options - solver options (optional) – .samples number of IS samples (default 1e4); .seed RNG seed for reproducibility (optional); .verbose print progress (default false); .qlen estimate the queue lengths too (default true). False

estimates ONLY G: the prefix-count coefficients are neither allocated nor accumulated and Q comes back zero. 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 read nothing else from it.

Returns:

G - IS estimate of the communicating-class normalizing constant G_C. lG - log(G). Q - (2 x R) IS estimate of the mean per-class queue length; Q(1,:) at

station 1, Q(2,:) = N - Q(1,:) at station 2.

Example (two OI/P&S queues, R=5, star swap graph):: mu1 = [1 2 0.5]; nb1 = {1,2,3,[1 3],[2 3]}; mu2 = [1 2]; nb2 = {1,2,[1 2],1,2}; rate = @(nb,mu,n) sum(mu(unique([nb{n>0}]))); mu = {@(n) rate(nb1,mu1,n), @(n) rate(nb2,mu2,n)}; H = [0 0 0 1 0; 0 0 0 0 1; 0 0 0 1 1; 0 0 0 0 0; 0 0 0 0 0]; [G,lG,Q] = pfqn_pas_is([1 1 1 3 3], mu, H, struct(‘samples’,1e5));

See also PFQN_PAS_NC, PFQN_NCOI, PAS_PLACEMENT, PAS_SWAP2ORDER.

pfqn_panaceald(L, N, Z, mu, terms)

PANACEA asymptotic expansion for load-dependent closed networks.

Parameters:

L – Service demand matrix (MxR).
N – Population vector (1xR).
Z – Think time matrix (DxR), summed over rows.
mu – Load-dependent rate matrix (Mx sum(N)).
terms – Number of terms in the normal-usage asymptotic series (1, 2, or 3; default 3), as in pfqn_panacea.

Returns:

Gn – Normalizing constant. lGn: Logarithm of normalizing constant.

[GN,LGN]=PFQN_PANACEALD(L,N,Z,MU,TERMS)

pfqn_oi_fnc(Phi, N, f, options)

[MUF, PSI, MU] = PFQN_OI_FNC(PHI, N, F, OPTIONS)

Order-independent (OI) generalization of the load-dependent functional server (FNC) of Casale, “On Single-Class Load-Dependent Normalizing Constant Equations”, QEST 2006 (Theorem 3, Corollary 1). It is the OI counterpart of PFQN_FNC.

Given the balance function Phi of an existing OI station in the model and a queue-dependent target f(n) (default f(n) = sum(n), the total occupancy), the routine builds a functional server: an auxiliary OI station whose balance function Psi satisfies the convolution identity

(Psi * Phi)(n) = (1 + f(n)) Phi(n), (*)

where * is the OI (population-lattice) convolution (Psi * Phi)(n) = sum_{0<=k<=n} Psi(k) Phi(n-k). By construction Psi(0)=1 (since f(0)=0). Inserting this server (with the same per-class demand as the existing station) into the network and letting G, G^{+} be the OI normalizing constants without and with it, (*) yields the FNC identity

E[f(n)] = G^{+}/G - 1,

i.e. the mean of the queue-dependent function is read off a ratio of normalizing constants, with no probabilities and no Little’s law. For f(n)=sum(n) this returns the exact total mean queue length of the station.

Construction (two steps):

Deconvolution of (*) for the FNC balance function, solved triangularly in column-major order (k<n precedes n):

Psi(n) = (1+f(n)) Phi(n) - sum_{0<=k<n} Psi(k) Phi(n-k).
Balanced-fairness inversion of Psi to the FNC rate function:
mu_f(n) = ( sum_{r: n_r>0} Psi(n-e_r) ) / Psi(n).

Step 2 alone is the scalar PFQN_FNC when R=1 and Phi is the trivial LI balance (Phi==1): there Psi==1 and mu_f==1, i.e. the FNC degenerates to an identical LI copy of the station, as in the QEST 2006 example.

As noted in that reference (Sec. 6), the FNC balance/rate may be signed or non-physical; this is immaterial because only the final normalizing-constant ratio is used, and it stays positive with its usual interpretation. Hence the functional server is best convolved through its balance function Psi rather than through the positivity-guarded rate peeling of PFQN_NCOI.

Parameters:

Phi - balance function of the existing OI station over the lattice, an – R-dimensional array of size (N_1+1) x … x (N_R+1) (Phi(n) at subscript n+1), or a flat column-major vector. Obtain it by the forward balanced-fairness recursion from the station rate function.
N - (1 x R) – full R-dimensional array (then N = size(Phi) - 1).
f - target queue-dependent function handle f (n), n a (1 x R) – vector, with f(0)=0 (default f = @(n) sum(n)).
options - solver options (optional, accepted for signature parity)

Returns:

muf - function handle muf(n) giving the FNC rate mu_f(n) (out-of-lattice: states return Inf). mu_f(0)=0; states with Psi(n)=0 return Inf.

Psi - R-dimensional array (lattice shape) with the FNC balance function. mu - R-dimensional array with the tabulated FNC rate mu_f.

See also PFQN_FNC, PFQN_NCOI.

pfqn_mvacld(L, N, Z, mu)

MVAC (Mean Value Analysis by Chain) for load-dependent networks.

Parameters:

L – Service demand matrix of the queue-length dependent centers (M x R).
N – Population vector (1 x R).
Z – Think time matrix of the infinite-server centers (1 x R or Mz x R).
mu – Load-dependent rate matrix (M x Nt), mu(j,n) = rate with n jobs.

Returns:

XN – Per-class throughput (1 x R). QN: Per-class mean queue-length at the centers (M x R). UN: Utilization of each center (M x 1), i.e. 1-P_j(0). CN: Per-class cycle time exclusive of think time (1 x R). pij: Marginal queue-length probabilities (M x (Nt+1)).

[XN,QN,UN,CN,PIJ] = PFQN_MVACLD(L,N,Z,MU)

Exact mean value analysis by chain (MVAC) of a closed multichain product-form queueing network that may contain queue-length dependent (QLD) service centers. This is the extension of Section V of Conway, de Souza e Silva and Lavenberg, “Mean Value Analysis by Chain of Product Form Queueing Networks”, IEEE Trans. Computers 38(3):432-442, 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 extra information comes almost for free, and PIJ is returned as a first-class output.

Notation follows PFQN_MVAC: j and i index centers (j = 1,…,J1 are the QLD centers of L, j = J1+1,…,J the IS centers of Z), k and 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, the recursion of Section V reads

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 of the algorithm 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 - (M x R)
N - (1 x R)
Z - (1 x R) or (Mz x R) – center. Defaults to zeros(1,R).
mu - (M x Nt) load-dependent rates, Nt >= sum(N); mu(j,n) – service rate of center j with n jobs present. mu(j,:) = 1 is a single-server fixed-rate queue, mu(j,n) = min(n,c) a c-server queue, mu(j,n) = n an infinite server. Defaults to ones(M,sum(N)), i.e. all centers SSFR, in which case the results agree with PFQN_MVAC.

Returns:

XN - (1 x R) per-class throughput at the reference station. QN - (M x R) per-class mean queue-length at the QLD centers. UN - (M x 1) utilization of each center, 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 - (1 x R) per-class cycle time exclusive of think time, 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 - (M x (sum(N)+1)) marginal queue-length probabilities,: pij(j,n+1) = P(n jobs at center j).

See also PFQN_MVAC, PFQN_MVALD, PFQN_DAC, PFQN_GLD.

pfqn_mvac(L, N, Z)

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

Parameters:

L – Service demand matrix of the single-server fixed-rate queues (M x R).
N – Population vector (1 x R).
Z – Think time matrix of the infinite-server centers (1 x R or Mz x R).

Returns:

XN – Per-class throughput (1 x R). QN: Per-class mean queue-length at the queues (M x R). UN: Per-class utilization at the queues (M x R). CN: Per-class residence time at the queues (M x R).

[XN,QN,UN,CN] = PFQN_MVAC(L,N,Z)

Exact mean value analysis by chain (MVAC) of a closed multichain product-form queueing network composed of single-server fixed-rate (SSFR) queues and infinite-server (IS) centers, as given in Conway, de Souza e Silva and Lavenberg, “Mean Value Analysis by Chain of Product Form Queueing Networks”, IEEE Trans. Computers 38(3):432-442, 1989.

Unlike the classic MVA recursion of PFQN_MVA, which recurs on the population vector and therefore costs O(prod(N+1)), MVAC recurs on the chains: each chain is reduced to single-customer chains and the removed chains are replaced by self-looping single-customer (SCSL) chains pinned at a service center. The multiplicity vector v = (v_1,…,v_J), where v_j is the number of SCSL chains at center j, indexes the recursion in place of the population vector. MVAC is therefore attractive for networks with few centers and many chains, where its cost grows only polynomially in the number of distinct chains, and it 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.

Throughout, j and i index service centers (j = 1,…,J1 are SSFR and j = J1+1,…,J are IS), k and l index single-customer chains, a_jk = theta_jk T_jk is the relative utilization of chain k at center j (i.e., its service demand), a_k = sum_j a_jk, and I_k = {v : sum_j v_j = K - k} with K the total number of single-customer chains. Writing L^k_j(v) for the mean number of customers at center j (SCSL customers excluded), L^k_{jl}(v) for the mean number of chain-l customers at center j, and lambda^k_k(v) for the throughput of chain k, all for the network with normalizing constant G_k(v), 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. 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 that visit only SSFR centers require instead 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 from the first execution. 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 rather than on K.

Parameters:

L - (M x R)
N - (1 x R)
Z - (1 x R) or (Mz x R) – one IS center. Defaults to zeros(1,R).

Returns:

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

See also PFQN_MVA, PFQN_RECAL, PFQN_CONV.

pfqn_mcub(L, N, Z)

Multiclass Composite Upper Bound (Kerola 1986) on per-class throughput for closed product-form networks.

Kerola’s composite bound method (Perf. Eval. 6:1-9, eqs. 10-16). Given per-class multiclass BJB lower bounds X_s^-, the residual- utilization composite upper bound is X_r <= min_k [1 - sum_{s!=r} X_s^- L_ks] / L_kr, computed at O(KR). Much tighter than per-class ABA at moderate load. (Named pfqn_mcub because pfqn_cub is the unrelated cubature NC method.)

Parameters:

L – Service demand matrix, station x class (M x R).
N – Population vector (1 x R).
Z – Think time vector (1 x R, default zeros).

Returns:

Xub – Composite upper throughput bound per class (1 x R). Xlb: Multiclass BJB lower throughput bound per class (1 x R, eq. 10).

pfqn_sdrcoeff(sdr)

Derived coefficients of a Krzesinski state-dependent routing structure.

Parameters:: sdr – State-dependent routing structure.
Returns:: c – Derived coefficient structure.

C = PFQN_SDRCOEFF(SDR)

Validates a state-dependent routing (SDR) structure and returns the derived coefficients of Krzesinski (1987), “Multiclass Queueing Networks with State-Dependent Routing”, Performance Evaluation 7:125-143, eqs. (11)-(14).

SDR describes the partition of the network into a subnetwork Q(V,V) subject to SDR and its complement Q(N-V,M-V), with Q(V,V) split into branches arranged in a hierarchy of nested subnetworks V_1 > V_2 > … > V_T. Its fields are, following the paper’s own indexing in which branch 1 is the complement M-V and the SDR branches are numbered 2..B:

sdr.entry station index of the entry center e of Q(V,V) sdr.departure station index of the departure center d of Q(V,V) sdr.branch 1xB cell, branch{b} = station indices of branch b (b>=2);

branch{1} is unused and holds the complement implicitly

sdr.entryOf 1xB, entryOf(b) = station index of the entry center e(b) sdr.departureOf 1xB, departureOf(b) = station index of d(b) sdr.level 1xB, level(b) = the unique t with B_b in V_t - V_{t+1} sdr.C 1xT, the coefficients C_t of eq. (11) sdr.d TxB, d(t,b) defined for 2<=b<=B and 1<=t<=level(b)

The returned C has fields

c.T, c.B, c.level, c.C, c.d as above; c.inA{t} branch indices b with level(b) >= t, i.e. the set A_t of

branches contained in the subnetwork Q(V_t,V_t)

c.Dtt(t) D_tt = sum_{b in A_t} d(t,b), eq. (14) c.Dprev(t) D_{t-1,t} = sum_{b in A_t} d(t-1,b), eq. (14), t>1 c.mmax(b) largest m_b with delta_{level(b),b}(m_b) >= 0 c.vmax(t) largest v_t with omega_tt(v_t) >= 0 and, for t>1, with

omega_{t-1,t}(v_t) >= 0

The bounds c.mmax and c.vmax are the population constraints that SDR imposes on branches and on subnetworks (Sec. 2.5); they are consequences of the routing coefficients, not independent inputs. A bound is Inf when the corresponding C_t is nonnegative, in which case SDR imposes no constraint.

pfqn_chow(L, N, Z, tol, maxiter, QN0, type, variant)

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 over (2.2)-(2.6) is run. Two estimators are given:

‘backward’ theta_ck = [Qhat_k(N - 1_c) - Qhat_k(N)] / Qhat_k(N) (2.16) ‘forward’ theta_ck = [Qhat_k(N) - Qhat_k(N + 1_c)] / Qhat_k(N + 1_c) (2.17)

Chow reports the forward form to be the more accurate of the two, so it is the default here. Setting every theta to zero recovers pfqn_lcp.

Parameters:

L – Service demand matrix (stations x classes).
N – Population vector.
Z – Think time vector.
tol – Tolerance for convergence.
maxiter – Maximum number of iterations.
QN0 – Initial guess for queue lengths.
type – Scheduling strategy type (default: PS).
variant – ‘forward’ (eq. 2.17, default) or ‘backward’ (eq. 2.16).

Returns:

XN – System throughput. QN: Mean queue lengths. UN: Utilization. RN: Residence times. it: Number of iterations performed.

pfqn_jdfun(nvec, jdscaling, classIdx)

AMVA joint-dependence function for non-product-form scaling.

Parameters:

nvec – Population state vector.
jdscaling – Cell array of joint-dependent scaling functions.
classIdx – Optional class index selecting eta_{i,r} (default: 1).

Returns:

r – Scaling factor vector for each station.

R = PFQN_JDFUN(NVEC, JDSCALING, CLASSIDX)

AMVA joint-dependence function. Returns, for every station i, the reciprocal of the joint-dependent scaling

eta_i(n_i1, …, n_iR)

evaluated at the per-class population vector NVEC(i,:), for class r=CLASSIDX.

JDSCALING{i} is a function handle of the joint per-class population vector at station i. It may return either

a scalar, i.e. a scaling eta_i(n) shared by every class (broadcast, as in the flagship min(ni(1),c)), or

a vector of length R, i.e. per-class scalings [eta_{i,1}(n), …, eta_{i,R}(n)], of which element CLASSIDX is taken (Sauer chain-dependent rate mu_{r,i}(n)).

Unlike PFQN_CDFUN (product-form beta_{i,r} depending on the own-class marginal n_{i,r}), eta may read the joint vector arbitrarily and is therefore NON-product-form: the AMVA result is an approximation with no exactness/uniqueness guarantee. The numerical evaluation matches PFQN_CDFUN; the distinction is semantic (product-form vs joint) and is carried by the separate sn.jdscaling field.

pfqn_dac(L, N, Z, mu)

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

Parameters:

L – Service demand matrix (MxR).
N – Population vector (1xR).
Z – Think time vector (1xR, default: zeros).
mu – Load-dependent rate matrix (MxNt, default: ones).

Returns:

Pjoint – Joint queue-length probabilities over the aggregate state space. states: Aggregate states, one per row of Pjoint. XN: Chain throughputs. QN: Mean queue lengths. UN: Station utilizations. CN: Cycle times. pi: Marginal queue-length probabilities.

[PJOINT,STATES,XN,QN,UN,CN,PI]=PFQN_DAC(L,N,Z,MU)

Distribution Analysis by Chain (DAC) for closed product-form queueing networks with single-server fixed-rate, infinite-server and queue-dependent service centers. Unlike MVA, RECAL or MVAC, the recursion returns the whole set of joint queue-length probabilities, which are required e.g. in availability modeling.

The recursion proceeds chain by chain 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.

Inputs:

L (MxR) service demand of chain r at station j N (1xR) number of customers of chain r Z (1xR) think time of chain r. 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 (MxNt) service rate of station j with n customers, Nt=sum(N).: Defaults to ones(M,Nt), i.e. single-server fixed rate. Use mu(j,:)=1:Nt for infinite server and mu(j,n)=min(n,c) for a c-server station.

Outputs:

Pjoint (SxR) probability of the aggregate state in the corresponding row: of STATES, S=nchoosek(Nt+J-1,J-1) with J the number of centers

states (SxJ) aggregate states, states(s,j) = customers at center j XN (1xR) throughput of chain r QN (MxR) mean number of chain-r customers at station j UN (Mx1) utilization of station j, i.e. 1-P_j(0) CN (1xR) cycle time of chain r, exclusive of think time pi (Mx(Nt+1)) pi(j,n+1) = marginal probability of n jobs at station j

Example (availability model, de Souza e Silva 1987, Section 3):

L = [5,0; 0,10; 2,1]; N = [1,3]; mu = [1,1,1,1; 1,2,2,2; 1,1,1,1]; [P,S] = pfqn_dac(L,N,[0,0],mu); AV = sum(P(S(:,1)==1 & S(:,2)>=1)); % system available

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

pfqn_cbh(L, N, Z, level)

Convolutional Bound Hierarchy (Dowdy, Eager, Gordon, Saxton 1984) for single-class closed product-form networks.

Level-level convolutional bound hierarchy on throughput. Fills column c = M-level of Buzen’s g-array with BJB-derived estimates (e_0=1, e_i=e_{i-1}/BOUND(i)), then convolves the remaining level servers exactly (g(n,m)=g(n,m-1)+L_m g(n-1,m)). Bounds tighten monotonically with level and equal the exact solution at level M (Dowdy et al. 1984). BJB upper fill yields the upper bound; BJB lower fill the lower bound.

Parameters:

L – Service demand vector (M x 1).
N – Population (scalar).
Z – Think time (scalar, default 0); folded as an extra IS column.
level – Number of exactly-convolved servers, 1..M (default 2).

Returns:

Xlo – Lower throughput bound. Xhi: Upper throughput bound.

pfqn_rgfmc(L, N, Z, options)

Multiclass Recursion by Generating Functions (RGF): the normalizing constant by iterated residues, with think times.

Exact normalizing constant of a closed MULTICLASS product-form network by eliminating one class at a time by residues, finishing in the single-class convolution of pfqn_rgf.

1. Harrison, S. Coury, “On the asymptotic behaviour of closed multiclass queueing networks”, Performance Evaluation 47:131-138, 2002, Thm 1, expresses the generating function of a q-class network in terms of those of (q-1)-class networks; P. G. Harrison, T. T. Lee, “A new recursive algorithm for computing generating functions in closed multi-class queueing networks”, IEEE MASCOTS 2004, eqs. (4)-(5), turns it into the RGF algorithm, bottoming out in a collection of single-class normalizing constants memoised by load vector (Sec. 3.4).

Parameters:

L – Service demand matrix (M x R).
N – Population vector (1 x R), nonnegative integers.
Z – Think time vector (1 x R). Default: zeros.
options – Struct with optional fields: .tol relative tolerance for calling two affine forms

proportional, hence one pole. Default: 1e-12.

.maxterms cap on residue terms carried between eliminations.
Default: 1e6.

.maxcancel nats of cancellation tolerated before refusing.
Default: 15.

Returns:

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

Think times are not in either paper:

Both write the generating function as the RATIONAL: H_q(M,X;z) = prod_i (1 - rho_i z)^-m_i,

with every node a load-independent single server. An infinite server multiplies this by the ENTIRE exp(sum_r Z_r z_r), and that breaks the step Thm 1 rests on: G_n(z’) = -sum_i r_i holds only because the residues of n(z)/d(z) sum to zero when deg d >= deg n + 2 (Bertozzi and McKenna, SIAM Review 35(2):239-268, 1993, fact (IV), p. 246), and an exponential numerator does not decay at infinity.

The delay is therefore carried by their own repair, eqs. (3.19)-(3.21) of the same paper: only the first k_r+1 Taylor coefficients of exp(Z_r z_r) can reach the coefficient of z_r^k_r, so replacing the exponential by that polynomial is EXACT, not an approximation, and leaves a rational integrand the residue calculus handles unchanged. The price is that the eliminated class’s population re-enters the term count, which is exactly the population-insensitivity Harrison-Lee Sec. 4 advertises; the class kept for the base case pays nothing.

Degeneracy:

Thm 1 assumes rho_iq ~= rho_lq and its Conclusion leaves the tied case open. Two affine forms name the SAME pole only when they are PROPORTIONAL, so fusing proportional forms into one factor of summed multiplicity disposes of the degeneracy with nothing else changed; a tie in the eliminated class alone merely leaves a form with no constant term, which the recursion carries.

Conditioning:

The elimination is exact in exact arithmetic but is an ALTERNATING sum over residues, so near-coincident loads over an eliminated class destroy significance. The worst cancellation ratio is tracked and the routine REFUSES past options.maxcancel rather than returning a confidently wrong lG. Everything else runs in the log domain, so no Poisson weight or binomial is ever formed as a naive ratio.

pfqn_cntol(N)

Chandy-Neuse population-scaled termination cutoff for approximate MVA.

Parameters:: N – Population vector.
Returns:: tol – Termination cutoff 1/(4000+16*sum(N)).

Published cutoff of the Linearizer termination test: K. M. Chandy, D. Neuse, “Linearizer: A Heuristic Algorithm for Queuing Network Models of Computing Systems”, Commun. ACM 25(2):126-134, 1982, p.129 and appendix. The iteration continues while

max_{i,r} |Q^I(i,r) - Q^{I-1}(i,r)| / N_r > 1/(4000 + 16*|N|),

|N| = sum(N). The paper motivates the scaling with |N|: at large populations removing one job changes the queue lengths very little, so a fixed cutoff would terminate the iteration prematurely. It also notes that the expression stays below 0.00025 even at very small populations.

The same expression is what LQNS uses as its termination test, set in the SchweitzerCommon constructor of libmva/src/mva.cc; that code carries no citation, and the paper above is its source.

Pass the result as the tol argument of pfqn_bs / pfqn_egflinearizer only if the plain cutoff is wanted with those functions’ own convergence metric. Passing tol = ‘cn’ (or NaN) instead selects BOTH this cutoff and the normalized-maximum metric of the paper, which is the published test.

pfqn_lldsingle(L, N, mu, options)

Exact normalizing constant for single-class limited load-dependent models.

Same recursion, same arithmetic and same result as pfqn_gldsingle, but with the auxiliary rate-offset axis truncated at the LIMITED LOAD-DEPENDENCE threshold instead of at the population. Unrolling the recursion of pfqn_gldsingle shows that its third index is an offset into station m’s rate function,

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) and the product collapses to (L_m/alpha_m(s_m))^j. Hence

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

and the N-s_m upper slices that pfqn_gldsingle computes are duplicates of one another. Capping the offset at s_m and reading g(m,n-1,min(t+1,s_m)) keeps every value it needs.

COST. O(N * sum_k s_k) time against O(M N^2) for pfqn_gldsingle, and O(N * max_k s_k) space against O(M N^2), the levels being rolled. On a multiserver model, where s_k is the server count, this is LINEAR in the population rather than quadratic. The two agree to the last bit, since the arithmetic performed is a subset of pfqn_gldsingle’s: the log-domain branch remains a sum of nonnegative terms and loses no digits to cancellation, unlike the closed form of pfqn_explicit_ld, which reaches the same asymptotics through an alternating sum.

There is no gain on a station whose rates never settle, an infinite server alpha(n)=n being the usual case: it gets s_k = N and costs what it costs in pfqn_gldsingle. The saving is over the OTHER stations, so a model carrying one delay among M queues drops from O(M N^2) to O(N^2 + N sum_k s_k).

The threshold is detected per station rather than declared, so an arbitrary rate matrix is accepted and simply yields s_k = N, at which point this is pfqn_gldsingle with its slices rolled. A MISSED tie only costs time; a FALSE tie would be a wrong answer, hence the strict tolerance, which follows pfqn_explicit_ld’s scan.

Parameters:

L – Service demand vector (Mx1).
N – Population (scalar).
mu – Load-dependent rate matrix (MxN), alpha_i(j) = mu(i,j).
options – Solver options.

Returns:

lG – Logarithm of normalizing constant. G: Normalizing constant. s: Detected per-station thresholds (Mx1), alpha_i(n)=alpha_i(s_i) for n>=s_i.

G=PFQN_LLDSINGLE(L,N,MU)

pfqn_qsa(L, N, Z, type, tol, maxiter, levels, QN0)

Queue-Shift Approximation (QSA) for closed product-form networks.

Schweitzer, Serazzi and Broglia (Tools’98, LNCS 1469, pp. 267-279) approximate the arrival-instant queue lengths through the absolute shift Y_ri(K) = 1 + Q_i(K-e_r) - Q_i(K) of the aggregate queue length, in place of the fractional deviations of Linearizer. The aggregate core problem (eq. 13a) is imposed at K, at every K-e_s and, in the three-level variant of eq. (16), at every K-e_s-e_t with the affine extrapolation of eq. (15).

The quintuple (16) is solved as a single system by the damped Newton method of Sect. 4. Successive substitution over the decomposed core problems is not used: it is unstable near saturation and converges to the degenerate root in which the bottleneck absorbs the whole population.

Parameters:

L – Service demand matrix (stations x classes).
N – Population vector.
Z – Think time vector.
type – Scheduling strategy per station; SchedStrategy.INF marks a delay centre.
tol – Residual tolerance (default: 1e-10).
maxiter – Maximum number of Newton iterations (default: 100).
levels – 2 for the two-level QSA of eq. (14), 3 for eq. (16) (default: 3).
QN0 – (M x R) queue lengths that warm-start the Bard-Schweitzer initialization; empty for the default cold start.

Returns:

Q – Mean queue lengths. U: Utilization. W: Residence times. C: Cycle times. X: System throughput. totiter: Newton iterations performed.

pfqn_ncoi(Z, N, mu, visits, options)

[G, LG, GTAB] = PFQN_NCOI(Z, N, MU, OPTIONS)

Normalizing constant for a closed product-form queueing network that comprises a single aggregated infinite-server (delay) node and an arbitrary number of order-independent (OI) / pass-and-swap stations with empty swap graph.

The OI stations are analyzed by the balanced-fairness recursion of Bonald & Proutiere (2003), “Insensitive bandwidth sharing in data networks”, 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 obtained by convolving the per-station balance functions with the multinomial delay factor F_Z(n) = prod_r Z_r^{n_r} / n_r!,

g_0(n) = F_Z(n), 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: both the balance functions and the convolution are tabulated over the count lattice 0 <= n <= N, never over orderings. That is legitimate exactly because an OI rate is permutation- invariant, so Phi(n) – itself the sum of the ordered-prefix weights prod_p 1/mu(c_1..c_p) over all orderings c of the multiset n – closes on the count vector. With a non-empty swap graph that closure fails and the microstate routine PFQN_PAS_NC must be used instead.

COST. With L = prod_r (N_r+1) the balance functions cost O(M R L) and the convolutions O(M sum_{n<=N} prod_r (n_r+1)) = O(M prod_r (N_r+1)(N_r+2)/2), i.e. the order of a load-dependent Buzen convolution: polynomial in the population for a fixed number of classes.

Parameters:

Z - (1 x R) – Z(r) = 1/sigma_r for a delay with per-class rate sigma_r.
N - (1 x R)
mu - cell array {1 x M} of function handles, one per OI station. Each – mu{m}(n) returns the total service rate of station m given the per-class occupancy (count) vector n (1 x R). For an order- independent station this rate depends only on the support of n (which classes are present), i.e. mu{m}(n) = sum of the capacities of the servers compatible with the classes present in n. A station state whose rate is non-positive is unreachable and is assigned a zero balance value. May be empty to model a pure delay network.
options - solver options (optional, currently unused; accepted for – signature parity with the other pfqn_* routines).

Returns:

G - Normalizing constant G(N). lG - log(G(N)). Gtab - (prod_r(N_r+1) x 1) the WHOLE lattice of normalizing constants,

Gtab(1 + sum(n .* strides)) = G(n) for every 0 <= n <= N, with strides = [1, cumprod(N(1:end-1)+1)]. The convolution produces this table anyway, so a caller that needs G at more than one population (throughputs, queue lengths, a fold of further stations) must take this output and index it, NOT re-call the routine per population – the latter costs a needless factor prod_r(N_r+1).

Example (IS + two OI stations, R classes):: oirate = @(n) sum(mu1(any(compat1(:, find(n>0)) ~= 0, 2))); oirate2 = @(n) sum(mu2(any(compat2(:, find(n>0)) ~= 0, 2))); G = pfqn_ncoi(1./sigma, N, {oirate, oirate2});

See also PFQN_PAS_NC, PFQN_OI_FNC, PFQN_OI_INSVC, PFQN_MVAOI.

pfqn_marie(L, N, Z, scv, varargin)

Marie’s iterative aggregation for closed networks with FCFS Coxian (non-exponential) service. Single-class exact-reducing; multiclass via QD-AMVA with class-dependent (cd) scaling.

Marie’s method (Marie 1979/1980): approximate mean performance of a closed queueing network with FCFS general (Coxian) service, via iterative aggregation-decomposition. Each station is analyzed in isolation as a lambda(n)/Cox/1 queue; the resulting conditional throughputs mu_i(n) drive a load-dependent aggregate solve, iterated to a fixed point. Single class (R=1): the aggregate is the exact LD product-form solve pfqn_mvald, and the method reduces to exact product form for exponential service. Multiple classes (R>1): the aggregate is QD-AMVA with class-dependent (cd) scaling beta_{i,r}(nvec) supplied by a multiclass Cox/1 isolation sub-model; exact product-form service (scv==1 with class-independent means) is dispatched to exact MVA, otherwise the result is a decomposition approximation.

Parameters:

L – Service demand matrix (M x R).
N – Population vector (1 x R).
Z – Think time vector (1 x R; total delay demand per class).
scv – Per-station per-class squared coefficient of variation (M x R). scv==1 exponential; scv<0.5 Erlang; scv>0.5 two-phase Coxian.
varargin – Optional trailing arguments, in order: tol, the convergence tolerance (default 1e-8); maxiter, the iteration cap (default 1000); nservers, the per-station server count (M x 1, default all 1), single-class only (multiserver multiclass isolation is not yet supported).

Returns:

X –

Throughput: single class M x 1 (per station); multiclass 1 x R: (per-class chain throughput, visits folded into L).

Q: Mean queue length (M x 1 single class, M x R multiclass). U: Utilization (same shape as Q). C: Residence time (same shape as Q). it: Iterations performed. mu: Converged LD data (M x N single class; cell of cd-scalings R>1).

pfqn_sens_validate()

PFQN_SENS_VALIDATE Numerically validate the queue-length moment / demand- derivative identities for closed product-form (BCMP) queueing networks.

Notation (single-server load-independent stations, plus optional delay Z):

Q_{i,r}(N) mean queue length of class r at station i, pop. N: -> pfqn_mva(L,N,Z), QN(i,r)
Q_{i,s}^{+k}(N) queue length of class s at ORIGINAL station i in the: network obtained by adding one replica of station k (a station with identical demands L(k,:)) -> pfqn_mva([L;L(k,:)],N,Z), read row i

D_{i,r} = L(i,r) service demand dQ_{k,r}/dD_{i,s} -> pfqn_sens(L,N,Z), dQ(k,r,pL(i,s)) N - 1_r population with one class-r job removed

The three “moment” identities and the three demand-derivative formulas from the prompt are each checked over all admissible index tuples (i,k,r,s) on a set of random closed networks. Analytic derivatives (pfqn_sens) are also cross-checked against central finite differences of pfqn_mva.

pfqn_nre(L, N, Z, alpha, options, vfix)

Normalizing constant via saddle-tilted Edgeworth (NRE) approximation.

Evaluates the Norlund-Rice integral form of the limited-load-dependent normalizing constant (Casale-Harrison-Ong, Perform. Eval. 152, 2021, Thm. 6) by steepest descent instead of a Laplace approximation on the untilted contour, as done by pfqn_nrl and pfqn_nrp. Two corrections are applied over those methods:

The integrand is invariant under t -> t + c*1, since h is homogeneous of degree sum(N) in the class variables and that degree cancels against exp(-1i*N*t). The redundant direction is quotiented out, so the integral is (R-1)-dimensional, not R-dimensional. pfqn_nrl and pfqn_nrp integrate over R dimensions and let the logit/probit Jacobian supply curvature along the null direction, which is an artifact of the change of variables rather than of the integrand. 2. The contour radii are tilted per class to the saddle point, i.e. to the X solving X_r*dlog(h)/dX_r = N_r, so that the origin is a stationary point of the phase. On the untilted contour X = 1 it is not, which is the leading source of bias in pfqn_nrl and pfqn_nrp.

A second-order Edgeworth term built from the third and fourth cumulants of the tilted distribution is then added, giving a relative error of O(1/sum(N)^2) instead of O(1) for a heuristic Gaussian fit.

All integrand evaluations are at real positive demands, so unlike pfqn_nrl and pfqn_nrp this method needs no complex arithmetic and runs entirely in the log domain through pfqn_lldsingle, whose cost is linear rather than quadratic in the population wherever the rates settle.

Cost is O(I*R^2 + R^4) evaluations of a single-class LLD normalizing constant, hence polynomial in the number of classes.

Parameters:

L – Service demand matrix (MxR).
N – Population vector (1xR).
Z – Think time vector (1xR or DxR).
alpha – Load-dependent rate matrix (Mx sum(N)).
options – Solver options.
vfix – Optional tilt to use instead of solving the saddle-point equation. 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; omit it for the standard estimator.

Returns:

lG – Logarithm of normalizing constant. lGs: Logarithm of the saddlepoint term alone, i.e. of lG with the

Edgeworth factor omitted, so that lG-lGs is the correction.

vsad: The tilt actually used, for reuse at a nearby population. G: Normalizing constant.

pfqn_minclasses(Usum, K, N)

Lower bound on the number of customer classes needed to explain a measured sum of device utilizations, from Dowdy et al. (1992), J. ACM 39(1), Section 4.7.

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, which is the point: it tells a clustering analysis how many classes it must at least find. An upper bound on R is meaningless (extra classes can always be introduced by splitting), so none is returned.

The paper’s example: K = 2 devices, N = 3 customers, measured sum_k U_k = 1.6. A single class admits at most 2N/(N+1) = 1.5, so the single-class assumption is unjustified and Rmin = 2.

Parameters:

Usum – Measured sum of device utilizations sum_k U_k (scalar).
K – Number of devices (scalar).
N – Total number of customers (scalar).

Returns:

Rmin –

Least R in 1..N with pfqn_usumbound(R,K,N) >= Usum; NaN when: Usum exceeds min(N,K) and so is unattainable by ANY class structure, which signals a measurement or bookkeeping error rather than a workload that needs more classes.

pfqn_ncldmx(lambda, D, N, Z, mu, S, varargin)

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

The closed-conditional normalizing constant of a mixed limited load-dependent (LLD) network equals a purely closed load-dependent normalizing constant in which every queueing station i carries the Bruell-Balbo-Afshari effective capacity rate mu_i^eff(n) = 1/EC_i(n), where EC is returned by pfqn_ldmx_ec and folds the open classes into the closed subnetwork. The open classes contribute the separable prefactor lGopen = sum_i log E_i(0), which reduces to -sum_i log(1-rho_i) in the load-independent limit.

Mean measures follow from that identification without ever enumerating the closed population lattice, which is what makes this the normalizing-constant counterpart of pfqn_mvaldmx rather than a rename of it:

closed throughputs are the ratios X_r = G(N-e_r)/G(N); - closed queue lengths are the conditional normalizing-constant recursion of the load-dependent closed network (pfqn_mushift / pfqn_fnc), applied to the effective-capacity rates; - open queue lengths are the Bruell-Balbo-Afshari sum Q_ir = lambda_r D_ir sum_n (n+1) EC_i(n+1) P_i(n) with its SATURATED TAIL FOLDED ONTO THE CLOSED MEAN. EC_i(n) is constant for n >= b_i, the level where the rate row stops growing, so writing EC_i(n) = EC_i^inf + delta_i(n) with delta_i(n)=0 for n >= b_i leaves Q_ir = lambda_r D_ir [ EC_i^inf (Q_i^closed + 1) + sum_{n=0}^{b_i-2} (n+1) delta_i(n+1) P_i(n) ] using sum_n P_i(n)=1 and sum_n n P_i(n)=Q_i^closed. Only the first b_i-1 marginal probabilities survive, and b_i is the number of servers, not the population: a single-server station needs none at all, and the formula collapses to the classical lambda_r D_ir (1+Q_i^closed)/(1-rho_i).

The marginals that remain are themselves normalizing-constant ratios, P_i(n) = sum_{|k|=n} F_i(k) G_{-i}(N-k) / G(N), evaluated as in solver_nc_margaggr.

Reproduces pfqn_mvaldmx to machine precision on multiserver, nonlinear-mu, think-time and multi-chain models; see test_pfqn_ncldmx.

Parameters:

lambda – Arrival rate vector (0 on closed classes).
D – Service demand matrix (MxR).
N – Population vector (Inf on open classes).
Z – Think time vector (closed classes).
mu – Load-dependent rate matrix (Mx>=sum(N_closed)).
S – Number of servers per station (currently informational).
varargin – Optional solver parameters forwarded to pfqn_ncld.

Returns:

lG – Logarithm of the closed-conditional normalizing constant. G: Closed-conditional normalizing constant (exp(lG)). lGopen: Logarithm of the open-class normalizing prefactor sum_i log E_i(0). XN: Throughputs (1xR): G(N-e_r)/G(N) on closed classes, lambda_r on open ones. QN: Mean queue lengths (MxR).

[LG,G,LGOPEN,XN,QN] = PFQN_NCLDMX(LAMBDA,D,N,Z,MU,S,VARARGIN)

pfqn_mvams_ilock(lambda, L, N, Z, mi, S, IL)

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

Parameters:

lambda – Arrival rate vector.
L – Service demand matrix.
N – Population vector.
Z – Think time vector.
mi – Queue replication factors (default: ones).
S – Number of servers per station (default: ones).
IL – Interlock matrix (R x R), see PFQN_MVA_ILOCK.

Returns:

XN – System throughput. QN: Mean queue lengths. UN: Utilization. CN: Residence times (M x R). lG: Always NaN.

[XN,QN,UN,CN,LOGG]=PFQN_MVAMS_ILOCK(LAMBDA,L,N,Z,MI,S,IL)

The interlock is defined only for closed single-server models, so this is the one shape accepted here; anything else is refused rather than served without the correction. Models with no interlock go to PFQN_MVAMS.

pfqn_ldmx_ec(lambda, D, mu)

Compute effective capacity terms for MVALDMX solver.

Parameters:

lambda – Arrival rate vector.
D – Service demand matrix.
mu – Load-dependent rate matrix.

Returns:

EC – Effective capacity matrix. E: E-function values. Eprime: E-prime function values. Lo: Open class load vector.

[EC,E,EPRIME,LO] = PFQN_MVALDMX_EC(LAMBDA,D,MU) Compute the effective capacity terms in MVALDMX Think times are not handled since this assumes limited load-dependence

pfqn_xzgsbup(L, N, Z)

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

Parameters:

L – Service demand vector.
N – Population.
Z – Think time.

Returns:

X – Upper bound on throughput.

pfqn_xzgsblow(L, N, Z)

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

Parameters:

L – Service demand vector.
N – Population.
Z – Think time.

Returns:

X – Lower bound on throughput.

pfqn_xzabaup(L, N, Z)

Upper asymptotic bound on throughput (Zahorjan-Balanced).

Parameters:

L – Service demand vector.
N – Population.
Z – Think time.

Returns:

XN – Upper bound on throughput.

pfqn_sdr(S, xi, N, sdr, alpha)

Exact product form of a multiclass network with state-dependent routing.

Parameters:

S – Mean service times per station and chain.
xi – Relative visit counts per station and chain.
N – Chain population vector.
sdr – State-dependent routing structure.
alpha – Optional load-dependent rate scalings.

Returns:

Q – Mean queue lengths per station and chain. X: Mean throughputs per station and chain. U: Mean utilizations per station and chain. R: Mean response times per station and chain. G: Normalizing constant. lG: Logarithm of the normalizing constant. prob: Stationary probability of every enumerated state. states: Enumerated state list.

[Q, X, U, R, G, LG, PROB, STATES] = PFQN_SDR(S, XI, N, SDR, ALPHA)

Exact evaluation of the product-form joint probability distribution of a multiclass closed queueing network with state-dependent routing, Krzesinski (1987), “Multiclass Queueing Networks with State-Dependent Routing”, Performance Evaluation 7:125-143, eq. (16):

P(n) = G^-1 prod_i f_i(n_i)
prod_t [Omega_{t-1,t}(v_t)/Omega_tt(v_t)] prod_{b in A_t - A_{t+1}} Delta_tb(m_b)

with f_i(n_i) = [n_i!/beta_i(n_i)] prod_j gamma_ij^{n_ij}/n_ij!, gamma_ij = xi_ij/mu_ij, beta_i(n) = alpha_i(n) beta_i(n-1), beta_i(0) = 1, and the cumulative coefficients

Omega_{t-1,t}(v) = omega_{t-1,t}(v-1) Omega_{t-1,t}(v-1), Omega(0) = 1, Omega_tt(v) = omega_tt(v-1) Omega_tt(v-1), Omega(0) = 1, Delta_tb(m) = delta_tb(m-1) Delta_tb(m-1), Delta(0) = 1.

This form is general in the branch topology: a branch may hold several interconnected centers. Only the MVA and convolution algorithm of the paper’s Section 4, implemented in PFQN_SDRMVA, is restricted to single-center branches. The normalizing constant here is obtained by summing the unnormalized weights over the whole reachable state space, which is exact for any branch topology but grows combinatorially with the populations.

Inputs:

S MxJ mean service times 1/mu_ij of a chain j customer at center i XI MxJ coefficients xi_ij of Section 3.2. For the entry and departure

centers of Q(V,V) and of every branch these are all equal and are NOT relative visit counts; use PFQN_SDRVISITS to obtain them

N 1xJ chain populations SDR state-dependent routing structure, see PFQN_SDRCOEFF ALPHA optional Mxmax(N) matrix of load-dependent rate scalings, with

ALPHA(i,k) = alpha_i(k) the rate multiplier when k customers are present at center i. Defaults to 1, a fixed-rate center. Use ALPHA(i,k) = k for an infinite-server center and min(k,c) for a c-server center

S and XI are required separately rather than as their product: under SDR the xi_ij are not visit ratios, so the per-center throughputs cannot be recovered from the demands alone.

Outputs Q, X, U and R are MxJ. X holds the per-center chain throughputs: T_ij = sum_n P(n) alpha_i(n_i) (n_ij/n_i) / S_ij,

U_ij = T_ij S_ij is the mean number of chain j customers in service, and R_ij = Q_ij / T_ij is the mean response time at the center.

pfqn_unique(L, mu, gamma)

Consolidate replicated stations into unique stations with multiplicity.

Parameters:

L – Service demand matrix (M x R).
mu – Load-dependent rate matrix (M x Ntot), optional.
gamma – Class-dependent service rate matrix (M x R), optional.

Returns:

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

PFQN_UNIQUE Consolidate replicated stations into unique stations with multiplicity

[L_UNIQUE, MU_UNIQUE, GAMMA_UNIQUE, MI, MAPPING] = PFQN_UNIQUE(L, MU, GAMMA)

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.

Input:: L - M x R demand matrix mu - M x Ntot load-dependent rate matrix (optional, pass [] if not used) gamma - M x R class-dependent service rate matrix (optional, pass [] if not used)
Output:: L_unique - M’ x R demand matrix with M’ <= M unique stations mu_unique - M’ x Ntot reduced load-dependent rates (empty if mu was empty) gamma_unique - M’ x R reduced class-dependent rates (empty if gamma was empty) mi - 1 x M’ multiplicity vector mapping - 1 x M vector mapping original station i to unique station index

pfqn_stdf_heur(L, N, Z, S, fcfsNodes, rates, tset)

Heuristic sojourn time distribution for multiserver FCFS nodes (McKenna 1987 variant).

Parameters:

L – Service demand matrix.
N – Population vector.
Z – Think time vector.
S – Number of servers per station.
fcfsNodes – Indices of FCFS nodes to analyze.
rates – Service rates (MxR matrix).
tset – Time points for CDF evaluation.

Returns:

RD – Cell array of sojourn time distributions {station,class}.

Heuristic sojourn time distribution analysis at multiserver FCFS nodes based on a variant of the method in J. McKenna 1987 JACM

pfqn_stdf(L, N, Z, S, fcfsNodes, rates, tset)

pfqn_scat(L, N, Z, type, tol, maxiter, QN0)

Neuse-Chandy SCAT approximate mean value analysis.

SCAT shares the Linearizer fixed point. It carries the mean queue lengths at the target population N and at the R reduced populations N-e_s, and corrects the Bard-Schweitzer proportionality assumption with the fraction difference

Delta(i,r,s) = Q(i,r|N-e_s)/(N-e_s)_r - Q(i,r|N)/N_r,

held fixed while an inner MVA fixed point is iterated. It differs from Linearizer in that this correction is refreshed ONCE: SCAT stops after the first pass, where Linearizer performs the fixed three passes of Chandy and Neuse (1982), Sec. 4. Cost is therefore about one third of Linearizer’s, and accuracy sits between Bard-Schweitzer (pfqn_bs, the Delta=0 special case) and Linearizer.

SCAT’s second departure from Linearizer, fitting a probability mass function centred on the mean queue length at queue-dependent centres instead of propagating the MVA distribution recursion (Krzesinski and Greyling 1984, Sec. 4), does not arise here: this entry point covers single-server and delay stations only, exactly as pfqn_linearizer does. That mass function is available separately as the ‘scat’ marginal rule of pfqn_ab_amva.

Reference: D. Neuse, K. M. Chandy, “SCAT: A Heuristic Algorithm for Queueing Network Models of Computing Systems”, ACM SIGMETRICS Perform. Eval. Rev. 10(3), 1981; K. M. Chandy, D. Neuse, “Linearizer: A Heuristic Algorithm for Queuing Network Models of Computing Systems”, Commun. ACM 25(2), 1982.

Parameters:

L – Service demand matrix (M x R).
N – Population vector (1 x R).
Z – Think time vector (1 x R) or matrix summed over rows.
type – Scheduling strategy per station; accepted for interface parity with pfqn_linearizer, but the residence-time recursion is discipline-independent.
tol – Convergence tolerance (default: 1e-8); ‘cn’ or NaN selects the Chandy-Neuse (1982) population-scaled termination test, see pfqn_cntol.
maxiter – Maximum inner iterations (default: 1000).
QN0 – (M x R) queue lengths that warm-start the Bard-Schweitzer initialization; empty for the default cold start.

Returns:

Q – Mean queue lengths. U: Utilization. W: Residence times. C: Cycle times. X: Class throughputs. totiter: Total iterations performed.

pfqn_scbgap(N, K, r, undominated)

Maximum relative throughput error of aggregating r customer classes, from Dowdy et al. (1992), J. ACM 39(1), Expressions (3)-(5).

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, and the error is then at most 50%; with r < N it is the partial-aggregation error of their Theorem 4.

The bound depends only on N, K and r, never on the demands, so it can be attached as a certified error bar to any result computed on merged chains – LINE merges classes into chains routinely through sn_get_demands_chain.

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.

Parameters:

N – Total number of customers (scalar). One customer per class, so N is also the number of classes before merging.
K – Number of queueing devices (scalar).
r – Number of classes merged into one (default N, full aggregation).
undominated – True to use the tighter Theorem-5 form, valid only when no class dominates, i.e. every customer’s total device demand is equal, and only for r <= K (default false).

Returns:

e – Maximum relative throughput error, in [0,1/2].

pfqn_sib(L, N, Z, level)

Srinivasan (1985) Successively Improving Bounds (SIB) on cycle time and throughput for single-class product-form closed networks.

Successively Improving Bounds (Srinivasan, “Successively Improving Bounds on Performance Measures for Product Form Queueing Networks”, IEEE ToC 1987 / TR 85-2). Closed-form hierarchy of upper/lower bounds on the cycle time W(N) and throughput X(N) of a single-class closed network of fixed-rate (and delay) stations, based on the MVA relation W(N)=L*(1+phi(N-1)) with phi(K)=sum_m rho_m Q_m(K). Level 1 is the closed form of Thm 2.1; higher levels use the S_i power sums (S_i=sum_m rho_m^i) via Thms 3.5 (upper) and 3.6 (lower), tightening monotonically toward exact. Bounds are always at least as tight as the Balanced Job Bounds. O(M) to compute; level n needs S_2..S_{n+2}.

Parameters:

L – Fixed-rate service demand vector (M x 1). Delay demand goes in Z.
N – Total population (scalar, N>=2).
Z – Think time (scalar; added to the cycle-time bracket).
level – Bound level (>=1, default 3). Higher = tighter, more S_i terms.

Returns:

Xlo – Lower bound on throughput X(N). Xhi: Upper bound on throughput X(N). Wlo: Lower bound on cycle time (residence + think) W(N). Whi: Upper bound on cycle time W(N).

pfqn_recal(L, N, Z, m0)

RECAL (REcursive CALculation) method for normalizing constant.

Parameters:

L – Service demand matrix.
N – Population vector.
Z – Think time vector (default: zeros).
m0 – Initial multiplicity vector (default: ones).

Returns:

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

[G,logG]=PFQN_RECAL(L,N,Z,M0)

pfqn_ble(L, N, Z)

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

Parameters:

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

Returns:

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

[GN,LGN]=PFQN_BLE(L,N,Z)

PFQN_BLE Asymptotic solution of closed product-form queueing networks by logistic expansion, corrected by the deficit that the expansion carries on the one model for which its reference measure is exact.

Cas17 Theorem 2 holds for eps >= eps_N > 0; the K(1+eps*N) self-looping populations are what 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 per Laplaced direction. The size of that bias follows from the control-variate reading of the closure: closing an integral with a reference measure whose exact value is known corrects the Laplace value by the factor the reference itself misses, and for a Dirichlet reference that factor is

-log L_alpha = sum_k r(alpha_k) - r(alpha_0), r(a) = gammaln(a) - ((a-1/2)*log(a) - a + log(2*pi)/2),

the Stirling remainder of its own parameters, with kappa = r(1) = 1-log(2*pi)/2. The reference here is not fitted but read off the balanced model, where the integrand is constant on the simplex and the uniform measure alpha = 1 is exactly right. That gives the two branch constants:

Z=0M*kappa - r(M), the M-1 simplex directions plus the Jacobian, with
the radial integral done exactly as Gamma(N+M);

Z>0M*kappa, the same simplex term plus the radial direction,
which is Laplaced in t=log(v) and contributes exactly +r(M), so the -r(M) cancels.

Both are exact deficits on the balanced model rather than fitted constants. Measured against exact convolution over random Z=0 models the residual after M*kappa-r(M) has median +0.006 nats, against +0.067 for the (M-1)*kappa used previously on this branch. The published expansion is NOT in error; see _kb/03-api-layer.md.

Input: L : MxR demand matrix. L(i,r) is the demand of class-r at queue i N : 1xR population vector. N(r) is the number of jobs in class r Z : 1xR think time vector. Z(r) is the total think time of class r

Output: Gn : estimated normalizing constant lGn: logarithm of Gn. If Gn exceeds the floating-point range, only lGn

will be correctly estimated.

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

pfqn_qzgbup(L, N, Z, i)

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

Parameters:

L – Service demand vector.
N – Population.
Z – Think time.
i – Station index.

Returns:

Qgb – Upper bound on mean queue length at station i.

pfqn_qzgblow(L, N, Z, i)

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

Parameters:

L – Service demand vector.
N – Population.
Z – Think time.
i – Station index.

Returns:

Qgb – Lower bound on mean queue length at station i.

pfqn_lld(L, N, mu, options)

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 turns 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 that pfqn_mushift applies.

pfqn_gld peels station M 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 the unmemoised binary recursion of pfqn_gld, which revisits the same states exponentially often. Without the saturation a memo would still be bounded, but by M * prod_r(N_r+1) * (|N|+1): the threshold is what replaces the population by the server count, exactly as in pfqn_lldsingle.

There is no gain on a station whose rates never settle, an infinite server alpha(n)=n being the usual case: it gets s_k = |N| and the memo keeps its population axis. The saving is over the OTHER stations.

Every shortcut of pfqn_gld is kept and evaluated on the same materialised arguments, so the two agree to the last bit rather than merely to a tolerance: the node reached at (m,n,j) sees the demand block L(1:m,:) and a rate block whose read range 1..sum(n) is entrywise the one pfqn_gld would have built there.

Parameters:

L – Service demand matrix (MxR).
N – Population vector (1xR).
mu – Load-dependent rate matrix (Mx sum(N)), alpha_i(j) = mu(i,j); default all ones.
options – Solver options.

Returns:

G – Normalizing constant. lG: Logarithm of the normalizing constant. s: Detected per-station thresholds (Mx1), alpha_i(n)=alpha_i(s_i) for n>=s_i.

[G,LG]=PFQN_LLD(L,N,MU,OPTIONS)

pfqn_propfair(L, N, Z)

Proportionally fair allocation approximation (Walton 2009).

Parameters:

L – Service demand matrix.
N – Population vector.
Z – Think time vector.

Returns:

G – Normalizing constant estimate. lG: Logarithm of normalizing constant. Xasy: Asymptotic throughput vector.

[G,LOG] = PFQN_PROPFAIR(L,N,Z)

pfqn_perm(A, m)

Permanent of a demand matrix, with optional column multiplicities.

Parameters:

A – Matrix whose permanent is required.
m – Multiplicity of each column of A (optional).

Returns:

val – Permanent value.

VAL = PFQN_PERM(A) permanent of the square matrix A VAL = PFQN_PERM(A, M) permanent of the matrix whose column J is column J of

A repeated M(J) times, so that SUM(M) == SIZE(A,1)

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. See matlab/util/perm.m and _kb/03-api-layer.md.

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

pfqn_procomom(L, N, Z, atol)

ProCoMoM algorithm for computing marginal queue-length probabilities.

Parameters:

L – Service demand matrix (M x R).
N – Population vector (1 x R).
Z – Think time vector (1 x R).
atol – Tolerance.

Returns:

Pr – Marginal probability matrix (M x sumN+1), Pr(k,j+1) = P(n_k = j). Q: Mean queue length vector (M x 1).

pfqn_lekt(L, N, Z)

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

Parameters:

L – Service demand matrix (MxR).
N – Population vector (1xR).
Z – Think time vector (1xR, default: zeros).

Returns:

Gn – Estimated normalizing constant. lGn: Logarithm of normalizing constant. route: ‘kt’ or ‘le’, the side the value was computed on.

[GN,LGN,ROUTE]=PFQN_LEKT(L,N,Z)

PFQN_LEKT The corrected logistic expansion and the corrected Knessl-Tier expansion are ONE estimator, evaluated in M-1 and in R dimensions; this routine computes it on whichever side is cheaper.

With a think time, pfqn_ble (LE plus M units of 1-log(2*pi)/2, one per Laplaced station direction) and pfqn_bkt (KT minus the Stirling remainder s(N_r) of every Laplaced class direction) are the same function of (L,N,Z): 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 on both sides. 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), with r the Stirling remainder of a Gamma direction; the common estimator is defined as the KT value, which is the better of the two on interior modes, 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).

ROUTE. The KT side solves an R-dimensional convex problem and an R x R determinant, the LE side an M-dimensional fixed point and an (M-1) x (M-1) one, so the KT side is taken when R <= M. It is also taken whenever a class self-loops (a single nonzero demand and no think time), since pfqn_kt extracts that class exactly where the logistic expansion only approximates it.

Input: L : MxR demand matrix. L(i,r) is the demand of class-r at queue i N : 1xR population vector. N(r) is the number of jobs in class r Z : 1xR think time vector. Z(r) is the total think time of class r

Output: Gn : estimated normalizing constant lGn : logarithm of Gn. If Gn exceeds the floating-point range, only lGn

will be correctly estimated.

route: ‘kt’ or ‘le’

References: G. Casale. Accelerating performance inference over closed systems by asymptotic methods. ACM SIGMETRICS 2017. C. Knessl, C. Tier. Asymptotic expansions for large closed queueing networks with multiple job classes. IEEE Trans. Computers, 41(4):480-488, 1992.

ljd_delinearize(idx, cutoffs)

NVEC = LJD_DELINEARIZE(IDX, CUTOFFS)

Inverse of LJD_LINEARIZE: recover the per-class population vector from its linearized index.

idx: 1-based linearized index cutoffs: [N1, N2, …, NK] - per-class cutoffs

Returns: [n1, n2, …, nK]

The forward map is idx = 1 + n1 + n2*(N1+1) + n3*(N1+1)*(N2+1) + …, i.e. a mixed-radix numeral with class k in radix (Nk+1), so the inverse is the digit-by-digit division that reads that numeral back.

pfqn_tay(L, N, Z, tol, maxiter, QN0)

Tay’s arrival-instant approximate mean value analysis.

Approximate MVA for closed multiclass product-form networks in which the arrival-instant queue lengths are estimated from the THROUGHPUT ELASTICITIES rather than from a population-shift heuristic (Tay 1987; presented as eqs. 4.8.2-1..3 of the Schweitzer-Serazzi-Broglia survey, where it is benchmarked against exact, Linearizer and Bard-Schweitzer on Tay’s Example 4).

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, so the cost is O(M R (R^3 + M R^2)) per sweep: more than Bard-Schweitzer, less than Linearizer’s R+1 auxiliary networks.

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.

Reference: Y. C. Tay, R. Suri, “Error bounds for performance prediction in queuing networks”, ACM TOCS 3(4), 1985; Y. C. Tay, “An approach to analyzing the behavior of some queueing networks”, Operations Research 40(S2), 1992; P. J. Schweitzer, G. Serazzi, M. Broglia, “A survey of bottleneck analysis in closed queueing networks”, Sec. 4.8.2.

Parameters:

L – Service demand matrix (M x R).
N – Population vector (1 x R).
Z – Think time vector (1 x R). Default: zeros.
tol – Convergence tolerance on the queue lengths. Default: 1e-6.
maxiter – Maximum number of iterations. Default: 1000.
QN0 – Initial guess for the queue lengths (M x R). Default: uniform.

Returns:

XN – Per-class throughputs (1 x R). QN: Mean queue lengths (M x R). UN: Utilizations (M x R). RN: Residence times (M x R). it: Number of iterations performed. QNarr: Arrival-instant queue lengths (M x R x R): QNarr(m,k,r) is the

class-k queue length at station m as seen by an arriving class-r job. These are the auxiliary quantities the method is tabulated on and they are NOT the queue lengths of the model re-solved at N - e_r, which is the same object only for an exact solution.

pfqn_nintmva(L, N, Z)

Mean value analysis at a nonintegral population (fractional-base aMVA).

Exact MVA recursion started from the FRACTIONAL base n_0 = N - floor(N) instead of from the empty network, giving mean performance measures of a single-class closed product-form network at a real-valued population (Dowdy and Gordon 1984, “aMVA”).

The recursion is the standard Reiser-Lavenberg one, R_i(n) = D_i (1 + Q_i(n-1)), X(n) = n/(Z + sum_i R_i(n)), Q_i(n) = X(n) R_i(n), stepped in unit increments from n = n_0 (where the arrival theorem term Q_i(n_0 - 1) 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, which is what a nonintegral degree of multiprogramming (a time-average over a measurement window) calls for.

Unlike pfqn_dnc this accepts a think time, since the delay enters the recursion and not a partial-fraction continuation. It is single-class: the multiclass recursion has no one-dimensional step. For fractional multiclass populations use pfqn_bs, which accepts them directly.

Reference: L. W. Dowdy, K. D. Gordon, “Algorithms for Nonintegral Degrees of Multiprogramming in Closed Queuing Networks”, Performance Evaluation 4(1):19-28, 1984.

Parameters:

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

Returns:

X – Throughput at population N. Q: Mean queue lengths (M x 1). U: Utilizations (M x 1). R: Residence times (M x 1).

pfqn_mvajd(Z, N, mu, Dli, visits, options)

[X, QJD, QLI, QDELAY, SJD] = PFQN_MVAJD(Z, N, MU, DLI, VISITS, OPTIONS)

Joint-dependent name of PFQN_MVAOI: the mean-value analysis of a closed network whose station rates read the whole per-class occupancy vector.

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 (the occupancy already committed at the bottom of station i), and never inspects the structure of mu_i. This is precisely 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. Order independence (mu_i constant on each support) is a modelling restriction, not an algorithmic one, so any joint-dependent scaling eta_i(n) (sn.jdscaling) is admissible.

Unlike the AMVA joint-dependence route (solver_amvald with sn.jdscaling, 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, at the cost of walking prod_r C(N_r+K+1,K+1) states with K joint-dependent stations.

Arguments and returns are exactly those of PFQN_MVAOI.

See also PFQN_MVAOI, PFQN_NCJD, PFQN_CLWJD, PFQN_MVAOI_MARG.

pfqn_bkt(L, N, Z)

Knessl-Tier expansion with the Stirling-remainder correction (BKT).

Parameters:

L – Service demand matrix (MxR).
N – Population vector (1xR).
Z – Think time vector (1xR, default: zeros).

Returns:

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

[GN,LGN]=PFQN_BKT(L,N,Z)

PFQN_BKT Knessl-Tier asymptotic expansion corrected for the Stirling remainder that steepest descent drops in each class direction.

pfqn_kt extracts N from the generating function of G by steepest descent on: F(u) = sum_r Z_r u_r - sum_k log(1-U_k) - sum_r N_r log u_r, U = L*u.

On the single-class, 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),

that is Stirling’s approximation of log(N!) in place of log(N!) itself. The expansion is therefore ABOVE the exact value by the Stirling remainder

s(N) = log(N!) - (N*log(N) - N + log(2*pi*N)/2)
= gammaln(N+1) - (N+1/2)*log(N) + N - log(2*pi)/2,

one term per Laplaced class direction, and BKT subtracts sum_r s(N_r).

s(N) is the class-direction analogue of the constant pfqn_ble adds per station direction: s(1) = 1 - log(2*pi)/2 = 0.0811 is that same constant, and s(N) = 1/(12*N) + O(N^-2) decays with the class population, so the correction matters most on lightly populated classes and on many-class models. Truncating it at 1/(12*N) loses an order of magnitude of accuracy, so the remainder is evaluated exactly from gammaln.

Only the classes that pfqn_kt actually Laplaces are corrected: a class with no jobs is dropped by pfqn_kt, and a self-looping class (one nonzero demand and no think time) has its coefficient extracted exactly, so neither contributes.

Measured over the 1562 models of the Cas17 dataset (Zenodo 546873, sec5.3.1, sigma=100), against exact convolution, the median absolute error in log G falls from 0.083 nats for KT to 1.4e-5 nats for BKT. On sec5.3.2 restricted to the models with a reference of known quality it falls from 0.64 to 0.02.

Input: L : MxR demand matrix. L(i,r) is the demand of class-r at queue i N : 1xR population vector. N(r) is the number of jobs in class r Z : 1xR think time vector. Z(r) is the total think time of class r

Output: Gn : estimated normalizing constant lGn: logarithm of Gn. If Gn exceeds the floating-point range, only lGn

will be correctly estimated.

References: C. Knessl, C. Tier. Asymptotic expansions for large closed queueing networks with multiple job classes. IEEE Trans. Computers, 41(4):480-488, 1992. G. Casale. Accelerating performance inference over closed systems by asymptotic methods. ACM SIGMETRICS 2017.

pfqn_hst(L, N, Z, ist)

Operational sensitivity of throughput to homogeneous-service-time (HST) violations, and the constrained worst case (Suri 1983).

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

The HST assumption states that the mean service time at station i does not depend on the queue length there. Suri (1983) perturbs it to S_i(n) = S_i (1 + a_n), one relative deviation a_n per queue-length level n, 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 the station utilization and the marginals taken from the product-form solution, 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 the deviations are not free: an operationally consistent perturbation must leave the observed mean service time unchanged, sum_n p_n a_n = 0 with p_n = P(n_i = n). 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 the ratio c_n/p_n exceeds a threshold, with at most one fractional coordinate. On the paper’s Figure 1 system it collapses 0.831 d to 0.102 d.

Reference: R. Suri, “Robustness of Queuing Network Formulas”, JACM 30(3):564-594, 1983 (eq. 3.11, Lemma 3.1, problem (P1)).

Parameters:

L – Service demand vector (M x 1) of the queueing stations.
N – Population (nonnegative integer scalar).
Z – Think time (scalar, default 0).
ist – Station index the HST perturbation is applied to (default: the bottleneck, argmax L).

Returns:

sens –

Struct with fields:: station the station index analysed; X product-form throughput; U utilization of that station; Q mean queue length there; Pgeq P(n_i >= k), k = 0..N; p P(n_i = k), k = 0..N; c sensitivity coefficients c_k, k = 1..N (eq. 3.11); total sum_k |c_k|, the unconstrained certificate per unit d; worst the (P1) optimum per unit d; astar the worst-case deviation profile a_k/d, k = 1..N.

pfqn_harel_bounds(rho, N, Z, maxUB)

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

Parameters:

rho – Relative utilizations (k x 1), all strictly positive.
N – Closed population, at least 1.
Z – Think time; must be zero.
maxUB – Largest extrapolation point; defaults to min(N,7).

Returns:

LB – Throughput lower bound at population N. UB: Upper bounds UB(n) for n = 2..maxUB (UB(1) unset). TH: Exact throughput TH(n) at population n, n = 1..maxUB.

[LB, UB, TH] = PFQN_HAREL_BOUNDS(RHO, N, Z, MAXUB)

Harel-Namn-Sturm throughput bounds of a single-class closed network of load-independent queues with relative utilizations rho.

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 already in SOLVER_BA: ‘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; that identity is the oracle the implementation is tested against. G is evaluated by the Newton-Girard recurrence n G(n) = sum_{i=1..n} A_i G(n-i) rather than from expanded polynomials in the power sums; the n <= 7 ceiling on the extrapolation point is kept from the reference implementation.

A nonzero think time is REFUSED rather than folded in: the bounds are derived for a network with no terminal population, so silently dropping Z would return a bound that does not bound.

See also PFQN_HAREL_LB, PFQN_HAREL_UB, PFQN_CA, PFQN_MVA, SOLVER_BA_ANALYZER.

pfqn_busyp_multiclass(alpha, mu, P, N, subnet, n, gamma, phi, tol, jobclass)

Mean busy period of order n for a subnetwork of a multichain network.

Parameters:

alpha – Relative arrival rates (JxR), one column per chain.
mu – Service rates (JxR), the chain-r rate at node j.
P – Routing matrix (JxJ), or a 1xR cell of per-chain matrices.
N – Population per chain (1xR), Inf entries for an open chain.
subnet – Indexes of the nodes forming the subnetwork.
n – Busy period order(s), 1 <= n <= sum(N).
gamma – External arrival rates (JxR), empty for a closed network.
phi – Load-dependent scaling (JxK), dimensionless; empty = single server.
tol – Relative tolerance of the open-network tail truncation.
jobclass – Chain whose own jobs are counted, or -1 to count every chain.

Returns:

b – Mean busy period duration(s), same size as n. lG: Log normalizing constants of the subnetwork over the lattice. lH: Log normalizing constants of the complement over the lattice.

[B,LG,LH] = PFQN_BUSYP_MULTICLASS(ALPHA, MU, P, N, SUBNET, N, GAMMA, PHI, TOL, JOBCLASS)

Multichain generalization of PFQN_BUSYP. The busy period of order n for a set of nodes I is the time from the instant a job entering I finds n-1 jobs in it up to the next instant when fewer than n remain, counting jobs of EVERY chain.

H. Daduna, “Busy Periods for Subnetworks in Stochastic Networks: Mean Value Analysis”, J. ACM 35(3), 1988, states Theorems 1 and 3 for a single chain and notes in Section 5 that they carry over to the whole product-form class. The proof uses only that the stationary law is product form and that the busy period is Keilson’s mean ergodic sojourn time on a level set, neither of which is single-chain, so replacing the scalar population by a per-chain vector m gives, for a closed network,

sum_{m : |m| >= n} G_I(m) H(N-m)

b(n,I) = ————————————————–
sum_{m : |m| = n-1} G_I(m) sum_r A_r(I) H(N-m-e_r)

with G_I and H the normalizing constants of the subnetwork and of its complement at a population VECTOR, and A_r(I) the chain-r arrival flow into I. The denominator is the exact chain-r flow across the cut: a chain-r departure from the complement at population k occurs at rate alpha_ir H(k-e_r)/H(k), and the H(k) cancels the state weight. At R=1 the inner sum holds the single term m=n-1 and H(N-m-e_1)=H(N-n), so the expression collapses to Theorem 1 exactly.

THE OPEN CASE NEEDS NO LATTICE. In an open product-form network the stations are independent and the total occupancy of node i is a function of the AGGREGATE load sum_r alpha_ir/mu_ir alone: summing the station function over the compositions of t collapses the multinomial to (sum_r rho_ir)^t. The open multichain problem is therefore the single-chain problem on aggregated demands, and this routine reduces it to PFQN_BUSYP rather than repeating it.

PER CLASS. With JOBCLASS = r the level set becomes {m_r >= n}, the jobs of chain r alone. Only chain-r arrivals move that level, so the flow sum loses its sum over r and the same two lattices serve every class.

A MIXED MODEL keeps the closed lattice with its OPEN dimensions TRUNCATED. The closed chains are conserved between the subnetwork and its complement, the open ones are not: the complement’s open count is free, so its open dimensions are summed out (Hbar below) and no e_r shift applies to an open chain, removing one job from an unbounded dimension leaving the same sum. The truncation is grown until the answer stops moving, and is the only approximation in that branch.

For a per-class query on a CLOSED chain of a mixed model there is an exact shortcut: marginalizing the open chains leaves a closed network with the demands deflated by 1/(1-rho_i^open). Deflate the RATES and not the visits, since A_r is built from the visit ratios and folding the deflation into alpha would scale the flow constant too.

pfqn_aghq(L, N, Z, q)

Adaptive Gauss-Hermite quadrature for the normalizing constant.

Parameters:

L – Service demand matrix (MxR).
N – Population vector (1xR).
Z – Think time vector (1xR).
q – Nodes per simplex direction (default: 3). q=1 reproduces pfqn_le.

Returns:

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

[GN,LGN]=PFQN_AGHQ(L,N,Z,Q)

pfqn_sqni(N, L, Z)

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

Parameters:

N – Population vector.
L – Service demand vector.
Z – Think time vector.

Returns:

Q – Mean queue lengths. U: Utilization. X: System throughput.

pfqn_sens_respt_validate()

Validation harness for pfqn_sens_respt, the FCFS sojourn-time moment analysis of Strelen (1990), Theorem 4.1.

Five references: A. brute-force enumeration. This is the strongest check because it shares none of Theorem 4.1’s algebra. The equilibrium product form is enumerated to get the exact arrival-theorem marginals p_i(j, N-1_l); the sojourn time conditioned on finding j jobs is known in closed form (Exp(mu) if j < b, otherwise an Erlang(j-b+1, b*mu) queueing delay plus an Exp(mu) service), so its moments are formed directly and mixed over j. Neither the coefficients a_(t,tau)(0) nor the recursion (4.2) enter, so the agreement tests both. B. the published table of Example 3.4 (continued) of the reference, which prints E(W_i) and sigma^2_(W_i) for the Kobayashi model. C. the internal identity W(i,l) = w_i(l)/V(i,l): the t = 1 case of (4.5) must reproduce the MVA residence time divided by the visit ratio, which is a completely different expression. D. pfqn_mva for the base measures in the single-server case. E. the single-job network, where an arriving job always finds an empty station, so W is exactly Exp(mu) and every moment is known in closed form.

Reference: J. C. Strelen, “Moment Analysis for Closed Queuing Networks and its Linearizer”, Performance Evaluation 11:127-142, 1990.

pfqn_sens_mvaldmx(lambda, D, N, Z, mu, S)

Exact queue-length variances and covariances for mixed open/closed product-form queueing networks with limited load dependence.

Exact second moments (variances and covariances) of the queue lengths of a mixed open/closed product-form queueing 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.

Reference: I. F. Akyildiz and J. C. Strelen, “Moment Analysis for Load-Dependent Mixed Product Form Queueing Networks”, IEEE Trans. Communications 39(6):828-832, 1991. The closed load-independent case reduces to E. de Souza e Silva and R. R. Muntz, IEEE Trans. Computers 37(9):1125-1129, 1988, which pfqn_sens_mva implements directly.

Parameters:

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

Returns:

mom –

A struct with the base measures and their second moments:

.X (1 x R), .Q (M x R), .U (M x R), .R (M x R) base measures, identical
to pfqn_mvaldmx(lambda,D,N,Z,mu,S).

.QCov (M x R x R) QCov(i,r,s) = Cov[n(i,r),n(i,s)], same-station. .QCovFull (M x R x M x R) QCovFull(i,r,j,s) = Cov[n(i,r),n(j,s)]. .QVar (M x R) QVar(i,r) = Var[n(i,r)]. .QTotVar (M x 1) QTotVar(i) = Var[sum_r n(i,r)]. .QCovAsym (scalar) 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.

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.

pfqn_sens_mva(L, N, Z, mi)

Exact per-station queue-length variances and covariances for closed product-form queueing networks, computed by an MVA-like moment recursion that does not require the full sensitivity Jacobian.

Exact second moments (variances and per-station covariances) of the queue lengths of a closed product-form (BCMP) queueing 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,

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 the reference below, 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) of the reference, 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.

Reference: E. de Souza e Silva and R. R. Muntz, “Simple Relationships Among Moments of Queue Lengths in Product Form Queueing Networks”, IEEE Trans. Computers 37(9):1125-1129, 1988 (Theorems 1-2 and Corollary 1). The underlying identity Cov[n(i,k),n(j,v)] = theta(i,k) * dQ(j,v)/dtheta(i,k) is the k=2 case of Theorem 1 of I. F. Akyildiz and J. C. Strelen, “Moment Analysis for Load-Dependent Mixed Product Form Queueing Networks”, IEEE Trans. Communications 39(6):828-832, 1991.

Parameters:

L – Service demand matrix (M x R), L(i,r) = visits_ir / rate_ir.
N – Population vector (1 x R).
Z – Think time vector (1 x R). Default: zeros.
mi – (Optional) Server multiplicity vector (1 x M). Default: ones.

Returns:

mom –

A struct with the base measures and their second moments:

.X (1 x R), .Q (M x R), .U (M x R), .R (M x R) base MVA measures,
identical to pfqn_mva(L,N,Z,mi).

.QCov (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 (M x R) QVar(i,r) = QCov(i,r,r) = Var[n(i,r)]. .QTotVar (M x 1) QTotVar(i) = Var[sum_r n(i,r)], the variance of the

total queue length at station i, i.e. sum_{r,s} QCov(i,r,s). This is Theorem 3 of the reference, obtained here without a capacity derivative.

.QCovAsym (scalar) 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.

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.

pfqn_sens_mom(L, N, Z, mi, groups)

Exact higher moments (up to order three) of the per-station total queue lengths of a closed product-form queueing network, by second-order differentiation of the MVA recursion.

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(1,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. pfqn_sens_mom_validate checks the per-class setting against brute force and against pfqn_sens_mva’s second moments, which it must reproduce exactly.

Reference: J. C. Strelen, “Moment Analysis for Closed Queuing Networks and its Linearizer”, Performance Evaluation 11:127-142, 1990, Theorems 2.1, 3.1, 3.2, 3.5 and equation (3.2).

Parameters:

L – Service demand matrix (M x R), L(i,r) = visits_ir / rate_ir.
N – Population vector (1 x R).
Z – Think time vector (1 x R). Default: zeros.
mi – (Optional) Server multiplicity vector (1 x M). Default: ones.
groups – (Optional) Class-to-group map (1 x R), a partition of the classes into G = max(groups) groups labelled 1..G. The moments returned are those of each group’s queue length at each station. Default: ones(1,R), i.e. one group holding every class, which is the per-station total.

Returns:

mom –

A struct with:

.X (1 x R), .Q (M x R), .U (M x R), .R (M x R) base MVA measures,
identical to pfqn_mva(L,N,Z,mi).

.m (M x G) m(i,g) = E[Q_(i,g)], the mean queue length of group g at
station i. With the default groups this is (M x 1), the station total.

.d2m (M x 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.

.Var (M x G) Var[Q_(i,g)]. .M2 (M x G) E[Q_(i,g)^2]. .M3 (M x G) E[Q_(i,g)^3]. .Skew (M x G) skewness of Q_(i,g). NaN where Var is zero. .Cov, .dm Cov((i,g),(j,g’)) = Cov[Q_(i,g),Q_(j,g’)] and the scaled

first derivative it comes from. Shaped (M x G x M x G) in general, but COLLAPSED to (M x M) in the default single-group case, where the group index carries no information and an (M x 1 x M x 1) array would only be awkward to index.

.CovAsym (scalar) raw asymmetry of Cov before symmetrization. The two
triangles are distinct expressions that must agree, so this is a live residual of the recursion; expect roundoff.

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.

pfqn_sens_linearizer(L, N, Z, tol, maxiter)

Approximate higher moments of the queue lengths of a closed product-form queueing network, by differentiating the Linearizer fixed point. Polynomial in the population, unlike the exact pfqn_sens_mom.

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]. pfqn_sens_linearizer_validate measures the error against the exact pfqn_sens_mom on models small enough for both, and asserts bands of that order rather than machine precision: this routine is an approximation and is expected to disagree with the exact answer.

Reference: J. C. Strelen, “Moment Analysis for Closed Queuing Networks and its Linearizer”, Performance Evaluation 11:127-142, 1990, Section 5, equations (5.1)-(5.8), the CORE-2 and LINEARIZER-2 algorithms. The delta definition follows the standard Chandy-Neuse Linearizer, v_i^(N-1_l’)(l) - v_i^(N)(l), which is what LINE’s pfqn_linearizer implements.

Parameters:

L – Service demand matrix (M x R), L(i,r) = visits_ir / rate_ir.
N – Population vector (1 x R).
Z – Think time vector (1 x R). Default: zeros.
tol – (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 – (Optional) Maximum CORE iterations. Default: 200.

Returns:

res –

A struct with:

.X (1 x R), .Q (M x R), .U (M x R), .W (M x R) approximate base measures. .m (M x 1) approximate E[Q_i], the total queue length at station i. .dm (M x M) dm(i,h) = x_h dm_i/dx_h, the scaled first derivative. .d2m (M x 1) d2m(i) = x_i^2 d^2m_i/dx_i^2. .Var (M x 1), .Cov (M x M), .M2 (M x 1), .M3 (M x 1), .Skew (M x 1)

the moments of (3.2), formed exactly as in pfqn_sens_mom but from the approximate derivatives.

.CovAsym (scalar) 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 (scalar) total CORE iterations performed.

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.

pfqn_schmidt_ext(D, N, S, sched)

Extended Schmidt MVA algorithm with queue-aware alpha corrections.

Parameters:

D – Service demand matrix (M x R).
N – Population vector (1 x R).
S – Number of servers per station (M x 1).
sched – Scheduling discipline per station.

Returns:

XN – System throughput. QN: Mean queue lengths. UN: Utilization. CN: Cycle times. T: Results table.

[XN,QN,UN,CN,T] = PFQN_SCHMIDT_EXT(D,N,S,SCHED)

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.

pfqn_respt_ps_moments(S, N, Z, method)

[W,W2,OUT] = PFQN_RESPT_PS_MOMENTS(S, N, Z, METHOD)

Sojourn-time moments at the processor-sharing station of the closed terminal-driven system of Mitra and Morrison (1983): 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 that paper:

‘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 4096 states and the asymptotic route otherwise.

The asymptotic route requires the normal-usage condition alpha > 0, where alpha = 1 - sum_r lambda_r/q_r with lambda_r = K(r)/Z(r) and q_r = 1/S(r), is the unutilized fraction of the CPU in the corresponding open system. Where it fails and the exact route is not affordable, the entry of W and W2 is NaN and OUT.method records ‘unavailable’; asking for ‘asymptotic’ explicitly in that regime is an error rather than a blank.

Input:: S - per-class mean service times at the PS station (1,R), positive N - per-class populations (1,R), non-negative integers Z - per-class mean think times (1,R), positive where N > 0 method - ‘auto’ (default), ‘exact’ or ‘asymptotic’
Output:: W - per-class mean sojourn times at the PS station (1,R) W2 - per-class second moments of the sojourn time (1,R) out - struct with fields method (1,R cell), c0, c1 (1,R, asymptotic route

only), alpha (1,R), nstates (1,R) and expansionParam (scalar Nexp)

A class with N(r) = 0 has no sojourn time and its entries are NaN.

Reference: D. Mitra, J. A. Morrison, “Asymptotic Expansions of Moments of the Waiting Time in Closed and Open Processor-Sharing Systems with Multiple Job Classes”, Adv. Appl. Prob. 15(4), 1983, Propositions 3 and 6.

See also: qsys_mm1_ps (the open counterpart, exact in closed form).

test_pfqn_manjunath()

TEST_PFQN_MANJUNATH Validate the Manjunath-Sikdar transform for product-form queueing networks against pfqn_ca and against brute-force enumeration.

THE TWO ORACLES DO NOT COME OUT OF THE IMPLEMENTATION. With no extra rows the transform computes the ordinary closed-network normalizing constant, for which pfqn_ca’s convolution recursion is an independent exact algorithm sharing no code; the two are algebraic identities for the same sum, so agreement to 1e-13 is the correct expectation and not a tolerance chosen to pass. With extra rows pfqn_ca has nothing to say, and the oracle becomes bcmp_enum below, which sums the BCMP product form over the enumerated state space and applies each row by direct comparison – the very enumeration the transform exists to avoid, so a coefficient-domain defect cannot hide behind a shared traversal.

pfqn_pbk(L, N, Z, k)

Iterative PB(k) proportional bounds (Eager-Sevcik) for single-class closed networks with delay.

PB(k): the Eager-Sevcik proportional (performance) bound at iteration count k, i.e. the level-k member of the Performance Bound Hierarchy computed at populations N, N-1, …, N-k (Casale-Muntz- Serazzi 2008 cite this as the iterative extension used in Tables 5,7). Reduces to the noniterative asymptotic bound at k=0 and tightens to exact as k->N. Backed by the validated PBH recursion (pfqn_pbh).

Parameters:

L – Service demand vector (M x 1).
N – Population (scalar).
Z – Think time (scalar, default 0).
k – Iteration count >= 0 (default 1).

Returns:

Xlo – Lower throughput bound. Xhi: Upper throughput bound.

pfqn_pas_nc(Z, N, mu, prec, options)

[G, LG] = PFQN_PAS_NC(Z, N, MU, PREC, OPTIONS)

Normalizing constant of a closed pass-and-swap (P&S) queueing network that comprises a single aggregated infinite-server (delay) node and an arbitrary number of order-independent (OI) / P&S stations, restricted to the recurrent communicating class selected by the placement order PREC.

With a non-empty swap graph the ordered-state chain is reducible (Comte & 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 that per-class constant 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 and does not collapse onto the count lattice. For the plain OI case (empty PREC, every ordering feasible) the count lattice does suffice and PFQN_NCOI computes the same G at far lower cost; use this routine only when a placement order is present, or as a microstate reference.

Method. Build station M’s chain head-first: appending class r at chain position k = sum(occ)+1 is admissible iff the placement order allows it (no class already placed at that station must come after r), and contributes the reciprocal OI prefix rate 1/mu_M(occ+e_r); the chain may be finalized only when occ is a placement-order ideal at full multiplicity, whereupon the recursion moves to station M-1. When every P&S station has been peeled the residual population sits at the delay node with the multinomial weight prod_r Z_r^{N_r}/N_r!. With PREC empty this is exactly the balanced-fairness convolution of Bonald & Proutiere (2003) evaluated ordering by ordering,

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

COST. The recursion visits one node per feasible ordered prefix, so with an empty PREC the node count is sum_{b<=N} C(|b|+M-1,M-1) * |b|!/prod_r b_r!, i.e. factorial in the total population sum(N). A placement order prunes the orderings (a total order leaves a single one per count split), which is what makes the microstate walk affordable in the P&S case.

Parameters:

Z - (1 x R) – Z(r) = 1/sigma_r for a delay with per-class rate sigma_r.
N - (1 x R)
mu - cell array {1 x M} of function handles, one per P&S station. Each – mu{m}(n) returns the total service rate of station m given the per-class occupancy (count) vector n (1 x R). For an OI station the rate depends only on the support of n (which classes are present), i.e. the sum of the capacities of the compatible servers. May be empty to model a pure delay network.
prec - placement order. Either a cell {1 x M} of (R x R) – matrices, one per station, or a single (R x 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). Empty or all-zero means no order: every ordering is feasible and G is the plain OI constant. NOTE the orientation: around a cycle 1->2->…->M->1 the chain of each downstream station is traversed in the opposite direction, so the downstream stations take the TRANSPOSE of the upstream order (prec = {P, P’} for a two-station cycle). Passing the same P to both stations of a cycle silently returns a smaller, wrong G.
options - solver options (optional, currently unused; accepted for – signature parity with the other pfqn_* routines).

Returns:

G - Normalizing constant G_C of the communicating class. lG - log(G_C).

Example (two P&S stations in a cycle with swap graph SWAP, no delay):: H = pas_swap2order(swap, {@(c) mu1, @(c) mu2}); P = pas_placement(H); G = pfqn_pas_nc([], N, {@(n) mu1, @(n) mu2}, {P, P’});

See also PFQN_NCOI, PFQN_PAS_IS, PAS_PLACEMENT, PAS_SWAP2ORDER.

pfqn_mva_ilock(L, N, Z, mi, IL)

Exact MVA recursion carrying the interlocked-flow correction.

Parameters:

L – Service demand matrix (M x R).
N – Population vector (1 x R).
Z – Think time vector (1 x R).
mi – (Optional) Server multiplicity vector (1 x M). Default: single servers.
IL – Interlock matrix (R x R). IL(r,s) is the share of the class-s queue that a class-r arrival cannot see, because that work was itself caused by the class-r request (Franks 1999, Ch. 4, Eq. 4.7).

Returns:

XN – System throughput (1 x R). QN: Mean queue length (M x R). UN: Utilization (M x R). CN: Residence time (M x R). lGN: Always NaN: the interlock leaves the model outside product form.

[XN,QN,UN,CN,LGN] = PFQN_MVA_ILOCK(L,N,Z,MI,IL)

Closed single-server models only. The correction replaces the arrival theorem term Q(n-1_s,i) by a per-class weighted sum sum_r ILw(s,r)*Qc(n-1_s,i,r), so the recursion has to carry per-class queue lengths that PFQN_MVA does not need. Pass an empty IL to PFQN_MVA instead.

The discounted arrival-instant queue is floored at the in-service component, as in lqns MVA::queueOnly_adjusted, so the correction damps itself out as a station saturates. That is a self-limiting guard, NOT a hard capacity test: sum_s XN(s)*L(i,s) <= mi(i) is still asserted nowhere. See git show 8bad654e7:_kb/log.md.

pfqn_panacea(L, N, Z, terms)

PANACEA (PAth-based Normal Approximation for Closed networks Estimation Algorithm).

Parameters:

L – Service demand matrix.
N – Population vector.
Z – Think time vector.
terms – 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:

Gn – Normalizing constant. lGn: Logarithm of normalizing constant.

[GN,LGN]=PFQN_PANACEA(L,N,Z,TERMS)

pfqn_oi_insvc(oirate, N, options)

[G, XI, PHI] = PFQN_OI_INSVC(OIRATE, N, OPTIONS)

Conditional mean number of IN-SERVICE jobs per class at an order-independent (OI) station, as a function of the per-class count vector n. This is the quantity underlying the LINE utilization convention at OI 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 that 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.toMarginal, PAS branch) and by LDES.

Because mu is permutation-invariant (a function of the count vector), the unnormalized weight of an ordering c of the multiset n factorizes over its prefixes as w(c) = prod_{p=1}^{|n|} 1/mu(n(c_1..c_p)), and the OI balance function is Phi(n) = sum_{orderings c of n} w(c), which obeys the standard balanced-fairness recursion (condition on the tail element c_{|n|}):

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

Conditioning the same way on the tail element 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 the sir-weighted balance Xi_r(n) = sum_{orderings c of n} 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 prod_r x_r^{n_r}, 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’s in-service mean then follows from the count marginal pM by 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 handle mu (n) – per-class count vector n (1 x R). mu(0) is taken as 0.
N - (1 x R)
options - solver options (optional, currently unused)

Returns:

g - (prod(N+1) x 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) x R) table with the sir-weighted balance Xi_r(n). Phi - (prod(N+1) x 1) table with the OI balance function Phi(n).

See also PFQN_OI_FNC, PFQN_NCOI, PFQN_MVAOI, PFQN_MVAOI_MARG.

pfqn_nrp(L, N, Z, alpha, options)

Normalizing constant via Norlund-Rice Probit (NRP) approximation.

Parameters:

L – Service demand matrix.
N – Population vector.
Z – Think time vector.
alpha – Load-dependent rate matrix.
options – Solver options.

Returns:

lG – Logarithm of normalizing constant.

pfqn_momlin(L, N, Z, tol, maxiter)

Moment linearizer: approximate mean queue lengths and their second moments (variance / covariance) for large closed product-form queueing networks.

Approximate first and second queue-length moments of a closed product-form network, scalable to large populations and many classes where exact MVA (exponential in the number of classes) and CoMoM (single-station) are infeasible.

Means are obtained from the Schweitzer-Bard AMVA fixed point. Second moments use the exact product-form identity Cov[n_{i,r},n_{j,s}] = D_{j,s} dQ_{i,r}/dD_{j,s}, with the demand derivatives obtained by analytically linearizing the AMVA fixed point (a “moment linearizer” in the sense of Strelen and Akyildiz, here developed per class). Both moments carry the AMVA approximation error and become exact only in the limits where Schweitzer-Bard is exact; for exact moments on tractable models use pfqn_sens (differentiated MVA / CoMoM).

Parameters:

L – Service demand matrix (M x R).
N – Closed population vector (1 x R).
Z – Think time vector (1 x R). Default: zeros.
tol – Convergence tolerance on the queue-length fixed point. Default 1e-8.
maxiter – Maximum iterations. Default 1000.

Returns:

Q – Mean queue length (M x R). X: Throughput per class (1 x R). U: Utilization (M x R). R: Residence time (M x R). QVar: Queue-length variance (M x R). QCov: Queue-length covariance tensor (M x R x M x R). dQ: Demand-derivative tensor, dQ(i,r,j,s) = dQ_{i,r}/dD_{j,s} (M x R x M x R).

pfqn_pam(L, N, Z, variant)

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

1. 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’ seed, then the last MVA step (eqs. 2.29-2.34) ‘pami’ PAMB, then throughputs scaled down wherever a centre would be driven past full utilization (eqs. 2.35 and step 4) ‘pamt’ seed at N - 1_i - 1_j, then the last TWO MVA steps, then the PAMI utilization capping (eqs. 2.36-2.40)

The seed spreads the whole class population over the queueing centres and ignores Z, exactly as published: PAM buys speed, not accuracy.

Parameters:

L – Service demand matrix (stations x classes).
N – Population vector.
Z – Think time vector.
variant – ‘pamb’ (default), ‘pami’ or ‘pamt’.

Returns:

XN – System throughput. QN: Mean queue lengths. UN: Utilization. RN: Residence times.

pfqn_looping(L, N, Z, tol, maxiter)

Eager Looping approximate MVA.

1. 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 (pfqn_pbh) starts from, so it carries a pair of bounds rather than a single fixed point.

It is built on the convolution identity of Zahorjan (1980)

Q_jk(N - 1_c) = [X_j^{+k}(N - 1_c) / X_j(N)] Q_jk(N),

with X_j^{+k}(N - 1_c) estimated from the level-0 multiple-class PBH upper bound B_j and X_j(N) from the optimistic response time R_j^(opt). A HEAP H_j is the class-j congestion that the current queue-length lower bounds have not yet 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, which are the largest and smallest delays one customer can inflict.

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.

Looping is a BOUNDING algorithm, not a point estimator: it returns the bracket X^(pess) of eq. (2.23) and its optimistic counterpart N_c/(Z_c + R_c^(opt)) built on eq. (2.24). It is reached through SolverBA as ‘looping.lower’/’looping.upper’.

Parameters:

L – Service demand matrix (stations x classes).
N – Population vector.
Z – Think time vector.
tol – Convergence tolerance (default: 1e-6).
maxiter – Maximum number of iterations (default: 1000).

Returns:

Xlo – Pessimistic (lower) throughput bound, one entry per class. Xup: Optimistic (upper) throughput bound, one entry per class. QN: Mean queue lengths on the pessimistic side, eq. (2.25). RN: Residence times, eq. (2.20). it: Number of iterations performed.

pfqn_manjunath(L, N, Z, A, b, sense)

Exact normalizing constant of a closed multiclass product-form network whose state space carries arbitrary linear integer constraints

This is the queueing-network half of the transform technique of Manjunath and Sikdar, of which lossn_manjunath is the loss-network half. The two solve the same problem – sum a product form over an irregular integer state space – but from opposite ends of the paper: lossn_manjunath implements Section 2.2, a set of ‘<=’ constraints over Poisson terms nu^n/n!, while this routine implements Section 3 together with Section 5.3, a MIXED set of ‘=’, ‘<=’ and ‘>’ constraints over the BCMP terms, where the population constraint of a closed network is itself one of the equalities.

Parameters:

L – Service demand of class r at queueing station i (MxR)
N – Population of class r (1xR nonnegative integers)
Z – Think time of class r at delay station k (MzxR), default zeros
A – Extra constraint coefficients on n(:) (Jx((M+Mz)*R), nonnegative integers), default empty
b – Extra constraint right-hand sides (Jx1 integers), default empty
sense – Row senses, char vector of ‘E’ (=), ‘L’ (<=), ‘G’ (>), default all ‘L’

Returns:

G – Normalizing constant lG: Logarithm of the normalizing constant peak: Peak number of live series coefficients, the realised cost stats: Per-class decomposition (1xR fields Q, X, U, think, blocked, delay). Needs one queueing station and one delay station outside the region

The model:

M queueing stations (FCFS, PS or LCFS, rows of L) and Mz delay stations (rows of Z) serve R closed classes with populations N. Writing n_ir for the number of class r jobs at station i and n_i = sum_r n_ir, the BCMP product form of Baskett-Chandy-Muntz-Palacios is

p(n) = (1/G) prod_{i queueing} n_i! prod_r L_ir^{n_ir}/n_ir!
prod_{i delay} prod_r Z_ir^{n_ir}/n_ir!

and G is its sum over the admissible set. Every state obeys the R population equalities sum_i n_ir = N_r; on top of those the caller may impose any number of further rows

sum_{i,r} A(j, i + S(r-1)) n_ir {=, <=, >} b(j), S = M + Mz,

i.e. A acts on n(:), the (M+Mz)-by-R occupancy matrix read column by column with the queueing stations first. With no extra rows the routine returns exactly the normalizing constant of pfqn_ca, which is the parity oracle used by the tests; with extra rows it answers a question no other routine in the pfqn family can, because the convolution and MVA recursions are built around the population constraint alone and have nowhere to carry a second one.

Why the generating function is a product, and why the n_i! disappears:

Marking class r by z_r and constraint row j by y_j, and abbreviating the monomial that one class r job at station i contributes as

u_ir = z_r prod_j y_j^{A(j, i + S(r-1))},

the sum over the occupancies of a single QUEUEING station is, by the multinomial theorem,

sum_{n_i.} n_i! prod_r (L_ir u_ir)^{n_ir}/n_ir!
= sum_k (sum_r L_ir u_ir)^k = 1 / (1 - sum_r L_ir u_ir),

so the n_i! that couples the classes at a queueing station is exactly what turns the station’s factor from an exponential into a geometric one. The paper reaches the same place through the Euler integral n! = int_0^inf e^{-t} t^n dt (Eqns 16-18), which is that geometric series evaluated; the closed form is used here because there is then no quadrature to discretize. A DELAY station has no n_i! and keeps its exponential, prod_r exp(Z_ir u_ir). Hence

F(z,y) = prod_{i=1}^{M} 1/(1 - sum_r L_ir u_ir)
prod_{k=1}^{Mz} prod_r exp(Z_kr u_kr)

and G is read off F as a coefficient: degree exactly N_r in z_r for every class, and for row j the degree dictated by its sense – exactly b_j for ‘=’, the sum of degrees 0 … b_j for ‘<=’ (the multiplier (y^{b+1}-1)/(y-1) of the paper’s Eqn 5, whose residue is that partial sum), and the complement of the latter for ‘>’ (Eqn 6).

Why this is a coefficient computation and not a quadrature:

The contour integrals of Eqn 9 all have their only pole at the origin, of order one more than the right-hand side, so each is a residue and hence a Taylor coefficient. The routine therefore never integrates: it carries F as a multivariate power series truncated at degree N_r in z_r and b_j in y_j. Truncation is exact because A is nonnegative – no monomial above the cut can be brought back down by a later factor.

Each queueing station is applied by SOLVING (1 - sum_r L_ir u_ir) x = ser rather than by expanding the geometric series, which is what keeps the cost at one pass. Every monomial of the operator carries z_r to a strictly higher power, so sweeping the lattice in increasing total class degree lets each coefficient read only coefficients already final: a Gauss-Seidel sweep whose result is the exact solve, not an iterate. A delay station has no such recurrence and is convolved with exp term by term, which is where the extra factor of the population in the cost comes from.

The elimination order is the memory bound:

Variable y_j is created when the first station its row touches is multiplied in and is discharged immediately after the last one, so peak memory is prod_r (N_r+1) times the product of (b_j+1) over the SIMULTANEOUSLY LIVE rows, not over all rows. A row that constrains one station therefore costs essentially nothing. The class dimensions are live throughout, so prod_r (N_r+1) is a floor on the cost – the same lattice pfqn_ca walks.

Scope:

Load-dependent and multiserver stations are NOT covered: their per-station term is not geometric, and while the paper admits an arbitrary f_i(n_i) in the single-class case (Section 2), the multiclass n_i! coupling used above then breaks. Use pfqn_gld or pfqn_conwayms for those. A and b must be integer valued and A nonnegative, since the residue argument counts whole units; a fractional entry is refused rather than rounded.

The per-class decomposition, and where the blocked jobs sit:

Asking for the fourth output STATS returns the whole solution of the constrained network, not just its normalizing constant. It requires the ONE configuration in which the truncated product form is the EXACT stationary law: a single queueing station inside the region and a SINGLE DELAY STATION OUTSIDE IT. Anything else is refused by name rather than answered wrongly.

Why that configuration is not a convenience: with one queueing station the state is the queue occupancy alone (the delay holds the complement) and every transition moves one job of one class by one unit, so the chain is a multidimensional birth-death process. That process is reversible, and Kelly’s truncation theorem then applies verbatim – restricting it to the coordinate-convex set A n <= b and renormalizing gives exactly the truncated product form. Add a second queueing station and the delay -> q1 -> q2 -> delay cycle destroys reversibility; truncation no longer preserves the product form, measured at 131% relative error on the stationary law of a 2-class, N = [2 2] instance. G and lG remain correct as a sum over the admissible set in every configuration; only the metrics are withheld.

Everything follows from two ratios of normalizing constants, both taken in the log domain so the internal rescaling cancels without being reconstructed:

X_r = G(N - e_r ; b - A(:,qcol_r)) / G(N ; b) P(n_qr = k) = G(N ; b, with the added row n_qr = k) / G(N ; b)

The first is the loss network’s g(C - A e_r) in another guise: removing one class r job from the queue leaves a state whose admission rule is shifted by that job’s own requirement column. The second is what an ‘=’ row is for.

A refused admission is a DELETED transition, so a blocked job never leaves the delay, and since the think time is exponential a held job is indistinguishable from one still thinking. The delay population carries both, and Little’s law separates them:

delay_r = N_r - Q_r (everything not at the queue) think_r = X_r Z_r (genuinely thinking) blocked_r = delay_r - think_r (held at the delay by the constraint)

This is NOT the WAITQ rule of SolverSSA/SolverCTMC/JMT, which moves a refused job out of the delay into a per-region FIFO counted at no station. That is a different chain, and on the reference instance below its class throughputs differ by 15%.

Examples

[G, lG] = pfqn_manjunath(L, N)
[G, lG] = pfqn_manjunath(L, N, Z)
[G, lG, peak] = pfqn_manjunath(L, N, Z, A, b, sense)
[G, lG, peak, stats] = pfqn_manjunath(L, N, Z, A, b, sense)

stats fields:: <table> <tr><th>Field<th>Description <tr><td>Q<td>Mean class r jobs at the queueing station <tr><td>X<td>Class r cycle throughput <tr><td>U<td>Class r utilization of the queueing station (X_r times its demand) <tr><td>think<td>Class r jobs genuinely thinking, X_r Z_r <tr><td>blocked<td>Class r jobs held at the delay by the constraint <tr><td>delay<td>Class r jobs at the delay, think + blocked </table>

References

D. Manjunath and B. Sikdar, Integral Expressions for the Numerical Evaluation of Product Form Expressions Over Irregular Multidimensional Integer Spaces. Sections 3 and 5.3.

pfqn_bklc(L, N, Z, method, tol, maxiter)

Birman-Kogan load concealment algorithm: multichain solved as single chain problems.

Birman and Kogan (Stochastic Models 8(3):543-563, 1992), Algorithm 2. The saddle point analysis of Corollary 2 shows that chain l may be solved on its own provided every station is 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; the paper’s contribution is the asymptotic argument that says when it is exact, and the extension to state dependent servers.

The single chain subproblem is solved either exactly (MVA) or by the uniform expansion of PFQN_BKUE.

Parameters:

L – Service demand matrix (stations x classes).
N – Population vector (1 x classes).
Z – Think time vector (default: zeros).
method – Single chain solver, ‘mva’ (default) or ‘ue’.
tol – Convergence tolerance on the throughputs (default: 1e-10).
maxiter – Maximum number of sweeps (default: 1000).

Returns:

X – Chain throughputs. Q: Mean queue lengths (stations x classes). U: Utilizations (stations x classes). it: Number of sweeps performed.

[X,Q,U,IT] = PFQN_BKLC(L,N,Z,METHOD,TOL,MAXITER)

Load concealment reduction of a multichain closed network.

pfqn_cyclet_ofree(v, mu, N, path, tset, options)

[f, F, MOM, OUT] = PFQN_CYCLET_OFREE(V, MU, N, PATH, TSET, OPTIONS)

Exact passage-time density, CDF and moments along an OVERTAKE-FREE PATH of a closed single-chain tree-like product-form network with population N.

V (1,M) visit ratios, MU (1,M) service rates, N the population, PATH the node list z = (z_1,…,z_m) of the overtake-free path with z_1 the root, TSET the time grid. PATH may instead be a cell array of paths, in which case OPTIONS.pathprob weights them and the outputs are the mixture; that is how a cycle time is assembled when the root branches.

Reference: P. G. Harrison and W. J. Knottenbelt, “Passage Time Distributions in Large Markov Chains”, 2002, Sec. 7.1, Theorems 1 and 2, after P. G. Harrison, J. Appl. Prob. 27, 1990 and H. Duduna, Adv. Appl. Prob. 14, 1982. The underlying sojourn-time result for overtake-free paths is F. Kelly and P. Pollett, Adv. Appl. Prob. 15, 1983.

THE ONE FACT THAT MAKES ALL THREE ROUTES WORK. Conditional on the path,

T | z = sum_{j in z} Erlang(u_{z_j} + 1, mu_{z_j})

with u distributed as the network’s equilibrium population vector AT N-1 (the arrival theorem). Hence the transform of Theorem 1 collapses to

L(s|z) = prod_{j in z} mu_j/(s+mu_j) * G(y(s), N-1) / G(x, N-1)

where x_i = v_i/mu_i and y_i(s) = x_i mu_i/(s+mu_i) on the path, x_i off it. One Buzen convolution per value of s.

OPTIONS.method selects the density route:

‘auto’ (default) ‘exact’ when the path rates are separated, else ‘lt’ ‘exact’ Theorem 2 in closed form; REQUIRES DISTINCT RATES on the path,

since its partial fractions divide by prod_{i~=j}(mu_i - mu_j)

‘lt’ the transform above inverted through api/lti: (OPTIONS.lti_method, default ‘euler’)

MOMENTS ARE NEVER TAKEN FROM THE DENSITY. They come from running the same Buzen convolution in the ring of truncated power series in s, so they are exact to machine precision, are unaffected by the time grid, and stay valid when the rates coincide and Theorem 2 does not apply.

NOTE ON THE PAPER. The inner sum of Theorem 2 reads (v_j t)^(c-i)/(c-i)! and that is CORRECT as printed, however odd the visit ratio looks against a time: substituting the service rate instead returns negative densities. Verified against a direct mixture-of-Erlangs oracle to 1e-15, and at the paper’s own N = 18 example against the transform route to 1e-11.

pfqn_lcp(L, N, Z, tol, maxiter, QN0, type)

Bard Large Customer Population (LCP) approximate MVA.

Bard, “Some extensions to multiclass queueing network analysis”, in Performance of Computer Systems, North-Holland, 1979. The first approximate MVA algorithm: it estimates the arrival-instant queue length 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.

Parameters:

L – Service demand matrix (stations x classes).
N – Population vector.
Z – Think time vector.
tol – Tolerance for convergence.
maxiter – Maximum number of iterations.
QN0 – Initial guess for queue lengths.
type – Scheduling strategy type (default: PS).

Returns:

XN – System throughput. QN: Mean queue lengths. UN: Utilization. RN: Residence times. it: Number of iterations performed.

pfqn_scb(L, N)

Dowdy-Carlson-Krantz-Tripathi (1992) single-class bounds on the performance of the multiclass system a single-class model aggregates.

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. L. W. Dowdy, B. M. Carlson, A. T. Krantz, S. K. Tripathi, “Single-Class Bounds of Multi-Class Queuing Networks”, J. ACM 39(1):188-213, 1992.

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 (demands weighted by the relative class throughputs) 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, hence X_R <= X_1*(N+m-1)/N. The single server capacity U_k,R <= 1 caps the same ratio at 1/(X_1*max(L)), and that cap is tight on the paper’s own worst case (m saturated devices, where D_k,1 = 1/X_R for every k), so both are applied.

Parameters:

L – Service demand vector of the single-class model (K x 1), queueing stations only. Delay stations are not admitted: Theorem 3 rests on the delay-free balanced-network throughput N/((N+m-1)D).
N – Total population (scalar, N >= 1).

Returns:

Xlo – Lower bound on multiclass total throughput X_R (= exact X_1). Xhi: Upper bound on X_R. Ulo: Lower bound on the per-device utilizations U_k,R (K x 1). Uhi: Upper bound on U_k,R (K x 1).

pfqn_bk(L, N, Z)

Birman-Kogan saddle point normalizing constant with bottleneck detection.

Birman and Kogan (Stochastic Models 8(3):543-563, 1992) evaluate the multichain partition function by the saddle point method applied to the Cauchy inversion of its generating function. 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, which is where the residue rather than the saddle carries the mass. The remaining stations are the paper’s large groups of identical stations and are exponentiated.

Parameters:

L – Service demand matrix (stations x classes).
N – Population vector (1 x classes).
Z – Think time vector (default: zeros).

Returns:

G – Normalizing constant. lG: Logarithm of the normalizing constant. X: Chain throughputs (the saddle point coordinates). U: Utilizations (stations x classes). A: Chains whose dedicated station is not saturated (eq. 29). B: Chains whose dedicated station is a bottleneck (eq. 30).

[G,LG,X,U,A,B] = PFQN_BK(L,N,Z)

Birman-Kogan saddle point expansion of the normalizing constant.

pfqn_xia(L, N, s)

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

Parameters:

L – Service demand vector (M x 1).
N – Closed population (scalar).
s – Server counts (M x 1).

Returns:

lGasy – Logarithm of the approximate normalizing constant.

LGASY = PFQN_XIA(L, N, S)

Xia’s asymptotic approximation of the normalizing constant of a load-dependent (multiserver) single-class closed network.

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

c being 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 complex here and NaN in the reference; NaN is returned so the refusal reads the same in every codebase. Only an infinite F (u == k exactly) is dropped: 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.

See also PFQN_GLD, PFQN_NCLD, PFQN_PANACEALD.

pfqn_rgf(L, N, Z)

Recursion by Generating Functions (RGF) for the normalizing constant of single-class closed product-form networks with replicated stations.

Exact normalizing constant of a single-class closed product-form network obtained by convolving the per-node generating-function sequences of Coury and Harrison (1997), Property 1, instead of the per-station Buzen recursion.

Each node contributes the coefficient sequence of its own generating function, and a GROUP of m stations sharing the same demand p is collapsed into the single negative-binomial sequence r(k) = C(k+m-1,k) p^k, i.e. the whole group costs one sequence rather than m convolution passes. The delay contributes the Poisson sequence r(k) = Z^k/k!. Convolving the G distinct sequences gives g(0..N) exactly, and lG = log g(N), X(n) = g(n-1)/g(n).

Cost O(G N^2) against Buzen’s O(M N), so RGF is the cheaper route precisely when the model is heavily replicated and the population is moderate (G N < M); it is otherwise kept for its exact group closed form. The whole recursion runs in the log domain, so no intermediate overflows or underflows are possible.

Reference: J. Coury, P. G. Harrison, “Asymptotic properties of queuing networks”, IEE Proc.-Comput. Digit. Tech. 144(5):247-254, 1997 (Property 1 and the five-sequence decomposition of Sec. 4).

Parameters:

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

Returns:

G – Normalizing constant. lG: Logarithm of the normalizing constant. lg: Logarithms of g(0), g(1), …, g(N) (1 x (N+1) vector).

pfqn_qdlin(L, N, Z, mu, nservers, tol, maxiter, wtol)

QD-LIN, the array-level twin of what SolverMVA computes for method=’qdlin’. The Linearizer of Chandy and Neuse (Commun. ACM 25(2), 1982) run inside the queue-dependent AMVA framework of Casale-Perez-Wang (IFIP PERFORMANCE 2015), so the load-dependent term g_k is evaluated at the CORRECTED arrival-instant queue rather than at the plain one.

THIS IS A TRANSCRIPTION OF solver_amvald.m TOGETHER WITH solver_amvald_forward.m, restricted to the domain a demand matrix describes: closed classes only, one chain per class, unit visits, PS queueing stations and one optional delay carrying Z. It is NOT an independent re-derivation, and it is not the Wang-Sevcik QDLIN of the same name in native Python and C++ before 2026-09-04, which was Bard-Schweitzer written out.

FOUR PROPERTIES OF THE REFERENCE ARE REPRODUCED DELIBERATELY:

THE GAMMA CORRECTION IS CLASS-AGGREGATE, IN SLICE 1. solver_amvald.m allocates the (K,M,K) per-class Linearizer array for qdlin but writes gamma(s,k) = sum_r Q_s(k,r)/(Nt-1) - sum_r Q(k,r)/Nt into it with two subscripts, which MATLAB linear-indexes to (s,k,1); slices 2..K stay zero while every reader indexes gamma per class. The correction that reaches the residence time is N_1*gamma(r,k,1) - [r==1]*gamma(r,k,1), which coincides with the queue-dependent AMVA form (Nt-1)*gamma_agg iff K == 1. method=’lin’ takes the per-class form instead. 2. A SINGLE-SERVER STATION STILL CARRIES A SOFTMIN TERM. The multiserver factor is pfqn_lldfun(1+arrival-instant total, [], nservers), whose softmin at c = 1 is not exactly 1, so qdlin does not reduce to a textbook single-server AMVA even when every station has one server. 3. THE WAIT FACTOR IS FLOORED AT wtol. This floor is LINE’s options.tol, a DIFFERENT knob from the convergence tolerance, and SolverMVA never sets it, so it stays at the lineDefaults 1e-4 while the fixed point converges to iter_tol 1e-6. MATLAB solver_amvald_forward.m does NOT carry it; native Python does, and removing it there was tried and reverted on 2026-09-04 because the unfloored python recursion diverges where MATLAB and C++ do not. See _kb/06-solver-catalog.md. 4. WHICH UTILIZATION IS REPORTED DEPENDS ON THE MODEL. The analyzer forwards the iterated Uchain to the deaggregation ONLY under lld, cd or jd scaling; with none of those the deaggregation recomputes T*S/c from the NOMINAL demand, and the two differ by the iteration residual.

MU AND NSERVERS ARE DIFFERENT MECHANISMS, unlike in pfqn_qdamva, which folds the multiserver curve into mu. Here mu is sn.lldscaling, an interpolated rate multiplier per station, and nservers is the server count feeding the softmin. A c-server station is nservers(k)=c, NOT a mu row of min(1:smax,c); the latter reproduces Queue.setLoadDependence, a different station.

Parameters:

L – (M x R) service demand matrix, queueing stations only.
N – (1 x R) population vector, finite.
Z – (1 x R) think time vector; a delay station carrying it is prepended to the station list when any entry is positive, exactly as the equivalent Network would hold one. Empty means no think time.
mu – (M x smax) load-dependent rate multipliers, sn.lldscaling; empty means none.
nservers – (M x 1) server counts; empty means one server everywhere.
tol – Convergence tolerance on the queue lengths (default 1e-6).
maxiter – Iteration budget (default 1000). The outer sweep and each inner sweep are capped at sqrt(maxiter) and the total number of forward evaluations at min(maxiter,10000), as in solver_amvald.
wtol – Floor on the wait factor (default 1e-4), LINE’s options.tol.

Returns:

Q – (M x R) mean queue lengths at the queueing stations. U: (M x R) per-class utilizations. R: (M x R) per-class residence times. X: (1 x R) per-class throughputs. C: (1 x R) per-class cycle times, think time included. iter: Number of forward evaluations performed.

[Q,U,R,X,C,ITER] = PFQN_QDLIN(L,N,Z,MU,NSERVERS,TOL,MAXITER,WTOL)

pfqn_usumbound(R, K, N)

Upper bound on the sum of device utilizations of a closed R-class network, from Dowdy et al. (1992), J. ACM 39(1), Expression (6).

Largest value the sum of device utilizations sum_k U_k,R can take in any closed product-form network with R classes, K devices and N customers (their Theorem 6): sum_k U_k,R <= (H-1) + (K-H+1)(N-H+1)/(K+N-2H+1), H = min(R,K). The bound is demand-free and nondecreasing in R, which is what makes it invertible into a lower bound on the number of necessary customer classes; see pfqn_minclasses. At R >= min(N,K) it reaches min(N,K), the trivial cap of one busy server per device.

The paper’s worked case is K = 2, N = 3, R = 1, giving 2N/(N+1) = 1.5: a measured sum of 1.6 then refutes the single-class assumption.

Parameters:

R – Number of single-customer classes (scalar, 1 <= R <= N).
K – Number of devices (scalar).
N – Total number of customers (scalar).

Returns:

Umax – Upper bound on sum_k U_k,R.

pfqn_ssd(L, N, Z, nservers)

Server-Station Disaggregation throughput bounds (Suri-Dallery 1986) for single-class closed networks with multiserver stations.

SSD multiserver bounds (SIGMETRICS 1986, Theorem 5). Each C_k-server station of loading L_k is bracketed by disaggregations: C_k balanced single-server stations of loading L_k/C_k (lower) and one single server of loading L_k/C_k (upper). With R_l=sum L_k, Y_l=max L_k/C_k, R_u=sum L_k/C_k, Y_u=R_u/K: X_l = N/(R_l+(N-1)Y_l) <= X(N) <= N/(R_u+(N-1)Y_u) = X_u, the upper bound taken jointly with the ABA bound min(N/R_l, C_b/L_b). O(K) cost, same order as BJB on single-server networks. With Z>0 the queueing terms carry the terminal-workload correction of Lazowska et al. 1984, Table 5.2: (N-1)Y_l/(1+Z/(N R_l)) on the lower bound and (N-1)Y_u/(1+Z/R_u) on the upper. Adding Z without it does not yield a bound.

Parameters:

L – Service demand vector (M x 1).
N – Population (scalar).
Z – Think time (scalar, default 0).
nservers – Per-station server counts C_k (M x 1, default all 1).

Returns:

Xlo – Lower throughput bound (Theorem 5). Xhi: Upper throughput bound (Theorem 5, joint with ABA).

pfqn_pbh(L, N, Z, level)

Performance Bound Hierarchy (Eager-Sevcik 1983) for single-class closed product-form networks.

Level-level Performance Bound Hierarchy throughput/queue bounds. i MVA steps from an ABA-initialized residence give nested optimistic/pessimistic bounds converging to exact MVA as level->N (Eager-Sevcik 1983, ACM TOCS 1(2):99-115, eqs. 5-13). Level 1 with Z=0 equals the BJB optimistic bound; level 2 with Z=0 equals the closed form 1 + S(N-1), S = sum L_k^2.

Parameters:

L – Service demand vector (M x 1).
N – Population (scalar).
Z – Think time (scalar, default 0).
level – Hierarchy level >= 0 (default 1); clamped to N.

Returns:

Xlo – Lower throughput bound (pessimistic residence). Xhi: Upper throughput bound (optimistic residence, joint with the

asymptotic bound and 1/max(L)).

Qlo: Per-station lower queue-length bound (M x 1), Little bracket. Qhi: Per-station upper queue-length bound (M x 1), Little bracket.

pfqn_ncjd(Z, N, mu, visits, options)

[G, LG, GTAB] = PFQN_NCJD(Z, N, MU, VISITS, OPTIONS)

Joint-dependent name of PFQN_NCOI: the balance-function convolution of a closed network whose station rates read the whole per-class occupancy vector.

The two names denote the SAME routine because the balanced-fairness recursion: Phi_i(0) = 1, 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. Any joint-dependent scaling eta_i(n) (sn.jdscaling) is therefore admissible here, with the usual proviso that the product form it induces is the balanced-fair one matched to that rate: Phi must stay positive for the result to be a distribution, which the routine enforces by zeroing the balance value of a non-positive rate.

Use PFQN_NCOI when the model is genuinely order independent and the name should say so; use PFQN_NCJD 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.

See also PFQN_NCOI, PFQN_MVAJD, PFQN_CLWJD, PFQN_PAS_NC.

pfqn_mcmc(L, N, Z, s, options)

Chen-O’Cinneide REGULARIZATION: throughputs and queue lengths of a closed multiclass product-form network, obtained by simulating a regularized network that has the same steady-state distribution.

Parameters:

L – Service demand matrix (stations x classes).
N – Closed population vector (1 x classes), finite and integer.
Z – Think time (1 x classes, or a matrix summed over its rows).
s – Server counts per station (stations x 1), Inf for an infinite-server station; [] means all stations single-server.
options – Solver options; .samples is the number of simulated service completions, .seed the RNG seed, .config.mcmc_batches the batch count and .config.mcmc_burnin the discarded warm-up fraction.

Returns:

X – Throughput estimate per class, X(r) = G(N-e_r)/G(N). Q: Mean queue length estimate per station and class. ci: Batch-means standard errors and two-sigma intervals for X and Q.

[X,Q,CI] = PFQN_MCMC(L,N,Z,S,OPTIONS)

Markov chain Monte Carlo estimator of the class throughputs X(r) = G(N-e_r)/G(N) and of the mean queue lengths Q(i,r) of a CLOSED multiclass product-form (BCMP, no type changes) network, by the REGULARIZATION algorithm of

W. Chen, C. A. O’Cinneide, “Towards a Polynomial-Time Randomized Algorithm for Closed Product-Form Networks”, ACM TOMACS 8(3):227-253, 1998.

The three steps of the paper are:

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

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 CI reports the paper’s two-sigma interval. 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 x R)
N - (1 x R)
Z - (1 x R) – rows; [] or zeros if the model has no delay.
s - (M x 1) – infinite server; [] (default) means all stations single-server.
options - solver options (optional) – .samples simulated service completions (1e5); .seed RNG seed for reproducibility (optional); .config.mcmc_batches batch count for the CIs (30); .config.mcmc_burnin discarded warm-up fraction (0.1).

Returns:

X - (1 x R) throughput estimates G(N-e_r)/G(N). Q - (M x R) mean queue lengths at the queueing stations; the delay

aggregate is not returned, the caller recovers it as Z.*X.

ci - struct with fields Xse, Xlo, Xhi (1 x R), Qse, Qlo, Qhi (M x R),: batches, samples and burnin. The intervals are two-sigma, as in the tables of the paper.

Example (the three single-server stations plus IS station of Example 5.1):: mu = [0.2 0.5 0.8]; sets = {[1 2 3],[1 2],[1 3],[2 3]}; L = zeros(3,4); Z = zeros(1,4); for c=1:4, L(sets{c},c) = 1./mu(sets{c})’; if c>1, Z(c) = 1/0.5; end, end X = pfqn_mcmc(L, 3*ones(1,4), Z, [], struct(‘samples’,1e5,’seed’,23000))

See also PFQN_NC, PFQN_MCI, PFQN_LS, PFQN_IS, PFQN_MVA.

pfqn_clust(L, N, Z, subnets, localclasses, inner, tol, maxiter)

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

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 (eq. 2.41) and the foreign classes into a per-centre utilization U_k (eq. 2.42),

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

Choosing the PE algorithm for every subnetwork reproduces global PE exactly, so the useful setting is Linearizer inside, PE outside: the cost then sits between pfqn_bs and pfqn_linearizer, which is the point of the method. The answer depends on the decomposition, which is an input.

When no decomposition is supplied the criterion of the paper is applied automatically: the cheap PAMB estimate (pfqn_pam) 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.

The name avoids pfqn_ca, which is the exact convolution algorithm.

Parameters:

L – Service demand matrix (stations x classes).
N – Population vector.
Z – Think time vector.
subnets – Cell array of station index vectors; their union must cover 1..M. Empty for the automatic decomposition above.
localclasses – Cell array of class index vectors, one per subnetwork, listing the local classes of that subnetwork. Empty for automatic.
inner – ‘lin’ (default) or ‘bs’, the algorithm run inside a subnetwork.
tol – Convergence tolerance (default: 1e-6).
maxiter – Maximum number of outer iterations (default: 1000).

Returns:

XN – System throughput. QN: Mean queue lengths. UN: Utilization. RN: Residence times. it: Number of outer iterations performed.

pfqn_bkue(L, N, Z)

Birman-Kogan uniform (van der Waerden) expansion for a single chain.

Birman and Kogan (Stochastic Models 8(3):543-563, 1992), Section 4. 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 van der Waerden uniform expansion keeps the pole and the saddle in one formula through the complementary error function, and so stays accurate on both sides of the crossing.

The published formula carries two typographical defects that the paper’s own Table 3 settles; see the notes in the code below.

Parameters:

L – Service demand vector (stations x 1), single class.
N – Population (scalar).
Z – Think time (default: 0).

Returns:

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

[G,LG] = PFQN_BKUE(L,N,Z)

Uniform expansion of the single chain normalizing constant.

pfqn_rd(L, N, Z, mu, options)

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

Parameters:

L – Service demand matrix.
N – Population vector.
Z – Think time vector.
mu – Load-dependent rate matrix.
options – Solver options.

Returns:

lGN – Logarithm of normalizing constant. Cgamma: Gamma correction factor.

pfqn_mvasjn(L, N, Z, scv, sjnset, V, options)

Mean value analysis of closed networks with shortest-job-next stations.

Approximate MVA for closed queueing networks in which a subset of the single-server stations schedules non-preemptively by shortest job next (SJN/SJF), the job size being known on arrival.

The SJN station is modelled by the conditional waiting time W(x,n) of a tagged customer whose service requirement is x, obtained from the arrival theorem as the sum of the residual life of the job in service, the work of the queued jobs that will be served before the tagged one, and the work of the jobs that overtake it while it waits (Kant 1992, eqs. 1-7):

W(x,n) = [ (1+CV^2) s U(n-1)/2 + X(n-1) phi(x,n-1) ] / [ 1 - X(n-1) theta(x) ], theta(x) = int_0^x t f(t) dt, phi(x,n) = int_0^x W(t,n) t f(t) dt, R(n) = s + int_0^inf W(x,n) f(x) dx.

The recursion is explicit: W(.,n) needs only phi(.,n-1), so it is carried alongside the population recursion of exact MVA. This is the unidirectional scheme of the reference and steps over the whole population lattice; pfqn_amvasjn is the fixed-point counterpart that trades the lattice for a Schweitzer closure on the same profile.

Two multiclass readings of “shortest job” are supported and selected by options.prio:

pooled (options.prio empty, the default): a job of any class with service requirement below x precedes the tagged job, so the sums over classes run over all of them. This is the direct multiclass reading of the derivation above and collapses to eq. (6) of the reference for a single class. - priority (options.prio a vector of distinct levels, 1 = highest): classes are non-preemptively prioritised at the SJN station and SJN applies only within a class, which is method A of Kant 1992, eq. (21). Method B of that paper, which evaluates the denominator at the non-integral population n - Q(n), is not implemented here; pfqn_nintmva provides the fractional-population recursion it needs.

The service time density at an SJN station is not an input: only its mean and squared coefficient of variation are, and the density is reconstructed by the two-moment branching-Erlang fit the reference prescribes (see sjn_fit). This makes theta(x) and the tail integrals closed form.

The x-integrals are evaluated on a uniform grid of options.ns subdivisions of [0, Lx] with Lx = options.Lfactor times the largest mean service time at the station, by composite Simpson, exactly as the reference does; W(.,n) is needed at the next population step so quadrature rules that sample at arbitrary abscissae cannot be used. Beyond Lx the profile is closed by the analytic tail W(x,n) = a_n - b_n exp(-c_n (x - Lx)) of eqs. (11)-(14).

SJN starves long jobs as the station saturates, and the arrival theorem then fails badly: when the fraction of work due to jobs no longer than x reaches one the recursion has no solution and the function throws ‘LINE:SjnStarvation’ rather than returning a value, which is the behaviour reported in the reference. Capping the utilization instead is not an option here: the cap rescales the conditional waiting time profile that the next population step reads back, so the correction compounds along the lattice and the recursion oscillates. pfqn_amvasjn solves the same profile by a fixed point and does cap, the iteration being self-consistent.

Reference: K. Kant, “MVA approximations for SJN scheduling”, Performance Evaluation 15(1):41-61, 1992.

Parameters:

L – Service demand matrix (M x R) of the queueing stations.
N – Population vector (1 x R).
Z – Think time vector (1 x R). Default: zeros.
scv – Squared coefficient of variation of the service times (M x R). Default: ones.
sjnset – Indices of the stations scheduling by SJN. Default: none.
V – Visit ratios (M x R), so that the per-visit service time is L./V. Default: ones.
options – Struct with fields ns (grid subdivisions, default 32), Lfactor (grid extent in mean service times, default 8) and prio (1 x R priority levels, default [] for the pooled reading).

Returns:

XN – System throughput (1 x R). QN: Mean queue length (M x R). UN: Utilization (M x R). CN: Residence time (M x R). WX: Struct array with the conditional waiting times at population N:

WX.station, WX.x (grid), WX.W (ns+1 x R) and WX.tail (R x 3 tail parameters a, b, c).

[XN,QN,UN,CN,WX] = PFQN_MVASJN(L,N,Z,SCV,SJNSET,V,OPTIONS)

pfqn_harel_ub(rho, N, n, Z)

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

Parameters:

rho – Relative utilizations (k x 1), all strictly positive.
N – Closed population, at least 1.
n – Extrapolation point, 2 <= n <= min(N,7).
Z – Think time; must be zero.

Returns:

UB – Throughput upper bound at population N.

UB = PFQN_HAREL_UB(RHO, N, N0, Z)

Upper bound alone of PFQN_HAREL_BOUNDS, 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)), A_i = sum_j rho_j^i,

from Harel, Namn and Sturm, “Simple bounds for closed queueing networks” (Queueing Systems 31, 1999). G is evaluated by the Newton-Girard recurrence n G(n) = sum_{i=1..n} A_i G(n-i); the n <= 7 ceiling is kept from the reference implementation. A nonzero think time is refused.

See also PFQN_HAREL_BOUNDS, PFQN_HAREL_LB.

pfqn_dnc(L, N)

Distinct-load Normalizing Constant (DNC) at a nonintegral population.

Normalizing constant and throughput of a single-class closed product-form network at a REAL-VALUED population, by partial-fraction inversion of the network generating function (Dowdy and Gordon 1984).

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 an analytic function of n. Evaluating it at a real n therefore interpolates the integral normalizing constants exactly (it reproduces them at every integer) and gives a smooth throughput curve X(N) = G(N-1)/G(N) through the integral points, rather than the rounding or linear interpolation the paper compares against. For all-distinct loads the coefficients have the closed form A_g = prod_{l~=g} x_g/(x_g - x_l), used directly; with repeated loads they are recovered from the M integral constants G(0..M-1), which determine them uniquely.

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 rejected here. Use pfqn_nintmva for nonintegral populations with a delay.

Reference: L. W. Dowdy, K. D. Gordon, “Algorithms for Nonintegral Degrees of Multiprogramming in Closed Queuing Networks”, Performance Evaluation 4(1):19-28, 1984.

Parameters:

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

Returns:

X – Throughput G(N-1)/G(N). G: Normalizing constant at population N. lG: Logarithm of the normalizing constant.

pfqn_busyp_clw(alpha, mu, P, N, subnet, n, gamma, isdelay, method)

Busy period of a subnetwork through normalizing-constant evaluations.

Parameters:

alpha – Relative arrival rates (JxR), one column per chain.
mu – Service rates (JxR), the chain-r rate at node j.
P – Routing matrix (JxJ), or a 1xR cell of per-chain matrices.
N – Population per chain (1xR), Inf entries for an open chain.
subnet – Indexes of the nodes forming the subnetwork.
n – Busy period order(s).
gamma – External arrival rates (JxR), empty for a closed network.
isdelay – Logical (Jx1), true at an infinite server.
method – Normalizing-constant method of the point evaluations (‘clw’).

Returns:

b – Mean busy period duration(s), same size as n.

B = PFQN_BUSYP_CLW(ALPHA, MU, P, N, SUBNET, N, GAMMA, ISDELAY, METHOD)

WHAT THIS BUYS OVER PFQN_BUSYP / PFQN_BUSYP_MULTICLASS. Those walk the whole population ladder (the whole lattice, multichain) because the numerator sums over {|m| >= n}. The complement of that set is the SHELLS |m| <= n-1, and summing the product form over the WHOLE lattice is the full network’s own normalizing constant, G_I and H convolving to it:

sum_{m : |m| >= n} G_I(m) H(N-m) = G(N) - sum_{m : |m| <= n-1} G_I(m) H(N-m)

so the busy period of order n needs only the n lowest shells plus ONE evaluation of G(N). The ordinary busy period n=1 collapses to three constants,

b(1,I) = [G(N) - H(N)] / sum_r A_r(I) H(N-e_r)

all of them point evaluations at or near the full population, which is what the normalizing-constant methods are built for. This routine calls CLW (Choudhury-Leung-Whitt, J. ACM 42, 1995, numerical inversion of the generating function); any method returning lG(N) can take its place. The cost stops depending on N: O(shells up to n-1) plus O(n*R) constant evaluations, against O(lattice) for the ladder routines.

THE OPEN CASE NEEDS NO INVERSION. The subnetwork’s constant sequence has generating function g(z) = prod_{i in I} f_i(z), and the tail is g(1) minus a partial sum,

b(n,I) = [g_I(1) - sum_{m=0}^{n-1} G_I(m)] / [G_I(n-1) C(I)],

with f_i(1) = 1/(1-rho_i) at a single server and exp(rho_i) at an infinite one. That removes the tail TRUNCATION of the ladder routine, not just its cost: the tail is exact.

ACCURACY. The numerator is a difference of two nearly equal quantities when the level set is unlikely, so the relative error grows with n: measured 6e-12 at n=1 against 2.7e-08 at n=N on a three-station closed model at N=20. Cost grows with n too, so the routine is most accurate where it is fastest.

SCOPE. CLW’s generating function covers single-server and infinite-server stations, so a general load-dependent scaling belongs to PFQN_BUSYP_MULTICLASS. The identity is for the AGGREGATE level set: a per-class one has complement {m_r <= n-1}, the whole lattice in the other chains, which buys nothing.

pfqn_xzabalow(L, N, Z)

Lower ABA (asymptotic bound analysis) bound on throughput. Not the Zahorjan-Balanced bound (that is pfqn_xzgsblow); single-class only.

Parameters:

L – Service demand vector.
N – Population.
Z – Think time.

Returns:

XN – Lower bound on throughput.

pfqn_jointmarg(n, L, N, infset, lGN, engine)

Joint probability of the per-station TOTAL queue lengths.

Parameters:

n – Per-station total queue lengths.
L – Service demand matrix, infinite servers included as rows.
N – Population vector.
infset – Rows of L that are infinite-server stations (optional).
lGN – Log normalizing constant at N (optional).
engine – Permanent engine (optional).

Returns:

pjoint – Joint probability of the total queue lengths. lpjoint: Its logarithm.

[PJOINT, LPJOINT] = PFQN_JOINTMARG(N, L, NPOP, INFSET, LGN, ENGINE)

Joint probability that station I holds N(I) jobs IN TOTAL, all classes summed out, in a closed multiclass product-form network:

P(n_1,…,n_M) = perm(A) / ( prod_r N_r! * prod_{j in INFSET} n_j! * G(N) )

with A the demand matrix whose column R is repeated N_r times and whose row I is repeated n_i times, so A is square of order sum(N). Unlike PFQN_JOINT, which takes the delay as a single aggregated row, every infinite-server station keeps its own row here and contributes its own 1/n_j! – the queueing stations contribute the n_i! that the permanent identity supplies, the infinite servers do not.

Input:

N - (M x 1) per-station total queue lengths, infinite servers: included; sum(N) must equal sum(NPOP)

L - (M x R) demand matrix, infinite-server rows included NPOP - (1 x R) per-class populations INFSET - indices of the rows of L that are infinite-server stations,

empty by default (every station is a queue)

LGN - log normalizing constant; computed with PFQN_CA when omitted,: aggregating the infinite-server rows into the think time (which is exact: the delay stations aggregate by the multinomial theorem, so G does not depend on how they are split)
ENGINE - ‘exact’ (default), ‘spm’, ‘bethe’, ‘heur’, ‘huberlaw’ or: ‘adapart’

Output:

PJOINT - joint probability LPJOINT - its logarithm, which survives populations the probability

itself underflows at

ENGINE ‘spm’ IS THE ONLY ONE THAT DOES NOT EXPAND THE MATRIX. The other approximate engines are handed A itself, of order SUM(N), where every column has multiplicity one and their accuracy decays with the population. PERM_SPM takes the ROW-replicated matrix with the class populations as its column multiplicities, which is the regime its expansion is asymptotically exact in: the Laplace integral has dimension R-1 whatever the population is, so the cost is independent of SUM(N) and the relative error is O((R-1)/min(N)). Measured on a 3-station 2-class model: 12.8% at N = (1,1), 4.2% at (3,3), 2.1% at (6,6), tracking 1/(8 min N_r). It degrades the other way round, when the CLASS COUNT grows at fixed population (2.7% at R = 2, 21% at R = 7, all at N_r = 3), because R-1 is the dimension being expanded in.

The bias is nearly constant across the lattice, so a caller sweeping the whole state space and RENORMALIZING to sum to one keeps far less of it: the total variation distance against the exact law is 5.0e-3 at N = (1,1), 8.4e-4 at (3,3) and 4.3e-4 at (5,5), better than ‘bethe’ and ‘heur’ at every population measured. Cost against ‘exact’, which collapses the repeated rows and columns and so grows with M and R rather than with SUM(N): comparable at R = 2, 12x faster at R = 4 and 16x at R = 7.

The identity holds for load-independent single-server queues plus infinite servers. Multiserver and load-dependent stations break the n_i! factor and are the caller’s responsibility to exclude.

ZERO ELEMENTS are safe under the exact engine and only under it: a station holding no jobs contributes no row, a class with no jobs contributes no column, a zero demand is an ordinary zero entry of A, and the permanent of the empty matrix is 1. The approximate engines are REFUSED on a matrix with a structural zero rather than having it floored at eps: Sinkhorn scaling needs full support, and the Bethe gap is a state-dependent lower bound that does not cancel when the estimates are normalized against each other.

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

pfqn_sens_mvaldmx_validate()

Validation harness for pfqn_sens_mvaldmx, the mixed load-dependent moment analysis of Akyildiz and Strelen (1991).

Five independent references, chosen so that every channel of the derivation is exercised by something that does not share its code:

pfqn_mvaldmx for the base measures X, Q, U, R. The primal must be reproduced entry by entry, otherwise the derivative is of the wrong function. B. central finite differences of pfqn_mvaldmx with respect to the demand-scaling parameter y(j,s). This checks the differentiated recursion itself, including the load-dependent channel dEC/dLo and the open-class channel dLo/dy of eq. (21), but does not check the identity Cov = d nbar / dy. C. pfqn_sens_mva in the closed load-independent limit. This checks the identity against the independently validated de Souza e Silva and Muntz recursion. D. brute-force enumeration of the product-form equilibrium distribution. Closed load-dependent models are enumerated exactly; mixed models are enumerated with the open populations truncated, which converges geometrically and is therefore checked at a looser tolerance. This is the only check that closes the loop on the identity in the mixed load-dependent case. E. symmetry of QCovFull. Cov[n(i,r),n(j,s)] and Cov[n(j,s),n(i,r)] are computed by differentiating two different classes’ equations, so their agreement is a nontrivial structural check.

pfqn_nrl(L, N, Z, alpha, options)

Normalizing constant via Norlund-Rice Logit (NRL) approximation.

Parameters:

L – Service demand matrix.
N – Population vector.
Z – Think time vector.
alpha – Load-dependent rate matrix.
options – Solver options.

Returns:

lG – Logarithm of normalizing constant.

pfqn_sens(L, N, Z, mi)

Exact analytic performance sensitivities for closed product-form queueing networks. Dispatches to a CoMoM-backed kernel for the single-station repairman model and to differentiated MVA otherwise.

Exact derivatives of the mean performance measures {X,Q,U,R} of a closed product-form (BCMP) queueing network with respect to the service demands L(i,r) and the think times Z(r). The derivatives are analytic (exact to machine precision), not finite differences.

Two exact kernels are dispatched transparently (local subfunctions): - sens_comom : Class-Oriented Method of Moments. Selected for the repairman model (a single single-server queue, M=1, plus a think-time delay with sum(Z)>0 and strictly positive demands). Polynomial in the number of classes R via normalizing-constant moment relations (queue-length covariances / replica moments). - sens_mva : forward-mode differentiation of the exact Reiser-Lavenberg MVA recursion. Used for every other model. Both kernels return the identical struct layout, so the choice is an internal optimization invisible to callers.

Reference: Z. Liu and P. Nain, “Sensitivity Results in Open, Closed and Mixed Product-Form Queueing Networks”, INRIA RR-1144, 1989; X.-R. Cao and D.-J. Ma, “Performance sensitivity formulae, algorithms and estimates for closed queueing networks with exponential servers”, Performance Evaluation 26:181-199, 1996; G. Casale, “CoMoM: Efficient Class-Oriented Evaluation of Multiclass Performance Models”, IEEE TSE 2011.

Parameters:

L – Service demand matrix (M x R), L(i,r) = visits_ir / rate_ir.
N – Population vector (1 x R).
Z – Think time vector (1 x R). Default: zeros.
mi – (Optional) Station residence multiplicity (1 x M). Default: ones.

Returns:

sens –

A struct with base measures and their Jacobians:

.X (1 x R), .Q (M x R), .U (M x R), .R (M x R) base MVA measures

(X system throughput per class, Q mean queue length, U utilization,
R residence time per visit-chain, i.e. CN of pfqn_mva).

.params 1 x P struct array describing each differentiation parameter,
fields .type (‘L’ or ‘Z’), .station (i, 0 for Z), .class (r).

.dX (R x P), .dQ (M x R x P), .dU (M x R x P), .dR (M x R x P)
derivatives of each base measure w.r.t. parameter p. For a ‘L’ parameter at (i,r) the derivative is d(.)/dL(i,r); for a ‘Z’ parameter at class r it is d(.)/dZ(r).

.QCov (M x R x M x R) QCov(i,r,j,s) = Cov[n(i,r),n(j,s)], the exact
queue-length covariance, a by-product of the Jacobian.

.QVar (M x R) QVar(i,r) = Var[n(i,r)]. .QTotVar (M x 1) QTotVar(i) = Var[sum_r n(i,r)]. .QCovAsym (scalar) roundoff-level residual of the moment recursion,

see pfqn_sens_mva.

Notes:

Mirrors pfqn_mva(L,N,Z,mi) exactly for the base measures (single-server or residence-multiplicity mi stations plus an infinite-server delay Z).
Derivatives w.r.t. a service rate mu(i,r) follow by the chain rule d(.)/dmu(i,r) = -(L(i,r)/mu(i,r)) * d(.)/dL(i,r).

pfqn_nc_sanitize(lambda, L, N, Z, atol)

Sanitize and preprocess network parameters for NC solvers.

Parameters:

lambda – Arrival rate vector.
L – Service demand matrix.
N – Population vector.
Z – Think time vector.
atol – Absolute tolerance.

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.

erase empty classes

pfqn_bjbk(L, N, Z, k)

Iterative BJB(k) balanced job bounds (Zahorjan et al.) for single-class closed networks with delay.

BJB(k): the iterative Balanced Job Bound at iteration count k (Casale-Muntz-Serazzi 2008, Tables 5,7). BJB(1) recovers the noniterative Balanced Job Bound; each additional iteration performs one exact MVA step from the balanced seed at N-k, and the bracket tightens to exact as k->N. Realized through the validated PBH recursion (pfqn_pbh), whose level-1 optimistic bound equals the noniterative BJB optimistic bound (Eager-Sevcik 1983, eq. 10).

Parameters:

L – Service demand vector (M x 1).
N – Population (scalar).
Z – Think time (scalar, default 0).
k – Iteration count >= 1 (default 1).

Returns:

Xlo – Lower throughput bound. Xhi: Upper throughput bound.

pfqn_mvamx(lambda, D, N, Z, mi)

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

Parameters:

lambda – Arrival rate vector.
D – Service demand matrix.
N – Population vector.
Z – Think time vector.
mi – Queue replication factors (default: ones).

Returns:

XN – System throughput. QN: Mean queue lengths. UN: Utilization. CN: Cycle times. lGN: Logarithm of normalizing constant.

[XN,QN,UN,CN,LGN] = PFQN_MVAMX(LAMBDA,D,N,Z, MI)

pfqn_mvams(lambda, L, N, Z, mi, S)

General-purpose MVA for mixed networks with multiserver nodes.

Parameters:

lambda – Arrival rate vector.
L – Service demand matrix.
N – Population vector.
Z – Think time vector.
mi – Queue replication factors (default: ones).
S – Number of servers per station (default: ones).

Returns:

XN – System throughput. QN: Mean queue lengths. UN: Utilization. CN: Residence times (M x R), as in PFQN_MVA. lG: Logarithm of normalizing constant.

[XN,QN,UN,CN,LOGG]=PFQN_MVAMS(LAMBDA,L,N,Z,MI,S)

Standard arrival theorem throughout. For the interlocked-flow correction of Franks (1999), Ch. 4, Eq. (4.7), call PFQN_MVAMS_ILOCK instead.

pfqn_comom(L, N, Z, atol)

CoMoM algorithm for computing the normalizing constant.

Parameters:

L – Service demand matrix.
N – Population vector.
Z – Think time vector.
atol – Tolerance.

Returns:

lG – Logarithm of the normalizing constant.

pfqn_clwoi(Z, N, mu, visits, options)

[G, LG] = PFQN_CLW_OI(Z, N, MU, VISITS, OPTIONS)

Normalizing constant of a closed product-form network made of an aggregated infinite-server (delay) node and an arbitrary number of ORDER-INDEPENDENT (OI) / pass-and-swap stations with empty swap graph, obtained by numerically inverting the multichain generating function with the lattice-Poisson algorithm of Choudhury, Leung and Whitt (J. ACM 42(5):935-970, 1995).

This is the OI counterpart of PFQN_CLW_LLD and the transform counterpart of the convolution routine PFQN_NCOI. Both return the same G(N); they differ in cost, see below.

GENERATING FUNCTION. The multichain generating function factorizes over the stations,

G(z) = exp( sum_r Z_r z_r ) prod_i F_i(z), F_i(z) = sum_n Phi_i(n) z^n,

with Phi_i the v-weighted balanced-fairness balance function of OI station i,

Phi_i(0) = 1, mu_i(n) Phi_i(n) = sum_{r: n_r>0} v_{i,r} Phi_i(n - e_r).

Unlike a load-dependent station, whose factor collapses to a function of the single argument sum_r rho_{ri} z_r, an OI station factor genuinely depends on the whole vector z, because mu_i(n) depends on the per-class occupancy n through its SUPPORT supp(n) = {r : n_r > 0}. It is nevertheless RATIONAL and available in closed form. Splitting the count lattice by support, on which mu_i(n) = mu_{i,S} is constant, and writing F_{i,S}(z) for the part of F_i carried by the states with supp(n) = S, the balance recursion gives

( 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 minus r}(z), F_{i,{}} = 1,

F_i(z) = sum_{S subseteq {1..R}} F_{i,S}(z), (*)

since removing one class-r job from a state of support S lands on support S when n_r >= 2 and on S minus r when n_r = 1. The singularities are the |S| hyperplanes sum_{r in S} v_{i,r} z_r = mu_{i,S}, one per support, in place of the single pole x = c_i of the load-dependent case. For R = 1 and constant rate c, (*) returns c/(c - v z), and a load-independent single-server queue (mu_{i,S} = 1 for all S) returns 1/(1 - sum_r v_{i,r} z_r), i.e. exactly the factor PFQN_CLW inverts.

INVERSION. G(N) is the coefficient of prod_r z_r^{N_r}, recovered by R nested one-dimensional lattice-Poisson inversions (CLW eq. 2.3) on contours of radius r_j = 10^{-gamma_j/(2 l_j N_j)}. The restrictive static scaling of CLW eqs. 5.41-5.46 is applied to the EXPANDED constraint matrix that lists one row per (station, nonempty support) pair, with unit-pole intensities v_{i,r}/mu_{i,S}: each such row is one singular hyperplane of (*), so the scaling keeps the whole contour inside the domain of analyticity exactly as the single-pole normalization does for PFQN_CLW_LLD. Recovery is in the log domain (CLW eq. 7.1).

SCOPE. The rates must be support-only, mu_i(n) = mu_i(supp(n)), which is the defining property of an OI station and what makes (*) a finite rational function. Every rate handle is therefore verified EXHAUSTIVELY on the count lattice 0 < n <= N before the inversion, at prod_r (N_r+1) evaluations per station (below the prod_r 2 l_r N_r contour points spent afterwards), and a state whose rate differs from that of its support is an error naming that state; a rate that varies inside a support is a general balanced-fairness station and must go to PFQN_NCOI. The violation is refused rather than warned-and-inverted because the inversion would otherwise return a plausible but wrong G(N) with no other symptom. A non-empty swap graph breaks the closure of Phi on the count vector altogether and requires the microstate routine PFQN_PAS_NC. Load-independent single-server queues need no special casing: pass them as OI stations with mu_i(n) = 1 and v_{i,r} = D_{i,r}.

COST. prod_r 2 l_r N_r contour points, each costing O(M R 2^R) for the M station transforms, against O(M prod_r (N_r+1)(N_r+2)/2) for the convolution of PFQN_NCOI. The inversion is therefore LINEAR rather than quadratic in each population and wins on large populations with few chains, while the 2^R per-point factor makes it lose as the number of chains grows. Unlike PFQN_NCOI it returns G at the single population N: throughputs need the R additional inversions at N - e_r.

Parameters:

Z - (1 x R) – Z(r) = 1/sigma_r for a delay with per-class rate sigma_r.
N - (1 x R)
mu - cell array {1 x M} of function handles, one per OI station. Each – mu{m}(n) returns the total service rate of station m at the per-class occupancy vector n (1 x R) and must depend on n only through its support. May be empty to model a pure delay network.
visits - (M x R) matrix, or {1 x M} cell of (1 x R) – class visit ratios v_{i,r} weighting the balance recursion. Default: unit visits.
options - struct with optional fields – .l (1 x R) inner lattice parameters l_j (roundoff control). .gamma (1 x R) aliasing parameters gamma_j (aliasing ~ 10^-gamma_j). Defaults follow CLW: l_1=1,g_1=11; l_2=l_3=2,g=13; l_j>=4=3,g=15.

Returns:

G - Normalizing constant G(N). Inf if it overflows the double range. lG - log(G(N)) (always finite when G > 0).

Example (delay + one OI station whose rate is the number of busy servers):: murate = @(n) sum(n > 0); G = pfqn_clwoi([1 2], [8 6], {murate});

See also PFQN_CLW_LLD, PFQN_CLW, PFQN_NCOI, PFQN_MVAOI, PFQN_PAS_NC.

pfqn_ab_amva(D, N, V, nservers, sched, fcfsSchmidt, marginalProbMethod)

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

Parameters:

D – Service time matrix (M x K).
N – Population vector (1 x K).
V – Visit ratio matrix (M x K).
nservers – Number of servers at each station (M x 1).
sched – Scheduling strategies for each station.
fcfsSchmidt – Whether to use Schmidt formula for FCFS stations.
marginalProbMethod – Method for marginal probability calculation (“ab” or “scat”).

Returns:

QN – Mean queue lengths. UN: Utilization. RN: Residence times. CN: Cycle times. XN: System throughput. totiter: Total number of iterations.

[QN,UN,RN,CN,XN,totiter] = PFQN_AB_AMVA(D,N,V,NSERVERS,SCHED,FCFSSCHMIDT,MARGINALPROBMETHOD)

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

pfqn_oi_is(N, mu, options)

[G, LG, Q] = PFQN_OI_IS(N, MU, OPTIONS)

Importance-sampling (IS) estimate of the normalizing constant of a closed two-station order-independent (OI) tandem. This is PFQN_PAS_IS with an EMPTY swap graph: it does not restrict the sampled orderings to a pass-and-swap communicating class but samples over ALL microstates (every ordering of the present jobs is feasible), thereby recovering the plain OI normalizing constant of PFQN_NCOI (which it estimates by Monte Carlo rather than exact balanced-fairness convolution).

Model. Two OI stations (1 = upstream, 2 = downstream) hold all N jobs (no delay). The ordered-state chain is irreducible, so the communicating class is the full set of orderings D = all permutations of the job multiset. Writing Phi_m for the balanced-fairness balance function of station m,

G = sum_{c in D} sum_{k=0}^{ell} Phi_1(c_{1..k}) Phi_2(c_{ell..k+1}),

with the ordered product Phi_m(q) = prod_{p=1}^{|q|} 1/mu_m(n(q_{1..p})), n(.) the per-class COUNT vector of the prefix (Comte-Dorsman product form: pi(c) = (1/G) prod_j 1/mu(c_1..c_j) prod_s sigma_s^{N_s-n_s}/(N_s-n_s)!). The argument is the prefix MULTISET, not its support: OI property P1 asks only that mu be permutation-invariant, i.e. a function of the counts. The two coincide for a compatibility rate that reads only which classes are present, and differ for any count-dependent OI rate – an INF station (mu(n) = sum_r n_r sigma_r) above all, which is exactly what SOLVER_NC_PAS_IS_ANALYZER feeds in as station 2 of a Delay + OI cycle.

Auto-normalized IS. Orderings c are drawn by placing, at each step, a uniformly random present class (no placement constraint); the draw probability p(c) is the product of the reciprocal branching factors. Then

G[xi] = E_{C~p}[ (sum_k xi(C,k) Phi_1(C_{1..k}) Phi_2(C_{ell..k+1})) / p(C) ],

and E[xi] = G[xi]/G[1] reuses the same samples (auto-normalized IS). Taking xi = number of class-r jobs in the prefix yields the class-r mean queue length at station 1.

Parameters:

N - (1 x R) closed population vector (the macrostate)
mu - cell {1 x 2} of function handles. mu{m} (n) – rank rate of station m for the per-class occupancy (count) vector n (1 x R). Permutation-invariant (OI property P1), but NOT necessarily a function of supp(n) alone. This is the svcRateFun stored on an OI station.
options - solver options (optional) – .samples number of IS samples (default 1e4); .seed RNG seed for reproducibility (optional); .verbose print progress (default false); .qlen estimate the queue lengths too (default true). False

estimates ONLY G: the prefix-count coefficients are neither allocated nor accumulated and Q comes back zero. The ordering is drawn from the same stream either way, so G is unchanged to the last bit.

Returns:

G - IS estimate of the OI normalizing constant (== PFQN_NCOI with Z=0). lG - log(G). Q - (2 x R) IS estimate of the mean per-class queue length; Q(1,:) at

station 1, Q(2,:) = N - Q(1,:) at station 2.

Example (two OI queues, R=5, no swap graph):: mu1 = [1 2 0.5]; nb1 = {1,2,3,[1 3],[2 3]}; mu2 = [1 2]; nb2 = {1,2,[1 2],1,2}; rate = @(nb,muv,n) sum(muv(unique([nb{n>0}]))); mu = {@(n) rate(nb1,mu1,n), @(n) rate(nb2,mu2,n)}; [G,lG,Q] = pfqn_oi_is([1 1 1 3 3], mu, struct(‘samples’,1e5));

See also PFQN_PAS_IS, PFQN_NCOI, PFQN_OI_FNC.

pfqn_mvaldmx(lambda, D, N, Z, mu, S)

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

Parameters:

lambda – Arrival rate vector.
D – Service demand matrix.
N – Population vector.
Z – Think time vector.
mu – Load-dependent rate matrix.
S – Number of servers per station.

Returns:

XN – System throughput. QN: Mean queue lengths. UN: Utilization. CN: Cycle times. lGN: Logarithm of normalizing constant. Pc: Marginal queue-length probabilities.

[XN,QN,UN,CN,lGN,Pc] = PFQN_MVALDMX(LAMBDA,D,N,Z,MU,S)

pfqn_mvaldms(lambda, D, N, Z, S)

Load-dependent MVA for multiserver mixed networks (wrapper for pfqn_mvaldmx).

Parameters:

lambda – Arrival rate vector.
D – Service demand matrix.
N – Population vector.
Z – Think time vector.
S – Number of servers per station.

Returns:

XN – System throughput. QN: Mean queue lengths. UN: Utilization (adjusted for multiservers). CN: Cycle times. lGN: Logarithm of normalizing constant.

[XN,QN,UN,CN] = PFQN_MVALDMS(LAMBDA,D,N,Z,S) Wrapper for pfqn_mvaldmx that adjusts utilizations to account for multiservers

pfqn_mvald(L, N, Z, mu, stabilize)

Exact MVA for load-dependent closed queueing networks.

Parameters:

L – Service demand matrix.
N – Population vector.
Z – Think time vector.
mu – Load-dependent rate matrix (MxNt).
stabilize – Force non-negative probabilities (default: true).

Returns:

XN – System throughput. QN: Mean queue lengths. UN: Utilization. CN: Cycle times. lGN: Logarithm of normalizing constant evolution. isNumStable: Numerical stability flag. pi: Marginal queue-length probabilities.

[XN,QN,UN,CN,LGN]=PFQN_MVALD(L,N,Z,MU)

pfqn_mva(L, N, Z, mi)

Exact Mean Value Analysis (MVA) for product-form queueing networks.

Parameters:

L – Service demand matrix (M x R).
N – Population vector (1 x R).
Z – Think time vector (1 x R).
mi – (Optional) Additive term of the residence-time recursion C(i,s)=L(i,s)*(mi(i)+Qarv), 1 for a queueing station. Default: ones. THIS IS NOT A SERVER COUNT: mi(i)=c INFLATES the residence time by c rather than adding c servers. For multiserver stations call PFQN_MVAMS(lambda,L,N,Z,mi,S), which passes S to the load-dependent recursion with mu(i,n)=min(n,S(i)).

Returns:

XN – System throughput (1 x R). QN: Mean queue length (M x R). UN: Utilization (M x R). CN: Residence time (M x R). lGN: Logarithm of the normalizing constant.

[XN,QN,UN,CN,LGN] = PFQN_MVA(L,N,Z,MI) [XN,QN,UN,CN] = pfqn_mva(L,N,Z,mi)

Standard arrival theorem. For the interlocked-flow correction of Franks (1999), Ch. 4, Eq. (4.7), call PFQN_MVA_ILOCK instead.

pfqn_mushift(mu, iset)

Shift load-dependent service rate vector by removing first element.

Parameters:

mu – Load-dependent rate matrix (MxN).
iset – Set of station indices to shift.

Returns:

mushifted – Shifted rate matrix (Mx(N-1)).

shifts the service rate vector

pfqn_mu_ms(N, m, c)

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

Parameters:

N – Maximum population.
m – Number of identical stations.
c – Number of servers per station.

Returns:

mu – Load-dependent service rate vector (1xN).

calculates the load-dependent rate of m identical c-server FCFS stations

pfqn_mmsample2(L, N, Z, samples)

Monte Carlo sampling for repairman models using McKenna-Mitra form.

Parameters:

L – Service demand vector.
N – Population vector.
Z – Think time vector.
samples – Number of samples.

Returns:

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

[G,LG] = PFQN_MMSAMPLE2(L,N,Z,SAMPLES)

pfqn_mmint2_gausslegendre(L, N, Z, m)

McKenna-Mitra integral with Gauss-Legendre quadrature.

Parameters:

L – Service demand vector.
N – Population vector.
Z – Think time vector.
m – Replication factor (default: 1).

Returns:

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

[G,LOGG] = PFQN_MMINT2_GAUSSLEGENDRE(L,N,Z,m)

Integrate McKenna-Mitra integral form with Gauss-Legendre in [0,1e6]

pfqn_mmint2_gausslaguerre(L, N, Z, m)

McKenna-Mitra integral with Gauss-Laguerre quadrature.

Parameters:

L – Service demand vector.
N – Population vector.
Z – Think time vector.
m – Replication factor (default: 1).

Returns:

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

[G,LOGG] = PFQN_MMINT2_GAUSSLAGUERRE(L,N,Z,m)

Integrate with Gauss-Laguerre

pfqn_mmint2(L, N, Z)

McKenna-Mitra integral form for repairman models using MATLAB integral.

Parameters:

L – Service demand vector.
N – Population vector.
Z – Think time vector.

Returns:

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

[G,LOGG] = PFQN_MMINT2(L,N,Z)

pfqn_mci(D, N, Z, I, variant)

Monte Carlo Integration (MCI) for normalizing constant.

Parameters:

D – Service demand matrix.
N – Population vector.
Z – Think time vector.
I – Number of samples (default: 1e5).
variant – MCI variant (‘mci’, ‘imci’, ‘amci’, ‘lhsmci’, ‘rm’; 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). Note that 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 – Normalizing constant estimate. lG: Logarithm of normalizing constant. lZ: Individual random sample log values.

[G,LG,LZ] = PFQN_MCI(D,N,Z,I,VARIANT)

Normalizing constant estimation via Monte Carlo Integration

Syntax: [G,lG,lZ] = pfqn_mci(D,N,Z,I,VARIANT) Input: D - demands (queues x classes) N - populations (1 x classes) Z - think times (1 x classes) I - samples VARIANT - ‘mci’, ‘imci’, ‘amci’, ‘lhsmci’, ‘rm’

Output: lG - estimate of logG lZ - individual random samples

Note: if the script returns a floating point range exception, double(log(mean(exp(sym(lZ))))) provides a better estimate of lG, but it is very time consuming due to the symbolic operations.

Implementation: Giuliano Casale (g.casale@imperial.ac.uk), 16-Aug-2013

pfqn_ls(L, N, Z, I)

Logistic sampling approximation for normalizing constant.

Parameters:

L – Service demand matrix (MxR).
N – Population vector (1xR).
Z – Think time vector (1xR).
I – Number of samples (default: 1e5).

Returns:

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

[GN,LGN]=PFQN_LS(L,N,Z,I)

pfqn_lldfun(n, lldscaling, nservers)

AMVA-QD load and queue-dependent scaling function.

Parameters:

n – Queue population vector.
lldscaling – Load-dependent scaling matrix.
nservers – Number of servers per station.

Returns:

r – Scaling factor vector.

R = PFQN_LLDFUN(N,MU,C)

pfqn_ld_is(L, N, Z, mu, options)

[G, LG] = PFQN_LD_IS(L, N, Z, MU, OPTIONS)

Importance-sampling (IS) estimate of the normalizing constant of a closed LOAD-DEPENDENT product-form queueing network. This is the 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 and each contributes 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, writing ell = sum(N) and letting a “cut vector” split 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 (n_1,…,n_S) 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=0}^{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 and is estimated by the sample mean.

Parameters:

L - (M x R)
N - (1 x R)
Z - (1 x R) aggregated think time (delay)
mu - load-dependent capacities. Either an (M x ell) – mu(i,k) the capacity of station i holding k jobs, or a cell {1 x M} of function handles mu{i}(k), or [] for the load-independent case mu(i,k) = 1 (see PFQN_IS).
options - solver options (optional) – .samples number of IS samples (default 1e4); .seed RNG seed for reproducibility (optional).

Returns:

G - IS estimate of the normalizing constant G(N). lG - log(G).

Example (2 queues + delay, load-dependent):: L = [0.5 0.3; 0.2 0.4]; N = [3 2]; Z = [1 1]; mu = [1 2 2 2 2; 1 1 1 1 1]; % station 1 is a 2-server queue [G,lG] = pfqn_ld_is(L, N, Z, mu, struct(‘samples’,1e5));

See also PFQN_IS, PFQN_OI_IS, PFQN_PAS_IS, PFQN_NCLD, PFQN_NC.

pfqn_linearizermx(lambda, L, N, Z, nservers, type, tol, maxiter, method, QN0)

Linearizer for mixed open/closed queueing networks.

Parameters:

lambda – Arrival rate vector (inf for closed classes).
L – Service demand matrix.
N – Population vector (inf for open classes).
Z – Think time vector.
nservers – Number of servers per station.
type – Scheduling strategy per station.
tol – Convergence tolerance (default: 1e-8).
maxiter – Maximum iterations (default: 1000).
method – Linearizer variant (‘lin’, ‘gflin’, ‘egflin’, default: ‘egflin’).
QN0 – (M x R) queue lengths that warm-start the Bard-Schweitzer initialization; empty for the default cold start.

Returns:

QN – Mean queue lengths. UN: Utilization. WN: Waiting times. CN: Cycle times. XN: System throughput. totiter: Total iterations.

function [Q,U,W,C,X,totiter] = PFQN_LINEARIZERMX(lambda,L,N,Z,nservers,type,tol,maxiter,method)

pfqn_linearizerms(L, N, Z, nservers, type, tol, maxiter, QN0)

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

Parameters:

L – Service demand matrix.
N – Population vector.
Z – Think time vector.
nservers – Number of servers per station.
type – Scheduling strategy per station (default: PS).
tol – Convergence tolerance (default: 1e-8).
maxiter – Maximum number of iterations (default: 1000).
QN0 – (M x R) queue lengths that warm-start the Bard-Schweitzer initialization; empty for the default cold start.

Returns:

Q – Mean queue lengths. U: Utilization. R: Residence times. C: Cycle times. X: System throughput. totiter: Total iterations performed.

Multiserver version of Krzesinski’s Linearizer as described in Conway 1989, Fast Approximate Solution of Queueing Networks with Multi-Server Chain- Dependent FCFS Queues. Some minor adjustments based on De Souza-Muntz’s description of the algorithm.

pfqn_linearizer(L, N, Z, type, tol, maxiter, QN0)

Linearizer approximation for single-server stations.

Parameters:

L – Service demand matrix.
N – Population vector.
Z – Think time vector.
type – Scheduling strategy type per station.
tol – Convergence tolerance; ‘cn’ or NaN selects the Chandy-Neuse (1982) population-scaled termination test, see pfqn_cntol.
maxiter – Maximum number of iterations.
QN0 – (M x R) queue lengths that warm-start the Bard-Schweitzer initialization; empty for the default cold start.

Returns:

Q – Mean queue lengths. U: Utilization. W: Waiting times. C: Cycle times. X: System throughput. totiter: Total iterations performed.

Single-server version of linearizer

pfqn_le(L, N, Z)

Logistic expansion (LE) asymptotic approximation for normalizing constant.

Parameters:

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

Returns:

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

[GN,LGN]=PFQN_LE(L,N,Z)

pfqn_lcfsqn_nc(alpha, beta, N)

Computes the normalizing constant for LCFS queueing networks

Parameters:

alpha – Service rates at LCFS station (1xR vector).
beta – Service rates at LCFS-PR station (1xR vector).
N – Population vector (default: ones(1,R)).

Returns:

G – Normalizing constant. Ax: Cell array of A matrices for each state.

[G,AX] = PFQN_LCFSQN_NC(ALPHA, BETA, N) Normalizing constant for multiclass LCFS queueing networks

This function computes the normalizing constant for a 2-station closed queueing network with:

Station 1: LCFS (Last-Come-First-Served, non-preemptive)

Station 2: LCFS-PR (LCFS with Preemption-Resume)

Parameters:

alpha - vector of inverse service rates at station 1 (LCFS) – alpha(r) = 1/mu(1,r) for class r
beta - vector of inverse service rates at station 2 (LCFS-PR) – beta(r) = 1/mu(2,r) for class r
N - population vector, N (r)

Returns:

G - normalizing constant Ax - cell array of A matrices for each state x=0:K

Reference:: G. Casale, “A family of multiclass LCFS queueing networks with order-dependent product-form solutions”, QUESTA 2026.

pfqn_lcfsqn_mva(alpha, beta, N)

Exact MVA for 2-station LCFS queueing networks.

Parameters:

alpha – Service rates at LCFS station (1xR vector).
beta – Service rates at LCFS-PR station (1xR vector).
N – Population vector (default: ones(1,R)).

Returns:

T – Throughput vector. Q: Mean queue lengths (2xR matrix). U: Utilization (2xR matrix). B: Back probability matrix (2xR).

[T,Q,U,B] = PFQN_LCFSQN_MVA(ALPHA, BETA, N) Mean Value Analysis for multiclass LCFS queueing networks

This function computes performance metrics for a 2-station closed queueing network with:

Station 1: LCFS (Last-Come-First-Served, non-preemptive)

Station 2: LCFS-PR (LCFS with Preemption-Resume)

Parameters:

alpha - vector of inverse service rates at station 1 (LCFS) – alpha(r) = 1/mu(1,r) for class r
beta - vector of inverse service rates at station 2 (LCFS-PR) – beta(r) = 1/mu(2,r) for class r
N - population vector, N (r) – (default: ones(1,R) - one job per class)

Returns:

T - throughput vector, T(r) = throughput of class r Q - queue length matrix, Q(i,r) = mean queue length at station i, class r U - utilization matrix, U(i,r) = utilization at station i, class r B - back probability matrix, B(i,r) = probability class r job is at

back of queue at station i

Note: This implementation uses log-space arithmetic to prevent numerical underflow. The results are mathematically exact (up to floating-point precision) - no approximations are made.

Reference:: G. Casale, “A family of multiclass LCFS queueing networks with order-dependent product-form solutions”, QUESTA 2026.

pfqn_lcfsqn_ca(alpha, beta, N)

Convolution algorithm for 2-station LCFS queueing networks

Parameters:

alpha – Service rates at LCFS station (1xR vector).
beta – Service rates at LCFS-PR station (1xR vector).
N – Population vector (default: ones(1,R)).

Returns:

G – Normalizing constant. V: Auxiliary normalization term.

[G,V] = PFQN_LCFSQN_CA(ALPHA, BETA, N) Convolution algorithm for multiclass LCFS queueing networks

This function computes the normalizing constant for a 2-station closed queueing network with:

Station 1: LCFS (Last-Come-First-Served, non-preemptive)

Station 2: LCFS-PR (LCFS with Preemption-Resume)

Parameters:

alpha - vector of inverse service rates at station 1 (LCFS) – alpha(r) = 1/mu(1,r) for class r
beta - vector of inverse service rates at station 2 (LCFS-PR) – beta(r) = 1/mu(2,r) for class r
N - population vector, N (r) – (default: ones(1,R) - one job per class)

Returns:

G - normalizing constant V - auxiliary normalization term

Reference:: G. Casale, “A family of multiclass LCFS queueing networks with order-dependent product-form solutions”, QUESTA 2026.

pfqn_mwrbb(V, S, N, Z, sched, prio)

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, as defined by Majumdar and 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 independent of the scheduling discipline. 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 the per-visit queueing delay bound d_kc^+ depends on the discipline at station k – FIFO (Theorem 1 / eq. 6, via Lemma 1), processor sharing (Lemma 2), preemptive priority (Lemma 3) and non-preemptive priority (Lemmas 4-5). The coupled inequalities are resolved by the interval- narrowing fixed point that reproduces the BNR-Prolog robust box bounds; for a single FIFO class it reduces to the Muntz-Wong asymptotic bounds.

Bounds are insensitive to the service-time distributions (only NBUE is assumed) and to routing dependencies; only the mean visits V, mean service demands S, populations N, think times Z, per-station disciplines and per- class priorities are required. The think time Z aggregates the pure-delay (infinite-server) stations; only queueing stations are passed in V,S.

Parameters:

V – (K x C) mean visits of class c at queueing station k
S – (K x C) mean service demand per visit of class c at station k
N – (1 x C) population of class c
Z – (1 x C) think time (pure delay) of class c
sched – (K x 1) discipline code per station: 0=FIFO (default),
prio – (1 x C) class priority, lower value = higher priority

Returns:

Xlo – (1 x C) lower bound on class throughput (Theorem 2) Xup: (1 x C) upper bound on class throughput (eqs. 2-3) Wlo: (K x C) per-visit residence time at station k for class c

Examples

[Xlo,Xup,Wlo] = pfqn_mwrbb(V,S,N,Z,sched,prio)

pfqn_kt(L, N, Z)

Knessl-Tier asymptotic expansion for normalizing constant.

Parameters:

L – Service demand matrix.
N – Population vector.
Z – Think time vector (default: zeros).

Returns:

G – Normalizing constant. lG: Logarithm of normalizing constant. X: System throughput. Q: Mean queue lengths.

Knessl-Tier asymptotic expansion, fixed to include the IS/think-time term (defect G18) and to evaluate the exponent at the exact saddle point.

Derivation. In LINE’s convention the generating function of G over the population vector is

sum_N G(N) prod_r u_r^N_r = exp(sum_r Z_r u_r) prod_k (1-sum_r L_kr u_r)^-1

so Cauchy extraction and steepest descent on: F(u) = sum_r Z_r u_r - sum_k log(1-U_k) - sum_r N_r log u_r, U_k = L(k,:)*u
give: G ~ (2 pi)^{-R/2} det(H)^{-1/2} exp(F(u*)) / prod_r u*_r H_rs = delta_rs N_r/u_r^2 + sum_k L_kr L_ks/(1-U_k)^2

where u* solves N_r = u_r Z_r + sum_k u_r L_kr/(1-U_k) (asymptotic MVA fixed point). This is Knessl-Tier Result 2 rewritten in LINE’s convention: their rho_k absorbs the think rate and the factor exp(sum_r Z_r u_r) is contained in Psi(1,y*) via the M!/(M-n)! Stirling terms. Stock pfqn_kt dropped the linear think term (it is absent from both the exponent and the Hessian) and evaluated the exponent at the AQL throughput.

pfqn_joint(n, L, N, Z, lGN)

Compute joint queue-length probability distribution.

Parameters:

n – Queue-length state vector (total or per-class).
L – Service demand matrix.
N – Population vector.
Z – Think time vector (optional).
lGN – Log normalizing constant at N (optional, computed if not provided).

Returns:

pjoint – Joint probability of state n.

pjoint = pfqn_joint(n,L,N,Z,lGN)

Compute the joint queue-length probability for vector n=(n_1,…,n_M) or (n_{11},…n_{M,R}), with M the number of queues, and R the number of classes. If there is a think time Z, then n is just the queue population. Optionally, the logarithm of the normalizing constant G at N can be passed as an input to reduce the computational cost.

Examples: * Total queue-lengths

pjoint=[]; for n=0:4 % queue 1

for z=0:4-n % think time
pjoint(end+1) = pfqn_joint([n;4-n-z],[10,2;5,4],[2,2],[91,92]);

end

end sum(pjoint)

Per-class queue-lengths pjoint=[]; for n1=0:4 % queue 1, class 1

for n2=0:3 % queue 1, class 2

for z1=0:4-n1 % think time, class 1

for z2=0:3-n2 % think time, class 2
pjoint(end+1) = pfqn_joint([n1,n2;4-n1-z1,3-n2-z2],[10,2;5,4],[4,3],[91,92]);

end

end

end

end sum(pjoint)

pfqn_grnmol(L, N)

Normalizing constant using Grundmann-Moeller quadrature.

Parameters:

L – Service demand matrix.
N – Population vector.

Returns:

G – Normalizing constant.

pfqn_is(L, N, Z, options)

[G, LG] = PFQN_IS(L, N, Z, OPTIONS)

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

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

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

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

which is unbiased for the exact constant of PFQN_NC.

Parameters:

L - (M x R)
N - (1 x R)
Z - (1 x R) aggregated think time (delay)
options - solver options (optional) – .samples number of IS samples (default 1e4); .seed RNG seed for reproducibility (optional).

Returns:

G - IS estimate of the normalizing constant G(N). lG - log(G).

Example (2 queues + delay):: L = [0.5 0.3; 0.2 0.4]; N = [3 2]; Z = [1 1]; [G,lG] = pfqn_is(L, N, Z, struct(‘samples’,1e5));

See also PFQN_LD_IS, PFQN_OI_IS, PFQN_PAS_IS, PFQN_NC.

pfqn_gld(L, N, mu, options)

Exact normalizing constant for load-dependent queueing networks.

Parameters:

L – Service demand matrix (MxR).
N – Population vector (1xR).
mu – Load-dependent rate matrix (Mx sum(N)).
options – Solver options.

Returns:

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

[G,LG]=PFQN_GLD(L,N,MU,OPTIONS)

pfqn_gflinearizer(L, N, Z, type, tol, maxiter, alpha, QN0)

Generalized fixed-point Linearizer with uniform scaling exponent.

Parameters:

L – Service demand matrix.
N – Population vector.
Z – Think time vector.
type – Scheduling strategy type per station.
tol – Convergence tolerance; ‘cn’ or NaN selects the Chandy-Neuse (1982) population-scaled termination test, see pfqn_cntol.
maxiter – Maximum number of iterations.
alpha – Uniform scaling exponent for all classes.
QN0 – (M x R) queue lengths that warm-start the Bard-Schweitzer initialization; empty for the default cold start.

Returns:

Q – Mean queue lengths. U: Utilization. W: Waiting times. C: Cycle times. X: System throughput. totiter: Total iterations performed.

pfqn_fnc(alpha, c)

Generate load-dependent rates for functional server model f(n)=n+c.

Parameters:

alpha – Rate parameters (Mx N matrix).
c – Constant offset parameter (default: auto-determined).

Returns:

mu – Load-dependent service rates. c: Determined offset constant.

generate rates for functional server f(n)=n+c

pfqn_expand(QN, UN, CN, mapping, M_original)

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

Parameters:

QN – Queue lengths from reduced model (M’ x R).
UN – Utilizations from reduced model (M’ x R).
CN – Cycle times from reduced model (M’ x R).
mapping – Mapping vector from pfqn_unique (1 x M), mapping(i) = unique station index.
M_original – Original number of stations M.

Returns:

QN_full – Queue lengths in original dimensions (M x R). UN_full: Utilizations in original dimensions (M x R). CN_full: Cycle times in original dimensions (M x R).

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

[QN_FULL, UN_FULL, CN_FULL] = PFQN_EXPAND(QN, UN, CN, MAPPING, M_ORIGINAL)

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

Input:: QN - M’ x R queue lengths from reduced model UN - M’ x R utilizations from reduced model CN - M’ x R cycle times from reduced model mapping - 1 x M vector from pfqn_unique (mapping(i) = unique station index) M_original - original number of stations M
Output:: QN_full - M x R queue lengths in original dimensions UN_full - M x R utilizations in original dimensions CN_full - M x R cycle times in original dimensions

pfqn_egflinearizer(L, N, Z, type, tol, maxiter, alpha, QN0, npasses)

Extended generalized fixed-point Linearizer approximation.

Parameters:

L – Service demand matrix.
N – Population vector.
Z – Think time vector.
type – Scheduling strategy type per station.
tol – Convergence tolerance (default: 1e-8); ‘cn’ or NaN selects the Chandy-Neuse (1982) population-scaled termination test, see pfqn_cntol.
maxiter – Maximum number of iterations (default: 1000).
alpha – Per-class scaling exponent vector.
QN0 – (M x R) queue lengths that warm-start the Bard-Schweitzer initialization; empty for the default cold start.
npasses – Number of Delta refresh passes (default 3, the Chandy-Neuse rule; pfqn_scat sets 1).

Returns:

Q – Mean queue lengths. U: Utilization. W: Waiting times. C: Cycle times. X: System throughput. totiter: Total iterations performed.

Single-server version of linearizer

pfqn_cub(L, N, Z, order, atol)

Cubature method for normalizing constant using Grundmann-Moeller rules.

Parameters:

L – Service demand matrix (MxR).
N – Population vector (1xR).
Z – Think time vector (1xR).
order – Degree of cubature rule (default: ceil((sum(N)-1)/2)).
atol – Absolute tolerance (default: 1e-8).

Returns:

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

[GN,LGN]=PFQN_CUB(L,N,Z,ORDER,TOL)

pfqn_conwayms(L, N, Z, nservers, type, tol, maxiter, QN0)

Multiserver Linearizer approximation (Conway 1989).

Parameters:

L – Service demand matrix.
N – Population vector.
Z – Think time vector.
nservers – Number of servers per station.
type – Scheduling strategy type per station (default: FCFS).
tol – Convergence tolerance (default: 1e-8).
maxiter – Maximum number of iterations (default: 1000).
QN0 – (M x R) queue lengths that warm-start the Bard-Schweitzer initialization; empty for the default cold start.

Returns:

Q – Mean queue lengths. U: Utilization. R: Residence times. C: Cycle times. X: System throughput. totiter: Total number of iterations.

Multiserver version of Linearizer as described in Conway 1989, Fast Approximate Solution of Queueing Networks with Multi-Server Chain- Dependent FCFS Queues

pfqn_ncld(L, N, Z, mu, varargin)

Normalizing constant for load-dependent closed networks.

Parameters:

L – Service demand matrix.
N – Population vector.
Z – Think time vector.
mu – Load-dependent rate matrix.
varargin – Optional solver parameters.

Returns:

lG – Logarithm of normalizing constant. G: Normalizing constant. method: Method used for computation.

[LGN,G,METHOD] = PFQN_NCLD(L,N,Z,VARARGIN)

pfqn_comomrm_orig(L, N, Z, atol)

Original CoMoM implementation for finite repairman model.

Parameters:

L – Service demand matrix.
N – Population vector.
Z – Think time vector.
atol – Absolute tolerance for numerical computations.

Returns:

lG – Logarithm of normalizing constant.

comom for a finite repairment model

pfqn_comomrm_ms(L, N, Z, m, S)

CoMoM for multiserver repairman model.

Parameters:

L – Service demand matrix.
N – Population vector.
Z – Think time vector.
m – Replication factor (default: 1).
S – Number of servers at queueing stations.

Returns:

G – Normalizing constant. lG: Logarithm of normalizing constant. prob: State probability distribution.

m: replication factor S: number of servers at the queueing stations

pfqn_comomrm_ld(L, N, Z, mu, options)

CoMoM for repairman model with load-dependent service rates.

Parameters:

L – Service demand matrix.
N – Population vector.
Z – Think time vector.
mu – Load-dependent rate matrix (MxNt matrix).
options – Solver options.

Returns:

G – Normalizing constant. lG: Logarithm of normalizing constant. prob: State probability distribution.

S: number of servers at the queueing stations

pfqn_comomrm(L, N, Z, m, atol)

CoMoM (Class-Oriented Method of Moments) for finite repairman model.

Parameters:

L – Service demand matrix.
N – Population vector.
Z – Think time vector.
m – Replication factor (default: 1).
atol – Absolute tolerance for numerical computations.

Returns:

lG – Logarithm of normalizing constant. lGbasis: Logarithm of basis functions.

comom for a finite repairment model

pfqn_lap(L, N, Z)

Laplace approximation for normalizing constant.

Parameters:

L – Service demand matrix.
N – Population vector.
Z – Think time vector.

Returns:

logI – Logarithm of normalizing constant approximation.

[LOGI] = PFQN_LAP(L,N,Z)

pfqn_nc(lambda, L, N, Z, varargin)

pfqn_ca(L, N, Z)

Convolution Algorithm for exact normalizing constant computation.

Parameters:

L – Service demand matrix.
N – Population vector.
Z – Think time vector.

Returns:

Gn – Normalizing constant. lGn: Logarithm of the normalizing constant.

pfqn_clwjd(Z, N, mu, visits, lcut, options)

[G, LG] = PFQN_CLWJD(Z, N, MU, VISITS, LCUT, OPTIONS)

Normalizing constant of a closed product-form network made of an aggregated infinite-server (delay) node and an arbitrary number of LIMITED JOINT-DEPENDENT (LJD) stations, obtained by numerically inverting the multichain generating function with the lattice-Poisson algorithm of Choudhury, Leung and Whitt (J. ACM 42(5):935-970, 1995).

This is the joint-dependent generalization of PFQN_CLWOI, which is the special case LCUT = 1. It stands to PFQN_CLWOI as PFQN_CLW_LLD stands to PFQN_CLW: a per-station cutoff beyond which the rate stops changing turns an infinite series into a rational function of the same denominators.

LIMITED JOINT DEPENDENCE. Station i has a rate mu_i(n) that reads the whole per-class occupancy vector n (a joint-dependent scaling, sn.jdscaling), but saturates coordinatewise: with a cutoff vector l_i = LCUT(i,:),

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

i.e. past l_{i,r} further class-r jobs no longer change the rate. The clipped vector t ranges over the finite box prod_r {0,…,l_{i,r}}. Order independence is l_i = 1 (t is then the support indicator); a multiserver station with c servers is l_i = c, since min(sum n, c) is a function of the clipped vector once every l_{i,r} >= c; a load-independent queue is l_i = 1 with a constant rate.

GENERATING FUNCTION. The multichain generating function factorizes over the stations,

G(z) = exp( sum_r Z_r z_r ) prod_i F_i(z), F_i(z) = sum_n Phi_i(n) z^n,

with Phi_i the v-weighted balance function of station i,

Phi_i(0) = 1, mu_i(n) Phi_i(n) = sum_{r: n_r>0} v_{i,r} Phi_i(n - e_r).

Splitting the count lattice by clipped region, on which mu_i is constant, and writing F_{i,t} for the part of F_i carried by the states with t_i(n) = t,

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

F_i(z) = sum_{t in box} F_{i,t}(z). (*)

The two sides differ because removing a class-r job leaves the region only where the coordinate is UNsaturated: for t_r < l_{i,r} the region fixes n_r = t_r, so n - e_r lands in t - e_r; for t_r = l_{i,r} the region is n_r >= l_{i,r}, so n - e_r lands in t (n_r > l) or in t - e_r (n_r = l), which is what puts v_{i,r} z_r on the left. Hence the singularities are the hyperplanes sum_{r in S} v_{i,r} z_r = mu_{i,t} over the SATURATED sets S = {r : t_r = l_{i,r}}: at most 2^R of them per station, however large the cutoffs are. Setting l_i = 1 reduces (*) to the support recursion of PFQN_CLWOI, and R = 1 reduces it to Bertozzi-McKenna eq. 2.19.

INVERSION. G(N) is the coefficient of prod_r z_r^{N_r}, recovered by R nested one-dimensional lattice-Poisson inversions (CLW eq. 2.3) on contours of radius r_j = 10^{-gamma_j/(2 l_j N_j)}. The restrictive static scaling of CLW eqs. 5.41-5.46 runs on the expanded constraint matrix listing one row per (station, saturated set S), carrying the unit-pole intensities v_{i,r}/min{mu_{i,t} : saturated set of t is S}: the smallest rate over the regions sharing a saturated set is the binding one. Rows dominated by a superset of no larger rate are dropped first. Recovery is in the log domain (CLW eq. 7.1).

SCOPE. The rate must be constant on each clipped region; this is checked on probe states and is an error otherwise. Any rate is admissible with LCUT = N (the default), the clipping being vacuous on the reachable lattice, at the cost given below. Rates that never saturate have no finite rational transform and are exactly the case where a large LCUT is mandatory.

COST. prod_r 2 l_r N_r contour points, each costing O(M R prod_r (LCUT(i,r)+1)) for the M station transforms, against O(M prod_r (N_r+1)(N_r+2)/2) for the convolution of PFQN_NCJD. The inversion is linear rather than quadratic in each population, but the per-point region box grows with the cutoff: with LCUT = N the box is the whole lattice and the convolution wins outright. The inversion pays off exactly when the joint dependence saturates early.

Parameters:

Z - (1 x R)
N - (1 x R)
mu - cell array {1 x M} of function handles, one per LJD station. Each – mu{m}(n) returns the total service rate of station m at the per-class occupancy vector n (1 x R). May be empty for a pure delay network.
visits - (M x R) matrix, or {1 x M} cell of (1 x R) – class visit ratios v_{i,r}. Default: unit visits.
lcut - (M x R) – l_{i,r} >= 1, or a scalar/row broadcast to every station. Entries are clipped to N_r, which is exact: a rate difference at n_r > N_r can only move coefficients with n_r > N_r, and every term reaching z^N has n_r <= N_r at every station. Default: N (no truncation).
options - struct with optional fields – .l (1 x R) inner lattice parameters l_j (roundoff control). .gamma (1 x R) aliasing parameters gamma_j (aliasing ~ 10^-gamma_j). Defaults follow CLW: l_1=1,g_1=11; l_2=l_3=2,g=13; l_j>=4=3,g=15.

Returns:

G - Normalizing constant G(N). Inf if it overflows the double range. lG - log(G(N)) (always finite when G > 0).

Example (delay + a station whose rate saturates at two jobs per class):: murate = @(n) 1 + sum(min(n, 2)); G = pfqn_clwjd([1 2], [8 6], {murate}, [], [2 2]);

See also PFQN_CLWOI, PFQN_CLW_LLD, PFQN_NCJD, PFQN_MVAJD, PFQN_NCOI.

pfqn_amvasjn(L, N, Z, scv, sjnset, V, options)

Approximate MVA of closed networks with shortest-job-next stations.

Fixed-point (Bard-Schweitzer) counterpart of pfqn_mvasjn for closed networks with non-preemptive shortest-job-next (SJN/SJF) stations.

pfqn_mvasjn carries the conditional waiting time profile W(x,n) over the whole population lattice, which costs prod(N+1) steps and rules the method out for large populations. The closure used here rests on the observation that

lam_k(n) W_k(x,n) f_k(x) dx

is the mean number of queued class-k customers whose service requirement lies in (x, x+dx), that is, the queue length resolved by job size. Schweitzer’s assumption is applied to that density rather than to its integral: removing one customer of class r scales the class-r size-resolved queue length by (N_r-1)/N_r and leaves the other classes unchanged. Integrating over x recovers the usual Schweitzer rule for the aggregate queue lengths, so the closure is the exact analogue of the one applied at the ordinary stations, and the two are used together consistently.

The unknowns are therefore the profiles W_r(x,N) on the quadrature grid together with the queue lengths, and they are found by successive substitution. Cost per iteration is O(M R ns), against the prod(N+1) M R ns of the exact recursion, and the population may be arbitrarily large. What is given up is the population dependence of the shape of W(x): the closure lets its level scale but keeps its shape fixed, whereas the true profile stiffens with the load because the denominator 1 - sum_k lam_k theta_k(x) sharpens. The error therefore concentrates at high utilization, where the SJN approximation is already at its weakest, and pfqn_mvasjn should be preferred whenever the lattice is affordable.

The response time equation, the two-moment branching-Erlang fit of the size distribution, the quadrature grid, the analytic tail and the pooled and priority multiclass readings are all shared with pfqn_mvasjn; see that function and the reference for their derivation.

Reference: K. Kant, “MVA approximations for SJN scheduling”, Performance Evaluation 15(1):41-61, 1992. The bidirectional use of the priority equations, of which this is the limiting form, is discussed in section 3.1 of that paper.

Parameters:

L – Service demand matrix (M x R) of the queueing stations.
N – Population vector (1 x R).
Z – Think time vector (1 x R). Default: zeros.
scv – Squared coefficient of variation of the service times (M x R). Default: ones.
sjnset – Indices of the stations scheduling by SJN. Default: none.
V – Visit ratios (M x R), so that the per-visit service time is L./V. Default: ones.
options – Struct with fields ns (grid subdivisions, default 32), Lfactor (grid extent in mean service times, default 8), prio (1 x R priority levels, default [] for the pooled reading), tol (default 1e-8) and iter_max (default 1000).

Returns:

XN – System throughput (1 x R). QN: Mean queue length (M x R). UN: Utilization (M x R). CN: Residence time (M x R). WX: Struct array with the converged conditional waiting times:

WX.station, WX.x (grid), WX.W (ns+1 x R) and WX.tail (R x 3 tail parameters a, b, c).

it: Number of iterations performed.

[XN,QN,UN,CN,WX,IT] = PFQN_AMVASJN(L,N,Z,SCV,SJNSET,V,OPTIONS)

mexify_pfqn

UNTITLED Generate static library pfqn_bs

See also CODER, CODER.CONFIG, CODER.TYPEOF, CODEGEN.

pfqn_busyp(alpha, mu, P, N, subnet, n, gamma, tol)

Mean busy period of order n for a subnetwork of a product-form network.

Parameters:

alpha – Relative arrival rates (1xJ), solution of the traffic equations.
mu – Load-dependent service rates (JxK), mu(j,k) with k jobs at node j, or a handle mu(j,kvec) when the rates do not saturate.
P – Routing matrix (JxJ).
N – Population (scalar, Inf for an open network).
subnet – Indexes of the nodes forming the subnetwork.
n – Busy period order(s), 1 <= n <= N.
gamma – External arrival rates (1xJ), empty for a closed network.
tol – Relative tolerance of the open-network tail truncation.

Returns:

b – Mean busy period duration(s), same size as n. lG: Log normalizing constants of the subnetwork, orders 0..K. lH: Log normalizing constants of the complement, orders 0..N.

[B,LG,LH] = PFQN_BUSYP(ALPHA, MU, P, N, SUBNET, N, GAMMA, TOL)

Mean duration of the busy period of order n for the subnetwork SUBNET, that 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.

Implements H. Daduna, “Busy Periods for Subnetworks in Stochastic Networks: Mean Value Analysis”, J. ACM 35(3), 1988, Theorem 1 (closed Gordon-Newell network) and Theorem 3 (open Jackson network). Both are evaluated in the log domain, which serves the same purpose as the ratio recursions of Corollaries 2 and 4, namely avoiding the overflow of the individual normalizing constants.

The paper is single-chain: ALPHA is the stochastic solution of x*P = x for a closed network and the solution of x = GAMMA + x*P for an open one, and every node is a state-dependent single-server FCFS station. By the insensitivity noted in the paper (Section 5) the result depends on the service processes only through the rates MU.

pfqn_bs(L, N, Z, tol, maxiter, QN0, type)

Bard-Schweitzer Approximate Mean Value Analysis (MVA).

Parameters:

L – Service demand matrix.
N – Population vector.
Z – Think time vector.
tol – Tolerance for convergence; ‘cn’ or NaN selects the Chandy-Neuse (1982) population-scaled termination test, see pfqn_cntol.
maxiter – Maximum number of iterations.
QN0 – Initial guess for queue lengths.
type – Scheduling strategy type (default: PS).

Returns:

XN – System throughput. QN: Mean queue lengths. UN: Utilization. RN: Residence times. it: Number of iterations performed.

[XN,QN,UN,RN]=PFQN_BS(L,N,Z,TOL,MAXITER,QN)

pfqn_aql(L, N, Z, TOL, MAXITER, QN0)

Approximate Queue Length (AQL) algorithm for product-form networks.

Parameters:

L – Service demand matrix.
N – Population vector.
Z – Think time vector.
TOL – Tolerance for convergence.
MAXITER – Maximum number of iterations.
QN0 – Initial guess for queue lengths.

Returns:

XN – System throughput. QN: Mean queue lengths. UN: Utilization. RN: Residence times. numIters: Number of iterations. AN: Average arrival rate at nodes.

pfqn_clw_lld(L, N, Z, mu, options)

Computes the normalization constant g(K) of a multichain closed product-form network with limited load-dependent (LLD) stations and (optionally) infinite-server delay by numerically inverting its p-dimensional generating function.

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

F_i(x) = sum_{n>=0} x^n / prod_{k=1}^n S_i(k),

and S_i(k) = mu(i,k) is the load-dependent rate scaling with k jobs at queue i. 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 (c_i = number of servers) and load-independent (F_i = 1/(1-x)) queues are special cases. Since g(K) depends on S_i(k) only for k <= sum(K), any general load-dependent input is truncated to LLD at l_i <= sum(K) without loss of exactness.

g(K) is the coefficient of prod_j z_j^{K_j}, recovered by p nested one-dimensional lattice-Poisson inversions (CLW, JACM 42(5):935-970, 1995, eq. 2.3) with a restrictive static scaling adapted from CLW eqs. 5.41-5.46: each queue is normalized by its pole c_i (unit-pole form, simple pole) and log-domain recovery (eq. 7.1) is applied.

Parameters:

L – (q’ x p) single-server relative traffic intensities, L(i,j)=rho_{ji}.
N – (1 x p) closed-chain population vector K.
Z – (1 x p) aggregate infinite-server relative intensities rho_{j0} (think-time term). Default: zeros(1,p).
mu – (q’ x n) load-dependent rate scalings, mu(i,k) = S_i(k); if fewer than sum(N) columns are given the last column is extended (LLD assumption). Default: ones (all queues load-independent).
options – struct with optional fields: .l (1 x p) inner lattice parameters l_j (roundoff control). .gamma (1 x p) aliasing parameters gamma_j (aliasing ~ 10^-gamma_j). Defaults follow CLW: l_1=1,g_1=11; l_2=l_3=2,g=13; l_j>=4=3,g=15.

Returns:

G – Normalization constant g(K). Inf if it overflows double range. lG: Natural logarithm of g(K) (always finite).

Scope: cost is prod_j 2 l_j K_j contour points, each requiring O(sum_i l_i) work, so the routine is practical for moderate populations and few chains. The numerator polynomials are evaluated in double precision; extreme LLD cutoffs (l_i > ~170 with large c_i) may overflow.

pfqn_cftp(L, N, S, nsamples, method)

Draws states exactly distributed according to the product-form stationary distribution of a closed single-class Jackson network with multiple servers, using monotone Coupling From The Past (Propp-Wilson) as proposed by Kijima and Matsui (WSC 2005, “Approximate/Perfect Samplers for Closed Jackson Networks”).

Parameters:

L – Service demand (visit ratio / service rate) vector, one entry per station, L(i) = theta_i/mu_i.
N – Total closed population K (scalar, single class).
S – Number of servers per station (default: 1; use Inf for a delay/infinite-server station).
nsamples – Number of independent samples to draw (default: 1).
method – ‘cftp’ for exact/perfect sampling (default) or ‘approx’ for the rapidly-mixing approximate sampler M_A.

Returns:

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

(‘approx’), returned as (nsamples x 1) for diagnostics.

[Q,X,T] = PFQN_CFTP(L,N,S,NSAMPLES,METHOD)

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

Reference: S. Kijima and T. Matsui, “Approximate/Perfect Samplers for Closed Jackson Networks”, Proc. Winter Simulation Conference, 2005.

Input: L - demands (stations x 1), L(i) = theta_i/mu_i N - total population K (scalar) S - servers per station (stations x 1), Inf for infinite server NSAMPLES - number of independent samples (default 1) METHOD - ‘cftp’ (exact, default) or ‘approx’ (rapidly mixing M_A)

Output: Q - empirical mean queue length (1 x stations) X - sampled states (nsamples x stations), each row sums to N T - coalescence horizon / mixing steps per sample (nsamples x 1)

ljd_linearize(nvec, cutoffs)

IDX = LJD_LINEARIZE(NVEC, CUTOFFS)

Convert per-class population vector to linearized index

nvec: [n1, n2, …, nK] - per-class populations cutoffs: [N1, N2, …, NK] - per-class cutoffs

Returns: 1-based linearized index

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

pfqn_mvaoi(Z, N, mu, Dli, visits, options)

[X, QOI, QLI, QDELAY, SOI] = PFQN_MVAOI(Z, N, MU, DLI, OPTIONS)

Mean-value analysis of a closed product-form queueing network composed 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) / pass-and-swap stations with empty swap graph. This is the mean-value counterpart of PFQN_NCOI and the marginal-distribution form PFQN_MVAOI_MARG: it returns the same exact per-class throughput and queue-lengths but WITHOUT computing any normalizing constant or joint marginal, using only mean quantities (throughputs, demands, queue-lengths) evaluated on shifted models. 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), whose “third form” (rate depending on the full per-class occupancy vector) is realized here, extended to MULTIPLE OI stations by carrying one rate-shift vector s_i per OI station i.

Throughout, r and s index job classes; i indexes OI stations; j indexes LI queues. For a single OI station and no LI queue the analysis recurs on the shift vector s_i (the OI occupancy already committed at the bottom of station i), Nn = N - s_i the jobs still to distribute:

Q^{(S)}(Nn) = sum_r U_r^{(S)}(Nn) ( e_r + Q^{(S+e_r@i)}(Nn - e_r) ), U_r^{(S)}(Nn) = D_r^{(S)}(Nn) X_r^{(S)}(Nn) (bottom-job utilization),

with the class-r OI demand and throughput satisfying: D_r^{(S)}(Nn) = (1/mu_i(s_i+e_r)) rho_{i,r}^{(S)}(Nn-e_r), Nn_r = 1, D_r^{(S)}(Nn) = [X_r^{(S)}(Nn-e_r)/X_r^{(S+e_r@i)}(Nn-e_r)] D_r^{(S)}(Nn-e_r), Nn_r >= 2, rho_{i,r}^{(S)}(M) = rho_{i,r}^{(S)}(M-e_s) X_s^{(S)}(M)/X_s^{(S+e_r@i)}(M), rho(0)=1, s ~= r,

and X_r^{(S)}(Nn) closed by population conservation. With K OI stations the shift becomes a K x R matrix S (row i = s_i); each OI station keeps its own D^i, rho^i and Q^i recursions driven by the common throughput X^{(S)}(Nn), and the conservation identity aggregates every station’s contribution:

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

where the LI queue Q^{(j)}_r = X_r D_{j,r} (1 + sum_s Q^{(j)}_s(Nn - e_r)) is the standard arrival-theorem term. States (S, Nn) are processed by increasing sum(Nn) so every reference lands at a strictly smaller free population.

Parameters:

Z - (1 x R)
N - (1 x R)
mu - cell array {mu_1,…,mu_K} of function handles; mu_i (n) – total service rate of station i for per-class occupancy n (1 x R). A bare function handle is accepted as the single-station shorthand.
Dli - (J x R) – (D_{j,r} = V_{j,r}/rate_{j,r}); empty or omitted when J = 0.
options - solver options (optional, currently unused)

Returns:

X - (1 x R) per-class throughput X_r = G(N-e_r)/G(N). Qoi - (K x R) per-class mean queue-length at each OI station (row i). Qli - (J x R) per-class mean queue-length at each LI queue (row j). Qdelay - (1 x R) per-class mean queue-length at the delay node (X.*Z). Soi - (K x R) per-class mean number of IN-SERVICE jobs at each OI station,

i.e. E[sir_r] with sir_r the count of 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 therefore obtained from the OI count marginal

pM_i(n|k) = (1/mu_i(n)) sum_r X_r(k) pM_i(n-e_r|k-e_r), pM_i(0|k) = 1 - sum_{n ~= 0} pM_i(n|k),

which is assembled here from the zero-shift throughputs X^{(0)}(k) already cached by the mean-value recursion above (no normalizing constant is formed). It is only computed when requested.

pfqn_gldsingle(L, N, mu, options)

Exact normalizing constant for single-class load-dependent models.

Parameters:

L – Service demand vector (Mx1).
N – Population (scalar).
mu – Load-dependent rate matrix (MxN).
options – Solver options.

Returns:

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

G=PFQN_GLDSINGLE(L,N,MU)

pfqn_conv(L, N, Z, cdscaling, options)

[G,LG] = PFQN_CONV(L, N, Z, CDSCALING, OPTIONS)

Multichain convolution algorithm for closed queueing networks with class-dependent service rates.

Implements the convolution algorithm of Sauer (1983), Section 5.2, “Computational Algorithms for State-Dependent Queueing Networks”, ACM TOCS, Vol. 1, No. 1, pp. 67-92.

The algorithm computes G(N) = (X_1 * X_2 * … * X_M)(N) where X_m(n) is the station factor at population vector n, and * denotes the multivariate discrete convolution:

A(n) = sum_{i: 0<=i<=n} B(i) * C(n-i)

For class-dependent stations, X_m(n) is computed recursively via Sauer eq. (40):

X_m(n) = (u_km / mu_km(n)) * X_m(n - e_k)

where mu_km(n) = (n_k/|n|) * beta_{m,k}(n) and beta is the DIMENSIONLESS class-dependent scaling of the service demand supplied by CDSCALING{m}: a handle of the per-class population vector n at station m, returning either a scalar (shared by every class) or a length-R vector. Equivalently

X_m(n) = (|n|/n_k) * (L(m,k)/beta_{m,k}(n)) * X_m(n - e_k),

which at beta = 1 is exactly the load-independent multinomial form, so a unit scaling means “no correction”. This is the same convention AMVA and CTMC use (effective service time ST/beta). Any saturation/cutoff is applied inside the handle.

For standard (load-independent) stations, X_m(n) reduces to the multinomial form and the convolution uses the efficient recurrence:

G_m(n) = G_{m-1}(n) + sum_r L(m,r) * G_m(n - e_r)

Parameters:

L - Service demand matrix (M x R)
N - Population vector (1 x R), must be finite (closed network)
Z - Think time vector (1 x R)
cdscaling - Cell array {M,1} of class-dependence handles beta_m (n) – empty entries denote load-independent stations
options - Solver options (optional)

Returns:

G - Normalizing constant G(N) lG - log(G(N))

pfqn_cdfun(nvec, cdscaling, classIdx)

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

Parameters:

nvec – Population state vector.
cdscaling – Cell array of class-dependent scaling functions.
classIdx – Optional class index selecting beta_{i,r} (default: 1).

Returns:

r – Scaling factor vector for each station.

R = PFQN_CDFUN(NVEC, CDSCALING, CLASSIDX)

AMVA-QD class-dependence function. 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=CLASSIDX.

CDSCALING{i} is a function handle 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

a vector of length R, i.e. the per-class scalings [beta_{i,1}(n), …, beta_{i,R}(n)], of which element CLASSIDX 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.

pas_swap2order(swap, listRate, N0)

H = PAS_SWAP2ORDER(SWAP, LISTRATE, N0)

Derive the GLOBAL placement-order DAG H of a closed two-station pass-and-swap (P&S) tandem 1->2->1 directly from its swap graph, for use with PFQN_PAS_IS.

With a non-empty swap graph the ordered-state chain is reducible (Comte & Dorsman, 2021, arXiv:2009.12299): the recurrent communicating class is the set of splits (c_{1..k}; c_{ell..k+1}) of the orderings c that are the linear extensions of a single placement partial order on the classes (their Prop.). PFQN_PAS_IS samples those orderings from H, so it needs exactly this global order – NOT the per-station placement orders consumed by the exact convolution PFQN_PAS_NC (which carry the same information one station at a time, transposed around the cycle, and are individually over-constrained for the single-order IS formulation). Feed H to PAS_PLACEMENT to obtain the precedence closure that PFQN_PAS_NC expects.

The placement order is a class-level property, independent of the per-class multiplicity, so it is extracted from the minimal single-job-per-class instance (N0 = ones(1,R)): enumerate the reachable communicating class from the all-in-queue-1 initial state; each reachable state (l1;l2) exposes the full ordering c = [l1, reverse(l2)] in D; then set H(i,j)=1 iff class i precedes class j in EVERY c in D (forced precedence). Incomparable pairs are left 0 (antichain). The result is valid for any population N.

Parameters:

swap - (R x R) – G(a,b)~=0 means class a chases class b. Empty/all-zero => plain OI: H = 0 (every ordering feasible).
listRate - cell {1 x 2} of OI microstate rate functions; listRate{m} (c) – the total service rate of queue m on the ordered prefix c. Used to prune zero-rate (non-head) completions.
N0 - (1 x R) minimal probing population; defaults to ones(1,R)

Returns:

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

See also PFQN_PAS_IS, PAS_PLACEMENT, PFQN_PAS_NC.

pas_placement(H)

[P, PLACEABLE] = PAS_PLACEMENT(H)

Placement-order logic of a pass-and-swap (P&S) / order-independent network with swap graph H. Isolates the check that decides which class orderings are feasible (adhere to the placement partial order) for PFQN_PAS_IS / PFQN_PAS_NC.

An ordering c = (c_1, …, c_ell) is FEASIBLE iff it is non-decreasing with respect to H, i.e. class a never appears before class b whenever H(b,a) ~= 0 (Comte & Dorsman, 2021, arXiv:2009.12299). Equivalently H(b,a) ~= 0 means b must be placed before a. Collecting these constraints and taking the transitive closure yields the precedence matrix

P(i,j) = 1 iff class i must be placed before class j,

so an ordering is feasible iff every class is placed only after all of its P-predecessors. H may be given as a mere Hasse diagram; the closure makes the full order explicit.

Parameters:

H - (R x R) swap-graph adjacency (H(b,a) ~= 0 forces b before a) – or all-zero graph yields P = 0 (no constraints, pure OI: every ordering feasible).

Returns:

P - (R x R) precedence closure; P(i,j)=1 iff i must precede j. placeable - function handle placeable(x) returning the row vector of class

indices that may be placed next given the remaining per-class count vector x (1 x R): those present classes with no remaining predecessor still to be placed. Used to enumerate/sample the feasible orderings and to check placement-order adherence.

See also PFQN_PAS_IS, PFQN_PAS_NC, PAS_SWAP2ORDER.

pfqn_clw(L, N, Z, m, options)

Computes the normalization constant 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 (Choudhury, Leung and Whitt, 1995).

The generating function of g(K) is (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 closed chains, i = 1..q’ indexes the distinct single-server queues with multiplicity m_i, rho_{ji} is the relative traffic intensity of chain j at queue i and rho_{j0} the aggregate relative traffic intensity of chain j at the infinite-server queues. g(K) is the coefficient of prod_j z_j^{K_j}, 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). The factors of (4.5) induce an interdependence graph on the chains (an edge whenever two chains share a queue). Removing a subset D disconnects the rest into components S_i(D), and the inversion dimension is |D| + max_i |S_i(D)| (eq. 3.2), minimized over D (eq. 3.3). With the D variables held on their contours the remaining factors split into groups with no variable in common, so each component is inverted separately and the results multiplied. The dimension reduction sets the inversion ORDER; the scaling is otherwise unchanged (Sec. 5.4). - EULER SUMMATION (Sec. 2.4, eqs. 2.20-2.22). The inner sum of (2.3) is nearly alternating, so for K_j > n+m it is replaced by the Euler sum of its first n+m+1 terms, applied twice (once for k >= 0, once for k < 0). Cost per chain drops from 2 l_j K_j to 2 l_j (n+m+1) evaluations, i.e. prod_j K_j becomes prod_j min(n+m+1, K_j) (eq. 2.26). The order m is doubled until the paper’s own estimate |E(m,n) - E(m,n+1)| falls below options.euler_tol, and the exact sum is taken once n+m reaches K_j, so the acceleration never costs accuracy: a fixed 32-term Euler sum is already 4e-4 nats off at K_j = 200 on Example 8.2.

Parameters:

L – (q’ x p) single-server relative traffic intensities, L(i,j)=rho_{ji}.
N – (1 x p) closed-chain population vector K.
Z – (1 x p) aggregate infinite-server relative intensities rho_{j0} (think-time term). Default: zeros(1,p).
m – (q’ x 1) queue multiplicities m_i. Default: ones(q’,1).
options – struct with optional fields: .l (1 x p) inner lattice parameters l_j (roundoff control),

indexed by chain. Default by inversion DEPTH: 1 at depth 1, 2 at depths 2-3, 3 deeper (Section 2.2, page 942).

.gamma (1 x p) aliasing parameters gamma_j (aliasing ~ 10^-gamma_j),
indexed by chain. Default by depth: 11, 13, 13, 15, …

.euler apply Euler summation where K_j > euler_n + euler_m.
Default true.

.euler_n number of terms summed exactly before averaging (n in
eq. 2.22). Default 11.

.euler_m starting order of the Euler averaging (m in eq. 2.22).
Default 20, so 32 terms per half-sum; doubled on demand.

.euler_tol relative tolerance on |E(m,n) - E(m,n+1)| below which the
Euler sum is accepted. Default 1e-10.

.euler_maxm largest Euler order reached by doubling. Default 160.
Beyond it the last estimate is returned rather than paying for the exact sum.

.beta (1 x 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 apply dimension reduction by decomposition. Default true. .dimred_maxd largest |D| examined when minimizing (3.3). Default 4.

Returns:

G – Normalization constant g(K). Inf if it overflows double range. lG: Natural logarithm of g(K) (always finite).

Validation: reproduces Table I (Example 8.1, p=1, every K up to 2e7) and Table II rows 1-7 (Example 8.2, p=4) to the 7 printed digits, and Table III (Example 8.3, p=11, which the reduction takes to dimension 2 and which is otherwise unreachable) to the same, its last four rows under the paper’s own scale tuning, options.beta(1) = 0.8 (0.95 on the last row; page 956). It matches exact convolution (pfqn_ca) to ~1e-9 wherever the magnitudes stay moderate. NOT reproduced: Table II rows 8-9 and Table IV, where the aliasing residual of the outer inversion exceeds the coefficient being extracted – raising gamma_1 by 7 moves the answer by exactly 7 decades, which is the signature of reading the residual rather than the coefficient. The paper computed those two tables with the model-specific ANALYTIC inner inversion of eq. (5.35) (“no l_j is involved for 2 <= j <= 11”), not with the general algorithm, so this is a limit of the scaling of Section 5 rather than of either acceleration: switching Euler summation off changes those rows in no digit at all.

cd_peak_scaling(beta, NK, K)

#ok<INUSD> CD_PEAK_SCALING Peak of a class-dependence handle over the population lattice

bmax = CD_PEAK_SCALING(beta, NK, K)

Peak of the class-dependence handle over the reachable population lattice 0 <= n(r) <= NK(r). The handle returns either a scalar (shared by every class) or a length-K vector, so the peak is taken over both the states and the classes: utilization is a per-station quantity, so the whole station shares one normalizer, as it does for max(lldscaling(ist,:)).

This is the single normalizer used to report utilization at stations with limited class dependence, U = T*S/bmax, so that every solver follows the same convention as solver_ncld does for lldscaling (U/max(lldscaling)).

laplaceapprox(h, x0)

Laplace approximation for multidimensional integrals.

Parameters:

h – Function handle to approximate integral of.
x0 – Point for Laplace approximation.

Returns:

I – Approximate integral value. H: Hessian matrix at x0. logI: Logarithm of integral value.

I = laplaceapprox(f,x0) approximates I=int f(x)dx by Laplace approximation at x0 example: I = laplaceapprox(@(x) prod(x),[0.5,0.5])

infradius_hnorm(x, L, N, alpha)

Helper function for infinite radius computation with normal CDF transformation.

Parameters:

x – Normal CDF transformation parameters.
L – Service demand matrix.
N – Population vector.
alpha – Load-dependent rate matrix.

Returns:

y – Evaluated function value for integration.

infradius_h(x, L, N, alpha)

Helper function for infinite radius computation with logistic transformation.

Parameters:

x – Logistic transformation parameters.
L – Service demand matrix.
N – Population vector.
alpha – Load-dependent rate matrix.

Returns:

y – Evaluated function value for integration.