A mathematical model for agroforestry optimization
Abstract: In the present article, we will describe some extensions of an agroforestry model that has been proposed
and computationally implemented in [7]. Our generalizations consist of the inclusion of two additional species
of tree, one culture, and a declaration of regeneration tours as variables definable by us as a parameter and the
weight allocation by rentability of the treeless soil utilization as well as an exhaustive exploration of the different
soil utilization scenarios in order to obtain the one who gives the best economic performance.
KeyWords: Optimal control, Mathematical modeling, Agroforestry, Environmental sustainability
Received: April 8, 2021. Revised: January 20, 2022. Accepted: February 10, 2022. Published: March 2, 2022.
1 Introduction
The problem of forest regeneration in lands whose
land use is shared between forest preservation, agri
culture, or cattle grazing has been one of the great
challenges of forest preservation. Such a situation
is the case of the Dehesas in Spain where private
landowners who make mixed use of the land are faced
with the dilemma of making sustainable use of their
lands in environmental terms and, on the other hand,
exploiting whatever resource is economically attrac
tive. In this last sense, various initiatives can be ob
served in Europe to promote the forestry development
of these multipleuse lands through government sub
sidies to stimulate the planting of trees. The develop
ment of wooded areas in these lands in private hands
can turn them into spaces for ecotourism or even de
pending on the type of tree planted; it could be used
to produce cork such as cork oaks or produce acorns
to feed pigs whose meat is highly valued for its qual
ity as would be the case of holm oaks. In the review
that we have carried out on this topic, we can men
tion the work carried out in the case of Dehesas in
Spain by the authors of [4, 5, 7]. In this article, we
are inspired by the work carried out by [7] where the
problem of economic and sustainable optimization of
pastures in the Andalusian region in southern Spain
is studied. In this article, he proposes a model of eco
nomic and agroforestry optimization where land use
for grazing, agriculture, and sowing of oaks and cork
oaks is combined. In this mixed use of the land that
they carry out in the article mentioned above, they
define some rules of rotation of both pasture crops
and tree renewal based on the lifetime of each tree
species and also taking into account the natural refor
estation specific to each species. The distribution of
wooded and nonwooded areas from which the opti
mization model is started is introduced as initial con
ditions. Once some of these parameters have been
defined, a linear optimization is carried out, taking
into account the restrictions imposed by the needs of
crop rotation as well as the search for obtaining eco
nomic profit. Of course, one of these restrictions is
not exceeding the limit of land available. We have
included the details of the model exposed in [7] as
well as a description of it in the section 9 of this ar
ticle for the purpose of completeness of the exposi
tion of our model. Based on the initial distribution of
surfaces as well as nonnegativity restrictions of the
variables involved, a MATLAB linprog toolbox for
linear optimization routine is called that produces the
decisions that will be made at each time step, and this
information feeds the dynamic model that describes
the behavior over time of the system from the appli
cation of the decisions calculated by the linear opti
mizer mentioned above. The time window used in
both the [7] model and ours is 0to 100 years. The so
lution to the forest optimization problem raised in [7]
is approached as an optimal control problem in dis
crete time and finite time horizon and that describes
the agroforestry decision processes by social and pri
vate agents in a pasture system. In this dynamic op
timization model, laws of motion of the states are de
fined during the time horizon, control variables and
conditions of nonnegativity of these, global restric
tion of resources, as well as artificial and natural re
generation restrictions of adult trees that are useful for
1CARLOS RODRÍGUEZ LUCATERO, 2MARCELO OLIVERA VILLAROEL,
3PAOLA OVANDO
1Departamento de Tecnologías de la Información
Universidad Autónoma Metropolitana Unidad Cuajimalpa, MEXICO
2Departamento de Teoría y Procesos del Diseño
Universidad Autónoma Metropolitana Unidad Cuajimalpa, MEXICO
3Research Staff Member Spanish Council for Scientific Research, SPAIN
WSEAS TRANSACTIONS on SYSTEMS and CONTROL
DOI: 10.37394/23203.2022.17.13
Carlos Rodríguez Lucatero,
Marcelo Olivera Villaroel, Paola Ovando
E-ISSN: 2224-2856
108
Volume 17, 2022
the definition of physical and legal restrictions in each
case study addressed in our research. Some additional
restrictions relative to the planting of holm oaks and
cork oaks on bare ground are added for each specific
case, which guides forest management to a steady
state. Similarly, an objective function is defined that
represents the variable of income at net present value
at a fixed horizon for each of the Dehesa farms stud
ied in [7]. The same authors of [7] define a gradi
ent vector of the objective function that allows them
to understand the operation of the proposed dynamic
optimization program, and from the expressions that
appear in said gradient vector, the linear nature of the
decision problem can be observed. In the model that
we have just described very quickly in the previous
paragraphs, a very complete analysis of the income
or earnings of a farm owner is made both in scenar
ios where there is public intervention through subsi
dies and incentives to friendly developments with the
environment, as in scenarios where there is no pub
lic intervention. In that sense, this work is very com
plete. However, the rules that are defined for the use
of the land for reforestation purposes or for agricul
ture are carried out by defining very specific rules
to the particular needs of each of the farms studied,
which seems to us to be a limitation of this model.
Another limitation of the model mentioned above is
that it only works with two tree species (holm oak and
cork oak). On the other hand, the optimization is car
ried out from a single scenario of initial distribution
of crops in bare soil. Due to these limitations in the
model proposed in [7], we propose a generalization
of the previous model by adding a greater number of
tree species, discarding the Adhoc rules of crop dis
tribution in bare soil and evaluating all possible crop
distributions on bare ground to choose the scenario
that maximizes income. The latter can generate de
generate cases in which the global resource restriction
is violated and that we eliminate after processing. On
the other hand, this exhaustive test of the possible ini
tial distributions can be costly in terms of calculation
time, which we try to alleviate by means of a random
sampling of the space of possible initial distributions
of bare soil. This model expands and refines [5] by
allowing a more flexible framework for land use re
placeability and tree species selection. This approxi
mation then allows us to examine in discrete time the
composition and distribution of managed forest, aban
doned forest, and nonforest uses at the farm level,
which characterizes private land holdings in the re
gion that is the empirical focus of our paper. This
type of models has been previously applied to ana
lyze efficient forest investment decisions through af
forestation programs [22] forest management invest
ments [5] and forest management decisions, such as
thinning [20],[21]. Optimal control techniques in dis
crete time have also been applied to solve farmland
decision problems, whereas the farm owner or holder
has a preference to conserve forest areas within the
farm [3],[19]. The contribution of our paper is as fol
lows. Our primary focus contributes to understand
ing how different initial land use distributions and en
vironmental conditions affect private actors’ invest
ment decisions in forest management and afforesta
tion. These changes in the distribution, in turn, throws
light on how these decisions affect the longerterm
provision of ecosystem services and (relatedly) for
est ecosystem asset values over time [24],[25]. That
is what our contribution in this article consists of,
which will be described in the following sections of
this work.
2 Optimal management of forest
systems based on environmental
economic behavior
In this article we develop a framework of economic
analysis for agroforestry systems by implementing of
a bioeconomic mathematical model based on a dy
namic optimization system inspired in the work de
veloped by [7] and [22]. For this reason we use as
object of evaluation of the optimization system, the
rent at present value in perpetuity of a typical farm.
The model allows to see the impact on the conser
vation of forest areas through the implementation of
policies of subsidies to forest products, the payment
of environmental services, the subsidies of the gen
eration of commercial plantations and the effects of
livestock control without restrictions in the systems
Forestry [8], being able to extrapolate to different
types of ecosystems as long as the level of analysis
is a farm type. The calculation of the annual income
of the farm type at present value considers an infinite
horizon time whereby it develops a process of nor
malization for generating an average farm. The op
timization model uses as methodology the resolution
of an optimal control problem in discrete time follow
ing the research line of [1],[3], [7], [19] considering
the decision of the agents that manage a typical farm
with forestry systems within their production system
and the preference to conserve forested area on the
farm. It is considered that the income of the farm in
cludes a function of the price of the land determined
by the age of the trees and the agricultural activi
ties developed on the farm, being a substitute of the
price estimated by the market of forest systems. It is
assumed that decisionmakers are maximizers of the
economic rent [4], [7]. In order to obtain the optimal
income of the farm analyzed, different land use op
tions are categorized in order to maximize the income
at present value of the integrated management of the
farm. In the model proposed in [7] it is assumed that
WSEAS TRANSACTIONS on SYSTEMS and CONTROL
DOI: 10.37394/23203.2022.17.13
Carlos Rodríguez Lucatero,
Marcelo Olivera Villaroel, Paola Ovando
E-ISSN: 2224-2856
109
the owner decides between the different land use op
tions from the initial distribution that faces as owner.
That is, it decides which should be the percentage of
land allocated to different forest and agricultural op
tions. Such decisions are taken in concordance with
the costs and longterm benefits that offer the farm
size. One of the possible options is the abandonment
of the farm in the long term. The model considers
that the different uses of the land are replaceable be
tween them, as long as the cost benefit of the chosen
activity pays the replenishment and opportunity costs
of the activities to be replaced, considering the sub
stitution between more than two forest species, thus
expanding the base model developed by [7] That is,
the objective function of the analysis framework is a
decision model where the renewal of the forest sys
tem associated with the farm, that has as options the
natural renewal with associated management, the ar
tificial renewal with the development of plantations
or the abandonment of the systems forestry with the
introduction of livestock systems without associated
forest management. Unlike [7] the model considers
nonlinear cost systems in the development of the in
come models of the type farms analyzed, consider
ing the substitution of more than two forest species
in the farm as well as pastures and agricultural sys
tems. Moreover, in contrast to the works of [11] and
[22] the farms abandonment is considered as a plau
sible option to the agroforestry management model.
Within the forest management process of the analysis
framework is considered a decision already made in
advance, optimal shifts according to the forest species
existing in each farm, that it is considered an exoge
nous decision to the process of optimal forest control,
solving the problem of [10]. The model’s objective
is to choose the type of renewal of forest systems as
sociated with the farm through the selection of opti
mal sequences of reforestation, forest plantations, and
natural regeneration of existing adult trees, abandon
ment of the forest or the implementation of new agri
cultural areas of pastures. The final goal is to offer
an optimal distribution of land use within the farm,
each time one of the cutting shifts of the existing for
est species is effective. The model takes into account
a typical farm considering the following restrictions
and assumptions to obtain the annual income of the
enterprise according to the type of renewal of the for
est systems associated with the farm:
1. The number of commercial forest species on the
farm
2. Production cycles of each forest species. The
oretical growth rate of the commercial forest
species
3. The land use of the farm including conserva
tion areas, agricultural (cereal crops, fruit trees,
grassland among others), forest management and
abandoned or dismantled areas
4. Age of the existing trees on the farm
5. Depending on the case, the management of ma
ture forests (steady state) or forest plantations in
their initial stages of growth can be considered.
6. Types of management of forest systems, if it is a
natural regeneration or a commercial plantation
7. Profits from agricultural activities are considered
part of the land use income
8. Livestock activities directly affect the processes
of natural regeneration
9. The costs and benefits of the production of
forestry and agricultural species in the long term
are considered, as well as the provision of other
types of services associated with forest systems
and the farm as recreational and hunting activi
ties.
10. Forestry activities include the collection of non
timber and timber products
11. Subsidies are considered by productive systems,
that is, forestry, agriculture and/or conservation
12. The costs and sale prices within the production
system are normalized to a reference year to
avoid short and medium term distortions such as
inflationary shocks or the exchange market. That
is, it is considered a constant price for the devel
opment of the analysis.
13. The model considers both positive and negative
price shocks
14. The model considers that the farm is taking mar
ket prices and cannot affect them in the short or
long term
15. Different discount rates are considered
In our model we call this initial distribution of land
use options scenario. Our model is a generalization of
the model proposed in [7] because we generate each
possible scenario or initial land distribution options
and for each possible scenario we apply the same kind
of decision optimization performed in [7]. We obtain
and record for each scenario its economic utility. Af
ter processing all the possible scenarios we select the
one that has the higher economic utility.
WSEAS TRANSACTIONS on SYSTEMS and CONTROL
DOI: 10.37394/23203.2022.17.13
Carlos Rodríguez Lucatero,
Marcelo Olivera Villaroel, Paola Ovando
E-ISSN: 2224-2856
110
3 Current dynamics of European
forest systems
The current dynamics of European forest systems
cannot be understood without the interference of hu
man activities in the forest environment; since the
composition and structure of the landscape is condi
tioned to the interactions of the existing conditions in
space and the use of space by humans [23]. The use
of European forest systems traditionally seeks to ef
fectively use the resources available in the environ
ment. Thus, the current distribution conditions of
agricultural activities with the efficient use of pastures
and land suitable for cultivation are associated with
the composition of forest systems for the extraction
of timber and nontimber resources; [12] generating
cultural landscapes with its own dynamics of biodi
versity and economic productivity [6] An example of
this, the cork oak forests and the combination of ce
real crops and pastures and hunting activities in the
region. The composition and structure of forests are
strongly associated with the profitability of the eco
nomic and of each region’s economic and social dy
namics so that the landscape units of these systems
are always in constant movement. This is evident in
the last 50 years with the demographic changes of ru
ral systems in Europe, where there is an abandonment
of rural activities and there is an ageing of population
structures, an example is the abandonment of farm
land, or the raising of sheep and goats [5], [15]. The
change in demographic composition has had an im
pact on the amount of energy contributed by humans
to their environment, leaving aside the agricultural
and forestry exploitation and with it a change in the
structure of the landscape [23]. Changes in this artifi
cial composition of the environment have many con
sequences among the most relevant: a loss of value
in the economic production of forest systems for their
owners and therefore an abandonment of preexisting
forest management; loss and spaces with very im
portant cultural, economic and social values; losses
of crops and pastures due to the increase in transi
tional spaces of thickets and groves, although this has
a lower erosion rate, higher water quality in rivers and
reservoirs and an unclear result on the gain and loss
of biodiversity according to the stage of forest suc
cession in which the abandoned space is located [17],
[18]. The longterm effect over outgoing is the gain
inhomogeneity of forest structure and a lower risk of
fires, although in the process of forest succession can
be generated colonization of shrub species, generalist
arboreal, and pyrophytes that make the landscape lose
its economic and environmental value and increase
their risk of fires [2] In this sense, the European Union
develops a series of public policies for forest manage
ment, which is summarized in Rural development pol
icy, Protection against fire and air pollution, Conser
vation of biodiversity, adaptation to climate change,
Competitiveness of forestry investigation. The ac
tions are on outgoing of these policies include sub
sidy programs for forest plantations, payment of en
vironmental services, incentives for the conservation
of nontimber extractions such as cork and acorn col
lection and nonextractive activities such as hunting
activities among others [9]. The changes described
both in the demographic composition of forest re
gions, local economic incentives and regions make
forest management a complex issue, which must be
faced by decisionmakers of public policies as well
as owners of farms with forestry systems throughout
the country. In this sense, develop a framework of
analysis that models the optimal behavior of a forest
state, based on the economic benefits of forest man
agement, is necessary to see the composition and dis
tribution of forest systems. In the case of the study,
the analysis will focus on optimizing the management
of a typical forest state. We develop an economic
analysis framework applied to agroforestry systems
through the implementation of a bioeconomic model
based on a dynamic optimization system. This type
of models have been previously applied to analyze
optimal forest investment decisions, though afforesta
tion programs [22], and or forest management invest
ments [5],[7]. A number of studies apply optimal con
trol techniques in discrete time to solve farmland de
cisions problems whereas the farm owner or holder
has a preference to conserve forest areas within the
farm [3]. This model expands and refines [5] optimal
control model allowing for a more flexible optimiza
tion framework in terms of land use replaceability and
species selection. More specifically, our model in
cludes two more forest species, considering in total
four forest species that can be selected out of eight
possible options at each farm. On the other hand, for
est species can be substituted at the same land plot, for
more profitable species as existing trees reach their
rotation age. This latter change respect to [5] model
is possible as our models allows for the competition
of the four forest species for land that is available for
afforestation/reforestation, which in turn accounts for
both nonforest land uses (i.e. treeless shrubs, pas
tures, crop lands and other) and abandoned forest.
The management of forest represents a complex prob
lem, since forest can be considered as renewable re
sources, at the time they can be exhausted in order
to satisfy other land use demands, such as food pro
duction or land development. A way to achieve an
equilibrium point between these opposites forces is
to propose models that include other forest benefits,
that are normally not considered in conventional eco
nomic analysis. To this end we use mathematical opti
mization tools based on the theory of optimal control.
WSEAS TRANSACTIONS on SYSTEMS and CONTROL
DOI: 10.37394/23203.2022.17.13
Carlos Rodríguez Lucatero,
Marcelo Olivera Villaroel, Paola Ovando
E-ISSN: 2224-2856
111
4 Dynamic optimization
We are going to model our forestry problem under the
methodology of dynamic optimization.
Our basic hypothesis is:
Uniqueness: although each of the phenomena an
alyzed ( land use on farms ) have similar behaviors
and react in a similar way to external and internal
stimuli as can be for instance, subsidies and taxes,
change in relative prices of forest, agricultural and
livestock products that determine the change in land
use on farms, but different departure point. Since con
ditions of distribution of land uses are very different
in each farm, being conditioned by geographic, cli
matic and historical factors, we consider that we can
apply linear optimization processes to our problem.
Our main goal is to optimize the incomes that are in
fluenced by the changes in land uses. There are two
main types of land use in the farms:
• Woodlands
• Treeless
Woodlands are classified by commercial species
that generate income to the owners of the farm.
Commercial species have two handling options:
• Natural character where species regenerate natu
rally, which mean that the species meets its cy
cle of natural reproduction and drive conditions
that allow their reproduction of the species in a
controlled manner, giving the conditions for the
prevalence of the species being replaced by sys
tems of secondary forests or thickets systems,
which are extremely expensive to replace other
land use. So, once the forest is replaced by scrub
systems there is no change in land use.
• The other option is artificial character where the
species are planted and managed in a controlled
manner within the farm
There are some conditions that have to be taken
into account within the process
• We can not have more land than the one provided
from the beginning.
• There are no negative soil quantities
• There are substitution between land uses , but this
substitution is between woodlands and treeless
uses depending on the profitability of each land
use.
• There is an additional limit given in handling land
use scrub , scrub arising from the abandonment
of trees that meet their natural cycle and are not
given the conditions for the prevalence of the
species being replaced by systems of secondary
forests or thickets systems , which are extremely
expensive to replace other land use. So , once the
forest is replaced by scrub systems no change in
land use.
5 Optimal control Model description
In this section we are going to describe the mathemat
ical tools used in the model proposed in [7] as well
as in our model. Both models are discretetime fi
nite temporal horizon optimal control models. Un
der this decision process, the agents try to maximize
the present value profit corresponding to the average
hectare of grassland and can adapt each year his forest
management schedule by means of artificial planta
tions of trees on treeless soil and by natural regenera
tion of existing adult trees. The Matlab code used for
implementing both models admit a maximal temporal
horizon lower than the critical regeneration age from
artificial planted blockhead whose value is T < 142
years. In the following subsections we will describe
the details of the mathematical model and the corre
sponding implementation in Matlab. For the sake of
completeness the datails of the model proposed in [7]
are included in the section 9 of the present article.
5.1 Control vectors of our model
The optimization process of the decisions taken by the
agents involved in the agroforestal process follow the
same procedure of the agents in model proposed in
[7], but we have added in our model two more control
vectors in order to include two new types of tree, the
Stone Pine tree and eucalyptus. The control vectors
are defined as follows:
zt=
ni,t +na
i,t
0
0
.
.
.
−ni,t
.
.
.
0
0
za,t =
ti,t
0
0
.
.
.
−na
i,t
.
.
.
0
0
(1)
vt=
ns,t +na
s,t
0
0
.
.
.
−ns,t
.
.
.
0
0
va,t =
ts,t
0
0
.
.
.
−na
s,t
.
.
.
0
0
(2)
WSEAS TRANSACTIONS on SYSTEMS and CONTROL
DOI: 10.37394/23203.2022.17.13
Carlos Rodríguez Lucatero,
Marcelo Olivera Villaroel, Paola Ovando
E-ISSN: 2224-2856
112
wt=
np,t +na
p,t
0
0
.
.
.
−np,t
.
.
.
0
0
wa,t =
tp,t
0
0
.
.
.
−na
p,t
.
.
.
0
0
(3)
rt=
ne,t +na
e,t
0
0
.
.
.
−ne,t
.
.
.
0
0
ra,t =
te,t
0
0
.
.
.
−na
e,t
.
.
.
0
0
(4)
Where np,t and ne,t represent de Stone Pine tree
and Eucalyptus surface naturally regenerated by
natural planted trees of such species, na
p,t and na
e,t
represent de Stone Pine tree and Eucalyptus surface
naturally regenerated by artifially planted trees
of such species, tp,t and te,t represent the surface
reforested with Stone Pine trees and Eucalyptus by
humans. The rest of the variables still represent the
surfaces mentioned in [7] model. The details of this
model can be consulted in the section 9 of the present
article.
Concerning the utilization of the deforested surface
of land the previous model due to [7], the author
imposed as a constraint just allowing one tree refor
esting species for each deforested type of surface of
a particular farm land. For avoiding these adhoc
rules of reforestation that depended on the particular
farm j. We eliminated the farm numbering unlike
what is done in [7] model where there are 4different
types of farm and in each case one or two species of
tree can be sown and in some cases only one type of
tree can be sown in a specific treeless terrain type.
See the equation 23 in section 9 of the present article
where we describe the precedent model of [7]. From
our point of view this was a very unnatural constraint
and we modified this criteria allowing to reforest
these surface by any species of tree using a weighted
asignation of surface to each species that depends on
their corresponding economic rentability.
ut=
((ξt−(ti,t +ts,t +tp,t +te,t)/4,
(ξt−(ti,t +ts,t +tp,t +te,t), ts,t)/4,
(ξt−(ti,t +ts,t +tp,t +te,t)/4,
(ξt−(ti,t +ts,t +tp,t +te,t)/4
(5)
Where ξtis a transition function of soil utilization and
is calculated as follows
ξt={x320,t−1+xa
320,t−1+s196,t +sa
194,t
+p320,t−1+pa
320,t−1+e320,t−1+ea
320,t−1
(6)
5.2 The state vector movement law of our
model
Introducing the new state variable as well as the new
control vectors the movement law equations that de
scribe the temporal evolution of the controlled state
variables are the following:
xt+1 =T xt+zt+1
xa,t+1 =T xa,t +za,t+1
pt+1 =Qpt+wt+1
pa,t+1 =Qpa,t +wa,t+1
et+1 =Ret+rt+1
ea,t+1 =Rea,t +ra,t+1
st+1 =Anst+vt+1
sa,t+1 =Aasa,t +va,t+1
dt+1 =dt+ut+1
=∀t∈[0, T −1], T < ∞
(7)
Where T, Anand Aahave the same meaning of the
one in the model proposed in [7] and where the matrix
Q, and Rare the transition matrix, corresponding to
stone pine and eucalyptus reforesting land and where
the elements are equal to zero excepted the elements
just under the main diagonal whose values are equal
to 1. The dimensions of Qand Rare 321 ×321.
5.3 The control variable constrains of our
model
The global resource restriction that establish that the
sizes of each farm stay constant, is now expressed as
follows:
Sj=i(T2)·(xt+xa,t) + i(T6)
·(pt+pa,t) + i(T4)·st
+i(T4)·sa,t +i(T8)·(et+ea,t)+
i(4) ·dt
∀j∈ {1,...,4}and ∀t∈ {1,...,T}
(8)
Where i(n)is a row vector of ones of dimension n
and Sjis the productive agricultural surface of farm
WSEAS TRANSACTIONS on SYSTEMS and CONTROL
DOI: 10.37394/23203.2022.17.13
Carlos Rodríguez Lucatero,
Marcelo Olivera Villaroel, Paola Ovando
E-ISSN: 2224-2856
113
jfor each j∈ {1, . . . 4}. Additionally some con
strains regarding the natural regeneration process are
imposed in such a way that the surface to be regen
erated cannot be greater than the surface of trees in
crital age. These constrains are expressed as follows:
ni,t ≤x249,t−1
na
i,t ≤xa
249,t−1
ns,t ≤s143,t−1
na
s,t ≤sa
141,t−1
np,t ≤p249,t−1
na
p,t ≤pa
249,t−1
ne,t ≤e249,t−1
na
e,t ≤ea
249,t−1
=∀t∈[0, T −1], T < ∞(9)
By other side, the author of [7] mentioned that given
the initial distribution of soil utilization, it can be ob
served that very big portions of land in each farm are
available for being reforested. Given this situation,
it can happen that the reforestation resources can be
exhausted in one period of decision. For avoiding
this case some constraints imposed on the number of
reforestation hectares should be compliant with the
spanish forestry and public laws. The corresponding
proposed upper bounds are represented by the follow
ing inequalities:
ni,t +na
i,t +ti,t ≤Sn
i,j
250
∀t∈[0, T −1], T < ∞
ns,t +na
s,t +ts,t ≤Sn
i,j
144
∀t∈[0, T −1], T < ∞
np,t +np
i,t +ti,t ≤Sn
i,j
250
∀t∈[0, T −1], T < ∞ ∧ j∈ {1,...,4}
ne,t +ne
i,t +ti,t ≤Sn
i,j
250
∀t∈[0, T −1], T < ∞ ∧ j∈ {1,...,4}
(10)
5.4 Control Optimization Objective
function of our model
The objective function of a control optimization
model represents the base on which a deciding agent
takes decisions. In the case of this agroforestry prob
lem the function whose value is maximized is a cubic
polynomial that evaluates the profit net present value
of the reforested land depending on the soil utiliza
tion. This function has the following mathematical
form:
fi=aiτ3+biτ2+ciτ+di∀τ∈ {0,...,320}Holm oak
fs=asτ3+bsτ2+csτ+ds∀τ∈ {0,...,196}blockhead
fp=apτ3+bpτ2+cpτ+dp∀τ∈ {0,...,320}stobe pine
fe=aeτ3+beτ2+ceτ+de∀τ∈ {0,...,320}eucalyptus
(11)
Based on these functions it can be calculated the land
price for hectare row vectors for each type of planta
tion and stored in the following vectors:
fi= (fi(0), . . . , fi(320))
fp= (fp(0), . . . , fp(320))
fs= (fs(0), . . . , fs(196))
fs,a = (fs(0), . . . , fs(196))
fd= (pct, pp, pm, po)
(12)
For enabling the agents to take decisions, it can
be calculated the level of annual return vectors using
diferent costing methods for each kind of soil utiliza
tion. The following row vectors store this informa
tion:
yi,n = (yi,n
0, . . . , yi,n
320)
yi,a = (yi,a
0, . . . , yi,a
320)
yp,n = (yp,n
0, . . . , yp,n
320)
yp,a = (yp,a
0, . . . , yp,a
320)
ys,n = (ys,n
0, . . . , ys,n
196)
ys,a = (ys,a
0, . . . , ys,a
194)
ye,n = (ye,n
0, . . . , ye,n
320)
ye,a = (ye,a
0, . . . , ye,a
320)
ys,d = (ys
ct, ys
p, ys
m, ys
o)
(13)
Assuming that the deciding agent maximizes the
finite horizon net present value of the profit associ
ated to the average hectare of the farm j, the objective
function gets the following form:
Fj=1
Sj{∑T
t=0 δt[yi,nxt
+yi,axa,t +ys,n st+ys,asa,t +yp,npt
+yp,apa,t +yq,d dt]
−[fi(x0+xa,0) + fss0+fs,asa,0
+fp(p0+pa,0) + fdd0]
+δT[fi(xT+xa,T ) + fssT
+fs,asa,T +fp(pT+pa,T ) + fddT]}
(14)
q=
i(ilex) ∀j
p(pine) ∀j
e(eucalyptus) ∀j
s(suber) ∀j
5.5 Our model control optimization
Finally the control optimization problem can be writ
ten as the maximization of the objective function as
follows:
max
{ni,t,na
i,t,ns,t ,na
s,t,np,t ,na
p,t,ti,t ,ts,t }t=T
t=1
Fj(15)
subject to: equations (25)(26)(27)(28)(29) and the
given intial conditions x0, xa,0, s0, sa,0, p0, pa,0, d0as
well as the time horizon T.
6 Matlab implementation of our
model
As was mentioned in section 1 our model differs from
the model proposed in [7] in two ways. The first is the
WSEAS TRANSACTIONS on SYSTEMS and CONTROL
DOI: 10.37394/23203.2022.17.13
Carlos Rodríguez Lucatero,
Marcelo Olivera Villaroel, Paola Ovando
E-ISSN: 2224-2856
114
introduction of an additional tree species, the stone
pines. The details of our mathematical model are de
scribed in section 5. We also mentioned at the final
part of the section 1 that our Matlab implementation
generalizes the model of [7] in the sense while this
model perform an optimization of the decision pro
cess just over one fixed scenario, our proposal first
generate all the possible scenarios, then process all
them in order to calculate for each one their present
value profit and finally extract the scenario with re
spect to the obtained present value profit. In this way
our model performs two kind of optimization, the first
is over the decision process in each scenario and the
second is on the selection of the scenario that has the
best economic performance. The steps followed by
our model can be expressed algorithmically as fol
lows:
1) Table_of_scenarios = Generation_of_scenarios;
2) for i=1 to rows(Table_of_scenarios)
3) decision_optimization(Table_of_scenarios(i));
4) end
5) Table_of_scenarios =
Sort_by_present_value(Table_of_scenarios);
6) Best_Scenario = Table_of_scenarios(last_row);
7) Output(Best_Scenario);
Each step of this process will be described in the
following subsections.
6.1 Definition and representation of a
scenario
In our model we have four tree species oak, block
head, stone pine, eucalyptus and four types of treeless
land. Associated to each tree species we have four
reforestation land when the trees become old. We
have as well four types of treeless land. In the model
proposed in [7] the land occupied by some type
of tree can only be reforested by the same kind of
tree. In our model we do not take into account this
constraint. Then we can reforest a treeland piece
that was occupied by oaks with stone pines. By
other side, if we reforested some treeland piece
with some type of tree, we cannot use this same
type for reforesting another reforestation treeland
piece. In what concerns the four types of utilizations
of treeless land, we can repeat the utilization type
in this kind of surfaces without restriction. Then
ascenario is composed by three parts. The first
part is associated with way we use the reforestation
land, the second part is associated with way we use
the treeless land, and third part is associated with
the economic performance of the utilization of the
land utilization of the reforestation land and the
treeless land. The first vector position is associated
with the oak reforestation land, the second vector
position with the blockhead reforestation land, the
third vector position is associated with stone pine
reforestetion land, the fourth vector position with the
eucalyptus reforestation land, etc. For representing a
scenario in Matlab we use a numerical vector with
nine places. The first four positions are used for the
reforestation land, the following four positions for
the treeless land, and the last position for the present
value economic performance. We assigned numbers
for codify each reforestation species as well as for
each treeless utilization as follows:
For reforestation land
species code
Oak 2
blockhead 4
stone pine 6
eucalyptus 8
For treeless land
use code
culture 9
Pastureland 10
scrub land 11
other 12
So an instance of scenario could be
4 6 8 2 10 10 9 9 1,000,000
This vector means that the oak reforestation land will
be occupied by blockhead, the blockhead reforesta
tion land will be occupied by stone pines, the stone
pine reforestation land will be occupied by culture,
the culture reforestation land will be occupied by
oaks, the scrub land as well by the pastureland, the
pastureland will be used for pastureland and the
culture land and the scrub land will be used for
culture. The last number means that the present value
profit of this scenario is 1,000,000 Euros.
6.2 Generation of the table of all possible
scenarios
Sometimes we have to face problems where an ex
haustive exploration of all the possible cases in the
state spaces is necessary and not only to calculate the
total number of possible configurations in such state
space. For example, looking for permutations of a
given set of objects. This kind of tasks are called by
some authors as enumeration problems but this term
can be mistaken as the calculation of the number of
possibilities and for this reason it is more appropri
ate to call them listing problems. Such is the kind of
problem that we need to face for the generation of all
possible scenarios of initial conditions in the land dis
tribution of the model. In this section the generation
of all the possible scenarios to be evaluated will be
WSEAS TRANSACTIONS on SYSTEMS and CONTROL
DOI: 10.37394/23203.2022.17.13
Carlos Rodríguez Lucatero,
Marcelo Olivera Villaroel, Paola Ovando
E-ISSN: 2224-2856
115
described. For this purpose we will have to gener
ate a table whose lines will be numerical vectors of
9positions where the 4first positions represent the
assignment of tree species to areas for artificial refor
estation, the next 4positions will correspond to the
assignment of types of crops to areas of nonforest use
and the last position will be used to record the evalua
tion of the cost at present value of said area allocation.
In the case of the available areas used for reforesta
tion, if one area is occupied by a species of tree, then
the other areas that remain available cannot be used to
reforest with the same species of tree, that is, in this
case we will be talking about permutations. no repe
tition. On the other hand, in the areas for agricultural
use, we can repeat crops, that is, we will be talking
about permutations with repetition. The number of
different ways of occupation of the reforestation land
can be calculated as the number of permutations with
out repetition of the four tree species.
The concepts of combinatorial analysis that we
will use to generate the different scenarios of land use
and reforestation are those of the rule of product, per
mutations without repetition and combinations with
repetition. We will define these concepts below that
are taken from [14].
Definition 1. (The combinatorial rule of product):
If an action of choice can be split into a first and a
second independent stages, and if there are mpos
sible outcomes for the first stage and if, for each of
these outcomes, there are npossible outcomes for the
second stage, then the total procedure can be carried
out, in the designated order, in m×nways.
This principle is sometimes called the principle of
choice. Based on the definition 1 we can be define the
following concepts.
Definition 2. (Variations): A variation of the kth
class of n elements is an ordered kelement group
formed from a set of nelements. The elements are
not repeated and depend on the order of the group’s
elements (therefore arranged).
Vk(n) = n×(n−1) ×(n−2)×
. . . ×(n−k+ 1) = n!
(n−k)!
(16)
Definition 3. (Permutation with repetitions): Same
as the definition 2 but when the objects can be re
peated. The calculation of nobjects in kpossible
positions in a linear array (using the definition1) is
P R(n, k) = n×n×n×. . . ×n
| {z }
k
=nk(17)
Details about such notions of combinatorial anal
ysis can be found in chapter 1 of book [14].
Thus, we can calculate the total number of possible
permutations without repetition of nobjects making
use of the formula of definition 2.
In our case n= 4 and k= 4. So the number of
possible land uses for reforestation would be V4(4) =
24.
In addition, the number of different ways of oc
cupation of the treeless land can be calculated as the
number of permutations with repetitions of nobjects
taken from kto kusing the formula of definition 3.
In our case n= 4 and k= 4. So the number of
possible land uses for non reforestation utilization of
the land is 44= 256. From all the above and applying
the product combinatorial principle, we have that the
number of possible scenarios and therefore the num
ber of rows in the table will be 24×256 = 6144. Once
we know the total number of possible scenarios that
correspond to the rows, we can proceed to the gen
eration of the 6144 rows of the table. Combinatorial
analysis allows us to be able to calculate the size that
a problem space can have. However, sometimes this
is not enough but it is necessary to explicitly obtain
the complete list of these possibilities. That is exactly
the problem that concerns us.
To do this, we first generate the list of all the per
mutations without repetition of the 4numbers that
correspond to the codes used to denote each species of
tree. For this, a Matlab program inspired by the Lex
icographic order permutation generation algorithm
that can be found in [16] (page 319), was imple
mented. The algorithm such as it is described in in
[16] (page 319) is
•L1.[Visit] Visit the permutation a1a2. . . an
•L2.[Find j] Set j←n−1. If aj≥aj+1 decrease
jby 1repeatedly until aj< aj+1. Terminate the
algorithm if j= 0.
•L3.[Increase aj] Set l←n.If aj≥aldecrease l
by 1repeatedly until aj< al. Then swap aj↔
al.
•L4.[Reverse aj+1 . . . an] Set k←j+ 1 and l←
n. Then while k < l interchange ak↔aland
set k←k+ 1, l ←l−1. Return to L1
The algorithm correction is assured as step L2 has
in the index jthe smallest index such that it has been
visited to this point of processing all possible permu
tations starting with the prefix a1a2. . . ajsuch that
the following permutation in lexicographic order will
make ajto have a bigger value, in L3 before the inter
change aj↔alwe have that aj+1 ≥. . . ≥al−1≥
al> aj≥al+1 ≥. . . ≥anand after the interchange
we have aj+1 ≥. . . ≥al−1≥aj> al≥al+1 ≥
. . . ≥an. According to D. Knuth in his book [16] this
WSEAS TRANSACTIONS on SYSTEMS and CONTROL
DOI: 10.37394/23203.2022.17.13
Carlos Rodríguez Lucatero,
Marcelo Olivera Villaroel, Paola Ovando
E-ISSN: 2224-2856
116
algorithm dates back to the 14th century in India, was
proposed by Narayana Pandita at that time and has
been rediscovered in other periods of the past such
as in the preface of the book Specimen Analyticum
by Lineis Curvis Ordinis by C.F. Rüdiger (Leipzig
1748). There are other more efficient algorithms to
solve the problem of listing all possible permutations
without repetition inspired by algorithms for obtain
ing ntuples of Gray codes one of which is the al
gorithm of Adjacent Interchanges whose detailed de
scription can be found at [16] page 320, as well as
another algorithm called Algorithm P(plain changes)
[16], page 321 which originates from the program
ming of the Church Steeples in England in the 17th
century (1653) (Cambridge Fortyeight, Ernest Mor
ris The History and Art of Change Ringing (1931),
R. Duckworth and F. Stedman, Tintinnalogia(1668)).
There are many algorithms for generating the list of
permutations without repetition, one of which can be
known as (Plain changes) or Algorithm P and that
can be consulted on page 322 of [16]. This algorithm
dates from the 17 century in England and was used by
workers in church steeples to program the melodies
in a more varied way from generating all the possible
permutations of the 5bells that the church towers had.
Regarding the generation of the list of the different
forms of occupation of the unfolded land, which we
represent in the positions 5to 8of our row vector of
9positions that make up each row of the aforemen
tioned table, the situation is much simpler when deal
ing with permutations with repetition. The generation
of this list of permutations with repetition in lexico
graphic order of 4tuples of 4possible positions is a
very simple algorithm since it would consist of nested
cycles in the following way
for i=1 to 4
for j=1 to 4
for k=1 to 4
for l=1 to 4
build-4-tuple(i,j,k,l);
Once the 44possible 4tuples of positions have
been generated, we choose the elements of the set
of codes associated with the ungrounded land uses,
which in this case are the values 9,10,11,12 and with
these values are filled the positions 5to 8of the row
vector of the table. For the purposes of our forestry
application it is sufficient to generate the list of pos
sible permutations in lexicographical order. However
it is interesting to mention that other forms of gener
ation of ntuples have been studied and are addressed
in chapter 7, subsection 7.2of the book [16]. Some
of these algorithms are related to Gray codes where it
is sought that the ntuples generated follow an order
such that between one ntuple and the next they only
differ by a single element as shown in the following
sequence of 4tuples in binary:
0000 0001 0011 0010
0110 0111 0101 0100
1100 1101 1111 1110
1010 1011 1001 1000
Such codes have application in the conversion of ana
log information to digital format and vice versa since
they minimize the errors that can be generated in said
format conversions. This type of coding was used
in color television broadcasting and pulse modulation
applications as well as in the analog transmission of
digital signals. Many variants of Gray codes and their
applications can be found in chapter 7of the book
[16].
6.3 Sorting the table and selection of the
best scenario
Once all the possible values of the 8first positions
of the row vector of the table have been generated,
we proceed to evaluate each of the scenarios corre
sponding to each row to fill the 9position of the row
that corresponds to the economic evaluation at present
value of said stage. It is important to mention that
some of these scenarios may violate the restriction of
the total land availability and that is why we submit
the table to a filter of these cases to eliminate them.
Once this purification has been carried out, we order
the table in ascending order with respect to the eco
nomic evaluation obtained. To do this, we apply a
standard algorithm for sorting by selection and after
this we will have in the last position of the table the
option or scenario whose economic evaluation was
the highest. Inglés
We should mention that the selection sort algo
rithm has a temporary performance in the worst case
of O(n2)where nis the number of rows to sort. This
algorithm is efficient enough to sort the table, how
ever for later versions we could use a more efficient
sorting method whose temporal performance is of the
order of O(nlog(n)) as is the case of MergeSort or
the Quicksort algorithm.
7 Results
We carry out our simulation by calling the dynamic
optimization routine for each row of our table of all
possible scenarios for assigning uses of land, both
wooded and bare. The time window used was T=
100, annual rate r= 0.03 and with these parameters
we call our simulator 6144 times that correspond to all
the possible scenarios that we generate. Once this is
done, we display the graphs corresponding to the best
scenario and the worst scenario for comparison pur
poses. The results are shown in the following graphs.
The starting land occupation (t= 0) for each call to
the dynamic optimization routine is:
WSEAS TRANSACTIONS on SYSTEMS and CONTROL
DOI: 10.37394/23203.2022.17.13
Carlos Rodríguez Lucatero,
Marcelo Olivera Villaroel, Paola Ovando
E-ISSN: 2224-2856
117
Figure 1: Worst and Best scenario. The blue line is the
total surface, black line represent Oak, red line rep
resent Stone pine, magenta line represent Blockhead,
yellow line represent Eucalyptus and green represent
pasture land.
Land Use Surface in Ha
Natural Oak tree 283
Artificial Oak tree 85
Natural blockhead tree 174
Artificial blockhead tree 88
Natural stone pine tree 195
Artificial stone pine tree 120
Natural Eucalyptus tree 45
Artificial Eucalyptus tree 530
culture 200
Pastureland 150
scrub land 80
other 50
After having carried out the evaluation of each sce
nario assigned to each row of the table, the results
shown are shown in the following figures.
8 Discussion and conclusions
In this article we propose a mathematical agroforestry
optimization model and from this we implement a
simulation in MATLAB. The model works with four
species of trees and four types of agricultural land
Figure 2: Screenshot of Best and Worst scenario ob
tained by the simulator
use and through the simulation we explore all pos
sible combinations of agricultural land use for refor
estation purposes to obtain the one that would be the
most economically profitable. The simulation of each
possible scenario took approximately three minutes
and exploring all the possibilities represented having
a laptop working for approximately a week. For this
reason, we found it necessary to take a significant
sample of all possible cases. Despite this, the time
to process all the cases in the sample took around 48
hours of continuous processing. For this reason, we
are considering the possibility that in future work we
will make use of parallel programming to try to re
duce processing time even more. Another possibility
for future work is to try to apply other combinatorial
optimization techniques to try to reduce the number
of cases to be explored.
Acknowledgments
Paola Ovando acknowledges the support of the
European Commission under the Marie Curie
IntraEuropean Fellowship Programme (PIEF2013
621940). Carlos Rodríguez Lucatero and Marcelo
Olivera Roel want to acknowledge the support given
by the Universidad Autónoma Metropolitana Unidad
Cuajimalpa.
WSEAS TRANSACTIONS on SYSTEMS and CONTROL
DOI: 10.37394/23203.2022.17.13
Carlos Rodríguez Lucatero,
Marcelo Olivera Villaroel, Paola Ovando
E-ISSN: 2224-2856
118
9 Appendix
9.1 State variables of David Martín Barroso
model
This model was presented in the Phd. Thesis of
Dr. David Martín Barroso in at the University Com
plutense of Madrid from Spain in the year 2003 [7].
The state variables of this model were
xt=
x0,t
.
.
.
x320,t
xa,t =
xa
0,t
.
.
.
xa
320,t
(18)
st=
s0,t
.
.
.
s197,t
sa,t =
sa
0,t
.
.
.
sa
197,t
(19)
dt=
d1,t
d2,t
d3,t
∀t∈ {0,1,...,T}, t < ∞(20)
Where xt, xa,t represent the surface amount of
natural and artificial generated Holm oak,st, sa,t
represent the surface amount of natural and artificial
generated blockhead,d1,t cereal cultivated surface,
d2,t pastureland, and d3,t scrub surface.
9.2 Control vectors of David Martín
Barroso model
The decisions taken by the agents involved in the
agroforestal process are represented by control vec
tor variables that influence the time evolution of the
state variables as usual in any optimization model.
This variables are called control vectors. There are
six control variables, three correspond to the regener
ation and reforestation decisions concerning the Holm
oak, and the other three to the blockhead. They are
mathematically represented as follows:
zt=
ni,t +na
i,t
0
0
.
.
.
−ni,t
.
.
.
0
0
za,t =
ti,t
0
0
.
.
.
−na
i,t
.
.
.
0
0
(21)
vt=
ns,t +na
s,t
0
0
.
.
.
−ns,t
.
.
.
0
0
va,t =
ts,t
0
0
.
.
.
−na
s,t
.
.
.
0
0
(22)
Where ni,t represent de Holm oak surface natu
rally regenerated by natural planted trees of such
species, na
i,t represent de Holm oak surface naturally
regenerated by artifially planted trees of such species,
ti,t represent the surface reforested with Holm oak
by humans. The variables ns,t, na
s,t and ts,t represent
the corresponding surfaces for the blockhead. The
negative entries correspond to the amount of trees
that have achieved their productive age. In what
concerns the dismasted land uses, the calculations
depend on the cases of study and because of that
the calculations are adhoc. In his thesis Dr. David
Martin Barroso [7] talks about four specific farms
(j= 1, . . . , 4). The corresponding vectors are:
ut=
(ξt,−ti,t,0) forj= 1
(0,−(ti,t +ts,t), ξt)for j= 2
(ξt,−(ti,t +ts,t),0) for j= 3
(0,(ξt−ti,t), ts,t )for j= 4
(23)
Where ξtis a transition function of soil utiliza
tion and is calculated as follows
ξt={x320,t−1+xa
320,t−1for j= 1
x320,t−1+xa
320,t−1+s196,t +sa
194,t ∀j= 1
(24)
9.3 The state vector movement law of David
Martín Barroso model
The movement law equations that describe the tem
poral evolution of the controled state variables are
the following.
xt+1 =T xt+zt+1
xa,t+1 =T xa,t +za,t+1
st+1 =Anst+vt+1
sa,t+1 =Aasa,t +va,t+1
dt+1 =dt+ut+1
=∀t∈ {0,...,T−1}, T < ∞
(25)
Where T, Anand Aaare transition matrix where all
the elements are equal to zero excepted the elements
just under the main diagonal whose values are equal
to 1. The dimensions of Tare 321 ×321,Anare
197 ×197 and Aaare 195 ×195.
WSEAS TRANSACTIONS on SYSTEMS and CONTROL
DOI: 10.37394/23203.2022.17.13
Carlos Rodríguez Lucatero,
Marcelo Olivera Villaroel, Paola Ovando
E-ISSN: 2224-2856
119
9.4 The control variable constrains of David
Martín Barroso model
In [7] they assume that the farm land owners don’t
change the size of their terrains by means of comer
cial transactions in such a way that the total surface
of each farm stay constant along the time. This global
resource constrain can be expresed by the following
equation
Sj=i(321)·(xt+xa,t)+i(197)·st+i(197)·sa,t+i(3)·dt
(26)
where Sjrepresent the total surface of the farm j∈
{1, . . . , 4}, for t∈ {1, . . . , T }, and i(n)is a row vec
tor of nones.
The author of the model in [7] also assumes some
non negativity constrains on the control variables ex
pressed mathematically as follows
ni,t, na
i,t, ns,t, na
s,t, ti,t, ts,t ≥0(27)
Additionally some constrains regarding the natural re
generation process are imposed in such a way that the
surface to be regenerated cannot be greater than the
surface of trees in crital age. That is expressed as fol
lows:
ni,t ≤x249,t−1
na
i,t ≤xa
249,t−1
ns,t ≤s143,t−1
na
s,t ≤sa
141,t−1
ti,t ≤d2,t if j= 1
ti,t +ts,t ≤d2,t if j= 2
ti,t +ts,t ≤d2,t if j= 3
(ts,t ≤d3,t)∧(ti,t ≤d2,t−1+ξt)if j= 4
=∀t∈[0, T −1], T < ∞
(28)
By other side, the author of [7] mentioned that given
the initial distribution of soil utilization, it can be ob
served that very big portions of land in each farm are
available for being reforested. Given this situation,
it can happen that the reforestation resources can be
exhausted in one period of decision. For avoiding
this case some constraints on the number reforesta
tion hectares should be imposed for being compliant
with the spanish forestry and public laws. The corre
sponding proposed upper bounds are represented by
the following inequalities:
ni,t +na
i,t +ti,t ≤Sn
i,j
250
∀t∈[0, T −1], T < ∞ ∧ j∈ {1,...,4}
ns,t +na
s,t +ts,t ≤Sn
i,j
144
∀t∈[0, T −1], T < ∞ ∧ j= 2,4
ns,t +na
s,t +ts,t ≤0.25·d2,0
144
∀t∈[0, T −1], T < ∞ ∧ j= 3
(29)
9.5 Objective function of David Martín
Barroso control optimization model
The objective function of a control optimization
model represents the base on which a deciding agent
takes decisions. In the case of this agroforestry prob
lem the function whose value is maximized is a cubic
polynomial that evaluates the profit net present value
of the reforested land depending on the soil utiliza
tion. This function has the following mathematical
form:
fi=aiτ3+biτ2+ciτ+di
∀τ∈ {0,...,320}Holm oak
fs=asτ3+bsτ2+csτ+ds
∀τ∈ {0,...,196}blockhead
(30)
Based on these functions it can be calculated the
land price for hectare row vectors for each type of
plantation and stored in the following vectors:
fi= (fi(0), . . . , fi(320))
fs= (fs(0), . . . , fs(196))
fs,a = (fs(0), . . . , fs(196))
fd= (pct, pp, pm)
(31)
For enabling the agents to take decisons, it can
be calculated the level of annual return vectors using
diferent costing methods for each kind of soil utiliza
tion. The following row vectors store this informa
tion:
yi,n = (yi,n
0, . . . , yi,n
320)
yi,a = (yi,a
0, . . . , yi,a
320)
ys,n = (ys,n
0, . . . , ys,n
196)
ys,a = (ys,a
0, . . . , ys,a
194)
ys,d = (ys
ct, ys
p, ys
m)
(32)
Assuming that the deciding agent maximizes the
finite horizon net present value of the profit associ
ated to the average hectare of the farm j, the objective
function gets the following form:
Fj=1
Sj{∑T
t=0 δt[yi,nxt
+yi,axa,t +ys,n st+ys,asa,t +yq,ddt]
−[fi(x0+xa,0) + fss0+fs,asa,0+fdd0]
+δT[fi(xT+xa,T ) + fssT+fs,asa,T +fddT]}
(33)
q={i(ilex) for j= 1,3
s(suber) for j= 2,4
9.6 David Martín Barroso control
optimization model
Finally the control optimization problem can be
written as the maximization of the objective function
as follows:
max
{ni,t,na
i,t,ns,t ,na
s,t,ti,t ,ts,t }t=T
t=1
Fj(34)
WSEAS TRANSACTIONS on SYSTEMS and CONTROL
DOI: 10.37394/23203.2022.17.13
Carlos Rodríguez Lucatero,
Marcelo Olivera Villaroel, Paola Ovando
E-ISSN: 2224-2856
120
subject to: equations (25)(26)(27)(28)(29) and the
given intial conditions x0, xa,0, s0, sa,0, d0as well as
the time horizon T.
9.7 State Variables of our model first
variation
Based on the model mentioned in the last subsections
our first variation on this model consited in adding
one additional tree species, the Stone Pine tree. For
this end we added the respective three decision adding
the corresponding control variable vectors. We also
have added more flexibility to the model of [7] by
generalizing with variables the corresponding reno
vation tours for ecah species of trees. Fisrt of all the
new set of state variables are represented mathemati
cally as follows:
xt=
x0,t
.
.
.
xT1,t
xa,t =
xa
0,t
.
.
.
xa
T2,t
(35)
st=
s0,t
.
.
.
sT3,t
sa,t =
sa
0,t
.
.
.
sa
T4,t
(36)
pt=
p0,t
.
.
.
pT5,t
pa,t =
pa
0,t
.
.
.
pa
T6,t
(37)
dt= d1,t
d2,t
d3,t !∀t∈ {0,1, . . . , T }, t < ∞(38)
References:
[1] Azqueta, D. Introducción a la economiía ambi
ental. McGrawHill Madrid 2002.
[2] Bielsa, I.; Pons, X.; Bunce, B. Agricul
tural abandonment in the north east
ern Iberian Peninsula: The use of ba
sic landscape metrics to support plan
ning. Journal of Environmental Planning
and Management 2005, 48(1), 85–102.
https://doi.org/10.1080/0964056042000308166
[3] Blot, J.; Hayek,N. Infinitehorizon optimal
control in the discretetime framework. New
York: Springer. (2014)
[4] Campos, P.; Ovando P.; Mesa B.; Oviedo J.
L. Environmental Income of Live stock Graz
ing on Privatelyowned Silvopastoral Farms in
Andalusia, Spain. Land Degra dation and De
velopment. 2016.
[5] Cerdá, E.; MartínBarroso, D. Optimal
control for forest management and conser
vation analysis in dehesa ecosystems. Euro
pean Journal of Operational Research (2013).
https://doi.org/10.1016/j.ejor.2012.12.010
[6] Crow, T. R., Host, G. E.; Mladenoff, D.
J. Ownership and ecosystem as sources
of spatial heterogeneity in a forested
landscape, Wisconsin, USA. Land
scape Ecology (1999), 14(5), 449–463.
https://doi.org/10.1023/A:1008084123874
[7] Martínez Barrosso D. Instrumentos de análi
sis económico en sistemas agroforestales. PhD
Thesis, Universidad Complutense de Madrid
2003.
[8] https://ec.europa.eu/web/forestry,2013
[9] http://eurnexo/legal
content/ES/TXT/?uri=LEGISSUM:l60040,
(revised January 2018), 2018
[10] Faustmann, M. The Faustmann
Model (part 1) Introduction to
Forestry,Forest Policy, and Economics.
web.archive.org/web/20111229211645/ 2011.
[11] Flichman G.; Allen T. Bioeconomic
modeling:Stateoftheart and key priori
ties. agris.fao.org 2015.
[12] GilTena, A.; Saura, S.; Brotons, L. Ef
fects of forest composition and structure
on bird species richness in a Mediter
ranean context: Implications for forest
ecosystem management. Forest Ecology
and Management (2007), 242(2–3), 470–476.
https://doi.org/10.1016/J.FORECO.2007.01.080
[13] GómezLimón, J.; De Lucío Fernández, J. V.
Changes in use and landscape preferences on
the agriculturallivestock landscapes of the
central Iberian Peninsula (Madrid, Spain).
Landscape and Urban Planning (1999), 44(4),
165–175. https://doi.org/10.1016/S0169
2046(99)000201
[14] Grimaldi R.P. Discrete and Combinatorial
Mathematics: An Applied Introduction Third
Edition AddisonWesley 1994.
[15] Kienast, F.; Frick, J.; van Strien, M. J.;
Hunziker, M. The Swiss Landscape Moni
toring Program A comprehensive indicator
set to measure landscape change. Eco
logical Modelling (2015), 295, 136–150.
https://doi.org/10.1016/j.ecolmodel.2014.08.008
WSEAS TRANSACTIONS on SYSTEMS and CONTROL
DOI: 10.37394/23203.2022.17.13
Carlos Rodríguez Lucatero,
Marcelo Olivera Villaroel, Paola Ovando
E-ISSN: 2224-2856
121
[16] Knuth D. The art of Computer Programming:
Combinatorial Algorithms, Part 1, Volume 4A,
AddisonWesley, 1997
[17] LasantaMartínez, T.; VicenteSerrano, S. M.;
CuadratPrats, J. M. Mountain Mediter
ranean landscape evolution caused by the
abandonment of traditional primary activities:
A study of the Spanish Central Pyrenees.
Applied Geography (2005), 25(1), 47–65.
https://doi.org/10.1016/j.apgeog.2004.11.001
[18] Lasanta, T.; Beguería, S.; GarcíaRuiz,
J. M. Geomorphic and Hydrological
Effects of Traditional Shifting Agricul
ture in a Mediterranean Mountain Area,
Central Spanish Pyrenees. Mountain Re
search and Development (2006), 26(2),
146–152. https://doi.org/10.1659/0276
4741(2006)26[146:GAHEOT]2.0.CO;2
[19] Newman D.H. The Optimal Forest Rotation: A
Discussion and Annotated Bibliography. Tech
nical Report USDAForest Service (1988).
[20] Pyy J., Ahtikoski A., Laitinen E., Siip
ilehto J. Introducing a NonStationary
Matrix Model for StandLevel Optimiza
tion, an EvenAged Pine (Pinus Sylvestris
L.) Stand in Finland forests, MDPI Open
Journal,http://www.mdpi.com/journal/forests,
Forests 2017, 8, 163; doi:10.3390/f8050163
[21] Pyy J., Ahtikoski A., Lapin A., Laitinen E.
Solution of Optimal Harvesting Problem by
Finite Difference Approximations of Size
Structured Population Model Mathematical
an Computational Applications, MDPI Open
Journal, http://www.mdpi.com/journal/mca,
Math. Comput. Appl. 2018, 23, 22;
doi:10.3390/mca23020022
[22] Ryan, M.;C. ODonoghe; Phillips, H. Mod
elling Financially Optimal Afforestation and
Forest Management Scenarios Using a Bio
economic Model, Open Journal of Forestry
(2016), 06(01), 19–38.
[23] VicenteSerrano, S. M.; Lasanta, T.; Cuadrat,
J. M.; Cuadrat, J. M. Transformaciones en
el paisaje del Pirineo como consecuencia del
abandono de las actividades económicas tradi
cionales. Pirineos (2000), 155(0), 111–133.
https://doi.org/10.3989/pirineos.2000.v155.91
[24] Ovando, P., Beguería, S., and Campos, P.
Carbon sequestration or water yield? The
effect of payments for ecosystem services
on forest management decisions in Mediter
ranean forests. Water Resources and Eco
nomics, Elsevier, Volume 28, October 2019,
100119,https://doi.org/10.1016/j.wre.2018.04.002
[25] Ovando, P., P. Campos, B. Mesa, A. Alvarez,
C. Fernández, J. Oviedo, A. Caparros, and B.
AlvarezFarizo Renta y capital de estudios de
caso de fincas agroforestales de Andalucía. In
P. Campos and P. Ovando (Eds.), Renta total
y capital de las fincas agroforestales de An
dalucía, pp. 156–445. Madrid: Memorias cien
tíficas de RECAMAN. Vol 4. Editorial CSIC.
Contribution of individual authors to
the creation of a scientific article
(ghostwriting policy)
Carlos Rodríguez Lucatero carried out the math
ematical modelling as well as the MATLAB im
plementation of the simulator. He also contributed
with the writing of the present article. Marcelo
Olivera Villarroel contributed in the modelling of
economic aspects as well as forestry aspects of the
mathematical model on which the simulator is based.
He also contributed with the writing of the present
article. Paola Ovando contributed in the modelling of
forestry aspects of the mathematical model as well as
with the forestry data on which the simulator is based.
Sources of funding for research
presented in a scientific article or
scientific article itself
Paola Ovando acknowledges the support of the
European Commission under the Marie Curie
IntraEuropean Fellowship Programme (PIEF2013
621940). Carlos Rodríguez Lucatero and Marcelo
Olivera Roel want to acknowledge the support given
by the Universidad Autónoma Metropolitana Unidad
Cuajimalpa.
Creative Commons Attribution
License 4.0 (Attribution 4.0
International , CC BY 4.0)
This article is published under the terms of the
Creative Commons Attribution License 4.0
https://creativecommons.org/licenses/by/4.0/deed.en_US
WSEAS TRANSACTIONS on SYSTEMS and CONTROL
DOI: 10.37394/23203.2022.17.13
Carlos Rodríguez Lucatero,
Marcelo Olivera Villaroel, Paola Ovando
E-ISSN: 2224-2856
122