solver_mva_qsys_analyzer.m

function [Q,U,R,T,C,X,lG,runtime,totiter,actualmethod] = solver_mva_qsys_analyzer(sn, options)
% [Q,U,R,T,C,X,LG,RUNTIME,ITER,ACTUALMETHOD] = SOLVER_MVA_QSYS_ANALYZER(QN, OPTIONS)
 
% Copyright (c) 2012-2026, Imperial College London
% All rights reserved.
 
T0=tic;
M = sn.nstations;
K = sn.nclasses;
Q = zeros(M,K); U = zeros(M,K);
R = zeros(M,K); T = zeros(M,K);
C = zeros(1,K); X = zeros(M,K);
totiter = 1;
 
method = options.method;
actualmethod = method; % resolved algorithm name, reported in the solver banner
source_ist = sn.nodeToStation(sn.nodetype == NodeType.Source);
queue_ist = sn.nodeToStation(sn.nodetype == NodeType.Queue);
lambda = sn.rates(source_ist)*sn.visits{source_ist}(sn.stationToStateful(queue_ist));
k = sn.nservers(queue_ist);
mu = sn.rates(queue_ist);
ca = sqrt(sn.scv(source_ist));
cs = sqrt(sn.scv(queue_ist));
 
line_debug('MVA qsys analyzer starting: method=%s, lambda=%g, mu=%g, k=%d', method, lambda, mu, k);
 
% Check for BMAP arrivals (batch Markovian)
if sn.procid(source_ist) == ProcessType.BMAP
    line_debug('BMAP arrival process detected');
 
    % Check if service is exponential
    if sn.procid(queue_ist) == ProcessType.EXP
        line_debug('Service is exponential, using MX/M/1 queue model');
        % sn.proc holds {D0, D1, D_batch1, ..., D_batchK} (JAR MatrixCell
        % layout); batch rates are weighted by the stationary vector of
        % the underlying CTMC, as in BMAP.getBatchRates
        proc_src = sn.proc{source_ist};
        if iscell(proc_src) && ~isempty(proc_src) && iscell(proc_src{1})
            proc_src = proc_src{1};
        end
        pie_b = ctmc_solve(proc_src{1} + proc_src{2});
        nbatch = length(proc_src) - 2;
        batch_rates = zeros(1, nbatch);
        for bsize = 1:nbatch
            batch_rates(bsize) = pie_b * proc_src{2+bsize} * ones(size(proc_src{1},1),1);
        end
        total_rate = sum(batch_rates);
        lambda_batch = total_rate; % batch event rate
        if total_rate > 0
            E_X = sum((1:nbatch) .* batch_rates) / total_rate;
            E_X2 = sum(((1:nbatch).^2) .* batch_rates) / total_rate;
        else
            E_X = 1;
            E_X2 = 1;
        end
 
        [W, Wq, U_mxm1, Q_mxm1] = qsys_mxm1(lambda_batch, mu, E_X, E_X2);
        R(queue_ist,1) = W * sn.visits{1}(sn.stationToStateful(queue_ist));
        C(1,1) = R(queue_ist,1);
        X(queue_ist,1) = lambda_batch * E_X;  % Job arrival rate
        U(queue_ist,1) = U_mxm1;
        T(source_ist,1) = lambda_batch * E_X;
        T(queue_ist,1) = lambda_batch * E_X;
        Q(queue_ist,1) = Q_mxm1;
        lG = 0;
        runtime=toc(T0);
        actualmethod = 'mxm1';
        return;
    end
end
 
if strcmpi(method,'exact')
    if ca == 1 && cs == 1 && k==1
        method = 'mm1';
        line_debug('Exact method selected: M/M/1 (ca=1, cs=1, k=1)');
    elseif ca == 1 && cs == 1 && k>1
        method = 'mmk';
        line_debug('Exact method selected: M/M/k (ca=1, cs=1, k=%d)', k);
    elseif ca == 1 && k==1
        method = 'mg1';
        line_debug('Exact method selected: M/G/1 (ca=1, k=1)');
    elseif cs == 1 && k==1
        method = 'gm1';
        line_debug('Exact method selected: G/M/1 (cs=1, k=1)');
    else
        line_error(mfilename,'MVA exact method unavailable for this model.');
    end
end
 
switch method
    case 'default'
        if ca == 1 && cs == 1 && k == 1
            method = 'mm1';
            line_debug('Default method: using M/M/1 exact solution\n');
        elseif ca == 1 && cs == 1 && k > 1
            method = 'mmk';
            line_debug('Default method: using M/M/k exact solution (k=%d)\n', k);
        elseif ca == 1 && k == 1
            method = 'mg1';
            line_debug('Default method: using M/G/1 exact solution\n');
        elseif cs == 1 && k == 1
            method = 'gm1';
            line_debug('Default method: using G/M/1 exact solution\n');
        elseif k > 1
            method = 'gigk';
            line_debug('Default method: using G/G/k approximation (k=%d)\n', k);
        else
            method = 'gig1.klb';
            line_debug('Default method: using G/G/1 KLB approximation\n');
        end
end
 
switch method
    case 'mm1'
        line_debug('Using M/M/1 exact solution');
        R = qsys_mm1(lambda,mu);
    case 'mmk'
        line_debug('Using M/M/k exact solution (k=%d)', k);
        R = qsys_mmk(lambda,mu,k);
    case {'rqna'}
        line_debug('Using RQNA (robust queueing) single-queue solution');
        arvMAP = sn.proc{source_ist}{1};
        rho1 = lambda/mu;
        IaFun1 = @(x) map_count_idc(arvMAP, x);
        [~, W1] = qsys_gig1_rq(rho1, mu, cs^2, IaFun1);
        R = W1 + 1/mu;
    case {'mg1', 'mgi1'}  % verified
        line_debug('Using M/G/1 exact solution');
        R = qsys_mg1(lambda,mu,cs);
    case {'gigk'}
        line_debug('Using G/G/k approximation (k=%d)', k);
        R = qsys_gigk_approx(lambda,mu,ca,cs,k);
    case {'gigk.kingman_approx'}
        line_debug('Using G/G/k Kingman approximation (k=%d)', k);
        R = qsys_gigk_approx_kingman(lambda,mu,ca,cs,k);
    case 'gig1.kingman'  % verified
        line_debug('Using G/G/1 Kingman upper bound');
        R = qsys_gig1_ubnd_kingman(lambda,mu,ca,cs);
    case 'gig1.heyman'
        line_debug('Using G/G/1 Heyman approximation');
        R = qsys_gig1_approx_heyman(lambda,mu,ca,cs);
    case {'gig1', 'gig1.allen'}
        line_debug('Using G/G/1 Allen-Cunneen approximation');
        R = qsys_gig1_approx_allencunneen(lambda,mu,ca,cs);
    case 'gig1.kobayashi'
        line_debug('Using G/G/1 Kobayashi approximation');
        R = qsys_gig1_approx_kobayashi(lambda,mu,ca,cs);
    case 'gig1.klb'
        line_debug('Using G/G/1 KLB approximation');
        R = qsys_gig1_approx_klb(lambda,mu,ca,cs);
    case 'gig1.marchal' % verified
        line_debug('Using G/G/1 Marchal approximation');
        R = qsys_gig1_approx_marchal(lambda,mu,ca,cs);
    case {'gm1', 'gim1'}
        line_debug('Using G/M/1 exact solution');
        mu = sn.rates(queue_ist);
        % Prefer the exact PH/M/1 sigma-root via the source's PH representation.
        % Falls back to the LST-based fzero if that's available, else to the
        % two-moment sigma-root fit of qsys_gg1.
        R = [];
        % A deterministic arrival has no exact finite PH -- sn.proc stores an
        % Erlang-n moment approximation, so the PH path returns an approximate
        % sigma (D/M/1 0.64477 instead of the exact 0.6275). Skip PH for Det and
        % use the exact LST sigma-root below (Erlang/HyperExp/Coxian/APH keep PH,
        % which is exact for them).
        if sn.procid(source_ist) ~= ProcessType.DET
            try
                phPair = sn.proc{source_ist}{1};
                if iscell(phPair) && numel(phPair) >= 2
                    D0_ph = phPair{1};
                    D1_ph = phPair{2};
                    if size(D0_ph, 1) == size(D0_ph, 2) && all(size(D0_ph) == size(D1_ph))
                        pie_src = map_pie({D0_ph, D1_ph});
                        res_phm1 = qsys_phm1(pie_src, D0_ph, mu);
                        R = res_phm1.meanSojournTime;
                    end
                end
            catch
                R = [];
            end
        end
        if isempty(R)
            try
                LA = @(s) sn.lst{source_ist}{1}(s);
                sigma = fzero(@(x) LA(mu-mu*x)-x, 0.5);
                R = qsys_gm1(sigma, mu);
            catch
                R = qsys_gg1(lambda, mu, ca^2, 1);
            end
        end
    otherwise
        line_error(mfilename,'Unsupported method for a model with 1 station and 1 class.');
end
actualmethod = method;
 
Rscalar = R;  % save scalar from qsys function
% lambda is the effective (per-visit) arrival rate to the queue = external
% source rate times the queue visit ratio; Rscalar is the per-visit sojourn
% time. For feedback models the queue visit ratio Vq > 1, so the per-visit
% and per-job (residence) quantities must be kept distinct: RespT/QLen are
% per-visit, while ResidT (=RespT*Vq downstream) and the system response C
% are per-job. Source throughput is the external arrival rate, not lambda.
Vq = sn.visits{1}(sn.stationToStateful(queue_ist));
srcRate = sn.rates(source_ist);
R = zeros(M,K);
R(queue_ist,1) = Rscalar;              % per-visit response time
C(1,1) = Rscalar * Vq;                 % system (per-job) response time
X(queue_ist,1) = srcRate;              % system throughput = external arrival rate
U(queue_ist,1) = lambda/mu/k;
T(source_ist,1) = srcRate;             % source throughput = external arrival rate
T(queue_ist,1) = lambda;               % queue throughput = effective arrival rate
Q(queue_ist,1) = lambda * Rscalar;     % Little's law at the queue (per-visit)
lG = 0;
runtime=toc(T0);
end