Amin Gheibi
email: [agheibi@scs.carleton.ca

Anil Maheshwari

Jörg-Rüdiger Sack
email: sack]@scs.carleton.ca

Minimizing Walking Length in Map Matching

In this paper, we propose a geometric algorithm for a map matching problem. More specifically, we are given a planar graph, H, with a straight-line embedding in a plane, a directed polygonal curve, T , and a distance value ε > 0. The task is to find a path, P , in H, and a parameterization of T , that minimize the sum of the length of walks on T and P whereby the distance between the entities moving along P and T is at most ε, at any time during the walks. It is allowed to walk forwards and backwards on T and edges of H. We propose an algorithm with O (mn (m + n) log(mn)) time complexity and O (mn (m + n)) space complexity, where m (n, respectively) is the number of edges of H (of T , respectively). As we show, the algorithm can be generalized to work also for weighted non-planar graphs within the same time and space complexities.

Introduction

Trajectory data are often obtained from global positioning system (GPS) devices. Such devices have accuracy limitations due to noise, sampling intervals, or poor signals (e.g., inside buildings) thus raw spatial trajectories tend not to be accurate. Under the assumption that the travel captured by the trajectory was following edges of a map (stored as a graph) the map matching problem arises. It asks to find a path on the map that "corresponds well" to the given trajectory. Map matching arises in different contexts and is a necessary step in preprocessing raw data before data mining [START_REF] Zheng | Trajectory Data Mining: An Overview[END_REF]. A variety of approaches have been used to solve the map matching problem (e.g. geometric, probabilistic methods, fuzzy logic, neural networks). In [START_REF] Chen | Map-matching algorithm for large-scale low-frequency floating car data[END_REF], Chen et al. discussed recent map matching algorithms when a trajectory is obtained from low-frequency GPS data of vehicles driving on a road network. Ruan et al. [START_REF] Ruan | A Method of Map Matching in Indoor Positioning[END_REF] studied indoor map matching technology based on personal motion states. In [START_REF] Asakura | A Pedestrian-oriented Map Matching Algorithm for Map Information Sharing Systems in Disaster Areas[END_REF], Asakura et al. proposed a pedestrian-oriented map matching algorithm in the context of disasters. In this context, refugees have battery-driven mobile GPS terminals and move to shelters at walking speed. They stated that in order to reduce battery consumption (which is vital in this context), they chose a geometric approach in which computation resources are less utilized when compared e.g., with probabilistic methods.

In this paper, we focus on geometric approaches. We assume that a map is given as a planar graph via a straight-line embedding in a plane. Therefore, a path in the graph corresponds to a polygonal curve in the plane. A trajectory is given as a directed polygonal curve from a starting point to an ending point. The objective is to find a path in the map which is most similar to the given trajectory.

To measure similarity, [START_REF] Zheng | Trajectory Data Mining: An Overview[END_REF] and [START_REF] Brakatsoulas | On map-matching vehicle tracking data[END_REF] observe that methods which consider global features of the input trajectories achieve more accurate results than local approaches. The Fréchet distance is a global similarity measure between curves, see e.g., the seminal paper [START_REF] Alt | Computing the Fréchet distance between two polygonal curves[END_REF]. Commonly, the Fréchet distance is illustrated as follows: Suppose a person wants to walk along one curve and his/her dog on another; the person is keeping the dog at a leash. Both person and dog walk, from starting point to ending points along their respective curves. The standard Fréchet distance is the minimum leash length required without either person or dog needing to backtrack. The weak Fréchet distance is a variant of the standard Fréchet distance in which backtracking on one or both curves is allowed. Alt and Godau [START_REF] Alt | Computing the Fréchet distance between two polygonal curves[END_REF] proposed algorithms to compute the standard and weak Fréchet distances in O(n 2 log n) time, where n is the maximum number of segments in the input polygonal curves. Har-Peled and Raichel [START_REF] Har-Peled | The Fréchet Distance Revisited and Extended[END_REF] showed that the weak Fréchet distance can be computed in quadratic time.

In [START_REF] Alt | Matching planar maps[END_REF], Alt et al. discussed the map matching problem set in the context of the standard Fréchet distance. I.e., their algorithm finds a path in the planar graph with minimum Fréchet distance to the given trajectory. The time complexity of their algorithm is O(mn log 2 mn) where m (n, respectively) is the number of edges in the input planar graph (the input polygonal curve, respectively). Brakatsoulas et al. [START_REF] Brakatsoulas | On map-matching vehicle tracking data[END_REF], extended the map matching algorithm of [START_REF] Alt | Matching planar maps[END_REF] for the weak Fréchet distance. The time complexity of their algorithm is O(mn log mn). In [START_REF] Chen | Approximate Map Matching with respect to the Frchet Distance[END_REF], Chen et al. proposed a (1 + ε)-approximation algorithm for the map matching problem when the similarity measure is the standard Fréchet distance and input model is more "realistic". They assumed that the input polygonal curve is cpacked and the input graph is φ-low density in R d (see Section 2 of [START_REF] Chen | Approximate Map Matching with respect to the Frchet Distance[END_REF]).

In [START_REF] Gheibi | Minimum Backward Fréchet Distance[END_REF], Gheibi et al. studied a natural optimization problem on the weak Fréchet distance, called the minimum backward Fréchet distance (MBFD) problem. There, the task is to determine a pair of walks for a given input leash length such that the total length of backtracking on both input polygonal curves is minimized. The cost of backtracking could represent, for example, the cost of moving against a flow, or the cost for a moving entity (e.g., a human, a humanoid robot) to move backwards because of the entity's physiology [START_REF] Flynn | Comparison of cardiopulmonary responses to forward and backward walking and running[END_REF]. They proposed an algorithm solving this problem within time complexity O(n 2 log n) and space complexity O(n 2 ), where n is the maximum number of segments in the polygonal curves. In [START_REF] Gheibi | Weighted Minimum Backward Fréchet Distance[END_REF], the weighted variant of the MBFD problem is solved in O(n 3 ) time. In this variant, each edge of the input polygonal curves has an associated non-negative weight to capture different costs for backward movement.

In this paper, we study the map matching problem when the similarity measure is the MBFD. More specifically, as input, we are given: a planar graph, H, with a straight-line embedding in a plane, a directed polygonal curve, T , and a distance ε > 0. As motions, both forward and backward motions along T and the edges of H are allowed. The objective is to find a path, P , in H, and a parameterization of T , that minimize the sum of the walk lengths along T and P while keeping a leash length of at most ε. We restrict the start and end point of P to be at a vertex of H. However, P may partially contain an edge of H. The difference between this problem setting and that optimization setting of [START_REF] Alt | Matching planar maps[END_REF] and [START_REF] Brakatsoulas | On map-matching vehicle tracking data[END_REF], is that here the total walking length along T and P is minimized while in the other settings the leash length is minimized. The optimization problem introduced in this paper, can also be used to track objects moving on road networks. To ensure high-quality tracking, the mobile tracker must remain within a distance of ε to the moving object, at all time. To minimize energy consumption, the tracker wants to minimizes walking distance. This type of scenario has been discussed in the context of wireless networks (see [START_REF] Bhuiyan | Local Area Prediction-Based Mobile Target Tracking in Wireless Sensor Networks[END_REF] and [START_REF] Vachhani | Continuous Spatio Temporal Tracking of Mobile Targets[END_REF]). Figure 1 shows an example of an embedding of a planar graph, H, in R 2 , a polygonal curve, T , and a length ε. The dog walks on T from T (0) to T (4) and the person chooses a path in H, from one vertex of H to another vertex. Two points, a and b, are determined on v 4 , v 5 and v 3 , v 4 respectively. A path in H, and a walk on T , that minimize the walking lengths on H and T , are as follows: the dog starts at T (0) and continues on T . The person starts at v 1 and walks on the edges v 1 , v 3 , v 3 , v 4 , and v 4 , v 5 . They move together in a forward direction until the dog reaches the end of the second segment of T and the person reaches the point a on v 4 , v 5 . Then, the dog continues to move forwards until the end of the third segment of T is reached, while the person moves backwards from a to v 4 and then to b. At the final step, they move forwards again together until the dog reaches the end of T and the person reaches v 5 . We show the path in the graph by the sequence of its vertices,

P * = [v 1 , v 3 , v 4 , a, v 4 , b, v 4 , v 5 ].
The structure of this paper is as follows. In Section 2, we discuss preliminaries and define the problem formally. In Section 3, we propose a polynomial time algorithm for the map matching problem introduced. Then, in Section 4, we develop an algorithm with improved time and space complexities. In Section 5, we sketch a solution to a weighted problem variant. Finally, in Section 6, we conclude the paper.

Preliminaries and Definitions

A geometric path in R 2 is a sequence of points in 2D Euclidean space, R 2 . A polygonal curve, or a discrete geometric path, is a geometric path, sampled by a finite sequence of points (called vertices), which are connected by line segments (called edges) in order. Let T : [0, n] → R 2 be a polygonal curve with n segments. A vertex of T is denoted by T (i), i = 0, . . . , n. Let H = V H , E H be a planar graph with a straight-line embedding in R 2 where V H (E H , respectively) is the set of vertices (edges, respectively) of H. In this paper, the geometric embedding of H is crucial and we simply refer to the straight-line embedding of the graph in R 2 as H. A path, P , in H, is a polygonal curve P : [0, 1] → H, such that P ⊂ H and P (0), P (1) ∈ V H .

A parameterization of a polygonal curve, T : [0, n] → R 2 , is a continuous function f : [0, 1] → [0, n], where f (0) = 0 and f (1) = n. Note that a parameterization is corresponding to a walk on T and the interval [0, 1] is representing time during the walk. Walking Length. Let f be a parameterization of a polygonal curve, T . Let D f ⊆ [0, 1] be the closure of the set of times in which f (t) is decreasing (i.e., the movement is backward). The walking length of T is defined by Formula 1, where . is the Euclidean norm and (.) is derivative.

L f (T ) := ||T || + 2 t∈D f ||(T (f (t))) ||dt (1)
Note that if f is monotone (i.e., there is no backward movements on T ), then

L f (T ) = ||T ||.
Problem Definition. Suppose H, T and a length, ε > 0, are given. The objective is to find a path in H and a parameterization of T such that sum of the length of P and the walking length of T is minimized (Formula 2). We consider only paths in the graph and parameterizations of T that guarantee to maintain the leash length at most ε, during the walks.

M ε (H, T ) := inf P ⊂H,f {||P || + L f (T )} (2)
Deformed free-space surface. The free-space diagram is a structure, used to decide whether the Fréchet distance between two polygonal curves is upper bounded by a given ε [START_REF] Alt | Computing the Fréchet distance between two polygonal curves[END_REF]. In [START_REF] Alt | Matching planar maps[END_REF], Alt et al. introduced a 3D structure, called free-space surface, to solve the decision version of their map matching problem.

Here, we use both free-space diagram and free-space surface. However, we modify them slightly, to fit our problem setting.

Let P : [0, 1] → H be a path in H with k + 1 vertices, [p 0 , p 1 , . . . , p k ]. The free-space diagram is the rectangle [0, 1]×[0, 1]
, partitioned into n columns and k rows. It consists of nk parameter cells C x,y , for x = 1, ..., n and y = 1, ..., k. Cell C x,y is the result of the product of two sub-intervals of [0, 1] that are mapped to edge ---------→ T (x -1)T (x) of T and edge ----→ p y-1 p y of P , respectively. We call a point (t

1 , t 2 ) ∈ [0, 1] 2 white if d(T (f (t 2 )), P (t 1 )) ≤ ε,
where d is the Euclidean distance; otherwise, we call it black. It has been shown that the set of all white points inside a cell C x,y is determined by the intersection of an ellipse with C x,y . This set is called the free-space region of that cell. The boundaries of a cell and its corresponding ellipse intersect at most eight times. These intersection points form at most four intervals of white points on the boundary of the cell (i.e., at most one interval per side of the cell). Note that two adjacent cells have the same interval on the shared side between the cells. The union of all cells' freespaces is the free-space (or white-space) of the diagram; it is denoted by W P . The complement of W P is the forbidden-space (or black-space) of the diagram and is denoted by B P . We stretch/compress the columns and rows of the freespace diagram, such that their widths and heights are equal to the lengths of the corresponding segments of T and P , respectively. The resulting diagram is called the deformed free-space diagram and is denoted by F ε (T, P ). In Figure 2, the free-space diagram F ε (T, P * ) is drawn, where

P * = [v 1 , v 3 , v 4 , a, v 4 , b, v 4 , v 5 ]
is a path in H, denoted as a sequence of its vertices.

Note that if the path P contains only a single vertex of H, v i ∈ V H , then F ε (T, P ) is a line segment and its length is equal to the Euclidean length of T . We call this 1D free-space diagram, F i , the deformed free-space line of v i . We denote the left endpoint of F i (i.e., the endpoint corresponding to T (0)) by s i and the right endpoint of F i (i.e., the endpoint corresponding to T (n)) by t i . If P contains only an edge, v i , v j ∈ E H , of H, then F ε (T, P ) has only one row. We call this row the deformed free-space face of v i , v j , and denote it by F j i . Note that F j i and F k j have F j in common. Therefore, gluing F j i and F k j along F j produces a conforming surface. Thus, we can construct the deformed freespace surface as follows. We first lay out the straight-line embedding of H in the xy-plane. For each edge v i , v j ∈ E H , we lay out F j i , orthogonal to the xy-plane, along z axis, such that F i (F j , respectively) is on top of v i (v j , respectively) and s i (s j , respectively) is in the xy-plane. Note that F j i is stretched along z axis from the plane z = 0 to the plane

z = T . Suppose Adj(v j ) is the set of all vertices v k ∈ V H such that v j , v k ∈ E H . We glue F j i to F k j along F j , where a b T P * v 3 v 4 v 4 v 4 v 5
Fig. 2. The free-space diagram Fε(T, P * ) is drawn. WP is the white area and BP is the gray area.

v k ∈ Adj(v j ). Also, we glue

F j i to F i h along F i , where v i ∈ Adj(v h ).
The result is a conforming 3D surface between two planes, z = 0 and z = T , called deformed free-space surface and is denoted by S = H × [0, T ]. Note that s i is on the plane z = 0 and t i is on the plane z = T , i = 1, . . . , |V H |. The union of the white-space (black-space, respectively) of all faces of S is called the whitesurface (black-surface, respectively) and is denoted by W (B, respectively). For the given planar graph H, the polygonal curve T , and the length ε in Figure 1, the corresponding deformed free-space surface is shown in Figure 3, from two points of views. Since the white-space of each cell of any F j i is convex, for simplicity, we just draw the white-space intervals on the boundary of the cells. In this figure, the red dashed polygonal curve is a path on the white-surface W, from s 1 to t 5 , that realizes

P * = [v 1 , v 3 , v 4 , a, v 4 , b, v 4 , v 5 ]
, in H, and a parameterization of T , that is an optimal solution to our problem setting. It intersects the following free-space faces sequentially:

F 3 1 , F 4 3 , F 5 4 , F 4 3 , F 5 4 .

Algorithm

In this section, we first transform the map matching problem to a shortest path problem on a weighted graph, G = V, E ; this yields a polynomial time algorithm. Before discussing the construction of G, we introduce a set of Steiner points on the boundary of the cells of S. Steiner Points. We position Steiner points so as to create intervals on the boundary of the cells of S. There are two types of intervals, Type 1 and Type 2. We classify the Steiner points based on the type of the intervals that they belong to. We denote the set of Type 1 (Type 2, respectively) Steiner points by S 1 (S 2 , respectively). The red dashed polygonal curve is a path on the white-surface that realizes an optimal solution to our problem setting.

Type 1. We say an interval is Type 1, if it lies completely in a plane, z = c, parallel to the xy-plane, where c is a constant. Each deformed free-space face, F j i , may have n + 1 Type 1 intervals, FI j i ( ), = 0, . . . , n, shared between its cells (where n is the number of edges in T ). For each interval FI j i ( ), we project the endpoints of FI j i ( ) orthogonally to all FI j i (k), k = . If the line segment from an endpoint of FI j i ( ) to its projection on FI j i (k) lies in the free-space of F j i and the projection point is not identical with an endpoint of FI j i (k), then we take the projection point as a Type 1 Steiner point (see Figure 4). The set of all Steiner points, obtained by the projections on F j i , for all v i , v j ∈ E H , is denoted by S 1 .

Type 2. We say an interval is Type 2 if it lies completely on a deformed freespace line. As we mentioned in Section 2, a plane z = c corresponds to a point on the given trajectory T . The intersection of z = c and S is an instance of H, denoted by H c . Note that some part (possibly empty) of H c is in W. Let z = h j be the corresponding plane of T (j), a vertex of T . The part of the deformed free-space surface, S, between the two parallel planes, z = h j-1 and z = h j , j = 1, . . . , n, corresponds to edge ---------→ T (j -1)T (j) of T . We denote this part of S by T j j-1 = H × [h j-1 h j ]. In T j j-1 , j = 1, . . . , n, there is at most one Type 2 interval per vertex v i ∈ V H . We denote these Type 2 intervals by T I j j-1 (i), i = 1, . . . , |V H |. Suppose z = c is the plane that is passing through an endpoint, p, of T I j j-1 (i). Let the intersection of z = c with T I j j-1 (k), k = i, be q k . Note that both p and q k are on the graph H c . Then, if q k is not an endpoint of T I j j-1 (k) and there is a path, from p to q k , in H c , that is in W, then q k is a Type 2 Steiner Point. The set of all Type 2 Steiner points is denoted by S 2 . An example is given in Figure 5. Suppose it is T j j-1 , for j = 1. In this example, there are four yellow intervals, T I j j-1 (1), T I j j-1 (3), T I j j-1 (5), and T I j j-1 [START_REF] Alt | Computing the Fréchet distance between two polygonal curves[END_REF]. The black points show the interval endpoints and red points show the Type 2 Steiner points. For simplicity, only two, out of eight planes, are drawn. The plane z 3 (z 6 , respectively) is passing through an endpoint of T I j j-1 (3) (T I j j-1 (6), respectively). The intersections of z 3 with T I j j-1 (5) and T I j j-1 (6) are Steiner points. However, the intersection of z 6 with T I j j-1 (3) is not a Steiner point. Constructing Graph. Now, we explain the construction of G = V, E . Recall that the white-surface (the white-space of S) is denoted by W. The vertices of W are the end points of the intervals on the boundary of the cells in S (at most 4 intervals may exist per cell). We denote the set of vertices of W by V W . The set of vertices, V , of G, is V = V W ∪ S 1 ∪ S 2 . Note that V contains all s i and t i if they are in W. Every two vertices, v 1 , v 2 ∈ V , that are on the boundary of a cell, are linked by two directed edges in E, from v 1 to v 2 , v 1 , v 2 , and vice versa, v 2 , v 1 . The weight of an edge e = v 1 , v 2 ∈ E, is its length in the L 1 metric, |e| 1 . Obtain an Optimal Solution. In order to have an optimal solution, at least one s i and one t j , i, j = 1, . . . , |V H |, must be in W. The main steps of the algorithm are as follows:

v i , v j
-Find all the vertices in V H that are in ε distance of T (0) (T (n), respectively), v i1 , . . . , v i k 1 (v j1 , . . . , v j k 2 , respectively). -Add an extra node, s , to G, and add k 1 extra directed edges, s , s i1 , . . . , s , s i k 1 , to E. The weight of these k 1 edges are set to zero. Analogously, add another extra node, t , to G, and add k 2 extra directed edges, t j1 , t , . . . , t j k 2 , t , to E. The weight of these k 2 edges are also set to zero. -Find a shortest path, from s to t , in G. Note that if there is no path from s to t in G, then there is no solution for the given leash length. -Remove s and t from the head and tail of the shortest path. The remaining path is from one s i to one t j . It gives an optimal solution to our problem setting.

Note that, a vertex of G (except s and t ) is also represented by a point in W. Therefore, the geometric embedding of a path, from one s i to one t j , in G, is constructed by connecting the consecutive vertices of the path in W by line segments.

Observation 1 Let Π be a path in the white-space, W, of a deformed free-space surface, S, from one s i to one t j . Π realizes a path, P : [0, 1] → H, in H, and a parameterization, f : [0, 1] → [0, n], of T , that maintain the leash length at most ε, for all t ∈ [0, 1].

Constructing a path in H. We can construct a path P in H, from the given path Π in W, as follows. As we mentioned earlier in this section, we have two types of intervals on the boundary of the cells, Type 1 and Type 2. A Type 1 interval lies completely on a plane, z = c, parallel to the xy-plane. A Type 2 interval lies completely on a deformed free-space line, F i . The path Π intersects a sequence of intervals (of both types). The path P in H is constructed by processing the intervals in this sequence. For each interval in this sequence, if it is Type 2 interval, on F i , then we append v i to the tail of P . If it is Type 1, then the intersection point, q, of Π and that interval, is appended to the tail of P , as a vertex of P . Note that q may not be a vertex of H. However, it is a point on an edge of H. At the end, we connect the consecutive vertices in P by straight line segments. Correctness. To establish the correctness, we use norms in two spaces: [START_REF] Zheng | Trajectory Data Mining: An Overview[END_REF] the Euclidean space of the embedding of the input graph and the polygonal curve, called the input space, (2) the deformed free-space surface, called the configuration space. In the input space, we denote the Euclidean length of a polygonal curve T by T . We also defined walking length of T , L f (T ), based on a parameterization f . Note that if f is a monotone parameterization, then L f (T ) = T . In configuration space, a path from an s i to a t j in W, is also denoted by its vertices, Π :

s i = p 1 , p 2 , . . . , p k = t j . The length, |.| 1 , of each segment of Π is calculated by the L 1 metric.
The length of a path, |Π| 1 , is the sum of the length of its segments. Lemma 1 is at the heart of the correctness proof. This section is concluded by a corollary to Lemma 1 and Observation 1, that is, in order to find a solution for our problem setting, it suffices to find a shortest path from s to t in G.

Lemma 1 For any path Π : s i = p 1 , p 2 , . . . , p k1 = t j in W, there is a path Π :

s i = p 1 , p 2 , . . . , p k2 = t j in W such that Π ⊂ G and |Π | 1 ≤ |Π| 1 .
Proof. The path Π intersects a sequence, SF , of deformed free-space faces. Every two consecutive faces in SF share a deformed free-space line. Therefore, we can unfold the free-space faces in the sequence, along the shared free-space lines. The result is a 2D free-space diagram, denoted by F ε (SF ). W.l.o.g., we can assume that F ε (SF ) is axis aligned in R 2 . The path Π is also unfolded into a 2D path in the white-space, W SF , of F ε (SF ). Note that unfolding does not change the length of a path. As an example, in Figure 6a, the result of unfolding the faces that are intersected by the red dashed polygonal curve in Figure 3, is shown.

Let Π opt : s i = q 1 , q 2 , . . . , q k3 = t j be a L 1 shortest path, from s i to t j , in W SF . Then, |Π opt | 1 ≤ |Π| 1 . To prove the lemma, it suffices to show that there is a path

Π ⊂ G = V, E , from s i to t j , in W SF , such that |Π | 1 = |Π opt | 1 .
We know that the vertices of Π opt are endpoints of some intervals on the boundary of the cells of F ε (SF ) [START_REF] Ghosh | An output-sensitive algorithm for computing visibility graphs[END_REF] (the well known rubber band property of shortest paths). Therefore, the vertices of Π opt are in V (i.e., the set of vertices of G). Thus, it is sufficient to show that for each edge, ---→ q i q i+1 , of Π opt , there is a path, π qiqi+1 , from

a b T s 1 t 5 a b q i q i+1 v 1 , v 3 v 3 , v 4 v 4 , v 5
q i to q i+1 , in G, that | ---→ q i q i+1 | 1 = |π qiqi+1 | 1 .
Two cases arise depending on whether ---→ q i q i+1 lies completely within a row (or a column) of ε (SF ), or not. We need a definition before discussing these cases. We assume that F ε (SF ) is an axis-aligned rectangle in a 2D Cartesian coordinate system, where the x-axis corresponds to T (Figure 6a). We say a path Π ∈ W SF is x-monotone (y-monotone, respectively), if any vertical (horizontal, respectively) line intersects it at most ones. Π is said to be xy-monotone, if it is both x-and y-monotone.

Case 1. In this case, ---→ q i q i+1 lies completely within a column (or a row) of F ε (SF ). Here, we discuss the case when it lies within a column (see Figure 6b); the arguments are analogous for case of a row. W.l.o.g we assume that ---→ q i q i+1 is xy-increasing. The other cases are symmetric. Edge ---→ q i q i+1 intersects a sequence of horizontal intervals, I z , within a column. They are sorted based on their y coordinates. We construct π qiqi+1 sequentially and always denote the last vertex appended to π qiqi+1 by π last . Initially, π qiqi+1 contains only q i and π last = q i . The sequence of intervals are processed sequentially. Suppose we processed interval I z and now we want to process I z+1 . We project orthogonally from π last to I z+1 . If the projection point exists (i.e., the perpendicular line from π last to I z+1 intersects I z+1 ), then append the projection point on I z+1 to π qiqi+1 and update π last . Otherwise, the closest endpoint of I z+1 to π last is appended to π qiqi+1 and we update π last . When all intervals, I z , have been processed, q i+1 is appended to π qiqi+1 .

Since the sorted list of intervals withing a column are traversed by π qiqi+1 sequentially, the path π qiqi+1 is y-monotone. Also, by construction, each vertex of π qiqi+1 either has the same x as its preceding vertex in π qiqi+1 (i.e., it is the result of the orthogonal projection) or its x is greater than its preceding vertex's x (since the orthogonal projection does not exist and ---→ q i q i+1 is xy-increasing inside the white-space). Therefore, the path π qiqi+1 is x-monotone. Thus, the path π qiqi+1 is xy-monotone. We know that the L 1 length of two xy-monotone paths that have the same starting and ending points, are equal. Therefore, | ---→

q i q i+1 | 1 = |π qiqi+1 | 1 .
Now, we prove that π qiqi+1 ⊂ G = V, E . It suffices to show that each vertex of π qiqi+1 is in V and between every two consecutive vertices of π qiqi+1 there is an edge in E. Each vertex of π qiqi+1 is either the result of the orthogonal projection or an endpoint of an interval. Therefore, each vertex is either a Steiner point or a vertex of the white-surface. In both cases, the vertex is in V . In addition, between every two consecutive vertices of π qiqi+1 there is an edge in E because every two consecutive vertices of π qiqi+1 lie on the boundary of a cell and, by the construction of G, all members of V that lie on the boundary of a cell are linked by edges in E.

Case 2. In [START_REF] Gheibi | Weighted Minimum Backward Fréchet Distance[END_REF], Section 4, Lemma 4, it is proved that if ---→ q i q i+1 does not lie completely within a row and within a column of F ε (SF ), then there is a xymonotone path π qiqi+1 , from q i to q i+1 , such that its edges lie completely within a row and within a column of F ε (SF ). For each edge of π qiqi+1 , we apply case 1. Then, we concatenate the resulting xy-monotone paths for edges of π qiqi+1 , to obtain π qiqi+1 . Since, xy-monotone paths for edges of π qiqi+1 are in G (as we proved in Case 1), the resulting path, π qiqi+1 , is a xy-monotone path in G.

Therefore, | ---→ q i q i+1 | 1 = |π qiqi+1 | 1 .
Corollary 1 For any pair of s i and t j , if t j is reachable from s i by a path in W, then there is a path from s i to t j , in G, that is a L 1 shortest path in W.

Corollary 2 A shortest path in G, from s to t , yields an optimal solution for our problem setting.

Proof. Let Π opt be a shortest path in G, from s to t . We remove s and t from the head and tail of Π opt . The result, Π ij , is a shortest path from s i to t j . Therefore, among all possible shortest paths Π k , for s k and t , k, = 1, . . . , |V H |, the pair (s i , t j ) has a shortest L 1 shortest path, Π ij . By Corollary 1, Π ij ⊂ G is a L 1 shortest path in W. Each point on Π ij is corresponding to a point, p, on H and a point, q, on T , such that the Euclidean distance of p and q is less than ε. By Observation 1, Π ij , is corresponding to a path, P , in H, from v i ∈ V H to v j ∈ V H , and a parameterization, f , of T . The summation of the Euclidean length of P , P , and the walking length of T , L f (T ), is equal to the L 1 length of Π ij . Since Π ij is a shortest L 1 shortest path, P and f minimize the matching cost, M ε (H, T ) (Equation 2).

Theorem 1 Let H be a planar graph with a straight-line embedding in a plane, T be a directed polygonal curve, and ε > 0 be a distance. A path, P : [0, 1] → H, between two vertices of H, and a parameterization, f , of T , that minimize the sum of the walking length of T and P , can be found in polynomial time and space. It is guaranteed that at any time t ∈ [0, 1], the Euclidean distance between P (t) and T (f (t)) is at most ε.

Proof. The correctness follows directly from Corollary 2. The deformed freespace surface, S, has O(mn) cells, where m (n, respectively) is the number of edges of H (T , respectively). Each cell of S has at most four intervals and at most O(m + n) Steiner points on its intervals. Therefore, the graph G has O (mn(m + n)) vertices and O mn(m + n) 2 edges (including the extra edges that connect s and t to the graph). In addition, it takes O n 2 time to compute all Type 1 Steiner points for each free-space face. Therefore, computing S 1 takes O mn 2 time. In order to compute Type 2 Steiner points, we use breadth first search for each interval endpoint to propagate the projection on the instance of the graph, H c , in the plane z = c. Therefore, computing S 2 takes O nm 2 time. At the end, it is possible to find a shortest path in G, from s to t , in O mn(m + n) 2 time, by using Dijkstra's algorithm. Therefore, both the total time and space complexities are O mn(m + n) 2 .

Improvement

In Section 3, we showed that the graph G = V, E contains a path that yields an optimal solution for our problem setting. The bottleneck in the time complexity of the algorithm in Section 3 is due to the number of edges of G. In this section, we construct a new graph G = V, E , such that |E | < |E| and it preserves the connectivity information of G. More precisely, if there is a path, from v i ∈ V to v j ∈ V , in G, then there is a path, from v i to v j , in G , with the same L 1 length.

Based on the construction of G, there are at most O(m + n) vertices in V (including the interval endpoints and Steiner points) on the boundary of each cell, C, of S. We connect these O(m + n) vertices by a linear number of edges, in E , as follows. The weight of each edge in E is equal to its L 1 -length. Let T , B, L, and R be the intervals on the top, bottom, left and right side of C, respectively. Suppose cell C is in a 2D Cartesian coordinate system and the vertices on each interval I ∈ {T, B, L, R} are sorted by x and y. Every two adjacent vertices, v i and v i+1 , on I, are linked by two directed edges, v i , v i+1 and v i+1 , v i (Figure 7). Every two of the eight interval endpoints are linked by two directed edges assuming they are not identical. A vertex v i on interval L (T , respectively), is linked by two directed edges to another vertex v i on R (B, respectively) if v i has the same y (x, respectively) coordinate as v i ; the two edges are denoted by v i , v i and v i , v i , respectively. By this approach, each vertex of G on the boundary of C is connected to a constant number of vertices of G on the boundary of C. It is now straightforward to prove the following lemma.

Lemma 2 Let v i and v j be two vertices, in V , on the boundary of a cell, C, of S. There is a path, from v i to v j , in G , that has the same L 1 length as the direct line segment between them.

Corollary 3 There is a path in G that realizes an optimal solution for our problem setting. T B L R Fig. 7. A cell of the free-space surface is drawn. The red solid line segments show the four intervals on the boundary of the cell. The arcs show the edges in E that connect every two adjacent vertices of G , on each interval. The dashed black line segments show the edges in E that connect a vertex with its orthogonal projection on the opposite side of the cell. The dash dotted blue line segments show some of the edges that connect endpoints of the intervals. For simplicity, we did not draw all of them.

Theorem 2 Let H be a planar graph with a straight-line embedding in a plane, T be a directed polygonal curve, and ε > 0 be a distance. A path, P : [0, 1] → H, between two vertices of H, and a parameterization, f , of T , that minimize the sum of the walking length of T and P , can be found in O(nm(n + m) log(nm)) time and O(nm(n + m)) space, where n (m, respectively) is the number of edges of T (H, respectively). It is guaranteed that at any time t ∈ [0, 1], the Euclidean distance between P (t) and T (f (t)) is at most ε.

Proof. The correctness follows directly from Corollary 3. The number of vertices and edges of G (and the total space complexity) is upper-bounded by O(nm(n+ m)). Using Dijkstra's algorithm, we find a shortest path in G , from s to t . Therefore, the time complexity of our algorithm is O(nm(n + m) log(nm)). Note that if there is no pair of (s k , t ), k, = 1, . . . , |V H |, in a connected component of G , then there is no feasible solution.

Weighted non-planar graphs

We assumed that the input graph H is planar. That makes the illustration of the algorithm easier since the faces of the free-space surface do not intersect except at the boundary of the faces. However, all lemmas and theorems, derived in Section 3 and 4, are proved without making an assumption regarding the planarity of H. Therefore, the algorithm proposed in this paper remains correct for any graph for which a straight-line embedding in a plane is provided (see [START_REF] Alt | Matching planar maps[END_REF], Section 2.7). In the embedding, the edges of the graph may intersect. Transition from one edge to another is allowed only at a vertex.

We also assumed that H is unweighted. Here, we sketch how the proposed algorithm can be generalized to also handle the problem instance, when H is weighted. Suppose that each edge of H has a non-negative, real weight. A weight could represent the cost of moving on the edge of the graph. The edges of the input polygonal curve T could also have weights capturing the costs of moving forwards and backwards. The objective is to find a path in H whose weighted walking length is minimized. In the weighted problem setting, inside each cell of the free-space surface, there are two weights, one corresponding to an edge of T and one corresponding to an edge of H. These weights are fixed inside the cell and do not change. Therefore, in the construction of G or G , instead of computing the L 1 length for each edge, e, we compute the orthogonal projections of e onto H and T . Then, we multiply the projection lengths with the corresponding weights, and the sum of these multiplications is the weight that we assign to e. The remaining parts of the algorithm remains the same and the time and space complexities do not change.

Conclusion

In this paper, we discussed a geometric algorithm for the map matching problem that minimizes the walking length. We established that this problem setting is dual to a weighted shortest path problem. Then, we proposed an algorithm with O (mn (m + n) log(mn)) time and O (mn (m + n)) space complexities, where m (n, respectively) is the number of edges of H (T , respectively). At the end, we discussed that the proposed algorithm is easily adaptable to handle weighted non-planar graphs. It is still open if we can improve the proposed algorithm further, for planar graphs. The main challenge here is the existence of cycles in the input graph and propagation through the cycles.

Fig. 1 .

 1 Fig. 1. An embedding of a planar graph, H, a polygonal curve, T , and a length ε are given. The path P * = [v1, v3, v4, a, v4, b, v4, v5], in H, is a part of a solution to the map matching problem instance. The edges of H that P * lies on, are illustrated in bold.

Fig. 3 .

 3 Fig.3. The free-space surface for the example of Figure1is drawn from two different viewpoints in 3D. The yellow line segments show the intervals on the cell boundaries. The red dashed polygonal curve is a path on the white-surface that realizes an optimal solution to our problem setting.

Fig. 4 .

 4 Fig. 4. The free-space face F j i is drawn. The endpoints of the intervals, FI j i ( ), are shown by points and the Type 1 Steiner points are shown by squares.

Fig. 5 .

 5 Fig. 5. An example of T j j-1 , for j = 1, is drawn. In this example, there are four intervals that are shown by yellow color. The black balls show the interval endpoints and red balls show the Type 2 Steiner points.

Fig. 6 .

 6 Fig. 6. a) The result of unfolding the sequence of deformed free-space faces that are intersected by the red dashed polygonal curve in Figure 3. It is a 2D free-space diagram, Fε(SF ). The red dashed polygonal curve is shown after unfolding. b) Illustration of case 1 in the proof of Lemma 1.

Acknowledgment

The authors would like to thank Carola Wenk for suggesting this topic and constructive comments, and Omid Gheibi for valuable discussions.


Thibault Gateau

Caroline P Carvalho

Mai-Huy Le

Frédéric Dehais

Considering Human's Non-Deterministic Behavior and his Availability State When Designing a Collaborative Human-Robots System

The objective of this study is to design a humanrobots system that takes into account the non-deterministic nature of the human operator's behavior. Such a system is implemented in a proof of concept scenario relying on a (MO)MDP decision framework that takes advantage of an eye-tracker device to estimate the cognitive availability of the human operator, and, some human operator's inputs to deduce where he is focusing his attention. An experiment was conducted with ten participants interacting with a team of autonomous vehicles in a Search & Rescue scenario. Our results demonstrate the advantages of considering the cognitive availability of a human operator in such a complex context and also the interest of using such a decisional framework that can formally integrate the non-deterministic outcomes which model the human behavior.

I. INTRODUCTION

Autonomous Air Vehicles (AAVs) systems are increasingly becoming more "autonomous" in achieving tasks in a variety of operational contexts (e.g. border security, inspection of hazardous areas, military operations). Despite these technical advances, the human is "still vital" [START_REF] Parasuraman | Humans: Still vital after all these years of automation[END_REF] to operate these automated systems, as stated by regulation authorities. Generally, the role of the human operator is to supervise the achievement of the mission, to take over during failure, and to make appropriate tactical and ethical decisions [START_REF] Galliott | Military Robots: Mapping the Moral Landscape[END_REF].

However, AAVs safety reports state that the human factors issues were involved in 80% of mishaps [START_REF] Williams | A summary of unmanned aircraft accident/incident data: Human factors implications[END_REF]. A deep analysis of these accidents reveals that they usually involve humanautomated systems interactions deficiencies rather than human errors per se. In fact, the developments of AAVs were generally achieved without taking into account the integration of the human operators [START_REF] Schurr | Impact of human advice on agent teams: A preliminary report[END_REF]: the human operator is commonly considered as a providential agent to take over when the artificial systems fail to adapt [START_REF] Schurr | Improving adjustable autonomy strategies for time-critical domains[END_REF]. The complexity of theses systems reduces the human operator's abilities to take immediate and appropriate actions [START_REF] Woods | Learning from automation surprises and" going sour" accidents[END_REF], especially after long periods of inactivity while monitoring AAVs by a user interface [START_REF] Boyer | Investigating mental workload changes in a long duration supervisory control task[END_REF], [START_REF] Durantin | Characterization of mind wandering using fNIRS[END_REF]. In addition, poor user interface design and careless implementation of authority sharing confuse the human operator when mental workload exceeds human capacity [START_REF] Durantin | Using near infrared spectroscopy and heart rate variability to detect mental overload[END_REF], resulting in poor human-system performance [START_REF] Tessier | Authority management and conflict solving in human-machine systems[END_REF]. One has to consider that this topic is getting critical as research efforts are currently oriented on the control of multiple AAVs by one (or few) human operator(s) [START_REF] Cummings | Automation architecture for single operator, multiple uav command and control[END_REF], [?]. 1 Université Fédérale de Toulouse Midi-Pyrennées, ISAE-SUPAERO, Institut Supérieur de l'Aéronautique et de l'Espace, 10, Av Edourd Belin -BP 54032 -31055 Toulouse Cedex 4 -FRANCE name.surname@isae-supaero.fr A promising avenue to deal with these issues is to consider that robot and human abilities are complementary and are likely to provide better performance when joined efficiently than when used separately. This approach, known as mixedinitiative [START_REF] Allen | Mixed-initiative interaction[END_REF] defines the role of the human and of the artificial agents according to their recognized skills. A key issue in designing such mixed-initiative interactions is to implement a decision system. This latter defines the role and the authority of human and artificial agents, while estimating the capabilities of the human, the robotic agent's status, and the achievement of the mission. Therefore, a first step is to implement methods and tools that can assess the "cognitive state" of the human operator. Many researches have been achieved in order to implement physiological or neurophysiological based inference systems to deduce cognitive limitations (eg. [START_REF] Gateau | Realtime state estimation in a flight simulator using fNIRS[END_REF], [START_REF] Régis | Formal detection of attentional tunneling in human operator-automation interactions[END_REF]) in order to dynamically adapt human-robot interactions [START_REF] Tessier | Authority management and conflict solving in human-machine systems[END_REF]. Despite these previous works, more research is still needed to improve their reliability due to intra and inter-individual differences. For instance, eyetracking techniques appear to be suitable in assessing the human operator's cognitive state in the context of a human-AAVs interaction. Desktop eye-trackers can be easily set up below the user's interface and allow the measurement of the operator's eye movement. Many metrics can be collected in real time (for a review see [START_REF] Holmqvist | Eye tracking: A comprehensive guide to methods and measures[END_REF]) to estimate cognitive impairment such as attentional tunneling [START_REF] Regis | Ocular metrics for detecting attentional tunnelling[END_REF], attentional confusion [START_REF] Dehais | Automation Surprise in Aviation: Real-Time Solutions[END_REF] or drowsiness [START_REF] Ji | Real-time eye, gaze, and face pose tracking for monitoring driver vigilance[END_REF].

The main idea of considering human operator's cognitive state is to merge such measurements with the states of the robots as well as the level of achievement of the mission. A decisional policy needs to be set up to manage the dynamics of the human-AAVs interactions, and possibly the AAVs' behavior. On one hand, decision-making systems can be governed by a plan [START_REF] Talamadupula | Coordination in human-robot teams using mental modeling and plan recognition[END_REF] to handle human-robot teaming scenarios using classical planning approaches. The drawback of such solutions is the impossibility to consider uncertainties on actions, partial observable states or potentially nondeterministic behavior of the human operator. On the other hand, the decision-making system can be driven by a policy resulting from the resolution of a (Partially Observable) Markov Decision Process ((PO)MDP) [START_REF] Smallwood | The optimal control of partially observable Markov processes over a finite horizon[END_REF], where the two types of agents, human or artificial, may face unexpected situations during the mission or having a partial observation of the state. These outcomes can be modeled as probabilistic effects of actions and as a probabilistic observation function respectively. In this sense, a framework has been proposed [START_REF] Taha | A pomdp framework for modelling human interaction with assistive robots[END_REF], which is able to adapt itself to the user's intention, getting feedback from the user in terms of satisfaction. A different way to manage interactions using POMDPs has been studied for assisting persons with dementia during hand-washing [START_REF] Hoey | Assisting persons with dementia during handwashing using a partially observable Markov decision process[END_REF]. [START_REF] Nikolaidis | Efficient model learning from joint-action demonstrations for human-robot collaborative tasks[END_REF] [START_REF] Nikolaidis | Efficient model learning from joint-action demonstrations for human-robot collaborative tasks[END_REF] proposed a MOMDP framework where the human is seen as a partially observable state variable in a joint action context involving a human and a robot that automatically learns the human's specific behavior. In another example [START_REF] Souza | MOMDPbased target search mission taking into account the human operator's cognitive state[END_REF], a decomposition -between the state vector of the robot considered as fully observable while the operator's cognitive state is considered as partially observable -is explored to model the mixed-initiative framework as a Mixed Observability Markov Decision Process (MOMDP) [START_REF] Ong | Planning under uncertainty for robotic tasks with mixed observability[END_REF]. However, this relevant work does not evaluate experimentally the human operator's behavior once integrated in such a system and the study only provides simulations results.

In this paper, we propose to model and to experimentally evaluate a scenario involving a team composed of AAVs and a human operator. The decisional policy is driven by a MOMDP model that aims to take into account the human's inputs and measures about his cognitive state in order to adapt the AAV's decision autonomy and so optimize the performance of the team. To the best of the authors' knowledge, such experimental evaluation has never been addressed yet. The main contribution of this work consists in:

(1) Including the human's cognitive state in the mission planning model, i.e considering the human's cognitive availability in the model to obtain a policy which drives the mixed-initiative process. The human cognitive's availability is derived from basic eye-tracking metrics and is also used in the on-line state estimation process, which provides appropriate actions for the whole system.

(2) Handling imperfections of the state transition function in the planning model, using an adapted decision model framework. The human operator's behavior is not deterministic. For instance, the reaction time can depend on the person solicited and on the context in which the person evolved. We propose therefore to take advantage of a decision framework that handles such uncertainties as probabilistic actions effects and as probabilistic observation function, and to apply such a decision framework to the system (AAVs and human operator). To this end, the (MO)MDP represents an elegant framework to manage such a decision problem.

(3) Building a proof of concept scenario tested in a well controlled experiment based on a collaborative mission including one human operator and AAVs.

The paper is organized as follows: first we present the materials and methods used in our study case. Next, the collaborative human-AAVs mission problem is formalized based on a MOMDP framework. After, the experimental protocol, metrics and the results obtained are presented. Finally we conclude this paper by discussing future works.

II. EXPERIMENTAL SCENARIO

Humans are generally confronted to concurrent tasks while performing system supervision. In this sense, the human operator was confronted to a dual-task paradigm in the scenario proposed. The first task, defined as Search & Rescue task, consisted of a classical search and rescue mission [START_REF] Parker | Distributed intelligence: Overview of the field and its application in multi-robot systems[END_REF]. The human operator had to collaborate with autonomous artificial agents to perform target identifications by means of an user interface as shown on the left side of Figure 1. The second task, defined as Short-Term Item Memorization task, consisted of a series of digits that the human operator had to memorize and report via the user interface as shown on the right side of Figure 1.

It was expected that the volunteers to have a lower performance when achieving both tasks than when achieving each of them separately.

A. Search & Rescue Task

Targets appeared on the simulation screen (the upper part of the AAVs monitoring interface on Fig. 1). A starting point base was defined for the team of AAVs that must: (1) visit target location; (2) take pictures of the target; (3) identify the target; (4) bring the target back to the base if it was a human or a material; (5) avoid bringing back targets that are neither human nor material.

Drones could request help from the human operator for target identification, according to the decisional policy computed (Cf. Section III-B.8). A request to the operator was materialized using a yellow bounding box visual effect on the target, in the simulation screen. For the human operator, the identification task was emulated by a dot counting, depending on the target's real nature (the lower part of Search & Rescue task interface -Fig. 1). The identification picture was available for a limited time on the screen, depending on the decisional policy, adding some time-pressure to the task [START_REF] Cummings | Automation architecture for single operator, multiple uav command and control[END_REF], and also to avoid instant target identification, subitizing behaviors were refrained by using more than six dots [START_REF] Mandler | Subitizing: an analysis of its component processes[END_REF].

B. Short-term Item Memorization Task

A key executive function when operating complex systems is the working memory [START_REF] Gateau | Realtime state estimation in a flight simulator using fNIRS[END_REF]. In this sense, we chose to design a well know task, the digit span task, to measure Working Memory (WM) storage and update capacity [START_REF] Miller | The magical number seven, plus or minus two: some limits on our capacity for processing information[END_REF]. The volunteers were instructed to retain digits displayed on half of a 21 inches gray screen (on the right part, resolution 1680x1050), via a script written in Python 2.7.3. Volunteers were asked to use a visual keypad to enter the sequence of seven digits using a mouse cursor, 500 ms after the last digit was presented. A visual keypad was preferred to refrain participant to look outside the screen while entering his response.

A series of digits (i.e. a trial) started with a 13 secondcross fixation. A variable stimulus onset asynchrony (SOA) of 1, 2, or 3 seconds extended the presentation of the fixation cross. Each digit was presented during 500 ms at a rate of 500 ms (i.e. an inter-stimulus interval of 500 ms). A trial ended after a 8 seconds response window. Therefore, a trial lasted 30 seconds in average. The complete task was composed of 20 trials and was designed to last in average 10 minutes.

The difficulty of such a task depends on the number of digits. This way, the difficulty was set to seven digits. In fact, the human WM is fundamentally limited to store seven (more or less two) items [START_REF] Miller | The magical number seven, plus or minus two: some limits on our capacity for processing information[END_REF]. Hence, volunteers were expected to get sufficiently involved into this WM task, without actually reaching mental overload when achieving this task alone [START_REF] Durantin | Using near infrared spectroscopy and heart rate variability to detect mental overload[END_REF].

III. COLLABORATIVE PROBLEM FORMALIZATION

A. MOMDP model formalism

The MOMDP model [START_REF] Ong | Planning under uncertainty for robotic tasks with mixed observability[END_REF] is an extension of the POMDP model [START_REF] Smallwood | The optimal control of partially observable Markov processes over a finite horizon[END_REF], which explores the fact that certain state variables are fully observable. Formally, a MOMDP is a tuple (X , Y, A, Ω, T X , T Y , Ω, R, γ, b 0 ), where: X (resp. Y) is the set of fully (resp. partially) observable state variables; A is the set of actions; Ω is the set of observations; T X (resp. T Y ) is the conditional probability transition function such that, T X (x, a, x , y) = p(x |x, y, a), ∀(x, a, x , y) ∈ X × A × X × Y (resp. T Y (x, y, a, x , y ) = p(y |x, y, a, x ), ∀(x, y, a, x , y

) ∈ X × Y × A × X × Y); O is the conditional probability observation function such that, O(o, a, x , y ) = p(o|x , y , a),∀(o, a, x , y ) ∈ Ω × A × X × Y; R : X × Y × A → R
is the reward function defined for a state-action pair; γ ∈ [0, 1[ is the discount factor; and, b 0 = (x 0 , b Y0 ) ∈ B is the initial probability distribution over states, where B is complete belief state space.

As the artificial decider agent does not have access to the true state of the world, it maintains a probability distribution over the partial observable state variables (initialized by b 0 = (x 0 , b Y0 )). And, after each action a executed, each next fully observable state x reached and each observation o perceived, the probability distribution must be updated using the Bayes' rule (η being a normalization constant): Note that the belief state b is a tuple formed by the couple

(x, b Y ), with b Y ∈ B Y (x). B Y (x) is a belief state space over Y conditioned by x: B Y (x) = {(x, b Y ), b Y ∈ B Y }. Notice that B Y (x) is a sub-space of the complete belief state space B, such that B = x∈X B Y (x).
The MOMDP solving can be seen as an optimization problem for finding a set of policies π x : B Y → A, which maximizes the criterion:

π * x ← arg max πx∈Π Eπ x ∞ t=0 γ t r((xt, b t Y ), π((xt, b t Y ))) b0) .
Such a criterion allows the definition of a value function and the use of dynamic programming to solve the problem. Then, the resulting set of policies π x : B Y → A, that can be extracted from the optimal value function, defines an optimal action which maximizes the expected gain for each reachable belief states (x, b Y ), considering that optimal actions would be performed for all possible successive belief states. For more details about MOMDP resolution, please refer to [START_REF] Ong | Planning under uncertainty for robotic tasks with mixed observability[END_REF].

B. MOMDP model for the collaborative mission

The collaborative mission consists in identifying targets, that can be either human, material, or none, and to bring to the base only relevant targets (Cf. Section II-A). The drones are able to sense the type of targets, but they have a less reliable probability of good classification than the one of the human operator. The human operator plays cooperatively with the drones in order to help them with a sensing task.

1) States, actions and observations: The MOMDP state space is composed by the cross-product between a set of state variables: the mission state of each drone i (drone i ), the input of the operator (operator-input), the availability of the human operator and the type of each target j (target j -type). For the following, let n d ∈ N * be the number of drones,

n t ∈ N * be the number of targets, i ∈ [[1; n d ]] and j ∈ [[1; n t ]].
Possible values for each discrete state variable are:

• dronei : { going-to-target1. We consider that the drones' position, the operator's input and the availability are fully observable state variables; and, the target's type is a partially observable state variable.

Therefore we have some discrete observation variables. The drones and the operator can observe the target's type. So, we have:

• o-dronei : {ohuman, omaterial, onone}, ∀i ∈ [[1; nd]]; • o-operator : {ohuman, omaterial, onone};
As actions are defined in a centralized way, the joint action set is composed by the cross product between drones' individual actions and user's interface possible actions. The user interface (i.e. ground station) action is used to relay a drone's requests to the human operator to look at a target (e.g. launching a yellow bounding box) in order to bring more reliable measures about the target type. So, we have: The transition function defines the dynamic of the model for each agent. We choose to define a clock of 1 second to trigger the model (ask for a new action). Therefore, all transition probabilities related to actions' effects are defined hereafter based on this clock.

2) Effects of drones' actions on drones' states: the action going-to-target j moves the drone i to the going-to-target j state with a probability of 1.0, except if the drone i is on the going-to-target j state. In this case, it moves the drone i to the at-target j state with a probability of 0.2, and remains on the going-to-target j state with 0.8. This helps us to indirectly model the drones' travel time for a movement that lasts approximately 5 seconds in our application. The action check does not change the state of drone i . It consists in sensing information about the target's type that the drone i is overflying. The action sample samples the considered target and changes the drone i 's state to going-to-base state with a probability of 1.0. The action go-to-base moves the drone i to the going-to-base state with a probability of 1.0, except if the drone i is already on the going-to-base state. In this case it moves the drone i to the base state with a probability of 0.20, and remains in going-to-base state with 0.8.

3) Effects of actions related to operator's request: the action show-target j depends on the availability state of the human operator. The probability effect of this action depends on the time taken by the human operator to answer the drones' request when available. To ask for the operator's attention, the ground station performs a visual effect to emphasize where the operator should look (e.g. by launching a yellow bounding box, see Fig. 1). Previous work has shown that the yellow bounding box visual effect that is used on the Search & Rescue task is very efficient, taking around 3 seconds in average to be seen [START_REF] Imbert | Attentional costs and failures in air traffic control notifications[END_REF]. Based on it, we consider that the operator would answer to the interface and so, changes its state to at-target j with a probability of 0.33 if available, and would remain on the previous state with the residual probability. The action show-other-task has no visual effects on the application. When the operator is considered as not available all actions have no effect.

4) Effects of drones' actions on target states: when drones are moving, target states do not change. But, when a drone samples a target, this target can change its type with an uniform probability over types. This allows us to model the fact that new targets can be integrated into the model, appearing anywhere on the application, once the actual one has been sampled. Also, when all the drones are on the base, all targets' types can change with an uniform probability. This trick allows us to exploit longer missions without reinitializing the planning problem.

5) The availability of the operator: the operator must perform the drones' monitoring task in priority. But, he can be requested when he is estimated available, e.g. when he has finished the other-task (Short-Term Item Memorization task) for instance. The other-task has a recurrent time duration of 30 seconds in average, with a 13 seconds cross fixation period which is considered as a break where the human operator could be estimated as available. Therefore, the human operator has a probability of 0.08 to switch from KO to OK, which is a worst-case probability, and vice-versa. drone. An informed observation can be done by a drone (resp. the human) only if it is over a target and only if it performs a check action (resp. only when he answers a show-target action request). The observation function is summarized in Table I. These probabilities were arbitrarily chosen for this study case. 7) Reward function: all moving actions have the same cost, 0.1 per drone, except if at least one drone is on the base (it is discounted from the cost). The check action has the cost of 5 per drone. We consider that the payoff related to a sample action is superior for human targets than for material ones (human 1000, material 300). If the drone samples a non relevant target it costs 10000 per target sampled. Note that if two drones sample the same target the payoff remains the same. It prevents multiple drones to sample the same target. The action sample is penalized if it is performed elsewhere than while flying a target. If a show-target request action is performed, it costs 2.5 if the operator is available, and 20 if he is not available.

8) Off-line policy generation: for this study case, and based on the parametric model before presented, we have generated the MOMDP model considering 2 drones and 3 targets. The model obtained was solved using the SARSOP algorithm [START_REF] Ong | Planning under uncertainty for robotic tasks with mixed observability[END_REF], with = 0.1 and with a memory limited to 6.5Gb (a computer with 8Gb RAM). The computed policy was then executed on the experimental platform described hereafter to handle the collaborative interaction between human and artificial agents.

IV. MATERIALS AND METHODS

A. Volunteers

Ten male volunteers (mean age: 33.0, STD= 8.1) completed the experiment. They have normal or corrected-tonormal vision, normal hearing, and no psychiatric disorders. After providing written informed consent, they were instructed to complete the experimental scenario (Cf. Section II).

B. Global System Architecture

The global system architecture is presented in Fig. 2. The interaction between the operator and the AAVs is ensured by the ATMOSPHEr1 interface [START_REF] Collart | Human-robot systems facing ethical conflicts: a preliminary experimental protocol[END_REF]. This interface can be seen as part of a ground station which is equipped with an eye-tracker sensor. 1) Eye-tracking Device: An eye-tracker (SMI RED250) was used to collect the human operator's eye gaze. Volunteers were asked to stay at a distance of around 60 cm from the screen. The calibration tolerance was set to be horizontally and vertically under 0.5 • . The refreshment rate of the eyetracker was set to 60Hz. We used SMI iView X T M Software Development Kit 2 to retrieve the volunteer's gaze position on the screen, and to export it through a Python TCP/IP interface that communicates with the Data Fusion module (Fig. 2). We used the participants' eye gaze as a straightforward indicator of their cognitive availability. In a nutshell, the discrete state variable availability (Cf. Section III-B.1) for the request of the decisional policy was set to OK when the participant gaze was on the left part of the screen (looking at Search & Rescue task), and KO otherwise (looking at Short-Term Item Memorization task or outside the screen, or if the gaze signal was lost).

C. Metrics 1) Performance evaluation on Search & Rescue task:

We evaluated the operator's Search & Rescue task performance with he ratio of the number of answered requests over the total number of requests asked by the policy. We were interested in how much the human operator was engaged on this task, that was defined as primary. We also logged the expected rewards of the human-robots system during experiments (MOMDP score given by the SARSOP policy simulator) as an indicator of policy execution efficiency (Cf. Section III-B.7).

2) Perfomance evaluation on Short-Term Item Memorization task: Each participant had the same number of digits series during each session. Therefore, we defined the performance criteria of Short-Term Item Memorization task 2 lyrawww.uvt.nl/ cenv/dci-lab/smi/iViewX.pdf as the number of digits series that has been fully retrieved, in the correct order.

D. Protocol's Sequence

First, the volunteer was trained on the Short-Term Item Memorization task alone (duration: 10 min), and then on the Search & Rescue task alone (duration: 10 min). Afterwards, the volunteer passed two times the dual task protocol:

• in one condition, the decisional system was using the eye gaze position measured by the eye-tracking device: we will call it the "with gaze monitoring" (with GM) condition (duration: 10 min); • in the other condition, the decisional system had no information about where the volunteer was looking, and it was considering him as "available": this was the "without gaze monitoring" (without GM) condition (duration: 10 min) and was used as a control condition.

The order between the with GM condition and the without GM condition was counter-balanced between the volunteers. The volunteers had no information about the type of condition they were achieving. At the end of each dual task condition, the volunteer completes a NASA-TLX form [START_REF] Hart | Development of NASA-TLX (task load index): Results of empirical and theoretical research[END_REF].

The typical total duration of a subject's session (informed consent approval, the two practice tasks, and the two dual tasks) was about one hour and fifteen minutes.

V. EXPERIMENTAL RESULTS

A. The Search & Rescue task's Scores

The Search & Rescue task's scores are summarized in Table II. One can observe that there is a tendency concerning a better request ratio score on with GM condition than on without GM condition. However, as the volunteers' request ratio score is not normally distributed, non-parametric Wilcoxon rank-sum tests were applied on median values for each condition. These tests reveal that ratio scores are not significantly different between the with GM condition, the without GM condition and the Search & Rescue task achieved alone. Based on this, we can not formally advance any conclusion with respect to this score. It is also worth noting that the variance was the largest in the without GM condition suggesting that the volunteers had more contrasted behavior regarding task prioritization. The strategy in the with GM condition proposed by the system led to less variability reflecting a more balanced behavior of human operators.

In the same Table, the average expected rewards were also summarized. Given that these expected rewards were averaged between 10 volunteers, this is obviously not enough to evaluate the performance of the MOMDP policy. Paired ttests did not reveal any significant difference concerning final expected reward (human-robots system) between the three conditions. Nevertheless, average expected rewards have the same order of magnitude. It means that even when the human operator was not available to help the AAVs (with GM condition), the robotic system could achieve the target recognition task alone.

B. Short-Term Item Memorization task's Scores

The Short-Term Item Memorization task's scores are summarized in Figure 3. For this scores, parametric two-tailed paired t-tests were applied on mean values for each condition. P-values were corrected for multiple comparisons (Bonferroni-Holm). Paired t-test results revealed that the volunteers perform significantly better on Short-Term Item Memorization task during the with GM condition than in the without GM condition (p corrected <0.01) (Fig. 3). They also performed significantly better during the Short-Term Item Memorization task alone training session than in the with GM condition (p corrected <0.001).

Without Gaze Monitoring

With

C. Volunteers' Subjective Feedback

A two tailed paired t-test was conducted on the NASA-TLX scores and revealed that the volunteers rated the with GM condition's workload statistically significantly lower than the without GM condition's workload (p<0.01). These results show the importance of taking into account a measure of the human self-estimated workload while interacting with robots and performing many tasks.

D. Synthesis

The behavioral measures reveal that the decision system including a cognitive state indicator (with GM condition) helps the participant to maintain a good performance on the primary drones monitoring task (Search & Rescue task), while improving the score on the secondary WM task (Short-Term Item Memorization task). Consistently, the subjective rating confirmed that with GM condition reduced the volunteer's workload. Our results confirm previous findings [START_REF] Parasuraman | Humans: Still vital after all these years of automation[END_REF] supporting that taking into account the human operator in the design of the system leads to better overall performance. In our context, the overall performance is the performance of the whole system that includes human-robots and all the tasks the human has to perform.

VI. CONCLUSION AND FUTURE WORKS

The objective of this study was to show that the (MO)MDP decision framework is suitable to design a human-robots interactive system. We implemented such a framework to handle the non-deterministic nature of the human operator's cognitive availability.

The results of this study supports the advantages of taking into consideration the cognitive availability of the human operator in such complex contexts and to integrate psychophysiological measurement: volunteers had better objective performances in the double task paradigm when their cognitive availability state was taken into account. The subjective measurements were also in line with these latter results as the volunteers experienced lower mental workload in the with GM condition than in the without GM one. Taken together these results demonstrate that the human is not a providential agent and that he can not take over or assist the robots at anytime.

The MOMDP formalism is suited to model probabilistic actions' effects, as the human's time to answer a request when solicited, and probabilistic observations as the misclassification rate when identifying targets. Moreover, such a formalism can be used to consider the human cognitive state as a partial observable state (see [START_REF] Nikolaidis | Efficient model learning from joint-action demonstrations for human-robot collaborative tasks[END_REF], [START_REF] Souza | MOMDPbased target search mission taking into account the human operator's cognitive state[END_REF] for instance). More interestingly, using an adapted decision model framework, such as a MOMDP, we can directly integrate confidence weight on the observation of the cognitive state we are obtaining. Though we admit that our metric to derive the cognitive availability was simple, this proof of concept demonstrates that it is easy to introduce uncertainties on observations concerning this state variable while maintaining the same framework. Therefore, we advocate for the suppleness of such a MOMDP formalism, that allows to easily extend the model to consider new state variables and observation variables. This approach paves the way to integrating more complex partial observable human cognitive states, such as mental workload, mind wandering state, and task engagement. Improvements can also be made concerning the model used (Cf. Section III): for example, the generation of the probability table (Cf. Section III-B.6) could be improved by using learned data. If probabilities tables are hard to be learned, one could consider an extension of the MOMDP model, the Qualitative Possibilistic MOMDPs [START_REF] Drougard | Qualitative possibilistic mixed-observable mdps[END_REF], where the possibility theory is applied instead of a probability one.

Fig. 1 .

 1 Fig. 1. Screen-shot of the participant's view of the dual task experimental setup.

  b o,a,x Y (y ) = η y ∈Y p(o|y , x , a)p(y |x, y, a, x )p(x |x, y, a)bY (y)

•

  a-dronei: {goto-target1...goto-targetnt, sample, check, base}, ∀i ∈ [[1; nd]]; • a-operator-request: {show-target1...show-targetn t , showother-task}.

Fig. 2 .

 2 Fig. 2. The complete system architecture. The human operator is able to interact with the Search & Rescue task and the Short-Term Item Memorization task interfaces, while an eye-tracker is monitoring his eyes position. The ATMOSPHEr interface sends to the Data Fusion module the current drones' state and observations, and the human operator's input (sensing task) with a frequency of 1Hz (the same clock used on the MOMDP belief update and action request too). The Data Fusion module sends to the MOMDP policy module the compiled couple (State, Observation). The MOMDP module updates the belief state and sends the next Action to the ATMOSPHEr interface.

Fig. 3 .

 3 Fig. 3. Mean number of fully correct digits on the three conditions (without GM, with GM and training on Short-Term Item Memorization task alone), on 20 trials. Red bars indicate standard errors and *: p corrected < 0.01, **: p corrected < 0.001.

  [START_REF] Woods | Learning from automation surprises and" going sour" accidents[END_REF] Observation function: the observation model is based on the fact that all recognition tasks performed by a human operator are still more reliable than when performed by a

	state/obs	ohuman	Drone (Operator) omaterial	onone
	human	0.65 (0.95) 0.30 (0.05) 0.05 (0.05)
	material 0.20 (0.05) 0.75 (0.95) 0.05 (0.05)
	none	0.05 (0.05) 0.05 (0.05) 0.90 (0.95)

TABLE I DRONES

 I (AND HUMAN) OBSERVATIONS PROBABILITY TABLE.

A Tiny Multi-agent system Oriented Simulator Platform for HMI Experiments -https://sourceforge.isae.fr/projects/atmospher/