*3.1. Preliminaries*

It is evident that if a *k*-out-of-*n* system's failure depends only on several its failed components, it coincides with the *k*-th order statistic from *n* iid rv *Ai* (*i* = 1, *n*) with a given cdf *<sup>A</sup>*(*t*). For simplicity, further we will denote order statistics *<sup>A</sup>*(1) ≤···≤ *<sup>A</sup>*(*k*) ≤ ··· ≤ *<sup>A</sup>*(*n*) of iid rv *Ai* (*i* = 1, *n*) by *Xi*, i.e., *Xi* = *<sup>A</sup>*(*i*) and *X*1 ≤ ··· ≤ *Xk* ≤ ··· ≤ *Xn*. Distributions of order statistics are well studied (see, for example, [30]), where it was shown that the joint probability density function (pdf) *fn*(*<sup>x</sup>*1, ... *xn*) of all order statistics *X*1 ≤ *X*2 ≤ ··· ≤ *Xn* from *n* iid rv *A*1, *A*2, ... , *An* with a given pdf *a*(*x*) has the following form:

$$f\_n(\mathbf{x}\_1, \mathbf{x}\_2, \dots, \mathbf{x}\_n) = n! a(\mathbf{x}\_1) a(\mathbf{x}\_2) \cdots a(\mathbf{x}\_n) \quad (\mathbf{x}\_1 \le \mathbf{x}\_2 \le \dots \le \mathbf{x}\_n).$$

By integration of this pdf with respect to last *n* − *k* variables one can simply find the joint pdf *fk*(*<sup>x</sup>*1, ... *xk*) of the first *k* order statistics *X*1 ≤ *X*2 ≤ ··· ≤ *Xk* from the *n* iid rv *Ai* (*i* = 1, *n*) in the domain *x*1 ≤ *x*2 ≤ ... ≤ *xk* in the form

$$f\_k(\mathbf{x}\_1, \mathbf{x}\_2, \dots, \mathbf{x}\_k) = \frac{n!}{(n-k)!} a(\mathbf{x}\_1) a(\mathbf{x}\_2) \dots a(\mathbf{x}\_k) (1 - A(\mathbf{x}\_k))^{n-k}.\tag{1}$$

However, if a failure of one of the system's components leads to the change in the residual lifetimes of all survived components, then their distributions are also changed.

#### *3.2. Transformation of Order Statistics*

Following the proposed model of the influence of components' failures on the residual lifetime of survivors, they are reduced by multiplying by some constant *ci* depending on the number of failed components. Denote by *Yi* (*i* = 1, *k*) the time of an *i*-th component failure under the conditions of increasing the load on survived components. To simplify the representation of these values in terms of order statistics *X*1 ≤ *X*2 ≤ ··· ≤ *Xn*, we introduce the following notations,

$$\mathbb{C}\_1 = (1 - c\_1), \quad \mathbb{C}\_2 = c\_1(1 - c\_2), \quad \dots, \quad \mathbb{C}\_{k-1} = c\_1 \cdots c\_{k-2}(1 - c\_{k-1}), \quad \mathbb{C}\_k = c\_1 \cdots c\_{k-1}.$$

In these notations, the following theorem holds.

**Theorem 1.** *The time to the considered system failure Yk is a linear function of order statistics of the following form:*

$$Y\_k = \mathbb{C}\_1 X\_1 + \mathbb{C}\_2 X\_2 + \dots + \mathbb{C}\_{k-1} X\_{k-1} + \mathbb{C}\_k X\_k. \tag{2}$$

**Proof.** To calculate the time of the system failure, we slightly expand the problem statement and calculate the successive time moments *Yi* (*i* = 1, *k*) of failures of the system's components under conditions of increasing load on the surviving components. To do that, we use a recursive procedure and denote by *X*(*j*) *i* the expected time moment of the *i*-th failure after the failure of the *j*-th component (*i* > *j*).

Thus, to start the induction, we have *X*(0) *i* = *Xi*. After the first failure of a component at time *Y*1 = *X*(0) 1 = *X*1 all residual lifetimes of surviving components that equal *Xi* − *X*1 for *i* > 1 decrease by a factor of *c*1, and therefore the expected failure times *X*(1) *i* for *i* > 1 take the form

$$X\_i^{(1)} = \quad X\_1^{(0)} + \mathbf{c}\_1 (X\_i^{(0)} - X\_1^{(0)}) = (1 - \mathbf{c}\_1)X\_1 + \mathbf{c}\_1 X\_i \quad i = \overline{2, k}.$$

Therefore *Y*2 = *X*(1) 2 = (1 − *<sup>c</sup>*1)*<sup>X</sup>*1 + *<sup>c</sup>*1*X*2.

Similarly, after the *j*-th failure at time *Yj* = *X*(*j*−<sup>1</sup>) *j* , the residual lifetimes *X*(*j*−<sup>1</sup>) *i* − *X*(*j*−<sup>1</sup>) *j* of all surviving components for all *i* > *j* decrease by a factor *cj*, 0 < *cj* < 1 (*j* = 1, *k*) and the expected failure times of components take the following form:

$$\begin{array}{rcl} X\_i^{(j)} &=& X\_i^{(j-1)}, \forall i \le j, \\ X\_i^{(j)} &=& X\_j^{(j-1)} + c\_j (X\_i^{(j-1)} - X\_j^{(j-1)}) = (1 - c\_j) X\_j^{(j-1)} + c\_j X\_i^{(j-1)}, \forall i > j. \end{array}$$

Thus, the expected failure times of the system components *Yj* (*j* = 1, *k*) under conditions of load redistribution equal to *Yj* = *X*(*j*−<sup>1</sup>) *j* (*j* = 1, *k*). Expressing *X*(*j*) *i* in terms of the original order statistics, we obtain the following expression for *i* > *j*:

$$\begin{split} X\_i^{(j)} &= \quad (1 - c\_1)X\_1 + c\_1(1 - c\_2)X\_2 + c\_1c\_2(1 - c\_3)X\_3 + \cdots + c\_1 \cdots c\_{j-1}(1 - c\_j)X\_j + c\_1 \cdots c\_j X\_i \\ &= \quad \sum\_{l=1}^j c\_1 \cdots c\_{l-1}(1 - c\_l)X\_l + c\_1 \cdots c\_j X\_i. \end{split} \tag{3}$$

Supposing that the last expression is true for a given *j* check it for all *i* > *j*:

$$\begin{aligned} X\_i^{(j+1)} &= (1 - c\_1)X\_1 + c\_1(1 - c\_2)X\_2 + c\_1c\_2(1 - c\_3)X\_3 + \dots + c\_1 \dots c\_{j-1}(1 - c\_j)X\_j \\ &+ \ c\_1 \cdot \dots \ c\_j(1 - c\_{j+1})X\_{j+1} + c\_1 \cdot \dots \cdot c\_j c\_{j+1} X\_i = \\ &= \sum\_{l=1}^j c\_1 \dotsm c\_{l-1}(1 - c\_l)X\_l + c\_1 \dotsm c\_j(1 - c\_{j+1})X\_{j+1} + c\_1 \dotsm c\_{j+1} X\_i = \\ &= \sum\_{l=1}^{j+1} c\_1 \dotsm c\_{l-1}(1 - c\_l)X\_l + c\_1 \dotsm c\_{j+1} X\_i. \end{aligned} \tag{4}$$

Hence, by the principle of mathematical induction, the equality (3) holds for any *j*. In terms of the original order statistics *Xi* (*i* = 1, *k*), we obtain for all *j* = 1, *k*:

$$X\_j = X\_j^{(j-1)} = (1 - c\_1)X\_1 + c\_1(1 - c\_2)X\_2 + \cdots + c\_1 \cdots c\_{j-2}(1 - c\_{j-1})X\_{j-1} + c\_1 \cdots c\_{j-1}X\_{j-1}$$

which, using the notation introduced earlier, leads to (2) for *j* = *k*, which completes the proof.

#### *3.3. Distribution of the System Failure Time*

Now move on to the calculation of the cdf *FYk* (*y*) of the system's time to failure *Yk* under the condition of redistribution of the load on the components. We will do that by taking into account expression (2) for the time of the system failure in terms of order statistics *Xi* and using Formula (1) for the joint distribution of the first *k* order statistics. Tosimplifytherepresentationofthiscdf weintroducethefollowingnotations,

$$\begin{array}{rcl}z\_0 &=& y\_\prime\\z\_i &=& z\_i(y; \mathbf{x}\_1, \dots, \mathbf{x}\_i) = \frac{y - \mathbf{C}\_1 \mathbf{x}\_1 + \mathbf{C}\_2 \mathbf{x}\_2 - \dots - \mathbf{C}\_i \mathbf{x}\_i}{\mathbf{C}\_{i+1}} \ (i = \overline{1, k-1}).\end{array} \tag{5}$$

With these notations the following theorem holds.

**Theorem 2.** *The distribution of the system's time to failure for y* > 0 *is*

$$F\_{Y\_k}(y) = \mathbf{P}\{Y\_k < y\}$$

$$= \frac{n!}{(n-k)!} \int\_0^{z\_0} a(\mathbf{x}\_1) d\mathbf{x}\_1 \int\_{\mathbf{x}\_1}^{z\_1} a(\mathbf{x}\_2) d\mathbf{x}\_2 \cdots \int\_{\mathbf{x}\_{k-1}}^{z\_{k-1}} a(\mathbf{x}\_k) (1 - A(\mathbf{x}\_k))^{n-k} d\mathbf{x}\_k. \tag{6}$$

**Proof.** According to Theorem 1 (see Formula (2)) the time *Yk* of the system failure is the linear function of the first *k* order statistics

$$Y\_k = \mathbb{C}\_1 X\_1 + \mathbb{C}\_2 X\_2 + \dots + \mathbb{C}\_{k-1} X\_{k-1} + \mathbb{C}\_k X\_k.$$

Therefore, for cdf *FYk* (*y*) of rv *Yk* in terms of pdf *fk*(*<sup>x</sup>*1, ... , *xk*) of the first *k* order statistics we obtain

$$\begin{aligned} F\_{Y\_k}(y) &= \mathbf{P}\{Y\_k < y\} \\ &= \mathbf{P}\{\mathbf{C}\_1 X\_1 + \mathbf{C}\_2 X\_2 + \dots + \mathbf{C}\_{k-1} X\_{k-1} + \mathbf{C}\_k X\_k < y\} \\ &= \int \dots \int f\_k(\mathbf{x}\_1, \mathbf{x}\_2, \dots, \mathbf{x}\_k) d\mathbf{x}\_1 \dots d\mathbf{x}\_{k\prime} \end{aligned} \tag{7}$$

where the integration domain is

$$D(\mathbf{x}\_1, \dots, \mathbf{x}\_k; y) = \{0 \le \mathbf{x}\_1 \le \dots \le \mathbf{x}\_k \colon \mathbb{C}\_1 \mathbf{x}\_1 + \mathbb{C}\_2 \mathbf{x}\_2 + \mathbb{C}\_3 \mathbf{x}\_3 + \dots + \mathbb{C}\_{k-1} \mathbf{x}\_{k-1} + \mathbb{C}\_k \mathbf{x}\_k \le y\}.$$

Let us represent this multidimensional integral as an iterated one. Taking into account that *x*1 ≤ *x*2 ≤···≤ *xk*, the integration domain can be transformed in the following way. For the last variable *xk* from the inequality

$$\mathbf{C}\_1 \mathbf{x}\_1 + \mathbf{C}\_2 \mathbf{x}\_2 + \mathbf{C}\_3 \mathbf{x}\_3 + \dots + \mathbf{C}\_{k-1} \mathbf{x}\_{k-1} + \mathbf{C}\_k \mathbf{x}\_k \le y\_\prime.$$

it follows that

$$\mathbf{x}\_{k} \le \frac{y - \mathbf{C}\_{1}\mathbf{x}\_{1} - \mathbf{C}\_{2}\mathbf{x}\_{2} - \mathbf{C}\_{3}\mathbf{x}\_{3} - \dots - \mathbf{C}\_{k-1}\mathbf{x}\_{k-1}}{\mathbf{c}\_{1} \cdots \mathbf{c}\_{k-1}} = z\_{k-1}(y; \mathbf{x}\_{1} \dots \mathbf{x}\_{k-1}).$$

Furthermore, taking into account that *xk*−1 ≤ *xk*, from the last inequality, it follows that

$$\mathbf{x}\_{k-1} \le \mathbf{x}\_k \le \frac{y - \mathbf{C}\_1 \mathbf{x}\_1 - \mathbf{C}\_2 \mathbf{x}\_2 - \mathbf{C}\_3 \mathbf{x}\_3 - \dots - \mathbf{C}\_{k-1} \mathbf{x}\_{k-1}}{\mathbf{c}\_1 \cdot \dots \cdot \mathbf{c}\_{k-1}}.$$

From this inequality with the simple algebra one can find

$$\mathbf{x}\_{k-1} \le \frac{y - \mathbf{C}\_1 \mathbf{x}\_1 - \mathbf{C}\_2 \mathbf{x}\_2 - \dots - \mathbf{C}\_{k-2} \mathbf{x}\_{k-2}}{\mathbf{c}\_1 \cdot \dots \cdot \mathbf{c}\_{k-2}} = z\_{k-2}(y; \mathbf{x}\_1, \dots, \mathbf{x}\_{k-2}).$$

Following in the same way we obtain for variable *x*2 the inequality

$$\begin{aligned} y &\geq \ (1 - c\_1)\mathbf{x}\_1 + c\_1(1 - c\_2)\mathbf{x}\_2 + c\_1 c\_2 \mathbf{x}\_3 \geq \\ &\geq \ (1 - c\_1)\mathbf{x}\_1 + c\_1(1 - c\_2)\mathbf{x}\_2 + c\_1 c\_2 \mathbf{x}\_2 = (1 - c\_1)\mathbf{x}\_1 + c\_1 \mathbf{x}\_2. \end{aligned}$$

from which it follows that

$$x\_2 \le \frac{y - (1 - c\_1)x\_1}{c\_1} \nu$$

and, at last,

$$y \ge (1 - c\_1)x\_1 + c\_1x\_1 = x\_1.$$

It means that 0 ≤ *x*1 ≤ *y*. This argumentation shows that the integration domain *<sup>D</sup>*(*<sup>x</sup>*1,..., *xk*; *y*) in terms of notations (5) can be represented as

$$D(\mathbf{x}\_1, \dots, \mathbf{x}\_k; y) = \{ \mathbf{x}\_{i-1} \le \mathbf{x}\_i \le z\_i(y; \mathbf{x}\_1, \dots, \mathbf{x}\_{i-1}) \mid (i = \overline{1, k}) \}.$$

Thus, using formula (1) for pdf *fk*(*<sup>x</sup>*1, ... , *xk*) for the first *k* order statistics and the above form of the integration domain, we can rewrite integral (7) for *y* ≥ 0 as

$$F\_{Y\_k}(y) = \frac{n!}{(n-k)!} \int\_0^y a(\mathbf{x}\_1) d\mathbf{x}\_1 \int\_{\mathbf{x}\_1}^{z\_1} a(\mathbf{x}\_2) d\mathbf{x}\_2 \cdots \int\_{\mathbf{x}\_{k-1}}^{z\_{k-1}} a(\mathbf{x}\_k) (1 - A(\mathbf{x}\_k))^{n-k} d\mathbf{x}\_k.$$

that ends the proof.

As a consequence of the theorem, the main system reliability characteristics can be calculated.

**Remark 1.** *Based on the distribution of the system's time to failure, any other system's reliability characteristics can be calculated, such as:*


#### *3.4. A Special Case: Exponential Distribution*

In a special case, when the system components' *Ai* (*i* = 1, *n*) lifetimes have exponential (*Exp*) distribution with a parameter *α* the integral (6) can be calculated analytically, but the calculations are rather cumbersome. We show it for the given value of *k* = 2. But for exponential distribution of the system components' lifetime, we propose another approach for the system lifetime distribution. It is based on the memoryless property of any exponentially distributed rv.

Denote by *Ti* the time interval between *i* −1-th and *i*-th components failures, *i* = 1, *k* − 1 (*T*0 = 0). Then due to the memoryless property of the exponential distribution the time to the *k*-th failure *Yk* is the sum

$$\chi\_k = T\_1 + T\_2 + \dots + T\_{k\nu}$$

of *k* independent exponentially distributed rv *Ti* with parameters

$$
\lambda\_1 = n\mathfrak{a}, \quad \lambda\_i = c\_1 c\_2 \cdots c\_{i-1} (n-i+1)\mathfrak{a} = \mathfrak{e}\_i (n-i+1)\mathfrak{a}, \; i = 2, k, \; i
$$

where for simplicity additional notations are used:

$$
\bar{c}\_{i} = \begin{cases} 1, & i = 1, \\ c\_{1} \cdot \cdots \cdot c\_{i-1}, & i = \overline{2, k}. \end{cases}
$$

The moment generating function (mgf) of the system's lifetime in this case has the following form:

$$\phi\_k(s) = \mathbb{E}\left[e^{-sY\_k}\right] = \prod\_{1 \le i \le k} \mathbb{E}\left[e^{-sT\_i}\right] = \prod\_{1 \le i \le k} \frac{\mathbb{E}\_i \lambda\_i}{s + \mathbb{E}\_i \lambda\_i}.$$

.

To apply the above theorem and the proposed approach, let us consider the simplest example of a *k*-out-of-*n* model with *k* = 2. In this case, suppose *c*1 = *c*. Thus, according to (1) the joint distribution of rv *<sup>X</sup>*(1), *<sup>X</sup>*(2) is

$$f\_2(\mathbf{x}\_1, \mathbf{x}\_2) = \frac{n!}{(n-2)!} a(\mathbf{x}\_1) a(\mathbf{x}\_2) (1 - A(\mathbf{x}\_2))^{n-2} = \frac{n!}{(n-2)!} a^2 e^{-a\mathbf{x}\_1} e^{-(n-1)a\mathbf{x}\_2}.$$

Calculate cdf *FY*2 (*y*) of the rv *Y*2 = (1 − *<sup>c</sup>*)*<sup>X</sup>*(1) + *cX*(2),

$$\begin{split} F\_{Y\_2}(y) &=& \mathbf{P} \{ (1-c)\mathbf{X}\_{(1)} + c\mathbf{X}\_{(2)} < y \} = \mathbf{P} \left\{ \mathbf{X}\_{(2)} < \frac{y - (1-c)\mathbf{X}\_{(1)}}{c} \right\} \\ &=& n(n-1)a^2 \int\_0^y e^{-ax\_1} dx\_1 \int\_{\mathbf{x}\_1}^c e^{-(n-1)ax\_2} dx\_2 \\ &=& 1 + \frac{n-1}{nc - (n-1)} e^{-ny} - \frac{nc}{nc - (n-1)} e^{-\frac{(n-1)a}{c}y} \end{split}$$

and therefore its pdf for *y* ≥ 0 is

$$f\_{Y\_2}(y) = \frac{n(n-1)\alpha}{nc - (n-1)} \left( e^{-\frac{(n-1)\alpha}{c}y} - e^{-\alpha xy} \right).$$

Please note that this result holds for *c* = (*n* − <sup>1</sup>)/*n* and in this case the distribution is a mixture of exponential distributions. The point *c* = (*n* − <sup>1</sup>)/*n* is a singular point for which cdf of the rv *Y*2 is the Erlang distribution,

$$F\_{Y\_2}(y) = 1 - e^{-n\alpha y} - n\alpha y e^{-n\alpha y}, \ y > 0,$$

with pdf

$$p\_{Y\_2}(y) = n^2 \lambda^2 y e^{-n\lambda y}, \ y > 0.$$

**Remark 2.** *The singularity in the calculation of the cdf of the system's lifetime arises because for some special values of the coefficient ci (here for c* = (*n* − 1)/*n) the moment generating function of the system's lifetime has multiple roots that leads to changing of the shape of distribution.*

With the help of another approach one can find mgf of the system's lifetime in the following form:

$$\phi\_2(s) = \frac{n(n-1)a^2}{s^2 + (2n-1)a s + n(n-1)a^2}$$

.

By expanding this expression into simple fractions, we find

$$\phi\_2(s) = \frac{n(n-1)\alpha}{s+n\alpha} - \frac{n(n-1)\alpha}{s+(n-1)\alpha}s$$

then, by calculating the inverse function, we obtain

$$f\_2(y) = n(n-1)a\left(e^{-(n-1)ay} - e^{-nay}\right),$$

which is the same as the result above for *c* = 1.

The analytical calculations of the reliability characteristics are not always possible. Nevertheless, their numerical analysis in the wide domain of initial data is possible. Therefore, in the next section a procedure for the numerical calculation of different reliability characteristics of the considered system will be proposed. Furthermore, in Section 5 this procedure will be used for the numerical analysis of the model with some examples.

#### **4. The General Calculation Procedure of the System Reliability Characteristics and Numerical Experiments**

Based on the results of the previous section, the general procedure for the problem solution can be implemented with the help of the following algorithm (Algorithm 1).

#### **Algorithm 1** : General algorithm for calculation of reliability function

**Beginning.** Determine: Integers *n*, *k*, real *ci* (*i* = 1, *k*), distribution *A*(*t*) of the system components' lifetime and its pdf.

**Step 1.** Taking into account that the system's failure moment *Yk* according to formula (2) equals

$$Y\_k = \mathbb{C}\_1 X\_1 + \mathbb{C}\_2 X\_2 + \dots + \mathbb{C}\_{k-1} X\_{k-1} + \mathbb{C}\_k X\_{k-1}$$

calculate the following,

$$\mathbb{C}\_{l} = \begin{cases} 1 - c\_{l\prime} & i = 1, \\ c\_{1} \cdots c\_{l-1} (1 - c\_{l}) & i = \overline{2, k - 1}, \\ c\_{1} \cdots c\_{k - 1} & i = k. \end{cases}$$

**Step 2.** Taking into account that according to formula (1), the joint distribution density of first *k* order statistics *X*1 ≤ *X*2 ≤···≤ *Xk* holds

$$f\_{X\_1 X\_2 \dots X\_k}(\mathbf{x}\_1, \mathbf{x}\_2, \dots, \mathbf{x}\_k) = \frac{n!}{(n-k)!} a(\mathbf{x}\_1) a(\mathbf{x}\_2) \dots a(\mathbf{x}\_k) (1 - A(\mathbf{x}\_k))^{n-k},$$

with which following to (6) calculate the reliability function

$$R(y) = 1 - F\_{\overline{\chi}\_{\mathbb{K}}}(y) = 1 - \frac{n!}{(n-k)!} \int\_0^y a(\chi\_1) d\chi\_1 \int\_{\chi\_1}^{z\_1} a(\chi\_2) d\chi\_2 \cdots \int\_{\chi\_{k-1}}^{z\_{k-1}} a(\chi\_k) (1 - A(\chi\_k))^{n-k} d\chi\_k$$

where the limits of integration are determined by the relation (5)

$$z\_0 = y, \quad z\_i = z\_i(y; \mathbf{x}\_1, \dots, \mathbf{x}\_i) = \frac{y - \mathbf{C}\_1 \mathbf{x}\_1 + \mathbf{C}\_2 \mathbf{x}\_2 - \dots - \mathbf{C}\_i \mathbf{x}\_i}{\mathbf{c}\_1 \mathbf{c}\_2 \dots \mathbf{c}\_i} \text{ ( $i = \overline{1, k-1}$ )}.$$

Find the values of the constants *ci* (singular points at which the denominator of the cdf *FYk* (*y*) turns into 0) for which the cdf changes its appearance.

**Step 3.** From the system reliability function *<sup>R</sup>*(*y*), calculate

– mean time to the system failure

$$\mu\_T = \mathbb{E}[\mathbb{Y}\_k] = \int\_0^\infty \mathbb{R}(y) dy;$$

– its variance

$$
\sigma\_T^2 = \mathbf{Var}[Y\_k] = \int\_0^\infty (y - \mu\_T)^2 f(y) dy, \quad \text{where} \quad f(y) = \frac{d}{dy} F\_{Y\_k}(y),
$$

*σ*

*μ*

.

*v* =

and coefficient of variation

**Stop.**

**Remark 3.** *The algorithm can also be used to solve other different problems, for example, to analyze the sensitivity of the system's reliability function and its characteristics to the shape of the lifetime distribution of the system's components.*

Furthermore, the algorithm will be applied to some examples.

#### **5. Numerical Experiments: 2-Out-Of-6 System**

According to Algorithm 1, we calculate the reliability function of a 2-out-of-6 system. Since such a system fails due to the failure of two components, we have only one constant that defines the decreasing residual lifetime of surviving components. Therefore, hereafter, we suppose *c*1 = *c*. Consider the Gnedenko–Weibull (*GW*) distribution as the lifetime distribution of the system's components, *A*(*t*) ∼ *GW<sup>θ</sup>*, *a* Γ(1 + *<sup>θ</sup>*−<sup>1</sup>), with the corresponding cdf

$$A(t) = 1 - \exp\left\{-\left(\frac{t\Gamma(1+\theta^{-1})}{a}\right)^{\theta}\right\}, t > 0\_{\prime\prime}$$

where


$$\bullet \qquad v = \frac{\sigma}{a} = a^{-1} \cdot \sqrt{\frac{\Gamma(1+2\cdot\theta^{-1})}{\Gamma(1+\theta^{-1})^2} - 1} \text{ is the coefficient of variation, } \frac{1}{a}$$

• *σ* is the standard deviation.

> Additionally, consider the Erlang (*Erl*) distribution, *A*(*t*) ∼ *Erl*(*l*, *θ*) with pdf

$$a(y) = \frac{\theta^l}{\Gamma(l)} y^{l-1} e^{-\theta y}, \ y > 0.$$

In this case, the distribution's parameters can be represented via the corresponding mean *a* and coefficient of variation *v* as follows,

$$l = v^{-2}, \quad \theta = (av^2)^{-1}.$$

For numerical experiments, we consider the reliability function and its characteristics of a 2-out-of-6 system for given distributions with a fixed mean *a* and different values of *v*. Thus, we can analyze the influence of the coefficient of variation of the repair time on the reliability characteristics of the system. In other words, investigate its sensitivity.

Suppose that the mean lifetime of the component *a* = 1. If *θ* = 1, *GW* and *Erl* distributions transform into the exponential one with the mean time *a* and the coefficient of variation *v* = 1. In this case, its reliability function is

$$R(t) = \frac{5e^{-6t} - 6c \cdot e^{-\frac{5t}{c}}}{5 - 6c}. \tag{8}$$

From Formula (8) it is clear that *c* = 56 leads to changing of the shape of distribution. Since calculating the coefficient *θ* for *GW* through the value of *v* is quite difficult, we define the parameter *θ* so that *v* ≈ 0.5. Moreover, if *θ* of *GW* takes non-integer values, it is not always possible to obtain a closed-form reliability function *R*(*t*) according to Algorithm 1 (the integrand takes a complex form). Therefore, define *θ* = 2, then, the coefficient of variation *v* = 0.5227. For *Erl* distribution, suppose that *v* = 0.5, which leads to *θ* = 4.

Suppose *c* = 0.1; 0.5; 1. Figure 3 illustrates the reliability function of the 2-out-of-6 system for different distributions, as well as *c* and *v*. Here, solid line means *v* = 1 and reliability function (8), dashed one is for *GW* with *v* = 0.5227 and dash-dotted is for *Erl* with *v* = 0.5. The legend of the figure denotes the color of line for different *c*.

**Figure 3.** Reliability function of a 2-out-of-6 system.

The figure shows that the higher reliability coincides with the lower value of *v*. The case *c* = 1 means the absence of load from the failed components to the surviving ones, thus this case corresponds to the highest reliability for different *v* compared to the values *c* < 1. Moreover, dependence of the reliability function curve on the shape of lifetime distribution is observed. On a small interval *y*, the system reliability as *A* ∼ *Erl* is higher than as *A* ∼ *GW* for each *c*, despite close values *v*. This may indicate the sensitivity of the reliability function not only to the shape of the lifetime distribution, but also to the corresponding value of the coefficient of variation.

According to the algorithm, we calculate other reliability characteristics of the 2-outof-6 system (Tables 1 and 2). These characteristics correspond to the system's reliability behavior, shown in Figure 3. The lower value of *v* leads to the higher value of the system lifetime expectation **<sup>E</sup>**[*<sup>Y</sup>*2], and the lower value of *c* leads to the lower value of **<sup>E</sup>**[*<sup>Y</sup>*2]. Moreover, as *v* ≈ 0.5 the relative error between the considered distributions is 14.11% for *c* = 0.1, 7.98% for *c* = 0.5 and 3.86% for *c* = 1.

**Table 1. <sup>E</sup>**[*<sup>Y</sup>*2] of a 2-out-of-6 system.


To distinguish coefficients of variation of the components and the whole system, denote them as *vcomp* and *vsys*, respectively. Thus, Table 2 shows the following. With a decrease in *c*, the coefficient of variation of the system *vsys* grows and tends to the value of the coefficient of variation of each system component *vcomp*. The increasing *vcomp* leads to the increasing *vsys* for all distributions and *c*. Thus, the coefficient of variation of the system *vsys* confirms that as *c* tends to 0 and *vcomp* tends to 1, variability with respect to the average lifetime of the system **<sup>E</sup>**[*<sup>Y</sup>*2] grows.

**Table 2.** *vsys* of a 2-out-of-6 system.


Table 1 showed that with increasing *c* and *v* ≈ 0.5, the mean system lifetime **<sup>E</sup>**[*<sup>Y</sup>*2] is very close. However, Figure 3 shows that over a small interval *y* with these *c* and *v*, the reliability of the system has significant differences. This leads to the study the quantiles of the system reliability. This measure shows how long the system will be reliable with a fixed probability. The quantiles *qγ* = *<sup>R</sup>*−<sup>1</sup>(*γ*) of the reliability function are presented in Figures 4–6. In all cases, red bullets correspond to *γ* = 0.99, whereas black bullets correspond to *γ* = 0.9.

All the values for quantiles *γ* = 0.999; 0.99; 0.9 are presented in Table 3 for different distributions. The values in the table show that for the presented quantiles *qγ*, the shape of the lifetime distribution of the system's components as well as its coefficient of variation play a critical role on the system's reliability. Therefore, for example, as *c* = 0.1 and *A* ∼ *Erl* a given reliability level 0.9 will last about 8 times longer than for *c* = 0.1 and *A* ∼ *Exp*. At that for *q*0.999, the difference for similar case is almost 40 times. As the coefficient *c* increases, this difference decreases for all values of the quantiles and lifetime distributions of the components. As *c* = 1 this difference is reduced by about two times. Thus, even as *c* = 1, which defines no changing in components' residual lifetimes, the influence not only of the lifetime distribution of the components but also its coefficient of variation on the reliability of the system is huge. This once again confirms the sensitivity of the reliability characteristics of the *k*-out-of-*n* system to the shape of the lifetime distribution and the coefficient of variation of system's components.

**Figure 4.** Reliability function with *v* = 1 and quantiles (*A* ∼ *Exp*). Red and black bullets are the points of intersection of the reliability function curves with fixed reliability levels of 0.99 and 0.9, respectively.

**Figure 5.** Reliability function with *v* = 0.5227 and quantiles (*A* ∼ *GW*). Red and black bullets are the points of intersection of the reliability function curves with fixed reliability levels of 0.99 and 0.9, respectively.

**Figure 6.** Reliability function with *v* = 0.5 and quantiles (*A* ∼ *Erl*). Red and black bullets are the points of intersection of the reliability function curves with fixed reliability levels of 0.99 and 0.9, respectively.

**Table 3.** Quantiles of reliability function *qγ*.


#### **6. Conclusions and the Further Investigations**

The reliability function of a new *k*-out-of-*n* : *F* model is investigated, under the new assumptions that the failures of its components lead to the increase in the load on the remaining ones and, consequently, to the change in their residual lifetimes. To model the situation, we proposed a novel approach based on the transformation of the order statistics of the system components' lifetimes , which is the new field of application of order statistics. An algorithm for calculation of the system's reliability function and its moments has been developed. Numerical experiments for the special case of the considered model based on the real-world systems have been carried out. The experiments show an essential sensitivity of the model reliability function and its moments to the shapes of the lifetime distributions of the system's components and their coefficient of variation.

Furthermore, it is proposed we extend this approach to the investigation of stationary characteristics of the model and consider its preventive maintenance, aiming to improve its reliability characteristics.

**Author Contributions:** Conceptualization, writing–original draft preparation, supervision, project administration, V.R.; validation, investigation, visualization, N.I.; writing–review and editing, methodology, D.K.; data curation, software, T.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This paper has been supported by the RUDN University Strategic Academic Leadership Program (recipients V.R, supervision and problem setting, N.I., visualization, D.K. writing–review and editing, T.M., analytic results). This paper has been partially funded by RFBR according to the research projects No.20-01-00575A (recipients V.R., conceptualization, and N.I., formal analysis) and RSF according to the research projects No.22-49-02023 (recipient N.I., validation, D.K. review and analytic results).

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

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** The authors express their gratitude to the Referees for the valuable suggestions, which improved the quality of the paper.

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