Combining Tabu Search and Genetic Algorithm

Kevin Boston and Pete Bettinger

ABSTRACT. The performance of two new heuristics that combine tabu search with a genetic crossover technique was examined for their usefulness in solving spatially constrained harvest scheduling problems. One heuristic utilized tabu search with 1-opt moves and a genetic crossover technique (TS/ GA), and the other utilized tabu search with both 1-opt and 2-opt moves and a genetic crossover technique (TS2/GA). The heuristics were tested on four problems. Three problems used the simple unit-restriction model (URM) to portray greenup constraints, allowing them to be solved with mathematical programming techniques and thus providing a benchmark to compare against the solutions produced by the heuristic techniques. These were considered hypothetical problems, since the datasets were grids, and the age classes were assigned with three different rules. The fourth problem uses an operational dataset and includes both a maximum opening size constraint, which limits the size of an individual opening, and a maximum average opening size constraint, which represents the greenup constraint contained in the American Forest and Paper Association (AF&PA) Sustainable Forestry Initiative (AF&PA 2000). The greenup constraints contained in these problems accurately portray the greenup constraints facing the forest industry in the United States. For the three hypothetical problems, the TS/GA technique found solutions with an objective function value between 96.6% and 99.1% of an estimated optimal value, and between 93.4% to 94.4% of a relaxed linear programming value. The TS2/GA found better solutions for all three hypothetical problems, with the objective function values between 98.2% and 99.7% of the estimated optimal value and the 96.3% to 97.2% of the relaxed linear programming value. Compared with a general tabu search technique that used 1-opt moves (TS), adding the genetic crossover technique resulted in a 2% increase in the objective function value, and adding the 2-opt intensification capability resulted in a further 1.5% improvement. A similar pattern was observed when the problem using the operational dataset was solved. The addition of the genetic crossover technique did not increase the time required to produce a solution to any particular problem, yet the addition of the 2-opt procedure did. TS2/GA required about twice as much time as did TS or TS/GA when applied to the three hypothetical datasets, and 67% more time when applied to the problem using the operational dataset. FOR. SCI. 48(1):35–46. Key Words: Harvest scheduling, tabu search, genetic algorithms, spatially constrained scheduling problems.

F

U NITED STATES and abroad are increasingly addressing spatial goals, through either regulatory or voluntary agreements. Today there is an increased awareness of the importance of OREST PLANNING EFFORTS IN THE

landscape pattern on wildlife populations, and spatially constrained harvest scheduling problems are now being solved that consider spatial habitat requirements. For example, some forest planning problems now include the

Kevin Boston, Carter Holt Harvey Fibre Solutions, Grayburn House, Leith Place, Tokoroa, New Zealand—Phone: 64-7-886-2760; E-mail: [email protected]. Pete Bettinger, Department of Forest Resources, Oregon State University, Corvallis, OR 97331—Phone: (541) 737-8549; E-mail: [email protected]. Manuscript received March 2, 2000. Accepted January 5, 2001.

Copyright © 2002 by the Society of American Foresters Forest Science 48(1) 2002

35

Downloaded from https://academic.oup.com/forestscience/article-abstract/48/1/35/4617418 by Humboldt State University Library user on 08 December 2018

Combining Tabu Search and Genetic Algorithm Heuristic Techniques to Solve Spatial Harvest Scheduling Problems

36

Forest Science 48(1) 2002

Random search techniques can be divided into three subcategories: (a) Monte Carlo integer programming, (b) simulated annealing, and (c) genetic algorithms. Monte Carlo integer programming is a sampling procedure that, for a single iteration of the procedure, randomly selects treatments for all units, typically subject to some user-defined constraints. The value of the objective function is then calculated, and if it is better than the best value previously found, the newly created solution replaces the best solution. This process is repeated for a user-specified number of iterations. Several studies have reported the results of using Monte Carlo integer programming (O’Hara et al. 1989, Nelson and Brodie 1990, Clements et al. 1990), generally noting that the results are within about 95% of optimal values. A single iteration of simulated annealing (SA) generally involves randomly selecting a unit and harvest timing to be added (or removed) from the solution. After selecting the unit and timing, the algorithm determines whether adding this would result in a feasible solution. If it would, and by doing so the objective function value would increase, the unit and its harvest timing are added to the solution. If it would lead to a lower objective function value, the following variable is calculated: Z = e (( current - best ) / temperature)

(1)

where current = the current objective function value (prior to adding the unit and harvest timing), best = best objective function value found during the search (so far), temperature = a pseudonym for length of time technique has run. If Z is greater than a randomly chosen number between 0 and 1, the unit and timing are added to the solution. The ability to avoid being trapped in a local optimum by using the random variable to accept an inferior solution is the main component to simulated annealing. There are other interpretations of SA where infeasible moves are allowed, but the objective function values are penalized to discourage their selection. The longer the SA technique runs, the cooler the “temperature” gets, and the less likely the technique will accept an inferior solution. Lockwood and Moore (1992), Dahlin and Sallnas (1993), and Van Duesen (1999) have used SA with some success in solving spatially constrained harvest scheduling problems. Threshold accepting and the great deluge algorithm are similar to simulated annealing, and have been shown to produce very good results as well in problems outside of forestry (Dueck and Scheuer 1990, Dueck 1993). The guiding premise of genetic algorithms (GA) is that simulating evolution via a computer algorithm can help solve complex problems (Glover et al. 1995). Using a biological analogy, a solution to a combinatorial problem is considered a chromosome in a biological organism. Each decision variable becomes a gene on a chromosome (i.e., a solution viewed sequentially by decision variables numbered 1 to n). The potential values for the decision variables are the alleles (Reeves 1993). New solutions are created through crossovers and random mutations. Mitchell (1995) describes the typical GA beginning with a randomly generated population of chromosomes. These represent


location of minimum cost corridors to connect similar wildlife habitat, while other problems include the arrangement of cover and forage to improve big game habitat. These types of problems quickly exceed the capability of traditional mathematical programming techniques, because often they are formulated as nonlinear 0-1 integer programming problems. Optimizing these problems involves tracking the interaction of management units, leading to a problem formulation that involves a large number of constraints and binary variables. Since solving these types of problems is formidable, no one particular solution approach has become generally acceptable (Murray and Church 1995). Spatially constrained forest plans can take on many forms, including those involving the timing and arrangement of harvests as well as the timing and arrangement of habitat patches. Here we are considering spatially constrained problems that involve the timing and arrangement of harvest activities. The development of constraints that prevent harvests of adjacent management units within a given window of time (greenup period) are commonly called “adjacency constraints” (McDill and Braze 2000). Murray and Church (1996) and McDill and Braze (2000) provide a thorough classification of adjacency constraints as well as an evaluation of the solution time and solution generation performance of a subset of the constraint groups. In the search to find an efficient method for solving spatially constrained harvest scheduling problems, many different approaches have been developed. Weintraub and Navon (1976) were among the first to solve the combined harvest scheduling and transportation planning problem as a 0-1 mixed-integer programming problem. Kirby et al. (1980) solved the combined harvest scheduling and transportation planning problem simultaneously by using a 0-1 mixedinteger programming technique. Hof and Joyce (1992), Hof et al. (1994), and Bevers and Hof (2000) all extend the use of integer programming to a variety of spatially constrained scheduling problems, yet the ability to apply these techniques over large land areas is uncertain. In spatially constrained problems with adjacency constraints [primarily using the unit restriction model described in Murray (1999)], Murray and Church (1996) and Snyder and ReVelle (1997) demonstrated ways to allow mixed-integer solvers to more efficiently find the optimal solution to large problems, and showed that improved formulations of the adjacency constraints can increase the size of a spatially constrained problem that can be optimally solved. Finally, Borges et al. (1999) and Hoganson et al. (1998) developed a dynamic programming technique to solve spatially constrained problems, by allowing the algorithm to move across the landscape in “strips,” and obtaining an optimal result for each strip while also acknowledging the results from analysis of previous strips as hard constraints. Decomposition techniques operate similarly, where the planning problem is divided into smaller problems that can be solved independently (Weintraub et al. 1994). In addition to these approaches, random search and hill-climbing heuristics have shown promise for solving spatially constrained harvest scheduling problems.

Murray and Church (1995) compared results generated by SA, MCIP, a random start hill-climbing or interchange algorithm, and TS, and showed that the interchange, SA, and TS techniques produced high quality solutions. They also found that there were no significant differences in the resulting final solutions of each technique, given the quality of the initial randomly generated solution. While no one approach could consistently identify superior solutions based on any one random initial solution, TS did produce the best solutions to the assumed problems. Boston and Bettinger (1999) compared MCIP, SA, and TS on four problems with between 3,000 and 5,000 0-1 integer variables. Simulated annealing produced the best result for three of the four problems, with solutions between 93% and 98% of the optimal solution. Tabu search found the best solution to one of the four problems. This brief review has shown that many different techniques have been developed to solve spatially constrained harvest scheduling problems, and the challenges to their implementation as well as the quality of their results vary considerably. The review also reveals that there have been few attempts to combine these techniques into a metaheuristic procedure. Glover et al. (1995) describe the potential for developing hybrid TS and GA procedures for solving combinatorial optimization problems. The hybrid algorithms use normal TS procedures combined with GAs for intensification and diversification of the search process. They believe that these hybrid algorithms may lead to better solutions to large combinatorial problems than solutions generated by the individual techniques alone. Colorni et al. (1998) and Costa (1995) provide examples of using a hybrid TS/GA technique to solve large, highly constrained scheduling problems outside of forestry. This article describes the development and testing of two heuristic techniques: (1) a hybrid TS and GA technique that uses 1-opt decision choices, and (2) a similar technique with 2-opt decision choices [see Bettinger et al. (1999) for a discussion of 1-opt and 2-opt decision choices]. Hypothetical datasets were used to compare the results generated by these two heuristic techniques to (1) an estimate of the optimal solution, (2) the upper bound on the solution obtained by relaxing the problem and solving it as a linear program, and (3) results generated by a general TS technique which uses 1opt decision choices. Finally, a comparison among the heuristics is made using an operational-sized dataset that combines both a maximum opening size constraint and a maximum average opening size constraint.

Methods We next describe the processes that are used to develop the tactical forest plans whose quality is assessed. A two-phase approach was used, where in the first phase a 30-yr strategic plan (using linear programming) was developed. Each time period in the 30 yr strategic plan was 1 yr in length, and the spatial restrictions were relaxed. In the second phase, the tactical forest plans were developed. Here, the volumes determined by the strategic plan were used as goals, and the spatial restrictions were incorporated as constraints. Forest Science 48(1) 2002

37


potential solutions to the problem. The fitness (objective function value) is calculated for each chromosome. Two solutions are selected, with the probability of selection based on each solution’s fitness value: the higher the fitness value, the greater the probability of being selected. A random crossover point is chosen for the two chromosomes, and they are split there and recombined (e.g., solution 1 is split into parts 1a and 1b, solution 2 into parts 2a and 2b, then 1a2b is created, as well as 1b2a). The two new chromosomes are mutated by randomly changing the values for some of the genes until the new population has been created. This process is repeated for a user-specified number of iterations (Mitchell 1995). Mullen and Butler (1997) developed a GA to solve a harvest scheduling problem with adjacency (greenup) constraints. They were able to solve planning problems with nearly 1,000 management units and 10 to 15 time periods and obtain objective function values within 98.3 to 99.8 of a relaxed LP value. Hill-climbing techniques include gradient search and tabu search (TS). These methods generally select units for harvest in a deterministic manner, where the choice that yields the largest improvement in the objective function is added to the solution. Hof and Joyce (1992) describe a gradient search method that illustrates the potential for solving complicated nonlinear harvest scheduling problems that creates a desired arrangement of the landscapes. Murray and Church (1995) also describe a hill-climbing routine applied to two unit-restriction adjacency problems. Tabu search differs from gradient search techniques because it allows inferior or infeasible changes to be made to a solution, thus forcing the search process to explore other regions of the solution space and avoid being trapped in a local optimum. Tabu search, which will be described in more detail below, has been shown to find the best solutions to a number of large combinatorial problems where the optimal solution is considered impossible to obtain (Glover et al. 1995). Within forestry, TS techniques have been used to solve harvest scheduling problems that included spatial cover and forage habitat requirements for big game (Bettinger et al. 1997) and water temperature and sediment constraints (Bettinger et al. 1998). The quality of heuristic solutions has recently received a fair amount of attention (Murray and Church 1995, Boston and Bettinger 1999) since the optimal solution to a planning problem is not guaranteed. In general, four methods can be used to evaluate the quality of solutions generated by a heuristic technique. The best method is to compare solutions against a known optimal solution, followed by comparing them against a solution where the integer constraints have been relaxed, providing an upper bound to the integer solution. The third method is to compare solutions with those produced by another heuristic technique. This will, however, allow only a relative comparison of better or worse solutions. The final method is to compare the solutions with an estimated optimal solution found using extreme value statistics, but this can result in misleading conclusions (Boston and Bettinger 1999).

N

maximize

T

ÂÂ n =1

( Revnt - Lcnt )Vnt Ant

t =1

(2)

where: N = number of harvest units; n = harvest unit; T = number of time periods; t = time period; Revnt = revenue per m3 for unit n harvested in time period t; Lcnt = logging cost per m3 for unit n harvested in time period t; Vnt = volume per hectare in unit n harvested during time period t; Ant = continuous variable representing the area of unit n harvested during time period t. All harvest units were assumed to be regenerated and available for harvest in later periods. A minimum harvest age was assumed for each timber stand, and an even-flow constraint was included for the three hypothetical problems: N

Â n =1

N

Vnt Ant -

ÂV

nt +1 Ant +1

n =1

=0

"t; t = 1,..., T - 1 (3)

Formulation of the Strategic Forest Plan for the Operational Dataset with Maximum Opening Size and Maximum Average Opening Size Constraints The strategic forest plan for the operational dataset with maximum opening size and maximum average opening size constraints was formulated similar to that for the three hypothetical datasets with URM greenup requirements. The objective function was similar to Equation (2), yet now certain product volumes are explicitly recognized: maximize

N

T

J

n =1

t =1

j =1

Â Â Â ( Rev

nt

- Lcnt )Vjnt Ant

(4)

where J = number of forest products; j = forest product. All harvest units were regenerated and available to harvest in later periods. In addition, the minimum harvest age was 19 yr, yet an even-flow constraint represented by Equation (3) was not used. Rather, the operational problem had three volume constraints, where the variation in the production of products (pulpwood, chip-and-sawlogs, and sawlogs) was constrained from one period to the next. Thus, individual product volumes could not change by more than a certain amount per time period, with deviationj representing that amount (5% for pulpwood and chip-andsaw logs, and 10% for sawlogs). Therefore, equations representing the upper and lower bounds of product production were used: 38


N

ÂV

N

jnt

n =1

X nt - (1 + deviation j )

ÂV

N

Â n =1

n =1

jnt +1 X nt +1

£0 (5)

" j, t; t = 1, 2,..., T - 1

N

V jnt X nt - (1 - deviation j )

ÂV n =1

jnt +1 X nt +1

£0 (6)

" j, t; t = 1, 2,..., T - 1

The volumes for each product j that were scheduled from this effort were used as the target volumes in the associated tactical plan, which had a 10 yr planning horizon and included the spatial constraints. Formulation of Tactical Problem for the Three Hypothetical Datasets with the Unit-Restriction Model Used to Account for Greenup Constraints Tactical plans for the three hypothetical datasets with greenup constraints were developed with a mathematical programming technique (integer programming) and three heuristic techniques: tabu search with 1-opt moves (TS), tabu search with 1-opt moves and a genetic crossover technique (TS/GA), and tabu search with 1-opt and 2-opt moves and a genetic crossover technique (TS2/GA). The heuristics are described shortly. The goal of tactical planning was to maximize net present value less a penalty for not achieving the volume goals, while accounting for spatial restrictions. The objective function for the problem is: maximize

N

T

n =1

t =1

Â Â ( Rev

T

-

nt

- Lcnt )Vnt Xnt

T

Â (Vp du ) - Â (Vp dl ) t

t =1

t

t

(7)

t

t =1

where Xnt = a binary (0,1) variable indicating whether unit n is harvested during time period t; Vpt = volume penalty per m3 during time period t; dut = positive deviation from volume goal during time period t; and dlt = negative deviation from volume goal during time period t. The volume goals we used were those that resulted from the strategic planning effort. Depending on the goals considered, it may be difficult to formulate a heuristic technique in terms familiar to traditional mathematical programming. A good example is where a problem has target volume goals, as this one does. Penalties were used to force the heuristic search process toward solutions that would emulate the strategic planning solution. Three constraints sets were included in the tactical plans. The first was a volume constraint, which was an accounting constraint used to sum the volume harvested in each period and determine the deviations from the volume goal: N

Â (V

nt Xnt ) -

n =1

dunt + dlnt = volume goalt

"t

(8)

The model creates a decision variable in each period for those harvest units with an age 19 yr or greater based on the


Formulation of the Strategic Forest Plan for the Three Hypothetical Datasets with the Unit-Restriction Model Used to Account for Greenup Constraints The strategic plan was formulated as a Model I linear programming problem with the goal of maximizing the discounted net present value of activities over thirty 1 yr periods. The strategic plan provides the target harvest goals for the tactical plan. The objective function consists of maximizing the net present value of revenue from harvests less logging costs

initial age of the unit. The second constraint, a singularity constraint, limits each unit to no more than one treatment during the planning horizon. T

Â

Xnt £ 1

"n

(9)

t =1

Un

Tm

ÂÂ k =1

Formulation of Tactical Problem for the Operational Dataset with Maximum Opening Size and Maximum Average Opening Size Constraints The tactical forest plan with maximum opening size and maximum average opening size constraints includes three product goals (sawlog, pulpwood, chip-and-saw logs) in each period as well as two greenup constraints. The objective function is nearly identical to the three previous data sets in that it maximizes the net present value less the deviations from the three product goals established by the strategic plan. J

( Xkm ) + un Xnt £ un

"n, t

maximize

(10)

T

nt

- Lc nt ) Vjnt Xnt

j =1 n =1 t =1

J

where Un = the set of adjacent neighbors to unit n and un is the number of elements in set Un; k = adjacent neighbors to unit n; Tm = the set of near-time periods; m = a near-time period; and Xkm = a binary (0,1) variable indicating whether unit k is harvested during near-time period m. There was no attempt to develop an improved constraint formulation to decrease the solution time. The goal was to generate solutions to the URM problems that would allow for a comparison of the results of our algorithm with the mathematical programming results. Using improved adjacency constraint formulations, such as the grid-packing algorithm described by Snyder and ReVelle (1997) or the clique-based algorithm described Murray and Church (1996), could result in improved solution times. For example, the nondominated clique selection approach has been shown to provide an LP bound that is superior to the formulations used here, with less computational effort, mainly because the nondominated clique selection approach can identify larger cliques, it does not prevent the identification of the true minimal set, and it contains more partial clique overlap (Murray and Church 1996). While this type of clique set provides a tighter problem structure and results in a tighter LP bound, the IP results are exactly the same among the techniques examined. Thus our justification for using constraints that some may feel are inferior is that we were only interested in generating IP bounds to several problems to allow for a comparison of the results with the heuristic techniques described here. Although the constraints are simple in shape, they are larger than most used in the research, and with multiple period exclusion they are closer to the problem size encountered by industrial organizations. Thus we felt these

N

Â Â Â ( Rev

-

T

ÂÂ

J

(Vp jt du jt ) -

j =1 t =1

T

Â Â (Vp

jt

dl jt )

(11)

j =1 t =1

The formulation includes similar volume and singularity constraints as used in the hypothetical data sets: N

Â (V

jnt

Xnt ) - dunt + dlnt = volume goal j , t

n =1

" j, t

(12)

T

ÂX t =1

nt

£ 1 "n

(13)

The main difference here is the formulation of the spatial constraints. The first greenup constraint is the maximum opening size constraint. The maximum opening size is the most common representation of actual greenup constraints used in forestry. Oregon has a 48.5 ha maximum opening constraint. In California, the maximum allowable opening constraint is 16.2 ha. Sweden also has an opening limitation of 20 ha. The maximum opening size constraint restricts the size of an opening that would result from the harvest of one or more units during the greenup exclusion period. For a 3 yr greenup constraint used in this model, we define a set (Tm) of near-time periods (mz, where m1 = t – 3, m2 = t – 2, m3 = t – 1, m4 = t, m5 = t + 1, m6 = t + 2, and m7 = t + 3; all mz ≥ 0 otherwise not in Tm; all mz £ T otherwise not in Tm). Therefore, an opening is not just the harvests that occur in time period t, but also includes those around each unit n that have occurred during the near time periods. The maximum opening size constraint is: Forest Science 48(1) 2002

39


The third set of constraints are the spatial restrictions for greenup, or more commonly called, adjacency considerations. We used a unit restriction method (URM) of controlling adjacency, wherein when one unit is in a regenerated state, all other adjacent units are restricted from being clearcut (Murray 1999). Our URM constraint held for 2 yr and was similar to the OAM formulation describe by McDill and Braze (2000). It can be described, however, in this fashion: for a 2 yr greenup constraint, we defined a set (Tm) of neartime periods (mz, where m1 = t – 2, m2 = t – 1, m3 = t, m4 = t + 1, m5 = t + 2; all mz ≥ 0 otherwise not in Tm; all mz £ T otherwise not in Tm). We then sought to control adjacent openings around a management unit being harvested that not only could occur in time period t, but also could occur during the near-time periods. The adjacency constraint is:

problems would be suitable for comparing the results of the heuristic techniques with a mathematical approach. To generate the IP solution, a Pentium II with 233 Megahertz computer with 128 Megabytes of memory was used to run CPLEX version 5 (Cplex Optimization Inc. 1998), the solver, for over 30 hr, until the computer’s available memory was exhausted. The result was the estimated optimal value, which was compared to the heuristic’s objective function values. The relaxed LP value that was solved as the initial step of the mixed-integer program was also reported to provide an upper bound on the objective function value.

ù È Un Tm ( Ak X km )ú + ( An Xnt ) Í úû ÍÎ k =1 m =1 (14) £ maximum opening size "n where Xnt = 1; t

ÂÂ

È Un Tm ù ( Ak X km )ú + ( A f X ft ) = O ft Í ÍÎ k =1 m =1 úû

ÂÂ

(15)

where Un = the set of adjacent neighbors to focal unit f; k = adjacent neighbors to focal unit f; Tm = the set of near time periods; m = a near time period; Ak = area of unit k; Xkm = 0,1 variable indicating whether unit k is harvested during near time period m; Af = area of focal unit f; and Xft = 0,1 variable indicating whether focal unit f is harvested during time period t. Because we are attempting to calculate the average opening size, we do not wish to count openings more than once. This miscounting could occur if we allow each unit (n) in an opening (composed of multiple units n) to be considered the “center.” Therefore, only one unit n can be delineated as the focal center f of the opening in any time period, and the total number of openings equals the number of focal centers of openings. Thus, F

Â

N

f £

f =1

Â

F

n

and

n =1

Â

f =1

N

X ft £

ÂX

nt

n =1

(16)

where F = the total number of openings; f = the focal center of an opening, or an opening. The average opening size for a set of openings in a time period t can then be constrained with the following equation: È F Í O ft Í f =1 Î

Â

ù ú / F £ average opening size " t ú û

(17)

where Oft = an opening centered around a focal unit f during time period t. For this problem, the maximum average opening used was 48 ha (120 ac), the minimum requirement of SFI. The size of the problem, approximately 13,000 0-1 variables, and the number of permutations required to describe all possible opening and average opening combina40


Heuristic Techniques Three heuristic techniques were developed to illustrate the increasing effectiveness of these solution processes when intensification and diversification routines are incorporated into a solution procedure; one is a standard heuristic, and two can be considered metaheuristics. They are, as mentioned earlier: tabu search with 1-opt moves (TS), tabu search with 1-opt moves and a genetic crossover technique (TS/GA), and tabu search with 1-opt and 2-opt moves and a genetic crossover technique (TS2/GA). All of the techniques include a 1opt tabu search technique. The 1-opt tabu search technique considers the entire neighborhood (the set of potential new solutions that can be created from the current solution given a change in the current solution) at each iteration and allows just one harvest unit to change its timing; thus it is a 1-opt tabu search algorithm (Bettinger et al. 1999). It has also been described by Murray and Church (1995), Bettinger et al. (1997), and Boston and Bettinger (1999) in detail. Two parameters commonly used to control a tabu search process are aspiration criteria and tabu memory. For all three heuristic techniques, the aspiration criterion was assumed to be the best objective function value, and the tabu memory was assumed to be 100 iterations of the process. At each iteration of a 1-opt tabu search process, one of five of the following events will occur while examining in succession each potential move in the 1-opt neighborhood: 1. The move is infeasible with respect to the constraints, and ignored. 2. The move is feasible, not tabu, and is the best one found during this iteration. The move becomes the candidate move. 3. The move is feasible, not tabu, and is not the best one found during this iteration. The move is ignored. 4. The move is feasible, tabu, the best one found during this iteration, yet the aspiration criterion is not met. The move is ignored. 5. The move is feasible, tabu, the best one found during this iteration, and the aspiration criterion is met. The move becomes the candidate move. After evaluating the entire neighborhood, the resulting candidate move is incorporated into the solution and the neighborhood and tabu list are updated. This describes a general tabu search process with 1-opt moves (TS). Tabu search can generally find good solutions to large combinatorial problems, but they tend to be concentrated in a small portion of the solution space (Glover et al. 1995). To help alleviate this potential problem, a diversification routine was also incorporated into the heuristics. The diversification routine works like this: at some point in time during the search process, those moves with the lowest frequency of being scheduled are forced into the solution, thus moving the search process to another part of the solution space. One


where Un = the set of adjacent neighbors, which share at least a point, to unit n; k = adjacent neighbors to unit n; Tm = the set of near-time periods; m = a near-time period; Ak = area of unit k; Xkm = 0,1 variable indicating whether unit k is harvested during near time period m; An = area of unit n; and Xnt = 0,1 variable indicating whether unit n is harvested during time period t. For this problem, the maximum opening size used was 91 ha (225 ac). A maximum average opening size greenup constraint is included because it is the greenup requirement contained in the American Forest and Paper Association (2000) Sustainable Forestry Initiative (SFI). It can be defined for each time period. If each opening is centered on a focal unit f during a time period t, we can define the size of the opening Oft as:

tions, prevent this problem from being solved optimally using 0-1 integer programming techniques.

(Figure 1). The goal here is to act as a bridge between two good solutions with the hope of finding superior solutions based on combining two local optimal solutions. The best feasible local optimal solutions that were found from the last two 1,000 iterations of the 1-opt tabu search routines became the parents for the mating. The chromosomes for the genetic crossover can be visualized as a column array containing a solution for the problem. Each unit is a gene on that chromosome, and the timing of each unit’s harvest is the allele. A random crossover point was determined,

Assign local optimums to parent1, parent2

Randomly calculate crossover point, K

I = I+1

Is iteration, I >K

Yes

No child1[I]= parent1[I] child2[I] =parent2[I]

Is adjacency violated in child1 or child2

child1[I]= uncut child2[I] =uncut

Yes

No

No

child1[I]= parent2[I] child2[I] =parent1[I]

Is iteration, I < number of units

Yes

No Calculate objective function

Save best and return to normal tabu search procedure

Figure 1. Detailed flowchart of genetic crossover technique. Forest Science 48(1) 2002

41


thousand iterations of the 1-opt TS technique were completed before the diversification routine was activated. The whole TS process (1,000 iterations of 1-opt TS plus diversification) was then repeated six times before a final solution was reported. Five hundred feasible final solutions were generated for three hypothetical datasets. The second heuristic technique (TS/GA) combined tabu search with a GA crossover technique that was activated after two sets of 1,000 iterations of the 1-opt TS procedure and one iteration of the diversification routine

Description of Data Three hypothetical datasets and one operational dataset were created to compare the tactical plans produced by the solution techniques described above. Each hypothetical dataset consisted of a 25 row by 25 column grid, and each cell in the grid represented a 30 ha logging unit. The age classes for each unit were created with a random number generator. In two of these datasets, the ages were uniformly distributed between 0 and 30 yr. The third hypothetical dataset had an aggregated distribution of age classes. In the first five columns of the grid, the ages were uniformly distributed between 10 and 15 yr. The next 6–10 columns in the grid had ages uniformly distributed between 15 and 20 yr. Columns 11 through 15 in the grid had ages uniformly distributed between 20 and 25 yr. Finally, the last columns in the grid, columns 16–25, had ages uniformly distributed between 25 and 30 yr. The yields for these management units were a function of the age of the unit, and represented potential yields of radiata pine (Pinus radiata) (Goulding 1995). The minimum harvest age is 19 yr. A delivered log price of $100 per m3 was assumed, as was an 8% discount rate applied to costs and revenues in the middle of each time period. A second uniform random number generator was used to calculate logging costs, which ranged 42


between $20 and $70 per m3 for each management unit. Finally, deviations from the volume goal (generated by the strategic plan) were penalized at a base rate of $100 per m3; the penalty was discounted to the present from the middle of the time period in which it occurred, using a 9% discount rate. After formulating the problem and examining the decision variables, we found that the resulting planning problems consisted of between 3,000 and 5,000 0–1 integer variables. We acknowledge that these are simple problem sets, but they allowed us to develop good estimates of the optimal solution (the best integer solution found after 30 hr of computation using a branch and bound algorithm), as well as an upper bound to the integer optimal solution. The operational dataset, from a forest products company in Georgia, represents a typical industrial ownership in the southeastern United States. The dataset contains loblolly pine (Pinus taeda) plantations, as well as a mixture of hardwood areas, which, for this analysis, are not managed for timber production. It contains approximately 1,300 potential logging units on about 15,900 ha, resulting in approximately 13,000 0–1 decision variables. Acquiring real databases for testing scheduling techniques is somewhat problematic. In this instance, the work that was performed was done in conjunction with a forest management organization that wished to remain anonymous. We appreciate the assistance they provided in loaning us their data, and in return we agreed not to show a map of the property. We have included, however, a table describing the age class distribution of the forest (Table 1). The minimum harvest age was assumed to be 19 yr. Yields were estimated with equations developed by Harrison and Borders (1996) for three types of wood products: sawlogs, chip-and-saw logs, and pulpwood. The price for each product was also used as the penalty value for deviations from the volume goal. These were $131 per m3 for sawlogs, $91 per m3 for chip-and-saw logs, and $31 per m3 for pulpwood. All penalty values were discounted by using an 8% real discount rate.

Results Hypothetical Datasets All of the solutions generated by the heuristic techniques were feasible for the hypothetical datasets using the URM adjacency constraints; however, the quality of the solutions varied with respect to an estimation of the exact solution (from integer programming) and the relaxed LP solution. For example, the 1-opt TS technique found solutions between 93.7% to 97.6% of the estimated optimal value, and between 93.9% to 94.4% of the relaxed LP value. While one might hope for a better similarity to the “optimal” values, these results were comparable to results reported by others for 1opt techniques. Murray and Church (1996) reported that their best 1-opt TS solutions were within 94% of a global optimal. Bettinger et al. (1998) reported 1-opt TS solutions that were between 85% and 98% of an estimated optimal solution. While these results may seem problematic and unconvincing, they do indicate that a standard 1-opt technique applied to the hypothetical datasets and associated problem produces results consistent with previous research.


and the two chromosomes were split, and then recombined to form two new solutions, or children. From unit 1 to the random crossover point, K, the solutions from parent 1 are transferred to child 1, and the solutions for parent 2 are transferred to child 2. From unit K + 1 to n, the solutions from parent 1 are transferred to child 2, and the solutions for parent 2 are transferred to child 1. If the greenup constraints are violated as a result of the genetic crossover, a unit that contributed to the problem is unscheduled for harvest. The result is two feasible solutions that combine elements from two local optimal solutions. The child with the highest objective function value becomes the starting point for another iteration of the 1-opt tabu search algorithm. The third heuristic technique (TS2/GA) incorporated a 2-opt intensification technique introduced by Bettinger et al. (1999) to the TS portion of the heuristic technique. The technique was similar to the TS/GA technique. However, after 1,000 1-opt iterations, the best solution found from the previous 1-opt solution was used as the starting point for 200 iterations of a 2-opt TS procedure, where the harvest timing of one unit was switched with that of another unit (Figure 2). The diversification routine was then again activated while maintaining feasibility. The process looped back though 1,000 iterations of the 1-opt TS process, then another 200 iterations of the 2-opt TS process, before applying the genetic crossover technique. The resulting child with the highest objective function value became the starting solution for another loop through the TS processes. In sum, 3 TS loops were performed (each with 2 sets of 1,000 1-opt iterations and 200 2-opt iterations) as well as 2 genetic crossovers. After all of these processes were completed, a final (best) feasible solution was reported. Five hundred feasible solutions were generated for the 3 datasets described below, and 5 for the operational dataset.

Generate random start; memory = 1


Iteration = iteration + 1

Core tabu search procedure

Is iteration < 1000

Yes

No Intensification tabu search procedure

Is iteration < 200

Yes

No

Is Yes memory = 1 or 3 or 5

Long term No diversification routine

Iteration = 0 memory = memory+1

No

Is Yes memory = 2 or 4 or 6

Genetic crossover routine

Iteration = 0 memory = memory+1

No Print results Figure 2. Flowchart of 2-opt tabu search/genetic crossover technique.


43

Land area (ha) 1,343 3,356 2,352 3,490 2,741 1,846 595 95 7 23 47

Figure 3. Comparison of 1-opt tabu search, TS/GA, TS2/GA, the optimal value, and the LP value for hypothetical dataset 1.

allowing Cplex to run for 30 hr and subsequently reporting the best result. Operational Dataset A similar pattern was found when the heuristic techniques were applied to the operational dataset using the maximum opening size constraint and the maximum average opening size constraint. Here, however, an exact solution to the problem was not generated; thus, the basis for comparison is the best 1-opt TS result. The TS/GA located solutions that were about 2% better than the 1-opt TS technique. An additional increase of 2.4% was realized with TS2/GA. These percentage increases may seem small, but they result in plans with approximately a $200,000 increase in the objective function value ($15 ha–1 or $37 ac–1), when each new technique (genetic crossover and 2-opt intensification) is added to the TS technique (Figure 6). With the operational dataset and associated planning problem, the most important goal was the development of chip-and-saw log volume. The TS2/GA had the lowest average deviation from the chip-and-saw volume goals (4.7%), while the TS/GA was slightly higher (4.9%) and the 1-opt TS was highest (5.3%). The TS2/GA has the lowest average deviation for the sawlog goals (31.4%), but the deviations

Percent of optimal IP value

The TS/GA technique found consistently better solutions than did the 1-opt TS technique (Figures 3–5). Here, the best solutions had objective function values between 96.6% to 99.1% of the estimated optimal value, and between 96.3% to 96.4% of the relaxed LP value. Thus the addition of the genetic crossover technique improves the objective function value by about 2% over the 1-opt TS procedure. Now we find that by enhancing the heuristic, we can develop even better results than when using a standard 1-opt technique. In addition, when incorporating the 2-opt TS intensification technique, we further improved solution values by approximately 1.5 to 2%. With TS2/GA, the objective function values were between 98.2% and 99.7% of the estimated optimal value, and between 96.3% and 97.2% of relaxed LP value. For the datasets used here, the potential gain in net present value was approximately $213 ha–1 ($86 ac –1). Not only did the TS2/ GA technique produce significantly better values, but with one dataset, the lowest value obtained from TS2/GA was better than the best value found with the 1-opt TS. The motivation for using these more advanced routines is the generation of more efficient solutions that will allow a forestry organization to produce higher returns to its shareholders. The time required to generate a solution to the tactical problem formulations was 2 min., on average, for the TS and TS/GA techniques, and 4 min. for the TS2/GA technique. The IP formulation of the tactical problem was solved by


44




Age class 0–4 5–9 10–14 15–19 20–24 25–29 30–34 35–39 40–44 45–49 50+

Percent of optimal IP value

Table 1. Age class distribution for operational data set.

were much higher than for chip-and-saw goals due to the lower levels of sawlog production, just one-tenth of the chipand-saw production, and the impact of the spatial constraints that prevented these stands from being harvested. The 1-opt TS had a 34.7% average deviation, while the TS/GA had a 36.2% average deviation. The TS/GA did have the lowest average deviation for the pulpwood goals (9%), while the TS2/GA (9.4%) and 1-opt TS (10.0%) had slightly higher deviations. Overall, the more refined search techniques produced less deviation from the goals defined by the strategic forest plan. The time required to generate a solution to the tactical problem formulations was 30 min., on average, for the TS and TS/GA techniques, and 45 min. for the TS2/GA technique. The combinatorial nature of the operational problem with both maximum and average opening size constraints prevented solving the problem using mixed-integer programming.

Discussion Tabu search is considered one of the most promising heuristic techniques for solving combinatorial optimization problems, since it is able to solve both linear and nonlinear optimization problems. This research has also shown that tabu search can produce high quality solutions to planning problems with simple adjacency constraints as well as solutions to operational problems with greenup constraints in the form of maximum opening size and maximum average opening size. The 1-opt TS technique demonstrated here is commonly used to solve difficult planning problems. We feel that it found solutions that were similar to those found by many other techniques currently being used. However, since a 1opt TS technique limits the search region to the current solution plus only one change, it may not sufficiently explore the solution space. The genetic crossover technique allowed the search process to move to other areas of the solution space, and improved the ability of tabu search to develop good solutions. Using the 1-opt TS technique to quickly move to a good solution increased the efficiency of the GA by not requiring it to begin with a randomly generated, inefficient set of solutions. The refined search neighborhood


45


Figure 6. Comparison of 1-opt tabu search, TS/GA, and TS2/GA for the operational dataset.

created with the 2-opt TS technique further improved the quality of the solutions by fine-tuning the results. In terms of actual computer processing time when searching for solutions, all three heuristic techniques required less solution time than did the integer programming technique. The 1-opt TS and the TS/GA techniques required less than 2 min. to generate a solution for each hypothetical dataset. The TS2/GA technique required approximately 4 min. to generate a solution for each hypothetical dataset. Although these times are significantly faster than the 30 hr we required to find solutions using 0-1 mathematical programming, a comparison between the two solutions is not a fair one, as new formulation techniques such as using cliques or grid packing algorithms for the URM have shown to improve the speed at which these problems can be solved with 0-1 integer programming algorithms. The goal was to provide a means to compare quality of the solutions for the simple URM problems before turning to the rather complex greenup problems that more accurately portray the constraints that the forest industry encounter. For the operational dataset, the 1-opt TS and TS/GA techniques required about 30 min. to generate a solution, while the TS2/GA required nearly 45 min. to generate a solution. In general, the more complicated TS/GA and TS2/ GA techniques required only minor modifications to the 1opt TS program. The main differences were that for the GA portion of the process, the best local optimum solution had to be maintained for all iterations of the heuristic technique, whereas with the 1-opt TS technique, only the best solution was tracked. The 2-opt TS technique required additional programming logic to make two simultaneous moves, and it significantly increased the time required to generate a solution. However, the quality of the solutions improved. The gain from the additional programming work may thus outweigh the additional work required to implement a metaheuristic such as TS2/GA. Based on computing time alone, one may still ask whether the use of a heuristic to solve these types of problems is necessary. Five-hundred runs of one heuristic require about 16–33 hr of computing time, which, when compared to the 30+ hr to generate a 0-1 mathematical programming solution, doesn’t look too bad. However, since after testing the heuristic, we could develop with a statistical analysis a relationship indicating that a solution within x% of the global optimum can be developed after only y runs of the heuristic. Some organizations, after assessing the trade-offs related to solution quality and scheduling algorithm time requirements, may be willing to take a risk and develop fewer heuristic solutions in order to arrive at a high quality answer in a timely manner. Thus the use of a heuristic may seem worthwhile. What may be of interest to forest planners is the amount of time required to develop the infrastructure (computer code) to be able to use heuristic techniques. Unfortunately, a direct comparison of the heuristic development time and that of developing a matrix to solve the problems with linear or integer programming is unavailable, since we did not factor this comparison into the project at its inception. We believe, however, that the development time may be comparable to that of developing a matrix for integer programming.

Conclusions

Literature Cited AMERICAN FOREST AND PAPER ASSOCIATION. 2000. The 1999 Sustainable Forestry Initiative (SFI)SM Standard. (http://afandpa.org/forestry/sfi/ sfi_standard.html). Accessed November 14, 2000. BETTINGER, P., K. BOSTON, AND J. SESSIONS. 1999. Intensifying a heuristic forest harvest scheduling procedure with paired attribute decision choices. Can. J. For. Res. 29:1784–1792. BETTINGER, P., J. SESSIONS, AND K. BOSTON. 1997. Using tabu search to schedule timber harvests subject to spatial wildlife goals for big game. Ecol. Model. 42:111–123.

DUECK, G. 1993. New optimization heuristics: The great deluge algorithm and record-to-record travel. J. Comput. Physics 104:86–92. DUECK, G., AND T. SCHEUER. 1990. Threshold accepting: A general purpose optimization algorithm appearing superior to simulated annealing. J. Comput. Physics 90:161–175. GLOVER, F., J.P. KELLY, AND M. LAGUNA. 1995. Genetic algorithm and tabu search: Hybrids for optimization. Comput. Oper. Res. 22:111–134. GOULDING, C.J. 1995. Growth and yield models. In 1995 Forestry Handbook, Hammond, D. (ed.). New Zealand Inst. of For., Christchurch. HOF, J.G., AND L.A. JOYCE. 1992. Spatial optimization of wildlife and timber in managed forest ecosystems. For. Sci. 38:489–508. HOF, J., M. BEVERS, L. JOYCE, AND B. KENT. 1994. An integer programming approach for spatially and temporally optimizing wildlife populations. For. Sci. 40:177–191. HARRISON, W.M., AND B.E. BORDERS. 1996. Yield prediction and growth projection for site-prepared loblolly pine plantations in the Carolinas, Georgia, Alabama, and Florida. PMRC Tech. Rep. 1996. Plantation Manage. Coop., D.B. Warnell School of For. Resour., Univ. of Georgia, Athens, GA. HOGANSON, H., AND J. BORGES. 1998. Using dynamic programming and overlapping subproblems to address adjacency in large harvest scheduling problems. For. Sci. 44:526–538. KIRBY, M.W., P. WONG, W.A. HAGER, AND M.E. HUDDLESTON. 1980. A guide to the integrated resource planning model. USDA Management Sciences Staff, Berkeley, CA. 212 p. LOCKWOOD, C., AND T. MOORE. 1992. Harvest scheduling with spatial constraints: A simulated annealing approach. Can. J. For. Res. 23:468–478. MCDILL, M.E., AND J. BRAZE. 2000. Comparing adjacency constraint formulations for randomly generated forest planning problems with four ageclass distributions. For. Sci. 46:423–436. MITCHELL, M. 1995. Genetic algorithms: An overview. Complexity 1:31–39. MULLEN, D.S., AND R.M. BUTLER. 1997. The design of genetic algorithm based spatially constrained timber harvest scheduling model. In Proc. of the seventh symp. on Systems analysis in forest resources, Vasievich, J.M., et. al. (eds.). USDA For. Serv., North Central Res. Sta., St. Paul, MN. MURRAY, A.T. 1999. Spatial restrictions in harvest scheduling. For. Sci. 45:45–52. MURRAY, A.T., AND R.L. CHURCH. 1995. Heuristic solution approaches to operational forest planning problems. OR Spektrum. 17:193–203.

BETTINGER, P., J. SESSIONS, AND K.N. JOHNSON. 1998. Ensuring the compatibility of aquatic habitat and commodity production goals in eastern Oregon with tabu search procedure. For. Sci. 44:96–112.

MURRAY, A.T., AND R.L. CHURCH. 1996. Analyzing cliques for imposing adjacency restrictions in forest models. For. Sci. 42:166–175.

BEVERS, M., AND J. HOF. 2000. Spatially optimizing wildlife habitat edge effects in forest management linear and mixed-integer programs. For. Sci. 45:232–248.

NELSON, J., AND J.D. BRODIE. 1990. Comparison of random search algorithm and mixed integer programming for solving area-based forest plans. Can. J. For. Res. 20:934–942.

BORGES, J.G., H.M. HOGANSON, AND D.W. ROSE. 1999. Combining a decomposition strategy with dynamic programming to solve spatially constrained forest management scheduling problems. For. Sci. 45:201–212.

O’HARA, A.J., B.H. FAALAND, AND B.B. BARE. 1989. Spatially constrained timber harvest scheduling. Can. J. For. Res. 19:715–724.

BOSTON, K., AND P. BETTINGER. 1999. An analysis of Monte Carlo integer programming, simulated annealing, and tabu search heuristics for solving spatial harvest scheduling problems. For. Sci. 45:292–301. CLEMENTS, S.E., P.L. DALLAIN, AND M.S. JAMNICK. 1990. An operational, spatially constrained harvest scheduling model. Can. J. For. Res. 20:1438– 1447. COLORNI, A., M. DORIGO, AND V. MANIEZZO. 1998. Metaheuristic for high school timetableing. Comput. Optim. Appl. 9:275–288. COSTA, D. 1995. An evolutionary tabu search algorithm and the NHL scheduling problem. Inform. 33:161–178. CPLEX OPTIMIZATION, INC. 1998. CPLEX version 5. Cplex Optimization, Inc., Incline Village, NV.

46


REEVES, C.R. 1993. Genetic algorithms. P. 151–188 in Modern heuristic techniques for combinatorial problems, Reeves, C.R. (ed.). Wiley, New York. SNYDER, S., AND C. REVELLE. 1997. Dynamic selection of harvests with adjacency restrictions: The SHARe model. For. Sci. 43:213–222. VAN DUESEN, P.C. 1999. Multiple solution harvest scheduling. Silva Fennica. 33(3):1–9. WEINTRAUB, A., AND D. NAVON. 1976. A forest management planning model integrating silvicultural and transportation activities. Manage. Sci. 22:1299–1309. WEINTRAUB, A., F. BARAHONA, AND R. EPSTEIN. 1994. A column generation algorithm for solving general forest planning problems with adjacency constraints. For. Sci. 40:142–161.


The heuristics described here are early attempts at combining some desirable characteristics (intensification, diversification) of tabu search and genetic algorithms to create a hybrid heuristic technique. The goal is to continue to explore and develop processes that may produce efficient solutions to spatially constrained harvest scheduling problems. With the addition of both the genetic crossover and 2-opt TS intensification techniques, the results showed an increase of almost 3.5% in the objective function values over a 1-opt TS technique for small hypothetical problems, and a 4.4% increase in objective function values for an operational problem. The genetic crossover technique alone resulted in no significant increase in the time required to generate a solution, but the addition of the 2-opt intensification technique resulted in a significant increase. The increase in objective function value (for example, $200,000) can easily be justified based on the additional cost of waiting for the solutions (15 min.). These additional solution-generation time requirements are also small when compared to the time required to gather and organize the data required to solve forest planning problems. These combined heuristics, however, can further be modified to solve planning problems with more complex spatial relationships such as those described by Bettinger et al. (1997, 1998, 1999). We suggest that others consider heuristics to both vary and intensify the search process when developing spatially constrained harvest scheduling plans.

DAHLIN, B., AND O. SALLNAS. 1993. Harvest scheduling under adjacency constraints—A case study from the Swedish sub-alpine region. Scand. J. For. Res. 8:281–290.