qsys_retrial_isstable.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_RETRIAL_ISSTABLE_H
#define LINE_API_QSYS_QSYS_RETRIAL_ISSTABLE_H
 
/**
 * @file
 * @ingroup api_qsys
 * Stability (positive recurrence) of the M/M/c/c retrial queue.
 *
 * Port of matlab/src/api/qsys/qsys_retrial_isstable.m. Poisson arrivals of
 * rate lambda_k, c exponential servers of rate mu_k, exponential retrials of
 * rate alpha_k (per orbiting customer under LINEAR, of the whole orbit under
 * CONSTANT), join probability p_k, persistence probability q_k, orbit
 * abandonment rate gamma_k and orbit capacity maxOrbit (-1 = unbounded).
 *
 * CONDITIONS:
 *   A finite orbit, gamma > 0, or q < 1 under LINEAR always gives a stable
 *   chain: the orbit then loses customers at a rate that grows with its size.
 *   LINEAR with q = 1 and gamma = 0: stable iff sum_k lambda_k p_k / mu_k < c
 *   (Falin and Templeton 1997, Sec. 2.1 for p = 1).
 *   CONSTANT, single class, gamma = 0: with x = B((lambda+alpha)/mu, c) the
 *   orbit drift is negative iff lambda p x < alpha (1-x) + alpha (1-q) x
 *   (Artalejo, Gomez-Corral and Neuts 2001 for p = q = 1).
 *
 * A multiclass queue under CONSTANT retrial, or with some classes self-
 * limiting and others not, has no known necessary-and-sufficient condition
 * and is refused (LINE:unsupported:Retrial).
 *
 * The margin is infinite when the orbit cannot grow without bound, so the
 * function needs a number type that holds an infinity (double or real).
 */
 
#include <cmath>
#include <cstddef>
#include <limits>
#include <string>
#include <vector>
 
#include "line/api/qsys/qsys_mmcc_retrial_fp.h"
#include "line/api/qsys/qsys_retrial_policy.h"
#include "line/num/number.h"
#include "line/util/error.h"
 
namespace line {
namespace qsys {
 
/** Return value of qsys_retrial_isstable, mirroring the three MATLAB outputs. */
template <class T>
struct RetrialStabilityResult {
    bool isStable = false;  ///< the orbit size is positive recurrent
    T margin;               ///< slack of the condition (positive iff stable; +Inf when the orbit cannot grow)
    std::string reason;     ///< the condition applied, as text
};
 
namespace retrial_isstable_detail {
 
/** A scalar broadcast to K entries, or a K-vector left as it is. */
template <class T>
std::vector<T> expand(const std::vector<T>& v, std::size_t K, const char* name, const T& def) {
    if (v.empty()) return std::vector<T>(K, def);
    if (v.size() == 1) return std::vector<T>(K, v[0]);
    if (v.size() != K)
        throw InputError(std::string("qsys_retrial_isstable: ") + name +
                         " must be a scalar or have one entry per class (" + std::to_string(K) + ").");
    return v;
}
 
}  // namespace retrial_isstable_detail
 
/**
 * @param lambda   arrival rate(s), one per class
 * @param mu       service rate(s) per server (scalar or one per class)
 * @param c        number of servers (no waiting room)
 * @param alpha    retrial rate(s); empty = 0
 * @param policy   LINEAR (default) or CONSTANT
 * @param p        join probability(ies), empty = 1
 * @param q        persistence probability(ies), empty = 1
 * @param gamma    orbit abandonment rate(s) per orbiting customer, empty = 0
 * @param maxOrbit orbit capacity, negative or infinite for unbounded (default -1)
 */
template <class T>
RetrialStabilityResult<T> qsys_retrial_isstable(const std::vector<T>& lambda, const std::vector<T>& mu,
                                                unsigned c, const std::vector<T>& alpha,
                                                RetrialPolicy policy = RetrialPolicy::LINEAR,
                                                const std::vector<T>& p = std::vector<T>(),
                                                const std::vector<T>& q = std::vector<T>(),
                                                const std::vector<T>& gamma = std::vector<T>(),
                                                double maxOrbit = -1.0) {
    static_assert(num_traits<T>::has_transcendental,
                  "qsys_retrial_isstable reports an infinite margin, which needs a type holding Inf");
    using retrial_isstable_detail::expand;
    const T zero = num_traits<T>::from_int(0);
    const T one = num_traits<T>::from_int(1);
    const T inf = num_traits<T>::from_double(std::numeric_limits<double>::infinity());
    const std::size_t K = lambda.size();
    if (K == 0) throw InputError("qsys_retrial_isstable: lambda must have at least one entry.");
    const std::vector<T> muv = expand(mu, K, "mu", one);
    const std::vector<T> pv = expand(p, K, "p", one);
    const std::vector<T> qv = expand(q, K, "q", one);
    const std::vector<T> gv = expand(gamma, K, "gamma", zero);
    const std::vector<T> av = expand(alpha, K, "alpha", zero);
    bool bad = c < 1;
    for (std::size_t k = 0; k < K; ++k)
        bad = bad || lambda[k] < zero || !(muv[k] > zero) || pv[k] < zero || pv[k] > one ||
              qv[k] < zero || qv[k] > one || gv[k] < zero;
    if (bad)
        throw InputError("qsys_retrial_isstable: Invalid input: rates must be positive, c >= 1, and p, q in [0, 1].");
 
    RetrialStabilityResult<T> out;
    if (std::isfinite(maxOrbit) && maxOrbit >= 0) {
        out.isStable = true;
        out.margin = inf;
        out.reason = "finite orbit: the chain is finite";
        return out;
    }
    std::vector<bool> selfLimiting(K);
    bool all = true, anyLim = false;
    for (std::size_t k = 0; k < K; ++k) {
        selfLimiting[k] = gv[k] > zero || (policy == RetrialPolicy::LINEAR && qv[k] < one) ||
                          pv[k] == zero || lambda[k] == zero;
        all = all && selfLimiting[k];
        anyLim = anyLim || selfLimiting[k];
    }
    if (all) {
        out.isStable = true;
        out.margin = inf;
        out.reason = "every orbit loses customers at a rate growing with its size";
        return out;
    }
    if (policy == RetrialPolicy::LINEAR) {
        if (anyLim)
            throw UnsupportedError(ErrorKind{"Retrial"},
                "qsys_retrial_isstable: No known stability condition for a multiclass retrial queue "
                "that mixes persistent classes with self-limiting ones (q < 1 or orbit impatience).");
        for (std::size_t k = 0; k < K; ++k)
            if (!(av[k] > zero)) {
                out.isStable = false;
                out.margin = T(zero - inf);
                out.reason = "LINEAR retrial with zero retrial rate: the orbit never empties";
                return out;
            }
        T load = zero;
        for (std::size_t k = 0; k < K; ++k) load = T(load + lambda[k] * pv[k] / muv[k]);
        out.margin = T(num_traits<T>::from_int(static_cast<long>(c)) - load);
        out.isStable = out.margin > zero;
        out.reason = "LINEAR retrial, persistent orbit: sum(lambda.*p./mu) < c";
        return out;
    }
    // CONSTANT
    if (K > 1)
        throw UnsupportedError(ErrorKind{"Retrial"},
            "qsys_retrial_isstable: No known necessary-and-sufficient stability condition for a "
            "multiclass retrial queue under the CONSTANT policy.");
    if (!(av[0] > zero)) {
        out.isStable = false;
        out.margin = T(zero - inf);
        out.reason = "CONSTANT retrial with zero retrial rate: the orbit never empties";
        return out;
    }
    const T x = erlang_b_recursive(T((lambda[0] + av[0]) / muv[0]), c);
    out.margin = T(av[0] * (one - x) + av[0] * (one - qv[0]) * x - lambda[0] * pv[0] * x);
    out.isStable = out.margin > zero;
    out.reason = "CONSTANT retrial: lambda*p*x < alpha*(1-x) + alpha*(1-q)*x, x = ErlangB((lambda+alpha)/mu, c)";
    return out;
}
 
}  // namespace qsys
}  // namespace line
 
#endif  // LINE_API_QSYS_QSYS_RETRIAL_ISSTABLE_H