solver_mam_qsys.m

function [QN, UN, RN, TN, piAgg, GM] = solver_mam_qsys(arrival, service, options)
% SOLVER_MAM_QSYS Unified solver for MAP/BMAP queueing systems with single server
%
% [QN, UN, RN, TN, PI, GM] = SOLVER_MAM_QSYS(ARRIVAL, SERVICE) solves a queueing
% system with Markovian arrival and service processes. Automatically detects
% whether arrivals or service are batched and applies the appropriate
% matrix-analytic method:
%   - BMAP/MAP/1: M/G/1 type analysis (batch arrivals, single service)
%   - MAP/BMAP/1: GI/M/1 type analysis (single arrivals, batch service)
%
% Input:
%   arrival - Arrival process: MAP object, BMAP object, or cell array
%             MAP: {D0, D1} or MAP object
%             BMAP: {D0, D1, D2, ..., DK} or BMAP object
%   service - Service process: MAP object, BMAP object, or cell array
%             MAP: {S0, S1} or MAP object
%             BMAP: {S0, S1, S2, ..., SK} or BMAP object
%   options - (optional) Options structure with fields:
%             .nMoments: number of queue length moments to compute (default: 3)
%
% Output:
%   QN    - Mean queue length E[N]
%   UN    - Server utilization rho
%   RN    - Mean response time E[R]
%   TN    - Throughput (arrival rate)
%   piAgg - Aggregated stationary probabilities
%   GM    - G matrix (M/G/1 type) or R matrix (GI/M/1 type)
%
% Queue Types:
%   BMAP/MAP/1 - Batch arrivals, single service (M/G/1 type)
%     Level increases by batch size (1, 2, 3, ...), decreases by 1
%   MAP/BMAP/1 - Single arrivals, batch service (GI/M/1 type)
%     Level increases by 1, decreases by batch size (1, 2, 3, ...)
%
% References:
%   Stathopoulos, Riska, Hua, Smirni: ETAQA algorithm
%   MAMSolver: www.cs.wm.edu/MAMSolver
%
% See also: MAP, BMAP
%
% Copyright (c) 2012-2026, Imperial College London
% All rights reserved.
 
%% Input validation and setup
if nargin < 3
    options = struct();
end
if ~isfield(options, 'nMoments')
    options.nMoments = 3;
end
 
%% Determine arrival type (MAP or BMAP)
[arrivalMatrices, arrivalIsBatch, ma] = parseProcess(arrival, 'arrival');
 
%% Determine service type (MAP or BMAP)
[serviceMatrices, serviceIsBatch, ms] = parseProcess(service, 'service');
 
%% Route to appropriate solver based on queue type
if arrivalIsBatch && ~serviceIsBatch
    % BMAP/MAP/1 queue - M/G/1 type analysis
    [QN, UN, RN, TN, piAgg, GM] = solveBMAPMAP1(arrivalMatrices, serviceMatrices, ma, ms, options);
elseif ~arrivalIsBatch && serviceIsBatch
    % MAP/BMAP/1 queue - GI/M/1 type analysis
    [QN, UN, RN, TN, piAgg, GM] = solveMAPBMAP1(arrivalMatrices, serviceMatrices, ma, ms, options);
elseif ~arrivalIsBatch && ~serviceIsBatch
    % MAP/MAP/1 queue - treat as BMAP/MAP/1 with batch size 1
    [QN, UN, RN, TN, piAgg, GM] = solveBMAPMAP1(arrivalMatrices, serviceMatrices, ma, ms, options);
else
    % BMAP/BMAP/1 - both batched, use finite CTMC truncation
    [QN, UN, RN, TN, piAgg, GM] = solveBMAPBMAP1(arrivalMatrices, serviceMatrices, ma, ms, options);
end
 
end
 
%% Helper function to parse arrival/service process
function [matrices, isBatch, nPhases] = parseProcess(proc, procType)
    if isa(proc, 'BMAP')
        K = proc.getMaxBatchSize();
        matrices = cell(1, K + 1);
        matrices{1} = proc.D(0);  % D0
        for k = 1:K
            matrices{k+1} = proc.D(1, k);  % Dk (batch size k)
        end
        isBatch = K > 1;
        nPhases = size(matrices{1}, 1);
    elseif isa(proc, 'MAP')
        matrices = {proc.D(0), proc.D(1)};
        isBatch = false;
        nPhases = size(matrices{1}, 1);
    elseif iscell(proc)
        matrices = proc;
        isBatch = length(proc) > 2;  % More than {D0, D1} means batch process
        nPhases = size(proc{1}, 1);
    else
        error('%s must be a MAP object, BMAP object, or cell array', procType);
    end
 
    % Validate dimensions
    for k = 1:length(matrices)
        if size(matrices{k}, 1) ~= nPhases || size(matrices{k}, 2) ~= nPhases
            error('All %s matrices must be %dx%d', procType, nPhases, nPhases);
        end
    end
end
 
%% BMAP/MAP/1 solver using M/G/1 type analysis
function [QN, UN, RN, TN, pi, G] = solveBMAPMAP1(D, S, ma, ms, options)
% Solves BMAP/MAP/1 using M/G/1 type matrix-analytic methods
%
% M/G/1 type structure (level can jump UP by any amount, DOWN by exactly 1):
%   A0 = I_ma \otimes S1           (service completion, level -1)
%   A1 = D0 \otimes I_ms + I_ma \otimes S0  (phase changes, level 0)
%   A_{k+1} = D_k \otimes I_ms     (batch arrival size k, level +k)
 
K = length(D) - 1;  % max batch size for arrivals
m = ma * ms;        % Combined phases per level
 
% Extract MAP service matrices
S0 = S{1};
S1 = S{2};
 
%% Construct M/G/1 type matrices using Kronecker products
I_ma = eye(ma);
I_ms = eye(ms);
 
A = zeros(m, m * (K + 2));
 
% A0 = I_ma \otimes S1 (service completion)
A0 = kron(I_ma, S1);
A(:, 1:m) = A0;
 
% A1 = D0 \otimes I_ms + I_ma \otimes S0 (phase changes only)
A1 = kron(D{1}, I_ms) + kron(I_ma, S0);
A(:, m+1:2*m) = A1;
 
% A_{k+1} = D_k \otimes I_ms for k >= 1 (batch arrivals)
for k = 1:K
    Ak = kron(D{k+1}, I_ms);
    A(:, (k+1)*m + 1 : (k+2)*m) = Ak;
end
 
% B = [B0, B1, B2, ..., BK] for level 0 (empty queue)
B = zeros(m, m * (K + 1));
 
% B0: At level 0, no service (add S1 back as self-loop)
B0 = kron(D{1}, I_ms) + kron(I_ma, S0 + S1);
B(:, 1:m) = B0;
 
% B_k = D_k \otimes I_ms for k >= 1 (batch arrivals from empty queue)
for k = 1:K
    Bk = kron(D{k+1}, I_ms);
    B(:, k*m + 1 : (k+1)*m) = Bk;
end
 
%% Verify stability condition
D0 = D{1};
D1_total = zeros(ma, ma);
for k = 1:K
    D1_total = D1_total + D{k+1};
end
 
% Stationary distribution of BMAP
pi_bmap = map_prob({D0, D1_total});
e_ma = ones(ma, 1);
 
% Total customer arrival rate: sum_k (k * π * Dk * e)
lambda_total = 0;
for k = 1:K
    rate_k = pi_bmap * D{k+1} * e_ma;
    lambda_total = lambda_total + k * rate_k;
end
 
% Compute service rate from MAP
pi_map = map_prob({S0, S1});
e_ms = ones(ms, 1);
mu = pi_map * S1 * e_ms;  % Service rate
 
% Utilization
rho = lambda_total / mu;
 
if rho >= 1
    warning('solver_mam_qsys:Unstable', ...
        'System is unstable (rho = %.4f >= 1). Results may be invalid.', rho);
end
 
%% Compute G matrix using MG1_G_ETAQA
G = MG1_G_ETAQA(A);
 
%% Compute stationary probabilities using MG1_pi_ETAQA
pi = MG1_pi_ETAQA(B, A, G);
 
%% Compute queue length moments using MG1_qlen_ETAQA
qlen_moments = zeros(1, options.nMoments);
for n = 1:options.nMoments
    qlen_moments(n) = MG1_qlen_ETAQA(B, A, pi, n);
end
 
%% Compute performance metrics
QN = qlen_moments(1);  % Mean queue length (E[N])
UN = rho;              % Utilization
TN = lambda_total;     % Throughput (total arrival rate)
RN = QN / TN;          % Mean response time (Little's law: E[R] = E[N] / λ)
 
end
 
%% MAP/BMAP/1 solver using GI/M/1 type analysis
function [QN, UN, RN, TN, piAgg, R] = solveMAPBMAP1(C, D, ma, ms, options)
% Solves MAP/BMAP/1 using GI/M/1 type matrix-analytic methods
%
% GI/M/1 type structure (level jumps UP by 1, DOWN by any amount):
%   A0 = C1 \otimes I_ms           (MAP arrival, level +1)
%   A1 = C0 \otimes I_ms + I_ma \otimes D0  (phase changes, level 0)
%   A_{k+1} = I_ma \otimes D_k     (batch size k service, level -k)
 
C0 = C{1};
C1 = C{2};
K = length(D) - 1;  % max batch size for service
m = ma * ms;        % Combined phases per level
 
%% Validate inputs
% Check MAP is valid generator
rowSumsC = sum(C0 + C1, 2);
if max(abs(rowSumsC)) > 1e-10
    error('MAP matrices C0 + C1 must have zero row sums');
end
 
% Check BMAP is valid generator
DTotal = D{1};
for k = 1:K
    DTotal = DTotal + D{k+1};
end
rowSumsD = sum(DTotal, 2);
if max(abs(rowSumsD)) > 1e-10
    error('BMAP matrices D0 + D1 + ... + DK must have zero row sums');
end
 
%% Compute arrival and service rates
% Arrival rate from MAP
piC = map_prob({C0, C1});
lambdaArr = piC * C1 * ones(ma, 1);
 
% Service rate from BMAP (total customers served per unit time)
piD = map_prob({D{1}, DTotal - D{1}});  % D0, sum of D1...DK
muTotal = 0;
for k = 1:K
    muTotal = muTotal + k * (piD * D{k+1} * ones(ms, 1));
end
 
% Utilization
rho = lambdaArr / muTotal;
 
if rho >= 1.0
    warning('solver_mam_qsys:Unstable', ...
        'System is unstable (rho = %.4f >= 1). Results may be invalid.', rho);
end
 
%% Construct A matrix for GI/M/1-type (repeating part)
A = zeros(m * (K + 2), m);
 
% A0: MAP arrival (level +1)
A0 = kron(C1, eye(ms));
A(1:m, :) = A0;
 
% A1: Phase changes only (level 0)
A1 = kron(C0, eye(ms)) + kron(eye(ma), D{1});
A(m+1:2*m, :) = A1;
 
% A_{k+1}: Batch service of size k (level -k)
for k = 1:K
    Ak = kron(eye(ma), D{k+1});
    A((k+1)*m+1:(k+2)*m, :) = Ak;
end
 
%% Construct B matrix for boundary levels
B = zeros(m * (K + 1) + m, m);
 
% B1 (level 0 -> level 0): Phase changes only, no service from empty
B1 = kron(C0, eye(ms)) + kron(eye(ma), D{1});
% Add service transitions that would go to negative levels back to level 0
for k = 1:K
    B1 = B1 + kron(eye(ma), D{k+1});
end
B(1:m, :) = B1;
 
% B2, B3, ..., B_{K+1}: Transitions to level 0 from levels 1, 2, ..., K
for j = 1:K
    Bj = zeros(m, m);
    for k = j:K
        Bj = Bj + kron(eye(ma), D{k+1});
    end
    B(m + (j-1)*m + 1 : m + j*m, :) = Bj;
end
 
%% Compute R matrix using GI/M/1 ETAQA
R = GIM1_R_ETAQA(A);
 
%% Compute stationary probabilities using GI/M/1 ETAQA
piAgg = GIM1_pi_ETAQA(B, A, R, 'Boundary', A0);
 
%% Compute performance metrics
% Mean queue length using ETAQA moments
QN = GIM1_qlen_ETAQA(B, A, R, piAgg, 1, 'Boundary', A0);
 
% Throughput equals arrival rate (in steady state)
TN = lambdaArr;
 
% Utilization
UN = rho;
 
% Mean response time via Little's Law
RN = QN / TN;
 
end
 
%% BMAP/BMAP/1 solver using finite CTMC truncation
function [QN, UN, RN, TN, piAgg, G] = solveBMAPBMAP1(D, S, ma, ms, options)
% Solves BMAP/BMAP/1 using finite CTMC truncation
%
% When both arrivals and service are batched, the process is neither
% M/G/1 type nor GI/M/1 type. We construct a finite truncated CTMC
% over the Kronecker product of arrival and service phases per level.
%
% State: (level, arrival_phase, service_phase)
%   level = number of customers in system (0, 1, ..., L)
%   arrival_phase = 1, ..., ma
%   service_phase = 1, ..., ms
 
Ka = length(D) - 1;  % max batch size for arrivals
Ks = length(S) - 1;  % max batch size for service
m = ma * ms;         % combined phases per level
 
%% Compute arrival and service rates for stability check
D1_total = zeros(ma, ma);
for k = 1:Ka
    D1_total = D1_total + D{k+1};
end
pi_bmap = map_prob({D{1}, D1_total});
e_ma = ones(ma, 1);
lambda_total = 0;
for k = 1:Ka
    lambda_total = lambda_total + k * (pi_bmap * D{k+1} * e_ma);
end
 
S1_total = zeros(ms, ms);
for k = 1:Ks
    S1_total = S1_total + S{k+1};
end
pi_smap = map_prob({S{1}, S1_total});
e_ms = ones(ms, 1);
mu_total = 0;
for k = 1:Ks
    mu_total = mu_total + k * (pi_smap * S{k+1} * e_ms);
end
 
rho = lambda_total / mu_total;
 
if rho >= 1
    warning('solver_mam_qsys:Unstable', ...
        'System is unstable (rho = %.4f >= 1). Results may be invalid.', rho);
end
 
%% Choose truncation level based on utilization
L = max(100, ceil(20 / (1 - min(rho, 0.99))));
L = min(L, 500);  % cap for memory
 
%% Build generator matrix Q as sparse (L+1)*m x (L+1)*m
N_states = (L + 1) * m;
I_ma = speye(ma);
I_ms = speye(ms);
 
% Pre-compute Kronecker products
arrKron = cell(1, Ka);
for k = 1:Ka
    arrKron{k} = kron(sparse(D{k+1}), I_ms);
end
svcKron = cell(1, Ks);
for k = 1:Ks
    svcKron{k} = kron(I_ma, sparse(S{k+1}));
end
phaseKron = kron(sparse(D{1}), I_ms) + kron(I_ma, sparse(S{1}));
 
% Estimate non-zeros and pre-allocate
nnz_est = (L+1) * m * m * (1 + Ka + Ks);
ii = zeros(nnz_est, 1);
jj = zeros(nnz_est, 1);
vv = zeros(nnz_est, 1);
cnt = 0;
 
for lev = 0:L
    base = lev * m;
 
    % Internal phase transitions: D_0 x I_ms + I_ma x S_0
    [ri, ci, vi] = find(phaseKron);
    n_entries = length(ri);
    ii(cnt+1:cnt+n_entries) = base + ri;
    jj(cnt+1:cnt+n_entries) = base + ci;
    vv(cnt+1:cnt+n_entries) = vi;
    cnt = cnt + n_entries;
 
    % Arrivals: batch size k, level -> min(level+k, L)
    for k = 1:Ka
        new_lev = min(lev + k, L);
        dest = new_lev * m;
        [ri, ci, vi] = find(arrKron{k});
        n_entries = length(ri);
        ii(cnt+1:cnt+n_entries) = base + ri;
        jj(cnt+1:cnt+n_entries) = dest + ci;
        vv(cnt+1:cnt+n_entries) = vi;
        cnt = cnt + n_entries;
    end
 
    if lev > 0
        % Service: batch size j, level -> max(level-j, 0)
        for j = 1:Ks
            new_lev = max(lev - j, 0);
            dest = new_lev * m;
            [ri, ci, vi] = find(svcKron{j});
            n_entries = length(ri);
            ii(cnt+1:cnt+n_entries) = base + ri;
            jj(cnt+1:cnt+n_entries) = dest + ci;
            vv(cnt+1:cnt+n_entries) = vi;
            cnt = cnt + n_entries;
        end
    else
        % Level 0: server idle, service completions stay at level 0
        for j = 1:Ks
            [ri, ci, vi] = find(svcKron{j});
            n_entries = length(ri);
            ii(cnt+1:cnt+n_entries) = base + ri;
            jj(cnt+1:cnt+n_entries) = base + ci;
            vv(cnt+1:cnt+n_entries) = vi;
            cnt = cnt + n_entries;
        end
    end
end
 
% Trim and assemble sparse Q
ii = ii(1:cnt);
jj = jj(1:cnt);
vv = vv(1:cnt);
Q = sparse(ii, jj, vv, N_states, N_states);
 
% Fix diagonal: Q(i,i) = -sum of off-diagonal entries in row i
d = full(sum(Q, 2));
Q = Q - spdiags(d, 0, N_states, N_states);
 
%% Solve for steady-state distribution
% Solve pi * Q = 0 with pi * e = 1
% Equivalent to Q' * pi' = 0
QT = Q';
% Replace last equation with normalization constraint
QT(N_states, :) = 1;
b = sparse(N_states, 1);
b(N_states) = 1;
pi = full(QT \ b)';
 
% Ensure non-negative
pi = max(pi, 0);
pi = pi / sum(pi);
 
%% Compute performance metrics
% Mean queue length
levels = 0:L;
QN = levels * reshape(sum(reshape(pi, m, L+1), 1), L+1, 1);
 
% Utilization and throughput
UN = rho;
TN = lambda_total;
 
% Mean response time via Little's law
RN = QN / TN;
 
% Aggregate probabilities by level
piAgg = sum(reshape(pi, m, L+1), 1);
 
G = [];  % No G matrix for truncation approach
 
end