qsys_mgck_sync.h

Go to the documentation of this file.

/*
 * Copyright (c) 2012-2026, QORE Lab, Imperial College London
 * All rights reserved.
 */
#ifndef LINE_API_QSYS_QSYS_MGCK_SYNC_H
#define LINE_API_QSYS_QSYS_MGCK_SYNC_H
 
/**
 * @file
 * @ingroup api_qsys
 * Synchronous approximation of the M/G/c/K queue with dual bounds (German,
 * EPEW 2025).
 *
 * Port of matlab/src/api/qsys/qsys_mgck_sync.m. The concurrent service clocks
 * of the c servers are replaced by synchronized ones, so that at most one
 * non-Markovian activity is enabled at any time and the process is solved
 * exactly by Markov renewal theory through its embedded Markov chain. Every
 * customer receives either its true service time or a shorter one (variant
 * "faster") or a longer one ("slower"), which brackets the M/G/c/K measures:
 * the faster model gives a lower bound on the mean number in system, on the
 * response time and on P(full), and an upper bound on P(idle) and the
 * throughput; the slower model the opposite bounds. The argument is
 * customer-wise and holds without losses; with blocking the bounds are
 * observed, not proven.
 *
 * MODEL WITHOUT EXPANSION (nphases = 0), state i = jobs in system. All jobs in
 * service share one clock. "faster": an arrival that finds a free server joins
 * the running batch and completes with it. "slower": such an arrival preempts
 * the clock, so the whole batch restarts with it. On expiry min(i,c) jobs
 * depart and the next min(i-c,c) start together. For c = 1 both variants are
 * the exact M/G/1/K queue.
 *
 * MODEL WITH EXPANSION (nphases = N >= 1, deterministic service only), state
 * (a, q, B): a = 1 when the owner server is busy, q the waiting jobs, B the
 * sorted multiset of phases completed by the other busy servers. While the
 * owner is busy a tick of length d/N advances every phase and the owner
 * completes after N ticks; a newcomer to another server starts at phase 0
 * ("faster") or -1 ("slower"). With the owner idle the ticks continue; an
 * arrival takes the owner role and preempts the tick, the other phases being
 * rounded up ("faster") or kept ("slower"). An arrival to the empty system
 * takes the owner role with probability 1/c.
 *
 * Embedded chain (German 2000, cascaded embedding for the expansion):
 *   P = I_E - diag^{-1}(Q_E) Q_E + Omega Delta + Psi Qbar
 *   C = -diag^{-1}(Q_E) + Psi,   pi = u C/(u C e),  u = u P
 * with Omega = I_G E[exp(Q_G S)] and Psi = I_G E[int_0^S exp(Q_G t) dt], both
 * evaluated by uniformization of Q_G at rate lambda from the service-law
 * coefficients of qsys_svclaw_coeffs.h.
 *
 * A truncated uniformization series (residual above 1e-6) warns and continues,
 * as the reference does.
 *
 * ARITHMETIC. Gated on num_traits<T>::has_transcendental: the coefficients need
 * exp and the series truncation is a tolerance.
 *
 * References:
 * [1] R. German, On Synchronous Approximations of Non-Markovian Processes with
 *     Dual Bounds, EPEW 2025.
 * [2] R. German, Performance Analysis of Communication Systems with
 *     Non-Markovian Stochastic Petri Nets, Wiley, 2000.
 */
 
#include <algorithm>
#include <cctype>
#include <cstddef>
#include <cstdio>
#include <iostream>
#include <map>
#include <string>
#include <vector>
 
#include "line/api/qsys/qsys_mapg1k.h"
#include "line/api/qsys/qsys_svclaw_coeffs.h"
#include "line/num/number.h"
#include "line/util/error.h"
#include "line/util/line_console.h"
#include "line/util/linalg.h"
#include "line/util/lu.h"
#include "line/util/matrix.h"
 
namespace line {
namespace qsys {
 
/** Return value of qsys_mgck_sync, mirroring the MATLAB result struct. */
template <class T>
struct MgckSyncResult {
    std::vector<T> queueLengthDist;  ///< P(n jobs in system), n = 0..K
    T p0;                            ///< P(idle)
    T pK;                            ///< P(full)
    T meanQueueLength;               ///< mean number in system
    T meanQueueLengthQ;              ///< mean number waiting
    T throughput;                    ///< lambda (1 - pK)
    T lossProbability;               ///< pK (PASTA)
    T meanSojournTime;               ///< meanQueueLength / throughput
    T meanWaitingTime;               ///< meanQueueLengthQ / throughput
    T utilization;                   ///< mean busy servers / c in the synchronized model
    std::string variant;             ///< "faster" or "slower"
    std::size_t nphases = 0;
    std::size_t stateSpaceSize = 0;
    std::string analyzer;
};
 
namespace mgck_sync_detail {
 
/** Row-masked copy: rows with keep[i] false are zeroed, i.e. diag(keep) M. */
template <class T>
Matrix<T> rowmask(const Matrix<T>& M, const std::vector<bool>& keep) {
    Matrix<T> R = M;
    const T zero = num_traits<T>::from_int(0);
    for (std::size_t i = 0; i < M.rows(); ++i)
        if (!keep[i])
            for (std::size_t j = 0; j < M.cols(); ++j) R(i, j) = zero;
    return R;
}
 
/** Column-masked copy: M diag(keep). */
template <class T>
Matrix<T> colmask(const Matrix<T>& M, const std::vector<bool>& keep) {
    Matrix<T> R = M;
    const T zero = num_traits<T>::from_int(0);
    for (std::size_t j = 0; j < M.cols(); ++j)
        if (!keep[j])
            for (std::size_t i = 0; i < M.rows(); ++i) R(i, j) = zero;
    return R;
}
 
template <class T>
Matrix<T> matadd(const Matrix<T>& A, const Matrix<T>& B) {
    Matrix<T> C = A;
    for (std::size_t i = 0; i < A.rows(); ++i)
        for (std::size_t j = 0; j < A.cols(); ++j) C(i, j) += B(i, j);
    return C;
}
 
/**
 * i_unif: Om = E[exp(QG S)], Ps = E[int_0^S exp(QG t) dt] from c_n at rate
 * theta, with d_n = (1 - sum_{k<=n} c_k)/theta.
 */
template <class T>
void unif(const Matrix<T>& QG, const std::vector<T>& cn, const T& theta, Matrix<T>& Om,
          Matrix<T>& Ps) {
    const std::size_t dim = QG.rows();
    const T zero = num_traits<T>::from_int(0);
    Matrix<T> P = eye<T>(dim);
    for (std::size_t i = 0; i < dim; ++i)
        for (std::size_t j = 0; j < dim; ++j) P(i, j) += QG(i, j) / theta;
    const std::size_t n = cn.size();
    // tailc(k) = sum_{m>=k} c_m; dn(k) = tailc(k+1)/theta, dn(last) = 0
    std::vector<T> dn(n, zero);
    {
        T tail = zero;
        for (std::size_t k = n; k-- > 0;) {
            dn[k] = tail / theta;
            tail += cn[k];
        }
    }
    Om = Matrix<T>(dim, dim, zero);
    Ps = Matrix<T>(dim, dim, zero);
    Matrix<T> Pn = eye<T>(dim);
    for (std::size_t k = 0; k < n; ++k) {
        for (std::size_t i = 0; i < dim; ++i)
            for (std::size_t j = 0; j < dim; ++j) {
                const T v = Pn(i, j);
                if (v == zero) continue;
                Om(i, j) += cn[k] * v;
                Ps(i, j) += dn[k] * v;
            }
        if (k + 1 < n) Pn = matmul(Pn, P);
    }
}
 
/** i_emc: stationary law of the embedded chain at regenerations, converted to time. */
template <class T>
std::vector<T> emc(const Matrix<T>& Q, const std::vector<bool>& inE, const Matrix<T>& OmDelta,
                   const Matrix<T>& Psi, const Matrix<T>& Qbar) {
    const std::size_t dim = Q.rows();
    const T zero = num_traits<T>::from_int(0), one = num_traits<T>::from_int(1);
    const Matrix<T> QE = rowmask(Q, inE);
    std::vector<T> dinv(dim, zero);
    for (std::size_t i = 0; i < dim; ++i)
        if (QE(i, i) != zero) dinv[i] = one / QE(i, i);
    Matrix<T> P = matadd(OmDelta, matmul(Psi, Qbar));
    for (std::size_t i = 0; i < dim; ++i) {
        if (inE[i]) P(i, i) += one;
        for (std::size_t j = 0; j < dim; ++j) P(i, j) -= dinv[i] * QE(i, j);
    }
    Matrix<T> C = Psi;
    for (std::size_t i = 0; i < dim; ++i) C(i, i) -= dinv[i];
    T rowdev = zero;
    for (std::size_t i = 0; i < dim; ++i) {
        T s = zero;
        for (std::size_t j = 0; j < dim; ++j) s += P(i, j);
        const T d = num_abs(T(s - one));
        if (d > rowdev) rowdev = d;
    }
    if (rowdev > num_traits<T>::from_double(1e-8))
        throw NumericError(ErrorKind{"NonConvergence"},
            "qsys_mgck_sync: embedded chain rows deviate from 1 by " +
            num_traits<T>::to_string(rowdev) + "; raise nmax");
    // A = (P - I)^T with its last row replaced by (C e)^T, b = e_dim
    Matrix<T> A(dim, dim, zero);
    for (std::size_t i = 0; i < dim; ++i)
        for (std::size_t j = 0; j < dim; ++j) A(i, j) = P(j, i) - (i == j ? one : zero);
    for (std::size_t j = 0; j < dim; ++j) {
        T s = zero;
        for (std::size_t k = 0; k < dim; ++k) s += C(j, k);
        A(dim - 1, j) = s;
    }
    std::vector<T> b(dim, zero);
    b[dim - 1] = one;
    const std::vector<T> u = solve(A, b);
    std::vector<T> pi = vecmul(u, C);
    T tot = zero;
    for (const T& v : pi) tot += v;
    for (T& v : pi) v /= tot;
    return pi;
}
 
template <class T>
std::vector<T> coefficients(const ServiceLaw<T>& svc, const T& lambda, const T& tol,
                            std::size_t nmaxCap) {
    // a truncated series warns and the solve continues, as in the reference
    return qsys_svclaw_coeffs(svc, lambda, tol, nmaxCap).cn;
}
 
/** i_plain: synchronous approximation without state expansion; state i <-> i jobs. */
template <class T>
std::vector<T> plain(const T& lambda, std::size_t c, std::size_t K, const ServiceLaw<T>& svc,
                     bool faster, const T& tol, std::size_t nmaxCap, std::vector<std::size_t>& level) {
    const std::size_t dim = K + 1;
    const T zero = num_traits<T>::from_int(0), one = num_traits<T>::from_int(1);
    Matrix<T> Q(dim, dim, zero), Qbar(dim, dim, zero), Delta(dim, dim, zero);
    Q(0, 0) = -lambda;
    Q(0, 1) = lambda;
    for (std::size_t i = 1; i <= K; ++i) {
        if (i < K) {
            Q(i, i) = -lambda;
            if (!faster && i < c) Qbar(i, i + 1) = lambda;  // joining a partial batch restarts its clock
            else Q(i, i + 1) = lambda;
        }
        Delta(i, i > c ? i - c : 0) = one;
    }
    std::vector<bool> inE(dim, false), inG(dim, true);
    inE[0] = true;
    inG[0] = false;
    const std::vector<T> cn = coefficients(svc, lambda, tol, nmaxCap);
    Matrix<T> Om, Ps;
    unif(rowmask(Q, inG), cn, lambda, Om, Ps);
    Om = rowmask(Om, inG);
    Ps = rowmask(Ps, inG);
    level.resize(dim);
    for (std::size_t i = 0; i < dim; ++i) level[i] = i;
    return emc(Q, inE, matmul(Om, Delta), Ps, Qbar);
}
 
/** A sparse transition list, summed into a dense matrix as MATLAB's sparse() does. */
template <class T>
struct Triplets {
    std::vector<std::size_t> from, to;
    std::vector<T> val;
    void add(std::size_t i, std::size_t j, const T& v) {
        from.push_back(i);
        to.push_back(j);
        val.push_back(v);
    }
    Matrix<T> dense(std::size_t dim) const {
        Matrix<T> M(dim, dim, num_traits<T>::from_int(0));
        for (std::size_t k = 0; k < val.size(); ++k) M(from[k], to[k]) += val[k];
        return M;
    }
};
 
/**
 * i_tick: advance every non-owner phase by one tick; a completed job frees its
 * server for a waiting one, which starts at phase 0 on the boundary.
 */
inline std::vector<int> tick(int a, int q, const std::vector<int>& B, int nph) {
    std::vector<int> Bn;
    int done = 0;
    for (int b : B) {
        if (b + 1 >= nph) ++done;
        else Bn.push_back(b + 1);
    }
    const int start = std::min(done, q);
    for (int k = 0; k < start; ++k) Bn.push_back(0);
    std::sort(Bn.begin(), Bn.end());
    std::vector<int> s;
    s.push_back(a);
    s.push_back(q - start);
    s.insert(s.end(), Bn.begin(), Bn.end());
    return s;
}
 
inline std::vector<int> make_state(int a, int q, const std::vector<int>& B) {
    std::vector<int> s;
    s.push_back(a);
    s.push_back(q);
    s.insert(s.end(), B.begin(), B.end());
    return s;
}
 
/** i_expanded: synchronous phase expansion; states are [a, q, sorted B phases]. */
template <class T>
std::vector<T> expanded(const T& lambda, std::size_t c, std::size_t K, const T& d, std::size_t nphU,
                        bool faster, const T& tol, std::size_t nmaxCap,
                        std::vector<std::size_t>& level) {
    const T zero = num_traits<T>::from_int(0), one = num_traits<T>::from_int(1);
    const int lo = faster ? 0 : -1;
    const int nph = static_cast<int>(nphU);
    const int ci = static_cast<int>(c), Ki = static_cast<int>(K);
    std::map<std::vector<int>, std::size_t> keys;
    std::vector<std::vector<int>> states;
    auto add = [&](const std::vector<int>& s) -> std::size_t {
        const auto it = keys.find(s);
        if (it != keys.end()) return it->second;
        states.push_back(s);
        const std::size_t id = states.size() - 1;
        keys[s] = id;
        return id;
    };
    add(make_state(0, 0, std::vector<int>()));
    // transitions: Q, Qbar, DBA (tick, owner busy), DA (owner completion), DB (tick, owner idle)
    Triplets<T> TQ, TQbar, TDBA, TDA, TDB;
    std::vector<bool> inA, inB;
    const T cT = num_traits<T>::from_int(static_cast<long>(c));
    for (std::size_t head = 0; head < states.size(); ++head) {
        const std::vector<int> s = states[head];  // copy: add() may reallocate
        const int a = s[0], q = s[1];
        const std::vector<int> B(s.begin() + 2, s.end());
        const int tot = a + q + static_cast<int>(B.size());
        inA.resize(states.size(), false);
        inB.resize(states.size(), false);
        if (tot == 0) {
            TQ.add(head, head, -lambda);
            TQ.add(head, add(make_state(1, 0, std::vector<int>())), lambda / cT);
            if (c > 1)
                TQ.add(head, add(make_state(0, 0, std::vector<int>(1, 0))),
                       lambda * (cT - one) / cT);
            inA[head] = false;
            inB[head] = false;
        } else if (a == 1) {
            inA[head] = true;
            inB[head] = false;
            if (tot < Ki) {
                std::vector<int> t;
                if (1 + static_cast<int>(B.size()) < ci) {
                    std::vector<int> Bn = B;
                    Bn.push_back(lo);
                    std::sort(Bn.begin(), Bn.end());
                    t = make_state(1, 0, Bn);
                } else {
                    t = make_state(1, q + 1, B);
                }
                TQ.add(head, head, -lambda);
                TQ.add(head, add(t), lambda);
            }
            TDBA.add(head, add(tick(1, q, B, nph)), one);
            TDA.add(head, add(q > 0 ? make_state(1, q - 1, B) : make_state(0, 0, B)), one);
        } else {
            inA[head] = false;
            inB[head] = true;
            if (tot < Ki) {
                std::vector<int> Bn;
                if (faster) {
                    for (int b : B)
                        if (b + 1 < nph) Bn.push_back(b + 1);  // a credited phase may complete the service
                } else {
                    Bn = B;
                }
                TQ.add(head, head, -lambda);
                TQbar.add(head, add(make_state(1, 0, Bn)), lambda);
            }
            TDB.add(head, add(tick(0, q, B, nph)), one);
        }
    }
    const std::size_t dim = states.size();
    inA.resize(dim, false);
    inB.resize(dim, false);
    std::vector<bool> inE(dim, false);
    for (std::size_t i = 0; i < dim; ++i) inE[i] = !inA[i] && !inB[i];
    const Matrix<T> Q = TQ.dense(dim), Qbar = TQbar.dense(dim);
    const Matrix<T> DBA = TDBA.dense(dim), DA = TDA.dense(dim), DB = TDB.dense(dim);
 
    const T theta = d / num_traits<T>::from_int(static_cast<long>(nphU));
    const std::vector<T> cn =
        coefficients(ServiceLaw<T>::deterministic(theta), lambda, tol, nmaxCap);
    Matrix<T> eA, iA;
    unif(rowmask(Q, inA), cn, lambda, eA, iA);
    const Matrix<T> OmTick = matmul(rowmask(eA, inA), DBA);  // one tick with the owner busy
    const Matrix<T> PsiTick = colmask(rowmask(iA, inA), inA);
    Matrix<T> PsiA = PsiTick;
    Matrix<T> W = OmTick;
    for (std::size_t j = 1; j < nphU; ++j) {
        PsiA = matadd(PsiA, matmul(W, PsiTick));
        W = matmul(W, OmTick);
    }
    const Matrix<T> OmA = matmul(W, DA);  // owner service: N ticks then its completion
    Matrix<T> eB, iB;
    unif(rowmask(Q, inB), cn, lambda, eB, iB);
    const Matrix<T> OmB = matmul(rowmask(eB, inB), DB);
    const Matrix<T> PsiB = colmask(rowmask(iB, inB), inB);
    level.resize(dim);
    for (std::size_t i = 0; i < dim; ++i)
        level[i] = static_cast<std::size_t>(states[i][0] + states[i][1]) + (states[i].size() - 2);
    (void)zero;
    return emc(Q, inE, matadd(OmA, OmB), matadd(PsiA, PsiB), Qbar);
}
 
}  // namespace mgck_sync_detail
 
/**
 * Synchronous approximation of the M/G/c/K queue with dual bounds.
 *
 * @param lambda  Poisson arrival rate (> 0)
 * @param c       number of identical FCFS servers (>= 1)
 * @param K       system capacity, jobs in service included (K >= c)
 * @param svc     service law, as for qsys_mapg1k
 * @param variant "faster" (lower bounds on the queue) or "slower" (upper bounds)
 * @param nphases 0 for the approximation without expansion; N >= 1 for the
 *                synchronous phase expansion, deterministic service only
 * @param tol     uniformization truncation tolerance
 * @param nmaxCap cap on the uniformization order
 */
template <class T>
MgckSyncResult<T> qsys_mgck_sync(const T& lambda, std::size_t c, std::size_t K,
                                 const ServiceLaw<T>& svc, const std::string& variant,
                                 std::size_t nphases, const T& tol, std::size_t nmaxCap) {
    static_assert(num_traits<T>::has_transcendental,
                  "qsys_mgck_sync requires transcendental arithmetic");
    const T zero = num_traits<T>::from_int(0), one = num_traits<T>::from_int(1);
    if (!(lambda > zero)) throw InputError("qsys_mgck_sync: lambda must be a positive real scalar");
    if (c < 1) throw InputError("qsys_mgck_sync: c must be a positive integer");
    if (K < c) throw InputError("qsys_mgck_sync: K must be an integer >= c");
    std::string v = variant;
    std::transform(v.begin(), v.end(), v.begin(),
                   [](unsigned char ch) { return static_cast<char>(std::tolower(ch)); });
    if (v != "faster" && v != "slower")
        throw InputError("qsys_mgck_sync: Unknown variant '" + variant +
                         "'; use 'faster' or 'slower'.");
    const bool faster = (v == "faster");
 
    std::vector<std::size_t> level;
    std::vector<T> pi;
    if (nphases == 0) {
        pi = mgck_sync_detail::plain(lambda, c, K, svc, faster, tol, nmaxCap, level);
    } else {
        if (svc.kind != ServiceKind::Deterministic)
            throw UnsupportedError(ErrorKind{"Distribution"},
                "qsys_mgck_sync: The synchronous phase expansion (nphases >= 1) requires a "
                "deterministic service time.");
        pi = mgck_sync_detail::expanded(lambda, c, K, svc.det, nphases, faster, tol, nmaxCap,
                                        level);
    }
 
    MgckSyncResult<T> r;
    r.queueLengthDist.assign(K + 1, zero);
    for (std::size_t s = 0; s < pi.size(); ++s) r.queueLengthDist[level[s]] += pi[s];
    T L = zero, Lq = zero, busy = zero;
    for (std::size_t n = 0; n <= K; ++n) {
        const T pn = r.queueLengthDist[n];
        L += num_traits<T>::from_int(static_cast<long>(n)) * pn;
        if (n > c) Lq += num_traits<T>::from_int(static_cast<long>(n - c)) * pn;
        busy += num_traits<T>::from_int(static_cast<long>(n < c ? n : c)) * pn;
    }
    r.p0 = r.queueLengthDist[0];
    r.pK = r.queueLengthDist[K];
    r.meanQueueLength = L;
    r.meanQueueLengthQ = Lq;
    r.throughput = lambda * (one - r.pK);
    r.lossProbability = r.pK;
    r.meanSojournTime = L / r.throughput;
    r.meanWaitingTime = Lq / r.throughput;
    r.utilization = busy / num_traits<T>::from_int(static_cast<long>(c));
    r.variant = v;
    r.nphases = nphases;
    r.stateSpaceSize = pi.size();
    if (nphases == 0)
        r.analyzer = "sync." + v + ":M/G/" + std::to_string(c) + "/" + std::to_string(K);
    else
        r.analyzer = "sync." + v + ".n" + std::to_string(nphases) + ":M/D/" + std::to_string(c) +
                     "/" + std::to_string(K);
    return r;
}
 
/** qsys_mgck_sync with the reference defaults tol = 1e-12, nmax = 200000. */
template <class T>
MgckSyncResult<T> qsys_mgck_sync(const T& lambda, std::size_t c, std::size_t K,
                                 const ServiceLaw<T>& svc, const std::string& variant,
                                 std::size_t nphases = 0) {
    return qsys_mgck_sync(lambda, c, K, svc, variant, nphases,
                          T(num_traits<T>::from_double(1e-12)), static_cast<std::size_t>(200000));
}
 
}  // namespace qsys
}  // namespace line
 
#endif  // LINE_API_QSYS_QSYS_MGCK_SYNC_H