First- and second-order level-dependent (multi-regime) Markovian fluid queues: the matrix-exponential building blocks of the stationary law. More...

#include <cmath>
#include <cstddef>
#include <vector>
#include "line/api/mam/mfq_ld_distr.h"
#include "line/api/mam/mfq_ld_mean.h"
#include "line/api/mam/mfq_solve.h"
#include "line/api/mam/qbd_r.h"
#include "line/num/number.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_solve.h:

Namespaces
namespace	line
namespace	line::mam

Enumerations
enum class	line::mam::FluidBoundary { line::mam::Reflective = 0 , line::mam::Absorbing = 1 }
	Boundary behaviour of one background state at a reflecting level. More...

Functions
template<class T>
LevelDependentFluidBlocks< T >	line::mam::mfq_ld_solve (const std::vector< Matrix< T > > &Q, const std::vector< Matrix< T > > &R, const std::vector< Matrix< T > > &S, const std::vector< T > &Thr, const std::vector< FluidBoundary > &boundaryL, const std::vector< FluidBoundary > &boundaryU, const std::vector< Matrix< T > > &Qt, const T &prec)
	Solve a first- or second-order level-dependent fluid queue.
template<class T>
LevelDependentFluidBlocks< T >	line::mam::mfq_ld_solve (const std::vector< Matrix< T > > &Q, const std::vector< Matrix< T > > &R, const std::vector< Matrix< T > > &S, const std::vector< T > &Thr)
	mfq_ld_solve with reflective boundaries, Qt = Q and the default prec = 1e-14.

Detailed Description

First- and second-order level-dependent (multi-regime) Markovian fluid queues: the matrix-exponential building blocks of the stationary law.

Port of matlab/src/api/mam/mfq_ld_solve.m and the BUTools SecondOrderLevelDependentFluidSolve it wraps. The generator, the drift and optionally the VARIANCE change at threshold levels, giving a piecewise homogeneous first- or second-order (Brownian) fluid queue. Setting every variance cell to zero reduces it to first order. The blocks it returns are exactly what mfq_ld_mean and mfq_ld_distr consume, so this closes that family rather than extending it.

METHOD. Per regime, the states with neither drift nor variance are censored out, and the remainder is split into three classes: positive drift with no variance, negative drift with no variance, and any state with variance. The density is anchored at both ends of the regime, forward from its lower threshold and backward from its upper one, and each direction's exponent comes from a QBD. The second-order terms are what force the QBD form: with a variance the balance equation is second order in the level, and multiplying through by a constant c chosen from the spectrum turns it into the matrix QUADRATIC of a discrete QBD, from whose R the exponent is recovered as K = (R - I) c. c is the smallest value that keeps the transformed triple substochastic, taken as the maximum over the drift states of -Q_ii/R_ii and over the variance states of the larger root of the discriminant, floored at 1.

WHICH QBD SOLVER, AND WHY NO NEW ONE WAS NEEDED. The reference calls QBD_CR(Bm, Lm, Fm), and the R it returns satisfies

Fm + R Lm + R^2 Bm = 0,

which is EXACTLY the contract of the port's existing qbd_R, F + R L + R^2 B = 0, under the direct mapping B = Bm, L = Lm, F = Fm. So the cyclic-reduction routine did not have to be transcribed.

That equation was established by MEASUREMENT, not by reading the argument names, and the first reading was wrong. QBD_CR's own error message rejects a triple whose sum is not "(sub)stochastic", which invites the conclusion that it wants the discrete-time form R = A0 + R A1 + R^2 A2 with A0 = Bm; the triples this caller builds sum to a GENERATOR, and the shift that reading implies (L = Lm - I) produces a completely different and wrong K. Feeding the actual regime-1 triple of a two-state first-order instance through MATLAB and evaluating all four candidate residuals settles it in one run: 2.2e-16 for the form above against 1.25, 1.5 and 2.75 for the others. The test asserts that residual directly on the port's own R, so the identity is pinned rather than inferred.

The boundary system then couples the K+1 point masses to the 2K density initial vectors through flux conservation at every threshold, the boundary conditions that reflective and absorbing states impose, and continuity of the second-order density across a threshold. One normalization row replaces the first flux equation.

REFERENCE DEFECT, reported and NOT worked around. mfq_ld_solve CRASHES ON ITS OWN DEFAULT ARGUMENTS whenever there is more than one regime:

Q = {[-2 2;1 -1],[-3 3;2 -2]}; R = {diag([1 -1]), diag([0.5 -2])};
S = {zeros(2), zeros(2)};  mfq_ld_solve(Q,R,S,[1 3])
Error: The logical indices contain a true value outside of the array
bounds.  SecondOrderLevelDependentFluidSolve line 202.

Line 24 defaults boundaryL to zeros(K,N), a K x N MATRIX, while the documented contract (line 8) and every use site want ONE ITEM PER BACKGROUND STATE, a length-N vector; line 202 then indexes ix = 1:N with that K x N logical mask. boundaryU inherits it. K = 1 survives only because zeros(1,N) happens to be the right shape. Passing explicit length-N vectors works and gives a correct answer. This port takes the boundary flags as length-N vectors, which is the documented contract, and defaults them to reflective, so it does not reproduce the crash.

ARITHMETIC. Templated on T and gated on num_traits<T>::has_transcendental: the QBD iteration is tolerance-terminated and expm is a Pade approximation. As in mfq_ld_distr, the ONLY eigenvalue computation is the branch test inside the normalization's integral of a matrix exponential, which selects between two algebraically equivalent formulas and never supplies a value that reaches the result, so Real instantiation is honest here. Contrast mfq_multiregime, where the Schur factors ARE the basis of the answer and double is the ceiling.

Definition in file mfq_ld_solve.h.

Namespaces

Enumerations

Functions

Detailed Description