solvers.FLD

fluid_qsys_horizon(options)

[BOOL, REASON, T0, T1] = FLUID_QSYS_HORIZON(OPTIONS)

@brief The integration window of a time-varying single-station fluid limit.

The time-varying limits (‘mol’, ‘mtginf’, ‘tvms’) report a TRAJECTORY, so an infinite or absent upper end of options.timespan leaves them nothing to report; the fluid analyzer applies the same rule to its own timespan.

A horizon is a solver OPTION and not a model feature, so the feature registry has no name for it and SolverFLD.supportsModelMethod has to ask this predicate directly. SOLVER_FLUID_QSYS_ANALYZER asks the same one on the solve path, which is what keeps the report and the run from disagreeing about whether a method can be asked for.

@param options solver options carrying the timespan @return bool true when the window is a finite non-empty interval @return reason the refusal, or ‘’ when BOOL is true @return t0 lower end of the window (0 when unset) @return t1 upper end of the window (meaningful only when BOOL is true)

solver_fluid_dae(sn, options)

[QN, UN, RN, TN, XVEC_IT, QNT, UNT, TNT, XVEC_T, T, ITERS, RUNTIME, MOMENTS] = SOLVER_FLUID_DAE(SN, OPTIONS)

Differential-algebraic formulation of the min-normal closure, backing options.method=’dae’.

SOLVER_FLUID_MOMENTS already solves a differential system (the mean) coupled to an algebraic one (the covariance). It solves them by SUCCESSIVE SUBSTITUTION: integrate the mean ODE to its fixed point at a held variance, solve the Lyapunov equation there, extract sigma2, repeat, up to 20 times and only to CoarseTol. This method states the same closure as one system and solves it as one system. Nothing about the closure changes – the drift, the rate factors and the Lyapunov equation are the same functions, taken unmodified from FLUID_MOMENT_TERMS and FLUID_LYAPUNOV – only the way the coupled equations are discharged.

Two modes, chosen by the horizon, because the DAE is genuinely different in each:

STEADY STATE (options.timespan(2) infinite, the usual case)

Solve the algebraic system

0 = D*r(x, sigma2) drift residual, nstate rows 0 = C*x - Nchain population conservation, one row

per closed chain

0 = sigma2 - sigmaOf(x,sigma2) closure consistency, one row per

closable station

simultaneously by damped Newton. Convergence is quadratic near the root instead of the linear rate of substitution, and the answer is converged to options.tol rather than to the outer CoarseTol.

THE FIXED POINT IS NOT FOUND BY INTEGRATING TO IT. Integrating a stable ODE until it stops moving is a poor way to solve f(x)=0: the cost is set by the slowest mode of the model rather than by the accuracy wanted, which is exactly why the stiff models are expensive here. A seed trajectory is still integrated once, cheaply and at the first-order closure, because Newton needs a starting point inside the basin; everything after that is algebraic.

TRANSIENT (finite horizon)

Integrate the index-1 DAE

d/dt x = D*r(x, sigma2(t)) differential, kept rows 0 = C*x - Nchain algebraic, one row per chain d/dt Sigma = A*Sigma + Sigma*A’ + Q differential

with a singular mass matrix. This is the ONLY fluid method other than ‘kp’ that produces a time-varying second moment, and unlike ‘minnormal’ – which evaluates the whole transient at the single STATIONARY variance – the variance here is the one the trajectory actually had at each instant.

WHAT THE ALGEBRAIC CONSTRAINT BUYS. Population conservation currently holds only to integrator tolerance: it is a consequence of the drift (the rows of D sum to zero on a closed chain), never an equation. Writing it as a constraint has two effects. It is enforced to solver tolerance rather than accumulated. And it is what makes the Newton system solvable at all: the drift Jacobian is singular along exactly the conserved directions – the same singularity FLUID_LYAPUNOV works around by projecting onto range(D) – so the constraint rows supply the missing rank instead of a pseudo-inverse hiding it.

WHY THE COVARIANCE IS NOT A NEWTON UNKNOWN. The obvious “solve everything at once” reading puts Sigma in the unknown vector. That is a trap: Sigma is nstate^2 entries, so the Jacobian is (nstate + nstate^2)^2 and the work is quartic – strictly worse than the cubic Lyapunov solves it would replace, and it would tighten the state cap rather than relax it. Sigma is LINEAR in itself for a held x, so it is eliminated by one Lyapunov solve per residual evaluation and only sigma2, M numbers, joins x in the unknown vector.

NO SIGMA2=0 SEED, SO NO KINK WORKAROUND. SOLVER_FLUID_MOMENTS must start its alternation at sigma2=0, where min(n,c) has no derivative, and a saturated model’s first-order fixed point lands on that kink by construction; it carries a two-sided-Jacobian probe to decide hyperbolicity side-independently there. The simultaneous solve never adopts sigma2=0 as an iterate, so that probe is not needed. It does NOT rescue a fixed point that genuinely sits on the kink at the converged variance – that is a real continuum of equilibria and no formulation removes it.

Parameters:

• sn - NetworkStruct

• options - solver options; options.tol sets the Newton tolerance, – options.config.dae_maxstate caps the simultaneous solve (default 100) and options.config.dae_maxcov caps the transient covariance (default 25)

Returns:

as SOLVER_FLUID_MOMENTS, plus MOMENTS.Sigmat / MOMENTS.QVart / MOMENTS.tvar on the transient route

See also SOLVER_FLUID_MOMENTS, FLUID_MOMENT_TERMS, FLUID_LYAPUNOV.

solver_fluid_analyzer(sn, options)

[QN, UN, RN, TN, CN, XN, T, QNT, UNT, TNT, XVEC, iter, aoiResults] = SOLVER_FLUID_ANALYZER(QN, OPTIONS)

fluid_petri_theta(terms, x, s2)

TH = FLUID_PETRI_THETA(TERMS, X, S2)

The closed enabling term of every transition mode, and its derivative.

A mode fires at a rate proportional to the number of BINDINGS the marking supports, min over its input arcs of m_a/w_a, capped at its server count and zeroed by any inhibitor arc that has reached its threshold. Under the Gaussian marginal that whole expression is

theta_j = ( prod_b Phi((thr_b - m_b)/sd_b) ) * E[ min_a(m_a/w_a), c_j ]

where the many-argument min is closed by FLUID_MINMULTI_CLOSURE and the inhibitor INDICATOR is closed by the normal CDF. Both collapse to their first-order form at zero variance – Phi becomes 1{m_b < thr_b} and the min closure becomes min() – so the mean-field limit is one code path, not two.

THE INHIBITOR GATE IS WHY A PETRI NET NEEDS A SMOOTHED CLOSURE AT ALL, quite apart from the accuracy argument. 1{m_b < thr_b} is a step, and a Newton solver has no derivative to descend on a step: the drift residual would be piecewise constant in m_b and the iterate would either sit still or chatter across the threshold. Phi((thr_b - m_b)/sd_b) is the same function smoothed by the marginal the closure already carries, so the Jacobian is defined everywhere the variance is positive.

THE VARIANCES ARE UNKNOWNS, NOT FUNCTIONS OF X. S2 holds one entry per TERMS.covPairs row, pinned by its own consistency row in the DAE, so the derivative below is with respect to the MEANS only – exactly how FLUID_DRIFT_JACOBIAN differentiates the queueing closure.

Parameters:

• terms - FLUID_PETRI_TERMS output

• x - state vector [marking; phase distributions]

• s2 - (npair x 1)

Returns:

th - struct with fields

theta - (nmodes x 1) closed enabling term dslot - (nmodes x 1) cell of coordinate indices dval - (nmodes x 1) cell of dTHETA/dX at those coordinates dep - (nmodes x 1) marking-dependent firing multiplier depslot - (nmodes x 1) cell of coordinate indices depval - (nmodes x 1) cell of dDEP/dX at those coordinates

See also FLUID_MINMULTI_CLOSURE, FLUID_PETRI_RATES, FLUID_PETRI_JACOBIAN.

fluid_petri_terms(sn, options)

TERMS = FLUID_PETRI_TERMS(SN, OPTIONS)

Event-based representation of the fluid marking process of a stochastic Petri net, in the contract SOLVER_FLUID_PETRI solves and FLUID_MOMENT_TERMS supplies for a queueing network.

A GSPN IS ALREADY A DENSITY-DEPENDENT MARKOV POPULATION PROCESS, which is the object the moment-closure family of SolverFLD is built on: the marking is the population, a transition mode is a reaction, its incidence column is the jump, and the rate law lambda*min(enabling degree, servers) is the same min() non-linearity the min-normal closure exists to smooth. Nothing about the closure changes here; only where the drift comes from.

dx/dt = D * r(x, Sigma, phi, mu)

THE STATE, x = [ m ; z ].

m(p,k) token mass of class k at place p. One coordinate per (place,

class) pair that some arc touches, that the initial marking loads, or that a Source feeds; a pair nothing reaches is dropped rather than carried as a null direction of the Newton system.

y(j,h) the number of mode-j servers running in phase h, for a mode whose

firing time has more than one phase. Their SUM is not free: the ENABLE synchronization latches it instantaneously to min(enabling degree, servers), so the latch is an ALGEBRAIC row with one free-sign unknown mu_j – the net rate at which servers start or stop – and the phase split evolves differentially.

CARRYING THE DISTRIBUTION INSTEAD OF THE COUNT LOOKS TIDIER AND IS WRONG. Writing y = theta(m)*z and integrating z removes mu, but the resulting z equation is missing the (dtheta/dt)/theta*(pie-z) term and carries a spurious factor theta, so it agrees with the count form only AT a fixed point. It also hides the coupling from the marking into the firing rate, which the covariance needs: on a closed cycle with Erlang(2) firing the count form reads 1.4564 against the exact 1.4531, the distribution form 1.4357.

There is no Source coordinate: an exogenous arrival is a CONSTANT-propensity event depositing one token, as it is in the NRM SPN runner, which is the correct exogenous Poisson event once the source pool is projected away (see FLUID_MOMENT_TERMS for the same projection in the queueing case). There is no Sink coordinate either: a firing arc into a sink is mass leaving the net.

THE EVENTS, one column of D each:

kind 1 a firing of mode j out of phase h into phase h’. The marking moves

by the incidence column post-pre; for a multi-phase mode the phase moves by D1(h,h’), which is the PH restart alpha and the MAP’s landing phase in the same expression, since a PH’s process pair has D1 = s*alpha.

kind 2 an internal phase change of mode j, D0(h,h’): the marking is

unchanged, so this column contributes nothing to the diffusion of the marking.

kind 3 an exogenous arrival, constant rate lambda*p, +1 token. kind 4 a firing of an IMMEDIATE mode, at the algebraic flow phi_j.

THE RATES:

theta_j(x) = g_inh_j(x) * E[ min_a( m_a/w_a ), c_j ] enabled bindings kind 1 = theta_j * z(j,h) * D1_j(h,h’) * g_dep_j(m) kind 2 = theta_j * z(j,h) * D0_j(h,h’) kind 3 = lambda * p kind 4 = phi_j

with the many-argument min closed by FLUID_MINMULTI_CLOSURE and the inhibitor indicator 1{m_b < thr} closed by the normal CDF Phi((thr-m_b)/sigma_b), which collapses to the hard indicator at sigma = 0 exactly as FLUID_MIN_CLOSURE does, so the first-order limit needs no second code path. For a single-phase mode z is absent and theta_j is the running-server count itself.

THE COVARIANCE LIVES ON THE MARKING ONLY. z is a distribution, so its fluctuations are O(1) rather than O(sqrt(N)) and the linear noise approximation does not apply to it at any scale – the same reason FLUID_MOMENT_TERMS projects the EXT source pool away. TERMS.covIdx therefore holds the marking coordinates and nothing else.

Parameters:

• sn - NetworkStruct holding Place and Transition nodes

• options - solver options

Returns:

terms - struct; see the field list assembled at the end of this function

See also SOLVER_FLUID_PETRI, FLUID_MINMULTI_CLOSURE, FLUID_MOMENT_TERMS.

fluid_petri_rates(terms, x, s2, phi, mu, th)

[R, TH] = FLUID_PETRI_RATES(TERMS, X, S2, PHI, MU) [R, TH] = FLUID_PETRI_RATES(TERMS, X, S2, PHI, MU, TH)

The rate of every event column of the fluid Petri drift.

kind 1 firing of mode j, phase h -> h’

single phase : rateBase * theta_j * dep_j multi phase : rateBase * y(j,h)

kind 2 internal phase change of mode j: rateBase * y(j,h) kind 3 exogenous arrival: rateBase, a constant kind 4 firing of an IMMEDIATE mode: PHI(j), an algebraic unknown kind 5 the server latch of a multi-phase mode: MU(j), a free-sign

algebraic unknown depositing at the firing process’s entry vector

THETA_J IS THE NUMBER OF RUNNING SERVERS a marking supports, closed by FLUID_PETRI_THETA. A SINGLE-PHASE mode holds no coordinate of its own, so its rate is theta_j directly. A MULTI-PHASE mode carries one coordinate per phase holding the servers sitting in it, and its rates are LINEAR in those: the marking reaches them through the latch row sum_h y(j,h) = theta_j, not through the rate.

THE MARKING-DEPENDENT MULTIPLIER APPLIES TO A FIRING ONLY, and only to a single-phase mode: SETFIRINGRATEDEPENDENCE accepts it for exponential timing alone, so there is no multi-phase clock for it to speed up.

Parameters:

• terms - FLUID_PETRI_TERMS output

• x - state vector

• s2 - the closure covariance entries

• phi - (nimm x 1)

• mu - (nlatch x 1)

• th - optional precomputed FLUID_PETRI_THETA result

Returns:

r - (nev x 1) rates th - the enabling terms, so a caller taking the Jacobian next reuses them

See also FLUID_PETRI_THETA, FLUID_PETRI_JACOBIAN, SOLVER_FLUID_PETRI.

fluid_petri_jacobian(terms, x, s2, phi, th)

[A, JR] = FLUID_PETRI_JACOBIAN(TERMS, X, S2, PHI) [A, JR] = FLUID_PETRI_JACOBIAN(TERMS, X, S2, PHI, TH)

Drift Jacobian of the fluid Petri net, A = D * dR/dX.

This is what the Lyapunov equation of the linear noise approximation is written about, so it has to be the derivative of the SAME rate vector FLUID_PETRI_RATES returns, not a finite-difference stand-in: a covariance solved about an inconsistent Jacobian is not the covariance of anything.

kind 1, single phase : rateBase * ( dTHETA*dep + THETA*dDEP ) kind 1/2, multi phase: rateBase * e_y(j,h), since the rate is LINEAR in the

servers sitting in that phase and reads no marking

kind 3 : zero, the arrival rate is a constant kind 4/5 : zero, an immediate flow and a server latch are

independent unknowns of the DAE, not functions of x

THE VARIANCES ARE HELD, as they are in FLUID_PETRI_THETA and in the queueing twin FLUID_DRIFT_JACOBIAN: the closure covariance is pinned by its own consistency row, so it does not move along a derivative with respect to x.

Parameters:

• terms - FLUID_PETRI_TERMS output

• x - state vector

• s2 - the closure covariance entries

• phi - immediate firing flows (unused in the derivative, kept for symmetry)

• th - optional precomputed FLUID_PETRI_THETA result

Returns:

A - (nstate x nstate) drift Jacobian Jr - (nev x nstate) rate Jacobian

See also FLUID_PETRI_RATES, FLUID_LYAPUNOV, FLUID_DRIFT_JACOBIAN.

fluid_petri_immediate(terms, x, imm)

IMM = FLUID_PETRI_IMMEDIATE(TERMS, X) IMM = FLUID_PETRI_IMMEDIATE(TERMS, X, IMM)

The active set of the IMMEDIATE transitions, and the equations that pin their flows.

AN IMMEDIATE TRANSITION HAS NO RATE, so it cannot be a column of the drift with a rate factor like every other event. Its fluid limit is a singular perturbation: the transition fires infinitely fast, so a marking that enables it cannot persist, and what survives in the limit is a FLOW. That flow is an algebraic unknown of the DAE, pinned by the constraint that the transition’s binding input place holds no mass:

phi_j >= 0, x_b = 0 for the coordinate b that binds mode j

which is the complementarity condition of the limit, and is the same object the capacity caps already solve as an active set – an inequality has no residual for a Newton solver, so the loop iterates over WHICH constraints bind and each pass is an equality-constrained solve.

TWO MODES DRAINING ONE PLACE NEED ONE MORE EQUATION THAN THE PIN GIVES, and the extra equation is the GSPN’s own conflict rule: among the enabled immediate modes of highest firing priority the branch is taken in proportion to the firing weights, so

phi_j * weight_l = phi_l * weight_j

for the members of one conflict group, and phi = 0 for a member of lower priority. That is exactly what the exact engines do at a vanishing marking (SOLVER_SSA_NRM’s SPNCOLLAPSE, STATE.AFTERGLOBALEVENT’s immediate branch), carried over to a continuous flow.

THE COUNT IS SQUARE BY CONSTRUCTION. Each active mode is assigned exactly one binding coordinate; the pinned coordinates are the image of that assignment, so V pins plus (F - V) ratio rows is F equations for F flows. A coordinate that happens to be empty without being anybody’s binding coordinate is not pinned – its own drift row determines it – and if the solve then pushes it negative, the caller REBINDS a mode to it and solves again, which is the active-set move.

Parameters:

• terms - FLUID_PETRI_TERMS output

• x - state to read the assignment off (the seed, or the last iterate)

• imm - a previous assignment to refresh; omitted to initialise

Returns:

imm - struct with fields

n - number of immediate modes active - (n x 1) logical, whether the mode carries a flow unknown bind - (n x 1) marking coordinate pinned for that mode, 0 if none pins - the distinct pinned coordinates rows - (nrow x 1) struct array describing each equation:

.kind ‘pin’ | ‘ratio’ | ‘zero’ .a, .b, .wa, .wb as the kind requires

See also SOLVER_FLUID_PETRI, FLUID_PETRI_TERMS.

fluid_petri_constraints(sn, terms)

CON = FLUID_PETRI_CONSTRAINTS(SN, TERMS)

Every finite place capacity as a linear constraint on the fluid marking.

Arow * x <= b

with one row per finite total capacity (SETCAPACITY on a Place) and one per finite per-class capacity (SETCLASSCAPACITY). This is the same family of rows the queueing DAE builds in FLUID_CAPACITY_CONSTRAINTS, restricted to what a Petri net can declare: a place has no scheduling, no region and no drop rule of its own, so there is exactly one gate and it is the same one the exact engines apply.

THE GATE IS A LOSS ON THE DEPOSIT. LINE loses the tokens a firing would push past a place’s capacity – the NRM’s APPLYPLACECAPS clamps the just-deposited slots after every firing, which is JMT’s and the CTMC’s semantics for a bounded place (an M/M/1/1 place at rho = 0.5 holds 1/3 of a token, not 1). The fluid analogue scales the DEPOSIT leg of every event that adds mass to the capped place, and leaves the removal leg alone: the firing still happens, and only the mass that does not fit is lost.

CONSERVATION AND LOSS CANNOT BOTH HOLD. A row that binds destroys mass, so any conserved quantity supported on the capped coordinates stops being conserved. SOLVER_FLUID_PETRI therefore drops those conservation rows for as long as the cap is active, which is the honest statement of what the model declares: a closed net whose place drops tokens is not closed while the cap binds. In practice this only arises where it should – an open net bounded to make it ergodic (the M/M/1/K place) has no conserved quantity to lose in the first place, since its arrivals already broke them.

Parameters:

• sn - NetworkStruct

• terms - FLUID_PETRI_TERMS output

Returns:

con - struct with fields

A - (ncon x nstate) constraint rows b - (ncon x 1) right-hand sides label - one description per row cover - (ncon x nstate) logical, the coordinates the row caps, i.e.

the deposits its multiplier throttles

See also SOLVER_FLUID_PETRI, FLUID_CAPACITY_CONSTRAINTS.

fluid_petri_conservation(terms)

CONS = FLUID_PETRI_CONSERVATION(TERMS)

The conserved quantities of a fluid Petri net, as equations.

Everywhere else in SolverFLD conservation is a CONSEQUENCE of the drift and therefore holds only to integrator tolerance. The DAE form states it, and for a Petri net the right statement is not “one row per closed chain” – a net has no chains – but the left null space of the jump matrix:

u’ D = 0 => u’ x is constant along every trajectory

On the marking coordinates those u are exactly the net’s P-INVARIANTS (an S-invariant is a non-negative left null vector of the incidence matrix, which is what D restricted to the places is). On a mode’s phase block the all-ones vector is one of them, and it is the statement that the phase coordinates are a DISTRIBUTION: sum_h z(j,h) = 1. Both come out of the same null space, so the phase normalisation needs no separate row and no separate mechanism.

AN OPEN NET LOSES THE ROWS ITS ARRIVALS BREAK, automatically and for the right reason: the arrival columns are part of D, so a u that an arrival moves is not in the null space and never becomes a constraint.

THE BASIS IS RATIONAL, not orthonormal. D is integral (arc multiplicities and unit phase moves), so NULL(.,’r’) returns exact rational rows through RREF, which keeps each constraint readable as a conservation statement about named places instead of an arbitrary orthogonal mixture of them.

Parameters:

terms - FLUID_PETRI_TERMS output

Returns:

cons - struct with fields

C - (ncons x nstate) conserved directions N - (ncons x 1) their value at the initial marking leak - max|C*D|, zero by construction; a positive value means the

basis and the jump matrix disagree

label - one description per row

See also SOLVER_FLUID_PETRI, FLUID_PETRI_TERMS, SPN_SINVARIANTS.

fluid_petri_applicable(sn, options)

[OK, REASON] = FLUID_PETRI_APPLICABLE(SN, OPTIONS)

Whether the fluid Petri route can answer this model, and why not when it cannot.

The route covers the class of nets the exact engines cover – Place, Transition, enabling, inhibiting and firing arcs, timed and immediate modes, single/k/infinite server firing, phase-type and MAP firing times, marking-dependent firing rates, bounded places, and a Source/Sink pair for an open net. What it does NOT cover is named here rather than discovered inside the solve, so a refused model says which declaration it was refused for.

A QUEUEING STATION IS THE ONE STRUCTURAL EXCLUSION. A net whose tokens also visit a Queue or a Delay is two formalisms at once, and LINE has no reference semantics for the hand-off: the Petri arcs and the routing matrix would each describe part of the movement and nothing pins how a token becomes a job. The same goes for a QUEUEING PLACE, whose embedded queue only SolverLDES simulates.

Parameters:

• sn - NetworkStruct

• options - solver options

Returns:

ok - true when the model can be solved by SOLVER_FLUID_PETRI reason - the refusal, empty when OK

See also SOLVER_FLUID_PETRI, FLUID_DAE_APPLICABLE.

fluid_minmulti_closure(mu, S, c)

[H, G, V] = FLUID_MINMULTI_CLOSURE(MU, S, C)

Min-normal moment closure of E[min(X_1,…,X_A,c)] for jointly normal X.

FLUID_MIN_CLOSURE closes the TWO-argument min() that a queueing station needs, min(n_i,c_i). A Petri-net transition mode needs the many-argument one: its enabling degree is

e(m) = min_a ( m_a / w_a )

over every input arc a, and the rate law then caps that at the mode’s server count. There is no closed form for the expectation of a min of more than two correlated normals, so this function uses the recursion of Clark (1961): the running min is replaced at each step by the normal with its exact first two moments, and the next argument is folded in with the exact bivariate formulas. Each step is FLUID_MIN_CLOSURE’s own expression, extended with the second moment and the cross-covariances Clark needs to carry the recursion forward:

th^2 = Var[Z] + Var[X] - 2*Cov[Z,X], al = (E[Z]-E[X])/th, p = Phi(-al) E[W] = E[Z]*p + E[X]*(1-p) - th*phi(al) E[W^2]= (E[Z]^2+Var[Z])*p + (E[X]^2+Var[X])*(1-p) - (E[Z]+E[X])*th*phi(al) Cov[W,Y] = Cov[Z,Y]*p + Cov[X,Y]*(1-p)

with W = min(Z,X). The deterministic cap C is folded in last, by FLUID_MIN_CLOSURE itself, so a single-arc mode reduces EXACTLY to the closure the queueing methods already use and no second code path exists for it.

THE RECURSION IS ORDER DEPENDENT, as Clark’s approximation always is: only the first two moments of the running min are kept, so folding the arcs in a different order gives a slightly different answer. The order here is the caller’s, i.e. increasing state coordinate, which is fixed by the layout and therefore reproducible.

THE VARIANCES ARE HELD. The derivative G is taken with respect to the MEANS only, exactly as FLUID_DRIFT_JACOBIAN differentiates the queueing closure: the covariance is a separate unknown of the DAE, pinned by its own consistency row, not a function of the mean along the Newton step.

Parameters:

• mu - (A x 1)

• S - (A x A)

• c - deterministic cap (the mode’s server count)

Returns:

h - E[min(X_1,…,X_A,c)] g - (A x 1) dH/dMU(a), the probability that arc a is the binding one v - Var[min(X_1,…,X_A)] before the cap, carried for the caller’s report

– Reference C. E. Clark, “The greatest of a finite set of random variables”, Operations Research 9(2):145-162, 1961.

See also FLUID_MIN_CLOSURE, FLUID_PETRI_TERMS, SOLVER_FLUID_PETRI.

fluid_share_closure(x, wv, C)

[S, DS, CN, DCN] = FLUID_SHARE_CLOSURE(X, WV, C)

Second-order closure of the capacity share of a sharing discipline. A DPS station gives coordinate j the fraction

S_j = w_j X_j / sum_m w_m X_m

of its capacity, and PS is the same expression with unit weights. The first-order fluid closure evaluates that ratio at the mean, which is not E[S_j]: the map is a ratio, so Jensen’s inequality biases it towards the coordinates carrying the LARGER weight. Writing u_j = w_j X_j and v = sum_m u_m, the delta method gives

E[u_j/v] = mu_j/v - Cov(u_j,v)/v^2 + mu_j*Var(v)/v^3 + O(sigma^3)

evaluated at the means. The correction is exactly capacity-conserving: sum_j Cov(u_j,v) = Cov(v,v) = Var(v), so the two correction terms cancel in the sum and sum_j S_j = 1 identically, as it must for a work-conserving discipline. That identity is the invariant to check on any change here.

The expansion is local and fails when v is small against its own standard deviation, where the exact expectation is a Cauchy-like integral with no finite mean. There a raw share can come out negative; it is clipped at zero and the survivors renormalised, which preserves the conservation identity.

THE EXPANSION CARRIES ITS OWN CONVERGENCE RATIO, and both corrections below are admitted only while that ratio is below one. Every second-order term here is a term of the series for E[1/v], whose successive terms are in the ratio

ratio = Var(v)/v^2

so the truncation is meaningful below 1 and the terms GROW above it. Nothing in the algebra notices: at a near-empty station the corrections are simply evaluated far outside the region where they mean anything, and come back larger than the quantity they correct. On layer 7 of test_LQN_13 (arbitraryMultiplicity) the station holds 0.5388 jobs with Var(N) = 2, so the ratio is 6.9, and Cov(S_j,N) came back at -1.86 where 0.71 is all any joint distribution can produce: S_j lies in [0,1], so its variance is at most 1/4 and Cauchy-Schwarz caps the covariance at sqrt(Var(N))/2. The drift that follows is not integrable – one mean solve took 287306 drift evaluations and 26 s against 774 and 0.03 s for the first-order one, and the moment fixed point runs up to ITER_MAX of those.

LOCAL_EXPANSION_WEIGHT therefore scales both corrections by a factor that is exactly one while the ratio is at most 1, falling smoothly to zero by 4 (a standard deviation twice the mean, where a non-negative v has no mass left near it), and is C^1 at both ends so the drift stays differentiable and the integrator keeps its step. Outside the region the closure degrades to the first-order share u/v, which is what SOLVER_FLUID_CLOSING answers with. Note the factor multiplies the SHARE correction and CN together, so the cancellation that makes the joint closure exact at an unsaturated station survives it identically: there S_j*N = X_j whatever the factor is.

Parameters:

• x - (n x 1)

• wv - (n x 1)

• C - (n x n) – first-order (plug-in) share

THE SHARE IS ONLY HALF OF THE RATE. What a station clears is S_j*psi(N), and the two factors are correlated through N, so E[S_j*psi] is not E[S_j]*E[psi]. CN returns the missing Cov(S_j, N), from which ODE_RATES_CLOSING_FACTORS builds the joint closure E[S_j*psi] = E[S_j]*E[psi] + psi’(n)*Cov(S_j,N). That term is what makes the product exact where min() is the identity: at an unsaturated station psi = N and psi’ = 1, and S_j*N = X_j identically, so the two second-order corrections must cancel and leave the plain mean. Separately closed they do not, and a five-server station holding 0.77 jobs reported a 12.4% queueing delay that cannot exist. CN sums to zero over j, because the shares sum to one at every point and Cov(1,N) = 0, so the joint closure is work-conserving for exactly the reason the plain share closure is.

Parameters:

• x - (n x 1)

• wv - (n x 1)

• C - (n x n) – first-order (plug-in) share

Returns:

s - (n x 1) expected shares, summing to one ds - (n x n) Jacobian ds_j/dx_m, with C held fixed cn - (n x 1) Cov(S_j, N) at the means, summing to zero dcn - (n x n) Jacobian dcn_j/dx_m, with C held fixed

See also ODE_RATES_CLOSING_FACTORS, FLUID_DRIFT_JACOBIAN, FLUID_MIN_CLOSURE.

fluid_capacity_closure(n, c, s2, lldrow, isInf)

[H, DH, D2H] = FLUID_CAPACITY_CLOSURE(N, C, S2, LLDROW, ISINF)

Moment closure of the station capacity term psi(X) and of its derivative.

Every scheduling branch of ODE_RATES_CLOSING_FACTORS scales the coordinates of a station by psi(n_i)/n_i, where psi is how much work the station clears at population n_i:

psi(n) = min(n,c) * alpha(n) at a queueing station psi(n) = n * alpha(n) at an infinite server

and alpha is the limited load-dependent scaling (1 when the station has none). This function returns E[psi(X)] for X ~ Normal(n, s2) together with d/dn E[psi(X)], which is what the drift and the Jacobian need. S2 = 0 gives psi(n) itself, so the first-order closure is the same code path.

Without load dependence the expectation is the closed form of FLUID_MIN_CLOSURE. With a tabulated alpha, FLUID_LLD_SCALING makes alpha piecewise linear on the integer lattice, so psi is piecewise QUADRATIC with breakpoints at the integers and at c, and the expectation is integrated segment by segment against the normal density using the truncated moments M0, M1, M2. The derivative is E[psi’(X)] by differentiation under the integral sign, valid because psi is Lipschitz.

The integration must be exact, NOT a fixed quadrature rule. Gauss-Hermite with fixed nodes applied to a piecewise-linear integrand does not smooth its kinks, it relocates them: the resulting estimate of E[psi] is itself piecewise linear in n, so its second derivative is zero almost everywhere and FLUID_REFINE_MEANFIELD silently returns a null correction. That is what the segment-wise closed form below avoids.

psi is extended by zero below n = 0, the only physically admissible continuation, since a station holding no jobs clears no work.

Parameters:

• n - mean population at the station

• c - number of servers (ignored when ISINF)

• s2 - population variance, 0 for the first-order closure

• lldrow - (1 x lldlimit)

• isInf - true at an infinite-server station

Returns:

h - E[psi(X)] dh - d/dn E[psi(X)] d2h - d2/dn2 E[psi(X)], which FLUID_DRIFT_JACOBIAN needs for the joint

product closure of the share and the capacity. psi is only piecewise smooth, so this is E[psi’’(X)] in the distributional sense: the segment-wise quadratic term PLUS an atom at every breakpoint where psi’ jumps, weighted by the normal density there. Dropping the atoms would report d2h = 0 for the pure min(), whose curvature is carried by nothing else.

See also FLUID_MIN_CLOSURE, FLUID_LLD_SCALING, ODE_RATES_CLOSING_FACTORS.

fluid_fixed_point_guard(ode_h, fpRate)

FCN = FLUID_FIXED_POINT_GUARD(ODE_H, FPRATE)

An ODE OutputFcn that halts a window whose state has SETTLED MID-FLIGHT.

WHY IT EXISTS. SOLVER_FLUID_ITERATION already short-circuits a window ENTERED at its fixed point, in closed form, because a stiff step controller handed a state it is already at cannot pick a step. That test is taken once, at the window’s first instant, and a window that reaches the fixed point AFTER its first step is left to grind out the rest of its span at the reciprocal of the FASTEST rate while the span runs to 10*iter/min_rate, set by the SLOWEST.

WHY IT IS AN OUTPUTFCN AND NOT AN EVENTS FUNCTION. The obvious slot is odeset(‘Events’), which costs nothing when it does not fire and is free of the conservation guard. It is the wrong one: the default accurateStiffOdeSolver is @lsoda_accurate_stiff, and LSODA_ODESOLVE implements the OutputFcn init/step/ done protocol but has NO Events support – so an Events-based check would work under @ode15s, do NOTHING on the default stiff path, and fail silently there. FLUID_OUTPUTFCN_CHAIN composes this guard with FLUID_CONSERVATION_GUARD so the one slot both integrators honour can carry both tests.

THE THRESHOLD IS ROUND-OFF, NOT THE SOLVER TOLERANCE, and it is the same GlobalConstants.Zero the entry test uses, deliberately: a state that merely satisfies the window loop’s driftTol (1e-4) is still moving, and cutting the window there was measured to shift results by 1.8e-5. Below Zero the drift is zero to double precision, so holding the state for the rest of the span is exact rather than approximate.

THE CALLER MUST STILL PAD THE TAIL. Halting through an OutputFcn returns a SHORT final time, which SOLVER_FLUID_ITERATION reads as the conservation guard and turns into ‘LINE:FluidNonHyperbolic’. LOCAL_SETTLE_TAIL re-tests the residual at the returned end state and, when it is settled, extends the trajectory to the requested instants at that state – so a settled window is indistinguishable from one that ran to its endpoint, and only a genuine conservation excursion still reads as short.

Parameters:

• ode_h - the drift, called as ode_h (t, y) – window is integrating, and AUTONOMOUS, which is what the caller’s arming condition (no RT_BREAKS) already guarantees

• fpRate - the slowest exit rate, normalising the residual to a rate ratio so – the test is scale-free in time as well as in population

Returns:

fcn - handle for odeset(‘OutputFcn’, …). Returns status 1 to halt.

See also SOLVER_FLUID_ITERATION, FLUID_CONSERVATION_GUARD, FLUID_OUTPUTFCN_CHAIN.

fluid_capacity_legs(terms, gates, stg, staged, active, x, sg, r, mult, stagedFlow)

[RUP, RIN, DS, DRAIN] = FLUID_CAPACITY_LEGS(TERMS, GATES, STG, STAGED, ACTIVE, X, SG, R, MULT, STAGEDFLOW)

The two legs of every event under the active caps, and the waiting rooms.

Shared by the steady-state residual and the transient right-hand side so that the two solve the SAME model and not two spellings of it. What differs between them is only what a STAGED cap’s multiplier means:

stagedFlow = false (steady state) – a drain RATE theta. The room mass is then

pinned by the cap (s = inflow/theta at the fixed point), which is what makes the algebraic system square, and the proportional split across a region’s rooms falls out of the common rate.

stagedFlow = true (transient) – the admitted FLOW itself. A rate cannot start

the constrained phase at all: at the instant the region fills the room is EMPTY, so theta*s is zero however large theta is, and holding the cap needs a finite admitted flow immediately. The flow is split across the room masses, or across their inflows while the rooms are still empty – and the two rules AGREE at a fixed point, where s_j is proportional to inflow_j, so the transient and the steady state describe one model.

A held or lost cap composes as a PRODUCT of fractions either way, which is what independent blocking gives and what keeps every active cap present in the Jacobian: a product has a live derivative in each factor. Several staged caps gating one room compose HARMONICALLY, 1/theta = sum 1/theta_f, because the waits a job serves in turn add; with one cap that is exactly theta_f.

Parameters:

• terms - the event representation from FLUID_MOMENT_TERMS

• gates - which cap throttles which event, from FLUID_CAPACITY_GATES

• stg - the waiting rooms, from FLUID_CAPACITY_STAGING

• staged - CON.STAGED, true where the blocked job waits in a room

• active - indices of the caps that currently bind

• x, sg - the state and the room populations

• r - the nominal rate of every event

• mult - one multiplier per ACTIVE cap

• stagedFlow - false for the steady state, true for the transient

Returns:

rup - the rate each event FIRES at rin - the rate mass LANDS at ds - the derivative of each room drain - each room’s total outflow

See also SOLVER_FLUID_DAE, FLUID_CAPACITY_GATES, FLUID_CAPACITY_STAGING.

solver_fluid_symodes(sn, options)

SYS = SOLVER_FLUID_SYMODES(SN, OPTIONS) Build a symbolic description of the mean-field ODE system integrated by SolverFLD for the ODE-based methods. Two representations are produced, mirroring the numerical solver code paths:

form = ‘W’: dx/dt = W’*theta(x) + lambda (methods: default, matrix, pnorm)

W (nstates x nstates) phase-transition rate matrix Alambda (nstates x 1) exogenous arrival rate vector theta_s = x_s * min(n_i, S_i)/n_i (‘min’ smoothing), or

= x_s * (1+(n_i/S_i)^p_i)^(-1/p_i) (‘pnorm’ smoothing)

with theta_s = 0 for states of Source (EXT) stations, and n_i the total mass at the station i owning state s.

form = ‘J’: dx/dt = J*r(x) (methods: closing, statedep, softmin)

J (nstates x nevents) stoichiometry (jump) matrix r_e(x) = coeff(e) * factor_e(x), where factor_e depends on the

scheduling strategy of the station owning the state variable that drives event e.

Factor types (form ‘J’):

‘lin’ x_v ‘min’ x_v * min(n_i, S_i)/n_i ‘ext1’ 1 - sum of x over phases 2..end of the class at the source ‘dpsmin’ x_v * min(n_i, S_i)/ntilde_i, weight w_ir folded into coeff;

ntilde_i = sum_r w_ir n_ir (normalized DPS weights)

‘dpspw’ piecewise: x_v if n_i <= S_i, else S_i*w_ir*x_v/ntilde_i ‘fcfsw’ x_v * min(n_i, S_i)/nhat_i, phase weight folded into coeff;

nhat_i = sum_u w_u x_u with w_u = -1/D0(k,k) (mean phase time)

‘fcfsws’ x_v * softmin(n_i, S_i; alpha)/nhat_i, phase weight in coeff

Copyright (c) 2012-2026, Imperial College London All rights reserved.

fluid_visited_pairs(sn, M, K)

VISITED = FLUID_VISITED_PAIRS(SN, M, K)

Mark the (station, class) pairs the model actually routes a job into, as an (M x K) logical read off the per-chain visit ratios SN.VISITS.

A fluid result cannot decide that question from the SIZE of QN or TN. Both carry a decaying remnant of the initial state, which is spread over pairs the class never reaches, and the remnant is whatever the integrator left behind when it stopped: measured at QN = 1.3e-12 and TN = 1.3e-13 on picard05 for test_CQN_Cox_CS_7, i.e. ABOVE GlobalConstants.Zero, so a threshold on them divides one remnant by the other and reports the station’s own service time as a response time. The visit ratios come from the routing solve instead, where an unrouted pair is zero to the last bits (2.7e-17 on that pair).

SN.VISITS{c} is indexed by STATEFUL node, hence the SN.STATIONTOSTATEFUL lookup; see _kb/04-networkstruct.md. A struct carrying no visit information at all decides nothing and every pair is reported visited.

Copyright (c) 2012-2026, Imperial College London All rights reserved.

fluid_capacity_constraints(sn, terms)

CON = FLUID_CAPACITY_CONSTRAINTS(SN, TERMS)

Every capacity limit in the model as a linear constraint on the fluid state.

TWO DECLARATIONS, ONE FAMILY OF ROWS. A finite capacity region caps a SET of stations jointly (addRegion); a station cap limits one station’s own buffer (setCapacity/setClassCapacity). LINE stores the region form four ways – a region-global job cap, a per-class job cap, a memory budget with per-class sizes, and an arbitrary linear constraint pair (A,b) – and the station form two, a total and a per-class buffer. All six are the same object once written against the state:

Arow * x <= b, Arow(s) = weight of the class of coordinate s,

zero outside the stations the limit covers

so the solver carries one mechanism rather than six. The global cap is that row with every weight one, a per-class cap is the row with one class weighted, the memory budget weights each class by its size, and the linear constraint pair supplies the weights directly.

WHAT DIFFERS IS NOT THE ROW BUT WHERE THE BLOCKED JOB GOES, which is decided per admission event in FLUID_CAPACITY_GATES rather than here.

WHAT IS REFUSED, AND WHY IT IS NOT A CONSTRAINT. Only a waiting queue or a drop is a constraint on this drift. Under WAITQ a blocked job waits and is admitted later, so the population is conserved and only the admission FLOW is throttled; under DROP the flow the cap will not take is discarded, which is what LINE does with an open arrival at a full buffer. BAS/BBS/RSRD give the upstream station a blocked-server state and the retrial rules add an orbit. Each of those changes the event set itself, so each needs a different drift rather than a constraint on this one, and is refused by name.

Parameters:

• sn - NetworkStruct, carrying nregions/region/regionmembers/regionrule/ – regionweight/regionsz/regionmaxmem/regionlincon and cap/classcap/ droprule

• terms - the event representation from FLUID_MOMENT_TERMS

Returns:

con - struct with fields

A (ncon x nstate) constraint rows in state space As (ncon x nstg) the same rows over the staging coordinates,

filled in by FLUID_CAPACITY_EXTEND

b (ncon x 1) right-hand sides region (ncon x 1) which region each row came from, 0 for a

station cap

station (ncon x 1) which station, 0 for a region cap class (ncon x 1) which class it limits, 0 when several staged (ncon x 1) true where the blocked job waits in a room label (ncon x 1) cell human-readable origin, for messages member (nregions x nstate) region membership, for the rooms coordClass (1 x nstate) the class of each coordinate

empty A/b when the model has no limit

See also SOLVER_FLUID_DAE, FLUID_CAPACITY_GATES, FLUID_CAPACITY_STAGING, FLUID_CAPACITY_EXTEND, FLUID_MOMENT_TERMS.

solver_fluid_kp(sn, options)

[QN,UN,RN,TN,XVEC_IT,QNT,UNT,TNT,XVEC_T,T,ITERS,RUNTIME,QVART,SIGMAT] = SOLVER_FLUID_KP(SN, OPTIONS)

Fluid and diffusion limits of the (MAP_t/Ph_t/inf)^N network of Y. M. Ko and J. Pender, “Diffusion limits for the (MAP_t/Ph_t/inf)^N queueing network”, Oper. Res. Lett. 45 (2017) 248-253.

The mean and the covariance of the limit are integrated jointly:

dq/dt = F(t,q) = A f(t,q) dSigma/dt = J Sigma + Sigma J’ + G, J = A df/dq, G = A diag(f) A’

with A the jump matrix whose column e is the jump vector of event e and f the event rate vector. G is exactly dH dH’ of Theorem 3.3, each independent Poisson term contributing l_e l_e’ f_e. Where f is affine in q – infinite-server stations and the arrival phase process – J does not depend on q and both equations close exactly, so for the (MAP_t/Ph_t/inf)^N case the mean and the covariance are exact rather than asymptotic. Finite-server stations are admitted through the usual fluid min(x,c) term, where the covariance degrades to a linear-noise approximation and a warning is emitted.

This method does NOT reuse the closing ODE. That formulation routes a departure from the source to the destination station and returns mass through the STATIONARY arrival-instant vector pie, replacing the D1’ operator by the rank-one map pie*(D1*e)’, i.e. by the PH renewal process with representation (pie, D0). Its stationary arrival rate is exact but its autocorrelation is gone, and a non-renewal arrival stream is the entire point of a MAP.

State layout, station-major, arrival phases before service phases:

u-block one per (EXT station, class): arrival MAP phase occupancy, sum 1 x-block one per (queueing station, class): fluid count in each service phase

MAPt, PHt and NHPP all reach the integrator as a piecewise-constant (D0, D1) schedule; an NHPP is the one-phase MAPt case, lifted in LOCAL_SCHEDULE.

Copyright (c) 2012-2026, Imperial College London All rights reserved.

solver_fluid_closing(sn, options)

[QN, UN, RN, TN, XVEC, QNt, UNt, TNt, XVEC_T, T, ITERS, RUNTIME] = SOLVER_FLUID_CLOSING(QN, OPTIONS)

solver_fld_cacheqn_analyzer(sn, options)

SOLVER_FLD_CACHEQN_ANALYZER Fluid solver for integrated caching-queueing networks

Delegates the decomposition-aggregation alternation between the isolated caches and the fluid ODE solution of the surrounding queueing network to da_cacheqn, until the cache arrival rates converge. RANDOM(m) (RR) and FIFO(m) use the refined mean field (cache_miss_rmf, 1/N-accurate; Gast15 Thm 1: pi_FIFO(m)=pi_RAND(m)); strict FIFO(m) uses its own position-resolved mean field (cache_miss_sfifo_rmf). LRU/HLRU/CLIMB/QLRU have no drift-based fluid model and are rejected (the FPI characteristic-time approximation is not a fluid method). Use SolverNC/SolverMVA or SolverLDES for those.

OPTIONS.METHOD SELECTS THE QUEUEING LAYER, not the cache one. ‘rmf’ solves the surrounding network with the first-order matrix method, which is the historical behaviour; ‘minnormal’ solves it with the second-order moment closure instead, so a cache model reaches the same E[min(X,c)] treatment as any other model and returns a covariance. The cache layer is the refined mean field either way – it has no first-order alternative here – and under ‘minnormal’ its own linear noise covariance is reported alongside, which is the second moment of the item occupancy.

MOMENTRESULTS is empty unless the moment closure ran. It then carries the queueing fields of SOLVER_FLUID_MOMENTS (Sigma, QVar, QStd, sigma2, stationBlock, classBlock) plus CACHE, a struct array with one entry per cache node: NODE, PI0 (per-item miss probability), SIGMA (occupancy covariance, item-major), PI0VAR (per-item miss-indicator variance) and MISSPROBVAR (per class, the delta-method variance of the miss probability that class sees).

ode_rate_base(sn, phi, Mu, PH, M, K, enabled, q_indices, rt, Kic, ~, all_jumps)

[RATEBASE, EVENTIDX] = GETRATEBASE(sn, PHI, MU, PH, M, K, MATCH, Q_INDICES, P, KIC, SCHED_ID, ALL_JUMPS)

fluid_refine_meanfield(x, sigma2, Sigma, terms, covblk, epsrel)

[V, INFO] = FLUID_REFINE_MEANFIELD(X, SIGMA2, SIGMA, TERMS, COVBLK, EPSREL)

Refined mean field correction of a fluid fixed point (Gast, POMACS 2017).

The mean-field fixed point x* is the leading term of an expansion of the true stationary mean in powers of the system size. The next term is obtained by carrying the second moment through the drift: writing A for the Jacobian at x* and B for its Hessian tensor, the correction V solves the linear system

A*V + (1/2) * sum_{j,k} Sigma_{jk} * d2F/dx_j dx_k = 0

with Sigma the stationary covariance from FLUID_LYAPUNOV. Because Sigma scales with the population, V is the O(1/N) term of the expansion written directly in job counts, so no explicit density rescaling is needed. The Hessian contraction is evaluated without ever forming the tensor: writing Sigma = sum_m lam_m*v_m*v_m’ by eigendecomposition,

sum_{jk} Sigma_{jk} d2F/dx_j dx_k = sum_m lam_m * d2F/dv_m^2

and each directional second derivative is one central second difference, so the cost is O(rank(Sigma)) drift evaluations rather than O(n^2).

The drift must be twice differentiable for this to mean anything. The first-order closure is only piecewise linear – its second derivative is zero away from the kink and a delta at it – so this function must be called on the Gaussian-closed drift, i.e. with the SIGMA2 that SOLVER_FLUID_MOMENTS converged to under options.method=’refined’. Passing SIGMA2 = 0 is rejected rather than silently returning zero.

Parameters:

• x - fluid fixed point

• sigma2 - (M x 1)

• Sigma - (n x n)

• terms - representation from FLUID_MOMENT_TERMS

• covblk - per-station covariance blocks closing the DPS share ratio

• epsrel - relative step of the second difference (default – 1e-4)

Returns:

V - (n x 1) correction to be added to x info - struct with fields rank, stepsize, residual

See also SOLVER_FLUID_MOMENTS, FLUID_LYAPUNOV, FLUID_MIN_CLOSURE.

fluid_lyapunov(A, Qdiff, D, tol)

[SIGMA, INFO] = FLUID_LYAPUNOV(A, QDIFF, D, TOL)

Stationary covariance of the linear noise approximation.

Around a fixed point x* of the fluid drift, the fluctuation process Z = (X - x*) obeys the linear stochastic differential equation dZ = A*Z*dt + sqrt(Qdiff)*dW, whose stationary covariance solves the Lyapunov equation

A*Sigma + Sigma*A’ + Qdiff = 0, Qdiff = D*diag(r(x*))*D’

A is singular whenever the model conserves population: every closed class contributes a left null vector, so the equation has no unique solution on the full state space. It does have one on the reachable subspace, which is exactly range(D): the state can only move along jump directions, so the fluctuation lives there and nowhere else. Both A = D*diag(rateBase)*G and Qdiff map into range(D) as well, so restricting to an orthonormal basis V of range(D) is an exact reduction, not an approximation, and the reduced Lyapunov equation is nonsingular whenever the fixed point is stable.

Parameters:

• A - (n x n)

• Qdiff - (n x n) diffusion matrix D*diag(r)

• D - (n x nevents)

• tol - stability margin; eigenvalues of the reduced A with real part – above -tol are reported as non-hyperbolic (default: sqrt(eps))

Returns:

Sigma - (n x n) stationary covariance, supported on range(D) info - struct with fields rank, maxRealEig, stable

See also SOLVER_FLUID_MOMENTS, FLUID_DRIFT_JACOBIAN.

solver_fluid_tbi_iteration(sn, N, Mu, Phi, PH, P, S, xvec_it, ydefault, slowrate, Tstart, max_time, options)

[XVEC_IT, XVEC_T, T, ITER] = SOLVER_FLUID_TBI_ITERATION(SN, N, MU, PHI, PH, P, S, XVEC_IT, YDEFAULT, SLOWRATE, TSTART, MAX_TIME, OPTIONS)

Trajectory-based iteration (TBI) for the transient fluid solution. The station set is partitioned into cells (tbi_partition). On each time segment, the IVP of every cell is solved with the state of the other cells frozen at the trajectory computed in the previous iteration (Gauss-Seidel or Jacobi waveform relaxation); iterations repeat until the trajectory sup-norm gap falls below options.config.tbi_tol. Cross-cell inflows are therefore evaluated on frozen trajectories, interior flows on the live cell state, matching the decomposed ODEs of the TBI method.

Reference: Sheldon, Tuncer, Casale, “TBI: Transient Hierarchical Modeling of Large-Scale Vehicle Sharing Systems”, IEEE T-ITS.

Copyright (c) 2012-2026, Imperial College London All rights reserved.

fluid_conservation_guard(phases, njobs, chains, tol)

FCN = FLUID_CONSERVATION_GUARD(PHASES, NJOBS, CHAINS, TOL)

An ODE OutputFcn that halts an integration whose state has left the model.

WHY IT EXISTS. The moment-closure drift can leave the simplex: on a station where min(n,c) is not the identity the Gaussian correction to the per-class share can drive a coordinate negative, the drift is conservative so another grows to match, and the trajectory runs away. odeset carries NonNegative over every coordinate, so the excursion is CLAMPED rather than reported – which INJECTS mass, collapses the step size, and leaves the window never returning. That is not a slow solve, it is a solve that does not terminate: one MATLAB suite run sat in test_CQN_Cox_CS_9 for 3h16m, and the 2026-08-27 run was killed after test11_interlock_lqnx had held the suite for 100 minutes, taking every block after it down with it. See _kb/06-solver-catalog.md.

THE TEST IS AN EXACT INVARIANT, not a heuristic bound on time or magnitude. The drift conserves the population of every CLOSED CHAIN exactly, so any deviation is the clamp injecting mass and nothing else. TOL is therefore a generous fraction of that population rather than a numerical tolerance: the integrator’s own error is ~1e-4 relative, while the documented excursion reaches 5.2e4 against a true population of 0.05. A closed model whose population has moved by TOL is no longer solving the model, whatever it is converging to.

THE CHAIN IS THE CONSERVED UNIT, NOT THE CLASS, and the difference is the whole correctness of this check. SN.NJOBS(k) is the population class k STARTS with; class switching then moves jobs between the classes of one chain, so only the chain total is invariant. Watching classes instead condemns every class-switching model out of hand – measured on cqn_twoclass_hyperl (313 of 447 accepted states), on init_state_ps (286 of 310) and on every one of the 162 fluid layers an LQN builds under the srvn.cs encoding, where the chain sum never moved at all. A cache model is the same story with the hit/miss classes.

A wall-clock budget would have caught the same hang and was rejected: it makes the answer depend on how busy the host is, so the same model would fall back on one machine and not on another. This invariant is deterministic.

Parameters:

• phases - (nstations,nclasses) phases per (station,class) – builds it from Mu; fixes the coordinate blocks of the state vector

• njobs - (1,nclasses) – its population is not conserved, and it is therefore not watched

• chains - (nchains,nclasses) – one chain per class, which is the safe reading rather than a guess: with no chain map there is no class switching to merge classes, so each class IS its own conserved unit

• tol - relative deviation that trips the guard (0.1 = 10%)

Returns:

fcn - handle for odeset(‘OutputFcn’, …). Returns status 1 to halt, which

the caller reads as a short final time and turns into ‘LINE:FluidNonHyperbolic’ so the existing fallback ladder answers.

See also SOLVER_FLUID_ITERATION, SOLVER_FLUID_MOMENTS, FLUID_LYAPUNOV.

fluid_capacity_gates(sn, terms, con)

GATES = FLUID_CAPACITY_GATES(SN, TERMS, CON)

Per cap and per event: is this event an admission the cap throttles, and WHERE DOES THE STOPPED MASS GO.

An admission is an event that pushes the constrained quantity UP, read off A*D rather than off the topology so that a per-class or memory-weighted row picks out its own admissions with no extra code.

THE THREE ANSWERS ARE THE MODEL, AND THEY ARE NOT INTERCHANGEABLE. LINE’s own semantics decides which one a cap gets, and the choice is visible in the answer – a job held upstream is still counted at that station, a lost job is counted nowhere and breaks flow balance across the cap on purpose, a staged job is counted in neither and is reported as blocked mass:

STAGED a finite capacity region under a waiting queue (CON.STAGED). The job

COMPLETES upstream service and waits outside the region, which is what JMT and LDES simulate; the upstream station empties exactly as it would with no region. See FLUID_CAPACITY_STAGING.

HELD a station buffer reached by a CLOSED class. State.arrivalIsLost

refuses to lose a closed job – population conservation is a defining invariant – and returns an empty successor instead, which DISABLES the upstream departure until room frees. The job is therefore still at the upstream station, in service as far as that station’s own metrics are concerned, so the fluid analogue scales the WHOLE event: both the removal upstream and the arrival.

LOSS a station buffer reached by an OPEN class. The same predicate loses

it: the external stream is memoryless, so a job that finds the buffer full simply never enters, the arrival event still fires (ArvR counts the offered job) and only the carried flow is admitted. The fluid analogue scales the ARRIVAL leg alone and destroys the difference.

Parameters:

• sn - NetworkStruct, read for njobs (open vs closed decides held vs loss)

• terms - the event representation from FLUID_MOMENT_TERMS

• con - the constraint set from FLUID_CAPACITY_CONSTRAINTS

Returns:

gates - struct with fields

gate (ncon x nevents) does this row throttle this event held (ncon x nevents) scale the WHOLE event: the job stays upstream loss (ncon x nevents) scale the ARRIVAL leg: the rest is destroyed Dn, DnExt, Dp (nstate x nevents) the jump matrix split in three, so

that an admission can remove mass upstream at one rate and deliver it at another. DN carries the removal at real stations, DNEXT the removal from the EXT source pool – which a LOST arrival must be returned to, since that coordinate is a normalisation and its drift row a real equation – and DP the arrival. Empty where no cap exists, which keeps the uncapped drift a single matrix product.

See also SOLVER_FLUID_DAE, FLUID_CAPACITY_CONSTRAINTS.

solver_fld_cacheqn_tran(sn, options, x0cell)

SOLVER_FLD_CACHEQN_TRAN Transient refined-mean-field cache trajectory.

Converges the per-class cache arrival rates with the same decomposition- aggregation alternation as the steady solver (da_cacheqn), then integrates the mean-field drift over OPTIONS.TIMESPAN to obtain the time-resolved per-class hit/miss probabilities of each cache. This is the transient counterpart of solver_fld_cacheqn_analyzer: the steady solver drives the drift to its fixed point, whereas here the same drift is integrated over the finite window from the supplied (or default) initial occupancy.

SN: NetworkStruct with at least one Cache node. OPTIONS: solver options; OPTIONS.TIMESPAN = [t0,t1] sets the window. X0CELL: optional cell(1,ncaches); X0CELL{c} seeds cache c occupancy

(flat DDPP state vector, n_items*(h+1)); [] uses the default all-items-outside-plus-first-m-in-list initial state.

Returns the time grid TCACHE, the per-cache per-class hit/miss probability trajectories HITPROB_T/MISSPROB_T (ncaches x nclasses x nt), and the cache node indices CACHES (rows ordered as find(sn.nodetype==Cache)).

Copyright (c) 2012-2026, Imperial College London All rights reserved.

fluid_is_single_queue(sn)

FLUID_IS_SINGLE_QUEUE Check if network is a valid single-queue topology for MFQ

[isSingleQueue, fluidInfo] = FLUID_IS_SINGLE_QUEUE(sn)

Validates that the network structure is suitable for the MFQ solver. MFQ requires: Source -> single Queue -> Sink topology with single server.

Parameters:

sn (struct) – Network structure from Network.getStruct()

Returns:

isSingleQueue (logical) – True if topology is valid for MFQ fluidInfo (struct): Contains topology information and error message if invalid

solver_mfq_prio(sn, options)

SOLVER_MFQ_PRIO Single-queue open priority model via the fluid priority queue.

[QN,UN,RN,TN,xvec_it,QNt,UNt,TNt,xvec_t,t,iters,runtime] = SOLVER_MFQ_PRIO(sn, options)

Analyzes a single-queue open system with class priorities using the fluid priority queue mfq_prio_queue (G. Horvath, “Efficient analysis of the MMAP[K]/PH[K]/1 priority queue”, EJOR 246(1):128-139, 2015). The per-class Markovian arrival processes are superposed into a joint background CTMC that modulates the per-class fluid input rates; the single server drains fluid at the constant rate d = mu, serving higher-priority fluid first (preemptive).

The fluid level is interpreted as the (per-class) queue length and the fluid sojourn time as the response time, consistent with solver_mfq. The method applies when arrivals are Markov-modulated (MAP) and the service rate is class-independent; simple-exponential arrivals or class-dependent service degenerate the fluid model and fall back to the matrix fluid method.

See also: solver_mfq(), mfq_prio_queue(), fluid_is_single_queue()

solver_fluid_aoi(sn, options)

SOLVER_FLUID_AOI Age of Information analysis using Markovian Fluid Queues

[QN, UN, RN, TN, xvec_it, QNt, UNt, TNt, xvec_t, t, iters, runtime, aoiResults] = SOLVER_FLUID_AOI(sn, options)

Analyzes Age of Information for single-queue open systems using the aoi-fluid MFQ solvers.

Supported systems:

• Capacity 1 (bufferless): PH/PH/1/1 or PH/PH/1/1* (preemptive)

• Capacity 2 (single-buffer): M/PH/1/2 or M/PH/1/2* (replacement)

Parameters:

• sn (struct) – Network structure from Network.getStruct()

• options (struct) – Solver options including: - config.aoi_preemption: Override preemption/replacement probability (0-1)

Returns:

QN, UN, RN, TN – Standard performance metrics xvec_it, QNt, UNt, TNt, xvec_t, t: Transient outputs (for compatibility) iters: Number of iterations (always 1) runtime: Execution time in seconds aoiResults (struct): Age of Information results:

.AoI_mean, .AoI_var: Mean and variance of AoI .PAoI_mean, .PAoI_var: Mean and variance of Peak AoI .AoI_g, .AoI_A, .AoI_h: Matrix exponential parameters for AoI CDF .PAoI_g, .PAoI_A, .PAoI_h: Matrix exponential parameters for Peak AoI CDF .systemType: ‘bufferless’ or ‘singlebuffer’ .preemption: Preemption/replacement probability used

References

aoi-fluid toolbox by Ozancan Dogan, Nail Akar, Eray Unsal Atay BSD 2-Clause License, 2020

See also: aoi_is_aoi(), aoi_extract_params(), solveBufferless, solveSingleBuffer

ode_statedep(x, Phi, Mu, PH, M, K, enabled, q_indices, rt, Kic, nservers, w, sched_id)

RATES = ODE_RATES_STATEDEP(X, M, K, Q_INDICES, KIC, NSERVERS, W, SCHED_ID)

ode_solve_stiff(ode_h, trange, y0, ode_opt, options)

[T, YT_E] = SOLVEODESTIFF(Y0)

ode_solve(ode_h, trange, y0, ode_opt, options)

[T, YT_E] = SOLVEODE(Y0)

fluid_dae_applicable(sn, options)

[OK, REASON] = FLUID_DAE_APPLICABLE(SN, OPTIONS)

Whether SOLVER_FLUID_DAE can answer this model, used by @SolverFLD/runAnalyzer to try ‘dae’ before dropping a declined ‘minnormal’ to a first-order method.

WHY THIS EXISTS SEPARATELY FROM FLUID_MINNORMAL_APPLICABLE. The two methods state the SAME closure and differ only in how the coupled equations are discharged, so a model ‘minnormal’ accepts is almost always one ‘dae’ accepts too. Almost: ‘dae’ carries a finite-difference Jacobian over the whole unknown vector rather than one Lyapunov solve, so its state cap is lower; it closes on the per-station variance only, so DPS and GPS are out; and it has no decomposition route, so a cache model is out. Those three are exactly the difference set, and naming them here keeps the fallback ladder from entering a rung that would refuse the model a moment later.

The test is STATIC, for the same reason FLUID_MINNORMAL_APPLICABLE is: a rung chosen by trial and rollback would make the reported method depend on a failed run. The one condition that cannot be static is the NON-HYPERBOLIC fixed point the ladder exists to route around – it exists only once the mean is solved – and ‘dae’ fails on it loudly, which is what moves the ladder to its last rung.

Every condition below mirrors a refusal SOLVER_FLUID_DAE would otherwise raise, so this function and those refusals must move together.

Parameters:

• sn - NetworkStruct, after sn_nonmarkov_toph

• options - solver options

Returns:

ok - true when ‘dae’ can be selected reason - one line naming the blocking feature, empty when OK

See also SOLVER_FLUID_DAE, FLUID_MINNORMAL_APPLICABLE.

fluid_moment_terms(sn, options)

TERMS = FLUID_MOMENT_TERMS(SN, OPTIONS)

Event-based representation of the fluid population process, as required by the moment-closure methods of SolverFLD.

The closing ODEs are a density-dependent Markov population process

dx/dt = F(x) = D * r(x), r_e(x) = rateBase(e) * g_e(x)

with D the jump matrix of ODE_JUMPS_NEW and g the rate-factor vector of ODE_RATES_CLOSING. SOLVER_FLUID_ODES discards D and r once it has composed the right-hand side, but the covariance equation of the linear noise approximation and the 1/N refinement need them separately: the diffusion matrix is D*diag(r(x))*D’, which cannot be recovered from F alone. This function rebuilds that representation from SN and returns it together with drift, rate and Jacobian handles that all take the closure variance as an explicit second argument.

Parameters:

• sn - NetworkStruct

• options - solver options

Returns:

terms - struct with fields D, rateBase, eventIdx, q_indices, Kic,

enabled, w, sched, S, M, K, nstate, stationBlock, classBlock, minExact, ratesFcn, driftFcn, jacFcn

See also SOLVER_FLUID_MOMENTS, FLUID_DRIFT_JACOBIAN, ODE_JUMPS_NEW.

ode_eliminate_immediate(all_jumps, rateBase, eventIdx, sn, options)

#ok<INUSL> [ALL_JUMPS_RED, RATEBASE_RED, EVENTIDX_RED, STATE_MAP, EMAP, ABSORB] = …

ODE_ELIMINATE_IMMEDIATE(ALL_JUMPS, RATEBASE, EVENTIDX, SN, OPTIONS)

Stochastic complementation of the IMMEDIATE coordinates of the fluid ODE.

A coordinate whose exit rate is GlobalConstants.Immediate (= 1/FineTol = 1e8) is not a fast coordinate, it is an INSTANTANEOUS one: the rate is LINE’s stand-in for infinity, written by SolverLN for the branch of an activity that takes no time (an entry called with probability y < 1 carries a second PH phase at InfRate entered with probability 1-y). Integrating it numerically is meaningless work – the mode relaxes 1e8 times faster than anything else in the model, so every integrator either crawls or steps over it – and it is what made a two-station LN layer take 145 s in the JAR and 70 s in MATLAB for an answer identical, to every digit, to the one the reduced system gives in 0.35 s. See _kb/06-solver-catalog.md.

THE REDUCTION IS EXACT, not an approximation. It is the ODE twin of CTMC_STOCHCOMP: the instantaneous coordinates F are absorbed into the timed ones S by the absorption probabilities of the embedded jump chain restricted to F, so the flow that would enter F is routed straight to where F would have sent it. On the layer above, the FCFS station’s immediate phase folds into the delay – Delay(0.3767) -> Queue(0.6105) with only 0.739012 of the departures entering the queue – whose fluid limit is Q = [2.192891786795962, 1.807108213330636], the number the stiff integration spends 145 s reaching.

WHY THIS IS A STRUCTURAL COMPOSITION AND NOT A GENERATOR ROUND TRIP. Every event of ODE_JUMPS_NEW is a single -1 at EVENTIDX and a single +1 at its destination, so a path through F composes to one event, -1 at the original source and +1 at the absorbing coordinate, that keeps the original source’s GATING. Rebuilding the events from a reduced generator instead (the shape this function had before) loses that identity, and with it the event ORDER that FLUID_MOMENT_TERMS reads throughputs off – which is the whole reason the moment-closure methods used to refuse the reduction outright. EMAP carries the identity across instead: EMAP(e,o) is the expected number of times the ORIGINAL event o fires per firing of the reduced event e, so a caller maps any per-event quantity with NEWATTR = EMAP * OLDATTR and gets an exact rate accounting. It is the identity when nothing is eliminated.

A COMPOSED EVENT CAN BE A DEPARTURE AT TWO STATIONS AT ONCE, which is why a single evIsDeparture flag cannot survive the composition: a job that leaves the delay, passes through the queue’s immediate phase and returns has completed at BOTH, and both throughputs must count it. EMAP gives it a row with weight on both original events, and the null jump it composes to (-1 and +1 on the same coordinate) correctly contributes nothing to the drift and nothing to the diffusion D*diag(r)*D’.

Parameters:

• all_jumps - [nstate x nevents] jump matrix of ODE_JUMPS_NEW

• rateBase - [nevents x 1] fixed part of each event rate

• eventIdx - [nevents x 1] source coordinate of each event, which is also – the coordinate whose occupancy gates it

• sn - NetworkStruct, unused, kept for the caller’s signature

• options - solver options; options.config.immediate_tol overrides the – detection threshold

Returns:

all_jumps_red - [nstate x nevents_red] jumps, in the ORIGINAL coordinate

layout, with the eliminated rows identically zero

rateBase_red - [nevents_red x 1] reduced rates eventIdx_red - [nevents_red x 1] source coordinates, original indexing state_map - the timed coordinates that survive, in increasing order Emap - [nevents_red x nevents] expected firings of each original

event per firing of each reduced one

absorb - [nstate x nstate] projector taking an initial condition to

the reduced coordinates: identity on the timed rows, the absorption distribution on the immediate ones. Mass parked on an eliminated coordinate would otherwise be frozen there for the whole integration, because nothing moves it any more

See also CTMC_STOCHCOMP, ODE_JUMPS_NEW, ODE_RATE_BASE, FLUID_MOMENT_TERMS.

generator_to_jumps(W)

[ALL_JUMPS, RATEBASE, EVENTIDX] = GENERATOR_TO_JUMPS(W)

Convert infinitesimal generator matrix to jump/rate representation

Input:

W: [n_states x n_states] infinitesimal generator matrix

W(i,j) = rate of transition from state i to state j (i~=j) W(i,i) = -sum of outgoing rates from state i

Output:

all_jumps: [n_states x n_transitions] matrix of state change vectors rateBase: [n_transitions x 1] vector of transition rates eventIdx: [n_transitions x 1] vector of source state indices

Copyright (c) 2012-2026, Imperial College London All rights reserved.

eliminate_immediate_matrix(W, sn, options)

[W_RED, STATE_MAP] = ELIMINATE_IMMEDIATE_MATRIX(W, SN, OPTIONS)

Eliminate immediate states from infinitesimal generator matrix W

Input:

W: [n_states x n_states] infinitesimal generator matrix sn: Model structure options: Solver options

Output:

W_red: Reduced generator matrix state_map: Mapping from reduced to original state indices

Copyright (c) 2012-2026, Imperial College London All rights reserved.

fluid_pstar(method, options, M)

[USEPNORM, PSTAR] = FLUID_PSTAR(METHOD, OPTIONS, M)

@brief The p-norm smoothing exponent in force for a matrix-family fluid drift.

Ruuskanen et al., PEVA 151 (2021), eq. (26)-(27) replace the hard min(n,c) of the matrix drift by ghat = (1 + (n/c)^p)^(-1/p), so p selects WHICH DRIFT is integrated rather than tuning one.

TWO WAYS IN, AND THE METHOD NAME IS ONE OF THEM. options.pstar or options.config.pstar switch the smoothing on under any matrix-family method, which is what a caller naming an exponent means. Failing that the method ‘pnorm’ supplies its own default of 20, the exponent whose behaviour matches softmin at alpha = 20 and the one SolverFluid.java and the C++ FluidOptions already default to. Without it the name ‘pnorm’ integrated the UNSMOOTHED matrix drift, so the label named a method that never ran.

One rule, two callers: SOLVER_FLUID_MATRIX integrates the drift this selects and SOLVER_FLUID_SYMODES exports it, so the exported system and the integrated one cannot drift apart.

@param method the requested method name, qualified (‘fluid.pnorm’) or not @param options solver options, read for pstar and config.pstar @param M number of stations, the length of the returned vector @return usePnorm true when the p-norm drift is the one to build @return pstar M-by-1 exponent per station, empty when USEPNORM is false

solver_fluid_petri(sn, options)

[QN, UN, RN, TN, XVEC_IT, QNT, UNT, TNT, XVEC_T, T, ITERS, RUNTIME, MOMENTS] = SOLVER_FLUID_PETRI(SN, OPTIONS)

Fluid (mean-field) analysis of a stochastic Petri net, as the differential- algebraic system of options.method=’dae’.

A GSPN IS A DENSITY-DEPENDENT MARKOV POPULATION PROCESS, so the min-normal closure SolverFLD already carries applies to it unchanged: the marking is the population, a transition mode is a reaction, its incidence column is the jump, and lambda*min(enabling degree, servers) is the same min() the closure exists to smooth. What this analyzer adds is not a new closure but the three things a Petri net needs that a queueing network does not, and each of them is a reason the DAE form is the only one that can host it:

CONSERVATION IS A P-INVARIANT, not a chain population. A net’s conserved quantities are the left null vectors of its jump matrix, and stating them as EQUATIONS is what makes them hold to solver tolerance instead of to integrator tolerance – and what supplies the rank the drift Jacobian is missing along exactly those directions.

AN IMMEDIATE TRANSITION IS AN ALGEBRAIC FLOW. In the fluid limit a zero-time transition fires infinitely fast, so a marking that enables it cannot persist and what survives is a flow phi >= 0 pinned by the constraint that the binding input place holds no mass. An ODE has nowhere to put that; the active-set loop this method already runs for capacity caps has exactly the right shape for it.

A BOUNDED PLACE IS A LINEAR INEQUALITY, throttling the deposit leg of every firing that would overfill it, which is the loss semantics the exact engines apply (the NRM’s APPLYPLACECAPS).

THE UNKNOWNS are u = [x; s2; phi; zeta]:

x the fluid state: the marking, followed by the servers each

multi-phase mode has running in each phase (see FLUID_PETRI_TERMS)

s2 the Sigma entries the closure reads, one per TERMS.covPairs row,

eliminated by one Lyapunov solve per residual evaluation so that the full covariance never enters the Newton vector

phi one firing flow per immediate mode mu one server-latch flow per multi-phase mode, free in sign zeta one admission fraction per binding place capacity

and the stacked residual is

0 = Dneg*r + gain(zeta).*(Dpos*r) drift, one row per coordinate 0 = C*x - N conservation, one row per invariant 0 = s2 - Sigma(covPairs) closure consistency 0 = x_b (pins) immediate: the binding place is empty 0 = phi_j*w_l - phi_l*w_j (ratios) immediate: the GSPN conflict rule 0 = sum_h y(j,h) - theta_j(x) the server latch of a multi-phase mode 0 = A_c*x - b_c each binding place capacity

solved by the same damped projected Newton the queueing DAE uses.

WHAT THE SECOND MOMENT IS. The covariance is the linear noise approximation of the WHOLE state, dSigma = A Sigma + Sigma A’ + D diag(r) D’ over the STOCHASTIC event columns only – the in-flight firings of a multi-phase mode are a population like any other, and dropping them severs the only path by which the marking reaches that mode’s firing rate. An immediate flow and a server latch carry no noise of their own: both are the limit of an infinitely fast mechanism whose fluctuation is slaved, so their columns are excluded and the fluctuation is reduced onto the manifold their constraints define (LOCAL_CLAMP_TANGENT) – ORTHOGONALLY for a latch and a capacity, and OBLIQUELY for an immediate pin, along the fast column itself, so that a deposit into the pinned place is forwarded rather than deleted.

Parameters:

• sn - NetworkStruct holding Place and Transition nodes

• options - solver options; options.tol sets the Newton tolerance and – options.config.dae_maxstate caps the simultaneous solve

Returns:

as SOLVER_FLUID_DAE, with MOMENTS.petri carrying the marking covariance, the per-mode firing flows and the conservation residuals

See also SOLVER_FLUID_DAE, FLUID_PETRI_TERMS, FLUID_PETRI_IMMEDIATE, FLUID_PETRI_CONSERVATION, FLUID_PETRI_CONSTRAINTS.

fluid_capacity_staging(terms, con)

STG = FLUID_CAPACITY_STAGING(TERMS, CON)

The waiting room outside a capped region, as fluid coordinates.

WHY THE CONSTRAINT ALONE IS NOT ENOUGH, FOR A REGION. Throttling the admission events of a region does hold its population at the cap, but it holds it by slowing the UPSTREAM STATION’S COMPLETIONS – an admission event is that station finishing a job – so the blocked mass piles up at a station it has already finished being served by. Where that station is a delay the error is visible as a broken Little’s law: measured against LDES the delay reported 12 jobs at a throughput of 1.6 and a think time of 1, a factor of 7.5 out, while the region population and the throughput were both right. A waiting queue means the job COMPLETES upstream service and then waits; it is somewhere else, and the model needs somewhere else to put it.

A STATION BUFFER GETS NO ROOM, and that is not an omission. There the job genuinely does stay where it was: LINE disables the upstream departure (State.arrivalIsLost) rather than moving the job out, so the blocked mass is still at the upstream station and still counted there. Giving it a room would move mass the reference keeps in place. See FLUID_CAPACITY_GATES.

So each capped region gains one coordinate per class, and every admission into it is split in two:

upstream -> staging at the nominal rate, untouched, so the upstream

station empties exactly as it would with no region

staging -> region at theta * s, the throttled leg, with theta the

algebraic unknown the constraint pins

The two jumps sum to the original one, so nothing about the event set changes except where the mass rests in between.

ONE ROOM PER (REGION, CLASS) AND ONE MULTIPLIER PER ROW, which is what lets two caps of one region bind at once – the case a single per-region throttle had to refuse. A room gated by several active rows drains at the HARMONIC composition of their rates, 1/theta = sum 1/theta_f, because the waits a job serves in turn add; with one row that is exactly theta_f, so the single-cap answer is unchanged. GATEDBY is what records which rows gate which room.

Parameters:

• terms - event representation from FLUID_MOMENT_TERMS

• con - constraint set from FLUID_CAPACITY_CONSTRAINTS, carrying per-region – membership and which rows stage

Returns:

stg - struct with fields

n number of staging coordinates region (n x 1) region each belongs to class (n x 1) class each carries idx (nregions x K) coordinate index, 0 where there is none adm (nevents x 1) is this event an admission into a region admRegion (nevents x 1) which region, 0 if not an admission admStage (nevents x 1) which staging coordinate, 0 if not gatedBy (ncon x n) which rows gate which room

See also FLUID_CAPACITY_CONSTRAINTS, FLUID_CAPACITY_EXTEND, SOLVER_FLUID_DAE.

fluid_min_closure(n, c, s2, vc, cov_nc)

[H, DH, D2H] = FLUID_MIN_CLOSURE(N, C, S2, VC, COV_NC)

Min-normal moment closure of E[min(X,Y)] for jointly normal X, Y.

The mean-field ODEs of SolverFLD close the moment hierarchy at first order, replacing E[min(X,Y)] by min(E[X],E[Y]). That closure is exact only where min() is locally linear, so its error peaks at the kink, i.e. exactly where the two means meet – rho ~ 1 at a queueing station, the “switch point” of the process-algebra literature. This function returns instead the expectation under a bivariate normal marginal, together with its derivative with respect to the first mean.

This is the min-normal closure of Guenther, Stefanek and Bradley (EPEW/UKPEW 2012, LNCS 7587:32-47, eq. 4), implemented in their general two-population form:

E[min(X,Y)] = E[X]*Phi((E[Y]-E[X])/th) + E[Y]*Phi((E[X]-E[Y])/th)

• th*phi((E[Y]-E[X])/th)

th = (Var[X] - 2*Cov[X,Y] + Var[Y])^(1/2)

and the derivative with respect to E[X] is Phi((E[Y]-E[X])/th), the probability that X is the smaller of the two. SolverFLD only ever needs the specialisation Y = c, the deterministic server count: VC = 0 and COV_NC = 0 give th = sqrt(S2) and the expression collapses to N - (N-C)*Phi(d) - sqrt(S2)*phi(d) with d = (N-C)/sqrt(S2), the truncated-normal form. The general arguments are kept so the function IS the published closure rather than one instance of it, and so a future state-dependent capacity (a population rather than a constant) needs no new derivation.

With TH = 0 the expressions collapse to min(n,c) and to the indicator 1{n < c}, so the first-order closure is recovered exactly and callers share a single code path. The derivative at the kink is taken as 0, the right derivative of min(), which is the convention already implied by the strict inequality test in ODE_RATES_CLOSING_FACTORS.

Parameters:

• n - mean of the first argument (scalar or vector)

• c - mean of the second argument (server count when deterministic)

• s2 - variance of the first argument

• vc - variance of the second argument (default 0, deterministic)

• cov_nc - covariance of the two arguments (default 0)

Returns:

h - E[min(X,Y)] dh - d/dE[X] E[min(X,Y)] = P(X < Y) under the normal marginal d2h - d2/dE[X]^2 E[min(X,Y)] = -phi(d)/th, the density of the kink. min()

is piecewise linear, so its second derivative is carried entirely by the atom at X = Y; smoothing over the normal marginal turns that atom into the density. Zero on the degenerate branch, where the closure is the first-order one and the kink is not smoothed at all.

See also FLUID_CAPACITY_CLOSURE, SOLVER_FLUID_MOMENTS.

solver_fluid_iteration(sn, N, Mu, Phi, PH, P, S, xvec_it, ydefault, slowrate, Tstart, max_time, options)

[XVEC_IT, XVEC_T, T, ITER] = SOLVER_FLUID_ITERATION(QN, N, MU, PHI, PH, P, S, YMEAN, YDEFAULT, SLOWRATE, TSTART, MAX_TIME, OPTIONS)

fluid_resolve_default_method(sn, options, model)

[METHOD, REASON] = FLUID_RESOLVE_DEFAULT_METHOD(SN, OPTIONS, MODEL)

Which concrete method options.method=’default’ stands for, for lang=’matlab’.

Preference order: ‘rmf’ for cache models, then ‘minnormal’ whenever FLUID_MINNORMAL_APPLICABLE accepts the model, then the historical choice of ‘closing’ for DPS and ‘matrix’ otherwise. The second-order closure dominates the first-order methods on every family measured against exact CTMC and is the only method that can represent GPS at all, so it is preferred wherever it applies.

This must be called BEFORE the feature gate, not from the dispatch switch: NetworkSolver.runAnalyzerChecks validates getMethodFeatureSet(options.method), and ‘default’ is not ‘minnormal’, so a GPS model would be rejected by the gate before ever reaching the resolution. Keeping the decision in one function is what stops the gate and the dispatch from disagreeing, which would silently send a GPS model to the matrix method.

Parameters:

• sn - NetworkStruct, after sn_nonmarkov_toph

• options - solver options

• model - the Network, optional. Only the binding-capacity test needs it, – because it reads the caps the USER set rather than the derived sn.classcap (see NetworkSolver.checkBindingCapacity); omitting it skips that branch, leaving the caller with the pre-existing preference order.

Returns:

method - the resolved method name reason - why ‘minnormal’ was declined, empty when it was selected

See also FLUID_MINNORMAL_APPLICABLE, SOLVER_FLUID_MOMENTS.

solver_fluid_matrix(sn, options)

[QN,UN,RN,TN,CN,RUNTIME] = SOLVER_FLUID_MATRIX(QN, OPTIONS)

solver_fluid_odes(sn, N, Mu, phi, PH, P, nservers, sched, schedparam, options)

[ODE_H,Q_INDICES,RT_BREAKS,ABSORB] = SOLVER_FLUID_ODES(sn, N, MU, PHI, PH, P, NSERVERS, SCHED, SCHEDPARAM)

RT_BREAKS are the instants at which the time-varying rate multiplier jumps, i.e. where the returned drift is DISCONTINUOUS in t. Empty unless an NHPP schedule is configured. The caller must integrate up to each of them and restart there rather than step across; see SOLVER_FLUID_RATEMULT.

solver_fluid_initsol(sn, options)

#ok<INUSD> [INIT_SOL, STATE] = SOLVER_FLUID_INITSOL(QN, OPTIONS) %#OK<INUSD>

solver_fluid(sn, options)

[QN,XVEC_IT,QNT,UNT,XVEC_T,T,ITERS,RUNTIME] = SOLVER_FLUID(QN, OPTIONS)

ode_softmin(x, Phi, Mu, PH, M, K, enabled, q_indices, rt, Kic, nservers, w, sched_id, alpha)

RATES = ODE_SOFTMIN(X, M, K, Q_INDICES, KIC, NSERVERS, W, SCHED_ID, ALPHA)

ode_rates_closing(x, M, K, enabled, q_indices, Kic, nservers, w, sched_id, rateBase, eventIdx, sigma2, lld, covblk)

RATES = ODE_RATES_BOTTQ(X, M, K, Q_INDICES, KIC, NSERVERS, W, STRATEGY, RATEBASE, EVENTIDX, SIGMA2, LLD, COVBLK)

ode_pnorm(x, Phi, Mu, PH, M, K, enabled, q_indices, rt, Kic, nservers, w, sched_id, pstar)

DX = ODE_PNORM(X, PHI, MU, PH, M, K, ENABLED, Q_INDICES, RT, KIC, NSERVERS, W, SCHED_ID, PSTAR)

ODE derivative function using p-norm smoothing for processor-share constraint. Based on Ruuskanen et al., PEVA 151 (2021).

This provides a smoother approximation than softmin, improving ODE stability for stiff problems.

@param x State vector: per (station, class) blocks of phase populations,

laid out at the offsets given by q_indices

@param Phi Cell array Phi{i}{c}(k), completion probability of phase k of

class c at station i

@param Mu Cell array Mu{i}{c}(k), phase-k service rate of class c at station i @param PH Cell array of the phase-type service representations per station

and class, read for the entry vector of a restarted service

@param M Number of stations @param K Number of classes @param enabled (M x K) flag marking the station-class pairs that are active @param q_indices (M x K) index of the first phase of class c at station i

within x

@param rt ((M*K) x (M*K)) routing table between station-class pairs @param Kic (M x K) number of phases of class c at station i @param nservers (M x 1) server count per station, Inf for infinite server @param w (M x K) scheduling weights, used by the share disciplines @param sched_id (M x 1) numeric scheduling strategy per station @param pstar Smoothing parameter vector (one per station) or scalar

Copyright (c) 2012-2026, Imperial College London All rights reserved.

ode_jumps_new(M, K, enabled, q_indices, P, Kic)

JUMPS = ODE_JUMPS_NEW(M, K, MATCH, Q_INDICES, P, KIC, STRATEGY)

fluid_method_refusal(sn, method, options, model)

[BOOL, REASON] = FLUID_METHOD_REFUSAL(SN, METHOD, OPTIONS, MODEL)

@brief Every structural rule of a fluid method that the feature registry cannot name, asked once for the run and once for the report.

A feature set can refuse a model only for HAVING a registered feature. The rules below are of the other kinds – a station count, a class count, a server count, “requires a cache”, “requires a patience law”, a finite horizon, a fork-join route, a binding buffer – and each of them used to live only inside an analyzer arm, where a report could not see it: SolverAUTO’s findSolver called ‘fluid.mol’ runnable on a Source -> Delay -> Sink model and ‘fluid.dae’ runnable on an Erlang-fed queue, and each run then stopped.

This function delegates to the per-family predicates, so that each rule has one home and two callers: @SolverFLD/runAnalyzer stops on it, and @SolverFLD/supportsModelMethod reports it.

single-station limits FLUID_QSYS_ADMITS, then FLUID_QSYS_HORIZON for the

three that report a trajectory

‘mfq’ (butools, aoi) no refusal: off the shape FLUID_MFQ_ADMITS decides

the name RESOLVES to ‘matrix’ (SolverFLD.resolveMethod), the documented fallback the analyzer arm keeps; the predicate is read here only for the AoI exemption

‘rmf’ no refusal either: with no Cache node it resolves

to ‘matrix’, the decomposition’s own network step

‘minnormal’, ‘refined’ FLUID_MINNORMAL_APPLICABLE (the guards of

SOLVER_FLUID_MOMENTS and FLUID_MOMENT_TERMS)

‘dae’ FLUID_PETRI_APPLICABLE on a Petri net, else

FLUID_DAE_APPLICABLE

‘diffusion’ single or infinite servers only, as

SOLVER_FLUID_DIFFUSION requires

every method FLUID_FORKJOIN_ADMITS, then the binding-capacity

gate (NetworkSolver.checkBindingCapacity)

THE CAPACITY GATE HAS THREE EXEMPTIONS, each because the buffer IS the model rather than a limit the drift ignores: ‘dae’ carries it as an algebraic constraint on the drift; ‘mol’ is stated for the Mt/G/s/0 loss system, where the server count is the buffer; and the age-of-information arm of ‘mfq’ is a bufferless or single-buffer queue by definition (AOI_IS_AOI), which the gate used to refuse before the arm could run. Nothing else in the FLD tree reads sn.cap or sn.classcap.

@param sn NetworkStruct, after sn_nonmarkov_toph, as the run sees it @param method the method name, any spelling SolverFLD.canonicalMethod takes @param options solver options (horizon, closure limits, rate schedules) @param model the Network, optional; only the capacity gate reads it, because

it tests the caps the USER set rather than the derived sn.classcap

@return bool true when METHOD may run on this model @return reason the refusal, or ‘’ when BOOL is true

fluid_mvn_rectangle(m, C, a, b, npoints)

[P, LOGP] = FLUID_MVN_RECTANGLE(M, C, A, B, NPOINTS)

Rectangle probability P(A <= Y <= B) for Y ~ Normal(M, C), the multivariate normal orthant/cell integral behind @SolverFLD/getProbAggr under the moment-closure methods.

The integral has no closed form beyond one dimension, so it is evaluated by the separation-of-variables transformation of Genz (1992): the Cholesky factor of C turns the rectangle into an iterated integral over the unit cube whose integrand is a product of normal-CDF differences, and the first coordinate is integrated exactly. The remaining cube is integrated with a DETERMINISTIC Richtmyer lattice rule, frac(k*sqrt(p_j)) over the first primes, averaged with its antithetic reflection. Determinism is required here, not merely convenient: the MATLAB, Java and Python twins must return the same number, and a randomized rule would make them agree only in distribution.

C may be SINGULAR, which is the common case: a closed population fixes the sum of the station coordinates, so the covariance of a station holding a whole class is rank deficient. A coordinate whose CONDITIONAL variance vanishes is not integrated; it is a hard constraint, contributing 1 when the conditional mean falls inside its interval and 0 otherwise.

Parameters:

• m - mean vector, d-by-1

• C - covariance matrix, d-by-d, symmetric positive SEMI-definite

• a - lower corner, d-by-1, -Inf allowed

• b - upper corner, d-by-1, +Inf allowed

• npoints - lattice points per antithetic pair (default 4096)

Returns:

p - the rectangle probability in [0,1] logp - log(p), -Inf when p is zero

See also GETPROBAGGR, SOLVER_FLUID_MOMENTS, FLUID_MIN_CLOSURE.

fluid_minnormal_applicable(sn, options)

[OK, REASON] = FLUID_MINNORMAL_APPLICABLE(SN, OPTIONS)

Whether SOLVER_FLUID_MOMENTS can answer this model, used by the ‘default’ method of SolverFLD to prefer ‘minnormal’ over ‘matrix’ when it applies.

The test is STATIC: it inspects the model and the options, never the solution, so a feature the model declares is decided here and only here. The one condition that cannot be static is a NON-HYPERBOLIC fluid fixed point (balanced bottlenecks, a saturated multiclass station, an overloaded open station): it exists only once the mean is solved. FLUID_LYAPUNOV detects it and raises ‘LINE:FluidNonHyperbolic’, and @SolverFLD/runAnalyzer switches a RESOLVED ‘minnormal’ to the first-order method on that identifier alone. An explicit options.method=’minnormal’ still fails loudly, so the closure never absorbs a real defect silently.

Every condition below mirrors a guard that SOLVER_FLUID_MOMENTS, FLUID_MOMENT_TERMS or @SolverFLD/runAnalyzer would otherwise raise, so this function and those guards must move together.

Parameters:

• sn - NetworkStruct, after sn_nonmarkov_toph

• options - solver options

Returns:

ok - true when ‘minnormal’ can be selected reason - one line naming the blocking feature, empty when OK

See also SOLVER_FLUID_MOMENTS, FLUID_MOMENT_TERMS.

fluid_lld_scaling(lldrow, n)

[A, DA] = FLUID_LLD_SCALING(LLDROW, N)

Limited load-dependent rate scaling at a CONTINUOUS population.

sn.lldscaling(i,:) tabulates the station rate multiplier at integer populations 1..lldlimit, and the discrete solvers read it as lldscaling(i, min(n, lldlimit)). The fluid state is continuous, so the table is read here by linear interpolation between consecutive entries, clamped to the first entry below n=1 and to the last entry above the table end, matching the clamping the CTMC already applies.

Parameters:

• lldrow - (1 x lldlimit) – empty means no load dependence and returns a = 1, da = 0

• n - population, scalar or vector, may be non-integer or negative

Returns:

a - interpolated scaling alpha(n) da - its derivative d alpha / d n, zero on the clamped tails

See also FLUID_CAPACITY_CLOSURE, ODE_RATES_CLOSING_FACTORS.

fluid_gps_share(xk, wk, vk)

[S, DS] = FLUID_GPS_SHARE(XK, WK, VK)

Expected capacity share of a GPS station under a normal marginal.

GPS divides the server by WEIGHT among the classes that are BACKLOGGED, and then equally among that class’s own jobs (State.afterEventStation, GPS branch: cir = min(nir,1), share = w_r/(w*cir’)). The share therefore depends on the backlog INDICATOR vector, not on the populations, and that is what makes GPS unreachable for a first-order closure: with continuous x_k > 0 every class is always backlogged, the indicator is identically one, and the share collapses to the constant w_r/sum_j w_j regardless of load. That constant is the heavy-traffic limit and is wrong everywhere else, so for GPS the second moment is not a correction, it is the entire mechanism.

The closure is an EXACT enumeration rather than an expansion. The share is piecewise CONSTANT over the 2^K backlog patterns, so

E[S_r] = sum_{A ni r} P(backlog set = A) * w_r / sum_{j in A} w_j

carries no truncation error once the pattern probabilities are given. Those come from the marginals, P(N_k >= 1) = Phi((x_k - 1/2)/sigma_k) with the continuity correction for an integer population, multiplied as if the backlogs were independent. That independence is the one approximation here and it is not innocuous: in a closed network the station coordinates are NEGATIVELY correlated through population conservation, so the exact treatment would need multivariate-normal orthant probabilities (closed form to K = 3, numerical beyond).

The empty pattern contributes zero share, so sum_r E[S_r] = 1 - P(all classes empty) rather than 1. That is deliberate and is how the idle server is represented: GPS is single-server, so the backlog indicator plays the role that min(n,c) plays at a PS station, and no separate capacity term is applied by the caller.

Unlike the DPS ratio closure the expansion is NOT perturbative in sigma: as sigma_k -> 0 the probability tends to a step at x_k = 1 and its derivative phi(.)/sigma_k diverges, so the Jacobian stiffens at low variance. VK = 0 falls back to the hard indicator with zero derivative, which is the correct mean-field starting point for the outer iteration.

Parameters:

• xk - (K x 1)

• wk - (K x 1)

• vk - (K x 1)

Returns:

s - (K x 1) expected capacity shares, summing to 1 - P(station empty) ds - (K x K) Jacobian ds_r/dx_m

See also ODE_RATES_CLOSING_FACTORS, FLUID_DRIFT_JACOBIAN, FLUID_SHARE_CLOSURE.

tbi_partition(sn, options)

CELLS = TBI_PARTITION(SN, OPTIONS)

Partition the station set into cells for trajectory-based iteration. Honors options.config.tbi_cells, a cell array of disjoint station index vectors covering 1:nstations. Otherwise stations are agglomerated greedily on the symmetrized station-level routing weights, targeting options.config.tbi_cellsize stations per cell (default 5), so that strongly coupled stations share a cell and connecting flows stay weak.

Copyright (c) 2012-2026, Imperial College London All rights reserved.

solver_fluid_ratemult(numEvents, M, K, enabled, q_indices, Kic, Mu, eventIdx, options)

[TGRID, MMAT, TBREAKS] = SOLVER_FLUID_RATEMULT(NUMEVENTS, M, K, ENABLED, Q_INDICES, KIC, MU, EVENTIDX, OPTIONS)

Build the time-varying per-event rate multiplier for the closing fluid ODE. Returns TGRID (1 x ngrid, strictly increasing) and MMAT (numEvents x ngrid); event e is scaled at time t by fluid_interpcols(TGRID, MMAT, t)(e). Returns [] when no time-varying source is configured, so the caller keeps the legacy autonomous closure.

TBREAKS lists the instants at which the multiplier JUMPS, i.e. the NHPP schedule’s own segment bounds. They are reported separately from TGRID because they are the only grid instants an integrator must not step across: the schedule is piecewise CONSTANT, so the drift is discontinuous there and nothing in the right-hand side tells a step controller where the jump is. SOLVER_FLUID_ITERATION makes each one an integration boundary. The other two sources carry no breaks: RATE_TRAJ and RATE_SCHED are sampled trajectories meant to be read as piecewise linear, which every integrator handles.

Two independent, composable sources are honoured (both reduce to a per-event multiplicative factor because the closing rate is rate = rateBase .* theta(x) and rateBase is linear in the station-class service/arrival rate):

1. options.config.rate_traj = {tgrid, Mmat} - a caller-supplied event multiplier matrix (used by the coupled LN layer transient).

2. options.config.nhpp_sched - a struct array of non-homogeneous (NHPP) source intensities, each with fields:

   .station station index in the sn station space (must be EXT/source) .class class index .nhpp the process handle exposing getRateAt(t) and, for grid

   construction, getBreakpoints/getPeriod/isCyclic

   The nominal (time-average) rate baked into rateBase for that station-class is Mu{station}{class}(1); the multiplier is getRateAt(t)/nominal, applied to every event sourced at that station and class (eventIdx == q_indices(station,class) + kic - 1).

Copyright (c) 2012-2026, Imperial College London All rights reserved.

fluid_hide_immediate(sn, options)

TF = FLUID_HIDE_IMMEDIATE(SN, OPTIONS)

Whether the fluid drift of this model should be built on the stochastic complement of its INSTANTANEOUS coordinates (see ODE_ELIMINATE_IMMEDIATE).

Every fluid route that builds its drift from the station/class/phase event set asks here rather than reading OPTIONS.CONFIG.HIDE_IMMEDIATE directly, so that the answer is the same one across ‘matrix’, ‘closing’, ‘statedep’, ‘tbi’, ‘minnormal’, ‘refined’ and ‘dae’. The flag defaults to TRUE for SolverFLD: a coordinate whose exit rate is GlobalConstants.Immediate is LINE’s stand-in for infinity, and integrating it is meaningless work no integrator does well.

THE STOCHASTIC PETRI NET ROUTE IS THE ONE EXCEPTION, and it is not a refusal. SOLVER_FLUID_PETRI carries immediate firings as ALGEBRAIC unknowns of an index-1 DAE (FLUID_PETRI_IMMEDIATE), which is a stronger treatment than absorbing them: it keeps the firing flow itself as a solved quantity rather than folding it into the timed events. That route never builds the event set this reduction acts on, so the answer here is simply FALSE and the flag is left alone rather than being turned into an error the way it once was.

Parameters:

• sn - NetworkStruct

• options - solver options

Returns:

tf - true when the immediate coordinates are to be complemented away

See also ODE_ELIMINATE_IMMEDIATE, SOLVER_FLUID_ODES, FLUID_MOMENT_TERMS.

fluid_capacity_extend(con, stg, terms)

CON = FLUID_CAPACITY_EXTEND(CON, STG, TERMS)

Extend every cap to the staging coordinates that hold mass INSIDE it.

A waiting room is outside the region it feeds, which is the whole point of it – but it is not outside every OTHER limit. Where two regions overlap, an admission into the inner one is an INTERNAL move of the outer one: the job leaves a station of the outer region, waits, and re-enters a station of the same outer region, never having left it. Counting only the state coordinates would take that mass out of the outer cap for as long as it waits, so the outer cap would be met on paper while the region actually held more; and the Newton system that results is inconsistent rather than merely inexact – two overlapping regions stalled at residual 5e-1 with the inner cap exceeded.

A room counts toward a row when the row weighs the room’s DESTINATION and also weighs every station that FEEDS it – that is exactly “the job was inside and stays inside”. A room fed from outside is a queue at the door and counts nowhere, as before.

Parameters:

• con - constraint set from FLUID_CAPACITY_CONSTRAINTS

• stg - staging from FLUID_CAPACITY_STAGING

• terms - the event representation from FLUID_MOMENT_TERMS

Returns:

con - the same struct with As (ncon x stg.n) filled in

See also FLUID_CAPACITY_CONSTRAINTS, FLUID_CAPACITY_STAGING.

fluid_drift_jacobian(x, M, K, enabled, q_indices, Kic, nservers, w, sched_id, rateBase, eventIdx, D, sigma2, lld, covblk)

[A, G] = FLUID_DRIFT_JACOBIAN(X, M, K, ENABLED, Q_INDICES, KIC, NSERVERS, W, SCHED_ID, RATEBASE, EVENTIDX, D, SIGMA2, LLD, COVBLK)

Analytic Jacobian of the fluid drift F(x) = D*(rateBase .* g(x)(eventIdx)), where g is the rate-factor vector built by ODE_RATES_CLOSING. The Jacobian drives the covariance (Lyapunov) equation of the linear noise approximation and the 1/N refinement, so it must mirror ODE_RATES_CLOSING branch by branch: any policy without a case there keeps g = x and contributes the identity here.

With sigma2 = 0 the derivative of the occupancy factor is the indicator of the unsaturated region, i.e. the a.e. derivative of the first-order closure. With sigma2 > 0 it is the smooth derivative of the Gaussian closure returned by FLUID_MIN_CLOSURE.

Parameters:

• x - state vector (phase-resolved populations)

• D - (n x nevents)

• sigma2 - (M x 1)

• remaining arguments as in ODE_RATES_CLOSING

Returns:

A - (n x n) Jacobian dF/dx G - (n x n) Jacobian dg/dx of the rate factors

See also ODE_RATES_CLOSING, FLUID_MIN_CLOSURE, FLUID_LYAPUNOV.

fluid_interpcols(tg, B, tt)

B = FLUID_INTERPCOLS(TG, B, TT) Clamped piecewise-linear interpolation of the columns of B at scalar time tt.

B is (nrows x ngrid) with column j sampled at time tg(j); tg is a strictly increasing grid. The result b is (nrows x 1). Times outside [tg(1),tg(end)] are clamped to the boundary columns (zero-order hold outside the grid). This is the shared time-varying-input evaluator used both by the TBI transient (cross-cell inflow drift) and by the time-varying rate-multiplier closure in solver_fluid_odes.

Copyright (c) 2012-2026, Imperial College London All rights reserved.

solver_mfq_aoi(sn, options)

SOLVER_MFQ_AOI Age of Information analysis using Markovian Fluid Queues

[QN, UN, RN, TN, xvec_it, QNt, UNt, TNt, xvec_t, t, iters, runtime, aoiResults] = SOLVER_MFQ_AOI(sn, options)

Analyzes Age of Information for single-queue open systems using the aoi-fluid MFQ solvers.

Supported systems:

• Capacity 1 (bufferless): PH/PH/1/1 or PH/PH/1/1* (preemptive)

• Capacity 2 (single-buffer): M/PH/1/2 or M/PH/1/2* (replacement)

Parameters:

• sn (struct) – Network structure from Network.getStruct()

• options (struct) – Solver options including: - config.aoi_preemption: Override preemption/replacement probability (0-1)

Returns:

QN, UN, RN, TN – Standard performance metrics xvec_it, QNt, UNt, TNt, xvec_t, t: Transient outputs (for compatibility) iters: Number of iterations (always 1) runtime: Execution time in seconds aoiResults (struct): Age of Information results:

.AoI_mean, .AoI_var: Mean and variance of AoI .PAoI_mean, .PAoI_var: Mean and variance of Peak AoI .AoI_g, .AoI_A, .AoI_h: Matrix exponential parameters for AoI CDF .PAoI_g, .PAoI_A, .PAoI_h: Matrix exponential parameters for Peak AoI CDF .systemType: ‘bufferless’ or ‘singlebuffer’ .preemption: Preemption/replacement probability used

References

aoi-fluid toolbox by Ozancan Dogan, Nail Akar, Eray Unsal Atay BSD 2-Clause License, 2020

See also: aoi_is_aoi(), aoi_extract_params(), solveBufferless, solveSingleBuffer

fluid_qsys_admits(sn, method)

[BOOL, REASON] = FLUID_QSYS_ADMITS(SN, METHOD)

@brief The shape a single-station fluid limit is stated for, asked as a gate.

‘ggisgi.fluid’, ‘ggingi.tga’, ‘tvms’, ‘mtginf’ and ‘mol’ are closed forms for ONE open class through ONE Source and ONE queueing station, not integrations of the network drift, so they cannot answer any other shape. Three of them are limits for a queue customers ABANDON and need a reneging patience law; three need a finite server count (the Erlang formula of ‘mol’, the sqrt(n) fluctuation of ‘tga’, the staffing of ‘tvms’). Every rule here used to be an inline error inside SOLVER_FLUID_QSYS_ANALYZER, invisible to a caller: a report offered ‘mol’ on a Source -> Delay -> Sink model and the run stopped on the infinite server count.

A time-varying SERVICE law is refused too: the analyzer reads the service ccdf off a stationary MAP pair and a schedule slot (NHPP/MAPt/PHt) is not one, so it silently fell back to an exponential of the same mean. The ARRIVAL side is where a schedule belongs, and SN_ARRIVAL_RATE_FUN reads it there.

The registry cannot name any of this (“exactly one station”, “requires a patience law”, “a finite server count”), hence the structural predicate. Called by SOLVER_FLUID_QSYS_ANALYZER, so the run stops on it, and by FLUID_METHOD_REFUSAL, so a report sees the same verdict. One predicate, two callers. The horizon rule of the three time-varying limits is a separate predicate, FLUID_QSYS_HORIZON, because it is a rule on the OPTIONS.

@param sn NetworkStruct of the model @param method the method name, any spelling SolverFLD.canonicalMethod takes @return bool true when METHOD may run on this model @return reason the refusal, or ‘’ when BOOL is true

fluid_dae_project(u, lb)

U = FLUID_DAE_PROJECT(U, LB)

An iterate projected onto its feasible box, the lower bound only.

WHY PROJECTED, AND NOT MERELY CLAMPED INSIDE THE RESIDUAL. The bounded unknowns of the DAE – a variance, an admission throttle, an immediate firing flow – are read through max(0,.) by the residual that uses them. Clamping there ALONE is a trap: once an iterate goes negative the residual stops depending on it, so the finite-difference column is exactly zero, the solver has no derivative to climb back on, and the unknown is pinned at the boundary for good. That is not hypothetical – it pinned the throttle at zero in every configuration of the capacity sweep, and the region then settled BELOW its cap instead of on it, while the residual reported the constraint as the only unmet equation.

Projecting the ITERATE instead keeps every evaluation inside the feasible box, where the forward difference across max(0,.) is live even exactly at zero, so the boundary can be left again.

Parameters:

• u - the iterate

• lb - lower bound per unknown, one entry per unknown; -Inf for a free one. – It is a VECTOR, never a count. The earlier form also accepted a scalar NFREE meaning “the first NFREE unknowns are free”, which is ambiguous the moment U holds a single unknown: a one-place, one-class net with no closure pair produced lb = -Inf, that was read as NFREE = -Inf, and the projection indexed U(-Inf+1:end).

Returns:

u - the projected iterate

See also FLUID_DAE_NEWTON, SOLVER_FLUID_DAE, SOLVER_FLUID_PETRI.

solver_fluid_moments(sn, options)

[QN, UN, RN, TN, XVEC_IT, QNT, UNT, TNT, XVEC_T, T, ITERS, RUNTIME, MOMENTS] = SOLVER_FLUID_MOMENTS(SN, OPTIONS)

Second-order moment-closure fluid analysis, backing options.method ‘minnormal’ and ‘refined’.

The default fluid methods close the moment hierarchy at first order: the drift of the mean depends on E[min(X_i,c_i)], which they replace by min(E[X_i],c_i). No second moment ever enters, so no variance is produced and the mean itself is biased wherever min() is not locally linear. The two methods here reinstate the second moment:

‘minnormal’ Min-normal closure (Guenther, Stefanek, Bradley). The drift

uses E[min(X_i,c_i)] under a normal marginal whose variance is produced by the covariance equation, so mean and covariance are solved self-consistently by fixed-point iteration. This corrects the mean, most visibly near rho = 1 where the first-order closure is worst.

‘refined’ Refined mean field (Gast). Adds the O(1/N) correction term of

the mean-field expansion to the MEAN-FIELD fixed point (not to the ‘minnormal’ one, which already resums it), computed by FLUID_REFINE_MEANFIELD.

All performance measures are read back from the same event representation that defines the drift, so throughputs balance flow at the fixed point under whichever closure was used.

Parameters:

• sn - NetworkStruct

• options - solver options; options.method selects the closure

Returns:

as SOLVER_FLUID_CLOSING, plus MOMENTS with fields Sigma (state-level covariance), QVar (station-class queue-length variance), QStd, sigma2 (per-station population variance), refinement (1/N correction, ‘refined’ only), outerIters

See also FLUID_MOMENT_TERMS, FLUID_LYAPUNOV, FLUID_REFINE_MEANFIELD.

solver_mfq(sn, options)

SOLVER_MFQ Solves single-queue open systems using BUTools MAPMAP1

[QN,UN,RN,TN,xvec_it,QNt,UNt,TNt,xvec_t,t,iters,runtime] = SOLVER_MFQ(sn, options)

Uses the BUTools library’s FluFluQueue function to analyze MAP/MAP/1 queues. Provides exact steady-state queue length and sojourn time moments for single-queue open systems with phase-type arrivals and service.

This solver applies only to single-queue topologies:

Source (external arrivals) -> Queue (finite server) -> Sink

Applicability:

• Single-queue open system

• Single-server (c=1) or infinite-server (c=Inf)

• No feedback loops

• Supported distributions: Exp, Erlang, HyperExp, Cox, APH, MAP, PH, MMDP

Parameters:

• sn (struct) – Network structure from Network.getStruct()

• options (struct) – Solver options including: - method: ‘mfq’ (should be set by dispatcher) - tol: Numerical tolerance (default: 1e-14)

Returns:

QN (MxK) – Mean queue lengths at each station UN (MxK): Utilizations at each station RN (MxK): Mean response times at each station TN (MxK): Throughputs at each station xvec_it (cell): Final state vector for compatibility QNt (MxK cell): Transient queue lengths (empty for MAPMAP1) UNt (MxK cell): Transient utilizations (empty for MAPMAP1) TNt (MxK cell): Transient throughputs (empty for MAPMAP1) xvec_t: Transient state vectors (empty for MAPMAP1) t: Time vector iters: Number of iterations (always 1 for MAPMAP1) runtime: Execution time in seconds

References

BUTools: Queueing and traffic modeling library https://github.com/ghorvath78/butools

See also: solver_fluid_analyzer(), fluid_is_single_queue(), FluFluQueue

fluid_outputfcn_chain(fcns)

FCN = FLUID_OUTPUTFCN_CHAIN(FCNS)

Compose several ODE OutputFcns into the ONE slot odeset provides.

WHY IT EXISTS. odeset carries a single ‘OutputFcn’, and the fluid solver has two independent reasons to halt a window: FLUID_CONSERVATION_GUARD, when the moment-closure drift has left the model, and FLUID_FIXED_POINT_GUARD, when the state has settled and the rest of the span is known in closed form. Neither is a special case of the other and both must run on the SAME integrator, so the slot holds a chain rather than a choice. Empty entries are dropped, so a caller may pass a slot that was never filled.

MATLAB’s protocol is honoured for every member: ‘init’ first, ‘’ after each accepted step, ‘done’ last. A halt is the disjunction – ANY member returning nonzero halts – but every member is still called on that step, because a member may be keeping state and skipping it would corrupt it. ‘done’ always reports 0: it is a teardown notification, not a decision point, and MATLAB ignores the status there.

Parameters:

fcns - cell array of OutputFcn handles, called as fcn (t, y, flag) – elements are ignored, so {[], guard} is a chain of one

Returns:

fcn - handle for odeset(‘OutputFcn’, …), or [] when nothing is left after

dropping the empties, which the caller may store back into the slot unchanged

See also SOLVER_FLUID_ITERATION, FLUID_CONSERVATION_GUARD, FLUID_FIXED_POINT_GUARD.

solver_fluid_qsys_analyzer(sn, options)

[Q,U,R,T,C,X,QT,UT,TT,RUNTIME,ACTUALMETHOD] = SOLVER_FLUID_QSYS_ANALYZER(SN, OPTIONS)

The single-station fluid limits: a Source -> Queue -> Sink model with one class, answered by a closed-form fluid or Gaussian limit rather than by integrating the network drift.

WHY THESE ARE FLUID METHODS AND NOT MVA ONES. Each depends on the service or patience law BEYOND ITS MEAN – the stationary point of the Liu-Whitt model is where the patience ccdf crosses 1/rho, the Mt/G/Inf mean is a convolution with the service ccdf – and each is the limit of a sequence of systems, not an approximation to a fixed one. That is the fluid solver’s contract.

METHODS

‘ggisgi.fluid’ - stationary point of the G/GI/s+GI fluid model

(Liu and Whitt, Operations Research 60(5), 2012)

‘ggingi.tga’ - truncated Gaussian approximation, which adds the

O(sqrt(n)) fluctuation around that point (Liu, Whitt and Yu, Naval Research Logistics 63(3), 2016)

‘tvms’ - the Gt/Mt/st+GI many-server fluid queue at CONSTANT

staffing (Liu and Whitt, INFORMS J. Computing 26(1), 2014)

‘mtginf’ - the exact Mt/G/Inf mean (Eick, Massey and Whitt,

Management Science 39(2), 1993)

‘mol’ - the modified-offered-load approximation for the finite

server count (Massey and Whitt, Ann. Appl. Prob. 4(4), 1994)

The transient tables Qt, Ut, Tt are returned for the three time-varying methods and are empty for the two stationary ones.

See also SOLVER_FLUID_ANALYZER, SN_ARRIVAL_RATE_FUN, SN_PATIENCE_HANDLES.

solver_fluid_mfq(sn, options)

SOLVER_FLUID_MFQ Solves single-queue open systems using BUTools MAPMAP1

[QN,UN,RN,TN,xvec_it,QNt,UNt,TNt,xvec_t,t,iters,runtime] = SOLVER_FLUID_MFQ(sn, options)

Uses the BUTools library’s FluFluQueue function to analyze MAP/MAP/1 queues. Provides exact steady-state queue length and sojourn time moments for single-queue open systems with phase-type arrivals and service.

This solver applies only to single-queue topologies:

Source (external arrivals) -> Queue (finite server) -> Sink

Applicability:

• Single-queue open system

• Single-server (c=1) or infinite-server (c=Inf)

• No feedback loops

• Supported distributions: Exp, Erlang, HyperExp, Cox, APH, MAP, PH, MMDP

Parameters:

• sn (struct) – Network structure from Network.getStruct()

• options (struct) – Solver options including: - method: ‘mfq’ (should be set by dispatcher) - tol: Numerical tolerance (default: 1e-14)

Returns:

QN (MxK) – Mean queue lengths at each station UN (MxK): Utilizations at each station RN (MxK): Mean response times at each station TN (MxK): Throughputs at each station xvec_it (cell): Final state vector for compatibility QNt (MxK cell): Transient queue lengths (empty for MAPMAP1) UNt (MxK cell): Transient utilizations (empty for MAPMAP1) TNt (MxK cell): Transient throughputs (empty for MAPMAP1) xvec_t: Transient state vectors (empty for MAPMAP1) t: Time vector iters: Number of iterations (always 1 for MAPMAP1) runtime: Execution time in seconds

References

BUTools: Queueing and traffic modeling library https://github.com/ghorvath78/butools

See also: solver_fluid_analyzer(), fluid_is_single_queue(), FluFluQueue

fluid_mfq_admits(sn)

[BOOL, REASON, ISAOI] = FLUID_MFQ_ADMITS(SN)

@brief Can the ‘mfq’ method (and its aliases ‘butools’, ‘aoi’) run here?

The Markovian fluid queue answers ONE open station, Source -> Queue -> Sink. SOLVER_FLUID_ANALYZER used to fall back to the matrix method on any other shape, and SOLVER_MFQ_PRIO did the same inside its own branch, so a caller’s ‘mfq’ label was answered by a different algorithm under that label. A listed name must run as itself, so the shapes are refused by name instead.

Three arms, tried in the analyzer’s own order:

the age-of-information shape (AOI_IS_AOI): one open class, one server, a

buffer of 1 or 2, FCFS/LCFS/LCFSPR. The finite buffer IS the model there, so ISAOI lets a caller exempt it from the binding-capacity gate, as ‘mol’ is exempt for the same reason.

the priority branch (SOLVER_MFQ_PRIO), taken when the classes carry

distinct priorities: a class-independent service rate, a MAP {D0,D1} arrival for every open class, at least one of them modulated. Each of those used to be a silent fallback inside the branch.

the plain branch (SOLVER_MFQ), which reads the arrival and service

processes of class 1 only, so it is stated for ONE open class.

Closed classes are refused by the feature set (SolverFLD.getMethodFeatureSet) rather than here: SOLVER_MFQ skips them and reports zeros.

Called by SOLVER_FLUID_ANALYZER, so the run stops on it, and by FLUID_METHOD_REFUSAL, so a report sees the same verdict before running. One predicate, two callers.

@param sn NetworkStruct of the model @return bool true when ‘mfq’ may run on this model @return reason the refusal, or ‘’ when BOOL is true @return isAoI true when the age-of-information arm is the one that runs

fluid_dae_newton(resid, u, tol, maxit, lb)

[U, IT, CONVERGED, RESNORM] = FLUID_DAE_NEWTON(RESID, U, TOL, MAXIT, LB)

Damped PROJECTED Newton with a finite-difference Jacobian and an Armijo backtrack on the residual norm, the solver behind the DAE formulation of the fluid closures.

The step is solved in least squares rather than by a square factorisation: the drift block is rank deficient by exactly the number of conserved quantities, and the conservation rows restore that rank, so the stacked system is consistent and overdetermined rather than square. The same holds for the Petri route, where the conserved quantities are the net’s P-invariants.

Parameters:

• resid - handle resid (u, quiet) – it cannot be evaluated at that iterate so the line search backs off

• u - starting iterate

• tol - convergence tolerance on the infinity norm of the residual

• maxit - maximum Newton steps

• lb - lower bound per unknown (-Inf where free) – the layout “the first NFREE unknowns are free, the rest are non-negative”; see FLUID_DAE_PROJECT

Returns:

u - the last iterate it - Newton steps taken converged - whether RESNORM fell below TOL resnorm - infinity norm of the residual at U

See also FLUID_DAE_FDJAC, FLUID_DAE_PROJECT, SOLVER_FLUID_DAE, SOLVER_FLUID_PETRI.

solver_fluid_diffusion(sn, options)

SOLVER_FLUID_DIFFUSION Diffusion approximation for closed multiclass BCMP networks

[QN,UN,RN,TN,XVEC_IT,QNT,UNT,TNT,XVEC_T,T,ITERS,RUNTIME] = SOLVER_FLUID_DIFFUSION(SN, OPTIONS)

Implements diffusion approximation using the Euler-Maruyama method for SDEs. This method models queue lengths as continuous variables with Brownian noise, providing stochastic extensions to deterministic fluid approximations.

Supported models:

• Closed multiclass queueing networks only

• Scheduling: PS, FCFS, INF (delay stations)

• Single-server or infinite-server stations

Options:

• iter_max: Number of simulation steps (default: 10000)

• timestep: Time step for Euler-Maruyama integration (default: 0.01)

Copyright (c) 2012-2026, Imperial College London All rights reserved.

solver_fluid_passage_time(sn, options)

RTRET = SOLVER_FLUID_PASSAGE_TIME(QN, OPTIONS)

solver_fluid_analyzer_jline(network, options)

[QN, UN, RN, TN, CN, XN, T, QNT, UNT, TNT, XVEC] = SOLVER_FLUID_ANALYZER_JLINE(QN, OPTIONS)

fluid_forkjoin_admits(sn, method)

[BOOL, REASON] = FLUID_FORKJOIN_ADMITS(SN, METHOD)

@brief Can METHOD run the fluid fork-join fixed point on this model?

A fork-join model is not integrated as one drift: the MMT transform replaces the fork by auxiliary OPEN classes arriving at a Source it adds (closed model or not, see ModelAdapter.mmt) and the answer is the fixed point of solving that mixed model repeatedly (@NetworkSolver/fjFixedPoint). So the inner solve is always a mixed network, and a method refused on one is refused here by name rather than by a failure inside the fixed point: ‘softmin’ and ‘statedep’ have no EXT branch, ‘refined’ is closed-only, and on an OPEN model the DAE form has no unknowns for the auxiliary classes, so its inner solve fails on the class count rather than returning a drift.

Called by @SolverFLD/runAnalyzer, so the run stops on it, and by @SolverFLD/supportsModelMethod, so a CALLER sees the same verdict before paying for the fixed point. One predicate, two callers.

@param sn NetworkStruct of the model @param method the concrete method name @return bool true when METHOD may run the fork-join fixed point here @return reason the refusal, or ‘’ when BOOL is true

fluid_dae_fdjac(resid, u, G)

J = FLUID_DAE_FDJAC(RESID, U, G)

Finite-difference Jacobian of a DAE residual at U, given its value G there.

Forward differences by default, and a BACKWARD difference for any column whose forward step lands where the residual cannot be evaluated – the closure raises on a non-hyperbolic fixed point, so a step that crosses into that region returns empty rather than a number. A column that fails both ways is left at zero and the least-squares step absorbs it.

Parameters:

• resid - handle resid (u, quiet) – be evaluated

• u - the iterate

• G - the residual at U

Returns:

J - (numel(G) x numel(u)) Jacobian

See also FLUID_DAE_NEWTON, SOLVER_FLUID_DAE, SOLVER_FLUID_PETRI.

ode_rates_closing_factors(x, M, K, enabled, q_indices, Kic, nservers, w, sched_id, sigma2, lld, covblk)

RATES = ODE_RATES_CLOSING_FACTORS(…) per-coordinate service share, before event indexing and before the constant rate factors are applied. Kept separate so the moment-closure methods can read the same service shares the ODE integrated (see SOLVER_FLUID_MOMENTS).

SIGMA2 is the per-station closure variance (0 or omitted: first-order closure). LLD is sn.lldscaling (empty: no load dependence). COVBLK is a per-station cell of coordinate covariance blocks, which closes the DPS capacity-share ratio at second order (see FLUID_SHARE_CLOSURE). All three leave the legacy code path bit-identical when absent, so the untouched methods are unaffected.