qsys_mg1_conservation.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_CONSERVATION_H
#define LINE_API_QSYS_QSYS_MG1_CONSERVATION_H
 
/**
 * @file
 * @ingroup api_qsys
 * Kleinrock's conservation law for work-conserving non-preemptive M/G/1
 * disciplines.
 *
 * Port of matlab/src/api/qsys/qsys_mg1_conservation.m. For any such
 * discipline sum_i rho_i Wq_i is invariant and equals
 *
 *   V = rho sum_i lambda_i E[S_i^2] / (2 (1 - rho)),
 *
 * so V/rho is the FCFS waiting time; given per-class waiting times W the
 * residual sum_i rho_i W_i - V checks a discipline against the law.
 *
 * ARITHMETIC: field; exact for T = Rational. A residual not asked for is NaN
 * in MATLAB, which the exact field cannot hold, so `hasResidual` flags it and
 * the res field is set only when W was given.
 */
 
#include <cstddef>
#include <limits>
#include <vector>
 
#include "line/api/qsys/qsys_res.h"
#include "line/num/number.h"
#include "line/util/error.h"
 
namespace line {
namespace qsys {
 
/** MATLAB's [V, res] outputs. */
template <class T>
struct Mg1ConservationResult {
    T V;                      ///< the invariant
    T fcfsWaitingTime;        ///< V/rho
    T residual;               ///< sum_i rho_i W_i - V, when W was given
    bool hasResidual = false; ///< false: MATLAB's residual = NaN
};
 
/**
 * @brief The M/G/1 conservation-law invariant.
 *
 * @param lambda per-class arrival rates
 * @param mu     per-class service rates
 * @param cs     per-class coefficients of variation of the service time
 * @param W      per-class mean waiting times to check; empty when not given
 */
template <class T>
Mg1ConservationResult<T> qsys_mg1_conservation(const std::vector<T>& lambda, const std::vector<T>& mu,
                                               const std::vector<T>& cs,
                                               const std::vector<T>& W = std::vector<T>()) {
    const T zero = num_traits<T>::from_int(0), one = num_traits<T>::from_int(1);
    const T two = num_traits<T>::from_int(2);
    const std::size_t R = lambda.size();
    bool bad = mu.size() != R || cs.size() != R;
    for (std::size_t i = 0; !bad && i < R; ++i)
        if (lambda[i] <= zero || mu[i] <= zero || cs[i] < zero) bad = true;
    if (bad)
        throw InputError(ErrorKind("Argument"),
                         "qsys_mg1_conservation: lambda, mu, cs must be equal-length, lambda and mu positive, "
                         "cs non-negative");
    std::vector<T> rho_i(R);
    T rho = zero, lES2 = zero;
    for (std::size_t i = 0; i < R; ++i) {
        rho_i[i] = T(lambda[i] / mu[i]);
        rho += rho_i[i];
        lES2 += T(lambda[i] * (one + cs[i] * cs[i]) / (mu[i] * mu[i]));
    }
    if (rho >= one) resdetail::throw_unstable("qsys_mg1_conservation", "rho", num_traits<T>::to_double(rho));
    Mg1ConservationResult<T> r;
    r.V = T(rho * lES2 / (two * (one - rho)));
    r.fcfsWaitingTime = T(r.V / rho);
    r.residual = zero;
    if (!W.empty()) {
        if (W.size() != R)
            throw InputError(ErrorKind("Argument"), "qsys_mg1_conservation: W must have one entry per class");
        T s = zero;
        for (std::size_t i = 0; i < R; ++i) s += T(rho_i[i] * W[i]);
        r.residual = T(s - r.V);
        r.hasResidual = true;
    }
    return r;
}
 
}  // namespace qsys
}  // namespace line
 
#endif  // LINE_API_QSYS_QSYS_MG1_CONSERVATION_H