qsys_engset.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_ENGSET_H
#define LINE_API_QSYS_QSYS_ENGSET_H
 
/**
 * @file
 * @ingroup api_qsys
 * Exact analysis of the Engset loss system M/G/c/c/K: K sources, each idle for
 * an Exp(alpha) time, c servers and no waiting room.
 *
 * Port of matlab/src/api/qsys/qsys_engset.m. The number of busy servers has
 * the truncated binomial law p_n ~ C(K, n) beta^n, n = 0..c, beta = alpha/mu,
 * which is insensitive to the service law. Time congestion is p_c; call
 * congestion (what an arriving request sees) is the time congestion of the
 * same system with K-1 sources.
 *
 * ARITHMETIC: double (gammaln, log-space normalization), 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 {
 
namespace engsetdetail {
/** The truncated binomial occupancy law, normalized in log space. */
inline std::vector<double> engset_dist(double beta, unsigned c, unsigned K) {
    using resdetail::gammaln;
    std::vector<double> logp(c + 1);
    for (unsigned n = 0; n <= c; ++n)
        logp[n] = gammaln(K + 1.0) - gammaln(n + 1.0) - gammaln(static_cast<double>(K) - n + 1.0) + n * std::log(beta);
    double mx = logp[0];
    for (double v : logp) mx = std::max(mx, v);
    std::vector<double> p(c + 1);
    double s = 0.0;
    for (unsigned n = 0; n <= c; ++n) s += (p[n] = std::exp(logp[n] - mx));
    for (double& v : p) v /= s;
    return p;
}
}  // namespace engsetdetail
 
/**
 * @brief Exact Engset M/G/c/c/K loss system.
 *
 * @param alpha request rate of one idle source
 * @param mu    service rate of one server
 * @param c     number of servers
 * @param K     number of sources, K >= c
 */
inline QsysRes qsys_engset(double alpha, double mu, unsigned c, unsigned K) {
    if (!(alpha > 0.0)) throw Error("qsys_engset: alpha must be a positive real scalar");
    if (!(mu > 0.0)) throw Error("qsys_engset: mu must be a positive real scalar");
    if (c < 1) throw Error("qsys_engset: c must be a positive integer");
    if (K < c) throw Error("qsys_engset: K must be an integer >= c");
    const double beta = alpha / mu;
    const std::vector<double> p = engsetdetail::engset_dist(beta, c, K);
    double callCong = 0.0;
    if (K > 1) {
        const std::vector<double> pK1 = engsetdetail::engset_dist(beta, std::min(c, K - 1), K - 1);
        if (c <= K - 1) callCong = pK1.back();
    }
    double EN = 0.0, EN2 = 0.0;
    for (unsigned n = 0; n <= c; ++n) {
        EN += n * p[n];
        EN2 += static_cast<double>(n) * n * p[n];
    }
    QsysRes r;
    r.set("timeCongestion", p.back());
    r.set("callCongestion", callCong);
    r.set("offeredRate", alpha * (K - EN));
    r.set("throughput", mu * EN);
    r.set("meanBusyServers", EN);
    r.set("varBusyServers", EN2 - EN * EN);
    r.set("utilization", EN / c);
    r.set("meanSojournTime", 1.0 / mu);
    r.set("queueLengthDist", p);
    r.analyzer = "exact:Engset M/G/" + std::to_string(c) + "/" + std::to_string(c) + "/" + std::to_string(K);
    return r;
}
 
}  // namespace qsys
}  // namespace line
 
#endif  // LINE_API_QSYS_QSYS_ENGSET_H