Optimal Planning of Wood Harvesting and Timber Supply in Russian Conditions

This paper describes an approach to the optimal planning of wood harvesting and timber supply for forest companies of Russia. Software and tools successfully used in other countries (e.g., Finland, Sweden, Canada, etc.) are not as effective in Russian conditions for a number of reasons. This calls for the development of an original approach to solve this problem with respect to Russia’s specific conditions. The main factors affecting the operation of wood harvesting companies in Russia were determined. The optimization problem was formulated taking into account all important features of wood harvesting specific to the country. The mathematical model of the problem was formulated and analyzed. An important requirement is that the solution algorithm should find high-quality plans within short computation times. The original problem was reduced to a block linear programming problem of large dimension, for which an effective numerical solution method was proposed. It is based on the multiplicative simplex method with column generation within Dantzig–Wolfe decomposition and uses heuristics to determine feasible solutions based on the branch and bound method. We tested the solution approach on real production data from a forest company in southern Karelia with a planning horizon up to a year. This case study involved 198 sites and 14 machines harvesting up to 200,000 cubic meters from an available stock volume of about 300,000 cubic meters. An increase in profit by 5% to 10% was observed, measured as revenue from the sale of products, net of harvesting and transportation costs.


Introduction
Currently, the active use of decision support systems based on mathematical modeling and operations research, along with geographical information systems (GIS), is considered the most effective way to increase production efficiency in the wood harvesting industry [1]. This approach minimizes the harvesting and transporting costs and leads to higher profits, and also improves other production and economic indicators without additional capital costs. Also, a new promising direction for better forest management focused on information and supporting economic, environmental, and sustainable decisions by using high technology sensing and analytical/digital tools, is the precision forestry. A detailed review of the contribution of the GIS, global navigation satellite systems (GNSS), and various kinds of sensors to forest operation can be found in [2].
Wood harvesting and timber transportation in Russia are typically conducted by wood-harvesting companies operating in leased forests. Forest companies have the opportunity to lease free forest areas from the state for up to 40 years, with possible further extension of the lease, subject to all difficult planning problems of a forestry company-the optimal planning of wood harvesting and timber supply, subject to aforementioned legislative and technological conditions and limitations specific for Russian forestry companies.
The paper is organized as follows. We describe the mathematical model and solution algorithm in Section 2. Section 3 contains numerical results of a case study. Section 4 considers some disadvantages of the approach. Section 5 concludes the paper.

Materials and Methods
The following main factors affect the operation of wood harvesting companies in Russia [18][19][20][21]: 1. Annual allowable cut (AAC): the annual volume of timber harvested in a leased forest area. Exceeding the AAC is strictly prohibited. It is regulated by state forestry authorities to ensure sustainable forest management; the volume of harvested wood should not exceed the estimated volume of the annual increase in the leased area. The AAC is assigned separately for production and protective forests, for coniferous and deciduous sites, for main and intermediate forest use, etc. All wood-harvesting companies strive to use the AAC to the maximum extent as far as their specific operational conditions allow.
2. The varying demand for wood harvesting products (assortments) in the operation locality: for example, demand for high-quality roundwood is generally higher than for low-quality firewood, etc. The current forest structure in Russia (natural mixed forests), in which wood-harvesting companies operate, quite often does not match the structure of customer demand. Thus, within the total supply volume there may be a shortage of some assortments and overproduction of others.
3. Accessibility of harvesting sites. One of the most serious problems of wood harvesting in Russia is the poor development of the forest road network; the average length of forest roads in the Russian Federation is 1.46 km per 1000 hectares, whereas in Western Europe and North America it ranges from 10 to 45 km. As the transporting cost accounts for a significant share of wood harvesting operation costs, it forces companies to determine the maximum distance from roundwood customers within which harvesting is economically feasible. For sites located beyond this distance, harvesting is unprofitable. If a new road must be constructed to access a harvesting site, its construction cost further reduces the cost-effective distance. Thus, forest areas that better match the demand structure may be economically unattractive.
4. Production facilities: the machinery, equipment, and labor employed by a wood harvesting company. They can be internal or external (hired).
Thus, to solve the problem of optimal planning of wood harvesting and timber supply on leased forest area, it is necessary to (1) determine the set of sites to be harvested and periods of harvesting on each one within the AAC and use available production facilities; and (2) determine, for each site in the plan, the composition of produced assortments and associate each cubic meter of each assortment to a customer to maximize profit.
At a tactical level, the typical planning horizon is 1 year. The solution to the problem requires such information as purchase prices for different types of assortments for different customers, harvesting and transporting costs, transporting distance, the accessibility of specific harvesting sites during the planning period, required cutting cycles, and the need for new road construction, among others.
More specifically, the main initial data for solving the problem include:

1.
A list of harvesting sites with locations, including connections to the road network, and necessary details for calculating the output volume of various assortments. From these sites, the algorithm will select the most profitable for harvesting, ceteris paribus.

2.
Consumer orders including the assortment type, delivery time, and price.

3.
Production facilities that may be involved in wood harvesting and transportation, including all necessary details, the main being productivity, cost, and location. 4.
Various technological limitations: AAC, required cutting cycles, machinery productivity, existing road network, seasonal accessibility of harvesting sites, road closing dates, etc. The initial data may be redundant, which is desirable for the problem solution quality. If the total volume of wood in a harvesting site equals AAC, the task will only be able to determine the harvesting periods and transportation routes. If the number of sites is greater, it enables choosing the most suitable ones. The situation is the same with orders and production facilities: if there is a choice, the plan will include orders and production facilities that result in maximal overall operational efficiency.
The main parameters of objects involved in the task are as follows: 1.
Planning parameters a. The planning period (start and end dates) b. The minimal period from clear-cutting the site (infrastructure cutting before road construction) to putting the road into operation (in days)

2.
Parameters related to harvesting sites a. The planned output of various assortments without specified length, calculated using the stock distribution by species and stand value class b. The felling method (clear-cutting, selective cutting, infrastructure cutting, thinning, etc.) c. Allowed operation periods (start and end dates), e.g., winter-only sites d. The earliest possible harvesting start date, if the site is adjacent to another site e. The obligatory indicator (0/1) that the site must be included in the plan (for whatever reason) f. Total operation cost, including harvesting, allotment, and road construction, taking into account the location of sand quarries g. The location on the road network map (latitude and longitude)

3.
Parameters of lease contracts a. AAC (in cubic meters) b. The maximum harvesting volume according to AAC separated by forest category (production, protective), felling method, and species (coniferous, deciduous) The goal is to find a set of sites the harvesting of which will yield the volume of wood matching customer orders by volume and structure, with the following optimization criteria and constraints. Parameters of the model: Planning period T G Time from clear-cutting to putting the road into operation V a AAC according to lease contract a ∈ A V ac Maximum total harvesting volume on sites complying with constraint c ∈ C for lease contract a ∈ A V k Maximum total harvesting volume on sites with felling method k ∈ K according to the specified productivity of wood harvesting machines Allowed supply period for order s ∈ S v s Required volume according to order s ∈ S for the supply period c s Price of wood for order s ∈ S r s ∈ Z + Priority of order Total volume of assortments on harvesting site d ∈ D T d Earliest possible harvesting start date for site d ∈ D according to required cutting cycles Total harvesting cost for site d ∈ D, including harvesting, allotment, and road construction, taking into account locations of sand quarries p h Productivity of wood harvesting machinery h ∈ H T e Operation start date for road section e ∈ E; L e Specified closing periods for road section e ∈ E : Let us build P as a set of all dates in the model: where {X} G = X, X + T G denotes the pair of dates, the initial date and the initial date offset by the period of time from clear-cutting to putting the road into operation. Then we reorder dates in P in ascending order, and P = {Θ t , t ∈ T}, where T is the index set of all time periods in the model with Θ t being the ending date of the appropriate period. In particular, • τ d , τ l1 , τ l2 , τ s1 , τ s2 , τ are the indices of dates T d , T l1 , T l2 , T s1 , T s2 , and T I . • ψ d and µ d are the indices of the start and end of harvesting periods for site d ∈ D, and µ g is the road construction completion date index (g ∈ G).
Using this set T, we append the list of model parameters: c dst is the transporting cost for one product unit from site d ∈ D to a customer demand point (terminal) for order s ∈ S during period t ∈ T, taking into account the operation start dates for road segments and closing dates, calculated based on the shortest distance from the harvesting site to the demand point for t ∈ T and specified transportation tariffs hr(t, h) is the number of hours of operation of forest machinery h ∈ H within period t ∈ T, taking into account its operation schedule (days per week and hours per day) To enable planning when the volume of harvested timber exceeds the volume of confirmed orders, we supplement the set of orders with a dummy (internal) order: S = S ∪ {o}, where o is the index of the internal order. For internal order, the following is true:  Unknown parameters: Indicator that site d ∈ D is included in the plan x dst Assortment volume produced at site d ∈ D according to order s ∈ S within period t ∈ T; z s Assortment volume delivered according to order s ∈ S.
With this notation, the mathematical model has the following formulation with two objective functions listed in order of priority: Objective function (2) maximizes the harvested volume for orders with higher priority: Coefficients p s , s ∈ S, accounting for the order priority, are determined as follows: Reorder the index set of orders S in ascending order of priority r s of order s ∈ S. Then calculate p s , s ∈ S so that p 1 = 1, p i = j<i p j v j + 1 to ensure that orders with higher priority are included in the optimal plan.
Objective function (3) subject to Constraints (4)- (14): Constraint (4) ensures that AAC is according to the lease contract: Constraint (5) ensures felling method limits: Constraint (6) Constraint (7) ensures matching consumer order volume: Constraint (8) respects harvesting machine productivity: d∈D s∈S Constraint (9) respects the volume of assortments available at harvesting sites: Constraints (10) ensure that harvesting is within the allowed time periods by setting to zero the volumes for time periods when harvesting is not allowed: prior to the allowed harvesting start date Constraint (11) ensures that obligatory sites are included in the plan: Constraint (12) ensures that required cutting cycles are respected: Constraints (13), with ψ d and µ d pointing to the start and end of harvesting at site d ∈ D, respectively, and µ g pointing to the adjacent road construction completion date, ensure continuity of harvesting operation at site d ∈ D and availability of adjacent road prior to the start of harvesting: Constraint (14) ensures non-negativity of volumes:
The large dimension is caused by constraint (9), reflecting the available volume of assortments on harvesting sites, and by constraints (12) and (13), accounting for the time parameters of harvesting sites and required cutting cycles.
The real production value of a wood-harvesting company (about 3000 harvesting sites per year and 50 different assortments) results in about 160,000 such constraints.
Note that constraints (10) on the equality of certain variables to zero do not increase the problem dimension. Similarly, constraint (11) does not affect the problem dimension, as it can be accounted for by appropriate correction of the right-hand side vector values.
One of the supporting optimization subproblems is to find the shortest paths on the road network graph. The following algorithms were tested: Floyd-Warshall algorithm, shortest path faster algorithm (SPFA), and Dijkstra's algorithm with a special pyramid data structure [25]. Expectedly, Dijkstra's algorithm with a pyramid showed the fastest calculation speed, as a typical forest road network has low density.

Solution Method
Applying linear programming methods to solving problem (1)- (14) is complicated by the following: The presence of two linear objective functions, ranked in order of priority, adds complexity to the problem [26][27][28]. It happens quite often that higher priority is assigned to a less profitable order of an important client, etc. These criteria will be satisfied as follows: In the first run, we solve problem (2), (4)- (14) by finding a set of columns that yield the optimum of objective function (2). Then, on this set of columns, we solve problem (3)- (14) to find the solution that yields the optimum of the objective function (3).

Boolean variables and nonlinear constraints
The multicriteria Boolean nonlinear programming problem is NP-hard, meaning that there are no effective exact solution algorithms of polynomial complexity [29]. The most common practical approaches to solving such problems are: 1.
Geoffrion relaxation with subsequent correction of the solution to relaxed problem [30] 2.

4.
Heuristic and metaheuristic algorithms (local search, taboo search, genetic algorithm, simulated annealing, etc.) [35,36] Following the Geoffrion scheme in general, the original problem is relaxed by excluding nonlinear constraints and discrete properties, with subsequent solution of the relaxed problem and its correction (in mathematical optimization and related fields, relaxation [30] is an approximation of a problem by a nearby problem that is easier to solve). However, the key point in the proposed approach is that nonlinear constraints and discrete properties are not removed from the model, but are taken into account using the column generation method described later.

Large dimension
The new relaxed problem (2), (15), (4)-(8), (16), (17) can be solved by linear programming methods. However, its large dimension in real production conditions (more than 100,000 constraints) does not allow the use of the standard simplex method.
Observe that the large dimension is mainly caused by constraint (16), reflecting the available volume of assortments at sites. However, the submatrix formed by this constraint has a block diagonal form [37], where each block corresponds to bounds on assortment volumes for a particular harvesting site. However, these blocks are independent (i.e., the intersection of their column sets is empty). Thus, we have the classical block linear programming problem [37] with coordinating the submatrix defined by constraints (4)-(8) and independent block submatrices, each of which corresponds to a harvesting site and has the form s∈S j t∈T x dst ≤ v dj , j ∈ N. However, this approach leads to a large number of block subproblems (equal to the number of harvesting sites), and thus, does not really help to tackle the large dimension. We propose combining the block problems of each harvesting site into other blocks as follows: Let U be the index set of blocks, u ∈ U. Harvesting sites with their adjacencies and other parameters of nonlinear constraints of the original problem are combined into one block. This makes it possible to consider nonlinear constraints (11)-(13) when solving the block subproblem using the column generation method, and not in the external coordinating problem. In practice, the number of blocks |U| can be estimated as √ |D|. Thus, optimization problem (2), (15), (4)-(8), (16), (17) can be formulated in the following block vector form: x u ≥ 0, u ∈ U.
The advantage of such representation lies in the possibility of applying Dantzig-Wolfe decomposition [38] with the column-generation method [39][40][41][42] to solve block LP problems. The method is based on transitioning to the set of extreme points α ω , ω ∈ Ω u of the subproblem corresponding to a block u ∈ U in the form of a convex hull x u = ω∈Ω u α ω λ ω , ω∈Ω u λ ω = 1. From now on, Ω u will be the index set of extreme points α ω satisfying constraint (20), and λ ω represents the unknown coefficients of the convex hull.
In this way, using the block representation of the new relaxed optimization problem (2), (15), (4)-(8), (16), (17) and Dantzig-Wolfe decomposition, we managed to substantially reduce the problem dimension and obtained problem (22)- (25), where (23) is related to the coordinating problem and (24) forms convex hull constraints for each block. In our practical case study presented in the Results section, the dimension was reduced from more than 100,000 constraints to about 150, with about 90 constraints of type (23) and about 60 of type (24).
We note that the matrix of problem (22)-(25) has a large number of columns, but there is no need to consider all extreme points for each block and solve problem (22)-(25) explicitly. Instead, it suffices to determine at each simplex method iteration just one point, which yields the maximum of the following LP problem, obtained from the optimality condition of the appropriate block subproblem: x u ≥ 0, where v is the vector of dual variables of the coordinating problem. Therefore, we obtain a supporting optimization problem (26)- (28) for the column generation method. It can be solved using the greedy algorithm [25], since all general constraints are taken into account in the coordinating problem. This allows us to account for nonlinear constraints of the original problem when solving the supporting problem in the framework of the column generation method. Thus, LP problem (22)-(25) is solved by using the column generation method with supporting problem (26)- (28) according to Dantzig-Wolfe decomposition.

Integer coefficients of the convex hull for the block LP problem
The optimal solution to problem (22)-(25) may include columns corresponding to different corner points for a certain block u ∈ U, which results in conflicting harvesting schemes. Thus, there is a need to choose a collaborative harvesting scheme for a block of harvesting sites. This is equivalent to choosing one corner point of the block problem. The corresponding constraint is: To account for this constraint, we propose a heuristic algorithm based on the branch and bound method within the standard scheme for solving integer linear programming problems. The algorithm has the following features:

•
The upper bound for an element is determined by the objective function value of the relaxed problem (22)-(25) excluding constraint (29). • A multilateral branching scheme is used, in which the element with the maximum upper bound is used as the branching element.

•
The branching of elements and solving of corresponding optimization problems are run in parallel. When the record is exceeded, it stops all running solutions whose upper bound is below the record. This significantly reduces the calculation time.
Since the solution tree in the branch and bound method can be quite deep, a heuristic was applied to save time: the calculation is stopped when a branch reaches certain depth n (the particular value is estimated in practice).
Thus, the original BNLP is solved by using linear programming algorithms with Dantzig-Wolfe decomposition within the column generation method, as well as the modified branch and bound algorithm to choose one corner point of the polyhedron for each block problem.

Implementation
The solution to the problem of optimal planning of wood harvesting and timber supply using the described method was implemented in Opti-Wood software, developed by the Russian IT company Opti-Soft together with researchers of Petrozavodsk State University [43] (Figure 1).
At the tactical level, Opti-Wood supports the planning of wood harvesting and timber supply, scheduling of wood harvesting and timber transport machines [10,11], planning of forest road construction and maintenance, scheduling of harvesting sites allocation, etc. At the operational level, the software supports the routing of timber trucks.
A core component of Opti-Wood is an optimization model that takes into account all important features of wood harvesting production in Russia and original optimization algorithms.
Another main component of Opti-Wood is the road network map with various objects attached to it: harvesting sites, intermediate storage areas, consumer demand points (storage), garages, terminals, quarries, etc. The spatial distribution of accessible harvesting sites (Figure 2a) and other objects (Figure 2b) strongly affects the calculation results. The accessibility of a harvesting site means both transport and organizational accessibility (i.e., the availability of information about the site needed for planning).
At the tactical level, Opti-Wood supports the planning of wood harvesting and timber supply, scheduling of wood harvesting and timber transport machines [10,11], planning of forest road construction and maintenance, scheduling of harvesting sites allocation, etc. At the operational level, the software supports the routing of timber trucks.

Results
We tested the methods and algorithms on a case study from a forest company in southern Karelia ( Figure 1) with a planning horizon up to a year. This case study involved 198 sites and 14 machines harvesting up to 200,000 cubic meters from an available stock volume of about 300,000 cubic meters. The calculation time on a personal computer (Intel Core i7-4790, 3.6 GHz CPU, 16 Gb RAM) in every case was under 5 minutes. Results are presented in Tables 1-3. 1 Observed profit at the forest company without using Opti-Wood software.

Results
We tested the methods and algorithms on a case study from a forest company in southern Karelia ( Figure 1) with a planning horizon up to a year. This case study involved 198 sites and 14 machines harvesting up to 200,000 cubic meters from an available stock volume of about 300,000 cubic meters. The calculation time on a personal computer (Intel Core i7-4790, 3.6 GHz CPU, 16 Gb RAM) in every case was under 5 min. Results are presented in Tables 1-3.   Therefore, the proposed numerical method and the algorithm implemented based on it raise the profit by 5% to 10% with given harvesting volumes over the annual planning horizon.
According to our modeling experience, the most significant factors for constructing an effective plan are the spatial locations of sites and consumers and the distribution of wood harvested on a site by product type, which depends on initial stand characteristics.
According to this optimal harvesting plan, other modules of Opti-Wood software were used to automatically calculate an optimal harvesting schedule, in which particular machinery is assigned to each harvesting site included in the plan for a precise operating period (Table 3, Figure 3). The objective of this problem was to minimize total machine relocation. A description of the applied heuristic and metaheuristic methods is provided in [10].  Therefore, the proposed numerical method and the algorithm implemented based on it raise the profit by 5% to 10% with given harvesting volumes over the annual planning horizon.
According to our modeling experience, the most significant factors for constructing an effective plan are the spatial locations of sites and consumers and the distribution of wood harvested on a site by product type, which depends on initial stand characteristics.
According to this optimal harvesting plan, other modules of Opti-Wood software were used to automatically calculate an optimal harvesting schedule, in which particular machinery is assigned to each harvesting site included in the plan for a precise operating period (Table 3, Figure 3). The objective of this problem was to minimize total machine relocation. A description of the applied heuristic and metaheuristic methods is provided in [10].

Discussion
The proposed approach allows the creation of cost-effective harvesting plans, taking into account specific features of wood harvesting in Russia. The practical feasibility of the created plans was confirmed by specialists of several forestry companies in Karelia and other regions of Russia.
As a disadvantage of the approach from the application standpoint, we note that it works with a limited list of harvesting sites and with fairly precise details, but not with the entire leased forest territory. Therefore, the solution is optimal only within the given list of sites. However, this approach

Discussion
The proposed approach allows the creation of cost-effective harvesting plans, taking into account specific features of wood harvesting in Russia. The practical feasibility of the created plans was confirmed by specialists of several forestry companies in Karelia and other regions of Russia.
As a disadvantage of the approach from the application standpoint, we note that it works with a limited list of harvesting sites and with fairly precise details, but not with the entire leased forest territory. Therefore, the solution is optimal only within the given list of sites. However, this approach is justified in practice, because most Russian companies do not have information about all of their leased forest territory with sufficient detail for adequate decision making. Moreover, this information is not always in digital form. Knowing precise details about harvesting sites is quite critical to ensure good quality of the plan: If incomplete or approximate initial data are supplied, the system would plan sales of nonexistent volumes of products or would not consume all harvested volume. Planning based on more detailed information is the next step toward more precise supply and demand matching.
So, this approach stimulates forest companies to verify information about harvesting sites before planning. Forest companies can use all available means of updating information. Whereas state forest inventory is carried out in a certain order, while waiting for their turn, companies often update information by themselves or by involving private contractors that use modern IT to collect and analyze stand data, e.g., data from satellite, airborne, unmanned aerial vehicles, global positioning systems, sensors, devices, analytical and communication tools [2], laser scanning and complex algorithms for processing raw data.
This work is usually performed gradually over several years. In addition, planning for the whole territory will significantly increase the number of objects, and consequently the calculation time.
Among other results, the proposed approach allows the identification of sites where harvesting would not be effective in current conditions. With such sites identified, forest management activities may be considered, e.g., thinning, etc. Automation and support of such planning may become one direction for further development of Opti-Wood software.
The achieved 5% to 10% increase in profit is generally in line with improvement rates reported by several other research teams [4]. It should be noted here that a one-to-one comparison of different approaches is difficult to undertake, since in every case results depend on particular production conditions and the quality of original "manual" planning. The models described in [3,[7][8][9] take into account many factors considered in the proposed model, however, these models fall short in accounting for Russian conditions as they omit such important factors as annual allowable cut, harvesting sites adjacency rules, and seasonal accessibility of harvesting sites.
As a computational disadvantage of this approach, we note the need to take into account integer convex hull coefficients λ ω , ω ∈ Ω u u ∈ U when using the modified branch and bound method. This leads to re-solving the original problem with an additional constraint at each step of the current set partition, which slows the calculation time. A possible remedy is the "warm" start, i.e., using the already found basis plan with a partial solution of the problem for new constraints. The effect of this modification on the solution quality of the main problem is another further research direction.

Conclusions
This paper describes the developed approach to optimizing wood harvesting planning and distribution of produced assortments between consumer demand points, taking into account a large number of legislative and technological conditions and limitations specific for Russian forestry companies.
A mathematical model of the problem is formulated and analyzed. This multicriteria Boolean nonlinear programming problem has a large dimension, two linear objective functions, and several nonlinear relationships between parameters.
The original problem was transformed to a block linear programming problem of large dimension. An effective numerical solution method is proposed based on a multiplicative simplex method with column generation within Dantzig-Wolfe decomposition, using heuristics to determine a feasible solution based on the branch and bound method. The proposed mathematical model and numerical methods for solving the described problem were implemented in the Opti-Wood software tool and tested on real production data from a forest company in southern Karelia. An increase in profit by 5% to 10% was observed, measured as revenue from the sale of products, net of harvesting and transportation costs.
Some disadvantages of the proposed approach are discussed, including the need to have a substantial amount of fairly detailed information about harvesting sites and stand characteristics, which is not always fully available. This lack of information is not specific to Russia, but can also be observed in other countries. Some managerial and technological measures are mentioned regarding how to identify and fill in the missing data.