Next Article in Journal
Natural Disturbance Dynamics Analysis for Ecosystem-Based Management—FORDISMAN
Previous Article in Journal
Medicinal and Aromatic Lamiaceae Plants in Greece: Linking Diversity and Distribution Patterns with Ecosystem Services
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Optimal Planning of Wood Harvesting and Timber Supply in Russian Conditions

1
Department of Applied Mathematics and Cybernetics, Petrozavodsk State University, Lenin av., 33, 185910 Petrozavodsk, Russia
2
Department of Transport and Production Machines and Equipment, Petrozavodsk State University, Lenin av., 33, 185910 Petrozavodsk, Russia
*
Author to whom correspondence should be addressed.
Forests 2020, 11(6), 662; https://doi.org/10.3390/f11060662
Submission received: 27 April 2020 / Revised: 5 June 2020 / Accepted: 8 June 2020 / Published: 10 June 2020
(This article belongs to the Section Wood Science and Forest Products)

Abstract

:
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.

1. 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 the rules and norms of sustainable forest management. Forest areas are leased on a competitive basis. Preference is given to companies that have all the necessary resources to conduct business and have proven themselves to be responsible forest users during previous periods. If multiple or gross violations of the principles of sustainable forest management are identified, the lease term may be terminated prematurely.
The supply chain with fully mechanized cut-to-length (CTL) harvesting is used most often. The assortments produced are classified according to their future use, with main groups being sawlogs, pulpwood, energy logs, and logging byproducts. Each group is divided into several subgroups according to their species, qualities, and dimensions. Sawlogs, pulpwood, and energy logs after harvesting are transported by truck or railroad either directly to mills or to intermediate storage areas. The remaining logging residue can be used for the production of bioenergy.
Wood harvesting in natural stands highly limits the possibility of changing the production volume of particular kinds of products (assortments), as their volume distribution varies according to the harvesting site. Thus, finding consumers for all kinds of harvested products is a challenge. The daily changes in the production spatial structure, natural and production conditions, and consumer demand and other factors necessitate an effective decision support tool for wood harvesting planning.
To date, many works have been published that describe approaches to solving various problems related to the optimal planning of wood harvesting [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15]. The majority of these approaches are based on operations research, because corresponding methods, algorithms, and procedures, when used correctly, can yield notable economic and social effects.
An overview and comparison of many works related to planning problems and solution methods can be found in [1]. Early stages were characterized by the use of linear programming (LP) models in forestry planning. In those models, forest management was accounted for in constraints [15], but spatial constraints have not yet been considered. However, such tasks in real production conditions have a very large dimension to be solved using the standard simplex method, because with CTL harvesting on hundreds of sites, dozens of product types are produced, and each is assigned to a specific consumer, the number of which is also measured in the dozens.
A model for harvesting plan optimization with seasonal storage and road network deployment over a 2- to 3-year planning horizon was presented in [3]. In that study, the mixed-integer programming (MIP) problem was solved by strengthening the formulation and using Lagrangean relaxation.
It should be noted that planning of wood harvesting cannot be separated from planning of transportation and road maintenance. A MIP model to solve this multi-element planning problem was proposed by Karlsson et al. [8].
Another important challenge is integrating trucks and other modes of transportation, specifically ships and trains. The development of a new decision support system for transportation planning in Swedish forestry is reported in [7].
Some models were developed to support forest planning taking into account only consumer demand volume, but not prices, harvesting, transportation, and costs [12].
Some authors proposed an optimization model that integrates the assignment of machines to harvest areas and schedules the harvest areas during the year for each machine. The assignment problem is solved in phase 1 and the scheduling problem is solved in phase 2 [4,9].
Nevertheless, in Russia, computer-based decision support tools that automate wood harvesting planning are still used rather seldom. Software and tools successfully used in other countries (e.g., Finland, Sweden, Canada, etc.) [1,2,3,4,5,6,7,8] are not so effective under Russian conditions for a number of reasons, including significant differences in legislation and rules governing wood harvesting due to the state ownership of forest lands, the prevalence of natural untreated mixed stands, the uneven distribution of logging operations throughout the year due to seasonal availability of forest areas, own standards for round wood, different categories and generally poor condition of roads, specific organizational structure of wood-harvesting companies, the mixture of different technologies and machinery, etc. [16,17]. The aim of the research described in this paper is to develop an approach to solve the first and one of the most 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.

2. 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:
  • 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.
  • Consumer orders including the assortment type, delivery time, and price.
  • Production facilities that may be involved in wood harvesting and transportation, including all necessary details, the main being productivity, cost, and location.
  • 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:
  • Planning parameters
    • The planning period (start and end dates)
    • The minimal period from clear-cutting the site (infrastructure cutting before road construction) to putting the road into operation (in days)
  • Parameters related to harvesting sites
    • The planned output of various assortments without specified length, calculated using the stock distribution by species and stand value class
    • The felling method (clear-cutting, selective cutting, infrastructure cutting, thinning, etc.)
    • Allowed operation periods (start and end dates), e.g., winter-only sites
    • The earliest possible harvesting start date, if the site is adjacent to another site
    • The obligatory indicator (0/1) that the site must be included in the plan (for whatever reason)
    • Total operation cost, including harvesting, allotment, and road construction, taking into account the location of sand quarries
    • The location on the road network map (latitude and longitude)
  • Parameters of lease contracts
    • AAC (in cubic meters)
    • The maximum harvesting volume according to AAC separated by forest category (production, protective), felling method, and species (coniferous, deciduous)
  • Parameters of intermediate storage
    • Maximum capacity (in cubic meters)
    • Closing date
    • Location on the road network map (latitude and longitude)
  • Parameters of consumer orders
    • Assortment with specific length (in cubic meters)
    • Required volume (in cubic meters)
    • Delivery period (start and end dates)
    • Order priority (integer) reflecting the importance of matching the volume in the order, determined based on the customer rating, and typically highest for contracted orders
    • Demand point on the road network map (latitude and longitude)
  • Parameters of wood harvesting machinery
    • Maximum harvesting volume for each felling method within the planning period (in cubic meters)
    • Productivity (in cubic meters per hour)
    • Number of working hours per day and working days per week
  • Transporting costs (tariffs)
    • Allowed minimum and maximum transporting distances (in km)
    • Transporting costs for each assortment for specified distance
  • Road network, determined by the graph G = <V, E>, where E is the set of road segments and V is the set of map points, and:
    • Each road segment is characterized by the maximum speed, closing periods (e.g., wintertime roads), and operation start date.
    • Map points may be garages, harvesting sites, intermediate storage areas, demand points, terminals, etc.
    • The accessibility, minimum travel time, and minimum distance between map points during a certain period of time are calculated from graph G.
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.
Optimization criteria
  • Maximum supply is for orders of highest priority.
  • Maximum profit is calculated as sales price minus harvesting and transportation costs.
Constraints
  • AAC is not exceeded.
  • Maximum capacity of intermediate storage is not exceeded.
  • Maximum productivity of harvesting machinery is not exceeded.
  • Volumes are delivered within specified periods.
  • Sites are harvested within allowed operation periods.
  • Intermediate storage closing dates are respected.
  • Harvesting at each site goes on continuously.
  • Assortment volumes at harvesting sites are not exceeded.
  • Harvesting sites with set obligatory indicators are included in the plan.
  • Required cutting cycles are respected.
  • Infrastructure sites must be harvested earlier than adjacent non-infrastructure sites.

2.1. Mathematical Model

Index sets:
DSet of sites allowed for harvesting, d D
ASet of lease contracts, a A
ISet of intermediate storage areas, i I
CSet of constraints on harvesting volume according to AAC by forest category, felling method, and species, c C
KSet of felling methods, k K ;
G D Set of infrastructure harvesting sites (for road construction), g G
Q D Set of harvesting sites with set obligatory indicator, q Q
D g D Set of non-infrastructure sites adjacent to infrastructure site, g G
D d D Set of non-infrastructure sites adjacent to harvesting site, d D \ G
D a D Set of harvesting sites corresponding to lease contract, a A
D a c D Set of harvesting sites complying with constraint c C for lease contract, a A
D k D Set of harvesting sites with felling method, k K
D i D Set of harvesting sites attached to intermediate storage, i I
NSet of assortments (without specifying the length), j N
SSet of consumer orders, s S
S j S Set of orders for assortment, j N
HSet of wood harvesting machinery, h H
ESet of road segments, e E
Parameters of the model:
T P = [ T 1 P , T 2 P ] Planning period
T G Time from clear-cutting to putting the road into operation
V a AAC according to lease contract a A
V a c 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
V i Maximum capacity of intermediate storage i I
T I Closing dates of intermediate storage areas
T s = [ T s 1 , T s 2 ] 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 + Priority of order s S
v d j Available volume of assortment j N on harvesting site d D
v d = j N v d j 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
L d Allowed operation periods at site d D : { [ T l 1 , T l 2 ] } ,   l L d
c d 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 : { [ T e 1 , T e 2 ] } ,   l L e
Let us build P as a set of all dates in the model:
P = { T 1 P } G { T 2 P } G { T I } G { d D { T d } G } { e E { T e } G } { s S { { T s 1 } G { T s 2 } G } } { d D , l L d { { T l 1 } G { T l 2 } G } } { e E , l L e { { T e 1 } G { T e 2 } G } } ,
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 ,   τ l 1 ,   τ l 2 ,   τ s 1 ,   τ s 2 ,   τ are the indices of dates T d ,   T l 1 ,   T l 2 ,   T s 1 ,   T s 2 , 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 d s t 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
h r ( 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:
  • v o = Potentially unlimited volume
  • c o = Penalty for harvesting assortments in excess of the ordered volume
  • c d o t = Penalty for transporting assortments in excess of the ordered volume
  • T o = T P The delivery period matches the planning period
  • r o = min s S \ { o } r s 1 Internal orders have the lowest priority
Unknown parameters:
  • y d { 0 , 1 } Indicator that site d D is included in the plan
  • x d s t 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:
s S p s z s max
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) maximizes the revenue from sale of products, net of harvesting, and transportation costs
s S c s z s d D s S t T c d s t x d s t d D c d y d max
subject to Constraints (4)–(14):
Constraint (4) ensures that AAC is according to the lease contract:
d D a s S t T x d s t V a ,   a A , d D a c s S t T x d s t V a c ,   a A ,   c C
Constraint (5) ensures felling method limits:
d D k s S t T x d s t V k ,   k K
Constraint (6) respects intermediate storage capacity:
d D i s S t T x d s t V i ,   i I
Constraint (7) ensures matching consumer order volume:
z s = d D t T x d s t ,   z s v s ,   s S
Constraint (8) respects harvesting machine productivity:
d D s S x d s t h H p h h r ( t , h ) ,   t T
Constraint (9) respects the volume of assortments available at harvesting sites:
s S j t T x d s t = y d v d j ,   d D ,   j N
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 T d ; outside the allowed operation periods [ T l 1 , T l 2 ] at site d D ; outside the allowed supply periods [ T s 1 , T s 2 ] for order s S ; and after the closing date for intermediate storage T I :
x d s t = 0 ,   t τ d ,   Θ τ d = T d ,   d D ,   s S x d s t = 0 ,   t { [ τ l 1 + 1 , τ l 2 ] } ,   Θ τ l 1 = T l 1 ,   Θ τ l 2 = T l 2 ,   l L d ,   d D ,   s S x d s t = 0 ,   t τ s 1 ,   t > τ s 2 ,   Θ τ s 1 = T s 1 ,   Θ τ s 2 = T s 2 ,   d D ,   s S x d s t = 0 ,   t > τ ,   Θ τ = T I ,   d D i ,   i I ,   s S
Constraint (11) ensures that obligatory sites are included in the plan:
y d = 1 ,   d Q
Constraint (12) ensures that required cutting cycles are respected:
y d + q D d y q = 1 ,   d D \ G y g = d D g y d ,   g G
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:
ψ d = arg min t T { x d s t > 0 } ,   d D ,   s S μ d = arg max t T { x d s t > 0 } ,   d D ,   s S s S x d s t > 0 ,   t [ ψ d , μ d ] ,   d D μ g < ψ d + 1 ,   g G ,   d D g
Constraint (14) ensures non-negativity of volumes:
y d { 0 , 1 } ,   x d s t 0 ,   z s 0

2.2. Analysis of the Resulting Optimization Problem

Problem (2)–(14) is a multicriteria Boolean nonlinear programming problem (BNLP) [22,23,24] of large dimension with linear objective functions (2) and (3) and several nonlinear relationships between parameters (constraints (12) and (13)).
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.

2.3. Solution Method

Applying linear programming methods to solving problem (1)–(14) is complicated by the following:
1. Multiple criteria
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).
2. 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:
  • Geoffrion relaxation with subsequent correction of the solution to relaxed problem [30]
  • Brute force algorithms of the "branch and border method" class [31,32,33]
  • Cutting-plane methods (Gomory’s cut) [32,34]
  • 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.
Another outcome of the column generation method (described later) would be that for y d ,   d D , the objective function (3) is equivalent to objective function (15):
s S c s z s d D s S t T ( c d s t c d v d ) x d s t max
and constraints (9) and (14) are equivalent to (16) and (17), respectively:
s S j t T x d s t v d j ,   d D ,   j N
x d s t 0 ,   z s 0
Constraints (10)–(13) and integer variables y d ,   d D will be handled with the column generation method. The new linear programming problem is (2), (15), (4)–(8), (16), (17).
3. 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 d s t v d j ,   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:
u U c u x u max ;
u U B u x u = b ;
A u x u = b u ,   u U ;
x u 0 ,   u U .
Objective function (18) corresponds to objective function (2) or (3). Constraint (19) corresponds to the coordinating problem (constraints (4)–(8)). Constraint (20) corresponds to constraint (16), and constraint (21) ensures nonnegativity of variables (following (17)).
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.
Next, problem (18)–(21) is reformulated using variables λ ω ,   ω Ω u   u U :
u U ω Ω u σ ω λ ω max ;
u U ω Ω u W ω λ ω = b ;
ω Ω u λ ω = 1 ,   u U ;
λ ω 0 ,   ω Ω u ,   u U ,
where W ω = B u α ω ,   σ ω = c u α ω ,   ω Ω u ,   u U .
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:
( c u v B u ) x u max ;
A u x u = b u ;
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.
4. 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:
λ ω { 0 , 1 } ,   ω Ω u ,   u U .
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.
  • Branching is based on dividing the currently allowed set of variables λ ω into two subsets, λ ω = 0 and λ ω = 1 , for some ω Ω u ,   u U such that 0 λ ω 1 .
  • 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.

2.4. 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).

3. 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 Table 1, Table 2 and Table 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].

4. 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.

5. 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.

Author Contributions

Conceptualization, A.S. (Anton Shabaev) and A.S. (Anton Sokolov); methodology, A.S. (Anton Shabaev), A.S. (Anton Sokolov), and A.U.; software, A.U. and D.P.; validation, A.S. (Anton Shabaev), A.S. (Anton Sokolov), A.U., and D.P.; formal analysis, A.S. (Anton Shabaev), A.S. (Anton Sokolov), A.U., and D.P.; investigation, A.S. (Anton Sokolov) and A.U.; resources, A.S. (Anton Shabaev); data curation, A.S. (Anton Shabaev) and A.S. (Anton Sokolov); writing—original draft preparation, A.S. (Anton Sokolov) and A.U.; writing—review and editing, A.S. (Anton Shabaev) and D.P.; visualization, A.S. (Anton Sokolov); supervision, A.S. (Anton Shabaev); project administration, A.S. (Anton Shabaev); funding acquisition, A.S. (Anton Shabaev). All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. D’Amours, S.; Rönnqvist, M.; Weintraub, A. Using Operational Research for Supply Chain Planning in the Forest Products Industry. INFOR 2008, 46, 265–281. [Google Scholar] [CrossRef]
  2. Picchio, R.; Proto, A.R.; Civitarese, V.; Di Marzio, N.; Latterini, F. Recent Contributions of Some Fields of the Electronics in Development of Forest Operations Technologies. Electronics 2019, 8, 1465. [Google Scholar] [CrossRef] [Green Version]
  3. Andalaft, N.; Andalaft, P.; Guignard, M.; Magendzo, A.; Wainer, A.; Weintraub, A. A problem of forest harvesting and road building solved through model strengthening and Lagrangean relaxation. Oper. Res. 2003, 51, 613–628. [Google Scholar] [CrossRef]
  4. Bredstrom, D.; Jonsson, P.; Rönnqvist, M. Annual planning of harvesting resources in the forest industry. Int. Trans. Oper. Res. 2010, 17, 155–177. [Google Scholar] [CrossRef]
  5. Chauhan, S.S.; Frayret, J.M.; LeBel, L. Supply network planning in the forest supply chain with bucking decisions anticipation. Ann. Oper. Res. 2011, 190, 93–115. [Google Scholar] [CrossRef]
  6. Flisberg, P.; Frisk, M.; Rönnqvist, M. Integrated harvest and logistic planning including road upgrading. Scand. J. For. Res. 2014, 29, 195–209. [Google Scholar] [CrossRef]
  7. Forsberg, M.; Frisk, M.; Rönnqvist, M. FlowOpt—A decision support tool for strategic and tactical transportation planning in forestry. Int. J. For. Eng. 2005, 16, 101–114. [Google Scholar] [CrossRef]
  8. Karlsson, J.; Rönnqvist, M.; Bergström, J. An optimization model for annual harvest planning. Can. J. For. Res. 2004, 34, 1747–1754. [Google Scholar] [CrossRef]
  9. Simonenkov, M.; Salminen, E.; Bacherikov, I. An optimization model for monthly timber flow planning. J. Energy Resour. Technol. 2016, 13, 1–29. [Google Scholar] [CrossRef]
  10. Shabaev, A.; Sokolov, A.; Urban, A.; Pyatin, D. An approach to the scheduling of wood harvesting machines. IOP Conf. Ser.: Environ. Earth Sci. 2019, 316, 012061. [Google Scholar] [CrossRef]
  11. Shabaev, A.; Sokolov, A.; Urban, A.; Pyatin, D. An approach to the optimal timber transport scheduling. E3S Web Conf. 2020, 164, 03019. [Google Scholar] [CrossRef]
  12. Beaudoin, D.; LeBel, L.; Frayret, J.-M. Tactical supply chain planning in the forest products industry through optimization and scenario-based analysis. Can. J. For. Res. 2007, 37, 128–140. [Google Scholar] [CrossRef] [Green Version]
  13. Carlsson, D.; Rönnqvist, M. Supply chain management in forestry—Case studies at Södra Cell AB. Eur. J. Oper. Res. 2005, 163, 589–616. [Google Scholar] [CrossRef]
  14. Dems, A.; Rousseau, L.-M.; Frayret, J.-M. Effects of deferent cut-to-length harvesting structures on the economic value of a wood procurement planning problem. Ann. Oper. Res. 2013, 232, 65–86. [Google Scholar] [CrossRef]
  15. Goltsev, V.; Tolonen, T.; Syunev, V.; Dahlin, B.; Gerasimov, Y. Wood harvesting and logistics in Russia—Focus on research and business opportunities: Final report of the research project. Work. Pap. Finn. For. Res. Inst. 2011, 210, 157. [Google Scholar]
  16. Sokolov, A.P.; Gerasimov, Y.Y. Functional Logistics of Wood Harvesting Enterprise; Publishing House of PetrSU: Petrozavodsk, Russia, 2013; p. 84. [Google Scholar]
  17. Gunn, E.A.; Rai, A.K. Modelling and decomposition for planning long-term forest harvesting in an integrated industry structure. Can. J. For. Res. 1987, 17, 1507–1518. [Google Scholar] [CrossRef]
  18. Bulatov, A.F.; Voronin, A.V.; Kuznetsov, V.A.; Pladov, V.A.; Shegelman, I.R. Mathematical Models and Methods of Planning and Management of a Regional Forestry Enterprise; Publishing House of PetrSU: Petrozavodsk, Russia, 2000; p. 227. [Google Scholar]
  19. Gerasimov, Y.; Karjalainen, T. Development of wood procurement in Northwest Russia: Round wood balance and unreported flows. Eur. J. For. Res. 2006, 125, 189–199. [Google Scholar] [CrossRef]
  20. Gerasimov, Y.; Sokolov, A.P.; Fjeld, D. Improving Cut-to-length Operations Management in Russian Logging Companies Using a New Decision Support System. Balt. For. 2013, 1, 89–105. [Google Scholar]
  21. Gerasimov, Y.; Sokolov, A.P.; Syunev, V.S. Development Trends and Future Prospects of Cut-to-length Machinery. Adv. Mat. Res. 2013, 705, 468–473. [Google Scholar] [CrossRef]
  22. Bazaraa, M.S.; Shetti, C.M. Nonlinear Programming Theory and Algorithms; Mir: Moscow, Russia, 1982; p. 584. [Google Scholar]
  23. Hadley, G. Nonlinear and Dynamic Programming; Mir: Moscow, Russia, 1967; p. 509. [Google Scholar]
  24. Gorodetsky, S.Y.; Grishagin, V.A. Nonlinear Programming and Multi-Extreme Optimization; UNN Publishing: Nizny Novgorod, Russia, 2007; p. 489. [Google Scholar]
  25. Cormen, T.H.; Leiserson, C.E.; Rivest, R.L. Introduction to Algorithms, 3rd ed.; The MIT Press: Cambridge, MA, USA, 2009; p. 1292. [Google Scholar]
  26. Emelianov, S.V.; Larichev, O.I. Multi-Criteria Decision-Making Methods; Znanie: Moscow, Russia, 1985; p. 32. [Google Scholar]
  27. Kogan, D.I. Discrete Multicriteria Distribution Type Problems; UNN Publishing: Nizny Novgorod, Russia, 1991; p. 82. [Google Scholar]
  28. Podinovsky, V.V. Multicriteria Problems with Criteria Ordered by the Importance. Automat. Remote Control 1976, 11, 118–127. [Google Scholar]
  29. Savage, J.E. The Complexity of Computing; Factorial: Moscow, Russia, 1998; p. 368. [Google Scholar]
  30. Golshtein, E.G. Theory of Duality in Mathematical Programming and Its Applications; Nauka: Moscow, Russia, 1971; p. 352. [Google Scholar]
  31. Land, A.H.; Doig, A.G. An automatic method of solving discrete programming problems. Econometrica 1960, 28, 497–520. [Google Scholar] [CrossRef]
  32. Korbut, A.A.; Finkelshtein, Y. Discrete Programming; Nauka: Moscow, Russia, 1969; p. 368. [Google Scholar]
  33. Germier, Y.B. Introduction to Operations Research Theory; Nauka: Moscow, Russia, 1971; p. 384. [Google Scholar]
  34. Mitchell, J.E. Integer Programming: Branch and Cut Algorithms, Encyclopedia of Optimization; Floudas, C., Pardalos, P., Eds.; Springer: Boston, MA, USA, 2008. [Google Scholar] [CrossRef]
  35. Aho, A.V.; Hopcroft, J.; Ullman, J.D. Data Structures and Algorithms; Williams Publishing: Moscow, Russia, 2000; p. 384. [Google Scholar]
  36. Panteleev, A.V. Metaheuristic Global Extremum Search Algorithms; MAI-PRINT: Moscow, Russia, 2009; p. 160. [Google Scholar]
  37. Golshtein, E.G.; Yudin, D.B. New trends in Linear Programming; Soviet Radio: Moscow, Russia, 1966; p. 524. [Google Scholar]
  38. Lasdon, L.S. Optimization Theory for Large Systems; Nauka: Moscow, Russia, 1975; p. 432. [Google Scholar]
  39. Gass, S. Linear Programming: Methods and Applications; Fizmatgiz: Moscow, Russia, 1961; p. 304. [Google Scholar]
  40. Karmanov, V.G. Mathematical Programming; Fizmatlit: Moscow, Russia, 2005; p. 264. [Google Scholar]
  41. Romanovsky, I.V. Algorithms for Extreme Problems Solving; Nauka: Moscow, Russia, 1977; p. 352. [Google Scholar]
  42. Romanovsky, I.V. Discrete Analysis; Nevsky Dialect: Saint-Petersburg, Russia, 2016; p. 352. [Google Scholar]
  43. Opti-Wood. Available online: https://www.opti-soft.ru/en/opti/wood (accessed on 27 April 2020).
Figure 1. Interface of Opti-Wood software.
Figure 1. Interface of Opti-Wood software.
Forests 11 00662 g001
Figure 2. Spatial distribution of (a) harvesting sites and (b) other objects.
Figure 2. Spatial distribution of (a) harvesting sites and (b) other objects.
Forests 11 00662 g002
Figure 3. Harvesting schedule visualization.
Figure 3. Harvesting schedule visualization.
Forests 11 00662 g003
Table 1. Plan calculation results.
Table 1. Plan calculation results.
PlanTotal Volume Harvested (1000 m3) Total Cost
(million rubles)
Total Revenue
(million rubles)
Profit
(rubles/m3)
Profit (Actual) 1
(rubles/m3)
Increase
(%)
1167.9237.2441.5121711535.5
2172254.6515.9152014316.2
3145.2196.4428.9143113109.2
4128.1217.7388.4133212288.5
5186.8304.2569.5142013058.9
1 Observed profit at the forest company without using Opti-Wood software.
Table 2. Main cost categories (million rubles).
Table 2. Main cost categories (million rubles).
PlanTotalHarvestingTransportationRoad Construction
1237.2108.272.856.2
2254.611185.854
3196.492.169.332.5
4217.782.364.369
5304.2121.295.883.9
Table 3. Machinery utilization.
Table 3. Machinery utilization.
PlanTotal Volume Harvested (1000 m3)Hours WorkedFleet Utilization RateTotal Time for Machinery Relocation (hours)Highest Possible Productivity (1000 m3)
1167.919,1370.8863.41190.8
217219,0930.8964.44193.2
3145.216,2250.96143.64151.2
4128.113,6560.8877.46145.4
5186.820,3390.9554.55196.6

Share and Cite

MDPI and ACS Style

Shabaev, A.; Sokolov, A.; Urban, A.; Pyatin, D. Optimal Planning of Wood Harvesting and Timber Supply in Russian Conditions. Forests 2020, 11, 662. https://doi.org/10.3390/f11060662

AMA Style

Shabaev A, Sokolov A, Urban A, Pyatin D. Optimal Planning of Wood Harvesting and Timber Supply in Russian Conditions. Forests. 2020; 11(6):662. https://doi.org/10.3390/f11060662

Chicago/Turabian Style

Shabaev, Anton, Anton Sokolov, Alexander Urban, and Dmitry Pyatin. 2020. "Optimal Planning of Wood Harvesting and Timber Supply in Russian Conditions" Forests 11, no. 6: 662. https://doi.org/10.3390/f11060662

APA Style

Shabaev, A., Sokolov, A., Urban, A., & Pyatin, D. (2020). Optimal Planning of Wood Harvesting and Timber Supply in Russian Conditions. Forests, 11(6), 662. https://doi.org/10.3390/f11060662

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop