api.mdd

mdd_rec_masked(mdds, g, mask)

G = MDD_REC_MASKED(MDDS, G_L) G = MDD_REC_MASKED(MDDS, G_L, MASK) Unnormalised mass of the masked subset of a reachable set held in a decision diagram, i.e. Algorithm 1 of the MDD-rec paper under a per-level restriction.

A product-form model has P(s) = (1/G) prod_k g_k(s_k) over its levels, and

G = sum_{s in S} prod_k g_k(s_k).

Summing state by state is exponential and numerically unstable. MDD-rec instead walks the diagram that already encodes S, accumulating the unnormalised mass of each node ONCE (Def. 4.4, Algorithm 1):

M(<l.p>) = sum_{v in S_l} g_l(v) * M(<l.p>[v]), M(TRUE)=1, M(FALSE)=0

so the cost is O(sum_l |nodes_l| * |S_l|) rather than O(|S|), and G = M(root).

THE MASK is how Sec. 5.3 of the paper computes measures. Restricting the sum at level l to a subset of its local values gives the unnormalised mass of the corresponding subset of S, so P(m_l = k) and P(e_j >= k) are the same recursion under a different mask rather than three separate algorithms.

– Input MDDS : MDD.toStruct of the reachable set (level 1 is the root) G_L : 1 x K cell, G_L{l}(v+1) is g_l(v), the per-level product-form factor MASK : 1 x K cell of logical row vectors; [] admits everything, which is the

plain MDD-rec of Algorithm 1 and returns the normalising constant

– Output G : the unnormalised mass of the masked subset

– Note The g_l themselves, and the test that the model has a product form at all, are the caller’s: the paper declares that out of scope (Sec. 3.2). Passing g_l that do not describe a product-form model returns a number that is not the normalising constant of anything, and nothing here can detect it.

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

See also MDD_REC, MDD_REC_MARGINAL, SPN_REC_ENABLED, MDD_REACHSET.

mdd_rec_marginal(mdds, g, l)

MASS = MDD_REC_MARGINAL(MDDS, G_L, L) Unnormalised masses of {s in S : s_L = k}, one per local value k of level L.

Divided by the normalising constant these are P(m_l = k) of Sec. 5.3 of the MDD-rec paper: the mean occupancy of a level is sum_k k * P(m_l = k), and its utilization 1 - P(m_l = 0).

– Input MDDS : MDD.toStruct of the reachable set G_L : 1 x K cell, G_L{l}(v+1) is g_l(v) L : level index, 1..K – Output MASS : 1 x domain(L) unnormalised masses, MASS(k+1) for local value 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. 5.3.

See also MDD_REC, MDD_REC_MASKED, SPN_METRICS.

mdd_rec(mdds, g)

G = MDD_REC(MDDS, G_L) Normalising constant of a product-form model whose reachable set is held in a decision diagram: G = sum_{s in S} prod_l g_l(s_l), by Algorithm 1 of the MDD-rec paper.

FORMALISM-AGNOSTIC. Nothing here knows what a level is: on the lattice sum_k s_k = n of a closed queueing network this collapses to Buzen’s convolution (paper Appendix B), and on an S-invariant reachable Petri net to the Coleman-Henderson-Taylor convolution, SPN_CONV (Sec. 5). Unlike either, it needs only that the reachable set be finite and encoded – no lattice, no S-invariant reachability.

– Input MDDS : MDD.toStruct of the reachable set G_L : 1 x K cell, G_L{l}(v+1) is g_l(v) – Output G : the normalising constant

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

See also MDD_REC_MASKED, MDD_REC_MARGINAL, SPN_METRICS, SPN_CONV.

mdd_reachset(domain, init, nextfun)

MDD = MDD_REACHSET(DOMAIN, INIT, NEXTFUN) Generate and store the reachability set into a quasi-reduced ordered MDD, using the decision diagram itself as the visited-state store, after A.S. Miner, G. Ciardo, “Efficient Reachability Set Generation and Storage Using Decision Diagrams”, ICATPN 1999, LNCS 1639, pp.6-25.

– Input DOMAIN : 1 x K vector of per-level local-state counts (values 0..domain(k)-1) INIT : 1 x K initial global state (0-based local values) NEXTFUN : function handle s -> T, returning an (m x K) matrix whose rows are

the successor states of s under the next-state function

– Output MDD : an MDD object holding every state reachable from INIT

– Remarks This is the basic (explicit-frontier) realisation: a breadth-first search enumerates successors while the MDD provides the O(K) membership test that replaces the usual explicit visited hash. The stored set lives entirely in the MDD (O(number of nodes) memory); only the transient BFS frontier is held explicitly. Symbolic image computation / saturation, which removes the explicit frontier too, is the natural next step but is out of scope here.

See also: MDD, mdd_closedqn().

mdd_mcd(mdds, desc, options)

OUT = MDD_MCD(MDDS, DESC) OUT = MDD_MCD(MDDS, DESC, OPTIONS) Miner-Ciardo-Donatelli approximate stationary analysis: solve a structured CTMC whose EXACT reachable state space is stored in a decision diagram, by building and iterating K level-CTMCs (a decision-diagram-guided aggregation), after A.S. Miner, G. Ciardo, S. Donatelli, “Using the exact state space of a Markov model to compute approximate stationary measures”, SIGMETRICS 2000.

The method never forms the |S|-state generator or probability vector. It keeps one CTMC per decision-diagram level k, over states M_k = {(p,i_k)} with p a level-k node and i_k a local state on a non-null arc, and iterates the coupled system to a fixed point. The single approximation (Eq. 5) is Pr{i_k | alpha} = Pr{i_k | p}: the local-state law at level k depends only on the node p, not the full path above it – justified by the exact reachability the diagram encodes. For product-form models the method is EXACT (paper Sec. 5), so on a single-class closed QN it reproduces SolverCTMC.

Orientation note: the paper indexes levels K (top/root) down to 1 (bottom/ terminal); the MDD class uses level 1 as the root. This routine works in the paper’s orientation, mapping paper level k to MDD level (K+1-k), i.e. to station (K+1-k).

– Input MDDS : struct from MDD.toStruct (the reachable set, MDD orientation) DESC : Kronecker rate descriptor from mdd_descriptor OPTIONS : struct, fields tol (1e-12), maxiter (500), verbose (false),

initpik (optional 1 x K cell of level warm-start vectors)

– Output OUT : struct with fields

QLen - 1 x K mean jobs per station X - 1 x K per-station throughput U - 1 x K utilisation pik - 1 x K cell, level-k stationary vectors over M_k Mrows - 1 x K cell, Mrows{k}(r,:) = [node local] of M_k row r levelSizes - 1 x K, |M_k| iters - fixed-point iterations performed

mdd_ps(mu, P, servers, N, options)

DESC = MDD_PS(MU, P, SERVERS, N, OPTIONS) Kronecker rate descriptor of a single-class closed queueing network whose stations are processor-sharing or infinite-server, with phase-type service.

Under processor sharing every job at a station is in service at once, each holding its own phase, so naming a single in-service phase (what MDD_DESCRIPTOR does, which is non-preemptive semantics) cannot represent the state. The local state here is instead the PER-PHASE COUNT vector v = (v_1,…,v_h), v_a jobs in phase a, with n = sum(v) jobs present. That is still a per-station quantity, so every event stays a product of per-level terms and the Kronecker form of Eq. 1 survives.

With one server shared by n jobs each job advances at rate 1/n, so from local state v with n = sum(v):

internal v -> v - e_a + e_b at v_a * D0(a,b) / n (a ~= b) departure v -> v - e_a at v_a * t(a) / n * P(i,j) arrival v -> v + e_b at pie(b)

An infinite-server (delay) station is the same without the 1/n scaling. For h = 1 the departure rate collapses to n*mu/n = mu at PS and to n*mu at IS, reproducing the usual single-server and delay rate laws.

The local domain is the number of compositions of 0..N over h phases, C(N+h,h), against 1+N*h for the non-preemptive encoding: the price of tracking every job’s phase rather than one.

– Input MU : 1 x K station service rates (ignored where OPTIONS.proc gives a law) P : K x K routing matrix (row-stochastic) SERVERS : 1 x K servers per station, 1 (PS) or Inf (IS); no other value has

a per-phase-count encoding here

N : closed population OPTIONS : struct, fields proc and pie as in MDD_DESCRIPTOR – Output DESC : descriptor for MDD_MCD, with K, N, domain, mu, servers, P,

nphases, valuemap, init, nextfun, events, and the population invariant through valuemap

mdd_profile(M, N): MDD_PROFILE(M, N) Profile the MDD reachability-set build for a single-class closed cyclic network with M stations and N jobs, printing the top hotspots by total time. Temporary developer harness (not part of the public API).

mdd_descriptor(mu, P, servers, N, options)

DESC = MDD_DESCRIPTOR(MU, P, SERVERS, N) DESC = MDD_DESCRIPTOR(MU, P, SERVERS, N, OPTIONS) Kronecker rate descriptor of a single-class closed queueing network, for the Miner-Ciardo-Donatelli approximate-aggregation solver (mdd_mcd), after A.S. Miner, G. Ciardo, S. Donatelli, “Using the exact state space of a Markov model to compute approximate stationary measures”, SIGMETRICS 2000.

The transition rate matrix is expressed compositionally as: R = sum_e ( kron_{k} W_k^e ) restricted to the reachable set,

with W_k^e[i_k,j_k] = lambda_k^e[i_k] * Prob_k^e(i_k,j_k) (Eq. 1). Each level k is a station and each event e is a completion at station a routed to b.

– Exponential stations The local state is the population alone:

W_a^e[i,i-1] = mu(a)*min(i,servers(a))*P(a,b) (i >= 1) – departure W_b^e[i,i+1] = 1 (i <= N-1) – arrival W_l^e = I (l ~= a,b)

– Phase-type stations The local state is the PAIR (population, phase of the job in service), encoded in one level rather than two. Splitting them does not work: on a completion routed into station b the phase at b restarts only when b was empty, a joint condition on b’s two components, which is not a product of per-level terms. Merging them keeps every event local. The encoding is

index 0 : station empty index 1 + (n-1)*h + (a-1) : n jobs present, job in service in phase a

so the domain is 1 + N*h and h = 1 reproduces the exponential encoding index = n exactly. With exit vector t = -D0*1 and entry law pie,

departure (n,a) -> (n-1,b) at t(a)*P(a,b)*pie(b) for n >= 2
(1,a) -> 0 at t(a)*P(a,b)

arrival 0 -> (1,b) at pie(b)
(m,c) -> (m+1,c) at 1 for m >= 1

internal (n,a) -> (n,b) at D0(a,b) for n >= 1, a ~= b

– Input MU : 1 x K station service rates. Entry i is ignored when station i is

given a phase-type law through OPTIONS.proc.

P : K x K routing matrix (row-stochastic) SERVERS : 1 x K servers per station (Inf for delay/IS) N : closed population OPTIONS : struct, fields

proc - 1 x K cell; entry i empty for an exponential station, or
a Markovian distribution object, or {D0,D1}

pie - 1 x K cell of entry laws; taken from the object, or the
stationary entry law of {D0,D1}, when omitted

sched- 1 x K cell of SchedStrategy values. Optional, and only
consulted to REJECT a phase-type law at a preemptive-resume station (see the restriction below). Supply it whenever the stations are not all non-preemptive, because the descriptor otherwise has no way to detect that case.

– Output DESC : struct with fields

K, N, domain, mu, servers, P nphases - 1 x K phases per station (1 when exponential) valuemap - 1 x K cell, valuemap{i}(idx+1) = jobs at station i init - initial local index per station nextfun - successor function over local indices, for mdd_reachset events - struct array, one per routing pair, with fields

a, b, lev (levels touched), W (matrices at lev)

– Restrictions A phase-type station must be single-server. With c > 1 or an infinite server the local state would have to count jobs per phase rather than name one phase, which is a different and much larger encoding.

A phase-type station must also be NON-preemptive. The composite level names the phase of the one job in service and restarts it at pie when the next job starts; under preemptive resume an arrival suspends that job and its phase has to be remembered, so the local state would need a stack of phases. This matters for LCFSPR, which is BCMP type 2 and stays product-form under general service: that insensitivity is real but is NOT reachable through this encoding. Exponential service is unaffected, preemption being immaterial by memorylessness. Pass OPTIONS.sched to have the case rejected rather than silently modelled as non-preemptive.

mdd_demo_mcd()

MDD_DEMO_MCD Demonstration and validation of the Miner-Ciardo-Donatelli approximate aggregation (mdd_mcd) on single-class closed exponential queueing networks. For product-form models the method is EXACT (SIGMETRICS 2000, Sec. 5), so it must reproduce the exact solution (mdd_closedqn / SolverCTMC) to machine precision while solving only K small level-CTMCs instead of the full chain.

mdd_closedqn(mu, P, servers, N, options)

OUT = MDD_CLOSEDQN(MU, P, SERVERS, N) OUT = MDD_CLOSEDQN(MU, P, SERVERS, N, OPTIONS) Solve a single-class closed exponential queueing network whose CTMC state space (reachable occupancy vectors) is stored in a Multi-valued Decision Diagram (MDD) instead of an explicit state list.

The reachable set is generated with MDD_REACHSET and the generator matrix is assembled using the MDD’s O(K) state-indexing (MDD/index), so no explicit (|S| x width) state matrix is ever materialised – the diagram is the store.

– Input MU : 1 x M vector of per-station exponential service rates P : M x M Markovian routing matrix (row-stochastic, irreducible) SERVERS : 1 x M number of servers per station (Inf for a delay/IS station) N : closed population (number of jobs) OPTIONS : (optional) struct; OPTIONS.verbose prints a storage summary,

OPTIONS.mdd reuses an already-built reachable set (as returned in OUT.mdd) instead of regenerating it

– Output OUT : struct with fields

mdd - the MDD holding the reachable occupancy set Q - sparse CTMC generator, rows aligned to MDD/index order pi - stationary distribution (1 x |S|) states - |S| x M occupancy matrix in MDD/index order QLen - 1 x M mean number of jobs per station U - 1 x M utilisation (busy servers / servers; mean busy for IS) X - 1 x M per-station throughput stats - MDD storage statistics (see MDD/stats) times - phase timings in seconds: reach (reachable-set build, 0

when OPTIONS.mdd is supplied), gen (generator assembly), solve (ctmc_solve), metrics (performance measures)

– Remarks For single-class exponential stations the aggregated (occupancy) chain is exact: the rate from n to n-e_i+e_j is mu(i)*min(n_i,servers(i))*P(i,j) for n_i>0. This matches LINE’s SolverCTMC on the same model.

See also: MDD, mdd_reachset(), solver_ctmc(), ctmc_solve.

mdd_bench_mcd(Ks, Ns, options)

T = MDD_BENCH_MCD(KS, NS) T = MDD_BENCH_MCD(KS, NS, OPTIONS) Cost profile of the Miner-Ciardo-Donatelli aggregation (mdd_mcd) against the exact CTMC solution of the same model, over a grid of station counts KS and populations NS on single-class closed cyclic exponential networks.

Three costs are separated because they are not interchangeable:

t_reachbuilding the exact reachable set into the MDD. Both methods need: it – mdd_mcd reads the diagram, so it can never beat this term.

t_exact : explicit generator assembly + ctmc_solve, given the diagram. t_mcd : descriptor + the K coupled level-CTMC fixed point, given the

diagram.

The crossover reported is therefore t_exact/t_mcd (the solve-phase speedup, the quantity the method is actually about) alongside the end-to-end ratio (t_reach+t_exact)/(t_reach+t_mcd), which is what a user observes.

Storage is compared as |S| (exact stationary vector) against sum_k |M_k| (the K level vectors that mdd_mcd stores instead).

– Input KS : vector of station counts (default [3 4 5 6]) NS : vector of populations (default [4 8 12 16 20]) OPTIONS : struct, fields

maxstates - skip the exact solve above this |S| (default 3e4) ctmcstates - also time LINE’s SolverCTMC below this |S|

(default 0 = never; SolverCTMC re-derives the state space itself and is the end-user baseline)

mcdopts - options struct forwarded to mdd_mcd reps - repetitions per phase, median reported (default 3) ctmcmethod - forced ctmc_solve method for the exact baseline;

‘direct’ keeps the sparse factorization above GMRES_MIN_STATES (6000), where the default path switches to GMRES and its cost stops tracking |S|

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

– Output T : table with one row per (K,N) case

Temporary developer harness (not part of the public API).

class MDD

Bases: handle

MDD - Quasi-reduced ordered Multi-valued Decision Diagram.

Compact symbolic store for a set of discrete states, after: A.S. Miner, G. Ciardo, “Efficient Reachability Set Generation and Storage Using Decision Diagrams”, ICATPN 1999, LNCS 1639, pp.6-25.

A global state is a K-tuple of local state values (one per level/ submodel), state(k) in {0,…,domain(k)-1}. The set is stored as a directed acyclic graph with K variable levels plus a terminal level: level 1 is the top (root), a node at level k has domain(k) outgoing arcs to level k+1 nodes, and a state belongs to the set iff its path of arcs reaches the TRUE terminal. Canonicity is enforced by a per-level unique table (no duplicate nodes) and by collapsing the all-FALSE node to the FALSE terminal. Storage is O(#nodes), typically O(K * #local-states), instead of O(|S|) as in an explicit state list.

Two constant terminals encode the boolean value of a completed path:: TERM_FALSE = 0 (empty subgraph / state not in set) TERM_TRUE = -1 (state in set)

Arcs of a level-k node hold ids of level-(k+1) nodes when k<K, or a terminal (TERM_TRUE/TERM_FALSE) when k==K.

Example

m = MDD([3 3 3]); % 3 levels, local values 0..2 m.insert([0 1 2]); m.insert([2 0 0]); m.member([0 1 2]) % -> true m.cardinality() % -> 2 m.index([2 0 0]) % -> 0-based rank in the set S = m.enumerate(); % all states, one per row, in index order

See also: mdd_reachset(), ctmc_ssg().

Constructor Summary

MDD(domain): MDD(domain) creates an empty set over the given per-level domains. domain is a 1 x K vector; state values at level k are 0..domain(k)-1.

Property Summary

K: number of variable levels

TERM_FALSE

empty subgraph

Type:: terminal node “0”

TERM_TRUE

completed path is accepted

Type:: terminal node “1”

domain: 1 x K, domain(k) = number of local states at level k

node: 1 x K cell, node{k} is (nnodes_k x domain(k)) of child ids

root: id of the top (level-1) node; TERM_FALSE for the empty set

uniq: 1 x K cell of dictionary(string->double), unique table per level

Method Summary

cardinality(): CARDINALITY() returns |S|, the number of stored states.

compact(): COMPACT() reclaims dead nodes left by the append-only build, rebuilding the level tables and unique tables so that only nodes reachable from the root remain. Membership/index/enumerate are unchanged. (A production MDD would reference-count instead and never accumulate dead nodes; this is the basic sweep.)

disp()

enumerate(): ENUMERATE() returns all stored states as rows, in index() order.

static fromStates(domain, S): MDD.fromStates(domain, S) builds an MDD from a matrix S whose rows are 0-based state tuples over the given per-level domains.

index(state): INDEX(state) returns the 0-based lexicographic rank of state among the stored set (level 1 most significant), or -1 if state is not in the set. This is a bijection S <-> {0,…,|S|-1}, so a generator matrix can be assembled without an explicit state list.

insert(state): INSERT(state) adds a K-tuple (row vector, 0-based values) to the set.

member(state): MEMBER(state) returns true iff state is in the set (O(K)).

stats(): STATS() returns a struct describing the storage of the current set. Only nodes reachable from the root are counted (the build is append-only, so superseded nodes may linger in the tables until COMPACT is called; see COMPACT):

numStates - |S| numNodes - reachable non-terminal nodes nodesPerLevel- 1 x K reachable node counts liveNodes - reachable nodes (= numNodes) tableNodes - nodes physically held in the tables (incl. dead) mddInts - integers in the reachable arc arrays (footprint) explicitInts - integers an explicit state list needs (|S|*K) compression - explicitInts / mddInts

toStruct(): TOSTRUCT() exports the diagram as plain arrays for downstream algorithms (e.g. mdd_mcd), trimmed to live rows. Fields:

K, domain, root - as the properties nnodes - 1 x K live node count per level node - 1 x K cell, node{k} is nnodes(k) x domain(k)

of child ids (level k+1 ids, or terminals at k=K)