qsys_mmcc.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_MMCC_H
#define LINE_API_QSYS_QSYS_MMCC_H
 
/**
 * @file
 * @ingroup api_qsys
 * Exact analysis of the M/M/c/c (Erlang loss) system.
 *
 * Port of matlab/src/api/qsys/qsys_mmcc.m. The blocking probability is the
 * Erlang B formula by its stable recursion B_k = A B_{k-1} / (k + A B_{k-1});
 * the occupancy is the truncated Poisson law, normalized in log space. The
 * result is insensitive to the service law beyond its mean.
 *
 * ARITHMETIC: double (gammaln, log-space normalization), as the reference.
 */
 
#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/M/c/c Erlang loss system.
 *
 * @param lambda arrival rate
 * @param mu     service rate of one server
 * @param c      number of servers
 */
inline QsysRes qsys_mmcc(double lambda, double mu, unsigned c) {
    using namespace resdetail;
    if (!(lambda > 0.0)) throw Error("qsys_mmcc: lambda must be a positive scalar");
    if (!(mu > 0.0)) throw Error("qsys_mmcc: mu must be a positive scalar");
    if (c < 1) throw Error("qsys_mmcc: c must be a positive integer");
    const double A = lambda / mu;
    double B = 1.0;
    for (unsigned k = 1; k <= c; ++k) B = A * B / (k + A * B);
    std::vector<double> logp(c + 1);
    for (unsigned n = 0; n <= c; ++n) logp[n] = n * std::log(A) - gammaln(n + 1.0);
    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));
    double EN2 = 0.0;
    for (unsigned n = 0; n <= c; ++n) {
        p[n] /= s;
        EN2 += static_cast<double>(n) * n * p[n];
    }
    const double EN = A * (1.0 - B);
    QsysRes r;
    r.set("blockingProb", B);
    r.set("throughput", lambda * (1.0 - B));
    r.set("lostRate", lambda * B);
    r.set("offeredTraffic", A);
    r.set("carriedTraffic", EN);
    r.set("overflowTraffic", A * B);
    r.set("meanQueueLength", EN);
    r.set("varQueueLength", EN2 - EN * EN);
    r.set("utilization", EN / c);
    r.set("meanSojournTime", 1.0 / mu);
    r.set("queueLengthDist", p);
    r.analyzer = "exact:M/M/" + std::to_string(c) + "/" + std::to_string(c);
    return r;
}
 
}  // namespace qsys
}  // namespace line
 
#endif  // LINE_API_QSYS_QSYS_MMCC_H