qsys_mg1_busy.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_MG1_BUSY_H
#define LINE_API_QSYS_QSYS_MG1_BUSY_H
 
/**
 * @file
 * @ingroup api_qsys
 * Busy-period measures of the M/G/1 queue.
 *
 * Port of matlab/src/api/qsys/qsys_mg1_busy.m. The busy period BP has
 *
 *   E[BP] = E[S]/(1-rho),   E[BP^2] = E[S^2]/(1-rho)^3,
 *   E[BP^3] = E[S^3]/(1-rho)^4 + 3 lambda E[S^2]^2/(1-rho)^5,
 *
 * the number served in it has mean 1/(1-rho) and variance
 * (rho(1-rho) + lambda^2 E[S^2])/(1-rho)^3, and given the service LST its own
 * LST solves Takacs' equation b(s) = B(s + lambda - lambda b(s)), here by the
 * reference's fixed-point iteration from b = 0.
 *
 * ARITHMETIC: double, as the reference.
 */
 
#include <cmath>
#include <functional>
#include <vector>
 
#include "line/api/qsys/qsys_res.h"
#include "line/util/error.h"
 
namespace line {
namespace qsys {
 
/** MATLAB's name-value options of qsys_mg1_busy; a negative ES3 means not given. */
struct Mg1BusyOptions {
    double ES3 = -1.0;                    ///< third moment of the service time
    std::function<double(double)> lst;    ///< the service LST at real arguments
    std::vector<double> sPoints;          ///< where to evaluate the busy-period LST
};
 
/**
 * @brief Busy-period measures of the M/G/1 queue.
 *
 * @param lambda arrival rate
 * @param mu     service rate
 * @param cs     coefficient of variation of the service time
 * @param o      ES3, lst and sPoints options
 */
inline QsysRes qsys_mg1_busy(double lambda, double mu, double cs, const Mg1BusyOptions& o = Mg1BusyOptions()) {
    if (!(lambda > 0.0) || !(mu > 0.0) || !(cs >= 0.0))
        throw InputError(ErrorKind("Argument"),
                         "qsys_mg1_busy: lambda and mu must be positive scalars and cs non-negative");
    const double rho = lambda / mu;
    if (rho >= 1.0) resdetail::throw_unstable("qsys_mg1_busy", "rho", rho);
    const double ES = 1.0 / mu;
    const double ES2 = (1.0 + cs * cs) / (mu * mu);
    const double EBP = ES / (1.0 - rho);
    const double EBP2 = ES2 / std::pow(1.0 - rho, 3);
    QsysRes r;
    r.set("meanBusyPeriod", EBP);
    r.set("varBusyPeriod", EBP2 - EBP * EBP);
    if (o.ES3 >= 0.0)
        r.set("thirdMomentBusyPeriod", o.ES3 / std::pow(1.0 - rho, 4) + 3.0 * lambda * ES2 * ES2 / std::pow(1.0 - rho, 5));
    r.set("meanServed", 1.0 / (1.0 - rho));
    r.set("varServed", (rho * (1.0 - rho) + lambda * lambda * ES2) / std::pow(1.0 - rho, 3));
    r.set("meanIdlePeriod", 1.0 / lambda);
    r.set("busyFraction", rho);
    r.set("busyPeriodRate", lambda * (1.0 - rho));
    if (!o.sPoints.empty()) {
        if (!o.lst)
            throw InputError(ErrorKind("Argument"), "qsys_mg1_busy: sPoints needs the service LST (option 'lst')");
        for (double s : o.sPoints)
            if (s < 0.0) throw InputError(ErrorKind("Argument"), "qsys_mg1_busy: sPoints must be non-negative");
        std::vector<double> v(o.sPoints.size());
        for (std::size_t k = 0; k < o.sPoints.size(); ++k) {
            const double s = o.sPoints[k];
            double b = 0.0;
            for (int it = 1; it <= 100000; ++it) {
                const double bn = o.lst(s + lambda - lambda * b);
                if (std::fabs(bn - b) <= 1e-14) {
                    b = bn;
                    break;
                }
                b = bn;
            }
            v[k] = b;
        }
        r.set("sPoints", o.sPoints);
        r.set("busyPeriodLst", v);
    }
    r.analyzer = "exact:M/G/1-busy";
    return r;
}
 
}  // namespace qsys
}  // namespace line
 
#endif  // LINE_API_QSYS_QSYS_MG1_BUSY_H