Multi-Objective Optimum Design and Maintenance of Safety Systems: An In-Depth Comparison Study Including Encoding and Scheduling Aspects with NSGA-II

: Maximising proﬁt is an important target for industries in a competitive world and it is possible to achieve this by improving the system availability. Engineers have employed many techniques to improve systems availability, such as adding redundant devices or scheduling maintenance strategies. However, the idea of using such techniques simultaneously has not received enough attention. The authors of the present paper recently studied the simultaneous optimisation of system design and maintenance strategy in order to achieve both maximum availability and minimum cost: the Non-dominated Sorting Genetic Algorithm II (NSGA-II) was coupled with Discrete Event Simulation in a real encoding environment in order to achieve a set of non-dominated solutions. In this work, that study is extended and a thorough exploration using the above-mentioned Multi-objective Evolutionary Algorithm is developed using an industrial case study, paying attention to the possible impact on solutions as a result of different encodings, parameter conﬁgurations and chromosome lengths, which affect the accuracy levels when scheduling preventive maintenance. Non-signiﬁcant differences were observed in the experimental results, which raises interesting conclusions regarding ﬂexibility in the preventive maintenance strategy.


Introduction
System Reliability (R(t)) can be defined as the probability of failure free operation under specified conditions over an intended period of time [1]. System Maintainability (M(t)) can be defined as the probability of being restored to a fully operational condition within a specific period of time [2]. These definitions lead to interest both in time taken for a system to failure (Time To Failure) and in time taken to repair the system (Time To Repair). System Availability (A(t)) can be defined as the fraction of the total time in which systems are able to perform their required function [2]. The concepts of Reliability and Maintainability are related to Availability in order to define the way in which the system is able to achieve the function for which it was designed, over a period of time. Therefore, whereas Reliability is a concept related to non-repairable systems, Availability is a concept related to repairable systems because it encompasses the complete failure-recovery cycle over the mission time. From above, it can be seen that repairable systems operate during a certain period of time (Time To Failure) until a failure occurs. After that, a period of time is needed to recover the system operating status (Time To Repair). This creates an interest in Time To Failure and Time To Repair, which can be modelled as random variables that can be represented by continuous probability distributions.
Optimisation through Evolutionary Algorithms is useful when complex problems have to be solved; that is, problems in which the number of potential solutions is usually high and whilst achieving the best solution is a daunting task and achieving the best solution can be extremely difficult, achieving at least a good solution (if not exactly the best) can be a manageable task [3]. Typical reliability optimisation problems have considered the possibility of maximising the system Availability. There are many techniques commonly used to improve it and the present paper focuses on two of them. On the one hand, modifying the structural design through redundancies improves the system Availability. A redundancy is a component added to a subsystem from a series-parallel configuration in order to increase the number of alternative paths [4]. On the other hand, an overall improvement of system availability is possible through preventive maintenance [5]. The unavailability of a system due to its failure can occur at any time, which necessitates a significant effort to recover the operating state. Conversely, a programmed shutdown to perform a preventive maintenance task represents a controlled situation with materials, spares and human teams available.
These techniques were jointly explored by the authors of the present paper preliminarily in [6], coupling Multi-objective EvolutionaryAlgorithm and Discrete Simulation and dealing with the joint optimisation of systems design (considering whether or not to include redundant devices) and their preventive maintenance strategy (choosing optimum periodical preventive maintenance times in relation to each system device). This allowed necessary preventive maintenance activities to be carried out. The system Availability and operation Cost were the objectives to maximise and minimise, respectively. A simulation approach was used in which each solution proposed by the Multi-objective Evolutionary Algorithm (using real encoding) was evaluated through Discrete Simulation; the technique used to build and modify (depending on the periodical preventive maintenance times) the Functionability Profile. This is a concept presented by Knezevic [7], which explains the history of the system varying among operating and recovery times. It is a powerful modelling technique, which permits analysis of complex systems with a reality more realistic representation of their behaviour. In the previous study, the performance of several configurations of the Multi-objective Evolutionary Algorithm Non-dominated Sorting Genetic Algorithm II (NSGA-II) [8] in a real encoding environment was explored. In the present paper an in-depth encoding comparative study is developed. Firstly, the real encoding case study is extended, and secondly, some binary encoding alternatives are explored, looking for possible advantages and disadvantages to encode this kind of real world problems. Moreover, an accuracy-level experiment for the preventive maintenance strategy is considered. The first part of the study determines the optimum periodical Time To Start a Preventive Maintenance activity measured in hours. Two more levels are explored in this study: days and weeks. There are preventive maintenance activities whose accuracy level in time can be of little importance. It may not be important to determine the exact instant for their development, being enough to define the day or the week. Therefore, the effect of several chromosome lengths is explored looking to improve the evolutionary process. Summarising the contributions of the present study: • In this work, seven encoding alternatives are thoroughly explored: Real encoding (with Simulated Binary Crossover), Binary encoding (with 1 point Crossover), Binary encoding (with 2 point Crossover), Binary encoding (with Uniform Crossover), Gray encoding (with 1 point Crossover), Gray encoding (with 2 point Crossover) and Gray encoding (with Uniform Crossover). Their performances are compared using the Hypervolume indicator and statistical significance tests. • Additionally, three accuracy levels on time are explored for the binary encoding; hours, days and weeks, in order to analyse the effect of chromosome length in the evolutionary search and final non-dominated set of solutions. Their performances are compared using the Hypervolume indicator and statistical significance tests. • The methodology is applied in an industrial test case, obtaining an improved non-dominated front of optimum solutions that could be considered as both a benchmark case and reference solution.
The paper is organised as follows. Section 2 explores the related literature. Section 3 shows an outline of the methodology. Section 4 presents the case study and Section 5 the description of experiments carried out (encodings and accuracy levels). In Section 6, results are shown and discussed, and finally, Section 7 introduces the conclusions.

Literature Review
Approaches dealing with reliability systems design (redundancy allocation) are presented first, and approaches dealing with preventive maintenance optimum design are presented secondly in this section. Finally, the third subsection deals with the simultaneous optimisation of reliability systems design and preventive maintenance.

Redundancy Allocation of Reliability Systems Design Optimisation
Improving the Reliability or Availability for series-parallel systems through redundancy allocation using Multi-objective Evolutionary Algorithms has been considered in several studies. Many authors have developed Genetic Algorithms-based approaches to solve the problem. Bussaca et al. [9] utilized a Multi-objective Genetic Algorithm to optimise the design of a safety system through redundancy allocation, considering components with constant failure rates. They employed profit during the mission time (formed by the profit from plant operation minus costs due to purchases and installations, repairs and penalties during downtime) and the reliability at mission time as objective functions. Marseguerra et al. [10] presented an approach that couples a Multi-objective Evolutionary Algorithm and Monte Carlo simulation for optimal networks design aimed at maximising the network reliability estimate and minimising its uncertainty. Tian and Zuo [11] proposed a multi-objective optimisation model for redundancy allocation for multi-state series-parallel systems. They used physical programming as an approach to optimise the system structure and a Genetic Algorithm to solve it. The objectives consisted of maximising the system performance and minimising its cost and weight. Huang et al. [12] proposed a niched Pareto Genetic Algorithm to solve reliability-based design problems aiming to achieve a high number of feasible solutions. They used reliability and cost as objective functions. Zoulfaghari et al. [13] presented a Mixed Integer Nonlinear Programming model to availability optimisation of a system taking into account both repairable and non-repairable components. They developed a Genetic Algorithm looking for maximum availability and minimum cost. Taboada et al. [14] introduced a Genetic Algorithm to solve multi-state optimisation design problems. The universal moment generating function was used to evaluate the reliability indices of the system. They used reliability, cost and weight as objective functions.
The use of the NSGA and NSGA-II algorithms has been extensive; some applications ranging between the years 2003-2021 are shown as follows. Taboaba et al. [15] presented two methods to reduce the size of the Pareto optimal set for multi-objective system-reliability design problems: on the one hand using a pseudo ranking scheme and on the other hand using data mining clustering techniques. To demonstrate the performance of the methods, they solved the redundancy allocation problem using the Non-dominated Sorting Genetic Algorithm (NSGA) to find the Pareto optimal solutions, and then the methods were successfully applied to reduce the Pareto set. They studied reliability, cost and weight as objective functions. However, the most widely used standard Genetic Algorithm is the second version of the NSGA Multi-objective Evolutionary Algorithm. Greiner et al. [16] introduced new safety systems multi-objective optimum design methodology (based on fault trees evaluated by the weight method) using not only the standard NSGA-II but also the Strength Pareto Evolutionary Genetic Algorithm (SPEA2) and the controlled elitist-NSGA-II, with minimum unavailability and cost criteria. Salazar et al. [17] developed a formulation to solve several optimal system design problems. They used the NSGA-II to achieve maximum reliability and minimum cost. Limbourg and Kochs [18] applied a specification method originating from software engineering named Feature Modelling and a NSGA-II with an external repository. They maximised the life distribution and minimised the system cost. Kumar et al. [19] proposed a multi-objective formulation of multi-level redundancy allocation optimisation problems and a methodology to solve them. They proposed a hierarchical Genetic Algorithm framework by introducing hierarchical genotype encoding for design variables. Not only the NSGA-II but also the SPEA2 Multi-objective Evolutionary Algorithm were used, considering reliability and cost as objective functions. Chambari et al. [20] modelled the redundancy allocation problem taking into account non-repairable series-parallel systems, non-constant failure rates and imperfect switching of redundant cold-standby components. They used NSGA-II as well as Multi-objective Particle Swarm Optimisation (MOPSO) to solve the problem with reliability and cost as objective functions. Safari [21] proposes a variant of the NSGA-II method to solve a novel formulation for redundancy allocation. He considered for the redundancy strategy both active and cold-standby redundancies with reliability and cost as objective functions. Ardakan et al. [22] solved the redundancy allocation by using mixed redundancies combining both active and cold-standby redundancies. Under this approach, the redundancy strategy is not predetermined. They studied reliability and cost as objective functions using the NSGA-II method. Ghorabaee et al. [23] considered reliability and cost as objective functions to solve the redundancies allocation problem using NSGA-II. They introduced modified methods of diversity preservation and constraints handling. Amiri and Khajeh [24] considered repairable components to solve the redundancy allocation problem in a series-parallel system. They used the NSGA-II method with availability and cost as objective functions. Jahromi and Feizabadi [25] presented a formulation for the redundancy allocation of non-homogeneous components considering reliability and cost as objective functions. The NSGA-II method was used to achieve the Pareto optimal front. Kayedpour et al. [26] developed an integrated algorithm to solve reliability design problems taking into account instantaneous availability, repairable components and a selection of configuration strategies (parallel, cold or warm) based on Markov processes and the NSGA-II method. As objective functions, they considered availability and cost. Sharifi et al. [27] presented a new multi-objective redundancy allocation problem for systems where the subsystems were considered weighted-k-out-of-n parallel. They used NSGA-II with reliability and cost as objective functions. Chambari et al. [28] proposed a bi-objective simulation-based optimisation model to redundancy allocation in series-parallel systems with homogeneous components to maximise the system reliability and minimise the cost and using NSGA-II. They considered optimal component types, the redundancy level, and the redundancy strategy (active, cold standby, mixed or K-mixed) with imperfect switching.
Other Multi-objective evolutionary and bio-inspired methods have been used in a lesser measure. Zhao et al. [29] optimised the design of series-parallel systems using the Ant Colony Algorithm in a multi-objective approach, considering reliability, cost and weight as objective functions. Chiang and Chen [30] proposed a Multi-objective Evolutionary Algorithm based on simulated annealing to solve the availability allocation and optimisation problem of a repairable series-parallel system. They applied their algorithm to two study cases presented in references [31] (with objective functions availability and cost) and [9]. Khalili-Damghani et al. [32] proposed a novel dynamic self-adaptive Multi-objective Particle Swarm Optimisation method to solve two-states reliability redundancy allocation problems. They contemplated reliability, cost and weight as objective functions. Jiansheng et al. [33] proposed an Artificial Bee Colony Algorithm to solve the redundancy allocation problem for repairable series-parallel systems where uncertainty in failure rates, repair rates and other relative coefficients involved were considered. They used availability and cost as objective functions. Samanta and Basu [34] proposed an Attraction-based Particle Swarm Optimisation to solve availability allocation problems for systems with repairable components, where they introduced fuzzy theory to manage uncertainties. They used availability and cost as objective functions.

Preventive Maintenance Strategy Optimisation
Improving the Reliability or Availability for series-parallel systems through the preventive maintenance strategy has been widely studied using Multi-objective Evolutionary Algorithms, too. Many authors used Genetic Algorithms as an optimisation method. Muñoz et al. [35] presented an approach based on Genetic Algorithms and focused on the global and constrained optimisation of surveillance and maintenance of components based on risk and cost criteria. Marseguerra et al. [36] coupled Genetic Algorithms and Monte Carlo simulation in order to optimise profit and availability. The Monte Carlo simulation was used to model system degradation while the Genetic Algorithm was used to determine the optimal degradation level beyond which preventive maintenance has to be performed. Gao et al. [37] studied the flexible job shop scheduling problem with availability constraints affecting maintenance tasks. They used a Genetic Algorithm looking for minimum makespan (time that elapses from the start of work to the end), time expended on a machine and the total time expended on all machines. Sánchez et al. [38] used Genetic Algorithms for the optimisation of testing and maintenance tasks with unavailability and cost as objective functions. They considered the epistemic uncertainty in relation to imperfect repairs. Wang and Pham [39] used a Genetic Algorithm to estimate the preventive maintenance interval allowing for imperfect repairs and the number of preventive maintenance activities before component needs to be replaced. They used availability and cost as objective functions. Ben Ali et al. [40] developed an elitist Genetic Algorithm to deal with the production and maintenance-scheduling problem, minimising makespan and cost. Gao et al. [5] studied preventive maintenance considering the dynamic interval for multi-component systems. They solved the problem using Genetic Algorithms with availability and cost as objective functions. An et al. [41] built a novel integrated optimisation model including the flexible job-shop scheduling problem to reduce energy consumption in the manufacturing sector. They considered degradation effects and imperfect maintenance. They proposed a Hybrid Multi-objective Evolutionary Algorithm taking into account the makespan, total tardiness, total production cost and total energy consumption as objective functions. Bressi et al. [42] proposed a methodology to minimise the present value of the life cycle maintenance costs and maximise the life cycle quality level of the railway track-bed considering different levels of reliability. They used a Genetic Algorithm to achieve optimal solutions.
The use of the NSGA-II algorithm and other non-dominated criterion-based approaches has been extensive with applications ranging across the years 2005-2021. Martorell et al. [43] proposed a methodology to take decisions and determine technical specifications and maintenance looking for increasing reliability, availability and maintainability. They used SPEA2 as an optimisation method. Oyarbide-Zubillaga et al. [44] coupled Discrete Event Simulation and Genetic Algorithms (NSGA-II) to determine the optimal frequency for preventive maintenance of systems under cost and profit criteria. Berrichi et al. [45] proposed a new method to deal with simultaneous production and maintenance scheduling. They used the Weighted Sum Genetic Algorithm (WSGA) and NSGA-II as optimisation methods to compare their performances. They worked with makespan and unavailability due to maintenance tasks as objective functions. Moradi et al. [46] studied simultaneous production and preventive maintenance scheduling to minimise the global time invested in production tasks and unavailability due to preventive maintenance activities. They used four Genetic Algorithms: NSGA-II, NRGA (Non-ranking Genetic Algorithm), CDRNSGA-II (NSGA-II with Composite Dispatching Rule and active scheduling) and CDRNRGA (NRGA with Composite Dispatching Rule and active scheduling). Hnaien and Yalaoui [47] considered a similar problem, minimising the makespan and the delay between the real and the theoretical maintenance frequency for two machines. They used NSGA-II and SPEA2, including two new versions based on the Johnson Algorithm. Wang and Liu [48] considered the optimisation of parallel-machine-scheduling integrated with two kinds of resources (machines and moulds) and preventive maintenance planning. They used makespan and availability as objective functions and NSGA-II as an optimisation method. Piasson et al. [49] proposed a model to solve the problem of optimising the reliability-centred maintenance planning of an electric power distribution system. They used NSGA-II to achieve the Pareto optimal front and, as objective functions, they minimised the cost due to maintenance activities and maximised the index of reliability of the whole system. Sheikhalishahi et al. [50] presented an open shop scheduling model that considers human error and preventive maintenance. They considered three objective functions: human error, maintenance and production factors. They used NSGA-II and SPEA2 as optimisation methods. As well as that, they used another Evolutionary Algorithm, the Multi-Objective Particle Swarm Optimisation (MOPSO) method. Boufellouh and Belkaid [51] proposed a bi-objective model that determines production scheduling, maintenance planning and resource supply rate decisions in order to minimise the makespan and total production costs (including total maintenance, resource consumption and resource inventory costs). They used NSGA-II and BOPSO (Bi-Objective Particle Swarm Optimization) as Evolutionary Optimisation Algorithms. Zang and Yang [52] proposed a multi-objective model of maintenance planning and resource allocation for wind farms using NSGA-II. They considered the implementation of maintenance tasks using the minimal total resources and at the minimal penalty cost.
Other Multi-objective Evolutionary methods have been used. Berrichi et al. [53] solved the joint production and preventive maintenance-scheduling problem using the Ant Colony Algorithm with availability and cost as objective functions. Suresh and Kumarappan [54] presented a model for the maintenance scheduling of generators using hybrid Improved Binary Particle Swarm Optimisation (IBPSO). As objective functions, they used a reduction in the loss of load probability and minimisation of the annual supply reserve ratio deviation for a power system. Li et al. [55] presented a novel Discrete Artificial Bee Colony Algorithm for the flexible job-shop scheduling problem considering maintenance activities. They used as objective functions the makespan, the total workload of machines and the workload of the critical machine.

Redundancy Allocation and Preventive Strategy Optimisation
Therefore, it is possible to improve the availability of repairable systems by dealing with their design and employing a preventive maintenance strategy. However, only a few works have been developed to look at the simultaneous optimisation of both from a multi-objective point of view. In Galván et al. [56], a methodology for Integrated Safety System Design and Maintenance Optimisation based on a bi-level evolutionary process was presented. While the inner loop is devoted to searching for the optimum maintenance strategy for a given design, the outer loop searches for the optimum system design. They used Genetic Algorithms as an optimisation method and cost and unavailability as objective functions. Okasha and Frangopol [57] considered the simultaneous optimisation of design and maintenance during the life cycle using Genetic Algorithms. They studied system reliability, redundancy and life-cycle cost as objective functions. Adjoul et al. [58] described a new approach to simultaneous optimisation of design and maintenance of multi-component industrial systems improving their performances with reliability and cost as objective functions. They used a two level Genetic Algorithm; the first optimises the design based on reliability and cost, and the second optimises the dynamic maintenance plan.
This work studies the simultaneous optimisation of design and preventive maintenance strategy coupling Multi-objective Evolutionary Algorithms and Discrete Simulation: a technique that has achieved good results in the Reliability field. Coupling Multi-objective Evolutionary Algorithms with Discrete Simulation has been studied, on the one hand, to supply redundancy allocation [10] and, on the other hand, to determine the preventive maintenance strategy [36,44]. Moreover, only a few works have been developed looking at design and corrective maintenance strategy simultaneously [59,60]. Nevertheless, to the knowledge of the authors of this work, coupling Multi-objective Evolutionary Algorithms and Discrete Simulation has not yet been explored for both joint optimisation of the design and preventive maintenance strategy as in the current study where, additionally, a variety of encoding schemes were considered in the analysis. In the studies cited in the present literature review, some used real encoding, while others used binary or integer encoding; however, not one of them studied the possible impact of such encoding schemes on performance. Moreover, an accuracy experiment is developed in the present paper, including the time unit needed to carry out the preventive maintenance activities. As shown in the above literature review, NSGA-II is one of the state-of-the-art Multi-objective Evolutionary Algorithms more commonly used to deal with redundancy allocation and scheduling preventive maintenance problems. Therefore, this method is thorough explored in the present paper.

Extracting Availability and Cost from Functionability Profiles
Availability can be computed by using the unconditional failure w(t) and repair v(t) intensities, as is explained in Ref. [2]. These are described by the following Equation (1), where f (t) is the failure density function of a system, t 0 f (t − u)v(u)du is the failure probability of the cited system in interval [0, t) having worked continuously since being repaired in [u, u + du) given that it was working at t = 0 and t 0 g(t − u)w(u)du is the repair probability in interval [0, t), given that it has been in the failed state since the previous failure in [u, u + du) and that it was working at t = 0. (1) When the system devices have exponential failure and repair intensities (constant failure and repair rates), it is not relatively straightforward to find its availability using the solutions of Equation (1). However, when the system devices do not have exponential failure and/or repair intensities, finding the system availability using Equation (1) is difficult, so a simulation approach may be more suitable. In this paper, the system availability is characterised in a simulation approach by using the system Functionability Profile, a concept introduced by Knezevic [7] and defined as the inherent capacity of systems to achieve the required function under specific features when they are used as specified. Speaking in general, all systems achieve their function at the beginning of their lives. However, irreversible changes occur over time and variations in system behaviour happen. The deviation of the variations in relation to the satisfactory features reveals the occurrence of system failure, which causes the transition from operating state to failure state. After failing, recovery activities (corrective maintenance) can recover its capacity to fulfil the required function when the system is repairable.
Additional tasks to maintain the operating status could be carried out. These are called preventive maintenance activities. These are generally less complex than a repair and should be done prior to failure. From the Functionability Profile point of view, the states of a repairable system fluctuate between operation and failure over the mission time. The shape of cited changes is called the Functionability Profile because it shows the states over the mission time. Therefore, Functionability Profiles depend on operation times (either Time To Failure or Time To Start a scheduled Preventive Maintenance activity) (t f 1 , t f 2 , ..., t f n ) and recovery times (either Time To Repair after failure or Time To Perform a scheduled Preventive Maintenance activity) (t r1 , t r2 , ..., t rn ). It is obvious that, after a period of operation, a period of recovery is needed.
In the present paper, the system Functionability Profile is built by using Discrete Event Simulation in a simulation approach as will be explained later. Once known, the system Functionability Profile will be able to compute the system Availability through the relation between the system operation times and the system recovery times. The system will be able to fulfil its purpose during t f times, so it is possible to evaluate its Availability at mission time by using Equation (2). (2) where: • n is the total number of operation times, Operation and recovery times are random variables so they may be treated statistically. They can be defined as probability density functions through their respective parameters. There are several databases such as OREDA [61], which supply the characteristic parameters for the referred functions, so operation and recovery times can be characterised for system devices. When systems are operating, earnings are generated due to the availability of the system. Conversely, when systems have to be recovered, economic cost is invested in order to regain the operation status. In this paper, the economic cost is a variable directly related to recovery times, which are related to corrective and preventive maintenance activities; quantities computed by Equation (3). where: • C is the system operation cost quantified in economic units, • q is the total number of corrective maintenance activities, • cc i is the cost due to the i-th corrective maintenance activity, • p is the total number of preventive maintenance activities, • cp j is the cost due to the j-th preventive maintenance activity.
In this work, costs derived from maintenance activities depend on fixed quantities per hour (corrective and preventive) so the global cost is directly related to the recovery times. Preventive maintenance activities are scheduled shutdowns so recovery times will be shorter and more economical than recovery times due to corrective maintenance activities (for reasons explained before, such as access to human personnel who are willing and/or trained, and/or the availability of spare parts). It will be necessary to carry out preventive maintenance activities to avoid long recovery times. These should ideally be done before the failure but as close as possible to it. Therefore, the basic Functionability Profile of the system (i.e., the system Functionability Profile, which represents the continuous cycle of failure-repair within the mission time, without considering preventive maintenance activities) should be modified by including preventive maintenance activities for system devices. This approach makes it possible to maximise the system Availability and minimise the costs due to recovery times.

Building Functionability Profiles to Evaluate the Objective Functions
It is necessary to characterise both the system Availability and the Cost from the system Functionability Profile in order to optimise the system design and preventive maintenance strategy. The Functionability Profiles of all system devices are built by using Discrete Event Simulation. Finally, the system Functionability Profile is built from these Functionability Profiles. With this purpose, it is necessary to have information about how to characterise operation Times To Failure (TF) and Times To Repair after failure (TR), which are related to the parameters of their probability density functions. The Functionability Profiles for system devices are built by generating random times, which are obtained from the respective probability density functions, both for Times To Failure (TF) and Times To Repair (TR). In order to modify the Functionability Profiles relating to the preventive maintenance activities, Times To Start a Preventive Maintenance (TP) have to be used. This information is supplied via each solution provided by the behaviour Algorithm (each individual of the population) which is used to build the Functionability Profile of each device through Discrete Simulation. Moreover, recovery times in relation to the preventive maintenance activities or Times To Perform a Preventive Maintenance activity (TRP) have to be introduced by generating random times within the limits previously fixed. The process is explained below: 1. System mission time must be defined and then, the process continues for all devices. 2. The device Functionability Profile (PF) must be initialised.

The Time To Start a Preventive Maintenance activity (TP) proposed by the Multi-objective
Evolutionary Algorithm is extracted from the individual (candidate solution) that is being evaluated and a Time To Perform a Preventive Maintenance activity (TRP) is randomly generated, within the limits previously fixed. 4. With reference to the failure probability density function related to the device, a Time To Failure (TF) is randomly generated, within the limits previously fixed. 5. If TP < TF, the preventive maintenance activity is performed before the failure. In this case, as many operating times units as TP units followed by as many recovery times units as TRP units are added to the device Functionability Profile. Each time unit represented in this way (both operating times and recovery times) is equivalent to one hour, day or week of real time, depending on the accuracy level chosen. 6. If TP > TF, the failure occurs before carrying out the preventive maintenance activity.
In this case, taking into consideration the repair probability density function related to the device, the Time To Repair after the failure (TR) is randomly generated, within the limits previously fixed. Then, as many operating times units as TF units followed by as many recovery times units as TR units are added to the device Functionability Profile. Each time unit represented in this way (both operating times and recovery times) is equivalent to one hour, day or week of real time, depending on the accuracy level chosen. 7. Steps 4 to 6 have to be repeated until the end of the device mission time. 8. Steps 2 to 7 have to be repeated until the construction of the Functionability Profiles of all devices. 9. After building the Functionability Profiles of the devices, the system Functionability Profile is built by referring to the series (AND) or parallel (OR) distribution of the system devices.
Once the system Functionability Profile is built, the values of the objective functions can be computed by using both Equation (2) (evaluating the system Availability in relation to the time in which the system is operating and being recovered) and Equation (3) (evaluating the system operation Cost depending on the cost of the time units in relation to the development of corrective or preventive maintenance).

Multi-Objective Optimisation Approach
The optimisation method used, the Non-dominated Sorting Genetic Algorithm II (NSGA-II) belongs to the Evolutionary Algorithms (EA) paradigm. These kind of methods use a population of individuals with a specific size. Each individual is a multidimensional vector, called a chromosome, representing a possible candidate solution to the problem, while the vector components are called genes or decision variables. NSGA-II uses the concept of Pareto dominance as the basis of its selection criterion. Extended information on Evolutionary Optimisation Algorithms can be found in Ref. [3] and related to Multi-objective Evolutionary Algorithms in Ref. [62]. In this work, each individual in the population consists of a string (of real numbers or integers, depending on the encoding used, with values between 0 and 1) in which the system design alternatives and the periodic Times To Start a Preventive Maintenance activity related to each device included in the system design is codified. Optimisation problems can be minimised or maximised for one or more of the objectives. In most cases, real world problems present various objectives in conflict with each other. These problems are called "Multi-objective" and their solutions arise from a solutions set which represents the best compromise among the objectives (Pareto optimal set) [62,63]. These kind of problems are described by Equation (4) (considering a minimisation problem) [3].
In problems defined in this way, the k functions have to be minimised simultaneously. In the present paper, the objective functions are, on the one hand, the system Availability (first objective function, which is maximised and which is mathematically expressed by Equation (2); maximising Availability is similar to minimise Unavailability) and, on the other hand, the operation Cost (second objective function, which is minimised and which is mathematically expressed by Equation (3)).

The Case Study
The case study is based on the case presented in Ref. [6]. It consists of simultaneously optimising the design and the preventive maintenance strategy for an industrial fluid injection system, which is composed of cut valves (V i ) and impulsion pumps (P i ), taking Availability and operation Cost as objective functions. The desired outcome is maximum Availability and minimum operation Cost. The higher the investment in preventive maintenance, the better the system Availability. Conversely, this policy implies the growth of unwanted Cost. The system scheme is shown in Figure 1. The data used in this work are shown in Table 1. Extended information regarding the parameters is supplied in Appendix A. The data were obtained from specific literature [61], expert judgement (based on professional experience from the Machinery & Reliability Institute (MRI), Alabama, USA) and mathematical relations. In this sense, the TR σ for valves and pumps were set in relation to the µ of their respective normal distribution functions and their TR min previously established. In relation to the TR max , it is known that 99.7% of the values of a normally distributed variable are included in the interval µ ± 3σ. The interval was extended to µ ± 4σ, taking into account anecdotal further values. As shown above, optimisation objectives consist of maximising the system Availability and minimising the operation Cost due unproductive system phases (both because the system is being repaired and because the system is being maintained). To do that: • It is necessary to establish the optimum period to perform a preventive maintenance activity for the system devices, and • It is necessary to decide whether to include redundant devices such as P2 and/or V4 by evaluating design alternatives. Including redundant devices will improve the system Availability but it will also increase the system operation Cost.

Description of the Experiments Carried Out
Two sets of experiment comparisons have been developed: first, comparing several encodings (real, binary and gray), and second, comparing several accuracies in the binary encoding. Finally, a description of the NSGA-II configuration is shown in the last subsection.

Comparing Encodings
From the optimisation point of view, it was explained before that the Evolutionary Algorithm (EA) uses a population of individuals called chromosomes, which represent possible solutions to the problem through their decision variables. The encoding of the system parameters is a crucial aspect of the algorithm. This has a significant influence on whether or not the algorithm works [3]. In the present paper, we intend to check whether there is a significant difference between the performances of the different encodings. Depending on the encoding type, each individual is codified as follows: • Real encoding: This is formed by strings of real numbers (with 0 as a minimum value and 1 as a maximum value) following the shape [ The presence of redundant devices, P2 and V4, is defined by B 1 and B 2 , respectively, and the optimum Time To Start a Preventive Maintenance activity in relation to each system device is represented by T 1 to T 7 . The values of the decision variables must be scaled and rounded, i.e., B 1 and B 2 are rounded to the nearest integer (0 implies that the respective device is not selected whereas 1 implies the opposite). T 1 to T 7 are scaled among the respective TP min and TP max (depending on the type of device) and rounded to the nearest integer using Equation (5) • Binary encoding: This is formed by strings of binary numbers that vary between 0 and 1, where the total number of bits is 103 and they are:  (7), where B represents the decimal value of the binary string T 18 to T 30 (e.g., if the values of the decision variables in binary encoding are 1 0 1 1 0 1 1 0 0 0 1 1 1, the value on a scale of 8192 steps will be 5831. If 5840 steps are scaled, the number achieved is 5831 × 0.712890625 ≈ 4157 steps. Therefore, the true Time To Start a Preventive Maintenance activity amounts to 4157 + 2920 = 7077 h).
5. T 89 to T 103 : These denote the Time To Start a Preventive Maintenance activity to the valve V7. The behaviour of its encoding is similar to the behaviour explained for the valve V1.
• Gray encoding: When binary encoding is used, close numbers could bring big scheme modifications (e.g., the code for 15 is 0 1 1 1 while the code for 16 is 1 0 0 0, which represents changes in four (all) bits). Conversely, very similar numbers can represent numbers that are very apart (e.g., the code for 0 is 0 0 0 0 while the code for 16 is 1 0 0 0). In addition to standard binary encoding, Gray encoding is used. Gray code is a binary code where the difference among neighbouring numbers always differs by one bit [64][65][66].
Therefore, the performance of real, standard binary and Gray encodings can be compared.

Comparing Accuracies
Once the encoding experiment had been developed, a second experiment was executed. This consisted of studying the possible impact of the size of the chromosome. In the encoding experiment, the hour was utilised by the chromosomes as a measure of time. In this case, based on the idea that the exact hour to develop a preventive maintenance task is not necessary, the day and the week were used as measures of time. Therefore, in these cases, the solution regarding preventive maintenance strategy consisted of supplying the Time To Start a Preventive Maintenance activity with the day or the week as a time unit, respectively. The consequence was a reduction in the size of the chromosome, which was applied to the binary encoding as follows: • Binary encoding-Days: This is formed by strings of binary numbers that vary between 0 and 1, where the total number of bits is 73 and they are:   • Binary encoding-Weeks: It is formed by strings of binary numbers that vary between 0 and 1, where the total number of bits is 54 and they are: 1. B 1 : This denotes the presence of the pump P2 in the system design (0 implies that the device is not selected whereas 1 implies the opposite).

NSGA-II Configuration
The parameters used to configure the NSGA-II evolutionary process are shown in Table 2. Depending on the encoding applied, specific parameters have to be set, which are described below: Each configuration was executed 10 times (for statistical purposes) with 10,000,000 evaluations used as the stopping criterion. Scale factors in relation to the value of the objective functions were used with the purpose of achieving a dispersed non-dominated front with the unit as maximum value. The values were obtained through a practical approach in which the values of the scale factors are extracted from the values of the objective functions when the optimisation process starts. This approach is based on the assumption that the values of the objective functions will improve over the evolutionary process. The scale factors were used as follows: • The scale factor used to compute the Cost was 1700 economic units. • The scale factor used to compute the system Unavailability was 0.003.
Finally, a two dimensional reference point is needed to compute the Hypervolume indicator. The cited point has to cover the values limited by the scale factors, which restricts the values of the objective functions to a maximum of one. The reference point is set to (2,2). The Software Platform PlatEMO [67] (programmed in MATLAB) is used to optimise the application case. The Design and Maintenance Strategy analysis software (as explained in Section 3) has been developed and implemented into the platform to solve the case study shown in Section 4.

Results
Due to the complexity of the problem, a general-purpose calculation cluster was used during the computational process. The cluster is composed of 28 calculation nodes and one access node. Each calculation node consists of two Intel Xeon E5645 Westmere-EP processors with 6 cores each and 48 GB of RAM memory, allowing 336 executions to be run simultaneously.
Once the results were obtained, valuable information emerged. For each analysed case, the following information is provided: Firstly, information related to the computational process is given with the purpose of showing the complexity of the problem and computational cost. It consists of the time taken for 10 executions of the nine configurations (three population sizes and three mutation rates) related to each analysed case. Secondly, the values of the main measures obtained for the final evaluation are shown. These measures are the Average, Median, Minimum Value, Maximum Value and Standard Deviation of the Hypervolume indicator [68] (HV). Thirdly, in order to establish the existence of significant differences among the performance of the analysed case, a rigorous statistical analysis is carried out. The Friedman's test allows significant differences among results obtained to be detected, and the null hypothesis (H 0 ) to be rejected in that case. The p-value is a useful datum, which represents the smallest significant value that can result in the rejection of H 0 . The p-value provides information about whether a statistical hypothesis test is significant (or not), and also indicates how significant the result is: The smaller the p-value (<0.05), the stronger the evidence against the null hypothesis. Finally, the Hypervolume is computed for the accumulated best non-dominated solutions obtained (the non-dominated front). These represent the best equilibrium solutions among the objectives and the computational procedure is described in Ref. [69].
Once the configurations had been ordered according to the Friedman's test values, one configuration of each analysed case was used for the final comparison taking two experiments into consideration: one looking at encodings and the other looking at to accuracy. In each case, additional information is given. The Hypervolume indicator average value evolution (in ten executions) is shown for each configuration. Moreover, box plots are given for the Hypervolume values distribution after the stopping criterion is met. In addition, the Friedman's test is used to detect significant differences among the performance of the configurations for each experiment. Finally, the accumulated best non-dominated solutions obtained (non-dominated front) are shown.
In Section 6.1, the results (of optimising the system that illustrate the case study) obtained after using the NSGA-II with different encodings (Real, Standard Binary and Gray) are presented; next, in Section 6.2 different accuracy levels of time (hours, days and weeks) in the conditions explained above are analysed and finally, in Section 6.3, the accumulated non-dominated designs are presented.

Encoding Experiment
The results of the Real encoding are presented first. Second, results of the binary encoding are shown, followed by the results of the Gray encoding. Finally, a comparison of the best performance cases of each codification is presented in Section 6.1.4.

Real Encoding
The results of using real encoding with simulated binary crossover are shown below. The computational time consumed is shown in Table 3. The average time represents the computational time regarding each one of ten executions and nine different configurations (real time consumed). The sequential time represents the computational time that would have been needed in the case of not using the cluster. The computational time shows the importance of using the cluster, which enables parallel processes. The relationship between method configurations (where N represents the population size and PrM the mutation probability) and identifiers is shown in Table 4. Moreover, statistical information in relation to the Hypervolume value at the end of the evolutionary process is shown (average, median, maximum, minimum and standard deviation, out of 10 independent executions). It is possible to conclude that the configuration with the identifier ID9 (with a population of 150 individuals and mutation probability of 1.5 gene per chromosome) presents the highest Hypervolume average value, the highest Hypervolume median value and the highest Hypervolume maximum value. The configuration with identifier ID3 (population of 150 individuals and mutation probability of 0.5 gene per chromosome) presents the highest Hypervolume minimum value and the configuration with identifier ID4 (population of 50 individuals and mutation probability of one gene per chromosome) presents the lowest Hypervolume standard deviation value. In order to establish the best behaviour amongst the configurations, a statistical significance hypothesis test was conducted. The average ranks computed through the Friedman's test and the p-value obtained (a value bigger than 0.05 implies that the null hypothesis cannot be rejected, suggesting that all configurations perform in a similar way) are shown in Table 4. It can be seen that the configuration with identifier ID4 (population of 50 individuals and mutation probability of one gene per chromosome) presents the best average rank (in order to maximise the Hypervolume, the average rank has to be as low as possible) so it was selected for the final comparison study among encoding configurations.
Finally, the best accumulated non-dominated solutions obtained through the final generation of the evolutionary process for all executions and all configurations were used to compute the accumulated Hypervolume (as described in Ref. [69]) whose value was 2.3943. As expected, the value is higher than 2.3307, the maximum value shown in Table 4.

Standard Binary Encoding
The results of using standard binary encoding with one point, two point and uniform crossover are shown below. The computational time consumed by each one is shown in Table 3. The relationship between method configurations and identifiers is shown in Table 5. Moreover, statistical information relating to the Hypervolume value at the final of the evolutionary process is shown.
For the binary encoding with one point crossover (B1PX), it is possible to conclude that the configuration with identifier ID8 (population of 100 individuals and mutation probability of 1. For the binary encoding with uniform crossover (BUX), it is possible to conclude that the configuration with identifier ID4 (population of 50 individuals and mutation probability of one gene per chromosome) presents the highest Hypervolume average value, the highest Hypervolume maximum value and the highest Hypervolume minimum value. The configuration with identifier ID2 (population of 100 individuals and mutation probability of 0.5 gene per chromosome) presents the highest Hypervolume median value and the configuration with identifier ID6 (population of 150 individuals and mutation probability of one gene per chromosome) presents the lowest Hypervolume standard deviation value. In order to establish the best behaviour amongst the configurations, a statistical significance hypothesis test was conducted. The average ranks computed through the Friedman's test and the p-value obtained (a value bigger than 0.05 implies that the null hypothesis cannot be rejected, suggesting that all configurations perform in a similar way) are shown in Table 5. It can be seen that the configuration with identifier ID8 (population of 100 individuals and mutation probability of 1.5 gene per chromosome) presents the best average rank for the binary encoding with one point crossover. It can be seen that the configuration with identifier ID2 (population of 100 individuals and mutation probability of 0.5 gene per chromosome) presents the best average rank both for the binary encoding with two point crossover and for the binary encoding with uniform crossover. These configurations were selected for the final comparison study among encoding configurations which is shown later. The best accumulated non-dominated solutions obtained through the final generation of the evolutionary process for all executions and all configurations were used to compute the accumulated Hypervolume whose values were 2.4142, 2.4298 and 2.3984 for the binary encoding with one point, two point and uniform crossover, respectively. As expected, the values are higher than 2.3394, 2.3627 and 2.3497, the maximum values shown in Table 5, respectively.

Gray Encoding
The results of using Gray encoding with one point, two point and uniform crossover are shown below. The computational time consumed by each one is shown in Table 3. The relationship between method configurations and identifiers is shown in Table 6. Moreover, statistical information in relation to the Hypervolume value at the end of the evolutionary process is shown. For the Gray encoding with one point crossover (G1PX), it is possible to conclude that the configuration with identifier ID2 (population of 100 individuals and mutation probability of 0.5 gene per chromosome) presents the highest Hypervolume average value, the highest Hypervolume median value and the highest Hypervolume maximum value. The configuration with identifier ID4 (population of 50 individuals and mutation probability of one gene per chromosome) presents both the highest Hypervolume minimum value and the lowest Hypervolume standard deviation value.
For the Gray encoding with two point crossover (G2PX), it is possible to conclude that the configuration with identifier ID1 (population of 50 individuals and mutation probability of 0.5 gene per chromosome) presents the highest Hypervolume average value, the highest Hypervolume median value and the highest Hypervolume minimum value. The configuration with identifier ID9 (population of 150 individuals and mutation probability of 1.5 gene per chromosome) presents the highest Hypervolume maximum value and the configuration with identifier ID5 (population of 100 individuals and mutation probability of one gene per chromosome) presents the lowest Hypervolume standard deviation value.
For the Gray encoding with uniform crossover (GUX), it is possible to conclude that the configuration with identifier ID4 (population of 50 individuals and mutation probability of one gene per chromosome) presents both the highest Hypervolume average value and the highest Hypervolume median value. The configuration with identifier ID5 (population of 100 individuals and mutation probability of 1 gene per chromosome) presents the highest Hypervolume maximum value, the configuration with identifier ID6 (population of 150 individuals and mutation probability of one gene per chromosome) presents the highest Hypervolume minimum value and the configuration with identifier ID2 (population of 100 individuals and mutation probability of 0.5 gene per chromosome) presents the lowest Hypervolume standard deviation value.
In order to establish the best behaviour amongst the configurations, a statistical significance hypothesis test was conducted. The average ranks computed through the Friedman's test and the p-value obtained (a value bigger than 0.05 implies that the null hypothesis cannot be rejected, suggesting that all configurations perform in a similar way) are shown in Table 6. It can be seen that the configuration with identifier ID5 (population of 100 individuals and mutation probability of 1 gene per chromosome) presents the best average rank for the Gray encoding with one point crossover. It can be seen that the configuration with identifier ID1 (population of 50 individuals and mutation probability of 0.5 gene per chromosome) presents the best average rank for the Gray encoding with two point crossover. It can be seen that the configuration with identifier ID4 (population of 50 individuals and mutation probability of one gene per chromosome) presents the best average rank for the Gray encoding with uniform crossover. These configurations were selected for the final comparison study among encoding configurations, which is shown later. The best accumulated non-dominated solutions obtained through the final generation of the evolutionary process for all executions and all configurations were used to compute the accumulated Hypervolume whose values were 2.4011, 2.3982 and 2.3829 for the Gray encoding with one point, two point and uniform crossover, respectively. As expected, the values are higher than 2.3556, 2.3364 and 2.3165, the maximum values shown in Table 6, respectively. The total computational time consumed is shown in Table 3. The computational cost shows the importance of using the cluster. Previously, configurations with the best average rank according to the Friedman's test were selected to be compared globally. These configurations are shown in Table 7. The Hypervolume average values evolution versus the evaluations number is shown in Figure 2a. The detail for the final evaluations (last million fitness function evaluations, from 9 to 10 million) is shown in Figure 2b. It can be seen that the configuration with identifier ID3 (with binary encoding and two point crossover, population of 100 individuals and mutation probability of 0.5 gene per chromosome) reaches the highest Hypervolume average value. Box plots of the Hypervolume values distribution at the end of the process are shown in Figure 3. Statistical information in relation to the Hypervolume value at the final of the evolutionary process is shown in Table 7. It can be seen that the configuration with identifier ID3 (with Binary encoding and two point crossover, population of 100 individuals and mutation probability of 0.5 gene per chromosome) presents both the highest Hypervolume average value and the highest Hypervolume median value, the configuration with identifier ID5 (with Gray encoding and one point crossover, population of 100 individuals and mutation probability of one gene per chromosome) presents the highest Hypervolume maximum value, the configuration with identifier ID6 (with Gray encoding and two point crossover, population of 50 individuals and mutation probability of 0.5 gene per chromosome) presents the highest Hypervolume minimum value and the configuration with identifier ID1 (with real encoding, population of 50 individuals and mutation probability of one gene per chromosome) presents the lowest Hypervolume standard deviation value. In order to establish whether one of the configurations performs better than any other, a statistical significance hypothesis test was conducted. The average ranks computed through the Friedman's test are shown in Table 7. It can be seen that the configuration with identifier ID3 (with binary encoding and two point crossover, population of 100 individuals and mutation probability of 0.5 gene per chromosome) presents the best average rank. However, the p-value computed (0.5979) implies that the null hypothesis (H 0 ) cannot be rejected (p-value > 0.05), so it is possible to conclude that, in the studied conditions, there is no configuration that performs better than another.
The best accumulated non-dominated solutions obtained for all encodings and configurations were used to compute the accumulated Hypervolume, whose value was 2.4553. As expected, the value is higher than 2.4298, the maximum accumulated value obtained after the evolutionary process for the standard binary encoding with two point crossover. This is shown in Table 8.

Accuracy Experiment
In the first experiment, a thorough comparison of the performances of encodings was developed using the hour as a time unit. Although non-significant differences among performances were found, the best average rank using the Friedman's test was presented by the standard binary encoding. For this reason, the results achieved for the standard binary encoding are used in the present experiment to compare the performance using the day and the week as time units.

Standard Binary Encoding (Days)
The results of using standard binary encoding with one point, two point and uniform crossover with the day as a time unit are shown below. The computational time consumed by each one is shown in Table 9. The relationship between method configurations and identifiers is shown in Table 10.  In order to establish the best behaviour amongst the range of configurations, a statistical significance hypothesis test was conducted. The average ranks computed through the Friedman's test and the p-value obtained (a value bigger than 0.05 implies that the null hypothesis cannot be rejected, suggesting that all configurations perform in a similar way) are shown in Table 10. It can be seen that the configuration with identifier ID6 (population of 150 individuals and mutation probability of one gene per chromosome) presents the best average rank for the binary encoding with one point crossover and the day as a time unit. Moreover, the p-value obtained of 0.0246 explains that the null Hypothesis can be rejected so, in this case, the configuration ID6 performs better than any other configuration. In order to find the concrete pairwise comparisons that produce differences, a statistical significance analysis using the Wilcoxon signed-rank test was run as explained by Benavoli et al. [70]. The results of using such a test, in which the configuration ID6 is compared with the rest of configurations, is shown in Table 11. It can be seen that configuration ID6 performs better than the configurations ID3, ID4 and ID8. Regarding the binary encoding with two point crossover and the day as a time unit, the same configuration (ID6) presents the best average rank. Finally, it can be seen that the configuration with identifier ID1 (population of 50 individuals and mutation probability of 0.5 gene per chromosome) presents the best average rank for the binary encoding with uniform crossover and the day as a time unit. In each case, the best configurations were selected for the final comparison study between the accuracy level encodings.
The best accumulated non-dominated solutions obtained through the final generation of the evolutionary process for all executions and all configurations were used to compute the accumulated Hypervolume whose values were 2.4007, 2.3942 and 2.3984 for the binary encoding with one point, two point and uniform crossover with the day as a time unit, respectively. As expected, the values are higher than 2.3430, 2.3372 and 2.3461, the maximum values shown in Table 10, respectively.

Standard Binary Encoding (Weeks)
The results of using standard binary encoding with one point, two point and uniform crossover and the week as a time unit are shown below. The computational time consumed by each one is shown in Table 9. The relationship between method configurations and identifiers is shown in Table 12. Statistical information in relation to the Hypervolume value at the final of the evolutionary process is also shown. For the binary encoding with one point crossover and the week as a time unit (B1PX-W), it is possible to conclude that the configuration with identifier ID1 (population of 50 individuals and mutation probability of 0.5 gene per chromosome) presents the highest Hypervolume average value, the highest Hypervolume median value and the highest Hypervolume maximum value. The configuration with identifier ID5 (population of 100 individuals and mutation probability of one gene per chromosome) presents both the highest Hypervolume minimum value and the lowest Hypervolume standard deviation value.
For the binary encoding with two point crossover and the week as a time unit (B2PX-W), it is possible to conclude that the configuration with identifier ID7 (population of 50 individuals and mutation probability of 1.5 gene per chromosome) presents both the highest Hypervolume average value and the highest Hypervolume maximum value. The configuration with identifier ID1 (population of 50 individuals and mutation probability of 0.5 gene per chromosome) presents both the highest Hypervolume median value and the highest Hypervolume minimum value. Finally, the configuration with identifier ID5 (population of 100 individuals and mutation probability of one gene per chromosome) presents the lowest Hypervolume standard deviation value.
For the binary encoding with uniform crossover and the week as a time unit (BUX-W), it is possible to conclude that the configuration with identifier ID8 (population of 100 individuals and mutation probability of 1.5 gene per chromosome) presents the highest Hypervolume average value, the highest Hypervolume median value and the highest Hypervolume minimum value. The configuration with identifier ID1 (population of 50 individuals and mutation probability of 0.5 gene per chromosome) presents the highest Hypervolume minimum value and the configuration with identifier ID2 (population of 100 individuals and mutation probability of 0.5 gene per chromosome) presents the lowest Hypervolume standard deviation value. In order to establish the best behaviour amongst configurations, a statistical significance hypothesis test was conducted. The average ranks computed through the Friedman's test and the p-value obtained (a value bigger than 0.05 implies that the null hypothesis cannot be rejected, suggesting that all configurations perform in a similar way) are shown in Table 12. It can be seen that the configuration with identifier ID5 (population of 100 individuals and mutation probability of one gene per chromosome) presents the best average rank for the binary encoding with one point crossover and the week as a time unit. It can be seen that the configuration with identifier ID1 (population of 50 individuals and mutation probability of 0.5 gene per chromosome) presents the best average rank for the binary encoding with two point crossover and the week as a time unit. It can be seen that the configuration with identifier ID8 (population of 100 individuals and mutation probability of 1.5 gene per chromosome) presents the best average rank for the binary encoding with uniform crossover and the week as a time unit. These configurations were selected for the final comparison study of the accuracy-level configurations.
The best accumulated non-dominated solutions obtained through the final generation of the evolutionary process for all executions and all configurations were used to compute the accumulated Hypervolume whose values were 2.4047, 2.4285 and 2.4037 for the binary encoding with one point, two point and uniform crossover with the week as a time unit, respectively. As expected, the values are higher than 2.3430, 2.3465 and 2.3336, the maximum values shown in Table 12, respectively.

Comparing Accuracy-Level Configurations
The global computational time consumed is shown in Table 9. The computational cost shows the importance of using the cluster. Previously, configurations with the best average rank according to the Friedman's test were selected to be globally compared. These configurations are shown in Table 13. The Hypervolume average values evolution versus the evaluations number is shown in Figure 4a. The detail for the final evaluations (last million fitness function evaluations, from 9 to 10 million) is shown in Figure 4b. It can be seen that the configuration with identifier ID2 (with Binary encoding, two point crossover, the hour as a time unit, population of 100 individuals and mutation probability of 0.5 gene per chromosome) reaches the highest Hypervolume average value.  Box plots of the Hypervolume values distribution at the end of the process are shown in Figure 5. They are ordered from left to right in relation to time units of hours (H), days (D) and weeks (W) and crossover types of one point (1PX), two point (2PX) and uniform crossover (UX). It can be seen that the medians are ordered from biggest to smallest for each group of crossover. The greater the accuracy the bigger the Hypervolume median value (it is greater for hours than for days and greater for days than for weeks). Statistical information in relation to the Hypervolume value at the end of the evolutionary process is shown in Table 13. It can be seen that the configuration with identifier ID2 (with binary encoding, two point crossover and the hour as a time unit, population of 100 individuals and mutation probability of 0.5 gene per chromosome) presents both the highest Hypervolume average value and the highest Hypervolume median value, the configuration with identifier ID9 (with binary encoding, uniform crossover and the week as a time unit, population of 100 individuals and mutation probability of 1.5 gene per chromosome) presents both the highest Hypervolume maximum value and the highest Hypervolume minimum value. The configuration with identifier ID7 (with binary encoding, one point crossover and the week as a time unit, population of 100 individuals and mutation probability of one gene per chromosome) presents the lowest Hypervolume standard deviation. In order to establish whether one of the configurations performs better than any other, a statistical significance hypothesis test was conducted. The average ranks computed through the Friedman's test are shown in Table 13. It can be seen that the configuration with identifier ID2 (with Binary encoding, two point crossover and the hour as a time unit, population of 100 individuals and mutation probability of 0.5 gene per chromosome) presents the best average rank. However, the p-value computed (0.4053) implies that the null hypothesis (H 0 ) cannot be rejected (p-value > 0.05), so it is possible to conclude that, in the studied conditions, there is no configuration that performs better than any other.
The best accumulated non-dominated solutions obtained were used to compute the accumulated Hypervolume, whose value was 2.4046. As expected, the value is higher than 2.4298, the maximum accumulated value obtained after the evolutionary process for the standard binary encoding with two point crossover and the hour as a time unit. This is shown in Table 14.

Accumulated Non-Dominated Set of Designs
The non-dominated solutions to the problem provided at the end of the evolutionary process for all executions, all configurations, all encodings and time units are shown in Figure 6. All optimum solutions belonging to the achieved non-dominated front are shown in Table 15. Unavailability (Q) is shown as a fraction, Cost is shown in economic units and the rest of the variables represent, for the respective devices, the optimum Time To Start a Preventive Maintenance activity with the hour, day or week as a time unit.  The solution with the lowest Cost (ID1) (823.38 economic units) represents the biggest Unavailability (0.002720). These values are followed by periodic optimum times (using the hour as a time unit in this case) measured from the moment in which the system mission time starts (Time To Perform a Preventive Maintenance activity (TRP) is not included). For solution ID1, it can be seen that periodic optimum Times To Start a Preventive Maintenance activity (TP) for devices P2 and V4 are not supplied. This is because the design alternative does not include such devices. The opposite case shows the biggest Cost (ID22) (1770.12 economic units) and the lowest Unavailability (0.000725). For solution ID22, periodic optimum Times To Start a Preventive Maintenance activity (TP) are supplied for all devices. This is because the design alternative includes devices P2 and V4. Other optimum solutions were found in these two solutions and can be seen in Table 15. The decision makers will need to decide which is the preferable design for them, taking into account their individual requirements.
Moreover, solutions were clustered in Figure 7 according to their final design. Solutions are shown from left to right and in ascending order in relation to the Cost from ID1 to ID22. The solutions contained in Cluster 1 (the solutions 1 to 2, see also Table 15) are the solutions in which non-redundant devices were included in the design. In this case, the system contains exclusively devices placed in series. These solutions present the lowest Cost and the biggest Unavailability. The solutions contained in Cluster 2 (the solutions 3 to 5, see also Table 15) are the solutions in which a redundant valve was included in the design as a parallel device. These solutions present a bigger Cost and a lower Unavailability than the solutions contained in Cluster 1. The solutions contained in Cluster 3 (the solutions 6 to 11, see also Table 15) are the solutions in which a redundant pump was included in the design as a parallel device. These solutions present a higher Cost and a lower Unavailability than the solutions contained in Clusters 1 and 2. Finally, the solutions contained in Cluster 4 (the solutions 12 to 22, see also Table 15) are the solutions in which both a redundant valve and a redundant pump were included in the design as parallel devices. These solutions present the biggest Cost and the lowest Unavailability. The best accumulated non-dominated solutions obtained were used to compute the accumulated Hypervolume, whose value was 2.4651. As expected, the value is higher than the rest of the maximum accumulated values obtained after the evolutionary process for the encoding experiment and the accuracy experiment. This is shown in Table 16. This value is also higher than the value obtained in [6] and could be considered as an actual benchmark value of the case study.

Conclusions
In the present paper, a methodology presented previously by the authors [6] was used. This consists of coupling a Multi-Objective Evolutionary Algorithm and Discrete Simulation in order to optimise simultaneously both the system design (based on redundancy allocation) and its preventive maintenance strategy (based on determining periodic preventive maintenance activities with regard to each device included in the design), whilst addressing the conflict between Availability and operational Cost. The Multi-Objective Evolutionary Algorithm gives rise to a population of individuals, each encoding one design alternative and its preventive maintenance strategy. Each individual represents a possible solution to the problem, which is then used to modify and evaluate the system Functionability Profile through Discrete Simulation. The individuals evolve generation after generation until the stopping criterion is reached. This process was applied to a technical system in a case study in which two experiments were developed: Firstly, an encoding experiment which consisted of comparing the performance of seven encoding types (real, standard binary with one point, two point and uniform crossover, and Gray with one point, two point and uniform crossover), and secondly, an accuracy level encoding which consisted of comparing the performance of using standard binary encoding with accuracy levels across a range of time unit (hours, days and weeks) with impact in the form of size of the chromosome (the smaller the time unit, the bigger the chromosome). The Multi-objective Evolutionary Algorithm used was NSGA-II and a set of optimum non-dominated solutions were obtained for all cases.
In conclusion, the use of the Multi-Objective Evolutionary Algorithm NSGA-II and Discrete Simulation to address the joint optimisation of the system design and its preventive maintenance strategy provides Availability-Cost-balanced solutions to the real world problem studied, where data based on specific bibliography, mathematical relations and field experience were used. The computational cost reveals the complexity of the process and the necessity of using a computing cluster, which allowed parallel executions.
With regard to the encoding experiment, the best ordered method by the Friedman test case based on the final Hypervolume indicator distributions was the two point crossover standard binary encoding, although no statistically significant differences were observed. With regard to the accuracy experiment, the best ordered method by the Friedman test case based on the final Hypervolume indicator distributions was the two point crossover standard binary encoding with hours as a time unit, although no statistically significant differences were observed. From the authors' point of view, an important conclusion arises from this experiment, which relates to flexibility regarding the time unit to schedule the preventive maintenance activities. Using the hour, the day or the week as a time unit does not affect significantly the performance of the configurations so, in the studied conditions, the preventive maintenance activities can be planned using weeks as a time unit. This allows a better range of time for planning than if days or hours are used as a time unit.
In addition, a higher benchmark value of the case study (in terms of Hypervolume indicator) is attained in this work as a reference.
As future work, the authors consider that these conclusions should be explored in greater depth extending the analysis to other real world problems in the reliability field, as well as comparing this analysis with other state of the art Multi-objective Evolutionary Algorithms.