<figure><img src="image_1.png" /></figure>

HAL
open science

# Designing compact, connected and gap-free reserves with systematic reserve site selection models

Adrien Brunel, Jérémy Omer, Antoine Gicquel, Sophie Lanco

► To cite this version:

Adrien Brunel, Jérémy Omer, Antoine Gicquel, Sophie Lanco. Designing compact, connected and gap-free reserves with systematic reserve site selection models. Applied Mathematical Modelling, 2024, 134, pp.307-323. 10.1016/j.apm.2024.06.001. hal-04414078v3

HAL Id: hal-04414078

https://hal.science/hal-04414078v3

Submitted on 21 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Designing compact, connected and gap-free reserves with systematic reserve site selection models

Adrien Brunel<sup>a,*</sup>, Jérémy Omer<sup>b</sup>, Antoine Gicquel<sup>b</sup>, Sophie Lanco<sup>a</sup>

## ARTICLE INFO

**Keywords:**

Reserve site selection

Integer linear programming

Connectivity

Compactness

Gap-free

Fragmentation

## ABSTRACT

Protected areas play a crucial role in current global policies to mitigate the erosion of biodiversity and systematic reserve site selection models are increasingly involved in their design. These models address the optimisation problem that seeks to cover spaces hosting biodiversity features with nature reserves at a minimum cost for human activities. To increase the likelihood of a successful implementation, reserves need to be spatially consistent. Widely used decision support tools such as Marxan and PrioritizR commonly enforce compactness indirectly by penalising the reserve perimeter in the objective function. Few other optimisation models explicitly consider spatial properties such as limited fragmentation, connectivity of selected sites, and buffer zones around them, etc. So far, no reserve site selection model can guarantee the production of a connected, compact, and gap-free reserve all at once. The impossibility of designing reserve solutions with desirable spatial properties using existing models makes it difficult to implement such solutions in the real world. Therefore, we propose a mixed-integer linear program to build a reserve that is connected, compact, and gap-free. To enforce these spatial attributes within a reserve site selection model, we used a multicommodity flow approach. We tested the computational feasibility of our model on generated instances and the real instance of Fernando de Noronha. The results indicate that a single model can be used to enforce compactness, connectivity, and the absence of gaps. Using this optimisation model, conservation practitioners can design reserve solutions with desirable spatial properties, thereby increasing the likelihood of a successful implementation.

## 1. Introduction

Biodiversity and habitats are threatened worldwide [1]. Building comprehensive networks of protected areas has become a popular conservation solution [2–4] and was shown to bring conservation benefits [5–7]. At sea, for instance, current political objectives are to cover 30% of the marine spaces under jurisdiction by 2030 with marine protected areas [8–10]. Similar concerns also exist on land [11,12]. Within this context, there is a strong demand in these spaces to find the best compromises between the protection of biodiversity and the sustainability of human uses. To address such problems, optimisation methods are commonly implemented [13–16]. The objective of these optimisation models is to identify the most cost-effective configuration of planning units

<footnote><sup>a</sup> UMR MARBEC (IRD, Univ. Montpellier, Ifremer, CNRS), 34200, Sète, France</footnote>

<footnote><sup>b</sup> Univ Rennes, INSA Rennes, CNRS, IRMAR - UMR 6625, F-35000, Rennes, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail addresses: adrien.brunel.pro@protonmail.com (A. Brunel), Jeremy.Omer@insa-rennes.fr (J. Omer), antoine.a.gicquel@inria.fr (A. Gicquel), sophie.lanco@ird.fr (S. Lanco).</footnote>
that collectively encompass a sufficient representation of each conservation feature. Once this representation is deemed sufficient, the conservation features are considered protected within the reserve. Since optimisation models can be customised at will by setting the appropriate constraints in the model, it is possible to produce a reserve that exhibits a set of desirable properties. Reserve design efforts focused on finding these appropriate constraints and explicitly address many conservation aspects (species representation, abundance/diversity trade-off, compactness, buffer zones, corridors, uncertain representation, sequential acquisition planning units, robustness to catastrophic events, etc.). Such methods are often embedded within a software, e.g. Marxan or PrioritizR [17–19]. They are designed to systematically select reserve sites and are used as decision support tools in real-world instances<fnref n="1" /> [20,21].

Nevertheless, solving these optimisation problems often results in the selection of scattered reserve sites. Yet, designing reserves that are compact, connected, and gap-free is usually needed for ecological, management and enforcement reasons [14,22–24]. A reserve is connected if one can move anywhere inside without having to leave it. A gap within a reserve is a zone outside the reserve one cannot leave without crossing the reserve. Compactness is not a binary concept, since a reserve is said more or less compact. It often indicates to what extent reserve sites are more or less aggregated. Several measures of compactness exist such as the area-to-perimeter ratio, the maximum distance between two reserve sites, the number of shared edges between reserve sites, etc. An illustration of these spatial attributes can be found in Fig. 1. Currently, the spatial attributes of reserves are poorly considered in decision support tools used for reserve selection. In the widely used decision support tools for reserve selection (e.g. Marxan, PrioritizR), the only spatial attribute explicitly addressed is the global compactness of a solution [17–19,25]. The compactness of a solution is enforced by directly penalising the overall perimeter of the reserve in the objective function of the optimisation problem addressed. Several issues arise with this approach. The linearisation of the perimeter expression involves the addition of many decision variables and constraints [26,27] which can be computationally expensive in an integer programming context. Also, this approach transforms the problem into a multi-objective problem where the cost of a solution and its perimeter are implicitly competing. In practice, the compactness multiplier is determined by trial and error until a solution meets the spatial requirements deemed satisfactory. This weakens the systematic nature of the reserve design approach, although a more systematic setting of the compactness multiplier is proposed in [28]. Improvements using both the reserve perimeter and area in the objective were proposed in [28] to enforce the compactness of the reserve. In the same line, a weighted combination of both compactness and connectivity measures are included in the objective and solved using metaheuristics in [29]. In any case, the connectivity of the reserve and the absence of gaps within it are not ensured, but rather, possibly emerging with the enforcement of the reserve compactness.

In an operation research context, several optimisation models were proposed to explicitly account for specific spatial properties [26,30–33]. For instance, optimisation models aim to design a reserve composed of a core area surrounded by a buffer zone [34–36]. It means every planning unit of the designated core area is always surrounded by other core planning units or buffer planning units. In particular, [36] proposed a model accounting for both the buffer zone and the connectivity of the core area. But these models do not necessarily ensure the resulting reserve is gap-free although such reserve can eventually emerge from them. A large family of models takes advantage of the modelling possibilities offered by the use of pairwise distances between candidate sites. Minimizing the sum of pairwise distances or the maximum distance between all reserve sites [37] favours compact reserves, but does not guarantee that the reserve is connected and gap-free. The same applies to models that constrain two distinct sites containing the same conservation feature to be closer than a predefined threshold distance [38]. Another large family of optimisation models takes advantage of graph theory [32,39–43], in particular to explicitly ensure the reserve connectivity. However, the site selection may still result in the inclusion of gaps within the reserve solution, which we define as a set of isolated sites not assigned to the reserve and entirely disconnected from the outside (i.e. surrounded by the reserve). A reserve perforated by gaps cannot be used in a large-scale reserve design. If gaps appear in a solution proposed by a decision support tool, they will either be arbitrarily incorporated into the reserve, artificially connected to the outside (in either case, this will often lead to the use of suboptimal solutions), or the provided solution will be ignored. Using models imposing connectivity and promoting compactness is likely to favour gap-free reserves, but this is not guaranteed. For instance, it may be necessary to design nature reserves around areas that cannot be included in the reserve, such as a harbour or a trade route. These areas cannot be enclosed by the reserve and must remain accessible from the outside. State-of-the-art models often provide a reserve solution with gaps in these cases. Consideration of gaps within reserves is rarely addressed in the literature. The absence of gaps in the reserve can be a *posteriori* achieved by iteratively searching a gap-free reserve among slightly suboptimal solutions [32]. This model does not *a priori* prevent gaps from being included within the reserve, but rather hopes such a solution exists even if the objective value is degraded. Such a procedure is interesting but does not guarantee to have the connected, compact and gap-free reserve with the best objective value. A model selecting “cellularly” convex reserve solutions (also in regular grids) that are thus connected and gap-free is given in [44]. This model could be used to avoid reserves with gaps, but the convexity requirement may neglect some admissible connected and gap-free solutions that are not cellularly convex.

In this work, we propose a model that guarantees by construction to provide the best connected, compact and gap-free reserve. Our optimisation model enforces the connectivity of both the reserve and the non-reserve areas, resulting in a connected and gap-free reserve. The overall compactness of the reserve is shaped by specifying a maximum radius or a maximum perimeter of the reserve. We show the reserve solutions provided by our approach on the real-world instance of Fernando de Noronha. We numerically assess the generality of the proposed approaches on several generated instances made of 300 or 500 planning units and 3 conservation features. This work improves the current models towards reserve site selection models that explicitly design reserves with desirable

<footnote n="1">¹ More case study examples can be found at https://marxansolutions.org/community/ and https://prioritizr.net/articles/publication_record.html.</footnote>
<figure><img src="image_4.png" /><figcaption>(a) Example of a scattered selection of reserve sites.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(c) Example of a connected selection of reserve sites made of one component and one gap within the selection.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(b) Example of a compact selection of reserve sites made of two components.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(d) Example of a connected and gap-free selection of reserve sites.</figcaption></figure>

**Fig. 1.** Illustration of the spatial characteristics of a reserve sites selection: compactness, connectivity and gap-free. A reserve is compact if its planning units are aggregated. A reserve is connected if one can move anywhere inside without having to leave it. A gap within a reserve is a zone outside the reserve you cannot leave without crossing the reserve.

spatial properties: compact, connected and gap-free. Therefore, the model we propose increases the chances of conservation science being successfully implemented by conservation practitioners.

## 2. Models

Here we present the integer linear program for reserve site selection that ensures the reserve is connected, compact and gap-free. Before all, we recall the general formulation of the reserve site selection problem. Then, the multicommodity flow approach using graph theory is presented for both the reserve and non-reserve areas. An illustration example of a reserve selection optimisation modelled using graph theory can be found in Fig. 2. Finally, a reduction of the problem is proposed, as well as an approach limiting the maximum radius of the graph of the reserve. The graph of the reserve is the graph induced by the planning units selected in the reserve solution.

### 2.1. General formulation of the reserve site selection problem

In a reserve site selection problem, the study area is discretised into a set of $J$ planning units within which a set of $I$ conservation features are distributed. The amount of conservation feature $i$ in the planning unit $j$ is a real positive number denoted $a_{ij}$. Each planning unit has a cost $c_j$ usually understood as the socio-economic cost associated with the closure of this unit. The decision is about whether to include the planning unit in the reserve. Consequently, we associate the decision variables $x_j$ with each planning unit $j$: $x_j = 1$ if a planning unit $j$ belongs to the reserve and $x_j = 0$ otherwise. One then seeks to find the least cost collection of planning units covering a sufficient amount for each conservation feature. The covering of a conservation feature $i$ is considered sufficient if it exceeds a user-defined level noted $t_i$.

$$ \sum_{j \in J} a_{ij} x_j \geq t_i, \forall i \in I \qquad (1) $$

Mathematically speaking, the general problem of reserve site selection is expressed as the following integer linear program $P_N$:
<figure><img src="image_8.png" /><figcaption>(a) Standard representation of the grid used to model a reserve site selection optimisation problem. Green planning units represent the selected reserve sites.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(b) Graph representation of the grid used to model connectivity constraints. Green nodes represent the selected reserve sites.</figcaption></figure>

<caption>**Fig. 2.** Example of the graph of a connected and gap-free reserve solution. Red arrows show the spanning tree of the reserve. The centre of the reserve tree is node 7. The radius of the graph of the reserve is 3 (reached by the path 7 → 6 → 5 → 9). Yellow arrows show the spanning tree of the non-reserve. The centre of the non-reserve tree is node α = 0 representing the outside area.</caption>

$$P_N : \left\{ \begin{array}{ll} \min_{x} & \sum_{j \in J} c_j x_j \\ \text{s.t.} & (1) \\ & x_j \in \{0, 1\} \quad \forall j \in J \end{array} \right.$$

In state-of-the-art models, the reserve perimeter is included in the objective function, to favour aggregated reserve solutions since a small perimeter involves a compact reserve. The perimeter is computed as the total length of the boundaries between reserved and non-reserved planning units. To model this, the length of the shared boundary between planning units $j_1$ and $j_2$ is denoted $b_{j_1j_2}$. A parameter $\beta$ is used to set the importance of compactness relative to the total cost of site selection. The quadratic expression of the perimeter can be linearised [45,27] by replacing the product $x_{j_1}x_{j_2}$ with the new binary decision variable $z_{j_1j_2}$ and add the following set of constraints:

$$\left\{ \begin{array}{ll} z_{j_1 j_2} - x_{j_1} \le 0 & \forall j_1 \in J \\ z_{j_1 j_2} - x_{j_2} \le 0 & \forall j_2 \in J \\ z_{j_1 j_2} - x_{j_1} - x_{j_2} \ge -1 & \forall j_1 \in J, \forall j_2 \in J \end{array} \right. \qquad (2)$$

Finally, the general formulation of the reserve site selection problem results in the following integer linear program $P_{N+COMP}$:

$$P_{N+COMP} : \left\{ \begin{array}{ll} \min_{x,z} & \sum_{j \in J} c_j x_j + \beta \sum_{j_1 \in J} \sum_{j_2 \in J} b_{j_1 j_2} (x_{j_1} - z_{j_1 j_2}) \\ \text{s.t.} & (1), (2) \\ & x_j, z_{j_1 j_2} \in \{0, 1\} \quad \forall j, j_1, j_2 \in J \end{array} \right.$$

This combinatorial optimisation problem is a minimum set cover problem known to be NP-hard [46]. It is a non-convex problem due to the binary nature of the decision variables. Yet, it can be expressed as an integer linear program and known solvers (like Gurobi or Cbc) can solve it for realistic instances in a reasonable time.

## 2.2. Connectivity of the reserve

The grid, resulting from the discretisation of the study area into planning units, is seen as a graph, where each planning unit $j \in J$ represents a node in the graph. The set of nodes is $J$. Planning units sharing an edge in the grid are considered neighbours and thus involve an edge $e = (e_1, e_2)$ in the graph between the nodes $e_1, e_2 \in J$. The set of edges is noted $E$. The corresponding directed edges $(e_1 \to e_2)$ and $(e_2 \to e_1)$ are called arcs. The set of arcs is noted $A$. We then use a multicommodity flow model developed in [47]. The idea is to constrain every node selected in the reserve to have a flow going from the source to the sink. The source is the commodity $k \in K$, i.e. a selected node, and the sink is the root node of the spanning tree. Therefore, we build a path connecting every selected nodes $k \in K$ and the root node which is constrained to belong to the reserve. The reserve is thus ensured to be connected. In this model, the set of commodities is $K = J$. The selection of the root node $j \in J$ of the spanning tree associated with the reserve is represented by the binary decision variable $r_j \in \{0, 1\}$. The selection of an arc $a \in A$ in the spanning tree associated with the reserve is represented by the binary decision variable $u_a \in \{0, 1\}$. The activation of the flow of commodity $k \in K$ between the source node $k$ and the sink node (i.e. root of the spanning tree) along the arc $a \in A$ is represented by the binary decision variable $f_a^k \in \{0, 1\}$. Let $V(n)$ be the set of neighbours nodes of node $n \in J$.
The selected arcs of the spanning tree must be between two nodes selected in the reserve:

$$ \left\{ \begin{array}{ll} u_a \le x_{a_1} & \forall a = (a_1, a_2) \in A \\ u_a \le x_{a_2} & \forall a = (a_1, a_2) \in A \end{array} \right. \qquad (3) $$

A maximum of one arc is activated by edge:

$$ u_{(e_1 \to e_2)} + u_{(e_2 \to e_1)} \le 1 \quad \forall e = (e_1, e_2) \in E \qquad (4) $$

The number of arcs in the tree is equal to the number of nodes minus 1 (prevent cycle formation):

$$ \sum_{a \in A} u_a = \sum_{j \in J} x_j - 1 \qquad (5) $$

The root of the tree must be in the reserve:

$$ r_j \le x_j \quad \forall j \in J \qquad (6) $$

There is only one root node for the spanning tree of the reserve:

$$ \sum_{j \in J} r_j \le 1 \qquad (7) $$

If the arc is not selected, all associated flow variables are set to 0:

$$ f_a^k \le u_a \quad \forall a \in A, \forall k \in K \qquad (8) $$

If the node is not selected, all the associated flow variables are set to 0:

$$ f_a^k \le x_k \quad \forall a \in A, \forall k \in K \qquad (9) $$

For commodity $k \in K$, the flow at the source node is 1, the flow at the sink node is 0, elsewhere for selected nodes, the flow entering is the same as the flow leaving the node.

$$ \left\{ \begin{array}{ll} \sum_{j \in V(n)} f_{(j \to n)}^k - \sum_{j \in V(n)} f_{(n \to j)}^k \le r_n & \forall k \in K, \forall n \in J \setminus \{k\} \\ \sum_{j \in V(n)} f_{(j \to n)}^k - \sum_{j \in V(n)} f_{(n \to j)}^k \ge 0 & \forall k \in K, \forall n \in J \setminus \{k\} \\ \sum_{j \in V(k)} f_{(k \to j)}^k - \sum_{j \in V(k)} f_{(j \to k)}^k = x_k - r_k & \forall k \in K \end{array} \right. \qquad (10) $$

Finally, the multicommodity flow model $P_{CON}$ for the reserve is:

$$ P_{CON} : \left\{ \begin{array}{ll} \min_{x,z,u,r,f} & \sum_{j \in J} c_j x_j + \beta \sum_{j_1 \in J} \sum_{j_2 \in J} b_{j_1 j_2} (x_{j_1} - z_{j_1 j_2}) \\ \text{s.t.} & (1) - (10) \\ & x_j, z_{j_1 j_2}, u_a, r_j, f_a^k \in \{0, 1\} \quad \forall j, j_1, j_2 \in J, \forall a \in A, \forall k \in K \end{array} \right. $$

## 2.3. Gap-free reserve

We apply the same multicommodity flow model to the non-reserve to have a connected non-reserve. A connected non-reserve implies that the reserve would not have gaps within it. Thus, the term $1 - x_j$ plays the role of the term $x_j$. We add a fictive node $\alpha$ in the graph representing the area outside the studied zone. Indeed, the non-reserve must be connected to the exterior area. Note that we fix the node $\alpha$ to be the root of the spanning tree of the non-reserve. The selection of an arc $a \in A$ in the spanning tree associated with the non-reserve is represented by the binary decision variable $v_a \in \{0, 1\}$. The activation of the flow of commodity $k \in K$ between the source node $k$ and the sink node $\alpha$ along the arc $a \in A$ is represented by the binary decision variable $g_a^k \in \{0, 1\}$. The set of edges and arcs associated with the fictive node $\alpha$ are respectively noted $E^f$ and $A^f$. Let $E^+ = E \cup E^f$, $A^+ = A \cup A^f$, and $J^+ = J \cup \{\alpha\}$.
The added set of constraints for the non-reserve is:

$$
\left\{
\begin{array}{ll}
x_a = 0 & \\
v_a \le 1 - x_{a1} & \forall a = (a_1 \to a_2) \in A^+ \\
v_a \le 1 - x_{a2} & \forall a = (a_1 \to a_2) \in A^+ \\
v_{(e_1 - e_2)} + v_{(e_2 - e_1)} \le 1 & \forall e = (e_1, e_2) \in E^+ \\
\sum_{a \in A^+} v_a = \sum_{j \in J^+} (1 - x_j) - 1 & \\
g_a^k \le v_a & \forall a \in A, \forall k \in K \\
g_a^k \le 1 - x_k & \forall a \in A, \forall k \in K \\
\sum_{j \in V(a)} g_{(a \to j)}^k - \sum_{j \in V(a)} g_{(j \to a)}^k \le 1 & \forall k \in K \\
\sum_{j \in V(n)} g_{(j \to n)}^k - \sum_{j \in V(n)} g_{(n \to j)}^k = 0 & \forall k \in K, \forall n \in J \setminus \{k\} \\
\sum_{j \in V(k)} g_{(j \to k)}^k - \sum_{j \in V(k)} g_{(k \to j)}^k = 1 - x_k & \forall k \in K \\
v_a, g_a^k \in \{0, 1\} & \forall a \in A^+, \forall k \in K
\end{array}
\right.
\qquad (11)
$$

By adding (11) to P<sub>CON</sub>, we get the integer linear program P<sub>CON+GF</sub> that ensure connected reserve solutions to be gap-free.

2.4. Compactness of the reserve

2.4.1. Maximum radius in the graph of the reserve

We want to avoid producing connected reserve solutions that spread across the entire study area. We thus impose the radius of the graph of the reserve to remain below a predefined threshold, denoted $R_{max}$ in the following. We define the radius of the reserve as the longest path starting from the source node and only composed of reserve nodes (illustration in Fig. 2). We have a double expectation with this additional constraint. First, we will produce more compact reserves and avoid cobweb shapes for the reserves by limiting the radius of the reserve. Secondly, by removing the nodes further than $R_{max}$ from consideration, we decrease the number of possible paths that satisfy the multicommodity flow constraints. This feature may counterbalance the addition of the maximum radius constraints and may even increase the solving speed.

Once the graph of an incumbent reserve solution is connected, we can define the centre and the radius of the graph of the reserve.
The centre is the selected node whose maximal distance from other selected nodes is the smallest. The radius is the maximum distance
in the graph between the centre and other selected nodes. Let $d(j_1, j_2)$ define the distance in the graph of the reserve between the
node $j_1 \in J$ and $j_2 \in J$. This distance corresponds to the shortest path in the reserve graph from node $j_1$ to node $j_2$. Note that the
global matrix of distances between all nodes of the grid was computed outside the solving procedure. All the selected nodes of the
incumbent connected reserve that are at a distance greater than $R_{max}$ from the centre are added to the set of commodities $K$. Then,
the following constraint is applied:

$$
\sum_{\substack{j \in J \\ d(j,k) \le R_{\max}}} \sum_{n \in V(j)} f_{j \to n}^k \le R_{\max} \quad \forall k \in K \quad (12)
$$

Finally, we impose the non-selection of nodes at a distance greater than $R_{max}$ from the root of the tree of the reserve:

$$
x_{j_1} \le 1 - r_{j_2} \quad \forall j_1 \in J, j_2 \in J, d(j_1, j_2) > R_{\text{max}} \qquad (13)
$$

By adding (12) and (13) to P<sub>CON+GF</sub>, we get the integer linear program that ensures connected and gap-free solutions to have a
radius smaller than R<sub>max</sub>.

2.4.2. Maximum perimeter of the reserve

As explained in Section 2.1, the compactness of a reserve in state-of-the-art models is enforced using a multi-objective approach
by penalising the reserve perimeter in the objective. Rather than that, we can keep a single objective formulation and specify a
maximum perimeter P<sub>max</sub> the reserve should not exceed. The associated constraint is:

$$
\sum_{j_1 \in J} \sum_{j_2 \in J} b_{j_1 j_2} (x_{j_1} - z_{j_1 j_2}) \leq P_{\text{max}} \qquad (14)
$$

By adding (12) and (13) and/or (14) to P<sub>CON+GF</sub>, we get the integer linear program that ensures connected and gap-free solutions to
have a perimeter smaller than P<sub>max</sub>. The models that in addition include constraints used to enforce the compactness of the reserve
will be named P<sub>CON+GF+COMP</sub>. The simple mention of compactness parameters, i.e. β, R<sub>max</sub>, P<sub>max</sub>, will remove any ambiguity
regarding our choice of constraints to enforce the compactness of the reserve in P<sub>CON+GF+COMP</sub>. If β is mentioned, it means we
include the penalty of the reserve perimeter in the objective and add the associated linearisation constraints (2) to the model. If
R<sub>max</sub> is mentioned, the constraints (12) and (13) are added to the model. If P<sub>max</sub> is mentioned, the constraint (14) is added to the
model.
<figure><img src="image_10.png" /><figcaption>Fig. 3. Chessboard overlay applied to the graph of the grid. Nodes are separated into two sets: black and white nodes. The multicommodity flow model is then only applied to the black nodes to decrease the size of the problem.</figcaption></figure>

## 2.5. Improvements of the model

### 2.5.1. Chessboard reduction

In a rectangular grid, if we want a connected and gap-free reserve, a node in a given binary state, i.e. selected or unselected, cannot be entirely surrounded by neighbouring nodes in the complementary state. The rectangular grid is thus assimilated to a chessboard, and the nodes are separated into two sets: black and white nodes (Fig. 3). This way, the 4 neighbouring nodes of a black node are white and vice versa. Let B be the set of black nodes and W the set of white nodes. We have $J = W \cup B$. In terms of constraints, we prevent white (respectively black) nodes of the grid in a given state from being surrounded by four black (respectively white) neighbours in the same state:

$$ \left\{ \begin{array}{l} x_j \le \sum_{i \in V(j)} x_i \quad \forall j \in J \\ 1 - x_j \le \sum_{i \in V(j)} (1 - x_i) \quad \forall j \in J \end{array} \right. \qquad (15) $$

Then, we apply the multicommodity flow model only to black nodes: the flow is now constrained to find a path from a black source node to the black nodes only. It means that the set of commodities is $K = B$ instead of $K = J$ in model $P_{CON}$ or $P_{CON+GF}$. This is the main motivation behind this chessboard reduction: we significantly reduce the number of expensive multicommodity flow constraints by only adding two constraints by node. This way, each node is either associated with a commodity or satisfies constraint (15) and has its neighbours associated with a commodity. In the following, the chessboard reduction is systematically applied.

### 2.5.2. Lazy constraints

Enforcing flow constraints for every node can be computationally challenging in multicommodity flow models. Lazy constraints are included in the model only if the incumbent solution does not satisfy them. Since the flow constraints (10) in the multicommodity flow model can be separated by commodities, these constraints are implemented as lazy ones. The motivation behind this choice is that a non-connected reserve is not a frequent case. We expect the iterative activation of the lazy constraints to be faster than considering the exhaustive set of flow constraints. In the following, the concerned constraints are systematically applied as lazy constraints.

A graph is connected if there is a path from any point to any other point in the graph. If a graph is not connected, the graph is made of two or more isolated connected subgraphs. We define the connected components of the graph as the maximum connected subgraphs, i.e. the connected subgraph composed of the maximum number of nodes. In our case, if the number of connected components of a reserve solution is strictly greater than 1, the reserve is not connected and we activate the flow constraints associated with a given commodity noted $k_1 \in K$. We define the return of a node as the ratio between the conservation feature's total amounts within a planning unit and its cost. The commodity $k_1$ is chosen as the node with the highest return among the nodes of each connected component of the incumbent reserve solution.

We do the same for the multicommodity flow model of the non-reserve. If the number of connected components of the non-reserve graph is greater than 1, the incumbent reserve solution has a gap within it, and we activate the flow constraints associated with a given commodity noted $k_2 \in K$. The commodity $k_2$ is chosen as the node with the lowest return among each connected component of the non-reserve. As above, a gap within the reserve is expected to be a rare case, so the lazy constraints allow a faster solving than including the exhaustive set of flow constraints (11).

The constraints (12) are also implemented as lazy constraints, and thus activated only if the incumbent reserve solution is spreading too much.
<figure><img src="image_11.png" /><figcaption>**Fig. 4.** Real-world data used for the reserve site selection numerical experiment. (A) Fishery-based cost layer in a continuous orange colour gradient. (B) Fish biomass conservation feature surrogate in a discrete purple colour gradient. (C) Continental shelf and (D) Shelf break habitat conservation feature surrogates in light and deep blue respectively. Opaque grey pixels are the planning units a priori excluded from the solution.</figcaption></figure>

## 3. Numerical experiments

We performed numerical experiments to validate, illustrate and then assess the numerical feasibility of the models described in Section 2. Our tests were performed on generated instances (3 conservation features, 300 or 500 planning units) and a real-world instance (3 conservation features, 756 planning units). The instances can be found at https://github.com/AdrienBrunel/data_generation. The experiments were realised on a personal computer (Intel Core i7-8850H CPU @ 2.60 GHz). The code used for the analyses of this work is open, free and available at https://github.com/AdrienBrunel/rssp-compact-connected-gapfree. We used Gurobi solver under a free academic license using the JuMP optimisation library [48] of Julia [49,50].

### 3.1. *Illustration on a real-world instance*

Fernando de Noronha is a small oceanic archipelago in the western tropical Atlantic, made up of 21 islands, islets and rocks with a total land area of 26 km<fnref n="2" />. It constitutes a genuine Brazilian natural and cultural heritage and a conservation showcase in Brazil. But it also faces many interests (e.g. tourism intensification, fisheries), which makes it an open laboratory for marine spatial planning. We summarise the main characteristics of the dataset below (see [51] for further details).

The geographical area was discretised according to a rectangular grid made of N=36×21=756 planning units with longitude and latitude respectively in [32.65° W, 32.30° W] and [3.95°S, 3.75°S] ranges. Planning units covering Fernando de Noronha land and harbour were a priori excluded from potential reserve site candidates. The considered conservation features are the fish biomass, the continental shelf and shelf break habitats. The cost layer was made of the fishing pressure intensity. Fig. 4 shows the input data involved in the case study.

For the first conservation feature, the fish biomass was estimated from *in situ* acoustic data [52,53]. Interpolating between sample data allowed producing a continuous distribution within the sampling area. Outside the surveyed area, values were set to 0, although the actual distribution was unknown. Then, ocean depth intervals were used as a surrogate for the two other conservation features: the continental shelf and shelf break habitats. Ocean bathymetry was obtained from GEBCO online platform.<fnref n="2" /> Finally, a segmentation model was applied to fishers' trajectories to derive the behavioural state for every GPS measure: fishing or travelling. This was then used to derive a quantitative proxy for the fishing pressure. The fishing activity proxy represents the cost vector in the objective function of the optimisation problem.

<footnote n="2">² General Bathymetric Chart of the Ocean. https://download.gebco.net/</footnote>
<table><caption>**Table 1**

Results of the numerical experiments for 36 × 21 planning units and 3 conservation features for the real-world instance of Fernando de Noronha. A summary of the characteristics of the reserve solutions is provided: computation time in seconds, reserve perimeter, total cost, radius of the reserve graph, number of connected components, number of gaps. Targets are the same for the three conservation features (50% or 70%). N = nominal, COMP = compactness, CON = connectivity, GF = gap-free.</caption><thead><tr><th>Targets</th><th>Model</th><th>Parameters</th><th>Time</th><th>Perimeter</th><th>Cost</th><th>Radius</th><th>Components</th><th>Gaps</th></tr></thead><tbody><tr><td>50%</td><td>N</td><td>β = 0</td><td>0.2</td><td>150</td><td>90.6</td><td>-</td><td>14</td><td>4</td></tr><tr><td>50%</td><td>N+COMP</td><td>β = 1</td><td>0.2</td><td>96</td><td>101.7</td><td>-</td><td>2</td><td>0</td></tr><tr><td>50%</td><td>CON+COMP</td><td>β = 1</td><td>19.5</td><td>92</td><td>106.0</td><td>16</td><td>1</td><td>0</td></tr><tr><td>50%</td><td>CON+COMP</td><td>β = 1, R<sub>max</sub> = 14</td><td>920.5</td><td>90</td><td>112.4</td><td>14</td><td>1</td><td>1</td></tr><tr><td>50%</td><td>CON+COMP+GF</td><td>β = 1, R<sub>max</sub> = 14</td><td>179.6</td><td>90</td><td>112.8</td><td>14</td><td>1</td><td>0</td></tr><tr><td>50%</td><td>CON+COMP+GF</td><td>P<sub>max</sub> = 80, R<sub>max</sub> = 14</td><td>366.1</td><td>80</td><td>125.6</td><td>14</td><td>1</td><td>0</td></tr><tr><td>70%</td><td>N</td><td>β = 0</td><td>0.0</td><td>156</td><td>200.2</td><td>-</td><td>7</td><td>7</td></tr><tr><td>70%</td><td>N+COMP</td><td>β = 1</td><td>0.1</td><td>108</td><td>216.8</td><td>-</td><td>2</td><td>2</td></tr><tr><td>70%</td><td>CON+COMP</td><td>β = 1</td><td>69.4</td><td>100</td><td>227.6</td><td>16</td><td>1</td><td>2</td></tr><tr><td>70%</td><td>CON+COMP+GF</td><td>β = 1</td><td>42.6</td><td>98</td><td>232.8</td><td>18</td><td>1</td><td>0</td></tr><tr><td>70%</td><td>CON+COMP+GF</td><td>β = 1, R<sub>max</sub> = 17</td><td>295.1</td><td>96</td><td>235.0</td><td>17</td><td>1</td><td>0</td></tr><tr><td>70%</td><td>CON+COMP+GF</td><td>β = 1, R<sub>max</sub> = 16</td><td>367.6</td><td>94</td><td>237.6</td><td>16</td><td>1</td><td>0</td></tr><tr><td>70%</td><td>CON+COMP+GF</td><td>P<sub>max</sub> = 90, R<sub>max</sub> = 16</td><td>164.2</td><td>90</td><td>243.6</td><td>16</td><td>1</td><td>0</td></tr></tbody></table>

Table 1 provides the characteristics of the reserve solutions computed using the models described in Section 2 on the real-world instance of Fernando de Noronha. The first observation is that the spatial coherence of a reserve is not guaranteed by state-of-the-art models. Fig. 5a and Fig. 6a show that the reserve site selection is significantly scattered for $β = 0$. Setting $β = 1$ in these models improved the global compactness of the reserve selection as illustrated in Fig. 5b and Fig. 6b but did not guarantee the connectivity of the reserve nor the absence of gaps within it (cf. Table 1). The state-of-the-art model with $β = 1$ (Fig. 6b) illustrates the problem with high covering demands and locked-out planning units: it naturally creates gaps by surrounding the locked-out planning units. When targets were set to 50%, obtaining a connected and gap-free reserve (cf. Fig. 5c) took 19.5 seconds. It required solving the model $P_{CON+COMP}$ with $β = 1$ since the reserve solution did not have any internal gap. However, when we increased the compactness demand ($β = 1$ and $R_{max} = 14$), a gap appeared within the solution. We removed this by solving the model $P_{CON+GF+COMP}$. We also observed that the solving of $P_{CON+GF+COMP}$ took less time than $P_{CON+COMP}$ (179.6 against 920.5 seconds). To obtain an even more compact reserve, we directly constrained the reserve perimeter to remain below 80 instead of the 90 of the reserve solution with $β = 1$ and $R_{max} = 14$. The connected, compact and gap-free solution (cf. Fig. 5d) was obtained in 366.1 seconds. When targets were set to 70%, solving the model $P_{CON+COMP}$ with $β = 1$ did not prevent the occurrence of gaps within the reserve (cf. Table 1). The compact, connected and gap-free reserve solution obtained by solving $P_{CON+GF+COMP}$ with $β = 1$ is shown in Fig. 6c. Also, the selection of isolated planning units is tolerated as soon as the perimeter involved contributed less to the objective than the selection cost. Unlike the state-of-the-art models’ solutions, our connected and gap-free reserve solutions left a path from the harbour of Fernando de Noronha to the outside area. It took only 42.6 seconds in this example. Then we increased our compactness demands by setting $P_{max} = 90$ and $R_{max} = 16$. The corresponding solution is shown in Fig. 6d. This reserve which looks more compact than Fig. 6c and is still connected and gap-free.

## 3.2. Feasibility assessment on generated instances

In this section, we tested our models for several generated instances to have a more accurate idea of the computation time needed and the extra cost involved to obtain connected and gap-free reserve solutions. We used a systematic way of building instances for our reserve site selection optimisation problems. The main principle was to build realistic spatial distributions for conservation features. To do so, the amount of a conservation feature within a planning unit was randomly drawn in a Gaussian distribution whose mean value decreases with the distance to the closest randomly drawn epicentres. If no epicentre was provided, the mean value depended on the distance to the locked-out planning units supposed to represent a shoreline. The standard deviation of the Gaussian distribution was set to 20% of the maximum value for every planning unit and conservation feature. The cost was uniformly drawn. More details can be found in Appendix C in [54]. The generation of instances is different from what is done in [32]. We wanted more realistic instances, closer to a real dataset, and where solutions were more likely to have gaps when targets were high. Detailed numerical results are exhaustively provided in Table 2 and Table 3 in Appendix A.1.

In any case, the nominal problems $P_N$ and $P_{N+COMP}$ are solved very fast, mostly under 1 second and 3.4 seconds at worst. However, the reserve solutions with $β = 0$ are very scattered, with many gaps, for all instances: 23.6 connected components and 11.0 gaps on average for 300 planning units; 41.5 connected components and 21.2 gaps on average for 500 planning units. The reserve solutions with $β = 1$ are less scattered, but still have several connected components and gaps in general: 2.7 connected components and 2.2 gaps on average for 300 planning units; 5.2 connected components and 3.4 gaps on average for 500 planning units.

The reserve solutions using our complete model $P_{CON+GF+COMP}$ with $β = 1$ provide by construction connected and gap-free reserves. For instances of 300 planning units, the mean computation time is 80.0 seconds with a standard deviation of 61.4 seconds, a minimum and maximum time of 15.8 and 229.6 seconds respectively. Obtaining a connected and gap-free reserve solution involves a mean relative extra cost of 2.7% (standard deviation of 1.8%, maximum of 5.7%) with respect to the state-of-the-art model with the same value of $β$. When we match the compactness demand, i.e. we constrain the perimeter to remain below the perimeter of the
<figure><img src="image_12.png" /><figcaption>**Fig. 5.** Reserve solutions of the real-world instance of Fernando de Noronha for several models. Conservation features targets are all set to 50%. Green planning units represent the reserve selection. Grey planning units are a priori excluded. N = nominal, COMP = compactness, CON = connectivity, GF = gap-free. Barplots at the top left of each panel show the coverage of each conservation feature (FB = fish biomass, CS = continental shelf, SB = shelf break) by the reserve solution. The target is represented by the red line.</figcaption></figure>

state-of-the-art reserve solution, the mean relative extra cost drops to 0.7% (standard deviation of 0.4%, maximum of 1.3%). The mean computation time for this model is 150.6 seconds. Once the reserve was compact, connected and gap-free, we evaluated the impact of an increase in compactness using the $R_{max}$ constraints. To do so, we set $R_{max}$ to the reserve radius obtained by solving $P_{CON+GF+COMP}$ with $\beta = 1$ minus 1, so the constraints were activated. For instances of 300 planning units, the mean computation time is 182.3 seconds with a standard deviation of 280.4 seconds, a minimum and maximum time of 21.1 and 928.7 seconds. The inclusion of $R_{max}$ constraints was sometimes associated with a decrease in computation time, sometimes with an increase, depending on the instances considered.

For instances of 500 planning units, the proof of optimality was not provided every time with a time limit set to 1000 seconds. However, the incumbent solution returned was still compact, connected and gap-free. For instances reaching the time limit, we considered a solving time of 1000 seconds and the characteristics of the incumbent suboptimal reserve solution. The mean computation time of the reserve solutions using our complete model $P_{CON+GF+COMP}$ with $\beta = 1$ was 480.0 seconds with a standard deviation of 382.8 seconds, a minimum and maximum time of 50.8 and 1000 seconds respectively. Obtaining a connected and gap-free reserve solution involved a mean relative extra cost of 3.4% (standard deviation of 1.7%, maximum of 6.9%) with respect to the state-of-the-art model with the same value of $\beta$. When we matched the compactness demand, i.e. we constrained the perimeter to remain below the perimeter of the state-of-the-art reserve solution, the mean relative extra cost dropped to 0.9% (standard deviation of 0.4%, maximum of 1.6%). For instances of 500 planning units, when a maximum radius was imposed to the reserve, it led to increased computational difficulties. Unlike instances of 300 planning units, the addition of the $R_{max}$ constraints systematically involved a greater computation time or an incumbent solution further from optimality than without $R_{max}$. For three instances, no solutions were found within the time limit.

## 4. Discussion

In this work, we presented a global integer linear program that designs compact, connected and gap-free reserves. The example of Fernando de Noronha showed our model is operational in real-world instances. The reserve solutions of our model are by construction connected and gap-free and therefore exceed by far the spatial quality when compared to state-of-the-art solutions as, to
<figure><img src="image_13.png" /><figcaption>**Fig. 6.** Reserve solutions of the real-world instance of Fernando de Noronha for several models. Conservation features targets are all set to 70%. Green planning units represent the reserve selection. Grey planning units are *a priori* excluded. N = nominal, COMP = compactness, CON = connectivity, GF = gap-free. Barplots at the top left of each panel show the coverage of each conservation feature (FB = fish biomass, CS = continental shelf, SB = shelf break) by the reserve solution. The target is represented by the red line.</figcaption></figure>

our knowledge, current reserve design literature cannot guarantee gap-free spatial requirements. Our results showed that obtaining connected and gap-free reserves only occasioned extra costs of approximately 3% on average. Our model involved a larger computation time which remains completely acceptable in practice. Indeed, the numerical experiments performed on generated instances showed that our model remains mostly under the time limit of 1000 seconds for instances composed of 500 planning units and 3 conservation features. As reserve selections often exhibit a “cobweb” pattern in models that aim to enforce connectivity [24,39], it was necessary to constrain the selection compactness. Compactness in our model is highly customizable since it can be explicitly constrained through a maximum perimeter or radius, or implicitly constrained through a compactness multiplier, as in state-of-the-art models. The results indicated that incorporating compactness constraints with the $\beta$ multiplier and associated constraints or the $P_{max}$ constraint can help to find a solution faster. Additionally, to maintain explicit reserve design and a single objective optimisation, it may be preferable to limit the maximum perimeter of the reserve instead of penalising it in the objective function. Finally, the code used for this work is free, open and available.

Instead of adopting an *a posteriori* approach that removes reserve solutions with gaps from the search space [32], we used an *a priori* approach introducing specific constraints within the model to build a connected non-reserve as advocated in the discussion of [32]. A multicommodity flow model and lazy constraints were used to enforce connectivity in the minimum set problem. Using a multicommodity flow in reserve site selection problems is not new and was already mentioned in [33] but is not a common approach [24,32,38,39]. However, applying a multicommodity flow model to both the reserve and non-reserve areas is new.

The primary limitation of our model is the restricted size of instances that can be solved within a reasonable timeframe. As expected, the optimisation problem we proposed remains computationally challenging. Models that consider spatial constraints are generally more demanding computationally, especially as the problem size increases [55]. The size of instances considered here is of the same order of magnitude as many examples found in the literature. More precisely, the number of planning units considered in this work is similar to other existing works: 100 planning units in [24], 131 in [37], 225 in [33], 324 planning units [38], 391 planning units in [39], 400 in [32]). The number of conservation features considered in other works is however an order of magnitude beyond. It is however difficult to strictly compare instance sizes since the instances were not generated using the same methodology.
Our results showed that obtaining a reserve that is compact, connected and gap-free involved a small increase in the site selection cost with respect to state-of-the-art models used in decision support tools such as Marxan and PrioritizR. The only price to pay is a greater computation time. Therefore, obtaining spatially consistent reserves is a more computational challenge rather than a question of socio-economic resources. On the other hand, computational difficulties can potentially be mitigated by providing feasible solutions before reaching the optimality proof. This approach would result in faster provision of compact, connected, and gap-free reserves but with a slightly worse objective value. While some may argue that authorizing suboptimal solutions is no different from using a metaheuristics approach, it is important to note that our model allows for explicit control over the spatial properties of the reserve solution, which can be finely tuned. The search for a suboptimal reserve solution that meets predefined spatial requirements differs from relying on the serendipity of a metaheuristic search.

This work did not consider the potential differences in the nature of locked-out planning units. Some planning units should be allowed to be crossed, while others should not. For instance, in marine reserve design, a planning unit consisting of land is locked-out and cannot be crossed. Conversely, a locked-out planning unit located in the harbour can be crossed. This difference has an impact on modelling as a path that guarantees connectivity in flow models cannot cross every locked-out planning unit. This aspect is not currently included in our model and may occasion further development in the future.

In this work, we focused on producing compact, connected and gap-free reserves which are by construction single-component. However, single-component reserves lack robustness against spatially-correlated risks (e.g. fire, epidemic, oil spill, invasive species). Robustness against spatially-correlated risks may be achieved by implementing backup coverage of conservation features. For instance, ensuring a sufficient representation of conservation features across multiple components of the reserve could fulfil this objective. Other optimisation models addressed this specific aspect [38,56–58], although these models were mainly based on presence/absence of features and should be adapted to non-binary data. Future research should build upon these different approaches to produce multiple-component reserves where each component is compact, connected and gap-free.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Data availability

Data will be made available on request.

## Acknowledgements

Adrien Brunel was funded by Région Occitanie and Créocean, Grant No. ALDOCT-000674. This work is a contribution to the PADDLE project (European Union's Horizon 2020 Research and Innovation program under Grant Agreement No. 734271), TRIATLAS project (European Union's Horizon 2020 Research and Innovation program under Grant Agreement No. 817578), to the Mixed International Laboratory TAPIOCA, and to the Young Team associated with IRD (JEAI) TABASCO. We also thank Tatiana Beltrão, Maria Ester Simoni and Paulo Travassos for the fishing data collection and processing.

## Appendix A

### A.1. Feasibility assessment on generated instances

In this section, we provide the detailed tables referred to in Section 3.2. Table 2 and Table 3 provide the characteristics of the reserve solutions computed using the models described in Section 2 for generated instances of respectively 300 and 500 planning units.

### A.2. Assessment of the compactness models

In this section, we aim to assess the difference in computation time needed to obtain a compact solution whether using constraints associated with the use of $\beta$ or $P_{max}$. To do so, we first solved $P_{CON+GF+COMP}$ with $\beta$ = 1. The perimeter of the reserve solution obtained was then used for $P_{max}$ when we solved $P_{CON+GF+COMP}$. As expected, we obtained the same solutions between the two models. For instances of 300 planning units, the mean computation time was 80.0 seconds for $P_{CON+GF+COMP}$ with $\beta$ and 74.8 seconds with $P_{CON+GF+COMP}$ with $P_{max}$. Then, Table 4 does not show a systematic trend between the two models since it sometimes took more time, sometimes less time, depending on the instance. However, for instances of 500 planning units, Table 5 shows a clear trend: models with constraints associated to the use of $\beta$ are solved faster than the models using the $P_{max}$ constraint for every instance. When the time limit was reached, the model using the constraints associated with $\beta$ provided a solution closer to optimality than the model with the $P_{max}$ constraint.
<table><caption>**Table 2**

Results for 10 generated instances of 20 × 15 planning units and 3 conservation features. A summary of the characteristics of the reserve solutions is provided: computation time in seconds, reserve perimeter, total cost, radius of the reserve graph, number of connected components, number of gaps. Conservation features targets are all set to 50%. N = nominal, COMP = compactness, CON = connectivity, GF = gap-free.</caption><thead><tr><th>Instance</th><th>Model</th><th>Parameters</th><th>Time</th><th>Perimeter</th><th>Cost</th><th>Radius</th><th>Components</th><th>Gaps</th></tr></thead><tbody><tr><td>1</td><td>N</td><td>β = 0</td><td>0.1</td><td>272</td><td>460.1</td><td>-</td><td>21</td><td>10</td></tr><tr><td>1</td><td>N+COMP</td><td>β = 1</td><td>0.4</td><td>106</td><td>544.3</td><td>-</td><td>2</td><td>3</td></tr><tr><td>1</td><td>CON+GF+COMP</td><td>β = 1</td><td>15.8</td><td>76</td><td>575.5</td><td>12</td><td>1</td><td>0</td></tr><tr><td>1</td><td>CON+GF+COMP</td><td>P<sub>max</sub> = 106</td><td>245.6</td><td>106</td><td>550.2</td><td>13</td><td>1</td><td>0</td></tr><tr><td>1</td><td>CON+GF+COMP</td><td>β = 1, R<sub>max</sub> = 11</td><td>24.9</td><td>82</td><td>571.4</td><td>11</td><td>1</td><td>0</td></tr><tr><td>2</td><td>N</td><td>β = 0</td><td>0.2</td><td>286</td><td>480.9</td><td>-</td><td>25</td><td>14</td></tr><tr><td>2</td><td>N+COMP</td><td>β = 1</td><td>0.4</td><td>96</td><td>580.2</td><td>-</td><td>1</td><td>3</td></tr><tr><td>2</td><td>CON+GF+COMP</td><td>β = 1</td><td>80.1</td><td>86</td><td>592.5</td><td>13</td><td>1</td><td>0</td></tr><tr><td>2</td><td>CON+GF+COMP</td><td>P<sub>max</sub> = 96</td><td>143.5</td><td>96</td><td>582.7</td><td>14</td><td>1</td><td>0</td></tr><tr><td>2</td><td>CON+GF+COMP</td><td>β = 1, R<sub>max</sub> = 12</td><td>30.9</td><td>88</td><td>590.9</td><td>12</td><td>1</td><td>0</td></tr><tr><td>3</td><td>N</td><td>β = 0</td><td>0.1</td><td>296</td><td>466.4</td><td>-</td><td>26</td><td>9</td></tr><tr><td>3</td><td>N+COMP</td><td>β = 1</td><td>0.2</td><td>112</td><td>541.2</td><td>-</td><td>2</td><td>4</td></tr><tr><td>3</td><td>CON+GF+COMP</td><td>β = 1</td><td>124.0</td><td>102</td><td>554.9</td><td>12</td><td>1</td><td>0</td></tr><tr><td>3</td><td>CON+GF+COMP</td><td>P<sub>max</sub> = 112</td><td>229.9</td><td>112</td><td>546.3</td><td>12</td><td>1</td><td>0</td></tr><tr><td>3</td><td>CON+GF+COMP</td><td>β = 1, R<sub>max</sub> = 11</td><td>21.1</td><td>86</td><td>572.3</td><td>11</td><td>1</td><td>0</td></tr><tr><td>7</td><td>N</td><td>β = 0</td><td>0.1</td><td>276</td><td>454.0</td><td>-</td><td>19</td><td>4</td></tr><tr><td>7</td><td>N+COMP</td><td>β = 1</td><td>0.1</td><td>142</td><td>527.9</td><td>-</td><td>4</td><td>1</td></tr><tr><td>7</td><td>CON+GF+COMP</td><td>β = 1</td><td>40.2</td><td>132</td><td>540.4</td><td>15</td><td>1</td><td>0</td></tr><tr><td>7</td><td>CON+GF+COMP</td><td>P<sub>max</sub> = 142</td><td>10.1</td><td>142</td><td>530.9</td><td>16</td><td>1</td><td>0</td></tr><tr><td>7</td><td>CON+GF+COMP</td><td>β = 1, R<sub>max</sub> = 14</td><td>101.3</td><td>130</td><td>543.0</td><td>14</td><td>1</td><td>0</td></tr><tr><td>12</td><td>N</td><td>β = 0</td><td>0.1</td><td>284</td><td>491.9</td><td>-</td><td>23</td><td>14</td></tr><tr><td>12</td><td>N</td><td>β = 1</td><td>0.1</td><td>112</td><td>572.1</td><td>-</td><td>2</td><td>2</td></tr><tr><td>12</td><td>CON+GF+COMP</td><td>β = 1</td><td>19.8</td><td>110</td><td>575.7</td><td>17</td><td>1</td><td>0</td></tr><tr><td>12</td><td>CON+GF+COMP</td><td>P<sub>max</sub> = 112</td><td>44.3</td><td>112</td><td>574.4</td><td>17</td><td>1</td><td>0</td></tr><tr><td>12</td><td>CON+GF+COMP</td><td>β = 1, R<sub>max</sub> = 16</td><td>345.5</td><td>106</td><td>580.0</td><td>16</td><td>1</td><td>0</td></tr><tr><td>13</td><td>N</td><td>β = 0</td><td>0.4</td><td>268</td><td>464.9</td><td>-</td><td>21</td><td>10</td></tr><tr><td>13</td><td>N+COMP</td><td>β = 1</td><td>0.5</td><td>106</td><td>555.0</td><td>-</td><td>1</td><td>2</td></tr><tr><td>13</td><td>CON+GF+COMP</td><td>β = 1</td><td>69.6</td><td>106</td><td>555.9</td><td>12</td><td>1</td><td>0</td></tr><tr><td>13</td><td>CON+GF+COMP</td><td>P<sub>max</sub> = 106</td><td>41.2</td><td>106</td><td>555.9</td><td>12</td><td>1</td><td>0</td></tr><tr><td>13</td><td>CON+GF+COMP</td><td>β = 1, R<sub>max</sub> = 11</td><td>55.3</td><td>92</td><td>570.9</td><td>11</td><td>1</td><td>0</td></tr><tr><td>14</td><td>N</td><td>β = 0</td><td>0.1</td><td>286</td><td>452.7</td><td>-</td><td>30</td><td>15</td></tr><tr><td>14</td><td>N+COMP</td><td>β = 1</td><td>0.2</td><td>120</td><td>524.3</td><td>-</td><td>2</td><td>4</td></tr><tr><td>14</td><td>CON+GF+COMP</td><td>β = 1</td><td>65.9</td><td>102</td><td>548.2</td><td>12</td><td>1</td><td>0</td></tr><tr><td>14</td><td>CON+GF+COMP</td><td>P<sub>max</sub> = 120</td><td>307.2</td><td>120</td><td>531.2</td><td>12</td><td>1</td><td>0</td></tr><tr><td>14</td><td>CON+GF+COMP</td><td>β = 1, R<sub>max</sub> = 11</td><td>104.8</td><td>112</td><td>539.8</td><td>11</td><td>1</td><td>0</td></tr><tr><td>15</td><td>N</td><td>β = 0</td><td>0.2</td><td>302</td><td>490.4</td><td>-</td><td>25</td><td>14</td></tr><tr><td>15</td><td>N+COMP</td><td>β = 1</td><td>0.2</td><td>148</td><td>565.4</td><td>-</td><td>6</td><td>1</td></tr><tr><td>15</td><td>CON+GF+COMP</td><td>β = 1</td><td>76.5</td><td>136</td><td>582.4</td><td>18</td><td>1</td><td>0</td></tr><tr><td>15</td><td>CON+GF+COMP</td><td>P<sub>max</sub> = 148</td><td>78.2</td><td>148</td><td>571.1</td><td>19</td><td>1</td><td>0</td></tr><tr><td>15</td><td>CON+GF+COMP</td><td>β = 1, R<sub>max</sub> = 17</td><td>928.7</td><td>138</td><td>582.1</td><td>16</td><td>1</td><td>0</td></tr><tr><td>16</td><td>N</td><td>β = 0</td><td>0.1</td><td>280</td><td>472.3</td><td>-</td><td>21</td><td>5</td></tr><tr><td>16</td><td>N+COMP</td><td>β = 1</td><td>0.7</td><td>124</td><td>567.6</td><td>-</td><td>4</td><td>0</td></tr><tr><td>16</td><td>CON+GF+COMP</td><td>β = 1</td><td>229.6</td><td>120</td><td>573.8</td><td>17</td><td>1</td><td>0</td></tr><tr><td>16</td><td>CON+GF+COMP</td><td>P<sub>max</sub> = 124</td><td>191.6</td><td>124</td><td>569.9</td><td>16</td><td>1</td><td>0</td></tr><tr><td>16</td><td>CON+GF+COMP</td><td>β = 1, R<sub>max</sub> = 16</td><td>175.3</td><td>118</td><td>575.8</td><td>16</td><td>1</td><td>0</td></tr><tr><td>17</td><td>N</td><td>β = 0</td><td>0.0</td><td>274</td><td>418.2</td><td>-</td><td>25</td><td>15</td></tr><tr><td>17</td><td>N+COMP</td><td>β = 1</td><td>1.0</td><td>114</td><td>477.8</td><td>-</td><td>3</td><td>2</td></tr><tr><td>17</td><td>CON+GF+COMP</td><td>β = 1</td><td>78.7</td><td>94</td><td>499.4</td><td>12</td><td>1</td><td>0</td></tr><tr><td>17</td><td>CON+GF+COMP</td><td>P<sub>max</sub> = 114</td><td>214.0</td><td>114</td><td>483.3</td><td>14</td><td>1</td><td>0</td></tr><tr><td>17</td><td>CON+GF+COMP</td><td>β = 1, R<sub>max</sub> = 11</td><td>35.3</td><td>84</td><td>515.6</td><td>11</td><td>1</td><td>0</td></tr></tbody></table>

## References

[1] IPBES, Global assessment report on biodiversity and ecosystem services of the intergovernmental science-policy platform on biodiversity and ecosystem services, Tech. Rep., Zenodo, version Number: 1 https://doi.org/10.5281/ZENODO.3831673, May 2019.

[2] P.C. Ticco, The use of marine protected areas to preserve and enhance marine biological diversity: a case study approach, Coast. Manage. 23 (4) (1995) 309–314, https://doi.org/10.1080/08920759509362274.

[3] M. Tundi Agardy, Advances in marine conservation: the role of marine protected areas, Trends Ecol. Evol. 9 (7) (1994) 267–270, https://doi.org/10.1016/0169-5347(94)90297-6.
<caption>Table 3</caption>

Results for 10 instances of 25 × 20 planning units and 3 conservation features. A summary of the characteristics of the reserve solutions is provided: computation time in seconds, reserve perimeter, total cost, radius of the reserve graph, number of connected components, number of gaps. Conservation features targets are all set to 50%. N = nominal, COMP = compactness, CON = connectivity, GF = gap-free. If the time limit is reached (TL = 1000 s), the optimality gap of the incumbent solution is given within brackets.

<table><thead><tr><th>Instance</th><th>Model</th><th>Parameters</th><th>Time</th><th>Perimeter</th><th>Cost</th><th>Radius</th><th>Components</th><th>Gaps</th></tr></thead><tbody><tr><td>4</td><td>N</td><td>β = 0</td><td>0.1</td><td>524</td><td>772.6</td><td>-</td><td>38</td><td>28</td></tr><tr><td>4</td><td>N+COMP</td><td>β = 1</td><td>0.4</td><td>248</td><td>925.0</td><td>-</td><td>3</td><td>7</td></tr><tr><td>4</td><td>CON+GF+COMP</td><td>β = 1</td><td>139.0</td><td>214</td><td>961.8</td><td>21</td><td>1</td><td>0</td></tr><tr><td>4</td><td>CON+GF+COMP</td><td>P<sub>max</sub> = 248</td><td>993.1</td><td>248</td><td>930.8</td><td>22</td><td>1</td><td>0</td></tr><tr><td>4</td><td>CON+GF+COMP</td><td>β = 1, R<sub>max</sub> = 20</td><td>284.3</td><td>204</td><td>973.4</td><td>20</td><td>1</td><td>0</td></tr><tr><td>5</td><td>N</td><td>β = 0</td><td>0.2</td><td>474</td><td>663.6</td><td>-</td><td>39</td><td>17</td></tr><tr><td>5</td><td>N+COMP</td><td>β = 1</td><td>0.8</td><td>224</td><td>778.9</td><td>-</td><td>7</td><td>4</td></tr><tr><td>5</td><td>CON+GF+COMP</td><td>β = 1</td><td>TL[0.5%]</td><td>222</td><td>792.8</td><td>21</td><td>1</td><td>0</td></tr><tr><td>5</td><td>CON+GF+COMP</td><td>P<sub>max</sub> = 224</td><td>TL[0.7%]</td><td>224</td><td>791.6</td><td>21</td><td>1</td><td>0</td></tr><tr><td>5</td><td>CON+GF+COMP</td><td>β = 1, R<sub>max</sub> = 20</td><td>TL[Infeasible]</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>6</td><td>N</td><td>β = 0</td><td>0.6</td><td>494</td><td>742.1</td><td>-</td><td>43</td><td>16</td></tr><tr><td>6</td><td>N+COMP</td><td>β = 1</td><td>3.4</td><td>248</td><td>882.8</td><td>-</td><td>6</td><td>2</td></tr><tr><td>6</td><td>CON+GF+COMP</td><td>β = 1</td><td>577.4</td><td>224</td><td>910.8</td><td>23</td><td>1</td><td>0</td></tr><tr><td>6</td><td>CON+GF+COMP</td><td>P<sub>max</sub> = 248</td><td>969.8</td><td>248</td><td>889.1</td><td>23</td><td>1</td><td>0</td></tr><tr><td>6</td><td>CON+GF+COMP</td><td>β = 1, R<sub>max</sub> = 22</td><td>TL[Infeasible]</td><td>-</td><td>-</td><td>-</td><td>1</td><td>0</td></tr><tr><td>8</td><td>N</td><td>β = 0</td><td>0.1</td><td>476</td><td>749.8</td><td>-</td><td>43</td><td>20</td></tr><tr><td>8</td><td>N+COMP</td><td>β = 1</td><td>0.2</td><td>196</td><td>880.8</td><td>-</td><td>5</td><td>3</td></tr><tr><td>8</td><td>CON+GF+COMP</td><td>β = 1</td><td>50.8</td><td>164</td><td>915.9</td><td>21</td><td>1</td><td>0</td></tr><tr><td>8</td><td>CON+GF+COMP</td><td>P<sub>max</sub> = 196</td><td>105.1</td><td>196</td><td>886.4</td><td>21</td><td>1</td><td>0</td></tr><tr><td>8</td><td>CON+GF+COMP</td><td>β = 1, R<sub>max</sub> = 20</td><td>179.7</td><td>166</td><td>914.1</td><td>20</td><td>1</td><td>0</td></tr><tr><td>9</td><td>N</td><td>β = 0</td><td>0.1</td><td>530</td><td>806.8</td><td>-</td><td>43</td><td>25</td></tr><tr><td>9</td><td>N+COMP</td><td>β = 1</td><td>1.4</td><td>228</td><td>957.5</td><td>-</td><td>3</td><td>4</td></tr><tr><td>9</td><td>CON+GF+COMP</td><td>β = 1</td><td>542.6</td><td>192</td><td>999.5</td><td>21</td><td>1</td><td>0</td></tr><tr><td>9</td><td>CON+GF+COMP</td><td>P<sub>max</sub> = 228</td><td>717.2</td><td>228</td><td>965.1</td><td>24</td><td>1</td><td>0</td></tr><tr><td>9</td><td>CON+GF+COMP</td><td>β = 1, R<sub>max</sub> = 20</td><td>TL[Infeasible]</td><td>-</td><td>-</td><td>-</td><td>1</td><td>0</td></tr><tr><td>10</td><td>N</td><td>β = 0</td><td>0.4</td><td>460</td><td>773.3</td><td>-</td><td>47</td><td>27</td></tr><tr><td>10</td><td>N+COMP</td><td>β = 1</td><td>2.6</td><td>194</td><td>912.6</td><td>-</td><td>3</td><td>3</td></tr><tr><td>10</td><td>CON+GF+COMP</td><td>β = 1</td><td>97.6</td><td>184</td><td>927.8</td><td>18</td><td>1</td><td>0</td></tr><tr><td>10</td><td>CON+GF+COMP</td><td>P<sub>max</sub> = 194</td><td>TL[0.1%]</td><td>194</td><td>919.1</td><td>20</td><td>1</td><td>0</td></tr><tr><td>10</td><td>CON+GF+COMP</td><td>β = 1, R<sub>max</sub> = 17</td><td>263.2</td><td>174</td><td>938.4</td><td>17</td><td>1</td><td>0</td></tr><tr><td>18</td><td>N</td><td>β = 0</td><td>0.2</td><td>494</td><td>771.8</td><td>-</td><td>41</td><td>17</td></tr><tr><td>18</td><td>N+COMP</td><td>β = 1</td><td>0.6</td><td>230</td><td>913.4</td><td>-</td><td>8</td><td>5</td></tr><tr><td>18</td><td>CON+GF+COMP</td><td>β = 1</td><td>TL[0.4%]</td><td>178</td><td>976.6</td><td>23</td><td>1</td><td>0</td></tr><tr><td>18</td><td>CON+GF+COMP</td><td>P<sub>max</sub> = 230</td><td>TL[0.5%]</td><td>230</td><td>927.0</td><td>22</td><td>1</td><td>0</td></tr><tr><td>18</td><td>CON+GF+COMP</td><td>β = 1, R<sub>max</sub> = 22</td><td>TL[4.1%]</td><td>222</td><td>976.8</td><td>21</td><td>1</td><td>0</td></tr><tr><td>19</td><td>N</td><td>β = 0</td><td>0.1</td><td>460</td><td>807.9</td><td>-</td><td>30</td><td>15</td></tr><tr><td>19</td><td>N+COMP</td><td>β = 1</td><td>0.5</td><td>210</td><td>940.9</td><td>-</td><td>6</td><td>2</td></tr><tr><td>19</td><td>CON+GF+COMP</td><td>β = 1</td><td>109.2</td><td>202</td><td>956.2</td><td>23</td><td>1</td><td>0</td></tr><tr><td>19</td><td>CON+GF+COMP</td><td>P<sub>max</sub> = 210</td><td>TL[0.3%]</td><td>210</td><td>951.5</td><td>22</td><td>1</td><td>0</td></tr><tr><td>19</td><td>CON+GF+COMP</td><td>β = 1, R<sub>max</sub> = 22</td><td>435.2</td><td>202</td><td>956.9</td><td>19</td><td>1</td><td>0</td></tr><tr><td>20</td><td>N</td><td>β = 0</td><td>0.1</td><td>506</td><td>800.6</td><td>-</td><td>50</td><td>26</td></tr><tr><td>20</td><td>N+COMP</td><td>β = 1</td><td>1.3</td><td>198</td><td>950.4</td><td>-</td><td>6</td><td>2</td></tr><tr><td>20</td><td>CON+GF+COMP</td><td>β = 1</td><td>895.1</td><td>184</td><td>969.5</td><td>20</td><td>1</td><td>0</td></tr><tr><td>20</td><td>CON+GF+COMP</td><td>P<sub>max</sub> = 198</td><td>TL[0.3%]</td><td>198</td><td>957.9</td><td>24</td><td>1</td><td>0</td></tr><tr><td>20</td><td>CON+GF+COMP</td><td>β = 1, R<sub>max</sub> = 19</td><td>TL[0.7%]</td><td>174</td><td>984.2</td><td>19</td><td>1</td><td>0</td></tr><tr><td>21</td><td>N</td><td>β = 0</td><td>0.3</td><td>470</td><td>736.0</td><td>-</td><td>41</td><td>21</td></tr><tr><td>21</td><td>N+COMP</td><td>β = 1</td><td>0.6</td><td>198</td><td>870.2</td><td>-</td><td>5</td><td>2</td></tr><tr><td>21</td><td>CON+GF+COMP</td><td>β = 1</td><td>388.3</td><td>166</td><td>905.8</td><td>18</td><td>1</td><td>0</td></tr><tr><td>21</td><td>CON+GF+COMP</td><td>P<sub>max</sub> = 198</td><td>616.7</td><td>198</td><td>875.6</td><td>19</td><td>1</td><td>0</td></tr><tr><td>21</td><td>CON+GF+COMP</td><td>β = 1, R<sub>max</sub> = 17</td><td>TL[0.2%]</td><td>162</td><td>911.4</td><td>17</td><td>1</td><td>0</td></tr></tbody></table>

[4] S. Le Saout, M. Hoffmann, Y. Shi, A. Hughes, C. Bernard, T.M. Brooks, B. Bertzky, S.H.M. Butchart, S.N. Stuart, T. Badman, A.S. Rodrigues, Protected areas and effective biodiversity conservation, Science 342 (6160) (2013) 803–805, https://doi.org/10.1126/science.1239268.

[5] J. Claudet, C.W. Osenberg, L. Benedetti-Cecchi, P. Domenici, J.-A. García-Charton, A. Pérez-Ruzafa, F. Badalamenti, J. Bayle-Sempere, A. Brito, F. Bulleri, J.-M. Culioili, M. Dimech, J.M. Falcón, I. Guala, M. Milazzo, J. Sánchez-Meca, P.J. Somerfield, B. Stobart, F. Vandeperre, C. Valle, S. Planes, Marine reserves: size and age do matter, Ecol. Lett. 11 (5) (2008) 481–489, https://doi.org/10.1111/j.1461-0248.2008.01166.x.

[6] S. Stolton, N. Dudley, Arguments for Protected Areas: Multiple Benefits for Conservation Use, earthscan edition, Earthscan, 2010.

[7] P. Liu, S. Jiang, L. Zhao, Y. Li, P. Zhang, L. Zhang, What are the benefits of strictly protected nature reserves? Rapid assessment of ecosystem service values in Wanglang, Nature Reserve, China, Ecosyst. Serv. 26 (2017) 70–78, https://doi.org/10.1016/j.ecoser.2017.05.014.
<table><caption>**Table 4**
Assessment of the computation time needed to enforce compactness whether using $\beta$ or $P_{max}$ constraints for 10 instances of 20 × 15 planning units and 3 conservation features. A summary of the characteristics of the reserve solutions is provided: computation time in seconds, reserve perimeter, total cost, radius of the reserve graph, number of connected components, number of gaps. Conservation features targets are all set to 50%.</caption><thead><tr><th>Instance</th><th>Model</th><th>Parameters</th><th>Time</th><th>Perimeter</th><th>Cost</th><th>Radius</th><th>Components</th><th>Gaps</th></tr></thead><tbody><tr><td>1</td><td>CON+GF+COMP</td><td>β = 1</td><td>15.8</td><td>76</td><td>575.5</td><td>12</td><td>1</td><td>0</td></tr><tr><td>1</td><td>CON+GF+COMP</td><td>P<sub>max</sub> = 76</td><td>12.0</td><td>76</td><td>575.5</td><td>12</td><td>1</td><td>0</td></tr><tr><td>2</td><td>CON+GF+COMP</td><td>β = 1</td><td>80.1</td><td>86</td><td>592.5</td><td>13</td><td>1</td><td>0</td></tr><tr><td>2</td><td>CON+GF+COMP</td><td>P<sub>max</sub> = 86</td><td>83.6</td><td>86</td><td>592.5</td><td>13</td><td>1</td><td>0</td></tr><tr><td>3</td><td>CON+GF+COMP</td><td>β = 1</td><td>124.0</td><td>102</td><td>554.9</td><td>12</td><td>1</td><td>0</td></tr><tr><td>3</td><td>CON+GF+COMP</td><td>P<sub>max</sub> = 102</td><td>47.4</td><td>102</td><td>554.9</td><td>12</td><td>1</td><td>0</td></tr><tr><td>7</td><td>CON+GF+COMP</td><td>β = 1</td><td>40.2</td><td>132</td><td>540.4</td><td>15</td><td>1</td><td>0</td></tr><tr><td>7</td><td>CON+GF+COMP</td><td>P<sub>max</sub> = 132</td><td>35.8</td><td>132</td><td>540.4</td><td>15</td><td>1</td><td>0</td></tr><tr><td>12</td><td>CON+GF+COMP</td><td>β = 1</td><td>19.8</td><td>110</td><td>575.7</td><td>17</td><td>1</td><td>0</td></tr><tr><td>12</td><td>CON+GF+COMP</td><td>P<sub>max</sub> = 110</td><td>118.0</td><td>110</td><td>575.7</td><td>17</td><td>1</td><td>0</td></tr><tr><td>13</td><td>CON+GF+COMP</td><td>β = 1</td><td>69.6</td><td>106</td><td>555.9</td><td>12</td><td>1</td><td>0</td></tr><tr><td>13</td><td>CON+GF+COMP</td><td>P<sub>max</sub> = 106</td><td>35.9</td><td>106</td><td>555.9</td><td>12</td><td>1</td><td>0</td></tr><tr><td>14</td><td>CON+GF+COMP</td><td>β = 1</td><td>65.9</td><td>102</td><td>548.2</td><td>12</td><td>1</td><td>0</td></tr><tr><td>14</td><td>CON+GF+COMP</td><td>P<sub>max</sub> = 102</td><td>108.6</td><td>102</td><td>548.2</td><td>12</td><td>1</td><td>0</td></tr><tr><td>15</td><td>CON+GF+COMP</td><td>β = 1</td><td>76.5</td><td>136</td><td>582.4</td><td>18</td><td>1</td><td>0</td></tr><tr><td>15</td><td>CON+GF+COMP</td><td>P<sub>max</sub> = 136</td><td>171.6</td><td>136</td><td>582.4</td><td>18</td><td>1</td><td>0</td></tr><tr><td>16</td><td>CON+GF+COMP</td><td>β = 1</td><td>229.6</td><td>120</td><td>573.8</td><td>17</td><td>1</td><td>0</td></tr><tr><td>16</td><td>CON+GF+COMP</td><td>P<sub>max</sub> = 120</td><td>98.5</td><td>120</td><td>573.8</td><td>17</td><td>1</td><td>0</td></tr><tr><td>17</td><td>CON+GF+COMP</td><td>β = 1</td><td>78.7</td><td>94</td><td>499.4</td><td>12</td><td>1</td><td>0</td></tr><tr><td>17</td><td>CON+GF+COMP</td><td>P<sub>max</sub> = 94</td><td>36.3</td><td>94</td><td>499.4</td><td>12</td><td>1</td><td>0</td></tr></tbody></table>

<table><caption>**Table 5**
Assessment of the computation time needed to enforce compactness whether using $\beta$ or $P_{max}$ constraints for 10 instances of 25 × 20 planning units and 3 conservation features. A summary of the characteristics of the reserve solutions is provided: computation time in seconds, reserve perimeter, total cost, radius of the reserve graph, number of connected components, number of gaps. Conservation features targets are all set to 50%. If the time limit is reached (TL=1000 s), the optimality gap of the incumbent solution is given within brackets.</caption><thead><tr><th>Instance</th><th>Model</th><th>Parameters</th><th>Time</th><th>Perimeter</th><th>Cost</th><th>Radius</th><th>Components</th><th>Gaps</th></tr></thead><tbody><tr><td>4</td><td>CON+GF+COMP</td><td>β = 1</td><td>139.0</td><td>214</td><td>961.8</td><td>21</td><td>1</td><td>0</td></tr><tr><td>4</td><td>CON+GF+COMP</td><td>P<sub>max</sub> = 214</td><td>438.9</td><td>214</td><td>961.8</td><td>21</td><td>1</td><td>0</td></tr><tr><td>5</td><td>CON+GF+COMP</td><td>β = 1</td><td>TL[0.5%]</td><td>222</td><td>792.8</td><td>21</td><td>1</td><td>0</td></tr><tr><td>5</td><td>CON+GF+COMP</td><td>P<sub>max</sub> = 222</td><td>TL[1.0%]</td><td>222</td><td>794.1</td><td>21</td><td>1</td><td>0</td></tr><tr><td>6</td><td>CON+GF+COMP</td><td>β = 1</td><td>577.4</td><td>224</td><td>910.8</td><td>23</td><td>1</td><td>0</td></tr><tr><td>6</td><td>CON+GF+COMP</td><td>P<sub>max</sub> = 224</td><td>630.1</td><td>224</td><td>910.8</td><td>23</td><td>1</td><td>0</td></tr><tr><td>8</td><td>CON+GF+COMP</td><td>β = 1</td><td>50.8</td><td>164</td><td>915.9</td><td>21</td><td>1</td><td>0</td></tr><tr><td>8</td><td>CON+GF+COMP</td><td>P<sub>max</sub> = 164</td><td>109.3</td><td>164</td><td>915.9</td><td>21</td><td>1</td><td>0</td></tr><tr><td>9</td><td>CON+GF+COMP</td><td>β = 1</td><td>542.6</td><td>192</td><td>999.5</td><td>21</td><td>1</td><td>0</td></tr><tr><td>9</td><td>CON+GF+COMP</td><td>P<sub>max</sub> = 192</td><td>717.2</td><td>192</td><td>999.5</td><td>21</td><td>1</td><td>0</td></tr><tr><td>10</td><td>CON+GF+COMP</td><td>β = 1</td><td>97.6</td><td>184</td><td>927.8</td><td>18</td><td>1</td><td>0</td></tr><tr><td>10</td><td>CON+GF+COMP</td><td>P<sub>max</sub> = 184</td><td>305.8</td><td>184</td><td>927.8</td><td>18</td><td>1</td><td>0</td></tr><tr><td>18</td><td>CON+GF+COMP</td><td>β = 1</td><td>TL[0.4%]</td><td>178</td><td>976.6</td><td>23</td><td>1</td><td>0</td></tr><tr><td>18</td><td>CON+GF+COMP</td><td>P<sub>max</sub> = 178</td><td>TL[0.1%]</td><td>178</td><td>974.8</td><td>23</td><td>1</td><td>0</td></tr><tr><td>19</td><td>CON+GF+COMP</td><td>β = 1</td><td>109.2</td><td>202</td><td>956.2</td><td>23</td><td>1</td><td>0</td></tr><tr><td>19</td><td>CON+GF+COMP</td><td>P<sub>max</sub> = 202</td><td>130.7</td><td>202</td><td>956.2</td><td>23</td><td>1</td><td>0</td></tr><tr><td>20</td><td>CON+GF+COMP</td><td>β = 1</td><td>895.1</td><td>184</td><td>969.5</td><td>20</td><td>1</td><td>0</td></tr><tr><td>20</td><td>CON+GF+COMP</td><td>P<sub>max</sub> = 184</td><td>930.6</td><td>184</td><td>969.5</td><td>20</td><td>1</td><td>0</td></tr><tr><td>21</td><td>CON+GF+COMP</td><td>β = 1</td><td>388.3</td><td>166</td><td>905.8</td><td>18</td><td>1</td><td>0</td></tr><tr><td>21</td><td>CON+GF+COMP</td><td>P<sub>max</sub> = 166</td><td>727.4</td><td>166</td><td>905.8</td><td>18</td><td>1</td><td>0</td></tr></tbody></table>
[8] IUCN, IUCN world parks congress 2014 bulletin, in: International Institute for Sustainable Development, vol. 89, 2014, p. 43.

[9] IUCN, Proceedings of the Members' Assembly, IUCN, Gland, Switzerland, 2016, p. 124, www.iucn.org/resources/publications.

[10] E. Commission, EU Biodiversity Strategy for 2030, Bringing nature into our lives, COM(2020) 380 finals, Brussels, 20.05.2020, Tech. Rep., 2020.

[11] J. Baillie, Y.-P. Zhang, Space for nature, Science 361 (6407) (2018) 1051, https://doi.org/10.1126/science.aau1397.

[12] E. Dinerstein, C. Vynne, E. Sala, A.R. Joshi, S. Fernando, T.E. Lovejoy, J. Mayorga, D. Olson, G.P. Asner, J.E.M. Baillie, N.D. Burgess, K. Burkart, R.F. Noss, Y.P. Zhang, A. Baccini, T. Birch, N. Hahn, L.N. Joppa, E. Wikramanayake, A global deal for nature: guiding principles, milestones, and targets, Sci. Adv. 5 (4) (2019) eaaw2869, https://doi.org/10.1126/sciadv.aaw2869.

[13] C. Margules, R.L. Pressey, Systematic conservation planning, Nature 405 (2000) 243–253.

[14] A. Moilanen, K.A. Wilson, H.P. Possingham (Eds.), Spatial Conservation Prioritization: Quantitative Methods and Computational Tools, Oxford Biology, Oxford University Press, Oxford; New York, 2009, oCLC: ocn277205285.

[15] A. Ando, J. Camm, S. Polasky, A. Solow, Species distributions, land values, and efficient conservation, Science 279 (5359) (1998) 2126–2128.

[16] R.R. Stewart, H.P. Possingham, Efficiency, costs and trade-offs in marine reserve system design, Environ. Model. Assess. 10 (3) (2005) 203–213, https://doi.org/10.1007/s10666-005-9001-y.

[17] I.R. Ball, H.P. Possingham, Marine Reserve Design using Spatially Explicit Annealing, Tech. Rep., The Great Barrier Reef Marine Park Authority, Mar. 2000.

[18] I.R. Ball, H.P. Possingham, M.E. Watts, Marxan and Relatives: Software for Spatial Conservation Prioritisation, Spatial Conservation Prioritisation: Quantitative Methods and Computational Tools, vol. 14, Oxford University Press, 2009, pp. 185–196, https://hdl.handle.net/1959.11/20240.

[19] J.O. Hanson, R. Schuster, N. Morrel, M. Strimas-Mackey, M.E. Watts, P. Arcese, J. Bennett, H.P. Possingham, prioritizr: Systematic Conservation Prioritization in R, https://CRAN.R-project.org/package=prioritizr, 2020.

[20] J. Flower, R. Ramdeen, A. Estep, L.R. Thomas, S. Francis, G. Goldberg, A.E. Johnson, W. McClintock, S.R. Mendes, K. Mengerink, M. O'Garro, L. Rogers, U. Zischka, S.E. Lester, Marine spatial planning on the Caribbean island of Montserrat: lessons for data-limited small islands, Conserv. Sci. Pract. 2 (4) (Apr. 2020), https://doi.org/10.1111/csp.2.158.

[21] L. Fernandes, J. Day, A. Lewis, S. Slegers, B. Kerrigan, D. Breen, D. Cameron, B. Jago, J. Hall, D. Lowe, J. Innes, J. Tanzer, V. Chadwick, L. Thompson, K. Gorman, M. Simmons, B. Barnett, K. Sampson, G. De'Ath, B. Mapstone, H. Marsh, H. Possingham, I. Ball, T. Ward, K. Dobbs, J. Aumend, D. Slater, K. Stapleton, Establishing representative no-take areas in the Great Barrier Reef: large-scale implementation of theory on marine protected areas, Conserv. Biol. 19 (6) (2005) 1733–1744, https://doi.org/10.1111/j.1523-1739.2005.00302.x.

[22] R.H. MacArthur, E.O. Wilson, The Theory of Island Biogeography, Princeton University Press, Princeton, 1967.

[23] J.M. Diamond, The island dilemma: lessons of modern biogeographic studies for the design of natural reserves, Biol. Conserv. 7 (2) (1975) 129–146, https://doi.org/10.1016/0006-3207(75)90052-X.

[24] A. Billionnet, Designing an optimal connected nature reserve, Appl. Math. Model. 36 (5) (2012) 2213–2223, https://doi.org/10.1016/j.apm.2011.08.002.

[25] M.E. Watts, I.R. Ball, R.S. Stewart, C.J. Klein, K. Wilson, C. Steinback, R. Lourival, L. Kircher, H.P. Possingham, Marxan with zones: software for optimal conservation based land- and sea-use zoning, Environ. Model. Softw. 24 (12) (2009) 1513–1521, https://doi.org/10.1016/j.envsoft.2009.06.005.

[26] A. Billionnet, Mathematical optimization ideas for biodiversity conservation, Eur. J. Oper. Res. (2013) 21.

[27] H.L. Beyer, Y. Dujardin, M.E. Watts, H.P. Possingham, Solving conservation planning problems with integer linear programming, Ecol. Model. 328 (2016) 14–22, https://doi.org/10.1016/j.ecolmodel.2016.02.005.

[28] M.D. McDonnell, H.P. Possingham, I.R. Ball, E.A. Cousins, Mathematical methods for spatially cohesive reserve design, Environ. Model. Assess. 7 (2) (2002) 107–114, https://doi.org/10.1023/A:1015649716111.

[29] D.J. Nalle, J.L. Arthur, J. Sessions, Designing compact and contiguous reserve networks with a hybrid heuristic algorithm, For. Sci. 48 (1) (2002) 10, https://doi.org/10.1093/forestscience/48.1.59.

[30] J.C. Williams, C.S. ReVelle, S.A. Levin, Using mathematical optimization models to design nature reserves, Front. Ecol. Environ. 2 (4) (2004) 98–105, https://doi.org/10.1890/1540-9295(2004)002[0098:UMOMTD]2.0.CO;2.

[31] J.C. Williams, C.S. ReVelle, S.A. Levin, Spatial attributes and reserve design models: a review, Environ. Model. Assess. 10 (2005) 19.

[32] A. Billionnet, Designing connected and compact nature reserves, Environ. Model. Assess. 21 (2) (2016) 211–219, https://doi.org/10.1007/s10666-015-9465-3.

[33] A. Billionnet, Designing Protected Area Networks: A Mathematical Programming Approach, edp sciences edition, Current Natural Sciences, Edp Sciences, 2021, oCLC: 9463862519.

[34] J.C. Williams, C.S. ReVelle, Reserve assemblage of critical areas: a zero-one programming approach, Eur. J. Oper. Res. 104 (3) (1998) 497–509, https://doi.org/10.1016/S0377-2217(97)00017-9.

[35] M.A. Clemens, C.S. ReVelle, J.C. Williams, Reserve design for species preservation, Eur. J. Oper. Res. (1999) 11.

[36] E. Álvarez Miranda, M. Goycoolea, I. Ljubić, M. Sinnl, The generalized reserve set covering problem with connectivity and buffer requirements, Eur. J. Oper. Res. 289 (3) (2021) 1013–1029, https://doi.org/10.1016/j.ejor.2019.07.017.

[37] H. Önal, R.A. Briers, Incorporating Spatial Criteria in Optimum Reserve Network Selection, Royal Society, Oct. 2002.

[38] J.C. Williams, Optimal reserve site selection with distance requirements, Comput. Oper. Res. 35 (2) (2008) 488–498, https://doi.org/10.1016/j.cor.2006.03.012.

[39] H. Önal, R.A. Briers, Optimal selection of a connected reserve network, Oper. Res. 54 (2) (2006) 379–388, http://www.jstor.org/stable/25146973.

[40] Y.-C. Wang, H. Önal, Designing connected nature reserve networks using a graph theory approach, Acta Ecol. Sin. 31 (5) (2011) 235–240, https://doi.org/10.1016/j.chnaes.2011.06.001.

[41] Y. Wang, H. Önal, Designing a connected nature reserve using a network flow theory approach, Acta Ecol. Sin. 33 (5) (2013) 253–259, https://doi.org/10.1016/j.chnaes.2013.07.004.

[42] N. Jafari, J. Hearne, A new method to solve the fully connected reserve network design problem, Eur. J. Oper. Res. 231 (1) (2013) 202–209, https://doi.org/10.1016/j.ejor.2013.05.015.

[43] T. Shirabe, A model of contiguity for spatial unit allocation, Geogr. Anal. 37 (1) (2005) 2–16, https://doi.org/10.1111/j.1538-4632.2005.00605.x.

[44] J.C. Williams, Convex land acquisition with zero-one programming, Environ. Plan. B, Plan. Des. 30 (2) (2003) 255–270, https://doi.org/10.1068/b12925.

[45] A. Billionnet, Optimisation Discrète, de la modélisation à la résolution par des logiciels de programmation mathématique, dunod Edition, Dunod, https://hal.archives-ouvertures.fr/hal-01125300, 2007.

[46] M.R. Garey, D.S. Johnson, Computers and Intractability: A Guide to the Theory of NP-Completeness, 27th edition, A Series of Books in the Mathematical Sciences, Freeman, New York [u.a.], 1979.

[47] T. Abdelmaguid, An efficient mixed integer programming model for the minimum spanning tree problem, Mathematics 6 (10) (2018) 183, https://doi.org/10.3390/math6100183, http://www.mdpi.com/2227-7390/6/10/183.

[48] I. Dunning, J. Huchette, M. Lubin, JuMP: a modeling language for mathematical optimization, SIAM Rev. 59 (2) (2017) 295–320, https://doi.org/10.1137/15M1020575, arXiv:1508.01982, http://arxiv.org/abs/1508.01982.

[49] J. Bezanson, S. Karpinski, V.B. Shah, A. Edelman, Julia: a fast dynamic language for technical computing, arXiv:1209.5145, http://arxiv.org/abs/1209.5145, Sep. 2012.

[50] J. Bezanson, A. Edelman, S. Karpinski, V.B. Shah, Julia: a fresh approach to numerical computing, arXiv:1411.1607, http://arxiv.org/abs/1411.1607, Jul. 2015.

[51] A. Brunel, J. Davret, B. Trouillet, N. Bez, J. Salvetat, A. Gicquel, S. Lanco Bertrand, Opening the black box of decision support tools in marine spatial planning: shedding light into reserve site selection algorithms for a balanced empowerment of stakeholders, Preprint HAL, https://hal.archives-ouvertures.fr/hal-03445922, 2021.
[52] A. Bertrand, FAROFA 3 cruise, RV TUBARÃO tigre publisher: sismer, https://doi.org/10.17600/18001381, https://campagnes.flotteoceanographique.fr/campagnes/18001381/, 2019.

[53] J. Salvetat, N. Bez, J. Habasque, A. Lebourges-Dhaussy, C. Lopes, G. Roudaut, M. Simier, P. Travassos, G. Vargas, A. Bertrand, Comprehensive spatial distribution of tropical fish assemblages from multifrequency acoustics and video fulfils the island mass effect framework, https://doi.org/10.1038/s41598-022-12409-9, https://www.nature.com/articles/s41598-022-12409-9, 2022.

[54] A. Brunel, J. Omer, S. Lanco Bertrand, Producing a diverse set of near-optimal reserve solutions with exact optimisation, Environ. Model. Assess. (Dec. 2022), https://doi.org/10.1007/s10666-022-09862-1, https://link.springer.com/10.1007/s10666-022-09862-1.

[55] Y. Wang, H. Önal, Q. Fang, How large spatially-explicit optimal reserve design models can we solve now? An exploration of current models’ computational efficiency, Nat. Conserv. 27 (2018) 17–34, https://doi.org/10.3897/natureconservation.27.21642, https://natureconservation.pensoft.net/articles.php?id=21642.

[56] S.A. Malcolm, C. ReVelle, Models for preserving species diversity with backup coverage, Environ. Model. Assess. 10 (2005) 99–105, https://doi.org/10.1007/s10666-004-5101-3.

[57] H.J. Albers, G.M. Busby, B. Hamaide, A.W. Ando, S. Polasky, Spatially-Correlated Risk in Nature Reserve Site Selection, PLoS ONE (2016).

[58] B. Hamaide, H.J. Albers, G. Busby, Backup coverage models in nature reserve site selection with spatial spread risk heterogeneity, Socio-Econ. Plan. Sci. 48 (2) (2014) 158–167, https://doi.org/10.1016/j.seps.2014.03.003.