- Email: [email protected]

Contents lists available at SciVerse ScienceDirect

Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa

A dynamic model of controlling invasive species İ. Esra Büyüktahtakın a,∗ , Zhuo Feng b , George Frisvold c , Ferenc Szidarovszky b , Aaryn Olsson d a

Department of Industrial and Manufacturing Engineering, Wichita State University, Wichita, KS 67260-0035, United States

b

Department of Systems and Industrial Engineering, University of Arizona, Tucson, AZ 85721-0020, United States

c

Department of Agricultural and Resource Economics, University of Arizona, Tucson, AZ 85721-0023, United States

d

School of Earth Sciences and Environmental Sustainability, Northern Arizona University, Flagstaff, AZ 86011-5694, United States

article

info

Article history: Received 16 September 2010 Received in revised form 29 June 2011 Accepted 11 August 2011 Keywords: Invasive species Biological invasion Non-native species Optimal control Integer programming

abstract A dynamic model of controlling invasive weeds is first developed which is a large scale, nonlinear 0-1 integer programming problem. This model is then applied for the case of control of the invasive grass, Pennisetum ciliare (buffelgrass), in the Arizona desert. The large size of the problem makes the application of direct optimization methods impossible, instead the most frequently suggested strategies were analyzed and their consequences compared. The model is more advanced and complex than those examined in earlier studies. Published by Elsevier Ltd

1. Introduction This study examines the spread and control of invasive weeds as a spatial-dynamic problem [1,2]. The spatial-dynamic framework is then applied to the management of Pennisetum ciliare (buffelgrass), an invasive fire-promoting African bunchgrass that spreads very rapidly across the Sonoran Desert. Buffelgrass forms dense stands, crowding out native species, thereby reducing species directly. Furthermore, it increases wildfire risk by promoting a grass/fire cycle in an ecosystem with little or no known natural fire history [3]. Details of the biological aspects of the buffelgrass invasion are given in our earlier paper [4]. Studies of invasive species control are often based on population growth without considering the important aspects of spatial heterogeneity. For example, damage often depends on the size, but not on the spatial distribution of the invasion, while costs are invariant to location [5,4,6]. The factors that can vary across the landscape are the potential for an invasive weed to become established, its carrying capacity, the cost of control and the damages it causes. Typically, damage induced by invasive weeds depends on a local population impact (e.g., through population density, level of competition, flammable biomass production, etc.) as well as the presence of valued or threatened resources at the site of infestation or nearby. The land manager’s problem is to minimize the total damage during a given time interval subject to budget and labor constraints. In our model, we consider multiple sources of spatial heterogeneity in the damage and cost functions. We utilize a gridded landscape in which each cell in the grid represents one acre of land. Damage caused by the buffelgrass invasion is defined as a weighted sum of damages to different resources, where the weights reflect management priorities. Opportunities to treat buffelgrass during the year are limited and, as such, management decisions prescribe to treat or not to treat a cell over a given period of time. Population in untreated cells grows according to the natural growth model; whereas, population decreases

∗

Corresponding author. E-mail address: [email protected] (İ.E Büyüktahtakın).

0898-1221/$ – see front matter. Published by Elsevier Ltd doi:10.1016/j.camwa.2011.08.037

İ.E Büyüktahtakın et al. / Computers and Mathematics with Applications 62 (2011) 3326–3333

3327

in treated cells based on the kill rate of the treatment. The objective function is the total discounted damage during the entire considered time interval. Most of the earlier works performed comparisons of specific strategies. Moody and Mack [7] and Martin et al. [8] compare the efficiency of different treatment variants on small, invasive weed populations over large, already established populations. Wadsworth et al. [9] compare alternative strategies based on proximity for human settlements, weed population size, age and spatial distribution similar to Stevens and Falk [3]. The benefits and costs of biological control program are compared in [10]. This paper makes the following contributions to the invasive species control literature. First, we develop a new dynamic programming model. No such optimization model was formulated earlier for buffelgrass treatment. In addition, our cost and damage functions are more realistic than those in earlier studies. Second, we assess the consequences of treatment on the growth function, which is also a new contribution. Moreover, we determine minimum labor requirements for creating defensible space for each valuable resource, which could improve labor allocation across areas. Because of the very large size of the problem, true optimization cannot be performed. Instead we follow the tradition of research comparing specific strategies; however, we consider additional factors such as reflecting future possible damage to the current time period which were not considered in earlier works. This paper develops as follows. In Section 2, the dynamic programming model will be introduced, and as a descriptive model, in Section 3 it will be used to assert and compare the consequences of the most commonly suggested treatment strategies. Section 4 will conclude the paper. 2. Mathematical model Let T years be the length of the considered time interval and t ∈ {0, 1, . . . , T } be any year. The area is a rectangular array of cells, each cell is identified by its coordinates (i, j), where i ∈ {1, 2, . . . , I }, and j ∈ {1, 2, . . . , J }. The management should decide which cells will be treated in each year. Therefore, the decision variables are the binary variables xi,j (t ) such that xi,j (t ) =

1, 0,

if cell (i, j) is treated in year t otherwise.

(1)

The objective function is the total damage caused by the weed population in the entire area during the entire time interval. The damage Dij (t ) in cell (i, j) depends on the resource at risk in cell (i, j), its proportion of the buffelgrass population Nij (t ) within that cell, and the management priority of the corresponding resource. So, the objective function has the form minimize

J − I − T −

Dij (t )

(2)

i=1 j=1 t =0

where Dij (t ) =

K −

dk Rijk Nij (t )

(3)

k =0

for damage weights dk affecting k ∈ {1, 2, . . . , K } resources (Rijk ) proportional to the population Nij (t ) in a cell (i, j) with Rijk ∈ {0, 1}. Rijk represents the spatially explicit quantification of a resource either as a binary (e.g., presence/absence of critical infrastructure) or continuous variable (e.g., biological diversity or probability density of an endangered species). Note that damage depends on the type of resources the manager wants to protect and, as such, different management priorities applied to the same landscape can lead to different estimates of damage (e.g., a resource manager tasked with conserving biodiversity would perceive the damage induced by the presence of the weed population differently than a manager tasked with preserving critical infrastructure). The state variables are defined as the actual population volume Nij (t ) of each cell at each year. The natural change in population density in two consecutive years can be formulated as Nij (t ) − Nij (t − 1) = g Nij (t − 1) , N ∼i∼j (t − 1) , Kij

(4)

where Kij is the carrying capacity in cell (i, j) and N ∼i∼j (t − 1) is the vector of the buffelgrass population densities in the eight cells surrounding cell (i, j) in year t − 1. Kij is the spatial representation of habitat suitability as predicted by a bioclimatic envelope (BEM) [11]. BEMs map the ecological niche in parameter space onto locations in physical space and provide an overall probability of each cell on the landscape being suitable for a species. Because many different methods and spatial covariates are involved in the development of BEMs, a discrete equation is not available. Rather, we refer to the habitat suitability as a static spatial layer. Kij is the suitability scaled by the cell area. The function g () has a logistic growth form, where population grows at an increasing rate at first, then at a decreasing rate as the population approaches the carrying capacity in the cell. Each cell receives propagules from plants within the cell and from the neighboring cells. The rate at

3328

İ.E Büyüktahtakın et al. / Computers and Mathematics with Applications 62 (2011) 3326–3333

which a cell receives propagules from its neighbors is described by an exponential decay function. In describing the effect of treatment on cell (i, j) we have two possibilities. Before treatment, the population is obtained from Eq. (3) as NijB (t ) = Nij (t − 1) + g Nij (t − 1) , N ∼i∼j (t − 1) , Kij .

If this value is less than the critical population

NijC ,

(5)

the treatment eradicates the buffelgrass from the cell. Otherwise, after

the treatment the population becomes NijB (t ) (1 − k), where k is the kill rate of herbicide treatment. In our study k = 0.9 is selected. So the state transition relations have the from NijB (t ) (1 − xij (t )),

Ni,j (t ) =

NijB

(t ) (1 − kxij (t )),

if NijB (t ) < NijC otherwise.

(6)

The first case shows that without treatment (xij (t ) = 0), Ni,j (t ) = NijB (t ) and if treatment is applied (xij (t ) = 1), then the cell population is eradicated. The second case can be interpreted similarly. If treatment is applied (xij (t ) = 1), then the population decreases according to the kill rate, otherwise it does not change. The land management faces budget and labor constraints. The budget constraint can be written as J I − −

Cij (t ) xij (t ) ≤ B (t )

(7)

i=1 j=1

where B(t ) is the annual budget in year t. The cost of treating cell (i, j) is assumed to increase linearly in the pre-treatment population NijB (t ), average cell slope sij and the distance dij of the cell from the closest road: Cij (t ) = c1 + c2 NijB (t ) + c3 sij + c4 dij

(8)

where the cost parameters c1 , c2 , c3 and c4 were estimated as (0.31, 2.91, 0.19, and 0.04), respectively using least squares regression based on recent records of actual treatments obtained from the Southern Arizona Buffelgrass Working Group [12]. These data are summarized from costs incurred by the University of Arizona, Saguaro National Park, and Arizona Department of Transportation. The labor constraint has a similar form, J I − −

Lij (t ) xij (t ) ≤ L (t )

(9)

i=1 j=1

where L(t ) is the maximum available labor in year t, and Lij (t ) is the labor requirement for treating cell (i, j) in year t. It is assumed that this is also a linear function of the pre-treatment population, average cell slope and the distance from the closest road: Lij (t ) = l1 + l2 NijB (t ) + l3 sij + l4 dij

(10)

where the values of parameters l1 , l2 , l3 and l4 are also estimated from actual treatment data obtained from the Southern Arizona Buffelgrass Working Group [12] and represent the needed treatment time (time/area treated) that increases with increasing slope and distance to road. The derived cost and damage functions were discussed with the members of the group, who approved them. A sensitivity analysis on the parameter values will be the subject of our future work. This model is a typical dynamic programming problem, so theoretically it is solvable. However in our case I = 40, J = 50, so we have 2000 cells. The length of the considered time interval is at least 10 years. So the number of state variables is 20,000 or more. We have 2000 decision variables in each year, so the number of binary decision variables is at least 20,000 implying that direct optimization of the nonlinear model is not practical. Instead, we will use the dynamic model as a descriptive tool to assess the consequences of a number of specific treatment strategies. 3. Analysis of specific strategies In this section, suggested treatment strategies will be considered and examined. 3.1. Preventing species establishment In this case, we study how much labor is required to prevent the establishment of buffelgrass after it first appears, and what the consequences of delays in initiating a treatment regime are. It is assumed that buffelgrass has initially discovered in 2.4% of the cells. We consider labor requirements given different start years until buffelgrass becomes completely eliminated. The strategy is to treat each infected cell. The results are shown in Fig. 1. There are neither labor nor budget constraints. It is clear that delay in the starting year significantly increases the labor requirement at the beginning, and then it decreases significantly in each consecutive year. The discounted total labor cost for each case is shown in Fig. 2. The horizontal line represents the starting year. The discount rate is chosen as 2% and 4% for the sake of comparison. Clearly, the total cost is a fast growing function of the starting year.

İ.E Büyüktahtakın et al. / Computers and Mathematics with Applications 62 (2011) 3326–3333

3329

30 27 400-hour work teams

24 21

Year 1

18

Year 3

15 Year 5

12 9

Year 9

6

Year 13

3 0 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 Year Fig. 1. Annual labor requirements with (400 h. work teams) to prevent buffelgrass establishment, with different starting years of control.

$700

($ thousands)

$600

r = 2%

$603

r = 4%

$487

$500

$422 $367

$400 $298

$300

$243

$200 $100 $-

$71 $64

$78 $70

Year 1

Year 3

$119 $105

Year 5 Year 9 Starting Year

Year 13

Year 17

Fig. 2. Discounted cumulative labor costs to prevent buffelgrass establishment ($ 1000) under 2% and 4% discount rate.

3.2. Decisions with binding resource constraints (A) Consider first a static rule by minimizing the total damage at a single year subject to only labor constraint. The objective function for time period t is to minimize the total damage: minimize

J I − −

Dij (t )

i=1 j=1

subject to the constraint J I − −

Lij (t ) xij (t ) ≤ L (t ) .

i =1 j =1

In Eq. (3), the data for the previous time period t − 1 are obtained from the optimal solution for that period. This approach starts at t = 1, when the initial data are given, then based on the left over resources we optimize for t = 2, and so on until the last period is considered. If all cells generating positive damages are treated and labor is still remaining, then we treat cells with the remaining labor to minimize the total buffelgrass population. The first step is to minimize function (2) with T = 1. In the second step, we minimize J I − −

Nij (1)

(11)

i =1 j =1

subject to the dynamic relations (5) and (6) with t = 1 with given Nij (0) value and the labor constraint (9), where L(1) is the remaining labor. (B) In the first step of the previous heuristic strategy, we can consider damages to the native vegetation, to saguaros and to buildings. The actual damage value Dij (1) depends on the priority of the management. In the second step, the objective

3330

İ.E Büyüktahtakın et al. / Computers and Mathematics with Applications 62 (2011) 3326–3333

1200

Damage Index

1000 800

(A) (B) w=0.5 (C) (D) w=0.5 No Treatment

600 400 200 0

1 4 7 10 13 16 19 22 25 28

Year Fig. 3a. Vegetation, L = 400, start year = 9.

1200

Damage Index

1000 (A) (B) w=0.5 (C) (D) w=0.5 No Treatment

800 600 400 200 0 1 3 5 7 9 1113151719212325 27 29

Year Fig. 3b. Vegetation, L = 2000, start year = 9.

function is a weighted average of the actual and potential future damages: J I − − (w Dij (1) + (1 − w)Dmax ) ij

(12)

i=1 j=1

where Dmax is the maximal potential damage buffelgrass can cause in cell (i, j), i.e. the damage of buffelgrass if no treatment ij is applied. (C) The Buffelgrass Working Group [12] recommends treating areas in at least three consecutive years. This is the strategy actually used in Arizona for treatment. If a cell with population near its carrying capacity is treated only once, then buffelgrass’ logistic growth will push the population back to the fast part of its growth path. This is the reason why repeated treatment is suggested. After the third treatment, we follow the strategy introduced in part (A). (D) Same as above with the only difference that after the third treatment, the strategy given in part (B) is applied. The static optimization problems with objective functions (2) with T = 1, (11) and (12) are standard 0–1 integer programming problems, which were solved by using the branch and bound method [13] and software ILOG CPLEX [14]. The above-described heuristic rules were applied for the three different resources at risk separately. In all cases, the results were obtained within a couple of seconds. Figs. 3a–3b show the damages for the above described strategies during a 20-year period if the treatment starts at year 9 and labor is constrained at 400 and 2000 h, when damages to only native vegetation are considered. In the initial years of treatment, all four strategies achieve comparable damage reduction. Prioritizing cells for treatment based on potential future damage as well as current damage (strategy B) outperforms control based solely on current damage (strategy A). This is the case with both low (Fig. 3a) and higher (Fig. 3b) labor availability. As might be expected, the improvement in damage reduction from this ‘‘potential damage weighting’’ trends upwards over time. While the difference in damage reduction does not increase monotonically in every year, it does trend upwards. Moreover strategy B dominates strategy A in every year except the first year of treatment. Again, this is to be expected given that strategy A minimizes current damage. The ‘‘treat three times’’ rule with potential damage weighting (strategy D) also outperforms the simple ‘‘treat three times’’ rule (strategy C) in every year except the initial treatment year. Again the improvement of D over C trends upwards over time, although it does not increase monotonically. Strategy D outperforms the other three in all years except the first year of treatment. Likewise, Strategy B ranks second in each year but year 9 (the first year of treatment). Treating three times is clearly an improvement when potential damage weighting is used (D dominates B). The improvement from treating three times over static optimization is more ambiguous, however. While strategy C dominates A with low labor availability (Fig. 3a), damage reductions are quite comparable until later years with larger labor endowments (Fig. 3b). Treating three times (Strategy C) emphasizes driving population (and

İ.E Büyüktahtakın et al. / Computers and Mathematics with Applications 62 (2011) 3326–3333

3331

1200

Damage Index

1000 800

(A) (B) w=0.5

600

(C) (D) w=0.5

400

No Treatment 200 0 1 3 5 7 9 11 13 15 17 19 21 23 25 27 29

Year Fig. 4a. Saguaro, L = 400, start year = 9.

1200

Damage Index

1000 800 (A) 600

(B) w=0.5 (C)

400

(D) w=0.5 No Treatment

200 0 1 3 5 7 9 11 13 15 17 19 21 23 25 27 29

Year Fig. 4b. Saguaro, L = 2000, start year = 9.

damage) down to zero in cells with established weed populations. This emphasis, however, is at the expense to treating new infestations, which may spread to cells with high-value resources to protect. Whether to focus on established populations or new outlying infestations is a long-debated issue in the invasive species management literature. Moody and Mack [7] is a classic reference in this literature. Two factors can combine account for the non-monotonicity in differences in damage. First, the logistic population growth curve implies that populations in a given cell (and damaged caused) may change little from one year to the next (if the cell has either very low or very dense populations) or change dramatically, (if the population is near the inflection point). Second, treatment reduces cell population by 90%. Again, populations can fall substantially if pre-treatment populations are large. Treating a cell a second or third time may contribute little to reducing current damage, especially compared to treating other cells with new populations, but driving an infestation from a particular area can have longer-term benefits. The superiority of treating three times over current damage minimization (C vs. A) appears sensitive to labor availability, Figs. 4a–4b show similar results if damages to saguaros are the main concern of the management. Again heuristic strategy (D) always outperformed the others, with strategy B as the second best. Treating three times improves damage reduction over current damage minimization with more severe labor constraints (Fig. 4a) and if potential damage weighting is used (D vs. C in both Figs. 4a and 4b). With large labor supplies (Fig. 4b), however, strategies A and C are comparable. Fig. 5 shows how strategy selection affects the risk on buildings. Since structures are usually on the periphery of the grid and initial infestation are relatively far from buildings, the damage to buildings is easy to reduce. With a labor constraint of 800 h per year, each strategy resulted in significant damage reduction. Again, strategy (D) was the best. The ranking of the other three strategies were in flux in early years, but in the longer term strategy (D) dominates strategy (B), which dominates strategy (C) dominating strategy (A) by year 20. We emphasize that a robust result was that strategy (D) that includes both Buffelgrass Working Group recommendations consistently dominated all other strategies. This was true for each resource requiring protection (vegetation, saguaros, buildings) and for each level of labor endowment. We also examined how the labor constraint affected the damage trajectories (Figs. 6a–6c). In each case, strategy (D) is followed and figures include the trajectories without any treatment in order to show the comparison of the optimal solutions to this worst case scenario. Fig. 6a shows damages to native vegetation under different labor constraints. In each case damage trends down, then reverses and trends up. This eventual upward trend occurs even with large labor supplies (L = 2000). Riparian vegetation is relatively dispersed throughout the grid. This makes it more difficult to create ‘‘defensible space’’ around a small number of cells. Fig. 6b shows that for saguaros, different labor

3332

İ.E Büyüktahtakın et al. / Computers and Mathematics with Applications 62 (2011) 3326–3333

120

Damage Index

100 80 (A) (B) w=0.5 (C) (D) w=0.5

60 40 20 0 1 3 5 7 9 11 13 15 17 19 21 23 25 27 29

Year Fig. 5. Buildings, L = 800, start year = 9.

1200

Damage Index

1000 800

L=0 L=400 L=800 L=1200 L=1600 L=2000

600 400 200 0

1 3 5 7 9 11 13 15 17 19 21 23 25 27 29

Year Fig. 6a. Riparian damage function with different annual labor budgets (strategy D).

1200

Damage Index

1000 800

L=0 L=400 L=800 L=1200 L=1600 L=2000

600 400 200 0 1 3 5 7 9 11 13 15 17 19 21 23 25 27 29

Year Fig. 6b. Sagurao damage function with different annual labor budgets (strategy D).

endowments can stabilize damage at different levels. For example, when L = 400, damage remains near 600 from years 23–29; when L = 1200, damage remains below 300 from years 23–29; when L = 1600, damage remains below 200 for years 23–29. Saguaros are less dispersed throughout the grid than riparian vegetation, but more dispersed than buildings. For labor endowments L ≥ 800, damage can be driven close to zero. There are three things to note. First, damage is not actually driven to zero. Labor is sufficient to create a ‘‘defensible space’’ around buildings. Populations of buffelgrass can continue to grow just beyond the fringe where they would create damage. This means that it is a continuing invasion threat even if it is not creating much current damage. Second, this implies that policies of creating defensible space need to be ‘‘long-term’’ policies with annual commitments of labor. It appears that defensible space can be created with L = 800 about as well as if L = 2000. Rather than applying L = 2000 to buildings in this area, our results suggest that many labor hours could be diverted to protecting other resources without increasing risk to buildings. Third, this implies that our approach can be used to determine minimum labor requirements for creating defensible space. This could improve labor allocation across areas. 4. Conclusions A dynamic model was developed to find optimal long-term treatment of buffelgrass in the Arizona desert. The model is a 0–1 nonlinear integer programming problem; however, the nonlinearity of the model and the huge number of decision

İ.E Büyüktahtakın et al. / Computers and Mathematics with Applications 62 (2011) 3326–3333

3333

1200

Damage Index

1000 L=0 L=400 L=800 L=1200 L=1600 L=2000

800 600 400 200 0 1 3 5 7 9 11 13 15 17 19 21 23 25 27 29

Year Fig. 6c. Building damage function with different annual labor budgets (strategy D).

variables make direct optimization impractical. Instead, the model was used as a descriptive tool to assess the consequences of given, most frequently suggested treatment strategies. In each case, only static optimization problems were solved. Four particular strategies were considered and compared. In all cases, heuristic strategy (D) outperformed the others. Heuristic strategy (D) mostly resembles the preliminary management recommendations of an inter-agency Buffelgrass Working Group [12]. In our future study, we will consider other strategy variants including multi-objective and game theoretical approaches and examine their effectiveness. Future research also includes the uncertain budget and invasion growth size for the control of the invasive species. References [1] M. Smith, J. Sanchirico, J.E. Wilen, The economics of spatial-dynamic processes: applications to renewable resources, Journal of Environmental Economics and Management 57 (2009) 104–121. [2] J.E. Wilen, Economics of spatial-dynamic processes, American Journal of Agricultural Economics 89 (2007) 1134–1144. [3] J. Stevens, D.A. Falk, Can buffelgrass invasions be controlled in the American southwest? Using invasion ecology theory to understand buffelgrass success and develop comprehensive restoration and management, Ecological Restoration 27 (2009) 417–427. [4] İ.E. Büyüktahtakın, F. Zhuo, A. Olsson, G. Frisvold, F. Szidarovszky, Invasive species control by a dynamic spatial process: an application to buffelgrass (Pennisetum ciliare) in Arizona, Journal of Environmental Management (2010) Under revision. [5] K. Burnett, B. Kaiser, J. Roumasset, Economic lessons from control efforts for an invasive species: Miconia calvescens in Hawaii, Journal of Forest Economics 13 (2007) 151–167. [6] D.I.S. Odom, O.J. Cacho, J.A. Sinden, G.R. Griffith, Policies for the management of weeds in natural ecosystems: the case of scotch broom (Cytisus scoparius, L.) in an Australian national park, Ecological Economics 44 (2003) 119–135. [7] M.E. Moody, R.N. Mack, Controlling the spread of plant invasions: the importance of nascent foci, Journal of Applied Ecology 25 (1988) 1009–1021. [8] B. Martin, D. Hanna, N. Korb, L. Frid, Decision Analysis of Alternative Invasive Weed Management Strategies for Three Montana Landscapes, Prepared by The Nature Conservancy of Montana, Helena MT and ESSA Technologies Ltd., Vancouver, BC, 2007, 34 pp. [9] R.A. Wadsworth, Y.C. Collingham, S.G. Willis, B. Huntley, P.E. Hulme, Simulating the spread and management of alien riparian weeds: are they out of control? Journal of Applied Ecology 37 (Suppl. 1) (2000) 28–38. [10] K.M. Jetter, J.M. DiTomaso, D.J. Drake, K.M. Klonsky, M.J. Pitcairn, D.A. Sumner, Biological control of yellow starthistle, in: D.A. Sumner (Ed.), Exotic Pests and Diseases: Biology and Economics for Biosecurity, Iowa State University Press, Ames, Iowa, 2003, pp. 225–241. [11] A.D. Olsson, Modeling the potential distribution of buffelgrass in southern Arizona: a comparison of logistic regression and artificial neural networks, Master’s Thesis, Prepared for the School of Renewable Natural Resources, University of Arizona, 2006. [12] A. Rogstad (Ed.), Southern Arizona Buffelgrass Strategic Plan: A Regional Guide for Control, Mitigation, and Restoration, USA Buffelgrass Working Group, Tucson, AZ, 2008. [13] G.L. Nemhauser, L.A. Wolsey, Integer and Combinatorial Optimization, John Wiley & Sons, New York, NY, 1988. [14] ILOG CPLEX, IBM ILOG CPLEX: High-performance mathematical programming engine, 2010. http://www-01.ibm.com/software/integration/ optimization/cplex/.