# Preventive Maintenance Interval Optimization for Continuous Multistate Systems

Authors: Ruiying Li; Xufeng Zhang
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

Preventive maintenance (PM), which is performed periodically on the system to lessen its failing probability, can effectively decrease the loss caused by the system breakdown or the performance degradation. The optimal PM interval has been well studied for both binary-state systems (BSSs) and discrete multistate systems (MSSs). However, in reality, the performance of many systems can change continuously, ranging from complete failure to perfect functioning. Considering such characteristics of systems, two types of performance-based measures, performance availability and probabilistic resilience, are addressed to quantify the system’s behaviour for continuous MSS. A Monte Carlo-based method is given to analyse the performance change process of the system, and an optimization framework is proposed to find the optimal PM interval with the considerations of per-unit-time cost, system breakdown rate, performance availability, and probabilistic resilience. A computer cluster is used as an example to illustrate the effectiveness of our proposed method.


---

## Body

## 1. Introduction

In industry, preventive maintenance (PM), a planned or scheduled maintenance activity, is an essential maintenance policy [1]. PM was proposed in 1940s, requiring engineers to take maintenance actions when the system is about to fail. In this way, the probability of system breakdown can be reduced or even avoided. Hence, how to determine the proper PM interval to avoid excessive or insufficient maintenance is a critical issue. With the practice of PM, maintenance activities changed from reactive to proactive. Later, condition-based maintenance (CBM) was proposed. Using such strategy, the optimal maintenance timing can be determined according to system’s current condition. However, some system’s performance cannot be easily monitored without interrupting the system’s operation, so PM is still widely applied in industry.The topic of how to determine an optimal PM interval for systems attracts a lot of researchers. The excessive maintenance may affect the normal operation of the system and result in high maintenance costs. On the contrary, the insufficient maintenance may cause frequent system breakdown, affecting the user’s experience. The optimal PM planning problem was firstly studied for binary-state systems (BSSs), which simplifies the system states into normal and failure [2, 3]. Group/block replacement models [4, 5] and opportunistic maintenance models [6] were proposed successively to find the optimal PM planning. However, if a system with multiple performance states is assumed as a BSS, its reliability is computed as the probability that their performance stayed above the threshold, and its performance changing process is then neglected. Considering that the system may have multiple operable states or operate at multiple performance levels, reliability theory then started to focus on discrete multistate systems (MSSs) [7], which has more than two but finite performance states. Levitin and Lisnianski [8] firstly studied the maintenance decision problem of discrete MSS and proposed the imperfect PM optimization framework in [9]. Then, many researchers proposed or applied several methods, such as the Markov process [10] and the universal generating function (UGF) [8, 9, 11], to find the optimal PM interval for discrete MSS. For example, Liu et al. [12] proposed a value-centric dynamic PM policy, and they determined which component to maintain and the maintenance level by maximizing the maintenance net value. However, in reality, the performance of many systems can change continuously, ranging from complete failure to perfect functioning. These systems are named as continuous MSS, or continuous-state systems briefly [13]. Obviously, binary-state-based or discrete multistate-based PM interval optimization frameworks are insufficient to describe the characteristics of the continuous MSS. For example, the Markov method is based on the models of system states and the state transitions, and UGF method is based on the possible system states and its corresponding probability. The infinite number of performance states makes the Markov and UGF methods difficult to apply for continuous MSS.The continuous changing performance of such systems requires advanced measures. Performance-based measures, e.g., performance availability and resilience, provide a new perspective to analyse such systems. Performance availability, an extension of traditional availability in the performance dimension, qualifies the average performance of the continuous MSS over its lifetime. Probabilistic resilience, a new measure to assess the system’s ability to withstand disturbances and return to a normal state quickly [14], focuses on the system performance after the disturbance occurs. Using such performance-based measures, Zhang and Li [15] compared different PM intervals of a discrete MSS. However, how to find the optimal PM interval and how to apply such method to continuous MSS was not addressed.To solve such problems, we propose a PM interval optimization method for continuous MSS. The contributions of this paper include the following:(1)
Two performance-based measures, performance availability and probabilistic resilience, are proposed to analyse and optimize the PM interval of continuous MSS with the consideration of its infinite number of performance states. Using such measures, the continuous MSS can be analysed more accurate, as its continuous changing performance is depicted.(2)
A general PM interval optimization framework is provided for continuous MSS. In the optimization model, we consider four types of decision elements, including the performance availability, the probabilistic resilience, the system breakdown rate, and the per-unit-time cost. Compared with the PM interval optimization methods for BSS or discrete MSS, the system’s continuous performance changing process is analysed by a Monte Carlo-based method in our framework.The remainder of the paper is organized as follows: Section2 describes our PM decision problem in detail. In Section 3 four key parameters, performance availability, probabilistic resilience, system breakdown rate, and per-unit-time cost, are proposed and applied as the decision elements in the PM interval optimization models for the continuous MSS. Monte Carlo method is used to compute the performance degradation and recovery process of the system and its components. In Section 4, the optimal PM interval is calculated for a computer cluster, and the optimal result is compared with the one obtained by enumeration to verify the effectiveness of our method. Finally, concluding remarks are provided in Section 5.


## 2. Problem Description

In this paper, we consider a continuous MSS which is composed of continuous multistate components. When a disturbance occurs on a component, its performance changes as Figure1 shows. There are two phases:(1)
Degradation phase: after the disturbance occurs on the component, its performance begins to decline. The component continuously operates in a degraded state until it breaks down or the PM activity starts.(2)
Recovery phase: the performance of the component begins to restore back after the corrective maintenance (CM) or preventive maintenance (PM) activity is taken.Figure 1
The resilience process of a continuous component/system.By performance monitoring and data fitting, functions of both the degradation and recovery processes of the components can be obtained. Cimellaro et al. [16] proposed three typical functions to describe the recovery behaviour of the MSS, including linear function, exponential function, and trigonometric function. These functions can also be used to describe the performance degradation process, see the study of Crk [17]. So equations (1)–(3) are applied in our paper to describe both the performance degradation and recovery processes:(1)ft=at−tsTp+b,(2)ft=a×exp−b×t−tsTp,(3)ft=a21+cosπbt−tsTp,where a and b are the constant values that can be calculated by curve fitting, ts is the start time of the corresponding process (i.e., the time that the component’s performance begins to decline or restore in the degradation and recovery process, respectively), and Tp is the time duration of the degradation or recovery process.As is well known, the performance of the system depends on that of its components. When the performance of components starts to degrade, that of the system may also decline. Then, the system will operate in a degradation state. When the system’s performance drops below the thresholdQL, CM activities will be taken on the components, and the system will restore back along with the components. Note that, if the system’s PM interval arrives before the system breaks down, PM activities are taken instead.In this paper, we aim to find an optimal PM interval for the continuous MMS composed of continuous multistate components. Our assumptions are as follows:(1)
The system after PM or CM is “as good as new” (i.e., perfect maintenance)(2)
There is only one maintenance channel, and only one component can be repaired each time(3)
The maintenance activity is performed immediately without administrative delay(4)
The disturbance occurrence time of components are independently and identically distributed(5)
Both the time duration of the performance degradation and recovery processes are also independently and identically distributed


## 3. Methods

In this section, we propose a PM interval optimization framework for continuous MSS based on the Monte Carlo method.


### 3.1. Decision Elements

Generally, the optimal PM interval is the time when the system is about to fail. As the performance changing is very essential for continuous MSS, we propose two types of performance-based measures, performance availability and probabilistic resilience, and apply them together with traditional PM decision elements, system breakdown rate, and per-unit-time cost, to find the optimal PM interval.


#### 3.1.1. Performance Availability

Availability is the proportion of time that a system is in a functioning condition for binary-state systems, and it measures the degree to which the system can be used. For MSS, the performance behaviour of the system also affects its availability. In Figure1, the solid line and the dotted line describe the real performance and the ideal performance of the system, respectively, and the performance availability can then be defined.Definition 1.
Performance availability is the ratio of the area under the system’s real performance curve to that under the ideal performance curve over time 0 toT.
It can be computed as(4)AP=∫0TQ1tdt∫0TQ0tdt,where Q0 and Q1 are the ideal and real performance functions of the system, respectively, and T is the operation cycle (e.g., 1 year). The performance availability is the average proportion of performance that a system behaves in [0, T]. One can see that the performance availability is an extension of the traditional availability in the performance dimension. For BSS, we have Q0 = 1, Q1 = 1 (if the system is operational or up), and Q1 = 0 (if the system fails or down). So our performance availability can be calculated as follows:(5)AP=uptimeoperation cycle,for BSS.


#### 3.1.2. Probabilistic Resilience

Resilience is also a performance-based measure, and it reflects the ability of the system to withstand disturbance and return to a normal state quickly. Under a given disturbance, a deterministic resilience measure can be defined.Definition 2.
Deterministic resilience is the ratio of the area beneath the performance curve within its maximum allowable recovery time after a given disturbance to the ideal one in the normal state.
This deterministic measure can be calculated as(6)ℝD=∫t0t0+TaQ1tdt∫t0t0+TaQ0tdt,where t0 is the disturbance occurrence time and Ta is system’s maximum allowable recovery time. ℝD presents the average performance of the system in t0,t0+Ta after the disturbance. Such measure is with clear physical meaning which reflects the “bounce back” ability of the system and can be used to compare the resilience of various systems on the same relative scale. Obviously, the system’s resilience is a random variable, and we define probabilistic resilience under a random disturbance.Definition 3.
Probabilistic resilience is the probability that the system’s resilience can satisfy users’ requirements under random disturbances.
The probabilistic resilience can be computed as(7)ℝP=PrℝD≥ℝD∗,where ℝD∗ is the resilience threshold defined by users. Such measure describes how well the system can meet its resilience requirements.


#### 3.1.3. System Breakdown Rate

System breakdown rate implies the probability that the system breaks down before the PM interval arrives. It can be defined as(8)η=PrSystem breaks down before PM interval arrives.


#### 3.1.4. Per-Unit-Time Cost

Cost is usually used in PM interval optimization to analyse the expense of the maintenance strategy. Considering the operation and maintenance process of the system, the cost includes six aspects: (1) operation cost; (2) maintenance labor cost; (3) maintenance material cost; (4) on-site maintenance cost; (5) system breakdown loss; and (6) system performance loss. Here, the operation cost reflects the resource consumption during the system operation; maintenance labor and material cost includes the labor and spare consumption and maintenance tool, equipment, and facility depreciation during the maintenance process; on-site maintenance cost refers to the cost of the preparation activity of each maintenance; system breakdown loss contains both the direct and indirect losses caused by the system breakdown; and system performance loss reflects the loss caused by the system performance degradation. Thus, we define the per-unit-time cost function as(9)CPUT=T∑i∫allsiCO,isifisidsi+CMMHC×tPM+tCM+∑iCMMC,iNi+CBMNM+CD×tD+CP∫0TQ0t−Q1tdtT,where T is the operation cycle (e.g., 1 year); CO,isi is the operation cost of component i under state si per unit time, fisi is the probability density function of state si; CMMHC is the maintenance labor cost per unit time, tPM and tCM are the total time of PM and CM over time [0, T], respectively; CMMC,i is the material cost for component i, and Ni is the number of times that the performance of component i degrades and needs to be repaired over time [0, T]; CBM is the on-site cost for one maintenance, and NM is the number of on-site maintenance over time [0, T]; CD is the system breakdown loss per unit time, and tD is the length of the system downtime over time [0, T]; and CP is the system performance loss per unit area in Figure 1.


### 3.2. Monte Carlo Simulation

For BSS and discrete MSS, Universal Generation Function and Markov process are often used in reliability analysis. However, for continuous MSS, it is impossible to model all the states and solve the corresponding differential equations. So the Monte Carlo method is used here to analyse the performance changing process of the continuous MSS under given PM intervals. As the system after PM or CM is supposed “as good as new,” we regard one simulation cycle as the time period from the time that the system is put into use to the time that its performance fully restores back after the maintenance activity (i.e, PM or CM) is taken. In Figure1, this time period is denoted as time [0, t2]. Multiple simulation cycles together form the operation cycle T. The randomness of both disturbances that the system suffered and the system’s response are then considered in each simulation cycle. In each simulation cycle, the simulation process is as follows:(1)
Sample and obtain the disturbance occurrence time and time duration of the performance degradation for each component according to their corresponding distributions.(2)
Sort the disturbance occurrence time in the ascending order, and denote them ast1, t2, t3....(3)
Letα = 0 and β = 0, where α is the number of the disturbances that have occurred and β is the number of disturbances in the current continuous disturbance behaviour. A continuous disturbance behaviour describes the situation that a new disturbance occurs on the system before the performance degradation of the previous one ends, and then these disturbances affect the system together. These continuous disturbances are denoted as a disturbance group.(4)
Letα = α + 1 and β = β + 1. Calculate the performance degradation end time of the current disturbance group tde by comparing it with the occurrence time of the α+1th disturbance.(a)
Iftde>tα+1, continue Step (4)(b)
Iftde>tα+1, the αth disturbance is the last one in the current disturbance group, go to Step (5)(5)
Calculate the performance degradation process of the current disturbance group, and determine whether the performance degradation in the current disturbance group will cause the system breakdown (i.e., the system performance drops below the thresholdQL). If so, record the time that the system breaks down as tQ.(6)
Comparetde and tQ with the PM interval TPM:(a)
If min{tde, tQ, TPM} = tQ, it means that the system breaks down in the current continuous disturbance process before the PM interval arrives. Go to Step (7) and start CM.(b)
If min{tde, tQ, TPM} = TPM, it means that the PM interval arrives in this continuous disturbance process before the system breaks downs. Go to Step (7) and start PM.(c)
If min{tde, tQ, TPM} = tde, it means that neither system breaks down nor PM interval arrives in this continuous disturbance process. Let β = 0, and go to Step (4).(7)
Sample and obtain the performance recovery time for each degraded component according to their distributions, and compute the recovery process of the system.From the Monte Carlo simulation, the system’s performance can be recorded every Δt time, and the decision elements can be estimated. For example, the performance availability can be estimated as(10)A^P=∑j=1N∑k=0te,j/Δt−1Q1,j,k+Q1,j,k+1∑j=1N∑k=0te,j/Δt−1Q0,j,k+Q0,j,k+1,where N is the number of iterations in the simulation, te,j is the system recovery end time in the jth simulation iteration, and Q0,j,k and Q1,j,k are the ideal and real system performance at time k·Δt in the jth simulation iteration, respectively. The performance-based probabilistic resilience can be estimated as(11)ℝ^P=NℝN,where N is the number of simulation iterations, Nℝ is the number of iterations whose deterministic resilience ℝ^D satisfies the resilience threshold ℝD∗ defined by users, and(12)ℝ^D,j=∑k=ts,j/Δtte,j/Δt−1Q1,j,k+Q1,j,k+1∑k=ts,j/Δtte,j/Δt−1Q0,j,k+Q0,j,k+1,where ts,j is the occurrence time of the first disturbance in the jth simulation iteration. The system breakdown rate can be estimated as(13)η^=NdN,where Nd is the number of the system breaks down (i.e., the number of CM) in the simulation. The per-unit-time cost can be estimated as(14)C^PUT=∑i=1nCO,isiΔtisi+CMMHC×∑j=1NtPM,j+tCM,j+∑i=1nCMMC,iNi+CBMNPM+NCM+CD×∑j=1NtD,j+CP∑j=1N∑k=ts,j/Δtte,j/Δt−1Q0,j,k−Q1,j,kΔt∑j=1Nte,j,where Δtisi is the time length of component i under state si in the whole simulation; tPM,j and tCM,j are the maintenance time of PM and CM in the jth iteration, respectively; n is the number of components in the system; NPM and NCM are the number of PM or CM in the whole simulation; and tD,j is the system downtime in the jth iteration.AfterN rounds of iterations, we can obtain N performance change processes of the system, and then the statistical values of the four decision elements under the given PM interval can be further calculated. The value of N depends on the requirements of simulation accuracy ε. According to the central limit theorem, the number of iterations can be computed as(15)N≥maxxZα/2σxεx2,where Zα/2 is the 1001−α/2th percentile of the standard normal distribution, 1−α is the confidence level (e.g., 1−α = 95%), and εx and σx are the accuracy requirement and mean square error for the parameter x (i.e., the four decision elements).


### 3.3. Optimization Model and Algorithm

Using the four decision elements mentioned above, we can build a PM interval optimization model for the continuous MSS. Depending on the problem, the previous PM interval optimization models can be divided into integer linear programming model [18], dynamic programming model [19], robust optimization model [20], and so on. For the continuous MSS, goal(s) and constraints in the optimization model can be determined according to the real problem. Hence, we provide two types of optimization models as examples. The first optimization model is a single-objective problem that minimizing the per-unit-time cost under the other three constraints as follows:(16)minCPUTs.t.η≤η∗AP≥AP∗ℝP≥ℝP∗,where AP∗, ℝP∗, and η∗ are the threshold of the performance availability, probabilistic resilience, and system breakdown rate, respectively. The second optimization model is a multiobjective optimization problem which minimizes the per-unit-time cost and the system breakdown rate while maximizes the performance availability and probabilistic resilience:(17)minCPUT,η,−AP,−ℝPTs.t.CPUT≤CPUT∗η≤η∗AP≥AP∗ℝP≥ℝP∗,where CPUT∗ is the threshold of the per-unit-time cost. Using the direct weighting method, this multiobjective optimization problem in equation (17) can be converted into a single one as(18)minf=k1CPUT+k2η−k3AP−k4ℝPs.t.CPUT≤CPUT∗η≤η∗AP≥AP∗ℝP≥ℝP∗,where f is the joint objective function and ki is the weighting factor for each objective. ki can be computed as ki=ki1×ki2, where ki1 reflects the importance of each objective, ∑ki1=1, and ki2 is used to adjust the effect of the difference in magnitude among objectives (e.g., ki2=1/∇fi2).To solve the PM interval optimization problem, three algorithms were usually applied, including genetic algorithms (GAs) [8, 9, 21], particle swarm optimization (PSO) [22], and simulated annealing (SA) [23]. Generally, the appropriate optimization algorithm is selected according to the type of the optimization model. For example, the simplex method can be used to solve linear optimization problems, golden section method and binary search method can be applied for nonlinear optimization problems with unimodal objective functions, and genetic algorithms, particle swarm optimization, and simulated annealing (SA) are applicable for nonlinear optimization problems with multimodal objective functions. In practice, the optimization algorithm should be determined according to the specific problem.


## 3.1. Decision Elements

Generally, the optimal PM interval is the time when the system is about to fail. As the performance changing is very essential for continuous MSS, we propose two types of performance-based measures, performance availability and probabilistic resilience, and apply them together with traditional PM decision elements, system breakdown rate, and per-unit-time cost, to find the optimal PM interval.


### 3.1.1. Performance Availability

Availability is the proportion of time that a system is in a functioning condition for binary-state systems, and it measures the degree to which the system can be used. For MSS, the performance behaviour of the system also affects its availability. In Figure1, the solid line and the dotted line describe the real performance and the ideal performance of the system, respectively, and the performance availability can then be defined.Definition 1.
Performance availability is the ratio of the area under the system’s real performance curve to that under the ideal performance curve over time 0 toT.
It can be computed as(4)AP=∫0TQ1tdt∫0TQ0tdt,where Q0 and Q1 are the ideal and real performance functions of the system, respectively, and T is the operation cycle (e.g., 1 year). The performance availability is the average proportion of performance that a system behaves in [0, T]. One can see that the performance availability is an extension of the traditional availability in the performance dimension. For BSS, we have Q0 = 1, Q1 = 1 (if the system is operational or up), and Q1 = 0 (if the system fails or down). So our performance availability can be calculated as follows:(5)AP=uptimeoperation cycle,for BSS.


### 3.1.2. Probabilistic Resilience

Resilience is also a performance-based measure, and it reflects the ability of the system to withstand disturbance and return to a normal state quickly. Under a given disturbance, a deterministic resilience measure can be defined.Definition 2.
Deterministic resilience is the ratio of the area beneath the performance curve within its maximum allowable recovery time after a given disturbance to the ideal one in the normal state.
This deterministic measure can be calculated as(6)ℝD=∫t0t0+TaQ1tdt∫t0t0+TaQ0tdt,where t0 is the disturbance occurrence time and Ta is system’s maximum allowable recovery time. ℝD presents the average performance of the system in t0,t0+Ta after the disturbance. Such measure is with clear physical meaning which reflects the “bounce back” ability of the system and can be used to compare the resilience of various systems on the same relative scale. Obviously, the system’s resilience is a random variable, and we define probabilistic resilience under a random disturbance.Definition 3.
Probabilistic resilience is the probability that the system’s resilience can satisfy users’ requirements under random disturbances.
The probabilistic resilience can be computed as(7)ℝP=PrℝD≥ℝD∗,where ℝD∗ is the resilience threshold defined by users. Such measure describes how well the system can meet its resilience requirements.


### 3.1.3. System Breakdown Rate

System breakdown rate implies the probability that the system breaks down before the PM interval arrives. It can be defined as(8)η=PrSystem breaks down before PM interval arrives.


### 3.1.4. Per-Unit-Time Cost

Cost is usually used in PM interval optimization to analyse the expense of the maintenance strategy. Considering the operation and maintenance process of the system, the cost includes six aspects: (1) operation cost; (2) maintenance labor cost; (3) maintenance material cost; (4) on-site maintenance cost; (5) system breakdown loss; and (6) system performance loss. Here, the operation cost reflects the resource consumption during the system operation; maintenance labor and material cost includes the labor and spare consumption and maintenance tool, equipment, and facility depreciation during the maintenance process; on-site maintenance cost refers to the cost of the preparation activity of each maintenance; system breakdown loss contains both the direct and indirect losses caused by the system breakdown; and system performance loss reflects the loss caused by the system performance degradation. Thus, we define the per-unit-time cost function as(9)CPUT=T∑i∫allsiCO,isifisidsi+CMMHC×tPM+tCM+∑iCMMC,iNi+CBMNM+CD×tD+CP∫0TQ0t−Q1tdtT,where T is the operation cycle (e.g., 1 year); CO,isi is the operation cost of component i under state si per unit time, fisi is the probability density function of state si; CMMHC is the maintenance labor cost per unit time, tPM and tCM are the total time of PM and CM over time [0, T], respectively; CMMC,i is the material cost for component i, and Ni is the number of times that the performance of component i degrades and needs to be repaired over time [0, T]; CBM is the on-site cost for one maintenance, and NM is the number of on-site maintenance over time [0, T]; CD is the system breakdown loss per unit time, and tD is the length of the system downtime over time [0, T]; and CP is the system performance loss per unit area in Figure 1.


## 3.1.1. Performance Availability

Availability is the proportion of time that a system is in a functioning condition for binary-state systems, and it measures the degree to which the system can be used. For MSS, the performance behaviour of the system also affects its availability. In Figure1, the solid line and the dotted line describe the real performance and the ideal performance of the system, respectively, and the performance availability can then be defined.Definition 1.
Performance availability is the ratio of the area under the system’s real performance curve to that under the ideal performance curve over time 0 toT.
It can be computed as(4)AP=∫0TQ1tdt∫0TQ0tdt,where Q0 and Q1 are the ideal and real performance functions of the system, respectively, and T is the operation cycle (e.g., 1 year). The performance availability is the average proportion of performance that a system behaves in [0, T]. One can see that the performance availability is an extension of the traditional availability in the performance dimension. For BSS, we have Q0 = 1, Q1 = 1 (if the system is operational or up), and Q1 = 0 (if the system fails or down). So our performance availability can be calculated as follows:(5)AP=uptimeoperation cycle,for BSS.


## 3.1.2. Probabilistic Resilience

Resilience is also a performance-based measure, and it reflects the ability of the system to withstand disturbance and return to a normal state quickly. Under a given disturbance, a deterministic resilience measure can be defined.Definition 2.
Deterministic resilience is the ratio of the area beneath the performance curve within its maximum allowable recovery time after a given disturbance to the ideal one in the normal state.
This deterministic measure can be calculated as(6)ℝD=∫t0t0+TaQ1tdt∫t0t0+TaQ0tdt,where t0 is the disturbance occurrence time and Ta is system’s maximum allowable recovery time. ℝD presents the average performance of the system in t0,t0+Ta after the disturbance. Such measure is with clear physical meaning which reflects the “bounce back” ability of the system and can be used to compare the resilience of various systems on the same relative scale. Obviously, the system’s resilience is a random variable, and we define probabilistic resilience under a random disturbance.Definition 3.
Probabilistic resilience is the probability that the system’s resilience can satisfy users’ requirements under random disturbances.
The probabilistic resilience can be computed as(7)ℝP=PrℝD≥ℝD∗,where ℝD∗ is the resilience threshold defined by users. Such measure describes how well the system can meet its resilience requirements.


## 3.1.3. System Breakdown Rate

System breakdown rate implies the probability that the system breaks down before the PM interval arrives. It can be defined as(8)η=PrSystem breaks down before PM interval arrives.


## 3.1.4. Per-Unit-Time Cost

Cost is usually used in PM interval optimization to analyse the expense of the maintenance strategy. Considering the operation and maintenance process of the system, the cost includes six aspects: (1) operation cost; (2) maintenance labor cost; (3) maintenance material cost; (4) on-site maintenance cost; (5) system breakdown loss; and (6) system performance loss. Here, the operation cost reflects the resource consumption during the system operation; maintenance labor and material cost includes the labor and spare consumption and maintenance tool, equipment, and facility depreciation during the maintenance process; on-site maintenance cost refers to the cost of the preparation activity of each maintenance; system breakdown loss contains both the direct and indirect losses caused by the system breakdown; and system performance loss reflects the loss caused by the system performance degradation. Thus, we define the per-unit-time cost function as(9)CPUT=T∑i∫allsiCO,isifisidsi+CMMHC×tPM+tCM+∑iCMMC,iNi+CBMNM+CD×tD+CP∫0TQ0t−Q1tdtT,where T is the operation cycle (e.g., 1 year); CO,isi is the operation cost of component i under state si per unit time, fisi is the probability density function of state si; CMMHC is the maintenance labor cost per unit time, tPM and tCM are the total time of PM and CM over time [0, T], respectively; CMMC,i is the material cost for component i, and Ni is the number of times that the performance of component i degrades and needs to be repaired over time [0, T]; CBM is the on-site cost for one maintenance, and NM is the number of on-site maintenance over time [0, T]; CD is the system breakdown loss per unit time, and tD is the length of the system downtime over time [0, T]; and CP is the system performance loss per unit area in Figure 1.


## 3.2. Monte Carlo Simulation

For BSS and discrete MSS, Universal Generation Function and Markov process are often used in reliability analysis. However, for continuous MSS, it is impossible to model all the states and solve the corresponding differential equations. So the Monte Carlo method is used here to analyse the performance changing process of the continuous MSS under given PM intervals. As the system after PM or CM is supposed “as good as new,” we regard one simulation cycle as the time period from the time that the system is put into use to the time that its performance fully restores back after the maintenance activity (i.e, PM or CM) is taken. In Figure1, this time period is denoted as time [0, t2]. Multiple simulation cycles together form the operation cycle T. The randomness of both disturbances that the system suffered and the system’s response are then considered in each simulation cycle. In each simulation cycle, the simulation process is as follows:(1)
Sample and obtain the disturbance occurrence time and time duration of the performance degradation for each component according to their corresponding distributions.(2)
Sort the disturbance occurrence time in the ascending order, and denote them ast1, t2, t3....(3)
Letα = 0 and β = 0, where α is the number of the disturbances that have occurred and β is the number of disturbances in the current continuous disturbance behaviour. A continuous disturbance behaviour describes the situation that a new disturbance occurs on the system before the performance degradation of the previous one ends, and then these disturbances affect the system together. These continuous disturbances are denoted as a disturbance group.(4)
Letα = α + 1 and β = β + 1. Calculate the performance degradation end time of the current disturbance group tde by comparing it with the occurrence time of the α+1th disturbance.(a)
Iftde>tα+1, continue Step (4)(b)
Iftde>tα+1, the αth disturbance is the last one in the current disturbance group, go to Step (5)(5)
Calculate the performance degradation process of the current disturbance group, and determine whether the performance degradation in the current disturbance group will cause the system breakdown (i.e., the system performance drops below the thresholdQL). If so, record the time that the system breaks down as tQ.(6)
Comparetde and tQ with the PM interval TPM:(a)
If min{tde, tQ, TPM} = tQ, it means that the system breaks down in the current continuous disturbance process before the PM interval arrives. Go to Step (7) and start CM.(b)
If min{tde, tQ, TPM} = TPM, it means that the PM interval arrives in this continuous disturbance process before the system breaks downs. Go to Step (7) and start PM.(c)
If min{tde, tQ, TPM} = tde, it means that neither system breaks down nor PM interval arrives in this continuous disturbance process. Let β = 0, and go to Step (4).(7)
Sample and obtain the performance recovery time for each degraded component according to their distributions, and compute the recovery process of the system.From the Monte Carlo simulation, the system’s performance can be recorded every Δt time, and the decision elements can be estimated. For example, the performance availability can be estimated as(10)A^P=∑j=1N∑k=0te,j/Δt−1Q1,j,k+Q1,j,k+1∑j=1N∑k=0te,j/Δt−1Q0,j,k+Q0,j,k+1,where N is the number of iterations in the simulation, te,j is the system recovery end time in the jth simulation iteration, and Q0,j,k and Q1,j,k are the ideal and real system performance at time k·Δt in the jth simulation iteration, respectively. The performance-based probabilistic resilience can be estimated as(11)ℝ^P=NℝN,where N is the number of simulation iterations, Nℝ is the number of iterations whose deterministic resilience ℝ^D satisfies the resilience threshold ℝD∗ defined by users, and(12)ℝ^D,j=∑k=ts,j/Δtte,j/Δt−1Q1,j,k+Q1,j,k+1∑k=ts,j/Δtte,j/Δt−1Q0,j,k+Q0,j,k+1,where ts,j is the occurrence time of the first disturbance in the jth simulation iteration. The system breakdown rate can be estimated as(13)η^=NdN,where Nd is the number of the system breaks down (i.e., the number of CM) in the simulation. The per-unit-time cost can be estimated as(14)C^PUT=∑i=1nCO,isiΔtisi+CMMHC×∑j=1NtPM,j+tCM,j+∑i=1nCMMC,iNi+CBMNPM+NCM+CD×∑j=1NtD,j+CP∑j=1N∑k=ts,j/Δtte,j/Δt−1Q0,j,k−Q1,j,kΔt∑j=1Nte,j,where Δtisi is the time length of component i under state si in the whole simulation; tPM,j and tCM,j are the maintenance time of PM and CM in the jth iteration, respectively; n is the number of components in the system; NPM and NCM are the number of PM or CM in the whole simulation; and tD,j is the system downtime in the jth iteration.AfterN rounds of iterations, we can obtain N performance change processes of the system, and then the statistical values of the four decision elements under the given PM interval can be further calculated. The value of N depends on the requirements of simulation accuracy ε. According to the central limit theorem, the number of iterations can be computed as(15)N≥maxxZα/2σxεx2,where Zα/2 is the 1001−α/2th percentile of the standard normal distribution, 1−α is the confidence level (e.g., 1−α = 95%), and εx and σx are the accuracy requirement and mean square error for the parameter x (i.e., the four decision elements).


## 3.3. Optimization Model and Algorithm

Using the four decision elements mentioned above, we can build a PM interval optimization model for the continuous MSS. Depending on the problem, the previous PM interval optimization models can be divided into integer linear programming model [18], dynamic programming model [19], robust optimization model [20], and so on. For the continuous MSS, goal(s) and constraints in the optimization model can be determined according to the real problem. Hence, we provide two types of optimization models as examples. The first optimization model is a single-objective problem that minimizing the per-unit-time cost under the other three constraints as follows:(16)minCPUTs.t.η≤η∗AP≥AP∗ℝP≥ℝP∗,where AP∗, ℝP∗, and η∗ are the threshold of the performance availability, probabilistic resilience, and system breakdown rate, respectively. The second optimization model is a multiobjective optimization problem which minimizes the per-unit-time cost and the system breakdown rate while maximizes the performance availability and probabilistic resilience:(17)minCPUT,η,−AP,−ℝPTs.t.CPUT≤CPUT∗η≤η∗AP≥AP∗ℝP≥ℝP∗,where CPUT∗ is the threshold of the per-unit-time cost. Using the direct weighting method, this multiobjective optimization problem in equation (17) can be converted into a single one as(18)minf=k1CPUT+k2η−k3AP−k4ℝPs.t.CPUT≤CPUT∗η≤η∗AP≥AP∗ℝP≥ℝP∗,where f is the joint objective function and ki is the weighting factor for each objective. ki can be computed as ki=ki1×ki2, where ki1 reflects the importance of each objective, ∑ki1=1, and ki2 is used to adjust the effect of the difference in magnitude among objectives (e.g., ki2=1/∇fi2).To solve the PM interval optimization problem, three algorithms were usually applied, including genetic algorithms (GAs) [8, 9, 21], particle swarm optimization (PSO) [22], and simulated annealing (SA) [23]. Generally, the appropriate optimization algorithm is selected according to the type of the optimization model. For example, the simplex method can be used to solve linear optimization problems, golden section method and binary search method can be applied for nonlinear optimization problems with unimodal objective functions, and genetic algorithms, particle swarm optimization, and simulated annealing (SA) are applicable for nonlinear optimization problems with multimodal objective functions. In practice, the optimization algorithm should be determined according to the specific problem.


## 4. A Numerical Example

A computer cluster contains 10 connected computers that work together, as shown in Figure2. Both the system and its components are with continuous multistates, and the load capacity of system (Q) is the sum of that of computers (qi). Under the normal state, the load capacity of the system and components are both normalized as 100%. The components may suffer various disturbances, and the occurrence probability of these disturbances follows the exponential distribution with parameter λ. Once a disturbance occurs on a computer, its load capacity degrades gradually following the exponential function in equation (2), and the time duration of the performance degradation follows exponential distribution with parameter τ. When the computer’s performance degrades below the threshold (qL = 50%), it stops working and waits for maintenance. The system’s load capacity degrades along with that of the computers, and it can work under degradation state until the system breaks down (i.e., the load capacity of the system drops below the system’s threshold (QL = 60%)) or the PM interval (T) arrivals. If the system breaks down before the PM interval arrives, CM starts. During CM, computers that suffered disturbance are repaired offline one after another. When the number of perfect computers reaches the threshold that can provide enough performance, the system restarts working. Then, the CM activity continues until all computers are repaired and the whole system returns back to normal. If the PM interval arrives before the system breaks down, PM starts and the degraded computers are repaired online one by one. The recovery process of computers also follows the exponential function in equation (2), and the maintenance time for PM and CM follows normal distributions with parameters (μPM, σPM2) and (μCM, σCM2), respectively. It is assumed that the maintenance strategy is first-in-first-out (i.e., first-degrade-first-repair), and there is only one maintenance channel. Our aim is to find the optimal PM interval for the continuous multistate computer cluster. The resilience-related and cost-related parameters are shown in Tables 1 and 2.Figure 2
A computer cluster.Table 1
Resilience-related parameters.
Parameter
Distribution
Parameter 1
Value
Parameter 2
Value
Disturbance/failure occurrence time
Exponential
λ
1/8760 h−1
—
—
Degradation time
Exponential
τ
1/100 h−1
—
—
PM time
Normal
μ
PM
24
σ
PM
2
CM time
Normal
μ
CM
24
σ
CM
2Table 2
Cost-related parameters.
Parameter
Value
Operation cost per unit time of a single computer (CO)
5 yuan/hour
Maintenance labor cost per unit time (CMMHC)
100 yuan/hour
Maintenance material cost of a single computer (CMMC)
3000 yuan
On-site cost per maintenance (CBM)
1000 yuan
System breakdown loss per unit time (CD)
10000 yuan/hour
System performance loss per unit area (CP)
100 yuan/area


### 4.1. Resilience Behaviour Analysis

We run the Monte Carlo simulation and obtained the resilience process of the computer cluster. During the operation, system may suffer both single and continuous disturbances. The difference is whether the performance degradation of the system under the disturbance ends before a new one occurs. In Figure3, the first disturbance is a single one, and the other three are continuous ones which forms a disturbance group. In this figure, the four red stars show the time that the disturbance occurred. One can see that the performance degradation caused by the last three disturbances is superimposed on each other.Figure 3
Single and continuous disturbance effect.Figure4 shows two types of typical resilience behaviours for the computer cluster, and the PM interval was set as 900 hours and 2000 hours, respectively. In Figure 4(a), four disturbances occurred in sequence, resulting in the gradual load capacity degradation of four computers in the system. The performance of the system degraded along with that of the computers, as shown in this figure. At time t = 900 hours, the PM interval arrived while the system was still working, and PM was then taken to restore the system. The four degraded computers were repaired one by one according to the disturbance occurrence sequence. As the PM activity was online, the system’s performance restored back gradually in this figure. In Figure 4(b), the PM interval was set as 2000 hours, and the load capacity of the system dropped below its threshold QL at time t = 1212 hours. Then, the offline CM was taken, and the four computers were completely repaired before they restarted working. So one can see that the system performance recovery process in Figure 4(b) presents a step-like performance recovery process.Figure 4
Typical resilience processes of the cluster system. (a) Recover by PM (PM interval = 900 hours). (b) Recover by CM (PM interval = 2000 hours).
(a)
(b)


### 4.2. PM Interval Optimization

Firstly, a multiobjective optimization model was built to find the optimal PM interval as follows:(19)minf=k1CPUT+k2η−k3AP−k4ℝPs.t.AP≥0.9ℝPℝD≥0.8≥0.9η≤0.1,where ki is the weighting factor for each decision element. According to the direct weighing method, ki=ki1×ki2. Here, K1 = [0.6, 0.2, 0.1, 0.1]T described the importance of the four elements, and K2 was determined according to the magnitude difference among these four elements as K2 = [0.0018, 66.4896, 35.28, 10.2043]T. Compared with the optimization model in equation (17), the constraint of the per-unit-time cost is omitted, as there is no budget limit in this problem.Secondly, we used the Monte Carlo method to analyse the behaviour of the cluster system under different PM intervals. We run the simulation for 1000 times at the PM interval 1000 hours and obtained the corresponding sample standard deviations for the four decision elements. LetεC^PUT=5, εη^=0.002, εA^P=0.002, εℝ^P=0.002, and 1−α = 95%, we finally determined the number of iterations as N = 2700 according to equation (15). Under such 2700 simulation iterations, the error can be computed as εCPUT′<4.3627, εAP′<0.0019, εη′<0.0018, and εℝ^P′<0.0017.Finally, optimization algorithms, the golden section method and genetic algorithm, were used to solve the optimization problem in equation (19). The specific algorithms are as follows:(1)
Golden-section search is a technique for finding an extremum of a unimodal function inside a specified interval [a, b]. The functional values of fx at two golden section points, c = a + 0.382(b − a) and d = a + 0.618(b − a), are calculated first. By comparing the functional values at a, c, d, and b, we can determine the extremum either in [a, d] or [c, b] using the characteristic of the unimodal function. In either case, a new narrower search interval can be obtained. Successively narrowing the range of values, the extremum value will finally be obtained within its accuracy requirement.(2)
Genetic algorithm is a metaheuristic inspired by the process of natural selection. The initial population of the first generation is generated randomly, and then bio-inspired operators such as mutation, crossover, and selection are used to search the extremum. In each generation, individual solutions with higher values of fitness function are more likely to be inherited by the next generation. Repeat the process for several generations, the genetic algorithm then has a tendency to converge towards the global optima of the problem. See Ref. [24] for details.Table3 compares the optimal results with the enumeration method. One can see that the golden section method is the most efficiency and accurate one.Table 3
Optimal results.
Method
Genetic algorithms
Golden section method
Enumeration method
Computation time (s)
332.09
49.08
1830.01
Joint objective function
−4.4749
−4.4723
−4.4735
Performance availability
0.9835
0.9980
0.9985
Probabilistic resilience
0.998
1.00
1
System breakdown rate
0.00
0.00
0.00
Per-unit-time cost (yuan/hour)
￥58.92
￥64.46
￥64.96
Optimal PM interval (h)
108.8549
98.8276
92


### 4.3. Discussion

Using the enumeration method, we can obtain the simulation results at different PM intervals (Figure5). One can see that both the performance availability and probabilistic resilience decrease along with the increase of the PM interval in Figures 5(a) and 5(b). As the PM is an online maintenance activity, the more frequent inspections and maintenances help reduce the system’s performance degradation, resulting in better performance availability and probabilistic resilience. Figure 5(c) shows an obvious result, i.e., the system breakdown rate increases along with the value of the PM interval. This is because it is more convenient to discover the abnormal behaviour of the system before the system breaks down with a smaller PM interval. In Figure 5(d), the per-unit-time cost drastically decreases at the very beginning and then increases gradually along with the increase of the PM interval. If the PM interval is too small, excessive PM activities are triggered, leading to large on-site cost. If the PM interval is too large, PM is insufficient, resulting in more system breakdown and larger performance loss. This also illustrates the importance of choosing a suitable PM interval. The result is also consistent with that of the related studies in binary-state systems and discrete MSS.Figure 5
The four decision elements of the system under different PM intervals. (a) Performance availability. (b) Probabilistic resilience. (c) System breakdown rate. (d) Per-unit-time cost (yuan/hour).
(a)
(b)
(c)
(d)Figure6 shows the result of the joint objective function in equation (19). It illustrates more clearly the necessity of determining the optimal PM interval.Figure 6
The joint objective function under different PM intervals.


## 4.1. Resilience Behaviour Analysis

We run the Monte Carlo simulation and obtained the resilience process of the computer cluster. During the operation, system may suffer both single and continuous disturbances. The difference is whether the performance degradation of the system under the disturbance ends before a new one occurs. In Figure3, the first disturbance is a single one, and the other three are continuous ones which forms a disturbance group. In this figure, the four red stars show the time that the disturbance occurred. One can see that the performance degradation caused by the last three disturbances is superimposed on each other.Figure 3
Single and continuous disturbance effect.Figure4 shows two types of typical resilience behaviours for the computer cluster, and the PM interval was set as 900 hours and 2000 hours, respectively. In Figure 4(a), four disturbances occurred in sequence, resulting in the gradual load capacity degradation of four computers in the system. The performance of the system degraded along with that of the computers, as shown in this figure. At time t = 900 hours, the PM interval arrived while the system was still working, and PM was then taken to restore the system. The four degraded computers were repaired one by one according to the disturbance occurrence sequence. As the PM activity was online, the system’s performance restored back gradually in this figure. In Figure 4(b), the PM interval was set as 2000 hours, and the load capacity of the system dropped below its threshold QL at time t = 1212 hours. Then, the offline CM was taken, and the four computers were completely repaired before they restarted working. So one can see that the system performance recovery process in Figure 4(b) presents a step-like performance recovery process.Figure 4
Typical resilience processes of the cluster system. (a) Recover by PM (PM interval = 900 hours). (b) Recover by CM (PM interval = 2000 hours).
(a)
(b)


## 4.2. PM Interval Optimization

Firstly, a multiobjective optimization model was built to find the optimal PM interval as follows:(19)minf=k1CPUT+k2η−k3AP−k4ℝPs.t.AP≥0.9ℝPℝD≥0.8≥0.9η≤0.1,where ki is the weighting factor for each decision element. According to the direct weighing method, ki=ki1×ki2. Here, K1 = [0.6, 0.2, 0.1, 0.1]T described the importance of the four elements, and K2 was determined according to the magnitude difference among these four elements as K2 = [0.0018, 66.4896, 35.28, 10.2043]T. Compared with the optimization model in equation (17), the constraint of the per-unit-time cost is omitted, as there is no budget limit in this problem.Secondly, we used the Monte Carlo method to analyse the behaviour of the cluster system under different PM intervals. We run the simulation for 1000 times at the PM interval 1000 hours and obtained the corresponding sample standard deviations for the four decision elements. LetεC^PUT=5, εη^=0.002, εA^P=0.002, εℝ^P=0.002, and 1−α = 95%, we finally determined the number of iterations as N = 2700 according to equation (15). Under such 2700 simulation iterations, the error can be computed as εCPUT′<4.3627, εAP′<0.0019, εη′<0.0018, and εℝ^P′<0.0017.Finally, optimization algorithms, the golden section method and genetic algorithm, were used to solve the optimization problem in equation (19). The specific algorithms are as follows:(1)
Golden-section search is a technique for finding an extremum of a unimodal function inside a specified interval [a, b]. The functional values of fx at two golden section points, c = a + 0.382(b − a) and d = a + 0.618(b − a), are calculated first. By comparing the functional values at a, c, d, and b, we can determine the extremum either in [a, d] or [c, b] using the characteristic of the unimodal function. In either case, a new narrower search interval can be obtained. Successively narrowing the range of values, the extremum value will finally be obtained within its accuracy requirement.(2)
Genetic algorithm is a metaheuristic inspired by the process of natural selection. The initial population of the first generation is generated randomly, and then bio-inspired operators such as mutation, crossover, and selection are used to search the extremum. In each generation, individual solutions with higher values of fitness function are more likely to be inherited by the next generation. Repeat the process for several generations, the genetic algorithm then has a tendency to converge towards the global optima of the problem. See Ref. [24] for details.Table3 compares the optimal results with the enumeration method. One can see that the golden section method is the most efficiency and accurate one.Table 3
Optimal results.
Method
Genetic algorithms
Golden section method
Enumeration method
Computation time (s)
332.09
49.08
1830.01
Joint objective function
−4.4749
−4.4723
−4.4735
Performance availability
0.9835
0.9980
0.9985
Probabilistic resilience
0.998
1.00
1
System breakdown rate
0.00
0.00
0.00
Per-unit-time cost (yuan/hour)
￥58.92
￥64.46
￥64.96
Optimal PM interval (h)
108.8549
98.8276
92


## 4.3. Discussion

Using the enumeration method, we can obtain the simulation results at different PM intervals (Figure5). One can see that both the performance availability and probabilistic resilience decrease along with the increase of the PM interval in Figures 5(a) and 5(b). As the PM is an online maintenance activity, the more frequent inspections and maintenances help reduce the system’s performance degradation, resulting in better performance availability and probabilistic resilience. Figure 5(c) shows an obvious result, i.e., the system breakdown rate increases along with the value of the PM interval. This is because it is more convenient to discover the abnormal behaviour of the system before the system breaks down with a smaller PM interval. In Figure 5(d), the per-unit-time cost drastically decreases at the very beginning and then increases gradually along with the increase of the PM interval. If the PM interval is too small, excessive PM activities are triggered, leading to large on-site cost. If the PM interval is too large, PM is insufficient, resulting in more system breakdown and larger performance loss. This also illustrates the importance of choosing a suitable PM interval. The result is also consistent with that of the related studies in binary-state systems and discrete MSS.Figure 5
The four decision elements of the system under different PM intervals. (a) Performance availability. (b) Probabilistic resilience. (c) System breakdown rate. (d) Per-unit-time cost (yuan/hour).
(a)
(b)
(c)
(d)Figure6 shows the result of the joint objective function in equation (19). It illustrates more clearly the necessity of determining the optimal PM interval.Figure 6
The joint objective function under different PM intervals.


## 5. Conclusions

In this paper, we provide a general PM interval optimization framework for continuous MSS. In order to analyse the continuous change behaviour of the system’s performance, two types of performance-based measures, performance availability and probabilistic resilience, are proposed. Monte Carlo-based simulation is used to obtain the performance data of the system. The optimization model and algorithm proposed in this paper are flexible and can be adjusted according to users’ requirements. A case of computer cluster is applied to demonstrate the use of our method, and a multiobjective optimization model is built and solved by both golden section method and genetic algorithm. Comparing with the enumeration results, the effectiveness of our optimization method is verified.In practice, before the system is put into use, predicted data (including the probability of disturbance occurrence, performance degradation function, recovery function, etc.) can be used to determine the optimal PM interval. During operation, operational data of the system can be combined to find more appropriate PM interval. It should be noted that the decision elements in the PM interval optimization model and the model itself are different for BSS, discrete MSS, and continuous MSS. The decision elements, model, and algorithm in the PM interval optimization problem should be determined according to the system characteristics. If inappropriate ones are applied, the result obtained is likely to deviate from the optimal value.In the future, we will further study the optimal PM interval for heterogeneous system with continuous multistate components.



---

*Source: 2942940-2020-01-20.xml*