qsys_mmc_prio_tardiness.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_MMC_PRIO_TARDINESS_H
#define LINE_API_QSYS_QSYS_MMC_PRIO_TARDINESS_H
 
/**
 * @file
 * @ingroup api_qsys
 * Exact stationary expected tardiness in a multiclass priority M/M/c queue
 * (Hafizoglu, Gel and Keskinocak, INFORMS J. Comput. 25(2), 2013).
 *
 * Port of matlab/src/api/qsys/qsys_mmc_prio_tardiness.m. Classes are listed in
 * priority order (class 1 highest), all share one exponential service rate mu,
 * FCFS within a class, under head-of-line ("np") or preemptive-resume ("pr")
 * priority. A class-i order that finds all servers busy waits q+1 busy periods
 * of an M/M/1 queue with the higher-priority rate and service rate c mu
 * (Theorem 2); under preemption the in-service states solve the tridiagonal
 * system of Theorem 4. Mixing over the stationary state seen on arrival gives
 * the response-time transform f(s); tardiness and its probability are the
 * inverses at d of (s E[T] - 1 + f(s))/s^2 and (1 - f(s))/s. Classes with no
 * higher-priority traffic use a closed form.
 *
 * ARITHMETIC: double. The result is a numerical Laplace inversion.
 */
 
#include <algorithm>
#include <cctype>
#include <cmath>
#include <complex>
#include <cstddef>
#include <limits>
#include <string>
#include <vector>
 
#include "line/api/lti/laplace_invert.h"
#include "line/api/qsys/qsys_res.h"
#include "line/util/error.h"
 
namespace line {
namespace qsys {
 
namespace tarddetail {
 
using Cplx = std::complex<double>;
 
inline bool is_preemptive(const std::string& type, const char* fn) {
    std::string t = type;
    std::transform(t.begin(), t.end(), t.begin(), [](unsigned char ch) { return static_cast<char>(std::tolower(ch)); });
    if (t == "np" || t == "hol" || t == "nonpreemptive" || t == "non-preemptive") return false;
    if (t == "pr" || t == "prpr" || t == "preemptive" || t == "fcfsprprio") return true;
    throw InputError(ErrorKind("Argument"), std::string(fn) + ": unknown priority type '" + type + "', use 'np' or 'pr'");
}
 
inline double erlang_c(unsigned c, double A) {
    double B = 1.0;
    for (unsigned k = 1; k <= c; ++k) B = A * B / (k + A * B);
    return c * B / (c - A * (1 - B));
}
 
/** LST of the M/M/1 busy period with rates lb and cmu, conjugate form. */
inline Cplx busy_lst(Cplx s, double lb, double cmu) {
    const double a = lb + cmu;
    const double b = 2.0 * std::sqrt(lb * cmu);
    return 2.0 * cmu / (s + a + std::sqrt(s + (a - b)) * std::sqrt(s + (a + b)));
}
 
inline Cplx ipow(Cplx z, unsigned k) {
    Cplx out(1.0, 0.0);
    while (k > 0) {
        if (k & 1u) out *= z;
        z *= z;
        k >>= 1;
    }
    return out;
}
 
/** Theorem 4 transforms f_0..f_{c-1}, by the tridiagonal (Thomas) solve. */
inline std::vector<Cplx> pr_state_lst(Cplx s, double lb, double mu, unsigned c) {
    std::vector<Cplx> diag(c), cp(c), dp(c), f(c);
    std::vector<double> sup(c, 0.0), sub(c, 0.0);
    for (unsigned v = 0; v < c; ++v) {
        diag[v] = s + (lb + (v + 1) * mu);
        if (v + 1 < c) sup[v] = -lb;
        if (v > 0) sub[v] = -static_cast<double>(v) * mu;
    }
    diag[c - 1] = s + (lb + c * mu) - lb * busy_lst(s, lb, c * mu);
    cp[0] = sup[0] / diag[0];
    dp[0] = mu / diag[0];
    for (unsigned v = 1; v < c; ++v) {
        const Cplx den = diag[v] - cp[v - 1] * sub[v];
        cp[v] = sup[v] / den;
        dp[v] = (mu - dp[v - 1] * sub[v]) / den;
    }
    f[c - 1] = dp[c - 1];
    for (int v = static_cast<int>(c) - 2; v >= 0; --v) f[v] = dp[v] - cp[v] * f[v + 1];
    return f;
}
 
/** Mean response times S_0..S_{c-1} of the Theorem 4 states. */
inline std::vector<double> pr_means(double lb, double mu, unsigned c) {
    std::vector<double> diag(c), sup(c, 0.0), sub(c, 0.0), rhs(c, 1.0), cp(c), dp(c), x(c);
    for (unsigned v = 0; v < c; ++v) {
        diag[v] = lb + (v + 1) * mu;
        if (v + 1 < c) sup[v] = -lb;
        if (v > 0) sub[v] = -static_cast<double>(v) * mu;
    }
    diag[c - 1] = c * mu;
    rhs[c - 1] = 1 + lb / (c * mu - lb);
    cp[0] = sup[0] / diag[0];
    dp[0] = rhs[0] / diag[0];
    for (unsigned v = 1; v < c; ++v) {
        const double den = diag[v] - cp[v - 1] * sub[v];
        cp[v] = sup[v] / den;
        dp[v] = (rhs[v] - dp[v - 1] * sub[v]) / den;
    }
    x[c - 1] = dp[c - 1];
    for (int v = static_cast<int>(c) - 2; v >= 0; --v) x[v] = dp[v] - cp[v] * x[v + 1];
    return x;
}
 
/** {S, Tard, P} of T = X + B Y, X ~ Exp(mu), Y ~ Exp(theta), B ~ Bernoulli(Pc). */
inline void exp_mix_tard(double mu, double theta, double Pc, double d, double& S, double& T, double& P) {
    S = 1.0 / mu + Pc / theta;
    if (d <= 0) { T = S - d; P = 1.0; return; }
    if (std::isinf(d)) { T = 0.0; P = 0.0; return; }
    const double delta = theta - mu;
    const double E1 = delta == 0 ? d : -std::expm1(-delta * d) / delta;
    const double emd = std::exp(-mu * d);
    T = (1 - Pc) * emd / mu + Pc * (emd / mu + emd * (1 + mu * E1) / theta);
    P = (1 - Pc) * emd + Pc * emd * (1 + mu * E1);
}
 
/** Tard and P by inverting (s S - 1 + F)/s^2 and (1 - F)/s at d. */
inline void invert_tard(const lti::LaplaceFn& F, double S, double d, lti::LaplaceMethod m, double& T, double& P) {
    if (d <= 0) { T = S - d; P = 1.0; return; }
    if (std::isinf(d)) { T = 0.0; P = 0.0; return; }
    T = lti::laplace_invert([&F, S](Cplx s) { return (s * S - 1.0 + F(s)) / (s * s); }, d, m);
    P = lti::laplace_invert([&F](Cplx s) { return (1.0 - F(s)) / s; }, d, m);
    T = std::min(std::max(T, 0.0), S);
    P = std::min(std::max(P, 0.0), 1.0);
}
 
/** Tard and P of T = Exp(mu) + Erlang(K, theta) by nonnegative Poisson series. */
inline void hypo_tard(unsigned K, double mu, double theta, double d, double S, double& T, double& P) {
    if (d <= 0) { T = S - d; P = 1.0; return; }
    if (std::isinf(d)) { T = 0.0; P = 0.0; return; }
    const double w = 1.0 - mu / theta;
    const std::size_t M = w <= 0 ? 0 : static_cast<std::size_t>(std::ceil(std::log(1e-18) / std::log(w)));
    const double y = theta * d;
    const std::size_t J = K + M + 2;
    std::vector<double> p(J), Qc(J + 1, 0.0);
    for (std::size_t j = 0; j < J; ++j) {
        p[j] = std::exp(-y + static_cast<double>(j) * std::log(y) - std::lgamma(static_cast<double>(j) + 1.0));
        Qc[j + 1] = Qc[j] + p[j];
    }
    P = Qc[K];
    T = 0.0;
    for (unsigned j = 0; j < K; ++j) T += (K - j) * p[j];
    double wm = 1.0;
    for (std::size_t m = 0; m <= M; ++m) {
        P += wm * p[K + m];
        T += wm * Qc[K + m + 1];
        wm *= w;
    }
    T /= theta;
    T = std::min(std::max(T, 0.0), S);
    P = std::min(std::max(P, 0.0), 1.0);
}
 
}  // namespace tarddetail
 
/** MATLAB's [Tard, Ptard, res] outputs. */
struct MmcPrioTardinessResult {
    std::vector<double> Tard;  ///< per-class E[(T_i - d_i)^+], class 1 first
    std::vector<double> Ptard; ///< per-class P(T_i > d_i)
    QsysRes res;               ///< meanSojournTime, tardiness, probTardy, dueDate, analyzer
};
 
/**
 * @brief Stationary per-class expected tardiness of a priority M/M/c queue.
 *
 * @param lambda per-class arrival rates, highest priority first
 * @param mu     common service rate
 * @param c      number of servers
 * @param d      one due date, or one per class (infinity = none)
 * @param type   "np" head-of-line or "pr" preemptive resume
 * @param method Laplace inversion method of lti::laplace_invert
 */
inline MmcPrioTardinessResult qsys_mmc_prio_tardiness(const std::vector<double>& lambda, double mu, unsigned c,
                                                      const std::vector<double>& d, const std::string& type = "np",
                                                      const std::string& method = "euler") {
    using namespace tarddetail;
    const char* fn = "qsys_mmc_prio_tardiness";
    const std::size_t R = lambda.size();
    bool bad = R == 0 || !(mu > 0);
    for (double l : lambda)
        if (!(l > 0)) bad = true;
    if (bad) throw InputError(ErrorKind("Argument"), std::string(fn) + ": lambda must be positive and mu a positive scalar");
    if (c < 1) throw InputError(ErrorKind("Argument"), std::string(fn) + ": c must be a positive integer");
    if (d.size() != 1 && d.size() != R)
        throw InputError(ErrorKind("Argument"), std::string(fn) + ": d must be a scalar or one due date per class");
    std::vector<double> dv(R);
    for (std::size_t i = 0; i < R; ++i) {
        dv[i] = d.size() == 1 ? d[0] : d[i];
        if (std::isnan(dv[i]))
            throw InputError(ErrorKind("Argument"), std::string(fn) + ": d must be a scalar or one due date per class");
    }
    const bool preemptive = is_preemptive(type, fn);
    const lti::LaplaceMethod m = lti::laplace_method(method);
    std::vector<double> Lam(R);
    double acc = 0;
    for (std::size_t i = 0; i < R; ++i) Lam[i] = (acc += lambda[i]);
    const double A = Lam[R - 1] / mu;
    if (A >= c) resdetail::throw_unstable(fn, "lambda/(c mu)", A / c);
    const double PcAll = erlang_c(c, A);
    MmcPrioTardinessResult out;
    out.Tard.assign(R, 0.0);
    out.Ptard.assign(R, 0.0);
    std::vector<double> S(R, 0.0);
    const double cmu = c * mu;
    for (std::size_t i = 0; i < R; ++i) {
        const double lb = i == 0 ? 0.0 : Lam[i] - lambda[i];
        const double Li = Lam[i];
        lti::LaplaceFn F;
        if (preemptive) {
            const double a = Li / mu, rho = Li / cmu;
            std::vector<double> pv(c + 1);
            pv[0] = 1.0;
            for (unsigned v = 1; v <= c; ++v) pv[v] = pv[v - 1] * a / v;
            double Z = pv[c] / (1 - rho);
            for (unsigned v = 0; v < c; ++v) Z += pv[v];
            for (double& x : pv) x /= Z;
            if (lb == 0) {
                exp_mix_tard(mu, cmu - Li, pv[c] / (1 - rho), dv[i], S[i], out.Tard[i], out.Ptard[i]);
                continue;
            }
            const std::vector<double> Sv = pr_means(lb, mu, c);
            double s0 = 0;
            for (unsigned v = 0; v < c; ++v) s0 += pv[v] * Sv[v];
            S[i] = s0 + pv[c] / (1 - rho) * Sv[c - 1] + pv[c] / ((cmu - lb) * (1 - rho) * (1 - rho));
            F = [lb, mu, c, cmu, rho, pv](Cplx s) {
                const std::vector<Cplx> f = pr_state_lst(s, lb, mu, c);
                const Cplx g = busy_lst(s, lb, cmu);
                Cplx o(0.0, 0.0);
                for (unsigned v = 0; v < c; ++v) o += pv[v] * f[v];
                return o + pv[c] * g * f[c - 1] / (1.0 - rho * g);
            };
        } else {
            const double sig = Li / cmu;
            if (lb == 0) {
                exp_mix_tard(mu, cmu - Li, PcAll, dv[i], S[i], out.Tard[i], out.Ptard[i]);
                continue;
            }
            S[i] = 1 / mu + PcAll / ((1 - sig) * (cmu - lb));
            F = [lb, mu, cmu, sig, PcAll](Cplx s) {
                const Cplx g = busy_lst(s, lb, cmu);
                return mu / (s + mu) * ((1 - PcAll) + PcAll * (1 - sig) * g / (1.0 - sig * g));
            };
        }
        invert_tard(F, S[i], dv[i], m, out.Tard[i], out.Ptard[i]);
    }
    const std::string tstr = preemptive ? "pr" : "np";
    out.res.set("meanSojournTime", S);
    out.res.set("tardiness", out.Tard);
    out.res.set("probTardy", out.Ptard);
    out.res.set("dueDate", dv);
    out.res.analyzer = "exact:M/M/" + std::to_string(c) + "-prio-" + tstr + "-tardiness";
    return out;
}
 
}  // namespace qsys
}  // namespace line
 
#endif  // LINE_API_QSYS_QSYS_MMC_PRIO_TARDINESS_H