# Health Indicator for Predictive Maintenance Based on Fuzzy Cognitive Maps, Grey Wolf, and K-Nearest Neighbors Algorithms

Authors: G. Mazzuto; S. Antomarioni; F. E. Ciarapica; M. Bevilacqua
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

An essential step in the implementation of predictive maintenance involves the health state analysis of productive equipment in order to provide company managers with performance and degradation indicators which help to predict component condition. In this paper, a supervised approach for health indicator calculation is provided combining the Grey Wolf Optimisation method, Swarm Intelligence algorithm, and Fuzzy Cognitive Maps. The k-neighbors algorithms is used to predict the Remaining Useful Life of an item, since, in addition to its simplicity, they produce good results in a large number of domains. The approach aims to solve the problem that frequently occurs in interpolation procedures: the approximation of functions belonging to a chosen class of functions of which we have no knowledge. The proposed algorithm allows maintenance managers to distinguish different degradation profiles in depth with a consequently more precise estimate of the Remaining Useful Life of an item and, in addition, an in-depth understanding of the degradation process. Specifically, in order to show its suitability for predictive maintenance, a dataset on NASA aircraft engines has been used and results have been compared to those obtained with a neural network approach. Results highlight how all of the degradation profiles, obtained using the proposed approach, are modelled in a more detailed manner, allowing one to significantly distinguish different situations. Moreover, the physical core speed and the corrected fan speed have been identified as the main critical factors to the engine degradation.


---

## Body

## 1. Introduction

Although predictive maintenance practices have existed for many years, only recently, thanks to the emerging Industry 4.0 technologies with increasingly reliable and affordable smart systems, it has become widely accessible [1]. It has several advantages, including machine life increase by 3–5%, reduced maintenance costs by up to 40%, and returns on investment up to 10 times [2].One of the most relevant steps in the prediction process is the choice of the best approach for the item behaviour assessment, such as data-driven or model-driven approach [3]. In particular, according to the platform developed by Patel et al. [4] for the application of Industry 4.0 principles to the industrial reality, the data-analytic layer is crucial to understand a plant functioning. Moreover, if properly designed, it allows users to identify the presence of invisible relations among data provided by the application layer [5]. It is also true that, according to the “no free lunch” theorems, a standard procedure for predictive maintenance does not exist. Still, it must be chosen among those that best suit the reality under analysis [6]. In any case, regardless of the adopted process, for a more accurate and optimal prediction, it is necessary to gather and analyse appropriately large amounts of data within a time frame [7, 8] with consequent problems deriving from the identification of the most accurate health indicators. The health of a system can be defined as the deviation or degradation of an item behaviour from its regular operating performance [9].The calculation of a suitable health indicator (HI) is fundamental to establish a link between the deviation or degradation of an item and its Remaining Useful Life (RUL). Thus, an accurate HI is a key for a more precise prediction tool, guaranteeing its reproducibility [10, 11]. This observation is the reason why many researchers focus their activity on this issue ranging from supervised and unsupervised algorithms [12, 13] to physical [14] and virtual [15] HIs.The HI assessment needs the monitored parameters provided by the physical data from sensors to be transformed into information represented as indicators. The potential benefits include not only the reduction of the quantity of data examined but also the maximisation of the useful information content [16].In this context, the proposed paper lays its foundations. An approach for HI definition and calculation is provided combining the Grey Wolf Optimisation (GWO) approach, belonging to the set of Swarm Intelligence algorithms, and Fuzzy Cognitive Maps (FCMs). Subsequently, the k-neighbors algorithms are used to predict the item RULs.The proposed approach, in comparison to previous studies presented in the literature, does not require knowledge about the gradients of the cost function and constrained functions, guaranteeing both reliable and robust performance and easy implementation. Moreover, it ensures extreme flexibility and adaptability to a given domain. It allows an in-depth understanding of a specific issue; thanks to the possibility of symbolically representing the relationships among all the involved variables.To present the methodology and analyse its performance, the rest of the paper is organised as follows. Section2 briefly describes the analysed literature on Swarm Intelligence algorithms and FCMs application to predictive maintenance. Section 3, divided into three sections, focuses on the explanation of the FCMs theory and GWO algorithm functioning. It then explains the proposed algorithm steps, underlining its benefits. Section 4 shows the results obtained using a dataset provided by NASA regarding the RUL prediction for aircraft engines and compares the results with an Artificial Neural Network approach. Conclusions have been drawn in Section 5.


## 2. Literature Review

As mentioned, research on predictive maintenance has grown in recent years due to the development of Industry 4.0 technology. Hence, to gather the most relevant contributions dealing with maintenance in general and FCMs and Swarm Intelligence applications in particular, a systematical approach has been adopted. The Scopus scientific database has been selected, considering that all the papers have an available full text written in English. All articles have been read to assess their relevance and pertinence to the theme developed in this study. In Table1, the combination of the keywords selected, the number of papers retrieved by Scopus, and the ones chosen for this literature review are reported.Table 1
Summary of the selected literary contributions.
Keywords# of papers# of relevant papers“Predictive maintenance” and “health indicator”2521“Predictive maintenance” and “Fuzzy Cognitive Map(s)” and “Swarm Intelligence”——“Predictive maintenance” and “Fuzzy Cognitive Map(s)”——“Maintenance” and “Fuzzy Cognitive Map(s)”4210“Predictive maintenance” and “Swarm Intelligence”21“Maintenance” and “Swarm Intelligence”1359“Predictive maintenance” and “Grey Wolf”1—“Maintenance” and “Grey Wolf”3310In recent literature, several contributions deal with the development of HIs aiming at predicting the need for maintenance interventions. For example, some authors propose the implementation of dashboards for the monitoring of the equipment health status in the semiconductor manufacturing industry [17, 18], while others focus on structural vibrations analysis [19] and RUL prediction [20, 21]. Various techniques and methodologies can be found in the literary contributions: for instance, Baraldi et al. [22] develop a differential evolution-based multiobjective model aiming at defining the health status of the system and adopting maintenance strategies; other authors, instead, apply artificial neural networks [23] or genetic algorithms [24] to model the health status of the system.To the best of the authors’ knowledge, there is no evidence of scientific papers dealing with predictive maintenance through the application of FCMs and Swarm Intelligence (SI) approaches. At the same time, a contribution can be found only through the application of SI methods. Li et al. [25], indeed, applied a multiclass relevant vector machine—optimised through the application of the SI dragonfly algorithm—to predict the failures of a diesel engine. Other SI applications to the maintenance field can be found in existing literature, for example, Zheng et al. [26] use the particle swarm optimisation to predict the performance degradation of aeroengines, considering aspects such as fuel consumption, rotor vibration, and thrust loss. A similar perspective is adopted by Hu et al. [27], in diagnosing the failures of a gearbox, through the particle swarm optimisation and the kernel extreme learning machine, and by Zhao and Liu [28] who solved the same class of problems through the rough set theory. Several further SI applications to the maintenance field instead focus on the definition of the maintenance scheduling [29–32]; R. [33].Going into detail regarding the GWO algorithm, some applications in the maintenance field can be found in the literature: the majority of them focus on the cost efficiency of the maintenance processes. For example, it is applied to optimise the design and maintenance of photovoltaic power plants [34] or to minimise maintenance costs of heat and power systems [35–38]. Kumar et al. [39] focus on both the reliability and the costs of a Space Shuttle, through the implementation of a multiobjective GWO. Dalla Vedova et al. [40], instead, compare different algorithms, among which the GWO is for the RUL estimation of an aircraft actuator, while Abdelghafar et al. [41] optimise a support vector machine through the GWO to improve the detection of satellite sensor failures. Some works focus on the scheduling through the implementation of the GWO Algorithm: it can be applied to solve job shop and maintenance scheduling problems [42] as well as to block flow shop scheduling, considering fuzzy processing times and dynamic maintenance strategies [43, 44].FCMs have proved to be useful tools in supporting the decision-making processes in the maintenance field. For instance, they can be applied to verify the impact of maintenance activities on a building’s energy efficiency [37] or to identify the factors affecting human reliability during the maintenance operations [45]. According to Gupta and Gandhi [46], data coming from maintenance work orders can be used to detect possible improvement areas in terms of component design. Dynamic risk modelling is also performed through the FMCs: in Lopez and Salmeron study [47], an FCM is built to assess the risk during the enterprise resource planning of maintenance processes, while in Jamshidi et al. [48], it is used to study the critical factors related to the maintenance outsourcing. Damage detection can also be performed through the FMCs. For instance, Senniappan et al. [49] propose an application for the early detection of damages in civil structures’ elements of support based on an FCM, modelling both the knowledge obtained from the domain experts and the existing literature. Instead, Lee et al. [50] use rule-based FCMs based on the experts’ knowledge and experience to identify the factors accelerating the deterioration of rubber components in order to predict the maintenance timing and structure a diagnostic process. In the work of Azadeh et al. [51], the FCM is used to assess which factors among cognitive and temporal ones have a more relevant impact on the execution of the maintenance interventions. Similarly, maintenance errors can be analysed through FCMs in order to highlight which are the most critical and repetitive ones and recommend modifications in the maintenance process or training [52]. Zhang et al. [53], instead, develop a robot dedicated to live maintenance whose behaviour is predicted through an FCM.According to the existing contributions, there is no evidence of the joint implementation of FCMs and GWO, even though both the methodologies have been successfully applied to the maintenance field. Among the benefits harboured by the GWO, its ability to work in a dynamic environment is one of the most useful in this application field. In parallel, the FCMs are useful for the qualitative simulation of a modelled system. To sum up, the benefits harboured by the joint implementation of the two techniques proposed in this research approach are the flexibility and adaptability, as well as the reliability and robustness of the performance.


## 3. The Research Approach

The general scheme of a predictive maintenance procedure proposed in this work is shown in Figure1 and described below.Figure 1
General scheme for a predictive maintenance procedure.(i)
Preprocessing Data. Preprocessing means the preparation of the dataset for analysis; it incorporates all the steps for dataset preparation. In this part of the process, it is essential to get as much information and indications as possible from the dataset.(ii)
Features’ Extraction. It is the step in which variables are selected and/or the amount of data to be processed is reduced, ensuring an accurate and complete description of the original dataset.(iii)
Splitting Data. This is an analytical step to understand how to train the machine learning system in the best way. As a matter of fact, within machine learning systems, there are two parts: the first is the training that, as the name may indicate, trains the course and teaches it how to act. After this step, the system is ready to perform what it has learned and to test if the training completed in the previous step was successful. This is done through the score or test. Given the significance that these two steps assume, it is of fundamental importance to understand the best way to divide the available data package in the right proportions.(iv)
Health Indicator Modelling. The sensor readings, reworked in the previous steps, are combined into a single parameter called health indicator to be used in the prediction of the adverse event.(v)
RUL Prediction. The RUL equipment is carried out in this work through the K-nearest neighbors classifier [54] and Weibull fitting [55].The core activities of this work are the HI definition and RUL Prediction. The innovative proposed methodology to develop these activities will be described in depth in Section3.1.


### 3.1. The HI Modelling

The proposed approach for the HI modelling is performed by the combined use of FCM and a Swarm Intelligence algorithm such as the GWO. Before describing the proposed approach, the FCM and GWO theories are briefly described in the following two sections.


#### 3.1.1. The FCMs’ Modelling

A cognitive map (CM) can be thought of as a concept map reflecting mental processing, comprised of collected information and several cognitive abstractions, individually filtered, about regarding physical phenomena and experiences [56]. Cognitive maps are visual representations of an individual’s mental model constructs, analogous to concept maps for representing human reasoning and knowledge or beliefs [7]. Thus, a generic problem is considered, and an expert panel of experts is formed for its in-depth analysis, since different individuals may face the same question differently. According to their area of expertise through fuzzy logic, they model collective FCM identifying concepts and relationships about regarding the considered problem. In particular, concepts, in number of N, are the FCM key elements that stand for the main characteristics of the abstract mental model for whichever complex system [57]. Once concepts are identified, experts are asked to assign a numerical value wij (the weight of the relation between concept ith and jth) for the W matrix, which represents the influence of concept Ci on concept Cj. According to equation (1), wij ranges in [–1, 1]. Specifically, wij = 0 indicates no causality between concepts, wij > 0 indicates causal Cj increases as Ci increases (or Cj decreases as Ci decreases), and wij < 0 shows causal decrease or negative causality (Cj decreases as Ci increases or Cj increases as Ci decreases):(1)FCM=w1,1…w1,N⋮⋱⋮wN,1…wN,N.Although many studies exist concerning the dynamical representation of an FCM, generally, the experts’ opinion aggregation of expert opinions for the collective weight matrix modelling is performed using the SUM method [58]. Then, overall linguistic weight is evaluated using the centre of gravity (COG) defuzzification method [59]. Some examples are presented by Bevilacqua et al. [7, 60, 61] and Stylios et al. [62] where a unique credibility value is assigned to each expert and a threshold function is used in the aggregation. On the contrary, a modification of the approach mentioned above has been provided by Stylios and Groumpos [63] and Stylios and Groumpos [64], introducing a corrective factor for the experts’ credibility evaluation. However, this approach does not take into consideration the fact that, in a complex multidisciplinary problem, most experts have in-depth knowledge of only parts of the problem and not the entire issue [65].Once the total weights’ matrix,W, has been designed, it is possible to analyse the system behaviour through simulations. Thus, if Ai defines the instantaneous value of concept Ci, its evolution over time can be evaluated computing the influence of the related concepts Cj on the specific concept Ci according to(2)Aik+1=fAik+∑j≠ij=1nAjkwi,j,where Aik + 1 is the value of concept Ci at simulation step k + 1 and Ajk is the value of concept Cj at simulation step k. Also, wij is the weight of the interconnection from concept Cj to concept Ci and f is an appropriate threshold function used to force the concept value to be monotonically mapped into a normalised range [66]. Other equations can be used in place of equation (2) as suggested by Mazzuto et al. [67] and Osoba and Kosko [68].An important topic in the FCM analysis is the indirect and total causal effect evaluation (Axelrod, 1976), whose knowledge allows an in-depth map analysis. The indirect effectIk of Ci concept on Cj concept can be defined as shown in(3)IkCi,Cj=minwCp,Cp+1.Ik is defined as the minimum numerical of the wij weight along a single path between concepts ith to jth. At the same time, the total causal effect T(Ci, Cj) (equation (4)) is the maximum of the indirect effect of concept Ci on concept Cj:(4)TCi,Cj=maxIkCp,Cp+1.According to Bevilacqua et al. [7], equation (3) can be described using the “weak ring in the chain” metaphor. Indeed, it is necessary for the identification of ATO identify concept concatenation as a chain where the weight wij is the hardiness of each chain ring. In the presence of a weak ring into the chain, it is not possible to consider it as a “resistant chain,” and its total hardiness is quantified with the hardiness of the weak ring. Therefore, once derived the value of hardiness is derived from by equation (3), and equation (4) allows defining the more resistant chain to be defined. Finally, the chain hardiness highlights the relevance of the first concept in the concatenations affecting the top event.In the proposed approach, the concepts of the FCM represent the working conditions of the component to be analysed, the sensor signals installed on the components, and the HI of the component. The FCM takes the advantage of the situation to identify the relationships among all the involved concepts in a matrix form to be used to calculate the health indicator for the RUL prediction.


#### 3.1.2. The GWO Algorithm

Mirjalili et al. [69] introduced the GWO, which mimics the hierarchy of leadership and the mechanism for hunting grey wolf packs in the wild. The algorithm divides the agents (grey wolves) into four different hierarchical categories called alpha (α), beta (ß), delta (δ), and omega (ω), in the descending order.Each hierarchy has different roles to find solutions, which in this case correspond to the prey. The leaders of the packs are the wolves called alphas. The alpha is primarily responsible for decisions about hunting, where to sleep, and so on. The alpha wolf is the dominant one, and the pack must follow his orders. The beta wolves identify the second level in the hierarchy. They are subordinate wolves that help the alpha in decision-making or other pack activities. Moreover, a beta wolf must not only respect the alpha but also command other lower level wolves.The lower level grey wolf is the omega. Theω plays the role of scapegoat, and it helps to satisfy the entire pack and maintain the dominant structure. Omega wolves must always submit themselves to all other dominant wolves. It may seem that the omega is not an essential individual in the pack, but it is also true that the entire pack faces internal struggles and problems if the omega is lost. If a wolf is not an alpha, beta, or omega, it is called a subordinate (or delta in some references). Delta wolves must submit themselves to alphas and betas, but dominate omegas. Scouts and hunters, for example, belong to this category. They are responsible for guarding the boundaries of the territory and warning the pack in case of danger. Hunters help the alphas and betas to hunt prey and provide food.To mathematically model the social hierarchy of wolves in the GWO design,α is therefore considered the most suitable (optimal) solution. Consequently, the second- and third-best solutions are named β and δ, respectively. The remaining candidate solutions are the ω ones.In the GWO algorithm,α, β, and δ wolves impose the rules of hunting and the ω ones follow them. In particular, the hunt is composed of three main phases such as (i) searching and chasing prey, (ii) surrounding and harassing the victim until it stops moving, and, finally, (iii) attacking the prey.After spotting the possible prey, the wolves begin to surround it and then move on to the attack. Equations (5) and (6) model mathematically encirclement behaviour:(5)D=C⋅xpt−xt,(6)xt+1=xpt−A⋅D,where D represents the difference between the position of the prey and the predator, t denotes the current iteration, xp specifies the location of the victim, and x indicates the wolf location. Equations (7) and (8) allow one to calculate the A and C values:(7)A=2⋅a⋅r1−a,(8)C=2⋅r2,where the components of a linearly decrease from 2 to 0 during each iteration and r1 and r2 are random arrays with ranging in [0 1], and they allow wolves to reach any position between the points, as illustrated in Figure 2.Figure 2
The pack hunting scheme [69].
(a)(b)As shown in Figure2(a), a wolf in position (X, Y) can update its location according to the prey’s position X∗,Y∗, and the same consideration is possible in 3D space (Figure 2(b)), or in n dimension space.It is assumed that alpha (best candidate solution), beta, and delta have a better knowledge of the potential position of the prey to simulate the hunting behaviour of wolves mathematically. Therefore, the first three best solutions are considered, and the other search agents (omega wolves) are obliged to update their positions according to the location of the best search agent [70].As mentioned above, wolves end the hunt by attacking their prey when it stops moving. If |A| < 1, the wolves begin the attack phase by moving towards the victim. Wolves look for prey mainly based on alpha, beta, and delta positions. In this phase of research (exploration), the wolves move away from each other to identify the different places of the prey (solutions). The vector A assumes values higher than one or less than −1 and forces the research agent to diverge from the victim. This emphasises the exploration and allows the GWO algorithm to search globally to find better prey. Thus, onceα, β, and δ wolves are identified, all of the members’ pack positions are updated according to(9)Dα=C1⋅xαt−xt,Dβ=C2⋅xβt−xt,Dδ=C3⋅xδt−xt,(10)x1t+1=xαt−A1⋅Dα,x2t+1=xβt−A2⋅Dβ,x3t+1=xδt−A3⋅Dδ,(11)xt+1=x1t+1+x2t+1+x2t+13.Figure3 describes the step to implement the GWO according to the mentioned equations.Figure 3
The Grey Wolf flowchart.The GWO has the advantage of having few parameters to initialise and be a flexible algorithm, so it can adapt to various practical engineering problems. Indeed, only the number of wolves in the pack (nPop) and the maximum number of iteration (MaxIt) must be initialised. In Figure 3, Iter is the current iteration. Moreover, the GWO can be easily implemented, and thanks to its hierarchical structure, which guarantees high accuracy in the solution.Although recently introduced, the GWO has been used in various fields of application. Das et al. [71] have tested the GWO to optimise the parameters of a PID controller used for speed control of a DC motor system. Komaki and Kayvanfar [72] proposed the application of GWO to program the optimal machining and assembly sequence to minimise the completion time. The results obtained with this algorithm were then compared with other methods. This comparison revealed that the GWO provided better performance. Nguyen et al. [73] used a multiobjective GWO to solve the problem of node location in a wireless sensor network. Various constraints were considered in the localisation model, including the limitation of spatial distance between nodes and the restriction of topology. The results of the simulations show significant improvements in terms of localisation accuracy and rate of convergence to the optimal solution, compared to those obtained with other methods. Song et al. [74] used GWO to estimate the parameters of Rayleigh Waves (a type of elastic surface wave).However, the research and development activities for this algorithm are still at an early stage [75]. As previously stated, the GWO has a strong exploration capacity, which can avoid convergence in excellent premises. This feature may lead the algorithm to slow convergence and indeed led us to try GWO to define the wij values of the FCM.


#### 3.1.3. The K-Nearest Neighbors Algorithm

An in-depth analysis of the k-nearest neighbors (KNN) algorithm allows underlining as it is simple and easy-to-implement supervised machine learning algorithm used to solve both classification and regression problems. Its functioning is based on the similarity of the characteristics: the closer an instance is to a data point, the more KNN will consider them similar [76].Once the HIs have been defined for each training unit, they can be used as models representing the degradation profile, from normal functioning to disruption. At this point, a set of modelsMi (with i = 1 to the number of items composing the training dataset) is available and usable to predict the RUL. Therefore, to find the most similar element, it is necessary to measure the distance between the model Mi and Y = y1, y2, ..., yr, which represents the HI of the test unit obtained through consecutive observations. The distance is calculated by the Euclidean distance (depending on the problem under examination) or by the mean value of the absolute residual (used in the proposed approach), as described by equations (11) and (12). Thus, the smaller the distance, the greater the similarity between the data point and the instance to be predicted:(12)dxi,xl=xi1−xl12+xi2−xl22+…+xip−xlp2,(13)dxi,xl=meanyi−yl,where yi is ith training model and yl is the l th testing one, and each of them is composed of xi and xl components.Then, the calculated distances are used as the argument to evaluate the similarity weight,swi,l, between the testing HI and all of the training ones considering(14)swi,l=exp−dxi,xl2.Once obtaining the similarity weights among the testing unit and training ones, it is possible to rank them in descending order and identify the number of similar unit SU as described in(15)SU=mink,N,where N is the number of training units and k is an initial fixed value. In particular, when k is small, it is limiting the prediction region, forcing the classifier to be “more blind” than the general distribution. On the contrary, a large k reduces the impact of the variance caused by a random error but runs the risk of ignoring small details that might be relevant. For the proposed approach, initially, k is fixed equal to 50. Thus, having similar units, it is possible, considering the relative end dates updated to the number of test samples, to fit a Weibull distribution to find the RUL.


#### 3.1.4. The Procedure for HI Assessment

Figure4 shows the framework used to esteem the aircraft engine HIs and, subsequently, to predict the RUL of the engine. The proposed algorithm can be classified within the condition monitoring techniques. It consists of a general framework and can be applied to any equipment. The dataset, both training and testing, are composed of the sensor readings of the considered items.Figure 4
The proposed framework approach.In the “Time Indicator Modelling” phase, a lifetime indicator (LTI) is defined. The sample number of each piece of equipment (that corresponds to the number of rows of the dataset) represents its life duration.The main idea of LTI is to model a degradation profile considering that, at the beginning of sampling, an item has the maximum reliability value (equal to 1), and when the disruption occurs, the item reliability has a minimum amount (equal to 0). The first value of LTI is equal to 1, and the last one is equal to zero, according to(16)TIm′=DURm−1DURm−2…0,(17)TImt'′=TImt'DURm,(18)LTImt=TImt′​′+1−TImt=1′​′,where DURm is the dataset length for the mth equipment. Each element of TIm indicates the remaining cycle times to the relative disruption at time t (then normalised in TIm”). Hence, each value of LTIm decreases from 1 to 0. LTIm represents the parameter to be esteemed and used in the algorithm for the HI estimation. Table 2 shows an example of LTI calculation.Table 2
An example of equations (16)–(18) application.
Steps descriptionEquation applicationSuppose machine 1 breaks after 10 sampling cycles; it follows thatDUR1=10, according to equation (16)TI1′=9876543210According to equation (17), the relative normalised valueTI1'′=0.90.80.70.60.50.40.30.20.10According to equation (18), the lifetime indicatorLTI1=10.90.80.70.60.50.40.30.20.1The “FCM Modelling” phase is the core of the proposed approach to identify the HI. Figure5 describes the iterative phase for the HI calculation reviewing the general GWO algorithm shown in Figure 3. In particular, the GWO algorithm is used for defining the weight of the relation between concept ith and jth (wij values) of the FCM matrix. The concepts of the FCM represent the working conditions of the equipment to be analysed, the sensor signals installed on the equipment, and the last concept which is the HI. In particular, since the purpose of the approach is the HI estimation using the FCM theory, the number of concepts (NC) to be used is equal to the number of reduced dataset variables’ number plus the HI (the algorithm output). This means that if the reduced dataset variables number is n, NC = n + 1.Figure 5
Framework of the proposed methodology.The iterative phase, shown in Figure5, is executed for each equipment belonging to the training dataset. In each iteration, the final α wolf position is assumed as the temporary FCMj and used as initial FCM for the next one. When the terminal equipment has been analysed, the relative FCM is considered as the optimal solution.The GWO algorithm is used to define thewij values of FCM. These values randomly in the range [−1 1] or [0 1], as required by the FCM theory, for each pack member. Then, equation (1) calculates the relative cost. By analysing Figure 5, it is possible to highlight how the lifetime indicator and the reduced dataset are used as input for the fitness function calculation.Assuming to have a reduced training dataset related toM equipment, with n main variables, it is divided in M reduced subdataset, each one related to a specific equipment. Thus, a reduced subdataset RDm (the sensors’ readings of equipment m with m = 1, 2, …, M) is available with the relative LTIm, in the form expressed in(19)RDm=Var1,1m…Varn,1m⋮⋱⋮Var1,fm…Varn,fm,LTIm=lti1…ltif.The termf identifies the instant in which the fault occurred; n represents the progressive number of the main reduced variables, and m is the number of considered equipment.The pack members’ position, obtained through the GWO application for a specific device, is given in the form of(20)FCMp,miter=w1,1…w1,NC⋮⋱⋮wNC,1…wNC,NC,where p is the pth member of the pack (p = 1, 2, …nPop), iter is the current iteration with iter ≤ MaxIt, and NC is the concepts number for the FCM algorithm defined before. Thus, FCMp,miter is the relative position of the pth pack member at the current iteration for the mth equipment.By analysing Figure5, it is possible to highlight the presence of two main “for loop.” The first one (external) is referred to the number of available equipment (M) and the second one (internal) to the maximum iteration number for GWO (MaxIt).At each iteration of the external loop (the iteration is equal to the equipment number in the dataset), the LTI related to a specific item is used as the benchmark for the positional cost calculation (if the iteration is equal to one, the LTI1 is examined). This means that the external loop has the objective of identifying the best FCMm for the mth equipment. In the internal loop, for the GWO application, all the pack members take a position within the domain space, updating it at each inner iteration. At the end of the inner loop, the best position for the considered item is identified (FCMm). The obtained FCMm is used as the initial position for the FCMm+1 identification (as long as m < M) to improve its accuracy.In the algorithm initialisation,FCM10 can be defined randomly if there is no knowledge of the involved equipment or a panel of experts cannot be established to model it, as described for the classical FCM design approach, according to the experience of each professional involved.As mentioned before, equation (1) evaluates the fitness cost value for each FCMp,miter. This step is the most critical in the whole algorithm, as highlighted by Mazzuto and Stylios [77]. Indeed, to calculate the cost connected to FCMp,miterr, it is necessary to consider all of the samples which make up RDm and LTIm. More accurately, if RDm has f samples (as described above) as well as LTIm, equation (1) has to be applied f times. Besides, since the number of iteration (k) in equation (1) depends on the function convergence or the fixed amount of repetition (FCMiter), the fitness cost evaluation requires an iteration number equal to f⋅FCMiter. Considering nPop wolves and a maximum iteration number for the GWO (MaxIt), the total number of iterations for the identification of the optimal FCM is equal to f⋅FCMiter⋅nPop⋅MaxIt.Considering thepth pack member, once its position is defined (FCMp,miter), its positional cost has to be calculated according to equation (1) for the best position identification.The main idea is to consider each sample in the reduced dataset for a specific item (RDm) as the initial array A0 at the specified iteration according to(21)Am0,j=Var1,jm…Varn,jm0,4j=1,2…f,m=1,2…M.The null value is because the output is considered within the set of FCM concepts, but it is the variable that needs to be taken into consideration.When the application of equation (1) reaches the convergent, Am∗,j will be described by(22)Am∗,j=Var1,jm^…Varn,jm^HIm∗,j,j=1,2…f,m=1,2…M,where Var1,jm^ is the convergent value for the variable Var1,jm and HIm∗,j is the esteemed output for the sample at time j for the mth engine.Once all of the samples inRDm have been processed, a final output array HIp,m∗ for the pth pack member, equation (23), will be available:(23)Hp,m∗=HIm∗,1HIm∗,2…HIm∗,f,f,m=1,2…M.Thus, the esteemed outputHIp,m∗ and the connected LTIm can be used to calculate the fitness cost value (Cp,miter), for the pth pack member and the mth item and at iteration iter, using the root mean squared error formula, as shown in(24)Cp,miter=∑j=1fLTImj−Hp,m∗,j2f.The root mean square error has been chosen because it describes efficiently how concentrated the data is around the line of best fit [78].Finally, once the optimal FCM to be used for the HI identification phases is identified, it is used to calculate the HIs in the “Training Health Indicator definitions” and “Testing Health Indicator definitions,” respectively.


## 3.1. The HI Modelling

The proposed approach for the HI modelling is performed by the combined use of FCM and a Swarm Intelligence algorithm such as the GWO. Before describing the proposed approach, the FCM and GWO theories are briefly described in the following two sections.


### 3.1.1. The FCMs’ Modelling

A cognitive map (CM) can be thought of as a concept map reflecting mental processing, comprised of collected information and several cognitive abstractions, individually filtered, about regarding physical phenomena and experiences [56]. Cognitive maps are visual representations of an individual’s mental model constructs, analogous to concept maps for representing human reasoning and knowledge or beliefs [7]. Thus, a generic problem is considered, and an expert panel of experts is formed for its in-depth analysis, since different individuals may face the same question differently. According to their area of expertise through fuzzy logic, they model collective FCM identifying concepts and relationships about regarding the considered problem. In particular, concepts, in number of N, are the FCM key elements that stand for the main characteristics of the abstract mental model for whichever complex system [57]. Once concepts are identified, experts are asked to assign a numerical value wij (the weight of the relation between concept ith and jth) for the W matrix, which represents the influence of concept Ci on concept Cj. According to equation (1), wij ranges in [–1, 1]. Specifically, wij = 0 indicates no causality between concepts, wij > 0 indicates causal Cj increases as Ci increases (or Cj decreases as Ci decreases), and wij < 0 shows causal decrease or negative causality (Cj decreases as Ci increases or Cj increases as Ci decreases):(1)FCM=w1,1…w1,N⋮⋱⋮wN,1…wN,N.Although many studies exist concerning the dynamical representation of an FCM, generally, the experts’ opinion aggregation of expert opinions for the collective weight matrix modelling is performed using the SUM method [58]. Then, overall linguistic weight is evaluated using the centre of gravity (COG) defuzzification method [59]. Some examples are presented by Bevilacqua et al. [7, 60, 61] and Stylios et al. [62] where a unique credibility value is assigned to each expert and a threshold function is used in the aggregation. On the contrary, a modification of the approach mentioned above has been provided by Stylios and Groumpos [63] and Stylios and Groumpos [64], introducing a corrective factor for the experts’ credibility evaluation. However, this approach does not take into consideration the fact that, in a complex multidisciplinary problem, most experts have in-depth knowledge of only parts of the problem and not the entire issue [65].Once the total weights’ matrix,W, has been designed, it is possible to analyse the system behaviour through simulations. Thus, if Ai defines the instantaneous value of concept Ci, its evolution over time can be evaluated computing the influence of the related concepts Cj on the specific concept Ci according to(2)Aik+1=fAik+∑j≠ij=1nAjkwi,j,where Aik + 1 is the value of concept Ci at simulation step k + 1 and Ajk is the value of concept Cj at simulation step k. Also, wij is the weight of the interconnection from concept Cj to concept Ci and f is an appropriate threshold function used to force the concept value to be monotonically mapped into a normalised range [66]. Other equations can be used in place of equation (2) as suggested by Mazzuto et al. [67] and Osoba and Kosko [68].An important topic in the FCM analysis is the indirect and total causal effect evaluation (Axelrod, 1976), whose knowledge allows an in-depth map analysis. The indirect effectIk of Ci concept on Cj concept can be defined as shown in(3)IkCi,Cj=minwCp,Cp+1.Ik is defined as the minimum numerical of the wij weight along a single path between concepts ith to jth. At the same time, the total causal effect T(Ci, Cj) (equation (4)) is the maximum of the indirect effect of concept Ci on concept Cj:(4)TCi,Cj=maxIkCp,Cp+1.According to Bevilacqua et al. [7], equation (3) can be described using the “weak ring in the chain” metaphor. Indeed, it is necessary for the identification of ATO identify concept concatenation as a chain where the weight wij is the hardiness of each chain ring. In the presence of a weak ring into the chain, it is not possible to consider it as a “resistant chain,” and its total hardiness is quantified with the hardiness of the weak ring. Therefore, once derived the value of hardiness is derived from by equation (3), and equation (4) allows defining the more resistant chain to be defined. Finally, the chain hardiness highlights the relevance of the first concept in the concatenations affecting the top event.In the proposed approach, the concepts of the FCM represent the working conditions of the component to be analysed, the sensor signals installed on the components, and the HI of the component. The FCM takes the advantage of the situation to identify the relationships among all the involved concepts in a matrix form to be used to calculate the health indicator for the RUL prediction.


### 3.1.2. The GWO Algorithm

Mirjalili et al. [69] introduced the GWO, which mimics the hierarchy of leadership and the mechanism for hunting grey wolf packs in the wild. The algorithm divides the agents (grey wolves) into four different hierarchical categories called alpha (α), beta (ß), delta (δ), and omega (ω), in the descending order.Each hierarchy has different roles to find solutions, which in this case correspond to the prey. The leaders of the packs are the wolves called alphas. The alpha is primarily responsible for decisions about hunting, where to sleep, and so on. The alpha wolf is the dominant one, and the pack must follow his orders. The beta wolves identify the second level in the hierarchy. They are subordinate wolves that help the alpha in decision-making or other pack activities. Moreover, a beta wolf must not only respect the alpha but also command other lower level wolves.The lower level grey wolf is the omega. Theω plays the role of scapegoat, and it helps to satisfy the entire pack and maintain the dominant structure. Omega wolves must always submit themselves to all other dominant wolves. It may seem that the omega is not an essential individual in the pack, but it is also true that the entire pack faces internal struggles and problems if the omega is lost. If a wolf is not an alpha, beta, or omega, it is called a subordinate (or delta in some references). Delta wolves must submit themselves to alphas and betas, but dominate omegas. Scouts and hunters, for example, belong to this category. They are responsible for guarding the boundaries of the territory and warning the pack in case of danger. Hunters help the alphas and betas to hunt prey and provide food.To mathematically model the social hierarchy of wolves in the GWO design,α is therefore considered the most suitable (optimal) solution. Consequently, the second- and third-best solutions are named β and δ, respectively. The remaining candidate solutions are the ω ones.In the GWO algorithm,α, β, and δ wolves impose the rules of hunting and the ω ones follow them. In particular, the hunt is composed of three main phases such as (i) searching and chasing prey, (ii) surrounding and harassing the victim until it stops moving, and, finally, (iii) attacking the prey.After spotting the possible prey, the wolves begin to surround it and then move on to the attack. Equations (5) and (6) model mathematically encirclement behaviour:(5)D=C⋅xpt−xt,(6)xt+1=xpt−A⋅D,where D represents the difference between the position of the prey and the predator, t denotes the current iteration, xp specifies the location of the victim, and x indicates the wolf location. Equations (7) and (8) allow one to calculate the A and C values:(7)A=2⋅a⋅r1−a,(8)C=2⋅r2,where the components of a linearly decrease from 2 to 0 during each iteration and r1 and r2 are random arrays with ranging in [0 1], and they allow wolves to reach any position between the points, as illustrated in Figure 2.Figure 2
The pack hunting scheme [69].
(a)(b)As shown in Figure2(a), a wolf in position (X, Y) can update its location according to the prey’s position X∗,Y∗, and the same consideration is possible in 3D space (Figure 2(b)), or in n dimension space.It is assumed that alpha (best candidate solution), beta, and delta have a better knowledge of the potential position of the prey to simulate the hunting behaviour of wolves mathematically. Therefore, the first three best solutions are considered, and the other search agents (omega wolves) are obliged to update their positions according to the location of the best search agent [70].As mentioned above, wolves end the hunt by attacking their prey when it stops moving. If |A| < 1, the wolves begin the attack phase by moving towards the victim. Wolves look for prey mainly based on alpha, beta, and delta positions. In this phase of research (exploration), the wolves move away from each other to identify the different places of the prey (solutions). The vector A assumes values higher than one or less than −1 and forces the research agent to diverge from the victim. This emphasises the exploration and allows the GWO algorithm to search globally to find better prey. Thus, onceα, β, and δ wolves are identified, all of the members’ pack positions are updated according to(9)Dα=C1⋅xαt−xt,Dβ=C2⋅xβt−xt,Dδ=C3⋅xδt−xt,(10)x1t+1=xαt−A1⋅Dα,x2t+1=xβt−A2⋅Dβ,x3t+1=xδt−A3⋅Dδ,(11)xt+1=x1t+1+x2t+1+x2t+13.Figure3 describes the step to implement the GWO according to the mentioned equations.Figure 3
The Grey Wolf flowchart.The GWO has the advantage of having few parameters to initialise and be a flexible algorithm, so it can adapt to various practical engineering problems. Indeed, only the number of wolves in the pack (nPop) and the maximum number of iteration (MaxIt) must be initialised. In Figure 3, Iter is the current iteration. Moreover, the GWO can be easily implemented, and thanks to its hierarchical structure, which guarantees high accuracy in the solution.Although recently introduced, the GWO has been used in various fields of application. Das et al. [71] have tested the GWO to optimise the parameters of a PID controller used for speed control of a DC motor system. Komaki and Kayvanfar [72] proposed the application of GWO to program the optimal machining and assembly sequence to minimise the completion time. The results obtained with this algorithm were then compared with other methods. This comparison revealed that the GWO provided better performance. Nguyen et al. [73] used a multiobjective GWO to solve the problem of node location in a wireless sensor network. Various constraints were considered in the localisation model, including the limitation of spatial distance between nodes and the restriction of topology. The results of the simulations show significant improvements in terms of localisation accuracy and rate of convergence to the optimal solution, compared to those obtained with other methods. Song et al. [74] used GWO to estimate the parameters of Rayleigh Waves (a type of elastic surface wave).However, the research and development activities for this algorithm are still at an early stage [75]. As previously stated, the GWO has a strong exploration capacity, which can avoid convergence in excellent premises. This feature may lead the algorithm to slow convergence and indeed led us to try GWO to define the wij values of the FCM.


### 3.1.3. The K-Nearest Neighbors Algorithm

An in-depth analysis of the k-nearest neighbors (KNN) algorithm allows underlining as it is simple and easy-to-implement supervised machine learning algorithm used to solve both classification and regression problems. Its functioning is based on the similarity of the characteristics: the closer an instance is to a data point, the more KNN will consider them similar [76].Once the HIs have been defined for each training unit, they can be used as models representing the degradation profile, from normal functioning to disruption. At this point, a set of modelsMi (with i = 1 to the number of items composing the training dataset) is available and usable to predict the RUL. Therefore, to find the most similar element, it is necessary to measure the distance between the model Mi and Y = y1, y2, ..., yr, which represents the HI of the test unit obtained through consecutive observations. The distance is calculated by the Euclidean distance (depending on the problem under examination) or by the mean value of the absolute residual (used in the proposed approach), as described by equations (11) and (12). Thus, the smaller the distance, the greater the similarity between the data point and the instance to be predicted:(12)dxi,xl=xi1−xl12+xi2−xl22+…+xip−xlp2,(13)dxi,xl=meanyi−yl,where yi is ith training model and yl is the l th testing one, and each of them is composed of xi and xl components.Then, the calculated distances are used as the argument to evaluate the similarity weight,swi,l, between the testing HI and all of the training ones considering(14)swi,l=exp−dxi,xl2.Once obtaining the similarity weights among the testing unit and training ones, it is possible to rank them in descending order and identify the number of similar unit SU as described in(15)SU=mink,N,where N is the number of training units and k is an initial fixed value. In particular, when k is small, it is limiting the prediction region, forcing the classifier to be “more blind” than the general distribution. On the contrary, a large k reduces the impact of the variance caused by a random error but runs the risk of ignoring small details that might be relevant. For the proposed approach, initially, k is fixed equal to 50. Thus, having similar units, it is possible, considering the relative end dates updated to the number of test samples, to fit a Weibull distribution to find the RUL.


### 3.1.4. The Procedure for HI Assessment

Figure4 shows the framework used to esteem the aircraft engine HIs and, subsequently, to predict the RUL of the engine. The proposed algorithm can be classified within the condition monitoring techniques. It consists of a general framework and can be applied to any equipment. The dataset, both training and testing, are composed of the sensor readings of the considered items.Figure 4
The proposed framework approach.In the “Time Indicator Modelling” phase, a lifetime indicator (LTI) is defined. The sample number of each piece of equipment (that corresponds to the number of rows of the dataset) represents its life duration.The main idea of LTI is to model a degradation profile considering that, at the beginning of sampling, an item has the maximum reliability value (equal to 1), and when the disruption occurs, the item reliability has a minimum amount (equal to 0). The first value of LTI is equal to 1, and the last one is equal to zero, according to(16)TIm′=DURm−1DURm−2…0,(17)TImt'′=TImt'DURm,(18)LTImt=TImt′​′+1−TImt=1′​′,where DURm is the dataset length for the mth equipment. Each element of TIm indicates the remaining cycle times to the relative disruption at time t (then normalised in TIm”). Hence, each value of LTIm decreases from 1 to 0. LTIm represents the parameter to be esteemed and used in the algorithm for the HI estimation. Table 2 shows an example of LTI calculation.Table 2
An example of equations (16)–(18) application.
Steps descriptionEquation applicationSuppose machine 1 breaks after 10 sampling cycles; it follows thatDUR1=10, according to equation (16)TI1′=9876543210According to equation (17), the relative normalised valueTI1'′=0.90.80.70.60.50.40.30.20.10According to equation (18), the lifetime indicatorLTI1=10.90.80.70.60.50.40.30.20.1The “FCM Modelling” phase is the core of the proposed approach to identify the HI. Figure5 describes the iterative phase for the HI calculation reviewing the general GWO algorithm shown in Figure 3. In particular, the GWO algorithm is used for defining the weight of the relation between concept ith and jth (wij values) of the FCM matrix. The concepts of the FCM represent the working conditions of the equipment to be analysed, the sensor signals installed on the equipment, and the last concept which is the HI. In particular, since the purpose of the approach is the HI estimation using the FCM theory, the number of concepts (NC) to be used is equal to the number of reduced dataset variables’ number plus the HI (the algorithm output). This means that if the reduced dataset variables number is n, NC = n + 1.Figure 5
Framework of the proposed methodology.The iterative phase, shown in Figure5, is executed for each equipment belonging to the training dataset. In each iteration, the final α wolf position is assumed as the temporary FCMj and used as initial FCM for the next one. When the terminal equipment has been analysed, the relative FCM is considered as the optimal solution.The GWO algorithm is used to define thewij values of FCM. These values randomly in the range [−1 1] or [0 1], as required by the FCM theory, for each pack member. Then, equation (1) calculates the relative cost. By analysing Figure 5, it is possible to highlight how the lifetime indicator and the reduced dataset are used as input for the fitness function calculation.Assuming to have a reduced training dataset related toM equipment, with n main variables, it is divided in M reduced subdataset, each one related to a specific equipment. Thus, a reduced subdataset RDm (the sensors’ readings of equipment m with m = 1, 2, …, M) is available with the relative LTIm, in the form expressed in(19)RDm=Var1,1m…Varn,1m⋮⋱⋮Var1,fm…Varn,fm,LTIm=lti1…ltif.The termf identifies the instant in which the fault occurred; n represents the progressive number of the main reduced variables, and m is the number of considered equipment.The pack members’ position, obtained through the GWO application for a specific device, is given in the form of(20)FCMp,miter=w1,1…w1,NC⋮⋱⋮wNC,1…wNC,NC,where p is the pth member of the pack (p = 1, 2, …nPop), iter is the current iteration with iter ≤ MaxIt, and NC is the concepts number for the FCM algorithm defined before. Thus, FCMp,miter is the relative position of the pth pack member at the current iteration for the mth equipment.By analysing Figure5, it is possible to highlight the presence of two main “for loop.” The first one (external) is referred to the number of available equipment (M) and the second one (internal) to the maximum iteration number for GWO (MaxIt).At each iteration of the external loop (the iteration is equal to the equipment number in the dataset), the LTI related to a specific item is used as the benchmark for the positional cost calculation (if the iteration is equal to one, the LTI1 is examined). This means that the external loop has the objective of identifying the best FCMm for the mth equipment. In the internal loop, for the GWO application, all the pack members take a position within the domain space, updating it at each inner iteration. At the end of the inner loop, the best position for the considered item is identified (FCMm). The obtained FCMm is used as the initial position for the FCMm+1 identification (as long as m < M) to improve its accuracy.In the algorithm initialisation,FCM10 can be defined randomly if there is no knowledge of the involved equipment or a panel of experts cannot be established to model it, as described for the classical FCM design approach, according to the experience of each professional involved.As mentioned before, equation (1) evaluates the fitness cost value for each FCMp,miter. This step is the most critical in the whole algorithm, as highlighted by Mazzuto and Stylios [77]. Indeed, to calculate the cost connected to FCMp,miterr, it is necessary to consider all of the samples which make up RDm and LTIm. More accurately, if RDm has f samples (as described above) as well as LTIm, equation (1) has to be applied f times. Besides, since the number of iteration (k) in equation (1) depends on the function convergence or the fixed amount of repetition (FCMiter), the fitness cost evaluation requires an iteration number equal to f⋅FCMiter. Considering nPop wolves and a maximum iteration number for the GWO (MaxIt), the total number of iterations for the identification of the optimal FCM is equal to f⋅FCMiter⋅nPop⋅MaxIt.Considering thepth pack member, once its position is defined (FCMp,miter), its positional cost has to be calculated according to equation (1) for the best position identification.The main idea is to consider each sample in the reduced dataset for a specific item (RDm) as the initial array A0 at the specified iteration according to(21)Am0,j=Var1,jm…Varn,jm0,4j=1,2…f,m=1,2…M.The null value is because the output is considered within the set of FCM concepts, but it is the variable that needs to be taken into consideration.When the application of equation (1) reaches the convergent, Am∗,j will be described by(22)Am∗,j=Var1,jm^…Varn,jm^HIm∗,j,j=1,2…f,m=1,2…M,where Var1,jm^ is the convergent value for the variable Var1,jm and HIm∗,j is the esteemed output for the sample at time j for the mth engine.Once all of the samples inRDm have been processed, a final output array HIp,m∗ for the pth pack member, equation (23), will be available:(23)Hp,m∗=HIm∗,1HIm∗,2…HIm∗,f,f,m=1,2…M.Thus, the esteemed outputHIp,m∗ and the connected LTIm can be used to calculate the fitness cost value (Cp,miter), for the pth pack member and the mth item and at iteration iter, using the root mean squared error formula, as shown in(24)Cp,miter=∑j=1fLTImj−Hp,m∗,j2f.The root mean square error has been chosen because it describes efficiently how concentrated the data is around the line of best fit [78].Finally, once the optimal FCM to be used for the HI identification phases is identified, it is used to calculate the HIs in the “Training Health Indicator definitions” and “Testing Health Indicator definitions,” respectively.


## 3.1.1. The FCMs’ Modelling

A cognitive map (CM) can be thought of as a concept map reflecting mental processing, comprised of collected information and several cognitive abstractions, individually filtered, about regarding physical phenomena and experiences [56]. Cognitive maps are visual representations of an individual’s mental model constructs, analogous to concept maps for representing human reasoning and knowledge or beliefs [7]. Thus, a generic problem is considered, and an expert panel of experts is formed for its in-depth analysis, since different individuals may face the same question differently. According to their area of expertise through fuzzy logic, they model collective FCM identifying concepts and relationships about regarding the considered problem. In particular, concepts, in number of N, are the FCM key elements that stand for the main characteristics of the abstract mental model for whichever complex system [57]. Once concepts are identified, experts are asked to assign a numerical value wij (the weight of the relation between concept ith and jth) for the W matrix, which represents the influence of concept Ci on concept Cj. According to equation (1), wij ranges in [–1, 1]. Specifically, wij = 0 indicates no causality between concepts, wij > 0 indicates causal Cj increases as Ci increases (or Cj decreases as Ci decreases), and wij < 0 shows causal decrease or negative causality (Cj decreases as Ci increases or Cj increases as Ci decreases):(1)FCM=w1,1…w1,N⋮⋱⋮wN,1…wN,N.Although many studies exist concerning the dynamical representation of an FCM, generally, the experts’ opinion aggregation of expert opinions for the collective weight matrix modelling is performed using the SUM method [58]. Then, overall linguistic weight is evaluated using the centre of gravity (COG) defuzzification method [59]. Some examples are presented by Bevilacqua et al. [7, 60, 61] and Stylios et al. [62] where a unique credibility value is assigned to each expert and a threshold function is used in the aggregation. On the contrary, a modification of the approach mentioned above has been provided by Stylios and Groumpos [63] and Stylios and Groumpos [64], introducing a corrective factor for the experts’ credibility evaluation. However, this approach does not take into consideration the fact that, in a complex multidisciplinary problem, most experts have in-depth knowledge of only parts of the problem and not the entire issue [65].Once the total weights’ matrix,W, has been designed, it is possible to analyse the system behaviour through simulations. Thus, if Ai defines the instantaneous value of concept Ci, its evolution over time can be evaluated computing the influence of the related concepts Cj on the specific concept Ci according to(2)Aik+1=fAik+∑j≠ij=1nAjkwi,j,where Aik + 1 is the value of concept Ci at simulation step k + 1 and Ajk is the value of concept Cj at simulation step k. Also, wij is the weight of the interconnection from concept Cj to concept Ci and f is an appropriate threshold function used to force the concept value to be monotonically mapped into a normalised range [66]. Other equations can be used in place of equation (2) as suggested by Mazzuto et al. [67] and Osoba and Kosko [68].An important topic in the FCM analysis is the indirect and total causal effect evaluation (Axelrod, 1976), whose knowledge allows an in-depth map analysis. The indirect effectIk of Ci concept on Cj concept can be defined as shown in(3)IkCi,Cj=minwCp,Cp+1.Ik is defined as the minimum numerical of the wij weight along a single path between concepts ith to jth. At the same time, the total causal effect T(Ci, Cj) (equation (4)) is the maximum of the indirect effect of concept Ci on concept Cj:(4)TCi,Cj=maxIkCp,Cp+1.According to Bevilacqua et al. [7], equation (3) can be described using the “weak ring in the chain” metaphor. Indeed, it is necessary for the identification of ATO identify concept concatenation as a chain where the weight wij is the hardiness of each chain ring. In the presence of a weak ring into the chain, it is not possible to consider it as a “resistant chain,” and its total hardiness is quantified with the hardiness of the weak ring. Therefore, once derived the value of hardiness is derived from by equation (3), and equation (4) allows defining the more resistant chain to be defined. Finally, the chain hardiness highlights the relevance of the first concept in the concatenations affecting the top event.In the proposed approach, the concepts of the FCM represent the working conditions of the component to be analysed, the sensor signals installed on the components, and the HI of the component. The FCM takes the advantage of the situation to identify the relationships among all the involved concepts in a matrix form to be used to calculate the health indicator for the RUL prediction.


## 3.1.2. The GWO Algorithm

Mirjalili et al. [69] introduced the GWO, which mimics the hierarchy of leadership and the mechanism for hunting grey wolf packs in the wild. The algorithm divides the agents (grey wolves) into four different hierarchical categories called alpha (α), beta (ß), delta (δ), and omega (ω), in the descending order.Each hierarchy has different roles to find solutions, which in this case correspond to the prey. The leaders of the packs are the wolves called alphas. The alpha is primarily responsible for decisions about hunting, where to sleep, and so on. The alpha wolf is the dominant one, and the pack must follow his orders. The beta wolves identify the second level in the hierarchy. They are subordinate wolves that help the alpha in decision-making or other pack activities. Moreover, a beta wolf must not only respect the alpha but also command other lower level wolves.The lower level grey wolf is the omega. Theω plays the role of scapegoat, and it helps to satisfy the entire pack and maintain the dominant structure. Omega wolves must always submit themselves to all other dominant wolves. It may seem that the omega is not an essential individual in the pack, but it is also true that the entire pack faces internal struggles and problems if the omega is lost. If a wolf is not an alpha, beta, or omega, it is called a subordinate (or delta in some references). Delta wolves must submit themselves to alphas and betas, but dominate omegas. Scouts and hunters, for example, belong to this category. They are responsible for guarding the boundaries of the territory and warning the pack in case of danger. Hunters help the alphas and betas to hunt prey and provide food.To mathematically model the social hierarchy of wolves in the GWO design,α is therefore considered the most suitable (optimal) solution. Consequently, the second- and third-best solutions are named β and δ, respectively. The remaining candidate solutions are the ω ones.In the GWO algorithm,α, β, and δ wolves impose the rules of hunting and the ω ones follow them. In particular, the hunt is composed of three main phases such as (i) searching and chasing prey, (ii) surrounding and harassing the victim until it stops moving, and, finally, (iii) attacking the prey.After spotting the possible prey, the wolves begin to surround it and then move on to the attack. Equations (5) and (6) model mathematically encirclement behaviour:(5)D=C⋅xpt−xt,(6)xt+1=xpt−A⋅D,where D represents the difference between the position of the prey and the predator, t denotes the current iteration, xp specifies the location of the victim, and x indicates the wolf location. Equations (7) and (8) allow one to calculate the A and C values:(7)A=2⋅a⋅r1−a,(8)C=2⋅r2,where the components of a linearly decrease from 2 to 0 during each iteration and r1 and r2 are random arrays with ranging in [0 1], and they allow wolves to reach any position between the points, as illustrated in Figure 2.Figure 2
The pack hunting scheme [69].
(a)(b)As shown in Figure2(a), a wolf in position (X, Y) can update its location according to the prey’s position X∗,Y∗, and the same consideration is possible in 3D space (Figure 2(b)), or in n dimension space.It is assumed that alpha (best candidate solution), beta, and delta have a better knowledge of the potential position of the prey to simulate the hunting behaviour of wolves mathematically. Therefore, the first three best solutions are considered, and the other search agents (omega wolves) are obliged to update their positions according to the location of the best search agent [70].As mentioned above, wolves end the hunt by attacking their prey when it stops moving. If |A| < 1, the wolves begin the attack phase by moving towards the victim. Wolves look for prey mainly based on alpha, beta, and delta positions. In this phase of research (exploration), the wolves move away from each other to identify the different places of the prey (solutions). The vector A assumes values higher than one or less than −1 and forces the research agent to diverge from the victim. This emphasises the exploration and allows the GWO algorithm to search globally to find better prey. Thus, onceα, β, and δ wolves are identified, all of the members’ pack positions are updated according to(9)Dα=C1⋅xαt−xt,Dβ=C2⋅xβt−xt,Dδ=C3⋅xδt−xt,(10)x1t+1=xαt−A1⋅Dα,x2t+1=xβt−A2⋅Dβ,x3t+1=xδt−A3⋅Dδ,(11)xt+1=x1t+1+x2t+1+x2t+13.Figure3 describes the step to implement the GWO according to the mentioned equations.Figure 3
The Grey Wolf flowchart.The GWO has the advantage of having few parameters to initialise and be a flexible algorithm, so it can adapt to various practical engineering problems. Indeed, only the number of wolves in the pack (nPop) and the maximum number of iteration (MaxIt) must be initialised. In Figure 3, Iter is the current iteration. Moreover, the GWO can be easily implemented, and thanks to its hierarchical structure, which guarantees high accuracy in the solution.Although recently introduced, the GWO has been used in various fields of application. Das et al. [71] have tested the GWO to optimise the parameters of a PID controller used for speed control of a DC motor system. Komaki and Kayvanfar [72] proposed the application of GWO to program the optimal machining and assembly sequence to minimise the completion time. The results obtained with this algorithm were then compared with other methods. This comparison revealed that the GWO provided better performance. Nguyen et al. [73] used a multiobjective GWO to solve the problem of node location in a wireless sensor network. Various constraints were considered in the localisation model, including the limitation of spatial distance between nodes and the restriction of topology. The results of the simulations show significant improvements in terms of localisation accuracy and rate of convergence to the optimal solution, compared to those obtained with other methods. Song et al. [74] used GWO to estimate the parameters of Rayleigh Waves (a type of elastic surface wave).However, the research and development activities for this algorithm are still at an early stage [75]. As previously stated, the GWO has a strong exploration capacity, which can avoid convergence in excellent premises. This feature may lead the algorithm to slow convergence and indeed led us to try GWO to define the wij values of the FCM.


## 3.1.3. The K-Nearest Neighbors Algorithm

An in-depth analysis of the k-nearest neighbors (KNN) algorithm allows underlining as it is simple and easy-to-implement supervised machine learning algorithm used to solve both classification and regression problems. Its functioning is based on the similarity of the characteristics: the closer an instance is to a data point, the more KNN will consider them similar [76].Once the HIs have been defined for each training unit, they can be used as models representing the degradation profile, from normal functioning to disruption. At this point, a set of modelsMi (with i = 1 to the number of items composing the training dataset) is available and usable to predict the RUL. Therefore, to find the most similar element, it is necessary to measure the distance between the model Mi and Y = y1, y2, ..., yr, which represents the HI of the test unit obtained through consecutive observations. The distance is calculated by the Euclidean distance (depending on the problem under examination) or by the mean value of the absolute residual (used in the proposed approach), as described by equations (11) and (12). Thus, the smaller the distance, the greater the similarity between the data point and the instance to be predicted:(12)dxi,xl=xi1−xl12+xi2−xl22+…+xip−xlp2,(13)dxi,xl=meanyi−yl,where yi is ith training model and yl is the l th testing one, and each of them is composed of xi and xl components.Then, the calculated distances are used as the argument to evaluate the similarity weight,swi,l, between the testing HI and all of the training ones considering(14)swi,l=exp−dxi,xl2.Once obtaining the similarity weights among the testing unit and training ones, it is possible to rank them in descending order and identify the number of similar unit SU as described in(15)SU=mink,N,where N is the number of training units and k is an initial fixed value. In particular, when k is small, it is limiting the prediction region, forcing the classifier to be “more blind” than the general distribution. On the contrary, a large k reduces the impact of the variance caused by a random error but runs the risk of ignoring small details that might be relevant. For the proposed approach, initially, k is fixed equal to 50. Thus, having similar units, it is possible, considering the relative end dates updated to the number of test samples, to fit a Weibull distribution to find the RUL.


## 3.1.4. The Procedure for HI Assessment

Figure4 shows the framework used to esteem the aircraft engine HIs and, subsequently, to predict the RUL of the engine. The proposed algorithm can be classified within the condition monitoring techniques. It consists of a general framework and can be applied to any equipment. The dataset, both training and testing, are composed of the sensor readings of the considered items.Figure 4
The proposed framework approach.In the “Time Indicator Modelling” phase, a lifetime indicator (LTI) is defined. The sample number of each piece of equipment (that corresponds to the number of rows of the dataset) represents its life duration.The main idea of LTI is to model a degradation profile considering that, at the beginning of sampling, an item has the maximum reliability value (equal to 1), and when the disruption occurs, the item reliability has a minimum amount (equal to 0). The first value of LTI is equal to 1, and the last one is equal to zero, according to(16)TIm′=DURm−1DURm−2…0,(17)TImt'′=TImt'DURm,(18)LTImt=TImt′​′+1−TImt=1′​′,where DURm is the dataset length for the mth equipment. Each element of TIm indicates the remaining cycle times to the relative disruption at time t (then normalised in TIm”). Hence, each value of LTIm decreases from 1 to 0. LTIm represents the parameter to be esteemed and used in the algorithm for the HI estimation. Table 2 shows an example of LTI calculation.Table 2
An example of equations (16)–(18) application.
Steps descriptionEquation applicationSuppose machine 1 breaks after 10 sampling cycles; it follows thatDUR1=10, according to equation (16)TI1′=9876543210According to equation (17), the relative normalised valueTI1'′=0.90.80.70.60.50.40.30.20.10According to equation (18), the lifetime indicatorLTI1=10.90.80.70.60.50.40.30.20.1The “FCM Modelling” phase is the core of the proposed approach to identify the HI. Figure5 describes the iterative phase for the HI calculation reviewing the general GWO algorithm shown in Figure 3. In particular, the GWO algorithm is used for defining the weight of the relation between concept ith and jth (wij values) of the FCM matrix. The concepts of the FCM represent the working conditions of the equipment to be analysed, the sensor signals installed on the equipment, and the last concept which is the HI. In particular, since the purpose of the approach is the HI estimation using the FCM theory, the number of concepts (NC) to be used is equal to the number of reduced dataset variables’ number plus the HI (the algorithm output). This means that if the reduced dataset variables number is n, NC = n + 1.Figure 5
Framework of the proposed methodology.The iterative phase, shown in Figure5, is executed for each equipment belonging to the training dataset. In each iteration, the final α wolf position is assumed as the temporary FCMj and used as initial FCM for the next one. When the terminal equipment has been analysed, the relative FCM is considered as the optimal solution.The GWO algorithm is used to define thewij values of FCM. These values randomly in the range [−1 1] or [0 1], as required by the FCM theory, for each pack member. Then, equation (1) calculates the relative cost. By analysing Figure 5, it is possible to highlight how the lifetime indicator and the reduced dataset are used as input for the fitness function calculation.Assuming to have a reduced training dataset related toM equipment, with n main variables, it is divided in M reduced subdataset, each one related to a specific equipment. Thus, a reduced subdataset RDm (the sensors’ readings of equipment m with m = 1, 2, …, M) is available with the relative LTIm, in the form expressed in(19)RDm=Var1,1m…Varn,1m⋮⋱⋮Var1,fm…Varn,fm,LTIm=lti1…ltif.The termf identifies the instant in which the fault occurred; n represents the progressive number of the main reduced variables, and m is the number of considered equipment.The pack members’ position, obtained through the GWO application for a specific device, is given in the form of(20)FCMp,miter=w1,1…w1,NC⋮⋱⋮wNC,1…wNC,NC,where p is the pth member of the pack (p = 1, 2, …nPop), iter is the current iteration with iter ≤ MaxIt, and NC is the concepts number for the FCM algorithm defined before. Thus, FCMp,miter is the relative position of the pth pack member at the current iteration for the mth equipment.By analysing Figure5, it is possible to highlight the presence of two main “for loop.” The first one (external) is referred to the number of available equipment (M) and the second one (internal) to the maximum iteration number for GWO (MaxIt).At each iteration of the external loop (the iteration is equal to the equipment number in the dataset), the LTI related to a specific item is used as the benchmark for the positional cost calculation (if the iteration is equal to one, the LTI1 is examined). This means that the external loop has the objective of identifying the best FCMm for the mth equipment. In the internal loop, for the GWO application, all the pack members take a position within the domain space, updating it at each inner iteration. At the end of the inner loop, the best position for the considered item is identified (FCMm). The obtained FCMm is used as the initial position for the FCMm+1 identification (as long as m < M) to improve its accuracy.In the algorithm initialisation,FCM10 can be defined randomly if there is no knowledge of the involved equipment or a panel of experts cannot be established to model it, as described for the classical FCM design approach, according to the experience of each professional involved.As mentioned before, equation (1) evaluates the fitness cost value for each FCMp,miter. This step is the most critical in the whole algorithm, as highlighted by Mazzuto and Stylios [77]. Indeed, to calculate the cost connected to FCMp,miterr, it is necessary to consider all of the samples which make up RDm and LTIm. More accurately, if RDm has f samples (as described above) as well as LTIm, equation (1) has to be applied f times. Besides, since the number of iteration (k) in equation (1) depends on the function convergence or the fixed amount of repetition (FCMiter), the fitness cost evaluation requires an iteration number equal to f⋅FCMiter. Considering nPop wolves and a maximum iteration number for the GWO (MaxIt), the total number of iterations for the identification of the optimal FCM is equal to f⋅FCMiter⋅nPop⋅MaxIt.Considering thepth pack member, once its position is defined (FCMp,miter), its positional cost has to be calculated according to equation (1) for the best position identification.The main idea is to consider each sample in the reduced dataset for a specific item (RDm) as the initial array A0 at the specified iteration according to(21)Am0,j=Var1,jm…Varn,jm0,4j=1,2…f,m=1,2…M.The null value is because the output is considered within the set of FCM concepts, but it is the variable that needs to be taken into consideration.When the application of equation (1) reaches the convergent, Am∗,j will be described by(22)Am∗,j=Var1,jm^…Varn,jm^HIm∗,j,j=1,2…f,m=1,2…M,where Var1,jm^ is the convergent value for the variable Var1,jm and HIm∗,j is the esteemed output for the sample at time j for the mth engine.Once all of the samples inRDm have been processed, a final output array HIp,m∗ for the pth pack member, equation (23), will be available:(23)Hp,m∗=HIm∗,1HIm∗,2…HIm∗,f,f,m=1,2…M.Thus, the esteemed outputHIp,m∗ and the connected LTIm can be used to calculate the fitness cost value (Cp,miter), for the pth pack member and the mth item and at iteration iter, using the root mean squared error formula, as shown in(24)Cp,miter=∑j=1fLTImj−Hp,m∗,j2f.The root mean square error has been chosen because it describes efficiently how concentrated the data is around the line of best fit [78].Finally, once the optimal FCM to be used for the HI identification phases is identified, it is used to calculate the HIs in the “Training Health Indicator definitions” and “Testing Health Indicator definitions,” respectively.


## 4. Research Approach Application

To explain the proposed approach and to test its accuracy for HI modelling, the Turbofan Engine Degradation Simulation Dataset has been used. It is available online on the NASA repository website (https://ti.arc.nasa.gov/tech/dash/groups/pcoe/prognostic-data-repository/, last access July 21, 2020).The aircraft gas turbine engine has an integrated control system, which consists of a fan-speed controller and a set of controllers and limiters. In particular, it includes three high-limit regulators aimed at preventing the engine from exceeding its designed parameters [79].Several categories of signals could be used, including temperature, pressure, speed, and air ratio to monitor the condition of the aircraft gas turbine engine. The dataset is composed of 21 sensors installed in the aircraft engine’s different components, allowing the health conditions of the aircraft engine to be monitored (see Figure6). An excerpt of the used dataset is shown in Table 3. To have a complete view of the dataset, it is possible to refer to Saxena et al. [79] and Xu et al. [80].Figure 6
Sensors implemented in the aircraft engine [80].Table 3
An excerpt of the training dataset referring to Engine 1.
Id engineCycleWorking conditionsSensors123123456789…1920211142.00.8100.0445.0549.71343.41112.93.95.7137.42211.98311.3…100.010.66.41220.00.7100.0491.2606.11477.61237.59.413.6332.12323.78713.6…100.024.414.71342.00.8100.0445.0549.01343.11117.13.95.7138.22211.98306.7…100.010.56.41442.00.8100.0445.0548.71341.21118.03.95.7138.02211.98312.4…100.010.56.41525.00.660.0462.5536.11255.21033.67.19.0174.81915.27994.9…84.914.08.71635.00.8100.0449.4554.81352.91117.05.58.0193.82222.88340.0…100.014.98.9170.00.0100.0518.7641.81583.51393.914.621.6552.52387.99050.5…100.038.923.51842.00.8100.0445.0549.11344.21110.83.95.7137.12211.98307.3…100.010.66.31942.00.8100.0445.0549.61342.91101.73.95.7138.02211.98307.8…100.010.66.311025.00.660.0462.5536.41251.91041.47.19.0174.71915.28005.8…84.914.38.611120.00.7100.0491.2606.91478.01233.19.413.6333.22323.78709.6…100.024.614.711235.00.8100.0449.4554.51366.01122.75.58.0193.72222.88337.5…100.014.78.911325.00.660.0462.5536.31257.81040.97.19.0174.51915.38000.1…84.914.48.611420.00.7100.0491.2607.31470.31242.49.413.6333.72323.78714.4…100.024.314.7………………………………………………The training dataset is made up of readings from 249 engines (for a total of 61249 rows and 26 columns), while the testing dataset is made up of data from 248 engines (for a total of 41214 rows and 26 columns). The approach evaluation has been carried out using Matlab 2019© installed on a Intel® Core™ i7-6700HQ CPU @ 2.60 GHz.The results of the proposed approaches have been compared to those obtained using an artificial neural network, due to the similarity between the Artificial Neural Network (ANN) and FCM, in order to evaluate the performance of these approaches. In light of this, to have comparable results, the initial dataset has been standardised according to the working conditions and then reduced through the trendability analysis [21] to guarantee the impartiality of the data suitability obtained with the two approaches. The reduced dataset has been used as input for the proposed approach and the ANN. More specifically, according to Figure 6, the number of reduced sensors is equal to 8, such as 2, 3, 4, 8, 9, 11, 13, and 17 (see Table 4).Table 4
An excerpt of the reduced and normalised training dataset referred to Engine 1.
Id engineCycleSensors23489111317110.435−1.385−1.282−0.608−0.877−1.207−0.386−0.32712−2.472−0.888−1.357−1.660−0.878−1.416−1.907−2.27713−1.137−1.432−0.825−0.357−1.126−1.311−0.464−0.92914−1.675−1.717−0.716−0.524−0.821−1.242−0.270−1.53115−1.655−0.935−1.763−0.911−1.441−0.861−1.032−1.27216−1.637−1.837−1.298−1.059−0.894−1.579−0.976−2.29317−1.296−0.753−1.167−1.106−0.674−1.691−1.305−0.97018−0.922−1.274−1.521−0.357−1.095−1.520−0.464−0.929190.155−1.472−2.530−0.440−1.066−2.076−0.347−1.531110−1.007−1.506−0.723−0.863−0.621−1.856−0.808−1.272111−0.879−0.826−1.855−1.466−1.102−1.698−1.433−1.087112−2.1320.199−0.637−1.020−1.037−1.543−1.012−1.063113−1.085−0.487−0.789−0.626−1.055−1.110−0.763−1.272114−0.013−1.985−0.805−1.368−0.836−1.487−1.338−0.492115−1.603−0.284−1.235−0.863−0.688−1.069−0.897−1.969116−2.315−1.459−0.353−1.514−0.689−1.769−1.291−2.277117−1.741−2.353−1.497−1.383−1.152−1.489−1.470−0.721118−0.125−1.193−1.179−0.273−1.001−1.172−0.425−0.929119−1.331−1.684−1.330−1.612−1.076−1.275−1.575−1.087120−1.705−0.936−1.622−1.612−0.877−1.275−1.670−1.087121−1.627−1.728−0.751−5.281−1.553−1.015−5.341−1.550122−0.348−1.623−1.684−1.514−1.505−1.557−1.670−1.087123−1.160−1.553−1.431−1.039−1.070−1.319−1.036−1.550124−1.707−1.503−1.862−0.721−1.158−1.027−0.763−2.666125−0.771−1.574−1.050−0.566−1.165−1.172−0.502−0.929…………………………These sensors will be the concepts for the realised FCM and the input for the ANN. Table5 shows the nomenclature used for each sensors, the concepts in FCM, and the input of the ANN.Table 5
Sensors nomenclature for FCM and ANN approaches.
Sensor 2Sensor 3Sensor 4Sensor 8Sensor 9Sensor 11Sensor 13Sensor 17FCM idC1C2C3C4C5C6C7C8ANN idI1I2I3I4I5I6I7I8


### 4.1. The Proposed Approach Results

Once the training dataset has been reduced, for each engine, according to equations (16)–(18), the relative LTI array has been calculated (see Figure 7) to be used as output in the proposed approach regarding the positional cost definition.Figure 7
LTI curves for the reduced training dataset.As far as the proposed approach is concerned, as mentioned in Section3.1.3, it can be initialised either using an FCM design referring to the experience of an expert panel or with a random matrix to be iteratively corrected. Due to the lack of availability of the experts concerning the aircraft engine knowhow, for the examined case study, a random initial FCM has been adopted.Since the training dataset is composed of 249 engines, the entire process has been carried out for 249 iterations during which the FCM obtained in the previous iteration is corrected. Figure8 shows the convergence curves during the algorithm iterations and highlights the final value of the last curve that shows the minimum root mean square error equal to 1.6117. Moreover, concerning the application of equation (2), the hyperbolic tangent function has been chosen as the threshold function f() with slope factor equal to 1. The maximum number of iterations for the positional cost calculus has been fixed equal to 50 and an additional threshold value, equal to 10–3, has been defined to potentially arrest the algorithm. The required training time to identify the final FCM has been calculated to be equal to 15 minutes due to the large numbers of samples composing the dataset.Figure 8
Convergence curves at each training iterations.Table6 shows the final wij values among the concepts of the optimal FCM and the output of the proposed algorithm. The last row shows all null values being the last concept C9, the HI, an output concept.Table 6
The identified FCM.
ConceptsC1C2C3C4C5C6C7C8C9C1—0.0900.016−0.0840.0150.6040.0890.1590.266C20.106—0.4050.509−0.2600.163−0.2550.069−0.056C30.0910.045—0.1730.0290.4200.3220.7490.022C40.1930.5950.016—0.2300.0980.0470.218−0.042C50.0400.4540.1640.198—0.3220.7730.648−0.998C60.0220.3090.1330.3460.391—0.077−0.892−0.276C70.3200.2930.1660.7701.0000.555—0.790−0.200C80.1710.122−0.069−0.1560.2390.1530.329—−0.062C9—————————Analysing Table6, it is possible to highlight the presence of some low values (less than 0.1). It would be possible to filter the final FCM so that these values could be considered null to facilitate the HI calculus. However, the additional filtering phase adds to the entire process a delay factor since the user should define the filter threshold value properly through specific algorithms, increasing the iteration time. Since the size of the FCM concepts set, for the examined case study, is not so big, the additional filter phase has been neglected.The final FCM can be graphically represented to evaluate all the concatenations among concepts, as shown in Figure9.Figure 9
The final FCM in the symbolic form.Once the final FCM has been obtained, the strength of the concepts involved can be analysed so as to identify the main relevant ones and, through the critical path analysis, identify those concepts that can be considered the leading causes of the degradation profile for each engine.Specifically, according to equation (3), Table 7 shows the Total Effects’ matrix where the (i, j) values of this matrix represent how much the concept indicated in the rows affects the concept indicated in the columns.Table 7
TEs matrix among all of the involved concepts.
ConceptsC1C2C3C4C5C6C7C8C9C1—0.3090.1330.3460.3910.6040.1590.6040.391C2−0.255—0.4050.509−0.2600.4050.3220.405−0.276C30.3200.309—0.3460.3910.4200.3290.7490.391C40.1930.5950.405—−0.2600.230−0.2550.230−0.260C50.3200.4540.4050.770—0.5550.7730.773−0.998C60.1930.3910.3090.3460.391—0.391−0.8920.391C70.3200.5950.2930.7701.0000.555—0.790−0.998C80.3200.2930.1660.3290.3290.3290.329—0.329C9—————————The most significant influence is related to the relationship between concept C7 (sensor number 13) and C1 (sensor number 2) with the strength value equal to1. This relation means that C7, the Corrected fan speed, is the main cause of the increase to C1 that is the Total temperature at the LPC outlet.Focusing on the degradation profile, the HI concept C9, it is possible to highlight that concept C7, jointly with concept C5 (sensor number 9,Physical core speed) have the most significant weight since the strength value is equal to -0.998. This mean that the fan speed increasing is the relevant cause of aircraft engine degradation.Analysing in depth all the critical paths (Table8) starting from each concept of the FCM and ending in the concept C9, it is evident how the strength of the relationship between C7 and C9 is not a direct influence. This is because C7 indirectly affects C9 through the influence on C5. Thus, C5 can be considered the most relevant cause of aircraft engine degradation. This could be an important consideration for a proper maintenance plan design.Table 8
Main paths affecting HI concept (C9).
Initial nodeEnd nodeTEC1C6C5C90.391C2C3C6C9−0.276C3C6C5C90.391C4C2C5C9−0.260C5C9−0.998C6C5C90.391C7C5C9−0.998C8C7C5C90.329Once the FCM has been analysed, it can be used to calculate the HIs for each engine in the training dataset (the first 50 HIs are shown in Figure10) and also for the testing dataset. In practical terms, the HI shape provides maintenance managers with the real RUL value.Figure 10
HIs’ curves for the first 50 engines of the training dataset through the proposed approach.


### 4.2. The Comparison between the Proposed Approach with ANN

Results of the proposed approach have been compared to those obtained using an Artificial Neural Network. The ANNs were chosen for the comparison as this methodology is one of the most used in literature for the evaluation of the HI and also due to the similarity with the FCM method. As far as the ANN is concerned, the best results have been obtained considering a two-level network composed of ten neurons, respectively. Referring to the Levenberg–Marquardt method [81] and using the mean squared error as a performance indicator, Figure 11 shows how the lower MSE is obtained at 60th epoch with a value ranging between [10−2 10−1].Figure 11
Best validation performance of the used ANN used.The LTI shown in Figure7 have been used to train the ANN in order to obtain the HI estimation, as reported in Figure 12.Figure 12
HIs’ curves for the first 50 engines of the training dataset through the ANN approach.The HIs defined for the engines in the testing dataset have been used for the RUL prediction using k-neighbors algorithm. Table9 shows an excerpt of the results. Specifically, engines are reported in the ascending order in terms of the FCM percentage error (%err FCM). The second column is the real RUL value for each machine (values provided by NASA).Table 9
An excerpt of estimated RUL for the testing dataset.
Id engineTrueRULEsFCMEsANN%Err FCM%Err ANN28382.194.731.09%14.13%9176178.11170.691.20%3.02%38987.89124.811.25%40.24%11176178.84219.611.61%24.78%45351.8450.542.19%4.64%18171167.05180.092.31%5.32%56462.4676.852.41%20.08%13333.8244.62.48%35.15%77071.93100.152.76%43.07%22163167.84160.832.97%1.33%13110113.36157.723.05%43.38%16109112.73155.483.42%42.64%21113117.74131.874.20%16.70%25111105.84148.764.64%34.02%84143.0340.104.95%2.20%63334.8637.355.64%13.18%247569.87105.276.84%40.36%277670.37101.517.41%33.57%267176.31101.377.48%42.77%174649.7657.998.17%26.07%153740.4834.039.40%8.03%122622.935.5411.90%36.69%233943.9555.1312.68%41.36%………………The estimated RUL values by FCM and ANN (columns 3 and 4, respectively) underline how the proposed approach performances are better than those obtained using ANN.Figures13 and 14 compare the Weibull distributions derived from the similarity approach based on FCM and ANN for HIs. These figures show, as an example, the case of Engine 1. For all other engines, similar results have been obtained.Figure 13
Engine 1, Probability Density Function, True RUL, and Estimated RUL with proposed approach (trueRUL = 33 and estRUL = 33.82).Figure 14
Engine 1, Probability Density Function, True RUL, and Estimated RUL with ANN (trueRUL = 33 and estRUL = 44.60).The RUL obtained by the ANN approach is affected by an overestimation with respect to the FCM one. Thus, few involved curves allow the reduction of the variability range for a likely prediction. This is a typical problem of ANN that has been overcome through the proposed method.The HIs defined using the proposed approach are more accurate and, in addition, the algorithm provides significant discrimination of all the considered aircraft engines (see Figure15). Thus, small variations in the sensor readings define quite distinct degradation profiles.Figure 15
HIs’ comparison using FCM and ANN for Engine 1.
(a)(b)


## 4.1. The Proposed Approach Results

Once the training dataset has been reduced, for each engine, according to equations (16)–(18), the relative LTI array has been calculated (see Figure 7) to be used as output in the proposed approach regarding the positional cost definition.Figure 7
LTI curves for the reduced training dataset.As far as the proposed approach is concerned, as mentioned in Section3.1.3, it can be initialised either using an FCM design referring to the experience of an expert panel or with a random matrix to be iteratively corrected. Due to the lack of availability of the experts concerning the aircraft engine knowhow, for the examined case study, a random initial FCM has been adopted.Since the training dataset is composed of 249 engines, the entire process has been carried out for 249 iterations during which the FCM obtained in the previous iteration is corrected. Figure8 shows the convergence curves during the algorithm iterations and highlights the final value of the last curve that shows the minimum root mean square error equal to 1.6117. Moreover, concerning the application of equation (2), the hyperbolic tangent function has been chosen as the threshold function f() with slope factor equal to 1. The maximum number of iterations for the positional cost calculus has been fixed equal to 50 and an additional threshold value, equal to 10–3, has been defined to potentially arrest the algorithm. The required training time to identify the final FCM has been calculated to be equal to 15 minutes due to the large numbers of samples composing the dataset.Figure 8
Convergence curves at each training iterations.Table6 shows the final wij values among the concepts of the optimal FCM and the output of the proposed algorithm. The last row shows all null values being the last concept C9, the HI, an output concept.Table 6
The identified FCM.
ConceptsC1C2C3C4C5C6C7C8C9C1—0.0900.016−0.0840.0150.6040.0890.1590.266C20.106—0.4050.509−0.2600.163−0.2550.069−0.056C30.0910.045—0.1730.0290.4200.3220.7490.022C40.1930.5950.016—0.2300.0980.0470.218−0.042C50.0400.4540.1640.198—0.3220.7730.648−0.998C60.0220.3090.1330.3460.391—0.077−0.892−0.276C70.3200.2930.1660.7701.0000.555—0.790−0.200C80.1710.122−0.069−0.1560.2390.1530.329—−0.062C9—————————Analysing Table6, it is possible to highlight the presence of some low values (less than 0.1). It would be possible to filter the final FCM so that these values could be considered null to facilitate the HI calculus. However, the additional filtering phase adds to the entire process a delay factor since the user should define the filter threshold value properly through specific algorithms, increasing the iteration time. Since the size of the FCM concepts set, for the examined case study, is not so big, the additional filter phase has been neglected.The final FCM can be graphically represented to evaluate all the concatenations among concepts, as shown in Figure9.Figure 9
The final FCM in the symbolic form.Once the final FCM has been obtained, the strength of the concepts involved can be analysed so as to identify the main relevant ones and, through the critical path analysis, identify those concepts that can be considered the leading causes of the degradation profile for each engine.Specifically, according to equation (3), Table 7 shows the Total Effects’ matrix where the (i, j) values of this matrix represent how much the concept indicated in the rows affects the concept indicated in the columns.Table 7
TEs matrix among all of the involved concepts.
ConceptsC1C2C3C4C5C6C7C8C9C1—0.3090.1330.3460.3910.6040.1590.6040.391C2−0.255—0.4050.509−0.2600.4050.3220.405−0.276C30.3200.309—0.3460.3910.4200.3290.7490.391C40.1930.5950.405—−0.2600.230−0.2550.230−0.260C50.3200.4540.4050.770—0.5550.7730.773−0.998C60.1930.3910.3090.3460.391—0.391−0.8920.391C70.3200.5950.2930.7701.0000.555—0.790−0.998C80.3200.2930.1660.3290.3290.3290.329—0.329C9—————————The most significant influence is related to the relationship between concept C7 (sensor number 13) and C1 (sensor number 2) with the strength value equal to1. This relation means that C7, the Corrected fan speed, is the main cause of the increase to C1 that is the Total temperature at the LPC outlet.Focusing on the degradation profile, the HI concept C9, it is possible to highlight that concept C7, jointly with concept C5 (sensor number 9,Physical core speed) have the most significant weight since the strength value is equal to -0.998. This mean that the fan speed increasing is the relevant cause of aircraft engine degradation.Analysing in depth all the critical paths (Table8) starting from each concept of the FCM and ending in the concept C9, it is evident how the strength of the relationship between C7 and C9 is not a direct influence. This is because C7 indirectly affects C9 through the influence on C5. Thus, C5 can be considered the most relevant cause of aircraft engine degradation. This could be an important consideration for a proper maintenance plan design.Table 8
Main paths affecting HI concept (C9).
Initial nodeEnd nodeTEC1C6C5C90.391C2C3C6C9−0.276C3C6C5C90.391C4C2C5C9−0.260C5C9−0.998C6C5C90.391C7C5C9−0.998C8C7C5C90.329Once the FCM has been analysed, it can be used to calculate the HIs for each engine in the training dataset (the first 50 HIs are shown in Figure10) and also for the testing dataset. In practical terms, the HI shape provides maintenance managers with the real RUL value.Figure 10
HIs’ curves for the first 50 engines of the training dataset through the proposed approach.


## 4.2. The Comparison between the Proposed Approach with ANN

Results of the proposed approach have been compared to those obtained using an Artificial Neural Network. The ANNs were chosen for the comparison as this methodology is one of the most used in literature for the evaluation of the HI and also due to the similarity with the FCM method. As far as the ANN is concerned, the best results have been obtained considering a two-level network composed of ten neurons, respectively. Referring to the Levenberg–Marquardt method [81] and using the mean squared error as a performance indicator, Figure 11 shows how the lower MSE is obtained at 60th epoch with a value ranging between [10−2 10−1].Figure 11
Best validation performance of the used ANN used.The LTI shown in Figure7 have been used to train the ANN in order to obtain the HI estimation, as reported in Figure 12.Figure 12
HIs’ curves for the first 50 engines of the training dataset through the ANN approach.The HIs defined for the engines in the testing dataset have been used for the RUL prediction using k-neighbors algorithm. Table9 shows an excerpt of the results. Specifically, engines are reported in the ascending order in terms of the FCM percentage error (%err FCM). The second column is the real RUL value for each machine (values provided by NASA).Table 9
An excerpt of estimated RUL for the testing dataset.
Id engineTrueRULEsFCMEsANN%Err FCM%Err ANN28382.194.731.09%14.13%9176178.11170.691.20%3.02%38987.89124.811.25%40.24%11176178.84219.611.61%24.78%45351.8450.542.19%4.64%18171167.05180.092.31%5.32%56462.4676.852.41%20.08%13333.8244.62.48%35.15%77071.93100.152.76%43.07%22163167.84160.832.97%1.33%13110113.36157.723.05%43.38%16109112.73155.483.42%42.64%21113117.74131.874.20%16.70%25111105.84148.764.64%34.02%84143.0340.104.95%2.20%63334.8637.355.64%13.18%247569.87105.276.84%40.36%277670.37101.517.41%33.57%267176.31101.377.48%42.77%174649.7657.998.17%26.07%153740.4834.039.40%8.03%122622.935.5411.90%36.69%233943.9555.1312.68%41.36%………………The estimated RUL values by FCM and ANN (columns 3 and 4, respectively) underline how the proposed approach performances are better than those obtained using ANN.Figures13 and 14 compare the Weibull distributions derived from the similarity approach based on FCM and ANN for HIs. These figures show, as an example, the case of Engine 1. For all other engines, similar results have been obtained.Figure 13
Engine 1, Probability Density Function, True RUL, and Estimated RUL with proposed approach (trueRUL = 33 and estRUL = 33.82).Figure 14
Engine 1, Probability Density Function, True RUL, and Estimated RUL with ANN (trueRUL = 33 and estRUL = 44.60).The RUL obtained by the ANN approach is affected by an overestimation with respect to the FCM one. Thus, few involved curves allow the reduction of the variability range for a likely prediction. This is a typical problem of ANN that has been overcome through the proposed method.The HIs defined using the proposed approach are more accurate and, in addition, the algorithm provides significant discrimination of all the considered aircraft engines (see Figure15). Thus, small variations in the sensor readings define quite distinct degradation profiles.Figure 15
HIs’ comparison using FCM and ANN for Engine 1.
(a)(b)


## 5. Discussion

The proposed approach has the ability to operate in a dynamic environment with no significant difference in the operation of the algorithm in steady state or dynamic mode, guaranteeing a reliable and robust performance together with an easy implementation.At the same time, it requires particular attention by users in defining all the involved parameters such as the size of the dataset, the number of agents to be used to find the final FCM, and the threshold values. Indeed, as discussed, the total number of iterations and therefore the total computational time to calculate HIs depends on them. However, the analysed case study has highlighted how this limitation can be overcome by applying, before the algorithm initialisation, a dataset reduction to minimise the involved variables number.Probably the most significant advantage of the swarm’s intelligence is its ability to operate in a dynamic environment. The swarm can continuously follow the path even for rapidly evolving optimisation. In principle, there is no significant difference in the operation of the algorithm in steady state or dynamic mode [82]. Moreover, these algorithms do not require knowledge, for example, about the gradients of the cost function and constrained functions. They guarantee reliable and robust performance together with an easy implementation [83]. Specifically, as far as the GWO is concerned, the most important one refers to the small number of parameters needed for its implementation and adjustment [84, 85].At the same time, regarding the utilisation of FCMs, among several advantages, the most important is their extreme flexibility and adaptability to a given domain, allowing qualitative simulation of a system once constructed. Furthermore, FCMs symbolically represent knowledge, converting the relations between the elements of a mental landscape to assess the impact of these elements [86, 87]. The use of FCMs demonstrates other additional benefits, including the use of fuzzy logic. Indeed, the fuzzy set theory allows the incorporation of uncertainty due to sparse and imprecise information [88]. A fuzzy value is a fuzzy representation of a specific property when it is not precisely known [89]. The fuzzy set theory and numbers are mainly used to quantify the grade to which a property can be connected with an object. It must not be confused with the concept of probability. Indeed, the causality among concepts is considered as a certainty, since the concept of causality is not used to try to identify or find relationships between factors such as structural equation model and/or Bayesian nets [90].


## 6. Conclusion

In this paper, an innovative supervised approach that combines a Swarm Intelligence algorithm, the GWO, and FCMs is proposed for HI analysis and calculation. This approach allows maintenance managers to predict the RUL of items through the use of k-neighbors algorithms as well as to have an in-depth understanding of the degradation process; thanks to the analysis of the main paths of concepts that affect the HI. In order to enhance the operating reliability and reduce maintenance costs, an integrated fault diagnosis and prognosis framework that analyse the machinery degradation process is necessary.In the proposed approach, the working conditions of the engines and the sensor signals installed on engines become the concepts of the FCM, while the GWO, a Swarm Intelligence algorithm, has been used for defining the connection weight among these concepts and the HI concepts.A dataset provided by NASA that concerns the data of aircraft engines has been used to test the proposed approach. The case used underlines a crucial aspect. Comparing the results with those obtained through neural networks, the proposed algorithm models, and all of the degradation profiles in a more detailed manner allows one to significantly distinguish different situations without imposing any specified mathematical functions. This consideration is reflected in fewer profiles that can be considered similar to the case in question and, consequently, give a more precise estimate of the RUL. Moreover, analysing the final FCM, the physical core speed and the corrected fan speed have been identified as the main critical factors to the engine degradation.Furthermore, the use of the FCM approach allows the user to be able to analyse in an intuitive way the relationships between the variables involved and thus have a greater understanding of the degradation process, which is impossible for an ANN. Indeed, in an ANN, the variables involved are the inputs for the system and the neurons concatenation has no meaning to understand the process. On the contrary, in an FCM, the variables are simultaneous inputs and “neurons,” so their concatenation gives more information about the process.The performance of the proposed approach has been demonstrated using a NASA dataset, but it can also be applicable to the other fault diagnosis and prognosis equipment. A wide range of experiments will be performed to investigate the robustness of the proposed method in our next step research. At the same time, it is evident how the proposed approach can be based not only on a feature reduction but also on the determination of the most useful items for the training phase. Indeed, considering all the variables involved for the algorithm application (number of wolves, maximum iteration number for the GWO, thresholds etc.), the total number of iterations for the identification of the optimal FCM can be very huge and time-consuming. For this reason, as further development is crucial to design a preliminary step to be used after feature extraction step.



---

*Source: 8832011-2021-02-17.xml*