"""
Quasi-Birth-Death (QBD) Process Utilities.
Native Python implementations for QBD matrix computations including
rate matrix R computation and QBD block construction.
Key algorithms:
qbd_R: Rate matrix via successive substitutions
qbd_R_logred: Rate matrix via logarithmic reduction
qbd_rg: Compute R and G matrices
References:
Original MATLAB: matlab/src/api/mam/qbd_*.m
Latouche & Ramaswami, "Introduction to Matrix Analytic Methods
in Stochastic Modeling", 1999
"""
import numpy as np
from numpy.linalg import LinAlgError
from scipy import linalg
from typing import Tuple, Optional
from dataclasses import dataclass
[docs]
@dataclass
class QBDResult:
"""Result of QBD analysis."""
R: np.ndarray # Rate matrix R
G: Optional[np.ndarray] = None # Rate matrix G
U: Optional[np.ndarray] = None # U matrix
eta: Optional[float] = None # Caudal characteristic
[docs]
def qbd_R(B: np.ndarray, L: np.ndarray, F: np.ndarray,
iter_max: int = 100000, tol: float = 1e-12) -> np.ndarray:
"""
Compute QBD rate matrix R using successive substitutions.
Solves the matrix quadratic equation:
R^2 * A_{-1} + R * A_0 + A_1 = 0
where A_{-1} = B, A_0 = L, A_1 = F.
Args:
B: Backward transition block A_{-1}
L: Local transition block A_0
F: Forward transition block A_1
iter_max: Maximum iterations (default: 100000)
tol: Convergence tolerance (default: 1e-12)
Returns:
Rate matrix R
References:
Original MATLAB: matlab/src/api/mam/qbd_R.m
"""
B = np.asarray(B, dtype=np.float64)
L = np.asarray(L, dtype=np.float64)
F = np.asarray(F, dtype=np.float64)
try:
L_inv = linalg.inv(L)
except LinAlgError:
L_inv = linalg.pinv(L)
Fil = F @ L_inv
BiL = B @ L_inv
R = -Fil
Rprime = -Fil - R @ R @ BiL
for _ in range(iter_max):
R = Rprime
Rprime = -Fil - R @ R @ BiL
if linalg.norm(R - Rprime, 1) <= tol:
break
return Rprime
[docs]
def qbd_R_logred(B: np.ndarray, L: np.ndarray, F: np.ndarray,
iter_max: int = 1000, tol: float = 1e-14) -> np.ndarray:
"""
Compute QBD rate matrix R using logarithmic reduction.
Uses the logarithmic reduction algorithm which has quadratic
convergence compared to linear convergence of successive substitutions.
Args:
B: Backward transition block A_{-1}
L: Local transition block A_0
F: Forward transition block A_1
iter_max: Maximum iterations (default: 1000)
tol: Convergence tolerance (default: 1e-14)
Returns:
Rate matrix R
References:
Original MATLAB: matlab/src/api/mam/qbd_R_logred.m
Latouche & Ramaswami, Ch. 8
"""
B = np.asarray(B, dtype=np.float64)
L = np.asarray(L, dtype=np.float64)
F = np.asarray(F, dtype=np.float64)
n = L.shape[0]
try:
L_inv = linalg.inv(-L)
except LinAlgError:
L_inv = linalg.pinv(-L)
A = L_inv @ F # A_1
C = L_inv @ B # A_{-1}
H = A.copy()
T = C.copy()
for _ in range(iter_max):
# Compute (I - AC - CA)^{-1}
I = np.eye(n)
M = I - A @ C - C @ A
try:
M_inv = linalg.inv(M)
except LinAlgError:
M_inv = linalg.pinv(M)
# Update
A_new = M_inv @ A @ A
C_new = M_inv @ C @ C
H_new = H + T @ M_inv @ A
if linalg.norm(H_new - H, 1) <= tol:
H = H_new
break
A = A_new
C = C_new
T = T @ M_inv @ (A + C)
H = H_new
# R = H * (-L)
R = H @ (-L)
return R
[docs]
def qbd_rg(B: np.ndarray, L: np.ndarray, F: np.ndarray,
method: str = 'logred', iter_max: int = 1000,
tol: float = 1e-14) -> QBDResult:
"""
Compute both R and G matrices for a QBD process.
G is the minimal non-negative solution to:
A_1 * G^2 + A_0 * G + A_{-1} = 0
R is the minimal non-negative solution to:
R^2 * A_{-1} + R * A_0 + A_1 = 0
Args:
B: Backward transition block A_{-1}
L: Local transition block A_0
F: Forward transition block A_1
method: 'logred' or 'successive' (default: 'logred')
iter_max: Maximum iterations
tol: Convergence tolerance
Returns:
QBDResult with R, G, U, and eta (caudal characteristic)
References:
Original MATLAB: matlab/src/api/mam/qbd_rg.m
"""
B = np.asarray(B, dtype=np.float64)
L = np.asarray(L, dtype=np.float64)
F = np.asarray(F, dtype=np.float64)
n = L.shape[0]
# Compute R
if method == 'logred':
R = qbd_R_logred(B, L, F, iter_max, tol)
else:
R = qbd_R(B, L, F, iter_max, tol)
# Compute G using the relation: G = A_{-1} * (R*A_{-1} + A_0)^{-1}
# Or alternatively iterate: G_{n+1} = -(A_0 + A_1*G_n^2)^{-1} * A_{-1}
try:
L_inv = linalg.inv(-L)
except LinAlgError:
L_inv = linalg.pinv(-L)
A = L_inv @ F
C = L_inv @ B
# Use logarithmic reduction for G
G = C.copy()
T = A.copy()
for _ in range(iter_max):
I = np.eye(n)
M = I - A @ C - C @ A
try:
M_inv = linalg.inv(M)
except LinAlgError:
M_inv = linalg.pinv(M)
A_new = M_inv @ A @ A
C_new = M_inv @ C @ C
G_new = G + T @ M_inv @ C
if linalg.norm(G_new - G, 1) <= tol:
G = G_new
break
A = A_new
C = C_new
T = T @ M_inv @ (A + C)
G = G_new
G = G @ (-L)
# Compute U = A_0 + A_1 * G
U = L + F @ G
# Caudal characteristic (spectral radius of R)
try:
eigvals = linalg.eigvals(R)
eta = np.max(np.abs(eigvals))
except:
eta = None
return QBDResult(R=R, G=G, U=U, eta=eta)
[docs]
def qbd_blocks_mapmap1(D0_arr: np.ndarray, D1_arr: np.ndarray,
D0_srv: np.ndarray, D1_srv: np.ndarray
) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
"""
Construct QBD blocks for a MAP/MAP/1 queue.
Builds the backward (B), local (L), and forward (F) transition
blocks for the QBD representation of a MAP/MAP/1 queue.
Args:
D0_arr: Arrival MAP D0 matrix
D1_arr: Arrival MAP D1 matrix
D0_srv: Service MAP D0 matrix
D1_srv: Service MAP D1 matrix
Returns:
Tuple of (B, L, F) QBD blocks
References:
Original MATLAB: matlab/src/api/mam/qbd_mapmap1.m
"""
D0_arr = np.asarray(D0_arr, dtype=np.float64)
D1_arr = np.asarray(D1_arr, dtype=np.float64)
D0_srv = np.asarray(D0_srv, dtype=np.float64)
D1_srv = np.asarray(D1_srv, dtype=np.float64)
na = D0_arr.shape[0]
ns = D0_srv.shape[0]
# Forward transitions (arrivals): F = D1_arr \otimes I_ns
F = np.kron(D1_arr, np.eye(ns))
# Backward transitions (service completions): B = I_na \otimes D1_srv
B = np.kron(np.eye(na), D1_srv)
# Local transitions: L = D0_arr \otimes I_ns + I_na \otimes D0_srv
L = np.kron(D0_arr, np.eye(ns)) + np.kron(np.eye(na), D0_srv)
return B, L, F
[docs]
def qbd_bmapbmap1(
MAPa: Tuple[np.ndarray, np.ndarray],
pbatcha: np.ndarray,
MAPs: Tuple[np.ndarray, np.ndarray]
) -> Tuple[np.ndarray, np.ndarray, list, np.ndarray, list]:
"""
Compute QBD blocks for a BMAP/BMAP/1 queue.
Constructs the QBD (Quasi-Birth-Death) transition blocks for a
BMAP/BMAP/1 queue with batch arrivals.
Args:
MAPa: Arrival process MAP as (D0, D1)
pbatcha: Probability distribution of batch sizes (array of length maxbatch)
MAPs: Service process MAP as (D0, D1)
Returns:
Tuple of (A0, A_1, A1_list, B0, B1_list) where:
A0: Local transition block
A_1: Downward transition block
A1_list: List of upward transition blocks for each batch size
B0: Initial boundary local block
B1_list: List of boundary upward blocks for each batch size
References:
Original MATLAB: matlab/src/api/mam/qbd_bmapbmap1.m
"""
D0_arr, D1_arr = MAPa
D0_srv, D1_srv = MAPs
D0_arr = np.asarray(D0_arr, dtype=np.float64)
D1_arr = np.asarray(D1_arr, dtype=np.float64)
D0_srv = np.asarray(D0_srv, dtype=np.float64)
D1_srv = np.asarray(D1_srv, dtype=np.float64)
pbatcha = np.asarray(pbatcha, dtype=np.float64)
na = D0_arr.shape[0]
ns = D0_srv.shape[0]
maxbatch = len(pbatcha)
# Build upward transition blocks for each batch size
A1_list = []
for b in range(maxbatch):
A1_b = np.kron(D1_arr * pbatcha[b], np.eye(ns))
A1_list.append(A1_b)
# Local transitions: A0 = D0_arr \otimes I_ns + I_na \otimes D0_srv (Kronecker sum)
A0 = np.kron(D0_arr, np.eye(ns)) + np.kron(np.eye(na), D0_srv)
# Downward transitions: A_1 = I_na \otimes D1_srv
A_1 = np.kron(np.eye(na), D1_srv)
# Boundary blocks (for level 0)
# MATLAB: B0 = krons(MAPa{1}, eye(ns)) = kron(D0_arr, I_ns) + kron(I_na, I_ns)
B0 = np.kron(D0_arr, np.eye(ns)) + np.eye(na * ns)
B1_list = []
for b in range(maxbatch):
B1_b = np.kron(D1_arr * pbatcha[b], np.eye(ns))
B1_list.append(B1_b)
return A0, A_1, A1_list, B0, B1_list
[docs]
def qbd_mapmap1(
MAPa: Tuple[np.ndarray, np.ndarray],
MAPs: Tuple[np.ndarray, np.ndarray],
util: Optional[float] = None
) -> Tuple[float, float, float, np.ndarray, np.ndarray, Optional[float],
np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray,
Tuple[np.ndarray, np.ndarray]]:
"""
Analyze a MAP/MAP/1 queue using QBD methods.
Solves a MAP/MAP/1 queue using Quasi-Birth-Death process methods,
computing throughput, queue length, utilization, and other metrics.
Args:
MAPa: Arrival process MAP as (D0, D1)
MAPs: Service process MAP as (D0, D1)
util: Optional target utilization to scale service rate
Returns:
Tuple of (XN, QN, UN, pqueue, R, eta, G, A_1, A0, A1, U, MAPs_scaled)
where:
XN: System throughput
QN: Mean queue length
UN: Utilization
pqueue: Queue length distribution
R: Rate matrix R
eta: Caudal characteristic (spectral radius of R)
G: Rate matrix G
A_1: Downward transition block
A0: Local transition block
A1: Upward transition block
U: Matrix U
MAPs_scaled: Scaled service process
References:
Original MATLAB: matlab/src/api/mam/qbd_mapmap1.m
"""
from ..mam import map_scale, map_lambda
D0_arr, D1_arr = MAPa
D0_srv, D1_srv = MAPs
D0_arr = np.asarray(D0_arr, dtype=np.float64)
D1_arr = np.asarray(D1_arr, dtype=np.float64)
D0_srv = np.asarray(D0_srv, dtype=np.float64)
D1_srv = np.asarray(D1_srv, dtype=np.float64)
na = D0_arr.shape[0]
ns = D0_srv.shape[0]
# Scale service process if target utilization provided
if util is not None:
lambda_a = map_lambda(D0_arr, D1_arr)
D0_srv, D1_srv = map_scale(D0_srv, D1_srv, util / lambda_a)
lambda_a = map_lambda(D0_arr, D1_arr)
lambda_s = map_lambda(D0_srv, D1_srv)
actual_util = lambda_a / lambda_s
# Build QBD blocks
A1 = np.kron(D1_arr, np.eye(ns)) # Forward (arrivals)
A0 = np.kron(D0_arr, np.eye(ns)) + np.kron(np.eye(na), D0_srv) # Local
A_1 = np.kron(np.eye(na), D1_srv) # Backward (services)
A0bar = np.kron(D0_arr, np.eye(ns)) # Boundary local
# Solve QBD using Cyclic Reduction (matching MATLAB: QBD_CR(A_1, A0, A1))
from line_solver.lib.thirdparty.smc import qbd_cr, qbd_pi
cr_result = qbd_cr(A_1, A0, A1)
G = cr_result['G']
R = cr_result['R']
U = cr_result['U']
# Compute caudal characteristic (eta = spectral radius of R)
eta_val = np.max(np.abs(linalg.eigvals(R)))
# Compute queue length distribution using QBD_pi
# MATLAB: pqueue = QBD_pi(A_1, A0bar, R, 'MaxNumComp', 1e2)
pi_flat = qbd_pi(A_1, A0bar, R, max_num_comp=100)
n_phases = na * ns
num_levels = len(pi_flat) // n_phases
pqueue = pi_flat.reshape(num_levels, n_phases)
# Retry with more components if needed (MATLAB line 75-77)
if np.sum(np.sum(pqueue[1:, :], axis=1)) < actual_util * 0.99:
pi_flat = qbd_pi(A_1, A0bar, R, max_num_comp=20000)
num_levels = len(pi_flat) // n_phases
pqueue = pi_flat.reshape(num_levels, n_phases)
# Compute performance measures (matching MATLAB lines 79-88)
if na == 1 and ns == 1:
UN = 1.0 - pqueue[0, 0]
QN = float(np.arange(pqueue.shape[0]) @ pqueue.flatten())
else:
UN = 1.0 - np.sum(pqueue[0, :])
QN = float(np.arange(pqueue.shape[0]) @ np.sum(pqueue, axis=1))
XN = lambda_a
MAPs_scaled = (D0_srv, D1_srv)
return XN, QN, UN, pqueue, R, eta_val, G, A_1, A0, A1, U, MAPs_scaled
[docs]
def qbd_raprap1(
RAPa: Tuple[np.ndarray, np.ndarray],
RAPs: Tuple[np.ndarray, np.ndarray],
util: Optional[float] = None
) -> Tuple[float, float, float, np.ndarray, np.ndarray, float,
np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
"""
Analyze a RAP/RAP/1 queue using QBD methods.
Solves a RAP/RAP/1 queue (Rational Arrival Process) using QBD methods,
computing throughput, queue length, utilization, and other metrics.
References:
N. G. Bean and B. F. Nielsen, "Quasi-Birth-and-Death Processes with
Rational Arrival Process Components", Stochastic Models, 26(3), 2010,
pp. 309-334. The analysis rests on the prediction-process
interpretation of a RAP due to Asmussen and Bladt, which is what
allows a QBD argument to be carried over to matrices that are not
nonnegative. The same prediction process underlies the
conditional-vector RAP sampler in RAP.sample.
Args:
RAPa: Arrival process RAP as (H0, H1)
RAPs: Service process RAP as (H0, H1)
util: Optional target utilization to scale service rate
Returns:
Tuple of (XN, QN, UN, pqueue, R, eta, G, B, L, F) where:
XN: System throughput
QN: Mean queue length
UN: Utilization
pqueue: Queue length distribution
R: Rate matrix R
eta: Caudal characteristic
G: Rate matrix G
B: Backward transition block
L: Local transition block
F: Forward transition block
References:
Original MATLAB: matlab/src/api/mam/qbd_raprap1.m
"""
from ..mam import map_scale, map_lambda
H0_arr, H1_arr = RAPa
H0_srv, H1_srv = RAPs
H0_arr = np.asarray(H0_arr, dtype=np.float64)
H1_arr = np.asarray(H1_arr, dtype=np.float64)
H0_srv = np.asarray(H0_srv, dtype=np.float64)
H1_srv = np.asarray(H1_srv, dtype=np.float64)
na = H0_arr.shape[0]
ns = H0_srv.shape[0]
# Scale service process if target utilization provided
if util is not None:
lambda_a = map_lambda(H0_arr, H1_arr)
H0_srv, H1_srv = map_scale(H0_srv, H1_srv, util / lambda_a)
lambda_a = map_lambda(H0_arr, H1_arr)
# QBD blocks of the RAP/RAP/1 queue: a level is the number in system, a
# phase is the (arrival,service) RAP phase pair, ARRIVAL phase major (the
# phase index is a*ns + s, the ordering produced by kron(RAPa, I_ns) and
# kron(I_na, RAPs)). The level rises on an arrival (F), falls on a service
# completion (B), and the two RAPs evolve independently between level
# changes (L, the Kronecker sum of the two hidden generators). At level 0
# the queue is empty, so no service completion can occur and only the
# arrival RAP evolves: the boundary local block is kron(C0, I) and the
# boundary up-block is F.
B_blk = np.kron(np.eye(na), H1_srv) # backward, service completions
L_blk = np.kron(H0_arr, np.eye(ns)) + np.kron(np.eye(na), H0_srv) # local
F_blk = np.kron(H1_arr, np.eye(ns)) # forward, arrivals
B1 = np.kron(H0_arr, np.eye(ns)) # boundary local
# Theorem 7 of Bean and Nielsen (2010), in qbd_rap: G from the quadratic
# matrix equation, U = L + F*G, R = F*inv(-U), and pi0 from the boundary
# equation pi0*(B1 + R*B) = 0 normalised by pi0*inv(I-R)*e = 1.
core = qbd_rap(F_blk, L_blk, B_blk, F_blk, B1, 0)
R = core.R
G = core.G
# Level series, truncated when the accumulated mass reaches 1-1e-10 or at
# 100 level vectors. qbd_rap returns the exact mean queue length in closed
# form, but this function keeps reporting the truncated series because
# pqueue is the documented return value and the three codebases must cut
# the tail at exactly the same point, otherwise the means disagree at the
# 1e-9 level.
n_phases = na * ns
max_num_comp = 100
levels = [core.pi0.reshape(1, n_phases)]
sumpi = float(np.sum(levels[0]))
numit = 1
while sumpi < 1.0 - 1e-10 and numit < 1 + max_num_comp:
levels.append(levels[numit - 1] @ R)
numit += 1
sumpi += float(np.sum(levels[numit - 1]))
num_levels = len(levels)
pqueue = np.vstack(levels)
# Compute performance measures matching MATLAB qbd_raprap1.m
eta = core.spr
if na == 1 and ns == 1:
UN = 1.0 - pqueue[0, 0]
else:
UN = 1.0 - np.sum(pqueue[0, :])
QN = float(np.arange(pqueue.shape[0]) @ np.sum(pqueue, axis=1))
XN = lambda_a
# Return B, L, F using original naming convention
B = B_blk
L = L_blk
F = F_blk
return XN, QN, UN, pqueue, R, eta, G, B, L, F
[docs]
@dataclass
class QbdRapResult:
"""
Result of the equilibrium analysis of a QBD with RAP components.
Attributes:
levelProb: Marginal level probabilities, levels 0..numLevels
QN: Mean queue length, computed exactly as pi0*R*inv(I-R)^2*e
R: Rate matrix R = A0*inv(-U)
G: Matrix G solving A0*G^2 + A1*G + A2 = 0
U: Matrix U = A1 + A0*G
spr: Spectral radius Sp(R); positive recurrent iff Sp(R) < 1
pqueue: (numLevels+1) x m array whose n-th row is the level vector pi_n
pi0: Level-0 vector pi_0, the boundary vector of Theorem 7
"""
levelProb: np.ndarray
QN: float
R: np.ndarray
G: np.ndarray
U: np.ndarray
spr: float
pqueue: np.ndarray
pi0: np.ndarray
[docs]
def qbd_rap_g(A0: np.ndarray, A1: np.ndarray, A2: np.ndarray,
block_scale: float) -> np.ndarray:
"""
Solve A0*G^2 + A1*G + A2 = 0 for the matrix G.
Uses the exact rank-one closed form when A2 has rank one, and otherwise
natural functional iteration as a warm start followed by Newton's method on
the Sylvester-form Jacobian. Never returns an unconverged iterate.
Args:
A0: Level-up block
A1: Local block
A2: Level-down block
block_scale: Scale of the blocks, used to set the residual tolerance
Returns:
The matrix G
Raises:
ValueError: if G cannot be computed to roundoff level
"""
m = A1.shape[0]
e = np.ones((m, 1))
res_tol = 1e-10 * block_scale
# Rank-one A2 = u*v admits the closed form G = e*v/(v*e), which is the case
# deliberately chosen in the example of Bean and Nielsen (2010) precisely so
# that G is available a priori. G is idempotent and, by conservativity,
# (A0+A1)*e = -A2*e = -u*(v*e), so (A0+A1)*e*v/(v*e) = -u*v = -A2.
_, s2, Vh2 = linalg.svd(A2)
if len(s2) > 1 and s2[0] > 0 and s2[1] <= 1e-10 * s2[0]:
v = Vh2[0, :].reshape(1, m)
ve = float((v @ e).item())
if abs(ve) < 1e-12 * np.max(np.abs(v)):
raise ValueError('A2 has rank one but its right factor v satisfies v*e = 0, '
'so the closed form G = e*v/(v*e) is undefined.')
G = e @ (v / ve)
res = float(linalg.norm(A0 @ G @ G + A1 @ G + A2, 'fro'))
if res > res_tol:
raise ValueError('The rank-one closed form for G leaves a residual '
'||A0*G^2 + A1*G + A2||_F = %g, which is above the roundoff '
'level %g.' % (res, res_tol))
return G
# General case. Natural functional iteration G <- inv(-A1)*(A2 + A0*G^2)
# gives a warm start; it is the standard Markovian iteration but has no
# convergence guarantee for blocks that are not nonnegative.
try:
mA1inv = linalg.inv(-A1)
except LinAlgError:
raise ValueError('The local block A1 is singular, the iteration for G cannot be started.')
G = np.zeros((m, m))
for _ in range(200):
Gnew = mA1inv @ (A2 + A0 @ G @ G)
if not np.all(np.isfinite(Gnew)):
break
step = float(linalg.norm(Gnew - G, 'fro'))
G = Gnew
if step <= 1e-14 * max(float(linalg.norm(G, 'fro')), 1.0):
break
if not np.all(np.isfinite(G)):
G = np.zeros((m, m))
# Newton's method on F(G) = A0*G^2 + A1*G + A2. The derivative in the
# direction H is (A0*G+A1)*H + A0*H*G, a Sylvester-type operator solved here
# through its Kronecker expansion, with column-major vectorisation.
Im = np.eye(m)
for _ in range(100):
res_mat = A0 @ G @ G + A1 @ G + A2
if linalg.norm(res_mat, 'fro') <= res_tol:
break
J = np.kron(Im, A0 @ G + A1) + np.kron(G.T, A0)
try:
h = linalg.solve(J, -res_mat.reshape(-1, order='F'))
except (LinAlgError, ValueError):
break
H = h.reshape(m, m, order='F')
G = G + H
if not np.all(np.isfinite(G)):
break
res = np.inf
ge_err = np.inf
if np.all(np.isfinite(G)):
res = float(linalg.norm(A0 @ G @ G + A1 @ G + A2, 'fro'))
ge_err = float(np.max(np.abs(G @ e - e)))
if not (res <= res_tol) or ge_err > 1e-8:
raise ValueError(
'Could not compute the matrix G for this QBD with RAP components: residual '
'||A0*G^2 + A1*G + A2||_F = %g against a tolerance of %g, and ||G*e-e||_inf = %g. '
'The blocks are not nonnegative, so neither the functional iteration nor '
"Newton's method is guaranteed to converge, and the justification of algorithms "
'for G in this setting is left as an open problem in Section 6 of N. G. Bean and '
'B. F. Nielsen, "Quasi-Birth-and-Death Processes with Rational Arrival Process '
'Components", Stochastic Models, 26(3), 2010, pp. 309-334. Supply a model with a '
'rank-one A2, for which G is available in closed form.' % (res, res_tol, ge_err))
return G
[docs]
def qbd_rap(A0: np.ndarray, A1: np.ndarray, A2: np.ndarray,
B0: Optional[np.ndarray] = None, B1: Optional[np.ndarray] = None,
numLevels: int = 20) -> QbdRapResult:
"""
Equilibrium analysis of a Quasi-Birth-and-Death process with Rational
Arrival Process (RAP) components.
The process is specified directly by its level-independent blocks
(A0,A1,A2) and its boundary blocks (B0,B1), where A0 drives level
increases, A2 drives level decreases and A1 the within-level evolution.
Unlike a Markovian QBD the blocks need not be nonnegative: they are only
required to be conservative, (A0+A1+A2)*e = 0, and to define a genuine RAP
through the prediction-process interpretation. This makes qbd_rap strictly
more general than qbd_raprap1, which builds a product-space QBD from two
INDEPENDENT RAPs; here the arrival process and the sequence of service
times may be driven from a shared phase space and therefore be
cross-correlated.
This is the block-level core of the RAP QBD family. qbd_raprap1 is the thin
wrapper over it that builds the product-space blocks of two independent
RAPs; callers with a coupled model must call qbd_rap directly because no
product form exists to factor out.
Algorithm (Theorem 7 of the reference):
1. Solve A0*G^2 + A1*G + A2 = 0 for G.
2. U = A1 + A0*G.
3. R = A0*inv(-U).
4. Find the row vector pihat0 with pihat0*(B1 + R*A2) = 0, pihat0*e = 1.
5. pi0 = K*pihat0 with K chosen so that pi0*inv(I-R)*e = 1.
6. pi_n = pi0*R^n, and the marginal level probability is pi_n*e.
The process is positive recurrent iff Sp(R) < 1 and step 4 has a solution.
Computation of G: the blocks are not nonnegative, so the probabilistic
iterations used for Markovian QBDs (logarithmic reduction, cyclic
reduction) carry no convergence guarantee here, and the paper explicitly
leaves the general case open ("The issue of justifying algorithms for the
evaluation of the matrix G for such processes has not been undertaken",
Section 6). See qbd_rap_g: the rank-one closed form is used when it
applies, otherwise functional iteration followed by Newton's method, and an
unconverged G is never returned.
References:
N. G. Bean and B. F. Nielsen, "Quasi-Birth-and-Death Processes with
Rational Arrival Process Components", Stochastic Models, 26(3), 2010,
pp. 309-334 (DTU technical report IMM-2007-20). The argument rests on
the prediction-process interpretation of a RAP due to Asmussen and
Bladt, which is what allows a QBD argument to be carried over to
matrices that are not nonnegative; the same prediction process
underlies the conditional-vector RAP sampler in RAP.sample.
Original MATLAB: matlab/src/api/mam/qbd_rap.m
Args:
A0: Level-up block (m x m)
A1: Local block (m x m)
A2: Level-down block (m x m)
B0: Boundary level-up block, defaults to A0
B1: Boundary local block, defaults to A1
numLevels: Highest level reported, default 20
Returns:
QbdRapResult with levelProb, QN, R, G, U, spr, pqueue and pi0
Raises:
ValueError: if the blocks are not conservative, if the process is not
positive recurrent, or if G cannot be computed
"""
A0 = np.asarray(A0, dtype=np.float64)
A1 = np.asarray(A1, dtype=np.float64)
A2 = np.asarray(A2, dtype=np.float64)
B0 = A0.copy() if B0 is None else np.asarray(B0, dtype=np.float64)
B1 = A1.copy() if B1 is None else np.asarray(B1, dtype=np.float64)
m = A1.shape[0]
for blk, name in ((A0, 'A0'), (A1, 'A1'), (A2, 'A2'), (B0, 'B0'), (B1, 'B1')):
if blk.shape != (m, m):
raise ValueError('All QBD blocks must be square and of the same order; '
'%s has shape %s.' % (name, blk.shape))
if numLevels < 0 or int(numLevels) != numLevels:
raise ValueError('numLevels must be a nonnegative integer.')
numLevels = int(numLevels)
e = np.ones((m, 1))
I = np.eye(m)
block_scale = max(1.0, float(linalg.norm(A0, 'fro')),
float(linalg.norm(A1, 'fro')), float(linalg.norm(A2, 'fro')))
# Conservativity of the repeating portion, (A0+A1+A2)*e = 0. This is the
# RAP analogue of the generator row-sum condition and every step below
# relies on it.
cons_a = float(np.max(np.abs((A0 + A1 + A2) @ e)))
if cons_a > 1e-8 * block_scale:
raise ValueError('The repeating blocks are not conservative: ||(A0+A1+A2)*e||_inf = %g. '
'A QBD with RAP components requires (A0+A1+A2)*e = 0.' % cons_a)
cons_b = float(np.max(np.abs((B0 + B1) @ e)))
if cons_b > 1e-8 * block_scale:
raise ValueError('The boundary blocks are not conservative: ||(B0+B1)*e||_inf = %g. '
'A QBD with RAP components requires (B0+B1)*e = 0 at level 0.' % cons_b)
# Step 1: matrix G.
G = qbd_rap_g(A0, A1, A2, block_scale)
# Steps 2 and 3: U and R.
U = A1 + A0 @ G
try:
R = A0 @ linalg.inv(-U)
except LinAlgError:
raise ValueError('The matrix U = A1 + A0*G is singular, R = A0*inv(-U) does not exist.')
# Corollary 8(i): positive recurrence.
spr = float(np.max(np.abs(linalg.eigvals(R))))
# The threshold carries a 1e-12 margin: at the null-recurrent boundary
# Sp(R) equals 1 in exact arithmetic but rounds to either side, and the
# three codebases must agree on rejecting it. Any model within 1e-12 of the
# boundary has an unbounded queue regardless.
if spr >= 1.0 - 1e-12:
raise ValueError('The process is not positive recurrent: Sp(R) = %.15g >= 1 '
'(Corollary 8 of Bean and Nielsen, 2010).' % spr)
# Step 4: boundary vector, pihat0*(B1 + R*A2) = 0 normalised to pihat0*e = 1.
# The left null space is extracted from the SVD of V' with a relative
# tolerance: V is only singular up to the accuracy with which R was
# computed, so an absolute tolerance is too strict here.
V = B1 + R @ A2
_, sv, Wh = linalg.svd(V.T)
null_tol = 1e-8 * max(float(sv[0]), 1.0)
if float(sv[-1]) > null_tol:
raise ValueError('The boundary equation x*(B1 + R*A2) = 0 has no nontrivial solution '
'(smallest singular value %g against tolerance %g), so the process is '
'not positive recurrent (Corollary 8(ii) of Bean and Nielsen, 2010).'
% (float(sv[-1]), null_tol))
if m > 1 and float(sv[-2]) <= null_tol:
raise ValueError('The boundary equation x*(B1 + R*A2) = 0 has a solution space of '
'dimension greater than one, the equilibrium vector is not unique.')
pihat0 = Wh[-1, :].reshape(1, m)
den = float((pihat0 @ e).item())
if abs(den) < 1e-12 * np.max(np.abs(pihat0)):
raise ValueError('The boundary vector cannot be normalised, x*e = 0.')
pihat0 = pihat0 / den
# Step 5: level-0 vector.
ImRinv = linalg.inv(I - R)
K = 1.0 / float((pihat0 @ ImRinv @ e).item())
pi0 = K * pihat0
# Consistency of the supplied boundary up-block: the level-0 balance
# equation pi0*B0 + pi1*A1 + pi2*A2 = 0 must hold with pi_n = pi0*R^n.
bal = pi0 @ B0 + pi0 @ R @ A1 + pi0 @ R @ R @ A2
bal_norm = float(np.max(np.abs(bal)))
if bal_norm > 1e-8 * block_scale * max(float(np.max(np.abs(pi0))), 1.0):
raise ValueError('The boundary block B0 is inconsistent with the repeating blocks: '
'||pi0*B0 + pi1*A1 + pi2*A2||_inf = %g. The level-0 balance equation '
'of Theorem 7 requires pi0*(B0-A0) = 0.' % bal_norm)
# Step 6: level vectors and marginal level distribution.
pqueue = np.zeros((numLevels + 1, m))
pin = pi0.copy()
for n in range(numLevels + 1):
pqueue[n, :] = pin[0, :]
pin = pin @ R
levelProb = pqueue @ e
levelProb = levelProb.reshape(-1)
# Exact mean queue length, sum_n n*pi0*R^n*e = pi0*R*inv(I-R)^2*e.
QN = float((pi0 @ R @ ImRinv @ ImRinv @ e).item())
return QbdRapResult(levelProb=levelProb, QN=QN, R=R, G=G, U=U, spr=spr, pqueue=pqueue,
pi0=pi0)
[docs]
def qbd_setupdelayoff(
lambda_val: float,
mu: float,
alpharate: float,
alphascv: float,
betarate: float,
betascv: float
) -> float:
"""
Analyze queue with setup delay and turn-off phases.
Performs queue-length analysis for a queueing system with setup
delay (warm-up) and turn-off periods using QBD methods.
The system operates as follows:
1. When empty and job arrives, server enters setup phase
2. After setup, server becomes active and serves jobs
3. When queue empties, server enters turn-off phase
4. After turn-off, server becomes idle
Args:
lambda_val: Arrival rate
mu: Service rate
alpharate: Rate of setup delay phase
alphascv: Squared coefficient of variation for setup delay
betarate: Rate of turn-off phase
betascv: Squared coefficient of variation for turn-off period
Returns:
Average queue length QN
References:
Original MATLAB: matlab/src/api/mam/qbd_setupdelayoff.m
"""
# Try to import AcyclicPHFromMeansAndSCVs, but it may not exist
try:
from line_solver.lib.thirdparty.butools.ph.baseph import AcyclicPHFromMeansAndSCVs
HAS_ACYCLIC_PH = True
except ImportError:
HAS_ACYCLIC_PH = False
AcyclicPHFromMeansAndSCVs = None
# Fit PH distributions for setup and turn-off phases
alpha_D0 = None
if HAS_ACYCLIC_PH:
try:
alpha_ph = AcyclicPHFromMeansAndSCVs([1.0 / alpharate], [alphascv])
alpha_D0 = alpha_ph[1] # Subgenerator matrix
alpha_D1 = -alpha_D0 @ np.ones((alpha_D0.shape[0], 1))
except:
alpha_D0 = None
if alpha_D0 is None:
# Fallback to Erlang approximation
na = max(1, int(round(1.0 / alphascv)))
rate_a = na * alpharate
alpha_D0 = np.diag([-rate_a] * na)
for i in range(na - 1):
alpha_D0[i, i + 1] = rate_a
alpha_D1 = np.zeros((na, 1))
alpha_D1[-1, 0] = rate_a
beta_D0 = None
if HAS_ACYCLIC_PH:
try:
beta_ph = AcyclicPHFromMeansAndSCVs([1.0 / betarate], [betascv])
beta_D0 = beta_ph[1]
beta_D1 = -beta_D0 @ np.ones((beta_D0.shape[0], 1))
except:
beta_D0 = None
if beta_D0 is None:
# Fallback to Erlang approximation
nb = max(1, int(round(1.0 / betascv)))
rate_b = nb * betarate
beta_D0 = np.diag([-rate_b] * nb)
for i in range(nb - 1):
beta_D0[i, i + 1] = rate_b
beta_D1 = np.zeros((nb, 1))
beta_D1[-1, 0] = rate_b
na = alpha_D0.shape[0]
nb = beta_D0.shape[0]
n = na + nb
# Build QBD blocks
# States: [setup phases (1..na), active + turn-off phases (na+1..n)]
F = np.zeros((n, n)) # Forward transitions (arrivals)
B = np.zeros((n, n)) # Backward transitions (service)
# Arrivals in setup phase
for i in range(na):
F[i, i] = lambda_val
# Arrivals in turn-off phase (go to active state)
for i in range(nb):
F[na + i, na] = lambda_val
F[na, na] = lambda_val
# Service completions (only from active state)
B[na, na] = mu
# Local transitions
L = np.zeros((n, n))
# Setup phase transitions
for i in range(na):
L[i, i] = alpha_D0[i, i] - lambda_val
if i < na - 1:
L[i, i + 1:na] = alpha_D0[i, i + 1:na]
else:
# Transition from last setup phase to active
L[na - 1, na] = -alpha_D0[na - 1, na - 1]
# Active state
L[na, na] = -mu - lambda_val
# Turn-off phase (only reachable from level 0)
for i in range(1, nb):
L[na + i, na + i] = -lambda_val
# Boundary block L0 for level 0
L0 = np.zeros((n, n))
# Setup phase at level 0
for i in range(na):
L0[i, i] = -lambda_val
# Turn-off phase at level 0
for i in range(nb):
L0[na + i, na + i] = beta_D0[i, i] - lambda_val if i < nb else -lambda_val
if i == nb - 1:
L0[na + i, 0] = -beta_D0[i, i] # Return to setup phase
elif i < nb - 1:
L0[na + i, na + i + 1] = -beta_D0[i, i]
# Compute R matrix using QBD_CR equivalent
R = qbd_R(B, L, F)
# Compute steady-state distribution using QBD_pi algorithm
# Follow MATLAB QBD_pi: convert to discrete time first
I = np.eye(n)
# Uniformization: find maximum exit rate from boundary block
lamb = max(-np.diag(L0))
if lamb <= 0:
lamb = 1.0
# Convert to discrete time stochastic matrices
B1_dt = L0 / lamb + I # Boundary local block
B0_dt = B / lamb # Backward transitions
# Compute stochastic matrix for level 0
stat_matrix = B1_dt + R @ B0_dt
# Find stationary distribution using stat() approach:
# Solve: K @ [A - I, e] = [0, ..., 0, 1]
# This is equivalent to: K @ (A - I) = 0 and K @ e = 1
e = np.ones((n, 1))
aug_matrix = np.hstack([stat_matrix - I, e])
y = np.zeros(n + 1)
y[-1] = 1.0
# Solve K @ aug_matrix = y using least squares (K = y @ pinv(aug_matrix))
try:
pi0 = linalg.lstsq(aug_matrix.T, y, cond=None)[0]
except:
pi0 = np.linalg.lstsq(aug_matrix.T, y, rcond=None)[0]
pi0 = np.abs(pi0) # Ensure non-negative
# Normalize using QBD normalization: pi @ (I-R)^{-1} @ 1 = 1
try:
temp = linalg.inv(I - R)
except:
temp = linalg.pinv(I - R)
norm_const = pi0 @ temp @ np.ones(n)
if norm_const > 0:
pi0 = pi0 / norm_const
# Build full probability vector pn following MATLAB QBD_pi
# Generate level probabilities until total mass approaches 1
max_num_comp = 500
pi_levels = [pi0]
sum_pi = np.sum(pi0)
numit = 1
while sum_pi < 1 - 1e-10 and numit < max_num_comp:
pi_next = pi_levels[-1] @ R
pi_levels.append(pi_next)
numit += 1
sum_pi += np.sum(pi_next)
# Concatenate all levels into a single vector (like MATLAB's reshape(pi', 1, []))
pn = np.concatenate(pi_levels)
# Mean queue length: the level probabilities are a flat vector of n phases
# per level, so level ni occupies pn[ni*n : (ni+1)*n], exactly n entries.
# Summing n+1 of them while advancing by n let each window reach into the
# next level, so one phase per level was counted twice under two different
# weights, and the terminating check then dropped the last level: the queue
# length came out high by up to 15%, worse the slower the setup, which is
# where the overlapped phases hold most mass. Level 0 is skipped, holding no
# jobs. Mirrors MATLAB qbd_setupdelayoff.m.
QN = 0.0
j = n
ni = 0
while j + n <= len(pn):
ni += 1
QN += ni * np.sum(pn[j:j + n])
j += n
return QN
__all__ = [
'QBDResult',
'qbd_R',
'qbd_R_logred',
'qbd_rg',
'qbd_blocks_mapmap1',
'qbd_bmapbmap1',
'qbd_mapmap1',
'qbd_raprap1',
'qbd_setupdelayoff',
]