Optimisation of Maintenance Policies Based on Right-Censored Failure Data Using a Semi-Markovian Approach

  Citation: Sánchez-Herguedas, A.; Mena-Nieto, A.; Rodrigo-Muñoz, F.; Villalba-Diez, J.; Ordieres-Meré, J. Optimisation of Maintenance Policies Based on Right-Censored Failure Data Using a Semi-Markovian Approach. Sensors 2022,22, 1432. https://doi.org/10.3390/s22041432 Academic Editor: Anastasios Doulamis Received: 29 December 2021 Accepted: 10 February 2022 Published: 13 February 2022 Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. Copyright: © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). sensors Article Optimisation of Maintenance Policies Based on Right-Censored Failure Data Using a Semi-Markovian Approach Antonio Sánchez-Herguedas 1, Angel Mena-Nieto 2, Francisco Rodrigo-Muñoz 3, Javier Villalba-Díez 4 and Joaquín Ordieres-Meré 5,* 1Department of Industrial Management, School of Engineering, University of Seville, Camino de los Descubrimientos s/n, 41092 Seville, Spain; [email protected] 2Department of Electrical and Thermal Engineering, Design, and Projects, School of Engineering, Campus El Carmen, University of Huelva, 21819 Huelva, Spain; [email protected] 3Department of Applied Mathematics II, School of Engineering, University of Seville, Camino de los Descubrimientos s/n, 41092 Seville, Spain; [email protected] 4Hochschule Heilbronn, Fakultät Management und Vertrieb, Campus Schwäbisch Hall, 74523 Schwäbisch Hall, Germany; javier[email protected] 5Escuela Técnica Superior de Ingenieros Industriales (ETSII), Universidad Politécnica de Madrid, 28006 Madrid, Spain *Correspondence: j.ordier[email protected]; Tel.: +34-910677107 Abstract: This paper exposes the existing problems for optimal industrial preventive maintenance intervals when decisions are made with right-censored data obtained from a network of sensors or other sources. A methodology based on the use of the z transform and a semi-Markovian approach is presented to solve these problems and obtain a much more consistent mathematical solution. This methodology is applied to a real case study of the maintenance of large marine engines of vessels dedicated to coastal surveillance in Spain to illustrate its usefulness. It is shown that the use of right-censored failure data significantly decreases the value of the optimal preventive interval calculated by the model. In addition, that optimal preventive interval increases as we consider older failure data. In sum, applying the proposed methodology, the maintenance manager can modify the preventive maintenance interval, obtaining a noticeable economic improvement. The results obtained are relevant, regardless of the number of data considered, provided that data are available with a duration of at least 75% of the value of the preventive interval. Keywords: maintenance interval; maintenance model; semi-Markov process; right-censored data; finite horizon; maintenance cost 1. Introduction The main goal of this work is to present a methodology that allows finding the optimal maintenance preventive interval when the failure data are right-censored. In particular, we are interested in finding out how the use of right-censored data affects the calculation of the optimal maintenance interval, considering the temporary maintenance schedule of maintenance interventions and the different types of costs incurred. This study is relevant, since these are the types of data that are generally available in the vast majority of companies, which the maintenance manager must use. When the maintenance engineer is faced with the problem of managing equipment subject to predetermined preventive maintenance, one of the tasks is to determine the preventive maintenance interval that economically optimises (other factors could also be optimised) the series of interventions that can be carried out on a physical asset over a specified period [ 1 ]. For this predetermined maintenance, the interventions can be of two types: corrective interventions originated by a failure of some component of the equipment (failure mode) and preventive interventions carried out after a certain number of hours, kilometres, or units produced. The former is produced randomly, and their number and Sensors 2022,22, 1432. https://doi.org/10.3390/s22041432 https://www.mdpi.com/journal/sensors Sensors 2022,22, 1432 2 of 17 location in time depend on the failure behaviour of the component with that particular failure mode. The second can be located in time, and the number of can be established along any considered time horizon. In those cases where the failure mode shows wear behaviour, these types of interventions are not independent of each other, since increasing the number of these interventions reduces the number of those due to component failures. However, the increase in the number of the second type of intervention implies an increase in the associated costs of preventive maintenance, but a reduction in the costs of corrective maintenance. Therefore, the optimal preventive maintenance interval that minimises maintenance costs (preventive and corrective) must be found. To find the optimal interval, the starting point is the existing status and the available information, which must reflect the behaviour in the event of failure, the costs of maintenance interventions, and the income generated by the operation of the equipment [ 2 ]. The information on the cost of each intervention is generated during the operation and maintenance phase of the equipment and is usually collected in the Computer-Aided Maintenance Management Systems (CMMS). Once the learning period for maintenance tasks is over, this information is usually relatively stable over time. In the same way, the information on the income from the operation of the equipment is usually collected in the Enterprise Resource Planning systems (ERPs) [ 3 ]. Information on behaviour in the event of failure can be obtained from various sources, the equipment manufacturer, the history of failures collected in the plant itself, or any other external database that presents operating conditions similar to that of the equipment under study. The information referring to equipment failures is usually right-censored due to the performance of preventive maintenance. On rare occasions, the maintenance manager has failure data where there is no right-censored data. This information, if it exists, is usually in the hands of manufacturers who have tested their equipment up to its failure point. The maintenance manager cannot bear the test cost, so the adopted policy usually follows the preventive interval set by the manufacturer. Such settings are frequently conditioned by economic interests or brand image. In the research carried out, an extensive data set was available; therefore, the knowledge based impact of right censoring can be identified. The literature on maintenance is vast, including Reliability Engineering [ 4 – 12 ], maintenance policies and modelling [ 13 – 20 ], and optimal preventive maintenance intervals [ 21 – 24 ]. However, fewer papers address the problem of uncertainty in lifetime distribution [ 25 – 28 ], either by analysing several time-based maintenance policies having uncertainty in the parameters of the lifetime distribution [ 29 ] or using right-censored data [ 30 ], and adopt the Markovian approach of transition between states (operational, preventive, and corrective maintenance intervention) [31–36]. The discrete-time and continuous-time Markov models satisfy the Markov property, i.e., the future only depends on the present, not on the past. Only the transitions between states and their respective probabilities are considered in the discrete-time model (Markov chain). The sojourn time in each state is irrelevant. In the continuous-time model, the sojourn time is a random variable with an exponential distribution due to the Markovian property. By contrast, in semi-Markov models, the sojourn time in each state does not follow an exponential distribution. It implies that semi-Markov models do not satisfy the Markovian property, which relaxes the constraints and improves the application value [ 37 ]. Therefore, successive transitions between states form a Markov chain, called an embedded Markov chain in the semi-Markov model. A wider review, which provides an interesting perspective, can be found in [7,38,39]. Additionally, other authors address the problem of using right-censored data by proposing alternative methods, for example, Li et al. [40] found the life distributions using a histogram-based technique to graphically obtain the maintenance interval. Mazzuchi and van Dorp [26] used the Newton–Raphson method and MLE to try to represent the remaining life of coupler knuckles in railway wagons. Taghipour and Banjevic [41] used the likelihood ratio test to check for trends in the failure data and the ME algorithm to find the parameters for trend analysis applied to pumps located in hospitals. Sensors 2022,22, 1432 3 of 17 However, to the best of our knowledge, very few studies provide rules and guidelines to modify the default maintenance interval optimally, and this paper helps to fill that gap. When the maintenance manager has censored data due to the implementation of preventive tasks, the proposed methodology represents an advance with respect to previous works such as [21–24,42]. This article differs from previous studies in providing the maintenance manager with a prediction of the length of the optimal preventive interval as well as a decision rule to improve the interval being used, based on the analysis of the results of various scenarios, and by using the methodology able to deal with right-censored data. This rule guides the manager in the direction of modification of the preventive interval to be used (increasing or decreasing it) and the intensity of the modification. In this way, the uncertainty induced by using censored data can be overcome. Therefore, to fill the identified gaps, two research questions are formulated: RQ1: How far is the obtained result from the optimal solution when the failure data are censored? RQ2: What strategy would help to move closer to the optimum based on censored data? The developed methodology guarantees that the solution found is optimal from a mathematical point of view, incorporating the time right-censored data, together with different types of costs and income. In addition, it applies a Semi–Markovian approach of transition between states [ 43 – 45 ], plus the z -transform, to find the optimal maintenance interval, which, as far as we know, is the first time it has been applied with right-censored data. It is also relevant that the developed methodology enables the assessment and comparison of alternative maintenance policies. The rest of the paper is organised as follows. Section 2provides the real data for this study, the research questions formulated, and the methodology applied to develop the mathematical optimisation model based on semi-Markov processes and z-transform. Section 3presents the results for the different cases studied. In Section 4, those results are highlighted, discussed, and compared. Finally, Section 5provides the conclusions drawn from the study. 2. Materials and Methods To adapt the preventive interval to the particular equipment conditions usage, the manager responsible for maintenance should use the information generated by the equipment itself. The failure information generated is right-censored and it is a challenge to determine the behaviour before the failure of the equipment. The method followed to find the optimal interval when the available failure data is right-censored is explained throughout the article. This process is graphically summarised in Figure 1. 2.1. Real Case—Data Selection and Information Processing This work analyses the behaviour of the O-rings located in the cooled crossover of a 12V diesel engine with 2 litres per cylinder. The crossover is a jacketed tube through which exhaust gases flow. The engine coolant cools this tube, and two O-rings are installed at each end of the tube to prevent coolant leakage. High temperatures affect the O-rings and are responsible for their degradation. In this experiment, failure data have been collected for several years, following the preventive replacement policy of O-rings every 4000 h (scenario A). During this time, it was observed that none of the O-rings had reached this limit, remaining at values very far from this value. This fact allows us to consider that the values of the observed failures can be considered as “not” censored. These failures are listed in Table 1. In addition, in some cases, the O-rings were replaced preventively because the last failure occurred close to the time of the preventive replacement. These are the censored values that are collected in Table 2. Sensors 2022,22, 1432 4 of 17 Design of the Semi-Markov Model for maintenance Obtaining the system of equations for average accumulated return Obtaining the formula of the average accumulated return by applying z transforms Obtaining the formula of optimal preventive interval, which maximizes the average accumulated return Section 2.3 Section 2.4 Section 2.5 Section 2.6 Obtaining the optimal preventive interval value t0 for two group of cases. Time to failure and preventive right-censored data Calculation of the observed function Weibull function adjustment Economic data and maintenance times Section 2.1 Section 2.2 Section 3 Figure 1. Graphical overview of the process. Table 1. Hours of operation of the O-rings until failure in Scenario A. Operation Hours to Failure 190 276 296 409 429 430 437 454 481 492 498 499 543 552 552 552 552 577 603 604 604 612 619 658 675 683 696 702 742 754 773 797 812 836 881 889 912 913 942 974 994 994 1014 1015 1024 1025 1041 1105 1183 1203 1211 1236 1238 1240 1249 1274 1295 1304 1312 1343 1345 1407 1413 1421 1442 1447 1486 1492 1542 1581 1601 1621 1630 1675 1735 1892 1926 1960 2006 2101 2242 2437 2450 Table 2. Hours of operation of the O-rings with censorship due to preventive maintenance policy in scenario A. Operation Hours to the Censorship 45 84 84 103 121 176 199 217 259 324 371 387 428 508 514 519 638 655 661 735 760 764 766 769 790 867 986 1006 1011 1111 1167 1188 1188 1395 1396 1505 1752 2016 Subsequently, the value of the preventive interval was modified, settling at 1000 h. In this case, the failure data were considerably reduced as many observations were rightcensored by the predetermined preventive maintenance performed at 1000 h (scenario B). This censorship radically changed the way of analysing the problem, since the results will depend less on the type of technique or method used and more on the type of censorship [ 46 ] and the extent of the censorship. Scenario B, corresponding to the preventive interval of 1000 h, is the one that is usually presented to the maintenance manager. Sensors 2022,22, 1432 5 of 17 Failure data and data censored by preventive corresponding to scenario B are, respectively, represented in Tables 3and 4. In Table 3, only seven of the nine failures that occurred during the period studied appear, since two of them, those corresponding to 171 h, did not correspond to the failure mode studied. Table 3. Hours of operation of the O-rings until failure in Scenario B. Operation Hours to Failure 242 480 522 633 663 839 845 Table 4contains 25 censored values and 87 values corresponding to a 1000-h preventive cycle completed. Table 4. Hours of operation of the O-rings with censorship due to preventive maintenance policy in scenario B. Operation Hours to the Censorship 93 93 128 128 128 128 155 161 220 220 240 240 337 367 380 380 467 467 478 485 485 520 758 829 829 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 2.2. Determination of the Failure Distribution Function Data from maintenance interventions are collected at the machine or through an automated procedure using a network of sensors. Its mathematical use for optimisation requires treating these data until obtaining information on the appropriate protocol that can be understood by the model proposed for the simulation. This treatment often supposes a loss of veracity. In other cases, it may happen that the information on which the treatment is started is not adequate. Achieving the most reliable information to the original exposed in the desired protocol is one of the objectives to be achieved in any optimisation process. In our case, the original starting data could be classified into two groups, those generated during the maintenance activities and those generated during the economic management of those maintenance activities. The first group included failure data and preventive activity data. These data are the hours of operation of the equipment in which preventive and corrective interventions are carried out and the duration of these interventions. In the first group, they were counted from the last repair or preventive change until the failure appeared or until a preventive task was carried out, and the equipment was inactive due to the intervention. The second group of data included the costs of each of the corrective and preventive interventions that were carried out and the income obtained from the use of the equipment. For the first data group, it was necessary to establish the method to determine the method of mathematically obtaining the failure distribution. Many examples can be found in the literature [ 47 , 48 ]. In our case, a Weibull probability distribution function was used, although other authors analyse different methods to estimate the parameters of the Weibull distribution [ 49 , 50 ]. Other authors have previously used other distributions. However, the Weibull distribution is the one that best adapts to the process of failure appearance in industrial assets. Furthermore, the Weibull distribution includes the exponential distribution and has the advantage of using two or three parameters instead of a single parameter of the Sensors 2022,22, 1432 6 of 17 exponential function. The exponential survival function has the hazard rate constant, while the Weibull survival distribution extends the exponential distribution to allow constant, increasing, or decreasing hazard rates. The estimation of the model parameters becomes a difficult task when censorship radically changes the available information and when the opinion of experts does not help to understand the system’s behaviour [51]. Many authors have estimated the parameters of the Weibull distribution for censored data through the Maximum Likelihood Estimation method (MLE), using algorithms for their numerical resolution such as the Maximised Expectation (EM) algorithm, see [ 52 – 54 ]. Other authors use and compare various methods such as the Synthetic Minority Oversampling TEchnique (SMOTE) [ 55 ], MLE, Least Square Estimation (LSE) [ 56 ], Weibull [ 57 ] probability plot, regression of the range of medians with the Benard approximation, and even combine methods, see [ 58 , 59 ]. Bayesian estimators [ 60 ] and other types of linear estimators [61] have also been used. At the same time, the repair and preventive replacement average times and the cost and income average of each type of intervention (corrective or preventive) are also calculated. This last type of data generally poses more difficulties to obtain than to calculate it. However, despite not being as obvious as the calculation of the previous means, the calculation of the failure distribution function can lead to significant management errors, mainly due to the lack of coherence of the starting data. When these failure data are right–censored, due to the performance of preventive maintenance, poor optimisation of the preventive interval will be achieved. This article focuses on the influence of the maintenance management data (first group) to determine the distribution function to perform the preventive interval optimisation calculations. In this research, a reduced number of data will be used, with which we will obtain an exact solution for the optimal preventive interval. If we had used a large number of data, as is the case with information collected online through sensors, the result would have been the same (exact solution). Nevertheless, in both cases, as will be shown, the results will be far from the optimal solution. The equipment was failing and was repaired, with few cases of preventive interventions (when the scheduled time according to predetermined maintenance was reached), which also lasted for a few hours of operation. This case represents a situation where the censorship was practically nonexistent, and the observed distribution function was very close to the real one. On the other hand, Tables 3and 4represent the failure and preventive data when a predetermined preventive threshold was established at 1000 h of operation. In this second case, the behaviour of the O-rings would not be known beyond 1000 h. In both cases, it does make sense to estimate the behaviour beyond the predetermined threshold or to analyse the differences between the optimum intervals in both cases. Therefore, a procedure to obtain estimations from the second case closer to those obtained for the first one can be distilled. Looking to gather more information, the following procedure was adopted: • From the failure data and preventive maintenance data collected in Tables 3and 4, the observed function was calculated applying the range of medians method and using the Benard approximation. • From the observed function, the theoretical Weibull function that fit the best was determined. Then, the least-squares method was used to obtain the theoretical function’s two or three parameters. • Subsequently, the theoretical function was used to optimise the preventive interval economically. In Figure 2, the procedures for the data in Tables 3and 4are summarised. The Weibull failure distribution function that best fit the observed function had two parameters, the shape parameter with value 1.88 and a scale parameter (characteristic life) with value 3603. These data were very different from the data obtained for the case of uncensored data and would give rise to a very different preventive interval. The detailed process for Sensors 2022,22, 1432 7 of 17 the data in Tables 1and 2can be verified in Sánchez-Herguedas et al. [42] . The Weibull functions adjusted to the uncensored data were, in that case, W(2.36, 1317, 0) and W(1.95, 1202, 117), where the first value corresponded to the shape parameter, the second to the scale parameter, and the third to the guaranteed life parameter. Figure 2. Weibull distribution function obtained from the observed function, which has been determined by Median Rank Regression (MRR) procedure and the Benard approximation. 2.3. Semi-Markov Maintenance Model with Returns for a Finite Period To study the behaviour of O-rings concerning failure, we need a mathematical model that reflects this behaviour, such as the semi–Markovian model of three states, which is used with the data that appear in Tables 1and 2. This model applies to the predetermined preventive maintenance of much industrial equipment, since it contains the two most representative maintenance states, the state of the equipment when corrective maintenance is performed after a failure (State S2) and the state of the equipment when preventive maintenance is performed after a predetermined number of operating hours (State S3). In these two cases analysed, the data correspond to preventive intervals of 4000 and 1000 h of operation, respectively. State S1 corresponds to the period in which the equipment is performing its required function. In this model, the equipment and others with the same characteristics evolve over time, changing the state according to the law of probability. When the system is in the operational state, the probability of equipment failure is distributed as a Weibull distribution. If the equipment fails, the system goes to the corrective state (S2) and develops a corrective intervention. Suppose the equipment does not fail after a time τ , the system transitions to the preventive state, developing a preventive intervention. Once these maintenance interventions (corrective and preventive) have been carried out, the system returns to its operational state. Both types of interventions have associated costs. On the other hand, during the operational state, the equipment generates income. Over time, these costs and income are accumulated in a variable called the average accumulated return, V(m) . The optimisation of this variable will allow the calculation of the optimal preventive interval, τ0 . These transitions between states and the accumulation of returns are graphically expressed in Figure 3. Figure 3. Transition process between states and accumulation of returns. Sensors 2022,22, 1432 8 of 17 2.4. Formulating a Difference Equation System for Average Accumulated Return The model that reflects the behaviour of the system is a semi–Markovian model. The variable to optimise is the average accumulated return in the successive transitions between states. In m transitions, starting from state i, the system accumulates m returns, which added to their respective signs (+ for income and − for costs) will constitute the accumulated return Qi(m) in m transitions from state i . For each value of i and m , the accumulated return is a random variable because once the initial state is determined, the next state is unpredictable, as well as the following states. Thus, in m transitions, the system can evolve in many ways. For this reason, it is interesting to know the average accumulated return instead of the accumulated return itself. Let vi(m) = E(Qi(m)) be the average accumulated return in m transitions when the system starts from the initial state i . To calculate vi(m) , a difference equation is constructed, separating the m transitions into two stages. The first stage is the one constituted by the transition from the initial state i to the next state j . As the second state can be any of the states of the system, the return vi( 1 ) in a single transition is a random variable that can reach the values ri1( 1 ) , ri2( 1 ) , · · · , rin( 1 ) , with the respective probabilities pi1 , pi2 , · · · , pin . The average return in that transition can be formulated as follows: vi(1) = n ∑ j=1 rij(1)·pij. (1) Now, the process continues with the remaining m− 1 transitions. Once the first transition has been made, the system is placed in a state called j , where j takes one of the values 1,2,3, · · · , n . The average accumulated return vj(m− 1 ) in the following transitions is a random variable that can reach the values v1(m− 1 ) , v2(m− 1 ) , · · · , vn(m− 1 ) , with probabilities of pj1 , pj2 , · · · , pjn , respectively, which remain constant throughout the m transitions, since the process is homogeneous. Therefore, the expected value of the return of the remaining m− 1 transitions can be formulated as: vj(m− 1 )) = ∑n q=1vj(m− 1 )·pjq . We conclude that the expected average return of the system in m transitions is calculated according to Equation (2): vi(m) = vi(1) + vj(m−1)·pj=vi(1) + n ∑ q=1 vj(m−1)·pjq (2) The vector V(m) is defined as V(m)=(v1(m) , v2(m) , · · · , vn(m))t , where the superscript tindicates the transpose in the matrix sense. This last equality can be written as a difference equation: V(m) = V(1) + P·V(m−1)(3) where P is the matrix of transition probabilities. The Equation (3) gives the average accumulated return in mtransitions from any possible starting state i. 2.5. Solving the System of Difference Equations by Applying Transforms Z Solving this difference equation requires the use of z transform and Laurent’s series. To reach its resolution, it is necessary to include the data and functions involved in the calculation: distribution functions of the time spent in each state, the returns for each state, and the matrices that describe the process. Time-related probabilistic functions for the tree states are as follows: •F(t)is the Cumulative Distribution Function (CDF) of equipment failures, and f(t) is the Probability Density Function (PDF). Based on the discussion in Section 2.3, we will use the three-parametric Weibull distribution. •G(tc) is the distribution function of the time the equipment remains under corrective maintenance and g(tc)is its probability density function. •H(tp) is the distribution function of the time that the equipment remains under preventive maintenance and h(tp)is its probability density function. Sensors 2022,22, 1432 9 of 17 Returns (costs for corrective and preventive states and income from operational state) are as follows: •R1, income per time unit that the system remains in State 1 (S1: Operational), 6 e/h. •R12, the cost of transition from State 1 to State 2, −4320 e. •R13, the cost of transition from State 1 to State 3, −1e. •R2 , the cost per time unit that the system remains in State 2 ( S2 : Corrective), − 95 e /h. •R21, cost of transition from State 2 to State 1, −620 EUR. •R3 , the cost per time unit that the system remains in State 3 (S3: Preventive), −82 e/h. •R31, the cost of transition from State 3 to State 1, −620 e. Matrices to describe the process are as follows: •P , the transition probability matrix between states ( pij is the probability of going from State ito j). •F , stay time matrix ( Fij is the average time in State i before the system goes to State j ). •R, returns matrix, where rij =Ri+Rij. With these data, making the z -transform of V(m+ 1 ) = V( 1 ) + P·V(m) , during the calculation, it is required to determine the matrix I−z−1Pto reach Equation (4). Z[V(m)] = 1 z−1V(1) + (I−z−1P)−1·V(1)(4) Developing the matrix (I−z−1P)−1 , the value of Z[V(m)] is reached and decomposed into simple fractions for each of its terms Z[v1(m)] , Z[v2(m)] , and Z[v3(m)] . Subsequently, using the appropriate Laurent expansion, the inverse z -transform of each one is calculated, obtaining the equations that describe the average accumulated returns in each transition for each of the starting states of the process, v1(m) , v2(m) , and v3(m) . The full development can be followed in Sánchez Herguedas et al. [62] . The average accumulated return in m transitions starting in the operating state is described by Equation (5). v1(m) =1 4h(2m+1+ (−1)m−1)·R1·Zτ 0t·f(t)dt +R12 ·F(τ) + (R1τ+R13)·(1−F(τ))+ (2m−1−(−1)m−1)·R2·Z∞ 0tc·g(tc)dtc+R21·F(τ)+ R3·Z∞ 0tp·h(tp)dtp+R31·(1−F(τ))i (5) The average accumulated return in m transitions starting in the corrective state is described by Equation (6). To deduce it, the same reasoning is followed as in the previous case . v2(m) =1 4h(2m−1−(−1)m−1)·R1·Zτ 0t·f(t)dt +R12 ·F(τ) + (R1·τ+R13)·(1−F(τ))+ (2m+1+ (−1)m−1)·R2·Z∞ 0tc·g(tc)dtc+R21+ (2m−3+ (−1)m−1)· R3·Z∞ 0tp·h(tp)dtp+R31−R2·Z∞ 0tc·g(tc)dtc+R21·(1−F(τ))i (6) The average accumulated return in m transitions starting in the preventive state is described by Equation (7). v3(m) =1 4h(2m−1−(−1)m−1)·R1·Zτ 0t·f(t)dt +R12 ·F(τ) + (R1·τ+R13)·(1−F(τ))+ (2m+1+ (−1)m−1)·R3·Z∞ 0tp·h(tp)dtp+R31+ (2m−3+ (−1)m−1)·R2·Z∞ 0tc·g(tc)dtc+R21−R3·Z∞ 0tp·h(tp)dtp+R31·F(τ)i (7) 2.6. Optimisation of the Averaged Accumulated Return When Starting from the Operational State Deriving and equalising to zero is possible to achieve the mathematical formula of the preventive interval that economically optimises the average accumulated return for Sensors 2022,22, 1432 16 of 17 18. Martín, M.G.; Álvarez, A.P.; Ordieres-Meré, J.; Villalba-Díez, J.; Morales-Alonso, G. New Business Models from Prescriptive Maintenance Strategies Aligned with Sustainable Development Goals. Sustainability 2020,13, 216. [CrossRef] 19. Ansari, F.; Glawar, R.; Nemeth, T. PriMa: A prescriptive maintenance model for cyber-physical production systems. Int. J. Comput. Integr. Manuf. 2019,32, 482–503. [CrossRef] 20. Młynarski, S.; Pilch, R.; Smolnik, M.; Szybka, J.; Wi ˛azania, G. A model of an adaptive strategy of preventive maintenance of complex technical objects. Eksploat. Niezawodn./Maint. Reliab. 2019,22, 35–41. [CrossRef] 21. White, J.R.; Widup, R. Factors to Consider When Determining Maintenance Intervals. IEEE Trans. Ind. Appl. 2014 ,50, 188–194. [CrossRef] 22. Coria, V.; Maximov, S.; Rivas-Dávalos, F.; Melchor, C.; Guardado, J. Analytical method for optimization of maintenance policy based on available system failure data. Reliab. Eng. Syst. Saf. 2015,135, 55–63. [CrossRef] 23. Zhao, X.; Cai, J.; Mizutani, S.; Nakagawa, T. Preventive replacement policies with time of operations, mission durations, minimal repairs and maintenance triggering approaches. J. Manuf. Syst. 2020,61, 819–829. [CrossRef] [PubMed] 24. Zhao, X.; Chen, M.; Nakagawa, T. Periodic replacement policies with shortage and excess costs. Ann. Oper. Res. 2020 . [CrossRef] 25. Bunea, C.; Bedford, T. The effect of model uncertainty on maintenance optimization. IEEE Trans. Reliabil. 2002 ,51, 486–493. [CrossRef] 26. Mazzuchi, T.A.; van Dorp, J.R. A Bayesian expert judgement model to determine lifetime distributions for maintenance optimisation. Struct. Infrastruct. Eng. 2012,8, 307–315. [CrossRef] 27. Zhu, J.; Wang, L.; Spiryagin, M. Control and decision strategy for a class of Markovian jump systems in failure prone manufacturing process. IET Control. Theory Appl. 2012,6, 1803–1811. [CrossRef] 28. van Dorp, J.R.; Mazzuchi, T.A. Three-point lifetime distribution elicitation for maintenance optimization in a bayesian context. In International Series in Operations Research & Management Science; Springer International Publishing: Cham, Switzerland, 2021; pp. 147–177. [CrossRef] 29. Fouladirad, M.; Paroissin, C.; Grall, A. Sensitivity of optimal replacement policies to lifetime parameter estimates. Eur. J. Oper. Res. 2018,266, 963–975. [CrossRef] 30. Wang, C.; Guo, J.; Shen, A. Sensitivity analysis of censoring data from component failure analysis and reliability evaluation for the aviation internet of things. Comput. Commun. 2020,157, 28–37. [CrossRef] 31. Hu, Q.; Yue, W. Optimal replacement of a system according toa semi-Markov decision process in a semi-Markov environment. Optim. Methods Softw. 2003,18, 181–196. [CrossRef] 32. Gualeni, P.; Perrera, F.; Raimondo, M.; Vairo, T. Accessibility for maintenance in the engine room: Development and application of a prediction tool for operational costs estimation. Ship Technol. Res. 2022. [CrossRef] 33. Wu, B.; Cui, L. Reliability analysis of periodically inspected systems with competing risks under Markovian environments. Comput. Ind. Eng. 2021,158, 107415. [CrossRef] 34. Tran, H.; Setunge, S.; Shi, L. Markov Chain–Based Inspection and Maintenance Model for Stormwater Pipes. J. Water Resour. Plan. Manag. 2021,147, 04021077. [CrossRef] 35. Wu, B.; Maya, B.I.G.; Limnios, N. Using semi-Markov chains to solve semi-Markov processes. Methodol. Comput. Appl. Probab. 2021,23, 1419–1431. [CrossRef] 36. Hu, J.; Shen, J.; Shen, L. Periodic preventive maintenance planning for systems working under a Markovian operating condition. Comput. Ind. Eng. 2020,142, 106291. [CrossRef] 37. Nie, R.; He, S.; Liu, F.; Luan, X. Sliding mode controller design for conic-type nonlinear semi-Markovian jumping systems of time-delayed Chua’s circuit. IEEE Trans. Syst. Man Cybern. Syst. 2019,51, 2467–2475. [CrossRef] 38. De Jonge, B.; Scarf, P.A. A review on maintenance optimization. Eur. J. Oper. Res. 2020,285, 805–824. [CrossRef] 39. Ruschel, E.; Santos, E.A.P.; Loures, E.d.F.R. Industrial maintenance decision-making: A systematic literature review. J. Manuf. Syst. 2017,45, 180–194. [CrossRef] 40. Li, D.; Li, Q.; Mingcheng, E.; Jiang, Z.; Ma, J. Failure analysis of coupler knuckle considering truncated and censored lifetime data. In Proceedings of the 2019 Prognostics and System Health Management Conference (PHM-Qingdao), Qingdao, China, 25–27 October 2019; IEEE: Piscataway, NJ, USA , 2019; pp. 1–5. 41. Taghipour, S.; Banjevic, D. Trend analysis of the power law process using Expectation–Maximization algorithm for data censored by inspection intervals. Reliab. Eng. Syst. Saf. 2011,96, 1340–1348. [CrossRef] 42. Sánchez-Herguedas, A.; Mena-Nieto, A.; Rodrigo-Muñoz, F. A new analytical method to optimise the preventive maintenance interval by using a semi-Markov process and z-transform with an application to marine diesel engines. Reliab. Eng. Syst. Saf. 2021,207, 107394. [CrossRef] 43. Wen, Y.; Rahman, M.F.; Xu, H.; Tseng, T.L.B. Recent Advances and Trends of Predictive Maintenance from Data-driven Machine Prognostics Perspective. Measurement 2021,187, 110276. [CrossRef] 44. Wu, B.; Cui, L.; Fang, C. Reliability analysis of semi-Markov systems with restriction on transition times. Reliab. Eng. Syst. Saf. 2019,190, 106516. [CrossRef] 45. Yi, H.; Cui, L.; Shen, J.; Li, Y. Stochastic properties and reliability measures of discrete-time semi-Markovian systems. Reliab. Eng. Syst. Saf. 2018,176, 162–173. [CrossRef] 46. Klein, J.P.; Moeschberger, M.L. Survival Analysis: Techniques for Censored and Truncated Data; Springer: Cham, Switzerland, 2003; Volume 1230. Sensors 2022,22, 1432 17 of 17 47. Shen, A.; Guo, J.; Wang, Z.; Jia, W. A novel reliability evaluation method on censored data. J. Mech. Sci. Technol. 2017 , 31, 1105–1117. [CrossRef] 48. Ahmed, E.A. Estimation of some lifetime parameters of generalized Gompertz distribution under progressively type-II censored data. Appl. Math. Model. 2015,39, 5567–5578. [CrossRef] 49. Teimouri, M.; Hoseini, S.M.; Nadarajah, S. Comparison of estimation methods for the Weibull distribution. Statistics 2013 , 47, 93–109. [CrossRef] 50. Genschel, U.; Meeker, W.Q. A Comparison of Maximum Likelihood and Median-Rank Regression for Weibull Estimation. Qual. Eng. 2010,22, 236–255. [CrossRef] 51. Ranjan, R.; Sen, R.; Upadhyay, S.K. Bayes analysis of some important lifetime models using MCMC based approaches when the observations are left truncated and right censored. Reliab. Eng. Syst. Saf. 2021,214, 107747. [CrossRef] 52. Balakrishnan, N.; Mitra, D. Left truncated and right censored Weibull data and likelihood inference with an illustration. Comput. Stat. Data Anal. 2012,56, 4011–4025. [CrossRef] 53. Ferreira, L.A.; Silva, J.L. Parameter estimation for Weibull distribution with right censored data using EM algorithm. Eksploat. Niezawodn./Maint. Reliab. 2017,19, 310–315. [CrossRef] 54. Joarder, A.; Krishna, H.; Kundu, D. Inferences on Weibull parameters with conventional type-I censoring. Comput. Stat. Data Anal. 2011,55, 1–11. [CrossRef] 55. Starling, J.K.; Mastrangelo, C.; Choe, Y. Improving Weibull distribution estimation for generalized Type I censored data using modified SMOTE. Reliab. Eng. Syst. Saf. 2021,211, 107505. [CrossRef] 56. Elmahdy, E.E. A new approach for Weibull modeling for reliability life data analysis. Appl. Math. Comput. 2015 ,250, 708–720. [CrossRef] 57. Zhang, L.; Xie, M.; Tang, L. Bias correction for the least squares estimator of Weibull shape parameter with complete and censored data. Reliab. Eng. Syst. Saf. 2006,91, 930–939. [CrossRef] 58. Musleh, R.M.; Helu, A. Estimation of the inverse Weibull distribution based on progressively censored data: Comparative study. Reliab. Eng. Syst. Saf. 2014,131, 216–227. [CrossRef] 59. Jia, X.; Wang, D.; Jiang, P.; Guo, B. Inference on the reliability of Weibull distribution with multiply Type-I censored data. Reliab. Eng. Syst. Saf. 2016,150, 171–181. [CrossRef] 60. Ducros, F.; Pamphile, P. Bayesian estimation of Weibull mixture in heavily censored data setting. Reliab. Eng. Syst. Saf. 2018 , 180, 453–462. [CrossRef] 61. Jabeen, R.; Ahmad, A.; Feroze, N.; Gilani, G.M. Estimation of location and scale parameters of Weibull distribution using generalized order statistics under type II singly and doubly censored data. Int. J. Adv. Sci. Technol. 2013,55, 67–80. 62. Sánchez Herguedas, A.; Crespo Márquez, A.; Rodrigo Muñoz, F. Optimizing preventive maintenance over a finite planning horizon in a semi-Markov framework. IMA J. Manag. Math. 2020,33, 75–99. [CrossRef]