An optimal control problem in forest management

An optimal control problem in forest management G. Corriga, D. Salimbeni, S. Sanna, and G. Usai lstituto di Elettrotecnica, Universit~ di Cagliari Piazza d'Armi, 09100 Cagliari, Italy (Received April 1987; revised August 1987)

An optimal time control problem arising in forest stand management is presented. The paper refers to a nonlinear dynamic model, expressed in terms of diameter classes, of the growth of a Mediterranean forest stand. The problem is tackled by resorting to the Pontryagin maximum principle and is thus reduced to a two-point boundary value problem. Being the Hamiltonian linear in the control, a singularity analysis is carried out which, taking into account the model structure, excludes the existence of singular extremals. A numerical example is finally presented.

Keywords: forest growth control, minimum time control Introduction The main goal pursued in the majority of forest management problems is the optimization of the revenue of the timber resources. ~-3 Numerous constraints, often ecological in nature, may limit the choice of optimum management strategy by prohibiting, for instance, the resort to certain harvesting techniques, sometimes very simple but which violate environmental constraints, such as clear-cut harvesting of mature trees and subsequent sowing or planting of new ones. Less frequent are those problems related to management of forest systems where the primary objective is to modify the forest structure, expressed in terms of relative incidence of the different tree diameter classes, so as to obtain in the shortest possible time a different structure more suited than the original one to the intended purpose of the forest resources. This is the case of artificial forest stands originally planted for some specific purpose, such as checking the advance of coastal dunes gradually invading agricultural land. In such situations, once the original goal has been attained, the forest structure is often ill-suited to commercial and/or tourist exploitation and must first be altered in order to ensure a desired future state, determined on the basis of forestry considerations, which will eventually allow a rational, even economic, exploitation. Forest stands of this type are usually particularly thickly grown with high density per unit area and, despite usually being even-aged, great diversity in tree size. Basically, optimum cutting rules need to be de-

328 Appl. Math. Modelling, 1988, Vol. 12, June

termined for the different diameter classes which achieve the desired stand structure in the least possible time. A prerequisite for addressing any optimization problem is the availability of appropriate mathematical models of forest growth. The paper presents a nonlinear model, expressed in terms of diameter classes, originally proposed by R. Leary 4 and adapted by the authors to suitably describe the dynamics of Mediterranean forest systems? The problem is tackled by resorting to the maximum principle of Pontryagin and hence is reduced to a twopoint boundary value problem. It is shown that for the model proposed the optimum cutting rule is always of the bang-bang type, where the determination of the extremal trajectories is particularly crucial. In the specific case, the method used 6 allows us to determine with a high degree of accuracy the switching times, which is a fundamental step for finding the optimum solution. The algorithm, implemented on an HP personal computer, has been applied to a three-diameter-class model.

Mathematical model of the forest stand The models examined herein have already been presented by the authors elsewhere. 5 They were originally proposed for uneven-aged forest stands. 4 Nevertheless, they can be applied to even-aged stands where competition between the components gives rise to different class diameters within the same stand. Each

© 1988ButterworthPublishers

Control problem in forest management: G. Corriga, D. Salimbeni, S. Sanna, G. Usai diameter class is characterized by a nominal diameter d~ which expresses the mean of the lower- and upperdiameter bound of the class. The different diameter classes have been identified in a certain initial instant to. In Refs. 4 and 5 two structurally similar models have been considered in which the state variables represent the sum of the diameters at breast height, referred to a unit area, of trees belonging to different diameter classes: Ni

xi(t) = x~ dij(t) j=

i = 1,2 . . . . . n

t > t,,

(I)

I

where d o. is the diameter of the jth tree in the ith diameter class, N; is the number of trees belonging to the diameter class d~, with d~ > d i ~, and n is the number of diameter classes in the stand. The first model, 2i(t) = aixi(t) exp

[1

- bi~xi(t) j=

i

"

- ui(t) _J

i = 1,2 . . . . . n

(2)

where th and b~ are constant parameters, greater than zero, which depend on the climatic-pedologic features and on the degree of interaction between the different diameter classes, is particularly appropriate for describing the growth dynamics of a forest system where the dominating effect of the taller over the smaller trees is prevalent. This is true, for instance, where direct solar radiation is the prime growth factor. The structure of the model makes explicit allowance for the fact that trees in each diameter class must only compete with those belonging to larger-diameter classes. The input uM), which specifically represents here the sum of tree diameters, in the ith class, cut in unit time, was explicitly introduced to enable identification of optimum cutting rules. In the case of bang-bang control, the problem consists in determining those periods, usually different for each diameter class, when the trees should not be cut (u~(t) = 0) and those when cutting should be done, cutting rate being known and constant (ui(t) = UiM). The second model, 2i(t) = a~x~(t) e x p

-b~

x~(t)

-uM)

i = 1,2 . . . . . n

(3)

whose parameters have the same meaning, expresses the fact that, similarly to the first model, growth rate of trees belonging to the ith diameter class is directly proportional to the number of trees in that class. However, it differs from the first model in that here the growth rate decreases exponentially, for any diameter class, in proportion to the total number of trees in the stand. This model is more appropriate for forest systems where the assumption that the trees belonging to the large-diameter classes do not have to compete with

those of the smaller ones, is unacceptable. This is true of certain types of forests in Mediterranean countries where there is an abundance of solar radiation, direct or indirect, throughout the vegetative phase (spring, summer, part of autumn). In this case it is not so much the competition for solar radiation, but rather that for water and nutrients in the soil by the roots of the trees that is the prime growth factor. Obviously, it is assumed that in this instance all the diameter classes compete with each other. This situation is represented by model (3). Since the study was concerned with forest dynamics in a Mediterranean environment, the mathematical model expressed by (3) is referred to in the sequel in the formulation of the optimum control problem, i.e., that of finding the cutting rule whereby the stand structure can be modified from an initial known state to a desired future state in the shortest possible time, as well as in the numerical example. Model (3) differs substantially from model (2) as to the properties of the solutions to the minimum time control problem. As shown in the fourth section, the model (3) structure excludes the existence of singular extremal controls, and hence the solution of the control problem is of the bang-bang type. Conversely, the model (2) structure allows for the existence of singular extremals, which, besides not having any practical interest, complicates the search for the optimum solution. Obviously, this does not exclude that under different environmental conditions (e.g., in north European climates) model (2) m a y b e more fitting. The coefficients ai and bi appearing in the two models can be calculated 5 once the growth data of the trees of the different diameter classes are known. In particular, yearly measurements of the diameters of a representative sample of the whole stand are required. Should such data not be available, they can often be determined a posteriori by cutting a sample portion of the stand and measuring the yearly growth rings (or using other growth measuring techniques applied in forestry). Obviously the validity of the models depends on the reliability of the growth data used for checking the models themselves as well as on the period covered by the data. Since the coefficients ai and b~ are determined on the basis of measurements taken over a long time, the models are not capable of representing seasonal growth but are more suited for describing stand growth over wider time horizons (several years).

Formulation of the optimum problem Having applied the maximum principle of Pontryagin, we can reduce the minimum time control problem-since, as shown in the next paragraph, the solution is of the bang-bang type and that no singular extremal control exists--to a two-point boundary value problem of the f o r m ( i = 1,2 . . . . . n) OH

~i = - -

c~Ai

,~i -

OH

OXi

Appl. Math. M o d e l l i n g , 1988, Vol. 12, J u n e

(4)

329

Control problem in forest management: G. Corriga, D. Salimbeni, S. Sanna, G. Usai xi(to) = Xio

Ai(to)= Bi

(5)

H(x(ts), A(ts), u(ts)) = 0

(6)

to <-t<-ts xi(ts) = xsi u~(t) = uiM

ifA~> 0

uAt) = 0

ifAi< 0

(7)

where the x / s are the components of state vector x the A,.'s are the components of the costate vector h. H = H(x(t), A(t), u(t)) is the Pontryagin state function (or Hamiltonian) X(to) denotes the known initial state, x(ts) the desired final state the B;'s, the components of the vector B, are the initial, unknown values of the costate variables A~ the relation H(t s) = 0, evaluated along an extremal trajectory, must be introduced as the (n + l)st condition necessary, along with the n boundary conditions referred to the final time, to find the n + 1 unknowns B~ and t S As is well known, the problem of finding the minimum time and the related control laws can be solved with iterative procedures, which consist of initially setting tentative values for t s and initial values Bi of the costate variable, then integrating the system (4) with the control laws (7) and verifying whether the final conditions (6) are fulfdled; i.e., compactly, we can write

..... x!l:? x: ..... ] L(B,ts) =

H(x(,s),A(,s),u(ts,,j = [':"]

(8)

L is a vector function with n + 1 components. Because equation (8) cannot usually be initially verified, we must resort to an appropriate algorithm which updates, at each iteration, the values of ts and B~. In the case in question, the algorithm due to Lastman 6 is used, where the corrections to the t~k' and B ~k~obtained from the kth iteration are expressed by

[

control problems in that an error in determining any one of them may lead to an error of the same or greater magnitude in the final trajectories, with the result that the final boundary conditions (8) cannot be met with the desired degree o f accuracy• In the application exemplified here, an algorithm based on the Newton-Raphson approach is used for finding the zero of a function.

Singularity analysis In solving optimal control problems where the Hamiltonian is linearly related to the control and the performance index does not contain the control, as in the case at hand, it typically happens that one or more components of the costate vector are identically zero in one or more finite time intervals. During such intervals the function H no longer explicitly depends on the corresponding component of the control vector, and as a result it does not provide any information on control extremals. In this event, the optimal control is bang-bang in certain time intervals, whereas in others it is expressed by appropriate time functions, called singular controls. 7"8 The aforementioned procedure, like other algorithms constructed for solving optimal control problems, usually fail where a singularity exists. Thus we want to know a priori whether singular solutions exist. Their existence depends on the structure of the model used to represent the dynamics of the system being controlled• In the case concerned, referring to model (3), we want to show that the problem does not admit extremal singular controls and that the optimum solution is therefore of the bang-bang type. To do this we only need to show that at least one of the necessary conditions for the existence of singular solutions is not satisfied: Setting

AB"'-] aixiexp

Ate-" J i

FOL(B''', t~") aK L O-B

• L(B"', t~") = - a t

OL(B"', t~")] -I Ots J

F[ 6t~" B"'-IJ

- bi

= fi

i = 1,2 . . . . . n

we can write the Hamiltonian function as (9)

and the scalar a r (0 < a r -< 1) is chosen such that

n

H = 1 + ~ a;(f; - ui)

(10)

i=l

IIL[B'" - ,~K,~B'", t~k' - ~KSt~'k']ll < IIE[B'", t~"] II which represents the convergence condition. The term OL/0B in (9) represents an (n + 1) x n matrix where the (i,j)th element is OLi/OBfi OL/Ots is instead a column matrix with OL~/Ots as the ith component. For the n switching functions, which in this case are the costate variables A~themselves, the problem arises of finding with very high accuracy the switching times. These parameters are critical in bang-bang optimal

330 Appl. Math. Modelling, 1988, Vol. 12, June

and it is identically zero along any extremai trajectory, since the control problem concerned has free final time. As stated in the foregoing, the singularity condition is verified if at least one costate variable is identically zero in a finite time interval. Another necessary condition, shown by B. S. G o h , 9 requires, in our case, that the Hessian matrix a2H] OXiOXj_]

i=1,2 ..... n j=l,2

..... n

Control problem in forest management: G. Corriga, D. Salimbeni, S. Sanna, G. Usai be positive semidefinite. This implies that all 2 x 2 determinants

I O2H

OXnOXs Ox,,OX~ Ox~

oq-I

s = 1,2 . . . . . n - 1

o2__H_H axe,

be greater than or equal to zero; i.e.,

Of course, such relations can only have any significance in the case of equality. It follows that if we assume at least one of the costate variables A; (i = 1, 2..... n) is zero, all terms of the type b~fi/x~ (i = 1, 2 . . . . . n) being different from zero, then all the other costate variables should also be zero. But this can never happen because we would get from (10) that H = 1, but the Hamiltonian function must be identically zero. This proves that no singular solutions exist for model (3).

The considerations expounded in the foregoing were applied to a problem of forest management using a mathematical model of the type expressed by (3), consisting of three diameter classes. The forest in question, a Pinus Pinea stand, was identified in an earlier work, s which used quasi-linearization techniques. The identification of a sample plot of 400 m 2 representative of the entire stand led to the following parameter matrix:

[a, b,]=/1.350 I-0.504-0.0.00457 003271 / a2

a3

b2 b3

L0.791

0.00313_1

The growth data used in identifying the model furnished the following actual values of the three diameter classes (expressing the sum of diameters of each class in centimeters): x(to) = [670

205

36] t

which gives the initial state. The final state, chosen on the basis of silvicultural considerations, can be expressed by the vector x(t/) = [540

180

1. Since the growth rates of the different diameter classes increase, according to model (3), with decreasing total biomass present in the stand, the minimum time control strategy leads to carrying out the required cutting at the initial instant. 2. For the same reasons, connected with the interaction between the different diameter classes, high cutting rates, which would involve a more rapid thinning out of the forest, allow the final desired state to be achieved in less time. In the case examined we have tf = 5 years and 5 months for u;M = 1200 cm/yr t s = 6 years and 4 months for u~M = 100 cm/yr

Numerical example

A=

In the case at hand uiM between 100 and 1200 cm/yr were considered. The above values as well as the values X(to) and x(tz) refer obviously to the area of 400 m 2. In the case of a forest of area S m 2 all these values should be multiplied by the ratio S/400. The solution to the problem is summarized in the diagrams of Figure 1, which show the curves for the two limit values taken for uiM. The values of u,-Mwithin the above range correspond to intermediate curves between the two limit curves for each diameter class. Some interesting qualitative considerations can be drawn from the results obtained:

60] T

The maximum value of the control uiM (i = 1, 2, 3) has been taken the same for all diameter classes. However, this does not impose any constraint because, if desired, different values can be assumed for the different classes. In reality, this is an indication of the " r e s o u r c e s " in terms of labour and machinery used in cutting and derives from precise technoeconomical considerations in relation to the different results obtainable with different cutting rates.

3. As a limit situation we can consider the one where cutting rate is infinite (impulsive control). This type of control (purely theoretical, however) would instantaneously transfer the initial known state to a new one, starting from which all the diameter classes would develop freely. This would result in achieving the desired final state in the absolute minimum time in the case at hand. 4. Finally, it is clear that should the cutting rate be too low, so that a nondecreasing pattern arises during the cutting period, even in just one of the classes to be thinned out, no solution to the problem of minimum time control exists. 700

Gm

5~

4u

~. 3m

zm

x 2 (t)

IN

x 3

....

,.,,,,,1

........ I

i,.I

......

iir..I

i~

..... 3

TIME

(t)

ii1.,,I,,

......... 4

I,, $

.........

I .....

,

6

[year:]

Figure 1 Optimal state trajectories in the computed solutions to the minimum time control problem with UiM = 100 (solid lines) and UiM = 1200 (dashed lines)

Appl. Math. Modelling, 1988, Vol. 12, June 331

Control problem in forest management: G. Corriga, D. Salimbeni, S. Sanna, G. Usai

Concluding remarks

References

A mathematical model is presented that represents the growth of a monospecies forest stand, expressed in terms of class diameters. The model allows for the specific interaction between trees of the different classes and is particularly suited for representing the dynamics of a Mediterranean forest stand. Model parameter identification was carried out by the authors, 5 resorting to a quasi-linearization technique and using growth data records. The problem of minimum time control, formulated for the model, allows us to solve forest management problems that require the alteration of the forest stand structure, expressed in terms of relative incidence of the different diameter classes, from an initial known state to a desired final state. It has been shown that for the model used the optimum cutting law is always of the bang-bang type, excluding the existence of singular extremal controls.

I

332 Appl. Math. Modelling, 1988, Vol. 12, June

2 3 4 5

6 7 8 9

Proceedings Workshop on Computer and Information Systems in Resource Management Decisions. U.S.D.A., Cleveland, Ohio. 1971 Field, D. B. Goal programming for forest management. Forest Science 1973. 19, 125-135 Adams, D. M. and Ek, A. R. Optimizing the management of uneven-aged forest stands. Can. J. For. Res. 1974, 4, 274-287 Leary. R. A. System identification principles in studies of forest dynamics. U.S.D.A. Forest Service, North Central For. Expt. Sin., St. Paul, Minn.. Res. Pap. NC-45 1970 Corriga, G., Salimbeni, D.. Sanna. S.. and Usai, G. Identification and validation of a forest dynamic model. Thirteenth lASTED International ConJ~,reme on Modelling and Simulation, Lugano. 1985 Lastman, G. J. A shooting method for solving two-point boundary-value problems arising from non-singular bang-bang optimal control problems. Int. J. Control 1978, 27, 513-524 Athans, M. and Falb, P. L. Optimal Control. McGraw-Hill, 1966 Bell, D. J., and Jacobson, D. H. Singular Optimal Control Problems. Academic Press, 1975 Goh, B. S. Necessary conditions for singular extremals involving multiple control variables. J. Siam Control 1966, 4, 716-731