solvers.CTMC

solver_ctmc_mdd_supports(sn, options)

[BOOL, REASON, P, KIND] = SOLVER_CTMC_MDD_SUPPORTS(SN, OPTIONS)

@brief Can the mdd decision-diagram method be asked for this model?

The model-shape gate of the mdd method, asked as a predicate rather than raised. SOLVER_CTMC_MDD_ANALYZER calls it before it builds anything, and SolverCTMC.supportsModelMethod calls it so that a CALLER (model.help, findSolver, SolverAUTO) sees the same verdict without paying for a run. One predicate with two callers is what stops the report and the analyzer from disagreeing about which models the method serves.

A STOCHASTIC PETRI NET IS EXEMPT from the shape rules: a Place model is read through SPN_MDD, which builds the reachable set and the Kronecker descriptor from the marking rather than from the (station,class) encoding, so neither the single-class rule nor the closed-population rule applies to it. The transient rule binds it all the same.

THE TWO DEEPER RULES ARE HERE TOO, and the analyzer reads their products off this predicate rather than recomputing them: P is the station-to-station routing chain of the single class, which must be stochastic (a completion must move the job to another station, so a Router or a leak refuses), and KIND is the local-state encoding the disciplines and service laws admit, ‘np’ (count, or count plus one non-preemptive phase) or ‘ps’ (per-phase counts, shared servers only). Both are decided from sn alone, so a caller is told about them at the same price as the class count. On a net, or on a refusal, P and KIND are empty.

@param sn NetworkStruct of the model @param options solver options (read for timespan); optional @return bool true when the mdd method may run @return reason the refusal, or ‘’ when BOOL is true @return P station-to-station routing matrix of the single class @return kind ‘np’ or ‘ps’, the local-state encoding to build

solver_ctmc_cftp_supports(sn, options)

[BOOL, REASON] = SOLVER_CTMC_CFTP_SUPPORTS(SN, OPTIONS)

@brief Can the cftp perfect sampler be asked for this model?

The model-class gate of the cftp method, asked as a predicate rather than raised. SOLVER_CTMC_CFTP calls it to refuse a model it cannot sample, and SolverCTMC.supportsModelMethod calls it so that a CALLER (model.help, findSolver, SolverAUTO) sees the same verdict before paying for a run. Keeping it in one place is the point: a second copy in the analyzer is how the report and the run drift into two different answers.

The sampler is exact only on the closed single-class product form its balance function encodes; anything else must be refused, not approximated. What the feature registry CAN name (open classes, non-product-form disciplines, non-exponential service, …) is also declared in SolverCTMC.getMethodFeatureSet. The structural rules the registry has no name for are the class count, the station count and the steady-state restriction. The rest is stated TWICE ON PURPOSE: SolverCTMC.supportsModelMethod asks this predicate BEFORE the feature gate, and the analyzer asks it INSTEAD of the feature gate, so a construct declared only in the featset would be sampled away by a run with options.enableChecks=false.

@param sn NetworkStruct of the model @param options solver options (read for timespan) @return bool true when the cftp sampler may run @return reason the refusal, or ‘’ when BOOL is true

solver_ctmc_transient_analyzer(sn, options)

[T,PIT,QNT,UNT,RNT,TNT,CNT,XNT,INFGEN,STATESPACE,STATESPACEAGGR,EVENTFILTRATION,RUNTIME,FNAME] = SOLVER_CTMC_TRANSIENT_ANALYZER(QN, OPTIONS)

solver_ctmc_avg_from_pi(sn, pivec, StateSpace, StateSpaceAggr, arvRates, depRates, options)

SOLVER_CTMC_AVG_FROM_PI Map a state distribution to mean performance metrics.

[QN,UN,RN,TN,CN,XN] = SOLVER_CTMC_AVG_FROM_PI(SN, PIVEC, STATESPACE,: STATESPACEAGGR, ARVRATES, DEPRATES)

Given an arbitrary probability vector PIVEC over the enumerated CTMC state space of SN (rows of STATESPACE / STATESPACEAGGR), returns the per-(station, class) mean queue length QN, utilization UN, response time RN, throughput TN, system response time CN and system throughput XN. The discipline-aware mapping is identical to the steady-state reduction performed by solver_ctmc_analyzer; it is factored here so that callers holding their own distribution (e.g. the SolverENV state-vector analyzer, which time-averages a transient distribution) can reuse it without re-solving for the stationary vector.

solver_ctmc_auxfilt_init(sn, Q)

DFILTAUX = SOLVER_CTMC_AUXFILT_INIT(SN, Q)

Allocate the derived START/PREEMPT filtrations: one all-zero sparse matrix of the shape of Q per (station, class). DFILTAUX.start{i,r}(s,ns) will hold the rate at which the transition s -> ns starts a class-r service at station i, and DFILTAUX.preempt{i,r}(s,ns) the rate at which it pushes a class-r job in service back into the buffer there.

solver_ctmc_gdfactor(sn, stateSpaceAggr, options)

GDFACTOR = SOLVER_CTMC_GDFACTOR(SN, STATESPACEAGGR, OPTIONS)

Tabulates the globally state-dependent rate scaling phi(n) declared through setGlobalDependence, one evaluation per CTMC state. Returns an (nstates x nstations*nclasses) matrix whose (s, i + (r-1)*M) entry is the scaling of class r at station i in state s, i.e. the column-major flattening of the (nstations x nclasses) matrix phi returns.

The handle is evaluated ONCE per state and never per transition: phi may be expensive (a bandwidth-sharing allocation solves a convex program per call), and within a state it is a constant that multiplies every rate at that state.

solver_ctmc_reward(sn, options)

SOLVER_CTMC_REWARD Compute steady-state and transient rewards

[STEADYSTATE, NAMES, STATESPACE, PI, V, T] = SOLVER_CTMC_REWARD(SN, OPTIONS)

Computes both steady-state expected rewards using the stationary distribution:

E[r] = sum_s pi(s) * r(s)

and transient value functions using uniformization with value iteration:

V^{k+1}(s) = r(s) + sum_{s’} P(s,s’) * V^k(s’)

INPUTS:

sn - NetworkStruct with .reward field populated options - Solver options with:

.rewardIterations - Number of iterations for transient (default 1000)

OUTPUTS:

steadyStateRewards - Vector of steady-state expected rewards [nRewards x 1] rewardNames - Cell array of reward names stateSpace - State space matrix [nStates x nDims] pi - Stationary distribution [1 x nStates] V - Cell array of value functions {nRewards x 1}

Each V{r} is [Tmax+1 x nStates] matrix

t - Time vector (scaled by uniformization rate)

solver_ctmc_mdd_analyzer(sn, options, model)

[QN,UN,RN,TN,CN,XN,INFO] = SOLVER_CTMC_MDD_ANALYZER(SN, OPTIONS) [QN,UN,RN,TN,CN,XN,INFO] = SOLVER_CTMC_MDD_ANALYZER(SN, OPTIONS, MODEL) Stationary analysis of a closed single-class network, or of a stochastic Petri net, whose CTMC state space is held in a decision diagram and solved by level 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.

This is the ‘mdd’ method of SolverCTMC. It never forms the |S|-state generator: the reachable set is stored in an MDD and K coupled level-CTMCs are iterated to a fixed point, so the memory cost is O(sum_k |M_k|) rather than O(|S|). The saving grows with the number of stations, and is negative at K=3, where the diagram compresses nothing.

– Exactness The single approximation is Pr{i_k | alpha} = Pr{i_k | p}. It is EXACT on product-form networks (paper Sec. 5), which covers exponential service under any work-conserving discipline and general service at PS or IS stations (BCMP types 2 and 3). It is an approximation otherwise, notably phase-type service at FCFS or LCFS, where errors of a fraction of a percent on the mean queue lengths have been observed.

– Petri nets Passing MODEL routes a net holding Places and Transitions through SPN_MDD instead of MDD_DESCRIPTOR: the levels are then (place, class) pairs plus one phase level per phase-type mode, and the measures come back per place. The approximation is the same Eq. 5 as for a queueing network, and it is exact on a product-form net, which SolverNC’s ‘rec’ method (SOLVER_NC_SPN_ANALYZER) solves exactly and far more cheaply – the aggregation earns its place on the nets that have NO product form.

A net carries no per-station service rate, so MDD_MCD returns only the level marginals. The token throughput is then assembled here from the mode rates and those marginals, under the SAME independence across levels that the aggregation already assumes: it is the method’s own approximation applied once more, not a second one layered on top.

– Output QN,UN,RN,TN : per station, mean queue length, utilization, response time and

throughput

CN,XN : per class, system response time and throughput INFO : struct with the diagram, the descriptor, the level sizes and

the iteration count

See also: mdd_mcd(), mdd_descriptor(), mdd_ps(), mdd_reachset(), spn_mdd(), SolverCTMC.

solver_ctmc_margaggr(sn, options)

[PNIR,PI,RUNTIME,FNAME] = SOLVER_CTMC_MARGAGGR(QN, OPTIONS)

solver_ctmc_marg(sn, options)

[PNIR,RUNTIME,FNAME] = SOLVER_CTMC_MARG(QN, OPTIONS)

solver_ctmc_jointaggr(sn, options)

[PNIR,RUNTIME,FNAME] = SOLVER_CTMC_JOINTAGGR(QN, OPTIONS)

solver_ctmc_joint(sn, options)

[PNIR,RUNTIME,FNAME] = SOLVER_CTMC_JOINT(QN, OPTIONS)

solver_ctmc_chain_transient(chain, pi0, timespan)

[PI_T, T] = SOLVER_CTMC_CHAIN_TRANSIENT(CHAIN, PI0, TIMESPAN)

Transient distribution of a user-supplied Markov chain over TIMESPAN = [t0,t1]. For a MarkovProcess the Kolmogorov forward equations are integrated from PI0; for a MarkovChain the distribution is advanced one step per unit of time, so the returned T holds the integer steps within TIMESPAN.

solver_ctmc_chain(chain, options)

[PI, INFGEN, STATESPACE, RUNTIME] = SOLVER_CTMC_CHAIN(CHAIN, OPTIONS)

Steady-state analysis of a user-supplied Markov chain, i.e. a MarkovProcess (CTMC, generator Q) or a MarkovChain (DTMC, transition matrix P). INFGEN is Q for a CTMC and the uniformized generator P-I for a DTMC, which carries the same stationary vector. STATESPACE is the state space attached to the chain, or the state indices when the chain carries none.

ctmc_stationary(Q, StateSpace, sn, options)

PI = CTMC_STATIONARY(Q, STATESPACE, SN, OPTIONS)

Single entry point for the stationary distribution of a CTMC generated from a NetworkStruct.

All the stationary mass of a reducible chain lives in its bottom strongly connected components, each weighted by the probability of being absorbed in it from the declared initial state; every other state is transient and carries zero. The block decomposition handles the irreducible case as the degenerate one BSCC / no transient states, so every CTMC solve goes through it and no dispatch can disagree with the algorithm about whether a chain is reducible.

See _kb/11-conventions-and-gotchas.md.

solver_ctmc_cftp(sn, options)

Draws iid states from the exact stationary distribution with pfqn_cftp and reduces them to the standard mean performance metrics.

Parameters:

sn – Network structure.
options – Solver options (method, samples, seed).

Returns:

QN – Mean queue length per station and class. UN: Utilization per station and class. RN: Response time per station and class. TN: Throughput per station and class. CN: System response time per class. XN: System throughput per class. Xs: Sampled states, one per row (samples x stations). Ts: Per-sample coalescence horizon or mixing steps. pAggr: Empirical probability of each distinct sampled state. SSq: Distinct sampled states, aligned with pAggr.

[QN,UN,RN,TN,CN,XN,XS,TS,PAGGR,SSQ] = SOLVER_CTMC_CFTP(SN, OPTIONS)

Perfect-sampling steady-state analysis of a closed single-class product-form network. States are drawn iid from the exact stationary distribution by monotone Coupling From The Past, so the estimator carries Monte Carlo error O(samples^(-1/2)) but never enumerates the state space.

Reference: S. Kijima and T. Matsui, “Approximate/Perfect Samplers for Closed Jackson Networks”, Proc. Winter Simulation Conference, 2005.

solver_ctmc_ratecomplement(D, nonimm, imm, Q12, Q22)

R = SOLVER_CTMC_RATECOMPLEMENT(D, NONIMM, IMM, Q12, Q22)

Long-run rate of an action, as seen from each tangible (non-vanishing) state, given the action’s rate filter D over the full state space.

Vanishing states are removed from the generator by stochastic complementation, so an action that fires only in vanishing states (a fork firing, a join departure, or the firing of an immediate SPN mode) would be lost if its rate were read off the tangible rows alone. The rate observed from tangible state s is the direct exit rate via the action plus the expected number of firings along the vanishing chain entered from s:

r = D(nonimm,:)*1 + Q12*(-Q22)^(-1)*(D(imm,:)*1)

where Q12 and Q22 are the tangible-to-vanishing and vanishing-to-vanishing blocks returned by ctmc_stochcomp.

ctmc_signal_lossy(sn, arvRates, probSysState, wset, isf)

LOSSY = CTMC_SIGNAL_LOSSY(SN, ARVRATES, PROBSYSSTATE, WSET, ISF)

Classes whose jobs can be annihilated by a G-network signal at the stateful node ISF, as a 1 x nclasses logical row. Such a job leaves the station without a service completion, so the arrival-based utilization estimator (offered load) is invalid there and only the departure-based (carried load) estimator UN = T * E[S] / c is meaningful.

A signal class R is active at ISF when its stationary arrival rate there is positive. A targeted signal (SN.SIGNALTARGET(R) >= 1) only removes that class; an untargeted one is class-agnostic and removes any non-signal class, matching State.afterEventStationSignal, MAM and LDES.

ctmc_signal_busy(sn, ind, ist, schedIst, SIst, StateSpace, istSpaceShift, wset, probSysState)

UNB = CTMC_SIGNAL_BUSY(SN, IND, IST, SCHEDIST, SIST, STATESPACE,: ISTSPACESHIFT, WSET, PROBSYSSTATE)

Exact per-class busy-server fraction at station IST, read off the enumerated state space as a 1 x nclasses row.

For a class that a G-network signal can annihilate, the departure-based estimator T*E[S]/c is exact only under exponential service: a job destroyed mid-service leaves behind busy time with no completion, so with phase-type service T*E[S]/c under-counts (M/Er2/1 with lambda+=0.5, lambda-=0.4 gives 0.34941 against a true 0.37696). The in-service occupancy below is exact for any service process.

PS-like disciplines share the servers among all resident jobs, so class k gets the weighted share n_k w_k / sum_j n_j w_j of the busy servers; the remaining disciplines expose the in-service indicator directly through State.toMarginal.

solver_ctmc_state_supports(sn, solverName)

[BOOL, REASON] = SOLVER_CTMC_STATE_SUPPORTS(SN, SOLVERNAME)

@brief Can the State machinery (lang/+State) serve the impatience laws: this model declares?

The rules of State.afterEventStation that a feature name cannot state, asked as a predicate rather than raised. The explicit-generator methods of SolverCTMC and BOTH engines of SolverSSA drive the same afterEvent code (the serial engine directly, the NRM through its own propensities, which read sn.retrialMu as a memoryless rate), so each solver used to carry a copy of these tests at the top of its analyzer (solver_ctmc.m, solver_ssa.m) and nothing above them asked: model.help offered ctmc.default and every ssa.* row on a model with a phase-type patience, and each then raised. ONE BODY, FOUR CALLERS: the two analyzers, which raise, and the two supportsModelMethod gates, which answer.

What it refuses, and why the registry cannot: Reneging, Balking and Retrial are DECLARED by both solvers, and the rules are about the LAW behind them, i.e. reneging with a non-exponential patience, balking outside QUEUE_LENGTH, a non-exponential retrial delay, a finite maxAttempts, a retrial station serving more than one class.

The multi-server DPS/GPS station afterEventStation also refuses is NOT a rule here: Queue.setNumServers refuses it at construction, so no built model carries one and the analyzer’s own raise is the only guard it needs.

@param sn NetworkStruct of the model @param solverName ‘SolverCTMC’ or ‘SolverSSA’, named in the refusal @return bool true when the model may run @return reason the refusal, or ‘’ when BOOL is true

solver_ctmc_analyzer(sn, options)

[QN,UN,RN,TN,CN,XN,INFGEN,STATESPACE,STATESPACEAGGR,EVENTFILTRATION,RUNTIME,FNAME,sn,AUXFILTRATION,STARTN,PREEMPTN] = SOLVER_CTMC_ANALYZER(sn, OPTIONS)

AUXFILTRATION carries the derived START/PREEMPT filtrations (see solver_ctmc); STARTN and PREEMPTN are the corresponding (station x class) rates pi*F*e: how often per unit time a class-r service starts at station i, and how often a class-r job in service is pushed back into the buffer there. For a lossless station with no in-service abandonment they satisfy STARTN == TN + PREEMPTN, which is the identity the tags exist to expose.

solver_ctmc_fcr_waitq(sn, options)

[SS,SSA,SSH,DFILT,SN,BASBLOCKQ,DFILTAUX]=SOLVER_CTMC_FCR_WAITQ(SN,OPTIONS)

DFILTAUX carries the derived START/PREEMPT filtrations, as in solver_ctmc. Here a single transition can start several services: the FIFO release cascade admits blocked jobs one after another, and each admission may take a server. The tags therefore accumulate along the cascade path and are written once, at the settled state, weighted by the same rate as the transition. Reachability-based state space and per-action rate filters for models with a finite capacity region (FCR) whose drop rule is WAITQ (waiting queue).

JMT WAITQ semantics (reference, mirrored by LDES): a job refused entry to a full region leaves the upstream station and waits in a per-region FIFO of (class, destination) tokens outside the region; after every transition that frees region capacity, tokens are released strictly in FIFO order (head-of- line: a stuck head blocks the queue) as long as the admission constraints (global cap, per-class caps, memory budget, linear constraints A*x<=b) permit; a fresh arrival that satisfies the constraints is admitted even if the FIFO is non-empty (it overtakes a head stuck on a different constraint). Blocked jobs are counted neither in the region occupancy nor in any station state, so station QLen excludes them, matching the JMT report convention.

True-BAS blocking between two stations is orthogonal to the region rule and is handled here as it is in the default generator: when a service completion at a BAS blocking station finds its destination unable to admit the job, the job is HELD AT THE SERVER (blocked marker set) rather than the departure being voided. Voiding it instead frees the server to re-serve the same job, which for exponential service is NOT equivalent to BAS: on release the held job enters the destination immediately, whereas a re-serving server must first draw a fresh completion, and throughput is understated. Those become-blocked arcs change the chain but are not departures, so they are accumulated separately in BASBLOCKQ and never enter DFILT.

The CTMC state is augmented as [h(1:nstateful), buf_1, …, buf_F] where h are the per-node hashed states and buf_f is the token FIFO of region f, padded with zeros to its maximum length. Classes whose region rule is DROP keep the transition-censoring behavior of the default generator (exact for memoryless sources, cross-validated against JMT).

solver_ctmc_auxrates(DfiltAux, probSysState, wset, M, K)

[STARTN, PREEMPTN] = SOLVER_CTMC_AUXRATES(DFILTAUX, PROBSYSSTATE, WSET, M, K)

Reduce the derived START/PREEMPT filtrations to (station x class) rates: StartN(i,r) = pi * F_start{i,r} * e, the long-run number of class-r service starts per unit time at station i, and likewise for preemptions. Summing the row of the filtration and weighting by the stationary probability is the same reduction the departure rates use, so the two are directly comparable: StartN == TN + PreemptN at a lossless station.

solver_ctmc_auxfilt_add(DfiltAux, sn, node, s, ns, w, starttag, preempttag, irow)

DFILTAUX = SOLVER_CTMC_AUXFILT_ADD(DFILTAUX, SN, NODE, S, NS, W, STARTTAG, PREEMPTTAG, IROW)

Accumulate the START/PREEMPT annotation of one successor row into the derived filtrations of the station behind NODE. W is the same weight the caller added to Dfilt{a}(S,NS), so the filtration integrates rate * count and pi*F*e is a rate of starts (or of preemptions) per unit time.

IROW selects the annotated successor row; a node that is not a station, or an event that tagged nothing, contributes nothing.

solver_ctmc(sn, options)

[Q,SS,SSQ,DFILT,ARVRATES,DEPRATES,QN,DFILTAUX]=SOLVER_CTMC(QN,OPTIONS)

DFILTAUX carries the two DERIVED event filtrations, DFILTAUX.start{i,r} and DFILTAUX.preempt{i,r}: the (state x state) rate at which a class-r job begins holding a server at station i, and the rate at which one is pushed back into the buffer there. They are deliberately NOT appended to DFILT: a START rides on the SAME arc as the ARV or DEP that causes it, so adding it to the eventFilt cell would double-count D1 in @SolverCTMC/sample.m (D0 = infGen - sum(eventFilt)) and break the pairing with sn.sync that getGenerator.m asserts. See _kb/06-solver-catalog.md.