LN: A Meta-solver for Layered Queueing Network Analysis

Abstract

We overview LN, a novel solver introduced in the LINE software package to analyze layered queueing network (LQN) models. The novelty of the LN solver lies in its capability to analyze LQNs with a user-defined combination of solution paradigms, including discrete-event and stochastic simulation, continuous-time Markov chain analysis (CTMC), normalizing constant evaluation (NC), matrix analytic methods (MAM), mean-field approximations (FLUID), and mean-value analysis (MVA). Being parametric in the solver used for each LQN layer, LN as a whole enables the efficient computation of advanced performance metrics such as marginal and joint state probabilities, response and passage time distributions, and transient measures, leveraging individual strengths of the supported solution paradigms. We discuss in particular recent developments added to NC, the default layer solver of LN, which significantly improve the solution of queueing network models obtained using loose layering of the LQN.

Appendices

A Proof of Theorem 1

For a homogeneous model, the CoMoM recurrence relation may be written as

$$ \boldsymbol{A} \varLambda (N) = \boldsymbol{B} \varLambda (N-1_R) $$

where

$$ \boldsymbol{A} = \begin{bmatrix} \boldsymbol{A}_{1,1} &{} \boldsymbol{A}_{1,2}\\ \boldsymbol{0} &{} \boldsymbol{A}_{2,2}\\ \end{bmatrix} \quad \boldsymbol{B} = \begin{bmatrix} \boldsymbol{B}_{1,1} &{} \boldsymbol{0}\\ \boldsymbol{B}_{2,1} &{} \boldsymbol{B}_{2,2}\\ \end{bmatrix} $$

Let $\boldsymbol{0}_{I,J}$ indicate a block of zeros of size $I\times J$. Defining $\tilde{D}=(D_1,\ldots ,D_R)^T$, we have

$$ \boldsymbol{A}_{1,1} = \begin{bmatrix} 1 &{} -\tilde{D}^T \\ \boldsymbol{0}_{R-1,1} &{} -m\mathop {\textrm{diag}}(\tilde{D}) \end{bmatrix} \quad \boldsymbol{A}_{1,2} = \begin{bmatrix} -1 &{} \boldsymbol{0}_{1,R-1} \\ \tilde{N} &{} -\mathop {\textrm{diag}}(\tilde{Z}) \end{bmatrix} \quad \boldsymbol{A}_{2,2} = N_R \boldsymbol{I} $$

$$ \boldsymbol{B}_{1,1} = D_R \boldsymbol{E}_{1,1} \quad \boldsymbol{B}_{2,1} = mD_R \boldsymbol{I} \quad \boldsymbol{B}_{2,2} = Z_R \boldsymbol{I} $$

The inverse of the block upper triagonal matrix A is now computed as

$$ \boldsymbol{A}^{-1} = \begin{bmatrix} \boldsymbol{A}_{1,1}^{-1} &{} ~~\boldsymbol{S} \boldsymbol{A}_{2,2}^{-1}\\ \boldsymbol{0}_{R-1,1} &{} ~~\boldsymbol{A}_{2,2}^{-1} \\ \end{bmatrix} $$

with $\boldsymbol{S}=-\boldsymbol{A}_{1,1}^{-1}\boldsymbol{A}_{1,2}$. Observe first that

$$ \boldsymbol{A}^{-1}_{1,1}= \begin{bmatrix} 1 &{} -m^{-1}e^T\\ \boldsymbol{0}_{R-1,1} &{} -m^{-1}\mathop {\textrm{diag}}(\tilde{D})^{-1} \\ \end{bmatrix} \quad \boldsymbol{A}^{-1}_{2,2}=N_R^{-1} \boldsymbol{I} $$

where $e^T=\tilde{D}^T\mathop {\textrm{diag}}(\tilde{D})^{-1}=(1,\ldots ,1)$. Thus

$$\begin{aligned} \boldsymbol{S}&= -\begin{bmatrix} 1 &{} -m^{-1}e^T\\ \boldsymbol{0}_{R-1,1} &{} -m^{-1}\mathop {\textrm{diag}}(\tilde{D})^{-1} \end{bmatrix} \begin{bmatrix} -1 &{} \boldsymbol{0}_{1,R-1} \\ \tilde{N} &{} -\mathop {\textrm{diag}}(\tilde{Z}) \end{bmatrix} \\ {}&= -m^{-1} \begin{bmatrix} - |\tilde{N}|-m &{} ~~\tilde{Z}^T \\ -\mathop {\textrm{diag}}(\tilde{D})^{-1} \tilde{N} &{} ~~\mathop {\textrm{diag}}(\tilde{D})^{-1}\mathop {\textrm{diag}}(\tilde{Z}) \end{bmatrix} \end{aligned}$$

Note that $\boldsymbol{A}_{2,2}$ is the only block that depends on $N_R$. We can therefore write

$$\begin{aligned} \boldsymbol{A}^{-1}\boldsymbol{B}&=\boldsymbol{F}_{1,R}+N_R^{-1}\boldsymbol{F}_{2,R} \end{aligned}$$

where

$$\begin{aligned} \boldsymbol{F}_{1,R}=\begin{bmatrix} \boldsymbol{A}^{-1}_{1,1} &{} \boldsymbol{0}\\ \boldsymbol{0} &{} \boldsymbol{0}\\ \end{bmatrix} \begin{bmatrix} D_R \boldsymbol{E}_{1,1} &{} \boldsymbol{0}\\ mD_R \boldsymbol{I} &{} Z_R \boldsymbol{I}\\ \end{bmatrix} =\begin{bmatrix} D_R \boldsymbol{E}_{1,1} &{} \boldsymbol{0}\\ \boldsymbol{0} &{} \boldsymbol{0}\\ \end{bmatrix} \end{aligned}$$

$$\begin{aligned} \boldsymbol{F}_{2,R}= \begin{bmatrix} \boldsymbol{0} &{} \boldsymbol{S}\\ \boldsymbol{0} &{} \boldsymbol{I}\\ \end{bmatrix} \begin{bmatrix} D_R \boldsymbol{E}_{1,1} &{} \boldsymbol{0}\\ mD_R \boldsymbol{I} &{} Z_R \boldsymbol{I}\\ \end{bmatrix} = \begin{bmatrix} mD_R \boldsymbol{S} &{} Z_R \boldsymbol{S}\\ mD_R \boldsymbol{I} &{} Z_R \boldsymbol{I}\\ \end{bmatrix} \end{aligned}$$

B Gaussian Quadratures

A Gauss-Laguerre quadrature of order K evaluates exponentially-weighted integrals by means of the approximation

$$\begin{aligned} \int _{x=0}^\infty e^{-x}f(x) dx \approx \sum _{k=1}^K w_k f(x_k) \end{aligned}$$

(4)

where $x_k$ denotes the k-th root of the Laguerre polynomial

$$ L_K(x) = \sum _{j=0}^K \left( {\begin{array}{c}K\\ i\end{array}}\right) \frac{(-1)^j}{j!} x^j $$

and with weights $w_k = {x_k} \left( {\left( k + 1\right) ^2 \left[ L_{k+1}\left( x_k\right) \right] ^2}\right) ^{-1}$.

Gauss-Legendre methods are similar but applicable to finite ranges [a, b]. Setting $a=0$ and large b they can also help evaluating the normalizing constant. Their main benefit is that nodes and weights do not incur the same floating-point range exceptions as observed instead for Gauss-Laguerre quadratures of large order. We point to [20] for further details on Gauss-Legendre methods.

C Proof of Theorem 2

For a homogeneous model with m identical single-server stations, the McKenna-Mitra integral takes the form

$$ G(m,N)=\frac{1}{\prod _{r=1}^R N_r!} \int _{u_1=0}^{+\infty } \cdots \int _{u_m=0}^{+\infty } e^{-(u_1+\ldots +u_m)} h(u_1+\ldots +u_m) du_1\cdots du_m $$

where $h(u)=\prod _{r=1}^R \left( Z_r + D_r u\right) ^{N_r}$. We note that the multidimensional integral may be interpreted as computing $E[h(U_1+\ldots +U_m)]$ for the i.i.d. exponential random variables $U_i\sim Exp(1)$. The result then readily follows after noting that $U_1+\ldots +U_m$ is Erlang-m distributed with density $f(u)=\frac{1}{(m-1)!}u^{m-1}e^{-u}$.

D Proof of Theorem 3

Let the entries of $\varPi (N)$ be indicated with $\widetilde{\pi }_N(n)$, $n=|N|,\ldots ,0$. A probabilistic population constraint holds for homogeneous models with $m=1$ [6, Thm. 6]

$$ N_R \widetilde{\pi }_N(n) = Z_R \widetilde{\pi }_{N-1_R}(n) + nD_R \widetilde{\pi }_{N-1_R}(n-1) $$

for all $n=1,\ldots ,|N|$ and where $\widetilde{\pi }_{N-1_R}(n-1)=0$ if $n=0$. With a load-dependent queueing station ($m=1$), the derivation in [6] generalizes with similar passages to the following form

$$\begin{aligned} N_R \widetilde{\pi }_N(n) = Z_R \widetilde{\pi }_{N-1_R}(n) + n\frac{D_R}{\mu (n)} \widetilde{\pi }_{N-1_R}(n-1) \end{aligned}$$

(5)

where $\mu (n)$ is the load-dependent scaling at the queueing station. Organizing (5) in matrix form, we get (3)

$$ \boldsymbol{T}_R = \begin{bmatrix} \quad Z_R\quad &{} |N|\dfrac{D_R}{\mu (|N|)} &{} 0 &{} \cdots &{} 0 \\ 0 &{} Z_R\quad \quad &{} (|N|-1)\dfrac{D_R}{\mu (|N|-1)} \quad &{} \cdots &{} 0 \\ \vdots &{} \vdots &{} \vdots &{} \ddots &{} \vdots \\ 0 &{} 0 &{} 0 &{} Z_R\quad &{} \dfrac{D_R}{\mu (1)} \\ 0 &{} 0 &{} 0 &{} 0 &{} Z_R \end{bmatrix} $$

As assumed, consider now an homogeneous layer, where there are m identical load-independent single-server stations. The proof follows by noting that, if $m>1$, the m queueing stations can be exactly replaced by a flow-equivalent server station with identical $D_1$,..., $D_R$ and [26]

$$ \mu (n)=\frac{n}{(n+m-1)} $$

The final expression for $\boldsymbol{T}_{R}$ follows after plugging the above expression for $\mu (n)$.

E Software Architecture Design

Figure 6 illustrates the key architectural elements of LINE, including the NetworkStruct data structure, and the Network, NetworkSolver, LayeredNetwork, EnsembleSolver and LayeredNetworkSolver classes.

Fig. 6.

Key architectural elements of LINE.

The Network object summarizes the model characteristics and acts as its persistence layer. The object is generated by the user either through a domain-specific language offered by LINE  [9] or via model-to-model transformations from other formats (e.g., JMT’s XML [2], PMIF [32]). Besides the model specification, a Network object can cache the model state space, its initial state, and retain information needed for the traffic equations in state-dependent models.

Each Network object is able to synthesize via the getStruct method a NetworkStruct data structure. The latter includes key model parameters, such as representations of service and arrival processes, job populations, and number of servers, among others. In addition, the data structure includes the routing table, the associated chains, and the average number of visits that each class pays to each node. NetworkStruct also offers indexing functions, that allow for example to differentiate between stations, where jobs can reside, and nodes, which are elements of the network traversed with zero service time (e.g., a fork).

The NetworkSolver object encodes a solver type, of which the aforementioned six LINE solvers are specific instances. The main role of this class is to ensure consistent computation of performance results, adopting identical conventions for reporting per-class and per-chain results to the end-user. Operational relationships are also applied by this object to derive certain performance metrics from the ones returned by the solvers, e.g. arrival rates from throughputs [14].

Each NetworkSolver object is equipped with a getAvgTable method that returns mean performance metrics for the model in a tabular format. The method invokes via the runAnalyzer method one of the solution methods offered by that solver, which operates solely on the NetworkStruct data structure. Model transformations that alter the model topology are conducted within runAnalyzer. An example is tagging a job class, which is used in response time distribution analysis.

The EnsembleSolver specifies the life-cycle for an iterative solution method that works on an ensemble of Network objects. This class allows to bind a NetworkSolver to each particular Network in the ensemble, applying consistently actions before, in-between, and after each iteration, and verifying convergence. It also harmonizes the presentation of ensemble-level results to the end-user. Besides LN, the Env solver is another example of EnsembleSolver, wherein iteration is used to analyze random environments [13].

The LN solver is a special instance of EnsembleSolver, operating on an ensemble consisting of the LQN layers. The LayeredNetwork class encompasses the objects that form an LQN, such as the Entry, Task, Host, and Activity classes.

The LayeredNetwork class offers a getEnsemble method that generates, and stores within the LayeredNetwork object, the ensemble of Network models, each mapping to a distinct LQN layer. Similarly to Network, this class also exposes a getStruct method that builds a static data structure of the LQN parameters.

The LN solver, implemented in the SolverLN class, is specified parametrically in terms of any of the LINE solvers, or a custom combination thereof. For example, the end-user may require to use LINE ’s simulators on layers that include non-Markovian service distributions (e.g., Pareto) and MVA otherwise.