Definitions

Let’s start with setting up the notation and the definitions. We will define failure probability density function (PDF), cumulative failure distribution(CDF) and the reliability function (\(R\)) as follows:

\[\begin{eqnarray} \rm{PDF}&:&f(t),\,t > 0 \label{eq:pdf}\\ \rm{CDF}&:&F(t)=\int_0^t d\tau f(\tau)\label{eq:cdf}\\ \rm{Reliability}&:&R(t)=1-F(t). \tag{1} \end{eqnarray}\] We will first prove that in both cases defined above, the reliability of a dual generator system is the same.

Run the first generator first until it dies

Define the time generator#1 dies as \(T_1=t_1\). The second generator will turn on at \(t_1\). At any time \(t\), the effective clock on generator#2 will show \(t-t_1\) (obviously with \(t> t_1\)). Therefore the probability that generator#2 will fail at \(t\) given generator#1 failed at \(t_1\) can be written as: \[\begin{eqnarray} f(T_2=t|T_1=t_1)=f(t-t_1)\Theta(t-t_1),\tag{2} \end{eqnarray}\] where \(\Theta\) is the unit step function ensuring \(t> t_1\). This is the conditional probability of generator#2 failing at \(t\) given that the first one failed at \(t_1\). To remove the condition, we need to multiply Eq. (2) with the probability of generator#1 to fail at \(t_1\) and integrate over \(t_1\):

\[\begin{eqnarray} f(T_2=t)&=&\int_0^\infty f(T_2=t|T_1=t_1)f(t_1)d t_1\nonumber\\ &=&\int_0^\infty f(t-t_1)\Theta(t-t_1) f(t_1)d t_1\nonumber\\ &=&\int_0^t f(t-t_1) f(t_1)d t_1,\tag{3} \end{eqnarray}\] which is the convolution integral of \(f(t)\) with itself. We can compute the integral once we specify \(f\), which we will do later.

Interleave generator runs

In this case, the workload is uniformly distributed over both generators until one of them fails. After the first one fails, the second one will continue alone until it also fails. This is actually a relatively simple calculation since the system failure time will be the sum of the individual generator failure times. In the context of random variables, we just need to compute the sum of two random variables, i.e. \(T=T_1+T_2\). The cumulative probability function of \(T\) is calculated as:

\[\begin{eqnarray} F_T(t)&=&P(T_1+T_2< t)\nonumber\\ &=&\int_{t_1+t_2<t}f_{T_1}(t_1)f_{T_2}(t_2) dt_1 dt_2\nonumber\\ &=&\int_{-\infty}^\infty\int_{-\infty}^{t-t_2}f_{T_2}(t_2)dt_2f_{T_1}(t_1) dt_1\nonumber\\ &=&\int_{-\infty}^\infty F_{T_2}(t-t_1)f_{T_1}(t_1) dt_1. \tag{4} \end{eqnarray}\]

The probability density function is the derivative of Eq. (4):

\[\begin{eqnarray} f_T(t)&=&\frac{d}{dt} F_T(t)=\int_{-\infty}^\infty f_{T_2}(t-t_1)f_{T_1}(t_1) dt_1\nonumber\\ &=&\int_{0}^t f(t-t_1)f(t_1) dt_1,\tag{5} \end{eqnarray}\] where the limits of the integral are truncated to the range where \(f\neq0\). Since Eqs. (3) and (5) are identical, we conclude that the system level failure distribution is identical in both cases.

Exponential distribution: \(\beta=1\)

This is a special case of Weibull distribution. \(f\) is parameterized by a single parameter \(\lambda\), which represents the failure rate: \[\begin{eqnarray} \rm{PDF}&:&f(t)=\lambda e^{-\lambda t}, t > 0\label{eq:pdfexp}\\ \rm{CDF}&:&F(t)=\int_0^t d\tau f(\tau)=1-e^{-\lambda t}\label{eq:cdfexp}\\ \rm{Reliability}&:&R(t)=1-F(t)=e^{-\lambda t}.\tag{6} \end{eqnarray}\]

From Eq. (5) we get:

\[\begin{eqnarray} f_T(t)&=&\int_{0}^t f(t-t_1)f(t_1) dt_1\nonumber\\ &=&\lambda^2 \int_{0}^t e^{-\lambda (t-t_1)}e^{-\lambda t_1} dt_1\nonumber\\ &=&\lambda^2e^{-\lambda t} \int_{0}^t dt_1=\lambda^2\, t \,e^{-\lambda t},\tag{7} \end{eqnarray}\] which is actually a \(\Gamma\) distribution (this observation will prove to be critical for generic Weibull case). The corresponding cumulative failure function is:

\[\begin{eqnarray} F_T(t)&=&\int_0^t d\tau f_T(\tau)=\lambda^2\int_{0}^td\tau\, \tau \,e^{-\lambda \tau}\nonumber\\ &=&-\lambda^2\frac{d}{d\lambda}\left[\int_{0}^td\tau\,e^{-\lambda \tau}\right]\nonumber\\ &=&\lambda^2\frac{d}{d\lambda}\left[\frac{e^{-\lambda t}-1}{\lambda}\right]\nonumber\\ &=&1-e^{-\lambda t}(1+\lambda t). \tag{8} \end{eqnarray}\]

And, finally, the reliability function reads: \[\begin{eqnarray} R(t)&=&1- F_T(t)=e^{-\lambda t}(1+\lambda t). \tag{9} \end{eqnarray}\]

Generic Weibull Distribution

Note that what made our computation super-easy was the fact that we enjoyed the homomorphism property of the exponential function: \(e^{-(t-t_1)}=e^{-t}e^{\,t_1}\), which resulted in neat cancellations. We won’t have that property for a generic Weibull distribution. This will make the convolution integral very hard to compute. Before referring to numerical computation, let’s push our luck a bit. In the case of exponential distribution for individual generators, we ended up with \(\Gamma\) distribution. Maybe, the generalized \(\Gamma\) function will be a good fit for generic Weibull case. Generalized \(\Gamma\) distribution comes with \(3\) parameters, which provides (more than) enough degrees of freedom to fit Weibull distribution. The functional form of \(\rm{gen}\Gamma\) is [1]:

\[\begin{eqnarray} \rm{gen}\Gamma(t)&=&\frac{(p/a^d)t^{d-1}e^{-(t/a)^p}}{\Gamma(d/p)}. \tag{10} \end{eqnarray}\] We are free to decide what \(p\), \(a\) and \(d\) would be so that this function reasonably represents the sum of two random variables. First consider \(p\), which appears in the exponent as the power of \(t\). From physical perspective, this represents the failure mode type. We would like to keep it same as \(\beta\), i.e. \(p=\beta\), \(\beta\) being the shape parameter of the generator Weibull distribution. As for the term \(t^{d-1}\): it is supposed to follow from an integration of two \(t^{\beta-1}\) terms multiplied together, that is \(t^{2\beta-1}\). So, let’s try \(d=2\beta\). One last parameter we have to fix is \(a\), and we can fix it in such a way that the mean value we will calculate using \(\rm{gen}\Gamma\) will be equal to the mean value of the sum of the two Weibull parameters.

The mean value of Weibull distribution, which is parameterized by scale parameter \(\alpha\) and shape parameter \(\beta\), is given by:

\[\begin{eqnarray} \mu &=&\alpha \Gamma(1+1/\beta),\tag{11} \end{eqnarray}\] which is per generator. For two generators, the mean value of the total time is

\[\begin{eqnarray} \mu_T &=&2\alpha \Gamma(1+1/\beta). \tag{12} \end{eqnarray}\]

On the \(\rm{gen}\Gamma\) side the mean value is

\[\begin{eqnarray} \mu_{\rm{gen}\Gamma} &=&a \frac{\Gamma[(d+1)/p]}{\Gamma[d/p]}=a\frac{\Gamma[2+1/\beta]}{\Gamma[2]}, \label{eq:gengmeantotal} \end{eqnarray}\] where we used \(p=\beta\) and \(d=2\beta\). We require these two to be the same, i.e., \(\mu_{\rm{gen}\Gamma} =\mu_T\). Solving this equality for \(a\), we get:

\[\begin{eqnarray} a &=&2\alpha \frac{\Gamma[1+1/\beta]}{\Gamma[2+1/\beta]} \tag{13}. \end{eqnarray}\]

So, the \(\rm{gen}\Gamma\) with the selected parameters becomes:

\[\begin{eqnarray} \rm{gen}\Gamma(t)&=&(\beta/a^{2\beta})t^{2 \beta-1}e^{-(t/a)^\beta}, \tag{14} \end{eqnarray}\] where \(a\) is defined in Eq. (13). This is the failure probability density \(f_T(t)\), and the cumulative failure probability, \(F_T(t)\) is computed as the integral, which results in \(\gamma(t)\), known as incomplete gamma function. For our purposes, we can compute the integral numerically.