api.npfqn

npfqn_dps_morrison(N, Z, S, w)

[Q,R,X,AUX] = NPFQN_DPS_MORRISON(N, Z, S, W)

Two-term heavy-usage asymptotic approximation for a closed queueing network with one infinite-server (think) station and one discriminatory processor-sharing (DPS) station, after

J.A. Morrison, “Asymptotic analysis of a large closed queueing network with discriminatory processor sharing”, Queueing Systems 9 (1991) 191-214.

The network is NOT product-form, so there is no normalizing constant here: the method expands the GENERATING FUNCTION of the balance equations. The substitution P(n) = <w,n> f(n) clears the DPS denominator and turns the balance recursion into a linear PDE with affine coefficients (Morrison eq. 2.5); rescaling z = 1 - xi/sqrt(N) and expanding in powers of N^(-1/2) leaves a degenerate leading operator whose kernel is the functions of the similarity variable eta, and the solvability condition along its characteristic gives an ODE for the amplitude (eq. 2.20). RESULT 1 (eq. 4.11) and RESULT 2 (eq. 4.17) are the two-term approximations returned here.

Input:: N - per-class populations (1,K), positive Z - per-class mean think time at the delay station (1,K), positive S - per-class mean service time at the DPS station (1,K), positive w - per-class DPS weights (1,K), positive
Output:: Q - mean number of class-k jobs at the DPS station (1,K) R - mean class-k sojourn time per visit to the DPS station (1,K) X - per-class throughput (1,K), by Little’s law on the think station aux - struct of intermediate quantities: the usage rho, the heavy-usage

parameter a, Morrison’s constants, the vector sigma of eq. (4.12), the W_m of eq. (3.23), and Qlead/Rlead, the leading-order (one-term) values, whose gap to Q/R measures the size of the correction.

Scaling. Morrison writes K_j = N b_j and lambda_j = N r_j g_j with N large and usage rho = sum_j b_j/g_j = 1 - a/sqrt(N). N is bookkeeping only: b, g and a all move with it and the approximation is invariant (verified numerically over four decades of N), so this function fixes N = 1, i.e. b = N_pop, g = Z./S, r = 1./Z and a = 1 - rho. The accuracy is governed by the PHYSICAL regime – large populations with rho near 1 – and not by any choice made here. rho > 1 (a < 0) is admissible: it is the saturated regime of Morrison’s appendix A, where the leading term agrees with the Mitra-Weiss fluid approximation.

npfqn_gtmtst_fluid(lambdaFuns, sFuns, muFuns, patienceCcdfs, P, T, varargin)

NPFQN_GTMTST_FLUID Time-varying open network of many-server fluid queues.

RESULT = NPFQN_GTMTST_FLUID(LAMBDAFUNS, SFUNS, MUFUNS, PATIENCECCDFS, P, T) solves an open network of M queues on [0,T]. Each queue is the Gt/Mt/st+GI fluid queue of QSYS_GTMTST_FLUID; the departure flow of queue i is routed to queue j with proportion P(i,j), and whatever is left leaves the network. LAMBDAFUNS, SFUNS, MUFUNS and PATIENCECCDFS are cell arrays with one entry per queue; P is an MxM substochastic matrix or a handle P(t) returning one.

THE NETWORK IS A FIXED POINT. The total arrival rate of queue j is

lambda_j(t) = lambda_j^0(t) + sum_i sigma_i(t) P_ij(t), sigma_i = mu_i B_i,

(eqs. 23-24) and sigma_i itself depends on lambda_i, because B_i does. The iteration starts from the external rates alone and adds one more traversal of the network per round, so the nth iterate is the fluid that has made n transitions through the network; the map is a monotone contraction, so the rates increase to the fixed point rather than oscillating.

Only the SERVICE COMPLETION flow is routed. Abandoning fluid leaves the network, which is what makes the traffic equations linear in sigma.

Options:: ‘dt’, DT - grid step, default T/2000 ‘B0’, VEC - initial fluid in service at each queue ‘w0’, VEC - initial boundary waiting time at each queue ‘tol’, TOL - sup-norm tolerance on the arrival-rate iteration, 1e-6 ‘maxIter’, K - cap on the iterations, default 100

Returns a struct with fields times, queues (a cell array of the per-queue structs of QSYS_GTMTST_FLUID), arrivalRates (the converged total rates, one row per queue), iterations and residual.

Example

fc = @(x) exp(-0.5*x); res = npfqn_gtmtst_fluid({@(t) 110+0*t, @(t) 0*t}, {100,80}, {1,1}, {fc,fc}, …

[0 1; 0 0], 60);

Reference: Y. Liu, W. Whitt (2014). Algorithms for time-varying networks of many-server fluid queues. INFORMS J. on Computing 26(1), 59-73.

See also QSYS_GTMTST_FLUID, QSYS_GGISGI_FLUID.

npfqn_feedback_elim(P, rho, varargin)

NPFQN_FEEDBACK_ELIM Near-immediate feedback elimination for RQNA.

RESULT = NPFQN_FEEDBACK_ELIM(P, RHO) computes, for each station of an open queueing network with routing matrix P and traffic intensities RHO, the probability that a departing customer returns to that station WITHOUT passing through a busier one, and the modified service description that eliminates that feedback.

WHY FEEDBACK BREAKS DECOMPOSITION. A parametric decomposition treats the arrival stream at each station as if it were renewal. Feedback destroys that badly: a customer that leaves a busy station and comes straight back arrives exactly when the station is busy, so the flow is strongly correlated with the queue it feeds. The fix is not to model the correlation but to REMOVE the feedback, by folding the repeated visits into the service time.

THE TRANSFORMATION. With feedback probability p, a customer is served a geometric number of times, so the effective service is S_p = sum_{i<=N} S_i with N geometric of mean 1/(1-p). Hence

effective mean service E[S]/(1-p),

effective service SCV p + (1-p)cs^2 (eq. 37 and the line after it),

fresh arrival rate lambda(1-p),

per-visit waiting time = (1-p) times the wait in the modified system.

The modified system has the SAME heavy-traffic limits for queue length, workload, waiting time and external departures, so the elimination is asymptotically exact rather than merely plausible.

NEAR-IMMEDIATE, NOT JUST IMMEDIATE. Feedback rarely returns a customer in one hop. What matters is whether it returns WITHOUT PASSING A BUSIER STATION: a detour through a station of lower traffic intensity is fast on the time scale of the busy station, so it behaves like immediate feedback.

Options:: ‘cs2’, VEC - service SCV per station; adds modifiedScv to the output ‘lambda’, VEC - arrival rate per station; adds modifiedRates ‘immediateOnly’, TF - keep only the self-loops P(i,i), i.e. immediate

feedback in the strict sense of Section 4.1

Returns a struct with fields feedbackProb (p-hat per station), visitInflation (1/(1-p), the mean visits per customer), modifiedScv, modifiedRates, modifiedRouting and reductionExact.

Example

% the three-station example of Dai et al. (1994), Section 6.1 of the paper P = [0 1 0; 0.5 0 0.5; 0 0.5 0]; res = npfqn_feedback_elim(P, [0.675 0.9 0.45]); res.feedbackProb % 0, 0.75, 0

Reference: W. Whitt, W. You (2022). A robust queueing network analyzer based on indices of dispersion. Naval Research Logistics 69(1), 36-56, Section 4.

See also NPFQN_TRAFFIC_IDC, SOLVER_RQNA.

mexify_npfqn

@brief MATLAB Coder script to generate MEX functions for npfqn_ module.

This script generates MEX (MATLAB Executable) versions of npfqn_ (Non-Product-Form Queueing Network) functions for improved performance.

Skipped functions (Coder-incompatible):

npfqn_nonexp_approx - Uses sn struct + SchedStrategy /: GlobalConstants OOP enums and dispatches on string method, none of which is supported by Coder.
npfqn_traffic_merge - Cell-array-of-MMAPs argument and dynamic: struct dispatching on config.merge / config.compress strings.
npfqn_traffic_merge_cs - Cell-array-of-MMAPs argument with: dynamic-typed config struct.
npfqn_traffic_split_cs - Builds variable-shape cell-of-cell output: (varargout) which Coder cannot infer.

See also CODER, CODER.CONFIG, CODER.TYPEOF, CODEGEN.

npfqn_traffic_rqt(lambda0, Gamma0, alpha0, F)

[LAMBDA,GAMMA,ALPHA] = NPFQN_TRAFFIC_RQT(LAMBDA0,GAMMA0,ALPHA0,F)

Effective arrival process perceived at each node of a single-class open queueing network under the Robust Queueing Theory (RQT) calculus. Solves the network characterization of Theorem 10 (Theorem 7 when all tail coefficients agree), which composes three operators: passage through a queue with adversarial servers leaves the uncertainty set unchanged (robust Burke, Theorem 4), superposition merges sets by Theorem 5, and thinning by a fraction f scales the rate by f and the variability by f^(-1/alpha) (Theorem 6). The resulting equations are

lambda_j = lambda0_j + sum_i lambda_i f_ij, Gamma_j = (1/lambda_j) [ 1{a0_j=ab_j} (lambda0_j Gamma0_j)^(p_j)

sum_i 1{ab_i=ab_j} (lambda_i Gamma_i)^(p_i) f_ij ]^(1/p_j),

with p_j = ab_j/(ab_j-1) and ab_j = min over the streams feeding j of their tail coefficients: the heaviest tail upstream dominates. Both are solved exactly rather than iteratively. The rate equations are the usual traffic equations, and in the variables z_j = (lambda_j Gamma_j)^(p_j) the variability equations are linear as well, so each is one linear system; ab is obtained by propagating the minimum to a fixed point.

Inputs:

LAMBDA0 - (J x 1) external arrival rate at each node, 0 where there is none GAMMA0 - (J x 1) variability parameter of each external arrival process,

which for a renewal stream is the interarrival standard deviation

ALPHA0 - (J x 1) tail coefficient in (1,2] of each external arrival process F - (J x J) routing probability matrix, F(i,j) = fraction of the jobs

leaving node i that go to node j (row sums <= 1)

Returns:: LAMBDA - (J x 1) effective arrival rate at each node GAMMA - (J x 1) effective variability parameter at each node ALPHA - (J x 1) effective tail coefficient at each node

Reference: C. Bandi, D. Bertsimas, N. Youssef (2015). Robust Queueing Theory. Operations Research 63(3), 676-700, Theorems 4-7 and 10.

npfqn_bnd_bpt(lambda0, mu, P, stationOf, c)

[ZLB, X, INFO] = NPFQN_BND_BPT(LAMBDA0, MU, P, STATIONOF, C)

First-order linear-programming relaxation of the achievable region of a multiclass open Markovian queueing network. Returns a LOWER bound ZLB on sum_r C(r)*x_r, where x_r is the mean sojourn time of class r, valid for EVERY non-idling scheduling policy.

A “class” here is a buffer with its own exponential service rate and its own Markovian routing, so a station that serves several customer types owns one class per type, and a customer type that visits a station twice owns two classes. The network is open: class r receives external Poisson arrivals at rate LAMBDA0(r) and, on completing service, becomes class r’ with probability P(r,r’) or leaves with the deficit probability 1 - sum_r’ P(r,r’).

METHOD. Uniformize the chain and let R(t) = sum_r f(r) n_r(t) for an arbitrary vector f. Writing the steady-state balance of E[R^2] gives an identity that is quadratic in f; since it holds for every f, the coefficient matrices of the two sides agree entrywise. The diagonal entries give one equation per class and the off-diagonal entries one per unordered pair, in the variables

x_r = E[T_r], the mean sojourn time of class r, I(r,l) = E[1{server sigma(r) busy with class r} * n_l], N(i,l) = E[1{server i idle} * n_l].

A third block states that the events “station i serves class r”, r in C_i, and “station i idle” are mutually exclusive and exhaustive, so their I and N terms sum to E[n_l] = lambda_l x_l. Minimizing C’x over this polyhedron is a relaxation of the true achievable region, hence a lower bound. This is the nonparametric derivation of the reference (obtained independently by Kumar and Kumar 1994); it dominates the parametric potential-function bound and needs O(K^2) variables and constraints rather than O(2^K) constraints.

EXACT ON M/M/1. With one class the LP reads mu*I11 - lambda^2*x = lambda and I11 + N11 = lambda*x with N11 >= 0, whence x >= 1/(mu-lambda) with equality.

NOT INCLUDED, DELIBERATELY. The inequality I(r,r) >= rho_r (a server busy with class r holds at least one class-r customer) is valid and would tighten the relaxation, but it is not part of the reference’s characterization, and reproducing the reference’s published bounds is the acceptance test here.

Inputs:

LAMBDA0 K x 1, external Poisson arrival rate into each class (0 if none) MU K x 1, exponential service rate of each class P K x K, P(r,r’) = probability class r becomes class r’ after

service; row sums must not exceed 1

STATIONOF K x 1, index in 1..M of the station serving each class C K x 1, objective weights (default: all ones)

Outputs:

ZLB lower bound on sum_r C(r)*x_r X K x 1, the x block of the LP optimizer. Only the objective value

is a bound; an individual X(r) is a vertex coordinate, not a bound on class r, unless C is the r-th unit vector

INFO struct with fields lambda (effective rates), rho (per-class: utilizations), rhoStation, exitflag, nvars, nrows

Reference: D. Bertsimas, I. Paschalidis, J. Tsitsiklis (1994). Optimization of multiclass queueing networks: polyhedral and nonlinear characterizations of achievable performance. Annals of Applied Probability 4(1), 43-75. See also D. Bertsimas (1995), Queueing Systems 21, 337-389, Theorem 9, which restates the same characterization.

npfqn_bnd_bgt(lambda, mu, sigma, J)

[QUB, INFO] = NPFQN_BND_BGT(LAMBDA, MU, SIGMA, J)

Piecewise-linear Lyapunov UPPER bound on the steady-state queue lengths of a multitype (deterministic-routing) multiclass Markovian queueing network, valid for EVERY work-conserving Markovian policy.

MODEL. J single-server stations; I customer types; type i arrives as a Poisson stream of rate LAMBDA(i) and passes through stages k = 1..numel(MU{i}), stage k being served at station SIGMA{i}(k) at exponential rate MU{i}(k). Class (i,k) is the buffer of type i at stage k; N = sum_i numel(MU{i}) is the number of classes.

METHOD. Solve the Down-Meyn global-stability linear program GLP[dm], eq. (25)-(28) of the reference, in the piecewise-linear Lyapunov function Phi(x) = max_j L^j’x:

L^j(i,1) lambda_i + mu(i,k) (L^j(i,k+1) - L^j(i,k)) + V_j <= -gamma
for (i,k) in station j

mu(i,k) (L^j(i,k+1) - L^j(i,k)) <= V_j for (i,k) not in j (1/(J-1)) sum_{j’ ~= j} L^j’(i,k) >= L^j(i,k) for (i,k) not in j L, V, gamma >= 0

with the convention L^j(i,Ji+1) = 0. A feasible solution with gamma > 0 certifies that EVERY work-conserving policy is stable, and a smoothed Phi is then a Lyapunov function with drift gamma/4 and an explicit exception parameter, which turns into the bound of the reference’s Theorem 4:

E[L^j’Q] <= 16 N J^2 (J-1) (Lmax+gamma)^3 / gamma^2

8 (Lmax + gamma/2)^2 / gamma =: U,

Lmax = max L^j(i,k), for every j. Since L >= 0 and Q >= 0, this gives E[Q(i,k)] <= U / max_j L^j(i,k) per class, and the geometric tail

P( L^j’Q >= B + 2(Lmax+gamma/2) m ) <= ((Lmax+gamma/2)/(Lmax+3gamma/4))^m.

THE RATES ARE RESCALED so that sum_i lambda_i + sum_{i,k} mu(i,k) = 1, the uniformization the reference imposes before Theorem 4. Queue lengths are counts and are unaffected by the time scale, so QUB is in the original units.

NORMALIZATION, WHICH THE REFERENCE LEAVES OPEN. GLP[dm] is homogeneous, and so is the bound: scaling (L,V,gamma) by t > 0 scales U by t and every denominator L^j(i,k) by t. This routine therefore fixes L^j(i,k) <= 1 and MAXIMIZES gamma, then breaks ties among gamma-optimal solutions by maximizing sum L – a degenerate optimum can otherwise leave some L^j(i,k) = 0 and report an infinite bound for a class for no reason.

THE BOUND IS LOOSE, and knowingly so: the exception parameter carries (Lmax+gamma)^3/gamma^2 and dominates as soon as J > 1. On M/M/1 (J=1, where that term vanishes) it is about 23x the exact mean queue length; on a two-station tandem it is three orders of magnitude above it. What is sharp is the STABILITY CERTIFICATE gamma > 0 and the geometric tail RATE; the constant in front of the tail is not.

Inputs:

LAMBDA I x 1, Poisson arrival rate of each type MU 1 x I cell, MU{i} is the vector of stage service rates of type i SIGMA 1 x I cell, SIGMA{i}(k) in 1..J is the station of stage k J number of stations (default: max over SIGMA)

Outputs:

QUB 1 x I cell, QUB{i}(k) upper bound on E[Q(i,k)]; Inf where the LP: optimum leaves max_j L^j(i,k) = 0
INFO struct with fields gamma, Lmax, L (J x N), V (J x 1), B (exception: parameter), U (the Theorem 4 bound on E[L^j’Q]), tailRatio (the geometric decay ratio), tailStep (2(Lmax+gamma/2)), rho (per-class load), rhoStation (J x 1), scale (the uniformization divisor), classType, classStage, classStation (N x 1 index maps)

Reference: D. Bertsimas, D. Gamarnik, J. N. Tsitsiklis (2001). Performance of multiclass Markovian queueing networks via piecewise linear Lyapunov functions. Annals of Applied Probability 11(4), 1384-1428, Section 5.1 (GLP[dm] of Down and Meyn 1997, and Theorem 4).

npfqn_traffic_split_rr(sn)

KRR = NPFQN_TRAFFIC_SPLIT_RR(SN)

Deterministic (round-robin) split degree of the departure stream of each station-class. KRR(i,r)=k>1 means that the class-r departures of station i are dispatched one-in-k by a round-robin node, so that a downstream flow carrying a fraction p of them is the k-fold convolution thinned with probability q=k*p and has SCV 1+p*(d2-k), against the Markovian 1+p*(d2-1). KRR(i,r)=1 marks an ordinary probabilistic split.

Only two topologies admit the deterministic rule: the station dispatches round-robin itself, or it feeds with probability one a router that does and whose pointer no other flow advances. Anything else falls back to k=1.

npfqn_traffic_idc(lambda0, P, c2a0, a0IdcFun, mu, cs2, sIdcFun, corrections)

ctx = NPFQN_TRAFFIC_IDC(lambda0,P,c2a0,a0IdcFun,mu,cs2,sIdcFun)

Traffic variability equations for the Robust Queueing Network Analyzer (RQNA) of W. Whitt and W. You (2018), “A Robust Queueing Network Analyzer Based on Indices of Dispersion”. Assembles and solves:

the limiting variability equations (eq. 42/44) for the asymptotic total-arrival variability parameters c2_{a,i} = I_{a,i}(Inf);

a solver (handle ctx.IaFun) for the time-dependent index of dispersion for counts (IDC) equations (eq. 40/43), returning I_{a,i}(t) for all internal arrival flows, using the default correction terms alpha_{i,j} (eq. 34) and beta_i (eqs. 38-39) and tuning function h(rho)=rho^2.

This models a single-class open queueing network of K single-server FCFS queues with Markovian routing P (P(i,j)=p_{i,j}).

Inputs (all K-vectors are column vectors, K = number of queues):: lambda0 : external arrival rate into each queue P : KxK routing matrix among queues c2a0 : asymptotic IDC (SCV) of each external arrival process a0IdcFun : handle, a0IdcFun(t) -> Kx1 external arrival IDC I_{a,0,i}(t) mu : service rate at each queue cs2 : service SCV c2_{s,i} sIdcFun : handle, sIdcFun(t) -> Kx1 service IDC I_{s,i}(t)
Output ctx: struct with fields: lambda,rho,Xi,c2a,c2d,c2aij,c2x and function handle IaFun(t) returning the Kx1 vector of total arrival IDCs I_{a,i}(t).

npfqn_sqd(sn, N, calibrationMode, serverBlockingTime, neighborMode, v1Policy, initialV1)

[X,Q,U,R] = NPFQN_BAS(SN, N, CALIBRATIONMODE, SERVERBLOCKINGTIME, NEIGHBORMODE, V1POLICY, INITIALV1)

Smith Queue Decomposition (SQD): approximate MVA for closed Blocking-After-Service (BAS) networks.

Solves a finite-buffer closed queueing network under Blocking-After-Service (manufacturing/transfer) blocking directly from its NetworkStruct. The method is an AMVA-style population recursion in which each finite-capacity station is described by a load-dependent effective service rate calibrated from an M/M/1/K blocking probability; downstream blocking is propagated through the effective routing between service stations. Delay (INF/EXT) stations are treated as infinite-capacity pure-delay nodes. Single-chain (chain-aggregated) demands only.

Returns per-station throughput X, queue length Q, utilization U, residence time R.

Originally contributed as SolverDBT by Avinash Bommareddy (Imperial College London FYP, 2026); refactored here into an sn-based API function.

npfqn_rqna_weight(t)

w = NPFQN_RQNA_WEIGHT(t) - Canonical RBM correlation weight function w*(t) used by the Robust Queueing Network Analyzer (RQNA).

w*(t) = 1 - (1 - c*(t))/(2 t), where c*(t) is the correlation function of the stationary version of canonical reflected Brownian motion (drift -1, diffusion coefficient 1),

c*(t) = 2(1 - 2t - t^2) Phi^c(sqrt(t)) + 2 sqrt(t) phi(sqrt(t)) (1 + t),

with Phi^c the standard-normal complementary cdf and phi its density. The weight is monotonically increasing with w*(0)=0 and w*(Inf)=1. Reference: W. Whitt and W. You (2018), “A Robust Queueing Network Analyzer Based on Indices of Dispersion”, eqs. (24)-(25).

Input: t - array of nonnegative time arguments Output: w - array of weights w*(t), same shape as t

npfqn_traffic_split_cs(MMAP, P, config)

Given a MMAP, produces a new array after split and class switching: P(r,(j-1)*R+s): prob that a class-r departure flows to destination j in class s out of R possible classes

empty = cellfun(@isempty, MMAP); MMAP(empty)=[]; P(:,empty)=[]; P(empty,:)=[];

npfqn_traffic_merge_cs(MMAP, prob, config): Given a cell array of n MMAPs with arrivals of R classes, produces a new MMAP after flow merging and class switching prob((i-1)*R+r,s): prob that a class-r arrival from MMAP i switches to class-s

npfqn_traffic_merge(MMAP, config): Given a cell array of n MMAPs with arrivals of R classes, produces a new MMAP after flow merging

npfqn_nonexp_approx(method, sn, ST, V, SCV, T, U, gamma, nservers): handler for non-exponential service and arrival processes