Non-dominated Sorting Genetic Algorithms for a Multi-objective Resource Constraint Project Scheduling Problem

1.1 Problem Statement

The resource constraint project scheduling problem (RCPSP) was proposed in 1969 by Pritsker et al. [30] for the first time. It was later proved to be NP-hard in the strong sense in 1983 by Blazewicz et al. [6]. In the RCPSP, we consider a scheduling problem with a set of non-preemptive jobs, J={0, 1, …, n, n+1}, submitted to a set of precedence relationships A. The jobs 0 and n+1 are dummy activities that represent the start and end of the project, respectively. The precedence graph, denoted by G=(J, A), can thus be defined with J and A. The jobs in process need given amounts of resources with limited capacities. A resource is called “renewable” if its capacity is brought to the maximum at the beginning of each period. Otherwise, a resource that is invested only once or progressively during the whole project is called “non-renewable.” The resource consumption cannot exceed the available resource quantities during each period of time. The basic RCPSP aims to minimize the makespan with respect to the precedence relationships and resource constraints.

Our earlier work is devoted to solve a bi-objective RCPSP [33]. In this paper, our problem is extended to an RCPSP with three criteria, namely the makespan, total tardiness of jobs, and workload balance. Moreover, in this work, we pay attention to resource allocation, which is a major problem that industry has to cope with. The non-dominated sorting genetic algorithms (NSGAII and NSGAIII) are firstly applied to solve our problem. Hence, instead of considering the Pareto dominance in the selection procedure, we integrate the Lorenz dominance and solve our problem with L-NSGAII [15] and then propose L-NSGAIII.

This paper is organized as follows: in Subsection 1.2, we present a specific literature review. Our three-criterion problem is presented in Section 2, in which we present a mathematical formulation. The NSGAII and NSGAIII are then applied to solve the multi-objective problem, as explained in Section 3. In Section 4, the principle of Lorenz dominance is explained as well as its application on the previously methods. In Section 5, we report the experimental results. Conclusion and perspectives are discussed in Section 6.

1.2 Literature Review

The RCPSP is firstly formulated with binary variables x_it that are equal to 1 if the job i starts at the moment t, and 0 otherwise, in the work of Pritsker et al. [30]. The completion time of a job i∈J can thus be written as ∑_0≤t≤Htx_it+p_i, where H is the scheduling horizon. Later, Christofides et al. [8] proposed a similar model by reformulating the precedence constraints. The work of Klein [22] defined two binary variables in order to describe the job status: s_it and f_it that take the value 1 if at time 0≤t≤H, the job i∈J has already started and finished, respectively. In Ref. [25], Mingozzi et al. proposed a formulation based on the idea of building sets of jobs that can be processed at the same time without violating the resource and precedence constraints. This formulation is proved to provide strong lower bounds on the makespan. The compatible job sets were considered by Alvarez-Valdés and Tamarit [3], and the authors proposed a mixed-integer linear programming formulation. Some other formulations are recently presented with innovative visions. Koné et al. [23] proposed an event-based formulation where an event can be the start/end of a job. The model of Bianco and Camaria [5] was based on the idea of job completion percentage: a job is only completed when the ratio reaches 100%. Other than the mathematical models, branch-and-bound methods were also applied by Brucker et al. [7], Christofides et al. [8], and Demeulemeester and Herroelen [14], to solve the RCPSP optimally. Some of these methods are known to be especially efficient, such as the work of Klein [22] and Brucker et al. [7]. However, due to the NP-hardness of the problem, exact methods only allow to solve instances with very limited number of jobs.

Meta-heuristics are also implemented by researchers: the genetic algorithm was applied by Debels and Vanhoucke [13] and Peteghem and Vanhoucke [29]. In Ref. [4], the authors used the Tabu search. Sirdey et al. [31] tackled the problems with the simulated annealing algorithm. In the work by Jia and Seo [19], a permutation-based artificial bee colony algorithm is chosen to find the solutions. A comparative study, including the methods mentioned above and particle swarm optimization, was carried out by Das and Acharyya [10].

However, unlike the RCPSP with a single criterion, the multi-objective RCPSP is not yet well documented. Only few works have been proposed to solve the multi-objective RCPSP with exact methods. A branch-and-bound-based procedure was proposed by Gutjahr [18] to solve a bi-criterion stochastic multi-mode RCPSP (MMRCPSP) under risk aversion, where the project makespan and the cost are to be optimized. We have considered the project makespan and the total job lateness in a bi-objective RCPSP, with the two-phase method (TPM) to find the optimal Pareto front [33].

More works have been proposed in the literature as approximated solutions. In Ref. [21], Kim et al. proposed a hybrid genetic algorithm with fuzzy logic controller in order to minimize the makespan and the total tardiness penalty. Simulated annealing and Tabu search were discussed in the works of Abbasi et al. [1] and Al-Fawzan and Haouari [2], respectively, where makespan and robustness were considered. The works of Vanucci et al. [32] focused on the bi-objective problem with makespan and total cost minimization in the case of MMRCPSP with NSGAII. The NSGAII can also be found in the work of Damak et al. [9], where an MMRCPSP was studied in order to minimize the makespan and the non-renewable resource cost. Khalili et al. [20] applied two genetic algorithms, namely multi-population genetic algorithm and two-phase sub-population genetic algorithm, to minimize the makespan and maximize the net present value. In the work by Gomes et al. [17], the minimization of makespan and total weighted job starting time was tackled and the efficiency of five meta-heuristic algorithms was compared. These methods were based on the multi-objective greedy randomized adaptive search procedure, the multi-objective variable neighborhood search, and the Pareto iterative local search.

Our works aim at solving the RCPSP in the multi-objective context. Our previous work [33] allowed finding the optimal Pareto front with the TPM for a bi-objective problem. In this paper, our work is extended by adding a third criterion, namely the workload balance level. The NSGAII is adapted to solve the three-criterion problem and compared to the recently proposed NSGAIII by Deb and Jain [11]. In addition to the Pareto niching mechanism, we also propose the L-NSGAII and L-NSGAIII for our problem by integrating the Lorenz dominance. For the sake of brevity, the makespan, total job tardiness, and workload balance are denoted by O₁, O₂, and O₃, respectively.

3.1 Chromosome and Initial Population

The initial population is constituted by N individuals, each of whom is identified by its chromosomes. Two chromosome parts are considered in our problem: a priority sequence of jobs, denoted by X₁, and a series of matching time lags, denoted by X₂. Regarding the example of Table 2, a whole chromosome is proposed in Figure 1 and denoted by χ.

Figure 1:

A Complete Chromosome with Nine Jobs.

For a given problem, X₁ is only valid if all the precedence relationships of jobs are respected. For instance, if a job i∈J is a direct or indirect predecessor of a job j∈J, then j cannot be ahead of i in the priority sequence.

For this purpose, we identify two sets of jobs: the available jobs (AJ) and the non-available jobs (NJ). At first, only job 0 is available so that AJ={0} and NJ=J\{0}. We choose the only job in AJ and put it in the first place of X₁. Next, considering that job 0 finished, we update the two job sets AJ and NJ: job 0 is firstly moved out of AJ; in the meantime, jobs whose predecessors are all finished are moved from NJ to AJ. In the example, jobs 1, 2, and 3 take job 0 as the only predecessor. Thus, after choosing job 0, they become available and we have AJ={1, 2, 3}. One job in AJ is then randomly chosen, for example, job 1. Next, we move job 1 out of AJ and move job 4 from NJ to AJ, as its only predecessor (job 1) is already chosen. We can thus update AJ={2, 3, 4}. Repeating this procedure explained above, we can obtain a valid sequence X₁, leading to a feasible solution.

Let us now focus on X₂, which is created specially for the workload balance criterion. This objective is measured by the gap between the highest and the lowest resource utilization. It is thus interesting for some jobs to wait before starting the operation so as to lower the eventual peak of resource consumptions. To cope with this problem, we introduce the second part of the chromosome, which contains a set of “non-negative time lags” between the actual and the earliest possible starting dates of jobs. For a given individual, if there exists at least one positive time lag, X₂ is said to be effective.

For a given chromosome, we consider each job one by one according to the priority sequence and assign a start time for each of them. A job i∈J can only start if three conditions are satisfied: firstly, all predecessors of i should already be finished; secondly, the time lag must be respected; and finally, there should be enough resources during periods when i is in process. With all conditions satisfied, the job starts at the earliest possible moment.

3.2 From Chromosome to Solution

Let us consider one single resource with periodic capacity q₁=5. Job durations and resource consumptions are defined in Table 2. In Figure 2, we provide the schedule obtained by X₁ of χ, where all jobs start as soon as possible and there is no time lag. This leads to a makespan of 6 and a workload balance of 3. Now let us add X₂ to complete the chromosome and assume that job 2 has a time lag that equals to 1. This means that job 2 cannot start earlier than t=1, and we find a schedule as shown in Figure 3, which proposes a higher makespan of 7 but a better workload balance of 2.

Figure 2:

Example of Solution.

Figure 3:

Schedule with TL₂=1.

Besides, it is worth noting that the priority sequence does not have to be the final schedule sequence due to the resource constraints. In the example of Figure 2, job 4 has the priority over job 3 in the chromosome but somehow starts later. As a matter of fact, given two jobs (i, j)∈J², if there is not any direct or indirect precedence relationship between them, even when i has the priority over j in X₁, it is possible that j starts before i. This case arises typically when the resource requirements of i and j are extremely high and low, respectively. Job i can only start when the resource availability is high; on the contrary, it is a lot easier for j to find a position in the partial schedule. As a result, it is possible to have two different chromosomes that lead to the same schedule, called “equivalent chromosomes.” In order to guarantee the diversity of solutions, equivalent chromosomes are not allowed in this step.

3.3 Crossover

Crossover allows creating diversified individuals. Both one-point and two-point crossovers are used in our program so as to explore more possibilities of chromosomes. The crossover points, denoted by CP₁ and CP₂, are randomly chosen among positions in [1, n−1], so that the chromosome can be divided into several parts. It is worth noting that CP₁ is by default in front of CP₂. Let us remember the example in Table 2 and denote FX₁, FX₂ as the chromosomes of the father and MX₁, MX₂ as those of the mother. As shown in Figure 4, the son, denoted by SX₁, SX₂, completely inherits the father’s chromosomes before CP₁ and after CP₂ (for one-point crossover, CP₂ is automatically set to the last position). For the remaining part, the son follows the mother’s sequence for the remaining jobs. However, considering the precedence relationships, a simple combination of partial chromosomes may lead to an invalid chromosome. To guarantee a feasible schedule in SX₁, for the positions between CP₁ and CP₂, jobs in FX₁ have to follow the same relative sequence in MX₁. In the meantime, it is important to keep the “pairing” relationship between X₁ and X₂. A job and its time lag should always be on the same position. The daughters (DX₁, DX₂) are created in the same way by reversing the roles of the father and the mother.

Figure 4:

Two-Point Crossover.

3.4 Mutation

Mutation allows furthering the diversity of chromosomes and can help step out of possible local optimum. Considering X₁, for a given individual, we choose firstly a mutation job, denoted by j_m. A time window, defined as between the end of the last predecessor and the beginning of the first successor, frames a free zone where j_m can start at any moment without violating the precedence relationships. We then choose a reference job randomly in this window and denote it as j_r. If j_m is before (after) j_r, then we put j_m right after (before) j_r and the rest of the jobs keep their relative order in the origin chromosome. Once again, we should pay attention to the pairing relationship between X₁ and X₂ and make sure that X₂ keeps us with X₁ in the mutation. Moreover, if the time lag of j_m was positive before the mutation, we set it to 0, and otherwise, we assign a positive time lag for it. In Figure 5, X₁, X₂ represents the original chromosome while X₁ʹ, X₂ʹ stands for that after the mutation.

Figure 5:

Mutation.

3.5 Selection

After the offspring, the children and parents are gathered into the same pool from which N individuals with be selected according to the Pareto dominance. Let us consider a problem with K objectives to minimize. For a given solutions s, denote O_k(s) its objective values for k=1, 2, …, K. For two solutions s₁ and s₂, if ∀k=1, 2, …, K, O_k(s₁)≤O_k and ∃i∈1, 2, …, K|O_i(s₁)<O_i(s₂), s₂ is said to be Pareto dominated by s₁ (s₁ π _P s₂). Hence, the two solutions s₁ and s₂ are mutually non-dominated if and only if s₁≡_P s₁ and s₂≡_P s₁. A set of mutually non-dominated solutions is on the same front.

During the selection, non-dominated fronts are established. Let us denote these fronts by {F₁, F₂, …, F_last}, where F₁ contains the best solutions while F_last is constituted by the worst ones in the given population. In most cases, F_last requires a selection within the non-dominated front. In this case, the NSGAII uses the crowding distance, denoted by D_c, to measure the solution distribution on the solution space. Solutions with higher values of D_c are selected as they are located on less visited zones. Readers can refer to Ref. [12] for more details.

As stated earlier, we have to deal with the equivalent chromosomes in our problem. For this purpose, we limit the size of each front F_i (i=1, 2, …) by a parameter δ, such that the number of solutions chosen on F_i should not exceed δ*N. This means that the selection within a front is not only proceeded for F_last, but may also be necessary for any other front. Generally speaking, the extreme solutions have the highest priority so that the exploration zone is as large as possible. However, all extreme solutions are not necessarily chosen. When lots of equivalent solutions are identified, choosing simply those with the highest D_c can limit the diversity of the population as well as the solution distribution on the search space. As a result, the equivalent solutions are gathered together as a group. The selection always follows the principle of crowding distance; however, for a group of equivalent solutions, only one is chosen at a time. Once all the groups are visited, we return to the head of the list and repeat the previous step until we have enough individuals.

3.6 Stopping Criteria and Equivalent Chromosomes

Two stopping criteria are considered in our program. If the first front does not vary for ϕ generations, we conclude the convergence of the population and the program is stopped. Otherwise, the program ends when the number of generation attains θ and we keep the last-found generation.

4.1 Definition

In a multi-objective problem with K criteria to minimize, given two solutions s₁ and s₂, s₂ is said to be Lorenz dominated by s₁ if the Lorenz vector of s₂ is Pareto dominated by that of s₁, denoted by s₁π_L s₂. In order to find the Lorenz vector of a solution s, we need to follow several steps. Firstly, the objective values are normalized to a value in [0, 1] by

where max O_k and min O_k are, respectively, the maximum and minimum values on the objective k∈1, …, K in the actual generation. The normalized objectives are then sorted in the decreasing order such that

The Lorenz vector can thus be found as

where OkL(s) stands for the k^th element in the Lorenz vector of s, computed by the cumulus of the first k elements in the sorted normalized vector. In Figure 6, we give an illustration for the bi-objective case. For a given solution s, the dominated area provided by L-dominance contains that provided by P-dominance. Besides, the L-dominated area is almost symmetric to the line O1L=O2L except for the symmetric point of s itself. This is due to the fact that values in the Lorenz vector is sorted in decreasing order. Finally, the Lorenz dominance can be established by comparing the Lorenz vectors according to the Pareto dominance: given two solution s_a and s_b, if the Lorenz vector of s_a dominates that of s_b in the Pareto sense, then s_a dominates s_b regarding the L-dominance.

Figure 6:

Pareto-dominated (Left) and Lorenz-dominated (Right) Area.

4.2 Numerical Example

For illustrative purposes, let us consider a bi-objective example with 10 solutions as in Table 3. The solutions are sorted with the Pareto dominance, as given in the column P-rank.

In order to implement the Lorenz dominance, the solutions are firstly normalized according to Eq. (16). Next, for each solution, their normalized objectives are sorted in decreasing order as described in Eq. (17). Thus, we can obtain the Lorenz vector OkL by cumulating the values of O′[k] as in Eq. (18). This transformation procedure is also illustrated in Figure 7. Finally, the solutions are re-ranked by comparing their Lorenz vectors according to the Pareto dominance.

Figure 7:

Numerical Example – Creation of the Lorenz Vectors from Original Solution Coordinates.

The L-dominance aims at giving priority to solutions that equally optimize each objective. In the example, s₅ is classified as rank 1 by the L-dominance as it proposes good performances on both O₁ and O₂. On the other hand, the solution s₁, which provides a good result on O₁ but a bad result on O₂, is poorly ranked by the L-dominance even if it belongs to the first Pareto front.

4.3 Adaptation on NSGAs

The Lorenz dominance is then integrated in our programs. The basic idea is to sort and select the solutions according to their Lorenz vectors instead of their original coordinates. For instance, at the beginning of the selection, we compute the Lorenz vectors before proceeding the non-dominated sorting. The Pareto fronts are thus built according to the Lorenz vectors of solutions. Moreover, while the selection within a non-dominated front is involved, the crowding distances (NSGAII), the reference points, and the hyperplane (NSGAIII) are calculated regarding the Lorenz vectors as well. Inheriting all the other procedures from the methods presented earlier, the L-NSGAII and L-NSGAIII are applied to solve our problem.

5.1 Data Generation and Experimental Design

As stated in our previous study [33], several parameters are used during the data generation in the literature, such as net complexity (NC) and resource strength. Inspired by these works and considering the proper characteristics of our problem, we proposed a new set of instances that allow evaluating our methods with adapted parameters, as described below:

• NC: ratio between the number of precedence relationships and the number of non-dummy jobs. A higher NC value leads to a more tightly constrained network with less possible schedules.

• Demand rate (DR): proportion of positive resource requirement q_ir with i∈J and r∈R.

• Demand factor (DF): parameter that allows generating resource requirements. For a given non-null requirement for i∈J and r∈R, q_ir is randomly chosen in [1, DF * q_r] according to uniform distribution.

• Delay proportion (DP): proportion of jobs whose due dates are set as their earliest finishing dates. The rest of the due dates are chosen randomly with uniform distribution in [EF_i+1, EF_i+p_i] for a given job i∈J.

With these parameters, we have generated instances with different combinations. Let us denote GY as the group of instances with Y non-dummy jobs.

As presented in previous sections, we solve this multi-objective problem with the NSGAII, L-NSGAII, NSGAIII, and L-NSGAIII (respectively denoted by M1, M2, M3, and M4). Instances with 30–150 jobs are considered with DP∈{0.3, 0.5, 0.7}, DR∈{0.4, 0.6, 0.8}, and DF∈{0.5, 0.7, 0.9}. For all the groups, we have set P_c=0.9, P_mut=0.05. For the sake of brevity, the subgroups are denoted by G30.γ.α.β, where γ, α, β indicate the position of the chosen value in the DP, DR, and DF sets, respectively. For instance, when DP=0.3, we have γ=0. Each subgroup contains 10 randomly created instances.

Our methods are evaluated and compared with two indicators – the hypervolume and the C-metric. The hypervolume is introduced by Zitzler et al. [34] and measures, for a Pareto front, the coverage level of its dominance zone. A value in (0, 1) is assigned to the evaluated front, where 1 stands for the perfect coverage and 0 represents extremely poor coverage. In the work of Fonseca et al. [16], the authors recently proposed a fast algorithm for computing the hypervolume and integrated the computation procedure especially for the three-objective case. We adapt this method as well as the program proposed by the authors to our problem. The C-metric allows comparing two Pareto fronts. It defines the proportion of solutions on one front that are dominated by at least one solution on the other front. A lower value for the C-metric stands for better performances.

5.2 Experimental Results

We observe firstly the performance of the four methods small instance G30 with NC=1. The population size is N=100 with the stopping criteria parameters θ=500 and ϕ=50. The experimental results are presented in Table 4, where the better results are highlighted in bold letters. As the computation times of the four methods are very close, they are not included in the table. All methods require 0.06 s on average for each generation.

Regarding the hypervolume, we can see that the NSGAII-based methods (M1 and M2) have better performances that the NSGAIII-based methods (M3 and M4). The NSGAII provides better solutions in 14 out of 27 subgroups, and is the best among the four methods. The NSGAIII-based methods give higher hypervolumes in only six groups among the 27 tested. These instances are characterized by medium resource requirements. In Ref. [33], we have shown that difficulty of instances arises with higher resource requirements. Furthermore, according to this metric, the Lorenz dominance does not seem to have strong improvements for the Pareto dominance on G30.

With the C-metric, the methods are studied in pairs. According to the comparison between M1 and M2, as well as that between M3 and M4, it is obvious that the Lorenz dominance improves the solution quality both for the NSGAII and NSGAIII approaches. Moreover, this improvement is stronger for the NSGAIII. By comparing M1 and M3, we find the same tendency as for the hypervolume, which means that the NSGAII outperforms the NSGAIII in most cases, especially in the easy and difficult instances. Finally, the comparison between M2 and M4 supports the previous statement that the Lorenz dominance has a much stronger improvement for the NSGAIII than the NSGAII.

Combining the two metrics, we find that the NSGAII holds its superior performance over the NSGAIII on easy and difficult instances, for both hypervolume and C-metric. Moreover, the Lorenz-based methods seem to be weaker. Nevertheless, while considering the C-metric, the improvements of Lorenz dominance are revealed. Even if it does not provide a higher hypervolume, it somehow generates more “truly” non-dominated solutions.

Moving on to medium-sized problems, for instance, the group G60 with NC=1.33, whose results are shown in Table 5. Once again, the bold letters represent the better results. In these experiments, we have set N=125, θ=700, and ϕ=70. The computation times of different methods are still very close, which is around 0.20 s. Compared to G30, NSGAIII-based methods start to show their advantages. Considering the hypervolume, M1 and M2 still provide better solutions for more instances; however, the proportion has fallen from 24 to 16 subgroups out of 27. In the meantime, the advantages of NSGAII are, on easy and difficult instances, reduced as well.

This tendency is proved by the comparison between M1, M3 and M2, M4 regarding the C-metric. It is shown that NSGAIII-based methods propose more interesting solutions. However, the improvements of Lorenz dominance are not as strong as in G30. The performance of Lorenz-based methods is slightly weaker than the Pareto-based methods. Within both the pairs M1/M2 and M3/M4, the better performances are rather evenly located and their behaviors on instances with moderate resource requirements are relatively close. However, it is worth noting that M4 still proposes better solutions on most instances than M2, while M3 is more efficient than M1 in around half of the cases. This meets our earlier remark that the Lorenz dominance has a stronger effect on NSGAIII than NSGAII.

Finally, we present the results on large instances by taking the example of G120. In Table 6 with NC=2, N=120, θ=1000, and ϕ=100 and better results marked in bold letters. The computation time for each generation varies from 0.45 s (M1) to 0.53 s (M4). The NSGAIII-based methods continue to outperform the NSGAII-based methods in 19 out of 27 subgroups on the hypervolume. The latter loses its superiority on both easy and difficult instances. The efficiency of NSGAIII-based methods is proved once more by comparing the C-metric of M1 and M3, as well as that of M2 and M4. Meanwhile, the improvements of Lorenz dominance rewind regarding the pairs M1, M2 and M3, M4. In Figure 8, where we provide an example comparing the NSGAIII and the L-NSGAIII, the advantage of Lorenz dominance is clearly shown. We also notice that the C-metrics of M2 and M4 are lower than those of M1 and M3, and the two Pareto fronts found by Lorenz-based methods are closer to each other than the other two methods. This means that the Lorenz dominance allows finding more stable solutions than the Pareto dominance.

Figure 8:

Example of Comparison Between NSGAIII and L-NSGAIII.

By observing the behaviors of different groups of instances, we notice that for small instances, NSGAII is the most robust method and the Lorenz dominance does not provide significant improvements. However, as the problem becomes more complicated, the NSGAIII-based methods show their advantage and take over the superiority. The Lorenz dominance always has an obvious effect of improvements considering the C-metric as it tends toward the ideal point, and this effect is stronger on NSGAIII than on NSGAII.