Journal Search Engine
Search Advanced Search Adode Reader(link)
Download PDF Export Citaion korean bibliography PMC previewer
ISSN : 1598-7248 (Print)
ISSN : 2234-6473 (Online)
Industrial Engineering & Management Systems Vol.16 No.3 pp.288-306

Solving a Multi-Objective Mathematical Model for a Multi-Skilled Project Scheduling Problem by Particle Swarm Optimization and Differential Evolution Algorithms

Parisa Shahnazari-Shahrezaei, Sina Zabihi, Reza Kia*
Department of Industrial Engineering, Firoozkooh Branch, Islamic Azad University, Firoozkooh, Iran
Corresponding Author,
November 22, 2016 April 29, 2017 May 8, 2017


A Multi-Skilled Project Scheduling Problem (MSPSP) that is an extension of a Multi-Mode Resource-Constrained Project Scheduling Problem (MM-RCPSP) has been generally addressed to schedule a project with staff members as resources. In MSPSP, each activity requires different specialties and each staff member has a known skill level in performing an activity. This causes to encounter a huge number of modes while performing activities of a project. This research focuses on a special type of MSPSP known as Multi-Objective Multi-Skilled Project Scheduling Problem (MOMSPSP) which incorporates some new objectives in the MSPSP and develops a multi-objective mixed-integer nonlinear programming (MINLP) model. The model is exactly solved for small-sized instances using CPLEX solver. To solve such a NP-hard problem for medium and large-sized instances, two efficient meta-heuristic algorithms based on Differential Evolution (DE) and Particle Swarm Optimization (PSO) are proposed. To evaluate the efficiency of the proposed algorithms, the results are compared with each other as well as to the optimal ones obtained by the CPLEX solver for small instances. Finally, the designed DE algorithm is identified as the superior proposed algorithm for solving the propounded MOMSPSP in terms of some performance metrics.



    Multi-skilled Project Scheduling Problem (MSPSP) is a combination of the classic Resource-Constrained Project Scheduling Problem (RCPSP) and Multi-Purpose Machine (MPM) Hartmann and Briskorn (2010) developed first by Neron and Boptista (2002). Besides the existing constraints in the classic RCPSP problem, they added a novel constraint to the problem inspired by MPMs and called the new problem MSPSP. In this problem it is assumed that the resources are staff members which are engaged in the project so that each staff member has more than one speciality, each activity requires some specialties with certain skill level of each specialty to be performed and every staff member can be allocated to at most one activity simulatnously. Since each staff member has different skill levels, many modes of this problem to perform the activities can be considered. Re- garding so many modes of this problem, it is impossible to solve the model even in medium sizes by the traditional methods proposed for the Multi-Mode Resource Constrained Project Scheduling Problem (MMRCPSP) which are usually based on exhaustive enumeration methods (Neron, 2002).

    Bellenguez and Neron (2004) presented a mathematical model for MSPSP by promoting MSPSP to multiskill mode and multi-level human resources and used lower bound method for minimizing the total project time. Corominas et al. (2006) studied the multi-task staff member’s allocation in MSPSP and stated a heuristic approach based on a sequence of dependent allocation problems. Kazemipoor et al. (2013) solved MSPSP problem applying scatter search and Tabu Search (TS). In another study, Kazemipoor et al. (2012) solved a new mixed-integer model for MSPSP using simulated annealing algorithm. Also, Kazemipoor et al. (2013) suggested a novel goal programming model for the multi-objective multi-skill portfolio scheduling problem and employed Differential Evolution (DE) and TS algorithms.

    Kadrou and Najid (2006) presented a practical production scheduling based on a RCPSP considering each activity needs a special skill to be done to minimize the project completion time. Pessan et al. (2007) used a branch and bound algorithm to solve a mathematical model for MSPSP assuming that the organization's activities can fail.

    Heimerl and Kolisch (2010) proposed a mixedinteger linear programming model for IT projects scheduling and internal and external human force resources allocation to the projects. They evaluated the model performance, solution time, and the gap between the solutions resulted from the exact and heuristic methods. Al- Anzi et al. (2010) focused on RCPSP considering human resources with various skills. Santos and Tereso (2011) suggested a mathematical model for solving multi-mode multilevel multi-skill RCPSP. Firat and Hurkens (2012) presented a mathematical model for the project scheduling problem with several heterogeneous resources assuming that each of the resources has various skills and each skill has diverse levels. Liu and Wang (2012) proposed a mathematical model with combination of single-skill and multi-skill work force in order to solve the Project Scheduling Problem (PSP) and applied a combination of constrained programming (CP) approach and scheduling rules-based heuristics.

    Yoosefzadeh and Tareghian (2013) developed an efficient hybrid method depending on a branch and bound with a Genetic Algorithm (GA) coupled with a new schedule generation scheme for solving a RCPSP. Montoya et al. (2014) extended a mathematical model as a combination of two problems, i.e., RCPSP and MPM and developed a novel method combining column generation approach and Branch and Price algorithm. Li and Womer (2009) proposed a Bender’s decomposition algorithm for solving RCPSP with multi-skill resources to minimize the human force-related costs. Dhib et al. (2012) presented new lower bounds for multi-skill project scheduling solution. Myszkowski et al. (2015) developed an ant colony optimization for a multi-skill RCPSP. Jia and Seo (2013) proposed an improved Particle Swarm Optimization (PSO) algorithm for a RCPSP. Khalili et al. (2013) proposed two meta-heuristic algorithms of multi-population and twophase sub-population GAs to simultaneously minimize the project makespan and maximize the project net present value (NPV) of a multi-objective RCPSP.

    Having reviewed the conducted studies in MSPSP, it is revealed that only minimization of the project time and costs of recruiting human force in the project have been focused on, while in the practical project environments, financial issues are very critical and determining. Thus, the incorporation of significant costs including idle time cost of non-allocated resources and inefficient utilization costs of resources allocated to activities at lower skill levels is the main motivation of this paper to present a multi-objective model with the intention of minimizing different project costs while reducing the project time concurrently.

    The objectives considered in the proposed mathematical model are as the following:

    • 1) Minimizing the project completion time: In project scheduling problems, project time minimization is always considered as one of the objectives.

    • 2) Minimizing the resources idle time: The resources considered in this study are staff members and the project cost is related to them. They are working on the project as a full-time manner. Also, they are paid fully on a regular basis during the project. Therefore, the idleness of staff members during the project is undesirable. Hence, the second objective in the proposed model is to minimize the staff members idle time during the whole project. In fact, through considering this objective, staff members are allocated in the whole project in such a manner that provides the opportunity for all of them to participate in the project as much as possible, so that no member is left idle. Since each member possesses different specialties with known skill levels, they have different salary. Therefore, the second objective can be defined as to minimize the staff member’s idleness cost.

    • 3) Minimizing the penalty of allocating the resources at lower skill levels: As mentioned earlier, it is assumed that the resources are staff members and each of them has some specialties with different skill levels. Also, the activities require specific skill levels (i.e., novice, standard and senior) in each specialty to be done. Furthermore, the staff member possessing higher skill level in a specialty is certainly capable of working at other lower skill levels in that specialty. Nevertheless, allocating a staff member to an activity requiring a skill level lower than the real skill level of that member causes inefficient utilization in the project regarding the staffs’ paid salaries. In this situation, the working time of the staff member is not utilized economically but the salary is fully paid based on the maximum skill level of the member. Consequently, the third objective is to minimize the penalty of not assigning the staff members at their maximum skill levels.

    Since the proposed MOMSPSP considering three aforementioned conflicting objectives is NP-hard, two meta-heuristic algorithms of PSO and DE are developed. Several numerical examples are also solved using the ε- constraint method as an exact one and the designed PSO and DE algorithms in order to investigate the efficiency of the developed algorithms.

    The remainder of this paper is as follows. Section 2 includes a novel mathematical model for MOMSPSP. DE and PSO algorithms are presented in Section 3. Section 4 deals with the computational results and validation of the proposed solution methods. Finally, Section 5 sums up the concluding remarks and future researches.


    In this section, in order to formulate the mathematical model for MOMSPSP, indices, parameters, and decision variables are firstly presented. Afterward, a mathematical model for the considered MOMSPSP regarding the following assumptions is defined.

    • 1. All required resources for the project are of human force type and are always available.

    • 2. An activity to be done needs various specialties with certain skill levels.

    • 3. The resources (i.e., staff members) possess different skill levels of diverse specialties, thus a staff member can be assigned to the activities based on his/her real skill level or a lower one regarding his/her specialty.

    • 4. The execution time of each activity is considered as an integer number.

    • 5. The staff members assigned to an activity with diverse specialties have to start and finish that activity simultaneously.

    • 6. Each staff member engaged in the project either has definitely a specialty or does not.

    • 7. Each activity has to be executed with no interruption.

    • 8. The time of setting up and allocating resources to the activities is considered negligible.

    • 9. Each staff member can only be assigned at one skill level of a specialty to one activity simultaneously. In other words, a staff member with a specified skill level of a specialty is not allowed to be assigned to more than one activity simultaneously.

    • 10. An activity may require a certain number of staff members in a specified skill level of a specialty to be executed.

    • 11. Each staff member can be assigned to diverse activities in different skill levels of different specialties subject to not being assigned to more than one activity simultaneously.

    • 12. The skill levels of each specialty are classified into three levels; namely, Novice, Standard and Senior.

    • 13. The personnel considered in this study are heterogeneous since they have different specialties with different skill levels as well different costs for idleness and inefficient allocation.

    To explain about personnel skills and specialties, we would say an example. Let us imagine a personal named Alex possessing 2 specialties, i.e., “project scheduling” and “computer programming.” He has skill level “senior” in speciality “project scheduling” and skill level “standard” in speciality “programming.”

    In addition, there is an activity named “project tracking” needs speciality “project scheduling” at skill level “novice.” As mentioned in the assumptions section, it is possible to assign activity “project tracking” to Alex as he has the required speciality “project scheduling.” However, in the case of assigning activity “project tracking” to Alex, a penalty would be incurred since Alex’s skill level “senior” in performing the assigned activity is higher than the required skill level “novice.”

    Indices and Sets

    • n   the index of the last activity

    • Ai   i∈{0,…n} the set of activities of the project

    • A0   the dummy start node of the project

    • An   the dummy end node of the project

    • (Ai, Aj) є E if there exists a precedence relation between Ai and Aj; IEMS-16-307_image1.gif

    • G = (A, E, P)   the precedence graph

    • K   number of specialties

    • L   number of levels per skill (1 is the first level and L is the highest level)

    • M   number of staff members

    • Ok   k ∈{0, …, K} set of specialties

    • O k l    k ∈{0, …, K}, l∈{0, …, L} set of skill levels for specialty k

      : For example, ( O 1 2 ) means capability of a staff member with specialty 1 at skill level 2)

    • Rm   m∈{0, …,M} set of staff members


    • Pi   i∈{0,…, n}Duration of activity Ai

    • t   ∈ {0, 1, … Tmax} Starting time of project activity, T max | = i = 1 n P i

    • Sml = l   m{0, …M}, k{0,…, K} the maximum skill level l at which staff member m Rm can perform specialty Ok

    • r m k l = 1    if staff member Rm can perform specialty Ok at level l; 0, otherwise

      It should be noted that every staff member at each specialty can be assigned to skill levels lower than his/ her maximum level, in other words if S m k = 1 , 1 1 , r m k i = 1 .

    • b i k l    The number of required staff members for performing specialty Ok at skill level l during executing activity Ai.

    • C1m   m∈{0, ⋯, M} idleness cost of staff member Rm per day.

    • C2m   m∈{0,⋯,M} Penalty cost of inefficient assigning staff member Rm at a skill level lower than his/her maximum skill level per day

    Decision Variables

    • x i t    1, if activity i starts at time position t; 0 otherwise

    • x i m t    1, if staff member m starts working on activity i at time position t; 0, otherwise

    • y i m k t    1, if staff member m with specialty k is assigned to activity i at skill level l; 0, otherwise

    The developed MOMSPSP is now formulated as a multi-objective mixed-integer nonlinear programming model:

    M i n Z 0 = C max = M i n t = 0 T max Z n t × t

    M i n Z 1 = M i n m = 1 M C 1 m × [ C max l = 1 L k = 1 K i = 1 N ( y i m k l × P i ) ]

    M i n Z 2 = M i n { i = 1 N p i × [ k = 1 K p i × ( m = 1 M { S m k × l = 1 L y i m k l l = 1 L y i m k l × l } × C 2 m ) ] }

    Subject to:    t = 0 T max ( Z i t × t ) + p i t = 0 T max ( Z j t × t ) i , j E

    t = 0 T max X i m t 1 i , m

    X i m t Z i t i , m , t

    X i m t + 1 Z i t + k = 1 K l = 1 L y i m k l i , m , t

    y i m k l r m k l i , m , k , l

    m = 1 M t = 0 T max X i m t = k = 1 K l = 1 L b i k l i

    m = 1 M y i m k l = b i k l i , k , l

    i = 0 N d = t P i + 1 t X i m d 1 m , t

    t = 0 T max X i m t = k = 1 K l = 1 L y i m k l i , m

    l = 1 L y i m k l 1 i , m , k

    t = 0 T max Z i t = 1 i

    X i m t , y i m k l , Z i t { 0 , 1 } i , m , t , k , l

    The model objective functions which are minimizing the project completion time, staff members idleness cost and inefficient allocation penalty of staff members to activities at skill levels lower than their maximum skill level are formulated as Equations (1) to (3), respectively. Equation (1) calculates the project completion time. As a matter of fact, Equation (2) results in minimizing the idleness cost of staff members. In other words, it minimizes the idleness times of staff members engaged in the project. In fact, this objective function is to maximize the utilization level of staff members.

    It is obvious that objective functions (1) and (2) (i.e., Z0 and Z1) are positively and highly correlated. It can be inferred that by improving objective (2), the other one will be improved. As a result, in most cases the optimal dominant front is unique. In addition, because of uniqueness of optimal solution and positive correlation between objective functions (1) and (2), multi-objective approaches will result in one solution (If they all obtain optimal solution). Thus, there is no necessity to implement Paretobased approach for obtaining optimal Pareto front of objective functions (1) and (2). Hence, objective function (1) can be ignored in the model in favor of objective function (2) to avoid the complexity of handling three objectives.

    Objective function (3) (i.e., Z2) leads to minimizing the penalty of not allocating staff members at their maximum skill levels. In fact, in this objective function, the idleness period of a staff member in the project is not important and the main concern is allocating the staff members based on their maximum skill levels. Consequently, objective functions (2) and (3) (i.e., Z1 and Z2) are optimized in a conflict manner.

    Equation (4) ensures complying precedence relations (activity i is the precedent for activity j). Equation (5) guarantees that staff member m can start performing activity i at only one time position.

    Inequality (6) describes that staff member m can start performing activity i at time position t if activity i is planned to be started at time position t. Similarly, Inequality (7) defines if activity i is planned to be started at time position t and staff member m is assigned to perform that activity with a determined skill level and specialty, then staff member m should start preforming activity i at time position t.

    Inequality (8) enforces that a staff member can be assigned to an activity if he/she has the required specialty at necessary skill level. Equation (9) guarantees that the total number of staff members allocated to one activity has to be equal to the number of staff members with diverse specialties at different skill levels required for executing that activity. Equation (10) ensures that the total number of staff members at specific skill level l of specialty k allocated to activity i must be equal to the number of staff members required at the considered skill level l of specialty k toexecute activity i.

    As mentioned earlier, preemption is not permitted in allocation of staff members at different skill levels of various specialties to the activities of project in the considered MSPSP. In fact, if a staff member is assigned to an activity, he/she is not allowed to get involved in another activity until it is completed. Inequality (11) ensures this rule. Equality (12) describes that if a staff member is assigned to an activity, he/she should start working on that activity at only one time position with a certain specialty at the required skill level. Inequality (13) ensures that a staff member with a certain specialty can be assigned to an activity at only one skill level. Equation (14) certifies that every activity has to be started at a certain time position of planning horizon. Set of Constraints (15) determine that all decision variables used in the model are binary.

    Since the MSPSP is NP-hard from the computational complexity viewpoint, it is impossible to solve the proposed mathematical model in real sizes using exact solution methods in spite of its linear nature. Hence, two efficient meta-heuristic-algorithms, namely DE and PSO are developed to solve such a hard problem.


    3.1.Generating the Initial Solutions

    Since the complexity class of the presented model is NP-hard, it cannot be solved in real sizes despite the fact it is linear. This matter causes DE and PSO meta-heuristic algorithms to be applied to solve the problem. Solution representation method, DE and PSO methodologies, and the proposed structure for these algorithms are discussed as follows.

    In this study, in order to generate the initial solutions to execute the activities sequentially, forward scheduling procedure (Briand, 2009) has been used with some modifications. Solution representation is made up of two parts. The first part is a one-dimensional matrix where the number of columns is equal to the number of activities used to represent the sequence of executing the activities.

    This part is depicted by an example in Table 1 where each element indicates the sequence of executing activities. In Table 1, the project activities are performed in a sequence based on the numbers in the second row. For instance, activity 3 is performed in sequence 2.

    The second part that as a one-dimensional matrix is used to represent the assignment of staff members to the skill levels of specialties required for each activity. This part is also illustrated by an example in Table 2 where the number of columns is equal to the number of skill levels of different specialties required to execute activity i.

    In this example, a project is assumed with 2 skill levels of 2 specialties and 5 staff members. As seen in Table 2, staff members 1 and 5 are assigned to activity i requiring 2 staff members at skill level 1 of specialty 1. Also, staff member 4 is assigned to activity i requiring 1 staff member in speciality 1 at skill level 2. Furthermore, staff member 2 is assigned to activity i requiring 1 staff member in speciality 2 at skill level 1. The skill level 2 of specialty 2 is not needed to perform activity i. In continuance, the propounded solution representation method is used to develop the proposed meta-heuristic algorithms.

    3.2.Proposed DE Algorithm

    DE algorithm is one of the computational methods for optimization of real functions with continuous values using evolutionary strategies presented first by Storn and Price (Price et al., 2006; Storn and Price, 1997). Similar to all evolutionary algorithms, two mutation and crossover operators are also employed in DE algorithm. In fact, the mutation operator generates a trial vector for each solution and crossover operator combines a pair of trial vectors to create new solutions (i.e., offsprings). Since the considered problem is a multi-objective one in discrete space, the proposed DE has to be updated to solve the problem (Qin et al., 2010).

    The general structure of the designed DE algorithm is depicted in Figure 1 and described as follows. First, an initial population with the size of N-pop is generated based on the forward scheduling procedure and then objective values of initial population are calculated. Repository set is considered as a store for non-dominated solutions with a capacity between 10 to (2×N-pop).

    The initial solutions are copied and stored in the Repository set based on non-dominated sorting rule and crowding distance metric defined in fast and elitist nondominated sorting genetic algorithm (NSGA-II) Deb et al. (Deb et al., 2002).

    The other components of the proposed DE including mutation operator, crossover operator, repair procedure, selection of candidate solutions strategy, population updating procedure and Repository set updating procedure are detailed out as follows.

    3.2.1.Mutation Operator

    In the DE algorithm developed in this research, mutation operator is used like the basic DE algorithm (Storn and Price, 1997) with some modification in order to search in the solution space (Qin et al., 2010). The mutation vector operates as follows.

    For both part of each feasible solution Si in the population, three indices p, j and k are generated randomly and their corresponding solutions (Sp, Sj, Sk) are called from the Repository set, unlike the basic DE algorithm where the corresponding solutions are recalled from the parent population. Then, a mutation vector is produced for each solution structure i using Equation (16).

    M i = s p + ( s j s k ) × F    i = 1 , 2 , , N p o p

    In Equation (18), subtraction and addition operators are defined element by element for the solutions. Also, F is a number smaller than 1 to control the population’s evolution rate.

    This vector is illustrated by an example in Figures 2(a)-3. The first parts of initial solutions are brought in Figures 2(a), Figures 2(b) and Figures 2(c). Regarding Equation (16) and F = 0.5, Mi is simply obtained as shown in Figure 3. Also, the second part can be similarly changed by the mutation vector.

    3.2.2.Crossover Operator

    After creating the mutation vector, it is necessary to determine a crossover operator for each part of the solutions. In the proposed algorithm, the crossover operator is applied as Equation (17).

    T i = { M i r i < C R S i o t h e r w i s e

    where CR is a random number in the interval [0, 1] and ri is a matrix whose elements are randomly generated in the in the interval [0, 1] for each Si and Mi. It is noted that Equation (17) is applied on both parts of the solution structure. Figure 4 illustrates the first part of feasible initial solution generated for a typical example. CR = 0.4 and ri in accordance with Figure 5 will result in Ti as seen in Figure 6.

    3.2.3.Repair Procedure

    Regarding previous s teps of DE, currently there is a solution Ti for each Si and its corresponding Mi that is used to convert solution Si to a feasible one by repair procedure. In fact, all solutions are repaired to become new feasible solutions as much as possible. Repairing in the proposed DE is performed as follows:

    3.2.3 .1.Repairing the First Part of Solutions

    In this section, with respect to precedence relations and existing Ti, the first part of new solution structure is generated. First, a new matrix is generated from matrix Ti where the activities are in the form of a permutation. Actually, the content of this new matrix is the position of ascended content of matrix Ti and is named as matrix Ri. Then, the project network is reviewed and activities that can be scheduled are considered as competitors. In this competition, the activity is winner (is scheduled in the current solution) whose corresponding cell number in matrix Ri is smaller. After scheduling the winner activity, there is still a set of competitors among which the winner is selected. This process continues until all project activities become feasible in terms of precedence relations. the Second Part of Solutions

    In this case, since there is no limitation of precedence relations, matrix Ti where the number of columns is equal to total skill levels of all specialties required for activity i is generated. Firstly, |Ti| is calculated and is named as matrix Bi. The elements of this matrix show the priority of skill levels of each speciality for allocation. Each element of Ti whose corresponding Bi is smaller has higher priority for allocation. This allocation continues until the staff members are assigned feasibly to all project activities based on their required specialties and skill levels. After repairing the both parts of a solution to make allocations feasible, the candidate solution is made based on forward scheduling procedure and all project activities are scheduled.

    3.2.4.Selection of Candidate Solutions Strategy

    After repairing the solutions generated by the crossover operator, fitness of solutions has to be evaluated. Since in the multi-objective optimization problems, it is impossible to compare directly the dominance of one solution to another one regarding only one objective function, it is inevitable to employ dominance concept and non-dominated solutions. Comparing the candidate solutions with parent solutions, three cases come up:

    • 1) The candidate solution dominates the parent solution and is substituted with the parent in the population.

    • 2) The parent solution dominates the candidate solution and the candidate solution is dismissed.

    • 3) Both solutions are considered as non-dominated ones in comparison with each other and in this case, the candidate solution is added to population.

    These steps are repeated until N-pop candidates are evaluated. Afterward, the population size will be between N-pop and 2×N-pop.

    3.2.5.Population Updating Procedure

    If the population size gets greater than N-pop, it is reduced to N-pop for the next generation. Truncation operation includes sorting the non-dominated solutions and evaluating the population of a Pareto front based on crowding distance metric. The superior N-pop solutions based on these metrics are transferred to the next generation as parent (Qin et al., 2010).

    3.2.6.Repository Set Updating Procedure

    The population obtained from the population updating procedure is added to the Repository set. Then, sorting the non-dominated solutions leads to 3 cases:

    • 1) The number of obtained non-dominated members (rank1) is less than the minimum capacity of Repository set. In this case, the rest of required members is supplied from rank2 similar to NSGA II. The obtained non-dominated members (rank1, rank2) remain in Repository set and the other dominated ones are dismissed.

    • 2) The number of obtained non-dominated members (rank1) is between the minimum and maximum capacity of Repository set. In this case, the obtained non-dominated members remain in Repository set and the other dominated ones are dismissed.

    • 3) The number of obtained non-dominated members (rank1) is greater than the maximum capacity of Repository set. In this case, 2×N-pop member’s which are superior in terms of crowding distance metric remain in Repository set and the rest of members are dismissed.

    In this stage, Repository set has been updated to use in the next iteration. Then, termination condition is checked. If it is satisfied, the algorithm stops and the obtained nondominated solutions existing in Repository set (rank1) are considered as final solutions.

    3.3.Proposed PSO Algorithm

    Particel swarm optimization is a random populationbased optimization algorithm presented by Eberhart and Kennedy (Eberhart and Kennedy, 1995). It starts by taking a swarm including N particles in a D-dimensional search space seeking a global optimal solution. The vectors Xi = (Xi1, …, XiD), Vi = (Vi1, …, ViD) and Pbesti = (Pbesti1, …, PbeastiD) are position, velocity and the best individual position for each particle, respectively. The position of each particle indicates a potential solution for the optimization problem. The fitness of a solution is evaluated through a pre-defined objective function. The best position visited by particle i is denoted by the best individual position vector Pbesti. Also, gbest = (gbest1, …, gbestD) is a vector for the best global positions visited in each direction by the whole swarm. For each particle, the current position in iteration t and the velocity in iteration t+1 determine the particle position in the next iteration t+1.

    Therefore, V i t and V i t + 1 stand for the velocities of particle i in iterations t and t+1, X i t and X i t + 1 are the positions of particle i in iterations t and t+1, P b e s t i t is the best position of particle i in iteration t and g b e s t i t is the best position among all particles in iteration t. C1 and C2 are the positive constants, r1 and r2 are the random numbers between 0 and 1 and ω is the inertia weight.

    Since MOMSPSP is a multi-objective problem with discrete space, the algorithm for solving the problems in this space has to be updated (dos Santos Coelho, 2009). In the designed PSO algorithm, the swarm size is expressed by N-particle. The initial swarm is generated based on the forward scheduling rule and the the objective functions for initial particles are calculated. Similar to DE algorithm, the Repository set has been taken as a non-dominated solutions storage whose minimum and maximum capacity are 10 and 2×N-particle. Also, the initial swarm solutions are stored in the Repository set based on the non-dominated solutions sorting procedure and crowding distance metric introduced by Deb et al. (2002). In Figure 7, PSO algorithm structure is presented to solve the MOMSPSP.

    In the developed PSO algorithm, updating the particle velocity and position, repairing the moved particles, evaluating the swarm fitness, updating Repository members, and updating the best individual position and the best global position are the main steps described as the following:

    3.3.1.Updating Particle Velocity and Position

    Due to the discrete nature of the designed MOMSPSP, a PSO standard coding strategy might not be used directly. Thus, in order to deal with the discrete structure of scheduling problems, increase algorithm efficiency and avoid being trapped in local optimums, a combination of PSO and DE operators have been employed to update the particles positions according to Equation (18) (Wu et al., 2011).

    i f r < 0.5    V i j t + 1 = W × V i j t + C 1 × r 1 × ( P b e s t i j t x i j t )      + C 2 × r 2 × ( g b e s t i j r x i j t )    x i j t + 1 = x i j t + V i j t + 1 E l s e    x i j t + 1 = p × o p t i t ( r r , j ) + ( 1 p )      × ( o p t i t ( r r , j ) x i j t )    V i j t + 1 = x i j t + 1 x i j t E n d

    where r and p are two random numbers between 0 and 1, and optit(rr) is the Repository set member chosen randomly.

    3.3.2.Repairing Moved Particles

    Considering a particle movement from position xti to xt+1i, for each initially generated solution sti, a solution st+1i is generated by PSO algorithm that has to be repaired in the case of infeasibility. In fact, repairing is done for all vectors of some newly-generated infeasible solutions. Repairing solution in the proposed PSO algorithm is carried out for both representation structures similar to DE algorithm described in subsection 3.2.3.

    3.3.3.Evaluating Swarm Fitness and Updating Repository Set

    In order to update the Repository set, all population members are added to the Repository set. Then, sorting non-dominated solutions is done similar to DE algorithm. In this step, the Repository set has been updated for the subsequent iteration. Now, the termination condition is checked. If it is not satisfied, the best individual position and the best global position of each particle have to be updated for the subsequent iteration.

    3.3.4.Updating the best Individual and Global Positions

    Since the number of the non-dominated solutions in the Pareto front is more than 1, each non-dominated solution can be gbest and transfer its position-related data to the current particles. According to the searching behavior of the particles in multi-objective problems, Pbest of a particle can be its current position. Thus, it is ineffective to use Pbest for guiding the particles towards new solutions. It also can be true for gbest. Therefore, to prevent from guiding the particles ineffectively, DE operators have been used with a slight modification Wu et al., 2011). As Equations (19) and (20) describe, the best individual position and the best global position are updated as it follows:

    i f α < 0.5   P b e s t i j t + 1 = x i j t E l s e   P b e s t i j t + 1 = x t ( r 1 , j ) + r a n d        × ( x t ( r 2 , j ) x t ( r 3 , j ) ) E n d

    i f β < 0.5    g b e s t i j t + 1 = x i j t E l s e    g b e s t i j t + 1 = o p t i t ( r 4 , j ) + r a n d        × ( o p t i t ( r 5 , j ) o p t i t ( r 6 , j ) ) E n d

    where rand, α and β are random numbers between 0 and 1. Also, x t ( r 1 ) , x t ( r 2 )  and  x t ( r 3 ) are the current swarm members selected randomly. Similarly, o p t t ( r 4 ) , o p t t ( r 5 )  and  o p t t ( r 6 ) stand for the Repository set members selected randomly.

    Having updated the best individual position and the best global position, these solutions have to be repaired according to the repairing solutions procedure. After repairing the best individual position and the best global position, the particles velocity and position have to be updated and the explained steps have to be repeated until reaching the termination condition, where the nondominated solutions (rank1) found in Repository set are considered as the final solutions for the developed MOMSPSP.


    4.1.Illustrative Numerical Example

    In order to verify the proposed model mathematically and illustrate the effects of incorporated features on its performance, a numerical example is solved by CPLEX solver executed on an Intel_ Quad_core 2.00 GHz Personal Computer with 4 GB RAM. The example is solved in a computational time equal to 39 minutes and 18 seconds. This example includes 7 activities where A0 is the dummy starting activity of the project and A6 is the dummy ending activity of the project, 4 staff members, 3 specialties at 3 skill levels. Table 3 represents the number of staff members with specialty k required for executing activity i at skill level l. For example, activity 4 needs 2 staff members with speciality 2 at skill level 1.

    Table 4 illustrates the penalty cost for executing an activity by a staff member at a skill level lower than his/her maximum skill level and idleness cost of each staff member per day.

    In Table 5, the capability of each staff member to do a speciality at a skill level is given. In Table 6, the maximum skill level of each staff member to do a speciality is presented.

    In Figure 8, precedence network of the project has been shown. In Figure 9, the obtained scheduling solution including the assignment of each staff member to each activity at a skill level and starting and ending times of each activity are represented.

    As seen in the Gantt chart depicted in Figure 9, to execute activity A1 started at time position 1 and finished at time position 10, two staff members including R2 with speciality 1 (O1) at skill level 2 ( O 1 2 ) and R4 with speciality 2 (O2) at skill level 1 ( O 2 1 ) have been assigned. Similarly, the assignment of the other staff members can be interpreted.

    Regarding the staff members assignment due to executing the project activities as depicted in Figure 9, it could be concluded that all model constraints including Equations (4)-(14) are observed. Also, given the optimal scheduling solution obtained by CPLEX solver and its display in Gantt chart, it is perceived that the assignment of staff members to project activities exactly follows the problem assumptions with no violation of model constraints.

    4.2.Parameters Setting for the Proposed Algorithms

    In order to set the parameters existing in the proposed algorithms, Taguchi design of experiments method has been used. In this step, the aim is to figure out the best PSO and DE algorithm parameters as the input variables to obtain the optimal solutions as the input variables. It is worth mentioning that in PSO and DE parameters setting, Mean Ideal Distance (MID) metric has been applied to determine the most suitable experiment. Taguchi method has been used to set PSO algorithm parameters including population size (N-pop), iteration number (N-iteration), individual learning coefficients (C1) and swarm learning coefficients (C2) and DE algorithm parameters including population size (N-pop), iteration number (N-iteration), combination probability (Cr) and control factor (F). Here, Taguchi method has been used for 4 factors at 3 levels for both algorithms. Table 7 presents the factors values at every level for both algorithms so that numbers 1, 2 and 3 stand for the levels of each factor.

    Table 8 and Table 9 show the orthogonal array L9 of 4 factors at 3 levels for the developed PSO and DE algo- rithms, respectively.

    A problem has been produced with random data and in order to achieve more reliable results, it has been tested 10 times. Thus, 10 results have been produced for each Taguchi test. As mentioned earlier, MID metric is used to determine the most appropriate test in the algorithms parameters setting.

    The results from Taguchi method is transformed into S/N rate. Figure 10 and Figure 11 depict average S/N ratio and average response ratio for each factor in PSO algorithm, respectively. Also, Table 10 shows the best level for these factors. Delta value in Table 11 indicates the average range of S/N ratio variations for each factor in PSO algorithm. A factor with a higher delta value achieves a higher rank. As seen in Figure 10 and Figure 11 and Delta values in Table 11, N-pop exerts the highest effect on PSO algorithm performance. C2, N-iteration and C1 are the other factors influencing PSO algorithm, orderly.

    Similarly, in DE parameters setti ng, the results obtained from Taguchi method are represented by Figure 12 and Figure 13 and Table 12 and Table 13. As it could be observed in Figure 12 and Figure 13, and also through Delta value in Table 13, N-pop has the maximum effect on DE algorithm, N-iteration, F and Cr are the other factors influencing DE algorithm, orderly.

    4.3.Computational Efficiency Analysis

    Since there is no test problems available in the Project Scheduling Problem Library (PSPLIB) for the presented MOMSPSP, ten numerical examples with different sizes have been produced randomly based on Table 14 to analyze the mathematical model performance and evaluate the developed algorithms efficiency.

    In this section, the efficiency of PSO and DE algorithms in solving the developed MOMSPSP is compared to CPLEX solver. PSO and DE algorithms have been coded using R2014a MATLAB software and the test problems have been run on an Intel_ Quad_core 2.00 GHz Personal Computer with 4 GB RAM. Comparing the solutions obtained for small-sized test problems by ε- constraint method and the developed algorithms is done to verify the proposed algorithms performance in reaching optimal solutions.

    4.3.1.Comparing the performance of ε-constraint method and the developed PSO and DE algorithms

    In this section, the presented model for the MOMSPSP is solved by ε-constraint method. The Pareto fronts obtained for problem No.1 solved by ε-constraint method are presented in Table 15 and depicted in Figure 14.

    Furthermore, in order to illustrate the consistency of optimal Pareto front obtained by ε-constraint method with Pareto fronts obtained by the developed algorithms, problem No.1 is solved by the proposed algorithms and the obtained Pareto fronts are represented in Table 16 and Table 17 and depicted in Figure 15 and Figure 16.

    Analyzing the Pareto fronts obtained by CPLEX solver and two meta-heuristic algorithms as depicted in Figures 14-16, it is concluded that all fronts have common or close solutions and are also consistent. That confirms the validity of the meta-heuristic algorithms developed for the proposed mathematical model in reaching optimal solutions for small-sized problems.

    4.3.2.Comparing the Performance of the Developed PSO and DE Algorithms

    Due to the stochastic nature of meta-heuristic algorithms, the produced test problems are solved 10 times by PSO and DE algorithms to obtain more reliable results. The computational times for solving small and large-sized problems by PSO and DE algorithms and CPLEX solver are given in Figure 17. Since the problem is NP-hard, it is impossible to solve the large-sized problems 4-10 with a reasonable time. Regarding Figure 17, it is revealed that in terms of computational times, DE algorithm has a slightly better performance than PSO algorithm. Also, the problem is strongly NP-hard and obtaining exact solutions is time-consuming and overwhelming as the problem size increase.

    To compare the results of PSO and DE algorithms on the test problems, SNS, MID, RAS and SM metrics (Nader et al., 2011) have been employed, which are briefly explained as follows.

    •Spread of Non-dominance Solution (SNS) Metric

    This metric calculates the span of Pareto front obtained by each algorithm. Obviously, the larger SNS is preferred. It is calculated by Equation (21) (Jolai et al., 2013):

    S N S = i = 1 n ( M I D C i ) 2 n 1

    •Mean Ideal Distance (MID) Metric

    Using this metric, the closeness distance between the obtained non-dominated solutions and the ideal point (0 , 0) are gained. This metric is estimated through Equation (22), where n stands for the number of obtained nondominated solutions and f1i and f2i are the 1st and 2nd objective function values for non-dominated solution i. The smaller this metric, the higher priority the algorithm has due to producing solutions with less distance from the ideal point (Behnamian et al., 2009).

    M I D = i = 1 n C i n ,     C i ( f 1 i ) 2 + ( f    i ) 2

    •The Rate of Achievement to Objectives Simultaneously (RAS) Metric

    This metric is seeking a state of equilibrium in optimizing conflicting objectives simultaneously and is calculated based on Equation (23) (Karimi et al., 2010) where f1i and f2i are the 1st and 2nd objective functions values for non-dominated solution i and n is the number of the obtained non-dominated solutions. The smaller this metric, the higher priority the algorithm has.

    R A S = i = 1 n ( f 1 i F i F i ) + ( f i F i F i ) n ,    F i = min ( f 1 i , f i )

    •Spacing Metric (SM)

    Employing this metric, the diversification uniformity of the non-dominated solutions is calculated according to Equation (24), where n represents the number of obtained non-dominated solutions and di is the minimum Euclidean distance between solution i and the solutions in the non-dominated solutions set calculated according to Equation (25). Also, d is the average of these distances. The smaller SM, the higher priority the algorithm has (Coello et al., 2007).

    S M = 1 ( n 1 ) i = 1 n ( d ¯ d i ) 2

    d i = min j ( | f 1 i f    j | + | f    i f    j | )    i , j = 1 , 2 , , n

    In the following, the results obtained for the performance metrics from each algorithm have been listed with Table 18 and Figure 18.

    Statistical analysis of the performance metrics obtained for PSO and DE algorithms is conducted in the following. First, the fit test (Stein, 1981) is done in Minitab 16 software to achieve the data distribution standardization. Since the metrics values resulted from solving ten test problems by two algorithms PSO and DE have no normal statistical distribution and the size of samples was less than 30, the Non-parametric Mann-Whitney test

    Siegel (1956) has been used to compare these metrics mean (Siegel, 1956). Then, four two-way hypothesis tests have been considered for examining the equality of mean values of MID, RAS, SM and SNS metrics derived by PSO and DE algorithms. In this study, the confidence level was taken as α = 0.05. Equations (26)-(29) indicate the two-way hypothesis tests for SNS, MID, RAS, and SM values, respectively. Mann-Whitney test results for 10 examples where each one has been run for 10 times are summarized as the following:(27)(28)

    H 0 : μ S N S ( D E ) = μ S N S ( P S O ) , H 1 : μ S N S ( D E ) μ S N S ( P S O )

    H 0 : μ M I D ( D E ) = μ M I D ( P S O ) , H 1 : μ M I D ( D E ) μ M I D ( P S O )

    H 0 : μ R A S ( D E ) = μ R A S ( P S O ) , H 1 : μ R A S ( D E ) μ R A S ( P S O )

    H 0 : μ S M ( D E ) = μ S M ( P S O ) , H 1 : μ S M ( D E ) μ S M ( P S O )

    Analyzing the P-value results achieved by Mann- Whitney test as shown in Table 19, it can be concluded that due to P-value being achieved more than α = 0.05 in these four tests, there is no reason to reject hypothesis H0 and there is no meaningful difference between the mean of four metrics obtained by PSO and DE algorithms. However, only regarding Figure 18 in a non-statistical analysis implies that totally DE algorithm has had slightly better performance according to SNS, MID and SM metrics and PSO algorithm is better based on RAS metric.


    In this study, a new multi-objective mathematical model in the MOMSPSP field has been presented to incorporate three objectives including: 1) minimizing the project completion time, 2) minimizing the resources idle time and 3) minimizing the penalty of inefficient allocating the resources. The proposed model has been solved optimally using CPLEX solver for few test examples in small sizes. Concerning high computational time to solve this problem optimally as a strongly NPhard one in real sizes, two efficient meta-heuristic algorithms PSO and DE have been developed. Analyzing the results achieved by the proposed meta-heuristic algorithms in solving the problem in different sizes draw the following conclusions:

    • Solving the developed mathematical model has created valid and feasible solutions for small-sized problems.

    • The Pareto fronts obtained by the proposed algorithms were conforming to the optimal Pareto fronts resulted from solving the model by ε – constraint method in small-sized problems.

    • With respect to the results of statistical tests, the proposed DE and PSO algorithms have almost equal results in terms of performance metrics.

    • Considering the computational times, DE algorithm has higher performance than PSO algorithm.

    Incorporating other features will be left to future research, such as:

    • Creating efficient lower bounds to be employed by exact methods for solving the developed MOMSPSP.

    • Considering the model parameters and variables as fuzzy numbers.

    • Combining MOMSPSP with the other existing problems in the project scheduling field such as MMRCPSP with crash times for activities.

    • Implementing the model for a case study with real data using the developed meta-heuristics

    This research did not receive any specific grant from funding agencies in the public, commercial, or not-forprofit sectors



    Structure of the proposed DE for MOMSPSP.


    Feasible solution.


    Mutation vector (M1) for solutions (S1, S2, S3).


    Feasible solution Si (parent).


    Matrix ri for the typical example.


    Matrix Ti for the typical example.


    Structure of the proposed PSO algorithm for MOMSPSP.


    Precedence network of the illustrative example.


    Gantt chart of staff member assignment and scheduling activities for the illustrative example.


    Average S/N ratio for each factor in PSO algorithm.


    Average response ratio for each factor in PSO algorithm.


    Average S/N ratio for each factor in DE algorithm.


    Average response ratio for each factor in DE algorithm.


    Two Pareto fronts obtained by ε-constraint method for problem No.1.


    Pareto front obtained by PSO algorithm for problem No.1.


    Pareto fronts obtained by DE algorithm for problem No.1.


    Comparing the performance of PSO, DE and CPLEX solver based on computational times.


    Comparing the SNS, MID, RAS and SM metrics of PSO and DE.


    An example for the sequence of executing activities in the solution representation

    An example for the assignment of staff members to activity i in the solution representation

    The number of staff members with different specialties required for executing each activity at different skill levels for the illustrative example

    Penalty cost of assigning a staff member to an activity at a skill level lower than his/her real skill level and idleness cost of each staff member in the illustrative example

    Capability of each staff member to do a speciality at a skill level in the illustrative example

    Maximum skill level of each staff member to do a speciality in the illustrative example

    Different levels for the parameters of PSO and DE algorithms

    Orthogonal array L9 for the parameters of PSO algorithm

    Orthogonal array L9 for the parameters of DE algorithm

    Best level for each factor in PSO algorithm

    Ranking the effects of each factor in PSO algorithm

    Best levels for each factor in DE algorithm

    Ranking the effects of each factor in DE algorithm

    Data for ten numerical examples

    The values obtained for two objective functions by ε -constraint method for problem No.1

    The values obtained for two objective functions by PSO algorithm for problem No.1

    The values obtained for two objective functions by DE algorithm for problem No.1

    Comparing the performance metrics of PSO and DE

    The Mann-Whitney test results achieved by PSO and DE algorithms for MID, RAS, SM and SNS values


    1. Al-Anzi F.S. , Al-Zame K. , Allahverdi A. (2010) Weighted multi-skill resources project scheduling , Journal of Software Engineering and Applications, Vol.3 (12) ; pp.1125
    2. Behnamian J. , Ghomi S.F. , Zandieh M. (2009) A multi-phase covering Pareto-optimal front method to multi-objective scheduling in a realistic hybrid flowshop using a hybrid metaheuristic , Expert Syst. Appl, Vol.36 (8) ; pp.11057-11069
    3. Bellenguez O. , Néron E. (2004) Lower bounds for the multi-skill project scheduling problem with hierarchical levels of skills , International Conference on the Practice and Theory of Automated Timetabling, ; pp.229-243
    4. Briand C. (2009) A new any-order schedule generation scheme for resource-constrained project scheduling , Oper. Res, Vol.43 (3) ; pp.297-308
    5. Coello C.A. , Lamont G.B. , Van Veldhuizen D.A. (2007) Evolutionary Algorithms For Solving Multi-Objective Problems 5, Springer,
    6. Corominas A. , Pastor R. , Rodríguez E. (2006) Rotational allocation of tasks to multifunctional workers in a service industry , Int. J. Prod. Econ, Vol.103 (1) ; pp.3-9
    7. Deb K. , Pratap A. , Agarwal S. , Meyarivan T.A. (2002) A fast and elitist multiobjective genetic algorithm: NSGA-II , IEEE Trans. Evol. Comput, Vol.6 (2) ; pp.182-197
    8. Dhib C. , Kooli A. , Soukhal A. , Néron E. (2012) Lower bounds for a multi-skill project scheduling problem , Operations Research Proceedings, Springer,
    9. dos Santos Coelho L. (2009) Multi-objective swarm intelligent systems: Theory & experiences, 261., Springer,
    10. Eberhart R. , Kennedy J. (1995) A new optimizer using particle swarm theory, Micro Machine and Human Science , Proceedings of the Sixth International Symposium on IEEE, ; pp.39-43
    11. Fırat M. , Hurkens C.A. (2012) An improved mipbased approach for a multi-skill workforce scheduling problem , J. Sched, Vol.15 (3) ; pp.363-380
    12. Hartmann S. , Briskorn D. (2010) A survey of variants and extensions of the resource-constrained project scheduling problem , Eur. J. Oper. Res, Vol.207 (1) ; pp.1-14
    13. Heimerl C. , Kolisch R. (2010) Scheduling and staffing multiple projects with a multi-skilled workforce , OR-Spectrum, Vol.32 (2) ; pp.343-368
    14. Jia Q. , Seo Y. (2013) An improved particle swarm optimization for the resource-constrained project scheduling problem , Int. J. Adv. Manuf. Technol, Vol.67 (9-12) ; pp.2627-2638
    15. Jolai F. , Asefi H. , Rabiee M. , Ramezani P. (2013) Bi-objective simulated annealing approaches for no-wait two-stage flexible flow shop scheduling problem , Sci. Iran, Vol.20 (3) ; pp.861-872
    16. Kadrou Y. , Najid N.M. (2006) A new heuristic to solve RCPSP with multiple execution modes and multi-skilled labor , Computational Engineering in Systems Applications, IMACS Multiconference on, IEEE,
    17. Karimi N. , Zandieh M. , Karamooz H.R. (2010) Bi-objective group scheduling in hybrid flexible flowshop: A multi-phase approach , Expert Syst. Appl, Vol.37 (6) ; pp.4024-4032
    18. Kazemipoor H. , Tavakkoli-Moghaddam R. (2013) Solving a novel multi- skilled project scheduling model by scatter search , S. Afr. J. Ind. Eng, Vol.24 (1) ; pp.121-131
    19. Kazemipoor H. , Tavakkoli-Moghaddam R. (2012) Solving a mixed-integer linear programming model for a multi-skilled project scheduling problem by simulated annealing , Management Science Letters, Vol.2 (2) ; pp.681-688
    20. Kazemipoor H. , Tavakkoli-Moghaddam R. (2013) A differential evolution algorithm to solve multi-skilled project portfolio scheduling problems , The International Journal of Advanced Manufacturing Technology, Vol.64 (5-8) ; pp.1099-1111
    21. Khalili S. , Najafi A.A. , Niaki S.T. (2013) Biobjective resource constrained project scheduling problem with makespan and net present value criteria: Two meta-heuristic algorithms , Int. J. Adv. Manuf. Technol, Vol.69 (1-4) ; pp.617-626
    22. Li H. , Womer K. (2009) Scheduling projects with multi-skilled personnel by a hybrid MILP/CP benders decomposition algorithm , J. Sched, Vol.12 (3) ; pp.281-298
    23. Liu S.S. , Wang C.J. (2012) Optimizing linear project scheduling with multi-skilled crews , Autom. Construct, Vol.24 (2) ; pp.16-23
    24. Montoya C. , Bellenguez-Morineau O. , Pinson E. , Rivreau D. (2014) Branch-and-price approach for the multi-skill project scheduling problem , Optim. Lett, Vol.8 (5) ; pp.1721-1734
    25. Myszkowski P.B. , Skowroński M.E. , Olech Ł.P. , Oślizło K. (2015) Hybrid ant colony optimization in solving multi-skill resource-constrained project scheduling problem , Soft Comput, Vol.19 (12) ; pp.3599-3619
    26. Naderi B. , Ghomi S.F. , Aminnayeri M. , Zandieh M. (2011) Scheduling open shops with parallel machines to minimize total completion time , J. Comput. Appl. Math, Vol.235 (5) ; pp.1275-1287
    27. Neron E. (2002) Lower bounds for the multi-skill project scheduling problem , Proceeding of the Eighth International Workshop on Project Management and Scheduling, ; pp.274-277
    28. Neron E. , Boptista D. (2002) Heuristics for the multi-skill project scheduling problem , International Symposium On Combinatorial Optimization,
    29. Pessan C. , Bellenguez-Morineau O. , Néron E. (2007) Multi-skill project scheduling problem and total productive maintenance , Proceedings of MISTA, ; pp.608-610
    30. Price K. , Storn R.M. , Lampinen J.A. (2006) Differential evolution: A practical approach to global optimization, Springer Science and Business Media,
    31. Qin H. , Zhou J. , Lu Y. , Wang Y. , Zhang Y. (2010) Multi-objective differential evolution with adaptive Cauchy mutation for short-term multi-objective optimal hydro-thermal scheduling , Energy Convers. Manage, Vol.51 (4) ; pp.788-794
    32. Santos M. A. , Tereso A. P. (2011) On the multi-mode, multi-skill resource constrained Project Scheduling Problem: A Software Application , Soft Computing in Industrial Applications, Springer,
    33. Siegel S. (1956) Nonparametric statistics for the behavioral sciences,
    34. Stein C.M. (1981) Estimation of the mean of a multivariate normal distribution , Ann. Stat, ; pp.1135-1151
    35. Storn R. , Price K. (1997) Differential evolution–a simple and efficient heuristic for global optimization over continuous spaces , J. Glob. Optim, Vol.11 (4) ; pp.341-359
    36. Wu D. , Ogawa M. , Suzuki Y. , Ogai H. , Kusaka J. (2011) Modified multi-objective particle swarm optimization: Application to optimization of diesel engine control parameter , SICE Journal of Control, Measurement, and System Integration, Vol.3 ; pp.315-323
    37. Yoosefzadeh H.R. , Tareghian H.R. (2013) Hybrid solution method for resource-constrained project scheduling problem using a new schedule generator , Int. J. Adv. Manuf. Technol, Vol.66 (5-8) ; pp.1171-1180