qsys_mmckk.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_MMCKK_H
#define LINE_API_QSYS_QSYS_MMCKK_H
 
/**
 * @file
 * @ingroup api_qsys
 * Exact analysis of the finite-source M/M/c/K/K (machine-repair) queue.
 *
 * Port of matlab/src/api/qsys/qsys_mmckk.m. K sources each operate for an
 * Exp(alpha) time and then queue for one of c Exp(mu) servers. The number at
 * the queue is a birth-death chain with
 *
 *   p_n ~ C(K, n) (alpha/mu)^n               n <= c
 *   p_n ~ C(K, n) (alpha/mu)^n n!/(c! c^(n-c))   n > c,
 *
 * normalized in log space. An arriving request sees Pi_n ~ (K - n) p_n
 * (n = 0..K-1); one that finds m-1 waiting waits an Erlang(m, c mu) time, which
 * gives the waiting and sojourn distributions in closed form.
 *
 * ARITHMETIC: double (gammaln, incomplete gamma), as the reference.
 */
 
#include <algorithm>
#include <cmath>
#include <functional>
#include <string>
#include <vector>
 
#include "line/api/qsys/qsys_res.h"
#include "line/util/error.h"
 
namespace line {
namespace qsys {
 
/**
 * @brief Exact M/M/c/K/K finite-source queue.
 *
 * @param alpha request rate of one operating source
 * @param mu    service rate of one server
 * @param c     number of servers
 * @param K     number of sources, K >= c
 * @param o     tPoints and percentiles options (nmax is not used)
 */
inline QsysRes qsys_mmckk(double alpha, double mu, unsigned c, unsigned K,
                          const QsysResOptions& o = QsysResOptions()) {
    using namespace resdetail;
    if (!(alpha > 0.0)) throw Error("qsys_mmckk: alpha must be a positive real scalar");
    if (!(mu > 0.0)) throw Error("qsys_mmckk: mu must be a positive real scalar");
    if (c < 1) throw Error("qsys_mmckk: c must be a positive integer");
    if (K < c) throw Error("qsys_mmckk: K must be an integer >= c");
    check_percentiles(o.percentiles, "qsys_mmckk");
 
    std::vector<double> logp(K + 1);
    for (unsigned n = 0; n <= K; ++n) {
        logp[n] = gammaln(K + 1.0) - gammaln(n + 1.0) - gammaln(static_cast<double>(K) - n + 1.0) +
                  n * std::log(alpha / mu);
        if (n > c) logp[n] += gammaln(n + 1.0) - gammaln(c + 1.0) - (static_cast<double>(n) - c) * std::log(c);
    }
    double mx = logp[0];
    for (double v : logp) mx = std::max(mx, v);
    std::vector<double> p(K + 1);
    double s = 0.0;
    for (unsigned n = 0; n <= K; ++n) s += (p[n] = std::exp(logp[n] - mx));
    for (double& v : p) v /= s;
 
    double L = 0.0, L2 = 0.0, Lq = 0.0;
    for (unsigned n = 0; n <= K; ++n) {
        L += n * p[n];
        L2 += static_cast<double>(n) * n * p[n];
        if (n > c) Lq += (static_cast<double>(n) - c) * p[n];
    }
    const double X = alpha * (K - L);
    std::vector<double> Pi(K);
    for (unsigned n = 0; n < K; ++n) Pi[n] = (static_cast<double>(K) - n) * p[n] / (K - L);
 
    const double theta = c * mu;
    std::vector<double> Piw, m;
    for (unsigned nw = c; nw + 1 <= K; ++nw) {
        Piw.push_back(Pi[nw]);
        m.push_back(static_cast<double>(nw - c + 1));
    }
    double probWait = 0.0, EW = 0.0, EW2 = 0.0;
    for (std::size_t i = 0; i < m.size(); ++i) {
        probWait += Piw[i];
        EW += Piw[i] * m[i];
        EW2 += Piw[i] * m[i] * (m[i] + 1.0);
    }
    EW /= theta;
    EW2 /= theta * theta;
 
    QsysRes r;
    r.set("meanQueueLength", L);
    r.set("varQueueLength", L2 - L * L);
    r.set("meanQueueLengthQ", Lq);
    r.set("meanWaitingTime", EW);
    r.set("varWaitingTime", EW2 - EW * EW);
    r.set("meanSojournTime", EW + 1.0 / mu);
    r.set("varSojournTime", EW2 - EW * EW + 1.0 / (mu * mu));
    r.set("probWait", probWait);
    r.set("condMeanWaitingTime", probWait > 0.0 ? EW / probWait : 0.0);
    r.set("throughput", X);
    r.set("utilization", X / (c * mu));
    r.set("meanOperating", K - L);
    r.set("queueLengthDist", p);
    r.set("arrivalDist", Pi);
 
    const std::function<double(double)> Wccdf = [Piw, m, theta](double t) {
        if (t < 0.0) return 1.0;
        double v = 0.0;
        for (std::size_t i = 0; i < m.size(); ++i) v += Piw[i] * gammainc_upper(theta * t, m[i]);
        return v;
    };
    const std::function<double(double)> Tccdf = [Pi, Piw, m, theta, mu, c](double t) {
        if (t < 0.0) return 1.0;
        double c0 = 0.0;
        for (std::size_t i = 0; i < std::min<std::size_t>(c, Pi.size()); ++i) c0 += Pi[i];
        double v = c0 * std::exp(-mu * t);
        if (c == 1) {
            for (std::size_t i = 0; i < m.size(); ++i) v += Piw[i] * gammainc_upper(mu * t, m[i] + 1.0);
        } else {
            const double d = theta - mu;
            for (std::size_t i = 0; i < m.size(); ++i) {
                const double lower = gammainc_lower(d * t, m[i]);
                const double cross =
                    lower > 0.0 ? std::exp(m[i] * std::log(theta / d) - mu * t + std::log(lower)) : 0.0;
                v += Piw[i] * (gammainc_upper(theta * t, m[i]) + cross);
            }
        }
        return v;
    };
    if (!o.tPoints.empty()) {
        r.set("tPoints", o.tPoints);
        r.set("waitingTimeCdf", cdf_of(Wccdf, o.tPoints));
        r.set("sojournTimeCdf", cdf_of(Tccdf, o.tPoints));
    }
    if (!o.percentiles.empty()) {
        const double hiT = (K + 50.0) / std::min(mu, theta);
        r.set("percentiles", o.percentiles);
        r.set("waitingTimePercentile", pct_of(Wccdf, o.percentiles, hiT));
        r.set("sojournTimePercentile", pct_of(Tccdf, o.percentiles, hiT));
    }
    r.analyzer = "exact:M/M/" + std::to_string(c) + "/" + std::to_string(K) + "/" + std::to_string(K);
    return r;
}
 
}  // namespace qsys
}  // namespace line
 
#endif  // LINE_API_QSYS_QSYS_MMCKK_H