qsys_mmc_prio_eta.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_ETA_H
#define LINE_API_QSYS_QSYS_MMC_PRIO_ETA_H
 
/**
 * @file
 * @ingroup api_qsys
 * Expected tardiness at arrival (ETA) in a multiclass priority M/M/c queue,
 * conditional on the state the arriving order observes (Hafizoglu, Gel and
 * Keskinocak, INFORMS J. Comput. 25(2), 2013).
 *
 * Port of matlab/src/api/qsys/qsys_mmc_prio_eta.m. For type "np", v is the
 * total number of orders in the system and q the number of class <= r orders
 * waiting (used when v >= c); for type "pr", v is the number of class <= r
 * orders in the system and q is ignored. Classes with no higher-priority
 * traffic use the closed form of Theorem 1; the others invert the transform of
 * Theorem 2 ("np") or Theorem 4 and eq. (17) ("pr").
 *
 * ARITHMETIC: double.
 */
 
#include <cmath>
#include <string>
#include <vector>
 
#include "line/api/qsys/qsys_mmc_prio_tardiness.h"
 
namespace line {
namespace qsys {
 
/** MATLAB's [Tard, Ptard, S] outputs, one entry of Tard and Ptard per due date. */
struct MmcPrioEtaResult {
    std::vector<double> Tard;  ///< E[(T - d)^+]
    std::vector<double> Ptard; ///< P(T > d)
    double S = 0.0;            ///< E[T]
};
 
/**
 * @brief Expected tardiness at arrival of a class-r order.
 *
 * @param lambda per-class arrival rates, highest priority first
 * @param mu     common service rate
 * @param c      number of servers
 * @param r      class of the arriving order, 1-based as in MATLAB
 * @param v      orders in the system ("np") or class <= r orders in the system ("pr")
 * @param q      class <= r orders waiting ("np" only)
 * @param d      due dates
 * @param type   "np" or "pr"
 * @param method Laplace inversion method
 */
inline MmcPrioEtaResult qsys_mmc_prio_eta(const std::vector<double>& lambda, double mu, unsigned c, unsigned r,
                                          unsigned v, unsigned q, const std::vector<double>& d,
                                          const std::string& type = "np", const std::string& method = "euler") {
    using namespace tarddetail;
    const char* fn = "qsys_mmc_prio_eta";
    bool bad = lambda.empty() || !(mu > 0);
    for (double l : lambda)
        if (!(l >= 0)) bad = true;
    if (bad) throw InputError(ErrorKind("Argument"), std::string(fn) + ": lambda must be nonnegative and mu a positive scalar");
    if (c < 1) throw InputError(ErrorKind("Argument"), std::string(fn) + ": c must be a positive integer");
    if (r < 1 || r > lambda.size()) throw InputError(ErrorKind("Argument"), std::string(fn) + ": r must be a class index of lambda");
    for (double x : d)
        if (std::isnan(x)) throw InputError(ErrorKind("Argument"), std::string(fn) + ": d must not be NaN");
    const bool preemptive = is_preemptive(type, fn);
    const lti::LaplaceMethod m = lti::laplace_method(method);
    double lb = 0;
    for (unsigned j = 0; j + 1 < r; ++j) lb += lambda[j];
    const double cmu = c * mu;
    if (lb >= cmu) resdetail::throw_unstable(fn, "higher-priority load", lb / cmu);
    unsigned K;
    if (preemptive)
        K = v + 1 > c ? v - c + 1 : 0;
    else
        K = v + 1 <= c ? 0 : q + 1;
    MmcPrioEtaResult out;
    out.Tard.assign(d.size(), 0.0);
    out.Ptard.assign(d.size(), 0.0);
    if (lb == 0 || (K == 0 && !preemptive)) {
        out.S = 1 / mu + K / (cmu - lb);
        for (std::size_t n = 0; n < d.size(); ++n) hypo_tard(K, mu, cmu, d[n], out.S, out.Tard[n], out.Ptard[n]);
        return out;
    }
    lti::LaplaceFn F;
    if (preemptive) {
        const std::vector<double> Sv = pr_means(lb, mu, c);
        if (v + 1 <= c) {
            out.S = Sv[v];
            F = [lb, mu, c, v](Cplx s) { return pr_state_lst(s, lb, mu, c)[v]; };
        } else {
            out.S = Sv[c - 1] + K / (cmu - lb);
            F = [lb, mu, c, cmu, K](Cplx s) { return pr_state_lst(s, lb, mu, c)[c - 1] * ipow(busy_lst(s, lb, cmu), K); };
        }
    } else {
        out.S = 1 / mu + K / (cmu - lb);
        F = [lb, mu, cmu, K](Cplx s) { return mu / (s + mu) * ipow(busy_lst(s, lb, cmu), K); };
    }
    for (std::size_t n = 0; n < d.size(); ++n) invert_tard(F, out.S, d[n], m, out.Tard[n], out.Ptard[n]);
    return out;
}
 
}  // namespace qsys
}  // namespace line
 
#endif  // LINE_API_QSYS_QSYS_MMC_PRIO_ETA_H