pfqn_mvaoi.m

function [X, Qoi, Qli, Qdelay, Soi] = pfqn_mvaoi(Z, N, mu, Dli, options)
% [X, QOI, QLI, QDELAY, SOI] = PFQN_MVAOI(Z, N, MU, DLI, OPTIONS)
%
% Mean-value analysis of a closed product-form queueing network composed of an
% aggregated infinite-server (delay) node, any number of load-independent (LI)
% single-server product-form queues, and any number of order-independent (OI) /
% pass-and-swap stations with empty swap graph. This is the mean-value
% counterpart of PFQN_OI_NC and the marginal-distribution form PFQN_MVAOI_MARG:
% it returns the same exact per-class throughput and queue-lengths but WITHOUT
% computing any normalizing constant or joint marginal, using only mean
% quantities (throughputs, demands, queue-lengths) evaluated on shifted models.
% It is the composition-dependent generalization of the Conditional MVA (CMVA)
% of Casale, "A Note on Stable Flow-Equivalent Aggregation in Closed Networks"
% (QUESTA 2009), whose "third form" (rate depending on the full per-class
% occupancy vector) is realized here, extended to MULTIPLE OI stations by
% carrying one rate-shift vector s_i per OI station i.
%
% Throughout, r and s index job classes; i indexes OI stations; j indexes LI
% queues. For a single OI station and no LI queue the analysis recurs on the
% shift vector s_i (the OI occupancy already committed at the bottom of station
% i), Nn = N - s_i the jobs still to distribute:
%   Q^{(S)}(Nn) = sum_r U_r^{(S)}(Nn) ( e_r + Q^{(S+e_r@i)}(Nn - e_r) ),
%   U_r^{(S)}(Nn) = D_r^{(S)}(Nn) X_r^{(S)}(Nn)   (bottom-job utilization),
% with the class-r OI demand and throughput satisfying
%   D_r^{(S)}(Nn) = (1/mu_i(s_i+e_r)) rho_{i,r}^{(S)}(Nn-e_r),                     Nn_r = 1,
%   D_r^{(S)}(Nn) = [X_r^{(S)}(Nn-e_r)/X_r^{(S+e_r@i)}(Nn-e_r)] D_r^{(S)}(Nn-e_r), Nn_r >= 2,
%   rho_{i,r}^{(S)}(M) = rho_{i,r}^{(S)}(M-e_s) X_s^{(S)}(M)/X_s^{(S+e_r@i)}(M),  rho(0)=1, s ~= r,
% and X_r^{(S)}(Nn) closed by population conservation. With K OI stations the
% shift becomes a K x R matrix S (row i = s_i); each OI station keeps its own
% D^i, rho^i and Q^i recursions driven by the common throughput X^{(S)}(Nn), and
% the conservation identity aggregates every station's contribution:
%   Nn_r = X_r Z_r + sum_j Q^{(j)}_r + sum_i Q^{(i)}_r,
% where the LI queue Q^{(j)}_r = X_r D_{j,r} (1 + sum_s Q^{(j)}_s(Nn - e_r)) is
% the standard arrival-theorem term. States (S, Nn) are processed by increasing
% sum(Nn) so every reference lands at a strictly smaller free population.
%
% Parameters:
%   Z   - (1 x R) think-time demand vector of the aggregated delay node.
%   N   - (1 x R) closed population vector, finite.
%   mu  - cell array {mu_1,...,mu_K} of function handles; mu_i(n) returns the OI
%         total service rate of station i for per-class occupancy n (1 x R). A
%         bare function handle is accepted as the single-station shorthand.
%   Dli - (J x R) per-class demand matrix of the LI single-server queues
%         (D_{j,r} = V_{j,r}/rate_{j,r}); empty or omitted when J = 0.
%   options - solver options (optional, currently unused).
%
% Returns:
%   X      - (1 x R) per-class throughput X_r = G(N-e_r)/G(N).
%   Qoi    - (K x R) per-class mean queue-length at each OI station (row i).
%   Qli    - (J x R) per-class mean queue-length at each LI queue (row j).
%   Qdelay - (1 x R) per-class mean queue-length at the delay node (X.*Z).
%   Soi    - (K x R) per-class mean number of IN-SERVICE jobs at each OI station,
%            i.e. E[sir_r] with sir_r the count of class-r jobs receiving a
%            strictly positive rank rate (see PFQN_OI_INSVC); the utilization of
%            OI station i is Soi(i,r)/c_i. Unlike X/Qoi/Qli, which are pure
%            mean-value quantities, Soi is a distributional statistic and is
%            therefore obtained from the OI count marginal
%              pM_i(n|k) = (1/mu_i(n)) sum_r X_r(k) pM_i(n-e_r|k-e_r),
%              pM_i(0|k) = 1 - sum_{n ~= 0} pM_i(n|k),
%            which is assembled here from the zero-shift throughputs X^{(0)}(k)
%            already cached by the mean-value recursion above (no normalizing
%            constant is formed). It is only computed when requested.
%
% Copyright (c) 2012-2026, Imperial College London
% All rights reserved.
 
if nargin < 4
    Dli = [];
end
if nargin < 5
    options = struct(); %#ok<NASGU>
end
if isa(mu, 'function_handle')
    mu = {mu};
end
if ~iscell(mu) || isempty(mu)
    line_error(mfilename, 'mu must be a (nonempty) cell of OI rate handles.');
end
K = numel(mu);
for i = 1:K
    if ~isa(mu{i}, 'function_handle')
        line_error(mfilename, 'each mu{i} must be a function handle mu_i(n).');
    end
end
 
R = numel(N);
N = round(N(:)');
Z = Z(:)';
if any(~isfinite(N))
    line_error(mfilename, 'pfqn_mvaoi requires finite (closed) populations.');
end
if isempty(Dli)
    Dli = zeros(0, R);
end
J = size(Dli, 1);
I = eye(R);
 
% Caches keyed by the string of [S(:).' , Nn]. Reachable states only.
Xc  = containers.Map('KeyType','char','ValueType','any');   % X^{(S)}(Nn)   1xR
Qlc = containers.Map('KeyType','char','ValueType','any');   % Qli^{(S)}(Nn) JxR
Dc  = cell(1, K);   % Dc{i}: D_i^{(S)}(Nn)  1xR
Qc  = cell(1, K);   % Qc{i}: Q_i^{(S)}(Nn)  1xR
Rc  = cell(1, K);   % Rc{i}: rho_i^{(S)}(M) 1xR (NaN until computed)
for i = 1:K
    Dc{i} = containers.Map('KeyType','char','ValueType','any');
    Qc{i} = containers.Map('KeyType','char','ValueType','any');
    Rc{i} = containers.Map('KeyType','char','ValueType','any');
end
 
zeroS = zeros(K, R);
k0 = statekey(zeroS, zeros(1,R));
Xc(k0)  = zeros(1, R);
Qlc(k0) = zeros(J, R);
for i = 1:K
    Dc{i}(k0) = zeros(1, R);
    Qc{i}(k0) = zeros(1, R);
end
 
% Enumerate all reachable (S, Nn) with, per class r, sum_i s_i,r + Nn_r <= N_r.
% A slack bucket (the last row) absorbs the remaining N_r, so every (K+2) x R
% matrix produced by sprod has column sums N: rows 1..K are the OI shift vectors
% s_1..s_K, row K+1 is the free vector Nn, and the last (slack) row is dropped.
% States are then processed by increasing sum(Nn).
states = {};
[si, nmat, SS, DD] = sprod(K+2, N);
while si >= 0
    states{end+1} = nmat(1:K+1, :); %#ok<AGROW>
    [si, nmat] = sprod(si, SS, DD);
end
sums = cellfun(@(nm) sum(nm(K+1, :)), states);
[~, ord] = sort(sums);
states = states(ord);
 
for p = 1:numel(states)
    nm = states{p};
    S  = nm(1:K, :);
    Nn = nm(K+1, :);
    key = statekey(S, Nn);
    if sum(Nn) == 0
        % No free jobs: throughput and all queue-ahead terms vanish. This must
        % be stored for EVERY shift S (not only S = 0) because Qsub references
        % (S + e_s@i, Nn - e_s) which reaches such states when Nn = e_s.
        Xc(key)  = zeros(1, R);
        Qlc(key) = zeros(J, R);
        for i = 1:K
            Dc{i}(key) = zeros(1, R);
            Qc{i}(key) = zeros(1, R);
        end
        continue
    end
 
    % Per-OI-station class demands D_i and queue-ahead vectors Qsub_i.
    Dt   = zeros(K, R);         % Dt(i,r)
    Qsub = cell(1, K);          % Qsub{i}(s,r) = Q_i^{(S+e_s@i)}(Nn - e_s), comp r
    for i = 1:K
        Qsub{i} = zeros(R, R);
        for r = 1:R
            if Nn(r) == 0
                continue
            end
            if Nn(r) == 1
                mur = mu{i}(shiftrow(S, i, I(r,:)));
                if mur > 0
                    Dt(i,r) = (1/mur) * rho_fn(i, r, S, oner(Nn,r), mu, Xc, Rc, I, R);
                end
            else
                Sp = S; Sp(i,:) = Sp(i,:) + I(r,:);
                xs  = Xc(statekey(S,  oner(Nn,r)));
                xs2 = Xc(statekey(Sp, oner(Nn,r)));
                if xs2(r) > 0
                    Dprev = Dc{i}(statekey(S, oner(Nn,r)));
                    Dt(i,r) = (xs(r)/xs2(r)) * Dprev(r);
                end
            end
        end
        for s = 1:R
            if Nn(s) > 0
                Ss = S; Ss(i,:) = Ss(i,:) + I(s,:);
                Qsub{i}(s,:) = Qc{i}(statekey(Ss, oner(Nn,s)));
            end
        end
    end
 
    % LI-queue arrival-theorem coefficients beta_j,r and aggregate per class.
    betaLI = zeros(J, R);       % beta_j,r = D_j,r (1 + sum_s Qli_j,s(Nn - e_r))
    for r = 1:R
        if Nn(r) == 0
            continue
        end
        Qsub_li = Qlc(statekey(S, oner(Nn,r)));     % J x R
        for j = 1:J
            betaLI(j,r) = Dli(j,r) * (1 + sum(Qsub_li(j,:)));
        end
    end
 
    % Population conservation  A X = Nn  over classes with Nn_r > 0.
    idx = find(Nn > 0);
    m = numel(idx);
    A = zeros(m, m);
    for a = 1:m
        r = idx(a);
        for b = 1:m
            s = idx(b);
            if s == r
                val = Z(r) + sum(betaLI(:,r));
                for i = 1:K
                    val = val + Dt(i,r) * (1 + Qsub{i}(r,r));
                end
                A(a,b) = val;
            else
                val = 0;
                for i = 1:K
                    val = val + Dt(i,s) * Qsub{i}(s,r);
                end
                A(a,b) = val;
            end
        end
    end
    % Throughput closure. The OI queue recurrence couples X_s (s~=r) via the
    % off-diagonal A_{r,s}, so the population-conservation identity A X = Nn is
    % not a per-class Little's-law ratio as in the canonical CMVA. Rather than
    % solving the linear system, we decouple it with the product-form throughput
    % ratio identity (valid at fixed shift S)
    %   X_s(S,Nn)/X_r(S,Nn) = X_s(S,Nn-e_r)/X_r(S,Nn-e_s),
    % which yields the scalar per-class recurrence
    %   X_r = Nn_r / ( A_{r,r} + sum_{s~=r} A_{r,s} X_s(S,Nn-e_r)/X_r(S,Nn-e_s) )
    % from same-shift one-job-less throughputs already in the cache.
    Xk = zeros(1, R);
    for a = 1:m
        r = idx(a);
        denom = A(a,a);
        for b = 1:m
            if b ~= a
                s = idx(b);
                Xner = Xc(statekey(S, oner(Nn,r)));   % X(S, Nn-e_r)
                Xnes = Xc(statekey(S, oner(Nn,s)));   % X(S, Nn-e_s)
                if Xnes(r) > 0
                    denom = denom + A(a,b) * (Xner(s) / Xnes(r));
                end
            end
        end
        if denom > 0
            Xk(r) = Nn(r) / denom;
        end
    end
 
    % Assemble OI-station queues Q_i and LI-queue queues Qli.
    Qk_li = zeros(J, R);
    for r = 1:R
        if Nn(r) == 0
            continue
        end
        Qk_li(:,r) = Xk(r) * betaLI(:,r);
    end
    for i = 1:K
        U = Dt(i,:) .* Xk;               % 1 x R bottom-job utilization
        Qi = zeros(1, R);
        for r = 1:R
            Qi(r) = U(r) + U * Qsub{i}(:,r);
        end
        Qc{i}(key) = Qi;
        Dc{i}(key) = Dt(i,:);
    end
    Xc(key)  = Xk;
    Qlc(key) = Qk_li;
end
 
keyN = statekey(zeroS, N);
X = Xc(keyN);
Qoi = zeros(K, R);
for i = 1:K
    Qoi(i,:) = Qc{i}(keyN);
end
Qli = Qlc(keyN);
Qdelay = X .* Z;
 
if nargout >= 5
    Soi = oi_insvc_means(Z, N, mu, Xc, zeroS, K, R);
end
end
 
% =========================================================================
function Soi = oi_insvc_means(Z, N, mu, Xc, zeroS, K, R)
% Mean number of in-service jobs per class at each OI station, E[sir_r], from
% the OI count marginal pM_i(n|k) built on the cached zero-shift throughputs
% X^{(0)}(k). Exact because in product form
% pM_i(n|k) = Phi_i(n) G_{-i}(k-n)/G(k) and X_r(k) = G(k-e_r)/G(k), so the
% recursion below reproduces the balanced-fairness identity for Phi_i.
shp = N + 1;
stride = ones(1, R);
for d = 2:R
    stride(d) = stride(d-1) * shp(d-1);
end
total = prod(shp);
subs = zeros(total, R);
for i = 1:total
    li = i - 1;
    for d = 1:R
        subs(i, d) = mod(li, shp(d));
        li = floor(li / shp(d));
    end
end
[~, ord] = sort(sum(subs, 2));      % process populations by increasing size
 
% Cache X^{(0)}(k) over the lattice.
Xk = zeros(total, R);
for i = 1:total
    Xk(i,:) = Xc(statekey(zeroS, subs(i,:)));
end
 
Soi = zeros(K, R);
for m = 1:K
    gm = pfqn_oi_insvc(mu{m}, N);        % E[sir_r | n] over the lattice
    muv = zeros(total, 1);
    for i = 1:total
        if sum(subs(i,:)) > 0
            muv(i) = mu{m}(subs(i,:));
        end
    end
    % pMv(a,b) = pM_m(n_a | k_b), filled for n_a <= k_b by increasing sum(k).
    pMv = zeros(total, total);
    pMv(1,1) = 1;                        % pM(0|0) = 1
    for bb = 1:total
        b = ord(bb);
        k = subs(b, :);
        if sum(k) == 0
            continue
        end
        acc0 = 0;
        for aa = 1:total
            a = ord(aa);
            n = subs(a, :);
            if sum(n) == 0 || any(n > k)
                continue
            end
            if muv(a) <= 0
                continue
            end
            acc = 0;
            for r = 1:R
                if n(r) > 0
                    acc = acc + Xk(b, r) * pMv(a - stride(r), b - stride(r));
                end
            end
            pMv(a, b) = acc / muv(a);
            acc0 = acc0 + pMv(a, b);
        end
        pMv(1, b) = 1 - acc0;            % empty-state probability by complement
    end
    idxN = 1 + sum(N .* stride);
    for r = 1:R
        Soi(m, r) = pMv(:, idxN)' * gm(:, r);
    end
end
end
 
% =========================================================================
% Helper functions
% =========================================================================
 
function rv = rho_fn(i, r, S, M, mu, Xc, Rc, I, R)
% rho_{i,r}^{(S)}(M): recursion on M (class r held fixed, M_r stays 0). NOTE all
% local names are r-prefixed to avoid corrupting the recursion via shared
% storage in nested calls.
rkey = statekey(S, M);
if isKey(Rc{i}, rkey)
    rrow = Rc{i}(rkey);
    if ~isnan(rrow(r))
        rv = rrow(r);
        return
    end
else
    rrow = nan(1, R);
end
if sum(M) == 0
    rv = 1.0;
    rrow(r) = rv;
    Rc{i}(rkey) = rrow;
    return
end
rs = -1;
for rss = 1:R
    if rss ~= r && M(rss) > 0
        rs = rss;
        break
    end
end
rxu  = Xc(statekey(S, M));
rSp  = S; rSp(i,:) = rSp(i,:) + I(r,:);
rxu2 = Xc(statekey(rSp, M));
rratio = 0.0;
if rxu2(rs) > 0
    rratio = rxu(rs) / rxu2(rs);
end
rv = rho_fn(i, r, S, oner(M, rs), mu, Xc, Rc, I, R) * rratio;
rrow(r) = rv;
Rc{i}(rkey) = rrow;
end
 
function row = shiftrow(S, i, e)
% Return the per-class occupancy row s_i + e for OI station i.
row = S(i,:) + e;
end
 
function v = oner(v, r)
% Decrement entry r of v by one.
v(r) = v(r) - 1;
end
 
function key = statekey(S, Nn)
% String key for the lattice state (S, Nn); S is K x R, Nn is 1 x R.
key = sprintf('%d_', [reshape(S.', 1, []), Nn]);
end