Multimode Preemptive Resource Investment Problem Subject to Due Dates for Activities : Formulation and Solution Procedure

The preemptive Multimode resource investment problem is investigated. The Objective is to minimize the total renewable/nonrenewable resource costs and earliness-tardiness costs by a given project deadline and due dates for activities. In this problem setting preemption is allowed with no setup cost or time. The project contains activities interrelated by finish-start type precedence relations with a time lag of zero, which require a set of renewable and nonrenewable resources. The problem formed in this way is an NP-hard. A mixed integer programming formulation is proposed for the problem and parameters tuned genetic algorithm (GA) is proposed to solve it. To evaluate the performance of the proposed algorithm, 120 test problems are used. Comparative statistical results reveal that the proposed GA is efficient and effective in terms of the objective function and computational times.


Introduction
The resource constrained project scheduling problem (RCPSP) is a known combinatorial optimization problem which is NP-hard [1].The objective of RCPSP is to minimize the makespan of the project while the renewable resources availabilities are considered given.In the literature there are several exact methods and heuristics that solve the RCPSP [2][3][4][5][6][7][8][9].
In the multimode version of the RCPSP (MRCPSP), each activity can be performed in one out of a set of modes, with a specific activity duration and resource requirements.The problem involves the selection of a mode for each activity and the determination of the activity start or finish times such that project makespan is minimized.However, the precedence and resource constraints should be satisfied.MRCPSP is also NP-hard as generalization of the RCPSP.In recent years several algorithms have been proposed to solve MRCPSP [10][11][12][13][14][15][16][17][18][19][20].
The resource investment problem (RIP) is a close variant of RCPSP.This problem consists of determination of the activities start times and renewable resources availabilities such that the total cost of the resources is minimized subject to a given project deadline.In RIP, project's makespan is not forgotten but it is controlled with a predetermined deadline as a constraint.This problem was introduced by Möhring [21].He showed that the problem is NP-hard.Rangaswamy [22] proposed a B&B for the RIP and applied it to the same instance set used by Demeulemeester [23].Drexl and Kimms proposed two lower bound procedures for the RIP based on Lagrangian relaxation and column generation methods [24].Other exact procedures are proposed by Demeulemeester [23] and Rodrigues and Yamashita [25].Yamashita et al. [26] proposed a multistart heuristic based on the scatter search.Shadrokh and Kianfar presented a genetic algorithm (GA) to solve the RIP when tardiness of the project is permitted with penalty [27].Ranjbar et al. developed a path relinking procedure and a genetic algorithm (GA) [28].Van Peteghem and Vanhoucke presented an artificial immune system algorithm for the problem [29].
The basic project scheduling problems assume that each activity, once started, will be executed until its completion.

Advances in Operations Research
Preemptive project scheduling problem refers to the scheduling problem which allows activities to be preempted at any time instance and restarted later on at no additional cost or time.In the case in which activities are preemptive, modeling and solving such a problem as a classical nonpreemptive RCPSP or RIP may lead to poor solutions.Such situations where preemption of an activity is beneficial or necessary are typical, for example, in industry processes where processing units, like reactors or filters, have to be cleaned after the completion of certain subactivities.In the textile industry, preemptions are inevitable.When the fabric type is changed on a machine, the wrap chain must be replaced which indicates the necessity of preemption.The literature on solution methods for the preemptive project scheduling problems is relatively scant.For the preemptive RCPSP, we refer to [30][31][32][33][34].For the preemptive MRCPSP, Buddhakulsomsiri and Kim proved that preemption is very effective to improve the optimal project makespan in the presence of resource vacations and temporary resource unavailability [35].Van Peteghem and Vanhoucke have proposed a genetic algorithm for the MRCPSP and its extension to the preempted case [36].
Minimization of the earliness-tardiness penalties of the project activities is a nonregular performance measure in just-in-time (JIT) environments.In this problem, activities have an individual due date with associated unit earliness and unit tardiness penalty.Traditionally, tardiness penalties due to delivery after a contractually arranged due date are considered.Tardiness leads to customer complaints, loss of reputation and profits, monetary penalties, or goodwill damages.Costs of earliness include extra storage requirements, idle times, implicitly incur opportunity costs, and storage costs due to insurance, theft, perishing, and bounded capital for the case of an activity is completed before the due date.Therefore, the minimization of the earliness-tardiness (E/T) costs is attractive in order to meet the requests of practice.For some recent solution procedures considering earlinesstardiness penalties we refer to [37][38][39][40].
There are some shortcomings in the basic RIP which is considered in this paper simultaneously.First, activities in the basic RIP are assumed nonpreemptive, while this assumption is not true in practice.Second, in the basic RIP, the determination of only renewable resources availabilities is investigated.Third, the basic RIP supposes single execution mode for activities.Last, due dates for activities are neglected.Therefore, the contribution of this paper is threefold: first, a mixed integer programming formulation is developed for the preemptive multimode RIP with due dates for the activities.We call this problem P-MRIP-WET.This model is not considered in the past literature.Second, a new efficient parameter-tuned GA is developed for the problem due to NP-hardness of the problem.Finally, the effectiveness of the proposed method to solve the P-MRIP-WET is analyzed statistically.
The reminder of the paper is organized as follows.Section 2 describes the problem P-MRIP-WET and mixed integer formulation for it.Section 3 explains the steps of proposed genetic algorithm (GA) to solve the problem.Section 4 contains the computational results and performance evaluation of proposed GA.Finally, Section 5 concludes the paper.

Problem Description
The preemptive multimode resource investment problem with weighted earliness-tardiness penalties (P-MRIP-WET) involves the scheduling of project activities on a set  of renewable resource types and a set   of nonrenewable resource types.Each activity  is performed in a mode   , which is chosen out of a set of   different execution modes, that is, with different durations and resource requirements.The duration of activity , when executed in mode   , is    .Each activity  in mode   requires      units of renewable resource type  ( = 1, . . ., ) during each time unit of its execution.For each renewable resource , the availability    is constant throughout the project horizon.Activity , executed in mode   , will also use  ] nonrenewable resource units (  = 1, . . .,   ) of the total available nonrenewable resource  ]   .In sequent, assume a project represented in AON format by a directed graph  = {,} where the set of nodes, , represents activities and the set of arcs, , represents finish-start precedence constraints with a time lag of zero.The preemptable activities are numbered from the dummy start activity 0 to the dummy end activity  + 1 and are topologically ordered; that is, each successor of an activity has a larger activity number than the activity itself.Also, we consider discrete time points for the preemption.As a rule, a preemptive problem is characterized by a complicated structure of its optimal solutions.When preemption is allowed at arbitrary times, the problem turns out to be intractable [41].In this situation, when the overall number of preemptions is unlimited and the set of admissible points for preemptions has continuous cardinality, we cannot utilize exact enumerative algorithms, unless a nontrivial preliminary analysis of properties of optimal solutions is performed.Such analysis reduces the original infinite set of possible points for preemption to a finite set.This reduction allows us to solve the problem by direct enumeration.This analysis is performed for some scheduling problems.Baptiste et al. prove that a wide class of preemptive scheduling models, including both machine and project scheduling models, have the "integer preemption property"; for any problem instance with integral input data, there exists an optimal schedule where all preemptions (as well as starting and completion times of jobs/activities) occur at integer time points.This conclusion is held for objective functions such as total weighted earliness-tardiness and total weighted number of late jobs [42,43].
The objective of the P-MRIP-WET is to schedule a number of activities, in order to minimize the total cost of the resources availabilities and earliness-tardiness penalties.A schedule  is defined by a vector of activity start times and is said to be feasible if all precedence relations, project deadline, and renewable and nonrenewable resource constraints are satisfied.However, execution modes   for activities, available resource capacities, and start times for activities have to be determined.We have the Notations section for P-MRIP-WET.
The variables     can only be defined over the time interval of the activity in question  ∈ [EST  , LFT  − 1].These limits are determined using the traditional forward and backward pass calculations considering duration of activity  based on the execution mode with lowest duration.The backward pass calculation is started from a fixed project deadline DL.In this paper, earliest finish time of dummy end activity EFT +1 multiplied by 1.3 is considered as project deadline.EFT +1 is computed using the traditional forward calculations considering duration of activity  based on execution mode with highest duration.
It is clear that an activity with duration of 0 is never in progress and thus does not have a corresponding decision variable which is set to 1.This problem, however, can be easily overcome: the dummy start and end activity are assigned a dummy mode with duration of 1. Also, the other parameters for dummy modes are assumed 0. All other activities with zero duration can be eliminated, provided that the corresponding precedence relations are adjusted appropriately.The resulting schedule may be transferred into a schedule for the original problem by removing the dummy start and end activity, and one time unit left shifting.
Using the above notation, the P-MRIP-WET can be mathematically formulated as follows: The objective in (1) is to minimize the total cost of the project resources and earliness-tardiness penalties.First part is weighted earliness-tardiness cost.Second part is costs related to renewable resources which are available and renewable from period to period.Only the total resource use at every time instant is constrained.Third part is costs related to nonrenewable resources which are available on a total project basis, with limited consumption availability for the entire project.Consider the following: Equation ( 2) guarantees that each activity  can be in progress with at most one mode in each time unit: Equation ( 3) ensures that each activity  should be in progress for    time units in assigned mode   : Equation ( 4) specifies that only one execution mode is allowed for each activity : Constraint set in (5) takes care of the renewable resources limitations: Constraint set in ( 6) takes care of the nonrenewable resources limitations: Constraint set in (7) computes the start time of first time unit of activity : Constraint set in (8) computes the start time of last time unit of activity : Equation ( 9) denotes the finish to start zero time lag precedence relations constraints: Equation ( 10) computes the earliness for each activity : Equation ( 11) computes the tardiness for each activity : Constraint in (12) guarantees that project deadline is not violated:

Proposed GA to Solve P-MRIP-WET
In this section we describe the genetic algorithm (GA) which is the well-known metaheuristic that has been successfully applied to a noticeable number of project scheduling problems [6,15,18,19,27,36,38], to solve P-MRIP-WET.In order to increase quality of the proposed GA, we implement two efficient local searches.We also consider results obtained from CPLEX software version 12.3 and randomly generated feasible solutions, to provide comparable computational efforts for the proposed GA.Solution representation, selection, and reproduction are the basic elements of GA.These elements must be well defined and adapted to a specific problem.Before the execution of the proposed GA all nonexecutable and inefficient modes can be omitted in order to reduce the search space [44].An execution mode   is called nonexecutable if its execution would violate the renewable resource constraints in any schedule.Also, a mode is called inefficient if there is another mode of the same activity with the same or higher duration and no more requirements for all resources.
3.1.Solution Representation.Two important representations for schedules are the random-key (RK) and the activity-list (AL) representation.It is deduced from experimental tests that procedures based on (AL) representation outperform the other procedures [45].In this study the AL representation is used to encode a project schedule and the parallel schedule generation scheme (PSGS) to decode the schedule representation to a schedule.We represent a feasible solution by an  +   +  +   element vector ().
The first  element represents execution modes for activities (mode list).In the sequel it is assumed that each original activity  with assigned mode   is replaced with    activities with unit duration.Each duration unit  = 1, . . .,    of an activity  is predecessor of duration unit (−1).Consequently, project new network has   = ∑  =1    activities with duration 1 and same resource requirements as original activities.Only duration unit  th   has fixed due date DD  .Also, without loss of generality, for duration units  = 1, . . ., (   − 1) of an activity  due dates are assumed 0 with no earliness-tardiness penalty.
The second   element represents a precedence-feasible permutation of activities (activity list).The third  element is renewable resource capacities (resource list).The last   element is nonrenewable resource capacities which is determined based on assigned modes in mode list.Consider the following: Having obtained a feasible solution represented by the vector described above, the starting times of all activities are defined by parallel SGS.The SGS determines how a feasible schedule is constructed by assigning a starting time to the activities with respect to the assigned modes and the resource availabilities.The PSGS sequentially iterates over the different decision points at which activities can be added to the schedule until a feasible complete schedule is obtained.At each decision point, the unscheduled activities whose predecessors have completed are considered (in the order of the activity list) and are scheduled with assigned modes (specified in the mode list) such that resource availabilities (specified in the resource list) are not violated.

Initial Generation.
The GA starts with producing the initial generation.As used in the literature [19], the initial generation can be created randomly.In proposed GA, each solution in initial generation is randomly created as follows.
First, we assign a randomly selected mode to each activity .
Second, we generate an empty   = ∑  =1    elements vector.To fill this vector an activity is randomly selected from set of activities that are not selected before and all their precedence are already met.This process is repeated until the activity list is full.Third, renewable resource  availability is randomly selected between lower and upper bound of    .Finally, the nonrenewable resource capacity   is computed based on assigned modes; that is, summation of resource requirement of the activities in the assigned modes for nonrenewable resource   is considered as  ]   .The capacity of the renewable resource type  cannot be less than    = Max  {Min   (     )}, while the upper bound of the capacity of the renewable resource type ,    , can be determined by unconstrained scheduling of activities at their earliest start times.

Reproduction.
Once the initial generation is created, the process continues in an iterative way in order to obtain the next high quality generation.The better solutions of the generation will have higher possibilities to survive for the reproduction.In proposed GA, the reproduction process consists of four phases.In the first one, the codified segments of two selected solutions are exchanged following the crossover operators.In the second phase, the selected solutions are randomly changed using the mutation operators.In the third phase, some improved solutions created by knowledge based local searches are inserted to the next generation.In the last phase, the current generation is sorted by best to worst fitness function.Then   percent top solutions from the current generation are copied to the next generation.In this way, these operators lead to the exploitation of good solutions and to the exploration of new areas of search space.
The selected solutions for crossover and mutation operators are called parents.The parents are selected based on their fitness values.We used the roulette wheel scheme, where each solution has a probability of selection that is directly proportional to its fitness value.The selected solutions are grouped in couples randomly.Then the crossover operators act with a probability   for each couple, and the mutation operators act on each child created by the crossover operators with a probability   .However,  LS percent of improved solutions created by two local searches are inserted to the next generation.In order to preserve the constant population size, we consider   = (1 −   −   −  LS ).

Crossover Operators.
Crossover operators, depend on the type of solution representation.In the proposed GA, one of the following operators randomly acts on the selected couple.
(i) Two integers V and  are randomly selected from set {1, . . ., }.Then the execution modes of original activities between V and  are crossed to form other two solutions.Subsequently, the activity list and resource capacities are updated.
(ii) Two integers  and  are randomly selected from set {1, . . ., }.Then the renewable resource availabilities  and  are crossed to form other two solutions.
(iii) Two integers  and  are randomly selected from set {1, . . .,   }.Then the activity list of parent solutions are divided in three parts and they are combined to form two news solutions (Figure 1).A, B, and C are three parts (subset of activities) of parent 1 and D, E, and F are three parts (subset of activities) of parent 2. First part (D) of child solution 1 is copied from parent solution 2. The next parts (a and b) of child 1 are taken from two first parts of parent solution 1 (A and B) in the same order.However, the activities that have already been added from the parent 2 may not be considered again.Then the nonrepetitive activities in C from parent solution 1 which have not appeared in F from parent solution 2 are added to part c of child 1.Finally, the activities that have not existed in child 1 are taken from the part F of parent 2 in the same order (f).The activity list of child 2 is produced similarly from parent 1 and parent 2.This operator guarantees the feasibility of the new children solutions because, as it is clear in Figure 1, the priority of activities in child 1 is D-a-b-c-f which preserves the priority of activities in parent 1 (A-B-C) and parent 2 (D-E-F).
Similarly, the priority of activities in child 2 is A-d-ef-c which preserves the priority of activities in parent 1 (A-B-C) and parent 2 (D-E-F).

Mutation Operators.
In the proposed GA, one of the following operators randomly acts on the selected child created with the crossover operations.
Then an execution mode different from the current mode of original activity   is assigned to it.Subsequently, the activity list and resource capacities are updated.
Then two predecessors of activity   are randomly selected and their positions in activity list are exchanged.Also, some activities between these positions are shifted to left or right such that feasibility of resulting solution is preserved.knowledge based local searches are employed on selected parents.
( , a different resource is randomly selected.This procedure is continued until the fitness of new solution is better than the current solution (and project deadline is not violated) or all resources are equal to lower (upper) bounds.The idea behind this operator is to balance between earliness-tardiness cost and the cost of the renewable resources.In doing so, a renewable resource is randomly selected and reduced (increased) with a probability which has direct (inverse) relation with the percentage of project activities which are on time or have earliness.Clearly, a low amount of resource availabilities leads to more tardiness and less earliness and resources cost.

(ii) Activity Based Local Search (LS2)
. This local search is applied if at least one activity with earliness exists in current solution.All activities with earliness are identified and their start times are calculated by parallel SGS.Let   and   denote the start time of first and last activities which have earliness, respectively.A random integer  is selected from set {  , . . .,   − 1}.Suppose   denotes the first time interval larger than  which is empty without any in-progress activity.Time interval [,  + 1] is left empty without any in-progress activity and all activities between  and   are shifted to right one unit time such that   is full.The fitness of new solution is evaluated.If it is better than the current solution and project deadline is not violated, it replaces the current one.

Parameters Tuning.
The performance of metaheuristics depends excessively on the value of their parameters.In this paper the Taguchi experimental design is used to tune the parameters of GA.Taguchi divides factors into controllable and noise factors and offers a set of orthogonal arrays for designing experiments of quality improvement.Although there is no direct control of noise factors, the Taguchi method determines the optimal level of controllable factors and minimizes the effect of noise [46].In the proposed GA, the factors that should be tuned are Pop, It,   ,   , and  LS , where Pop is the population size, It is the number of iterations (generations), and   ,   , and  LS are the crossover, mutation, and local search probabilities.A set of 27 randomly generated problems with 20 nondummy activities are used for parameter tuning.Using MINITAB software version 16, based on a L27 orthogonal array design, the average S/N ratio is obtained at each level (Figure 2).Different levels of these factors and the optimal levels (in bold) of Pop, It,   ,   , and  LS are shown in Table 1.

Stopping Criterion.
The procedure is continued until a predetermined number of generations are produced.We obtained good results by indexing the number of generations to the size of the problem, that is, use of the small number of generations for small problems and large number of generations for larger problems.So, the tuned values of Pop and It (90, 30) are used for problems with about 20 activities.However, these values are increased proportionally for larger sizes.

Performance Evaluation
4.1.Test Problems.In order to validate the proposed GA algorithm, a set of 120 problems was generated by the generator ProGen developed by Drexl et al. [47], using the parameters given in Table 2.
For each combination of problem parameters 5 problems were generated.The deadlines were generated in the same The due dates were generated in the same way as Vanhoucke et al. [37].First, a maximum due date was obtained by multiplying the critical path length by 1.5 for each project.Then, random numbers generated between 1 and maximum due date.The numbers are sorted and assigned to the activities in increasing order.

Experimental Results
. The proposed GA were coded in Borland C++ 5.02 and executed on a personal computer with an Intel Core2Dou, 2.5 GHz processor, and 3 GB memory.Table 3 present the computational results of the proposed algorithm.For problems with 20 subactivities, the results are compared with the optimal solutions obtained by CPLEX.However, for problems with 50 and 80 subactivities which CPLEX is unable to solve, the results are compared with the best randomly generated solution (BRGS).Proposed GA ARD-GA: average relative deviation percentages for the GA; ARD-BRGS: average relative deviation percentages for the BRGS.
Relative deviation (RD) percentage for each problem is obtained by the following formula: where  is the value of objective function and  * is the optimal solution obtained by CPLEX.In the cases for which CPLEX is not able to solve the problem optimally,  * is the best obtained solution by the GA or the best randomly generated solution.Table 3 reveals that when the number of subactivities is equal to 20, all 40 problems can be solved optimally by CPLEX 12.3 within 1000 second.It is noticeable that when the number of subactivities is equal to 20, the mathematical model contains averagely 57 variables and 139 constraints.Average CPU-time for CPLEX indicates that when the number of resource types is increased the complexity of problem is increased.Also, Table 3 shows that when the number of subactivities is greater than 20, while the CPLEX is unable to solve the problem, there is a solution by the proposed GA in a satisfying CPU time.ARD for the GA algorithm are not high.This means that proposed GA gives robust solutions.Also, NPMG for the GA indicates that too many of executions of problems with 20 subactivities reach to the optimum solution.Also, Table 3 shows that best randomly generated solutions (BRGS) conquest GA only in 1.75% runs (21 out of 1200 runs).Also, Table 3 indicates that solutions obtained by BRGS are not reliable measured by ARD.Consequently, these results confirm the advantages of the proposed GA to solve the P-MRIP-WET.

Sensitivity Analysis on Local Searches.
In order to evaluate the effect of the knowledge based local searches on speed and quality of GA a sensitivity analysis is done.The results are demonstrated in Figures 3 and 4. Figure 3 shows that, by simultaneously using LS1 and LS2, the proposed GA has best quality measured by objective function.While, by only using LS1 and LS2, objective function will be 1.79% and 3.57% worse, respectively.However, applying GA without any local searches leads to solutions about 6.94% worse than simultaneously using LS1 and LS2.It is clear from Figure 4, by simultaneously using LS1 and LS2, that the proposed GA has lowest speed measured by CPU time.By only using LS1 speed of algorithm reduces slightly, while by only using LS2 algorithm speeds up considerably (about 18%).Also, applying GA without any local searches leads to about 55% more speed than simultaneously using LS1 and LS2.

Summary and Conclusions
In this paper, we attempted to solve the preemptive multimode resource availability cost problem with due dates for the activities (P-MRIP-WET).The objective is to schedule the activities in order to minimize the total cost of the resources availabilities and earliness-tardiness penalties subject to the precedence constraints and a fixed deadline on project.This problem has not been studied ever before.The problem was described with an integer programming model, and then the genetic algorithm (GA) was proposed to solve it.However, two knowledge based local searches are proposed to improve the evolution speed of GA.The parameters of proposed GA are tuned based on Taguchi experimental design.The performance of the proposed algorithm on 120 test problems was compared with the results of the CPLEX and best randomly generated solutions (BRGS).From the computation results, we could clearly see that the proposed GA could efficiently solve the project scheduling problem.

Figure 2 :
Figure 2: The mean S/N ratio plot for each level of the factors.

Figure 3 :Figure 4 :
Figure 3: The effect of local searches on quality of GA.
1 ,   2 , . . .,    ) , ( 1 ,  2 , . . ., i) Resource Based Local Search (LS1).The percentage of project activities which are on time or have earliness (denoted by ) is computed.A random real number (denoted by ) is selected from interval (0, 1).Then, an integer  is randomly selected from set {1, . . ., }.If  ≤  the availability of renewable resource  is reduced (else it is increased) one unit.If the current availability of renewable resource  is equal to lower (upper) bound of

Table 1 :
Factors and factors levels.

Table 2 :
The parameter settings for the problem set.

Table 3 :
Comparison of the results obtained by the CPLEX, GA, and BRGS.
Number of renewable resource(s), index by    : Number of nonrenewable resource(s), index by     : Set of execution modes for activity ,  ∈  |  |: Number of execution modes for activity , index by      : Duration of activity  in mode   ,  ∈ ,   ∈        : Resource requirement of activity  in mode   for renewable resource type ,  = 1, . . ., ,  ∈ ,   ∈    ]     : Resource requirement of activity  in mode   for nonrenewable resource type   ,   = 1, . . .,   ,  ∈ ,   ∈   DD  : Due date of activity ,  ∈  DL: Project deadline   : Per unit earliness cost of activity ,  ∈    : Per unit tardiness cost of activity ,  ∈     : Unitcostofrenewableresourcetype,  = 1, . . .,   ]   : Unit cost of nonrenewable resource type   ,   = 1, . . .,   EST  : Earliest start time of activity  (computed by assuming infinite resource capacities) LFT  : Latest finish time of activity  (computed by assuming infinite resource capacities) : Objective function (total cost of the resources availabilities and earliness-tardiness penalties)    : Constant availability of renewable resource type  throughout the project horizon,  = 1, . . .,  (integer decision variable)  ]   : Totalavailabilityofnonrenewableresource type   ,   = 1, . . .,   (integer decision variable)   : Earliness of activity  (integer decision variable)   : Tardiness of activity  (integer decision variable)   Min : Start time of first time unit of activity  (integer decision variable)   Max : Start time of last time unit of activity  (integer decision variable)     : 1, if activity  in mode   is in progress at time interval [,  + 1], 0, otherwise (binary decision variable)    : 1, if activity  is executed in mode   , 0, otherwise (binary decision variable).