qsys_mg1_exceptional.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_EXCEPTIONAL_H
#define LINE_API_QSYS_QSYS_MG1_EXCEPTIONAL_H
 
/**
 * @file
 * @ingroup api_qsys
 * M/G/1 with exceptional first service: the job that starts a busy period has
 * service law Bf, every other job law B (Welch, 1964).
 *
 * Port of matlab/src/api/qsys/qsys_mg1_exceptional.m. With x = P(empty) =
 * (1 - lambda E[B]) / (1 + lambda E[Bf] - lambda E[B]), the fraction of time
 * serving first jobs is rhof = lambda x E[Bf] and other jobs rho = lambda (1-x) E[B];
 *
 *   Wq = (rhof E[Bf^2]/(2E[Bf]) + rho E[B^2]/(2E[B])) / (1 - lambda E[B]),
 *   W  = Wq + x E[Bf] + (1 - x) E[B].
 *
 * ARITHMETIC: field; exact for T = Rational.
 */
 
#include "line/api/qsys/qsys_res.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 Mg1ExceptionalResult {
    T W;             ///< mean sojourn time
    T Wq;            ///< mean waiting time
    T U;             ///< utilization rhof + rho
    T Q;             ///< mean number in system
    QsysResT<T> res; ///< meanWaitingTime ... rhoOther, analyzer
};
 
/**
 * @brief M/G/1 with exceptional first service.
 *
 * @param lambda arrival rate
 * @param mu     service rate of ordinary jobs
 * @param cs     coefficient of variation of the ordinary service time
 * @param muf    service rate of the job opening a busy period
 * @param csf    coefficient of variation of the exceptional service time
 */
template <class T>
Mg1ExceptionalResult<T> qsys_mg1_exceptional(const T& lambda, const T& mu, const T& cs, const T& muf,
                                             const T& csf) {
    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 || muf <= zero || cs < zero || csf < zero)
        throw InputError(ErrorKind("Argument"),
                         "qsys_mg1_exceptional: lambda, mu, muf must be positive scalars and cs, csf non-negative");
    const T EB = T(one / mu), EB2 = T((one + cs * cs) / (mu * mu));
    const T EBf = T(one / muf), EBf2 = T((one + csf * csf) / (muf * muf));
    if (T(lambda * EB) >= one)
        resdetail::throw_unstable("qsys_mg1_exceptional", "lambda E[B]", num_traits<T>::to_double(T(lambda * EB)));
    const T x = T((one - lambda * EB) / (one + lambda * EBf - lambda * EB));
    const T rhof = T(lambda * x * EBf);
    const T rho = T(lambda * (one - x) * EB);
    Mg1ExceptionalResult<T> r;
    r.Wq = T((rhof * EBf2 / (two * EBf) + rho * EB2 / (two * EB)) / (one - lambda * EB));
    r.W = T(r.Wq + x * EBf + (one - x) * EB);
    r.U = T(rhof + 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", r.U);
    r.res.set("probEmpty", x);
    r.res.set("rhoFirst", rhof);
    r.res.set("rhoOther", rho);
    r.res.analyzer = "exact:M/G/1-exceptional";
    return r;
}
 
}  // namespace qsys
}  // namespace line
 
#endif  // LINE_API_QSYS_QSYS_MG1_EXCEPTIONAL_H