qsys_mg1kk.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_MG1KK_H
#define LINE_API_QSYS_QSYS_MG1KK_H
 
/**
 * @file
 * @ingroup api_qsys
 * Exact analysis of the finite-source M/G/1/K/K (machine-repair) queue.
 *
 * Port of matlab/src/api/qsys/qsys_mg1kk.m. K sources each operate for an
 * Exp(alpha) time and then queue FCFS at one server of general service law B.
 * Takacs (1962) gives the idle probability from the service LST at the
 * points k alpha,
 *
 *   P0 = 1 / (1 + K alpha E[S] sum_{n=0}^{K-1} C(K-1, n) prod_{k=1}^n (1 - B(k alpha))/B(k alpha)),
 *
 * whence X = (1 - P0)/E[S] and E[T] = K/X - 1/alpha. The queue-length law
 * comes from the chain embedded at departures, whose transition probabilities
 * are alternating sums of B(k alpha); MATLAB's check that its P0 agrees with
 * Takacs' is kept: when the alternating sums cancel it warns and returns the
 * distribution (and the variance) as NaN, the means being the Takacs sum.
 *
 * ARITHMETIC: double, as the reference. MATLAB's overdetermined `A \ b`
 * (least squares on the consistent system [P' - I; 1] pi = [0; 1]) is solved
 * here by LU on the square system that replaces one of the K dependent
 * balance rows by the normalization; the solution of a consistent system is
 * unique, so the two agree to rounding.
 */
 
#include <algorithm>
#include <cmath>
#include <cstdio>
#include <functional>
#include <iostream>
#include <limits>
#include <string>
#include <vector>
 
#include "line/api/qsys/qsys_res.h"
#include "line/lang/distribution.h"
#include "line/util/error.h"
#include "line/util/line_console.h"
#include "line/util/lu.h"
#include "line/util/matrix.h"
 
namespace line {
namespace qsys {
 
/**
 * @brief Exact M/G/1/K/K finite-source queue, service given by its LST.
 *
 * @param alpha request rate of one operating source
 * @param K     number of sources
 * @param lst   service-time LST at real arguments
 * @param ES    mean service time
 */
inline QsysRes qsys_mg1kk(double alpha, unsigned K, const std::function<double(double)>& lst, double ES) {
    using resdetail::nchoosek_r;
    if (!(alpha > 0.0)) throw Error("qsys_mg1kk: alpha must be a positive real scalar");
    if (K < 1) throw Error("qsys_mg1kk: K must be a positive integer");
    if (!lst) throw Error("qsys_mg1kk: svc must be a Distribution or an LST function handle");
    if (!(ES > 0.0)) throw Error("qsys_mg1kk: with an LST handle the mean service time ES must be given");
 
    std::vector<double> Bs(K + 1);
    for (unsigned k = 0; k <= K; ++k) Bs[k] = lst(k * alpha);
    double Bn = 1.0, S = 1.0;
    for (unsigned nn = 1; nn + 1 <= K; ++nn) {
        Bn = Bn * (1.0 - Bs[nn]) / Bs[nn];
        S += nchoosek_r(K - 1.0, nn) * Bn;
    }
    const double P0 = 1.0 / (1.0 + K * alpha * ES * S);
    const double X = (1.0 - P0) / ES;
    const double ET = K / X - 1.0 / alpha;
    const double EW = ET - ES;
 
    Matrix<double> P(K, K, 0.0);
    for (unsigned j = 0; j < K; ++j) {
        const unsigned nop = j == 0 ? K - 1 : K - j;
        const unsigned base = j == 0 ? 0 : j - 1;
        for (unsigned i = 0; i <= nop; ++i) {
            double acc = 0.0;
            for (unsigned l = 0; l <= i; ++l)
                acc += ((l % 2) ? -1.0 : 1.0) * nchoosek_r(i, l) * Bs[nop - i + l];
            P(j, base + i) += nchoosek_r(nop, i) * acc;
        }
    }
    for (unsigned j = 0; j < K; ++j) {
        double rs = 0.0;
        for (unsigned i = 0; i < K; ++i) {
            P(j, i) = std::max(P(j, i), 0.0);
            rs += P(j, i);
        }
        for (unsigned i = 0; i < K; ++i) P(j, i) /= rs;
    }
    // [P' - I] with its last balance row replaced by the normalization
    Matrix<double> A(K, K, 0.0);
    std::vector<double> b(K, 0.0);
    for (unsigned r = 0; r + 1 < K; ++r)
        for (unsigned cc = 0; cc < K; ++cc) A(r, cc) = P(cc, r) - (r == cc ? 1.0 : 0.0);
    for (unsigned cc = 0; cc < K; ++cc) A(K - 1, cc) = 1.0;
    b[K - 1] = 1.0;
    std::vector<double> Pi = solve(A, b);
 
    std::vector<double> pn(K + 1, 0.0);
    double mass = 0.0;
    for (unsigned n = 0; n < K; ++n) {
        pn[n] = X * Pi[n] / ((static_cast<double>(K) - n) * alpha);
        mass += pn[n];
    }
    pn[K] = 1.0 - mass;
    if (std::fabs(pn[0] - P0) > 1e-6 * std::max(P0, 1e-12) || pn[K] < -1e-9) {
        char buf[384];
        std::snprintf(buf, sizeof buf,
                      "qsys_mg1kk: the embedded-chain distribution disagrees with the Takacs P_0 (%g vs %g): "
                      "the alternating arrival sums lost accuracy at K=%u; queueLengthDist, arrivalDist "
                      "and varQueueLength are returned as NaN, the means are unaffected.",
                      pn[0], P0, K);
        // the means above are the Takacs sum of positive terms and do not
        // read the chain, so only the distribution is withdrawn
        if (util::LineConsole::get_verbose() != util::VerboseLevel::SILENT)
            std::cerr << "[LINE] Warning: " << buf << std::endl;
        std::fill(pn.begin(), pn.end(), std::numeric_limits<double>::quiet_NaN());
        std::fill(Pi.begin(), Pi.end(), std::numeric_limits<double>::quiet_NaN());
    } else {
        pn[K] = std::max(pn[K], 0.0);
    }
    double EN = 0.0, EN2 = 0.0;
    for (unsigned n = 0; n <= K; ++n) {
        EN += n * pn[n];
        EN2 += static_cast<double>(n) * n * pn[n];
    }
    QsysRes r;
    r.set("meanQueueLength", X * ET);
    r.set("varQueueLength", EN2 - EN * EN);
    r.set("meanQueueLengthQ", X * EW);
    r.set("meanWaitingTime", EW);
    r.set("meanSojournTime", ET);
    r.set("throughput", X);
    r.set("utilization", 1.0 - P0);
    r.set("meanOperating", K - X * ET);
    r.set("queueLengthDist", pn);
    r.set("arrivalDist", Pi);
    r.analyzer = "exact:M/G/1/" + std::to_string(K) + "/" + std::to_string(K);
    return r;
}
 
/**
 * @brief Exact M/G/1/K/K finite-source queue, service given as a Distribution
 *        (its LST and mean, as MATLAB's svc.evalLST / svc.getMean).
 */
inline QsysRes qsys_mg1kk(double alpha, unsigned K, const lang::Distrib<double>& svc) {
    const lang::Distrib<double> d = svc;
    return qsys_mg1kk(alpha, K, [d](double s) { return lang::dist_lst(d, s); }, lang::dist_moment(d, 1));
}
 
}  // namespace qsys
}  // namespace line
 
#endif  // LINE_API_QSYS_QSYS_MG1KK_H