qsys_mg1_setup.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_SETUP_H
#define LINE_API_QSYS_QSYS_MG1_SETUP_H
 
/**
 * @file
 * @ingroup api_qsys
 * M/G/1 with setup times: the server switches off when empty and needs a setup
 * time U before serving the job that finds it off.
 *
 * Port of matlab/src/api/qsys/qsys_mg1_setup.m. By the decomposition of
 * Fuhrmann and Cooper, the waiting time is the M/G/1 waiting time plus the
 * equilibrium excess of the setup period,
 *
 *   Wq = lambda E[S^2] / (2 (1 - rho)) + (2 E[U] + lambda E[U^2]) / (2 (1 + lambda E[U])).
 *
 * ARITHMETIC: field. Only the squared coefficient of variation appears, so the
 * computation is exact for T = Rational.
 */
 
#include <string>
 
#include "line/api/qsys/qsys_res.h"
#include "line/lang/distribution.h"
#include "line/num/number.h"
#include "line/util/error.h"
 
namespace line {
namespace qsys {
 
/** MATLAB's [W, Wq, U, Q, res] outputs. */
template <class T>
struct Mg1SetupResult {
    T W;             ///< mean sojourn time
    T Wq;            ///< mean waiting time
    T U;             ///< utilization rho
    T Q;             ///< mean number in system lambda W
    QsysResT<T> res; ///< meanWaitingTime ... probSetup, analyzer
};
 
/**
 * @brief M/G/1 with setup times, setup given by its first two moments.
 *
 * @param lambda arrival rate
 * @param mu     service rate
 * @param cs     coefficient of variation of the service time
 * @param EU     mean setup time
 * @param EU2    second moment of the setup time
 */
template <class T>
Mg1SetupResult<T> qsys_mg1_setup(const T& lambda, const T& mu, const T& cs, const T& EU, const T& EU2) {
    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);
    if (lambda <= zero || mu <= zero || cs < zero)
        throw InputError(ErrorKind("Argument"),
                         "qsys_mg1_setup: lambda and mu must be positive scalars and cs non-negative");
    const T EUsq = T(EU * EU);
    const T tol = T(num_traits<T>::from_double(1e-12) * (EUsq > one ? EUsq : one));
    if (EU < zero || EU2 < T(EUsq - tol))
        throw InputError(ErrorKind("Argument"),
                         "qsys_mg1_setup: setup moments must satisfy EU >= 0 and EU2 >= EU^2");
    const T rho = T(lambda / mu);
    if (rho >= one) resdetail::throw_unstable("qsys_mg1_setup", "rho", num_traits<T>::to_double(rho));
    const T ES2 = T((one + cs * cs) / (mu * mu));
    Mg1SetupResult<T> r;
    r.Wq = T(lambda * ES2 / (two * (one - rho)) + (two * EU + lambda * EU2) / (two * (one + lambda * EU)));
    r.W = T(r.Wq + one / mu);
    r.U = rho;
    r.Q = T(lambda * r.W);
    r.res.set("meanWaitingTime", r.Wq);
    r.res.set("meanSojournTime", r.W);
    r.res.set("meanQueueLength", r.Q);
    r.res.set("meanQueueLengthQ", T(lambda * r.Wq));
    r.res.set("utilization", rho);
    r.res.set("probIdle", T((one - rho) / (one + lambda * EU)));
    r.res.set("probSetup", T((one - rho) * lambda * EU / (one + lambda * EU)));
    r.res.analyzer = "exact:M/G/1-setup";
    return r;
}
 
/**
 * @brief M/G/1 with setup times, setup given as a Distribution (MATLAB passes
 *        it in place of EU and reads its mean and variance).
 */
template <class T>
Mg1SetupResult<T> qsys_mg1_setup(const T& lambda, const T& mu, const T& cs, const lang::Distrib<T>& U) {
    const T EU = lang::dist_moment(U, 1);
    return qsys_mg1_setup(lambda, mu, cs, EU, T(lang::dist_moment(U, 2)));
}
 
}  // namespace qsys
}  // namespace line
 
#endif  // LINE_API_QSYS_QSYS_MG1_SETUP_H