A method for obtaining the preventive maintenance interval in the absence of failure time data

Ek s p l o a t a c j a i NiEzawodNosc – Ma i N t E N a N c E a N d REliabi lit y Vo l . 24, No. 3, 2022 564 Eksploatacja i Niezawodnosc – Maintenance and Reliability Volume 24 (2022), Issue 3 journal homepage: http://www.ein.org.pl Indexed by: (*) Corresponding author. E-mail addresses: 1. Background The maritime transport sector must face two significant challenges: the ecological transformation helping to comply with international agreements to reduce emissions and the digital transformation derived from the development of Industry 4.0. In both cases, maintenance management can play a crucial role. Regarding the first of them, maritime transport was responsible for more than 2.89% (1,076 “million tonnes”) of the world’s total greenhouse gas emissions in 2018, but it will continue to increase in the future due to the globalization of markets and to the growth of international trade [18]. In maritime transport, components installed on ships are subject to hours of operation where maintenance interventions are interspersed to ensure the desired operating condition or restore it to that condition in case of failure. The maintenance manager must ensure that their equipment provides optimal system reliability, availability and security values, keeping maintenance costs as low as possible to help maintain profitability for the shipowner’s business [2]. One of the ways to reduce greenhouse gas emissions and other polluting gases caused by ships is to improve their maintenance operations through their life cycle. The maintenance manager usually does not modify the preventive intervals that the equipment manufacturer has designed to reduce the failure. Conditions of use and maintenance often change from design conditions. In these cases, continuing using the manufacturer's preventive intervals can lead to non-optimal management situations. This article proposes a new method to calculate the preventive interval when the hours of failure of the assets are unavailable. Two scenarios were created to test the effectiveness and usefulness of this new method, one without the failure hours and the other with the failure hours corresponding to a bypass valve installed in the engine of a maritime transport surveillance vessel. In an easy and fast way, the proposed method allows the maintenance manager to calculate the preventive interval of equipment that does not have installed an instrument for measuring operating hours installed. Highlights Abstract A new method to find the preventive interval opti-• mizes the asset's costs and incomes. This method contains low mathematical complex-• ity and rapid execution. The Weibull distribution and a semi-Markovian • model are used to compute the returns. Use of failure data when the time values for such • failures are unknown. Apply to those assets that do not have an instru-• ment that measures the operating hours. A method for obtaining the preventive maintenance interval in the absence of failure time data Antonio Sánchez-Herguedas a, Ángel Mena-Nieto b,*, Francisco Rodrigo-Muñoz c a University of Seville, Department of Industrial Management I, School of Engineering, Camino de los Descubrimientos s/n, 41092 Seville, Spain b University of Huelva, Department of Electrical and Thermal Engineering, Design and Projects, School of Engineering, Centre for Advanced Studies in Physics, Mathematics and Computation, Campus El Carmen, 21071 Huelva, Spain c University of Seville, Department of Applied Mathematics II, School of Engineering, Camino de los Descubrimientos s/n, 41092 Seville, Spain Sánchez-Herguedas A, Mena-Nieto A, Rodrigo-Muñoz F. A method for obtaining the preventive maintenance interval in the absence of failure time data. Eksploatacja i Niezawodnosc – Maintenance and Reliability 2022; 24 (3): 564–573, http://doi.org/10.17531/ ein.2022.3.17. Article citation info: ship maintenance, preventive interval, maintenance model, semi-Markovian process, incomplete data, Weibull function, maintenance costs. Keywords This is an open access article under the CC BY license (https://creativecommons.org/licenses/by/4.0/) Therefore, any improvement in the maintenance of ships will contribute to reducing their emissions of pollutants and greenhouse gases. Moreover, as in other economic sectors, the current digital transformation is changing how maintenance is managed [38], [32]. New maintenance strategies are increasing the reliability of assets during their lifecycle, reducing the frequency of preventive maintenance [9]. Maintenance is key to the development of Industry 4.0 [7]. The massive use of intelligent sensors allows an increasing amount of data to be obtained, which must be analysed efficiently and effectively to support decision-making and increasingly complex technological systems management [30]. Focusing now on the case of ship maintenance, sometimes, the maintenance manager finds himself with assets that suffer repetitive failures even when the maintenance strategy set by the manufacturer has been followed. In these cases, the operating and maintenance conditions of these assets probably differ from the conditions for which they were designed. The modification of the maintenance intervals A. Sanchez-Herguedas (ORCID: 0000-0001-5135-3250): [email protected], A. Mena-Nieto (ORCID: 0000-0002-0828-0612): [email protected], F. Rodrigo-Munoz (ORCID: 0000-0001-5350-9378): [email protected] Ek s p l o a t a c j a i Ni E z a w o d N o s c – Ma i N t E N a N c E a N d REliabilit y Vo l . 24, No. 3, 2022 565 according to each asset’s particular conditions of use must be dealt with by the person in charge of maintenance [45]. Available information is an essential factor to consider. To make an accurate analysis it is not only necessary to know the number of failures over time, and the time at which failures occur. It is also necessary to consider other factors, such as the number of cycles of preventive interventions performed, or the product usage that characterise the ageing of the asset. Recently, methods using this information have been proposed, based on collected parameters and machine learning. However, when this data is not available, only the number of failures that have occurred in a time interval, we are faced a different scenario [20, 21]. This work shows how to optimize the periodic preventive maintenance interval of a ship engine exhaust component when the available information on the failure mode understudy is minimal. This kind of problem with a scarcity of data is a more common type of problem than one might think in the professional practice of maintenance engineers in many facilities, equipment and fleets of land or sea transport. Despite the current proliferation of sensors in companies, there are situations, especially in small companies and industries, where, for economic or other reasons, they do not record sufficient historical failure data on the assets they own, so they must manage the maintenance of those assets and save the maintenance costs with little or no data, beyond that provided by their manufacturers. There are not so many studies that advise maintenance engineers on what maintenance interval to adopt when there is very little historical data on equipment failures. They need to figure out how to predict when the next failure will occur and how it will occur. This challenge is more difficult when there is no historical failure data. In other words, they still need to know how to proceed with very little data regarding asset failures in the data-driven methods era. However, the scientific literature on how to address this problem is scarcer than would be desirable, so this paper responds to that gap. A sample of recent work is present below, which has proposed mathematical models to address situations with incomplete data. Valis et al. [40] present new single and multiple error state-space models to model and predict the reliability of water mains from patchy and sparse failure data. Nobakhti et al. [27] propose a hybrid approach using fault tree analysis and the Mamdani fuzzy inference to obtain the reliability response as a function of a few frequently operating pressure and temperature, which avoid the lack of historical data. Andrzejczak et al. [4] study the problem of the lack of data in new technical facilities, presenting a method, with a Bayesian approach, for estimating the probability distribution of the lifetime for these assets based on expert assessments of three parameters characterizing the expected lifetime of these assets. Furthermore, they show some practical applications using the Weibull distribution. Yun-Fei et al. [26] establish a hidden semi-Markovian model that clusters the incomplete degradation data. The model allows predictions for the remaining useful life and is applied to bearing failure data. Zhang et al. [47] developed a Markov chain Monte Carlo method to perform multiple imputations for incomplete correlated ordinal data using the multivariate probit model and simulations to compare the performance of their method with other methods for various missing data scenarios. Li et al. [23] tackle fault detection and diagnosis from incomplete data. As a preliminary step, they propose the adjacent information recovery filter to recover the missing data from sensors. This filter considers the time series adjacency information through a hidden Markov model. Yuguang et al. [46] present a method based on a Markov chain model and generalized projection non-negative matrix factorization to detect and diagnose faults in industrial processes. This method is applicable when there are missing data in incomplete measurements since they often are correlated with some of the available variables. Aguirre- Salado et al. [1] analyse the maximum intensity levels of earthquakes with incomplete data on the coast of southern Mexico using a random censorship approach. They use a flexible semi-parametric Bayesian approach, whose parameters are estimated through the Markov chain Monte Carlo. In the work of Lupton and Allwood [25], a general procedure based on a Bayesian approach is developed for its application to material flow analysis. Using Markov Chain Monte Carlo simulations, the procedure uses incomplete or missing data to map global steel production. Liu et al. [24] apply the delay time theory to assess the reliability of a system with sufficient inspection data but insufficient failure data. First, they develop an optimization model for individual components that minimizes maintenance, failure, and downtime costs over the component lifecycle and then extend it to multi-component serial systems, applying it to a locomotive. Yamany and Abraham [43] develop and validate a non-homogeneous Markovian pavement performance model that improves on previously proposed probabilistic pavement performance models because it allows addressing cases where there are no historical preventative maintenance data. Finally, one of the most recent works related to our research is Zhao et al. [49]. Using a study case from the petrochemical industry, they apply a Bayesian framework to update unknown parameters in a Wiener degradation model to overcome the practical challenges of information sparsity. The latest trend to address the lack of data is to use artificial intelligence techniques; see, for example, references [3], [14] and [50]. The main objective of this paper is to optimize the periodic preventive maintenance interval of a ship engine exhaust component when the available information on the failure mode understudy is minimal. For that purpose, the results obtained under two scenarios are compared. In the first scenario A, only information regarding the number of failures that have occurred during a period of time is available, while for the second scenario B, the hours of operation at which the failures have occurred are also available. The active component to be analysed presents a failure mode due to the accumulation of combustion residues on the sealing surface and the axis of the second turbo actuation bypass valve. The failure in the closure causes the engine malfunction, emitting exhaust gases with a high content of polluting particles, which forces the equipment to stop for cleaning and replacement, and therefore the cessation of activity and its incomes. This article uses the Benard’s approximation as an estimator to determine the observed failure distribution function [28], However, in each scenario, a different method is used to determine the hours at which the failures occur. Scenario A proposes a method based on the total hours of operation and the number of failures counted. In scenario B, the failure hours are taken as problem data. The process continues in parallel for the two scenarios by calculating the Weibull cumulative distribution [39]. Some authors use other methods to calculate this distribution [4]. In both cases, the curve is fitted, minimizing the mean square error (MSE) [12]. Once the distribution function has been theorized, it is necessary to use a mathematical model that simulates the evolution over time of the activity of the element under study (the bypass valve). Traditionally, Markov processes have been used to represent the evolution over time of industrial components [16], since they allow different states of activity to be established (operation, corrective, preventive, etc.) through which the element can go through its life cycle. These models are very useful for establishing behaviour because they demarcate history from the future and establish a law of transition probabilities between states. However, these models use constant failure rates over time [41]. Carrying out preventive maintenance tasks on the equipment is justified because they present increasing failure rates. It forces the use of semi-Markovian type models that preserve the advantages of Markov processes and allow the introduction of variable failure rates [22], [42] and [31]. To calculate the preventive interval, we use a method developed from a semi-Markovian process with three possible states (operational, corrective and preventive), using the direct costs of corrective and preventive maintenance tasks and other costs associated with changes in statements, and the income obtained from the use of the asset [34] . The rest of the document is organized as follows: Section 2 contains the description of the method followed to reach the value of the preventive interval. The entire section is divided into five parts. In Ek s p l o a t a c j a i NiEzawodNosc – Ma i N t E N a N c E a N d REliabi lit y Vo l . 24, No. 3, 2022 566 Section 2.1, the failure data and the administrative information regarding the preventive, corrective intervention and operation are selected, particularly for a case of an actual ship engine. In Section 2.2, the procedure for determining the failure hours when these are not available is applied, and the observed distribution function is calculated. In Section 2.3, the estimation of the theoretical distribution functions (Weibull functions) is carried out. In Section 2.4, the three-state semi-Markovian model that governs the behaviour over time of the failure behaviour of the asset under study (bypass valve) is presented. In Section 2.5, the method for calculating the average accumulated return for each transition between states is shown. In Section 2.6, the mathematical formula that provides the value of the preventive interval that optimizes the average accumulated return is shown. Section 3 presents the results for each scenario of the preventive interval for the different transitions. In Section 4, the results obtained are discussed, and the method is proposed as a tool to help the maintenance manager calculate the preventive interval when the hours of failure of the asset are not available. Finally, Section 5 presents the conclusions of this work. 2. Material and methods Due to price, most physical assets do not have a measuring instrument that allows a chronology of the events that happen to them (a clock or hour meter). For this reason, Computer-Aided Maintenance Management Systems (CAMMS) lack this information. When trying to use failure data to make calculations and optimizations, this can be an inconvenience or one of the aggravating factors that prevent satisfactory results [19], [29], [37] and [48]. For the development of this article, the failure data will be used together with its failure appearance times, comparing the preventive interval obtained after applying the methodology with the interval obtained when only the failure data without appearance hours are known. A procedure will be used to obtain a series of failure hours in scenario A and, on the other hand, the Benard’s approximation formula will be used to calculate the observed failure distribution function, which is common to both scenarios. It is then unified for the two scenarios from this point to reach the two theoretical distribution functions, Weibull functions [6], one for each scenario. From here, a semi-Markovian model is established that simulates the evolution in time, accounting for the transitions between the states [15]. The model allows one to determine the average accumulated return for each of the transitions. For this, it requires information of an administrative nature: costs of maintenance interventions and income from the operation of the equipment. This average accumulated return depends on the preventive interval, so it is possible to obtain the preventive interval that optimizes the average accumulated return. The objective of our work is to compare the size of the optimal preventive interval for the two established scenarios and to verify that the method followed when we only have the number of failures is a method whose results are very close to the results achieved when the same number of failures are also used. This article is based on [34] , where other distributions (Exponential, Log-normal and Normal) have been analysed, concluding that the one that best fits the failure times is the Weibull distribution. The use of the Weibull distribution for the analysis of failure times is now widespread. See [5], [6] and [44]. In the pioneering paper [xx], failure data is analysed for the first time. The author determines that the best fitting distribution is the exponential distribution. However, this is a particular case of the Weibull distribution which, because it has more parameters, is capable of covering more general situations. Moreover, the exponential distribution, due to its lack of memory, does not allow preventive maintenance, only corrective maintenance, which would render this article meaningless. 2.1. Real case. Data selection and process information The physical asset under study is a bypass valve installed in marine diesel engines. Its mission is to control the flow of exhaust gases, determining the volumetric configuration where they will expand, in an initial configuration made up of one turbo or in a second, made up of the two turbos. This valve is essential if one wants a quick increase in rpm and therefore the power delivered. Its failure produces a significant emission of exhaust gases. During a period of 15 months and two weeks, valve failure data corresponding to 5 engines have been collected. A total of 16 bugs have been recorded, while the total hours that the five engines have worked is 26,400 hours. As the patrol boats have an hour meter, it has been possible to record the exact times in which the failures have occurred. The failure times already ordered from lowest to highest are presented in Table 1. From the information of the maintenance interventions carried out, it has been possible to establish that the duration of the corrective repair task to solve the fault has an average of 8 hours, while the average of the preventive task carried out to avoid the failure is 5 hours. The technical personnel who carry out both interventions have different categories, so the cost of the corrective intervention is €95/hour, while that of the preventive intervention is €82/hour. The cost borne by the customer after the failure occurs is valued at €3,270, while the cost borne after deciding to carry out the preventive intervention is €1. The costs of spare parts used in both interventions coincide at €360. The owner’s income for the use of the equipment is €5/hour. 2.2. Estimation of the observed failure distribution function The observed failure distribution function is the function that is estimated from the collected failure data. To develop it, the Benard’s approximation is used as an estimator, ( ) ( ) 0.3 / 0.4 i FiN=−+ , where i represents the order of the failure and N the total number of failures (in this case 16) . The failure distribution function observed is the same for the two scenarios described, since only the number of failures is involved in its calculation. The information available in scenario A is the 16 failures in the 26,400 hours and the information for scenario B is shown in Table 1. The application of Benard’s approximation requires establishing an order in the failures. Place the youngest bug in order 1 and the oldest bug in order 16. Table 2 shows these values. In order to determine the theoretical distribution function and the observed failure distribution function, it is necessary to know the number of hours that the asset was operating before the failure occurred. For scenario B, this information is collected in Table 1. However, this table of values must be constructed for scenario A based on the available information. The procedure to follow is as follows. From the equipment’s total operating hours (26,400 hours), the • average occurrence of failures is calculated, 1,650 hours, dividing 26,400 hours by the 16 failures. The interval between failures is established, the average appear-• ance of failures divided by the number of failures, giving it a value of 1650/16 = 104 hours. The failures occurrence is distributed around the mean, with each • failure occurrence increasing or decreasing by the amount of interval between failures. If the number of failures is even, the two central events (failures • 8 and 9) are placed, equidistant from the mean separated by the interval. The values obtained are shown in Table 3. In this case, since it is an even number of failures, failure 8 corresponds to the value 1.650-52 = 1.598. Failure 9 corresponds to the value 1.650+52 = 1.702. 2.3. Estimation of the theoretical distribution functions From the failure data for each scenario, Table 1 and Table 3, and the observed failure distribution function data, Table 2, the theoretical failure distribution function that best fits must be estimated for each scenario for the pair (hours of failure, function observed). In the case Ek s p l o a t a c j a i Ni E z a w o d N o s c – Ma i N t E N a N c E a N d REliabilit y Vo l . 24, No. 3, 2022 567 of physical assets, it is usually adjusted to the Weibull distribution function [5]. Scenario A We first try to fit the observed function to the two-parameter Weibull αβ ,� () , using columns 5 and 4 of Table 4 (columns 1, 3 and 2 correspond to Table 2 and Table 3). To do this, we represent the points of the pair (hours of failure, observed function) in the graph on the left of Fig. 1, using the logarithmic scales, ( ) ( ) l n 1 / 1 — i ln F on the vertical axis and i lnt on the horizontal axis. Excel allows you to plot the trend line that best fits the points using the method of least squares. The form parameter α of the Weibull function coincides with the slope of the trend line, while the scale parameter β corresponds to the negative exponential of the quotient between the ordinate at the origin and the slope [34]. Table 1. An ordered list of hours of failures in the valves of the five engines during a 26,400 hour-period. Scenario B Failure, i 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 Fi 0.0427 0.1037 0.1646 0.2256 0.2866 0.3476 0.4085 0.4695 0.5305 0.5915 0.6524 0.7134 0.7744 0.8354 0.8963 0.9573 Table 2. Values of the observed failure distribution function obtained from Benard's approximation. Scenarios A and B Failure, i 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 Fi 0.0427 0.1037 0.1646 0.2256 0.2866 0.3476 0.4085 0.4695 0.5305 0.5915 0.6524 0.7134 0.7744 0.8354 0.8963 0.9573 Table 3. List of failure hours built for scenario A. Failure 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 Hours 870 974 1,078 1,182 1,286 1,390 1,494 1,598 1,702 1,806 1,910 2,014 2,118 2,222 2,32 2,430 Table 4. Data for the construction of the graphs in Fig. 1 Failure Failure hours, tiObserved function, Filn tiln ln (1/(1−Fi (t))) ln(ti−γ) 1 870 0.04268 6.768493 −3.132225 5.799093 2 974 0.10366 6.881411 −2.212435 6.073045 3 1,078 0.16463 6.982863 −1.715435 6.287859 4 1,182 0.22561 7.074963 −1.363831 6.464588 5 1,286 0.28659 7.159292 −1.085620 6.614726 6 1,390 0.34756 7.237059 −0.850883 6.745236 7 1,494 0.40854 7.309212 −0.644061 6.860664 8 1,598 0.46951 7.376508 −0.455772 6.964136 9 1,702 0.53049 7.439559 −0.279633 7.057898 10 1,806 0.59146 7.498870 −0.110737 7.143618 11 1,910 0.65244 7.554859 0.055260 7.222566 12 2,014 0.71341 7.607878 0.222919 7.295735 13 2,118 0.77439 7.658228 0.398070 7.363914 14 2,222 0.83537 7.706163 0.590023 7.427739 15 2,326 0.89634 7.751905 0.818304 7.487734 16 2,430 0.95732 7.795647 1.148658 7.544332 Fig. 1. Graphic representation of failure hours and observed function in logarithmic coordinates. Trend lines and curves. Scenario A Ek s p l o a t a c j a i NiEzawodNosc – Ma i N t E N a N c E a N d REliabi lit y Vo l . 24, No. 3, 2022 568 From the fitted line obtained for the linear trend ( 3.5706 26.828yx=− ), the values of α =357. and β =1 832, . However, the fitted curve of order 2 yxx=− +− () 0 9314 17 177 76 429 2 ... shows that a guaranteed life γ (location parameter), can exist since the coefficient of 2 x is negative, graph on left of Fig. 1. To find the value of γ, we must move the origin ()ti− γ until the coefficient of 2 x is zero. The points of the new pair (failure hours minus the guaranteed life, observed function) are re-plotted on the graph on the right of Fig. 1. The logarithmic scales ln ( ) ( ) ln 1 / 1 — i ln F and ln ln ( ) i tγ− are now used, and the values in columns 5 and 6 of Table 4. Values are given to γ until the coefficient of 2 x is zero, and the fitting curve becomes a straight line yxx=− +− () 0 0003 2 1824 15 578 2 ... .From this line, the 2.18α= and 1,266β= , values are obtained, which, with the adjusted value for 540γ= , are the three parameters of the Weibull function. Scenario B In this scenario, we know the times at which the failures occurred, Table 1. These data are collected in column 2 of Table 5 and modify the values of columns 4 and 6 that vary concerning those of Table 4. Following the same procedure described for scenario A, the graphs in Fig. 2 are constructed from the data in Table 5. Again values are given to γ until the coefficient of 2 x is zero, and the fitting curve is flattened. From this straight line, the values of 2.12α= . and 1,028β= are obtained, which together with the adjusted value for 756γ= constitute the three parameters of the Weibull function for scenario B. 2.4. Semi-Markovian model of evolution of the system in time The calculation of the preventive interval of the bypass valve requires developing a mathematical model that reflects the behaviour of the asset over time. This model should allow the evolution of the asset between three states (operational S1, corrective S2 and preventive S3). You can use stochastic discrete or continuous-time processes; an introductory class between them is the Markov process [17]. This process meets the Markovian property, “the future only depends on the present, not the past”. In the discrete-time model (chain of Markov), only transitions between states and their probabilities are taken into account. The permanence time in each state is not relevant. Due to the Markovian property, the permanence time is a random variable with exponential distribution in the continuous-time model. However, the physical asset object of study presents an increase in the failure rate as the operating time increases. This forces one to discard Markovian models and opt for semi-Markovian processes. This type of model dif- Fig. 2. Graphic representation of failure hours and observed function in logarithmic coordinates. Trend lines and curves. Scenario B Table 5. Data for the construction of the graphs of Fig. 2 Failure Failure hours, tiObserved function, Fi ln ti ln ln Ft i 11/− () () () ln ti ()− γ 1 991 0.04268 6.898715 −3.132225 5.456175 2 1,082 0.10366 6.986566 −2.212435 5.784440 3 1,315 0.16463 7.181592 −1.715435 6.324717 4 1,342 0.22561 7.201916 −1.363831 6.371954 5 1,405 0.28659 7.247793 −1.085620 6.474199 6 1,515 0.34756 7.323171 −0.850883 6.630947 7 1,520 0.40854 7.326466 −0.644061 6.637520 8 1,570 0.46951 7.358831 −0.455772 6.700977 9 1,592 0.53049 7.372746 −0.279633 6.727671 10 1,635 0.59146 7.399398 −0.110737 6.777874 11 1,769 0.65244 7.478170 0.055260 6.919881 12 1,797 0.71341 7.493874 0.222919 6.947168 13 1,837 0.77439 7.515889 0.398070 6.984901 14 2,177 0.83537 7.685703 0.590023 7.258553 15 2,433 0.89634 7.796880 0.818304 7.424285 16 2,536 0.95732 7.838343 1.148658 7.483919 Ek s p l o a t a c j a i Ni E z a w o d N o s c – Ma i N t E N a N c E a N d REliabilit y Vo l . 24, No. 3, 2022 569 fers from Markovians in that the time of permanence in each state does not follow an exponential distribution [13], which implies that semi- Markovian models do not meet the Markovian property. However, the successive transitions between states form a chain of Markov, called the Markov chain embedded in the Semi-Markovian model [8]. The model evolves over time, Fig. 3, accumulating returns (income in S1 and costs in S2 and S3) due to the times of permanence in each state and transitions between states. To this end, three types of square order matrices are developed, the transition probability matrix between states [ ] P , the matrix of permanence times in each state [ ] Q and the matrix of returns by permanence in each state and the transition to the next [ ] R . In this case, the non-null elements of the matrix [ ] P are 12 13 21 31 , , pppyp . The non-null elements of the matrix [ ] Q are 12 13 21 31 , , q q qyq . The non-null elements of the matrix [ ] R are 12 13 21 31 , , rrryr . If the asset fails before reaching the time of the preventive interval τ , the asset passes to the corrective state. After a time 21 q , the asset returns to the operating state. If during the operation, time τ is reached, the asset passes to the preventive state. After a time 31 q , the asset returns to the operating state. The element of the matrix [ ] P , 12 p , is the probability of failing the asset before reaching preventive maintenance. It is defined by the value that reaches the failure distribution function in time τ . The element 13 p is the probability of achieving preventive maintenance. It is defined as 12 1p− . The elements 21 31 pyp take value 1 since the asset always returns to the operating state after the corrective or preventive task. The 12 q element is the average time that the asset remains operative before failing. The 13 q element is the average time that the asset state before moving on to the preventive state 13 13 1 13 rqRR=⋅+ . It is also composed of two concepts, the income per operating time 1 Rτ⋅ and the cost motivated by the activation of the preventive 13 R . The 21 r element is the average return that demands the asset when it remains in a corrective state before moving to the operational state 21 21 2 21 rqRR=⋅+ . It comprises two concepts, the cost per repair time 21 2 qR⋅ . (being 2 R the cost per hour of repair). and the cost motivated by activating the operating state 21 R . The element 31 r , is the average return demanded by the asset when it remains in the preventive state before going to the operating state 31 31 3 31 rqRR=⋅+ . It is also composed of two concepts, the cost per time dedicated to the preventive 31 3 qR⋅ . (being 3 R the cost per hour of preventive task) and the cost motivated by activating the operating state 31 R . 2.5. Calculation of the average accumulated return When the valve is operating, income is obtained, but when it is subjected to corrective tasks after a failure or preventive tasks after a certain time operating, expenses occur. The return in each transition is a random variable, so we cannot calculate it. However, it is possible to calculate the average accumulated return ( ) i vm in m transitions, starting from state i . In the first transition, assuming that we start from state S1, the average return in the first transition can be expressed as follows ( ) 3 111 1 1 · jj j vrp = =∑ . Once the asset is found in the j (S2 or S3) state, you can perform 1m− transitions and will accumulate returns in all of them. The average accumulated return in those 1m− transitions is a random variable, denoted by ( ) 1 j vm− , which can take values ( ) ( ) ( ) 123 1 , 1 , 1vmvmvm−−− with probabilities 123 , , jjj ppp that remain constant throughout the 1m− transitions, since the process is homogeneous. Its value can be calculated as: ( ) 3 1 1 · j ij j vmp = − ∑ . The average accumulated return in m transitions can be expressed as the sum of the returns of the first transition and the remaining 1m− transitions, according to the expression: ( ) ( ) ( ) 3 11 1 1 1 j ij j vmv vmp = =+ −⋅ ∑ . Similarly, expressions for ( ) 2 vm and ( ) 3 vm can be calculated [35]. It is a system of difference equations [36]. This system of equations can be resolved by z-transform [11], obtaining the equation Eq. ((1)): Fig. 3. Description of the transition process between states and the accumulation of returns associated with permanence and transition between states remains operative before passing to preventive, that is, the preventive interval τ . The 21 q element is the average time that the asset remains in a corrective state. We give this element the value of the average of the distribution function of the repair times, ( ) c Gt . The 31 q element is the average time that the asset remains in a preventive state. We consider this value the average distribution function of preventive intervention times, ( ) p Ht . The value taken by the 12 q element is conditioned that time τ is not reached and expressed in the form: ( ) ( ) 12 0 1 · q t f t dt F τ τ =∫ . The 12 r element is the average return that delivers the asset when it remains in an operating state before moving to the corrective state 12 12 1 12 rqRR=⋅+ . It is composed of two concepts, the operating time income 12 1 qR⋅ . (being 1 R the income per hour of operation) and the cost motivated by the appearance of the failure 12 R . The 13 r element is the average return that delivers the asset when it remains in an operating The equation obtained depends on m and τ , and uses the distribution functions of failure times ( ) Ft , repair time ( ) c Gt and preventive time ( ) p Ht and their respective density functions ( ) ft , ( )) c gt and ( ) p ht . In the same way, the average accumulated returns are calculated starting from the S2 and S3 states. 2.6. Calculation of the preventive interval Starting from the equation obtained for ( ) 1 vm , the derivative is calculated with respect to τ and it equals zero ( ) ( ) 1 /0dv m dτ= . Then one proceeds to replace the failure distribution function F(t) by the Weibull function of three parameters and the other two distribution functions of time in the corrective and preventive states due to its average values [33]. The next step is cleared, the value of the preventive interval 0 τ that maximizes the average accumulated return, ( ) 1 vm , for each transition m , obtaining la Eq. (2): vm mRtf tdtRFRRF m 1 1 1 0 12 113 1 42111 () =++− () () () + () ++ − () − ∫ τ ττ τ () (()        +−−− () () () +         () −∞ ∫ 2111 2 0 21 mRtgtdtRF m ccc � τ ++ () +        − () ()                 ∞ ∫ Rtht dt RF ppp3 0 31 1� τ (1) Ek s p l o a t a c j a i NiEzawodNosc – Ma i N t E N a N c E a N d REliabi lit y Vo l . 24, No. 3, 2022 570 τγ β α αα 0 11 12 13 1 122121 211 211 − () =⋅ − −+ −−− () ++− () + − − − R RR m m Rq R m m−−− () Rq R 33131 (2) This Eq. (2) depends solely on the data of the case to be analysed and the transition m . 3. Results To compare the two scenarios, it is necessary to define the data required by Eq. (2). All the data are the same for both scenarios except the distribution functions of the failures times. Both functions were calculated in section 2.3. For scenario A a Weibull distribution is used (2.18, 1,266, 540), a Weibull (2.12, 1,028, 756) is used for Scenario B. The other data needed to find the optimal interval are set out in Table 6. The results found for the first ten transitions in scenario A are collected in Fig. 4. The results found for the first ten transitions on Scenario B are collected in Fig. 5. 4. Discussion At the beginning of the article, the possibility of establishing an acceptable value for the preventive interval was raised when the asset failure hours were not available. The only available information focused on the number of failures that occur in a certain number of hours of operation. In order to make a comparison, two scenarios were established. The first scenario, A, describes the proposed problem. A second scenario B includes the data corresponding to the asset failure hours under study. In both cases, the same failure mode is studied, omitting failure times for scenario A which are included for scenario B. Scenario A proposes a model for the distribution of failure hours, while in scenario B, the distribution of failure hours is taken directly from the starting data. The procedure followed for both scenarios is the same, except for the construction phase of the observed failure distribution function. This phase ends similarly for both scenarios, applying Benard’s approximation to obtaining the observed failure distribution function. The theoretical distribution functions are obtained; it is a tri-para- metric Weibull in both scenarios. The semi-Markovian model calculates the preventive interval that maximizes the average accumulated return in the following phase. Scenario A obtains an optimal preventive interval with a value 1,592 A τ= hours, with a medium accumulated return at the end of 10 transitions of ( ) 110 22,564 € A v= . In scenario B an optimal preventive interval is obtained with a value 1,475 B τ= hours, with an average accumulated return at the end of 10 transitions of ( ) 110 23,913 € B v= . Analysing these results, it is observed that the preventive interval in scenario A is only 8% higher than the preventive interval in scenario B. On the other hand, the average accumulated return after ten transitions decreases by 5.6%, supposing an excellent approximation. As close if we speak in terms of preventive maintenance intervals, these values indicate that the procedure used for obtaining schedules values of failures allows considerably good results to set up the value of the preventive interval when we do not have the failure times data. This methodology is applicable when failures are uniformly distributed around the mean operating time between failures (MTOBF). The proposed dispersion foresees the existence of an optimal preventive interval. Both values are known to maintenance engineers for equipment showing wear and tear. For those where the failure occurs randomly, the calculation of the preventive interval does not make sense. To analyse the effect of the dispersion of the data, we can compare the results obtained by varying the shape parameter of the function Weibull. Table 7 shows the relationship between α and τ0 if the other parameters and the data are held constant. As expected, for values of α ≥ 2 the values of τ0 approach a limit value around 1,500 hours. On the other hand, when failures occur randomly, the preventive interval increases considerably. These values correspond to this case, if other data are used, the results will be different, but there is always a limit value and the preventive interval will increase as the randomness of the failures increases (decrease α). This is important when only the number of Table 6. Data needed to find the preventive interval τ0 Failure Time Distribution Function Distrib. Function Repair Times Distrib. Function Preventive Times Weibull (α, β, γ) Normal (μ1, σ1)Normal (μ2, σ2) α β γ μ1μ2 2.18 / 2.12 1,267 / 1,028 540 / 756 8 5 Operational Returns Corrective Returns Preventive Returns R1 (€/hour) R12 (€) R13 (€) R2 (€/hour) R21 (€) R3 (€/hour) R31 (€) 5.0 −3,270 −1.0 −95.0 −360 −82.0 −360 Fig. 4. Graphic representation and values of the evolution of the preventive interval for home transition. Scenario A Table 7. Variation of the optimal preventive interval when changing the dispersion of failures (α, shape parameter of the Weibull) α1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 4.0 τ0 (hour) 2,755 2,012 1,764 1.649 1.587 1,551 1,529 1,516 1,507 1,503 1.502 Ek s p l o a t a c j a i Ni E z a w o d N o s c – Ma i N t E N a N c E a N d REliabilit y Vo l . 24, No. 3, 2022 571 failures is available. The maintenance engineer must estimate the random component of the failures. Finally, the results obtained with the method presented in the paper have been compared with the utilization of Monte Carlo simulation. A dynamic simulation model has been built using continuous time stochastic simulation [10]. States, transitions, sojourn times and returns have been defined. The variability is introduced by three seeds, for the preventive time, for the corrective time and for the random number generating failures according to the Weibull distribution. Difference equations for average accumulated return are used in the model and the simulation time step considered is one hour. 200 simulations were done per scenario A ( 1,592 hoursτ= ), and scenario B ( 1,475 hoursτ= ). A total number of 30,000 hours of simulation were considered. The results from this simulation are shown in Table 8 and can be compared with the ones presented in the last column from the analytical method. In Fig. 6, a graph containing the sensitivity results for scenario A is presented. The selected variable is returned after 10 transitions, and the confidence bounds illustrate the percentile values for the total 200 simulations over the time horizon simulated (30,000 hrs). The mean value plot is included in red (ending in 22,564, as presented in Table 8). 5. Conclusions The study’s objective presented in this article is to develop a procedure to find the preventive interval that optimizes the returns obtained by the operation and maintenance of an asset. This method has been designed to be applied when the time values for when such failures occur are unknown. This is very typical of those assets that do not have an instrument that measures the operating hours of the asset. The only requirement is to know the number of failures that have occurred during a period of time. The procedure requires the establishment of a distribution of possible failure hours and, on the other hand, determining the observed failure distribution function using Bernard’s approximation and theoretical distribution functions. A semi-Markovian model, which evolves over time through transitions between three states, is then used to compute the average accumulated return. From there, the preventive interval that maximizes the average accumulated return is calculated. To verify that the procedure followed to establish the distribution of the failure hours is correct, the same procedure has been applied to a second scenario where the failure hours are known. The result has been very similar, which validates this procedure as a method to reach the preventive interval when the hours of failure of the assets are not available. The calculation of the preventive interval for maintenance managers is a task that is often not carried out, and the value imposed by the manufacturer is trusted, even if the operating and maintenance conditions differ from those of design, even more so in the case where data on failure hours is not even available. The described procedure is easy to apply and within reach of any maintenance manager due to its low mathematical complexity and its rapid execution. This procedure can become a tool to help asset management. Acknowledgments This research was funded by the “Consejería de Transformación Económica, Industria, Conocimiento y Universidades de la Junta de Andalucía” and by the European Regional Development Fund (FEDER), within the framework of the FEDER of Andalusia 2014-2020 Project Reference UHU-202031, and by the grant Ref. RTI2018-094614-B-I00 (SMASHING Project) into the “Programa Estatal de I+D+i Orientada a los Retos de la Sociedad” funded by MCIN/AEI/ 10.13039/497 501100011033. Fig. 5. Graphic representation and values of the evolution of the preventive interval for each transition. Scenario B Table 8. Monte Carlo simulation. Returns values v1 (10) for scenarios A and B Variable Count Min Max Mean Median StDev Analytical Method Returs for m=10 Scenario A, τ=1,592 200 3,246 39,069 22,564 22,677 7.696 22,564 Scenario B, τ=1,475 200 8,318 33,999 23,913 24,124 5.758 23,913 Fig. 6. Sample sensitivity results for the scenario A. Confidence bounds in colours and mean value plot in red Ek s p l o a t a c j a i NiEzawodNosc – Ma i N t E N a N c E a N d REliabi lit y Vo l . 24, No. 3, 2022 572 REFERENCES 1. Aguirre-Salado A I, Vaquera-Huerta H, Aguirre-Salado C A et al. Facing missing observations in data-A new approach for estimating strength of earthquakes on the Pacific coast of southern Mexico using random censoring. Applied Sciences (Switzerland) 2019; 9(14): 1–14, https://doi.org/10.3390/app9142863. 2. Ait Allal A, Mansouri K, El Had K, Youssfi M. World Class Ship Maintenance and Reliability for a Sustainable Shipping Industry Competitiveness. In Kacprzyk J, Balas VE, Ezziyyani M (eds): Advanced Intelligent Systems for Sustainable Development (AI2SD’2020), Cham, Springer International Publishing: 2022: 1–10. 3. Amruthnath N, Gupta T. Fault class prediction in unsupervised learning using model-based clustering approach. 2018 International Conference on Information and Computer Technologies, ICICT 2018., IEEE: 2018: 5–12, https://doi.org/10.1109/INFOCT.2018.8356831. 4. Andrzejczak K, Bukowski L. A method for estimating the probability distribution of the lifetime for new technical equipment based on expert judgement. Eksploatacja i Niezawodnosc – Maintenance and Reliability 2021; 23(4): 757–769, https://doi.org/10.17531/ein.2021.4.18. 5. Assis E M, Figueirôa Filho C L S, Lima G C et al. Comparison between maintenance policies based on q-Weibull and Weibull models. International Journal of Quality and Reliability Management 2022; 39(1): 258–279, https://doi.org/10.1108/IJQRM-09-2019-0283. 6. Aydi W, Alduais F S. Estimating weibull parameters using least squares and multilayer perceptron vs. bayes estimation. Computers, Materials and Continua 2022; 71(2): 4033–4050, https://doi.org/10.32604/cmc.2022.023119. 7. Bajic B, Rikalovic A, Suzic N, Piuri V. Industry 4.0 Implementation Challenges and Opportunities: A Managerial Perspective. IEEE Systems Journal 2021; 15(1): 546–559, https://doi.org/10.1109/JSYST.2020.3023041. 8. Chen Y, Liu Y, Jiang T. Optimal maintenance strategy for multi-state systems with single maintenance capacity and arbitrarily distributed maintenance time. Reliability Engineering and System Safety 2021; 211(February): 107576, https://doi.org/10.1016/j.ress.2021.107576. 9. Crespo Márquez A. Benefits of Digital Transformation for Maintenance Management Systems. Market Trends. In Crespo Márquez A (ed): Digital Maintenance Management. Guiding Digital Transformation in Maintenance, Cham, Springer International Publishing: 2022: 3–11, https://doi.org/10.1007/978-3-030-97660-6_1. 10. Crespo Márquez A. Dynamic modelling for supply chain management: Dealing with front-end, back-end and integration issues. Springer London: 2010. doi:10.1007/978-1-84882-681-6, https://doi.org/10.1007/978-1-84882-681-6. 11. Cui Z, Lee C, Liu Y. Single-transform formulas for pricing Asian options in a general approximation framework under Markov processes. European Journal of Operational Research 2018; 266(3): 1134–1139, https://doi.org/10.1016/j.ejor.2017.10.049. 12. Deviren Saygin S, Erpul G. Modeling aggregate size distribution of eroded sediment resulting from rain-splash and raindrop impacted flow processes. International Journal of Sediment Research 2019; 34(2): 166–177, https://doi.org/10.1016/j.ijsrc.2018.10.004. 13. Dhawalikar M N, Mariappan V, Srividhya P K, Kurtikar V. Multi-state failure phenomenon and analysis using semi-Markov model. International Journal of Quality and Reliability Management 2018; 35(9): 2080–2091, https://doi.org/10.1108/IJQRM-01-2016-0001. 14. Fernandes M, Corchado J M, Marreiros G. Machine learning techniques applied to mechanical fault diagnosis and fault prognosis in the context of real industrial manufacturing use-cases: a systematic literature review. Applied Intelligence 2022. doi:10.1007/s10489-022- 03344-3, https://doi.org/10.1007/s10489-022-03344-3. 15. Girtler J, Rudnicki J. The matter of decision-making control over operation processes of marine power plant systems with the use of their models in the form of semi-markov decision-making processes. Polish Maritime Research 2021; 28(1): 116–126, https://doi.org/10.2478/ pomr-2021-0011. 16. Hu J, Shen J, Shen L. Periodic preventive maintenance planning for systems working under a Markovian operating condition. Computers and Industrial Engineering 2020; 142: 106291, https://doi.org/10.1016/j.cie.2020.106291. 17. Hu J, Xu A, Li B, Liao H. Condition-based maintenance planning for multi-state systems under time-varying environmental conditions. Computers and Industrial Engineering 2021. doi:10.1016/j.cie.2021.107380, https://doi.org/10.1016/j.cie.2021.107380. 18. International Maritime Organization (IMO). Fourth IMO Greenhouse Gas Study 2020. London, 2021. 19. Jing H, Yao P, Song L et al. Analysis on optimization decision of pavement maintenance mode based on analytic hierarchy process. DYNA 2021; 96(6): 605–611, https://doi.org/10.6036/10332. 20. Khan M M, Tse P W, Trappey A J C. Development of a novel methodology for remaining useful life prediction of industrial slurry pumps in the absence of run to failure data. Sensors 2021; 21(24): 22, https://doi.org/10.3390/s21248420. 21. Kim D, Lee S, Kim D. An applicable predictive maintenance framework for the absence of run-to-failure data. Applied Sciences (Switzerland) 2021; 11(5180): 17, https://doi.org/10.3390/app11115180. 22. Kumar G, Varghese J P. Optimum preventive maintenance policy for a mechanical system using Semi-Markov method and Golden section technique. 2018 IEEE International Conference on Industrial Engineering and Engineering Management (IEEM), Bangkok, 2018: 232–236, https://doi.org/10.1109/IEEM.2018.8607343. 23. Li D, Zhou Y, Hu G, Spanos C J. Handling Incomplete Sensor Measurements in Fault Detection and Diagnosis for Building HVAC Systems. IEEE Transactions on Automation Science and Engineering 2020; 17(2): 833–846, https://doi.org/10.1109/TASE.2019.2948101. 24. Liu G, Chen S, Jin H, Liu S. Optimum opportunistic maintenance schedule incorporating delay time theory with imperfect maintenance. Reliability Engineering and System Safety 2021; 213: 107668, https://doi.org/10.1016/j.ress.2021.107668. 25. Lupton R C, Allwood J M. Incremental Material Flow Analysis with Bayesian Inference. Journal of Industrial Ecology 2018; 22(6): 1352– 1364, https://doi.org/10.1111/jiec.12698. 26. Ma Y F, Jia X, Hu Q et al. A New State Recognition and Prognosis Method Based on a Sparse Representation Feature and the Hidden Semi- Markov Model. IEEE Access 2020; 8: 119405–119420, https://doi.org/10.1109/ACCESS.2020.3005165. 27. Nobakhti A, Raissi S, Khalili K, Soltani R. Dynamic reliability assessment of a complex recovery system using fault tree, fuzzy inference and discrete event simulation. Eksploatacja i Niezawodnosc – Maintenance and Reliability 2021; 23(4): 593–604, https://doi.org/doi. org/10.17531/ein.2021.4.1. 28. Oh G. Fatigue fracture and probabilistic assessments of a cone and pipe welded structure of stainless steels. Fatigue and Fracture of Engineering Materials and Structures 2022; 45(1): 40–54, https://doi.org/10.1111/ffe.13581. 29. Parra C, Viveros P, Kristjanpoller F et al. Audit and diagnosis in asset management and maintenance applied in the electrical industry. DYNA 2021; 96(3): 238, https://doi.org/10.6036/10037. 30. Pech M, Vrchota J, Bednář J. Predictive maintenance and intelligent sensors in smart factory: Review. Sensors 2021; 21(4): 1–39, https://doi.