Stationary density and distribution of a first- or second-order level-dependent (multi-regime) Markovian fluid queue, evaluated at requested fluid levels. More...

#include <cmath>
#include <complex>
#include <cstddef>
#include <vector>
#include "line/api/mam/mfq_ld_mean.h"
#include "line/num/number.h"
#include "line/util/eig.h"
#include "line/util/error.h"
#include "line/util/expm.h"
#include "line/util/linalg.h"
#include "line/util/lu.h"
#include "line/util/matrix.h"

Include dependency graph for mfq_ld_distr.h:

Namespaces
namespace	line
namespace	line::mam

Enumerations
enum class	line::mam::FluidDistrKind { line::mam::Pdf , line::mam::Pdfd , line::mam::Cdf , line::mam::Cdfm }
	Which functional of the stationary level law to evaluate. More...

Functions
template<class T>
std::vector< std::vector< T > >	line::mam::mfq_ld_distr (const LevelDependentFluidBlocks< T > &b, FluidDistrKind what, const std::vector< T > &points)
	Stationary density or distribution of a level-dependent fluid queue.

Detailed Description

Stationary density and distribution of a first- or second-order level-dependent (multi-regime) Markovian fluid queue, evaluated at requested fluid levels.

Port of matlab/src/api/mam/mfq_ld_distr.m and the BUTools LevelDependentFluidStationaryDistr it wraps. Its inputs are the matrix-exponential building blocks that mfq_ld_solve returns, the same LevelDependentFluidBlocks that mfq_ld_mean consumes, so the two are interchangeable consumers of one solve.

THE DENSITY. Over regime k, that is over the level interval [T(k), T(k+1)], the stationary density is anchored at both ends,

pi_k(x) = iniF_k exp(KF_k (x - T(k))) cloF_k
        + iniB_k exp(KB_k (T(k+1) - x)) cloB_k,

with point masses at the K+1 thresholds. Four quantities are offered: Pdf the per-state density Pdfd its derivative Cdf P(X < p), so a point mass sitting exactly at p is EXCLUDED Cdfm P(X <= p), so that mass is INCLUDED The Cdf/Cdfm distinction is not cosmetic here: a level-dependent fluid queue puts genuine atoms at its thresholds, so the two differ by a finite amount at every threshold and by nothing anywhere else.

INTEGRATING A SINGULAR EXPONENT. The cumulative forms need int_0^L exp(M u) du. When M is non-singular that is (-M)^-1 (I - exp(M L)); when M has a zero eigenvalue – which is the normal case for one of the two directions, not an edge case – the inverse does not exist and the reference deflates it instead. With l and r the left and right null vectors of M normalized so that l r = 1,

int_0^L exp(M u) du = (-(M - r l))^-1 (I - exp((M - r l) L))
                    + r l (L + exp(-L) - 1),

the rank-one shift moving the zero eigenvalue to -1 and its contribution being added back in closed form. The null vectors come from BUTools CRPSolve, which is a LINEAR SOLVE and not an eigenvector computation: it replaces the first column of M by ones and solves, which is legitimate because M has zero row sums by construction. That is worth stating because it is the reason this function does NOT need an eigendecomposition for its arithmetic.

WHERE THE EIGENVALUES DO ENTER, AND WHY Real IS STILL HONEST. The reference chooses which of KF, KB to deflate by comparing min |eig(KF)| against min |eig(KB)|, i.e. by asking which is closer to singular. That comparison is the ONLY use of an eigendecomposition in the whole function, and its result is a BRANCH SELECTION, a discrete choice between two formulas for the same integral. No eigenvalue flows into any returned number. The port therefore converts the two matrices to double for the comparison alone, exactly as util/eig.h instructs its callers to do, and carries out the integral itself in T. At Real50 the returned values are genuinely Real50-accurate; what is computed in double is which of two algebraically equivalent routes to take. If the build has no LAPACK, eig_values refuses by name, and this function refuses with it rather than guessing the branch.

ARITHMETIC. Gated on num_traits<T>::has_transcendental: expm is a tolerance-controlled Pade approximation in any arithmetic, and the deflation formula evaluates exp(-L). See the paragraph above for why the LAPACK dependency does not reduce this to double.

Definition in file mfq_ld_distr.h.

Namespaces

Enumerations

Functions

Detailed Description