Next Article in Journal
Self-Management Portfolio System with Adaptive Association Mining: A Practical Application on Taiwan Stock Market
Next Article in Special Issue
Estimation of the Optimal Threshold Policy in a Queue with Heterogeneous Servers Using a Heuristic Solution and Artificial Neural Networks
Previous Article in Journal
The Higher-Order of Adaptive Lasso and Elastic Net Methods for Classification on High Dimensional Data
Previous Article in Special Issue
Graph Metrics for Network Robustness—A Survey
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Queuing-Inventory Models with MAP Demands and Random Replenishment Opportunities

by
Srinivas R. Chakravarthy
1,* and
B. Madhu Rao
2
1
Departments of Industrial and Manufacturing Engineering & Mathematics, Kettering University, Flint, MI 48504, USA
2
Department of Business Systems and Analytics, School of Business Administration, Stetson University, Deland, FL 32723, USA
*
Author to whom correspondence should be addressed.
Mathematics 2021, 9(10), 1092; https://doi.org/10.3390/math9101092
Submission received: 18 April 2021 / Revised: 30 April 2021 / Accepted: 10 May 2021 / Published: 12 May 2021
(This article belongs to the Special Issue Distributed Computer and Communication Networks)

Abstract

:
Combining the study of queuing with inventory is very common and such systems are referred to as queuing-inventory systems in the literature. These systems occur naturally in practice and have been studied extensively in the literature. The inventory systems considered in the literature generally include ( s , S ) -type. However, in this paper we look at opportunistic-type inventory replenishment in which there is an independent point process that is used to model events that are called opportunistic for replenishing inventory. When an opportunity (to replenish) occurs, a probabilistic rule that depends on the inventory level is used to determine whether to avail it or not. Assuming that the customers arrive according to a Markovian arrival process, the demands for inventory occur in batches of varying size, the demands require random service times that are modeled using a continuous-time phase-type distribution, and the point process for the opportunistic replenishment is a Poisson process, we apply matrix-analytic methods to study two of such models. In one of the models, the customers are lost when at arrivals there is no inventory and in the other model, the customers can enter into the system even if the inventory is zero but the server has to be busy at that moment. However, the customers are lost at arrivals when the server is idle with zero inventory or at service completion epochs that leave the inventory to be zero. Illustrative numerical examples are presented, and some possible future work is highlighted.

1. Introduction

Models for inventory management under uncertainty have been studied extensively. The two key questions of when and how many to order have been addressed under a variety of factors such as the nature of inventory review (continuous or periodic), order quantity (fixed or variable), lead time for an order to be fulfilled (negligible, constant or random), nature of demand (deterministic or random), and other factors to optimize a function of various costs such as ordering, carrying inventory, lost sales, etc.
Most models assume a single supplier and fixed cost of replenishment. Some models incorporate the availability of random opportunities for replenishment which may lower system costs due to reduced unit cost and/or ordering cost. We refer to them as opportunistic replenishment. Friend [1] studied systems with special opportunities occurring according to a Poisson process. These opportunities, which may be exercised only at the instant of their occurrence, are offered for the same unit price but at reduced or zero ordering cost. Hurter and Kaminsky [2] extend this model to systems where the special opportunities may be exercised during a random period. Silver, Robb, and Rahnama [3] developed an efficient heuristic for the Hurter and Kaminsky [2] model. Gurnani [4] considered periodic review systems and used game theory to study the economics of placing a special order between reviews if a discounted sale is offered. Moinzadeh [5] considered systems with constant demand rate, zero lead time, and special discounted opportunities occurring at exponentially distributed intervals, and obtained optimal order quantity at regular price when inventory reaches zero, and the threshold and order quantity at discounted price. Feng and Sun [6] suggest modifying the ( s , S ) policy to a four-parameter system (threshold and order quantity for regular and discount purchases) and proposed a bisection search procedure to determine the optimal values of the four parameters. Goh and Sharafali [7] consider the model in [1,2] and incorporate the policy of passing the cost reduction due to special purchases downstream to increase demand. Chaouch [8] considers the model in [1,2] when special and regular prices are valid over alternating exponential intervals. Tajbakhsh and Zolfaghari [9] consider systems with discounts offered at exponential intervals with the discount price given by a discrete random variable and develop optimal order quantities at each special price and further extended the model to the case of multiple suppliers. Karimi-Nasab and Konstantaras [10] consider a system with constant demand, periodic random review intervals (uniformly distributed or exponential subject a maximum and minimum) and random special sale offers, and determine maximum inventory level for regular purchases, and a higher maximum for special purchases. Den Boer and Zwart [11] consider a system where the management makes simultaneous decisions on whether to take advantage of the special discounted offer, and the selling price of a unit similar to [7]. All the above references assume that the lead time for receiving a special replenishment is negligible.
In all the above models, it is assumed that the special opportunities for replenishment are always considered as a supplement to the normal replenishment process. In many situations such as drugstores, groceries, small supermarkets, etc., the suppliers visit the retailers at random (but frequent) intervals to offer special sales. This raises the possibility that for some systems it may be more economical to manage inventory solely based on special offers. This is particularly attractive for non-critical items (e.g., canned goods, generic medication) where stock outs are not critical. Special replenishment opportunities offer lower unit cost and/or reduced ordering cost, and are usually available without delay.
Inventory models discussed above assume that the time to process a demand is negligible. In many situations, completing a customer’s service requires other resources, in addition to the item(s) in inventory (e.g., requiring a pharmacist to process the sale for prescription medication). Such inventory systems are referred to as queuing-inventory (QI) models or inventory models with positive service time, and have received a great deal of attention recently. Research in queuing-inventory systems may be classified based on features such as the nature of customer arrival process, service time distribution, server vacation, service discipline, customer behavior, inventory review policy, replenishment policy, stock-out assumptions, perishability of units, lead time for replenishment, among others. We refer the readers to Krishnamurthy, Shajin, and Narayanan [12] for a general review of queuing-inventory models, to Choi et al. [13] for a survey of QI models with phase-type service times, and to Melikova and Shahmaliyeva [14] for a discussion of QI models for perishable items. Recent generalizations of QI models are due to Chakravarthy, Maity and Gupta [15], who studied a QI model in which the service is carried out in batches, Chakravathy [16], who studied a QI model in which the demand occurs in batches, and Chakravarthy and Rumyantsev [17], who study a QI system with batch demand, Markovian Arrival Process (MAP), and phase-type service time and exponential replenishment lead times.
In this paper, we explore queuing-inventory systems which solely use special opportunities for replenishment. This research is at the interface of the two areas of traditional inventory systems, and the queuing-inventory systems. Most traditional inventory models consider the demand to occur one unit at a time, at either constant or exponential intervals. Our models extend this to consider demand that can occur in batches of random size, with arrivals following a very general process which allows for a broad class of inter-arrival times and autocorrelation between inter-arrival times. For both types of inventory systems, it is a significant departure to manage inventory exclusively based on random replenishment opportunities. The objective of this study was to understand the behavior of such systems and compare them with the traditional ( s , S ) -type inventory management systems.

2. General Description

The customer arrivals are modeled using Neuts’ versatile Markovian point process, now referred to as batch Markovian arrival process ( B M A P ). If customers arrive one at a time, the process is referred to as Markovian arrival process or briefly M A P . While the customers arrive singly to the service facility, but their demands for inventory items are in varying sizes. One of the major advantages of using M A P to model the arrival process is to capture any correlation present in successive arrivals. The M A P is completely described by two parameter matrices, say, ( E 0 , E 1 ) of dimension m. While, E 0 governs transitions within the underlying Markov process with generator E = E 0 + E 1 , the transitions governed by E 1 correspond to the arrivals of customers. We assume that E is an irreducible generator. M A P and B M A P have been extensively used in stochastic model and there is a vast literature available on these. We refer the reader to [18,19,20,21,22,23,24,25,26,27,28,29] for details on M A P and B M A P and other key references.
Each customer demands a random number of items from inventory before being processed by a single server. At the instant of an arrival, the available inventory is reduced by the amount demanded by the customer. The number of remaining items in inventory is referred to as available to promise or ATP in the business community. In this paper, we always refer to ATP when referring to inventory. Physical inventory is not relevant because once allocated to a customer, a unit is no longer available.
When the inventory level is positive, but less than the number of units demanded by a customer, the customer’s demand is partially satisfied and the inventory level is set to zero. Arrivals to the system when inventory level is zero, are denied entry and are lost. The time to service a customer has a phase-type or P H -distribution and is independent of the number of units demanded. For details on P H -distributions and their usefulness in stochastic modeling, we refer the reader to [29,30,31].
Opportunities to replenish the inventory occur according to a Poisson process. We assume that the lead time associated with a replenishment is zero. Replenishment decisions are based on two parameters K and L ( 0 L < K < ) and for convenience we refer to the policy as a ( K , L ) policy. If the inventory at the time a replenishment opportunity becomes available is equal to i, a replenishment order to bring the inventory level to K is placed with probability 1 if i L , and with with probability a i L if L < i < K .
In this paper, we consider two types of models, referred to as Opportunistic Model 1 and Opportunistic Model 2, and these are described in the following sections.
Any queuing-inventory model with positive service time can be looked at two different ways according to these two models to see which one will be better from either customer or management or from both points of view. Therefore, there is a trade-off between the two models and a comparison of these would benefit practicing managers to make a decision on the choice of the models.
We use the classical matrix-analytic method (see, e.g., [30]) to formulate the associated stochastic process with matrix-geometric steady-state probability vector. Performance measures of interest are obtained in terms of the steady-state probabilities.
The structure of the present paper is as follows. Models 1 and 2 are analyzed in Section 3 and Section 4, respectively. Some selected system performance measures are presented in Section 3.2 and Section 4.2. In Section 5, we numerically study the key performance measures of both models in steady-state as well as present a cost analysis to compare the two models. Some concluding remarks are noted in Section 6. The following notation is used in the paper.
  • The subscript i , i = 1 , 2 refers to the model under consideration.
  • The symbol stands for the transpose notation.
  • The notations ⊗ and ⊕ stand for the Kronecker product and Kronecker sum, respectively (see, e.g., [32,33]).
  • e = ( 1 , 1 , , 1 ) , whose dimension should be clear in the context. Where more clarity is needed, the dimension will be mentioned, e.g., e ( m ) is a column vector of 1s of dimension m.
  • e i = ( 0 , 0 , , 1 , 0 , , 0 ) , where 1 is in the ith position.
  • I denotes an identity matrix, whose dimension is dictated by the context.
  • Δ ( E 1 , , E r ) denotes a diagonal matrix with diagonal (possibly block) entries given by E 1 through E r . In the context where this notation is used, it will be clear whether the entities are scalars or vectors or matrices.

3. Opportunistic Model 1

In this first model, we consider the scenario where the customers arrive to a service facility consisting of a single server. The arrivals are modeled using Neuts’ versatile Markovian point process, now referred to as batch Markovian arrival process ( B M A P ). While the customers arrive singly to the service facility, their demands for inventory items are in varying sizes. The M A P parameter matrices describing the arrivals of customers are ( E 0 , E 1 ) of dimension m. The customers’ demands are in batches of varying sizes. Suppose that M denotes the demand size. We assume that M has the following discrete distribution on finite support,
P ( M = i ) = p i , 1 i N , i = 1 N p i = 1 .
In order to determine the arrival rate of the customers to the system, we first define the invariant vector of the generator E = E 0 + E 1 that governs the underlying Markov process for the arrivals. Suppose that θ is the invariant probability vector of E. That is,
θ E = 0 , θ e = 1 .
Then, the arrival rate of the customers is given by
λ = θ E 1 e .
Let μ M to be the mean of the demand size of each customer. That is,
μ M = i = 1 N i p i .
In this paper, we assume that the demands of the customers are either fully met or partially met or not met at all, depending on, respectively, if the demands in the inventory are fully adequate, partially adequate, or the inventory level is zero. When the inventory level is zero at the time of the arrival of a customer, the arriving customer will be lost.
The size of the inventory is assumed to be finite with a capacity to hold at most K items. As described and motivated earlier, the main purpose of this paper was to incorporate opportunistic replenishment into stochastic models dealing with queues with inventory and positive services. The opportunistic events occur according to a Poisson process with parameter γ . Therefore, when an opportunity occurs for replenishing the inventory, the system has the ability to bring the current inventory level to the maximum level through a given probability structure. Suppose that the inventory level at the time a replenishment opportunity is i. Then, the opportunity will be availed to bring the inventory level to K with probability 1 if i L ; however, if i > L , then with probability a i L , 0 a i L < 1 , the opportunity will be availed to bring the inventory to K, for L < i K 1 .
When a customer enters the system (either with fully or partially met demands), the inventory level will be reduced by that customer’s demand size. Note that the only way an arriving customer will be lost is when the inventory level is zero at that moment. However, customers have to wait until their demands are processed with a positive service time, which is assumed to be random. The services follow a phase-type ( P H ) -distribution with representation ( β , S ) of order n. Note that the mean service rate, say, μ , is given by
μ = β ( S ) 1 e 1 .
For use in sequel, we define S 0 = S e .
Suppose we define the quantities:
  • J 1 ( t ) to be the number of customers in the system (including one in service),
  • J 2 ( t ) to be the level of the inventory,
  • J 3 ( t ) to be the phase of the service (if the server is idle, this will not be defined),
  • J 4 ( t ) to be the phase of the arrival process,
at time t. The Markov process { ( J 1 ( t ) , J 2 ( t ) , J 3 ( t ) , J 4 ( t ) ) : t 0 } is a quasi-birth-and-death process ( Q B D ) with state space given by
Ω 1 = { ( 0 , j , r ) : 0 j K , 1 r m } { ( ( i , j , k , r ) : 0 j K , 1 k n , 1 r m , i 1 } .
We now define the levels based on the set of states defined above. Suppose that 0 = { ( 0 , j , r ) : 0 j K , 1 r m } and i = { ( i , j , k , r ) : 0 j K , 1 k n , 1 r m } , i 1 . This way, we notice that 0 corresponds to an idle system while i corresponds to the system having i customers with one in service, and the inventory, and the arrival process are in various states. For use in the sequel, we define a few additional notations.
  • p ^ i , probabilities of demands greater than i, are computed as
    p ^ i = k = i + 1 N p k , 1 i N 1 .
  • P, a square matrix of dimension K + 1 is defined as
    P = 0 0 0 0 0 0 1 0 0 0 0 0 p ^ 1 p 1 0 0 0 0 p ^ 2 p 2 p 1 0 0 0 p ^ N 1 p N 1 p N 2 0 0 0 0 p N p N 1 0 0 0 0 0 p N 0 0 0 0 0 0 p 2 p 1 0 .
    A quick look at P indicates that the customers’ demands are taken into consideration at the time of their arrivals. The first column of P justifies that the structure due to the customers’ demands are met partially. It is easy to verify that
    P e = ( 0 , 1 , , 1 ) .
The generator, Q 1 , of the Q B D -process for Opportunistic Model 1 is given by (note that blank spaces correspond to the entries are zero)
Q 1 = B P β E 1 I S 0 I A 1 A 0 A 2 A 1 A 0 A 2 A 1 A 0 .
For a better display of the matrices, we take
F = E γ I , F 0 = E 0 γ I , F r = E 0 a r γ I , 1 r K L 1 , S E = ( S E ) γ I , S E 0 = ( S E 0 ) γ I , S E r = ( S E 0 ) a r γ I , 1 r K L 1 .
The matrices appearing in Q 1 are as follows.
Mathematics 09 01092 i001
Mathematics 09 01092 i002
and
A 0 = P I E 1 , A 2 = I S 0 β I .
Suppose that η = ( η 0 , , η K ) denotes the invariant probability vector of A = A 0 + A 1 + A 2 , that is, η A = 0 and η e = 1 . The following lemma, which is intuitively clear, is very useful in numerical computation as well as for probabilistic interpretation.
Lemma 1.
The vector η is such that
i = 0 K η i = ( μ β ( S 1 ) θ ) .
Proof. 
Setting A = A 0 + A 1 + A 2 , S ˜ = S + S 0 β , and noting that A ( e I I ) = [ e ( S ˜ E ) ] , we see that
η A = 0 η A ( e I I ) = 0 η [ e ( S ˜ E ) ] = 0 ,
which implies that
η ( e I ) [ S ˜ E ] = 0 .
On noting that the the invariant vectors of S ˜ and E, are, respectively, given by μ β ( S 1 ) and θ , the stated result follows from Equation (16).
The computation of η is done by exploiting the structure of the matrices appearing in A. From the definition of η and the equations associated with the invariant vector, we get
η 0 = j = 1 N η j p j , 0 ( I E 1 ) [ γ I ( S ˜ E ) ] 1 , η k = j = k + 1 N + k η k p j , k ( I E 1 ) [ γ I ( S ˜ E 0 ) ] 1 , 1 k L , η k = j = k + 1 N + k η k p j , k ( I E 1 ) [ a k L γ I ( S ˜ E 0 ) ] 1 , L + 1 k K 1 , η K = γ [ j = 0 L η k + j = L + 1 K 1 a j L η k ] ( ( S ˜ E 0 ) ] 1 .
The following lemma, which again is intuitively clear, gives necessary and sufficient condition for the stability of the model under study. □
Lemma 2.
The queuing-inventory model with opportunistic replenishment with the generator given in (9) is stable if and only if
λ η 0 ( e E 1 e ) < μ .
Proof. 
From the celebrated Neuts’ ergodicity criterion (see, e.g., [30]) for stability of Q B D -process: π A 0 e < π A 2 e , the stated result follows immediately after some elementary matrix manipulations.
Letting x = ( x 0 , x 1 , x 2 , ) , with the vector x 0 of dimension ( K + 1 ) m and each of the vectors x i , i 1 , of dimension ( K + 1 ) m n , be the steady-state probability vector of Q ( 1 ) process under consideration, we notice the following:
x Q 1 = 0 , x e = 1 .
The following theorem is a direct consequence of Neuts’ result on Q B D -process (see, e.g., [30]) and we state it here for the sake of completeness. □
Theorem 1.
Under the stability condition given in (18), the steady-state probability vector, x, of Q 1 is of matrix-geometric type:
x i = x 1 R i 1 , i 1 ,
where the (rate) matrix R is the minimal non-negative solution to:
R 2 A 2 + R A 1 + A 0 = 0 .
and the vectors x 0 and x 1 are obtained by solving
x 0 B + x 1 ( I S 0 I ) = 0 , x 0 ( P β E 1 ) + x 1 ( A 1 + R A 2 ) = 0 ,
subject to:
x 0 e + x 1 ( I R ) 1 e = 1 .

3.1. Computation of R

Without any additional special structure to the rate matrix, R, it is pertinent that we compute it numerically. Of course, any sparsity in the coefficient matrices seen in the matrix-quadratic equation given in (21) should be fully exploited (see, e.g., [20,21,30,34]), especially, when m , n , L , and K are sufficiently large. In addition, the fact that R A 2 e A 0 e , will help to identify if there are any zero rows in R further simplifying the computational steps. For example, let V = R 2 and partition R into matrices of dimension m n as follows.
R = R 0 , 0 R 0 , 1 R 0 , K R 0 , 0 R 0 , 1 R 0 , K R K , 0 R K , 1 R K , K .
Noting that for the model under study R 0 , j , 0 j K = 0 , we can exploit that observation. Now, the equation in (21) with the help of the partition of R can be rewritten for numerical implementation as follows, noting that, for 1 i K ,
V i , 0 ( S 0 β I ) + R i , 0 ( S E γ I ) + p i , 0 ( I E 1 ) = 0 , V i , j ( S 0 β I ) + R i , j ( S E 0 γ I ) + p i , j ( I E 1 ) = 0 , 1 j L , V i , j ( S 0 β I ) + R i , j ( S E 0 ) + p i , j ( I E 1 ) = 0 , L + 1 j K 1 , V i , K ( S 0 β I ) + R i , K ( S E 0 ) + γ R i , L + p i , K ( I E 1 ) = 0 .
If m and n are very large, one can further exploit the coefficient matrices seen in (25). For example, the Kronecker products and sums appearing in the Equation (25) can be exploited so as to deal only with matrices of dimension m. Suppose that
ζ i , j = lim t P ( J 2 ( t ) = i , J 4 ( t ) = j ) = , 0 i K , 1 j m .
Recall that J 2 ( t ) , and J 4 ( t ) are, respectively, track of the number in the inventory and the phase of the arrival process. Let ζ = ( ζ 0 , ζ 1 , , ζ K ) be such that ζ i , 0 i K , of dimension m gives the (marginal) probability vector of seeing i in the inventory with the phase of the arrival process being in various phases in steady-state.
Lemma 3.
The (marginal) probability vector ζ is independent of the service time distribution.
Proof. 
Observe that the steady-state equations given in (36) when rewritten yield
x 0 B + x 1 ( I S 0 I ) = 0 , x 0 ( P β E 1 ) + x 1 A 1 + x 2 A 2 = 0 , x i 1 A 0 + x i A 1 + x i + 1 A 2 = 0 , i 2 .
Suppose we define H = B + P E 1 , a square matrix of dimension ( K + 1 ) m n , and partition it into blocks of dimension m n × m n as
H = H 0 , 0 H 0 , 1 H 0 , K H 1 , 0 H 1 , 1 H 1 , K H K , 0 H K , 1 H K , K .
From the definitions of B (see Equation (11)) and P (see Equation (7)), it is easy to identify the (block) elements of H.
It is easy to verify that A ( I e I ) is given by
e H 0 , 0 e H 0 , 1 e H 0 , K e H 1 , 0 e H 1 , 1 e H 1 , K e H K , 0 e H K , 1 e H K , K .
By post-multiplying the second and third equations in (27) by I e I and adding the resulting equations with the first equation in (27), we get
x 0 H + i = 1 x i A ( I e I ) = 0 ,
which, on observing that
ζ = x 0 + i = 1 x i ( I e I ) ,
implies that
ζ H = 0 .
The above equation yields the stated result as the matrix H is independent of the service time distribution. □
Remark 1.
Lemma 3 is intuitively clear. This is due to the fact that the inventory level is decreased at arrival times and the service time plays no role.

3.2. Selected System Measures in Steady-State

For qualitative evaluation of the models presented in this paper, we look at the system performance measures in the following table. For Opportunistic Model 1, the measures along with their formulas are as follows.
1.
Server idle probability: ν 1 = x 0 e .
2.
Probability of idle server with positive inventory: ν 1 I = x 0 e x 0 , 0 e .
3.
Percent of server idle time with positive inventory: ν 1 * = ν 1 I / ν 1
4.
Mean number of customers in the system: μ 1 = x 1 ( I R ) 2 e .
5.
Variance of the number of customers in the system:
σ 1 2 = 2 x 1 ( I R ) 3 e μ 1 ( 1 + μ 1 ) .
6.
Probability of customer loss (at arrivals): ϑ 1 = 1 λ i = 0 x i , 0 ( e E 1 e ) .
7.
Mean inventory level: μ ^ 1 = j = 1 K j x 0 , j e + i = 1 x i , j e .
8.
Variance of the inventory level: σ ^ 1 2 = j = 1 K j 2 x 0 , j e + i = 1 x i , j e μ ^ 1 2 .
9.
Mean cycle time of replenishment:
κ 1 = γ j = 0 L x 0 , j e + j = L + 1 K 1 a j L i = 1 x i , j e 1
10.
Mean replenishment quantity:
Γ 1 = γ j = 0 L ( K j ) x 0 , j e + i = 1 x i , j e + j = L + 1 K 1 ( K j ) a j L x 0 , j e + i = 1 x i , j e .
11.
Probability of procuring an order when a replenishment opportunity arises:
ξ 1 = j = 0 L x 0 , j e + i = 1 x i , j e + j = L + 1 K 1 a j L x 0 , j e + i = 1 x i , j e .

4. Opportunistic Model 2

In the Opportunistic Model 1, we assumed that the customers are lost at arrivals when the inventory level is zero. However, in Opportunistic Model 2, we will admit customers into the system even if the inventory is zero as long as the server is busy serving. If the server is idle with zero inventory, then the arriving customers will be lost. In this model, we assume that the inventory level is reduced by the amount the customers’ demand is either fully or partially met just before the service starts (or equivalently just after a service completion). Thus, there is a possibility that the customers may be lost at service completion epochs when the inventory level is zero. All the customers waiting at the time of a service completion which results in zero inventory will be lost. This is due to the fact that a service requires at least one inventory. Arrivals, services, and the probabilistic structure for availing the opportunistic events are all the same as in the Opportunistic Model 1.
The Opportunistic Model 2 is of the G I / M / 1 -type [30]. Suppose that we define
  • J ^ 1 ( t ) to be the number of customers in the system,
  • J ^ 2 ( t ) to be the level of the inventory,
  • J ^ 3 ( t ) to be the phase of the service (if the server is idle, this will not be defined),
  • J ^ 4 ( t ) to be the phase of the arrival process,
at time t. The stochastic process { ( J ^ 1 ( t ) , J ^ 2 ( t ) , J ^ 3 ( t ) , J ^ 4 ( t ) ) : t 0 } is a Markov process with the state space given by
Ω 2 = { ( 0 , j , r ) : 0 j K , 1 r m } { ( ( i , j , k , r ) : 0 j K , 1 j n , 1 r m , i 1 } .
We define the set of states similar to Opportunistic Model 1. The generator, Q 2 , of the system governing Model 2 is of the form
Q 2 = B P β E 1 I S 0 I G 1 G 0 C G 2 G 1 G 0 C G 2 G 1 G 0 .
The matrix B is as given (11) and the matrices C , G 0 , G 1 , and G 2 appearing in Q 2 are as follows.
C = e 1 ( K + 1 ) e 1 ( K + 1 ) S 0 I , G 0 = I D 1 , G 2 = P S 0 β I ,
and
Mathematics 09 01092 i003
where S E r , 0 r K L 1 are as defined in (10). Note that the level 0 can be reached from any other level and hence this queuing model can be thought as a catastrophic model, and hence is always stable. In addition, observe that when γ the model becomes unstable when λ μ .
Suppose that y, partitioned as y = ( y 0 , y 1 , y 2 , ) , satisfies
y Q 2 = 0 , y e = 1 .
The following theorem is a direct consequence of Neuts’ result (see, e.g., [30]) and for the sake of completeness, we will register it here.
Theorem 2.
Assuming that γ is finite, the steady-state probability vector, y, of Q 2 is of matrix-geometric type and is given by:
y i = y 1 R ^ i 1 , i 1 ,
where the (rate) matrix R ^ is the minimal non-negative solution to:
R ^ 2 G 2 + R ^ G 1 + G 0 = 0 .
and the vectors y 0 and y 1 are obtained by solving
y 0 B + y 1 [ ( I S 0 I ) + R ^ ( I R ^ ) 1 C ] = 0 , y 0 ( P β E 1 ) + y 1 [ G 1 + R ^ G 2 ] = 0 ,
subject to
y 0 e + y 1 ( I R ^ ) 1 e = 1 .
The rate matrix R ^ can be efficiently computed by exploiting the structure of the coefficient matrices and the steps are similar to the computation of R seen in the Opportunistic Model 1 and hence will be omitted. Again, the exploitation of the structure of the coefficient matrices will result in dealing with matrices of smaller dimensions.

4.1. Computation of y 0 and y 1

Due to the special structure of the coefficient matrices appearing in the solution of the vectors, y 0 and y 1 , it is worth to exploit them as follows. First, we partition these vectors as follows.
y 0 = ( y 0 , 0 , , y 0 , K ) , y 1 = ( y 1 , 0 , , y 1 , K ) .
Note that the vectors, y 0 , j , 0 j K , have dimension m, whereas the vectors, y 1 , j , 0 j K , have dimension m n .
y 0 , 0 = y 1 ( I R ^ ) 1 S 0 I 0 ( γ I E ) 1 , y 0 , j = k = 1 n y 1 , j , k S k 0 ( γ I E 0 ) 1 , 1 j L , y 0 , j = k = 1 n y 1 , j , k S k 0 ( a j L γ I E 0 ) 1 , L + 1 j K 1 , y 0 , K = γ k = 0 L y 0 , k + γ k = L + 1 K 1 a k L y 0 , k + k = 1 n y 1 , j , k S k 0 S k 0 ( E 0 ) 1 .
y 1 , j = k = 0 K p k , j ( β E 1 ) + k = 0 K y 1 , k D k , j [ γ I ( S E 0 ) ] 1 , 0 j L , y 1 , j = k = 0 K p k , j ( β E 1 ) + k = 0 K y 1 , k D k , j [ a j L γ I ( S E 0 ) ] 1 , L + 1 j K 1 , y 1 , K = γ k = 0 L y 1 , k + γ k = L + 1 K 1 a k L y 1 , k + k = 0 K y 1 , j , k D k , j [ ( S E 0 ) ] 1 ,
and the normalizing constant is given by
k = 0 K y 0 e + k = 0 K y 1 b = 1 ,
where the block entries, D j , k , 0 j , k K , of the matrix D and the vector b are defined as
b = ( I R ^ ) 1 e , R ^ G 2 = D .
The vector b = ( b 0 , , b K ) is obtained as follows. Suppose that ( I R ^ ) b = e ( ( K + 1 ) m n ) , then we have
b i = ( I R ^ i , i ) 1 e ( m n ) + j = 0 , j i K R ^ i , j b .
The block entries of D are computed as follows
D j , k = k = 0 K p k , j R ^ i , k ( S 0 β I ) , 0 j , k K .
Note that the quantity
( I R ^ ) 1 S 0 I 0 = ( I R ^ ) 1 S 0 e 1 ( m ) S 0 e 2 ( m ) S 0 e m ( m ) 0 0 0 ,
can be obtained by solving for d j , 1 j m , which satisfies
( I R ^ ) d j = S 0 e j ( m ) 0 ,
and one can use the same type of equations seen in (44) by replacing e ( m n ) with
S 0 e j ( m ) 0 .

4.2. Selected System Measures in Steady-State

For qualitative evaluation of the models presented in this paper, we look at the system performance measures in the following table. For the Opportunistic Model 1, the measures along with their formulas are as follows.
1.
Server idle probability: ν 2 = y 0 e .
2.
Probability of idle server with positive inventory: ν 2 I = y 0 e y 0 , 0 e .
3.
Percent of server idle time with positive inventory: ν 2 * = ν 2 I / ν 2 .
4.
Mean number of customers in the system: μ 2 = y 1 ( I R ^ ) 2 e .
5.
Variance of the number of customers in the system:
σ 2 2 = 2 y 1 ( I R ^ ) 3 e μ 2 ( 1 + μ 2 ) .
6.
Probability of customer loss at arrivals:
ϑ 2 = 1 λ y 0 , 0 ( e E 1 e ) + i = 1 ( i 1 ) y i , 0 ( S 0 e ) .
7.
Mean inventory level: μ ^ 2 = j = 1 K j y 0 , j e + i = 1 y i , j e .
8.
Variance of the inventory level: σ ^ 2 2 = j = 1 K j 2 y 0 , j e + i = 1 y i , j e μ ^ 2 2 .
9.
Mean cycle time of replenishment:
κ 2 = γ j = 0 L y 0 , j e + j = L + 1 K 1 a j L i = 1 y i , j e 1 .
10.
Mean replenishment quantity:
Γ 2 = γ j = 0 L ( K j ) y 0 , j e + i = 1 y i , j e + j = L + 1 K 1 ( K j ) a j L y 0 , j e + i = 1 y i , j e .
11.
Probability of procuring an order when a replenishment opportunity arises:
ξ 2 = j = 0 L y 0 , j e + i = 1 y i , j e + j = L + 1 K 1 a j L y 0 , j e + i = 1 y i , j e .
Note that the loss probability, ϑ 2 , has two components, the first is for the loss at arrivals ( ϑ 2 a ) and the second one is for the loss at a service completion due to zero inventory ( ϑ 2 d ) , which makes the server become idle and hence it cannot serve any customer.

5. Illustrative Numerical Examples

First, we describe the experimental setup. We start with considering the dependence of relative difference of performance measures across various batch size distributions.
We define the following arrival processes (along with their parameters including the standard deviations):
ERA
Erlang distributed inter-arrival times with density λ ˜ k x k 1 e λ ˜ x / ( k 1 ) ! , λ ˜ = 4 , and k = 4 (standard deviation 0.5774 ).
HEA
Hyperexponential inter-arrival times with density i = 1 k p i λ ˜ i e λ ˜ i x , k = 4 , λ ˜ i = ( 63.1 ) · 10 i , i = 1 , 2 , 3 , and p = ( 0.6 , 0.25 , 0.10 , 0.05 ) (standard deviation 4.9629 ).
PCA
Markov arrival process ( E 0 , E 1 ) with positive correlation 0.4637 (standard deviation 1.3153 ), where
E 0 = 1.05 1.05 0 0 1.05 0 0 0 10.5 , E 1 = 0 0 0 1.035 0 0.0105 0.105 0 10.395 .
The parameters of distributions are normalized so as to obtain a unit arrival rate. The graphs of the probability density function of the three inter-arrival times are plotted in Figure 1 and the joint probability function of the M A P with positive correlation is plotted in Figure 2.
We also use the following service time distributions:
ERS
Erlang distributed service times having density μ ˜ k x k 1 e μ ˜ x / ( k 1 ) ! with k = 3 .
EXS
Exponentially distributed service times of rate μ .
HES
Hyperexponentially distributed service times with density i = 1 k p i μ ˜ i e μ ˜ i x , k = 3 , p = ( 0.7 , 0.25 , 0.05 ) and μ ˜ i = ( 8.2 μ ) · 10 i , i = 1 , , 3 .
The service rate μ will be normalized so as to arrive at a given service rate. Obviously, the above three P H -distributions are qualitatively different covering a wide range applicable in practice. The plot of these three distributions is given in Figure 3.
For our illustrative examples, here we take the batch size distribution to be uniform in { 1 , , N } and the batch size maximum to be N, which is set at N = 7 . It should be pointed out that we did consider other batch size distributions like constant, truncated Poisson, and truncated geometric. For the range of parameters considered, we observed that numerical results with other batch size distribution are similar to those with uniform distribution.
When a replenishment opportunity occurs, an order is placed with probability 1 when the inventory level is at most L, and with probability a i , when the inventory level is i , L + 1 < i K 1 . Order quantity is always determined to bring the inventory level to K. We considered three types of probability functions for a i , namely, constant, linearly decreasing, and non-linearly decreasing for i = L + 1 i K 1 . We observed that the effect of the type of probability function used for a i on the system performance is minimal. In many cases, the corresponding results for the three probability functions are equal within the fourth decimal.
Queuing-inventory systems consist of two interacting subsystems, namely the inventory subsystem (ISS) and the customer subsystem (CSS).
  • In ISS, inventory is consumed and replenished as needed and is characterized by the parameters K , L , γ and the distribution of the time between two opportunities for replenishment. The arrival process to CSS impacts ISS through the demand for items. μ ^ 1 [ μ ^ 2 ], σ ^ 1 [ σ ^ 2 ], Γ 1 [ Γ 2 ], ξ 1 [ ξ 2 ], κ 1 [ κ 2 ] represent the measures of performance for ISS for Model 1 [Model 2].
  • In CSS, customers arrive, receive service (plus items from inventory), and depart and are characterized by the arrival and service processes. Some arrivals are lost due to lack of inventory at the time of arrival. In Model 2, customers may also be lost at a service completion epoch due to lack of inventory at that moment. Customer loss is impacted by the availability of inventory in ISS. μ 1 [ μ 2 ] , σ 1 [ σ 2 ] , ν 1 [ ν 2 ] , ν 1 * [ ν 2 * ] , ϑ 1 [ ϑ 2 a , ϑ 2 d ] represent measures of performance for CSS for Model 1 [Model 2].
Inventory levels control the interaction between the two subsystems. If the inventory levels are high [low], customer loss and the interaction between ISS and CSS will be low [high]. In the limiting case when infinite inventory is maintained, interaction between the two subsystems decreases to zero and Model 2 converges to Model 1. In this paper, we consider ϑ 1 [ ϑ 2 ] and μ 1 [ μ 2 ] as good measures for describing the extent of interaction between ISS and CSS for the two models.
Table 1 and Table 2 respectively summarize the results for CSS and ISS for a broad range of parameter values. These tables represent a subset of the numerical results used in the following qualitative summary of the effect of various parameters on the performance of two subsystems for the two models. Table 3 compares the system measures for the two models. In these tables, the abbreviations Arr. and Ser. are used to identify the arrival process and service time distribution.

5.1. Impact of Service Rate and Service Time Distribution

In Model 1, customers do not enter the system when inventory level is zero, and when a customer enters the system the inventory level is immediately reduced to reflect the new customer’s demand. As a consequence, the service rate and service time distribution do not impact ISS. They do impact ISS in Model 2 because customers who arrive when inventory level is zero wait in line hoping for a replenishment during their wait. However, this impact is minimal because the likelihood of a replenishment during their wait is typically very small.
As it might be expected, the performance measures for CSS improved with increasing service rate in both models. The overall impact of service time variability is to degrade the performance of the system, with HES generating the strongest impact. Specifically, increasing service rate increased the probability of an idle server ( ν 1 , ν 2 ) and decreased the number of customers in the system ( μ 1 , μ 2 ). Increasing variability of service time distribution, on the other hand, increased the mean ( μ 1 , μ 2 ), the standard deviation ( σ 1 , σ 2 ), and the coefficient of variation of the number of customers in the system, and increased the percent of time the server is idle with positive inventory ( ν 1 * , ν 2 * ) . The probability of customer loss in Model 1 ( ν 1 ) was unaffected by the service time variability, but ν 2 increased marginally. We also observed that when the service rate is less than the arrival rate, ν 1 < ϑ 1 and ν 2 < ϑ 2 ; and when service rate is greater than the arrival rate, ν 1 > ϑ 1 and ν 2 > ϑ 2 .

5.2. Effect of Arrival Process

Arrival process impacts both the ISS and the CSS and in general, increasing variability of arrivals had the effect of degrading the performance of the system. It is interesting to note that for systems with both a highly variable arrival process and service time distribution, the combined effect is more attenuated than the individual effects.
Systems with highly variable arrival processes (HEA and PCA relative to ERA) had higher probability of a customer loss ( ϑ 1 and ϑ 2 ), and higher mean ( μ 1 , μ 2 ), standard deviation ( σ 1 , σ 2 ), and the coefficient of variation of the number of customers in the system. Similarly, HEA and PCA produced higher mean ( μ ^ 1 , μ ^ 2 ), standard deviation ( σ ^ 1 , σ ^ 2 ) but lower coefficient of variation of the inventory levels. Mean cycle times ( ξ 1 , ξ 2 ) were lower but the order quantities ( Γ 1 , Γ 2 ) remained relatively unaffected. HEA and PCA increased the probability that the server is idle ( ν 1 , ν 2 ) as well as the percent of time the server is idle with inventory of units on hand ( ν 1 * , ν 2 * ), significantly degrading the system performance. The impact of HEA was more pronounced than PCA.
Figure 4 and Figure 5 summarize the effect of the arrival and service processes on the steady-state marginal distributions of the number of customers in the system. In these figures, X and Y axes are truncated at much smaller values than required by the results in order to highlight the differences in the probability distributions.
Figure 4 and Figure 5 indicate that increasing variability in the system (either for service times or for the arrival process) shifts the probability from states with smaller number in the system to the states with higher number in the system, thus creating a longer tail for the distribution. Model 2, in general, has a slightly longer tail relative to Model 1.

5.3. Impact of γ

γ is not under the control of the management but understanding its effect on the system is essential to choose K and L efficient operation of the system. We considered values of γ =1/10, 1/15, 1/20, and 1/25 (with λ = 1 ) in the numerical analysis.
For a given K and L, larger values of γ resulted in more frequent usage of opportunistic events to procure replenishment or smaller order cycle times, κ 1 and κ 2 ), and hence smaller size replenishments, Γ 1 and Γ 2 , resulting in a net increase in the mean inventory levels μ ^ 1 and μ ^ 2 ). This observation should be interpreted with caution because it states the effect of larger γ , keeping K and L at the same level. Larger γ indicates more frequent opportunities for replenishment permitting smaller K and L to achieve the same service level. Larger γ , in addition to increasing the means, μ ^ 1 and μ ^ 2 , also increased the standard deviations, σ ^ 1 and σ ^ 2 , but with a net decrease in the coefficients of variation of the inventory levels. Larger γ also had a similar effect on the numbers of customers in the system by increasing the means, μ 1 and μ 2 , and the standard deviations, σ 1 and σ 2 , but resulted in a net decrease in the coefficients of variation. Further discussion on the combined effect of changes in γ with other parameters such as K and L is presented in the following.
Inventory system with larger values of γ , indicating frequent replenishment opportunities, will permit maintaining smaller inventories to achieve the same level of service.

5.4. Impact of K and L

In the inventory subsystem, as K increases, with all other parameters remaining the same, more items are procured when opportunistic events occur leading to larger inventories on average. Larger inventories lead to larger cycle times ( κ 1 and κ 2 ) because of an increased likelihood of skipping current opportunities to wait for future ones. The net effect is an increase in the means μ ^ 1 and μ ^ 2 and also an increase in the standard deviations σ ^ 1 and σ ^ 2 , but a reduction in the coefficient of variation of the inventory levels. Changes in L had a slightly different, and smaller, effect on the inventory subsystem than corresponding changes in K. As L increases with all other parameters remaining the same, opportunistic events are availed more often (smaller κ 1 and κ 2 ) and in smaller quantities. The net result is once again an increase in the means, μ ^ 1 and μ ^ 2 , and an increase in the standard deviations, σ ^ 1 and σ ^ 2 ), but a reduction in the coefficient of variation, of the numbers in the system.
In the customer subsystem, increasing K or L keeping the rest of the parameters the same, led to an increase in the means, μ 1 and μ 2 , and an increase in the standard deviations, σ 1 and σ 2 , but a reduction in the coefficient of variation, of the number of customers in the system. Increasing K or L also resulted in a decrease the customer loss probability ( ϑ 1 , ϑ 2 ) , server idle probability ( ν 1 , ν 2 ), and an increase in the percent of server idle time with positive inventory ( ν 1 * , ν 2 * ). Effect of changes in K are much stronger than corresponding changes in L.

5.5. Cost Analysis

A system that can provide the desired service level at the minimum cost is the goal of any optimization. We consider three costs associated with the system, namely, cost of carrying inventory, cost of placing an order, and cost of lost customers/demand. A fourth item, the cost of keeping customers waiting, is not considered in this study because it was not considered particularly relevant. Let c 1 , c 2 , and c 3 denote the costs of carrying a unit inventory per unit time, cost of placing an order, and the average cost of a lost customer. Z 1 and Z 2 , the total cost per unit time for the two models can be expressed as follows. The system cost per unit time is then given by
Z 1 = c 1 μ ^ 1 + c 2 κ 1 + c 3 ϑ 1 Z 2 = c 1 μ ^ 2 + c 2 κ 2 + c 3 ( ϑ 2 a + ϑ 2 d )
The first two components in Z 1 and Z 2 represent the cost of the inventory subsystem, and the third item represents the cost of the customer subsystem. The goal of management is to choose the values of K and L that minimize the total cost. Using values of c 1 = 0.25, c 2 = 100, and c 3 = 10, a summary of the cost for a select set of parameter values is presented in Table 4, where the minimum cost values for each combination of K and L are highlighted in bold. In order to find the best combination of values of K and L, a more refined search than in Table 4 needs to be performed.

5.6. Comparing Model 1 and Model 2

Table 3 compares the steady behavior of the system under the two models. For a given set of parameter values, Model 1 performed better than Model 2 in terms of the mean number of customers in the system, mean inventory level, and probability of an idle server. For Model 2, the loss of demand can occur either at arrival, or at service completion epochs. ϑ 2 is always larger than ϑ 1 , but the difference is significantly higher for the case of HES. This high loss rate may cause higher idle probability. Higher variability of the arrival process (as in the case of HEA or PSA) appears to mitigate the effect. For Model 2, the proportion of customers lost at departure instants varied from a low of 9.01% to a high of 68.4%. We notice that these percentages are significantly larger for HES, HEA, and PCA.

5.7. Comparing ( K , L ) System with ( s , S ) System

In this subsection, we compare the characteristics of the system considered in this paper (the ( K , L ) system) with the ( s , S ) system considered in Chakravarthy and Rumyantsev [17]. In the ( K , L ) system, the replenishment opportunities occur randomly at exponentially distributed intervals with no lead time for delivery. In the ( s , S ) system in [17], orders are placed as inventory is depleted but the items are received after an exponentially distributed lead time. In other words, the random lead time in the ( s , S ) system is replaced by the random interval between replenishment opportunities in the ( K L ) system. For effective comparison, we use exactly the same parameters for both systems and set K = S and L = s , and the mean lead time for delivery in the ( s , S ) system equal to the mean time between two replenishment opportunities in the ( K , L ) system. Table 5 presents a summary comparison for γ = 0.1, μ = 1.1, K = S = 60 , L = s = 20 for both Model 1 and Model 2, for all arrival and service distributions considered in Section 3.
The summary indicates that for the combinations of parameters considered, ( K , L ) system has a smaller server idle probability, smaller customer loss probability but on all other counts ( s , S ) system had better operational performance for both Model 1 and Model 2. Considering the customer subsystem, ( K , L ) model has larger mean and standard deviation of the number of customers in the system relative to the ( s , S ) system. In the inventory subsystem, ( K , L ) model has a larger cycle time, and a larger mean and standard deviation of the inventory level relative to the ( s , S ) system. This indicates that orders are placed less often in the ( K , L ) system, and more units are ordered each time an order is placed, resulting in larger mean inventory level with greater variability. The overall conclusion from these observations is that the additional uncertainty introduced in the system by the random supply process adversely affects the system performance.
Cost computations are not performed because the results depend significantly on the relative values of c 1 , c 2 , and c 3 . Furthermore, direct comparison of total costs could be misleading because one has to optimize the two systems separately and then compare the two optimal solutions to see which performs better. Even that is not very meaningful, because one does not choose between a ( K , L ) system or an ( s , S ) system.

6. Concluding Remarks

In this paper, we considered a queuing-inventory system where the replenishment opportunities occur at random intervals at which time the decisions are made with regard to the procurement and the quantities. We proposed an opportunistic system involving two parameters, L and K, similar to the well-known ( s , S ) inventory management system. The cut-off point, L, between 0 and K, largely determines whether to avail the opportunistic replenishment with certainty or with a probabilistic rule. The two models considered differ in the way the customers enter the system or are getting lost. In Model 1, the customers are lost when the inventory level is zero (irrespective of whether the server is idle or busy). In Model 2, the customers are allowed to enter even when the inventory level is zero but as long as the server is busy serving. Further, the customers may be lost at a service completion due to lack of inventory to start a new service.
The analysis of the results from an extensive set of numerical experimentation under various scenarios of the parameters indicates that with all other parameters remaining the same, Model 1 appears to have better performance characteristics in terms of the number of customers in the system, inventory levels in the system, and the likelihood of losing a customer.
Variability in the overall system behavior can be due to the variability in the arrival process or the service time distribution. In general, increasing the variability of either the input process or the service time distribution, led to an increase in overall variability in the system behavior. When both arrival process and service time distribution are highly variable, the combined effect is more attenuated than the individual effects.
With the tools developed in this paper and the understanding of the system behavior presented in this paper, it is possible to conduct an efficient search for the value of the key parameters K, and L that minimizes the total system cost. While γ , the rate at which replenishment opportunities arise, is not under the control of the management, one might surmise that an optimal inventory system with larger values of γ will have smaller cost because of the availability of frequent replenishment opportunities making it possible to achieve the same level of service with smaller inventories.
A comparison with the ( s , S ) system studied by Chakravarthy and Rumyantsev [17] revealed that, the ( s , S ) system has better operating characteristics in terms of number in the system and inventory levels, apparently due to the random supply process. This could be mitigated by the possibly having a lower purchase cost per unit. This aspect is not studied in this paper. In addition, it would be interesting to compare the system studied here with the ( s , S ) system in which the inventory level is always brought to level S, irrespective of when the replenishment occurs.
A possible extension to this study could be to investigate systems with both regular channels of purchase (as in the ( s , S ) system) and the opportunistic replenishment opportunities (as in the ( K , L ) system) with possibly lower per unit cost of purchase in the context of M A P arrivals of demands and phase-type services.
In this paper, we explored two queuing-inventory models where the replenishment opportunities occur at random. Such models are appropriate for non-critical items (e.g., canned goods, generic medication) where stock outs for short periods of time are not serious. Special replenishment opportunities typically offer lower unit cost with little or no ordering cost, and units are usually available without delay. The results of this paper and the detailed discussion of the system behavior, offer a useful tool to compare the traditional inventory systems (e.g., ( s , S ) system) with the opportunistic models to determine the most effective method to manage inventories under a more general scenario of using M A P for the demand process and phase-type distributions for the processing of the demands.

Author Contributions

Conceptualization: S.R.C. and B.M.R.; methodology: S.R.C.; software: S.R.C.; validation: S.R.C. and B.M.R.; Formal analysis: S.R.C. and B.M.R.; investigation: S.R.C. and B.M.R.; writing—original draft preparation, S.R.C. and B.M.R.; writing—review and editing, S.R.C. and B.M.R.; visualization, S.R.C. and B.M.R. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Data Availability Statement

Not applicable.

Acknowledgments

The authors gratefully acknowledge the many suggestions by the anonymous reviewers which significantly improved the presentation of the paper.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Friend, J.K.; Kaminsky, F.C. Stock Control with Random Opportunities for Replenishment. J. Oper. Res. Soc. 1960, 11, 130–136. [Google Scholar] [CrossRef]
  2. Hurter, A.P.; Kaminsky, F.C. Inventory Control with a Randomly Available Discount Purchase Price. Oper. Res. Q. 1968, 19, 433–444. [Google Scholar] [CrossRef]
  3. Silver, E.A.; Robb, D.J.; Rahnama, M.R. Random Opportunities for Reduced Cost Replenishments. IIE Trans. 1993, 25, 111–120. [Google Scholar] [CrossRef]
  4. Gurnani, H. Optimal ordering policies in inventory systems with random demand and random deal offerings. Eur. J. Oper. Res. 1996, 95, 299–312. [Google Scholar] [CrossRef]
  5. Moinzadeh, K. Replenishment and Stocking Policies for Inventory Systems with Random Deal Offerings. Manag. Sci. 1997, 43, 334–342. [Google Scholar] [CrossRef]
  6. Feng, Y.; Sun, J. Computing the optimal replenishment policy for inventory systems with random discount opportunities. Oper. Res. 2001, 49, 790–795. [Google Scholar] [CrossRef]
  7. Goh, M.; Sharafali, M. Price-dependent inventory models with discount offers at random times. Prod. Oper. Manag. 2002, 11, 139–156. [Google Scholar] [CrossRef]
  8. Chaouch, B.A. Inventory control and periodic price discounting campaigns. Nav. Res. Logist. 2007, 54, 94–108. [Google Scholar] [CrossRef]
  9. Tajbakhsh, M.M.; Lee, C.-G.; Zolfaghari, S. An inventory model with random discount offerings. Omega 2011, 39, 710–718. [Google Scholar] [CrossRef]
  10. Karimi-Nasab, M.; Konstantaras, I. An inventory control model with stochastic review interval and special sale offer. Eur. J. Oper. Res. 2013, 227, 81–87. [Google Scholar] [CrossRef]
  11. Den Boer, A.; Perry, O.; Zwart, B. Dynamic pricing policies for an inventory model with random windows of opportunities. Nav. Res. Logist. 2018, 65, 660–675. [Google Scholar] [CrossRef] [Green Version]
  12. Krishnamoorthy, A.; Shajin, D.; Narayanan, V. Inventory with positive service time: A survey. In Advanced Trends in Queueing Theory; Mathematics and Statistics Series; ISTE &amp Wiley: London, UK, 2019. [Google Scholar]
  13. Choi, K.H.; Yoon, B.K.; Moon, S.A. Queueing Inventory Systems with Phase-type Service Distributions: A Literature Review. Ind. Eng. Manag. Syst. 2019, 18, 330–339. [Google Scholar] [CrossRef]
  14. Melikova, A. Z.; Shahmaliyeva, M. O. Markov Models of Inventory Management Systems with a Positive Service Time. J. Comput. Syst. Sci. Int. 2018, 57, 766–783. [Google Scholar] [CrossRef]
  15. Chakravarthy, S.R.; Maity, A.; Gupta, U.C. Modeling and Analysis of Bulk Service Queues with an Inventory under (s,S) Policy. Ann. Operat. Res. 2017, 258, 263–283. [Google Scholar] [CrossRef]
  16. Chakravarthy, S.R. Queueing-inventory models with batch demands and positive service times. Math. Game Theory Appl. 2019, 11, 95–120. [Google Scholar]
  17. Chakravarthy, S.R.; Rumyantsev, A. Analytical and simulation studies of queueing-inventory models with MAP demands in batches and positive phase type services. Simul. Model. Pract. Theory 2020, 103, 102092. [Google Scholar] [CrossRef]
  18. Alfa, A.S. Queueing Theory for Telecommunications; Springer Science+Business Media, LLC: Berlin, Germany, 2010. [Google Scholar]
  19. Artalejo, J.R.; Gomez-Correl, A.; He, Q.M. Markovian arrivals in stochastic modelling: A survey and some new results. SORT 2010, 34, 101–144. [Google Scholar]
  20. Chakravarthy, S.R. The batch Markovian arrival process: A review and future work. In Advances in Probability Theory and Stochastic Processes; Krishnamoorthy, A., Ed.; Notable Publications Inc.: Trenton, NJ, USA, 2001; pp. 21–39. [Google Scholar]
  21. Chakravarthy, S.R. Markovian Arrival Processes. In Wiley Encyclopedia of Operations Research and Management Science; John Wiley & Sons: New York, NY, USA, 2010. [Google Scholar]
  22. Chakravarthy, S.R. Matrix-Analytic Queueing Models. In An Introduction to Queueing Theory, 2nd ed.; Narayan Bhat, U., Ed.; Birkhauser; Springer Science + Business Media: New York, NY, USA, 2015; Chapter 8. [Google Scholar]
  23. He, Q.-M. Fundamentals of Matrix-Analytic Methods; Springer: New York, NY, USA, 2014. [Google Scholar]
  24. Lucantoni, D.M.; Meier-Hellstern, K.S.; Neuts, M.F. A single-server queue with server vacations and a class of nonrenewal arrival processes. Adv. Appl. Probl. 1990, 22, 676–705. [Google Scholar] [CrossRef]
  25. Lucantoni, D.M. New results on the single server queue with a batch Markovian arrival process. Stoch. Model. 1991, 7, 1–46. [Google Scholar] [CrossRef]
  26. Neuts, M.F. A versatile Markovian point process. J. Appl. Prob. 1979, 16, 764–779. [Google Scholar] [CrossRef]
  27. Neuts, M.F. Structured Stochastic Matrices of M/G/1 Type and Their Applications; Marcel Dekker, Inc.: New York, NY, USA, 1989. [Google Scholar]
  28. Neuts, M.F. Models based on the Markovian arrival process. IEICE Trans. Commun. 1992, E75B, 1255–1265. [Google Scholar]
  29. Neuts, M.F. Algorithmic Probability: A Collection of Problems; Chapman and Hall: New York, NY, USA, 1995. [Google Scholar]
  30. Neuts, M.F. Matrix-Geometric Solutions in Stochastic Models: An Algorithmic Approach; The Johns Hopkins University Press: Baltimore, MD, USA, 1981. [Google Scholar]
  31. Neuts, M.F. Probability distributions of phase type. In Liber Amicorum Prof. Emeritus H. Florin; Department of Mathematics, University of Louvain: Ottignies-Louvain-la-Neuve, Belgium, 1975; pp. 173–206. [Google Scholar]
  32. Marcus, M.; Minc, H. A Survey of Matrix Theory and Matrix Inequalities; Courier Corporation: Chelmsford, MA, USA, 1992. [Google Scholar]
  33. Steeb, W.-H.; Hardy, Y. Matrix Calculus and Kronecker Product: A Practical Approach to Linear and Multilinear Algebra, 2nd ed.; World Scientific Publishing Company: Singapore, 2011. [Google Scholar]
  34. Latouche, G.; Ramaswami, V. Introduction to Matrix Analytic Methods in Stochastic Modeling; ASA-SIAM: Philadelphia, PA, USA, 1999. [Google Scholar]
Figure 1. Probability density functions of the inter-arrival times.
Figure 1. Probability density functions of the inter-arrival times.
Mathematics 09 01092 g001
Figure 2. Joint Probability Density Function of M A P P C .
Figure 2. Joint Probability Density Function of M A P P C .
Mathematics 09 01092 g002
Figure 3. Probability density function of the service times.
Figure 3. Probability density function of the service times.
Mathematics 09 01092 g003
Figure 4. Marginal distributions of number of customers in the system (Model 1).
Figure 4. Marginal distributions of number of customers in the system (Model 1).
Mathematics 09 01092 g004
Figure 5. Marginal distributions of number of customers in the system (Model 2).
Figure 5. Marginal distributions of number of customers in the system (Model 2).
Mathematics 09 01092 g005
Table 1. Customer subsystem ( μ = 1.1).
Table 1. Customer subsystem ( μ = 1.1).
Model 1Model 2
Arr.Ser. γ K L μ 1 σ 1 ν 1 ν 1 * ϑ 1 μ 2 σ 2 ν 2 ν 2 * ϑ 2 a ϑ 2 d
ERAERS0.0550200.8001.2430.5940.1460.5530.8511.3870.6050.1400.5200.046
ERAERS0.0550300.8411.2840.5810.1510.5390.8961.4330.5910.1450.5050.045
ERAERS0.0560200.9291.3620.5520.1620.5070.9781.5000.5640.1560.4760.044
ERAERS0.0560300.9821.4120.5360.1690.4901.0331.5540.5480.1610.4600.043
ERAERS0.150201.3071.5360.3940.2950.3331.3561.6830.4140.2750.3000.056
ERAERS0.150301.4121.6180.3680.3140.3051.4651.7700.3880.2920.2750.052
ERAERS0.160201.4671.6570.3550.3250.2901.5041.7890.3770.3020.2630.051
ERAERS0.160301.5901.7540.3280.3480.2611.6321.8890.3490.3230.2360.048
ERAHES0.0550203.5987.5540.5940.3010.5544.49310.9160.6460.2820.4610.149
ERAHES0.0550303.8998.0400.5810.3100.5394.82711.4570.6330.2910.4470.150
ERAHES0.0560204.5539.0000.5510.3350.5065.16011.8270.6130.3110.4200.154
ERAHES0.0560304.9769.6480.5360.3460.4905.58012.4760.5980.3210.4040.154
ERAHES0.150208.94214.6320.3940.5220.3336.78513.3280.5040.4740.2630.192
ERAHES0.1503010.23316.2460.3680.5490.3057.42814.1500.4830.4940.2420.189
ERAHES0.1602010.93517.0850.3550.5650.2907.59014.2880.4740.5090.2310.191
ERAHES0.1603012.65019.1470.3280.5960.2618.36815.2600.4520.5320.2100.187
HEAERS0.0550201.8543.7550.7020.5910.6724.1039.8680.7290.5270.3680.334
HEAERS0.0550301.9403.8660.6950.5970.6644.41110.3660.7170.5340.3550.334
HEAERS0.0560202.4334.6330.6640.6180.6304.89811.0880.6990.5440.3370.332
HEAERS0.0560302.5604.7900.6550.6270.6205.30511.7150.6850.5520.3230.330
HEAERS0.150203.1405.3780.5890.8110.5486.31012.7240.6300.7070.2020.391
HEAERS0.150303.3635.6370.5760.8210.5337.10413.8250.6060.7180.1860.381
HEAERS0.160204.1356.7370.5460.8360.5007.41614.2060.5970.7230.1790.378
HEAERS0.160304.4697.1210.5300.8480.4838.44115.5770.5700.7360.1620.365
HEAHES0.0550202.7246.1660.7020.5840.6726.33117.0110.7380.5300.3960.315
HEAHES0.0550302.8636.3850.6950.5900.6646.75617.7560.7280.5360.3850.316
HEAHES0.0560203.6087.6590.6640.6110.6307.35118.5280.7090.5490.3610.319
HEAHES0.0560303.8187.9750.6550.6190.6207.89719.4370.6970.5560.3480.318
HEAHES0.150205.20810.0940.5890.7810.5488.91220.3010.6540.7020.2400.379
HEAHES0.150305.61010.6560.5760.7900.5339.73521.5060.6360.7100.2260.374
HEAHES0.160206.86312.5350.5460.8060.50110.19221.9850.6230.7190.2120.374
HEAHES0.160307.46813.3480.5300.8170.48311.22623.4410.6040.7290.1960.368
PCAERS0.0550200.7622.3370.6730.4730.6402.37212.7370.6820.4620.4190.230
PCAERS0.0550300.7982.4070.6600.4980.6262.64113.8430.6690.4860.3960.239
PCAERS0.0560200.9783.0940.6450.5170.6103.02015.3350.6580.5030.3750.248
PCAERS0.0560301.0283.1950.6320.5460.5953.40816.8160.6430.5310.3490.257
PCAERS0.150201.5534.6130.5530.7270.5084.36219.5150.5740.7010.2180.311
PCAERS0.150301.6544.8250.5370.7680.4905.18622.2250.5550.7400.1880.320
PCAERS0.160202.0926.2370.5290.7650.4815.42723.0060.5550.7360.1870.321
PCAERS0.160302.2386.5430.5130.8060.4646.57326.4740.5360.7740.1580.328
PCAHES0.0550202.1325.2850.6730.5250.6405.14222.4480.7010.4960.4320.239
PCAHES0.0550302.2905.5520.6600.5490.6265.57423.6850.6880.5160.4130.244
PCAHES0.0560202.6166.3160.6450.5680.6106.03024.9150.6780.5340.3900.256
PCAHES0.0560302.8116.6340.6320.5940.5956.58926.4490.6640.5560.3690.261
PCAHES0.150204.2489.2390.5530.7620.5087.71228.4190.6070.7120.2540.314
PCAHES0.150304.6069.7770.5370.7970.4908.54530.4790.5900.7420.2300.319
PCAHES0.160205.19511.1930.5290.7960.4818.87031.2870.5880.7420.2220.325
PCAHES0.160305.62511.8380.5130.8300.4649.93033.8250.5710.7720.1990.329
Table 2. Inventory subsystem— μ = 1.1, ERS.
Table 2. Inventory subsystem— μ = 1.1, ERS.
Model 1Model 2
Arr. γ K L μ ^ 1 σ ^ 1 Γ 1 ξ 1 κ 1 μ ^ 2 σ ^ 2 Γ 2 ξ 2 κ 2
ERA0.05502011.69216.32147.0550.73527.22812.39816.64246.8970.71827.861
ERA0.05503012.41016.72744.5360.80424.87813.15517.03444.2600.78925.342
ERA0.05602015.42119.91856.7460.67729.53716.25320.24756.5780.66030.296
ERA0.05603016.40320.33354.2400.73527.21517.22320.63254.0030.72027.791
ERA0.10502018.14217.41644.6100.58317.15818.84917.48644.3790.56617.678
ERA0.10503019.86417.75240.3560.67314.85520.58317.77840.0150.65815.206
ERA0.10602022.99520.64454.0130.51519.43323.72720.71553.7880.50020.003
ERA0.10603025.10120.82049.8220.58317.15525.83920.81849.5090.56917.592
HEA0.05502022.05922.42347.8420.52937.78123.09621.21346.3140.49940.081
HEA0.05503022.78622.72646.1490.56335.50223.94321.33143.3730.55735.894
HEA0.05602026.85526.41357.5320.49940.04928.05824.90855.9050.46243.287
HEA0.05603027.82026.71455.8480.52937.78229.09024.93253.0560.51139.182
HEA0.10502028.90921.53946.5400.37626.57228.91019.60343.7570.36327.541
HEA0.10503030.17821.72044.1340.41024.36830.31119.38339.3220.43123.199
HEA0.10602034.64925.15055.9540.34828.75234.48422.87352.9920.32830.487
HEA0.10603036.26025.21953.5750.37626.57336.08022.41948.7080.38226.208
PCA0.05502016.74717.68445.5210.61532.50817.32217.64445.1900.60333.193
PCA0.05503018.18418.11741.8380.69728.70118.78818.02241.3350.68829.058
PCA0.05602021.26921.10355.1110.55336.14921.97921.00054.7030.53937.141
PCA0.05603023.03021.43651.5620.61532.51023.77821.25550.9430.60433.102
PCA0.10502023.34617.18842.4580.45422.05523.96516.85541.8080.44022.751
PCA0.10503025.98017.20636.9850.54118.48026.61816.73636.1010.53418.762
PCA0.10602028.54720.16751.7330.39425.40829.26219.73050.9260.37926.462
PCA0.10603031.49219.90546.5400.45422.05532.21119.31145.4460.44422.564
Table 3. Comparing Model 1 and Model 2
Table 3. Comparing Model 1 and Model 2
Model 1 ( μ = 1.1, γ = 0.1, K = 50)
Arr.Ser. L μ 1 σ 1 μ ^ 1 σ ^ 1 Γ 1 κ 1 ν 1 ν 1 * ϑ 1 ϑ 1 *
ERAERS201.3071.53618.14217.41644.61017.1580.3940.2950.3330.000
ERAERS301.4121.61819.86417.75240.35614.8550.3680.3140.3050.000
ERAEXS201.8442.39118.14317.41544.60917.1580.3940.3450.3330.000
ERAEXS302.0232.56619.86717.75140.35314.8550.3680.3660.3050.000
ERAHES208.94214.63218.14117.41544.60817.1650.3940.5220.3330.000
ERAHES3010.23316.24619.86517.75040.34814.8610.3680.5490.3050.000
HEAERS203.1405.37828.90921.53946.54026.5720.5890.8110.5480.000
HEAERS303.3635.63730.17821.72044.13424.3680.5760.8210.5330.000
HEAEXS203.2445.60528.90921.53946.54026.5700.5890.8090.5480.000
HEAEXS303.4765.87830.17821.72044.13424.3660.5760.8190.5330.000
HEAHES205.20810.09428.90921.53946.54026.5730.5890.7810.5480.000
HEAHES305.61010.65630.17721.72044.13424.3670.5760.7900.5330.000
PCAERS201.5534.61323.34617.18842.45822.0550.5530.7270.5080.000
PCAERS301.6544.82525.98017.20636.98518.4800.5370.7680.4900.000
PCAEXS201.6814.77123.34517.19042.46022.0500.5530.7320.5080.000
PCAEXS301.7934.99125.98217.20536.98318.4790.5370.7720.4900.000
PCAHES204.2489.23923.34517.18942.45822.0550.5530.7620.5080.000
PCAHES304.6069.77725.97817.20736.98718.4800.5370.7970.4900.000
ERAERS201.3561.68318.84917.48644.37917.6780.4140.2750.3560.156
ERAERS301.4651.77020.58317.77840.01515.2060.3880.2920.3270.160
ERAEXS201.9162.67819.33417.62544.31218.0460.4280.3140.3710.210
ERAEXS302.1012.86021.11217.92539.93215.5090.4020.3320.3420.217
ERAHES206.78513.32822.40818.54144.22820.7770.5040.4740.4540.422
ERAHES307.42814.15024.55119.02040.07118.0150.4830.4940.4310.439
HEAERS206.31012.72428.91019.60343.75727.5410.6300.7070.5930.659
HEAERS307.10413.82530.31119.38339.32223.1990.6060.7180.5670.673
HEAEXS206.48013.19029.02519.63543.82627.7330.6320.7060.5950.654
HEAEXS307.28014.30630.44619.43739.45723.4270.6090.7170.5700.668
HEAHES208.91220.30130.27919.90744.55030.0390.6540.7020.6190.612
HEAHES309.73521.50631.91219.89640.90226.2000.6360.7100.6000.624
PCAERS204.36219.51523.96516.85541.80822.7510.5740.7010.5290.587
PCAERS305.18622.22526.61816.73636.10118.7620.5550.7400.5090.630
PCAEXS204.56820.08724.06716.89141.82322.8560.5760.7020.5330.584
PCAEXS305.39222.75026.71916.77436.13418.8510.5570.7400.5120.625
PCAHES207.71228.41925.60717.46042.14824.9120.6070.7120.5680.552
PCAHES308.54530.47928.30717.48936.85820.8220.5900.7420.5490.581
Table 4. Total cost summary γ = 0.1, μ = 1.1.
Table 4. Total cost summary γ = 0.1, μ = 1.1.
Service Dist.ERSEXSHES
Arr. L / K 405060405060405060
Model 1
ERA1013.08512.94613.12513.08412.94213.12313.08012.94113.123
ERA2013.98213.69513.79913.98113.69513.79713.97913.69113.796
ERA3015.30414.75114.71815.30314.75014.71715.30114.74514.708
HEA1015.45315.99916.66715.45316.00016.66915.45316.00016.668
HEA2015.92516.47017.14515.92616.47117.14615.92616.47017.146
HEA3016.44416.98217.66216.44516.98317.66316.44416.98317.663
PCA1014.16114.44514.97614.16414.44614.98114.16214.44414.979
PCA2015.35515.44815.88615.35615.45115.88815.35515.44915.885
PCA3017.12516.80717.04617.12716.81017.05017.12516.80817.047
Model 2
ERA1013.25813.18413.41213.36313.33313.59813.91314.12014.596
ERA2014.16913.92814.07714.28114.08414.26014.87114.95715.364
ERA3015.54114.99314.97915.65615.14415.16116.15816.00016.289
HEA1015.62216.26217.03415.63516.28317.06215.77816.50717.358
HEA2016.28016.78617.46516.29216.81217.50216.42317.08917.901
HEA3017.26317.55718.10617.25817.57718.14317.22917.79518.550
PCA1014.30614.66415.27714.33314.70315.32314.61815.09615.803
PCA2015.53315.68116.17415.55615.72016.22815.78716.09316.710
PCA3017.36017.07217.34717.37517.10717.40117.42817.37317.817
Table 5. Comparing the ( K , L ) system with the ( s , S ) system.
Table 5. Comparing the ( K , L ) system with the ( s , S ) system.
( K , L ) System ( γ = 0.1, μ = 1.1, K = 60, L = 20) ( s , S ) System ( γ = 0.1, μ = 1.1, S = 60, s = 20)
Model 1
Arr.Ser. ϑ 1 ν 1 μ 1 σ 1 μ ^ 1 σ ^ 1 κ 1 ϑ 1 ν 1 μ 2 σ 2 μ 2 σ ^ 2 κ 1
ERAERS0.3550.2901.4671.65722.99520.64419.4330.4230.3651.2091.47515.66016.17916.221
ERAEXS0.3550.2902.1162.64422.99920.64219.4390.4230.3651.6862.26715.66316.17816.222
ERAHES0.3550.29010.93517.08522.98620.64619.4370.4230.3657.78413.21415.65716.17916.223
HEAERS0.5460.5004.1356.73734.64925.15028.7520.6290.5922.4384.41424.48119.09825.437
HEAEXS0.5460.5014.2727.01434.65025.14928.7490.6290.5922.5204.60124.48119.09825.436
HEAHES0.5460.5016.86312.53534.64925.15028.7510.6290.5924.0188.25924.48119.09825.437
PCAERS0.5290.4812.0926.23728.54720.16725.4080.5680.5251.2203.46821.74816.97221.571
PCAEXS0.5290.4822.2396.40928.55020.16325.4110.5680.5251.3383.62021.74916.97021.571
PCAHES0.5290.4815.19511.19328.54720.16325.4150.5680.5253.6777.96821.74716.97221.572
Model 2
Arr.Ser. ϑ 2 ν 2 μ 2 σ 2 μ ^ 2 σ ^ 2 κ 2 ϑ 2 ν 2 μ 2 σ 2 μ 2 σ ^ 2 κ 2
ERAERS0.3770.3151.5041.78923.72720.71520.0030.4410.3851.2711.63716.35716.73016.748
ERAEXS0.3910.3302.1512.88324.27320.83820.4450.4530.3981.7962.60016.84316.94617.113
ERAHES0.4740.4217.59014.28827.67321.70823.6190.5240.4766.41812.95919.84618.27019.642
HEAERS0.5970.5567.41614.20634.48422.87330.4870.6510.6165.76612.09925.21818.10926.769
HEAEXS0.5990.5597.59414.67434.62322.90930.6900.6530.6185.93112.56725.30718.15326.928
HEAHES0.6230.58610.19221.98536.11823.21033.1740.6750.6438.19919.48826.32618.62928.909
PCAERS0.5550.5085.42723.00629.26219.73026.4620.5830.5404.02418.49422.36317.22622.283
PCAEXS0.5570.5125.64123.52729.38019.76826.5650.5850.5434.22019.05822.43617.20722.384
PCAHES0.5880.5478.87031.28731.06920.38728.7990.6170.5787.24627.33223.75517.80024.279
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Chakravarthy, S.R.; Rao, B.M. Queuing-Inventory Models with MAP Demands and Random Replenishment Opportunities. Mathematics 2021, 9, 1092. https://doi.org/10.3390/math9101092

AMA Style

Chakravarthy SR, Rao BM. Queuing-Inventory Models with MAP Demands and Random Replenishment Opportunities. Mathematics. 2021; 9(10):1092. https://doi.org/10.3390/math9101092

Chicago/Turabian Style

Chakravarthy, Srinivas R., and B. Madhu Rao. 2021. "Queuing-Inventory Models with MAP Demands and Random Replenishment Opportunities" Mathematics 9, no. 10: 1092. https://doi.org/10.3390/math9101092

APA Style

Chakravarthy, S. R., & Rao, B. M. (2021). Queuing-Inventory Models with MAP Demands and Random Replenishment Opportunities. Mathematics, 9(10), 1092. https://doi.org/10.3390/math9101092

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop