INTRODUCTION

Average productivity of stands and its spatial distribution, considering the displacement of harvest operation team and the transport of logs, represent two of the main factors that contribute to the forest harvest performance. Therefore, the harvest schedule generation should consider these factors directly or indirectly.

The road network should be an integrated part of the forest harvest planning, what makes it a very complex task. The computational procedures of harvest optimization integrated to the road network can reduce the budget and the environmental impact caused by such activity. These mathematical tools allow the production of scenarios that represents the problem at hand, providing important information to support the decision-makings (^{MCDILL, 2014}).

^{Kirby et al. (1986}) were one the first authors to describe an integer formulation called Integrated Resource Planning Models (IRPM) of forest management, the class of operational research models that simultaneously integrates the problem of harvest and the configuration of the forest road network. ^{Guignard et al. (1998}) and ^{Andalaft et al. (2002}) described scientific approaches of IRPM formulation by the exact method through the Mixed Integer Linear Problems (MILP), using the Branch-and-Bound algorithm.

The shortest path algorithms, such as the dynamic programming algorithm called Floyd-Warshall (1962), have been used in formulation of IRPM models through MILP models. In these cases, the models are formulated including flow restrictions generated by the shortest path algorithms (^{AKAY et al. 2013}).

Based on the IRPM formulations, it is possible to measure the influence of several parameters related to forest yield, harvest dispersion and production unit cost. From this assumption, the following questions may arise: How much can these models influence timber production cost? What is the volumetric production in such scenarios? What are the computation tools to model such problems?

The purpose of this study is to apply decision making aid tools on the forest harvest planning in Eucalyptus spp. and Pinus spp. stands. The specific objectives were: 1) Minimize the production costs resulting from harvesting operations and forest road maintenance. 2) Simulate forest harvest scenarios with the adoption of different minimal harvest ages of Eucalyptus spp. and Pinus spp. stands. 3) Compare them according to the following performance parameters: harvest income, used road network and the production unit cost. 4) Check the influence of Floyd-Warshall shortest path algorithm in the performance parameters.

MATERIAL AND METHODS

Study area

The data considered for the present paper were obtained from a forest project of a Brazilian forest company located in the state of Paraná, southern Brazil. The total planted area covers 23,330 hectares, containing 2,455 stands of *Pinus* spp. and Eucalyptus *spp*. The main purpose of the stands is pulpwood production, the management regime is clear cut followed by replanting and the rotation age is of 8 and 16 years for Eucalyptus *spp*. and *Pinus* spp., respectively. Additionally, there are also older stands due to spatial location and other secondary production goals.

Table 1 presents the initial structure of the forest stands by age class considering the frequency of Pinus spp. and Eucalyptus *spp*. stands, respectively. This forest was used to validate the proposed mathematical model, but their applications are feasible in other forests too.

Cost definitions

The operational costs covered two activities: harvesting and road network maintenance. The harvest cost of each stand was calculated per tons of produced timber, ranging from 8 up to 16 R$/ton, including the procedures of cutting, processing, dragging and loading.

The road network was composed of unpaved roads (main and secondary) and paved roads. However, the road building cost depended only on the class of unpaved road. Fixed costs of road maintenance were established for two unpaved road classes: main road class and secondary road class, whose costs were equal to R$15,000/km and R$5.000/km, respectively. The construction of new road sections, in this study, was not considered.

Generation of flow routes

To establish fixed values for road network maintenance is not a viable condition because with one segment of main road it is possible to access many stands, resulting in combinatorial possibilities of options. Based on this assumption, different paths from the stands to the paved road were simulated based on the creation of graphs representing such components.

After the graphs buildings, the Floyd-Warshall minimal path algorithm was applied to generate the array of paths containing the sequences of road sections that represents the minimal path to be made from the stand up to the paved roads.

The algorithm was implemented using the software Visual Basic ^{(r)} for Applications (VBA), available in Microsoft^{(r)} Excel^{(r)} 2013. Following the structure presented by ^{Cormen et al. (1990}), where: is the weight of shortest path *i → j* with all intermediate vertices in {1, 2,... , *n*}. When k = 0, the path from the vertex *i* to the vertex *j* without intermediate vertices with numeration bigger than 0 is absolutely no intermediate vertex.

Considering that for any path all the intermediate vertices are in the set {1, 2,... , *n*}, the array D^{(n)}=(d_{ij}
^{(n))} provides the expected answer d_{ij}
^{(n)}
*= (i, j) for all i, j V*.

The resultant array of paths was used to formulate IRPM models of forest management.

Mathematical formulation

The mathematical model was formulated aiming to minimize the costs involved in the harvest and forest road maintenance subject to the production targets. A Integer Linear Program (ILP) model was proposed basing on the formulation of the Model Type I presented by ^{Johnson e Scheurman (1977}), for the annual harvest scheduling. Where *A*: total number of road sections; *U*: total number of stands; *T*: total number of annual planning periods; *G*: total number of graphs; *P*: total number of products; *Y*
_{ijk}: binary decision variable of harvesting of the stand *i*, in period *j,* contained on the graph *k*; *C*
_{ijk}: harvest cost of the stand *i*, in period *j,* contained on the graph *k; Z _{tjk}
*: binary decision variable of the road section

*t*, in period

*j*, contained on the graph

*k;*R

_{ijk}: maintenance cost of the road section

*i*, in period

*j,*contained on the graph

*k; I*: set of road sections that make up the path

_{tjk}*t*, in period

*j*, contained on the graph

*K*; V

*: amount of timber product*

_{pij}*p*, in stand

*i,*on the period

*j*; L

*: lower limit of timber production according to the annual demand of each product; L*

_{pi}*: upper limit of timber production according to the annual demand of each product.*

_{ps}The Objective Function was:

The set expression (2) requires each stand to be harvested only one time over the planning horizon. The constraints defined by (3) and (4) were formulated to ensure that production is within annual desired limits.

Aiming to compose the forest road network to be built by period, the constraints in (5) were formulated by the generation of array of paths from the Floyd-Warshall algorithm. Figure 1 presents the illustration of sample road network map and some paths.

The constraints in (6) impose that the decision variables assume only binary values.

Scenario optimization

The annual production goal defined by the company was 684 thousand and 964 thousand tons of timber for *Eucalyptus* spp. and *Pinus* spp. species, respectively. The growth predictions were determined following the production curves supplied by the company in function of the site, specie and age.

According to preliminary tests, each optimization scenario was simulated obeying 2% deviation margin of the timber production goals. The planning horizon corresponded to five years.

The specification of each scenario was first based on the variations of minimal harvest ages. In the Scenario I, we named the variation of the minimal harvest age as precocious age, in which the lower limits were equal to 4 and 6 years old for *Eucalyptus* spp. and *Pinus* spp., respectively. Regarding to Scenario II, we named its variation as semiprecocious age, whose lower limits were 5 and 10 years old, respectively. Finally, Scenario III represented the normal variation with lower limits equal to 6 and 12 years old, respectively. We called simulated scenarios the scenarios that contained all the model constraints.

Secondly, to enable new comparisons, three other scenarios without constraints (5) were formulated, corresponding to Scenario I* presenting precocious harvest ages, Scenario II* presenting semiprecocious harvest ages and Scenario III* considering normal harvest ages. We called the new scenarios group as witness scenarios.

After the formulations, GUROBI OPTIMIZER 6.0.4. was used to solve them and the processing limit time defined was 2 hours (7,200 seconds).

RESULTS AND DISCUTIONS

The optimization problems of forest planning developed in this study are considered complex due to the large number of constraints and decision variables involved in the formulation of ILP models, becoming in many cases impossible to obtain the optimal integer solution within the processing limit time defined. However, it was possible to obtain results very close to the best bounds. Table 2 presents the size of the MILP models, the objective function value and the distance in percentage of the best bound for the six scenarios. The best bound represents a value in which there is mathematical guarantee that the objective function value is not lower than it. Residual Gap is the relative difference between best known integer solution and best bound.

When compared to the witness scenarios, the simulated scenarios present highest number of variables and constraints, due to the presence of constraints (5).

The variations of the minimal harvest ages change the number of variables and constraints in inverse proportions, because the number of available stands with ages greater than to the established minimal ages decrease with increasing of minimal harvest ages that consequently reflect on decreasing of the variables and constraints number. Specifically, Scenarios I and I* encompasses the largest number of harvest opportunities due to its lower harvest ages (precocious), besides presenting the highest number of variables and constraints, when compared to the other scenarios group.

That number of variables and constraints in the simulated scenarios are considered high in relation to literature. As example, ^{Andalaft et al. (2002}) integrated the roads in forest harvesting considering 2-5 years planning horizon and annual periods, in which they formulated MILP models with more than 12,000 variables and 7,500 constraints, which resulted in solutions very close to the best integer solution.

Since the production varied ± 2% in the scenarios, the harvest area is inversely proportional to the productivity of the harvest scheduled stands. Figure 2 presents the cumulative harvest area by period and species for each scenario.

With the same minimal harvest ages, the witness scenarios, on average, presented harvest area 1.85% lower than the one obtained in the simulated scenarios, reflecting the non-utilization of constraints (5).

Among the simulated scenarios, Scenario III presented the lower total area for harvesting with 12,900 ha, a value 4.9% lower than Scenario II and 5.3% lower than the Scenario I, which corresponds to a discrepancy equal to 632 ha and 686 ha, respectively.

Eucalyptus spp. and Pinus spp. productivities were affected by the decrease of the minimal harvest age. However, by looking the harvest areas and objective function values of the scenarios, the production increase opposes the cost reduction. So by allowing harvest of young stands, with the consequent production decreases, the cost also decreases.

The road network distance resulted from the scenarios define the harvest dispersion in the evaluated periods. Figure 3 presents the sum of the road network distances built by period for each scenario.

As we expected, due to unintegrated models to the road network, the witness scenarios presented the road network with the longest distances, presenting a total average equal to 267.3 km, i.e. 60.2% higher than the simulated scenarios (106.1) km.

Among the simulated scenarios, Scenario I, which considers lower minimal harvest ages, presented the road network with lower distance, equals to 98 km, i.e value 1.2% lower than the Scenario II and 20% lower than the Scenario III, what corresponds to a discrepancy of 1.3 km and 20 km, respectively. Therefore, the larger road network distances in the Scenario III was consequence of the increasing of the minimal harvest ages that decreases the number of available harvest stands, resulting in increased road network density.

When the road network distance is decreased by planning period, the transportation costs and environmental impact tend to reduce. These are important aspects as highlighted by ^{Malinovsky (2010)}, where an average value of R$8,231.00/km was obtained for restoring roads with gravel. Besides, ^{Shan et al. (2009}) consider the building and restoration of forest roads an important issue due to its high potential environmental impact on the several involved ecosystems.

Table 3 present the total cost corresponding to the sum of costs of harvesting and road restoration, as well as the unit cost of timber production (R$/ton), which was calculated relating total cost of the projects (R$) and total production (ton).

On the simulated scenarios, the road maintenance cost decreases an average of R$1,174,382, i.e 60.4% of the witness scenarios. In the other hand, we observed increase of the harvest cost in the simulated scenarios, an average of R$959,144, i.e. 6.7% in relation to witness scenarios. This average difference was not enough to increase the production unit cost significantly on simulated scenarios.

Of the simulated scenarios, Scenario I presented the lowest unit cost with a value of R$9.69/ton, while for Scenario II and Scenario III the values were equal to R$9.74/ton and R$9.75/ton, respectively. For Scenario I, the decreasing of the minimal harvest ages increased the number of available stands for harvesting and consequently the possibilities of harvesting. Therefore, the optimal schedule planning becomes less restricted, resulting in lower unit production cost.

The road investment changes according to managed forest. The generation of harvest schedule guidelines requires integrated simulations like the ones presented in this paper and also other techniques, according to ^{Karlsson et. al (2006)}.

The heuristics techniques were used to design of road network for timber harvesting operations. ^{Weintraub (2000}) described two heuristic models: OPTIMED to decide road building options, and PLANEX, to locate access forest roads. ^{Akay (2004}) used Simulating Annealing to guide the search for minimizing the costs of construction, transportation, and maintenance of a forest road network. The results indicated that total road cost was reduced about 25%, but the main advantage is the short processing times of the technique.

CONCLUSIONS

The proposed method to generate flow routes of production applying the Floyd-Warshall algorithm is effective to formulate the IRPM models of forest management by means of ILP.

The formulated ILP models are efficient to minimize the production unit cost, considering harvest and forest road maintenance costs. On simulated scenarios, the unit production cost decreases an average of R$0.6, i.e 5.8% lower than the witness scenarios.

In the generated scenarios, the decreasing of the minimal harvest age reduces the production unit cost. In the evaluated scenarios such reduction is up to R$0.16/ton (1.5% less then).

The obtained integer solutions using ILP models are very close to the best possible solution, therefore, approximation methods are not required. This suggests that the methodology can be applied to forest harvesting integrated to the road network.