pfqn_mvablo.h

Go to the documentation of this file.

/*
 * Copyright (c) 2012-2026, QORE Lab, Imperial College London
 * All rights reserved.
 */
#ifndef LINE_API_PFQN_PFQN_MVABLO_H
#define LINE_API_PFQN_PFQN_MVABLO_H
 
/**
 * @file
 * @ingroup api_pfqn
 * MVABLO: approximate MVA for closed networks with Blocking-After-Service (BAS).
 *
 * Templated port of matlab/src/api/pfqn/pfqn_mvablo.m, i.e. MVABLO of
 * I. F. Akyildiz, "Mean Value Analysis for Blocking Queueing Networks",
 * IEEE TSE 14(4):418-428, 1988 (Sec. III-IV).
 *
 * Closed single-class network of M stations with exponential service. A queue is
 * a single server with total capacity K[i] (waiting plus in service). Blocking is
 * Blocking-After-Service (BAS): a job completing service at j whose next station
 * i is full stays at j, and j's server stops, until i releases a job. Delay (infinite-server) stations are a LINE extension: no
 * queueing term, never full, but they can be blocked.
 *
 * The population recursion is exact MVA until a mean queue length exceeds its
 * capacity, Q[i] > K[i]+tol. The pass is then repeated with station i marked full
 * (its arriving job's own "+1" term is dropped, eq. 15) and each upstream station
 * j charged the blocking time S[i]*P(j,i)*e[j]/e[i] (eq. 14). Passes repeat until
 * no capacity is exceeded; the full marks and blocking times carry over to the
 * next population, as in the paper.
 *
 * Results follow the paper, so a blocked job is counted at the SOURCE, whereas
 * LINE's CTMC/JMT/LDES count it at the full destination.
 *
 * ARITHMETIC. Only field operations and comparisons, so every T is admitted,
 * Rational included. Capacities are doubles because an unbounded one is +inf.
 */
 
#include <cmath>
#include <cstddef>
#include <limits>
#include <string>
#include <utility>
#include <vector>
 
#include "line/num/number.h"
#include "line/util/error.h"
#include "line/util/matrix.h"
 
namespace line {
namespace pfqn {
 
/** Return value of pfqn_mvablo, mirroring [X,Q,U,R,BT,iter,converged]; vectors are (M). */
template <class T>
struct MvabloResult {
    std::vector<T> X;   ///< throughput of each station
    std::vector<T> Q;   ///< mean number of jobs, blocked jobs at the source
    std::vector<T> U;   ///< utilization X.*S, blocked time excluded
    std::vector<T> R;   ///< mean residence time per visit, blocking time included
    Matrix<T> BT;       ///< (M x M) per-visit blocking time at j caused by a full i
    int iter = 0;       ///< total number of passes over all populations
    bool converged = true;  ///< false if a population exhausted iter_max passes (deadlock)
};
 
/**
 * Runs MVABLO.
 *
 * @param N        population (positive integer)
 * @param S        (M) mean service time per visit
 * @param P        (M x M) station routing matrix, row-stochastic and irreducible
 * @param K        (M) capacity, +inf for unbounded; ignored at delays
 * @param isdelay  (M) true at infinite-server stations, empty for none
 * @param tol      slack on the capacity test Q[i] > K[i]+tol
 * @param iter_max maximum passes at one population
 */
template <class T>
MvabloResult<T> pfqn_mvablo(long N, const std::vector<T>& S, const Matrix<T>& P,
                            const std::vector<double>& K, std::vector<bool> isdelay,
                            double tol = 1e-4, int iter_max = 1000) {
    using nt = num_traits<T>;
    const std::size_t M = S.size();
    const T zero = nt::from_int(0);
    const T one = nt::from_int(1);
    if (isdelay.empty()) isdelay.assign(M, false);
    const std::string caller = "pfqn_mvablo: ";
    if (N <= 0)
        throw InputError(ErrorKind{"Argument"}, caller + "MVABLO requires a positive integer population N.");
    if (P.rows() != M || P.cols() != M || K.size() != M || isdelay.size() != M)
        throw InputError(ErrorKind{"Argument"},
            caller + "MVABLO requires S, K and isdelay of length M and an M x M routing matrix P.");
    double sumK = 0.0;
    bool anyDelay = false;
    for (std::size_t i = 0; i < M; ++i) {
        const double s = nt::to_double(S[i]);
        if (!(s >= 0.0) || !std::isfinite(s))
            throw InputError(ErrorKind{"Argument"}, caller + "MVABLO requires finite nonnegative service times.");
        if (!isdelay[i] && !(K[i] >= 1.0))
            throw InputError(ErrorKind{"Argument"},
                caller + "MVABLO requires every queue capacity K(i) to be at least 1.");
        sumK += K[i];
        anyDelay = anyDelay || isdelay[i];
        double row = 0.0;
        for (std::size_t j = 0; j < M; ++j) {
            const double p = nt::to_double(P(i, j));
            if (p < 0.0)
                throw InputError(ErrorKind{"Argument"}, caller + "MVABLO requires a row-stochastic routing matrix P.");
            row += p;
        }
        if (std::fabs(row - 1.0) > 1e-8)
            throw InputError(ErrorKind{"Argument"}, caller + "MVABLO requires a row-stochastic routing matrix P.");
    }
    std::vector<std::vector<bool>> reach(M, std::vector<bool>(M, false));
    for (std::size_t i = 0; i < M; ++i)
        for (std::size_t j = 0; j < M; ++j) reach[i][j] = i == j || P(i, j) > zero;
    const int squarings = static_cast<int>(std::ceil(std::log2(static_cast<double>(M < 2 ? 2 : M)))) + 1;
    for (int h = 0; h < squarings; ++h) {
        std::vector<std::vector<bool>> next(M, std::vector<bool>(M, false));
        for (std::size_t i = 0; i < M; ++i)
            for (std::size_t j = 0; j < M; ++j) {
                bool r = reach[i][j];
                for (std::size_t l = 0; l < M && !r; ++l) r = reach[i][l] && reach[l][j];
                next[i][j] = r;
            }
        reach.swap(next);
    }
    for (std::size_t i = 0; i < M; ++i)
        for (std::size_t j = 0; j < M; ++j)
            if (!reach[i][j])
                throw InputError(ErrorKind{"Argument"},
                    caller + "MVABLO requires a routing matrix that is irreducible on the stations.");
    if (!anyDelay && static_cast<double>(N) >= sumK)
        throw UnsupportedError(ErrorKind{"Structure"},
            caller + "MVABLO requires N < sum(K) = " + std::to_string(sumK) +
            " to rule out a deadlock of the blocking network; N = " + std::to_string(N) + ".");
 
    // visit ratios e = e*P, normalized at station 1; every output is invariant to the scale.
    // A = P' - I with row 1 replaced by e_1, solved by Gaussian elimination with partial pivoting.
    Matrix<T> A(M, M, zero);
    std::vector<T> e(M, zero);
    for (std::size_t i = 0; i < M; ++i)
        for (std::size_t j = 0; j < M; ++j) A(i, j) = T(P(j, i) - (i == j ? one : zero));
    for (std::size_t j = 0; j < M; ++j) A(0, j) = j == 0 ? one : zero;
    e[0] = one;
    for (std::size_t c = 0; c < M; ++c) {
        std::size_t piv = c;
        for (std::size_t r = c + 1; r < M; ++r)
            if (std::fabs(nt::to_double(A(r, c))) > std::fabs(nt::to_double(A(piv, c)))) piv = r;
        if (piv != c) {
            for (std::size_t j = 0; j < M; ++j) std::swap(A(c, j), A(piv, j));
            std::swap(e[c], e[piv]);
        }
        for (std::size_t r = 0; r < M; ++r) {
            if (r == c || A(r, c) == zero) continue;
            const T f = T(A(r, c) / A(c, c));
            for (std::size_t j = c; j < M; ++j) A(r, j) = T(A(r, j) - f * A(c, j));
            e[r] = T(e[r] - f * e[c]);
        }
    }
    for (std::size_t i = 0; i < M; ++i) e[i] = T(e[i] / A(i, i));
 
    std::vector<T> q(M, zero), z(M, one), t(M, zero);
    Matrix<T> bt(M, M, zero);
    MvabloResult<T> res;
    T lambda = zero;
    for (long n = 1; n <= N; ++n) {
        const std::vector<T> qprev = q;
        int passes = 0;
        while (true) {
            ++passes;
            ++res.iter;
            T denom = zero;
            for (std::size_t i = 0; i < M; ++i) {
                T bsum = zero;
                for (std::size_t j = 0; j < M; ++j) bsum = T(bsum + bt(i, j));
                t[i] = isdelay[i] ? T(S[i] + bsum) : T(S[i] * (z[i] + qprev[i]) + bsum);
                denom = T(denom + e[i] * t[i]);
            }
            lambda = T(nt::from_int(n) / denom);
            std::vector<bool> full(M, false);
            bool anyFull = false;
            for (std::size_t i = 0; i < M; ++i) {
                q[i] = T(e[i] * lambda * t[i]);
                full[i] = !isdelay[i] && std::isfinite(K[i]) && q[i] > nt::from_double(K[i] + tol);
                anyFull = anyFull || full[i];
            }
            if (!anyFull) break;
            if (passes >= iter_max) {
                res.converged = false;
                break;
            }
            for (std::size_t i = 0; i < M; ++i) {
                if (!full[i]) continue;
                z[i] = zero;
                for (std::size_t j = 0; j < M; ++j)
                    if (j != i && P(j, i) > zero) bt(j, i) = T(bt(j, i) + S[i] * P(j, i) * e[j] / e[i]);
            }
        }
    }
    res.X.assign(M, zero);
    res.U.assign(M, zero);
    for (std::size_t i = 0; i < M; ++i) {
        res.X[i] = T(e[i] * lambda);
        res.U[i] = T(res.X[i] * S[i]);
    }
    res.Q = q;
    res.R = t;
    res.BT = bt;
    return res;
}
 
}  // namespace pfqn
}  // namespace line
 
#endif  // LINE_API_PFQN_PFQN_MVABLO_H