Using Mixed Integer Goal Programming in Final Yield Harvest Planning: A Case Study from the Mediterranean Region of Turkey

: A mixed integer goal programming model is developed to address the regeneration planning problems of even-aged forests in the Mediterranean region of Turkey. The unique aspect of the goal programming formulation is to minimize deviations in scheduled wood product volumes and the size of harvest areas within each time period, as these are important goals for the management area. About 98% of the forests in Turkey are considered even-aged, and 2% are uneven-aged. Therefore, an age class method is used for the planning of even-aged forests. For the areas where this method is applied, reaching the optimal age class structure is the ﬁrst priority. This involves implementing ﬁnal harvests (clearcuts) to regenerate an amount of forest area into each age class. To meet the local market’s needs, forest enterprises also require the ﬁnal yield to be fairly equal each year. Further, it is desired that the harvest area (regeneration area) is relatively equal each year, to address operational considerations. A linear goal programming model is developed to address the problem. The minimization of deviations from both the harvest area and harvest volume targets are incorporated as goals in the objective function of the model. Several scenarios are solved using the extended version of Lingo 16. A scenario with weights of 0.8 for area and 0.2 for volume produces the best results. Here, the total deviation for 20 years is 3.8 ha in area and 2889 m 3 in volume. In the actual regeneration plan, the area deviation for 10 years is 54.72 ha (6.2% of total regeneration area), and the volume deviation is 20,472 m 3 (9.8% of harvest volume). The model described through this study can be developed further and integrated into forest management planning software and processes used for the planning of even-aged forests in the Mediterranean region.


Introduction
During the last two decades, Turkey has experienced considerable economic, environmental, and social developments which have increased the demands and expectations for forest resources. Until recently, forests have been predominantly the source of wood production, however, now, they are viewed as a source of non-wood forest products and provider of ecological and sociocultural services. In line with these developments, Turkey has changed its classical forest management planning model and implemented a functional ecosystem-based planning approach [1]. Modern harvest planning is needed not only to ensure that a plan of action provided to foresters is efficient from an economic perspective, but also that it recognizes as many of the quantifiable management issues as possible to forest management objectives that relate to commodity production, wood flow, economic returns, and perhaps other objectives [37][38][39][40][41][42][43][44][45][46]. More recently, spatial considerations have been incorporated into a GP problem. For example, [14] developed mixed integer GP models for thinning block designation based on the distance between scheduled blocks, and [47] and [48] developed processes for aggregated forest harvesting activities on the basis of their geographic location.
As a problem-solving methodology, GP is very useful in understanding contrasts in outcomes from the use of different assumptions in forest management planning processes [49], yet in one review, only 14 papers were located describing the use of GP in forestry from the onset of its use through about 2015 [50]. Many GP problem formulations use decision variables that are assigned continuous real values or that utilize data structures that preclude the ability to control the harvest of adjacent stands, e.g., [51]; relatively few examples, e.g., [7,14,52] of GP in forestry have utilized integer decision variables and attempted to control the timing and placement of harvests. Others have developed pairwise adjacency constraints to maximize age contrast between stands [51] and maximize clustering of stands within defined time periods [47]. Furthermore, while other research has proposed alternative GP formulations, we feel the type of integrated forest management planning issue described here has yet to be described in the forest management literature.
The objective of this study is to develop and assess a mixed integer multi-objective GP model that addresses spatial forest planning problems involving time-space arrangements of regeneration sites. Previous work of ours in this area involved creating thinning blocks using a mixed integer goal programming process [14]. One hypothesis for the current study is that GP can address a type of integrated forest management planning issue described above, which has yet to be assessed in the forest management literature. Another hypothesis is that through the GP modeling effort, more efficient final yield harvest forest plans can be developed, as compared to actual forest plans that have been developed through other means. Along these lines, a linear regeneration model is developed. The model aims to regenerate an equal amount of forest area every year to schedule an equal amount of product. This includes regenerating neighboring sub-compartments to final harvests that are smaller than 5.0 ha in size in the same year. In accordance with legislation associated with the management of the study area forest, areas applying a final harvest cannot exceed 25.0 ha during the greening (green-up) period. With this study, we also introduce the use of the greening (green-up) period approach in the final yield harvest planning of Turkish pine (Pinus brutia Ten.) forests in the region.
Without using a decision support system such as the one we propose, it is difficult to schedule the regeneration of an equal amount of area each year and to avoid regenerating forest areas larger than the allowable clear-cut size during the greening period. In practice, the regeneration activities are conducted by local foresters using detailed silviculture plans produced from forest management plans. In the best situation, area deviation may be more than 5% of the total regeneration area each year, and the volume deviation may be around 10% of the harvest volume each year. However, these deviations are acceptable in Turkish forest management. Avoiding clearing larger areas is crucial for ecological and environmental issues such as soil conservation, wildlife protection, and biodiversity conservation.

Study Area
The study area is located within the Akoren Planning Unit, in the southern part of Turkey ( Figure 1). The study area contains 5380 ha of coniferous forests, where Turkish pine is the dominant tree species. Other coniferous tree species include Anatolian black pine (Pinus nigra Arn. subsp. pallasiana (Lamb.) Holmboe var. pallasiana), fir (Abies cilicica Carr.), cedar (Cedrus libani A. Rich.) and juniper (Juniperus spp.). The growing stock of the forest is 559,006 m 3 and the annual increment is 22,422 m 3 . The main management objective of this coniferous forest is pine wood production. Since even-aged forest management is applied, every year the same amount of forest area is being regenerated and the same amount of wood flow is desired from the regeneration cuts (clearcuts).
Since even-aged forest management is applied, every year the same amount of forest area is being regenerated and the same amount of wood flow is desired from the regeneration cuts (clearcuts).

Study Area Data
While Turkish pine is the dominant tree species in the study area, Anatolian black pine stands (36 ha) and hardwood-coniferous mixed stands (27 ha) are also present. Only 5366.5 ha of productive forest (crown closure is more than 10%) of the planning unit is considered in this research, and 13.3 ha of degraded forest is ignored. The majority of the forest (55%) is considered site class II (medium). In this study, the site class differences have been neglected. It should be noted that the thinning operations in Turkish forestry are conducted at the compartment level and the regeneration operations are conducted at the stand (sub-compartment) level.
According to the current forest management plan for the study area [53], the annual allowable cut from final harvests is 23,350 m 3 , arising from 1777 ha over 20 years (88.85 ha per year). Additionally, the number of sub-compartments (stands) planned for regeneration is 166. According to the Communiqué on Technical Principles of Silvicultural Applications published by the Turkish Forest Service [54], final harvest areas shall not be larger than 25 ha in even-aged forests that have a production function. Therefore, stands larger than 25 ha were subdivided in ArcMap. For this study, the green-up period was assumed to be five years. During this time, it was not desired to regenerate the two adjacent compartments that comprise together an area larger than 25 ha. Likewise, to avoid regeneration activities in very small areas, it is assumed that adjacent stands smaller than 5 ha could be harvested together in the same year.

Problem Formulation
For the regeneration (final harvest) problem, the deviations from the harvest area (regeneration area) and harvest volume (final yield harvest) targets should be minimized. The foresters also want to cut each stand as a whole when they are entered, therefore, the decision for a stand is discrete (harvest or do not harvest). Furthermore, adjacent final harvests need to be spatially and temporally controlled. The nature of the management situation suggests that mixed integer multi-objective GP would be applicable to solve this problem.
The objective function of the linear regeneration model minimizes the deviations from wood volume and regeneration area targets over a 20-year time horizon. The constraints of the model are area constraints that prevent work in both small and very large areas during the greening (green-up) period. The mathematical formulation of the model is as follows:

Study Area Data
While Turkish pine is the dominant tree species in the study area, Anatolian black pine stands (36 ha) and hardwood-coniferous mixed stands (27 ha) are also present. Only 5366.5 ha of productive forest (crown closure is more than 10%) of the planning unit is considered in this research, and 13.3 ha of degraded forest is ignored. The majority of the forest (55%) is considered site class II (medium). In this study, the site class differences have been neglected. It should be noted that the thinning operations in Turkish forestry are conducted at the compartment level and the regeneration operations are conducted at the stand (sub-compartment) level.
According to the current forest management plan for the study area [53], the annual allowable cut from final harvests is 23,350 m 3 , arising from 1777 ha over 20 years (88.85 ha per year). Additionally, the number of sub-compartments (stands) planned for regeneration is 166. According to the Communiqué on Technical Principles of Silvicultural Applications published by the Turkish Forest Service [54], final harvest areas shall not be larger than 25 ha in even-aged forests that have a production function. Therefore, stands larger than 25 ha were subdivided in ArcMap. For this study, the green-up period was assumed to be five years. During this time, it was not desired to regenerate the two adjacent compartments that comprise together an area larger than 25 ha. Likewise, to avoid regeneration activities in very small areas, it is assumed that adjacent stands smaller than 5 ha could be harvested together in the same year.

Problem Formulation
For the regeneration (final harvest) problem, the deviations from the harvest area (regeneration area) and harvest volume (final yield harvest) targets should be minimized. The foresters also want to cut each stand as a whole when they are entered, therefore, the decision for a stand is discrete (harvest or do not harvest). Furthermore, adjacent final harvests need to be spatially and temporally controlled. The nature of the management situation suggests that mixed integer multi-objective GP would be applicable to solve this problem.
The objective function of the linear regeneration model minimizes the deviations from wood volume and regeneration area targets over a 20-year time horizon. The constraints of the model are area constraints that prevent work in both small and very large areas during the greening (green-up) period. The mathematical formulation of the model is as follows: where, i years (1, 20) j objectives: 1 = area scheduled for harvest, 2 = volume scheduled for harvest w − ij weight for negative deviations in objective j, year i w + ij weight for positive deviations in objective j, year i d − ij negative deviation in objective j, year i d + ij positive deviation in objective j, year i k regeneration stands X ki binary (0, 1) decision variable representing the harvest of stand k during year i A ki area available for harvest in stand k during year i V ki volume available for harvest in stand k during year i X mi binary (0, 1) decision variable representing the harvest of stand m during year i Pr pi pairs of adjacent stands AC i total area scheduled for harvest in year i AT i area target in year i VC i total volume scheduled for harvest in year i VT i volume target in year i Equation (1) is the objective function, which minimizes the deviations from harvest area and harvest volume targets, which were 88.85 ha and 23,350 m 3 per year. Equation (2), together with Equation (9), forces a stand to be scheduled for harvest. Equation (3) represents the accounting rows to add up area scheduled for harvest. Equation (4) represents the accounting rows to add up the volume scheduled for harvest. Equation (5), together with Equation (10), forces adjacent stands k and m, that are smaller than 5 ha, to be scheduled for harvest in same year. Equation (6) ensures adjacent stands k and m not be scheduled for harvest during the greening period (5 years). Twenty equations, represented by Equation (7), determine deviations in area scheduled for harvest. Twenty equations, represented by Equation (8), determine deviations in the volume scheduled for harvest.
The problem was solved using the extended version of Lingo 16.0 [55]. This version of Lingo allows an unlimited number of continuous value variables, integer variables, and nonlinear variables. Table 1 displays the total number of variables and constraints of our model. The problem was solved using different weights for the deviations in harvest volume and harvest area that ranged from 0.0 to 1.0, in 10 percent intervals, which created 11 scenarios (Table 2). These weights were selected in order to observe the sensitivity of changes in the volume and harvest area (regeneration area) outcomes when emphasis on the objective function elements was altered. The software was operated on a PC with a 2.60 GHz Intel ® Core™ i7 processor and 16 GB of RAM. A maximum of about 500 million iterations was allowed for the regeneration model. The computing times by scenarios can be found in Table 3. Table 1. Total number of variables and constraints used in the models.

Variables and Constraints Linear Regeneration Model
Variables 3620 Nonlinear variables 0 Integer variables 3500 Constraints 2767 Nonlinear constraints 0

Results from the Regeneration Model
Using the regeneration model, Lingo was able to produce a feasible solution to the stated regeneration problem for each of the 11 management scenarios. Scenarios 1 and 11 are the reference scenarios for the regeneration area and harvest volume, respectively, because they represent the greatest weights applied to the area deviations (scenario 1) and the volume deviations (scenario 11). A total of 210 h and 4 min and 59 s was required to run all 11 scenarios. The longest computing time was 35 h and 39 min and 2 s (scenario 4).
In the results, the total deviation from the target value for the area reference scenario (scenario 1) was only 3.8 ha for 20 years and the total deviation from the target value for the volume reference scenario (scenario 11) was 1025 m 3 . When the area and volume deviations were evaluated together, scenario 3 provided the best results. When the area and volume deviations were evaluated individually, respectively, scenarios 1, 3, 2, and 5, for area and 8, 11, 9, and 6 scenarios for volume provided better results (Table 4). It is interesting that scenario 8, with weights of 0.3 for area and 0.7 for volume, produced a 874 m 3 deviation in scheduled volume, which is better than the reference scenario (scenario 11). We allowed a maximum of about 500 million iterations to solve the problem. If we allow more iterations, the volume reference scenario would produce a better result. The actual area scheduled for regeneration for all 11 scenarios is displayed in Table 5 and the deviations from the area target (88.85 ha) in Table 6. However, in addition to the results of area reference scenario (scenario 1), the results of scenario 3, 2, and 5, which gave better results respectively, are presented in Figure 2.     Since the stands represent real forest data, the sizes of stands are all different. Therefore, a solution showing an exact adherence to the target (i.e., with 0 deviation from the target) is likely impossible. The area scheduled for harvest under the area reference scenario (scenario 1) deviated from the target by only 3.8 ha. The deviations here ranged from −0.95 (19th year) to +0.65 (1st year) ha per year. When the area and volume deviations were evaluated together, the total deviation was only 3.8 ha for Forests 2020, 11, 744 9 of 17 20 years for scenario 3, which provided the best result. The deviations from the target value range from −0.55 ha (4th year) to +0.55 ha (9th year). In other words, the maximum deviation from the area target value was 0.55 ha (0.6% of the annual regeneration area). In Turkish forest management, this deviation is very well acceptable.
As stated earlier, the area to be regenerated annually in the study area is 88.85 ha (1777 ha/20 years). If local foresters decide to implement the results of scenario 3, they annually will regenerate a forest area ranging from 88.3 ha to 89.4. This means that the maximum deviation will be only 0.55 ha, which is 0.6 per cent of the annual regeneration area. The regeneration operations in Turkey are very costly and it is not easy to find experienced forestry workers to employ each year to conduct harvesting and regeneration activities. If this deviation increases, this means the forest enterprise regenerates forests in different area sizes each year. This will affect the annual budgeting issues, operational efficiency, the needs of local timber markets, and the employment opportunities for experienced workers. More importantly, natural regeneration requires considerable attention. Quality seed collected from the same forest, intensive vegetation and litter clearing within the harvest area, and light tillage of the ground may be needed. After regeneration activities are completed, other forestry operations should be conducted such as seedling counting, vegetation management, and seedling maintenance. The success of all these operations is related to the annual regeneration area size. If the annual deviation from this target is minimal, the forest enterprise will be much more organized to conduct these regeneration operations.

The Results for Volume
The actual volume scheduled for harvest for all 11 scenarios is displayed in Table 7 and the deviations from the volume target (23,350 m 3 ) in Table 8. In addition to the results of volume reference scenario (scenario 11), the results of scenario 8, 9, and 6, which gave better results respectively, are presented in Figure 3.   Similar to the stand area issue, since the stands represent real forest data, the sizes of the stands are all different. Thus, the volume available in a stand is different. Therefore, a solution showing an exact adherence to the target (i.e., with 0 total deviation from the target volume) is likely impossible. For the volume reference scenario (scenario 11), the total deviation from the volume target was 1025 m 3 over 20 years. The deviations here ranged from −158 m 3 (13th year) to +166 m 3 (10th year). The total deviation for 20 years was 2889 m 3 under scenario 3, which provided the best result if we evaluate area and volume goals together. The deviations in volume scheduled in scenario 3 ranged from −1.047 m 3 (1st year) to 276 m 3 (20th year), and the ratio of deviations between the target volume and the scheduled volume varied between −4.48% and +1.18% per year. Again, this deviation is also  Similar to the stand area issue, since the stands represent real forest data, the sizes of the stands are all different. Thus, the volume available in a stand is different. Therefore, a solution showing an exact adherence to the target (i.e., with 0 total deviation from the target volume) is likely impossible. For the volume reference scenario (scenario 11), the total deviation from the volume target was 1025 m 3 over 20 years. The deviations here ranged from −158 m 3 (13th year) to +166 m 3 (10th year). The total deviation for 20 years was 2889 m 3 under scenario 3, which provided the best result if we evaluate area and volume goals together. The deviations in volume scheduled in scenario 3 ranged from −1.047 m 3 (1st year) to 276 m 3 (20th year), and the ratio of deviations between the target volume and the scheduled volume varied between −4.48% and +1.18% per year. Again, this deviation is also acceptable in Turkish forest management. In order to better illustrate the results of the linear regeneration model, a map was prepared for scenario 3, which again seemed to be the best solution achieved (Figure 4).  As noted earlier, the 23,350 m 3 annual final yield was set as the volume target. If local foresters decide to use the results of scenario 3, they will harvest wood ranging from 22,303 m 3 to 23,626 m 3 each year. This means the maximum deviation will be 1047 m 3 in the first year and the minimum deviation will be only 3 m 3 in the third year. If they prefer to use the results of scenario 8, the total volume deviation for 20 years will be only 874 m 3 . In this case, the regeneration area deviation will be 34.70 ha for 20 years, ranging from 4.05 ha (1st year) to −5.65 ha (16th year). Since this study doesn't evaluate the economic results of the management scenarios, we are not recommending the local foresters which scenario they would decide to use. However, the largest income source of forest enterprises in Turkey are timber sale revenues. Furthermore, around 90 per cent of the timber they produce comes from the final yield harvest (regeneration). In this case, minimizing the volume As noted earlier, the 23,350 m 3 annual final yield was set as the volume target. If local foresters decide to use the results of scenario 3, they will harvest wood ranging from 22,303 m 3 to 23,626 m 3 each year. This means the maximum deviation will be 1047 m 3 in the first year and the minimum deviation will be only 3 m 3 in the third year. If they prefer to use the results of scenario 8, the total volume deviation for 20 years will be only 874 m 3 . In this case, the regeneration area deviation will be 34.70 ha for 20 years, ranging from 4.05 ha (1st year) to −5.65 ha (16th year). Since this study doesn't evaluate the economic results of the management scenarios, we are not recommending the local foresters which scenario they would decide to use. However, the largest income source of forest enterprises in Turkey are timber sale revenues. Furthermore, around 90 per cent of the timber they produce comes from the final yield harvest (regeneration). In this case, minimizing the volume deviation is important to meet the needs of the markets in a sustainable manner.

Comparison of the Actual Regeneration Plan Data with the Results of the Regeneration Model
The actual regeneration plan [56], or detailed silviculture plan, for the study area indicated that over 10 years 881.8 ha of forests would be regenerated, and 209,600 m 3 of wood material would be produced (Table 9). This plan was developed using less sophisticated mathematical methods but took into account the growth of the forest. Given this information, we assumed that an area of 88.18 ha (881.8 ha/10 years) and a volume of 20.960 m 3 (209,600 m 3 /10 years) would be scheduled each year, since finer detail was unavailable from the plan. Using scenario 3 of the regeneration model, the total deviation for the first 10 years was 2.1 ha in the area (0.24% of total 10-year regeneration area) and 1959 m 3 in volume (0.84% of total 10-year harvest volume). Using the guidance from the actual regeneration plan, the area deviation would be 54.72 ha (6.2% of total regeneration area) from the assumed area target, and the volume deviation would be 20,472 m 3 (9.8% of harvest volume) from the assumed volume target (Table 10). This comparison suggests some efficiencies in the regeneration program can be realized using the mathematical programming approach.

Discussion
As all of the forests in Turkey are managed using plans prepared in consideration of basic principles such as sustainability, economics, productivity, multi-purpose utilization, protection of biological diversity, aesthetics and protection of other landscape values, and carbon balance, the use of appropriate techniques to achieve management objectives is important. Final yield and intermediate yield harvest plans are included in these management plans. The regeneration activities envisaged in the management plan are carried out according to detailed silviculture plans. In the long run, in order to harvest an equal amount of wood every year, it is necessary to cut a relatively equal area of forest area that will provide this each year. This is only possible in forests that are optimal in terms of their age class distribution. The way to meet these requirements is to establish an equal area of new stands every year and to manage them to the end of the rotation length. This approach, often considered as a regulated forest, should be taken into consideration when preparing the detailed silviculture plans. Furthermore, in accordance with legislation, areas applying a final harvest cannot exceed 25 ha during the greening (green-up) period. It seems impossible to regenerate an equal amount of area each year and to avoid regenerating forest areas larger than 25 ha during the greening period without using decision support systems. We have demonstrated in this case that a mixed integer GP model can be developed to provide efficient plans which can be used as guidance for implementation of activities by field foresters.
With a total of 3620 variables, the linear regeneration model presents a difficult problem. Both this situation and the use of 500 million iterations have greatly extended the Lingo's computing time. It took a total of 210 h and 4 min and 59 s to run all 11 scenarios. The longest computing time was 35 h and 39 min and 2 s (scenario 4). One of the factors directly affecting the computing time is undoubtedly the number of adjacency constraints, therefore one of the drawbacks of the model was that its long computing time needed to provide proper results when the adjacency constraints are numerous. While there is an inverse relationship between the clear-cut area size and the number of constraints, there is a linear relationship between the green-up period and the number of constraints. Conversely, as the green-up period increases, so generally does the number of constraints. Others recommend that the clear-cut size area be as large as possible and that the green up period is as short as possible to reduce the number of adjacency constraints [57].
We also explored the sensitivity of the model to produce the goals desired, by assessing how the results (harvest volume and regeneration area) changed when the weights applied to the deviational values varied. Managers seeking to locate forest planning solutions may desire that two (or more) goals have the same relative weight, and therefore one challenge to the planner is to determine the normalized weights that are appropriate for the problem [14]. A major issue of debate among the researchers using G concerns the use of normalization techniques to overcome incommensurability. Incommensurability in a weighted GP occurs when deviational variables are represented by different units (ha, m 3 , km, etc.). The summation of different units may cause an unintentional bias towards the objectives that have a larger magnitude. To overcome this difficulty, it may be necessary to divide each objective by a constant pertaining to that objective, to ensure that all objectives have roughly the same magnitude. Such a constant is known as a normalization constant. There are several normalization methods, each with its own normalization constant [58]. Some of them are percentage normalization [59], Euclidean normalization [60,61], summation normalization [62], and zero-one normalization [63]. However, these normalization methods cannot guarantee that the achieved objectives are consistent with their goals [64]. Therefore, with regard to the type of forestry problem addressed in this research, this area of investigation is still open.

Conclusions
In this study, we explained how one might use mixed integer GP in the Turkish forest management planning system to address final harvest issues. This study introduced a new approach to regeneration planning. This study provided an enhancement on previous research in this area as well. The previous research proposed an approach in intermediate harvest planning and stand tending block designation using GP models. The compartments were grouped together to create blocks by minimizing deviation between scheduled compartments. It was the first application of GP in intermediate harvest planning in Turkish forestry. The present study applied a GP model in final yield harvest planning, which is the other main aspect of forest management plans. This study can also be improved by integrating these different models into one comprehensive model to address both problems at the same time.
The biggest disadvantage of this study is that it does not take into account natural disturbances such as forest fires, insect and disease outbreaks, and drought and windthrow. This is not unusual, since most forest planning models ignore these issues, due to their high level of uncertainty. However, the process described has some major advantages over the current approaches to develop forest management plans in Turkey. Since the regeneration operations in Turkey are very costly and locating experienced forestry workers to employ to conduct harvesting and regeneration activities is not easy, using this type of model will help forest enterprises to regenerate their forests successfully. Some certainty and regularity of annual production levels can alleviate wood flow fluctuations and help forest enterprises increase their forestry operation efficiency. Further, they will be better positioned to meet the needs of local markets in a sustainable manner by regulating final yield harvest, a major source of income. Should decision-makers consider that this study be useful for forest planning and contribute to Turkish forestry, it would be of great benefit to develop a special matrix generator to solve these (and similar) problems. Finally, the approach can be tested in different planning units and can be employed world-wide in the planning of even-aged forests, primarily pine forests.