Optimal Fishery Modelling Population Dynamics with Optimal Harvest Strategies

(1)

Optimal Fishery

Modelling Population Dynamics with Optimal Harvest Strategies

Degree Project in Engineering Physics, First Level, SA114X Department of Mathematics, division of Numerical Analysis

Mats Barkman Supervisor: Anders Szepessy

(2)

Abstract

Determining sustainable harvest strategies are of interest in or-der to avoid overexploitation of ecosystems. In this report the dynamics of fish populations affected by fishery is studied. The aim is to find an optimal catch policy under which fishery shall be performed in order to obtain maximum catch without the fish population facing extinction. By formulating this as an optimal control problem, Pontryagin’s principle can be used in order to find optimal harvest strategies. The numerical results show that for a spatially independent model of a single population, there is an optimal population size which gives maximum yield. For interde-pendent species, interference by fishery may change the qualitative behaviour of the dynamical system. The majority of the studied dynamical systems had a stationary optimal strategy.

(3)

Sammanfattning

För att undvika överexploatering krävs h˚allbara metoder för

ut-nyttjande av ekosystem. H¨ar studeras hur fiske p˚averkar

dy-namiken hos populationer av fisk. M˚alet med studien ¨ar att hitta

en optimal strategi s˚a att f˚angsten maximeras under villkoret att

fiskpopulationen ej utrotas. Med hj¨alp av optimal styrteori och

Pontryagin’s princip, kan s˚adana strategier erh˚allas. Numerisk

analys av det optimala styrproblemet visade att det f¨or oberoende

fiskarter existerar en optimal populationsstorlek. F¨or beroende

fiskarter, kunde dynamiken av ekosystemet kvalitativt f¨or¨andras.

I majoriteten av de studerade exemplen var den optimala strategin tidsoberoende.

(4)

1 Introduction

Commercial fishery relies on the ability of the fish to repopulate. Since interference with ecosystems changes its dynamics, uncontrolled harvesting may result in overexploitation. Therefore it is of interest to find sustainable strategies yielding large catch.

The dynamics of the fish population may be dependent of various ecological variables, such as the population size of natural predators, climate change and access of food. In order to determine optimal fishery policies, a suitable model for the ecosystem must be obtained.

Apart from complex dynamical system of the ecosystem, defining optimal fishery is a complex problem in bioeconomics. In order to obtain a suitable definition, both ecological and economical viewpoints need to be taken into consideration. In this study optimal result is considered as the outcome yielding the highest mean value of catch over a given time period, without the fish population facing extinction. The simplicity of this definition makes mathematical analysis possible but does not necessary imply trivial results. A

simi-lar definition of optimality has been used to justify marine reserves1_{, studying a stationary}

spatially dependent model [Neu03]. In this report a non-stationary spatially independent model is considered. The optimal harvest strategy will be analysed for various ecosystems, containing multiple interdependent fish species.

2 Background

2.1 Optimal Control Problem

The task of finding the optimal fishery strategy can be formulated as an optimal control

problem. Consider a variable X(t) ∈ Rn _{that satisfies the ordinary differential equation}

˙

X(t) = dX(t)

dt = f (X(t), α(t)), 0 < t < T, X(0) = X0, (1)

where t is time, the control function α ∈ A and where A = {α : [0, T ] → A} is the set of

allowed control functions, the flux f : Rn_{× A → R}n _{and the initial condition X}

0 ∈ Rn.

Let X = {X : [0, T ] → Rn_{} be the set of possible X. The optimal control problem is to}

find an extreme value of the functional J : X × A → R J =

Z T

0

h(X(t), α(t))dt + g(X(T )), (2)

with given cost functions h : Rn_{× [0, T ] → R and g : R}n _{→ R and flux f and initial}

condition X0. That is, find a control for all times that optimizes Equation (2) with

Equation (1) as a constraint.

2.2 Symplectic Maps and Hamiltonian System

Differential equations of motion are imperative of understanding the time evolution of dynamical systems. Lagrange formulated a method of finding these differential equations.

Introduce general coordinates q = (q1, . . . , qn)T and its time derivatives ˙q which together

describe the state of a system with n degrees of freedom. The kinetic energy T (q, ˙q) and

1

(6)

potential energy V (q, ˙q) can be expressed in terms of a given state (q, ˙q). By introducing the Lagrangian

L(q, ˙q) := T (q, ˙q) − V (q, ˙q),

using calculus of variation and Hamilton’s principle give Euler-Lagrange’s equations [Apa12] d dt ∂L ∂ ˙qi = ∂L ∂qi , i = 1, . . . , n. 2.2.1 Hamiltonian Systems

By introducing the Hamiltonian

H(p, q) := pTq − L(q, ˙q),

where p is the momenta

pi :=

∂L ∂ ˙qi

, i = 1, . . . , n,

Euler-Lagrange’s equations can be reduced to a system of first order. Euler-Lagrange’s equations are equivalent to

( ˙ qi(t) = ∂piH(p(t), q(t)) ˙ pi(t) = −∂qiH(p(t), q(t)) , i = 1, . . . n.

The Hamiltonian is interpreted as the energy of the system. By introducing z = (p, q)T

the system of differential equations can be described as

˙z(t) = S−1· ∇H(z(t)), S = 0 I

−I 0

, (3)

where I is the n × n identity matrix [EHW02].

2.2.2 Symplectic Maps

A property of Hamiltonian systems is that its flow is a symplectic map. The definition is presented below.

Definition: A linear mapping A is symplectic if

ω(Aξ, Aη) = ω(ξ, η), ∀ξ, η ∈ R2n,

where ω(ξ, η) is the sum of the oriented areas of the projection of P on the planes (pi, qi)

and P = {tξ + sη|0 ≤ t, s ≤ 1}. The statement is equivalent to ω(ξ, η) = ξT_{Sη, where S}

is the symplectic matrix presented in Equation (3).

Similarly the definition of a symplectic differentiable map is the following.

Definition: Let U ⊂ R2n _{be an open subset. A differentiable map g : U → R}2n _is

called symplectic if the Jacobian of g is symplectic everywhere.

One property of symplectic differentiable maps is that the area of a manifold is preserved under the flow of the differentiable map. Since Hamiltonian system are symplectic, the area of a section in phase space in preserved under a Hamiltonian flow [EHW02].

(7)

2.3 Hamilton-Jacobi-Bellman Equation and Pontryagin’s Principle

The optimal control problem stated in Section 2.1 can reformulated using dynamic

pro-gramming. Introduce the value function u : Rn_{× [0, T ] → R}

u(x, t) := sup α∈A Z T t h(X(t), α(t))dt + g(X(T )) , X(t) = x.

Since u(·, T ) = g(X(T )) is known, the dynamic programming approach is to iterate

back-wards in time in a manner such that the global extrema of u(X0, 0) is obtained by finding

the optimal path [CMS+_{10]. This leads to the Hamilton-Jacobi-Bellman partial differential}

equation (HJB equation)

∂tu(x, t) + H(∂xu(x, t), x) = 0, (4)

where H : Rn_{× R}n_{→ R is the Hamiltonian}

H := max

α∈A(∂xu(x, t)f (x, α) + h(x, α)).

The characteristics of the HJB equation is a Hamiltonian system. By introducing λ(t) :=

∂xu(x, t) the HJB equation gives

( ˙_{X(t) = ∂}_λ_{H(X(t), λ(t))} _{X(0) = X}₀

˙λ(t) = −∂XH(X(t), λ(t)) λ(T ) = g0(X(T ))

, (5)

which is Pontryagin’s principle. However the Hamiltonian is in general not differentiable.

In order to obtain well posedness it is necessary that (∂λH, −∂XH) are Lipschitz

contin-uous2_{. In addition, Pontryagin’s principle does not guarantee that a global maximum is}

found [CMS+_10].

3 Dynamics of Fish Populations

3.1 The Logistic Equation

Modelling biological systems can be very complex, since in an ecosystem, many natural parameters may affect the population growth. One model introduced 1838 by Verhulst [GQ38] is the logistic equation for population growth

˙ X(t) = F (X(t)) = rX(t) 1 −X(t) K . (6)

Here X is the population size of the fish species. The intrinsic growth rate r > 0 and the carrying capacity K > 0 are natural constants for each species. If X is larger than K, the population size will decrease and if X is less than K, the population size will increase. The dynamical system has one stable fixed point corresponding to X = K.

The solution to Equation (6) can be found analytically as

X(t) = X0

ert

K + X0(ert− 1)

,

2

Note that Lipschitz continuous flux is needed in any differential equation in order to provide well posedness.

(8)

and is presented in the Figure 1 for different initial conditions X0. Time t 0 1 2 3 4 5 6 7 8 9 10 F is h P op u lat ion X 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

Figure 1: Solutions to Equation (6) with r = K = 1 with different initial conditions

X0 ∈ [0, 2]. The red lines show the fixed points X = 0 and X = K.

3.2 Lotka-Volterra Model of Competition

The dynamics of populations was expanded further, taking into account interaction be-tween species, by Volterra, Lotka and Gause [GFG35]. Firstly, consider two interacting

species with population sizes X1, X2, intrinsic growth rates r1, r2 and carrying capacities

K1, K2. The dynamics is proposed to be governed by

˙ X1(t) = r1X1(t) 1 −X1(t) + w12X2(t) K1 , ˙ X2(t) = r2X2(t) 1 −X2(t) + w21X1(t) K2 ,

where w12 and w21 are positive weights. The species compete for the same natural

re-sources, for example food and habitat. Therefore each interaction between the species decreases the growth of the populations. However, the two species may be affected differ-ently. Therefore the qualitative behaviour of the dynamical system may vary [Cla76].

(9)

X1 X2

(a) Competitive coexistence.

X1 X2 (b) Competitive exclusion. X1 X2 (c) Competitive exclusion.

Figure 2: Illustration of the qualitative behaviour of the dynamics of the Lotka-Volterra

model of competition. The dashed line corresponds to ˙X1 = 0 and the dot-dashed line to

˙

X2 = 0. The black arrows represent the vector field and the blue lines some trajectories.

If the interaction between the species is weak, the dynamical system has a stable fixed

point (where the lines corresponding to ˙X1 = 0 and ˙X2 = 0 intersect in Figure 2a). That

is, the dynamics allows for a competitive coexistence.

However, competitive exclusion occurs when the interaction is stronger. If both species are noticeably affected by each other, the dynamical system has two stable fixed points, corresponding to one or the other species facing extinction, as presented in Figure 2b. If one of the species is significantly more affected than the other, there is only one stable fixed point, as in Figure 2c.

3.3 Generalised Model of Competition

In Section 3.2 the model was restricted to two competing species. For n competing species

Xi, where i = 1, 2, . . . , n, the generalisation of the two dimensional Lotka-Volterra

dynam-ical system is ˙ Xi(t) = Fi(X(t)) = riXi(t)  1 − 1 Ki n X j=1 wijXj(t)  , i = 1, 2, . . . n, (7)

(10)

where wij is related to the affect the species j has on species i. The weights can be

collectively represented using the weight matrix W

W =      1 w12 · · · w1n w21 1 w2n .. . . .. wn1 1      .

For convenient notation, the vector X = (X1, X2, . . . , Xn)T is introduced in order two

write the dynamical system as

˙

X(t) = F (X(t)), (8)

where F : Rn_{→ R}n_{and F}

i(X) is given in Equation (7). The dynamical system can exhibit

various properties depending on the parameters and the value of n. For higher dimensions the behaviour of the dynamical can be complex. For example,

W =     1 1.09 1.52 0 0 1 0.44 1.36 2.33 0 1 0.47 1.21 0.51 0.35 1     ,     r1 r2 r3 r4     =     1 0.72 1.53 1.27     ,     K1 K2 K3 K4     =     1 1 1 1     , (9)

which is one of the lowest dimensional dynamical systems that exhibit chaotic behaviour

[VWA+_{06]. Its numerical solution is presented in Figure 3.}

Figure 3: Four-dimensional dynamical system in Equation (9) with chaotic behaviour.

The initial condition is X(0) = X0 = (0.2, 0.2, 0.2, 0, 2)T.

3.4 Fishery

If fishery by man is accounted for in the dynamics of the fish population, Equation (8) is modified accordingly

˙

X(t) = f (X(t), E(t)) = F (X(t)) − C(X(t), E(t)), (10)

where E = (E1, E2, . . . , En)T ∈ A is the control and C : Rn× A → Rn is the catch. The

value of C(X(t), E(t)) is the rate at which fishery is performed at time t. One suitable model for the catch is the bilinear relationship

(11)

introduced by Schaeffer [Sch54]. The control Ei can be considered as the effort put into

fishery of species i. The effort can for example be the time spent fishing, the amount of active fishing vessels or the efficiency of the fishing techniques. The model infers that an increase in either effort or population size yields a larger catch.

Consider first the case when n = 1 and the effort E is constant. In this case, the so-lution to Equation (10) can be found analytically

X(t) =    X0 ke krt k+X0/K(ekrt−1), k = 1 − E r 6= 0 rt K + 1 X0 −1 , k = 0 . (11)

Figure 4 illustrates the stationary solutions to Equation (10) and their stability when considering a single fish species.

X

F(X) C(X, E)

Figure 4: Note that X∗ _{= K 1 −}E

r marked as a dashed line is a stable fixed point of

the system ˙X(t) = F (X(t), E) − C(X(t), E) when E is constant, assuming that E < r.

If E > r there is no stable fixed point, since it would imply that C(X, E) > F (X) for all X 6= 0.

Let C(X∗_{, E) be the equilibrium catch rate,}

C(X∗, E) = EK 1 −E r , E < r.

It is easily seen that the maximum equilibrium catch rate is achieved when E = E∗ = r/2

and the equilibrium population size X∗ _{= K/2. Geometrically this corresponds to the line}

C(X, E) intersecting F (X) at its maximum point. The maximum equilibrium catch rate is referred to as the Maximum Sustainable Yield (MSY) [Cla76].

(12)

4 Analysis of Optimal Control Problem

4.1 Formulation of Optimal Control Problem

Primarily the functions for the optimal control problem in Equations (1)-(2) have to be specified. The flux f is given by Equation (10). The aim is to find an effort such that the functional J

J =

Z T

0

h(X, E)dt + g(X(T )),

is maximized. Since the mean value of the catch rate is to be maximized h is chosen as

h(X, E) = 1 T n X i=1 Ci(X(t), E(t)).

The cost function g must be specified. Since a large population size at the time horizon T is beneficial, it is suitable to choose g as an increasing function. However, exceeding the

initial population size X0 shall not be noticeably more rewarding than reaching the initial

population size. Consider first the case of only one species. The cost function g is then chosen as g(X) = D arctan 10X X0 ,

where D ∈ R is some constant. A reasonable value of D is derived from that the reward from g(X(T )) and from the fishery should be of the same magnitude. Assume that X(t) =

X0 for all times and the equation

1 T Z T 0 C(X0, E0)dt = D arctan 10X0 X0 | {z } ≈π/2 ≈ π 2D,

with E0 such that F (X0) − C(X0, E0) = 0. That is,

E0 = r 1 −X0 K . This gives a plausible value of the constant D.

D = 2r π X0 1 −X0 K .

In order to generalise for n competing species, the cost function g is chosen to g(X) = n X i=1 2ri π 1 −X0,i Ki arctan 10 Xi X0,i ,

where X0,i is the initial population size of species i.

However, the condition X(T ) = X0 could possibly replace the cost function g. The

mathematical interpretation of this is that g is an impulse function at X = X0, making it

infinitely rewarding terminating at X0.

Reaching the initial population size at the time horizon might be impossible for some dynamical systems. Not specifying g is possible in the one dimensional case, as long as

(13)

some conditions are met. The effort E can only take values between 0 and the maximum

effort Emax. If the fish initial population size is larger than the carrying capacity of the

species, it is impossible to reach the initial size at the time horizon, since ˙X(t) < 0 for

all times even though no fishery is performed. On the other hand, if the maximum effort can be performed without the fish population decreasing in size, the initial population size cannot be reached at the time horizon. In conclusion, the optimal control problem can be

formulated using the functional J1 or the functional J2 where

J1= 1 T Z T 0 C(X(t), E(t))dt, X(T ) = X0, (12) and J2= n X i=1 1 T Z T 0 Ci(X(t), E(t))dt + 2ri π 1 −X0,i Ki arctan 10 Xi X0,i . (13)

However a solution satisfying X(T ) = X0is not possible for all dynamical systems. Finding

necessary conditions that allow for such a solution for systems with arbitrary number of populations is difficult. Therefore Equation (12) is only studied in the case of one fish species.

4.1.1 Defining New Variables

By introducing new dimensionless variables the amount of constants in the optimal control problem can be reduced. The reduction is only beneficial when considering one fish species. Therefore the following rescaling is only performed:

x = X K, x0 = X0 K , t 0 = rt, T0 = rT, α = E r, αmax= Emax r .

This modifies the constraint given in Equation (10) accordingly

dx(t0₎

dt0 = x(t

0_{)(1 − x(t}0_{)) − α(t}0_)x(t0_).

The functional J1 becomes

rK 1 T0 Z T0 0 x(t0)α(t0)dt0 ! .

Since the term rK is constant, without loss of generality J1 is chosen as

J1= 1 T0 Z T0 0 x(t0)α(t0)dt0.

Performing the rescaling of the variables gives that only two constants, the maximum value of the control and the time horizon are needed to describe the problem. However, in order to keep notation not too prolix Equations (14)-(15) will be used in the optimal control problem in the case of one fish species. The functional is

J1 = 1 T Z T 0 x(t)α(t)dt, (14)

and the constraint is

˙x(t) = f (x(t), α(t)) = x(t)(1 − x(t)) − α(t)x(t), x(0) = x(T ) = x0, (15)

where α(t) ∈ [0, αmax]. In the case of competing species, the optimal control problem is to

find the efforts Ei(t) ∈ [0, Emax,i] that maximizes Equation (13) with Equation (10) and

(14)

4.2 Application of Pontryagin’s Principle

The Pontryagin’s principle (Equation (5)) will now be applied to derive the Hamiltonian dynamical system of the fish species X and co-states λ. Provided the Hamiltonian system, the optimization problem stated in Section 4 can be solved. Firstly the derivation of the Hamiltonian system using Equations (14) and (15) will be shown. The Hamiltonian H is

H(x, λ) = max

α∈A(λf (x, α) + h(x, α)) = maxα∈A

λx(1 − x) − αx 1

T − λ

.

The linearity in the control α results in the optimal control α?(λ)

α?(λ) =

(

αmax, λT < 1

0, λT > 1 (16)

Since α? is discontinuous, the Hamiltonian H will not be C1. The flux in the Hamiltonian

system will therefore not be Lipschitz continuous and well posedness is not obtained. By adjusting the cost function h using the continuation principle

h(x, α) = αx1 − δα T , δ > 0, gives α?(λ) as α?(λ) =      αmax, λT < 1 − 2δαmax 1−λT 2δ , 1 − 2δαmax< λT < 1 0, λT > 1 ,

which is Lipschitz continuous and well posedness is obtained. The process of converting the ill-posed problem into a well posed is called regularisation. As δ approaches zero, the initial optimal control is obtained.

λT -3 -2 -1 0 1 2 3 α⋆ ( λ) / αm a x 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 δ=1 δ=0.5 δ=0.25 δ=0.125 δ=0.0625

Figure 5: As δ approaches zero, the regularized optimal control α?(λ) converges to the

initial optimal control.

Given the optimal control, the Hamiltonian system describing the dynamics is ( ˙x(t) = x(t)(1 − x(t)) − α?(λ(t))x(t), ˙λ(t) = λ(t)(2x(t) − 1) + α?(λ(t)) λ(t) − _T1 + δα?(λ(t)) 2 T , x(0) = x(T ) = x0. (17)

(15)

For n competing species, the functional J2 in Equation (13) and the flux in Equation (10)

are used. Similarly as in the one dimensional case, the continuation principle is used to obtain well posedness redefining the cost function h as

h(X, E) = n X i=1 EiXi 1 − δEi T , δ > 0.

This gives the optimal control E?(λ) = (E?,1(λ1), . . . , E?,n(λn))T where

E?,i(λi) =     

Emax,i, λiT < 1 − 2δEmax,i 1−λiT 2δ , 1 − 2δEmax,i < λiT < 1 0, 1 < λiT . The Hamiltonian H is H(X, λ) = λTf (X, E?(λ)) + h(X, E?(λ)) = = n X j=1 λj(Fj(X) − E?,j(λj)Xj) + E?,j(λj)Xj 1 − δE?,j(λj) T

The Hamiltonian system of X and λ is consequentially determined by ˙ Xi(t) = ∂λiH(X(t), λ(t)) = Fi(X(t)) − Ci(X(t), E?(λ(t)) = = riXi(t)  1 − 1 Ki n X j=1 wijXj(t)  − E?,i(λi(t))Xi(t), i = 1, 2, . . . n, and3 ˙ λi(t) = −∂XiH(X(t), λ(t)) = . . . = = −riλi(t)  1 − 1 Ki n X j=1 wijXj(t)  + n X j=1 rjwjiλj(t) Xj(t) Kj + +E?,i(λi(t)) λi(t) − 1 T + δE?,i(λi(t)) 2 T , i = 1, 2, . . . n,

with boundary conditions

Xi(0) = X0,i, λi(T ) = ∂Xig(X(T )), i = 1, 2, . . . n.

5 Numerical Algorithms

5.1 Symplectic Euler Method

The continuous time evolution of the characteristics of the Hamiltonian is transformed into a time discrete evolution. That is, the solution is determined by the values of the fish populations and the co-states at N points in time. These points are collected in the time vector t

t = (0, ∆t, 2∆t, . . . , (N − 1)∆t)T, ∆t = T

N − 1.

3

(16)

The values of the fish population sizes and co-states at these points in time are represented by the vectors X and λ

X = (X₀T, X₁T, . . . , xT_{N −1})T, Xi ≈ X(ti)

λ = (λT₀, λT₁, . . . , λT_{N −1})T, λi≈ λ(ti).

Note that Xi ∈ Rn here is a vector representing the population sizes of all fish species

at time ti and shall not be misinterpreted as the population size of species i. Similarly

λi ∈ Rn represents all the co-states at time ti.

In order to preserve the symplectic structure of the Hamiltonian system, the symplec-tic Euler method [EHW02] is used to approximate the derivatives.

(

Xi+1 = Xi+ ∆t · ∇λH(Xi, λi+1)

λi+1 = λi− ∆t · ∇XH(Xi, λi+1)

, i = 0, 1, . . . , N − 2. (18)

5.2 Iterative Procedure of Equation Solving

In order to find X and λ they must satisfy the iterated map in Equation (18). Therefore

an equation M (X, λ) = 0 can be constructed, where M : RnN_{× R}nN _{→ R}2nN_{. Its solution}

X, λ satisfies Equation (18) and the boundary conditions. The function M is                            M1 = X0− X0 (M2, M3, . . . , M2+n)T = X1− X0− ∆t · ∇λH(X0, λ1) (M3+n, M3+n+1, . . . , M3+2n)T = λ1− λ0+ ∆t · ∇XH(X0, λ1) .. . (M2nN −3n+1, . . . , M2nN −2n)T = XN −1− XN −2− ∆t · ∇λH(XN −2, λN −1) (M2nN −2n+1, . . . , M2nN −n)T = λN −1− λN −2+ ∆t · ∇XH(XN −2, λN −1) (M2nN −n+1, . . . , M2nN)T = λN −1− ∇g(XN −1) or XN −1− X0 ,

where the last row depends on whether the condition at the time horizon is specified by cost function g or that the initial population size is reached.

Once the function M has been defined, the equation M (X, λ) = 0 is solved numerically using the MATLAB method fsolve from the Optimization Toolbox [MATa].

5.2.1 Regularisation Iteration

Since the original optimal control is obtained as δ → 0, the equation M (X, λ) = 0 should

be solved for some sufficiently small δmin. In order to obtain the solution for small values

of δ, the equation M (X, λ) = 0 is solved for some initial δ0. The solution is used as initial

guess for solving M (X, λ) = 0 for another δ1 < δ0. The procedure is iterated until the

desired limit δmin is reached. The process of finding the solution to the optimal control

(17)

Solve equation M (X, λ) = 0 Initial guess X0_{, λ}0 δi+1= δi/2 else Solution X, λ found if δ0 δi> δmin

Use Xi_{, λ}i_{as initial guess}

Figure 6: Flowchart of the procedure to find the solution to the optimal control problem using regularisation iteration.

5.2.2 Time Horizon Iteration

An initial guess is always needed in the iteration process. Obtaining suitable initial guesses is important for the solution to converge. For small time horizons T , it is reasonable that deviations from the initial population sizes are small. Hence a suitable choice as initial

guess for small time horizons should be close to X0. However, for larger time horizons,

the behaviour may deviate from the initial condition. Therefore finding a suitable initial guess is more difficult for larger time horizons.

In order to obtain a initial guess for large time horizons, the solution to the equation M (X, λ) = 0 for some smaller time horizon may be used. Therefore the following itera-tion process can be used:

1. Solve M (X, λ) = 0 for some small time horizon T0, for which an initial guess is easily

obtained.

2. Use the solution as initial guess to the equation M (X, λ) = 0 for some larger time

horizon Tstart+ ∆T where ∆T > 0 is some constant.

3. Iterate the procedure of adding ∆T to the current time horizon and equation solving until the desired time horizon T is reached.

The time horizon iteration is depicted in Figure 7.

Solve equation M (X, λ) = 0 Initial guess X0_{, λ}0 else Solution X, λ found if Ti+1= Ti+ ∆T T0 Ti< T

Use Xi_{, λ}i_{as initial guess}

Figure 7: Flowchart of the procedure to find the solution to the optimal control problem using time horizon iteration.

(18)

For some systems, the regularisation iteration is sufficient in order to obtain a numerical solution. For more complex system with large time horizons, the regularisation iteration can be combined with the time horizon iteration, as depicted in Figure 8.

Regularisation iteration for some T0 and δmin,1 Regularisation iteration for Time horizon iteration with step ∆T until T is reached some T and δmin,2

Figure 8: Iterative procedure combining regularisation and time horizon iteration.

6 Results

6.1 Single Fish Species

Firstly the case of a single species is studied. The numerical result is presented in Figure

9 for some different combinations of αmax and x0 using the functional and the boundary

condition in Equation (12). For a single fish species, the numerical iteration process only involves regularisation iteration.

t 0 1 2 3 4 5 6 7 8 9 10 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 x α (a) αmax = 1, x0= 0.7 t 0 1 2 3 4 5 6 7 8 9 10 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 x α (b) αmax = 1, x0= 0.2 t 0 1 2 3 4 5 6 7 8 9 10 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 x α (c) αmax = 0.4, x0= 0.7 t 0 1 2 3 4 5 6 7 8 9 10 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 x α (d) αmax = 2, x0= 0.2

Figure 9: Numerical solution with the number of point N = 150, T = 10 and δmin = 2−10

using Equation (12).

(19)

Figures 9a, 9b and 9d the optimal strategy is stationary corresponding to α = 1/2 and x = 1/2. This result coincides with the analytical solution corresponding to maximum sustainable yield discussed in Section 3.4. However, the optimal state can not be reached

if αmax< 1/2 as presented in Figure 9c.

In conclusion, the effort α is

α(t) = ( ϕ+(t), x0 < 1/2 ϕ−(t), x0 > 1/2 , where ϕ+(t) =      0, 0 < t < tr1 1/2, tr1 < t < T − td1 αmax, T − td1< t < T , and ϕ−(t) =      αmax, 0 < t < tr2 1/2, tr2 < t < T − td2 0, T − td2< t < T ,

provided that αmax > 1/2. The time tr is the time needed to reach the equilibrium

population size and td the time needed to go from the equilibrium size to the initial

population size. Since the effort is piecewise constant, Equation (11) can be used to find

the times tr and td analytically.

tr1 = ln 1 − x0 x0 , td1= (₁ kln x0 k−x0 · (2k − 1) , k = 1 − αmax6= 0 1 x0 − 2, αmax= 1 , tr2= (₁ kln k−x0 x0 · 1 2k−1 , k = 1 − αmax6= 0 2 −_x1 0, αmax= 1 , td2= ln x0 1 − x0 .

Using these expressions, the critical time horizon Tc= tr+tdcan be found. The equilibrium

population size is possible to reach if T > Tc.

Tc=    ±ln 1−x0 x0 + 1 x0 − 2 , αmax= 1 ±ln 1−x0 x0 + 1 kln x0 k−x0 · (2k − 1) , k = 1 − αmax6= 0 , +x0<1/2 −x0>1/2.

(20)

t 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 x α tr T₋td (a) T = Tc+ 1 t 0 0.5 1 1.5 2 2.5 3 3.5 4 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 x α tr T₋td (b) T = Tc t 0 0.5 1 1.5 2 2.5 3 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 x α (c) T = Tc− 1

Figure 10: Numerical solution for different T with N = 150, αmax = 1, x0 = 0.2 and

δmin= 2−13.

The functional in Equation (13) is also used in the case of one fish species. The result is presented in Figure 11. t 0 1 2 3 4 5 6 7 8 9 10 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 X 1 E 1 (a) αmax = 1, x0= 0.7 t 0 1 2 3 4 5 6 7 8 9 10 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 X 1 E 1 (b) αmax = 1, x0= 0.2

Figure 11: Numerical solution with the number of point N = 150, T = 10 and δmin = 2−10

using the functional in Equation (13).

(21)

com-pared. The only difference in the functionals is the cost function related to the fish population size at the time horizon. Hence difference is only expected in the result close to the time horizon. The results in Figure 9b and Figure 11b coincide fairly well. However in Figure 9a and Figure 11a there is an evident difference. The effort goes to zero close

to the time horizon in Figure 9a in order to meet the condition x(T ) = x0. In Figure 11a

the effort goes to its maximum value, which results in a drop in the population size. This difference may be explained by the reward from fishing being noticeably larger than the reward from cost function g.

6.2 Two Competing Species

The numerical algorithm is also applied to some dynamical systems involving two interact-ing fish species. Three examples will be considered, one of each class presented in Figure 2. For all examples, both regularisation and time horizon iteration is performed (see Figure 8). In Appendix B, the analytical stationary population sizes yielding maximum catch are derived and compared to the numerical result.

6.2.1 Competitive Coexistence Example

The dynamical system is determined by

W = 1 0.5 0.4 1 , r1 r2 =2 1 , K1 K2 =1 1 , (19)

and Equation (7). If no fishery is performed, this is an example of competitive coexistence as in Figure 2a. The optimal fishery strategy is determined under the conditions

Emax,1 Emax,2 =1 1 , X0,1 X0,2 =0.5 0.5 , T = 10. (20)

The numerical calculation was conducted for two sets of parameters. Both calculations shared the parameters,

N = 50, T0 = 5, ∆T = 1, δmin,2 = 2−10. (21)

For the first set δmin,1 = 2−4 was used and for the other δmin,1 = 2−7 was used. The

(22)

t 0 1 2 3 4 5 6 7 8 9 10 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 X₁ X₂

(a) Time evolution if no fishery is performed.

t 0 1 2 3 4 5 6 7 8 9 10 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 X₁ E₁ X₂ E₂

(b) Optimal fishery with δmin,1= 2−4.

t 0 1 2 3 4 5 6 7 8 9 10 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 X₁ E₁ X₂ E₂

(c) Optimal fishery with δmin,1= 2−7.

Figure 12: Numerical study of the dynamical system determined by Equations (19)-(21) with functional in Equation (13).

The numerical algorithm for both values of δmin,1 converges to same solution (see

Figures 12b and 12c). The efforts and fish population sizes stabilize and reach an optimal stationary solution. In Figure 12a the stationary solution corresponds to a larger value of

X2 than X1. On the other hand, if fishery is performed as in Figures 12b and 12c, then

X1 is larger than X2.

6.2.2 Competitive Exclusion Example 1

Consider the dynamical system determined by

W =1 3 3 1 , r1 r2 =1 1 , K1 K2 =1 1 . (22)

If no fishery is performed, the dynamical system has two stable fixed points, corresponding to one or the other fish species facing extinction, as in Figure 2b. The case with

Emax,1 Emax,2 =1 1 , X0,1 X0,2 =0.45 0.5 , T = 10, (23)

is studied. The parameters used in the algorithms are the same as in the example of competitive coexistence in Section 6.2.1. The numerical result is presented in Figures 13b and 13c.

(23)

t 0 1 2 3 4 5 6 7 8 9 10 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 X₁ X₂

t 0 1 2 3 4 5 6 7 8 9 10 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 X₁ E₁ X₂ E₂

Figure 13: Numerical study of the dynamical system determined by Equations (21), (22) and (23) with functional in Equation (13).

In Figure 13b the optimal fishery policy is oscillatory, and not stationary. Fishery is performed in such a manner that only one fish species is fished at the time, allowing for the other one to regain in population size. The reward is evaluated numerically using MATLAB function trapz [MATb] to approximately 0.57 using the functional in Equation (13).

The optimal strategy in Figure 13c is neither stationary nor oscillatory. In this case the reward is numerically calculated to approximately 0.60. This indicates that the fish-ery strategy presented in Figure 13c is more optimal than the strategy in Figure 13b. However, the difference in reward is small. The numerical approximations in the algo-rithms may cause errors larger than the difference. The difference in reward from the two strategies is not a contradiction. As mentioned in Section 2.3, Pontryagin’s principle does not guarantee a global maximum. Hence both strategies can give local maximums of the reward.

Finally the same example is considered with an extended time horizon T = 15. The numerically calculated optimal strategies are presented in Figure 14.

(24)

t 0 5 10 15 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 X₁ E₁ X₂ E₂

(a) Optimal fishery with δmin,1 = 2−4.

t 0 5 10 15 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 X₁ E₁ X₂ E₂

Figure 14: Numerical study of the dynamical system determined by Equations (21), (22) and (23) but with T = 15 and functional in Equation (13).

The oscillatory behaviour in Figure 13b is carried on to the solution for the larger time horizon (see Figure 14a). However, the frequency of the oscillations has decreased. There-fore it cannot be concluded that there is a certain oscillation frequency yielding optimal catch for larger time horizons. Numerical approximation of the reward related to Figure 14a is calculated to 0.57.

The qualitative behaviour of the strategies in Figures 13c and 14b are very similar. The optimal strategy in Figure 14b is retrieved by rescaling the strategy in Figure 13c such that it fits the larger time horizon. The reward related to Figure 14b is calculated numerically to approximately 0.60.

6.2.3 Competitive Exclusion Example 2

Consider the dynamical system determined by the following values of the weights, carrying capacities and intrinsic growth rates

W = 1 1 0.5 1 , r1 r2 =1 1 , K1 K2 =1 1 . (24)

The dynamical system without fishery is an example of competitive exclusion (see Figure

2c). In this case the fish population X1 faces extinction and X2 survives as in Figure 15a.

Consider fishery with Emax,1 Emax,2 =1 1 , X0,1 X0,2 =0.5 0.5 , T = 10. (25)

The parameters used in the numerical algorithms are same parameters as in Section 6.2.1. The numerical result is presented in Figures 15b and 15c.

(25)

t 0 1 2 3 4 5 6 7 8 9 10 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 X₁ X₂

0 1 2 3 4 5 6 7 8 9 10 t 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 X₁ E₁ X₂ E₂

The numerical result presented in Figures 15b and 15c coincide. Similarly to the result in Section 6.2.1, the efforts and fish population sizes reach an optimal stationary solution. That means that maximum catch rate is achieved when fishery is performed in a manner such that both fish species survive. Hence the optimal strategy preserves biological diversity.

6.3 Four Competing Species

Numerical study of optimal fishery for the chaotic ecosystem of four interdependent fish species in Equation (9) is also performed. The maximal efforts, initial population sizes and the time horizon are

    Emax,1 Emax,2 Emax,3 Emax,4     =     1 1 1 1     ,     X0,1 X0,2 X0,3 X0,4     =     0.2 0.2 0.2 0.2     , T = 80. (26)

The optimal fishery strategy is found numerically using both regularisation and time horizon iterations. The parameters used are

N = 100, T0 = 5, ∆T = 0.5, δmin,1= 2−4, δmin,2= 2−9. (27)

(26)

t 0 10 20 30 40 50 60 70 80 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 X1 X2 X3 X4

t 0 10 20 30 40 50 60 70 80 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 X1 E1 X2 E2 X3 E3 X4 E4 (b) Optimal fishery.

The behaviour of the dynamical system in the case of no fishery is very complex, and exhibit aperiodic oscillations. However, Figure 16b infers that the optimal strategy is to perform fishery in such a manner that the population sizes stabilize.

7 Discussion

In this study, fishery was modelled by introducing the concept of fishery effort. The effort put into fishery can be considered as actions yielding larger catch. For example the time spent fishing or the number of active fishery vessel. When considering multiple fish species, there was an effort restricted to each species with an associated maximal effort. However in practice considering a maximal effort for each species might not be an accurate model. Fishing vessels could be equipped with multiple fishing tools, allowing them to vary which fish species to harvest. In this case the accumulated effort put into fishery would have an upper bound. The optimization problem could then be formulated as how to distribute the effort among the different species.

The results indicate that the effort put into fishery of different species may vary. There might be practical issues with this strategy, since controlling the harvest of each species separately may be difficult. If the fishing techniques of each species were different, the species-specific efforts would be possible. However, if this is not the case, more applicable models for the effort must be introduced. For example an effort which is shared among all species could be introduced. The shared effort model would imply that whenever fishery is performed, all fish species are harvested, but the amount of yield of each species could still be dependent of the current population sizes.

Apart from modelling the effort, defining optimal result is very complex. In this study, optimal result was considered as maximal catch over a given time period. This model has many drawbacks. For example, it does not take into account the fish price on the market. In addition, it allows for periods during which no fishery is performed. During these pe-riods, the amount of products from fishery will drop on the market, creating substantial gaps between supply and demand. Therefore, a more applicable definition of optimality must be introduced. For example large variance in the yield over time could be punished.

(27)

Another way would be to give the effort some lower boundary, which would correspond to the minimal effort needed to maintain stable market prices.

In this study the aim has been to find an optimal fishery strategy, assuming that all fishermen follow the same policy. When acting as one unit, the optimal outcome maxi-mizes collective yield. However, in reality this is not the case. Commercial fishery is often open access, and there are many individuals trying to maximize their own yield. Consider the scenario where there are a given number of individuals harvesting a certain region, and that the fish population is such that the accumulated catch is maximal. If another in-dividual A starts harvesting in the same region, the fish population would decrease which reduces the accumulated catch. This might not seem optimal, but from the viewpoint of the A, the catch has increased, since A did not fish in the region previously. This dilemma when individuals act according to their own self-interest negatively affects the common situation is referred to as Tragedy of the commons. Hence even though there is an optimal strategy for the accumulated fishery, it might be unstable.

The dynamics of the fish population has in this study assumed to only be dependent of time. However, in reality, fish do not stay in the same habitat for all times, but migrate to different regions. Hence a more accurate model would consider fish population density, and coordinates describing the position of the fish. The logistic equation would then be expanded into a partial differential equation. A spatially dependent but time independent model has been studied to justify marine reserves [Neu03]. Using the partial differential equation could then possibly result in the position of the marine reserves to be dependent on time.

For a single fish species, the result indicated that the optimal strategy coincide with the maximum sustainable yield (MSY). For most of the examples considered with inter-dependent fish species, the population sizes stabilized to some stationary solution. This includes the example of four interdependent fish species which exhibit chaotic behaviour if no fishery is performed. In Appendix B, it is shown that for the examples of two inter-acting species with stabilized strategy, this stationary numerical solutions coincide with MSY. In the case of competitive exclusion in Section 6.2.2, the optimal solution did not stabilize. However, it is important to keep in mind that only special cases of the dynamical have been studied. Therefore the results do not necessarily hold for the general dynamical systems.

The numerical solution of the competitive exclusion example in Section 6.2.2 converged to two different solutions. This is possible since Pontryagin’s principle only finds local maximum. In order to obtain a global maximum, the HJB equation (Equation (4)) should be solved. However, solving the HJB equation for high dimensional n is computationally

demanding [CMS+_10].

8 Conclusion

In conclusion, Pontryagin’s principle can be used in order to obtain optimal fishery strate-gies. By using a regularization method and iterative algorithms, the optimal policies can be found for various dynamical systems and for large time horizons. The majority of the studied examples had a stationary optimal strategy. However, altering the iteration pro-cess could result in convergence to different solutions. This is possible since Pontryagin’s

(28)

principle only finds local maximums. Future studies could involve solving the HJB equa-tion, which guarantees global maximum. In addition more accurate dynamical systems and applicable definitions of optimality would be of interest to examine.

Acknowledgement

I would like to express my greatest gratitude to my supervisor Anders Szepessy for having devoted his time and interest in the project.

(29)

References

[Apa12] N. Apazidis. Mekanik II. Studentlitteratur, 2nd edition, 2012.

[Cla76] C. W. Clark. Mathematical Bioeconomics: The Optimal Management of

Re-newable Resources. John Wiley & Sons, 1976.

[CMS+_{10] J. Carlsson, K. Moon, A. Szepessy, R. Tempone, and G. Zouraris. Stochastic}

differential equations: Models and numerics. February 2010.

[EHW02] C. Lubich E. Hairer and G. Wanner. Geometric Numerical Integration.

Springer Series in Computational Mathematics. Springer, 2nd edition, 2002.

[GFG35] A. A. Witt G. F. Gause. Behavior of mixed populations and the problem of

natural selection. The American Naturalist, 69(725):596–609, 1935.

[GQ38] J. G. Garnier and A. Qu´etelet. Correspondance math´ematique et physique,

volume 10. L’Observatoire de Bruxelles, 1838.

[MATa] Equation solving algorithms. http://se.mathworks.com/help/optim/ug/

equation-solving-algorithms.html. Accessed: 2016-04-23, 21:20.

[MATb] trapz. http://se.mathworks.com/help/matlab/ref/trapz.html#bua4lsr.

Accessed: 2016-05-22, 13:06.

[Neu03] M. G. Neubert. Marine reserves and optimal harvesting. Ecology Letters,

6(843-849), 2003.

[Sch54] M.B. Schaeffer. Some aspects of the dynamics of populations important to

the management of the commercial marine fisheries. Inter-American Tropical Tuna Comission Bulletin, 1(2):23–56, 1954.

[VWA+_{06] J. A. Vano, J. C. Wildenberg, M. B. Anderson, J. K. Noel, and J. C. Sprott.}

Chaos in low-dimensional lotka–volterra models of competition. Nonlinearity, 19:2391–2404, September 2006.

(30)

Appendix A

In this section the derivation of ˙λ(t) from the Hamiltonian is presented in the general case for n species.

˙

λi(t) = −∂XiH(X(t), λ(t)),

where the Hamiltonian H is H(X, λ) = n X j=1 λj(Fj(X) − E?,j(λj)Xj) + E?,j(λj)Xj 1 − δE?,j(λj) T . Differentiation gives ∂XiH(X, λ) = n X j=1   λj∂XiFj(X) + E?,j(λj) ∂XiXj | {z } δij 1 − δE?,j(λj) T − λj   = = n X j=1 λj∂XiFj(X) + E?,i(λi) 1 T − λi − δ TE?,i(λi) 2_. ₍₂₈₎ Equation (7) gives ∂XiFj(X) = ∂Xi rjXj 1 − 1 Kj n X k=1 wjkXk !! = = rj   ∂XiXj | {z } δij 1 − 1 Kj n X k=1 wjkXk ! −Xj Kj X wjk∂XiXk | {z } δik   = = rj δij 1 − 1 Kj n X k=1 wjkXk ! − Xj Kj wji ! . Therefore n X j=1 λj∂XiFj(X) = n X j=1 rjλj δij 1 − 1 Kj n X k=1 wjkXk ! −Xj Kj wji ! = = riλi 1 − 1 Ki n X k=1 wikXk ! − n X j=1 rjλj Xj Kj wji. (29)

Using Equations (28) and (29) gives ˙ λi(t) = −riλi(t) 1 − 1 Ki n X k=1 wikXk(t) ! + n X j=1 rjwjiλj(t) Xj(t) Kj + +E?,i(λi(t)) λi(t) − 1 T + δ TE?,i(λi(t)) 2_.

(31)

Appendix B

Here the stationary solutions maximizing yield are analysed in the case of two fish species. The analytic solution is lastly compared to the numerical result presented in Section 6.2. The dynamical system is

dXi(t) dt = Fi(X(t)) − Ci(X(t), E(t)), where F1(X) = r1X1 1 −X1+ w12X2 K1 , C1(X, E) = E1X1, and F2(X) = r2X2 1 −w21X1+ X2 K2 , C2(X, E) = E2X2.

The stationary solutions correspond to Fi(X) − Ci(X, E) = 0 for i = 1, 2. For each i,

the curve defined by Fi(X) − Ci(X, E) = 0 projected onto the X1X2-plane is a line. The

equation for each line is

1 −E1 r1 −X1+ w12X2 K1 = 0, and 1 −E2 r2 −w21X1+ X2 K2 = 0. For convenience the following variables are defined,

x1= X1 K1 , x2 = X2 K2 , α1= E1 r1 , α2= E2 r2 , β12= K2 K1 w12, β21= K1 K2 w21. (30)

The equation for the lines are therefore (

x1= 1 − α1− β12x2,

x2= 1 − α2− β21x1

. (31)

The goal is to find the reduced efforts α1 and α2such that the total catch Ctot = C1+C2is

maximized. The optimization problem can also be formulated with x1and x2as variables.

Ctot(x1, x2) = E1X1+ E2X2 = {Equation (30)} = r1K1α1x1+ r2K2α2x2=

= {Equation (31)} = r1K1(x1− x21− β12x1x2) + r2K2(x2− x22− β21x1x2).

Optimization using ∇Ctot(x1, x2) = 0 gives

0 = ∂x1Ctot(x1, x2) = r1K1(1 − 2x1− β12x2) − r2K2β21x2,

0 = ∂x2Ctot(x1, x2) = −r1K1β21x1+ r2K2(1 − 2x2− β21x1).

Thus a system of linear equation can be solved in order to find x1 and x2 that maximizes

stationary yield. 2a c c 2b x1 x2 =a b , where a = r1K1, b = r2K2, c = r1K1β12+ r2K2β21. The solution is x1 = 2ab − bc 4ab − c2, x2= 2ab − ac 4ab − c2. (32)

Note that the optimal fish population sizes are equal if a = b, that is r1K1 = r2K2. The

(32)

The analytical stationary optimal solution can be compared with the numerical result in Section 6.2. Introduce

X₁stat, X₂stat, E₁stat, Estat₂ ,

as the optimal stationary population sizes and efforts. Given the dynamical system, the optimal solution can be calculated using Equations (30)-(32). Figures 17-19 show the comparison with the numerical result.

t 0 1 2 3 4 5 6 7 8 9 10 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 X₁ E₁ X₂ E₂ Xstat 1 Estat 1 Xstat 2 Estat 2

Figure 17: Comparison between numerical result of the competitive coexistence example in Section 6.2.1 and analytical stationary solution. The numerical result coincides with the optimal stationary solution.

Figure 17 infers that the numerical solution approaches the optimal analytical station-ary solution. t 0 1 2 3 4 5 6 7 8 9 10 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 X₁ E₁ X₂ E 2 Xstat 1 Estat 1 Xstat 2 Estat 2

(a) Optimal fishery with δmin,1 = 2−4.

t 0 1 2 3 4 5 6 7 8 9 10 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 X 1 E 1 X 2 E 2 Xstat 1 Estat 1 Xstat 2 Estat 2

Figure 18: Comparison between numerical result of the competitive exclusion example in Section 6.2.2 and analytical stationary solution.

In Figure 18a the population sizes oscillate around the analytical stationary solution. The mean efforts in the oscillatory region coincide with the stationary efforts. Numerical integration shows that the oscillatory solution yields larger reward than the stationary. In Figure 18b there is no clear connection between the numerical result and the analytical stationary solution.

(33)

0 1 2 3 4 5 6 7 8 9 10 t 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 X 1 E 1 X 2 E 2 Xstat 1 Estat 1 Xstat 2 Estat 2

Figure 19: Comparison between numerical result of the competitive exclusion example in Section 6.2.3 and analytical stationary solution. The numerical result coincide with the optimal stationary solution.

Figure 19 shows that the numerical solution stabilizes to the analytical stationary solution.

Optimal Fishery Modelling Population Dynamics with Optimal Harvest Strategies

Optimal Fishery

Modelling Population Dynamics with Optimal Harvest Strategies

Contents

1

Introduction

2

Background

3

Dynamics of Fish Populations

4

Analysis of Optimal Control Problem

5

Numerical Algorithms

6

Results

7

Discussion

8

Conclusion

Acknowledgement

References

Appendix A

Appendix B