ISSN : 2234-6473 (Online)

**Industrial Engineering & Management Systems Vol.12 No.3 pp.281-287**

DOI : https://doi.org/10.7232/iems.2013.12.3.281

# A Roots Method in GI/PH/1 Queueing Model and Its Application

### Abstract

### 1. INTRODUCTION

Efficient numerical solutions of queues have attracted attention of many researchers in computational probability over the past few years. As a result, several methods of inversion of transforms have been developed. Since transform methods have been considered numerically difficult to implement, they started using the matrixgeometric method (MGM) for analyzing various queueing models such as GI/M/1 and other related models. Therefore, many researchers have focused on getting the exact square matrix R that is nonnegative solution to matrix- quadratic equation as fast as possible.

This paper applies the roots method which was first introduced by Chaudhry et al. (2012) in dealing with the model GI/C-MSP/1 that a transform method can be exploited to get the unknown constants and the steady-state probability. In this connection, see also Chaudhry et al. (1990), Chaudhry and Templeton (1983) and references therein. The roots method first accurately calculates the unknown constants and then uses these constants to get probability distributions of the number in the system at pre-arrival, arbitrary, and post-departure epochs as well as the sojourn-time distribution and other performance measures.

We also investigate the quality of products before making a decision to buy some important goods such as electronics and oil. The airplane oil inspection is conducted to prevent an accident like a drop or a breakdown. Since airplane oil inspection has also several steps, sometimes airplane oil inspection has caused oil train congestion. The growing oil train congestion can lead to an increase in waiting time and running short of fuel. Therefore, it is important to analyze the phenomenon with more accuracy because it can create a deep bottleneck for the rapid and safe logistics support.

The system under consideration has a special feature in which interarrival time of oil occurs by a reorder. It is important to maintain a constant inventory level of oil for continuous flights in the Republic of Korea Air Force (ROKAF). To maintain a constant level of oil, it is necessary to reorder oil at a right time. A reorder is mainly affected by the flight sortie with scheduling. The interarrival time from the actual data follows a normal distribution with a small deviation in this system. We also presented a higher description model for the airplane oil inspection procedure using phase-type rather than exponential distribution. Since the inspection is conducted by skilled airmen whose maintaining cost is quite high, an air base has just one inspection team. For this reason, we used the GI/PH/1 queueing model. The obtained results offer fundamental data which support safe and economic operations. In addition, the model may be applicable to various other fields as well as oil train congestion problems.

We introduce a method to analyze oil train congestion problem using the roots method in GI/PH/1 queueing system. We also present an oil train congestion analysis model with complicated inspection steps. This paper gives information about the mean number of oil trains and their sojourn time in one of ROKAF’s bases.

This paper is organized as follows. In Section 2, we briefly review the related literature. In Section 3, we describe the procedures of airplane oil inspection in ROKAF. The numerical analysis uses the roots method to get the steady-state probability at three epochs from the associated characteristic equation of the vector-generating function of system-length distribution. In Section 4, we deal with the transition rate matrix and numerical results on the probability of the system-length distribution at several epochs as well as performance measures.

### 2. LITERATURE REVIEW

Many researchers have theoretically analyzed several GI/PH/1 queueing models with various types of services as well as arrivals. However, very little seems to have been done on the application of GI/PH/1 queueing models to the real world. In addition, it was very hard to find researches in GI/PH/1 queueing model using the roots method.

Ramaswami and Latouche (1989) and Neuts (1981) discussed a GI/PH/1 queueing model using the MGM. Ramaswami and Lucantoni (1988) derived recursive relations for computing the higher moments of the stationary waiting time distribution in a stable GI/PH/1 queue. Neuts (1981) indicated stationary probability distributions of the waiting time of an arriving customer and of the virtual waiting time, as well as that of the time-in-system which may be computed by the numerical integration of highly structured, finite systems of differential equations with constant coefficients. Grassmann (1982) complements the work of Neuts (1981) by giving an efficient method based on randomization for finding the transition matrix of the imbedded Markov chain. Kao (1991) also computed steadystate probabilities of GI/PH/1 queueing model using the state reduction procedure. Baum and Breuer (2006) suggested applying Foster’s criteria that are used to obtain necessary and sufficient conditions ensuring the recurrence and ergodicity of an embedded Markov chain to GI/PH/1 queueing system. Sengupta (1990) and Asmussen (1992) showed that the waiting time has a matrix exponential representation, presenting algorithms for the GI/PH/1 queue using a nonlinear matrix equation which can be solved using an iterative technique. Asmussen et al. (1996) showed that the waiting-time distribution in a GI/PH/1 queue is also phase-type. Kao (1996) compared alternative approaches to computing the R matrices of GI/PH/1 queues. The alternative approaches are composed of state reduction algorithm described in Kao (1991), two algorithms proposed by Latouche (1993) and an algorithm using the notion of matrix exponential proposed by Sengupta (1990). Kim (2006) analyzed an asymptotic behavior of the loss probability for the GI/PH/1/K queue as K tends to infinity when the traffic intensity ρ is strictly less than one. Dudin et al. (2004) studied GI/PH/1 queueing system with a BMAP flaw of negative customers. Chakravarthy (1992) showed efficient algorithmic procedures for the steady-state analysis of a GI/PH/1 queueing model with group services. Dudin and Klimenok (2003) considered an optimal admission control of customers into a GI/PH/1 queueing system. Guo et al. (2007) discussed the queue system GI/PH/1 with repairable service, where the lifetime and the repair time of the service station are both PH random variables. Sengupta (1989) extended the results for the discrete problem to the continuous case, providing the waiting time and queue length distributions of the GI/PH/1 queue with firstin firstout discipline. In addition, he showed that the method is also applicable in solving an inventory problem in which the supply and demand are subject to seasonal fluctuations.

### 3. MODEL DESCRIPTION

#### 3.1 GI/PH/1

We considered the GI/PH/1 system with first-come first-served discipline. The system features arrivals generated by a renewal process with the interarrival time distributed by a distribution function F(x) and the service times distributed with PH(τ , T) distribution of phase m∈N. PH(τ , T) defined by the interval that Markov chain (MC) having the following infinitesimal generator enters into absorbing state where the initial state is determined by the m-state vector τ . The infinitesimal generator Q for the service process is

where a is the vector of size m. Its entries define transition rates of MC from transient states to absorbing state. The m×m dimensional matrix T is always nonsingular and invertible. Further, it is clear that Te + a = 0, where e = (1, …, 1)^{T}, since each row of Q sums to zero. For more details on phase-type distribution, see Latouche and Ramaswami (1999) and Neuts (1981).

With this description, the inspection process can be presented as a two-dimensional Markov process, {(N(t), J(t)), t ≥ 0}, where N(t) and J(t) represent the number of customers (oil trains) in the system and the service step at time t, respectively. If the initial vector is τ , the generator matrix of the service phase process S is given by the matrix

Based on the above explanation, we denote the transition probability matrix of the number of customers in the system atarrival instants by

The general structure of is determined by two considerations. First, between two arrival instants, the number of users in the system can increase by at most one. Second, as long as the system is not empty, the change of the number of users between arrival instants T_{n}, T_{n+1 }does not depend on the number of users in the system at time T_{n}. The entry (A_{k})_{i,j} of the matrix A_{k} indicates the probability that between two consecutive arrivals instants T_{n} and T_{n+1}, k customers are served and the phase of service changes from i to j. The entry (B_{k})_{i,j} of the matrix Bk indicates the probability that between two consecutive arrivals instants T_{n} and T_{n+1}, at least k+1 customers are served and the phase of service changes from i to j. Then, A_{n}, B_{n} can be given by the following equations.

where F(t) is the distribution function of the interarrival time of a customer and the (i, j)th entry of the matrix P(n,t) represents the probability that n customers are served and the phase of the service process changes from i to j during time interval (0, t].

#### 3.2 Analysis of the Model

For the analysis of the model described in the previous section, we use the result from Chaudhry et al. (2012). We summarize the procedures in brief. Let π^{-} = (π^{-}(0), π^{-}(1), π^{-}(2), …, π^{-}(n)), n ≥ 0, be the stationary distribution of number of customers in the system at a pre-arrival epochs, and let π^{-*}(z) = be the vector generating function (VGF) of the distribution. Using , we have

Multiplying (6) and (7) by z^{0} and z^{n}, respectively, umming them, and using VGF , we finally get

where adj[M] is the adjoint of a square matrix M, det[M] is the determinant of M, and I_{m} is the identity matrix of order m×m. To evaluate the vector π^{-}(n), n≥0 in (8), we first determine an analytic expression for each component of π_{j}^{-*}(z). As each component π_{j}^{-*}(z) is convergent in |z| ≤ 1, the zeros of det[I_{m}-zA(z^{-1})] whose absolute value is less than or equal to one must be the zeros of the numerator of each component (8). So, the equation of det[I_{m}z-A(z)] = 0 has exactly m roots (γ_{i}) inside the unit circle |z| = 1. We assume that < 1 for the stability of the model and all γ_{i} (1 ≤ i ≤ m) are distinct. Applying partialfraction method used by Chaudhry and Templeton (1983, p. 229), we have

Where k_{ij} are constants to be determined. Expanding both sides of (9) in series at the point z = 0, we get the following result for system-length distribution at pre-arrival epoch.

and hence

Using results of Markov renewal theory and semi- Markov processes, we get system-length distribution at arbitrary epoch. Let Ω_{n} and Φ_{n} denote the matrices of order m×m whose (i, j)th element represents the limiting probability that n customers are served during an elapsed interarrival time of the arrival process and the service process is in phase j at the end of n-th service completion, given that there were exactly n and at least (n+1) customers in the system respectively, and the state of the service process was i. This leads to

The mean system length can be obtained as Eq. (14).

From Little’s law, we also have mean sojourn time was Eq. (15).

Also, using level-crossing arguments given in Chaudhry and Templeton (1983, p. 299), we get system-length distribution at post-departures as follows:

#### 3.3 The Procedure of Airplane Oil Inspection

In this paper, we describe the procedure of airplane oil inspection receiving oil from trains in one of the ROKAF’s bases. Since poor quality oil could cause accidents, the airplane oil inspection is one of the most important works in the ROKAF. The oil is supplied usually 3?5 times a week depending on attrition volume, and most of the inspections are done manually. Oil tank capacity is limited because to build new oil tanks is costly. Due to the shortage of fuel tanks, usage (consumption) increases lead to the increase in the number of reorders. Increase of the frequency of reorder could cause oil train congestion in inspection process and further delay in lead time. Therefore, we have to decide on the reorder point considering further delay caused by oil train congestion.

The procedure of airplane oil inspection is described in Figure 1. Once oil arrives, there is an examination such as sealed, cover state check and water examination with the naked eye. If the first step is good, it has to proceed to the second step called “lower layer sampling test” composed of precipitation reaction and freezing prevention. Third step is an “all-layer test” consisting of six inspection items such as conductivity and filtering test. Finally, they get a receiving test. We assumed that the inspection went into the first step without any queueing if other steps fail. Phase-type distributions cannot only be extended with little extra complication, but also explain the complex phenomenon correctly such as inspection procedures.

**Figure 1.** The procedure of airplane oil inspection.

### 4. CASE STUDY AND NUMERICAL RESULTS

#### 4.1 Interarrival Time

In this study, the intervals of trains transporting oil are distributed by a normal distribution with a small deviation. A small deviation is caused by oil’s lead time. Table 1 gives statistics for interarrival time. It is necessary to evaluate more suitable distribution when squared error is lower and corresponding p-value is higher than 0.05.

**Table 1.** Results of the interarrival time analysis

However, this distribution is difficult to apply for the queueing model. Therefore, we use a phase-type distribution to fit the observed interarrival time and compute the mean waiting time in the queue with this fitted interarrival time distributions. To fit a phase-type distribution with the observed interarrival time, we use the expectationmaximization (EM) algorithm for phase-type distributions. The EM algorithm is an iterative method for maximum likelihood estimation. For more details on the EM algorithm, see Asmussen et al. (1996) and Dempster et al. (1977).

There is a slight difference between real data being approximated by a normal distribution and phase-type distribution as may be seen in Figure 2.

**Figure 2.** The comparison of interarrival time between normal and phase-type distribution.

We could get the estimated phase-type distributions 3×3 matrix for interarrival time. The more the number of phase is, the smaller χ^{2} value is. Since χ^{2} is almost the same after the 3rd phase, we approximated interarrival time by 3×3 phase-type distributions matrix for normal interarrival time.

#### 4.2 Service Time

Table 2 gives the statistical results on the assumption of the service time whether it follows phase-type or exponential distribution. If service time is assumed to have exponential distribution, null hypothesis is rejected because p-value is below 0.05. Therefore, it cannot be assumed that service time is an exponential distribution. Meanwhile, in case service time is assumed to have phasetype distribution, we can adopt the null hypothesis because p-value is over 0.05. We conclude that it is better to assume that the service time follows phase-type rather than exponential distribution.

**Table 2.** Results of the service time analysis

There are many differences between real data being approximated by exponential distributions and phase-type distributions as may be seen in Figure 3. We can see the phase-type distribution fits the data better than exponential distribution as in Figure 3.

**Figure 3.** The comparison of fitting result for service time between phase-type and exponential distributions.

On analyzing service time, we assumed that every inspection step has the exponential distribution with rates μ_{1}, … , μ_{4}. The s_{1}, … , s_{4} and f_{1}, … , f_{4} are test success and failure probabilities of each step, respectively. We present details of the transition rate probability, PH(τ, T) of service time in Table 3.

**Table 3.** The transition rate probability of service time

We analyzed system-length distribution at several epochs on the basis of the result equation in Section 3.2. Some outputs of this model are shown in Tables 4–6, where the results are given for the GI/PH/1 queue. In addition, Table 7 lists remarkable differences between the models using exponential service time and phase-type service time.

**Table 4.** System-length distribution at pre-arrival epoch

**Table 5.** System-length distribution at arbitrary epoch

**Table 6.** System-length distribution at post-departure epoch

**Table 7.** The comparison of performance measures with GI/PH/1 and GI/M/1

In case a manager makes a decision for order intervals based on a model with exponential distribution, he or she makes a mistake to increase the inventory level based on the false information. Through this study, we found that the application of wrong distributions can cause different customer’s interarrival time.

We calculated the performance measures of a GI/M/1 system using the roots method suggested by Chaudhry and Templeton (1983, p. 116) for the model GI^{r}/M/1 queue, where r is the exact batch size.

### 5. CONCLUSION AND FUTURE SCOPE

In this paper, we applied GI/PH/1 queueing model to the oil train congestion problem and analyzed severalepochs probabilities of our model by applying the roots method suggested by Chaudhry for the model GI/CMSP/ 1. Since the procedure of airplane oil inspection is composed of several service phases, it has caused oil train congestion sometimes. We assumed oil inspection time as phase-type distribution which can present more detailed information. We used a simple procedure to evaluate the steady-state system-length distribution at several epochs.

Our analysis of the interarrival time and inspection (service) time of oil supply from the actual data which is given from one of the ROKAF’s bases indicate that the interarrival time of oil follows a normal distribution with a small deviation, and the service time follows phase-type distribution, which was first introduced by Neuts to deal with the shortfalls of exponential distributions.

### Reference

2.Asmussen, S., Nerman, O., and Olsson, M. (1996), Fitting phase-type distributions via the EM algorithm, Scandinavian Journal of Statistics, 23(4), 419-441.

3.Baum, D. and Breuer, L. (2006), Applying Foster's criteria to a GI/PH/1 queueing system, Cybernetics and Systems Analysis, 42(3), 433-439.

4.Chakravarthy, S. (1992), A finite capacity GI/PH/1 queue with group services, Naval Research Logistics, 39(3), 345-357.

5.Chaudhry, M. L. and Templeton, J. G. C. (1983), A First Course in Bulk Queues, John Wiley & Sons, New York, NY.

6.Chaudhry, M. L., Harris, C. M., and Marchal, W. G. (1990), Robustness of root finding in single-server queueing models, INFORMS Journal on Computing, 2(3), 273-286.

7.Chaudhry, M. L., Samanta, S. K., and Pacheco, A. (2012), Analytically explicit results for the GI/C-MSP/1/∞ queueing system using roots, Probability in the Engineering and Informational Sciences, 26(2), 221- 244.

8.Dempster, A. P., Laird, N. M., and Rubin, D. B. (1977), Maximum likelihood from incomplete data via the EM algorithm, Journal of the Royal Statistical Society: Series B, 39(1), 1-38.

9.Dudin, A. N. and Klimenok, V. I. (2003), Optimal admission control in a queueing system with heterogeneous traffic, Operations Research Letters, 31(2), 108- 118.

10.Dudin, A. N., Kim, C. S., and Semyonova, O. V. (2004), An optimal multithreshold control for the input flow of the GI/PH/1 queueing system with a BMAP flow of negative customers, Avtomatika i Telemekhanika, 9, 71-84.

11.Grassmann, W. K. (1982), The GI/PH/1 queue: a method to find the transition matrix, INFOR, 20, 144-156.

12.Guo, M. M., Tian, N. S., and Liu, A. Y. (2007), Queue system GI/PH/1 with repairable service station, Operations Research and Management Science, 16(5), 69-74.

13.Kao, E. P. C. (1991), Using state reduction for computing steady state probabilities of queues of GI/PH/1 types, ORSA Journal on Computing, 3(3), 231-240.

14.Kao, E. P. C. (1996), A comparison of alternative approaches for numerical solutions of GI/PH/1 queues, INFORMS Journal on Computing, 8(1), 74-85.

15.Kim, J. S. (2006), Asymptotic analysis of the loss probability in the GI/PH/1/K queue, Journal of Applied Mathematics & Computing, 22(1), 273-283.

16.Latouche, G. (1993), Algorithms for infinite Markov chains with repeating columns. In Meyer, C. D. and Plemmons, R. D. (eds.), Linear Algebra, Markov Chains and Queueing Models, Springer, New York, NY, 231-265.

17.Latouche, G. and Ramaswami, V. (1999), Introduction to Matrix Analytic Methods in Stochastic Modeling, Society for Industrial and Applied Mathematics, Philadelphia, PA.

18.Neuts, M. F. (1981a), Matrix-Geometric Solutions in Stochastic Models: An Algorithmic Approach, Johns Hopkins University Press, Baltimore, MD.

19.Neuts, M. F. (1981b), Stationary waiting-time distributions in the GI/PH/1 queue, Journal of Applications Probability, 18(4), 901-912.

20.Ramaswami, V. and Latouche, G. (1989), An experimental evaluation of the matrix-geometric method for the GI/PH/1 queue, Communications in Statistics: Stochastic Models, 5(4), 629-667.

21.Ramaswami, V. and Lucantoni, D. M. (1988), Moments of the stationary waiting time in the GI/PH/1 queue, Journal of Applications Probability, 25(3), 636-641.

22.Sengupta, B. (1989), Markov processes whose steady state distribution is matrix-exponential with an application to the GI/PH/1 queue, Advances in Applied Probability, 21(1), 159-180.

23.Sengupta, B. (1990), The semi-Markovian queue: theory and applications, Communications in Statistics: Stochastic Models, 6(3), 383-413.