# Optimal Harvesting Effort for Nonlinear Predictive Control Model for a Single Species Fishery

Authors: Messaoud Bounkhel; Lotfi Tadj
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

We use nonlinear model predictive
control to find the optimal harvesting effort of a renewable resource system
with a nonlinear state equation that maximizes a nonlinear profit function.
A solution approach is proposed
and discussed and satisfactory numerical illustrations are provided.


---

## Body

## 1. Introduction

Consumption of the world’s natural resources is increasing at a disturbing rate. The United Nations Environment Programme (UNEP) warned that current voracious consumption of resources cannot be sustained.Unlike petroleum, oil, copper, and gold, fish are renewable resources. However, more people are eating fish than ever before, and fish stocks are declining alarmingly. Aquaculture is failing to fill the gaps between the supply and the demand, for lack of a better management, as reported by a recent Food and Agriculture Organization (FAO) review.Better management of fisheries in the high seas, conservation of the biodiversity of ecosystems and species related to it, and reduction of illegal catch of popular and consumed worldwide fish are required to reverse the negative trends threatening fish and the ocean environment on which they depend. Collective actions at all levels and extensive cooperation optimizing the use of depleted resources are needed to help the world abandon the race for fish and adopt an ecosystemic approach that is crucial to ensure the health and future productivity of these key marine ecosystems.Clark [1] summarizes the current worldwide crisis in marine fisheries as “too many boats chasing too few fish.”Since the earliest models of Gordon [2] and Schaefer [3], the fishery resource has received a lot of attention. The basic model of renewable resource exploitation is
(1)dxtdt=Gx-ht,x0=x0,
where the state variable x(t) denotes the biomass of the fish stock at time t, x0 is the initial population of fish, G(x) is the net logistic growth rate of the population biomass, and h(t) is the harvest of the resource stock at time t. The traditional logistic growth rate assumed by most researchers is the logistic model
(2)Gx=rxt1-xtK,
where r>0 is the intrinsic growth rate and K>0 is the carrying capacity of the environment. Introducing the control variable E(t), the harvesting effort at time t, the harvest of the resource, and the harvesting effort are related through thecatch-effort relation
(3)ht=qEtxt,
and q is called thecatchability coefficient. The net economic revenue is
(4)Rt=pqxt-cEt,
where p is the selling price per unit of biomass and c is the cost of harvesting per unit of effort.Ganguly and Chaudhuri [4] use the catch-rate function
(5)ht=qEtxtaEt+bxt,
where a and b are positive constants. Assuming that an external agency regulates the fishery by imposing a suitable tax per unit biomass of landed fish, the fishing effort is taken as a dynamic variable depending on the capital invested in the fishery
(6)Et=αQt,0≤α≤1,dQtdt=It-γQt,
and I(t) is the gross investment rate at time t, Q(t) is the amount of capital invested in the fishery at time t, and γ is the rate of depreciation of capital. Optimal control theory is used to obtain the optimal harvesting policy.Fan and Wang [5] consider the logistic model (2) and assume that r and K are both periodic functions with respect to t:
(7)Gx=rtxt1-xtKt.
They apply qualitative methods and optimal control methods to determine the optimal harvesting policy, including the optimal harvesting time-spectrum.Peng [6] considers the same model as Ganguly and Chaudhuri [4] with the Gompertz law of growth
(8)Gx=rxtln⁡xtK,
instead of the logistic model (2). He also uses optimal control theory to obtain the optimal harvesting policy.Joshi et al. [7] analyze a spatial extension of the dynamic Gordon-Schaefer bioeconomic model. They consider a stock density u(x,t) that diffuses and is advected within a smooth, bounded habitat Rn for a finite length of time T. The diffusion coefficients are aij for i,j=1,…,n, the advection coefficients are bi for i=1,…,n, the stock grows at a rate f(u) that depends only on the local stock density, and the harvest rate h(x,t) is proportional to both the stock density and the effort density:
(9)ut=fu-hu+∑i,j=1naijx,tuxixj+∑i=1nbix,tuxi.
The subscripts represent partial derivatives. Joshi et al. prove the existence of an optimal control and derive the necessary conditions that an optimal control must satisfy.E. Braverman and L. Braverman [8] study the optimal harvesting strategy for populations whose dynamics is described by reaction-diffusion equations. They consider three production functions: logistic, Gilpin-Ayala, and Gompertz type. They investigate the maximum yield for both continuous and impulsive models.In Halkos and Papageorgiou [9], the variations of the renewable resource stock evolve according to the differential equation
(10)dxtdt=fxt-qt,
while the discounted net economic revenue is given by
(11)R=e-ρtpt-cxtqt,
where ρ>0. They use optimal control theory to derive the conditions under which the renewable resource harvesting model achieves a unique steady-state equilibrium.Duncan et al. [10] use a nonlinear dynamic discount rate δ(t) instead of the linear discount rate ρt, to reflect the fact that the utility derived from a harvest may be worth more at one point in time than another. The discount factor e-ρt in (11) will be replaced with the more general one given by
(12)Dt=exp⁡-∫0tδτdτ.Optimal control theory has been extensively used to determine the optimal harvesting policy for renewable resources such as fish stocks. Not only in the extensions of the basic model that we have described above, but also in more integrated models which involve two or more species, structured models, a population of consumers, predator-prey models, reserve-unreserve areas, and so forth.Our intention in this paper is to use a different approach, model predictive control (MPC). Model predictive control for linear constrained systems provides excellent control solutions both theoretically and practically. Many systems, however, such as in renewable resources, are inherently nonlinear. This motivates the use of nonlinear model predictive control. Basically, in model predictive control an optimal control problem is solved for the current system state. MPC is based on an iterative process over finite horizon. At timet0 the current state is sampled and a cost minimizing control strategy is computed for a relatively short prediction time T. Specifically, an online computation is used to explore state trajectories from the current state until the end of the prediction interval [t0,t0+T]. Only the first step of the optimal control is implemented; then the state is sampled again at time t0+T. The computations are repeated starting from the current state, yielding a new control and new predicted state path. We mention here that the online nonlinear optimization step is always an issue, and thus suboptimal MPC algorithms with online linearization and quadratic optimization have been used in process control. For more details, see, for example, Ellis et al. [11] and Simon [12] and the references therein.MPC is an advanced method of process control that has been successfully used in the process industries, especially in chemical processes; see, for instance, Goodwin et al. [13], Qin and Badgwell [14], and the references therein. Among some recent references, we cite Petersen and Jorgensen [15] who use MPC to maximize profit of the fermentation process, which is a widely used process in production of many foods, beverages, and pharmaceuticals. del Favero et al. [16] report the first wearable artificial pancreas outpatient study based on MPC and investigate specifically its ability to control postprandial glucose, one of the major challenges in glucose control. Sun et al. [17] study trajectory generation for a mothership that tows a drogue using a flexible cable. The optimal trajectory for the towed cable system with tension constraints is generated using MPC. Karamanakos et al. [18] present an MPC approach for dc-dc boost converters. Berkemeier et al. [19] use nonlinear MPC for vehicle control and explore whether straight-forward application results in computations take too long for real-time use. Methods for speeding up the computations are discussed. Pytel and Kozak [20] deal with the effective predictive control algorithm for the gas turbine on the base of a mathematical model obtained from measured I/O data. Finally, Bréchet et al. [21] use MPC to study the change of the atmospheric temperature within the next 150 years.Clark [22] argues that when the fisherman operates his vessel so as to obtain the largest possible income from each day’s fishing, the daily effort cost c(E) is assumed to be nonlinear with increasing marginal cost c′(E). In agreement with this argument, we define an objective in which the cost, and therefore the profit to maximize, is nonlinear. Since the dynamics of the problem are already nonlinear, the result is a highly nonlinear formulation. To avoid the computational burden in the online optimization phase of nonlinear model predictive control (NMPC), approximate solutions are sought. The approximations employed are good enough and largely compensate for the extra effort required to reach optimal solutions.The model is formulated in Section2 and solved in Section 3. Illustrative examples are presented in Section 4. Section 5 presents a conclusion and future research directions.


## 2. Model Formulation

LetH>0 and consider on the planning interval [0,H] a renewable resource, such as a fishery, whose population dynamics are governed by a state equation of the type
(13)ddtBt=FBt-ht,B0=B0,
where the state variable B(t) represents the population biomass at time t, F(B) is the logistic growth rate function, and h(t) is the harvesting rate at time t.We use the logistic growth rate function(14)FB=rBt1-BtK,
where r>0 is the intrinsic growth rate and the nonnegative integer K is the carrying capacity of the fish population.We also use the rate of harvesth(t) based on the constant “catch-per-unit-effort” and usually considered in fishery models
(15)ht=qEtBt,
where q>0 is the catchability coefficient and the control variable E(t) is the effort of harvesting at time t. Thus, the differential equation (13) becomes
(16)ddtBt=rBt1-BtK-qEtBt.
Finally, let t0∈[0,H] and consider the prediction interval [t0,t0+T], where T≪H. To write the objective function for our nonlinear model predictive control model, we denote by p>0 the constant price per unit biomass and by c>0 the constant cost of harvesting per unit biomass. Assuming a nonlinear (quadratic) cost (see Clark [22]), we seek to maximize the profit
(17)Jt0,E=k2Bt0+T+∫t0t0+TpqBτEτ-cEτ2dτ,
where k∈(0,c] represents the salvage value of the ending state.An NMPC approach is used in the next section to determine the control variable at timet that maximizes the profit function (17) subject to the state equation (16).


## 3. Model Solution

Different techniques have been proposed in the literature to speed up the calculation of the optimal control variable of the problem stated above. We use an approximate calculation of the integral in the objective function (17) as follows. Put
(18)Fτ=pqBτEτ-cEτ2.
Now divide the time interval [t0,t0+T] into m subintervals of equal length h=T/m; then use the trapezoid formula for m intervals. The objective function (17) becomes
(19)Jt0,E≃k2Bt0+mh+h2Ft0+2∑i=1m-1Ft0+ih+Ft0+mh,
where, as in (18), we have
(20)Ft0+ih=pqBt0+ihEt0+ih-cEt0+ih2.
In order to calculate the sum that appears in (19), we write the linear approximation of the state variable B(t), which in conjunction with the state equation (16) yields
(21)Bt0+ih≃Bt0+ihdBtdt≃C1t0,i-C2t0,iEt0,
where
(22)C1t0,i=Bt01+ihr1-Bt0K,C2t0,i=ihqBt0.
Now, substitute (21) in (20) to get
(23)Ft0+ih≃pqC1t0,i-C2t0,iEt0Et0+ih-cEt0+ih2.
We thus obtain the second term in the right-hand side of objective function (19) as
(24)h2Ft0+2∑i=1m-1Ft0+ih+Ft0+mh≃hpq2Bt0Et0-ch2Et02+hpq∑i=1m-1C1t0,iEt0+ihhhhhhhh-hpqEt0∑i=1m-1C2t0,iEt0+ihhhhhhhh-ch∑i=1m-1Et0+ih2+hpq2C1t0,mEt0+mhhhhhhhh-hpq2C2t0,mEt0Et0+mhhhhhhhh-ch2Et0+mh2.
Also, the first term in the right-hand side of objective function (19) is given by
(25)k2Bt0+mh=k2C1t0,m-C2t0,mEt0.
Combining (24) and (25), and after some simple computations, the objective function (19) becomes
(26)Jt0,E≃Mt0-Gt0⊤Et0+Et0⊤Qt0Et0,
where M(t0)≔(k/2)C1(t0,m) is independent of the control, E(t0) and G(t0) are the (m+1)-tuple vectors defined by
(27)Et0≔Et0,Et0+h,hhhEt0+2h,…,Et0+mh⊤;Gt0≔-hpq12Bt0-k2hpqC2t0,m,C1t0,1,hhhhhhhC1t0,2,…,C1t0,m-1,12C1t0,m⊤,
and Q(t0) is a square matrix of order (m+1):(28)Qt0≔-ch12pqcC2t0,1pqcC2t0,2⋯pqcC2t0,m-1pq2cC2t0,m010⋯00001⋯00⋮⋮⋮⋮⋮⋮000⋯10000⋯012.Clearly, its inverse Q(t0)-1 exists (since c>0 and h>0) and(29)Qt0-1≔-1ch2-2pqcC2t0,1-2pqcC2t0,2⋯-2pqcC2t0,m-1-2pqcC2t0,m010⋯00001⋯00⋮⋮⋮⋮⋮⋮000⋯10000⋯02.Note that the matrix Q(t0) is negative definite. The global mmaximum of J(t0,·) is reached at E(t0) such that
(30)Et0=12Qt0-1Gt0.
In receding horizon, we obtain E(t0) using the formula
(31)Et0=1,0,0,…,0·Et0
and hence
(32)Et0=pqc212Bt0-k2hpqC2t0,m∑i=1m-1C1t0,iC2t0,ihhhhhh-2pqc∑i=1m-1C1t0,iC2t0,ihhhhhh-∑i=1m-1C1t0,iC2t0,ipqcC1t0,mC2t0,m.First Term in (32).Consider
(33)212Bt0-k2hpqC2t0,m=Bt0-khpqC2t0,m=Bt0-khpqmqhBt0=1-kmpBt0.Second Term in (32). Let α=∑i=1m-1‍i=m(m-1)/2 and β=∑i=1m-1‍i2=m(m-1)(2m-1)/6:
(34)2pqc∑i=1m-1C1t0,iC2t0,i=2pqc∑i=1m-1Bt01+ihr1-Bt0KihqBt0=2pqcBt02hqα+βhr1-Bt0K=2hpq2cBt02α+βhr-βhrKBt0=2hpq2cα+βhrBt02-βhrKBt03.Third Term in (32).Consider
(35)pqcC1t0,mC2t0,m=pqcmBt0hqBt01+mhr1-Bt0K=hmpq2cBt021+hmr-hmrKBt0=hmpq2c1+hmrBt02-hmrKBt03.
Therefore (32) becomes
(36)Et0=pqc1-kmpBt0hmpq2chhhhhhhh-2hpq2cα+βhrBt02-βhrKBt03hhhhhhhh-hmpq2c1+hmrBt02-hmrKBt03.
Since the choice of t0 was arbitrary in [0,H] we deduce the general relationship between the optimal control E and the state B over the whole horizon interval [0,H] as follows:
(37)E·=B·b1+b2B·+b3B·2,
where
(38)b1≔pqc1-kmp,b2≔-pqc2hpq2α+βhrc+hmpq21+hmrc=-hp2q3c22α+βhr+m1+hmr,b3≔pqc2hpq2cβhrK+hmpq2chmrK=h2p2q3rc2K2β+m2.Substituting (37) in the state equation (16) we obtain the following differential equation:
(39)ddtBt=rBt1-BtK-qBt2b1+b2Bt+b3Bt2,
which will be solved numerically over the interval [t0,t0+T] with the initial condition B(t0)=Bt0 to get the solution B over [t0,t0+T] for any t0∈[0,H]. Once B is found we replace its value in E to find the optimal value of the harvesting effort.Since the harvesting effortE must be nonnegative due to the real life assumptions on the model, we use the following maximum formula, usually used in optimal control theory (see, for instance, Sethi and Thompson [23]), for the optimal control E*:
(40)E*·=max⁡0,E·.


## 4. Illustrative Examples

We provide in this section simulation examples to show different types of solutions that can be obtained using the results obtained in the previous section. For a given time horizonH=10, we take for all the following simulations an instant time t0=2 and a prediction horizon T=1.25.Case 1.
Consider the following values of parameters:c=0.1, r=0.6, K=100, q=0.01, k=0.1, p=5, h=0.05, m=25, and B0=0.2. Figure 1 shows the variations of the optimal state and control variables. To obtain these graphs, first, the differential equation (39) is solved numerically using the mathematical package MATLAB. The result is given in the form of the left graph in Figure 1. The biomass of the fish stock increases monotonically, starting from the initial state Bt0=0.2. Then, using (37), the optimal harvesting effort is obtained as the right graph in Figure 1. It is also a monotonically increasing function of time. Finally, the optimal objective profit at any time t0 is evaluated by substituting the expression of the optimal control (30) in (26) and we obtain
(41)J*t0=Mt0-14Gt0⊤Qt0-1Gt0.Figure 1
Variations ofB* and E* as functions of t on [t0,t0+T] in Case 1.Case 2.
We keep all the parameters constant as in Case1, except for the growth coefficient which is taken to be r=1.9 instead of 0.6. The results are depicted in Figure 2. We note that the optimal biomass of the fish stock is again increasing as in Case 1; however it reaches a much higher level since the growth rate is much higher. The optimal harvesting effort increases at first to reach a maximum effort and then decreases to reach the value 0 at time t=2.983 and remains 0 on the interval [2.983,3.25]. The reason the harvesting effort is 0 is that (37) is negative on this interval.Figure 2
Variations ofB* and E* as functions of t on [t0,t0+T] in Case 2.Case 3.
Again, we keep all the parameters constants as in Case1, except for the initial biomass which is taken to be Bt0=2 instead of 0.2. The results are plotted in Figure 3. The biomass increases monotonically starting from the initial value 2. It is much higher than the biomass of Case 1. The resulting harvesting effort is monotonically decreasing during the prediction interval [t0,t0+T]=[2,3.25].Figure 3
Variations ofB* and E* as functions of t on [t0,t0+T] in Case 3.Case 4.
Once again, we keep all the parameters constants as in Case1, except for the growth coefficient which is taken to be r=0.0004 instead of 0.6 and for the constant price per unit biomass p=20 instead of 5. The results are plotted in Figure 4. Contrarily to the previous cases, the biomass in this case is decreasing monotonically during the prediction interval.Figure 4
Variations ofB* and E* as functions of t on [t0,t0+T] in Case 4.


## 5. Conclusion

Optimal control theory has been plentifully used to determine at timet=0 the optimal harvesting effort on a planning interval [0,H]. In contrast, starting at time t0=0 with the current state B(0)=B0, NMPC determines the optimal harvesting effort on the prediction interval [t0,t0+T]=[0,T]. Then at time t0=T with the current state B(t0)=B(T), the optimal harvesting effort is determined on the prediction interval [t0,t0+T]=[T,2T]. This process is repeated over and over until time H. The determination of the optimal control on the time intervals [t0,t0+T] is the step that requires the most calculations. We overcome difficulty by using a judicious approximation of the objective function calculation. The results obtained are easily implementable as shown in Section 4.The method described in this paper is quite robust and we propose to further experiment it on more complex models. For example, the catch-rate function (5) could be used instead of the function (3); the intrinsic growth of rate r and the carrying capacity of the fish population could be periodic functions; the Gompertz law of growth (8) could be used instead of the logistic growth rate (2); and/or a discounted cost function with either a constant or dynamic discount factor could be used.



---

*Source: 367593-2015-02-12.xml*