Sojourn-time distribution of the M/G/1 processor-sharing queue. More...

#include <algorithm>
#include <cmath>
#include <complex>
#include <cstddef>
#include <functional>
#include <limits>
#include <vector>
#include "line/api/pfqn/pfqn_asympt_common.h"
#include "line/num/number.h"
#include "line/util/error.h"
#include "line/util/expm.h"
#include "line/util/lu.h"
#include "line/util/matrix.h"

Include dependency graph for qsys_mg1_ps.h:

Classes
struct	line::qsys::Mg1PsOptions
	Options of qsys_mg1_ps, mirroring the reference's name-value pairs. More...
struct	line::qsys::Mg1PsResult
	Everything qsys_mg1_ps returns. More...

Namespaces
namespace	line
namespace	line::qsys

Typedefs
using	line::qsys::Cplx = std::complex<double>

Functions
Mg1PsResult	line::qsys::qsys_mg1_ps (double lambda, const std::vector< double > &alpha, const Matrix< double > &Tmat, const Mg1PsOptions &opt)
	Sojourn-time distribution of the M/G/1 processor-sharing queue.

Detailed Description

Sojourn-time distribution of the M/G/1 processor-sharing queue.

Port of matlab/src/api/qsys/qsys_mg1_ps.m (twins Qsys_mg1_ps.java, python/line_solver/api/qsys/mg1ps.py). MATLAB is the reference.

Jobs arrive Poisson at rate lambda at one egalitarian processor-sharing server whose service requirement has LST bhat(tau) and mean m1. Writing V(x) for the sojourn of a tagged job of requirement x and rho = lambda m1 < 1, Ott (1984) and Yashkov (1983) give

E[exp(-s V(x))] = (1 - rho) / D(s,x),

D(s,x) being the inverse Laplace transform, evaluated at x, of

f(tau;s) = [ (1-rho) tau^2 - (1-rho) lambda (1-bhat(tau)) tau
             + s rho tau - s lambda (1-bhat(tau)) ]
           / [ tau^2 (tau - s - lambda (1-bhat(tau))) ].

The transform is exact but IMPLICIT: f has to be inverted in tau. For phase-type service f is a proper rational function of tau, the double pole at the origin cancels, and D(s,x) comes out in closed form as a finite sum of residues – so the M/PH/1-PS queue, hence every service law that can be fitted phase-type, is exactly solvable. For a transform supplied as a callback, f is inverted numerically on a Bromwich contour placed to the right of the dominant singularity tau*(s), the unique root of tau = s + lambda (1 - bhat(tau)) in the right half plane, which the fixed point of that equation reaches at geometric rate rho.

THE CONDITIONAL SOJOURN IS ATOMIC ON THE LATTICE t = (k+1) x, and that is not a numerical artifact. Processor sharing gives every job in the system the same service rate, so if the k jobs present on arrival all outlive the tagged job and nothing arrives, the sojourn is exactly (k+1) x. The k = 0 atom, (1-rho) exp(-lambda x), is the probability of finding the system empty and sharing it with nobody, and it is the only one that stays exact for general service. It is removed before inverting in s; the remaining atoms make the conditional CDF JUMP and leave no density, so the density is reported as NaN on the lattice rather than as the finite garbage a smooth inversion returns there.

TWO SUBSTITUTIONS FOR MATLAB, both self-contained rather than approximate:

roots(P) is computed by Durand-Kerner rather than by the eigenvalues of the companion matrix, so this header does not require LAPACK. The polynomial has degree n+1 in the phase count, which is small, and the iteration converges to machine precision on it; the roots are a SET, and the residue sum that consumes them does not depend on their order.
The Golub-Welsch call behind the panelled quadrature is the tree's own gauss_legendre, which computes the same Legendre nodes by Newton. The reference notes that the rule is fixed so that every codebase shares it; the nodes agree to machine precision, so they do.

ARITHMETIC: double. The inversion, the root finding and the contour are all inherently floating point, and the reference's own tolerances are absolute in double.

Definition in file qsys_mg1_ps.h.

Classes

Namespaces

Typedefs

Functions

Detailed Description