qsys_mmcc_retrial_fp.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_MMCC_RETRIAL_FP_H
#define LINE_API_QSYS_QSYS_MMCC_RETRIAL_FP_H
 
/**
 * @file
 * @ingroup api_qsys
 * Fixed-point approximation for the M/M/c/c retrial queue.
 *
 * Port of matlab/src/api/qsys/qsys_mmcc_retrial_fp.m. Blocked customers join an
 * orbit and retry; under the assumption that the retrial rate is small relative
 * to the service rate the superposition of fresh and retrial arrivals is
 * approximated by a Poisson stream of rate lambda + r, with r the solution of
 *
 *     r = (lambda + r) B((lambda + r)/mu, c),
 *
 * B being Erlang's loss formula. Cohen (1957); Phung-Duc, "Retrial Queueing
 * Models: A Survey on Theory and Applications" (2019), eq. (1).
 *
 * With non-persistent customers (Phung-Duc 2019, Sec. 3.2), a blocked fresh
 * customer joins the orbit with probability p and a blocked retrial customer
 * returns to it with probability q. The orbit then sends retrials at the rate
 *
 *     r = lambda p B / (1 - q B),   B = B((lambda + r)/mu, c),
 *
 * iterated as r <- lambda p B + q B r, which is Cohen's iteration when
 * p = q = 1 (so the iteration count is unchanged there). A fresh customer is
 * lost with probability B(1-p) + B p (1-q) B/(1 - q B), zero when p = q = 1.
 *
 * The reference evaluates B by the rational recursion
 * B_k = a B_{k-1}/(k + a B_{k-1}), which is exactly what line::lossn::erlang_b
 * declines to do (it goes through log/exp and says so). This port keeps the
 * reference's recursion, so B itself is a rational function of a and needs no
 * transcendental arithmetic. The static_assert is nonetheless present because
 * the OUTER iteration is a fixed point tested against a tolerance: it converges
 * geometrically but not in a finite number of field operations, and at Rational
 * the iterates would grow without bound in representation size while never
 * meeting an exact stopping rule.
 *
 * The iteration is a monotone increasing map of r started from 0 and bounded
 * above by lambda/(1-B) at the fixed point, so it converges from below and the
 * successive-difference test is a genuine stopping criterion rather than a
 * heuristic. The reference does NOT report non-convergence; when maxiter is
 * exhausted it silently returns the last iterate, and the port reports the
 * iteration count so the caller can tell the two apart, exactly as the third
 * MATLAB output niter does.
 */
 
#include <cstddef>
#include <string>
 
#include "line/num/number.h"
#include "line/util/error.h"
 
namespace line {
namespace qsys {
 
/** Return value of qsys_mmcc_retrial_fp, mirroring the four MATLAB outputs. */
template <class T>
struct MmccRetrialFpResult {
    T blockingProbability;  ///< B((lambda + r)/mu, c)
    T retrialRate;          ///< r, the extra arrival rate contributed by the orbit
    T lossProbability;      ///< probability that a fresh customer is eventually lost
    std::size_t iterations;  ///< iterations performed, = maxiter when unconverged
    bool converged;          ///< whether the successive-difference test was met
};
 
/**
 * Erlang's loss formula B(a, c) by the numerically stable rational recursion
 * B_0 = 1, B_k = a B_{k-1}/(k + a B_{k-1}). Exact in any field: no logs, no
 * factorials, no cancellation.
 */
template <class T>
T erlang_b_recursive(const T& a, unsigned c) {
    T b = num_traits<T>::from_int(1);
    for (unsigned i = 1; i <= c; ++i) b = a * b / (num_traits<T>::from_int(static_cast<long>(i)) + a * b);
    return b;
}
 
/**
 * M/M/c/c with retrials by the Cohen fixed point.
 *
 * @param lambda  fresh arrival rate
 * @param mu      service rate of one server
 * @param c       number of servers, which is also the capacity
 * @param tol     stopping tolerance on |r_{k+1} - r_k|
 * @param maxiter iteration cap
 * @param p       probability that a blocked fresh customer joins the orbit (default 1)
 * @param q       probability that a blocked retrial customer stays in the orbit (default 1)
 */
template <class T>
MmccRetrialFpResult<T> qsys_mmcc_retrial_fp(const T& lambda, const T& mu, unsigned c, const T& tol,
                                            std::size_t maxiter,
                                            const T& p = num_traits<T>::from_int(1),
                                            const T& q = num_traits<T>::from_int(1)) {
    static_assert(num_traits<T>::has_transcendental,
                  "qsys_mmcc_retrial_fp is a tolerance-terminated fixed point");
    const T zero = num_traits<T>::from_int(0);
    if (lambda <= zero) throw InputError("qsys_mmcc_retrial_fp: arrival rate must be positive");
    if (mu <= zero) throw InputError("qsys_mmcc_retrial_fp: service rate must be positive");
    if (c == 0) throw InputError("qsys_mmcc_retrial_fp: at least one server is required");
    if (tol <= zero) throw InputError("qsys_mmcc_retrial_fp: tolerance must be positive");
    if (maxiter == 0) throw InputError("qsys_mmcc_retrial_fp: maxiter must be positive");
    const T one = num_traits<T>::from_int(1);
    if (p < zero || p > one || q < zero || q > one)
        throw InputError("qsys_mmcc_retrial_fp: The join and persistence probabilities P and Q must be in [0, 1].");
    if (!(q < one) && p > zero &&
        !(lambda * p / mu < num_traits<T>::from_int(static_cast<long>(c)))) {
        // a persistent orbit grows without bound: no fixed point exists
        throw InputError(ErrorKind("Unstable"),
            "qsys_mmcc_retrial_fp: unstable retrial queue: lambda*p/mu = " +
            std::to_string(num_traits<T>::to_double(T(lambda * p / mu))) + " >= c = " +
            std::to_string(c) + ", the persistent orbit grows without bound");
    }
 
    T r = zero;
    std::size_t it = 0;
    bool converged = false;
    for (; it < maxiter;) {
        ++it;
        const T a = (lambda + r) / mu;
        const T b = erlang_b_recursive(a, c);
        // r = lambda*p*B + q*B*r, the same iteration as Cohen's r = (lambda + r)*B,
        // to which it reduces when p = q = 1
        const T rnew = lambda * p * b + q * b * r;
        const T gap = num_abs(T(rnew - r));
        r = rnew;
        if (gap < tol) {
            converged = true;
            break;
        }
    }
 
    MmccRetrialFpResult<T> out;
    out.retrialRate = r;
    out.blockingProbability = erlang_b_recursive(T((lambda + r) / mu), c);
    const T B = out.blockingProbability;
    out.lossProbability = T(B * (one - p) + B * p * (one - q) * B / (one - q * B));
    out.iterations = it;
    out.converged = converged;
    return out;
}
 
/** qsys_mmcc_retrial_fp with the reference defaults tol = 1e-10, maxiter = 10000. */
template <class T>
MmccRetrialFpResult<T> qsys_mmcc_retrial_fp(const T& lambda, const T& mu, unsigned c) {
    return qsys_mmcc_retrial_fp(lambda, mu, c, T(num_traits<T>::from_double(1e-10)),
                                static_cast<std::size_t>(10000));
}
 
}  // namespace qsys
}  // namespace line
 
#endif  // LINE_API_QSYS_QSYS_MMCC_RETRIAL_FP_H