Sojourn time distribution in a MAP/M/1 processor-sharing queue. More...

#include <cmath>
#include <cstddef>
#include <vector>
#include "line/api/mam/map_moment.h"
#include "line/api/mam/qbd_r.h"
#include "line/api/mc/ctmc_solve.h"
#include "line/num/number.h"
#include "line/util/error.h"
#include "line/util/linalg.h"
#include "line/util/matrix.h"

Include dependency graph for map_m1ps.h:

Classes
struct	line::mam::MapM1psResult< T >
	What the two MAP/M/1-PS sojourn entry points return. More...

Namespaces
namespace	line
namespace	line::mam

Functions
template<class T>
T	line::mam::map_compute_R_residual (const Matrix< T > &C, const Matrix< T > &D, const T &mu, const Matrix< T > &R)
	Residual \|\|D + R (C - mu I) + mu R^2\|\|_inf of the MAP/M/1 rate equation.
template<class T>
Matrix< T >	line::mam::map_compute_R (const Matrix< T > &C, const Matrix< T > &D, const T &mu, unsigned max_iter, const T &tol)
	Rate matrix R of a MAP/M/1 queue, the minimal nonnegative solution of D + R (C - mu I) + mu R^2 = 0, by the iteration R <- -D (C - mu I + mu R)^-1 (map_compute_R.m).
template<class T>
Matrix< T >	line::mam::map_compute_R (const Matrix< T > &C, const Matrix< T > &D, const T &mu)
	map_compute_R with the reference defaults, 1000 iterations and tolerance 1e-10.
template<class T>
Matrix< T >	line::mam::map_compute_R_quadratic (const Matrix< T > &C, const Matrix< T > &D, const T &mu, unsigned max_iter, const T &tol)
	The same R by the other splitting, R <- (D + mu R^2) (mu I - C)^-1, warm started at -D (C - mu I)^-1 and with the scalar case solved in closed form (the private compute_R_matrix of map_m1ps_cdfrespt.m).
template<class T>
Matrix< T >	line::mam::map_compute_R_quadratic (const Matrix< T > &C, const Matrix< T > &D, const T &mu)
	map_compute_R_quadratic with the reference defaults, 5000 iterations, 1e-10.
template<class T>
std::vector< std::vector< std::vector< T > > >	line::mam::map_m1ps_h_recursive (const Matrix< T > &C, const Matrix< T > &D, const T &mu, std::size_t N, std::size_t K)
	The vectors h_{n,k} of Theorem 1 (map_m1ps_h_recursive.m).
template<class T>
MapM1psResult< T >	line::mam::map_m1ps_sojourn (const Matrix< T > &C, const Matrix< T > &D, const T &mu, const std::vector< T > &x, const T &epsilon, const T &epsilon_prime)
	Complementary sojourn time distribution of a MAP/M/1-PS queue (map_m1ps_sojourn.m).
template<class T>
MapM1psResult< T >	line::mam::map_m1ps_sojourn (const Matrix< T > &C, const Matrix< T > &D, const T &mu, const std::vector< T > &x)
	map_m1ps_sojourn with the reference defaults, epsilon 1e-11 and 1e-10.
template<class T>
MapM1psResult< T >	line::mam::map_m1ps_cdfrespt (const Matrix< T > &C, const Matrix< T > &D, const T &mu, const std::vector< T > &x, const T &epsilon, const T &epsilon_prime)
	Complementary sojourn time distribution of a MAP/M/1-PS queue by the spectral-radius truncation (map_m1ps_cdfrespt.m).
template<class T>
MapM1psResult< T >	line::mam::map_m1ps_cdfrespt (const Matrix< T > &C, const Matrix< T > &D, const T &mu, const std::vector< T > &x)
	map_m1ps_cdfrespt with the reference defaults, epsilon 1e-11 and 1e-10.

Detailed Description

Sojourn time distribution in a MAP/M/1 processor-sharing queue.

Templated port of matlab/src/api/mam/map_compute_R.m, map_m1ps_h_recursive.m, map_m1ps_sojourn.m and map_m1ps_cdfrespt.m, which implement Theorem 1 of H. Masuyama and T. Takine, "Sojourn time distribution in a MAP/M/1 processor-sharing queue", Operations Research Letters 31(6), 2003, 406-412.

The arrival process is the MAP (C, D) – C carries the hidden transitions, D the arrivals – and service is exponential of rate mu shared equally, so with n jobs present each is served at rate mu/n. The queue length is a QBD whose rate matrix R is the minimal nonnegative solution of

D + R (C - mu I) + mu R^2 = 0,

and the complementary sojourn time distribution is, by uniformization at theta + mu with theta = max_i |C_ii|,

W^c(x) = (1/lambda) sum_n pi_0 R^n D sum_k e^-(theta+mu)x
         ((theta+mu) x)^k / k! h_{n,k},

with pi_0 = pi (I - R) and the vectors h_{n,k} from the recursion

h_{n,0}   = e
h_{n,k+1} = [ n mu/(n+1) h_{n-1,k} + (theta I + C) h_{n,k}
              + D h_{n+1,k} ] / (theta + mu),   h_{-1,k} = 0.

ARITHMETIC.

map_m1ps_h_recursive is a FINITE recursion in the entries of C and D, with theta a maximum of absolute diagonal entries, so it instantiates at every arithmetic including Rational and returns exact fractions.
map_compute_R is a fixed-point iteration driven to a tolerance and is gated on num_traits<T>::has_transcendental, for the same reason as qbd_R (see qbd_r.h).
map_m1ps_sojourn and map_m1ps_cdfrespt need the Poisson weights e^-a a^k / k! and are gated as well.

THE TWO ENTRY POINTS ARE NOT THE SAME FUNCTION, despite identical signatures and identical documentation in the reference. They differ in three ways, all reproduced here:

The queue-length truncation. map_m1ps_sojourn finds the smallest N with (1/lambda) sum_{n<=N} pi_0 R^n D e > 1 - epsilon, scanning n = 0..1000, and falls back to N = 100 when the scan never gets there. map_m1ps_cdfrespt instead estimates N from the spectral radius of R, N = ceil(log(epsilon (1 - sp)) / log(sp)), clamps it to [10, 10000] and then truncates AGAIN at run time as soon as ||pi_0 R^n D||_inf drops below epsilon/100.
The stationary vector. map_m1ps_sojourn solves the (M+1) x M overdetermined system [Q; e'] x = [0; 1] in the least-squares sense, whereas map_m1ps_cdfrespt replaces the LAST ROW of Q' by e' and solves the resulting square system. Both give the stationary vector of an irreducible generator; the port uses the square solve of ctmc_solve for both, which is the same vector.
R itself. map_m1ps_sojourn calls map_compute_R, iterating R <- -D (C - mu I + mu R)^-1. map_m1ps_cdfrespt has a PRIVATE compute_R_matrix that iterates the different splitting R <- (D + mu R^2) (mu I - C)^-1 from the warm start -D (C - mu I)^-1, with 5000 rather than 1000 iterations, and for M = 1 solves the scalar quadratic in closed form. The two fixed points coincide – both are the minimal nonnegative solution – so the two are exposed here as map_compute_R and map_compute_R_quadratic and the tests check they agree to the residual of the defining equation.

REFERENCE DEFECT. The second output of both functions, W_bar_n, is documented as "the conditional complementary distribution for customers finding n customers in the system", but both compute it as sum(sum_k)/M – the arithmetic MEAN OVER PHASES of the uniformized h-weighted sum, with no reference to the phase distribution at an arrival and no normalization by the probability of finding n customers. It is not a conditional distribution: on the M/M/1-PS instance of the tests (M = 1, lambda = 0.8, mu = 1) the n = 0 curve at x = 0 is 1.0 and DECREASES correctly, but for M > 1 the phases are weighted uniformly rather than by pie, so it is not a probability of anything. The port returns it under the name w_bar_n_unweighted to make the meaning explicit, and computes it identically so a caller comparing against MATLAB sees the same numbers.

Definition in file map_m1ps.h.

Classes

Namespaces

Functions

Detailed Description