qbd_rap.h File Reference

Quasi-birth-death processes with rational arrival process components, and the RAP/RAP/1 queue built on top of them. More...

#include <cmath>
#include <cstddef>
#include <string>
#include <vector>
#include "line/api/mam/map_moment.h"
#include "line/api/mam/map_transform.h"
#include "line/api/mam/mmap_lambda.h"
#include "line/api/mam/qbd_r.h"
#include "line/num/number.h"
#include "line/util/eig.h"
#include "line/util/error.h"
#include "line/util/linalg.h"
#include "line/util/matrix.h"

Include dependency graph for qbd_rap.h:

Go to the source code of this file.

Classes
struct	line::mam::QbdRapResult< T >
	Everything qbd_rap returns. More...
struct	line::mam::QbdRapRap1Result< T >
	Everything qbd_raprap1 returns. More...

Namespaces
namespace	line
namespace	line::mam

Functions
template<class T>
QbdRapResult< T >	line::mam::qbd_rap (const Matrix< T > &A0, const Matrix< T > &A1, const Matrix< T > &A2, const Matrix< T > &B0, const Matrix< T > &B1, std::size_t numLevels)
	Equilibrium analysis of a QBD with RAP components (qbd_rap.m).
template<class T>
QbdRapResult< T >	line::mam::qbd_rap (const Matrix< T > &A0, const Matrix< T > &A1, const Matrix< T > &A2)
	qbd_rap with the reference defaults, B0 = A0, B1 = A1 and 20 levels.
template<class T>
QbdRapRap1Result< T >	line::mam::qbd_raprap1 (const Map< T > &arrival, const Map< T > &service_in, const T &util)
	RAP/RAP/1 queue (qbd_raprap1.m).
template<class T>
QbdRapRap1Result< T >	line::mam::qbd_raprap1 (const Map< T > &arrival, const Map< T > &service)
	qbd_raprap1 without rescaling the service process.

Detailed Description

Quasi-birth-death processes with rational arrival process components, and the RAP/RAP/1 queue built on top of them.

Templated port of matlab/src/api/mam/qbd_rap.m (including its local qbd_rap_g) and matlab/src/api/mam/qbd_raprap1.m, following N. G. Bean and B. F. Nielsen, "Quasi-Birth-and-Death Processes with Rational Arrival Process Components", Stochastic Models 26(3), 2010, 309-334. The equilibrium construction is their Theorem 7 and the stability test their Corollary 8.

The process is given by its repeating blocks (A0, A1, A2) – A0 up, A1 local, A2 down – and its boundary blocks (B0, B1). Unlike a Markovian QBD the blocks need not be nonnegative; they are only required to be conservative, (A0 + A1 + A2) e = 0 and (B0 + B1) e = 0. The analysis rests on the prediction-process interpretation of a RAP, which is what lets a QBD argument survive the loss of nonnegativity.

Theorem 7, step by step:

1. G solves A0 G^2 + A1 G + A2 = 0.
2. U = A1 + A0 G.
3. R = A0 (-U)^-1.
4. pihat0 (B1 + R A2) = 0 with pihat0 e = 1.
5. pi_0 = K pihat0 with K chosen so pi_0 (I - R)^-1 e = 1.
6. pi_n = pi_0 R^n. Positive recurrence holds iff Sp(R) < 1 and step 4 has a unique solution.

COMPUTING G. The blocks are not nonnegative, so logarithmic and cyclic reduction carry no convergence guarantee, and the paper leaves the general case explicitly open (Section 6). Two paths, exactly as in the reference:

• A2 of rank one, A2 = u v: then G = e v / (v e) solves the equation in closed form. Conservativity gives (A0 + A1) e = -A2 e = -u (v e), and G is idempotent, so A0 G^2 + A1 G = (A0 + A1) e v/(v e) = -u v = -A2. This is the case of the paper's own example.
• otherwise, natural functional iteration G <- (-A1)^-1 (A2 + A0 G^2) as a warm start, then Newton on the Sylvester-form Jacobian (A0 G + A1) H + A0 H G = -(A0 G^2 + A1 G + A2), solved through its Kronecker expansion (I (x) (A0 G + A1) + G^T (x) A0) vec(H). An unconverged G is never returned: the residual and the constraint G e = e are both checked and a failure raises NumericError carrying both numbers.

DIVERGENCES FROM THE REFERENCE, all in how a quantity is EXTRACTED rather than in what it is, and all tested:

• the right factor v of a rank-one A2 is taken as the row of A2 with the largest infinity norm instead of the top right singular vector. G depends on v only through v/(v e), and every nonzero row of a rank-one matrix is a scalar multiple of v, so the two agree exactly; this keeps the step inside the templated arithmetic instead of routing it through a double-precision SVD. The rank test itself still uses the singular values (util/eig.h), matching the reference's sv(2) <= 1e-10 sv(1).
• the boundary vector of step 4 is obtained from the linear system x V = 0, sum(x) = 1 (qbd_detail::statvec) rather than from the last right singular vector of V^T. The solution is unique up to scale precisely when the reference's own second-smallest-singular-value test passes, so the accept/reject decision is unchanged, but the vector is computed at the working precision instead of in double.
• rcond(-U) is replaced by the EXACT reciprocal 1-norm condition number 1 / (||X||_1 ||X^-1||_1). MATLAB's rcond only estimates that quantity.

ARITHMETIC. Gated on num_traits<T>::has_transcendental: the general G is a fixed-point iteration plus Newton driven to a tolerance, and Sp(R) is the modulus of an eigenvalue, which is algebraic and not rational. Sp(R) is computed by converting R to double and calling LAPACK (util/eig.h), so at Real<D> the STABILITY GATE is only double-accurate; every returned quantity is computed at the full working precision. The gate is a comparison against 1 - 1e-12 and any model that close to the null-recurrent boundary has an unbounded queue anyway, which is why the precision loss is confined there.

Definition in file qbd_rap.h.