1.INTRODUCTION
Multiskilled Project Scheduling Problem (MSPSP) is a combination of the classic ResourceConstrained Project Scheduling Problem (RCPSP) and MultiPurpose 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 MultiMode 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 multilevel human resources and used lower bound method for minimizing the total project time. Corominas et al. (2006) studied the multitask 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 mixedinteger model for MSPSP using simulated annealing algorithm. Also, Kazemipoor et al. (2013) suggested a novel goal programming model for the multiobjective multiskill 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 multimode multilevel multiskill 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 singleskill and multiskill work force in order to solve the Project Scheduling Problem (PSP) and applied a combination of constrained programming (CP) approach and scheduling rulesbased 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 multiskill resources to minimize the human forcerelated costs. Dhib et al. (2012) presented new lower bounds for multiskill project scheduling solution. Myszkowski et al. (2015) developed an ant colony optimization for a multiskill RCPSP. Jia and Seo (2013) proposed an improved Particle Swarm Optimization (PSO) algorithm for a RCPSP. Khalili et al. (2013) proposed two metaheuristic algorithms of multipopulation and twophase subpopulation GAs to simultaneously minimize the project makespan and maximize the project net present value (NPV) of a multiobjective 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 nonallocated resources and inefficient utilization costs of resources allocated to activities at lower skill levels is the main motivation of this paper to present a multiobjective 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 fulltime 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 NPhard, two metaheuristic 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.
2.THE PROPOSED MATHEMATICAL MODEL FOR MOMSPSP
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

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

A_{0} the dummy start node of the project

A_{n} the dummy end node of the project

(A_{i}, A_{j}) є E if there exists a precedence relation between A_{i} and A_{j};

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

O_{k} k ∈{0, …, K} set of specialties

${O}_{k}^{l}$ k ∈{0, …, K}, l∈{0, …, L} set of skill levels for specialty k
: For example, (${\text{O}}_{\text{1}}^{\text{2}}$) means capability of a staff member with specialty 1 at skill level 2)

R_{m} m∈{0, …,M} set of staff members
Parameters

P_{i} i∈{0,…, n}Duration of activity A_{i}

t ∈ {0, 1, … T_{max}} Starting time of project activity, ${T}_{\mathrm{max}}={\displaystyle {\sum}_{i=1}^{n}{P}_{i}}$

S_{ml} = l m{0, …M}, k{0,…, K} the maximum skill level l at which staff member m R_{m} can perform specialty O_{k}

${r}_{mk}^{l}=1$ if staff member R_{m} can perform specialty O_{k} 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}_{mk}=1,\text{\hspace{0.17em}}\forall \text{\hspace{0.17em}}{1}^{\prime}\le 1,\text{\hspace{0.17em}}{r}_{mk}^{i}=1$.

${b}_{ik}^{l}$ The number of required staff members for performing specialty O_{k} at skill level l during executing activity A_{i}.

C_{1m} m∈{0, ⋯, M} idleness cost of staff member R_{m} per day.

C_{2m} m∈{0,⋯,M} Penalty cost of inefficient assigning staff member R_{m} 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}_{im}^{t}$ 1, if staff member m starts working on activity i at time position t; 0, otherwise

${y}_{imk}^{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 multiobjective mixedinteger nonlinear programming model:
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., Z_{0} and Z_{1}) 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), multiobjective 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., Z_{2}) 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., Z_{1} and Z_{2}) 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 NPhard 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 metaheuristicalgorithms, namely DE and PSO are developed to solve such a hard problem.
3.METAHEURISTIC SOLUTION ALGORITHMS
3.1.Generating the Initial Solutions
Since the complexity class of the presented model is NPhard, it cannot be solved in real sizes despite the fact it is linear. This matter causes DE and PSO metaheuristic 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 onedimensional 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 onedimensional 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 metaheuristic 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 multiobjective 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 Npop 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 nondominated solutions with a capacity between 10 to (2×Npop).
The initial solutions are copied and stored in the Repository set based on nondominated sorting rule and crowding distance metric defined in fast and elitist nondominated sorting genetic algorithm (NSGAII) 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 (S_{p}, S_{j}, S_{k}) 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).
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).
where C_{R} is a random number in the interval [0, 1] and r_{i} is a matrix whose elements are randomly generated in the in the interval [0, 1] for each S_{i} and M_{i}. 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. C_{R} = 0.4 and r_{i} in accordance with Figure 5 will result in T_{i} as seen in Figure 6.
3.2.3.Repair Procedure
Regarding previous s teps of DE, currently there is a solution T_{i} for each S_{i} and its corresponding M_{i} that is used to convert solution S_{i} 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 T_{i}, the first part of new solution structure is generated. First, a new matrix is generated from matrix T_{i} 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 T_{i} and is named as matrix R_{i}. 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 R_{i} 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.
3.2.3.2.Repairing the Second Part of Solutions
In this case, since there is no limitation of precedence relations, matrix T_{i} where the number of columns is equal to total skill levels of all specialties required for activity i is generated. Firstly, T_{i} is calculated and is named as matrix B_{i}. The elements of this matrix show the priority of skill levels of each speciality for allocation. Each element of T_{i} whose corresponding B_{i} 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 multiobjective 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 nondominated 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 nondominated ones in comparison with each other and in this case, the candidate solution is added to population.
These steps are repeated until Npop candidates are evaluated. Afterward, the population size will be between Npop and 2×Npop.
3.2.5.Population Updating Procedure
If the population size gets greater than Npop, it is reduced to Npop for the next generation. Truncation operation includes sorting the nondominated solutions and evaluating the population of a Pareto front based on crowding distance metric. The superior Npop 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 nondominated solutions leads to 3 cases:

1) The number of obtained nondominated members (rank_{1}) is less than the minimum capacity of Repository set. In this case, the rest of required members is supplied from rank_{2} similar to NSGA II. The obtained nondominated members (rank_{1}, rank_{2}) remain in Repository set and the other dominated ones are dismissed.

2) The number of obtained nondominated members (rank_{1}) is between the minimum and maximum capacity of Repository set. In this case, the obtained nondominated members remain in Repository set and the other dominated ones are dismissed.

3) The number of obtained nondominated members (rank_{1}) is greater than the maximum capacity of Repository set. In this case, 2×Npop 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 (rank_{1}) 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 Ddimensional search space seeking a global optimal solution. The vectors X_{i} = (X_{i}_{1}, …, X_{iD}), V_{i} = (V_{i}_{1}, …, V_{iD}) and Pbest_{i} = (Pbest_{i}_{1}, …, Pbeast_{iD}) 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 predefined objective function. The best position visited by particle i is denoted by the best individual position vector Pbest_{i}. Also, gbest = (gbest_{1}, …, gbest_{D}) 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, $Pbes{t}_{i}^{t}$ is the best position of particle i in iteration t and $gbes{t}_{i}^{t}$ is the best position among all particles in iteration t. C_{1} and C_{2} are the positive constants, r_{1} and r_{2} are the random numbers between 0 and 1 and ω is the inertia weight.
Since MOMSPSP is a multiobjective 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 Nparticle. 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 nondominated solutions storage whose minimum and maximum capacity are 10 and 2×Nparticle. Also, the initial swarm solutions are stored in the Repository set based on the nondominated 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).
where r and p are two random numbers between 0 and 1, and opti^{t}(rr) is the Repository set member chosen randomly.
3.3.2.Repairing Moved Particles
Considering a particle movement from position x^{t}_{i} to x^{t+1}_{i}, for each initially generated solution s^{t}_{i}, a solution s^{t+1}_{i} 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 newlygenerated 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 nondominated 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 nondominated solutions in the Pareto front is more than 1, each nondominated solution can be g_{best} and transfer its positionrelated data to the current particles. According to the searching behavior of the particles in multiobjective problems, P_{best} of a particle can be its current position. Thus, it is ineffective to use P_{best} for guiding the particles towards new solutions. It also can be true for g_{best}. 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:
where rand, α and β are random numbers between 0 and 1. Also, ${x}^{t}\left({r}_{1}\right),\text{\hspace{0.17em}}{x}^{t}\left({r}_{2}\right)\text{\hspace{0.17em}and\hspace{0.17em}}{x}^{t}\left({r}_{3}\right)$ are the current swarm members selected randomly. Similarly, $op{t}^{t}\left({r}_{4}\right),\text{\hspace{0.17em}}op{t}^{t}\left({r}_{5}\right)\text{\hspace{0.17em}and\hspace{0.17em}}op{t}^{t}\left({r}_{6}\right)$ 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 (rank_{1}) found in Repository set are considered as the final solutions for the developed MOMSPSP.
4.COMPUTATIONAL RESULTS
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 A_{1} started at time position 1 and finished at time position 10, two staff members including R_{2} with speciality 1 (O_{1}) at skill level 2 (${O}_{1}^{2}$) and R_{4} with speciality 2 (O_{2}) 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 (Npop), iteration number (Niteration), individual learning coefficients (C_{1}) and swarm learning coefficients (C_{2}) and DE algorithm parameters including population size (Npop), iteration number (Niteration), combination probability (C_{r}) 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, Npop exerts the highest effect on PSO algorithm performance. C_{2}, Niteration and C_{1} 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, Npop has the maximum effect on DE algorithm, Niteration, F and C_{r} 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 smallsized 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 metaheuristic algorithms as depicted in Figures 1416, it is concluded that all fronts have common or close solutions and are also consistent. That confirms the validity of the metaheuristic algorithms developed for the proposed mathematical model in reaching optimal solutions for smallsized problems.
4.3.2.Comparing the Performance of the Developed PSO and DE Algorithms
Due to the stochastic nature of metaheuristic 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 largesized problems by PSO and DE algorithms and CPLEX solver are given in Figure 17. Since the problem is NPhard, it is impossible to solve the largesized problems 410 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 NPhard and obtaining exact solutions is timeconsuming 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 Nondominance 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):
•Mean Ideal Distance (MID) Metric
Using this metric, the closeness distance between the obtained nondominated 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 f_{1i} and f_{2i} are the 1st and 2nd objective function values for nondominated 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).
•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 f_{1i} and f_{2i} are the 1st and 2nd objective functions values for nondominated solution i and n is the number of the obtained nondominated solutions. The smaller this metric, the higher priority the algorithm has.
•Spacing Metric (SM)
Employing this metric, the diversification uniformity of the nondominated solutions is calculated according to Equation (24), where n represents the number of obtained nondominated solutions and d_{i} is the minimum Euclidean distance between solution i and the solutions in the nondominated 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).
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 Nonparametric MannWhitney test
Siegel (1956) has been used to compare these metrics mean (Siegel, 1956). Then, four twoway 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 twoway hypothesis tests for SNS, MID, RAS, and SM values, respectively. MannWhitney test results for 10 examples where each one has been run for 10 times are summarized as the following:(27)(28)
Analyzing the Pvalue results achieved by Mann Whitney test as shown in Table 19, it can be concluded that due to Pvalue 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 nonstatistical 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.
5.CONCLUSIONS
In this study, a new multiobjective 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 metaheuristic algorithms PSO and DE have been developed. Analyzing the results achieved by the proposed metaheuristic algorithms in solving the problem in different sizes draw the following conclusions:

Solving the developed mathematical model has created valid and feasible solutions for smallsized 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 smallsized 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 metaheuristics
This research did not receive any specific grant from funding agencies in the public, commercial, or notforprofit sectors