getMomentStationTable.m

function [MomentStationTable, mom] = getMomentStationTable(self, order)
% GETMOMENTSTATIONTABLE Exact higher moments of the total queue length.
%
% [MOMENTSTATIONTABLE, MOM] = GETMOMENTSTATIONTABLE(SELF) returns a table with
% one row per Station giving the moments of the TOTAL queue length at that
% station, Q_i = sum_r n(i,r):
%   QLen, QLenVar, QLenSCV
%
% [..] = GETMOMENTSTATIONTABLE(SELF, ORDER) selects which moment orders to
% report. ORDER is a set: a scalar k is read as 1:k, "everything up to order
% k"; an explicit vector selects exactly those orders.
%   1        the mean only:             QLen
%   2        (default) mean and second moment: QLen, QLenVar, QLenSCV
%   3        also adds QLenM3 and QLenSkew
%   [1 3]    the mean and the third moment, without the variance
% Order 1 contributes QLen, order 2 contributes QLenVar and QLenSCV, order 3
% contributes QLenM3 and QLenSkew.
%
% Why this is a separate table from getMomentTable. The moments beyond the
% second are generated by differentiating with respect to x_i, the reciprocal
% of the capacity of station i, which scales the service times of ALL classes
% at that station at once. That parameter therefore produces moments of the
% station total, not of any one class: there is no per-class third moment to
% report, and inventing one by splitting the total would be fiction. The
% per-class second moments, which do exist, are in getMomentTable. The two are
% consistent: Var[Q_i] here equals the sum of getMomentTable's per-class
% covariances at station i over all class pairs.
%
% The ALGORITHM is chosen by the solver's method, set at construction, not by
% an argument here: a method is a property of the solver object, so passing one
% per call would let a single solver answer with two different algorithms.
%   SolverMVA(model)                   -> pfqn_sens_mom, exact, but it walks the
%                                         whole population lattice at a cost of
%                                         prod(N+1), so it is unusable once the
%                                         populations are large.
%   SolverMVA(model,'method','lin')    -> pfqn_sens_linearizer, approximate and
%                                         polynomial-time. The reference reports
%                                         relative errors below 2.1% on E[Q],
%                                         4.1% on E[Q^2] and 6.2% on E[Q^3].
% Any Linearizer-family method ('lin', 'amva.lin', 'egflin', 'gflin') takes the
% approximate path; every other method takes the exact one.
%
% Restricted to closed models. Mixed and open second moments are available per
% class from getMomentTable; the higher moments of the reference are stated for
% closed networks only.
%
% MOM is the underlying pfqn_sens_mom / pfqn_sens_linearizer struct, which also
% carries the cross-station covariance matrix .Cov that this table does not
% show.
%
% Reference: J. C. Strelen, "Moment Analysis for Closed Queuing Networks and
% its Linearizer", Performance Evaluation 11:127-142, 1990, Theorem 3.1 and
% equation (3.2).
%
% See also: getMomentTable, getAvgTable, getSensitivityTable.
 
if nargin < 2 || isempty(order)
    order = 2;
end
order = validateMomentOrder(order, 3);
% The algorithm is a property of the solver, not of this call: it comes from
% the method set at construction, e.g. SolverMVA(model,'method','lin'). Taking
% a per-call method argument here would have given the same solver object two
% disagreeing methods.
method = self.getOptions.method;
 
sn = self.model.getStruct();
R = sn.nclasses;
N = sn.njobs;
 
[~, D, Np, Z, ~, Ssrv, ~] = sn_get_product_form_params(sn);
queueIndices = find(sn.nodetype == NodeType.Queue);
Mq = numel(queueIndices);
Ztot = sum(Z, 1);
 
% The moments rest on the product-form identity Cov = L dQ/dL, which is a
% theorem about the product form: it follows from a d/da log G = E[n], which
% needs the state distribution to be exponential-family in the demands. Outside
% product form, L dQ/dL remains computable and is simply NOT a covariance, so
% differentiating a non-product-form solver would return a confident wrong
% number. Refuse rather than do that; exact moments of such a model come from
% SolverCTMC's stationary distribution, or from LDES via setReward.
if ~sn_has_product_form(sn)
    line_error(mfilename, 'getMomentStationTable requires a product-form model: the moment identity Cov[n,n] = L dQ/dL holds only under product form, so no correct value exists here. Use SolverCTMC (exact distribution) or SolverLDES with setReward for the moments of a non-product-form model.');
end
if any(isinf(N))
    line_error(mfilename, 'getMomentStationTable supports closed models only. Per-class second moments of an open or mixed model are available from getMomentTable.');
end
if any(Ssrv > 1)
    line_error(mfilename, 'getMomentStationTable supports single-server stations only: the higher-moment recursion of the reference is stated for load-independent stations. Per-class second moments of a multiserver model are available from getMomentTable.');
end
 
% Two algorithms exist for these moments, so the solver's method selects
% between them: any Linearizer-family method ('lin', 'amva.lin', 'egflin',
% 'gflin') takes the approximate polynomial-time path, and every other method
% takes the exact one. The solver has already rejected methods it does not
% support, so no further validation is needed here.
if isLinearizerMethod(method)
    mom = pfqn_sens_linearizer(D, Np, Ztot);
elseif isExactMvaMethod(method)
    mom = pfqn_sens_mom(D, Np, Ztot);
else
    % An approximate method with no hand-differentiated counterpart. The
    % identity does not care HOW the mean queue lengths were obtained, so the
    % solver itself is used as a mean-value oracle and differentiated
    % numerically. This is exactly the move the Linearizer makes analytically
    % (Strelen Sec. 5: differentiate the approximate fixed point, then apply the
    % exact identity), generalized to any product-form method. The moments
    % inherit the accuracy of that method's means.
    mom = momentsByFiniteDifference(self, sn, queueIndices);
end
 
Station = {};
QLen = []; QLenVar = []; QLenSCV = []; QLenM3 = []; QLenSkew = [];
for ist = 1:Mq
    if all(D(ist, :) <= 0)
        continue;   % no class visits this station
    end
    Station{end+1, 1} = sn.nodenames{queueIndices(ist)}; %#ok<AGROW>
    QLen(end+1, 1) = mom.m(ist);            %#ok<AGROW>
    QLenVar(end+1, 1) = mom.Var(ist);       %#ok<AGROW>
    if mom.m(ist) > 0
        QLenSCV(end+1, 1) = mom.Var(ist) / mom.m(ist)^2; %#ok<AGROW>
    else
        QLenSCV(end+1, 1) = NaN;            %#ok<AGROW>
    end
    QLenM3(end+1, 1) = mom.M3(ist);         %#ok<AGROW>
    QLenSkew(end+1, 1) = mom.Skew(ist);     %#ok<AGROW>
end
 
vars = {Station};
names = {'Station'};
if any(order == 1)
    vars{end+1} = QLen;     names{end+1} = 'QLen';
end
if any(order == 2)
    vars{end+1} = QLenVar;  names{end+1} = 'QLenVar';
    vars{end+1} = QLenSCV;  names{end+1} = 'QLenSCV';
end
if any(order == 3)
    vars{end+1} = QLenM3;   names{end+1} = 'QLenM3';
    vars{end+1} = QLenSkew; names{end+1} = 'QLenSkew';
end
MomentStationTable = table(vars{:}, 'VariableNames', names);
end
 
% =========================================================================
function tf = isExactMvaMethod(method)
% Methods whose means are the exact MVA recursion, so pfqn_sens_mom's analytic
% derivatives apply directly.
tf = any(strcmpi(method, {'default', 'mva', 'exact'}));
end
 
% =========================================================================
function mom = momentsByFiniteDifference(self, sn, queueIndices)
% Moments of the per-station totals from ANY solver and method, by central
% differences of that method's OWN mean queue lengths.
%
% The identity does not care how the means were obtained, so the solver's own
% method is used as a mean-value oracle. The parameter is y_i, a scaling of
% station i's whole demand column, which is Strelen's x_i; at y = 1 the
% y-derivatives are the scaled x-derivatives that (3.2) asks for. Because
% D(i,r) = visits(i,r)/rate(i,r), scaling station i's rates by 1/f scales its
% whole demand column by f, which is exactly the column perturbation wanted.
% NetworkStruct is a plain struct (see lang/NetworkStruct.m, a function, not a
% classdef), so copying it to perturb the rates cannot disturb the caller's sn.
%
% The oracle re-runs THIS solver, so every method of every product-form solver
% is covered, including the normalizing-constant methods of SolverNC (comom, ca,
% le, ...), which no hand-differentiated implementation reaches.
%
% Cost: 2*M extra solves. Accuracy: the moments inherit the accuracy of the
% method's means, and the second derivative inherits the usual h^2 truncation.
M = numel(queueIndices);
h = 1e-4;
qst = sn.nodeToStation(queueIndices);
m0 = solveTotals(self, sn, qst);
dm = zeros(M, M); d2m = zeros(M, 1);
for hcol = 1:M
    mp = solveTotals(self, scaleDemandColumn(sn, qst(hcol), 1+h), qst);
    mm = solveTotals(self, scaleDemandColumn(sn, qst(hcol), 1-h), qst);
    dm(:, hcol) = (mp - mm) / (2*h);
    d2m(hcol) = (mp(hcol) - 2*m0(hcol) + mm(hcol)) / h^2;
end
mom = packFiniteDifference(m0, dm, d2m);
end
 
% =========================================================================
function sn2 = scaleDemandColumn(sn, station, factor)
% Scale one station's whole demand column by FACTOR, via its service rates.
sn2 = sn;
sn2.rates(station, :) = sn.rates(station, :) / factor;
end
 
% =========================================================================
function m = solveTotals(self, sn, qst)
% Station totals under the solver's OWN method, whatever that is.
QN = solveMeansForStruct(self, sn);
m = sum(QN(qst, :), 2);
end
 
% =========================================================================
function mom = packFiniteDifference(m, dm, d2m)
% (3.2), applied to numerically obtained derivatives.
M = numel(m);
mom.m = m;
mom.dm = dm;
mom.d2m = d2m;
Cov = (dm + dm.') / 2;
mom.CovAsym = max(max(abs(dm - dm.')));
mom.Cov = Cov;
Var = zeros(M,1); M2 = zeros(M,1); M3 = zeros(M,1); Skew = zeros(M,1);
for i = 1:M
    Var(i) = dm(i,i);
    M2(i) = dm(i,i) + m(i)^2;
    M3(i) = d2m(i) + (1 + 3*m(i))*dm(i,i) + m(i)^3;
    mu3 = M3(i) - 3*m(i)*M2(i) + 2*m(i)^3;
    if Var(i) > 0
        Skew(i) = mu3 / Var(i)^1.5;
    else
        Skew(i) = NaN;
    end
end
mom.Var = Var; mom.M2 = M2; mom.M3 = M3; mom.Skew = Skew;
end
 
% =========================================================================
function tf = isLinearizerMethod(method)
% True for the Linearizer family of solver methods. Everything else maps to the
% exact recursion: the moment analysis has only these two algorithms, and the
% exact one is the right default for a method with no approximate counterpart.
tf = any(strcmpi(method, {'lin', 'amva.lin', 'egflin', 'gflin'}));
end
 
% =========================================================================
function order = validateMomentOrder(order, maxorder)
% ORDER is a set of moment orders. A scalar k is shorthand for 1:k, so that
% getMomentStationTable(2) means "up to the second moment" and not "the second
% moment alone"; a vector of two or more entries is taken literally.
%
% Consequence of MATLAB's isscalar: a one-element vector IS a scalar, so [2]
% takes the 1:k path and yields [1 2]. "The second moment alone" is therefore
% not expressible, which is deliberate: a variance with no mean beside it is not
% a useful table, and [2 3] remains available for the higher orders.
%
% Non-integers are rejected on BOTH paths. Rounding them silently would accept
% [1 2.5] as [1 3], i.e. answer a question that was not asked.
if ~isnumeric(order) || isempty(order) || any(~isfinite(order(:)))
    line_error(mfilename, sprintf('order must be an integer in 1..%d, or a vector of such integers.', maxorder));
end
if any(order(:) ~= round(order(:))) || any(order(:) < 1) || any(order(:) > maxorder)
    line_error(mfilename, sprintf('order must be an integer in 1..%d, or a vector of such integers.', maxorder));
end
if isscalar(order)
    order = 1:order;
    return;
end
order = unique(order(:)');
end