## Abstract

Government agencies lack robust modelling tools to manage the spread of invasive alien species (IAS). In this paper, we combine optimal control and simulation methods with biological invasion spread theory to estimate the type of optimal policy and switching point of control efforts against a spreading IAS. We employ information-gap (info-gap) theory to assess how the optimal solutions differ from a policy that is most robustly immune to unacceptable outcomes. The model is applied to the potential invasion of the Colorado potato beetle in the UK. Under no uncertainty, we demonstrate that for many of the parameter combinations the optimal control policy corresponds to slowing down the invasion. The info-gap analysis shows that eradication policies identified as optimal under no uncertainty are robustly the best policies even under severe uncertainty, i.e. even if they are likely to turn into slowing down policies. We also show that the control of satellite colonies, if identified as optimal under no uncertainty, will also be a robust slowing down policy for IAS that can spread by long distance dispersal even for relatively ineffective control measures. The results suggest that agencies adopt management strategies that are robustly optimal, despite the severe uncertainties they face.

## 1. Introduction

The introduction of invasive alien species (IAS) is one of the main causes of the loss of global biodiversity. IAS can cause the extinction of vulnerable native species through predation, grazing, competition and habitat alteration (Mack 2000). In addition, IAS impose enormous costs on human health, agriculture, forestry, fisheries and water use, utilities, buildings and natural areas (US Congress Office of Technology Assessment 1993).

The IAS invasion process can be divided into three stages, namely entry, establishment and spread. Depending on which stage an invasion is at, different management decisions can be taken by government agencies responsible for managing IAS, i.e. decisions can be designed to prevent, eradicate, contain, slow down and/or accept the invasion (Sharov & Liebhold 1998; Leung *et al*. 2002). Identifying the optimal management policy decision and the point at which one should switch between different management decisions is a very complex task and the development of modelling tools to assist government agencies would be of great benefit.

Great insight has been gained into the bioeconomics of IAS management in recent years. Analytical models have been devoted to the optimal allocation of resources for preventative measures (Horan *et al*. 2002) or after establishment of an IAS in order to determine when eradication is the optimum policy (Eiswerth & Van Kooten 2002; Olson & Roy 2002). Other approaches that integrate the invasion stages have focused on assessing the optimal trade-off between exclusion and control efforts (Leung *et al*. 2002; Kim *et al*. 2006; Finnoff *et al*. 2007).

These modelling approaches have largely concentrated on IAS population dynamics instead of using theoretical spread models for IAS (a good review is provided by Hastings 1996). Instead, demographic models are in some cases employed as substitutes for spread models (e.g. logistic growth model). However, demographic models alone are unlikely to provide accurate predictions of invasion spread rates because, in order to relate population growth to spread velocity, it is necessary to take into account the spatial dispersal patterns of the invader (Higgins & Richardson 1996). Some notable exceptions of bioeconomic models that consider the dispersal patterns of the invader are those that incorporate the spread predictions of reaction–diffusion (R–D) models (Fisher 1937; Skellman 1951) and models that incorporate IAS performing long-distance dispersal, e.g. stratified diffusion models (Shigesada *et al*. 1995), into the management of insect and plant invasions (e.g. Sharov & Liebhold 1998; Sharov 2004; Cacho *et al*. 2008; Hyder *et al*. 2008).

In this paper, we combine optimal control methods (Pontryagin *et al*. 1962; Sethi & Thompson 2000) and simulation methods (Moody & Mack 1988; Taylor & Hastings 2004) with biological invasion spread theory (Fisher 1937; Skellman 1951; Shigesada *et al*. 1995) to estimate the optimal management policy against a spreading IAS. The effect of severe uncertainty on the decision making of the management of IAS spread is analysed using information-gap (info-gap) theory (Ben-Haim 2006). We specifically aim to answer the question: does the optimal IAS invasion management policy under no uncertainty differ from a policy that robustly protects us from unacceptable outcomes? The model is applied to the potential invasion by Colorado potato beetle (CB) *Leptinotarsa decemlineata* (Say) (Insecta: Coleoptera: Chrysomelidae) in the UK. CB is one of the most devastating pests of potato and since 1877 the UK has adopted a successful policy of preventing CB entry and eradication of any breeding colonies (Bartlett 1980; Waage *et al*. 2005; see electronic supplementary material). As a result, the study offers insight into the robustness of the optimal control policies and time at which to switch management efforts against an IAS invasion.

## 2. Methods

### 2.1. Basic model: spread by reaction–diffusion

The basic model is inspired by the works of Sharov & Liebhold (1998) and the integration of optimal control methods with epidemiological theory for disease spread management (Rowthorn *et al*. in press). We consider an already established IAS that is spreading following an R–D model (Fisher 1937; Skellman 1951). R–D models are partial differential equations where random diffusion in a homogeneous environment is assumed. The main parameters are *ϵ*, the intrinsic rate of population growth, and *d*, the diffusivity of the population. The solution of the R–D model is
2.1by which spread is predicted to follow a continuous expansion at an asymptotically constant radial velocity represented by *c*. A homogeneous landscape implies that *c* is constant in every direction, leading to a circular (or fraction of a circle if physical barriers occur) invasion front that is centred at the initial establishment point. The IAS invasion generates damages *D*(*x*) that are assumed to follow a linear relationship with the area invaded (Parker 1999)
2.2where *D** is the unit cost of damage caused by the IAS per unit of area invaded at the average population abundance and is assumed to be constant; *x* is the radius of the invaded area and in an optimal control theory context is called the ‘state variable’ because it describes the size (state) of the invasion; and *k* denotes the proportion of the circular invasion front that can spread without physical barriers.

The government agency can decide upon controlling the invasion using a moving barrier zone. A barrier zone is defined as the area bordering the expansion front of the invasion where management activities are carried out with the aim of reducing the velocity or even to lead to eradication of the invasion. For example, moving barrier zones were employed for the eradication of the boll weevil (*Anthonomous grandis*) in the USA (Sharov 2004). The barrier zone leads to costs *R*(*u,x*) which are proportional to the length of the invasion front (2π*x*/*k*) times the unit cost of control of an infested unit of area (*p*_{R}, which encompasses the unit cost of detection and control activities and is assumed constant) and the desired reduction in spread velocity (*u*):
2.3where *θ* is the rate of effectiveness of the control measures. In the context of optimal control theory, *u* is called the ‘control variable’ because it is the variable under the discretionary choice of the agency, i.e. the agency can choose its level to influence the state of the invasion.

#### 2.1.1. Optimal control under constrained control resources in a finite landscape

We assume that the government agency has limited resources and can at any one time accomplish only a reduction in the invasion spread velocity of up to *u*_{max}. We further assume that the total area susceptible to be invaded (susceptible range) can be described by a circle or fraction of a circle. The problem for the government agency is to choose the optimal time path for *u* (control variable) in order to minimize the net present value (NPV) of the total overall costs caused by the invasion and its management
2.4subject to
2.5
2.6
2.7
2.8where equation (2.4) is the objective function where *T* = time horizon, *r* = discount rate; equation (2.5) is the equation of motion of the radius of the invasion; equation (2.6) is the restriction of non-negativity of the control variable and the maximum value that *u* can take (*u*_{max}); equation (2.7) is the requirement of non-negativity of the state variable *x* and the constraint by which *x* cannot be bigger than the radius of the maximum susceptible range; and equation (2.8) is the initial boundary condition that reflects that at the moment of discovery the invasion has a size of *x* = *x*_{0} due to undetected spread. In addition, *c* = 0 for *x* = 0 (eradication) and *x* = *x*_{max} (total susceptible range is invaded).

We employ optimal control theory (Pontryagin maximum principle; Pontryagin *et al*. 1962) to solve the problem. The maximum principle aims at identifying the optimal path in time of the control variable in order to minimize the objective function. The constraints in the control variable and the spread dynamics of the IAS subject to the barrier zone are introduced in the optimization problem by means of the co-state variables. The mathematical vehicle by which the co-state variables take part in the optimal control problem is called the Hamiltonian function (Chiang 1992). The constraints in the state variable were introduced into the Hamiltonian using the indirect adjoining method (Sethi & Thompson 2000) and extensive numerical methods were employed to verify and obtain the solutions of the problem (see appendix A for details).

### 2.2. Extended model: spread by stratified diffusion

The basic model was extended to incorporate the possibility that the IAS can spread by long-distance dispersal events. A spatially implicit discrete-time spread model representing an IAS spreading by stratified diffusion (Shigesada *et al*. 1995) is employed. The model is composed of a system of difference equations (appendix A.5). The initial main colony grows following an R–D model (equation (2.1)) and generates new migrating individuals at a rate *λ* that is proportional to its area. The migrating individuals can establish and generate a new colony (‘satellite colony’). The model assumes complete spatial randomness on the location where the colonies establish (Bogich *et al*. 2008). The location of satellite colonies follows a scattered colony model by which the colonies are assumed to not coalesce with the other colonies in the time horizon considered (Shigesada *et al*. 1995).

The probability of establishment (*p*_{e}) is assumed equal to the density of the host in the landscape. We assume that there are no Allee effects (reduced survival of new colonies with a small number of individuals due, for example, to the difficulty of finding a mating partner) and the number of individuals capable of forming a colony arriving at the same location in the same time period is assumed equal to one. Once established, the satellite colony grows following also an R–D model and generates migrating individuals at a rate *λ*.

We assume that the government agency prioritizes the control of satellite colonies and uses the remainder of the allocated expenditure for the control of the initial main colony (identified as the most effective spread velocity reduction by Moody & Mack 1988). Instead of using a constraint on the maximum spread velocity reduction, we set a financial budget constraint for control measures. For simplicity, satellite colonies are assumed to be detected once they reach a threshold area (*A*_{d}).

We relax the assumption that the marginal cost of management is constant and equal to the unit cost of control of the IAS (*p*_{R}). We consider now the total management costs to be composed of searching costs (*σ*) and control costs (*p*_{R}*A*_{R}, where *A*_{R} is the area from which the IAS has been removed). Searching activities follow diminishing marginal returns that become relevant in the case of eradication campaigns that reduce the size of the invasion. This represents the increasing difficulty of finding the remaining IAS when the density of the population of the IAS is low. We use the following expression to represent the marginal cost of management (*c*′(*A*_{inv}), where *A*_{inv} is the area invaded by the IAS; modifying Burnett *et al*. 2007):
and
where *p*_{S} is the unit cost of searching a unit of area susceptible of invasion and *ω* is the radius of a security buffer zone around the initial size of the invasion (*x*_{0}) at discovery.

We incorporated the potential failure of the eradication campaign (by missing the control of some individuals) and the post-eradication costs, i.e. monitoring to certify that the IAS has been completely eradicated. Failure to completely eradicate leads to the re-emergence of colonies that would need to be eradicated in the next time periods. The probability of not achieving complete eradication was transformed into expected years of extended eradication campaigns. The control activities derived from the failure to eradicate are assumed to be applied to an area as big as the initial invaded area upon discovery plus a security buffer zone. The costs due to post-eradication activities consist of searching costs for 3 years starting from the year that the last colony re-emerged. We employed numerical methods to estimate the optimal time paths of control in the case of the extended model (see appendix A.4).

### 2.3. Model validation

The successful eradication campaign against CB in the UK in 1976 (see electronic supplementary material) was employed to validate whether the model would identify as optimal a policy that pursues eradication in these circumstances. Both basic and extended models identified a policy of eradication as the optimal policy. The main reason is that the area invaded was small at discovery. We considered the models valid conceptually since they are based on established biological invasion spread theory and optimization techniques.

### 2.4. Information-gap theory: robust decision making using severe uncertainty

Government agencies face severe uncertainty when making decisions regarding the control of IAS invasions. We adopt an info-gap approach (Ben-Haim 2006) to assess the robustness of the optimal control policy. Info-gap theory was developed by Ben-Haim (2006) to assist decision making in situations where the information about some areas of the system modelled is highly deficient. It seeks to obtain robust management decisions that fulfil a performance requirement for the widest uncertainty range in the model. Info-gap is especially suitable when probabilistic models of uncertainty are unreliable, inappropriate or unavailable (Regan *et al*. 2005). As an advantage over other probabilistic methods of uncertainty modelling, info-gap does not require one to make underlying assumptions about the distributions of uncertainty of the parameters. Info-gap is also concerned with the mapping of the propagation of uncertainty in the model into uncertainty in the decision-making domain (Ben-Haim 2004). We consider info-gap more adequate than traditional methods to deal with uncertainty for the case of management of the spread of IAS. The reasons are that we concern ourselves with providing government agencies with an analysis that does not require knowledge about the uncertainty distributions or ranges for the model parameters (this is very useful given the gaps in knowledge the agencies need to deal with) and our objective is to try to assess how different is the optimal policy obtained under no uncertainty from the most robust policy under severe uncertainty.

Three main components are required for info-gap analysis: (i) a mathematical process model, (ii) a performance requirement, and (iii) a model of uncertainty.

(i) The mathematical process model is an abstract representation of the reality. It condenses what the analyst considers to be the fundamental processes of the system. In our cases, the process models are the above-described basic and extended models.

(ii) The performance requirement of a decision (policy) is the level of an indicator of performance that is calculated by the process model. In our case, the performance indicator is the NPV of the total costs owing to the invasion and its control. The objective is to reduce the NPV of total costs as much as possible, but there is a level of the indicator that is considered unacceptably high, i.e. our minimum aspiration. We consider unacceptable to spend more resources for the control of the IAS than the costs of living with the IAS for the given time horizon under a policy of total acceptance of the invasion and assuming no uncertainty (baseline performance criterion). In addition, we consider one exigent and one modest performance criteria that are a decrease and increase of 50 per cent of the baseline performance criterion, respectively. We set the baseline performance criterion as
2.9by which the NPV of the total costs (TC) is less than or equal to the NPV of costs of living with the IAS when no control is used against the IAS invasion (TC_{Q}_{=0}).

(iii) The model for uncertainty expresses what is unknown about the parameters in the process model. It is an unbounded family of nested sets of possible values. Each set corresponds to a degree of knowledge deficiency according to the level of nesting (Ben-Haim 2004). We centre our analysis on the uncertainty of the IAS spread velocity (*c*), the effectiveness of the control measures (*θ*) and the probability of failure to achieve a complete eradication (*f*_{e}). The corresponding info-gap models are expressed as the sets: *U*_{c}(*α*, *c*^{be}), *U*_{θ}(*α*, *θ*^{be}) and *U*_{fe}(*α*, *f*_{e}^{be}). *α* is the info-gap between what is known and what needs to be known for an ideal solution and *c*^{be}, *θ*^{be} and *f*_{e}^{be} are our best estimates of the model parameters. The greater *α* is, the greater the range of possible variation. The value of *α* is unknown and unbounded and expresses the idea that possibilities expand as the info-gap grows, imbuing *α* with its meaning of ‘horizon of uncertainty’ (Ben-Haim 2004). We adopt an interval model of uncertainty where the fractional deviation of the values of the uncertain parameters *c*, *θ* and *f*_{e} from their best estimates is no greater than *α* (e.g. Regan *et al*. 2005; Nicholson & Possingham 2007; table 1):

Hence, given a horizon of uncertainty *α*, the value of the uncertain parameters are in the intervals

We only consider values that are lower than *θ*^{be} and higher than *f*_{e}^{be} and *c*^{be} to represent the most unfavourable conditions.

Info-gap theory identifies as the best policy the one that is most robustly satisfying (Ben-Haim, 2006), i.e. the goal is not to minimize the NPV of total costs but to maximize the reliability of an acceptable outcome. The most robust policy will be the one that presents most immunity to unacceptable outcomes. Once we have defined the process models, the performance criterion and the model of uncertainty, we can use them to estimate the robustness of the model. We employ a robustness function :
2.12where equation (2.12) states that is equal to the maximum value of α, in such a way that the maximum of the NPV of total costs fulfils the performance criterion (equation (2.9)) given *policy* option *i* and uncertainty in the parameters *c*, *θ* and *f*_{e} (that is modelled as a family of nested and unbounded sets that increases with *α*). We varied the parameters *c*, *θ* and *f*_{e} simultaneously at the same relative rate (e.g. Nicholson & Possingham 2007).

The types of policies compared using info-gap analysis were (i) the solution of the basic and extended models under no uncertainty, (ii) the agency perseveres on eradication during the time horizon, (iii) total acceptance of the invasion, and (iv) the agency controls for two years and ‘learns’ about the effectiveness of the control; if the invaded area is reduced (effective control), the agency will continue controlling for the invasion, and if the invaded area is not reduced, the agency accepts the invasion.

## 3. Results

### 3.1. Basic model: optimal policy and time to switch management policies

The type of optimal policy obtained by the basic model was the so-called ‘bang-bang’ optimal control by which it is optimal to apply maximum resources to control the invasion at the beginning and then switch to acceptance at some determined point in time. The optimal control policy was framed around the following points in time: the time to switch from management to acceptance that is the solution of the maximum principle without constraints in the state variable (*τ*) and the time when eradication is achieved (*t*_{erad}) or the total susceptible range is invaded (*t*_{x}_{max}) (see appendix A: equations (A 11) and (A 14) are solved to obtain *t*_{erad} or *t*_{x}_{max} and *τ*, respectively). These points in time depend on the characteristics of the government agency, IAS and initial conditions (see the electronic supplementary material for a sensitivity analysis of the switching time to model parameters and appendix A.2, where equation (A 14) relates model parameters to the switching time).

The types of optimal control policies are: (i) slow down the invasion and then switch to a policy of acceptance; (ii) we control the invasion during the entire time period and not the entire susceptible area is occupied; (iii) we accept the invasion without any attempt to control it; (iv) we control until an undetermined point in time between zero and *t*_{x}_{max} and we do not control when the entire susceptible range is occupied; (v) we control until total eradication (*t*_{erad}) and we stop controlling when eradication is achieved; and (vi) is analogous to case (i) except that the invasion invades the entire susceptible range during the acceptance policy. In cases (i), (ii) and (iii), the results were verified using extensive numerical methods for a variety of parameter values (see appendix A.4). In cases (iv), (v) and (vi), the state constraints and the control variable constraints are binding within the same time horizon. For these cases, no analytical solution could be found. Instead, numerical methods were used to obtain the switching point (see appendix A.4).

### 3.2. Basic and extended models: optimal type of policy under no uncertainty

The most common optimal management policy for both models and parameter ranges was one of slowing down the CB invasion (figure 1). The rationale behind a slowing down policy is the delay of the invasion so that no large areas can be occupied at the beginning of the time horizon. Because invasion size is directly related to impact costs, this delay implies the avoidance of considerable costs. Slowing down policies are especially attractive if there are cost-effective means to delay the invasion.

Both the basic and extended models identified eradication as the optimal policy for IAS invasions with the following characteristics: the invasions are discovered at a small size (figure 1*b*,*d*), the agency has high resources to control the invasion, the IAS spreads at low velocities, it is not costly to control and the IAS leads to high economic impacts on the invaded area (figure 1*a*,*c*). In the case of IAS that can spread by long-distance dispersal, the range of the parameters for which eradication was optimal was much smaller (figure 1, compare (*b*) with (*d*)). Total acceptance, in contrast, was optimal for fast spreading IAS that are very costly to control (figure 1*a*,*c*) and have already invaded a large area at the moment of discovery (figure 1*b*,*d*). Unexpectedly, total acceptance of IAS spreading by long-distance dispersal occurred for high budgets allocated for control and large areas invaded. The reason is that large budgets allow both for the eradication of the satellite colonies and the control of the main initial colony. Whereas the eradication of the satellite colonies is a cost-effective slowing down measure, the reduction of the main colony is not cost-effective when it is very large because no significant reductions in its spread velocity can be achieved.

### 3.3. Robustness of the model under severe uncertainty

We used the info-gap model (equations (2.10) and (2.11)) and the basic and extended models to evaluate the robustness function (equation (2.12)). The robustness curves for the four types of policies considered are demonstrated in figure 2. They show the NPV of total costs for all values of the horizon of uncertainty *α* between 0 and 1. For *α* = 0, i.e. there is no uncertainty regarding the parameters *c*, *θ* and *f*_{e}, the NPV of total costs is lowest for the solution provided by the basic and extended models (‘no uncertainty’ curves in figure 2).

In the case of spread by reaction–diffusion (figure 2*a*), the optimal policy under no uncertainty was one of eradication. As *α* increases the NPV of total costs increases faster for the solution provided by the basic model (no uncertainty, figure 2*a*) than for the ‘eradication policy’ (perseverance on eradication) and ‘learning agency’ (keep controlling as long as the invasion size is reduced) curves. Thus, the latter two policies represent the most robust policies for values of *α* below 0.2. The reason is that both policies persevere on eradication, whereas the no uncertainty curve does not attain a complete eradication under uncertainty and turns into a slowing down policy. When *α* reaches approximately 0.5 in the case of R–D spread, the learning agency realizes that control cannot attain a reduction in the invasion and accepts the invasion. This decision causes an increase of the NPV of total costs of the learning agency curve. The reason is that the avoided costs from slowing down the invasion were not seized by this type of agency. This leaves the eradication policy as the most robust policy. The eradication policy curve turns, for values of *α* > 0.5, into a slowing down policy and is still consistently the most robust curve for the three performance criteria considered. The robustness of the eradication policy is approximately 0.5, 0.6 and 0.7 for the exigent, baseline and modest performance criteria, respectively. The interpretation of a robustness of 0.7 is that all the parameters considered for the info-gap analysis can vary from their original values up to a fraction of 0.7 without leading to an NPV of total costs greater than the performance criteria. It is only for values of *α* greater than approximately 0.8 that the eradication policy curve presents higher NPV of total costs than the other policies.

For the case in which CB can perform long-distance dispersal, slowing down the invasion is the optimal policy if there is no uncertainty. Figure 2*b* shows that the optimal policy under no uncertainty is also the most robust policy under severe uncertainty for values of *α* very close to 1. The reason is that control measures directed to new satellite colonies considerably slow down the invasion. Even if such measures present very low effectiveness (*α* close to 1), the NPV of total costs would be much lower than a policy of total acceptance.

The consideration of different performance criteria (baseline, modest and exigent in figure 2) did not affect the selection of the most robust policy. However, if we had considered a very modest performance criterion, e.g. NPV of total costs of £12 million in the case of reaction–diffusion (figure 2*a*) and £68 million in the case of stratified diffusion (figure 2*b*), the policy of total acceptance would have been the most robust policy (with robustness equal to 1). We can see that there is a trade-off between how exigent is the aspiration level (low NPV of total costs) and how large is the immunity to desirable outcomes. Remarkably, as we increase our aspiration level, the types of policies identified as optimal by the models if there is no uncertainty are consistently the most robust policies under severe uncertainty.

## 4. Discussion

We have presented simple, yet general, optimal control and simulation models to identify the switching point for the management of a spreading invasion. The methods developed show how a policy of slowing down the invasion is in many cases the optimal approach. Furthermore, the results overturn the intuition that eradication campaigns that are optimal under no uncertainty are questionable if the probability of failure of the eradication campaign is high when considering uncertainty. They also contradict the intuition that agencies might not be able to best manage the invasion because of the severe uncertainty faced.

Under no uncertainty, we showed that eradication was only optimal for low initial sizes of invasion (Sharov 2004) and when the control measures taken were able to reduce the spreading velocity (figure 1*b*,*d*); this was particularly the case when the IAS was able to spread by long distance dispersal mechanisms (figure 1*d*). Eradication was also preferred for low velocity of spread of the invader (Cacho *et al*. 2008; figure 1*a*). However, the eminently most common optimal control policy for many of the parameter combinations was a policy where control switched to acceptance of the invasion within the time horizon, showing how, even if eradication was not feasible, slowing down the spread until a certain point in time was optimal (Sharov & Liebhold 1998).

Info-gap theory (Ben-Haim 2006) was employed to assess how different were the model solutions under no uncertainty from robust solutions under severe uncertainty, i.e. policies that would provide the greatest immunity to unacceptable outcomes. The types of policies identified as optimal under no uncertainty were, unexpectedly, also the most robust policies.

This result has implications for management: agencies should not be deterred from carrying out an eradication campaign even if there is severe uncertainty about the spread capacity of the invader, the effectiveness of the measures available and the probability of not achieving complete eradication. We show that even if the eradication campaign is unsuccessful, the control efforts will have served as very beneficial slowing down activities as long as the eradication campaign was initially identified as optimal under no uncertainty.

The economic reason behind the robustness of eradication policies (identified as optimal under no uncertainty) turning into slow down policies is that the agency, when eradicating, deals with an invasion that only represents a relatively small fraction of all the potential area invaded. The control activities are cost-effective for the control in small areas because they lead to high reductions of spread velocity, e.g. the same reduction of spread velocity in a small colony will imply a much smaller barrier zone in terms of area treated than in a large colony. It will be when the invasion becomes very extensive that control will not be cost-effective and there will be a switch to acceptance.

At times, however, optimal policies under no uncertainty will not be the most robust policies. This could happen when there is a high difference in benefits between an eradication policy and a slowing down policy. For instance, if the IAS was a quarantine pest, and should its presence cause the loss of important export markets, eradication would produce a sharp increase in total benefits owing to the re-establishment of export markets (Fraser *et al*. 2006). This nonlinearity might cause eradication to be the economically optimal policy under no uncertainty but not the most robust policy under severe uncertainty, as the efforts to pursue eradication would be of little benefit if the eradication was not finally achieved.

Identification of a slowing down policy (but not an eradication policy) as optimal under no uncertainty does not imply that starting an eradication campaign that might turn into a slowing down policy is a robust strategy. This could be due to high differences in the costs of eradication and slowing down campaigns. For instance, costly detection campaigns being carried out in inaccessible areas or high post-eradication costs. For the case of CB spreading by R–D, however, the consideration of these costs did not indicate a policy of eradication to be not optimal or to be less robust than a policy of total acceptance. This remained the case even when the CB invasion could never be totally eradicated (figure 2*a*, ‘never eradicate’ scenario).

Agencies that can learn from the effectiveness of the control campaigns in the initial years and give up on control when the invasion cannot be reduced were shown to perform less robustly than agencies that are consistent on eradication and slowing down campaigns. In reality, agencies will not give up on control if they perceive that a cost-effective slowing down of the invasion is being achieved. The normal sequence of management performed by an agency is, if considered optimal, to attempt to accomplish eradication and if it is not possible, slow down the invasion until the area invaded means slowing down activities are not cost-effective. The results of the models suggest that the normal sequence of management of the agencies is actually very close to the most robustly optimal approach available despite the severe gaps in the information they have.

In the case of CB potential invasion in the UK, the current policy of eradication was identified as the optimal and most robust policy by the model provided that discovery occurs for small invaded areas. If an invasion breached eradication campaigns and reached areas that made eradication unfeasible, a policy of slowing down the CB invasion by control of new satellite colonies would be the most robust policy.

Further improvements of the model could be brought about by relaxing the assumption of a homogeneous landscape. This could be achieved by adopting a spatially explicit simulation approach (Gilligan *et al*. 2007). In this case, more flexible spread models such as metapopulation models (Rowthorn *et al*. in press), cellular automata or individual-based models could be considered. In addition, agencies deal with multiple IAS simultaneously, little is known about the optimal management and the effect of prioritizing strategies, sometimes in response to international obligations.

In this paper, a combination of optimal control and simulation methods, biological invasion spread theory and info-gap theory was presented to estimate the optimal policy and switching point of invasion management campaigns. The models can be easily applied to other IAS (e.g. information to parametrize the model for several IAS can be found in Andow *et al*. (1990) and Waage *et al*. (2005)), although comprehensive information for many IAS might be lacking. This reflects the importance of compiling and sharing data on historical eradication campaigns and spread events at the international level. The models represent useful tools for preliminary exploration of the optimal policy against a spreading IAS given a set of biological and economic parameters. The integration of the dispersal patterns of the invader in the bioeconomic modelling of IAS invasions is strongly recommended. This integration together with info-gap analysis will help us to construct robust tools for the management of IAS invasions.

## Acknowledgements

We thank Christos Gavriel and two anonymous referees for their helpful comments that greatly improved this paper. We thank Emily Nicholson for introducing us to info-gap theory. The research was funded by a grant of the UK DEFRA (NB 53 9005) and the Rural Economy and Land Use Program (RES-229-25-0005). The views expressed herein are exclusively our own and do not represent the views of DEFRA.

## Appendix A

### A.1. Application of Pontryagin maximum principle

Taking into account the constraints, the resulting current value (Chiang 1992) Lagrangian–Hamiltonian equation is (the minimization problem was transformed into a maximization problem by multiplying the objective function by minus one)

A 1Applying the Pontryagin maximum principle (Pontryagin *et al*. 1962) and the indirect adjoining method for pure inequality state constraints (Sethi & Thompson 2000), the following set of conditions can be obtained:

Equation (A 2) indicates that the optimal control *u**(*t*) must maximize the Lagrangian–Hamiltonian for all *t* within the time horizon considered; equation (A 3) is the equation of motion for *x*; equation (A 4) is the equation of motion of the co-state variable λ_{c} modified for the current value Hamiltonian; equation (A 5) is the transversality conditions for a vertical terminal line at *t* = *T*; equations (A 6) and (A 7) are the conditions due to the constrained state variable (equation (2.7)). The complementary-slackness conditions state that η_{1} and η_{2}, the Lagrangian multipliers, will be zero unless *x* = 0 and *x*=*x*_{max}, respectively (the state constraints become binding).

Given that *L*_{c} is linear in the control variable *u*, we obtain a bang-bang solution for *u* (Chiang 1992). ∂*L*_{c}/∂*u* is called the switching function and is referred to as *σ*. To maximize *L*_{c}, the boundary solution *u** = 0 (acceptance of invasion) should be chosen if *σ* is negative and *u** = *u*_{max} will be chosen if *σ* is positive. Only if *σ* = 0 for a positive interval of time does the Lagrangian–Hamiltonian not depend on *u* and we obtain a singular sub-arc. The optimal control is described as
A 8where
A 9If there is a singular sub-arc (*σ* = 0) *u** (0<*u**<*u*_{max}), equation (A 9) indicates that the marginal avoided cost of reducing the size of the invasion (λ_{c}) must equal the marginal costs that led to such reduction. If there is no singular sub-arc, the optimal control contains only the extreme levels of control and there will be as many switches (from *u** = *u*_{max} to *u** = 0 or vice versa) as the number of roots that *σ* has. Because it was not possible to check for singular sub-arcs analytically, we employed numerical methods instead (see appendix A.4).

### A.2. Unconstrained solution

We initially attempt to solve the maximum principle as an unconstrained problem assuming that the constraints in the state variable are not binding for all *t* (equations (A 6) and (A 7)). We need to determine the roots in the switching function *σ*. In this case, the sign of *σ* depends on the co-state variable *λ*_{c}. We proceed to investigate the form of *λ*_{c}. Equation (A 4) can be expressed as

We proposed as optimal path of *u* one in which *u* equals *u*_{max} for all *t* and evaluated the state and a co-state solution.

Since *u* is constant and equal to *u*_{max}, we can integrate equation (A 3) and apply the boundary condition (2.8) to obtain

Substituting equation (A 11) into equation (A 10) and setting *u*=*u*_{max}, we can solve equation (A 10) as an ordinary differential equation
A 12where *a* is an integration constant that is defined by applying the boundary conditions (A 5). Rearranging terms, we obtain

Then, substituting equations (A 11) and (A 13) into equation (A 9) and rearranging terms, the expression of *σ* results:
A 14where
The switching function *σ* has only one root for *t* = *τ*. We employed numerical methods to obtain *τ* (we used the function *FindRoot* in Mathematica, Wolfram Research Inc. 2005). The solution obtained is to apply maximum control efforts until *t* = *τ* and then accept the invasion from *τ* to *T*. Extensive numerical analysis for a variety of parameter values confirmed that this solution was a maximum (see appendix A.4).

We define the junction points in time when the state variable become binding as: *t*_{erad}, time of eradication, and *t*_{x}_{max}, time of total invasion. There are three cases when the pure state variable inequality constraints are not binding for all *t* ∈ [0,*T*] (eradication or total invasion does not occur, i.e. *η*_{1} = 0 and *η*_{2} = 0) and the solution of the unconstrained problem is an optimum
A 15

These are cases (i), (ii) and (iii) of §3.1.

### A.3. Constrained solution

In the cases when the junction points *t*_{erad} or *t*_{x}_{max} ∈ [0,*T*], the state constraints become binding. Both the control variable and the state variable constraints are binding in the same time horizon and it was not possible to solve the co-state variables analytically. Hence, the switching point between control and no control before reaching the junction point (state variable constraint becomes binding) could not be obtained analytically. Once the junction point has been reached, by complementary slackness we know that (equation (A 6) or (A 7)): (*c* − *θu*) = 0; since, by definition, *c* = 0 when *x* = 0 or *x*=*x*_{max}, *u* has to be zero whilst the state constraint is binding. We used numerical methods to obtain the switching point (see appendix A.4).

### A.4. Numerical methods

A discrete time-step simulation model that represented the problem was developed. The control levels were discretized so that the model could be solved using a genetic algorithm (Goldberg 1989) with @RiskOptimizer (Palisade-Corporation 2006). A genetic algorithm is a numerical optimization method inspired from evolutionary biology. A computer simulation is performed where a population of abstract representations of candidate solutions of the optimization problem (chromosomes) evolves to better solutions according to a fitness criterion (Goldberg 1989). Our fitness criterion was that chromosomes leading to lower mean of the NPV of total costs were the fittest ones. We used an initial population of optimal paths for *u* of 500, a crossover rate of 0.5 and a mutation rate of 0.1 (Palisade-Corporation 2006). The model was used to verify the conclusions of the analytical approach and both models presented quantitative agreement for the parameter sets that led to cases (i), (ii) and (iii). This procedure was used to derive the switching point for cases (iv), (v) and (vi) (see §3) and to obtain the optimal path in the extended model. The genetic algorithm did not identify singular sub-arcs in any of the solutions for the parameters considered.

### A.5. Extended model: spread by stratified diffusion

The model is represented by the system of difference equations
where *A*_{0}, *A*_{1}, … , *A*_{T} represents the area invaded when *t* = 0, 1, … *T*; and *N*_{t} is the total number of colonies at time period *t*. Each difference equation expresses the increases in size of the existing colonies following R–D and the generation of new satellite colonies minus the area treated against the IAS.

## Footnotes

- Received July 2, 2009.
- Accepted August 7, 2009.

- © 2009 The Royal Society