Li et al.: MAP Parameter Estimation of QBD Queueing Systems with Utilization Data

Nomenclature↩︎

Acronyms and Abbreviations

Notations

1 Introduction↩︎

Parameter estimation of queueing systems is a long-standing topic, which has been extensively studied in past decades [1]–[4]. Typically, a queueing system can be symbolized as \(A/B/S/K\), where \(A\) represents the arrival distribution, \(B\) represents the service distribution, \(S\) is the number of service nodes, and \(K\) is the capacity of the queueing system. For example, a queueing model where the arrival process follows a Poisson distribution, the service time follows an exponential distribution, the number of service nodes is 1, and the queue length is \(K\) can be represented as \(M/M/1/K\). For a Poisson process, the renewal process requires a sequence of independent and identically distributed non-negative random variables formed between inter-arrival times.

A Markovian arrival process (MAP) is a generalization of the Poisson process with dependence between the inter-arrival times and non-exponential inter-arrival time distribution [5]. MAP was first proposed by Neuts [6] and has been widely used to analyze mathematically stochastic behaviors such as reliability, teletraffic, and performance evaluation [7], [8]. MAP is defined as a specific continuous-time Markov chain (CTMC) to represent the arrivals of a queueing system. A MAP consists of phase and level processes with corresponding discrete state spaces, respectively. The phase process is represented by a CTMC, which reflects the dynamic process of the internal state. The level process represents the number of events, which can be a counting process such as a Poisson process. MAP usually has two main properties. (i) MAP is the most versatile stochastic counting process, and it contains many arrival processes, such as the Poisson process, Markov-modulated Poisson process (MMPP) [9], and batch MAP (BMAP) [10]. In addition, MAP is known to be dense [11] for any stochastic point process, which can approximate complex stochastic counting processes. (ii) MAP can represent the time correlation in arrival streams. The commonly used arrival processes, such as the homogeneous Poisson process (HPP), often assume that the inter-arrivals are independent. However, this assumption is unrealistic in practice since the inter-arrival times are usually correlated, such as the long-range dependence problem in communication networks [12]. In actual Ethernet traffic, long-range dependence problems are caused by self-similar processes. Therefore, MAP-based queueing system modeling has become significant.

Generally, the parameters (e.g., arrival and service rates) can be estimated by observing the behavior of a queueing system during a specified time interval. For example, all the above literature collected inter-arrival times to estimate the parameters of the arrival stream. To date, moment-based and likelihood-based approaches are the two most widely used approaches to parameter estimation for MAP. The moment-based approach determines the MAP parameters to fit the theoretical moments to empirical moments from observed data [13], [14]. The likelihood-based approach aims to find the MAP parameters from the empirical data when their likelihood is maximum [15].

However, collecting and recording inter-arrival times is costly and unrealistic in real-world situations. Additional data that record the available information of queueing systems (e.g., waiting times and queue lengths) are easier to obtain that are used for parameter estimation [16], [17]. However, such data are difficult to collect in a computer system from observation of a queue. Utilization data is defined as the ratio of the time interval during which the system is busy to a fixed monitoring time. In a computer system, utilization data such as CPU utilization is one of the most commonly used statistics to monitor CPU behavior during task execution. Most operating systems have the function to calculate CPU utilization by default. Unfortunately, utilization data belongs to incomplete data that only can be collected during observable time intervals, while the information in the unobservable periods is missing. Although utilization data is easy to obtain, exact inter-arrivals and waiting times are not known by observation. To the best of our knowledge, no studies have been conducted to estimate MAP parameters from utilization data due to the limitations of traditional parameter estimation approaches.

An expectation-maximization (EM) algorithm is an iterative way to compute maximum likelihood estimation (MLE) from incomplete data [18]–[20]. Therefore, an EM algorithm is used to estimate the MAP parameters from utilization data. To further generalize the proposed approach, this study seeks to estimate the parameters of a MAP quasi-birth-death (QBD) queueing system. The main contributions of this study are summarized below:

• Data collection of utilization data: unlike traditional observable data, the arrival process is unknown from utilization data. Additionally, utilization data contains unobservable and observable intervals, and only utilization with the observable interval can be collected.

• Generalization of a queueing system via MAP and QBD: MAP is a generalization of the Poisson process to represent more complex bursts and associated traffic streams. QBD is a generalization of the birth-death process and is a powerful tool for modeling many stochastic phenomena.

• An EM algorithm-based parameter estimation approach: unlike the traditional moment-based and MLE-based approaches, a nontrivial EM algorithm can compute the MLE for utilization data containing unobservable intervals.

2 Related Works↩︎

A temporal point process is a stochastic process used to model a time series of binary events that occur at random intervals in continuous time [21]. A continuous-time stochastic process can be modeled by a CTMC with a discrete state space.CTMC models have applications in many fields, such as traffic modeling [22], [23], earthquake prediction [24], and performance evaluation [25], [26].

2.1 Parameter Estimation of MAP↩︎

MAP is one of the most flexible stochastic processes and is defined as a specific CTMC. Since MAP can approximate any point process, it is often used to model general arrival and service processes in a queueing system [27]–[30]. Bruneo et al. [27] proposed a virtualized system with a regenerative policy to evaluate the performance of the arrivals following the Poisson process and MMPP. Experiments show that the proposed workload-based policy is superior to the time-based strategy. Zheng et al. [28] considered the workload-based and time-based policies for arrival streams following MAP. They also analyzed the loss probability of the transactions and the average response time of the estimated parameters based on Markovian arrivals. Klimenok et al. [29] estimated the parameters of a queueing system in which the arrival stream and service time followed MAP and exponential distribution, respectively. They assumed that the system has two servers and that the buffer size is infinite. Vygovskaya et al. [30] estimated the parameters of the arrival process with a simple variant of MAP (i.e., MMPP) as the waiting times. Their queueing system can be represented as \(MMPP/M/2\). Note that the collected time series data are observable.

2.2 Parameter Estimation for Observable Data↩︎

Generally, moment-based and likelihood-based estimations are the two major approaches to fit the observed data and MAP parameter estimation. The moment-based estimation approach can reduce computational costs compared with the likelihood-based estimation.

In moment estimation, the parameters of MAP can be estimated by matching the observed data between the sample moments and the population moments [31]. Heffes et al. [32] provided an explicit formula to estimate the parameters of the two-state MMPP (i.e., MMPP(2)) from empirical moments by observing the number of arrivals. Anderson et al. [33] estimated the superposition parameters of the two-state MAP (i.e., MAP(2)) based on moments by collecting the number of arrivals of the switched Poisson process. Yoshihara et al. [14] fitted the superposition of MMPP(2) using a moment estimation approach for self-similar traffic. Mitchell et al. developed an approximate model to estimate the parameters of the \(G/G/1/N\) queueing system by observing inter-arrival times. Additionally, Telek et al. [34] considered a model from the moment inter-arrival time distribution to estimate the parameters of MAP.

The likelihood-based estimation is used to find the optimal MAP parameters by maximizing the likelihood function from a series of observations [35], [36]. Although MLE is a well-known method for parameter estimation of a queueing system, it has significant limitations for estimating the parameters of MAP. MAP consists of many parameters, which require large matrix operations. Therefore, MLE cannot work well with an increasing number of phases in MAP. Baum et al. [37] proposed the forward-backward algorithm for MAP in the statistical estimation of the probabilistic functions for Markov chains, which is the earliest paper for the EM algorithm. Deng et al. [38] used an EM algorithm to convert the MMPP from continuous time to a discrete-time domain, obtaining an MLE containing the model parameters. Breuer [39] designed a specification of the classical EM algorithm for the communication systems of MAP and BMAP. Roberts et al. [40] developed a scaling forward-backward algorithm to improve Rydeń’s EM algorithm to estimate the parameters of MMPP. Furthermore, Buchholz [41], [42] proposed a two-step EM algorithm for fitting the PH distribution and HMM parameters based on real traffic observations. Horvat́h et al. [43] presented a two-step MAP fitting method. The first step fitted the PH distribution with the inter-arrival times. The second step approximated the lag correlation values. Okamura et al. [44] improved the two-step fitting algorithm for MAP with inter-arrival time data.

2.3 Parameter Estimation for Utilization Data↩︎

Almost all the above studies estimated the parameters of queueing systems from observable data such as arrival times and waiting times. However, such parameter estimation approaches cannot apply to group and utilization data. Due to the greater complexity of these data types, their parameters are difficult to estimate. Okamura et al. [18] proposed an EM algorithm for estimating the parameters of MAP for group data. Group data is defined as a group of arrival times, with each group being one bin. In each observed time interval, the number of arrival times is collected. They proposed a novel EM algorithm to calculate the expected log-likelihood function (LLF) based on the group data and then maximize the LLF to estimate the parameters. Li et al. [45] proposed the MLE approach and used the utilization data to estimate the parameters of the queueing system \(M_t/M/1/K\), where the arrival process followed a nonhomogeneous Poisson process (NHPP). Since the arrival intensity of the NHPP is a continuous function of time \(t\), parameter estimation is difficult. To solve the problems, they approximated the NHPP with a piecewise constant function that converted the NHPP into a series of HPPs.

This study presents an EM algorithm to estimate the parameters of a QBD queueing system with arrivals following MAP from utilization data. Compared to previous work [45], the queueing system that follows MAP arrivals is more general. A nontrivial EM algorithm is proposed to make the estimation approach work well on the queueing system. Moreover, QBD is a generalization of the birth-death process and is a powerful tool for modeling stochastic processes. The proposed EM algorithm can estimate the MAP parameters for QBD queueing systems with utilization data.

3 Preliminaries↩︎

In this section, some preliminary knowledge on MAP, MMPP, and QBD process are presented, which is necessary to understand this work.

3.1 Markovian Arrival Process↩︎

MAP is a remarkably versatile modeling tool in point process theory. Usually, the arrival rate of MAP can be governed by CTMC. Formally, let \(\boldsymbol{D}_0\) denote the infinitesimal generator of the underlying CTMC where no arrivals occur, and \(\boldsymbol{D}_1\) denote a rate matrix that triggers a change of state of the CTMC. Then, the infinitesimal generator matrix of the CTMC is defined by \(\boldsymbol{Q}_{MAP}\) as follows: \[\begin{align} \boldsymbol{Q}_{MAP} = \begin{pmatrix} \boldsymbol{D}_0&\boldsymbol{D}_1\\ &\boldsymbol{D}_0&\boldsymbol{D}_1\\ &&\ddots&\ddots\\ &&&\boldsymbol{D}_0&\boldsymbol{D}_1\\ \end{pmatrix}. \end{align}\]

A \(m\)-state MAP can be denoted as MAP(\(m\)). Let \(q_{i,j(i \neq j)}\) denote the transition rate from phase \(i\) to phase \(j\), where \(q_{i,i} = \sum_{j=1,j \neq i}^m q_{i,j} + \sum_{j=1}^m \lambda_{i,j}\). Then, the \(\boldsymbol{D}_0\) and \(\boldsymbol{D}_1\) matrices can be represented as follows: \[\begin{align} \boldsymbol{D}_0 &= \begin{pmatrix} -q_{1,1} & q_{1,2} & \cdots & q_{1,m}\\ q_{2,1} & -q_{2,2} & \cdots & q_{2,m}\\ \vdots & \vdots & \ddots & \vdots\\ q_{m,1}\;&q_{m,2} & \cdots & -q_{m,m} \\ \end{pmatrix}, \label{eq:d0} \end{align}\tag{1}\] \[\begin{align} \boldsymbol{D}_1 &= \begin{pmatrix} \lambda_{1,1} & \lambda_{1,2} & \cdots & \lambda_{1,m} \\ \lambda_{2,1} & \lambda_{2,2} & \cdots & \lambda_{2,m} \\ \vdots & \vdots & \ddots & \vdots\\ \lambda_{m,1} & \lambda_{m,2} & \cdots & \lambda_{m,m} \\ \end{pmatrix}. \end{align}\]

Let \(\{N(t); t \geq 0\}\) and \(\{J(t); t \geq 0\}\) represent the cumulative number of arrivals during the time interval \([0, t)\) and the phase at time \(t\), respectively. In this study, \(N(t)\) and \(J(t)\) are called level and phase, respectively. The infinitesimal generator of the phase process \(J(t)\) can be represented by \(\boldsymbol{D}_0+\boldsymbol{D}_1\) with the following properties: \[\begin{align} (\boldsymbol{D}_0+\boldsymbol{D}_1)\boldsymbol{1}=\boldsymbol{0} \end{align}\] where \(\boldsymbol{1}\) and \(\boldsymbol{0}\) are two column vectors with all elements of 1 and 0, respectively. For a MAP(\(m\)), let \(\boldsymbol{\pi}_s\) be the stationary probability vector and \(\boldsymbol{\pi}_s=(\pi_{s_1},\pi_{s_2},...,\pi_{s_m})\). Then, \[\begin{align} \boldsymbol{\pi}_s(\boldsymbol{D}_0+\boldsymbol{D}_1)=\boldsymbol{0}, \quad \boldsymbol{\pi}_s \boldsymbol{1}=1. \end{align}\] Define the fundamental arrival rate as \(\tilde{\lambda}\). \(\tilde{\lambda}\) can be represented by the average number of arrivals, which is given by \[\begin{align} \tilde{\lambda}=\boldsymbol{\pi}_s\boldsymbol{D}_1\boldsymbol{1}. \end{align}\] Define the matrix \(\boldsymbol{P}_k(t)\) whose \((i, j)\)-element is given by \[\begin{align} [\boldsymbol{P}_k(t)]_{i,j}=P(N(t)=k, J(t)=j | N(0)=0, J(0)=i). \end{align}\] Then, the differential-difference equations can be computed by \[\begin{align} \frac{d}{dt}\boldsymbol{P}_0(t)=\boldsymbol{P}_0(t)\boldsymbol{D}_0, \end{align}\] \[\begin{align} \frac{d}{dt}\boldsymbol{P}_k(t)=\boldsymbol{P}_k(t)\boldsymbol{D}_0 + \boldsymbol{P}_{k-1}(t)\boldsymbol{D}_1, \quad k=1,2,3,.... \end{align}\]

MMPP is a doubly stochastic process where the arrival rate of a Poisson process is modulated by a CTMC, hence the name. The MMPP can be identified as a special case of MAP, where \(\boldsymbol{D}_0\) and \(\boldsymbol{D}_1\) of the MMPP are defined by Eq. (1 ) and a diagonal matrix, whose representation is \[\begin{align} \boldsymbol{D}_1=\text{diag} \{\lambda_1,...,\lambda_m\} \end{align}\] In a word, the main difference between MAP and MMPP is whether the phase change after an arrival.

3.2 Quasi-birth-death Process↩︎

In general, a QBD is a special case of infinite-state CTMCs. A continuous-time QBD Markov process can be defined by a two-dimensional process \(\{N(t), J(t)\}\). The \(level\) \(N(t)\) may have a finite or infinite number of states, \(K\). Assume that there is only one server in a queueing system with the finite capacity \(K (\geq 1)\). Also, the arrivals of the queueing system are served according to the FCFS discipline. In the case of no arrival in the system, the transaction process for an arrival starts immediately. In the case that there are less than \(K\) arrivals in the system, an arrival waits to be served until all previous arrivals completing service in the system. In the case that there are \(K\) arrivals in the system, the next arrivals are refused. Then the infinitesimal generator matrix \(\boldsymbol{Q}_{QBD}\) of the Markov process can take the \((K+1)\)-by-\((K+1)\) block tridiagonal form as follows: \[\begin{align} \boldsymbol{Q}_{QBD} = \left( \begin{array}{c|ccccc} \boldsymbol{B}_0 & \boldsymbol{A}_0 & & & & \\ \hline \boldsymbol{A}_2 & \boldsymbol{A}_1 & \boldsymbol{A}_0 & & & \\ & \ddots & \ddots & \ddots & & \\ & & \boldsymbol{A}_2 & \boldsymbol{A}_1 & \boldsymbol{A}_0 & \\ & & & & \boldsymbol{A}_2 & \boldsymbol{A}_1 +\boldsymbol{A}_0 \\ \end{array}\right) \label{eq:qbd} \end{align}\tag{2}\] where \(\boldsymbol{B}_0\) indicates the case of state transitions at idle period, \(\boldsymbol{A}_0\) represents the state transitions with arrivals, \(\boldsymbol{A}_1\) is the state transitions with no arrival, and \(\boldsymbol{A}_2\) expresses the state transitions during the job service process.

For example, consider a specific queueing system, whose arrival process follows a MAP, and the service time follows an exponential distribution with the service rate \(\mu\). The queueing system can be indicated by the symbol \(MAP/M/1/K\), and the infinitesimal generator matrix for the queueing system can be expressed by: \[\begin{align} \boldsymbol{Q}_{QBD} = \begin{pmatrix} \boldsymbol{D}_0 & \boldsymbol{D}_1 & \\ \mu \boldsymbol{I}& \boldsymbol{D}_0-\mu \boldsymbol{I}& \boldsymbol{D}_1 & \\ & \ddots & \ddots & \\ & & \mu \boldsymbol{I}& \boldsymbol{D}_0+\boldsymbol{D}_1-\mu \boldsymbol{I}\\ \end{pmatrix} \label{eq:q95qbd} \end{align}\tag{3}\]

3.3 Brief Overview of the EM Algorithm↩︎

The EM algorithm is considered as an iterative method to calculate the MLE from incomplete data [18]–[20]. Define \({\mathcal{D}}\) and \({\mathcal{U}}\) as the observable data and unobservable data, respectively. Then, a set of model parameters \(\boldsymbol{\theta}\) are estimated according to the given observable data \({\mathcal{D}}\). In general, the EM algorithm consists of two steps:

where \(E_{{\mathcal{U}}}\) is the expectation operator for the unobservable data \({\mathcal{U}}\). Eq. (4 ) provides an updated formula for the parameters. That is, the estimated parameters of the M-step are used in the next E-step. The parameters \(\boldsymbol{\theta}\) are updated iteratively with the two steps until the LLF or parameters converge.

4 Parameter Esitmation for MAP↩︎

In this section, we consider the parameter estimation for a QBD process from utilization data. More precisely, an MLE method via the EM algorithm is considered to estimate MAP parameters.

4.1 Utilization Data Format↩︎

Here, we consider an EM algorithm with utilization data. Utilization data is a kind of time-series data, which is usually defined as the time fraction of busy time over the time length of observable periods. To estimate MAP parameters only from the utilization data, the following assumptions are given.

(1)

Utilization data can be monitored at every time interval.

(2)

Each time interval consists of two successive periods: unobservable time interval and observable time interval.

(3)

There is only one state change at most in every observable time interval. That is, a state change from busy (idle) to idle (busy).

The first two assumptions are reasonable according to the definition of utilization data. Also, the third assumption is reasonable that at most only one state change can occur in every observable period when the monitoring time of the utilization data is sufficiently small.

Formally, we define a series of samples of the utilization data in a monitoring time as \({\mathcal{D}}=(u_1, u_2, ..., u_N)\), where \(u_n\) is the \(n\)th utilization sample of the observable period, and \(0 \le u_n \le 1\). Fig. 1 demonstrates an possible behavior of system state. In the figure, let \(t_u\) and \(t_o\) represent the two successive unobservable time interval and observable time interval, and \(t_o \ll t_u\). Notice that the numbers of samples are counted only in the observable time interval, while the numbers of samples in the unobservable time intervals are unknown. In other words, the sample data in the unobservable time interval are missing. Besides, let \(B_{t_o}^{(n)}\) and \(I_{t_o}^{(n)}\) be the busy time and idle time for the \(n\)th observable time interval. Then, the utilization is formulated by \[\begin{align} u_n=\frac{B_{t_o}^{(n)}}{B_{t_o}^{(n)}+I_{t_o}^{(n)}} \end{align}\]

4.2 A QBD Queueing System with Markovian Arrivals↩︎

Consider a QBD queueing system, in which the jobs arrive at the system following MAP. To formulate the MLE via the EM algorithm, the infinitesimal generator \(\boldsymbol{Q}_{QBD}\) in Eq. (2 ) are divided into four sub-matrices according to the two solid lines: \[\begin{align} &\boldsymbol{Q}_{00}=\boldsymbol{B}_0=\boldsymbol{D}_0, \\ &\boldsymbol{Q}_{01}=\left(\boldsymbol{A}_0, \boldsymbol{O}, ..., \boldsymbol{O}\right)=\left(\boldsymbol{D}_1, \boldsymbol{O}, ..., \boldsymbol{O}\right),\\ &\boldsymbol{Q}_{10}^T=\left(\boldsymbol{A}_2, \boldsymbol{O}, ..., \boldsymbol{O}\right), \end{align}\] \[\begin{align} \boldsymbol{Q}_{11} &= \begin{pmatrix} \boldsymbol{A}_1 & \boldsymbol{A}_0 & & \\ \boldsymbol{A}_2 & \boldsymbol{A}_1 & \boldsymbol{A}_0 & \\ & \ddots & \ddots & \\ & \boldsymbol{A}_2 & \boldsymbol{A}_1 & \boldsymbol{A}_0 \\ & & \boldsymbol{A}_2 & \boldsymbol{A}_1+\boldsymbol{A}_0 \\ \end{pmatrix} \end{align}\] where \(\boldsymbol{O}\) represents \(m \times m\) zero matrix, and \(\boldsymbol{Q}_{10}^T\) indicates the matrix transpose operation of \(\boldsymbol{Q}_{10}\). Moreover, since \(\boldsymbol{Q}\) is a \((K+1)\)-by-\((K+1)\) block matrix, and each block is a \(m \times m\) matrix, the dimensions of matrix \(\boldsymbol{Q}_{00}, \boldsymbol{Q}_{01}, \boldsymbol{Q}_{10}\) and \(\boldsymbol{Q}_{11}\) are \(m \times m, m \times mK, mK \times m\) and \(mK \times mK\), respectively.

Here, \(\boldsymbol{Q}_{00}\) indicates that the system is an idle state without the occurrence of state change between idle and busy, while \(\boldsymbol{Q}_{01}, \boldsymbol{Q}_{10}\) and \(\boldsymbol{Q}_{11}\) represent that the system state changes from idle to busy, busy to idle, and busy to busy, respectively.

4.3 M-Step Formulas↩︎

Consider an \(m\)-state MAP with utilization data \({\mathcal{D}}\). Also, for the sake of notational convenience, we define the following several unobservable variables and observable variables.

\(B_i\)

an indicator random variable for the event that the pahse is \(i\) at the initial time \(t=0\).

\(Z_i^{[n,l]}\)

cumulative sojourn time that the number of jobs is \(l\) in phase \(i\) during the \(n\)th unobservable time interval.

\(N_{i,j}^{[n,l]}\)

the number of phase transitions from phase \(i\) to phase \(j\) that the number of jobs is \(l\) during the \(n\)th unobservable time interval.

\(A_{i,j}^{[n,l]}\)

the number of arrivals leading to phase transitions from phase \(i\) to phase \(j\) that the number of jobs is \(l\) during the \(n\)th unobservable time interval.

\(S_{i,j}^{[n,l]}\)

the number of services leading to phase transitions from phase \(i\) to phase \(j\) that the number of jobs is \(l\) during the \(n\)th unobservable time interval.

\({\tilde{Z}}_i^{[n,l]}\)

cumulative sojourn time that the number of jobs is \(l\) in phase \(i\) during the \(n\)th observable time interval.

\({\tilde{N}}_{i,j}^{[n,l]}\)

the number of phase transitions from phase \(i\) to phase \(j\) that the number of jobs is \(l\) during the \(n\)th observable time interval.

\({\tilde{A}}_{i,j}^{[n,l]}\)

the number of arrivals leading to phase transitions from phase \(i\) to phase \(j\) that the number of jobs is \(l\) during the \(n\)th observable time interval.

\({\tilde{S}}_{i,j}^{[n,l]}\)

the number of services leading to phase transitions from phase \(i\) to phase \(j\) that the number of jobs is \(l\) during the \(n\)th observable time interval.

Define a vector of parameters \(\boldsymbol{\theta}:=\{\pi_{i}, q_{i,j}, \lambda_{i,j}, \mu_{i,j}\}\), the unobservable variables \({\mathcal{U}}:=\{B_i, Z_i^{[n,l]}, N_{i,j}^{[n,l]}, A_{i,j}^{[n,l]}, S_{i,j}^{[n,l]}\}\) and the observable variables \({\mathcal{D}}:=\{{\tilde{Z}}_i^{[n,l]}, {\tilde{N}}_{i,j}^{[n,l]}, {\tilde{A}}_{i,j}^{[n,l]}, {\tilde{S}}_{i,j}^{[n,l]}\}\), for \(i,j=1,2,...,m\) and \(n=1,..,N\). Then the MLEs of the parameters can be obtained from Eq. (4 ): \[\begin{align} \label{eq:pi} \pi_i:\text{E}[B_i|{\mathcal{D}}] \end{align}\tag{5}\] \[\begin{align} \label{eq:q} q_{i,j}:=\frac{\sum_{n=1}^{N}\text{E}\left[N_{i,j}^{[n,l]}+{\tilde{N}}_{i,j}^{[n,l]}|{\mathcal{D}}\right]}{\sum_{n=1}^{N}\text{E}\left[Z_i^{[n,l]}+{\tilde{Z}}_i^{[n,l]}|{\mathcal{D}}\right]},\quad i\neq j \end{align}\tag{6}\] \[\begin{align} \label{eq:lam} \lambda_{i,j}:=\frac{\sum_{n=1}^{N}\text{E}\left[A_{i,j}^{[n,l]}+{\tilde{A}}_{i,j}^{[n,l]}\right]}{\sum_{n=1}^{N}\text{E}\left[Z_i^{[n,l]}+{\tilde{Z}}_i^{[n,l]}|{\mathcal{D}}\right]} \end{align}\tag{7}\] \[\begin{align} \label{eq:mu} \mu_{i,j}:=\frac{\sum_{n=1}^{N}\text{E}\left[S_{i,j}^{[n,l]}+{\tilde{S}}_{i,j}^{[n,l]}|{\mathcal{D}}\right]}{\sum_{n=1}^{N}\text{E}\left[Z_i^{[n,l]}+{\tilde{Z}}_i^{[n,l]}|{\mathcal{D}}\right]} \end{align}\tag{8}\]

Note that the subscript of the expectation operations is omitted for the sake of simplicity. In addition, all the expected values in the above Eq. (5 )-(8 ) are calculated based on the previous M-step. Also, the EM algorithm provides an initial guess of MAP parameters at the initial step.

4.4 E-Step Formulas↩︎

In the E-step, the analytical forms of the expected values in Eq. (5 )-(8 ) are derived. Define \(T\) as the length of monitoring time interval, and \(T=t_u+t_o\). Then the cumulative time sequence for utilization data can be indicated by \(s_0=0<s_1<...<s_N\), i.e., \(s_n=n \times T, (n=1,2,...,N)\). Let \({\mathcal{A}}\) be the following event: \[\begin{align} {\mathcal{A}}_n=\{N(s_n^+) - N(s_n^-) = a_n\} \end{align}\] where \(s_n^-\) and \(s_n^+\) represent the left limit and right limit, and \[\begin{align} &\{N(s_n^-)=x, N(s_n^+)=y\} \nonumber\\ &=\lim_{\Delta t \to +0}\{N(s_n-\Delta t)=x, N(s_n+\Delta t)=y\} \end{align}\] Then the forward, backward and overall events can be indicated by \({\mathcal{F}}_n={\mathcal{A}}_1...{\mathcal{A}}_n\), \({\mathcal{B}}_n={\mathcal{A}}_n...{\mathcal{A}}_N\) and \({\mathcal{O}}={\mathcal{A}}_1...{\mathcal{A}}_N\), respectively. For simplicity, let \(P(A)\) represent the probability of an indicator random variable \(A\).

Define \(\boldsymbol{f}(n), \tilde{\boldsymbol{f}}(n)\) and \(\tilde{\boldsymbol{f}^{\prime}}(n)\) as \(1 \times m(K+1)\) row vectors, who represent the probabilities (likelihoods) for the forward events in the \(n\)th period. Also, define \(\boldsymbol{b}(n), \tilde{\boldsymbol{b}}(n)\) and \(\tilde{\boldsymbol{b}^{\prime}}(n)\) as \(m(K+1) \times 1\) column vectors, who represent the probabilities (likelihoods) for the backward events in the \(n\)th period. Specifically, \(\boldsymbol{f}(n), \boldsymbol{b}(n)\) and \(\tilde{\boldsymbol{f}}(n), \tilde{\boldsymbol{b}}(n)\)represent the probability vectors at the beginning of the unobservable time interval and observable time interval, respectively. \(\tilde{\boldsymbol{f}^{\prime}}(n), \tilde{\boldsymbol{b}^{\prime}}(n)\) are the probability vectors of the system state change between busy and idle in the obserable time interval.

For a better understanding, fig. 2 shows some examples to explain the probability vectors of forward and backward events.

Notice that each of the six vectors is partitioned by \(levels\) into subvectors. In other words, each vector has \(1 \times (K+1)\) block subvector with dimension \(1 \times m\). Depending on whether the system state is idle or busy, the six vectors can be re-divided into two parts, which are demonstrated as follows: events. \[\begin{align} \label{eq:fb01} \begin{aligned} &\boldsymbol{f}(n)=(\boldsymbol{f}_0(n), \boldsymbol{f}_1(n)), \quad \boldsymbol{b}(n)=(\boldsymbol{b}_0(n), \boldsymbol{b}_1(n)) \\ &\boldsymbol{\tilde{\boldsymbol{f}}}(n)=(\tilde{\boldsymbol{f}}_0(n), \tilde{\boldsymbol{f}}_1(n)), \quad \boldsymbol{\tilde{b}}(n)=(\tilde{\boldsymbol{b}}_0(n), \tilde{\boldsymbol{b}}_1(n)) \\ &\boldsymbol{\tilde{f}}^{\prime} (n)=(\tilde{\boldsymbol{f}}_0^{\prime} (n), \tilde{\boldsymbol{f}}_1^{\prime} (n)), \quad \boldsymbol{\tilde{b}}^{\prime} (n)=(\tilde{\boldsymbol{b}}_0^{\prime} (n), \tilde{\boldsymbol{b}}_1^{\prime} (n)) \end{aligned} \end{align}\tag{9}\] where the two elements of every vector represent the probabilities when the utilization state being idle and busy in the \(n\)th monitoring period, and the dimensions are \(1 \times m\) and \(1 \times mK\), respectively.

Consider the indicator random variable \({\mathcal{O}}\), we can formulate the Eq. (5 ) as \[\begin{align} \label{eq:pi95i} \pi_i:=\frac{\text{E}\left[B_i{\mathcal{O}}\right]}{P({\mathcal{O}})}=\frac{\pi_i [\boldsymbol{b}(1)]_i}{\boldsymbol{\pi} \boldsymbol{b}(1)} \end{align}\tag{10}\] where \(\boldsymbol{\pi}\) be the initial probability vector, and \(\boldsymbol{\pi}=(\pi_1,\pi_2,...,\pi_m)\) with \(\sum_{i=1}^{m}\pi_i=1\).

Next, according to the Eq. (6 ), we consider the expected values \(\text{E}\left[Z_i^{[n,l]}|{\mathcal{D}}\right], \text{E}\left[\tilde{Z}_i^{[n,l]}|{\mathcal{D}}\right], \text{E}\left[N_{i,j}^{[n,l]}|{\mathcal{D}}\right]\) and \(\text{E}\left[\tilde{N}_{i,j}^{[n,l]}|{\mathcal{D}}\right]\), respectively. For the calculation of \(\text{E}\left[Z_i^{[n,l]}|{\mathcal{D}}\right]\), since \(\text{E}\left[Z_i^{[n,l]}|{\mathcal{D}}\right]= \text{E}\left[Z_i^{[n,l]}{\mathcal{O}}\right]/P({\mathcal{O}})\), the subsequent analysis treats only \(\text{E}\left[Z_i^{[n,l]}{\mathcal{O}}\right]\), which can be formulated by \[\begin{align} &\text{E}\left[Z_i^{[n,l]}{\mathcal{O}}\right]=\nonumber\\ &\int_{0}^{t_u} \left[\boldsymbol{f}(n)\exp\left(\boldsymbol{Q}s\right)\right]_{(l,i)}\left[\exp\left(\boldsymbol{Q}(t_u-s)\right)\tilde{\boldsymbol{b}}(n)\right]_{(l,i)}ds \end{align}\] Note that \([\cdot]_{(l,i)}\) indicates the \(i\)th element of the \(l\)th block in the block vector \([\cdot]\).

Similarly, \(\text{E}\left[\tilde{Z}_i^{[n,l]}|{\mathcal{D}}\right]\) can be derived by \[\begin{align} \text{E}\left[\tilde{Z}_i^{[n,l]}|{\mathcal{D}}\right]=\frac{\text{E}\left[\tilde{Z}_i^{[n,l]}{\mathcal{O}}\right]}{P({\mathcal{O}})} \end{align}\] For \(\text{E}\left[\tilde{Z}_i^{[n,l]}{\mathcal{O}}\right]\), the calculations are given by the two equations as follows according to the number of jobs \(l\).

if \(l=0\): \[\begin{align} &\text{E}\left[\tilde{Z}_i^{[n,l]}{\mathcal{O}}\right]= \int_{0}^{(1-u_n)t_o} \left[\tilde{\boldsymbol{f}}_0(n)\exp\left(\boldsymbol{Q}_{00} s\right)\right]_i\nonumber\\ &\times \left[\exp\left(\boldsymbol{Q}_{00}\left((1-u_n)t_o-s\right)\right)\boldsymbol{Q}_{01}\tilde{\boldsymbol{b}}_1^{\prime}(n)\right]_i ds \nonumber\\ &+ \int_{0}^{(1-u_n)t_o}\left[\tilde{\boldsymbol{f}}_1^{\prime}(n)\boldsymbol{Q}_{10}\exp\left(\boldsymbol{Q}_{00} s\right)\right]_i \nonumber\\ &\times \left[\exp\left(\boldsymbol{Q}_{00}\left((1-u_n)t_o-s\right)\right){\boldsymbol{b}}_0(n+1)\right]_i ds \end{align}\] where \([\cdot]_i\) indicates the \(i\)th element of vector \([\cdot]\).

if \(l>0\): \[\begin{align} &\text{E}\left[\tilde{Z}_i^{[n,l]}{\mathcal{O}}\right]=\int_{0}^{u_n t_o} \left[\tilde{\boldsymbol{f}}_1(n)\exp\left(\boldsymbol{Q}_{11} s\right)\right]_{(l,i)}\nonumber\\ &\times \left[\exp\left(\boldsymbol{Q}_{11}(u_n t_o-s)\right)\boldsymbol{Q}_{10}\tilde{\boldsymbol{b}}_0^{\prime}(n)\right]_{(l,i)}ds \nonumber\\ &+ \int_{0}^{u_n t_o}\left[\tilde{\boldsymbol{f}}_0^{\prime}(n)\boldsymbol{Q}_{01}\exp\left(\boldsymbol{Q}_{11} s\right)\right]_{(l,i)}\nonumber\\ &\times \left[\exp\left(\boldsymbol{Q}_{11}(u_n t_o-s)\right){\boldsymbol{b}}_1(n+1)\right]_{(l,i)}ds \end{align}\]

Similarly, \(\text{E}\left[N_{i,j}^{[n,l]}|{\mathcal{D}}\right]\) can be derived by \[\begin{align} \text{E}\left[N_{i,j}^{[n,l]}|{\mathcal{D}}\right]=\frac{\text{E}\left[N_{i,j}^{[n,l]}{\mathcal{O}}\right]}{P({\mathcal{O}})} \end{align}\] \[\begin{align} &\text{E}\left[N_{i,j}^{[n,l]}{\mathcal{O}}\right]= \int_{0}^{t_u} \left[\boldsymbol{f}(n) \exp\left(\boldsymbol{Q}s\right)\right]_{(l,i)}[\boldsymbol{B}_0]_{i,j}\nonumber\\ &\times \left[\exp\left(\boldsymbol{Q}(t_u-s)\right)\tilde{\boldsymbol{b}}(n)\right]_{(l,j)}ds \end{align}\] where \([\cdot]_{i,j}\) represents the \((i,j)\)th element of the matrix \([\cdot]\).

In addition, \(\text{E}\left[\tilde{N}_{i,j}^{[n,l]}\right]\) can be derived by \[\begin{align} \text{E}\left[\tilde{N}_{i,j}^{[n,l]}|{\mathcal{D}}\right]=\frac{\text{E}\left[\tilde{N}_{i,j}^{[n,l]}{\mathcal{O}}\right]}{P({\mathcal{O}})} \end{align}\] According to the number of jobs \(l\), the computation of \(\text{E}\left[\tilde{N}_{i,j}^{[n,l]}{\mathcal{O}}\right]\) are as follows.

if \(l=0\): \[\begin{align} &\text{E}\left[\tilde{N}_{i,j}^{[n,l]}{\mathcal{O}}\right]= \int_{0}^{(1-u_n)t_o} \left[\tilde{\boldsymbol{f}}_0(n) \exp\left(\boldsymbol{Q}_{00} s\right)\right]_i [\boldsymbol{B}_0]_{i,j}\nonumber\\ &\times \left[\exp\left(\boldsymbol{Q}_{00}\left((1-u_n)t_o-s\right)\right)\boldsymbol{Q}_{01}\tilde{\boldsymbol{b}}_1^{\prime}(n)\right]_j ds \nonumber\\ &+ \int_{0}^{(1-u_n)t_o}\left[\tilde{\boldsymbol{f}}_1^{\prime}(n)\boldsymbol{Q}_{10}\exp\left(\boldsymbol{Q}_{00} s\right)\right]_i [\boldsymbol{B}_0]_{i,j}\nonumber\\ &\times \left[\exp\left(\boldsymbol{Q}_{00}\left((1-u_n)t_o-s\right)\right){b}_0(n+1)\right]_j ds \end{align}\]

if \(l>0\): \[\begin{align} &\text{E}\left[\tilde{N}_{i,j}^{[n,l]}{\mathcal{O}}\right]= \int_{0}^{u_n t_o} \left[\tilde{\boldsymbol{f}}_1(n) \exp\left(\boldsymbol{Q}_{11} s\right)\right]_{(l,i)}[\boldsymbol{B}_0]_{i,j}\nonumber\\ &\times \left[\exp\left(\boldsymbol{Q}_{11}(u_n t_o-s)\right)\boldsymbol{Q}_{10}\tilde{\boldsymbol{b}}_0^{\prime}(n)\right]_{(l,j)}ds \nonumber\\ &+ \int_{0}^{u_n t_o}\left[\tilde{\boldsymbol{f}}_0^{\prime}(n)\boldsymbol{Q}_{01}\exp\left(\boldsymbol{Q}_{11} s\right)\right]_{(l,i)}[\boldsymbol{B}_0]_{i,j}\nonumber\\ &\times \left[\exp\left(\boldsymbol{Q}_{11}(u_n t_o-s)\right){\boldsymbol{b}}_1(n+1)\right]_{(l,j)}ds \end{align}\]

Next, according to the Eq. (7 ), we consider the expected values \(\text{E}\left[A_{i,j}^{[n,l]}|{\mathcal{D}}\right]\) and \(\text{E}\left[\tilde{A}_{i,j}^{[n,l]}|{\mathcal{D}}\right]\), respectively. \(\text{E}\left[A_{i,j}^{[n,l]}|{\mathcal{D}}\right]\) can be derived by \[\begin{align} \text{E}\left[A_{i,j}^{[n,l]}|{\mathcal{D}}\right]=\frac{\text{E}\left[A_{i,j}^{[n,l]}{\mathcal{O}}\right]}{P({\mathcal{O}})} \end{align}\] \[\begin{align} &\text{E}\left[A_{i,j}^{[n,l]}{\mathcal{O}}\right]= \int_{0}^{t_u} \left[\boldsymbol{f}(n) \exp\left(\boldsymbol{Q}s\right)\right]_{(l,i)}[\boldsymbol{A}_0]_{i,j}\nonumber\\ &\times \left[\exp\left(\boldsymbol{Q}(t_u-s)\right)\tilde{\boldsymbol{b}}(n)\right]_{(l+1,j)}ds \end{align}\]

Similarly, \(\text{E}\left[\tilde{A}_{i,j}^{[n,l]}\right]\) is derived by \[\begin{align} \text{E}\left[\tilde{A}_{i,j}^{[n,l]}|{\mathcal{D}}\right]=\frac{\text{E}\left[\tilde{A}_{i,j}^{[n,l]}{\mathcal{O}}\right]}{P({\mathcal{O}})} \end{align}\]

if \(l=0\): \[\begin{align} &\text{E}\left[\tilde{A}_{i,j}^{[n,l]}{\mathcal{O}}\right]= \int_{0}^{(1-u_n)t_o} \left[\tilde{\boldsymbol{f}}_0(n) \exp\left(\boldsymbol{Q}_{00} s\right)\right]_i [\boldsymbol{A}_0]_{i,j}\nonumber\\ &\times \left[\exp\left(\boldsymbol{Q}_{00}\left((1-u_n)t_o-s\right)\right)\boldsymbol{Q}_{01}\tilde{\boldsymbol{b}}_1^{\prime}(n)\right]_j ds \end{align}\]

if \(l>0\): \[\begin{align} &\text{E}\left[\tilde{A}_{i,j}^{[n,l]}{\mathcal{O}}\right]= \int_{0}^{u_n t_o} \left[\tilde{\boldsymbol{f}}_1(n)\exp\left(\boldsymbol{Q}_{11} s\right)\right]_{(l,i)}[\boldsymbol{A}_0]_{i,j}\nonumber\\ &\times \left[\exp\left(\boldsymbol{Q}_{11}(u_n t_o-s)\right)\boldsymbol{Q}_{10}\tilde{b}_0^{\prime}(n)\right]_{(l+1,j)}ds \nonumber\\ &+ \int_{0}^{u_n t_o}\left[\tilde{\boldsymbol{f}}_0^{\prime}(n)\boldsymbol{Q}_{01}\exp\left(\boldsymbol{Q}_{11} s\right)\right]_{(l,i)}[\boldsymbol{A}_0]_{i,j}\nonumber\\ &\times \left[\exp\left(\boldsymbol{Q}_{11}(u_n t_o-s)\right){b}_1(n+1)\right]_{(l+1,j)}ds \end{align}\]

Finally, we calculate \(\text{E}\left[S_{i,j}^{[n,l]}|{\mathcal{D}}\right]\) and \(\text{E}\left[\tilde{S}_{i,j}^{[n,l]}|{\mathcal{D}}\right]\) according to the Eq. (8 ). \(\text{E}\left[S_{i,j}^{[n,l]}|{\mathcal{D}}\right]\) is derived by: \[\begin{align} \text{E}\left[S_{i,j}^{[n,l]}|{\mathcal{D}}\right]=\frac{\text{E}\left[S_{i,j}^{[n,l]}{\mathcal{O}}\right]}{P({\mathcal{O}})} \end{align}\] Then, \(\text{E}\left[S_{i,j}^{[n,l]}{\mathcal{O}}\right]\) can be expressed by \[\begin{align} &\text{E}\left[S_{i,j}^{[n,l]}{\mathcal{O}}\right]= \int_{0}^{t_u} \left[\boldsymbol{f}(n)\exp\left(\boldsymbol{Q}s\right)\right]_{(l,i)}[\boldsymbol{A}_2]_{i,j}\nonumber\\ &\times \left[\exp\left(\boldsymbol{Q}(t_u-s)\right)\tilde{\boldsymbol{b}}(n)\right]_{(l-1,j)}ds \end{align}\]

Similarly, \(\text{E}\left[\tilde{S}_{i,j}^{[n,l]}|{\mathcal{D}}\right]\) derived by \[\begin{align} \text{E}\left[\tilde{S}_{i,j}^{[n,l]}|{\mathcal{D}}\right]=\frac{\text{E}\left[\tilde{S}_{i,j}^{[n,l]}{\mathcal{O}}\right]}{P({\mathcal{O}})} \end{align}\] According to the number of jobs \(l\), the equation \(\text{E}\left[\tilde{S}_{i,j}^{[n,l]}{\mathcal{O}}\right]\) can be divided into two parts according to \(l\):

if \(l=1\): \[\begin{align} &\text{E}\left[\tilde{S}_{i,j}^{[n,l]}{\mathcal{O}}\right]= \int_{0}^{(1-u_n)t_o} \left[\tilde{\boldsymbol{f}}_1(n)\exp\left(\boldsymbol{Q}_{11} s\right)\right]_{(l,i)} [\boldsymbol{A}_2]_{i,j}\nonumber\\ &\times \left[\exp\left(\boldsymbol{Q}_{11}\left((1-u_n)t_o-s\right)\right)\boldsymbol{Q}_{10}\tilde{\boldsymbol{b}}_0^{\prime}(n)\right]_{(l-1,j)} ds \end{align}\]

if \(l>1\): \[\begin{align} \label{eq:S95ij} &\text{E}\left[\tilde{S}_{i,j}^{[n,l]}{\mathcal{O}}\right]= \int_{0}^{u_n t_o} \left[\tilde{\boldsymbol{f}}_1(n)\exp\left(\boldsymbol{Q}_{11} s\right)\right]_{(l,i)}[\boldsymbol{A}_2]_{i,j}\nonumber\\ &\times \left[\exp\left(\boldsymbol{Q}_{11}(u_n t_o-s)\right)\boldsymbol{Q}_{10}\tilde{\boldsymbol{b}}_0^{\prime}(n)\right]_{(l-1,j)}ds \nonumber\\ &+ \int_{0}^{u_n t_o}\left[\tilde{\boldsymbol{f}}_0^{\prime}(n)\boldsymbol{Q}_{01}\exp\left(\boldsymbol{Q}_{11} s\right)\right]_{(l,i)}[\boldsymbol{A}_2]_{i,j}\nonumber\\ &\times \left[\exp\left(\boldsymbol{Q}_{11}(u_n t_o-s)\right){\boldsymbol{b}}_1(n+1)\right]_{(l-1,i)}ds \end{align}\tag{11}\]

4.5 Computation of E-step Formulas↩︎

The E-step of MAP estimation with utilization data requires the computation of the probabilities of the forward and backward events and their convolutions.

According to fig. 2, we first derive the likelihoods for forward events in the \(n\)th period. Since the state transitions in unobservable intervals can not be observed, \(\tilde{\boldsymbol{f}}(n)\) is computed by \[\begin{align} \label{eq:hatf} \tilde{\boldsymbol{f}}(n)=\boldsymbol{f}(n)\exp\left(\boldsymbol{Q}t_o\right) \end{align}\tag{12}\] In the \(n\)th observable interval, according to the idle state and busy state of the system, the computation of the likelihood \(\tilde{\boldsymbol{f}}^{\prime}(n)\) are divided into two parts, which are expressed as: \[\begin{align} \tilde{\boldsymbol{f}}_0^{\prime} (n)=\tilde{\boldsymbol{f}}_0(n) \exp\left(\boldsymbol{Q}_{00}(1-u_n) t_o\right) \end{align}\] \[\begin{align} \tilde{\boldsymbol{f}}_1^{\prime} (n)=\tilde{\boldsymbol{f}}_1(n) \exp\left(\boldsymbol{Q}_{11}u_n t_o\right) \end{align}\] Also, the instantaneous transition probability \(\boldsymbol{f}(n+1)\) can be calculated based on \(\tilde{\boldsymbol{f}}^{\prime} (n)\). According to Eq. 9 , \(\boldsymbol{f}_0(n+1)\) and \(\boldsymbol{f}_1(n+1)\) are shwon as follows:

if \(u_n=0\): \[\begin{align} \boldsymbol{f}_0(n+1)=\tilde{\boldsymbol{f}}_0^{\prime}(n) \end{align}\] \[\begin{align} \boldsymbol{f}_1(n+1)=\boldsymbol{0}_{1\times m(K+1)} \end{align}\]

if \(0<u_n<1\): \[\begin{align} \boldsymbol{f}_0(n+1)=\tilde{\boldsymbol{f}}_1^{\prime}(n) \boldsymbol{Q}_{10}\exp\left(\boldsymbol{Q}_{00}(1-u_n) t_o\right) \end{align}\] \[\begin{align} \boldsymbol{f}_1(n+1)=\tilde{\boldsymbol{f}}_0^{\prime}(n) \boldsymbol{Q}_{01} \exp\left(\boldsymbol{Q}_{11} u_n t_o\right) \end{align}\]

if \(u_n=1\): \[\begin{align} \boldsymbol{f}_0(n+1)=\boldsymbol{0}_{1\times m(K+1)} \end{align}\] \[\begin{align} \label{eq:f1} \boldsymbol{f}_1(n+1)=\tilde{\boldsymbol{f}}_1^{\prime}(n) \end{align}\tag{13}\] Then, the convolution of the \(n\)th forward likelihood \(\boldsymbol{f}(n)\) can be calculated according to the Eq. (12 )-(13 ).

Similarly, the computations of backward likelihoods in \(n\)th period can be expressed as follows. \[\begin{align} \label{eq:b} \boldsymbol{b}(n)=\exp\left(\boldsymbol{Q}t_o\right)\tilde{\boldsymbol{b}}(n) \end{align}\tag{14}\] \[\begin{align} \tilde{\boldsymbol{b}}_0^{\prime} (n)=\exp\left(\boldsymbol{Q}_{00}(1-u_n) t_o\right) \boldsymbol{b}_0(n+1) \end{align}\] \[\begin{align} \tilde{\boldsymbol{b}}_1^{\prime} (n)=\exp\left(\boldsymbol{Q}_{11}u_n t_o\right) \boldsymbol{b}_1(n+1) \end{align}\]

if \(0<u_n<1\): \[\begin{align} \tilde{\boldsymbol{b}}_0(n)=\exp\left(\boldsymbol{Q}_{00}(1-u_n) t_o\right)\boldsymbol{Q}_{01}\tilde{\boldsymbol{b}}_1^{\prime}(n) \end{align}\] \[\begin{align} \tilde{\boldsymbol{b}}_1(n)=\exp\left(\boldsymbol{Q}_{11}u_n t_o\right)\boldsymbol{Q}_{10}\tilde{\boldsymbol{b}}_0^{\prime}(n) \end{align}\]

if \(u_n=0\): \[\begin{align} \tilde{\boldsymbol{b}}_0(n)=\tilde{\boldsymbol{b}}_0^{\prime}(n) \end{align}\] \[\begin{align} \tilde{\boldsymbol{b}}_1(n)=\boldsymbol{0}_{m(K+1)\times 1} \end{align}\]

if \(u_n=1\): \[\begin{align} \tilde{\boldsymbol{b}}_0(n)=\boldsymbol{0}_{m(K+1)\times 1} \end{align}\] \[\begin{align} \label{eq:hatb} \tilde{\boldsymbol{b}}_1(n)=\tilde{\boldsymbol{b}}_1^{\prime}(n) \end{align}\tag{15}\] Then, \(\boldsymbol{b}(1)\) can be derived according to the convolution equations above of Eq.(14 )-(15 ).

In general, the computations of \(\boldsymbol{f}(n)\) and \(\boldsymbol{b}(1)\) take too much cost with the number of phases \(m\) and the capacity \(K\) increasing. To reduce the computation, \(\boldsymbol{f}(n)\) and \(\boldsymbol{b}(1)\) can be computed by using uniformization technique [46]. Let \(q\) be a positive constant, and \(q\) is larger than the maximum value of the absolute diagonal element of \(\boldsymbol{D}_0\). Then, the equation for the uniformization is given by \[\begin{align} \exp \left(\boldsymbol{D}_0t\right)=\sum_{z=0}^{\infty} e^{-qt}\frac{(qt)^z}{z!}(\boldsymbol{I}+\boldsymbol{D}_0/q) \end{align}\] where \(\boldsymbol{I}\) is the \(m\times m\) identity matrix. Since there exists the sum operation of the above equation, whose upper limit is \(\infty\), we usually truncate it by a certain point in the practical computation. In general, a simple and effective way to truncate the upper limit is based on Poisson probability mass function.

Now, we summarize the EM algorithm to estimate the parameters of the MAP with utilization data.

In step 5, we need to determine the converge condition to out the iteration process. A simple method of the condition to stop the iteration is based on the difference between the two successive likelihoods or the relative difference between the two successive likelihoods. The coverage condition is formulated as \[\begin{align} |\text{LLF}(\boldsymbol{\theta}^\prime)-\text{LLF}(\boldsymbol{\theta})| < \epsilon_a,~~\frac{|\text{LLF}(\boldsymbol{\theta}^\prime)-\text{LLF}(\boldsymbol{\theta})|}{|\text{LLF}(\boldsymbol{\theta})|}< \epsilon_r \end{align}\] where \(\boldsymbol{\theta}^\prime\) and \(\boldsymbol{\theta}\) represent the parameters at the current iteration and the previous iteration. Notice that in this paper, we use \(P({\mathcal{O}})=\boldsymbol{\pi} \boldsymbol{b}(1)\) as the LLF.

4.6 Optimal Model Selection↩︎

Consider a \(m\)-state MAP. The value \(m\) determines the number of phases of a MAP. In general, the accuracy of the model improves with the number of phases \(m\) increases. However, a large \(m\) sometimes may cause the overfitting problem, so that the accuracy of fitting gets worse. In our previous work [45], we used AIC [47] to select the optimal model for \(NHPP/M/1/K\) queueing systems. And the experiments were verified that AIC can work well in selecting optimal queueing models. In this paper, we also use the AIC, which can quantify the goodness of fit of the model. The formula below definites the AIC: \[\begin{align} \text{AIC}=-2\text{LLF}+2\left(\#~\text{of free parameters}\right) \end{align}\]

In the case of general MAP with \(m\) phases, the number of free parameters in \(\boldsymbol{\pi}, \boldsymbol{D}_0\) and \(\boldsymbol{D}_1\) is \(2m^2-1\). On the other hand, there are \(m^2\) free parameters in an MMPP. Then, the optimal number of phases can be computed by determining the smallest value of the AIC.

5 Conclusion↩︎

This paper investigated the parameter estimation problem for MAP-driven QBD queueing systems using utilization data. Unlike conventional approaches that rely on inter-arrival times, waiting times, or queue-length observations, utilization data only records the proportion of busy time within each observable monitoring interval. As a result, key information such as arrival epochs, service completions, phase transitions, and queue-length trajectories is partially or completely unobserved. To address this incomplete-data setting, we formulated the queueing system as a QBD process and developed an EM algorithm for maximum likelihood estimation of the MAP and service parameters.

The proposed method explicitly separates observable and unobservable intervals and derives the expected sufficient statistics required in the E-step, including cumulative sojourn times, phase transitions, arrivals, and service completions. These quantities are then used in the M-step to update the model parameters iteratively. In addition, AIC-based model selection is incorporated to determine the number of MAP phases, providing a practical way to balance fitting accuracy and model complexity.

The main significance of this study is that it enables MAP parameter estimation from utilization data, which is much easier to obtain in practical computer systems than detailed event-level observations. By combining MAP, QBD modeling, and EM-based inference, the proposed framework extends queueing parameter estimation to a more realistic monitoring scenario. Future work will include more extensive validation with real system traces, robustness analysis under different monitoring resolutions, and extensions to more general service-time distributions and multi-server queueing systems.

References↩︎