A PARALLEL SIMULATED ANNEALING ( PSA ) FOR SOLVING PROJECT SCHEDULING PROBLEM WITH DISCOUNTED CASH FLOW POLICY IN PRICING STRATEGY OF THE PROJECT SUPPLIERS

Original scientific paper Resource-constrained project scheduling problem is known as a NP-Hard problem in literature. In this research, discounted cash flow policy is suggested for the resources-constrained project scheduling problem for the first time while in classic models, it has been assumed that price of the required resources is fixed for performing the activities and resources can be prepared only with one price rate in the market. Goal of this problem is to determine optimal starting time of the project activities considering precedence constraints and the available resources such that the project completion time can be minimized. In order to solve the proposed model, a hybrid algorithm based on two algorithms i.e. genetic and simulated annealing has been suggested. In this method, genetic algorithm has been designed as the main framework of the proposed method and simulated annealing method as a new operator and in order to improve local search of the main algorithm. Since values of the parameters have considerable effect on efficiency of these algorithms, therefore, a new statistical approach based on the stepwise regression has been presented to set the proposed algorithms parameters. Results of the calculations show high efficiency of proposed algorithm in terms of solution time and optimal solutions.


Introduction
One of the applicable research problems in management and programming of the projects is the resources-constrained project scheduling which attracted attention of many researchers.This problem is applied broadly in fields of management of the industrial projects, construction engineering, software development as well as production scheduling [1,2].Main goal of this problem is to schedule a project in order to minimize completion time of the entire project such that precedence constraint and the available resources are observed.It is assumed that the project activities are independent of each other and completion time of each one of the activities is dependent on the resources allocated to that activity.The available resources are constrained and preemption of the activities is not permissible.
Relationship between the spent resources for the activities and completion time of activities is one of the factors which form the project scheduling problem.Extent to which the resources are effective on the required time for completion of the activities can be presented as time function in terms of resources.It is evident that no comprehensive and equal form can be considered for structure and behavior of these functions but the factors such as nature of the available activity, methods and technologies, access to more resources etc. form time function of each activity.Different functions have been presented in the literature for showing relationship between time and resources including continuous linear, discrete, convex nonlinear and concave nonlinear functions.Special functions in resource -constrained project scheduling problem are multimode functions in which any activity can be programmed in one of many options in terms of relationship between completion time and the required resources.In this problem, determination of optimal time of activities and the allocated resources to them results in a balance between completion time and total project cost.The goal is to determine allocation of the states optimality to the activities such that the precedence constraints and the available resources are not violated and completion time of the project can be minimized.Blazevic et al. [1] show that the projects scheduling problem as a general case of the job-shop production belongs to NP-Hard optimization problems class for which it is difficult to calculate good approximate answers.On the other hand, it is impossible to find an efficient algorithm to solve the optimality of samples with large dimensions from this problem and at reasonable solution time.Some of the useful papers which review projects scheduling problem solving resource-constrained techniques have been presented by Heroen et al. [3].Like most complex optimization problems, the solution approaches presented in the literature include board scope of the precise approaches, heuristic methods and metaheuristic algorithms.Many precise approaches have been presented in the literature for these problems [4÷6].Some researchers mentioned their research field to calculate a lower bound for the mentioned problem.
Brooker and Nest [7] suggest lower bound for the problem based on the given value since project completion time.A hybrid approach has been presented by Demasi et al. [8] based on constrained programming and integer programming [8].Most heuristic algorithms for the intended problem use priority rules to design the answer constructive methods.These methods start with a subset of activities and schedule a new activity at any stage until all activities are programmed.Bokter [9] studied 21 heuristic scheduling methods to solve the resource-constrained project scheduling problem and finally mentioned their combination as the new heuristic algorithm.In addition, Bokter has suggested more precise heuristic algorithm using critical path method and its calculations.In recent years, Notes et al. [11] studied agent based algorithms for the mentioned problem and Lova et al. [12] presented methods based on the priority rules to solve that.Most heuristic algorithms which have been presented use two serial and parallel productions methods of scheduling.
A comprehensive review of these two approaches has been presented to solve the project scheduling problems.The last class of the presented solution approaches for the project scheduling problem included application of the metaheuristic algorithms.These methods start with an initial answer and search for the optimal answer with a set of answers.
Different metaheuristic algorithms have been suggested in the literature for the projects scheduling problem.Genetic algorithm is one of the known methods which have been presented by [14÷16].Simulated annealing [17], tabu search [18], algorithms based on particle swarm optimization [19,20] are the other metaheuristic methods which have been suggested in the literature to solve this problem.Also Dalfard and Ranjbar [21] used simulated annealing algorithm for Multiprojects scheduling with resource constraints and Priority rules.Renato and Silva [22] developed a hybrid method for a project scheduling problem.The other new researched are included Pieter Leyman and Vanhoucke that they developed a new scheduling technique for the resource-constrained project scheduling problem with discounted cash flows [23].Joseph and Venkateshan develped an integer programming formulation for the project scheduling problem with irregular time-cost tradeoffs [24].
As mentioned above, all previous researches on the resource -constrained project scheduling considered fixed price for providing any required resources in their proposed models while in some projects, the project managers face discount request for purchasing the resources.In this case, the project managers and programmers should evaluate and decide allocation of resources to the activities based on the available offer.Since the available models have considered only fixed price for purchasing the resources from suppliers, therefore, they cannot satisfy such conditions.
In this research, we suggest Discounted cash flow policy with resource-constrained projects scheduling problem in which effect of different policies of the suppliers on the project scheduling is studied.In order to solve the proposed problem in this research, a hybrid approach based on genetic algorithm has been presented in which a new operator has been used as local search.Main body of the algorithm is composed of genetic algorithm operators and simulated annealing algorithm has been used.At the end, a new method based on stepwise regression has been used to set the algorithm parameters in order to improve performance of the proposed approach.
When we face many independent variables, it is very important to determine suitable combination of these parameters for presence in the model.This subject will be more important when the used algorithm is affected by many parameters.Calculation results show suitable performance of the proposed approach versus other approaches.
Section 2 relates to definition of the problem and its mathematical formulation.The proposed solution based on genetic algorithm and new hybrid operator is presented in section 3. Results of tests design and comparative analysis of the proposed method are shown in section 4. At the end, section 5 includes conclusion.

Problem mathematical formulation
In this research, we consider a resource -constrained project scheduling problem in which completion time of the activities is a multi-mode discrete function.This problem is defined as follows.The project is shown as an activity on arrow (AOA) and has 2 independent activities.In order to show precedence relation between activities of the project, an acyclic graph G is used.
Successor activities can start when all precedence activities have been completed.We show completion time of the i th activity with   , k-type resource which is allocated to this activity with   decision variable and completion time of activity with   ( = 1, 2, … ., ).It has been also assumed that the number of required resources equals  and   resource is accessible in all programming horizon.Two kinds of constraint are defined for activities: (1) precedence constraint which results in the i th activity starts when all of precedence activities of   are completed and (2) constrained use of the available resources which guarantees that total resource k consumed by the activities does not exceed accessible resource   .Basic model of the project scheduling problem with the obtained resources is presented.

min 𝐹𝐹 𝑛𝑛+1
(1) where target function is minimization of the project completion time which has been mentioned by the 2+1 th virtual activity completion time.Series of the first constraints (1.1) ensure that the precedence relations are established between activities of the project while the second type constraints (2.1) apply the resources constraints.In order to determine   function in terms of the allocated resources to the i th activity, it has been assumed that each one of the activities can be programmed in one of two cases of  1 and  2 .Case M 1 indicates allocation of the highest value of resources to the activity and consequently the shortest completion time and case  2 indicates allocation of the lowest value of resources to the activity and consequently the longest completion time.In order to model this function, we assume that d im indicates completion time for the i th activity and  , specifies k resource allocated to the i th activity while this activity is programmed in   , ( = 1, 2).In this case, main problem (1) can be reformulated as an integer programming problem.

min 𝐹𝐹 𝑛𝑛+1
(2) Goal of the above model ( 2) is to determine optimal starting time of the project activities based on optimal level of the resources allocated to them so that total time of project can be minimized.In the classic models of the project scheduling problem which was described above , the resources pricing is fixed based on the price policy and for this reason, issue of cost was not considered in these models.In this section, stepwise discount policy is defined based on different purchased resources (all unit discounts) in order to study effect of different pricing policies on the projects scheduling. (  ) shows unit price of the k type resource.

𝑆𝑆(𝑅𝑅
Based on this pricing policy, total purchasing cost of  unit from  ℎ resource i.e. (  ) will be as follows by calculating cumulative rate of Eq. (3).
Using Eq. ( 4), the project scheduling problem model is suggested as follows when the discount is suggested as pricing policy of the suppliers.
In constraint (5.2),  is total accessible budget for supply of the required resources of the project and function   is calculated though Eq. ( 4).

Parallel simulated annealing algorithm
In this section, we present a genetic algorithm based on combination of the genetic algorithms and simulated annealing (SA) to solve the resource-constrained project scheduling model and discount type pricing policy.This method uses genetic algorithm as the main framework of the proposed method and simulated annealing as the local search.In the proposed method, after the answer space was searched by the genetic algorithm using the primary population, local search algorithm improved the answer and the resulting answer returns to the genetic algorithm to start the next iteration.
Genetic algorithm is involved in local optimality like other algorithms based on random search [25].This specification is obtained due to myopic behavior of these algorithms but the variable neighborhood search is used as a metaheuristic algorithm which was considered recently to improve many methods based on local search [26].
This algorithm has been based on systematic and concurrent search of some types of neighborhood structures.Interesting specifications of this algorithm result in its increasing use in different problems.Flexibility and easy execution for different problems, use of many neighborhood search structures and low probability of its involvement in local optimality are among the interesting specifications of this metaheuristic algorithm.
In this research, we use simulated annealing algorithm as a new operator in genetic algorithm in order to improve quality of the answers resulting from genetic algorithm and reduce the local optimality.As mentioned above, simulated annealing algorithm has removed this weakness of metaheuristic algorithms due to use of several structures of neighborhood search.

Main framework of the proposed algorithm
Genetic algorithm which was introduced and expanded by Holland [26] is one of the random search methods which are used to solve the complex optimization problems.This algorithm starts with some initial answers which have been randomly produced and produces the next generations by applying the probable operators.Genetic algorithm directs the population to the optimal answer by the crossover and mutation operators.Although it is very difficult to design a genetic algorithm which directs the answer to optimality, in a genetic algorithm, a very important issue is evaluation of the answers quality by the fitness function.Answers with better values of fitness function have more potential to produce the next generation.Genetic process continues for many times until a predetermined condition is realized.
In recent years, many metaheuristic algorithms have been used to solve the project scheduling problem.In this research, we use genetic algorithm as main framework of the proposed solution.After production of the primary population and selection of parents, crossover and mutation operators are applied randomly.Answer obtained from this reproduction is exposed to the third operator called hybrid operator.This new operator receives the reproduced child as the input and does local search based on the simulated annealing method to improve the answer.This improved answer is transferred to the next generation.This procedure prevents from involvement in local optimality while improving quality of the resulting answer and also reduces convergence period of the algorithm.

Answer representation 𝑀𝑀
Since constraint of the accessible resources relates to the entire programming horizon, a bimodal allocation vector ( 1 ,  2 ) corresponding to all activities of the project is suitable for representing answers of the resource -constrained project scheduling problem (chromosome) and discount policy in pricing.For example, for a network with four activities, answer of V={ 1 ,  2 } indicates allocation of mode  1 to activities 1, 2 and 4 as well as allocation of mode  2 to activity 3.In order to simplify answers representation, we represent  1 with code 0 and  2 with code 1 (Fig. 1).An important point in this representation is feasibility of the answers which is analyzed in the next sections.

Producing the primary population
The primary answers play very important role in access to the ultimate answers with better quality.In recent years, many researchers have discussed how quality of the primary answers is effective on quality of the ultimate answers.In this regard, we try to improve quality of the input answers to the algorithm by presenting a heuristic algorithm.
Step 1: we define list of the programmed activities as  and put  = ф Step 2: calculate reduction in completion time of the activities for each resource increase unit in case of change from  1 to  2 for all ′ member works: Step 3: put ′ member works in ascending order in terms of   .
Step 4: for the first activity of step 3, put   = 2 , exclude the activity from list ′ and add it to list L.
Based on the presented heuristic method, the available resources are allocated only to the activities which have the highest reduction in completion time for each additional resource unit.The answer obtained from this algorithm is included as one of the primary population answers and other population will replace by random production of the answer.

Evaluating answers
In this section, the produced answers are evaluated by the algorithm target function (completion time of project).As mentioned before, important point in evaluation of the answers is their feasibility.For this purpose, different feasibility tests have been used as fine function in literature.For one given answer , value of () shows feasibility of .
Vector  is corresponding to a feasible answer if and only if () = 0.The following fitness function is suggested using this fine function.
where  an upper is bound for the project completion time and is calculated by adding the longest possible time for performing the activities.

Crossover and mutation operators
Crossover reproduces children as an operator by changing some parts of the parent's chromosomes.Children have some traits of their parents and for this reason; they can be better or worse than their parents.In this research, uniform crossover where a probability vector is used to reproduce children is applied.Size of the vector equals the number of project activities.If the probability is smaller than 0,5, its corresponding element will be transferred from the first parent to the child, otherwise, the second parent is selected.In addition, in order to increase variety of the answers, mutation operator is applied to make random changes in chromosomes.In this research, random mutation is used in which a chromosome element is randomly selected and its state changes ( 2 →  1 ,  1 →  2 ) .

Hybrid operator: simulated annealing
As mentioned before, a new operator is added to the genetic algorithm based on simulated annealing as the hybrid operator in order to improve answers resulting from genetic algorithm and reduce local optimality probability.Simulated annealing algorithm starts with an initial feasible answer for the problem and considers it as current solution.At any temperature, we study neighborhoods of the current solution with iterations.After production of neighbor answer and calculation of its target function, its target function is compared with current target function.Weak answers are accepted with regard to Boltzmann function in which ΔE is variation in target function for going from one answer to neighbor answer and   is current temperature.Therefore, acceptance of a worse movement is a function of two M M M M 0 0 1 0 temperatures of system and variations of target function.The lower the system temperature, the less probable the weak movement is accepted.After completion of local access, the current solution is compared with the best answer of algorithm and it is replaced if required.
According to temperature reduction model, temperature is reduced and the next stage of algorithm starts.Strong point of simulated annealing algorithm is escape from optimal local point [28].On the contrary, one of the weaknesses is high volume of computer calculation to reach very good answers.The mentioned problems occur when the mentioned algorithm is not populated oriented and starts searching with one point of answer space.One of the solutions to this problem is that some parallel simulated annealing processes are used for searching answer space.There are different ways for parallelizing simulated annealing algorithms.Hirayasu et al suggested one of the most efficient parallelizing methods [29].In their methods, searching process starts with some points of answer space and crossover operator is used in genetic algorithm for production of good answers after some steps.In this article, a similar method is used for solution.
In this method,  simulated annealing process ) ..., , , ( starts searching answer space from stochastic answers.These operations start with common initial temperature for all processes.In each process, search operations are performed until reaching balance according to simulated annealing algorithm.These operations are performed in each one of the parallel processes and a neighbor answer is produced from the current answer and target function value is calculated.If neighbor answer is better than the current answer, this answer replaces current answer; otherwise, neighbor answer is accepted with probability of p according to .In this relation, Δf is a target function variation and  is current temperature of the process.Neighborhood production process continues until reaching balance.After all  processes reached balance, n child is produced from n current answer.Then n better answer is considered as current answers out of total n current answer and n child answer.Temperature reduction process is performed and in case the final condition is not met, simulated annealing processes start their activity with the current answers.

Temperature reduction mechanism
In this research, temperature reduction model has been considered as linear.In linear model, temperature decreases on the basis of stage number.On the other hand, temperature reduction rate is adjusted during execution of algorithm.Temperature reduction rule and movement toward system cooling are as cooling function.
) where α indicates temperature reduction factor and a constant lower than 1 and parameter t controls frequencies of system reduction.

Algorithm stoppage criterion
The algorithm has been designed in such a manner that when algorithm reaches final temperature, algorithm will be ended and the last target optimizing answers will emerge as answer in output.Fig. 2 shows the used parallel simulated annealing algorithm diagram in this research.

Experimental results
In this section, performance of the proposed hybrid approach is evaluated.For this purpose, standard examples are required.Since the proposed model is presented for the first time in the literature, a random production approach is used to give examples in small and large dimensions.

Data production for the standard problems
In order to evaluate the proposed approach, different classes of problems are produced in this section for the resource -constrained project scheduling problem in discount pricing policy which is presented by the supplier.Data required for examples of this problem are classified as the following factors: Number of Activity (n): two options are considered for design of the number of activities  = {40, 80}.
Precedence Relations: in order to design a network with n activity, first, a random number which indicates the number of precedence relations is produced.Then two activities are randomly selected for each relation.The first number is the precedence activity number and the second number is the successor activity number.
Upper and Lower Bound for Activities (  1 and   2 ): in order to produce resources constraints for activities, a uniform distribution is used in interval (100, 200) for lower bound and uniform distribution is used in interval (300, 400) for the upper bound.
Upper and Lower Bound for Performing Activities (  1 and   2 ): it is assumed that the lower bound for performing activities is in interval (15,30) and its upper bound is in interval (30, 45) and as uniform distribution.
Discount function: the number of points of failure as well as unit price of the resources in points of failure is randomly selected.
Accessible budget (M): for the available budget, two options are assumed:  = {5000, 7500}.In order to produce the standard examples, resourceconstrained scheduling problem and resources price discount, we should only execute steps 1 to 6 based on the examples.

Regulating parameters of the proposed algorithm
In recent years, it was mentioned in the literature that efficiency of the solution algorithms based on random search is considerably affected by different factors.One of these factors is parameters of these algorithms [24].A genetic algorithm cannot reach desirable quality of the ultimate answer if its parameters have not been set correctly.For this reason, we try to obtain the most suitable values from parameters of the proposed algorithm in order to have access to the desirable answers.In this regard, we suggest stepwise regression approach for the first time in literature to set a metaheuristic algorithm.When we face many independent variables, it is very important to determine suitable combination of these parameters to be present in the model and estimate dependent variable.This will be more important when the used algorithm has been designed as a combination of the metaheuristic algorithms and consequently there are more parameters of hybrid method and it will be more important to set parameters.Stepwise regression is a stable statistical tool to select the best combination of model variables which can be done with forward selection and backward elimination.Forward selection starts with an empty set of variables and adds the variable to the model which is most effective on the model step by step while backward elimination starts with a set of total variables.Stepwise regression approach is a combination of two approaches of forward selection and backward elimination.In fact, this approach adds variables with the highest effect to the model and eliminates the variables with the lowest effect from the model.Our proposed algorithm includes parameters of population size, generation number, crossover rate, mutation rate and the number of simulated annealing iterations.The proposed approach for parameters setting of the algorithm includes the following two stages: 1-Determining parameters effective on performance of the algorithm and regression relation on these parameters 2-Optimizing the obtained regression relation in the first stage to determine their suitable values.
Tab. 1 ranges the algorithm parameters for the large and small dimensions problems.
In order to facilitate calculation of stepwise regression, we used Stepwise graphical user interface of MATLAB software.Results of these calculations are shown in Figs. 3 and 4 for large and small examples.For each important parameter in the model, regression coefficients which have been calculated through least squares and their confidence level of 90 % are marked with blue color.For unimportant parameters which have not been included in the model, a red curve indicates their regression coefficients in case they enter the model and the confidence level is 90 %.In these diagrams,  1 indicates population size,  2 indicates generation number,  3 indicates crossover rate,  4 indicates mutation rate and  5 indicates number of neighborhood search iterations.As it is observed based on the obtained results, parameters of generation number, crossover rate and mutation rate are effective for small dimensions of the problem as well as parameters of population size and generation number for the large dimensions of the problem (marked with blue color).In order to reconfirm effectiveness of these parameters, added variable plots are drawn in Figs. 5 and 6 for each one of the model parameters in both dimensions of the problem.
These plots are used to recognize important variables in multivariate regression functions and their non-uniform behaviors indicate effectiveness of the variable on regression model.As observed above, results obtained from the drawn plots confirm results of the stepwise regression.Behavior of unimportant parameters in these plots is uniform and behavior of effective parameters is linear with nonzero slope which confirms the previous results.In order to determine suitable values of parameters, two optimization problems with the above target function and with constraints of table 1 are executed for problems with large and small dimensions.

Comparing results of stepwise regression and Taguchi method
In this section, we evaluate results obtained from stepwise regression approach based on proposed Taguchi method.This method which was presented for the first time by Mr. Taguchi in 1960 is one of the common methods for parameter setting of the metaheuristic algorithms which has been used in the literature in recent years.Taguchi has presented an evaluation as signal to noise ratios for his proposed method.
In this method, target functions are classified into three main groups of the smaller the Better, The Larger the Better and The Nominal value -is Expected each using special formula to calculate / ratios.Since goal of this research is to minimize completion time of the total project, therefore, we use the smaller the Better function in this research.For this purpose, the following formula is used to evaluate target function in which the larger the function, the better the quality of answers and the shorter the completion time of the project. ) Based on levels of the parameters which were considered in the previous section to set parameters of the proposed hybrid algorithm, suitable values for parameters of the algorithm were calculated which are summarized in Tab. 2.

Results of simulation
In order to evaluate the answers obtained from the proposed hybrid algorithm, answers resulting from LINGO 8.0 Optimization Software are used.In fact, LINGO is able to solve nonlinear optimization problems including integer and continuous variables using branch and bound algorithm.In order to produce standard problems, 10 examples with small dimensions (n = 40) and 10 examples with large dimensions (n = 80) are produced using the described procedure in section 4.1.For each example, LINGO model was executed and results were recorded.
On the other hand, due to random nature of the proposed genetic approach in section 3, the obtained answers are not equal in several iterations for a special example.For this reason, for each example, the proposed algorithm was executed for 30 times and the obtained results were recorded.In order to improve performance of the algorithm, population size for all problems equals to 100, generation number equals to 250, crossover and hybrid operator equals to 0,98 and mutation operator probability equals to 0,05.Results obtained from calculations are given in Tab. 3.
As results of Tab. 3 show, the proposed method in 13 problems is able to find the optimal answer.Difference in values of target function for the proposed algorithms was only 0,4358 % compared with branch and bound approach which shows good reliability and effectiveness of the proposed method.However, the Taguchi method reaches optimal answer for 4 examples and error of this method is averagely 0,5972 compared with the optimal answer.The proposed approach has produced better answer in 12 examples of Taguchi method while Taguchi approach has produced better answer than stepwise regression method only in 4 examples.

Conclusion
In this research, the resource-constrained project scheduling problem has been combined with discounted cash flow policy as the pricing strategy of the project suppliers for the first time in the literature.It has been assumed that the resource -constrained project scheduling model with supply chain of the project has been evaluated and studied to enable the managers and programmers to make the best decision about optimal use of budget of the project in case they face suggestion of the suppliers to change price of the resources step by step.
After mathematical formulation of the problem, a new genetic method has been suggested based on combination of genetic algorithm and simulated annealing algorithm.This method uses genetic algorithm as the main body of method and simulated annealing algorithm as a new operator.This new operator is beside crossover and mutation operators and tries to apply all possible improvements by applying some neighborhood search structures through the local search.
Then the improved answer is returned into the main algorithm to form the next generation.Use of this operator prevents premature convergence of the algorithm and involvement in local optimality while improving quality of the answers.Application of stepwise regression approach was suggested for the first time in the literature to set parameters of the metaheuristic algorithms in this paper.At the end, some standard problems in large and small dimensions are produced and results of a comparative analysis shows suitable performance of the proposed method well in order to evaluate efficiency of the proposed approach.The results showed that our proposed algorithm is too fast and it obtained optimal answers or near optimal answers in a few seconds.

Figure 1
Figure 1 Answer representation

Figure 2
Figure 2 Parallel simulated annealing algorithm diagram

Figure 3 Figure 4 Figure 5
Figure 3 Results of stepwise regression calculation for the small dimensions of examples

Figure 6
Figure 6 Added variable plot for the large dimensions of problems After this stage, regression function of the project completion time is extracted from examples for the small and large dimensions based on regression coefficients of the effective parameters. +1 = 203,642 − 0,9234  2 − 1,4480 3 − 1,5570 4  +1 = 52,577 − 1,4683  1 − 0,6623

Table 1
Algorithms parameters bounds

Table 2
Results of Taguchi method for parameters of the algorithm