pfqn_oqnrc.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_OQNRC_H
#define LINE_API_PFQN_PFQN_OQNRC_H
 
/**
 * @file
 * @ingroup api_pfqn
 * OQNRC: approximate MVA for open queueing networks with restricted capacity.
 *
 * Port of matlab/src/api/pfqn/pfqn_oqnrc.m, i.e. the algorithm of
 * R. Buitenhek, G.-J. van Houtum, H. Zijm, "Approximate MVA algorithms for open
 * queueing networks with restricted capacity", University of Twente, 1997.
 *
 * An open network whose inner system holds at most Ncap[g] jobs of the classes
 * in token pool g. Blocked arrivals wait outside (WAITQ, the paper) or are lost
 * (DROP). Each pool is a closed token chain with a load-dependent
 * synchronization station, nu(1) = lam_g/q_g, q_g = 1 - lam_g/TH_g^C(N_g)
 * (q_g = 1 under DROP), solved by the approximate MVA recursion of the paper's
 * Appendix over the population lattice. FCFS multi-server stations use the
 * residual-service approximation, PS and INF are exact.
 *
 * ARITHMETIC. double only: the recursion is an approximation and the WAITQ
 * external queue needs a geometric tail.
 */
 
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <cstdio>
#include <string>
#include <vector>
 
#include "line/util/error.h"
#include "line/util/matrix.h"
 
namespace line {
namespace pfqn {
 
/** Result of pfqn_oqnrc. Station matrices are M x R, class vectors length R. */
struct OqnrcResult {
    Matrix<double> Q, T, U, R;   ///< queue length, throughput, utilization, response time per visit
    std::vector<double> X;       ///< carried throughput of each class
    std::vector<double> Qext;    ///< mean class-r arrivals waiting outside (WAITQ)
    std::vector<double> Ploss;   ///< loss probability of each class (DROP)
    int iter = 0;                ///< multi-pool fixed point iterations
};
 
namespace oqnrc_detail {
 
/** Approximate MVA over the population lattice 0..Npop (pool 0 fastest). */
struct Lattice {
    const std::vector<double>& lamg;
    const std::vector<std::vector<double>>& D;
    const std::vector<std::vector<double>>& V;
    const std::vector<std::vector<double>>& ES;
    const std::vector<std::vector<double>>& ESrem;
    const Matrix<double>& S;
    const std::vector<int>& pool;
    const std::vector<double>& c;
    const std::vector<bool>& isdelay;
    const std::vector<bool>& isps;
    std::vector<std::vector<double>> TH;
    std::vector<std::size_t> stride;
    std::vector<std::vector<double>> ETr;
 
    void run(const std::vector<int>& Npop, const std::vector<bool>& hasSync, const std::vector<double>& q,
             bool wantET) {
        const std::size_t M = D.size();
        const std::size_t G = Npop.size();
        const std::size_t Rn = pool.size();
        stride.assign(G, 0);
        std::size_t P = 1;
        for (std::size_t g = 0; g < G; ++g) {
            stride[g] = P;
            P *= static_cast<std::size_t>(Npop[g] + 1);
        }
        int cmax = 1;
        for (std::size_t m = 0; m < M; ++m)
            if (!isdelay[m]) cmax = std::max(cmax, static_cast<int>(c[m]));
        const std::size_t cm1 = static_cast<std::size_t>(cmax);
        // flat storage: pm[(idx*M + m)*cmax + l], Qh/ELh[(idx*M + m)*G + h]
        std::vector<double> pm(P * M * cm1, 0.0), Qh(P * M * G, 0.0), ELh(P * M * G, 0.0);
        for (std::size_t m = 0; m < M; ++m) pm[m * cm1] = 1.0;
        std::vector<double> sp0(P * G, 1.0), sQ(P * G, 0.0), sL(P * G, 0.0);
        TH.assign(P, std::vector<double>(G, 0.0));
        ETr.assign(M, std::vector<double>(Rn, 0.0));
        auto qh = [&](std::size_t i, std::size_t m, std::size_t h) -> double& { return Qh[(i * M + m) * G + h]; };
        auto el = [&](std::size_t i, std::size_t m, std::size_t h) -> double& { return ELh[(i * M + m) * G + h]; };
        auto pp = [&](std::size_t i, std::size_t m, std::size_t l) -> double& { return pm[(i * M + m) * cm1 + l]; };
        std::vector<int> n(G, 0);
        for (std::size_t idx = 1; idx < P; ++idx) {
            std::size_t rr = idx;
            for (std::size_t gg = G; gg-- > 0;) {
                n[gg] = static_cast<int>(rr / stride[gg]);
                rr -= static_cast<std::size_t>(n[gg]) * stride[gg];
            }
            for (std::size_t g = 0; g < G; ++g) {
                if (n[g] == 0) continue;
                const std::size_t prev = idx - stride[g];
                double den = 0.0;
                for (std::size_t m = 0; m < M; ++m) {
                    if (V[m][g] <= 0) continue;
                    double base = 0.0, fac = 1.0;
                    if (isdelay[m]) {
                        // exact: no queueing
                    } else if (isps[m]) {
                        double tot = 0;
                        for (std::size_t h = 0; h < G; ++h) tot += qh(prev, m, h) + el(prev, m, h);
                        fac = 1 + tot / c[m];
                    } else {
                        for (std::size_t h = 0; h < G; ++h)
                            base += qh(prev, m, h) * ESrem[m][h] + el(prev, m, h) * ES[m][h] / c[m];
                    }
                    den += V[m][g] * (base + fac * ES[m][g]);
                    if (wantET && idx == P - 1)
                        for (std::size_t r = 0; r < Rn; ++r)
                            if (static_cast<std::size_t>(pool[r]) == g) ETr[m][r] = base + fac * S(m, r);
                }
                if (hasSync[g])
                    den += (sL[prev * G + g] + sQ[prev * G + g] + sp0[prev * G + g] * q[g]) / lamg[g];
                TH[idx][g] = n[g] / den;
            }
            for (std::size_t m = 0; m < M; ++m) {
                if (isdelay[m]) continue;
                const int cm = static_cast<int>(c[m]);
                for (std::size_t h = 0; h < G; ++h) {
                    if (n[h] == 0 || D[m][h] <= 0) continue;
                    const std::size_t prev = idx - stride[h];
                    double Qp = 0, ELp = 0;
                    for (std::size_t k = 0; k < G; ++k) {
                        Qp += qh(prev, m, k);
                        ELp += el(prev, m, k);
                    }
                    const double coef = D[m][h] * TH[idx][h];
                    qh(idx, m, h) = coef / cm * (Qp + pp(prev, m, static_cast<std::size_t>(cm - 1)));
                    el(idx, m, h) = coef / cm * (ELp + Qp);
                    for (int l = 1; l < cm; ++l)
                        pp(idx, m, static_cast<std::size_t>(l)) += coef / l * pp(prev, m, static_cast<std::size_t>(l - 1));
                }
                double s = 0;
                for (int l = 1; l < cm; ++l) s += pp(idx, m, static_cast<std::size_t>(l));
                for (std::size_t h = 0; h < G; ++h) s += qh(idx, m, h);
                pp(idx, m, 0) = std::max(0.0, 1 - s);
            }
            for (std::size_t g = 0; g < G; ++g) {
                if (hasSync[g] && n[g] > 0) {
                    const std::size_t prev = idx - stride[g];
                    const double ETs = (sL[prev * G + g] + sQ[prev * G + g] + sp0[prev * G + g] * q[g]) / lamg[g];
                    sQ[idx * G + g] =
                        std::min(1.0, TH[idx][g] / lamg[g] * (sQ[prev * G + g] + sp0[prev * G + g] * q[g]));
                    sL[idx * G + g] = TH[idx][g] * ETs;
                    sp0[idx * G + g] = 1 - sQ[idx * G + g];
                }
            }
        }
    }
};
 
}  // namespace oqnrc_detail
 
/**
 * Runs OQNRC.
 *
 * @param lam      external arrival rate of each class (length R)
 * @param A        arrival rate of class r at station m (M x R)
 * @param S        mean service time (M x R)
 * @param SCV      squared coefficient of variation of the service time (M x R)
 * @param nservers servers per station, +inf at a delay (length M)
 * @param isps     true at processor-sharing stations (length M)
 * @param pool     token pool of each class, 0-based (length R)
 * @param Ncap     tokens of each pool (length G)
 * @param waitq    true when blocked arrivals wait, false when they are lost
 * @param tol      multi-pool fixed point tolerance
 * @param iter_max multi-pool fixed point iteration cap
 */
inline OqnrcResult pfqn_oqnrc(const std::vector<double>& lam, const Matrix<double>& A, const Matrix<double>& S,
                              const Matrix<double>& SCV, const std::vector<double>& nservers,
                              const std::vector<bool>& isps, const std::vector<int>& pool,
                              const std::vector<int>& Ncap, bool waitq, double tol, int iter_max) {
    const std::size_t M = A.rows();
    const std::size_t R = A.cols();
    const std::size_t G = Ncap.size();
    for (std::size_t g = 0; g < G; ++g)
        if (Ncap[g] < 1)
            throw InputError(ErrorKind{"Argument"},
                             "pfqn_oqnrc: Every token pool needs a positive integer number of tokens.");
    std::vector<double> lamg(G, 0.0);
    for (std::size_t r = 0; r < R; ++r) lamg[static_cast<std::size_t>(pool[r])] += lam[r];
    for (std::size_t g = 0; g < G; ++g)
        if (!(lamg[g] > 0))
            throw InputError(ErrorKind{"Argument"}, "pfqn_oqnrc: Every token pool needs a positive arrival rate.");
    std::vector<bool> isdelay(M, false);
    for (std::size_t m = 0; m < M; ++m) isdelay[m] = std::isinf(nservers[m]);
    std::vector<std::vector<double>> D(M, std::vector<double>(G, 0.0)), V = D, ES = D, ESrem = D;
    for (std::size_t g = 0; g < G; ++g)
        for (std::size_t m = 0; m < M; ++m) {
            double asum = 0, as = 0, as2 = 0;
            for (std::size_t r = 0; r < R; ++r) {
                if (static_cast<std::size_t>(pool[r]) != g) continue;
                const double a = A(m, r), s = S(m, r);
                asum += a;
                as += a * s;
                as2 += a * (1 + SCV(m, r)) * s * s;
            }
            if (asum <= 0) continue;
            V[m][g] = asum / lamg[g];
            D[m][g] = as / lamg[g];
            ES[m][g] = D[m][g] / V[m][g];
            if (!isdelay[m] && ES[m][g] > 0) {
                const double c = nservers[m];
                const double ES2 = as2 / asum;
                ESrem[m][g] = (c - 1) / (c + 1) * ES[m][g] / c + 2 / (c + 1) * (1 / c) * ES2 / (2 * ES[m][g]);
            }
        }
    oqnrc_detail::Lattice cs{lamg, D, V, ES, ESrem, S, pool, nservers, isdelay, isps, {}, {}, {}};
 
    std::vector<double> q(G, 1.0);
    std::vector<std::vector<double>> THC(G);
    int it = 0;
    if (waitq) {
        for (it = 1; it <= iter_max; ++it) {
            std::vector<double> qnew = q;
            for (std::size_t g = 0; g < G; ++g) {
                std::vector<bool> hasSync(G, true);
                hasSync[g] = false;
                cs.run(Ncap, hasSync, q, false);
                std::size_t base = 0;
                for (std::size_t h = 0; h < G; ++h)
                    if (h != g) base += static_cast<std::size_t>(Ncap[h]) * cs.stride[h];
                THC[g].assign(static_cast<std::size_t>(Ncap[g]), 0.0);
                for (int k = 1; k <= Ncap[g]; ++k)
                    THC[g][static_cast<std::size_t>(k - 1)] = cs.TH[base + static_cast<std::size_t>(k) * cs.stride[g]][g];
                const double thN = THC[g].back();
                if (lamg[g] >= thN * (1 - 1e-12)) {
                    char buf[320];
                    std::snprintf(buf, sizeof(buf),
                                  "pfqn_oqnrc: Token pool %zu is unstable: its arrival rate %g is not below the "
                                  "throughput %g of the region holding N=%d tokens, so the queue of waiting "
                                  "arrivals grows without bound.",
                                  g + 1, lamg[g], thN, Ncap[g]);
                    throw NumericError(ErrorKind{"Instability"}, buf);
                }
                qnew[g] = 1 - lamg[g] / thN;
            }
            double delta = 0;
            for (std::size_t g = 0; g < G; ++g) delta = std::max(delta, std::fabs(qnew[g] - q[g]));
            q = qnew;
            if (G == 1 || delta < tol) break;
        }
        if (it > iter_max) it = iter_max;
    }
    cs.run(Ncap, std::vector<bool>(G, true), q, true);
    const std::vector<double> THN = cs.TH.back();
 
    std::vector<double> Xg(G, 0.0), Qg(G, 0.0), Pg(G, 0.0);
    for (std::size_t g = 0; g < G; ++g) {
        if (waitq) {
            Xg[g] = lamg[g];
            double pik = 1.0, sum = 0.0;
            for (int k = 0; k < Ncap[g]; ++k) {
                sum += pik;
                pik *= lamg[g] / THC[g][static_cast<std::size_t>(k)];
            }
            const double rho = lamg[g] / THC[g].back();
            const double Z = sum + pik / (1 - rho);
            Qg[g] = pik * rho / ((1 - rho) * (1 - rho)) / Z;
        } else {
            Xg[g] = THN[g];
            Pg[g] = std::max(0.0, 1 - THN[g] / lamg[g]);
        }
    }
    OqnrcResult res;
    res.Q = Matrix<double>(M, R, 0.0);
    res.T = Matrix<double>(M, R, 0.0);
    res.U = Matrix<double>(M, R, 0.0);
    res.R = Matrix<double>(M, R, 0.0);
    res.X.assign(R, 0.0);
    res.Qext.assign(R, 0.0);
    res.Ploss.assign(R, 0.0);
    for (std::size_t r = 0; r < R; ++r) {
        const std::size_t g = static_cast<std::size_t>(pool[r]);
        for (std::size_t m = 0; m < M; ++m) {
            const double t = A(m, r) * Xg[g] / lamg[g];
            const double qq = THN[g] * A(m, r) / lamg[g] * cs.ETr[m][r];
            res.T(m, r) = t;
            res.Q(m, r) = qq;
            res.U(m, r) = t * S(m, r) / (isdelay[m] ? 1.0 : nservers[m]);
            res.R(m, r) = t > 0 ? qq / t : 0.0;
        }
        res.X[r] = lam[r] * Xg[g] / lamg[g];
        res.Qext[r] = Qg[g] * lam[r] / lamg[g];
        res.Ploss[r] = Pg[g];
    }
    res.iter = it;
    return res;
}
 
}  // namespace pfqn
}  // namespace line
 
#endif  // LINE_API_PFQN_PFQN_OQNRC_H