Systematic Conservation Planning for Sustainable Land-use Policies: A Constrained Partitioning Approach to Reserve Selection and Design.

1 Introduction

In the context of the global biodiversity crisis, it is urgent to strengthen the conservation of natural habitats through the establishment of nature reserves [Prendergast et al., 1993; Haddad et al., 2015]. Accordingly, two of the United Nations Sustainable goals have been focused on the conservation of marine and terrestrial habitats. However, to be efficient the selection and the design of reserves must be based on a systematic approach [Margules and Pressey, 2000], and sustainable land use policies must promote socioeconomic development and nature conservation in a balanced way. On top of that, the design of protected buffer zones surrounding sensitive areas is an important aspect, as it can mitigate negative edge-effects and contribute to reducing fragmentation, by fostering recolonization and habitat restoration [Harris, 1988; Fahrig, 2003]. Buffer zones are, for instance, an important element of UNESCO World Heritage Sites [Feilden and Jokilehto, 1998] and Man and the Biosphere reserves [Batisse, 1982]. More recently, the IUCN protected areas management categories provided a more comprehensive framework by promoting the partitioning of the space into several levels of protection, which can be nested [Dudley, 2008].

∗Contact Author

In conservation biology, this concern lies in the framework of systematic conservation planning and formalizes as the reserve selection and design problem, which also lies in the more recent framework of computational sustainability [Gomes, 2009]. The reserve selection and design problem aims at partitioning the geographical space into at least two regions: one dedicated to habitats and biodiversity conservation, the other for socioeconomic development. However, effective strategies usually involve more regions with several nested levels of protection. Each region is defined by a combination of coverage and spatial constraints, and some other constraints such as the buffer zone, that can involve several regions. Finally, optimization objectives can be defined, such as minimizing the cost of a region or maximizing the coverage of certain features. Figure 1 depicts an example with a grid partitioned into three regions.

Many reserve selection and design models have been devised, e.g. [ReVelle et al., 2002; Williams et al., 2005; Sarkar, 2012; Dilkina et al., 2017]. Usually, these methods rely on ad-hoc heuristics, metaheuristics, or mixedinteger linear programming (MILP). More recently, CP models for and reserve selection and design has been published [Bessière et al., 2015; Justeau-Allaire et al., 2018]. While the first is focused on wildlife corridor design, the second is, to the best of our knowledge, the only generic CP model for this problem which combines both covering and spatial aspects. Sadly, in this model the search is focused on a single region, the reserve. As a consequence, socioeconomic constraints cannot be expressed on the out-reserve area. In addition, no more than two regions can be defined, and the buffer zone constraint is lacking. As a matter of fact, this constraint has, to our knowledge, only been modeled by some MILP approaches [Williams et al., 2005; Billionnet, 2013] in a local fashion for three-regions configurations (core area, buffer zone and out-reserve area). This approach has some limitations since it does not account for the reciprocity between the regions (e.g. a buffer zone can exist without a core area).

The current models are limited in their flexibility, as each address a specific subset of variants of the general problem. Because these subsets are usually different, it is also difficult to provide systematic model comparisons. For instance, to the best of our knowledge no existing MILP model is able to tackle problems with more than three regions. Some heuristics and metaheuristics are able to, e.g. Marxan with zones [Watts et al., 2009], but they don't provide strict and explicit control over spatial attributes. Finally, no existing model implements a complete buffer zone constraint. However, it is of great interest for managers to have the possibility of seamlessly considering those aspects.

Although rarely used in this context, CP is a good candidate for devising a more generic model for reserve selection and design. In fact, CP provides both flexibility and expressiveness, as it allows seamless integration of complex and heterogeneous constraints into a single model. Moreover, CP is a complete approach that is able to provide satisfiability and optimality proofs. In this paper, we show how to encode a CP model that allows the definition of an arbitrary number of regions, on which any constraint can be seamlessly applied. On top of that, we provide a complete and generic formulation of the buffer zone constraint, which can be reused to compose more complex spatial constraints. Finally, we experiment our model on a real-world dataset addressing the problem of forest fragmentation in the south of New Caledonia, a biodiversity hotspot located in the South Pacific (this dataset was already used in [Justeau-Allaire et al., 2018]). Through this use case, we show that the genericity provided by our model allows addressing problems that were not possible to tackle until now.

2 ProblemDescription

2.1 The Grid

The problem applies in a discretized geographical space. Several types of tessellation are possible (e.g. square grid, hexagonal grid, irregular grid). In accordance to available data, we only consider the regular square grid. However, the methods described in this paper can easily be transposed to other types of tessellation. In the context of reserve selection and design, a grid cell is called a site (the terms parcel and planning unit are also used in the literature).

Definition 1 (2D regular square grid). A M × N regular square grid S is a tessellation of the 2-dimensional space into |S| = M×N unit squares (called sites in our context), M and N are respectively the number of rows and columns. A site is uniquely identified by its zero-based matrix coordinates xi (row) and xj (column) or by its flattened index x = Nxi +xj (this indexing is independent of the tessellation).

Let S be a regular square grid.

Definition 2 (ω-connected neighborhood). We denote by ωconnected neighborhood any function Γω : S ↦→P(S) (with P the power set) that associates to x ∈ S a set Γω(x) ⊆ S representing the neighbors of x, according to the label ω.

Given that, we derive some specific neighborhoods that will be useful in the rest of the paper (some of them are illustrated in Figure 2), ∀x ∈ S, with x ≡ (xi, xj):

Definition 3 (Connected component). A set cc ⊆ S of connected sites (according to a neighborhood definition) is called a connected component (abbreviated CC). Note: This is an extension of the classical graph definition of connected component. If the neighborhood is not symmetric (e.g. representing environmental flows), the definition of strongly connected component must be used.

Definition 4 (Region). A set R ⊆ S associated to a single land-use policy is called a region. A non-empty region is composed of one or several CCs. The set of CCs associated to a region R is denoted by cc(R).

2.2 The Features

In the context of reserve selection and design, a feature corresponds to a characteristic of the geographical space that can be spatially represented with a numerical value for each site. A feature can represent biodiversity (e.g. species, habitats), or socioeconomic values (e.g. exploitable land, customary area). Three data types can describe a feature: binary data (e.g. presence of exploitable land), quantitative data (e.g. species abundance) and probabilistic data (e.g. species distribution model, SDM, representing either a probability of presence or a habitat suitability index). We denote the value associated to a feature f in the site x by vf (x).

2.3 Towards a Partitioning Formulation

In [Justeau-Allaire et al., 2018], the problem is stated as follows: given a M × N grid S, find R ⊆ S such that a set of constraints C are satisfied by R. In this formulation, the constraints are organized into two categories, coverage and spatial constraint. They can be formalized as follows.

Coverage Constraints

Let R be a region and F a set of features:

Constraint A (Covered features). R is a cover of F if every feature of F is present in at least one site of R. In this context, a feature f is considered to be covered by a site x if and only if vf (x) ≥ 1 that is ∀f ∈ F, ∃x ∈ R, vf (x) ≥ 1.

Constraint B (α-covered features). The constraint holds if and only if every feature of F has a probability of at least α to lie in R: ∀f ∈ F, ∏x∈S (1 - vf (x)) ≤ 1 - α.

Constraint C (k-redundant features). The constraint holds if and only if every feature of F is present in at least k site of R: ∀f ∈ F, ∃X ⊆ R, |X | ≥ k ∧ ∀x ∈ X , vf (x) ≥ 1.

Spatial Constraints

Let R be a region:

Constraint D (Number of CCs, aka Number of reserves). The constraint holds if and only if the number of CCs of R is bounded: minNbCC ≤ |cc(R)| ≤ maxNbCC .

Constraint E (Region size, aka Reserve System Area). The constraint holds if and only if the size of the region is bounded: minSize ≤ |R| ≤ maxSize.

Constraint F (CCs size, aka Reserve areas). The constraint holds if and only if the smallest (respectively largest) CC of R contains at least minSizeCC (respectively maxSizeCC) sites: ∀C ∈ cc(R), minSizeCC ≤ |C| ≤ maxSizeCC .

We suggest here a more generic formulation of the problem: given a M ×N grid S, find a partitioning of S into n regions {R0, ...,Rn-1} such that each region Ru satisfies a given set of constraints Cu ⊆ {A, . . . , F}. Using this formulation, any constraint in the previous catalog can be seamlessly applied to any region.

3 AGeneric CPModel

In this section we introduce a generic CP model associated with the partitioning formulation of the reserve selection and design problem.

3.1 The Base Model

Decision Variables

To each site x ∈ S we associate an integer variable ρx ∈ [0, n[. If x lies in Ru then ρx = u. An instantiation of these variables defines de facto a partitioning of the grid into (at most) n regions: ∀x ∈ S, ρx ∈ [0, n[.

Set Variables

Set variables are an abstraction providing an efficient, expressive and compact way to solve combinatorial problems through set-based modeling. The domain of a set variable X is a set interval [X,X], with X and X two sets (respectively the lower and upper bounds). Given that, an instantiation of X is a subset of X, such that X is a subset of X [Gervet, 1995]: X ∈ [X, X] ⇔ X ⊆ X ⊆ X.

Each region is represented by a set variable Ru ∈ [∅, P(S)] that is channeled with the decision variables such that ρx = u if and only if x ∈ Ru. This channeling ensures that the sets are all disjoint and that they form a partition of S: ∀u ∈ [0, n[, Ru ∈ [∅, S], ρx = u ⇔ x ∈ Ru.

Graph Variables

Similarly to set variables, graph variables are an abstraction providing an efficient and expressive way to model combinatorial problems with graphs. A graph variable G is defined by a graph interval [G, G] (with G and G two graphs, respectively the lower and upper bounds), such that an instantiation of G is a subgraph of G and G is a subgraph of G [Dooms, 2006; Fages, 2014]: G ∈ [G, G] ⇔ G ⊆ G ⊆ G.

First, we define the full spatial graph GS = (S, ES), where a vertex is associated with each site of the grid S, and such that there is an edge between two vertices if and only if they are adjacent in the grid, ES = {(x, y) | y ∈ Γ4(x)}. In the scope of this paper we represent adjacency with the fourconnected neighborhood Γ4, but any other definition could be used. An illustration of GS is provided in Figure 3. Then, similarly to the model defined in [Justeau-Allaire et al., 2018], to each region Ru we associate a graph variable Gu = (Ru,Eu). Gu is the subgraph of GS induced by Ru. These variables will be used to define connectivity and size constraints on the regions and their CCs. Each graph Gu has the empty graph as lower bound, and the full spatial graph GS as upper bound. Formally, for all u ∈ [0, n[, Gu = (Ru,Eu) ⊆ GS, such that Eu = {(x, y) | (x, y) ∈ Ru 2 ∧ y ∈ Γ4(x)}. An illustration of {G0,G1,G2} for a 3-regions partitioning is provided in Figure 3.

User Constraints

Finally, any constraint defined in 2.3 can be seamlessly applied to any region. For more details on how to apply them, refer to [Justeau-Allaire et al., 2018]. On top of that, the genericity provided by the partitioning perspective as well as the expressiveness provided by set and graph variables allows the modeling of more complex constraints. We detail, in the following, the buffer zone constraint.

3.2 The Buffer Zone Constraint

As mentioned in the introduction, the buffer zone constraint is of great interest for managers. Here, we provide a set-based generic formulation of the constraint. A buffer zone is an area separating the periphery of two areas. We first introduce the notion of generalized neighborhood.

Definition 5 (Generalized neighborhood). Let R be a region, and Γω a neighborhood definition. The generalized neighborhood Γω(R) of R is the union of the neighborhood of every site in R: Γω(R) = ∪x∈R Γω(x).

Constraint G (Buffer zone constraint). Let Γω be a neighborhood definition, Ru and Rv be two regions, and B a third region intended to be a buffer zone between Ru and Rv. The buffer zone constraint buffer[Γω](Ru,Rv,B) holds if and only if:

Consistency of the Buffer Zone Constraint

We now study the consistency of the buffer zone constraint. To this end, we rely on the definitions and results introduced by [Walsh, 2003] on set/multiset constraints. In particular, we rely on the definition of bound consistency (BC) and on the following one: a constraint decomposition is a normal form if and only if decomposing constraints are at most ternary and of the form X ⊆ Y,X = Y ∪ Z,X = Y ∩ Z,X = Y - Z, X ̸= Y, |X| = I, occ(I, X) = m or occ(m, X) = I, where X, Y and Z are set/multiset variables, I an integer variable and m an integer. The constraint occ(m, X) = I (respectively occ(I, X) = m) holds if and only if I (respectively m) equals the number of occurrences of m (respectively I) in X.

Proposition 1. The constraint N = Γω(R) (with N and R set variables) can be decomposed into a normal form.

Proof. The decomposition is based on additional variables N R x , for each site x:

(2)

These additional variables N R x are constrained by the following decomposition:

Then, ∪x∈S N R x can be decomposed into a normal formwith |S| -1 ternary union constraints and |S| -2 intermediary set variables: M1 = N R 0 ∪N R 1 ∧M2 =M1 ∪N R 2 ∧ ... ∧N = M|S|-2 ∪ N R |S|-1 .

Corollary 1. BC on N = Γω(R) is equivalent to BC on the decomposed normal form.

Proof. It is straightforward to decompose N = Γω(R) into a normal form by Proposition 1. Moreover, N = Γω(R) does not contain a repeated occurrence of variables, thus, BC on N = Γω(R) is equivalent to BC on the decomposed normal form (see [Walsh, 2003]).

Corollary 2. The constraint buffer[Γω](Ru,Rv,B) can be decomposed into a normal form. Consequently, BC on buffer[Γω](Ru,Rv,B) is equivalent to BC on the decomposed normal form.

Proof. The buffer constraint as defined in (1) can be decomposed as following:

Then, Corollary 1's exact same reasoning applies.

Time Complexity of the Buffer Zone Constraint

Finally, we study the time complexity of the buffer zone constraint filtering. Once more, we rely on [Walsh, 2003] which provides both filtering rules to enforce BC on the decomposed normal form of the constraint, and the worst-case time complexity associated with such a filtering: it is at most O(enm2 ) where e is the number of constraints, n the number of variables and m the maximum cardinality of the set variables.

Proposition 2. Enforcing Bound Consistency on buffer [Γω](Ru, Rv, B) can be done in O(|S|4 ).

Proof. According to the proofs of Proposition 1 and Corollary 2, the number of constraints and the number of additional variables in the normal form are in O(|S|). Since the cardinality of all the set variables is bounded by |S|, the complexity is in O(|S|4 ).

3.3 Extending the Model

We now consider a set of operational scenarios to demonstrate how expressive our CP model is, and how it can be extended and adapted to complex requirements.

Width of the Buffer Zone

The buffer constraint allows great control over the spatial attributes of the buffer zone through the neighborhood definition. A good example consists in controlling the width of the buffer zone. To do so, we introduce an alternative version of any neighborhood which integrates the notion of width, as illustrated in Figure 4.

Definition 6 (k-wide neighborhood). Let S be a regular square grid and Γω a neighborhood definition. The k-wide neighborhood of a site x ∈ S, denoted by Γkω (x), is defined by the following recursion:

Nesting Several Protection Levels

Another interesting application of the buffer zone constraint is the spatial nesting of regions representing several levels of protection. Such a configuration, thereby, can be desired to design a landscape where the level of protection increases gradually as the habitat gets more sensitive, as illustrated in Figure 5. The IUCN protected area management categories provide good guidelines to define the policies associated with such configurations [Dudley, 2008]. In our CP model, expressing a nesting is easy. Let Γω be a neighborhood definition, {R0, ...,Rn-1} a partition, and {R0, ...,Rm-1} be the regions to be nested. We assume that R0 is the core region and Rm-1 the periphery region. The nesting can be expressed as follows: ∀i ∈ [0, m[, buffer [Γω](Ri, S \ (Ri ∪Ri+1), Ri+1). The term S \ (Ri ∪ Ri+1) represents the area outside Ri and Ri+1, the latter thus being a buffer separating Ri from this outside area, as Russian nested dolls.

4 Use Case

New Caledonia is a large archipelago located in the South Pacific and the smallest biodiversity hotspot in the world. New Caledonian terrestrial flora, notably, is distinguished by a high rate of endemism (one of the highest in the world), and the presence of a large number of ancient lineages. However, New Caledonian forests are, as are most tropical forests, endangered with surface loss and fragmentation. Regarding this, Ibanez et al. [2017] conducted a study in a 60km2 sensitive area located in the South of the main Island, "Grande Terre". In this section, we rely on the dataset from this study and the Species Distribution Models (SDMs) produced from [Pouteau et al., 2015; Schmitt et al., 2017]. The area is tessellated into a 46 × 75 square grid (a site is 1.7ha) were 5431 trees were identified among 97 communities, over 88 forest fragments. This area harbors 223 tree species for which 173 SDMs were produced. The 50 species without SDM were arbitrarily considered as endangered. According to the recent literature, this problem can be considered as a large one [Wang et al., 2018]. Our model was implemented 1 and ran with Choco-solver and its Choco-graph extension (for graph variables) [Prud'homme et al., 2017], on a Linux laptop with (Intel Core i5-5200U CPU 2.20GHz×4, 8GB RAM). We ran optimization scenarios under a time limit of 4h, focusing on the ability to find solutions for many constraint configurations rather than optimality proof. Numbers of solutions found and solving times for the first and best solutions found are provided in Table 1 (none were proven optimal).

1 The source code is available on GitHub as an open source project: https://github.com/dimitri-justeau/choco-reserve.

Figure 6: Use case scenarios best solutions mappings.

4.1 Original Scenario (SC1)

First, we show that our model is able to tackle the use case as defined by Justeau-Allaire et al.. This use case consists in highlighting a partitioning of the study area into two regions: R0 the reserve and R1 the out-reserve area. The reserve is the only constrained region: it must be composed of at most two CCs (Constraint D) with a surface area of at least 300ha (Constraint F). The total area of the reserve must not exceed 1000ha (Constraint E). Moreover, each endangered species must have at least one known occurrence in the reserve (Constraint A), every other species must be covered with a minimum probability of 0.8 (Constraint B), and at least 340ha of forest area must be covered (Constraint C). Finally, a set of sites that are lakes or mining areas cannot belong to the reserve. This constraint was added by restricting the domain of R0. We tried to minimize the total size of the reserve (SC1, Figure 6a). The solving time of our approach is comparable to [Justeau-Allaire et al., 2018], the best solution found within the time limit is even slightly better. If such a scenario provides useful insights, the produced delineation suffers from a major limitation: 99.3% of the selected sites are on the edge (i.e. adjacent to the out-reserve area).

4.2 Extended Scenarios (SC2 and SC3)

We show how the previous use case can be extended to a more realistic one including buffer zone constraints. To this end, we rely on Ibanez et al. conclusions on edgeeffects and their conservation implications in the studied area: "[...] the surrounding vegetation including secondary forest at the edge and the vegetation matrix should also be protected to promote the long process of forest extension and subsequently reduce edge-effects [...]". Accordingly, we refined the use case by defining three regions: R0 the core, R1 the buffer zone and R2 the out-reserve area, the protected area being R0 ∪ R1. Pursuant with Ibanez et al. results which suggest that the spatial extent of the edge-effect applies within the first 100m to 300m, two scenarios were considered: an optimistic one (SC2, Figure 6b) with a 130-180 m wide buffer zone buffer[Γ8](R0,R2,R1), and a pessimistic one (SC3, Figure 6c) with a 260-360 m wide buffer zone buffer[Γ28 ](R0,R2,R1) 2 . In these scenarios, the region size constraint was applied to the protected area and became conflicting with the number of CCs constraint. Relying on conservation scientists' recommendations, we relaxed it and defined another optimization objective: minimize |cc(R1)|, that is minimize the fragmentation.

In SC2, we could satisfy the constraints with 14 CCs and 37% of the protected area located in the core (R0). In SC3, the constraints could not be satisfied, as the core cannot cover 200 sites of forest without exceeding the total area constraint. We thus relaxed the forest covering constraint to the protected area. This relaxation lead to a 10 CCs solution with only 10.6% of the protected area located in the core.

2 Width is variable because of the square grid: diagonal is √2 larger than horizontal/vertical.

4.3 Final Scenario (SC4)

Finally, we suggest a trade-off between SC2 and SC3 with SC4 (Figure 6d). In this final scenario we defined four regions: R0 the core, R1 the inner buffer, R2 the outer buffer and R3 the out-reserve area, the protected area being R0 ∪ R1 ∪ R2. These three regions were defined as nested: buffer [Γ8](R0, R2 ∪R3, R1) ∧ buffer [Γ8](R1, R0 ∪R3, R2). The species coverage constraints were still restricted to the core, however, the forest coverage constraint was only relaxed to the core and the inner buffer (R0 ∪R1). According to conservation scientists' feedbacks on SC3 results, we changed the optimization objective to: maximize the core area. The best solution found is composed of 9 CCs with 22.6% of the protected area located in the core, which is significantly better than SC3.

5 Conclusion

In this paper, we introduced a generic CP model that is able to tackle a high variety of reserve selection and design problems by providing high levels of flexibility and expressiveness. It is the first approach to allow the definition of an arbitrary number of regions on top of which any coverage or spatial constraint can be explicitly expressed. In addition, we provided the first CP formulation of the buffer zone constraint, which is compatible with any neighborhood definition in the tessellated geographical space and can be reused to compose more complex spatial constraints. Moreover, we provided insights on the consistency associated with the constraint, as well as on its worst-time complexity. Relying on a use case based on a real-world dataset, we showed how our model is able to support systematic conservation planning through a progressive and exploratory process. Through diverse scenarios, we highlighted useful insights for managers and conservation scientists. In particular, we showed how the buffer zone constraint can be composed to prospect more complex conservation scenarios. On top of that, our implementation showed its ability to quickly find solutions to the decision problem (cf. Table 1), demonstrating its potential for exploring many scenarios.

To conclude, our constrained partitioning approach for reserve selection and design provides the basis of a generic and exploratory decision support tool for systematic conservation planning and computational sustainability. It now remains to work closely with conservation scientists and managers to refine it and integrate it in decisional processes, in order to move towards more sustainable land-use policies. Providing proofs of optimality is, on the other hand, a prospect for technical future work.