Performance of a Synchronisation Station with Abandonment

Abstract

The paper presents a Markovian queueing model for assessing the performance of synchronisation between stations in a production system. The system at hand consists of K distinct buffers, each buffer storing an item that is needed for the next production stage. Departures are immediate when all items are present. Due to the presence of multiple buffers, there is no reasonably fast way to calculate the stationary distribution of the Markov chain. Therefore, we focused on the series expansion of the stationary distribution in terms of the arrival rate. We provide a fast algorithm for calculating these terms. Comparing our results with stochastic simulation, we show that the expansion approach converges to the simulated values for a wide range of arrival rates.

1. Introduction

The synchronisation of material flows is a key element in supply chain management and plays a pivotal role in ensuring efficiency in various industries including automotive [1] and electronics [2,3]. By coordinating different stages of a production process, resources such as machines, materials, and labour can be utilised more efficiently, which, in turn, leads to an increase in overall productivity. A synchronised production process allows for faster and more-predictable lead times. Nevertheless, due to the presence of uncertainty in production times at the different stages of the production process, some buffering may be required to streamline the process. For example, Borodin [4] and Romero-Silva and Hurtado-Hernàndez [5] note that the supply of components to assembly operations is subject to variability and uncertainty, as well as the task of assembling these components itself. In the absence of buffers, the subsequent production stages are prone to blocking, which reduces overall production efficiency. Such blocking is either caused by a lack of input materials or by the absence of storage at the output.

In this paper, we introduce a Markovian queueing model that captures the essential features and dynamics of synchronisation under uncertainty. The model at hand consists of multiple queues in parallel, each queue representing a distinct product inventory that is used to mitigate the blocking of the adjacent production stage. The different inventories are tied to distinct independent production stages, while the subsequent stage requires all items to be present. We, here, in particular, focus on the case where intermediate storage in the inventories is limited in time. The products cannot be used beyond a certain time limit, which is captured by abandonments in the different queues. Hence, the model at hand involves multiple queues with abandonments and synchronous departures from the different queues. From a mathematical perspective, the synchronisation system at hand, therefore, shares some properties with the well-investigated fork–join queues [6,7,8,9], as well as with queues with paired service [10,11,12]. In a fork–join queue, a job is split up between servers (forking) and, then, again, combined in a later stage (joining). Hence, in a fork–join queue, there is synchronisation between arrivals, which gives rise to a different type of dynamic. Queues with paired service, also referred to as assembly-like queues or kitting processes, also assume that the arrival processes are independent. However, departures are not immediate, as pairing requires service.

The synchronisation model at hand with only two buffers further corresponds to a so-called double-ended queue [13,14]. As service is immediate, one of the queues is empty, and the synchronisation station can be modelled by an integer-valued Markov chain. Positive states correspond to items in one queue, while negative states correspond to items in the other queue. Such double-ended queues are, e.g., motivated by taxi stands, where a trip starts once the taxi is present, as well as its passengers.

The main contribution of the present paper involves extending the analysis of the double-ended queue to Markovian models for the synchronisation of multiple streams. When the number of streams is limited and the provided buffer capacity per stream is limited, the corresponding Markov chain has a limited number of states, and the steady-state solution of the Markov chain can be readily obtained by solving the linear system of balance equations. As the state space of the synchronisation station is multidimensional, the model suffers from the so-called curse of dimensionality: the state space of the Markov chain grows exponentially in the number of dimensions, which makes the computational cost of directly solving the balance equations prohibitive for practical applications. A second contribution, therefore, lies in devising an efficient algorithm for numerically solving the Markovian synchronisation model. We, here, rely on a series expansion approach. Instead of calculating the stationary transition probabilities, the series expansion approach studies the stationary transition probabilities as functions of a system parameter (the arrival rate $\lambda$ in the present setting) and calculates the terms in the series expansion of these functions. At first sight, this does not mitigate the computational burden of the large state space, as there are as many equations for each term in the series expansions as there are equations in the balance equations. However, if the chain at $\lambda =0$ only allows for transitions in one direction (either upwards or downwards), the transition matrix is a triangular matrix, and the computational complexity for a chain with M states is reduced to at most $O\left(M^2\right)$. Moreover, in many practical applications, the number of possible transitions from a state is much smaller than the size of the state space, so that the numerical complexity for the calculation of each term is $O\left(M\right)$ and the complexity for the calculation of the first N terms of the expansion is only $O\left(MN\right)$.

The presented series expansion technique is known by various terms in the literature, including Markov chain perturbation, light-traffic approximation, and the power series algorithm. We refer to the surveys on series expansion techniques [15,16] and the references therein. Despite the lack of a fixed nomenclature, perturbation methods primarily originate from the necessity to analyse how a system’s performance is affected by changes in specific parameters. Singular perturbations, which disrupt the class-based structure of the unperturbed chain, have been a significant focus in research, as seen in works such as [17,18,19]. When the perturbed parameter relates to the arrival rate of a queue, the series expansion is a light-traffic approximation. Lastly, the power series algorithm transforms a Markov chain of interest into a collection of Markov chains indexed by a potentially artificially introduced variable $\rho$. This transformation is chosen to allow for a straightforward computation of perturbations around $\rho =0$, while as $\rho$ approaches 1, the original Markov chain is regained. Hence, the series expansion can be used for approximating solutions of the original Markov chain. This holds true if the radius of convergence of the series expansion exceeds $\rho =1$ [20,21,22,23]. The methodology in the present setting is a light-traffic approximation, though the numerical results reveal that the results are accurate for a wide range of arrival rates. The approach is primarily algorithmic, with a focus on a recursion for calculating many terms in the series expansion, rather than on closed-form expressions for the first few terms. Such an algorithmic approach previously allowed for analysing queueing systems with coupled service [11,12], retrial queues with heterogeneous retrials [24], and epidemic processes in finite populations [25]. In order to improve convergence, we further rely on convergence acceleration techniques. Such techniques not only accelerate the convergence in the region of convergence, but also extend the region of convergence [26,27,28]. We, here, in particular, apply Wynn’s $ϵ$-algorithm to improve convergence. Wynn’s $ϵ$-algorithm is an efficient implementation of Shank’s transformation, which is the best all-purpose acceleration method according to [29]. Under non-stringent assumptions on the sequences, the transformation converges faster to the limit of the original series in the region of convergence and converges to the analytic continuation of diverging sequences [30].

The remainder of this paper is organised as follows. In the next section, we introduce the Markovian modelling assumptions and find the corresponding set of balance equations. We study the synchronisation station in the “light-traffic” regime in Section 3. In contrast to a direct solution, we show that the computational cost of calculating the terms in the light traffic approximation is considerably lower. We, then, extend the modelling assumptions by the inclusion of a Markovian environment variable, which modulates the arrival and abandonment rates in Section 4. We illustrate our approach by a set of numerical experiments in Section 5. Finally, the conclusions are drawn in Section 6.

2. Mathematical Model

We considered a synchronisation station with N finite or infinite capacity queues, each queue storing a particular part. Parts arrive at the nth queue in accordance with a Poisson process with rate ${\lambda }_n>0$, $n\in \{1,\dots ,N\}$. Let $C_n\in \mathbb{N}\setminus \left\{0\right\}\cup \{\infty \}$ denote the (possibly infinite) capacity of the nth queue. Whenever all queues are non-empty, a container with all parts immediately leaves the synchronisation station. Note that, at any point in time, at least one of the queues will be empty. If this is not the case, a container is ready and departs immediately. Apart from container departures, parts may also leave prior to being packed in a container. To be more precise, each part has a finite maximal residence time, which is exponentially distributed. Let $1/{\alpha }_n$ denote the average maximal residence time of parts in the nth queue. In other words, each part in the nth queue abandons the queue with rate ${\alpha }_n>0$.

Considering the modelling assumptions outlined above, the Markovian state of the synchronisation station is completely determined by the number of parts in each queue. Let $X_n\left(t\right)\in \{0,\dots ,C_k\}\doteq {\mathcal{C}}_k$ denote the number of parts in the nth queue at time t, and let $\mathit{X}\left(t\right)=[X_1\left(t\right),\dots ,X_N\left(t\right)]$ denote the vector of queue contents at time t. As one of the queues is empty for sure, the state space of the Markov process $\mathit{X}\left(t\right)$ is given by

$$\mathcal{C}=({\mathcal{C}}_1\times \dots C_n)\setminus ({\mathcal{C}}_1^0\times \dots \times {\mathcal{C}}_N^0),$$

with ${\mathcal{C}}_k^0={\mathcal{C}}_k\setminus \left\{0\right\}=\{1,\dots ,C_k\}$. Moreover, it is easy to verify that this Markov process is ergodic. First, the Markov chain is irreducible, as state $\mathit{X}=\mathbf{0}$ can be reached by abandonments from any state and any state can be reached from state $\mathbf{0}$ by arrivals. Secondly, comparing the system at hand with the (ergodic) system without container departures ensures that state zero is recurrent. Indeed, in the absence of container departures, departures in the queues stem from abandonments, and each queue can be modelled by an infinite server queue.

We focus on the stationary distribution of the Markov chain. For $\mathit{x}\in \mathcal{C}$, let $\pi \left(\mathit{x}\right)$ denote the stationary probability to be in state $\mathit{x}=[x_1,x_2,\dots ,x_N]$:

$$\pi \left(\mathit{x}\right)=\underset{t\to \infty }{lim}\mathsf{P}[\mathit{X}\left(t\right)=\mathit{x}].$$

Moreover, to simplify the notation, it is convenient to set $\pi \left(\mathit{x}\right)=0$ for $\mathit{x}\ne \mathcal{C}$. State transitions are either triggered by part abandonments or by new part arrivals. In state $\mathit{x}$, parts in the nth queue are abandoned with rate $x_n{\alpha }_n$, the new state being $\mathit{x}-{\mathit{e}}_n$. Here, ${\mathit{e}}_n$ is a vector of all zeros, apart from its nth element, which is one. For $x_n<C_n$, a new part arrives at the nth queue with rate ${\lambda }_n$. If ${\mathit{x}}_n+{\mathit{e}}_n\in \mathcal{C}$, the arrival does not complete a container and the state after the arrival is ${\mathit{x}}_n+{\mathit{e}}_n$. If ${\mathit{x}}_n+{\mathit{e}}_n\notin \mathcal{C}$, the new arrival completes a container, such that there is an immediate departure of this container. Therefore, the state after the arrival is ${\mathit{x}}_n+{\mathit{e}}_n-\mathit{e}$, where $\mathit{e}$ is a vector of ones. In view of these observations, we find the following set of balance equations:

$$\begin{array}{c}\pi \left(\mathit{x}\right)\sum _{n=1}^N({\lambda }_n{\mathsf{1}}_{\{x_n<C_n\}}+{\alpha }_n{x}_n)\hfill \\ \hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}=\sum _{n=1}^N\pi (\mathit{x}+{\mathit{e}}_n){\alpha }_n(x_n+1)+\pi (\mathit{x}-{\mathit{e}}_n){\lambda }_n+\pi (\mathit{x}-{\mathit{e}}_n+\mathit{e}){\lambda }_n,\end{array}$$

(1)

for $\mathit{x}\in \mathcal{C}$. Here, ${\mathsf{1}}_{\{·\}}$ is the indicator function, which evaluates to 1 if its argument is true and to 0 if this is not the case. Note that the expression above holds both for finite and infinite $C_n$. For infinite $C_n$, the condition $x_n<C_n$ trivially holds for all $x_n$. Also, the assumption $\pi \left(\mathit{x}\right)=0$ for $\mathit{x}\notin \mathcal{C}$ greatly simplifies the balance equations as we can omit the indicators that the states $\mathit{x}+{\mathit{e}}_n$, $\mathit{x}-{\mathit{e}}_n$ and $\mathit{x}-{\mathit{e}}_n+\mathit{e}$ are part of the state space.

3. Light-Traffic Analysis

For small N and $C_n$, the balance Equation (1) is readily solved by direct methods like, e.g., Gaussian elimination. However, already for moderate N and $C_n$, the size of the state space of the Markov process makes a direct calculation of the stationary probabilities computationally demanding. Therefore, we focus on the synchronisation system in the light-traffic regime. Our approach can approximate the exact distribution with high accuracy and low computational effort, whereby the computational effort for the calculation of the Tth term in the series expansion is $O(Nmin{(T,C)}^N)$ with $C={max}_k{C}_k$, cf. infra.

To this end, we set ${\lambda }_n={\kappa }_n\lambda$ and introduce the Maclaurin series expansion in $\lambda$ of the solution of the balance equations:

$$\pi \left(\mathit{x}\right)=\sum _{k=0}^{\infty }{\pi }_k\left(\mathit{x}\right){\lambda }^k.$$

Before proceeding, we note that, for finite $C_k$, it is easy to check—by applying Cramer’s rule, that one explicitly expresses $\pi \left(\mathit{x}\right)$ as a fraction of two determinants; that the stationary probabilities are rational functions of $\lambda$. Since rational functions have no singularities other than poles in the extended complex plane and there is no pole at $\lambda =0$, we conclude that the stationary probabilities are analytic functions in a region around $\lambda =0$. This, then, justifies the series expansion. However, analyticity is not necessarily required. In the following numerical examples, we rely on the $ϵ$ method for convergence acceleration. This method is known to extend the region of convergence and even works when the radius of convergence is 0 [30].

Substituting the series expansions above in the balance equations and isolating the terms in ${\lambda }^k$ leads to

$$\begin{array}{c}{\pi }_k\left(\mathit{x}\right)={\left(\sum _{n=1}^N{\alpha }_n{x}_n\right)}^{-1}(-{\pi }_{k-1}\left(\mathit{x}\right)\sum _{n=1}^N{\kappa }_n{\mathsf{1}}_{\{x_n<C\}}\hfill \\ \hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+\sum _{n=1}^N{\pi }_k(\mathit{x}+{\mathit{e}}_n){\alpha }_n(x_n+1)+{\pi }_{k-1}(\mathit{x}-{\mathit{e}}_n){\kappa }_n+{\pi }_{k-1}(\mathit{x}-{\mathit{e}}_n+\mathit{e}){\kappa }_n),\end{array}$$

(2)

for $k\in \{1,2,\dots \}$ and for $\mathit{x}\in \mathcal{C}\setminus \left\{\mathbf{0}\right\}$. For $\mathit{x}=\mathbf{0}$, the normalisation condition:

$$\sum _{\mathit{x}\in \mathcal{C}}\pi \left(\mathit{x}\right)=\sum _{\mathit{x}\in \mathcal{C}}\sum _{k=0}^{\infty }{\pi }_k\left(\mathit{x}\right){\lambda }^k=1$$

yields

$${\pi }_k\left(\mathbf{0}\right)=-\sum _{\mathit{x}\in \mathcal{C}\setminus \left\{\mathbf{0}\right\}}{\pi }_k\left(\mathit{x}\right),$$

(3)

$k\in \{1,2,\dots \}$. Moreover, for $k=0$, one easily finds

$${\pi }_0\left(\mathit{x}\right)={\mathsf{1}}_{\{\mathit{x}=\mathbf{0}\}}.$$

The latter result is also intuitively clear. In the absence of part arrivals, all part queues are empty.

In contrast to the balance Equation (1), the system of Equations (2) and (3) is easily solved. Indeed, close inspection of Equations (2) and (3) reveals that one can solve these equations in reverse lexicographical order: to calculate ${\pi }_k\left(\mathit{x}\right)$ by (2) and (3), one only needs $(k-1)$ the-order terms and terms in ${\pi }_k\left(\mathit{y}\right)$ for $\mathit{y}$ lexicographically larger than $\mathit{x}$. In addition, by the so-called k-events rule, one can show that ${\pi }_k\left(\mathit{x}\right)=0$ if the state $\mathit{x}$ cannot be reached from state $\mathbf{0}$ with k arrivals (and possibly any number of other transitions). In this case, we have ${\pi }_k\left(\mathit{x}\right)=0$ if the total queue size exceeds k: ${\sum }_n{x}_n>k$. Summarising, we find that one needs to calculate at most $O(min{(T,C)}^N)$ unknowns for the Tth-order terms, each term requiring $O\left(N\right)$ additions. Hence, the calculation of the Tth-order term is of order $O(Nmin{(T,C)}^N)$. Note that the factor $min(T,C)$ follows from the n-events rule used above. We have ${\pi }_T\left(\mathit{x}\right)=0$ and state that it cannot be reached from the empty queueing state with T arrivals.

Once the terms in the series expansions of the stationary probabilities have been calculated, various performance measures of interest can be studied. The Kth order expansion of the ℓth moment of the nth queue equals

$$\mathsf{E}\left[X_n^{𝓁}\right]\approx \sum _{k=0}^K\left(\sum _{\mathit{x}\in \mathcal{C}}{\pi }_k\left(\mathit{x}\right)x_n^{𝓁}\right){\lambda }^k.$$

(4)

Similar expressions are readily obtained for the cross-moments of the different queue content as well.

The Kth-order expansion of the abandonment rate ${\nu }_n$ and loss probability ${\eta }_n$ in the nth queue further equal

$${\nu }_n\approx \sum _{k=0}^K\left(\sum _{\mathit{x}\in \mathcal{C}}{\pi }_k\left(\mathit{x}\right)x_n{\alpha }_n\right){\lambda }^k,{\eta }_n\approx \sum _{k=0}^K\left(\sum _{\mathit{x}\in \mathcal{C}}{\pi }_k\left(\mathit{x}\right){\mathsf{1}}_{\{x_n=C_n\}}\right){\lambda }^k.$$

Indeed, each of the $x_n$ parts in the nth queue is abandoned with rate ${\alpha }_n$, while there is a loss if there is an arrival in a full queue. Finally, the Kth-order expansion of the (simultaneous) departure rate $\mu$ from the queues equals

$$\mu \approx \sum _{k=0}^K\left(\sum _{\mathit{x}\in \mathcal{C}}{\pi }_k\left(\mathit{x}\right)\sum _{n=1}^N{\kappa }_n{\mathsf{1}}_{\{x_n=0,x_m>0,\forall m\ne n\}}\right){\lambda }^{k+1}.$$

Here, we used the observation that we have a departure if there is an arrival in an empty queue, and the other queues or not empty.

Remark 1.

The expressions of the performance measures show that the terms in the series expansions of the performance measures can be calculated once the corresponding terms of the stationary probabilities are obtained. This observation allows for a considerable reduction of the memory requirements of the algorithm. Indeed, while close inspection of (2) shows that the nth terms need to be preserved during the calculation of the $(n+1)$th terms, once these calculations are complete and once the nth terms of the performance measures are calculated, the nth terms are no longer required, and the storage can be freed or overwritten. Summarising, the storage requirement is $O(Nmin{(T,C)}^N)$, which only depends on the number of terms by the truncation of the state space, which follows from the n-events rule (cf. supra).

4. Markovian Environment

We now extend our findings to a synchronisation station, where the arrival rates in the different queues, as well as the abandonment rates from the different queues depend on an exogenous Markov process with a finite state space. More precisely, let $Y\left(t\right)$ be an ergodic Markov process with finite state space $\mathcal{M}$ and with generator matrix $B={\left[{\beta }_{ij}\right]}_{i,j\in \mathcal{M}}$. Here, ${\beta }_{ij}$ denotes the transition rate from state i to state j ($i,j\in \mathcal{M}$ with $i\ne j$), while the diagonal elements are equal to

$${\beta }_{ii}\doteq -\sum _{j\in \mathcal{M}\setminus \left\{i\right\}}{\beta }_{ij},$$

as usual (for $i\in \mathcal{M}$). For further use, let $\mathit{\nu }$ denote the invariant distribution of this Markov process. That is, $\mathit{\nu }$ solves $\mathit{\nu }B=0$ and $\mathit{\nu }{\mathit{e}}^{\prime }=1$ or, equivalently,

$$\mathit{\nu }={(B+{\mathit{e}}^{\prime }\mathit{e})}^{-1}\mathit{e}.$$

When the exogenous state is $j\in \mathcal{M}$, new parts arrive at the nth queue with rate ${\lambda }_n^j>0$, while parts in the nth queue are abandoned with rate ${\alpha }_n^j>0$. Retaining the notation of the preceding section, it is now easily verified that the process $\left(\mathit{X}\right(t),Y(t\left)\right)$ constitutes an ergodic Markov process with state space $\mathcal{C}\times \mathcal{M}$. For $(\mathit{x},j)\in \mathcal{C}\times \mathcal{M}$, let

$$\pi (\mathit{x},j)=\underset{t\to \infty }{lim}\mathsf{P}[\mathit{X}\left(t\right)=\mathit{x},Y\left(t\right)=j]$$

denote the limiting distribution of the Markov process.

Remark 2.

For the ease of exposition, we assumed positive arrival and abandonment rates for all states of the environment. This is not strictly required. In the finite setting, a positive arrival and abandonment rate in one of the states suffices. In the infinite setting, it suffices to check that each queue is ergodic in the absence of container departures.

Assuming $\pi (\mathit{x},j)=0$ for $(\mathit{x},j)\notin \mathcal{C}\times \mathcal{M}$, the balance equations then read

$$\begin{array}{c}\pi (\mathit{x},j)\left(\sum _{n=1}^N({\lambda }_n^j{\mathsf{1}}_{\{x_n<C_n\}}+{\alpha }_n^j{x}_n)+\sum _{k\in \mathcal{M}\setminus \left\{j\right\}}{\beta }_{jk}\right)=\sum _{i\in \mathcal{M}\setminus \left\{j\right\}}\pi (\mathit{x},i){\beta }_{ij}\hfill \\ \hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+\sum _{n=1}^N\left(\pi (\mathit{x}+{\mathit{e}}_n,j){\alpha }_n^j(x_n+1)+\pi (\mathit{x}-{\mathit{e}}_n,j){\lambda }_n^j+\pi (\mathit{x}-{\mathit{e}}_n+\mathit{e},j){\lambda }_n^j\right),\end{array}$$

(5)

for $\mathit{x}\in \mathcal{C}$ and $j\in \mathcal{M}$. For convenience, let $\mathit{\pi }\left(\mathit{x}\right)={\left[\pi (\mathit{x},j)\right]}_{j\in \mathcal{M}}$ denote the row vector that collects the stationary probabilities for a fixed $\mathit{x}\in \mathcal{C}$. Moreover, let $A_n$ and ${\Lambda }_n$ denote diagonal matrices with jth diagonal entry ${\alpha }_n^j$ and ${\lambda }_n^j$, respectively. The balance equations can then be rewritten as follows

$$\begin{array}{c}\mathit{\pi }\left(\mathit{x}\right)\sum _{n=1}^N({\Lambda }_n{\mathsf{1}}_{\{x_n<C_n\}}+A_n{x}_n)\hfill \\ \hfill \hspace{1em}\hspace{1em}\hspace{1em}=\mathit{\pi }\left(\mathit{x}\right)B+\sum _{n=1}^N\left(\mathit{\pi }(\mathit{x}+{\mathit{e}}_n)A_n(x_n+1)+\mathit{\pi }(\mathit{x}-{\mathit{e}}_n){\Lambda }_n+\mathit{\pi }(\mathit{x}-{\mathit{e}}_n+\mathit{e}){\Lambda }_n\right).\end{array}$$

(6)

We now express ${\Lambda }_k$ in terms of a single rate $\lambda$, ${\Lambda }_k=\lambda K_k$, and as for the synchronisation station without an exogenous environment, we introduce the Maclaurin series expansion in $\lambda$ of the solution of the balance equations:

$$\mathit{\pi }\left(\mathit{x}\right)=\sum _{k=0}^{\infty }{\mathit{\pi }}_k\left(\mathit{x}\right){\lambda }^k.$$

Plugging the series expansion above in (4) and comparing terms in equal powers of $\lambda$ then yield

$${\mathit{\pi }}_0\left(\mathit{x}\right)=\mathbf{0},$$

(7)

and

$$\begin{array}{c}{\mathit{\pi }}_k\left(\mathit{x}\right)=(-{\mathit{\pi }}_{k-1}\left(\mathit{x}\right)\sum _{n=1}^N{K}_n{\mathsf{1}}_{\{x_n<C_n\}}+\sum _{n=1}^N{\mathit{\pi }}_k(\mathit{x}+{\mathit{e}}_n)A_n(x_n+1)\hfill \\ \hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+\sum _{n=1}^N{\mathit{\pi }}_{k-1}(\mathit{x}-{\mathit{e}}_n)K_n+{\mathit{\pi }}_{k-1}(\mathit{x}-{\mathit{e}}_n+\mathit{e})K_n){\left(\sum _{n=1}^N{A}_n{x}_n-B\right)}^{-1},\end{array}$$

(8)

for $k=1,2,\dots$ and $\mathit{x}\in \mathcal{C}\setminus \left\{\mathbf{0}\right\}$. Note that Equation (7) simply states that there are no parts in any of the queues in absence of arrivals.

For $\mathit{x}=\mathbf{0}$, we further find

$$\begin{array}{c}{\mathit{\pi }}_0\left(\mathbf{0}\right)B=\mathbf{0},\hfill \\ {\mathit{\pi }}_k\left(\mathbf{0}\right)B={\mathit{\pi }}_{k-1}\left(\mathbf{0}\right)\sum _{n=1}^N{K}_n-\sum _{n=1}^N\left({\mathit{\pi }}_k\left({\mathit{e}}_n\right)A_n+{\mathit{\pi }}_{k-1}(\mathit{e}-{\mathit{e}}_n)K_n\right).\hfill \end{array}$$

(9)

Combining these expressions with the normalisation conditions:

$${\mathit{\pi }}_0\left(\mathbf{0}\right){\mathit{e}}^{\prime }=1,\sum _{\mathit{x}\in \mathcal{C}}{\mathit{\pi }}_k\left(\mathit{x}\right){\mathit{e}}^{\prime }=0,$$

then yields ${\mathit{\pi }}_0\left(\mathbf{0}\right)=\mathit{\nu }$ and

$${\mathit{\pi }}_k\left(\mathbf{0}\right)=\left({\mathit{\pi }}_{k-1}\left(\mathbf{0}\right)\sum _{n=1}^N{K}_n-\sum _{n=1}^N\left({\mathit{\pi }}_k\left({\mathit{e}}_n\right)A_n+{\mathit{\pi }}_{k-1}(\mathit{e}-{\mathit{e}}_n)K_n\right)+\kappa \mathit{e}\right){(B+{\mathit{e}}^{\prime }\mathit{e})}^{-1},$$

with

$$\kappa ={\mathit{\pi }}_k\left(\mathbf{0}\right){\mathit{e}}^{\prime }=-\sum _{\mathit{x}\in \mathcal{C}\setminus \left\{\mathbf{0}\right\}}{\mathit{\pi }}_k\left(\mathit{x}\right){\mathit{e}}^{\prime }.$$

Note that we cannot simply obtain ${\mathit{\pi }}_k\left(\mathbf{0}\right)$ from Equation (9) as B is not invertible.

Again, the terms in the series expansion can be easily determined in reverse lexicographical order, where, by the k-events rule, one can show that ${\mathit{\pi }}_k\left(\mathit{x}\right)=\mathbf{0}$ if the state $\mathit{x}$ cannot be reached from state $\mathbf{0}$ with k arrivals (and possibly any number of other transitions). The scalar operations of the preceding section are now replaced by matrix operations. Hence, following the same reasoning as in the preceding section, we find that the numerical complexity for calculating the Tth order term is $O(Nmin(T,C)M^3)$, where M denotes the number of environment states. Once the series expansion of the stationary vectors $\mathit{\pi }\left(\mathit{x}\right)$ is obtained, series expansions for various performance measures of interest are readily obtained. The expression (4) for the ℓth moment of the queue size remains valid, with $\pi \left(\mathit{x}\right)={\sum }_{j\in \mathcal{M}}\pi (\mathit{x},j)$. The Kth order expansion of the departure rate now reads

$$\mu \approx \sum _{k=0}^K\left(\sum _{\mathit{x}\in \mathcal{C}}\sum _{j\in \mathcal{M}}{\pi }_k\left(\mathit{x}\right)\sum _{n=1}^N{\kappa }_n^j{\mathsf{1}}_{\{x_n=0,x_m>0,\forall m\ne n\}}\right){\lambda }^{k+1}.$$

5. Numerical Examples

We now illustrate our numerical approach by a set of examples. To assess the accuracy of the series approximations, we compare the results that were obtained by the analytic approach with stochastic simulation. We first demonstrate the accuracy for the synchronisation system without the random environment and, then, extend our examples by including an environment variable that allows for including burstiness in the arrival processes of one of the queues.

Preliminary experimentation showed that the series expansions only converge to the correct value for small $\lambda$ values, which is an indication that the region of convergence of the series expansion is limited. Luckily, convergence acceleration techniques allow for extending the region of convergence [26,27,28]. We, here, implement Wynn’s $ϵ$-algorithm to speed up convergence, as well as to extend the region of convergence. To this end, for a performance measure of interest, introduce the sequence of partial sums $S_n={\sum }_{k=0}^n{s}_k{\lambda }^k$. Here, the sequence $s_k$ corresponds to the coefficients in the expansion of the performance measure. Expressions for these terms are given in Section 3 and are readily calculated once the terms in the expansions of the stationary probabilities are determined. The $ϵ$-algorithm then approximates the limit $n\to \infty$ by ${ϵ}_k^{\left(0\right)}$, for k even (odd indices are intermediate values), recursively defined as follows:

$$\begin{array}{cc}\hfill {ϵ}_{-1}^{\left(n\right)}& =0,{ϵ}_0^{\left(n\right)}=S_n,\hfill \\ \hfill {ϵ}_{r+1}^{\left(m\right)}& ={ϵ}_{r-1}^{(m+1)}+\frac{1}{{ϵ}_r^{(m+1)}-{ϵ}_r^{\left(m\right)}},\hfill \end{array}$$

for $k,n\in \{0,1,2,\dots \}$. Our implementation includes a simple protection against division by zero, as well as the simple stopping rule of [31] (p. 213).

5.1. Basic Model

We first consider a synchronisation station with five queues, with equal arrival rates $\lambda$ in all queues and with capacities $C_1=C_2=10$, $C_3=15$, $C_4=20$, and $C_5=25$. The abandonment rates in the different queues equal ${\alpha }_1=0.25\gamma$, ${\alpha }_2=0.225\gamma$, ${\alpha }_3=0.2\gamma$, ${\alpha }_4=0.175\gamma$, and ${\alpha }_5=0.15\gamma$. Here, $\gamma$ is a scale factor for the abandonment rates. This factor can be used to study how the abandonment rates affect performance, with higher $\gamma$ corresponding to more abandonments. With these queue capacities, there are approximately ${10}^6$ states. Figure 1a,b depict the mean and variance of the queue sizes versus the arrival rate $\lambda$, respectively. The scale parameter is fixed to $\gamma =1$ for all curves. The different curves are calculated by applying the $ϵ$-algorithm on the 30th-order series expansion. Not all terms in the series expansion are necessarily used though, the expansion order being an upper bound for the number of terms in the stopping rule of [31]. The markers were obtained by means of stochastic simulation and are included to verify the accuracy of the series expansion approach. The simulation results were obtained by means of the replication–deletion approach. We simulated 100 replications of 5 × ${10}^6$ state transitions. The obtained $99\%$ confidence intervals are smaller than $1\%$ of the simulated value and are, therefore, omitted from the plots. The figures show that the correspondence between the simulation results and the series expansion approach is excellent for the complete range of the $\lambda$-values. We would like to point out that the $ϵ$-algorithm is key for the accuracy of the results. Without convergence acceleration, the region of convergence of the series expansion is limited, and the results are only accurate for very small arrival rates ($\lambda <0.05$). We excluded the curves without acceleration from the figures to enhance clarity.

To evaluate the effect of the abandonment rates on performance, Figure 2 depicts the mean total queue content and the departure rate from the queues vs. the arrival rate $\lambda$ for different values of the scale factor $\gamma$, as indicated. Again, values from the (30th-order) series expansion approach were compared with the equivalent values obtained by stochastic simulation. As in Figure 1, the simulation results were obtained by means of the replication–deletion approach with 100 replications of 5 × ${10}^6$ state transitions. The simulation study showed that the combination of the series expansion approach and the $ϵ$-algorithm leads to very accurate results. More arrivals translate into larger queue sizes, while more abandonments imply lower queue sizes. Analogously, the departure rate increases with the arrival rate and decreases with the abandonment rate. Recall that the departure rate only includes departures by synchronisation and excludes departures by abandonments. Ideally, the departure rate should equal the arrival rate, so that all arrivals can leave the queues by synchronisation. There are two impediments to such a perfect operation of the synchronisation station. First, some part arrivals may not enter the synchronisation station due to a lack of queue capacity. Secondly, parts may be abandoned, prior to synchronisation. The latter effect is dominant in the present example, especially for higher abandonment rates (that is, for higher $\gamma$).

5.2. Markovian Environment

So far, the queueing model did not include the Markovian environment variable of Section 4. We now consider a model with 10 queues, all with capacity $C_n=2$ and all with abandonment rate ${\alpha }_n=0.1$, $n\in \{1,\dots ,N\}$. The arrivals in all but the first queue stem from independent Poisson processes with common rate ${\lambda }_n=1$. In contrast, the part arrivals in the first queue stem from an interrupted Poisson process. More precisely, we introduce a random environment with two states, say on and off, and there are only arrivals in the on state. Let a denote the transition rate from the on to the off state, and let b denote the transition rate from the off state to the on state. The fraction of time $\sigma$ the Markov environment is in the on state and the length $\kappa$ of an on and off period then equal

$$\sigma =\frac{b}{a+b},\kappa =\frac{1}{a}+\frac{1}{b}.$$

The pair $(\sigma ,\kappa )$ uniquely determines the transition rates a and b. We further assume that the (Poisson) arrival rate in the on state equals ${\widehat{\lambda }}_1=1/\sigma$, such that the average arrival rate in the first queue is 1 as well. One expects that smaller $\sigma$ and larger $\kappa$ affect performance negatively. For smaller $\sigma$, the arrivals are more concentrated, while for larger $\kappa$, the periods with and without arrivals are longer. This means that the queue can build up for a longer time.

Figure 3a shows the mean queue size as a function of the arrival rate for different values of $\sigma$ and fixed $\kappa =20$. Only the mean queue size of the first two queues is shown. Due to symmetry, the mean queue size of the remaining queues is identical to the mean queue size of the second queue. The curves were obtained by the 80th-order series expansion approach, complemented with the $ϵ$-algorithm. The markers again represent stochastic simulation results. The simulation results were obtained by means of the replication–deletion approach with 100 replications of ${10}^9$ state transitions. The upper and lower bounds of the obtained $99\%$ confidence intervals are visually indiscernible and, therefore, omitted. The correspondence between the analytic and simulation results is again good, though, by the additional complexity of the environment, a higher order expansion is needed to obtain accurate results. Recall that the arrivals in queue 1 are more concentrated in time for lower values of $\sigma$. This leads to a lower mean queue size for the first queue and a higher mean queue size for the second queue. The concentration of arrivals leads to considerable loss in the first queue, followed by periods without arrivals. The loss explains why there are fewer parts in this queue on average. On the other hand, the absence of parts in the first queue implies that there are fewer synchronised departures, which, in turn, means that there are more parts in the other queues. This latter observation is confirmed by Figure 3b, which shows the departure rate $\mu$ for different $\kappa$ and $\sigma$. Both an increase of $\kappa$ and a decrease $\sigma$ lead to a decrease of the departure rate. Also, Figure 3b compares the analytic results (80th-order expansion) with the simulation results. In contrast to the preceding figures, the accuracy is not as good, especially for $\kappa =40$ and $\sigma =0.25$, which is the example where the arrivals are concentrated over long periods, followed by long periods without arrivals.

To conclude this numerical investigation, we briefly summarise some observations on the applicability and accuracy of the series expansion approach. Foremost, we found that the method works well for a large number of parameter sets, though the number of terms in the expansion that is required to obtain accurate results does depend on the specific parameter values. In general, there is no direct relation between the size of the state space and the required order. We, however, observed that the presence of the Markovian environment most often required higher-order expansions. As the calculation of the series expansion is linear in the order of the expansion, computational limitations did not prevent the use of higher-order expansions. However, we found that Wynn’s algorithm may encounter numerical instability when dealing with higher-order series expansions for specific $\lambda$-values, so there is a practical limit on the order.

6. Conclusions

In this study, we introduced a Markovian queueing model to evaluate the performance of synchronisation among stations within a complex production system featuring K distinct buffers. These buffers store parts necessary for advancing to the subsequent production stages, and departures occur instantly upon the availability of all required items. The intricacies arising from multiple buffers render the calculation of the stationary distribution of the Markov chain a challenging task. To address this challenge, we applied a series expansion approach where the arrival rate was used as the parameter of the expansion.

Our primary contribution lies in the development of a rapid algorithm tailored to the computation of terms within the series expansion. By employing this algorithm, we mitigated the computational challenges associated with the calculation of the stationary distribution. Our numerical experiments confirmed that the proposed algorithm enables efficient and accurate computation of the stationary distribution. By comparing with stochastic simulation results, the accuracy was observed consistently across a broad spectrum of arrival rates, highlighting the robustness and versatility of our methodology in capturing the dynamic behaviour of the production system.