In the study that follows, we consider a dataset consisting of the x and y-coordinates, calculated from the GPS longitudes and latitudes, representing the recorded position of the robot at time t: D = x t , y t , t , t ∈ T .

These coordinates span six different disjoint geographical patches within the whole vineyard, as illustrated in Figure 2, that have been recorded in a period of time T leading to the maps reported in Figure 3 that reflects the robot’s trajectory.

Maintenance operations are also known and properly identified in the provided dataset. Thus, the dataset consists of a collection of n discrete, finite and compact two-dimensional trajectories S 1 , … , S n .

We are interested in extracting the geometrical and topological features of the trajectories in D across different scales. For that purpose, we introduce the so-called Rips filtration. We construct a Rips complex from simplexes of varying dimensions that are generalizations of triangles of varying dimensions. More specifically, a d-simplex is the smallest convex set of d + 1 points, x ₀, … , x _d where x ₁−x ₀, … , x _d −x ₀ are linearly independent, as illustrated in Figure 4. The so-called abstract simplicial complex is a finite collection of sets that is closed under the subset relation, i.e., if a ∈ A and b ⊂ a, then b ∈ A.

Let S be a trajectory, defined from a finite compact set of points in R 2 , and ϵ ≥ 0. The Rips complex of S at scale ϵ, R ϵ ( S ) , is the abstract simplicial complex consisting of all subsets of diameter up to ϵ: R ϵ S ≔ σ ⊂ S | diam σ ≤ ϵ , where the diameter of a set of points is the maximum distance between any two points in the set.

Geometrically, we can construct the Rips complex by considering balls of radius ϵ 2 , centred at each point in S . Whenever d balls have pairwise intersections, we add a d−1 dimensional simplex. An example of Rips complex is given in Figure 5.

A filtration of a simplicial complex K is a nested sequence of subcomplexes starting at the empty set and ending with the full simplicial complex ∅ ⊂ K 0 ⊂ ⋯ ⊂ K .

By varying the value of the scale parameter ϵ, from ϵ _min = 0 to ϵ max = diam ( S ) we get a family of nested Rips complexes known as the Rips filtration.

In order to have a more exhaustive view on how the features are changing across different scales, the appearance and disappearance of each feature within the filtration is tracked and coded into the homology groups H k ( S ) , where k is the homology dimension. The elements of a Homology Group H k ( S ) are classes of chain of simplexes (“packets”) in the Rips complex. The use of homology groups allows us to perform algebraic operations over the simplicial elements. The homology group H 0 ( S ) represents the vertices, while the homology group H 1 ( S ) represents the cycles (loops) formed in the simplicial complex. Since our data is in R 2 we are only interested in k = 0 and k = 1.

Given a homology group, we can now define how to track the appearance of the features across different scales, by defining the homology group at a scale ϵ, H k ϵ ( S ) . It represents the classes of simplexes as described previously, but taken from R ϵ ( S ) . That is, the elements of R ϵ ( S ) with a filtration value lower than ϵ. This approach is known as the persistent homology. It allows to quantify the appearance and disappearance of the features across the different scales (discretized by considering m values related to ϵ _j , j = 0, … , m):

• For H 0 ( S ) , the birth scale of all vertices is set to zero, while the death scale is the filtration value at which the vertex has been joined to another one by a segment.

• For H 1 ( S ) , the birth scale of a cycle is the filtration value at which a loop has been formed, while the death scale is the filtration value at which the interior of the loop has been covered.

We can formalize this as follows:

• The birth scale b _γ of the feature γ

The persistence of the features throughout the scales can then be represented by the so-called persistence barcode of S . It is a histogram, where the bar associated to each feature starts at the birth scale and ends at the death scale.

An example of persistent homology computation is given with the rips complex in Figure 6, and the associated barcode in Figure 7. A loop is formed at ϵ = 0.9 (birth) and then covered at ϵ = 1.8 (death). It is represented by the red bar.

A more compact representation of the features persistence is the persistence diagram of S , defined from PD S = b γ , d γ : γ ∈ H k , where b _γ and d _γ are the birth and death scales associated to the feature γ. In what follows, in the trajectories analysis, we only consider one-dimensional features, i.e., k = 1.

The persistence diagram associated with the Rips complex shown in Figure 6 is given in Figure 8. An equivalent representation of the persistence diagram consists in the so-called life-time diagram of S , which is constructed by means of a bijective transformation T (a, b) = (a, b−a), acting over PD ( S ) , that is, LT S ≔ a , b − a ∈ R 2 : a , b ∈ PD S .

In order to use the persistence features in a machine learning approach, we construct the so-called persistent image of S . First, observe that LT ( S ) is a finite set of p points, LT S = a 1 , b 1 − a 1 , … , a p , b p − a p , and such that b ₁−a ₁ ≤ b ₂−a ₂ ≤ …, ≤ b _p −a _p . Then, consider a non-negative weighting function given by w : LT S → 0,1 a i , b i − a i ↦ w a i , b i − a i = b i − a i b p − a p , for 1 ≤ i ≤ p .

Finally, we fix M, a natural number, and take a bivariate normal distribution g _u (x, y) centered at each point u ∈ LT ( S ) with a variance σ I 2 = b p − a p M I 2 I 2 is the 2 × 2 identity matrix). A persistence kernel is then defined according to: ρ S : R 2 → R x , y ↦ ρ S x , y = ∑ u ∈ LT S w u g u x , y .

We associate to a robot trajectory S ∈ R 2 a matrix in R M × M as follows: let δ > 0 be a non-negative, small enough real number, and then consider a squared region Ω S , δ = [ a , b ] × [ c , d ] ⊂ R 2 , covering the support of ρ S ( x , y ) up to a certain precision δ, such that ∬ Ω S , δ ρ S x , y d x d y ≥ 1 − δ .

Then, we consider two uniform partitions of the intervals a = p 0 ≤ p 1 ≤ … , ≤ p M = b  and  c = q 0 ≤ q 1 ≤ … , ≤ q M = d .

Finally, we express Ω S , δ from Ω S , δ = ⋃ i = 0 M − 1 ⋃ j = 0 M − 1 p i , p i + 1 × q j , q j + 1 = ⋃ i = 0 M − 1 ⋃ j = 0 M − 1 P i j .

The persistence image of S associated with the partition P = { P i j } is then described by the R M × M matrix with elements: P I S , M , P , δ i j = ∬ P i j ρ S x , y d x d y  for  0 ≤ i , j ≤ M − 1 .

An example of persistence computation for a given trajectory is given in Figure 9.

Consider two data sets S u and S v representing two trajectories. A matching between two persistence diagrams, PD ( S u ) and PD ( S v ) , is a map ψ, that reads: ψ : PD S u → PD S v , such that ∀ γ = ( b , d ) ∈ PD ( S u ) , ψ γ = ψ 1 b , ψ 2 d = b ′ , d ′ ∈ PD S v .

The map ψ associates each feature from PD ( S u ) to a feature from PD ( S v ) . The optimal matching between PD ( S u ) and PD ( S v ) is a matching ψ ̂ ψ ̂ : PD S u → PD S v , minimizing the transport cost C to move the features from PD ( S u ) to PD ( S v ) : C min = ∑ γ ∈ PD S u ‖ γ − ψ ̂ γ ‖ 2 = ∑ b , d ∈ PD S u ‖ b − ψ ̂ 1 b , d − ψ ̂ 2 d ‖ 2 = ∑ b , d ∈ PD k S u b − ψ ̂ 1 b 2 + d − ψ ̂ 2 d 2 .

Then, to measure the degree of similarity between two trajectories S u and S v we consider the Wasserstein distance (Villani, 2006; Peyré and Cuturi, 2019) between PD ( S u ) and PD ( S v ) W PD S u , PD S v = ∑ b , d ∈ PD S u b − ψ ̂ 1 b 2 + d − ψ ̂ 2 d 2 , where ψ ̂ is the optimal matching between PD ( S u ) and PD ( S v ) .

An example of matching between the persistence diagrams of two trajectories is given in Figure 10.

Consider now a collection S 1 … S n of trajectories with their associated diagrams PD 1 … PD n .

Since the space of persistence diagrams equipped with the Wasserstein distance, the Wasserstein space, is not a linear space, the notion of barycentres (Agueh and Carlier, 2011) can be extended for the persistence diagrams using the so-called Frechet mean (Turner et al., 2014), which always exists in the context of averaging finitely many diagrams.

The Frechet mean of PD 1 … PD n is any diagram minimizing the map E : μ ↦ ∑ i = 1 n W μ , PD i 2 .

The computation of the barycentre μ has proven to be challenging, and multiple approaches can be used, such as the Sinkhorn algorithm (Cuturi and Doucet, 2014). We will use the one based on the Hungarian algorithm presented in Turner et al. (2014) and consider Partial Optimal Matchings (Divol and Lacombe, 2020), as the diagrams may not be of the same size. In this case, points from the diagonal are matched with the remaining (exceeding) points.

In our case, we estimate the barycentres of a finite family of persistence diagrams, taking a Lagrangian approach by tracking the individual points of the diagrams. Given a collection PD 1 … PD n of persistence diagrams, we proceed as follows:

1) Initialize the estimation μ of the barycenter at a certain diagram μ = PD i 0 .

4. If μ ̂ minimizes E , return μ ̂ . Otherwise, update μ = μ ̂ and go back to 2.

An example of a barycentre of three persistence diagrams is given in Figure 11.

Image classification is a procedure that is used to automatically categorize images into classes by assigning to each image a label representative of its class. A supervised classification algorithm requires a training sample for each class, that is, a collection of data points whose class of interest is known. Labels are assigned to each class of interest. The classification problem applied to a new observation (data) is thus based on how close a new point is to each training sample. The Euclidean distance is the most common metrics used in low-dimensional datasets. The training samples are representative of the known classes of interest to the analyst. In order to classify the persistence images, we considered the logistic regression algorithm.

Consider a training set X i i = 1 n of flattened persistence images, i.e., M × M-component vectors, computed from a set S i i = 1 n of trajectories as described earlier. Associated is a list ( Y i ) i = 1 n of binary labels {0, 1}, describing whether an image X i is in the interest set or not.

The training of the L 2 -penalized logistic regression binary classifier is then the minimization of a cost function as described in the following optimization problem: min ω , c 1 2 ω T ω + C ∑ i = 0 p log exp Y i X i T ω + c + 1 .

Here ω are the weights we optimize over, c a Bernoulli mean vector of the weights, and C an inverse regularization parameter. Once trained, the model is evaluated on a unseen set of flattened persistence images. The metric used for the model evaluation is the Accuracy Score defined in the next section.

Evaluating a classification model consists of determining how often labels are correctly or wrongly predicted for the testing samples. In other words, it is counting how many times a sample is correctly or wrongly labelled into a particular class. We distinguish four qualities:

• TP (True Positive): the correct prediction of a sample into a class;

These quantities are involved in the definition of the model performances estimator, the AccuracyScore (A). It is giver by the ratio of the number of correct predictions over the number of all samples, expressed by A = T P + T N T P + T N + F P + F N

We first want to predict whether a robot is in a certain patch. For that purpose we choose one parcel as a target, and train a classification model as described in Section 2.6. The complete dataset consists of daily trajectories for 240 days. For each day a persistence image is computed, which will then be used as input for the model (a sample is depicted in Figure 9). The samples are labelled according to the target patch (patch #3): 1 if the robot is in the target patch (114 samples), and 0 otherwise (126 samples). The dataset is split into 65% for training and 35% for testing. The proposed classifier achieves an 80% accuracy score in predicting the patch at which the robot is, based on the persistence images.

Then, we consider daily trajectories in the same patch (patch #3), consisting of 50 samples. For each day, a persistence image is computed, that will be used as input in the classifier. The periods considered here are the ones in between two consecutive maintenance operations of the robot. The samples are labelled 0 if they are associated to a day before the maintenance date (25 samples), 1 otherwise (25 samples). The dataset is split into 65% for training and 35% for testing. The model achieves a 90% accuracy score predicting the period associated to the sampled trajectories. The model high accuracy proves that the topological descriptors have enough information about the pathways to allow detecting patterns related to maintenance events, fact that could be used for predictive maintenance purposes.

Figure 14 depicts the Wasserstein distance between the persistence diagrams for consecutive daily trajectories, with the maintenance operation emphasized in red, whereas Figure 15 shows the barycentres of each period between consecutive maintenance operations. As it can be noticed from the persistence images in Figure 15, maintenance operations affect the topology of the trajectory, as it was expected from the fact that classification performs successfully as just reported.

To better support our hypothesis about the effect of maintenance on the trajectory topology, we consider the first operation interval, the one before the first maintenance, that correspond to the first persistence image in Figure 15 (left), and divide it in two parts with identical length. Then, the associated barycentres in both half intervals are obtained. Both are represented in Figure 16.

As it can be noticed, both of them resemble very much to the one associated to the whole interval (the first picture in Figure 15), with a Wasserstein distance of 23.5 and 24.1 (computed on the associated diagrams). Conversely, the distance of both to the second period is much higher (46.5 and 26.1).

Finally, we can compute the distance of the 5 latter periods to the first one (Figure 15) and we have: 36.7, 31.0, 34.2, 33.3, 35.2, so significantly higher than the first period compared to its own two halves.

These results support again our assumption on the effect of maintenance on the trajectory topology.