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

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.

Messaoud Bounkhel 1 and Lotfi Tadj 2

Academic Editor:Shihua Li

1, Department of Mathematics, College of Science, King Saud University, Riyadh 11451, Saudi Arabia
2, Department of Information Systems and Decision Sciences, Fairleigh Dickinson University, Vancouver, BC, V6B 2P6, Canada

Received 8 October 2014; Revised 7 January 2015; Accepted 27 January 2015; 12 February 2015

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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 [figure omitted; refer to PDF] where the state variable [figure omitted; refer to PDF] denotes the biomass of the fish stock at time [figure omitted; refer to PDF] , [figure omitted; refer to PDF] is the initial population of fish, [figure omitted; refer to PDF] is the net logistic growth rate of the population biomass, and [figure omitted; refer to PDF] is the harvest of the resource stock at time [figure omitted; refer to PDF] . The traditional logistic growth rate assumed by most researchers is the logistic model [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is the intrinsic growth rate and [figure omitted; refer to PDF] is the carrying capacity of the environment. Introducing the control variable [figure omitted; refer to PDF] , the harvesting effort at time [figure omitted; refer to PDF] , the harvest of the resource, and the harvesting effort are related through the catch-effort relation [figure omitted; refer to PDF] and [figure omitted; refer to PDF] is called the catchability coefficient . The net economic revenue is [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is the selling price per unit of biomass and [figure omitted; refer to PDF] is the cost of harvesting per unit of effort.

Ganguly and Chaudhuri [4] use the catch-rate function [figure omitted; refer to PDF] where [figure omitted; refer to PDF] and [figure omitted; refer to PDF] 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 [figure omitted; refer to PDF] and [figure omitted; refer to PDF] is the gross investment rate at time [figure omitted; refer to PDF] , [figure omitted; refer to PDF] is the amount of capital invested in the fishery at time [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] 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 [figure omitted; refer to PDF] and [figure omitted; refer to PDF] are both periodic functions with respect to [figure omitted; refer to PDF] : [figure omitted; refer to PDF] 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 [figure omitted; refer to PDF] 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 [figure omitted; refer to PDF] that diffuses and is advected within a smooth, bounded habitat [figure omitted; refer to PDF] for a finite length of time [figure omitted; refer to PDF] . The diffusion coefficients are [figure omitted; refer to PDF] for [figure omitted; refer to PDF] , the advection coefficients are [figure omitted; refer to PDF] for [figure omitted; refer to PDF] , the stock grows at a rate [figure omitted; refer to PDF] that depends only on the local stock density, and the harvest rate [figure omitted; refer to PDF] is proportional to both the stock density and the effort density: [figure omitted; refer to PDF] 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 [figure omitted; refer to PDF] while the discounted net economic revenue is given by [figure omitted; refer to PDF] where [figure omitted; refer to PDF] . 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 [figure omitted; refer to PDF] instead of the linear discount rate [figure omitted; refer to PDF] , 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 [figure omitted; refer to PDF] in (11) will be replaced with the more general one given by [figure omitted; refer to PDF]

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 time [figure omitted; refer to PDF] the current state is sampled and a cost minimizing control strategy is computed for a relatively short prediction time [figure omitted; refer to PDF] . Specifically, an online computation is used to explore state trajectories from the current state until the end of the prediction interval [figure omitted; refer to PDF] . Only the first step of the optimal control is implemented; then the state is sampled again at time [figure omitted; refer to PDF] . 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, Brechet 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 [figure omitted; refer to PDF] is assumed to be nonlinear with increasing marginal cost [figure omitted; refer to PDF] . 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 Section 2 and solved in Section 3. Illustrative examples are presented in Section 4. Section 5 presents a conclusion and future research directions.

2. Model Formulation

Let [figure omitted; refer to PDF] and consider on the planning interval [figure omitted; refer to PDF] a renewable resource, such as a fishery, whose population dynamics are governed by a state equation of the type [figure omitted; refer to PDF] where the state variable [figure omitted; refer to PDF] represents the population biomass at time [figure omitted; refer to PDF] , [figure omitted; refer to PDF] is the logistic growth rate function, and [figure omitted; refer to PDF] is the harvesting rate at time [figure omitted; refer to PDF] .

We use the logistic growth rate function [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is the intrinsic growth rate and the nonnegative integer [figure omitted; refer to PDF] is the carrying capacity of the fish population.

We also use the rate of harvest [figure omitted; refer to PDF] based on the constant "catch-per-unit-effort" and usually considered in fishery models [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is the catchability coefficient and the control variable [figure omitted; refer to PDF] is the effort of harvesting at time [figure omitted; refer to PDF] . Thus, the differential equation (13) becomes [figure omitted; refer to PDF] Finally, let [figure omitted; refer to PDF] and consider the prediction interval [figure omitted; refer to PDF] , where [figure omitted; refer to PDF] . To write the objective function for our nonlinear model predictive control model, we denote by [figure omitted; refer to PDF] the constant price per unit biomass and by [figure omitted; refer to PDF] the constant cost of harvesting per unit biomass. Assuming a nonlinear (quadratic) cost (see Clark [22]), we seek to maximize the profit [figure omitted; refer to PDF] where [figure omitted; refer to PDF] represents the salvage value of the ending state.

An NMPC approach is used in the next section to determine the control variable at time [figure omitted; refer to PDF] 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 [figure omitted; refer to PDF] Now divide the time interval [figure omitted; refer to PDF] into [figure omitted; refer to PDF] subintervals of equal length [figure omitted; refer to PDF] ; then use the trapezoid formula for [figure omitted; refer to PDF] intervals. The objective function (17) becomes [figure omitted; refer to PDF] where, as in (18), we have [figure omitted; refer to PDF] In order to calculate the sum that appears in (19), we write the linear approximation of the state variable [figure omitted; refer to PDF] , which in conjunction with the state equation (16) yields [figure omitted; refer to PDF] where [figure omitted; refer to PDF] Now, substitute (21) in (20) to get [figure omitted; refer to PDF] We thus obtain the second term in the right-hand side of objective function (19) as [figure omitted; refer to PDF] Also, the first term in the right-hand side of objective function (19) is given by [figure omitted; refer to PDF] Combining (24) and (25), and after some simple computations, the objective function (19) becomes [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is independent of the control, [figure omitted; refer to PDF] and [figure omitted; refer to PDF] are the [figure omitted; refer to PDF] -tuple vectors defined by [figure omitted; refer to PDF] and [figure omitted; refer to PDF] is a square matrix of order [figure omitted; refer to PDF] : [figure omitted; refer to PDF] Clearly, its inverse [figure omitted; refer to PDF] exists (since [figure omitted; refer to PDF] and [figure omitted; refer to PDF] ) and [figure omitted; refer to PDF] Note that the matrix [figure omitted; refer to PDF] is negative definite. The global mmaximum of [figure omitted; refer to PDF] is reached at [figure omitted; refer to PDF] such that [figure omitted; refer to PDF] In receding horizon, we obtain [figure omitted; refer to PDF] using the formula [figure omitted; refer to PDF] and hence [figure omitted; refer to PDF] First Term in (32). Consider [figure omitted; refer to PDF] Second Term in (32) . Let [figure omitted; refer to PDF] and [figure omitted; refer to PDF] : [figure omitted; refer to PDF] Third Term in (32). Consider [figure omitted; refer to PDF] Therefore (32) becomes [figure omitted; refer to PDF] Since the choice of [figure omitted; refer to PDF] was arbitrary in [figure omitted; refer to PDF] we deduce the general relationship between the optimal control [figure omitted; refer to PDF] and the state [figure omitted; refer to PDF] over the whole horizon interval [figure omitted; refer to PDF] as follows: [figure omitted; refer to PDF] where [figure omitted; refer to PDF]

Substituting (37) in the state equation (16) we obtain the following differential equation: [figure omitted; refer to PDF] which will be solved numerically over the interval [figure omitted; refer to PDF] with the initial condition [figure omitted; refer to PDF] to get the solution [figure omitted; refer to PDF] over [figure omitted; refer to PDF] for any [figure omitted; refer to PDF] . Once [figure omitted; refer to PDF] is found we replace its value in [figure omitted; refer to PDF] to find the optimal value of the harvesting effort.

Since the harvesting effort [figure omitted; refer to PDF] 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 [figure omitted; refer to PDF] : [figure omitted; refer to PDF]

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 horizon [figure omitted; refer to PDF] , we take for all the following simulations an instant time [figure omitted; refer to PDF] and a prediction horizon [figure omitted; refer to PDF] .

Case 1.

Consider the following values of parameters: [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] . 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 [figure omitted; refer to PDF] . 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 [figure omitted; refer to PDF] is evaluated by substituting the expression of the optimal control (30) in (26) and we obtain [figure omitted; refer to PDF]

Figure 1: Variations of [figure omitted; refer to PDF] and [figure omitted; refer to PDF] as functions of [figure omitted; refer to PDF] on [figure omitted; refer to PDF] in Case 1.

[figure omitted; refer to PDF]

Case 2.

We keep all the parameters constant as in Case 1, except for the growth coefficient which is taken to be [figure omitted; refer to PDF] instead of [figure omitted; refer to PDF] . 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 [figure omitted; refer to PDF] and remains 0 on the interval [figure omitted; refer to PDF] . The reason the harvesting effort is 0 is that (37) is negative on this interval.

Figure 2: Variations of [figure omitted; refer to PDF] and [figure omitted; refer to PDF] as functions of [figure omitted; refer to PDF] on [figure omitted; refer to PDF] in Case 2.

[figure omitted; refer to PDF]

Case 3.

Again, we keep all the parameters constants as in Case 1, except for the initial biomass which is taken to be [figure omitted; refer to PDF] instead of [figure omitted; refer to PDF] . 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 [figure omitted; refer to PDF] .

Figure 3: Variations of [figure omitted; refer to PDF] and [figure omitted; refer to PDF] as functions of [figure omitted; refer to PDF] on [figure omitted; refer to PDF] in Case 3.

[figure omitted; refer to PDF]

Case 4.

Once again, we keep all the parameters constants as in Case 1, except for the growth coefficient which is taken to be [figure omitted; refer to PDF] instead of [figure omitted; refer to PDF] and for the constant price per unit biomass [figure omitted; refer to PDF] instead of [figure omitted; refer to PDF] . 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 of [figure omitted; refer to PDF] and [figure omitted; refer to PDF] as functions of [figure omitted; refer to PDF] on [figure omitted; refer to PDF] in Case 4.

[figure omitted; refer to PDF]

5. Conclusion

Optimal control theory has been plentifully used to determine at time [figure omitted; refer to PDF] the optimal harvesting effort on a planning interval [figure omitted; refer to PDF] . In contrast, starting at time [figure omitted; refer to PDF] with the current state [figure omitted; refer to PDF] , NMPC determines the optimal harvesting effort on the prediction interval [figure omitted; refer to PDF] . Then at time [figure omitted; refer to PDF] with the current state [figure omitted; refer to PDF] , the optimal harvesting effort is determined on the prediction interval [figure omitted; refer to PDF] . This process is repeated over and over until time [figure omitted; refer to PDF] . The determination of the optimal control on the time intervals [figure omitted; refer to PDF] 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 [figure omitted; refer to PDF] 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.

Acknowledgments

The authors would like to thank the referees for carefully reading the paper and making relevant suggestions to improve its final version. The authors extend their appreciations to the Deanship of Scientific Research at King Saud University for funding the work through the research group Project no. RGP-024.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

References

[1] C. W. Clark The Worldwide Crisis in Fisheries: Economic Modelsand Human Behavior , Cambridge University Press, Cambridge, UK, 2006.

[2] H. S. Gordon, "The economic theory of a common-property resource: the fishery," Journal of Political Economy , vol. 62, no. 2, pp. 124-142, 1954.

[3] M. B. Schaefer, "Some aspects of the dynamics of population important to the management of commercial marine fisheries," Bulletin of the Inter-American Tropical Tuna Commision , vol. 1, no. 2, pp. 25-56, 1954.

[4] S. Ganguly, K. S. Chaudhuri, "Regulation of a single-species fishery by taxation," Ecological Modelling , vol. 82, no. 1, pp. 51-60, 1995.

[5] M. Fan, K. Wang, "Optimal harvesting policy for single population with periodic coefficients," Mathematical Biosciences , vol. 152, no. 2, pp. 165-177, 1998.

[6] T. Peng, "A bioeconomic model of continuous harvesting for a single species fishery," International Journal of Pure and Applied Mathematics , vol. 48, no. 1, pp. 21-31, 2008.

[7] H. R. Joshi, G. E. Herrera, S. Lenhart, M. G. Neubert, "Optimal dynamic harvest of a mobile renewable resource," Natural Resource Modeling , vol. 22, no. 2, pp. 322-343, 2009.

[8] E. Braverman, L. Braverman, "Optimal harvesting of diffusive models in a nonhomogeneous environment," Nonlinear Analysis. Theory, Methods & Applications , vol. 71, no. 12, pp. e2173-e2181, 2009.

[9] G. E. Halkos, G. Papageorgiou, "Dynamic optimization in natural resources management," Journal of Environment Management and Tourism , vol. 2, no. 2, pp. 92-97, 2010.

[10] S. Duncan, C. Hepburn, A. Papachristodoulou, "Optimal harvesting of fish stocks under a time-varying discount rate," Journal of Theoretical Biology , vol. 269, no. 1, pp. 166-173, 2011.

[11] M. Ellis, H. Durand, P. D. Christofides, "A tutorial review of economic model predictive control methods," Journal of Process Control , vol. 24, no. 8, pp. 1156-1178, 2014.

[12] D. Simon http://liu.divaportal.org/smash/get/diva2:690771/FULLTEXT01.pdf Model predictive control in flight control design stability and reference tracking [Ph.D. thesis] , Department of Electrical Engineering, Linköping University, Linköping, Sweden, 2014.

[13] G. C. Goodwin, S. F. Graebe, M. E. Salgado Control System Design , Prentice Hall, Upper Saddle River, NJ, USA, 2001.

[14] S. J. Qin, T. A. Badgwell, "A survey of industrial model predictive control technology," Control Engineering Practice , vol. 11, no. 7, pp. 733-764, 2003.

[15] L. N. Petersen, J. B. Jorgensen, "Real-time economic optimization for a fermentation process using model predictive control," in Proceedings of the European Control Conference (ECC '14), pp. 1831-1836, IEEE, Strasbourg, France, June 2014.

[16] S. del Favero, D. Bruttomesso, F. Di Palma, G. Lanzola, R. Visentin, A. Filippi, R. Scotton, C. Toffanin, M. Messori, S. Scarpellini, P. Keith-Hynes, B. P. Kovatchev, J. H. DeVries, E. Renard, L. Magni, A. Avogaro, C. Cobelli, "First use of model predictive control in outpatient wearable artificial pancreas," Diabetes Care , vol. 37, no. 5, pp. 1212-1215, 2014.

[17] L. Sun, J. D. Hedengren, R. W. Beard, "Optimal trajectory generation using model predictive control for aerially towed cable systems," Journal of Guidance, Control, and Dynamics , vol. 37, no. 2, pp. 525-539, 2014.

[18] P. Karamanakos, T. Geyer, S. Manias, "Direct voltage control of DC-DC boost converters using enumeration-based model predictive control," IEEE Transactions on Power Electronics , vol. 29, no. 2, pp. 968-978, 2014.

[19] M. D. Berkemeier, S. Perez, D. Bevly, "On the suitabilityof nonlinear model predictive control for unmanned ground vehicles," in American Control Conference (ACC '14), pp. 4605-4610, Portland, Ore, USA, June 2014.

[20] A. Pytel, S. Kozak, "Modelling and effective predictive control of gas turbine process," in Proceedings of the 15th International Carpathian Control Conference (ICCC '14), pp. 469-474, IEEE, Velke Karlovice, Czech Republic, May 2014.

[21] T. Brechet, C. Camacho, V. M. Veliov, "Model predictive control, the economy, and the issue of global warming," Annals of Operations Research , vol. 220, pp. 25-48, 2014.

[22] C. W. Clark Mathematical Bioeconomics: The Mathematics of Conservation , John Wiley & Sons, Hoboken, NJ, USA, 2010., 3rd.

[23] S. P. Sethi, G. L. Thompson Optimal Control Theory: Applications to Management Science and Economics , Kluwer Academic Publishers, Dodrecht, The Netherlands, 2000., 2nd.