qsys_mh21.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_MH21_H
#define LINE_API_QSYS_QSYS_MH21_H
 
/**
 * @file
 * @ingroup api_qsys
 * Exact analysis of the M/H2/1 queue: hyperexponential service with branch
 * probabilities (p1, 1-p1) and rates (mu1, mu2).
 *
 * Port of matlab/src/api/qsys/qsys_mh21.m. The probability generating
 * function of the number in system is rational with two poles z1 < z2 outside
 * the unit disk, the roots of a1 a2 z^2 - (a1 + a2 + a1 a2) z + 1 + a1 + a2 - rho
 * (a_i = lambda/mu_i), so P_n = C1 z1^-n + C2 z2^-n, and the sojourn and
 * waiting times are mixtures of two exponentials of rates lambda (z_i - 1).
 *
 * ARITHMETIC: double, as the reference.
 */
 
#include <algorithm>
#include <cmath>
#include <functional>
#include <vector>
 
#include "line/api/qsys/qsys_mg1.h"
#include "line/api/qsys/qsys_res.h"
#include "line/util/error.h"
 
namespace line {
namespace qsys {
 
/**
 * @brief Exact M/H2/1 queue.
 *
 * @param lambda arrival rate
 * @param p1     probability of the first service phase
 * @param mu1    rate of the first phase
 * @param mu2    rate of the second phase
 * @param o      tPoints, percentiles and nmax options
 */
inline QsysRes qsys_mh21(double lambda, double p1, double mu1, double mu2,
                         const QsysResOptions& o = QsysResOptions()) {
    using namespace resdetail;
    check_percentiles(o.percentiles, "qsys_mh21");
    if (!(lambda > 0.0) || !(mu1 > 0.0) || !(mu2 > 0.0))
        throw Error("qsys_mh21: lambda, mu1 and mu2 must be positive scalars");
    if (!(p1 >= 0.0) || p1 > 1.0) throw Error("qsys_mh21: p1 must lie in [0, 1]");
    const double p2 = 1.0 - p1;
    const double rho = lambda * (p1 / mu1 + p2 / mu2);
    if (rho >= 1.0) throw_unstable("qsys_mh21", "rho", rho);
    const double a1 = lambda / mu1, a2 = lambda / mu2;
    // roots([a1 a2, -(a1 + a2 + a1 a2), 1 + a1 + a2 - rho]), real parts sorted
    const double qa = a1 * a2, qb = -(a1 + a2 + a1 * a2), qc = 1.0 + a1 + a2 - rho;
    double z1, z2;
    const double disc = qb * qb - 4.0 * qa * qc;
    if (disc >= 0.0) {
        const double sq = std::sqrt(disc);
        const double qq = -0.5 * (qb + (qb >= 0.0 ? sq : -sq));
        z1 = qq / qa;
        z2 = qc / qq;
        if (z1 > z2) std::swap(z1, z2);
    } else {
        z1 = z2 = -qb / (2.0 * qa);
    }
    if (std::fabs(z1 - z2) < 1e-10 * std::max(1.0, std::fabs(z2)))
        throw NumericError(ErrorKind{"Root"},
                           "qsys_mh21: the two roots coincide; use qsys_mg1 or qsys_mm1 for this service law");
    const double C1 = (z1 - 1.0) * (1.0 - rho * z2) / (z1 - z2);
    const double C2 = (z2 - 1.0) * (1.0 - rho * z1) / (z2 - z1);
    const double mua = lambda * (z1 - 1.0), mub = lambda * (z2 - 1.0);
    const double pia = C1 * z1 / (z1 - 1.0), pib = C2 * z2 / (z2 - 1.0);
    auto hp = [&](double s) {
        return 1.0 - lambda * (p1 * mu1 / ((mu1 + s) * (mu1 + s)) + p2 * mu2 / ((mu2 + s) * (mu2 + s)));
    };
    const double wa = -(1.0 - rho) / hp(-mua), wb = -(1.0 - rho) / hp(-mub);
 
    const double ES = p1 / mu1 + p2 / mu2;
    const double ES2 = 2.0 * (p1 / (mu1 * mu1) + p2 / (mu2 * mu2));
    const double ES3 = 6.0 * (p1 / (mu1 * mu1 * mu1) + p2 / (mu2 * mu2 * mu2));
    QsysRes res = qsys_mg1_res(lambda, 1.0 / ES, std::sqrt(std::max(ES2 / (ES * ES) - 1.0, 0.0)), ES3);
 
    long nmax = o.nmax;
    if (nmax < 0) nmax = std::max(1L, static_cast<long>(std::ceil(std::log(1e-14) / std::log(1.0 / z1))));
    std::vector<double> P(static_cast<std::size_t>(nmax + 1));
    for (long n = 0; n <= nmax; ++n)
        P[static_cast<std::size_t>(n)] = C1 * std::pow(z1, -static_cast<double>(n)) + C2 * std::pow(z2, -static_cast<double>(n));
    res.set("queueLengthDist", P);
    res.set("roots", std::vector<double>{z1, z2});
    res.set("sojournRates", std::vector<double>{mua, mub});
    res.set("sojournWeights", std::vector<double>{pia, pib});
    res.set("waitingWeights", std::vector<double>{wa, wb});
    const std::function<double(double)> Wccdf = [=](double t) {
        return t < 0.0 ? 1.0 : wa * std::exp(-mua * t) + wb * std::exp(-mub * t);
    };
    const std::function<double(double)> Tccdf = [=](double t) {
        return t < 0.0 ? 1.0 : pia * std::exp(-mua * t) + pib * std::exp(-mub * t);
    };
    if (!o.tPoints.empty()) {
        res.set("tPoints", o.tPoints);
        res.set("waitingTimeCdf", cdf_of(Wccdf, o.tPoints));
        res.set("sojournTimeCdf", cdf_of(Tccdf, o.tPoints));
    }
    if (!o.percentiles.empty()) {
        const double hi = 50.0 / std::min(mua, mub);
        res.set("percentiles", o.percentiles);
        res.set("waitingTimePercentile", pct_of(Wccdf, o.percentiles, hi));
        res.set("sojournTimePercentile", pct_of(Tccdf, o.percentiles, hi));
    }
    res.analyzer = "exact:M/H2/1";
    return res;
}
 
}  // namespace qsys
}  // namespace line
 
#endif  // LINE_API_QSYS_QSYS_MH21_H