Spatially-explicit optimal reserve design models select best sites from a set of candidate sites to assemble nature reserves to protect species (or habitats) and these reserves display certain spatial attributes which are desirable for species. These models are formulated with linear 0–1 programming and solved using standard optimisation software, but they were run on different platforms, resulting in discrepant or even conflicting messages with regard to their computational efficiency. A fair and accurate comparison of the convenience of these models would be important for conservation planners who use these models. In this article, we considered eight models presented in literature and tested their computational efficiency using randomly generated data sets containing up to 2000 sites. We focused on reserve contiguity and compactness which are considered crucial to species persistence. Our results showed that two of these models, namely Williams (2002) and Önal et al. (2016), stand out as the most efficient models. We also found that the relative efficiency of these models depends on the scope of analysis. Specifically, the Williams (2002) model solves more of the test problems when contiguity is the only spatial attribute and a large subset of the candidate sites needs to be selected. When compactness is considered also, the Önal et al. (2016) model generally performs better. Large scale models are found to be difficult to solve in a reasonable period of time. We discussed factors that may affect those models’ computational efficiency, including model size, share of selected sites, model structure and input data. These results provide useful insight and guidance to conservation practitioners and researchers who focus on spatial aspects and work with large-scale data sets.

nature reserve design, mixed integer programming, computational efficiency, contiguity, compactness, spatial optimisation

Anthropogenic activities have caused tremendous impacts on ecosystems all over the world. These impacts include but are not limited to environmental pollution, habitat loss and fragmentation, invasion of exotic species and climate change (e.g. Aplet and McKinley 2017). These impacts led to an increased rate of extinction of species and impaired the services that ecosystems have been providing to humans. Establishing nature conservation reserves have been adopted across the globe as a direct way to restore ecosystems and protect species. However, resources devoted to this purpose such as land and budget have always been scarce and they often face conflicting demands from other sectors of society such as housing and industry. With the aim of designing ecologically effective and economically efficient nature reserves comes the science of reserve design (Kingsland 2002).

The reserve design problem is defined as selecting a subset from a larger set of candidate sites to assemble nature conservation reserves to provide adequate protection to targeted species or habitats. This problem has been studied using various approaches, including site scoring (e.g. Pressey and Nicholls 1989), gap analysis (e.g. Scott et al. 1993), heuristics (e.g. Pressey et al. 1997) and formal optimisation (e.g. Possingham et al. 2000, Williams et al. 2004). Gap analysis and heuristics have been widely used both in academic research and conservation practices, primarily because they are intuitive, easy to use and computationally convenient. Some popular nature reserve design software such as Marxan and C-plan use heuristics as their major solution algorithm (Ball et al. 2009, Sarkar 2012). However, it has long been demonstrated that these techniques find optimal solutions only occasionally and often generate solutions which significantly deviate from the optimal ones (Önal 2003, Williams et al. 2005). This means that they may not be allocating scarce conservation resources in the best possible manner. The problem can also be formulated as an optimisation model, specifically as a linear mixed integer programme (MIP), to determine the best selection of sites while meeting the conservation goals. The resulting model can be solved using MIP solvers such as CPLEX and GUROBI integrated in off-the-shelf mathematical modelling software such as GAMS and AMPL (AMPL 2016, GAMS 2016). Finding the best possible solution helps guarantee efficient allocation of conservation resources. Due to this advantage, optimal reserve design models have received substantial attention since the 1990s.

The early optimisation formulations of the reserve design problem adapted the set covering and maximal covering frameworks well known in the operations research literature (Cocks and Baird 1989, Underhill 1994, Camm et al. 1996, Church et al. 1996, Ando et al. 1998). These initial formulations neglected spatial attributes of the reserve, usually ending up with a solution where selected sites are highly scattered over the landscape and far apart from each other. A reserve with such a spatial configuration is not very meaningful in practice because both species survival and reserve management benefit from a spatially coherent reserve. The second generation optimum reserve design models incorporated spatial attributes to improve the coherence and functionality of the reserves (Williams et al. 2005). Amongst the spatial attributes discussed in literature, the most important ones are spatial contiguity (or connectivity) and compactness of the selected sites. Spatial contiguity means that any two sites in the reserve can be linked by a path of selected sites. Compactness is a more difficult concept to define, but, in general, it means that the selected sites must be close to each other (or grouped together) to the greatest possible extent. Some studies used the total distance between pairs of selected sites and the boundary length of the reserve as proxy measures for compactness. These two spatial attributes have been successfully modelled using graph theory, network flow theory and mathematical techniques (see Williams et al. 2005, Moilanen et al. 2009, Billionnet 2013 for reviews of those studies). Incorporating these spatial attributes into the optimisation frameworks requires more complex formulations, however, including a large number of additional variables and constraints. Thus, the reserve selection problem gets computationally more difficult when spatial configuration is of concern. For instance, a set covering formulation involving thousands of candidate sites could be solved in seconds (Rodrigues and Gaston 2002), but a compactness model involving about 700 sites could not be solved to optimality within 10 hours (Fischer and Church 2003).

Several modelling approaches have been introduced in the reserve design literature to incorporate spatial aspects in site selection. Due to their differences in structure and size, those models differ in their computational performance. Moreover, individual models were run on different platforms with different CPU speeds and solved using different MIP solvers and, therefore, it is not clear how those models compare to each other in terms of their computational efficiency. A fair and accurate comparison of the convenience of these models would be important for on-the-ground conservation practitioners. This is the main motivation for this article. For this purpose, we selected eight optimisation models incorporating (or could incorporate easily by proper modification) contiguity and compactness and conducted computational experiments with these models by running them on the same computer and using the same data sets that were generated randomly. We aimed to answer two questions: i) how do these models perform in solving contiguity and compactness problems with different sizes and ii) what are the likely factors affecting their computational efficiencies?

We carried out the tests on a Lenovo desktop computer with Intel Pentium (R) CPU of 2.8 gigahertz and RAM of 8.0 gigabyte and the operation system was Windows 64-bit. As the MIP solver, GUROBI 5.0 integrated in GAMS 23.9 was used.

Based on the computational experiments we conducted, we identified four important factors that affect the computational efficiency of MIP models in nature reserve design problems with contiguity and compactness considerations: 1) Model size, 2) The share of selected sites, 3) Model structure and 4) Input data. We discuss these factors below.

Model size is the most important factor affecting the computational efficiency of these models. As the model size (especially the number of binary variables) increases, solution time may increase substantially. Figure 1 depicts an example using the Williams (2002) model. As the number of candidate sites n increased from 100 to 1500, which gradually increased the number of variables from 6,091 to 18,831 and the number of constraints from 5,772 to 13,142, the solution time increased exponentially from 3 seconds to nearly 400 seconds. This is largely due to the number of nodes solved during the branch and bound procedure used by most MIP solvers. In general more nodes have to be solved when the model has a larger size.

Figure 1.

Relationship between processing time and the number of candidate sites (Williams (2002) model was used. population coefficient​α​=2.0. See text for detailed explanation of​α​).

Our test runs showed that the solution times also depend significantly on the share of selected sites in the optimal solution, which in turn depends on the specification of the population coefficient (​α​). This can be seen clearly, for instance, in the results obtained with the Shirabe (2005), Conrad et al. (2012) and Jafari and Hearne (2013) models for n=100 (see Table 2) and in the Williams (2002) model results for n=200 (see Table 3). Generally, these models follow a fast-slow-fast pattern, that is, as the share of selected sites increases, the solution time increases first, reaches a peak (where the share of selected sites is about 20–30 percent of all the candidate sites) and then decreases. This may be explained by the number of binary variables that have to be fixed at the values of 0 or 1 during the solution process. When the share of selected sites is small (or large), a large portion of binary variables will be fixed at 0 (or 1), leaving fewer combinations (a smaller branch and bound tree) that need to be searched for improving the objective function value (A brief description of the branch and bound algorithm can be found in Önal 2003). Similar results have been reported earlier by Williams (2002).

The Önal et al. (2016) model seems to be an exception. As noted above, it showed a relatively steady pattern in solution times for all population coefficient settings. The solution time tends to decrease only when a very small number of sites are to be selected.

Model structure is the third important factor that affects computational performance of the models we tested. Although the number of variables and constraints in a model may be substantial, due to the algebraic model structure, a large number of those variables and constraints may actually be redundant. Those redundant elements are eliminated during the preliminary solution procedure (Presolve in some GAMS solvers such as CPLEX and GUROBI) and the model sizes are much smaller after the Presolve. This is one of the main reasons for the superior computational advantage of the Önal et al. (2016) model. This situation was also observed in empirical applications of the Önal et al. (2016) model.

In addition, some algebraic forms may be integer friendly and some may not. For example, the graph-theoretic contiguity model introduced by Önal and Briers (2006) has the following constraint:

​​ (10)

where X_ij and U_j are binary variables as defined before and N_j is the set of j’s neighbour sites. This constraint states that if site j is selected, then at most four arcs can direct to it from its neighbour sites. The computational efficiency of the model can be improved substantially by modifying this constraint as in constraints (3) and (4). A comparison of the original Önal and Briers (2006) model and the modified model is depicted in Figure 2.

Figure 2.

Significance of model structure in affecting the model’s computational efficiency (number of sites n=36. The size of population to be covered equals to​α​times n.).

It is clear that replacing constraint (10) with constraints (3) and (4) slashes the solution time. We performed a paired samples t-test for each value of ​α​ to test the significance of the difference between these two models’ processing times. The result reported a p value less than 0.05 for each ​α​ value, indicating that the difference between the processing times is statistically significant. This result should be attributed to the fact that a constraint like (10) is integer unfriendly while the constraints (3) and (4) improve the model’s computational efficiency by helping construct a unimodular constraint matrix (ReVelle 1993, Williams 2002). Indeed, having a unimodular constraint matrix is one of the main reasons for the Williams (2002) model’s good efficiency.

Input data such as land price and species distribution across candidate sites are also crucial factors that can affect the model’s computational efficiency. Figure 3 depicts the processing times used by the models of Önal et al. (2016) and Williams (2002) when solving a 100-site problem with population coefficient ​α​=0.5. The problem size remained unchanged across the 50 runs while species population and site costs in each run were varied.

Figure 3.

Computational efficiency of the Önal et al. (2016) and Williams (2002) models in solving a 100-site problem with various species distributions and site costs data.

The shortest time used by the Önal et al. (2016) model is only about 1.2 seconds, while the longest time is 131.2 seconds. For the Williams (2002) model, the shortest and longest times are 0.3 seconds and 84.5 seconds, respectively. These results indicate that alterations in the input data may increase or decrease the solution time substantially even if the model size remains the same. This can also be explained by the workings of the branch and bound procedure. One data set may yield a ‘good’ integer solution in an early iteration and eliminate a large portion of the branch and bound tree; therefore, a fewer number of nodes needs to be solved (or just the opposite may occur). This seems to have happened in runs 4, 17, 33 and 44 with the Önal et al. (2016) model and in runs 1, 4, 7, 29, 39 and 44 with the Williams (2002) model.

We should note that none of these factors alone can determine the performance of a MIP model, instead they affect the model’s performance together. Solution techniques have been designed to speed up the computational process, such as supplying a cut value (Önal 2003, Cerdeira et al. 2005, Carvajal et al. 2013) or using a hybrid two-phase solution method (Conrad et al. 2012). We did not explore the effects of these techniques on the reserve design models tested here.

We have compared the computational efficiencies of eight representative optimal reserve design models that explicitly incorporate contiguity and compactness considerations. We identified the models presented by Williams (2002) and Önal et al. (2016) as the most efficient models. They both solved a 100-site contiguity problem easily, but encountered difficulty in solving the 200-site and larger contiguity problems. When compactness is also included, the Williams (2002) model had difficulty in solving a 100-site problem, while the Önal et al. (2016) model could solve a 200-site problem in a reasonable processing time (about 20 minutes). The Williams (2002) model could solve contiguity problems and some contiguity and compactness problems with up to 2000 sites when a large fraction of the candidate sites (such as 75 percent or more) needs to be selected. The Önal et al. (2016) model displayed a steady performance and outperformed the Williams (2002) model when both contiguity and compactness are considered, but still the model could not be solved within 1 hour when 500 sites are involved. We identified model size, share of selected sites, model structure and input data as the four important factors affecting the computational efficiency of spatially-explicit optimal reserve design models. When developing such a model, one should try to define as few variables and equations as possible and to construct the model in such a way that its constraint matrix is unimodular (ReVelle 1993). These findings may be useful to conservation planners who select and use models in reserve design or management practices. This work may also provide guidance to academic researchers who try to develop computationally efficient optimal reserve design models.

The current nature reserves, whether terrestrial or marine, may not be adequate in providing sufficient protection to biodiversity. For instance, only over 1.89% of the marine protected areas (MPAs) worldwide is covered by exclusively no-take MPAs, far from the Convention on Biological Diversity’s (CBD) Aichi Target 11 of 10% MPA coverage by 2020 and even further from the 30% target recommended by the IUCN World Parks Congress 2014 (IUCN 2018). It is almost certain that the current systems of nature reserves should be expanded. Acquiring more protected areas, however, would have to be done in a world where the human economy struggles to grow and the ecosystems and nonhuman species demands protection. Indeed, there is a conflict between economic growth and wildlife conservation (Czech 2014). Various disciplines such as ecological economics and environmental philosophy discussed this conflict and offered principles and guidelines to mediate the conflict and achieve sustainability. The science of nature reserve design or, more specifically, the optimal reserve design models, can contribute too. These models aim to solve practical conservation problems by incorporating ecological and socioeconomic considerations and help best prioritise conservation effort by guaranteeing optimal solutions. With the remarkable improvements in optimisation methods and computational power, it is reasonable to anticipate that spatial optimisation models will likely play a more significant role in the future than it does today in analysing and acquiring new protected areas.

The authors thank Dr. Falk Huettmann for his comments which greatly improve the quality of this manuscript. This work was supported by National Natural Science Foundation of China [71202096] and CREES Project [ILLU 05–0361].