Product form of a stochastic Petri net: decide whether one exists and derive the per-level factors g_l that mdd_rec and spn_metrics take as input. More...

#include <algorithm>
#include <cmath>
#include <cstddef>
#include <limits>
#include <map>
#include <string>
#include <vector>
#include "line/api/mc/ctmc_solve.h"
#include "line/api/spn/spn_mdd.h"
#include "line/util/error.h"
#include "line/util/lstsq.h"
#include "line/util/matrix.h"

Include dependency graph for spn_pf.h:

Classes
struct	line::spn::SpnPfOptions
	Options of the product-form derivation. More...
struct	line::spn::SpnPfResult< T >
	The product form, and the certificate that it is one. More...

Namespaces
namespace	line
namespace	line::spn

Functions
template<class T>
SpnPfResult< T >	line::spn::spn_pf (const qn::NetworkStruct< T > &sn, const SpnPfOptions &options=SpnPfOptions())
	Derive the product form of a stochastic Petri net.

Detailed Description

Product form of a stochastic Petri net: decide whether one exists and derive 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 equals 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 tested for here:

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 through a marking-dependent firing rate 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.

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, reached in a form that is unique whichever least-squares primitive a codebase carries: solve (rows rows^T) w = rhs and set x = rows^T w. Any two solutions w of that system give the SAME rows^T w, so the answer does not depend on how the rank-deficient solve breaks its tie.

ARITHMETIC. The logarithm and the exponential are unavoidable here – y is the exponential of a least-squares solution – so the whole derivation needs a transcendental field and is refused under exact arithmetic by name. What consumes the g_l afterwards, mdd_rec and spn_metrics, stays rational.

References: 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

Definition in file spn_pf.h.

Classes

Namespaces

Functions

Detailed Description