api.mam

ldqbd_mphc(D0, D1, alpha, c, arrRate, sf)

LDQBD_MPHC Level-dependent QBD blocks of an M/PH/c queue, exact in the phases.

[Q0, Q1, Q2] = LDQBD_MPHC(D0, D1, ALPHA, C, ARRRATE) returns the block-tridiagonal generator of a c-server FCFS queue whose service is the Markovian process (D0, D1) restarted from ALPHA, with a level-dependent arrival rate. Level n is the number of jobs at the station; the coordinate INSIDE a level is the MULTISET of the phases the min(n,c) busy servers sit in, so the chain is exact for phase-type service at any number of servers.

[Q0, Q1, Q2] = LDQBD_MPHC(…, SF) additionally scales the station’s total service rate at level n by SF(n) (load dependence). Each busy server then runs at SF(n)/min(n,c) of its nominal speed, so the aggregate capacity is SF(n) times nominal and SF(n) = min(n,c) reproduces the unscaled queue exactly. Pass [] for no scaling.

@par Parameters: <table> <tr><th>Name<th>Description <tr><td>D0<td>service sub-generator (p x p): phase changes without completion <tr><td>D1<td>service completion block (p x p); D1 = (-D0*1)*ALPHA for a PH <tr><td>ALPHA<td>1 x p vector a server starts each new job in <tr><td>C<td>number of identical servers (>= 1; capped at the top level) <tr><td>ARRRATE<td>1 x (Nlev+1); ARRRATE(n+1) is the arrival rate out of level n <tr><td>SF<td>optional 1 x Nlev total-service-rate multiplier per level </table>

@par Returns: <table> <tr><th>Name<th>Description <tr><td>Q0<td>{Nlev x 1}; Q0{n+1} is the upward block, level n -> n+1 <tr><td>Q1<td>{(Nlev+1) x 1}; Q1{n+1} is the local block of level n <tr><td>Q2<td>{Nlev x 1}; Q2{n} is the downward block, level n -> n-1 </table>

WHY THE MULTISET. The collapsed alternative – one PH process run at min(n,c) times its speed – gets the aggregate rate right but forgets which phase each busy server is in, which is not a detail: it makes the c servers behave like one fast server whose remaining work is a single phase-type variable. Tracking counts rather than an ordered tuple costs NCHOOSEK(min(n,c)+p-1, p-1) states per level instead of p^min(n,c), because identical servers are exchangeable.

Level sizes therefore GROW over the boundary levels 0..c and repeat above them: level 0 is the single empty configuration, and Q0{1}, Q2{1} are the rectangular blocks that join a level to a differently sized neighbour. LDQBD, LDQBD_R and LDQBD_PI all accept that heterogeneity.

References

S. Asmussen and J.R. Moller, “Calculation of the steady state waiting time distribution in GI/PH/c and MAP/PH/c queues”, Queueing Systems 37(1):9-29, 2001. M. F. Neuts, “Matrix-geometric solutions in stochastic models”, Johns Hopkins University Press, 1981.

ph_multisets(p, k)

PH_MULTISETS Configurations of k identical servers over p service phases.

M = PH_MULTISETS(P, K) returns the compositions of K into P nonnegative parts, one per row: M(r,i) is the number of the K busy servers sitting in phase i. There are NCHOOSEK(K+P-1, P-1) rows, the multiset count of Asmussen and Moller (2001): identical servers are exchangeable, so only the phase COUNTS carry information and the ordered space of size P^K collapses onto this one.

The order is fixed and shared by every caller, so a configuration index means the same thing in each of them: the first part descends. K = 1 therefore yields the identity rows e_1 … e_p in phase order, which is what makes the c = 1 case of LDQBD_MPHC coincide with the plain phase indexing.

Examples

ph_multisets(2, 0) -> [0 0] ph_multisets(3, 1) -> [1 0 0; 0 1 0; 0 0 1] ph_multisets(2, 2) -> [2 0; 1 1; 0 2]

References

S. Asmussen and J.R. Moller, “Calculation of the steady state waiting time distribution in GI/PH/c and MAP/PH/c queues”, Queueing Systems 37(1):9-29, 2001.

map_m1ps_sojourn(C, D, mu, x, varargin)

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

W_bar = MAP_M1PS_SOJOURN(C, D, mu, x) computes the complementary distribution function of the sojourn time in a MAP/M/1 processor-sharing queue at the points specified in x.

[W_bar, W_bar_n] = MAP_M1PS_SOJOURN(C, D, mu, x) also returns the conditional complementary distributions W_bar_n{i} for customers finding i-1 customers in the system on arrival.

[…] = MAP_M1PS_SOJOURN(C, D, mu, x, ‘Param’, Value) specifies optional :Parameters: * ‘Epsilon’ - Truncation parameter for queue length (default – 1e-11)

‘EpsilonPrime’ - Truncation parameter for uniformization (default – 1e-10)

‘Verbose’ - Display computation progress (default – false)

Input:: C - M x M matrix governing MAP transitions without arrivals D - M x M matrix governing MAP transitions with arrivals mu - Service rate (scalar, mu > 0) x - Vector of time points at which to evaluate W_bar(x) = Pr[W > x]
Output:: W_bar - Vector of same size as x, containing Pr[W > x] W_bar_n - Cell array of conditional distributions (optional)

The processor-sharing (PS) discipline shares the server equally among all customers. When n customers are present, each receives service at rate 1/n.

The algorithm implements Theorem 1 from:: Masuyama, H., & Takine, T. (2003). Sojourn time distribution in a MAP/M/1 processor-sharing queue. Operations Research Letters, 31(6), 406-412.

Example

% M/M/1-PS queue with lambda=0.8, mu=1 lambda = 0.8; mu = 1; C = -lambda; D = lambda; % Poisson arrivals x = linspace(0, 10, 100); W_bar = map_m1ps_sojourn(C, D, mu, x); plot(x, W_bar); xlabel(‘x’); ylabel(‘Pr[W > x]’); title(‘Sojourn time distribution for M/M/1-PS’);

See also: MAP_M1PS_H_RECURSIVE, MAP_COMPUTE_R

ldqbd_pi(R, Q0, Q1, Q2, varargin)

Computes the stationary distribution for a level-dependent QBD

This function implements Algorithm 3 from the Phung-Duc et al. paper for computing the stationary distribution of level-dependent QBD processes.

The algorithm uses: 1. Boundary equation: pi_0 * (Q1^(0) + R^(1)*Q2^(1)) = 0 2. Forward recursion: pi_n = pi_{n-1} * R^(n) 3. Normalize: sum(pi) = 1

Parameters:

R – Cell array {R^(1), R^(2), …, R^(N)} - rate matrices from ldqbd_R
Q0 – Cell array {Q0^(0), Q0^(1), …, Q0^(N-1)} - upward transitions
Q1 – Cell array {Q1^(0), Q1^(1), …, Q1^(N)} - local transitions
Q2 – Cell array {Q2^(1), Q2^(2), …, Q2^(N)} - downward transitions
options – Optional struct with fields: verbose (default false)

Returns:

pi – Row vector [pi_0, pi_1, …, pi_N] - stationary distribution pi_cell: (Optional) Cell array with full phase-level distributions

Examples

pi = ldqbd_pi(R, Q0, Q1, Q2)
pi = ldqbd_pi(R, Q0, Q1, Q2, options)

References

T. Phung-Duc, H. Masuyama, S. Kasahara, Y. Takahashi, “A Simple Algorithm for the Rate Matrices of Level-Dependent QBD Processes”, QTNA 2010.

map2_fit_idc(e1, e2, e3, I)

Fits a second-order MAP matching the first three moments and the asymptotic index of dispersion

A MAP(2) has a geometrically decaying autocorrelation, so its index of dispersion obeys

I = SCV + (SCV-1)*g2/(1-g2)

as reported in Section 5.2.2 of Casale, Mi, Cherkasova and Smirni, IEEE Trans. Soft. Eng. 37(5), 2011. The relation is inverted in closed form as g2 = (I-SCV)/(I-1) and the resulting decay rate is passed to map2_fit, which is the explicit inverse characterization of Heindl, Horvath and Gross. A third moment outside the feasible region is replaced by its lower limit (3/2)*e2^2/e1, the largest heavy-tail decay a MAP(2) admits.

The paper returns an exponential whenever SCV <= 1 or I < SCV, on the grounds that burstiness is then negligible. The rule does more than avoid an infeasible fit and must not be relaxed: a flow-equivalent server whose service is exponential and load dependent is exact for a product-form subnetwork by Norton’s theorem, whereas any MAP(2) fitted to the marginal inter-departure statistics is not, because the departure stream of the subnetwork is not independent of the rest of the model. Fitting the sub-exponential SCV of a non-bursty aggregate was measured to cost up to 2.2% of throughput on a three-station exponential network that the exponential fallback reproduces exactly.

Parameters:

e1 – Mean inter-arrival time
e2 – Second moment of the inter-arrival times
e3 – Third moment of the inter-arrival times
I – Asymptotic index of dispersion

Returns:

MAP – Fitted process in the form {D0,D1} status: 0 all four descriptors matched, 1 exponential as

See also

map2_fit, map_idc, fes_map_moments()

Examples

[MAP,status] = map2_fit_idc(e1,e2,e3,I)

dmap_to_dph(DMAP)

Discrete phase-type law underlying a renewal DMAP

Inverts dph_to_dmap: with D1 = a*alpha of rank one, alpha is recovered by normalizing any nonzero row of D1 and A is D0. Errors on a DMAP that does not renew, since no discrete phase-type law then describes its interevent times.

Parameters:: DMAP – Cell {D0, D1} of a renewal DMAP
Returns:: alpha – 1xm initial phase probability row vector A: mxm substochastic transient matrix

Examples

[alpha, A] = dmap_to_dph(DMAP)

dmap_super(A, B)

Superposition of two discrete-time arrival streams

Superposing two DMAPs on a slotted time scale is NOT a DMAP: both streams can fire in the same slot, so the merged stream carries batches. With A = {A_0,…,A_p} and B = {B_0,…,B_q} the merged batch matrices are E_k = sum_{i+j=k} kron(A_i, B_j), which is a discrete batch MAP of order (p+q). Folding E_2 into E_1 would conserve neither the arrival rate nor the slot in which the work appears, so the batch dimension is kept and the downstream station is solved as an M/G/1-type chain instead of a QBD.

Parameters:

A – Cell {A_0, A_1, …} of batch matrices (DMAP is {D0,D1})
B – Cell {B_0, B_1, …} of batch matrices

Returns:

E – Cell {E_0, E_1, …} of merged batch matrices

Examples

E = dmap_super(A, B)

dmap_compress_batch(B, maxOrder)

Reduces the order of a discrete batch arrival stream

Superposing slotted streams multiplies phase dimensions, so a decomposition over several stations has to bound them. The reduction keeps the two features that the downstream M/G/1-type solve consumes: the law of the time between NONEMPTY slots, matched to three moments by dmap_compress, and the stationary batch-size distribution conditional on a nonempty slot, kept exactly. The reduced stream is B_k = q_k * D1, B_0 = D0, whose event rate equals the original one by construction.

Parameters:

B – Cell {B_0, B_1, …} of batch matrices
maxOrder – Largest phase order left uncompressed

Returns:

Bc – Cell {B_0, B_1, …} of bounded phase order

Examples

Bc = dmap_compress_batch(B, maxOrder)

qbd_setupdelayoff(lambda, mu, alpharate, alphascv, betarate, betascv)

Analyzes queue length for system with setup and turn-off phases

This function performs queue-length analysis for a queueing system with setup delay and turn-off periods using QBD methods.

Parameters:

lambda – 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:

QN – Average queue length

Examples

QN = qbd_setupdelayoff(lambda, mu, alpharate, alphascv, betarate, betascv)

mmap_compress(MMAP, config)

Compresses a Marked MAP into a smaller representation

This function compresses an MMAP (Marked Markovian Arrival Process) into a smaller representation using various compression methods including mixture fitting, MAMAP2, and M3PP approaches.

Parameters:

MMAP – Original Marked Markovian Arrival Process
config – (Optional) Configuration struct with compression method

Returns:

MMAP – Compressed MMAP

Examples

MMAP = mmap_compress(MMAP)
MMAP = mmap_compress(MMAP, config)

mfq_multiregime(Q, R, Qt, Rt, T, pdfpoints, cdfpoints)

Solves a multi-regime feedback Markovian fluid queue.

Thin LINE wrapper around the multiregime routine implementing the method of H. E. Kankaya and N. Akar, “Solving Multi-Regime Feedback Fluid Queues”. The generator and drift rates are regime dependent, with separate boundary behaviour (feedback) generators/rates at each threshold.

Parameters:

Q – Cell of (N,N) generators, one per regime k=1..K
R – Cell of (N,N) diagonal drift-rate matrices per regime
Qt – Cell of boundary generators, k=0..K
Rt – Cell of boundary drift-rate matrices, k=0..K
T – Vector of regime thresholds (length K)
pdfpoints – Levels at which to evaluate density / density derivative
cdfpoints – Levels at which to evaluate the distribution function

Examples

[pdf,pdfd,cdf,cdfm] = mfq_multiregime(Q,R,Qt,Rt,T,pdfpoints,cdfpoints)

Returns:: pdf, pdfd (density derivative), cdf P(X<p), and cdfm P(X<=p).

mfq_ld_distr(masses, iniF, KF, cloF, iniB, KB, cloB, T, what, points)

Stationary fluid-level distribution from mfq_ld_solve output.

Thin LINE wrapper around LevelDependentFluidStationaryDistr. Evaluates the stationary density or distribution of a first/second-order level-dependent fluid queue at the requested points, using the matrix-exponential building blocks returned by mfq_ld_solve.

Parameters:

masses,iniF..cloB – Building blocks returned by mfq_ld_solve
T – Vector of regime thresholds (length K)
what – ‘pdf’, ‘pdfd’ (density derivative), ‘cdf’ P(X<p), or ‘cdfm’ P(X<=p)
points – Fluid levels at which to evaluate

Examples

res = mfq_ld_distr(masses,iniF,KF,cloF,iniB,KB,cloB,T,what,points)

Returns:: res – (numel(points), N) matrix with the per-state values at each point.

qbd_mapmap1(MAPa, MAPs, util)

Analyzes MAP/MAP/1 queue using Quasi-Birth-Death process

This function solves a MAP/MAP/1 queue using QBD (Quasi-Birth-Death) methods, computing throughput, queue length, utilization, and other performance metrics.

Parameters:

MAPa – Arrival process (MAP)
MAPs – Service process (MAP)
util – (Optional) Target utilization to scale service rate

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 A_1: Downward transition block A0: Local transition block A1: Upward transition block U: Matrix U MAPs: Scaled service process

Examples

[XN, QN, UN, pqueue, R, eta, G, A_1, A0, A1, U, MAPs] = qbd_mapmap1(MAPa, MAPs)
[XN, QN, UN, pqueue, R, eta, G, A_1, A0, A1, U, MAPs] = qbd_mapmap1(MAPa, MAPs, util)

[XN,QN,UN,PQUEUE,R,ETA]=QBD_MAPMAP1(MAPA,MAPS,UTIL)

qbd_depproc_etaqa_ps(MAPa, MAPs, n)

Builds a MAP approximation of the departure process for a MAP/MAP/1 PS queue

This function constructs a finite-state MAP {D0, D1} representation of the departure process from a MAP/MAP/1 queue with Processor Sharing (PS) discipline, using ETAQA truncation at QBD level n.

Compared to the FCFS variant (qbd_depproc_etaqa), the PS discipline splits service completions at level j into a departure component B*(1/j) and an internal transition component B*(1-1/j), reflecting the rate-dependent sharing of the server among j jobs. The tail approximation at level n uses Bbar and Bhat matrices weighted by 1/n and (n-1)/n respectively.

Parameters:

MAPa – Arrival process in MAP format {D0, D1}
MAPs – Service process in MAP format {D0, D1}
n – Truncation level (number of QBD levels to represent explicitly)

Returns:

D – Departure process in MAP format {D0, D1}

Examples

D = qbd_depproc_etaqa_ps(MAPa, MAPs, n)

qbd_depproc_etaqa(MAPa, MAPs, n)

Builds a MAP approximation of the departure process for a MAP/MAP/1 FCFS queue

This function constructs a finite-state MAP {D0, D1} representation of the departure process from a MAP/MAP/1 queue with FCFS discipline, using ETAQA truncation at QBD level n. The resulting MAP captures both the SCV (squared coefficient of variation) and autocorrelation structure of inter-departure times.

The QBD process is formed with forward matrix F = kron(MAPa{2}, I_ns), local matrix L = krons(MAPa{1}, MAPs{1}), and backward matrix B = kron(I_na, MAPs{2}). Levels 0..n-1 are represented explicitly, and level n uses a tail approximation via the G matrix.

Parameters:

MAPa – Arrival process in MAP format {D0, D1}
MAPs – Service process in MAP format {D0, D1}
n – Truncation level (number of QBD levels to represent explicitly)

Returns:

D – Departure process in MAP format {D0, D1}

Examples

D = qbd_depproc_etaqa(MAPa, MAPs, n)

qbd_bmapbmap1(MAPa, pbatcha, MAPs)

Constructs QBD blocks for batch arrival/service queue

This function computes the QBD (Quasi-Birth-Death) transition blocks for a BMAP/BMAP/1 queue with batch arrivals and batch service.

Parameters:

MAPa – Arrival process (MAP)
pbatcha – Probability distribution of batch sizes
MAPs – Service process (MAP)

Examples

qbd_bmapbmap1(MAPa, pbatcha, MAPs)

qbd_R_logred(B, L, F, iter_max)

Computes QBD rate matrix R via logarithmic reduction

This function computes the rate matrix R for a Quasi-Birth-Death (QBD) process using the logarithmic reduction method.

Parameters:

B – Backward transition block A_(-1)
L – Local transition block A_0
F – Forward transition block A_1
iter_max – (Optional) Maximum iterations (default: 100000)

Returns:

R – Rate matrix R

Examples

R = qbd_R_logred(B, L, F)
R = qbd_R_logred(B, L, F, iter_max)

Logarithmic reduction method

qbd_R(B, L, F, iter_max)

Computes QBD rate matrix R via successive substitutions

This function computes the rate matrix R for a Quasi-Birth-Death (QBD) process using the successive substitutions method.

Parameters:

B – Backward transition block A_(-1)
L – Local transition block A_0
F – Forward transition block A_1
iter_max – (Optional) Maximum iterations (default: 100000)

Returns:

R – Rate matrix R

Examples

R = qbd_R(B, L, F)
R = qbd_R(B, L, F, iter_max)

Successive substitutions method

mmdp_isfeasible(Q, R)

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

bool = MMDP_ISFEASIBLE(Q, R)

Checks if the given matrices define a valid Markov-Modulated Deterministic Process (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 (matrix) – n×n generator matrix
R (matrix) – n×n diagonal rate matrix

Returns:

bool (logical) – True if (Q, R) defines a valid MMDP

map_m1ps_h_recursive(C, D, mu, N, K)

MAP_M1PS_H_RECURSIVE Recursive computation of h_{n,k} coefficients

h = MAP_M1PS_H_RECURSIVE(C, D, mu, N, K) computes the h_{n,k} vectors used in the sojourn time distribution for a MAP/M/1-PS queue.

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

Input:

C - M x M matrix governing MAP transitions without arrivals D - M x M matrix governing MAP transitions with arrivals mu - Service rate (scalar) N - Maximum value of n to compute (determines rows) K - Maximum value of k to compute (determines columns)

Output:

h - Cell array of size (N+1) x (K+1): h{n+1, k+1} contains the M x 1 vector h_{n,k} (indices shifted by 1 for MATLAB 1-based indexing)

Reference:

Masuyama, H., & Takine, T. (2003). Sojourn time distribution in a MAP/M/1 processor-sharing queue. Operations Research Letters, 31(6), 406-412.

See also: MAP_M1PS_SOJOURN, MAP_COMPUTE_R

map_jointpdf_derivative(MAP, iset)

Computes partial derivative at zero of MAP joint PDF

This function computes the partial derivative at 0 of a MAP’s joint probability density function. Based on A. Horvath et al. “A Joint Moments Based Analysis of Networks of MAP/MAP/1 Queues”.

Parameters:

MAP – Markovian Arrival Process {D0, D1}
iset – Vector of indices defining the partial derivative

Returns:

gamma – Calculated derivative value

Examples

gamma = map_jointpdf_derivative(MAP, iset)

partial derivative at 0 of a MAP’s joint PDF A. Horvath et al. A Joint Moments Based Analysis of Networks of MAP/MAP/1 Queues

map_compute_R(C, D, mu)

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

R = MAP_COMPUTE_R(C, D, mu) computes the matrix R, which is the minimal nonnegative solution of the matrix equation:

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

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

Input:: C - M x M matrix governing MAP transitions without arrivals D - M x M matrix governing MAP transitions with arrivals mu - Service rate (scalar)
Output:: R - M x M rate matrix (minimal nonnegative solution)
Reference:: Masuyama, H., & Takine, T. (2003). Sojourn time distribution in a MAP/M/1 processor-sharing queue. Operations Research Letters, 31(6), 406-412.

See also: MAP_M1PS_SOJOURN, MAP_M1PS_H_RECURSIVE

map_ccdf_derivative(MAP, i)

Computes derivative at zero of MAP complementary CDF

This function computes the derivative at 0 of a MAP’s Complementary Cumulative Distribution Function (CCDF). Based on A. Horvath et al. “A Joint Moments Based Analysis of Networks of MAP/MAP/1 Queues”.

Parameters:

MAP – Markovian Arrival Process {D0, D1}
i – Order of the derivative

Returns:

nu – Calculated derivative value

Examples

nu = map_ccdf_derivative(MAP, i)

derivative at 0 of a MAP CCDF A. Horvath et al. A Joint Moments Based Analysis of Networks of MAP/MAP/1 Queues

ldqbd_R(Q0, Q1, Q2, varargin)

Computes all rate matrices R^(n) for a level-dependent QBD

This function implements the backward recursion algorithm for computing rate matrices of level-dependent QBD processes with heterogeneous dimensions.

For a level-dependent QBD with levels 0, 1, …, N: R^(N) = Q0^(N-1) * (-Q1^(N))^{-1} R^(n) = Q0^(n-1) * (-Q1^(n) - R^(n+1) * Q2^(n+1))^{-1} for n = N-1,…,1

Dimensions: R^(n) is (states at level n-1) x (states at level n) Q0^(n-1) is (states at level n-1) x (states at level n) Q1^(n) is (states at level n) x (states at level n) Q2^(n) is (states at level n) x (states at level n-1)

Parameters:

Q0 – Cell array {Q0^(0), Q0^(1), …, Q0^(N-1)} - upward transitions
Q1 – Cell array {Q1^(0), Q1^(1), …, Q1^(N)} - local transitions
Q2 – Cell array {Q2^(1), Q2^(2), …, Q2^(N)} - downward transitions
options – Optional struct with fields: verbose (default false)

Returns:

R – Cell array {R^(1), R^(2), …, R^(N)} - rate matrices

Examples

R = ldqbd_R(Q0, Q1, Q2)
R = ldqbd_R(Q0, Q1, Q2, options)

References

T. Phung-Duc, H. Masuyama, S. Kasahara, Y. Takahashi, “A Simple Algorithm for the Rate Matrices of Level-Dependent QBD Processes”, QTNA 2010.

ldqbd(Q0, Q1, Q2, varargin)

Solves level-dependent QBD processes using matrix continued fractions

This function implements Algorithm 1 from “A Simple Algorithm for the Rate Matrices of Level-Dependent QBD Processes” by Phung-Duc, Masuyama, Kasahara, and Takahashi (2010), QTNA Conference.

It computes the rate matrices R^(n) for an ergodic level-dependent QBD process with a finite number of levels, and optionally the stationary distribution.

For a level-dependent QBD with levels 0, 1, …, N, the infinitesimal generator has block-tridiagonal structure:

Q^(0)_1 Q^(0)_0 O O … O Q^(1)_2 Q^(1)_1 Q^(1)_0 O … O O Q^(2)_2 Q^(2)_1 Q^(2)_0 … O … … … … … … O O O … Q^(N)_2 Q^(N)_1

where Q_0^(n) are upward transitions (level n to n+1), Q_1^(n) are local transitions (within level n), and Q_2^(n) are downward transitions (level n to n-1).

Parameters:

Q0 – Cell array {Q0^(0), Q0^(1), …, Q0^(N-1)} - upward transitions
Q1 – Cell array {Q1^(0), Q1^(1), …, Q1^(N)} - local transitions
Q2 – Cell array {Q2^(1), Q2^(2), …, Q2^(N)} - downward transitions
options – Optional struct with fields: epsilon (default 1e-10),

Returns:

R – Cell array {R^(1), R^(2), …, R^(N)} - rate matrices pi: Row vector [pi_0, pi_1, …, pi_N] - stationary distribution

See also HYPEREXP, MAP_FIT.

mexify_mam

@brief MATLAB Coder script to generate MEX functions for mam_ module.

This script generates MEX (MATLAB Executable) versions of MAM (Matrix Analytic Methods) functions for improved performance.

Skipped functions (Coder-incompatible):: ldqbd, ldqbd_R, ldqbd_pi - cell array arguments qbd_rg, qbd_mapmap1, qbd_raprap1 - cell array arguments qbd_bmapbmap1 - cell array arguments map_ccdf_derivative - cell/struct MAP argument map_jointpdf_derivative - cell/struct MAP argument map_m1ps_h_recursive - cell return values map_m1ps_sojourn - function handle validators map_m1ps_cdfrespt - function handles + nested functions mmap_compress - struct/cell arguments

See also CODER, CODER.CONFIG, CODER.TYPEOF, CODEGEN.

qbd_depproc_jointmom(MAPa, MAPs, iset)

Computes joint moments E[X_0^i * X_1^j] of consecutive inter-departure times

Given arrival and service MAPs for a MAP/MAP/1 queue, this function computes joint moments of consecutive inter-departure times using the QBD structure. The initial vector is constructed from the stationary distribution at departure epochs via ETAQA.

The departure SCV and lag-1 ACF can be obtained as: E1 = qbd_depproc_jointmom(MAPa, MAPs, [1,0]) % E[X] E2 = qbd_depproc_jointmom(MAPa, MAPs, [2,0]) % E[X^2] E11 = qbd_depproc_jointmom(MAPa, MAPs, [1,1]) % E[X_0*X_1] SCV = (E2 - E1^2) / E1^2 ACF = (E11 - E1^2) / (E2 - E1^2)

Parameters:

MAPa – Arrival process in MAP format {D0, D1}
MAPs – Service process in MAP format {D0, D1}
iset – Matrix of moment orders [i1,j1; i2,j2; …] where each

Returns:

JM – Vector of joint moments, one per row of iset

Examples

JM = qbd_depproc_jointmom(MAPa, MAPs, iset)

mg1_dt_queue(ARV, SVC, options)

Discrete-time single-server queue with batch DMAP arrivals

Solves the DBMAP/DMAP/1 queue on a slotted time scale under the late arrival system with delayed access (LAS-DA): within a slot the service completion resolves first, arrivals are appended at the end of the slot and cannot enter service before the next slot, and the level is read after both. This is the convention of Q_DT_MAP_MAP_1, whose QBD blocks this function reproduces for a single arrival per slot, and of the LDES slotted engine.

The chain is M/G/1-type because a slot may deliver a batch: with arrival matrices A_k and service pair (S0,S1), A^(-1) = kron(A_0,S1), A^(k) = kron(A_k,S0) + kron(A_{k+1},S1), B^(k) = kron(A_k,I), the boundary row holding the empty system, where no service runs.

Parameters:

ARV – Cell {A_0, A_1, …} of arrival batch matrices
SVC – Cell {S0, S1} of the service DMAP
options – (Optional) reads config.space_max and config.dt_maxlevel

Returns:

QN – Mean number in system at slot boundaries (LAS-DA) UN: Utilization, i.e. fraction of slots with the server busy TN: Throughput in jobs per slot ql: Queue length pmf, ql(i) = Prob[i-1 in system] DEP: Departure process as a DMAP {D0,D1} on the truncated chain

Examples

[QN, UN, TN, ql] = mg1_dt_queue(ARV, SVC)
[QN, UN, TN, ql, DEP] = mg1_dt_queue(ARV, SVC, options)

dph_to_dmap(alpha, A)

Renewal DMAP of a discrete phase-type law

Maps (alpha, A) to {D0, D1} with D0 = A and D1 = a*alpha, a = e - A*e. The resulting DMAP renews the phase at every event, so its interevent times are i.i.d. copies of the DPH and D0+D1 is stochastic by construction.

Parameters:

alpha – 1xm initial phase probability row vector
A – mxm substochastic transient matrix

Returns:

DMAP – Cell {D0, D1}

Examples

DMAP = dph_to_dmap(alpha, A)

dph_from_dist(procType, meanSlots, scv)

Exact discrete phase-type representation of a lattice-valued law

Returns (alpha, A) with P[X=k] = alpha*A^(k-1)*a, a = e - A*e, k = 1,2,…, measuring X in slots. The three families below are represented EXACTLY, not moment-matched: Geometric by its single-phase chain, Det(k) by the k-step chain, DiscreteUniform by the hazard chain h_j = 1/(hi-j+1) on j >= lo. A family outside this set is an error rather than a fit, because a fitted surrogate would silently leave the lattice the caller is relying on.

Parameters:

procType – ProcessType id (GEOMETRIC, DET, DUNIFORM)
meanSlots – Mean of the law expressed in slots
scv – Squared coefficient of variation (used by DUNIFORM only)

Returns:

alpha – 1xm initial phase probability row vector A: mxm substochastic transient matrix

Examples

[alpha, A] = dph_from_dist(procType, meanSlots, scv)

dmap_thin(A, p)

Bernoulli thinning of a discrete-time batch arrival stream

Routes each event of a discrete batch MAP to the branch with probability p, independently across events. Since a slot may carry a batch of n events, the number routed is Binomial(n,p) and the thinned matrices are B_k = sum_{n>=k} nchoosek(n,k) p^k (1-p)^(n-k) A_n. The phase process is untouched, so the result is exact for Bernoulli (PROB/RAND) routing.

Parameters:

A – Cell {A_0, A_1, …} of batch matrices
p – Routing probability in [0,1]

Returns:

B – Cell {B_0, B_1, …} of thinned batch matrices

Examples

B = dmap_thin(A, p)

dmap_lambda(A)

Event rate per slot of a discrete-time batch arrival stream

Returns pi * sum_k k*A_k * e with pi the stationary vector of the phase chain sum_k A_k. This counts EVENTS per slot, not slots with at least one event, so it is the quantity that Little’s law consumes downstream. For a plain DMAP {D0,D1} it reduces to pi*D1*e.

Parameters:: A – Cell {A_0, A_1, …} of batch matrices
Returns:: lambda – Mean number of events per slot

Examples

lambda = dmap_lambda(A)

dmap_is_renewal(DMAP)

Tests whether a DMAP is a renewal process

A DMAP renews at every event exactly when D1 has rank one, i.e. D1 = a*alpha with a the absorption vector of D0. Interevent times are then i.i.d. and the process is the DMAP form of a discrete phase-type law, which is what lets the caller use the cheaper Q_DT_PH_PH_1 entry point.

Parameters:: DMAP – Cell {D0, D1}
Returns:: bool – True when D1 has rank one

Examples

bool = dmap_is_renewal(DMAP)

dmap_compress(DMAP, maxOrder)

Reduces the order of a DMAP by matching interevent moments

Leaves the process untouched while its order stays within maxOrder, and otherwise replaces it by the two-phase discrete phase-type law with the same first three interevent moments (BuTools DPH2From3Moments), read back as a renewal DMAP. Correlation is NOT preserved, which mirrors what the continuous-time ‘mixture.order1’ compression does, and it is the reason the multi-station discrete-time path is an approximation. When the moment triple is outside the DPH(2) region the fallback keeps the exact mean with a Geometric, so the arrival rate of the decomposition is conserved in every branch.

Parameters:

DMAP – Cell {D0, D1}
maxOrder – Largest order left uncompressed

Returns:

DMAPc – Cell {D0, D1} of order at most max(maxOrder,2)

Examples

DMAPc = dmap_compress(DMAP, maxOrder)

qbd_raprap1(RAPa, RAPs, util)

Analyzes RAP/RAP/1 queue using Quasi-Birth-Death process

This function solves a RAP/RAP/1 queue (Rational Arrival Process) using QBD methods, computing throughput, queue length, utilization, and other performance metrics.

The two RAPs are INDEPENDENT of each other, so the QBD phase space is the product of the two phase spaces and the blocks factor as Kronecker products. This function is the thin product-space wrapper; the analysis itself is the block-level core qbd_rap.m, which solves an arbitrary QBD with RAP components. A model whose arrival process and sequence of service times share a phase space, and are therefore cross-correlated, has no such product structure and must call qbd_rap directly.

Parameters:

RAPa – Arrival process (RAP)
RAPs – Service process (RAP)
util – (Optional) Target utilization to scale service rate

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

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.m.

Phase ordering:: The QBD phase is the pair (arrival phase, service phase) laid out with the ARRIVAL phase major and the service phase minor, i.e. the phase index is (a-1)*ns + s. That is the ordering produced by kron(RAPa, eye(ns)) and kron(eye(na), RAPs), and pqueue is indexed by it downstream in solver_mam_basic, so the two Kronecker factors must not be swapped.

Examples

[XN, QN, UN, pqueue, R, eta, G, B, L, F] = qbd_raprap1(RAPa, RAPs)
[XN, QN, UN, pqueue, R, eta, G, B, L, F] = qbd_raprap1(RAPa, RAPs, util)

[XN,QN,UN,PQUEUE,R,ETA]=QBD_RAPRAP1(RAPA,RAPS,UTIL)

qbd_rap(A0, A1, A2, B0, B1, numLevels)

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.

Parameters:

A0 – Level-up block (m x m)
A1 – Local block (m x m)
A2 – Level-down block (m x m)
B0 – (Optional) boundary level-up block, default A0
B1 – (Optional) boundary local block, default A1
numLevels – (Optional) highest level reported, default 20

Returns:

levelProb – Row vector of marginal level probabilities, levels 0..numLevels QN: Mean queue length (exact, pi0*R*inv(I-R)^2*e) R: Rate matrix R G: Matrix G solving A0*G^2 + A1*G + A2 = 0 U: Matrix U = A1 + A0*G spr: Spectral radius Sp(R) pqueue: (numLevels+1) x m matrix of the vectors pi_n pi0: Level-0 vector pi_0, the boundary vector of Theorem 7

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 equilibrium construction below is their Theorem 7 and the stability test is their Corollary 8. 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.m.

Algorithm (Theorem 7):

Solve A0*G^2 + A1*G + A2 = 0 for G.
U = A1 + A0*G.
R = A0*inv(-U).
Find the row vector pihat0 with pihat0*(B1 + R*A2) = 0, pihat0*e = 1.
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). Two paths are therefore taken: - If A2 has rank one, A2 = u*v, then G = e*v/(v*e) solves the equation

exactly. This is immediate from conservativity: G is idempotent and (A0+A1)*e = -A2*e = -u*(v*e), so (A0+A1)*e*v/(v*e) = -u*v = -A2. This is the case covered by the paper’s example, and the residual is checked.

Otherwise the quadratic matrix equation is solved numerically by natural functional iteration G <- inv(-A1)*(A2 + A0*G^2) used as a warm start, followed by Newton’s method on the Sylvester-form Jacobian, (A0*G+A1)*H + A0*H*G = -(A0*G^2 + A1*G + A2), solved through its Kronecker expansion. If the residual does not reach roundoff level, or the iterate does not satisfy the stochastic-analogue constraint G*e = e, an error is raised rather than returning an unconverged G.

Examples

[levelProb, QN, R, G, U, spr, pqueue, pi0] = qbd_rap(A0, A1, A2)
[levelProb, QN, R, G, U, spr, pqueue, pi0] = qbd_rap(A0, A1, A2, B0, B1)
[levelProb, QN, R, G, U, spr, pqueue, pi0] = qbd_rap(A0, A1, A2, B0, B1, numLevels)

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, in which arrivals and
services share a phase space, must call qbd_rap directly because no
product form exists to factor out.

g  = 0.25;
A1 = [-1 0 0; -2/3 -1 1; 2/3 -1 -1];
Da = [14/5 -9/10 -9/10; 26/15 -8/15 -8/15; 58/15 -19/15 -19/15];
Ds = [1; 2/3; 4/3]*[3 -1 -1];
[levelProb, QN] = qbd_rap(g*Da, A1, (1-g)*Ds, g*Da, g*A1, 8);
% levelProb reproduces Table 1 of Bean and Nielsen (2010)

[LEVELPROB,QN,R,G,U,SPR,PQUEUE,PI0]=QBD_RAP(A0,A1,A2,B0,B1,NUMLEVELS)

map_m1ps_cdfrespt(C, D, mu, x, varargin)

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

W_bar = MAP_M1PS_CDFRESPT(C, D, mu, x) computes the complementary distribution function of the sojourn time in a MAP/M/1 processor-sharing queue at the points specified in x.

[W_bar, W_bar_n] = MAP_M1PS_CDFRESPT(C, D, mu, x) also returns the conditional complementary distributions W_bar_n{i} for customers finding i-1 customers in the system on arrival.

[…] = MAP_M1PS_CDFRESPT(C, D, mu, x, ‘Param’, Value) specifies optional :Parameters: * ‘Epsilon’ - Truncation parameter for queue length (default – 1e-11)

‘EpsilonPrime’ - Truncation parameter for uniformization (default – 1e-10)

‘Verbose’ - Display computation progress (default – false)

Input:: C - M x M matrix governing MAP transitions without arrivals D - M x M matrix governing MAP transitions with arrivals mu - Service rate (scalar, mu > 0) x - Vector of time points at which to evaluate W_bar(x) = Pr[W > x]
Output:: W_bar - Vector of same size as x, containing Pr[W > x] W_bar_n - Cell array of conditional distributions (optional)

The processor-sharing (PS) discipline shares the server equally among all customers. When n customers are present, each receives service at rate 1/n.

The algorithm implements Theorem 1 from:: Masuyama, H., & Takine, T. (2003). Sojourn time distribution in a MAP/M/1 processor-sharing queue. Operations Research Letters, 31(6), 406-412.

Example

% M/M/1-PS queue with lambda=0.8, mu=1 lambda = 0.8; mu = 1; C = -lambda; D = lambda; % Poisson arrivals x = linspace(0, 10, 100); W_bar = map_m1ps_cdfrespt(C, D, mu, x); plot(x, W_bar); xlabel(‘x’); ylabel(‘Pr[W > x]’); title(‘Sojourn time distribution for M/M/1-PS’);

qbd_fundmat(B, L, F, matrices, precision, maxNumIt)

QBD_FUNDMAT Cyclic-reduction fundamental matrices G (and optionally R) of a QBD.

Operates directly on the raw level blocks (B,L,F) of a homogeneous QBD and supports complex-valued L (as required when L is shifted by -s*I in a Laplace-domain transient analysis). This differs from QBD_RG, which builds the blocks from MAP representations and is only exercised at real argument.

VARARGOUT = QBD_FUNDMAT(B, L, F, MATRICES) returns the matrices named by the character codes in MATRICES (default ‘G’), in that order. Supported codes: ‘G’, ‘R’, ‘GR’, ‘RG’.

The G matrix is obtained by cyclic reduction (Bini-Meini logarithmic reduction); R is recovered from G as R = Fm*(I-(Lm+Fm*G))^-1.

mfq_sojourn(Q, Rin, Rout, Q0, transToPH)

Sojourn-time distribution of a fluid queue as an ME/PH representation.

Thin LINE wrapper around FluidQueueSTD. Returns the distribution of the time a fluid drop spends in a Markov-modulated fluid queue with input rate matrix Rin and output (service) rate matrix Rout, modulated by generator Q, as a matrix-exponential (ME) or phase-type (PH) representation (alpha, A).

Parameters:

Q – (N,N) generator of the background Markov chain
Rin – (N,N) diagonal input fluid-rate matrix
Rout – (N,N) diagonal output (service) fluid-rate matrix
Q0 – (Optional) level-0 generator (default: Q)
transToPH – (Optional) true to return a PH (else ME) representation

Examples

[alpha,A] = mfq_sojourn(Q,Rin,Rout)
[alpha,A] = mfq_sojourn(Q,Rin,Rout,Q0,transToPH)

Returns:: alpha, A – ME/PH representation of the sojourn-time distribution.

mfq_prio_queue(Q, R, d, varargin)

Analyzes an MMAP[K]/PH[K]/1-type fluid priority queue.

Thin LINE wrapper around the BUTools-family routine FluidPrioQueue, implementing the method of G. Horvath, “Efficient analysis of the MMAP[K]/PH[K]/1 priority queue”, EJOR 246(1):128-139, 2015. A background Markov chain with generator Q modulates the per-class fluid input rates (matrix R, one row per priority class) and the fluid is drained at the constant service rate d, higher-priority fluid first.

Parameters:

Q – (N,N) generator of the modulating Markov chain
R – (K,N) per-class fluid input rates in the background states
d – Constant fluid service rate (scalar, positive)
… – Measure/option pairs: ‘flMoms’,’flDistr’,’stMoms’,’stDistr’,

Examples

varargout = mfq_prio_queue(Q, R, d, ...)
[flM] = mfq_prio_queue(Q, R, d, 'flMoms', n)
[cdf] = mfq_prio_queue(Q, R, d, 'stDistr', points, 'erlMaxOrder', L)

Returns:: One output per requested performance measure; each column corresponds to a priority class.

mfq_ld_solve(varargin)

Matrix-exponential solution of a multi-regime Markovian fluid queue.

Thin LINE wrapper around SecondOrderLevelDependentFluidSolve. The generator, drift and (optionally) variance change at threshold fluid levels, yielding a piecewise-homogeneous first- or second-order (Brownian) fluid queue. Setting the variance cells S to zero reduces the model to first order. The returned matrix-exponential building blocks are consumed by mfq_ld_distr / mfq_ld_mean.

Parameters:

Q – Cell of (N,N) generators, one per regime
R – Cell of (N,N) diagonal drift-rate matrices per regime
S – Cell of (N,N) diagonal variance matrices per regime (0 = first order)
T – Vector of regime thresholds (length K)
boundaryL/U – (Optional) per-state boundary flag (0 reflective, 1 absorbing)
Qt – (Optional) cell of boundary generators
prec – (Optional) numerical precision (default 1e-14)

Examples

[masses,iniF,KF,cloF,iniB,KB,cloB] = mfq_ld_solve(Q,R,S,T)
[...] = mfq_ld_solve(Q,R,S,T,boundaryL,boundaryU,Qt,prec)

Returns:: masses (K+1 point-mass vectors) and per-regime forward/backward matrix-exponential parameters iniF,KF,cloF,iniB,KB,cloB.

mfq_ld_mean(masses, iniF, KF, cloF, iniB, KB, cloB, T)

Stationary mean fluid level from mfq_ld_solve output.

Thin LINE wrapper around LevelDependentFluidStationaryMean. Returns the scalar mean fluid level E[X] of a first/second-order level-dependent fluid queue in closed form, using the matrix-exponential building blocks returned by mfq_ld_solve.

Parameters:

masses,iniF..cloB – Building blocks returned by mfq_ld_solve
T – Vector of regime thresholds (length K)

Examples

res = mfq_ld_mean(masses,iniF,KF,cloF,iniB,KB,cloB,T)

Returns:: res – scalar mean fluid level E[X].

mfq_fluflu_sojourn(Qin, Rin, Qout, Rout, srv0stop, transToPH)

Sojourn-time distribution of a fluid queue with fluid-modulated service.

Thin LINE wrapper around FluFluSTD. Returns the sojourn-time distribution of a fluid queue in which both the arrival and the service processes are Markov-modulated fluid flows, as a matrix-exponential (ME) or phase-type (PH) representation (alpha, A).

Parameters:

Qin – (Na,Na) generator of the arrival-modulating chain
Rin – (Na,Na) diagonal arrival fluid-rate matrix
Qout – (Ns,Ns) generator of the service-modulating chain
Rout – (Ns,Ns) diagonal service fluid-rate matrix
srv0stop – true if service stops when the server fluid level hits zero
transToPH – (Optional) true to return a PH (else ME) representation

Examples

[alpha,A] = mfq_fluflu_sojourn(Qin,Rin,Qout,Rout,srv0stop)
[alpha,A] = mfq_fluflu_sojourn(Qin,Rin,Qout,Rout,srv0stop,transToPH)

Returns:: alpha, A – ME/PH representation of the sojourn-time distribution.

qbd_rg(MAPa, MAPs, util)

Computes QBD transition matrices for MAP/MAP/1 queue

This function computes the QBD (Quasi-Birth-Death) matrices R and G along with transition blocks B, L, F for a MAP/MAP/1 queue.

Parameters:

MAPa – Markovian Arrival Process for arrivals
MAPs – Markovian Arrival Process for service
util – (Optional) Target utilization for scaling service rate

Returns:

R – Rate matrix R (forward transitions) G: Rate matrix G (backward transitions) B: Backward transition matrix L: Local transition matrix F: Forward transition matrix U: Utilization

Examples

[R, G, B, L, F, U] = qbd_rg(MAPa, MAPs)
[R, G, B, L, F, U] = qbd_rg(MAPa, MAPs, util)

[XN,QN,UN,PQUEUE,R,ETA]=QBD_MAPMAP1(MAPA,MAPS,UTIL)