Complex electromechanical systems (CEMSs) have been used in industry, civil machinery, aerospace, and so on (Mi et al., 2016; Wang et al., 2022). Reliability synthesis and prediction of CEMSs have long been critical issues within the field of reliability engineering. CEMSs are usually composed of optics, machinery, and electricity (Tang et al., 2021). Some CEMSs have very complex structures, some have huge volumes, some have long service life, and some are expensive (Han et al., 2021). Thus, reliability synthesis and prediction of CEMSs only relying on life testing are limited by multiple factors such as test scheme, test site, test time, test cost, and test conditions. However, at the stage of the project demonstration and initial design of CEMSs, reliability synthesis and prediction of CEMSs are necessary (Fuqiu et al., 2020). Therefore, it challenges the realization of reliability synthesis and prediction of CEMSs.

Since Buehler proposed a reliability synthesis method for a series system with two binomial components (Buehler, 1957), scholars from classical, Bayes, and fiducial schools have conducted a lot of research on reliability synthesis and prediction of systems. The classical school put forward the maximum likelihood estimation (MLE) method, modified maximum likelihood estimation (MML) method, sequential reduction (SR) method, combined MML and SR method (CMSR), and so on (Yu et al., 2013). The MLE method is only applicable to large sample tests, and the failure distribution is an unbounded symmetric normal distribution. The MML method is only applicable to binomial (pass–fail) component systems (Fo and Xiong, 2009). The disadvantage of the SR method is that successive compression of test data leads to information loss, resulting in large estimation variance and a more conservative confidence limit of system reliability (Zhu and J Sh, 1990). The algorithm of the CMSR method is complex, and it is not suitable for the system with zero failure number of unit test data. The Bayes method only solves the reliability interval estimation of binomial or exponential approximation for complex systems (Guo and Wilson, 2013). The fiducial method has been proved unsuitable for system reliability evaluation (Zhou and Weng, 1990). Therefore, it is difficult to use the abovementioned methods to synthesize and predict the CEMSs’ reliability, and reliability synthesis and prediction of CEMSs are still problems in reliability engineering (Peng et al., 2013). Although the Monte Carlo method can easily deal with the reliability synthesis problem of systems with different unit distributions (Yeh et al., 2010), the simple Monte Carlo method, bootstrap method, double Monte Carlo method, and asymptotic method introduced in many studies can only solve the lower confidence limit of system reliability, and the interval estimation problem of system life distribution types and parameters are unsolved (Meng and Wen-Tao, 2021). Despite decades of research, there is still no general method to realize the reliability synthesis and prediction of CEMSs (Kovacs et al., 2019).

It is well known that a CEMS is composed of several units (subsystems), and the system reliability depends on the reliability of its constituent units (subsystems) (Negi and Singh, 2015), (Liu et al., 2015). When the constituent units’ life of CEMS follows different distributions, it is difficult to synthesize and predict the CEMS reliability from units’ reliability data information (Weiyan et al., 2009). Hence, some scholars have conducted a lot of correlational research. For example, Guo et al. (2014) analyzed some theories and methods for system reliability synthesis, and the new idea was proposed for the problem aiming at specific engineering backgrounds. Considering that some CEMSs contain outsourced components, Sun et al. (2018) computed system reliability and component importance measures. According to the equal-principle first and second moment of reliability, Yu et al. (2013) investigated Bayes reliability confidence limit for a series-parallel system consisting of different distribution units. Graves and Hamada (2016) evaluated the likelihood for simultaneous failure time data when monitoring was stopped. In view of the lack of reliability data information of CEMSs, Wilson et al. (2006) and Yuan et al. (2019) discussed how to make full use of limited reliability data information of CMESs to synthesize and predict system reliability. Although scholars have discussed the reliability synthesis and prediction of CEMS from many aspects, there are still challenges by many factors such as difficulty in establishing a reliability model for CMES, lack of reliability data information, high test cost, and so on. It is necessary to be studied further.

In order to solve the above problems and realize the reliability synthesis and prediction of CEMS, a numerical analysis method based on the computer simulation and the Monte Carlo (MC) method is proposed in this article. The remainder of this article is organized as follows: In Section 2, we propose a method to obtain the life simulation values for the constituent units of CEMS when the life distribution types and distribution parameters of constituent units are known. In Section 3, we take a method to obtain the system life simulation value from the units’ life simulation values. One system life simulation value represents a reliability test of CEMS. Computer simulation can be realized instead of life test when there are enough system life simulation values. In Section 4, we carry out the initial selection of CMES life distribution types and goodness-of-fit test. In Section 5, a case study of a CMES is presented to demonstrate the effectiveness of our proposed method. Finally, conclusions are made in Section 6.

It is assumed that a CEMS is composed of n different life distribution units, and the life distribution and the appropriate parameters of each constituent unit are known. The life simulation values for n units are generated using the random variable simulation conversion equation and through Monte Carlo simulation (Yoshida and Akiyama, 2011; Wang et al., 2012; Zhang et al., 2014). According to the common distributions in engineering, such as the exponential distribution, Weibull distribution, normal distribution, logarithmic normal distribution, extreme value distribution, Gamma distribution etc., the respective random variable simulation conversion equations g ( ϑ i ) are deduced and listed in Table 1. If a unit is subject to any other distribution, its random variable simulation conversion equation can be solved similarly using the inverse function method.

In order to ensure that all units are independent of each other, the pseudo-random number is computer-generated for each unit and subject to 0-1uniform distribution. ϑ i = R N D ( 1 )       i = 1,2 , ⋯ , n .

According to the life distribution and parameters of each constituent unit, the life simulation value t i for the corresponding unit can be obtained through logical operation. t i = g ( ϑ i )       i = 1,2 , ⋯ , n .

In this way, a set of life simulation values of n units can be obtained, and the first simulation values are completed. The m groups of life simulation values are obtained after m cycles.

For the CEMS composed of n units, the life simulation value for each unit can be obtained through each simulation, and they are denoted by t j 1 , t j 2 , ⋯ , t j n . According to the reliability logical relationship between the CEMS and each component unit, the life simulation values for the CMES are obtained on the minimal path set method (Cancela et al., 2013; Schallert, 2014). S = ∪ h = 1 q S h = ∪ h = 1 q [ ∩ x i , x j , ⋯ ∈ S h ( x i , x j , ⋯ ) ] , (1) where S denotes the CEMS, q denotes the number of the minimum path sets in the CEMS, h = 1,2 , ⋯ , q , and S h = ( x i , x j , ⋯ ) , x i , x j … is all the units in a minimum path set.

Therefore, the life simulation value t S for a CEMS can be obtained through logical operation after each simulation. t S = max 1 < h < q ( t S h ) = max 1 < h < q [ min t i , t j … ∈ S h ( t i , t j , … ) ] , (2) where t S is the life simulation value for the CEMS, t i , t j … is the life simulation value for each unit in the minimum path set, and t S h is the life simulation value for the hth minimum path set.

The minimal path sets of the typical system are shown in Figure 1, and the system has four minimum path sets: ( x 1 , x 5 , x 4 ) , ( x 2 , x 5 , x 3 ) . The life simulation values for the system are obtained by the following equation (Hong-Bo and Guo, 2009): t S = max [ min ( t 1 , t 2 ) , min ( t 3 , t 4 ) , min ( t 1 , t 5 , t 4 ) , min ( t 3 , t 5 , t 2 ) ] . (3)

In general, the CEMS is usually considered to be composed of several simple serial or parallel subsystems, and the following equations are established: t S k = min t i , t j … ∈ S k ( t i , t j , … ) , (4) t S l = max t i , t j … ∈ S l ( t i , t j , … ) , (5) where t S k is the life simulation value for the serial subsystem and t S l is the life simulation value for the parallel subsystem.

For any CEMS, the simulation value array T ( r × n ) for the units can be obtained by carrying out r simulations over n component units. T ( r × n ) = | t 11 , t 21 , ⋯ , t n 1 t 12 , t 22 , ⋯ , t n 2 ⋮ t 1 r , t 2 r , ⋯ , t n r | . (6)

Each simulation can yield a set of life simulation values for the units, and one life simulation value for the system can be obtained by Eq. 2. If r simulations are performed, r life simulation values can be obtained, namely, t S 1 , t S 2 , ... , t S r , and t _Sr is equivalent to the full life test over r CEMSs. If r is large enough, the life distribution types for CEMSs can be counted and deduced from t S 1 , t S 2 , ... , t S r . Therefore, the probability estimation for the reliability measures of the CEMS is taken.

In order to determine the life distribution types for the system, this thesis makes use of the probability graph estimation method and goodness-of-fit test method. First, the probability graph for common distributions is designed and constructed on the computer, and t S 1 , t S 2 , ... , t S r is drawn and fit on the probability graph. Then, the residual sum of squares under each distribution is calculated, and the distribution whose residual sum of squares is minimal is selected as the initial distribution. Finally, the Pearson χ 2 goodness-of-fit test is performed (Praks and Gono, 2011; Aguwa and Sadiku, 2012), and the life distribution type for the system is determined after the hypothesis test is accepted.

An MW wind turbine unit is mainly composed of impeller, gearbox, generator, yaw system, pitch system, brake system, lubrication system, electrical system, and frequency converter in serial connection, and its reliability block diagram (RBD) is shown in Figure 2.

The given life distributions and parameters for all components in the wind power generator unit are listed in Table 3 (Tavner et al., 2007; Spinato et al., 2009; Guo et al., 2012).

According to the relevant parameters of each unit listed in Table 3, the life simulation value for each unit is generated using the simulation conversion equation in Table 1. We set r = 500 , and the 500 life simulation values for units are obtained as listed in Table 4 (excerpt). 500 life simulation values for the system are obtained through logical operation.

1) Initial selection of the life distribution types for the wind turbine.

The characteristic data in Table 5 are processed using probability graph estimation to get the residual sum of squares under different distributions, Q j , as listed in Table 6.

It is shown in Table 6 that the residual sum of squares for the two-parameter Weibull distribution was the least, and it can be used as the initial distribution.

2) Goodness-of-fit test for the initial distribution.

According to the probability graph estimation, m = 1.1925 , η = 13207 h , the Pearson χ 2 test statistics is obtained from Eqs 12 and 13: χ 2 = 58.7058.

Given the confidence probability γ = 0.9 , χ 0.9 2 ( 47 ) is obtained from Eq. 15. χ 0.9 2 ( 47 ) = 60.0312.

We can obtain χ 2 less than χ 0.9 2 ( 47 ) . χ 2 < χ 0.9 2 ( 47 ) .

The hypothesis test is passed at a high probability, and then the life data of this wind turbine follows the two-parameter Weibull distribution, the distribution parameters, the point estimation and the interval estimation of reliability measures can be further solved.

We set t ¯ S , 1 , t ¯ S , 2 , ... , t ¯ S , p as a complete sample, the likelihood function for the two-factor Weibull distribution can be expressed as L ( m , η ) = ∏ k = 1 p m η ( t ¯ S , k η ) m − 1 ⁡ exp [ − ( t ¯ S , k η ) m ] . (17)

The logarithm is taken from ln ⁡ L ( m , η ) = p ⁡ ln m η + ( m - 1 ) ∑ k = 1 p ln t ¯ S , k η − ∑ k = 1 p ( t ¯ S , k η ) m . (18)

The constraint conditions are set as ∂ ln ⁡ L ∂ m = 0 ,       ∂ ln ⁡ L ∂ η = 0.

The likelihood equations are established and arranged as { ∑ k = 1 p ( t ¯ S , k ) m ⁡ ln t ¯ S , k ∑ k = 1 p ( t ¯ S , k ) m − 1 m = 1 p ∑ k = 1 p ln t ¯ S , k η m = 1 p ∑ k = 1 p ( t ¯ S , k ) m , (19) t ¯ S , 1 , t ¯ S , 2 , ⋯ , t ¯ S , 50 are substituted into Eq. 19 to solve the maximum likelihood estimation for the parameter m, η . { m ^ = 1.2427 η ^ = 13114 .

Therefore, T ¯ S = η ^ Γ ( 1 + 1 m ^ ) , (20) where T ¯ S is the maximum likelihood estimation for the average life of the wind turbine.

As m ^ > 1 , gamma function can be approximated as Γ ( 1 + 1 m ^ ) ≈ 1 − 0.5748646 1 m ^ + 0.9512363 ( 1 m ^ ) 2 − 0.6998588 ( 1 m ^ ) 3 − 0.6998588 ( 1 m ^ ) 3 + 0.4245549 ( 1 m ^ ) 4 − 0.1010678 ( 1 m ^ ) 5 . (21)

If m ^ ≤ 1 , the calculation can be simplified using the recursive equation and T ¯ S can be calculated as follows: T ¯ S = 12287.02 h .

With the given t value, the maximum likelihood estimation for the reliability function R ( t ) is R ^ S ( t ) = exp [ − ( t η ^ ) m ^ ] . (22)

Therefore, the maximum likelihood estimation for the failure rate function λ ( t ) is λ ^ S ( t ) = m ^ η ^ ( t η ^ ) m ^ − 1 . (23)

If the confidence probability is set as 1 − α , the confidence interval for the distribution parameter m, η is { m L = m ^ exp ( − u 1 − α 2 1.0490 p − 1 ) m U = m ^ exp ( u 1 − α 2 1.0490 p − 1 ) , (24) { η L = η ^ exp ( − u 1 − α 2 1.0810 m ^ p − 1 ) η U = η ^ exp ( u 1 − α 2 1.0810 m ^ p − 1 ) , (25) where u 1 − α 2 is the quantile of standard normal distribution, and it can be calculated by Eq. 16.

To simplify the calculation, the shape parameter m can be taken as its maximum likelihood estimation, and it is taken as m ^ = 1.2427 .

When the given confidence probabilities 1 − α are 0.9, 0.8, and 0.7, respectively, the confidence interval for the distribution parameter m, η , and lower confidence limit for the average life of the wind turbine are shown in Table 7. ( R S ( t ) ) L = exp [ − ( t η L ) m ^ ] . (26)

The upper confidence limit ( 1 − α ) for its failure rate is approximated as ( λ S ( t ) ) U = m ^ η L ( t η L ) m ^ − 1 / h , (27) when t is 1000, 1500, 2000, 2500, and 3000h, and the confidence probabilities 1 − α are 0.9, 0.8, and 0.7, respectively, the reliability lower confidence limit ( 1 − α ) and failure rate upper confidence limit ( 1 − α ) of wind turbine are shown in Figure 3 and Figure 4, respectively.

It can be seen from Figure 3 and Figure 4 that all lower limits of reliability ( R S ( t ) ) L show a downward trend, while all failure rates show an upward trend under different confidence probabilities.

For the CEMS with many limiting factors, this article presented a method of replacing life test with computer analog simulation for CEMS. The life simulation values of system constituent units of CEMS could be obtained by their life distribution types and relevant distribution parameters. The life simulation values of CEMS were obtained according to the reliability logic relationship and reliability logic block diagram of the CEMS, and the simulation was instead of the life test of CEMS when the simulation times were enough. The method proposed in this thesis was of great significance to save test time and test cost, especially for the CEMSs that could not carry out reliability life test.

The method proposed in this thesis could solve the problems of system reliability synthesis and prediction of CEMS with different distributions of units. When the “pyramid model” was used for system reliability level by level synthesis and prediction. Compared with the traditional life test method of CMES, it could save a lot of test cost, test time, test site, and so on. This approach could be applied to the early program design, prototype development, trial production, and other stages of production. In engineering applications, the CMES reliability logic block diagram could be drawn by computer, and all the smallest path sets could be obtained by the node traversal optimization algorithm. Therefore, it has obvious application value in engineering.

Finally, the proposed method was applied to the wind turbine, and the reliability lower confidence limit ( 1 − α ) and failure rate upper confidence limit ( 1 − α ) of wind turbine were calculated with given the life time t. The method is also applicable to other CEMSs and can be utilized to provide guidance for system design, maintenance planning, and so on. This paper provides an effective and flexible method for reliability synthesis and prediction of CEMS, which can be easily implemented in engineering practices.