Transient Response of Homogenous and Nonhomogenous Bernoulli Production Lines

Abstract

The transient response of production systems is of significant importance especially if present advancements in Digital Twinning technology are taken into account. While the steady-state response enables long-term strategic decision making, the transient response enables more detailed simulation concerning aspects like production losses and preventive maintenance. This is especially relevant if nonhomogenous aspects of production systems are taken into account. An analytical and approximative solution to the problem of the transient response of homogenous and nonhomogenous Bernoulli production systems is developed in this paper based on the eigendecomposition of transition matrices, the eigenvalue problem, and the finite-state method. In particular, sub-resonant and resonant nonhomogeneous production lines are introduced for the first time. Also, the most significant key performance indicators are developed as functions of the time elapsed from the first cycle. Finally, the relationship between the number of eigenvalues and the accuracy of the results is inspected by employing a sensitivity analysis. The presented theoretical framework was employed in the case of a wood processing facility to present the potential application of the theory in the case of long- and short-term management of production systems.

1. Introduction

Production systems play a central role in the modern manufacturing of goods and products distributed all over the world. In an economic sense, manufacturing stands for a path to creating and developing new jobs, incomes, and sophisticated supply chains. It has significant potential to tackle environmental and social issues by introducing sustainable, energy-efficient, and cost-effective production concepts as well as suitable incomes for workers of different qualifications [1,2]. Also, production systems are often recognized as the driving force of the Digital Twinning concept, especially since it was first introduced in the case of manufacturing systems and associated lifecycle management of products [3]. However, these concepts would remain a mere list of hypothetical ideas if sophisticated mathematical models were not employed. Indeed, the mathematical background of production system engineering and the associated modeling approaches are the key enablers of frameworks such as Digital Twinning, factories for the future, Overall Equipment Efficiency (OEE), Lean production, Industry 4.0, the Internet of Things, Augmented and Virtual Reality, and others [4].

The existing mathematical models employ different approaches to evaluate the key performance indicators (KPIs) of production systems and address design issues. Some of them include Markov processes, queuing models, stochastic automata networks, Petri nets, and process and diodic algebra used to evaluate configurations such as serial, assembly, splitting, re-entrant, and rework lines as well as job shops, flexible manufacturing cells, or other specific systems [5,6,7]. Amongst these, Markov chains and simulations based on queuing models are usually addressed in the majority of publications and applications. Hence, potentially significant financial benefits, a wide application range of production system engineering, and the large community of scientists and professionals employing it in a variety of cases justify the constant requirement for deepening our understanding and knowledge of mathematical modeling, especially concerning Markov chains as one of the most representative cases.

The application of the Markovian framework in the case of production systems was considered in 1962 for the first time by Sevast’yanov [8] in the case of a rather simple production system composed of two machines and one buffer placed between them. At that time, a quite complex analytical solution to the steady-state response of the system based on integral equations was formulated, leading to a conclusion that further generalizations were simply elusive. This fact motivated intensive research on approximative approaches to the problem, including more realistic arrangements of production systems with more than two machines. First, the decomposition method was presented by Gershwin [9] as an approximative procedure to evaluate the key performance indicators of serial lines. Then, an asymptotic analysis technique (or the aggregation procedure) was introduced by Lim et al. [10] and successfully applied later on in many cases, including primarily the steady-state response of serial, splitting, or assembly production lines composed of machines of Bernoulli, exponential, or non-exponential reliability [11]. Finally, an analytical solution to the problem has been developed in [12], followed by the newly introduced finite-state method, e.g., [13,14]. All three approximative methodologies, the decomposition method, the aggregation procedure, and the finite-state method, were successfully applied in the steady-state case of production systems.

The transient response has been seldom reported in the present literature. A few examples include only the application of the aggregation procedure in the case of Bernoulli production lines [15]. The main reasons for that can be found in the demanding theoretical framework as well as in the computational burden associated with the transient analysis of large-scale and transition-rich stochastic systems leading to a lack of analytical results that could be used to validate approximative approaches. Nevertheless, the importance of the production system’s transient response is reflected in the ability to evaluate the production losses due to the warm-up period and to schedule the production in a way that minimizes such losses [16]. Also, such an approach is a prerequisite to the development and support of preventive maintenance planning and raw material management as it enables the analysis of nonhomogenous (reliability of machines is a function of time) production systems [17,18]. Therefore, transient analysis embraces all aspects of the production dynamics and may be considered a step toward a more realistic description of production systems, especially in the context of Digital Twinning.

Hence, the main purpose of this paper is to evaluate the transient response of production systems including the homogenous and nonhomogenous reliability of machines. First, the analytical solution to the transient response of production systems will be developed using a direct solution of balance equations and the mode superposition method. Then, the same problem will be considered in the framework of the finite-state method (FSM) to address dimensionality issues of the problem including computation time and data storage demands. In all cases, the key performance indicators will be presented as a function of time elapsed from the production initialization.

The remainder of the paper is structured as follows: A brief literature review is presented in the next subsection. The analytical solution to the problem is presented in Section 2, including the approximative approach as well. Illustrative examples of homogenous and nonhomogenous serial production lines are presented in detail in Section 3, including implications for splitting and assembly systems as well. Finally, the main conclusions and prospects for further research are outlined in Section 4.

Brief Literature Review

The available literature related to the transient behavior of production systems is sparse and only a few references can be pointed out as relevant to production system engineering problems. Some of the earliest works on similar topics include the Markov drift model of a fluid movement between producers and consumers mutually coupled by a buffer [19], production losses of a single machine of exponential reliability [20], and a continuous model of two machines with finite intermediary buffer [21]. The problem of production losses due to the transient effects of serial production lines with Bernoulli machines was first considered in [16] using the aggregation procedure as an approximative approach. The problem was solved by employing the second-largest eigenvalue and pre-exponential factors, yielding the production rate and the work-in-process as functions of the time slot. In that way, production losses could be evaluated concerning the steady-state values obtained using the standard aggregation procedure [11]. The key assumption of this approach is the application of the pseudo-modal superposition method that employs only the first two eigenvalues of the transition matrix while neglecting others. However, this may be misleading as other eigenvalues may be of the same order of magnitude as the second one and may in that way impact the dynamics of the system’s transition toward its steady-state response. Based on these facts, transient analysis of production systems was pointed out as a still open area, especially regarding complex systems [22].

The same problem with the Bernoulli serial production line was addressed in [15], including transient aspects of the consumption rate as well as probabilities of blockage and starvation. All of the key performance indicators were considered in that case using the exact evaluation and the aggregation-based algorithm relying on the second-largest eigenvalue only. Besides that, the aggregation approach to the transient effects of the Bernoulli serial production lines was applied in several other cases including serial lines with two geometric machines [23], float-based systems with two Bernoulli machines [24], and dairy filling and packing systems [25]. Also, the same effects were considered in the case of serial lines with geometric machines [26], closed production lines with Bernoulli machines [27], serial lines with perishable products [28,29], and assembly systems with Bernoulli machines [30].

The governing idea behind the presented approaches to the transient analysis of production lines is the possibility of evaluating the key performance indicators more realistically by taking into account the initial state of the system as well as its movement through the system space toward steady-state behavior. However, the vast majority of considered problems are assumed to be homogenous Markov chains, while omitting nonhomogenous cases except for, according to the authors’ best knowledge, the research presented in [15]. But, in reality, the reliability of machines is a function of the time elapsed from the last scheduled maintenance. Hence, nonhomogenous aspects of the problem can contribute to the development of more realistic mathematical models reflecting both the productivity and maintainability aspects of production systems [31].

Given that, the main goal of this paper is to develop a mathematical model of Bernoulli production lines capable of evaluating transient performance by employing homogenous and nonhomogenous Markov chains. First, the analytical model will be presented using the homogenous and nonhomogenous formulation of the generalized transition matrix [12], direct solution of balance equations, and modal superposition approach. Second, the finite-state method [32,33] will be further extended to the case of transient evaluation of homogenous and nonhomogenous systems to cope with well-known CPU (Central Processing Unit) and memory storage issues of large-scale and transition-rich stochastic systems. It is expected that this new approach will enable the integration of a more realistic Markovian framework within continuously developing Digital Twinning platforms, yielding more reliable predictive analytics as well as maintenance scheduling conditioned to lesser production losses.

2. Mathematical Modeling

As outlined in the Introduction section, the mathematical modeling of Bernoulli production lines will take into account analytical and approximate approaches. The importance of the analytical solution, although computationally demanding, is reflected in the possibility of validating the approximative approaches. Also, it will enable quantification of the relationship between the number of modes taken into account within the modal superposition method and the accuracy of the results. This is of theoretical and practical interest as the second-largest eigenvalue assumption may not be sufficient to represent the dynamic behavior of a random system. Following this conclusion, the extent of the frequency band which has to be taken into account can be defined.

2.1. Analytical Solution

Consider a serial production line composed of M machines m_i, i = 1, 2, …, M where each one is of the Bernoulli reliability, and M − 1 buffers b_i, i = 1, 2, …, M − 1, each placed between two adjacent machines. Assume that all machines have identical cycle times and that the first machine can not be starved (infinite material source capacity) nor the last one can be blocked (infinite product sink capacity). Also, assume that the status of each machine is determined independently of others at the beginning of each cycle and that the system obeys the mass conservation law. This kind of production line is usually referred to as a Bernoulli production line [11] and is associated with the state space spanning over all possible combinations of buffer occupancies. The dynamics of this kind of random system can be mathematically described through a transition matrix as a Markov chain. Hence, given the initial state of the system, ${\left\{\pi \right\}}_0$, and the transition matrix $\left[P\right]$, the state of the system in the cycle n, ${\left\{\pi \right\}}_n$, can be determined as [34,35]

$${\left\{\pi \right\}}_n={\left[P\right]}^n{\left\{\pi \right\}}_0,$$

(1)

or alternatively as

$${\left\{\pi \right\}}_n=\left[P\right]{\left\{\pi \right\}}_{n-1}, n\ge 1.$$

(2)

Given the transition matrix, [P], Equations (1) and (2) can be solved directly for each cycle n. Indeed, the available formulations of transition matrices in the case of serial, splitting, and assembly Bernoulli production lines [12,13,15] stand as the key enablers of a direct approach to the problem. However, for complexity reasons, the direct solution of Equation (1) or (2) quickly becomes rather cumbersome, especially in the case of large-scale, dense, and transition-rich systems. Also, we are often interested in the relationship between the system properties reflected in the transition matrix and the key performance indicators. Hence, a different path such as the eigendecomposition and the modal superposition method may be of interest.

The eigendecomposition of a transition matrix, [P], yields

$$\left[P\right]=\left[U\right]\left[\Lambda \right]{\left[U\right]}^{-1},$$

(3)

where [U] is a matrix whose columns are the right eigenvectors of [P], and [Λ] is the diagonal matrix with corresponding eigenvalues λ₁, λ₂, …, λ_d, $d={\displaystyle \prod _{e=1}^{M-1}\left(N_e+1\right)}$. If the right eigenvectors ${\left\{u\right\}}_i$ are normalized using the left eigenvectors ${〈v〉}_i$ in a way that a relationship ${〈u〉}_i{\left\{v\right\}}_i=1$, i = 1, 2, …, d, holds, Equation (3) can be simplified into

$$\left[P\right]=\left[U\right]\left[\Lambda \right]{\left[V\right]}^{\mathrm{T}},$$

(4)

using the identity ${\left[U\right]}^{-1}={\left[V\right]}^{\mathrm{T}}$, where [V] is a matrix whose rows are the left eigenvectors of [P]. Consequently, introducing Equation (4) into Equation (2), the following relationship holds

$${\left\{\pi \right\}}_n=\left[U\right]\left[\Lambda \right]{\left[V\right]}^{\mathrm{T}}{\left\{\pi \right\}}_{n-1}, n\ge 1,$$

(5)

which can be presented more conveniently in the form of a sum as

$${\left\{\pi \right\}}_n={\displaystyle \sum _{i=1}^d{\lambda }_i}{\left\{u\right\}}_i{〈v〉}_i{\left\{\pi \right\}}_{n-1},$$

(6)

also known as the modal superposition approach. Each summand in Equation (6) represents a probability distribution vector associated with the ith eigenvalue. Consequently, each key performance indicator associated with the nth cycle may be presented in the form of a sum, where each summand will correspond to a particular eigenvalue.

A very important concept justifying the application of Markov chains in the case of production system engineering is the ergodic theorem ensuring that the distribution ${\left\{\pi \right\}}_n$ is unique for each n and limiting when $n\to \infty$. In the homogenous case, a Markov chain is ergodic if it is irreducible (all system states communicate with each other) and aperiodic (at least one diagonal element of a transition matrix is larger than 0). Provided that it remains irreducible and aperiodic, the nonhomogenous Markov chain must meet additional criteria on weak and strong ergodicity to comply with standard concepts of production system engineering. The weak ergodicity criterion is based on Dobrushin’s ergodic coefficient, δ,

$$\delta =\frac{1}{2}\underset{i,j}{\sup }{\displaystyle \sum _k\left|p_{ik}\left(n\right)-p_{jk}\left(n\right)\right|},$$

(7)

representing the supremum of the difference between matrix rows. The weak condition is met if $\underset{n\to \infty }{\lim }\delta =0$, that is if rows of the transition matrix become more and more similar with each cycle. Also, the condition

$${\displaystyle \sum _{s=0}^{\infty }\left(1-\delta \left[P\left(n_s,n_{s+1}\right)\right]\right)}=\infty ,$$

(8)

must hold for the family of transition matrices $\left[P\left(n_s,n_{s+1}\right)\right]$ where n_s is a strictly increasing sequence of integers [34]. Each weakly ergodic Markov chain is also strongly ergodic if a unique steady-state distribution exists in each cycle n, which is by definition valid for all aperiodic and irreducible random systems. As demonstrated earlier in [12,13,14,15], the transition matrices associated with serial, splitting, and assembly Bernoulli production lines stand for irreducible and aperiodic random processes and have unique and limiting distributions referred to as the steady-state probability distributions. Since in that way the conditions in Equations (7) and (8) are met, Equation (6) can be extended to the nonhomogenous case as

$${\left\{\pi \right\}}_n={\displaystyle \sum _{i=1}^d{\lambda }_i\left(m\right)}{\left\{u\left(m\right)\right\}}_i{〈v\left(m\right)〉}_i{\left\{\pi \right\}}_{n-1}, n,m\ge 1,$$

(9)

where the eigenvalues ${\lambda }_i\left(m\right)$ and right and left eigenvectors, ${\left\{u\left(m\right)\right\}}_i$ and ${〈v\left(m\right)〉}_i$, are functions of the second-order cycle, m, representing a nonhomogenous transition matrix [P(m)] with entries p_ij(m). Given the standard assumptions of the Bernoulli production lines, a ratio m/n may take two distinct values: m/n > 1 and m/n = 1, which can be referred to as the sub-resonant and resonant nonhomogenous Bernoulli production lines. The over-resonant case, m/n < 1, implies variability of machine reliability during cycle time and for that reason does not meet the set of standard assumptions. Hence, it can be addressed only as a reduction to the resonant case if the cycle time is condensed to some subset of production activities.

Finally, by observing Equation (9), the problem of the transient response reduces to the eigenvalue problem employed to calculate eigenvalues and the associated right and left eigenvectors of a transition matrix for each cycle n, respectively, second-order cycle m. Both eigenvalues and eigenvectors can be real numbers and may also appear as complex conjugate pairs [34,36,37], while the first (and the largest) eigenvalue always takes unit value. Thus, given the transition matrix, they may be calculated using some of the well-known procedures like the Power or QR method [38]. However, for dimensionality issues, it may easily become inconvenient especially when realistic production systems are considered. For these reasons, an approximative procedure to evaluate the transient response of nonhomogenous Bernoulli production lines will be presented and validated against corresponding analytical solutions. To do that, the transient response of a nonhomogeneous Bernoulli production line composed of two machines and one buffer will be considered and the results will be used as a basis for an approximative solution using the finite-state method.

2.2. Two-Machine One-Buffer Bernoulli Production Line

Consider a serial production line composed of two machines of Bernoulli reliability and one buffer of capacity N. Then, the associated transition matrix, [P(m)], takes a well-known form [11]

$$\left[\begin{array}{cccccccccc}1-p_1& \left(1-p_1\right)p_2& 0& 0& 0& \dots & & \dots & 0& 0\\ 0& D& \left(1-p_1\right)p_2& 0& \ddots & & & & ⋮& ⋮\\ ⋮& p_1\left(1-p_2\right)& D& \left(1-p_1\right)p_2& \ddots & & & & & \\ & 0& p_1\left(1-p_2\right)& D& \ddots & & & & & \\ & ⋮& 0& \left(1-p_1\right)p_2& \ddots & & & & & \\ 0& & ⋮& 0& \ddots & & & & & \\ p_1& & & ⋮& \ddots & & & \ddots & 0& \\ 0& & & & & & & \ddots & \left(1-p_1\right)p_2& 0\\ ⋮& & & & & & & \ddots & D& p_1\left(1-p_2\right)\\ 0& 0& 0& 0& 0& \dots & & \dots & p_1\left(1-p_2\right)& p_1{p}_2+\left(1-p_2\right)\end{array}\right]$$

(10)

where $D=\left(1-p_1\right)\left(1-p_2\right)+p_1{p}_2$. Note that the reliabilities p₁ and p₂ are functions of the second-order cycle, m; however, for clarity reasons, this notation was omitted from Equation (10). The right eigenvector ${\left\{u\left(m\right)\right\}}_i$ corresponding to eigenvalue λ_i and second-order cycle m equals

$${{\left\{u\left(m\right)\right\}}_i}^{\mathrm{T}}=〈\begin{array}{ccccc}{P}_{i,0}& P_{i,1}& P_{i,2}& \cdots & P_{i,N}\end{array}〉,$$

(11)

where

$$\begin{array}{c}{P}_{i,1}=\frac{-\left(1-p_1-{\lambda }_i\right)P_{i,0}}{\left(1-p_1\right)p_2},\\ P_{i,2}=\frac{-p_1{P}_{i,0}-\left(D-{\lambda }_i\right)P_{i,1}}{\left(1-p_1\right)p_2},\\ P_{i,j}=\frac{-p_1\left(1-p_2\right)P_{i,j-2}-\left(D-{\lambda }_i\right)P_{i,j-1}}{\left(1-p_1\right)p_2}, j=3,4,\dots ,N.\end{array}$$

(12)

Similarly, the left eigenvector ${〈v\left(m\right)〉}_i$ takes the form

$${〈v\left(m\right)〉}_i=〈\begin{array}{ccccc}{Q}_{i,0}& Q_{i,1}& Q_{i,2}& \cdots & Q_{i,N}\end{array}〉,$$

(13)

where

$$\begin{array}{c}{Q}_{i,0}=1,\\ Q_{i,1}=\frac{-\left(1-p_1-{\lambda }_i\right)Q_{i,0}}{p_1},\\ Q_{i,j}=\frac{-\left(1-p_1\right)p_2{Q}_{i,j-2}-\left(D-{\lambda }_i\right)Q_{i,j-1}}{p_1\left(1-p_2\right)}, j=2,3,\dots ,N.\end{array}$$

(14)

The remaining unknowns ${\lambda }_i$ and $P_{i,0}$ can be determined using additional conditions. In the first case, the last row of the transition matrix, Equation (10), may be employed since the system of equations is overdetermined. Hence, it follows that the unknown eigenvalue can be calculated iteratively using a relationship

$$P_{i,N}=\frac{-p_1\left(1-p_2\right)}{p_1{p}_2+\left(1-p_2\right)-{\lambda }_i}{P}_{i,N-1}.$$

(15)

In the second case, the unknown probability $P_{i,0}$ equals

$$P_{i,0}=\frac{1}{{〈u〉}_i{\left\{v\right\}}_i},$$

(16)

since the eigendecomposition of a transition matrix is based on normalization through the left eigenvectors. It may be easily verified that in the case when i = 1 and ${\lambda }_1=1$, Equations (11)–(16) result in a left eigenvector taking a unit value and the right one corresponding to the steady-state probability distribution of a system at hand, as presented in [11]. Also, note that according to Equation (15), ${\lambda }_i\in ℝ$. This relationship is valid only for a line composed of two machines and one buffer. However, this is not necessarily the case for nonsymmetric matrices with real entries such as transition matrices associated with lines composed of more than two machines. In these more complex cases, eigenvalues may take real or complex values, with the latter appearing in complex conjugate pairs.

2.3. The Finite-State Method

The finite-state method is based on the discretization of the system’s state space using finite-state elements and the associated proportionality property of element-level probability distributions, resulting in a possibility to assume the independence of events at each buffer and to reconstruct the probability distribution at the system level, ${\left\{\pi \right\}}_n$. This approach has so far been successfully applied and validated against rigorous analytical results in cases of the steady-state response of serial [14], splitting [11], and assembly [15] Bernoulli production lines. Here, as a further generalization, it will be extended to the transient behavior of Bernoulli production lines assuming both homogenous and nonhomogenous definitions of machine reliabilities, i.e., of the transition matrix.

The finite-state method employs a set of finite-state elements, each composed of only two machines and one intermediate buffer. Depending on the particular arrangement of the production line, each finite-state element may be characterized as upstream or downstream conditioned to the position of the least reliable machine, m_m. In the first case, each upstream element is composed of machine m_e, buffer b_e, and machine m_m where e < m, while in the second case, each downstream element consists of machine m_m, buffer b_e, and machine m_e+1, where e ≥ m. Consequently, in the upstream case reliabilities p₁ and p₂, Equations (11)–(16) are identical to p_e, e < m, respectively, p_m, while in the downstream case, they are identical to p_m, respectively, p_e+1, e ≥ m.

Hence, the unknown eigenvalues, ${\lambda }_i^e$, and eigenvectors, ${\left({\left\{u\left(m\right)\right\}}_i^e\right)}^{\mathrm{T}}={〈\begin{array}{ccccc}{P}_{i,0}& P_{i,1}& P_{i,2}& \cdots & P_{i,N}\end{array}〉}^e$ and ${〈v\left(m\right)〉}_i^e={〈\begin{array}{ccccc}{Q}_{i,0}& Q_{i,1}& Q_{i,2}& \cdots & Q_{i,N}\end{array}〉}^e$, at the level of finite-state element e, e = 1, 2, …, M − 1 can be determined using Equations (11)–(16). Given the eigenvalues ${\lambda }_i^e$ and the independence of events assumption, the system-level eigenvalues, ${\lambda }_i$, equal

$${\lambda }_i={\lambda }_j^{e=1}{\lambda }_k^{e=2}\dots {\lambda }_l^{e=M-1},$$

(17)

where $i=0,1,\dots ,d$, $d={\displaystyle \prod _{e=1}^{M-1}\left(N_e+1\right)}$, and $j,k,\dots ,l=0,1,\dots ,N_e$. Similarly, the element f of the system-level eigenvectors, $u{\left(m\right)}_{i,f}$ and $v{\left(m\right)}_{i,f}$, associated with eigenvalue ${\lambda }_i$, Equation (17), equals

$$\begin{array}{l}{u}_{i,f}\left(m\right)=P_{i,j}^{e=1}{P}_{i,k}^{e=2}\dots P_{i,l}^{e=M-1},\\ v_{i,f}\left(m\right)=Q_{i,j}^{e=1}{Q}_{i,k}^{e=2}\dots Q_{i,l}^{e=M-1},\end{array}$$

(18)

where $i,f=0,1,\dots ,d$, $d={\displaystyle \prod _{e=1}^{M-1}\left(N_e+1\right)}$, and $j,k,\dots ,l=0,1,\dots ,N_e$. Note that each element f of the right and left eigenvector corresponds to a unique combination of indices j, k, …, l reflecting the occupancy level of each buffer. All combinations may thus be represented as a Cartesian product

$$\begin{gathered} \left\{\left(j,k,\dots ,l\right)=(h_1,h_2,\dots ,h_{M-1})\right\}=\left\{j=h_1,h_1\in 0,1,\dots ,N_1\right\}\times \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\left\{k=h_2,h_2\in 0,1,\dots ,N_2\right\}\times \dots \times \left\{l=h_{M-1},h_{M-1}\in 0,1,\dots ,N_{M-1}\right\}, \end{gathered}$$

(19)

where $h_e, e=0,1,\dots ,M-1$ is the occupancy level of the eth buffer.

2.4. Key Performance Indicators

The key performance indicators may be determined once the system-level eigenvalues and corresponding eigenvectors as well as the probability distribution vector are known. Although many different indicators can be employed in realistic cases, here we will focus only on the production rate of the sth machine, PR_s, the work-in-process at the eth buffer, WIP_e, and the probabilities of blockage and starvation of the sth machine, BL_s and ST_s, as generic key performance indicators. Other indicators, such as the expected lead time, expected energy consumption, expected production costs, expected profitability, or efficiency, are considered case-specific. Hence, they will be omitted from this paper. However, once the probability distribution vector is defined, their evaluation is quite straightforward.

The production rate, PR_s, stands for the expected number of products delivered by the machine s. Usually, the production rates of the last machine in the serial or assembly line or of the last machine in each branch of the splitting system are of interest. Let us, for simplicity, first consider the simplest case of the serial production line. In that case, the production rate of the machine M equals

$$P{R}_M\left(n,m\right)=P\left(\left\{b_{M-1} \mathrm{is} \mathrm{not} \mathrm{empty}\right\}\cap \left\{m_M \mathrm{is} \mathrm{not} \mathrm{down}\right\}\right).$$

(20)

By employing Equations (9) and (18), Equation (20) yields

$$P{R}_M\left(n,m\right)=p_M\left(1-{\displaystyle \sum _{h_1=0}^{N_1}{\displaystyle \sum _{h_2=0}^{N_2}\dots {\displaystyle \sum _{h_{M-2}=0}^{N_{M-2}}{\displaystyle \sum _{i=1}^d{\lambda }_i{P}_{i,h_1}^{e=1}{P}_{i,h_2}^{e=2}\dots P_{i,h_{M-2}}^{e=M-2}{P}_{i,0}^{e=M-1}{〈v〉}_i{\left\{\pi \right\}}_{n-1}}}}}\right),$$

(21)

that for proportionality reasons [32] can be simplified into

$$P{R}_M\left(n,m\right)=p_M\left(1-{\displaystyle \sum _{i=1}^d{\lambda }_i{P}_{i,0}^{e=M-1}{〈v〉}_i{\left\{\pi \right\}}_{n-1}}\right).$$

(22)

By the same logic, the production rate associated with machine s equals

$$P{R}_s\left(n,m\right)=p_s\left(1-{\displaystyle \sum _{i=1}^d{\lambda }_i{P}_{i,0}^{e=s-1}{〈v〉}_i{\left\{\pi \right\}}_{n-1}}\right), s=2,3,\dots ,M.$$

(23)

The work-in-process, WIP_e, states the average number of parts contained at the buffer b_e and can be expressed for each finite-state element, e, as

$$WI{P}_e\left(n,m\right)={\displaystyle \sum _{h_1=0}^{N_1}{\displaystyle \sum _{h_2=0}^{N_2}\dots {\displaystyle \sum _{h_{M-1}=0}^{N_{M-1}}{h}_e{P}_{h_1{h}_2\dots h_{M-1}}}}}.$$

(24)

Hence, by employing Equations (9) and (18), Equation (24) can be presented as

$$WI{P}_e\left(n,m\right)={\displaystyle \sum _{h_1=0}^{N_1}{\displaystyle \sum _{h_2=0}^{N_2}\dots {\displaystyle \sum _{h_{M-1}=0}^{N_{M-1}}{\displaystyle \sum _{i=1}^d{\lambda }_i{h}_e{P}_{i,h_1}^{e=1}{P}_{i,h_2}^{e=2}\dots P_{i,h_{M-1}}^{e=M-1}{〈v〉}_i{\left\{\pi \right\}}_{n-1}}}}},$$

(25)

which for proportionality reasons simplifies into

$$WI{P}_e\left(n,m\right)={\displaystyle \sum _{h_e=0}^{N_e}{\displaystyle \sum _{i=1}^d{\lambda }_i{h}_e{P}_{i,h_e}^e{〈v〉}_i{\left\{\pi \right\}}_{n-1}}}, e=1,2,\dots M-1.$$

(26)

Further, the probabilities of blockage and starvation, BL_s and ST_s, stand for probabilities of events

$$\begin{array}{c}B{L}_s=P\left(\left\{m_s \mathrm{is} \mathrm{up}\right\}\cap \left\{b_s \mathrm{is} \mathrm{full}\right\}\cap \left\{m_{s+1} \mathrm{is} \mathrm{down}\right\}\right.\\ \cup \left.\left\{m_s \mathrm{is} \mathrm{up}\right\}\cap \left\{b_s \mathrm{is} \mathrm{full}\right\}\cap \left\{m_{s+1} \mathrm{is} \mathrm{blocked}\right\}\right), s=1,2,\dots ,M-1,\\ S{T}_s=P\left(\left\{m_s \mathrm{is} \mathrm{up}\right\}\cap \left\{b_{s-1} \mathrm{is} \mathrm{empty}\right\}\right), s=2,3,\dots ,M.\end{array}$$

(27)

Again, by employing Equations (9) and (18), Equation (27) after simplification yields

$$\begin{array}{c}B{L}_s\left(n,m\right)=p_s\left(1-p_{s+1}+B{L}_{s+1}\right){\displaystyle \sum _{i=1}^d{\lambda }_i{P}_{i,N_e}^{e=s}{〈v〉}_i{\left\{\pi \right\}}_{n-1}}, s=1,2,\dots ,M-1; B{L}_M=0,\\ S{T}_s\left(n,m\right)=p_s{\displaystyle \sum _{i=1}^d{\lambda }_i{P}_{i,0}^{e=s-1}{〈v〉}_i{\left\{\pi \right\}}_{n-1}}, s=2,3,\dots ,M; S{T}_1=0.\end{array}$$

(28)

It is also of interest to consider the key performance indicators in a special case when the transition matrix is homogenous and when the initial distribution ${\left\{\pi \right\}}_0$ states that all buffers are empty. In that case, since $Q_{i,1}=1$, ${\left({\left\{\pi \right\}}_0\right)}^{\mathrm{T}}=〈\begin{array}{cccc}1& 0& \dots & 0\end{array}〉$, and ${〈v〉}_i{\left\{\pi \right\}}_0=1$, Equations (23), (26), and (28) can be simplified into

$$\begin{array}{c}P{R}_s\left(n\right)=p_s\left(1-{\displaystyle \sum _{i=1}^d{\lambda }_i^n{P}_{i,0}^{e=s-1}}\right), s=2,3,\dots ,M,\\ WI{P}_e\left(n\right)={\displaystyle \sum _{h_e=0}^{N_e}{\displaystyle \sum _{i=1}^d{\lambda }_i^n{h}_e{P}_{i,h_e}^e}}, e=1,2,\dots ,M-1,\\ B{L}_s\left(n\right)=p_s\left(1-p_{s+1}+B{L}_{s+1}\right){\displaystyle \sum _{i=1}^d{\lambda }_i^n{P}_{i,N_e}^{e=s}}, s=1,2,\dots ,M-1; B{L}_M=0,\\ S{T}_s\left(n\right)=p_s{\displaystyle \sum _{i=1}^d{\lambda }_i^n{P}_{i,0}^{e=s-1}}, s=2,3,\dots ,M; S{T}_1=0.\end{array}$$

(29)

In addition, note that when i = 1, λ_i takes a unit value and Equation (28) becomes identical to the key performance indicators valid in the steady-state case [32]. Finally, it may seem that the key performance indicators presented in Equation (28) are not dependent on the initial distribution ${\left\{\pi \right\}}_0$. However, such a conclusion may be misleading. Although it does not appear directly in Equation (28), it will indirectly influence the key performance indicators through the gradual buildup of available system states. In other words, as the cycles pass one after another more and more system states will be included in the state space until its full dimension d is achieved.

2.5. KPI Evaluation Algorithm

As already pointed out, the main goal of the presented analytical and approximative procedures is to evaluate KPIs associated with different production systems. Hence, a KPI evaluation algorithm was developed in a general-purpose and imperative programming language using Equations (20)–(28) with the probability distribution vector, $\left\{\pi \right\}$, determined according to Equation (9) in the case of an analytical approach and Equation (18) in the case of an approximative approach. In both cases, the eigendecomposition is employed to cope with dimensionality issues of the eigenvalue problem. Additionally, the finite-state method is employed to bypass the issue of the transition matrix formulation using a set of two-machine one-buffer elements and the independence of events assumption. The developed algorithm and model assumptions are summarized in Figure 1.

3. Application of the Developed Theory

The introduced theory has been applied in benchmark cases of Bernoulli serial lines composed of five machines and four buffers, each with a capacity of 5, to validate the approximative approach against the analytical solution to the problem. Three different cases were considered, namely the transient response of a homogenous serial Bernoulli production line, as well as the transient response of the sub-resonant and resonant nonhomogenous Bernoulli production lines. The considered input data are outlined in

Table 1. The input data for benchmark cases, N₁ = N₂ = N₃ = N₄ = 5.

Reliability of Machine mi		Homogenous		Nonhomogenous
		/		m/n > 1		m/n = 1
		m	Δpi	m	Δpi	m	Δpi
p1	0.4	/	/	10n	0.01	n	0.001
p2	0.5	/	/	10n	0.01	n	0.001
p3	0.6	/	/	10n	0.01	n	0.001
p4	0.7	/	/	10n	0.01	n	0.001
p5	0.8	/	/	10n	0.01	n	0.001

including reliabilities of machines, p_i, second-order cycles, m, and machine deterioration rates, Δp_i. In all cases, the initial state of the system was set to 0 at each buffer, while planned maintenance of all machines was assumed between cycles 100 and 101 for nonhomogenous lines. Each case was first evaluated using the analytical approach (AN), Equation (2), to obtain the probability distribution vector for each cycle, n. In that way, the key performance indicators were determined using Equations (20), (24), and (27). Similarly, the probability distribution vectors and the associated key performance indicators were calculated using the finite-state method (FSM), Equations (18), (22), (26), and (28). A comparison between the obtained results is presented in Figure 2, Figure 3, Figure 4 and Figure 5 to justify the application of the finite-state method in more complex cases involving significantly larger system states.

A comparison between Re(λ) calculated using the AN and FSM approach in the case of a homogenous serial Bernoulli line is presented in Figure 2. Quite good agreement between both methods may be noted, especially in the range of the most influential system states for d < 500, while some minor discrepancies are present for d > 1100 since, by definition, Equation (17), the FSM yields only positive and real-valued eigenvalues based on eigenvalues of each finite-state element. This fact, as well as the absence of the imaginary part of the eigenvalue, does not impact the accuracy of the FSM since Im(λ) << Re(λ) and Im({p_n}) << Re({p_n}).

A detailed comparison between KPIs (the production rate, PR, the work-in-process, WIP_i, the probability of blockage, BL_i, and the probability of starvation, ST_i) obtained using the analytical (AN) and the approximative approach (FSM) is presented in Figure 3, Figure 4 and Figure 5 in cases of homogenous, sub-resonant nonhomogenous, and resonant nonhomogenous serial Bernoulli production lines, Table 1. In all cases, quite good agreement between the obtained results may be noted except in the case of the earliest phase of the system operation where FSM demonstrates a more pronounced gradient of KPIs concerning the cycle number, n. This effect is present since a gradual increase in the available system states results in a rather coarse discretization of the state space. It, however, diminishes once all available states become reachable.

In the case of a homogenous Bernoulli production line, all KPIs asymptotically approach the steady-state values which, in cases of the production rate and the probability of starvation, equal $P{R}_{\mathrm{steady}\text{-}\mathrm{state}}=\min \left\{p_i\right\}$, i = 1, 2, …, M, respectively, $S{T}_{i, \mathrm{steady}\text{-}\mathrm{state}}=p_i-P{R}_{\mathrm{steady}\text{-}\mathrm{state}}$, i = 2, 3, …, M, while the work-in-process and the probability of blockage converge to a constant value given the ascending arrangement of the machine reliabilities. Also, the transport delay related to the time required for input to reach output may be noted in the case of all KPIs, especially in the case of the production rate, Figure 3a, taking zero value for n < 4, or more generally for n < M − 1 [11]. On the other hand, the transient response of both sub-resonant and resonant nonhomogenous Bernoulli production lines, Figure 4 and Figure 5, closely resembles the dynamics of the machine deterioration rates. Thus, in the first case, all KPIs result in stepwise functions, Figure 4, while in the second case, they become continuous except for the stepwise increase due to scheduled maintenance, Figure 5. In addition, the transportation delay is present in both cases for all KPIs, similar to the homogenous case.

As already pointed out in the Introduction section, an important issue related to the transient response evaluation of production systems is the number of eigenvalues taken into account within the modal superposition method. Its reduction may decrease the CPU requirements; however, the associated impact on the accuracy of the results is not known. Presently, the only known reference to this issue may be found in [11], where the second largest eigenvalue is associated with the duration of transient behavior as a function of the system properties. Therefore, a sensitivity analysis concerning the number of eigenvalues and the accuracy of the results was performed for the resonant nonhomogenous Bernoulli serial line as the most representative case. First, the number of eigenvalues was decreased to five by excluding the smallest eigenvalue at the level of finite-state elements. In that way, the scale of the problem was reduced to d = 625 (recall that d stands for the total number of eigenvalues at the level of the production system). Similarly, the sensitivity analysis was performed by taking into account four, three, two, and only the largest eigenvalue for each finite-state element leading to the problem scales d = 256, d = 81, d = 16, and d = 1. In all cases, a comparison against the full-scale problem including six eigenvalues per element and d = 1296 was enabled. Surprisingly, the reduction in the problem scale did not impact the final results at all, except for the case d = 1 where negligible disagreement between the results was detected during the first 10 to 50 cycles, depending on the particular KPI, Figure 6, Figure 7, Figure 8 and Figure 9. For clarity reasons, Figure 6, Figure 7, Figure 8 and Figure 9 bring cases d = 1296, d = 16, and d = 1, while all other cases remain identical to the referent one. Also, Figure 6, Figure 7, Figure 8 and Figure 9 present only cycles associated with some influence of the problem scale reduction, while in all other cycles each KPI remains the same as the referent one, i.e., it is not affected by the problem scale reduction. Consequently, a significant reduction in the problem scale may be introduced for practical problems without impacting the accuracy of the KPI evaluations. In the most extreme case, the problem of the transient response may be even reduced to d = 1 by taking into account only steady-state probability distributions, $\left\{u_1^e\right\}$ and $\left\{v_1^e\right\}$, and the associated eigenvalues, ${\lambda }_1^e=1$, for each finite-state element.

3.1. Implications for Splitting and Assembly Systems

The developed theory has so far been presented in the case of serial production lines. However, practical engineering problems often involve production systems composed of splitting and assembly lines as well. Consequently, the transient response of such systems is of great interest and it may be addressed using the finite-state method and the outlined theory. This conclusion may be justified by the fact that the application of the finite-state framework in the case of splitting [13] and assembly [15] lines boils down to a set of mutually communicating serial lines representing primary and secondary material flows. In the first case, this communication is realized through splitting factors, while in the second case it is governed by states of buffers placed immediately before machines involved in merging operations. Hence, once formulated, the finite-state elements of the splitting or assembly production line are combined according to the proportionality property into a system-leveled eigenvector and employed further on in Equation (9) to determine the associated probability distribution vector, {p}_n, in cycle n. Finally, all conclusions related to the transient response of serial Bernoulli production lines are also valid in the splitting and assembly case as well.

3.2. Application Case

The introduced transient response evaluation of the Bernoulli production system was applied in the case of a wood processing facility producing flooring elements, Figure 10. The steady-state response of the same system was already addressed in [13], including factory-floor acquisition of machine reliabilities, storage capacities, splitting rates, and reliability of the acquired data. Hence, only its brief description will be provided here, while focusing on the associated transient behavior and its comparison to the steady-state values. Additionally, transient properties of production systems may heavily impact (especially in terms of financial losses) facilities with large series of products. It is, thus, of interest to demonstrate the application of the derived theory in the case of a wood processing facility as one representative of high-volume production.

The production facility consists of seven machines where operations of log debarking, band sawing, circular sawing, and multiple rip sawing are performed. The reliability of each machine, the capacity of buffers, and splitting rates were acquired through the factory-floor measurement,

Table 2. The reliability of machines, the capacity of buffers, and splitting rates [13].

	Wood Processing Facility
Reliability of machines	p1	p2	p3	p4	p5	p6	p7
	1	0.956	0.985	0.992	0.962	0.979	0.985
Capacity of buffers	N1	N2	N3	N4	N5	N6
	13	6	4	6	11	6
Rates	r1	r2	r3	r4	r5	r6
	0.4	0.6	0.8	0.2	0.7	0.3

. The final product of the facility is 2-m-long and 50-mm-thick wooden beams produced at the machine m₇, while all other products including sawn boards of thickness less than 50 mm, beams of thickness less than 50 mm, and beams of length less than 2 m are stored immediately after operations m₃, m₆, and m₇ for further processing. The transient response of the system at hand was evaluated for a period of 2200 cycles starting with empty buffers. In that way, a typical working shift of the facility was simulated using a nonhomogenous Bernoulli splitting production line. The finite-state method was employed and the deterioration rate, Δp_i, of 0.01 per 100 cycles was assumed for each machine. Only the largest eigenvalue and the associated eigenvector were employed for each finite-state element. The obtained KPIs are presented In Figure 11, Figure 12 and Figure 13 as functions of the cycle number, n, including production rates of machines m₃, m₆, and m₇, work-in-process contained at buffers b₂, b₃, b₄, b₅, and b₆, and probabilities of starvation of machines m₄, m₅, m₆, and m₇. As probabilities of blockage take values well below 10⁻⁴, we are omitting them from the presentation of the results.

The obtained results, Figure 11, Figure 12 and Figure 13, resemble the behavior of the wood processing facility during one working shift. The stepwise nature of all KPIs may be noted as a result of the nonhomogenous effects and machine deterioration rates. Also, all KPIs associated with the 100th cycle remain equal to KPIs obtained for the steady-state and homogenous response of the same system,

Table 3. Comparison between KPIs: steady-state response of the homogenous splitting line, transient response of the nonhomogenous splitting line, n = 100, and the factory floor data.

KPIs	Nonhomogenous Case, n = 100	Steady-State, Homogenous Case [13]	Factory-Floor Data [13]
PR3	0.382	0.382	n/a
PR6	0.115	0.114	n/a
PR7A	0.321	0.321	0.331
PR7B	0.138	0.138	0.176
WIP2	0.391	0.392	1.093
WIP3	0.585	0.585	0.506
WIP4	0.630	0.630	0.909
WIP5	0.118	0.117	0.695
WIP6	0.472	0.472	0.283
ST3	0.603	0.603	n/a
ST4	0.418	0.418	n/a
ST5	0.388	0.388	n/a
ST6	0.864	0.865	n/a
ST7	0.526	0.526	n/a

. The reliability of the approach may also be noted if its results are compared with the data acquired at the factory floor [13]. This is especially valid in the case of PR_7A and PR_7B, with some discrepancies related to the work-in-process. These discrepancies are dominantly caused by the measurement uncertainty conditioned to the factory floor environment. Given the properties of the system, production rates and probabilities of starvation decrease with the increase in the cycle number, Figure 11 and Figure 13, while work-in-process increases simultaneously, Figure 12. The decrease and increase in KPIs closely obey the linear relationship between machine deterioration rates and the number of cycles with gradients dominated by splitting rates. The most interesting KPI of the system at hand is PR_7A, which may be pointed out as quite sensitive to the nonhomogenous effects of the production system leading to an expected production loss of about 73 pieces per working shift. By considering its long-term effects, this loss could be decreased if preventive maintenance is introduced to the system at the point when the accumulated financial loss meets the maintenance cost. Otherwise, it will be decreased only partly.

4. Conclusions

Mathematical modeling has an important role in the development of production system engineering, especially when Digital Twinning, the cost-effectiveness, and energy efficiency of manufacturing processes are taken into account. These can be considered within the Markovian framework using analytical and approximative approaches. Both of them were addressed in this paper. First, the analytical solution to the problem was developed as a direct solution to the balance equations by employing the eigendecomposition of homogenous and nonhomogenous transition matrices. Second, the same problem was solved using the finite-state method to tackle issues such as CPU and memory storage requirements associated with large-scale and transition-rich stochastic systems. A good agreement between the obtained results was proven in both cases.

In addition to that, a sensitivity analysis regarding the relationship between the number of eigenvalues and the accuracy of the results was performed pointing to the possibility of additional CPU and memory storage reduction by employing only the first two eigenvalues per finite-state element, or even by employing the steady state solution only. Finally, the developed theory was applied in the case of the transient response of a wood processing facility producing flooring elements. The mathematical model of the facility was developed as a sub-resonant nonhomogenous Bernoulli splitting production line, and all of the relevant KPIs were evaluated for each of the 2200 cycles representing one working shift. The assumed machine deterioration rate resulted in an expected loss of about 73 pieces of final product per one working shift. Hence, it was suggested to develop a maintenance strategy that would minimize these losses concerning maintenance costs.

The transient response of production systems is a challenging topic that should be further extended to cases including machine reliabilities defined using geometric, exponential, or non-exponential probability distributions as well as by assuming asynchronous relationships between machine cycles. This implies a broader incorporation of nonhomogenous aspects into the modeling and analysis of production systems. In addition, nonhomogenous properties of the system could be extended to cases including limited material sources or sinks to simulate the behavior of systems more realistically, including possible failures in the raw material supply or final product delivery. The outlined theory could also be applied effectively in the development of different preventive maintenance strategies by considering different aspects of cost–benefit analysis, especially in the context of Digital Twinning. Finally, it is well-known that stochastic modeling is relevant to a variety of scientific areas including not only mass, job-shop, or project-based production but also finance, physics, communication systems, biology, ecosystem, and others. Hence, the presented research may be also applied beyond production system engineering including the diversity of dense, transition-rich, homogenous and nonhomogenous systems.