The MAP/D/1 FCFS queue: deterministic service of length s fed by a Markovian arrival process. More...

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

Classes
struct	line::qsys::MapD1Result< T >
	Return value of qsys_mapd1, mirroring the MATLAB result struct. More...

Namespaces
namespace	line
namespace	line::qsys

Functions
template<class T>
MapD1Result< T >	line::qsys::qsys_mapd1 (const mam::Map< T > &arrival, const T &s, std::size_t dist_size, unsigned max_arrivals, std::size_t max_levels, const T &tol)
	MAP/D/1 by the exact embedded M/G/1-type chain.
template<class T>
MapD1Result< T >	line::qsys::qsys_mapd1 (const mam::Map< T > &arrival, const T &s)
	qsys_mapd1 with 100 materialized levels, an arrival-count cap of 4096, a level cap of 20000 and tolerance 1e-14.

Detailed Description

The MAP/D/1 FCFS queue: deterministic service of length s fed by a Markovian arrival process.

ALGORITHM, AND HOW IT DIFFERS FROM THE MATLAB REFERENCE. matlab/src/api/qsys/qsys_mapd1.m delegates to qsys_mapdc, which calls Q-MAM's Q_CT_MAP_D_C, which is not transcribed here. This port computes the SAME quantities exactly, by the embedded-chain route, and NOT by an Erlang-k approximation: a deterministic service time has no phase-type representation, so no QBD can carry it, but an M/G/1-type Markov chain can.

The construction is in four exact steps.

MAP counting probabilities over one service. P_k(u) is the matrix whose (i,j) entry is P(k arrivals in [0,u], phase j at u | phase i at 0). They satisfy dP_k/du = P_k D0 + P_{k-1} D1, so the row (P_0(s) ... P_K(s)) is the top block row of exp(M s) with M the (K+1)-block bidiagonal matrix carrying D0 on the diagonal and D1 on the superdiagonal. One expm gives all of them. K is grown until the missing mass max_i (1 - sum_k P_k(s)e) is below the tolerance, so the truncation is a controlled quantity and not a modelling choice.
Time-in-level integrals. I_k = integral_0^s P_k(u) du follows from the same differential equation without a second expm: integrating gives P_k(s) - P_k(0) = I_k D0 + I_{k-1} D1, hence I_0 = (P_0(s) - I) D0^-1 and I_k = (P_k(s) - I_{k-1} D1) D0^-1. D0 is a non-singular sub-generator.
The embedded chain at departure epochs. With the level equal to the number left behind by a departure and the phase equal to the arrival phase, the chain is skip-free to the left with blocks A_k = P_k(s) from a busy level and B_k = (-D0)^-1 D1 A_k from the empty level (wait for the next arrival, then serve). Its stationary vector comes from Ramaswami's recursion: with G the minimal solution of G = sum_k A_k G^k, and the partial sums Ahat_i = sum_{k>=i} A_k G^(k-i), Bhat_i = sum_{k>=i} B_k G^(k-i),
```
x_0 (B_0 + Bhat_1 (I - Ahat_1)^-1 A_0) = x_0,
x_n = [ x_0 Bhat_n + sum_{k=1}^{n-1} x_k Ahat_{n-k+1} ] (I - Ahat_1)^-1.
```

The time-stationary distribution. The departure-epoch vector x is NOT the time-stationary one unless the arrivals are Poisson, so the two are not interchanged here. Averaging over an inter-departure cycle with the integrals of step 2,
```
p_0   = lambda x_0 (-D0)^-1 e,
p_n   = lambda [ y_0 I_{n-1} e + sum_{m=1}^{n} x_m I_{n-m} e ],  n >= 1,
```
with y_0 = x_0 (-D0)^-1 D1 the phase at the arrival that ends an idle period. Since I_k e sums to s e over k, this construction satisfies sum_n p_n = lambda (E[idle] + s) = 1 and p_0 = 1 - lambda s identically, which the test file asserts. The mean follows in closed form from sum_k I_k and sum_k k I_k without materializing the levels.

meanWaitingTime is (L - rho)/lambda by Little's law and meanSojournTime is that plus s.

MEASURED AGREEMENT (MATLAB R2025a, T = double). meanQueueLength agrees with the reference; meanWaitingTime does not, and the reference is the one that is wrong – see the defect note below.

M/D/1 collapse, qsys_mapd1(D0 = [-2], D1 = [2], s = 1/3), rho = 2/3. Textbook L = rho + rho^2/(2(1-rho)) = 4/3 and Wq = rho s/(2(1-rho)) = 1/3. The port returns 1.33333333333307 and 0.333333333333203, i.e. the textbook values to 2.0e-13 and 4.0e-13. MATLAB's meanQueueLength is 1.333333330345585, itself 2.2e-9 below the textbook value (its own maxNumComp truncation), so the port and MATLAB differ by 2.24e-9.
Correlated MMPP2 arrivals D0 = [-2.5 0.2; 0.1 -0.7], D1 = diag(2.3, 0.6) (lambda = 7/6), s = 0.4, rho = 7/15: MATLAB meanQueueLength 1.118300079597477, port 1.11830008358428, relative difference 3.57e-9.
Erlang-2 arrivals D0 = [-4 4; 0 -4], D1 = [0 0; 4 0] (lambda = 2), s = 0.3, rho = 0.6: MATLAB meanQueueLength 0.7758216464528508, port 0.775821646757718, relative difference 3.93e-10.

Both non-Poisson values were confirmed independently of Q-MAM by the Erlang-k limit through LINE's own MATLAB qbd_mapmap1: replacing the deterministic service by an Erlang-k of the same mean and Richardson- extrapolating the O(1/k) convergence from k = 40 and k = 80 gives 1.11832240806 and 0.775808057317, which match the port to 2.0e-5 and 1.7e-5, the residual of the extrapolation itself. The port additionally satisfies sum_n p_n = 1 and p_0 = 1 - rho to 1e-15 on all three instances, both being identities of the construction rather than fitted quantities.

At T = Real50 the port reproduces its own double results to 2e-14, so the double values above are not precision-limited.

DEFECT IN THE MATLAB REFERENCE (reported, not fixed here). qsys_mapdc, and therefore qsys_mapd1, computes meanWaitingTime as a left-rectangle sum of the survival function of the Q-MAM waiting-time CDF with step s/numSteps and numSteps defaulting to 1. Two consequences: (a) even for Poisson arrivals the default is a one-point quadrature. On the M/D/1 instance above it returns 0.4444444443040794 against the exact 1/3, a +33% error, converging as O(1/numSteps): 0.3350694 at numSteps = 64 and 0.3334418 at numSteps = 1024. (b) for non-Poisson arrivals it does not converge to the right value at all. On the correlated instance it converges to 0.3540778458263887 while Little's law applied to its own (correct) meanQueueLength gives 0.5585429253692660; on the Erlang-2 instance it converges to 0.1428343846295046 against 0.0879108232264254. Errors of -37% and +62%, in opposite directions, so this is not a quadrature artefact. The reference's meanQueueLength is correct in all three cases – that is what the Erlang-k cross-check above establishes – so this port takes Little's law as the definition of meanWaitingTime and does not reproduce the reference's waiting-time numbers. Reproduction: run qsys_mapd1([-2],[2],1/3,'numSteps',N) for N = 1, 64, 1024 and compare meanWaitingTime against 1/3.

ARITHMETIC. Gated on num_traits<T>::has_transcendental for two independent reasons: step 1 calls expm, a scaling-and-squaring Pade approximation that is tolerance-controlled and cannot be exact in any arithmetic, and step 3 computes G by a fixed-point iteration that does not terminate in a finite number of field operations. Steps 2 and 4 are finite exact matrix algebra given P_k(s) and G, and add no error of their own.

Definition in file qsys_mapd1.h.

Classes

Namespaces

Functions

Detailed Description