solvers.SSA

ssa_warmupfrac(sn, options)

WF = SSA_WARMUPFRAC(SN, OPTIONS)

@brief Effective warmup (initial-transient) discard fraction of an SSA run @param sn Network structure @param options Solver options; options.config.warmupfrac may be unset ([]) @return wf Fraction of the requested samples (firings) whose tallies are discarded

An explicit options.config.warmupfrac is honoured as given, 0 included. When it is unset the rule is the same in all four codebases: a STEADY-STATE run (infinite options.timespan end) of a model with at least one OPEN class drops the first 20% of its samples, the LDES default, and any other run drops none. An open model starts EMPTY, far from its stationary load, so on a heavily loaded queue the first stretch of the path holds fewer jobs than the stationary mean and biases QN low. A closed model starts at a feasible state of its fixed population and the transient is short against the run, so it keeps every sample. A transient run (finite timespan) and a sample path keep every sample, since the initial transient is what they are asked for.

ssa_reply_held(sn, tranSysState)

HELD = SSA_REPLY_HELD(SN, TRANSYSSTATE) Time-averaged servers held for a pending REPLY signal, per (station, calling class).

A job of a class r with SN.SYNCREPLY(r) >= 0 leaves its caller for the callee but KEEPS its server there until the REPLY returns. The hold is counted in the reply block at the tail of the caller’s local state (State.replyBlockInfo), not in its buffer or servers, so State.toMarginal does not see it. SolverCTMC adds the same block mean to QLen and Util (solver_ctmc_analyzer); this is its sample-path twin, read from the recorded trajectory: TRANSYSSTATE{1} is the cumulative time and TRANSYSSTATE{1+isf} the state of stateful node isf held over each step. HELD is all-zero for a model without synchronous calls.

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

ssa_nrm_refusal(sn)

[OK, REASON] = SSA_NRM_REFUSAL(SN) Whether an EXPLICIT options.method=’nrm’ can answer this model.

NOT SSA_NRM_ELIGIBLE, and the difference is the whole point of having two predicates. SSA_NRM_ELIGIBLE asks whether the NRM should be PREFERRED, and it is deliberately wide: SOLVER_SSA_ANALYZER consults it on the ‘default’ and ‘parallel’ paths to pick an engine. This one asks what a GATE must ask, which is whether the name the caller typed will produce an answer.

The two differ because the explicit ‘nrm’ arm of SOLVER_SSA_ANALYZER FALLS BACK to the serial engine, with a warning, for reneging patience, balking, phase-type service, global dependence and cache retrieval. A model carrying any of those still gets a correct answer under ‘nrm’, so refusing it here would withdraw a row the solver honours – the report would be wrong in the opposite direction.

What is left are the TWO conditions that raise rather than fall back: a scheduling discipline the reaction network has no form for (solver_ssa_analyzer_nrm.m:24, ‘solver_ssa_analyzer_nrm:UnsupportedPolicy’), and a marking-dependent firing rate (SSA_FIRINGDEP_REFUSAL), which the explicit arm refuses by name as C++ does, so that ‘nrm’ never answers with the nominal rate. ‘default’ sends the latter to the serial engine instead.

The C++ port gates on the wider set ON PURPOSE and is not diverging: its NRM has no fallback arm at all, so every one of the six conditions raises there. Each gate states its own engine’s reach.

ssa_nrm_guards(sn)

G = SSA_NRM_GUARDS(SN) The per-feature tests the NRM engine’s dispatch and gate both read.

ONE BODY, THREE CALLERS. SOLVER_SSA_ANALYZER reads the fields INDIVIDUALLY, because its explicit ‘nrm’ arm falls back to the serial engine with a warning that names the offending construct, so it has to know WHICH test failed. SSA_NRM_ELIGIBLE reads them together to answer “should the NRM be preferred”. SSA_NRM_REFUSAL reads only sched and firingdep, because those are the two conditions the explicit arm does not fall back on.

These were local functions of SOLVER_SSA_ANALYZER until 2026-09-05, which is why the eligibility test could not be asked from a gate. Moving them here is what lets all three callers share one definition instead of three copies.

Fields (all true = the NRM can take it):

sched every station’s discipline has a reaction form phase non-exponential service only where every job present is in service renege no station renegs with non-exponential patience balk no balking strategy outside QUEUE_LENGTH gd no global (Whittle) dependence cache every cache configuration (always true; see below) fcr every finite capacity region rule (always true; see below) spn a Petri net the NRM’s reaction builder can read (see below) immfeed no EFFECTIVE immediate feedback (sn_has_immfeed: declared AND a

self-loop in sn.rt, so a declaration nothing can exercise is inert)

reply no synchronous call (REPLY signal) holds a server anywhere sdr no state-dependent (SDR) routing, which only the serial engine reads firingdep no marking-dependent firing rate (setFiringRateDependence), which

the reaction builder has no form for (SSA_FIRINGDEP_REFUSAL)

ssa_nrm_eligible(sn)

[OK, REASON] = SSA_NRM_ELIGIBLE(SN) Whether the next-reaction-method engine can run this model AT ALL.

THIS IS THE PREFERENCE QUESTION, not the gate question. SOLVER_SSA_ANALYZER asks it to decide whether the ‘default’ and ‘parallel’ paths should run the NRM instead of the serial engine. A caller asking whether an EXPLICIT options.method=’nrm’ will produce an answer must ask SSA_NRM_REFUSAL instead, which is narrower: the explicit arm FALLS BACK to the serial engine, with a warning, for every one of these conditions except the scheduling one.

The tests themselves live in SSA_NRM_GUARDS, which SOLVER_SSA_ANALYZER also reads field by field so its fallback warnings can name the construct. One body, three callers.

The NRM simulates open and closed models alike – it only lacks Fork/Join node handling – so the model-class exclusion is Fork/Join, not the INF/PS-only sn_is_population_model.

ssa_mlmc_feature_set()

SSA_MLMC_FEATURE_SET Population-network envelope; structural checks run in the adapter.

ssa_imm_tput_balance(sn, QN, RN, TN, CN, XN)

[TN, RN, CN, XN] = SSA_IMM_TPUT_BALANCE(SN, QN, RN, TN, CN, XN)

Throughput of the station-class pairs whose service is Immediate, taken from flow balance rather than from the sample path.

Such a station holds a job for zero time, so its departure rate IS its arrival rate. Both SSA engines estimate TN as the time integral of the departure rate (NRM: the propensity; serial: the empirical state distribution times the departure rates), which at an Immediate station is pure noise: each visit lasts ~Exp(Immediate) and contributes rate*dt ~ Exp(1) instead of one departure. On the Tutorial 6 cache model this left Client Tput 2.9944 against the ArvR 2.9214 its feeders imply. TN is therefore replaced by the arrival rate given by the same flow balance that reports ArvR (sn_get_arvr_from_tput: the timed stations’ integrals, the measured cache hit/miss split, routing through stateless nodes), sweeping until a chain of Immediate stations settles. The reference throughput XN, the response times RN of the replaced pairs and the cycle times CN follow. Skipped under JSQ/SQ routing, whose split is state dependent and is not in sn.rt.

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

ssa_fluid_refusal(sn, method)

[OK, REASON] = SSA_FLUID_REFUSAL(SN, METHOD) Whether SolverSSA can simulate this Fluid Stochastic Petri Net under METHOD.

A FLUID NET RUNS ON THE NRM PETRI RUNNER ONLY. Its levels move continuously between firings, so the simulator is a piecewise-deterministic Markov process: the runner keeps one clock per timed mode over the marking, and the levels, their guard crossings and their bound hits ride beside those clocks as scheduled non-firing events (SOLVER_SSA_NRM_SPN). The serial engine walks the State arms of the queueing stations and has no level to move, so ‘serial’ is refused by name and ‘default’/’parallel’ run the NRM. A net the NRM’s reaction builder cannot read (a non-exponential firing) has no engine left, and is refused with the reason the builder would give. A marking-dependent firing or flow rate is read by the runner, segment by segment.

Parameters:

• sn - NetworkStruct

• method - options.method

Returns:

ok - true when the net holds no fluid, or holds fluid and METHOD can run it reason - the refusal, empty when OK

See also SN_FLUID_REFUSAL, SSA_FIRINGDEP_REFUSAL, SSA_NRM_GUARDS.

ssa_firingdep_refusal(sn)

[OK, REASON] = SSA_FIRINGDEP_REFUSAL(SN) Whether the NRM engine of SolverSSA can answer a Petri net whose firing rates depend on the marking (Transition.setFiringRateDependence).

IT CANNOT, AND SAYS SO BY NAME. The NRM builds one constant-propensity reaction per timed mode, so the g(marking) multiplier has no place in it and answering with the nominal rate would be a wrong number under the caller’s own method name. The SERIAL engine does apply it: State.afterGlobalEvent scales the mode’s firing rate by firingdep{m}(ep_space) on the state it is leaving, the same code and the same marking matrix SolverCTMC uses, so the refusal is confined to the NRM (as in C++, ssa_check_firingdep).

ONE PREDICATE, FOUR CALLERS: SSA_NRM_GUARDS stores it as g.firingdep, so ‘default’ and ‘parallel’ do not prefer the NRM on such a net and run the serial engine; SSA_NRM_REFUSAL answers with it so model.help does not offer ssa.nrm; solver_ssa_analyzer raises with it under an explicit ‘nrm’, and the NRM reaction builder does on the enableChecks=false path.

A FLUID NET IS EXEMPT. Its runner (SPNFLUIDRUN) recomputes every propensity at each segment, and the marking is constant within a segment, so it applies g(marking) exactly, on the same input-place marking the serial engine reads.

solver_ssa_reachability(sn, options)

[SSQ,SSH,QN] = SOLVER_SSA_REACHABILITY(QN,OPTIONS)

solver_ssa_nrm_space_analyzer(sn, options)

SOLVER_SSA_ANALYZER_NRM Performance indices from SSA/NRM simulation This variant runs the next‑reaction–method SSA on the network SN and returns mean throughput (XN), utilisation (UN), queue length (QN), response time (RN), throughput per station (TN) and class residence time (CN). Results are based on the empirical state probabilities obtained from SOLVER_SSA_NRM. * Only stations whose scheduling policy is INF, EXT or PS are supported. Any other policy triggers an error. Cache nodes are not supported in this simplified variant. *

solver_ssa_nrm_space(sn, options)

SOLVER_SSA_NRM_SPACE Steady‑state analysis via the Next‑Reaction Method (SSA)

[PI, SSQ, ARVRATES, DEPRATES, SN] = SOLVER_SSA_NRM_SPACE(SN, OPTIONS) runs a stochastic simulation of the queueing network described in SN for OPTIONS.samples reaction firings using Gibson & Bruck’s Next‑Reaction Method. During the run it:

• observes every distinct global state visited and the time spent in it;

• accumulates the per‑state propensities of all enabled reactions, including self‑loops (service completions routed back to the same queue).

Outputs

PI – 1×S vector of empirical steady‑state probabilities

(sojourn‑time fractions) for the S unique states;

SSQ – S×(M·R) matrix listing those states row‑by‑row in the

flattened (station, class) order;

DEPRATES – S×(M·R) matrix of total departure rates from each queue

in the corresponding state;

SN – (Possibly updated) network structure.

See also NEXT_REACTION_METHOD.

solver_ssa_nrm_OMR(sn, options)

SOLVER_SSA_NRM Steady‑state analysis via the Next‑Reaction Method (SSA)

[PI, SSQ, ARVRATES, DEPRATES, SN] = SOLVER_SSA_NRM(SN, OPTIONS) runs a stochastic simulation of the queueing network described in SN for OPTIONS.samples reaction firings using Gibson & Bruck’s Next‑Reaction Method. During the run it:

• observes every distinct global state visited and the time spent in it;

• accumulates the per‑state propensities of all enabled reactions, including self‑loops (service completions routed back to the same queue).

Outputs

PI – 1×S vector of empirical steady‑state probabilities

(sojourn‑time fractions) for the S unique states;

SSQ – S×(M·R) matrix listing those states row‑by‑row in the

flattened (station, class) order;

DEPRATES – S×(M·R) matrix of total departure rates from each queue

in the corresponding state;

SN – (Possibly updated) network structure.

See also NEXT_REACTION_METHOD.

solver_ssa_nrm_OM2R2(sn, options)

SOLVER_SSA_NRM Steady‑state analysis via Gibson & Bruck’s Next‑Reaction Method (SSA) as modified in Anderson’s, THE JOURNAL OF CHEMICAL PHYSICS 127, 214107, 2007.

solver_ssa_gdfactor(sn, cur_state)

GDNOW = SOLVER_SSA_GDFACTOR(SN, CUR_STATE)

Evaluates the globally state-dependent rate scaling phi(n) declared through setGlobalDependence on the CURRENT sample-path state, returning the (nstations x nclasses) matrix of scalings.

SolverCTMC tabulates phi once per state of the enumerated space; a simulator instead has one state at a time, so the same factorization applies with the table collapsed to a single row. Within a state phi is a constant, so it multiplies the rate of every station service event at that state.

solver_ssa_findenabled(sn, node_a, enabled_next_states, cur_state, outprob_a, event_a, class_a, isSimulation, node_p, local, outprob_p, event_p, class_p, sync, gsync, depRatesSamples, samples_collected, arvRatesSamples, last_node_a, last_node_p, eventCache, startRatesSamples, preemptRatesSamples)

The START/PREEMPT accumulators and per-transition tags mirror the inlined enumeration in solver_ssa.m: this is the SAME enumeration, kept as a function, so tagging only one of the two would make the counters depend on the use_inline performance flag.

solver_ssa_analyzer_serial(sn, init_state, options, isHashed)

[XN,UN,QN,RN,TN,CN,TRANSYSSTATE,TRANSYNC,SN,STARTN,PREEMPTN,TRANSTARTTAG,TRANPREEMPTTAG]=SOLVER_SSA_ANALYZER_SERIAL(SN, OPTIONS)

STARTN and PREEMPTN are the (station x class) service-start and preemption rates, estimated exactly as TN is: the time average of the rate the enabled transitions carry. TRANSTARTTAG/TRANPREEMPTTAG are the per-step tags of the sampled path, used by @SolverSSA/sample.m to emit START and PREEMPT events.

solver_ssa_analyzer_parallel(sn, init_state, laboptions)

[XN,UN,QN,RN,TN,CN]=SOLVER_SSA_ANALYZER_PARALLEL(SN, INIT_STATE, LABOPTIONS)

Worker-count-invariant parallel SSA.

The simulation budget is split into a FIXED number of independent replications R (laboptions.config.nreplicas). Replication r simulates ceil(samples/R) events and is seeded deterministically with (seed+r-1): inside a parfor body spmdIndex resolves to 1 in solver_ssa, so the effective seed is exactly laboptions.seed+r-1, independent of which worker executes the iteration. The R replications are distributed over whatever workers the parallel pool provides (and run on the client if no pool exists), and the per-replication estimates are averaged.

Because replication r always uses the same seed and sample budget no matter how many workers are available, the returned averages depend only on (seed, samples, R) and are invariant to the worker count. This is the key difference from the previous spmd implementation, whose result varied with the pool size (it both divided the budget by, and seeded labs from, the runtime number of labs).

ssa_mlmc(problem, T, options)

SSA_MLMC Estimate a vector terminal reward of a Markov jump process. problem.x0: d-by-1 counts; problem.nu: d-by-K fixed increments; problem.rates(x): K-by-1 nonnegative rates, defined also on auxiliary leap states; problem.observe(x): vector reward, defined on those states. Exact-corrected split coupling: Anderson & Higham, arXiv:1107.2181. Pilot paths are discarded and production counts are fixed before sampling. Intervals are pointwise normal approximations, not certified error bounds.

ssa_mlmc_problem(model, rewardName)

SSA_MLMC_PROBLEM Compile exponential population networks into jump channels. Auxiliary leap populations are signed; rate laws use their positive part. No population is clipped, so every linear channel invariant is preserved. Rewards must therefore be defined on signed populations as well.

solver_ssa_analyzer(sn, options)

[QN,UN,RN,TN,CN,XN,RUNTIME,…,STARTN,PREEMPTN,TRANSTARTTAG,TRANPREEMPTTAG] = SOLVER_SSA_ANALYZER(SN, OPTIONS)

STARTN/PREEMPTN are the derived (station x class) service-start and preemption rates, and TRANSTARTTAG/TRANPREEMPTTAG the per-step tags of the sampled path. Both back-ends fill them, so the counters do not depend on which one the dispatch picks. Wrapper that selects the most suitable SSA performance-analysis back-end.

If every station uses scheduling policy INF, EXT, or PS (and the network has no cache nodes), the faster Next-Reaction-Method analyser

-> solver_ssa_analyzer_nrm

is invoked. Otherwise the original serial / parallel analysers are used.

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

solver_ssa_analyzer_nrm(sn, options)

SOLVER_SSA_ANALYZER_NRM Performance indices from SSA/NRM simulation This variant runs the next‑reaction–method SSA on the network SN and returns mean throughput (XN), utilisation (UN), queue length (QN), response time (RN), throughput per station (TN) and class residence time (CN). Results are based on the empirical state probabilities obtained from SOLVER_SSA_NRM. * Only the scheduling policies listed in ALLOWEDSCHED below are supported. Any other policy triggers an error. Cache nodes are not supported in this simplified variant. *

solver_ssa(sn, init_state, options, eventCache)

[PI,SSQ,ARVRATES,DEPRATES,TRANSYSSTATE,QN,DLYRATES,STARTRATES,PREEMPTRATES,TRANSTARTTAG,TRANPREEMPTTAG]=SOLVER_SSA(QN,OPTIONS)

STARTRATES and PREEMPTRATES are the derived START/PREEMPT rates per unique state, laid out like ARVRATES/DEPRATES: how fast the transitions enabled in that state start a class-r service at a stateful node, or push a class-r job in service there back into the buffer.

TRANSTARTTAG{e} and TRANPREEMPTTAG{e} are the tags of the transition that actually FIRED at step e, as rows [statefulIndex class]. The trace needs them explicitly: @SolverSSA/sample.m rebuilds events from sn.sync{tranSync(e)}, which by construction carries no derived tag.

class EventCache

EventCache Event caching utilities for SSA solver optimization

EventCache provides static methods to create and manage event caching for the Stochastic State-space Analysis (SSA) solver. Event caching improves performance by storing computed events to avoid redundant calculations during simulation.

@brief Event caching system for SSA solver performance optimization

The cache stores mappings of state transitions and their associated events, enabling faster simulation execution by reusing previously computed results.

Method Summary

static create(enabled, sn)

CREATE Create an event cache container

@brief Creates a dictionary for event caching if enabled @param enabled Boolean flag to enable/disable caching @param sn Network structure containing model information @return eventCache dictionary instance or empty array

solver_ssa_nrm(sn, options)

SOLVER_SSA_NRM Steady‑state analysis via the Next‑Reaction Method (SSA)

[QN, UN, RN, TN, CN, XN, LG, SN] = SOLVER_SSA_NRM(SN, OPTIONS) runs a stochastic simulation of the queueing network described in SN for OPTIONS.samples reaction firings using Gibson & Bruck’s Next‑Reaction Method. During the run it:

• computes performance metrics directly during simulation;

• returns standard queueing performance measures.

Outputs

QN – M×K matrix of mean queue lengths UN – M×K matrix of utilizations RN – M×K matrix of response times TN – M×K matrix of throughputs CN – 1×K vector of cycle times XN – 1×K vector of system throughputs LG – Logarithm of normalizing constant (not computed) SN – (Possibly updated) network structure. STARTN – M×K rate at which a class-r service starts at station i PREEMPTN – M×K rate at which a class-r job in service is displaced

See also SOLVER_SSA_NRM_SPACE, NEXT_REACTION_METHOD_DIRECT.