qbd_setupdelayoff_mmk.h

Go to the documentation of this file.

/**
 * @file qbd_setupdelayoff_mmk.h
 * @brief Multi-server queue with setup and delayed-off servers.
 *
 * Templated port of matlab/src/api/mam/qbd_setupdelayoff_mmk.m.
 */
 
#ifndef LINE_API_MAM_QBD_SETUPDELAYOFF_MMK_H
#define LINE_API_MAM_QBD_SETUPDELAYOFF_MMK_H
 
#include <cmath>
#include <cstddef>
#include <functional>
#include <map>
#include <string>
#include <vector>
 
#include "line/api/mam/qbd_setupdelayoff.h"
#include "line/api/mc/ctmc_solve.h"
#include "line/lang/lang_types.h"
#include "line/num/number.h"
#include "line/util/error.h"
#include "line/util/matrix.h"
 
namespace line {
namespace mam {
 
/** Mean number present and throughput of the setup/delayed-off pool. */
template <class T>
struct SetupDelayoffMMk {
    T QN;  ///< mean number at the station: queued, waiting for a setup or in service
    T XN;  ///< throughput of the station
};
 
namespace detail {
 
/**
 * Stationary (QN, XN) of the chain of qbd_setupdelayoff_mmk with arrivals at
 * rate lam(n) while n < nmax jobs are present; PTOP receives P(n = nmax).
 */
template <class T>
SetupDelayoffMMk<T> setupdelayoff_mmk_solve(int k, int nmax, const std::function<T(int)>& lam,
                                            const T& mu, const T& alpharate, const T& alphascv,
                                            const T& betarate, const T& betascv, double& ptop) {
    const T zero = num_traits<T>::from_int(0);
    const Matrix<T> Ta = coxian_phase_subgen(alpharate, alphascv);
    const int na = static_cast<int>(Ta.rows());
    std::vector<T> ta(na, zero);
    for (int i = 0; i < na; ++i) {
        T s = zero;
        for (int j = 0; j < na; ++j) s += Ta(i, j);
        ta[i] = T(-s);
    }
    const Matrix<T> Tb = coxian_phase_subgen(betarate, betascv);
    const int nb = static_cast<int>(Tb.rows());
    std::vector<T> tb(nb, zero);
    for (int i = 0; i < nb; ++i) {
        T s = zero;
        for (int j = 0; j < nb; ++j) s += Tb(i, j);
        tb[i] = T(-s);
    }
 
    // state = [q, b, a, s_1..s_na, d_1..d_nb]; off = k - b - a - sum(s) - sum(d)
    const int L = 3 + na + nb;
    std::vector<std::vector<int>> states;
    std::map<std::vector<int>, std::size_t> index;
    states.push_back(std::vector<int>(L, 0));
    index[states[0]] = 0;
    struct Arc {
        std::size_t from, to;
        T rate;
    };
    std::vector<Arc> arcs;
    std::vector<std::pair<std::vector<int>, T>> out;
    auto emit = [&](const std::vector<int>& st, const T& rate, int i1, int dv1, int i2, int dv2) {
        std::vector<int> x = st;
        x[i1] += dv1;
        if (i2 >= 0) x[i2] += dv2;
        out.emplace_back(x, rate);
    };
    for (std::size_t head = 0; head < states.size(); ++head) {
        const std::vector<int> st = states[head];
        out.clear();
        const int q = st[0], b = st[1], a = st[2];
        int ssum = 0, dsum = 0;
        for (int i = 0; i < na; ++i) ssum += st[3 + i];
        for (int j = 0; j < nb; ++j) dsum += st[3 + na + j];
        const int off = k - b - a - ssum - dsum;
        const int n = q + b;
        const T l = n < nmax ? lam(n) : zero;
        if (num_traits<T>::to_double(l) > 0) {
            if (a > 0) {
                emit(st, l, 2, -1, 1, 1);
            } else if (dsum > 0) {
                for (int j = 0; j < nb; ++j)
                    if (st[3 + na + j] > 0)
                        emit(st, T(l * num_traits<T>::from_int(st[3 + na + j]) / num_traits<T>::from_int(dsum)),
                             3 + na + j, -1, 1, 1);
            } else if (off > 0) {
                emit(st, l, 0, 1, 3, 1);
            } else {
                emit(st, l, 0, 1, -1, 0);
            }
        }
        if (b > 0) {
            const T rb = T(num_traits<T>::from_int(b) * mu);
            if (q > 0) emit(st, rb, 0, -1, -1, 0);
            else emit(st, rb, 1, -1, 3 + na, 1);
        }
        for (int i = 0; i < na; ++i) {
            const int si = st[3 + i];
            if (si <= 0) continue;
            for (int i2 = 0; i2 < na; ++i2)
                if (i2 != i && num_traits<T>::to_double(Ta(i, i2)) > 0)
                    emit(st, T(num_traits<T>::from_int(si) * Ta(i, i2)), 3 + i, -1, 3 + i2, 1);
            if (num_traits<T>::to_double(ta[i]) > 0) {
                // a finished setup takes the head of the queue, or stays ACTIVE and idle
                std::vector<int> x = st;
                x[3 + i] -= 1;
                if (q > 0) {
                    x[0] -= 1;
                    x[1] += 1;
                } else {
                    x[2] += 1;
                }
                out.emplace_back(x, T(num_traits<T>::from_int(si) * ta[i]));
            }
        }
        for (int j = 0; j < nb; ++j) {
            const int dj = st[3 + na + j];
            if (dj <= 0) continue;
            for (int j2 = 0; j2 < nb; ++j2)
                if (j2 != j && num_traits<T>::to_double(Tb(j, j2)) > 0)
                    emit(st, T(num_traits<T>::from_int(dj) * Tb(j, j2)), 3 + na + j, -1, 3 + na + j2, 1);
            if (num_traits<T>::to_double(tb[j]) > 0)
                emit(st, T(num_traits<T>::from_int(dj) * tb[j]), 3 + na + j, -1, -1, 0);
        }
        for (const auto& o : out) {
            auto it = index.find(o.first);
            std::size_t to;
            if (it == index.end()) {
                to = states.size();
                index[o.first] = to;
                states.push_back(o.first);
            } else {
                to = it->second;
            }
            arcs.push_back(Arc{head, to, o.second});
        }
    }
    const std::size_t m = states.size();
    Matrix<T> Q(m, m, zero);
    for (const Arc& arc : arcs) {
        Q(arc.from, arc.to) += arc.rate;
        Q(arc.from, arc.from) -= arc.rate;
    }
    const std::vector<T> pi = mc::ctmc_solve(Q);
    double total = 0.0;
    for (std::size_t i = 0; i < m; ++i) total += std::max(0.0, num_traits<T>::to_double(pi[i]));
    SetupDelayoffMMk<T> res;
    res.QN = zero;
    res.XN = zero;
    ptop = 0.0;
    if (total <= 0) return res;
    const T norm = num_traits<T>::from_double(total);
    T QN = zero, B = zero;
    for (std::size_t i = 0; i < m; ++i) {
        if (num_traits<T>::to_double(pi[i]) <= 0) continue;
        const T p = T(pi[i] / norm);
        QN += T(num_traits<T>::from_int(states[i][0] + states[i][1]) * p);
        B += T(num_traits<T>::from_int(states[i][1]) * p);
        if (states[i][0] + states[i][1] == nmax) ptop += num_traits<T>::to_double(p);
    }
    res.QN = QN;
    res.XN = T(mu * B);
    return res;
}
 
}  // namespace detail
 
/**
 * Mean number present and throughput of a k-SERVER queue whose servers each pay
 * a setup when woken from off and count down a delay-off when idle.
 *
 * The M/M/k/setup/delayedoff model of Gandhi, Harchol-Balter and Adan
 * (Performance Evaluation 67(11), 2010), solved exactly as a finite CTMC over its
 * REACHABLE states, under the admission rules of the LDES engine
 * (Solver_ssj_ln.admitWithSetup), which are what a SetupTask means:
 *
 * - an arrival takes an ACTIVE idle server if one exists, else one still
 *   counting down its delay-off (the countdown is cancelled and it pays
 *   nothing), else it queues and wakes ONE off server, if any is off;
 * - a server that completes takes the head of the queue, or starts its delay-off
 *   when the queue is empty, and powers off when that expires;
 * - a setup is never cancelled: a server finishing it takes the head of the
 *   queue, or stays ACTIVE and idle with no countdown when the queue is empty.
 *
 * A state is (queued, busy, active-idle, setup phase counts, delay-off phase
 * counts), the remaining servers being off. Both phases are taken in canonical
 * Coxian form, entered at phase 1 for every SCV; servers in a phase are
 * exchangeable, so an arrival that catches a delay-off server catches one in
 * phase j with probability proportional to the count there, which is exact
 * whenever the delay-off is exponential. With k = 1 the chain is the one
 * `qbd_setupdelayoff_closed` solves.
 *
 * With N finite the population is closed and ARV is the complementary delay Z,
 * so lambda(n) = (N-n)/Z. With N infinite ARV is the Poisson rate and the level
 * is truncated where the probability of the top level falls below
 * GlobalConstants::Zero, doubling the truncation until it does; the queue must
 * then be stable, arv < k*mu.
 *
 * ARITHMETIC. Gated on num_traits<T>::has_transcendental, like the closed twin:
 * the Coxian fit takes a square root.
 *
 * @param k number of servers
 * @param N closed population, or infinity for an open queue
 * @param arv complementary delay Z when N is finite, arrival rate when it is not
 * @param mu service rate of one server
 * @param alpharate rate of the setup phase
 * @param alphascv SCV of the setup phase
 * @param betarate rate of the delay-off phase
 * @param betascv SCV of the delay-off phase
 */
template <class T>
SetupDelayoffMMk<T> qbd_setupdelayoff_mmk(int k, double N, const T& arv, const T& mu,
                                          const T& alpharate, const T& alphascv,
                                          const T& betarate, const T& betascv) {
    static_assert(num_traits<T>::has_transcendental,
                  "qbd_setupdelayoff_mmk requires transcendental arithmetic");
    const T zero = num_traits<T>::from_int(0);
    if (k <= 0 || num_traits<T>::to_double(mu) <= 0)
        throw InputError("qbd_setupdelayoff_mmk needs k >= 1 servers and a positive service rate");
    SetupDelayoffMMk<T> res;
    res.QN = zero;
    res.XN = zero;
    double ptop = 0.0;
    if (std::isfinite(N)) {
        const int pop = static_cast<int>(std::llround(N));
        if (pop <= 0) return res;
        T Z = arv;
        if (num_traits<T>::to_double(Z) < lang::GlobalConstants::FineTol)
            Z = num_traits<T>::from_double(lang::GlobalConstants::FineTol);
        const std::function<T(int)> lam = [pop, Z](int n) {
            return T(num_traits<T>::from_int(pop - n) / Z);
        };
        return detail::setupdelayoff_mmk_solve<T>(k, pop, lam, mu, alpharate, alphascv, betarate,
                                                  betascv, ptop);
    }
    if (num_traits<T>::to_double(arv) <= 0) return res;
    if (num_traits<T>::to_double(arv) >= k * num_traits<T>::to_double(mu))
        throw InputError("qbd_setupdelayoff_mmk: the open queue is unstable, arrival rate " +
                         std::to_string(num_traits<T>::to_double(arv)) + " against capacity " +
                         std::to_string(k * num_traits<T>::to_double(mu)));
    const T rate = arv;
    const std::function<T(int)> lam = [rate](int) { return rate; };
    int nmax = 2 * k + 32;
    while (true) {
        res = detail::setupdelayoff_mmk_solve<T>(k, nmax, lam, mu, alpharate, alphascv, betarate,
                                                 betascv, ptop);
        if (ptop < lang::GlobalConstants::Zero) return res;
        nmax *= 2;
    }
}
 
}  // namespace mam
}  // namespace line
 
#endif  // LINE_API_MAM_QBD_SETUPDELAYOFF_MMK_H