Marie's iterative aggregation-decomposition for closed networks with FCFS general (Coxian) service. More...

#include <cmath>
#include <cstddef>
#include <limits>
#include <vector>
#include "line/api/mc/ctmc_solve.h"
#include "line/api/pfqn/pfqn_mva.h"
#include "line/api/pfqn/pfqn_mvams.h"
#include "line/num/number.h"
#include "line/util/error.h"
#include "line/util/matrix.h"

Include dependency graph for pfqn_marie.h:

Classes
struct	line::pfqn::MarieCoxFit< T >
	Coxian phase representation: phase rates and per-phase completion probabilities. More...
struct	line::pfqn::MarieCdScaling< T >
	Class-dependent scaling of one station, tabulated on the integer population box. More...
struct	line::pfqn::MarieResult< T >
	Result of pfqn_marie, mirroring the six MATLAB outputs. More...

Namespaces
namespace	line
namespace	line::pfqn

Functions
template<class T>
MarieCoxFit< T >	line::pfqn::marie_cox_fit (const T &mean, const T &scv)
	Closed-form Coxian fit of a mean and an SCV (matlab/src/lang/processes/Coxian.m, fitMeanAndSCV), with the branch thresholds at CoarseTol = 1e-3.
template<class T>
std::vector< T >	line::pfqn::marie_cd_eval (const MarieCdScaling< T > &cd, const std::vector< T > &nv)
	Evaluate a class-dependent scaling at a real-valued population vector (cdscale_eval in the reference): the interpolated ratio, guarded against a non-positive or non-finite value and clamped to [1e-3, 1e3].
template<class T>
MarieResult< T >	line::pfqn::pfqn_marie (const Matrix< T > &L, const std::vector< int > &N, const std::vector< T > &Z, const Matrix< T > &scv, double tol, int maxiter, const std::vector< int > &nservers)
	Marie's method for a closed network with FCFS Coxian service.
template<class T>
MarieResult< T >	line::pfqn::pfqn_marie (const Matrix< T > &L, const std::vector< int > &N, const std::vector< T > &Z, const Matrix< T > &scv)
	Reference defaults: tol 1e-8, maxiter 1000, single server everywhere.

Detailed Description

Marie's iterative aggregation-decomposition for closed networks with FCFS general (Coxian) service.

Templated port of matlab/src/api/pfqn/pfqn_marie.m. There is no JAR counterpart, so MATLAB is the only reference.

SINGLE CLASS (R = 1). Each station's service is given a Coxian phase representation matching its mean L(i) and its squared coefficient of variation scv(i). The aggregate model is the exact load-dependent product form pfqn_mvald driven by a multiplier lattice mu(i,n), from which the marginal queue-length distribution P_i(n) at the full population is read off. Flow balance across the n <-> n+1 cut of that birth-death marginal gives the complementary arrival rate seen by station i,

lambda_i(n) = (mu(i,n+1)/L(i)) P_i(n+1) / P_i(n), n = 0, ..., N-1,

and the lambda(n)/Cox/1(-m) isolation chain – states (n,k) with n present and the head job in phase k, phase rates scaled by min(n,m) for m servers – is solved for its stationary distribution. Its conditional departure rate

mu_i(n) = sum_k p(n,k) rate_k phi_k min(n,m) / sum_k p(n,k)

is converted back to a multiplier (multiplied by L(i)) and fed to the next aggregate solve. The iteration stops when max|mu_new - mu| < tol. For exponential service (scv == 1) the isolation chain is the M/M/1(-m) queue, mu_i(n) = min(n,m)/L(i), so the multiplier lattice is the initial one, the first iteration already reproduces exact product form and the loop exits on the second pass.

MULTIPLE CLASSES (R > 1). Exponential and class-independent demands at every station is genuine BCMP FCFS and is dispatched to exact pfqn_mva. Otherwise the aggregate is a Schweitzer-style multiclass AMVA in which the class-r demand at station i is divided by a class-dependent scaling beta_{i,r}(nvec) = muCox_{i,r}(nvec) / muExp_{i,r}(nvec), the ratio of the conditional class-r throughput of a multiclass Cox/1 FCFS isolation chain to that of the same chain with exponential service of the same means. beta == 1 therefore recovers plain FCFS AMVA and carries only the non-exponential correction. The isolation chain is fed the aggregate per-class throughput (Baynat-Dallery isolation) and the outer loop iterates to a fixed point on X. beta is tabulated on the integer population box and read at the real-valued arrival-instant populations by multilinear interpolation; the reference guards a non-positive or non-finite ratio by falling back to 1 and clamps the result to [1e-3, 1e3], and both are reproduced.

COXIAN FIT. Coxian.fitMeanAndSCV is closed form and is reproduced here as marie_cox_fit, with tol = GlobalConstants.CoarseTol = 1e-3 (matlab/lineStart.m):

|scv - 1| <= tol exponential, n = 1, mu = [1/mean], phi = [1] 0.5 + tol < scv < 1 - tol hypoexponential, n = 2, mu = 2/mean/(1 +- sqrt(2 scv - 1)), phi = [0,1] scv <= 0.5 + tol Erlang, n = ceil(1/scv), mu = (n/mean) 1, phi = e_n scv > 1 + tol Coxian-2, mu = [2/mean, 1/(scv mean)], phi = [1 - 1/(2 scv), 1]

The exponential, hypoexponential and Coxian-2 branches match the requested mean and scv exactly. The Erlang branch does NOT: ceil(1/scv) is an integer, so the fitted scv is 1/ceil(1/scv) <= scv, with equality only when 1/scv is an integer. That is the reference's behaviour and is reproduced, not improved. The phase count also goes through a double-valued ceiling, which is the one place the fit is not a pure field computation.

DIVERGENCE from the reference, deliberate and documented. MATLAB solves the isolation chains as the overdetermined least-squares system [Q'; ones] p = [0; 1] via backslash. This port calls mc::ctmc_solve, which replaces one balance equation by the normalization and solves the resulting square system, and which additionally splits a reducible generator into its weakly connected components. The two are the same computation whenever the isolation chain is irreducible, which is the case whenever every complementary arrival rate lambda_i(n) is positive – the only regime in which the reference's own least-squares answer is a probability vector at all.

Arithmetic: INEXACT BY CONSTRUCTION, and additionally gated on has_transcendental. The method is a fixed-point decomposition stopped on a tolerance, the Coxian fit of an scv in (0.5, 1) needs a square root that has no exact rational counterpart, and the multiclass path reports throughputs from an approximate MVA. It is therefore unavailable at T = Rational.

Definition in file pfqn_marie.h.

Classes

Namespaces

Functions

Detailed Description