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

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.¹ 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. Galantai² describes the theory of the Newton’s method. The ideas have been extended in many directions. One such case is Wierzbicki,³ who developed a method for quadratic approximations based on augmented Lagrangian functions for non convex nonlinear programming problems. Another case is Powel,⁴ who created a new algorithm for unconstrained optimization models, leading to quadratic approximations via interpolations. Li⁵ made adaptive quadratic approximations to be used in structural and topology optimization and Lee et al.⁶ 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. Lohmander^7-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?

Lohmander¹⁰ 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. Lohmander⁷ contains several optimization methods and typical solutions to adaptive forest management problems.

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.

$b_i(0)=b_{i0}\forall i$ (1)

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

$d_i(t)=2\sqrt{\left(\frac{b_i(t)}{\pi }\right)}$ (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)\in \left\{0,1\right\}\forall i,t$ (3)

${\displaystyle \sum _{t=0}^{T}u(i,t)}\le 1\forall i$ (4)

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

$s_m(i)=\{\begin{array}{cc}1& ifS(i)=m\\ 0& ifS(i)\ne m\end{array}\forall i,m\in \left\{1,.,N\right\}$ (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,¹ 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.,¹¹ 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.¹² 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.¹³

$b_i(t+1)=\{\begin{array}{cc}{b}_i(t)+G\left(b_i(t),S(i),L(i,t)\right)& foru(i,t)=0\\ 0& foru(i,t)=1\end{array}$ (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={\displaystyle \sum _{t=0}^T{e}^{-rt}{\displaystyle \sum _{i=1}^{I}u(i,t)\left(P(S(i),t)V(S(i),b_i(t))-C(S(i),b_i(t),t))\right)}}$ (7)

$e^{-rt}$ is the discounting factor of period t with rate of interest r.$P(S(i),t)$ denotes price per cubic meter of specie i sin period t.$P(S(i),t)$ t is a stationary variable which is stochastic at $t-n,\forall n>0$ .$E\left(P(S(i),t)\right)$ is the expected value of $P(S(i),t)$ at $t-n,\forall n>0$ . In the two species case, if trees i and j belong to different species, then $P(S(i),t)$ and $P(S(j),t)$ have correlation $\rho$ .$\left(-1\le \rho \le 1\right).V(S(i),b_i(t))$ is the volume of tree as a function of the species and the basal area.$C(S(i),b_i(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”D_i(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.${\alpha }_m$ is the value of the limit diameter D_i(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.

$D_i(t)={\displaystyle \sum _{m=1}^N{\alpha }_m{s}_m(i)}+{\alpha }_P\left(\frac{P(S(i),t)-E\left(P(S(i),t)\right)}{E\left(P(S(i),t)\right)}\right)+{\alpha }_{L}L(t,i)$ (8)

In case the diameter $d_i(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 ${\alpha }_P$ and ${\alpha }_L$ may be different for different species.)

$u(i,t)=\{\begin{array}{cc}0& if\left(d_i(t)-D_i(t)\right)\left(1-{\displaystyle \sum _{t=0}^{t-1}u(i,t)}\right)\le 0\\ 1& if\left(d_i(t)-D_i(t)\right)\left(1-{\displaystyle \sum _{t=0}^{t-1}u(i,t)}\right)>0\end{array}\forall i,t\in \left\{0,\kappa ,2\kappa ,\mathrm{...},n\kappa \right\}$ (9)

$\kappa$ denotes the harvest decision interval. This is an integer,$\kappa \ge 1$ .n is the total number of harvest decision intervals and $n\kappa =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\left({\displaystyle \sum _{t=0}^T{e}^{-rt}{\displaystyle \sum _{i=1}^{I}u(i,t)\left(P(S(i),t)V(S(i),b_i(t))-C(S(i),b_i(t),t)\right)}}\right)$ (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 $\left({\alpha }_{1_0},{\alpha }_{2_0},{\alpha }_{P_0},{\alpha }_{L_0}\right)$ of the optimal values of the parameters $\left({\alpha }_1,{\alpha }_2,{\alpha }_P,{\alpha }_L\right)$ .

Create a number, W, of alternative parameter combinations, w , such that stochastic variables $\left({\epsilon }_1,{\epsilon }_2,{\epsilon }_P,{\epsilon }_L\right)$ are added to $\left({\alpha }_{1_0},{\alpha }_{2_0},{\alpha }_{P_0},{\alpha }_{L_0}\right)$ . The probability density functions of these stochastic variables are defined with consideration of the interesting parameter ranges.$\left(d_{j_{\lambda }},d_{j_{\gamma }}\right)$ denote the lower ($\lambda$ ) and upper ($\gamma$ ) bounds of${\epsilon }_j$, $j\in \left\{1,2,P,L\right\}$ .

$d_{1_{\lambda }}<{\epsilon }_1<d_{1_{\gamma }};d_{2_{\lambda }}<{\epsilon }_2<d_{2_{\gamma }};d_{P_{\lambda }}<{\epsilon }_P<d_{P_{\gamma }};d_{L_{\lambda }}<{\epsilon }_L<d_{L_{\gamma }}$ . Normally, this means that $d_{1_{\lambda }}<0<d_{1_{\gamma }};d_{2_{\lambda }}<0<d_{2_{\gamma }};d_{P_{\lambda }}<0<d_{P_{\gamma }};d_{L_{\lambda }}<0<d_{L_{\gamma }}$ . Let ${\epsilon }_{j_w}$ denote the value of ${\epsilon }_j$ in parameter combination w.$\left({\alpha }_{1_w},{\alpha }_{2_w},{\alpha }_{P_w},{\alpha }_{L_w}\right)=\left({\alpha }_{1_0}+{\epsilon }_{1_w},{\alpha }_{2_0}+{\epsilon }_{2_w},{\alpha }_{P_0}+{\epsilon }_{P_w},{\alpha }_{L_0}+{\epsilon }_{L_w}\right)$ .

For each random parameter combination w, estimate the value of $E(Z)=E\left(Z\left({\alpha }_{1_w},{\alpha }_{2_w},{\alpha }_{P_w},{\alpha }_{L_w},\mathrm{...}\right)\right)$ .

Make a quadratic approximation $\Phi =\Phi \left({\alpha }_1,{\alpha }_2,{\alpha }_P,{\alpha }_L\right)$ of the function $E(Z)=E\left(Z\left({\alpha }_1,{\alpha }_2,{\alpha }_P,{\alpha }_L,\mathrm{...}\right)\right)$ 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.

$\Phi =\Phi \left({\alpha }_1,{\alpha }_2,{\alpha }_P,{\alpha }_L\right)$ (11)

$\begin{array}{l}\Phi =\Phi \left({\alpha }_1,{\alpha }_2,{\alpha }_P,{\alpha }_L\right)\\ =k_1{\alpha }_1+k_2{\alpha }_2+k_P{\alpha }_P+k_L{\alpha }_L+\\ +k_{11}{\alpha }_1{}^2+k_{22}{\alpha }_2{}^2+k_{PP}{\alpha }_P{}^2+k_{LL}{\alpha }_L{}^2+\\ +k_{12}{\alpha }_1{\alpha }_2+k_{1P}{\alpha }_1{\alpha }_P+k_{1L}{\alpha }_1{\alpha }_L+k_{2P}{\alpha }_2{\alpha }_P+k_{2L}{\alpha }_2{\alpha }_L+k_{PL}{\alpha }_P{\alpha }_L\end{array}$ (12)

$\Phi ({\alpha }_1,{\alpha }_2,{\alpha }_P,{\alpha }_L)\approx E(Z({\alpha }_1,{\alpha }_2,{\alpha }_P,{\alpha }_L,\mathrm{...}))$ (13)

$\underset{k_0,k_1,\mathrm{...},k_{PL}}{\min }\Psi =E\left[{\left(\Phi \left({\alpha }_1,{\alpha }_2,{\alpha }_P,{\alpha }_L;k_1,\mathrm{...},k_{PL}\right)-E(Z({\alpha }_1,{\alpha }_2,{\alpha }_P,{\alpha }_L,\mathrm{...}))\right)}^2\right]$ (14)

Use the quadratic approximation $\Phi =\Phi \left({\alpha }_1,{\alpha }_2,{\alpha }_P,{\alpha }_L\right)$ to determine the approximately optimal values of the parameters ${\alpha }_1,{\alpha }_2,{\alpha }_P,{\alpha }_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).

$\{\begin{array}{c}\frac{d\Phi }{d{\alpha }_1}=k_1+2k_{11}{\alpha }_1+k_{12}{\alpha }_2+k_{1P}{\alpha }_P+k_{1L}{\alpha }_L=0\\ \frac{d\Phi }{d{\alpha }_2}=k_2+k_{12}{\alpha }_1+2k_{22}{\alpha }_2+k_{2P}{\alpha }_P+k_{2L}{\alpha }_L=0\\ \frac{d\Phi }{d{\alpha }_P}=k_P+k_{1P}{\alpha }_1+k_{2P}{\alpha }_2+2k_{PP}{\alpha }_P+k_{PL}{\alpha }_L=0\\ \frac{d\Phi }{d{\alpha }_L}=k_L+k_{1L}{\alpha }_1+k_{2L}{\alpha }_2+k_{PL}{\alpha }_P+2k_{LL}{\alpha }_L=0\end{array}$ (15)

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

$\left[\begin{array}{cccc}2k_{11}& k_{12}& k_{1P}& k_{1L}\\ k_{12}& 2k_{22}& k_{2P}& k_{2L}\\ k_{1P}& k_{2P}& 2k_{PP}& k_{PL}\\ k_{1L}& k_{2L}& k_{PL}& 2k_{LL}\end{array}\right]\left[\begin{array}{c}{\alpha }_1\\ {\alpha }_2\\ {\alpha }_P\\ {\alpha }_L\end{array}\right]=\left[\begin{array}{c}-k_1\\ -k_2\\ -k_P\\ -k_L\end{array}\right]$ (16)

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

$\left|M\right|=\left|\begin{array}{cccc}2k_{11}& k_{12}& k_{1P}& k_{1L}\\ k_{12}& 2k_{22}& k_{2P}& k_{2L}\\ k_{1P}& k_{2P}& 2k_{PP}& k_{PL}\\ k_{1L}& k_{2L}& k_{PL}& 2k_{LL}\end{array}\right|$ (17)

$\begin{array}{l}\left|2k_{11}\right|<0,\left|\begin{array}{cc}2k_{11}& k_{12}\\ k_{12}& 2k_{22}\end{array}\right|>0,\\ \left|\begin{array}{ccc}2k_{11}& k_{12}& k_{1P}\\ k_{12}& 2k_{22}& k_{2P}\\ k_{1P}& k_{2P}& 2k_{PP}\end{array}\right|<0,\left|\begin{array}{cccc}2k_{11}& k_{12}& k_{1P}& k_{1L}\\ k_{12}& 2k_{22}& k_{2P}& k_{2L}\\ k_{1P}& k_{2P}& 2k_{PP}& k_{PL}\\ k_{1L}& k_{2L}& k_{PL}& 2k_{LL}\end{array}\right|>0\end{array}$ (18)

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

${\alpha }_1=\frac{\left|\begin{array}{cccc}-k_1& k_{12}& k_{1P}& k_{1L}\\ -k_2& 2k_{22}& k_{2P}& k_{2L}\\ -k_P& k_{2P}& 2k_{PP}& k_{PL}\\ -k_L& k_{2L}& k_{PL}& 2k_{LL}\end{array}\right|}{\left|M\right|}$ (19)

${\alpha }_2=\frac{\left|\begin{array}{cccc}2k_{11}& -k_1& k_{1P}& k_{1L}\\ k_{12}& -k_2& k_{2P}& k_{2L}\\ k_{1P}& -k_P& 2k_{PP}& k_{PL}\\ k_{1L}& -k_L& k_{PL}& 2k_{LL}\end{array}\right|}{\left|M\right|}$ (20)

${\alpha }_P=\frac{\left|\begin{array}{cccc}2k_{11}& k_{12}& -k_1& k_{1L}\\ k_{12}& 2k_{22}& -k_2& k_{2L}\\ k_{1P}& k_{2P}& -k_P& k_{PL}\\ k_{1L}& k_{2L}& -k_L& 2k_{LL}\end{array}\right|}{\left|M\right|}$ (21)

${\alpha }_L=\frac{\left|\begin{array}{cccc}2k_{11}& k_{12}& k_{1P}& -k_1\\ k_{12}& 2k_{22}& k_{2P}& -k_2\\ k_{1P}& k_{2P}& 2k_{PP}& -k_P\\ k_{1L}& k_{2L}& k_{PL}& -k_L\end{array}\right|}{\left|M\right|}$ (22)

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.

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

The author declares there is no conflict of interest.

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.