1. INTRODUCTION
Project management techniques have evolved and been practiced for half a century. A large part of efforts in the area of development of project management concepts has been dedicated to generating scheduling models (Sabela, 2020;Hartani et al., 2021;Barzamini et al., 2022;Zabihi et al., 2012). As such, project scheduling is recognized as one of the fundamental pillars of project management (Tirkolaee et al., 2019). Among the existing scheduling problems, the resourceconstrained project scheduling problem (RCPSP) has been introduced by various researchers as one of the most significant and challenging operational problems (Koulinas et al., 2014). In fact, it is the most general scheduling problem, and job shop scheduling, flow shop scheduling, and other scheduling problems are all a subset of this issue (Kurilova, 2021). The main objective of RCPSP is to allocate a set of constrained resources to activities by adhering to prerequisite relationships to achieve various and specific goals. Nevertheless, the most common objective of these methods is to minimize the project completion time (Pahlevan et al., 2021;Golmohammadi et al., 2016;Alwreikat and Rjoub, 2020;Khan et al., 2020). Various objective functions considered for RCPSPs divide the problems into different branches, one of which is the classic RCPSP. The primary goal of this technique is to reduce the implementation time of a project by finding a proper sequence of project activities in a way that precedence and latency constraints of the project and the different types of resource constraints available are satisfied simultaneously (Susilawati, 2021;Goli and Malmir, 2020;Shannaq et al., 2020;Seyedhosseini et al., 2016). This could be recognized as an NPhard optimization problem (Bukata et al., 2015). Proposing more efficient solutions for such problems has been taken into account by several researchers in the past few years (Valls et al., 2005;Rjoub et al., 2017).
Different methods have been proposed in the literature for solving RCPSP, including exact solution generation methods and algorithms. One of these algorithms is the branch and bound, which has been presented as a proper technique for solving RCPSP and achieving an optimal solution (Aravindhan et al., 2021). Given the classification of RCPSP as an NPhard problem in terms of computational complexity, these methods can only solve examples with less than 60 activities (Koulinas et al., 2014). Therefore, heuristic and metaheuristic algorithms are the only solutions for solving RCPSP with a high number of activities (Valls et al., 2005). Several algorithms have been proposed for these problems; for instance, Kharat and Shelar (2021), presented a solution based on priority rules. In addition, Nonobe and Ibaraki (2002) and Klein (2000) developed techniques based on tabu search algorithms. In another study, Bouleimen and Lecocq (2003) applied a simulated annealing algorithm, and Susminingsih et al. (2021) and Fleszar and Hindi (2004) created neighborhood search techniques. Among the populationbased algorithms, the genetics algorithm was applied in (Qazani et al., 2021;Valls et al., 2008), whereas scatter search algorithm and ant colony algorithm were used by Debels et al. (2006) and Harmoko (2021). Other algorithms used to solve RCPSP include artificial neural network algorithm (Sembekov et al., 2021), bee colony algorithm (Ziarati et al., 2011), and frog jumping steps algorithm (Fang and Wang, 2012). In addition, many researchers have used a hybrid approach to solving RCPSP. Wang and Fang (2012), and Susminingsih et al. (2011), developed techniques with a hybrid approach.
The present study presents a novel technique to solve classic RCPSP using a forest optimization algorithm (FOA) and gradual simulated annealing. The remainder of the article is structured as follows: Section 2 describes the mathematical model of the problem and section 3 introduces the FOA. In addition, section 4 briefly explains the gradual simulated annealing, and section 5 presents the proposed method. Furthermore, section 6 evaluates the performance of the proposed algorithm. Finally, section 7 concludes and presents suggestions for future studies.
2. STATEMENT OF THE PROBLEM
If we consider the activities of a project in a set such as N and the prerequisite relationships between the activities in a set such as A, RCPSP can be represented by a directed acyclic graph (G = (N, A). Activities are numbered from 0 to 1 + n, in which activities 0 and 1 + n are dummy activities and indicate the beginning and end of the project, respectively. The implementation time of each activity, such as I, is shown by di, and d0=dn+1=0. In addition, there are K types of renewable resources, and there are RK units of the kth resource. Each activity, such as i, requires a certain amount of k resource to be implemented, which is shown by rik and r0k = rn+1 = 0. The start time and the end time of each activity, such as i, is shown by Si and Fi, respectively. In addition, P(t) presents activities running at t time, and the T variable is equal to the longest project implementation time.
Before presenting the model, the variables used in it are described:

A_{i} : A set of all risk responses;

R : A set of all risks;

W : A set of project activities;

M : A pair of responses eliminating each other;

M : A set of paired responses with positive effects on each other;

B : Total budget to provide project risk responses;

W_{k} : kth activity;

R_{j} : Response to jth risk;

A_{i} : Response to ith risk;

C_{ij} : Cost of implementing ith response to jth risk

${S}_{j}^{k}$ : Estimating days of delay in case of jth risk during kth activity;

${q}_{j}^{k}$ : Estimating quality loss in case of jth risk during kth activity;

e_{ij} : Estimating risk reduction rate in terms of ith response to jth risk;

${S}_{ij}^{k}$ : Estimating decrease in working days of kth activity in terms of ith response to jth risk;

${q}_{ij}^{k}$ : Estimating improved quality of kth activity in terms of ith response to jth risk;

ε^{k} : Time between end of kth activity and its subsequent activity;

δ^{k} : Acceptable quality of kth activity without affecting that of its subsequent ones;

T_{Max} : Maximum allowable project delay;

Q_{Max} : Maximum allowable project quality loss;

X_{ij} : Zeroone decision variable: whenever ith response is chosen for jth risk, it is equal to one, and otherwise zero.
The first objective function of the zeroone goalprogramming model maximizes the effect of risk response (Objective Function 1); in other words, this objective function takes an optimal response to reduce the negative impact of risks as much as possible. However, the main point not mentioned in previous research is that the project manager does not aim to mitigate negative effects irrespective of response costs. This means that the minimum cost to deal with risks must be considered (Objective Function 2) along with lowering negative effects, so this generates a dualobjective optimization problem. The model is illustrated below.
In this model, constraint 3 ensures that the total cost of providing risk responses does not exceed the total budget. Inserting Max_{j} in this constraint is also to avoid recounting the responses. As well, constraint 4 shows that each activity must be fulfilled within its allowable time and not cause delays in the subsequent activities. Constraint 5 is the same as before, except that it confirms that reducing the quality of one activity does not lower that of subsequent ones and result in an unacceptable decline in project quality. Moreover, constraint 6 indicates that the final activity ends with an allowable delay and the project is not delayed in an unreasonable manner. Constraint 7 also determines the standard of the final project quality and does not reduce it. As well, Constraint 8 is utilized for mutually exclusive responses. That is, if one of them is employed, the other response should not be chosen. Constraint 9 is further applied for responses, one of which must be selected. Furthermore, Constraint 10 maintains that if one of these two responses is chosen, the other one must be subsequently selected. The last constraint also shows the problem decision variable.
3. FOREST OPTIMIZATION ALGORITHM (FOA)
FOA is deemed one of the newest evolutionary algorithms introduced by Ghaemi and Derakhshani (2014). The algorithm has been inspired by various seed dispersal methods of trees in forests, such as local seed dispersal and longdistance seed dispersal (Cain et al., 2000).
3.1 Construction of Primary Forest Trees
Similar to other evolutionary algorithms, FOA starts with an initial population of trees, where each tree can be a solution for the problem. In addition to the values of the variables, a tree also has an age, which is initially considered zero. In addition, each tree is shown in the form of an array of 1 × (1 + N_{var}), where N_{var} shows the number of variables in the problem.
3.2 Local Seeding
The operator acts by choosing a cell from the tree for each zeroaged tree and selecting a random number, such as r, from the range of (−Δx and Δx). Afterwards, the value of the house is added to the r number. Δx is a value smaller than the upper bound of the variables. The process is repeated at a frequency of local seeding changes (LSC). In fact, LSC shows the number of trees that must be generated from each zeroaged tree at the local seeding stage. After performing the above operation, the age variable of newly formed trees is set to zero, and the variable increases by one unit for other trees in the forest.
3.3 Forest Size Restriction
Two variables are used to restrict the trees in the forest. The first variable is life time, which determines the maximum age of each tree while the second variable is area limit, which is equal to the maximum number of trees that can exist in the forest. After local dispersal, the trees that have reached their life time are eliminated from the forest and added to the set of candidates, which encompasses the trees eliminated from the forest. Afterwards, if the number of trees left in the forest is more than the area limit variable, the trees in the forest are sorted based on the value of their fit function and some trees are selected from the best trees to the number of area limit for survival of the forest and the rest of the trees are added to the candidate set.
3.4 Global Seeding
In this step, a number of trees in the set of candidates are first selected, and the number is determined by the transfer rate variable, which represents the percentage of trees in the set of candidates. Afterwards, a certain number of cells are selected from each array (tree), and the value of the cell in the tree is replaced by a new value of the permissible interval for problem variables. The number of cells selected from each tree is determined by another variable of the algorithm known as global seeding changes (GSC). After the mentioned operations, the variable of tree age is set at zero and added to the forest.
3.5 Maintenance and Updating the Best Tree
After each iteration of the algorithm and performing the above steps, the best tree is selected in terms of the value of the fit function and its age value is set to zero so that it can perform local dispersal operations in the next step.
3.6 Termination Conditions
Similar to other metaheuristic algorithms, three termination conditions can be considered for forest optimization algorithm: first, a certain number of algorithm iteration; second, no change is observed in the value of the fit function of the best tree during several iterations, and third, the condition for achieving a certain level of accuracy.
4. GRADUAL SIMULATED ANNEALING
Simulated annealing was developed by Christofides et al. (1987) as an alternative to local search and a possible general method for solving optimization problems. This technique is practiced to reach a properly sorted solid matter and minimize its energy. This technique involves exposing the substance to high temperature and then gradually reducing the temperature. The main idea in this method is based on accepting the worse solutions in some conditions in order to escape the local optimum. Simulated annealing is so named because of its analogy to the process of physical. The gradual simulated annealing process, which reduces energy in solid, was defined by Agarwal et al. (2011) as follows: at each stage, the atom is slightly displaced, which leads to a change in the energy system shown by ΔE. If ΔE≤0, the displacement of two atoms will be accepted and the displaced solid structure or atom will be used as the starting point for the next step. On the other hand, if ΔE>0 , the process will be dealt with in a probabilistic manner, meaning that the possibility of acceptance of the solid structure will be determined using the equation below, where T is the temperature degree and k_{b} is the Boltzmann constant. To solve an optimization problem, the simulated annealing algorithm starts with an initial solution and moves to the neighboring solutions in a loop. The algorithm will select the neighboring solution (i.e., moves the solution) if it is more efficient than the current solution. Otherwise, the algorithm accepts the solution as the current solution with the possibility of exp(−ΔE/T) . In this equation, ΔE is the difference between the objective function of the current solution and the neighboring solution, and T is a parameter known as temperature. The temperature is gradually reduced after performing several iterations at each temperature. The temperature is set very high at the initial stages so that there is a higher probability of accepting worse solutions. There will be a lower chance of accepting worse solutions in the final stages by a gradual decrease of temperature, and the algorithm will converge toward a proper solution.
6. PROPOSED METHOD TO SOLVE RCPSP USING FOREST OPTIMIZATION ALGORITHM
Since forest optimization algorithm has been introduced to solve continuous problems and RCPSP is a discrete problem, some changes should be made, such as redefining LSC and GLC steps in the initial version of the algorithm to adapt it to the problem conditions in order to solve the RCPSP using the forest optimization algorithm. In addition, the gradual simulated annealing mechanism is used in the present article to improve the performance of the proposed algorithm and increase the possibility of transferring fit trees to the next stage at each step of the algorithm.
6.1 Structure of Trees
Each tree is shown with an array of 1 × (1 + N_{var}), where N_{var} is the number of project activities plus two dummy activities at the beginning and end of the process, and the other cell of the array is related to the variable of tree age. Each tree encompasses a series of project activities in a way that precedence and latency constraints of the project network are observed. In other words, each tree involves a doable sequence for performing project activities. To generate each tree, one activity is randomly selected from activities, the prerequisites for which are placed in the current home, for each cell of the array, and the process continues until reaching the final activity.
6.2 Implementation Duration Calculation Method
In order to estimate the duration of the project for each tree, we start from the array and identify activities that could be performed simultaneously (i.e., their prerequisites have been performed completely). These activities are implemented simultaneously if allowed by resource constraints. Otherwise, the activity with the shortest implementation time remained is selected among the running activities and ended, followed by reassessing the mentioned operation at the current time. The process continues until reaching the n+1 activity; evidently, the fit value of each tree equates to the start time of the n+1 process.
6.3 Simulation of Local Seed Dispersal Operator
In order to simulate the local seed dispersal operator, two elements of the array are first selected and the current activity in the cell with a larger number is transferred to a cell with a smaller number if the precedence and latency constraints of the project network are not violated. Afterwards, the activities between the two selected cells are transferred one unit forward. Otherwise, the process is repeated until the process is made possible.
6.4 Simulation of Global Seed Dispersal Operator
In order to simulate the global seed dispersal operator, a cell of the tree array is selected first, and the process of generating a doable solution is continued from that cell. Since one cell of the array is randomly selected from the candidate activities, there is a high possibility of generating new sequences by this technique. In the proposed method, the GSC variable equates to the maximum number of selected cells. If the value is chosen to be high, there is a lower possibility of changes in the structure of the selected tree. On the other hand, if the value is chosen to be very low, the selected tree will be extremely changed after the global seed dispersal process. Therefore, it is crucial to select an appropriate value for the parameter in the algorithm implementation process. Following the generation of a new tree with the method explained above, the gradual simulated annealing mechanism is used to select between the tree chosen from the candidate list and the changed tree to be added to the forest in the next algorithm stages. To this end, we first calculate the time (fit) of the tree changed in the previous stage. Afterwards, the selected tree will be added to the forest and its age will be set at zero if its time is shorter than the time of the tree selected from the list of candidates. Otherwise, an appropriate solution will be selected with the possibility of ${e}^{\frac{df}{{T}_{0}}}$ and if the probability is larger than a uniform random number between zero and one.
7. COMPUTATIONAL RESULTS
The examples existing in PSPLIB, which are known as reference examples fir RSPSP, are used to evaluate the performance of the proposed method. This set includes 90, 60, 30 and 120activity examples, and the 30, 60, and 120activity examples have been used in previous studies to assess the performance of various algorithms. While an optimal solution exists for 30activity examples, the lower bound scheduling, which is obtained by the classic critical path method (CPM), is used as the criterion for assessing the deviation of algorithms. In this article, 20 examples are selected randomly from the examples existing for each mentioned group. All of these examples have been mentioned in (Cain et al., 2000). The proposed method is coded by C# programming language and run in Visual studio 2013 programming space on a computer with 2.5 GHz Intel Core i5 processor and 4 GB of RAM with Windows 8.1 operating system. Adjusting the variables of any evolutionary algorithm plays an important role in its implementation. The variables are initialized for the proposed method following the implementation of various experiments, as follows: the number of initial forest trees and remaining trees in the forest at each stage (area limit) equal to six, GSC variable value equal to onefifth of the number of activities per problem, LSC variable value equal to five, transfer rate variable value equal to four, maximum age of each tree (life time) equal to 5, and cooling rate equal to 0.85. In addition, the termination condition of the algorithm, is equal to producing a certain number of scheduling (trees), similar to previous studies. In this article, this number is equal to 5000.
Regarding 30activity examples, the proposed method is able to achieve optimal solutions for all 20 examples selected in a rapid manner and solve the examples without deviation. Table 1 shows the lower bound of scheduling for 60activity examples and the solutions obtained in this regard. As observed, the proposed method is able to solve the examples of this category with very little difference compared to the low bound scheduling.
Table 2 presents the lower bound of scheduling for examples with 120 activities and solutions obtained by the proposed method. The solutions of the presented method show the suitable applicability of the method based on the high number of activities and the mean deviation of the existing methods, which is equal to 34.384. Figure 1 compares the percentage of deviation in two 60 and 120activity examples.
According to Figure 1, the deviation percentage of different activities is higher when a project with 120 activities is assessed, compared to a project with 60 activities. In other words, increasing the number of activities has a direct impact on the degree of project deviation, which confirms the fact that the larger the size of the project, the higher the necessity of resource consumption optimization in this regard.
8. CONCLUSION
The present study suggested a novel technique for RCPSP based on a forest optimization algorithm. Since the primary forest optimization algorithm is an evolutionary algorithm for solving continuous optimization problems, a new method was proposed in this paper by applying changes and redefining its leading operators to solve the project scheduling problem with classic resource constraints. In addition, tree selection was optimized for the following stages of the algorithm by using the gradual simulated annealing mechanism in the global seed dispersal stage. According to computational results, the proposed method had an extremely suitable performance regarding solving the problem. Given the suitability of the proposed method in solving RCPSP, it should be used to solve other problems related to project scheduling, including multimode resourceconstrained project scheduling problems and resourceconstrained project scheduling problems with minimal and maximal time lags in future studies.