Application of Stochastic Optimization to Optimal Preventive Maintenance Problem

Petr Volf

Faculty of Science, Humanities and Education, The Technical University of Liberec, Liberec, Czech Republic.

DOI: 10.4236/jamp.2021.910157                

Abstract

The contribution deals with the optimization of a sequential preventive maintenance schedule of a technical device. We are given an initial time-to-failure probability distribution, model of changes of this distribution after maintenance actions, as well as the costs of maintenance, of a device acquisition, and of the impact of failure. The maintenance timing and, eventually, its extent, are the matter of optimization. The objective of the contribution is two-fold: first, to formulate a proper (random) objective function evaluating the lifetime of the maintained device relatively to maintenance costs; second, to propose a numerical method searching for a maintenance policy optimizing selected characteristics of this objective function. The method is based on the MCMC random search combined with simulated annealing. It is also shown that such a method is rather universal for different problem specifications. The approach will be illustrated on an artificial example dealing with accelerated lifetime after each maintenance action.

Keywords

Reliability, Preventive Maintenance, MCMC Algorithms, Simulated Annealing, Stochastic Optimization, Accelerated Lifetime Model

1. Introduction

The maintenance objective is as a rule to increase the lifetime of maintained device and to reduce the maintenance costs simultaneously. As these two requirements are in fact contradictory, one has to consider a criterion combining them, using for instance relative cost to lifetime unit, or, on the contrary, the proportion of effective lifetime to costs needed for the device acquisition and maintenance, not neglecting costs caused by eventual device failure.

The main contribution of the present study is the following: firstly, a quite general model of optimal preventive maintenance problem is formulated, including a rather universal scheme of sequential maintenance actions. Then, a procedure of numerical solution is proposed. It is based on the MCMC (Markov Chain Monte Carlo) method solving the stochastic optimization problem via controlled simulation of possible maintenance scenarios. Both these features of the approach are rather original (to author’s knowledge). The way of solution is demonstrated on artificial examples.

There exists a number of works dealing with the optimal maintenance problems from different points of view. Some are closely connected with concrete applications; others are more theoretical; however, most of them can serve at least as an inspiration for solutions of real problems. It should be also the case of the present paper. In order to put this to a wider research context, let us here present a brief selection of references relevant to the theme.

The paper [1] reviews computational approaches to maintenance scheduling. Several major categories are distinguished. One group of approaches is those utilizing traditional optimization-based techniques, such as dynamic programming, yielding exact optimal solution. Another type of approach, taking into account stochastic aspects of the problem, deals with heuristic methods, including simulated annealing, genetic algorithm, or guided randomized search methods. A similar overview of both the state of the art and directions of development is given in [2] . Both categorizations concern prevailingly to methods of solution; however, first, the problem description (model) has to be specified, concerning the lifetime distribution, possible actions, their consequences as well as costs.

Sometimes, the schedule is given (allowing for certain tolerance); the objective is to optimize distribution and allocation of sources and people to carry out all required actions at proper time, see for instance [3] . In the present paper, on the contrary, the time schedule and, eventually, the extent of maintenance are to be optimized, while their impact and costs are given (though depending, naturally, on maintenance actions). Schematically, the setting is sketched in Figure 1 below. As the lifetime (survival) of a device is random, characterized by its probability distribution, we deal with the problem of stochastic optimization. In general cases, it could be rather difficult to find one optimal solution. Sometimes a unique solution does not exist, or the surface of objective function is flat around its optimum. Thus, the intensive simulation (under selected strategy) is a common way of how to find at least sufficiently sub-optimal solution (see also [4] ).

More complicated problem arises in connection with optimal maintenance in complex systems. A crucial task there is to describe the dependence among system components and their mutual impact. Such a description could be based on the

Figure 1. Scheme of repeated maintenance actions.

7FTA analysis, which should contain also a model of dynamic reaction to events (defects, maintenances) of the system elements. One example of such an approach could be found in [5] . A nice literature overview concerning maintenance decision objective functions is presented in [6] . The authors also discuss the gap between maintenance impact and costs, leading to a multi-criteria optimization problem. A useful overview of problems connected with maintenance optimization is given also in [7] , it can serve as an inspiration to further deeper study of this area and offers also a number of relevant sources.

Another set of models studied in optimal maintenance literature, see for instance [8] , considers also inspections of actual system state. Then the maintenance is optimized in a sequential manner conditionally on the current prediction of system reliability. In the framework of our setting, inspections could be included as a special maintenance action, then the optimal maintenance policy has to be recomputed sequentially, when a new information on actual lifetime distribution is provided.

It is necessary to mention here also the models of imperfect repairs (we shall return to this theme in more detail in Section 3). The models are often based on the reduction of the hazard rate, either directly or indirectly (by shifting the virtual age of the system). If the state of the system is characterized by its degradation, the repair degree can be connected with the reduction of the degradation level. However, here other problems arise. First, how the degradation should be characterized and evaluated. Actually, the degradation could be connected directly with the growth of cumulated hazard rate, i.e. with one of basic characteristics of survival time probability distribution. There are numerous other possibilities. One of standard approaches consists in describing the process of degradation and repairs as a discrete (Markov or hidden Markov) chain on a finite set of states, see for instance [9] . However, the second problem is how to connect the state improvement with maintenance costs. This has to be settled as well.

In the present study, the specification of device degradation is not solved, as, in fact, it could be as well characterized by cumulated hazard rate, i.e. given by lifetime probability distribution. The impact of a maintenance is then modeled via certain change of this distribution. The organization of the paper is the following: In Section 2 a sketch of general model of maintenance sequence and its impact is provided and the objective function of stochastic optimization problem is formulated. Then, Section 3 recalls several models of non-complete maintenance and its impact to future lifetime. One of them, namely the model of accelerated ageing, is then described in more details. Its particular version is then specified to be used in concluding example. In Section 4, the method of numerical solution is described. It utilizes the random search via the Metropolis algorithm combined with simulated annealing. It is, in Section 5, applied to a solution of the optimal maintenance in an artificial example. Finally, the conclusion discusses how the approach can be adapted to different more complicated variants of the optimal maintenance problem.

2. Problem Formulation

Figure 1 shows a general scheme of a sequence of maintenance actions. The symbols used there have the following sense:

Notation:

$T_j,j=1,2,\cdots$ —times of maintenance actions, to be optimized, $T_0=0$.

$T_F$ —time of failure (random), it depends on maintenances performed before its occurrence.

$X_j,j=1,2,\cdots$ —survival times after j-th maintenance (random variables),

$X_0$ —initial time to failure (random variable).

$F_j,f_j,j=0,1,2,\cdots$ —cumulative distribution functions and densities of $X_j$.

${\bar{F}}_j=1-F_j$ —corresponding survival functions.

$C_j,j=1,2,\cdots$ —cost of j-th maintenance, $C_F$ —cost of failure,

$C_A$ —cost of device acquisition, $C_0=C_A+C_F$, C—total costs.

${\tau }_j=T_j-T_{j-1},j=1,2,\cdots$ —time periods between maintenances.

$Z=T_F/C$ —random objective function.

$M_j,j=1,2,\cdots$ —events, maintenance actions, F—event, failure.

While the initial distribution of $X_0$, i.e. $F_0$, is assumed to be known, other distributions $F_j,j=1,2,\cdots$ depend on maintenances, that is why the model of such a dependence should be specified. Moreover, in a general case also the costs of actions may depend on their extent. Further, as the problem is stochastic, we have to select an objective function, i.e. a corresponding random variable and its representative characteristics, e.g. its expectation, or its certain quantile. Let us first assume that all components and parameters of the scheme are fixed. Then it is possible, at least formally, to evaluate distribution of following random variables, which might be of the interest:

1) Distribution of the time to failure T_F: In each interval $\left(T_{j-1},T_j\right),j=1,2,\cdots$ its distribution is given (conditionally, provided the device survives $T_{j-1}$ ) by a probability density function $f_{j-1}\left(t\right)$ with $t\in \left(\mathrm{0,}{\tau }_j\right)$, and, simultaneously, the probability that the device survives over $T_j$ equals ${\bar{F}}_{j-1}\left({\tau }_j\right)$. Let us assume that there are exactly K maintenance times, fixed in advance (the device can but need not survive to them). Hence, when densities $f_j\left(t\right)$ are specified, we are able to derive distribution of lifetime $T_F$. In particular, its expectation $E{T}_F$ then equals:

$\begin{array}{l}E{T}_F={\displaystyle {\int }_0^{{\tau }_1}t{f}_0\left(t\right)\text{d}t}\\ +{\bar{F}}_0\left({\tau }_1\right)\left[{\displaystyle {\int }_0^{{\tau }_2}t{f}_1\left(t\right)\text{d}t}+T_1+{\bar{F}}_1\left({\tau }_2\right)\left[\cdots +{\bar{F}}_{K-1}\left({\tau }_K\right)\left[E{X}_K+T_K\right]\cdots \right]\right].\end{array}$ (1)

2) Distribution of costs, C: It is discrete; if the failure occurs in interval $\left(T_j,T_{j+1}\right),j=0,1,2,\cdots ,K-1$, then the total costs are $C={\displaystyle {\sum }_{m=0}^j}{C}_m$, after the last $T_K$ the costs equal $C={\displaystyle {\sum }_{m=0}^K}{C}_m$. Corresponding probabilities are then ${\displaystyle {\prod }_{m=0}^{j-1}}{\bar{F}}_m\left({\tau }_{m+1}\right)\cdot F_j\left({\tau }_{j+1}\right)$, while ${\displaystyle {\prod }_{m=0}^{K-1}}{\bar{F}}_m\left({\tau }_{m+1}\right)$, is the probability of survival over the last $T_K$.

3) Distribution of random variable $Z=T_F/C$, the time-to-failure to costs- to-failure. Again, when the model is specified, there is no problem to derive this distribution from above. As it was already said, we might be interested in maximization of its certain quantile, however, in the present paper, we shall search for parameters maximizing its expectation. For given all model parameters, the expectation EZ has a form quite similar to (1), each term of (1) has to be divided by corresponding costs shown in point 2. For instance, when the number of planned maintenances is just K = 2, then

$EZ=\frac{1}{C_0}{\displaystyle {\int }_0^{{\tau }_1}t{f}_0\left(t\right)\text{d}t}+{\bar{F}}_0\left({\tau }_1\right)\left[\frac{{\displaystyle {\int }_0^{{\tau }_2}t{f}_1\left(t\right)\text{d}t}+{\tau }_1}{C_0+C_1}+{\bar{F}}_1\left({\tau }_2\right)\left[\frac{E{X}_2+{\tau }_1+{\tau }_2}{C_0+C_1+C_2}\right]\right].$ (2)

It is seen that, in general case, the solution of the optimization task could be a rather demanding problem. Though the structure is sequential, it is not possible to solve the problem parts separately. For instance, it is rather difficult to find an optimal ${\tau }_2$ as a function of ${\tau }_1$. On the other hand, as shown above, for given all variables and parameters, the objective function can be evaluated. Hence, optimal parameters can be obtained by a convenient method of (randomized, possibly) search.

2.1. Sketch of Optimization Procedure

While a numerical example is solved in detail in Section 5, here the main features of the solution approach are described. The objective is to find the configuration of maintenance times $\left\{T_j,j=1,2,\cdots \right\}$ (and, eventually, the maintenance extent parameters) maximizing EZ. In the numerical example in Section 5 the corresponding parameter (namely the parameter of time acceleration, $\alpha$ ) is fixed, in a more general case such parameters could also be the subject of optimal search, then, naturally, the costs of maintenances have to depend on them, too.

The procedure starts from a reasonably chosen initial configuration and is searching for the best one, comparing their values EZ. It means that for each considered configuration corresponding EZ has to be computed. The way of search depends on selected search method. There exists a number of random search methods, each consisting, as a rule, of repeating two sub-steps: First, a new configuration is proposed, second, it is accepted or not. For instance, genetic algorithms can propose a new configuration as a combination of already accepted, they store a whole set of best results. On the contrary, the method used here is based on the Metropolis algorithm. The method creates a (Markov) chain of configurations, at each step both the proposal distribution and the acceptance probability depend just on the actual state and proposed new state. The algorithm accepts also non-optimal configurations, which helps the procedure to escape from the regions of local extremes. In order to force the procedure to converge to the global extreme, the method is combined with simulated annealing. Other details of the method are given in Section 4. As a rule, such a global search method should be started several times, in order to be more sure that really the global extreme has been approached.

When the search ends (and we hope that the optimal solution is achieved), it could be of interest to evaluate corresponding distributions of optimized variable Z, also of T_F and C, eventually their other characteristics, e.g. selected quantiles.

2.2. Illustrative Example

The objective of the following example is just to illustrate the form of results, i.e. the evaluation of mentioned distributions and also the maintenance times optimization, for a rather simple special case. Let the initial distribution of X₀ be Weibull with parameters a = 100 (scale), b = 2 (shape), i.e. its characteristics are $E{X}_0~89$, $\text{std}{X}_0~46$. Further, let us for simplicity assume that this distribution remains the same after all maintenances (i.e. the maintenances are “complete”), their maximal number is K = 10, inter-maintenance times are also the same, ${\tau }_j=55,j=1,2,3,\cdots ,10$. Further, $C_0=10$, all other $C_j=1,j=1,2,\cdots ,10$. Figure 2 shows distributions of all three variables, i.e. $T_F\mathrm{,}C\mathrm{,}Z$. Obtained characteristics of these variables (namely their expectations and standard deviations) then are: $E{T}_F=191.11$, $\text{std}{T}_F=179.86$, $EC=11.98$, $\text{std}C=3.36$, $EZ=10.86$, $\text{std}Z=7.36$, $\text{median}Z=9.66$.

In fact, $\tau =55$ is approximately the optimal inter-maintenance interval in this simplified setting. For comparison, let us consider the same case as above, except that the number of maintenances is limited to only 3, on the other hand we optimized all 3 inter-maintenance times (using already the method of random

Figure 2. Cumulative distribution function of time-to-failure T (above), the dotted curve there is the c.d.f. of Weibull(100, 2) r.v. X₀. The distribution of total costs C is in the middle, cumulative distribution function of variable Z is below.

search described later in Section 4, however assuming just complete repairs at maintenance times). It was found that the maximum of EZ was attained approximately for ${\tau }_1,{\tau }_2,{\tau }_3=41,50,66$ with $EZ=12.30$, while median $Z=11.94$.

2.3. A Variant Model with Finite Time Horizon

The time of model introduced in (1) has no upper bound, assuming that each finite lifetime can be exceeded with positive probability. While it holds theoretically, in practice it is not realistic. Therefore let us now assume that there is an upper time limit T_R, given as a rule by the manufacturer, after which the device is innovated completely (or exchanged for new). The model is then slightly different. For instance, let us consider again fixed number K of planned maintenances between 0 and T_R. Let us recall that C_A is the cost of the device acquisition, C_F the cost of failure consequences, $C_0=C_A+C_F$. Now, after a convenient re-arrangement, the expression for EZ is

$\begin{array}{c}EZ=\frac{{\displaystyle {\int }_0^{{\tau }_1}t{f}_0\left(t\right)\text{d}t}}{C_0}+{\displaystyle \underset{m=1}{\overset{K}{\sum }}\left[{\displaystyle \underset{j=1}{\overset{m}{\prod }}{\bar{F}}_{j-1}\left({\tau }_j\right)}\cdot \frac{{\displaystyle {\int }_0^{{\tau }_{m+1}}t{f}_m\left(t\right)\text{d}t}+T_m}{{\displaystyle {\sum }_{i=0}^m{C}_i}}\right]}\\ +\frac{{\displaystyle {\prod }_{j=1}^{K+1}{\bar{F}}_{j-1}\left({\tau }_j\right)}\cdot T_R}{{\displaystyle {\sum }_{i=1}^m{C}_i}+C_A},\end{array}$ (3)

where ${\tau }_j=T_j-T_{j-1}$ for $j=1,2,\cdots ,K$, ${\tau }_{K+1}=T_R-T_K$ is the time remaining to $T_R$ after the last maintenance action. Thus, in general, except the maintenance times $T_j,j=1,2,\cdots ,K$, also their number K should be optimized.

There are many model variants. For instance, one just slightly different from the setting considered here assumes the renewal after a given number of maintenances, similarly like in [10] . It means that K is given, at T_K the device is innovated, however the maintenance times T_m should be found optimally as well. Expression (3) is adapted easily to this case. The way of solution could be also quite similar to our solution here.

In the numerical example above it was assumed that after a maintenance the device is as new. This is evidently unrealistic in many instances. The model (3), on the contrary, is more universal and considers changed distributions. That is why the next section is devoted to a brief presentation of commonly used models taking into account the maintenance impact.

3. Models of Partial Repairs

The renewal (or the “complete repair”) means that the device is repaired fully (or exchanged for a new one) and that, consequently, the successive random variables—times to failure—are distributed identically and independently. On the contrary, the minimal repair means that the repair has no innovative effect (we can imagine just switching an electronic device off and on, to restart it).

There are several natural ways how the models of complete repairs can be extended to repairs incomplete. Some authors, e.g. [11] , speak in this context on “Generalized renewal processes” meaning that after the maintenance and other interventions the distribution of time-to-failure is changed. One of basic contributions is in [12] : Let $F_0$ be the distribution function of the time to failure of a new system. Assume that a maintenance reduces the virtual age. It could be modeled directly by switching from the actual age $V$ to $\delta \cdot V$, with a parameter $\delta \in \left[\mathrm{0,1}\right]$ characterizing the degree of repair. Hence, the remaining lifetime is from this moment described by the distribution function $F_1\left(t\right)=F_0\left(t+\delta V\right)/{\bar{F}}_0\left(\delta V\right),t\ge 0$.

In [12] , several sub-models of imperfect repairs are then specified. Denote by ${\delta }_n$ the degree of the n-th repair, $T_n$ the time survived after the last repair. Then in Model I the n-th repair cannot remove old damages, $V_n=V_{n-1}+{\delta }_n\cdot T_n$. On the contrary, the Model II allows for a reduction of the whole virtual age, namely $V_n={\delta }_n\cdot \left(V_{n-1}+T_n\right)$. Notice that special cases contain the complete repair model with ${\delta }_n=0$, minimal repair model, ${\delta }_n=1$, and frequently used variant with constant degree $\delta$. Naturally, there are many other approaches, e.g. considering a randomized degree of repair, the regressed degree (based on the system history), accelerated virtual ageing, change of hazard rate (i.e. its reduction as a rule) etc. A nice overview of models of partial maintenance can be found e.g. in [13] .

Accelerated Failure Time Model

In the rest of this paper, the specification of the model presented schematically in Figure 1 will use the Accelerated Failure Time (AFT) model assuming that after an event the subjective time speed is changed, as a rule accelerated. Formally, the maintenance returns the virtual age back to zero, however, the speed of ageing is then higher, by an accelerating parameter $\alpha >1$. In general, this parameter value, $\alpha ={\alpha }_j$, can depend on maintenance action $M_j$. Then, when the distribution function of the variable $X_{j-1}$ is $F_{j-1}\left(t\right)$, after the next j-th maintenance $M_j$ the distribution function of $X_j$ is $F_j\left(t\right)=F_{j-1}\left({\alpha }_j\cdot t\right)$ and its time starts from zero anew. Hence, for probability densities it means that $f_j\left(t\right)={\alpha }_j\cdot f_{j-1}\left({\alpha }_j\cdot t\right)$.

In the sequel, it will be assumed that the initial lifetime distribution of $X_0$ is given and that the acceleration parameter $\alpha$ is constant and known. It means that also the distributions of variables $X_j,j=1,2,\cdots$ are known, they do not depend on maintenance timing. If we further assume that all maintenance costs are the same, $C_j=C_1,j=1,2,\cdots$, the problem description seems to be rather simple. Nevertheless, still the maintenance schedule optimization is not less difficult and needs a convenient numerical procedure.

4. Method of Solution

The method used here will be based on the Markov Chain Monte Carlo approach (MCMC, see for instance [14] ) generating a chain of possibly sub-optimal solutions, and on its combination with the simulated annealing reducing variability of accepted sequence of solution and “forcing” it to approach the optimum. In this framework, two variants of optimization can be considered. First, optimal times $T_j$ are found for several fixed values K, then the optimal K from them is selected. The second method should be able to change K during the process of search. In the framework of the MCMC algorithms, the method is known as the reversible jumps MCMC. As the structure of corresponding algorithm is then more complicated, we shall use the first approach.

Hence, let us assume that K is fixed, we search for optimal vector of maintenance times $T=\left(T_1\mathrm{,}\cdots \mathrm{,}{T}_K\right)$ from the subset $\mathcal{T}=\left\{0<T_1<T_2<\cdots <T_K<T_R\right\}\subset R^K$. The model sketch in Figure 1 as well as its description in (3) suggest that a reasonable method could consists in a sequential alternating of times $T_j$. We shall follow the Metropolis algorithm. It starts from an initial configuration of $T$, e.g. selected equidistantly inside $\left(\mathrm{0,}{T}_R\right)$. At each step, one new component $T_j^{new}$ is proposed, selected randomly in $\left(T_{j-1}\mathrm{,}{T}_{j+1}\right)$, and accepted with probability

$P=\min \left\{1,\exp \left(E{Z}^{new}-E{Z}^{old}\right)\right\}.$ (4)

Hence we can control the distribution proposing new candidates, the most rough choice is the uniform distribution in $\left(T_{j-1}\mathrm{,}{T}_{j+1}\right)$, a finer selection could be beta distribution, or another distribution limited to this interval, with the last accepted value $T_j^{old}$ as its parameter (constrained Gauss distribution could be used as well). In every case, such a procedure innovating one component after another creates a Markov sequence of sets of times $T^{\left(n\right)}$, with limit distribution given by a probability density function proportional to $\exp \left(EZ\left(T\right)\right)$. The rule (4) accepts, with positive probability, also non-optimal values, therefore it is able to jump from a local maximum region. However, the objective is to find the global maximum of EZ. The method of simulated annealing suggests, instead of acceptance probability (4), to use its modification, namely

$P=\min \left\{\mathrm{1,}\exp \left(\frac{E{Z}^{new}-E{Z}^{old}}{Temp}\right)\right\}\mathrm{.}$ (5)

Then, the acceptance rule is controlled by a variable $Temp>0$ (“temperature”), it is seen that when $Temp$ is close to zero, acceptance probability of non-optimal solutions (i.e. the cases $E{Z}^{new}<E{Z}^{old}$ ) is reduced and closer to zero as well. The success of simulated annealing application depends on the selection of rate of convergence of $Temp$ to zero. For instance, $Temp\left(n\right)~1/\log \left(n+1\right)$, with n denoting the iteration loop number, suffices theoretically. However, for practical use it tends to zero too quickly, there does not exist a universal and exact rule how to adapt this rate to actual case.

5. Numerical Example of Optimal Maintenance

For numerical illustration we selected again the distribution of $X_0~\text{Weibull}\left(a=100,b=2\right)$, final time of renewal was set to $T_R=300$. Notice that without any maintenance $P\left(X_0>300\right)~0.00001$. Further, the acceleration parameter was $\alpha =1.1$. Thus, unlike in the illustrative example at the end of Section 2, the distribution changed, in the AFT model manner, after each action. The objective was to find an optimal scheduling of maintenances, i.e. times $T_1\mathrm{,}\cdots \mathrm{,}{T}_K$ (or, equivalently, inter-maintenance periods ${\tau }_1\mathrm{,}\cdots \mathrm{,}{\tau }_K$ ), and also optimal number K maximizing EZ, where random variable Z was now the proportion of $T=\min \left(T_F\mathrm{,}{T}_R\right)$ to costs spent to T, including $C_0=C_A+C_F$ if failure occurred before $T_R$, or just $C_A$ if device survived to $T_R$. The expression for EZ is written in (3), the method of randomized search for optimal solution is sketched roughly above. As it was said there, optimal times $T_j$ were found for different fixed K, then the best K was selected.

In the first experiment the costs were set to $C_A=5,C_j=1$ for maintenances $j=1,2,\cdots ,K$, the cost of failure impact $C_F=5$, too. Initial location of points $T_k,k=1,2,\cdots ,K$ was regular inside $\left(\mathrm{0,}{T}_R\right)$. Then, the iterative search followed the Metropolis algorithm described above, with N = 10,000 loops. In each loop the innovation of each point $T_k,k=1,2,\cdots ,K$ was proposed sequentially, candidates were selected randomly uniformly from an interval inside $\left(T_{k-1}\mathrm{,}{T}_{k+1}\right)$ (where $T_0=0$, $T_{K+1}=T_R$ ). It means that computation contained together $N\cdot K$ steps. New point $T_k$ was accepted along (5), with variable $Temp\left(n\right)=\exp \left(-c\cdot {\left(n-1\right)}^2\right)$. Here $n=\mathrm{1,}\cdots \mathrm{,}N$ was the number of loop, c was such that final $Temp\left(N\right)=0.001$ (i.e. $c~7\times {10}^{-8}$ ). Notice that such a choice yields an “S-curve” decreasing from 1 to 0 (for instance $Temp\left(5000\right)~0.19$, in fact it is the survival function of Weibull distribution). Its shape allows for rough allocation of optimal configuration in the first part of iteration, making it more accurate in the second part. Simultaneously, also the width of proposal interval was decreased, which resulted in reduction of number of rejected proposals.

The best result was found for $K=5$. Optimal times of maintenance were $T_1,\cdots ,T_5=53.9,108.3,162.7,2145.1,262.3$. Achieved maximal $EZ=11.7$. Figure 3 shows the distribution of costs and the distribution of variable Z in optimal case. Then, Figure 4 presents, in its upper part, the trace of search for maximal value of EZ, while its lower part shows how location of points $T_k$ was concentrated around their optimal values. Both graphs display values from the second part of iteration, i.e. last 5000 loops, 25,000 iterations. It must be said that the shape of objective function around its optimum was rather flat, maximal values of EZ for selections $K=4$ or $K=6$ were just slightly smaller (by several %), as well as the result for quite regular inter-maintenance periods.

In order to explore a solution sensitivity, namely the change of optimal solution for different costs, in the second case we set $C_F=10,C_j=1$ for maintenances, while $C_A=5$ remained unchanged. Now, in order to reduce the probability of rather expensive failure event, optimal number of maintenances was higher, namely $K=6$, optimal times then $T_1,\cdots ,T_6=55.00,106.82,156.84,205.22,241.47,274.54$, maximum of $EZ=9.17$.

On the contrary, the third case assumed that there were no additional failure expenses. We selected $C_F=0$, $C_A=C_0=10$ and $C_j=1$ for maintenances.

Figure 3. Probability distribution of total costs C (upper plot), circles correspond to costs with failure, the star sign at value 10 to costs of all maintenances plus C_A (i.e. without failure). Below is the cumulative distribution function of variable Z.

Figure 4. Upper plot shows the sequence of EZ, lower plot the histogram of accepted locations of times T_k. Both from the second half of iterations.

As expected, optimal number of maintenances was lower, namely $K=4$, optimal maintenance times were $T_1,\cdots ,T_4=47.45,99.21,152.91,206.92$, and achieved maximum of EZ was 10.93.

As shown in Figure 3, for given K and maintenance times $T$ it is possible to derive the distribution function of variable Z, therefore the task of maximization of its another characteristics (e.g. of a selected quantile) can be solved rather straightforwardly. Recall that the p-quantile of a random variable Z, $Q_Z\left(p\right)$, fulfils $P\left(Z\ge Q_Z\left(p\right)\right)\ge 1-p$. In other words, it means that Z is larger than its quantile value with guaranteed probability. Let us consider again the 1-st case, i.e. the case with $C_A=C_F=5$. The maximal quantile problem was solved with the aid of the same algorithm as the search for maximal expectation. The result for maximal median of Z was not far from optimal solution for EZ, with $\text{max}\left(\text{med}Z\right)=11.239$ achieved now for $K=6$. Then, for instance, the 25% quantile of Z was maximized for $K=7$, optimal $T=30.27,54.69,76.25,101.49,127.07,148.10,187.61$ and achieved max $Q_Z\left(0.25\right)=7.239$.

As the last example, let us consider a variant where two types of preventive maintenance are available. Let the parameters of the case be the same as above, namely the initial distribution being Weibull(100, 2), $C_F=C_A=5$, $T_R=300$, and let there be a choice between two actions, the first accelerating the lifetime by ${\alpha }_1=1.1$, with cost $C_1=1$ (i.e. as above), the second cheaper, with cost $C_2=0.5$, however with the lifetime accelerated more, with ${\alpha }_2=1.2$. The method of random search via the Metropolis algorithm with simulated annealing, proposing also types of actions, yielded an optimal solution consisting of $K=6$ maintenances in $\left(\mathrm{0,}{T}_R\right)$ and combining both maintenance types. Namely, the procedure found a solution with $K=6$, optimal sequence of times 51.83, 105.21, 161.21, 212.10, 251.53, 279.43, corresponding maintenance action types were 2, 2, 1, 1, 2, 2, achieved $EZ=11.91$. This is more than in the case with just 1-st action (results are above), and also when just the 2-nd type is used (in this case, maximal $EZ=11.42$ was achieved, for $K=7$ ).

6. Conclusions

In the present paper, a rather general model of sequential maintenance actions, including their effects and costs, was sketched, and the problem of stochastic optimization of action times was formulated. Further, the method of numerical solution, in fact also rather general, was proposed. Different approaches to maintenance impact characterization were then described and discussed. Finally, in order to demonstrate the way of solution, the model was specified, the method applied and illustrated on an artificial example.

The approach proposed in the present contribution can easily be adapted to different formulations of the maintenance problem. The modification can concern to following aspects (among others):

1) Modeling the reliability via different probabilistic distributions of time-to- failure. For reliability of one unit, the standard choices are the Weibull (exponential as well), or the log-normal, gamma, Gumbel and Gauss distribution; sometimes their combination is applied. However, even the empirical distribution based on lifetime data directly, i.e. yielding a stepwise reliability function, can be used.

2) Impact of maintenance action is, in general, described via a change of failure rate. Two ways on how this could be formulated are presented in Section 3. Naturally, other descriptions are possible.

3) Maintenance actions scheduling. Standardly, regular inter-maintenance periods are considered (as in [15] ), in such a way this part of problem is simplified and one can concentrate to other aspects of optimal solutions (e.g. the extent of maintenance). It has to be said that the solution assuming constant periods could be, as a rule, close to the optimal solution without this assumption. Still, even a simple example in subsection 2.2 shows that there is a difference.

4) Formal objective of optimization. Its choice depends on preferences. In the present contribution, the proportion of time-to-failure to spent costs is considered, its mean value taken as a criterion. Use of selected quantile is discussed, too. A variant approach considers multiple criteria (see, for instance, [16] comparing several scenarios). As a rule, in the end, it is necessary to select weights of each criterion and to optimize a combination of them.

5) Optimization with constraints (for instance concerning the budget or other sources for maintenance actions). This is another, rather practical approach to multi-objective optimization, where an optimal solution is found in restricted space of “strategies”. Most of randomized optimization procedures (including the one presented here) can control its search to stay inside prescribed space (the simplest way consists in rejection of inadmissible configurations).

6) Maintenance of a system of components, the theme of many studies. The key point of presented approach (see Figure 1) is the knowledge (i.e. a good model) of the time-to-failure probabilistic distribution of individual components before and after each maintenance action. Then, the distribution for the whole system failure can be recomputed, provided its proper scheme is given. Again, optimal times and also the type of actions should be selected, by the action type we mean here what component (and, eventually, to which degree) is to be maintained.

In fact, in a system with parallel structure (including “K of N” system), the failures of individual components should be taken into account rather than the failure of the whole system. The system is still available after one component failure; however, its capacity is reduced during its repair, which causes additional costs. Further, the number of redundant components could also be a matter of optimal choice (cf. [17] ).

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Oke, S.A. (2004) Maintenance Scheduling: Description, Status and Future Directions. The South African Journal of Industrial Engineering, 15, 101-117. https://doi.org/10.7166/15-1-241
[2]	Sharma, A., Yadava, G.S. and Deshmukh, S.G. (2011) A Literature Review and Future Perspectives on Maintenance Optimization. Journal of Quality in Maintenance Engineering, 17, 5-25. https://doi.org/10.1108/13552511111116222
[3]	Senra, P., Lopes, I. and Oliveira, J.A. (2017) Supporting Maintenance Scheduling: A Case Study. Procedia Manufacturing, 11, 2123-2130. https://doi.org/10.1016/j.promfg.2017.07.342
[4]	Lahiani, N., El Mhamedi, A., Hani, Y. and Triki, A. (2014) Simulation Based Optimization Approach to Solve a Maintenance Process Problem. Proceedings of 2014 International Conference on Control, Decision and Information Technologies (CoDIT), Metz, 3-5 November 2014, 146-151. https://doi.org/10.1109/CoDIT.2014.6996884
[5]	Turan, O., Lazakis, I., Judah, S. and Incecik, A. (2011) Investigating the Reliability and Criticality of the Maintenance Characteristics of a Diving Support Vessel. Quality and Reliability Engineering International, 27, 931-946. https://doi.org/10.1002/qre.1182
[6]	Hatsey, N.H. and Birkie, S.E. (2021) Total Cost Optimization of Submersible Irrigation Pump Maintenance Using Simulation. Journal of Quality in Maintenance Engineering, 27, 187-202. https://doi.org/10.1108/JQME-08-2018-0064
[7]	Van Horenbeek, A., Pintelon, L. and Muchiri, P. (2010) Maintenance Optimization Models and Criteria. International Journal of System Assurance Engineering and Management, 1, 189-200. https://doi.org/10.1007/s13198-011-0045-x
[8]	Wu, L. and Zhou, Q. (2021) Adaptive Sequential Predictive Maintenance Policy with Nonperiodic Inspection for Hard Failures. Quality and Reliability Engineering International, 37, 1173-1185. https://doi.org/10.1002/qre.2788
[9]	Kim, M.J. and Makis, V. (2010) A Maintenance Model with Minimal and General Repair. IMA Journal of Management Mathematics, 21, 385-393. https://doi.org/10.1093/imaman/dpn016
[10]	Castro, I.T. (2009) A Model of Imperfect Preventive Maintenance with Dependent Failure Modes. European Journal of Operational Research, 196, 217-224. https://doi.org/10.1016/j.ejor.2008.02.042
[11]	Tanwar, M., Rai, R.N. and Bolia, N. (2014) Imperfect Repair Modeling Using Kijima type Generalized Renewal Process. Reliability Engineering and System Safety, 124, 24-31. https://doi.org/10.1016/j.ress.2013.10.007
[12]	Kijima, M. (1989) Some Results for Repairable Systems with General Repair. Journal of Applied Probability, 26, 89-102. https://doi.org/10.2307/3214319
[13]	Pham, H. and Wang, H. (1996) Imperfect Maintenance. European Journal of Operational Research, 94, 425-438. https://doi.org/10.1016/S0377-2217(96)00099-9
[14]	Gamerman, D. and Lopes, H.F. (2006) Markov Chain Monte Carlo, Stochastic Simulation for Bayesian Inference (2nd edition). Chapman and Hall, London. https://doi.org/10.1201/9781482296426
[15]	Coria, V.H., Maximov, S., Rivas-Dvalos, F., Melchor, C.L. and Guardado, J.L. (2015) Analytical Method for Optimization of Maintenance Policy Based on Available System Failure Data. Reliability Engineering and System Safety, 135, 55-63. https://doi.org/10.1016/j.ress.2014.11.003
[16]	Linnéusson, G., Ng, A.H.C. and Aslam, T. (2018) Quantitative Analysis of a Conceptual System Dynamics Maintenance Performance Mode Using Multi-Objective Optimisation. Journal of Simulation, 12, 171-189. https://doi.org/10.1080/17477778.2018.1467849
[17]	Mellal, M.A., Zio, E. and Williams, E.J. (2020) Cost Minimization of Repairable Systems Subject to Availability Constraints Using Efficient Cuckoo Optimization Algorithm. Quality and Reliability Engineering International, 36, 1098-1110. https://doi.org/10.1002/qre.2617