Petri net analysis of a queueing inventory system with orbital search by the server

1 Introduction

Stochastic management inventory systems have many practical applications and have been studied extensively in the literature (see [1–4]). These systems play an important role in real-life situations (in manufacture, warehouse, supply chains, etc.). The complexity of stochastic inventory systems depends on various characteristics:

• Single product inventory and multi-product inventory

• Sources (finite, infinite)

• Arrival demand (batch arrival, generally distributed, etc.)

• Service time and lead time (positive/instantaneous exponential/non-exponential, etc.)

• Ordering policies: ( s , S ) , ( s , Q ) , etc.

• Retrial, vacation, breakdown, etc.

In most works on the inventory systems, the authors assume that the population of demands is very large. However, in several practical situations, the number of demands who access the system is finite [2,12–15]. Thus, when a source is free at time moment t , it may generate primary demands during interval ( t , t + h ) with probability λ h + 0 ( h ) . In [2], a continuous review ( s , S ) inventory system with retrial, service facility, and finite sources of demands is analyzed.

Several research articles have dealt with the inventory systems with positive or zero service time/lead time. Artalejo et al. studied the retrial inventory systems with immediate service [5]. The notion of inventory model with positive service time was initiated by Sigman and Simchi-Levi [16]. Shajin and Krishnamoorthy obtained the product-form solution for M / M / 1 retrial queueing inventory system with positive service time [17]. An extensive survey on inventory systems with positive service time could be found in [18].

The lead time may actually be uncertain in duration due to variability in shipping times, material availability, and supplier processing times. The inventory systems with instantaneous lead time were studied by Liu and Yang [19], Berman and Sapna [20], etc. The previous works have extended to include exponential lead time (see [5,21,22]). The inventory models with positive lead time are complex to analyze. Still more complex are the models in which the lead time has a general distribution [23–25], the reason for which is that there have not been significant studies on retrial inventory systems with arbitrary distributed lead time.

The applicability of Petri nets for stochastic analysis of inventory queueing system has received few attention [15,26]. In [15], the authors have analyzed priority multi-server retrial inventory queues with Markovian arrival process generated by single or dual sources and exponential service times where both queueing and orbit spaces are assumed finite using the generalized stochastic Petri net (GSPN) formalism. In [26], Chen et al. studied a finite source inventory system using the batch deterministic and stochastic Petri nets (SPNs).

The dynamic [27] evolution of inventory systems proceeds from one discrete state to another at arbitrary moments in time. The steady-state of stochastic inventory systems can be obtained either by discrete event simulation or by stochastic process theory. Most works on quantitative evaluation of these systems have addressed models with exponentially distributed durations, which always satisfy the Markov property. In this case, the underlying stochastic process of the model is a continuous time Markov chain (CTMC), and standard numerical algorithms can be used to compute the state probabilities. From these systems, the interesting performance measures of the system can then be derived. For multi-dimensional Markov chains underlying to the inventory systems, an appropriate numbering of the states yields a repetitive structure infinitesimal generator matrix, which can be formulated as birth-and-death process, quasi-birth-and-death process [28,29], M / G / 1 -type, etc. This repetitive structure of the infinitesimal generator matrix allows the calculation of the state probabilities by approaches known as the matrix-geometric method, matrix analytic method [30], etc.

However, in several application contexts, the inventory systems are characterized by deterministic timers or non-exponential durations (e.g., arrival pattern, the service pattern of servers, lead time, and lifetime). In this case, the underlying stochastic process of these systems falls in more complex classes of the so-called non-Markovian processes but solution may still be viable if the model satisfies the Markov property at specific times (regeneration points). The traditional approaches to solve non-Markovian models use phase-type expansions, apply the method of supplementary variables [31], or construct an embedded Markov chain [32]. The supplementary variable method is known to be originated by Cox [33] and has become one of the most frequently used approaches for both the continuous and discrete-time queuing systems. The framework of this method in stochastic models leads to the notion of the partial differential equations, ordinary differential equations, and integral equations [34]. The supplementary variable method is known as an efficient way of deriving the steady-state solution for the stochastic process underlying a Markov regenerative stochastic Petri net (MRSPN) formalism [31].

This article is an extension of the work [35] presented in the conference “IEEE International Conference on Recent Advances in Mathematics and Informatics (ICRAMI’2021),” where we have presented an approach for modeling and analyzing the inventory system with finite sources of demands, retrial demands, positive service time, exponential lead time, ( s , S ) replenishment policy, and demand search from the orbit. Immediately after a service completion, the server with probability p makes an instantaneous search for an orbital demand for the next service or remains idle with the complimentary probability ( 1 − p ) . The notion of orbital search for inventory system was introduced with the hope that it would decrease the length of server idle period and not neglected the unsatisfied (orbital) demands. We gave an extensive analysis of this system using GSPN formalism. This high-level modeling formalism allows us to generate the underlying diagram state and the CTMC of the studied system. However, the state space of CTMC underlying our GSPN increases exponentially as a function of the source size, the maximum inventory level, and the inventory level threshold. So, for real applications, the models may have a very large state space (state space explosion). Hence, by exploiting the repetitive structure of the transition rate matrix of our GSPN, we develop a recursive algorithm for automatically calculating the steady-state probability vector and the performance indices, without generating the reachability graph or the CTMC.

In this extended work, we added some numerical results to original work [35] in order to highly illustrate the influence of the system parameters on some performance characteristics. Numerically, we investigate the sensitivity analysis of the search probability, the maximum inventory level, the inventory level thresholds, and the number of sources over the total expected cost rate. Most of the aforementioned works assume that the lead time is instantaneous or exponentially distributed. However, this basic lead time assumption is far from the real situation. For this, we assume that the lead time is generally distributed. So the underlying process to the inventory system considered is a Markov regenerative process. Under this assumption, we suggested to use the MRSPN formalism in order to establish the appropriate model for our system. This formalism allows us to capture easily the interaction between the components of this system, to better understand its behavior, to generate its underlying diagram state, and to compute its performances. To obtain the quantitative analysis of this MRSPN model, we used the supplementary variable method. In addition, we established an algorithm that comprises steady-state solution and computes performance measures. The numerical results for performance measures and the total expected cost rate are presented under the assumption that the lead time follows the following family of distributions: exponential, Erlang, hypoexponential, and hyperexponential.

2 SPNs

SPNs are defined by associating a random variable called firing time for transitions [36]. A random firing time elapses after a transition is enabled until it fires. Many stochastic variants of SPN have been proposed, which can be used as high-modeling formalism for performances evaluation and reliability of systems, like GSPN [37], deterministic and SPNs, MRSPN [38]. According to the type of distribution allowed for the firing time, the SPN may correspond to a wide range of stochastic processes.

The GSPN class contains two types of transitions: immediate transitions and exponential transitions. The marking process ( ℳ t ) t ≥ 0 generated by the GSPN class is a semi-Markovian process; this process includes two types of markings: tangible markings and vanishing markings. The quantitative analysis of the GSPN can be summarized as follows:

• Generating the reachability graph of the GSPN,

• Eliminating the vanishing markings in order to obtain the reduced reachability graph,

• Converting the reduced reachability graph to CTMC with state space H (set of tangible markings),

• Constructing the transition rate matrix A of the CTMC,

• Computing the steady-state probability vector v of the CTMC by solving the linear system of balance equations:

  v A = 0 , v e = 1 ,

  where e is a column vector of 1’s of appropriate dimension and 0 is a vector of zeros of appropriate dimension.

• Computing the characteristics of the GSPN model (mean number of customers in places, throughput of transitions, probability of event, etc.).

• The generator matrix A contains all exponentially distributed-state transitions that do not preempt a GEN transition.

• The exponential-state transitions that preempt a GEN transition are denoted by the generator matrix A ¯ .

• The probability that the firing of a GEN transition leads to state “j,” given that the transition is fired in state i, is represented by the matrix of branching probabilities Δ .

• Step 01. Obtain the two matrices Ω and Ψ :

  Ω = ∑ g ∈ G I GEN ∫ 0 ∞ e A GEN x d F GEN ( x ) , Ψ = ∑ g ∈ G I GEN ∫ 0 ∞ e A GEN x [ 1 − F GEN ( x ) ] d x ,

  where F GEN ( x ) is the firing time distribution of the GEN transition g , G denotes the set of all GEN transitions of the MRSPN, A GEN is the generator matrix of the subordinated CTMC of the GEN transition g , I GEN is the diagonal matrix whose ith element is equal to one if i ∈ G and equal to zero otherwise.

• Step 02. Compute the vector x by solving the linear system:

  x S = 0 , x L e = 1 ,

  where, S = A EXP + Ω Δ + Ψ A ¯ − I GEN is the generator matrix of a CTMC, we refer to it as the embedded CTMC, L = I EXP + Ψ is the conversion matrix, I EXP is the diagonal matrix whose ith element is equal to one if the state i is exponential and equal to zero otherwise and A EXP filtered matrix is defined by: A EXP = I EXP A .

• Step 03. Compute the steady-state probability distributions v of the MRSPN by the formula:

  v = x L .

3 The proposed GSPN model

We consider a continuous review inventory system with finite sources of demands ( 2 ≤ N < ∞ ) , retrial demands, positive service time, positive lead time, ( s , S ) replenishment policy, and demand search from the orbit. The following assumptions and notations are considered:

• The arrival process of primary demands is quasi-random with parameter λ ( > 0 ) .

• The maximum inventory level is S units. The policy used is ( s , S ) , in which an order is placed for a quantity (“ Q = S − s > s + 1 ”) units up to S whenever the inventory level falls to the threshold s or below.

• The lead time distribution is exponential with parameter α ( > 0 ) .

• When the inventory level is zero or the server is occupied, any arriving primary demand joins a virtual room called orbit. These orbiting demands compete for their demands according to an exponential distribution with parameter γ ( > 0 ) . We consider the constant retrial policy (i.e., the retrial rate is independent of the number of customers in the orbit).

• Immediately after a service completion, the server takes a demand from the orbit with probability p , when there are demands in the orbit and the stock is not empty, and with probability ( 1 − p ) , the server remains idle.

• We suppose that the service time following exponential distribution with rate μ ( > 0 ) and the search time is negligible.

Figure 1

GSPN model of the inventory system with finite sources of demands, retrial demands, service time, lead time, and demand search from the orbit.

Our model contains six places p i , i = 1 , 6 ¯ (noted by circles) and nine transitions t i , i = 1 , 9 ¯ . The white rectangular boxes represent the exponential transitions ( t 1 , t 3 , t 7 , and t 9 ), and the thin bars represent the immediate transitions ( t 2 , t 4 , t 5 , t 6 , and t 8 ). The interpretation of the places and transitions of our model is explained in Table 1.

The markings (ordinary states) of our GSPN are given by:

and its initial marking is M 0 = ( N , 0 , 0 , S , 0 , 0 ) . According to the equation “ # p 1 + # p 3 + # p 6 = N ,” we deduce that the vectors (micro-states) M i = ( # p 3 , # p 4 , # p 6 ) provide all information needed for states description of this GSPN model. So, the state space of this GSPN is defined as:

and the total number of its states is equal to “ N ( 2 S + 1 ) + 2 .”

The marking process ( ℳ t ) t ≥ 0 underlying the GSPN depicted in Figure 1 is isomorphic to a three-dimensional CTMC that expresses server status ( # p 3 = 0 if the server is idle and # p 3 = 1 if the server is busy), # p 4 the number of demands in the orbit, and # p 6 the inventory level. We arranged the state space of this CTMC as:

where 00 = ≪ 00 k ≫ , k = 0 , … , N , 1j = ≪ 1 j k ≫ , j = 1 , … , S , k = 0 , … , N − 1 , 0j = ≪ 0 j k ≫ , j = 1 , … , S , k = 0 , … , N − 1 , 0Q = ≪ 0 Q k ≫ , k = 0 , … , N .

The infinitesimal generator A can be expressed in the block form:

where

The dimensions of the entries (sub-matrices) of A are given in Table 2.

4 The steady-state solution

The GSPN model given in Figure 1 is bounded and admits the initial marking M 0 as home state. Then, the steady-state probability distribution v = v ≪ i , j , k ≫ of this GSPN exists and is unique. The elements v ≪ i , j , k ≫ are given by:

Using the block-structured matrix A , the vector v can be represented by:

where

The linear system v A = 0 yields the following global balance equations:

• v 00 A 0 + v 11 B 0 = 0 .

• For j = 1 , … , s ,

  v 1j A 1 + v 0j D + v 1j+1 E = 0 , v 0j A 2 + v 1j+1 B = 0 .

• For j = s + 1 , … , Q − 1 ,

  (1) v 1j A 3 + v 0j D + v 1j+1 E = 0 , v 0j A 4 + v 1j+1 B = 0 .

• For j = Q ,

  v 1j A 3 + v 0j D Q + v 1j+1 E = 0 , v 00 C Q + v 0j A Q + v 1j+1 B Q = 0 .

• For j = Q + 1 , … , S − 1 ,

  v 1j-Q C + v 1j A 3 + v 0j D + v 1j+1 E = 0 , v 0j-Q C + v 0j A 4 + v 1j+1 B = 0 .

• For j = S ,

  v 1s C + v 1S A 3 + v 0S D = 0 , v 0s C + v 0S A 4 = 0 .

where the matrices ℋ ( 0 , j ) and ℋ ( 1 , j ) are given by:

and

And I is an identity matrix.

The vector v 0S can be obtained by solving the two equations:

5 The proposed MRSPN model

In this section, we consider that the lead time is of general distribution with cumulative distribution function F GEN ( x ) . Thus, the suggested model is depicted in Figure 2. The black rectangular box represents the general transition t 7 of lead time.

Figure 2

MRSPN model of the inventory system with finite sources of demands, retrial demands, service time, lead time and demand search from the orbit.

The matrix A is given by:

where

The other matrices: A 3 , A 4 , A Q , B , B 0 , B Q , D , D Q , and E are the same to those obtained for the previous GSPN model.

The matrix of branching probabilities Δ is given by:

The matrix A ¯ is a zero matrix (i.e., there is not an exponential-state transition that preempts a GEN transition t 7 ).

6 Performance indices

The performance measures of our inventory systems and the total expected cost rate in the steady-state can be defined as follows:

• The mean inventory level, n l ,

  n l = Q v ≪ 0 Q N ≫ + ∑ j = 1 S j ∑ k = 0 N − 1 [ v ≪ 0 j k ≫ + v ≪ 1 j k ≫ ] .

• The expected number of customers in the orbit, n o ,

  n o = N v ≪ 0 Q N ≫ + ∑ k = 1 N k v ≪ 00 k ≫ + ∑ j = 1 S ∑ k = 1 N − 1 k [ v ≪ 0 j k ≫ + v ≪ 1 j k ≫ ] .

• The expected reorder level, n r ,

  n r = ∑ k = 0 N − 1 μ v ≪ 1 s + 1 k ≫ .

• The probability that server is busy, n B ,

  n B = ∑ j = 1 S ∑ k = 0 N − 1 v ≪ 1 j k ≫ .

• The effective search rate, n ES ,

  n ES = ∑ j = 2 S ∑ k = 1 N − 1 p μ v ≪ 1 j k ≫ .

• The expected number of successful retrials, n ESR ,

  n ESR = γ [ v ≪ 0 Q N ≫ + ∑ j = 1 S ∑ k = 1 N − 1 v ≪ 0 j k ≫ ] .

• The probability that inventory level is zero, p 0 ,

  p 0 = ∑ k = 0 N v ≪ 00 k ≫ .

• The probability that inventory level is greater than s , p s ,

  p s = v ≪ 0 Q N ≫ + ∑ j = s + 1 S ∑ k = 0 N − 1 [ v ≪ 0 j k ≫ + v ≪ 1 j k ≫ ] .

• The mean generation of primary calls, λ e ,

  λ e = ∑ k = 0 N − 1 ( N − k ) v ≪ 00 k ≫ + ∑ j = 1 S ∑ k = 0 N − 1 ( N − k ) v ≪ 0 j k ≫ + ∑ j = 1 S ∑ k = 0 N − 2 ( N − k − 1 ) v ≪ 1 j k ≫ .

• The mean waiting time, ω ,

  ω = n o λ e .

• The mean response time, ϖ ,

  ϖ = n o + n B λ e .

• The total expected cost rate T C in the steady-state is given by:

  T C = C h Q v ≪ 0 Q N ≫ + ∑ j = 1 S j ∑ k = 0 N − 1 [ v ≪ 0 j k ≫ + v ≪ 1 j k ≫ ] + C s ∑ k = 0 N − 1 μ v ≪ 1 s + 1 k ≫ + C w N v ≪ 0 Q N ≫ + ∑ k = 1 N k v ≪ 00 k ≫ + ∑ j = 1 S ∑ k = 1 N − 1 k [ v ≪ 0 j k ≫ + v ≪ 1 j k ≫ ] ,

  where C h , C s , and C w are, respectively, the inventory carrying cost per unit item per unit time, setup cost per order, and the waiting cost for an orbiting demand per unit time.

7 Numerical application

We construct two programs based on the formulas obtained in previous sections in order to compute numerically the probability distributions of the two models (GSPN and MRSPN). We show the influence of the system parameters on the system performance measures and the total expected cost rate. The numerical results for:

1. the GSPN model are given in Tables 3, 4, 5 and illustrated in Figures 3, 4, 5, and 6.

2. the MRSPN model is given in Table 6.

Figure 3

Effect of arrival rate λ and the search probability p on the mean response time ϖ ; N = 12 , S = 15 , s = 5 , λ = 0.1 , … , 1.5 , μ = 1.2 , γ = 0.15 , α = 1 , and p = 0 , … , 1 .

Figure 4

Effect of arrival rate λ and the retrial rate γ on the mean response time ϖ ; N = 16 , S = 14 , s = 6 , λ = 0.15 , … , 3 , μ = 2 , γ = 0.15 , … , 1.5 , α = 1 , and p = 0.4 .

Figure 5

Effect of arrival rate λ and the retrial rate γ on the total expected cost rate T C ; N = 20 , S = 25 , s = 8 , λ = 1 , … , 5 , μ = 3.5 , γ = 0.1 , … , 2 , α = 2 , p = 0.3 , C h = 2.8 , C s = 10 , and C w = 17 .

Figure 6

Effect of the replenishment rate α and the service time rate μ on the total expected cost rate T C ; N = 18 , S = 25 , s = 10 , λ = 1.5 , μ = 1 , … , 6 , γ = 0.5 , α = 0.5 , … , 4 , p = 0.6 , C h = 0.2 , C s = 7 , and C w = 16 .