A Modified PSO Algorithm for Minimizing the Total Costs of Resources in MRCPSP

We introduce a multimode resource-constrained project scheduling problem with finish-to-start precedence relations among project activities, considering renewable and nonrenewable resource costs. We assume that renewable resources are rented and are not available in all periods of time of the project. In other words, there is a mandated ready date as well as a due date for each renewable resource type so that no resource is used before its ready date. However, the resources are permitted to be used after their due dates by paying penalty costs. The objective is to minimize the total costs of both renewable and nonrenewable resource usage. This problem is called multimode resource-constrained project scheduling problem with minimization of total weighted resource tardiness penalty cost MRCPSP-TWRTPC , where, for each activity, both renewable and nonrenewable resource requirements depend on activity mode. For this problem, we present a metaheuristic algorithm based on a modified Particle Swarm Optimization PSO approach introduced by Tchomté and Gourgand which uses a modified rule for the displacement of particles. We present a prioritization rule for activities and several improvement and local search methods. Experimental results reveal the effectiveness and efficiency of the proposed algorithm for the problem in question.


Introduction
The resource-constrained project scheduling problem RCPSP is the scheduling of project activities subject to precedence relations as well as renewable resource constraints with the objective of minimizing the makespan of the project.Each nonpreemptive activity in RCPSP can be done in a single mode.For more information on RCPSP and solution methods, we refer to Demeulemeester and Herroelen 1 .In the multimode RCPSP MRCPSP , a set of problem, the objective function is minimization of total costs of renewable and nonrenewable resources.
MRCPSP-TWRTPC is a generalization of the RCPSP problem, and, considering the NP-hardness of RCPSP 12, 13 , the MRCPSP-TWRTPC problem is NP-hard as well, and hence, metaheuristic method is the practical approach.In the remainder of this paper, we introduce a metaheuristic-based PSO algorithm for solving this problem.PSO, in its present form, dates back to 1990s; however, in this short period, PSO has shown good performance in a variety of application domains, particularly in the constrained optimization problems.Many researchers studied PSO widely and proposed several modifications.In this paper, we use a modified PSO algorithm developed and used by Tchomté and Gourgand for solving RCPSP efficiently 14 .
The rest of this paper is organized as follows.In the next section, MRCPSP-TWRTPC is described in detail and is formulated in a mathematical model.In Section 3, a description of the PSO algorithm and its modifications is presented, and in Section 4 an algorithm based on modified PSO introduced by Tchomté and Gourgand 14 is explained.Section 5 is the experimental analysis.Finally, Section 6 concludes the work.

Problem Description
In MRCPSP-TWRTPC, a project is to be scheduled in order to minimize its total costs.Resources available for completing project activities can be classified as either renewable or nonrenewable.Activity j may have a number of execution modes M j .Each activity mode specifies the activity duration and the activity requirements for the certain amount of renewable and nonrenewable resources.Each type of limited renewable resource is rented for a fixed time interval, starting from its ready time and ending with its due date, and is not available before its ready time but can be used after its due date with tardiness penalty cost.Nonrenewable resources are not limited.All activities are ready at the beginning of the project, and no preemption is permitted.If an activity is started under a specific mode, the activity mode cannot be changed.Activity j executed in mode m has duration d jm and requires r jmk units of renewable resource k and n jmk units of nonrenewable resource k.The project network is depicted by an activity on node AON representation with finish-to-start precedence relations and zero time lag.Dummy activities 1 and n correspond to start and completion of the project.The list of activities is topologically numbered; that is, each predecessor of every activity has a smaller number than the number of activity itself.Also, we define the earliest and latest start time of activity j by EST j and LST j , respectively.EST j s and LST j s are computed by CPM forward and backward passes using the mode with shortest duration for each activity and assigning LST n LFT n T , where T is an upper bound for project makespan determined by any valid method, such as the sum of the longest duration of entire project activities plus the ready times of renewable resources.Consequently, each activity j can only be performed in time period EST j , LST j .
We define problem parameters as follows: We also define the decision variables as follows: x jmτ ⎧ ⎨ ⎩ 1, if activity j is started under mode m in period τ, 0, otherwise, , where CP k is the release time of resource k by the project.
The mixed integer programming model for this problem can be formulated as follows:

2.10
In the above model, objective function 2.2 is project cost minimization in which the first and second terms are total costs of using nonrenewable resources and total penalty costs of renewable resources tardiness, respectively.Constraint set 2.3 ensures that each activity j is started under one of its modes within its specified start time periods, that is, EST j , LST j .Constraint set 2.4 forces precedence relationship between activities.Constrain 2.5 limit renewable resource usage.According to constraint 2.6 , renewable resources cannot be used before their ready times and their tardiness periods are determined by constraint 2.7 .Finally, constraint sets 2.8 , 2.9 , and 2.10 are nonfunctional ones.

Particle Swarm Optimization
The PSO algorithm is relatively recent, evolutionary, and population-based metaheuristic, originally developed by Kennedy and Eberhart 15 and redefined by Shi and Eberhart 16 .In spite of the early stages of the PSO method, it has found broad applications in combinatorial and constrained optimization domains and is currently a main research topic.PSO inspired from the social behavior of natural swarms exploits a swarm of particles for search space that are updated from iteration to iteration.Each particle corresponds to a candidate solution assessed by the objective function of the problem in question and is considered as a point in an n-dimension space.The status of a particle is represented by its position and velocity 15 .The n-dimensional position for particle i at iteration t can be represented as X i,t x t i1 , x t i2 , . . ., x t in .Similarly, the velocity is also an n-dimension vector for particle i at iteration t which can be denoted as V i,t v t i1 , v t i2 , . . ., v t in .Using the fitness evaluation, each particle remembers the best position it has perceived so far, referred to as P i,t , and the best position of all particles, the best position of the best particle in the entire swarm, referred to as G t .The position vector of each particle at iteration t is updated using 3.1 , and the particle moves to its new position: X i,t 1 X i,t V i,t .

3.1
Velocity vector is also updated by 3.2 This vector is a function of three components: previous velocity of the particle, the best experience of the particle, also, the entire swarm's best experiences up to the current iteration which are called inertia, cognition part, and social part, respectively 16 .
Updating process continues until the termination criterion is met which usually is the maximal number of generations, processing time, or the best particle position of the whole swarm that cannot improve further after a predefined number of generations.
In 3.2 , r 2 and r 3 are real random numbers with uniform distribution which are usually selected from the interval 0, 1 .c 2 and c 3 are known constants as learning factors, showing the significance of local and global best experiences, respectively.Also, c 1 is defined as a positive inertia weight which was first introduced by Shi and Eberhart 17 .This parameter can be specific for each particle 18 .Liu et al. 19 and Shi and Eberhart 20 introduced time-decreasing inertia weight.
The PSO parameters analyses have been the subject of several researches.For instance, Tchomté and Gourgand 14 determined some conditions for parameters to ensure that each particle converges to some equilibrium point after enough number of iterations.
Although PSO has been originally designed for solving continuous problems, it can be used for solving discrete problems as well.Different techniques have been designed to use this algorithm for combinatorial optimization problems such as the ones introduced by Jarboui et al. 8 , Clerc 21 , and Kennedy and Eberhart 22 .
In this paper, we use the modified PSO approach which was introduced by Tchomté and Gourgand 14 as an extension of PSO that integrates a new displacement rule of the particles.The computational results of their algorithm showed that their PSO algorithm outperformed all state-of-the-art algorithms in solving RCPSP, and this is the reason for selecting their approach for our problem.We describe this modified PSO method in the following.
A metaheuristic algorithm should be able to explore search space effectively and efficiently.PSO algorithm should be intelligent enough to both intensively explore regions of the search space with high-quality solutions and to diversely move to unexplored regions of the search space.These two techniques that were introduced by Glover and Laguna 23 are known as intensification and diversification methods, respectively.Tchomté and Gourgand 14 analyzed particle trajectories and modified particle position updating rules.The idea originated from the PSO applications in which particles basically move from their current positions toward the best local and global positions P i,t , G t , but the particles do not get close enough to P i,t and G t .As a result, diversification is performed well, but intensification is not.Consequently, Tchomté and Gourgand 14 proposed a new particle displacement rule to improve the intensification process by letting each particle visit two positions S i,t and T i,t before moving from current position, X i,t , to the next position X i,t 1 .First, the inertia has influences on the position by making the particle move from X i,t to S i,t .Then the cognition part moves the particle to T i,t and finally under social part affect the particle to reach its new position, X i,t 1 , at the next iteration.Adapted from Tchomté and Gourgand 14 , particle displacement in the classical PSO and this modified PSO has been shown in Figures 1 a and 1 b , respectively.For this purpose, the particle position updating rule is Then, in each time step t, velocity vector is updated as follows:

3.5
Tchomté and Gourgand 14 showed that the necessary conditions for coefficients so that the particles converge to the equilibriums are satisfying 3.6 plus 3.7 or 3.6 plus 3.8 : 3.8

Modified PSO for MRCPSP-TWRTPC
In this section, we present a modified PSO algorithm, using the approach of Tchomté and Gourgand 14 , for solving MRCPSP-TWRTPC.Algorithm 1 shows the pseudocode.In this algorithm the ith particle position at iteration t is represented by the n-dimensional vector

X i,t
x t i1 , x t i2 , . . ., x t in in which x t ij , j 1, 2, 3, . . ., n is the mode assignment to the jth activity at iteration t and is an integer in the interval 1, M j .A feasible schedule of the project is constructed from each X i,t .For this purpose, first the activities are prioritized, see Section 4.3.Then, using single-pass parallel schedule generation scheme, the activities are scheduled, see Section 4.4.Certain local search and improvement procedures are applied on this solution to reach a better schedule, see Sections 4.5, 4.6, 4.7, and 4.8.
Each particle's fitness value is determined by calculating total cost of the final schedule.Then, if necessary, local and/or global best positions are updated, see Section 4.9.If termination criterion is not met, particle positions and velocity vectors are updated by 3.3 and 3.4 , respectively, for the next iteration, see Section 4.10.Different parts of this algorithm are described in more details as follows.

Preprocessing
Sprecher et al. 24 introduced several preprocessing rules in order to reduce feasible space of MRCPSP.Later, these rules have been used in other articles such as Lova et al. 6 , Peteghem and Van Vanhouck 10 , and Hartmann and Briskorn 25 .Considering the similarities between MRCPSP-TWRTPC and MRCPSP, we apply two of these rules to our proposed problem.One is the nonexecutable mode elimination rule for an activity.For a nonexecutable mode, the amount of the resource needed for executing the activity is more than the resource availability.Another method is inefficient mode elimination method.A given mode is inefficient for an activity if there is another mode for which the activity duration is less, and that activity can be accomplished with less total amount of both renewable and nonrenewable resources.Therefore, activities modes are analyzed one by one and nonexecutable and inefficient modes are deleted.

Generating Initial Particle Swarm
Initial particle swarm is generated randomly.Here, component j, j 1, . . ., n, of either position or velocity vectors for each particle is generated randomly from the ranges 1, M j and −M j , M j , respectively, as there is no nonrenewable resource constraint in the problem; moreover, all nonexecutable modes have been deleted, initial mode assignments are feasible, and no modification is needed.

Activity Priority for Scheduling
In order to generate a solution in MRCPSP-TWRTPC, two issues are to be decided: activities mode assignment and scheduling of activities.By specifying mode assignment for a solution, the cost of nonrenewable resources is determined and fixed.Then, scheduling of activities is performed with the objective of minimizing the total cost of renewable resource tardiness penalties.Therefore, the priorities of activities for scheduling are determined by renewable resources.Our procedure for this purpose is as follows.
First, we define the set of activities which need renewable resource k as the activity set of resource k ASR k .Each activity in this set may have immediate or nonimmediate predecessors that may not be a member of this set.We define the set of these predecessors which are not members of ASR k as activities predecessors of resource k APR k .Then, the pairs of ASR k and APR k , k 1, . . ., R, are prioritized by index k using the heuristic that activities in ASR k and APR k for the resource which has more potential of causing tardiness penalty should receive higher priority of being scheduled.To access the potential of the kth resource tardiness penalty cost, we note that this penalty cost is equal to P k × max{0, CP k − d k }, where CP k is the release time of resource k in the project, and hence P k × CP k − d k is a good measure for prioritizing the resources for this purpose.Now the question is how to access CP k without knowing the schedule.We use the following procedure to access CP k , k 1, . . ., R.
Since no activity can start sooner than the ready time of all the resources it needs, in order to take into account the resource ready time, we add one dummy node k for each resource k, k 1, . . ., R, to the project network.Each dummy activity R is single mode with no resource requirement, and the duration of r k .All these dummy activities can be scheduled at time zero.Then, for any activity j which needs renewable resource k, we set dummy activity k as one of its predecessors.So execution of activity j is not possible before the time r k .Note that dummy activity k is a member of APR k .The length of the critical path of subnetwork k is denoted by CP k and is considered as a relative measure of CP k for k 1, . . ., R. CP k is computed using CPM method.After computation of CP k for all resources, the parts of activity sets ASR k and APR k , k 1, . . ., R, are prioritizing using the value of P k × CP k − d k .
In our procedure for prioritizing activities for scheduling, we give more priority to the activities in APR k than activities in ASR k , since activities in ASR k cannot be scheduled unless activities in APR k are scheduled.
Finally, for each resource k, it is necessary to prioritize the activities belonging to each set of ASR k and APR k .We notice when a set of activities using a renewable resource are to be scheduled, we actually deal with an RCPSP, because activities under specified modes are to be scheduled in order to execute within the shortest possible time.Therefore, in order to prioritize activities of a set of ASR k or APR k , we apply one of the most efficient heuristics to scheduling of activities in RCPSP.Lova et al. 6 compared a number of most efficient heuristics for prioritization of activities in RCPSP and found out that activities prioritization based on nondecreasing order of the sum of the latest start and finish time LSTLFT gained the best results among single-pass heuristics.This has been among the best multipass methods as well.Multipass methods need much more computation times than single-pass methods, but they usually result in negligible improvement of the solution; hence, we use the LSTLFT method with single-pass scheduling in order to prioritize activities of ASR k and APR k .

Scheduling Activities
We use parallel schedule generation scheme for scheduling activities, since it fits well with delay local search which we will use.For more details on parallel schedule generation scheme see Section 6.3.2.1.2 of Demeulemeester and Herroelen 1 .
Parallel schedule generation scheme repeats over the separate decision points at which activities can be added to the schedule.A decision point in time horizon corresponds with the completion times of already scheduled activities; thus, at most n decision points need to be considered.At each decision point, the unscheduled activities whose predecessors have been accomplished are taken into consideration in the order of priority list and are scheduled if no resource conflict exists at that time instant.
In this problem, renewable resources have ready time and availability of them differs before and after these times and some new activities may be able to be scheduled after these times.Therefore, we consider ready times in addition to the completion times of activities for choosing new point in time horizon.A new point in time horizon is the closest point to the current point among the ready times of renewable resources and the completion time of scheduled activities.

Delay Local Search
This local search method was used by Chen et al. 26 for the RCPSP problem to escape from local minimum.This method performs like the mutation operator in genetic algorithm and can delay scheduling each activity in spite of its priority to let other activities be scheduled sooner and some resources asked by selected activity will be retained for other activities.
In order to use resources more efficiently, scheduling some activities are delayed in spite of their priority.So other activities can be scheduled sooner.If these selected activities are not delayed, they could delay other project activities for a rather long time.Therefore, each activity is delayed if q ≤ q 0 , where q is randomly selected form uniform distribution 0, 1 and q 0 0 < q 0 < 1 is the probability of delay and called "delay rate."Considering the efficiency of this delay local search in shortening project makespan shown in Chen et al. 26 , we apply this method to scheduling activities.

Mode Assignment Modification-Part I
After scheduling activities, the current schedule may have a set of activities with positive free floats.We call this set FFA.For each j ∈ FFA, it may be possible to change the mode of activity j and reschedule it, using its free float, such that its finish time is not increased and hence no other activity is affected.The mode change and rescheduling of activity j can reduce nonrenewable resource costs, but it can also change CP k , release time of resource k, k 1, . . ., R. The change of CP k can change renewable resource k tardiness, k max{0, CP k −d k }, and its cost, p k l k .If we set the availability of resource k, k 1, . . ., R, for the periods after max{CP k , d k } equal to zero and then reschedule activity j, we are sure that this scheduling is not going to increase renewable resource tardiness penalties.Considering the above points, we define "mode assignment modification-part I" as follows and use it as a local search procedure in our algorithm.
3 For each j ∈ FFA, as we go through forward pass, consider the least nonrenewable resource cost mode of activity j.If it is not its current mode, reschedule activity j in this mode using its free float and considering resource constraints, without increasing its finish time.If this rescheduling is not possible, check the next least nonrenewable resource cost mode of activity j.

Local Left Shift Improvement
Mode assignment modification-part I can reduce the renewable resource requirements of the project for certain periods under the resulted schedule.This should be expected since usually the mode with less nonrenewable resources has longer duration and requires less renewable resources per period.Hence, the possibility of local left shift of certain activities exists.Therefore, we perform the standard local left shift method after the mode assignment modification-part I.
After performing local left shift, we might be able to modify some of the activities mode assignment again by mode assignment modification-part I. Hence, these two procedures are used one after another until no improvement can be achieved in the schedule.

Mode Assignment Modification-Part II
We consider the set of project activities with no successors and call this set NSA.The direct predecessor activities of dummy activity n make NSA.For any j ∈ NSA if we change the mode of activity j and reschedule it, the schedule of no other activity changes and the value of cost change of the project is Δc ΔNRC j R k 1 p k Δl k , where ΔNRC j is the change of nonrenewable resource cost of mode for activity j and Δl k is the change of kth resource tardiness.Since activity j has no successor, Δl k can be computed easily.If the value of change effect on total cost is negative, the mode change for activity j ∈ NSA is justifiable.Considering above points, we define the following local search procedure as "mode assignment modificationpart II" and use it in our algorithm.
1 Let NSA {j | activity j is the direct predecessor of dummy activity n}.
2 For each j ∈ NSA, consider some modes for activity j, in which nonrenewable resources requirement cost is less and this cost saving is more than the probable increase in the penalty cost of renewable resources.Compute Δc for each of them.
Considering renewable resource constraints, if the least Δc is negative, its corresponding mode replaces the current mode of activity j.

Updating the Local and Global Best Solutions
As mentioned earlier, the PSO algorithm stores the best local solution obtained by each particle and the best global solution obtained by the entire swarm.Therefore, after evaluating all the particles of the swarm at iteration, the best local solution for each particle up to the current iteration is compared with the current solution of the particle.If the current solution has lower total cost, the best local solution of the particle is replaced by the current solution of the particle.Similarly, the best global solution is updated.

Updating Particles Position and Velocity
To update position of the particles to generate new solutions for the next iteration, firstly, particles velocity is updated using 3.4 .In this process each element of the velocity vector which is more than M i is changed to M i and each element of the velocity vector which is less than −M i is changed to −M i .Then, particles positions are updated using 3.3 .Similar to the procedure for updating each particle velocity, each element of the new position vector which is more than M i is changed to M i and each element of the new position vector which is less than 1 is changed to 1.As the elements of position vector determine the mode assignment of activities and should be integer, each noninteger element is rounded to the closest integer.Each test instance was run for all 3 5 × 2 486 permutations of parameter values, a total of 48,600 problems.The tuned values of the parameters are q 0 0.4, P 60, c 1 0.2, c 2 0.5, c 3 1.

Algorithm Validity
As our research is new, we could not find any solved problem in the literature.So we developed some instance problems whose optimum objective function values are in hand.Then we solved these instances with our algorithm and tested the results.In order to generate these instances, we used sample problems introduced in Section 5.1; however, we modified the due dates of renewable resources as follows.We generated a random feasible schedule for each instance, after assigning the least nonrenewable cost mode to each activity of the project.In this schedule, we determined the release time of each renewable resource.Subsequently, we set the due date for each renewable resource equal to its release time; hence, this schedule has zero tardiness penalty cost and its objective function value is equal to the cost of nonrenewable resources.As we assigned the least nonrenewable cost modes to the activities, this schedule is optimal and the optimum value of its objective function is available.
All sample problems were modified with the above procedure and solved with the PSO algorithm.The termination criterion was generation of 600 schedules.In order to assess the validity of the algorithm, d, the percent deviation of the objective function value from optimum, computed for each test problem solved, where d 100 × Z p − Z opt /Z opt , Z p is the objective function value of the best solution achieved by the PSO algorithm and Z opt , the optimal objective function value of the instance.Table 2 shows the average and standard deviation of d for each problem set.The low values of average and standard deviation of d reveal good performance of the algorithm with small CPU time, although we do not have any standard for comparison.
Table 5 shows the effect of deleting mode assignment improvement method-part II.We can see that in most cases the performance of the algorithm deteriorates, and in all cases the average CPU time increases remarkably.In the following, we explain the reason for increasing the average CPU time.

Conclusions
In this paper, we introduced MRCPSP-TWRTPC problem as a resource-oriented cost minimization project scheduling problem considering both renewable and nonrenewable resource costs.We formulated and mathematically modeled this problem as mixed integer programming model and discussed its NP-hardness.Subsequently, we developed a metaheuristic algorithm to tackle the proposed project scheduling problem.We briefly reviewed the applications of the PSO algorithm for solving combinatorial and constrained optimization problems.Thereafter, we applied a modified PSO algorithm including modified updating rules for particles velocity and position.In order to generate feasible schedules, we used the PSO algorithm for activity mode assignment and developed a novel heuristic technique to prioritize activities for parallel scheduling scheme.Two improvement heuristics, delay local search and local left shift, in line with two mode assignment modification methods, were implemented to improve the solutions.The computational results revealed proper algorithm robustness in solving different instances especially with high number of iterations.Also, the validity analysis showed small deviations from the optimal solutions for the test instances in reasonable solving time.Finally, we assessed two improvement methods used in our algorithm to demonstrate their good performance.

Figure 1 :
Figure 1: Particle displacement in the swarm: a classical PSO, b modified PSO.