Reliability Assessment via Combining Data from Similar Systems

1. Introduction

Reliability assessment serves as a critical element in equipment test and evaluation, forming an integral component of general quality characteristic assessment. This evaluation dimension maintains inherent connections with maintainability, supportability, and availability analyses, collectively determining operational suitability of equipment systems [1]. While reliability metrics such as Mean Time Between Failure (MTBF), mission reliability, and mean flight hours between failure vary across systems under test (SUT), their assessment remains fundamental to system validation.

Traditional reliability assessment methodologies typically require extensive failure datasets, presenting significant challenges in operational testing contexts. Unlike dedicated reliability trials, operational tests prioritize mission-specific evaluations within constrained timeframes and resource allocations [2]. Consequently, the reliability assessment objectives in operational testing have shifted from comprehensive performance mapping to maximizing evaluation effectiveness under specific mission [3]. This paradigm necessitates developing optimized test strategies that balance data availability with assessment rigor, particularly crucial for time-sensitive procurement programs facing compressed testing schedules [4].

Advances in data collection and storage technologies have expanded potential data sources for reliability assessment beyond conventional field tests. Current approaches increasingly incorporate historical failure data from analogous systems, performance metrics, and expert knowledge to complement operational test data. Among these supplementary sources, data from similar systems proves particularly valuable due to inherent design commonalities in modern equipment development. Modular architectures and evolutionary design practices often create structural, functional, and operational similarities between SUTs and previous-generation systems, both within product families and across different lineages [5].

Effective utilization of similar system data for reliability assessment necessitates addressing two critical challenges: (1) quantifying similarity metrics between reference systems and the SUT, and (2) developing appropriate reliability modeling frameworks. Current similarity quantification methods predominantly employ statistical techniques like grey correlation analysis [6], Pearson’s correlation coefficient [7], Mahalanobis distance [8], Euclidean distance [9], and included angle cosine [10]. However, these approaches require substantial datasets and demonstrate limited efficacy under sparse data conditions. Inheritance factor methodologies attempt to identify common system elements through expert-driven value assignment [11], though subjective bias remains a concern. Regarding model development, the similar product method can elucidate the influence of environmental and derating factors on system failure rates, thereby establishing a correlation between the failure rates of similar systems [12]. However, this approach does not allow for the full utilization of the system’s reliability data. By integrating the system development mechanism and physical parameter analysis, the life equivalence method can equalize the failure time of similar systems [13]. It is essential that the evaluator possess a comprehensive understanding of the system’s operational mechanism and physical characteristics. The application of fuzzy theory [14] and evidence theory [15] introduce non-probabilistic theory into the data analysis process of similar systems. Bayesian theory [16] is a widely utilized approach in the data analysis of similar systems. It employs a prior distribution to express the information of similar systems, combines the data of the SUT via the inheritance factor, or combines different prior distributions with the data model of the SUT, respectively, and then combines the evaluation results to form a comprehensive conclusion [17]. However, the calculation process of this method contains some expressions of complex forms that are difficult to solve, which presents a challenge in practical applications.

To address these limitations in operational testing scenarios, this paper puts forth a methodology for the combination of data from similar subsystems with the aim of conducting a reliability assessment for subsystem under test (SSUT). The methodology encompasses three phases: (1) Mission decomposition using hierarchical method to identify Mission Essential Functions (MEFs) and corresponding subsystems, (2) Phylogenetic Tree (PT) representation of subsystem Bill of Materials (BOM) with environmental considerations, followed by Mixed Integer Programming (MIP)-based similarity factor computation, and (3) Bayesian reliability modeling incorporating similarity covariates, with posterior distribution sampling via Markov Chain Monte Carlo (MCMC) techniques. The No-U-Turn Sampler (NUTS) algorithm enables efficient sampling from the posterior distributions. The following is a presentation of three innovations contained in this article.

(1) The impact of subsystem configuration and external circumstances on the reliability of basic components is taken into account in the similarity factor calculation methodology. PT is employed to illustrate the subsystem structure, the process of calculating the distance between subsystems is subjected to rigorous analysis, and MIP is utilized to facilitate the expeditious resolution of the similarity factor.

(2) A reasonable reliability assessment model is established, incorporating similarity factors, failure data, and regression coefficient priors. The NUTS-enhanced MCMC sampling ensures robust posterior estimation.

(3) Implements system mission decomposition to establish precise data combination targets, aligning with operational test requirements.

The subsequent sections are organized as follows: Section 2 presents the mission decomposition framework for identifying subsystems requiring data combination. Section 3 details the calculation process of similarity factor, including the construction of PT representation, the calculation of distance, and the application of MIP. Section 4 describes the reliability assessment model. Section 5 serves to validate the proposed methodology and demonstrate its superiority over existing techniques through the presentation of two illustrative examples. Section 6 briefly summarizes the entire paper and draws some conclusions.

2. Hierarchical Decomposition of the SUT’s Mission

In operational test, SUT is required to perform the designated test activities in accordance with the mission profile, which delineates the SUT’s tasks, activities, frequencies, periods, environmental conditions, and operational contexts at each stage. Additionally, the mission profile specifies the operational activities and environmental conditions pertinent to specific tasks, organized according to a timeline [18]. Therefore, the SUT in this context can be considered as a Phased-Mission-System (PMS), which may have different functional structures in different mission phases, requiring different subsystems to perform different MEFs, where MEF refers to the minimum set of functions that must be performed to accomplish the specified operational tasks [19]. The SUT in operational test is given one or more missions to describe the critical tasks that must be accomplished to meet the user’s needs. The aim of the mission is explicit, while the specific ways of accomplishing the mission need not be explicit. For example, a certain type of carrier-based aircraft is required to provide close air fire support for a ship formation. If there is insufficient test data for the SUT and other data sources must be found, it is not possible to measure similarity for the entire system, because only certain subsystems or components may be similar between the SUT and the comparison system. As an illustration, the reliability of the 57-mm gun on the Littoral Combat Ship (LCS) can be assessed by combining the 57-mm gun test data from the National Security Cutter (NSC). Despite the fact that it may have been modified to fit different platforms, there are still considerable areas of commonality. The mission of the SUT is decomposed hierarchically into MEFs, which are then mapped to the subsystems. By comparing the mission profile with the decomposition structure, the subsystems involved in operational test can be clearly presented, thus allowing the similarity factors between different subsystems to be targeted. The tree structure of mission decomposition is illustrated in Figure 1.

As illustrated in Figure 1, the mission of SUT is structured in accordance with the mission$\to$task$\to$subtask$\to$metatask$\to$MEF hierarchy, with a corresponding relationship between the MEFs and the subsystems. Furthermore, subtasks may also exist at multiple levels, and metatask is the task that do not necessitate further decomposition at a specific level of granularity. In operational tests, due to the characteristics of the PMS, the object of data combination is subsystem.

3. Similarity Factor Calculation

This section provides a detailed account of the methodology employed in calculating the similarity factor. It begins by outlining the process of establishing a PT representation of the subsystem. This is followed by a description of the approach used to calculate distance, taking into account the external circumstance. The subsequent step is to present the MIP solution for calculating distance and deriving the similarity factor.

3.1. Establishment of PT Representation

The calculation of similarity factor requires the establishment of a suitable representation model to quantify the relevant elements. In addition to its application in the field of reliability assessment, similarity analysis has also been the subject of considerable interest among experts and scholars in the domains of design and manufacturing. Group technology [20] represents a traditional methodology for the analysis of product similarity. However, it necessitates a considerable investment of effort to ensure the continual refinement of the coding system. Similarity analysis based on functional features has also been developed in the field of design retrieval [21]. The challenge in applying it lies in the representation of codes for identifying functional features and the calculation of similarity based on functional feature ranking. Following years of practice and investigation, it has become evident that reliability is contingent upon the design, is refined during production, and is manifested in use [22]. The product structure can be conceptualized as a form of data that is shaped during the extended design and manufacturing process, reflecting the composition and interrelationships of the product’s internal components. Consequently, the product structure is inextricably linked to the reliability level, and research into its similarities has also garnered significant interest [23].

BOM was originally proposed by Dr. Orlicky [24] as a form of product data structure subordinate to the Material Requirements Planning (MRP) system commonly used for production planning and inventory control. BOM is a list of components, raw materials and their quantities required to produce the final product, which can be represented by an unordered root tree. Literatures [25,26,27] employ BOM-based similarity metrics to ascertain the similarity between products. However, factors such as material composition, test circumstance, and operator level may influence product reliability data and thus warrant further investigation.

Figure 2 illustrates the BOM tree structure of two three-level subsystems, comprising eight and seven distinct basic components (numbered 8,9,5,2,6,10,11,12), respectively. The number adjacent to each component indicates the quantity required to form a parent component. For instance, component 1 is the parent of component 4 in T1, and three component 4 are necessary to produce one component 1. Figure 2 is merely an illustrative example, and the actual scenario may encompass a greater number of components and hierarchies.

A BOM tree can be transformed into a PT, thereby facilitating the calculation of similarity factors. PT, also referred to as molecular evolutionary tree, is a method utilized in bioinformatics for the classification of a group of species that have descended from a common ancestor, as well as for the study of their evolutionary history [28]. Figure 3 depicts two PTs corresponding to Figure 2. The arbitrarily ordered numbers “①②③④⑤” indicate nodes, and the number of components associated with each node is shown in parentheses adjacent to the basic components.

3.2. Distance Calculation Principle

Robinson-Foulds (RF) distance is a widely utilized metric for comparing PTs [29]. However, it is unable to account for quantitative differences between nodes and components. This limitation can be addressed by performing a bijective mapping between the subsystems [26]. Taking the PTs in Figure 3 as an example, each node in T1 is mapped to the most similar node in T2, and this process is repeated for each node in T2. The resulting mapping distances are then summarized in Table 1 and Table 2.

In Table 1, column 1 contains all the nodes in T1, and column 2 contains the nodes in T2 that are most similar to the nodes in column 1. This can be determined visually in simple PTs. Column 3 shows the components that are present in T1 but not in T2. Column 4 enumerates the components that are more in T1 than in T2, among those that are shared by the corresponding nodes. Column 5 contains the hierarchical level differences between the components contained in columns 3~4 and the corresponding nodes. For example, if node ① is marked at level 0, then components 6, 10, and 11 are, in turn, marked at level 2, 3, and 3, respectively. Column 6 represents the number of components in column 3 and those in column 4 in T1 over T2, respectively. Column 7 denotes the mapping distance of the corresponding node, which is the summation resulting from dividing the corresponding elements in columns 6 and 5. The meanings of the columns in Table 2 are analogous to those expressed above.

As evidenced by Table 1 and Table 2, the distance between corresponding nodes is determined by the discrepancy in components and the structural characteristics of the node, irrespective of the designation of the intermediate node. This is a significant attribute of PT. For example, the node ⑤ in Figure 3 corresponds to component 7 in Figure 2. However, the components that comprise it in the two subsystems are distinct, which may introduce ambiguity in other similarity analysis methods. The question of whether the numbering of component 7 in the two subsystems should be consistent or not will yield disparate results. The method presented in this paper will remain unaffected by this discrepancy. In addition, the basic distance between different components is 1, which needs to be multiplied by the quantity difference and divided by the hierarchical level difference. For different components, an increase in quantity will correspondingly increase the distance between nodes, while the level difference is equivalent to the weight, and the smaller the difference, the greater the weight. In accordance with the data presented in Table 1 and Table 2, the distance between the subsystems can be calculated as the sum of all distances in both mapping directions, resulting in a total distance of 21.66.

In addition to being determined by intrinsic attributes such as structure, material, function, design principles, and working principles, similarity factor may also be affected by external factors such as the working circumstance, the operator level, and interference from other subsystems within the SUT. These external factors can be considered as the external circumstance. In the PT, it can be posited that the basic components with the same number have essentially the same intrinsic attributes. Consequently, the same components operating in disparate external circumstance can be regarded as representing different working modes of the components. Literature [23] presents a methodology for calculating the product similarity of components in disparate operational modes, which entails initially calculating the similarity factor of the components across various operational modes and subsequently multiplying the number of components in the compared products by the aforementioned similarity factor to derive the equivalent number of components. Nevertheless, when the number of components in the compared product is greater than that in the original product, it results in a situation where the similarity between the two products is greater than it would be if the working mode were not taken into account. This is evidently illogical. In this paper, it is assumed that a similar factor ${\rho }_i\left(0<{\rho }_i\le 1\right)$ has been obtained for the influence of the external circumstance on the reliability of two components numbered i in two PTs. The less quantity is multiplied by ${\rho }_i$ and regarded as the equivalent quantity. If the quantities are equal, either one is multiplied by ${\rho }_i$. Take Figure 3 as an example. Assuming that $\left({\rho }_8,{\rho }_9,{\rho }_5,{\rho }_6,{\rho }_{11},{\rho }_{12}\right)=\left(0.95,0.85,0.78,0.89,0.91,0.88\right)$, the aforementioned method can be employed to modify Figure 3 into Figure 4. At this time, the distance is 36.065, which shows that the distance between subsystems will increase when considering the influence difference of the external circumstance.

The value of factor ${\rho }_i$ affects the final distance and should therefore be carefully considered in practical applications. The subjective experience of engineering experts is an important type of information, and it is a credible practice to assign a value to it after analysis and discussion.

3.3. MIP Solution of Distance and Determination of Similarity Factor

The preceding section 3.2 presents the fundamental principle for calculating the distance between subsystems, taking into account the external circumstance. In operational test, the SSUT is typically more complex, making it challenging to identify the most similar node through observation when performing node mapping, particularly when the external circumstance is taken into account. In such case, the distance can be solved using MIP.

3.3.1. Matrix Encoding Corresponding to PT

MIP is an extension of linear programming, wherein decision variables encompass both continuous and discrete types [30]. Prior to solving the distance using MIP, the PT must be converted to matrix form. The conversion process is outlined below:

step 1: Create matrices ${\mathit{M}}_1={\left[m_{ij}\right]}_{p\times q}$ and ${\mathit{M}}_2={\left[w_{bj}\right]}_{k\times q}$ corresponding to two PTs. The $i$-th ($b$-th) row of ${\mathit{M}}_1$(${\mathit{M}}_2$) corresponds to the node $i$($b$) of the PT. $p$($k$) is the total number of nodes, while $q$ denotes the total number of component numbers. All component numbers should be arranged in ascending order, and the values assigned to the matrix ${\mathit{M}}_1$ should be determined in accordance with the following Equation:

$$\left\{\begin{array}{l}{m}_{ij}=0 \mathrm{node} i \mathrm{does} \mathrm{not} \mathrm{contain} \mathrm{the} j-\mathrm{th} \mathrm{component} \\ m_{ij}=1 \mathrm{else} \end{array}\right.$$

(1)

Subsequently, the matrix ${\mathit{M}}_2$ should be assigned values in a similar manner.

step 2: In the matrix ${\mathit{M}}_1$ resulting from step 1, the element $m_{ij}$ is to be replaced with $m_{ij}/l_{ij}$, where $l_{ij}$ is the hierarchical level difference between the j-th component and node i. The same replacement process is to be applied to the elements in matrix ${\mathit{M}}_2$.

step 3: For matrix ${\mathit{M}}_1$ from step 2, multiply the elements in the j-th column by the equivalent quantity of the j-th component in the corresponding subsystem (taking into account the external circumstance). The same process is to be repeated for the elements in matrix ${\mathit{M}}_2$.

PTs depicted in Figure 4 can be converted into two matrices by following the aforementioned steps:

$$M_1=\left[\begin{array}{ccccccccc}2& 4.5& 4.45& 11.4& 2.55& 4/3& 20/3& 10/3& 0\\ 0& 9& 0& 17.1& 3.825& 0& 0& 0& 0\\ 0& 0& 0& 34.2& 7.65& 0& 0& 0& 0\\ 0& 0& 8.9& 0& 0& 2& 10& 5& 0\\ 0& 0& 0& 0& 0& 4& 20& 10& 0\end{array}\right]$$

(2)

$$M_2=\left[\begin{array}{ccccccccc}0& 2.34& 5& 12& 3& 0& 18.2/3& 8.8/3& 4/3\\ 0& 4.68& 0& 18& 4.5& 0& 0& 0& 0\\ 0& 0& 0& 36& 9& 0& 0& 0& 0\\ 0& 0& 10& 0& 0& 0& 9.1& 4.4& 2\\ 0& 0& 0& 0& 0& 0& 18.2& 8.8& 4\end{array}\right]$$

(3)

3.3.2. MIP Construction and Similarity Factor Calculation

Construct the following MIP optimization model based on matrices ${\mathit{M}}_1$ and ${\mathit{M}}_2$:

$$\begin{array}{l}\mathrm{Minimize} {\displaystyle \sum _{i=1}^p{\displaystyle \sum _{b=1}^k{x}_{ib}}}+{\displaystyle \sum _{b=1}^k{\displaystyle \sum _{i=1}^p{y}_{bi}}}\\ s.t.\left\{\begin{array}{l}{x}_{ib}+\mathrm{C}\left(1-z_{ib}\right)\ge {\displaystyle \sum _{j=1\left|m_{ij}\ge w_{bj}\right.}^q\left(m_{ij}-w_{bj}\right),\forall i=1,\mathrm{...},p,\forall b=1,\mathrm{...},k}\\ y_{bi}+\mathrm{C}\left(1-e_{bi}\right)\ge {\displaystyle \sum _{j=1\left|w_{bj}\ge m_{ij}\right.}^q\left(w_{bj}-m_{ij}\right),\forall b=1,\mathrm{...},k,\forall i=1,\mathrm{...},p}\\ {\displaystyle \sum _{b=1}^k{z}_{ib}=1,\forall i=1,\mathrm{...},p}\\ {\displaystyle \sum _{i=1}^p{e}_{bi}=1,\forall b=1,\mathrm{...},k}\end{array}\right.\end{array}$$

(4)

where the meanings of $p,k,q,m_{ij}$ and $w_{bj}$ are the same as in section 3.3.1. If node i is mapped to node b, the value of $z_{ib}$ is 1, otherwise it is 0. If node b is mapped to node i, the value of $e_{bi}$ is 1, otherwise it is 0. $x_{ib}$ indicates the distance from node i to node b, and $y_{bi}$ indicates the distance from node b to node i. The objective function is the sum of the distances of all mappings, i.e. the distance between the two subsystems, in Equation (4). The first two inequalities within the constraint define $x_{ib}$ and $y_{bi}$, respectively, while the decision variables $z_{ib}$ and $e_{bi}$ in conjunction with constant C, which takes a large value, are used to ensure that $x_{ib}$ and $y_{bi}$ are assigned 0 when the distance between two nodes that are not mapped to each other is to be computed. The final two equations of the constraint serve to guarantee that each node is mapped to a single node. The optimization problem shown in Equation (4) can be expeditiously resolved during the programming implementation via the Python language by importing the function ‘LpProblem’ from the module ‘pulp’ [31]. Upon feeding matrices ${\mathit{M}}_1$ and ${\mathit{M}}_2$ into the MIP model, the resulting distance between the two subsystems is 36.071. The similarity factor between them is:

$$s={\displaystyle \frac{D_{\max }-D}{D_{\max }}}$$

(5)

where $D_{max}$ represents the maximum possible distance between two subsystems, which is equal to the sum of all elements of matrices ${\mathit{M}}_1$ and ${\mathit{M}}_2$. This represents the situation when all nodes do not match at all. $D$ is the result of the distance between the subsystems. For matrices ${\mathit{M}}_1$ and ${\mathit{M}}_2$, D_max = 329.26, s = 0.89.

4. Reliability Assessment Model

Assuming that there are I similar subsystems for a given SSUT, collect test data for these I+1 subsystems and establish a data model $\left(t_i,{\delta }_i,s_i\right)$, where $i$=1,$\cdots$,n, and n is the data size. $t_i$ is the observation time for a subsystem, ${\delta }_i$ is an indicator that takes 1 if $t_i$ is a failure data and takes 0 if $t_i$ is a right-censored data. $s_i$ is the similarity factor of the subsystem to which the data belongs, and takes 1 when the data belongs to the SSUT.

As the Weibull distribution is a prevalent tool in the field of reliability engineering and is intimately associated with other distribution types, including the exponential distribution and the Gamma distribution [32], this paper employs a two-parameter Weibull distribution to model the failure time, utilizing the following probability density function (p.d.f.):

$$f(t\left|m,\lambda \right.)=m{t}^{m-1}\exp \left(\lambda -\exp (\lambda )t^m\right)$$

(6)

where m>0 is the shape parameter, $\lambda =\mathit{log}(\eta )$ and $\eta$>0 is the scale parameter, t>0 is the failure time. The reliability function and the mathematical expectation of failure time are

$$R(t\left|m,\lambda \right.)=\exp \left(-\exp (\lambda )t^m\right)$$

(7)

$$\mathrm{E}(\mathrm{T})=\exp (-{\displaystyle \frac{\lambda }{m}})\mathsf{\Gamma }\left(1+{\displaystyle \frac{1}{m}}\right)$$

(8)

where $\mathsf{\Gamma }\left(\cdot \right)$ is the Gamma function. While there are some discrepancies in the failure time distribution parameters between similar subsystems, their test data can be regarded as originating from the same population at a specific significance level, and there should also be some kind of correlation between the distribution parameters. In this paper, we employ the corresponding trick of Weibull survival regression to consider the similarity factor as a covariate, thereby expressing the correlation as follow:

$$\left\{\begin{array}{l}{m}_i=\exp \left({\alpha }_1+{\beta }_1{s}_i\right)\\ {\lambda }_i={\alpha }_2+{\beta }_2{s}_i\end{array}\right.$$

(9)

Subsequently, the likelihood of parameters α₁, β₁, α₂ and β₂ is derived as follow:

$$\begin{array}{l}L({\alpha }_1,{\beta }_1,{\alpha }_2,{\beta }_2\left|D\right.)={\displaystyle \prod _{i=1}^n\left[\exp \left({\alpha }_1+{\beta }_1{s}_i\right){t_i}^{\exp \left({\alpha }_1+{\beta }_1{s}_i\right)-1}\right.}\hfill \\ {\left.\cdot \exp \left({\alpha }_2+{\beta }_2{s}_i-\exp ({\alpha }_2+{\beta }_2{s}_i){t_i}^{\exp \left({\alpha }_1+{\beta }_1{s}_i\right)}\right)\right]}^{{\delta }_i}\\ \cdot {\left[\exp \left(-\exp ({\alpha }_2+{\beta }_2{s}_i){t_i}^{\exp \left({\alpha }_1+{\beta }_1{s}_i\right)}\right)\right]}^{1-{\delta }_i}\end{array}$$

(10)

where D=$\left(t_i,{\delta }_i,s_i\right)$. The prior distributions of parameters α₁, β₁, α₂ and β₂ can all be assumed to be normal distributions, specifically N(μ₁, σ₁²), N(μ₂, σ₂²), N(μ₃, σ₃²) and N(μ₄, σ₄²). In accordance with Bayesian theory, the joint posterior distribution of (α₁, β₁, α₂, β₂) is given by:

$$\begin{array}{c}\pi \left({\alpha }_1,{\beta }_1,{\alpha }_2,{\beta }_2\left|D\right.\right)\propto L({\alpha }_1,{\beta }_1,{\alpha }_2,{\beta }_2\left|D\right.)\cdot \varphi ({\displaystyle \frac{{\alpha }_1-{\mu }_1}{{\sigma }_1}})\cdot \varphi ({\displaystyle \frac{{\beta }_1-{\mu }_2}{{\sigma }_2}})\\ \cdot \varphi ({\displaystyle \frac{{\alpha }_2-{\mu }_3}{{\sigma }_3}})\cdot \varphi ({\displaystyle \frac{{\beta }_2-{\mu }_4}{{\sigma }_4}})\end{array}$$

(11)

where $\varphi$ is the standard normal distribution p.d.f.. The complex form of Equation (11) makes it difficult to give a closed expression, which can be sampled from the joint posterior distribution using MCMC. Specifically, No-U-Turn Sampler (NUTS) is employed to execute the sampling process. NUTS is an extension of Hamilton Monte Carlo (HMC) algorithm [33]. It consists of a sequence of steps based on first-order gradient information, which can circumvent the random wandering behavior and sensitivity to related parameters exhibited by many alternative MCMC methods. Furthermore, it can converge more rapidly to high-dimensional target distribution. Bayes point estimates for the parameters α₁, β₁, α₂ and β₂ can be readily obtained from the posterior samples.

It is of greater significance that Bayes point estimates and credible intervals for reliability and failure time can be obtained from the posterior sample. The specific steps are as follows:

step 1: Let i = 1.

step 2: Use NUTS to draw samples from the posterior distribution shown in Equation (11).

step 3: Convert samples ${\alpha }_1^i,{\beta }_1^i,{\alpha }_2^i$ and ${\beta }_2^i$ to samples R_i and E_i using Equations (7) ~ (9), and let i = i+1.

step 4: Repeat steps 2 ~ 3 until i = N to obtain samples {R_i, i = 1,$\cdots$,N} and {E_i, i = 1,$\cdots$,N}.

Arrange the elements in {R_i, i = 1,$\cdots$,N} in ascending order as R₁≤R₂≤$\cdots$R_N，and do the same for {E_i, i = 1,$\cdots$,N}. The resulting Bayes point estimates of both reliability and failure time are

$$\widehat{R}={\displaystyle \frac{1}{N}}{\displaystyle \sum _{i=1}^N{R}_i}, \widehat{\mathrm{E}}\left(\mathrm{T}\right)={\displaystyle \frac{1}{N}}{\displaystyle \sum _{i=1}^N{\mathrm{E}}_i}$$

(12)

The credible intervals for reliability and failure time under credible level (1-α) are

$$\left[R_{\alpha N/2},R_{N-\alpha N/2}\right],\left[{\mathrm{E}}_{\alpha N/2},{\mathrm{E}}_{N-\alpha N/2}\right]$$

(13)

5. Numerical Experiments

5.1. Application of the Proposed Method

This subsection provides an application of the proposed method, utilizing the operational test of a particular type of equipment as a case study. The mission of this equipment is assumed to be breaching opponent’s defensive position, which is decomposed hierarchically in accordance with the method described in Section 2, and the simplified tree structure is illustrated in Figure 5.

The lower-level decompositions of the “ Buildup” and “Assault “ tasks are simplified in Figure 5 because this subsection employs the “Hydraulic subsystem” under the “Firepower preparation” task branch as an illustrative example. It is comprised of hydraulic pump, control component, and hydraulic connector and pipe. Additionally, five distinct variants of the equipment are available, each exhibiting variations in the hydraulic subsystem. The hydraulic subsystems are designated as H1~H6, with H1 representing the SSUT. Their BOM tree structures are illustrated in Figure 6. Some components in the BOMs are simplified. Some identical or inconveniently disclosed components are removed, and the different components are highlighted to maintain the complexity of the method verification process. However, it does not indicate that the proposed method is incapable of handling more complex hierarchies. The six hydraulic subsystems contain a total of 21 different components, each of which is assigned a distinct number, as illustrated in Table 3. The PT corresponding to each hydraulic system is depicted in Figure 7.

Let H1 be paired with H2 ~ H6 respectively, and calculate the similarity factors according to the method described in Section 3. The distances between H1 and H2 ~ H6 can be calculated as: 15.67, 21.33, 12.83, 15.50, 42.50. The maximum possible distances are 212.33, 225.33, 228.83, 220.50, 217.50. Substituting the value of D_max as 228.83 into Equation (5) yields the similarity factors, which are 0.93, 0.91, 0.95, 0.94 and 0.81, respectively. The result can also be visually perceived from Figure 7. Compared with H1, H6 is observed to lack components “2” and “13” and to have additional “20” and “21”. This discrepancy is more pronounced than that observed in the other subsystems, providing further evidence that this methodology is sound.

Following a period of deliberation, the engineering experts arrived at a consensus regarding the external circumstance impact on the same components, which is expressed in matrix form as follow:

$$\left[\begin{array}{ccccccccccccccc}0.86& 0.92& 0.95& 0.78& 0.95& 0.88& 0.82& 0.98& 0.74& 0.88& 0.77& 0& 0.75& 0.89& 0.91\\ 0.78& 0.75& 0.83& 0.76& 0.99& 0.89& 0.98& 0.94& 0.72& 0.97& 0.78& 0& 0.77& 0.93& 0.88\\ 0.87& 0.85& 0.88& 0.84& 0.83& 0.86& 0.92& 0.88& 0.75& 0.94& 0.94& 0.71& 0.78& 0.77& 0.98\\ 0.88& 0& 0.98& 0.91& 0.72& 0.92& 0.94& 0.93& 0.88& 0.79& 0.96& 0.73& 0.98& 0.83& 0.99\\ 0.99& 0& 0.89& 0.91& 0.88& 0.96& 0.82& 0.88& 0.76& 0.93& 0.81& 0.71& 0& 0.80& 0.75\end{array}\right]$$

(14)

The element in the i-th row and the j-th column of the matrix represents the similarity factor between two components numbered j in the hydraulic subsystems H(i+1) and H1 with respect to the influence of the external circumstance. In considering the impact of the external circumstance on reliability, the distances between H1 and H2 ~ H6 can be solved as: 29.67, 29.05, 44.11, 45.94, 59.57. The similarity factors are 0.87, 0.84, 0.81, 0.80, 0.76, respectively.

Figure 8 illustrates the data pertaining to the six hydraulic subsystems. “×” represent the failure time data, “•” represent the right-censored data, and subsystem number 1 corresponds to SSUT. Based on engineering experience, it can be inferred that the data follow a Weibull distribution.

NUTS is employed for sampling from the joint posterior distribution as shown in Equation (11). The prior distributions N(μ₁, σ₁²), N(μ₂, σ₂²), N(μ₃, σ₃²) and N(μ₄, σ₄²)are all taken as N(0, 100²), indicating that four fairly diffuse priors are employed. The joint posterior distribution will be primarily influenced by the test data illustrated in Figure 8. Figure 9 depicts the posterior distributions of the parameters, accompanied by the posterior mean and the highest density interval (HDI) [34] for each parameter. Figure 10 illustrates the reliability function of the SSUT over time, accompanied by the associated credible interval. The estimate of the expected failure time is 170.45h, with a 95% credible interval of [126.50h, 234.04h]. It should be noted that only one failure time is available for the SSUT in Figure 8. Consequently, the maximum likelihood estimation method will no longer be applicable in the case where the failure time follows a Weibull distribution.

5.2 Comparison of Methods

The calculation method of similarity factor is compared with the Jaccard similarity measure. The Jaccard coefficient J_ij between products i and j is calculated using the following equation:

$$J_{ij}={\displaystyle \frac{N_s}{N_s+N_b+N_c}}$$

(15)

where N_s is the number of components shared between products i and j. N_b is the number of components belonging to i but not j, and N_c is the number of components belonging to j but not i. J_ij is expressed only as the ratio between the number of shared components and the total number of components, regardless of the relationship between components. The Jaccard coefficients between H1 and H2~H6 are calculated using Equation (15), yielding values of 0.85, 0.81, 0.89, 0.87 and 0.65, respectively, in the absence of external circumstance. The results of the two methods indicate a consistent ability to distinguish similarities between H1 and H2~H6. However, a notable discrepancy is observed in the Jaccard coefficient between H6 and H1, which is too small and may result in the exclusion of H6 data during data collection. As illustrated in Figure 7, H6 and H1 exhibit a consistent hierarchy. Discrepancies exist in the constituent components 2, 15, 20 and 21. The weights of them differ depending on the hierarchical level. The Jaccard coefficient does not distinguish between them well, and exaggerates the weights of 2 and 20.

Simulation is used to verify the accuracy of the reliability assessment method in section 4. Suppose that for SSUT, there are two similar subsystems. The similarity factors between subsystems 1~2 and SSUT are set to 0.78 and 0.95, respectively. The parameters in Equation (9) are set as α₁=0.36, β₁=0.45, α₂=1.2 and β₂=0.67, and the failure time distribution parameters of the three subsystems are then solved. Subsequently, the failure time data for the subsystems are generated, with SSUT collecting one data point and subsystems 1~2 collecting six and eight data points, respectively. Finally, the failure time distribution parameters of SSUT are estimated according to the method described in section 4, with the prior distributions of α₁, β₁, α₂ and β₂ remaining N(0, 100²).

For comparison, the method described in the literature [16] is also employed for the estimation of the distribution parameters. Let θ = (m, λ), and the prior obtained from the data of the similar subsystem i is ${\pi }_i\left(\mathit{\theta }\right)$, i=1,2. Subsequently, the priors are combined equally into a single prior distribution $\pi \left(\mathit{\theta }\right)=0.5{\sum }_{i=1}^2\left[s_i{\pi }_i\left(\mathit{\theta }\right)+1-s_i\right]$, where s_i is the “inheritance factor” and is equivalent to the similarity factor. The part inside the symbol “[]” of the equation is a mixture prior distribution, which contains a informative prior ${\pi }_i\left(\mathit{\theta }\right)$ and a uninformative prior represented by a uniform distribution. In accordance with Bayesian theory, the posterior distribution $\pi \left(\mathit{\theta }\left|D\right.\right)$ can be derived from the test data D of SSUT and $\pi \left(\mathit{\theta }\right)$, which can be used to estimate $\mathit{\theta }$.

Both the proposed method and the existing method need to employ NUTS for the purpose of sampling from $\pi \left(\mathit{\theta }\left|D\right.\right)$. Given that the status of SSUT and subsystems 1~2 can be interchanged, the parameters of subsystem 1~2 can also be estimated using the same method. To eliminate the influence of random effects, each method is repeated 100 times, and relative deviation between the estimated value and the true value is then subjected to statistical analysis. The results are presented in Table 4. It is evident that the absolute value of relative deviation of the proposed method is less than that of the existing method, indicating that it has a greater accuracy. The primary reason for this discrepancy is that, although the existing method integrates engineering data and expert insights into the determination of inheritance factor, it lacks the utilization of test data to align the relationship between inheritance factor and failure time distribution parameters. In contrast, the proposed method can address this shortcoming, thereby conferring upon the reliability assessment model a degree of learning capacity.

6. Conclusions

In operational testing scenarios, data from similar systems can effectively complement the test data of the SUT. This study proposes a novel reliability assessment framework that combines data from SSUT through systematic data combination. The mission of SUT is initially decomposed to identify SSUT that requires data combination operation. Subsequently, BOM trees of similar subsystems are matched to PTs. MIP is then employed to determine the distance between subsystems, thereby obtaining the similarity factors. The impact of the external circumstance on subsystem components is taken into account during the solution process. Ultimately, a reliability assessment model is constructed, and Bayesian theory and NUTS are utilized to obtain the reliability assessment results. Numerical experiments demonstrate the framework’s practical applicability and superior performance compared to conventional methods. Three conclusions were drawn as follows:

(1) Operational test constraints significantly challenge standalone reliability assessment. Data combination from similar systems emerges as an effective mitigation strategy, particularly beneficial for systems with limited test data availability.

(2) System reliability demonstrates strong dependence on structural configuration and operational environment. Matching the BOM tree of a subsystem to a PT enables MIP to rapidly and accurately solve the similar factors.

(3) It is essential to ascertain the relationship between the similarity factor and failure time distribution parameters via the underlying system mechanisms, and should also be considered for learning from test data.

This work constitutes a critical component of our ongoing research program for operational reliability evaluation. Future efforts will focus on hierarchical reliability aggregation from subsystem to system level, ultimately determining mission reliability through comprehensive verification.