Introduction to population dynamics and resource exploitation
EE. Canuto, D. Mazza, Introduction to population dynamics and resource exploitation January 29, 2021 1
Introduction to population dynamics and resource exploitation
Enrico Canuto and Daniele Mazza Former Faculty, Politecnico di Torino [email protected] [email protected] Version 1.1, September 2019
Table of Contents
Introduction to population dynamics and resource exploitation ................................................ 1
Abstract ...................................................................................................................................... 2
Acronyms ................................................................................................................................... 3 Introduction ........................................................................................................................ 3 Single-species growth and mortality: the logistic equation ............................................... 7
Logistic equation and free response ............................................................................ 7
Parameter estimation of the logistic rate ..................................................................... 9
Application to US crude oil production ..................................................................... 12 Two-species competition: the Lotka-Volterra equation ................................................... 16
The Lotka-Volterra equation ..................................................................................... 16
Equation derivation and meaning ....................................................................... 16
Equilibrium points and conservation equation ................................................... 17
Other proofs that trajectories are closed ............................................................. 18
Average population ............................................................................................ 19
Perturbation equation ......................................................................................... 20
Graphical plots ................................................................................................... 20
The data that suggested V. Volterra equation .................................................... 21
Orbit period................................................................................................................ 22
Introduction ........................................................................................................ 22
S-D. Shih derivation ........................................................................................... 22 . Canuto, D. Mazza, Introduction to population dynamics and resource exploitation January 29, 2021 2
Estimation of the LV equation parameters ................................................................ 26
Method ............................................................................................................... 26
Monte Carlo test ................................................................................................. 28
Discrete time LV equation ......................................................................................... 31
Forward Euler method ........................................................................................ 31
Stabilization and frequency tuning ..................................................................... 31
Application of LV equation: snowshoe hare and Canadian lynx data ...................... 33
Application of LV equation: resource exploitation ................................................... 37
A limiting form of the LV equation ................................................................... 37
A first-order equation of the capital stock driven by production ....................... 40
Application to experimental data: California gold rush ..................................... 42 Riferences ......................................................................................................................... 44 Appendix .......................................................................................................................... 46
Poincaré maps ............................................................................................................ 46
Abstract
The paper was suggested by a brief note of the second author about the application of the Hubbert’s curve to predict decay of resource exploitation. A further suggestion came from the interpretation of the Hubbert’s curve in terms of a specific Lotka-Volterra (LV) equation. The link with population dynamics was obvious as logistic function and LV equation were proposed within the demography science field. Mathematical population dynamics has a history of about two centuries. The first principle and model of population dynamics can be regarded the exponential law of T. R. Malthus. In the XIX century, the Malthusian demographic model was first refined to include mortality rate by B. Gompertz. In the early XIX century the model was further refined by P-F. Verhulst by introducing the standard logistic function. The previous models only concern the population of a single species. In the early XX century, the American demographer A.J. Lotka and the Italian mathematician V. Volterra proposed a pair of state equations which describe the population dynamics of two competing species, the predator and the prey. This paper is concerned with the single and two-species fundamental equations: the logistic and LV equation. The paper starts with the generalized logistic equation whose free response is derived together with equilibrium points and stability properties. The parameter estimation of the logistic function is applied to the raw data of the US crude oil production. The . Canuto, D. Mazza, Introduction to population dynamics and resource exploitation January 29, 2021 3 paper proceeds with the Lotka-Volterra equation of the competition between two species, with the goal of applying it to resource exploitation. At the end, a limiting version of the LV equation is studied since it describes a competition model between the production rate of exploited resources and the relevant capital stock employed in the exploitation.
Acronyms
AS asymptotically stable, asymptotic stability DT discrete-time LV Lotka-Volterra LTI linear time-invariant RF recovery factor US United States of America Introduction
The paper was suggested by a brief note of the D. Mazza about the application of the Hubbert’s curve (in other terms, the derivative of the logistic function) to predict decay of resource exploitation [1]. A further suggestion came from the interpretation of the Hubbert’s curve in terms of a specific Lotka-Volterra (LV) equation by U. Bardi and A. Lavacchi [2]. The link with population dynamics was obvious as logistic function and LV equation were proposed within the demography science field. Mathematical population dynamics has a history of about two centuries. The first principle and model of population dynamics can be regarded the exponential law [3] of T. R. Malthus [1766-1834], which is the free response of the state equation , Vol dx t bx t x t xdt . (1) In (1),
Vol denotes the unit of the population volume x , usually equal to the number of individuals or to their concentration in a region, and b is the growth rate in Vol/unit , where unit is a generic time unit. In the XIX century, the Malthusian demographic model was first refined [4] to include mortality rate (the decay rate of resource exploitation) by B. Gompertz [1779-1865] The corresponding state equation is We will use the symbol [Vol] to denote the population volume measuring unit. An alternative symbol would be [Count]. . Canuto, D. Mazza, Introduction to population dynamics and resource exploitation January 29, 2021 4 max 0 0 ln , Vol dx t xbx t x t xdt x , (2) where max lim t x x t is the asymptotic limit and b is proportional to the initial growth rate. The free response, known as the Gompertz’s function, is given by exp ln / 1 exp x t x b x bt . (3) By recalling the algebraic limit maxmax nn xx nn x x , (4) replacement of max ln / x x in (2) with the LHS of (4), provides the generalized logistic state equation ([5], also Richards’ equation, from F.J. Richards who proposed the equation in 1959) n dx t b x x t x t xdt n x . (5) The free response of (5) for n , max max / 2 x t x and t , holds maxmax maxmax exp1 tanh2 21 exp b t tx xx t b t t . (6) The function in (6) was proposed in the early XIX century [6] by P-F. Verhulst [1804-1849] as an alternative refinement of (1) and is known as the standard logistic function , not to be confused with the symmetric logistic rate x t . The term ‘logistic’ was created without explanation by P-F. Verhulst, presumably on the Ancient Greek λογῐστῐκή (the art of computing, a branch of the Ancient Greek Mathematics). For max t t the population rate x t grows and for max t t decays to zero. The symmetric logistic rate x t is also known as Hubbert’s curve (and so it will be referred here) in the field of Earth’s resource exploitation, since it was proposed by the American geologist M. King Hubbert in 1956 [7], thus establishing a link between population dynamics and resource exploitation. Gompertz’s and logistic equations (2) and (5) are autonomous state equations [17], since their solution only depends on initial conditions, and not on external perturbations and events. To obtain a finite asymptotic population limit max lim t x x t , an intrinsic mortality term (known also as self-competition), i.e. the negative quadratic term in (5), is subtracted from+ the growth term, i.e. the positive and linear term in (5). Gompertz’s and Richards’ equations only concern the population of a single species. In the early XX century, the American demographer A.J. Lotka (1880-1949) in 1920 [8] and the . Canuto, D. Mazza, Introduction to population dynamics and resource exploitation January 29, 2021 5 Italian mathematician V. Volterra (1860-1940) in 1926 [9] proposed a pair of state equations which describe the population dynamics of two competing species, the predator whose volume is denoted by x and the prey whose volume is x : predator: , 0prey: , 0 x t b a x t x t x xx t b a x t x t x x , (7) The coefficients in (7) are positive. Each of the two equations follows the same template of the logistic equation (5) for n , but the mortality term (the negative term in the second equation of (7)) of the prey and the growth term of the predator (the positive term in the first equation of (7)) are made extrinsic since they depend on the population volume of both species, thus accounting for their competition. The LV equations (7) have been extended to include intrinsic mortality terms as in (5): predator: , 0prey: , 0 x t b a x t a x t x t x xx t b a x t a x t x t x x , (8) with positive coefficients. Equation (8) can be extended to a finite number n of competing species. By denoting the volume vector with x , and using the element-wise product symbol , the dim n LV equation becomes , 0 t A t t x b x x x x , (9) where b is the vector of the intrinsic birth (positive) or death (negative) rates and the n n matrix A with components ij a accounts for relationships between species. The diagonal term ii a is usually negative and prevents the species to grow indefinitely. Two species , i j x x are competing when ij a and ji a . A species i x is a predator at the expense of j x , when ij a and ji a . The two-dimensional equation (8) can be viewed as the linear case of the generic predator-prey equation [10] predator: , , 0prey: , , 0 x t b x c x x x x xx t x b x c x x x x x , (10) where b is the intrinsic predator mortality, b x is the per capita growth rate of the prey in absence of predator, , 0 c x x is the per capita consumption rate of the predator, known as the trophic function (from the assumed Ancient Greek τρόφικος, relative to τροφή, food, nourishment). , c x x determines, through the efficiency , the mortality rate of the prey. A common assumption as in (7) is that the trophic function is only dependent on the prey volume, i.e. c c x . Based on experimental data, R. Arditi and L. R. Ginzburg . Canuto, D. Mazza, Introduction to population dynamics and resource exploitation January 29, 2021 6 suggested in [10] that actually c depends on the ratio between prey and predator volumes, i.e. / c c x x , and that c is a concave monotonic function tending to as asymptote. Since the formalization in 1972 [11] by the British theoretical biologist J. Maynard Smith [1920-2004] (following a verbal statement of the American geneticist G.R. Price [1922-1975]) of the evolutionary stable strategies - a Nash equilibrium of the game theory which is evolutionary stable, in the sense that only natural selection protects the population evolution from small external perturbations - population dynamics has been complemented by the theory of games in order to model and explain population evolution within the evolutionary game theory [12]. This paper is just concerned with the single and two-species fundamental equations: the logistic and LV equation. Section 2 is devoted to the generalized logistic equation in (5). In Section 2.1, the free response (or logistic function) is derived together with equilibrium points and stability properties. The parameter estimation of the logistic function in Section 2.2 will be restricted to the standard logistic function and will be applied in Section 2.3 to the raw data of the US crude oil production from 1860 to 2018 [13]. Section 3 is devoted to the Lotka-Volterra (LV) equation of the competition between two species, with the goal of applying it to resource exploitation as suggested by [2]. To this end, the well-known LV equation (7) is studied in Section 3.1 pointing out the properties of the periodic free response and of the relevant closed trajectories around the so-called critical point. Section 3.2, which can be omitted at a first reading, derives, with the help of the literature [14], a method for computing the period of any LV closed trajectory. Section 3.3 provides a three-step method for estimating the six parameters of the LV equation (7), namely the four parameters , , / , / b b a b a b of the equation and initial conditions. The first step aims to estimate the four parameters from explicit equations based on population measurements. At this stage, initial conditions are assumed to be equal to initial measurements. The second step refines the previous estimates by minimizing an error function. However, any two slightly different LV equations, although periodic and stable around their critical point, will drift from each other like two sine functions tuned on slightly different frequencies. To remove drift, the estimated equation must be tuned on the measurement data by means of a dynamic feedback, which is driven by the model error between measurements and the reconstructed response. This refinement is done by a third estimation step which is based on a stable discrete-time version of the LV equation, the topic being treated in Section 3.4. The section will point out that, although stable around the same critical equilibrium of the continuous-time equation, also the DT version will produce a drifting but bounded response. The drift can be removed by a dynamic feedback as already mentioned. Section 3.5 will apply the estimation procedure of Section 3.3 to the well-known ‘lynx and hare data’ collected around the early XX century in the northern Canadian forests [16], although recent experimental data suggest that the LV model . Canuto, D. Mazza, Introduction to population dynamics and resource exploitation January 29, 2021 7 is too simple. Finally, Section 3.6 studies the limiting version of the LV equation, which, as suggested in [2], describes a competition model between the production rate of exploited resources and the relevant capital stock employed in the exploitation. It is shown that the capital stock is the solution of a linear time-invariant (LTI) equation driven by the production rate. This allows to design a simple method for estimating equation parameters, namely, the capital obsolescence rate and the transformation coefficient between resources and capital. The method is applied to the sparse and uncertain data of the California gold rush, [18], [19], with some interesting result. Single-species growth and mortality: the logistic equation
Logistic equation and free response
Let us consider the generalized logistic equation in (5) and apply, under n , the variable change max / 1 n z t x t x , (11) and the relevant differential identity n x t x tdx dx dzdz n x x x x nz . (12) Starting from (5) and by applying (11) and (12), we obtain, after some manipulations, a new equation as follows: n x t x t x tb x t xx n x xx tx z tz n bn x tz x z t b z t z t z . (13) Free response.
The solution of the last equation in (13), which is LTI, is the following exp exp1 1 tt z t b t t x b b t dz t t , (14) where exp t t b t t and the change of variable in (11) provides the free response of (5) . Canuto, D. Mazza, Introduction to population dynamics and resource exploitation January 29, 2021 8 max 1/max 00 nn xx t x t tx . (15) The response (15) reduces to the standard logistic function in (6) for n and / 2 x x . Equilibrium points and stability.
The last equation in (13) has a single natural equilibrium z , which is asymptotically stable (AS), since (14) provides lim 1 t z t and shows that the eigenvalue of the LTI perturbation equation is b . Unlikely, the original equation (5) (the first equation in (13)) has two equilibrium points, x x x . (16) The former equilibrium in (16), which corresponds to the divergent lim t z t , is unstable, since the corresponding perturbation equation with / 0 b n and x x holds , 0 bx t x i x xn . (17) The second equilibrium max x x corresponds to z and is left to the reader to prove that is AS under the LTI perturbation equation of (5) and to find the relevant eigenvalues. Generalized Hubbert’s curve.
Last but not least, we derive from (15) a new expression of the logistic rate x t in (5), which holds: maxmax 00 1/ 1max 00
11 1 n nn xx t txx t b x t tx . (18) The maximum value max x and the argument * arg max x x x x can be derived from the derivative identity / 0 dx x dx , which from (5) holds: max 1/*max max1/ 1maxmax n nn dx x b xndx n xx xxx x nx bx n . (19) The time argument max arg max t t x t is obtained from the free response (15) and holds . Canuto, D. Mazza, Introduction to population dynamics and resource exploitation January 29, 2021 9 *max 0 0 t t x x . (20) In the standard logistic case, corresponding to the Hubbert’s curve, we obtain * 0 maxmax max / 2/ 4 x x xx bx . (21) Replacement of n x x n in (18), and of t with max t provides the generalized and standard logistic rate expressions max max1/ 1maxmax max 2
1 @ 11 nmqx x n t tx t b n t tx t tx t b nt t . (22) Figure 1 shows three profiles of the normalized logistic rate max ( ) / x t x with max
10 unit t and
1/ unit b . The unitary exponent n separates two types of skew-symmetric profiles: 1) fast rising profiles with n , 2) fast decaying profiles with n . Figure 1. Normalized logistic rate for different n . Parameter estimation of the logistic rate
We will restrict to the standard logistic rate in (22) (the second equation), which is a function of three parameters collected in max max , , x b t p , where max x is the asymptotic population volume, b is the growth rate coefficient and t is the time instant of peak logistic rate. The functions in (6) and (22) suggest that knowledge of max t allows to estimate two parameters of p , namely . Canuto, D. Mazza, Introduction to population dynamics and resource exploitation January 29, 2021 10 max maxmax max max max
22 / 4 / x x tb x t x t x x . (23) The critical step is the estimation of max t , i.e. of the time instant when max max x t x and the second time derivative x t becomes zero. From (22) (second equation), the second derivative holds max max2max max3max t t t tx t x b x tt t . (24) When experimental data of the logistic rate x t are just increasing, max t cannot be directly estimated from the raw data, but one can assume that either the last measurement corresponds to the maximum value of the second derivative or the relevant time instant can be obtained from the first derivative of the experimental data, after a convenient interpolation. Let us denote the this time instant with inf t (it is the inflection point of the increasing logistic rate): it is the instant when the second derivative of the logistic rate, that is x t , becomes equal to zero. The third derivative holds
22 max max3max max 4max t t t tx t x b t t t t , (25) and by setting the derivative to zero, provides the interval max inf t t as follows x t t t b b . (26) From (26) we can obtain the value of the logistic rate: inf max / 6 x t bx (27) from the identity inf max t t . Since identities (26) and (27) are just two equations with the three unknowns of max max , , x b t p , we need a third identity, possibly only involving experimental data of x t . For instance, we can search the interval half max t t such that half inf / 2 x t x t and half inf t t . The relevant identity provides the second-degree equation
22 2 2 half max
10 1 0, =exp b t t , (28) whose solution, in the case of half max t t , holds max half t t b b . (29) . Canuto, D. Mazza, Introduction to population dynamics and resource exploitation January 29, 2021 11 By combining (26), (27) and (29), and by assuming to know inf arg max t t x t and half t from the identity half inf / 2 x t x t , one obtains a first-trial estimate p of the parameters in p , namely inf halfmax infmax inf b t tx y t bt t b , (30) where the concave mark denotes measurements, the convex mark estimates, and notations have been simplified by the identity y t x t . We assume also that the raw data i i y t x t , i N , have been interpolated to obtain uniform sampling times i t iT . Given the m -th step parameter estimate , 0,1,... m m p , the sampled measurements y i x i of the logistic rate and the estimated logistic rate , , y i m x i m p p , the estimation RMS error can be constructed as follows Ni s k y i y i mN p , (31) and can be progressively minimized with respect to p . The final RMS error is indicated by s , and it can be normalized by the mean value y i defined by Ni y y iN . (32) When the covariance of the measurement error is assumed to be known, the Cramer-Rao bound (the lower bound of the covariance matrix [17]) P p of the parameter vector p , can be estimated from the gradient , , / t y t g p p p of the logistic rate in (22) (second row), computed at p p (the parameter estimate) and at the sampling sequence , 0,..., 1 i t iT i N . The gradient holds t t t ty tt y t b t t bx b t t t t pg p pp . (33) If we assume a statistical independent measurement error , y i y i y i p with known variance y i , the estimated covariance is computed as
11 20 , ,
N Tyi
P i i i p g p g p , (34) where the variance pk of a single parameter , 1, 2, 3 k p k holds pk kk P p . (35) . Canuto, D. Mazza, Introduction to population dynamics and resource exploitation January 29, 2021 12 We will use also the normalized standard deviation (std) / pk kk k P p p . The covariance in (34) can be a-posteriori checked though a Monte Carlo run [17]. Application to US crude oil production
We apply the estimation method in Section 2.2 to a well-known set of raw data, the US annual production of crude oil from 1860 to 2019 [13], i.e. from the first successful use of a drilling rig in a well of Pennsylvania. The well was just drilled for producing oil by E.L. Drake and associates on August 27, 1959. The success of the Drake’s well quickly led to oil drilling in other locations of United States and to a first wave of investments in oil drilling, refining and marketing. The principal product of oil was kerosene which replaced whale oil for illumination purposes. Before Drake’s success, oil was an accidental byproduct of wells drilled for salt brine. Generally speaking, petroleum reservoirs are permeable rock formations with minute porous spaces that contain trapped crude oil, natural gas or condensate (a low-density mixture of hydrocarbon liquids present in natural gas fields). The accumulated oil or gas is prevented to move upwards by an impermeable cap rock. Under favorable conditions, usually at the early stage of exploitation (the so-called primary recovery), oil production is made possible by drilling a well that produces a pressure difference between well and oil reservoir, thus allowing hydrocarbons to flow toward the well bore. The recovery factor (RF) expresses the amount of crude oil (or gas) which has been extracted as a percentage of the estimated hydrocarbon content of the reservoir. The primary recovery RF varies between 5 and 10%. After the primary recovery, pressure depletion alone becomes insufficient and injection of water or gas is required to keep the reservoir pressure at a convenient level for crude oil recovery. This process is known as the secondary recovery, which allows to raise the RF to a range between 15 and 40%. This implies that reservoirs are abandoned after the secondary recovery with, on the average, two thirds of their hydrocarbon content still left on the ground. The production raise in Figure 2 until the seventies of the XX century and the subsequent decay indicate the progressive exploitation and cessation of the secondary recovery in the US wells. In the last twenty years, corresponding to the new increase of US crude oil production in Figure 2, a pair of unconventional exploitation processes have been implemented. The first one, referred to as tertiary recovery or end-of-reservoir process, aims to increase the RF of the abandoned wells to about the 60% of the estimated hydrocarbon content. The second process addresses rocks with low-permeability, like shale rocks (organic rich fine-grain sedimentary rocks), in which case hydraulic fracturing can stimulate the flow. It should be remarked that one of the first source of mineral oil, since prehistoric times, was oil shale, and oil shale industry, started in France in 1837 and steadily grew since before World War I, but declined, after World War II, due to accessibility of the Middle East crude oil. . Canuto, D. Mazza, Introduction to population dynamics and resource exploitation January 29, 2021 13
Raw data and spline interpolation are shown in Figure 2. The interpolation time unit is 0.1 year. The ordinate is
10 m /y ,which was obtained from the conversion of the original thousand US barrels per day. Figure 2.
Interpolated and raw data of US crude oil production. An immediate inspection of Figure 2 reveals that around the year 2000 a different kind of US crude oil production was implemented as explained above, which implies that at least two different logistic production rates y t x t and y t x t must be summed up, namely y t y t y t , (36) where y reached the maximum value around the year 1980, whereas y is still increasing. A cascade estimation of the parameters of , y t y t has been done. The parameter vector , , x b t p of y has been estimated by restricting the raw data to end,1 i t t (from inspection) and the parameter , , x b t p of y from the residuals of , y i y i y i p , but restricted to start,2 i t t . The initial parameter estimates p and p have been obtained with the methods of Section 2.2. More specifically p has been computed from (23) and the estimate of max,1 t from data inspection. Instead, p has been computed from (30). The final parameter estimates have been obtained by minimizing the estimation RMS error in (31). The MATLAB fminsearch(.) function has been employed. . Canuto, D. Mazza, Introduction to population dynamics and resource exploitation January 29, 2021 14 Figure 3.
Raw data, estimated logistic rate profile and estimation residuals. Figure 3 compares interpolated raw data y i with the estimated profile y of the total function in (36). The estimation error profile y i is also shown. Estimated parameters and their first-step estimated values are reported in Table 1. Table 1. US crude oil production estimation No Parameter Symbol Unit Value Comments First logistic rate y
1 Limit of the cumulative production max,1 x
10 m b y max,1 t y 1976.5 Same 4 Limit of the cumulative production max,1 x
10 m b y max,1 t y 1976.2 Same 7 Estimation RMS error s
10 m /y 16.6 end,1 t t
8 Fractional RMS error / s y fraction 0.079 same . Canuto, D. Mazza, Introduction to population dynamics and resource exploitation January 29, 2021 15
9 Normalized parameter std pk fraction <0.005
1, 2, 3 k (each parameter) Second logistic rate y
10 Limit of the cumulative production max,2 x
10 m b y max,2 t y 2022 Same 13 Limit of the cumulative production max,2 x
10 m b y max,2 t y 2020.5 Same 16 Estimation RMS error s
10 m /y 25.2 start,2 i t t
17 Normalized parameter std pk fraction <0.1
1, 2, 3 k (each parameter) Total logistic profile 178 Estimation RMS error s
10 m /y 17.7 19 Fractional RMS error / s y fraction 0.076 Initial and final estimates of the first logistic profile y look very close since fairly all the profile is available from raw data. Something different occurs to the second profile y , since the rate peak time in Table 1, row 14, has been anticipated, by the RMS error minimization, of two years, with respect to the first-step estimate in Table 1, row 11, obtained from (30) and the guess inf,2 t . The difference is due to the intermediate peak at about the 2015 epoch, which likely hides an intermediate logistic profile. Only future production data may confirm or contradict the prediction. The a priori fractional parameter standard deviation pk has been computed from (35), by assuming a uniform measurement error std, which has been set equal to the RMS of the estimation error in Table 1, rows 7 and 16. The last assumption should be conservative since . Canuto, D. Mazza, Introduction to population dynamics and resource exploitation January 29, 2021 16 the estimation error, as one can see from Figure 3 (the green line), looks more a model error than a measurement error. Two-species competition: the Lotka-Volterra equation
The Lotka-Volterra equation
Equation derivation and meaning
The well-known Lotka-Volterra (LV) equation of the predator-prey competition was independently proposed by the American demographer A.J. Lotka (1880-1949) in 1920 and by the Italian mathematician V. Volterra (1860-1940) in 1926. A.J. Lotka proposed the model to describe a hypothetical chemical reaction in which chemical concentrations oscillate [1]. V. Volterra proposed the model to explain the increase of predator fish (and the corresponding decrease of prey fish) in the Adriatic Sea after the World War I [9]. If we refer to predator-prey species competition, only two species are considered, predator k and prey k . Their population volume or size [Vol] defines the state variables, which combine in the vector , x x x . Their per capita growth/decay rate is defined by / ln / , 1, 2, 0 k k k k k g t x t x t d x dt k x t . (37) Let us assume a proportional decay (predator) and growth (prey) rate
0, 0 t g g A t g x x x b x , (38) where , , 0 k b b b b [Vol/s] is the intrinsic growth/decay vector, decay for predator and growth for the prey, i.e. the growth/decay rate in absence of competition. The matrix A [(Vol/s)/Vol], the interaction matrix, describes the effect of the population volume on the per capita growth. The resulting autonomous equation , 0 t t t x g x x x x , (39) where denotes the Hadamard element-wise product, is known as the Kolmogorov’s predator-prey equation. In the LV equation, only mutual interactions are accounted for, which leads to ij aA aa . (40) Let us remark that A in (40) has the same form of the state matrix of an undamped oscillator. By combining (37), (38) and (40), the LV state equation reads as , 0 , 0 t A t t t x b x x x x x . (41) The equation components are as follows: . Canuto, D. Mazza, Introduction to population dynamics and resource exploitation January 29, 2021 17 predator: , 0 , 0prey: , 0 , 0 x t b a x t x t x x x tx t b a x t x t x x x t , (42) where all the coefficients and variables are nonnegative. The equation block-diagram is shown in Figure 4. The three loops, the decay and growth loops around the predator and prey state variables, respectively, and the overall oscillating loop from predator to prey are shown. The dashed lines may account for environmental actions. x t x t a b x x b a Predator b Prey Decay loopOscillating loopGrowth loop
Figure 4.
Block-diagram of the LV equation (42).
Equilibrium points and conservation equation
Equation (42) has a pair of equilibrium points: the zero equilibrium
0, 0 x and the critical (equilibrium) point , / , / x x b a b a x . (43) To proceed to the analysis, (42) is rewritten in terms of the normalized population volumes , / , / z z x x x x z , (44) which transformation provide the following nonlinear equation z t bt tz t bt z t x z t x z z z z zx z . (45) The LV equation (45) cannot be integrated in closed-form, but a constant integral is found by multiplying the first equation in (45) by b z z , the second one by b z z and by summing their products to get the following identity
1 / 1 / 0 b z z z b z z z , (46) which is equivalent to the differential equation . Canuto, D. Mazza, Introduction to population dynamics and resource exploitation January 29, 2021 18 ln ln 0 d b z z b z zdt . (47) In turn, the differential equation (47) can be integrated to provide the constant function ln lnln ln 0, 0
V b z t z t b z t z t Vb z z b z z t z z , (48) which is clearly positive. Equation (48) is a conservation equation , which is typical of oscillatory systems: at any time, the integral V t V z z is conserved. This implies that the trajectories of (45), lying in the positive quadrant t z Z = of the state space, are periodic. Prey population increases when predator population decreases and vice versa.
Other proofs that trajectories are closed
The gradient and the Hessian matrix of V z V b z b zV b z b z zz zz , (49) show that [1,1] V z is the unique minimum value of V z in the positive quadrant, and all the trajectories are contour curves of V z around V b b E z . This implies that U V V z z z , i.e. that ln 1 ln 10, 0
U b z t z t b z t z tU E U zz z , (50) is a candidate Lyapunov function and E can be defined as the ‘orbit energy’ in [Vol/s] units. The energy has been set to zero in the critical point. The time derivative has been already computed in (46) and holds zero
1 / 1 / 0
U t b z z z b z z z z , (51) which confirms that trajectories are closed around the critical point x . A second proof comes from the use of Poincaré map (see Section 5.1). To this end, let us rewrite V V z z with the help of (48), as follows ln / ln / b z z b z z b z z b z z , (52) and let us select the line , 1, 2 k k z z k , as the return surface k . With the notation / k k k z z , the identity (52) simplifies to . Canuto, D. Mazza, Introduction to population dynamics and resource exploitation January 29, 2021 19 ln 1 j k j j z , (53) and shows that the solution j j j z z always exists and indicates that the return point coincides with the initial point as for closed trajectories, in agreement with the Poincaré map P z z . The unit value j implies that the crossing point is unique and is just the initial point. In other terms, the trajectory is tangent to the line k . Intermediate crossing points exist and depend on the value of j . Under j , the intermediate crossing point holds , k j k j z z z , whereas under j it holds , k j k j z z z . A significant property is that the intermediate crossing point component j k z does not depend on k z and thus on k . A third proof comes by converting U z in (50) into a Hamiltonian function . The nonlinear state transformation ln exp , 1, 2 k k k k s z s z k changes U z into a Hamiltonian function as follows exp 1 exp 1 0 H b s s t b s t s t s . (54) To prove that H s is Hamiltonian, let us transform the state equations (45) in terms of , s s s : z t b z z s t b sz t b z z s t b sHs t s ssHs t s ss s s , (55) where the last two state equations constitute a Hamiltonian system. Since the Hamiltonian in (54) does not explicitly depend on time, the energy conservation law states that the Hamiltonian is equal to the ‘constant energy’ of the orbit, namely E H s . Average population
A further meaning of the critical point , 0 x x x is that of providing the average population volume along any closed trajectory. In fact, let us consider a generic trajectory of period T such that T z z and from (45) the equation / 0ln , 0/ 0 z z bd t tz z bdt z z z z z . (56) Integration of (56) along the interval T provides . Canuto, D. Mazza, Introduction to population dynamics and resource exploitation January 29, 2021 20
10 02 T bT d Tb z z z z , (57) and consequently, as expected, T T d dT T z z x x . (58)
Perturbation equation
The perturbation equation about the critical point is expected to be oscillatory. Let us define the perturbation
1, 1 z z z z z of the normalized state z in (44). Equation (45) can be rewritten as z t bt tz t b z z z z . (59) The LTI part of the state matrix in (59) has a pair of imaginary eigenvalues LV LV j b b j j f , which do not depend on the prey-predator interaction coefficients
12 21 , a a . The relevant period 1 / LV LV
P f is just the period of the small orbits close to the critical point, since the orbit period will depend on the initial state z encoded in the orbit energy E defined in (50). Graphical plots
Figure 5 shows the time profile of the pair (predator, prey) by using the , x x x variables of equation (7) and starting from the following initial condition: x x x . (60) The equation coefficients and the population average have values: A b x . (61) . Canuto, D. Mazza, Introduction to population dynamics and resource exploitation January 29, 2021 21 Figure 5.
Left: time profiles of the largest trajectory on the right figure. Right: LV trajectories about the critical point (equal to the population average) [4,8].
The data that suggested V. Volterra equation
The weakness of the LV equation is of being autonomous , thus incapable of accounting for external perturbations. In fact, the LV equation suggested by V. Volterra [9], could not per se explain, as the output of a cause-effect phenomenon, the predator fish fraction (mainly of the clade Selachimorpha, like sharks) increase in Northern Adriatic Sea during and after the World War I, since the cause should have been related to war and specifically to a diminished fishing activity. The record of the predator fish fraction over the whole caught fish per year [16] is shown in Figure 6. Figure 6.
Raw and interpolated data of the predator fish fraction.\ The available data are scarce and moreover only the predator fish fraction is available, which implies that the LV equation (7) should be rewritten as follows. Let
X x x the total amount of caught fish per year, and / , 1, 2 k k x X k the fraction of predator and prey. It is straightforward to rewrite (7) as follows predator: 1 , 0total: , 0 t b t t a X t t tX t t X t X X , (62) where the intrinsic rate perturbation t and the initial fish population X are not available. A similar equation holds for the prey prey: 1 , 0 t b t t a X t t t , (63) where, now, the rate perturbation t due to time varying fishing, affects b with the opposite sign of b in (62). In conclusion, fishing reduction implies t , and consequently an increase of the prey intrinsic growth rate, now b t , and a reduction of the predator . Canuto, D. Mazza, Introduction to population dynamics and resource exploitation January 29, 2021 22 intrinsic mortality rate, now b t . Scarce and incomplete data prevent further considerations and numerical estimation. Orbit period
Introduction
This section may be omitted at a first reading. We have already noticed that the perturbation period LV P b b is just the period of the small orbits about the critical point / , / b a b a x . The issue is to find the expression of a generic period P E as a function of the initial conditions, which are encoded in the orbit energy E . Some equivalent expressions can be found in the literature as reported by S-D. Shih in [14], who in addition derived the same expression as that of F. Rothe [15], but in a simpler way. S-D. Shih derivation
We report the steps of the S-D. Shih‘s derivation. Let us begin from U z in (50), which is set equal to the ‘orbit energy’ U E z and then exponentiated as follows exp ln 1 / ln 1 exp /exp / exp / 1 /exp b b z z b b z z E bz z E b b bz z . (64) As a second step, let us solve for z by separating the variables as follows /1 2 2 2 21 1 1 1 exp / / exp / ,exp , 0 b b h z z z E E b z Eh z z z z , (65) where h z is shown in Figure 7, z being a generic nonnegative real, and E is the energy at the critical point. Since h z is not one-to-one invertible, it must be written as the sum of two invertible functions h z h z h z , which are defined as follows exp , 0 1exp , 1 h z z z zh z z z z . (66) . Canuto, D. Mazza, Introduction to population dynamics and resource exploitation January 29, 2021 23 Figure 7.
The components of ( ) h z . The range of , 0 h z z is bounded and holds h z e . With the help of the notation k k g u h u , we can formally write , , 0 1, , 1 z g z E zz g z E z , (67) where , z E has been defined in (65). Equations in (67) need to be completed with their range z z z , which is obtained from the inequalities /2 2 2 22 21 b b z E z z E E bE z zb . (68) The second row inequalities have two solutions
1, 1 z z , which coalesce into for z z , and are defined by exp 1 /exp 1 / z E g E bz E g E b . (69) By replacing exp / exp 1 / z z E b into (64) one finds that the pair , z z corresponds to the orbit points , P P with coordinates z z and z z , as Figure 8 shows. The same holds for the pair , z E z E and the points , P P with coordinates ,1 z z and , ,1 z z . This implies that the period computation can be written as the sum of time-integrals along the four orbit branches mod 1,4 , , 0,1, 2, 3 k k P P k , once the relation between the time differential dt and the state differential either dz (in this case) or dz has been found. . Canuto, D. Mazza, Introduction to population dynamics and resource exploitation January 29, 2021 24 Figure 8.
The clockwise normalized orbit starting from the red point. The time differentials are obtained from (45) by replacing either z with , 0,1 k k z g z k in the second equation or z with , 0,1 k k z g z k in the first equation, which provides , 11 , k kk k dzdt z b g z Edzdt z b g z E . (70) Using the second differential, the orbit period P E in [s] units, is the solution of the following integral z Ez E dzP E b zg z E g z Edzb zg z E g z E , (71) where g applies to z , g to z and the minus sign to dz . Following [14], the above integral can be converted into a convolution integral through the substitution exp exp 1 / z z b , (72) which entrains the following transformations exp 1 / , 1, 1, exp 1 k z g b g z g zEz E b , (73) and the new integration limits and differential . Canuto, D. Mazza, Introduction to population dynamics and resource exploitation January 29, 2021 25 z z z Ez z ddz b z (74) Replacement of (73) and (74) in (71) and some manipulations allow us to rewrite (71) as the following convolution integral ELV
EP E G G db b , (75) where the dimensionless function to be integrated holds G g g . (76) In the numerical integration of (75) written as NiLV i E iEP E G Gb bEi i N , (77) one should pay attention that G is singular for , the first-order approximation being lim 2 / G . (78) This implies that, unlike in (77), the integration step should become smaller for and E . Figure 9 shows a typical argument of the summation in (77) with a fractional energy quantization / 0.001 E . The fractional error between the period computed by (77) and the period estimated by simulation is about . A fractional error larger than the quantization is likely due to the uniform quantization of the numerical integration. . Canuto, D. Mazza, Introduction to population dynamics and resource exploitation January 29, 2021 26 Figure 9.
Typical argument of (77) and the average (the computed period).
Estimation of the LV equation parameters
Method
Estimation of LV parameters from experimental data is organized into three steps. 1)
Fist we perform the estimation of the population averages , x x and of the rate coefficients , b b from experimental data , , 0,..., 1 x i x i i N , where uniform sampling, i.e. i t iT , is assumed. If necessary, the samples can be interpolated through splines, especially in view of the population rate computation , x i x i and of the uniform sampling. Let us collect the previous estimates and the initial conditions in the vector
10 1 12 1 2 20 2 21 2 1 x b a b x x b a b x p . (79) 2) The second step performs a refinement of the previous estimate through the progressive minimization of an error functional , , , ; m m m x i x i x i x i m p , m , between raw (interpolated) data , x i x i and simulated samples , m m x i x i , where the subscript m denotes a generic optimization step, until some convergence criterion is met. The simulated samples derive from the integration of (42) based on the estimated parameter vector , 0 m m p . The error functional is defined as the RMS of the normalized error between raw and simulated samples as follows.
21 20 1 ;1 N k kmm i k k x i x i mN x p . (80) The normalization is fixed and is done with respect to the raw data average , x x obtained at the first step. The MATLAB fminsearch(.) function will be employed, as already done in Section 2.3. 3) The last step interprets the final residual errors ; , 1, 2 k k k x i x i x i k p of the second step as the result of an external inputs k u i to the autonomous LV equation (42), which therefore is rewritten as predator: , 0 , 0prey: , 0 , 0 x t b a x t x t u t x x x tx t b a x t x t u t x x x t . (81) To this end, a DT state equation of (81) is necessary as it will be done in Section 3.4. secondary, k u i must be estimated as a dynamic feedback driven by k x i as shown in the same Section. . Canuto, D. Mazza, Introduction to population dynamics and resource exploitation January 29, 2021 27 The first-step estimation equations are obtained as follows. Let us rewrite equation (42) at the uniform measurement times i N / 1 // 1 /1 1,
N Ni i x i x i b x i xx i x i b x i xx x i x x iN N , (82) where the round mark denotes measurement and estimate, and the nonzero average k x is estimated from the measurements. If the measurement noise k x i is zero mean, the estimate , 1, 2 k x k is unbiased. Since the long-term averages of the velocity and of k k x x i are zero, the measurement record must be split into two segments , k k S S defined by ; , , ; , k i k k k i k k t x i x i i i t x i x i i i
S S , (83) where / 0 i P T accounts for the estimation transient of k x i , P is the estimated period and is assumed. The end time / i i mP T selects an integer number of periods. The rate coefficients are estimated by separately computing the ratios between the partial sums over the positive and negative segments as follows / /1 / 1 // /1 / 1 / i ii ii ii i i i i it ti it ti i i it ti it t x t x t x t x tb x t x x t xx t x t x t x tb x t x x t x S SS SS SS S , (84) where 1 k k accounts for the uncertainty of two terms in (84), and in absence of any other information 1 / 2 k k . The above estimates can be proved to be unbiased, when measurement errors are zero mean. The lower bound of the covariance T P p p p p p E E E of the estimate p at the second step, may be predicted by the Cramèr-Rao bound P p [17], where p is the true parameter and P p holds
11 10 , ,
N Ti
P G i S i G i p p p . (85) In (85), diag , S i i i is the covariance of the measurement error , i x x i x of the population measurements , i x x i x , under assumption of normality and statistical independence. The matrix , G t p is the Jacobian matrix of the measurement model, i.e. of the free response of the LV equation (81). The columns of the Jacobian matrix . Canuto, D. Mazza, Introduction to population dynamics and resource exploitation January 29, 2021 28 , ,, x t x tG t p pp p p (86) are the solution of the linear time-varying state equations
0, 1, 0, 0, 0, 0, 00, 1, , 0, 0, 0 0, 0, 0, 0,1, 0, 00, 0, 0, 0,1, x t xx xx t a b a xt xx t xx xx t a b a xt x t gp p p pg gp p p pg , (87) which are driven by , t x x x and by the initial conditions. Monte Carlo test
The first estimation step corresponding in (82) and (84) has been checked by a Monte Carlo trial of m M runs [17]. The first, second and third estimation steps will be applied in Section 3.5. The discrete-time state equation of the third estimation step will be explained in Section 3.4. In the Monte Carlo test, the estimation of the initial conditions
10 20 , x x is assumed to be equal to the first measurement, i.e.
10 20 1 2 , 0 , 0 x x x x . (88) The measurement noise has been designed to be low frequency such to modulate the amplitude of the periodic LV response as in Figure 10, left: the noise cutoff frequency 0.13 Hz, 1, 2 wk f k is smaller than the response frequency 0.4 Hz LV f . Figure 10, left, shows the time profile of the deviation k k i x x t , which is clearly disturbed by the low-frequency measurement noise. Figure 10, right, shows the time profile of volume rate error i x t . The volume rate reconstruction is essential for the first-step intrinsic rate estimation in (84). . Canuto, D. Mazza, Introduction to population dynamics and resource exploitation January 29, 2021 29 Figure 10.
Left: deviation from the average population. Right: volume rate estimation error. Table 2 shows test parameters and statistical results. They refer to normalized estimates, i.e. the estimated parameter k p is divided by the true parameter value k p , which implies unitary expected value: / 1 k k p p E , (89) where , , , p b p x p b p x . Table 2 shows Monte Carlo bias k and standard deviation k of each parameter, defined by
10 210
M kk km kM kk km k p mM pp mM p . (90) Table 2. Monte Carlo test parameters and statistics No Parameter Symbol Unit Value Comments Volume rate filter and measurement noise parameters 1 Estimated period P s 2.5 (0.4 Hz) 2 Measurement time unit T s 0.02 3 Volume rate filter eigenvalues , Discrete time 4 Number of used periods 20 5 Low-frequency noise standard deviation , w w Vol
6 Cutoff frequency , w w f f Hz . Canuto, D. Mazza, Introduction to population dynamics and resource exploitation January 29, 2021 30
7 Measurement noise fractional bias / , / w w x x fraction
To avoid negative measurements Monte Carlo test 8 Number of runs M
100 9 Predator intrinsic rate b , fraction Predator average population x , fraction
11 Prey intrinsic rate b , fraction
12 Prey average population x , fraction Monte Carlo statistics in Table 2, rows 9 and 11, shows that the first-step estimation of the intrinsic rates is biased. This is due to the measurement error bias in Table 2, row 7. Figure 11 shows the quantile-quantile (Q-Q) plot of the normalized estimate / , 1, 2 k k b m b k , of the intrinsic rate coefficients versus the Gaussian variables / k g m M defined by ; , / k k k F g m M , (91) where ; , / k k k
F g m M is the cumulative Gaussian distribution defined by the estimation mean k and standard deviation k in (90). The Monte Carlo estimates fit the Gaussian distribution. . Canuto, D. Mazza, Introduction to population dynamics and resource exploitation January 29, 2021 31 Figure 11.
Monte Carlo estimates of the normalized intrinsic rate coefficients versus Gaussian variable.
Discrete time LV equation
Forward Euler method
The simplest discretization method is the forward Euler method which amounts to replace the state derivative in (45) with t T tt T z zz . (92) The corresponding DT equation from (45) holds z i Tbi i iz i Tb z z z z z z . (93) and has the same equilibrium points of the CT equation, namely
0, 0 z and the critical point z . The stability of the trajectories around the critical point depends on T as follows. To this end, let us consider the perturbation z z and the corresponding LTI equation: z i bi i T iz i b z z z x x . (94) The characteristic equation of the LTI part of the state matrix in (94) is the following T b b T b b , (95) and shows that the DT equation is always unstable for T . Stabilization and frequency tuning
When a DT state equation is employed as a state predictor for prediction, control or identification purposes, open-loop stability is not essential, since the state predictor is the result of a feedback stabilization of the DT state equation itself, the feedback being driven by the sampled measurements of the real dynamic system (alternatively, of an accurate simulator). As a matter of fact, weaker feedback corrections will be applied to more accurate DT equations. In the case of nonlinear state equations, no general rule exists for building an exact DT equation, i.e. such that, given the time unit T and the initial condition z , the solution i z of (93) is the same as the sampled solution iT z of (45) for i . To reach this objective, we proceed in two steps. . Canuto, D. Mazza, Introduction to population dynamics and resource exploitation January 29, 2021 32 Stabilization of the perturbed LTI equation. We find the sampled-data equation of the LTI part of (59), which amounts to cos 1 sin1 sin cos 1
LV LVLV LV bT Tbi i i I A T ib T Tb z z z z . (96) We replace the LTI matrix in (93) with the matrix
A T in (96) in order to guarantee the same equilibrium points of (93) and the perturbation stability of (96). We should remark that the free response of (96) around the critical point is only quasi-periodic , since in general the period 2 /
LV LV P imposed by the eigenvalues, is not a multiple of the time unit T . The DT nonlinear equation holds
01 , 00, 1, 2
DT DT DT DT DTDT k DT k k z ii i A T iz ix i z i x k z z z z z z . (97) 2)
Frequency tuning.
In general, equation (97), though stable, drifts in a bounded way (see Figure 12) from the CT equation (45), like a pair of sine functions tuned on slightly different frequencies. In other words, the error equation between CT and DT is Lyapunov unstable but Lagrange stable [17]. The drift is eliminated by adding to (97) an input i u which is synthesized as a dynamic feedback driven by the model error m i DT i t i x x x . (98) The simplest dynamic feedback is first-order as follows u u m u uu m i i L ii i L i x x x x xu x x . (99) The gain matrices diag , L l l and diag , L l l are designed to stabilize the unstable error equation, which can be written as
01 2 i i DTDT i i i t t i X ii X i t X x x x x h x x x u x xx z x , (100) where h is the perturbation to be cancelled. Figure 12 shows the open-loop bounded drift of the error equation (100) affecting the model error m i x in (98). . Canuto, D. Mazza, Introduction to population dynamics and resource exploitation January 29, 2021 33 Figure 12.
Bounded drifting error between CT equation and stable DT equation of the LV dynamics. Figure 13 shows the LF components of i u , i.e. the state vector u i x in (99), in charge of cancelling h in (100) in the unnoisy (left) and noisy case (right). The closed loop frequency BW is 1 Hz. Figure 13. Left: periodic feedback without measurement noise. Right: random feedback in the noisy case.
Application of LV equation: snowshoe hare and Canadian lynx data
The parameter estimation method of Section 3.3 will be applied to a short-time period of the well-known ‘snowshoe hare and Canadian lynx population data’ collected in the northern Canadian forests from 1845 to 1933 [16]. The recorded population data were actually the number of furs caught by the Hudson Bay Company, and an underlying assumption of the data analysis is that the number of caught furs was proportional to the actual population of hares and lynxes. The raw data are limited to 20 years from 1900 to 1920 and denote the yearly average population. The raw data are firstly interpolated with a cubic spline to create a smooth profile, . Canuto, D. Mazza, Introduction to population dynamics and resource exploitation January 29, 2021 34 as it allows population rate computation. N interpolation points have been assigned to each year, which implies a time unit 0.1 y s T . The initial time t has been set to zero, i.e. t . The data uncertainty due to measuring errors is not known. The time profiles of interpolated and raw data are reported in Figure 14. The population volume is expressed in kVol=1000 individuals. Figure 14. Spline interpolated and raw data (marked by crosses and circles). The first and second step estimated parameters are reported in Table 3. The first-step estimated value is in bracket. Table 3.
Hare and Lynx case parameter estimation No Parameter Symbol Unit Value Comments 1 Lynx initial population x kVol 4.04 (4.00) 2 Hare initial population x kVol 36.22 (30.00) 3 Lynx rate coefficient b kVol/y 0.92 (0.73) 4 Hare rate coefficient b kVol/y 0.48 (0.55) 5 Lynx average population x kVol 19.48 (20.86) 6 Hare average population x kVol 33.78 (34.35) 7 LV period LV P y 9.48 (9.88) . Canuto, D. Mazza, Introduction to population dynamics and resource exploitation January 29, 2021 35
8 Mean energy E kVol/y 0.38 (0.29) 9 Residual RMS, lynx x Vol/Vol 0.14 (0.62) Fractional error 10 Residual RMS, hare x Vol/Vol 0.13 (0.47) same 11 Error correlation (angle) (rad) -0.24 (1.81) The second-step estimated time profile is compared to the interpolated raw data in Figure 15, left. The estimation residuals, normalized by the population averages, are shown in Figure 15, right. Hare and lynx residuals are close to be each other orthogonal, as proved by the small correlation in Table 3, row 11. The estimated profiles show a significant deviation from the interpolated raw data during the low population period of both hares and lynxes from year 5 to 13. We can only imagine an external cause, like for instance, a higher proportion of caught furs with respect to the actual population than otherwise, or a more complex competition dynamics as proved in [20] and [21]. Figure 15. Left: estimated population time profile compared to raw data. Right: fractional residuals. The deviation the estimated profile from raw data is shown also by the plot, in Figure 16, of the energy E defined in (50). The deviation from year 5 to year 13 is evident and justifies the higher energy of the estimated model (black line). . Canuto, D. Mazza, Introduction to population dynamics and resource exploitation January 29, 2021 36 Figure 16.
Energy of the interpolated raw data (red) and the mean value (dashed red) compared with the energy of the estimated population profiles (black). The third estimation step provides the input perturbation k u t in (81) as the combination of low (LF)-and high-frequency (HF) components. The corresponding DT equation is as follows: uk uk k uk ukk uk k x i x i w i x xu i x i w i , (101) where the feedback signal , 0,1 jk w i j is proportional to the model error ˆ k k k x i x i x i , and ˆ k x i is the state variable of the DT LV equation to be explained below. The resulting PI (proportional and integrative) feedback driven by k x i has the form uk uk k k uk ukk uk k k x i x i L x i x xu i x i L x i , (102) where , k k L L are feedback gains fixing the eigenvalues (and the frequency BW) of the closed-loop system made by the DT LV equation and (102). As we will see below, the accuracy of the DT LV equation depends on the time unit T , which has been fixed smaller than the interpolation unit s T , i.e. 0.2 0.02 y s T T . The closed-loop BW c f has been fixed close to the interpolation Nyquist frequency, i.e. 5 Hz c f . The estimated DT input signals , 1, 2 k u i k , in [KVol/y] units, are shown in Figure 17. The residual normalized error / k k x x becomes negligible with respect to residuals in Figure 15, right The residula RMS is less than 0.001 Vol/Vol, a value to be compared with Table 3, rows 9 and 10. . Canuto, D. Mazza, Introduction to population dynamics and resource exploitation January 29, 2021 37 Figure 17.
Input signals (population rates) to LV equation to explain the residuals of the second step. The time profiles in Figure 17 appear to be dominated by 2-year period fluctuations, which are visible also in Figure 15, right, and to a less extent in Figure 15, left. Likely, such fluctuations may be fit by a short-time, time-varying dynamics: cycle explanation and fitting is beyond the paper aim. The correlation between lynx and hare profiles is about -0.4, a bit larger in magnitude than the correlation in Table 3, rows 11. The application of the two-species LV equation (7) to the lynx-hare competition has been judged too simple, from experimental data as in [20] and [21]. In essence, the hare is influenced by many predators other than lynx, and the lynx is primarily influenced by the snowshoe hare.
Application of LV equation: resource exploitation
A limiting form of the LV equation
In [2] a limiting form of the LV equations has been applied to resource exploitation and specifically to justify, under simple assumptions, and to generate the well-known bell-shaped Hubbert’s curve, which has been obtained from the logistic equation in Section 2.1. The key difference from the logistic equation is the assumption of two competing species, the resource, playing the role of the ‘prey’, and the human capital, employed in exploitation, playing the role of ‘predator’. Both resource and capital interact. The key assumption of the mathematical model is that the resource either cannot reproduce or the reproducing rate is negligible (non-renewable resource). Both capital and resource current amount (volume, stock) can be measured in currency units, energy units or mass units. Let us rewrite (42) with the external perturbation k u t , already introduced in (81): capital: , 0 , 0resource: , 0 , 0 x t b a x t x t u t x x x tx t b a x t x t u t x x x t , (103) . Canuto, D. Mazza, Introduction to population dynamics and resource exploitation January 29, 2021 38 where the same unit [Vol] applies to the pair of state variables. The assumption that resource cannot reproduce by itself implies b . (104) By assuming zero perturbation k u t , the equilibrium points of (103) and (104) are as follows:
0, 0 x x x . (105) In other terms, any resource value is an equilibrium point. Among such points, the critical point is defined by
0, / c c c x x b a x . It is left to the reader to find the forced equilibrium points under k k u t u . The local stability of the equilibrium points in (105) can be inferred from the perturbation equation of the state variable x x x and perturbation input , u u u , namely
11 1 2 2 021 2 c a b x xt t ta x x x u x x . (106) We remark that (106) is still nonlinear since perturbations from zero equilibrium cannot be negative. The zero equilibrium
0, 0 x x x is clearly stable since x . The stability of the nonzero equilibrium
0, 0 x x x splits in two cases: 1) Instability, when / 1 c x x . In the equality case, (106) becomes the series of two integrators. 2) Stability (not asymptotical stability), when / 1 c x x . Moreover, since x t , the initial resource perturbation can only decrease as shown by the free response exp , sgn anyexp t x t a t ax t x a x x a d . (107) The free response shape under k u t of the nonlinear equation (103) can be inferred by solving separately the two equations as follows
21 10 1 0 22 20 21 1 200 tt xx t x b t dt xx t x a x d x . (108) The resource volume x t always decreases to the limiting value x in (110) below, whereas the capital volume increases for x t x and then decreases to zero. As a result, the capital time profile x t follows a bell-shaped profile. The non-negative resource production p t , defined by . Canuto, D. Mazza, Introduction to population dynamics and resource exploitation January 29, 2021 39
21 1 2 p t x t a x t x t , (109) follows the bell shape too, since x t is a monotonic profile. Bell-shaped profiles can be considered as the limit of the LV periodic profiles, when the period LV P b b , which occurs in this case, since b has been assumed. The resource volume x t decreases to a limiting value x , since the integral of the capital x t in (108) is finite. From (108), x holds exp x a x d x . (110) A further variable of interest as pointed out in [2] is the ‘return on investment’ (ROI) / unit q t p t x t which from (109) is proportional to the resource volume
21 2 q t a x t , (111) and, therefore, is always decreasing. As a further property, let us denote with t the time when x t x . From (108) we find the equivalent of (58) to the present aperiodic solution, namely: t x d x b at . (112) Typical time profiles are in Figure 18, left. The counter-clockwise phase diagram production-capital is in Figure 18, right. The equation parameters are
310 20 1 12 21 , , , , 0.5 10 ,3.0,1.0,1.0, 0.5 x x b a a . (113) Figure 18. Left: Time profiles of capital stock, resource stock, production and ROI. Right: production-capital phase diagram. Since as shown by (118), the capital stock x can be made function of the production p , a phase diagram as in Figure 18, right, can be plotted. The diagram, as indicated by the initial arrow, is counter-clockwise, which implies that the internal area has negative sign. The internal area is the result of the integral . Canuto, D. Mazza, Introduction to population dynamics and resource exploitation January 29, 2021 40 p t p tp p V t x p p d x p dp , (114) where t is defined by p t p . To show that the integral in (114) is always negative let us define two other time instants: c t is defined by c x t x and corresponds to arg max c t t x t . (115) The second time instant p t is defined by arg max p t t p t . (116) The integral in (114) is always negative because p c t t as one can check from Figure 18, right. To prove it, let us set to zero the production time derivative which provides p p p p t x x px p t b x t x x t xa , (117) and proves that p c t t since p p t . In summary, model (103) with b , tells us that the capital stock is more efficiently employed when the resource stock is somewhat larger than the critical value c x , the limiting epoch being p t . The same meaning is given by the ROI q t . Though very simple and maybe questionable, equations (103) and (104) imply that the equilibrium capital x is zero as in (105). On the contrary, a nonzero equilibrium capital / x a b would appear for b as in the classical LV equation (42) or by assuming a nonzero mean perturbation u t u . This implies that renewable resources push a nonzero steady-state human capital! A first-order equation of the capital stock driven by production
If production p t x t is known instead of the resource stock x t , the state equation (103) becomes first-order and LTI as follows , 0 , 0/ x t b x t a p t u t x x x ta a a , (118) where a expresses how resources are transformed into capital. The solution of (118) is explicit as follows exp exp t x t b t x b t a p u d . (119) The corresponding sampled-data equation in the difference form, to be used for the parameter estimation, holds . Canuto, D. Mazza, Introduction to population dynamics and resource exploitation January 29, 2021 41 x i x i x i x i p i x xab T b Tb , (120) where u i has been assumed. The available measurements are given by , 0,..., 1 x i x i x i i Np i p i p i , (121) where , x i p i denote the measurement noise. Parameter estimation assumes (120) and (121), but exploits the orthogonality of the pair , / x x T to implement a well-conditioned Gramian matrix. To this end, the measurement equation is written as p i x i x i T ib bT Ta a b T , (122) where i p i x i x i T is the measurement error, and where the parameter conversion from , to , b a holds ln 1 / / , / b T T a b . (123) Also in presence of a zero-mean measurement error, the estimate , γ is biased, because of the noisy regression variables , 1 / x i x i T . Since the bias depends on the covariance of the measurement error, which is a second-order term in the expansion of γ , we can assume γ γ E . Given the covariance , , T P γγ E of the estimation error γ γ γ , the error covariance of the pair , b a can be computed from the following expansion in terms of the components of γ : / / / / b b b b Tb ba a a a aT (124) Table 4 and Figure 19 show estimation results of the model in (113). The measured variables are the capital stock and the production. The measurement error is a coloured noise with a driving noise pair with standard deviation , x p and cutoff frequency , x p f f as reported in Table 4. The last row of Table 4 reports the relative error RMS of the production estimation, defined by
211 00
Np N ii p i p iNp i . (125) . Canuto, D. Mazza, Introduction to population dynamics and resource exploitation January 29, 2021 42 Table 4.
Test case parameters No Parameter Symbol Unit Value Comments 1 Estimated interval H unit 10 2 Measurement time unit T unit 0.05 3 Low-frequency noise standard deviation , x p Vol (capital stock, production) 4 Cutoff frequency , x p f f
5 Capital obsolescence rate b fraction <0.02 6 Resource to capital transformation a fraction <0.03 7 Production estimation relative error (RMS) p fraction 0.09 Figure 19. True production, estimation and estimation error.
Application to experimental data: California gold rush
In [2], several fitting examples from experimental data have been proposed. We only consider the scarce data of the California gold rush accounted for in [18] and [19]. The gold rush (1848-1955) began on January 24, 1848, when flakes of gold were found in the American River by J.W. Marshall, a carpenter working to build a water-powered saw mill, the . Canuto, D. Mazza, Introduction to population dynamics and resource exploitation January 29, 2021 43
Sutter’s Mill, near Coloma, California. In the March 1848, the discovery was announced through the streets of San Francisco by the newspaper publisher and merchant S. Brennan, who holding a vial of gold, walked along the streets shouting ‘Gold, gold, gold from the American River’. The US-wide spread news attracted from 1848 to 1855 about 300000 people from the rest of United States and abroad. Production and capital stock raw data are available from [18] and [19]. They are rather scarce, sparse and uncertain. Gold production p i is measured in MUSdollar/y=10 USdollar/y . Capital stock x i is assumed to be proportional to the estimated number of gold prospectors [kVol]. Production and capital raw data have been spline interpolated with a time unit of T , beginning in 1843 and assuming zero production from 1843 to 1847. Raw and interpolated data are assumed to satisfy the DT equation (120) which has been converted into (122) for estimation purposes. No measurement error is known. The least-squares solution of (122) is reported in Table 5. The capital obsolescence was found to be , which means a reduction of about 50% per year. The resource to capital transformation was estimated to be about 1.6 prospector per 1000 US dollars, which means that in average each prospector produced about 625 US dollars per year. By assuming a price of about 0.75 US dollars each gram of gold, the average prospector produced less than 1 kg of gold per year. Table 5. California gold rush parameters No Parameter Symbol Unit Value Comments 1 Estimation interval H y 30 1943-1973 2 Interpolation time unit T y 0.1 3 Capital obsolescence rate b a Vol/USDollar p
10 USDollar/y 4.9 6 Relative estimation error / p p fraction 0.16 Notwithstanding the sparse and uncertain data, the interpolated capital profile (the blue curve in Figure 20, left) seems capable of fitting the interpolated gold production except during the peak production around the 1852 year and during the stabilized production set up after the gold . Canuto, D. Mazza, Introduction to population dynamics and resource exploitation January 29, 2021 44 rush, from about 1865. The peak production around the 1852 year may be the result of improvements of the gold-recovery techniques from placer mining to hydraulic mining. In fact, by assuming the smooth capital curve in Figure 20, left, a production peak higher than the estimated profile should imply that other capitals than prospectors were employed. Figure 20. Left: Capital and production raw data and estimated production. Right: production raw and interpolated data, estimated profile and estimation error. By assuming the measurement error variance to be of the order of p in Table 5, row 5, the fractional standard deviation of , b a was estimated from (124) to be less than 10%. As expected, the capital time profile lags the production profile of about two years, if measured between the peak values of the estimated profiles. Riferences [1]
D. Mazza, La curva sigmoide e il grafico di Hubbert (in Italian), February 2019, available from https://resources-and-reserves.blogspot.com/. [2]
U. Bardi and A. Lavacchi, A simple interpretation of Hubbert’s model of resource exploitation,
Energies,
T. R. Malthus,
An Essay on the Principle of Population As It Affects the Future Improvement of Society, with Remarks on the Speculations of Mr. Goodwin, M. Condorcet and Other Writers (1 ed.).
London, J. Johnson in St Paul's Church-yard, 1798. [4]
B. Gompertz, On the Nature of the Function Expressive of the Law of Human Mortality, and on a New Mode of Determining the Value of Life Contingencies . Philosophical Transactions of the Royal Society of London,
Vol. 115, 1825, pp. 513–585. [5]
J.F. Richards, A flexible growth function for empirical use,
Journal of Experimental Botany , Vol. 1, No. 10, 1959, pp. 290-310. . Canuto, D. Mazza, Introduction to population dynamics and resource exploitation January 29, 2021 45 [6]
P-F. Verhulst, Notice sur la loi que la population poursuit dans son accroissement,
Correspondance Mathématique et Physique,
Vol. 10, 1838, pp. 113–121. [7]
M. King Hubbert,
Nuclear Energy and the Fossil Fuels , Drilling and Production Practice , American Petroleum Institute & Shell Development Co. Publication No. 95, [8]
A.J. Lotka, Analytical note on certain rhythmic relations in organic systems,
Proc. Nat. Acad ., 6 (1920), 410-415. [9]
V. Volterra, Fluctuations in the abundance of a species considered mathematically,
Nature
118 (1926), 558-560. [10]
R. Arditi and L. R. Ginzburg, Coupling in predator-prey dynamics: ratio-dependence,
Journal of Theoretical Biology,
Vol. 139, 1988, pp. 311-326. [11]
J. Maynard Smith, Game Theory and The Evolution of Fighting , in On Evolution,
Edinburgh University Press, Cambridge, UK,1972. [12]
J. Maynard Smith,
Evolution and the Theory of Games,
Cambridge University Press, 1982. [13]
S-D. Shih, The period of a Lotka-Volterra system,
Taiwanese Journal of Mathematics, vol. 1, No.4, 1997, pp. 451-470. [15]
F. Rothe, Thermodynamics of the Volterra model, in J.P. Aubin, D. Saari and K. Sigmund (eds)
Dynamics of Macrosystems. Lecture Notes in Economics and Mathematical Systems , vol 257, pp. 49-62, 1985, Springer, Berlin, Heidelberg. [16]
J.M. Mahaffy, Math 636 – Mathematical modelling, Fall Semester 2010, Lotka-Volterra models, San Diego State University, 2000, available from https://jmahaffy.sdsu.edu/courses/f09/math636/lectures/lotka/qualde2.html. [17]
E. Canuto, C. Novara, L. Massotti, D. Carlucci, C. Perez Montenegro,
Spacecraft dynamics and control: the embedded model control approach.
Butterworth-Heinemann (Elsevier), 2018. [18]
M.J. Rohrbough,
Days of gold: the California gold rush and the American Nation,
University of California Press, 1997, Berkeley, CA. [19]
R.W. Paul,
California gold,
University of Nebraska Press, 1947, Lincoln, NE. . Canuto, D. Mazza, Introduction to population dynamics and resource exploitation January 29, 2021 46 [20]
N.C. Stenseth, W. Falck, O.N. Bjørnstadt and C.J. Krebs, Population regulation in snowshoe hare and Canadian lynx: asymmetric food web configurations between hare and lynx,
Proc. Natl. Acad. Sci. USA, Ecology,
Vol. 94, May 1997, pp. 5147-5152. [21]
R. Tyson, S. Haines and K.E. Hodges, Modeling the Canada Lynx and the snowshoe hare population: the role of specialist predators,
Theoretical Ecology,
Vol. 3, 2010, pp. 97-111. Appendix
Poincaré maps
Given an autonomous state equation of order n , 0 , dim t t n x f x x x x , (126) a way to find the properties of the trajectories in the phase space n X is to define a surface X of dimension n . Given an initial state x corresponding to a trajectory intersection, defined by f x n x where n x is the normal of the surface tangent Figure 21plane in x , find the trajectory intersections , ,... x x with the surface itself at times ... t t , such that sgn sgn , 0 k k k f x n x f x n x , (127) i.e. surface crossing occurs in the same direction (see Figure 21). The Poincaré or return map : P X X is the map that relates a generic initial state to the next return state , k k k k P t t x x . (128) The map does not exist if the trajectory does not return to the surface. There may exists trajectories which does not cross or just touch the surface. By excluding pathological trajectories, if the trajectory is close (orbit) we have k k k k P x x x x (129) and the time interval k k k T t t is the period of the orbit. In this case we may just restrict to the initial state x and to the first return state x x . A typical return surface is a n dimensional hyperplane defined by constant j j x t x x X, . (130) . Canuto, D. Mazza, Introduction to population dynamics and resource exploitation January 29, 2021 47 x x x x x x x x Figure 21.