Submit manuscript...
eISSN: 2574-8092

International Robotics & Automation Journal

Research Article Volume 5 Issue 2

Control function optimization for stochastic continuous cover forest management

Peter Lohmander

Optimal Solutions in cooperation with Linnaeus University, Sweden

Correspondence: Peter Lohmander, Optimal Solutions in cooperation with Linnaeus University, Sweden

Received: February 13, 2019 | Published: April 18, 2019

Citation: Lohmander P. Control function optimization for stochastic continuous cover forest management. Int Rob Auto J. 2019;5(2):85-89. DOI: 10.15406/iratj.2019.05.00178

Download PDF

Abstract

Economically optimal management of continuous cover forests can be obtained via a new approach to adaptive control function optimization. We maximize our objective function, the expected present value, with consideration of stochastic prices, timber quality variations and dynamically changing spatial competition. The parameters of the control function are optimized via the first order optimum conditions based on a multivariate polynomial approximation of the objective function. The second order maximum conditions are investigated and used to determine relevant parameter ranges. The procedure is described and optimal results are derived for general function multi species CCF forests. A numerically specified model with empirical data showed that if the stochastic price variations are not considered when the harvest decisions are taken, the expected present value is reduced by 23%.

Keywords: economic optimization, multi species forests, forest management, stochastic processes, adaptive optimal control

Introduction

Economic forest management is an interesting area from a methodological point of view. Several dynamic and stochastic processes should be considered. Market prices are very important to optimal decisions, but often change rapidly and cannot be perfectly predicted. The central question is this: What is the best way to sequentially update the information and adaptively determine the management decisions?

If we are interested in economic results and market adapted forest management decisions, we cannot ignore prices, costs and other economically important parameters. It is necessary to consider the degree of predictability of future values of these parameters. Furthermore, in a dynamic world, the economically optimal levels of adjustments of management decisions to changes in prices, that are not perfectly predictable, are important. In stochastic markets, production capacity levels, stock policies and flexibility are important to the expected profitability.

Forest management decisions can be taken at many different levels. When the level of detail increases, the number of partial decision problems increases almost without bound. In continuous cover forestry, CCF, we may consider the management of each tree as a decision problem. Should we harvest this tree now or wait until some future point in time? Furthermore, these decision problems at the tree level are not independent. If one tree is harvested now, the available space increases for other trees in the neighbourhood to continue growing.

With large numbers of problem dimensions, stochastic parameter changes, large numbers of nonlinearities and adaptive decisions, the problem structure makes it difficult or even impossible to utilize standardized linear and nonlinear programming methods from the field of operations research. In order to give relevant solutions to real world problems, it is necessary to let the model contain the relevant structure with respect to how different parts of the analyzed system are connected and influence each other. One way to do this is to use stochastic simulation.

In the present paper, stochastic simulation will be used as a part of an adaptive control function optimization procedure. This approach has been developed by Lohmander.1 With stochastic simulation as a subroutine, it is possible to search the best way to control the system to reach the most desirable solution, in case the following procedure is utilized: First, a stochastic simulation model of the complete system under analysis is developed. The adaptive control of this system is defined via a specified control function.

General theoretical principles in the field of analysis can be used to define the functional form of the adaptive control function to be used in the system. The exact values of the optimal parameters of the control function are still unknown. Next, the complete system is simulated with a large number of alternative control function parameter value combinations. Thereafter, multidimensional regression analysis is used to determine an approximating function that gives the expected objective function value of the system as a nonlinear function of the control function parameters. Then, we maximize the value of the approximating function. From the first order optimum conditions, the optimal parameter values of the control function are determined. In case the approximating function is quadratic, the equation system of first order optimum conditions is linear. Then, the optimum is usually unique and it is possible that the approximating function can be shown to be strictly concave. If that is the case, the derived control function parameters give a maximum that is globally unique.

In case the approximating function is not quadratic, but for instance cubic, the analysis is more complicated. Then, the equation system of first order optimum conditions is not linear. Still, if the equation system only contains a limited number of nonlinearities, the solutions may be calculated via elimination and analytical methods of quadratic, cubic or quartic equations.

In such cases, it may be found that the approximating function is strictly concave in some region(s). Then, it may be possible to show that one of the solutions of the first order optimum conditions gives an optimum that is also a locally unique maximum. In some cases, it may be possible to show that we also have a unique global maximum. In this paper, this method will be utilized to derive optimal adaptive control functions in forest management.

The idea to use approximations, in particular quadratic approximations of the functions to be optimized, is not new. Sir Isaac Newton developed the approach denoted Newton’s method. Newton’s method is usually described in terms of root-finding, but it can also be understood as maximizing a local quadratic approximation to the objective function. Galantai2 describes the theory of the Newton’s method. The ideas have been extended in many directions. One such case is Wierzbicki,3 who developed a method for quadratic approximations based on augmented Lagrangian functions for non convex nonlinear programming problems. Another case is Powel,4 who created a new algorithm for unconstrained optimization models, leading to quadratic approximations via interpolations. Li5 made adaptive quadratic approximations to be used in structural and topology optimization and Lee et al.6 created an algorithm that gave local quadratic approximations for penalized optimization problems.

Optimization problems with many dimensions, nonlinearities and stochastic disturbances are common in most sectors of the economies. Lohmander7-9 presents alternative optimization methods to handle such situations in a rational way.

In this paper, we focus on CCF, continuous cover forestry. Initially, there are a large number of spatially distributed trees of different sizes. The central question is this: What is the best way to adaptively control and manage such a forest?

Lohmander10 shows that there is considerable option values associated with mixed forests. Single species forests give fewer options to adjust production to possible stochastic events. For instance, prices of different species may change. Then, it is valuable to have the option to adjust the harvest activities to these changing market conditions. Furthermore, some species may be negatively affected by pests, insects or large animals. Some species may not produce well in case the climate changes or if pollution increases. In these cases, it is valuable to be able to adaptively adjust forest production. This can easily be done if we already have several species growing in the forest. In some cases, it is possible to calculate the expected present value of forestry, conditional on the initial species mix. Lohmander7 contains several optimization methods and typical solutions to adaptive forest management problems.

Methods

In the present study, we develop and describe a general analytical and numerical method to handle management decision problems of this type: We want to optimize the harvest decisions over time. We want to maximize our objective function, the expected present value. The prices of the different species are stochastic. The problem is solved using an adaptive control function. The parameters of the control function are optimized via the first order optimum conditions of an approximation of the expected objective function. The second order maximum conditions are investigated. The expected objective function is estimated via Monte Carlo simulation.

In this section, a general procedure is given for a multi species forest with trees of different sizes. Now, we consider a mixed species forest. Initially, there is a large number of spatially distributed trees of different sizes and species. We want to optimize the harvest decisions over time. We want to maximize our objective function, the expected present value. The prices of the different species are stochastic.

bi(0)=bi0ibi(0)=bi0i (1)

bi(t) is the basal area of tree number i(at height 1.3 meters above ground) at time (t)(from t0=0t0=0 ). Each period normally represents one year but other time intervals are sometimes more relevant. The initial condition is bio.di(t) is the diameter of tree i at time t at height 1.3 meters above ground.

di(t)=2(bi(t)π)di(t)=2(bi(t)π) (2)

u(i,t) represents the control decision. If u(i,t)=1, the tree i is harvested in period t. Otherwise,u(i,t)=0.

u(i,t){0,1}i,tu(i,t){0,1}i,t (3)

Tt=0u(i,t)1iTt=0u(i,t)1i (4)

bi(t) develops according to a discrete time process. The increment, growth,, is a function of the basal area bi(t), the species s(i) and the competition,L(i,t).S(i){1,.,N}S(i){1,.,N} where N is the total number of species. is a ”species dummy variable” defined in (5).

sm(i)={1ifS(i)=m0ifS(i)mi,m{1,.,N}sm(i)={1ifS(i)=m0ifS(i)mi,m{1,.,N} (5)

The level of competition of relevance to tree i at time t is denoted L(i,t). In different studies,L(i,t) has been specified in different ways. Usually, L(i,t) is a strictly increasing function of the size (for instance basal area) of neighbour trees. Furthermore, neighbour trees that are close to tree influence the value of L(i,t) more than what more distant trees do. In Lohmander,1 L(i,t) is the total basal area of neighbour trees per hectare within a circle of radius 10 meters, where tree   represents the center of the circle. In Lohmander et al.,11 L(i,t) is a nonlinear function of the properties of the competitors; L(i,t) decreases with the distance to the competitors and increases with the size of the competitors.

The general function G(.), for L(i,t)=0, has been defined and presented by Lohmander.12 Lohmander also derived a differential equation consistent with G(.) and the dynamic properties of the basal area development. Empirical estimations of the parameters of G(.) with variations of L(i,t) have been performed for forests with different tree species in Iran by Hatami et al.13

bi(t+1)={bi(t)+G(bi(t),S(i),L(i,t))foru(i,t)=00foru(i,t)=1bi(t+1)={bi(t)+G(bi(t),S(i),L(i,t))foru(i,t)=00foru(i,t)=1 (6)

The present value of all profits is denoted Z. This is the discounted net value of all harvests. Hence it is a function of all harvest decisions, the rate of interest, the prices of the different species, the volumes of trees at different points in time and the harvest costs.

Z=Tt=0ertIi=1u(i,t)(P(S(i),t)V(S(i),bi(t))C(S(i),bi(t),t)))Z=Tt=0ertIi=1u(i,t)(P(S(i),t)V(S(i),bi(t))C(S(i),bi(t),t))) (7)

ertert is the discounting factor of period t with rate of interest r.P(S(i),t)P(S(i),t) denotes price per cubic meter of specie i sin period t.P(S(i),t)P(S(i),t) t is a stationary variable which is stochastic at tn,n>0tn,n>0 .E(P(S(i),t))E(P(S(i),t)) is the expected value of P(S(i),t)P(S(i),t) at tn,n>0tn,n>0 . In the two species case, if trees i and j belong to different species, then P(S(i),t)P(S(i),t) and P(S(j),t)P(S(j),t) have correlation ρρ .(1ρ1).V(S(i),bi(t))(1ρ1).V(S(i),bi(t)) is the volume of tree as a function of the species and the basal area.C(S(i),bi(t),t))C(S(i),bi(t),t)) denotes the harvest cost of tree i. This cost is a function of species, basal area and time.

This problem is highly stochastic, multidimensional and nonlinear. Furthermore, it contains a large number of integer variables. It is necessary to define a reasonable type of adaptive control function that can be used to handle the many control decisions in a way that takes the stochastic prices and competition between trees into account. Then the parameters of the control function may be optimized. For this purpose, the following rule is suggested: First, we calculate the ”limit diameter”Di(t) of tree i at time t. The limit diameter is a function of the tree species index, the relative deviation of the price from the expected level and the competition.αmαm is the value of the limit diameter Di(t) if the species is m, and at the same time, in case the price is a the expected level and L(i,t)=0.

Di(t)=Nm=1αmsm(i)+αP(P(S(i),t)E(P(S(i),t))E(P(S(i),t)))+αLL(t,i)Di(t)=Nm=1αmsm(i)+αP(P(S(i),t)E(P(S(i),t))E(P(S(i),t)))+αLL(t,i) (8)

In case the diameter di(t)di(t) is larger than the limit diameter, we instantly harvest. Otherwise we wait at least one more period before we harvest the tree. (The particular functional form (8) can be generalized in several ways. For instance, the values of the parameters αPαP and αLαL may be different for different species.)

u(i,t)={0if(di(t)Di(t))(1t1t=0u(i,t))01if(di(t)Di(t))(1t1t=0u(i,t))>0i,t{0,κ,2κ,...,nκ}u(i,t)=⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎪ ⎪ ⎪ ⎪ ⎪ ⎪0if(di(t)Di(t))(1t1t=0u(i,t))01if(di(t)Di(t))(1t1t=0u(i,t))>0i,t{0,κ,2κ,...,nκ} (9)

κκ denotes the harvest decision interval. This is an integer,κ1κ1 .n is the total number of harvest decision intervals and nκ=Tnκ=T , where T is the total planning horizon.

Z is a function of many things, including the stochastic price outcomes. When the control decisions are optimized, we are interested in the expected value of Z, namely E(Z), which is defined in (10). We may estimate E(Z) for given parameter values via Monte Carlo simulation. The average value of Z is determined based on a large number of random outcomes of the stochastic prices of the different species. The correlations may be estimated from real price series and Cholesky factorization can be used to generate the correlated price series. It is a good idea to store all of the simulated price series and to use the same set of simulated price series in every step of the control function parameter optimizations.

E(Z)=E(Tt=0ertIi=1u(i,t)(P(S(i),t)V(S(i),bi(t))C(S(i),bi(t),t)))E(Z)=E(Tt=0ertIi=1u(i,t)(P(S(i),t)V(S(i),bi(t))C(S(i),bi(t),t))) (10)

In the following derivations,N is assumed to take the value 2, which is a typical case in real applications. The procedure can easily be extended to other values of N.

Procedure:

Make an initial guess (α10,α20,αP0,αL0)(α10,α20,αP0,αL0) of the optimal values of the parameters (α1,α2,αP,αL)(α1,α2,αP,αL) .

Create a number, W, of alternative parameter combinations, w , such that stochastic variables (ε1,ε2,εP,εL)(ε1,ε2,εP,εL) are added to (α10,α20,αP0,αL0)(α10,α20,αP0,αL0) . The probability density functions of these stochastic variables are defined with consideration of the interesting parameter ranges.(djλ,djγ)(djλ,djγ) denote the lower (λλ ) and upper (γγ ) bounds ofεjεj, j{1,2,P,L}j{1,2,P,L} .

d1λ<ε1<d1γ;d2λ<ε2<d2γ;dPλ<εP<dPγ;dLλ<εL<dLγd1λ<ε1<d1γ;d2λ<ε2<d2γ;dPλ<εP<dPγ;dLλ<εL<dLγ . Normally, this means that d1λ<0<d1γ;d2λ<0<d2γ;dPλ<0<dPγ;dLλ<0<dLγd1λ<0<d1γ;d2λ<0<d2γ;dPλ<0<dPγ;dLλ<0<dLγ . Let εjwεjw denote the value of εjεj in parameter combination w.(α1w,α2w,αPw,αLw)=(α10+ε1w,α20+ε2w,αP0+εPw,αL0+εLw)(α1w,α2w,αPw,αLw)=(α10+ε1w,α20+ε2w,αP0+εPw,αL0+εLw) .

For each random parameter combination w, estimate the value of E(Z)=E(Z(α1w,α2w,αPw,αLw,...))E(Z)=E(Z(α1w,α2w,αPw,αLw,...)) .

Make a quadratic approximation Φ=Φ(α1,α2,αP,αL)Φ=Φ(α1,α2,αP,αL) of the function E(Z)=E(Z(α1,α2,αP,αL,...))E(Z)=E(Z(α1,α2,αP,αL,...)) according to the lines suggested by equations (11) – (14), via OLS, the ordinary least squares method.

It is important to be aware that cubic approximations or other functional forms may sometimes be more relevant.

Φ=Φ(α1,α2,αP,αL)Φ=Φ(α1,α2,αP,αL) (11)

Φ=Φ(α1,α2,αP,αL)=k1α1+k2α2+kPαP+kLαL++k11α12+k22α22+kPPαP2+kLLαL2++k12α1α2+k1Pα1αP+k1Lα1αL+k2Pα2αP+k2Lα2αL+kPLαPαLΦ=Φ(α1,α2,αP,αL)=k1α1+k2α2+kPαP+kLαL++k11α12+k22α22+kPPαP2+kLLαL2++k12α1α2+k1Pα1αP+k1Lα1αL+k2Pα2αP+k2Lα2αL+kPLαPαL (12)

Φ(α1,α2,αP,αL)E(Z(α1,α2,αP,αL,...))Φ(α1,α2,αP,αL)E(Z(α1,α2,αP,αL,...)) (13)

mink0,k1,...,kPLΨ=E[(Φ(α1,α2,αP,αL;k1,...,kPL)E(Z(α1,α2,αP,αL,...)))2] (14)

Use the quadratic approximation Φ=Φ(α1,α2,αP,αL) to determine the approximately optimal values of the parameters α1,α2,αP,αL . The approximate optimal values can be used as new initial conditions, and the approximation process can continue any number of iterations until the solution is considered satisfactory.

The first order optimum conditions are found in (15).

{dΦdα1=k1+2k11α1+k12α2+k1PαP+k1LαL=0dΦdα2=k2+k12α1+2k22α2+k2PαP+k2LαL=0dΦdαP=kP+k1Pα1+k2Pα2+2kPPαP+kPLαL=0dΦdαL=kL+k1Lα1+k2Lα2+kPLαP+2kLLαL=0 (15)

We may determined the optimal parameter values via (16).

[2k11k12k1Pk1Lk122k22k2Pk2Lk1Pk2P2kPPkPLk1Lk2LkPL2kLL][α1α2αPαL]=[k1k2kPkL] (16)

|M| is presented in (17) and the second order maximum conditions are found in (18).

|M|=|2k11k12k1Pk1Lk122k22k2Pk2Lk1Pk2P2kPPkPLk1Lk2LkPL2kLL| (17)

|2k11|<0,|2k11k12k122k22|>0,|2k11k12k1Pk122k22k2Pk1Pk2P2kPP|<0,|2k11k12k1Pk1Lk122k22k2Pk2Lk1Pk2P2kPPkPLk1Lk2LkPL2kLL|>0 (18)

The optimal parameter values are obtained via (19)–(22).

α1=|k1k12k1Pk1Lk22k22k2Pk2LkPk2P2kPPkPLkLk2LkPL2kLL||M| (19)

α2=|2k11k1k1Pk1Lk12k2k2Pk2Lk1PkP2kPPkPLk1LkLkPL2kLL||M| (20)

αP=|2k11k12k1k1Lk122k22k2k2Lk1Pk2PkPkPLk1Lk2LkL2kLL||M| (21)

αL=|2k11k12k1Pk1k122k22k2Pk2k1Pk2P2kPPkPk1Lk2LkPLkL||M| (22)

Conclusion

It is possible to find optimal solutions also if the management problems have large numbers of integer variables, nonlinearities and stochastic processes. The introduced and tested methods are quite general and can be applied to many other kinds of problems in other sectors. The present approach makes it possible to determine optimal adaptive control rules and to estimate the economic values of mixed forests with trees in many size classes and of many species. With traditional forest management planning methods, the market price variations, locally relevant competition information, multi species management options and variations in timber quality are not considered in the optimal way. It is important to make market adapted harvest decisions. A numerically specified model with empirical data, Lohmander et al (2017), showed that if the stochastic price variations are not considered when the harvest decisions are taken, the expected present value is reduced by 23%. As a result, the economic values of optimally managed forests are underestimated via traditional calculation methods.

Acknowledgments

The author appreciates funding from FORMAS via Linnaeus University in Sweden.

Conflicts of interest

The author declares there is no conflict of interest.

References

  1. Lohmander P. Optimal Stochastic Dynamic Control of Spatially Distributed Interdependent 
    Production Units. In: Cao BY, editor. IWDS 2016: Fuzzy Information and Engineering and Decision. 2018. p. 115–122.
  2. Galantai A. The theory of Newton’s method, Journal of Computational and Applied Mathematics. 2000;124(1–2):25–44.
  3. Wierzbicki AP. A quadratic approximation method based on augmented Lagrangian functions for nonconvex nonlinear programming problems. International Institute for Applied Systems Analysis. 1978. 54 p.
  4. Powel MJD. UOBYQA: unconstrained optimization by quadratic approximation. Mathematical Programming. 2002;92(3):555–582.
  5. Li L, Khandelwal K. An adaptive quadratic approximation for structural and topology optimization. Computers & Structures. 2015;151:130–147.
  6. Lee S, Kwon S, Kim Y. A modified local quadratic approximation algorithm for penalized optimization problems. Computational Statistics and Data Analysis. 2016;94:275–286.
  7. Lohmander P. Adaptive optimization of forest management in a stochastic world. In: Weintraub A, editor. Handbook of Operations Research in Natural Resources. 2007. p. 525–544.
  8. Lohmander P. ICMDS 2016 Conference report. Fuzzy Information and Engineering. 2017;9(2):269–270.
  9. Lohmander P. Applications and Mathematical Modeling in Operations Research. In: Cao BY, editor. IWDS 2016: Fuzzy Information and Engineering and Decision. 2018. p. 46–53.
  10. Lohmander, P. Optimal sequential forestry decisions under risk. Annals of Operations Research. 2000;95(1–4):217–228.
  11. Lohmander P, Olsson JO, Fagerberg N, et al. High resolution adaptive optimization of continuous cover spruce forest management in southern Sweden. Symposium on Systems Analysis in Forest Resources. 2017. p. 37–38.
  12. Lohmander P. A General Dynamic Function for the Basal Area of Individual Trees Derived from a Production Theoretically Motivated Autonomous Differential Equation. Iranian Journal of Management Studies. 2017;10(4):917–928.
  13. Hatami N, Lohmander P, Moayeri MH, et al. A basal area increment model for individual trees in mixed continuous cover forests in Iranian Caspian forests. Journal of Forestry Research. 2018. p 1–8.
Creative Commons Attribution License

©2019 Lohmander. This is an open access article distributed under the terms of the, which permits unrestricted use, distribution, and build upon your work non-commercially.