Port of ode_eliminate_immediate.m, eliminate_immediate_matrix.m and ode_solve_stiff.m: the two answers to an IMMEDIATE transition. More...

#include <algorithm>
#include <cmath>
#include <cstddef>
#include <exception>
#include <functional>
#include <string>
#include <type_traits>
#include <vector>
#include "line/api/mc/dtmc_solve.h"
#include "line/lang/qn/network_struct.h"
#include "line/util/linalg.h"
#include "line/solvers/fluid/fluid_odes.h"
#include "line/util/error.h"
#include "line/util/lsoda.h"
#include "line/util/matrix.h"
#include "line/util/ode.h"

Include dependency graph for fluid_stiff.h:

Classes
struct	line::fluid::FluidImmediateResult
	What an elimination attempt produced. More...
struct	line::fluid::FluidImmediateMatrix< T >
	What the matrix-level elimination produced. More...
struct	line::fluid::FluidStiffOptions
	Controls for the stiff arm. More...

Namespaces
namespace	line
namespace	line::fluid

Functions
double	line::fluid::fluid_immediate_transition_tol ()
	The reference's two thresholds, which do NOT agree and are not meant to.
double	line::fluid::fluid_immediate_state_tol ()
template<class T, class Opt>
bool	line::fluid::fluid_hide_immediate (const qn::NetworkStruct< T > &sn, const Opt &opt)
	Stochastic complementation of the INSTANTANEOUS coordinates of a fluid drift, the twin of ode_eliminate_immediate.m.
FluidImmediateResult	line::fluid::fluid_eliminate_immediate (const FluidOdeSystem &sys, double imm_tol=fluid_immediate_transition_tol())
template<class T>
FluidImmediateResult	line::fluid::fluid_eliminate_immediate (const qn::NetworkStruct< T > &sn, const FluidOdeSystem &sys, double imm_tol=fluid_immediate_transition_tol())
	The same, at the reference's own signature, which carries sn.
template<class T>
FluidImmediateMatrix< T >	line::fluid::fluid_eliminate_immediate_matrix (const Matrix< T > &W, double imm_tol=fluid_immediate_state_tol())
	Port of eliminate_immediate_matrix.m: the same elimination on a generator that is already assembled.
OdeSolution< double >	line::fluid::fluid_ode_solve_stiff (const std::function< void(double, const double , double )> &f, double t0, double t1, const std::vector< double > &y0, const FluidStiffOptions &opt=FluidStiffOptions())
	Port of ode_solve_stiff.m.
LsodaOptions	line::fluid::fluid_lsoda (const LsodaOptions &lopt)
	The step budget LSODA gets on a fluid leg before the stiff arm takes over.
std::vector< double >	line::fluid::fluid_integrate_leg (const std::function< void(double, const double , double )> &f, double t0, double t1, const std::vector< double > &y0, const LsodaOptions &lopt)
	One integration leg, with the reference's retry on a failed solve.
LsodaSolution	line::fluid::fluid_integrate_grid (const std::function< void(double, const double , double )> &f, const std::vector< double > &y0, const std::vector< double > &grid, const LsodaOptions &lopt)
	The same retry over a whole output grid, for the callers that ask LSODA for a trajectory rather than an endpoint.

Detailed Description

Port of ode_eliminate_immediate.m, eliminate_immediate_matrix.m and ode_solve_stiff.m: the two answers to an IMMEDIATE transition.

WHERE THE STIFFNESS COMES FROM. An Immediate distribution fires at GlobalConstants.Immediate, 1e8, while the rest of the model runs at rates of order one. The drift then has a mode with time constant 1e-8 alongside modes with time constant 1, and an explicit step controller is pinned to the fastest one for the WHOLE integration, long after that mode has died: the run either crawls or goes unstable. That is stiffness, and it is a property of the equations rather than of the integrator.

THE TWO ANSWERS ARE NOT EQUALLY GOOD, AND THE REFERENCE PREFERS THE FIRST. Eliminating the immediate transitions removes the fast mode from the system, so what is left is not stiff at all. Falling back to a stiff integrator keeps the fast mode and pays an implicit solve to stay stable across it. The first is ALGEBRA and the second is NUMERICS, which is why the first is exact:

ELIMINATION IS EXACT. The stochastic complement of a generator over a retained set is the generator of the process WATCHED ONLY ON THAT SET, S = Q11 + Q12 (-Q22)^-1 Q21. Its stationary law is the original's, conditioned on the retained set and renormalized. It is an identity, not an approximation, and it holds at any rate ratio. The one place a gap appears is when the answer is compared against the UNREDUCED system: the immediate states hold stationary mass of order 1/Immediate, about 1e-8, which the reduced system does not carry and the conditioning divides out.

THE STIFF INTEGRATOR ONLY COPES. It keeps the fast mode and controls the error on it, so its answer carries the integrator's tolerance and nothing better.

A CAVEAT THE REFERENCE DOES NOT STATE. Reconstructing a generator from the jump/rate representation reads the drift as dx/dt = x W, which is what the fluid drift IS wherever the state-dependent factor g(x) is x itself: an infinite server, or any station holding fewer jobs than it has servers. At a SATURATED station g is min(n_i, S_i)/n_i times x and the reconstructed W is not the drift there, so the complement is exact for the linear part only. That is the reference's behaviour, reproduced rather than corrected, and it is why the elimination is applied to the transitions and not to the metrics.

WHAT WAS ALREADY HERE, AND IS NOT REBUILT. mc::ctmc_stochcomp is the complement, with the shared LU factorization the reference gets from backslash. ode_rosenbrock4 in util/ode.h is the stiff integrator: a four-stage L-stable Rosenbrock method, the family MATLAB's ode23s belongs to, with an embedded estimate for step control. And the fluid path's ordinary integrator is LSODA, which already switches itself from Adams to BDF when it detects stiffness; the point of this file is therefore not "a stiff solver exists" but that the stiffness can be removed before anyone integrates.

Definition in file fluid_stiff.h.

Classes

Namespaces

Functions

Detailed Description