qsys_mgc_ps.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_MGC_PS_H
#define LINE_API_QSYS_QSYS_MGC_PS_H
 
/**
 * @file
 * @ingroup api_qsys
 * Exact analysis of the M/G/c processor-sharing queue.
 *
 * Port of matlab/src/api/qsys/qsys_mgc_ps.m. The queue-length law is
 * insensitive and equals that of M/M/c: with A = lambda/mu, a = A/c and the
 * Erlang C probability C,
 *
 *   E[N] = C a/(1 - a) + A,   Var[N] = a C (1 + a - a C)/(1 - a)^2 + A (1 + C),
 *
 * and E[T] = E[N]/lambda by Little's law. The distribution is the M/M/c law,
 * normalized in log space and truncated where its geometric tail falls below
 * 1e-14.
 *
 * ARITHMETIC: double, as the reference.
 */
 
#include <algorithm>
#include <cmath>
#include <string>
#include <vector>
 
#include "line/api/qsys/qsys_res.h"
#include "line/util/error.h"
 
namespace line {
namespace qsys {
 
/**
 * @brief Exact M/G/c-PS queue.
 *
 * @param lambda arrival rate
 * @param mu     service rate of one server
 * @param c      number of servers
 */
inline QsysRes qsys_mgc_ps(double lambda, double mu, unsigned c) {
    using namespace resdetail;
    if (!(lambda > 0.0)) throw Error("qsys_mgc_ps: lambda must be a positive scalar");
    if (!(mu > 0.0)) throw Error("qsys_mgc_ps: mu must be a positive scalar");
    if (c < 1) throw Error("qsys_mgc_ps: c must be a positive integer");
    const double A = lambda / mu;
    const double a = A / c;
    if (a >= 1.0) throw_unstable("qsys_mgc_ps", "lambda/(c mu)", a);
    double B = 1.0;
    for (unsigned k = 1; k <= c; ++k) B = A * B / (k + A * B);
    const double C = c * B / (c - A * (1.0 - B));
    const double EQ = C * a / (1.0 - a);
    const double VQ = a * C * (1.0 + a - a * C) / ((1.0 - a) * (1.0 - a));
    const double EN = EQ + A;
    const double VN = VQ + A * (1.0 + C);
    std::vector<double> logp(c + 1);
    for (unsigned n = 0; n <= c; ++n) logp[n] = n * std::log(A) - gammaln(n + 1.0);
    const double logpc = logp[c];
    const double ntail = std::max(0.0, std::ceil(std::log(1e-14 * (1.0 - a)) / std::log(a)));
    for (long j = 1; j <= static_cast<long>(ntail); ++j) logp.push_back(logpc + j * std::log(a));
    std::vector<double> lognorm(logp.begin(), logp.begin() + c);
    lognorm.push_back(logpc - std::log(1.0 - a));
    double mx = lognorm[0];
    for (double v : lognorm) mx = std::max(mx, v);
    double s = 0.0;
    for (double v : lognorm) s += std::exp(v - mx);
    const double logP0 = -(mx + std::log(s));
    std::vector<double> pn(logp.size());
    for (std::size_t i = 0; i < logp.size(); ++i) pn[i] = std::exp(logp[i] + logP0);
    QsysRes r;
    r.set("meanQueueLength", EN);
    r.set("varQueueLength", VN);
    r.set("meanSojournTime", EN / lambda);
    r.set("throughput", lambda);
    r.set("utilization", a);
    r.set("queueLengthDist", pn);
    r.set("tailRatio", a);
    r.analyzer = "exact:M/G/" + std::to_string(c) + "-PS";
    return r;
}
 
}  // namespace qsys
}  // namespace line
 
#endif  // LINE_API_QSYS_QSYS_MGC_PS_H