Multi-regime FEEDBACK Markovian fluid queue: density, density derivative and distribution of the fluid level. More...

#include <cmath>
#include <cstddef>
#include <vector>
#include "line/api/mam/mfq_solve.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_multiregime.h:

Classes
struct	line::mam::MultiRegimeResult
	Return value of mfq_multiregime: one row of N per-state values per point. More...

Namespaces
namespace	line
namespace	line::mam

Functions
MultiRegimeResult	line::mam::mfq_multiregime (const std::vector< Matrix< double > > &Q, const std::vector< std::vector< double > > &R, const std::vector< Matrix< double > > &Qt, const std::vector< std::vector< double > > &Rt, const std::vector< double > &Thr, const std::vector< double > &pdfpoints, const std::vector< double > &cdfpoints)
	Multi-regime feedback fluid queue.

Detailed Description

Multi-regime FEEDBACK Markovian fluid queue: density, density derivative and distribution of the fluid level.

Port of matlab/src/api/mam/mfq_multiregime.m and the BUTools multiregime it wraps, which implements H. E. Kankaya and N. Akar, "Solving Multi-Regime Feedback Fluid Queues". The generator and the drift rates are regime dependent, and separate FEEDBACK generators and rates govern the behaviour at each threshold, which is what distinguishes this from the level-dependent model behind mfq_ld_solve: there the thresholds only separate regimes, here the process can behave differently while sitting exactly on one.

METHOD. Per regime, the zero-drift states are censored out and the remaining generator is rescaled by the drifts, giving A = Qbar diag(1/R). The spectrum of A splits into zero, negative and positive parts, and the three corresponding invariant subspaces carry the constant, the decaying-upward and the decaying-downward components of the density. A block-triangularizing similarity Y is built from the ordered real Schur form of A plus two Sylvester solves, after which the density in regime k is

pi_k(x) = a0_k L0_k + an_k exp(An_k (x - T_k)) Ln_k

ap_k exp(-Ap_k (T_{k+1} - x)) Lp_k.

The unknown coefficients, together with the point masses at the K+1 thresholds, solve one linear system assembled from the reference's equations (8)-(16): flow balance at each boundary, the boundary conditions that the feedback rates impose on each state, and one normalization row that replaces the first balance equation.

WHY THE SCHUR FORM, AND WHY NOT EIGENVECTORS. The eigenvector basis exists only when A is diagonalizable and is complex whenever the spectrum is; substituting it here would be a different algorithm wearing the same name. The real Schur basis is orthogonal and real for every real A, and the ordering by eigenvalue sign class is what isolates the three subspaces. That is why line/util/eig.h grew schur_decomposition and schur_reorder for this function.

DOUBLE ONLY, AND HONESTLY SO. Unlike mfq_ld_distr – where the only eigenvalues merely SELECT a branch and never reach a returned number – here the Schur factors Z and T are the basis in which An, Ap, L0, Ln and Lp are expressed, so they enter every value the function returns. util/eig.h is double-only by design (the eigenvalues of a rational matrix are algebraic, not rational, and LAPACK cannot supply a multiprecision QR), so this function takes and returns Matrix<double> rather than being templated. Instantiating it at Real50 would advertise a precision it cannot deliver, since everything downstream of the Schur step would be carrying double-accurate inputs. It is registered as {Double} for that reason and no other. When the build has no LAPACK it refuses by name through schur_decomposition rather than guessing a basis.

COMPARING THE OUTPUT AGAINST MATLAB. The Cdfm column at level zero is an EXACT ZERO in this port and comes back from MATLAB as -6.26e-17, its solve having rounded. Any comparison of that column must therefore be ABSOLUTE, not relative: a relative test against a value that is exactly zero fails on a difference of one ulp and reports a defect that is not there. The same holds for any state carrying no mass at a threshold. This is stated here rather than only in the test so that a later reader does not "tidy" the absolute comparison into a relative one and manufacture a failure.

THE SYLVESTER EQUATIONS ARE NOT THE OBSTACLE. A X + X B = C vectorizes to (I kron A + B^T kron I) vec(X) = vec(C), one ordinary linear solve at these block sizes, needing no Bartels-Stewart and no Schur form of its own. The reference calls MATLAB's sylvester, which is Bartels-Stewart; the two compute the same X, and at the sizes reached here the direct solve is not the slower one.

Definition in file mfq_multiregime.h.

Classes

Namespaces

Functions

Detailed Description