Matrix-Analytic Methods

QBD processes and matrix-analytic solutions.

The mam module implements matrix-analytic methods for analyzing queues with structured Markov chains, including QBD, MAP, RAP, and G/M/1 and M/G/1 type processes.

Key function categories:

Matrix-Analytic Methods (MAM) for MAP/PH distributions.

Native Python implementations for analyzing Markovian Arrival Processes (MAPs), Phase-Type (PH) distributions, and related matrix-analytic methods.

Key algorithms:

map_piq: CTMC steady-state of MAP map_pie: Embedded DTMC steady-state map_lambda: Arrival rate computation map_mean, map_var, map_scv: Moment computations solver_mam_map_bmap_1: MAP/BMAP/1 queue solver using GI/M/1-type ETAQA solver_mam_bmap_map_1: BMAP/MAP/1 queue solver using M/G/1-type ETAQA

map_infgen(D0, D1)[source]

Compute the infinitesimal generator of a MAP.

The generator Q = D0 + D1 represents the underlying CTMC.

Parameters:

  • D0 (numpy.ndarray) – Hidden transition matrix (non-arrival transitions)

  • D1 (numpy.ndarray) – Visible transition matrix (arrival transitions)

Returns:

Infinitesimal generator matrix Q = D0 + D1

Return type:

numpy.ndarray

map_piq(D0, D1=None)[source]

Compute steady-state distribution of the underlying CTMC of a MAP.

Solves πQ = 0 where Q = D0 + D1 is the generator.

Parameters:

  • D0 (numpy.ndarray) – Hidden transition matrix, or stacked [D0, D1] if D1 is None

  • D1 (numpy.ndarray) – Visible transition matrix (optional)

Returns:

Steady-state probability vector π

Return type:

numpy.ndarray

map_pie(D0, D1=None)[source]

Compute equilibrium distribution of embedded DTMC.

The embedded DTMC has transition matrix P = (-D0)^{-1} * D1. Its steady-state is π_e = π * D1 / (π * D1 * e).

Parameters:

  • D0 (numpy.ndarray) – Hidden transition matrix, or stacked [D0, D1] if D1 is None

  • D1 (numpy.ndarray) – Visible transition matrix (optional)

Returns:

Equilibrium distribution of embedded DTMC

Return type:

numpy.ndarray

map_lambda(D0, D1=None)[source]

Compute the arrival rate (λ) of a MAP.

The arrival rate is λ = π * D1 * e where π is the steady-state and e is the column vector of ones.

Parameters:

  • D0 (numpy.ndarray) – Hidden transition matrix, or stacked [D0, D1] if D1 is None

  • D1 (numpy.ndarray) – Visible transition matrix (optional)

Returns:

Arrival rate λ

Return type:

float

map_mean(D0, D1=None)[source]

Compute mean inter-arrival time of a MAP.

The mean is 1/λ where λ is the arrival rate.

Parameters:

  • D0 (numpy.ndarray) – Hidden transition matrix, or stacked [D0, D1] if D1 is None

  • D1 (numpy.ndarray) – Visible transition matrix (optional)

Returns:

Mean inter-arrival time

Return type:

float

map_var(D0, D1=None)[source]

Compute variance of inter-arrival times of a MAP.

Var[X] = E[X²] - E[X]²

Uses map_moment for consistency with Kotlin implementation.

Parameters:

  • D0 (numpy.ndarray) – Hidden transition matrix, or stacked [D0, D1] if D1 is None

  • D1 (numpy.ndarray) – Visible transition matrix (optional)

Returns:

Variance of inter-arrival times

Return type:

float

map_scv(D0, D1=None)[source]

Compute squared coefficient of variation (SCV) of a MAP.

SCV = Var[X] / E[X]² = (E[X²] - E[X]²) / E[X]²

Parameters:

  • D0 (numpy.ndarray) – Hidden transition matrix, or stacked [D0, D1] if D1 is None

  • D1 (numpy.ndarray) – Visible transition matrix (optional)

Returns:

Squared coefficient of variation

Return type:

float

map_moment(D0, D1, k)[source]

Compute the k-th moment of inter-arrival time distribution.

E[X^k] = k! * π_e * (-D0)^{-k} * e

where π_e is the embedded DTMC steady-state.

Parameters:
Returns:

k-th raw moment

Return type:

float

map_scale(D0, D1, factor)[source]

Scale a MAP by a given factor.

Scaling changes the time scale: D0’ = factor * D0, D1’ = factor * D1.

Parameters:
Returns:

Tuple of scaled (D0’, D1’)

Return type:

Tuple[numpy.ndarray, numpy.ndarray]

map_normalize(D0, D1)[source]

Normalize a MAP to have unit mean inter-arrival time.

Parameters:
Returns:

Tuple of normalized (D0’, D1’)

Return type:

Tuple[numpy.ndarray, numpy.ndarray]

map_isfeasible(D0, D1, tolerance=1e-10)[source]

Check if (D0, D1) form a valid MAP.

A valid MAP requires: - D0 has non-positive diagonal and non-negative off-diagonal - D1 has non-negative elements - D0 + D1 is a valid generator (row sums = 0)

Parameters:
Returns:

True if valid MAP

Return type:

bool

exp_map(lambda_rate)[source]

Create a MAP representation of an exponential distribution.

Parameters:

lambda_rate (float) – Rate parameter (λ > 0)

Returns:

Tuple of (D0, D1) matrices representing Exp(λ)

Return type:

Tuple[numpy.ndarray, numpy.ndarray]

erlang_map(k, lambda_rate)[source]

Create a MAP representation of an Erlang-k distribution.

Parameters:

  • k (int) – Number of phases (k >= 1)

  • lambda_rate (float) – Overall rate parameter

Returns:

Tuple of (D0, D1) matrices representing Erlang(k, k*λ)

Return type:

Tuple[numpy.ndarray, numpy.ndarray]

hyperexp_map(probs, rates)[source]

Create a MAP representation of a hyperexponential distribution.

Parameters:
Returns:

Tuple of (D0, D1) matrices representing hyperexponential

Return type:

Tuple[numpy.ndarray, numpy.ndarray]

map_exponential(mean)[source]

Create a MAP representation of an exponential distribution (MATLAB-style).

This is a MATLAB-compatible wrapper that takes mean instead of rate.

Parameters:

mean (float) – Mean inter-arrival time (= 1/λ)

Returns:

Tuple of (D0, D1) matrices representing Exp(1/mean)

Return type:

Tuple[numpy.ndarray, numpy.ndarray]

Examples

>>> D0, D1 = map_exponential(2)  # Poisson process with rate λ=0.5
map_erlang(mean, k)[source]

Create a MAP representation of an Erlang-k distribution (MATLAB-style).

This is a MATLAB-compatible wrapper that takes mean as first argument.

Parameters:

  • mean (float) – Mean inter-arrival time

  • k (int) – Number of phases

Returns:

Tuple of (D0, D1) matrices representing Erlang-k with given mean

Return type:

Tuple[numpy.ndarray, numpy.ndarray]

Examples

>>> D0, D1 = map_erlang(2, 3)  # Erlang-3 with mean 2
map_hyperexp(probs, means)[source]

Create a MAP representation of a hyperexponential distribution (MATLAB-style).

This is a MATLAB-compatible wrapper that takes means instead of rates.

Parameters:

  • probs (numpy.ndarray) – Probability vector for choosing each phase

  • means (numpy.ndarray) – Mean values for each phase (1/rates)

Returns:

Tuple of (D0, D1) matrices representing hyperexponential

Return type:

Tuple[numpy.ndarray, numpy.ndarray]

map_gamma(mean, n)[source]

Create a MAP approximation of a Gamma distribution.

Uses Erlang-n with matching mean as an approximation.

Parameters:

  • mean (float) – Mean of the distribution

  • n (int) – Number of phases (shape parameter)

Returns:

Tuple of (D0, D1) matrices

Return type:

Tuple[numpy.ndarray, numpy.ndarray]

map_sumind(maps)[source]

Compute the sum of independent MAPs.

Creates a MAP representing the sum (concatenation) of independent random variables represented by the input MAPs.

Parameters:

maps (list) – List of MAPs, each as (D0, D1) tuple

Returns:

Tuple of (D0, D1) matrices representing the sum

Return type:

Tuple[numpy.ndarray, numpy.ndarray]

Examples

>>> # Sum of exponential and Erlang-2
>>> MAP1 = map_exponential(1.0)
>>> MAP2 = map_erlang(1.0, 2)
>>> D0, D1 = map_sumind([MAP1, MAP2])
map_cdf(D0, D1, points)[source]

Compute cumulative distribution function of inter-arrival times.

F(t) = 1 - π_e * exp(D0*t) * e

Parameters:
Returns:

CDF values at specified points

Return type:

numpy.ndarray

Examples

>>> map_cdf(D0, D1, 1.0)  # Returns P(T <= 1)
>>> map_cdf(D0, D1, [1.0, 5.0])  # Returns [P(T<=1), P(T<=5)]
map_pdf(D0, D1, points)[source]

Compute probability density function of inter-arrival times.

f(t) = π_e * exp(D0*t) * (-D0) * e

Parameters:
Returns:

PDF values at specified points

Return type:

numpy.ndarray

Examples

>>> map_pdf(D0, D1, [0.5, 1.0, 2.0])
map_sample(D0, D1, n_samples, rng=None)[source]

Generate random samples from a MAP distribution.

Parameters:
Returns:

Array of inter-arrival times

Return type:

numpy.ndarray

map_skew(D0, D1)[source]

Compute skewness of inter-arrival times.

Parameters:
Returns:

Skewness of inter-arrival times

Return type:

float

map_kurt(D0, D1)[source]

Compute kurtosis of inter-arrival times.

Parameters:
Returns:

Kurtosis of inter-arrival times

Return type:

float

map_acf(D0, D1, lags=1)[source]

Compute autocorrelation coefficients of inter-arrival times.

Parameters:
Returns:

Array of autocorrelation coefficients at specified lags

Return type:

numpy.ndarray

Examples

>>> map_acf(D0, D1)  # lag-1 autocorrelation
>>> map_acf(D0, D1, np.arange(1, 11))  # first 10 autocorrelations
map_acfc(D0, D1, kset, u)[source]

Compute autocorrelation of counting process at given lags.

Parameters:
Returns:

Autocorrelation coefficients at specified lags

Return type:

numpy.ndarray

map_idc(D0, D1)[source]

Compute the asymptotic index of dispersion.

I = SCV * (1 + 2 * sum_{k=1}^{inf} rho_k)

where SCV is the squared coefficient of variation and rho_k is the lag-k autocorrelation coefficient.

Parameters:
Returns:

Asymptotic index of dispersion

Return type:

float

map_count_mean(D0, D1, t)[source]

Compute mean of counting process at resolution t.

Parameters:
Returns:

Mean arrivals in (0, t]

Return type:

numpy.ndarray

map_count_var(D0, D1, t)[source]

Compute variance of counting process at resolution t.

Parameters:
Returns:

Variance of arrivals in (0, t]

Return type:

numpy.ndarray

Reference:

He and Neuts, “Markov chains with marked transitions”, 1998

map_varcount(D0, D1, tset)[source]

Compute variance of counting process (alternative implementation).

Parameters:
Returns:

Variance at each time point

Return type:

numpy.ndarray

map_count_moment(D0, D1, t, orders)[source]

Compute power moments of counts at resolution t.

Uses numerical differentiation of the moment generating function.

Parameters:
Returns:

Power moments of counts

Return type:

numpy.ndarray

map_mmpp2(mean, scv, skew=-1, acf1=-1)[source]

Fit an MMPP(2) as a MAP.

Parameters:

  • mean (float) – Mean inter-arrival time

  • scv (float) – Squared coefficient of variation

  • skew (float) – Skewness (-1 for automatic minimization)

  • acf1 (float) – Lag-1 autocorrelation (-1 for maximum feasible)

Returns:

Tuple of (D0, D1) matrices

Return type:

Tuple[numpy.ndarray, numpy.ndarray]

Examples

>>> D0, D1 = map_mmpp2(1, 2, -1, 0.2)  # Minimal skewness, ACF=0.2
>>> D0, D1 = map_mmpp2(1, 2, -1, -1)   # Minimal skewness, max ACF
map_gamma2(D0, D1)[source]

Compute the second largest eigenvalue of embedded DTMC.

This is the autocorrelation decay rate.

Parameters:
Returns:

Second largest eigenvalue (gamma_2)

Return type:

float

map_rand(k=2)[source]

Generate a random MAP of order k.

Parameters:

k (int) – Order of the MAP (default: 2)

Returns:

Tuple of (D0, D1) matrices

Return type:

Tuple[numpy.ndarray, numpy.ndarray]

map_randn(k, mu=(1.0, 1.0), sigma=(0.5, 0.5))[source]

Generate a random MAP with normally distributed elements.

Parameters:
Returns:

Tuple of (D0, D1) matrices

Return type:

Tuple[numpy.ndarray, numpy.ndarray]

map_renewal(D0, D1)[source]

Remove all correlations from a MAP, creating a renewal process.

Parameters:
Returns:

Renewal MAP with same CDF but no correlations

Return type:

Tuple[numpy.ndarray, numpy.ndarray]

map_embedded(D0, D1)[source]

Compute embedded discrete-time transition matrix.

P = (-D0)^{-1} * D1

Parameters:
Returns:

Transition matrix of embedded DTMC

Return type:

numpy.ndarray

map_sum(D0, D1, n)[source]

Create MAP for sum of n IID random variables.

Parameters:
Returns:

Tuple of (D0_new, D1_new) for the sum

Return type:

Tuple[numpy.ndarray, numpy.ndarray]

map_super(D0_a, D1_a, D0_b, D1_b)[source]

Create superposition of two MAPs.

Parameters:
Returns:

Superposed MAP (D0, D1)

Return type:

Tuple[numpy.ndarray, numpy.ndarray]

map_mixture(alpha, maps)[source]

Create probabilistic mixture of MAPs.

Parameters:

  • alpha (numpy.ndarray) – Probability vector for choosing each MAP

  • maps (list) – List of MAPs as (D0, D1) tuples

Returns:

Mixture MAP (D0, D1)

Return type:

Tuple[numpy.ndarray, numpy.ndarray]

map_max(D0_a, D1_a, D0_b, D1_b)[source]

Create MAP for max(X, Y) where X ~ MAP_A and Y ~ MAP_B.

Parameters:
Returns:

MAP for the maximum (D0, D1)

Return type:

Tuple[numpy.ndarray, numpy.ndarray]

map_timereverse(D0, D1)[source]

Compute time-reversed MAP.

Parameters:
Returns:

Time-reversed MAP (D0_r, D1_r)

Return type:

Tuple[numpy.ndarray, numpy.ndarray]

map_mark(D0, D1, prob)[source]

Mark arrivals from a MAP according to given probabilities.

Parameters:
Returns:

MMAP as list [D0, D1_class1, D1_class2, …]

Return type:

list

map_stochcomp(D0, D1, retain_idx)[source]

Stochastic complementation to reduce MAP order.

Parameters:
Returns:

Reduced MAP (D0_new, D1_new)

Return type:

Tuple[numpy.ndarray, numpy.ndarray]

map_kpc(maps)[source]

Kronecker product composition of MAPs.

Parameters:

maps (list) – List of MAPs as (D0, D1) tuples

Returns:

Composed MAP (D0, D1)

Return type:

Tuple[numpy.ndarray, numpy.ndarray]

map_bernstein(f, n=20)[source]

Convert distribution to MAP via Bernstein approximation.

Parameters:

  • f – PDF function handle

  • n (int) – Number of phases (default: 20)

Returns:

MAP representation (D0, D1)

Return type:

Tuple[numpy.ndarray, numpy.ndarray]

map_pntiter(D0, D1, na, t, M=None)[source]

Compute probability of na arrivals in interval [0, t].

Uses iterative bisection method (Neuts and Li).

Parameters:

  • D0 (numpy.ndarray) – Hidden transition matrix

  • D1 (numpy.ndarray) – Visible transition matrix

  • na (int) – Number of arrivals

  • t (float) – Time interval length

  • M (int | None) – Bisection depth (auto-computed if None)

Returns:

Matrix of probabilities

Return type:

numpy.ndarray

map_pntquad(D0, D1, na, t)[source]

Compute probability of na arrivals using ODE solver.

Parameters:
Returns:

Matrix P_n(t)

Return type:

numpy.ndarray

map2_fit(e1, e2, e3=-1.0, g2=0.0)[source]

Fit a MAP(2) distribution to moments and autocorrelation.

Based on: A. Heindl, G. Horvath, K. Gross “Explicit inverse characterization of acyclic MAPs of second order”

Parameters:

  • e1 (float) – First moment E[X]

  • e2 (float) – Second moment E[X^2]

  • e3 (float) – Third moment E[X^3], or -1 for automatic selection, -2 for minimum, -3 for maximum, or negative fraction for interpolation

  • g2 (float) – Autocorrelation decay rate (gamma_2)

Returns:

  • MAP: Tuple (D0, D1) if successful, None if failed

  • error_code: 0 for success, >0 for various errors

Return type:

Tuple of (MAP, error_code) where

Error codes:

0: Success 10: Mean out of bounds 20: Correlated exponential 30: h2 out of bounds 40: h3 out of bounds 51-54: g2 out of bounds

map_joint(D0, D1, a, i)[source]

Compute joint moments of a MAP.

E[(X_{a1})^{i1} * (X_{a1+a2})^{i2} * …]

Parameters:
Returns:

Joint moment

Return type:

float

map_issym(D0, D1=None)[source]

Check if MAP contains symbolic elements.

In Python, this always returns False as we use numpy arrays.

Parameters:
Returns:

False (Python arrays are always numeric)

Return type:

bool

map_feastol()[source]

Get tolerance for feasibility checks.

Returns:

Tolerance exponent (10^(-feastol))

Return type:

int

map_feasblock(E1, E2, E3, G2, opt='')[source]

Fit the most similar feasible MAP(2).

Parameters:

  • E1 (float) – First moment (mean)

  • E2 (float) – Second moment (or SCV if opt=’scv’)

  • E3 (float) – Third moment

  • G2 (float) – Autocorrelation decay rate

  • opt (str) – ‘scv’ if E2 is SCV instead of second moment

Returns:

Feasible MAP (D0, D1)

Return type:

Tuple[numpy.ndarray, numpy.ndarray]

map_block(E1, E2, E3, G2, opt='')[source]

Construct a MAP(2) from moments and autocorrelation.

Parameters:

  • E1 (float) – First moment

  • E2 (float) – Second moment (or SCV if opt=’scv’)

  • E3 (float) – Third moment

  • G2 (float) – Autocorrelation decay rate

  • opt (str) – ‘scv’ if E2 is SCV

Returns:

MAP (D0, D1)

Return type:

Tuple[numpy.ndarray, numpy.ndarray]

map_largemap()[source]

Get threshold for “large” MAP where exact computation is expensive.

Returns:

Return type:

Order threshold (default

mmap_infgen(D0, D_list)[source]

Compute the infinitesimal generator of an MMAP.

Parameters:

  • D0 (numpy.ndarray) – Hidden transition matrix (non-arrival transitions)

  • D_list (List[numpy.ndarray]) – List of arrival marking matrices [D1, D2, …, Dm]

Returns:

Infinitesimal generator Q = D0 + D1 + D2 + … + Dm

Return type:

numpy.ndarray

mmap_normalize(D0, D_list)[source]

Normalize an MMAP to ensure -diag(D0) = rowsums(D1 + D2 + … + Dm).

This ensures D0 represents true hidden transitions (negative diagonal).

Parameters:
Returns:

Tuple of (normalized_D0, normalized_D_list)

Return type:

Tuple[numpy.ndarray, List[numpy.ndarray]]

mmap_super_safe(mmap_list, maxorder=1000, method='default')[source]

Safely superpose multiple MMAPs into a single MMAP.

Combines multiple MMAPs while ensuring the order of the resulting MMAP does not exceed maxorder. If combining would exceed maxorder, alternative methods are used (fitting to exponential or simplified MAP).

Parameters:

  • mmap_list (List[Tuple[numpy.ndarray, List[numpy.ndarray]]]) – List of (D0, D_list) tuples, one per MMAP

  • maxorder (int) – Maximum allowed order (state space size) for result

  • method (str) – Combination method (“default” or “match”)

Returns:

Tuple of (D0_super, D_list_super) for the superposed MMAP

Return type:

Tuple[numpy.ndarray, List[numpy.ndarray]]

Algorithm:

  1. Sort MMAPs by squared coefficient of variation (SCV) of unmarked process

  2. Iteratively combine: start with simplest, add others in order

  3. If Kronecker product would exceed maxorder, use simpler approximations

mmap_mark(D0, D_list, prob)[source]

Reclassify arrivals in an MMAP according to a probability matrix.

Converts an MMAP with K classes to one with R classes based on a KxR probability matrix that describes how arrivals are reclassified.

Parameters:

  • D0 (numpy.ndarray) – Hidden transition matrix

  • D_list (List[numpy.ndarray]) – List of K arrival marking matrices [D1, D2, …, DK]

  • prob (numpy.ndarray) – K x R probability matrix where prob[k,r] is the probability that a class-k arrival in the original MMAP becomes a class-r arrival in the new MMAP

Returns:

Tuple of (D0_new, D_list_new) where D_list_new has R matrices

Return type:

Tuple[numpy.ndarray, List[numpy.ndarray]]

Example

If prob = [[0.8, 0.2], [0.3, 0.7]], a 2-class MMAP becomes a 2-class MMAP where 80% of class-1 arrivals are now class-1, 20% become class-2, etc.

mmap_scale(D0, D_list, M, max_iter=30)[source]

Scale the mean inter-arrival times of an MMAP.

Adjusts all matrices to achieve specified mean inter-arrival times.

Parameters:

  • D0 (numpy.ndarray) – Hidden transition matrix

  • D_list (List[numpy.ndarray]) – List of arrival marking matrices

  • M (float | numpy.ndarray) – Desired mean inter-arrival times. Can be: - Single float: uniform scaling of all classes - Array of K values: individual scaling per class

  • max_iter (int) – Maximum iterations for optimization (only for vector M)

Returns:

Tuple of (D0_scaled, D_list_scaled)

Return type:

Tuple[numpy.ndarray, List[numpy.ndarray]]

Algorithm:

For single M: Scale all matrices uniformly by ratio = M_old / M For vector M: Use iterative coordinate descent to find scaling factors

mmap_hide(D0, D_list, types)[source]

Hide (remove) specified arrival classes from an MMAP.

The hidden classes are set to zero matrices, effectively removing them from observation while maintaining the underlying stochastic process.

Parameters:
Returns:

Tuple of (D0_hidden, D_list_hidden) where specified classes are zero

Return type:

Tuple[numpy.ndarray, List[numpy.ndarray]]

Example

mmap_hide(D0, [D1, D2, D3], types=[1]) hides class 2 (index 1)

mmap_compress(D0, D_list, method='default', target_order=None)[source]

Compress an MMAP using various approximation methods.

Reduces the state space of an MMAP while preserving key statistical properties (moments, inter-arrival time distribution, etc.).

Parameters:

  • D0 (numpy.ndarray) – Hidden transition matrix

  • D_list (List[numpy.ndarray]) – List of arrival marking matrices

  • method (str) – Compression method. Options: - “default”, “mixture”: Mixture fitting (order 1) - “exponential”: Fit to single-state Poisson

  • target_order (int | None) – Target state space size (overrides method if specified)

Returns:

Tuple of (D0_compressed, D_list_compressed)

Return type:

Tuple[numpy.ndarray, List[numpy.ndarray]]

Notes

The “default” mixture method fits each class independently to MMPP(2) if possible, then combines them. This is fast and often effective.

More advanced methods (m3pp, mamap2) can be added as extensions.

mmap_exponential(lambda_rate, nclasses=None)[source]

Create an exponential MMAP (Poisson arrival process).

Parameters:

  • lambda_rate (float | numpy.ndarray) – Arrival rate (scalar) or rate per class (vector)

  • nclasses (int | None) – Number of classes (if lambda_rate is scalar)

Returns:

Tuple of (D0, D_list) representing exponential MMAP

Return type:

Tuple[numpy.ndarray, List[numpy.ndarray]]

mmap_issym(mmap)[source]

Check if an MMAP contains symbolic elements.

In Python/NumPy, we don’t have built-in symbolic support like MATLAB, so this always returns False for numpy arrays.

Parameters:

mmap (List[numpy.ndarray]) – List of matrices [D0, D1, D2, …, Dm]

Returns:

False (symbolic computation not supported in numpy)

Return type:

bool

mmap_isfeasible(mmap, tol=None)[source]

Check whether an MMAP is feasible up to the given tolerance.

Checks: - Elements are real (no imaginary parts) - D0 + D1 rows sum to zero (generator property) - Diagonal of D0 is negative - Off-diagonal of D0 is non-negative - All D1c matrices are non-negative - D1 = D11 + D12 + … + D1C

Parameters:

  • mmap (List[numpy.ndarray]) – List of matrices [D0, D1, D2, …, Dm]

  • tol (float) – Tolerance for feasibility checks (default: 1e-6)

Returns:

True if MMAP is feasible, False otherwise

Return type:

bool

mmap_lambda(mmap)[source]

Compute the arrival rate of each class in an MMAP.

Parameters:

mmap (List[numpy.ndarray]) – List of matrices [D0, D1, D2, …, Dm]

Returns:

Array of arrival rates, one per class

Return type:

numpy.ndarray

mmap_count_lambda(mmap)[source]

Compute the arrival rate of the counting process for each class.

The rate for class k is λ_k = π * D_k * e where π is the steady-state of the underlying CTMC.

Parameters:

mmap (List[numpy.ndarray]) – List of matrices [D0, D1, D2, …, Dm]

Returns:

Array of arrival rates, one per class

Return type:

numpy.ndarray

mmap_pie(mmap)[source]

Compute the stationary probability of the DTMC embedded at restart instants after an arrival of each class.

Parameters:

mmap (List[numpy.ndarray]) – List of matrices [D0, D1, D2, …, Dm]

Returns:

Matrix of shape (m, n) where row c is the steady-state for class c

Return type:

numpy.ndarray

mmap_pc(mmap)[source]

Compute the arrival probabilities of each class for the given MMAP.

Parameters:

mmap (List[numpy.ndarray]) – List of matrices [D0, D1, D2, …, Dm]

Returns:

Array of arrival probabilities, one per class

Return type:

numpy.ndarray

mmap_embedded(mmap)[source]

Compute the embedded DTMC transition matrices for each class.

For class k, the embedded matrix is Pc_k = (-D0)^{-1} * D_{2+k}

Parameters:

mmap (List[numpy.ndarray]) – List of matrices [D0, D1, D2, …, Dm]

Returns:

List of embedded DTMC transition matrices, one per class

Return type:

List[numpy.ndarray]

mmap_sample(mmap, n_samples, pi=None, seed=None)[source]

Generate samples from a Marked MAP.

Parameters:

  • mmap (List[numpy.ndarray]) – List of matrices [D0, D1, D2, …, Dm]

  • n_samples (int) – Number of random samples to collect

  • pi (numpy.ndarray) – Initial probability (if None, uses steady-state)

  • seed (int) – Random seed for reproducibility

Returns:

  • T: Array of inter-arrival times

  • A: Array of class labels (1-indexed for MATLAB compatibility)

Return type:

Tuple of (T, A) where

mmap_rand(order, classes)[source]

Generate a random MMAP with given order and number of classes.

Parameters:

  • order (int) – Number of states

  • classes (int) – Number of arrival classes

Returns:

List of matrices [D0, D1, D21, D22, …, D2K]

Return type:

List[numpy.ndarray]

mmap_timereverse(mmap)[source]

Compute the time-reversed MMAP.

Parameters:

mmap (List[numpy.ndarray]) – List of matrices [D0, D1, D2, …, Dm]

Returns:

Time-reversed MMAP

Return type:

List[numpy.ndarray]

mmap_sum(mmap, n)[source]

Create an MMAP representing the sum of n identical MMAPs (n-fold sum).

This creates an MMAP whose inter-arrival times are distributed as the sum of n inter-arrival times from the original MMAP.

Parameters:
Returns:

Summed MMAP

Return type:

List[numpy.ndarray]

mmap_super(mmap_a, mmap_b=None, opt='default')[source]

Superpose two or more MMAPs.

Parameters:

  • mmap_a (List[numpy.ndarray]) – First MMAP (or list of MMAPs if mmap_b is None)

  • mmap_b (List[numpy.ndarray]) – Second MMAP (optional)

  • opt (str) – “default” - each class in MMAPa and MMAPb is distinct in result “match” - class c in both maps is mapped to class c in result

Returns:

Superposed MMAP

Return type:

List[numpy.ndarray]

mmap_mixture(alpha, maps)[source]

Create a probabilistic mixture of MAPs.

Each MAP in the list is selected with probability alpha[i].

Parameters:

  • alpha (numpy.ndarray) – Array of mixing probabilities

  • maps (List) – List of MAPs (each MAP is [D0, D1])

Returns:

Mixed MMAP with len(maps) classes

Return type:

List[numpy.ndarray]

mmap_max(mmap_a, mmap_b, k)[source]

Create MMAP for the maximum of arrivals from two synchronized MMAPs.

This models a synchronization queue where arrivals from both sources must be paired before release.

Parameters:
Returns:

MMAP for the synchronized (maximum) process

Return type:

List[numpy.ndarray]

mmap_maps(mmap)[source]

Extract K MAPs, one for each class of the MMAP[K] process.

For class k, the MAP is {D0 + D1 - D_{2+k}, D_{2+k}}

Parameters:

mmap (List[numpy.ndarray]) – List of matrices [D0, D1, D21, D22, …, D2K]

Returns:

List of MAPs, one per class

Return type:

List[Tuple[numpy.ndarray, numpy.ndarray]]

mmap_count_mean(mmap, t)[source]

Compute the mean of the counting process at resolution t.

Parameters:

  • mmap (List[numpy.ndarray]) – List of matrices [D0, D1, D2, …, Dm]

  • t (float) – Time period for counting

Returns:

Array of mean counts, one per class

Return type:

numpy.ndarray

mmap_count_var(mmap, t)[source]

Compute the variance of the counting process at resolution t.

Parameters:

  • mmap (List[numpy.ndarray]) – List of matrices [D0, D1, D2, …, Dm]

  • t (float) – Time period for counting

Returns:

Array of variances, one per class

Return type:

numpy.ndarray

mmap_count_idc(mmap, t)[source]

Compute the per-class Index of Dispersion of Counts at resolution t.

IDC = Var[N(t)] / E[N(t)]

Parameters:
Returns:

Array of IDC values, one per class

Return type:

numpy.ndarray

mmap_count_mcov(mmap, t)[source]

Compute the count covariance between each pair of classes at time scale t.

Parameters:
Returns:

m x m covariance matrix

Return type:

numpy.ndarray

mmap_idc(mmap)[source]

Compute the asymptotic Index of Dispersion of Counts for each class.

This is the limit of IDC(t) as t -> infinity.

Parameters:

mmap (List[numpy.ndarray]) – List of matrices [D0, D1, D2, …, Dm]

Returns:

Array of asymptotic IDC values, one per class

Return type:

numpy.ndarray

mmap_sigma(mmap)[source]

Compute one-step class transition probabilities.

p_{i,j} = P(C_k = j | C_{k-1} = i)

Parameters:

mmap (List[numpy.ndarray]) – List of matrices [D0, D1, D2, …, Dm]

Returns:

C x C matrix of transition probabilities

Return type:

numpy.ndarray

mmap_sigma2(mmap)[source]

Compute two-step class transition probabilities.

p_{i,j,h} = P(C_k = h | C_{k-1} = j, C_{k-2} = i)

Parameters:

mmap (List[numpy.ndarray]) – List of matrices [D0, D1, D2, …, Dm]

Returns:

C x C x C 3D array of transition probabilities

Return type:

numpy.ndarray

mmap_forward_moment(mmap, orders, norm=True)[source]

Compute the theoretical forward moments of an MMAP.

Forward moments are E[X^k | previous arrival was class c].

Parameters:

  • mmap (List[numpy.ndarray]) – List of matrices [D0, D1, D2, …, Dm]

  • orders (numpy.ndarray) – Array of moment orders to compute

  • norm (bool) – If True, normalize by class probability

Returns:

Matrix of shape (C, len(orders)) where element (c, k) is the order-k forward moment for class c

Return type:

numpy.ndarray

mmap_backward_moment(mmap, orders, norm=True)[source]

Compute the theoretical backward moments of an MMAP.

Backward moments are E[X^k | next arrival will be class c].

Parameters:

  • mmap (List[numpy.ndarray]) – List of matrices [D0, D1, D2, …, Dm]

  • orders (numpy.ndarray) – Array of moment orders to compute

  • norm (bool) – If True, normalize by class probability

Returns:

Matrix of shape (C, len(orders)) where element (c, k) is the order-k backward moment for class c

Return type:

numpy.ndarray

mmap_cross_moment(mmap, k)[source]

Compute the k-th order moment of inter-arrival times between class pairs.

Computes E[X^k | previous = class i, next = class j] for all i, j.

Parameters:

  • mmap (List[numpy.ndarray]) – List of matrices [D0, D1, D2, …, Dm]

  • k (int) – Order of the moment

Returns:

C x C matrix where element (i, j) is the k-th moment for transitions i->j

Return type:

numpy.ndarray

mmap_modulate(P, HT, MMAP)[source]

Modulate MMAPs in continuous time according to holding time distributions.

Creates an MMAP that switches between different MMAPs according to a continuous-time Markov modulated process.

Parameters:

  • P (numpy.ndarray) – J x J transition probability matrix between states

  • HT (List) – List of phase-type distributions for holding time in each state

  • MMAP (List) – List of MMAPs, one for each state

Returns:

Modulated MMAP

Return type:

List[numpy.ndarray]

ldqbd(Q0, Q1, Q2, options=None)[source]

Solve a level-dependent QBD process.

Parameters:

  • Q0 (List[numpy.ndarray]) – List of upward transition matrices [Q0^(0), Q0^(1), …, Q0^(N-1)] Q0^(n) has shape (states_n, states_{n+1})

  • Q1 (List[numpy.ndarray]) – List of local transition matrices [Q1^(0), Q1^(1), …, Q1^(N)] Q1^(n) has shape (states_n, states_n)

  • Q2 (List[numpy.ndarray]) – List of downward transition matrices [Q2^(1), Q2^(2), …, Q2^(N)] Q2^(n) has shape (states_n, states_{n-1})

  • options (LdqbdOptions | None) – LdqbdOptions instance (uses defaults if None)

Returns:

LdqbdResult containing rate matrices R and stationary distribution pi

Return type:

LdqbdResult

Algorithm:

  1. Compute rate matrices backward using continued fraction recursion

  2. Compute stationary distribution forward using R matrices

class LdqbdResult(R, pi)[source]

Bases: object

Result of LDQBD solver.

R

List of rate matrices R^(1), R^(2), …, R^(N)

Type:

List[numpy.ndarray]

pi

Stationary distribution vector [pi_0, pi_1, …, pi_N]

Type:

numpy.ndarray

R: List[numpy.ndarray]
pi: numpy.ndarray
class LdqbdOptions(epsilon=1e-10, max_iter=1000, verbose=False)[source]

Bases: object

Options for LDQBD solver.

epsilon

Convergence tolerance (default: 1e-10)

Type:

float

max_iter

Maximum number of iterations (default: 1000)

Type:

int

verbose

Print debug information (default: False)

Type:

bool

epsilon: float = 1e-10
max_iter: int = 1000
verbose: bool = False
class MAPBMAP1Result(mean_queue_length, utilization, mean_response_time, throughput, pi, R, mean_batch_size)[source]

Bases: object

Result of MAP/BMAP/1 queue analysis.

mean_queue_length: float

Mean queue length E[N]

utilization: float

Server utilization rho

mean_response_time: float

Mean response time E[R]

throughput: float

Throughput (arrival rate)

pi: numpy.ndarray

Aggregated stationary probabilities [pi0, pi1, piStar]

R: numpy.ndarray

R matrix

mean_batch_size: float

Mean batch size of service

solver_mam_map_bmap_1(C0, C1, D)[source]

Solve a MAP/BMAP/1 queue using GI/M/1-type matrix-analytic methods.

The MAP is specified by matrices (C0, C1) where: - C0: transitions without arrivals - C1: transitions triggering arrivals

The BMAP for service is specified by matrices {D0, D1, D2, …, DK} where: - D0: transitions without service completions - Dk: transitions triggering batch service of k customers (k >= 1)

The GI/M/1-type structure for MAP/BMAP/1 is: ```

B1 A0 0 0 … B2 A1 A0 0 …

Q = B3 A2 A1 A0 …

…

```

Where:

A0 = C1 otimes I_ms (MAP arrival, level +1) A1 = C0 otimes I_ms + I_ma otimes D0 (phase changes, level 0) A_{k+1} = I_ma otimes D_k (batch size k service, level -k)

Parameters:
Returns:

MAPBMAP1Result with performance metrics

Return type:

MAPBMAP1Result

class BMAPMAP1Result(mean_queue_length, utilization, mean_response_time, throughput, pi, G, mean_batch_size)[source]

Bases: object

Result of BMAP/MAP/1 queue analysis.

mean_queue_length: float

Mean queue length E[N]

utilization: float

Server utilization rho

mean_response_time: float

Mean response time E[R]

throughput: float

Throughput (total customer arrival rate)

pi: numpy.ndarray

Aggregated stationary probabilities [pi0, pi1, piStar]

G: numpy.ndarray

G matrix

mean_batch_size: float

Mean batch size of arrivals

solver_mam_bmap_map_1(D, S0, S1)[source]

Solve a BMAP/MAP/1 queue using M/G/1-type matrix-analytic methods.

The BMAP is specified by matrices {D0, D1, D2, …, DK} where: - D0: transitions without arrivals - Dk: transitions triggering batch size k arrivals (k >= 1)

The MAP for service is specified by matrices (S0, S1) where: - S0: transitions without service completions - S1: transitions triggering service completions

The M/G/1-type structure for BMAP/MAP/1 is: ```

A0 = I_ma otimes S1 (service completion, level -1) A1 = D0 otimes I_ms + I_ma otimes S0 (phase changes, level 0) A_{k+1} = D_k otimes I_ms (batch arrival size k, level +k)

```

Parameters:

  • D (List[numpy.ndarray]) – BMAP matrices as list [D0, D1, D2, …, DK]

  • S0 (numpy.ndarray) – MAP service matrix for transitions without service completions

  • S1 (numpy.ndarray) – MAP service matrix for service completions

Returns:

BMAPMAP1Result with performance metrics

Return type:

BMAPMAP1Result

class QBDResult(R, G=None, U=None, eta=None)[source]

Bases: object

Result of QBD analysis.

G: numpy.ndarray | None = None
U: numpy.ndarray | None = None
eta: float | None = None
R: numpy.ndarray
qbd_R(B, L, F, iter_max=100000, tol=1e-12)[source]

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.

Parameters:

  • B (numpy.ndarray) – Backward transition block A_{-1}

  • L (numpy.ndarray) – Local transition block A_0

  • F (numpy.ndarray) – Forward transition block A_1

  • iter_max (int) – Maximum iterations (default: 100000)

  • tol (float) – Convergence tolerance (default: 1e-12)

Returns:

Rate matrix R

Return type:

numpy.ndarray

References

Original MATLAB: matlab/src/api/mam/qbd_R.m

qbd_R_logred(B, L, F, iter_max=1000, tol=1e-14)[source]

Compute QBD rate matrix R using logarithmic reduction.

Uses the logarithmic reduction algorithm which has quadratic convergence compared to linear convergence of successive substitutions.

Parameters:

  • B (numpy.ndarray) – Backward transition block A_{-1}

  • L (numpy.ndarray) – Local transition block A_0

  • F (numpy.ndarray) – Forward transition block A_1

  • iter_max (int) – Maximum iterations (default: 1000)

  • tol (float) – Convergence tolerance (default: 1e-14)

Returns:

Rate matrix R

Return type:

numpy.ndarray

References

Original MATLAB: matlab/src/api/mam/qbd_R_logred.m Latouche & Ramaswami, Ch. 8

qbd_rg(B, L, F, method='logred', iter_max=1000, tol=1e-14)[source]

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

Parameters:

  • B (numpy.ndarray) – Backward transition block A_{-1}

  • L (numpy.ndarray) – Local transition block A_0

  • F (numpy.ndarray) – Forward transition block A_1

  • method (str) – ‘logred’ or ‘successive’ (default: ‘logred’)

  • iter_max (int) – Maximum iterations

  • tol (float) – Convergence tolerance

Returns:

QBDResult with R, G, U, and eta (caudal characteristic)

Return type:

QBDResult

References

Original MATLAB: matlab/src/api/mam/qbd_rg.m

qbd_blocks_mapmap1(D0_arr, D1_arr, D0_srv, D1_srv)[source]

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.

Parameters:
Returns:

Tuple of (B, L, F) QBD blocks

Return type:

Tuple[numpy.ndarray, numpy.ndarray, numpy.ndarray]

References

Original MATLAB: matlab/src/api/mam/qbd_mapmap1.m

qbd_bmapbmap1(MAPa, pbatcha, MAPs)[source]

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.

Parameters:
Returns:

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

Return type:

Tuple of (A0, A_1, A1_list, B0, B1_list) where

References

Original MATLAB: matlab/src/api/mam/qbd_bmapbmap1.m

qbd_mapmap1(MAPa, MAPs, util=None)[source]

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.

Parameters:
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

Return type:

Tuple[float, float, float, numpy.ndarray, numpy.ndarray, float | None, numpy.ndarray, numpy.ndarray, numpy.ndarray, numpy.ndarray, numpy.ndarray, Tuple[numpy.ndarray, numpy.ndarray]]

References

Original MATLAB: matlab/src/api/mam/qbd_mapmap1.m

qbd_raprap1(RAPa, RAPs, util=None)[source]

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.

Note: This is a simplified implementation. For full RAP analysis, specialized algorithms like those in SMCSolver are needed.

Parameters:
Returns:

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

Return type:

Tuple of (XN, QN, UN, pqueue, R, eta, G, B, L, F) where

References

Original MATLAB: matlab/src/api/mam/qbd_raprap1.m

qbd_setupdelayoff(lambda_val, mu, alpharate, alphascv, betarate, betascv)[source]

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

Parameters:

  • lambda_val (float) – Arrival rate

  • mu (float) – Service rate

  • alpharate (float) – Rate of setup delay phase

  • alphascv (float) – Squared coefficient of variation for setup delay

  • betarate (float) – Rate of turn-off phase

  • betascv (float) – Squared coefficient of variation for turn-off period

Returns:

Average queue length QN

Return type:

float

References

Original MATLAB: matlab/src/api/mam/qbd_setupdelayoff.m

map_ccdf_derivative(MAP, i)[source]

Compute derivative at zero of a MAP’s complementary CDF.

Calculates the i-th derivative at t=0 of the complementary cumulative distribution function (CCDF) of the inter-arrival time distribution.

Formula: ν_i = π_e * D_0^i * e

where π_e is the embedded stationary vector and e is the column vector of ones.

Parameters:

  • MAP (List[numpy.ndarray]) – Markovian Arrival Process as [D0, D1]

  • i (int) – Order of the derivative

Returns:

Value of the i-th derivative at zero

Return type:

float

References

Original MATLAB: matlab/src/api/mam/map_ccdf_derivative.m A. Horvath et al., “A Joint Moments Based Analysis of Networks of MAP/MAP/1 Queues”

map_jointpdf_derivative(MAP, iset)[source]

Compute partial derivative at zero of a MAP’s joint PDF.

Calculates the partial derivative at t=0 of the joint probability density function of consecutive inter-arrival times.

For index set {i_1, i_2, …, i_k}: γ = π_e * D_0^{i_1} * D_1 * D_0^{i_2} * D_1 * … * D_0^{i_k} * D_1 * e

Parameters:
Returns:

Value of the partial derivative at zero

Return type:

float

References

Original MATLAB: matlab/src/api/mam/map_jointpdf_derivative.m A. Horvath et al., “A Joint Moments Based Analysis of Networks of MAP/MAP/1 Queues”

map_factorial_moment(MAP, k)[source]

Compute the k-th factorial moment of a MAP.

The k-th factorial moment is computed using derivatives of the inter-arrival time distribution.

Parameters:

  • MAP (List[numpy.ndarray]) – Markovian Arrival Process as [D0, D1]

  • k (int) – Order of the factorial moment

Returns:

k-th factorial moment

Return type:

float

References

Based on MAP moment formulas from matrix-analytic methods

map_joint_moment(MAP, k, l)[source]

Compute the (k,l)-th joint moment of consecutive inter-arrival times.

E[X_n^k * X_{n+1}^l] for a MAP with inter-arrival times X_n.

Parameters:

  • MAP (List[numpy.ndarray]) – Markovian Arrival Process as [D0, D1]

  • k (int) – Power for first inter-arrival time

  • l (int) – Power for second inter-arrival time

Returns:

Joint moment E[X_n^k * X_{n+1}^l]

Return type:

float

References

Based on joint moment formulas for MAPs

map_m1ps_cdf_respt(C, D, mu, x, epsilon=1e-11, epsilon_prime=1e-10)[source]

Compute complementary sojourn time CDF for MAP/M/1-PS queue.

Based on:

Masuyama, H., & Takine, T. (2003). Sojourn time distribution in a MAP/M/1 processor-sharing queue.

Parameters:

  • C (numpy.ndarray | list) – MAP C matrix (transitions without arrivals).

  • D (numpy.ndarray | list) – MAP D matrix (transitions with arrivals).

  • mu (float) – Service rate (must be > 0).

  • x (numpy.ndarray | list) – Time points for CDF evaluation.

  • epsilon (float) – Queue length truncation parameter.

  • epsilon_prime (float) – Uniformization truncation parameter.

Returns:

Complementary CDF values W_bar(x) = Pr[W > x] at each point in x.

Return type:

numpy.ndarray

Example

>>> # Exponential arrivals (Poisson process with rate lambda=0.5)
>>> C = np.array([[-0.5]])
>>> D = np.array([[0.5]])
>>> mu = 1.0
>>> x = np.array([0.0, 0.5, 1.0, 2.0, 5.0])
>>> cdf = map_m1ps_cdf_respt(C, D, mu, x)
map_compute_R(C, D, mu)[source]

Compute rate matrix R for MAP/M/1 queue.

Computes the minimal nonnegative solution of the matrix equation:

D + R(C - mu*I) + mu*R^2 = 0

This matrix is used in the analysis of MAP/M/1 queues based on quasi-birth-death processes.

Parameters:

  • C (numpy.ndarray | list) – M x M matrix governing MAP transitions without arrivals

  • D (numpy.ndarray | list) – M x M matrix governing MAP transitions with arrivals

  • mu (float) – Service rate (scalar)

Returns:

M x M rate matrix (minimal nonnegative solution)

Return type:

numpy.ndarray

References

Original MATLAB: matlab/src/api/map/map_compute_R.m Masuyama, H., & Takine, T. (2003). Sojourn time distribution in a MAP/M/1 processor-sharing queue.

map_m1ps_h_recursive(C, D, mu, N, K)[source]

Recursive computation of h_{n,k} coefficients for MAP/M/1-PS.

The h_{n,k} vectors satisfy the recursion (Theorem 1 in paper):

h_{n,0} = e (vector of ones), for n = 0, 1, … h_{n,k+1} = 1/(theta+mu) * [n*mu/(n+1) * h_{n-1,k} + (theta*I + C) * h_{n,k}

  • D * h_{n+1,k}]

where h_{-1,k} = 0 for all k

Parameters:

  • C (numpy.ndarray | list) – M x M matrix governing MAP transitions without arrivals

  • D (numpy.ndarray | list) – M x M matrix governing MAP transitions with arrivals

  • mu (float) – Service rate (scalar)

  • N (int) – Maximum value of n to compute (determines rows)

  • K (int) – Maximum value of k to compute (determines columns)

Returns:

2D list h[n][k] of (M,1) arrays containing h_{n,k}

Return type:

list

References

Original MATLAB: matlab/src/api/map/map_m1ps_h_recursive.m Masuyama, H., & Takine, T. (2003).

map_m1ps_sojourn(C, D, mu, x, epsilon=1e-11, epsilon_prime=1e-10)[source]

Compute sojourn time distribution in MAP/M/1-PS queue.

Alias for map_m1ps_cdf_respt.

Parameters:

  • C (numpy.ndarray | list) – MAP C matrix (transitions without arrivals)

  • D (numpy.ndarray | list) – MAP D matrix (transitions with arrivals)

  • mu (float) – Service rate (must be > 0)

  • x (numpy.ndarray | list) – Time points for CDF evaluation

  • epsilon (float) – Queue length truncation parameter

  • epsilon_prime (float) – Uniformization truncation parameter

Returns:

Complementary CDF values W_bar(x) = Pr[W > x]

Return type:

numpy.ndarray

References

Original MATLAB: matlab/src/api/map/map_m1ps_sojourn.m

mmdp_isfeasible(Q, R, tol=1e-10)[source]

Check if (Q, R) defines a valid MMDP.

Requirements: - Q must be a valid generator (square, row sums = 0, proper signs) - R must be diagonal with non-negative entries - Q and R must have compatible dimensions

Parameters:

  • Q (numpy.ndarray) – Generator matrix (n x n)

  • R (numpy.ndarray) – Rate matrix (n x n diagonal) or rate vector (n,)

  • tol (float) – Numerical tolerance for validation

Returns:

True if (Q, R) defines a valid MMDP, False otherwise

Return type:

bool

Examples

>>> Q = np.array([[-0.5, 0.5], [0.3, -0.3]])
>>> R = np.array([[2.0, 0], [0, 5.0]])
>>> mmdp_isfeasible(Q, R)
True