Really Ageing Systems Undergoing a Discrete Maintenance Optimization

: In general, a complex system is composed of different components that are usually subject to a maintenance policy. We take into account systems containing components that are under both preventive and corrective maintenance. Preventive maintenance is considered as a failure-based preventive maintenance model, where full renewal is realized after the occurrence of every n th failure. It offers an imperfect corrective maintenance model, where each repair deteriorates the component or system lifetime, the probability distribution of which gradually changes via increasing failure rates. The reliability mathematics for unavailability quantiﬁcation is demonstrated in the paper. The renewal process model, involving failure-based preventive maintenance, arises from the new corresponding renewal cycle, which is designated a real ageing process. Imperfect corrective maintenance results in an unwanted rise in the unavailability function, which can be rectiﬁed by a properly selected failure-based preventive maintenance policy; i.e., replacement of a properly selected component respecting both cost and unavailability after the occurrence of the n th failure. The number n is considered a decision variable, whereas cost is an objective function in the optimization process. The paper describes a new method for ﬁnding an optimal failure-based preventive maintenance policy for a system respecting a given reliability constraint. The decision variable n is optimally selected for each component from a set of possible realistic maintenance modes. We focus on the discrete maintenance model, where each component is realized in one or several maintenance mode(s). The ﬁxed value of the decision variable determines a single maintenance mode, as well as the cost of the mode. The optimization process for a system is demanding in terms of computing time because, if the system contains k components, all having three maintenance modes, we need to evaluate 3 k maintenance conﬁgurations. The discrete maintenance optimization is shown with two systems adopted from the literature.


Introduction
Various maintenance strategies have recently been intensively studied and developed to improve the reliability, availability and usability of relevant industrial systems. Unexpected failures may cause dangerous situations for human lives, unplanned production outages, etc. This is why systems must be protected against them. System reliability can be significantly improved by applying the optimal maintenance strategy. One can read in the work of [1] that maintenance is no longer a necessary evil and that production companies should invest in maintenance to maximize revenues. To find the best and most suitable strategy for every unit or subsystem, preferred economic decisions have to be made that make it possible to achieve a profit. Therefore, evaluation of the performances of different strategies is often used ( [1,2]). We can distinguish between the different maintenance strategies that can be used: corrective maintenance (CM) and preventive maintenance (PM) strategies. CM strategies, often known as restoration or repair strategies, are only launched when a failure occurs and the system is broken. The system is then returned to a functioning state by applying a maintenance action. PM strategies aim to prevent the system from undergoing undesired breakdowns. PM is usually carried out on operating systems and it reduces ageing processes, which means that the probability of system failures is decreased. Maintenance modelling is a dynamically developing recent scientific discipline, and we do not aim to present a general overview of references on maintenance here. However, several papers can be discussed, some of them dealing with CM, others with PM. For example, authors of [3] mention that there are several different types of PM and CM actions depending on the degree of restoration, including perfect maintenance action, which restores the system to an as-good-as-new function state; and imperfect maintenance action, which can have several restoration levels between perfect maintenance and minimal maintenance, with the latter restoring the system to the same state it was in just before the failure occurred, with the same failure rate that it had just before the minimal maintenance action.
A repair action can gradually restore a system to its initial level, with the system being returned step-by-step to the operating state, resulting in perfect CM. This approach is in detail introduced for example in [4] and the authors call it a gradual CM repair strategy. In general, we can say that imperfect maintenance is a kind of intervention where the system is returned to someplace between being as good as new and as bad as it used to be. However, for any maintenance intervention, it is necessary to establish the level of imperfect maintenance. Authors in [5] introduced such imperfect restoration when they assumed that the system age is affected by the kind of maintenance intervention. Apart from the system age, the failure rate may also become worse due to maintenance actions ( [6]). Both age shortening and failure rate adaptation are assumed in the so-called hybrid model, which seems to be a more realistic model (see more in [7,8]). Whereas the former authors addressed the optimal maintenance decision for binary systems to maximize the reliability of the next mission under imperfect maintenance, the latter introduced the imperfect maintenance model, where fixed maintenance action corrects a system functionality to any state between minimal repair and full renewal; i.e., perfect repair. For extensive discussions regarding imperfect maintenance models, readers may refer to authors in [9,10].
In this article, we use the imperfect CM process, which degrades the system lifetime at any CM intervention due to the growing failure rate. Gradual changes in the failure rate result in corresponding changes in the lifetime probability distribution. The basic aim of this article is to describe a realistic failure-based PM model that undergoes a process of discrete optimization. PM is considered a failure-based preventive maintenance model (FBM), where full renewal is realized at the occurrence of every nth failure. This model is related to the imperfect CM process, an overview of which is provided in [11]. For instance, a comparable strategy, where a unit is replaced at the nth failure and (n − 1) previous failures are repaired with minimal repair, was proposed in the cost-focused study [12]. However, the author supposed that the failure rate would not be violated after carrying out minimal repair. Stochastic models that describe the failure pattern of repairable units subject to minimal maintenance are discussed in [13]. The unit can be replaced at time T or at the nth failure, whichever occurs first, and n can be minimized in the context with both repair and replacement costs, as is developed in [14]. We use the imperfect CM model, where each repair action deteriorates the component or system lifetime, the probability distribution of which is gradually changed via the increasing failure rate. This model is of particular interest for experts analyzing components of power distribution networks ( [15]), the reliability of which will be studied by the authors of this article in the near future.
Innovative reliability mathematics for the unavailability quantification of components is presented in this article. The new renewal process model, involving the FBM, is designated as a real ageing process. The imperfect CM results in an unwanted rise in the unavailability function, which can be rectified with a properly selected FBM process. This means that the renewal of a component starts with the occurrence of the nth failure. The number n is considered a decision variable, whereas cost is an objective function in the optimization process.
This article describes a new method that can be used to find an optimal FBM strategy to solve a particular optimization problem while respecting a given reliability constraint. The above-mentioned decision variable determining the different maintenance modes of a system component is optimally selected from a set of possible realistic maintenance modes. Thus, the discrete maintenance model is considered, where each component can work in one or several maintenance modes. The fixed value of the decision variable n determines one maintenance mode of the component, which predetermines both the unavailability course and cost. Different maintenance modes of system components result in different system configurations, each having a specific unavailability course, as well as cost. The optimization process is demanding in terms of computing time because a complex system can have many maintenance configurations. The discrete maintenance optimization is demonstrated with two systems adopted from the literature.

Optimization Problem
Each optimization problem works on the presumption that an objective function f (x) that varies in a given range must be optimized; i.e., either maximized or minimized, constrained by several restrictions imposed on the decision variables. The optimization problem in this article can be formulated using the following objective function f (x), which represents the minimum cost: where x = (x 1 , . . . , x k ) ∈ R k is a decision variable, and k is the number of system components each having the decision variable x i = n i , which can be optimized if needed. Each component undergoes a real ageing process, including imperfect CM, until the occurrence of the n i -th failure, which starts its restoration (renewal). In most cases, both f (x) and U S (x) are complicated, either linear or nonlinear, functions of the decision variable vector: x = (x 1 , . . . , x k ) = (n 1 , . . . n k ) that constitute parameters for which optimal values must be found.

Discrete Maintenance Model
Maintenance optimization can be classified in different ways. A recent thorough classification can be found in [16]. Concerning the optimization outcome and decision variables, the discrete maintenance model can be classified as a process that finds optimized parameter values defining a single maintenance strategy selected a priori; e.g., in this paper, the type of action performed (repair, replacement). There are different optimization approaches that take into account the previously selected decision variables. The methodology used in this paper can be included among the mathematical approaches in which the optimization problem is formulated utilizing mathematical equations, which are then solved using differential calculus to identify optimal parameters for the maintenance strategy.
We introduced a discrete maintenance model for real multi-component systems with non-identical components for the first time in [17], in which systems with repairable components and latent failures were optimized using changeable periods of inspection as a decision variable. Real complex systems are composed of a finite number of repairable and maintained components. Each component can be operated in different discrete maintenance modes. A maintenance mode of the i-th component is determined by a prescribed value of the decision variable (x i = n i ) that influences the maintenance cost of the mode. Given that a system is composed of k components and each component has four maintenance modes in general, we have to investigate 4 k maintenance configurations of the system. Any system configuration can be described by a maximal system unavailability U S (x) and a total cost C S , which is usually obtained as a sum of the costs of all component modes forming the configuration. The optimal system configuration is detected under requirements (1) and (2). This maintenance model is defined in this article as a discrete maintenance model. Similar models of maintenance optimization have been used in other publications. For example, the work of authors in [18][19][20] addresses the optimization problem under maintenance policies including two decision variables. One of them is the maximum number of failures before the system undergoes a perfect restoration and the second is the inspection interval used to detect hidden failures. In [21] authors consider only one component system subject to two types of age-dependent failure. Catastrophic failures are detected through periodic inspections, whereas minor failures are followed by minor repairs. The system is preventively replaced either at an optimal multiple of the inspection time or after the n-th minor failure, whichever comes first. Both parameters are decision variables. The cost per unit of time is an objective function. The objective is to obtain a cost-minimizing policy. Authors in [22] also present a maintenance model for a system subject to two types of unrevealed failures: minor and catastrophic. The system is replaced at the occurrence of the n-th minor failure, a catastrophic failure, or due to working age, etc.

The Method for Unavailability Calculation for a Complex System
The basic methodology, including algorithms, for unavailability calculation for a complex system with maintenance was developed in [23,24]. The system structure and its functionality are described through the use of a directed acyclic graph (AG). An AG contains nodes and edges. The system function or non-function state constitutes the highest node, which is at the top of the AG. Subsystems (components) are described through internal (terminal) nodes, all of which are interconnected by edges. An AG cannot contain feedback loops because it is acyclic. Terminal nodes characterize the stochastic behavior of the input components of a system, which means that each terminal node must be provided with information about the probability distribution of its lifetime, as well as maintenance characteristics and parameters. Stemming from this information, the unavailability function of each terminal node is further computed through the renewal process model described in the following section. The final system unavailability function is obtained through the unavailability functions of all terminal nodes. Other details concerning the computing algorithm that calculates the system unavailability function from component (i.e., terminal node) unavailability functions are described in our previous research work in ([17], pp. 86-87).

Unavailability Analysis of Terminal Nodes
The reliability mathematics used in this section results partially from the work of authors in [25] and partially from other articles presented in the references ( [26][27][28][29][30]). The mathematics used in these sources was developed to a large extent to consider the new renewal cycle of a terminal node, as introduced below. A complex maintained system consists of particular components that, in the context of the AG system structure, are denoted terminal nodes. In this section, the unavailability function of a terminal node is investigated. The renewal cycle of a terminal node starts at time t = 0. The first failure occurs at the time X 1 and CM of the node starts immediately. This CM action takes the time Y 1 . The next failure occurs after the time X 2 after the end of the repair of the first failure elapses. In this way, the renewal cycle continues until the nth failure occurs. This is followed by the FBM replacement of the node. When this is completed, the renewal cycle ends, as is demonstrated in Figure 1, where the length of the renewal cycle is T. We will call the progress over one renewal cycle a real ageing process. This evolution is characterized by gradually changing probability distributions of random variables X 1 , X 2 , . . . X n due to different degradation processes caused by CM. h failure occurs. This is followed by the FBM replacement of the node. When this is completed, the renewal cycle ends, as is demonstrated in Figure 1, where the length of the renewal cycle is T. We will call the progress over one renewal cycle a real ageing process. This evolution is characterized by gradually changing probability distributions of random variables X1, X2, … Xn due to different degradation processes caused by CM.
First, the unavailability function of the terminal node at a given time in the first renewal cycle is determined. Figure 1. The first renewal cycle. The time from the beginning of the renewal cycle to the occurrence of the second failure is a random variable * = + + , etc. The renewal cycle T is terminated by the n-th failure and followed by the replacement. The replacement time is the random variable .
The random variables , … , , , … , are supposed to be independent. The terminal node is out of service at time only if it is in the process of being repaired (or replaced) after the i-th failure. This happens only when 0 ≤ * ≤ ≤ * + for some ∈ {1, … , }.
It is well-known that Let us denote the probability that the terminal node is unavailable at time due to its CM or FBM after the -th failure ( ). It follows that The probability that the node is out of order at time in the first renewal cycle (let us denote it ( )) fulfils the following equation: The length of the -th renewal cycle is the random variable and the -th renewal cycle is the time interval ( , ⟩ where = 0 and is the value of the random variable = ∑ , as is shown in Figure 2.
The renewal cycle T is terminated by the n-th failure and followed by the replacement. The replacement time is the random variable Y n .
First, the unavailability function of the terminal node at a given time t in the first renewal cycle is determined.
The terminal node is out of service at time t only if it is in the process of being repaired (or replaced) after the i-th failure. This happens only when It is well-known that Let us denote the probability that the terminal node is unavailable at time t due to its CM or FBM after the i-th failure P i (t). It follows that The probability that the node is out of order at time t in the first renewal cycle (let us denote it u 1 (t)) fulfils the following equation: The length of the k-th renewal cycle is the random variable T k and the k-th renewal cycle is the time interval ( s k−1 , s k where s 0 = 0 and s k is the value of the random variable S k = ∑ k i=1 T i , as is shown in Figure 2. The unavailability function of a terminal node describes the probability that the node is out of order at a given time t. It can be determined as follows. Let us denote this probability ( ). It is obvious that where ( ) is the probability that the node is out of order at time , which belongs to the -th renewal cycle ( , ⟩. The unavailability function of a terminal node describes the probability that the node is out of order at a given time t. It can be determined as follows. Let us denote this probability u(t). It is obvious that where u k (t) is the probability that the node is out of order at time t, which belongs to the k-th renewal cycle ( s k−1 , s k . The value u 1 (t) is given by Formula (3). The values u k (t) for k = 2, 3, . . . can be determined in the following way, assuming that the FBM restores the terminal node to an as-good-as-new state. It holds that . . and, in general, for k ≥ 2 (see Figure 3): Under the assumption that the FBM restores the node to an as-good-as-new state, it also holds that The unavailability function of a terminal node describes the probability that the node is out of order at a given time t. It can be determined as follows. Let us denote this probability ( ). It is obvious that where ( ) is the probability that the node is out of order at time , which belongs to the -th renewal cycle ( , ⟩. The value ( ) is given by Formula (3). The values ( ) for = 2, 3, … can be determined in the following way, assuming that the FBM restores the terminal node to an as-good-as-new state. It holds that ( ) = ( − ), ( ) = ( − ), … and, in general, for ≥ 2 (see Figure 3): where is the cumulative distribution function of the random variable = ∑ . Under the assumption that the FBM restores the node to an as-good-as-new state, it also holds that = = for all , ∈ . This means that is the ( − 1)-fold convolution of the cumulative distribution function , where the random variable = ∑ + ∑ is the length of the renewal cycle. Hence
Equations (4) and (5) then imply: Equations (4) and (5) then imply: where Hence, the function G(s) is the solution of the integral equation Thus, knowing G(s), the required unavailability function u(t) can be determined as a solution of the integral Equation (6).
The starting point for the development of this methodology lies in our previous work introduced in [31], where dormant systems under inspection were investigated. The main difference is in the renewal cycle. A dormant system has a different renewal cycle because failures are not detected immediately but only at special inspection times.

Cost Model
The cost model of a system configuration can be derived by adding up all the contributions arising from both the repair and replacement processes of a mode over all of the system components. A maintenance mode of a component has two main cost contributions: the cost of FBM, given by the replacements of the component during a mission time T M that depends on the decision variable of the component; and the cost of the imperfect repair process, which further depends on the mean number of failures during the mission time T M and the CM parameters. In practical situations, the cost contributions result from a database for the year and give an average yearly cost for the system configurations in a monitored period. In the remainder of this article, the cost will be computed in non-identified cost units based on the summation principle.
To obtain the cost of one system configuration, we simply add up the costs of all maintenance modes of all system components. The mean cost of one maintenance mode of the j-th component C T M (j) can be computed as follows: where C CM (j) . . . CM cost = cost of one repair action for the j-th component C T M (j) . . . mean cost of one maintenance mode of the j-th component It is worth noting that the mean maintenance cost of the j-th component C T M (j) depends on the decision variable n j , not only directly, as follows from Formula (7), but also indirectly via n R (j) (see Formula (8) for MTTF(j)). The total cost of one system configuration C S is obtained by summing up these contributions, described by Formula (7), over all of the system components k: A similar principle for the computation of the total maintenance cost for the whole of a system was used in [32].

One Maintained Component System Adopted from the Literature
First, we take into account one component system that is under both PM and CM. PM is considered as FBM, where full renewal is realized at the occurrence of every nth failure. The imperfect CM model causes a real ageing process, where each CM intervention deteriorates the system lifetime at an increasing failure rate. Proceeding from the previously derived reliability mathematics, the system unavailability function can be determined from Formula (6). The imperfect CM results in an unwanted rise in the unavailability function that can be reduced by the properly selected FBM process, meaning that renewal of the system starts with the occurrence of the n-th failure. After the renewal, the system is restarted to an as-good-as-new state. The number n is considered here as a changing decision variable permitting optimization. The exact specification of the decision variable determines a system configuration that is connected to a specific cost computed according to Formula (9), where k = 1.
The discrete maintenance optimization method is shown for a system adopted from [25] where a stochastic alternating renewal process model is derived but a different renewal cycle is considered. The distribution function of the first failure time X 1 is a Weibull distribution with the shape parameter β = 2 and scale parameter α = 600 days. Imperfect CM is characterized by a random repair time with rectangular distribution in an interval of <12, 16> days, so that the mean time for CM is 14 days. We further suppose that the replacement time in the FBM model is deterministic-the renewal duration is 7 days; i.e., it is shorter than the CM time since it may be scheduled beforehand, which is given by the fact that the repair team has at its disposal information about (n − 1)-th failure.
The real ageing process can be realized in the following way. If a failure occurs, it is followed by a standard CM action, which recovers the health of the system to some extent, but its failure rate is always worse when compared with the system health before failure. We presume that, following the first system failure and the CM action, the growth of the failure rate can be estimated by the quotient q a , which ranges from 1 to 1.5 and by which the failure rate is multiplied. Worsening of the failure rate will continue after the second failure and will be followed by a repair time, etc., until the time of the FBM intervention; i.e., the time of the n-th failure. Provided that the first failure time X 1 follows a Weibull distribution, the failure rate can be expressed in the following way: X 2 also has a Weibull distribution failure rate, which is worsened as follows: Mathematics 2022, 10, 2865 9 of 17 and the worsening continues with each subsequent failure in the same way. In the time X n−1 , the failure rate can be expressed as: λ(t) = β α 2 .q a n−2 .t (12) and the FBM is started at the time X n , which launches a renewal (replacement). Thus, the system progressively deteriorates, and its lifetime distribution is modified after each CM action until an FBM time comes. This is in accordance with real practice because each failure followed by a CM action exposes the system to shocks that accumulate, and the system becomes increasingly worse. For example, the level of deterioration after the CM of circuit breakers, components of a power distribution network, has been intensively and personally discussed with experts and authors ( [15]), who stated that each CM action makes the initial technical health of the component worse by approximately 10-25%. That is why we selected the following computing experiments with an appropriate value for the quotient q a = 1.25, which means the limit worsening of the failure rate by 25% after each CM intervention. The difference between real and theoretical ageing is demonstrated in Figure 4, which illustrates the time-dependent unavailability function of the system u(t), computed according to Formula (6), with a mission time of 4000 days for two system modes:

1.
A real ageing mode without FBM; i.e., n = ∞ and ageing quotient q a = 1.25, so that each failure is followed by CM and the subsequent system lifetime has an increasingly worse failure rate, in accordance with Formulas (10)-(12); 2.
A traditional ageing mode, where ageing is due to an increasing Weibull failure rate. This ageing process can be denoted as theoretical ageing, where each failure is followed by a standard CM intervention that makes the system as good as new, such that the system is replaced (n = 1).
Mathematics 2022, 10, x FOR PEER REVIEW 10 of 18 In the first mode, we can see a rapid rise in the unavailability function by the end of mission time, whereas, in the second mode, the unavailability is stabilized shortly after 500 days, which is close to the mean lifetime of the system (531.7 days). These two curves are limiting curves for optimization. The unavailability growth of the real ageing mode can be reduced using a properly selected FBM; i.e., with a properly selected value for the decision variable n respecting a selected unavailability restriction.
Let us take into account that the maximal permissible value of US(x) is U0 = 0.04. It is In the first mode, we can see a rapid rise in the unavailability function by the end of mission time, whereas, in the second mode, the unavailability is stabilized shortly after 500 days, which is close to the mean lifetime of the system (531.7 days). These two curves are limiting curves for optimization. The unavailability growth of the real ageing mode can be reduced using a properly selected FBM; i.e., with a properly selected value for the decision variable n respecting a selected unavailability restriction.
Let us take into account that the maximal permissible value of U S (x) is U 0 = 0.04. It is necessary to find an optimal system mode for the objective function given by Equation (1); i.e., we look for an optimal value for the decision variable n. To solve this optimization problem, we must provide data related to the maintenance cost: the cost of one FBM action is C R = 12 and the cost of one CM action is C CM = 6, which is in good agreement with practical experience because FBM entails the replacement of the old system with a new one, which is more expensive than repairing it. The results of our optimization process are shown in Table 1, which gives the clear conclusion that the optimal value for the decision variable n = 5 and the corresponding minimal cost based on Formula (9) is 59.97 units. Figure 5 shows the dependence of unavailability on time in the optimal mode respecting the given restriction, as well as its comparison with all of the computed system modes.    The results in Table 1 show that the anticipated gradual growth in unavailability depends on the increasing value of the decision variable n. On the other hand, the cost shows a decreasing trend up to the value of the decision variable n = 5 (except for n = 4) and an increasing trend for the greater values of n. This outcome can easily be explained: the increasing value of the decision variable n results in an increasing number of failures per mission time n R , but the increase in n is faster. For small values of n, it holds true that n R > n, which means that FBM replacements have a direct influence on cost. When n approaches n R , after a certain value of n is reached, the number of FBM replacements is constant and equal to 1, such that the total cost increases successively only at the expense of the repair cost. If n exceeds n R (n > n R ), FBM replacement is not achieved; thus, it does not affect the total cost. For example, the decision variable n = 2 produces n R = 7.7, such that three FBM replacements influence the total cost, whereas for n = 5, . . . , 9, n R ranges from 9.0 to 10.85, which means that all of these modes have only one FBM replacement and the total cost only changes due to repairs.
We would like to remark that the cost optimization resulting from Formula (7) can be significantly influenced by other parameters-for example, C R , C CM -as well as the mission time T M ( [21]), which has a decisive effect on n R . These parameters are fixed in the analysis and the only decision variable is the number of failures until replacement n.

Unavailability Quantification and Discrete Maintenance Optimization of a Four-Component System
The second system selected for optimization using the discrete maintenance model based on FBM is a series-parallel system, which consists of two subsystems connected in series. Each subsystem consists of two components connected in parallel, which is demonstrated in Figure 6. The system has been used many times by different authors to demonstrate an imperfect maintenance model; for example, in [7,32], etc. The system parameters and maintenance characteristics are shown in Table 2. The distribution function of the first failure time X 1 of all components is a Weibull distribution, with the shape parameter β = 2 and scale parameter α = 1500 h for both components of the first subsystem, and α = 2000 h for both components of the second subsystem. Imperfect CM is characterized by random repair times with a rectangular distribution and a CM MTTR = 300 h for both components of the first subsystem and an MTTR = 200 h for both components of the second subsystem. The replacement time in the FBM model is deterministic, with durations of 75 h for both components of the first subsystem and 50 h for both components of the second subsystem. Similarly to the first example, the replacement time for all components is shorter than the CM time since it can be scheduled beforehand. The imperfect CM model is characterized by a real ageing process, as described above. We presume that the real ageing mode of all components is characterized by the same ageing quotient q a = 1.25, so that each component failure (prior to the nth) is followed by CM and the subsequent component lifetime has an increasingly worse failure rate, in accordance with Formulas (10)- (12). The system is investigated with the mission time T M = 8000 h.

Unavailability Quantification and Discrete Maintenance Optimization of a Four-Component System
The second system selected for optimization using the discrete maintenance model based on FBM is a series-parallel system, which consists of two subsystems connected in series. Each subsystem consists of two components connected in parallel, which is demonstrated in Figure 6. The system has been used many times by different authors to demonstrate an imperfect maintenance model; for example, in [7,32], etc. The system parameters and maintenance characteristics are shown in Table 2. The distribution function of the first failure time X1 of all components is a Weibull distribution, with the shape parameter β = 2 and scale parameter α = 1500 h for both components of the first subsystem, and α = 2000 h for both components of the second subsystem. Imperfect CM is characterized by random repair times with a rectangular distribution and a CM MTTR = 300 h for both components of the first subsystem and an MTTR = 200 h for both components of the second subsystem. The replacement time in the FBM model is deterministic, with durations of 75 h for both components of the first subsystem and 50 h for both components of the second subsystem. Similarly to the first example, the replacement time for all components is shorter than the CM time since it can be scheduled beforehand. The imperfect CM model is characterized by a real ageing process, as described above. We presume that the real ageing mode of all components is characterized by the same ageing quotient qa = 1.25, so that each component failure (prior to the nth) is followed by CM and the subsequent component lifetime has an increasingly worse failure rate, in accordance with Formulas (10)- (12). The system is investigated with the mission time TM = 8000 h. Table 2 shows the FBM costs CR and CM cost CCM of all components.    First, the unavailability functions of all components (terminal nodes) were computed according to Formula (6). Then, the system from Figure 6 was transformed into a corresponding AG using the methodology described in [17]. Thereafter, the system unavailability function was calculated using the AG. The difference between real and theoretical ageing is shown in Figure 7, which illustrates the time-dependent unavailability function of the system u(t) with the mission time of 8000 h for two system configurations: 1.
All components (terminal nodes) are in a real ageing mode without FBM; i.e., n = ∞, each failure is followed by CM, and the subsequent lifetime has a correspondingly increased failure rate.

2.
All components are in the traditional ageing mode (theoretical), where ageing is due to an increasing Weibull failure rate and each failure is followed by CM, which brings the system to an as-good-as-new state, such that components are replaced (n = 1).
Mathematics 2022, 10, x FOR PEER REVIEW 13 of 18 2. All components are in the traditional ageing mode (theoretical), where ageing is due to an increasing Weibull failure rate and each failure is followed by CM, which brings the system to an as-good-as-new state, such that components are replaced (n = 1). Not surprisingly, we see a similar effect as was identified for the previous one-component system. In the first system configuration, we could see a rapid rise in the unavailability function, whereas, in the second configuration, the unavailability function peaks at 0.0476 for about 1710 h, which is close to the mean lifetimes of both components of the second subsystem (1772.4 h), and thereafter becomes stabilized. These two curves are, again, limiting curves for optimization. The unavailability growth of the real ageing configuration can be reduced through properly selected FBM for all components. Let us take into account the unavailability restriction that states that the maximal permissible value of US(x) is U0 = 0.08. It is necessary to obtain an optimal system configuration for the objective function given by Equation (1) by applying our discrete maintenance model. We search for optimal component modes; i.e., values of the decision variable n for all four components. In the preliminary calculations, we found that the optimal choice of decision variables ranges from n = 6 to n = 8. As our system is composed of four components and each component has three maintenance modes, we have to investigate 3 4 = 81 maintenance Not surprisingly, we see a similar effect as was identified for the previous onecomponent system. In the first system configuration, we could see a rapid rise in the unavailability function, whereas, in the second configuration, the unavailability function peaks at 0.0476 for about 1710 h, which is close to the mean lifetimes of both components of the second subsystem (1772.4 h), and thereafter becomes stabilized. These two curves are, again, limiting curves for optimization. The unavailability growth of the real ageing configuration can be reduced through properly selected FBM for all components. Let us take into account the unavailability restriction that states that the maximal permissible value of U S (x) is U 0 = 0.08. It is necessary to obtain an optimal system configuration for the objective function given by Equation (1) by applying our discrete maintenance model.
We search for optimal component modes; i.e., values of the decision variable n for all four components. In the preliminary calculations, we found that the optimal choice of decision variables ranges from n = 6 to n = 8. As our system is composed of four components and each component has three maintenance modes, we have to investigate 3 4 = 81 maintenance configurations of the system. The results of this optimization process are shown in Table 3. Table 3. Results for optimization process for decision variable n in a four-component system. The sequence of the index digits in Table 3 corresponds to the numbering of the components (2,2), (2,1), (1,2), and (1,1) and determines the component maintenance modes; i.e., values for the decision variable n that range from 6 to 8 (the corresponding digits range from 1 to 3). For example, the sequence index at row number 4 describes such a maintenance configuration, where component (2,2) has the decision variable n = 8 (denoted by the number 3 in the sequence), whereas all the components (2,1), (1,2), and (1,1) have a maintenance mode with the decision variable n = 6 (all denoted by the number 1 in the sequence).
The values for maximal system unavailability U S (x) in Table 3 are ordered in ascending order so that it can be seen that the last value of U S (x), which is below the restriction U 0 = 0.08, is in row number 46 and the minimal cost of the first 46 calculations computed according to Formula (9) is 125.23, which is also found in row number 46. This represents the optimal system configuration with the following optimal values for the decision variable n: n = 6 for components (2,1) and (2,2) and n = 7 for components (1,1) and (1,2). A comparison of the unavailability function of the optimal FBM and the unavailability limiting curves is shown in Figure 8.  Table 3 are caused by frequent FBM replacements-see often-repeated index 1, which means a minimal value for the decision variable n = 6. T effect is particularly noticeable in the last two columns of the index sequence, which c respond to components (1,1) and (1,2). Frequent FBM replacements produce excessive tal costs that vary between 130-140 units. On the other hand, in the last section of the Ta 3 indexes, 2 and 3 are most common (see again the indexes corresponding to compone (1,1) and (1,2)), which implies system configurations with higher maintenance modes; component modes with higher decision variables n = 7 and 8. These values for the decis variable for components (1,1) and (1,2) exceed the mean number of failures per miss time nR such that FBM replacement is not achieved in most cases, and a relevant par the total cost of these configurations is at the expense of repair costs, which are lower th replacement costs ( Table 2). As a consequence of this, the decreasing total cost mos varies between 126-131 units, whereas the maximal system unavailability US(x) increa to a value of 0.0923.
One more observation can be added. It was mentioned above that components (1 and (1,2) have a higher impact on the decision-making process related to the optimizat problem. One more reason for this is that the mean time to failure of these two com nents is less than the mean time to failure of the remaining two components (comp 1329 versus 1772 h). As the first block in the series-parallel system in Figure 6 has grea unavailability than the second block, it is natural that it is of greater importance to so the optimization problem related to all systems. FBM replacements affect both unavailability and cost. We can further see in Table 3 that the cost of all configurations ranges from 125.2 to 140.6 cost units. Low unavailability values in the first section of Table 3 are caused by frequent FBM replacements-see the often-repeated index 1, which means a minimal value for the decision variable n = 6. This effect is particularly noticeable in the last two columns of the index sequence, which correspond to components (1,1) and (1,2). Frequent FBM replacements produce excessive total costs that vary between 130-140 units. On the other hand, in the last section of the Table 3 indexes, 2 and 3 are most common (see again the indexes corresponding to components (1,1) and (1,2)), which implies system configurations with higher maintenance modes; i.e., component modes with higher decision variables n = 7 and 8. These values for the decision variable for components (1,1) and (1,2) exceed the mean number of failures per mission time n R such that FBM replacement is not achieved in most cases, and a relevant part of the total cost of these configurations is at the expense of repair costs, which are lower than replacement costs ( Table 2). As a consequence of this, the decreasing total cost mostly varies between 126-131 units, whereas the maximal system unavailability U S (x) increases to a value of 0.0923.
One more observation can be added. It was mentioned above that components (1,1) and (1,2) have a higher impact on the decision-making process related to the optimization problem. One more reason for this is that the mean time to failure of these two components is less than the mean time to failure of the remaining two components (compare 1329 versus 1772 h). As the first block in the series-parallel system in Figure 6 has greater unavailability than the second block, it is natural that it is of greater importance to solve the optimization problem related to all systems.

Conclusions
In this paper, we introduced Weibull-based ageing systems that undergo discrete maintenance optimization. The Weibull-based ageing process is considered imperfect CM, where each failure and follow-up repair degrade the system to some extent. After the occurrence of the n-th failure, where n can be determined for each component as an optimal value for the decision variable of the optimization process, FBM is launched and the component is replaced with a new one. The optimization is realized in a context with minimal system costs and a prescribed unavailability restriction. The corresponding reliability mathematics for the unavailability quantification of terminal nodes was derived in this article because this unavailability is used as the main input in system unavailability quantification using AGs as the system representation. Although the renewal process model that originates from the new renewal cycle was only developed in this article, in general, we can conclude that it is a limited case of the model developed in [31] where the length of the inspection period approaches zero. A cost model respecting the imperfect CM process with FBM was further introduced because it is indispensable for discrete maintenance optimization.
Numerical experiments showed that the discrete maintenance optimization method is a viable method to make an optimal decision for FBM-replacement of system components that undergo real ageing. Although the computing process may in some cases be heavy on computing time (see the 81 system configurations of the four component systems), none of the computations in this paper exceeded 4 min. All these computing experiments, including both the development of the algorithm for the unavailability function u(t) and the discrete maintenance optimization, were numerically realized with the high-performance language MATLAB on computing equipment with the following characteristics: Intel (R) Core™ i7-3770 CPU @ 3.4 GHz and 3.9 GHz, 8.00 GB RAM.
Our future research work will be based on our original achievements resulting from our cooperation with power industry experts ( [17]). It will involve maintenance optimization of power networks with a particular focus on circuit breakers.

Conflicts of Interest:
The authors declare no conflict of interest.

CM
corrective maintenance PM preventive maintenance FBM failure-based preventive maintenance AG directed acyclic graph C S total maintenance cost of a system configuration f (x) = min C S objective function U(x,t) instantaneous time-dependent unavailability function U S (x) maximal system unavailability within a mission time T M U 0 a specified limitation of U S (maximal permissible value) x = (x 1 , . . . , x k ) ∈ R k decision variable n i number of failures of the i-th component, which starts its renewal k number of system components X 1 the time from the beginning of the renewal cycle to the first failure X i the time from the occurrence of the i-th failure to the end of its CM repair (for i ∈ {1, . . . , n − 1}) Y i the time from the occurrence of the i-th failure to the end of its CM repair (for i ∈ {1, . . . , n − 1}) Y n the time from the occurrence of the n-th failure to the end of the FBM replacement X * i = X 1 + Y 1 + · · · + the time from the beginning of the renewal cycle to the occurrence X (i−1) + Y (i−1) + X i of the i-th failure, X * 1 = X 1 , X * 2 = X 1 + Y 1 + X 2 , etc. F X (t) = P(X ≤ t) the cumulative distribution function of the random variable X F X (t) = 1 − F X (t) the reliability function of the random variable X