api.spn

spn_lpbnd(sn, options)

BND = SPN_LPBND(SN) BND = SPN_LPBND(SN, OPTIONS) Linear-programming bounds on the mean marking and the throughputs of a stochastic timed Petri net.

The stationary chain is relaxed to a MOMENT POLYTOPE: the uniformized evolution equation is written for E[X_p], E[X_p^2] and E[X_p1 X_p2], which gives linear equalities among the mean marking x, the enabling probabilities q and the products y(p,t) = E[X_p e_t]; behavioural and probabilistic inequalities are added on top; and every reported measure is then obtained by minimising and maximising its linear form over that polytope. Any stationary point of the true chain satisfies every row, so the two optima BRACKET the exact value whatever the polytope leaves out.

This is the Petri-net sibling of the QRF bounds in SolverBA: same technique, a different index space, and a LINEAR objective, so there is no stationary point to escape from and the answer is a property of the model alone.

VARIABLES, over place levels l = 1..L and modes e = 1..E:

x(l) E[X_l], mean tokens >= 0 q(e) P(mode e enabled) in [0,1] th(e) throughput of mode e >= 0 u(e) state-equation firing counts >= 0 y(l,e) E[X_l e_e] >= 0 (Markovian only)

u is EXISTENTIAL and is not reported: E[X] is a convex combination of reachable markings, each of which is m0 + C h for some nonnegative integer h, so the mean satisfies m0 + C u for some nonnegative real u. It is not a mean firing count and has no steady-state value.

LEVELS ARE (place, class) PAIRS, PLACE-MAJOR, level (pp-1)*R + k, the same coordinates SPN_MDD, SPN_SINVARIANTS and SPN_CONV use. MODES are (transition, mode) pairs in node order. Both index spaces are returned in BND.

THE TOKEN COUNTS CREATED BY A FIRING ARE DETERMINISTIC IN LINE, which removes a whole branch of the reference: it allows sigma_{t,p}(n) to be random and splits the covariance family into an independent case (its eq. 7) and a selective one (its eq. 8). SETFIRINGOUTCOME takes an integer weight, so E[sigma^2] = sigma^2 and E[sigma_p1 sigma_p2] = sigma_p1 sigma_p2 hold exactly and eq. (7) is the correct form. Eq. (8) has no LINE model behind it and is deliberately absent.

– Input SN : a NetworkStruct holding Places and Transitions OPTIONS : struct, all fields optional

markovian true (default) uses the second-moment, covariance and Little’s

law families, which need exponential firing times; false drops them and the whole y block, leaving the operational bound, which needs only a mean firing time and so admits any phase-type law. A law that is not phase-type at all is refused by both: SN carries such a mode as its own parameter list rather than a (D0,D1) pair, and no mean can be read off it without a per-distribution table

assumelive false (default). true adds the two liveness rows, which are

valid only on a live net; see the note below

init initial tokens per place level, place-major; omitted reads the

declared marking through State.initialOccupancy

tol LP feasibility slack added to the inequality sides, default 0 verbose print the polytope size

– Output BND : struct with fields

places place node indices, in level order levelname 1 x L cell, “Place.Class” modes 1 x E struct array with fields trans, mode, enab, inhib,

fire, rate

tokens 2 x L, row 1 the minimum and row 2 the maximum of x(l) placeTput 2 x L, bracket on sum_e pi_e(l) mu_e q(e), tokens removed

from level l per unit time

modeTput 2 x E, bracket on th(e) modeUtil 2 x E, bracket on q(e) bound 1 x L, the a priori per-level bound B_l read off the

P-invariants, Inf where none constrains the level

nplacelevels, nclasses, markovian, nvars, nrows

LIVENESS IS OFF BY DEFAULT, AND THAT IS DELIBERATE. The reference’s two liveness rows (sum_t q_t >= 1 and x_p <= sum_t y_{p,t}) hold only on a live net, and liveness is not something this function can cheaply certify – an inhibitor arc alone is enough to deadlock a net that looks well formed. A bound that silently assumed it would be wrong rather than loose on exactly the models where a bound is most wanted, so the rows are available under OPTIONS.assumelive and absent otherwise.

WHAT THE ROWS ARE WORTH, MEASURED. They are the whole of the lower side. On the reference’s own Table 2 (its Fig. 2b production line, five rate vectors) ASSUMELIVE reproduces its published l.b. column to four decimals – 1.1653 against 1.165, 1.8288 against 1.829, 1.5814 against 1.581, 1.3592 against 1.359, 1.3497 against 1.350 – while without them the Markovian lower bound collapses onto the OPERATIONAL one (0.9302, 1.4815, 1.1111, 1.1110 against that column’s 0.930, 1.481, 1.111, 1.111) on four of the five. The upper side needs neither row and matches the published u.b.2 either way.

– Reference Z. Liu, “Performance Analysis of Stochastic Timed Petri Nets Using Linear Programming Approach”, IEEE Trans. Software Engineering 24(11), 1998, 1014-1030. The constraint families are its Table 1, p. 1022; the bracket statement is its Theorem 3, p. 1021.

See also SPN_SINVARIANTS, SPN_METRICS, SPN_MDD, SOLVER_BA_SPNLP_ANALYZER.

spn_pf(model, options)

PF = SPN_PF(MODEL) PF = SPN_PF(MODEL, OPTIONS) Decide whether a stochastic Petri net has a product-form equilibrium distribution and, when it has, return the per-level factors g_l that MDD_REC and SPN_METRICS take as input.

THIS IS THE PART THE MDD-REC PAPER DECLARES OUT OF SCOPE (FGCS Sec. 3.2). Every other function in api/spn/ receives the g_l already formed; this one derives them from the net, which is what lets a solver reach them.

– The theory, in one paragraph Write I(t), O(t) for the input and output vectors of mode t and lambda_t for its rate constant. Henderson-Taylor and Coleman-Henderson-Taylor show that a net whose firing rate has the form

r_t(m) = lambda_t psi(m - I(t)) / psi(m), m >= I(t)

has invariant measure pi(m) = psi(m) prod_l y_l^{m_l} whenever the positive vector y satisfies COMPLEX BALANCE: reading the distinct vectors appearing as some I(t) or O(t) as the COMPLEXES of the net, the flow into every complex must equal the flow out of it,

sum_{t : O(t)=v} lambda_t y^{I(t)} = ( sum_{t : I(t)=v} lambda_t ) y^v.

Two choices of psi are realisable in LINE’s own rate law, and they are the two this function tests for:

psi = 1 r_t = lambda_t, the rate of a SINGLE-SERVER mode.

pi(m) = prod_l y_l^{m_l}, so g_l(k) = y_l^k.

psi = prod_l 1/m_l! r_t = lambda_t prod_l m_l!/(m_l-I_l)!, MASS ACTION,

reached in LINE through Transition.setFiringRateDependence or, for a mode drawing one token from one place, by infinite-server semantics. pi(m) = prod_l y_l^{m_l}/m_l!, so g_l(k) = y_l^k/k!.

Which one holds is not guessed from the model API: the effective rate LINE would use, lambda_t min(enabling degree, servers) g(m), is EVALUATED at every reachable marking and compared against both laws. A net that matches neither under one common psi is refused by name, never approximated.

– Solving for y Complex balance reads A_lambda Psi(y) = 0 with A_lambda the Laplacian of the weighted digraph on complexes and Psi(y)_v = y^v. That Laplacian is the TRANSPOSED GENERATOR of a Markov chain that hops from complex to complex at the rate of the mode joining them, so its kernel on one linkage class is that chain’s stationary distribution and CTMC_SOLVE returns it – strictly positive exactly when the class is strongly connected, which is weak reversibility. With that positive vector kappa in hand y follows from the LINEAR system in x = log y

(v - v0) x = log kappa_v - log kappa_v0, v, v0 in the same linkage class,

solved in minimum norm. Feinberg’s Deficiency Zero Theorem says this system is consistent for every choice of rate constants when the net is weakly reversible and its deficiency c - l - s is zero, which is why those two numbers are reported; but consistency is CHECKED rather than assumed, so a net of positive deficiency whose particular rates still admit a complex-balanced point is accepted on the evidence.

– The gauge, and why the minimum-norm solution is the canonical one Complex balance fixes y only up to y -> y .* exp(u) for any u orthogonal to the stoichiometric subspace S. Such a shift multiplies pi(m) by exp(u’m), which is CONSTANT on one compatibility class, so every reported measure is invariant under it – but the normalising constant G itself is not, it scales by that constant. A gauge must therefore be FIXED, or the four codebases would report four different G on the same net. The one fixed here is x in the row space of the constraint matrix, i.e. the minimum-norm solution, and it is reached in a form that is unique whichever least-squares primitive a codebase carries: solve (rows*rows’)w = rhs and set x = rows’*w. Any two solutions w of that system give the SAME rows’*w, so the answer does not depend on how the rank-deficient solve breaks its tie.

– Input MODEL : a Network holding Places and Transitions OPTIONS : struct, fields

bound - per-place-level token bound, passed to SPN_MDD tol - relative tolerance of the rate-law and complex-balance

checks (default 1e-9)

verbose - print the certificate (default false)

– Output PF : struct with fields

g - 1 x L cell, g{l}(k+1) = g_l(k), ready for MDD_REC y - 1 x L positive vector solving complex balance kind - ‘geometric’ or ‘massaction’, the psi that was found complexes - C x L matrix of the distinct complexes deficiency - c - l - s linkage - number of linkage classes srank - rank of the stoichiometric subspace weaklyreversible - true when every linkage class is strongly connected residual - relative complex-balance residual at y mdds, info - the reachable set and metadata SPN_MDD returned

– Reference J. L. Coleman, W. Henderson, P. G. Taylor, “Product form equilibrium distributions and a convolution algorithm for stochastic Petri nets”, Performance Evaluation 26(3), 1996. M. Feinberg, “Complex balancing in general kinetic systems”, Arch. Rational Mech. Anal. 49, 1972. D. F. Anderson, G. Craciun, T. G. Kurtz, “Product-form stationary distributions for deficiency zero chemical reaction networks”, Bull. Math. Biol. 72, 2010.

See also MDD_REC, SPN_METRICS, SPN_MDD, SPN_CONV.

spn_sinvariants(sn, init)

INV = SPN_SINVARIANTS(SN) INV = SPN_SINVARIANTS(SN, INIT) Minimal-support S-invariants (P-invariants) of a stochastic Petri net, and the load vector V = S m0.

An S-invariant is a non-negative left null vector of the incidence matrix, U’ C = 0, so U’ m is conserved by every firing. The minimal-support ones form a basis of all of them and are what the convolution algorithm SPN_CONV decomposes the reachability set along; SPN_MDD uses a single positive invariant for a much weaker purpose, to bound each place a priori.

FARKAS’ ALGORITHM, on [C | I]: for each transition column in turn, keep the rows that already annihilate it and add, for every pair of rows of opposite sign in it, the positive combination that cancels it; then drop every row whose support strictly contains another’s, which is what leaves the minimal supports. Rows are kept in integer arithmetic and divided by their gcd, so a multiplicity is never lost to rounding and two invariants that differ only by a positive scale are the same row.

ARC MULTIPLICITIES MUST BE INTEGRAL. A fractional arc has no Petri-net meaning and would make the gcd normalisation and the ILP-free convolution both wrong, so it is refused rather than rounded.

Levels are the (place, class) pairs of SPN_MDD, place-major, so the invariants come out in the coordinates the decision diagram and SPN_CONV both use.

– Input SN : a NetworkStruct holding Places and Transitions INIT : initial tokens per place level, place-major; omitted or [] takes them

from the reference station of each closed class, as SPN_MDD does

– Output INV : struct with fields places (node indices), S (one row per invariant, one

column per place level), V = S m0 and m0

– Reference S. Balsamo, A. Marin, I. Stojic, “Computation of the normalising constant for product-form models of distributed systems with synchronisation”, Future Generation Computer Systems 111 (2020) 475-490, Sec. 3.1.

See also SPN_CONV, SPN_MDD.

spn_rec_enabled(mdds, g, mde, nplacelevels)

EN = SPN_REC_ENABLED(MDDS, G_L, MDE, NPLACELEVELS) Enabling-degree distribution of one mode of a product-form stochastic Petri net, by the masked MDD-rec recursion.

The enabling degree of a mode in marking m is

e(m) = min_{l : I_l > 0} floor(m_l / I_l),

zero when any inhibitor threshold is met. P(e >= k) is therefore the mass of the marking subset in which EVERY input level holds at least k*I_l tokens and no inhibitor fires, which is a per-level restriction and so exactly what MDD_REC_MASKED computes: the paper’s second modified recurrence is the same walk under a different mask, not a second algorithm.

– Input MDDS : MDD.toStruct of the reachable set built by SPN_MDD G_L : 1 x K cell of per-level product-form factors MDE : one entry of INFO.modes as returned by SPN_MDD NPLACELEVELS : how many leading levels are place levels – Output EN : struct with fields

ge - ge(k+1) is the mass of {e >= k}; ge(1) is the whole set eq - eq(k+1) is the mass of {e == k}, i.e. ge(k+1)-ge(k+2) maxDegree - the largest enabling degree the place bounds permit, E_j

The masses are UNNORMALISED, as in the paper; divide by G from MDD_REC for probabilities. SPN_METRICS does that and turns them into the transition measures.

– Reference S. Balsamo, A. Marin, I. Stojic, “Computation of the normalising constant for product-form models of distributed systems with synchronisation”, Future Generation Computer Systems 111 (2020) 475-490, Sec. 5.3.

See also MDD_REC_MASKED, SPN_METRICS, SPN_MDD.

spn_metrics(mdds, g, info)

MET = SPN_METRICS(MDDS, G_L, INFO) Stationary measures of a product-form stochastic Petri net from the MDD-rec masses.

n(P_j) = sum_k k P(m_j = k) mean tokens u(P_j) = 1 - P(m_j = 0) place utilization u(T_j) = P(e_j >= 1) transition utilization x(T_j) = sum_k min(k, c_j) W(T_j) P(e_j = k) throughput x(P_j) = sum_T I_j(T) x(T) tokens removed per unit time

ONE DEVIATION FROM THE PAPER’S x(T_j), AND IT IS A GENERALISATION. The paper writes x(T_j) = sum_k k W(T_j) P(e_j = k), which is INFINITE-SERVER firing semantics – every enabling set fires in parallel. LINE’s own rate law is min(enabling degree, nmodeservers) * W(T), so c_j above is the mode’s server count: c_j = 1 recovers single-server semantics, x = W(T) P(e >= 1), and c_j = Inf recovers the paper’s formula exactly. Using the paper’s form for a single-server mode would report a throughput that grows with the token population of a net whose transition can only fire one set at a time.

The measures come out of ONE reachable set and ONE set of g_l, so they are mutually consistent by construction: no per-measure fixed point, no iteration.

– Input MDDS : MDD.toStruct of the reachable set built by SPN_MDD G_L : 1 x K cell of per-level product-form factors g_l(v) INFO : the metadata SPN_MDD returned alongside the diagram – Output MET : struct with fields G, tokens, placeUtil, placeTput, modeUtil,

modeTput and marginal, the per-level P(m_l = k)

– Reference S. Balsamo, A. Marin, I. Stojic, “Computation of the normalising constant for product-form models of distributed systems with synchronisation”, Future Generation Computer Systems 111 (2020) 475-490, Sec. 3.1 and Sec. 5.3.

See also MDD_REC, MDD_REC_MARGINAL, SPN_REC_ENABLED, SPN_MDD.

spn_mdd(model, options)

[MDDS, DESC, INFO] = SPN_MDD(MODEL) [MDDS, DESC, INFO] = SPN_MDD(MODEL, OPTIONS) Decision-diagram reachable set and Kronecker rate descriptor of a stochastic Petri net, so that MDD_MCD can analyse it.

Levels are of two kinds:

place levelsone per (Place, class) pair, holding a token count. A

multiclass net therefore has P*R of them, ordered place-major, so level (p-1)*R+k is class k in place p.

phase levelsone per mode whose firing time has more than one phase,

holding the phase the running server occupies.

The rate structure factorises exactly under single-server firing semantics: a mode fires at a constant rate whenever every input level holds its enabling multiplicity and no inhibitor level has reached its threshold, so

W_l^e[i, i + fire(l) - enab(l)] = 1 for enab(l) <= i < inhib(l)

at every place level. A phase-type mode contributes two event families on its phase level, the internal phase changes D0 (marking unchanged) and the firings D1 (marking moved), each gated by the same per-level enabling indicators. Both are products of per-level terms, which is what Eq. 1 of the paper requires.

– Input MODEL : a Network holding Places and Transitions OPTIONS : struct, fields

bound - per-place-level token bound (scalar or 1 x P*R);

inferred from a place invariant when omitted

phmemory - ‘exact’ (default) or ‘resume’, see the note below descriptor - build the Kronecker rate descriptor (default true).

Pass false for the MDD-rec route, which reads only the reachable set: see the note below

verbose - print the net summary (default false)

– Output MDDS : MDD.toStruct of the reachable set DESC : Kronecker descriptor for MDD_MCD, with fields K, domain, events,

levelkind (1 place, 2 phase) and, when one exists, invariant = struct(‘weights’, w, ‘value’, v) with zero weight on the phase levels

INFOstruct with fields places, placenames, classnames, levelname,

levelkind, modes, init, mdd, phaseof

– Phase-type firing and the memory policy LINE discards a running server’s phase when its mode becomes disabled (State/afterGlobalEvent zeroes the server phases on the disable action), i.e. preemptive repeat. Resetting a mode’s phase is then triggered by a JOINT condition on the place levels, which is not a product of per-level terms and has no Kronecker form. What this descriptor encodes is preemptive resume: a disabled mode’s phase freezes and continues when the mode is re-enabled. The two policies coincide exactly when a mode is never disabled while running, so reachability records, for free, whether any phase-type mode was ever found disabled. OPTIONS.phmemory=’exact’ (the default) errors when one was; ‘resume’ proceeds deliberately with the resume semantics.

– Reachable-set-only mode (OPTIONS.descriptor = false) MDD_REC and SPN_METRICS need the reachable set and the metadata, never the rate descriptor: the product form supplies the rates. So the restrictions below that exist only because a Kronecker form must factorise per level – marking-dependent firing rates and multi-server modes drawing from several places – are lifted, and DESC comes back carrying K, domain, levelkind and the invariant but no events. In exchange the firing times must be EXPONENTIAL: a phase level is a descriptor device, and a product-form marking process is memoryless in the marking alone. Immediate modes stay refused, since a vanishing marking carries no probability and would have to be eliminated from the reachable set first.

– Other restrictions (each is an error, never a silent approximation) No immediate transitions (they make vanishing states, which must be eliminated before a Kronecker rate descriptor exists) and no marking-dependent firing rates. A multi-server mode is accepted only when its enabling touches ONE level, because the enabling degree min_l floor(m(l)/enab(l)) is otherwise not a product of per-level terms.

See also: mdd_mcd(), mdd_descriptor(), mdd_reachset().

spn_conv(S, V, g)

G = SPN_CONV(S, V, G_L) G = SPN_CONV(INV, G_L) Normalising constant of an S-invariant reachable product-form stochastic Petri net, by convolution over the invariant load vector.

With S the minimal-support S-invariant matrix and V = S m0 the load vector, the reachability set of an S-INVARIANT REACHABLE net is exactly {m >= 0 : S m = V}, and conditioning on the marking of one place partitions it (Lemma 5.1 of the MDD-rec paper). Writing G_j(W) for the mass of the markings supported on the first j places with S m = W,

G_0(W) = [W == 0], G_j(W) = sum_i g_j(i) G_{j-1}(W - i S_j),

and G = G_n(V). On a net whose only invariant is “the tokens are conserved” this is Buzen’s convolution for a closed queueing network, one place per station.

NO ILP IS SOLVED. Coleman-Henderson-Taylor obtain the marking set M_p(P’,W) from the feasibility of an integer program so that the sum skips the terms that contribute nothing. Here the sum simply runs over i whose residual W - i S_j stays non-negative and the recursion returns zero on an infeasible residual, which gives the same value: the ILP is an optimisation of the enumeration, not part of the definition. Memoising on (j, W) keeps the walk over the reachable residuals rather than over all of them.

S-INVARIANT REACHABILITY IS NOT CHECKED, and cannot be cheaply: no algorithm is known that decides it without generating the reachability set. On a net that fails it, {m : S m = V} is strictly larger than the reachable set and this returns a normalising constant over unreachable markings too, which is why MDD_REC – which walks the reachable set itself – is the general algorithm and this one the special case. Compare the two on a new net before trusting this one on it.

– Input S : S(i,p), the minimal-support S-invariants, one row per invariant; or the

struct SPN_SINVARIANTS returned, in which case V holds the factors

V : the load vector S m0, one entry per invariant G_L : 1 x n cell, G_L{p}(i+1) is g_p(i), the product-form factor of i tokens

in place level p; its length bounds the marking of that level

– Output G : the normalising constant

– Reference J. Coleman, W. Henderson, P. Taylor, “Product form equilibrium distributions and a convolution algorithm for stochastic Petri nets”, Performance Evaluation 26(3), 1996, 159-180; the point of comparison for MDD-rec in S. Balsamo, A. Marin, I. Stojic, Future Generation Computer Systems 111 (2020) 475-490, Sec. 5.1.

See also MDD_REC, SPN_SINVARIANTS, SPN_MDD.

spn_bench_mdd(OPTIONS)

OUT = SPN_BENCH_MDD() OUT = SPN_BENCH_MDD(OPTIONS) Cost profile of the decision-diagram aggregation on a SCALABLE stochastic Petri net, the flexible-manufacturing-system (FMS) family of G. Ciardo, K. Trivedi, “A decomposition approach for stochastic reward net models”, Perform. Eval. 18(1), 1993, which is the standard scaling benchmark for decision-diagram state-space methods.

The net has three machine lines, one job class each, that take a part from a common pool, process it and return it, plus an assembly stage that consumes one finished part of every class at once. The reachable set grows steeply in the part count N while the number of levels stays fixed, which is the regime the aggregation is built for: |S| grows polynomially of degree 3R in N and sum_k |M_k| does not.

– Input OPTIONS : struct, fields

parts - vector of part counts to sweep (default 1:5) maxstates- skip a row whose exact solve would exceed this many

states (default 200000); the aggregation still runs

verbose - print the table as it is produced (default true)

– Output OUT : struct array, one entry per part count, with fields parts, nstates,

levels, sumMk, compression, iters, t_reach, t_solve, noAggregation and, when the exact solve was run, t_exact and maxerr

– Remarks The exact reference is the explicit generator assembled through MDD/index, which is why MAXSTATES gates it: the aggregation is the cheap side and the point of the benchmark is to watch the two diverge.

See also: spn_mdd(), mdd_mcd(), mdd_bench_mcd().