problem  Input: Question: optimisation_program min

$\textnormal{\@name}\left\\{\begin{tabular}[]{rlll}{\@sense}&&\\\
{s.t.}\end{tabular}\right.$

# Optimizing the ecological connectivity of landscapes with generalized flow
models and preprocessing

François Hamonic  Aix-Marseille Univ, CNRS, Université de Toulon, LIS,
Marseille, France  Aix-Marseille Univ, CNRS, Univ Avignon, IRD, IMBE,
Marseille, France  Cécile Albert  Aix-Marseille Univ, CNRS, Univ Avignon, IRD,
IMBE, Marseille, France  Basile Couëtoux  Aix-Marseille Univ, CNRS, Université
de Toulon, LIS, Marseille, France  Yann Vaxès  Aix-Marseille Univ, CNRS,
Université de Toulon, LIS, Marseille, France

###### Abstract

In this article we consider the problem of optimizing the landscape
connectivity under a budget constraint by improving habitat areas and
ecological corridors between them. We model this problem as a discrete
optimization problem over graphs in which vertices represent the habitat areas
and arcs represent a probability of connection between two areas that depend
on the quality of the respective corridor. We propose a new generalized flow
model that improves existing models for this problem. Then, following the
approach of Catanzaro et al. [7] for the robust shortest path problem, we
design an improved preprocessing algorithm that reduces the size of the graphs
on which we compute generalized flows. Reported numerical experiments
highlight the benefits of both contributions that allow to solve larger
instances of the problem. These experiments also show that several variants of
greedy algorithms perform relatively well in practice while they return
arbitrary bad solutions in the worst case.

Keywords : Combinatorial optimization, environment and climate change,
landscape connectivity, network flow, shortest path computation, mixed integer
linear program.

##  1 Introduction

###  1.1 Context and ecological motivation

Habitat loss is a major cause of the rapid decline of biodiversity [6].
Further than reducing the available resources, it also increases the
discontinuities among small habitat areas called patches, this is a phenomenon
known under the name habitat fragmentation [12]. While habitat loss tends to
reduce the size of populations of animals and plants, habitat fragmentation
also makes it harder for organisms to move around in landscapes. This
decreases the access to resources and the gene flow among populations.
Landscape connectivity, defined as the degree to which the landscape
facilitates the movement of organisms between habitat patches [26], then
becomes of major importance for biodiversity and its conservation. Accounting
for landscape connectivity in restoration or conservation plans thus appears
as a key solution to maximize the return on investment of the scarce financial
support devoted to biodiversity conservation.

Graph-theoretical approaches are useful in modelling habitat connectivity
[27]. Indeed, a landscape can be viewed as a directed graph in which vertices
are the habitat patches and each arc indicates a way for individuals to travel
from one patch to another. The weight of a vertex represents its quality –
patch area is often used as a surrogate for quality – and the weight of an arc
mesures the difficulty for an individual to make the corresponding travel.
This quantity is often approximated by a function of the border to border
distance between the patches. Interestingly, this approach can be used for a
variety of ecological systems like terrestrial (patches of forests in an
agricultural area, networks of lakes or wetlands), riverine (patches are
segments of river than can be separated by human constructions like dams that
prevent fishes’ movement) or marine (patches can be reefs that are connected
by flows of larvae transported by currents). With this formalism, ecologists
have developed many connectivity indicators [16, 23, 14] that aim to quantify
the quality of a landscape with respect to the connections between its habitat
patches.

###  1.2 Indicators of landscape connectivity

Among the proposed indicators, the Probability of Connectivity (PC) [23] and
its derivative the Equivalent Connected Area (ECA) [25] have received
encouraging empirical support [18, 3, 24]. Given a graph $G=(V,A)$ with
probability on edges $(\pi_{a})_{a\in A}$ and weights on vertices
$(w_{v})_{v\in V}$,

$ECA(G)=\sqrt{\sum_{s\in V}\sum_{t\in V}\left(w_{s}w_{t}\Pi_{st}\right)}$

where $\Pi_{st}$ is the probability of connection from patch $s$ to patch $t$

$\Pi_{st}=\max_{P:\textnormal{$st$-path}}\prod_{a\in P}\pi_{a}\,.$

The following probabilistic analysis explains why this indicator has been
called Equivalent Connected Area in [25]. Let $\cal W$ be the area of a
rectangle containing the landscape under study. We consider a stochastic
process that consists in choosing two points $p$ and $q$ uniformly at random
in the rectangle. The indicator PC is the expected value of a random variable
equal to $0$ if either $p$ or $q$ does not belong to a patch and $\Pi_{st}$ if
$p$ belongs to $s$ and $q$ belongs to $t$ (recall that $\Pi_{st}=1$ if $s=t$).
Let $w_{u}$ denote the area of the patch $u.$ Since the probability that $p$
belongs to $u$ is $w_{u}/{\cal W}$ and the events $p\in s$ and $q\in t$ are
independent, by linearity of expectation, $PC$ can be expressed as follows:

$PC(G)=\frac{\sum_{s,t\in V}w_{s}w_{t}\Pi_{st}}{{\cal
W}^{2}}=\frac{ECA(G)^{2}}{{\cal W}^{2}}\,.$

(a) ecological landscape [21]

(b) graph representation

Figure 1: Modelisation of an ecological landscape by a graph.

The equivalent connected area of a landscape is the area of a single patch
whose PC value is equal the PC value of the original landscape. If the area of
the patches and the landscape are normalized to make $\cal W$ equal to $1$
then PC is the square of ECA. Therefore, optimizing PC and optimizing ECA are
equivalent problems. ECA is often considered by researchers interested in
landscape connectivity because it represents an area, a more concrete quantity
than the expected value of a random variable.

###  1.3 Optimizing ECA under a budget constraint

One key question for which landscape connectivity indicators have been used in
the last years is to identify which elements of the landscape (habitat
patches, corridors) should be preserved from destruction or restored in order
to maintain a well-connected network of habitat under a given budget
constraint [2]. This translates into identifying the set of vertices or arcs
that optimally maintains a good level of connectivity. Many mathematical
programming models have been introduced in the literature to help decision-
makers to protect biodiversity. For a review of these models, we refer to the
monography [4], the survey article [5] and the references therein. Here, we
consider the budget-constrained ECA optimization problem (BC-ECA-Opt) that
looks for the best combination of arcs among a set $\Phi$ that could be
restored or upgraded in order to maximize the ECA in the landscape within a
limited budget:

Input: a graph $G=(V,A)$ with weights on vertices $(w_{v})_{v\in V}$ and
probabilities on the arcs $(\pi_{a})_{a\in A}$, $\Phi\subseteq A$ a subset of
arcs with improved probabilities $(\pi^{\prime}_{a})_{a\in\Phi}$ and costs
$(c_{a})_{a\in\Phi}$, and a budget $B\in\mathbb{N}$.

Ouput: A subset $S\subseteq\Phi$ such that $\sum_{a\in S}c_{a}\leq B$
maximizing $ECA(G(S))$ where $G(S)$ is the graph $G$ where we replace the
value of $\pi_{a}$ by $\pi^{\prime}_{a}$ for all $a\in S$.

Interestingly, this problem addresses both restoration and conservation cases.
The restoration case starts from the current landscape and aims to restore a
set of elements among the different feasible options. The conservation case
starts with a landscape in which we search for the elements to protect among
those that will be altered or destroyed if nothing is done.

###  1.4 Previous works

Until now, ecologists have mostly tackled this problem by ranking each
conservation or restoration option by its independent contribution to ECA,
i.e. the amount by which ECA varies if the option is purchased alone. Such an
approach overlooks the cumulative effects of the decisions made like
unnecessary redundancies or potential synergistic effects, e.g. improving two
consecutive corridors results in a greater increase in ECA than the sum of the
increases achieved by improving each corridor independently. This could lead
to solutions that are more expensive or less beneficial to ECA than an optimal
solution. Some studies have tried to overcome this limitation by considering
tuples of options [20, 17]. In [20], the authors show that the brute force
approach rapidly becomes impractical for landscape with more than 20 patches
of habitat. Few studies have explored the search for an optimal solution but,
in most cases, the underlying graph was acyclic (river dendritic networks).
For instance, a polynomial time approximation scheme has been proposed when
the underlying graph is a tree. Indeed, [28] describes a dynamical programming
algorithm with rounding that computes, for any $\epsilon,$ a
$(1-\epsilon)$-approximated solution in time $n^{8}/\epsilon$ where $n$ is the
number of nodes of the tree. More recently, [29] has introduced a XOR sampling
method based on a mixed integer formulation. To our knowledge, this is the
unique linear programming formulation of BC-ECA-Opt already proposed. Solving
optimally the problem with their mixed integer formulation does not scale to
landscapes with few hundreds patches.

To compute $ECA(G)$ we have to compute $\Pi_{st}$ for every pair of vertices
$s,t\in V.$ If we consider the length function $l$ defined by
$l_{a}=-\log\pi_{a}$ for $a\in A$ then the distance function $d$ induced by
$l$ on $G$ verifies $d(s,t)=-\log(\Pi_{st})$ for $s,t\in V$. Thus our problem
is closely related to solving shortest path problems on a graph with different
length functions, namely one length function for each restoration plan (or
scenario). A preprocessing step has been proposed by Catanzaro et al. [7] to
address such problems. It consists in identifying a subset of arcs that can
either be removed or contracted to reduce the size of the graph considered.
More formally, following [7], we call an arc $(u,v)$ $t$-strong if it belongs
to a shortest path from $u$ to $t$ in every possible scenario. Symetrically,
we call an arc $(u,v)$ $t$-useless if there is no scenario such that $(u,v)$
is on a shortest path from $u$ to $t.$ Useless arcs were called $0$-persistent
in [7]. Since we need to specify a target vertex $t$ for which $(u,v)$ is
useless, we didn’t adopt the same terminology. In [7], given an arc $(u,v)$
and a vertex $t$, the authors give a sufficient (but not necessary) condition
for $(u,v)$ to be $t$-strong and use it to design a $O(|A|+|V|\log|V|)$
algorithm that decides, in most cases, whether $(u,v)$ is $t$-strong. The
authors also propose a $O(|A|+|V|\log|V|)$ algorithm to test whether a given
arc $(u,v)$ is $t$-useless.

###  1.5 Contribution

The contribution of this article is twofold. First, we propose a new mixed
integer linear formulation which is more compact than the one of [29]. Then,
we design a new preprocessing algorithm that computes $t$-strong and
$t$-useless arcs for all $t\in V.$

Our new mixed integer formulation is based on a generalized flow formulation
(see for instance [1]) instead of a standard network flow formulation. This
leads to two improvements. Firstly, our formulation has a linear objective
function as opposed to the model of [29] that was using a piecewise constant
approximation and additional binary variables to handle non-linearity.
Secondly, the new formulation aggregates into a single generalized flow the
contribution to the connectivity of several source/sink pairs having the same
source whereas the previous model treated every pair separately. More
precisely, our model uses $O(|V||A|)$ flow variables with $O(|V|^{2})$
constraints, and $|\Phi|$ integer variables, while the model of [29] was using
of $O(|V|^{2}|A|)$ flow variables with $O(|V|^{3})$ constraints, and
$O(|\Phi|+|V|^{2})$ integer variables. Recall that $\Phi$ is the set of
elements that are threatened or could be restored.

The other contribution of this article is a preprocessing step that speeds-up
the resolution of the problem by reducing the size of the directed graph on
which a generalized flow to a particular sink $t$ have to be computed. Our
algorithm improves the one of [7] in two ways. First, we replace the
sufficient condition defined in [7] by a necessary and sufficient condition.
This allows us to compute the entire set of $t$-strong arcs instead of a
subset of them. We also improve the time complexity by a factor $|V|$ as we
only need to run our algorithm on each arc instead of each pair of arc and
vertex.

###  1.6 Organization of the article

Section 2 is devoted to the description of our mixed integer formulation of
the problem. In Section 3 we first explain why the removal of $t$-useless arcs
and the contraction of $t$-strong arcs does not modify the objective function
of any restoration/conservation plan. Then, we design and analyze a
$O(|A|+|V|\log|V|)$ algorithm that computes the set of all vertices $t$ for
which a given arc $(u,v)$ is $t$-strong or the set of all vertices $t$ for
which $(u,v)$ is $t$-useless. In section 4, we describe several greedy
algorithms for BC-ECA-Opt and provide some instances where they compute
solutions far from being optimal. Section 5 compares our optimization approach
to simple greedy algorithms in terms of running times and quality of the
solutions found on a set of experimental cases. These experiments show that
the preprocessing is very effective and that the greedy approach performs
quite well on these instances.

##  2 An improved MIP formulation for BC-ECA-Opt

We decompose $ECA(G)$ as $\sum_{t\in V}w_{t}f_{t}$ where $f^{t}=\sum_{s\in
V}w_{s}\cdot\Pi_{st}$. We first show that $f^{t}$ can be expressed as the
maximum quantity of generalized flow that can be sent to $t$ across the
network if, for every patch $s\in V$, $w_{s}$ units of flow are available at
$s$ and appropriate multipliers are chosen for each arc of the network. Recall
that a generalized flow differs from a standard flow by the fact that each arc
$a$ has a multiplier $\pi_{a}$ such that the quantity of flow leaving arc $a$
is equal to the quantity of flow $\phi_{a}$ entering in $a$ multiplied by
$\pi_{a}$ (see [1] for an introduction to network flows). For each vertex
$u\in V,$ let $\delta^{\text{out}}_{u}$ be the set of arcs leaving $u$ and
$\delta^{\text{in}}_{u}$ the set of arcs entering $u$. As explained below, the
linear program $\mathcal{(P)}$ has $f^{t}$ as optimum value.

max

Constraints (A1) require that the total quantity of flow leaving $u$ is at
most $w_{u}$ plus the total quantify of flow entering $u$, i.e. the quantity
of flow available in vertex $u$. Constraint (A2) requires that $z$ is equal to
the total quantity of flow entering $t$ plus $w_{t}$ minus the total quantity
of flow leaving $t$. Finally, constraints (A3) state that each arc carries a
non negative quantity of flow from its source to its sink.

######  Lemma 1.

Any optimal solution of $\mathcal{(P)}$ is obtained by sending, for every
vertex $s\in V-\\{t\\},$ $w_{s}$ units of flow from $s$ to $t$ along most
reliable $st$-paths.

###### Proof.

First notice that the quantity of flow arriving at $t$ when $w_{s}$ units flow
are sent from $s$ to $t$ along an $st$-path $P$ is $w_{s}$ times the
probability of path $P.$ Let $\phi^{\prime}$ be an optimal solution of
$\mathcal{(P)}$ maximizing the quantity of flow routed along a path which is
not a most reliable path. Suppose, by contradiction, that $\phi^{\prime}$
sends $\epsilon>0$ units of flow along an $st$-path $P^{\prime}$ whose
probability is smaller than the probability of a most reliable $st$-path $P.$
Let $\phi$ be the flow obtained from $\phi^{\prime}$ by decreasing the flow
sent on path $P^{\prime}$ by $\epsilon$ and increasing the flow sent on path
$P$ by $\epsilon.$ Since the probability of $P$ is larger than the probability
of $P^{\prime},$ $\phi$ sends more flow to $t$ than $\phi^{\prime},$ leading
to a contradiction with the choice of $\phi^{\prime}.$ ∎

######  Corollary 1.

The optimal value of $\mathcal{(P)}$ is $f^{t}=\sum_{s\in
V}w_{s}\cdot\Pi_{st}$.

###### Proof.

Since the objective is to maximize $z,$ no flow leaves $t$ in any optimal
solution of $\mathcal{(P)},$ i.e.
$\sum\limits_{a\in\delta^{\text{out}}_{t}}\phi_{a}=0.$ Constraint (A2) ensures
that $z$ is the quantity of flow received by $t$ plus $w_{t}.$ By Lemma 1,
there exits an optimal solution $\phi$ such that every vertex $s$ distinct
from $t$ send $w_{s}$ units of flow on a most reliable $st$-path. Hence, for
every vertex $s$ distinct from $t$, the flow received by $t$ from $s$ is
$w_{s}\cdot\Pi_{st}$ and thus the value of $z$ is $\sum_{s\in
V}w_{s}\cdot\Pi_{st}.$ ∎

With this linear programming definition of $f^{t},$ we can now present our MIP
program of BC-ECA-Opt. For each arc $a=(u,v)\in\Phi,$ we add another arc
$a^{\prime}=(u,v)$ of probability $\pi_{a^{\prime}}=\pi^{\prime}_{a}$ that can
be viewed as an improved copy of arc $a$. This improved copy can be used only
if the improvement of arc $a$ is purchased. For each arc $a\in\Phi,$ $x_{a}$
is equal to one if the improvement of arc $a$ is purchased and zero otherwise.
We denote by $\Psi$ the set of improved copies of arcs in $\Phi.$ In the
following MIP program, $\delta^{\text{out}}_{u}$ and $\delta^{\text{in}}_{u}$
are defined with respect to the set of arcs $A^{\prime}=A\cup\Psi.$

max

Constraints of (B1) and (B2) are simply constraints of (A1) and (A2) for all
possible target vertex $t$. For each arc $a\in\Phi,$ the big-M constraints
(B3) ensure that if the improvement of arc $a$ is not purchased, i.e.
$x_{a}=0,$ then the flow on arc $a^{\prime}$ is null. The constant $M_{a}$ is
an upper bound of the flow value on arc $a.$ We could simply take
$M_{a}=\sum_{u\in V}w_{u}$ for all $a\in\Phi$ but more precise estimations are
possible for a better linear relaxation and thus a faster resolution of the
MIP program. Constraint (B4) ensures that the total cost of the improvements
is at most $B.$ This MIP program has $O(|V|(|A|+|\Phi|))$ flow variables,
$|\Phi|$ binary variables and $O(|V|^{2}+|V||\Phi|)$ constraints.

###  2.1 Extension to patch improvements

Here, we explain how to extend our model to the version with patch
improvements where it is possible to increase the weight of a vertex $u$ from
$w_{u}$ to $w^{+}_{u}$ at cost $c_{u}.$ For that, we process all the
occurrences of $w_{u}$ as follows. Let $y_{u}$ be a binary variable equal to
$1$ if the vertex $u$ is improved and $0$ otherwise. We add a term
$y_{u}c_{u}$ in the budget constraint for every vertex $u$ in the set $W$ of
vertices that can be improved. When $w_{u}$ appears as an additive constant,
we simply replace it by $w_{u}+y_{u}(w^{+}_{u}-w_{u})$. Note that $w_{u}$ only
appears as a coefficient in the objective function in the form $w_{u}f_{t}$.
In this case, to avoid a quadratic term, we use a standard McCormick
linearization [13]. We replace the product $w_{u}f_{t}$ by
$w_{u}f_{t}+(w^{+}_{u}-w_{u})f_{t}^{\prime}$ where $f_{t}^{\prime}$ is a new
variable that is equal to $f_{t}$ if $y_{u}=1$ and $0$ otherwise. To achieve
this values of $f_{t}^{\prime}$, for all $u\in V,$ we add the constraints
$f_{t}^{\prime}\leq f_{t}$ and $f_{t}^{\prime}\leq y_{u}M$ where $M$ is larger
than any values of $f_{t}$. As it is a maximization program and
$f_{t}^{\prime}$ appears with a positive coefficient in the objective
function, the first constraint guaranties that $f_{t}^{\prime}=f_{t}$ if
$y_{u}=1$ in any optimal solution and the second constraint guaranties that
$f_{t}^{\prime}=0$ if $y_{u}=0$.

##  3 Preprocessing

The size of the mixed integer programming formulation of BC-ECA-Opt given in
Section 2 grows quadratically with the size of the graph that represents the
landscape. In this section, we describe a preprocessing step that reduces the
size of this graph. For that, we adopt the approach used by Catanzaro et al.
[7] for the robust shortest path problem. We introduce a notion of strongness
of an arc with respect to a target vertex $t\in V$. We call a solution
$x\in\\{0,1\\}^{\Phi}$ of BC-ECA-Opt a scenario and say that the distances are
computed under scenario $x$ when the length of every $a\in A$ is $l^{-}_{a}$
if $x_{a}=1$ and $l_{a}$ otherwise. We denote $d_{x}(s,t)$ the distance
between $s$ and $t$ when the arc lengths are set according to the scenario
$x.$ Following [7], an arc $(u,v)$ is said to be $t$-strong if, for every
scenario $x\in\\{0,1\\}^{\Phi}$, $(u,v)$ belongs to a shortest path from $u$
to $t,$ i.e. $d_{x}(u,t)=d_{x}(u,v)+d_{x}(v,t)$ for every
$x\in\\{0,1\\}^{\Phi}$. An arc $(u,v)$ is said to be $t$-useless if, for every
scenario $x\in\\{0,1\\}^{\Phi}$, arc $(u,v)$ does not belong to any shortest
$ut$-path when arc lengths are set according to the scenario $x.$ Useless arcs
were called $0$-persistent in [7] but here we need to specify the target $t$.
We denote by $S(t)$ the set of arcs $(u,v)$ such that $(u,v)$ is $t$-strong
and by $W(t)$ the set of arcs $(u,v)$ such that $(u,v)$ is $t$-useless. In
[7], given an arc $(u,v)$ and a vertex $t$, the authors identify a sufficient
(but not necessary) condition for $(u,v)\in S(t)$ and use it to design a
$O(|A|+|V|\log|V|)$ algorithm that can, in most cases, identify if $(u,v)\in
S(t)$. The authors also propose a $O(|A|+|V|\log|V|)$ algorithm to test
whether a given arc $(u,v)$ belongs to $W(s)$ or not. In Section 3.1, given an
arc $(u,v)$ of $G$, we show how to adapt the Dijkstra’s algorithm to compute
in $O(|A|+|V|\log|V|)$ the set of all vertices $t$ such that $(u,v)\in S(t)$
or the set of all vertices $t$ such that $(u,v)\in W(t)$. This improves the
results of [7] by providing a necessary and sufficient condition for $(u,v)$
to be $t$-strong, i.e. we compute the entire set $S(t)$ while the algorithm of
[7] computes a subset of $S(t).$ It also reduces the time complexity for
computing $S(t)$ and $W(t)$ for all $t$ by a factor $|V|$ as we only need to
run the algorithms on each arc instead of each pair of arc and vertex. Then,
we explain how the knowledge of $S(t)$ and $W(t)$ for all $t$ can be used to
define a smaller equivalent instance of the problem.

###  3.1 Computing $S(t)$ for all $t$

Given a scenario $x\in\\{0,1\\}^{\Phi},$ the fiber $F_{x}(u,v)$ of arc
$(u,v)\in A$ is the set of vertex $t$ such that $(u,v)$ belongs to a shortest
path from $u$ to $t$ when arc lengths are set according to $x,$ i.e.

$F_{x}(u,v)=\\{t\in V:d_{x}(u,t)=l_{x}(u,v)+d_{x}(v,t)\\}$

Let $F(u,v)$ be the intersection of the fibers of $(u,v)\in A$ over all
possible scenarios $x\in\\{0,1\\}^{\Phi},$ i.e.
$F(u,v):=\bigcap_{x}F_{x}(u,v).$ By definition, an arc $(u,v)$ is $t$-strong
if $t$ belongs to $F_{x}(u,v)$ for every scenario $x$ i.e.

$S(t)=\left\\{(u,v):t\in F(u,v)\right\\}$

In order to compute every $S(t)$, we first compute $F(u,v)$ for every arc
$(u,v)\in A$ and then we transpose the representation to get $S(t)$ for every
vertex $t\in V$. Let $y\in\\{0,1\\}^{\Phi}$ be the following scenario :

$y_{wt}=\left\\{\begin{array}[]{ll}0&w\in F(u,v)\textbf{or}(w,t)=(u,v)\\\
1&w\notin F(u,v)\end{array}\right.$ (1)

######  Lemma 2.

The intersection $F(u,v)$ of the fibers of $(u,v)\in A$ over all possible
scenarios is the fiber of $(u,v)$ under the scenario $y,$ i.e.
$F(u,v)=F_{y}(u,v).$

###### Proof.

Since the inclusion $F(u,v)\subseteq F_{y}(u,v)$ is obvious, it suffices to
prove that $F_{y}(u,v)\subseteq F_{x}(u,v)$ for any scenario
$x\in\\{0,1\\}^{\Phi}.$ By way of contradiction, suppose that $x$ is a
scenario such that $F_{y}(u,v)\setminus F_{x}(u,v)$ contains a vertex $t$ at
minimum distance from $u$ in the scenario $y.$ Let $P$ be a shortest $ut$-path
containing $(u,v)$ under the scenario $y.$ For every vertex $w$ of $P,$
$d_{y}(u,w)\leq d_{y}(u,t)$ and $w\in F_{y}(u,v).$ Hence, by the choice of
$t$, $w$ belongs to $F_{x}(u,v)$ for every scenario $x.$ Therefore $w$ belongs
to $F(u,v)$ and the length of every arc of $P$ is set to its upper bound in
the scenario $y,$ i.e. $l_{y}(P)=l(P).$ Now, let $Q$ be a shortest $ut$-path
in the scenario $x.$ If the path $Q$ contains a vertex $z\in F(u,v)$ then
$d_{x}(u,t)=d_{x}(u,z)+d_{x}(z,t)=l_{x}(u,v)+d_{x}(v,z)+d_{x}(z,t)=l_{x}(u,v)+d_{x}(v,t),$
a contradiction with $t\notin F_{x}(u,v).$ Therefore, the vertices of $Q$ do
not belong to $F(u,v)$ and thus their lengths are set to their lower bound in
$y,$ i.e. $l(Q)=l^{-}(Q).$ We deduce that $l_{x}(P)\leq l_{y}(P)\leq
l_{y}(Q)\leq l_{x}(Q)$, which contradicts $t\notin F_{x}(u,v).$ ∎

We describe a $O(|A|+|V|\log|V|)$ time algorithm that, given an arc $(u,v)$,
computes simultaneously the scenario $y$ defined by (1) and the fiber
$F_{y}(u,v)$ of arc $(u,v)$ with respect to this scenario. The algorithm is an
adaptation of the Dijkstra’s shortest path algorithm that assigns colors to
vertices. We prove that it colors a vertex $w$ in blue if $w$ belongs to
$F_{y}(u,v)$ and in red otherwise. At each step, we consider a subset of
vertices $S\subseteq V$ whose colors have been already computed. Before the
first iteration, $S=\\{u\\}$ and the length of every arc leaving $u$ is set to
its lower bound except the length of $(u,v)$ which is set to its upper bound
according to scenario $y.$ At each step, since the color of every vertex of
$S$ is known, the length, under scenario $y,$ of every arc $(w,t)$ with $w\in
S$ is also known. Therefore, it is possible to find a vertex $t\in V-S$ at
minimum distance from $u,$ under scenario $y,$ in the subgraph
$G[S\cup\\{t\\}]$ induced by $S\cup\\{t\\}.$ Following Dijkstra’s algorithm
analysis, we know that the distance under scenario $y$ from $u$ to $t$ in
$G[S\cup\\{t\\}]$ is in fact the distance between $u$ and $t$ in $G$. This
allows us to determine if there exists a shortest path from $u$ to $t$ in the
scenario $y$ passing via $(u,v)$ and to color the vertex $t$ accordingly. We
end-up the iteration with a new vertex $t$ whose color is known and that can
be added to $S$ before starting the next iteration. The algorithm terminates
when $S=V.$

Input : $G=(V,A,l,l^{-})$, $(u,v)\in A$

Output : $\\{t\in V:(u,v)\text{is}t\text{-strong}\\}$

foreach  _$(u,w)\in\delta^{\text{out}}_{u}\setminus\\{(u,v)\\}$_ do

$d(w)\leftarrow l^{-}_{uw}$

$\gamma(w)\leftarrow red$

$d(v)\leftarrow l_{uv}$; $\gamma(v)\leftarrow blue$

$S\leftarrow\\{u\\}$; $\gamma(u)\leftarrow red$

while  _$S\neq V$_ do

Pick $t\in V-S$ with smallest $d(t)$ breaking tie by choosing a vertex $t$
such that $\gamma(t)=blue$ if it exists

if  _$\gamma(t)$ is $blue$_ then

foreach  _$(t,w)\in\delta^{\text{out}}_{t}$ such that $d(w)\geq d(t)+l_{tw}$_
do

$d(w)\leftarrow d(t)+l_{tw}$

$\gamma(w)\leftarrow blue$

else

foreach  _$(t,w)\in\delta^{\text{out}}_{t}$ such that $d(w)>d(t)+l^{-}_{tw}$_
do

$d(w)\leftarrow d(t)+l^{-}_{tw}$

$\gamma(w)\leftarrow red$

$S\leftarrow S\cup\\{t\\}$

return $\\{t\in V:\gamma(t)=blue\\}$

Algorithm 1 Computes the set $S$ of vertices $t$ such that $(u,v)$ is
$t$-strong

For every vertex $w\in V-S$, the estimated distance $d(w)$ is the length of a
$uw$-path under the scenario $y$ in the subgraph $G[S\cup\\{w\\}].$ An arc $a$
is blue if $a=uv$ or if its origin is blue, the other arcs are red, i.e. $a$
is blue if $y_{a}=0$ and red if $y_{a}=1.$ The estimated color $\gamma(w)$ of
$w$ is blue if there exists a blue $uw$-path of length $d(w)$ in
$G[S\cup\\{w\\}]$ and red otherwise. The correctness of Algorithm 1 follows
from the following Lemma.

######  Lemma 3.

For every vertex $t\in S,$ $\gamma(t)$ is blue if and only if $t\in
F_{y}(u,v).$

###### Proof.

We proceed by induction on the number of vertices of $S.$ When $S=\\{u\\},$
the property is verified. Now, suppose the property true before the insertion
in $S$ of the vertex $t$ such that $d(t)$ is minimum. By induction hypothesis,
the lengths of arcs having their source in $S$ are set according to $y.$
Therefore, following Dijkstra’s algorithm analysis, we deduce that $d(t)$ is
the length of the shortest path from $u$ to $t$ in the graph $G$ under
scenario $y.$ If $t$ has been colored blue then there exists a blue vertex
$w\in S$ such that $d(t)=d(w)+l_{wt}.$ By induction hypothesis $w\in
F_{y}(u,v)$ and there exists a shortest $ut$-path under scenario $y$
containing $(u,v),$ i.e. $t\in F_{y}(u,v).$ Now, suppose that $t$ has been
colored in red. By contradiction, assume there exists a shortest $ut$-path
under scenario $y$ that contains $(u,v).$ This path cannot contain a vertex
outside $S$ except $t$ because otherwise, since arc lengths are non-negative,
its length according to $y$ would be greater than $d(t)$ by the choice of $t.$
Therefore, the predecessor $w$ of $t$ in this path belongs to $S.$ Since $w$
belongs to a shortest $ut$-path passing via $(u,v),$ by induction, it was
colored blue. But in this case, there exists a blue vertex $w$ such that
$d(w)+l(w,t)=d(t),$ and $t$ was colored blue as well, a contradiction. ∎

We are now ready to state the main result of this section.

######  Proposition 1.

Given a graph $G=(V,A),$ two arc-length functions $l^{-}$ and $l$ such that
$0\leq l^{-}_{a}\leq l_{a}$ for every arc $a\in A$, and an arc $(u,v)$,
Algorithm 1 computes in $O(|A|+|V|\log|V|)$ the set of vertex $t$ such that
$(u,v)$ is $t$-strong.

###### Proof.

The lemma 2 shows that given an arc $(u,v)$ the set of vertices $S$ such that
$t\in S$ if and only if $(u,v)$ is $t$-strong is a fiber for a specific
scenario. The lemma 3 shows that the algorithm compute this fiber. Therefore
the Algorithm 1 is correct. Analogously to the Dijkstra’s algorithm, Algorithm
1 can be implemented to run in $O(|A|+|V|\log|V|)$ by using a Fibonacci heap
as priority queue. ∎

###  3.2 Computing $W(t)$ for all $t$

An adaptation of Algorithm 1 can compute, given an arc $(u,v)$ the set of
vertices $W$ such that $t\in W$ if and only if $(u,v)$ is $t$-usless. Before
describing this adaptation, we explain how the two problems are related. For
that, we first introduced a strengthening of the notion of strongness. We will
say that an arc $(u,v)$ is strictly $t$-strong if it belongs to all shortest
$ut$-paths for every scenario $x\in\\{0,1\\}^{\Phi}.$ Recall that
$t$-strongness requires only the existence of a shortest $ut$-path passing via
$(u,v)$ for every scenario $x\in\\{0,1\\}^{\Phi}.$ Algorithm 1 can be easily
adapted to compute for every arc $(u,v)$ the set of vertex $t$ such that
$(u,v)$ is strictly $t$-strong. It suffices to change the way, the algorithms
breaks tie between a red and a blue path and the choice of $t$ in case of tie.
Namely, in the first internal loop, the condition for coloring $w$ in blue
becomes $d(w)>d(t)+l_{tw}$ while the condition for coloring $w$ in red in the
second internal loop becomes $d(w)\geq d(t)+l^{-}_{tw}.$ Moreover, when we
choose $t$ such that $d(t)$ is minimal, we break tie by choosing a red vertex
if it exists. Clearly, these small changes exclude the existence of a red path
of length $d(w)$ between $u$ and a blue vertex $w$. Therefore, $(u,v)$ belongs
to every path of length $d(w)$ and $(u,v)$ is strictly $w$-strong whenever $w$
is blue. We call the resulting algorithm the strict version of Algorithm 1.
The next step is to extend the notion of strict strongness to a subset of arcs
having the same source. For any vertex $u\in V,$ a subset
$\Gamma\subseteq\Gamma_{u}^{+}$ of arcs is strictly $t$-strong if, for every
scenario $x\in\\{0,1\\}^{\Phi},$ all shortest $ut$-paths intersect $\Gamma.$
By definition, an arc $(u,v)$ is $t$-useless if and only if
$\Gamma_{u}^{+}\setminus\\{(u,v)\\}$ is strictly $t$-strong. Indeed, every
$ut$-path avoiding $(u,v)$ intersects $\Gamma_{u}^{+}\setminus\\{(u,v)\\}$
and, conversely, $(u,v)$ belongs to every $ut$-path avoiding
$\Gamma_{u}^{+}\setminus\\{(u,v)\\}$. Hence, computing the set of vertex $t$
such that $(u,v)$ is $t$-useless amounts to compute the set of vertex $t$ such
that $\Gamma_{u}^{+}\setminus\\{(u,v)\\}$ is strictly $t$-strong. An algorithm
that computes this set of vertices can be obtained from the strict version of
Algorithm 1 by modifying only the initialization step: arc $(u,v)$ is colored
in red and its length is set to $l^{-}_{uv}$ while arcs of
$\delta^{\text{out}}_{u}\setminus\\{(u,v)\\}$ are colored in blue and their
lengths are set to their upper bound. A correctness proof very similar to the
one of Algorithm 1 (and that we will not repeat) shows that a vertex is
colored blue by Algorithm 2 if and only if $(u,v)$ is $t$-useless. Since the
two algorithms have clearly the same time complexity, we conclude this section
with the following result.

Input : $G=(V,A,l,l^{-})$, $(u,v)\in A$

Output : $\\{t\in V:(u,v)\text{is}t\text{-useless}\\}$

$d(v)\leftarrow l^{-}_{uv}$; $\gamma(v)\leftarrow red$

foreach  _$(u,w)\in\delta^{\text{out}}_{u}\setminus\\{(u,v)\\}$_ do

$d(w)\leftarrow l_{uw}$

$\gamma(w)\leftarrow blue$

$S\leftarrow\\{u\\}$; $\gamma(u)\leftarrow blue$

while  _$S\neq V$_ do

Pick $t\in V-S$ with smallest $d(t)$ breaking tie by choosing a vertex $t$
such that $\gamma(t)=red$ if it exists

if  _$\gamma(t)$ is $blue$_ then

foreach  _$(t,w)\in\delta^{\text{out}}_{t}$ such that $d(w)>d(t)+l_{tw}$_ do

$d(w)\leftarrow d(t)+l_{tw}$

$\gamma(w)\leftarrow blue$

else

foreach  _$(t,w)\in\delta^{\text{out}}_{t}$ such that $d(w)\geq
d(t)+l^{-}_{tw}$_ do

$d(w)\leftarrow d(t)+l^{-}_{tw}$

$\gamma(w)\leftarrow red$

$S\leftarrow S\cup\\{t\\}$

return $\\{t\in V:\gamma(t)=blue\\}$

Algorithm 2 Computes the set of vertex $t$ such that $(u,v)$ is $t$-useless

######  Proposition 2.

Given a graph $G=(V,A),$ two arc-length functions $l^{-}$ and $l$ such that
$0\leq l^{-}_{a}\leq l_{a}$ for every arc $a\in A$, and an arc $(u,v)$,
Algorithm 2 computes in $O(|A|+|V|\log|V|)$ the set of vertex $t$ such that
$(u,v)$ is $t$-useless.

###  3.3 Operations to reduce the size of the graph

Removal of an arc from $W(t)$. Let $(u,v)\in W(t)$, since for all scenarios
$x\in\\{0,1\\}^{\Phi}$, $(u,v)$ does not belong to any shortest path from $u$
to $t,$ its removal does not affect the distance from any vertex to $t$. It is
clear, from the definition of ECA, that the removal of $(u,v)$ does not affect
the contribution of $t$ to ECA. Let $f^{x}_{t}(G):=\sum_{s\in
V}w_{s}w_{t}\Pi^{x}_{st}$ be the contribution to ECA of all the pairs having
sink $t$ in the graph $G$ when the probability of connection $\Pi^{x}_{st}$ is
computed under the scenario $x.$

######  Lemma 4.

Let $(u,v)\in W(t)$ and let $G^{\prime}$ be a graph obtained from $G$ by
removing arc $(u,v).$ Then, for all scenario $x\in\\{0,1\\}^{\Phi},$ it holds
that $f^{x}_{t}(G)=f^{x}_{t}(G^{\prime}).$

###### Proof.

For every scenario $x\in\\{0,1\\}^{\Phi}$, $(u,v)$ does not belong to any
shortest path from $u$ to $t.$ Therefore, the removal of $(u,v)$ cannot affect
the probability of connection $\Pi_{ut}$ for any vertex $t$, Thus
$f^{x}_{t}(G)=f^{x}_{t}(G^{\prime}).$ ∎

Contraction of an arc in $S(t)$. Now, assume that $(u,v)\notin\Phi,(u,v)\in
S(t)$. The contraction of $(u,v)$ consists in replacing every arc
$(w,u)\in\delta^{\text{in}}_{u}$ by an arc $(w,v)$ of length
$l^{\prime}_{wv}=l_{wu}+l_{uv}$ and by removing the vertex $u$ and all its
outgoing arcs. The weight of $u$ in $G$ is moved to the weight of $v$ in
$G^{\prime}.$ Namely, the weight of $v$ in the new graph is
$w^{\prime}_{v}=w_{v}+w_{u}\exp(-l_{uv}).$ Let $G^{\prime}$ be the graph
obtained from $G$ by contracting $(u,v).$ The next lemma establishes that the
contribution of $t$ to ECA in $G$ is equal to its contribution in
$G^{\prime}.$

(a)  (b)

Figure 2: (a) A graph $G$ before contraction of an arc $(u,v).$ (b) A graph
$G^{\prime}$ obtained from $G$ by contracting arc $(u,v)$. The weight
$w^{\prime}_{v}$ of $v$ in $G^{\prime}$ is equal to
$w_{v}+w_{u}\exp(-l_{uv}).$

######  Lemma 5.

Let $(u,v)\in S(t)$ and let $G^{\prime}$ be the graph obtained from $G$ by
contracting arc $(u,v)$ and modifying accordingly the weight of $w_{v}.$ For
every scenario $x\in\\{0,1\\}^{\Phi},$ $f^{x}_{t}(G)=f^{x}_{t}(G^{\prime}).$

###### Proof.

Let $s$ be a vertex of $G$ and $x$ be a scenario in $\\{0,1\\}^{\Phi}$. We
denote $f^{x}_{st}(G)$ the contribution to ECA of the pair $st$ with scenario
$x\in\\{0,1\\}^{\Phi}$ in $G.$ If $s\notin\\{u,v\\},$ it is easy to check
that, for any $x\in\\{0,1\\}^{\Phi}$ the length of the shortest path from $s$
to $t$ in $G$ is equal to the length of the shortest path from $s$ to $t$ in
$G^{\prime}.$ Moreover, the weight of $s$ is the same in $G$ and in
$G^{\prime}.$ Therefore, by the definition of ECA,
$f^{x}_{st}(G)=f^{x}_{st}(G^{\prime}),$ for any $x\in\\{0,1\\}^{\Phi}.$ On the
other hand, the contributions of the pairs $ut$ and $vt$ in $G$ sum to the
contribution of $v$ in $G^{\prime}.$ Indeed,

$\displaystyle f^{x}_{ut}(G)+f^{x}_{vt}(G)$ $\displaystyle=$ $\displaystyle
w_{u}\exp(-d_{x}(u,t))+w_{v}\exp(-d_{x}(v,t))$ $\displaystyle=$ $\displaystyle
w_{u}\exp(-(l_{uv}+d_{x}(v,t)))+w_{v}\exp(-d_{x}(v,t))$ $\displaystyle=$
$\displaystyle(w_{v}+w_{u}\exp(-l_{uv}))\exp(-d_{x}(v,t))$ $\displaystyle=$
$\displaystyle w^{\prime}(v)\exp(-d_{x}(v,t))$ $\displaystyle=$ $\displaystyle
f^{x}_{vt}(G^{\prime}).$

We conclude that, for any $x\in\\{0,1\\}^{\Phi},$ $f^{x}_{t}(G)=\sum_{s\in
V}f^{x}_{st}(G)=f^{x}_{ut}(G)+f^{x}_{vt}(G)+\sum_{s\in
V-\\{u,v\\}}f^{x}_{st}(G)=f^{x}_{vt}(G^{\prime})+\sum_{s\in
V-\\{u,v\\}}f^{x}_{st}(G^{\prime})=f^{x}_{t}(G^{\prime}).$ ∎

Graph reduction. Let $G_{t}$ be the graph obtained from $G$ by deleting every
arc of $W(t)$ and contracting every arc of $S(t)$. By induction the
contribution of $t$ is preserved. The preprocessing step consists in computing
the graph $G_{t}$ for every vertex $t.$ For that, we compute $F(u,v)$ for
every arc $(u,v)$ in $O(|A|\cdot(|A|+|V|\log|V|)).$ Then we transpose the
representation to obtain $S(t)$ for each vertex $t.$ Similarly, we compute
$W(t)$ for each vertex $t$ within the same complexity. The experiments of
Section 5 show that, when the number of arcs that can be protected is much
smaller than the total number of arcs, replacing $G$ by $G_{t}$ in the
construction of the mixed integer program reduces significantly the size of
the MIP formulation and the running times.

##  4 Greedy algorithms and their limits

In the lack of an efficient method to compute an optimal solution of BC-ECA-
Opt, ecologists often use greedy algorithms to compute sub-optimal solutions
[11]. In this section, we present four commonly used greedy algorithms and
highlight their pathological cases. In section 5.3, we compare the quality of
the solutions obtained with these greedy algorithms with the optimal solution
on four case studies.

The Incremental Greedy (IG) algorithm starts from the graph with no improved
element. At each step $i$, the algorithm selects the element $e$ with the
greatest ratio $\delta_{e}^{i}/c_{e}$ until no more element fits in the
budget. Here, $\delta_{e}^{i}$ denotes the difference between the value of ECA
with and without the improvement of the element $e$ at the step $i.$ As usual,
$c_{e}$ is the cost of improving the element $e$. The element $e$ can be
either an arc or a vertex.

The Decremental Greedy (DG) algorithm, similar to the Zonation algorithm [15],
algorithm starts from the graph with all improvements performed and
iteratively removes the improvement of the element $e$ with the smallest ratio
$\delta_{e}^{i}/c_{e}$. DG finishes with incremental steps to ensure there is
no free budget left. These algorithms perform at most $|\Phi|$ steps and at
each step $i$ need to compute $\delta_{e}^{i}$ for each element $e$. It is
easy to implement IG in $O(|V|^{3}+|\Phi|^{2}\cdot|V|^{2})$ by using an all
pair shortest path algorithm to compute in $O(|V|^{3})$ the initial distance
matrix and then by performing $|\Phi|$ steps in which the computation of
$\delta^{i}_{e}$ for each arc $e\in\Phi$ takes $O(|V|^{2}).$ Indeed, when we
decrease the length of an edge we can update the distance matrix in
$O(|V|^{2})$. For the implementation of DG, when we increase the length of an
edge we cannot update the distance matrix as easily as for IG. We can
recompute in $|V|^{3}$ the whole distance matrix for computing each
$\delta^{i}_{e}$ and the algorithm runs in $O(|\Phi|^{2}\cdot|V|^{3})$.
Dynamically updating shortest path lengths would improve computational
complexity [8]. These complexities are already too large for the practical
instances handled by ecologists which can have few thousands of patches. Most
of the studies using the PC or ECA indicators use simpler algorithms that we
call Static Increasing (SI) and Static Decreasing (SD). These algorithms are
variants of the greedy algorithms which do not recompute the ratio
$\delta_{e}/c_{e}$ of each element $e$ at each step and thus are faster but do
not account for cumulative effects nor redundancies.

Below, we provide instances on which IG and DG performs poorly compared to an
optimal solution. On these instances, it is easy to check that the solutions
returned by SI and SD are not better than the solutions returned by IG and DG.
In the following instances, all arcs have a probability $1$ if improved and
$0$ otherwise and have unitary costs. Recall that a spider is a tree
consisting of several paths glued together on a central vertex (Fig. 3a, 4a
and 5a).

Bad case for IG: The graph is a spider with $2k$ branches: $k$ long branches
with two edges, an intermediate node of weight $0$ and a leaf node of weight
$1$, and $k$ short branches consisting of a single edge with a leaf of very
small weight $\epsilon>0$, see Fig.2 (a). All branches are connected to a
central node of weight $1$. IG performs poorly on this instance. Indeed, IG is
tricked into purchasing short branches with very small ECA improvement because
purchasing an edge of a long branch alone does not increase ECA at all. An
optimal solution results in larger value of ECA by improving pairs of arcs of
long branches.

In this case, IG does not perform well while DG finds an optimal solution
computed by the MIP solver except when the budget is $1.$ In this case, the
reverse occurs: DG performs badly while IG is optimal. Indeed, DG realizes
that the budget is not sufficient to improve two arcs of a long branch only
after removing the improvements of all short branches.

(a)

(b)

Figure 3: An instance on which Incremental Greedy fails. (a) the graph of the
IG bad case with $k=4$, (b) ratio of the increase in ECA between the solutions
returned by IG and DG and an optimal solution for several budgets.

Bad case for DG: The graph is obtained from a star with $k+1$ branches by
replacing one branch by a path of length $k.$ The central node and all leaves
except the leaf of the path have weight $1.$ The leaf of the path has weight
$1+\epsilon$. The internal nodes of the path have weight $0$. DG performs
poorly on this instance because it removes one by one the branches of the star
for which $\delta_{e}/c_{e}=1$ before removing an edge of the path for which
$\delta_{e}/c_{e}=1+\epsilon.$ When the budget is at least $2,$ an optimal
solution removes all the edges of the path before removing an edge of another
branch.

(a)

(b)

Figure 4: An instance on which Decremental Greedy fails. (a) the graph of the
DG bad case with $k=5$, (b) ratio of the increase in ECA between the solutions
returned by IG and DG and an optimal solution for several budgets.

Bad case for IG and DG: The graph is a spider with $k+1$ branches. All
branches except one are paths of length 2 with an internal node of weight $0$
and a leaf of weight $1.$ The last branch is a path of length $2k$ with
internal nodes of weight $\epsilon>0$ and a leaf of weight $1+\epsilon$. All
branches intersect in a central node of weight $1$. In this case, both
Incremental and Decremental Greedy fail. On one hand, IG selects the edges of
the path of length $2k$ one by one and does not realize that by taking two
edges of a short branch it could improve much more ECA. On the other hand, DG
removes first the edges of the short branch because the weight of leaf of a
long branche is $1+\epsilon$ while the weight of the leaf of a short branch is
$1.$ Hence, DG and IG return the same law quality solution.

(a)

(b)

Figure 5: An instance on which both Incremental and Decremental Greedy
algorithms fail. (a) the graph of the IG and DG bad case with $k=5$, (b) ratio
of the increase in ECA between the solutions returned by IG and DG and an
optimal solution for several budgets.

The case of Figure 5 illustrates the fact that IG and DG do not provide any
constant approximation guarantee (even for trees), i.e. for any constant
$0<\alpha<1$ there exists an instance of BC-ECA-Opt such that $ALG<\alpha OPT$
where ALG is the value of ECA for the best solution among those returned by IG
and DG and $OPT$ is the value of ECA for an optimal solution.

##  5 Numerical experiments

In this section, we report on our computational experiments in order to
demonstrate the added benefit of our MIP formulation and preprocessing step.
We performed the numerical experiments on a desktop computer equipped with an
Intel(R) Core(TM) i7-8700k 4.8 gigahertz and 32 gigabytes of memory and
running Manjaro Linux release 21.2.4. We implemented our model as well as the
preprocessing and greedy algorithms in C++17 using Gurobi Optimizer [10]
version 9.1.1 with default settings for solving MIP formulations, the graph
library LEMON [9] version 1.3.1 for managing graph algorithms, and the library
TBB [19] version 2020.3 for multithreading the preprocessing and greedy
algorithms. Code is available at https://gitlab.lis-
lab.fr/francois.hamonic/landscape_opt_networks_submission.

###  5.1 Instances

Below, we briefly describe the case studies on which we conduct experiments.

Case study 1 consists in identifying among a set of 15 dams present on the
Aude river (France) those that need to be equipped with fish passes in order
to restore the river connectivity for trouts [22]. The graph is a tree of $45$
vertices with $88$ arcs of which $30$ represent damns and can be improved by
increasing their probability from $0$ to $0.8$.

Case study 2 consists in identifying the remnant forest patches that need to
be preserved from deforestation in the Montreal neighborhood (Canada) to
guaranty habitat connectivity for the wood frog [2]. The graph is a planar
graph of $598$ vertices and $989$ arcs whose $260$ vertices can be improved by
increasing their quality and $80$ arcs can be improved by increasing their
probability from $0$ to $1$.

Case study 3 consists in identifying street sections in which planting trees
can improve the connectivity of the urban canopy for the European red squirrel
in the city of Aix-en-Provence. The graph is a triangular grid of $6186$
vertices and $27818$ arcs where $47$ street sections, each with an average of
$90$ arcs, can be improved by increasing the probability of each arc $a$ from
$\pi_{a}$ to $\pi_{a}^{1/6}$.

Case study 4 consists in identifying wastelands that need to be preserved from
artificialization to maintain connectivity among urban parks and the
surrounding natural massifs in the city of Marseille for songbirds (e.g.
Eurasian blackcap). The baseline graph is a near complete graph of $297$ nodes
and $25024$ arcs, of which $100$ represent wastelands and can be improved by
increasing their probability from $0$ to $1$.

###  5.2 Scalability and benefits of the preprocessing

In this section we address the scalability of our approach and the added
benefits of our preprocessing step. For this purpose we execute our method on
about one hundred instances obtained from the four case studies by varying the
budget.

Figure 6: Execution times on the four case studies as a function of the budget
(missing points correspond to instances that do not finish within 10 hours)

For the Aude and Montreal cases, the preprocessing reduces the resolution time
by about 10 times (Fig. 6). Without preprocessing, the Aix and Marseille
instances are not solved by the optimizer within 10 hours, whereas with
preprocessing they become solvable in about half an hour and half a minute
respectively. In most unfinished instances, the optimizer reaches the optimal
solution but is not able to complete the exploration of the search space in
the allotted time.

case | MIP | preprocessed MIP
---|---|---
#var | #const | time | #var | #const | p. time | time
Aude | 4061 | 2551 | 120 ms | 1069 | 1055 | 3 ms | 20 ms
Montreal | 830530 | 262445 | 4 mins | 318848 | 167153 | 0.26 s | 19 s
Aix | 1748708 | 624841 | – | 555124 | 295010 |  $3$ s | 1600 s
Marseille | 4949825 | 78410 | – | 41676 | 22465 |  $0.9$ s | 7 s
Table 1: Comparison of the MIP and the preprocessed MIP according to the
number of variables (#var), the number of constraints (#const), the
preprocessing time (p. time) and the average computation time (time).

The preprocessing represents a small portion of the total computation time for
all case studies (Table. 1). The number of variables of the model is reduced
by about $75\%$ in the Aude case, $60\%$ in the Montreal case, $70\%$ in the
Aix case and $99\%$ in the Marseille case. This last number is explained by
the fact that the Marseille graph is near complete and a large proportion of
its arcs are $t$-useless for some vertex $t$. For constraints, the reduction
is $60\%$ for the Aude case, $33\%$ for the Montreal one, $53\%$ for the Aix
case and about $70\%$ for the Marseille case.

#wasteland | MIP | preprocessed MIP | DG
---|---|---|---
#var | #const | time | #var | #const | time | time
20 | 345775 | 18410 | 12 s | 10716 | 7197 |  $<1$ s |  $<1$ s
50 | 717901 | 36410 | 2 mins | 34613 | 21485 | 2 s | 7 s
80 | 1306597 | 59810 | 32 mins | 76281 | 42039 | 13 s | 30 s
110 |  | - |  | 132225 | 66759 | 1 min | 1 min 30 s
140 |  | - |  | 207999 | 96600 | 3 mins | 3 mins
170 |  | - |  | 308355 | 132701 | 28 mins | 7 mins
Table 2: Comparison of the MIP, the preprocessed MIP and DG according to the
number of variables (#var), the number of constraints (#const) and the time
(on average with 20 different budget values) it takes to solve the Marseille
instance with different numbers of wastelands

We see in Table 2 that the computation time of the MIP without preprocessing
increases very quickly. It takes more than 30 minutes on average for instances
with 80+ wastelands whereas the preprocessed one can be solved in less than 30
minutes with up to 170 wastelands. This is due to the preprocessing step that
significantly reduces the time required to solve the linear relaxation by
reducing the number of variables, constraints and non-zero entries of the
mixed integer program. The preprocessed MIP is faster than the greedy
algorithm for instances with at most 140 wastelands (the preprocessing step
was not used for greedy algorithms).

Finally, we run our preprocessing on 400 randomly generated instances from red
a model of the landscape around Montreal for hares of $8733$ vertices and
$18422$ arcs [2] to study the impact of the preprocessing on the MIP
formulation size. For building these instances, we take 20 connected subgraphs
of 500 nodes and for each graph we create 20 instances by randomly picking a
percentage of arcs whose probability could be increased from $\pi$ to
$\sqrt{\pi}$.

Figure 7: Box plots showing the percentage of constraints, variables and non-
zero entries removed by the preprocessing as a function of the percentage of
improvable arcs. The red line is the median, the dashed green line is the
mean, the box represents the values between the 25th and 75th percentiles and
the whiskers the min and max values.

Figure 7 shows a box plot ofthe percentage of reduction of the number of
constraints, variables and non-zero entries of the MIP formulation with
respect to the percentage of arcs that could be improved. The preprocessing
removes almost all the elements of the MIP formulation when the number of
improvable arcs arrives close to zero. This reduction decreases with the
number of arcs that can be improved. When $20\%$ of the arcs can be improved,
the preprocessing removes on average $80\%$ and at least $70\%$ of the model’s
variables, on average $60\%$ and at least $45\%$ of the model’s constraints,
and on average $75\%$ and at least $65\%$ of the model’s non-zero entries.
Even when $100\%$ of the arcs can be improved, the preprocessing reduces on
average by $35\%$ the model’s size. Since the gap between the 25th and 75th
percentiles does not exceed $17\%,$ the reduction seems to be robust. These
results look consistent with the one of Table 1. Indeed, in the case of Aix,
about $15\%$ of the arcs can be improved and our preprocessing reduces the
number of variables by $70\%$ and the number of constraints by $65\%.$

###  5.3 Quality of the solutions

Figure 8: Percentage gain in ECA achieved by the solutions of the different algorithms compared to the optimal solution for different budget values.  | IL | DL | IG | DG
---|---|---|---|---
min. | avg. | min. | avg. | min. | avg. | min. | avg.
Aude | 80.3 % | 94 % | 14.7 % | 83.9 % | 82.8 % | 94 % | 88.1 % | 97.2 %
Montreal | 92 % | 99.2 % | 97.7 % | 99.6 % | 98.4 % | 99.8 % | 98.4 % | 99.8 %
Aix | 81.6 % | 96.2 % | 65.2 % | 95.6 % | 81.6 % | 98.7 % | 54.1 % | 97.4 %
Marseille | 97.8 % | 99.5 % | 95 % | 99.3 % | 97.8 % | 99.6 % | 97.8 % | 99.6 %
Table 3: Minimum and average optimilaty ratio for each algorithm and case
study.

For each of the case studies and each of the four algorithms, there is at
least one budget value for which the quality of the solution is significantly
lower than the quality of the optimal solution, the greatest departures being
observed at lower budget values (Fig. 8, Table. 3). Greedy versions of
incremental and decremental algorithms perform on average better than their
static counterpart (Table. 3). The minimum and average optimality ratio in the
Aude and Aix cases is lower than in the other cases, for all algorithms
(Table. 3). Static and greedy algorithms are generally quite close to the
optimal solution ($5\%$ lower on average). However, all algorithms, whether
static or greedy, incremental or decremental, provide poor quality solutions
for some budget values (Table. 3).

##  6 Conclusion

This article introduces a new MIP formulation for BC-ECA-Opt and shows that
this formulation allows to optimally solve instances having up to 150 habitat
patches while previous formulations, such as those described in [29] are
limited to 30 patches. The preprocessing step reduces significantly the size
of the graphs on which a generalized flow has to be computed, thus enabling to
scale up to even larger instances. We showed that this preprocessing step
allows to greatly reduce the MIP formulation size and that its benefits
increase when the proportion of arcs whose lengths can change decrease. This
allows to tackle instances up to 300 habitat patches.

The optimum solutions obtained experimentallly are compared to the ones
returned by several greedy algorithms. Interestingly, we found that greedy
algorithms perform well in practice despite the arbitrary bad cases we
spotted. Therefore, greedy algorithms remains a reasonable choice for
instances too large to be solved optimally by an MIP solver. Our next goals
will be to experiment our approach on other practical instances of the problem
arising from different ecological contexts.

On the theoretical side, we would like to investigate the problem from the
point of view of an approximation algorithm. For instance, is it possible to
find reasonable assumptions under which greedy algorithms are guaranteed to
return a solution whose ECA value is at least a constant fraction of the
optimal ECA? If these assumptions are fulfilled by the real instances that we
considered, this would explain our experimental observations. Moreover, since
a polynomial time approximation scheme has been given in the case of trees
[28], it could also be interesting to know for which larger classes of graphs
constant factor approximation algorithms for this problem exist. Another
interesting question is to determine whether good solutions could be obtained
by decomposing geographically the problem, by solving independently a
subproblem for each region and then by reassembling the solutions. In this
case, a notion of fairness could help to allocate the budget among the regions
so that each region can enhance its own internal connectivity keeping a part
of the budget to enhance the connectivity between the regions. Since ECA is
based on the equivalence between a landscape and a patch, such a multilevel
optimization approach looks promising.

##  7 Acknowledgments

We are grateful to the referees for a careful reading and many useful comments
and suggestions. The research on this paper was supported by Région Sud
Provence-Alpes-Côte d’Azur, Natural Solutions and the ANR project DISTANCIA
(ANR-17-CE40-0015). The Aix case study belongs to the Baum program
(Biodiversity Urban Development Morphology) supported by the PUCA, the OFB and
the DGALN. For stimulating exchanges on the case studies, we also thank:
Patrick Bayle, Simon Blanchet, Andrew Gonzalez, Maria Dumitru, Jérôme Prunier,
Bronwyn Rayfield, Benoit Romeyer and Keoni Saint-Pé.

## References

  * [1] Ravindra K. Ahuja, Thomas L. Magnanti and James B. Orlin  “Network Flows: Theory, Algorithms, and Applications”  USA: Prentice-Hall, Inc., 1993
  * [2] Cécile H. Albert, Bronwyn Rayfield, Maria Dumitru and Andrew Gonzalez  “Applying network theory to prioritize multispecies habitat networks that are robust to climate and land-use change”  In _Conservation Biology_ 31, 2017, pp. 1383–1396
  * [3] Marcelo Awade, Danilo Boscolo and Jean Paul Metzger  “Using binary and probabilistic habitat availability indices derived from graph theory to model bird occurrence in fragmented forests”  In _Landscape Ecology_ 27 Springer, 2012, pp. 185–198
  * [4] A. Billionnet  “Designing protected area networks”  EDP Sciences, Paris, France, 2021
  * [5] Alain Billionnet  “Mathematical optimization ideas for biodiversity conservation”  In _European Journal of Operational Research_ 231.3, 2013, pp. 514–534  DOI: https://doi.org/10.1016/j.ejor.2013.03.025
  * [6] E.. Brondizio, J. Settele, S. Díaz and H.. Ngo  “Global assessment report on biodiversity and ecosystem services of the Intergovernmental Science- Policy Platform on Biodiversity and Ecosystem Services”  IPBES, Bonn, Germany, 2019
  * [7] Daniele Catanzaro, Martine Labbé and Martha Salazar-Neumann  “Reduction approaches for robust shortest path problems”  In _Computers & OR_ 38, 2011, pp. 1610–1619
  * [8] Camil Demetrescu and Giuseppe F Italiano  “A new approach to dynamic all pairs shortest paths”  In _Journal of the ACM (JACM)_ 51.6 ACM New York, NY, USA, 2004, pp. 968–992
  * [9] Balázs Dezs, Alpár Jüttner and Péter Kovács  “LEMON - an Open Source C++ Graph Template Library”  In _Electron. Notes Theor. Comput. Sci._ 264.5 Elsevier Science Publishers B. V., 2011, pp. 23–45
  * [10] LLC Gurobi Optimization  “Gurobi Optimizer Reference Manual”, 2022  URL: https://www.gurobi.com
  * [11] Jeffrey O Hanson, Richard Schuster, Matthew Strimas-Mackey and Joseph R Bennett  “Optimality in prioritizing conservation projects”  In _Methods in Ecology and Evolution_ 10.10 Wiley Online Library, 2019, pp. 1655–1663
  * [12] Jochen AG Jaeger et al.  “Landscape fragmentation in Europe”  Publications office of the European Environmental Agency, Luxembourg, 2011, pp. 157–198
  * [13] Garth P. McCormick  “Computability of global solutions to factorable nonconvex programs: Part I - Convex underestimating problems”  In _Mathematical Programing_ 10.1, 1976, pp. 147–175  DOI: 10.1007/BF01580665
  * [14] Brad H. McRae, Brett G. Dickson, Timothy H. Keitt and Viral B. Shah  “USING CIRCUIT THEORY TO MODEL CONNECTIVITY IN ECOLOGY, EVOLUTION, AND CONSERVATION”  In _Ecology_ 89, 2008, pp. 2712–2724
  * [15] Atte Moilanen et al.  “Prioritizing multiple-use landscapes for conservation: methods for large multi-species planning problems”  In _Proceedings of the Royal Society B: Biological Sciences_ 272.1575 The Royal Society London, 2005, pp. 1885–1891
  * [16] Lucía Pascual-Hortal and Santiago Saura  “Comparison and development of new graph-based landscape connectivity indices: Towards the priorization of habitat patches and corridors for conservation”  In _Landscape Ecology_ 21 Springer, 2006, pp. 959–967
  * [17] Juliana Pereira, Santiago Saura and Ferenc Jordán  “Single-node vs. multi-node centrality in landscape graph analysis: Key habitat patches and their protection for 20 bird species in NE Spain”  In _Methods in Ecology and Evolution_ 8, 2017, pp. 1458–1467
  * [18] Miguel Pereira, Pedro Segurado and Nuno Neves  “Using spatial network structure in landscape management and planning: A case study with pond turtles”  In _Landscape and Urban Planning_ 100, 2011, pp. 67–76
  * [19] Chuck Pheatt  “Intel® threading building blocks”  In _Journal of Computing Sciences in Colleges_ 23.4 Consortium for Computing Sciences in Colleges, 2008, pp. 298–298
  * [20] Lidón Rubio, Örjan Bodin, Lluís Brotons and Santiago Saura  “Connectivity conservation priorities for individual patches evaluated in the present landscape: How durable and effective are they in the long term?”  In _Ecography_ 38, 2015, pp. 782–791
  * [21] Deborah A. Rudnick et al.  “The role of landscape connectivity in planning and implementing conservation and restoration priorities”  In _Issues in Ecology_ 16, 2012, pp. 1–20
  * [22] Keoni Saint-Pé  “In situ quantification of brown trout movements”, 2019
  * [23] Santiago Saura and Lucía Pascual-Hortal  “A new habitat availability index to integrate connectivity in landscape conservation planning: Comparison with existing indices and application to a case study”  In _Landscape and Urban Planning_ 83, 2007, pp. 91–103
  * [24] Santiago Saura and Josep Torne  “Conefor Sensinode 2.2: a software package for quantifying the importance of habitat patches for landscape connectivity”  In _Environmental modelling & software_ 24 Elsevier, 2009, pp. 135–139
  * [25] Santiago Saura, Christine Estreguil, Coralie Mouton and Mónica Rodríguez-Freire  “Network Analysis to Assess Landscape Connectivity Trends: Application to European Forests (1990-2000)”  In _Ecological Indicators_ 11, 2011, pp. 407–416
  * [26] Philip D. Taylor, Lenore Fahrig, Kringen Henein and Gray Merriam  “Connectivity Is a Vital Element of Landscape Structure”  In _Oikos_ 68 [Nordic Society Oikos, Wiley], 1993, pp. 571–573
  * [27] Dean Urban and Timothy Keitt  “LANDSCAPE CONNECTIVITY: A GRAPH-THEORETIC PERSPECTIVE”  In _Ecology_ 82, 2001, pp. 1205–1218
  * [28] Xiaojian Wu, Daniel R Sheldon and Shlomo Zilberstein  “Stochastic Network Design in Bidirected Trees”  In _Advances in Neural Information Processing Systems_ 27 Curran Associates, Inc., 2014, pp. 882–890
  * [29] Yexiang Xue et al.  “Dynamic Optimization of Landscape Connectivity Embedding Spatial-capture-recapture Information”  In _31st AAAI Conference on Artificial Intelligence_ 31, 2017, pp. 4552–4558