Optimal control of harvest timing in discrete population models

Harvest plays an important role in management decisions, from fisheries to pest control. Discrete‐time models enable us to explore the importance of timing of management decisions, including the order of events of particular actions. We derive novel mechanistic models featuring explicit within‐season harvest timing and level. We explore optimization of within‐season harvest level and timing through optimal control of these population models. With a fixed harvest level, harvest timing is taken as the control. Then both harvest timing and level are used as controls. We maximize an objective functional which includes management goals of maximizing yield, maximizing stock, and minimizing costs associated with both harvest intensity and harvest timing. While standard models with compensatory population dynamics predict it is best to harvest as early as possible in the season, we find instances where harvesting later in the season is optimal. Furthermore, we discover interesting oscillations in the population size, which would be unexpected in the model without time‐varying controls.

, Kellogg et al. (1988), Tang et al. (2010), and Xu et al. (2005). However, when trying to identify the best harvest time, optimal control has only been used for continuous-time population models as far as we know. With continuous-time models, all the events representing biological mechanisms are happening at the same time, and thus it is difficult to represent the seasonal structure in the population dynamics, such as recruitment taking place only in certain, relatively short periods of time within a season. The reason for this lack of harvest timing controls in discrete-time population models probably lies in the difficulty of finding examples of discrete-time models taking into account the moment of harvest intervention.
Discrete-time population models give the expected state of a population N t+1 in time t + 1 as a function of the state of the population N t , N N g N t = ( ), = 1,2,…, (1) where g is the per-capita population growth function. Several authors have modified the basic model (1) to build discrete-time harvesting models in two limit situations. More specifically, when the intervention takes places exactly at the beginning or at the end of the period of time separating two consecutive seasons. The idea is simple and we illustrate it with constant effort harvesting. In such a case, a proportion of the population γ (0, 1) ∈ is removed. So, depending on whether the intervention occurs at the beginning or the end of the season one arrives, respectively, at model (3) Models of the types (2) and (3) have been very useful to understand how harvesting affects the stability or the persistence of a population, see, for example, Carmona and Franco (2011), Liz (2010), Liz and Hilker (2014), Liz and Ruiz-Herrera (2012), and Schreiber (2001). However, it is not clear which one of these models, if any, should be used when the intervention occurs neither at the beginning nor at end of the season but somewhere in between. Seno (2008Seno ( , 2010 proposed the only existing discrete-time model we know of that explicitly describes harvest at any specific time in the season, thus overcoming the previous limitation of models (2) and (3). The novel idea of Seno consists of introducing a new parameter θ that measures the percentage of the season elapsed until harvesting takes place. With this new parameter, Seno builds a new population growth function as a convex combination of the population growth functions of (2) and (3). See Section 2.3 for the explicit equation of the model.
Seno's approach has the advantage of providing the perhaps simplest model so far to understand the effect of harvest time in a population. Using it, several authors have studied how changes in the parameter θ affect the population size, stability, or persistence; for example, see Cid et al. (2014), Franco et al. (2017, Liz (2017), and Wang et al. (2017). Moreover, it gives rise to interesting mathematical problems about the global properties of the limit dynamical models that are inherited by the convex combination (Franco et al., , 2020. Nevertheless, Seno's approach has also drawbacks. The model lacks a mechanistic derivation from basic biological principles. Furthermore, as a consequence of this heuristic derivation, the yield associated with the intervention is not defined, making the model unsuitable for the use of optimal control techniques. To resolve this issue, we introduce novel discrete-time harvesting models that take into account the time of intervention but are derived from basic biological principles and have a clearly defined yield. We consider two such models, one for population dynamics of Beverton-Holt type and the other one of Ricker type (Beverton & Holt, 1993;Ricker, 1954). These novel models allow us to study the problem of harvest timing and its level using optimal control techniques in discrete-time systems, which is the goal of this paper. The remainder is organized as follows. In Section 2, we present the novel models and compare them to Seno's approach. Then, in Section 3, we apply them to numerical optimal control techniques considering different scenarios. Finally, in Section 4, we discuss the results obtained.

| VARIABLE HARVEST TIMING MODELS
We will derive (and briefly analyze) two mechanistic models of harvest timing, one a Beverton-Holt-type model and the other a Ricker-type model. The two models differ in the way they account for continuous within-season natural mortality. Other than that, they are based on identical assumptions, which we will outline in the following. Finally, we will briefly compare the two models to the Seno approach.

| Derivation
Consider a single-species population and let N t be the population size measured at time t, where t = 0, 1, 2, … We refer to the time period between the consecutive steps t and t + 1 as season. The season length depends on the ecology of a population. For example, in studying deer it might be reasonable for a time step to take 1 year. However, in studying an organism with a shorter or longer generation time a time step might take 2 weeks or several years. In our model, we assume that the season has unit length. The time-discrete censuses N t measure population sizes between-season. They are determined by processes within-season. We will use the notation n τ ( ) for the within-season population size, where τ [0, 1] ∈ is the time within a season. Figure 1 illustrates the relationships between N , n, t, and τ.
Within-season, we assume that the population undergoes three events: density-dependent mortality continually during the season; proportional harvest at a time θ [0, 1] t ∈ within the season; and a density-independent recruitment pulse. Table 1 gives an overview of the processes within-season and at which time they occur. Here, we assume that population census takes place immediately after recruitment. This may apply to species where breeding individuals or their offspring are easy to observe and measure, for example, when eggs or simply nests can be counted or individuals become more sedentary during breeding. By contrast, other species may be easier to census before reproduction, for example, when individuals are large in body size or do not hide in seclusive breeding sites. In this case, the model formulation would need to change the order of recruitment and census.
Let us now go through a mathematical description of the processes in the course of a season and eventually arrive at a single difference equation that maps the population size from one season to the next. We start with an initial value of the within-season population size that is given by the preceding census, that is, , there is continuous density-dependent mortality, which can be described by the differential equation where β > 0 is a mortality parameter. For the ease of exposition, we ignore densityindependent mortality. The solution of (5) to the initial condition (4) is there is an instantaneous harvesting event with a yield Y t proportional to the current population size where γ (0, 1) ∈ is the proportion of the population that is being harvested. Immediately after the harvest, the population size is In the time period , there is still continuous density-dependent mortality, which is again described by (5). Its solution to the initial condition (7) gives the population size At time τ = 1 − , there is a recruitment pulse, giving a population size n b n (1 ) = (1 ) + − at the immediately following census. Parameter b is the recruitment factor, and we assume that recruitment always increases population size, that is, b > 1. The population size in the next season is given by N n Though the process of mortality occurs continuously, we can represent the population dynamics by the recurrence equation (8). It represents an iterative formula for our mechanistic model with variable harvest timing within the season. Due to its structural similarity with the Beverton-Holt map, we refer to (8) as the Beverton-Holt-type model.

| Some analysis
We will first study the impact of harvest timing on the population size and yield. We will then determine the fixed points and their stability for the case of constant harvest times that do not change between seasons. Finally, we consider the special cases of harvesting at the very beginning and end of the season, to gain some more intuition about the structure of Equation (8).
First, to understand how the timing of harvest affects the population and yield, note that the next-season population size of the model (8) is a strictly decreasing function of θ t , because This means that harvesting at θ = 0 t will result in the highest next-season population size, while harvesting any time later in the season will result in a diminished next-season population size. Therefore, if we were to consider only the management goal of maximizing next-season population size we would always harvest at the beginning of the season, θ = 0 t . This result corroborates existing modeling insights from the literature (Kokko & Lindström, 1998). The explanation is that harvesting individuals that would have later died of natural causes has no direct effect on the population size, but releases density-dependent intraspecific competition in the remaining population. The earlier an individual is harvested, the more the remaining population benefits because of reduced negative density dependence.
The yield is also a strictly decreasing function of harvest timing, because Therefore, to maximize yield we would choose to harvest at the beginning of the season. Any later season harvest would result in a lower yield. The explanation is straightforward. When we harvest a fixed proportion of a population that declines with later harvest timing, we catch more the earlier we harvest. Second, if we set the harvest timing equal to a constant, we can find equilibria for our model. Let θ θ = t , a constant, for all t. Then there are two equilibria, namely, N * = 0 0 corresponding to extinction and The nontrivial equilibrium is always asymptotically stable when it exists, that is, for b > γ 1 1 − . In this case, the trivial equilibrium is unstable. Otherwise, the trivial equilibrium is asymptotically stable, as the recruitment factor is too small to replace the losses due to harvesting.
Third, consider the limit cases of (8), which are for θ 1 t ≡ . Equations (9) and (10) consist of three functions composed in a different order. The three functions are the mortality function f The order in which these functions are composed tells us the order of events for each equation. Equation (9) for harvesting right at the beginning of the season can be written as Clearly, the order of events is the first harvest, then mortality, and finally recruitment. Equation (10) for harvesting at the end of the season can be written as This means that Equation (10) does not distinguish the order between the events of harvesting and recruitment, since they are both assumed to be density-independent and are assumed to occur consecutively. Recruitment may occur before harvest or harvest may occur before recruitment, and the expression would be the same. Mortality must have occurred first. Since the order of recruitment and harvest would affect yield, it is important that we remain consistent with our assumption that recruitment occurs last in the season during the remainder of our analysis. The order of events for Equation (10), then, is mortality, then harvest, and then recruitment.

| Derivation
We now make a slightly different assumption on the within-season mortality that will result in a Ricker-type model. The assumptions on all the other processes are identical to the ones of the Beverton-Holt-type model.
Here, we assume that the within-season mortality can be described by That is, the per-capita mortality is given by μN − t with μ > 0 being a parameter. The value of N t is the population size at the beginning of the season and determines the per-capita mortality, which remains constant for a given season. For example, this assumption may be applicable to species for which the survival of an individual during the season is mainly determined by its ability to acquire food resources at the beginning of the season; in the remainder of the season it will then live from the energy reserve or food stock accumulated at the beginning of the season. For example, think of species that, before the winter season, collect and store food that they eat over winter. The resource availability at the beginning of the season is set by the population size, assuming that the resources are shared evenly between individuals as is typical for scramble competition. Another example is species, like, salmon that cannibalizes on young conspecifics or eggs, since the mortality rate can be assumed proportional to the early-season population size when individuals are young (Brauer & Castillo-Chavez, 2012).
By contrast, Equation (5) of the Beverton-Holt-type model assumes a per-capita mortality that is proportional to the current population size throughout the season. This can be viewed as a form of contest competition, where individuals continually compete for resources. The assumptions behind Equations (11) and (5) are made to be identical with the ones in Kokko and Lindström (1998). Now let us go again through the processes as in Table 1. For the time period τ θ (0 , ) t + − ∈ , we have to solve (11) with the initial condition (4), resulting in Taking the remaining population as the initial condition for the time interval τ θ ( , 1 ) t + − ∈ , applying again the continuous mortality (11), and multiplying by a recruitment factor b > 1 at τ = 1 − finally gives This is identical to the Ricker map with r b γ = ln( (1 − )) and K r μ = ∕ .

| Some analysis
The recurrence equation (13) for the population size is independent of the harvest timing, which is a somewhat unexpected result and in contrast to the Beverton-Holt-type model (8).
The dynamics of Equation (13) is well known. There are two equilibria, namely, N * = 0 0 corresponding to extinction and N K * = 1 . The nontrivial equilibrium can be asymptotically stable or undergo a sequence of period-doubling bifurcations to chaos (May & Oster, 1976). N * 1 is positive for b γ (1 − ) > 1, in which case N * 0 is unstable. If N * 1 does not exist, N * 0 is asymptotically stable. Interestingly, the yield (12) does depend on the harvest timing. It strictly decreases as a function of θ t , because These results can be explained as follows. The population size is independent of the harvest timing because the constant mortality is set by the conditions at the beginning of the season and thus remains density-independent throughout the season. That the yield decreases with later harvest timing is a straightforward effect of catching a fixed proportion of a declining population.

| Seno's approach
Seno's approach accounts for harvest timing in discrete-time models in a heuristic way. It was proposed and first used in Seno (2008Seno ( , 2010. For a population with a per-capita growth function g N ( ) t in the absence of harvesting, the formulation reads harvest after reproduction harvest before reproduction (15) Equation (15) is a convex combination of two order-of-events cases, namely, (a) harvest after reproduction and (b) harvest before reproduction. They correspond to the limit situations (2) and (3). The weight coefficient of the convex combination is interpreted as the within-season harvest time θ [0, 1] ∈ . Though the time of harvest and the proportion of harvest are known, the yield is not specified, a topic on which Seno makes no comment. We know what the population size is at time t and t + 1, but the population size at any intermediate time is unknown. One might speculate that the yield could be approximated by a convex combination of the yields for the two order-of-events cases. Yet, it is also unclear when reproduction occurs. One term assumes reproduction occurs before harvest, while the other describes reproduction occurring after harvest.
In the following, we will compare the population sizes predicted by the Beverton-Holt-type model and the Ricker-type model with the corresponding Seno models. First, consider the percapita population growth g N , which corresponds to Beverton-Holt dynamics. This actually matches our Beverton-Holt-type model (8), as we obtain exactly this expression for θ = 0 t and for θ = 1 t . Now, plugging g N ( ) t into the Seno model (15) and comparing it with the next-season population size of (8), one can show that That is, the population sizes in the Seno model are greater than or equal to the population sizes in the mechanistically derived model. So the Seno model will never underestimate, but could overestimate population sizes.
To study the differences between the Beverton-Holt-type model and the corresponding Seno model numerically, we used Latin Hypercube Sampling (LHS) to explore the multidimensional parameter space (Marino et al., 2008). We built an LHS matrix with 100,000 possible parameter sets. We chose the ranges N . By θ here we mean a constant, that is, θ θ = t for each t. The large upper bound for b is motivated by species, such as gypsy moths (Lymantria dispar), who oviposit 100-1000 eggs per female (Doane & McManus, 1981), and coho salmon (Oncorhynchus kisutch), who have up to several thousand eggs in their spawning nests (Fleming & Gross, 1990). Note that such high recruitment factors are followed by an interval of mortality in our two models. First we assumed each parameter was uniformly distributed on these ranges, and then in our second comparison we used log scale sampling for b. To compare the time-series outputs of the two models, we used to measure the relative differences, where S t indicates a value in the time series given by Seno's model and N t a value in the time series given by our mechanistically derived model. The relative differences for all of the parameter sets are very small, with the overwhelming majority of the values virtually identical to zero (see Table 2). We conclude that though the two models (15) and (8) look like they would act on θ t rather differently, in practice this does not seem to be the case.
Second, consider the per-capita population growth g N be ( ) = t μx − , which corresponds to Ricker dynamics. Plugging this into the Seno model, we obtain population sizes that do depend on θ t in contrast to our mechanistically derived Ricker-type model (13). Now, let us look at the limit case θ = 0 t for the Seno model, which reads This does not match our mechanistic model (13). However, the limit case of the Seno model for − t which is identical to (13). One of the limit cases for Seno's model thus matches the Ricker-type mechanistic derivation. We might conclude that the latter is somehow less general than the corresponding Seno model. Or we might conclude that the Seno model for these growth dynamics adds unnecessary complexity for the biological situation. For more insights into the Seno model with Ricker dynamics, see . We performed a similar analysis also for the Ricker dynamics, using an LHS matrix with again 100,000 possible parameter sets. We chose the ranges and T [5, 25] ∈ . The Seno model differs more from the mechanistic derivation (13) for the Ricker dynamics than for the Beverton-Holt dynamics. But for most of the parameter space explored the differences are still quite small, see Table 2.

| OPTIMAL CONTROL OF HARVEST INTENSITY AND TIMING IN BOTH MODELS
For our mechanistic models we consider possible components which a management goal might have and some simple ways of incorporating them into our objective functional. We do not consider optimal control of the Seno model because the yield is unknown, even though it is often an important factor in management decisions. For a specific application, the objective functional might be quite different from what we consider here. For instance, pest control programs will aim to minimize population size and the invoked costs. In the context of exploitation as in fisheries or hunting, managers might try to maximize population size, using a weighted sum of population size at each time step, where T is the final time. The dynamics of N t can be found in Equation (8) for the Beverton-Holt-type model or (13) for the Ricker-type model. The coefficients A t are weights for the population for times t T < , while nonzero A T indicates population conservation at the final time. In some applications, managers may worry only about the final state by setting A = 0 t for t T < . Note that we will use the notation A to mean the vector (A A A , , …, T

2
). Alternately managers may want to maximize yield, perhaps to maximize revenue, where B t is revenue per unit and Y t denotes the yield at time t. Note that Y t is a function of N t given by (6) or (12). We will use the notation B to mean the vector (B B B , , …, T 0 2 −1 ). We will see in our results some interplay between A T and B. The population N 0 is given, but the yield between t = 0 and 1 is dependent on the control time θ 0 , thus the indices on the sums in (17) and (18) start at different times. Managers may also want to minimize cost in some way. The objective functional we wish to maximize as a function of harvest timing is given by where δ (0, 1] ∈ is a discount factor, which ensures that funds obtained earlier are worth more due to being able to invest them. We assume 2 , which indicates that the cost of harvest is minimized in the middle of the season. We used a harvest cost that is lowest in the middle of the season to illustrate the results of costs that were minimized during a particular time within the season rather than the beginning or the end. Qualitatively similar results can be obtained with harvest cost having a minimum at other intermediate points in the season. For a specific application, a more complicated cost function (such as a cubic) may be required to account for important variations of costs within the season. This could be because of labor seasonality, because of target species behavior, such as migration or seasonal aggregation, or other reasons. For example, in Norway red deer (Cerus elaphus) may be less costly to hunt during late fall and winter because they migrate from remote areas high in the mountainous interior to coastal lowlands which are more accessible to hunters (Loe et al., 2016). During the winter they also occupy a smaller range than during the spring and summer, allowing hunters to increase their efficiency. Also, it may not be best to harvest at θ 0 = because many species take part of the season to mature and harvesting juveniles is undesirable.
There can be an interplay in balancing the three parts of the objective functional (19). The term with A i coefficients will be maximized for the Beverton-Holt-type model at θ 0 = , because population size is decreasing as a function of θ t . For the Ricker-type model, this term is independent of θ t . For both models, the term with B i coefficients will be maximized at θ 0 = because yields for both models are decreasing as a function of θ t . Therefore, if we considered a cost that was minimized at the beginning of the season the optimal control problem would be trivial. If harvest timing is minimized at the beginning or end of the season (θ 0 = or 1) the harvest timing question reduces to an exploration of optimal order of events, as seen in Bodine et al. (2012), Zhong and Lenhart (2013), Kokko and Lindström (1998), Lundberg and Lundberg (2013), Jonzen and Lundberg (1999), and Ratikainen et al. (2008). Note that units of the A t , B t , and C i , i = 1, 2, 3, are adjusted so that the terms in J can be interpreted in units of a currency.
For our two population models and using a finite number of time steps, the states, controls, and the objective functional values are uniformly bounded, and thus it can easily be shown that there exists an optimal control (θ γ *, *) in the set

| Optimal control of harvest timing only in our Beverton-Holt-type model
We will now present numerical results and explain the optimization technique for controlling the harvest timing only. For our optimization scheme to find the maximum of our objective functional with respect to the control of harvest timing with harvest intensity held constant, we use MATLAB functions fmincon and MultiStart. Because fmincon is a local minimizer, we use MultiStart to search our parameter space for the global minimum. We choose 100 start points for MultiStart, which then runs fmincon from each of those start points. Thus MultiStart in conjunction with fmincon can provide substantial coverage of the parameter space and achieve a close approximation of the global optimum. Now, to perform optimal control on harvest timing while holding harvest intensity constant, we find θ* by maximizing θ J γ ( , ) over θ as a control and with γ fixed. One can view θ* as a function of γ. Figure 2 shows how the maximum value of the components of the vector θ* varies over the range of γ. Other parameters are held at the baseline values That is, for intermediate values of the harvest intensity, the optimal timing is at the beginning of the season. Small or large harvest intensities move the optimal timing later into the season. For extremely small or large intensities, the optimal timing converges to 1 2 , where an explanation for this could be the following. On the one hand, very large harvest intensities increase the importance of costs, causing costs to dominate the choice of θ. On the other hand, with very small harvest intensities, the observed dynamics suggest that the effect of timing on yields diminishes faster in harvest intensity than the effect of timing on costs.

| Optimal control of our Beverton-Holt-type model
Now we perform optimal control of our Beverton-Holt-type model with two controls: θ and γ, using direct optimization of the objective functional (19) with MultiStart selecting start points (θ γ , ) to use in fmincon. We now describe our notation for the displayed results and for the comparisons with beginning-and end-of-season harvest. Each time-series graph in the following figures shows the optimal control in black squares and the optimal state of population size in blue circles. Each of the controls is graphed on the x-axis with its position on the x-axis marking the timing of the harvest. For example, if the harvest is set for θ 0 = , then the controls will appear at t T = 0, 1, …, − 1; but if the controls are set for θ 0.5 = , then the markers for the controls will appear at t T = 0.5, 1.5, …, − 0.5. The value for the objective functional at the optimal θ is given in the captions, labeled θ γ J ( *, *). For comparison, the values of the objective functional for θ 0 = and 1 are also given, for which we have chosen the notation γ J 0 ( , *) and γ J 1 ( , *), where  0 = (0, …, 0) T ∈ and  1 = (1, …, 1) T ∈ . All the results below are for the vectors A 1 = and B 1 = . Using a thorough search with MultiStart, there was no evidence of nonuniqueness of optimal controls in these numerical results.
Our set of baseline parameters for these simulations is With these parameter values the without-harvest population equilibrium is 50. With optimal control as in Figure 3a, it is best to harvest always at the beginning of the season, meaning θ 0 * = . Harvest intensity begins below 0.2 in the first season, allowing the population to reach a higher level. Once the population reaches this level, harvesting can be maintained at an intensity of above 0.4 without a drop in population size. Harvest intensity is higher during the final harvest (above 0.6) because the final population size is not weighted higher than intermediate population sizes, so it is best to get a higher final yield. The values of the objective functional are θ γ γ J J 0 ( *, *) = 596.00 = ( , *) and γ J 1 ( , *) = 282.31; that is, harvesting at the beginning of the season gives a 53% higher objective functional value than harvesting at the end of the season. F I G U R E 2 Optimal control of harvest timing for fixed values of the harvest intensity γ. In black we have the relative difference of objective functional compared to harvesting at the beginning of the season, * In Figure 3b we keep the parameters from above except for b = 2, C = 1 1 , and A = 0, meaning recruitment has decreased from the baseline of b = 6 and cost associated with timing has increased from the baseline of C = 0.1 1 . This produces cycles of harvest timing and harvest intensity combinations. Harvest intensity is very high (about 0.9 in the first cycle and about 0.5 in the remaining cycles) for one season and then very low (close to 0) for the following season. On the season with the high intensity, harvest timing is at the beginning of the season. This leads to a drop in population size in the next season. Then for the low-intensity harvest season, harvest timing is at midseason. This leads to an increase in population size in the next season. Overall, this produces cycles in the managed population where the population size drops low, then recovers over the low-harvest midseason event and then drops low again. In this case, the cycle in harvest intensity seems to be the crucial part of the maximization process; this is suggested by the objective functional values for combining exclusively beginning-of-season harvest timing with the same cycle of harvest intensity being similar to that for the cyclic harvest timing. In other words θ γ J ( *, *) differs from γ J 0 ( , *) by only 2.4%. The values of the objective functional are θ γ J ( *, *) = 42.00, γ J 0 ( , *) = 41.00, and γ J 1 ( , *) = 10.71. Figure 3c was generated with the baseline parameters except for C = 50 1 (increased from C = 0.1 1 ). Again, we observe a cyclic harvest pattern: for one season it is best to harvest hardly at all (intensity close to 0) with midseason harvest while for the next season it is best to harvest early in the season with high intensity (about 0.6). This demonstrates that the cycles in population size and harvest control pattern can emerge from increasing C 1 , that is, from harvesting at the season beginning or end being more costly. Interestingly, even though the early-season harvest intensities during the cycle are higher than those for smaller values of C 1 (Figure 3b), the population size oscillates around higher levels than in Figure 3b. This is due to the initial harvest in the first season being nil rather than large.
3.3 | Optimal control of harvest timing and harvest intensity in our Ricker-type model We perform optimal control of both γ and θ using the objective functional (19) with our Rickertype model (13) determining the state values and the corresponding yield (12). We again use direct optimization of the objective functional with fmincon and MultiStart in MATLAB. Our set of baseline parameters for these simulations with the Ricker-type model is These initial parameters emphasize yield and cost but not population size. Later we will explore setting A 1 = . Figure 4a was generated with these baseline parameters. In this case there is a cycle with one high-intensity early-season harvest followed by two very-low-intensity midseason harvests. The high-intensity harvests are all greater than 0.8, with the final harvest having intensity close to 1 causing the population to near extinction. The values for the objective functional are θ γ J ( *, *) = 83.79, γ J 0 ( , *) = 83.64, and γ J 1 ( , *) = 6.60, with θ γ J ( *, *) and γ J 0 ( , *) differing by 0.18%. With the baseline parameters we see cycles in the population which is again caused by timevarying harvest intensities and harvest timings. Note that these cycles occur for parameter values where the population without harvesting would not normally experience cycles, that is, where the Ricker map exhibits a stable equilibrium. Figure 4b was generated with several parameters differing from the baseline: μ = 1, b = 5, N = 1

| DISCUSSION
To investigate the effects of harvest timing, we carefully built two mechanistic models with the features of time-varying level and moment of the harvest. With the order of events being mortality-harvest-mortality-recruitment, the mortality was density-dependent and continuous, while recruitment and harvest were both density-independent and discrete in time. We performed optimal control of two discrete models of harvest timing we developed, a novel application of the theory.
We compared our models with the Seno model, the only discrete-time model of explicit within-season harvest timing in the literature. We found surprisingly little difference between the Beverton-Holt-type mechanistic model and the corresponding Seno model. Though the Seno model always predicts larger next-season size than our model, for large parts of parameter space they provide virtually identical time series. Despite this similarity in outputs, our model offers the benefit of known order of events and, critically, known yield. The Ricker-type mechanistic model was identical to one of the limit cases for the corresponding Seno model.
However, the difference between our Ricker-type model and Seno's corresponding model was nonetheless small over large parts of parameter space. We found that there are some cases where midseason harvest is preferable to beginning-or end-of-season harvest. This is a novel result, based on taking into account harvesting costs. If harvesting costs are ignored, existing theory predicts largest population sizes at the beginning of the season (Kokko & Lindström, 1998). A strong Allee effect population size has been reported to lead to nonmonotonic relationships between equilibrium population size and harvest timing, with population size initially decreasing and later increasing-but it is unclear whether they could grow larger than at the beginning of the season Cid et al. (2014). In general, whether controlling only harvest timing or both harvest timing and harvest intensity, midseason harvest tends to occur when yield and harvest intensity are low, or when cost associated with harvest timing is high. Low yield might result from low-intensity harvesting, or from the depletion of stock due to high-intensity harvesting, or as an intermediate season (or two) in a cycle between high-intensity and high-yield seasons. Indeed, under some circumstances midseason harvest is not only preferable to beginning-or end-of-season harvest, but also the only option that will result in a profit rather than a loss.
When only controlling harvest intensity, we did not find any circumstances under which harvesting midseason was preferable to harvesting at the beginning of the season (see Burton, 2020, pp. 58-65). This is because all the terms in the objective functional are either maximized at θ 0 = or do not depend on θ at all. The lowest difference between the midseason option and the beginning-of-season option came when yield was low due to low recruitment.
When exploring control of both harvest timing and intensity for the Beverton-Holt-type mechanistic model, we found fewer variations in our results than when exploring the Rickertype mechanistic model. In all of the situations explored in controlling both intensity and timing for Beverton-Holt-type model, the final harvest had much higher intensity than intermediate harvests. This was due to the importance placed on yield, because even placing importance on stock conservation resulted in high-intensity final harvest. When no importance was placed on stock conservation, the final harvest depleted the stock.
We observed cycles in two cases, both with high cost placed on harvest timing (see Figure 3b,c). We hypothesize that the cycles are caused by an alternating focus on maximizing yield (via a high harvest intensity early in the season when the population size is large) in one season and on minimizing harvest-related costs (via harvesting close to midseason and as little as possible). As the latter allows the population to recover, there is sufficient stock again in the next season to generate revenue from intense harvesting, and the cycle repeats. Hence, the oscillations may be based on the fact that yield, population size, and harvest-related costs contribute to the objective functional, which can be maximized by a rotational pattern in time. This suggests that a key ingredient for the cycles to occur is the possibility that both the harvest intensity and harvest timing can vary from season to season, that is, that both controls are time-varying.
The oscillations are remarkable because we do not know of another harvesting model that is able to destabilize Beverton-Holt population dynamics (except for a harvest control rule inducing a discontinuity in the dynamical system; Lois-Prados & Hilker, 2021). As cycles have important implications for population management (Barraquand et al., 2017), many empirical and theoretical studies have addressed the impact of harvesting on cycles (Anderson et al., 2008;Beddington & May, 1977;Dattani et al., 2011;Hilker & Liz, 2020;Hilker & Westerhoff, 2006;May et al., 1978;Milner-Gulland & Mace, 1998;Sah et al., 2013;Shelton & Mangel, 2011). Here, we have demonstrated for the first time that cycles can emerge as the optimal strategy when harvest timing is taken into account.
For the Ricker-type mechanistic model we observe a broader variety in results, including some interesting cycles. The prominence of cycles among these optimal control results likely is related to our choice of a Ricker-type mortality term. Since the Ricker model is known to exhibit cycles, we expected that cycles would be a possibility. Yet, we emphasize that the cycles reported occur for parameter values where the population dynamics in the absence of harvesting would be stable. Hence, the cycles are driven by the time-varying controls.
In many situations, control of harvest intensity seems to be much more important than control of harvest timing. Indeed, only with higher cost associated with harvest timing does γ J 0 ( , *) differ significantly from θ γ J ( *, *). When the time horizon is increased, several types of things may happen. The pattern of the optimal controls might not change at all (see Burton, 2020, pp. 107-116). The pattern might completely switch to a different one. The pattern may initially be the same, but then break down. A pattern may establish where before there was none. In general, over longer time periods the timing of harvest seems to be relatively more important, as seen in the increased relative difference between γ J 0 ( , *) and θ γ J ( *, *). Such optimal control results can give interesting dynamics, especially unusual oscillations in the population levels which would not be found in such models without time-varying timing and level of the harvest.
In this paper, we have considered seasonally reproducing species that have a recruitment pulse in a short period of the season and decline continuously throughout the remainder of the season. More continuous or repeated forms of recruitment during the season would require changes to the model structure. For instance, experimental removals of songbirds, grouse, waders, and raptors have been quickly replaced by 'floaters,' that is, migrating, nonbreeding adults without territory (Newton, 1992); this is likely to compensate density-dependent effects due to harvest timing. Also, strong Allee effects in a population have been shown in mathematical models to have reverse effects such that later harvest timing can be beneficial for populations (Cid et al., 2014). The current models assume recruitment directly into the harvestable population. In the future we plan to study the effect of harvest timing on agestructured populations. We also plan to develop more mechanistic models with different types of processes for reproduction and mortality and explore optimal control of both harvest timing and intensity for those models.