On Reliability of a Double Redundant Renewable System with a Generally Distributed Life and Repair Times

Abstract

The paper provides reliability analysis of a cold double redundant renewable system assuming that both life-time and repair time distributions are arbitrary. The proposed approach is based on the theory of decomposable semi-regenerative processes. We derive the Laplace–Stieltjes transform of two main reliability measures like the distribution of the time between failures and the time to the first failure. The transforms are used to calculate corresponding mean times. It is further derived in closed form the time-dependent and time stationary state probabilities in terms of the Laplace transforms. Numerical results illustrate the effect of the type of distributions as well as their parameters on the derived reliability and probabilistic measures.

1. Introduction

The use of parallel or redundant units is an often recommended way of increasing the system reliability. In the last few decades, many authors have studied the redundant systems under different sets of assumptions about the system’s attributes and level of generality, the distribution of life, and repair times when one or more units are standby. In this paper, we formulate the model where both of distributions are assumed to be arbitrary. Such models are interesting both from theoretical and practical points of view. Theoretically, these studies are related to the creation and development of new mathematical methods such as regenerative methods, especially methods of the Decomposable Semi-Regenerative Processes (DSRP), and Markovization methods via introduction of the supplementary variables. In practical terms, the tasks under consideration form a basis for analyzing the reliability of different complex systems with generally distributed times.

Smith’s idea about regenerative stochastic processes [1] found many generalizations and applications. This approach still helps to solve many applied problems since it reduces the complexity of the task. Instead of a whole time period, the process can be analyzed within a separate regenerative period of time. However, if the process describing the system behavior is complicated enough in order to be exhaustively analyzed on the regenerative period, even in this case, different generalizations of the regenerative idea are possible. One possible generalization is to combine Smith’s idea with a semi-Markov approach [2], which seems quite appropriate for real models. The result of such combination is a semi-regenerative process. In Rykov and Yastrbrnrtsky [3], it appears under the name “Regenerative processes with several types of regeneration points”. Klimov [4] considered it as a “semi-Markov process with additional trajectories”. In accordance with Jakod [5] and Nummeline [6], in this paper, we will call such process a semi-regenerative process (SRP).

The next possible generalization is to construct within the basic regeneration cycle other regeneration points called embedded regeneration points. Sometimes, in a particular model, one can find some embedded regeneration points and to represent the process distribution in a separate regeneration cycle in terms of its distribution in embedded regeneration cycles of the second level. Initially, such processes appeared in [7,8] under the name “Regenerative processes with embedded regeneration cycles” or “Decomposable Semi-Regenerative Processes” (DSRP). Regenerative processes with embedded regeneration cycles are widely used applications in operations research. In [7], regenerative processes have been applied to study priority queueing systems. In [9,10], the one-server queueing system with infinite queue, recurrent input flow, and general service time distribution ($GI|GI|1,\infty$—system in Kendall’s notation [11]) has been studied in terms of these processes (see also [12]). In paper [13], the authors investigated a model of a polling system by means of the DSRP. Another approach related to the study of complex systems (proposed in Belyaev [14]) contains the introduction of so-called supplementary variables. According to this method, the original model can be described by a continuous-time Markov chain (Markovization method). The implementation of this method to study reliability models is proposed by Rykov [15], Mazumdar [16], and Mine and Asakura [17].

The redundancy of some technical system means the duplication of critical components resulting in increasing reliability of the system. The double redundant systems are well known and were studied and applied under different assumptions like repair maintenance, preventive maintenance policy, exponential distributions for life and repair times, and so on; see, e.g., Gaver [18], Gnedenko et al. [19], Mine and Asakura [20], and Srinivasan [21]. In this paper, we study a homogeneous cold double redundant system with arbitrary life- and repair time distributions. In such a system, the cold standby redundancy means that the redundant component does not fail in a standby mode. The closed form solutions for the main reliability characteristics have been derived in terms of an integral transformation like a Laplace transform. The proposed analysis is based on the DSRP, the stochastic equations for the corresponding random variables (RVs), and the theory of regenerative processes [1].

The rest of the paper is organized as follows. In Section 2, the problem set, some assumptions, and notations are introduced. In Section 3, the Laplace–Stieltjes transforms (LSTs) of the regeneration period (the system life cycle), and the Laplace transforms (LTs) of the times to the first failure and between system failures are calculated. In Section 4 and Section 5, the time-dependent system state probabilities (TDSSPs) and stationary system state probabilities (SSSPs) are calculated. Finally, using some numerical examples, we demonstrate the effect of type of distributions and their parameter on the reliability and probabilistic measures of the system in Section 6. The paper ends with a conclusion.

2. The Problem Set: Assumptions and Notations

Consider a homogeneous cold double redundant renewable system with generally distributed life- and repair times, which, according to modified Kendall’s notations [11], will be denoted as <$G{I}_2/GI/1$>. The system consists of two identical units which can be in two possible states: operational and failed. The system fails when both of units are in a failed state. For the renewable system, different strategies of its renovation are possible. In this paper, we consider a strategy when a new system cycle begins at the point in time when one unit starts working while the other unit failed and is being repaired.

Denote by $A_i$ $(i=1,2,\dots )$, life times of the system units, by $B_i$ $(i=1,2,\dots )$, its repair times, and suppose that all these random variables (RVs) are mutually independent and identically distributed. Thus, denote by $A\left(t\right)=\mathbb{P}[A_i\le t]$ and $B\left(t\right)=\mathbb{P}[B_i\le t]$ the corresponding cumulative distribution functions (CDFs). Denote also by A and B the RVs with the same CDFs as $A_i$ and $B_i$, respectively. Suppose that the instantaneous failures and repairs are impossible and their mean times are finite:

$$A\left(0\right)=B\left(0\right)=0,a={\int }_0^{\infty }(1-A\left(x\right))dx<\infty ,b={\int }_0^{\infty }(1-B\left(x\right))dx<\infty .$$

Denote by E the set of system states. We assume that $E=\{i=0,1,2\}$, where i stands for the number of failed units, and we introduce a random process $J=\left\{J\right(t),t\ge 0\},$ where

$$J\left(t\right)=\mathrm{number} \mathrm{of} \mathrm{failed} \mathrm{units} \mathrm{at} \mathrm{time}t.$$

A simple state transition diagram is illustrated in Figure 1.

Denote by F the time between system failures, i.e., the period of time the system begins to operate again up to the next system failure. We remind that the system is again in an operational state when operation of one unit and repair of another one begin simultaneously. Denote by $F_1$ the time to first system failure (the time from the epoch when both units are in operational state up to the system failure) (see Figure 2). The functions $F\left(t\right)=\mathbb{P}[F\le t]$ and $F_1\left(t\right)=\mathbb{P}[F_1\le t]$ are corresponding CDFs.

We are interesting in calculation of the reliability function

$$R\left(t\right)=\mathbb{P}[F_1>t]=1-F_1\left(t\right),$$

distribution of time to the first failure $F_1\left(t\right)$, time-dependent system state probabilities (TDSSPs)

$${\pi }_j\left(t\right)=\mathbb{P}[J\left(t\right)=j],j=0,1,2,$$

and stationary system state probabilities (SSSPs)

$${\pi }_j=\underset{t\to \infty }{lim}{\pi }_j\left(t\right)\equiv \underset{t\to \infty }{lim}\mathbb{P}[J\left(t\right)=j],j=0,1,2,$$

as well as the availability coefficient

$$K_{\mathrm{av}.}={\pi }_0+{\pi }_1=1-{\pi }_2.$$

The following notations will be used next.

• The Laplace–Stiltjes transforms (LSTs) of their CDFs are denoted as:

  $$\begin{array}{ccc}\hfill \tilde{a}\left(s\right)& =& \mathbb{E}\left[e^{-sA}\right]=\underset{0}{\overset{\infty }{\int }}{e}^{-sx}dA\left(x\right),\tilde{b}\left(s\right)=\mathbb{E}\left[e^{-sB}\right]=\underset{0}{\overset{\infty }{\int }}{e}^{-sx}dB\left(x\right),\hfill \\ \hfill \tilde{f}\left(s\right)& =& \mathbb{E}\left[e^{-sF}\right]=\underset{0}{\overset{\infty }{\int }}{e}^{-sx}dF\left(x\right),{\tilde{f}}_1\left(s\right)=\mathbb{E}\left[e^{-s{F}^{\left(1\right)}}\right]=\underset{0}{\overset{\infty }{\int }}{e}^{-sx}d{F}_1\left(x\right).\hfill \end{array}$$
• In the case when the repair completion comes before the failure of other unit and when an operative unit fails before the repair completion of other unit, the distributions are equal, respectively,

  $$\begin{array}{c}\hfill a_B\left(t\right)={\int }_0^{t}B\left(x\right)dA\left(x\right)=\mathbb{P}[B\le A\le t],b_A\left(t\right)={\int }_0^{t}A\left(x\right)dB\left(x\right)=\mathbb{P}[A\le B\le t].\end{array}$$

  Applying the Laplace transform (LT) for $a_B\left(t\right)$ and $b_A\left(t\right)$, we get

  $$\begin{array}{cc}\hfill & {\tilde{a}}_B\left(s\right)={\int }_0^{\infty }{e}^{-st}d{a}_B\left(t\right)={\int }_0^{\infty }{e}^{-sx}B\left(x\right)dA\left(x\right),\hfill \\ \hfill & {\tilde{b}}_A\left(s\right)={\int }_0^{\infty }{e}^{-st}d{b}_A\left(t\right)={\int }_0^{\infty }{e}^{-sx}A\left(x\right)dB\left(x\right).\hfill \end{array}$$

  (1)
• Their derivatives with sign minus in point zero are denoted by:

  $$\begin{array}{ccc}\hfill a_B& =& -\frac{d}{ds}{\tilde{a}}_B\left(s\right){|}_{s=0}={\int }_0^{\infty }xB\left(x\right)dA\left(x\right),\hfill \\ \hfill b_A& =& -\frac{d}{ds}{\tilde{b}}_A\left(s\right){|}_{s=0}={\int }_0^{\infty }xA\left(x\right)dB\left(x\right).\hfill \end{array}$$

  (2)
• The probabilities $\mathbb{P}[B\le A]$ and $\mathbb{P}[B\ge A]$ are associated with transformations through the following ratios:

  $$\begin{array}{ccc}\hfill {\tilde{a}}_B\left(0\right)=\underset{0}{\overset{\infty }{\int }}B\left(x\right)dA\left(x\right)& =& \mathbb{P}[B\le A]\equiv p,\hfill \\ \hfill {\tilde{b}}_A\left(0\right)=\underset{0}{\overset{\infty }{\int }}A\left(x\right)dB\left(x\right)& =& \mathbb{P}[B>A]\equiv q=1-p.\hfill \end{array}$$
• Note the property of transformations (1)

  $${\tilde{a}}_{1-B}\left(s\right)=\tilde{a}\left(s\right)-{\tilde{a}}_B\left(s\right),{\tilde{b}}_{1-A}\left(s\right)=\tilde{b}\left(s\right)-{\tilde{b}}_A\left(s\right).$$

  (3)

3. Reliability Function

Process J is a regenerative one. A trajectory of this process is illustrated in Figure 2.

Here, F stands for the time between system failures, i.e., the period of time the system begins to operate again (state 1) up to the next system failure (state 2). The variable G specifies the length of a regenerative cycle according to a defined renovation policy. The following lemma holds for the LSTs of the time F between failures and the time to the first failure $F_1$.

Lemma 1.

The LST $\tilde{f}\left(s\right)=\mathbb{E}\left[e^{-sF}\right]$ of the time F between failures and the LST ${\tilde{f}}_1\left(s\right)=\mathbb{E}\left[e^{-s{F}_1}\right]$ of the time to the first failure $F_1$ are of the form:

$$\tilde{f}\left(s\right)=\frac{\tilde{a}\left(s\right)-{\tilde{a}}_B\left(s\right)}{1-{\tilde{a}}_B\left(s\right)},{\tilde{f}}_1\left(s\right)=\tilde{a}\left(s\right)\frac{\tilde{a}\left(s\right)-{\tilde{a}}_B\left(s\right)}{1-{\tilde{a}}_B\left(s\right)}.$$

(4)

Proof. 

The life time of the system F is a time between two successive failures. After a failure in state 0, the system goes to state 1 where there can be two events: either the recovery of a unit in state 1 in time B and subsequent transition again to state 0 or a failure of a second unit in time A will happen. After a failure in state 1, the system comes back to state 1 in time B. Hence, as it can be seen from Figure 2, the system lifetime F satisfies the following stochastic equation:

$$F=\left\{\begin{array}{cc}A+F\hfill & \mathrm{if}B<A,\hfill \\ A\hfill & \mathrm{if}B>A.\hfill \end{array}\right.$$

(5)

The time to the first system failure $F_1$ satisfies

$$F_1=A+F.$$

(6)

Applying LT to the Equation (5) and taking into account notations (1), one can obtain

$$\begin{array}{cc}\hfill \tilde{f}\left(s\right)& =\mathbb{E}\left[e^{-sF}\right]=\underset{0}{\overset{\infty }{\int }}{e}^{-st}dF\left(t\right)=\hfill \\ \hfill & =\underset{0}{\overset{\infty }{\int }}{e}^{-sx}\left[B\left(x\right)\tilde{f}\left(s\right)+(1-B\left(x\right))\right]dA\left(x\right)=\hfill \\ \hfill & =\tilde{f}\left(s\right){\tilde{a}}_B\left(s\right)+{\tilde{a}}_{1-B}\left(s\right).\hfill \end{array}$$

(7)

Due to the last equation, one can find that the LST of the system lifetime $\tilde{f}\left(s\right)$ is given by the first relation in of (4). The second equation of lemma directly follows from stochastic Equation (6). □

Remark 1.

Let ${\tilde{\phi }}_i\left(s\right)$, $i=0,1,2$, be the LST of the distribution of time to the first failure of the system, starting from state i at $t=0$. In previous statement ${\tilde{f}}_1\left(s\right)={\tilde{\phi }}_0\left(s\right)$. The following equations hold:

$$\begin{array}{cc}\hfill & {\tilde{\phi }}_0\left(s\right)={\tilde{q}}_{0,1}\left(s\right){\tilde{\phi }}_1\left(s\right),\hfill \\ \hfill & {\tilde{\phi }}_1\left(s\right)={\tilde{q}}_{1,1}\left(s\right){\tilde{\phi }}_1\left(s\right)+{\tilde{q}}_{1,2}\left(s\right).\hfill \end{array}$$

Here, ${\tilde{q}}_{i,j}\left(s\right)$ stands for the LST of the transition time distribution $q_{i,j}\left(t\right)$ from state i to state j. Obviously, we have ${\tilde{q}}_{0,1}\left(s\right)=\tilde{a}\left(s\right)$, ${\tilde{q}}_{1,1}\left(s\right)={\tilde{a}}_B\left(s\right)$ and ${\tilde{q}}_{1,2}\left(s\right)={\tilde{a}}_{1-B}\left(s\right)$. Expressing ${\tilde{\phi }}_0\left(s\right)$, we get the relation

$${\tilde{\phi }}_0\left(s\right)=\frac{{\tilde{q}}_{0,1}\left(s\right){\tilde{q}}_{1,2}\left(s\right)}{1-{\tilde{q}}_{1,1}\left(s\right)},$$

which coincide with the relation (4) for the first time to failure.

The main result of this section is the following theorem:

Theorem 1.

The LT ${\tilde{R}}_1\left(s\right)$ of the system reliability function $R_1\left(t\right)=1-F_1\left(t\right)$ is

$${\tilde{R}}_1\left(s\right)=\frac{(1-\tilde{a}\left(s\right))(1+\tilde{a}\left(s\right)-{\tilde{a}}_B\left(s\right))}{s(1-{\tilde{a}}_B\left(s\right))}.$$

(8)

Proof. 

Taking into account that the LT of any CDF connected with its LST through a relation

$${\tilde{F}}_1\left(s\right)=\underset{0}{\overset{\infty }{\int }}{e}^{-st}{F}_1\left(t\right)dt=\frac{1}{s}\underset{0}{\overset{\infty }{\int }}{e}^{-st}d{F}_1\left(t\right)=\frac{1}{s}{\tilde{f}}_1\left(s\right),$$

the LT ${\tilde{R}}_1\left(s\right)$ of the reliability function $R_1\left(t\right)=1-F_1\left(t\right)$ can be derived using (4) in form

$$\begin{array}{cc}\hfill {\tilde{R}}_1\left(s\right)& =\frac{1}{s}-{\tilde{F}}_1\left(s\right)=\frac{1}{s}(1-{\tilde{f}}_1\left(s\right))=\frac{1-{\tilde{a}}_B\left(s\right)-{\tilde{a}}^2\left(s\right)+\tilde{a}\left(s\right){\tilde{a}}_B\left(s\right)}{s(1-{\tilde{a}}_B\left(s\right))}\hfill \\ \hfill & =\frac{(1-\tilde{a}\left(s\right))(1+\tilde{a}\left(s\right)-{\tilde{a}}_B\left(s\right))}{s(1-{\tilde{a}}_B\left(s\right))},\hfill \end{array}$$

which ends the proof. □

Corollary 1.

The LT $\tilde{R}\left(s\right)$ of the survival probability $R\left(t\right)=1-F\left(t\right)$ is of the form

$$\begin{array}{c}\hfill \tilde{R}\left(s\right)=\frac{\tilde{R}\left(s\right)}{1+\tilde{a}\left(s\right)-{\tilde{a}}_B\left(s\right)}=\frac{1-\tilde{a}\left(s\right)}{s(1-{\tilde{a}}_B\left(s\right))},\end{array}$$

(9)

which follows directly from (4) and (8).

Corollary 2.

The mean system life time between failures $\mathbb{E}\left[F\right]$ and mean time to the first failure $\mathbb{E}\left[F_1\right]$ are of the form:

$$\mathbb{E}\left[F\right]=\frac{a}{q},\mathbb{E}\left[F_1\right]=a+\frac{a}{q}.$$

(10)

Proof. 

The first relation follows from the equality

$$\mathbb{E}\left[F\right]=\tilde{R}\left(0\right)=\underset{s\to 0}{lim}\frac{1-\tilde{a}\left(s\right)}{s(1-{\tilde{a}}_B\left(s\right))}=\frac{a}{1-{\tilde{a}}_B\left(0\right)}=\frac{a}{q},$$

and the second one follows directly from the relation (6). □

4. Time Dependent System State Probabilities

4.1. Definition

For calculation of the time dependent system state probabilities (TDSSP), we use a renewal theory. In our case, the process J is a regenerative process with a delay (see Figure 2) and its regeneration times are

$$S_0=0,S_1=A_1,S_2=S_1+G_1,\dots ,S_{k+1}=S_k+G_k,\dots .$$

Here, regeneration cycles $G_i(i=1,2,\dots )$ are the time intervals between two successive returns of the process J into state 1 after a system failure, when one of the system units begins to operate and another one is being repaired. Thus, if the system functioning begins when both units are in operational state, the first regeneration period differs from the others and has the duration $A_1$. The state probabilities of the process J in any time t

$${\pi }_j\left(t\right)=\mathbb{P}\{J\left(t\right)=j\},j=0,1,2$$

can be represented in terms of its distribution over separate regeneration period

$${\pi }_j^{\left(1\right)}\left(t\right)=\mathbb{P}[J\left(t\right)=j,t<G],j=0,1,2,$$

(11)

and the renewal function $H\left(t\right)$ can be defined as follows:

$${\pi }_i\left(t\right)={\pi }_i^{(1,1)}\left(t\right)+\underset{0}{\overset{t}{\int }}dH\left(u\right){\pi }_j^{\left(1\right)}(t-u).$$

(12)

The distribution of the process J on this period is of the form

$${\pi }_j^{(1,1)}\left(t\right)=\mathbb{P}[J\left(t\right)=j,t<A_1]={\delta }_{j0}(1-A\left(t\right)).$$

(13)

In terms of the CDF $G\left(t\right)=\mathbb{P}[G_i\le t]$ of RVs $G_i$, a renewal function $H\left(t\right)$ is determined as follows:

$$H\left(t\right)=\sum _{k\ge 1}\mathbb{P}\left[\left(\sum _{1\le i\le k}{G}_i\right)\le t\right]=\sum _{k\ge 1}{G}^{*k}\left(t\right).$$

(14)

On the basis of the foregoing, the next lemma can be proved. Consider, first of all, the regeneration cycle distribution.

Lemma 2.

The LST of the regeneration cycle is of the form:

$$\tilde{g}\left(s\right)=\mathbb{E}\left[e^{-sG}\right]=\frac{{\tilde{b}}_A\left(s\right)}{1-{\tilde{a}}_B\left(s\right)}.$$

(15)

Proof. 

The regeneration cycle is a time between two successive visits of state 1 where two events can occur: either the recovery of a unit in state 1 in time B and subsequent transition to state 0 or a failure of a second unit in time A happens. Figure 2 shows that RV G satisfies the following stochastic equation due to independence of RVs A and B:

$$G=\left\{\begin{array}{cc}A+G\hfill & \mathrm{if}A>B,\hfill \\ B\hfill & \mathrm{if}A\le B.\hfill \end{array}\right.$$

(16)

Applying LST to this stochastic equation, we obtain

$$\begin{array}{ccc}\hfill \tilde{g}\left(s\right)& =\hfill & \mathbb{E}\left[e^{-sG}\right]=\underset{0}{\overset{\infty }{\int }}{e}^{-st}dG\left(t\right)\hfill \\ & =\hfill & \underset{0}{\overset{\infty }{\int }}dA\left(x\right)\left[B\left(x\right)e^{-sx}\tilde{g}\left(s\right)+{\int }_{y>x}{e}^{-sy}dB\left(y\right)\right]=\hfill \\ & =\hfill & \tilde{g}\left(s\right)\underset{0}{\overset{\infty }{\int }}{e}^{-sx}B\left(x\right)dA\left(x\right)+{\int }_0^{\infty }{e}^{-sy}A\left(y\right)dB\left(y\right)=\hfill \\ & =\hfill & \tilde{g}\left(s\right){\tilde{a}}_B\left(s\right)+{\tilde{b}}_A\left(s\right),\hfill \end{array}$$

(17)

which implies the expression (15) for the LST of the regeneration period. □

As a corollary, we obtain the mean length of the regeneration period.

Corollary 3.

The mean length of the regeneration period can be defined as:

$$\mathbb{E}\left[G\right]=\frac{a_B+b_A}{q}.$$

(18)

Proof. 

The proof is obtained by differentiation

$$\mathbb{E}\left[G\right]=-\frac{d}{ds}\tilde{g}\left(s\right){|}_{s=0}=-\frac{{\tilde{b}}_A^{\prime }\left(0\right)(1-{\tilde{a}}_B\left(0\right))+{\tilde{a}}_B^{\prime }{\tilde{b}}_A\left(0\right)}{{(1-{\tilde{a}}_B\left(0\right))}^2}.$$

This expression leads to expression (18) due to properties of ${\tilde{a}}_B\left(0\right)$ and ${\tilde{b}}_A\left(0\right)$. □

Lemma 3.

The LST of the system renewal function is given by:

$$\tilde{h}\left(s\right)=\frac{{\tilde{b}}_A\left(s\right)}{1-({\tilde{a}}_B\left(s\right)+{\tilde{b}}_A\left(s\right))}.$$

(19)

Proof. 

From the renewal theory, it is well known (and follows from (14)) that the LST of the renewal function $H\left(t\right)$ is defined as $\tilde{h}\left(s\right)=\tilde{g}\left(s\right){(1-\tilde{g}\left(s\right))}^{-1}$. Thus, substitution to this expression of (15) for $\tilde{g}\left(s\right)$ leads to (19). □

Theorem 2.

The LTs ${\tilde{\pi }}_j\left(s\right)$ of the TDSSPs ${\pi }_j\left(t\right)$ $(j=0,1,2)$ are of the form:

$${\tilde{\pi }}_j\left(s\right)={\delta }_{j0}\frac{1-\tilde{a}\left(s\right)}{s}+\frac{{\tilde{b}}_a\left(s\right)}{1-({\tilde{a}}_B\left(s\right)+{\tilde{b}}_A\left(s\right))}{\tilde{\pi }}_j^{\left(1\right)}\left(s\right),j=0,1,2,$$

(20)

where ${\tilde{\pi }}_j^{\left(1\right)}\left(s\right)(j=0,1,2)$ are the LTs of the TDSSPs ${\pi }_j^{\left(1\right)}\left(t\right)$ $(j=0,1,2)$ in a separate regeneration period. These probabilities will be calculated in the next subsection.

Proof. 

The form of Equations (12) and (14) prompts to calculate the process state probabilities in terms of their LTs and the renewal function LST

$${\tilde{\pi }}_j\left(s\right)={\int }_0^{\infty }{e}^{-st}{\pi }_j\left(t\right)dt,\tilde{h}\left(s\right)={\int }_0^{\infty }{e}^{-st}dH\left(t\right).$$

Applying LT to Equation (12) and taking into account Equation (13), one can obtain

$${\tilde{\pi }}_j\left(s\right)={\delta }_{j0}\frac{1-\tilde{a}\left(s\right)}{s}+\tilde{h}\left(s\right){\tilde{\pi }}_j^{\left(1\right)}\left(s\right).$$

(21)

Substitution into this equality of the expression for $\tilde{h}\left(s\right)$ (19) leads to (20) and completes the proof. □

4.2. The State Probabilities in a Separate Regeneration Period

Now, we calculate the TDSSPs in a separate regeneration period. The probability ${\pi }_2^{\left(1\right)}\left(t\right)$ can be calculated very easily in the main level regeneration period.

Lemma 4.

The LT of the second TDSSP ${\pi }_2^{\left(1\right)}\left(t\right)$ in the main regeneration period is given by

$${\tilde{\pi }}_2^{\left(1\right)}\left(s\right)=\frac{\tilde{a}\left(s\right)-({\tilde{a}}_B\left(s\right)+{\tilde{b}}_A\left(s\right))}{s(1-{\tilde{a}}_B\left(s\right))}.$$

(22)

Proof. 

Due to representation of the regeneration period by Formula (16), the event $\left\{J\right(t)=2,t<G\}$ occurs if and only if either the event $\{A_1\le t\le B_1\}$ occurs or the events $\{t>u=A_1>B_1\}$ and $\left\{J\right(t-u)=2,t-u<G\}$ occur. Thus, as it is shown in Figure 2, for the probability ${\pi }_2^{\left(1\right)}\left(t\right)$, it holds that

$${\pi }_2^{\left(1\right)}\left(t\right)=\mathbb{P}[J\left(t\right)=2,t<G]=\mathbb{P}[A\le t<B]+\underset{0}{\overset{t}{\int }}dA\left(u\right)B\left(u\right){\pi }_2^{\left(1\right)}(t-u).$$

From this equation for the LST, it follows

$${\tilde{\pi }}_2^{\left(1\right)}\left(s\right)=\underset{0}{\overset{\infty }{\int }}{e}^{-st}A\left(t\right)(1-B\left(t\right))dt+{\tilde{a}}_B\left(s\right){\pi }_2^{\left(1\right)}\left(s\right)$$

and

$${\tilde{\pi }}_2^{\left(1\right)}\left(s\right)=\frac{1}{1-{\tilde{a}}_B\left(s\right)}\underset{0}{\overset{\infty }{\int }}{e}^{-st}A\left(t\right)(1-B\left(t\right))dt.$$

Calculating LTs of the integral in these expressions using partial integration, we obtain

$$\begin{array}{ccc}\hfill {\tilde{\pi }}_2^{\left(1\right)}\left(s\right)& =\hfill & \frac{1}{1-{\tilde{a}}_B\left(s\right)}\underset{0}{\overset{\infty }{\int }}{e}^{-st}A\left(t\right)(1-B\left(t\right))dt=\hfill \\ & =\hfill & -\frac{1}{s(1-{\tilde{a}}_B\left(s\right))}{e}^{-st}A\left(t\right)(1-B\left(t\right))dt{|}_0^{\infty }+\hfill \\ & +\hfill & \frac{1}{s(1-{\tilde{a}}_B\left(s\right))}\underset{0}{\overset{\infty }{\int }}{e}^{-st}[dA\left(t\right)(1-B\left(t\right))-A\left(t\right)dB\left(t\right)]=\hfill \\ & =\hfill & \frac{\tilde{a}\left(s\right)-({\tilde{a}}_B\left(s\right)+{\tilde{b}}_A\left(s\right))}{s(1-{\tilde{a}}_B\left(s\right))}\hfill \end{array}$$

that completes the proof. □

Since the calculation of the probabilities ${\pi }_j^{\left(1\right)}(j=0,1)$ is not a trivial task, we intend to apply the theory of DSRP [3,8,12]. According to this theory, the further analysis can be described as follows. As it was shown in the Equation (12), the process distribution can be represented in terms of its distribution in embedded regeneration cycles and embedded renewal function. This approach is useful if the behavior of the process is too complex in order to derive a state distribution directly, but the embedded regeneration epochs are available. At known regeneration epochs $S_k^{\left(1\right)}$, the distribution in the basic regeneration period (regeneration period of the first level) ${\pi }_i^{\left(1\right)}\left(t\right)$ similar to the Equation (12) can be represented in terms of distributions in embedded regeneration cycles (regeneration cycles of the second level) ${\pi }_j^{\left(2\right)}\left(t\right)$ and embedded renewal function $H^{\left(1\right)}\left(t\right)$ in the following way:

$${\pi }_j^{\left(1\right)}\left(t\right)={\pi }_j^{\left(2\right)}\left(t\right)+\underset{0}{\overset{\infty }{\int }}d{H}^{\left(1\right)}\left(u\right){\pi }_j^{\left(2\right)}(t-u),j=0,1,$$

(23)

where

$${\pi }_j^{\left(2\right)}\left(t\right)=\mathbb{P}\{J\left(t\right)=j,t<G^{\left(2\right)}\},j=0,1$$

is the process state distribution in a separate regeneration period of the second level. In this case as embedded regeneration epochs, we choose the random number $\nu =min\{n:A_n<B_n\}$ of time epochs

$$S_1^{\left(1\right)}=A_1{1}_{\{A_1>B_1\}},S_2^{\left(1\right)}=S_1^{\left(1\right)}+A_2{1}_{\{A_1>B_1,A_2>B_2\}},\dots$$

up to time, when the event $\{A_n\le B_n\}$ occurs for the first time. It means that the time interval $G^{\left(1\right)}$ between embedded regeneration epochs has a distribution $G^{\left(1\right)}\left(t\right)=A\left(t\right)$ and these epochs in the time interval F are determined by Equation (5). According to the theory of DSRP, the embedded renewal function $H^{\left(1\right)}\left(t\right)$ satisfies the equation

$$H^{\left(1\right)}\left(t\right)=A\left(t\right)+\underset{0}{\overset{t}{\int }}d{H}^{\left(1\right)}\left(u\right)A(t-u)-F\left(t\right),$$

(24)

where $F\left(t\right)$ is a CDF of the time between system failures determined by its LST (4).

Similar to the basic case, the solution of Equations (23) and (24) can be represented in terms of their LTs and LSTs:

$${\tilde{\pi }}_j^{\left(2\right)}\left(s\right)=\underset{0}{\overset{\infty }{\int }}{e}^{-st}{\pi }_j^{\left(2\right)}\left(t\right)dt,{\tilde{h}}^{\left(1\right)}\left(s\right)=\underset{0}{\overset{\infty }{\int }}{e}^{-st}d{H}^{\left(1\right)}\left(t\right).$$

The next lemma specifies connections between process distributions in the first and in the second level regeneration cycles in terms of their LTs.

Lemma 5.

The TDSSPs in the first and in the second level regeneration cycles in terms of the LTs satisfy the relation

$${\tilde{\pi }}_j^{\left(1\right)}\left(s\right)={\tilde{\pi }}_j^{\left(2\right)}\left(s\right)\frac{{\tilde{a}}_B\left(s\right)}{1-{\tilde{a}}_B\left(s\right)}j=0,1.$$

(25)

Proof. 

Applying LT to Equation (23), we get

$${\tilde{\pi }}_j^{\left(1\right)}\left(s\right)=(1+{\tilde{h}}^{\left(1\right)}\left(s\right)){\tilde{\pi }}_j^{\left(2\right)}\left(s\right),$$

(26)

where, due to (24), the LST ${\tilde{h}}^{\left(1\right)}\left(s\right)$ of the embedded renewal function $H^{\left(1\right)}\left(t\right)$ is of the form

$${\tilde{h}}^{\left(1\right)}\left(s\right)=\tilde{a}\left(s\right)+{\tilde{h}}^{\left(1\right)}\left(s\right)\tilde{a}\left(s\right)-\tilde{f}\left(s\right),$$

which leads in turn to

$${\tilde{h}}^{\left(1\right)}\left(s\right)=\frac{\tilde{a}\left(s\right)-\tilde{f}\left(s\right)}{1-\tilde{a}\left(s\right)}.$$

Substitution of this relation to the Equation (26) and taking into account expressions for $\tilde{f}\left(s\right)$ from (4) implies the result (22) that completes the proof. □

We have to calculate now only ${\tilde{\pi }}_j^{\left(2\right)}\left(s\right)(j=0,1)$.

Lemma 6.

In notations (4), the LTs of the second level system state probabilities are:

$$\begin{array}{ccc}\hfill {\tilde{\pi }}_0^{\left(2\right)}\left(s\right)& =& \frac{1}{s}(\tilde{a}\left(s\right)-({\tilde{a}}_B\left(s\right)+{\tilde{b}}_A\left(s\right)));\hfill \\ \hfill {\tilde{\pi }}_1^{\left(2\right)}\left(s\right)& =& \frac{1}{s}[1-(\tilde{a}\left(s\right)+\tilde{b}\left(s\right))+({\tilde{a}}_B\left(s\right)+{\tilde{b}}_A\left(s\right))].\hfill \end{array}$$

(27)

Proof. 

First, consider appropriate expressions for ${\tilde{\pi }}_j^{\left(2\right)}\left(t\right)(j=0,1)$. From Figure 2, it follows that

• The event $\{J\left(t\right)=0,t<G^{\left(1\right)}\}$ occurs if and only if $\{B\le t<A\}$;
• The event $\{J\left(t\right)=1,t<G^{\left(1\right)}\}$ occurs if and only if $\{t<B\le A\}$ or $\{t<A\le B\}$;

Hence, the appropriate probabilities are

$$\begin{array}{ccc}\hfill {\pi }_0^{\left(2\right)}\left(t\right)& =\hfill & \mathbb{P}[B<t<A]=B\left(t\right)(1-A(t\left)\right),\hfill \\ \hfill {\pi }_1^{\left(2\right)}\left(t\right)& =\hfill & \mathbb{P}[t<B<A]+\mathbb{P}[t<A<B]=\hfill \\ & =\hfill & \underset{t}{\overset{\infty }{\int }}(1-A\left(u\right))dB\left(u\right)+\underset{t}{\overset{\infty }{\int }}(1-B\left(u\right))dA\left(u\right).\hfill \end{array}$$

(28)

Calculating LTs of these expressions by partial integration in terms of notations (1) gives

$$\begin{array}{ccc}\hfill {\tilde{\pi }}_0^{\left(2\right)}\left(s\right)& =\hfill & \underset{0}{\overset{\infty }{\int }}{e}^{-st}B\left(t\right)(1-A\left(t\right))dt=\hfill \\ & =\hfill & \frac{1}{s}\underset{0}{\overset{\infty }{\int }}{e}^{-st}[dB\left(t\right)(1-A\left(t\right))-B\left(t\right)dA\left(t\right)]=\hfill \\ & =\hfill & \frac{1}{s}(\tilde{b}\left(s\right)-({\tilde{a}}_B\left(s\right)+{\tilde{b}}_A\left(s\right)));\hfill \\ \hfill {\tilde{\pi }}_1^{\left(2\right)}\left(s\right)& =\hfill & \underset{0}{\overset{\infty }{\int }}{e}^{-st}\underset{t}{\overset{\infty }{\int }}\left[(1-A(u\left)\right)dB\left(u\right)+(1-B(u\left)\right)dA\left(u\right)\right]=\hfill \\ & =\hfill & \underset{0}{\overset{\infty }{\int }}\underset{0}{\overset{u}{\int }}{e}^{-st}dt\left[(1-A(u\left)\right)dB\left(u\right)+(1-B(u\left)\right)dA\left(u\right)\right]=\hfill \\ & =\hfill & \frac{1}{s}\underset{0}{\overset{\infty }{\int }}(1-e^{-su})\left[(1-A(u\left)\right)dB\left(u\right)+(1-B(u\left)\right)dA\left(u\right)\right]=\hfill \\ & =\hfill & \frac{1}{s}[1-(\tilde{a}\left(s\right)+\tilde{b}\left(s\right))+({\tilde{a}}_B\left(s\right)+{\tilde{b}}_A\left(s\right))]\hfill \end{array}$$

that ends the proof. □

By combining all obtained results for LTs of the TDSSPs, one can get the following statement.

Theorem 3.

The LTs of the TDSSPs are of the form:

$$\begin{array}{ccc}\hfill {\tilde{\pi }}_0\left(s\right)& =\hfill & \frac{1-\tilde{a}\left(s\right)}{s}+\frac{{\tilde{b}}_A\left(s\right)}{s(1-{\tilde{a}}_B\left(s\right))}\times \frac{\tilde{b}\left(s\right)-({\tilde{a}}_B\left(s\right)+{\tilde{b}}_A\left(s\right))}{1-({\tilde{a}}_B\left(s\right)+{\tilde{b}}_A\left(s\right))};\hfill \\ \hfill {\tilde{\pi }}_1\left(s\right)& =\hfill & \frac{{\tilde{b}}_A\left(s\right)}{s(1-{\tilde{a}}_B\left(s\right))}\times \frac{1-(\tilde{a}\left(s\right)+\tilde{b}\left(s\right))+({\tilde{a}}_B\left(s\right)+{\tilde{b}}_A\left(s\right))}{1-({\tilde{a}}_B\left(s\right)+{\tilde{b}}_A\left(s\right))};\hfill \\ \hfill {\pi }_2\left(s\right)& =\hfill & \frac{{\tilde{b}}_A\left(s\right)}{s(1-{\tilde{a}}_B\left(s\right))}\times \frac{\tilde{a}\left(s\right)-({\tilde{a}}_B\left(s\right)+{\tilde{b}}_A\left(s\right))}{1-({\tilde{a}}_B\left(s\right)+{\tilde{b}}_A\left(s\right))}.\hfill \end{array}$$

(29)

Proof. 

For $j=0,1$, the proof is presented as a series of equalities:

$$\begin{array}{ccc}\hfill {\tilde{\pi }}_j\left(s\right)& =\hfill & {\delta }_{j0}\frac{1-\tilde{a}\left(s\right)}{s}+\tilde{h}\left(s\right){\pi }_j^{\left(1\right)}\left(s\right)=\hfill \\ & =\hfill & {\delta }_{j0}\frac{1-\tilde{a}\left(s\right)}{s}+\frac{\tilde{g}\left(s\right)}{1-\tilde{g}\left(s\right)}(1+{\tilde{h}}_j^{\left(1\right)}\left(s\right)){\pi }_j^{\left(2\right)}\left(s\right)=\hfill \\ & =\hfill & {\delta }_{j0}\frac{1-\tilde{a}\left(s\right)}{s}+\frac{{\tilde{b}}_A\left(s\right)}{1-({\tilde{a}}_B\left(s\right)+{\tilde{b}}_A\left(s\right))}\frac{\tilde{a}\left(s\right)-\tilde{f}\left(s\right)}{1-\tilde{a}\left(s\right)}{\pi }_j^{\left(2\right)}\left(s\right)=\hfill \\ & =\hfill & {\delta }_{j0}\frac{1-\tilde{a}\left(s\right)}{s}+\frac{{\tilde{b}}_A\left(s\right)}{1-({\tilde{a}}_B\left(s\right)+{\tilde{b}}_A\left(s\right))}\frac{1}{1-{\tilde{a}}_B\left(s\right)}{\pi }_j^{\left(2\right)}\left(s\right).\hfill \end{array}$$

Substitution of expressions (27) for ${\tilde{\pi }}_j^{\left(2\right)}\left(s\right)$ leads to the first two formulas from (29). For the last one, it holds that

$${\tilde{\pi }}_2\left(s\right)=\tilde{h}\left(s\right){\tilde{\pi }}_2^{\left(1\right)}\left(s\right)=\frac{\tilde{g}\left(s\right)}{1-\tilde{g}\left(s\right)}{\tilde{\pi }}_2^{\left(1\right)}\left(s\right)=\frac{{\tilde{b}}_A\left(s\right)}{1-({\tilde{a}}_B\left(s\right)+{\tilde{b}}_A\left(s\right))}{\tilde{\pi }}_2^{\left(1\right)}\left(s\right).$$

Substitution into this formula of expression (22) for ${\tilde{\pi }}_2^{\left(1\right)}\left(s\right)$ leads to the last formula from (29). □

5. Stationary State Probabilities

To calculate the stationary state probabilities (SSSPs), we use a Tauber theorem. Due to (29), we get

$${\pi }_j=\underset{t\to \infty }{lim}{\pi }_j\left(t\right)=\underset{s\to 0}{lim}s{\tilde{\pi }}_j\left(s\right).$$

(30)

By means of this formula, we can prove the next statement.

Theorem 4.

The system state stationary probabilities are:

$$\begin{array}{ccc}\hfill {\pi }_0& =& 1-\frac{b}{a_B+b_A},\hfill \\ \hfill {\pi }_1& =& \frac{a+b}{a_B+b_A}-1,\hfill \\ \hfill {\pi }_2& =& 1-\frac{a}{a_B+b_A}.\hfill \end{array}$$

(31)

Proof. 

According to the limit relation (30), we obtain

$$\begin{array}{ccc}\hfill {\pi }_0& =\hfill & \underset{s\to 0}{lim}\left[1-\tilde{a}\left(s\right)+\frac{{\tilde{b}}_A\left(s\right)}{1-{\tilde{a}}_B\left(s\right)}\times \frac{\tilde{b}\left(s\right)-({\tilde{a}}_B\left(s\right)+{\tilde{b}}_A\left(s\right))}{1-({\tilde{a}}_B\left(s\right)+{\tilde{b}}_A\left(s\right))}\right]=\hfill \\ & =\hfill & \underset{s\to 0}{lim}\frac{\tilde{b}\left(s\right)-({\tilde{a}}_B\left(s\right)+{\tilde{b}}_A\left(s\right))}{1-({\tilde{a}}_B\left(s\right)+{\tilde{b}}_A\left(s\right))};\hfill \\ \hfill {\pi }_1& =\hfill & \underset{s\to 0}{lim}\left[\frac{{\tilde{b}}_A\left(s\right)}{1-{\tilde{a}}_B\left(s\right)}\times \frac{1-(\tilde{a}\left(s\right)+\tilde{b}\left(s\right))+({\tilde{a}}_B\left(s\right)+{\tilde{b}}_A\left(s\right))}{1-({\tilde{a}}_B\left(s\right)+{\tilde{b}}_A\left(s\right))}\right];\hfill \\ & =\hfill & \underset{s\to 0}{lim}\frac{1-(\tilde{a}\left(s\right)+\tilde{b}\left(s\right))+({\tilde{a}}_B\left(s\right)+{\tilde{b}}_A\left(s\right))}{1-({\tilde{a}}_B\left(s\right)+{\tilde{b}}_A\left(s\right))};\hfill \\ \hfill {\pi }_2& =\hfill & \underset{s\to 0}{lim}\frac{{\tilde{b}}_A\left(s\right)}{1-{\tilde{a}}_B\left(s\right)}\times \frac{\tilde{a}\left(s\right)-({\tilde{a}}_B\left(s\right)+{\tilde{b}}_A\left(s\right))}{1-({\tilde{a}}_B\left(s\right)+{\tilde{b}}_A\left(s\right))}=\hfill \\ & =\hfill & \underset{s\to 0}{lim}\frac{\tilde{a}\left(s\right)-({\tilde{a}}_B\left(s\right)+{\tilde{b}}_A\left(s\right))}{1-({\tilde{a}}_B\left(s\right)+{\tilde{b}}_A\left(s\right))}.\hfill \end{array}$$

Using L’Hopital’s rule, we get the next

$$\begin{array}{ccc}\hfill {\pi }_0& =\hfill & \frac{{\tilde{b}}^{\prime }\left(0\right)-({\tilde{a}}_B^{\prime }\left(0\right)+{\tilde{b}}_A^{\prime }\left(0\right))}{-({\tilde{a}}_B^{\prime }\left(0\right)+{\tilde{b}}_A^{\prime }\left(0\right))}=1-\frac{b}{a_B+b_A};\hfill \\ \hfill {\pi }_1& =\hfill & \frac{-(\tilde{a}\left(0\right)+{\tilde{b}}^{\prime }\left(0\right))+({\tilde{a}}_B^{\prime }\left(0\right)+{\tilde{b}}_A^{\prime }\left(0\right))}{-({\tilde{a}}_B^{\prime }\left(0\right)+{\tilde{b}}_A^{\prime }\left(0\right))}=\frac{a+b}{a_B+b_A}-1;\hfill \\ \hfill {\pi }_2& =\hfill & \frac{{\tilde{a}}^{\prime }\left(0\right)-({\tilde{a}}_B^{\prime }\left(0\right)+{\tilde{b}}_A^{\prime }\left(0\right))}{-({\tilde{a}}_B^{\prime }\left(0\right)+{\tilde{b}}_A^{\prime }\left(0\right))}=1-\frac{a}{a_B+b_A}\hfill \end{array}$$

that completes the proof. □

Example 1.

As an example, consider a Markov model $<M_2|M|1>$, when

$$A\left(t\right)=1-e^{-\alpha t},B\left(t\right)=1-e^{-\beta t}.$$

For the SSSPs, we get in terms of parameter $\rho =\frac{\alpha }{\beta }$:

$$a_B+b_A=\left[\frac{1}{\alpha }+\frac{1}{\beta }-\frac{1}{\alpha +\beta }\right]=\frac{{\alpha }^2+\alpha \beta +{\beta }^2}{\alpha \beta (\alpha +\beta )}$$

Now, for ${\pi }_j(j=0,1,2)$, it holds that

$$\begin{array}{ccc}\hfill {\pi }_0& =\hfill & 1-\frac{b}{a_b+b_A}=1-\frac{\alpha (\alpha +\beta )}{{\alpha }^2+\alpha \beta +{\beta }^2}=\frac{1}{1+\rho +{\rho }^2};\hfill \\ \hfill {\pi }_1& =\hfill & \frac{a+b}{a_B+b_A}-1=\frac{\alpha \beta }{{\alpha }^2+\alpha \beta +{\beta }^2}=\frac{\rho }{1+\rho +{\rho }^2};\hfill \\ \hfill {\pi }_2& =\hfill & 1-\frac{a}{a_B+b_A}=\frac{{\alpha }^2}{{\alpha }^2+\alpha \beta +{\beta }^2}=\frac{{\rho }^2}{1+\rho +{\rho }^2}.\hfill \end{array}$$

This result coincides with those calculated by direct approach using Birth and Death process for the Markov case.

6. Numerical Analysis

In this section, we study the sensitivity of reliability function, TDSSPs and SSSPs to types of distribution functions and their parameters. The life- and repair-times will follow exponential, PH-type, Pareto and log-normal distributions. For valid comparison of the results, we choose the parameters of these distributions in such a way that the mean and the variance of the life and repair times always coincide with an exponential version, i.e., $\mathbb{E}\left[A\right]=a=\frac{1}{\alpha }$ and $\mathbb{E}\left[B\right]=b=\frac{1}{\beta }$. Now, we give a short description of the considered life-time distributions $A\left(t\right)$ with appropriately selected parameters. The repair-time distributions $B\left(t\right)$ are defined by analogy.

The PH-type distribution $\mathcal{PH}$ with representation $(\eta ,M)$,

$$A\left(t\right)=1-\eta e^{Mt}\mathbf{1},\mathbb{E}\left[A\right]=a=-\eta M^{-1}\mathbf{1}=\frac{1}{\alpha },\mathbb{V}\left[A\right]=2\eta M^{-2}\mathbf{1}-a^2=\frac{1}{{\alpha }^2},$$

where $\eta =(0.7,0.1,0.2)$ and

$$\begin{array}{cc}\hfill & M=\left(\begin{array}{ccc}-0.80& 0.19& 0.50\\ 3.35& -3.50& 0.10\\ 0.80& 0.60& -1.50\end{array}\right)\mathrm{for}a=10,M=\left(\begin{array}{ccc}-2.80& 1.14& 1.00\\ 13.39& -15.00& 1.50\\ 2.80& 1.60& -4.50\end{array}\right)\mathrm{for}a=2,\hfill \\ \hfill & M=\left(\begin{array}{ccc}-8.00& 1.94& 5.00\\ 33.48& -35.00& 1.00\\ 8.00& 6.00& -15.00\end{array}\right)\mathrm{for}a=1,M=\left(\begin{array}{ccc}16.00& 1.03& 10.00\\ 25.43& -35.00& 2.00\\ 16.00& 10.00& -30.00\end{array}\right)\mathrm{for}a=0.2,\hfill \\ \hfill & M=\left(\begin{array}{ccc}-35.00& 18.81& 5.00\\ 94.70& -110.00& 10.00\\ 10.00& 60.00& -80.00\end{array}\right)\mathrm{for}b=0.1.\hfill \end{array}$$

The Pareto distribution $\mathcal{PAR}$ with parameters $(k,t_{\min })$,

$$A\left(t\right)=1-{\left(\frac{t_{\min }}{t}\right)}^k,t\ge t_{\min },k=1+\sqrt{2},t_{\min }=\frac{\sqrt{2}}{\alpha (1+\sqrt{2})}.$$

The log-normal distribution $\mathcal{LN}$ with parameters $(\mu ,\sigma )$,

$$A\left(t\right)=\Phi \left(\frac{ln\left(x\right)-\mu }{\sigma }\right),\mu =ln\frac{1}{\alpha }-\frac{ln\left(2\right)}{2},\sigma =\sqrt{ln\left(2\right)}.$$

Now, we fix the parameters of the repair-time distributions to guarantee that $\mathbb{E}\left[B\right]=b=0.1$ and $\mathbb{V}\left[B\right]=0.01$. In Figure 3, we demonstrate the reliability $R_1\left(t\right)$ and survival $R\left(t\right)$ functions obtained by numerical inversion with Euler algorithm of the corresponding LTs (8) and (9) in case $A\sim \mathcal{PH}({\eta }_1,M_1)$ and $B\sim \mathcal{PH}({\eta }_2,M_2)$. Obviously, for $a\ge 1$, we observe that the functions exhibit a heavier tail. Note that these graphics are indistinguishable from those obtained for exponential distributions.

In case it is not feasible to evaluate Laplace–Stieltjes transforms in explicit form (4), different approaches can be used instead to recover the distribution function via moments. We have tried the maximum entropy solution and simple fitting by a mixture of exponential distributions. Both of these methods work and give reasonable approximations. For details about maximum entropy solution, the reader is referred to Falin et al. [22]. To illustrate this methodology, we limit ourselves to analyzing the length of the system life time between failures F, when the first, e.g., four moments $\mathbb{E}\left[F^k\right]={(-1)}^k\frac{d^k}{d{s}^k}\tilde{f}\left(s\right){|}_{s=0}$, $1\le k\le m=4$, are available,

$$\begin{array}{cc}\hfill & \mathbb{E}\left[F\right]=\frac{a}{q},\hfill \\ \hfill & \mathbb{E}\left[F^2\right]=\frac{2a·a_B+q{a}^{\left(2\right)}}{q^2},\hfill \\ \hfill & \mathbb{E}\left[F^3\right]=\frac{3a(q{a}_B^{\left(2\right)}+2a_B^2)+q(3a^{\left(2\right)}{a}_B+q{a}^{\left(3\right)})}{q^3},\hfill \\ \hfill & \mathbb{E}\left[F^4\right]=\frac{4a(6a_B^2+6q{a}_B{a}_B^{\left(2\right)}+q^2{a}_B^{\left(3\right)})+q(12a^{\left(2\right)}{a}_B^2+4q{a}^{\left(3\right)}{a}_B+q(6a^{\left(2\right)}{a}_B^{\left(2\right)}+q{a}^{\left(4\right)}))}{q^4},\hfill \end{array}$$

where $a^{\left(k\right)}=\mathbb{E}\left[A^k\right]$, $a_B^{\left(k\right)}={\int }_0^{\infty }{x}^{k}B\left(x\right)dA\left(x\right)$. From the maximum entropy formalism, it is known that the maximum entropy solution is of the form

$$\begin{array}{c}\hfill \widehat{f}\left(t\right):=\widehat{f}(t;m)=e^{-c_0-{\sum }_{k=1}^m{c}_k{t}^k},t>0\end{array}$$

(32)

with $c_0=ln{\int }_0^1{e}^{-{\sum }_{k=1}^m{c}_k{t}^k}dt$. The Lagrangian coefficients $c_k$, $1\le k\le m$, can be calculated as a solution of the following minimization problem:

$$\begin{array}{c}\hfill g(c_1,\dots ,c_m)=ln{\int }_0^{\infty }{e}^{-{\sum }_{k=1}^m{c}_k(t^k-\mathbb{E}\left[F^k\right])}dt.\end{array}$$

To solve this problem and obtain optimal $c_k^{*}$, we use Nelder and Mead’s algorithm. This direct search method doesn’t use derivatives and hence can be used also in case of an almost singular Hessian of $g(c_1,\dots ,c_m)$.

In Figure 4, we plot the maximum entropy solution $\widehat{R}\left(t\right)$ for the function $R\left(t\right)$ in the case when RV A und B either simultaneously $A\sim \mathcal{LN}({\mu }_1,{\sigma }_1)$ and $B\sim \mathcal{LN}({\mu }_2,{\sigma }_2)$ or $A\sim \mathcal{PH}({\eta }_1,M_1)$ and $B\sim \mathcal{PH}({\eta }_2,M_2)$. We see that, in the log-normal case, the functions exhibit heavier tails by increasing of the mean value as. Figure 3a and Figure 4b can be used to compare the maximum entropy solution with a true function in the case of a PH-type distribution. We may notice no sufficient difference in graphics of functions.

Figure 5 and Figure 6 provide the TDSSPs for the PH-type distributed life- and repair-times for different values of a. It should be noticed that the state probability function ${\pi }_0\left(t\right)$ exhibits a concave structure with a corresponding minimum which is shifted to the left as a decreases. The function ${\pi }_0\left(t\right)$ in Figure 5a is shown only up to its minimum but further it increases as well and converges to a stationary value. In addition, we can see how fast the function ${\pi }_i\left(t\right)$ converges to a stationary value ${\pi }_i$ as a decreases. We strongly believe that the proposed above methodology of distribution recovery via moments will also be helpful for evaluation of the TDSSPs, when there are no closed forms of LSTs for life- and repair-time distributions.

The results of SSSP and mean lifetime calculation for different types of distribution pairs $A\left(t\right)$ and $B\left(t\right)$ are summarized in

Table 1. ${\pi }_i$, $i=0,1,2$, $\mathbb{E}\left[F\right]$ and $\mathbb{E}\left[F_1\right]$ versus $A\left(t\right)$ and $B\left(t\right)$.

$\mathit{A}\left(\mathit{t}\right)$, $\mathit{B}\left(\mathit{t}\right)$	${\mathit{\pi }}_0$	${\mathit{\pi }}_1$	${\mathit{\pi }}_2$	$\mathbb{E}\left[\mathit{F}\right]$	$\mathbb{E}\left[{\mathit{F}}_1\right]$
$\mathcal{E}$	0.9009	0.0900	0.0090	11.0000	12.0000
$\mathcal{PH}$	0.9009	0.0901	0.0090	11.0051	12.0051
$\mathcal{PAR}$	0.9001	0.0989	0.0010	519.0910	520.0910
$\mathcal{LN}$	0.9003	0.0966	0.0031	39.5976	40.5976

,

Table 2. ${\pi }_i$, $i=0,1,2$, $\mathbb{E}\left[F\right]$ and $\mathbb{E}\left[F_1\right]$, $A\sim \mathcal{E}\left(\alpha \right)$ versus $B\left(t\right)$.

$\mathit{B}\left(\mathit{t}\right)$	${\mathit{\pi }}_0$	${\mathit{\pi }}_1$	${\mathit{\pi }}_2$	$\mathbb{E}\left[\mathit{F}\right]$	$\mathbb{E}\left[{\mathit{F}}_1\right]$
$\mathcal{PH}$	0.9009	0.0901	0.0090	11.0000	12.0000
$\mathcal{PAR}$	0.9007	0.0919	0.0073	10.7972	11.7972
$\mathcal{LN}$	0.9009	0.0903	0.0088	10.9739	11.9739

,

Table 3. ${\pi }_i$, $i=0,1,2$, $\mathbb{E}\left[F\right]$ and $\mathbb{E}\left[F_1\right]$, $A\sim \mathcal{PH}({\eta }_1,M_1)$ versus $B\left(t\right)$.

$\mathit{B}\left(\mathit{t}\right)$	${\mathit{\pi }}_0$	${\mathit{\pi }}_1$	${\mathit{\pi }}_2$	$\mathbb{E}\left[\mathit{F}\right]$	$\mathbb{E}\left[{\mathit{F}}_1\right]$
$\mathcal{E}$	0.9009	0.0901	0.0090	11.0050	12.0050
$\mathcal{PAR}$	0.9007	0.0919	0.0073	11.0050	12.0050
$\mathcal{LN}$	0.9009	0.0903	0.0088	11.0050	12.0050

,

Table 4. ${\pi }_i$, $i=0,1,2$, $\mathbb{E}\left[F\right]$ and $\mathbb{E}\left[F_1\right]$, $A\sim \mathcal{PAR}(k,t_{\min })$ versus $B\left(t\right)$.

$\mathit{B}\left(\mathit{t}\right)$	${\mathit{\pi }}_0$	${\mathit{\pi }}_1$	${\mathit{\pi }}_2$	$\mathbb{E}\left[\mathit{F}\right]$	$\mathbb{E}\left[{\mathit{F}}_1\right]$
$\mathcal{E}$	0.9000	0.0993	0.0001	1297.1700	1298.1700
$\mathcal{PH}$	0.9000	0.0999	0.0001	1299.5000	1300.5000
$\mathcal{LN}$	0.9001	0.0994	0.0005	452.3640	453.3640

and

Table 5. Case 5: ${\pi }_i$, $i=0,1,2$, $\mathbb{E}\left[F\right]$ and $\mathbb{E}\left[F_1\right]$, $A\sim \mathcal{LN}(\mu ,\sigma )$ versus $B\left(t\right)$.

$\mathit{B}\left(\mathit{t}\right)$	${\mathit{\pi }}_0$	${\mathit{\pi }}_1$	${\mathit{\pi }}_2$	$\mathbb{E}\left[\mathit{F}\right]$	$\mathbb{E}\left[{\mathit{F}}_1\right]$
$\mathcal{E}$	0.9003	0.0969	0.0028	35.4011	36.4011
$\mathcal{PH}$	0.9003	0.0969	0.0028	35.3988	36.3988
$\mathcal{PAR}$	0.9002	0.0973	0.0024	59.4631	60.4631

. In all cases, the parameters were selected for mean and variance value constraints, i.e., $a=\mathbb{E}\left[A\right]=\frac{1}{\alpha }$, $b=\mathbb{E}\left[B\right]=\frac{1}{\beta }$, $\mathbb{V}\left[A\right]=\frac{1}{{\alpha }^2}$ and $\mathbb{V}\left[A\right]=\frac{1}{{\beta }^2}$. In experiments demonstrated in Table 1, the life- and repair time distributions have the same distributions but with different parameters fitting to $a=1$ and $b=0.1$. High mean values for the time to failure and the time between failures are a consequence of the Pareto distribution property that is defined for $t\ge t_{\min }$.

The exponential and PH-type distributions generate almost indistinguishable results. We also notice that the distributions with heavier tails decrease the probability of two operational units.

Table 2, Table 3, Table 4 and Table 5 present calculations for the systems with different pairs of life- and repair times distributions. The results for exponential and PH-type life type distributions in Table 2 and Table 3 are very similar. In addition, there are no sufficient differences in obtained results between different repair time distributions.

As it is shown in Table 4 and Table 5, the heavy-tailed Pareto and log-normal life-time distributions are a cause for reduction of the system failure probability and increase the mean time to failure and the mean time between failures. Hence, we notice that reliability and probabilistic characteristics of the system are more sensitive to changing the form of distribution functions when they exhibit heavier tails.

7. Conclusions

We provide probabilistic and reliability analysis of a cold double redundant renewable system with generally distributed life and repair times. We derive in closed-form the reliability function, the time-dependent as well as stationary system state distributions using LSTs and LTs. Numerical results illustrate sensitivity and effect of distribution functions and their parameters on reliability and probabilistic measures of the system. The presented system can be made preventive maintainable. We believe that the corresponding maintenance policy can be effective and leads to an increase in the system reliability.