ANALYZING LOAD SHARING SYSTEM RELIABILITY: A MODIFIED WEIBULL DISTRIBUTION APPROACH Текст научной статьи по специальности «Электротехника, электронная техника, информационные технологии»

ANALYZING LOAD SHARING SYSTEM RELIABILITY: A MODIFIED WEIBULL DISTRIBUTION APPROACH

Santosh S. Sutar

Yashwantrao Chavan School of Rural Development, Shivaji University, Kolhapur - 416 004, India,

sutarssantosh@gmail.com •

Chandrakant G. Gardi*

Department of Statistics, Punyashlok Ahilyadevi Holkar Solapur University,

Solapur, 413255, India, chandrakant.gardi@gmail.com

Somnath D. Pawar

Department of Statistics, Shivaji University, Kolhapur - 416 004, India. sdpawar23@gmail.com

* Corresponding author Abstract

Load sharing systems have the ability to distribute the workload among its components. For analyzing a two component parallel load sharing system, the accelerated failure time (AFT) based model with component lifetimes as linear failure rate distribution have been recently proposed in the literature. In the present study, the component lifetimes are assumed to follow a modified Weibull distribution, which is the generalization of many standard lifetime distributions such as exponential, Weibull, Rayleigh, and linear failure rate. The use of modified Weibull distribution leads to a new family ofbivariate distributions for ordered random variables. We have also looked into the associated inference techniques for the proposed model. In order to evaluate the effectiveness of the suggested estimating approaches, we conducted a simulation study. In order to provide a practical application and better understanding, we carefully examine a dataset related to motors.

Keywords: accelerated failure time model, conditional distribution, load sharing, modified Weibull distribution, order statistics

1. Introduction

Load sharing systems are characterized by their ability to distribute the workload among multiple components, such that if one component fails, the remaining components bear the additional workload. This can either increase or decrease the load on each surviving component. Load sharing systems have been extensively investigated in various engineering domains, such as software and hardware reliability, power plants, computing workload analysis, and fiber composites (Wang et al. [1]).

Liu [2] presents various instances that illustrate the concept of load sharing systems. These include scenarios like electric generators distributing an electrical load within a power plant, CPUs operating in a multiprocessor computer system, cables supporting a suspension bridge, bolts

fastening a wheel assembly onto a truck, and valves or pumps functioning in a hydraulic system. When any of these components fail, the remaining components must bear an additional load, which can elevate their failure rates. In an intriguing study by Drummond et al. [3] conducted on vertebrate species, it was observed that when a litter mate dies due to food shortage, the surviving offsprings receive a larger portion of the available food supply, leading to improved growth. This finding highlights how the failure of one individual can inadvertently benefit the surviving members. Furthermore, in the realm of software testing, the detection of a fault can uncover previously undetected critical faults. This demonstrates that a component's failure can facilitate the discovery of other hidden issues, thereby enhancing the overall reliability of the system. These examples collectively illustrate that when a component fails, it can actually enhance the remaining components' remaining lifespan, resulting in a higher growth rate for the surviving components.

Daniels [4] conducted the first study on the phenomenon of load sharing and load sharing systems. A thorough analysis of load sharing systems up till 2009 is present in Dewan and Naik-Nimbalkar [5]. Deshpande et al. [6], Park [7], Singh and Gupta [8], Park [9], Gurov and Utkin [10], Sutar and Naik-Nimbalkar [11]-[12], Krivtsov et al. [13], Wang et al. [1] and Sutar and Naik-Nimbalkar [14] have all published studies and modeled the load sharing phenomenon since then.

The study of load-sharing systems with a k-out-of-m configuration was suggested by Sutar and Naik-Nimbalkar [11] with modelling strategy based on the accelerated failure time (AFT) model. They concentrated on a particular configuration, a parallel load sharing system consisting of two components with baseline as the linear failure rate distribution. The associated inference techniques were also covered by the researchers. The distributions used there in for the ordered random variables are a subset of a larger family of distributions known as sequential order statistics. Kamps [15] first described this family of distributions and Cramer and Kamps [16] further developed them.

This study utilizes the load sharing model based on accelerated failure time (AFT), proposed by Sutar and Naik-Nimbalkar [11], to examine the load sharing phenomenon within a parallel system consisting of two components. We adopt a modified Weibull distribution (MWD) as the baseline distribution for the components in the system, which is characterized by three parameters: Ai, A2, and A3. The introduction of this three-parameter MWD was done by Sarhan and Zaindin [17].

It is important to highlight that the three-parameter MWD provides a comprehensive representation of various distributions, including exponential, Weibull, Rayleigh and linear failure rate. Thus, the MWD serves as a versatile baseline distribution for the component lifetime in any load sharing system. The subsequent sections of this paper are organized as follows.

In Section 2, we address the AFT-based load-sharing model for a parallel load sharing system consisting of two components and with a modified Weibull distribution for component lifetimes. In Section 3, the inference procedures are thoroughly examined, while Section 4 focuses on the simulation study. Section 5 demonstrates an application using real data, and the last section presents the concluding remarks.

2. Proposed AFT based Load sharing model

We investigate a parallel system consisting of two components. The cumulative distribution function (c.d.f.) of the components follows a MWD characterized by three parameters: A1, A2, and A3. The probability density function (p.d.f.), survival function (s.f.), and hazard rate function of the MWD with parameters A1, A2, and A3 are provided as follows.

f (u) = A + A2A3uA3-1) exp { - (\1 u + A2ma^|, u > 0, A1, A2, A3 > 0, A1 + A2 > 0 (1)

F(u) = exp {- (Aiu + A2uAs) j ,u > 0,A1,A2,A3 > 0,Ai + A2 > 0 h(u) = A + A2A3uA3-1) , u > 0, A1, A2, A3 > 0, A1 + A2 > 0

The three-parameter modified Weibull distribution, denoted as MWD(A1, A2, A3), serves as a generalization of the following distributions.

(a) It represents the Exponential distribution, ED(A1), when A2 is set to 0 and A3 is finite.

(b) It encompasses the Weibull distribution, WD(A2, A3), when A1 is set to 0.

(c) It corresponds to the Rayleigh distribution, RD(A2), when A3 is set to 2 and A1 is set to 0.

(d) It encompasses the Linear failure rate, LFR(A1, A2), when A3 is set to 2.

For more details on MWD(A1, A2, A3), one can refer Sarhan and Zaindin [17] and references cited therein.

The load sharing behavior observed in a system comprising two components is captured by the AFT model, which was introduced by Sutar and Naik-Nimbalkar [11]. We denote the lifetimes of the two components in the system as V1 and V2. These lifetimes are considered independent and identically distributed random variables. The baseline densities of V1 and V2 are denoted as f1 (■) and f2(■), respectively, while their corresponding baseline survival functions are denoted as F1 (■) and F2(■). Let X = min(V1, V2) denote time of the first failure and Y = max(V;L, V2) denote the time of the second failure or the system failure time. Consequently, the marginal density of the first failure can be expressed as follows.

g(x) = (2A1 + 2A2A3xA3-1) exp { - (2A1 x + 2A2xA3) }, x > 0, A1 > 0, A2, A3 > 0, A1 + A2 > 0.

(2)

It is worth mentioning that the distribution of the first failure is identical to the baseline distribution, which is a modified Weibull distribution with parameters (2A1,2A2, A3). Following the AFT load sharing model, the conditional density of variable Y given that X = x, as well as the joint density of the variables X and Y, can be expressed in the following manner.

g(y|x) = { $ + ^y^ } exp { - A (y " x) " £ (yA3 " xA3)} , (3)

0 < x < y < to, $ > 0, A1 > 0, A2 > 0, A3 > 0, A1 + A2 > 0,

g(x, y) = 2 (A1 + A2A3xA3-1){ $ +

x exp { - $ (y - x) - (yA - xA3) - (2A1 x + 2A2xA3 ^ , (4)

0 < x < y < to, $ > 0, A1 > 0, A2 > 0, A3 > 0, A1 + A2 > 0, It is important to note that when the parameter $ is equal to 1, the joint density described in equation (4) simplifies to the joint density of independent random variables X and Y. Essentially, when $ = 1, it indicates no load sharing effect, and hence the occurrences of the two failures (first and second) are independent of each other. The parameter $ is referred to as the load sharing parameter in this context. In the following section, we will delve into the inference procedures associated with this concept.

3. Inference procedures

In this section, we examine different methods for estimating the unknown parameters and introduce a testing procedure to assess the presence of the load sharing effect.

3.1. Direct estimation procedure

The complete data, denoted as (x, y) = {(xi, yi) : xi < yi; i = 1,2,..., n}, represents the set of observations. The log-likelihood function based on this complete data can be expressed as follows.

n / \ n , \

log L = n log 2 + £ log (Ai + A2A3xf3-1) + £ log (AipA3 + A2A3pyf3 -1) - n(A3 + 1) log p

i=1 i=1 A n A n n n

- A^T £ (yi - Xi) - A3 £ (yf3 - xf3) - 2A1 £ Xi - 2A2 £ xf3. p i=1 P3 i=1 i=1 i=1

We observe that, the log-likelihood equations d L = 0, ddAgL = 0, L = 0 and L = 0 do not have explicit solutions for the parameters A1, A2, A3, p. The mathematical expressions for the score functions, specifically d dpgL = 0, d^g-L = 0, ^AgL = 0, and L = 0, can be found in Appendix (A).

In the following subsection, we outline a two-step approach for determining the values of the unknown parameters A1, A2, A3, and p.

3.2. Two-step parameter estimation procedure

The process of estimating the values of A1, A2, A3, and p has been conducted using a two-step methodology.

Step 1. We observe the first failure, X, and estimate baseline parameters, namely, A1, A2 and A3 by using the MCEM procedure proposed by Sutar [18].

Step 2. In order to estimate the load sharing parameter p, we utilize the conditional distribution of Y given X = x, as expressed in equation (3). The estimates of A1, A2, and A3 obtained in Step 1 are then substituted into that equation to perform the estimation. We refer to this estimation process as a two-step estimation procedure, and the subsequent subsections outline these two steps in detail.

3.2.1 Estimation of A1, A2 and A3 (Step 1)

It is worth noting that the distribution of the first failure, X, is also a modified Weibull distribution (MWD) with parameters 2A1, 2A2, and A3. Let us denote 2A1 = 71, 2A2 = 72, thus distribution of X is MWD with parameters 71, y2 and A3. We use the MCEM algorithm proposed by Sutar [18] for finding the estimates of 71, y2 and A3. To implement the proposed MCEM algorithm, we take two independent random variables U1 and U2, which has exponential (71) and Weibull (72, A3) distributions, with their respective survival functions as exp(-71M1) and exp(-y2uf3). Let 71, 72 and A3 be the MLEs of 71, 72 and A3 obtained through MCEM algorithm, then the MLEs of A1, A2 and A3 can be obtained as A1 = Yr, A2 = Y2 and A3.

3.2.2 Estimation of p (Step 2)

To estimate the load sharing parameter p, we utilize the conditional distribution of Y given X = x as described in equation (3). In this study, two methods are proposed for estimating p, which are discussed as follows.

Method I : It can be noted, the conditional distribution of Y given X = x as truncated MWD with parameters p = 01 (say), p^ = d2 (say), A3 = 93 (say) truncated at X = x. Furthermore, the conditional distribution mentioned is equivalent to the distribution of the minimum value between

two independent random variables, denoted as Wi and W2. Specifically, Wi follows a truncated exponential distribution, which is truncated below x and has a parameter 01. On the other hand, W2 follows a truncated Weibull distribution, also truncated below x, with parameters 02 and 03. The survival functions of W1 and W2 are exp{-01(w1 — x)} and exp{-02(w03 — x03)}, respectively. Let us consider the complete data for i = 1,2,..., n as a set of 2n independent random variables denoted as (W1i, W2i). The random variables W1i represent truncated exponential distributions, truncated below xi, with a parameter 01, while W2i represent truncated Weibull distributions, truncated at xi, with parameters (02, 03). Additionally, we define Z2i as the minimum value between W1i and W2i. Consequently, Z2i follows a truncated MWD (Minimum of Weibull and Exponential Distribution) distribution, characterized by a probability density function (p.d.f)

g(z2i) = (01 + 0203z23—1) exp {—01 (z2i — xi) — 02(z23 — x03)} ,

0 < xi < z2i < TO, fi > 0,01 > 0,02 > 0,03 > 0,01 + 02 > 0. We can regard the observed values y = (y1, y2,..., yn) as corresponding to the values of Z2 =

(Z2^ Z22,..., Z2n).

The joint density of W1 and W2 given X = x(= (x1, x2,..., xn)) can be written as

n { }

g(w1, w2|x) = {01,02,03}n n w03—1ex^ — 01(w1i — xi) — 02(w03 — x03 )} , (5)

i=1

The log-likelihood can be expressed in the following manner.

n n n

log L = n log (010203) + 03 E log(u2i) — 01 E uu — 02 E ^3.

i=1 i=1 i=1

In order to perform the E step, it is necessary to calculate the conditional expectation of Ec [log L|Z2]. This can be represented as follows.

Ec [log L|Z2] = n log (0*0203) + 03Ec

E log(U2i)|Z2

— 0i Ec

E Uu |Z2

i=1

- 07* E,

i=1

0,

E U3IZ2

i=1

(6)

Remark 1. For i equal to 1, 2, and 3, the variables 0i and 0* represent the values of 0i at the current iteration and the next iteration of the MCEM (Monte Carlo Expectation-Maximization)

( p)

algorithm. Specifically, if 0i = 0y represents the estimated value of 0i at the p-th iteration,

and 0* = 0(p+1) represents the estimated value of 0i at the (p + 1)-th iteration, then 0i and 0* respectively denote the values of 0i at the p-th and (p + 1)-th iterations of the MCEM algorithm.

The conditional density of W11 given X = x and Z21 = z21 is a mixed probability density function (p.d.f.) and can be expressed as follows.

0 { }

g(w1|x, z21) = (-1 0 —1){I(w1=z21)+ 0203z231—1 exp { —01(w1 — z21)} Wz^)} .

( 01 + 0203z21 )

where, Ia(■) is indicator function defined on set A. The details regarding the same are Appendix (B). Thus, the conditional expectation of W11 given X and Z21 can be obtained as

E(W1|X, Z21) = J w1g(w1 |x,z21)dw1

0 f f TO --1-) |Z21 + 0203Z23—1 exp {01Z21} Jz w1 exp {—01w1} dw^ .

01 + 0203z21

03—1'

n

By using the result,

£ w1exp{-01 wi}dwi = Z21 + 1) g? Z21},

we get

E(Wi|X,Z2i) = 7-(z2I + 0203Z2l-iexp{glZ21 }(fli Z21 + 1)Z21} }

(01 + 0203z|-1) I fl1 J

1 + { Z21 - (01 + 0203z23-1 ) ^ = K(Z21 )(say),

and hence

(n \ n n n I

E Wii ix, Z2 = n + E Z2i - E (-

i=i J 0i i=i i=i (01 + 0203Z2I-1) Thus, given {Z2i = yi, i = 1,2,..., n} and {Xi = xi, i = 1,2,..., n}, we get

in \ n n n I

E E Wii|x,y = n + Eyi - E 1

Vi=1 V 01 i=1 i=1 (01 + 0203y03-1

Likewise, the conditional density function of W21 given X = x, Z21 = z21, and the corresponding conditional expectation can be expressed as follows.

g(w2|x,z21) = 0203 -1\ {z231-11(w2=z21) + 01w23-1 e{-02(W203 -z21)}I(w2>z21) }

(^01 + 0203z21 J ^ '

and

Ec {wg|x Z21} = + (^fl^L) r wfl +*-1e-

1 ' («1 + 0203Z2i(ft + 0203Z22'i-^ ^

After simplification, we get,

E {W3| XZ } = °2°3 Z1+9 3+1 + W " ^ r (u + 0 Z3) ^ e-udu

E^W2iIX/ ^ = (01 +0203 Z2i"1) + (01 + 0203 Z2i"1) J0 ^ + ^J 6 du

Let

fli

V + 2Z

fl3

T = j (u + 02Z23) "3 e-udu = E

where V has an exponential distribution with mean 1. By employing the Monte Carlo technique, we can replace T with a Monte Carlo sum, which is given as follows.

T = E

, fli I m . fli

V + 02= mm E (v, + 02Z23l) fl3

where, {v1, v2,..., vm } is random sample of size m (sufficiently large) from exponential distribution with mean 1. Thus, we get

fli

763 + 03-1 a a-%

W . |XZ = 6263Z23 I ~ E (v + e2ZO^

i=1 "2i ^^ j i=i j ^fli + fl2fl3z|3-^ m (01 + 8283Z2?-1) 7=1 V 7 2i ;

9

and hence {Z2i = yi, i = 1,2,..., n}, that is Z2 = y, we get

62%:

0* + 03 — 1

03—y E (vj + 02 y0

Ec\ EWin*,y = E i ( ' "•" 0 1 y + D11_

\i=1 -J i=1 1 (01 + M^03—M m 01 + 02%p— 1 j=1

By applying similar arguments to calculate Ec ^En=1 log (W2i) |x,y^, we obtain

eJ E log (W2i) lx, y = E

0203 log yi

i=1

1 \ (01 + 0203y03—^

+ E

1

3-

1t (£ fa + 02y03) — lOg0^.

1 [ m02 (01 + 0203y

To carry out the M-step, we obtain the following expression from equation (6).

dEc [log L|X, Z2 ] = n_E d0* 0* c

dEc [logL|X,Z2] _ n

- — — — Er

d02 02 c

E W1i |X, Z

i=1

E WnX, Z2

i=1

dEc [log L|X, Z2] = n_E

d0* 0* c

E log(W2i )|X, Z2

i=1

2*

dEc [En=1 W2||X,Z2]

d0 *

From (7) and (8), we get

dEc [logL|X,Z2 d 0*

dEc [log L|X, Z] d0 *

0 ^ 0*

0 ^ 02

Ec [En=1 W1i |X, Z2]' n

Ec

En=1 W2i3|X, Z2

(7)

(8) (9)

(10) (11)

Based on the expression (5), it can be inferred that t(W1) = En=1 W1i serves as the sufficient statistic for 01. In order to carry out the M-step, we equate the sufficient statistic to its expectation, which takes the following form in our scenario.

E[W1 ] = 0

(12)

(p)

The EM iterations alternate between expressions (12) and (10). Let 1 represent the estimate of 01 at the p-th iteration step of the MCEM algorithm. The updated estimate 0(p+1) is determined by the following equation.

01 p+1) = ^ + -

1 > 0(p) n

E (yi—xi) — E

i=1

i=1 (0( p)+02 p) 03 p) y03p)—1

(13)

To estimate 03, we substitute equation (11) into equation (9), resulting in a nonlinear equation in terms of 3. This equation can be expressed as follows.

dEc

dEc

log L|x y

log W2,3 |x,y

d0 *

03* Ec

E (log(W2i)) |x,y

i=1

B0i

n

En=1 W2? |x,y

n

n

n

1

1

n

(r)

We employ the Newton-Raphson iterative method to estimate 0|. Let 03 denote the estimate of 03

(r+1)

at the r-th iteration of the Newton-Raphson method. The updated estimate 3 is determined by the following equation. { }

dEc [log L|x,y]

i(r+1)

-

dt!')

d2 Ec [log L|x,y]

d( e<r))2

(14)

The detailed expressions used in equation (14) can be found in Appendix (C). The process should be iterated until the convergence criterion is satisfied. The value of 03 at the (p + 1)-th iteration

of MCEM, denoted as 03

0.

(p+1)

, represents the stabilized value of 3 obtained through the

Newton-Raphson method. Once we get the estimate of 03, we can obtain 02 at (p + 1)th iteration by substituting (14) in (11). That is 02 p+1) is given by

,(p+1)

n

En

'¿=1

0(p+i) W2? y

(15)

Thus, we get the estimates 0 = (A1/p), 02 = (A2/pA3) of (A1/p) and (A2/pA3) respectively. Now, we have the two estimates of p as p01 = (A1/01) and p2 = (A2/02)1/A3.Here, A1, A2, and A3 represent the estimates of A1, A2, and A3, respectively, obtained in Step 1. The estimate of p is obtained by taking the average of p1 and p2. This estimation method is referred to as the 'average method' for estimating p.

Method II : The likelihood, based on the conditional density of Y given X = x, can be expressed as follows.

L = n ,1 + A2M

L=ü\ ß + ßX3

A3-1'

A1

A2

exH E (yi - Xi) — ^A23 E (yA3 - xA3)

ß

¿=1

Then the log-likelihood can be written as

1 L f, A1 + A2A3yi

log L=SSI

A3-1

IA3 ¿=1

A2

- A1 E (yi - Xi) E (yf3 - xf3)

I

i=1

Then we have,

d

A1 + A2A3 y,3~

n I ß^ ßA3+1

A1

IA3 ,=1

A2 A3

dß log L iSi A^ I A2A3yA3-M 1 ß2i=Ti ^ IA3+1 i=1'

rl + £E (yi - Xi) -SE (yA3 - xA3)

a1 I

i +

IA3

d2 n

^L=- fi

2 A3 - 1

A1 + A2 a2 y.-3

I2 ^ ßA3+1

2A1 + A2 A3(A3 + 1)y-3 1

I3 + ßA3+2

+ , ,1 + A2A3 y-3 1

+ ^ ß2 + ßA3+1

,1 + A2A3yi 3 I + ßA3

-^E (yi - xi) - A2ApAA+2+1) E (yA3 - xA3).

To estimate the load sharing parameter p, we substitute the maximum likelihood estimates (MLEs) A1, A2, and A3 of A1, A2, and A3, respectively, into the aforementioned expressions. The Newton-Raphson iteration method is then employed to obtain the estimate of p. Let p(m) be the estimate of p at mth iteration. The estimate of p at (m + 1)th iteration is given by

ß(m+1) = ß(m) -

dß log L) l(A1,A2,

A3)

di21^ L)1 (-1,-2,-3)

3

2

n

n

n

2

We termed this procedure to estimate fi as the 'iteration method'. The subsequent section presents the test procedure used to assess the presence of the load sharing effect.

3.3. Test Procedure

To test the load sharing effect, we set up the null hypothesis H0 stating that the failure of a

component does not affect the survival components, and the alternative hypothesis H1 stating that there exists a load sharing phenomenon. Specifically, we express the null hypothesis as H0 : fi = 1, indicating no load sharing effect, and the alternative hypothesis as H1 : fi = 1, indicating the presence of a load sharing effect. To test these hypotheses, we employ a score type test, as used by Sutar and Naik-Nimbalkar [11]. The test statistic follows an asymptotic distribution with 1 degree of freedom.

In the subsequent section, we present a simulation study to evaluate and compare the performance of two estimation methods: the average method and the iterative method. The simulation study aims to assess the accuracy and efficiency of these methods under various scenarios and conditions. By conducting simulations and analyzing the results, we can gain insights into the strengths and limitations of each method and make informed decisions about their suitability for practical applications.

4. Simulation study

In this section, we performed a simulation study to assess the performance of the proposed estimation procedure in estimating unknown parameters. We generated a total of 10,000 samples from the joint density described in equation (4) for different combinations of sample sizes (n) and parameter values. This allowed us to examine the behavior and accuracy of the estimation procedure under various scenarios and conditions.

For sample sizes of n = 20, 30, 50, and 100, we considered different parameter combinations, namely (A1, A2, A3, fi) as (1,2,0.5,0.5), (1,2,0.5,1), (1,2,0.5,1.5), (1,2,1,0.5), (1,2,1,1), (1,2,1,1.5), (1,2,2,0.5), (1,2,2,1), (1,2,2,1.5), (2,2,0.5,0.5), (2,2,0.5,1), (2,2,0.5,1.5), (2,2,1,0.5), (2,2,1,1), (2,2,1,1.5), (2,2,2,0.5), (2,2,2,1), and (2,2,2,1.5).

The average estimates of the unknown parameters (A1, A2, A3, fi) obtained through Method-I (Two-step Procedure), denoted as (A1, A2, A3, fi), along with their corresponding standard errors (SE), i.e., SE(A1), SE(A2), SE(A3), and SE(fi), are presented in Table 1 and Table 2. The simulation results reveal a clear pattern: as the sample size grows larger, the standard errors exhibit a decreasing trend. This indicates that larger sample sizes lead to enhanced precision in estimating the parameters, implying that the estimates become more accurate and reliable.

We also conducted a simulation study for the iterative method. We generated 10,000 samples with sizes n = 30,50, and 100 from the joint density given in equation (4). We considered different parameter combinations as (1,2,0.5,0.5), (1,2,0.5,1), (1,2,0.5,1.5), (1,2,1,0.5), (1,2,1,1), (1,2,1,1.5), (1,2,2,0.5), (1,2,2,1), (1,2,2,1.5).

The estimates of the unknown parameters (A1, A2, A3, fi), where the estimate of fi obtained through the Method-II (iterative method), along with their corresponding standard errors (SEA), SE(A2), SE(A3), and SE(fi)), are reported in Table 3. We observed that compared to the estimates obtained by the average method, the estimates obtained by the iterative method had higher standard errors and tended to be overestimated.

Same phenomenon is observed for the simulation for study corresponding to the parameter combination (A1, A2, A3, fi) as (2,2,0.5,0.5), (2,2,0.5,1), (2,2,0.5,1.5), (2,2,1,0.5), (2,2,1,1), (2,2,1,1.5), (2,2,2,0.5), (2,2,2,1) and (2,2,2,1.5). Hence, we decided not to report the simulation results corresponding to these parameter combination.

Table 1: The average estimates of (A1, Ai, A3, fi) obtained through the two-step procedure.

n A1 A2 A3 fi A1 A2 A3 J SE(A1) SE(A2) SE(A13) SE(/J)

30 1 2 0.5 0.5 1.6347 3.2042 1.4232 0.7157 0.0798 0.0854 0.0693 0.0516

50 1 2 0.5 0.5 1.5741 2.5862 1.0873 0.6217 0.0759 0.0871 0.0669 0.0463

100 1 2 0.5 0.5 1.2432 2.3124 0.7554 0.5482 0.0748 0.0854 0.0675 0.0454

30 1 2 0.5 1 1.6865 3.256 1.8661 1.4245 0.0981 0.0847 0.0775 0.0611

50 1 2 0.5 1 1.4885 2.5789 1.0832 1.3029 0.0868 0.1001 0.0798 0.0579

100 1 2 0.5 1 1.1941 2.2879 0.7286 1.1229 0.0782 0.0861 0.0692 0.0461

30 1 2 0.5 1.5 1.8624 3.2677 1.9431 1.8289 0.1031 0.0962 0.0885 0.0721

50 1 2 0.5 1.5 1.5765 2.4989 0.9867 1.7110 0.0887 0.1042 0.0875 0.0598

100 1 2 0.5 1.5 1.2299 2.3093 0.7589 1.5603 0.0798 0.0954 0.0781 0.0476

30 1 2 1 0.5 1.8921 3.2776 2.4102 0.7372 0.0986 0.0865 0.0703 0.0627

50 1 2 1 0.5 1.5132 2.5867 1.5305 0.6134 0.0764 0.0967 0.0686 0.0511

100 1 2 1 0.5 1.2389 2.3123 1.2397 0.5623 0.0831 0.0876 0.0779 0.0467

30 1 2 1 1 1.7868 3.2682 2.3405 1.4321 0.1031 0.0881 0.0832 0.0684

50 1 2 1 1 1.5105 2.6105 1.5193 1.3139 0.0872 0.0989 0.0794 0.0673

100 1 2 1 1 1.1962 2.2961 1.2204 1.1283 0.0864 0.0872 0.0689 0.0551

30 1 2 1 1.5 1.8611 3.2692 2.3952 1.8382 0.1084 0.1051 0.0967 0.0685

50 1 2 1 1.5 1.5902 2.5156 1.4734 1.7102 0.0972 0.1098 0.0935 0.0637

100 1 2 1 1.5 1.2346 2.3203 1.2087 1.5589 0.0876 0.0971 0.0798 0.0472

30 1 2 2 0.5 1.9614 3.2658 3.4231 0.7267 0.0889 0.0847 0.0704 0.0542

50 1 2 2 0.5 1.5837 2.5991 2.5237 0.6079 0.0769 0.0885 0.0679 0.0476

100 1 2 2 0.5 1.2437 2.2984 2.2389 0.5472 0.0768 0.0853 0.0672 0.0462

30 1 2 2 1 1.8773 3.2674 3.4212 1.4298 0.0974 0.0869 0.0773 0.0616

50 1 2 2 1 1.4932 2.5813 2.5326 1.2998 0.0867 0.0983 0.0795 0.0568

100 1 2 2 1 1.2167 2.2916 2.2193 1.1183 0.0792 0.0851 0.0693 0.0463

30 1 2 2 1.5 1.8823 3.2593 3.4261 1.8672 0.1006 0.0975 0.0869 0.0578

50 1 2 2 1.5 1.5824 2.5139 2.4934 1.6672 0.0891 0.1092 0.0879 0.0573

100 1 2 2 1.5 1.2305 2.3027 2.1979 1.5723 0.07967 0.0945 0.0797 0.0493

Table 2: The average estimates of (A1, A2, A3, fi) obtained through the two-step procedure.

n Ai A2 A3 ß Ai a2 A3 ß SE(A1) SE(A2) SE(A13) SE(ß)

30 2 2 0.5 0.5 2.9587 3.2681 1.9672 0.7361 0.0893 0.0842 0.0669 0.0578

50 2 2 0.5 0.5 2.5831 2.5951 1.0851 0.6138 0.0765 0.0879 0.0677 0.0456

100 2 2 0.5 0.5 2.2360 2.3013 0.7542 0.5479 0.0754 0.0851 0.0679 0.0442

30 2 2 0.5 1 2.8476 3.2821 1.8567 1.4310 0.0981 0.0843 0.0774 0.0621

50 2 2 0.5 1 2.4881 2.5813 1.0829 1.3021 0.0870 0.0995 0.0792 0.0582

100 2 2 0.5 1 2.1933 2.2929 0.7300 1.1232 0.0798 0.0856 0.0686 0.0450

30 2 2 0.5 1.5 2.8621 3.2713 1.9413 1.8348 0.1116 0.0961 0.0873 0.0635

50 2 2 0.5 1.5 2.5855 2.5018 0.9964 1.7027 0.0877 0.1030 0.0862 0.0604

100 2 2 0.5 1.5 2.2320 2.3058 0.7542 1.5590 0.0802 0.0934 0.0770 0.0482

30 2 2 1 0.5 2.9745 3.2799 2.3941 0.7416 0.0971 0.0845 0.0689 0.0618

50 2 2 1 0.5 2.5941 2.5979 1.5238 0.6156 0.0771 0.0962 0.0690 0.0504

100 2 2 1 0.5 2.2468 2.3015 1.2416 0.5601 0.0815 0.0860 0.0758 0.0485

30 2 2 1 1 2.8891 3.2801 2.3407 1.4392 0.1028 0.0871 0.0817 0.0671

50 2 2 1 1 2.4932 2.5918 1.5262 1.3124 0.0881 0.0998 0.0811 0.0656

100 2 2 1 1 2.2003 2.3056 1.2185 1.1273 0.0841 0.0887 0.0694 0.0529

30 2 2 1 1.5 2.8601 3.2713 2.4006 1.8425 0.1061 0.1027 0.0946 0.0684

50 2 2 1 1.5 2.5888 2.5139 1.4866 1.7073 0.0926 0.1115 0.0915 0.0630

100 2 2 1 1.5 2.2387 2.3117 1.2032 1.5622 0.0857 0.0965 0.0831 0.0490

30 2 2 2 0.5 2.9601 3.2589 3.4098 0.7244 0.0862 0.0818 0.0671 0.0512

50 2 2 2 0.5 2.5785 2.5853 2.5164 0.6021 0.0741 0.0861 0.0664 0.0447

100 2 2 2 0.5 2.2351 2.2937 2.2336 0.5423 0.0727 0.0815 0.0639 0.0413

30 2 2 2 1 2.8744 3.2701 3.4189 1.4251 0.0948 0.0841 0.0751 0.0591

50 2 2 2 1 2.4821 2.5784 2.5261 1.2982 0.0836 0.0971 0.0783 0.0555

100 2 2 2 1 2.1927 2.2891 2.2164 1.1153 0.0774 0.0816 0.0683 0.0436

30 2 2 2 1.5 2.8513 3.2587 3.4228 1.8294 0.0987 0.0939 0.0849 0.0611

50 2 2 2 1.5 2.5751 2.4992 2.4839 1.6982 0.0838 0.1018 0.0839 0.0572

100 2 2 2 1.5 2.2194 2.2943 2.1926 1.5468 0.0782 0.0896 0.0744 0.0451

Table 3: The average estimates o/(Aj, A2, A3, p) obtained through the iterative method.

n A1 A2 A3 ß A1 A2 A3 SE(Ai) se(A2) se(A13) SE(ß)

30 1 2 0.5 0.5 1.9887 3.3081 1.9873 0.8154 0.0893 0.0902 0.0678 0.0603

50 1 2 0.5 0.5 1.6912 2.6332 1.0923 0.7385 0.0782 0.0893 0.0691 0.0472

100 1 2 0.5 0.5 1.3497 2.3213 0.7743 0.6293 0.0772 0.0869 0.0695 0.0449

30 1 2 0.5 1 1.9534 3.3109 1.9576 1.4734 0.1023 0.0953 0.0867 0.0703

50 1 2 0.5 1 1.5934 2.6156 1.0921 1.3921 0.0899 0.1012 0.0823 0.0612

100 1 2 0.5 1 1.3173 2.4679 0.7591 1.1581 0.0823 0.0897 0.0728 0.0499

30 1 2 0.5 1.5 1.9727 3.3278 1.9825 1.9568 0.1792 0.1034 0.0957 0.0684

50 1 2 0.5 1.5 1.7455 2.5357 0.9764 1.8627 0.0893 0.1125 0.0897 0.0692

100 1 2 0.5 1.5 1.5671 2.3513 0.7934 1.6193 0.0934 0.1045 0.0842 0.0502

30 1 2 1 0.5 1.9834 3.2895 2.4231 0.8245 0.1034 0.0957 0.0725 0.0769

50 1 2 1 0.5 1.6761 2.6281 1.5482 0.7372 0.0821 0.0993 0.0723 0.0584

100 1 2 1 0.5 1.3182 2.3756 1.2949 0.6429 0.0902 0.0931 0.0784 0.0521

30 1 2 1 1 1.9233 3.2956 2.3756 1.5334 0.1342 0.0913 0.0882 0.0705

50 1 2 1 1 1.5421 2.6492 1.5849 1.4294 0.0917 0.1034 0.0942 0.0736

100 1 2 1 1 1.2951 2.3682 1.2735 1.2661 0.0879 0.0921 0.0704 0.0569

30 1 2 1 1.5 1.9349 3.2937 2.4623 1.9173 0.1236 0.1412 0.1034 0.0756

50 1 2 1 1.5 1.5923 2.5634 1.4954 1.8728 0.1054 0.1532 0.1031 0.0713

100 1 2 1 1.5 1.3681 2.3542 1.2348 1.6212 0.0886 0.0993 0.0902 0.0534

30 1 2 2 0.5 1.9789 3.2782 3.4267 0.8178 0.0921 0.0941 0.0714 0.0589

50 1 2 2 0.5 1.5845 2.5936 2.5372 0.6912 0.0797 0.0901 0.0686 0.0484

100 1 2 2 0.5 1.3383 2.3417 2.2756 0.6389 0.0810 0.0889 0.0725 0.0467

30 1 2 2 1 1.8972 3.2852 3.4462 1.5343 0.1034 0.0872 0.0792 0.0602

50 1 2 2 1 1.4939 2.5789 2.5319 1.4014 0.0913 0.0989 0.0810 0.0579

100 1 2 2 1 1.3214 2.2973 2.2682 1.2314 0.0824 0.0901 0.0713 0.0498

30 1 2 2 1.5 1.9344 3.2604 3.4610 1.9282 0.1083 0.0991 0.0892 0.0659

50 1 2 2 1.5 1.6952 2.5021 2.4934 1.8317 0.0884 0.1153 0.0897 0.0604

100 1 2 2 1.5 1.3156 2.3023 2.2103 1.7116 0.0816 0.0927 0.0829 0.0523

5. Illustration

In this section, we have applied the AFT based load sharing model and estimation procedures to motor data obtained from Reliability Edge Home [19]. The dataset consists of 18 systems, each consisting of two motors operating continuously in parallel. The failure times of both motors, along with their identification labels A and B, were recorded.

Our objectives were twofold. Firstly, we aimed to determine whether the modified Weibull distribution (MWD) is an appropriate baseline distribution for modeling the lifetimes of both components. Secondly, we aimed to test whether there exists a load sharing phenomenon, where the failure of one motor affects the working of the other.

To assess the appropriateness of the MWD as the baseline distribution, we conducted a

Kolmogorov-Smirnov type test, which confirmed its suitability. However, it should be noted that the test was conservative due to the estimation of unknown parameters. We utilized a two-step estimation procedure, with the estimation of ft being conducted using the 'average method'. The estimated values of A1, A2, A3, and ft were found to be 0.0028, 2.08168 x 1016, 31.6118, and 1.9847, respectively.

To investigate the presence of load sharing among the motor failure times, we employed a score-type test proposed by Sutar and Naik-Nimbalkar [11]. The computed test statistic value was 19.564, which surpassed the critical values at both the 1% and 5% significance levels. Consequently, we can infer that the failure of one motor has a significant impact on the lifetime of the other. This finding supports the existence of a load sharing phenomenon, where the failures of individual components influence the performance of the remaining components in the system. This conclusion is further supported by the estimated value of ft being 1.9847 (significatly different than 1), suggesting that these 18 parallel systems exhibit load sharing or a load sharing phenomenon among the component failures.

6. Conclusions

In our study, we focused on a two-component parallel load sharing system and utilized the accelerated failure time based load sharing model to capture the load sharing behavior observed in this system. We chose the modified Weibull distribution as the baseline distribution for the component lifetime. We proposed two procedures for estimating the model parameters within this framework and also discussed a test procedure for assessing the presence of load sharing in such systems. Furthermore, we conducted a simulation study to evaluate the performance of the proposed estimation procedures, which demonstrated satisfactory results. To illustrate the practical applicability of the load sharing system, we analyzed a specific dataset. It is worth mentioning that the modeling and analysis of load sharing phenomena can be extended to more complex systems, such as a k-out-of-m system.

Acknowledgment

The Santosh S. Sutar is grateful for the financial support provided by Department of Science and Technology (DST), Science and Engineering Research Board (SERB), Government of India under a Core Research Grant (File No. CRG/2021/005672).

References

[1] Wang, D., Jiang, C. and Park, C. (2019). Reliability analysis of load-sharing systems with memory. Lifetime Data Analysis, 25:341-360.

[2] Liu, H. (1998). Reliability of a Load-Sharing k-out-of-n:G System: Non-iid Components with Arbitrary Distributions. I.E.E.E. Transactions on Reliability, 47:279-284.

[3] Drummond, H., Vazquez, E., Sanchez-Colon, S., Martinez-Gomez, M., Hudson, R. (2000). Competition for milk in the domestic rabbit: survivors benefit from littermate deaths. Ethology, 106:511-526.

[4] Daniels, H. E. (1945). The statistical theory of the strength of bundles of threads. I. Proceedings of the Royal Society London A, 183:404-435.

[5] Dewan I. and Naik-Nimbalkar U. V. (2010). Load-Sharing Systems. Wiley Encyclopedia of Operations, Research and Management Science, https://doi.org/10.1002/9780470400531.eorms0476.

[6] Deshpande, J. V., Dewan I. and Naik-Nimbalkar, U. V. (2010). A family of distributions to model load sharing systems. J. Statist. Plann. Infer., 140:1441-1451.

[7] Park, C. (2010). Parameter estimation for the reliability of load-sharing systems. IIE Trans, 42:753-765.

[8] Singh B. and Gupta, P. K. (2012). Load-sharing system model and its application to the real data set. Mathematics and Computers in Simulation 82:1615-1629.

Park, C. (2013). Parameter estimation from load-sharing system data using the expectation maximization algorithm. HE Trans, 45:147-163.

Gurov S. V. and Utkin, L. V. (2014). A load-share reliability model under the changeable piecewise smooth load. /onrnal o/Qnality and Reliability Engineering, 1-11. Sutar S. S. and Naik-Nimbalkar, U. V. (2014). Accelerated Failure Time Models For Load Sharing Systems. I.E.E.E. Transactions on Reliability, 63:706-714.

Sutar S. S. and Naik-Nimbalkar, U. V. (2016). A model for k-out-of-m load sharing systems. Commnnications in Statistics - Theory and Methods, 45(20):5946-5960.

Krivtsov, V., Amari S. and Gurevich V. (2017). Load sharing in series configuration. Quality and Reliability Engineering International, 1-12.

Sutar S. S. and Naik-Nimbalkar, U. V. (2019). A load share model for non-identical components of a k-out-of-m system. Applied Mathematical Modelling, 72:486-498. Kamps, U. (1995). A concept of generalized order statistics. /. Statist. Plann. Infer. 48:1-23. Cramer, E. and Kamps, U. (2003). Marginal distributions of sequential and generalized order statistics. Metrika, 58:293-310.

Sarhan A. M. and Zaindin M. (2009). Modified Weibull distribution. Applied Science, 11:123136.

Sutar S. S. (2014). Parameter estimation of the modified Weibull distribution using Monte Carlo Expectation Maximization algorithm. Model Assisted Statistics and Applications, 11:171-178.

Reliability Edge Home, (2003), Quarter 2, 3, (http://www.reliasoftindia.com/newsletter/ 4q2002/qalt.htm).

and

Appendix (A): Expressions involved in Score functions

dlogL _ " (A1A3ßA3 1 + A2a3yA3 ^ _ n(A3 +1) Ai \h( _ N dß k (A1 ßA3 + A2 A3ßyA3-1) ß + ß2 HVl Xl)

A2A3 A3 A3n n

ßA+T L (Vr - x3 ) =

i=1

d log L = » ßA3-1 + » 1

dA1 ¿=1 (A1 ßA3 + A2 A3 ßyA3-1) i=1 (A1 + A2 A3 xA3-1)

1 n n

- ß L(yi- x) - 2L xi = 0 ß i=1 i=1

d log L " ßA3vA3-1 + f A3 xf3-1

L A ,1 1 „„.A3-1) + L '

dA2 i=1 (A1 ßA3 + A2 A3 ßyA3-1) i=1 (A1 + A2 A3 xA3-1)

A n n

- A3 L(yA - xA3) - 2A2 L xA3 = 0,

ß 3 i=1 i=1

d log L " (A1 ßA3 logß + A2ßyA3-1 + A2A3ßyA3-1 logyi) dA3 = i=1 (A1 ßA3 + A2 A3 ßyA3-1)

+ £ A2^ + A2A3) - n _ A2 £ /1( ß <

i=1 (A1 + A2A3 xA3-1J i=A ßJ vß,

- A2 t (j) A3 log (I) - 2A2 ExA3 log(xi) = 0.

(16)

Appendix (B): The information pertaining to the conditional density of

W11 given X = x and Z21 = z2i

Consider

G(wL,z21 |X = x) = P(Z21 > z21, W11 > w1 |X = x) = f?w1|x=x (max(Z2l, wI))fw2 |

X=x (Z21)

= exp j -0l [max(z21, wL) - x] - 02 (z2l - x03) j .

Due to ordering in z21 and wL, we have following three cases-

1. z21 > wL i.e. max(z21, wL ) = z21.

2. z21 < wL i.e. max(z21, wL ) = wL.

3. z21 = wL i.e. max(z21, wL) = z21 or wL. When, z21 > wL we have,

G(wL, z21 |X = x) = exp { -0l (z21 - x) - 02(z2i - x03 )} .

Thus, we have

d2 -

g(wL,z21|X = x) = -—-— G(z21, wL|x) = 0, z21 > wL > 0. dz2idWi

When, z21 < wL we have,

G(wL,z21 |X = x) = exp j -01(w1 - x) - 02(z2l - x03) j , wL > z21 > 0,

and hence

32

g(wi,z2i|X = x) = G(z2i, wi|x) , Wi > Z21 > 0.

That is

g(wL, z211X = x) = 0l0203z2l-1 exp { -0l (wL - x) - 02 (z2l - x03 )} , wL > z21 > 0. When, z21 = wL we have,

G(wL,z21|X = x) = exp{-0l(z21 - x) - 02(z| - x03)} , wL = z21 > 0. Thus, we have

d2

g(wi,z2i|X = x) = G(z2i, wi|x) , wi = Z21 > 0

= 0Lexp j -0l(z21 - x) - 02(z| - x03)} j ,z21 = wL , wL = z21 > 0.

Thus, by combining all the above cases, we can write the joint density as

g(wi,z2i|X = x) = 010203z2l-1 exp{-01 (wi - x) - 02(z| - x03)}I(z2i = Wi)

+01 exp{-01 (z2i - x) - 02(z2i - x03)}I(z2i < Wi ). Hence, the conditional density of WL1 given X = x, Z21 = z21 can be obtained as g(wi, z2i |X = x)

g(wi|x, Z21)

g(Z2i|X = X)

-1-—T {7(wi=Z2i) + 020 3z23i-1exp {-0l(wl - Z21)} i(wi>Z2i)}

Ql + 0203z2i-

Appendix (C): Expressions involved in (14)

Ec

d2 Ec

with

log y

■ g(r)

W2i3 |x, y

d(g 3r))2

?3r))2

n

Ec

(r)

W2? ^y

( g(r+1)

3E^Ln=1 W3 |x, y

dg,

(r)

E

m £ 3r+1)+%-1

log(yi)

1 (g1 + g2g3yg3-1)

and

(r+1)

2r) yg3) g3r)

g{r) g2r) - ^ (vj + g2r) y\3) i=1 h mg? (g 1+g2g3yg3-1)

+

i=1 j=1

e3r+1) e3r+1)

g(r) g2r) - g3r) (vj+g2r) #) g3r) log fa+g2r) y^)

mg3r) (g1 + g2g3yf -1

i g(r+1) 1

d E^Ln=1 |x,y|

d( g 3r))2

E

i=1

g 2g3j/ g3r+1)+g3-1 (log)yi ))2

(g 1 + g2g3yg3-1

g^

>)>) g(r)

g^

nm w ^ (logg2r^2 (vj+g2r)g3r)

+ ¿L iL

i=1 j=1

(r+1)

m

3-1

^r)a(r) g(r)

n m 1 2 3

-EE

i=1 j=1

g^

(log(g«)) (Vj + g2r)yf) g3r) log (Vj + g2r)yf)

m

(r+1)

(r)

3-1

(r+1)

-EE

i=1 j=1

g(r)g2r) g3r) (vj+g2r)y?) g3r) log ^+g«yg3)

m(g3r)) (g1+g2g3y

3-1

(r+1)

(r+1)

+

i=1 j=1

g(r)g2r) g3r) (v+g2r)y?) g3r) (log ^+g2r)y?)y

m

(r)

3-1

n

2

3

3

nm

3

2