Integration of electromagnetism with multi-objective evolutionary algorithms for RCPSP

Accepted Manuscript

Integration of electromagnetism with multi-objective evolutionary algorithms for RCPSP Jing Xiao, Zhou Wu, Xi-Xi Hong, Jian-Chao Tang, Yong Tang PII: DOI: Reference:

S0377-2217(15)00985-6 10.1016/j.ejor.2015.10.059 EOR 13343

To appear in:

European Journal of Operational Research

Received date: Revised date: Accepted date:

31 January 2014 22 July 2015 26 October 2015

Please cite this article as: Jing Xiao, Zhou Wu, Xi-Xi Hong, Jian-Chao Tang, Yong Tang, Integration of electromagnetism with multi-objective evolutionary algorithms for RCPSP, European Journal of Operational Research (2015), doi: 10.1016/j.ejor.2015.10.059

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Highlights • Electromagnetism (EM) is integrated with multi-objective evolutionary algorithms (MOEAs). • Performance of EM integration to NSGA-II, SPEA2 and MOEA/D for RCPSP is compared. • Employment of EM to NSGA-II is the most effective.

CR IP T

• EM does not improve the performance of MOEA/D on RCPSP.

AC

CE

PT

ED

M

AN US

• Explanation of why EM heuristics plays different roles in the studied algorithms is given.

1

ACCEPTED MANUSCRIPT

Integration of electromagnetism with multi-objective evolutionary algorithms for RCPSP Jing Xiao∗, Zhou Wu, Xi-Xi Hong, Jian-Chao Tang, Yong Tang

CR IP T

School of Computer Science, South China Normal University, Guangzhou 510631, China

Abstract

M

AN US

As one of the most challenging combinatorial optimization problems in scheduling, the resourceconstrained project scheduling problem (RCPSP) has attracted numerous scholars’ interest resulting in considerable research in the past few decades. However, most of these papers focused on the single objective RCPSP; only a few papers concentrated on the multi-objective resource-constrained project scheduling problems (MORCPSP). Inspired by a procedure called electromagnetism (EM), which can help a generic population-based evolutionary search algorithm to obtain good results for single objective RCPSP, in this paper we attempt to extend EM and integrate it into three reputable state-of-the-art multi-objective evolutionary algorithms (MOEAs) i.e. non-dominated sorting based multi-objective evolutionary algorithm (NSGA-II), strength Pareto evolutionary algorithm (SPEA2) and multi-objective evolutionary algorithm based on decomposition (MOEA/D), for MORCPSP. We aim to optimize makespan and total tardiness. Empirical analysis based on standard benchmark datasets are conducted by comparing the versions of integrating EM to NSGAII, SPEA2 and MOEA/D with the original algorithms without EM. The results demonstrate that EM can improve the performance of NSGA-II and SPEA2, especially for NSGA-II.

1. Introduction

PT

ED

Keywords: Multiple objective programming, Heuristics, Resource-constrained project scheduling, Electromagnetism, Evolutionary algorithms

AC

CE

The research on the classical resource-constrained project scheduling problem (RCPSP) has been widely expanded over the past several decades. According to the reviews provided by Herroelen et al. (1999, chap. 1), Brucker et al. (1999), Icmeli et al. (1993), Kolisch and Padman ¨ (2001) and OZdamar and Ulusoy (1995), the RCPSP can be stated as follows. A set of activities N numbered from 0 to n is to be scheduled without preemption. Activity i has a duration di and needs some resources which are limited during every moment to be completed. We assume that the renewable resource set is R and the availability of each resource type k, k ∈ R, which is a constant ak throughout the project horizon. Each activity i needs rik units of resource k during each period of its duration, i ∈ N, k ∈ R. The dummy start and end activities 0 and n have zero duration ∗ Corresponding author at: School of Computer Science, South China Normal University, Guangzhou 510631, China. Tel.: +86 13926287616 Email addresses: [email protected] (Jing Xiao), [email protected] (Zhou Wu), [email protected] (Xi-Xi Hong), [email protected] (Jian-Chao Tang), [email protected] (Yong Tang)

Preprint submitted to European Journal of Operational Research

October 30, 2015

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

and do not require any resource. The precedence relation among the activities is that one activity can not be started until all its predecessors have been finished. A solution for RCPSP is feasible if and only if the precedence and resource constraints are satisfied. In general, a solution of RCPSP is represented as a schedule or a list of start times s = (s0 , · · · , sn ) that implies a corresponding finishing times f = (f0 , · · · , fn ). Obviously, as every manager usually wants to finish projects as quickly as possible with minimum cost and high quality, RCPSP is a multi-objective optimization problem. Slowi´ nski (1981) is the first to describe the multi-objective resource-constrained project scheduling problem (MORCPSP) framework and list all kinds of objectives. After the further research on MORCPSP, the primary objectives include makespan, activity tardiness, net present value (NPV), resource investment (RI), robustness. MORCPSP is strongly NP-hard. As one kind of the most effective algorithms used for solving NP-hard problems, metaheuristic algorithms have been developed and proved to be successful in practice. Certainly many multiobjective metaheuristic approaches named MOEAs (multi-objective evolutionary algorithms) have been developed rapidly and applied for multi-objective optimization problems widely. Fonseca et al. (1993) proposed multi-objective genetic algorithm (MOGA); Srinivas and Deb (1994) proposed non-dominated sorting genetic algorithm (NSGA) and Horn et al. (1994) proposed niched Pareto genetic algorithm (NPGA). These algorithms were generally classified as the first generation of MOEAs in which selecting individuals based on Pareto ranks and maintaining population diversity by fitness sharing were the common features. The second generation of MOEAs based on elitism strategy has been proposed successively since 1999. Zitzler and Thiele (1999) presented Strength Pareto Evolutionary Algorithm (SPEA) and the improved version SPEA2 was described in Zitzler et al. (2001). Pareto Archived Evolution Strategy(PAES) was proposed by Knowles and Corne (2000). Pareto Envelope-Based Selection Algorithm (PESA) and its improved version of PESA-II were introduced in Corne et al. (2000) and Corne et al. (2001) respectively. Erickson et al. (2001, chap. 48) advanced NPGA2, and the famous NSGA-II was proposed by Deb et al. (2002). On the other hand, Zhou et al. (2007) presented Regularity Model Based Multi-objective Estimation of Distribution Algorithm (RM-MEDA), and Zhang and Li (2007) introduced Multiobjective Evolutionary Algorithm Based on Decomposition (MOEA/D) by combining traditional mathematical programming with evolutionary algorithms. Even though the MOEAs of the second generation have been proved to be successful in handling multi-objective optimization problems with high performance, only a few MOEAs such as NSGA-II have been used for dealing with MORCPSP. Meanwhile, many excellent hybrid algorithms have been presented and applied for classical single objective RCPSP successfully. As the electromagnetism (EM) heuristics has been applied for many unconstrained global optimization problem, Debels et al. (2006) extended the electromagnetism heuristics and integrated it into an scatter search (SS) frame for single objective RCPSP. By utilizing the principle that the EM force generated by two solutions would guide the search direction, scatter search integrated with EM outperformed other state-of-the-art heuristics for single objective RCPSP in this literature. However, to the best of our knowledge, no paper was found in literature integrating EM heuristics into MOEAs to solve MORCPSP. Inspired by the method presented by Debels, we modify the EM heuristics to make it suitable for MORCPSP and combine it with several state-of-the-art MOEAs i.e. NSGA-II, SPEA2 and MOEA/D in this paper. The paper is organized as follows: in Section 2, a literature review is presented. In Section 3, some basic concepts concerning MORCPSP are introduced. Then the representation and evaluation measurements for MORCPSP are defined in Section 4. In Section 5, we show how the EM 3

ACCEPTED MANUSCRIPT

methodology can be modified to be used for MORCPSP. Brief introductions of NSGA-II, SPEA2, MOEA/D and the method of integrating EM heuristics into the three algorithms are given in Section 6. Section 7 shows experimental results and comparison analysis. Finally, some concluding remarks are given in Section 8. 2. Literature review

AC

CE

PT

ED

M

AN US

CR IP T

Slowi´ nski (1981) presented the first multi-objective RCPSP framework that applied multiobjective linear programming to a MORCPSP with splittable activities, multiple modes, renewable, non-renewable and doubly constrained resources. It is notable that this procedure belonged to single objective optimization process actually because it could only optimize one objective each time. Furthermore, Slowi´ nski et al. (1994) proposed a decision support system to optimize several objectives including project completion time, smoothness of the resource profile, total resource consumption, weighted resource consumption, weighted flow time and net present value by using parallel priority rules, simulated annealing, branch and bound. Hapke et al. (1997) considered a multi-mode project scheduling problem under multi-category resource constraints with fuzzy time parameters of activities. A metaheuristic procedure Pareto simulated annealing (PSA) was presented to generate approximation solutions. Hapke et al. (1999, chap. 16), Pan and Yeh (2003, chap. 145) conducted the similar research. These three papers employed PSA frequently, and were partial to the fuzzy version of MORCPSP. Hapke et al. (1998) proposed a two-stage interactive search over a non-dominated solution space to handle a multiple-criteria PSP with multiple modes, renewable, non-renewable and doubly constrained resources. In order to optimize makespan, resource utilization smoothness, maximum lateness, NPV and project cost, in the first stage the approach applied PSA to generate a large representative sample of approximately nondominated schedules. An iterative search over the sample based on the discrete version of the Light Beam Search (LBS) procedure was organized in the secondary stage. Nabrzyski and W¸eglarz (1995) described a decision support system that could generate a set of feasible schedules using tabu search (TS) for individual criteria. In Nabrzyski and W¸eglarz (1999, chap. 17), they presented a knowledge-based multi-objective project scheduling system that handled a class of non-preemptive scheduling problems. The system made decisions during search. The optimization procedure was separated into several steps, and the trade-off solutions generated by each step guided the next step. These two models proposed by Nabrzyski considered the similar objective functions to Hapke et al. (1998), but difficult to obtain enough non-dominated solutions. In order to minimize the makespan, the ”weighted” lateness of activities and the violation of resource constraints, Viana and Pinho de Sousa (2000) applied multi-objective versions of simulated annealing (MOPSA) and tabu search (MOTS) to RCPSP. In consequence of the absence of resource constraints, the procedure generated different schedules; it seemed difficult to obtain excellent solutions for medium or large instances in various performance measurements. Al-Fawzan and Haouari (2005) took two objectives of makespan and robustness (maximize the sum of the free slack of activities) into account. A MOTS that ran a single-objective TS with a different linearly aggregated function several times was implemented for this model. The approach could obtain good solutions in terms of makespan while one of the drawbacks was that it would probably miss the solutions with good performance for robustness but not good enough at makespan. Abbasi et al. (2006) aggregated two objectives of minimizing makespan and maximizing robustness into a linear objective and applied SA to solve it. 4

ACCEPTED MANUSCRIPT

3. MORCPSP

PT

ED

M

AN US

CR IP T

Most of the aforementioned papers applied SA and TS to solve MORCPSP, but failed to obtain enough Pareto solutions. Apart from the methods motioned above, other algorithms have been used to solve MORCPSP. Kazemi and Tavakkoli-Moghaddan (2008) presented a mathematical model for MORCPSP with discounted cash flows that aimed at minimizing makespan and maximizing net present value. Due to the computational complexity of MORCPSP, the software proposed for solving the model only could handle the test instances with small number of activities and resources. NSGA-II was applied to this model for the same test instances as well. Aboutalebi et al. (2012) compared the performance of NSGA-II and Multi-objective Particle Swarm Optimization (MOPSO) for the bi-objective RCPSP. The adopted objectives were the minimization of project makespan and maximization of net present value. The computation result demonstrated the superior performance of NSGA-II in terms of metric C and maximum spread metric. Ballest´ın and Blanco (2011) focused on the theoretical study and gave some proof for MOPSP with regular objective functions(ROFs). Comparisons among NSGA-II, SPEA2 and PSA in several performance measures were presented. However, the representation of a solution was complicated and the measures to evaluate the quality of a solution obtained by algorithms were intricate because of all kinds of artificial parameters. Wang et al. (2013) presented a Pareto-archived estimationof-distribution algorithm (PAEDA) for the multi-objective resource-constrained project scheduling problem with makespan and resource investment criteria. The numerical experiment reuslts showed that the PAEDA outperformed the existing methods. Gomes et al. (2014) implemented five multi-objective metaheuristic algorithms, based on multi-objective GRASP (MOG), multiobjective variable neighborhood search (MOVNS), a MOG using NVS as local search (GMOVNS), a MOVNS with an intensification procedure (MOVNS I) and Pareto iterated local search (PILS), for MORCPSP to optimize makespan and total weighted start time of the activities. Four performance measurements: distance metric, hypervolume indicator, epsilon metric and error ratio were used for comparisons. A detailed experimental analysis indicated MOVNS as the most efficient one in terms of the distance metric, MOG and MOVNS to be superior in terms of the hypervolume and epsilon metrics and PILS as the best in terms of the error ratio. From all the reviews, the research on MORCPSP is still limited and there are only a few papers which applied MOEAs to MORCPSP. No paper is found in literature combining popular MOEAs i.e, NSGA-II, SPEA2 and MOEA/D with other heuristic approaches for dealing with MORCPSP.

CE

A multi-objective optimization problem (MOP) can be stated as follows: minimize F (x) = (f1 (x), · · · , fm (x)) subject to x ∈ Ω

(1)

∀i = 1, 2, · · · , m, fi (xa ) ≤ fi (xb ) ∧

(2)

AC

where Ω is the decision (variable) space, Rm is the objective space, and F : Ω → Rm consists of m objective functions. If the variable space is continuous, the MOP is considered as a multiobjective continuous optimization problem and the MOP is termed a multi-objective combinatorial optimization problem when the variable is discrete. Very often, no solution in Ω can minimize all the objectives simultaneously since the objectives in (1) contradict with each other. As a decision maker, the tradeoffs among the objectives have been taken into consideration usually. The best tradeoffs among the objectives can be defined in terms of Pareto optimality. Let xa , xb in Ω, xa dominates xb , (xa ≺ xb ) if and only if ∃j = 1, 2, · · · , m, fi (xa ) < fi (xb ) 5

ACCEPTED MANUSCRIPT

A solution x∗ ∈ Ω is Pareto optimal (non-dominated) if and only if ¬∃x ∈ Ω x ≺ x∗

(3)

AN US

CR IP T

F (x∗ ) is then called a Pareto optimal (objective) vector. The set of all Pareto optimal solutions (non-dominated solutions) is called Pareto set (PS) and the set of all Pareto optimal objective vectors is called Pareto front (PF), referring to Miettinen (1999). In MORCPSP, the precedence relationship among activities is illustrated by a directed acyclic graph and every node in the graph has a series of attributes associated with time and resource constraints. There are many objective functions of interest for managers. Some of them can be listed as follows in Table 1 by referring to Ballest´ın and Blanco (2011). Some problem notations are necessary for formal description: ddi and rdi denote the due date and release date of activity i respectively; wi is the weight associated with cost or reward; α is the discount rate. In the final formula in Table 1, the slack of i on s is the amount of time that activity i can slip without delaying the start of the next activities while maintaining resource feasibility. All the objectives aim to be minimized but N P V and robustness should be maximized. In this paper, we focus on the makespan fn and total tardiness (T W T without considering the weights).

PT

ED

Description makespan sum of start times weighted sum of start times maximum tardiness total amount of tardy activities total weighted tardiness total weighted flow time net present value robustness

M

Table 1: The descriptions of objective functions

Definition fn (the P finish time of end activity) SST = Pni=0 si W ST = ni=0 wi si M T = max{0, fi − ddi , i = 1, · · · , n} T AT A =P |{fi > ddi , i = 1, · · · , n}| n TWT = P i=0 wi max{0, fi − ddi } n TWFT = i (fi − rdi ) Pn i=0 w−αf i w e N P V = i=0 i Pn i=0 slack(i, s)

CE

4. Encoding for RCPSP

AC

Once the fundamental definitions of MORCPSP have been established, we discuss the representation scheme and evaluation procedure of solutions for RCPSP. It is well known that a solution for RCPSP is a schedule, but most metaheuristic algorithms for RCPSP do not operate on a schedule directly. It is convenient and effective when a solution is represented in some particular ways which can be transformed into a schedule using a scheduled generation scheme (SGS). According to Brucker et al. (1999), even though there are many representations for schedules of heuristics for the RCPSP, random-key (RK) and activity-list (AL) are the most important ones. The corresponding recombination methods generally include traditional one point or two point crossover and other proven effective conventional crossover such as parameterized uniform crossover(Mendes et al. (2009)). Due to the disadvantage that a single schedule can be represented by different activity lists or random keys, Debels et al. (2006) modified the standard random-key (SRK) representation in order to guarantee that each RK corresponded to a unique schedule. Taking a project depicted in Figure 1, for example, we can give a brief illustration for SRK. Assuming that activity i only 6

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

need one kind of renewable resource, and the amount of resource units required for activity i is ri . The duration and due date of activity i are di and ddi respectively; a feasible schedule could be depicted as Figure 2 if each type of resource availability is 2 (we omit the dummy start and end activity, horizontal axis represents the time and vertical axis represents the resource). According to the starting time of each activity, we assign a priority for each activity. The activity with the earliest starting time receives the highest priority denoted by 1, so we can obtain a RK vector x = [1; 5; 3; 2; 4; 6; 7; 8; 9] for the schedule in Figure 2. Because the priorities of the activities which have the same starting time can be exchanged without affecting the associated schedule, a schedule has multiple representations. A unique standardized random key (SRK) representation for schedule in Figure 2 can be gained when we attribute the lowest priority to all activities starting at the same time such as 1 for activity 1 and activity 4, 4 for activity 2 and activity 5. By doing this for x, we finally obtain the unique SRK vector x0 = [1; 4; 3; 1; 4; 6; 7; 8; 9]. For every new priority vector created in the algorithm, it is convenient to transform it into SRK form by using serial SGS. We can convert it to a SRK vector in accordance with the starting time obtained by serial SGS. Before we use SGS to generate SRK vector, we can get an initial random key chromosome vector RK. In the initial step, we can generate initial chromosomes by topological sorting. The precedence relationship between each task is depicted with an acyclic graph G which includes n vertices. Assuming that Idpv denotes the amount number of the precursor tasks of vertex v and Sucv represents all the successor tasks of vertex v, we can obtain an initial random key chromosome vector RK as follows. ∀v ∈ G, ∝v = Idpv ; s = {}, i = 1, prio = 1; Repeat until i = n V = {v ∈ / s |∝v = 0} Push all v ∈ V onto a stack in an arbitrary order; Pop a task w from the stack, i = i + 1; s = s ∪ {w}; RK[w] = prio, prio = prio + 1; ∀v ∈ Sucw , ∝v =∝v −1; Then, we can build a schedule from any one chromosome by using SGS procedure. The SGS procedure is briefly described as follows: P rew : the precursor tasks of task w; Sucw : the successor tasks of task w; sw : start time of task w; fw : finish time of task w; rwk : request quantity resource k of task w; Rk0 (t): current quantity of resource k at time t; ∀v ∈ G, ∝v = Idpv ; f easibletask = {}; For i = 1 to i = n V = {v ∈ / feasibletask |∝v = 0} Insert all v ∈ V into a queue q in ascendent order with RK[v]; Delete the head task w from queue q; if P rew = φ then sw = 0 else sw = maxj∈P rew {fj }; sw = min{t | rwk ≤ Rk0 (τ ) ∀k ∈ K, τ ∈ [t, t + dw ], t ≥ sw }; 7

ACCEPTED MANUSCRIPT

di ri 0 0

0

i

ddi 1 1

2

4

2 1

0

2 1

1

3

2

5

2 1

2

4

3 1

5

2

4

6 3

1

8 2 5

2

7

10 r

0

10

8

10

2

0

4

0 4 1

9

1

12 0

1

2

3

4

6

5 5

6

7

7 8

8

9

9 10 11 12 13 14 15 16 17 18 t

Figure 2: A schedule for the example project

CR IP T

Figure 1: Project example

2

3

1

AN US

fw = sw + dw ; Update resource profile; feasibletask = feasibletask ∪{w}; ∀v ∈ Sucw , ∝v =∝v −1; With the above SGS procedure, we can obtain the SRK chromosome when we give the tasks the same random keys if the tasks start at the same time. As we have known that a feasible schedule could be obtained by using serial SGS, in order ¨ to improve the solution quality, Li and Willis (1992), Ozdamar and Ulusoy (1996) presented an iterative scheduling technique. The principle contains two stages as follows.

ED

M

Stage 1. Backward shift: When a schedule for example as Figure 2 has been obtained by serial SGS, we consider the end time of activities as priority rule by assigning the highest priority to the activity with maximum end time. Let the maximum finish time be a start point, according to the above priority rule and the precedence relationships between activities, we can schedule the activities in the direction from right to left and obtain a new schedule. The new schedule is depicted in Figure 3. Stage 2. Forward shift: In reverse, we rank the priorities in descending order by the starting time of activities generated in Stage 1. With this order we schedule the activities from left to right to obtain another schedule illustrated by Figure 4.

AC

CE

PT

Assuming that dd = [2; 5; 5; 2; 4; 8; 10; 10; 12] is the due date set of this example project (see Figure 1), according to the objective functions, the makespan and total tardiness of the schedule depicted in Figure 2 are 18 and 23 respectively. The makespan and total tardiness of another schedule depicted in Figure 4 are equal to 15 and 21. It is obvious to show the improvement of the schedule by using backward and forward shifts. It has also been proved that the schedule makespan obtained by backward and forward shifts is never higher than the makespan of the schedule without shifting. In fact, we reduce the makespan by shifting each activity as much as possible to the right without affecting the project completion time and then shifting them to left as much as possible in r

r

2

2

1 0

1

2

3

5

3

1

2 4

5

6

7

8

6

2

9

8 4

6 4

3

1

5

1

7

9 10 11 12 13 14 15 16 17 18 t

0

Figure 3: A schedule generated by backward shift

1

2

3

4

5

6

7

8

7 8

9

9 10 11 12 13 14 15 16 17 18 t

Figure 4: A schedule generated by forward shift

8

ACCEPTED MANUSCRIPT

a similar way. Attributing to the left shift, the starting time of most activities will be minimized, of course the total tardiness of a schedule will be reduce too. In this paper, we iteratively perform backward and forward shifts until no further improvement can be found. 5. Extending EM heuristics for MORCPSP

CE

PT

ED

M

AN US

CR IP T

Before introducing the EM heuristics, we illustrate how to use traditional two-point crossover operator to recombine offspring. If p1 and p2 are two parents, we can get two offspring c1 and c2 with two-point crossover operation as follows. Randomly generate two integers γ1 and γ2 , 1 < γ1 < γ2 < n. For i ∈ [1, γ1 ] and i ∈ [γ2 + 1, n], c1 [i] = p1 [i], c2 [i] = p2 [i]; For i ∈ [γ1 + 1, γ2 ], c1 [i] = p2 [i], c2 [i] = p1 [i]. Consider as an example p1 = [1; 3; 6; | 1; 5; 4; | 7; 8; 9] p2 = [1; 1; 3; | 6; 4; 4; | 8; 7; 8] then c1 = [1; 3; 6; | 6; 4; 4; | 7; 8; 9] c2 = [1; 1; 3; | 1; 5; 4; | 8; 7; 8] Finally, we also should build schedules and get SRK chromosomes by using SGS on offspring. Electromagnetism (EM) heuristics for unconstrained global optimization problems is proposed by Birbil and Fang (2003). Each point represented a solution has a charge related to an objective function value in a multi-objective solution space by analogy with electromagnetism. Two arbitrary points in a population of solutions will generate a force with each other. In an evolutionary search algorithm, each point (solution) in the population will exert attraction or repulsion on other points with the force which is proportional to the product of the charges and inversely proportional to the distance between the points. The principle of EM heuristics to guide the search is that a point with a relatively good objective function value attracts the other ones and the points with inferior objective values repel the other population members. Based on the research of Birbil and Fang (2003), for RCPSP, Debels et al. (2006) modified and integrated the EM heuristics into SS procedure for new solution recombination. The main idea of his modification is as follows. Let xi and xj be two feasible solutions of RCPSP, which are represented by SRK. A charge qij that depends on the relative difference in objective function value between xi and xj can be defined as: f (xi ) − f (xj ) qij = (4) f (xworst ) − f (xbest )

AC

where xworst and xbest are the worst and best solutions in population. If f (xi ) < f (xj ), i.e. xi is better than xj , qij is negative and xj repels xi . The opposite, i.e. xj attracts xi when xi is worse than xj . The force exerted on solution xi by solution xj can be defined as: Fij = (xj − xi ) · qij

(5)

It is notable that xi and xj are vectors, so Fij is a vector as well. Finally, the search can move from solution xi to xi + Fij in the direction of xj . For the RCPSP, forces are exerted in changing in the priority of each activity. Two boundaries pmin ∈ [1; n − 1] and pmax ∈ [2; n] (n denotes the total amount of activities) are selected randomly and the SRK values between pmin and pmax are updated according to the force exerted in these 9

ACCEPTED MANUSCRIPT

dimensions. Be different from the activities with SRK values between pmin and pmax , the activities with SRK values lower than pmin or higher than pmax preserve the priorities structure from the original solution through subtracting a large constant from all SRK values lower than pmin and adding a large constant to all those larger than pmax . The last step is transforming the resulting RK vector into SRK format. A hybrid two-point/EM crossover operator can be described briefly as follows:

CR IP T

1. SRK < pmin : a large constant value is subtracted from the priority. 2. pmin ≤ SRK ≤ pmax : a force calculated by EM is added to the priority value. 3. SRK > pmax : a large constant value is added to the priority value.

Due to the definition of formula (4) considers only single objective, we extend it for MOP as: qij =

m X k=1

ϕk

fk (xi ) − fk (xj ) fkworst − fkbest

(6)

M

AN US

where m is the number of objectives, fk (xi ) is the k-th objective value of xi , fkworst and fkbest are the worst and best values of k-th objective, ϕk is a constant factor. We can obtain the worst and best values for each objective in the initial population and define them as the f worst and f best . In the next generations, we update the f worst and f best by obtaining the objective values worse than current f worst and better than current f best . Similar to single objective case, qij is positive and xj attracts xi in its direction if xi is dominated by xj . On the contrary, xj repels xi when xi dominates xj . For example, xj will guide the move direction of xi with the attractive force Fij when xj dominates xi in Figure 5.

2

6 4

3

1

ij

j

r 2

2

ED

i

5

1 0

1

2

3

4

5

6

7

7 8

9

9 10 11 12 13 14 15 16 17 18 t

8

PT

1

Figure 7: A schedule for the example project (associated with vector x2 )

Figure 5: xj attracts xi

CE

r 2

6

1 1

5

AC

4 2

3

0

1

2

3

4

5

6

7

8

r 2

1 7

8

4 2

9

9 10 11 12 13 14 15 16 17 18 19 t

0

1

2

5

3

1

3

4

5

6 6

7

8

7

8

9

9 10 11 12 13 14 15 16 17 18 t

Figure 6: A schedule for the example project (associated Figure 8: The resulting schedule after two-point/EM with vector x1 ) crossover operator(x00 )

We will demonstrate how to implement hybrid two-point/EM crossover operator by the case of example project depicted in Figure 1 in Section 3. Let x1 = [1; 3; 6; 1; 5; 4; 7; 8; 9] and x2 = [1; 1; 3; 6; 4; 4; 8; 7; 8] are two feasible SRK vectors associated with two schedules depicted in Figure 6 and Figure 7. According to the due date set dd = [2; 5; 5; 2; 4; 8; 10; 10; 12], two bi-objective solutions f (x1 ) = (19, 22) and f (x2 ) = (15, 21) are able to be calculated using the objective 10

ACCEPTED MANUSCRIPT

CR IP T

functions presented in Section 3. The first objective is makespan and total tardiness value is the second objective. We assume that 22 and 15 are the worst and best values of makespan, and the worst and best values of total tardiness equal to 25 and 18, so that the charge q12 = 5/7 ≈ 0.71. According to F12 = q12 (x2 − x1 ), force F12 on x1 can be calculated. In Table 2 (an adaptation of Debels et al. (2006)), the components of x1 between [pmin ; pmax ] are bolded and we can obtain a new 0 where F 0 is represented as a vector and calculated by updating F schedule x01 with x01 = x1 + F12 12 12 with substracting or adding a constant according to the hybrid two-point/EM crossover operator. Finally, x01 can be converted to its SRK form x001 by using SGS. The schedule corresponded to x001 is depicted in Figure 8, with f (x001 ) = (18, 21). This hybrid two-point/EM crossover operator on x1 not only preserves the part of SRK vector larger than pmax and lower than pmin but also results in the middle part between pmin and pmax attracted by x2 . Meanwhile, we can reproduce another offspring that copies the lowest and highest priorities from x2 , with a middle part repelled by x1 . The hybrid two point/EM crossover operator illustrated above will be used for generating offsprings instead of the traditional two-point crossover in the recombination steps of three MOEAs introduced in Section 6.

2 3 1 -1.42 -1.42 1.58 3

3 6 3 -2.13 -2.13 3.87 4

4 1 6 3.55 -10 -9 1

M

1 1 1 0 -10 -9 1

5 5 4 -0.71 -0.71 4.29 5

6 4 4 0 0 4 5

7 7 8 0.71 0.71 7.71 7

8 8 7 -0.71 +10 18 8

9 9 8 -0.71 +10 19 9

ED

Activities x1 x2 F12 0 F12 x01 x001

AN US

Table 2: Illustration of the execution of hybrid two-point/EM crossover operator

6. MOEAs integrated with EM heuristics

AC

CE

PT

In this section we will introduce three multi-objective evolutionary algorithms which are applied to solve MORCPSP. On the other hand, we also demonstrate how to integrate EM heuristics with MOEAs in the pivotal steps. NSGA-II and SPEA2 are two most popular algorithms which have been widely used for many MOPs successfully and effectively. By combining mathematical programming with evolutionary algorithm, MOEA/D decomposes a multi-objective optimization problem into a number of scalar optimization subproblems and optimizes them simultaneously. It has been shown that MOEA/D receives a fairly good effect in continuous MOP and discrete MOP. 6.1. Integration of EM with NSGA-II for RCPSP Firstly, the outline of the NSGA-II is as follows:

Input: N (population size) Output: N S (non-dominated set) Step 1: Create an initial population P0 of size N and an empty population Q0 . Set t=0 //initialization S Step 2: Rt = Pt Qt //combine parent and children population Step 3: F = fast-non-dominated-sort(Rt ) // F = (F1 , F2 , · · · ) all non-dominated fronts of Rt 11

ACCEPTED MANUSCRIPT

Step 4: Pt+1 = ∅ and i = 1 Until Pt+1 < N //till the parent population is filled crowding-distance-assignment (Fi ) //calculate crowding distance in Fi S Pt+1 = Pt+1 Fi //include i-th non-dominated front in the parent pop

CR IP T

5: sort(Pt+1 , n ) //sort in descending order using n ) 6: Pt+1 = Pt+1 [0 : N ] //choose the first N elements of Pt+1 7: If stopping condition is satisfied, N S is set of the non-dominated solutions in Pt+1 . Stop. 8: Qt+1 = make-new-pop (Pt+1 ) //using selection, traditional crossover(or hybrid twopoint/EM crossover) and mutation to create a new population Qt+1 Step 9: t = t + 1. Go to Step 2. Step Step Step Step

CE

PT

ED

M

AN US

Fast non-dominated sorting and crowding distance assignment are the primary parts of NSGAII. The population is classified into a number of layers by dominance ranks. In the first layer, F1 is consisted of non-dominated solutions in population. Each solution in the second layer F2 should be dominated by one solution in F1 at least. Similarly all solutions can be ranked and each solution does not dominate or be dominated by any other in the same layer. The new parent population Pt+1 is formed by adding solutions from the first front till the size exceeds N . Thereafter, the solutions of the last accepted front are sorted according to n , then the first N solutions are picked. The partial order (n ) can be established by the dominance rank and crowding distance of each solution. Given two individuals i and j, we have i n j if (irank < jrank ) OR (irank = jrank AND idistance > jdistance ). Individual i has a lower dominance rank when individual i dominates individual j. The partial order is defined as that the individual with lower dominance rank will be selected when the dominance ranks of two individuals are different. We prefer to choose the individual with larger crowding distance when two individuals have the same equal dominance rank, i.e. they do not dominate each other. To calculate the crowding distance, the population is sorted according to each objective value in ascending order. For each objective, the solutions with the smallest and largest values are assigned an infinite distance to guarantee that boundary points are always selected. All the other intermediate solutions are assigned a distance value equal to the absolute normalized difference in the objective value of two adjacent solutions. The overall crowding distance of a solution is the summation of distance values of each objective. If both two solutions belong to the same layer, the one with larger crowding distance will be selected. In Step 8, we use hybrid two-point/EM crossover operator (see Section 5) to generate new offsprings for new population when we combine EM heuristics with NSGA-II.

AC

6.2. Integration of EM with SPEA2 for RCPSP As another classical algorithm, the main frame of SPEA2 can be described as follows:

Input: N (population size) N (archive size) Output: N S (non-dominated set) Step 1: Generate an initial population P0 and empty archive p0 . Set t = 0; // initialization Step 2: Calculate fitness values of individuals in Pt and Pt ; // fitness assignment Step 3: Copy all non-dominated individuals in Pt and Pt to Pt+1 . If the size of Pt+1 exceeds N then reduce Pt+1 by truncation. Otherwise if the size of Pt+1 is less than N then fill Pt+1 with dominated individuals in Pt and Pt // environmental selection 12

ACCEPTED MANUSCRIPT

CR IP T

Step 4: If stopping criterion is satisfied then set NS to the set of decision vectors represented by the non-dominated individuals in Pt+1 . Stop. // termination Step 5: Perform binary tournament selection with replacement on Pt+1 in order to fill the mating pool.// mating selection Step 6: Apply recombination and mutation operators to the mating pool and set Pt+1 to the resulting population. Increment generation counter t = t + 1 and go to Step2.//variation, we can apply hybrid two-point/EM crossover for recombination here

AN US

In SPEA2, each individual i in population Pt and archive Pt will be assigned a strength S(i) representing the number of solutions it dominates and a density estimation D(i) representing the distance to its nearest k-th nearest data point. The fitness of individual i used for Step 3 to fill Pt+1 is the summation of strengths of solutions that dominates i, plus density estimation D(i). For the truncation operator, the k-th nearest distance is adopted: the individual which has the minimum distance to another individual is truncated at each stage; if there are several individuals with minimum distance, consider the second smallest distances and so forth. We apply hybrid two-point/EM algorithm to generate new solutions instead of the traditional two-point crossover in Step 6 when EM heuristics is integrated into SPEA2. 6.3. Integration of EM with MOEA/D for RCPSP

CE

PT

ED

M

Unlike NSGA-II and SPEA2, MOEA/D does not treat a MOP as a whole. MOEA/D explicitly decomposes the MOP into N scalar optimization subproblems which can be solved simultaneously by evolving a population of solutions. Let λ1 , · · · , λN be a set of uniform weight vector and z ∗ be the reference point. These parameters are used for defining the scalar functions of subproblems. With different decomposition approaches, the scalar functions are different. In MOEA/D, each weight vector λi has a neighborhood consisted of its several closest weight vectors in λ1 , · · · , λN . The neighborhood of weight vector λi decides the neighborhood of i-th subproblem. The neighborhood of each subproblem is used for generating new solutions. The population is composed of the best solution called current solution found so far for each subproblem. At each generation t, MOEA/D maintains a population of N solutions corresponding to N subproblems and z = (z1 , · · · , zm )T (a set of the best value found so far for each objective). F V is the objective function value of a solution. In the same way, the hybrid two-point/EM crossover operator will be used for recombination if we apply EM heuristics for MOEA/D in the part of reproduction in Step 2.1. The algorithm works as follows:

AC

Input: N the number of the subproblems considered in MOEA/D A uniform spread of N weight vector: λ1 , · · · , λN ; T the number of the weight vectors in the neighborhood of each weight vector. Output: EP Step 1: Initialization. Step 1.1: Set EP = ∅ Step 1.2: Compute the Euclidean distance between any two weight vectors and then work out the T closest weight vectors to each weight vector. For each i = 1, · · · , N , set B(i) = i1 , · · · , iT , where λi1 , · · · , λiT are the closest weight vectors to λi . Step 1.3: Generate an initial population x1 , · · · , xN randomly. Set F V i = F (xi ). Step 1.4: Initialize by a problem-specific method. 13

ACCEPTED MANUSCRIPT

Step 2: Update. For i = 1, · · · , N do

CR IP T

Step 2.1: Reproduction: Randomly select two indexes k, l from B(i), and then generate a new solution y from xk and xl by using recombination operators. // recombination operators include traditional crossover(or hybrid two-point/EM crossover) and mutation Step 2.2: Update of z: For each j = 1, · · · , m, if zj < fj (y) then set zj = fj (y). Step 2.3: Update of Neighboring Solutions: For each index j ∈ B(i), if y is better than xj for j-th subproblem, the set xj = y and F V j = F (y). Step 2.4: Update of EP : Remove from EP all the solutions dominated by F (y). Add F (y) to EP if no solution in EP dominates F (y). Step 3: If stopping criteria is satisfied, stop and output EP . Otherwise, go to Step 2.

AN US

7. Computational experiments

We implement NSGA-II, SPEA2, MOEA/D and the EM integration versions of these three algorithms in C for Windows. The performance of these algorithms are compared and analyzed too. As test instances, we choose 160 instances from the standard j30,j60,j90,j120 set, each set contains 40 instances, introduced in Kolisch et al. (1995) for the RCPSP and borrow the due date generated by Ballest´ın et al. (2006, chap. 4) where the procedure to generate the due date in dataset is presented.

ED

M

7.1. Performance measures In our experiments, the following performance measures are adopted to evaluate the solution sets produced by different evolutionary algorithms.

PT

(i) Set Coverage (C-Metric), which computes the percentage of solutions of one set dominated by solutions of another set. Let A and B be two sets of solutions. The set coverage is defined as: |{u ∈ B | ∃v ∈ A : v ≺ u}| C(A, B) = × 100% (7) |B|

AC

CE

The value C(A, B) = 100% means that all solutions of B are dominated by solutions of A and C(A, B) = 0 implies no solution of B is dominated by solution of A. It is notable that C(A, B) is not necessarily equal to 100% − C(B, A). (ii) Distance from the PF (D-Metirc) is a value representing how ”far” an approximation set is from the pareto front. Assuming P ∗ be a set of uniformly distributed points along the PF. Let A be a set of solutions, the average distance from A to P ∗ is defined as: D(A, p∗ ) =

∀v∈ A d(v, p∗ ) |A|

(8)

where d(v, p∗ ) is the minimum Euclidean distance between v and the points in P ∗ . If P ∗ is large enough to represent the PF very well, set A must be very close to the PF when D(A, p∗ ) is low enough. In our experiments, we do not know the actual PF. The non-dominated solutions obtained from all algorithms are used as the reference set P ∗ for each instance. On the other hand, we incorporate a simple objective normalization technique into D-metric measure since 14

ACCEPTED MANUSCRIPT

the difference between two objectives is large in the MORCPSP. A simple normalization method is to replace objective fi (i = 1, · · · , m) by fi =

fi − zimin zimax − zimin

(9)

CR IP T

where zimax and zimin is the maximum and minimum value of i-th objective in reference set P ∗. 7.2. Results and discussion

ED

M

AN US

In our experimental studies, for standard j30,j60,j90,j120 set, the number of iterations is 50,75,100,125,the population size is 100,200,300,400. ϕk for every objectives in formula (6) is equal to 1.0. We use Tchebycheff approach and set neighborhood size T to 10 for MOEA/D. The mutation rate of 0.06 is considered for all algorithms. 10 runs are performed for each algorithm on each instance. Table 3 presents the average mean of D-Metric values of the final solutions for each dataset compared to the original algorithms without EM. The better D-Metric values for each dataset is highlighted in bold. This table reveals that both NSGA-II EM (NSGA-II integrated with EM heuristics) and SPEA2 EM (SPEA2 integrated with EM heuristics) obtain lower D-Metric values than original NSGA-II and SPEA2 for j120 dataset, and higher D-Metric values for j30,j60 datasets. MOEA/D EM (MOEA/D integrated with EM heuristics) does not outperform MOEA/D for all datasets. In order to evaluate the convergence of all algorithms, Figure 9 and Figure 10 present the change of average D-Metric values along with the increasing number of schedules on some test instances. We find that MOEAs integrated with EM heuristics converges faster than their original versions in most j90,j120 cases. As the number of activities grows, NSGA-II EM outperforms other algorithms. Table 3: Average D-Metric values of the solutions from standard datasets found by all algorithms

NGSA-II EM

NSGA-II

PT

Dataset

0.1907 0.7130 2.6945 0.3490

0.1152 0.6359 2.9409 0.4538

SPEA2

MOEA/D EM

MOEA/D

0.2264 0.9355 3.7074 0.4335

0.1141 0.9092 3.4522 0.4947

0.4534 1.1111 4.0784 0.8140

0.4330 1.0228 4.0384 0.8073

CE

j30 j60 j90 j120

SPEA2 EM

AC

Figure 11 shows the potential non-dominated solutions obtained by 10 runs of each algorithm for the instances of dataset j120 on Figure 10. RPF is the reference pareto front set and it consists of all the solutions using the six algorithms studied in this paper for that instance for each run. From Figure 10, we can observe that almost all solutions gained by NSGA-II are dominated by solutions founded by NSGA-II EM. SPEA2 EM works out a majority of solutions which dominate the solutions obtained by SPEA2, but most of the solutions gained by MOEA/D and MOEA/D EM are dominated by each other. It means that the EM heuristics can not help to improve the performance of MOEA/D. Figure 12,13 shows that in terms of C-Metric, the final solutions obtained by NSGA-II are better than NSGA-II EM in most j30,j60 cases,but in most j90,j120 cases,NSGA-II EM is better. SPEA2 EM gets similar results but MOEA/D EM can not receive better performance than MOEA/D. 15

ACCEPTED MANUSCRIPT

0.6

NSGA-II

NSGA-II_EM

CR IP T

NSGA-II

NSGA-II_EM

1.2

SPEA2

SPEA2

SPEA2_EM

1.0

average D-metric value

MOEA/D MOEA/D_EM

0.4

0.2

MOEA/D MOEA/D_EM

0.8

0.6

AN US

average D-metric value

SPEA2_EM

0.4

0.2

0.0 0

1

2

3

4

5

0

2

3

M

NSGA-II_EM SPEA2 SPEA2_EM

ED

average D-metric value

MOEA/D MOEA/D_EM

PT

1.5

0.5

0

CE

1.0

4

number of schedules

AC

4

12

16

MOEA/D MOEA/D_EM

3

2

1

0 8

x10

NSGA-II

5

SPEA2_EM

2.0

6 3

SPEA2

2.5

5

(b) j3045 10

NSGA-II

3.0

4

number of schedules

NSGA-II_EM

average D-metric value

1

x10

(a) j309 3 3.5

6 3

number of schedules

0

4

8

3

x10

(c) j609 8

number of schedules

(d) j6045 9

Figure 9: The average D-Metric values of all algorithms for some test instances

16

12

16 3

x10

ACCEPTED MANUSCRIPT

14 13

SPEA2

12

SPEA2_EM

6

average D-metric value

MOEA/D MOEA/D_EM

5 4 3

NSGA-II_EM SPEA2

SPEA2_EM

11

MOEA/D

10

MOEA/D_EM

9 8 7 6 5 4

2 1

NSGA-II

AN US

average D-metric value

7

NSGA-II_EM

CR IP T

NSGA-II

8

3

0

6

12

18

24

2

30

0

6

12

18

24

30

3

x10

number of schedules

1.1

j1207_1

j12039_1

NSGA-II

NSGA-II

NSGA-II_EM

1.0

NSGA-II_EM 0.5

SPEA2

SPEA2

SEPA2_EM

SEPA2_EM

ED

0.9

average D-metric value

MOEA/D MOEA/D_EM

0.8

0.7

0.6

PT

average D-metric value

x10

(b) j9045 2

M

(a) j909 10

3

number of schedules

0.5

0.4

0.2 0

CE

0.3

10

20

MOEA/D_EM 0.4

0.3

0.2

30

number of schedules

AC

MOEA/D

40

50

0

10

20

30

40

50 3

3

number of schedules

x10

(c) j1207 1

(d) j12039 1

Figure 10: The average D-Metric values of all algorithms for some test instances

17

x10

ACCEPTED MANUSCRIPT

j12039_1 j1207_1

120

RPF

160

100

140

60

120 100

40

80

20

60

106

108

110 f

1

RPF SPEA2

M 104

106

108

ED

20

110

f

j1207_1

160

112

102

103

104

105

106

RPF SPEA2 SPEA2_EM

120

AC

80 60 95

114

98

99

100

101

102

103

104

105

106

1

RPF

180

MOEA/D

MOEA/D MOEA/D_EM

MOEA/D_EM

160

f2

140 120 100 80

20

60 106

97

j12039_1

RPF

40

104

96

f

CE

140

0

101

1

j12039_1

1

PT

180

60

100

100

40

80

99

180

f2

f2 60

100

98

140

80

120

97

160

100

0

96

SPEA2_EM

120

f2

95

114

f

j1207_1

140

112

AN US

104

CR IP T

f2

80

f2

NSGA-II NSGA-II_EM

NSGA-II_EM

0

RPF

180

NSGA-II

108

110 f

112

95

114

96

97

98

99

100

101 f

1

102

103

104

105

106

1

Figure 11: Non-dominated solutions obtained by 10 runs of each algorithm on j1207 1(left) and j12039 1(right)

18

ACCEPTED MANUSCRIPT

C(NSGA-II_EM, NSGA-II)

C(NSGA-II_EM, NSGA-II)

C(NSGA-II, NSGA-II_EM)

%

%

80

90

75 70

80 Average C-metric value

60 55 50 45 40 35 30 25 20 15 10

70 60 50 40

CR IP T

65

Average C-metric value

C(NSGA-II, NSGA-II_EM)

100

85

30 20 10

5 0

0

-5 0

4

8

12

16

20

24

28

32

36

40

44

0

4

8

Instance

20

%

28

32

36

40

44

C(SPEA2, SPEA2_EM)

AN US

85

75 70 65 60 55 50 45 40 35 30

M

25 20 15 10 5

Average C-metric value

100

80

80

60

40

20

0

4

8

12

16

20

ED

0

0 -5

24

28

32

36

40

0

44

4

8

12

16

20

%

PT

C(MOEA/D, MOEA/D_EM)

AC

40

20

4

8

40

44

70 60 50 40 30 20

0

0

36

80

Average C-metric value

60

32

C(MOEA/D, MOEA/D_EM)

90

CE

80

28

C(MOEA/D_EM, MOEA/D)

C(MOEA/D_EM, MOEA/D)

%

24

Instance

Instance

Average C-metric value

24

C(SPEA2_EM, SPEA2)

C(SPEA2, SPEA2_EM)

90

Average C-metric value

16

Instance

C(SPEA2_EM, SPEA2)

%

12

12

16

20

24

28

32

36

40

44

Instance

10

0

4

8

12

16

20

24

28

32

36

40

Instance

Figure 12: Average set coverage between algorithms integrated with EM and their original versions for j30(left) and j60(right)

19

44

ACCEPTED MANUSCRIPT

C(NSGA-II_EM, NSGA-II)

C (N S G A -II_ E M , N S G A -II) C (N S G A -II, N S G A -II_ E M )

%

C(NSGA-II, NSGA-II_EM)

1 0 0

%

9 0

100

7 0

60

40

20

6 0 5 0 4 0 3 0

CR IP T

80

A v e r a g e C - m e tr ic v a lu e

Average C-metric value

8 0

2 0 1 0 0

0

0

0

4

8

12

16

20

24

28

32

36

40

4

8

44

1 6

2 0

2 4

C(SPEA2_EM, SPEA2)

3 2

3 6

4 0

C (S P E A 2 _ E M , S P E A 2 ) C (S P E A 2 , S P E A 2 _ E M )

AN US

%

C(SPEA2, SPEA2_EM)

%

2 8

In s ta n c e

Instance

1 0 0

100

60

M

40

20

0 0

4

8

12

16

20

24

28

32

36

40

A v e r a g e C - m e tr ic v a lu e

8 0

80

ED

Average C-metric value

1 2

6 0

4 0

2 0

0

44

4

8

1 2

1 6

2 0

2 4

2 8

3 2

3 6

4 0

In s ta n c e

Instance

C (M O E A /D _ E M , M O E A /D ) C (M O E A /D , M O E A /D _ E M )

PT

C(MOEA/D_EM, MOEA/D)

%

C(MOEA/D, MOEA/D_EM)

100

CE

90

70

AC

60

8 0

A v e r a g e C - m e tr ic v a lu e

80 Average C-metric value

% 9 0

50 40

7 0

6 0

5 0

30 20

4 0

0

4

8

12

16

20

24

28

32

36

40

44

0

4

8

1 2

1 6

2 0

2 4

2 8

3 2

3 6

4 0

In s ta n c e

Instance

Figure 13: Average set coverage between algorithms integrated with EM and their original versions for j90(left) and j120(right)

20

ACCEPTED MANUSCRIPT

7.3. Experiments given the same running time

CR IP T

From the above experiments, integrating EM with NSGA II and SPEA2 has a better performance in standard j120 set. Table 4 shows the running time and we can see that integrating EM does not bring computing burden. We also compare NSGA-II and SPEA2 with or without EM given the same running time. Table 5 shows the D-Metric values of dataset j120 with the same running time. Lower D-Metric is highlighted in bold. The results show that NSGA II EM is better than NSGA II in most instances, while SPEA2 EM performs slightly better than SPEA2. Figure 14 shows the C-Metric values. We can get the same conclusion as we discussed D-Metric values in Table 5. Table 4: Average running time of the algorithms for j120 (s)

NGSA-II EM

NSGA-II

SPEA2 EM

SPEA2

j1201 1 j1203 1 j1203 4 j1207 1 j1207 2 j1207 3 j1208 2 j1208 4 j1209 1 j12013 1 j12013 3 j12013 5 j12014 1 j12014 4 j12017 1 j12017 4 j12018 1 j12018 2 j12019 4 j12019 6 j12020 1 j12021 6 j12022 6 j12028 2 j12029 5 j12033 1 j12034 4 j12034 6 j12035 2 j12035 3 j12037 6 j12038 1 j12039 1 j12039 2 j12040 1 j12040 2 j12049 3 j12053 3 j12055 4 j12058 6

53.79 50.65 47.51 60.50 59.92 59.15 60.60 57.49 55.93 63.29 61.74 53.92 57.55 56.67 66.48 60.65 61.29 60.21 55.49 53.62 54.07 49.58 47.25 54.32 53.06 59.46 50.78 53.60 54.89 51.52 66.21 58.15 56.53 56.73 51.56 55.60 53.61 54.92 49.91 57.03

55.32 52.53 49.92 65.22 63.27 62.40 64.18 63.30 58.19 66.99 64.87 57.78 59.22 59.13 70.72 66.18 63.03 62.78 59.77 56.67 54.04 51.28 48.16 58.19 56.57 60.58 52.83 56.33 57.98 53.92 73.11 59.34 57.50 60.66 50.96 55.62 56.91 59.39 52.61 62.24

96.44 100.61 142.73 100.70 100.23 105.44 98.32 98.73 112.67 98.89 102.18 101.47 104.99 96.38 102.65 99.48 101.60 100.74 92.28 98.16 95.53 152.11 116.19 105.56 105.80 97.77 101.94 108.28 98.88 101.00 98.14 91.00 97.04 93.45 99.99 101.70 99.95 98.70 110.20 94.02

105.08 114.05 165.15 110.97 109.06 110.37 108.75 109.17 140.77 110.43 111.10 116.45 116.58 106.99 112.10 111.73 114.33 114.63 106.98 111.56 106.90 212.70 151.91 127.23 116.57 106.73 113.40 124.37 111.59 114.53 105.66 100.88 108.21 105.45 113.47 115.40 110.52 114.84 121.29 105.24

M

ED

PT CE AC

AN US

Instance

21

ACCEPTED MANUSCRIPT

Table 5: D-Metric values of the solutions from j120 with the same running time

NGSA-II EM

NSGA-II

SPEA2 EM

SPEA2

j1201 1 j1203 1 j1203 4 j1207 1 j1207 2 j1207 3 j1208 2 j1208 4 j1209 1 j12013 1 j12013 3 j12013 5 j12014 1 j12014 4 j12017 1 j12017 4 j12018 1 j12018 2 j12019 4 j12019 6 j12020 1 j12021 6 j12022 6 j12028 2 j12029 5 j12033 1 j12034 4 j12034 6 j12035 2 j12035 3 j12037 6 j12038 1 j12039 1 j12039 2 j12040 1 j12040 2 j12049 3 j12053 3 j12055 4 j12058 6

0.1811 0.1004 0.0926 0.2264 0.1740 0.1179 0.2782 0.2758 1.6862 0.3461 0.3315 0.6778 0.4556 0.2811 0.3680 0.3659 0.9113 0.7030 0.2996 0.2901 0.3298 0.0196 0.0381 0.2182 1.5082 0.2214 0.2078 0.1709 0.2197 0.2696 0.1854 0.2663 0.1808 0.1802 0.4996 0.2797 0.1629 0.2835 0.4114 0.2253

0.2147 0.2228 0.0244 0.3447 0.2378 0.1676 0.3806 0.3518 1.7419 0.2897 0.3951 0.7182 0.4797 0.2989 0.3961 0.3827 0.9026 0.9177 0.3328 0.2742 0.3571 0.0020 0.0216 0.1863 1.4806 0.2160 0.2454 0.1291 0.2176 0.3095 0.2261 0.3011 0.1958 0.2150 0.5528 0.3031 0.1634 0.2887 0.3504 0.2638

0.1659 0.1573 0.2264 0.2696 0.2194 0.1413 0.3648 0.3384 2.3424 0.3739 0.4667 0.8311 0.6389 0.3052 0.4224 0.4128 0.8568 0.8588 0.3433 0.3646 0.3930 0.0242 0.0420 0.2753 1.9834 0.2543 0.2908 0.1843 0.2682 0.3025 0.2162 0.2735 0.1935 0.2146 0.7534 0.4197 0.2161 0.3935 0.4799 0.2567

0.1621 0.1573 0.1114 0.3458 0.2689 0.1420 0.3906 0.3188 2.0191 0.3370 0.4506 0.8437 0.5836 0.3287 0.4476 0.4341 0.8235 0.9041 0.3491 0.3500 0.3918 0.0082 0.0468 0.2722 1.8186 0.2234 0.3191 0.1650 0.2539 0.3047 0.2431 0.3006 0.2119 0.2125 0.5975 0.3933 0.2295 0.4188 0.5037 0.2868

AN US

M

ED

PT

CR IP T

Instance

AC

CE

From the experiment results showed above, obviously, in j90,j120 set, combing EM heuristics with NSGA-II and SPEA2 can receive performance improvement, especially for NSGA-II. NSGAII EM outperforms NSGA-II in terms of solution quality and convergence, so does SPEA2 EM. Comparing with MOEA/D, MOEA/D EM does not have a better effect. The reason why EM plays different roles in these three algorithms is derived from the different individual selection strategy for crossover operator. In Step 8 of NSGA-II, the whole population is used for selection and crossover. For SPEA2, the mating pool used for recombination and mutation is constructed with the individuals selected from the archive population. According to the principle of EM, the larger the difference between two individuals used for crossover is, the more obvious the effect is. Because EM will guide the search direction in terms of the objective values, EM can hardly work when the two individuals used for crossover have similar objective values. In NSGA-II, all individuals in population take part in crossover operator, so the diversity and differences in individuals can be guaranteed for EM. SPEA2 only selects individuals from archive population for crossover, so EM 22

ACCEPTED MANUSCRIPT

C(NSGA-II_EM, NSGA-II)

C(SPEA2_EM, SPEA2)

C(NSGA-II, NSGA-II_EM)

%

100

100

60

40

20

80

60

40

20

0 0

4

8

12

16

20

24

28

32

36

40

44

0

0

4

CR IP T

80 Average C-metric value

Average C-metric value

C(SPEA2, SPEA2_EM)

%

8

12

16

20

24

28

32

36

40

Instance

Instance

AN US

Figure 14: Average set coverage between algorithms integrated with EM and their original versions for j120 set with the same running time

ED

M

heuristics integrated into SPEA2 can not work as well as in NSGA-II due to the fact that individuals in archive are elitist and the similarities among them reduce the effectiveness of the EM heuristics. Nevertheless, the approach to generate new solutions in MOEA/D is quiet different from NSGA-II and SPEA2. In Step 2.1 of MOEA/D, a new solution is generated from two neighbors which are similar by using crossover and mutation operator. Furthermore, updating the neighborhood with this new solution will make the neighbor individuals more similar in Step 2.3. So it may be the important reason why EM could not improve the performance of MOEA/D. To sum up, NSGA-II EM is probably the most effective algorithm for MORCPSP among the algorithms we studied.

PT

8. Conclusions

AC

CE

In this paper, we have studied if the electromagnetism (EM) heuristic procedure could be integrated into multi-objective evolutionary algorithms (MOEAs) for the multi-objective resourceconstrained project scheduling problem (MORCPSP), one of the most challenging combinatorial optimization problems in scheduling. Not only the approach of extending EM heuristics to multiobjective optimization but also how to combine it with MOEAs framework for MORCPSP have been introduced. We have conducted experiments on integrating EM heuristics with three state-ofthe-art MOEAs including NSGA-II, SPEA2 and MOEA/D, and compared with their corresponding original versions for 160 standard test instances. The computational results show that integrating EM heuristics with NSGA-II and SPEA2, especially for NSGA-II, could improve various aspects of performance, and it may not suitable to combine EM with MOEA/D. We also give a brief explanation that why EM plays different roles in these three algorithms. Finally , we can get a conclusion that NSGA-II integrated with EM heuristics may be the most effective algorithm for MORCPSP among the algorithms we studied. The further research will consist of doing research on the other objectives mentioned in this paper and developing exact algorithms capable of calculating the real Pareto front for MORCPSP. 23

44

ACCEPTED MANUSCRIPT

Acknowledgments We would like to thank the anonymous reviewers to improve the quality of this paper. This work was partially supported by the National Natural Science Foundation of China (NSFC) projects No.61202296, 71102146 and 61272067, the National High-Technology Research and Development Program (”863”Program) of China under Grand No.2013AA01A212, the Natural Science Foundation of Guangdong Province project No.S2012030006242.

CR IP T

References

AC

CE

PT

ED

M

AN US

Abbasi, B., Shadrokh, S., Arkat, J., 2006. Bi-objective resource-constrained project scheduling with robustness and makespan criteria. Applied Mathematics and Computation 180 (1), 146–152. Aboutalebi, R., Najafi, A., Ghorashi, B., 2012. Solving multi-mode resource-constrained project scheduling problem using two multi-objective evolutionary algorithms. Afr J Bus Manage 6 (11), 4057–65. Al-Fawzan, M. A., Haouari, M., 2005. A bi-objective model for robust resource-constrained project scheduling. International Journal of Production Economics 96 (2), 175–187. Ballest´ın, F., Blanco, R., 2011. Theoretical and practical fundamentals for multi-objective optimisation in resourceconstrained project scheduling problems. Computers & Operations Research 38 (1), 51–62. Ballest´ın, F., Valls, V., Quintanilla, S., 2006. Due dates and rcpsp. In: J´ ozefowska, J., W¸eglarz, J. (Eds.), Perspectives in Modern Project Scheduling. Vol. 92 of International Series in Operations Research & Management Science. Springer US, Ch. 4, pp. 79–104. ˙ Fang, S.-C., 2003. An electromagnetism-like mechanism for global optimization. Journal of global Birbil, S ¸ . I., optimization 25 (3), 263–282. Brucker, P., Drexl, A., M¨ ohring, R., Neumann, K., Pesch, E., 1999. Resource-constrained project scheduling: Notation, classification, models, and methods. European Journal of Operational Research 112 (1), 3–41. Corne, D. W., Jerram, N. R., Knowles, J. D., Oates, M. J., J, M., 2001. Pesa-ii: Region-based selection in evolutionary multiobjective optimization. In: Proceedings of the Genetic and Evolutionary Computation Conference (GECCO2001). Morgan Kaufmann Publishers, pp. 283–290. Corne, D. W., Knowles, J. D., Oates, M. J., 2000. The pareto envelope-based selection algorithm for multiobjective optimization. In: Parallel Problem Solving from Nature PPSN VI. Springer, pp. 839–848. Deb, K., Pratap, A., Agarwal, S., Meyarivan, T., 2002. A fast and elitist multiobjective genetic algorithm: Nsga-ii. Evolutionary Computation, IEEE Transactions on 6 (2), 182–197. Debels, D., De Reyck, B., Leus, R., Vanhoucke, M., 2006. A hybrid scatter search/electromagnetism meta-heuristic for project scheduling. European Journal of Operational Research 169 (2), 638–653. Erickson, M., Mayer, A., Horn, J., 2001. The niched pareto genetic algorithm 2 applied to the design of groundwater remediation systems. In: Zitzler, E., Thiele, L., Deb, K., Coello Coello, C., Corne, D. (Eds.), Evolutionary MultiCriterion Optimization. Vol. 1993 of Lecture Notes in Computer Science. Springer Berlin Heidelberg, Ch. 48, pp. 681–695. Fonseca, C. M., Fleming, P. J., et al., 1993. Genetic algorithms for multiobjective optimization: Formulation, discussion and generalization. In: ICGA. Vol. 93. pp. 416–423. Gomes, H. C., de Assis das Neves, F., Souza, M. J. F., 2014. Multi-objective metaheuristic algorithms for the resource-constrained project scheduling problem with precedence relations. Computers & Operations Research 44 (0), 92 – 104. Hapke, M., Jaszkiewicz, A., Slowinski, R., 1997. Fuzzy project scheduling with multiple criteria. In: Fuzzy Systems, 1997., Proceedings of the Sixth IEEE International Conference on. Vol. 3. IEEE, pp. 1277–1282. Hapke, M., Jaszkiewicz, A., Slowi´ nski, R., 1998. Interactive analysis of multiple-criteria project scheduling problems. European Journal of Operational Research 107 (2), 315–324. Hapke, M., Jaszkiewicz, A., Slowi´ nski, R., 1999. Fuzzy multi-mode resource-constrained project scheduling with multiple objectives. In: W¸eglarz, J. (Ed.), Project Scheduling. Vol. 14 of International Series in Operations Research & Management Science. Springer US, Ch. 16, pp. 353–380. Herroelen, W., Demeulemeester, E., Reyck, B., 1999. A classification scheme for project scheduling. In: W¸eglarz, J. (Ed.), Project Scheduling. Vol. 14 of International Series in Operations Research & Management Science. Springer US, Ch. 1, pp. 1–26.

24

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

Horn, J., Nafpliotis, N., Goldberg, D. E., 1994. A niched pareto genetic algorithm for multiobjective optimization. In: Evolutionary Computation, 1994. IEEE World Congress on Computational Intelligence., Proceedings of the First IEEE Conference on. IEEE, pp. 82–87. Icmeli, O., Erenguc, S. S., Zappe, C. J., 1993. Project scheduling problems: A survey. International Journal of Operations & Production Management 13 (11), 80–91. Kazemi, F., Tavakkoli-Moghaddan, R., 2008. Solving a multi-objective multi-mode resource-constrained project scheduling with discounted cash flows. In: Proceedings of the 6th international management conference. Knowles, J. D., Corne, D. W., 2000. Approximating the nondominated front using the pareto archived evolution strategy. Evolutionary computation 8 (2), 149–172. Kolisch, R., Padman, R., 2001. An integrated survey of deterministic project scheduling. Omega 29 (3), 249–272. Kolisch, R., Sprecher, A., Drexl, A., 1995. Characterization and generation of a general class of resource-constrained project scheduling problems. Management science 41 (10), 1693–1703. Li, K., Willis, R., 1992. An iterative scheduling technique for resource-constrained project scheduling. European Journal of Operational Research 56 (3), 370–379. Mendes, J., Goncalves, J., Resende, M., 2009. A random key based genetic algorithm for the resource constrained project scheduling problem. Computers & Operations Research 36 (1), 92–109. Miettinen, K., 1999. Nonlinear multiobjective optimization. Vol. 12 of International Series in Operations Research & Management Science. Kluwer Academic. Nabrzyski, J., W¸eglarz, J., 1995. On an expert system with tabu search for multiobjective project scheduling. In: Emerging Technologies and Factory Automation, 1995. ETFA ’95, Proceedings., 1995 INRIA/IEEE Symposium on. Vol. 3. pp. 87–94. Nabrzyski, J., W¸eglarz, J., 1999. Knowledge-based multiobjective project scheduling problems. In: W¸eglarz, J. (Ed.), Project Scheduling. Vol. 14 of International Series in Operations Research & Management Science. Springer US, Ch. 17, pp. 383–411. ¨ OZdamar, L., Ulusoy, G., 1995. A survey on the resource-constrained project scheduling problem. IIE Transactions 27 (5), 574–586. ¨ Ozdamar, L., Ulusoy, G., 1996. A note on an iterative forward/backward scheduling technique with reference to a procedure by li and willis. European Journal of Operational Research 89 (2), 400–407. Pan, H., Yeh, C.-H., 2003. A metaheuristic approach to fuzzy project scheduling. In: Palade, V., Howlett, R., Jain, L. (Eds.), Knowledge-Based Intelligent Information and Engineering Systems. Vol. 2773 of Lecture Notes in Computer Science. Springer Berlin Heidelberg, Ch. 145, pp. 1081–1087. Slowi´ nski, R., 1981. Multiobjective network scheduling with efficient use of renewable and nonrenewable resources. European Journal of Operational Research 7 (3), 265–273. Slowi´ nski, R., Soniewicki, B., W¸eglarz, J., 1994. Dss for multiobjective project scheduling. European Journal of Operational Research 79 (2), 220–229. Srinivas, N., Deb, K., 1994. Multiobjective optimization using nondominated sorting in genetic algorithms. Evolutionary computation 2 (3), 221–248. Viana, A., Pinho de Sousa, J., 2000. Using metaheuristics in multiobjective resource constrained project scheduling. European Journal of Operational Research 120 (2), 359–374. Wang, L., Fang, C., Mu, C.-D., Liu, M., Aug 2013. A pareto-archived estimation-of-distribution algorithm for multiobjective resource-constrained project scheduling problem. Engineering Management, IEEE Transactions on 60 (3), 617–626. Zhang, Q., Li, H., 2007. Moea/d: A multiobjective evolutionary algorithm based on decomposition. Evolutionary Computation, IEEE Transactions on 11 (6), 712–731. Zhou, A., Zhang, Q., Jin, Y., Sendhoff, B., Tsang, E., 2007. Global multiobjective optimization via estimation of distribution algorithm with biased initialization and crossover. In: Proceedings of the 9th Annual Conference on Genetic and Evolutionary Computation. GECCO ’07. ACM, pp. 617–623. Zitzler, E., Laumanns, M., Thiele, L., 2001. Spea2: Improving the strength pareto evolutionary algorithm. Technical Report 103, Computer Engineering and Networks Laboratory (TIK), (ETH) Zurich, Switzerland. Zitzler, E., Thiele, L., 1999. Multiobjective evolutionary algorithms: A comparative case study and the strength pareto approach. Evolutionary Computation, IEEE Transactions on 3 (4), 257–271.

25