A joint replenishment production-inventory model as an MMAP[K]/PH[K]/1 queue

Abstract In this paper we analyze a continuous review finite capacity production-inventory system with two products in inventory. With stochastic order quantities and time between orders, the model reflects a supply chain that operates in an environment with high levels of volatility. The inventory is replenished using an independent order-up-to (s, S) policy or a can-order (s, c, S) joint replenishment policy in which the endogenously determined lead times drive the parameters of the replenishment policy. The production facility is modeled as a multi-type MMAP[K]/PH[K]/1 queue in which there are K possible inventory positions when the order is placed and the age process of the busy queue has matrix-exponential distribution. We characterize the system and determine the steady state distribution using matrix analytic methods. Using numerical methods we obtain the inventory parameters that minimize the total ordering and inventory related costs. We present numerical comparisons of independent and joint replenishment policies with varying lead times, order quantities, and cost reductions. We further demonstrate the interplay between the two products in terms of lead times, order quantities and costs.


Introduction
Inventory management has been studied for decades, but volatile market conditions have increased the complexity of modelling and analysis, as supply chain factors such as production rate, demand rate, and lead time become more variable.While much of the research in inventory management deals with single product analysis, companies need to make decisions on many aspects of inventory management, and the replenishment of multiple items is a common operational concern [6] .
We study a continuous review production-inventory system with a single production line, two products in inventory, stochastic order quantities, stochastic time between orders, and where the objective is to determine the parameters of the inventory policy for each product that minimizes the total ordering and inventory related costs for both products.
We consider two policies of replenishing each product's inventory: an order-up-to (s, S) policy in which each product is replenished independently (Arrow et al. [1] ); and a can-order joint replenishment (s, c, S) policy (Balintfy [5] ), where if the inventory position of one product reaches its reorder point s or below, all other products whose inventory position is below their can-order level c, are included in the same replenishment order and any order placed raises the inventory position of that product to the order-up-to level S. Considerable savings may be realised by such coordinated replenishments (Bastos et al. [8] ) and an important factor making this approach cost effective is the joint set-up cost structure.However, the interaction between items makes determination of optimal parameters a difficult task.
The lead time is the time between the moment an order is placed and the moment it is received in inventory and thus consists of a waiting time in queue (if the system is busy), a setup time and a production time.The lead time depends on the way orders are placed and the production process, and influences inventory levels and inventory related costs.
Much of the inventory literature treats lead times as a fixed constant or as an exogenous variable with a given probability distribution, which means that the time required to deliver an order is assumed to be independent of the size of the current order and independent of the lead time of previous orders (Kulkarni and Yan [22] , Zipkin [40] ).This approach is justified when both production and inventory are decoupled through a large inventory at the production; if the owner of the production system guarantees a fixed delivery date; or if transportation lead times are much longer than production lead times, as in these scenarios the inventory policy does not impact significantly on lead times.
When multiple products are produced on the same production line, the order process of each individual product loading the production queue, impacts the lead times for all products ordered and all inventory levels are affected.Ignoring this impact may lead to suboptimal inventory parameters and higher costs.For example, in the context of a fixed order period (r, S) policy, it was shown that relying on exogenous lead times can result in a severe underestimation of the required inventory on hand (Van Houdt and Pérez [39] ).
In an integrated production-inventory setting, the interaction between the inventory control system and the production system is explicitly modelled: inventory influences production by initiating orders, and production influences inventory by completing and delivering orders to inventory.Accordingly, the lead times are endogenously determined by the production system: as the orders of each product load the production line, they determine the lead times, hence influence the time to replenish each product in inventory, and impact the holding and backlogging costs.Thus there is correlation between demand and lead times: high demand depletes inventory more, creating a larger order to replenish inventory, which in turn takes longer to produce.In order to characterize the retailer's inventory behavior in a production-inventory system in an exact way, we include this correlation in our analysis.
Much queueing theory literature assumes the arrival processes of customers are independent processes.The introduction of the Markov arrival process with marked transitions, MMAP[K] (He [17] ), allowed modelling of complicated queues with multiple types of customers and correlated arrival times (Artalejo [2] ).Further, the application of standard matrix analytic methods has allowed the development of efficient algorithms for computing relevant performance measures (for example, He [15] and Horvath [18] ).
In this paper, the production facility is modelled as a MMAP[K]/PH[K]/1 queue (He [16] , Van Houdt and Blondia [38] ).Queues of this type have a single server, correlated interarrival times and K types of customers, where each customer type may have a different phase-type service requirement.In our case, the customer types correspond to the possible inventory positions when the order is placed.
The state space for the process records information on the type of order, the inventory position when the order was placed, and the current state of production or setup.The age process that observes this queue during the busy period keeps track of the age of the order in production, the current production phase and has matrix-exponential distribution (Sengupta [31,32] ).To enable efficient computation of the inventory level distribution, we construct a fluid queue to enable computation of required matrices (Dziel et al. [13] ) and from these we can determine the costminimizing order quantity of the replenishment orders.
In our model the setup, change-over and production times are stochastic and follow a phase-type distribution [24] .The advantage of a phase-type distribution is that its Markovian nature allows for an exact analysis and performance evaluation.
The contribution of this work is to construct an exact evaluation for the can-order joint replenishment policy while taking endogenous lead times into account, and hence to study the complex interplay between mutually dependent product inventory levels.We provide illustrations of the interaction between the inventory properties for the two product case, in terms of lead times, order quantities, and cost reductions from a can-order joint replenishment policy as compared with an independent order-up-to policy.The analysis is applicable to more products but leads to state space explosion.
The structure of this paper is as follows: in Section 3 we define the productioninventory model and then, in Section 4, define the associated Markov process characterizing the production-inventory system.From the matrix exponential form of the steady state distribution, we construct the associated fluid process in Section 5.In Section 6 we present results required for the computation of the expected costs, outline the numerical method, and present some individual and joint replenishment outcomes, illustrating the impact of different parameter values.Concluding remarks are given in Section 7.

Background
Inventory analyses are generally aimed at minimizing ordering and holding costs while satisfying demand and there is a vast literature addressing differing assumptions on how customer orders and production systems are managed.We mention just a few variations.For example, inventory control policies assuming independent and identically distributed inter-event times have received most attention, but a recent study of production policies to minimise expected holding and backlog costs, while considering correlated inter-arrival and processing times, showed that effective production control policies should take correlations in service and demand into account (Dizbin [12] ).We use a model that assumes customer orders are held until the item is produced.Much research seeks ordering policies to maintain the inventory position at a constant level, but a recent analysis illustrates opportunities for performance improvement by using the current inventory level to set a dynamic target for inventoryin-transit, and to place orders that follow this target (Stolyar [36] ).Other work considers situations where a customer may have a ready alternative source and sales are lost (Faaland [14] ), and situations where the supplier has access to an emergency supply to raise the inventory level (Barron [7] ).
In the context of an inventory managing multiple items, the joint replenishment problem deals with the prospect of saving resources through coordinated replenishments, by considering interdependencies amongst product orders involving a single supplier.Studies of joint replenishment have led to many models, solution methods and applications [8] .As the complexity of many variations of the joint replenishment problem has been determined to be NP-hard [11] , studies focus on approximate and heuristic methods.
Replenishment models cover both deterministic demand types, for examples see [8] , and stochastic demand types, for example see [6,21,22] .With the stochastic demand models, there are two main replenishment policies: (i) the periodic replenishment policy that considers a Poisson distributed demand, in which all products follow the same replenishment periodic interval, and each one is replaced to a pre-determined inventory level [4] .(ii) the continuous review can-order policy, in which if any item reaches a mustorder point, a replenishment order is made, and if the other items have inventory levels below a can-order point, they are included in the same order [5,20,30] .On stochastic demand models based on the can-order policyi, solution methods have focused mostly on determining optimal values for re-order and can-order parameters in relation to the selected policy [8] .
Markov processes are valuable as approximate processes, due to their versatility in matching key statistical properties of the demand process [10,26] .Phase-type distributions [27] are often used, and in which their matrix structure supports efficient computational implementation [7] .When general distributions are appropriate for modeling the demand, phase-type distributions can be taken into account in a natural way, as any non-negative continuous distribution of the probability can be approximated with a phase-type distribution [3] .
Prior studies of models incorporating endogeneous lead times include: early work on assemble-to-order production-inventory systems with stochastic lead times, in which one must allow for correlated demands by jointly managing inventories and production capacities across several items -that work used matrix geometric methods for exact analysis to evaluate order-fulfilment performance measures (Song, Hu, Liu [34] ); performance analysis of a single product production-inventory system with periodic demands (Boute et al. [9] ); and more recent work studying a single product dual-source inventory system where the lead times at both sources are stochastic and endogeneous, namely sojourn times in a queueing system -where the authors establish a new approach to finding an optimal policy (Song et al. [35] ), in part drawing on results on endogeneous stochastic supply networks (Song and Zipkin [33] ).
A key element influencing the complexity and usefulness of any Markov model is the choice of the state space and transition function.In our model, the state representation records information on the type of order, the inventory position when the order was placed, and the current state of production or setup, and the analysis takes advantage of matrix analytic and fluid queue results applicable to the associated age process.An approach similar to ours was used by Noblesse et al. [28] for the single item problem, but the analysis in the current paper allows for systems with batch demand arrivals, variable order quantities, and non-exponential production times.Related approaches in less complex settings than that studied here have been used by Horvath and Van Houdt [19] , Otten et al. [29] , Van Houdt [37] .

Model description
Consider a production-inventory system with j ∈ {1, 2} products, where demands arrive according to a compound Poisson process with rate λ (j) .The sizes of the batch demands are i.i.d. and follow a general discrete distribution with maximum demand size m j .Let d x denote the probability of a demand of size x for product j.If inventory is insufficient to fulfill demand, unmet demand is backlogged.
The inventory of product j is replenished according to a (s j , c j , S j ) policy, where c j ≥ s j : if the inventory position of product j reaches its reorder point s j , product ℓ ̸ = j will also place an order if its inventory position is at or below the can-order level c ℓ , otherwise it will not join the order.Setting c j = s j , the (s j , c j , S j ) policy reduces to an (s j , S j ) policy, therefore our analysis holds for both joint and independent replenishment policies.The maximum order size is o (j) = S j − s j + m j − 1.
Orders are sent to the production queue and produced on a first-come-first-served basis by a single machine which produces the units of each order sequentially.Each order requires a major setup time and a production time per unit.If a joint order is placed, only one major setup time is required and an additional minor change-over time.Typically, the minor change-over time will be smaller than the major setup time, although the model does not require this assumption.
When both the setup is complete and the last unit of the order is produced, the order is replenished in inventory.If a joint order is placed, the order is only delivered in inventory when the production for both products is completed.We assume the unit production time is identical for both products and follows phase-type distribution P H(γ p , U p ), that is, with subgenerator matrix U p and initial probability vector γ p of order n p .The major setup time and minor change-over time are assumed to both follow phase-type distributions with P H(γ + , V + ) of order n + , and P H(γ − , V − ) of order n − , respectively.
Denote by η j the expected number of orders placed per unit time of only product j and η c j the expected number of joint orders per unit time initiated by product j.An inventory holding cost h j is incurred per unit per time and Φ + j is the expected number of units in inventory of product j at a random point in time.A backlogging cost p j is incurred per unit per time and Φ − j denotes the expected number of units of product j backlogged at a random point in time.The expected total cost per unit of time for both products is where a major fixed order cost k is incurred per order placed and a minor fixed order cost k j if product j is included in the order.For ease of reference, a summary of key notation is provided in Table 1.
The objective is to determine the inventory parameters (s j , c j , S j ), j ∈ {1, 2}, that minimize E[C].We next describe the Markov process of the underlying model, and key theoretical results are set out in the following two sections, enabling the iterative procedure used for the numerical solutions discussed in Section 6.3.

Markov process of the two product system
Below we describe the state representation and transitions (Section 4.1).For the busy process we use a state space construct similar to that in He [15] and Noblesse et al. [28] .Then we present the transition rate and density matrices (Sections 4.2 and 4.3), and introduce the product utilisation rate, ρ p (Section 4.4) that is used in the inventory level distribution that is derived in Section 5.3.

Process description
When the server is idle, we use as state space I = {(i, j) : where each element of I characterizes an inventory position.Clearly I has dimension (S 1 −s 1 )(S 2 −s 2 ).A Markov process (I t : t ∈ R + ) then records the inventory depletions due to arrivals of demands.
When the server is busy, we represent the system using the Markov process (X t , L t : t ∈ R + ).The first component X t ∈ R + represents the age, i.e., the elapsed time since the order in service was placed, and the second component L t ∈ P represents the state of the order in service at time t.This busy process describes the transitions in the MMAP[K]/PH[K]/1 queue that models the production facility (He [17] ).
For an order in production, the state space P is defined as P = P 1 ∪ P 2 , where the states in P j describe the order when initiated by product j, i.e., the product that reached its reorder point, and we set P j = P j,indiv ∪ P j,joint , partitioned according to whether the order in production is an individual or joint order.
For an order initiated by product 1, define where i j represents the inventory position of product j when the order in production at time t was placed, i.e., at time t − X t , and the third component, z, represents the progress of the production: i.e., it is either in major setup (in a state in z ∈ {1, . . ., n + }), or in production with r − 1 remaining items to be produced to complete the order and the server is currently in state k ∈ {1, . . ., n p }.So z = n + +(r −1)n p +k. Accordingly, z 1,indiv = n + + o (1) n p .Define In contrast to P 1,indiv , here we do not need to keep track of the inventory position i 2 as it will increase to S 2 irrespective of the value of i 2 .We note that the range of z is now larger as the order also involves at most S 2 − s 2 − 1 items of product 2 and the server can also be in one of the states of the minor change-over time.If z ∈ {1, . . ., n − }, the server is in state z of the minor change-over time, if z ∈ {n − + 1, . . ., n − + n + } a major setup is ongoing and the remaining z values correspond to the production of an item.
When the inventory position of product 1, exceeds c 1 , the order only includes product 2, while for inventory positions up to c 1 , the order is joint with product 1.So, for an order initiated by product 2, define where z 2,indiv = n + + o (2) n p and We have defined the model from the perspective of product 1 and thus the sets P 1,indiv and P 2,indiv are structurally different.This representation provides sufficient information to derive the inventory distributions, but enables smaller matrix calculations than if we retained the same structure for P 2,indiv as P 1,indiv .
We also note that to obtain a Markov process it is not necessary to keep track of the inventory position, i 1 of product 1 (as the inventory position increases to S 1 irrespective of the value of i 1 ), but we have included this in the representation as it can then be used in deriving the inventory level distribution of product 1 in Section 5.3, which in turn is used in our iterative procedure to compute the expected cost in Section 6.3.
The process evolves as follows: as the production of an order continues, X t increases linearly over time and a downward jump in X t occurs when an order completes service.Three types of transitions can occur from state (X t , L t ) = (x, i): • A transition to state (x, j) for j ̸ = i ∈ P when the state of production/ changeover/setup changes but the same order remains in service.Denote by (A 0 ) i,j as the rate of this transition.• A transition to a state (y, j) with y ∈ [x − u, x), for 0 < u < x when an order completes service and there is a subsequent order in the queue enters service.When the inter-arrival time of the subsequent order is (at most) u, time units, the order has waited at least x − u time units, illustrated in Figure 1.Denote by A i,j (u) the transition rate, with dA i,j (u) its corresponding density function.
The matrix A(u) of rates A i,j (u) takes into account the correlation between j (which represents the new order size) and u (the inter-arrival time between the new order and the previous one).• A transition to state (0, j) when the order completes service and the queue is empty upon service completion.This occurs with rate ∞ x dA i,j (u).The matrices A 0 and dA(u) are described in the next sections.

Transition rate matrix A
where ++ is an order |P j | matrix that contains the rates when the order was initiated by product j.Then where I m1 records the inventory position of product 1 the moment the order was placed, the matrices F ++,joint and F ++,indiv describe the joint and individual orders initiated by product 1, respectively, and ⊗ denotes the Kronecker product.Then, is an order z 1,joint matrix and describes the transitions in production and setup phases and where the vectors u p = −U p e and v + = −V + e denote the completion rates of the unit production and the major setup time, respectively.Define, F an order (S 2 − c 2 )z 1,indiv matrix.Since there is no minor setup time, we record the inventory position of product 2 with I S2−c2 , the moment the order was placed.Then we obtain a zero matrix, except for the S 2 − c 2 sub-matrices on the diagonal where each of sub-matrix defines the transitions in production and setup and specifying i 2 .
Similarly, define where is of order (S 1 − c 1 )z 2,indiv and is of order (c 1 − s 1 )z 2,joint .In both joint and individual orders, we keep track of the inventory position of product 1 at the moment the order was placed.

Density matrix dA(u)
The density function of the rate matrix A(u) can be expressed as where F +− describes the transitions from the states in P prior to service completion to the inventory positions of both products immediately after the order was placed.
The matrix exponential e F−−u gives the density of the inter-arrival time of u, where F −− describes the evolution of the inventory position until the subsequent order is placed.The matrix F −+ describes the transitions to the initial production state of the subsequent order beginning at t − X t + u.

Production completion of the order in service (F +− )
Denote by α j the vector of the steady state probabilities of the inventory position immediately after a replenishment, then α j = (1, 0, . . ., 0) with length (S j − s j ).Let +− describes the transitions when the order was initiated by j.Recalling that completion occurs after a major setup time in an individual order, or after a minor change over time of a joint order, we have with dimensions |P| × (S 1 − s 1 )(S 2 − s 2 ) and where ν − , 0, . . ., 0) ⊺ a vector of length z j,joint , where v ⊺ + denotes the transpose of v + .The vector e m1 records the inventory level of product 1 during production as a new order is about to begin.The upper sub-matrix in Equation (11)  describes an individual order of product 1, then the inventory position of product 1 is raised to S 1 and we record the inventory position of product 2 with I S2−c2 .The lower sub-matrix describes a joint order, in which case the inventory positions of both products are raised to their respective S j .Similarly, where I S1−c1 is used to keep track of the inventory position of product 1 when only product 2 is replenished, and e c1−s1 records the inventory position of product 1 in a joint order.

Inventory position until subsequent order (F −− )
The matrix F −− is of order (S 1 − s 1 )(S 2 − s 2 ) and describes the transition rates in the inventory positions of both products due to demand arrivals until the reorder point of one of the products is reached.Then . . .λ (1) d (1) S1−s1−1 0 −λ (1) . . .λ (1) d The time u between subsequent orders is then equivalent to the time it takes for the process to decrease from the inventory positions just after an order was placed until absorption in one of the reorder points (or below), thus its density is given by e F−−u .
4.3.3.Start of production of the subsequent order (F −+ ) −+ ] where F −+ records the transitions when the order was initiated by product j and which has dimensions where is a column vector η k of length S j − s j and contains the rates of demand arrivals for product 1 which result in an order of size S 1 − s 1 + k.The matrix Γ (1) k determines the state in which production starts if the order size of product 1 equals S 1 − s 1 + k and if the inventory position of product 2 is S 2 −i 2 +1 at the moment when product 1 initiates an order, for i 2 = 1, . . ., S 2 − s 2 .If only product 1 is replenished, Γ k also records the inventory position of product 2, i 2 .Thus Γ (1)   k is an (S 2 − s 2 ) × |P 1 | matrix with all its entries equal to zero, except for n p entries given by γ p .For i 2 = 1, . . ., S 2 − c 2 the inventory position of product 2 remains above c 2 and the n p nonzero entries of row i 2 of Γ (1)   k are on the columns corresponding to the states (k, we have a joint order and the n p nonzero entries on row i 2 of Γ (1)   k are on the columns corresponding to the states (k, n − s +n + s +z) of P 1,joint with z = (S 1 −s 1 +k+i 2 −2)n p +1 to (S 1 − s 1 + k + i 2 − 1)n p as the joint order has size Similarly, F −+ defines the transition rates when product 2 reaches its reorder point, F (2) k an (S 1 − s 1 ) × |P 2 | matrix with all its entries equal to zero, except that each row contains n p entries equal to γ p .More specifically, the non-zero entries on row i 1 ∈ {1, . . ., S 1 − c 1 } appear on the columns corresponding to the states (i

Utilisation of the production system
Define ρ p as the production utilisation rate, that is, excluding the setup times, then where γ p (−U p ) −1 e is the mean unit production time and is the mean demand size for product j.The overall utilisation ρ is then determined by the number of orders placed, where −γ + V −1 + e and −γ − V −1 − e are the mean major setup and minor change-over times, respectively.

Steady state analysis
We consider the form of steady state of the Markov process in Section 5.1 and construct the associated fluid queue in Section 5.2.In Section 5.3 we obtain an expression for the inventory position and the number of demand arrivals since an order was placed and combine the results to determine the inventory level distribution.

Preliminaries
For x > 0, j ∈ P, denote the steady state density of (X t , L t ) as and the corresponding vector δ(x) = {δ j (x), j ∈ P}.If ρ < 1 the process (X t , L t ) has matrix exponential form (Sengupta [31,32] ) and there exists a matrix T of order |P| such that where the matrix T is the smallest non-negative solution to and δ(0) = θ(−T ), where θ is the stationary vector of See Van Houdt [37] for details.
Using the matrices A 0 and dA(x), the matrix T can be computed iteratively, T n+1 = A 0 + ∞ 0 e Tnx dA(x), with T 0 = 0 (Sengupta [31] ).As this method results in linear convergence and is therefore impractical under high loads, in the next section we present a more efficient approach.

Reduction to a fluid queue
To employ a method for computing T that converges quadratically we construct a fluid queue (Dzial et al. [13] ).In order to obtain a fluid queue, the fluid process must be skip-free in both directions so we replace the immediate downward jumps in X t by intervals of the appropriate length during which the level decreases linearly.Then there are |P| phases in which the fluid increases and |I| artificial phases in which the fluid decreases and the rate matrix of the underlying continuous-time Markov chain is The matrix F ++ = A 0 contains the rates at which the phase changes while the fluid increases (the same order remains in service), F +− contains the transitions from an increase to a decrease in fluid (a service completion occurs), F −+ contains the transitions from a decrease to an increase in fluid (a new order enters service) and F −− contains the transition rates while the fluid decreases (demands arrive depleting the inventory, but no product has reached its order point).
When the fluid queue hits level zero the phase continues to evolve according to F −− until an event part of F −+ occurs and the fluid level becomes positive again.
Censoring the queue only on the periods where the level increases we obtain the Markov process of Section 4, in which the age always increases, unless there is a downward jump.Due to Equation ( 9), a downward jump to level zero occurs with rate ∞ x dA i,j (u) = F +− e F−−x (−F −− ) (−1) F +− i,j where the part (−F −− ) (−1) F +− captures the evolution of the phase of the fluid queue at zero.
If we take the expression for the steady state of a fluid queue (Latouche [23] ) and observe the queue only when the level is increasing, the steady state has the form given in Equation (20) with and δ(0) = θ(−T ), where θ is the stationary vector of and matrix Ψ, the first return probabilities to the initial level, is the minimal nonnegative solution to an algebraic Riccati equation To reduce the computation time of θ, we employ the power method (Van Houdt and Perez [39] ).Further, as F ++ has block diagonal form, we can apply the algorithm of Meini [25] , which exploits this structure to reduce the memory and time complexity.
Observing the phase process only during the intervals of time in which the fluid level is decreasing, that is, when the production facility is idle, the rate matrix is from which we obtain the steady state probabilities of the idle production facility, θ, such that θD = 0, θe = 1.

Inventory level distribution of product 1
For i ≥ 0, denote by ϕ i the probability the inventory level is S 1 − i.To derive this distribution, we separately consider the situations in which the server is busy and where it is idle.Firstly when the server is busy, define q k,n as the probability that the inventory position of product 1 was S 1 −k at the moment when the order in production was placed and n demand arrivals of product 1 have depleted inventory since.The corresponding vector q n = {q k,n , k = 0, 1, . . ., o (1) } is partitioned such that where the entries refer to individual orders placed by product 2, joint orders initiated by product 2, and orders initiated by product 1 (joint or individual), respectively.During the time x this order has spent in the system, the probability that the inventory is further depleted by the arrival of n demands for product 1 is (1/n!) λ (1) x n e −(λ (1) x) . Let λ (1) x n n! e −λ (1) x dx = δ(0) λ (1) n λ (1) then q n = p n (x) I m1 ⊗ e |P1|/m1 , (I c1−s1 ⊗ e z 2,joint ) , (I S1−c1 ⊗ e z 2,indiv ) .We can then write where 1 − ρ is the probability of a idle server in which the inventory level of product 1 is at least s 1 + 1 and the idle probabilities g i are given by where θ is derived from the fluid queue construction in Equation ( 26) and where When the server is busy, with probability ρ, ϕ i is obtained by q k,n and the n-fold convolution of the demand size distribution, ks . ( Then the expected number of units in inventory and expected number of units backlogged are Denote by L the lead time and define ζ(x) as its steady state density, where η = F +− e, and with corresponding distribution function then the expected lead time of a random order is E[L] = −θT −1 η(θη) −1 .Conditioning on the order if product 1 is included, the product dependent expected lead time is where using the expression for F +− in Equation ( 12), and where • denotes the Hadamard product.

Expected total cost
For the computation of the expected total cost, we need to characterize the time between orders.In Section 6.1 we consider the time, η j , between individual orders, and similarly in Section 6.2 for the time, η c j , between joint orders.In Section 6.3 we present the iterative procedure employed to derive the policy parameters for both products and give numerical illustrations.

Expected number of individual orders per unit time
The time between two subsequent individual orders has distribution P H( π j , F (j) −− ) of order (S 1 − s 1 )(S 2 − s 2 ) and the time to place the subsequent order is equivalent to the time to reach absorption.The expected time between subsequent individual orders of product j is given by η The matrix F −− records the transitions in inventory position until the subsequent replenishment, thus only a decrease in inventory positions is possible.However in F (j) −− , an increase in inventory positions is possible when product i ̸ = j places an order.Then for j = 1, where the second term refers to an increase in inventory positions due to a joint order initiated by product 1, and the third term refers to a (joint or individual) order placed for product 2.The vector π 1 contains the probabilities of starting in a given state following an individual order of product 1 and is the stationary distribution of , as an individual order for product 1 can only occur if the inventory position of product 2 is at least c 2 + 1.Similarly, F and π 2 is the stationary distribution of

Expected number of joint orders per unit time
The time between two subsequent joint orders initiated by product 1 has distribution e.As the inventory position of both products after placing a joint order is (S 1 , S 2 ), the initial state vector is α 1 ⊗α 2 and where the second term accounts for the individual orders placed of product 1, and the third term refers for the orders placed by product 2.The time between two subsequent joint orders initiated by product 2 is

Iterative solution and numerical examples
To obtain the inventory parameters that minimize E[C] in Equation ( 1), we use an iterative procedure to obtain successive approximations for both products using neighbor search and steepest descent algorithms.Given a set of initial inventory parameters for product 1, we determine the optimal inventory parameters for product 2. Then given this inventory policy for product 2, we determine the inventory parameters for product 1 that minimize E[C] and repeat until convergence is reached.We note that to compute the cost for both products simultaneously would involve matrices with significantly larger dimensions and hence be less tractable.We present results for six experiments for different values of the ordering costs and minor change-over times, as these drive the order quantities, (S j − s j ) and (S j − c j ).The parameter values are chosen to highlight the impact of changing the key replenishment drivers and costs, being the ordering costs and setup times.The holding and backlogging costs represent a service level of 90%.Parameter values common to all experiments are listed in Table 2 and the different values are listed in Table 3.The unit production time, major setup and minor change-over times are each exponentially distributed and the demand size distribution is zero-truncated binomial.
We demonstrate the iterative procedure for deriving parameters for the case of an independent replenishment (s j , S j ) policy and discuss the interplay between the quantities of interest.We then illustrate the cost reduction by using a joint replenishment policy compared to independent policy.
Obtaining cost minimizing parameters for individual replenishment.For Experiment 1, we outline the solution for (s * j , S * j ).The iterations are listed in Table 4 with an initial guess of S 1 − s 1 = 5, from which we obtain (s * 1 , S * 1 ) = (16, 38) and (s * 2 , S * 2 ) = (11, 27).To understand how the quantities of interest change throughout the iteration, Figures 2-4 illustrate how ρ, E[C] and s j evolve.
Iteration 1.For the initial order quantity S 1 − s 1 = 5, the small size leads to multiple setups, congested production, long lead times and high costs.The solid line Figure 2(a) illustrates how ρ decreases as the order quantity of product 2 increases.The cost-minimizing reorder point s 2 as a function of the order quantity is illustrated in Figure 3(a) and Figure 4(a) illustrates the impact on E[C].Due to the small order of product 1, a large order of product 2 is optimal, as this tempers the utilisation rate and lead times.Iteration 2. Given S 2 −s 2 = 52, the solid lines of Figures 2,3,4(b) illustrate the impact of the order quantity of product 1 on the utilisation rate ρ, optimal reorder point s 1 and expected total costs E[C], respectively.Increasing the order quantity of product 1 from 5 to 19, leads to more favorable lead times (due to a reduction in setups) and lower total system costs.Iteration 3. Given S 1 −s 1 = 19, the dashed lines of Figures 2,3 The results for the all experiments are listed in Table 5.As expected, the minor change-over time has no impact as no joint orders are placed.Different values of k and k j have no impact as long as k + k j remains the same for product j, which explains why the result of experiments 1 − 4 are the same for the (s j , S j ) policy.In experiments 5 and 6 k + k j = 0 and the optimal order quantities are only slightly smaller compared to experiments 1 − 4 due to the impact of the order quantities on setup times and lead times.
Comparison of independent and joint replenishments.Following a similar iterative procedure, we present the solutions for (s j , c j , S j ) in Table 5.In Table 6 we give the cost reductions of (s, c, S) compared to the (s, S) policy, listed alongside the ratios of the fixed ordering costs k/k j and the mean setup time vs. the mean change-over time.As expected, E[C] is lower in the joint replenishment scenarios.
In Experiment 6, due to the small change-over time, we find the highest gains with a cost reduction of 10.73%, in which the benefit is mostly seen from the reduction in setup time by placing joint orders.We see the benefit of joint replenishment is not only influenced by the ratio of k/k j , but also by the ratio of the joint major setup time and the changeover time.Comparing the individual and joint policies we find that the mean number of orders placed reduces in our set of experiments by more than 14%.
We demonstrate that this cost discrepancy increases in the ratio of k and k j and in the ratio of the major setup time and the minor change-over time.The impact of a can order policy on lead times is unclear.Although the number of setups reduces, lead times may not necessarily go down due to a batching time, i.e., waiting time after production until the joint order is completed.

Concluding remarks
The analysis of sufficiently detailed production-inventory models is challenging when accounting for the number of variable quantities in the supply-chain process.We consider a continuous review finite capacity production-inventory system with two products in inventory.As the order process of each individual product impacts the utilisation rate of the production system and thus the lead times for all products, the inventory levels are mutually dependent.
The contribution of this work is to construct an exact evaluation for the can-order joint replenishment policy while taking endogenous lead times into account and hence, using numerical examples, to study the complex interplay between product inventory levels, adding to the relatively limited studies of multi-item inventory analysis.The analysis in the current paper allows for systems with batch demand arrivals, variable order quantities, and non-exponential production times.Modelling the production facility as a multi-type queue, with correlated arrivals and type dependent service, the matrix exponential form of the age process of the busy queue is used to compute the inventory level distribution.For a given set of input values, our analysis generates parameters to specify a minimal cost can-order policy, a well-known heuristic for the joint replenishment problem.
We use numerical methods to obtain inventory parameters that minimize total ordering and holding costs and provide examples of the interaction between the inventory properties for the two product case, in terms of lead times, order quantities, and cost reductions from a can-order joint replenishment policy as compared with an independent order-up-to policy.We demonstrate scenarios where there is decreased cost of a can-order (s, c, S) joint replenishment policy as compared with an independent order-up-to (s, S) policy.Our analysis also can apply to a setting with multiple retailers of a single product that is produced on the same production line (Zipkin [41] ) and joint replenishment in transport applications (Padilla Tinoco et al. [30] ).A limitation is the computation time, which grows rapidly with the number of products, the maximum demand size or the order quantity.The analysis is applicable to more products but further work includes establishing optimality of the solution obtained numerically, and development of heuristics to manage combinatorial state space explosion.The analysis can be extended to a system with MAP arrival processes, however this will increase the dimension of the matrices involved (by a factor equal to the product of the number of states of both MAPs).Notation Definition λ (j) demand arrival rate for product j m j maximum batch demand size d (j) x probability of a demand of size x o (j) maximum order size P H(γp, Up) unit production time of order np P H(γ + , V + ) major setup time of order n + P H(γ − , V − ) minor change-over time of order n − η j expected number of orders placed per unit time of only j η c j expected number of joint orders per unit time placed by j Φ + j expected number of units in inventory of j Φ − j expected number of units of j h j inventory holding cost per unit per time p j backlogging cost per unit per time k major fixed order cost per order k j minor fixed order cost if j is included in the order

Collated Author's Response
We would like to thank the reviewers and the editor for their constructive and detailed feedback and the opportunity to resubmit.We have performed a thorough revision of the manuscript, improving the presentation of the paper.We address each specific comment from the editor and each reviewer in the tables below, and also provide the latexdiff comparison from the the previous submission.
Editor comments Response R1 recommends acceptance, subject to some minor revision.The referee mentions that the paper is not easy to follow, but does not have any good idea how to improve it, and therefore leave it as is.R2 is also positive about the technical contribution of the paper.He/she too complains about the paper's writing and exposition.This referee recommends a major revision to improve the readability, clarity, and exposition, and provides some detailed suggestions.R2 feels that the paper will be acceptable after the paper's clarity and exposition has improved to satisfactory level.I would like to thank the referees for their valuable feedback.I concur with the referees assessment of the paper.Given that both reports are consistent in nature, and the remaining weakness is on writing and exposition, I recommend a minor revision.But please take this revision seriously, your contribution can be appreciated and impactful only when the readers can understand it.

Reviewer 1 comments
Response 1.It seems to me that the paper is mathematically correct.2. The main contribution of the paper is the introduction of a Markov process for an inventory model with two types of products.The introduction of the Markov process, and a Markov modulated fluid flow process, is highly technical and creative.Using the stationary distribution of the Markov process, a number of quantities for the inventory model are obtained, which makes it possible to find a good inventory control policy.3. It is not easy for me to follow the details for the introduction of the Markov process.On the other hand, I am not able to suggest a better way to do it.Thus, the presentation of the paper, especially about the Markov process, is acceptable to me.In summary, the methodology utilized in this paper is innovative and the analysis is mathematically correct.I recommend the paper be accepted for publication in Stochastic Models.Minor editorial comments.1) Page 7, lines 17-18: Elaborate English.
We have added a sentence on the minor change-over time.
2) Page 7, line 43 and line 47: Why is "where" used?Is it better to use "and"?
We have made this suggested change.4) Page 8, lines 17 -24: I do not fully understand the meaning of this sentence.I think it is related to the introduction of the Markov process.For example, the sets P 1,indiv and P 2,indiv are structurally different.This paragraph has been rewritten to improve the explanation and a further comment has been added in Section 4.1 as to why they are structurally different.5) Page 8, line 46: Is it better to define the "state space I" by (i, j), s 1 + 1 <= i <= S 1 , 2 + 1 <= j <= S 2 ?
This improvement has been made.6) Page 10, lines 36 -43: It seems to me that this comment is true for the joint order case (for both i 1 and i 2 ).P 1,indiv has i 1 and i 2 , but P 2,indiv does not have i 2 .(I think it is related to page 8, lines 17-24.) We have expanded the statement to improve the clarity.
The vector e m1 was defined, but the paragraph has been rearranged to improve clarity.8) Page 14, line 43: Is the notation "T " for matrix transpose defined?Definition added.We have added further detail of the boundary behaviour on page 13. 10) Page 21, equation ( 28): Please make sure that the elements in q n are arranged in consistent with that in equation (26).
The ordering is consistent.We agree with the reviewer and have modified the title of Section 6.3.16) Page 29: Any comments on whether or not the analysis can be extended to a system with MAP arrival processes for the two types of products?
We have added a comment on the extension to MAP arrival processes in the conclusion -noting that this will increase the dimension of the matrices involved (by a factor equal to the product of the number of states of both MAPs).

Continued on next page
Reviewer 1 comments Response 17) Page 29: Any comments on whether or not the analysis can be extended to a system with different unit production times?
We note that in principle this could be done, noting that the number of phases in the set P 2,joint increases quite significantly as we need to keep track of both order sizes, instead of just the sum of the two.So the size of for instance the F ++,joint matrix would involve the product of o (1) and (S 2 − s 2 − 1) instead of its sum.

Reviewer 2 comments
Response The paper analyzes a continuous-review finite capacity production-inventory system with two products in inventory.Two replenishment policies are considered: the order-up-to (s, S) policy and the can-order (s, c, S) policy.The main tool is the matrix geometric method.Their aim is to derive the inventory parameters that minimize the total ordering and inventory related costs.Numerical examples, comparison and insights are also provided.The model seems to be important, practical and can be applied in many cases.However, some serious changes should be made before acceptance.The Introduction and the description Sections are poorly written, and thus, it makes the manuscript difficult to understand.The introduction isn't clear, cumbersome, lacks necessary details, explanations and recent studies.The analysis should be restructured.Numerical examples, observations and insights should be expanded.Also, the language needs some improvements; I recommend the use of a technical editor.I recommend a major revision; Following, I detail my comments.Introduction.Motivation, Examples, Contribution.The authors should motivate their model and add real-world examples of the model.There are no real-world examples or applications of the two policies.In fact, the policies are shortly described rather than extendedly.The contribution of the paper is not clear.Although there is a list of good papers, they are not described, and their contribution is not clear.As a result, the contribution and the uniqueness of the paper is not clear.In addition, no main results, observations, conclusions and insights are provided in the Introduction.
Section 1 has been reworked to emphasise the motivation and contribution, and added comments on real world interest.

Figure 1 .
Figure 1.: Sample path of X t .Upon production completion of the n th order at age x, the subsequent order in service has already spent at least x − u time units in the system.

9 )
Page 19, Section 5.2 [old reference Section 4.1] : It is not clear how the border x = 0 (zero fluid level) is defined.How does the underlying Markov chain changes when the fluid level enters or leaves zero.

11 )
Page 21, equation (30): Is there a short justification for it?Section 5.3 has been rearranged to improve the clarity of the construction of the inventory level distribution.12) Page 22, equation(35): Please check the correctness of the column vector with elements {1, 0, 1} on the right hand side.This has been checked.13) Page 23, line 13, the title of Subsection 5.1: Should "per unit time" be added to it?This modification has been made.14) Page 24, line 45, the tile of Subsection 5.2: Should "per unit time" be added to it?This modification has been made.15) Page 25, line 30, the title of Section 6: It might be better to add "and numerical examples" to it.

Table 1 .
: Summary of key notation

Table 2 .
: Parameter values common to all numerical experiments

Table 3 .
: Ordering costs and setup times for each experiment

Table 5 .
: Optimal replenishment parameters for each experiment