*Article* **Reliability Simulation of Two Component Warm-Standby System with Repair, Switching, and Back-Switching Failures under Three Aging Assumptions**

**Kiril Tenekedjiev 1,2,\*, Simon Cooley 1, Boyan Mednikarov 2, Guixin Fan <sup>1</sup> and Natalia Nikolova 1,2**


**Abstract:** We analyze the influence of repair on a two-component warm-standby system with switching and back-switching failures. The repair of the primary component follows a minimal process, i.e., it experiences full aging during the repair. The backup component operates only while the primary component is being repaired, but it can also fail in standby, in which case there will be no repair for the backup component (as there is no indication of the failure). Four types of system failures are investigated: both components fail to operate in a different order or one of two types of switching failures occur. The reliability behavior of the system is investigated under three different aging assumptions for the backup component during warm-standby: full aging, no aging, and partial aging. Four failure and repair distributions determine the reliability behavior of the system. We analyzed two cases—in the First Case, we utilized constant failure rate distributions. In the Second Case, we applied the more realistic time-dependent failure rates. We used three methods to identify the reliability characteristics of the system: analytical, numerical, and simulational. The analytical approach is limited and only viable for constant failure rate distributions i.e., the First Case. The numerical method integrates simultaneous Algebraic Differential Equations. It produces a solution in the First Case under any type of aging, and in the Second Case but only under the assumption of full aging in warm-standby. On the other hand, the developed simulation algorithms produce solutions for any set of distributions (i.e., the First Case and the Second Case) under any of the three aging assumptions for the backup component in standby. The simulation solution is quantitively verified by comparison with the other two methods, and qualitatively verified by comparing the solutions under the three aging assumptions. It is numerically proven that the full aging and no aging solutions could serve as bounds of the partial aging case even when the precise mechanism of partial aging is unknown.

**Keywords:** state probability functions; partial aging in standby; Monte Carlo simulation; qualitative and quantitative verification of simulation model

### **1. Introduction**

Assessing the reliability of a system is a key engineering task that has economic and safety implications. Having a better understanding of failure/repair rates of system components is a tool to design highly reliable systems and conduct repair operations at adequate cost levels while complying with adequate and reasonable maintenance schedules. A common approach for improving the reliability is to provide redundancy for excessively failing components. The redundant components may operate simultaneously in a sense that the system will never fail if at least one of the parallel components operates [1]. Another possibility is to design a "*k* out of *n*" configuration where all *n* components are in operation and the system is not failing if at least *k* of them operate properly [2]. However, the standby

**Citation:** Tenekedjiev, K.; Cooley, S.; Mednikarov, B.; Fan, G.; Nikolova, N. Reliability Simulation of Two Component Warm-Standby System with Repair, Switching, and Back-Switching Failures under Three Aging Assumptions. *Mathematics* **2021**, *9*, 2547. https://doi.org/ 10.3390/math9202547

Academic Editors: Maria Luminit,a Scutaru and Catalin I. Pruncu

Received: 23 August 2021 Accepted: 6 October 2021 Published: 11 October 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

arrangement is the simplest, cheapest and the most utilized one; the system operates with some of its components (called primary) whereas the redundant components (called backup components) are in standby, but when a primary component fails, they take its place [3]. In this paper, we will focus on a two-component system with a standby arrangement where the backup components may fail either while in standby, or during operation after some imperfect switching mechanism has put those online. The switching mechanism can be continuous type when it actively monitors the primary component and makes its own decisions, but it can malfunction at any time [4]. However, in this paper, we will treat exclusively the widespread mechanisms that can fail only on demand when the switching is needed [5]. According to the failure intensities of the backup component, such systems are classified as cold-standby, warm-standby, and hot-standby [6]. In the hot-standby system, the intensity of failures of the backup component is the same during standby and operation, whereas in the cold-standby, there is no failure in standby. We will concentrate on the two-component warm-standby systems where there are failures of the backup component in standby, but with smaller intensity compared to its operational mode.

Additionally, the reliability of a system with backup components depends on the way of aging of the backup components while in standby. Previous works have identified three types of aging of the backup component during warm-standby: full aging, no aging and partial aging [7]. The full aging assumption means that the component changes its failure/repair rate during standby as if it is operational. Under the no aging assumption, the component does not change its failure/repair rate during standby. The partial aging assumption models the intermediate situation where the backup component experiences some wear during standby, but at a slower rate than if operational.

If some components of the system are deemed repairable, the system can be brought to its full operational capacity by replacing parts or by making adjustments [8]. In most works, the focus is on single-component repairable systems under various repair activities. A detailed discussion on how such tasks can be approached with modern statistical tools is offered in [9]. There are different types of repairs that can be adopted depending on objectives. The first possibility is the so-called perfect repair (a.k.a. as-good-as-new (AGAN) repair), where the primary component fails and it is replaced or restored to its original or good-as-new condition [10]. Minimal repair restores the device to the condition it was in immediately before the failure [11] (pp. 226–227). There may also be intermediate types of repairs (e.g., the partial perfect repair procedures mentioned in [8]). In the current work, the focus is on the case of minimal repair of the special type worse-than-old (WTO) [12]. The assumption is that during this repair, the non-repaired elements of the primary unit age as if the latter was operational.

In this paper, we will investigate the effects that adding repair and back-switching failures to a two-component warm-standby system with switching mechanism has on the reliability of the system. Our goal is to analyze how this affects the system reliability under different aging assumptions in standby. In such a system, the primary component begins operation, and when it fails, the system will try to activate the backup component, but a switching failure is possible. When the backup component is operational, the primary component undergoes minimal repair. If the latter finishes before the backup component fails, the system will try to activate back the primary component, but again a back-switching failure is possible. However, it is possible that the backup component will fail in standby. In that case, there will be no repair for the backup component since there will be no indication of the failure. The system is considered to have failed when either both components fail to operate at any given time, or when, after primary component failure, the switching mechanism fails on demand to switch the system operation to the secondary component (switching failure), or when, after a successful repair, the switching mechanism fails on demand to switch the system operation back to the primary component (back-switching failure). The primary component undergoes minimal repair, i.e., the primary component experiences full aging during the repair. The reliability behavior of the system will be investigated under three different aging assumptions for the backup component during

standby: full aging, no aging and partial aging. Only mechanical aging of the components will be considered, which excludes any influence by some software aging (for discussion of the latter topic see [13,14]).

The focus of our investigation would be a two-component warm-standby system with repair, switching, and back-switching failures, denoted as 2SBRSBF. Our study concentrates on the characteristics of the uptime of the 2SBRSBF. We only consider repairs of the failed primary component of a working system. The repair of the failed system, which relates to downtime, is an important component of system availability, but it is outside of the scope of our paper (for elaborate case study of data center availability using Markovian modeling, see [14]). In [15], the causes of system and component failures were classified as technological failures, natural disaster failures, and man-made disasters (e.g., terrorism). In our study, we will consider only the technological failure of 2SBRSBF since the other two types tend to cause dependent component failures, which is outside of the scope of our study. Here, the standby mode of the 2SBRSBF is defined as a situation, where the primary component is working properly, and the backup component is fully operational, but its failure will not affect the normal operation of the system at this moment (in [16] such component configuration is classified as "active/cold-standby").

The reliability behavior of the 2SBRSBF depends on four distributions: the failure and the repair distributions of the primary component, the failure distributions of the backup component in operation and in standby. We will analyze two cases for those distributions. In the First Case, all distributions will be with constant failure/repair rates. In the Second Case, the more realistic time-dependent failure/repair rates will be applied.

We will use three methods to identify the reliability characteristics of the 2SBRSBF: analytical, numerical, and simulational. The analytical approach is applicable for the First Case distributions. We will develop novel analytical solutions for the state probability functions in the case of exponential distributions. The numerical method creates and integrates simultaneous Ordinary Differential Equations (ODEs) for 2SBRSBF. This method is applicable for any set of First Case distributions and for Second Case distributions under the assumption of full aging in standby. However, there are no simultaneous ODEs that describe the behavior of 2SBRSBF with time-dependent distributions (i.e., Second Case) under no aging or partial aging assumptions in standby. To facilitate the simulational solution, we will introduce a novel method to generate failure times of the backup component in standby under the assumptions of full aging, no aging, or partial aging. Using this method, we will modify and generalize the algorithm from [17] to simulate the behavior of 2SBRSBF and to calculate its most important reliability characteristics. That algorithm will produce a novel simulation solution for any set of distributions (i.e., the First Case and the Second Case) under any of the three aging assumptions for the backup component in standby. The proposed algorithm will be validated quantitively by comparing with the analytical and with the numerical solutions (if those exist) as well as quantitatively by comparing with the full aging results.

In what follows, Section 2 summarizes the state-of-the-art in the field and outlines the contributions of our paper. Section 3 will setup the problem for reliability characteristics assessment of a 2SBRSBF function. In Section 4, we present a novel analytical solution of the formulated problem in the case of distributions with constant failure/repair rate. A numerical solution will be identified in Section 5 where a system of four simultaneous deferential algebraic equations will model the 2SBRSBF in the case of full aging of the backup component during standby. In Section 6, the same problem will be solved with simulation which can be used with any distributions under three different assumptions about the aging mechanism of the backup component. Section 7 contains the results of three numerical examples, where we validate the proposed simulational algorithm quantitatively (by comparing with the analytical and the numerical solutions when those exist) and qualitatively (by checking whether the effects of no aging and partial aging correspond to the logically expected ones). Section 8 concludes the paper.

#### **2. Related Works and Contributions of the Paper**

Although the publications about warm-backup system reliability are growing recently, they are rare in comparison with reliability studies of cold-backup and hot-backup system, since the realistic models of the former tend to be more elaborate [18]. In [19] (pp. 113–115), the analytical solution for two-component warm-standby system with switching failure (2SBSF) was developed. The switching mechanism fails on demand. The failure distributions were considered exponential. Hence, no aging effect was taken into account. Explicit formulae were derived for the reliability of the system and for all state probability functions. In [20] (pp. 167–170), a model of a two-component warm-standby system (2SB) with arbitrary failure distributions was proposed. Although no particular simulation algorithm was developed, general advice was given on how to acquire the reliability function and the state probability functions using Monte Carlo simulation and how to deal with different aging assumptions. In [21], a model of a 2SB system was proposed under general standby, which generalizes the three special cases of warm-, cold- and hot-standby. The failure distributions can be arbitrary. The aging effects are accounted for using a pre-specified virtual aging function. An integral equation, connecting the failure rates and the virtual aging function with the reliability of the system was proposed. In [6], these results from [21] were generalized to solve the problem of allocation of redundancy that includes two independent and one generalized standby component. The reliability and the state probability functions of a generic two-component standby system under full aging, no aging, and partial aging were identified with a simulation algorithm in [22] using arbitrary failure distributions. That solution is verified with analytical and numerical special cases. The results from this work were expanded in [17] to model the 2SBSF, but some numerical problems connected with random variate generation and arbitrary failure rate calculations were resolved.

The majority of the above models consider aging effects, but none of them has repairs.

A 14-states model of two dissimilar warm-standby subsystems in series with repair were discussed in [23]. The failure distributions are exponential, and the system is with constant repair rates. The type of repair is AGAN. The problem of aging is not considered. Some analytical steady-state characteristics of the system are provided using Laplace transforms. Those characteristics for two-component warm-standby system with repair (2SBR) can be obtained as special cases from the results in the paper. The work [4] performs reliability analysis for a two-component warm-standby system with repair and switching failures (2SBRSF). The failure and repair rates are constant. The switching mechanism is of continuous type and has its own failure distribution. This leads to a possibility of repairing the failed backup unit while the primary component operates. All failure and repair distributions are exponential. Any failure of the switch leads to system failure. The repairs are AGAN and no aging is considered. The system has 10 states. The reliability and the state probability functions of the system were identified with a numerical algorithm as a solution of an ODE system. Another interesting two-identical-component standby system is given in [24]. The type of standby is difficult to determine since the failure in standby mode is deterministic and happens after surpassing a pre-specified time. The failure distribution of the operating unit is exponential, but the repair rate is arbitrary. There is no switching failure, but the switching mechanism inspects the failed standby unit and decides whether to replace it or to repair it. No aging is considered in this model. Some steady-state measures of reliability are obtained using semi-Markov models. In [25], the authors propose a system with *m* identical components working in parallel with *s* components in warm-standby. The system includes a service station that can also fail and be repaired. There are no switching failures, and all failure and repair distributions of the components are exponential. The failure and repair distributions of the service station are also exponential. The repairs are AGAN, and no aging is considered. The reliability and the state probability functions of the system were approximated using symbolic computer software. The work [18] presents a system of *n* components in series with one component in warm-standby. There are neither switching failures nor aging

considerations. The failure distributions are exponential, but the repair distributions are arbitrary. The system is also subjected to non-repairable failures. Some reliability and availability steady-state characteristics of the system are derived using Laplace transforms. In [26], the authors discuss a three identical component warm-standby system. Initially, the primary component is working, and the other two are in standby. The failure of the operating unit and the repairs are with random distributions, however in standby there is constant failure rate. The repairs are AGAN, there are no switching failures, and no aging is considered. An integral equation, connecting the failure rates with the reliability of the system is proposed.

The models with repair discussed above do not consider any aging effects.

In Table 1, we summarized seven characteristics for each of the above-discussed 12 papers plus the current work. The information in Table 1 highlights the novelty of our work against the discussed state-of-the-art studies in the literature. The contributions of our study can be outlined as follows:



**Table 1.** Overview of the state-of-the-art publications in the warm-standby area.

#### **3. States, Transition Rates, and Distributions**

The dynamics of a 2SBRSBF system can be determined by its transition between several possible states [27]. The 2SBRSBF has four major states, but State 4 (where the 2SBRSBF system is not operational) is subdivided into 4 substates, called types.

In State 1, the primary component operates, the backup component is fully operational but is in standby. Sooner or later, one of the two components will fail:


In State 2 sooner or later either the primary component will be repaired, or the backup component will fail. Then one of the following two events will occur:


In State 3, sooner or later, the primary component will fail and there will be no operational backup component to take over. The system will transit to State 4 where the 2SBRSBF system is not operational (type *d* system failure).

The State 4, where the 2SBRSBF system is not operational, is irreversible in our model regardless of the type of the system failure.

The described system is partially observable since we will not know whether the system is in State 1 or in State 3, but State 4 and State 2 are observable. At the same time, 2SBRSBF is controllable by three trivial event-driven decisions: (a) when the primary component fails, attempt to move to State 2, by switching to the backup unit; (b) when the backup unit is in operation, start repairing the primary component; (c) when the primary component is repaired, attempt to move to State 1, by back-switching to the primary unit.

The state function *Pg*(*t*) (for *g* = 1,2,3,4) measures the probability of the 2SBRSBF to be in State *g* at time *t* (for *t* ≥ 0). Since the system will be in one state and in one state only at any non-negative time moment *t*, then:

$$P\_1(t) + P\_2(t) + P\_3(t) + P\_4(t) = 1,\text{ for } t \in [0, \infty) \tag{1}$$

The 2SBRSBF system starts in fully operational mode so initially it will be in State 1:

$$P\_1(0) = 1 \text{ and } P\_2(0) = P\_3(0) = P\_4(0) = 0 \tag{2}$$

If the four state functions are identified, then the 2SBRSBF system is quantitatively described and we can calculate all its reliability characteristics. The reliability of the system is the sum of the first three state probabilities (i.e., the probability not to be in State 4):

$$R\_{sys}(t) = P\_1(t) + P\_2(t) + P\_3(t) = 1 - P\_4(t), \text{ for } t \in [0, \infty) \tag{3}$$

The mean time to failure (MTTF) of the 2SBRSBF system can be calculated as:

$$MTTF\_{\text{sys}} = \int\_0^\infty R\_{\text{sys}}(t)dt\tag{4}$$

The time for which the reliability of the system will be *α* is known as *α*-design life (*tdes,α*). It can be identified as the unique solution of Equation (5) in the domain *tdes*,*<sup>α</sup>* ∈ (0, ∞):

$$R\_{\text{sys}}(t\_{\text{des},\mathfrak{a}}) = \mathfrak{a}, \text{ for } \mathfrak{a} \in (0,1) \tag{5}$$

The median (*Mediansys*), the B1 life (*B*1*\_life*), the B10-life (*B*10*\_life*), and the interquartile range (*IQRsys*) of the 2SBRSBF system reliability can be easily estimated using Equation (5) respectively as *tdes*,0.5, *tdes*,0.99, *tdes*,0.9, and *tdes*,0.25 − *tdes*,0.75 [20] (pp. 87–88).

To identify the four required state functions of the 2SBRSBF system, we need to know:


Each of the four PDFs, *fk*(*t*), (for *k* = 1, 2, 3, 4) can be transformed into four alternative forms: a cumulative distribution function (CDF), *Fk*(*t*), a failure/repair rate, *λk*(*t*) (as shown in [28]), a complementary CDF, or *Rk*(*t*), and an inverse CDF, i.e., *F*−<sup>1</sup> *<sup>k</sup>* (*p*). The five forms *fk*(*t*), *Fk*(*t*), *λk*(*t*), *Rk*(*t*), and *F*−<sup>1</sup> *<sup>k</sup>* (*p*) contain the same information and are equivalent. In the ideal world the domain of the first four functions and the range of the last one will be *t* ∈ [0, ∞) where *t* can be interpreted as time. However, this is not always the case—those failure and repair distributions are based on information about the behavior of the components. The first step is to summarize the available information in several nodes of the CDF. If the reliability information is in the form of fully observed or multiply sensor data, then we can produce an empirical distribution, using either the Kaplan-Meier product limit estimator method [29] (see the function *ecdf.m* in [30], which embodies the method) or the invertible ECDF estimator with maximum count of nodes [31], or any other modern method. If the information is in the form of expert knowledge, then we can extract subjective quantiles using the triple bisection method [32] as described in [33]. The second step is to fit a parametric distribution of some type to the nodes of the CDF identified in the first step. The work [20] (p. 399) gives several reasons to use parametric distributions rather than empirical ones, with the most important one being that empirical distributions can only be trusted at the beginning of the failure/repair process. Regardless of the method utilized to identify the parameters in the second step (least square, maximum likelihood estimation, Bayesian estimation, etc.), it is quite possible that some of the derived parametric distributions would have substantial support for negative values of the argument *t*. For purely pragmatic reasons, we assume that for each *k*, we are given only procedures to calculate *fk*(*t*), *Fk*(*t*), and *F*−<sup>1</sup> *<sup>k</sup>* (*p*). Such numerical procedures exist in almost any software package. For example, the Statistics and Machine Learning Toolbox in MATLAB contains the *pdf.m*, *cdf.m*, and *icdf.m* which calculate the PDF, the CDF, and the inverse CDF values for any distribution object created by the *makedist.m* [30]. The latter can choose a wide variety of parametrical 1D distributions with arbitrary specified parameters. Unluckily, some of those parametrical distributions are defined over the whole real axis (e.g., the normal distribution, or the extreme value distribution). Traditionally, no numerical procedures are given for estimating the values of *λk*(*t*) and *Rk*(*t*), which have to be approximated using *fk*(*t*), *Fk*(*t*), *F*−<sup>1</sup> *<sup>k</sup>* (*p*). In this paper, any of the procedures *fk*(*t*), *Fk*(*t*), *Rk*(*t*), *λk*(*t*), *F*−<sup>1</sup> *<sup>k</sup>* (*p*) will be called the *k*th original distribution since the five of them describe in alternative form the uncertainty of a real continuous variable:

$$\begin{cases} \text{ (a) } f\_k(t), \text{ for } k = 1, 2, 3, 4 \text{ with Domain } t \in ( -\infty, +\infty) \\ \qquad \text{ (b) } F\_k(t) = \int\_0^t f\_k(t)dt, \text{ for } k = 1, 2, 3, 4 \text{ with Domain } t \in ( -\infty, +\infty) \\ \qquad \quad \text{ (c) } \lambda\_k(t) = f\_k(t) / [1 - F\_k(t)], \text{ for } k = 1, 2, 3, 4 \text{ with Domain } t \in ( -\infty, +\infty) \\ \qquad \text{ (d) } R\_k(t) = 1 - F\_k(t), \text{ for } k = 1, 2, 3, 4 \text{ with Domain } t \in ( -\infty, +\infty) \\ \qquad \text{ (e) } F\_k^{-1}(p), \text{ for } k = 1, 2, 3, 4 \text{ with Domain } p \in [0, 1] \end{cases} \tag{6}$$

Here, *Rk,*(*t*) from Equation (6) a is aka original reliability/repair function when the real argument *t* is non-negative and can be interpreted as time. In our problem, the argument *t* would be most often the time (or other suitable non-negative variable, e.g., mileage), so we will use the original distribution in Equation (6) a–e to approximate their truncated versions which take the form of conditional distributions provided that the failure/repair has not happened till time 0:

$$\begin{cases} \text{ (a)} \, f\_{k,trm}(t) = f\_k(t|0) = f\_k(t)/R\_k(0), \text{for} = 1, 2, 3, 4 \text{with } \text{Domain} \in [0, +\infty) \\ \text{ (b)} \, F\_{k,trm}(t) = F\_k(t|0) = 1 - R\_k(t)/R\_k(0), \text{for} = 1, 2, 3, 4 \text{with } \text{Domain} \in [0, +\infty) \\ \text{ (c)} \, \lambda\_{k,trm}(t) = \lambda\_k(t|0) = f\_{k,trm}(t)/[1 - F\_{k,trm}(t)], \text{for} = 1, 2, 3, 4 \text{with } \text{Domain} \in [0, +\infty) \\ \text{ (d)} \, R\_{k,trm}(t) = R\_k(t|0) = 1 - F\_{k,trm}(t), \text{for} = 1, 2, 3, 4 \text{with } \text{Domain} \in [0, +\infty) \\ \text{ (e)} \, F\_{k,trm}^{-1}(p) = F\_k^{-1}(p|0), \text{for} = 1, 2, 3, 4 \text{with } \text{Domain} \in [0, 1] \end{cases} \tag{7}$$

In this paper, any of the functions *fk,trun*(*t*), *Fk,trun*(*t*), *Rk,trun*(*t*), *λk,trun*(*t*), *F*−<sup>1</sup> *<sup>k</sup>*,*trun*(*p*) will be called the *k*th truncated distribution, since the five of them describe in alternative forms the uncertainty of a real non-negative continuous variable which can be interpreted as time. The *Rk,trun*(*t*) from Equation (7) d is aka truncated reliability/repair function. Let us concentrate on the 2SBRSBF system at time *t*:


• The rate for transitioning between State 2 and State 4 (type *c* system failure) will depend on *P*2(*t*) and on the conditional failure distribution *f* 2(*τ|tage*) (failure density of the backup component in operation, given that it has not failed till time *tage*). Here *tage* is the equivalent aging of the backup component in operation. It depends on the type of aging and possibly on the backup component history of utilization (alternating between operational and standby modes).

The four state functions of 2SBRSBF system can be identified using computer simulation in the above setup for any set of distributions and aging assumptions during standby. However, for verification purposes, two alternative solution methods can be developed for some special cases of the 2SBRSBF system. This approach was successfully applied in [34] for verification of a novel simulation-based optimization algorithm used in redundancy allocati on problems using Markovian models as special cases.

If we have a set of First Case distributions, then all state transitions will depend on the absolute densities, rather than from conditional ones. The reason is that the exponential distributions have no memory, and hence any aging assumptions are irrelevant. Then the probabilities for transitioning between the states depend only on the current state of the system, but not on the history describing how the system turns out to be in the current state. This means that the 2SBRSBF system with First Case distributions degenerates to a Markov model [27] (more precisely to a partially observable Markov decision process [35]). Such Markov model can be conveniently visualized with the Rate Diagram (RD) [20] (pp. 155–170) shown in Figure 1a. Using that RD, we will derive an analytical solution for the four state probability functions of the 2SBRSBF system with First Case distributions.

**Figure 1.** Rate diagram for 2SBRSBF with: (**a**) First Case distributions; (**b**) Second Case distributions with full aging in standby.

If we have a system with full aging assumption during standby, then the equivalent aging of the backup component in operation rate *tage*, described above in the transitioning between State 2 and State 4 (type *c* system failure) will be simply the current time *t*. The reason is that the backup component is assumed to age during standby in the same fashion as in operation, which shows that any failure of the backup component during operation will behave like a first failure at time *t*. This means that the 2SBRSBF system with full aging assumption degenerates to a semi-Markov model where the transition probabilities depend not only on the current state but also on the current time [36]. The semi-Markov model can be conveniently visualized with the Generalized Rate Diagram (GRD) shown in Figure 1b [37] (pp. 521–526). Using the GRD, we can describe the 2SBRSBF system with simultaneous ODEs. This is possible because the failure/repair rate of any distribution, *F*(*t*), at time *t\** coincides with the failure/repair rate of the conditional distribution *F*(*τ|T*) at the same time *t\**, if *T* ≤ *t\**. This trivial fact is proven in Appendix A. The derived Cauchi problem can be solved numerically. Obviously, such solution exists also for the First Case distribution, which will allow the comparison of the analytical and the numerical solutions.

Neither the analytical, nor the numerical solutions can be derived for the cases of the Second Case distribution under the assumptions of no aging and partial aging since general aging effects cannot be described by any Markovian or semi-Markovian model and there is no system of ODE which fully quantifies the reliability behavior of 2SBRSBF unless when the primary component is subjected to full aging in standby (see [14,36]).

#### **4. Analytical Solution**

This solution is applicable only for First Case distributions, where the failure/repair rates are constant. The rate diagram in Figure 1a can be represented as a system of three ODEs from Equation (8) about the first three state probability functions [38]:

$$\begin{cases} \frac{dP\_1}{dt}(t) = -(\lambda\_1 + \lambda\_3)P\_1(t) + (1 - p\_r)\lambda\_4 P\_2(t) \\ \frac{dP\_2}{dt}(t) = \left(1 - p\_f\right)\lambda\_1 P\_1(t) - (\lambda\_4 + \lambda\_2)P\_2(t) \\ \frac{dP\_3}{dt}(t) = \lambda\_3 P\_1(t) - \lambda\_1 P\_3(t) \end{cases} \tag{8}$$

The initial conditions are given in Equation (2). After solving Equation (8), the last probability function, *P*4(*t*), can be estimated from Equation (3) as the complement to 1 of the other state probability functions. The analytical solution of 2SBRSBF with First Case distributions can be described as: "set the constants from Equation (9) and form the state probability functions from Equation (10)" (see Appendix B for the proof).

$$\begin{array}{l} K = \left(\lambda\_1 + \lambda\_2 + \lambda\_3 + \lambda\_4\right) / 2; \mathbf{C} = \left(\lambda\_1 + \lambda\_3\right) \left(\lambda\_2 + \lambda\_4\right) - \left(1 - p\_f\right) \left(1 - p\_r\right) \lambda\_1 \lambda\_4\\ s\_1 = -K + \sqrt{K^2 - \overline{C}}; s\_2 = -K - \sqrt{K^2 - \overline{C}}\\ A\_1 = \frac{s\_1 + \lambda\_2 + \lambda\_4}{s\_1 - s\_2}; B\_1 = \frac{s\_2 + \lambda\_2 + \lambda\_4}{s\_2 - s\_1}; A\_2 = \frac{\left(1 - p\_f\right)\lambda\_1}{s\_1 - s\_2}; B\_2 = \frac{\left(1 - p\_f\right)\lambda\_1}{s\_2 - s\_1}\\ A\_3 = \frac{\lambda\_3(s\_1 + \lambda\_2 + \lambda\_4)}{(s\_1 - s\_2)(\lambda\_1 + s\_1)}; B\_3 = \frac{\lambda\_3(s\_2 + \lambda\_2 + \lambda\_4)}{(s\_2 - s\_1)(\lambda\_1 + s\_2)}; C\_3 = \frac{\lambda\_3(\lambda\_1 + \lambda\_2 + \lambda\_4)}{(\lambda\_1 - s\_1)(\lambda\_1 - s\_2)} \end{array} \tag{9}$$

$$\begin{aligned} \text{Domain}: t \in [0, \infty) \\ \begin{cases} P\_1(t) = A\_1 e^{s\_1 t} - B\_1 e^{s\_1 t} \\ P\_2(t) = A\_2 e^{s\_1 t} - B\_2 e^{s\_2 t} \\ P\_3(t) = A\_3 e^{s\_1 t} - B\_3 e^{s\_2 t} + C\_3 e^{-\lambda\_1 t} \\ P\_4(t) = 1 - (A\_1 + A\_2 + A\_3) e^{s\_1 t} + (B\_1 + B\_2 + B\_3) e^{s\_2 t} - C\_3 e^{-\lambda\_1 t} \end{cases} \end{aligned} \tag{10}$$

The reliability of the system from Equation (11) and its MTTF from Equation (12) are derived as special cases of Equations (3) and (4):

$$R\_{sys}(t) = (A\_1 + A\_2 + A\_3)e^{s\_1t} - (B\_1 + B\_2 + B\_3)e^{s\_2t} + C\_3e^{-\lambda\_1 t}, \text{ for } t \in [0, \infty) \tag{11}$$

$$MTTF\_{sys} = -\left(A\_1 + A\_2 + A\_3\right)/s\_1 + \left(B\_1 + B\_2 + B\_3\right)/s\_2 + C\_3/\lambda\_1\tag{12}$$

#### **5. Numerical Solution**

This solution is applicable for any Second Case distribution with full aging of the backup component in standby and for any First Case distribution. The GRD in Figure 1b can be represented as a system of four simultaneous DAEs from Equation (16) about the four state probability functions, *Pg*(*t*) (*g* = 1,2,3,4). The system of DAEs will be numerically integrated from 0 to *tend*, where the latter will be selected sufficiently large, so *Rsys*(*tend*) ≈ 0 (< 0 .01). The main numerical difficulty in solving Equation (16) is to advise a procedure for stable approximation of the failure/repair rates, *λk*(*t*) (*k* = 1,2,3,4), at any *t* ∈ [0, *tend*]. That problem is far from trivial since sometimes *Fk*(*t*) is so close to 1, that the denominator of Equation (7) turns into 0. For each of the four distributions, using the

original inverse CDF function, we can calculate the time *tλ*,*k*, where the denominator of Equation (7) equals to 100 times the machine epsilon ():

$$t\_{\lambda,k} = F\_k^{-1}(1 - 100\varepsilon), \text{ for } k = 1, 2, 3, 4 \tag{13}$$

The approximated failure/repair rate, *λk,a*(*t*) (*k* = 1,2,3,4) equals to Equation (7) if its denominator is greater than 100 or equals the failure/repair rate at *tλ*,*<sup>k</sup>* otherwise:

$$\lambda\_{k,\mathfrak{a}} = \begin{cases} f\_k(t) / \left[1 - F\_k(t)\right] & , \, t \in \left[0, t\_{\lambda,k}\right] \\ f\_k(t\_{\lambda,k}) / \left[1 - F\_k(t\_{\lambda,k})\right] & , \, t \in \left(t\_{\lambda,k}, \infty\right) \end{cases} \text{ where } k = 1, 2, 3, 4 \tag{14}$$

Equation (14) produces numerically stable approximations of the failure/repair rates at any non-negative time not greater than *tend*. This is true even when a distribution is truncated which means that *Fk*(0) > 0 and only its part in the non-negative domain has to be used. Then, according to Appendix A, the value of the failure/repair rate for any non-negative time will be the same as that of the non-truncated distribution since the truncated distribution can be represented as a conditional nontruncated one:

$$F\_{k,trun}(t) = F\_k(t|T\_0 = 0) = 1 - \left[1 - F\_k(t)\right] / \left[1 - F\_k(0)\right], \ t \ge 0\tag{15}$$

Now, we can write the DAE system corresponding to Figure 1b:

$$\begin{cases} \frac{dP\_1}{dt}(t) = -[\lambda\_{1,a}(t) + \lambda\_{3,a}(t)]P\_1(t) + (1 - p\_r)\lambda\_{4,a}(t)P\_2(t) \\ \frac{dP\_2}{dt}(t) = \left(1 - p\_f\right)\lambda\_{1,a}(t)P\_1(t) - [\lambda\_{4,a}(t) + \lambda\_{2,a}(t)]P\_2(t) \\ \frac{dP\_3}{dt}(t) = \lambda\_{3,a}(t)P\_1(t) - \lambda\_{1,a}(t)P\_3(t) \\ 0 = P\_1(t) + P\_2(t) + P\_2(t) + P\_2(t) - 1 \end{cases} \tag{16}$$

The dependent variables can be organized in a 4D vector: <sup>→</sup> *y* (*t*) = [*P*1(*t*), *P*2(*t*), *P*3(*t*), *P*4(*t*)]*T*. The DAE from Equation (16) is semi-explicit with differential index 1. It has a singular constant mass matrix:

$$\mathcal{M}\begin{pmatrix}t,\vec{\boldsymbol{y}}\end{pmatrix} = \begin{pmatrix}1 & 0 & 0 & 0\\0 & 1 & 0 & 0\\0 & 0 & 1 & 0\\0 & 0 & 0 & 0\end{pmatrix} \tag{17}$$

The Jacobian matrix of the RHS of Equation (16) depends only on the time *t*:

$$J\begin{pmatrix} -\lambda\_{1,a}(t) - \lambda\_{3,a}(t) & (1 - p\_r)\lambda\_{4,a}(t) & 0 & 0\\ (1 - p\_r)\lambda\_{1,a}(t) & -\lambda\_{4,a}(t) - \lambda\_{2,a}(t) & 0 & 0\\ \lambda\_{3,a}(t) & 0 & -\lambda\_{1,a}(t) & 0\\ 1 & 1 & 1 & 1 \end{pmatrix} \tag{18}$$

The initial conditions given in Equation (2) together with Equations (16) and (17) form a Cauchi problem:

$$\mathbf{M}\begin{pmatrix}t\ \stackrel{\rightarrow}{y}\end{pmatrix}\stackrel{\rightarrow}{y}^{t}(t) = \stackrel{\rightarrow}{f}\begin{pmatrix}t, \stackrel{\rightarrow}{y}(t)\end{pmatrix}\text{ with }\stackrel{\rightarrow}{y}\_{\text{ini}} = \stackrel{\rightarrow}{y}(0) = \left[P\_1(0), P\_2(0), P\_3(0), P\_4(0)\right]^T = \left[1, 0, 0, 0\right]^T\tag{19}$$

Here, <sup>→</sup> *y* (*t*) = [*dP*1(*t*)/*dt*, *dP*2(*t*)/*dt*, *dP*3(*t*)/*dt*, *dP*4(*t*)/*dt*] *<sup>T</sup>*, *M t*, → *y* is the mass ma-

trix (17), and the 4D <sup>→</sup> *f t*, → *y* (*t*) is the RHS of Equation (16). The problem from Equation (19) can be numerically integrated (e.g., using *ode15s.m* from MATLAB [39]) at 2000 evenly distributed time points from 0 to *tend*:

$$t\_i = (i - 1)t\_{end} / 1999\text{ , for } i = 1, 2, \dots, 2000\tag{20}$$

The reliability function and the*MTTFsys* can be calculated approximatingEquations (3) and (4)

$$R\_{sys}(t\_i) = 1 - P\_4(t\_i) \text{, for } i = 1, 2, \dots, 2000 \tag{21}$$

$$MTTF\_{sys} = \left[ R\_{sys}(t\_1) + R\_{sys}(t\_{2000}) + 2\sum\_{i=2}^{1999} R\_{sys}(t\_i) \right] t\_{end} / 1999 \tag{22}$$

#### **6. Simulation Solution**

as:

This solution is applicable for any set of distribution (First Case or Second Case) and for any type of aging of the backup component in standby (full aging, no aging, or partial aging). Any simulation uses multiple pseudo-realities to study the system in question. The information from each generated pseudo-reality will be kept in an EC, whose definition and properties will be discussed in Section 6.1. In Section 6.2 we will concentrate on the development of specific functions generating random time intervals for the 2SBRSBF system. Those functions will be used in Section 6.3 where an algorithm will be developed to generate a random EC describing the 2SBRSBF system. In Section 6.4 we will extract the information in the generated ECs to calculate the state probability functions and the rest of the reliability characteristics of a 2SBRSBF system.

#### *6.1. Definition and Properties of the Event Chains for 2SBRSBF*

In the simulational solution, we generate a large count *N* of pseudo-realities in which we observe the behavior of the 2SBRSBF system from time 0 to system failure or to time *tend* whichever comes first. As in the numerical solution (described in Section 5) the constant *tend* is selected sufficiently large, so *Rsys*(*tend*)< 0.01. The pseudo-realities are described with the ECs introduced in [22] where the EC of the *j*th pseudo-reality is defined as the set:

$$EC\_{\hat{\jmath}} = \left\{ \left[ timepsr\_{\hat{\jmath}}(k), statepsr\_{\hat{\jmath}}(k) \right] \, | \, k = 1, 2, \dots, q\_{\hat{\jmath}} \right\} \tag{23}$$

The notation in Equation (23) shows that the *j*th pseudo-reality contains *qj* state transitions (called events) where the *k*th consecutive event which happened at time *timepsrj*(*k*) is a transition to state/substate *timepsrj*(*k*). The latter is coded either with 1, 2, and 3 respectively for State 1, State 2, and State 3, or with 40, 41, 42, and 43 respectively for system failure type *b*, type *a*, type *c*, and type *d* (all of them denoting State 4). Any EC for a 2SBRSBF system should have the following properties:

**p1)** It contains at least one event: *qj* ≥ 1.

**p2)** The events happen at strictly increasing times: *timepsrj*(*k*)< *timepsrj*(*k +* 1) for *k* = 1,2, . . . ,(*qj* − 1).

**p3)** The initial event is at time zero: *timepsrj*(1) = 0.

**p4)** The final event happens before *tend*: *timepsrj*(*qj*) < *tend*.

**p5)** The simulation starts with fully operational system: *statepsrj*(1) = 1.

**p6)** Whenever a system failure is observed the simulation ends: if *statepsrj*(*b*) > 3, then *qj* = *b*.

**p7)** Whenever the State 3 is observed either it is the last event, or the next event is the system failure type *d*: if *statepsrj*(*b*) = 3, then either *qj* = *b*, or *qj* = (*b +* 1) and *statepsrj*(*qj*) = 43. **p8)** The State 3 and the State 4 (in all its substates) can happen only once: #[*statepsrj*(*k*) = 3] ≤ 1,

#[*statepsrj*(*k*) = 40] ≤ 1, #[*statepsrj*(*k*) = 41] ≤ 1, #[*statepsrj*(*k*) = 42] ≤ 1, #[*statepsrj*(*k*) = 43] ≤ 1.

**p9)** The State 1 and State 2 alternate in the beginning of the EC including to the *h*th event and neither one happens later: *statepsrj*(*k*) = 1 if and only if *k* is odd and *k* ≤ *h*, whereas *statepsrj*(*k*) = 2 if and only if *k* is even and *k* ≤ *h.*

**p10)** There could be maximum two events after *h*: *h* ≤ *qj* ≤ (*h* + 2).

**p11)** If there are events after the *h*th one, they are either a transition to State 3 or a transition to State 4 (in all its substates): *statepsrj*(*k*) ≥ 3 for all *k* > *h* and *k* ≤ *qj.*

**p12)** The State 2 can be observed only on an even position and the previous event is always a transition to State 1: if *statepsrj*(*b*) = 2, then *b* is even and *statepsrj*(*b* − 1) = 1.

**p13)** The State 3 can be observed only on an even position and the previous event is always a transition to State 1: if *statepsrj*(*b*) = 3, then *b* is even and *statepsrj*(*b* − 1) = 1.

The formulated EC properties will facilitate the generation of time-period variates presented in Section 6.2. The algorithm described in Section 6.3 will generate ECs with the formulated EC properties. The latter will be used in Section 6.4 to prove the methods for extracting reliability information from the generated set of ECs for 2SBRSBF system.

#### *6.2. Generating Times Periods Using Conditional Distributions from 2SBRSBF*

As discussed in Section 3, to simulate an EC of a 2SBRSBF system we need to generate random time-periods complying with the conditional failure distributions *f* 1(*τ|t*), *f* 3(*τ|t*), and *f* 2(*τ|tage*) and with the conditional repair distribution *f* 4(*τ|t*), where *t* and *tage* are non-negative values.

We do not know which of the four original distributions, *fk*(*t*) (*k* = 1,2,3,4), are defined only in the non-negative domain and which are defined in the entire real axes so we need to substitute them with their truncated distributions, *fk*,*trunc*(*t*) = *fk*(*t|*0) for *k* = 1,2,3,4. Noting that if the first condition is met, then *fk*,*trunc*(*t*) = *fk*(*t|*0) = *fk*(*t*) (*k* = 1,2,3,4), and we can safely work only with truncated distributions. So, strictly speaking, we need to generate time-period variates from the conditional truncated distributions *f* 1,*trun*(*τ|t*), *f* 3,*trun*(*τ|t*), *f* 2,*trun*(*τ|tage*), and *f* 4,*trun* (*τ|t*). However, for any *k* it is true that:

$$f\_{k,t\text{run}}(\mathbf{r}|t) = \frac{f\_{k,t\text{run}}(\mathbf{r}+t)}{R\_{k,t\text{run}}(t)} = \frac{f\_k(\mathbf{r}+t|0)}{R\_k(t|0)} = \frac{f\_k(\mathbf{r}+t+0)/R\_k(0)}{R\_k(t+0)/R\_k(0)} = \frac{f\_k(\mathbf{r}+t)}{R\_k(t)} = f\_k(\mathbf{r}|t) \tag{24}$$

According to Equation (24) the conditional truncated distributions coincide with the conditional original distributions. In case *t* and *tage* are known entities we can generate random time-period variates as special cases of the Practical Indirect Sampling Method from Conditional CDF (PISMCF) [17] where the algorithm is motivated, formalized, illustrated, and proven. On its basis we can define a three-attribute procedure, PISMCF(.), which generates numerically stable random time interval variate, Δ*τ*, from a given conditional CDF, *F*(*t*|*Tsurv*), where *Tsurv* is a non-negative real number representing the time of survival:

$$
\Delta \tau = \text{PISMCF}\left(F(.), F^{-1}(.), T\_{surv}\right) \tag{25}
$$

In Equation (25), *F*(.) is the unconditional CDF which can express *F*(*t*|*Tsurv*) using Equation (26):

$$1 - F(t|T\_{surv}) = \frac{1 - F(t + T\_{surv})}{1 - F(T\_{surv})} \tag{26}$$

The second argument, *F*–1(.), of the PISMCF procedure from Equation (25) being the inverse CDF, can be used to estimate the time *t<sup>λ</sup>* where the denominator of Equation (26) is 100 machine epsilons ():

$$t\_{\lambda} = F^{-1}(1 - 100\varepsilon) \tag{27}$$

In short, the algorithm for estimating Equation (25) is: (a) Calculate *t<sup>λ</sup>* using Equation (27); (b) if *Tsurv* < *tλ*, then set *Tcut* = *Tsurv*, else set *Tcut* = *tλ*; (c) Generate *RD* as a uniformly distributed variate in the unit interval (0,1); (d) estimate *pRD* = 1 − *RD* [1 − *F*(*Tcut*)]; (e) Set Δ*τ* = *F*−1(*pRD*).

Let us assume that while performing the simulation of the *j*th pseudo-reality for the 2SBRSBF system we have observed only the first *kcur* state events. The simulation probably will continue and therefore, the EC is yet incomplete:

$$EC\_{j}^{inc} = \left\{ \left[ timepsr\_{j}(k), statepsr\_{j}(k) \right] \middle| k = 1, 2, \cdots, k\_{cur} \right\} \tag{28}$$

Then, the current state of the system is *scur* = *statesprj*(*kcur*) and the simulational time is *Tcur* = *timesprj*(*kcur*) < *tend* (see EC property p3). The incomplete EC in Equation (28) is never empty since *kcur* ≥ 1 (see EC properties p1, p3, and p5).

If *scur* > 3, we do not need to generate any time-period variates since it shows a system failure, i.e., end of the simulation in the *j*th pseudo-reality (see EC properties p8, p11, and p6).

If *scur* is 1, we need to generate two possible time-period variates: the time to failure of the primary unit, Δ*τ*1, *f p*, and the time to failure in standby of the backup unit, Δ*τ*1, *f b*. Using Equation (25):

$$
\Delta \tau\_{1,fp} = \text{PISMCF}\left(F\_1(.), F\_1^{-1}(.), T\_{cur}\right) \tag{29}
$$

$$
\Delta \tau\_{1,fb} = \text{PISMCF}\left(F\_3(.), F\_3^{-1}(.), T\_{cur}\right) \tag{30}
$$

If *scur* is 3, we do not need to generate any time-period variate since the possible time to failure of the primary unit is known to be <sup>Δ</sup>*τ*1, *f p* <sup>−</sup> <sup>Δ</sup>*τ*1, *f b* , where Δ*τ*1, *f p* and Δ*τ*1, *f b* are generated in the previous State 1 (see EC properties p9 and p13).

If *scur* > 3, we do not need to generate any time-period variate since we have observed a system failure of some type which means that the simulation in the *j*th pseudo-reality should stop and therefore *qj* = *kcur* (see EC properties p6, p8, and p11).

If *scur* is 2, we need to generate two possible time-period variates: the time to repair of the primary unit, Δ*τ*2,*rp*, and the time to failure in operation of the backup unit, Δ*τ*2, *f b*. Using Equation (25):

$$
\Delta \tau\_{2,rp} = \text{PISMCF}\left(F\_4(.), F\_4^{-1}(.), T\_{cur}\right) \tag{31}
$$

$$
\Delta \tau\_{2,fb} = \text{PISMCF}\left(F\_2(.), F\_2^{-1}(.), t\_{4\text{gc}}\right) \tag{32}
$$

If the 2SBRSBF operates with First Case distributions, the equivalent age, *tage*, of the backup unit when it starts operation at time *Tcur* is rather irrelevant since *F*2(.) is the CDF of an exponential distribution. Then we can compare the state probability functions derived by the simulational solution with the same acquired, on one hand, from numerical solution with the DAE system from Equation (16) according to the RD in Figure 1b and on the other hand with the analytical solution from Equations (9)–(12) according to the RD in Figure 1a.

If, however, the 2SBRSBF operates with Second Case distributions, then in order to use Equation (32), we have to determine the equivalent age, *tage*, at time *Tcur*. Since we need Δ*τ*2, *f b* only when the system is State 2, it follows that *kcur* is even (see EC property p9). From the beginning of the *j*th pseudo-reality up to time *Tcur*, the backup component has been in standby (*kcur*/2) times when the primary component was operating till its failure (see EC property p9). Up to *Tcur*, the backup component has never failed in standby when the primary component was in operation, i.e., during the compound time interval with overall positive length *Tsb* (see EC properties p6, p8, and p9). The latter time length can be defined using Equation (28) as:

$$T\_{sb} = \sum\_{i=1}^{k\_{\rm cur}/2} \left[ timepsr\_{\hat{\jmath}}(2i) - timepsr\_{\hat{\jmath}}(2i-1) \right] \tag{33}$$

On the other hand, the backup component has been in operation (*kcur*/2 − 1) times, when the primary component was in successful repair (see EC property p9). Up to *Tcur*, the backup component has never failed during operation when the primary component was successfully repaired, i.e., during the compound time interval with overall non-negative length *Toper* time (see EC properties p7, p8, and p9). The latter time length can be estimated by noting that up to *Tcur* the 2SBRSBF system is either in State 1 or in State 2:

$$T\_{\rm optr} = T\_{\rm cur} - T\_{\rm sb} \tag{34}$$

The non-negative value of *tage* will be the sum of the backup component operation time, *Toper*, with the equivalent operating time *Tequ* with which the backup component would age during the standby time *Tsb*:

$$t\_{\rm agc} = T\_{\rm oper} + T\_{\rm equ} \tag{35}$$

The equivalent operating time *Tequ* depends on aging in standby mechanism under which the 2SBRSBF system functions. There are three alternative assumptions for the nature of this aging in standby mechanism: full aging, no aging, or partial aging.

The full aging assumption accepts that the aging of the backup component during standby is the same as during operation (see Figure 2a):

$$T\_{equ} = T\_{sb} \Rightarrow t\_{agv} = T\_{oper} + T\_{sb} = T\_{cur} - T\_{sb} + T\_{sb} = T\_{cur} \tag{36}$$

**Figure 2.** *Cont*.

**Figure 2.** Identification of the equivalent aging time for different aging assumptions: (**a**) under full aging; (**b**) under no aging; (**c**) under partial aging.

Under the full aging assumption, a 2SBRSBF can be described with the DAE system from Equation (16) according to the RD in Figure 1b, which allows us to acquire numerical solution for Second Case distributions. The numerical solution can be compared with simulational state probability functions.

The no aging assumption accepts that the backup component during standby never ages (see Figure 2b):

$$T\_{\rm equ} = 0 \Rightarrow t\_{\rm q\%} = T\_{\rm optr} + 0 = T\_{\rm cur} - T\_{\rm sb} \tag{37}$$

Under the no aging assumption, a 2SBRSBF cannot be described with a DAE system since no RD adequately reflects the reliability behavior of the 2SBRSBF. For Second Case distributions, the only possible solution is the simulational one.

The partial aging assumption accepts that the backup component in standby ages to the same reliability as the backup component in operation during the equivalent operating time *Tequ sb* (see Figure 2c):

$$F\_{2,turn} \left( T\_{sb}^{cpu} \right) = F\_{3,turn} (T\_{sb}) \tag{38}$$

Equation (38) in simplified form was firstly proposed in [22], where it was successfully tested for Two-Component Standby System with Failures in Standby. In a real 2SBRSBF system, the failures of the backup component will be more frequent during operation than during standby which means that *F*2*,trun*(*t*) ≥ *F*3*,trun*(*t*) for any non-negative time *t.* This inequality, together with Equation (38), assures that practically always *Tequ* ∈ [0, *Tsb*]. Applying Equation (15) twice to Equation (38) we get:

$$T\_{\rm cqu} = \begin{cases} F\_2^{-1}(p^{\rm cqu}) & \text{if } p^{\rm cqu} < 1 - 100\varepsilon \\ t\_{\lambda, 2} & \text{if } p^{\rm cqu} \ge 1 - 100\varepsilon \end{cases},\\ \text{where } p^{\rm cqu} = 1 - \frac{1 - F\_2(0)}{1 - F\_3(0)}[1 - F\_3(T\_{\rm sb})] \tag{39}$$

In Equation (39), *tλ*,2 is calculated with Equation (13) for *k* = 2, therefore it uses the ideas in Equation (14) for stable approximation of the equivalent operating time, *Tequ*, at any *Tcur* ∈ [0, *tend*] for arbitrary incomplete EC from Equation (23) describing the behavior of a 2SBRSBF.

Under the partial aging assumption, the 2SBRSBF cannot be described with a DAE system since no RD adequately reflects the reliability behavior of the 2SBRSBF similarly to the no aging assumption. Again, for Second Case distributions, the only possible solution is the simulational one.

We combined Equations (33)–(37) into a six-attribute procedure, TAGEASS(.), which gives numerically stable estimates for the equivalent age, *tage*, of the backup unit under any of the three aging assumptions:

$$t\_{\rm açc} = \text{TAGEASS}\left(F\_2(.), F\_2^{-1}(.), F\_3(.), EC\_j^{\text{inc}}, T\_{\rm cur}, Flag\_A\right) \tag{40}$$

In Equation (40), the variable *FlagA* is 1, 2 or 3, respectively when the 2SBRSBF operates under the full aging, no aging, or partial aging assumptions. Then Equation (40) can be estimated using Algorithm 1.

**Algorithm 1** Equivalent Age Estimation in the *j*th Pseudo-Reality for 2SBRSBF

	- 5.1) Calculate the positive constant, *tλ*,2, using (13) with *k* = 2.
	- 5.2) Calculate the probability, *p*equ, using the second part of (39).
	- 5.3) Calculate the equivalent operating time, *T*equ, the first part of (39)

*6.3. Event Chain Generation for 2SBRSBF*

After developing the procedures for random time-period generation in Section 6.2, we may simulate an EC for the *j*th pseudo-reality of a 2SBRSBF system which satisfies all EC properties defined in Section 6.1.

The following is given:


It is easy to demonstrate that any EC generated by Algorithm 2 satisfies all EC properties formulated in Section 6.1.

The event chain for the *j*th pseudo-reality, *ECj* can be calculated using Algorithm 2.

#### **Algorithm 2** Generation of the Event Chain for the *j*th Pseudo-Reality of 2SBRSBF

	- 1.1) Set, *Tcur* =0 (the current system time is zero)
	- 1.2) Set, *kcur* = 1 (the current count of events is one)
	- 1.3) Set, *timepsrj*(*kcur*) = *Tcur* (the time of the first event is zero)
	- 1.4) Set, *statepsrj*(*kcur*) = 1 (the system starts from State 1)
	- 2.1) Estimate, Δ*τ*1, *f p* = PISMCF *F*1(.), *F*−<sup>1</sup> <sup>1</sup> (.), *Tcur* (the possible time to failure of the primary unit).
	- 2.2) Estimate, Δ*τ*1, *f b* = PISMCF *F*3(.), *F*−<sup>1</sup> <sup>3</sup> (.), *Tcur* (the possible time to standby failure of the backup unit).
		- 2.3) If *tend* ≤ *Tcur* + Δ*τ*1, *f p* and *tend* ≤ *Tcur* + Δ*τ*1, *f b* (the end of simulation comes first), then:
			- 2.3.1) Set *qj* = *kcur* (the last event count in *ECj*)
			- 2.3.2) Set *ECj* = *ECinc <sup>j</sup>* (the final *ECj*)
			- 2.3.3) Stop the Algorithm
		- 2.4) If Δ*τ*1, *f p* ≤ Δ*τ*1, *f b* (the primary unit is failing first), then:
			- 2.4.1) Set *kcur*= *kcur* +1 (new event)
			- 2.4.2) Set *Tcur*= *Tcur* + Δ*τ*2, *f b* (new current system time)
			- 2.4.3) Set *timepsrj*(*kcur*) = *Tcur* (the time of the new event)
			- 2.4.4) Generate *RN* as an evenly distributed number in the unit interval (check which is the new state) 2.4.4.1) If *RN* > *pf*, then *statepsrj*(*kcur*) = 2 (i.e., no switching failure, move to State 2) A. If *RN* ≤ *pf*, then *statepsrj*(*kcur*) = 41 (i.e., switching failure, move to State 4, type *a*)
		- 2.5) If Δ*τ*1, *f p* > Δ*τ*1, *f b* (the backup unit is failing first), then:
			- 2.5.1) Set *kcur*= *kcur* + 1 (new event)
			- 2.5.2) Set *Tcur*= *Tcur* + Δ*τ*1, *f b* (new current system time)
			- 2.5.3) Set *timepsrj*(*kcur*) = *Tcur* (the time of the new event)
			- 2.5.4) Set *statepsrj*(*kcur*) = 3 (move to State 2)
	- 3.1) Estimate, Δ*τ*2,*rp* = PISMCF *F*4(.), *F*−<sup>1</sup> <sup>4</sup> (.), *Tcur* (the possible time to repair of the primary unit).
	- 3.2) Estimate *tage* = TAGEASS *F*2(.), *F*−<sup>1</sup> <sup>2</sup> (.), *F*3(.), *ECinc <sup>j</sup>* , *Tcur*, *FlagA* (the equivalent age of the backup unit)
	- 3.3) Estimate, Δ*τ*2, *f b* = PISMCF *F*2(.), *F*−<sup>1</sup> <sup>2</sup> (.), *tage* (the possible time to operational failure of the backup unit).
	- 3.4) If *tend* ≤ *Tcur* + Δ*τ*2,*rp* and *tend* ≤ *Tcur* + Δ*τ*2, *f b* (the end of simulation comes first), then:
		- 3.4.1) Set *qj* = *kcur* (the last event count in *ECj*)
		- 3.4.2) Set *ECj* = *ECinc <sup>j</sup>* (the final *ECj*)
		- 3.4.3.) Stop the Algorithm
	- 3.5) If Δ*τ*2,*rp* ≤ Δ*τ*2, *f b* (the primary unit is repaired first), then:
		- 3.5.1) Set *kcur* = *kcur* + 1 (new event)
		- 3.5.2) Set *Tcur* = *Tcur* + Δ*τ*2,*rp* (new current system time)
		- 3.5.3) Set, *timepsrj*(*kcur*) = *Tcur* (the time of the new event)
		- 3.5.4) Generate *RN* as an evenly distributed number in the unit interval (check which is the new state) 3.5.4.1) If *RN* > *pr*, then *statepsrj*(*kcur*) =1 (no back-switching failure, move to State 1)
		- 3.5.4.2) If *RN* ≤ *pr*, then *statepsrj*(*kcur*) = 40 (back-switching failure, move to State 4, type *b*)
	- 3.6) If Δ*τ*2,*rp* > Δ*τ*2, *f b* (the backup unit is failing first), then:
		- 3.6.1) Set, *kcur* = *kcur* + 1 (new event)
		- 3.6.2) Set, *Tcur* = *Tcur* + Δ*τ*2, *f b* (new current system time)
		- 3.6.3) Set, *timepsrj*(*kcur*) = *Tcur* (the time of the new event)
		- 3.6.4) Set, *statepsrj*(*kcur*) =42 (move to State 4, type *c*)
	- 4.1) If *tend* ≤ *Tcur* + Δ*τ*1, *f p* − Δ*τ*1, *f b* (the end of simulation comes first), then:
		- 4.1.1) Set, *qj*= *kcur* (the last event count in *ECj*)
		- 4.1.2) Set, *ECj* = *ECinc <sup>j</sup>* (the final *ECj*)
		- 4.1.3) Stop the Algorithm
	- 4.2) If *tend* > *Tcur* + Δ*τ*1, *f p* − Δ*τ*1, *f b* (the primary unit is failing first), then:
		- 4.2.1) Set, *kcur*= *kcur*+1 (new event)
		- 4.2.2) Set, *Tcur*= *Tcur*+Δ*τ*1, *f p* − Δ*τ*1, *f b* (new current system time)
		- 4.2.3) Set, *timepsrj*(*kcur*) = *Tcur* (the time of the new event)
		- 4.2.4) Set, *statepsrj*(*kcur*) = 43 (switching failure, move to State 4, type *d*)
	- 5.1) Set, *qj* = *kcur* (the last event count in *ECj*)
		- 5.2) Set, *ECj* = *ECinc <sup>j</sup>* (the final *ECj*)
		- 5.3) Stop the Algorithm

### *6.4. Extracting Reliability Information from the Simulated ECs*

Let *N* be a large positive integer representing the count of the randomly simulated pseudo-realities. Using Algorithm 2, we can simulate *ECj*, for *j* = 1,2, ... , *N*. In this section, we will extract the reliability information from the simulated ECs, approach which is the essence of any Monte Carlo simulation [37] (pp. 290–294).

Let us calculate the state probability functions at the 2000 evenly distributed time points from 0 to *tend* given in Equation (20). For a given *ECj* we can estimate the state, *Sti,j*, at each of the time points *ti*:

$$St\_{i,j} = \begin{cases} statepsr\_j(k) \text{ if } timepsr\_j(k) \le t\_i < timepsr\_j(k+1) \text{, for } k < q\_j\\ statepsr\_j(q\_j) \text{ if } timepsr\_j(q\_j) \le t\_i \le t\_{end} \end{cases}, \text{ where } \begin{array}{l} i = 1, 2, \dots, 2000\\ j = 1, 2, \dots, N \end{array} \tag{41}$$

From Equation (41) it is easy to estimate the values of the first three state probability functions at the time point, *ti*:

$$P\_{\mathcal{S}}(t\_i) = \frac{1}{N} \# \{ \mathcal{S}\_{i,j} = \mathcal{g} \, | \, j = 1, 2, \dots, N \}, \text{ where } \mathcal{g} = 1, 2, 3 \text{ and } i = 1, 2, \dots, 2000 \tag{42}$$

In Equation (42) the # *Si*,*<sup>j</sup>* = *g <sup>j</sup>* <sup>=</sup> 1, 2, . . . , *<sup>N</sup>* is the count of all states at the time point *ti* which are equal to *g*.

The fourth state probability function can be estimated using Equation (1) as:

$$P\_4(t\_i) = 1 - P\_3(t\_i) - P\_2(t\_i) - P\_3(t\_i) \text{ , for } i = 1, 2, \dots, 2000 \tag{43}$$

The reliability function and the*MTTFsys* can be approximated with Equations (21) and (22). According to the ES property p1, the reliability in Equation (22) has decreasing nodes:

$$R\_{sys}(t\_i) \ge R\_{sys}(t\_{i+1}) \text{ , for } i = 1, 2, \dots, 1999 \tag{44}$$

One way to identify the *α*-design life, *tdes,<sup>α</sup>* for given *α* is to transform the nodes, *ti*, *Rsys*(*ti*) *<sup>i</sup>* <sup>=</sup> 1, 2, . . . , 2000 , of the system reliability from Equation (22) into strictly decreasing purged nodes !*t pu <sup>i</sup>* , *<sup>R</sup>pu sys t pu i* " *i* = 1, 2, . . . , *npu*# where:

$$R\_{sys}^{pu} \left( t\_i^{pu} \right) > R\_{sys}^{pu} \left( t\_i^{pu} \right), \text{ for } i = 1, 2, \dots, n^{pu} \tag{45}$$

Such a purging procedure is proposed in [17], where the algorithm is motivated, formalized, illustrated, and proven. In short, it runs in the steps summarized in Algorithm 3.

**Algorithm 3** Purging Algorithm

(a) Identify the time of the first purged node *t pu* <sup>1</sup> , *<sup>R</sup>pu sys t pu* 1 = 1 as the greatest *ti* for which *Rsys*(*ti*) = 1;


Having the strictly decreasing purged system reliability function, we can identify the *α*-design life, *tdes,<sup>α</sup>* for any *α* ∈ ! *Rpu sys t pu npu* , 1" :

$$t\_{\rm des,a} = t\_i^{pu} + \left[ R\_{\rm sys}^{pu} \left( t\_i^{pu} \right) - a \right] \frac{t\_{i+1}^{pu} - t\_i^{pu}}{R\_{\rm sys}^{pu} \left( t\_i^{pu} \right) - R\_{\rm sys}^{pu} \left( t\_{i+1}^{pu} \right)}, \text{ for } R\_{\rm sys}^{pu} \left( t\_i^{pu} \right) \ge a > R\_{\rm sys}^{pu} \left( t\_{i+1}^{pu} \right) \tag{46}$$

As discussed in Section 3, the reliability numerical characteristics *Mediansys*, *B*1*\_life*, *B*10*\_life*, and *IQRsys* can be estimated as *tdes*,0.5, *tdes*,0.99, *tdes*,0.9, and *tdes*,0.25 − *tdes*,0.75 respectively by applying Equation (46) five times.

The simulational solution is universal and exists even when the numerical and analytical solutions are impossible. Even when the numerical and the analytical solutions exist, the simulational solution can provide richer reliability information.

For example, it is obvious that the 2SBRSBF system will have 100% chance to ever be in the State 1. It is also clear that if tend is correctly selected, then the 2SBRSBF system will have more than 99% chance to ever be in the State 4. However, it is interesting to know the chance, *P*2*,ever*, for the 2SBRSBF system to ever be in the State 2, since that probability will help us plan the resources needed for the repair of the primary unit. Similarly, the chance, *P*3*,ever*, for the 2SBRSBF system to ever be in the State 3 is important, since that will show us the prevalence of the failure in standby of the backup unit. So, for a given 2SBRSBF system, we can estimate the chances, *Pg,ever*, for *g* = 1,2,3:

$$P\_{\text{g.cvar}} = \frac{100}{N} \# \{ \exists i, \text{ that } \mathbb{S}\_{i,j} = \emptyset | j = 1, 2, \dots, N \}, \text{ where } \mathbb{g} = 1, 2, 3 \tag{47}$$

In Equation (47), # ∃*i*, that *Si*,*<sup>j</sup>* = *g <sup>j</sup>* <sup>=</sup> 1, 2, . . . , *<sup>N</sup>* is the count of pseudo-realities in which State *g* can be found at least once. Similarly, for a given 2SBRSBF system we can estimate the chance, *P*4*,ever* as:

$$P\_{4, \text{rev}} = \frac{100}{N} \# \left( \exists i\_{\text{\\_that } S\_{i,j}} > 3 \, | \, j = 1, 2, \dots, N \right) \tag{48}$$

In Equation (48), # ∃*i*, that *Si*,*<sup>j</sup>* > 3 *<sup>j</sup>* <sup>=</sup> 1, 2, . . . , *<sup>N</sup>* is the count of pseudo-realities in which State 4 (system failure) can be found at least once.

As another example for reliability information, which can be acquired neither with the numerical, nor with the analytical solution, can be found in the four conditional chances, *Pcond <sup>g</sup>*,*ever* (for *g* = 40, 41, 42, 43), of the 2SBRSBF system to have respectively type *b*, type *a*, type *c*, or type *d* system failure, provided that system has failed:

$$P\_{\text{g, rev}}^{\text{cond}} = 100 \frac{\#\left(\mathcal{S}\_{i, \eta\_j} = \mathcal{g} \, \middle| \, j = 1, 2, \dots, N\right)}{N P\_{4, \text{rev}} / 100}, \text{ where } \mathcal{g} = 40, 41, 42, 43 \tag{49}$$

The information in Equation (49) allows to identify the types of system failures which dominate the 2SBRSBF system. That knowledge will increase the efficiency of the reliability improvement measures. Equations (42), (47)–(49) use the frequentist interpretation of probability [40] (pp. 42–43).

Knowing how to simulate an EC for the *j*th pseudo-reality of a 2SBRSBF system, allows us to develop the simulational solution of a given 2SBRSBF system. We have the following given:


The proposed algorithm in [17] uses simulation to find the reliability characteristics of a two-component standby systems with switching failures and aging in standby. The simulational solution for 2SBRSBF system can be obtained through a generalization of that algorithm, which is formalized as Algorithm 4 below.

With the formulation of Algorithm 4 the universal simulational solution for a 2SBRSBF system is complete.

#### **Algorithm 4** Simulational Solution of a 2SBRSBF System


#### **7. Illustrative Examples**

#### *7.1. Examples Setup*

We shall analyze three Illustrative Examples. In all of them, the probability for switching failure is *pf* = 0.12, whereas the probability for back-switching failure is *pr* = 0.03. The ratio between those values is plausible for the following reasons. If the switching is successful, it means that the switching device operated properly. Then a back-switching failure is less probable since it will be demanded shortly afterwards (the repair time of the primary component is much smaller than its failure time).

In Example 1, any of the four original distributions has a constant failure/repair rate *λ<sup>k</sup>* shown in Table 2 (for *k* = 1,2,3,4). The PDFs of the original exponential distributions are:

$$f\_k(t) = \lambda\_k e^{-\lambda\_k t}, \text{ for } t \ge 0 \text{ where } k = 1, 2, 3, 4 \tag{50}$$

**Table 2.** Description of the original distributions in Example 1.


The PDFs, the reliability/repair functions, and the failure/repair rates of the truncated distributions from Equation (50) are plotted in Figure 3. Example 1 will illustrate the behavior of the 2SBRSBF system with First Case distributions. Here, the original and the truncated distributions coincide.

In Example 2, the original distributions are as follows:

(1) a Rayleigh distribution with shape parameter *b*<sup>1</sup> [41] for the failures of the primary component:

$$f\_1(t) = (t/b\_1)e^{-0.5(t/b\_1)^2}, \text{ for } t \ge 0\tag{51}$$

(2) a normal distribution with mean value *μ*<sup>2</sup> h and standard deviation *σ*<sup>2</sup> h [42] for the failures of the backup component in operation:

$$f\_2(t) = \frac{1}{\sqrt{2\pi}\sigma\_2} e^{-0.5(t-\mu\_2)^2/(\sigma\_2)^2}, \text{ for } t \in \left(-\infty, +\infty\right) \tag{52}$$

(3) a Weibull distribution with a scale parameter *θ*<sup>3</sup> h and a shape parameter *β*<sup>3</sup> [43] for the failures of the backup component in standby:

$$f\_3(t) = \frac{\beta\_3}{\theta\_3} \left(\frac{t}{\theta\_3}\right)^{\beta\_3 - 1} e^{-(t/\theta\_3)^{\beta\_3}}, \text{ for } t \ge 0 \tag{53}$$

(4) a lognormal distribution with median time *tmed*,4 h and shape parameter *s*<sup>4</sup> [44] for the repairs of the primary component:

$$f\_4(t) = \frac{1}{\sqrt{2\pi} \text{s}\_4 t} e^{-0.5 \ln^2 \left(t/t\_{mol,4}\right)/\left(s\_4\right)^2}, \text{ for } t \ge 0 \tag{54}$$

The original distribution Example 2 are described in Table 3. The PDFs, the reliability/repair functions, and the failure/repair rates of the truncated distributions from Equations (51)–(54) are plotted in Figure 4. Example 2 will illustrate the behavior of the 2SBRSBF system with Second Case distributions where the failures of the backup component in operation have an Increasing Failure Rate (IFR). Such a typical situation can occur when the operational failure is caused mainly by high wearing in the backup component [11] (pp. 73–75). Here, the original and the truncated distributions coincide except for the *f* 2(*t*) and *f* 2,*trunc*(*t*).

**Figure 3.** The truncated distributions in Example 1. The three failure distributions are shown in section (**a**–**c**), whereas the repair distribution is shown in section (**d**–**f**). The reliability/repair functions, the PDFs and the failure/repair rates are given respectively in the first (sections (**a**,**d**)), the second (section (**b**,**e**)), and the third row (sections (**c**,**f**)).

**Table 3.** Description of the original distributions in Example 2.

**Figure 4.** The truncated distributions in Example 2. The three failure distributions are shown in section (**a**–**c**), whereas the repair distribution is shown in section (**d**–**f**). The reliability/repair functions, the PDFs and the failure/repair rates are given respectively in the first (sections (**a**,**d**)), the second (section (**b**,**e**)), and the third row (sections (**c**,**f**)).

In Example 3 the distributions are the same as in Example 2, except for the second type, which changes to:

2) a lognormal distribution with median time *tmed*,2 h and shape parameter *s*<sup>2</sup> for the failures of the backup component in operation:

$$f\_2(t) = \frac{1}{\sqrt{2\pi}s\_2t}e^{-0.5\ln^2\left(t/t\_{mol2}\right)/\left(s\_2\right)^2}, \text{ for } t \ge 0\tag{55}$$

The original distribution Example 3 are described in Table 4. The PDFs, the reliability/repair functions, and the failure/repair rates of the truncated distributions from Equations (51), (53)–(55) are plotted in Figure 5. Example 3 will illustrate the behavior of the 2SBRSBF system with Second Case distributions where the failures of the backup component in operation have a Decreasing Failure Rate (DFR). Such an atypical situation can occur when the operational failure is caused mainly by high child mortality in the backup component [11] (pp. 73–75). Here, the original and the truncated distributions coincide.

**Table 4.** Description of the original distributions in Example 3.

**Figure 5.** The truncated distributions in Example 3. The three failure distributions are shown in section (**a**–**c**), whereas the repair distribution is shown in section (**d**–**f**). The reliability/repair functions, the PDFs and the failure/repair rates are given respectively in the first (sections (**a**,**d**)), the second (section (**b**,**e**)), and the third row (sections (**c**,**f**)).

#### *7.2. Example 1 Solution*

Since in Example 1, we are dealing with First Case distributions, the type of aging has no effect on the reliability performance of the 2SBRSBF system. The simulation solution was obtained by Algorithm 4 with *N* = 10,000 pseudo-realities for time from 0 to *tend* = 20,000 h. Four typical pseudo-realities are shown in Table 5 where the different types of system failures are demonstrated. The four state probability functions are shown in Figure 6a–d, respectively. The system reliability function is depicted in Figure 7. The simulation reliability at *tend* was negligible (as required *R*sys(20,000) = 0.0024 < 0.01) which justifies the selection of *tend*. Important simulation numerical characteristics of the 2SBRSBF reliability can be found in Table 6. The chances of some events of interest (described in Section 6.4) can be found in Table 7. It is revealing so see that the backup component has approximately 69% chance to endure failure in standby (State 3). Another useful fact is that the switching failures (Type *a*) are more frequent than the backup component failures in operation (Type *c*) (17% vs. 11% conditional chance). That fact suggests that it is easier to improve the reliability by upgrading the switching mechanism than by upgrading the backup unit.


**Table 5.** Four typical pseudo-realities from Example 1.

**Figure 6.** State probability functions for Example 1 (with states 1 through 4 given in sections (**a**–**d**) respectively) from the analytical, numerical and simulation solution.


**Table 6.** Reliability characteristics of the 2SBRSBF from Example 1.

**Table 7.** Chances for events of interest in % for Example 1.


The simulation results were verified by comparison with the precise analytical solution (see Section 4). According to Table 6, the precise analytical MTTF is 4282 h, whereas the simulational MTTF is estimated as 4294 h, which contains less than 0.3% error.

Also, the simulational results were verified by comparison with the numerical solution (see Section 5), which, as seen from Figures 6 and 7, produced undistinguishable curves from the simulational state probabilities and the simulational reliability function. According to Table 6, the numerical MTTF is 4280 h, whereas the simulational MTTF is estimated as 4294 h. The numerical solution for Example 1 (as well as in Examples 2 and 3) was derived by solving the index-1 DAE system described in Section 5 with the MATLAB multistep procedure *ode15s.m*. The software successfully integrated the DAE system from 0 to *tend* = 20,000 h using variable-step method of variable order from 1 to 5 [45].

As seen from Figures 6 and 7, the analytical and the numerical solutions produce undistinguishable curves from the simulational state probabilities and the simulational reliability function. The observed overlap is an essential part of the verification of the presented simulation algorithm: in the case of exponential distribution, the model is Markovian, where the analytical, the numerical, and the simulational solutions should practically coincide.

#### *7.3. Example 2 Solution*

Since Example 2 deals with Second Case distributions, the type of aging has an effect on the reliability performance of the 2SBRSBF system. Three simulational solutions were obtained by repeatedly using Algorithm 4 with *N* = 10,000 pseudo-realities for the three aging assumptions: full aging, no aging and patrial aging of the backup component in standby. Each of those solutions was estimated for time from 0 to *tend* = 8000 h. The three sets of four state probability functions are shown in Figure 8a–d, respectively. The three system reliability functions are depicted in Figure 9. The simulational reliabilities at *tend* were negligible and much lower than 0.01 (for full aging-*R*sys(8000) = 0; for no aging-*R*sys(8000) = 0.0025; for partial aging-*R*sys(8000) = 0.0003) which justifies the selection of *tend*.

**Figure 8.** State probability functions for Example 2 (with states 1 through 4 given in sections (**a**–**d**) respectively) under the three aging assumptions.

**Figure 9.** Reliability functions for Example 2 under the three aging assumptions.

Important simulational numerical characteristics of the 2SBRSBF reliabilities can be found in Table 8 for the three types of aging. The chances of some events of interest (described in Section 6.4) can be found in Table 9 for each of the three aging assumptions. It is revealing to see that the backup component has between 31% and 41% chance to endure failure in standby (State 3) depending on the aging model. An interesting dynamic is observed in the conditional chances of observing the different types of failure. At full aging, the backup component failures during primary repair (Type *c*) have more than 50% chance, whereas the primary component failures after failure in standby (Type *d*) constitute only around 30% of the system failures. At no aging, the backup failures in operation are less likely and, therefore, the primary component failures after backup failure in standby (Type *d*) are more frequent than the backup component failures during primary repair (Type *c*) (41% vs. 36% conditional chance). At the same time, Type *c* and Type *d* system failures are marginally the same at partial aging of the backup component in standby (37% vs. 42% conditional chance). Those facts suggest that to increase the reliability of the 2SBRSBF it is of paramount importance correctly to identify the aging mechanism of the backup unit during standby.


**Table 8.** Reliability characteristics of the 2SBRSBF from Example 2 under the three aging assumptions.


**Table 9.** Chances in % for events of interest for Example 2 under the three aging assumptions.

In Example 2, the distribution of the backup component failures in operation has an IFR (see the blue line in Figure 5c), indicating that the wear out is the most likely reason for those failures. This is by far the most widespread case in the engineering practice where the backup component operates at the rear end of the bathtub curve [46]. Then, the severity of the aging should increase the failure incidence of the backup component in operation and subsequently should decrease the reliability. As expected, the system reliability function is the best at no-aging and worst at full aging (see Figure 9 for 1500–5500 h). The MTTF increases from 2837 h at full aging, through 3242 h at partial aging, to 3457 h at no aging, which corresponds to substantial 21% improvement. Similar behavior can be observed in the median, B10 life, and at the B1 life (see Table 8). Another expected result is that the state probability functions for partial aging are between the state probability functions for no aging and full aging (see Figure 8). The real distinction between the three curves can be seen in State 2 probability function (Figure 8b) which is very sensitive to the aging mode. The observed forms of the State 2 probability functions are justifiable since the severity of aging increases the incidence of failure of the operational backup unit, which moves the system to State 4 and decreases the probability of the 2SBRSBF to be in State 2. All the above can serve as a qualitative validation of Algorithm 4 for simulating the reliability behavior of the 2SBRSBF system.

Also, the simulational results were quantitatively verified by comparison with the numerical solution (as described in Section 7.2), which, as seen from Figures 8 and 9, produced undistinguishable curves from the simulational state probabilities and the simulational reliability function in the case of full aging of the backup component during standby. This overlap is an important result: under the full-aging assumption the model is semi-Markovian where the numerical, and the simulational solutions should practically coincide. According to Table 8, the numerical MTTF and the simulational MTTF at full aging are estimated to be equal (2837 h). Note that the analytical solution is impossible to be derived in Example 2 since the failure/repair rates are not constant.

#### *7.4. Example 3 Solution*

Since Example 3 deals with Second Case distributions, similarly to Example 2, the type of aging has effect on the reliability performance of the 2SBRSBF system. Three simulational solutions were obtained by repeatedly using Algorithm 4 with *N* = 10,000 pseudo-realities for the three aging assumptions: full aging, no aging and patrial aging of the backup component in standby. Each of those solutions was estimated for time from 0 to *tend* = 12,000 h. The three sets of four state probability functions are shown in Figure 10a–d, respectively. The three system reliability functions are depicted in Figure 11. The simulational reliabilities at *tend* were negligible and lower than 0.01 (for full aging-*R*sys(12000) = 0.0062; for

no aging-*R*sys(12,000) = 0.0011; for partial aging-*R*sys(12,000) = 0.002) which justifies the selection of *tend*.

**Figure 10.** State probability functions for Example 3 (with states 1 through 4 given in sections (**a**–**d**) respectively) under the three aging assumptions.

**Figure 11.** Reliability functions for Example 3 under the three aging assumptions.

Important simulational numerical characteristics of the 2SBRSBF reliabilities can be found in Table 10 for the three types of aging. The chances of some events of interest (described in Section 6.4) can be found in Table 11 for each of the three aging assumptions.


B1 life 4.870 × 10+2 h 5.031 × <sup>10</sup>+2 h 4.999 × <sup>10</sup>+2 <sup>h</sup> Mean value (Analytical) NA NA NA Mean value (Numerical) 3.653 × 10+3 <sup>h</sup> NA NA

**Table 10.** Reliability characteristics of the 2SBRSBF from Example 3 under the three aging assumptions.

**Table 11.** Chances in % for events of interest for Example 3 under the three aging assumptions.


In Example 3, the distribution of the backup component failures in operation has an DFR (see the blue line in Figure 5c), indicating that the child mortality is the most likely reason for those failures. This is a very rare case in the engineering practice where the backup component operates at the front end of the bathtub curve. Then, the severity of the aging should decrease the failure incidence of the backup component in operation and subsequently should increase the reliability. As expected, the system reliability function is the worst at no-aging and best at full aging (see Figure 11 for 2000–8000 h). The MTTF increases from 3139 h at no aging, through 3187 h at partial aging, to 3625 h at full aging, which corresponds to noticeable 16% improvement. Similar behavior can be observed in the median, B10 life, and at the B1 life (see Table 10). Another expected result is that the state probability functions for partial aging are between the state probability functions for no aging and full aging (see Figure 10). The real distinction between the three curves can be seen in State 2 probability function (Figure 10b) which is very sensitive to the aging mode. The observed forms of the State 2 probability functions are justifiable since the severity of aging decreases the incidence of failure of the operational backup unit, which moves the system to State 4 and increases the probability of the 2SBRSBF to be in State 2. All the above can serve as a qualitative validation of Algorithm 4 for simulating the reliability behavior of the 2SBRSBF system.

A partial overlap between the no aging simulation solution and the partial aging simulation solution can be spotted in Figures 10 and 11. The same can also be observed in Figures 8 and 9 to a lesser extent. Those partial overlaps reflect the fact that for almost all realistic distribution sets, under the applied method, the solution of partial aging is much closer to the solution with no aging assumption than to the solution with full aging assumption.

Again, the simulational results were quantitatively verified by comparison with the numerical solution (as described in Section 7.2) which as seen from Figures 10 and 11 produced undistinguishable curves from the simulational state probabilities and the simulational reliability function in the case of full aging of the backup component during standby (for comment on the observed overlap see Section 7.3). According to Table 10, the numerical MTTF and the simulational MTTF at full aging are estimated to be virtually equal (3653 h vs. 3652 h, respectively). Note that the analytical solution is impossible to be derived in Example 2 since the failure/repair rates are not constant.

#### **8. Conclusions**

In this paper, we investigated the reliability effect of introducing a primary component minimal repair in a two-component standby system with switching failures and aging in warm-standby. A novel analytical solution was derived for distributions with constant failure/repair rates. Under a full aging assumption of the backup component during standby, an index-1 DAE system of four simultaneous equations with constant mass singular matrix was proposed and solved to numerically approximate the state probability functions and system reliability. A universal simulational algorithm was designed to solve the 2SBFSR system under three types of aging. That algorithm generates pseudo-realities with ECs, which satisfy the newly formulated EC properties for the 2SBFSR system. Novel function to assess the equivalent age of the backup component under arbitrary aging mechanisms was proposed and utilized during the EC generation. The system has a stable operation with any type of distribution. There is a significant practical benefit in the ability of the user to write their own distribution functions, which reflect several modes of failure during operation, several modes of failure during warm-standby, and several modes of repair.

Three numerical examples were elaborated to validate quantitatively and qualitatively the simulational solution. To model the 2SBRSBF system with partial aging in standby, we assumed that that the backup component in standby ages to the same reliability as the backup component in operation. That is a logical and plausible hypothesis that allows to produce a tractable aging model whose results can be treated as best estimate. Even if the real aging mechanism is different the numerical examples show that the partial aging results always will be bounded by the full aging and the no-aging results. That fact allows the designers and the maintenance staff to correctly assess the effect of alternative measures aiming at improving the system reliability even if the precise aging in standby mechanism is known.

Although our model may look too specific and simplified, it is easily scalable. The demonstrated methodology can easily be applied to multiple-component warm-standby system with random configuration. We have not given such an example for purely volume constraints in this work. Any different aging assumptions can be incorporated by modifying Algorithm 1 (hence the function TAGEASS). All aspects and elements of such a multi-component warm-standby system can be found in 2SBRSBF. In such a way, our model is suitable for applications in industrial systems, manufacturing, design of ship electrical and propulsion systems, power plants, etc.

As a direction for future studies, we may study the ways to adapt our procedures to the case of perfect repair [10] and intermediate repair [8], as this work only analyzed the case of minimal repair.

**Author Contributions:** Conceptualization, K.T. and N.N.; methodology, K.T., N.N. and S.C.; software, K.T. and S.C.; validation, S.C. and G.F.; investigation, K.T., B.M. and G.F.; data curation, B.M. and S.C.; writing—original draft, K.T., S.C. and B.M.; writing—review and editing, N.N., G.F. and B.M.; visualization, K.T., N.N. and S.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data sharing not applicable. No new data were created or analyzed in this study. Data sharing is not applicable to this article.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **Appendix A**

**Given:** Let the random variable *T* be the time to failure (or repair) of a component. Also, let the random event *A*(*T*0) be that the component is operational at, or has not been repaired up to, time *T*0. In fact, *T* is the deterministic time *T*<sup>0</sup> plus the random time period till the next event (failure or repair). This definition of *T* is true only for Appendix A. Then:


**Prove:** The unconditional and the conditional failure (repair) rate are equal for any *t*∗ ≥ *T*0, i.e., *λ*(*t*∗) = *λcond*(*t* ∗ −*T*0|*T*0), for *t*∗ ∈ [*T*0, ∞).

**Proof.** The unconditional *R*(*t*) and *f*(*t*) are given in Figure A1a,c. The relationship between these functions is:

$$f(t) = \frac{dF(t)}{dt} = \frac{d[1 - R(t)]}{dt} = -\frac{dR(t)}{dt} \text{ for } t \in [0, \infty) \tag{A1}$$

Similarly, the conditional *Rcond*(*τ*|*T*0) and *fcond*(*τ*|*T*0) are given in Figure A1b,d. The relationship between these functions is:

$$f\_{\rm cond}(\boldsymbol{\tau}|T\_0) = \frac{d F\_{\rm cond}(\boldsymbol{\tau}|T\_0)}{d\boldsymbol{\tau}} = \frac{d \left[1 - R\_{\rm cond}(\boldsymbol{\tau}|T\_0)\right]}{d\boldsymbol{\tau}} = -\frac{d R\_{\rm cond}(\boldsymbol{\tau}|T\_0)}{d\boldsymbol{\tau}} \text{ for } \boldsymbol{\tau} = (t - T\_0) \in [0, \infty) \tag{A2}$$

According to [11] (p. 72), the value of the conditional *Rcond*(*τ*|*T*0) can be expressed as the ratio of two unconditional values of *R*(*t*):

$$R\_{cond}(\tau|T\_0) = \frac{R(\tau + T\_0)}{R(T\_0)} \text{ for } \tau = (t - T\_0) \in [0, \infty) \tag{A3}$$

The interdependency between Figure A1a,b illustrates Equation (A3). The constant *R*(*T*0) is the height of the red vertical line in Figure A1a.

**Figure A1.** A generic distribution described by: (**a**) unconditional reliability; (**b**) conditional reliability; (**c**) unconditional density; (**d**) conditional density.

Let us take the first derivative about *τ* from Equation (A3) and multiply both sides by negative 1. Then,

$$-\frac{d\mathcal{R}\_{cond}(\mathbf{r}|T\_0)}{d\tau} = -\frac{d}{d\tau}\frac{\mathcal{R}(\tau+T\_0)}{\mathcal{R}(T\_0)}\text{ for}\tau = (t-T\_0)\in[0,\infty)\tag{A4}$$

Let us simplify the RHS of Equation (A4) using Equation (A1) and utilizing that *τ* = (*t* − *T*0):

$$\begin{array}{rcl}-\frac{d}{d\tau}\frac{R(\tau+T\_0)}{R(T\_0)} &=& -\frac{1}{R(T\_0)}\frac{dR(\tau+T\_0)}{d\tau} = -\frac{1}{R(T\_0)}\frac{dR(\tau+T\_0)}{d(\tau+T\_0)}\frac{d(\tau+T\_0)}{d\tau} \\ &=& -\frac{1}{R(T\_0)}\frac{dR(t)}{dt}\frac{dt}{d\tau} = \frac{1}{R(T\_0)}f(t)\frac{d(\tau+T\_0)}{d\tau} \\ &=& \frac{f(\tau+T\_0)}{R(T\_0)}(1) = \frac{f(\tau+T\_0)}{R(T\_0)}\end{array} \tag{A.5}$$

According to Equation (A2), the LHS of Equation (A4) is *fcond*(*τ*|*T*0). Then, from Equations (A4) and (A5) it follows that:

$$f\_{cond}(\tau|T\_0) = \frac{f(\tau + T\_0)}{R(T\_0)} \text{ for } \tau = (t - T\_0) \in [0, \infty) \tag{A6}$$

The interdependency between Figure A1c,d illustrates Equation (A6). The constant *R*(*T*0) is the area of the green patch in Figure A1c, since from Equation (A1) it follows that *R*(*T*0) = ∞ *f*(*t*)*dt*.

The conditional failure (repair) rate of *T* if *A*(*T*0) can be transformed using Equations (A3) and (A6):

$$\lambda\_{\text{cond}}(\mathbf{r}|T\_0) = \frac{f\_{\text{cond}}(\mathbf{r}|T\_0)}{R\_{\text{cond}}(\mathbf{r}|T\_0)} = \frac{f(\mathbf{r} + T\_0)}{R(T\_0)} \div \frac{R(\mathbf{r} + T\_0)}{R(T\_0)} = \frac{f(\mathbf{r} + T\_0)}{R(\mathbf{r} + T\_0)} \text{ for } \mathbf{r} = (t - T\_0) \in [0, \infty) \tag{A7}$$

Let's select a time point *t*∗ ≥ *T*0. The unconditional failure (repair) rate of *T* at time *t\** is:

$$
\lambda(t\*) = \frac{f(t\*)}{R(t\*)}\tag{A8}
$$

The nominator and the denominator in Equation (A8) are respectively the heights of the blue lines in Figure A1a,c. The conditional time *τ* is simply the time *t* delayed with *T*<sup>0</sup> (i.e., *t* = *τ* + *T*0).

From here, the relative time moment *τ*∗ which coincides with time *t\** is:

$$
\tau \ast = t \ast - T\_0 \tag{A9}
$$

Equation (A9) is illustrated by the transition from Figure A1a to Figure A1b, and in the transition from Figure A1c to Figure A1d.

The value of *λcond*(*τ*|*T*0) at relative time point *τ*∗ can be easily calculated from Equation (A7) utilizing Equations (A8) and (A9):

$$\begin{array}{rcl} \lambda\_{\text{cond}}(t \ast -T\_0 | T\_0) &=& \frac{f\_{\text{cond}}(t \ast -T | T\_0)}{\overline{R}\_{\text{cond}}(t \ast -T | T\_0)} = \lambda\_{\text{cond}}(\tau \ast | T\_0) \\ &=& \frac{f(\tau \ast +T\_0)}{\overline{R}(\tau \ast +T\_0)} = \frac{f(t \ast)}{\overline{R}(t \ast)} = \lambda(t \ast) \end{array}, \text{ for} \mathbf{f} \ast \in [T\_0, \infty) \tag{A10}$$


### **Appendix B**

*T*0

**Given:** Let *λ*1, *λ*2, *λ*3, and *λ*<sup>4</sup> be real positive constants, whereas *pr* and *pf* are real positive constants less than 1. The real functions *P*1(*t*), *P*2(*t*), and *P*3(*t*) are defined in the Domain *t* ∈ [0, ∞) and satisfy the system from Equation (A11) of three simultaneous ordinary differential equations. The initial conditions of the functions are given in Equation (A12).

$$\begin{cases} \frac{dP\_1}{dt}(t) = -(\lambda\_1 + \lambda\_3)P\_1(t) + (1 - p\_r)\lambda\_4 P\_2(t) \\ \frac{dP\_2}{dt}(t) = \left(1 - p\_f\right)\lambda\_1 P\_1(t) - (\lambda\_4 + \lambda\_2)P\_2(t) \\ \frac{dP\_3}{dt}(t) = \lambda\_3 P\_1(t) - \lambda\_1 P\_3(t) \end{cases} \tag{A11}$$

$$P\_1(0) = 1,\\ P\_2(0) = 0,\\ P\_3(0) = 0 \tag{A12}$$

**Find:**


$$\text{(c)}\quad \text{The quantity } MTTF\_{sys} = \int\_0^\infty R\_{sys}(t)dt.$$

#### **Solution:**

(a) Taking Laplace transformation [47] (pp. 331–335) of the three equations in Equation (A11) yields a system of three algebraic equations about the Laplace transforms *Y*1(*s*), *Y*2(*s*), and *Y*3(*s*) of the functions *P*1(*t*), *P*2(*t*), and *P*3(*t*), where *s* is a complex number known as frequency:

$$\begin{cases} sY\_1(s) - P\_1(0) = -(\lambda\_1 + \lambda\_3)Y\_1(s) + (1 - p\_r)\lambda\_4 Y\_2(s) \\ sY\_2(s) - P\_2(0) = (1 - p\_r)\lambda\_1 Y\_1(s) - (\lambda\_4 + \lambda\_2)Y\_2(s) \\ sY\_3(s) - P\_3(0) = \lambda\_3 Y\_1(s) - \lambda\_1 Y\_3(s) \end{cases} \tag{A13}$$

Substituting Equation (A12) into Equation (A13) and simplifying gives:

$$\begin{cases} (s + \lambda\_1 + \lambda\_3)\mathbf{Y}\_1(s) - (1 - p\_r)\lambda\_4 \mathbf{Y}\_2(s) = 1 \\ - \left(1 - p\_f\right)\lambda\_1 \mathbf{Y}\_1(s) + (s + \lambda\_2 + \lambda\_4)\mathbf{Y}\_2(s) = 0 \\ - \lambda\_3 \mathbf{Y}\_1(s) + (\lambda\_1 + s)\mathbf{Y}\_3(s) = 0 \end{cases} \tag{A14}$$

The first two equations in Equation (A14) can be solved for *Y*1(*s*), *Y*2(*s*) using the Cramer's rule [48]:

$$Y\_1(s) = \frac{s + \lambda\_2 + \lambda\_4}{(s + \lambda\_1 + \lambda\_3)(s + \lambda\_2 + \lambda\_4) - (1 - p\_r)\left(1 - p\_f\right)\lambda\_1\lambda\_4} \tag{A15}$$

$$Y\_2(s) = \frac{\left(1 - p\_f\right)\lambda\_1}{(s + \lambda\_1 + \lambda\_3)(s + \lambda\_2 + \lambda\_4) - (1 - p\_r)\left(1 - p\_f\right)\lambda\_1\lambda\_4} \tag{A16}$$

The denominator in both Equations (A15) and (A16) is a quadratic polynomial with real coefficients 1, *K*, and *C*:

$$(s + \lambda\_1 + \lambda\_3)(s + \lambda\_2 + \lambda\_4) - (1 - p\_\mathcal{I})\left(1 - p\_\mathcal{f}\right)\lambda\_1\lambda\_4 = s^2 + 2\mathcal{K}s + \mathcal{C} \tag{A17}$$

where the real constants *K* and *C* are:

$$\begin{aligned} K &= (\lambda\_1 + \lambda\_2 + \lambda\_3 + \lambda\_4)/2 \\ C &= (\lambda\_1 + \lambda\_3)(\lambda\_2 + \lambda\_4) - \left(1 - p\_f\right)(1 - p\_r)\lambda\_1\lambda\_4 \end{aligned} \tag{A18}$$

We will prove that the discriminant, Δ, of the quadratic polynomial Equation (A17) is always positive:

$$\begin{split} \Delta = (2k)^2 - 4(1)\mathbb{C} &= \left[2(\lambda\_1 + \lambda\_2 + \lambda\_3 + \lambda\_4)/2\right]^2 - 4\left[(\lambda\_1 + \lambda\_3)(\lambda\_2 + \lambda\_4) - \left(1 - p\_f\right)(1 - p\_r)\lambda\_1\lambda\_4\right] \\ &= \left(\lambda\_1 + \lambda\_2 + \lambda\_3 + \lambda\_4\right)^2 - 4(\lambda\_1 + \lambda\_3)(\lambda\_2 + \lambda\_4) + 4\left(1 - p\_f\right)(1 - p\_r)\lambda\_1\lambda\_4 \\ &\quad \left(\lambda\_1 + \lambda\_2 + \lambda\_3 + \lambda\_4\right)^2 - 4(\lambda\_1 + \lambda\_3)(\lambda\_2 + \lambda\_4) + 4(1 - 1)(1 - p\_r)\lambda\_1\lambda\_4 \\ &= \left(\lambda\_1 + \lambda\_2 + \lambda\_3 + \lambda\_4\right)^2 - 4(\lambda\_1 + \lambda\_3)(\lambda\_2 + \lambda\_4) = \left[(\lambda\_1 + \lambda\_3) + (\lambda\_2 + \lambda\_4)\right]^2 - 4(\lambda\_1 + \lambda\_3)(\lambda\_2 + \lambda\_4) \\ &= \left(\lambda\_1 + \lambda\_3\right)^2 + \left(\lambda\_2 + \lambda\_4\right)^2 + 2(\lambda\_1 + \lambda\_3)(\lambda\_2 + \lambda\_4) - 4(\lambda\_1 + \lambda\_3)(\lambda\_2 + \lambda\_4) \\ &= \left(\lambda\_1 + \lambda\_3\right)^2 + \left(\lambda\_2 + \lambda\_4\right)^2 - 2(\lambda\_1 + \lambda\_3)(\lambda\_2 + \lambda\_4) = \left[(\lambda\_1 + \lambda\_3) - (\lambda\_2 + \lambda\_4)\right]^2 \ge 0 \\ &\Rightarrow \Delta > 0 \end{split} \tag{A1.9}$$

In Equation (A19) we used that 4 1 − *pf* (<sup>1</sup> <sup>−</sup> *pr*)*λ*1*λ*<sup>4</sup> <sup>&</sup>gt; 0 since 1 − *pf* > 0, (1 − *pr*) > 0, *λ*<sup>1</sup> > 0, and *λ*<sup>4</sup> > 0. From Equation (A19) it follows that the roots *s*<sup>1</sup> the *s2* of the quadratic polynomial Equation (A19) are always real and different:

$$s\_{1,2} = \left(-2K \pm \sqrt{\Delta}\right) / 2 = \left(-2K \pm \sqrt{4K^2 - 4C}\right) / 2 = -K \pm \sqrt{K^2 - C} \tag{A20}$$

In Equation (A20) we assume that *<sup>s</sup>*<sup>1</sup> <sup>&</sup>gt; *<sup>s</sup>*<sup>2</sup> (i.e., *<sup>s</sup>*<sup>1</sup> <sup>=</sup> <sup>−</sup>*<sup>K</sup>* <sup>+</sup> <sup>√</sup> √ *K*<sup>2</sup> − *C* and *s*<sup>2</sup> = −*K* − *K*<sup>2</sup> − *C*). It can easily be seen that the constants *<sup>s</sup>*<sup>1</sup> the *<sup>s</sup>*<sup>2</sup> are always negative. Using the quadratic factorization formula together with Equation (A17) the denominator in both Equations (A15) and (A16) can be factored to:

$$\mathbf{s}^2 + 2\mathbf{K}\mathbf{s} + \mathbf{C} = \mathbf{1}(\mathbf{s} - \mathbf{s}\_1)(\mathbf{s} - \mathbf{s}\_2) = (\mathbf{s} - \mathbf{s}\_1)(\mathbf{s} - \mathbf{s}\_2) \tag{A21}$$

From Equations (A15)–(A17), and (A21), *Y*1(*s*), *Y*2(*s*) can be simplified to:

$$Y\_1(s) = \frac{s + \lambda\_2 + \lambda\_4}{(s - s\_1)(s - s\_2)}\tag{A22}$$

$$Y\_2(s) = \frac{\left(1 - p\_f\right)\lambda\_1}{(s - s\_1)(s - s\_2)}\tag{A23}$$

Substituting Equation (A22) in the third equation of Equation (A14) we can find *Y*3(*s*):

$$Y\_3(s) = \frac{\lambda\_3 Y\_1(s)}{(\lambda\_1 + s)} = \frac{\lambda\_3 (s + \lambda\_2 + \lambda\_4)}{(s - s\_1)(s - s\_2)(\lambda\_1 + s)}\tag{A24}$$

The identified *Y*1(*s*), *Y*2(*s*), and *Y*3(*s*) are rational fractions according to Equations (A22)–(A24). To facilitate the inverse Laplace transform, those rational fractions can be subjected to a partial fraction decomposition [49] (pp. 533–540):

$$Y\_1(s) = \frac{s + \lambda\_2 + \lambda\_4}{(s - s\_1)(s - s\_2)} = \frac{A\_1}{(s - s\_1)} + \frac{B\_1}{(s - s\_2)}\tag{A25}$$

The constants *A*<sup>1</sup> and *B*<sup>1</sup> in Equation (A25) are:

$$A\_1 = \frac{s\_1 + \lambda\_2 + \lambda\_4}{s\_1 - s\_2} \text{ and } B\_1 = \frac{s\_2 + \lambda\_2 + \lambda\_4}{s\_2 - s\_1} \tag{A26}$$

$$Y\_2(s) = \frac{\left(1 - p\_f\right)\lambda\_1}{(s - s\_1)(s - s\_2)} = \frac{A\_2}{(s - s\_1)} + \frac{B\_2}{(s - s\_2)}\tag{A27}$$

The constants *A*<sup>2</sup> and *B*<sup>2</sup> in Equation (A27) are:

$$A\_2 = \frac{\left(1 - p\_f\right)\lambda\_1}{s\_1 - s\_2} \text{ and } B\_2 = \frac{\left(1 - p\_f\right)\lambda\_1}{s\_2 - s\_1} \tag{A28}$$

$$Y\_3(s) = \frac{\lambda\_3(s + \lambda\_2 + \lambda\_4)}{(s - s\_1)(s - s\_2)(\lambda\_1 + s)} = \frac{A\_3}{(s - s\_1)} + \frac{B\_3}{(s - s\_2)} + \frac{C\_3}{(\lambda\_1 + s)}\tag{A29}$$

The constants *A*3, *B*3, and *C*<sup>3</sup> in Equation (A29) are:

$$A\_3 = \frac{\lambda\_3(s\_1 + \lambda\_2 + \lambda\_4)}{(s\_1 - s\_2)(\lambda\_1 + s\_1)}, B\_3 = \frac{\lambda\_3(s\_2 + \lambda\_2 + \lambda\_4)}{(s\_2 - s\_1)(\lambda\_1 + s\_2)}, \text{ and } C\_3 = \frac{\lambda\_3(\lambda\_1 + \lambda\_2 + \lambda\_4)}{(\lambda\_1 - s\_1)(\lambda\_1 - s\_2)}\tag{A30}$$

Now, we can apply the inverse Laplace transform over Equations (A25), (A27), and (A29) and find the solutions *P*1(*t*), *P*2(*t*), and *P*3(*t*) of the stated initial-value problem:

$$\begin{cases} \text{Domain}: t \in [0, \infty) \\ P\_1(t) = A\_1 e^{s\_1 t} - B\_1 e^{s\_2 t} \\ P\_2(t) = A\_2 e^{s\_1 t} - B\_2 e^{s\_2 t} \\ P\_3(t) = A\_3 e^{s\_1 t} - B\_3 e^{s\_2 t} + C\_3 e^{-\lambda\_1 t} \end{cases} \tag{A31}$$

(b) Using the Equation (A31), the required functions can be simplified to:

$$\begin{array}{lcl} & \text{Domain :} t \in [0, \infty) \\ P\_{4}(t) & = 1 - P\_{1}(t) - P\_{2}(t) - P\_{3}(t) \\ & = 1 - \left( A\_{1}e^{s\_{1}t} - B\_{1}e^{s\_{2}t} \right) - \left( A\_{2}e^{s\_{1}t} - B\_{2}e^{s\_{2}t} \right) - \left( A\_{3}e^{s\_{1}t} - B\_{3}e^{s\_{2}t} + C\_{3}e^{-\lambda\_{1}t} \right) \\ & = 1 - A\_{1}e^{s\_{1}t} + B\_{1}e^{s\_{2}t} - A\_{2}e^{s\_{1}t} + B\_{2}e^{s\_{2}t} - A\_{3}e^{s\_{1}t} + B\_{3}e^{s\_{2}t} - C\_{3}e^{-\lambda\_{1}t} \\ & = 1 - \left( A\_{1} + A\_{2} + A\_{3} \right)e^{s\_{1}t} + \left( B\_{1} + B\_{2} + B\_{3} \right)e^{s\_{2}t} - C\_{3}e^{-\lambda\_{1}t} \end{array} \tag{A32}$$

$$\begin{array}{lcl} & \text{Domain:} \, t \in [0, \infty) \\ \, R\_{\text{sys}}(t) & = 1 - P\_4(t) \\ & = 1 - \left[ 1 - \left( A\_1 e^{s\_1 t} - B\_1 e^{s\_2 t} \right) - \left( A\_2 e^{s\_1 t} - B\_2 e^{s\_2 t} \right) - \left( A\_3 e^{s\_1 t} - B\_3 e^{s\_2 t} + C\_3 e^{-\lambda\_1 t} \right) \right] \\ & = 1 - 1 + A\_1 e^{s\_1 t} - B\_1 e^{s\_2 t} + A\_2 e^{s\_1 t} - B\_2 e^{s\_2 t} + A\_3 e^{s\_1 t} - B\_3 e^{s\_2 t} + C\_3 e^{-\lambda\_1 t} \\ & = (A\_1 + A\_2 + A\_3)e^{s\_1 t} - (B\_1 + B\_2 + B\_3)e^{s\_2 t} + C\_3 e^{-\lambda\_1 t} \end{array} \tag{A33}$$

(c) The required improper integral for *MTTFsys* when the integrand is given by Equation (A33) can be calculated using the following formula:

$$\int\_0^\infty e^{-at}dt = \frac{1}{a}\text{ where }a > 0\tag{A34}$$

Then,

$$\begin{aligned} MTTF\_{sys} &= \bigcap\_{0}^{\infty} R\_{sys}(t)dt = \bigcap\_{0}^{\infty} (A\_1 + A\_2 + A\_3)e^{\varepsilon\_1 t} - (B\_1 + B\_2 + B\_3)e^{\varepsilon\_2 t} + \mathcal{C}\_3 e^{-\lambda\_1 t}dt \\ &= (A\_1 + A\_2 + A\_3) \int\_0^{\infty} t^{s\_1 t} dt - (B\_1 + B\_2 + B\_3) \int\_0^{\varepsilon\_2 t} t^{s\_2 t} dt + \mathcal{C}\_3 \int\_0^{\varepsilon} e^{-\lambda\_1 t} dt \\ &= -(A\_1 + A\_2 + A\_3)/s\_1 + (B\_1 + B\_2 + B\_3)/s\_2 + \mathcal{C}\_3/\lambda\_1 \end{aligned} \tag{A35}$$

In the derivation shown in Equation (A35) we applied Equation (A34) three times since *s*<sup>1</sup> < 0, *s*<sup>2</sup> < 0, and (–*λ*1) < 0.

#### **References**

