api.qsys

qsys_mtgs0_mol(lambdaFun, serviceCcdf, ES, s, tvals, varargin)

QSYS_MTGS0_MOL Modified-offered-load approximation for a time-varying system.

RESULT = QSYS_MTGS0_MOL(LAMBDAFUN, SERVICECCDF, ES, S, TVALS) approximates the blocking probability of the Mt/G/S/0 loss system at the times TVALS, with arrival rate LAMBDAFUN, service ccdf SERVICECCDF of mean ES.

THE ONE IDEA. A stationary loss system with offered load a blocks with probability B(s,a). In a time-varying system the question is WHICH LOAD to put in that formula. The pointwise stationary approximation (PSA) uses the instantaneous one, lambda(t)ES. The modified offered load (MOL) uses the offered load of the corresponding INFINITE-SERVER system,

m(t) = ES E[lambda(t - Se)] = int_0^Inf lambda(t-x)P(S>x)dx,

which is EXACT for that system and therefore carries the time lag and the smoothing the finite-server system also has. MOL is then B(s,m(t)). The difference between the two is precisely the lag: PSA peaks when the arrival rate peaks, MOL peaks later, and the real system peaks later too.

WHAT TO EXPECT. Measured against the exact time-varying birth-death chain on a sinusoidal rate, MOL cuts the mean RELATIVE error roughly threefold (0.13 against 0.44 at s = 100) because it gets the phase right. It does not always win on ABSOLUTE error, which is dominated by the peak of the cycle where both are weakest. Under constant input MOL is exact.

Options: ‘delay’ (use Erlang C, i.e. the delay probability of an Mt/M/s queue rather than the blocking probability of a loss system), plus every option of QSYS_MTGINF, which computes the offered load.

Returns a struct with fields times, offeredLoad, instantLoad, probBlockMOL, probBlockPSA, meanBusyMOL and arrivalRate.

Reference: W. A. Massey, W. Whitt (1994). An analysis of the modified offered load approximation for the nonstationary Erlang loss model. Annals of Applied Probability 4(4), 1145-1160; W. Whitt (1991). The pointwise stationary approximation for Mt/Mt/s queues is asymptotically correct as the rates increase. Management Science 37(3), 307-314.

See also QSYS_MTGINF, QSYS_MMK_QED, LOSSN_ERLANGFP.

qsys_mtginf(lambdaFun, serviceCcdf, ES, tvals, varargin)

QSYS_MTGINF Exact time-varying analysis of the Mt/G/Inf queue.

RESULT = QSYS_MTGINF(LAMBDAFUN, SERVICECCDF, ES, TVALS) evaluates the infinite-server queue with a non-homogeneous Poisson arrival rate LAMBDAFUN and iid service times with complementary cdf SERVICECCDF and mean ES, at the times TVALS.

THE RESULT IS EXACT, not an approximation. With infinitely many servers customers never interact, so the model is a Poisson random measure and the number in system at time t is POISSON with mean

m(t) = E[ int_{t-S}^{t} lambda(u) du ] = ES * E[lambda(t - Se)]

= int_0^Inf lambda(t-x) P(S > x) dx

where Se is the STATIONARY-EXCESS (equilibrium) law of the service time, with density P(S>x)/ES. Because the law is Poisson, the variance equals the mean and every quantile follows from it.

THE PHYSICS. Reading m(t) as ES*E[lambda(t-Se)] says the time-varying load is the stationary load ES*lambda(t) subjected to a TIME LAG and a SPACE SHIFT: to first order m(t) ~ ES*lambda(t - E[Se]) with E[Se] = E[S^2]/(2*ES), so peak congestion LAGS peak arrival rate, and by more than the mean service time when the service law is variable. The pointwise stationary approximation ES*lambda(t) is the zeroth-order term of the same expansion, which is exactly why it misses the lag.

Options:

‘startTime’, T0 - the system started empty at T0; the default -Inf assumes

the arrival rate has been running forever

‘ES2’, M2 - the second moment of the service time, which adds the lag

E[Se] and the first-order lag approximation to the output

‘servicePdf’, G - the service density, used for the exact departure rate;

without it the departure rate comes from the flow balance m’(t) = lambda(t) - delta(t) by a central difference

‘tol’, TOL - service-tail cut for the age integral, default 1e-12 ‘panels’, N - Simpson panels for that integral, default 4000 ‘maxAge’, A - cap on the age integrated over, default 1e12

Returns a struct with fields:

times - the requested times meanNumber - m(t), the Poisson mean varNumber - equal to meanNumber, the law being Poisson arrivalRate - lambda(t) departureRate - delta(t) = E[lambda(t-S)] offeredLoadPSA - ES*lambda(t), the pointwise stationary approximation meanLag - E[Se], when ES2 is given lagApproximation - ES*lambda(t-E[Se]), when ES2 is given

Example

% sinusoidal arrivals, exponential service of rate 2 res = qsys_mtginf(@(t) 10+5*sin(t), @(x) exp(-2*x), 1/2, linspace(0,2*pi,50));

Reference: S. G. Eick, W. A. Massey, W. Whitt (1993). The physics of the Mt/G/infinity queue. Operations Research 41(4), 731-742.

See also QSYS_MGINF, QSYS_ERLANGA, QSYS_GGISGI_FLUID.

qsys_mmk_qed_staffing(lambda, mu, target, varargin)

QSYS_MMK_QED_STAFFING Square-root staffing of the M/M/s queue.

RESULT = QSYS_MMK_QED_STAFFING(LAMBDA, MU, TARGET) returns the smallest number of servers whose Halfin-Whitt delay probability is at most TARGET, together with the quality-of-service parameter behind it.

THE RULE. Invert alpha(beta) = TARGET for the server slack beta, then staff

s = ceil( a + beta sqrt(a) ), a = lambda/mu,

the SQUARE-ROOT STAFFING rule: the base a erlangs of work plus a safety cushion that grows only as the square root of the load. Doubling the load needs only sqrt(2) times the cushion, which is why large service systems can be both highly utilized and responsive, and why the QED regime is the one large systems are actually run in.

RESULT = QSYS_MMK_QED_STAFFING(…, ‘criterion’, C) chooses what TARGET means:

‘delay’ - P(W > 0) <= TARGET (the default) ‘meanwait’ - E[W] <= TARGET, in time units ‘servicelevel’ - P(W <= TARGET.deadline) >= TARGET.level, TARGET being a

struct with those two fields. In the QED regime the wait is exponential with rate s*mu - lambda given that it is positive, so P(W > t) = alpha(beta) exp(-(s*mu-lambda) t).

RESULT = QSYS_MMK_QED_STAFFING(…, ‘exact’, true) then walks s up or down until the EXACT Erlang C measure meets the target, starting from the square-root answer. That costs an O(s) Erlang C evaluation per step and is what to use when the answer must be defensible rather than asymptotic.

Returns a struct with fields:

numServers - the recommended s beta - the server slack achieved, (s-a)/sqrt(s) betaTarget - the slack the target asks for, before rounding s up offeredLoad - a = lambda/mu probDelay - the QED delay probability at the recommended s meanWait - the QED mean wait at the recommended s serviceLevel - P(W <= deadline) at the recommended s, for the

‘servicelevel’ criterion

exactUsed - whether the exact Erlang C refinement was applied

Example

res = qsys_mmk_qed_staffing(1000, 1, 0.2); % 20% delayed at most res = qsys_mmk_qed_staffing(1000, 1, struct(‘deadline’,0.02,’level’,0.8), …

‘criterion’,’servicelevel’);

Reference: S. Halfin, W. Whitt (1981). Heavy-traffic limits for queues with many exponential servers. Operations Research 29(3), 567-588. The staffing see also W. Whitt (2007), What you should know about queueing models to set staffing requirements in service systems, Naval Research Logistics 54(5), 476-484.

See also QSYS_MMK_QED, QSYS_MMK_QED_ALPHA, QSYS_MMK.

qsys_mmk_qed_alpha(beta)

ALPHA = QSYS_MMK_QED_ALPHA(BETA)

The Halfin-Whitt delay-probability function

alpha(beta) = [ 1 + beta Phi(beta)/phi(beta) ]^(-1), beta > 0,

with phi and Phi the standard normal density and cdf. It is the limit of the Erlang C delay probability of the M/M/s queue as s -> Inf with beta = (1-rho)sqrt(s) held fixed, and it decreases strictly from 1 at beta = 0 to 0 as beta -> Inf, which is what makes it invertible for staffing.

BETA may be an array. Non-positive entries return 1: with no server slack every arrival is delayed.

Evaluated as phi/(phi + beta*Phi) rather than as the reciprocal of 1 + beta*Phi/phi. The two are the same function, but the quotient Phi/phi overflows once phi underflows (beta beyond about 38), whereas this form degrades to 0/(0+beta) = 0, which is the correct limit.

Reference: S. Halfin, W. Whitt (1981). Heavy-traffic limits for queues with many exponential servers. Operations Research 29(3), 567-588.

See also QSYS_MMK_QED, QSYS_MMK_QED_STAFFING.

qsys_mmk_qed(lambda, mu, s)

QSYS_MMK_QED Halfin-Whitt QED approximation for the M/M/s queue.

RESULT = QSYS_MMK_QED(LAMBDA, MU, S) returns the quality-and-efficiency-driven (QED) approximation for the M/M/S queue with arrival rate LAMBDA and service rate MU per server.

THE REGIME. Let S grow with the offered load a = LAMBDA/MU so that the SERVER SLACK stays of order sqrt(S), i.e.

beta = (1 - rho) sqrt(S) = (S - a)/sqrt(S)

is held fixed. Halfin and Whitt proved that the delay probability then has the non-degenerate limit

alpha(beta) = [ 1 + beta Phi(beta)/phi(beta) ]^(-1)

with phi and Phi the standard normal density and cdf. That is the whole point of the regime: servers are busy a fraction 1 - beta/sqrt(S) of the time, so efficiency tends to 1, and yet the delay probability tends to a constant strictly between 0 and 1, so quality does not collapse. Neither the underloaded regime (alpha -> 0) nor the overloaded one (alpha -> 1) has that property.

WHY USE IT WHEN M/M/s IS EXACT. Erlang C needs a sum of S terms of the form a^j/j!, which overflows in double precision well before the thousands of servers a large contact centre or a datacentre thread pool has; alpha(beta) is three transcendental calls at any S. The approximation is also the object that the staffing rule inverts, see QSYS_MMK_QED_STAFFING.

Returns a struct with fields:

offeredLoad - a = lambda/mu, in erlangs trafficIntensity - rho = a/s beta - the QED server-slack parameter (s-a)/sqrt(s) probDelay - alpha(beta), the probability an arrival waits meanWaitDelayed - E[W | W > 0] = 1/(s*mu - lambda), exact for M/M/s meanWait - E[W] = alpha(beta)/(s*mu - lambda) meanQueueLength - E[Q] = lambda E[W], customers waiting meanNumber - E[N] = a + E[Q] utilization - rho

ACCURACY. The error is O(1/sqrt(s)) at fixed beta: at s = 100 and beta = 0.5 the exact Erlang C is 0.5065 against 0.5045 here, and at s = 100000 it is 0.50461 against 0.50454.

An overloaded model, beta <= 0, has no QED limit; probDelay is then reported as 1 and the waiting-time fields as Inf, which is what the M/M/s queue does.

Example

res = qsys_mmk_qed(990, 1, 1000); % 1000 servers, 99% loaded res.probDelay

Reference: S. Halfin, W. Whitt (1981). Heavy-traffic limits for queues with many exponential servers. Operations Research 29(3), 567-588.

See also QSYS_MMK, QSYS_MMK_QED_STAFFING, QSYS_ERLANGA.

qsys_mgisrgi_whitt(lambda, mu, s, r, patience, varargin)

QSYS_MGISRGI_WHITT Engineering solution of the M/GI/s/r+GI queue.

RESULT = QSYS_MGISRGI_WHITT(LAMBDA, MU, S, R, PATIENCE) computes every standard steady-state measure of the call-center model M/GI/s/r+GI: Poisson arrivals at rate LAMBDA, iid general service times of mean 1/MU, S servers, R extra waiting spaces and iid patience times with a general distribution.

THE TWO APPROXIMATIONS. The general patience law is replaced by STATE- DEPENDENT Markovian abandonment: a customer who is jth from the end of a queue abandons at rate delta_j = h(j/lambda), where h = f/(1-F) is the patience hazard rate, because a customer in that position has been waiting for about j/lambda (eq. 3.3). The total abandonment rate with k waiting is then Delta_k = sum_{j<=k} delta_j (eq. 3.4). The general service law is replaced by an exponential of the same mean (Section 5), which is accurate here because with many servers and non-negligible abandonment the model behaves like a loss system, where the service law is insensitive beyond its mean. What is left is the Markovian M/M/s/r+M(n) model, solved as a birth-and-death process.

WHAT THE PATIENCE LAW CONTRIBUTES. Only the hazard function NEAR THE ORIGIN, not the mean and not the tail: waits are O(1/sqrt(s)) in the many-server regime, so a customer either abandons early or never. That is the paper’s main modelling insight and the reason a two-moment patience fit is not enough.

PATIENCE accepts three forms:

scalar THETA - exponential patience of rate THETA, h(t) = THETA. The

approximations are then EXACT and the model is Erlang A (see QSYS_ERLANGA).

function handle - the hazard rate h(t), used as in eq. (3.3). struct(‘ccdf’,G) - the complementary cdf G(t) = 1-F(t), used through the

integrated form Delta_k = -log G(k/lambda) of eq. (3.6), which is the variant to use when the density is not smooth.

RESULT = QSYS_MGISRGI_WHITT(…, ‘wPoints’, T) also returns the waiting-time cdfs at the times T, obtained by numerically inverting the transforms of eqs. (7.22)-(7.23) and (7.32)-(7.33) with the Abate-Whitt EULER algorithm. Other options: ‘maxQueue’ (truncation level used when R is Inf, default 100000), ‘tol’ (relative tail tolerance for that truncation, default 1e-14), ‘invMethod’ and ‘invN’ (inversion method and node count, default ‘euler’, 41).

Returns a struct with fields:

queueLengthDist - P(N = k) for k = 0..s+r, N the number in system probLoss - P(an arrival is blocked) = p_{s+r}, zero when R is Inf probNoWait - P(W = 0) among entering customers probServed - P(S), an entering customer is eventually served probAbandon - P(A) = 1 - P(S) meanNumber - E[N] varNumber - Var[N] meanQueueLength - E[Q] varQueueLength - Var[Q], Q = (N-s)^+ meanWaitServed - E[W|S] varWaitServed - Var[W|S] meanWaitAbandon - E[W|A] varWaitAbandon - Var[W|A] meanWait - E[W] over entering customers, zeros included secondMomentWait - E[W^2] over entering customers utilization - E[min(N,s)]/s, the fraction of servers busy throughput - rate of served customers, lambda(1-P_loss)P(S) abandonRate - rate of abandoning customers, lambda(1-P_loss)P(A) abandonRates - delta_j, j = 1..r totalAbandonRates - Delta_k, k = 0..r waitPoints - the times T, when requested cdfWaitServed - P(W <= t | S) cdfWaitAbandon - P(W <= t | A) cdfWait - P(W <= t) over entering customers

ACCURACY. Exact for M/M/s/r+M. Elsewhere the paper reports errors of a few percent against simulation, degrading as the service SCV moves away from 1.

Example

% M/M/100/200+M with mean patience 1, from Table 1 of the paper res = qsys_mgisrgi_whitt(102, 1/10, 100, 200, 1);

Reference: W. Whitt (2005). Engineering solution of a basic call-center model. Management Science 51(2), 221-235.

See also QSYS_ERLANGA, QSYS_GGSGI_FLUID, LAPLACE_INVERT_EULER.

qsys_maxima_twomoment(n, mean_, cs2, varargin)

QSYS_MAXIMA_TWOMOMENT Two-moment approximation for the maximum of n variables.

RESULT = QSYS_MAXIMA_TWOMOMENT(N, MEAN, CS2) approximates the mean of the maximum of N iid non-negative variables with the given mean and squared coefficient of variation.

THE SHAPE OF THE ANSWER. For a law with an exponential-like tail the maximum of N samples grows like c~^2 (log N + …): doubling N adds a constant, it does not scale the answer. What the two moments buy is the SLOPE c~^2 of that logarithm and an offset eta:

x_n(q) = c~^2 [log(n eta) - log log(1/q)], E[M_n] = c~^2 [log(n eta) + gamma]

with, for CS2 >= 1, c~^2 = CS2 and eta = (CS2+1)/(2 CS2^2) from the H2 representative, and for CS2 < 1 the shifted-exponential representative c~^2 = sqrt(CS2), eta = exp((1-sqrt(CS2))/sqrt(CS2)).

WHEN NOT TO USE IT. The extreme-value form needs N past a threshold n* ~ CS2/q, because with a highly variable law most of the N samples come from the short component and only about N p of them can contend for the maximum. Measured against exact maxima, the closed form is within a few percent for N >= 100 at CS2 = 4 and 16, and useless at N = 10 for CS2 = 16 – which is exactly what n* predicts. RESULT.reliable reports the test.

AND WHEN TWO MOMENTS ARE NOT ENOUGH. Below CS2 = 1 the maximum is genuinely family-dependent: an Erlang and a shifted exponential with the same two moments have maxima that differ by tens of percent and diverge as N grows, because their tails decay at different rates. Measured on Erlang-4, both the closed form and the fitted shifted exponential are 15-23% high, and they agree with each other, so the gap is the model’s, not the arithmetic’s.

Options: ‘q’ (a quantile level in (0,1); the mean is returned without it), ‘exactFitted’ (also compute the maximum exactly from the fitted representative by integrating 1-F^n, which is the paper’s other recommendation; default true).

Returns a struct with fields value, slope, eta, threshold, reliable, family and exactFittedValue.

Reference: C. Crow, D. Goldberg, W. Whitt (2007). Two-moment approximations for maxima. Operations Research 55(3), 532-548.

See also FJ_RMAX, FJ_XMAX_ERLANG.

qsys_gtmtst_fluid(lambdaFun, sFun, muFun, patienceCcdf, T, varargin)

QSYS_GTMTST_FLUID The Gt/Mt/st+GI many-server fluid queue.

RESULT = QSYS_GTMTST_FLUID(LAMBDAFUN, SFUN, MUFUN, PATIENCECCDF, T) solves the time-varying many-server fluid queue on [0,T]: arrival rate LAMBDAFUN(t), staffing SFUN(t), exponential service at rate MUFUN(t), general patience with complementary cdf PATIENCECCDF, unlimited waiting room.

THE MODEL ALTERNATES BETWEEN TWO REGIMES, and the whole algorithm is the bookkeeping of that alternation:

UNDERLOADED the queue is empty and every arrival enters service at once, so

the system is the infinite-server fluid model and B obeys B’(t) = lambda(t) - mu(t)B(t) (eq. 18 in its Mt form). It ends when B reaches s while lambda exceeds the rate Gamma(t) = s’(t) + s(t)mu(t) at which capacity frees up (15).

OVERLOADED every server is busy, B(t) = s(t), fluid enters service at

exactly Gamma(t), and the queue is described by its BOUNDARY WAITING TIME w(t), the age of the oldest fluid still waiting. Content of age x is what arrived x ago and has not abandoned, q(t,x) = lambda(t-x)F^c(x), and the boundary moves by the delay differential equation (21)

w’(t) = 1 - Gamma(t) / [lambda(t-w(t)) F^c(w(t))].

It ends when w returns to 0 with lambda no longer above Gamma (14).

WHY w AND NOT Q. The queue content is a functional of w, but not the other way round: two systems with the same Q and different age profiles abandon at different rates. Tracking the boundary keeps the age profile exact, which is what makes a general patience law admissible at all.

Options:

‘dt’, DT - grid step, default T/2000 ‘B0’, B - fluid in service at time 0, default 0 ‘w0’, W - boundary waiting time at time 0, default 0 ‘sPrime’, SP - s’(t); differentiated from SFUN numerically when absent ‘patiencePdf’, F - the patience density, for the abandonment rate;

differenced from the ccdf when absent

‘lambdaPast’, L - the arrival rate before time 0, needed only when the

queue starts non-empty

Returns a struct on the grid with fields times, regime (1 overloaded, 0 underloaded), B, Q, X = B+Q, w, v (potential waiting time), sigma (service completion rate), alpha (abandonment rate), utilization, arrivalRate, staffing and capacityRate (Gamma).

Example

fc = @(x) exp(-0.5*x); res = qsys_gtmtst_fluid(@(t) 100+30*sin(t), 100, 1, fc, 40);

Reference: Y. Liu, W. Whitt (2012). The Gt/GI/st+GI many-server fluid queue. Queueing Systems 71, 405-444; Y. Liu, W. Whitt (2014). Algorithms for time-varying networks of many-server fluid queues. INFORMS J. on Computing 26(1), 59-73.

See also QSYS_GGISGI_FLUID, QSYS_MTGINF, NPFQN_GTMTST_FLUID.

qsys_gig1_bnds_extremal(lambda, mu, ca, cs, varargin)

QSYS_GIG1_BNDS_EXTREMAL Extremal two-moment bounds for the GI/GI/1 queue.

RESULT = QSYS_GIG1_BNDS_EXTREMAL(LAMBDA, MU, CA, CS) returns the tightest known interval for the mean steady-state waiting time of a GI/GI/1 queue known ONLY through the first two moments of its interarrival and service times, together with the classical bounds it improves on.

WHAT THE INTERVAL MEANS. Two moments do not determine E[W]; they determine a SET of possible values, and the width of that set is the honest uncertainty in any two-moment approximation. The extremal distributions are the ones that attain its ends:

• the lower end is attained by deterministic interarrival times and a three-point service law concentrated on multiples of that interval, and has the closed form rho((1+cs^2)rho - 1)^+ / (2(1-rho)) (eq. 2.12);

• the upper end is attained (asymptotically) by TWO-POINT laws: an interarrival law with an atom at 0, and a service law whose upper atom runs off to infinity while its probability vanishes. Making an interarrival time larger only empties the queue once, but making a service time larger delays every customer behind it, which is why the two ends look so different.

HOW THE UPPER END IS COMPUTED. Chen and Whitt reduce that limit to a D(1/p)/RS(D(rho),p)/1 model with p = 1/(1+ca^2) and RS a geometric random sum, then evaluate its mean waiting time by Spitzer’s identity with the negative binomial pmf (their Algorithm 1). The sum is truncated in both indices, so this bound is a numerical limit, not a formula. The closed-form companion (eq. 3.4) uses the D/M/1 root delta = exp(-(1-delta)/rho) and is within about 1% of it.

Options:

‘K’, K - truncation of the negative binomial value, default 4000 ‘N’, N - truncation of the random-walk length, default 2000 ‘skipTight’, TF - skip the O(K*N) tight bound and return the closed forms

only, default false

Returns a struct whose waiting-time fields are TIMES IN QUEUE (add 1/MU for a response time), all for the given LAMBDA:

trafficIntensity - rho = lambda/mu lowerBound - the tight lower bound, eq. (2.12) upperBound - the conjectured tight upper bound, eq. (3.2) by Algorithm 1 upperBoundClosed - the closed-form upper bound, eq. (3.4) upperBoundDaley - Daley’s bound, eq. (2.7) upperBoundKingman- Kingman’s bound, eq. (2.6) heavyTraffic - the heavy-traffic approximation, eq. (2.9) delta - the D/M/1 root behind upperBoundClosed relativeWidth - (upperBound-lowerBound)/upperBound, the fraction of the

answer that two moments genuinely leave undetermined

Example

res = qsys_gig1_bnds_extremal(0.5, 1, 2, 2); % rho = 0.5, ca^2 = cs^2 = 4 [res.lowerBound res.upperBound res.upperBoundKingman] % 0.750 3.470 5.000

Reference: Y. Chen, W. Whitt (2020). Algorithms for the upper bound mean waiting time in the GI/GI/1 queue. Queueing Systems 94, 327-356. The lower bound is classical, restated there as eq. (2.12); Kingman (1962) and Daley (1977) are the two established upper bounds.

See also QSYS_GIG1_LBND, QSYS_GIG1_UBND_KINGMAN, QSYS_GIG1_RQ.

qsys_ggnm_diffusion(lambda, mu, n, m, ca, cs, varargin)

QSYS_GGNM_DIFFUSION Diffusion approximation for the G/GI/n/m queue.

RESULT = QSYS_GGNM_DIFFUSION(LAMBDA, MU, N, M, CA, CS) approximates the steady state of a queue with a general arrival process of rate LAMBDA and variability CA, iid general service of mean 1/MU and variability CS, N servers and M extra waiting spaces (M = Inf for an unbounded queue).

THE APPROXIMATION IS ONE DIFFUSION WITH TWO REGIONS. Below the staffing level the queue behaves like an infinite-server system, whose limit is NORMAL with variance-to-mean ratio the ASYMPTOTIC PEAKEDNESS

z = 1 + (ca^2 - 1) omega_G, omega_G = int G^c(x)^2 dx / int G^c(x) dx

(eqs. 1.6-1.7); above it the queue behaves like a single-server queue, whose limit is EXPONENTIAL with variability v = (ca^2 + cs^2)/2 (eq. 3.7). The steady-state law is a normal piece spliced to an exponential piece and every measure below is an integral of that density (eq. 3.14).

WHAT z SAYS. The service-time distribution enters the delay probability ONLY through omega_G, which is 1 for deterministic service, 1/2 for exponential, and falls toward 0 as service gets more variable. So at ca^2 = 1 the delay probability does not depend on the service law at all (z = 1), which is the long-standing M/GI/n-by-M/M/n approximation; away from ca^2 = 1 it does, and this is how much.

With M = Inf the delay probability is alpha(beta/sqrt(z)) for the Halfin-Whitt function alpha (eq. 3.10), so this generalizes QSYS_MMK_QED.

Options:

‘serviceCcdf’, G - G^c(x) = P(S > x); the exponential of rate MU by default ‘tol’, TOL - service-tail cut for the peakedness integral, 1e-12 ‘panels’, P - Simpson panels for it, default 4000

Returns a struct with fields beta (the QED server slack), gamma (the scaled waiting room), peakedness (z), peakednessWeight (omega_G), variability (v), probDelay, probBlock, meanQueueLength, meanNumber, meanWait, utilization, throughput and trafficIntensity.

Example

res = qsys_ggnm_diffusion(95, 1, 100, 20, 1.5, 0.8);

Reference: W. Whitt (2004). A diffusion approximation for the G/GI/n/m queue. Operations Research 52(6), 922-941.

See also QSYS_MMK_QED, QSYS_GIGK_APPROX_WHITT, QSYS_ERLANGA.

qsys_ggisgi_fluid(lambda, mu, s, patienceCcdf, varargin)

QSYS_GGISGI_FLUID Steady state of the G/GI/s+GI fluid model.

RESULT = QSYS_GGISGI_FLUID(LAMBDA, MU, S, PATIENCECCDF) computes the unique steady state of the deterministic fluid model of a multiserver queue with customer abandonment: arrival rate LAMBDA, S servers each of rate MU, general service times and general patience times whose complementary cdf is the function handle PATIENCECCDF, F^c(t) = P(patience > t).

THE MODEL. Scale the content by S and let S grow. Customers become quanta of fluid but their sojourns do not shrink, so the ages survive the limit: the state is the density b(x) of fluid that has been IN SERVICE for time x and the density q(x) of fluid that has been WAITING for time x. With rho = lambda/(s*mu),

rho <= 1 b(x) = rho G^c(x), q = 0, no abandonment, no wait rho > 1 b(x) = G^c(x), q(x) = rho F^c(x) on [0,w] and 0 beyond

where the queue boundary w solves F^c(w) = 1/rho (eq. 3.6). That one equation carries the whole overloaded regime: fluid that survives its patience for w enters service, so the surviving fraction F^c(w) must equal the fraction 1/rho the servers can absorb.

WHAT THE DISTRIBUTIONS CONTRIBUTE (Corollary 3.1). The rates and the number in service depend on G and F only through their means. The wait w, the queue content and its age profile depend on F BEYOND its mean but on G only through its mean. The age profile in service depends on G beyond its mean. Neither the number of servers nor anything about the arrival process beyond its rate appears at all, which is why this model says nothing about the QED regime and everything about the overloaded one.

Options:

‘servingCcdf’, G - the service-time ccdf, needed only for the in-service

age density; defaults to exponential of rate MU

‘agePoints’, X - ages at which to return the two densities ‘tol’, T - bisection tolerance for w, default 1e-12 ‘maxTime’, TMAX - largest age searched for w, default grows automatically

Returns a struct with fields:

regime - ‘underloaded’, ‘balanced’ or ‘overloaded’ trafficIntensity - rho = lambda/(s*mu) offeredWait - w, the wait of every customer who is served, 0 unless overloaded meanWait - E[W] over all customers, = int_0^w F^c(t)dt = m_a F_e(w) meanWaitServed - w again, the fluid wait being deterministic meanWaitAbandon - E[patience | patience <= w] probAbandon - 1 - 1/rho when overloaded, 0 otherwise meanQueueLength - Q = lambda * meanWait, in customers meanNumberInService - B = min(lambda/mu, s), in customers meanNumber - B + Q utilization - min(rho,1) throughput - min(lambda, s*mu) abandonRate - lambda - throughput agePoints, serviceAgeDensity, queueAgeDensity - the densities per server at

the requested ages, b(x) and q(x)

ACCURACY. This is the s -> Inf limit, so it is an approximation at finite S that improves with S and with the overload. Whitt (2004) reports it as crude at s = 100, rho = 1.02 and good at s = 100, rho = 1.10.

Example

% 100 agents, 10% overload, exponential patience of mean 5 res = qsys_ggisgi_fluid(110, 1, 100, @(t) exp(-t/5));

Reference: W. Whitt (2006). Fluid models for multiserver queues with abandonments. Operations Research 54(1), 37-54, Theorem 3.1 and Corollary 3.2.

See also QSYS_ERLANGA, QSYS_MGISRGI_WHITT.

qsys_ggingi_tga(lambda, mu, n, ca, cs, patienceCcdf, varargin)

QSYS_GGINGI_TGA Truncated Gaussian approximation for the G/GI/n+GI queue.

RESULT = QSYS_GGINGI_TGA(LAMBDA, MU, N, CA, CS, PATIENCECCDF) approximates the steady state of a heavily-loaded multiserver queue with abandonment: a general stationary arrival process of rate LAMBDA and variability CA, iid general service of mean 1/MU and variability CS, N servers, unlimited waiting room and iid general patience with complementary cdf PATIENCECCDF.

THE APPROXIMATION IS A FLUID CENTRE PLUS A GAUSSIAN FLUCTUATION, TRUNCATED. In the efficiency-driven regime (rho > 1 held fixed as N grows) the fluid limit gives the centre – all servers busy, waiting time w = F^-1(1-1/rho), queue Q = lambda int_0^w F^c – and the many-server central limit theorem gives a NORMAL fluctuation of order sqrt(N) around it. Adding the two directly can produce negative queues and negative waits, so both are TRUNCATED at zero, which is what makes the formulas usable down to moderate overload; the paper reports good accuracy for rho > 1.02 and abandonment rates below 2.

Three independent sources of variability enter separately, which is what lets the exponential-service formula be generalized: the service law appears only as the factor (cs+1)rho in sigma_W^2 (eq. 24), which is 2rho at cs = 1.

An UNDERLOADED model (rho <= 1) has no queue in the limit; the number in system is then normal with the infinite-server variance, whose variance-to-mean ratio is the asymptotic peakedness of QSYS_GGNM_DIFFUSION.

Options:

‘patiencePdf’, F - the patience density; differenced from the ccdf when absent ‘serviceCcdf’, G - G^c(x) = P(S > x), used only in the underloaded branch

Returns a struct with fields regime, trafficIntensity, fluidWait, fluidQueueLength, meanWait, varWait, meanQueueLength, varQueueLength, meanNumberInService, meanNumber, probDelay, probAbandon, sigmaW and sigmaX.

Example

fc = @(x) exp(-0.5*x); fp = @(x) 0.5*exp(-0.5*x); res = qsys_ggingi_tga(120, 1, 100, 1, 1, fc, ‘patiencePdf’, fp);

Reference: Y. Liu, W. Whitt, Y. Yu (2016). Approximations for heavily-loaded G/GI/n+GI queues. Naval Research Logistics 63(3), 187-217.

See also QSYS_GGISGI_FLUID, QSYS_ERLANGA, QSYS_GGNM_DIFFUSION.

qsys_erlanga(lambda, mu, theta, s, r, varargin)

QSYS_ERLANGA Exact analysis of the Erlang A model M/M/s/r+M.

RESULT = QSYS_ERLANGA(LAMBDA, MU, THETA, S) analyzes the Erlang A queue: Poisson arrivals at rate LAMBDA, exponential service of rate MU at each of S servers, and exponential patience of rate THETA, so a waiting customer abandons after an exponential time of mean 1/THETA. The waiting room is infinite; THETA > 0 makes the model ergodic at every load, including LAMBDA above S*MU.

RESULT = QSYS_ERLANGA(LAMBDA, MU, THETA, S, R) allows only R extra waiting spaces, so an arrival finding S+R customers is blocked and lost. R = Inf is the default.

The number in system is the birth-and-death process with birth rate LAMBDA and death rate min(k,S)*MU + (k-S)^+ *THETA, so every measure below is EXACT: this is the special case in which the state-dependent Markovian approximation of QSYS_MGISRGI_WHITT reproduces the model rather than approximating it (eq. 7.12 of the reference). THETA = 0 recovers M/M/s/r without abandonment, and then a finite R is required whenever LAMBDA >= S*MU.

Accepts and returns exactly what QSYS_MGISRGI_WHITT does, including the ‘wPoints’ option for the waiting-time cdfs.

Example

res = qsys_erlanga(102, 1/10, 1/10, 100); res.probAbandon % fraction of arrivals that abandon

Reference: W. Whitt (2005). Engineering solution of a basic call-center model. Management Science 51(2), 221-235, Section 7 and eq. (7.12). The model itself is due to C. Palm (1937, 1957).

See also QSYS_MGISRGI_WHITT, QSYS_MMK, QSYS_MMCK.

qsys_mmapgk1(MMAP, svc, varargin)

QSYS_MMAPGK1 Per-type waiting times of the MMAP[K]/G[K]/1 FCFS queue.

RESULT = QSYS_MMAPGK1(MMAP, SVC) analyzes a single-server FCFS queue with:

MMAP - cell {D0, D1, …, DK}: a marked Markovian arrival process of order

ma whose type-k arrivals carry the block Dk

SVC - cell of K service laws, one per type. Each entry is a LINE

Distribution object, a cell {D0s, D1s} holding a phase-type representation, or a struct with fields

.lst function handle s -> E[exp(-s S)], admitting COMPLEX s .moments raw moments [E[S], E[S^2], …], at least numWMoms+1 of

them

The struct form is what a solver hands in: sn.lst already carries the transform of the ORIGINAL law, while sn.proc carries only its phase-type fit, and it is the original that this analysis needs. The service laws may be GENERAL and need not share a family across types.

RESULT = QSYS_MMAPGK1(…, ‘wPoints’, T) also evaluates the per-type waiting time CDF at the times T by numerical transform inversion. RESULT = QSYS_MMAPGK1(…, ‘numWMoms’, N) returns N moments per type (default 3). RESULT = QSYS_MMAPGK1(…, ‘quadNodes’, Q) sets the quadrature order used for a service law with no closed-form matrix transform (default 120).

THE METHOD, which is He’s, theorem for theorem. FCFS makes the actual waiting time of a customer the WORKLOAD it finds on arrival, so everything follows from the joint transform of workload and arrival phase, f(s)_j = E[exp(-s V) 1{phase = j}], which by He’s Theorem 4.1 (eq. 4.6) satisfies

f(s) [ s I + D0 + sum_k Dk gk(s) ] = s v0, (*)

with v0 the idle-phase vector, his y0. The unknown v0 needs NO search for the roots of the determinant: the matrix U solving

U = D0 + sum_k Dk Fk(U), Fk(U) = int_0^inf exp(U t) dFk(t),

is his eq. (4.4), the generator of the underlying Markov process obtained by EXCISING the busy periods, and eq. (4.5) with Theorem 4.2 give y0 Q = 0 and y0 e = 1 - rho, i.e.

v0 = (1 - rho) pi_U.

The same vector is what the analyticity of (*) forces: for every left eigenpair (w, u) of U one has w [D0 + sum_k Dk gk(-u) + (-u) I] = 0, so the roots of (*) in the closed right half plane are exactly s = -u over the spectrum of U, and imposing v0 r_i = 0 at each right null vector reproduces the stationary vector to 2.5e-13. The stationary route is the one taken, as it needs no complex eigenvector and no rule for telling the structural root at the origin from a genuine one.

The per-type actual waiting time is the workload seen by a type-k arrival, biased by that type’s own arrival block, which is his Theorem 5.1 eq. (5.1) summed over the post-arrival phase:

E[exp(-s Wk)] = f(s) Dk e / lambda_k.

SCOPE. He allows an arrival to be a BATCH carrying a sequence of types, and his Theorem 5.3 then multiplies the transform by prod_{i<n} f*_{h_i}(s), the service of the customers ahead of the tagged one WITHIN its own batch. This function covers the single-customer-per-arrival case, his Special case 3.3, where that product is empty – which is exactly the MMAP convention LINE carries, {D0, D1, D^(1), …, D^(K)} with one customer per epoch.

Returns a struct with fields:

lambda - per-type arrival rates (1 x K) utilization - rho = sum_k lambda_k E[S_k] idleVector - v0, the idle-phase vector, summing to 1 - rho waitLST - function handle s -> row vector of E[exp(-s Wk)] waitMoments - (K x numWMoms) per-type waiting time moments meanWaitingTime - per-type mean waiting time (1 x K) meanSojournTime - per-type mean sojourn time (1 x K) meanQueueLength - E[N], by Little over all types waitCDF - (K x numel(wPoints)) per-type waiting time CDF waitPoints - the requested points analyzer - name of the analyzer used

Reference:

Qi-Ming He, “The versatility of MMAP[K] and the MMAP[K]/G[K]/1 queue”, Queueing Systems 38(4):397-418, 2001.

See also qsys_mapg1(), qsys_mapphc(), MMAPPH1FCFS

qsys_mapphc(D0, D1, alpha, S, c, varargin)

QSYS_MAPPHC Analyzes a MAP/PH/c FCFS queue exactly.

RESULT = QSYS_MAPPHC(D0, D1, ALPHA, S, C) analyzes a c-server FCFS queue with:

D0, D1 - arrival MAP of order ma ALPHA - PH service initial vector (1 x ms) S - PH service sub-generator (ms x ms) C - number of servers (identical, so the service law is shared)

RESULT = QSYS_MAPPHC(…, ‘maxNumComp’, N) caps the number of queue length probabilities returned (default 500). RESULT = QSYS_MAPPHC(…, ‘numWMoms’, K) returns K waiting-time moments (default 3). RESULT = QSYS_MAPPHC(…, ‘wPoints’, T) evaluates the waiting-time CCDF at the times T.

THE STATE SPACE. With c identical servers the server identities carry no information, so the service phases are held as a MULTISET: a configuration is n = (n_1..n_ms) with sum(n) = k servers busy in phase i. There are nchoosek(ms+k-1,k) of them, the count of Asmussen and Moller (2001), against ms^k for the ordered space. Levels 0..c-1 are the boundary (level = servers busy), levels >= c repeat and carry the queue.

THE WAITING TIME. An arrival that finds j customers waiting ahead of it waits for j+1 service completions, so Wq is the (j+1)-st event time of the configuration MAP (Lc, Cdep) started at the arrival-epoch configuration. Folding the matrix-geometric level distribution over j gives

G’(t) = G(t) Lj + R G(t) Cj, G(0) = (I-R)^-1 kron(D1,I) / lambda, P(Wq > t) = pi_c G(t) e,

a LINEAR matrix ODE, so Wq is matrix-exponential. Its transform obeys the generalized Sylvester equation g(sI-Lj) - R g Cj = G(0), and every moment reuses that one operator with a different right-hand side.

Returns a struct with fields:

meanQueueLength - E[N], number in system meanWaitingTime - E[Wq], time in queue meanSojournTime - E[W] = E[Wq] + E[service] utilization - rho = lambda E[service] / c, per server queueLengthDist - P(N = n), n = 0, 1, … waitingTimeMoments - E[Wq^k], k = 1..numWMoms waitingTimeCCDF - P(Wq > t) at the requested wPoints probWait - P(Wq > 0), an arrival finds every server busy phaseCount - nchoosek(ms+c-1,c), the repeating configuration count analyzer - name of the analyzer used

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. D.P. Gaver, P.A. Jacobs, G. Latouche, “Finite birth-and-death models in randomly changing environments”, Adv. Appl. Probab. 16:715-731, 1984.

See also qsys_mapmc(), qsys_mapph1(), qsys_phmc(), ph_multisets()

mexify_qsys

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

This script generates MEX (MATLAB Executable) versions of qsys_ functions for improved performance. All functions in this module are pure numerical computations (scalar or vector input/output).

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

qsys_tandem_ub_ciucu(x, lst, p, mu, dlst)

QSYS_TANDEM_UB_CIUCU Tail bounds for a GI/Hn/1 -> ./Hn/1 tandem.

RESULT = QSYS_TANDEM_UB_CIUCU(X, LST, P, MU) returns polynomial-exponential upper bounds on the tails of the end-to-end waiting time W and sojourn time S of a tandem of two single-server FCFS stations, fed by a renewal arrival process and serving the same hyperexponential law at both stations. X is the vector of thresholds at which the tails are bounded, LST is a function handle evaluating the interarrival Laplace-Stieltjes transform E[e^{-s X}] for s >= 0, and P, MU are the phase probabilities and rates of the service law Y, Z ~ sum_i P(i) Exp(MU(i)); a scalar P = 1 gives exponential service.

RESULT = QSYS_TANDEM_UB_CIUCU(X, LST, P, MU, DLST) also takes a handle evaluating E[X e^{-s X}], i.e. minus the derivative of LST. Without it that derivative is obtained by a Richardson-extrapolated central difference, which costs about four extra LST evaluations and loses roughly four digits.

The bounds are those of Ciucu and Mehri: with theta the positive root of E[e^{theta (Y-X)}] = 1 (Lundberg/Kingman exponent of the first station, which by identical service is also the exponent of the tandem) and alpha = E[X e^{-theta X}], the test function

gamma(u,v) = 1{0<=u<=v} [1 - A e^{-theta u} - (B + C u + D v) e^{-theta v}]

is made to satisfy the integral inequality of their Theorem 1(b) by the five sufficient conditions of their Lemma 4, which fix

A = 1, C = theta sum_i p_i/(mu_i-theta) / sum_i p_i mu_i/(mu_i-theta)^2, D = (-C E[U e^{theta V}]/E[V e^{theta V}]) v 0, U = Y-X, V = Z-X, B = C (1/mu_1 - alpha E[e^{theta Z}]) if D = 0,

= (C+D)/(mu_1-theta) - theta/mu_1 if D > 0,

with mu_1 the smallest service rate. Their Corollary 2 then gives

P(S > x) <= sum_i p_i { e^{-mu_i x}

• mu_i/(mu_i-theta) (A+B) (e^{-theta x} - e^{-mu_i x})

• mu_i/(mu_i-theta)^2 (C+D) (((mu_i-theta)x-1) e^{-theta x}

  • e^{-mu_i x}) }

and the corresponding closed form for W when the service is exponential. E[V e^{theta V}] is positive at any stable load, so D is always well defined: h(s) = E[e^{s(Z-X)}] is convex with h(0) = h(theta) = 1, hence h’(theta) > 0.

The two exponentials mix a polynomial of degree one in x, which is what lets the bound follow the concave bend of the tail on a linear-log scale where a purely exponential bound cannot. In the M/M/1 -> ./M/1 case the five inequalities hold as equalities, so gamma is the exact joint distribution and both bounds are exact:

P(W > x) = (1 - 2 theta^2/(mu(mu+theta)) + x(mu-theta)theta/(mu+theta)) e^{-theta x} P(S > x) = (1 + theta x) e^{-theta x} .

Away from it the bound stays sharp: against an exact CTMC reference for the Erlang(2)/M/1 -> ./M/1 tandem it is within 2% at P(S>x) = 1e-2 and within 0.6% at 5e-10, with the correct asymptotic slope theta^2/(mu(1-alpha mu)). Accuracy degrades with service variability, to about a factor of two at CV(Y) = 2.

Returns a struct with fields:

S - upper bound on P(S > x), one entry per threshold, capped at 1 W - upper bound on P(W > x), NaN unless the service is exponential theta - the tail decay rate, positive root of E[e^{theta (Y-X)}] = 1 alpha - E[X e^{-theta X}] A,B,C,D- the coefficients of gamma fixed by Lemma 4 analyzer - ‘qsys_tandem_ub_ciucu’

Example, a D/M/1 -> ./M/1 tandem at utilization 3/4 with unit service rate:

r = qsys_tandem_ub_ciucu([5 10], @(s) exp(-s*4/3), 1, 1, @(s) (4/3)*exp(-s*4/3));

Reference: F. Ciucu, S. Mehri, “On the Distribution of Sojourn Times in Tandem Queues”, Proc. ACM Meas. Anal. Comput. Syst. 9(2), Article 27, 2025 (ACM SIGMETRICS 2025). Registered in .citations() as ‘tandemub’.

qsys_mm1_lindley(lambda, mu, Wn, mmax)

QSYS_MM1_LINDLEY Conditional waiting-time moments of the M/M/1 Lindley recursion.

RESULT = QSYS_MM1_LINDLEY(LAMBDA, MU, WN) returns the exact conditional mean and variance of the waiting time of customer n+1 in an FCFS M/M/1 queue with arrival rate LAMBDA and service rate MU, given that customer n waited WN. WN may be a vector, in which case every output is evaluated elementwise.

RESULT = QSYS_MM1_LINDLEY(LAMBDA, MU, WN, MMAX) also returns the conditional raw moments of orders 1 to MMAX.

This is one step of Lindley’s recursion W_{n+1} = max(W_n + S_n - A_n, 0) with A_n ~ Exp(LAMBDA) and S_n ~ Exp(MU). Unlike every other qsys_* function, the quantities here are conditional on the current state rather than stationary, so they are defined and finite for any load, including LAMBDA >= MU.

The m-th conditional moment is

E[W_{n+1}^m | W_n] = lambda mu/(lambda+mu) [ S + T ], S = sum_{k=0}^{m} C(m,k) W_n^k (m-k)! / mu^(m-k+1), T = (-1)^m e^{-lambda W_n} (Gamma(m+1,-lambda W_n) - m!) / lambda^(m+1),

where the density of S_n - A_n is the asymmetric Laplace density lambda mu/(lambda+mu) times e^{-mu x} for x > 0 and e^{lambda x} for x < 0. Because m+1 is a positive integer, the upper incomplete gamma function admits the finite form Gamma(m+1,x) = m! e^{-x} sum_{k=0}^{m} x^k/k!, valid at the negative argument -lambda W_n needed here. Substituting it cancels the growing exponential and leaves the numerically stable

T = (-1)^m m! ( sum_{k=0}^{m} (-lambda W_n)^k/k! - e^{-lambda W_n} ) / lambda^(m+1),

which is what this function evaluates. No incomplete gamma routine is needed.

The mean is returned from the equivalent explicit form

E[W_{n+1} | W_n] = W_n + (lambda-mu)/(lambda mu)

• mu e^{-lambda W_n} / (lambda (lambda+mu)),

and the variance as the second moment less the squared mean.

Returns a struct with fields:

mean - Conditional mean E[W_{n+1} | W_n], same size as WN var - Conditional variance Var[W_{n+1} | W_n], same size as WN moments - numel(WN) x MMAX conditional raw moments, orders 1 to MMAX mmax - Highest moment order computed analyzer - Identifier string

Examples

r = qsys_mm1_lindley(0.8, 1.0, 2.0); r.mean % 1.8902 r.var % 1.7016 qsys_mm1_lindley(0.8, 1, 0).mean % 0.4444, an empty-queue step

Reference: S. Palomo, J. Pender, “Learning the Tandem Network Lindley Recursion”, Proc. Winter Simulation Conference, 2021, theorem 1 and corollary 2. Verified against 4e6 Monte Carlo replications to 5e-4 relative error for m = 1, 2, 3.

See also QSYS_HH1_LINDLEY, QSYS_TANDEM_LINDLEY, QSYS_MM1_TANDEM_LINDLEY, QSYS_MM1

Copyright (c) 2012-2026, Imperial College London All rights reserved.

qsys_gigk_rqt_gamma(rho, mu, Gamma_a, sigma_s, k, alpha_a, regime)

[GAMMA_S,THETA]=QSYS_GIGK_RQT_GAMMA(RHO,MU,GAMMA_A,SIGMA_S,K,ALPHA_A,REGIME)

Service variability parameter of the Robust Queueing Theory (RQT) framework, obtained from the first two moments by the adaptation of Section 7.1,

Gamma_s = (2 (theta0 + theta1 sigma_s^2/k + theta2 Gamma_a^2 rho^2 k))^((a-1)/a)

• Gamma_a k^((a-1)/a),

where (theta0,theta1,theta2) are regressed so that the worst-case system time of Theorem 3 approximates the MEAN system time of the corresponding stochastic queue. The arrival side needs no adaptation: Gamma_a = sigma_a for an external renewal stream, and the network traffic equations (NPFQN_TRAFFIC_RQT) carry it to the internal streams. Since the last term cancels Gamma_a at alpha=2, the adaptation acts on the sum Gamma_a + Gamma_s/k^(1/alpha) that Theorem 3 reads.

THE FACTOR 2 IS NOT IN THE PRINTED FORMULA and is restored here. Section 7.1 states that the functional form is motivated by Kingman’s bound, which the alpha=2 bound of Theorem 3 reproduces when (Gamma_a+Gamma_s)^2 = 2(sigma_a^2 + sigma_s^2); the published (theta0,theta1,theta2) are all near unity, i.e. they are corrections to that bound rather than a substitute for its factor 2. Dropping the factor puts M/M/1 about 40% BELOW its exact mean system time at rho=0.9 and M/M/3 about 45% below, which contradicts the errors of at most 9.5% that Tables 2-3 report; restoring it gives +4.7% and +8.8%, inside that envelope. See _kb/06-solver-catalog.md.

CAUTION: even so the form is not dimensionally homogeneous, since theta0 is an additive constant on a scale of variances, so it is only valid in the time unit the regression was run in. This routine therefore evaluates it in units of the mean service time, 1/mu = 1, and converts the result back. Do not call it with a time scale of your own choosing.

Inputs:

RHO - Traffic intensity lambda/(k*mu) MU - Service rate of each server, which sets the time unit GAMMA_A - Variability parameter of the arrival uncertainty set SIGMA_S - Standard deviation of the service time K - Number of servers (default 1) ALPHA_A - Effective arrival tail coefficient in (1,2] (default 2) REGIME - Adaptation regime, Table 1: ‘independent’ (default, service

distribution unknown), ‘normal’ or ‘pareto’

Returns:

GAMMA_S - Variability parameter of the service uncertainty set THETA - The (theta0,theta1,theta2) triple used

Reference: C. Bandi, D. Bertsimas, N. Youssef (2015). Robust Queueing Theory. Operations Research 63(3), 676-700, Section 7.1 and Table 1.

qsys_gigk_rqt(lambda, mu, Gamma_a, Gamma_s, k, alpha_a, alpha_s)

[W,RHOHAT,SWORST]=QSYS_GIGK_RQT(LAMBDA,MU,GAMMA_A,GAMMA_S,K,ALPHA_A,ALPHA_S)

Robust Queueing Theory (RQT) worst-case system time of a G/G/k FCFS queue. The arrival and service processes are not described by distributions but by the polyhedral uncertainty sets

U^a = { T : (sum_{i=k+1}^n T_i - (n-k)/lambda)/(n-k)^(1/alpha_a) >= -Gamma_a } U^s = { X : (sum_{i=k}^n X_i - (n-k+1)/mu)/(n-k+1)^(1/alpha_s) <= Gamma_s }

whose shape follows the (generalized) central limit theorem: alpha=2 is the finite-variance regime, alpha in (1,2) the heavy-tailed one. The performance analysis is then a worst-case optimization rather than an expectation.

The returned W is the closed-form bound of Theorem 3 (Theorem 8 when the two tail coefficients differ, with alphabar = min(alpha_a,alpha_s)),

W <= (alphabar-1)/alphabar^(alphabar/(alphabar-1))

• lambda^(1/(alphabar-1)) (Gamma_a + Gamma_s/k^(1/alphabar))^(alphabar/(alphabar-1)) / (1-rho)^(1/(alphabar-1)) + k/lambda,

which for k=1 reduces to Theorem 2 and, at alphabar=2, to the Kingman-like form (lambda/4)(Gamma_a+Gamma_s)^2/(1-rho) + 1/lambda. SWORST returns instead the exact worst case over the uncertainty sets, eq. (45), i.e. the supremum over the integer x = nu-j+1 >= 1 of

x/mu + Gamma_s x^(1/alpha_s) - k(x-1)/lambda + Gamma_a (k(x-1))^(1/alpha_a),

a one-dimensional problem. W >= SWORST by construction, and the two agree closely in heavy traffic. The arrival deviation ADDS to the worst case, since the adversary shortens the interarrival times: the sign printed in eq. (12) is easily misread as a subtraction of the whole arrival bracket, and reading it that way puts SWORST an order of magnitude below W.

Note that W is a SYSTEM time (waiting plus service), not a waiting time, and that its additive term is k/lambda rather than the mean service time 1/mu.

Inputs:

LAMBDA - Arrival rate MU - Service rate of each server GAMMA_A - Variability parameter of the arrival uncertainty set GAMMA_S - Variability parameter of the service uncertainty set K - Number of servers (default 1) ALPHA_A - Arrival tail coefficient in (1,2] (default 2) ALPHA_S - Service tail coefficient in (1,2] (default 2)

Returns:

W - Closed-form bound on the system time (Theorem 3 / Theorem 8) RHOHAT - Modified utilization (so that M/M/1 formulas still hold) SWORST - Exact worst-case system time over the uncertainty sets, eq. (45)

Reference: C. Bandi, D. Bertsimas, N. Youssef (2015). Robust Queueing Theory. Operations Research 63(3), 676-700.

qsys_gig1_rqt(lambda, mu, Gamma_a, Gamma_s, alpha_a, alpha_s)

[W,RHOHAT,SWORST]=QSYS_GIG1_RQT(LAMBDA,MU,GAMMA_A,GAMMA_S,ALPHA_A,ALPHA_S)

Robust Queueing Theory (RQT) worst-case system time of a G/G/1 FCFS queue, the single-server case of QSYS_GIGK_RQT. The closed-form bound is Theorem 2,

W <= (alpha-1)/alpha^(alpha/(alpha-1)) * lambda^(1/(alpha-1))

(Gamma_a+Gamma_s)^(alpha/(alpha-1)) / (1-rho)^(1/(alpha-1)) + 1/lambda,

and SWORST is the exact worst case over the uncertainty sets, eq. (12).

Inputs:

LAMBDA - Arrival rate MU - Service rate GAMMA_A - Variability parameter of the arrival uncertainty set GAMMA_S - Variability parameter of the service uncertainty set ALPHA_A - Arrival tail coefficient in (1,2] (default 2) ALPHA_S - Service tail coefficient in (1,2] (default 2)

Returns:

W - Closed-form bound on the system time (Theorem 2) RHOHAT - Modified utilization (so that M/M/1 formulas still hold) SWORST - Exact worst-case system time over the uncertainty sets, eq. (12)

Reference: C. Bandi, D. Bertsimas, N. Youssef (2015). Robust Queueing Theory. Operations Research 63(3), 676-700.

qsys_mg1_ps(lambda, svc, svcparam, varargin)

RESULT = QSYS_MG1_PS(LAMBDA, ALPHA, T, …) RESULT = QSYS_MG1_PS(LAMBDA, BLST, M1, …)

Sojourn time distribution of the M/G/1 processor-sharing queue.

Jobs arrive in a Poisson stream of rate LAMBDA at a single egalitarian processor-sharing server whose service requirement has Laplace-Stieltjes transform bhat(tau) and mean M1. Writing V(x) for the sojourn time of a tagged job of service requirement x and RHO = LAMBDA*M1 < 1, Ott (1984) and Yashkov (1983) express the conditional transform as

E[exp(-s V(x))] = (1-RHO) / D(s,x),

where D(s,x) is the inverse Laplace transform, evaluated at x, of

f(tau;s) = [ (1-RHO)*tau^2 - (1-RHO)*LAMBDA*(1-bhat(tau))*tau

• s*RHO*tau - s*LAMBDA*(1-bhat(tau)) ]

/ [ tau^2 * (tau - s - LAMBDA*(1-bhat(tau))) ].

The transform is exact but implicit, since f must be inverted in tau. For phase-type service f(tau;s) is a proper rational function of tau, the double pole at tau = 0 cancels, and D(s,x) is obtained in closed form as a finite sum of residues (or, for repeated poles, from the matrix exponential of the companion realization). This makes the M/PH/1-PS queue, hence every service law LINE can fit with a phase-type distribution, exactly solvable. For a service transform supplied as a function handle, f is inverted in tau numerically on a Bromwich contour placed to the right of the dominant singularity tau*(s), the unique root of tau = s + LAMBDA*(1-bhat(tau)) in the right half plane, which the fixed-point iteration of that equation reaches at geometric rate RHO.

The conditional sojourn time is atomic on the lattice t = (k+1)*x, for k = 0,1,2,…: processor sharing gives every job in the system the same amount of work, so if the k jobs present on arrival all outlive the tagged job and no arrival intervenes, the sojourn is exactly (k+1)*x. For exponential service the masses are A_k = (1-RHO)*RHO^k*exp(-k*MU*x)*exp(-LAMBDA*(k+1)*x), the k = 0 term being the probability of finding the system empty and sharing it with nobody, which is the only one that stays exact for general service. The k = 0 atom is removed before inverting in s; the remaining atoms make cdfCond jump and leave no density, so pdfCond is NaN on the lattice.

Input:

lambda - Poisson arrival rate, positive scalar svc - phase-type initial probability vector alpha (1,n), or a

function handle bhat(tau) returning the service LST, which must accept complex arguments

svcparam - phase-type subgenerator T (n,n) when svc is a vector, or the

mean service time m1 when svc is a function handle

Optional name-value pairs:

‘x’ - service requirements to condition on (default []) ‘s’ - transform arguments at which to tabulate the LST (default []) ‘t’ - times at which to evaluate the sojourn time distribution

(default [])

‘nterms’ - function evaluations per numerical Laplace inversion, odd

(default 41)

‘pdf’ - service density handle, needed to remove the conditioning

when svc is a function handle (default [], filled in automatically on the phase-type path)

Output (struct):

rho - utilization LAMBDA*M1 m1, m2 - first two moments of the service requirement (m2 is NaN

when only a transform handle is supplied)

lstCond - handle (s,x) -> E[exp(-s V(x))] lstExcess - handle (s,x) -> E[exp(-s (V(x)-x))], bounded at large s lstUncond - handle s -> E[exp(-s V)], by quadrature over the density dominantRoot - handle s -> tau*(s) x, s - the requested grids lstCondVal - (numel(x),numel(s)) values of lstCond lstUncondVal - (1,numel(s)) values of lstUncond atomCond - (1,numel(x)) atom (1-RHO)*exp(-LAMBDA*x) at t = x atomUncond - (1-RHO)*bhat(LAMBDA), mass of the unshared jobs meanCond - (1,numel(x)) exact conditional mean x/(1-RHO) m2Cond - (1,numel(x)) conditional second moment, from the

transform derivatives

varCond - (1,numel(x)) conditional variance meanUncond - exact unconditional mean M1/(1-RHO) m2Uncond - unconditional second moment varUncond - unconditional variance t - the requested time grid pdfCond - (numel(x),numel(t)) density of V(x), NaN on the lattice cdfCond - (numel(x),numel(t)) P(V(x) <= t), atom included pdfUncond - (1,numel(t)) density of V cdfUncond - (1,numel(t)) P(V <= t)

References

T. J. Ott, “The sojourn-time distribution in the M/G/1 queue with processor sharing”, J. Appl. Prob. 21(2), 1984, pp. 360-378. S. F. Yashkov, “A derivation of response time distribution for an M/G/1 processor-sharing queue”, Probl. Contr. Inform. Theory 12, 1983, pp. 133-148. Q. Zhen, C. Knessl, “Asymptotic expansions for the sojourn time distribution in the M/G/1-PS queue”, Math. Meth. Oper. Res. 74, 2011, equations (2.2)-(2.5).

Copyright (c) 2012-2026, Imperial College London All rights reserved.

qsys_tandem_lindley(A, S, W0)

QSYS_TANDEM_LINDLEY Tandem network Lindley recursion on a sample path.

RESULT = QSYS_TANDEM_LINDLEY(A, S) propagates the waiting times of a series of K single-server FCFS stations in tandem, driven by the primitives of the sample path. A is the length-N vector of interarrival times at the first station, A(n) separating customers n and n+1, and S is the N-by-K matrix of service times, S(n,k) being the service time of customer n at station k. All waiting times start from zero.

RESULT = QSYS_TANDEM_LINDLEY(A, S, W0) starts customer 1 from the given length-K vector of waiting times instead of from an empty network.

At the first station this is Lindley’s recursion,

W(n+1,1) = max(W(n,1) + S(n,1) - A(n), 0).

Downstream the interarrival time is not a primitive: the arrival epoch of customer n at station k is its departure epoch from station k-1, so the interarrival time at station k is the interdeparture time upstream. Writing G(n,k) for the interarrival time at station k between customers n and n+1, with G(n,1) = A(n), the exact interdeparture identity is

G(n,k+1) = G(n,k) + W(n+1,k) - W(n,k) + S(n+1,k) - S(n,k),

equivalently and more transparently

G(n,k+1) = max(G(n,k) - W(n,k) - S(n,k), 0) + S(n+1,k),

an idle period at station k followed by the next customer’s service there. The recursion at station k is then

W(n+1,k) = max(W(n,k) + S(n,k) - G(n,k), 0).

Nothing here is distributional, so the recursion is exact for arbitrary interarrival and service times, dependent or not, and is the reference a simulated tandem sample path can be checked against directly. It reproduces a direct event-driven tandem simulation to 1e-12 over four stations.

Note that proposition 1 of the reference states this identity without the S(n+1,k) - S(n,k) term, which makes it wrong as a sample-path identity: the omitted difference has mean zero, so the mean interdeparture time survives, but individual waiting times do not. Implementing it as published gives station-1 waiting times that are correct and downstream ones that are not, by up to several mean service times. The form above is used instead.

Returns a struct with fields:

W - N-by-K waiting times, W(n,k) for customer n at station k G - N-by-K interarrival times, G(n,k) between customers n and n+1 at

station k, so G(:,1) is A; the last row is NaN, there being no customer N+1 to separate from

T - N-by-K sojourn times, W + S departure- N-by-K departure epochs of each customer from each station analyzer - Identifier string

Examples

A = exprnd(1/0.8, 1000, 1); S = exprnd(1, 1000, 2); r = qsys_tandem_lindley(A, S); mean(r.W(:,1)) % upstream mean wait mean(r.W(:,2)) % downstream mean wait

Reference: S. Palomo, J. Pender, “Learning the Tandem Network Lindley Recursion”, Proc. Winter Simulation Conference, 2021, equations 2 and 3 and proposition 1, the last corrected as described above; D. V. Lindley, “The Theory of Queues with a Single Server”, Proc. Camb. Phil. Soc. 48, 1952.

See also QSYS_MM1_TANDEM_LINDLEY, QSYS_MM1_LINDLEY, QSYS_HH1_LINDLEY

Copyright (c) 2012-2026, Imperial College London All rights reserved.

qsys_mm1_tandem_lindley(lambda, mu1, mu2, Wk, Wk1)

QSYS_MM1_TANDEM_LINDLEY Conditional downstream waiting time in an M/M/1 tandem.

RESULT = QSYS_MM1_TANDEM_LINDLEY(LAMBDA, MU1, MU2, WK, WK1) returns the conditional mean waiting time of customer n+1 at the downstream station of a two-station single-server tandem queue, given that customer n waited WK at the upstream station and WK1 at the downstream one. LAMBDA is the external arrival rate at the upstream station, MU1 and MU2 the two service rates. WK and WK1 are scalars or arrays of the same size.

The point of the tandem recursion is that the interarrival time at the downstream station is the interdeparture time upstream, not an independent draw. With A ~ Exp(LAMBDA) the interarrival time upstream and S1, S1’ the service times upstream of customers n and n+1, that interdeparture time is

D = max(A - WK - S1, 0) + S1’,

an idle period followed by the next service, and the downstream Lindley step is W2_{n+1} = (WK1 + S2 - D)^+ with S2 ~ Exp(MU2) independent of D.

Because A is exponential, max(A - WK - S1, 0) is zero with probability 1-q and Exp(LAMBDA) with probability

q = e^{-LAMBDA WK} MU1/(LAMBDA+MU1) = P(the upstream server goes idle),

so D is either Exp(MU1) or the sum of Exp(MU1) and Exp(LAMBDA). Averaging the downstream step over both cases needs only two elementary transforms of

g(d) = E[(WK1 + S2 - d)^+] = WK1 - d + 1/MU2 for d <= WK1,

= e^{-MU2 (d-WK1)}/MU2 for d > WK1,

namely J(c) = int_0^inf e^{-cu} g(u) du and Jw(c) = int_0^inf u e^{-cu} g(u) du, both closed form, giving

E[W2_{n+1} | WK, WK1] = (1-q) MU1 J(MU1) + q C, C = LAMBDA MU1 (J(MU1) - J(LAMBDA))/(LAMBDA-MU1) if LAMBDA ~= MU1,

= MU1^2 Jw(MU1) if LAMBDA == MU1.

As WK grows the upstream server never idles, q vanishes, and the mean tends to MU1 J(MU1) = E[g(S1’)], as it must.

Two caveats, both inherited from the reference and both quantified here.

First, this is exact for the step taken in isolation, that is when the conditioning pair is independent of the four primitives that drive the step. In a running tandem it is not: the downstream wait WK1 was itself determined by an interdeparture time containing S1, so conditioning on (WK, WK1) is not conditioning on a Markov state of the tandem. Measured against a 4e6-customer simulation of the real tandem at LAMBDA = 0.8, MU1 = MU2 = 1, the formula is within 0.4% to 1.3% away from the empty state and 7% at WK = WK1 = 0, where the entanglement is strongest. Treat it as exact for one isolated step and as a good approximation in a running tandem.

Second, this closed form was derived here rather than transcribed from the reference’s theorem 4, because that theorem rests on its proposition 1, which omits a service-time difference and so does not describe a tandem queue; see QSYS_TANDEM_LINDLEY. The two differ: at LAMBDA = 0.8, MU1 = MU2 = 1 and WK = WK1 = 0 the published route gives 0.3016 against 0.3457 here, the latter matching simulation of the step to 6e-4 relative error.

As in the reference, the upstream interarrival time is taken to be Exp(LAMBDA), which by Burke’s theorem is also the stationary interdeparture law, so the same formula is applied at any pair of consecutive stations of a longer M/M/1 tandem, with the caveat above compounding.

Returns a struct with fields:

mean - Conditional mean downstream waiting time, size of WK interdepMean - Conditional mean interdeparture time E[D | WK] = 1/MU1 + q/LAMBDA idleProb - q, the probability the upstream server idles, size of WK analyzer - Identifier string

Examples

r = qsys_mm1_tandem_lindley(0.8, 1.0, 1.0, 2.0, 3.0); r.mean % 2.9061 qsys_mm1_tandem_lindley(0.8, 1, 1, 0, 0).mean % 0.3457, both stations empty

Reference: S. Palomo, J. Pender, “Learning the Tandem Network Lindley Recursion”, Proc. Winter Simulation Conference, 2021, proposition 1 and theorem 4, corrected as described above.

See also QSYS_TANDEM_LINDLEY, QSYS_MM1_LINDLEY, QSYS_HH1_LINDLEY

Copyright (c) 2012-2026, Imperial College London All rights reserved.

qsys_lindley_moment(lambda, mu, Wn, m)

QSYS_LINDLEY_MOMENT One conditional Lindley moment for exponential primitives.

VALUE = QSYS_LINDLEY_MOMENT(LAMBDA, MU, WN, M) returns E[max(WN + S - A, 0)^M] with A ~ Exp(LAMBDA) and S ~ Exp(MU), for the vector of current waiting times WN and the integer order M >= 1.

This is the kernel shared by QSYS_MM1_LINDLEY, which calls it once per moment order, and QSYS_HH1_LINDLEY, which mixes it over the arrival and service phases. See QSYS_MM1_LINDLEY for the derivation and for why the upper incomplete gamma function reduces to a finite sum here.

See also QSYS_MM1_LINDLEY, QSYS_HH1_LINDLEY

Copyright (c) 2012-2026, Imperial College London All rights reserved.

qsys_hh1_lindley(lambda, pa, mu, ps, Wn, mmax)

QSYS_HH1_LINDLEY Conditional waiting-time moments of the Hl/Hn/1 Lindley recursion.

RESULT = QSYS_HH1_LINDLEY(LAMBDA, PA, MU, PS, WN) returns the exact conditional mean and variance of the waiting time of customer n+1 in an FCFS queue with hyperexponential interarrival and service times, given that customer n waited WN. LAMBDA and PA are the arrival phase rates and their probabilities, MU and PS the service phase rates and their probabilities. WN may be a vector.

RESULT = QSYS_HH1_LINDLEY(LAMBDA, PA, MU, PS, WN, MMAX) also returns the conditional raw moments of orders 1 to MMAX.

Hyperexponential primitives are mixtures of exponentials, so conditioning on the arrival phase i and the service phase j reduces one Lindley step to the M/M/1 step of QSYS_MM1_LINDLEY at rates LAMBDA(i) and MU(j), and the conditional moment is the corresponding mixture

E[W_{n+1}^m | W_n] = sum_i sum_j PA(i) PS(j) E_{ij}[W_{n+1}^m | W_n].

Phases are drawn independently for each customer, which is what makes the mixture exact rather than an approximation; a Markov-modulated arrival stream would not decompose this way.

Note that the variance is not the corresponding mixture of the per-phase variances, because the phase is itself random: it is recovered here from the first two mixed raw moments, which adds the between-phase spread of the means.

Returns a struct with fields:

mean - Conditional mean E[W_{n+1} | W_n], same size as WN var - Conditional variance Var[W_{n+1} | W_n], same size as WN moments - numel(WN) x MMAX conditional raw moments, orders 1 to MMAX mmax - Highest moment order computed analyzer - Identifier string

Examples

r = qsys_hh1_lindley([0.5 2], [0.4 0.6], [1 4], [0.7 0.3], 1.0); r.mean % 1.0484

% a degenerate mixture reproduces the M/M/1 result qsys_hh1_lindley(0.8, 1, 1.0, 1, 2.0).mean % 1.8902

Reference: S. Palomo, J. Pender, “Learning the Tandem Network Lindley Recursion”, Proc. Winter Simulation Conference, 2021, theorem 3. Verified against 4e6 Monte Carlo replications to 2e-3 relative error for m = 1, 2.

See also QSYS_MM1_LINDLEY, QSYS_TANDEM_LINDLEY, QSYS_MM1_TANDEM_LINDLEY

Copyright (c) 2012-2026, Imperial College London All rights reserved.

qsys_mm1_ps(lambda, mu)

[W,W2,ALPHA] = QSYS_MM1_PS(LAMBDA, MU)

Exact sojourn-time moments of the multiclass M/M/1-PS queue.

Class j arrives in a Poisson stream of rate LAMBDA(j) and requires an exponential amount of service with rate MU(j). The processor is shared equally by all jobs in service, so the class of a job affects its sojourn time both through its own service rate and through the mix of rates of the jobs it shares the processor with. With ALPHA = 1 - sum_j LAMBDA(j)/MU(j) the unutilized fraction of the processor, the moments of the sojourn time W_r of a tagged class-r job are

E[W_r] = 1/(ALPHA*MU(r)) E[W_r^2] = 2/(ALPHA*MU(r))^2 * [1 - sum_j lambda_j (mu_j-mu_r)/(mu_j(mu_j+mu_r))]

/ [1 - sum_j lambda_j/(mu_j+mu_r)]

which is equation (7) of Mitra and Morrison (1983). Both are exact, not asymptotic: the open system is the N -> infinity limit of the closed terminal-driven system whose moments that paper expands in 1/N, and the leading term of the expansion is exact in the limit. For a single class the second moment reduces to the classical 4/(mu^2 (1-rho)^2 (2-rho)) of Coffman, Muntz and Trotter (1970).

Input:

lambda - per-class Poisson arrival rates (1,R), non-negative mu - per-class exponential service rates (1,R), positive

Output:

W - per-class mean sojourn times (1,R) W2 - per-class second moments of the sojourn time (1,R) alpha - unutilized fraction of the processor, 1 - sum_j lambda_j/mu_j

Reference: D. Mitra, J. A. Morrison, “Asymptotic Expansions of Moments of the Waiting Time in Closed and Open Processor-Sharing Systems with Multiple Job Classes”, Adv. Appl. Prob. 15(4), 1983, equation (7).

Copyright (c) 2012-2026, Imperial College London All rights reserved.

qsys_mapg1k(D0, D1, svc, K, varargin)

RESULT = QSYS_MAPG1K(D0, D1, SVC, K)

Exact analysis of a MAP/G/1/K queue with tail drop: Markovian arrivals, arbitrary service time distribution F, and a finite buffer of K packets (the position held by the packet in transmission included).

Unlike QSYS_MAPG1, the service time is NOT fitted to a phase-type distribution: F enters exactly, through the functionals A_m and Q_m evaluated by uniformization of the arrival MAP. Unlike QSYS_MG1K_LOSS, which embeds the same way but assumes Poisson input, arrivals may be a general MAP, so flows with equal rate but different interarrival variability or autocorrelation are told apart.

D0, D1 - MAP parameter matrices (M x M), D0 + D1 an irreducible generator SVC - service time descriptor, a struct with field ‘type’:

‘gamma’fields alpha (shape), theta (scale). Covers Exp

(alpha=1) and Erlang (alpha integer).

‘det’ : field d (constant service time) ‘ph’ : fields alpha (1 x p row), T (p x p subgenerator) ‘density’ : field pdf (handle), optional field tmax

K - buffer size in packets, K >= 1

RESULT = QSYS_MAPG1K(…, ‘tol’, TOL) sets the uniformization truncation tolerance (default 1e-12). RESULT = QSYS_MAPG1K(…, ‘nmax’, N) caps the uniformization order.

Returns a struct with fields:

p0 - stationary probability of an empty buffer pK - stationary probability of a full buffer lossProbability - loss ratio of the aggregate arrival stream, 1-T/lambda throughput - aggregate throughput [pkts/s] lambda - aggregate arrival rate of the MAP meanServiceTime - S = E[service time] utilization - 1 - p0 rho - offered load lambda*S nmax - uniformization order used sigma - stationary law of the embedded chain, K*M entries pKvec - 1 x M, P(buffer full, phase j), summing to pK p0vec - 1 x M, P(buffer empty, phase j), summing to p0 plevel - 1 x (K+1), time-stationary P(level = l), l = 0..K meanQueueLength - E[number in system], sum_l l*plevel(l+1)

Method. The chain embedded at departure epochs is used, in the state (n,j): n = 0..K-1 packets left behind by a departure, j = MAP phase. With A_m the matrix of “m arrivals during a service, phase i -> j”,

n >= 1: n’ = n-1+min(m, K-n), via A_m, overflow sum_{m>=K-n} A_m n == 0: the phase first jumps by (-D0)^{-1}*D1 (the idle period ends at

an arrival), and the service then proceeds as from n = 1.

Its stationary law sigma gives, by Markov renewal reward,

E[cycle] = S + sum_j sigma(0,j)*idle_j, idle = (-D0)^{-1}*e T = 1/E[cycle], p0 = (sum_j sigma(0,j)*idle_j)/E[cycle] = 1 - T*S pK = E[time at level K per cycle]/E[cycle], from Q_m,

where Q_m is the expected time within a service with exactly m arrivals so far. Time-stationary p0 and pK follow, so no PASTA assumption is needed on the MAP side.

This is not the transform solution of Theorem 1 of [1], which is stated in terms of a sequence R_m obeying R(z) = z*(A(z)-z*I)^{-1}. That sequence grows geometrically, at a rate set by the smallest zero of det(A(z)-z*I), while the quantity extracted from it stays O(s) as s -> 0+; cond(G(s)) therefore grows like that ratio^K and crosses the double-precision ceiling near K = 20 for the flows of [1] under gamma service with CV = 2, and near K = 10 under constant service, where A_0 = exp(D0*d) has entries O(1e-8). Reference [1] evaluates its formulae in arbitrary precision, so the restriction is invisible there. The embedded chain used here has every entry a probability or a time and is stable for any K and any F.

TEST (M/M/1/5, exact loss 0.04812030):

r=qsys_mapg1k(-2,2,struct(‘type’,’gamma’,’alpha’,1,’theta’,1/3),5); r.pK

References: [1] Chydzinski, A. Per-Flow Throughput of a FIFO Buffer. Applied System

Innovation 2026, 9, 112.

[2] Niu, Z.; Cooper, R.B. Transform-Free Analysis of M/G/1/K and Related

Queues. Mathematics of Operations Research 1993, 18, 486-510.

See also QSYS_MAPG1K_PERFLOW, QSYS_MAPG1, QSYS_MG1K_LOSS, QSYS_MMCK.

qsys_ldps_workload(lambda, B, alpha, N, t, ngrid)

[F,t,p] = QSYS_LDPS_WORKLOAD(lambda,B,alpha,N,t) - Stationary distribution of the quantity of work in a single-stage load-dependent processor sharing station with Poisson arrivals and blocking.

Assumed model (Cohen, 1979, Sect. 9; the model of Sect. 7 with one stage):

• a single service stage fed by a Poisson arrival stream of rate lambda;

• blocking capacity N: a request arriving when N requests are already present is lost and leaves no trace on the state;

• generalized processor sharing: when x requests are present each of them accrues service at rate f(x), so the stage completes work at total rate x*f(x). This function is parametrized by the LINE load-dependent total rate scaling alpha(x)=x*f(x), i.e. the argument of setLoadDependence at a PS station;

• required service times are i.i.d. with absolutely continuous distribution B and finite mean beta.

The result is insensitive to B beyond its shape only through the equilibrium residual distribution Psi below; the state probabilities p depend on B only through beta.

Let psi_t denote the total amount of service still to be given to the requests present at time t. Cohen (1979) eqs. (9.1)-(9.3) give

Pr{psi_t < psi} = sum_{h=0}^{N} p_h Psi^{h*}(psi), p_h = (rho^h/h!) phi(h) / sum_{k=0}^{N} (rho^k/k!) phi(k), rho=lambda*beta phi(h) = 1/prod_{k=1}^{h} f(k), phi(0)=1, Psi(psi) = int_0^psi (1-B(v))/beta dv,

with Psi^{h*} the h-fold convolution of Psi and Psi^{0*} degenerate at zero. Substituting f(k)=alpha(k)/k the factorial cancels, leaving

p_h propto rho^h / prod_{k=1}^{h} alpha(k),

which is the familiar load-dependent birth-death form. Psi is the equilibrium (residual life) distribution of B, so psi_t is a mixture of h-fold convolutions of residual service times, with an atom p_0 at zero.

Inputs:

lambda : Poisson arrival rate (finite, positive) B : Distribution object for the required service time (continuous,

finite positive mean)

alpha : rate scaling alpha(n)=n*f(n) for n=1..N (finite, positive) N : blocking capacity (finite positive integer) t : optional grid at which the CDF is returned. Default: an

automatically sized grid covering the bulk of the distribution.

ngridoptional number of points of the internal uniform quadrature

grid on which the convolutions are formed. Default 2001.

Outputs:

FF(j) = Pr{psi_t <= t(j)}. Note F(1)=p(1)=Pr{psi_t=0} when t(1)=0,

since the workload has an atom at zero of size p_0.

t : grid at which F is reported p : (1,N+1) vector, p(h+1)=Pr{x_t=h}, the stationary number in system

This is the model of Cohen (1979) Sect. 9 only. It is not the weighted GPS/DPS discipline of SchedStrategy.GPS, whose per-class weights this formula does not represent.

Reference: J.W. Cohen, “The multiple phase service network with generalized processor sharing”, Acta Informatica 12, 245-284 (1979), Sect. 9.

qsys_bmapm1(D, mu, varargin)

QSYS_BMAPM1 Analyzes a BMAP/M/1 queue by the matrix-analytic (M/G/1-type) method.

RESULT = QSYS_BMAPM1(D, MU) analyzes a single-server queue fed by a batch Markovian arrival process and with exponential service of rate MU.

Inputs:

D - cell array {D0, D1, …, DK} of BMAP matrices. D0 carries the hidden

transitions, Dk (k >= 1) the transitions that release a batch of k customers.

MU - exponential service rate.

Optional parameters:

‘Uniformization’ - uniformization constant q used to randomize the

generator into a discrete-time M/G/1-type chain. It must dominate every total outflow rate; by default it is chosen as max_i(-D0(i,i)) + mu, rounded up.

‘MaxIter’ - maximum functional iterations for G (default 10000) ‘Tolerance’ - convergence tolerance for G (default 1e-12) ‘MaxLevel’ - level truncation used for the queue-length distribution

(default: adaptive, see qsys_bmapphnn_retrial)

‘TailTolerance’ - relative truncation target for the level distribution

(default 1e-10)

Beyond the usual performance measures the result exposes the intermediate matrix-analytic quantities themselves, so that the algorithm can be inspected and taught rather than only its output:

theta - stationary vector of the BMAP phase process, sum_k D_k lambda - mean arrival rate, theta * sum_k k*D_k * e rho - offered load lambda/mu q - uniformization constant actually used A0, A1, Bk - randomized blocks: A0 = (mu/q)I is a service completion

(level down by one), A1 = (1/q)(D0 - mu*I) + I keeps the level, Bk{k} = (1/q)D_k raises it by k

B0 - boundary local block (1/q)D0 + I, used at level 0 where no

service can complete

A - A0 + A1 + sum_k Bk{k}, the phase process of the chain alpha - stationary vector of A G - minimal non-negative solution of

G = A0 + A1*G + sum_k Bk{k}*G^(k+1)

drift - alpha*(sum_k k*Bk{k})*e - alpha*A0*e. The queue is stable

iff this is strictly negative

decayRate - geometric decay rate of the level probabilities, measured as

the limiting ratio pi_(n+1)/pi_n. Reported rather than derived from a spectral convention so that it is unambiguous

levelProb - level probabilities pi_n as rows (level 0 first) pi0 - probability the system is empty (equals 1-rho exactly)

Example

% Example 6.4 of Bolch et al. D0 = [-2, 1/2; 1/3, -3]; D1 = [1/4, 1/2; 1/3, 1]; D2 = [1/4, 1/2; 1, 1/3]; result = qsys_bmapm1({D0, D1, D2}, 11);

See also qsys_mapm1(), qsys_mapph1(), qsys_bmapphnn_retrial()

Copyright (c) 2012-2026, Imperial College London All rights reserved.

qsys_mmapg1k(D0, D1c, svc, K, varargin)

RESULT = QSYS_MMAPG1K(D0, D1C, SVC, K)

Exact per-class throughput and loss ratio of an MMAP[K]/G/1/K queue with tail drop: marked Markovian arrivals, arbitrary service time distribution F common to all classes, and a finite buffer of K packets (the position held by the packet in transmission included).

Two classes of equal arrival rate but different interarrival variability or autocorrelation receive different loss ratios, which is the effect that motivates [1]. Aggregate-only finite-buffer analyses cannot express it: they return a single blocking probability p and set T_k = lambda_k*(1-p), making the loss ratio identical across classes by construction.

D0 - M x M hidden transition matrix of the arrival MMAP D1C - 1 x R cell array, D1C{k} = M x M arrival matrix of class k.

D0 + sum_k D1C{k} must be an irreducible generator.

SVC - service time descriptor, see QSYS_MAPG1K K - buffer size in packets, K >= 1

Options are forwarded to QSYS_MAPG1K (‘tol’, ‘nmax’).

Returns a struct with fields:

throughput(k) - throughput of class k [pkts/s] lossRatio(k) - loss ratio of class k, in [0,1] lambda(k) - arrival rate of class k [pkts/s] lambdaAggregate, throughputAggregate, lossAggregate p0, pK - empty/full buffer probabilities pKvec - 1 x M, P(buffer full, phase j) plevel - 1 x (K+1), time-stationary P(level = l) meanQueueLength - E[number in system] meanServiceTime, rho, utilization

Method. The aggregate MAP {D0, sum_k D1C{k}} drives QSYS_MAPG1K, whose embedded chain returns the joint law of buffer level and MAP phase. A class-k arrival leaves phase i at rate (D1C{k}*e)_i, so the rate of class-k arrivals that meet a full buffer is pKvec*D1C{k}*e, and

lambda_k = pi*D1C{k}*e, L_k = (pKvec*D1C{k}*e)/lambda_k,

with pi the stationary phase law of the aggregate. This is exact: no independence between classes is assumed and no PASTA argument is used, the phase resolution of pKvec doing the work instead.

Relation to [1]. Reference [1] instead keeps one flow exact and replaces the rest by a Poisson stream of the same rate, invoking Palm-Khinchin, and repeats that once per flow. That approximation is needed only when the joint arrival process is unavailable, its exact form costing prod_n M_n phases for N independent flows. When the joint MMAP is already at hand, as it is inside a solver that propagates MMAPs between stations, the phase resolution above is both exact and cheaper. Use QSYS_MAPG1K_PERFLOW for the setting of [1], where flows are given as N separate MAPs.

Assumes a single server and a service law that is iid and independent of class. Per-class service makes the departure rate depend on which class holds the server, which this model does not represent.

References: [1] Chydzinski, A. Per-Flow Throughput of a FIFO Buffer. Applied System

Innovation 2026, 9, 112.

See also QSYS_MAPG1K, QSYS_MAPG1K_PERFLOW, QSYS_MG1K_LOSS.

qsys_mapg1k_perflow(MAPS, svc, K, varargin)

RESULT = QSYS_MAPG1K_PERFLOW(MAPS, SVC, K)

Per-flow throughput and loss ratio of a FIFO buffer with tail drop that is fed by N flows of arbitrary, mutually different statistical character.

Flow n is described by its own MAP, so two flows may share an arrival rate and still differ in the shape and autocorrelation of their interarrival times. The buffer holds K packets including the one in transmission, and the transmission time follows an arbitrary distribution F.

MAPS - 1 x N cell array, MAPS{n} = {D0n, D1n}, the MAP of flow n. The

modulating orders M_n may differ from flow to flow.

SVC - service time descriptor, see QSYS_MAPG1K K - buffer size in packets, K >= 1

Options are forwarded to QSYS_MAPG1K (‘tol’, ‘nmax’).

Returns a struct with fields:

throughput(n) - throughput of flow n [pkts/s] lossRatio(n) - loss ratio of flow n, in [0,1] lambda(n) - arrival rate of flow n [pkts/s] lambdaAggregate - sum_n lambda(n) throughputAggregate- sum_n throughput(n) lossAggregate - aggregate loss ratio, sum_n L(n)*lambda(n)/lambda p0(n), pK(n) - empty/full buffer probabilities of the n-th model rho - offered load lambdaAggregate * E[S]

Method. The exact model of N flows would need a Markov chain tracking the modulating state of every flow jointly with the buffer occupancy, hence prod_n M_n * (K+1) states; [1] notes this is already out of reach at N=10, M_n=3, K=10 (over 4e11 transition matrix entries). Instead, one model per flow is solved: flow n is kept exactly as MAP_n, while the other N-1 flows are replaced by a single Poisson stream of rate lambda - lambda_n. That substitution is justified by the Palm-Khinchin limiting theorem on the superposition of many point processes, so it is an approximation that improves as N grows, and it is applied N times, once per flow, so no flow is ever the one being Poissonized when its own throughput is computed. Superposing MAP_n with the Poisson background yields the MAP

D0 = D0n - lambdaBar_n*I, D1 = D1n + lambdaBar_n*I, ([1], eq. 5)

which is passed to QSYS_MAPG1K. Each solve costs O((K*M)^3) for the embedded chain, so the whole sweep is O(N*(K*M)^3) against the O(M^(3N)*K^3) of the exact joint model: linear rather than exponential in the flow count.

Accuracy. [1] reports errors against simulation of the exact model below about 8% for N >= 9 with K >= 20, falling to 2.1% at K=50 and 0.5% at K=100, and to 1.2% at N=900. Errors are largest when flows are few, highly variable, and the buffer is small.

TEST (Table 2 of [1], K=20, rho=1, gamma service with CV=2; the nine flows are defined in eq. (55)-(63)). Tables 2-9 of [1] all reproduce to within 0.07 pkts/s, the precision to which they are published:

throughput ~ [89.8 179.6 269.4 85.2 169.3 253.4 56.9 110.8 166.6] pkts/s

References: [1] Chydzinski, A. Per-Flow Throughput of a FIFO Buffer. Applied System

Innovation 2026, 9, 112. Theorem 1.

See also QSYS_MAPG1K, QSYS_MG1K_LOSS, QSYS_MAPG1.

qsys_mm1_dps(lambda, mu, w, tol, maxCutoff)

[T,RHO] = QSYS_MM1_DPS(LAMBDA, MU, W, TOL, MAXCUTOFF)

Numerically exact M/M/1 Discriminatory Processor Sharing (DPS) queue.

Solves the multiclass DPS continuous-time Markov chain on the per-class population vector (n_1..n_K): arrivals lambda(k), class-k completion rate mu(k)*n_k*w(k)/sum_j n_j*w(j). The state space is truncated at a total population level chosen from the geometric tail bound and doubled until the per-class mean counts are stable, so the result is exact to solver precision and conserves the M/M/1 total for equal service rates by construction.

Input:

lambda - per-class Poisson arrival rates (1,K) mu - per-class exponential service rates (1,K) w - per-class DPS weights (1,K), positive tol - convergence tolerance on the mean counts (default 1e-10) maxCutoff - hard bound on the truncation level (default 2048)

Output:

T - per-class mean response times (1,K) via Little’s law rho - total utilization sum_k lambda(k)/mu(k)

Copyright (c) 2012-2026, Imperial College London All rights reserved.

qsys_phph1(alpha, T, beta, S, varargin)

QSYS_PHPH1 Analyzes a PH/PH/1 queue using BUTools MMAPPH1FCFS.

RESULT = QSYS_PHPH1(ALPHA, T, BETA, S) analyzes a PH/PH/1 queue with:

ALPHA - Arrival PH initial probability vector (1 x n) T - Arrival PH generator matrix (n x n) BETA - Service PH initial probability vector (1 x m) S - Service PH generator matrix (m x m)

The arrival PH is converted to an equivalent MAP representation.

RESULT = QSYS_PHPH1(…, ‘numQLMoms’, K) computes K queue length moments RESULT = QSYS_PHPH1(…, ‘numQLProbs’, N) computes N queue length probs RESULT = QSYS_PHPH1(…, ‘numSTMoms’, K) computes K sojourn time moments

Returns a struct with fields:

meanQueueLength - Mean number of customers in system meanWaitingTime - Mean waiting time in queue meanSojournTime - Mean sojourn time (waiting + service) utilization - Server utilization queueLengthDist - Queue length distribution P(Q=n) queueLengthMoments- Raw moments of queue length sojournTimeMoments- Raw moments of sojourn time analyzer - Name of analyzer used

See also MMAPPH1FCFS, qsys_mapph1(), qsys_mapmap1()

qsys_phmc(alpha, T, mu, c, varargin)

QSYS_PHMC Exact PH/M/c (Neuts’ matrix-geometric).

Solves the GI/M/c QBD via R^2*A2 + R*A1 + A0 = 0, where

A0 = (-T*1)*alpha, A1 = T - c*mu*I, A2 = c*mu*I.

pi_n = pi_c * R^(n-c) for n >= c. Boundary states pi_0..pi_c are obtained from balance equations + normalization sum_{n<c} pi_n*1 + pi_c*(I-R)1 = 1.

Inputs:

ALPHA - PH entry probability row vector (1 x k) T - PH sub-generator (k x k) MU - Exponential service rate per server (positive scalar) C - Number of servers (positive integer)

Optional name-value:

‘maxIter’ - max iterations for R fixed-point (default 50000) ‘tol’ - convergence tolerance (default 1e-14)

RESULT struct fields:

meanQueueLength, meanWaitingQueue, meanWaitingTime, meanSojournTime, utilization, analyzer

qsys_phm1(alpha, T, mu)

QSYS_PHM1 Exact PH/M/1 (GI/M/1 with phase-type inter-arrivals).

RESULT = QSYS_PHM1(ALPHA, T, MU) solves the GI/M/1 sigma-root

sigma = psi_A(mu*(1 - sigma)), psi_A(s) = ALPHA*(s*I - T)(-T*1)

and returns time-average performance metrics.

Inputs:

ALPHA - PH entry probability row vector (1 x k) T - PH sub-generator (k x k) MU - Exponential service rate (positive scalar)

RESULT struct fields:

meanQueueLength - L = rho/(1 - sigma) meanWaitingQueue - Lq = rho*sigma/(1 - sigma) meanWaitingTime - Wq = Lq/lambda meanSojournTime - W = Wq + 1/MU utilization - rho = lambda/MU sigma - GI/M/1 root in (0, 1) analyzer - identifier string

See also qsys_gm1(), qsys_dmc()

qsys_mxm1(lambda_batch, mu, E_X, E_X2_or_Var_X, varargin)

QSYS_MXM1 - MX/M/1 queue with batch arrivals

[W, Wq, U, Q] = QSYS_MXM1(LAMBDA_BATCH, MU, E_X, E_X2) [W, Wq, U, Q] = QSYS_MXM1(LAMBDA_BATCH, MU, BATCH_SIZES, PMF)

Analytical solution for MX/M/1 queueing system with batch Markovian arrivals and exponential service.

Input formats:

1. Moment-based: qsys_mxm1(lambda_batch, mu, E_X, E_X2) - lambda_batch: Batch arrival rate - mu: Service rate - E_X: Mean batch size - E_X2: Second moment of batch size

2. PMF-based: qsys_mxm1(lambda_batch, mu, batch_sizes, pmf) - lambda_batch: Batch arrival rate - mu: Service rate - batch_sizes: Array of batch sizes (e.g., [1, 2, 4, 8]) - pmf: Probability mass function for batch sizes

3. Variance-based: qsys_mxm1(lambda_batch, mu, E_X, Var_X, ‘variance’) - lambda_batch: Batch arrival rate - mu: Service rate - E_X: Mean batch size - Var_X: Variance of batch size - ‘variance’: Flag to indicate variance mode

Outputs:

W - Mean time in system Wq - Mean waiting time in queue U - Server utilization Q - Mean queue length (including service)

The formula accounts for both queueing delay and internal batch delay:

Wq = rho/(mu*(1-rho)) + (E[X²] - E[X])/(2*mu*E[X]*(1-rho))

^M/M/1 term ^Internal batch delay

Copyright (c) 2012-2026, Imperial College London All rights reserved.

qsys_mmk(lambda, mu, k)

W=QSYS_MMK(LAMBDA,MU,K)

qsys_mmck(lambda, mu, c, K)

QSYS_MMCK Exact closed-form analysis of an M/M/c/K queue.

RESULT = QSYS_MMCK(LAMBDA, MU, C, K) analyzes an M/M/c/K queue with:

LAMBDA - Poisson arrival rate (positive scalar) MU - Exponential service rate per server (positive scalar) C - Number of servers (positive integer) K - System capacity, total jobs allowed (integer, K >= C)

Stationary distribution (Erlang-B/C truncated form):

a = lambda/mu, rho = a/c p_n = (a^n / n!) * p_0 for 0 <= n <= c p_n = (a^c / c!) * rho^(n-c) * p_0 for c <= n <= K p_0 = 1 / [ sum_{n=0}^{c-1} a^n/n! + (a^c/c!) * sum_{n=0}^{K-c} rho^n ]

Returns a struct with fields:

meanQueueLength - Mean number of jobs in system, L meanQueueLengthQ - Mean number waiting in queue, Lq meanWaitingTime - Mean waiting time in queue, Wq (Little) meanSojournTime - Mean sojourn time, W = Wq + 1/mu utilization - Per-server utilization, lambda_eff/(c*mu) throughput - Effective throughput, lambda*(1 - p_K) lossProbability - Blocking probability, p_K queueLengthDist - Distribution as 1x(K+1) row: [p_0, p_1, …, p_K] analyzer - Identifier string

Examples

r = qsys_mmck(1, 1, 2, 4); % M/M/2/4 with rho = 0.5 r.meanQueueLength % 1.1304 r.lossProbability % 0.0435

See also QSYS_MM1K_LOSS, QSYS_MMK, QSYS_MAPDC

Copyright (c) 2012-2026, Imperial College London All rights reserved.

qsys_mmcc_retrial_fp(lambda, mu, c, tol, maxiter)

[blocProb, r, niter] = QSYS_MMCC_RETRIAL_FP(LAMBDA, MU, C, TOL, MAXITER)

Fixed-point approximation for M/M/c/c retrial queues.

Customers arrive at rate LAMBDA to a system with C servers, each with service rate MU. Blocked customers join an orbit and retry. Under the assumption that the retrial rate is small relative to the service rate, the total arrival flow (fresh + retrial) is approximated by a Poisson process with rate LAMBDA + r, where r satisfies the fixed-point equation:

r = (lambda + r) * B(lambda/mu + r/mu, c)

and B(a, c) is the Erlang-B blocking probability for offered load a and c servers.

INPUT:

lambda : arrival rate mu : service rate per server c : number of servers (= capacity, no waiting room) tol : convergence tolerance (default: 1e-10) maxiter : maximum iterations (default: 10000)

OUTPUT:

blocProb : blocking probability (fraction of arrivals lost or retried) r : additional arrival rate due to retrials niter : number of iterations to converge

REFERENCE:

Cohen (1957), fixed-point approximation for M/M/c/c retrial queues. Phung-Duc, “Retrial Queueing Models: A Survey on Theory and Applications”, 2019, Eq. (1).

qsys_mm1k_loss(lambda, mu, K)

Niu-Cooper, Transform-Free Analysis of M/G/1/K and Related Queues, Mathematics of Operations Research Vol. 18, No. 2 (May, 1993), pp. 486-510 (25 pages) TEST: mu=3; lambda=2; K=5; [sigma,rho,lossprob]=qsys_mg1k_loss(lambda,@(t)mu.*exp(-mu.*t),K)

qsys_mm1(lambda, mu)

W=QSYS_MM1(LAMBDA,MU)

qsys_mginf(lambda, mu, varargin)

[L,Lq,W,Wq,p0]=QSYS_MGINF(LAMBDA,MU) [L,Lq,W,Wq,p0,pk]=QSYS_MGINF(LAMBDA,MU,K)

Exact solution for M/G/∞ queue (infinite servers). Performance is independent of the service time distribution shape (G). The number of customers in the system follows a Poisson distribution.

Input:

lambda - arrival rate mu - service rate (mean service time = 1/mu) K (optional) - state for probability computation

Output:

L - average number of customers in system Lq - average number of customers in queue (always 0) W - average time in system (= 1/mu) Wq - average waiting time in queue (always 0) p0 - probability of empty system pk (optional) - probability of exactly k customers in system

qsys_mg1k_loss_mgs(lambda, mu, mu_scv, K)

10. MacGregor Smith - Optimal Design and Performance Modelling of M/G/1/K Queueing Systems

qsys_mapdc(D0, D1, s, c, varargin)

QSYS_MAPDC Analyzes a MAP/D/c queue using Q-MAM.

RESULT = QSYS_MAPDC(D0, D1, S, C) analyzes a MAP/D/c queue with:

D0 - MAP hidden transition matrix (n x n) D1 - MAP arrival transition matrix (n x n) S - Deterministic service time (positive scalar) C - Number of servers (positive integer)

RESULT = QSYS_MAPDC(…, ‘maxNumComp’, N) sets max queue length components (default 1000) RESULT = QSYS_MAPDC(…, ‘numSteps’, K) sets waiting time distribution granularity (default 1)

Returns a struct with fields:

meanQueueLength - Mean number of customers in system meanWaitingTime - Mean waiting time in queue meanSojournTime - Mean sojourn time (waiting + service) utilization - Server utilization (per server) queueLengthDist - Queue length distribution P(Q=n) waitingTimeDist - Waiting time CDF at discrete points analyzer - Name of analyzer used

The waiting time distribution is evaluated at points {0, s/numSteps, 2*s/numSteps, …} where s is the deterministic service time and numSteps controls the granularity.

See also Q_CT_MAP_D_C, qsys_mapmc(), qsys_mapd1()

qsys_mg1k_loss(lambda, svc_density, K)

[LOSSPROB,RHO]=QSYS_MG1K_LOSS(LAMBDA,SVC_DENSITY,K)

Exact M/G/1/K loss probability via the Markov chain embedded at service-start epochs (transform-free analysis in the spirit of Niu-Cooper).

State: number of customers waiting in the queue immediately after a service start, q in {0,…,K-2} (capacity K includes the job in service; just after a departure at most K-1 jobs remain, one of which enters service). With a_j = P(j Poisson arrivals during a service time):

q=0if no arrival occurs during the service the system empties and

the next service starts with the next arrival (q’=0), so both a_0 and a_1 lead to q’=0 and j>=2 arrivals lead to q’=j-1;

q>=1: q’ = q-1+j, with arrivals beyond the free capacity lost

(aggregated in the last column).

The loss probability follows from the renewal-reward argument

E[cycle] = E[S] + sigma_0*a_0/lambda, lambda_eff = 1/E[cycle], P_loss = 1 - lambda_eff/lambda = 1 - 1/(rho + sigma_0*a_0)

where sigma is the stationary distribution at service-start epochs.

TEST: mu=3; lambda=2; K=5; [lossprob_mg1k,rho]=qsys_mg1k_loss(lambda,@(t)mu.*exp(-mu.*t),K) (exact M/M/1/5 value: 0.04812030)

Reference: Niu, Cooper. Transform-Free Analysis of M/G/1/K and Related Queues. Mathematics of Operations Research 18(2), 1993, 486-510.

qsys_gg1(lambda, mu, ca2, cs2)

[W,RHOHAT]=QSYS_GG1(LAMBDA,MU,CA2,CS2) analyzes a G/G/1 queue.

Uses exact methods for special cases (M/M/1, M/G/1, G/M/1) and Allen-Cunneen approximation for the general case. In the G/M/1 case, the interarrival-time distribution is fitted from (LAMBDA,CA2) by a two-moment renewal process (H2 with balanced means for CA2>1, mixed Erlang for CA2<1) and sigma is the root of sigma = A*(mu*(1-sigma)), with A* the interarrival-time LST.

Inputs:

LAMBDA - Arrival rate MU - Service rate CA2 - Squared coefficient of variation of inter-arrival time CS2 - Squared coefficient of variation of service time

Returns:

W - Average time in system (response time) RHOHAT - Modified utilization (so that M/M/1 formulas still hold)

qsys_mg1_srpt(lambda, mu, cs)

QSYS_MG1_SRPT Compute mean response time for M/G/1/SRPT queue

[W, RHO] = QSYS_MG1_SRPT(LAMBDA, MU, CS) computes the mean response time for each job class in an M/G/1 queue with Shortest Remaining Processing Time (SRPT) scheduling.

SRPT is a size-based policy: it always serves the job with the smallest remaining processing time, preempting whenever a shorter job arrives. The class-conditional response time is obtained from the Schrage-Miller formula (Eqs (1)-(3) of Bansal-Harchol-Balter, SIGMETRICS 2003, citing Schrage-Miller 1966). For a job of size x:

E[T(x)] = E[W(x)] + E[R(x)] E[W(x)] = lambda*(m2(x) + x^2*(1-F(x))) / (2*(1-rho(x))^2) (waiting) E[R(x)] = integral_0^x dt/(1-rho(t)) (residence)

where f(t) is the overall (mixture) job-size density, F(t) its CDF,

rho(x) = lambda * integral_0^x t*f(t) dt (load from jobs of size <= x) m2(x) = integral_0^x t^2 f(t) dt

The per-class mean is E[T_r] = integral_0^inf E[T(x)] f_r(x) dx, with f_r the class-r size density and f = sum_r (lambda_r/lambda) f_r the mixture. Because E[T(x)] depends only on the job size (SRPT is size-based, not class-based), this integral is exact. The integrals are evaluated by cumulative trapezoidal quadrature on a common grid.

Each class is represented by a job-size distribution matched to its (mean=1/mu_r, scv=cs_r^2): exponential when cs_r=1, a two-phase balanced hyperexponential when cs_r>1, and a Tijms mixture of Erlang-(k-1)/Erlang-k when cs_r<1. For the fully exponential case this reproduces the exact M/M/1/SRPT hyperexponential-mixture result.

Parameters:

• lambda – Vector of arrival rates per class

• mu – Vector of service rates per class

• cs – Vector of coefficients of variation per class (cs=1 for exponential)

Returns:

W – Vector of mean response times per class (original class order) rho : System load measure Q/(1+Q) with Q = sum(lambda.*W)

References

• L. E. Schrage and L. W. Miller, “The queue M/G/1 with the shortest remaining processing time discipline”, Operations Research, 14:670-684, 1966.

• N. Bansal and M. Harchol-Balter, “Analysis of SRPT scheduling: investigating unfairness”, SIGMETRICS 2001, Sec. 4, Eqs (1)-(3).

Copyright (c) 2012-2026, Imperial College London All rights reserved.

qsys_mg1_setf(lambda, mu, cs)

QSYS_MG1_SETF Compute mean response time for M/G/1/SETF queue

[W, RHO] = QSYS_MG1_SETF(LAMBDA, MU, CS) computes the mean response time for each job class in an M/G/1 queue with SETF (Shortest Elapsed Time First) scheduling.

SETF is the non-preemptive version of FB/LAS (Feedback/Least Attained Service). Under SETF: - Jobs are ordered by their attained service (elapsed processing time) - The job with the least attained service has highest priority - However, once a job begins service, it runs to completion (non-preemptive)

The difference from FB/LAS is that arriving jobs with less attained service must wait until the currently serving job completes, rather than preempting it.

For SETF, the mean response time follows a modified FB formula that accounts for the non-preemptive nature. The formula includes an additional residual service time term due to the non-preemptive blocking:

E[T(x)]^SETF = E[T(x)]^FB + E[R] / (1 - rho_x)

where E[R] is the mean residual service time and rho_x is the truncated load.

Parameters:

• lambda – Vector of arrival rates per class

• mu – Vector of service rates per class

• cs – Vector of coefficients of variation per class

Returns:

W – Vector of mean response times per class rho : Overall system utilization (modified for Little’s law)

References

• M. Nuyens and A. Wierman, “The Foreground-Background queue: A survey”, Performance Evaluation, 2008.

• A. Wierman and M. Harchol-Balter, “Classifying scheduling policies with respect to unfairness in an M/GI/1”, SIGMETRICS 2003.

Copyright (c) 2012-2026, Imperial College London All rights reserved.

qsys_mg1_psjf(lambda, mu, cs)

QSYS_MG1_PSJF Compute mean response time for M/G/1/PSJF queue

[W, RHO] = QSYS_MG1_PSJF(LAMBDA, MU, CS) computes the mean response time for each job class in an M/G/1 queue with Preemptive Shortest Job First (PSJF) scheduling.

Under PSJF, priority is based on a job’s original size (not remaining size). Jobs with smaller original sizes always preempt jobs with larger sizes.

For PSJF, the mean response time for a job of size x is given by (Section 3.2 of Wierman-Harchol-Balter 2003):

E[T(x)]^PSJF = (lambda * integral_0^x t^2*f(t)dt) / (2*(1-rho(x))^2)

• x / (1 - rho(x))

where:

• rho(x) = lambda * integral_0^x t*f(t)dt (truncated load)

• f(t) is the service time density (mixture of exponentials)

For exponential service, the expected response time for class k is computed by integrating E[T(x)] over the exponential distribution of class k sizes.

CLASSIFICATION (Wierman-Harchol-Balter 2003):

PSJF is “Always Unfair” - some job size is treated unfairly under all loads and all service distributions.

Parameters:

• lambda – Vector of arrival rates per class

• mu – Vector of service rates per class

• cs – Vector of coefficients of variation per class (cs=1 for exponential)

Returns:

W – Vector of mean response times per class rho : Overall system utilization (modified for Little’s law)

References

• A. Wierman and M. Harchol-Balter, “Classifying scheduling policies with respect to unfairness in an M/GI/1”, SIGMETRICS 2003, Section 3.2.

• L. Kleinrock, “Queueing Systems, Volume II: Computer Applications”, Wiley, 1976.

Copyright (c) 2012-2026, Imperial College London All rights reserved.

qsys_mg1_prio(lambda, mu, cs)

W=QSYS_MG1_PRIO(LAMBDA,MU,CS)

Analyzes an M/G/1 queueing system with non-preemptive (Head-of-Line) priorities.

This function computes per-class mean waiting times for multiple priority classes using the Pollaczek-Khinchine formula extended for priorities.

For K priority classes (class 1 = highest priority), the waiting time for class k is:

W_k = (B_0 / (1 - sum_{i=1}^{k-1} rho_i)) * (1 / (1 - sum_{i=1}^{k} rho_i))

Where:

• rho_i = lambda_i / mu_i (utilization of class i)

• B_0 = sum_{i=1}^{K} lambda_i * (1/mu_i^2) * (1 + cs_i^2)

• cs_i = coefficient of variation for class i service time

Parameters:

• lambda – Vector of arrival rates per priority class (class 1 = highest)

• mu – Vector of service rates per priority class

• cs – Vector of coefficients of variation per priority class

Returns:

W – Vector of mean waiting times per priority class rho : Overall system utilization (scalar)

Example

% Two priority classes lambda = [0.2; 0.3]; % Arrival rates mu = [1.0; 1.0]; % Service rates cs = [0.5; 1.0]; % Coefficients of variation [W, rho] = qsys_mg1_prio(lambda, mu, cs);

REFERENCE:

Kleinrock, L., “Queueing Systems, Volume I: Theory”, Wiley, 1975, Section 3.5

qsys_mg1_lrpt(lambda, mu, cs)

QSYS_MG1_LRPT Compute mean response time for M/G/1/LRPT queue

[W, RHO] = QSYS_MG1_LRPT(LAMBDA, MU, CS) computes the mean response time for each job class in an M/G/1 queue with Longest Remaining Processing Time (LRPT) scheduling.

Under LRPT, the job with the longest remaining processing time receives exclusive service. This is a remaining-size based policy that favors large jobs.

For LRPT, the slowdown for a job of size x is given by (Section 3.2 of Wierman-Harchol-Balter 2003):

E[S(x)]^LRPT = 1/(1-rho) + lambda*E[X^2]/(2*x*(1-rho)^2) E[T(x)]^LRPT = x * E[S(x)]^LRPT

For exponential service, the expected response time for class k is computed by integrating E[T(x)] over the exponential distribution of class k sizes.

CLASSIFICATION (Wierman-Harchol-Balter 2003):

LRPT is “Always Unfair” - it favors large jobs at the expense of small jobs.

Parameters:

• lambda – Vector of arrival rates per class

• mu – Vector of service rates per class

• cs – Vector of coefficients of variation per class (cs=1 for exponential)

Returns:

W – Vector of mean response times per class rho : Overall system utilization (modified for Little’s law)

References

• A. Wierman and M. Harchol-Balter, “Classifying scheduling policies with respect to unfairness in an M/GI/1”, SIGMETRICS 2003, Section 3.2.

Copyright (c) 2012-2026, Imperial College London All rights reserved.

qsys_mg1_fb(lambda, mu, cs)

QSYS_MG1_FB Compute mean response time for M/G/1/FB (Feedback/LAS) queue

[W, RHO] = QSYS_MG1_FB(LAMBDA, MU, CS) computes the mean response time for each job class in an M/G/1 queue with Feedback (FB) scheduling, also known as Least Attained Service (LAS) or Shortest Elapsed Time (SET).

Under FB/LAS, the job with the least attained service (smallest age) receives priority. This is an age-based policy where priority depends on how much service a job has received, not its original or remaining size.

For FB, the mean response time for a job of size x is given by (Section 3.3 of Wierman-Harchol-Balter 2003):

E[T(x)]^FB = (lambda * integral_0^x t*F_bar(t)dt) / (1-rho_x)^2

• x / (1 - rho_x)

where:

• rho_x = lambda * integral_0^x F_bar(t)dt (load from jobs completing

  service before reaching age x)

• F_bar(t) = 1 - F(t) (complementary CDF / survival function)

For exponential service, the expected response time for class k is computed by integrating E[T(x)] over the exponential distribution of class k sizes.

CLASSIFICATION (Wierman-Harchol-Balter 2003):

FB is “Always Unfair” - some job size is treated unfairly under all loads and all service distributions. However, FB approximates SRPT for heavy-tailed distributions and is practical since job sizes need not be known in advance.

Parameters:

• lambda – Vector of arrival rates per class

• mu – Vector of service rates per class

• cs – Vector of coefficients of variation per class (cs=1 for exponential)

Returns:

W – Vector of mean response times per class rho : Overall system utilization (modified for Little’s law)

References

• A. Wierman and M. Harchol-Balter, “Classifying scheduling policies with respect to unfairness in an M/GI/1”, SIGMETRICS 2003, Section 3.3.

• L. Kleinrock, “Queueing Systems, Volume II: Computer Applications”, Wiley, 1976.

Copyright (c) 2012-2026, Imperial College London All rights reserved.

qsys_mg1(lambda, mu, cs)

W=QSYS_MG1(LAMBDA,MU,CS)

qsys_mapph1(D0, D1, sigma, S, varargin)

QSYS_MAPPH1 Analyzes a MAP/PH/1 queue using BUTools MMAPPH1FCFS.

RESULT = QSYS_MAPPH1(D0, D1, SIGMA, S) analyzes a MAP/PH/1 queue with:

D0 - MAP hidden transition matrix (n x n) D1 - MAP arrival transition matrix (n x n) SIGMA - PH service initial probability vector (1 x m) S - PH service generator matrix (m x m)

RESULT = QSYS_MAPPH1(…, ‘numQLMoms’, K) computes K queue length moments RESULT = QSYS_MAPPH1(…, ‘numQLProbs’, N) computes N queue length probs RESULT = QSYS_MAPPH1(…, ‘numSTMoms’, K) computes K sojourn time moments

Returns a struct with fields:

meanQueueLength - Mean number of customers in system meanWaitingTime - Mean waiting time in queue meanSojournTime - Mean sojourn time (waiting + service) utilization - Server utilization queueLengthDist - Queue length distribution P(Q=n) queueLengthMoments- Raw moments of queue length sojournTimeMoments- Raw moments of sojourn time analyzer - Name of analyzer used

See also MMAPPH1FCFS, qsys_mapmap1(), qsys_phph1()

qsys_mapmap1(C0, C1, D0, D1, varargin)

QSYS_MAPMAP1 Analyzes a MAP/MAP/1 queue using BUTools MMAPPH1FCFS.

RESULT = QSYS_MAPMAP1(C0, C1, D0, D1) analyzes a MAP/MAP/1 queue with:

C0 - Arrival MAP hidden transition matrix (n x n) C1 - Arrival MAP arrival transition matrix (n x n) D0 - Service MAP hidden transition matrix (m x m) D1 - Service MAP observable transition matrix (m x m)

The service MAP is converted to an equivalent PH representation.

RESULT = QSYS_MAPMAP1(…, ‘numQLMoms’, K) computes K queue length moments RESULT = QSYS_MAPMAP1(…, ‘numQLProbs’, N) computes N queue length probs RESULT = QSYS_MAPMAP1(…, ‘numSTMoms’, K) computes K sojourn time moments

Returns a struct with fields:

meanQueueLength - Mean number of customers in system meanWaitingTime - Mean waiting time in queue meanSojournTime - Mean sojourn time (waiting + service) utilization - Server utilization queueLengthDist - Queue length distribution P(Q=n) queueLengthMoments- Raw moments of queue length sojournTimeMoments- Raw moments of sojourn time analyzer - Name of analyzer used

See also MMAPPH1FCFS, qsys_mapph1(), qsys_phph1()

qsys_mapm1(D0, D1, mu, varargin)

QSYS_MAPM1 Analyzes a MAP/M/1 queue using Q-MAM.

RESULT = QSYS_MAPM1(D0, D1, MU) analyzes a MAP/M/1 queue with:

D0 - MAP hidden transition matrix (n x n) D1 - MAP arrival transition matrix (n x n) MU - Exponential service rate

This is a convenience wrapper for QSYS_MAPMC with c=1.

RESULT = QSYS_MAPM1(…, ‘maxNumComp’, N) sets max queue length probs (default 500)

Returns a struct with fields:

meanQueueLength - Mean number of customers in system meanWaitingTime - Mean waiting time in queue meanSojournTime - Mean sojourn time (waiting + service) utilization - Server utilization queueLengthDist - Queue length distribution P(Q=n) waitingTimePH - Struct with alpha and T for waiting time PH analyzer - Name of analyzer used

See also QSYS_MAPMC, qsys_mapph1()

qsys_mapg1(D0, D1, serviceMoments, varargin)

QSYS_MAPG1 Analyzes a MAP/G/1 queue using BUTools MMAPPH1FCFS.

RESULT = QSYS_MAPG1(D0, D1, SERVICEMOMENTS) analyzes a MAP/G/1 queue with:

D0 - MAP hidden transition matrix (n x n) D1 - MAP arrival transition matrix (n x n) SERVICEMOMENTS - First k raw moments of service time [E[S], E[S^2], …]

(k = 2 or 3 for best accuracy)

The general service time is fitted to a Phase-Type distribution using moment matching before analysis.

RESULT = QSYS_MAPG1(…, ‘numQLMoms’, K) computes K queue length moments RESULT = QSYS_MAPG1(…, ‘numQLProbs’, N) computes N queue length probs RESULT = QSYS_MAPG1(…, ‘numSTMoms’, K) computes K sojourn time moments

Returns a struct with fields:

meanQueueLength - Mean number of customers in system meanWaitingTime - Mean waiting time in queue meanSojournTime - Mean sojourn time (waiting + service) utilization - Server utilization queueLengthDist - Queue length distribution P(Q=n) queueLengthMoments- Raw moments of queue length sojournTimeMoments- Raw moments of sojourn time analyzer - Name of analyzer used

See also MMAPPH1FCFS, APHFrom3Moments, qsys_mapph1()

qsys_mapmc(D0, D1, mu, c, varargin)

QSYS_MAPMC Analyzes a MAP/M/c queue using Q-MAM.

RESULT = QSYS_MAPMC(D0, D1, MU, C) analyzes a MAP/M/c queue with:

D0 - MAP hidden transition matrix (n x n) D1 - MAP arrival transition matrix (n x n) MU - Exponential service rate C - Number of servers

RESULT = QSYS_MAPMC(…, ‘maxNumComp’, N) sets max queue length probs (default 500)

Returns a struct with fields:

meanQueueLength - Mean number of customers in system meanWaitingTime - Mean waiting time in queue meanSojournTime - Mean sojourn time (waiting + service) utilization - Server utilization (per server) queueLengthDist - Queue length distribution P(Q=n) waitingTimePH - Struct with alpha and T for waiting time PH analyzer - Name of analyzer used

See also Q_CT_MAP_M_C, qsys_mapph1(), qsys_mapm1()

qsys_mapd1(D0, D1, s, varargin)

QSYS_MAPD1 Analyzes a MAP/D/1 queue using Q-MAM.

RESULT = QSYS_MAPD1(D0, D1, S) analyzes a MAP/D/1 queue with:

D0 - MAP hidden transition matrix (n x n) D1 - MAP arrival transition matrix (n x n) S - Deterministic service time (positive scalar)

RESULT = QSYS_MAPD1(…, ‘maxNumComp’, N) sets max queue length components (default 1000) RESULT = QSYS_MAPD1(…, ‘numSteps’, K) sets waiting time distribution granularity (default 1)

This is a convenience wrapper for qsys_mapdc with c=1.

See also qsys_mapdc(), Q_CT_MAP_D_C, qsys_mapm1()

qsys_is_retrial(sn)

Checks if network is a valid BMAP/PH/N/N bufferless retrial queue

Validates that the network structure matches the requirements for the BMAP/PH/N/N retrial queue solver: - Single bufferless queue (capacity == number of servers) - Retrial drop strategy configured - BMAP/MAP arrival process at source - PH/Exp service at queue - Open class model

Based on: Dudin et al., “Analysis of BMAP/PH/N-Type Queueing System with Flexible Retrials Admission Control”, Mathematics 2025, 13(9), 1434.

Parameters:

sn – Network structure

Returns:

isRetrial – True if network is valid BMAP/PH/N/N retrial topology retInfo: Struct with parameters for the retrial solver

Examples

[isRetrial, retInfo] = qsys_is_retrial(sn)

qsys_gm1(sigma, mu)

W=QSYS_GM1(SIGMA,MU)

qsys_gigk_approx_whitt(lambda, mu, ca, cs, k)

[W,RHOHAT]=QSYS_GIGK_APPROX_WHITT(LAMBDA,MU,CA,CS,K)

GI/G/k approximation of Whitt (1993), eqs. (2.16)-(2.25):

Wq = phi(rho,ca^2,cs^2,k) * ((ca^2+cs^2)/2) * Wq(M/M/k)

where phi interpolates the Cosmetatos M/D/k (phi1) and D/M/k (phi3) correction factors. Exact for M/M/k; reduces to the Cosmetatos M/D/k approximation for cs=0. Implements eq. (2.25) as printed, which was validated here against the paper’s Tables 5-7 (New column).

Inputs:

LAMBDA - Arrival rate MU - Service rate per server CA - Coefficient of variation of inter-arrival time CS - Coefficient of variation of service time K - Number of servers

Returns:

W - Average time in system (response time) RHOHAT - Modified utilization (so that M/M/1 formulas still hold)

Reference: Whitt, W. (1993). Approximations for the GI/G/m queue. Production and Operations Management 2(2), 114-161.

qsys_gigk_approx_kingman(lambda, mu, ca, cs, k)

W=QSYS_GIG1_UBND_KINGMAN(LAMBDA,MU,CA,CS,K)

qsys_gigk_approx_cosmetatos(lambda, mu, ca, cs, k)

[W,RHOHAT]=QSYS_GIGK_APPROX_COSMETATOS(LAMBDA,MU,CA,CS,K)

GI/G/k approximation by interpolation of the M/M/k, M/D/k and D/M/k queues (Cosmetatos 1982; Page 1982):

Wq = [ca^2*cs^2 + ca^2*(1-cs^2)*phi1/2 + (1-ca^2)*cs^2*phi3/2]*Wq(M/M/k)

where phi1 and phi3 are the Cosmetatos (1975) correction factors for M/D/k and D/M/k, with the safeguards of Whitt (1993). The D/D/k corner has Wq=0. The interpolation requires ca^2<=1 and cs^2<=1; outside this region the Lee-Longton scaling Wq = ((ca^2+cs^2)/2)*Wq(M/M/k) is used.

Inputs:

LAMBDA - Arrival rate MU - Service rate per server CA - Coefficient of variation of inter-arrival time CS - Coefficient of variation of service time K - Number of servers

Returns:

W - Average time in system (response time) RHOHAT - Modified utilization (so that M/M/1 formulas still hold)

References: Cosmetatos, G.P. (1975). Approximate explicit formulae for the average queueing time in the processes (M/D/r) and (D/M/r). INFOR 13, 328-331. Page, E. (1982). Tables of waiting times for M/M/n, M/D/n and D/M/n and their use to give approximate waiting times in more general queues. J. Opl. Res. Soc. 33, 453-473.

qsys_gigk_approx(lambda, mu, ca, cs, k)

W=QSYS_GIGK_APPROX(LAMBDA,MU,CA,CS,K)

qsys_gig1_ubnd_kingman(lambda, mu, ca, cs)

[W,RHOHAT]=QSYS_GIG1_UBND_KINGMAN(LAMBDA,MU,CA,CS)

Kingman’s upper bound on the mean waiting time of a G/G/1 queue:

Wq <= lambda*(sa^2+ss^2)/(2*(1-rho)), sa^2=ca^2/lambda^2, ss^2=cs^2/mu^2

The returned W adds the mean service time, so it upper-bounds the mean response time (time in system).

Reference: Kingman, J.F.C. (1962). Some inequalities for the queue GI/G/1. Biometrika 49(3/4), 315-324.

qsys_gig1_lbnd(lambda, mu, ca, cs)

[W,RHOHAT]=QSYS_GIG1_LBND(LAMBDA,MU,CA,CS)

Computes fundamental theoretical lower bounds for G/G/1 queues. These are the minimum possible values that performance measures cannot fall below for any realization of the arrival and service processes.

Inputs:

LAMBDA - Arrival rate MU - Service rate CA - Coefficient of variation of inter-arrival time CS - Coefficient of variation of service time

Returns:

W - Lower bound on average time in system (= 1/mu) RHOHAT - Modified utilization (so that M/M/1 formulas still hold)

qsys_gig1_rq(rho, mu, cs2, IaFun)

[Z,W,Q,X] = QSYS_GIG1_RQ(rho,mu,cs2,IaFun) - Robust Queueing (RQ) approximation for a single G/GI/1 queue partially characterized by its arrival rate, index of dispersion for counts (IDC) and the first two moments of the service time. Implements the mean steady-state workload

Z* = sup_{x>=0} { -(1-rho) x + sqrt( 2 rho x (I_a(x) + c2_s) / mu ) }

and the derived steady-state performance measures. Reference: W. Whitt and W. You (2018), “A Robust Queueing Network Analyzer Based on Indices of Dispersion”, eqs. (13),(16)-(18).

Inputs:

rho : traffic intensity lambda/mu (0<rho<1) mu : service rate cs2 : service SCV c2_s IaFun : handle, IaFun(x) -> arrival IDC I_a(x) at time argument x>0

Outputs:

Z : mean steady-state workload E[Z] W : mean steady-state waiting time E[W] Q : mean steady-state queue length E[Q] (number waiting + in service) X : mean number in system (= Q here for single class), kept for interface

qsys_gig1_approx_myskja2(lambda, mu, ca, cs, q0, qa)

[W,RHOHAT]=QSYS_GIG1_APPROX_MYSKJA2(LAMBDA,MU,CA,CS,Q0,QA)

Myskja’s enhanced third-moment approximation of the mean response time (time in system). For ca=1 the interpolation parameter theta is a 0/0 form, so the exact M/G/1 result is returned instead.

qa = third relative moment E[X^3]/6/E[X]^3, X=inter-arrival time r.v. q0 = lowest value of the relative third moment for a given mean and SCV

qsys_gig1_approx_marchal(lambda, mu, ca, cs)

W=QSYS_GIG1_APPROX_MARCHAL(LAMBDA,MU,CA,CS)

qsys_gig1_approx_kobayashi(lambda, mu, ca, cs)

W=QSYS_GIG1_APPROX_KOBAYASHI(LAMBDA,MU,CA,CS)

qsys_gig1_approx_klb(lambda, mu, ca, cs)

[W,rhohat]=QSYS_GIG1_APPROX_KLB(LAMBDA,MU,CA,CS)

qsys_gig1_approx_myskja(lambda, mu, ca, cs, q0, qa)

[W,RHOHAT]=QSYS_GIG1_APPROX_MYSKJA(LAMBDA,MU,CA,CS,Q0,QA)

Myskja’s third-moment approximation of the mean waiting time

Wq = rho/(2*mu*(1-rho))*((1+cs^2)+(q0/qa)^(1/rho-rho)*(1/rho)*(ca^2-1))

exact for M/G/1 (ca=1). The returned W adds the mean service time (response time, time in system).

qa = third relative moment E[X^3]/6/E[X]^3, X=inter-arrival time r.v. q0 = lowest value of the relative third moment for a given mean and SCV

qsys_gig1_approx_heyman(lambda, mu, ca, cs)

W=QSYS_GIG1_APPROX_HEYMAN(LAMBDA,MU,CA,CS)

qsys_gig1_approx_kimura(lambda, mu, ca, cs)

[W,RHOHAT]=QSYS_GIG1_APPROX_KIMURA(LAMBDA,MU,CA,CS)

Kimura’s diffusion-interpolation approximation of the mean waiting time

Wq = rho*(ca^2+cs^2)/(mu*(1-rho)*(1+ca^2))

exact for M/M/1 and M/G/1. The returned W adds the mean service time (response time, time in system).

Reference: Kimura, T. (1986). A two-moment approximation for the mean waiting time in the GI/G/s queue. Management Science 32(6), 751-763.

qsys_gig1_approx_allencunneen(lambda, mu, ca, cs)

W=QSYS_GIG1_APPROX_ALLENCUNNEEN(LAMBDA,MU,CA,CS)

qsys_gig1_approx_gelenbe(lambda, mu, ca, cs)

[W,RHOHAT]=QSYS_GIG1_APPROX_GELENBE(LAMBDA,MU,CA,CS)

Gelenbe’s diffusion approximation with instantaneous-return boundary:

p(0) = 1-rho, p(n) = rho*(1-rhat)*rhat^(n-1), n>=1 rhat = exp(-2*(1-rho)/(rho*ca^2+cs^2))

hence E[N] = rho/(1-rhat) and the mean response time (time in system) W = E[N]/lambda = 1/(mu*(1-rhat)).

Reference: Gelenbe, E. (1975). On approximate computer system models. Journal of the ACM 22(2), 261-269.

qsys_dmc(lambda_arr, mu, c, varargin)

QSYS_DMC Analyzes a D/M/c queue (deterministic interarrivals, exp service).

RESULT = QSYS_DMC(LAMBDA, MU, C) returns time-average performance metrics by embedding the system at arrival epochs and integrating over the inter-arrival cycle. The death-only sub-generator A[m,m-1] = min(m,c)*MU encodes service completions; X_{n+1} = expm(A*s)[X_n+1, :].

Optional name-value parameters:

‘truncation’ - state-space truncation (default auto) ‘quadSteps’ - trapezoidal-rule steps over a cycle (default 200)

RESULT struct fields:

meanQueueLength - time-average E[N] meanWaitingQueue - time-average Lq = E[(N-c)+] meanWaitingTime - Wq = Lq / lambda meanSojournTime - Wq + 1/MU utilization - rho = lambda / (c*mu) analyzer - identifier

See also qsys_mapdc(), qsys_mdc_crommelin

qsys_bmapphnn_retrial(D, beta, S, N, alpha, gamma, p, R, varargin)

QSYS_BMAPPHNN_RETRIAL Analyzes a BMAP/PH/N/N bufferless retrial queue.

RESULT = QSYS_BMAPPHNN_RETRIAL(D, BETA, S, N, ALPHA, GAMMA, P, R) analyzes a BMAP/PH/N/N bufferless retrial queueing system with admission control.

This implements the algorithm from: Dudin et al., “Analysis of BMAP/PH/N-Type Queueing System with Flexible Retrials Admission Control”, Mathematics 2025, 13(9), 1434.

Inputs:

D - Cell array {D0, D1, …, DK} of BMAP matrices

D0: hidden transition matrix (V x V) D1, …, DK: arrival matrices for batches of size 1, …, K

BETA - PH service initial probability vector (1 x M) S - PH service subgenerator matrix (M x M) N - Number of servers (also capacity, hence bufferless) ALPHA - Retrial rate per customer in orbit GAMMA - Impatience (abandonment) rate per customer in orbit P - Probability of batch rejection when not enough servers R - Admission threshold (scalar or 1 x V vector per BMAP state)

When n > R(nu), arriving customers go to orbit

Optional parameters:

‘MaxLevel’ - Fixed orbit truncation level. When empty or non-positive

(default) the level is chosen adaptively: it is doubled until the mass retained at the top level contributes less than ‘TailTolerance’ of the mean orbit length. A fixed level disables the adaptive refinement.

‘Tolerance’ - Convergence tolerance (default: 1e-10) ‘TailTolerance’- Relative orbit-truncation error target (default: 1e-6) ‘MaxDim’ - Cap on the total generator dimension explored by the

adaptive refinement (default: 2e5)

‘MaxBlockSize’ - Cap on the per-level block size V*d (default: 5000).

Exceeding it is an error: the phase-type service order and the server count make the level block intractable.

‘Verbose’ - Print progress messages (default: false)

Returns a struct with fields:

L_orbit - Mean number of customers in orbit N_server - Mean number of busy servers L_system - Mean number in system (orbit + servers) Utilization - Server utilization (N_server / N) Throughput - System throughput P_idle - Probability all servers are idle P_empty_orbit - Probability orbit is empty P_empty_system - Probability system is empty (idle and empty orbit) pi - Stationary distribution (levels x Vd) truncLevel - Truncation level used truncError - Relative orbit-truncation error estimate at truncLevel analyzer - Name of analyzer used

Example

% M/M/3/3 retrial queue (exponential arrivals and service) D = {-2.0, 2.0}; % Exp(2) arrivals beta = 1; S = -1; % Exp(1) service N = 3; alpha = 0.5; % Retrial rate gamma = 0; % No impatience p = 0; % No batch rejection R = 2; % Admission threshold result = qsys_bmapphnn_retrial(D, beta, S, N, alpha, gamma, p, R);

See also qsys_mapph1(), qsys_is_retrial()

Copyright (c) 2012-2026, Imperial College London All rights reserved.