# Model Predictive Control of the Harvesting Effort of a Sustainable Seafood with a Nonlinear State Equation

Authors: Messaoud Bounkhel; Lotfi Tadj; Ramdane Hedjar
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

This paper is concerned with the nonlinear model predictive control of the harvesting effort of a renewable resource system, namely, sustainable seafood (farmed fish), with a nonlinear state equation. A solution approach is proposed and discussed, and numerical illustrations are provided.


---

## Body

## 1. Introduction

A lot of the resources we use daily cannot be recovered once they are gone. This type of resources is called nonrenewable resources. For the reason of the overconsumption of nonrenewable resources due to the exponential evolution of the industry, we have to turn to another type of resources which can replenish with the passage of time, either through biological reproduction or other naturally recurring processes. This type of resources is called renewable resources, such as biomass resources (fish, trees, food crops, algae, agricultural and forestry byproducts, and even methane fumes from landfills), water, wind, and solar. Renewable resources may be the source of many types of power for renewable energy. However, if the exploitation rate or consumption rate of the renewable resource exceeds its renewal rate, the renewal will not be ensured. Then, we could run out of renewable resources, and for this reason, we should turn to sustainable fishing. Solutions to this problem involve several plans of action not only by governments but also by consumers and communities; everyone must play a larger role in trying to achieve better management of renewable resources. Mathematicians can also participate in this effort by mathematically modeling the phenomenon as an optimal control problem and solving it for the optimal harvesting rate which will ensure sustainable fishing.The traditional optimal control theory has been used extensively to effectively and efficiently manage the fishery resource. A large number of papers are available. For an account on the theory and applications, the reader is referred to the books by Clark [1, 2].Model predictive control (MPC) is usually formulated as the repeated solution of an optimal control problem. Based on measurements obtained at timet0, the controller determines the input such that a predetermined objective function is minimized over a control horizon t0,t0+T. The optimal input is implemented only until the next sampling instant t0+T. Using the new system state at the new sampling time, the procedure is repeated. Allgöwer et al. [3] summarizes a standard NMPC scheme work as follows:(1)
Obtain estimates of the states of the system(2)
Calculate an optimal input minimizing the desired cost function over the prediction horizon using the system model for prediction(3)
Implement the first part of the optimal input until the next sampling instant(4)
Continue with (2)It is customary to distinguish between linear MPC (LMPC) used when the system dynamics are linear models and its extension, nonlinear MPC (NMPC) used when the system dynamics are nonlinear models.MPC has received a great interest during the last decades in the process industries, especially in chemical processes. Goodwin et al. [4] reported more than 2000 applications of MPC. There are many excellent reviews of MPC. We cite the most recent one, Qin and Badgwell [5], which cites the others, and also reports more than 4,500 applications.In the present paper, we aim at using the nonlinear model predictive control, an optimal control approach that we have not seen applied in the context of nonrenewable resources. We consider a resource biomass subject to harvesting. The evolution of the resource biomass is described by a nonlinear differential equation. The optimal harvesting effort is obtained in the closed form. The model considered is nonlinear, and hence, NMPC is used to study a tracking problem. We hope to contribute to increase the knowledge about NMPC to stimulate new research in this promising field.We formulate the model in Section2 and develop the solution procedure in Section 3. Simulation results are presented in Section 4. Section 5 discusses the importance of our results, and Section 6 concludes the paper.


## 2. Model Formulation

We consider a resource biomass of densityBt and denote by Et the harvesting effort at time t. Denoting by q>0 the catchability coefficient, by r>0 the growth rate of the resource, and by K>0 its carrying capacity, the dynamics of the resource during the planning horizon 0,H are governed by the usual ordinary differential equation (see, for instance, Dubey et al. [6] and the references therein): (1)ddtBt=rBt1−BtK−qEtBt,B0=B0.Define the targetsB^ (goal density of resources biomass) and E^ (goal harvesting effort) and introduce the shifted variables:(2)ΔBt=Bt−B^,ΔEt=Et−E^.It is desired to have the state variableBt converge toward its goal B^ and the control variable Et toward its goal E^. Note that, in order to have ΔE⟶0, one has to choose the goal harvesting effort:(3)E^=rq1−B^K.Penaltiesℓ and p are incurred for each variable to deviate from its goal so that at any instant of time t0∈0,H and given a prediction horizon of length T>0 (with T≪H), the controller seeks to determine the harvesting rate over t0,t0+T that minimizes the following objective function:(4)Jt0,E=k2ΔBt0+T2+∫t0t0+Tℓ2ΔBτ2+p2ΔEτ2dτ,where k>0 represents the penalty of the terminal state.


## 3. Model Analysis

It is well known that the main drawback of MPC is that the computational burden involved in the optimization stage is fairly expensive. This has limited its practical use for nonlinear systems. Various techniques have been used to alleviate this burden. In practice, successive linearizations which result in simple quadratic optimization are used, for example, Lawrynczuk [7].It is widely accepted that approximate solutions are preferred to exact ones to remove the computational encumbrance from the online implementation. Our approach consists in using the trapezoidal rule to approximate the definite integral in (4). To this end, put(5)Fτ=ℓ2ΔBτ2+p2ΔEτ2,and divide the prediction horizon t0,t0+T into m subintervals of equal length h=T/m. Using the trapezoid formula for m intervals, the objective function (4) becomes(6)Jt0,E≃k2ΔBt0+mh2+h2Ft0+2∑i=1m−1Ft0+ih+Ft0+mh,where, as in (5),(7)Ft0+ih=ℓ2ΔBt0+ih2+p2ΔEt0+ih2.To deal with the termΔBt0+ih2 in (7), we write the linear approximation of the state variable ΔB as follows:(8)ΔBt0+ih≃ΔBt0+ihdBtdt≃ΔBt0+ihBt0r1−Bt0K−qE^−ihqΔEt0Bt0.For simplicity, let(9)A1t0:=hBt0r1−Bt0K−qE^,A2t0:=hqBt0.Then,(10)ΔBt0+ih≃ΔBt0+iA1t0−iA2t0ΔEt0.Taking the square, we have(11)ΔBt0+ih2≃ΔBt0+iA1t0−iA2t0ΔEt02≃C1t0,i−C2t0,iA2t0ΔEt0+i2A2t02ΔEt02,where(12)C1t0,i:=ΔBt02+2iA1t0ΔBt0+i2A1t02,C2t0,i:=2iΔBt0+2i2A1t0.Let(13)C1t0=ℓ∑i=1m−1C1t0,i,C2t0=ℓA2t0∑i=1m−1C2t0,i,C3t0=ℓA2t02∑i=1m−1i2=δA2t02,δ:=ℓmm−12m−16.After some lengthy computations, we put the objective function (6) under its final form:(14)Jt0,E=M1t0−Gt0TΔEt0+ΔEt0Tℚt0ΔEt0,where (15)M1t0:=hℓ4ΔBt02+h2C1t0+hℓ+2k4C1t0,m,ΔEt0:=ΔEt0,ΔEt0+h,ΔEt0+2h,…,ΔEt0+mhm+1×1T,Gt0:=h2C2t0+hℓ+2k4C2t0,mA2t0,0,…,0m+1×1T,ℚt0:=diagph4+h2C3t0+hℓ+2km24A2t02,ph2,…,ph2,ph4m+1×m+1.Note thatM1t0 is independent of the control and that the diagonal matrix ℚt0 is positive definite and invertible provided p≠0. We derive the optimal control in the following two cases:Case 1 (p > 0).
Clearly, the minimum ofJt0 is reached at ΔE∗t0 such that(16)ΔE∗t0=12ℚt0−1Gt0.
As mentioned in Introduction, in the receding horizon, only the first component of the vectorΔE∗t0 is implemented until the next sampling instant and then the whole process is repeated. From (16), we obtain ΔEt0 explicitly as(17)ΔEt0=1,0,0,⋯,0⋅ΔE∗t0=q−δ1Bt0+δ2Bt02−δ3Bt03p+δ4Bt02,where(18)δ1:=mhℓm+2kB^,δ2:=mhℓm+2k+hrK2δh+m2hℓ+2kB^,δ3:=hrK2δh+m2hℓ+2k,δ4:=2δ+hℓ+2km2qhq.
Since our previous analysis is valid for anyt0∈0,H, we substitute the expression E⋅=E^+q−δ1B⋅+δ2B⋅2−δ3B⋅3/p+δ4B⋅2 in the state equation (1) to obtain the following differential equation:(19)ddtBt=r−qE^Bt−rKBt2+q2δ1Bt2−δ2Bt3+δ3Bt4p+δ4Bt2,which can be solved numerically on the interval t0,t0+T with the initial value Bt0=Bt0, for any t0∈0,H.Case 2 (p = 0).
In this case, the matrixℚt0−1 does not exist; however, the objective function has the following simpler form:(20)Jt0,E=M1t0−Gt0ΔEt0+Qt0ΔEt02,where(21)M1t0:=rh4ΔBt02+h2C1t0+rh+2k4C1t0,m,Gt0:=h2C2t0+rh+2k4C2t0,mA2t0,Qt0:=h2C3t0+rh+2km24A2t02>0.
In this case, the minimum ofJt0 is reached at(22)ΔEt0=Gt02Qt0,which yields after computations:(23)Et0=E^+qδ2δ4−qδ1δ4Bt0−qδ3δ4Bt0.
Similarly to the previous case, the state equation is used to find the state variableBt over the interval t0,t0+T as follows:(24)ddtBt=r−qE^Bt−rKBt2−qBtE^+qδ2δ4−qδ1δ4Bt−qδ3δ4Bt=q2δ1δ4+r−2qE^−q2δ2δ4Bt+q2δ3δ4−rKBt2=a+bBt+cBt2,where(25)a=q2δ1δ4>0,b=r−2qE^−q2δ2δ4,c=q2δ3δ4−rK.
Observe thatc=2δrhq−1/K2δ+m2qrh+2k<0 and so c≠0 which means that the previous differential equation is a Riccati differential equation and has an explicit solution. To find that solution, we compute Δ=b2−4ac, and we find out that Δ>0, and hence, the solution has the form(26)Bt=w1+e2wt+k0c1−e2wt+k0−b2c,where w=Δ/2 and k0=1/2lnB0+b/2c−w/c/B0+b/2c+w/c.


## 4. Numerical Simulations

We implement the results obtained in the previous section for the casep>0, in which the differential equations are solved numerically. The case p=0 is much easier due to the explicit form of the solutions of the differential equations.


### 4.1. Base Case

For a given time horizonH=240, we take for all the following simulations an instant time t0=0 and a prediction horizon of length T=60.Consider the following values of the parameters:p=1, r=0.6, ℓ=1, k=1; K=100, q=0.1, B^=10, E^=r/q1−B^/K=1.35, m=1500, h=T/m=0.04; and B0=B0=40. Figures 1 and 2 show that the solutions Bt (the blue continuous curve) and Et (the blue discontinuous curves) converge to their corresponding targets B^ (the multicolored line) and E^ (the golden line) over the interval t0,t0+T, respectively. Figure 3 shows that the objective function goes to zero, as we expected, at the end of the prediction interval. We point out that the value of J is computed at any time t0∈0,H by substituting the expression of the optimal control (16) in equation (14), and we obtain(27)J∗t0=M1t0−14Gt0Tℚt0−1Gt0.Figure 1
Variations ofBt and B^ as functions of time (t).Figure 2
Variations ofEt and E^ as functions of time (t).Figure 3
Variations ofJt as a function of time (t).In this case, we foundJ∗t0=J∗0=239.13×106.


### 4.2. Sensitivity Analysis

For the study of sensitivity of the model, we varied three parametersr,k,andℓ, and for the other parameters, the same study can be conducted. Consider the following values of the parameters:(28)r∈0.1,0.2,0.3,0.4,H=60,t0=0,T=6,ℓ=1,k=1,K=100,q=0.1,p=1,B^=10,m=1500,B0=B0=40,k∈1,4,12,20,H=60,t0=0,T=6,r=0.4,p=1,ℓ=1,K=100,q=0.1,B^=10,E^=rq1−B^K=1.35,m=1500,B0=B0=40,ℓ∈1,10,20,30,40,50,60H=60,t0=0,T=6,k=1,r=0.4,p=1,K=100,q=0.1,B^=10,E^=rq1−B^K=1.35,m=1500,B0=B0=40.Our numerical results are presented in Figures4–6. Figure 4 shows that the objective function increases when the parameter r increases from 0.1 to 0.4. Similarly, Figures 5 and 6 show that J increases when the values of the parameters k and ℓ increase, respectively, from 1 to 20 and 1 to 60. These results are exactly what we expected.Figure 4
Variations ofJt for r=0.1,0.2,0.3,0.4.Figure 5
Variations ofJt for k=1,4,12,20.Figure 6
Variations ofJt for ℓ=1,10,20,30,40,50,60.


## 4.1. Base Case

For a given time horizonH=240, we take for all the following simulations an instant time t0=0 and a prediction horizon of length T=60.Consider the following values of the parameters:p=1, r=0.6, ℓ=1, k=1; K=100, q=0.1, B^=10, E^=r/q1−B^/K=1.35, m=1500, h=T/m=0.04; and B0=B0=40. Figures 1 and 2 show that the solutions Bt (the blue continuous curve) and Et (the blue discontinuous curves) converge to their corresponding targets B^ (the multicolored line) and E^ (the golden line) over the interval t0,t0+T, respectively. Figure 3 shows that the objective function goes to zero, as we expected, at the end of the prediction interval. We point out that the value of J is computed at any time t0∈0,H by substituting the expression of the optimal control (16) in equation (14), and we obtain(27)J∗t0=M1t0−14Gt0Tℚt0−1Gt0.Figure 1
Variations ofBt and B^ as functions of time (t).Figure 2
Variations ofEt and E^ as functions of time (t).Figure 3
Variations ofJt as a function of time (t).In this case, we foundJ∗t0=J∗0=239.13×106.


## 4.2. Sensitivity Analysis

For the study of sensitivity of the model, we varied three parametersr,k,andℓ, and for the other parameters, the same study can be conducted. Consider the following values of the parameters:(28)r∈0.1,0.2,0.3,0.4,H=60,t0=0,T=6,ℓ=1,k=1,K=100,q=0.1,p=1,B^=10,m=1500,B0=B0=40,k∈1,4,12,20,H=60,t0=0,T=6,r=0.4,p=1,ℓ=1,K=100,q=0.1,B^=10,E^=rq1−B^K=1.35,m=1500,B0=B0=40,ℓ∈1,10,20,30,40,50,60H=60,t0=0,T=6,k=1,r=0.4,p=1,K=100,q=0.1,B^=10,E^=rq1−B^K=1.35,m=1500,B0=B0=40.Our numerical results are presented in Figures4–6. Figure 4 shows that the objective function increases when the parameter r increases from 0.1 to 0.4. Similarly, Figures 5 and 6 show that J increases when the values of the parameters k and ℓ increase, respectively, from 1 to 20 and 1 to 60. These results are exactly what we expected.Figure 4
Variations ofJt for r=0.1,0.2,0.3,0.4.Figure 5
Variations ofJt for k=1,4,12,20.Figure 6
Variations ofJt for ℓ=1,10,20,30,40,50,60.


## 5. Discussion

The growing human needs for more food and energy have led to increased exploitation of several biological resources. Humans today extract and use around 50% more natural resources than only 30 years ago, at about 70 billion tons of raw materials a year. On the contrary, there is a global concern to protect the ecosystem at large. A sustainable resource use is required: high resource consumption is not a requirement for high quality of life.To avoid overexploitation, and hence possibly extinction, of a resource, policies have to be devised to restrict the amount of fish caught. Many alternatives have been proposed, e.g., Behringer and Upmann [8]:(i)
The fisher undertakes several fishing journeys within the fixed time horizon(ii)
The fisher is granted a residual or salvage value for the stock of the resource remaining at the end of the fishing period(iii)
The establishment of harvesting quotas or catch shares(iv)
The specification of a minimum mesh sizeWe propose in this paper the alternative of specifying a goal density of resources biomassB^. This could be, for example, a fraction of the carrying capacity K. This will automatically generate a goal harvesting effort E^ by (3). A regulating agency imposes penalties when either the density of resources biomass or harvesting effort deviates from its goal. This approach of specifying a target is widely used in other areas such as engineering (see, for example, Udwadia [9] and the references therein, Çimen and Banks [10], Tan et al. [11], and Hedjar and Bounkhel [12]) and production planning (see Hedjar et al. [13], Hedjar et al. [14], Tadj et al. [15], and Bounkhel et al. [16]).A novel and effective use of the MPC method is successfully demonstrated in the numerical simulation section. Contrary to the maximum principle which determines at time zero the optimal harvesting rate for the whole planning interval0,H, MPC achieves a better and more flexible control by determining the optimal harvesting rate at any time t0∈0,H on successive small prediction intervals t0,t0+T, until the end of the planning horizon. On each small interval, the state and control variables are tracked toward their respective goal.


## 6. Summary and Future Research

We have considered in this paper a resource biomass subject to harvesting. It is sought to maintain the density of the resource biomass and the harvesting effort at some desired levels. Since the dynamics of the system are nonlinear, a nonlinear model predictive approach is applied. The optimal harvesting effort is obtained in closed form. When the penalty for the harvesting effort to deviate from the desired harvesting rate is zero, the optimal density of the resource biomass is obtained explicitly too. When this penalty is not zero, the optimal density of the resource biomass is obtained as the numerical solution of a differential equation. In this case, we have provided some numerical example to show how the optimal state and control variables are obtained. We have also calculated the optimal objective function value and numerically shown the effect of the penalties on this optimal cost.This first application of nonlinear model predictive control to renewable resources can be further investigated in different ways. We propose the following possible extensions that may be worth investigating:(i)
The catchability coefficientq is defined as the average proportion of a stock that is taken by each unit of fishing effort. We have assumed here that it is constant. However, q may be dynamic (increase or decrease over time) due to many factors such as changes in fishing methods or technology or gain of experience and skill over time. Alternatively, q may be state dependent instead of time dependent.(ii)
When writing the dynamics of the model, we have assumed in (1) that the harvesting rate function at any time t is given by ht=qEtBt. Other harvesting rate functions such as ht=qEtBt/aEt+bBt and ht=qEtBt/1+aEt+bBt, where a and b are positive constants, have been taken by other researchers.(iii)
We have mentioned in this paper that the goal density of the resource biomassB^ should be a fraction of the environment capacity K. An optimal estimation of B^ could be an interesting research question.



---

*Source: 5823465-2019-08-08.xml*