Queue with Batch Service, Batch Size Determined by Emergency Customers Sinu Lal T. S.1 , A. Krishnamoorthy1 , and V. C. Joshua1 Department of Mathematics, CMS College Kottayam 686001, Kerala, India {sinulal, krishnamoorthy, vcjoshua}@cmscollege.ac.in http://www.cmscollege.ac.in Abstract. We consider a queueing system offering batch service for het- erogeneous customers. Two types of customers called ordinary customers and emergency customers arrive to the system according to a marked Markovian arrival process of order 2. There is an infinite queue and finite buffer in front of the station. Ordinary customers arrive to the infinite queue and emergency customers to the finite buffer. Service is provided in batches with maximum batch size K. If the number of emergency cus- tomers is greater than or equal to K, the first K of them will be taken together into service and if it is less than K, ordinary customers are taken along with emergency customers so as to maximize the batch size. In the absence of emergency customers, if there are at least K ordinary customers, then first K of them will be served next. For a batch of size i, 1 ≤ i ≤ K, the service time follows a phase type distribution with representation P H(αi , Si ). The model possesses the characteristics of a vacation queueing model. The system is analyzed using Matrix analytic Method. Performance characteristics are derived. The model is numeri- cally illustrated with suitable example. Keywords: Batch Service, Marked Markovian Arrival Process, Emer- gency Customer, Matrix Analytic Methods, Phase Type Distributions 1 Introduction In this paper we study the mathematical model of a courier delivery system with emergency arrivals and bulk service. The model has found applications in priority based signal transmission systems as well. Numerous real life problems can be modelled as queues with bulk service and a rich literature on bulk queues is available. Delivery systems of online marketing, courier services, transportation vehicles and airline services etc. are typical examples for bulk service systems. Each of these systems possesses its own special features. In delivering items, the emergency orders are usually served with priority. Sometimes instantaneous service is called for the emergency customers due to urgency (e.g. delivery of highly perishable items like samples of body fluids for diagnostic service, animal sperm, medical essentials upon order placing etc.). It is effective to incorporate threshold policies in bulk queues to maximize revenue generation and a great exploration for mathematical results is possible in this direction. Copyright © 2020 for this paper by its authors. Use permitted under Creative Commons License Attribution 4.0 International (CC BY 4.0). One of the early works in bulk queues is due to M. F. Neuts [13]. In this work, the customers are taken together to the service station until queue length reaches a specified level (L). If the queue grow beyond L and reaches K all the customers up to K are taken together, the customers beyond K must wait in queue. H. Gold and P. Tran-Gia [6] studied an M/G[a, b]/1−S queueing system. The main motivation for this model is the manufacturing systems with batch serving stations (machines for computer components and chip production). M. L. Chaudhry and U. C. Gupta [5] describe an M/Ga,b /1/N queue and queue is studied using supplementary variables and embedded Markov chain techniques. A finite buffer M/G/1 queue with general bulk service rule and single vacation is studied by U. C. Gupta and K. Sidkar [8] and in this model batch size is restricted to a range of values, and when the queue length falls below the infimum, the server goes for a vacation. W. B. Powell [16] studies a general class of vehicle dispatching strategies for bulk arrivals, bulk service queue. D. J. Van De Rzee, A. Van Haeten and P. C. Schuur [22] describe control strategies for reading the cost of multi server batch operation systems. S. Kuppa and G. Dattatreya [10] describe frame ag- gregation in unsaturated WLAN’s with finite buffers and in this model method for aggregation of more than one data form for transmission is described as an application for bulk queues. A. Banerjee and U. C. Gupta [1] analyze a single server finite buffer queue with Poisson arrival and bulk service and batches are arbitrarily distributed and depends the size of the current batch in service and the aim is to reduce congestion in the system. J. Baetens, B. Steyaert, D. Claeys, and H. Bruneel [2] analyze a discrete time BM AP/Gl,c /1 queue in which service time of a batch depends the current batch size and a timing mechanism to reduce waiting time is also proposed. M. Yu and A. S. Alfa [23] describe an algorithm for computing the queue length distribution at various time epochs in a DM AP/Ga,b /1/N queue with batch size dependent service times and the impact of correlation on system characteristic is also discussed. [11,7,4] and [2] give live descriptions of batch service queues with dependent batch size. A. Sikadar and S. K. Samantha [19] studied a bulk service queue with service vacation and vacation starts when all the customers are exhaustively served. In [17] G. V. Krishna Reddy, R. Nadarajan and P. R. Kandasamy considered a bulk service system with heterogeneous arrivals and bulk service is provided to the class of customers with low priority. In [2] J. Baetens, B. Steyaert, D. Claeys, and H. Bruneel analyze a two class batch service queueing model with variable server capacity and batch size is determined by number of consecutive same-class customers. In [3] J. Baetens, B. Steyaert, D. Claeys, and H. Bruneel analyze delay of a random customer in a two class batch service queueing model with variable service capacity and batch size is determined by the length of the sequence of same class customers. An excellent survey by S. Sasikala, K. Indhira on bulk service queueing models is demonstrated in [18]. A phase type distribution may be defined as the distribution of the time until absorption takes place in a Markov process with a finite state space and a single absorption state defined over nonnegative real line. A phase type distribution with transient states {1, 2, . . . , n} and an absorbing state n + 1 is represented by a two tuple of the form (α, T ), where α is the probability vector of length n according to which the process  the initial state from {1, 2, . . . , n} and T  selects T T0 is an n × n matrix such that generates the process, given the column 0 0 vector T 0 satisfies the condition T e+T 0 = 0, where e is the vector of ones. (α, T ) is called the representation of the phase type distribution. The distribution F of time until the chain gets absorbed into the state n + 1 is given by F (x) = 1 − αexp(T x)e, x ≥ 0. The set of all phase type distributions is a dense subset of the set of all distributions on the non-negative real line and hence it is a best tool to approximate any arbitrary distribution in this set. For more descriptions on phase type distributions see [14]. A Marked Markovian arrival process (MMAP) is a stochastic point process with heterogeneous arrivals in discrete or continuous time. MMAP may be de- scribed as follows: Let C be set of indices which describes different types of customers. Let Nh (t) be the number of arrivals of type h in [0, t] such that Nh (0) = 0, ∀h ∈ C. Consider the set of nonnegative matrices {Dh : h ∈ C}. Let D0 is a matrix with nonnegative P off-diagonal elements and negative diag- onal elements, and D = D0 + h∈C Dh be an infinitesimal generator of order m, and {I(t) : t ≥ 0} be a continuous time Markov chain defined by D. Then (D0 , Dh , h ∈ C) defines an MMAP {Nh (t), h ∈ C, I(t), t ≥ 0}. The Markov chain {I(t) : t ≥ 0} is called the underlying Markov chain of {Nh (t), h ∈ C, I(t), t ≥ 0}. For description of MMAP model see [9]. Analysis of queues using matrix analytic method can be seen in [20,21] Upcoming sections are arranged as follows: Section 2 describes the mathe- matical model stability condition and stationary distribution is also obtained in this section. Performance characteristics are included in Section 3. Service time analysis is done in Section 4. Waiting time analysis is given in Section 5. The model is numerically illustrated in Section 6. Section 7 concludes the work. 2 Mathematical Model The customers arriving to the system are of two types, namely ordinary cus- tomers and emergency customers. The arrival is according to a MMAP deter- mined by the matrices D0 , D1 and D2 . D0 gives the rate of transitions in the underlying process without an arrival, D1 and D2 gives the rate of transitions for ordinary and emergency arrivals respectively. The matrix D = D0 + D1 + D2 is an infinitesimal generator of the underlying process. If θ is the steady state distribution of D then λi = θDi e, i = 1, 2 is the fundamental rates of ordinary and emergency arrivals. If the covariance of number of ordinary customers n1 (t) and that of emer- gency customers n2 (t) is given by 2 X cov(n1 (t), n2 (t)) = −(2λ1 λ2 + θ( Dk (D − eθ)−1 D3−k )e)t k=1 2 X + θ( Dk (D − eθ)−1 exp(DtI (D − eθ)−1 )D3−k ))e. (1) k=1 The service discipline is as follows: service is provided in batches with maximum batch size K. If the number of emergency customers is greater than or equal to K, the first K of them will be taken together into service and if it is less than K, ordinary customers are taken along with emergency customers so as to maximize the batch size. In the absence of emergency customers, if there are at least K ordinary customers, then first K of them will be served next. Service time for a batch of size j is distributed with a phase type representation P H(αj , Sj ). The following are system descriptors at the time t. • N1 (t) - Length of the infinite queue. • N2 (t) - Number of customers in finite buffer. • I(t) - Status of the server. • S(t) - Phase of the service. • a(t) - Phase of MMAP. The server status is defined as follows.  0, if the server is idle I(t) = i, if the sever is busy with i customers, 1 ≤ i ≤ K. We define Y(t) = {N1 (t), N2 (t), I(t), S(t), a(t)}. Then {Y(t), t ≥ 0} is a continuous time Markov process on the state space S = ∪i≥0 L(i). For each i < K, L(i) = L1 (i) ∪ L2 (i), where L1 (i) = {(i, 0, l), 1 ≤ l ≤ a} and L2 (i) = {(i, i1 , 1, j, k), 0 ≤ i1 ≤ M, 1 ≤ k ≤ a, 1 ≤ j ≤ r}. L1 contains the states in which there are no priority customers and the server is idle and L2 (i) corresponds to states in which server is active and in this case buffer can be empty or non empty. For i ≥ K, L(i) = {(i, i0 , 1, j, k), 0 ≤ i0 ≤ M, 1 ≤ k ≤ a, 1 ≤ j ≤ r}. The states in all the levels greater than or equal to K will only contain states in which server is busy. The infinitesimal generator Q of Y takes the form   A00 A01  A10 A11 A12     A20 A21 A22     . .. ..  ..   . .    AK−10 A K−11 A K−12    Q= # ,  A1 AK1 A0  # A2 AK+11 A0     .. .. ..   . . .      #   A K A 2K−11 A 0  ..   .. .. . . .  A00 = D0 , A01 = D1 α1 ⊗ D2 , j X Sj0 ⊗ Ia , 1 ≤ j ≤ K,  Aj0 = er ⊗ i=0 i X  Ai1 = diag D0 , ( Si ) ⊕ D0 , 1 ≤ i ≤ K − 1, j=1  Ai1 = diag D0 , SK ⊕ D0 , i > K,   D1 αi ⊗ D2 Ai2 = , 1 ≤ i ≤ K − 1, Ir ⊗ (D1 + D2 )   αK ⊗ (D1 + D2 )  AK2 = , Ai2 = Ir ⊗ (D1 + D2 ) , i > K, Ir ⊗ (D1 + D2 ) A# 0 0  l = er ⊗ (SK ⊗ Ia ) er ⊗ (αl ⊗ (SK ⊗ Ia )) , K + 1 ≤ l ≤ 2K, A# 0 0  l = er ⊗ (SK ⊗ Ia ) er ⊗ (αl ⊗ (SK ⊗ Ia )) , l > 2K. Note that the Kronecker sum ⊕ and product ⊗ are used. The non zero blocks in the above matrix arise due to the transitions described below. From To Transition rate (i,0,0,l), 0 ≤ i < K − 1 (i,i’,0,k) D0 (l, k) (i,0,0,l), i ≥ K − 1 (i,1,j,k) D1 (l, k) (h,0,0,l), h ≥ 0 (h,1,j,k) αhj D2 (l, k) (i,0,g,j,k), 1 ≤ i ≤ K − 1 (i,0,g,0,k) Sg0 (j) (i,i’,g,j,k), 1 ≤ k, k ∗ ≤ a (i, i0 , g, j, k ∗ ) U (k, k ∗ ), U = (Sg ) ⊕ D0 (i,i’,g,j,k), ≤ k, k ∗ ≤ a (i, 1, j, k ∗ ) U (k, k ∗ ), U = Si ⊕ D0 , 1 j0 (i,1,j,k), K + 1 ≤ i ≤ 2K − 1 (0, 0, k) SK j0 (i,1,j,k), K + 1 ≤ i ≤ 2K − 1 (0,0,k) Sm ∗ ∗ 0 (i,1,j,k),K + 1 ≤ i ≤ 2K − 1, (i, 1, j, k ) V (k, k ), V = αi−K ⊗ (SK ⊗ Ia ) ∗ 1 ≤ k, k ≤ a (i,1,j,k), 0 ≤ i ≤ K − 2 (i, 1, j, k ∗ ) Ω(k, k ∗ ), Ω = Ir ⊗ (D1 + D2 ) (i,1,j,k),0 ≤ i ≤ K − 2 (i, 0, k ∗ ) Ω 0 (k, k ∗ ), Ω 0 = Ir ⊗ (D0 ) (i,1,j,k),i > K − 2 (i, 1, j, k ) Ω(k, k ∗ ), Ω = Ir ⊗ (D1 + D2 ) ∗ Clearly in the infinitesimal generator Q there are non-zero blocks A∗j , 1 ≤ j ≤ K, due to which Q loses its QBD structure. In order to gain a QBD structure for the matrix Q, the levels are redefined by merging them appropriately. The levels lK + 1 to (l + 1)K, l = 0, 1, 2, . . . , are merged together, and modified generator Q0 is described below. After merging of cells the modified form of Q take the structure   A00 A01  B10 B11 B12    0 Q =  B 20 B 21 B 0 .    B 2 B 1 B 0   .. .. .. . . . A00 = D0 , A01 = (D1 α1 ⊗ D2 ), B10 = [A∗1 A∗2 . . . A∗K ]T , B11 = diag(A11 , A21 , . . . , AK1 ) + (Ω0 + diag(A12 , A22 , . . . , Ak−12 )), B21 = diag(AK+11 , AK+21 , . . . , A2K1 ) + (Ω0 + diag(A12 , A22 , . . . , Ak−12 )), B1 = diag(A2K+11 , A2K+21 , . . . , A3K1 ) + (Ω0 + diag(A12 , A22 , . . . , Ak−12 )), B12 = diag(diag(D0 , Ir ⊗ (D1 + D2 ), . . . , diag(D0 , Ir ⊗ (D1 + D2 )), B2 = diag(A# # # 1 , A2 , . . . , AK ), B0 = diag(Ir ⊗ (D1 + D2 ), . . . , Ir ⊗ (D1 + D2 )). Clearly Q0 posses a QBD structure. The matrix B = B0 + B1 + B2 is an infinitesimal generator and the steady state distribution of B is obtained below. B = IK ⊗ diag(V1 , V2 , . . . , Var ), Vi = e ⊗ αi ⊗ (S30 ⊗ Ia + SK ⊗ Ia + Ir ⊗ D), 1 ≤ i ≤ ar. Let π = (π1 , π2 , . . . πK ) be the corresponding steady state probability vector, πi = (πi1 , πi2 , . . . , πiK ), i ≥ 1. π may be obtained as the solution to the system πB = 0 and πeKar = 1 where eKar is a column vector of ones having length Kar. B is a singular matrix as it can have at most Kar − 1 linearly independent rows. Hence πB = 0 has a non trivial solution π, this solution can be normalized to satisfy the condition πeKar = 0. Theorem 1. The system is stable if and only if K X K X πi Ir ⊗ (D1 + D2 ) < πi (e ⊗ αi ⊗ (S30 ⊗ Ia )). i=1 i=1 Proof. The system is stable if and only if πB2 e < πB0 e, see [15]. From matrices defined in Lemma 2 XK πB2 e = πi Ir ⊗ (D1 + D2 ), i=1 K X πB0 e = πi (e ⊗ αi ⊗ (S30 ⊗ Ia )). i=1 t u The stationary distribution of the system process is obtained as follows. Un- der the assumption of the stability condition, the steady state probability dis- tribution exists. Let y = (y0 , y1 , y2 , . . . ) be the steady state probability vector of the Markov Chain Y, where yi = (yi0 , yi1 , . . . , yiM ). Then y is the unique solution to the system of equations yQ = 0 and ye = 1. Then each component yi ia a vector of length Kar. From yQ0 = 0 and ye = 1, we have the system of equations y0 A00 + y1 B10 = 0 y0 A01 + y1 B11 + y2 B20 = 0 y1 B12 + y2 B21 + y3 B2 = 0 y1 B0 + y3 B1 + y4 B2 = 0 .. . Now from Matrix analytic method, yc+i = yc Ri , i = 0, 1, 2 . . .,where R is the minimal nonnegative solution of matrix quadratic equation R2 A2 +RA1 +A0 = 0. R is computed algorithmically, using the logarithmic reduction algorithm [12]. 3 Performance Characteristics of the System • Expected number of ordinary customers in the system. ∞ X EN1 = iyi e. i=0 • Expected number of emergency customers. ∞ X X M EN2 = jyi e. i=0 j=0 • Expected rate of departure from the system. K−1 X ∞ l+2K−1 X X Er = yi Ai e + yi Ai# e. i=0 l=0 i=l+K • Probability that server is idle with i customers i < K customers in the queue. K−1 X i 0 Pidle = yi0 e. i=0 4 Expected Service Time of a Customer The expected service time of a customer is described as follows. The service process can be considered as a continuous time Markov process on the finite state space {#} ∪ {(j, l, m), 1 ≤ j ≤ K, 1 ≤ l ≤ r, 1 ≤ m ≤ a} and {#} is the absorbing state.   0 0 Qs = , Σ Σ0   A11 A12   A21 A22   . .  0 T Σ=  .. ..  , Σ = A10 A20 . . . AK0 .    AK−11 AK−12  AK1 The service time of an arbitrary customer follows phase type distribution with 1 representation (β, Σ), where β = (β1 , β2 , . . . , βK , βK+1 ), βi = (K+1)ar αi ⊗ eT 1 P K and βK = ar i=1 αir+1 . Expected service time of an arbitrary customer is given by Est = −βΣ −1 e. 5 Expected Waiting Time of an Emergency Customer When at Least K Ordinary Customers are in Waiting The expected waiting time of an emergency customer is described as follows. Consider the process H = {(τ (t), I(t), s(t)), t ≥ 0}. Then H is clearly an ir- reducible continuous time stochastic process defined on a finite state space {(i, j, k), 1 ≤ i ≤ M, 1 ≤ j ≤ K, 1 ≤ k ≤ r} ∪ {∗}, where ∗ represents the absorption state and all other states are transient. The infinitesimal generator of this process is as shown below. C C0   QH = , 0 0 where  ∗  CM CM −K ∗   CM −1 CM −(K+1)    ..   .  C= ∗,    CK C0   ..   .  C1 Ch = diag(S1 , S2 , . . . , SK ), 1 ≤ h ≤ M, 0 S10 ⊗ αK   0 ...  .. . .  . . S20 ⊗ αK  Cg∗ =    .. ..  , 0 ≤ g ≤ M − K,   . .   0  SK−1 ⊗ αK  0 0 0 C 0 = (CM 0 0 , CM 0 T −1 , . . . , C1 ) , Ci0 = (0, 0, . . . , SK 0 ⊗ αK ), 0 ≤ i ≤ M. The waiting time of an emergency customer follows a phase type distribution 1 with representation (γ, C), where γ = (γ1 , γ2 , . . . , γK ) and γi = Kr α i ⊗ eT . Expected waiting time of an emergency customer is Eew = −γC −1 e. 6 Numerical Example We illustrate the model by considering a system with K = 2, the service times are exponentially distributed with parameters µ1 and µ2 . We take arrival rate of ordinary customers as λ and arrival rate of emergency customers as γ. The matrices defining the MMAP process takes the form D0 = (−λ − γ), D1 = (λ), D2 = (γ). Fig. 1. Variation in queue length with respect arrival rates. Fig. 2. Variation in number of customers in queue with variation in µ1 ans µ2 . Fig. 3. Idleness with λ. In Figure 1 variation in expected number of ordinary customers in the system, (EN1 ) with respect to the variations in arrival rates of ordinary and emergency customers is depicted. It is observed from Figure 1 that accumulation of ordinary customers increases in the queue as arrival rate of emergency customers increases. Figure 2 describes the variation in (EN1 ) with respect to the service rates µ1 and µ2 . Figure 3 shows that the probability that system is in idle state decreases with the increase of rate arrival of the ordinary customers. 7 Conclusion In this paper we have studied a queue with batch service and batch size de- pends the number of customers in the buffer. The model is observed in a courier delivering system with two kinds of arrivals. The strategy for batch size determi- nation is designed for reducing system cost. The model is analyzed using matrix analytic method and is illustrated with a numerical example. 8 Acknowledgment Sinu Lal T. S. thanks University Grants Commission(UGC) of India for UGC- Junior Research Fellowship(Roll no.-432693,June 2017). References 1. Banerjee, A., Gupta, U. C., Sikdar, K.: Analysis of finite buffer bulk-arrival bulk- service queue with variable service capacity and batch size-dependent service, Inter- national Journal of Mathematics in Operational Research, 5, (3), 358-386 (2013). 2. Baetens, J., Steyaert, B., Claeys, D., Bruneel, H.: System occupancy of a two class batch service queue with class dependent variable server capacity, International Conference on Analytical and Stochastic Modeling Techniques and Applications, 32-44, Springer (2016). 3. Baetens, J., Steyaert, B., Claeys, D., Bruneel, H.: Delay analysis of a two class batch service queue with class dependent variable server capacity. Mathematical Methods of Operations Research,88(1), 37-57 (2018). 4. Baetens, J., Steyaert, B., Claeys, D., Bruneel, H.: System occupancy in a multi- class batch-service queueing system with limited variable service capacity. Annals of Operations Research, 1-24 (2019). 5. Chaudhry, M. L., Gupta, U. C.: Modelling and analysis of M/Ga,b /1/n queue a simple alternative approach, Queueing Systems, 31(1-2), 95-100 (1999). 6. Gold, H., Tran-Gia, P.: Performance analysis of a batch service queue arising out of manufacturing system modelling. Queueing Systems, 14(3-4), 413-426 (1993). 7. Gupta, G., Banerjee, A.: On M/G(a;b) /1/n queue with batch size and queue length dependent service. International Conference on Mathematics and Computing, 249- 262, Springer (2018). 8. Gupta, U. C., Sikdar, K.: The fnite-buffer m/g/1 queue with general bulk-service rule and single vacation,Performance Evaluation, 57(2), 199-219 (2004). 9. He, Q-M.: Fundamentals of matrix-analytic methods, Vol. 365, Springer, New York (2014). 10. Kuppa, S., Dattatreya, G.: Modeling and analysis of frame aggregation in unsatu- rated wlans with finite buffer stations, IEEE International Conference on Commu- nications, 3, 967-972, IEEE (2006). 11. Maity, A., Gupta, U. C.: Analysis and optimal control of a queue with infinite buffer under batch size dependent versatile bulk-service rule. Oper. research, 52(3), 472-489 (2015). 12. Latouche, G., Ramaswami, V.: Introduction to matrix analytic methods in stochas- tic modeling, Vol. 5, Siam (1999). 13. Neuts, M. F.: A general class of bulk queues with Poisson input. The Annals of Mathematical Statistics, 38(3), 759–770, JSTOR (1967). 14. Neuts, M. F.: Probability distributions of phase type. Liber Amicorum Prof. Emer- itus H. Florin (1975). 15. Neuts, M. F.: Matrix-geometric solutions in stochastic models: an algorithmic ap- proach. Courier Corporation (1994). 16. Powell, W. B.: Analysis of vehicle holding and cancellation strategies in bulk arrival, bulk service queues. Transportation Science. 19(4), 352–377, INFORMS (1985). 17. Reddy, G. K., Nadarajan, R., Kandasamy, P.: A nonpreemptive priority multiserver queueing system with general bulk service and heterogeneous arrivals. Computers and operations research, 20(4), 447-453 (1993). 18. Sasikala, S., Indhira, K.: Bulk service queueing models survey. International Jour- nal of Pure and Applied Mathematics 106(6), 43-56 (2016). 19. Sikdar, K., Samanta, S.: Analysis of a finite buffer variable batch service queue with batch markovian arrival process and servers vacation. Oper. research, 53(3), 553-583 (2016). 20. Sinu Lal, T. S., Krishnamoorthy, A., Joshua, V. C.: A Multiserver Tandem Queue with a Specialist Server Operating with a Vacation Strategy. Automation and Re- mote Control, 81, 760-773 (2020). 21. Sinu Lal, T. S., Krishnamoorthy, A., Joshua, V. C. September. A Queueing In- ventory System with Search and Match-An Organ Transplantation Model. Interna- tional Conference on Distributed Computer and Communication Networks, 273-287, Springer, Cham (2019). 22. Van De Rzee, D. J., Van Harten, A., Schuur, P. C.: Dynamic job assignment heuristics for multi-server batch operations-a cost based approach. International Journal of Production Research, 35(11), 3063-3094 (1997). 23. Yu, M., Alfa, A. S.: Algorithm for computing the queue length distribution at various time epochs in DM AP/G(1,a,b) /1/n queue with batch size dependent service time. European Journal of Operational Research, 244(1), 227-239 (2015).