Automatic Preference Based Multi-objective Evolutionary Algorithm on Vehicle Fleet Maintenance Scheduling Optimization
Yali Wang, Steffen Limmer, Markus Olhofer, Michael Emmerich, Thomas Baeck
AAutomatic Preference Based Multi-objectiveEvolutionary Algorithm on Vehicle Fleet MaintenanceScheduling Optimization
Yali Wang a , Steffen Limmer b , Markus Olhofer b , Michael Emmerich a , ThomasB¨ack a a Leiden Institute of Advanced Computer Science, Leiden University, 2333 CA Leiden, TheNetherlands b Honda Research Institute Europe GmbH, 63073 Offenbach, Germany
Abstract
A preference based multi-objective evolutionary algorithm is proposed for gener-ating solutions in an automatically detected knee point region. It is named Au-tomatic Preference based DI-MOEA (AP-DI-MOEA) where DI-MOEA standsfor Diversity-Indicator based Multi-Objective Evolutionary Algorithm). AP-DI-MOEA has two main characteristics: firstly, it generates the preference regionautomatically during the optimization; secondly, it concentrates the solutionset in this preference region. Moreover, the real-world vehicle fleet maintenancescheduling optimization (VFMSO) problem is formulated, and a customizedmulti-objective evolutionary algorithm (MOEA) is proposed to optimize main-tenance schedules of vehicle fleets based on the predicted failure distributionof the components of cars. Furthermore, the customized MOEA for VFMSOis combined with AP-DI-MOEA to find maintenance schedules in the auto-matically generated preference region. Experimental results on multi-objectivebenchmark problems and our three-objective real-world application problemsshow that the newly proposed algorithm can generate the preference regionaccurately and that it can obtain better solutions in the preference region. Es- ∗ Corresponding author
Email addresses: [email protected] (Yali Wang), [email protected] (Steffen Limmer),
[email protected] (MarkusOlhofer), [email protected] (Michael Emmerich), [email protected] (Thomas B¨ack)
Preprint submitted to Journal of L A TEX Templates January 26, 2021 a r X i v : . [ c s . N E ] J a n ecially, in many cases, under the same budget, the Pareto optimal solutionsobtained by AP-DI-MOEA dominate solutions obtained by MOEAs that pursuethe entire Pareto front. Keywords:
DI-MOEA, evolutionary algorithm, multi-objective optimization,preference, scheduling optimization
1. Introduction
The fleet maintenance scheduling optimization (VFMSO) problem was ini-tially proposed in [1] due to the increasing demand by companies, corporations,and organizations of all sorts which rely on vehicle fleets to deliver productsand services and need to maintain vehicles for safety reasons. In the problem, avehicle fleet, such as a taxi fleet, a bus fleet, etc., can be maintained in multipleseparate workshops according to a maintenance schedule. To be specific, eachworkshop has its own capacity and ability, meaning that on the one hand, eachworkshop has its own team and each team can work on only one car at the sametime; on the other hand, each workshop is limited to the maintenance of thespecific component(s) due to restrictions in the equipment or skill level of thestaff. The maintenance schedule is optimized for each component based on itsremaining useful lifetime (RUL) which has been predicted through predictiveapproaches or models [2]. Furthermore, the cost and time which are needed bydifferent workshops are considered because it is possible that the maintenance ofthe same component produces different costs and workloads when the operationis performed in different workshops. The VFMSO problem is essential becauseit can not only ensure the safety of vehicles for use; at the same time, it canlead to low maintenance costs and longer lives for vehicles as well.To enhance the approach in [1], to be specific, to handle the uncertaintyin the problem and apply it on new application scenarios, in this paper, weimprove it from the following two aspects:1. There exists a lot of uncertainty when we use the predicted RUL for eachcomponent as its due date, because no matter how accurate the predictive2odel is, it is still possible that the component will break on other dates:before the due date or later. Therefore, instead of only the RUL, we decideto involve the RUL probability distribution as the foundation to assign themaintenance time in the scheduling optimization.2. The VFMSO problem usually leads to a large and complex solution space,however, finding the most preferred solution is the ultimate goal. To thisend, AP-DI-MOEA (Automatic Preference based DI-MOEA) is developedbased on the framework of DI-MOEA (Diversity-Indicator based Multi-Objective Evolutionary Algorithm) [3]. The new algorithm can generatethe preference region automatically and find solutions with a more fine-grained resolution in the preference region.This paper is organized as follows. Section 2 formulates the enhancedVFMSO problem. A literature review on preference based optimization is pro-vided in Section 3. The customized multi-objective evolutionary algorithm forthe enhanced VFMSO is introduced in Section 4, and in Section 5, we explainAP-DI-MOEA. Section 6 presents and discusses experiments and their results.Lastly, Section 7 concludes the paper and outlines directions for future work.
2. Problem Formulation
For a car fleet operated by an operator, the components of cars (e.g., springs,brakes, tires or the engine) can fail and should be maintained regularly. Someseparate workshops are available for the maintenance of the car fleet, and therepair time and maintenance cost are known for each component in each work-shop. Beside the time and cost for repairing the car component, a fixed set-upcost and set-up time are also considered for each visit of a car to a workshop,which correspond to the cost and time required for the preparation of the main-tenance operation.The enhanced VFMSO problem addressed in this paper is defined as follows:1. There are n cars C = { C , C , · · · , C n } and m workshops W = { W , W , · · · , W m } . 3. Each car C i comprises l i components to be maintained for i = 1 , · · · , n .3. For each component O ij ( j = 1 , · · · , l i ), i.e., the j th component of car C i ,there is a set of workshops capable of repairing it. The set of workshopsis represented by W ij which is a subset of W .4. The processing time for maintaining component O ij in workshop W k ispredefined and denoted by p ijk .5. The cost for maintaining component O ij in workshop W k is predefinedand denoted by q ijk .6. The set-up time of car C i in workshop W k is predefined and denoted by x ik .7. The set-up cost of car C i in workshop W k is predefined and denoted by y ik .8. The number of teams in workshop W k is predefined and denoted by z k .9. The previous repair time of component O ij is recorded and denoted by L ij .At the same time, the following assumptions are made:1. All workshops and teams are available at the time of the optimization andassumed to be continuously available.2. All the components are independent from each other.3. Times required for transport of cars from/to workshops are included inthe maintenance time and cost of cars, and the set-up time.4. Environmental changes (such as car accidents) are not considered here.5. There are no precedence constraints among the components of differentcars. Cars are maintained on a first-come-first-served basis.6. Each team can only work on one operation at a time and an operation,once started, must run to completion.7. No operation can start before the completion of the previous operation.Two constraints are considered in the problem. As mentioned earlier, eachworkshop can only repair specific components, and this is the first constraint.4nother constraint is that the maintenance periods of different operations forthe same car should not overlap. It is obviously wrong if two overlapping main-tenance operations of a car are assigned to different workshops because one carcannot be in two different workshops at the same time. If two overlapping main-tenance operations of a car are assigned to the same workshop, it is not correcteither because these two maintenance operations should be grouped together asone operation in this case. The grouping strategy will be explained in Section4. Three objectives are assumed to be relevant for the vehicle fleet operator,which are the total workload, total cost and expected number of failures. Ina multi-objective optimization problem, the objectives typically are conflicting,i.e., achieving the optimal value for one objective requires some compromise onother objectives. In our problem, the fact that a faster maintenance usuallyis more expensive leads to the conflict between the first two objectives. Theexpected number of failures counts the times when the vehicles are broken onthe road. Here, the expected value is used because the actual value is unknownat the time of the optimization due to uncertainties in the predictions. Whenthe expected number of failures is large, less maintenance tasks are performed,therefore, the workload and cost can drop.Let T k denote the sum of the times spent for all operations that are processedin workshop W k ; M i the sum of all costs spent for all maintenance operationsof car C i ; F ij the number of failures of component O ij . Three objectives can bedefined as:Minimize the total workload: f = m (cid:88) k =1 T k (1)Minimize the total cost: f = n (cid:88) i =1 M i (2)Minimize the expected number of failures: f = n (cid:88) i =1 l i (cid:88) j =1 E ( F ij ) (3)5 . LITERATURE REVIEW Multi-objective scheduling optimization is a major topic in the research ofmanufacturing systems. Its fundamental task is to organize work and work-loads to achieve comprehensive optimization in multiple aspects, such as theprocessing time, processing cost and production safety, by deploying resources,setting maintenance time and processing sequence. In past decades, this issuehas received a great deal of interest and research in different fields, such asscheduling of charging/discharging for electric vehicles [4]; scheduling in cloudcomputing [5]; scheduling of crude oil operations [6]; scheduling in the manufac-turing industry to reduce carbon emissions [7]; scheduling medical treatmentsfor resident patients in a hospital [8]; scheduling for Internet service providers[9], and so on.As a typical workshop style, the flexible job shop scheduling problem (FJSP)is an essential branch of production planning problems. The FJSP consists ofa set of independent jobs to be processed on multiple machines, and each jobcontains several operations with a predetermined order. It is assumed thateach operation must be processed in specified processing time on a specific ma-chine out of multiple alternatives. The problem has been extensively studiedin the literature (for example, [10], [11], [12]). The FJSP is the research ba-sis of the maintenance scheduling optimization problem and many real-worldproblems extend the standard FJSP by adding specific features. [13] considersFJSP-PPF (process plan flexibility), where jobs can have alternative processplans. It is assumed that the process plans are known in advance and thatthey are represented by linear precedence relationships. Because only one ofthe alternative plans has to be adopted for each job, the FJSP-PPF deals withnot only routing and sequencing sub-problems, but also the process plan se-lection sub-problem. In this paper, a mixed-integer linear programming modelis developed for the FJSP-PPF. In [14], a mathematical model and a geneticalgorithm are proposed to handle the feature of overlapping in operations. It isassumed that a lot which contains a batch of identical items is transferred from6ne machine to the next only when all items in the lot have completed theirprocessing, therefore, sublots are transferred from one machine to the next forprocessing without waiting for the entire lot to be processed at the predecessormachine, meaning that starting a successor operation of job is not necessary tofinish of its predecessor completely. Three features are considered in [15], whichare (1) job priority; (2) parallel operations: some operations can be processedsimultaneously; (3) sequence flexibility: the sequence of some operations can beexchanged. A mixed integer liner programming formulation (MILP) model isestablished to formulate the problem and an improved differential evolution al-gorithm is designed. Because of unexpected events occurring in most of the realmanufacturing systems, there is a new type of scheduling problem known as thedynamic scheduling problem. This type of problem considers random machinebreakdowns, adding new machine, new job arrival, job cancellation, changingprocessing time, rush order, rework or quality problem, due date changing, etc.Corresponding works on the FJSP include [16], [17], [18], [19]. Compared withthe standard FJSP, our VFMSO problem has some special properties: (1) flex-ible sequence: the sequence of the components is not predefined, but mainlyinfluenced by the RUL probability distribution. (2) multiple problem param-eters: besides the processing time, other problem parameters like the mainte-nance cost, set-up time, set-up cost, repair teams, etc, also have impacts on theresult.Our real-world problem, like many other multi-objective optimization prob-lems, can lead to a large objective space. However, finding a well-distributedset of solutions on the Pareto front requests a large population size and com-putational effort. Therefore, instead of spreading a limited size of individualsacross the entire Pareto front, we decide to only focus on a part of the Paretofront, to be specific, the search for solutions will be only guided towards thepreference region which, in our algorithm, is determined by the knee point. Ithas been argued in the literature that knee points are most interesting solutions,naturally preferred solutions and most likely the optimal choice of the decisionmaker (DM) [20, 21, 22, 23]. 7he knee point is a point for which a small improvement in any objectivewould lead to a large deterioration in at least one other objective. In the lastdecade, several methods have been presented to identify knee points or kneeregions. Das [20] refers the point where the Pareto surface “bulges” the mostas the knee point, and this point corresponds to the farthest solution from theconvex hull of individual minima which is the minima of the single objectivefunctions. Zitzler [24] defines (cid:15) -dominance: a solution a is said to (cid:15) -dominate asolution b if and only if f i ( a ) + (cid:15) ≥ f i ( b ) ∀ i = 1 , ..., m where m is the numberof objectives. A solution with a higher (cid:15) -dominance value with respect to theother solutions in the Pareto front approximation, is a solution having highertrade-offs and in this definition corresponds to a knee point. The authors of[25] propose to calculate the density of solutions projected onto the hyperplaneconstructed by the extreme points of the non-dominated solutions, then identifythe knee regions based on the solution density.Different algorithms of applying knee points in MOEA have also been pro-posed. Branke [23] modifies the second criterion in NSGA-II [26], and replacesthe crowding distance by either an angle-based measure or a utility-based mea-sure. The angle-based method calculates the angle between an individual andits two neighbors in the objective space. The smaller the angle, the more clearlythe individual can be classified as a knee point. However, this method can onlybe used for two objective problems. In the utility-based method, a marginal util-ity function is suggested to approximate the angle-based measure in the case ofmore than two objectives. The larger the external angle between a solution andits neighbors, the larger the gain in terms of linear utility obtained from sub-stituting the neighbors with the solution of interest. However, the utility-basedmeasure is not suited for finding knees in concave regions of the Pareto front.Rachmawati [27, 28] proposes a knee-based MOEA which computes a trans-formation of original objective values based on a weighted sum niching approach.The extent and the density of coverage of the knee regions are controllable bythe parameters for the niche strength and pool size. The strategy is susceptibleto the loss of less pronounced knee regions.8ch¨utze [29] investigates two strategies for the approximation of knees ofbi-objective optimization problems with stochastic search algorithms. Severalnew definitions for identifying knee points and knee regions for bi-objectiveoptimization problems has been suggested in [30] and the possibility of applyingthem has also been discussed.Besides the knee points, the reference points, which are normally provided bythe DM, have also been used to find a set of solutions near reference points. Deb[31] proposes an MOEA, called R-NSGA-II, by which a set of Pareto optimalsolutions near a supplied set of reference points can be found. The dominancerelation together with a modified crowding distance operator is used in thismethodology. For all solutions of the population, the distances to all referencepoints are calculated and ranked. The lowest rank (over all reference points) ofa solution is used as its crowding distance. Besides, a parameter (cid:15) is used tocontrol the spread of obtained solutions. Bechikh proposes KR-NSGA-II [32] byextending R-NSGA-II. Instead of obtaining the reference points from the DM, inKR-NSGA-II, the knee points are used as mobile reference points and the searchof the algorithm was guided towards these points. The number of knee pointsof the optimization problem is needed as prior information in KR-NSGA-II.Gaudrie [33] uses the projection (intersection in case of a continuous front) ofthe closest non-dominated point on the line connecting the estimated ideal andnadir points as default preference. Conditional Gaussian process simulationsare performed to create possible Pareto fronts, each of which defines a samplefor the ideal and the nadir point, and the estimated ideal and nadir are themedians of the samples.Rachmawati and Srinivasan [34] evaluate the worthiness of each non-dominatedsolution in terms of compromise between the objectives. The local maxima isthen identified as potential knee solutions and the linear weighted-sums of theoriginal objective functions are optimized to guide solutions toward the kneeregions.Another idea of incorporating preference information into evolutionary multi-objective optimization is proposed in [35]. They combine the fitness function9nd an achievement scalarizing function containing the reference point. In thisapproach, the preference information is given in the form of a reference pointand an indicator-based evolutionary algorithm IBEA [36] is modified by embed-ding the preference information into the indicator. Various further preferencebased MOEAs have been suggested, e.g., [37, 38, 39].In our proposed algorithm, i.e., AP-DI-MOEA, we adopt the method from[20] to identify the knee point, design the preference region based on the kneepoint, and guide the search towards the preference region. The advantages of ouralgorithm are: (1) no prior knowledge is used in identifying the knee point andknee region; (2) the preference region is generated automatically and narroweddown step by step to benefit its accuracy; (3) our strategy cannot only handle bi-objective optimization problems, but also tri- and many-objective problems; (4)although we integrate the strategy with DI-MOEA, it may be integrated withany standard MOEAs (such as NSGA-II [26], SMS-EMOA [40] and others);(5) the proposed algorithm is capable of finding preferred solutions for multi-objective optimization problems with linear, convex, concave Pareto fronts anddiscrete problems.
4. Customized Algorithm for Vehicle Fleet Maintenance SchedulingOptimization
For our real-world VFMSO problem, we first define the execution windowfor each component based on its predicted RUL probability distribution whichis assumed to be a normal distribution. The execution window suggests that themaintenance of the component can only start at a time spot inside the window.The mean ( µ ) and standard deviation ( σ ) of the predicted RUL probabilitydistribution determine the interval of the execution window, which is definedas: [ µ − × σ , µ +2 × σ ]. The interval is chosen relatively long because 95% of thevalues are within two standard deviations of the mean, therefore, maintenancebefore or after the interval hardly makes sense.After the determination of the execution window, the following two special10trategies have been taken to improve the process of scheduling optimization: • Grouping components. • Obtaining the penalty cost and expected number of failures by MonteCarlo simulation.Lastly, evolutionary algorithm (EA) is chosen to solve this real-world applica-tion problem due to its powerful characteristics of robustness and flexibility tocapture global solutions of complex combinatorial optimization problems. More-over, EAs are well suited to solve multi-objective optimization problems due totheir ability to approximate the entire Pareto front in a single run.
It would be troublesome and also a waste of time and effort to send a car toworkshops repeatedly in a short period of time to repair different components.In our algorithm, since each component has its execution window for its main-tenance, it is possible to combine the maintenance of several components to onevisit if their execution windows overlap. Especially, by grouping the mainte-nance of multiple components into one maintenance operation, the set-up costand set-up time are charged only once for the complete group of components.
Figure 1: Possible groups for a car with eight components.
Figure 1 represents the execution windows of eight components of a car.The overlap of the execution windows shows the possibility of grouping thesecomponents. Combining components can only be effective if there is a common11verlap of the execution window for components from the same car, and thestarting time of the group operation must lie within the common overlap. Inthis example, component c c c Within the execution window of a component, an arbitrary time can bechosen as the starting time for maintaining the component. The maintenancetime of each component should be as close as possible to its real due date,because: • Performing the maintenance too early results in higher maintenance costsin the long term, because more maintenance tasks have to be done. • The risk of breaking down on the road will increase if the maintenancedate is too late.Therefore, we use Monte Carlo simulation to simulate the “real” due dates foreach component. In our experiments, stability can be achieved at a few hundredsamples, in our case, 1000 samples of the due date are generated in the executionwindow of each component according to its predicted RUL probability distri-bution (see Section IV). Figure 2 shows an example of the execution windowevolved from the predicted RUL probability distribution of a component. After1000 sampled due dates are generated in the execution window, the scheduledmaintenance date of the component is compared with these samples one by one,and each comparison can lead to three situations. Let us use d vij to denote the v th due date sample of component O ij ; and D ij the scheduled maintenance dateof component O ij . Three possibilities after the comparison are:Case 1) D ij < d vij The scheduled maintenance date is earlier than the sample (or the “real” duedate) means that the component will be maintained before it is broken. Inthis case, its useful life between the maintenance date and the due date will be12asted. Therefore, a corresponding penalty cost is imposed to reflect the waste.To calculate the penalty cost, a linear penalty function is suggested based onthe following assumptions: • If a component is maintained when it is new or the previous maintenancehas just completed, the penalty cost would be the full cost of maintainingit, which is c + s : the maintenance cost of the component and the set-upcost of the car; • If a component is maintained at exactly its due date, the penalty costwould be 0.
Figure 2: Execution window of a component.
Assume d vij is “Sampled Due date” in Figure 2, and D ij is “Maintenancedate a”, in this case, D ij is earlier than d vij . The penalty cost of “Maintenancedate a” for “Sampled Due date” would be the vertical dotted line above “Main-tenance date a”.Case 2) D ij > d vij The scheduled maintenance date is later than the sample means that the main-tenance date is too late and the defect occurs on the use. Still, d vij is “SampledDue date” in Figure 2, but the scheduled maintenance date D ij is “Maintenancedate b”. In this case, D ij is later than d vij , and the vehicle will break down onthe road. In our algorithm, the number of failures will be increased by one.Case 3) D ij = d vij To solve our application problem with an EA, there are several basic issueswe need to deal with, such as, how to represent an individual or solution inthe population (Chromosome Encoding); how to take these chromosomes into aprocess of evolution (Genotype-Phenotype Mapping); how to create variationsof solutions in each iteration (Genetic Operators). Details of these topics aregiven in the following subsections.
In our algorithm, a three-vector chromosome (Figure 3) is proposed to rep-resent an individual, and the three vectors are: • Group structure vector: the group structures of components. • Starting time vector: the starting times of operations. • Workshop assignment vector: the workshops for operations.The group structure vector gives the information of which components are inthe same group, it is initialized by randomly picking a feasible group structurefor each car (check the details in [1]). The generation of the starting time vectorshould be later than the generation of the group structure vector because the14 igure 3: Three-vector chromosome. starting time of each operation is determined by the execution window which isthe entire execution window of the component for a single-component operationor the execution window intersection for a group operation. A time spot israndomly selected from the execution window or execution window intersectionfor each operation in order to initialize the starting time vector.A workshop is considered as “several workshops” based on its capacity (thenumber of teams). By this way, the schedule for each workshop team can beachieved from the solution. For example, consider that two workshops havethree and four repairing teams respectively. Then, group operations can berandomly assigned to seven “workshops”, the former three and the latter fourrepresent the corresponding teams in two workshops.
To use the power of EAs to obtain a better population, we need to evalu-ate each chromosome and give the better ones higher probabilities to produceoffspring. This is done by genotype-phenotype mapping or decoding the chro-mosome. In our problem, it is to convert an individual into a feasible schedule tocalculate the objectives and constraints which represent the relative superiorityof a chromosome. The genotype-phenotype mapping can be easily achieved inour algorithm because the group structure, the starting time and the workshopteam of the operations can be acquired directly from each individual. Whenconverting an individual into a schedule, it is possible that the processing timesof two or more operations assigned to the same workshop team are overlappingsince the starting time of each operation is decided in the starting time vector.In this situation, the principle of first-come-first-served is followed: the startingtime and processing time of the earlier started operation remain the same; thestarting time of the later started operation is delayed until the completion of the15revious operation; the processing time of the later started operation remainsthe same; while, an extra waiting time is added to the later started operationas a penalty because the vehicle waits in the workshop for the maintenance.
In accordance with the problem and its encoding, specific crossover andmutation operators have been designed for our problem (check the details in [1]).Both operators are applied separately to the three parts of the chromosome.For the group structure vector, multi-point crossover can be used as crossoveroperator and the number of cutting points depends on the length of the vec-tor. The same cutting points can be applied to the starting time vector whenperforming crossover. However, the change on the group structure vector as aconsequence of the crossover may result in the invalidity of genes in the start-ing time vector because it is possible that the group members and executionwindow intersections have changed due to the new group structure. Therefore,when performing the crossover on the starting time vector, the starting times ofall operations should be checked based on the new group structure and a newstarting time is produced randomly from the correct intersection in the casethat the starting time of an operation is invalid. The multi-point crossover canbe applied to the workshop assignment vector as well.The mutation operator alters one or more gene values in a chromosome.Similarly, the mutation should be operated on the group structure vector firstdue to its impact on the starting time vector; the starting time of operationsshould be checked and corrected after the mutation is done on the group struc-ture vector. Afterwards, several gene values can be altered in the staring timevector and workshop assignment vector to generate a new individual.
5. Proposed Preference based Algorithm
As the number of objectives and decision variables increases, the numberof non-dominated solutions tends to grow exponentially [41]. This brings morechallenges on achieving efficiently a solution set with satisfactory convergence16nd diversity. At the same time, a huge number of solutions is needed to approx-imate the entire Pareto front. However, a big population means more computa-tional time and resources. To overcome these difficulties, we propose an auto-matic preference based MOEA, which can generate the preference region or theregion of interest (ROI) automatically and find non-dominated solutions in thepreference region instead of the entire Pareto front. The automatic preferencebased MOEA is developed based on the framework of DI-MOEA (Diversity-Indicator based Multi-Objective Evolutionary Algorithm) [3]. We call our newalgorithm AP-DI-MOEA.DI-MOEA is an indicator-based MOEA, it has shown to be competitive toother MOEAs on common multi-objective benchmark problems. Moreover, it isinvariant to the shape of the Pareto front and can achieve evenly spread Paretofront approximations. DI-MOEA adopts a hybrid selection scheme: • The ( µ + µ ) generational selection operator is used when the parent pop-ulation can be layered into multiple dominance ranks. The intention is toaccelerate convergence until all solutions are non-dominated. • The ( µ + 1) steady state selection operator is adopted in the case that allsolutions in the parent population are mutually non-dominated and thediversity is the main selection criterion to achieve a uniform distributionof the solutions on the Pareto front.DI-MOEA employs non-dominated sorting as the first ranking criterion; thediversity indicator, i.e., the Euclidean distance based geometric mean gap indi-cator, as the second, diversity-based ranking criterion to guide the search. Twovariants of DI-MOEA, denoted as DI-1 and DI-2, exist, which use the crowdingdistance and diversity indicator, respectively, as the second criteria in the ( µ + µ ) generational selection operator. While, to ensure the uniformity of thefinal solution set, the diversity indicator is used by both variants in the ( µ +1) steady state selection operator. Analogously, two variants of AP-DI-MOEA,i.e., AP-DI-1 and AP-DI-2, are derived from the two variants of DI-MOEA.17he workings of AP-DI-MOEA are outlined in Algorithm 1. (Exceedance of) Enum P is a predefined condition (In our algorithm,
Enum P is the numberof evaluations.) to divide the algorithm into two phases: learning phase anddecision phase. In the learning phase, the algorithm explores the possible area ofPareto optimal solutions and finds the rough approximations of the Pareto front.In the decision phase, the algorithm identifies the preference region and findspreferred solutions. When the algorithm starts running and satisfies
Enum P at some moment, the first preference region will be generated and
Enum P willbe updated for determining a new future moment when the preference regionneeds to be updated. The process of updating
Enum P continues until the end.The first
Enum P is a boundary line. Before it is satisfied, AP-DI-MOEA runsexactly like DI-MOEA to approximate the whole Pareto front; while, after itis satisfied, the preference region is generated automatically and AP-DI-MOEAfinds solutions focusing on the preference region. The subsequent values of
Enum P define the later moments to update the preference region step by step,eventually, a precise ROI with a proper size can be achieved.The first/new preference region is formed based on the population at themoment when the condition of
Enum P is satisfied, especially the knee pointof the population. Algorithm 2 gives the details of line 41 in Algorithm 1,it introduces the steps of finding the knee point of a non-dominated solutionset and constituting a hypercube shaped preference region according to theknee point. Figure 4 also gives an illustration of finding the knee point inbi-dimensional space. Firstly, the upper quartile objective values (line 13 inAlgorithm 2) in the solution set are used as a boundary to define outliers andsolutions outside this boundary are removed (line 16-20 in Algorithm 2). Theextreme solutions (the solutions with the maximum value in one objective) (line21 in Algorithm 2) are then found inside the boundary and a hyperplane isformed based on the extreme solutions. In a bi-dimensional space (Figure 4),the hyperplane is only a line connecting two extreme solutions. According to thenumbers of points below and above the hyperplane (line 22 - 23 in Algorithm2), the shape of the solution set can be roughly perceived. We will distinguish18 lgorithm 1
AP-DI-MOEA-1 P ← init (); // Initialize random population2: popsize ← | P | ;3: existRegion ← false; //mark whether there exists a preference region or not4: Enum C ← the number of evaluations so far;5: Enum P ← the number of evaluations before generating a (new) preference region;6: ( R , ..., R ‘ ) ← Partition P into fronts with increasing dominance ranks; //Non-dominated sorting7: for each i ∈ { , . . . , ‘ } do
8: calculate crowding distance for all solutions on R i ;9: end for t ← while Stop criterion not satisfied() do P t +1 ← ∅ ;13: if ( ‘ t > || t == 0) then
14: // ( µ + µ ) generational selection operator15: Q t ← Gen ( P t ); // Generate offspring with the size of popsize by recombination and mutation16: Evaluate Q t ;17: M t +1 = P t ∪ Q t // Combine offspring and parent population18: ( R , ..., R ‘ t +1 ) ← Partition M t +1 into fronts with increasing dominance ranks; //Non-dominated sorting19: i ← while | P t +1 | < popsize do P t +1 ← P t +1 ∪ R i ;22: i ← i + 1;23: end while if ( | P t +1 | > popsize ) then n ← | P t +1 | − popsize while n > do if ( existRegion == false) then
28: calculate crowding distance for all solutions on R i − ;29: remove the least contributor to rank and crowding distance in R i − from P t +1 ;30: else
31: calculate crowding distance and Euclidean distance for all solutions R i − ;32: remove the the least contributor to rank, crowding distance and Euclidean distance in R i − from P t +1 ;33: end if n ← n − end while end if else
38: // ( µ + 1) steady state selection operator39: if ( Enum C > Enum P ) then existRegion ← true;41: calculate P region ; //generate a (new) preference region42: update
Enum P ; //
Enum P ← the next number of evaluations to update the preference region43: end if q ← Gen ( P t ); // Generate only one offspring by variation45: Evaluate q ;46: P t +1 ← P t ∪ { q } ;47: ( R , ..., R ‘ t +1 ) ← Partition P t +1 into fronts with increasing dominance ranks; //Non-dominated sorting48: for each i ∈ { , . . . , ‘ t +1 } do if ( existRegion == true) then
50: calculate diversity indicator contribution and Euclidean distance for all solutions on R i ;51: else
52: calculate diversity indicator contribution for all solutions on R i ;53: end if end for if ( existRegion == true) then
56: remove the last solution from P t +1 based on rank, diversity indicator contribution and Euclidean distance;57: else
58: remove the last solution from P t +1 based on rank and diversity indicator contribution;59: end if end if t ← t + 1;62: l t ← the number of fronts of P t ;63: end while lgorithm 2 Finding the knee point and defining the preference region. n ← the number of objectives;2: P t ← current population;3: (cid:15) ; //parameter ( >
0) for distinguishing convex/concave shape;4: popsize ← | P t | ; //population size5: Declare Q [ n ]; //upper quartile objective values of P t
6: Declare L [ n ]; //worst objective values of P t
7: Declare knee [ n ]; //knee point of P t
8: Declare
P region [ n ]; //preference region of P t
9: Declare
Expoints [ n ][ n ]; //extreme points (single-objectives)10: foundknee ← false ;11: for each i ∈ { , . . . , n } do
12: sort( P t ) by the i th objective in ascending order;13: Q [ i ] ← P t .get index( × popsize ).get obj( i ); //upper quartile value of the i th objective14: L [ i ] ← P t .get index( popsize ).get obj( i ); //the largest (worst) value of the i th objective15: end for for all solution s ∈ P t do if s .get obj( i = 1 , ..., n ) > Q [ i ] then
18: remove s from P t ;19: end if end for Expoints [ (cid:5) ][ (cid:5) ] ← extreme points in P t ;22: num a ← the number of points in concave region of hyperplane formed by Expoints [ (cid:5) ][ (cid:5) ];23: num v ← | P t | − num a ; //the number of points in convex region24: if ( num v − num a > (cid:15) ) then
25: //roughly convex shape26: remove solutions in concave region from P t ;27: else if ( num a − num v > (cid:15) ) then
28: //roughly concave shape29: remove solutions in convex region from P t ;30: else
31: //roughly linear shape32: for all solution s ∈ P t do
33: calculate hypervolume of s with reference point L [ (cid:5) ];34: end for knee [ (cid:5) ] ← solution with the largest hypervolume value;36: foundknee ← true ;37: end if if ( foundknee == false ) then for all solution s ∈ P t do
40: calculate distance between s and hyperplane formed by Expoints [ (cid:5) ][ (cid:5) ];41: end for knee [ (cid:5) ] ← solution with the largest distance;43: end if for each i ∈ { , . . . , n } do P region [ i ] ← knee [ i ] + ( L [ i ] − knee [ i ]) × end for igure 4: Finding the knee point in bi-dimensional space. between “convex” and “concave” regions. Points in the convex ( concave ) region are dominating (dominated by) at least one point in the hyperplane spannedby the extreme points. However, when the number of the points in the convexregion and the number of points in the concave region is close enough, it impliesthat the shape of the current solution set is almost linear. This occurs bothwhen the true Pareto front is linear and when the solution set is converged verywell in a small area of the Pareto front. A parameter (cid:15) then is used to representthe closeness and it is a small number decided by the size of the solution set.In the case that the shape of the current solution set is (almost) linear, thesolution with the largest hypervolume value with regards to the worst objectivevector (line 14 in Algorithm 2) is adopted as the knee point (line 32 - 36 inAlgorithm 2). While, under the condition that the shape of the current solutionset is convex or concave, the solution in the convex or concave region with thelargest Euclidean distance to the hyperplane is chosen as the knee point (line39 - 42 in Algorithm 2). After the knee point is found, the preference regioncan be determined based on the knee point by the following formula: P region [ i ] = knee [ i ] + ( L [ i ] − knee [ i ]) ×
85% (4)21et i denotes the i th objective, as in Algorithm 2, L [ i ] is the worst valueof the i th objective in the population, knee [ i ] is the i th objective value of theknee point and P region [ i ] is the upper bound of the i th objective. W.l.o.g. weassume the objectives are to be minimized and the lower bound of preferenceregion is the origin point. According to the formula, we can see that the firstpreference region is relatively large (roughly 85% of the entire Pareto front).With the increase in the number of iteration, the preference region will beupdated and becomes smaller and smaller because every preference region picks85% of the current Pareto front. Eventually, we want the preference regioncan reach a proper range, say, 15% of the initial Pareto front. The process ofnarrowing down the preference region step by step can benefit the accuracy ofthe preference region.In the interest of clarity, Algorithm 1 only shows the workings of AP-DI-1, the workings of AP-DI-2 can be obtained by replacing crowding distancewith the diversity indicator contribution. In the ( µ + µ ) generational selectionoperator (line 14 - 36 in Algorithm 1), when there is no preference region,the second ranking criteria (the crowding distance for AP-DI-1; the diversityindicator for AP-DI-2) for all solutions on the last front are calculated and thepopulation will be truncated based on non dominated sorting and the secondranking criteria (line 28 - 29 in Algorithm 1). While, if a preference regionalready exists, both the second ranking criteria and Euclidean distance to theknee point for all solutions on the last front are calculated and the populationwill be truncated based on first non dominated sorting, then the second rankingcriteria, lastly, Euclidean distance to the knee point (line 31 - 32 in Algorithm1). In the ( µ + 1) steady state selection operator (line 38 - 59 in Algorithm 1),firstly, the value of Enum P is compared with the current number of evaluationsto determine if a (new) preference region should be generated. When it is timeto do so, the preference region is generated through Algorithm 2 (line 41 inAlgorithm 1), at the same time, the value of
Enum P is updated to the nextmoment when the preference region is to be updated (line 42 in Algorithm 1).There are different strategies to assign the values of
Enum P . In our algorithm,22e divide the whole computing budget into two parts, the first half is used tofind an initial entire Pareto front approximation, and the second half is used toupdate the preference region and find solutions in the preference region. Assumethe total computing budget is
Enum T (the number of evaluations), then thefirst value of
Enum P is × Enum T . Due to the reason that we expect afinal preference region with a size of around 15% of the initial entire Paretofront and each new preference region takes 85% of the current Pareto front,according to the formula: 0 . ≈ .
14, the value of
Enum P can be updatedby the following formula:
Enum P = Enum P + (
Enum T / /
12 (5)Another half of budget can be divided into 12 partial-budgets and a newpreference region is constituted after each partial-budget. In the end, the finalpreference region is achieved and solutions focusing on this preference region areobtained. For the rest part of the ( µ + 1) steady state selection operator, like-wise, when there is a preference region, three ranking criteria (1. non-dominatedsorting; 2. diversity indicator; 3. the Euclidean distance to the knee point) worktogether to achieve a well-converged and well-distributed set of Pareto optimalsolutions in the preference region.
6. Experimental Results
In this section, simulations are conducted to demonstrate the performanceof proposed algorithms on both benchmark problems and our real-world appli-cation problems. All experiments are implemented based on the MOEA Frame-work ( ), which is a Java-based framework formulti-objective optimization.For the two variants of AP-DI-MOEA: AP-DI-1 and AP-DI-2, their per-formances have been compared with DI-MOEA: DI-1, DI-2 and NSGA-III [42].23e compare our algorithm with NSGA-III because NSGA-III is a representativestate-of-the-art evolutionary multi-objective algorithm and it is very powerfulto handle problems with non-linear characteristics. For bi-objective benchmarkproblems, algorithms are tested on ZDT1 and ZDT2 with 30 variables. Forthree objective benchmark problems, DTLZ1 with 7 variables and DTLZ2 with12 variables are tested. For the real-world application problem of VFMSO,experiments have been conducted on two instances with different sizes. Theconfigurations of the two instances, such as the predicted RUL probability distri-bution, the processing time and maintenance cost of each component, the set-uptime and cost of each car, are made available on http://moda.liacs.nl . Onevery problem, we run each algorithm 30 times with different seeds, while thesame 30 different seeds are used for all algorithms. All the experiments are per-formed with a population size of 100; and for bi-objective problems, experimentsare run with a budget of 22000 (objective function) evaluations, DTLZ threeobjective problems with a budget of 120000 evaluations, the VFMSO problemswith a budget of 1200000 evaluations. This setting is chosen to be more realisticin the light of the applications in scheduling that we ultimately want to solve.
Bi-objective problems are optimized with a total budget of 22000 evaluations,when the number of evaluations reaches 10000 times, the first preference regionis generated, then after every 1200 evaluations, the preference region will beupdated. Figure 5 shows the Pareto front approximations from a typical run onZDT1 (left column) and ZDT2 (right column). The graphs on the upper row areobtained from DI-1 and AP-DI-1, while the graphs on the lower row are fromDI-2 and AP-DI-2. In each graph, the entire Pareto front approximation fromDI-MOEA and the preferred solutions from AP-DI-MOEA (or
AP solutions )are presented, at the same time, the preference region of AP-DI-MOEA is alsoshown by the gray area.Besides the visualization of the Pareto fronts, we also compute the knee pointof the entire final Pareto front approximation from DI-MOEA via the strategy24 a) ZDT1 (b) ZDT2Figure 5: Pareto front approximation on ZDT1 and ZDT2. • If the knee point from DI-MOEA is in the preference region achieved byits derived AP-DI-MOEA; • If the knee point from DI-MOEA is dominated by or dominating APsolutions; or if it is a non-dominated solution (mutually non-dominatedwith all AP solutions).Table 1 shows the results of 30 runs. For ZDT1 problem, all 30 knee pointsfrom DI-1 and DI-2 are in the preference regions from AP-DI-1 and AP-DI-2 re-spectively; in all these knee points, 10 from DI-1 and 7 from DI-2 are dominatedby AP solutions. For ZDT2 problem, most knee points are not in correspondingpreference regions, but for those in the preference regions, almost all of themare dominated by AP solutions. Please note that when a knee point from DI-MOEA is outside of the preference region from AP-DI-MOEA, it is not possiblethat it can dominate any AP solutions because all AP solutions are in the pref-erence region and only solutions in the left side of the gray area can dominateAP solutions.
Table 1: Space and dominance relation of knee point from DI-MOEA and AP solutions onZDT problems.
Problem ZDT1 ZDT2Algorithm DI-1/ DI-2/ DI-1/ DI-2/AP-DI-1 AP-DI-2 AP-DI-1 AP-DI-2In Incomparable 20 23 1 1preference Dominated 10 7 9 9region Dominating 0 0 0 0Outside Incomparable 0 0 20 20p-region Dominated 0 0 0 0We also perform the same comparison between AP-DI-MOEA and NSGA-III, the results are shown in Table 2. For ZDT1 problem, all knee points fromNSGA-III are in the preference regions from AP-DI-MOEA. Some of these kneepoints dominate AP solutions. For ZDT2 problem, most knee points from26SGA-III are not in the preference regions and these knee points are incompa-rable with AP solutions. For the knee points in the preference regions, all threedominating relations with AP solutions appear. For both problems, when theknee point from NSGA-III is dominating AP solutions, it only dominates oneAP solution.
Table 2: Space and dominance relation of knee point from NSGA-III and AP solutions onZDT problems.
Problem ZDT1 ZDT2Algorithm NSGA-III/ NSGA-III/ NSGA-III/ NSGA-III/AP-DI-1 AP-DI-2 AP-DI-1 AP-DI-2In Incomparable 14 19 3 1preference Dominated 0 0 2 3region Dominating 16 11 4 6Outside Incomparable 0 0 21 20p-region Dominated 0 0 0 0Instead of spreading the population across the entire Pareto front, we onlyfocus on the preference region. To ensure that our algorithm can guide the searchtowards the preference region and the achieved solution set is distributed acrossthe preference region, we compare the performance of AP-DI-MOEA, DI-MOEAand NSGA-III in the preference region. For each Pareto front approximationfrom DI-MOEA and NSGA-III, the solutions in the corresponding preferenceregion from AP-DI-MOEA are picked, and we compare these solutions with APsolutions through the hypervolume indicator. The point formed by the largestobjective values over all solutions in the preference region is adopted as thereference point when calculating the hypervolume indicator. It has been foundthat all hypervolume values of new solution sets from DI-MOEA and NSGA-IIIin the preference region are worse than the hypervolume values of the solutionsets from AP-DI-MOEA, which proves that the mechanism indeed works inpractice. Figure 6 shows box plots of the distribution of hypervolume indicatorsover 30 runs. 27 a) ZDT1 (b) ZDT2Figure 6: Boxplots comparing the hypervolume values on ZDT1 and ZDT2.
DTLZ1 and DTLZ2 are chosen as three objective benchmark problems toinvestigate our algorithms. They are performed with a total budget of 120000fitness evaluations, when the evaluation reaches 60000 times, the first preferenceregion is formed, then after every 5000 evaluations, the preference region isupdated. Figure 7 shows the Pareto front approximations from a typical runon DTLZ1 (left column) and DTLZ2 (right column). The upper graphs areobtained from DI-1 and AP-DI-1, while the lower graphs are from DI-2 andAP-DI-2. In each graph, the Pareto front approximations from DI-MOEA andcorresponding AP-DI-MOEA are given. Since the target region is actually anaxis aligned box, the obtained knee region (i.e., the intersection of the axisaligned box with the Pareto front) has an inverted triangle shape for these twobenchmark problems.Table 3 shows the space and dominance relation of the knee point fromDI-MOEA and the solution set from AP-DI-MOEA over 30 runs. For DTLZ128 a) DTLZ1 3 objective problem (b) DTLZ2 3 objective problemFigure 7: Pareto front approximation on DTLZ1 and DTLZ2.
Table 3: Space and dominance relation of knee point from DI-MOEA and AP solutions onDTLZ problems.
Problem DTLZ1 DTLZ2Algorithm DI-1/ DI-2/ DI-1/ DI-2/AP-DI-1 AP-DI-2 AP-DI-1 AP-DI-2In Incomparable 29 27 10 13preference Dominated 0 0 1 0region Dominating 0 0 0 0Outside Incomparable 1 3 19 17p-region Dominated 0 0 0 0AP-DI-1 and AP-DI-2 have also been compared with NSGA-III in the sameway. Table 4 shows the comparison result. For DTLZ1, the average numberof solutions from NSGA-III in the corresponding preference regions from AP-DI-MOEA is six. Still, almost all knee solutions from NSGA-III are in thepreference region. For DTLZ2, the average number of solutions from NSGA-III in the corresponding preference region from AP-DI-MOEA is less than one,30hile, in more than half of 30 runs, the knee points from NSGA-III are still inthe preference region. To some extent, it can be concluded that the preferenceregions from AP-DI-MOEA are accurate. It can also be observed that AP-DI-1 behaves better than AP-DI-2 on DTLZ2, because two knee points fromNSGA-III dominate the solutions from AP-DI-2.
Table 4: Space and dominance relation of knee point from NSGA-III and AP solutions onDTLZ problems.
Problem DTLZ1 DTLZ2Algorithm NSGA-III/ NSGA-III/ NSGA-III/ NSGA-III/AP-DI-1 AP-DI-2 AP-DI-1 AP-DI-2In Incomparable 30 29 14 17preference Dominated 0 0 1 1region Dominating 0 0 0 2Outside Incomparable 0 1 15 10p-region Dominated 0 0 0 0Similarly, we pick from DI-MOEA and NSGA-III solutions which are inthe corresponding preference region of AP-DI-MOEA, and the hypervolumeindicator value is compared between these solutions and AP solutions. It hasbeen found that all hypervolume values of solutions from AP-DI-MOEA arebetter than those of solutions from DI-MOEA and NSGA-III. The left columnof Figure 8 shows box plots of the distribution of hypervolume values over 30runs on DTLZ1, and the right column shows the hypervolume comparison onDTLZ2.In our experiments, we decide half of the total budget is used to find aninitial Pareto front because it turned out to be a good compromise: half budgetfor the initial Pareto front and another half budget for the solutions focusingon the preference region. We also run experiments using 25% and 75% of thetotal budget for the initial Pareto front. Figure 9 presents the entire Paretofront from DI-MOEA and the Pareto front from AP-DI-MOEA with differentbudgets for the initial Pareto front. The left two images are on DTLZ1 andthe right two images are on DTLZ2. The uppper two images are from DI-1and AP-DI-1; the lower two images are from DI-2 and AP-DI-2. In the legend31 a) DTLZ1 (b) DTLZ2Figure 8: Boxplots comparing the hypervolume values on DTLZ1 and DTLZ2. labels, 50%, 25% and 75% indicate the budgets which are utilized to find theinitial entire Pareto front. It can be observed that the preference region fromAP-DI-MOEA with 50% of budget are located on a better position than with25% and 75% budgets, and the position of the preference region from AP-DI-MOEA with 50% of budget is more stable. Therefore, in our algorithm, 50% ofbudget is used before the generation of preference region.
The budget of 1200000 evaluations has been used on the real-world applica-tion problems, and 600000 of them are for the initial Pareto front. After that,the preference region is updated after every 50000 evaluations. The VFMSOproblem has been tested with different sizes. Figure 10 shows Pareto front ap-proximations of a problem with 20 cars and 3 workshops (V1), and each carcontains 13 components: one engine, four springs, four brakes and four tires[43]. It can be observed that AP-DI-1 and AP-DI-2 can zoom in the entirePareto front and find solutions in the preference region, at the same time, both32 a) DTLZ1 3 objective problem (b) DTLZ2 3 objective problemFigure 9: Pareto front approximation by different budgets generating initial Pareto front. (a) DI-1 & AP-DI-1 (b) DI-2 & AP-DI-2Figure 10: Pareto front approximation on VFMSO problem with 20 cars and 3 workshops.
Table 5 gives the space and dominance relation of knee points from DI-MOEA and solutions from AP-DI-MOEA on these two VFMSO problems. Forboth problems, only few knee points from DI-MOEA are in the preference re-gions of AP-DI-MOEA, and the main reason is that the Pareto front of AP-DI-MOEA converges better than that of DI-MOEA, in some cases, the Pareto frontof DI-MOEA cannot even reach the corresponding preference region. More im-portantly, it can be observed that most knee points from DI-MOEA, no matter34 a) DI-1 & AP-DI-1 (b) DI-2 & AP-DI-2Figure 11: Pareto front approximation on VFMSO problem with 30 cars and 5 workshops.(a) V1 (b) V2Figure 12: Pareto front approximation on VFMSO problem by DI-MOEA, AP-DI-MOEA andNSGA-III.
Table 5: Space and dominance relation of knee point from DI-MOEA and AP solutions onV1 and V2.
Problem V1 V2Algorithm DI-1/ DI-2/ DI-1/ DI-2/AP-DI-1 AP-DI-2 AP-DI-1 AP-DI-2In Incomparable 0 0 0 0preference Dominated 9 7 9 6region Dominating 0 0 0 0Outside Incomparable 4 9 3 3p-region Dominated 17 14 18 21Table 6 gives the space and dominance relation of knee points from NSGA-III and AP solutions. For both problems, again, most knee points from NSGA-III are not in the preference regions of AP-DI-MOEA. Some knee points fromNSGA-III are dominated by AP solutions and most of them are incomparablewith AP solutions.
Table 6: Space and dominance relation of knee point from NSGA-III and AP solutions on V1and V2.
Problem V1 V2Algorithm NSGA-III/ NSGA-III/ NSGA-III/ NSGA-III/AP-DI-1 AP-DI-2 AP-DI-1 AP-DI-2In Incomparable 0 0 0 1preference Dominated 0 1 3 2region Dominating 0 0 1 1Outside Incomparable 23 24 21 18p-region Dominated 7 5 5 8
7. CONCLUSIONS
In this paper, a preference based multi-objective evolutionary algorithm,AP-DI-MOEA, is proposed. In the absence of explicitly provided preferences,36he knee region is usually treated as the region of interest or preference region.Given this, AP-DI-MOEA can generate the knee region automatically and canfind solutions with a more fine-grained resolution in the knee region. This hasbeen demonstrated on the bi-objective problems ZDT1 and ZDT2, and the threeobjective problems DTLZ1 and DTLZ2. In the benchmark, the new approachwas also proven to perform better than NSGA-III which was included in thebenchmark as a state-of-the-art reference algorithm.The research for the preference based algorithm was originally motivated bya real-world optimization problem, namely, Vehicle Fleet Maintenance Schedul-ing Optimization (VFMSO), which is described in this paper in a new formu-lation as a three objective discrete optimization problem. A customized set ofoperators (initialization, recombination, and mutation) is proposed for a multi-objective evolutionary algorithm with a selection strategy based on DI-MOEAand, respectively, AP-DI-MOEA. The experimental results of AP-DI-MOEA ontwo real-world application problem instances of different scales show that thenewly proposed algorithm can generate preference regions automatically and it(in both cases) finds clearly better and more concentrated solution sets in thepreference region than DI-MOEA. For completeness, it was also tested againstNSGA-III and a better approximation in the preference region was observed byAP-DI-MOEA .Since our real-world VFMSO problem is our core issue to be solved, and itsPareto front is convex, we did not consider problems with an irregular shape. Itwould be an interesting question how to adapt the algorithm to problems withmore irregular shapes. Besides, the proposed approach requires a definition ofknee points. Future work will provide a more detailed comparison of differentvariants of methods to generate knee points, as they are briefly introduced inSection 3. In the application of maintenance scheduling, it will also be importantto integrate robustness and uncertainty in the problem definition. It is desirableto generate schedules that are robust within a reasonable range of disruptionsand uncertainties such as machine breakdowns and processing time variability.37 cknowledgment
This work is part of the research programme Smart Industry SI2016 withproject name CIMPLO and project number 15465, which is (partly) financedby the Netherlands Organisation for Scientific Research (NWO).