Ali El

Olivier Brun
email: brun@laas.fr

A Game-Theoretic Algorithm for the Joint Routing and VNF Placement Problem

Keywords: Service Function Chain, Virtual Network Function, Slice as a Service, Game Theory, Optimization

Network Function Virtualization (NFV) simplifies the deployment of network services by leveraging virtualization technologies to make the management of network functions more flexible and cost efficient. The deployment of these services requires the allocation of Virtual Network Function -Forwarding Graphs (VNF-FGs), which implies placing and chaining VNFs according to the requests of VNF-FGs. In this paper, we consider the offline allocation of VNF-FGs problem to improve resource utilization and reduce total costs. We focus on how VNF-FG demands are routed so as to optimize resource utilization without adding capacity to the infrastructure. Given a non-linear cost function associated to each network resource, we formulate the problem as a non-linear single-path routing problem in an extended graph. Then, we propose to adapt a single-path routing heuristic algorithm inspired from game theory to solve it. We show that this algorithm converges and establishes its approximation ratio in a number of cases. Experimental results obtained for different network topologies and different cost functions show that this algorithm provides very good quality solutions with substantially lower computing times compared to the optimal solution.

Introduction

Network Function Virtualization (NFV) opens up the possibility to replace hardwired middleboxes implementing specific network functions (e.g., firewalls or Web proxies) by Virtual Network Functions (VNFs) running in data centers and operated as cloud services. In addition to reducing the cost and complexity of the network, this approach enables to deploy and compose VNFs so as to create network slices on demand. Network slices are isolated logical networks coexisting simultaneously on the same Physical Substrate Network (PSN), each one being tailored to the needs of a specific application. The VNFs composing the network slice as well as the virtual links established between them are described by a VNF Forwarding Graph (VNF-FG). A VNF-FG, according to ETSI, is a directed graph of an ordered set of VNFs that compose an end-to-end service [START_REF] Nfv | Network Functions Virtualisation (NFV): Architectural framework[END_REF].

The logical service graph described by a VNF-FG is decoupled from the network and has to be mapped to existing PSN capabilities [START_REF] Herrera | Resource allocation in NFV: A comprehensive survey[END_REF]. This implies not only the fulfillment of slice requirements in terms of Quality of Service (QoS), but also considering the constraints of the underlying infrastructure. Two types of decisions have to be made: (a) where to run the VNFs, and (b) how to interconnect them in the physical network, taking into account specific VNF ordering requirements. They are usually made so as to optimize a certain performance objective, which may be related for instance to the resource utilization of the PSN, while satisfying a number of technical constraints (e.g., some VNFs should be placed in the vicinity of end users). This problem, which is known as the VNF Forwarding Graph Embedding (VNF-FGE) problem, is an extension of the Virtual Network Embedding problem [START_REF] Fischer | Virtual network embedding: A survey[END_REF][START_REF] Laghrissi | A survey on the placement of virtual resources and virtual network functions[END_REF] which is known to be NP-hard [START_REF] Amaldi | On the computational complexity of the virtual network embedding problem[END_REF]. The VNF-FGE problem has been studied in the online context, where requests to create network slices are submitted one after the other [START_REF] Soualah | A green VNF-FG embedding algorithm[END_REF][START_REF] Mechtri | VNF placement and chaining in distributed cloud[END_REF][START_REF] Carpio | Replication of virtual network functions: Optimizing link utilization and resource costs[END_REF][START_REF] Quang | Single and multi-domain adaptive allocation algorithms for VNF forwarding graph embedding[END_REF][START_REF] Amine | Shortening the deployment time of SFCs by adaptively querying resource providers[END_REF], and in the offline context where multiple VNF-FGs have to be embedded in the physical network at the same time [START_REF] Luizelli | Piecing together the NFV provisioning puzzle: Efficient placement and chaining of virtual network functions[END_REF][START_REF] Tastevin | A graph approach to placement of service functions chains[END_REF][START_REF] Zhong | Cost-aware service function chaining with reliability guarantees in NFV-enabled inter-dc network[END_REF][START_REF] Guo | Cost-aware placement and chaining of service function chain with VNF instance sharing[END_REF].

In this work, we study the offline joint routing and VNF placement problem with non-linear cost functions to reduce total costs and improve resource utilization under physical infrastructure capacity constraints. In particular, we consider a physical infrastructure with known transmission and capacity resources used to set up logical service requests. As these requests arrive one after the other in time, without knowing neither the future requests nor the logical networks that will be removed in the future, the individual placement and routing decisions made are necessarily sub-optimal in hindsight, which may result in poor resource utilisation. Under this model, we focus in this paper on the re-optimization of mapping existing logical networks onto the physical infrastructure. This will allow for a better load distribution which in turn will increase the chance to accommodate new logical networks without having to add additional infrastructure resources capacity.

Despite the previous efforts on the VNF-FGE problem, little considered the re-optimization of resources utilization on the physical infrastructure with non-linear cost functions. In [START_REF] Luizelli | Piecing together the NFV provisioning puzzle: Efficient placement and chaining of virtual network functions[END_REF], the authors aimed at minimizing the number of VNF instances mapped on the infrastructure to improve resource utilization. However, the work ignored the routing cost from connecting the VNFs of the VNF-FG request. Inspired by the previous work to share the same VNF by multiple VNF-FGs, the work in [START_REF] Tastevin | A graph approach to placement of service functions chains[END_REF] and [START_REF] Zhong | Cost-aware service function chaining with reliability guarantees in NFV-enabled inter-dc network[END_REF] studied the offline VNF-FGE problem to reduce both deployment and forwarding traffic costs. In [START_REF] Tastevin | A graph approach to placement of service functions chains[END_REF], the authors proposed a graph-based heuristic inspired from the Viterbi algorithm to reduce the deployment cost of requested VNF-FG by reducing the number of PoPs hosting the VNF instances. However, this solution results in an unbalanced resource utilization. In [START_REF] Zhong | Cost-aware service function chaining with reliability guarantees in NFV-enabled inter-dc network[END_REF], the authors focused at minimizing the cost of routing and deploying the requested VNFs, while guaranteeing a certain level of reliability. The work however, does not consider an efficient use of resources for future requests. The closest to our work is [START_REF] Guo | Cost-aware placement and chaining of service function chain with VNF instance sharing[END_REF] that focuses on minimizing both resource utilization and total placement costs. Despite the similarities with our work, the model considered is different from ours. First, the authors considered an infrastructure with only PoP nodes that are used to host VNFs. Hence, their work do not consider the routing problem. Second, the authors considered the possibility to deploy only one type of VNF in a PoP. In addition to the above limitations, none of the previous works considered nonlinear cost functions which are essential if one wishes to express that the cost per unit capacity of a resource grows with its utilization rate.

In the setting that we consider, we are given a physical network whose links have known transmission capacities and in which some nodes provide known computing resources. A nonlinear cost function is associated to each network resource and allows to measure its congestion level. A set of traffic demands flow through the network, each traffic demand being characterised by its source, destination, volume and the sequence of VNFs it must flow through. The volume of a traffic demand may change after each VNF, due to, for example, encryption or data aggregation operations. An initial solution specifies, for each traffic demand, on which compute nodes the VNFs are executed and the path taken by the demand in the network between these VNFs. We show that the VNF placement and chaining problem can be formulated as a non-linear single-path routing problem in an extended graph. To solve this problem, we adapt a single-path routing algorithm inspired from Game Theory which was proposed in [START_REF] Brun | A penalized best-response algorithm for non-linear singlepath routing problems[END_REF]. We show that this algorithm converges and establishes its approximation ratio in a number of cases. Experimental results obtained for different network topologies and different cost functions show that this algorithm provides very good quality solutions with a rather modest computation time.

The paper is organized as follows. We describe the problem addressed in this paper and introduce our main notations in Section 2. In Section 3, we introduce the concept of expanded network. The problem is then formulated as a non-linear single path routing problem in the expanded network in Section 4. We present the proposed algorithm in Section 5. Section 6 is devoted to the performance evaluation of this algorithm using diverse network topologies and cost functions. Finally, some conclusions are drawn in Section 7.

Problem Statement

We are given a physical network represented by a directed graph G = (V, E), where V is the set of nodes and E is the set of links. A subset D ⊆ V of the nodes have the required compute and storage resources to host VNF. These nodes will be called the function nodes in the following. They are used to execute one or more VNFs out of a set F of K VNFs, where function f ∈ F is available at function nodes v ∈ D(f ) ⊆ D. A function node corresponds to a logical entity representing reserved computing and storage capacity in a data center. Multiple VNFs may be deployed at the same function node, and a VNF may be replicated at different function nodes. The transmission links of the network and the function nodes are the critical resources in our problem and in the following we shall denote by R = E ∪ D the set of these resources.

We are also given a set W of traffic demands. Each traffic demand w ∈ W is characterized by its source node s w , its destination node t w and the sequence F w = f w 1 , f w 2 , . . . , f w Kw of network functions it must flow through in that prescribed order. The sequence is called the Service Chain associated to the traffic demand. The processing path of demand w may therefore be decomposed in different stages, where stage 0 corresponds to the transmission of the original traffic volume from the source to VNF f w 1 (excluded), stage 1 represents the processing at VNF f w 1 followed by the transmission to VNF f w 2 , etc. The volume of a traffic demand may change at each stage due to, e.g., packet encryption/decryption or to data aggregation operations. Moreover, the processing capacity required per unit of traffic may be different for one VNF to another. We therefore define b w k,r as the capacity required from resource r for processing demand w at stage k. If r is a transmission link e ∈ E, then this capacity just represents the volume of traffic send from VNF f w k to VNF f w k+1 for 0 < k < K(w), and the volume of traffic from s w (resp. VNF f w K(w) ) to VNF f w 1 (resp. t w ) for k = 0 (resp. = K(w)). This volume of traffic should be the same for all links e and is expressed in bit/s. Similarly, if r is a function node v ∈ D(f w k ), the capacity b w k,r represents the (fractional) number of cores required for processing the traffic of demand w by VNF f w k . Our notations are summarized in Table 1.

Figure 1 provides a simple example of the setting considered in the paper, in which there are two traffic demands. The service chain for the traffic demand

s 1 -t 1 is F 1 = (f 1 , f 2 ), whereas the service chain for the traffic demand s 2 -t 2 is F 2 = (f 1 ). The volume of traffic demand s 1 -t 1 is b 1 0 = 20
units of traffic and the volume of traffic demand s 2 -t 2 is 10 units. To simplify, we assume that the volumes do not change along the paths from the sources to the destinations. We also assume that VNF f 1 requires 1.5 cores per unit of traffic, whereas VNF f 2 requires 1.0 core per unit of traffic. The capacity of all links are taken to be c e = 60 units of traffic, and the processing capacity of the function nodes are c D = 70 and c E = 40. The figure illustrates two possible deployment scenarios. In the first scenario (see Fig. 1a), all instances of all VNFs are deployed at node D, whereas in the second scenario (see Fig. 1b) the instances of VNF 1 required for processing the traffic demand s 1 -t 1 are executed at node E. We analyze below each scenario, assuming that the cost function for resource r is yr cr-yr , where y r is the traffic volume on resource r and c r is the capacity of resource r:

(a) Scenario 1: In this scenario, the utilization rates of the network links are as follows: 1 3 for one link, 1 6 for two links and 1 2 for four links. This yields a total link cost of 1

× 1/3 1-1/3 + s 2 b 1 0 = 20 b 2 0 = 10 f 1 s 1 t 1 C B D E f 1 , f 2 t 2 A
(a) Deployment scenario 1 in which all VNFs are executed at node D.

s 2 b 1 0 = 20 b 2 0 = 10 f 1 s 1 t 1 C B D E f 1 , f 2 t 2 A
(b) Deployment scenario 2 in which instances of VNF 1 are executed at node E for the traffic demand s1 -t1.

Figure 1: A simple example network , in which we assume that VNF f 1 may be executed at nodes D and E, whereas VNF f 2 is only available at node D.

2 × 1/6 1-1/6 + 4 × 1/2
1-1/2 = 4.9. However, the total capacity required at node D is 65 cores ((20 + 10) × 1.5 for VNF f 1 , and 20 × 1 for VNF f 2 ), yielding a total cost for function nodes equals to 65 70-65 = 13. Hence the total network cost for this deployment scenario is 17.9.

(a) Scenario 2: In this scenario, the utilization rates of the network links are as follows: 1 3 for seven links, 1 6 for three links and 1 2 for three other links. The total link cost is therefore 7

× 1/3 1-1/3 + 3 × 1/6 1-1/6 + 3 × 1/2 1-1/2 = 7.1.
The computing capacity required at node D is 35 cores: 10 × 1.5 for VNF f 1 , and 20 × 1 for VNF f 2 . The computing capacity required at node E is 20 × 1.5 = 30 cores. The total computing cost of function nodes is hence 35 70-35 + 30 40-30 = 4. We conclude that the total network cost in this scenario is 11.1, which shows that it achieves a better resource utilization.

The problem at hands amounts to choosing a single path through the network for each traffic demand such that the function nodes are visited in the prescribed order and so that a certain network cost function is minimized. We shall shortly discuss the network cost function to be optimized, but we would first like to emphasize that a first difficulty is related to the modelling of precedence constraints on the order in which the function nodes must be visited by each traffic request. The key observation is that the overall path from s w to t w for traffic demand w can be viewed as a collection of segments, the first segment connecting s w to one of the function nodes hosting f w 1 through some intermediate nodes, the second segment connecting the end node of the first segment with one of the function nodes hosting f w 2 , etc, until the destination node t w is reached. These segments are connected by decision points where we choose which nodes will compute the functions in F w . As observed in [START_REF] Nguyen | Routing via functions in virtual networks: The curse of choices[END_REF], it is possible to construct an expanded network to model both the segment construction and the function selection at the nodes. We shall discuss the construction of this expanded network in Section 3. In Section 4, we formulate the single-path routing problem considered in this paper using the concept of expanded network.

Expanded Network

The expanded network G * is made up of K layers. Each layer G (k) = (V (k) , E (k) ), k = 0, . . . , K, is a copy of of the original network G = (V, E). We shall use the notation v (k) (resp. e (k) ) to refer to node v ∈ V (resp. link e ∈ E) at layer k. We use each layer to find one segment of the total path for each traffic demand w. Consecutive layers are connected by artificial edges which are used to model the function selection at the nodes. More precisely, for each function node v ∈ D and each layer k = 0, . . . , K -1, there is an artificial directed edge from node v (k) to node v (k+1) . The interpretation is that a path for traffic demand w goes through this artificial link if to refer to the artificial edge (v (k) , v (k+1) ) and we shall denote the set of artificial edges between layers k = 0, . . . , K -1 and k + 1 by

f1 f1 f1, f2 s (0) 1 s (0) 2 t (0) 1 t (0) 2 A (0) C (0) B (0) D (0) E (0) s (1) 1 s (1) 2 t (1) 1 t (1) 2 A (1) C (1) B (1) D (1) E (1) s (2) 1 s (2) 2 t (2) 1 t (2) 2 A (2) C (2) B (2) D (2) E (2) f1, f2
L (k) = ℓ (k) v : v ∈ D . In summary, the set of nodes of the expanded network G * is V * = K k=0 V (k) ,
whereas its set of edges is

E * = K k=0 E (k) ∪ K-1 k=0 L (k) .
The construction of the expanded network is shown in Figure 2 for the example network of Figure 1.

In the following, we say that an artificial edge ℓ

(k) v ∈ L (k) is feasible for demand w if and only if v ∈ D(f w k )
, that is, it connects two consecutive copies of a node v hosting function f w k . A feasible path π w ⊂ E * for the traffic demand w ∈ W is then defined as a path from node s (0) w to node t (Kw) w in the expanded network such that all artificial edges in the path are feasible for demand w. Figure 2 shows two feasible paths for the demand from node s 1 to node t 1 . Note that the artificial edges appearing in a path specify the function nodes at which the VNFs are executed, while the other edges in the path specify the path segments used to reach these function nodes. The feasibility of a path ensures that the network functions are executed in the prescribed order.

Given a set π ⊂ E * of edges in the expanded network and a link e ∈ E of the original network, we define the constant δ e,k π as 1 if e (k) ∈ π, and 0 otherwise. Similarly, for any set π ⊆ E * and any function node v ∈ D, we define the constant δ v,k π as 1 if the artificial edge ℓ

(k)
v ∈ π, and 0 otherwise. In the following, the notation δ r,k π will refer to one of the above-defined constants, depending on whether resource r ∈ R is a transmission link e or a function node v. We shall say for short that resource r appears in path π at layer k of the expanded network and write (r, k) ∈ π when δ r,k π = 1. Similarly, we write r ∈ π if δ r,k π = 1 for some k.

Single-Path Routing Problem

In the following, we shall assume that we are given a set Π w of candidate feasible paths through the expanded network for each demand w ∈ W. A feasible solution for the joint routing and VNF placement problem is then defined as a vector π = (π w ) ∈ Π, where π w is the path assigned to traffic demand w and Π = w∈W Π w . Given such a feasible solution π, the traffic sent by demand w on resource r ∈ R of the original network is:

y w r (π w ) = K k=0 δ r,k πw b w k,r , (1)
from which it follows that the traffic flowing on resource r is y r (π) = w∈W y w r (π w ). We assume that to each resource r ∈ R of the physical network is associated a non-decreasing function ϕ r : R + → R + and that the cost of this resource has the form y r ϕ r (y r ). The function ϕ r () may be interpreted as the cost per unit of capacity of resource r, and it may depend on the total traffic flowing on that resource. A typical example is the Kleinrock function ϕ r (y r ) = 1/(c r -y r ), where c r is a constant representing the (transmission or processing) capacity of resource r. This function assumes that the total traffic on a resource is smaller than its capacity and that the cost per unit of capacity grows boundlessly as the former approaches the latter. Other examples are the linear function ϕ r (y r ) = y r /c r or the quadratic function ϕ r (y r ) = y 2 r /c 2 r , which are valid even for y r ≥ c r . In the above examples, the cost per unit capacity grows as the traffic y r on the resource increases, in contrast with linear models in which ϕ r is assumed to be a constant function. We emphasize that the functions ϕ r () may differ for different types of resources or even from one resource to the other.

The goal is to find a feasible solution π ∈ Π that minimizes the network cost F (π) = r∈R y r (π) ϕ r (y r (π)). Formally, the problem is as follows:

minimize

F (π) = r∈R y r (π) ϕ r (y r (π)) (OPT)
subject to:

π w ∈ Π w , w ∈ W. (2)
It is worth mentioning that problem (OPT) can be cast as a 0-1 mathematical programming problem by introducing the following binary variables:

x w,π = 1 if demand w is routed along path π in G * , 0 otherwise, (3)
so that it becomes equivalent to: minimize r∈R y r ϕ r (y r ) subject to:

y r = w∈W π∈Πw K k=0 δ r,k π b w k,r
x w,π , for all r ∈ R,

x w,π = 1, for all w ∈ W,

x w,π ∈ {0, 1}, for all π ∈ Π w and w ∈ W.

Except when the functions ϕ r () are constant functions, the above problem belongs to the class of non-linear mathematical programs with integer variables, a class of mathematical programs which are known to be extremely hard to solve. This motivates the development of an efficient approximation method for solving problem (OPT).

Penalized Best-Response Algorithm

As discussed in Section 4, the joint routing and VNF placement problem can be cast as a non-linear single path routing problem in the expanded network. The problem differs however significantly from traditional single-path routing problems as on one hand the same resource appears at different layers of the expanded network and on the other hand the volume of a traffic demand may change from one layer to the other. Nevertheless the heuristic algorithm that we propose in this section is directly inspired from an algorithm proposed in [START_REF] Brun | A penalized best-response algorithm for non-linear singlepath routing problems[END_REF] for solving such problems. The idea of the algorithm is to view the traffic demands as the players of a non-cooperative game in which each each player independently optimizes its own objective function. Starting from an arbitrary initial feasible solution, the algorithm then mimics the best-response dynamics of the game, that is, the players take turns in some order to adapt their strategy based on the most recent known strategy of the others. The players objective functions are designed in such a way that (a) the best-response dynamics converges to a Nash equilibrium of the game in a finite number of steps, (b) all optimal solutions of problem (OPT) are Nash equilibria of the game, and (c) an upper bound on the approximation ratio of the algorithm can be proven in some cases.

Let us think of the traffic demands as the players of the game. The strategy of player w ∈ W is the path π w it chooses in the set Π w , and a strategy profile of the game is a feasible solution to problem (OPT), that is, a vector π ∈ Π. Given the strategy of the other players π -w = (π u ) u̸ =w , we shall assume that the player w seeks to solve the following problem:

minimize π∈Πw c w (π, π -w ) = f w (π, π -w ) + p w (π, π -w ), (OPT-w)
where the value f w (π, π -w ) associated to path π by player w reflects the cost of this path, whereas the term p w (π, π -w ) is a penalty term measuring the impact of player w's choice on other players. More precisely, we shall assume that the selfish objective function of player w is:

f w (π, π -w ) = r∈π y w r (π) ϕ r y -w r + y w r (π) , (7)
where y -w r = u̸ =w y u r (π u ) is the total traffic sent over resource r by all the other players u ̸ = w. The term f w (π, π -w ) thus represents the cost of path π for player w, given the fixed strategies π -w of the other players.

In order to define the penalty term p w (π, π -w ), assume that player w chooses path π ∈ Π w and consider a resource r appearing in the path π at some layer. Then, for all players u ̸ = w, the cost of this resource increases by:

y u r (π u ) ϕ r y -w r + y w r (π) -ϕ r y -w r .
As a consequence, we define the penalty term p w (π, π -w ) as follows: With ( 7) and ( 8), the objective function c w (π, π -w ) = f w (π, π -w ) + p w (π, π -w ) of player w is perfectly defined. Given the strategies π -w of the other players, the path π minimizing c w (π, π -w ) in problem (OPT-w) is known as the best response of player w.

The pseudocode of our heuristic is given in Algorithm 1. Note that it is fully similar to the algorithm in [START_REF] Brun | A penalized best-response algorithm for non-linear singlepath routing problems[END_REF], the only difference being in the way the traffic flowing on a resource is computed. The algorithm starts from an initial feasible solution π (0) . At each iteration n, the players update their strategies in a given order by computing their best responses to the strategies of the others (lines 3 -5). Note that a player w deviates from its strategy π n) . The algorithm stops when no player can decrease its cost by unilaterally deviating from its strategy, that is, π (n+1) = π (n) , which means that a Nash equilibrium has been reached (line 7).

(n) w at iteration n to a new strategy π ′ if and only if c w π ′ , π (n) -w < c w π (
Theorem 1 below states the main properties of Algorithm 1. The proofs of properties (a), (b) and (c) are straightforward adaptations of those of similar results in [START_REF] Brun | A penalized best-response algorithm for non-linear singlepath routing problems[END_REF]. The proof of property (d) is notably more involved, though it follows the same lines as the proof of Theorem 3 in [START_REF] Brun | A penalized best-response algorithm for non-linear singlepath routing problems[END_REF]. The main difference is that we have to handle the fact that the volume of a traffic demand can change along its path. Note that property (c) follows from property (d) by taking d = 0. All proofs are omitted due to the lack of space but can be found in [START_REF] Amine | Offline Re-Optimization of VNF Placement Decisions for Existing Network Slices -A Game-Theoretic Algorithm[END_REF]. (c) It computes an optimal solution of problem (OPT) for linear resource costs, that is, when the functions ϕ r () are constant functions ϕ r (y) = a r .

(d) Its approximation ratio is 2

1 d+1 -1 -(d+1)
when the functions ϕ r () are polynomials of degree

d, that is, ϕ r (y) = d n=0 a r,n y n .
Algorithm 1 is optimal for linear resource costs and returns a solution whose cost is at most 5.83× the cost of an optimal solution for quadratic resource costs. Unfortunately, the approximation ratio of the algorithm for polynomial resource costs degrades quickly with the degree of the polynomial (it is already equals to 56.9 for d = 2). Our numerical results suggest however that the heuristic provides close-to-optimal solutions in most cases.

The proof of Theorem 1 shows that the game that we consider is a potential game, which implies that the worst-case complexity of Algorithm 1 is exponential in the number of traffic demands [START_REF] Fabrikant | The complexity of pure nash equilibria[END_REF]. However, as we will see in Section 6 below, its convergence in practice is much faster than suggested by this worst-case result.

Performance Evaluation

In this section, we experimentally evaluate the performance of the penalized best-response algorithm. We first describe the different types of objective functions used to solve the penalized best-response algorithm and the procedure for generating random instances in sections 6.1 and 6.2. Then, we present in section 6.3 the numerical results for the different topologies evaluated under the different types of objective functions.

Cost Functions

Given the cost per-unit of capacity of resource r function ϕ r , we denote by Φ r = y r ϕ r (y r ) the cost function of resource r. We consider four types of cost functions: linear, piece-wise linear, quadratic and the Kleinrock Function (M/M/1).

1. Linear cost function: the per-unit cost does not depend on the utilization rate of the resource r (ϕ r (y r ) = a, where a is a constant). Here, the objective function is linear with the amount of traffic which is solved using Gurobi 9.0 as the ILP solver. We note that in this case the solution may over utilize the resources capacity. We choose a = 5, and we have:

Φ r = 5 × y r (9)
2. Piece-Wise Linear cost function: we consider the following increasing per-unit cost function:

ϕ r (y r ) =      3 if 0 ≤ yr cr ≤ 1 4 , 5 -1 2 cr yr if 1 4 < yr cr ≤ 3 4 , 10 -17 4 cr yr if 3 4 < yr cr . (10)
where c r is the capacity of resource r. This function expresses that it is cheap to utilize a resource with a small utilization rate, whereas, as the load approaches c r , it becomes more expensive. Problem (OTP) can then be formulated as a Mixed-Integer Linear Problem and solved using Gurobi 9.0 as the MILP solver:

minimize r∈R Φ r
subject to:

y r = w∈W π∈Πw K k=0 δ r,k π b w k,r
x w,π , for all r ∈ R,

π∈Πw x w,π = 1, for all w ∈ W, (11)
Φ r ≥ ay r -bc r , for all r ∈ R and (a, b) ∈ C,

x w,π ∈ {0, 1}, for all π ∈ Π w and w ∈ W.

y r ≥ 0, Φ r ≥ 0, for all r ∈ R, (14)
where C = {(3, 0), (

Quadratic cost function: The per-unit cost function is linear (ϕ r (y r ) = ay r ). This results in a quadratic objective function which is solved using Gurobi 9.0 as the Quadratic Integer Program solver. Similar to the previous function, the more is the utilization rate of a resource, the more expensive it becomes.

Φ r = y 2 r /c 2 r (16)
4. Kleinrock cost function (M/M/1): with this function, it becomes hard to find the optimal solution. However, it guarantees that the utilization of a resource r does not exceed its capacity c r . Φ r = y r c r -y r

We formulate problem (OPT) as a Bilinear Program, and we solve it with Gurobi 9.0 Bilinear solver. minimize r∈R Φ r subject to:

y r = w∈W π∈Πw K k=0 δ r,k π b w k,r
x w,π , for all r ∈ R,

π∈Πw x w,π = 1, for all w ∈ W, (18)
y r -Φ r c r + Φ r y r = 0, for all r ∈ R,

x w,π ∈ {0, 1}, for all π ∈ Π w and w ∈ W.

y r ≥ 0, Φ r ≥ 0, for all r ∈ R, (21)

Generation of Random Instances

In order to evaluate the performance of the penalized best-response algorithm, we use five network topologies (see Table 2) collected from the IEEE literature and from The Internet Topology Zoo [START_REF] Knight | The internet topology zoo[END_REF]. As illustrated in Figure 3 for the NSFNET topology, we choose three source nodes (s 1 , s 2 and s 3 ) and three destination nodes (t 1 , t 2 and t 3 ) in each topology. The original topology is also expanded by adding three function nodes v 1 , v 2 and v 3 and connecting them to the other nodes. We consider three VNFs f 1 , f 2 and f 3 and assume that each one requires one core per unit of traffic. Function node v 1 can host VNFs f 1 and f 2 , whereas function node v 2 (resp. v 3 ) can host VNFs f 1 and f 3 (resp. f 2 and f 3 ). As shown in Figure 3, each VNF is thus available at two different locations.

We then generate 100 random instances for each network topology. A random instance is composed of 25 traffic demands, which are randomly generated as follows. The source and destination nodes of a traffic demand are chosen randomly in the sets {s 1 , s 2 , s 3 } and {t 1 , t 2 , t 3 }, respectively, according to uniform distributions. The volume of a traffic demand is drawn from a uniform distribution in the interval [START_REF] Nfv | Network Functions Virtualisation (NFV): Architectural framework[END_REF][START_REF] Amaldi | On the computational complexity of the virtual network embedding problem[END_REF], and, to simplify, we assume that it does not change along the path from the source to the destination. Finally, the processing path of a traffic demand in chosen uniformly at random in the set

{(f 1 , f 2 ), (f 1 , f 3 ), (f 2 , f 3 ), (f 1 , f 2 , f 3 )}.
Each traffic demand must thus go through two or three VNFs, which are available at multiple locations. The possible paths for a traffic demand w are generated as follows. We first compute two possible paths between the source node s w and each location where the first VNF f w 1 is available by solving a 2-shortest path problem (assuming unit weights for the edges of the expanded network). We then apply the same procedure for computing two possible paths between each possible location of the VNF f w 1 and each possible location of f w 2 , etc. The final set of candidate paths Π w is obtained by connecting the different segments forming the overall path from s w to t w .

The initial solution is obtained by choosing randomly a path for each traffic demand. We assume that all links and all function nodes have the same capacity. This capacity is computed for each random instance so as to obtain a network congestion rate of 0.83.

We report in Table 3 the number of candidate paths for a traffic demand |Π w | calculated for ) to choose a processing path size of two (three, resp.). In the case of CESNET topology and under all the generated instances, we could not find two candidate paths for each VNF option, hence the reduced size of Π w .

The final set of candidate paths calculated for each traffic demand is used to find the optimal solution with Gurobi solver. However, when running the penalized best-response algorithm, we consider only the k ′ -shortest paths from the previously obtained set. For instance, a traffic demand in the AARNET topology has on average 52 different paths that will be used to find the optimal solution. For k ′ = 10, only the 10 shortest paths will be considered in the penalized best-response algorithm.

Numerical Results

We report in Tables 4 and5 the performance of the penalized Best Response (BR) algorithm for different cost functions under distinct topologies. We set k ′ = 10, and we evaluate the algorithm performance by calculating the relative gap between the penalized BR solution and the optimal solution from Gurobi. Following the obtained results, the penalized BR algorithm provides very good quality solutions for the different cost functions with impressive computation time. Its average relative gap to the optimal solution is under 5%, and it is almost negligible for small topologies such as CESNET that has a result value of at most 0.07%. The average execution time of the algorithm depends on the size of the topology and on the cost function used. It is also directly proportional to the number of paths considered for each demand. Nevertheless except for the PWL cost function where the MILP solver is known to be efficient and quick at solving, the solver is significantly slower for the other more complex cost functions. In the case of the quadratic cost function, the maximum relative error of the BR algorithm over the 100 instances is 3.725%. However, if we look at the computing times in Table 5, we observe that the BR algorithm is at most 16 times faster and at worst 6.3 times faster than the solver. It is interesting to note that for the M/M/1 cost function, the solver has reached its time limit that has been set to 5 minutes, while the BR algorithm takes less than 10 seconds to execute with a solution as close as 1.43% to the solver. While PWL approximations can sometimes be considered for non-linear functions, it is not evident how to choose the intervals to build an approximation for non-linear cost functions.

The box plot in Figure 4 compares the distribution of the 100 instances generated for each scenario. The line that divides the box into two parts marks the median (Q2) of the data. The upper quartile (Q3) value is represented by the upper box line. 75% of the values fall below Q3. The lower quartile (Q1) value is the lower box line. 25% of the values fall below Q1. Hence, the values inside the box represent 50% of the data. The end of the lower and upper whiskers represents the minimum and maximum values, respectively. We first note that there are no outliers in the data samples which means that all the values of the data fall within the normal range. We observe that the values obtained with the quadratic cost function are more dispersed than the values obtained from the other cost functions.

In order to evaluate the choice of k ′ , we report in Table 6 the average relative gap and the computing times over different values of k ′ . We evaluate the results with the quadratic cost function (similar results were observed for the other cost functions). As expected, the optimality gap decreases when more routing paths are added to the set of candidate paths. However, the decrease of this gap is insignificant for k ′ > 10. In the case of the CESNET topology and since on average there are only 5 possible paths for a given demand (refer to Table 3), the gap optimality remains constant for k ′ ≥ 5. On the other hand, increasing the routing paths increases the running time of the BR algorithm. In this work, we found that the best trade-off between quality solution and running time is for k ′ = 10.

Conclusion

In this paper, we have considered the offline VNF placement and chaining problem, assuming that a non-linear cost function is associated to each network resources. With respect to previous works, the main originality of the considered model is that the cost per unit capacity of a resource is not constant, instead it grows with its utilization rate, which is an essential feature to achieve a better load distribution. We have formulated the problem as a single-path routing problem in an extended network and adapted an existing game-theoretic algorithm to solve it. Our numerical results suggest that the algorithm provides near-optimal solutions in substantially lower computing times than the solver, in particular for highly non-linear cost functions.

Acknowledgment

This work was supported by the French DGA project ONSET.

Figure 2 :

 2 Figure 2: Expanded network for the example of Figure 1. The artificial edges between consecutive layers are shown with dotted lines. Two feasible paths for the demand from node s 1 to node t 1 are shown, with blue-coloured edges for the first one and with red-coloured edges for the other one (common edges are purple-coloured).

Algorithm 1 n ← n + 1 7 :Theorem 1 .

 171 Penalized best-responseRequire: π (0)1: until π (n+1) = π (n) . 8: return π (n) Algorithm 1 has the following properties:(a) It converges in a finite number of steps.(b) If it ever reaches a global optimum of problem (OPT), it returns this optimal solution.

Figure 3 :

 3 Figure 3: NSFNET topology which has been expanded with function nodes v 1 , v 2 and v 3 . The nodes s 1 , s 2 and s 3 are the source nodes and the nodes t 1 , t 2 and t 3 are the destination nodes of the traffic demands.

Figure 4 :

 4 Figure 4: Box plot of the relative gap to the optimal solution with k ′ = 10.

Table 1 :

 1 Notations.Capacity required by demand w from resource r at stage k.

	G = (V, E)	Physical network, where V is the set of nodes and E is the set of links.
	D	Set of function nodes.
	R	Set of network resources (transmission links and function nodes).
	F (resp. F w ) Set of K network functions (resp. service chain for demand w).
	D(f )	Set of function nodes where VNF f can be executed.
	W	Set of traffic demands.
	s w (resp. t w ) Source (resp. destination) node of demand w.
	K w	Number of VNFs in the service chain F w .
	f w k b w k,r	k th VNF in the service chain F w .

Table 2 :

 2 Topologies : number of nodes and links.

	Topology	# Nodes # Links
	AARNET	19	46
	ARPANET	25	56
	NSFNET	13	30
	IBM	18	48
	CESNET	10	18
	s2	s1	t1
				f2, f3
	v2		t2	v3
	f1, f3
			s3
	t3		v1 f1, f2

Table 3 :

 3 Number of candidate paths Π w . The max (min, resp.) value is obtained when the size of the processing path is three (two, resp.) resulting in |Π w | = 27 (|Π w | = 2 5 , resp.). We remind the reader that we have k = 2 candidate paths for each VNF option. The average value is obtained over 100 problem instances where traffic demands have a probability of 3/4 (1/4, resp.

	Topology Max Min Average
	AARNET	128	32	52
	ARPANET	128	32	47
	NSFNET	128	32	48
	IBM	128	32	51
	CESNET	8	4	5
	each studied topology.

Table 4 :

 4 Relative gap (% ± std) to the optimal solution with k ′ = 10.

	Topology	PWL		Quadratic	M/M/1
	AARNET	0.71 ±0.65	1.81 ±0.85	1.43 ±1.76
	ARPANET	1.15 ±0.5		3.725 ±1.55	1.5 ±0.79
	NSFNET	0.42 ±0.49	1.67 ±1.12	0.31 ±0.32
		IBM	1.7 ±0.8		2.85 ±0.98	0.53 ±0.71
	CESNET	0.051 ±0.003 0.047 ±0.008 0.070 ±0.056
		8 10					PWL function Quadratic function M/M/1
	Relative gap (%)	4 6
		2
		0	AARNET	ARPANET	NSFNET	IBM	CESNET

Table 5 :

 5 Average computing times (seconds±std) for the 100 problem instances with k ′ = 10.

	Function Costs	AARNET BR Gurobi	ARPANET BR Gurobi	NSFNET BR Gurobi	BR	IBM Gurobi	CESNET BR Gurobi
	PWL	7.1±0	1.2	8.09±0	1.72	4.93±0	0.72	6.83±0	0.9	1.96±0	0.07
	Quadratic	7.89±0 126±0 8.89±0 124±0	6.3±0	40±0	8.27±0 109±0 2.33±0	0.3
	M/M/1	8.17±0	300	9.84±0	300	7.52±0	300	9.17		300	2.79±0	300

Table 6 :

 6 Average relative gap and computing times for the 100 problem instances with the quadratic cost function.

	Topology	gap	k ′ = 5 runtime	k ′ = 10 gap runtime	k ′ = 20 gap runtime	k ′ = 32 gap runtime
	AARNET	5.2%	3.8 sec	1.8%	7.9 sec	1.4%	17 sec	1.4% 26.7 sec
	ARPANET 5.5%	4 sec	3.7%	8.9 sec	3%	19 sec	3%	30 sec
	NSFNET	1.8%	2.9 sec	1.7%	6.3 sec	1.6%	12 sec	1.6%	18 sec
	IBM	5.88% 3.5 sec 2.85% 8.2 sec 2.77% 12.8 sec 2.75%	30 sec
	CESNET	0.05% 1.9 sec 0.05% 2.3 sec 0.05% 2.4 sec 0.05% 2.4 sec


Valentin Lauret
email: valentin.lauret@ens-lyon.fr

Hélène Labach

Matthieu Authier

Olivier Gimenez

Running head: Integrated single-visit occupancy models

Keywords: Bottlenose dolphins, Ecological monitoring, Integrated species distribution models, Multimethod, Occupancy models, Single-visit

A major challenge in statistical ecology consists of integrating knowledge from different datasets to produce robust ecological indicators. To estimate species distribution, occupancy models are a flexible framework that can accommodate several datasets obtained from different sampling methods. However, repeating visits at sampling sites is a prerequisite for using standard occupancy models. Occupancy models were recently developed to analyze detection/non-detection data collected during a single visit. To date, single-visit occupancy models have never been used to integrate several different datasets. Here, we showcase an approach that combines two datasets into an integrated single-visit occupancy model. As a case study, we estimated the distribution of common bottlenose dolphin (Tursiops truncatus) over the North-western Mediterranean Sea by combining 24,624 km of aerial surveys and 21,464 km of at-sea monitoring. We compared the outputs of single-vs. repeated-visit occupancy models into integrated occupancy models. Integrated models allowed a better sampling coverage of the targeted population, which provided a better precision for occupancy estimates than occupancy models using datasets in isolation. Overall, single-and repeated-visit integrated occupancy models produced similar inference about the distribution of bottlenose dolphins. We suggest that single-visit occupancy models open promising perspectives for the use of existing ecological datasets.

Introduction

In large-scale ecological analysis, several parallel monitoring programs are often carried out to collect ecological data [START_REF] Zipkin | Synthesizing multiple data types for biological conservation using integrated population models[END_REF]. Ecological monitoring programs are conducted by organizations operating across different time scales, geographic scales and funding initiatives [START_REF] Lindenmayer | The science and application of ecological monitoring[END_REF]. A major challenge is integrating knowledge from different monitoring programs to produce robust ecological indicators that may be used to inform decision-making [START_REF] Fletcher | A practical guide for combining data to model species distributions[END_REF][START_REF] Zipkin | Addressing data integration challenges to link ecological processes across scales[END_REF]. Recently, modelling tools have emerged to combine multiple data sources to estimate species distributions and Integrated models refer to the approaches that combine different data sources [START_REF] Miller | The recent past and promising future for data integration methods to estimate species' distributions[END_REF][START_REF] Isaac | Data Integration for Large-Scale Models of Species Distributions[END_REF]. The main purpose of integrated models is to improve the accuracy of ecological indicators [START_REF] Fletcher | A practical guide for combining data to model species distributions[END_REF][START_REF] Zipkin | Innovations in data integration for modeling populations[END_REF]. Species distributed over large areas could particularly benefit from integrated models because they allow for a global coverage of species occurrence by combining different data sources collected at different spatial scales [START_REF] Miller | The recent past and promising future for data integration methods to estimate species' distributions[END_REF]. To estimate species distribution in the face of uncertainties inherent to data collection, occupancy models are commonly used statistical tools [START_REF] Mackenzie | Estimating site occupancy rates when detection probabilities are less than one[END_REF]. Occupancy models have been developed to estimate species distribution while accounting for false negatives in the observation process [START_REF] Mackenzie | Estimating site occupancy rates when detection probabilities are less than one[END_REF]. Estimating occupancy when species detection is not perfect requires performing repeated visits to a set of sites to assess the detection probability [START_REF] Mackenzie | Occupancy estimation and modeling: inferring patterns and dynamics of species[END_REF]. However, repeating visits is sometimes unfeasible due to associated costs and logistical issues. In this context, two relevant developments of occupancy models have been recently proposed. First, integrated occupancy models combine data from different monitoring programs to improve the estimation of species distribution [START_REF] Nichols | Multi-scale occupancy estimation and modelling using multiple detection methods[END_REF][START_REF] Miller | The recent past and promising future for data integration methods to estimate species' distributions[END_REF][START_REF] Fletcher | A practical guide for combining data to model species distributions[END_REF]. Second, [START_REF] Lele | Dealing with detection error in site occupancy surveys: what can we do with a single survey[END_REF] used occupancy models to estimate species distribution and detectability while having only one visit at the sampling site, i.e. hereafter single-visit occupancy models. An increasing number of studies suggest that under certain conditions, single-visit models produce robust estimates of occupancy without repeating visits at the sampling sites [START_REF] Lele | Dealing with detection error in site occupancy surveys: what can we do with a single survey[END_REF][START_REF] Sólymos | Revisiting resource selection probability functions and single-visit methods: clarification and extensions[END_REF][START_REF] Peach | Single-visit dynamic occupancy models: an approach to account for imperfect detection with Atlas data[END_REF]. Besides, single-visit occupancy offers the possibility to work with existing datasets that did not carry out repeated visits, which is relevant to population biology and management. In this paper, we develop an integrated approach that combines two single-visit occupancy models and demonstrate that combining several datasets into integrated single-visit occupancy models leads to accurate ecological parameter estimation. We also investigate the performance of single-visit vs. repeated-visit occupancy models. As a case study, we focused on the distribution of Bottlenose dolphins (Tursiops truncatus) in the North-Western Mediterranean Sea. We combined aerial surveys and at-sea monitoring into integrated occupancy models and we compared the outputs of integrated occupancy models to occupancy models using each dataset in isolation. Last, we discuss the advantages of integrated singlevisit occupancy models to deal with existing ecological monitoring programs.

Methods

MODEL DESCRIPTION

Latent ecological process

Occupancy models estimate spatial distribution while accounting for imperfect species detection [START_REF] Mackenzie | Estimating site occupancy rates when detection probabilities are less than one[END_REF]. The formulation of occupancy models as state-space models allows distinguishing the latent ecological state process (i.e. species present or absent at a gridcell) from the detection process [START_REF] Royle | A Bayesian state-space formulation of dynamic occupancy models[END_REF]. We denote zi the latent occupancy of grid-cell i (z = 1, presence; z = 0, absence). We assume zi is drawn from a Bernoulli distribution with Ψi the probability that the species is present at grid-cell i:

zi ~ Bernoulli(Ψi)
We modelled Ψ as a function of some environmental covariate on a logit scale, say habitat. logit(Ψi) = b0 + b1 habitati where parameters b0, and b1 are to be estimated.

Repeated-visit observation process

In standard occupancy designs, each grid-cell is visited J times to estimate the detection probability. We denote yi,j (yi,j = 0, no detection ; yi,j = 1, detection) the observations corresponding to the data collected at grid-cell i during visit j (j =1,..,J). Repeating visits at a grid-cell allows estimating species detectability, with pi,j being the probability of detecting the species at visit j given it is present at grid-cell i: yi,j | zi ~ Bernoulli(zi pi,j)

Single-visit observation process

The difference with repeated-visit occupancy models lies in the number of sampling occasions which is J = 1 in single-visit occupancy models. The j subscript is dropped and we denoted yi the observation corresponding to the data collected at site i. Subsequently, pi is the probability of detecting the species during the single visit given it is present at site i: yi | zi ~ Bernoulli(zi pi). Single-visit occupancy models require certain conditions to be fulfilled for estimating detection probabilities reliably. First, different continuous covariates should be used to estimate detection and occupancy probabilities [START_REF] Lele | Dealing with detection error in site occupancy surveys: what can we do with a single survey[END_REF][START_REF] Peach | Single-visit dynamic occupancy models: an approach to account for imperfect detection with Atlas data[END_REF]. Second, the number of detections may affect the estimation of occupancy in the case of rare or ubiquitous species [START_REF] Peach | Single-visit dynamic occupancy models: an approach to account for imperfect detection with Atlas data[END_REF]. Third, the use of inappropriate link functions to model the detection process may lead to model misspecification and biased interpretation (e.g. loglink and scaled logit link function on detection, [START_REF] Knape | Estimates from non-replicated population surveys rely on critical assumptions[END_REF]. However, most often, the logit link function is used for detection, which makes the single-visit approach valid [START_REF] Sólymos | Revisiting resource selection probability functions and single-visit methods: clarification and extensions[END_REF]. Despite these concerns, simulation studies have showed that situations where single-visit occupancy models fail are rare (Sólymos andLele 2016, Peach et al. 2017) and, in practice, the conditions for a valid application of single-visit occupancy models are often fulfilled [START_REF] Sólymos | Revisiting resource selection probability functions and single-visit methods: clarification and extensions[END_REF]. We detailed the modeling assumptions of singlevisit occupancy models in Appendix S4. Because the number of detections is an important condition to accurately estimate single-visit occupancy parameters [START_REF] Peach | Single-visit dynamic occupancy models: an approach to account for imperfect detection with Atlas data[END_REF], we expect that integrated approaches will be beneficial to single-visit occupancy modelling by increasing the number of detections (true occupancy) available.

Integrated occupancy models

We developed an integrated occupancy model using data from two independent monitoring programs, say A and B. The state process driving the latent occupancy state of site i, zi, remains unchanged and is drawn from a Bernoulli distribution with probability y, which is modeled as a function of environmental covariates. The observation of the targeted species at site i during occasion j may take four values with yi,j = 0 for no detection, yi,j = 1 for detection in dataset A, yi,j = 2 for detection in dataset B, and yi,j = 3 for detection in both datasets A and B. For convenience, we drop the subscripts in the notation as the formulation of the integrated observation process is identical whether we consider single-visit occupancy (i.e. J = 1) or repeated-visit occupancy (J > 1). Assuming that detection methods are independent, the observation process can be written using detection probability by the monitoring program A (pA) and detection probability by the monitoring program B (pB):

𝑦|𝑧~𝑀𝑢𝑙𝑡𝑖𝑛𝑜𝑚𝑖𝑎𝑙(1, 𝑧𝜋) with 𝜋 = [𝑝 ! , 𝑝 " , 𝑝 # , 𝑝 $ ] = [𝑝𝑟(𝑦 = 0), 𝑝𝑟(𝑦 = 1), 𝑝𝑟(𝑦 = 2), 𝑝𝑟(𝑦 = 3)] 𝜋 = [1 -𝑝 % -𝑝 & + 𝑝 % 𝑝 & , 𝑝 % (1 -𝑝 & ), 𝑝 & (1 -𝑝 % ), 𝑝 % 𝑝 & ]
We modeled monitoring-specific detection probabilities as functions of the sampling effort of each monitoring program:

logit(pA) = a0A + a1A log(SeffA) logit(pB) = a0B + a1B log(SeffB).
where the parameters a0A, a1A, a0B, and a1B are to be estimated. For example, if we assume that the detection history at site i is yi = {2,0,1,2} over J = 4 sampling occasions, i.e. the species was detected by monitoring program B only at sampling occasions j = 1 and j = 4, then went undetected at j = 2, and was detected by monitoring program A only at j = 3, then for singlevisit integrated occupancy we consider yi = {3} because both monitoring programs detected the species at site i. We ran a simulation study comparing the performance of single-vs. repeated-visit occupancy over different scenarios affecting occupancy, and detection probabilities (Appendix S1).

BOTTLENOSE DOLPHINS CASE STUDY

We aimed at estimating bottlenose dolphin (Tursiops truncatus) distribution in an area of 255,000 km 2 covering the North-Western Mediterranean. The protected status of this species within the French seas led to the development of specific programs to monitor Mediterranean bottlenose dolphins within the implementation of the European Marine Strategy Framework Directive (2008/56/EC; MSFD), which involve estimating common bottlenose dolphin distribution. We considered two large-scale monitoring programs about bottlenose dolphins.

We divided the study area in 4,356 contiguous pixel/grid-cells creating a 5'x5' Mardsen grid (WGS 84) that we used for all the occupancy models we considered. We used data from at-sea surveys over 21,464 km of the French continental shelf (456 grid-cells sampled, 10.46% of the total number of grid-cells). Observers performed monitoring aboard small sailing and motor boats to locate and photo-identify bottlenose dolphins all year long between 2013 and 2015 [START_REF] Labach | Distribution and abundance of bottlenose dolphin over the French Mediterranean continental shelf[END_REF]. At-sea surveys detected 129 distinct bottlenose dolphin groups located in 89 different grid-cells. At-sea surveys did not include planned repeated visits, some grid-cells have been visited once, and others have been visited 50 times. Then, using repeated-visits occupancy models to analyze the at-sea monitoring requires considering only the grid-cells sampled multiple times and hence to drop the data collected in grid-cells sampled only once.

Single-visit models enable us to include all data, even data collected in grid-cells that were surveyed only once, which make at-sea a relevant candidate for single-visit model implementation. Besides, we considered data collected during aerial line-transects covering 24,624 km of the French Exclusive Economic Zone (EEZ), targeting marine megafauna, and following a distance sampling protocol. The survey sampled 1336 grid-cells (i.e. 30.67% of the total number of grid-cells). Aerial surveys produced 130 distinct bottlenose dolphin detections located in 87 grid-cells. Sampling effort for aerial surveys was homogeneous over the study area with three or four replicates per line-transect between November 2011 and August 2012 [START_REF] Laran | Seasonal distribution and abundance of cetaceans within French waters-Part I: The North-Western Mediterranean, including the Pelagos sanctuary[END_REF]. Because we used occupancy models, we only considered detection/nodetection data, which lead to a binary 0/1 dataset. Hence, multiple sightings detected in the same groups were coded as 1. Thus, we obtain the two aerial and at-sea detection/no-detection datasets that we analyzed with occupancy models. An important assumption of single-season occupancy models is that the latent ecological state of a grid-cell (the zi's) remains unchanged between the repeated visits (MacKenzie 2006). When monitoring highly mobile species, such as cetaceans, the closure assumption is likely to be violated because individuals can move into and out of the sampling grid-cell. The size of the grid-cells is much lower than dolphins' range of activity. If individuals' movement in and out of the sampling units is random, then the occupancy estimator is unbiased [START_REF] Kendall | Relaxing the closure assumption in occupancy models: staggered arrival and departure times[END_REF]). However, it is unlikely the case for bottlenose dolphins because their use of space is driven by ecological and environmental factors, and occupied locations are used only temporarily by individuals (MacKenzie 2006; [START_REF] Neilson | Animal movement affects interpretation of occupancy models from camera-trap surveys of unmarked animals[END_REF]. Closure assumption is crucial to the interpretation of occupancy model's parameters. In cases where this assumption is known to be violated, the parameter is usually interpreted as the probability that a location is used by the species as opposed to probability of species presence. In this situation, the occupancy estimator Ψi represents the probability that grid-cell i is used by the target species [START_REF] Kendall | Relaxing the closure assumption in occupancy models: staggered arrival and departure times[END_REF], being interpreted as space-use by bottlenose dolphins. Occupancy and space-use refer to distinct ecological concepts. Occupancy describes the species home range that can be defined as the geographic range of occurrence, while space-use refers to the usage made by individuals of the different components of the home range (e.g. feeding locations, migratory routes, [START_REF] Johnson | The Comparison of Usage and Availability Measurements for Evaluating Resource Preference[END_REF]. Then, both single-visit and repeated-visits occupancy models infer the probability that a particular grid-cell is used by the species. The detection probability now accounts for both the probability of detecting the species given that the species is available for sampling, and the probability that the species is using the grid-cell during the sampling, reflecting that the species might occupy the grid-cell but not during the sampling occasion (MacKenzie 2006). As stated above, single-visit occupancy relaxed the closure assumption because the inference of the detection probability does not require site closure between the repeated visits. However, the interpretation of the occupancy parameter is still space-use in the case of bottlenose dolphins because our single-visit data is collected during multiple years and dolphin are expected to move in and out the sampling unit area during the sampling period.

Because at-sea and aerial surveys were performed during different years, we considered them as independent. In 2018, recent Mediterranean scale aerial monitoring program sampled French Mediterranean following the same line-transect protocol as the aerial dataset we analyzed (ACCOBAMS Survey Initative, 2018). Preliminary and unpublished results from the 2018 program estimated similar common bottlenose dolphin distribution to that of 2011-2012. Then, we assumed that space-use remained unchanged during the monitoring period (i.e. 2011 to 2015). Besides, we neglected the seasonal variation in the bottlenose dolphin space-use in this case study. Concerning the ecological process, we used two environmental covariates to estimate the space-use of bottlenose dolphins: i) bathymetry, which is expected to have a positive effect on bottlenose dolphins' occurrence [START_REF] Bearzi | Ecology and conservation of common bottlenose dolphins Tursiops truncatus in the Mediterranean Sea[END_REF][START_REF] Labach | Distribution and abundance of bottlenose dolphin over the French Mediterranean continental shelf[END_REF], and ii) sea surface temperature (SST, AQUA MODIS | NASA 2019, https://neo.sci.gsfc.nasa.gov/), which is locally related to dolphins' prey abundance and hence expected to affect local distribution of bottlenose dolphins [START_REF] Bearzi | Ecology and conservation of common bottlenose dolphins Tursiops truncatus in the Mediterranean Sea[END_REF]). We extracted average SST between 2011 and 2015 value in each grid-cell, making SST a cell-specific covariate. Similarly, bathymetry had a single value per grid-cell. We checked for correlation between the two covariates and the Pearson coefficient was < 0.3. Then, we modelled Ψ as a function of bathymetry, SST, and the interaction between bathymetry and SST on a logit scale:

logit(Ψi) = b0 + b1 bathymetryi + b2 SSTi + b3 bathymetryi SSTi
Regarding the observation process, we calculated the transect length (in km) prospected by each monitoring protocol within each grid-cell during a time period. Sampling effort was therefore a grid-cell-specific and time-specific covariate; SeffA refers to the sampling effort of the aerial monitoring program while SeffS refers to the sampling effort of the at-sea monitoring program.

We modeled monitoring-specific detection probabilities as functions of the relevant sampling effort: logit(pa) = a0 + a1 log(SeffA) logit(ps) = a0 + a1 log(SeffS).

Regarding the repeated-visit occupancy models, we divided the detection/non-detection datasets into four sampling occasions (J = 4): winter (January, February, March), spring (April, May, June), summer (July, August, September), autumn (October, November, December). For the single-visit occupancy models, we considered the entire monitoring program in a single occasion. For example, let us assume that the detection history at site i is yi= {0,1,1,0} in repeated-visit occupancy, i.e. the species was detected at sampling occasions j = 2 and j = 3, and went undetected at j = 1, and j = 4, then for single-visit occupancy we have yi= {1}. In addition, the single-visit sampling effort in a grid-cell was equal to the sum of the sampling effort over the 4 sampling occasions of the repeated-visit occupancy model.

Performances of integrated models

To assess the added value of combining aerial and at-sea datasets into integrated single-visit occupancy models, we analyzed 3 datasets: i) the aerial dataset, ii) the at-sea dataset, and iii) the two datasets together into an integrated occupancy model. For each of these datasets, we applied repeated-visit and single-visit occupancy models. Besides the case study, we also carried out a simulation study to test for the performances of integrated occupancy models (Appendix S2). In Appendix S5, we go through a worked example of the likelihood functions for single-visit, repeated-visit, integrated repeated-visit, and integrated single-visit occupancy models. In Appendix S4, we listed the modeling assumptions required to run the different occupancy models.

Bayesian implementation

We ran all models with three Markov Chain Monte Carlo chains with 100,000 iterations each in JAGS [START_REF] Plummer | JAGS: A program for analysis of Bayesian graphical models using Gibbs sampling[END_REF] called from R (R Core Team, v 3.2.5 2019) using the r2jags package [START_REF] Su | R2jags: Using R to Run[END_REF]. We checked for convergence calculating the R-hat parameter [START_REF] Gelman | Bayesian data analysis[END_REF]) and reported posterior means and 95% credible intervals (CI) for each regression coefficient of covariates affecting space-use probability (Fig. 1). Hereafter, we considered effect size of a covariate as the estimate of its regression coefficient. We discussed the effect of a covariate whenever the 95% CI of its associated parameter did not include 0. From covariates' effect size, we calculated the predicted space-use by bottlenose dolphins (i.e. Ψ, Fig. 2). We reported maps of standard deviation of Ψ (Fig. 2B). On the maps, we displayed mean and standard deviation of Ψ for coastal and pelagic seas according to a 500m deep boundary that corresponds to the separation of continental shelf from the abysses. Data and codes are available on

Results

All models produced similar predictions of space used by bottlenose dolphins (Fig. 2).

The 95% CI of SST, and of the interaction between SST and bathymetry included 0 in all models (Fig. 1). The probability of space-use increased with decreasing bathymetry for all models (Fig. 1). Bathymetry ranges from altitude of 0 m to -3,488 m deep, hence a positive influence of bathymetry referred to a preference for a high seafloor (e.g. 0-200m depth).

Overall, maps showed greater probabilities of space-use on the continental shelf (mean Ψ = 0.76 SD ± 0.17) than on the high seas (mean Ψ = 0.40 SD ± 0.15), although magnitudes of Ψ were different between models (Fig. 2). Bathymetry posterior means were highest for at-sea occupancy (although the 95% CI of effect size included 0), which resulted in models using only at-sea survey data predicting the highest contrast between the continental shelf and the highseas. Bathymetry effect size was the lowest for aerial occupancy while maps from integrated occupancy models displayed moderate contrast of space-use between shelf and pelagic waters (Fig. 2). Single-visit occupancy models exhibited similar covariates estimates to those of repeated-visit occupancy models (Fig. 1). For aerial occupancy, we noticed similar space-use prediction between single-and repeated-visit (Fig. 2A). For at-sea, predicted space-use probabilities were different between single-visit and repeated-visit occupancy models (Fig. 2).

When considering the covariates' effect size (Fig. 1), the widths of the 95% CI were not smaller for integrated occupancy than for occupancy models using datasets in isolation.

However, when looking at the standard deviation of the predicted probability of space-use, integrated occupancy models had a better precision than aerial or at-sea occupancy models separately, (Fig. 2B). The use of integrated single-visit occupancy models also improved precision in predicted space-use compared to single-visit occupancy built from aerial and atsea datasets separately (Fig. 2B). Inspecting the simulation results, we found that 1) integrated occupancy models produced more precise estimates of covariates effect size than occupancy models fitted to a single dataset (Appendix S2), and 2) single-visit occupancy models produced similar results to repeated-visit occupancy models (Appendix S1).

Discussion

Integrated single-visit occupancy models provide reliable ecological inference

Ecological inference from integrated occupancy models lied within the range of the estimates obtained with each dataset separately (Fig. 1). Across all occupancy models, the effects of environmental covariates were similar and consistent with previous studies.

Bottlenose dolphins were more likely to use shallower seas [START_REF] Bearzi | Ecology and conservation of common bottlenose dolphins Tursiops truncatus in the Mediterranean Sea[END_REF][START_REF] Labach | Distribution and abundance of bottlenose dolphin over the French Mediterranean continental shelf[END_REF], and depth had a stronger effect than SST on the use of space by bottlenose dolphins [START_REF] Torres | Fine-scale habitat modeling of a top marine predator: do prey data improve predictive capacity[END_REF]). However, we found variations among models in the estimation of the probability of space-use by dolphins (Fig. 1). In particular, at-sea occupancy models predicted that dolphins make little use of the pelagic seas compared to the continental shelf, while aerial occupancy models predict more homogeneous space-use between coasts and pelagic seas.

Aerial surveys detected several dolphin groups in the high depths while at-sea surveys detected none. Detecting offshores groups tempered the preference for low-depth seafloors in aerial occupancy models (Appendix S6). Besides, we recommend caution in interpreting predicted maps of space-use as predicted space-use was sensitive to the mean value of covariate effect size. Therefore, depth being the only covariate that affect space-use probability, maps of predicted space-use were mostly driven by bathymetry effect size, and did not account for precision associated with space-use prediction. Because depth posterior mean was similar between occupancy models, differences between predicted space-use maps do not provide a relevant illustration to compare occupancy models performances, nor they reflect the uncertainty associated with the occupancy models' estimates. To study the benefits of singlevisit and integrated occupancy models to accommodate existing ecological datasets, we emphasize standard deviation maps and the credible intervals of covariates effect size (Fig. 1-2B). Integrated occupancy models had a better precision in space-use than models using aerial or at-sea surveys separately (Fig. 2). This result was supported by our simulation study which demonstrates the better performance of integrated occupancy models at estimating covariate effect size compared to occupancy models from a single dataset (Appendix S2). Single-visit occupancy models gave similar estimates to those obtained with repeated-visit occupancy models, although repeated-visit occupancy models exhibited better precision (Fig. 1-2B), as well as in our simulations (Appendix S1). In the bottlenose dolphins case study, we considered two existing monitoring programs that were not initially designed for occupancy modeling. In the at-sea monitoring, repeated line-transects were not implemented, nor the high depths were sampled, which made at-sea occupancy unlikely to exhibit precise estimates at our spatial extent. The two datasets exhibit complementary characteristics. While aerial surveys covered a larger spatial scale than at-sea surveys, at-sea surveys exhibited a better detection rate.

Detection probability was greater for at-sea surveys (p = 0.18 SD ± 0.04) than for aerial surveys (p = 0.10 SD ± 0.03). Regarding the aerial dataset, the number of occurrences was low despite the important coverage of the monitoring design (i.e. bottlenose dolphins were detected in 6.5%

of sampled grid-cells), which might hinder the implementation of single-visit occupancy when the number of occurrences is less than 10% of the sampling units [START_REF] Peach | Single-visit dynamic occupancy models: an approach to account for imperfect detection with Atlas data[END_REF]. However, the at-sea dataset had occurrences in 19.5% of sampled units. Using integrated occupancy models enables to combine low-frequency occurrence data like the aerial dataset with another dataset to increase the amount of information about the ecological state process and helps mitigating the issue of low number of occurrences.

Ecological implications and perspectives

Overall, we illustrate that: i) Integrating datasets into occupancy models improves the precision of space-use estimates, and ii) Single-visit occupancy models can reliably accommodate the lack of repeated visits that occurs frequently. Integrated occupancy models produced more reliable estimates than occupancy models using datasets in isolation in both the bottlenose dolphin data analyzes and the simulations. Our finding on the bottlenose dolphins case study is a good illustration of the well-known benefit of combining datasets into integrated species distribution models to increase precision in ecological inference [START_REF] Fletcher | A practical guide for combining data to model species distributions[END_REF].

Some advanced developments of occupancy models allow combining datasets to estimate occupancy parameters at multiple spatial scales [START_REF] Nichols | Multi-scale occupancy estimation and modelling using multiple detection methods[END_REF][START_REF] Pavlacky | Hierarchical multi-scale occupancy estimation for monitoring wildlife populations[END_REF].

Besides, integrated occupancy modeling has also been used to evaluate ecological monitoring programs prior to their implementation (e.g., comparing capabilities of different detection devices, Otto & Roloff 2011;Haynes et al. 2013). Here, we emphasize the benefit of considering integrated methods combined with single-visit occupancy modeling after data collection. When the species of interest either occurs over a large spatial scale or is a highly mobile species (such as bottlenose dolphins), considering multiple sampling methods is effective to monitor the entire population making the most of each device [START_REF] Zipkin | Synthesizing multiple data types for biological conservation using integrated population models[END_REF]. In particular, integrating a large volume of data, such as those that can be leveraged through citizen-science programs or with dedicated NGOs over the years can make the most of ecological monitoring programs for the furthering of many applied situations [START_REF] Zipkin | Innovations in data integration for modeling populations[END_REF]). However, caution should be taken as integrating data is not always beneficial and requires additional modelling assumptions according to the particularity of each dataset ot include (Dupont et al., 2019;Farr et al., 2020;[START_REF] Fletcher | A practical guide for combining data to model species distributions[END_REF][START_REF] Simmonds | Is more data always better? A simulation study of benefits and limitations of integrated distribution models[END_REF].

Although repeated-visit occupancy models remain statistically more precise, there are benefits in using single-visit occupancy models. The ability of single-visit occupancy to relax the closure assumption is appealing, because this assumption is often incompatible with the behavior of mobile species and for numerous monitoring programs of animal populations [START_REF] Rota | Occupancy estimation and the closure assumption[END_REF][START_REF] Issaris | Occupancy estimation of marine species: dealing with imperfect detectability[END_REF][START_REF] Sólymos | Revisiting resource selection probability functions and single-visit methods: clarification and extensions[END_REF][START_REF] Lele | Dealing with detection error in site occupancy surveys: what can we do with a single survey[END_REF][START_REF] Kendall | Relaxing the closure assumption in occupancy models: staggered arrival and departure times[END_REF].

In this study, the closure assumption is unlikely to be valid for bottlenose dolphins over the time span of the two monitoring programs, because dolphins obviously would not remain into the same grid-cell. Besides, when financial or logistical costs are limited, implementing a single-visit monitoring design could provide robust ecological inference while explicitly accounting for imperfect species detection [START_REF] Lele | Dealing with detection error in site occupancy surveys: what can we do with a single survey[END_REF][START_REF] Dénes | Biome-scale signatures of land-use change on raptor abundance: insights from single-visit detectionbased models[END_REF]. Overall, increasing quantity and types of biodiversity data are becoming available [START_REF] Isaac | Data Integration for Large-Scale Models of Species Distributions[END_REF].

Numerous monitoring programs do not rely on protocols implementing repeated visits like, e.g., historical monitoring programs, or citizen science programs (Tingley andBeissinger 2009, Zipkin and[START_REF] Zipkin | Synthesizing multiple data types for biological conservation using integrated population models[END_REF]. Then, using single-visit occupancy models helps making efficient use of available data, which is of great interest in many ecological applications [START_REF] Nichols | Monitoring for conservation[END_REF]Williams 2006, Sólymos and[START_REF] Sólymos | Revisiting resource selection probability functions and single-visit methods: clarification and extensions[END_REF]. In this context, [START_REF] Miller | The recent past and promising future for data integration methods to estimate species' distributions[END_REF] encouraged further developments of methods mixing standardized and non-standardized datasets. To illustrate, we built an integrated occupancy model mixing repeated-visit occupancy models for aerial surveys and single-visit occupancy models for at-sea surveys (Appendix S3). One could also extend integrated occupancy models to more than two datasets. However, caution should be taken when integrating datasets, as combining different sources of information does not always outperform the analysis of single datasets in isolation (Lele andAllen 2006, Simmonds et al. 2020). The flexibility of occupancy models provided a relevant framework to combine monitoring programs and to accommodate different types of data collection. Integrated and single-visit occupancy models contribute to widen the scope of possibilities. We emphasize the usefulness of both integrated and single-visit approaches to deal with existing datasets. We anticipate that their combination into integrated single-visit approaches will be of most interest for many parties in ecological research.

Figure captions

Figure 1: Effect size of bathymetry, sea surface temperature (SST), and interaction between SST and bathymetry on the probability Ψ that a site is used by Bottlenose dolphins (Tursiops truncatus). The posterior mean is provided with the associated 95% credible interval. "SV" refers to single-visit occupancy models, "RV" to repeated visit occupancy models, and "IOM" stands for integrated occupancy models, in which aerial surveys and at-sea surveys are combined. Estimates are given on the logit scale.

Figure 2 :

 2 Figure 2: A. Probability of predicted space-use by Bottlenose dolphins (Tursiops

  Figure 1

Acknowledgements

The French Ministry in charge of the environment (Ministère de la Transition Energetique et Solidaire) and the French Office for Biodiversity (OFB) funded the project SAMM. The PELAGIS observatory, with the help of the OFB, designed, coordinated and conducted the survey. We are grateful to all financial partners of the GDEGeM project. We warmly thank technical and scientific participants of GDEGeM. We also thank three anonymous reviewers for their insightful comments and suggestions.

Literature cited

ACCOBAMS Survey Initiative. 2018. ASI 2018 Summer Survey | Accobams. https://accobams.org/main-activites/accobams-survey-initiative-2/asi-preliminaryresults/.

data will be stored on the following permanent repository https://github.com/valentinlauret/IntegratedSingleVisitOccupancy

Data S1, and on GitHub at https://github.com/valentinlauret/IntegratedSingleVisitOccupancy.