"""
Modified-offered-load and pointwise-stationary approximations for time-varying
multiserver systems.
Native Python twin of matlab/src/api/qsys/qsys_mtgs0_mol.m, implementing the
approximation analyzed by W. A. Massey and 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, and the pointwise stationary
approximation of W. Whitt (1991), Management Science 37(3), 307-314.
"""
from typing import Any, Callable, Dict, Optional, Sequence
import numpy as np
from .mtginf import qsys_mtginf
[docs]
def erlang_b(s: int, a: float) -> float:
"""
Erlang B blocking probability with ``s`` servers and offered load ``a``, by
the recursion ``B_j = a B_{j-1}/(j + a B_{j-1})``, which never forms
``a^s/s!`` and so never overflows.
Args:
s: number of servers
a: offered load in erlangs
Returns:
The probability that all servers are busy.
"""
b = 1.0
for j in range(1, int(s) + 1):
b = a * b / (j + a * b)
return b
[docs]
def erlang_c(s: int, a: float) -> float:
"""
Erlang C delay probability with ``s`` servers and offered load ``a``, from
the same recursion; 1 when the load saturates the servers.
Args:
s: number of servers
a: offered load in erlangs
Returns:
The probability that an arrival waits.
"""
if a >= s:
return 1.0
b = erlang_b(s, a)
rho = a / s
return b / (1.0 - rho * (1.0 - b))
[docs]
def qsys_mtgs0_mol(lambdaFun: Callable[[Any], Any], serviceCcdf: Callable[[Any], Any],
ES: float, s: int, tvals: Sequence[float],
startTime: float = -np.inf, delay: bool = False,
ES2: Optional[float] = None, **kwargs) -> Dict[str, Any]:
"""
Modified-offered-load (MOL) and pointwise-stationary (PSA) approximations for
a time-varying multiserver system.
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. PSA uses the instantaneous one, ``lambda(t)E[S]``.
MOL uses the offered load of the corresponding INFINITE-SERVER system,
m(t) = E[S] E[lambda(t - S_e)] = 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 that 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.
WHY IT WORKS. The blocking system differs from the infinite-server one only
in what happens at the ceiling, and the ceiling does not change the AGE
structure of the load much when blocking is not extreme. That is why the
approximation is asymptotically correct in the many-server regime and
degrades when blocking is heavy.
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: that is dominated by the peak of the cycle,
where both approximations are weakest. Under constant input MOL is exact,
reducing to the stationary Erlang formula.
Args:
lambdaFun: the arrival rate; must accept arguments in the past when
``startTime`` is infinite
serviceCcdf: G^c(x) = P(S > x)
ES: the mean service time
s: number of servers
tvals: times at which to evaluate
startTime: time the system started empty; -inf assumes an infinite past
delay: use Erlang C rather than Erlang B, i.e. approximate the DELAY
probability of an Mt/M/s queue rather than the blocking probability
of an Mt/G/s/0 loss system
ES2: second moment of the service time, passed through for the time lag
**kwargs: passed to :func:`qsys_mtginf`
Returns:
Dict with ``times``, ``offeredLoad`` (m(t)), ``instantLoad``
(lambda(t)E[S]), ``probBlockMOL``, ``probBlockPSA``, ``meanBusyMOL``
(the carried load ``m(t)(1-B)`` for the loss model), and, when ``ES2``
is given, ``meanLag``.
References:
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.
"""
s = int(round(s))
if s < 1:
raise ValueError('The number of servers s must be at least 1.')
inf_server = qsys_mtginf(lambdaFun, serviceCcdf, ES, tvals, startTime=startTime, ES2=ES2,
**kwargs)
t = inf_server['times']
m = inf_server['meanNumber']
inst = inf_server['offeredLoadPSA']
f = erlang_c if delay else erlang_b
mol = np.array([f(s, float(v)) for v in m])
psa = np.array([f(s, float(v)) for v in inst])
result: Dict[str, Any] = {
'times': t,
'offeredLoad': m,
'instantLoad': inst,
'probBlockMOL': mol,
'probBlockPSA': psa,
'meanBusyMOL': m * (1.0 - mol) if not delay else np.minimum(m, float(s)),
'arrivalRate': inf_server['arrivalRate'],
}
if 'meanLag' in inf_server:
result['meanLag'] = inf_server['meanLag']
return result