Solving conservation planning problems with integer linear programming (appendices) HAWTHORNE L. BEYER1,3 , YANN DUJARDIN2 , MATT WATTS1 & HUGH P. POSSINGHAM1 1
ARC Centre of Excellence for Environmental Decisions, Centre for Biodiversity & Conservation Science, University of Queensland, Brisbane QLD 4072, Australia 2 CSIRO Ecosystem Sciences, Ecosciences Precinct, Dutton Park, QLD 4102, Australia 3 Corresponding author email:
[email protected]
A
Summary
Section B provides a detailed case study of how to implement the Marxan with Zones (Watts et al., 2009) objective functions in an integer linear programming (ILP) framework using the linearisation techniques described in the main text. Section C is a detailed description of implementing ILP optimisation problems using R and Gurobi, including documented code for the simulations reported in the main text.
B
Linearisation of the Marxan With Zones reserve selection optimisation problem
A variation on the Marxan problem in which each planning unit can take on more than two states is called ‘Marxan with Zones’ (Watts et al., 2009). This functionality allows several additional types of problems to be addressed depending on what the stratification of states represent (management actions, zoning types, etc). For example, Klein et al. (2010) use this approach to assign marine planning units to one of five zones that regulated the conservation status and level of fishing allowed, ranging from complete conservation protection to unrestricted commercial fishery. Marxan With Zones solves the following optimisation problem:
1
min
N X Z X i=1 k=1
s.t.
Z X
cik xik +
Z X Z X k=1 m=1
bkm
N X N X
xik xjm vijkm
i=1 j=1
xik = 1, i ∈ N
k=1 N X Z X
(1) aij pjk xij ≥ Tj , j ∈ N
i=1 k=1 N X
aij xik ≥ Ujk , j ∈ N, k ∈ K
i=1
xi ∈ {0, 1}, i ∈ N where N and Z refer to the number of planning units and zones (actions) respectively, xik is a binary decision variable that is 1 when unit i is assigned to zone k and 0 otherwise. In the first part of the objective function cik is the cost of planning unit i when assigned to zone k. Thus, the costs associated with selecting a planning unit can vary among the different zone types and planning units. The second part of the objective function provides a way of enforcing or facilitating the assignment of zones between pairs of units (often neighbouring units). For example, conflicting zones in neighbouring units could be discouraged by assigning a higher penalty in the array v. It is also possible to specify v to facilitate greater clumping of units of a given zone type by assigning a cost to same-type pairs that is less than the cost of any different-type pairs. Although v is large (a four dimensional array), it may contain many zeros if costs are only assigned to neighbouring units, making it more straightforward to linearise this problem (because we can ignore any instances of xik xjm for which vijkm = 0). As with the Marxan objective function (Eqn 4), parameter matrix bkm can be adjusted by the decision maker to control the overall strength of the penalties. The first constraint ensures that each planning unit is assigned to exactly one zone. In the second constraint aij represents the amount or value of feature j in unit i and pjk is the proportion of that feature that is preserved when the unit is assigned to zone k. Thus, summed over all planning units and zones, this constraint ensures that minimum representation targets (Tj ) are met. The final constraint ensures that targets (U) for the amounts of feature j are represented by specific zone types. For example, if the targets refer to species abundances, and the zones refer to multiple land use types, then this last constraint could be used to ensure that at least 50% of the abundance was represented by fully protected conservation areas. As with Marxan, the constraints in the Marxan with Zones problem are incorporated into the objective function in the form of a shortfall penalty function (Watts et al., 2009) that allows solutions to be found even when all targets are not met. When all targets are met, however, the shortfall penalty is 0. For the purposes of implementing Marxan with Zones as an ILP problem we choose not to merge the constraints with the objective function as this allows for greater flexibility in customising the problem. The Marxan with Zones problem can be linearised in a similar way as the linearisation of the Marxan problem. The first component of the objective function is linear with respect to the decision variables x and requires no modification. The second component of the objective function contains the quadratic expression xik xjm , which can be linearised by replacing it with a 2
new decision variable zijkm and implementing the following additional constraints:
zijkm − xik − xjm ≤ −1 This constraint ensures that zijkm = 1 when xik = xjm = 1, and that zij = 0 when either xik = 0 or xjm = 0. Note that this constraint differs from the two constraints needed to linearise the Marxan objective function because the sign of the quadratic term differs. We can also express the objective function using set notation, which makes it clearer that the penalty term is only evaluated for neighbours:
min
N X Z X i=1 k=1
cik xik +
Z X Z X k=1 m=1
bkm
X
xik xjm vijkm
(2)
(i,j)∈E
where E defines the set of all neighbouring planning units. In the worst case scenario the linearisation of the Marxan with Zones problem could result in approximately N Z + N 2 Z 2 decision variables and N 2 Z 2 constraints in addition to the structural constraints defining the targets. However, the linearisation terms can be omitted whenever vijkm = 0, which typically applies to all non-connected planning units. In fact, in most applications the array v is sparse. For example, in the case of a hexagon grid a planning unit will have at most 6 neighbours, thereby adding less than 6N Z decision variables (note that planning units on the edges of the grid have fewer than 6 neighbours) and 6N Z constraints, which is a fraction of the size of the worst case scenario. Nevertheless, the dimension of the problem can still increase rapidly with N and Z. A number of customisations of this problem are possible. It may be appropriate to allow no zone assignment to some units in cases where only a subset of unit and zone combinations PZ are required to meet targets. This can be achieved by changing the first constraint to k=1 xik ≤ 1, thereby relaxing the constraint that forces all planning units to be assigned to exactly one zone. Another problem is forcing some units to be assigned particular zones. Even though there is no ‘decision’ to be made for such units it is often necessary to include them in optimisation problems because they can affect zone assignments to other units via the second component of the optimisation function. For example, in a planning problem it may be unrealistic to transform urban areas into any other land use type, so they would be fixed as urban. But if a penalty has been implemented for placing protected areas adjacent to urban area these planning units would need to be included in the problem. Although it is possible to force zone assignments by assigning large costs to alternative zones in the cost matrix cik , a more direct and definitive approach is to implement a constraint xik = 1 for each ik corresponding to existing urban planning units. We note that such constraints can also be implemented in Marxan with Zones (see Watts et al., 2008). Similarly, to prevent the assignment of incompatible zones in adjacent planning units (e.g. a dredge spoil dump next to a marine protected area) two constraints can be implemented:
3
xik + xjm ≤ 1, i ∈ N, j ∈ N xim + xjk ≤ 1, i ∈ N, j ∈ N where k and m are the incompatible zones. In Marxan With Zones these relationships will be enforced less explicitly by means of the values that are assigned to the cost matrix cik and the penalty array vijkm . Indeed, adding additional constraints to enforce these effects is detrimental to processing times as it increases the dimension size of the problem. The benefit of using the constraint approach, however, is that it is definitive. Finding suitable values for c and v to prevent undesirable solutions from arising often involves a process of trial and error, particularly when parameter b is also being adjusted as this affects the penalties. The decision maker must inspect the solutions to determine whether undesirable assignments have occurred and this can be difficult to do effectively for very large problems. Implemented as a constraint, we are certain that these undesirable assignments are precluded, thereby removing the reliance on users to appropriately tune the problem formulation.
C
Implementation of linear programming methods to conservation planning
Here, we explain how to solve ILP problems using R and the commercial optimisation software Gurobi. This includes code for simulating 10,000 planning units with cost and species data for 10 species (Section 2), using this data to solve a simple, linear ILP reserve selection problem (Section 3), using this data to solve the non-linear Marxan objective function (Section 4), and some general advice for solving large ILP problems (Section 5). We recommend that readers understand how the simpler reserve selection problem is solved (Section 3) before moving on to the Marxan problem (Section 4), which is considerably more complex. We note that copying and pasting code from a PDF format can sometimes cause issues with how quotation marks are translated into plain text. If you experience problems you may wish to first look at whether replacing the single and double quotation marks resolves the problem.
C.1
Alternative linear programming software
Alternatives to Gurobi include the commercial software CPLEX (IBM, Inc., 2009), MatLab (MATLAB, 2014) and Mathematica (Mathematica, 2014), and the open source software lpSolve (Berkelaar et al., 2005). We have found that processing times with lpSolve can be considerably longer than Gurobi. Furthermore, the R interface to lpSolve (the ‘lpSolve’ R package) lacks some important features, including the ability to accept constraint matrices in the sparse format and setting to determine when the solver exits (a time limit or gap statistic). Thus, only the smaller ILP problems we evaluate could be implemented though the lpSolve R interface and the more complex problems would have to be implemented through the lpSolve programming interface (the API) or possibly the command line interface. We do not explore this functionality here. 4
C.2
Generating spatial datasets
The following R code was used to simulate input data for the case study, including 10 raster (matrix) datasets representing the value of each planning unit (one raster cell) to each of 10 species of conservation concern, and another raster representing the cost of selecting a planning unit for inclusion in the reserve system. To run this code it is necessary to first install the RandomFields library in R.
# generates simulated spatial (raster) datasets # if not installed, the RandomFields library can be installed using this command: # install.packages(’RandomFields’) library(RandomFields) # 1. Species data representing the value of each planning unit # (raster cell) to each species # user defined parameters: # set the number of rows and columns of the raster: nr
