Optimal Fishery with Coastal Catch

In many spatial resource models it is assumed that the agent is able to determine the harvesting activity over the complete spatial domain. However, agents frequently have only access to a resource at particular locations at which the moving biomass, such as fish or game, may be caught or hunted. To analyse this problem, we set up a simple optimal control model of boundary harvesting. Using the Pontryagin’s Maximum Principle we derive the associated canonical system, and numerically compute canonical steady states and optimal time dependent paths, and characterise the optimal control and the associated stock of the resource. Finally, we extend our model to a predator-prey model of the Lotka-Volterra type, and show how the presence of two species enriches the results of our basic model. For both models we illustrate the dependence of the optimal steady states and the optimal paths on the cost parameters.


Introduction
Optimal control theory is an important tool to design optimal harvesting strategies in the management of natural resources. While the early work only considers the temporal dimension, see, e. g., the familiar monographs of Conrad and Clark (1987), Conrad (2010) and Clark (2010), more recent work captures the spatial dimension as well. The early models in spatial resource economics feature discrete patches, where at each location of the resource the stock evolves according to an ordinary differential equation (ODE). Migration of the biomass is then modelled as entry and exit of the biomass from one location to the other. Notably, Sanchirico and Wilen (1999) set up a bioeconomic model with a finite number of patches with migration of the biomass and reallocation of effort between patches, and demonstrate how biological and economic features jointly determine the equilibrium distribution of the biomass and the harvesting effort over space and time.
However, in many cases the continuous process of migration is more adequately described by partial differential equations (PDEs) characterising the spread or diffusion of the resource in the domain. 1 Thus, a branch of spatial resource economics has developed, which originates from theoretical biology and applied mathematics and makes use of PDEs to model the migration of biomass, such as fish or wildlife, and which becomes increasingly popular. For example, Bai and Wang (2005) consider a spatially distributed species with movement, characterised by a spatially non-homogeneous Gilpin-Ayala diffusive equation (with Neumann boundary condition), and investigate the optimal harvesting policy with spatially variable harvesting effort. In this way they are able to generalise the spatially homogeneous approach (see, e. g., Fan and Wang, 1998;Clark, 2010). Similarly, characterisations of optimal controls are established for many models so that control theory can now readily applied in these cases when the state variables (stocks of resources) are governed by PDEs. For example, necessary conditions of the maximum principle for infinite horizon intertemporal optimization problems with diffusion can be found in Brock and Xepapadeas (2008). Recently, Bressan et al. (2013) consider a resource problem where harvesting activity must be allocated over a bounded connected domain, while at each time total harvesting activity is bounded by a capacity constraint; assuming that profit is linear in the harvesting rate, these authors demonstrate the existence of a constant profit-maximizing harvesting rate.
In parallel with this, comparative statics, that is analyses of the effects of parameter values on state variables, optimal controls and the maximised objective function, have made a substantial progress in recent years. For example, Montero (2001), utilizing the optimal harvesting model previously studied by Cañada et al. (1998) and Montero (2000), 2 investigates the effects of an increase in the spatial domain on the maximised payoff. Montero shows that the effect of the domain depends on the type of the boundary condition: for a Dirichlet boundary conditions-which models a lethal frontier or hostile shore-the optimal benefit increases as the domain becomes larger. However, under a Neumann boundary condition-which models that the species cannot pass the boundary-this conclusion does not necessarily hold. It is therefore crucial which type of a boundary condition applies to appropriately describe the bioeconomic situation.
A major field of application of control theory in this context are the benefits and costs of marine reserves (regions of no-harvest) in the management of fisheries. For example, in a discrete space model Sanchirico and Wilen (2005) demonstrate that the creation of marine reserves may constitute an economically optimal policy. Analogously, Neubert (2003) considers a resource (with logistic growth) on a continuous, finite one-dimensional spatial domain with continuous diffusion. Assuming Dirichlet boundary conditions, he solves for the spatial distribution of fishing effort that maximizes the steady state yield, and shows that the establishment of marine reserves is always part of an optimal harvesting policy, with the size of the reserve depending on the size of the spatial domain. Ding and Lenhart (2009) extend Neubert's model to a multidimensional spatial domain, and show that also in this case marine reserves may continue to be part of the optimal strategy of the management of fisheries. Joshi et al. (2009) extend the work of Neubert (2003) and Ding and Lenhart (2009) and characterize the spatio-temporal distribution of harvest effort which maximizes the present value of harvest (in the absence of harvesting cost).
Spatial bioeconomic models are also suitable to describe the behaviour of interacting species such as inhibitor-activator or predator-prey systems which allow for movement (diffusion) of the species. For instance, Leung (1995) investigates spatially distributed profitmaximizing harvesting policies in a predator-prey model with Lotka-Volterra type of interaction and diffusion. 3 Subsequently, optimal control problems in a diffusive predator-prey system are studied by Fister (1997Fister ( , 2001, Fister and Lenhart (2006) and Chang and Wei (2012). 4 While Fister (1997) considers harvesting of both the prey and the predator population (on the full spatial domain), the economic agent in Fister (2001) is interested in harvesting the predator only and merely appreciates the existing stock of the prey. In the latter model the author introduces a boundary control which allows for controlling the migration of the populations across the boundary (e. g., a fence, mesh size of a net, filter). 5 Brock andXepapadeas (2005, 2008) examine conditions for the emergence of spatial heterogeneity of the resource and (harvesting) effort concentration. Based on the linearisation around homogeneous states, they find a mechanism capable of creating spatial heterogeneity, which they call optimal diffusion induced instability, and which is similar to the so-called Turing mechanism in uncontrolled systems (Turing, 1952;Murray, 2003). Subsequently, Xepapadeas (2010) explores this instability in a model where optimizing agents generate optimal agglomerations, and Brock and Xepapadeas (2010) study interacting economic-ecological systems with diffusion of two interacting resources with inhibitor-activator characteristics 3 The standard predator-prey model with Lotka-Volterra type of interaction can be found, for example, in Clark (2010, Section 8.2). For a brief survey of single and multi-species population models with diffusion see Okubo and Levin (2001) and Murray (2003). 4 The corresponding simpler case of only one species and one control is studied by Stojanovic (1991); Leung and Stojanovic (1993); He et al. (1994He et al. ( , 1995. 5 Belyakov and Veliov (2014) also investigate an optimal harvesting problem with an age-structured population of fish (but without spatial dimension) and selective fishing where only fish of prescribed size is harvested; a similar problem is dealt with by Quaas et al. (2013).
(reaction-diffusion system). This work has been extended by Grass and Uecker (2015) and Uecker (2016) by following the linear instabilities to compute numerical bifurcation diagrams of heterogeneous steady states, identifying potentially optimal ones, and, moreover, computing the optimal controls leading to such patterned optimal states. Both, Brock and Xepapadeas (2010) and Uecker (2016), also compare the competitive equilibrium, where a myopic economic agent is located at each point in space, with the social optimum, where a foresighted social planner has control over all economic activities, find that the social optimum may be significantly higher than the competitive equilibrium, and identify the co-states (or shadow prices) of the social optimum as optimal taxes in the competitive model.
Here we set up and analyse yet another class of diffusive optimal control problems: the case where the agent only has access to a resource at particular locations at which the moving biomass, such as fish, may be caught or harvested. This limited access may be the consequence of geographic or physical constraints that render a widely spread harvesting activity utterly impossible, or may result from legal restrictions, such as the presence of reserve areas or of external property rights, that prohibit harvesting outside designated areas. (For example, fishing may be restricted to certain places at the shore of a lake.) We begin our analysis of this type of a problem by setting up a simple one-species infinitetime-horizon optimal control problem for a fishery model with diffusion and boundary catch, and then extend it to a two species predator-prey model. Due to boundary harvesting in our approach, the no-flux or Neumann boundary conditions of previous models need to be modified to take into account the take-out of fish precisely at the boundary, i. e., at the coast; this renders the boundary condition at the point of take out of the Robin type.
For the scalar model in a one dimensional space, some first conclusions about candidates for optimal steady states, i. e., optimal steady harvesting rates at the boundary and corresponding distributions of fish in the domain, can be obtained from a phase plane analysis of the state equation. However, to study optimal time-dependent paths, and to extend the analysis to predator-prey models or other generalizations, we need to resort to numerical methods. Thus, we use the Pontryagin maximum principle to derive the so called canonical systems for both models, i. e., the necessary first order optimality conditions. The associated stationary problems are systems of non-linear elliptic PDEs where the controls determine the boundary conditions, and in order to compute their solutions, i. e., the so called canonical steady states, we use the continuation and bifurcation software pde2path (see Uecker et al., 2014). In a second step, using the add-on package p2poc (see Uecker, 2015) we identify the optimal steady states and compute their canonical paths, similar to the applications in Grass and Uecker (2015) and Uecker (2016).
The method can easily be generalized to more complicated boundary control problems, and as an example we also consider a two-species reaction-diffusion system of Lotka-Volterra type. It turns out that for both classes of models, i. e., the scalar one and the interacting species model, a "moderate harvesting policy" is optimal, see sections 3 and 4 for the precise results. This is quite intuitive as excessive fishing leads to a drastic diminution of the stock and thus impairs the conditions for future yield; while a moderate fishing activity forgoes present profits, but saves some of the stock for later catch and growth. Thus, in both the one-species and the two-species case a balanced fishing path is optimal. We complement our analyses by computing the evolution of the respective shadow prices and their spatial distributions in order to enhance the understanding for our findings. Finally, we show how our results depend on the cost of fishing effort: while the qualitative results are quite robust over a large domain of parameter values, polar specifications may seriously, though still in an intuitive way, affect the qualitative features.
Beyond the specific insights into optimal boundary harvesting achieved here, our analysis also makes a methodological contribution: while, according to our knowledge, problems of optimal boundary control with PDE constraints have not yet been utilised in the economic literature, our analysis demonstrates that economic problems of this type actually existand may be solved. While the presence of boundary control problems in resource and environmental economics is apparent from this paper, formally similar problems arise in other areas of economics. Possible examples are the control of the proliferation of an infectious disease (epidemic) as a public health protection policy; the control of the flow of information in problems of advertising, information management or crowd motion and herding; the launch of new products when demand depends on consumer (the diffusion of) experience etc. Thus, this paper should also help advance the consideration of boundary control methods in other areas of economics.
The rest of the paper is organized as follows. In Section 2 we set up our basic one-species fishery model with boundary catch, provide a phase plane analysis of the constraints for the stationary case and derive the canonical system. In Section 3 we describe our numerical method, and present our results for the one species model; while the results of the two-species model are presented in Section 4. Finally, in Section 5 we conclude by discussing our results and indicating possible directions for future research.

A Basic Model
Consider a fishery problem where harvesting (fishing) can be done on the boundary of an area populated by some species of fish. For example, think of a fisherman catching fish from the shore. For ease of tractability, we consider a one-dimensional space represented by the interval Ω := (0, l x ). Fishing can be done at location x = 0-the position of the fishermanonly. Let v = v(x, t) be the biomass of fish at location x ∈ Ω at time t ∈ T ≡ [0, ∞). The catch depends on the available biomass of fish v and on the harvesting effort k of the fisherman. We specify the catch (or harvest) as a standard Cobb-Douglas function, The fisher is interested in maximizing the profit from his fishing activity. Let p > 0 denote the market price of one unit of fish, and c > 0 the (constant) per unit cost of harvesting effort. Since fish is a non-durable good, the catch is offered at the market immediately when it is realised. Thus, we model the instantaneous profit from harvesting as The total growth of the stock at a given location x is governed by net growth of the biomass and movement of fish. The growth of the biomass (net of mortality) is assumed to follow bistable growth, (2) f is initially convex and then concave with β > 0 representing the minimum viable population; thus, the growth function exhibits critical depensation (see Conrad andClark 1987, p. 63 andDa Lara andDoyen 2008, p. 18f). The movement of fish is modelled as diffusion, i. e., by a term proportional to ∆v, where ∆ denotes the Laplace operator with respect to x. 6 Thus, the biomass of fish evolves according to the following system of differential equations and boundary conditions (BC) where D is the diffusion coefficient, n denotes the exterior normal to the boundary ∂Ω, and the normal derivative on the left-hand side of equations (3b) and (3c) is defined as Equation (3a) describes the total change in biomass resulting from autonomous growth and diffusion. It states that the rate of change of the state variable, i. e., the concentration of the resource, at a given point in space, is determined by the growth function f , which represents the biological growth process of the resource when left unimpaired, and by the movement (or dispersion) of the stock described by the term D∆v, which reflects the standard assumption that the flux of the stock is supposed to be proportional to the gradient of the size of the stock, with the movement of the stock taking place from places of high towards those of low concentration. In the one-dimensional case considered here we may assume that D = 1, 7 but for conceptual clarity and generalizations we keep D, for the moment. The fact that harvesting takes place at the left boundary of Ω motivates equation (3b), which represents the flux boundary condition at the left. It captures the fact that there is a negative input at x = 0: the negative of the value of the spatial derivative of the stock at x = 0, e. g., g > 0 amounts to less fish at x = 0 than at some x > 0 close to x = 0. Since this spatial differential is induced by the harvest h, we specify g as for some γ > 0. Hence, equation (3b) represents the idea that a larger harvest at x = 0, i. e., a larger take out of fish at the left, increases the differential in stocks between x = 0 and its (right) neighbourhood x > 0. Correspondingly, equation (3c) captures the idea that there is no fishing at the right boundary, and it thus represents the zero-flux (or Neumann) boundary condition at the right.
In equation (3e), γ can be thought of as the inverse of the replacement flux of fish: when γ is high (low) a given amount of fishing leads to a large (small) differential in stocks near to the location of the fisher, and this differential is due to a slow (fast) replacement of fish due to diffusion. Thus, in a simple setting γ could be chosen proportional to 1/D, but we keep it as an independent parameter, essentially to possibly model special conditions for the replacement fluxes at the boundary.
Finally, given the instantaneous profit J c , the agent seeks to maximize the total discounted profits We thus have an optimal control (OC) problem with the partial differential equation (PDE) constraints (3a)-(3e) and the boundary control k.

Phase plane analysis of the constraint for the stationary case
To get an intuition for the constraints in (3), we first briefly sketch a phase plane analysis of the constraint in the stationary case, i. e., for where we set D = 1, and as a shorthand notation write v ≡ ∂ x v and v ≡ ∂ 2 x v. This is a Hamiltonian system with conserved energy = 0, and hence the orbits are level lines of E, see Figure 1(a). In particular, E has a saddle-point at 0), see the outermost curve in Figure 1(b). This homoclinic decomposes the phase plane into an inner region containing the minimum V m = (β, 0) of E and periodic orbits V per , and an outer region, which does not contain bounded solutions. The linearisation of (6) at V m is which yields the eigenvalues λ 1,2 = ± (β − 1)β = ±i β − β 2 . For example, for β = 0.6 this yields P min = 2π/ Im(λ 1 ) ≈ 12.83 as the infimum period of the periodic solutions. Moreover, the period P (V per ) grows with the amplitude of V per and P (V per ) → ∞ as V per → V hom . In Figure 1(b) and (c) we plot a number of initial conditions (v (0), v (0)) that yield admissible orbits, i. e., orbits with v (l x ) = 0. Since we are later interested in the case of harvesting at x = 0 we restrict our attention to v (0) > 0. For l x = 20 these orbits come in two continuous families S 1 (black) and S 2 (blue). The orbits from S 1 are monotonous in v and are close to the stable manifold of (1, 0), i. e., the homoclinic V hom , while those from S 2 make a "loop" around V = (β, 0). the orbits of (6) are the level lines of E. (b) Initial conditions leading to orbits fulfilling the right BC v (l x ) = 0, l x = 20, and associated orbits; the black dots yield monotonous v(·), while the blue dots give orbits with a loop around v = β. These initial conditions were obtained from integrating (6) backward in "time" Thus we have found two families S 1,2 of steady states fulfilling (3a) and (3c), and for reasonable choices of parameters p, c and γ we may expect that for all v ∈ S := S 1 ∪ S 2 there exists an effort level k > 0 such that the boundary condition (3b) on the control is also fulfilled. Now, the crucial question is: which of these steady states are optimal, i. e., maximize J c in S? Moreover, we are interested in the solution of the intertemporal OC problem (3). In order to characterise the optimal fishing effort, we next derive the necessary first order optimality conditions, known as the canonical system.

Derivation and discussion of the canonical system
The so called canonical system (CS) formalism, also known as Pontryagin's Maximum Principle, 8 yields first order necessary optimality conditions. For our case where harvesting effort k is a boundary control, we follow Tröltzsch (2010, Section 3.1) and consider the Lagrangian where λ : Ω × T → R is the Lagrange multiplier for the PDE constraint, also called co-state, which may be interpreted as the shadow price of the biomass at location x at time t. Using integration by parts in x we have and integration by parts in t yields where, since we may restrict our analysis to bounded state and co-state variables v and λ, we used the so-called transversality condition 9 Evaluating the stock of fish by its shadow price λ, equation (11) specifies that the present value of the existing biomass living in Ω converges to zero as we consider the very distant future. Thus, using the BCs (3b) and (3c), i. e., D∂ n v| x=0 = −g and D∂ n v| x=lx = 0, we obtain The first variation of L with respect to v, applied to a test-function φ ∈ C ∞ (Q) with φ(·, 0) = 0, yields Therefore, by density of C ∞ 0 (Q) in L 2 (Q), Q = Ω×T , and by density of ∂ n C ∞ (Ω) in L 2 (∂Ω), the condition ∂ v Lφ = 0 yields ρλ − ∂ t λ − ∆λ − ∂ v f (v) = 0 and the boundary conditions D∂ n λ − ∂ v J c + λ∂ v g = 0 and D∂ n λ| x=lx = 0. Thus, the CS is where k is obtained from k(t) = argmax k L(v(·, t), λ(·, t), k). In the absence of control constraints, 10 the condition ∂ k L = 0 yields System (13) summarizes the necessary first order optimality conditions for the optimal control problem (3). In particular, by equation (13e) the optimal effort is determined so that the value of the marginal product of effort equals its cost. This is because the term ∂ k h = (1−α) (v/k) α represents the marginal product of effort in harvesting, with 1−α being the productivity of effort (or more precisely, the elasticity of the marginal product of effort in harvesting), and p − γλ represents the total value of one unit of biomass (fish) caught.
In the first place, this value equals its market price p, however since each unit can only be caught once, future opportunities are impaired by any take out. More precisely, the extent to which a take out affects future catch depends on the influx or the replacement of fish at the coast (at the boundary), measured by γ. If γ is large, the replacement flux of fish is low, so that the recovery of the stock takes time; on the contrary, if γ is low (and possibly zero in the limit) the replacement rate is high, so that the stock recovers quickly (and possibly instantaneously). Since λ equals the shadow price of fish, the term γλ represents the future reduction in the benefit of the catch due to today's take out of one unit of the biomass-and this value must be subtracted from the market price of fish. We want to solve system (13) for k on the infinite time horizon t ∈ [0, ∞), and thus at first might want to think of (13) as an initial value problem. However, equation (13a) provides initial data for only half the variables, while equation (13b) represents backward diffusion, implying that (13) is not an initial value problem. Instead, we proceed similar to Grass et al. (2008, Chapter 7), Grass and Uecker (2015) and Uecker (2015Uecker ( , 2016: 11 Letting u ≡ (v, λ), we write (13) as where B(u) = (B 1 (u), B 2 (u)) = 0 encodes the boundary conditions at the left and right boundaries, respectively, and where we generally suppress the dependence on parameters.
Next we restrict condition (11) further to lim t→∞ u(t) =û, whereû is a canonical steady state (CSS), i. e., a steady state of the canonical system, and thus a solution of In a first step we then numerically compute the CSSs. Some first hints for candidates of an CSS can be obtained from the phase plane analysis of the stationary state, equation (4) in section 2.1. However, this did not take into account the optimization problem, and hence neither the co-state λ nor the left boundary condition B 1 (u) = 0, which involves the control k, which via condition (13e) is a function of u. As (14) is a non-linear elliptic system, at a given set of parameters (ρ, α, p, c, γ, β) we may expect multiple CSSs,û = (v,λ), with (generically) different values J c (v), or J(v) = 1 ρ J c (v), for different CSSs; and the question arises which of these CSSs maximizes J c amongst the CSSs. In a second step, we want to find optimal steady states (OSSs), and their associated optimal paths. That is, given some initial state v 0 , we want to compute canonical paths (CPs) connecting v 0 to some CSSû, and then compare their respective values.

Canonical steady states
A standard method to get numerical insight into the (possibly non-unique) solutions u of a non-linear equation such as (15), and their dependence on parameters ψ, is as follows: one first seeks a solution u 0 (ψ 0 ) at some parameters values ψ 0 . This first step is sometimes easy, for instance if there is a trivial solution u = 0, at least for some particular parameter values, and sometimes hard. In the latter case, one typically uses a guess u 0 (ψ 0 ), and then tries to improve this guess iteratively to a solution u(ψ 0 ), for instance by the so called Newton method. Then, given a first solution u(ψ 0 ), there are standard methods of so called numerical continuation (see Keller, 1977 andDoedel, 2007) to continue the solution in one or more parameters, i. e., to find a nearby solution u(ψ 0 + δ), where δ is a small change of parameters. This can be repeated to obtain continuous branches of solutions. On these branches there may be bifurcation points, at which other branches of (qualitatively different) solutions bifurcate, and there are numerical methods to detect these bifurcations and compute the bifurcating branches.
pde2path is a Matlab package to do such a numerical continuation and bifurcation analysis for PDEs of type (14) (and also for more general PDEs) over one-, two-and threedimensional domains. Moreover, in the package p2poc (see Uecker, 2015) it has been linked with the boundary-value-problem solver TOM (see Mazzia and Sgura, 2002) to compute time-dependent canonical paths; and, for instance, Uecker (2016) applies this to a vegetation system with distributed controls, which leads to a four component reaction diffusion system with many different solution branches and associated bifurcations. The situation here turns out to be more simple though, and in particular we do not find bifurcations. Nevertheless, it is still useful to study the parameter dependence of the solutions, i. e., the continuation of the solutions of system (14) (and of the generalization (26) below) in parameters.
In our numerical simulations of (14) we set D = 1 and l x = 20 as in section 2.1, and hence a priori know that for a CSSû = (v,λ),v must belong to either the family S 1 of monotonous steady states, or the family S 2 of steady states with "loops". Thus, given parameter values of (ρ, α, p, c, β, γ), the steady-state of the canonical system (15) selects a discrete set of CSSs, some corresponding to the family S 1 , and some to S 2 . In addition to D = 1 and l x = 20, we now fix (ρ, α, p, β, γ) = (0.01, 0.3, 1, 0.6, 0.5), as a base parameter set for the canonical system (13), and take the costs of fishing effort c ∈ (0, 2), which is the economically the most significant variable, as our primary continuation parameter. Figure 2(a) shows two (obtained by different initial guesses for (v, λ)) branches of CSSs as functions of the effort cost c, i. e., so called bifurcation diagrams (or here rather continuation diagrams). We display the instantaneous profit J c , the biomass of fish v, the co-state variable λ, harvesting effort k and the resulting harvest h along these branches. Figure 2(b) displays example plots of v and λ at the specific parameter values c = 0.1 and c = 2, together with the corresponding numbers for k, h and J c . 12 We begin with the discussion of the monotone solution, the black graphs in Figure 2(a) and the plots in Figure 2(b1). Figure 2(a) shows that an increase in the harvesting cost c renders harvesting effort more and more unattractive: with an increase in c effort k is reduced implying a reduction in harvest h, and thus an increase in the stock v; in parallel, higher cost reduce the value of the stock λ. Thus, as expected, the resulting profit falls with higher values of c. Focusing on the cases c = 0.1 and c = 2, Figure 2(b1) depicts the spatial distribution of v and λ. If c is small (here c = 0.1), harvesting effort is high (h = 0.1285); and since harvesting can only be done at the left boundary, the value of the resource is highest at x = 0 and is monotonously decreasing to zero for more remote populations. For high harvesting cost, harvesting becomes almost unprofitable implying that the value of the resource is close to zero even at the boundary x = 0. Since harvesting effort is low, the take out is marginal, and thus the stock is distributed almost equally on the domain Ω = (0, 20).
The blue branches in (a) and the associated solutions in (b2) show different behaviour. Here, an increase in c induces a (slight) increase in the harvesting effort k and thus a higher catch h. Clearly, this behaviour fosters the reduction in profit resulting from higher cost. 13 Clearly, such an solution is economically implausible and thus irrelevant. This can also be seen by inspection of the co-state variable λ(x), which even becomes negative in parts of the domain, and is yet negative at x = 0 for c = 2, see right part of Figure 2 (b2); moreover, λ(x) no longer approaches zero for large values of x. For these reasons, the second branch of a CSS in clearly not optimal.

The saddle point property, and canonical paths
In Section 3.1 we identified CSSs of system (3), i. e., branches of solutions of the canonical system (15). Now, given an initial state v 0 we want to determine whether or not there exist canonical paths (CPs) connecting v 0 to some CSSû. That is, we are interested in time dependent solutions t → u(t) of (14) such that This is a connecting orbits problem, not an initial value problem, as only the first component u 1 | t=0 =v 0 is fixed, while u 2 = λ| t=0 and hence the control k(0) is free. Thus, different situations may arise: 1. There is a unique CP connecting v 0 to one CSSû.
2. There is a unique CSSû which can be reached from v 0 , but different CPs to do so. 3. Different CSSsû (1) ,û (2) , . . . can be reached from v 0 , and for each target there may be more than one CP; 4. There is no CSSû which can be reached from v 0 . If, given v 0 , there is more than one CP, we can compare the respective values of J for those paths, and decide which one is optimal. Put differently, we can also consider a given CSŜ u and ask from which v 0 it can be reached by a suitably chosen CP. In particular, a CSSû that can be reached from allnearby v 0 and such that the associated CP maximizes J(v 0 , ·) is called a locally stable optimal steady state (OSS), while a CSS which can be reached from all v 0 that are admissible, i. e., here all v 0 ≥ 0 (pointwise), and such that the CPs maximize J(v 0 , ·), is called a globally stable OSS.
In general, all of the alternatives 1.-4. above can occur in a given system, and the domains of attraction of different locally stable OSSs are separated by so called Skiba manifolds; see, for example, Grass et al. (2008) for various ODE applications; and Grass and Uecker (2015) and Uecker (2016) for some PDE examples. However, here the situation turns out to be much more straightforward. Numerically, we proceed as follows. Given the spatial discretization of G(u) = 0 with 2n degrees of freedom, i. e., u ∈ R 2n , (14) becomes a coupled system of 2n ODEs, which with a slight abuse of notation, we again write as Here M ∈ R 2n×2n is the mass matrix of FEM (finite element method) approximation. We choose a truncation time T and approximate (17) by where E s (û) is the stable eigenspace ofû for the linearisation M d dtũ = −∂ u G(û)ũ of (18). At t = 0 we already have the boundary conditions v t=0 = v 0 for the states. Then, in order to obtain a well-defined two point boundary value problem in time we need dim E s (û) = n.
Since the eigenvalues of the linearisation are always symmetric around ρ/2 (see Grass and Uecker, 2015, Appendix A) we always have dim E s (û) ≤ n. The number d(û) = n−dim E s (û) is called the defect ofû, a CSSû with d(û) > 0 is called defective, and if d(û) = 0, thenû has the so called saddle point property (SPP). Clearly, these are the only CSSs such that for general v 0 close toû we may expect a solution for the connecting orbits problem (18a), (18b). See Grass and Uecker (2015) for further comments on the significance of the SPP (19) on the discrete level, and its (mesh-independent) meaning for the canonical system as a PDE; and Uecker (2015) for algorithmic details how to implement (18b), and how to find CPs connecting some v 0 toû by a continuation process in the initial states. The CSSs from branch b1 are the only CSSs with the SPP, and hence are globally stable OSSs, while all other CSSsû have a defect d(û) > 0. Figure 3(a) exemplarily shows a CP starting from v 0 , as given in b2/pt0, connecting to b1/pt0; and Figure 3(b) shows the control k and the instantaneous profit J c (t) along this path, together with the control and profit of the initial CSS. Of course, we could as well start with any v 0 , but the idea of the simulations in Figure 3 is that for some reason, e. g., historic, the fisherman starts in some (suboptimal) CSS at t = 0. The CPs then show exactly how to control the system to govern it to the OSS in an optimal way. The strategy to go from b2/pt0 to b1/pt0 is as follows: initially the effort k is reduced substantially below k 0 , leading to a drastic decrease in instantaneous profit J c . The initially low harvesting activity allows the population v to grow, first near the boundary x = 0, and eventually for all x. This policy allows the stock to recover and then the fisherman to increase effort above k 0 of the CSS b1/pt0. For t > 50 we have essentially reached the target CSS b1/pt0, and as a consequence of controlling the system from the starting CSS to the target CSS the total discounted profit increases by about 8.6%.
Another natural situation is depicted in Figure 3(c,d), where we start from the homogeneous steady state v ≡ 1 of the uncontrolled system (with zero flux BC), which corresponds to the case where before t = 0 there is no fishing at all, and the stock of the biomass has thus reached its maximum. Here the optimal strategy is to begin with intense harvesting and then to monotonously decrease harvesting effort to reach the target value. Due to the fact that the initially unimpaired resource provides excellent conditions for harvesting, this state allows to begin with a harvesting rate exceeding the CSS b1/pt0, and total discounted profit surpasses the profit of the CSS (J = 4.596 compared with J 1 = 4.056).
As already indicated in section 3.1, also in other parameter regimes the CSSsû * with monotonous v are the only CSSs with the SPP and maximize J c amongst the CSSs. Similarly, it turns out thatû * always is the unique globally stable CSS, and the CPs to reachû * from some v 0 are quite similar to the examples in Figure 3.

A predator-prey system
The scalar model of Section 2 can be greatly generalized. Here we consider a standard Lotka-Volterra system for two species: the prey (v 1 ) and the predator (v 2 ) in the form with diffusion constants d j , and self damping parameter of the prey β > 0. This system can more compactly be written as it follows that is the unique steady state of the ODE system d dt v = f (v) in the first quadrant, and is globally stable. Similarly, using the functional it follows that for d 1,2 sufficiently large, V * is the unique steady state of (20) with zero flux BC, and is globally stable, see, e. g., Hastings (1978). Analogous to (3) we consider a boundary fishing problem for the Lotka-Volterra system (20), and introduce J and controls via To apply Pontryagin's Maximum Principle we introduce the co-states λ 1,2 and the Lagrangian Integration by parts in x and t now yields, similar to section 2.2, Then ∂ v L = 0 yields the evolution and the BCs of the co-states (combining with (20), to have it all together) and ∂ k L = 0 yields As above, we first compute canonical steady states, i. e., we start with the stationary version G(u) = 0, u = (v, λ) of system (26), again on the domain Ω = (0, 20). We choose the base parameter set (β, d 1 , d 2 , γ 1 , γ 2 , ρ, α 1 , α 2 , p 1 , p 2 ) = (0.6, 1, 10, 0.1, 0.1, 0.03, 0.4, 0.4, 20, 10), and consider the costs (c 1 , c 2 ) as our continuation parameters, starting with c 1 = c 2 = 0.1. According to specification (28), we assume that the predator species moves faster than the prey, d 2 = 10 > 1 = d 1 , and that the market price of the prey exceeds the price of the predator, p 1 = 20 > 10 = p 2 , so that, disregarding the interaction of both species, the fisher is interested in catching the prey rather than the predator. Yet in view of the interaction of the species, the fisher may consider catching the predator as well in order to "protect" the prey from being eaten by the former.
To find a CSS we use initial guesses of the form with parameters δ ∈ (0, 1) and s close to 1, and some variations of (29). But if such an initial guess of this form yields convergence to a CSS, then (for the base parameters (28)) this convergence always leads to the same CSS. Note that a complete graphical phase plane analysis similar to Section 2.1 is no longer possible, and thus it is not clear how many CSSs may exist for (26). However, given that for d 1,2 sufficiently large and zero flux BCs V * is the unique steady state of the Lotka-Volterra system (20), it appears reasonable to expect that CSSs of system (26) are unique for the given parameter regime. (Note that we actually assumed rather large values for the diffusion parameters, (d 1 , d 2 ) = (1, 10), relative to the size of the domain.) Figure 4 depicts the CSSs for the parameter specification (28), and their dependence on the cost parameters (c 1 , c 2 ). 14 Parts (a) and (b) show relevant quantities at the left boundary as functions of c 1 and c 2 , respectively, while part (c) shows the spatial shape of the CSSs for selected values of (c 1 , c 2 ). As expected, an increase in c i (i = 1, 2) leads to a reduction in effort k i and the associated harvest h i , and thus brings forth a recovery of the stock v i . In addition, an increase in c i also imposes an indirect effect on v 3−i resulting from the interaction between both species. Consider an increase in the effort cost c 1 . Clearly, there is no direct effect of c 1 on the fishing effort k 2 as the cost of this activity as well as the market price p 2 are unaffected. However, since an increase in c 1 results in less fishing effort k 1 and thus in a higher stock v 1 , the living conditions of the predator species improve. Accordingly, the stock v 2 tends to increase, but since v 2 and k 2 are complements in the fishing technology, effort k 2 can be reduced with the catch h 2 still going up. Naturally, with higher cost c 1 the value of the the prey λ 1 falls, but the indirect effect, just explained, renders the value of the predator λ 2 to go up. The direct and the indirect effect of an increase in c 1 are also reflected in the profit terms. As expected, J 1 is a decreasing function of c 1 , but the induced interaction effects between both species make J 2 to increase with c 1 . Since the direct effect dominates, total profit J falls with higher cost.
Due to the biologically asymmetric situation of both species, the effect of an increase in c 2 has somewhat different indirect effects. In this case, higher effort cost c 2 lead to a reduction in effort k 2 and thus to increase in the stock v 2 . But with an increase in the predator population the prey population becomes more threatened-rendering its population to decline; and for that reason the fishing activity k 1 is reduced. In this way, since the increase in the harvesting cost c 2 acts as a protection of the predator against being fished, the growth of that population exerts a negative effect on the prospects of fishing for the fisherman. Accordingly, with the population of the predator becoming sufficiently large, the value of a unit of this species becomes negative. This explains why λ 2 is decreasing in c 2 and why it becomes (quickly) negative as the stock v 2 rises.
The qualitative structure of the spatial distribution of both species (and of the associated shadow prices) is quite robust with respect to changes in the costs c 1 , c 2 . Inspecting Figures 4(c1), (c2) and (c3), we infer that by catching the predator species the fisherman makes sure that its stock is kept low at the coast (left boundary) so as to safeguard the prey there. In fact, the stock of the prey reaches its maximum close to the coast, but harvesting at the coast causes the stock to decrease drastically there (unless c 1 is very high, see Figure 4(c2)). In any case, fish located close to the coast is more worthy than at more distant locations, where it is inaccessible for the fisherman, i. e., λ 1 and λ 2 are both decreasing in x. There is one exception, though: when, as explained above, c 2 is sufficiently large, such that it is very expensive to catch the predator species, this stock may swell until it interferes with the prospects of the fisherman to catch the prey. In this case, the value of the predator is negative, λ 2 < 0, and because the damage caused by the predator species is the larger the closer it gets to the shore, λ 2 has its minimum directly at the shore.
In Figure 5 we illustrate the transition dynamics from the unique steady state (22) (with no fishing) to the CSS. Setting (c 1 , c 2 ) = (0.1, 0.1) in Figure 5(a) and (c 1 , c 2 ) = (20, 0.1) in Figure 5(b), these figures can directly be compared with Figure 4(c1) and Figure 4(c2), respectively. When the fishing cost of both species are low, (c 1 , c 2 ) = (0.1, 0.1), the transition to the CSS is accomplished by extensive fishing of both species at an initial phase, with fishing intensities decreasing from high towards low values. Thus, there is a some initial overfishing (because immediate profits are more desirable than those in the future due to the discounting), followed by a recovery phase during which the stocks of the CSS are reached from below, and during which the fisher increases the fishing intensity to that of the CSS. 15 If fishing the prey is very costly, (c 1 , c 2 ) = (10, 0.1), the CSS is characterised by a low effort level k 1 and hence by a high stock of the prey, cf. Fig. 4(c2) and 5(b). Similarly, along the CP leading to the CSS the main harvest is on the predator, with fishing effort being reduced over time, leading to the gradual increase of the prey at the left boundary-and then on the complete domain. However, some harvesting activity on the prey still takes place. Finally, for (c 1 , c 2 ) = (0.1, 10), the roles are basically reversed, subject to the indirect effects resulting from the implicit protection of the predator species by a high effort cost c 2 , as explained above.

Discussion and extensions
To the best of our knowledge, in the context of economics this is the first detailed numerical analysis of infinite time horizon optimal control problems with PDE constraints and a boundary control. We set up a one-species and a two-species fishery model, compute the respective canonical steady states, single out the optimal steady states and explore canonical paths connecting some arbitrary initial state to the optimal steady state. In both cases, the results appear natural and intuitive, and they are remarkably robust with respect to changes in the cost parameters in both optimal control problems.
In particular, we illustrate how the optimal fishing strategy acknowledges the dynamics of growth and movement of fish in a bounded domain (e. g., a lake). It is a crucial feature of our model that fishing is restricted to take place at the boundary of the domain (e. g., the shore) only. The optimal policy then compromises between immediate and future yield taking into account that a higher stock left in the sea may favour future growth of the resource, while a larger take out of fish may initiate the influx of new fish, due to movement taking place from highly to sparsely populated areas.
Things become slightly more intricate in a two-species Lotka-Volterra model. In this case, the interaction between both species, the prey and the predator, has to be taken into account as well. This interaction provides an additional incentive to catch the predator species in order to protect the prey for the purpose of own take out. In particular, we illustrate how this policy depends on the relative effort cost of fishing either species. Also, for both models we diagrammatically explore the optimal policy paths starting from some suboptimal state and connecting to the optimal steady state (canonical paths).
Our approach may be extended to more complicated and realistic models, presumably yielding more complicated and less intuitive results. One obvious generalization would be to generalise the one-dimensional problem by, for instance, considering either spatially dependent (PDE) coefficients, or more complicated multi-species fishery models. A second class of generalizations would be to extend the spatial domain to two dimensions. The spatial domain could either take the form of a shallow lake (modelled as a two-dimensional domain), or of the transversal positions in a shallow river (with two horizontal directions), or of the depth dependence of fishing at a fixed position in a river with one horizontal dimension. In all of these cases, we then have a boundary control function defined on a given part of the boundary, rather than a scalar control. Additionally we may consider advection modelling the transport of fish by flow. Finally, the model (either with one or multiple species) could be used to model competition between two or more agents (fishermen). Even in the basic model of a one-dimensional spatial space, the introduction of a second agent positioned, for example, at the other end of the one-dimensional lake would allow formulating an interesting non-cooperative dynamic harvesting game (viz. a differential game) in fishing activities. Such a model features similarities with, and thus presumably an applicability to, competition models in industrial organization where in the framework of horizontal product differentiation and mobile consumers, demand may grow (or decrease) over time. Therefore, the interaction between two players in the presence of a mobile renewable resource (or non-constant mobile demand) will represent an interesting direction for future research.