solvers.AG

ag_inf_supports(sn)

[BOOL, REASON] = AG_INF_SUPPORTS(SN)

@brief Can the RCAT decomposition serve the infinite-server stations of: this model?

EVERY RCAT COMPONENT IS A SINGLE-SERVER QBD: build_rcat (solver_ag) gives each (station, class) pair a level-INDEPENDENT completion block kron(I, D1^s) and never reads sn.sched or sn.nservers, so a Delay, or any station scheduled INF, is decomposed as an M/PH/1 queue. That is exact only when the component never holds more than one job, i.e. when the class is closed with a population of one: then the level is 0 or 1 and a single-server chain IS the infinite-server chain. With two jobs or more, or with an open class, the component serves at mu where the station serves at n*mu, and the answer is that of a different model.

A self-looping class is exempt: it never leaves its reference station, and solver_ag reports it from sn.njobs directly (QN = N, UN = N at an INF station), which is exact.

ONE PREDICATE, TWO CALLERS: SolverAG.supportsModelMethod answers with it and solver_ag raises with it, so model.help cannot offer what the run refuses. The registry cannot state the rule because ‘Delay’ and ‘SchedStrategy_INF’ ARE declared; what decides is the population behind them.

@param sn NetworkStruct of the model @return bool true when every INF station may be decomposed @return reason the refusal, or ‘’ when BOOL is true

solver_ag(sn, options)

SOLVER_AG RCAT agent-based methods for SolverAG

[QN, UN, RN, TN, CN, XN, ITER] = SOLVER_AG(SN, OPTIONS)

Uses RCAT (Reversed Compound Agent Theorem) to find product-form solutions for queueing networks.

Each (station, class) pair becomes an isolated component, and the components are coupled only through the reversed rates of the synchronizing actions. A component is a QBD whose LEVEL is the queue length and whose PHASE is the pair (arrival phase, service phase), laid out in the Kronecker order of qbd_mapmap1: an arrival moves the level up carrying kron(D1^a, I), a service completion moves it down carrying kron(I, D1^s), the busy levels evolve under krons(D0^a, D0^s) and level zero under kron(D0^a, I), because no server is running there. With exponential processes every block is 1x1 and the QBD collapses to the scalar birth-death chain this analyzer used before, entry for entry.

'inap' - Iterative Numerical Approximation Procedure (default, fast)

'inapplus' - Improved INAP with weighted rates (no normalization)

'inapinf' - INAP with matrix-geometric solution of the isolated open: components (no state-space truncation), per Marin, Rota Bulo, Balsamo, “A Numerical Algorithm for the Decomposition of Cooperating Structured Markov Processes”, MASCOTS 2012.

'exact' - Not available (see solver_ag_autocat, unreachable)

ag_solve_component(Qk, mph, nlev, lvl)

AG_SOLVE_COMPONENT Stationary vector of one isolated component.

PI = AG_SOLVE_COMPONENT(QK, MPH, NLEV, LVL)

Shared by every execution backend and by the cluster fallback, so a remote agent and a local one are solved by the same code.

A component with a single phase per level is the birth-death chain the analyzer has always built, and the ratio recursion is both exact and stable there; a phase-expanded component is block tridiagonal instead, and the matrix analogue of that recursion (linear level reduction) keeps the same stability at the 100-level truncation, where a null-space solve is already ill-conditioned. Anything that reaches beyond the neighbouring level – a catastrophe, a batch removal – is neither, and falls back to ctmc_solve.

ag_is_tridiagonal(Q): IS_TRIDIAGONAL Check if a matrix is tridiagonal

ag_is_block_tridiagonal(Q, lvl): IS_BLOCK_TRIDIAGONAL True when every transition of Q stays within the neighbouring level, LVL being the level index of each state.

ag_birth_death_solve(Q)

BIRTH_DEATH_SOLVE Solve equilibrium of a birth-death (tridiagonal) CTMC

For a birth-death chain with birth rate lambda_n = Q(n, n+1) and death rate mu_n = Q(n, n-1), the equilibrium is computed using the recursion pi(n) = pi(n-1) * lambda(n-1) / mu(n).

This is numerically stable and avoids the ill-conditioned linear system that plagues null-space methods for large state spaces.

solver_ag_autocat(R, AP, varargin)

SOLVER_AG_AUTOCAT Native MATLAB implementation of autocat

[X, PI, Q, STATS] = SOLVER_AG_AUTOCAT(R, AP, OPTIONS)

Searches for RCAT product-form solutions using various relaxation methods. Native MATLAB implementation without YALMIP dependency.

Input:

R - Cell array of rate matrices AP - Action-Process mapping (A x 2) OPTIONS (name-value pairs):

‘relaxation’ - Relaxation method (default: ‘auto’)
‘lpr’ - LP relaxation with McCormick envelopes ‘tlpr0inf’ - Tightened LP with 0-level + infinity ‘tlpr1inf’ - Tightened LP with 1-level + infinity ‘tlprUinf’ - Tightened LP with U-level + infinity ‘ens’ - Exact nonlinear system (fmincon) ‘qcp’ - Quadratically constrained (fmincon) ‘ma’ - Mean Approximation (first moment matching) ‘va’ - Variance Approximation (slack variables) ‘minres’ - Minimum Residual (L1-norm of RC3) ‘inap’ - Iterative fixed-point approximation ‘tma1inf’ - Mean Approx + Tightened LP with 1-level + infinity ‘zpr’ - Zero potential relaxation with cutting planes ‘tzpr0inf’ - Tightened zero potential (0-level) ‘tzpr1inf’ - Tightened zero potential (1-level) ‘auto’ - Automatic selection

‘policy’ - Bound selection policy (default: ‘vol’)
‘vol’, ‘lu’, ‘l’, ‘u’

‘maxiter’ - Maximum iterations (default: 100) ‘tol’ - Convergence tolerance (default: 1e-4) ‘verbose’ - Verbosity level 0/1/2 (default: 1)

Output:

x - Action rate parameters pi - Equilibrium distributions Q - Generator matrices stats - Statistics struct

solver_ag_analyzer(sn, options)

[QN,UN,RN,TN,CN,XN,RUNTIME,METHOD,TOTITER,PERCRESULTS] = SOLVER_AG_ANALYZER(SN, OPTIONS)

Analyzer for the agent-based (RCAT) solver. Every (station, class) pair is an isolated component and the components are coupled only through the reversed rates of the synchronizing actions, so the whole analysis is the fixed point over that scalar vector; see solver_ag.

ag_stat_vector(C)

Stationary vector of a generator C, allowing a reducible one.

Replaces the first balance equation by the normalization, which is the equation it is redundant with (the columns of a generator sum to zero), and solves the resulting square system. Unlike a null-space solve this stays well posed when the chain is reducible with ONE closed class, which the level-0 chain of a phase-expanded component routinely is: a phase-type restarts in the support of alpha, so every service phase outside that support is unreachable once the queue has emptied at least once.

ag_qbd_finite_solve(Q, m, nlev)

QBD_FINITE_SOLVE Stationary vector of a finite block-tridiagonal generator.

Linear level reduction: censor the chain level by level from the top,: C(nlev-1) = B(nlev-1), C(n) = B(n) + F(n) * (-C(n+1))^-1 * D(n+1),

with B, F and D the diagonal, up and down blocks. C(0) is the generator of the chain censored on level 0, so pi_0 is its stationary vector and the rest follows from pi_(n+1) = pi_n F(n) (-C(n+1))^-1. This is the block form of BIRTH_DEATH_SOLVE and reduces to it entry for entry when m == 1.

ag_exec_resolve(options)

EXEC = AG_EXEC_RESOLVE(OPTIONS)

Resolve the execution backend of the RCAT fixed point into a struct the sweep can act on without re-reading options.

THE BACKEND CHANGES WHO EVALUATES AN AGENT, NEVER WHAT IT EVALUATES TO. Agent k’s generator is Q_k(x) = L_k + sum_{c passive at k} x_c Pb_c, which depends on the other agents only through the scalar reversed rates x. The sweep is Jacobi – every x_a is read off the PREVIOUS sweep’s stationary vectors, then all agents re-solve – so the agent order is immaterial and a parallel or distributed sweep produces the SAME iterates as the serial one.

HOW FAR THAT SURVIVES FLOATING POINT DEPENDS ON WHO RUNS THE AGENT SOLVE. ‘parallel’ is bit-identical to ‘serial’: the same code in the same process, differing only in an order that does not matter. ‘cluster’ is bit-identical only when the worker runs the same implementation as the coordinator – the wire is exact, since JSON round-trips a double without loss, but the stationary vector comes back from the WORKER’s solve, so a MATLAB coordinator driving a Java ag-worker agrees to a few ulp rather than bit for bit. That is the ordinary cross-codebase difference, not a protocol defect, but a cluster run must not regenerate a golden that a same-language run will later read.

Fields:: mode ‘serial’ | ‘parallel’ (alias ‘para’) | ‘cluster’ nworkers ‘parallel’: pool size, 0 = MATLAB’s default pool endpoints ‘cluster’: cell array of ‘host:port’ worker addresses timeout ‘cluster’: seconds to wait on a worker before solving its

agents locally instead

ag_cluster_sweep(x, Aa, Pb, L, ACT, PSV, numProcesses, A, N, meta)

[PI, Q] = AG_CLUSTER_SWEEP(X, AA, PB, L, ACT, PSV, NUMPROCESSES, A, N, META)

One sweep of the reversed-rate fixed point with the agents solved on remote ag-worker processes.

THE ONLY THING THAT CROSSES THE WIRE PER SWEEP IS X, one double per action, and back come the agents’ stationary vectors. The generator itself is rebuilt here rather than shipped, because assembling it is the cheap half (O(N^2)) and solving it is the expensive half (O(N^3)); see ag_agent_generator.

A worker that is missing, slow or broken is not fatal: its agents are solved locally through the same agent path, so the result is the run’s result either way and only the wall clock changes. That is a property of the decomposition, not a fallback bolted on – an agent depends on the rest of the model only through x, so anyone holding x can solve it.

ag_cluster_session(exec, agents, reset)

SESSION = AG_CLUSTER_SESSION(EXEC, AGENTS, RESET)

Open (or reuse) the connections to the ag-worker processes and assign the agents to them once, up front.

WHY THE ASSIGNMENT IS SEPARATE FROM THE SWEEP. An agent’s static data – its local rate matrix and the passive/active matrices of the actions it takes part in – does not change across the fixed point; only the reversed rates x do, and those are one double per action. Shipping the matrices once and then exchanging x per sweep is what makes remote execution worth doing; shipping them every sweep would put an O(N^2) payload on the wire to save an O(N^3) solve and lose most of the benefit.

The partition is round-robin over the worker list in agent index order, so it is a pure function of (number of agents, number of workers): a rerun assigns the same agents to the same workers, which is what keeps a distributed run reproducible.

A worker that cannot be reached is NOT fatal. Its slot is left empty and ag_cluster_sweep solves those agents locally, because any agent can be solved anywhere given x. A cluster run therefore degrades to a slower run, never to a wrong one or a failed one.

ag_cluster_agents(Aa, Pb, L, ACT, PSV, numProcesses, A, N, meta)

AGENTS = AG_CLUSTER_AGENTS(AA, PB, L, ACT, PSV, NUMPROCESSES, A, N, META)

The STATIC description of every agent, i.e. everything a worker needs that does not change across the fixed point. Only the reversed rates x change, so this is built once per solve and shipped once per worker.

Matrices travel as [row, col, value] triplets with 1-based indices, because the passive and active matrices of an action are nearly empty (an action touches one level transition) and a dense N-by-N payload would be mostly zeros. The receiving worker rebuilds the same matrices and therefore the same generator, which is what makes a remote agent’s answer identical to a local one rather than merely close.

ag_blk(n, m): Row/column range of level N (0-based) in a component with M phases per level.

ag_agent_generator(k, x, Aa, Pb, L, ACT, PSV, A, N)

AG_AGENT_GENERATOR Agent k’s generator at the current reversed rates.

QK = AG_AGENT_GENERATOR(K, X, AA, PB, L, ACT, PSV, A, N)

SPLIT OUT FROM THE SOLVE BECAUSE THE TWO HALVES COST DIFFERENT ORDERS. Assembling the generator is O(N^2) and solving it is O(N^3), so the ‘cluster’ backend ships only the stationary vector back and rebuilds the generator on the coordinator: sending an N-by-N matrix per agent per sweep to save the cheaper half would spend more on the wire than it saves on the worker.