Julien Alexandre
email: alexandre@ensta-paristech.fr

Alexandre Sandretto

Olivier Chapoutot
email: chapoutot@ensta-paristech.fr

Mullier
email: mullier@ensta-paristech.fr

Formal Verification of Robotic Behaviors in Presence of Bounded Uncertainties

Robotic behaviors are mainly described by differential equations. Those mathematical models are usually not precise enough because of inaccurately known parameters or model simplifications. Nevertheless, robots are often used in critical contexts as medical or military fields. So, uncertainties in mathematical models have to be taken into account in order to produce reliable and safe analysis results. A framework based on interval analysis is proposed to safely verify and analyze robotic behaviors with bounded uncertainties. It follows an interval constraint programming approach, combined with validated numerical integration methods to deal with differential equations. A case study on robust path planning is presented to emphasize the efficiency of the complete framework.

Introduction

Designing control algorithms of robots is a challenging activity in order to guarantee the success of the mission and taking into account safety constraints. Moreover, control architecture is made of several layers strongly interacting with each other which makes the design complex. Indeed, several kinds of algorithms ranging from Proportional-Integral-Derivative controllers (PID controllers) to neural networks are used and combined to achieve the missions of robots.

In control theory, mathematical models of the controlled element are essential to define control algorithms. Nevertheless, such models are usually a simplification of the real world. Mainly, accurate models are expensive to define, or even impossible, due to the presence of unknown, or hardly measurable, values of the involved parameters (e.g., the inertial moment). Mathematical models represent possible behaviors of robots in function of the control algorithms. So they play a critical role in the design or in the verification process of robots.

The framework given by constraint programming [START_REF] Rueher | Solving continuous constraint systems[END_REF] is powerful enough to solve many problems with numerous constraints. Unfortunately, taking into account the dynamic of robots is not easy in this framework. In contrary, validated numerical integration methods are mature techniques to simulate, with bounded uncertainties, dynamical systems described by differential equations [START_REF] Nedialko | Validated solutions of initial value problems for ordinary differential equations[END_REF][START_REF] Rauh | A novel interval arithmetic approach for solving differential-algebraic equations with ValEncIA-IVP[END_REF][START_REF] Alexandre | Validated Explicit and Implicit Runge-Kutta Methods[END_REF][START_REF] Alexandre | Validated Simulation of Differential Algebraic Equations with Runge-Kutta Methods[END_REF]. While such methods produce rigorous enclosures of the trajectories, they can not consider constraints on the dynamics without major modifications. Our goal is to bridge the gap between these two approaches.

In this article, we propose an approach to analyze or verify control algorithms of robots based on mathematical models of their behaviors. An extension of constraint programming framework with differential equations is proposed, hence our models can take into account bounded uncertainties on the parameters offering guarantees on robot behaviors, while taking into account temporal constraints as safety. We extend current work on this field, such as [START_REF] Janssen | Multistep filtering operators for ordinary differential equations[END_REF][START_REF] Cruz | Constraint reasoning over differential equations[END_REF][START_REF] Goldsztejn | Including ordinary differential equations based constraints in the standard CP framework[END_REF], by considering more complex constraints involving time.

Formal verification of robotic application is not new and several previous work address this problem. Nevertheless, they usually rely on temporal logic specification as Linear Temporal Logic (LTL) to prove the correctness. For example, this approach is used in [START_REF] Hoxha | Planning in dynamic environments through temporal logic monitoring[END_REF][START_REF] Liu | Synthesis of reactive switching protocols from temporal logic specifications[END_REF]. Satisfiability Modulo Theories techniques (SMT techniques), which are described in [START_REF] Eggers | SAT modulo ODE: A direct SAT approach to hybrid systems[END_REF][START_REF] Gao | dReal: An SMT solver for nonlinear theories over the reals[END_REF], can also be used to prove safety or even stability of hybrid systems. Our contribution differs from these works by adopting constraint satisfaction problem formulation in order to easily specify the requirements.

The paper is organized as follows. In Section 2, the basics of interval analysis, constraint satisfaction problem and guaranteed numerical integration are provided. The main contribution is presented in sections 3 and 4. A case study is described in Section 5 before concluding in Section 6.

Preliminary Notions

A presentation of the main mathematical tools is given in this section. First, basics of interval analysis is provided in Section 2.1. Second, Section 2.2 gives a quick presentation of constraint satisfaction problem. Finally, a short introduction of validated numerical integration is presented in Section 2.3.

Interval Analysis

The simplest and most common way to represent and manipulate sets of values is interval arithmetic (see [START_REF] Moore | Interval Analysis[END_REF]). An interval [x i ] = [x i , x i ] defines the set of reals x i such that x i ≤ x i ≤ x i . IR denotes the set of all intervals over reals. The size or the width of [x i ] is denoted by w([

x i ]) = x i -x i .
Interval arithmetic extends to IR elementary functions over R. For instance, the interval sum, i.e.,

[x 1 ] + [x 2 ] = [x 1 + x 2 , x 1 + x 2 ]
, encloses the image of the sum function over its arguments.

An interval vector or a box [x] ∈ IR n , is a Cartesian product of n intervals. The enclosing property basically defines what is called an interval extension or an inclusion function.

Definition 1 (Inclusion function) Consider a function

f : R n → R m , then [f ] : IR n → IR m is said to be an extension of f to intervals if ∀[x] ∈ IR n , [f ]([x]) ⊇ {f (x), x ∈ [x]} .
It is possible to define inclusion functions for all elementary functions such as ×, ÷, sin, cos, exp, and so on. The natural inclusion function is the simplest to obtain: all occurrences of the real variables are replaced by their interval counterpart and all arithmetic operations are evaluated using interval arithmetic. More sophisticated inclusion functions such as the centered form, or the Taylor inclusion function may also be used (see [START_REF] Jaulin | Applied Interval Analysis[END_REF] for more details). In the first example of division, the result is the interval containing all the real numbers because denominator contains 0.

Example 1 (Interval arithmetic) A few examples of arithmetic operations between interval values are given

As an example of inclusion function, we consider a function p defined by p(x, y) = xy + x .

The associated natural inclusion function is

[p]([x], [y]) = [x][y] + [x],
in which variables, constants and arithmetic operations have been replaced by its interval counterpart. And so

p([0, 1], [0, 1]) = [0, 2] ⊆ {p(x, y) | x, y ∈ [0, 1]} = [0, 2].

Numerical Constraint Satisfaction Problems

In this section, we introduce the numerical constraint satisfaction problem (NCSP) formalism, following the description given in [START_REF] Rueher | Solving continuous constraint systems[END_REF], and present some basics on constraint programming. A NCSP (V, D, C) is defined as follows:

• V := {v 1 , . . . , v n } is a finite set of variables which can also be represented by the vector v; • C := {c 1 , . . . , c m } is a set of constraints of the form c i (v) ≡ f i (v) = 0 or c i (v) ≡ g i (v) 0, with f i : R n → R, g i : R n → R for 1 i m. Constraints C are interpreted as a conjunction of equalities and inequalities.

• D := {[v 1 ], . . . , [v n ]}
The solution of a NCSP is a valuation of v ranging in [v] and satisfying the constraints C.

The approach of NSCP is both powerful to address complex problems (NP-hard problem with numerical issues, even in critical applications) and simple in the definition of a solving framework [START_REF] Lebbah | Accelerating filtering techniques for numeric CSPs[END_REF][START_REF] Benhamou | Clp (intervals) revisited[END_REF]. Indeed, the classical algorithm to solve a NCSP is the branch-and-prune which needs only an evaluation of the constraints and an initial domain for variables. If this algorithm is sufficient for many problems, some improvements have been achieved and the algorithms based on contractors have emerged [START_REF] Chabert | Contractor programming[END_REF]. The branch-and-contract algorithm consists in two main steps i) the contraction (or filtering) of one variable and the propagation to the others till a fixed point is reached, then ii) the bisection of the domain of one variable in order to obtain two problems, easier to solve. A more detailed description follows.

Contraction

A filtering algorithm or contractor is used in a NCSP solver to reduce the domain of variables till a fixed point (or close to), by respecting the local consistencies. A contractor Ctc can be defined with the help of constraint programming, analysis or algebra, but it must satisfy three properties:

• Ctc(D) ⊆ D: the contractance,

• Ctc cannot remove any solution: it is conservative,

• D ⊆ D ⇒ Ctc(D ) ⊆ Ctc(D): the monotonicity.
There are many contractor operators defined in the literature, among them the most notable are:

• (Forward-Backward contractor) By considering only one constraint, this method computes the interval enclosure of a node in the tree of constraint operations with the children domains (the forward evaluation), then refines the enclosure of a node in terms of parents domain (the backward propagation). For example, from the constraint x + y = z, this contractor refines initial domains • (Newton contractor) This contractor, based on first order Taylor interval extension:

[f ]([x]) = f (x * ) + [J f ]([x])([x] -x * ) with x * ∈ [x]
, has the property of

-if 0 ∈ [f ]([x]), -then [x] k+1 = [x] k ∩ x * -[J f ]([x] k ) -1 f (x * ) is a tighter inclusion of the solution of f (x) = 0.
Some other contractors based on Newton have been defined such as Krawczyk [START_REF] Jaulin | Applied Interval Analysis[END_REF].

Propagation

If a variable domain has been reduced, the reduction is propagated to all the constraints involving that variable, which allows to narrow the other variable domains. This process is repeated till a fixed point is reached.

Branch-and-Prune

A Branch-and-Prune algorithm consists on alternatively branching and pruning to produce two sub-pavings L and S, with L the boxes too small to be bisected and S the solution boxes. We are then sure that all solutions are included in L ∪ S and that every point in S is solution.

Literally, this algorithm browses a list of boxes W, initially W is set with a vector [x] made of the elements of D, and for each one i) Prune: the NCSP is evaluated (or contracted) on the current box, if it is a solution it is added to S, otherwise ii) Branch: if the box is large enough, it is bisected and the two boxes resulting are added into W, otherwise the box is added to L.

Example 2 An example of the problems that the previously presented tools can solve is taken from [START_REF] Lhomme | Consistency techniques for numeric CSPs[END_REF]. The considered NCSP is defined such that

• V = {x, y, z, t} • D = {[x] = [0, 1000], [y] = [0, 1000], [z] = [0, 3.1416], [t] = [0, 3.1416]} • C = {xy + t -2z = 4; x sin(z) + y cos(t) = 0; x -y + cos 2 (z) = sin 2 (t); xyz = 2t}
To solve this CSP we use a Branch-and-Prune algorithm with the Forward-Backward contractor and a propagation algorithm. We obtain one solution ([1.999, 2.001], [1.999, 2.001], [1.57, 1.571], [3.14159, 3.1416]). The algorithm needs 6 bisections.

Validated Numerical Integration Methods

A differential equation is a mathematical relation used to define an unknown function by its derivatives. These derivatives usually represent the temporal evolution of a physical quantity. Because the majority of systems is defined by the evolution of a given quantity, differential equations are used in engineering, physics, chemistry, biology or economics as examples. Mathematically, differential equations have no explicit solutions, except for few particular cases. Nevertheless, the solution can be numerically approximated with the help of integration schemes such as Taylor series [START_REF] Nedialko | Validated solutions of initial value problems for ordinary differential equations[END_REF] or Runge-Kutta methods [START_REF] Alexandre | Validated Simulation of Differential Algebraic Equations with Runge-Kutta Methods[END_REF][START_REF] Alexandre | Validated Explicit and Implicit Runge-Kutta Methods[END_REF].

In the following, we consider a generic parametric differential equation as an interval initial value problem (IIVP) defined by

     ẏ = F (t, y, x, p, u) 0 = G(t, y, x, p, u) y(0) ∈ Y 0 , x(0) ∈ X 0 , p ∈ P, u ∈ U, t ∈ [0, t end ] , (1)
with

F : R × R n × R m × R r × R s → R n and G : R × R n × R m × R r × R s → R m .
The variable y of dimension n is the differential variable while the variable x is an algebraic variable of dimension m with an initial condition y(0) ∈ Y 0 ⊆ R n and x(0) ∈ X 0 ⊆ R n . In other words, differential-algebraic equations (DAE) of index 1 are considered, and in the case of m = 0, this differential equation simplifies to an ordinary differential equation (ODE). Note that usually, the initial values of algebraic variable x are computed by numerical algorithms used to solve DAE but we consider it fixed here for simplicity. Variable p ∈ P ⊆ R r stands for parameters of dimension r and variable u ∈ U ⊆ R s stands for a control vector of dimension s. We assume standard hypotheses on F and G to guarantee the existence and uniqueness of the solution of such problem.

A validated simulation of a differential equation consists in a discretization of time, such that t 0 • • • t end , and a computation of enclosures of the set of states of the system y 0 , . . . , y end , by the help of a guaranteed integration scheme. In details, a guaranteed integration scheme is made of • an integration method Φ(F, G, y j , t j , h), starting from an initial value y j at time t j and a finite time horizon h (the step-size), producing an approximation y j+1 at time t j+1 = t j + h, of the exact solution y(t j+1 ; y j ), i.e., y(t j+1 ; y j ) ≈ Φ(F, G, y j , t j , h);

• a truncation error function lte Φ (F, G, y j , t j , h), such that y(t j+1 ; y j ) = Φ(F, G, y j , t j , h)+lte Φ (F, G, y j , t j , h).

Basically, a validated numerical integration method is based on a numerical integration scheme such as Taylor series [START_REF] Nedialko | Validated solutions of initial value problems for ordinary differential equations[END_REF] or Runge-Kutta methods [START_REF] Alexandre | Validated Simulation of Differential Algebraic Equations with Runge-Kutta Methods[END_REF][START_REF] Alexandre | Validated Explicit and Implicit Runge-Kutta Methods[END_REF] which is extended with interval analysis tools to bound the truncation error, i.e., the distance between the exact and the numerical solutions. Mainly, such methods work in two stages at each integration step, starting from an enclosure [y j ] y(t j ; y 0 ) at time t j of the exact solution, we proceed by: i. a computation of an a priori enclosure [ y j+1 ] of the solution y(t; y 0 ) for all t in the time interval [t j , t j+1 ]. This stage allows one to prove the existence and the uniqueness of the solution.

ii. a computation of a tightening of state variable [y j+1 ] y(t j+1 ; y 0 ) at time t j+1 using [ y j+1 ] to bound the term lte Φ (F, G, y j , t j , h).

Sometimes, additive constraints can be considered, coming from a mechanical constraint, an energy conservation, a control law, or whatever [START_REF] Alexandre | Validated Simulation of Differential Algebraic Equations with Runge-Kutta Methods[END_REF]. It is important to understand that these constraints are linked to the system, valid all the time, and thus semantically different from a constraint coming from a satisfaction problem which are valid only at particular time instants. Finally a constrained parametric IIVP is considered where

         ẏ = F (t, y, x, p, u) 0 = G(t, y, x, p, u) 0 = H(t, y, p, u) y(0) ∈ Y 0 , x(0) ∈ X 0 , p ∈ P, u ∈ U, t ∈ [0, t end ] , (2)
H : R × R n × R r × R s → R c .
Methods defined in [START_REF] Alexandre | Validated Simulation of Differential Algebraic Equations with Runge-Kutta Methods[END_REF][START_REF] Alexandre | Validated Explicit and Implicit Runge-Kutta Methods[END_REF] can take into account such kind of constraints by adapting the first stage of validated numerical integration methods, that is the computation of the a priori enclosure.

A validated simulation starts with the interval enclosures [y(0)], [x(0)], [p] and [u] of respectively, Y 0 , X 0 , P, and U. It produces two lists of boxes:

• the list of discretization time steps. {[y 0 ], . . . , [y end ]};

• the list of a priori enclosures:

{[ y 0 ], . . . , [ y end ]}.
Based on these lists, two functions depending on time can be defined

R : R → IR n t → [y] (3)
with {y(t; y 0 ) :

∀y 0 ⊆ [y 0 ]} ⊆ [y],
and

R : IR → IR n t, t → [ y] (4)
with {y(t; y 0 ) :

∀y 0 ∈ [y 0 ] ∧ ∀t ∈ [t, t]} ⊆ [ y].
Function R, defined in (3), is obtained by new applications of validated integration method starting from [y k ] at t k , and finishing at t with t k < t < t k+1 . Function R, defined in (4), is obtained with the union of [ y k ] with k = a, . . . , b and t a < t < t < t b . These functions are then strictly conservative.

More abstractly, the functions R and R define two interval enclosures of the solution function of differential equations defined in Equation (2).

Example 3 We consider an Initial Value Problem based on the following constrained differential algebraic equation

                       F :      ẏ1 = -y 3 y 2 -(1 + y 3 )x 1 ẏ2 = y 3 y 1 -(1 + y 3 )x 2 ẏ3 = 1 G : 0 = (y 1 -x 2 )/5 -cos(y 2 3 /2) 0 = (y 2 -x 1 )/5 -sin(y 2 3 /2) H : 0 = x 2 1 + x 2 2 -1 (5)
and the initial values y(0) = (5, 1, 0) and x(0) = (-1, 0). Interesting variables are y 1 and y 2 , the states of the system. A picture of the result of validated simulation is given in Figure 1 and The problem given by Equation (5) has a known exact solution which is

y(t) =      y 1 = sin(t) + 5 cos(t 2 /2) y 2 = cos(t) + 5 sin(t 2 /2) y 3 = t .
If we compute the previous enclosures with the exact solution (by combining interval arithmetic with bisection method to obtain sharp results), we find: The enclosure property of R(t) and R(t) is then verified.

CSP with Differentials Constraints

Constraint Satisfaction Differential Problems (CSDP) extending the definition of NCSP have been defined in [START_REF] Cruz | Constraint reasoning over differential equations[END_REF] by adding new variables to represent time derivative, i.e., dy dt and a new kind of constraints of the form dy dt = f (y, t). The main drawback of CSDPs is that temporal properties cannot be encoded because the time variable t is not a variable of the problem. So, the property "for all t in [0, 10], the system does not collide with an obstacle" cannot be represented as a CSDP in an obvious manner.

The work [START_REF] Goldsztejn | Including ordinary differential equations based constraints in the standard CP framework[END_REF] considers differential equations as an implicit definition of temporal functions. In consequences, considering only solution function Φ of differential equations makes an easy embedding of differential equations into NCSP framework. Unfortunately, NCSPs naturally express existential problems, i.e., problems of the form ∃x ∈

[x], n i=1 f i (x) = 0 ∧ m i=1 g i (x)
0. While this approach is appealing to simply embed differential equations into constraint programming framework, it is also unable to express constraints where the time variable is a variable of the problem.

In both previous work [START_REF] Cruz | Constraint reasoning over differential equations[END_REF] and [START_REF] Goldsztejn | Including ordinary differential equations based constraints in the standard CP framework[END_REF], the main drawback is the lack of quantification on variables. A proposition of extension of these approaches to cope with this limitation is the main contribution of the paper.

In order to be able to define a CSDP which can answer a question of the type "Will the system collide with an obstacle before ten seconds?", we propose a new definition of CSDP which mainly introduces quantification over variables.

Definition 2 (QCSDP) Let S be a differential system as defined in (2) and t end ∈ R + the time limit of the problem handled. A QCSDP is a CSP defined by:

• a set of variables V including at least the time variable t, a vector y 0 of dimension n representing the initial state-space of S, variables associated to the parameters p, variables representing the control input u of S. We represent these variables by the vector v;

• an initial domain D containing at least [0, t end ] for time t, the initial domain Y 0 of initial condition y 0 , U for control input u, and P for parameters p;

• a set of constraints C = {c 1 , . . . , c } composed of predicates over sets, that is, constraints of the form:

c i (v) ≡ Qv ∈ D i .f i (v) A, ∀1 i e,
with Q ∈ {∃, ∀}, f i : ℘(R |V| ) → ℘(R q ) stands for non-linear arithmetic expressions defined over variables v and solution of differential system S, y(t; y 0 , p, u) ≡ y(v). The operator denotes {⊂ , ⊃ , ∩ =∅ , ∩ =∅ }, i.e., the set-based operations, and the set A ⊂ R q where q > 0.

Note that ∩ =∅ stands for the predicate f (v) ∩ A = ∅ and ∩ =∅ stands for the predicate f (v) ∩ A = ∅. ℘(E) denotes the power-set of the set E, while |E| denotes the cardinality of E.

Note that in Definition 2, constraints are expressed with set-based operator instead of equalities and inequalities to simplify the notations despite of the complete equivalence with equalities and inequalities. Moreover, set-based operations will emphasize, in Section 4, one important difficulty due to the representation of sets with respect to the correctness of the proposed approach.

As the QCSDP over sets is not computable in general, a set abstraction has to be used. Moreover, solving arbitrary quantified constraints over reals is a hard problem hence we focus on certain classes of problems related to the robotic domain. In particular, three main classes are considered. Problem 1 deals with the design of robot platform or control algorithms and for example is treated with a CSP approach in [START_REF] Sandretto | Appropriate design guided by simulation: An hovercraft application[END_REF] or in [START_REF] Sandretto | Tuning PI controller in non-linear uncertain closed-loop systems with interval analysis[END_REF]. Problem 2 deals with the design of robust control algorithms. It is the case study of Section 5 and it has been considered in the context of non-linear sampled switch systems in [START_REF] Le Coent | Control of Nonlinear Switched Systems Based on Validated Simulation[END_REF]. Problem 3 consists to find the subset of initial states such that there exists at least one evolution "viable", in the sense that at each time, the state of the evolution remains confined to the viability kernel [START_REF] Aubin | A survey of viability theory[END_REF]. The problem of viability kernel computation has been considered by the help of our method in [START_REF] Monnet | Viability kernel computation based on interval methods[END_REF].

Problem 1 (Appropriate design) A parametric dynamics of a robot described by a differential system is considered, the parameters represent some possible choices in the characteristics of the robot such as its dimension, the power of the engine and so on. An appropriate design problem aims at finding a value of the parameters such as a set of functional requirements is fulfill. In other words, QCSDP problems are made of constraints over variables p, u and y 0 based on quantification ∃p, ∀u, ∀y 0 , and some temporal constraints (see Section 4).

Problem 2 (Control synthesis) A parametric dynamics of a robot described by a differential system is considered, the parameters represent perturbation. A control synthesis problem aims at finding value of the control u such that some requirements are fulfill taking into account all possible perturbations. In other words, QCSDP problems are made of constraints over variables p, u and y 0 based on quantification ∃u, ∀p, ∀y 0 , and some temporal constraints (see Section 4).

Problem 3 (Viability kernel) A parametric dynamics of a robot described by a differential system is considered, the parameters represent perturbation. A viability kernel problem aims at finding initial values of the system such that there exists at least one evolution remaining into the viability kernel taking into account all possible perturbations. In other words, QCSDP problems are made of constraints over variables p, u and y 0 based on quantification ∃y 0 , ∀p, ∃u, and some temporal constraints (see Section 4).

Note that it is also possible to add a cost function in the Definition 2, and hence consider constraint optimization problems.

Reasoning Under Set Abstraction

In Definition 2, QCSDPs are defined over sets of values and hence they are not computable. Abstraction with boxes is considered in this section. By denoting α() the box abstraction, the constraints of QCSDP presented above becomes:

α(c i (v)) ≡ α(Qv ∈ D i .f i (v) A) (6)
≡ Qv ∈ [d i ].α(f i (v) A) (7)
≡ Qv ∈ [d i ].[f i ](v)α( A) . (8)
From Equation [START_REF] Janssen | Multistep filtering operators for ordinary differential equations[END_REF] to Equation [START_REF] Cruz | Constraint reasoning over differential equations[END_REF], abstraction consists on an over-approximation of the domain of variables. Then, from Equation [START_REF] Cruz | Constraint reasoning over differential equations[END_REF] to Equation ( 8), the interval enclosure is used as an abstraction of function f . It means that only enclosure of the solution of differential equations is considered, even if an inner approximation could be available. Indeed, enclosure is mainly favored because the approach presented in this paper is oriented to safety problems and not to liveness. Abstraction α( A) is not distributive in a safe manner. Applied to the set A, the function α ∈ {Hull, Int} stands for an interval abstraction function which maps a set A to a set of boxes. Hull(A) stands for the interval over-approximation of the set A while Int(A) stands for an interval representation (a single interval or a paving) of an under-approximation of the set A. By an abuse of notation, in Definition 3 will denote the abstraction of the set-theoretic operator given in Definition 2 but it is important to keep in mind that the implementation is strongly dependent on the abstraction of the set A.

The possible operations are gathered in the Table 1. Due to the abstraction, all these operators are weak. Nevertheless, some results after abstraction imply the original logic formula. In Table 1, question mark "?" means that neither true nor false implies a guaranteed result on the original logic formula. On the contrary, true means that if the abstract formula is valid, then the original one is valid, i.e., [f ](v) ⊂ Int(A) ⇒ f (v) ⊂ A. In the same way, false means that if the abstract formula is not valid, then the original one is not valid, i.e., ¬(

[f ](v) ∩ =∅ Hull(A)) ⇒ ¬(f (v) ∩ =∅ A). Abstraction function α( A)
can be distribute into α( )α(A) by following carefully the results presented in Table 1.

Finally, we propose the Box-QCSDP formalism, which is a QCSP where sets are represented by boxes and is formally defined in Definition 3. It mainly introduces two abstractions: the abstraction of the solution of differential systems as defined in (2); and the abstraction of constraints c i , in particular, the abstraction of the set A.

The first and major abstraction allows us to manage with the differential system in a CSDP. It is performed by a validated simulation method as described in 2.3. This approach provides to us two abstraction functions guaranteed to enclose all the behaviors of the differential system.

The second abstraction is done to be able to handle the sets by an outer or an inner approximation, depending on the associated logic operation. In order to correctly translate a satisfaction problem, these capabilities offered by QCSDP have to be used in the right way.

Definition 3 (Box-QCSDP) Let S be a differential system as defined in (2) and t end ∈ R + the time limit of the problem handled. A Box-QCSDP is defined by:

• a set of variables V including at least the time variable t, a vector y 0 of dimension n representing the initial state-space of S, variables associated to the parameters p, variables representing the control input u of S. We represent these variables by the vector v;

• an initial box domain [d] containing at least [0, t end ] for time t, the initial box [y 0 ] of initial condition y 0 , a box [u] for control input u, and a box [p] for parameters p;

• a set of constraints C = {c 1 , . . . , c e } composed of predicates over sets, that is, constraints of the form

c i ≡ Qv ∈ [d i ].[f i ](v)α ( A) , ∀1 i e,
with Q ∈ {∃, ∀}, [f i ] : IR |V| → IR q stands for an inclusion function of non-linear arithmetic expressions defined over variables v and an interval enclosure of the solution of differential system S, [y](t; [y 0 ], [p], [u]). The set theory based operations are expressed by ∈ {⊂, ⊃, ∩ =∅ , ∩ =∅ } and a set A ⊂ R q where q > 0. The function α( A) maps a part of logic operations into abstract ones, following Table 1.

Definition 3 is very dependent of the implementation. The following remarks give some insight to produce safe results.

∀t ∈ [0, t end ], R(t) ⊆ Int(A) In A at τ R(τ ) ⊆ Int(A) Has crossed A ∃t ∈ [0, t end ], R(t) ∩ Hull(A) = ∅ Go out A ∃t ∈ [0, t end ], R(t) ∩ Hull(A) = ∅ Has reached A R(t end ) ∩ Hull(A) = ∅ Finish in A R(t end ) ⊂ Int(A)
Remark 2 Solutions of differential equations need to be pre-computed to obtain our abstraction functions R(t) and R([t]). If the parameters or control inputs involved in the definition of the differential equation is reduced by the CSDP solver, then the abstraction functions need to be computed again. It is a part of the propagation algorithm of a CSP.

QCSDP, see Definition 2, allows one to describe a complex satisfaction problem in order to reason about a system defined by a differential equation. Nevertheless, a definition of the temporal constraints has not been done yet, it is the purpose of the end of this section. We consider the main interesting properties for robotic applications, by providing verbal operators in the sense of "system has the property of ", such as "stay in till".

By convention, these Boolean operations provide a guaranteed true answer. It means that "true" is proven while "false" means "uncertain" due to the abstraction of boxes. In table 2, the different available verbal operators and their Box-QCSDP translation are gathered. Despite this list of specifications is limited it emphasizes the main difficulty to reason on abstraction of the behaviors of robots.

From table 2 and due to the abstraction, "Go out" is not equivalent to ¬("Stay in"). Operators "Has reached" and "Has crossed" are poor tests (true also means uncertain). They are equivalent but the validity of "Has reached" at the end of a simulation implies a notion of stability. On the other hand, the negation of these latter operators are strong and are useful to speed up a Box-QCSDP solver.

Following the same generalization than for Constraint Optimization Problems based on CSP, a cost function can be added to Box-QCSDP to perform optimization. Definition 4 (Optimization on Box-QCSDP) A Box-Quantified Constraint Optimization Differential Problem involves in addition of a Box-QCSDP an objective function Q(v 1 , . . . , v n ) that has to be maximized or minimized over the set of all feasible solutions.

The details on algorithms to solve optimization problems is out of the scope of this paper, but there are the same than the ones used for CSP [START_REF] Hansen | Global Optimization Using Interval Analysis[END_REF].

Box-QCSDP and this list of specifications have been implemented in a library named DynIBEX In Box-RRT algorithm, Step 3 is a Box-QSCDP. In particular, there are several constraints that must be taken into account while the path is generated. First, the path generation works on a finite dimension map and so trajectories must stay inside this map. Obviously, a constraint expressing the avoidance of obstacles shall be considered. An other constraint has to be defined to detect when a trajectory reaches the targeted region of the state-space. In consequence, the path generation is translated into Box-QCSDP which is defined in Equation [START_REF] Liu | Synthesis of reactive switching protocols from temporal logic specifications[END_REF]. Note that this Box-QCSDP has to be solved at each iteration of the Box-RRT algorithm until a path has been found.

In Equation [START_REF] Liu | Synthesis of reactive switching protocols from temporal logic specifications[END_REF], the constraint c 1 states that the trajectory [car](t) from [s nearest ] to [s new ] has to stay inside the map Map. The constraint c 2 states that the trajectory [car](t) shall not intersect an obstacle and finally the constraint c 3 states that a path is found if one state of trajectory is included into [s f ].

An example of the exploration of the state-space produced by the Box-RRT algorithm is given in Figure 3. The red squares stand for obstacles while the blue square represents [s f ]. The small black square in the center of the figure is [s 0 ] while the black square inside [s f ] is the final state of the trajectory reaching the target.

Conclusion

In this paper, we presented a framework for the formal verification of robotic behaviors in presence of uncertainties. Based on the constraint satisfaction problem formalism and by the help of interval tools such as interval arithmetic and the validated integration of differential equations, a method to verify the safety, in a guaranteed manner, has been developed. Our approach mainly exploits box abstraction of uncertainties, sets, functions and differential equations. Box abstraction is a very powerful approach to manage with uncertainties and, more particularly, with differential equations. Nevertheless, this abstraction introduces another difficulty, which is the handling of logic formula, because some logic operators become weak. We carefully analyzed available operators and a new formalism called Box-QCSDP has been finally presented. A well known problem, the robot path planning, has been studied under the formalism of Box-QCSDP. An algorithm based on Rapid-exploring Random Tree has been defined including constraints (obstacles and objectives) and solved with our method. The resulting path is guaranteed by construction, considering the robot dynamics, which is an interesting improvement of classical path planners. To conclude, we presented in this paper a new formalism that allows one to solve complex problems linking dynamical systems and temporal logic, in a validated way.

[- 2

 2 , 5] + [-8, 12] = [-10, 17] [-10, 17] -[-8, 12] = [-10, 17] + [-12, 8] = [-22, 25] [-10, 17] -[-2, 5] = [-15,

  [x], [y] and [z] from a forward evaluation [z] = [z] ∩ ([x] + [y]), and from two backward evaluations [x] = [x] ∩ ([z] -[y]) and [y] = [y] ∩ ([z] -[x]).

Figure 1 :

 1 Figure 1: Results of validated simulation: y 2 w.r.t. y 1 .

Figure 2 .

 2 Figure Results of validated simulation: y 2 w.r.t. time.

•

  y(1.1) = ([5.003, 5.003]; [3.297, 3.297]; [1.1]) • y([0.7, 0.8]) = ([[5.464, 5.494][1.977, 2.269]; [0.7, 0.8])

Remark 1

 1 If a validated simulation cannot reach t end and fails at t = t fail then R(t) = [-∞, ∞], for all t t fail . Hence the enclosure of [y](t; [y 0 ], [p], [u]) is always correct and respects the contractor definition.

V

  = {t, x 0 , y 0 , θ 0 , v, δ} ,D = {[t] = [0, 1], [x 0 ] = [0.0, 0.3], [y 0 ] = [0.0, 0.3], [θ 0 ] = [0.0]} , C = { c 1 = ∀s 0 ∈ [s 0 ], ∀t ∈ [t] • [car](t; s 0 ) ⊂ Map, c 2 = ∀s 0 ∈ [s 0 ], ∀t ∈ [t], ∀s a ∈ [L a ] • [car](t; s 0 ) ∩ =∅ s a , c 3 = ∀s 0 ∈ [s 0 ], ∃t ∈ [t] • [car](t; s 0 ) ⊂ [s f ] } .

Figure 3 :

 3 Figure 3: Results of Box-RRT

  is a set of intervals such that [v i ] contains all possible values of v i . It can be represented by a box [v] gathering all [v i ];

Table 1

 1

		: Available operations
			α(A)
			Int(A) Hull(A)
		⊂	true	?
	α(f ) [f ]	⊃ ∩ =∅	false ?	? true
		∩ =∅	?	false

  1 [START_REF] Alexandre | DynBEX: a Differential Constraint Library for Studying Dynamical Systems[END_REF]. map, denoted by Map, represented by a paving, i.e., a set of non-overlapping boxes and from a list of obstacles, denoted by L a , Box-RRT looks for a safe path between [s 0 ] and [s f ], avoiding all obstacles in L a and staying in Map. To do this, a tree T of possible paths is generated from a random walk inside the state-space represented by a paving. Box-RRT is an iterative algorithm which mainly follows three steps until[s f ] is reached i. pick a random state [s random ], ii. find the closest state of [s random ] in the tree T named [s nearest ], iii. pick a control u and compute [s new ] which is the final state of the trajectory starting from [s nearest ] with the control u. If this trajectory reaches an obstacle, an other state is chosen. Otherwise if [s new ] ⊂ [s f ] then a path is found and the algorithm terminates. Otherwise, the target [s f ] has not been reached and [s new ] is a possible safe state being on the path to reach the target [s f ]. Hence, [s new ] is added as a new node in T and the exploration of the state-space continues.

http://perso.ensta-paristech.fr/ ~chapoutot/dynibex/

Acknowledgment

The authors are grateful to François Pessaux for his pertinent comments on the early versions of this work.

This research benefited from the support of the "Chair Complex Systems Engineering -Ecole Polytechnique, THALES, DGA, FX, DASSAULT AVIATION, DCNS Research, ENSTA ParisTech, Télécom ParisTech, Fondation ParisTech and FDO ENSTA". This research is also partially funded by DGA MRIS.

Case Study: Robot Path Planning

We present an application of Box-QCSDP based on the algorithm Box-RRT defined in [START_REF] Pepy | Reliable robust path planning with application to mobile robots[END_REF] which aims at producing a robust path for a car-like vehicle while avoiding obstacles.

The dynamic model of the car is defined by

The state of the car is defined by the car position (x, y) and its heading θ. L is the distance between the front and rear wheels. There are two control inputs v for the longitudinal speed and δ ∈ [δ min , δ max ] for the steering angle. The noise element ε = [-0.001, 0.001] stands for the steering angle imprecision. A second source of uncertainties is the initial position of the vehicle. We denote by car(t; s 0 ) the solution of Equation ( 9) and by [car](t; s 0 ) its interval enclosure produce by guaranteed integration methods. We denote by s 0 = (x 0 , y 0 , θ 0 ) and [s 0 ] = [x 0 ] × [y 0 ] × [θ 0 ] the state variables and the initial conditions. We denote by [s f ] the targeted state-space. Box-RRT algorithm extends Rapid Random Tree (RRT) algorithm with interval analysis tools. A brief recall of the Box-RRT algorithm is given before introducing our formalization with Box-QCSDP. From a


Thomas Caissard

David Coeurjolly

Jacques-Olivier Lachaud

Tristan Roussillon

Heat kernel Laplace-Beltrami operator on digital surfaces

 to digital surfaces embedded in 3D. The core of the method relies on an accurate estimation of measures associated to digital surface elements. We experimentally evaluate the interest of this operator for digital geometry processing tasks.

Introduction

Objectives In geometry processing, Partial Differential Equations (PDEs) containing Laplace-Beltrami operator arise in many applications such as surface fairing, mesh smoothing, mesh parametrization, remeshing, mesh compression, feature extraction or shape matching (see [START_REF] Lévy | Spectral Mesh Processing[END_REF] for an extensive survey). On digital surfaces, few digital Laplace-Beltrami operators has been proposed and none has been evaluated in terms of multigrid convergence (convergence of the operator toward the continuous one in digitization of smooth manifolds on grid with decreasing gridstep).

Contributions In this article, we propose a discrete Laplace-Beltrami operator on digital surfaces (boundaries of subsets of Z 2 embedded in 3D). This new operator adapts Belkin et al. [START_REF] Belkin | Discrete laplace operator on meshed surfaces[END_REF] on our specific data. The method uses an accurate estimation of areas associated with digital surface elements. This estimation is achieved through a convergent digital normal estimator described in [START_REF] Coeurjolly | Multigrid convergent principal curvature estimators in digital geometry[END_REF]. We show experimental convergence of our operator but also that none of the existing approaches adapted to digital surfaces achieves such convergence. Finally, we illustrate the interest of the discretized laplacian on digital surface geometry processing.

Related works

The Laplace-Beltrami operator being a second order differential operator (divergence of the function gradient, see Sect. 2), a discrete calculus framework is required to define such operator on embedded combinatorial structures such as meshes or digital surfaces. First works on discrete calculus may be found in the Regge Calculus [START_REF] Regge | General relativity without coordinates[END_REF] for quantum physics, where tetrahedra in combination with edge lengths are used. Works on geometric acquisition devices and models drove studies toward calculus working on meshes and more generally on simplicial complexes. Early works include a definition of the Laplace-Beltrami operator using the classical cotangent formula [START_REF] Pinkall | Computing discrete minimal surfaces and their conjugates[END_REF] for solving the problem of minimal surfaces, which is an analog of the standard finite element method [START_REF] Lévy | Spectral Mesh Processing[END_REF]. Exact calculus generalizing the cotangent discretization in 2D based on finite elements [START_REF] Polthier | Identifying vector field singularities using a discrete Hodge decomposition[END_REF] emerged from the German school but with a restriction to triangular complexes.

In a more generic discrete calculus perspective, the Discrete Exterior Calculus (DEC) framework was then developed in the computational mathematics and geometry processing community. Another more recent formulation of the DEC comes from Hirani's thesis [START_REF] Hirani | Discrete exterior calculus[END_REF] and later by the monograph [START_REF] Desbrun | Discrete exterior calculus[END_REF]. On triangular meshes, DEC based Laplace-Beltrami operator and the cotangent based one coincide.

In [START_REF] Hildebrandt | On the convergence of metric and geometric properties of polyhedral surfaces[END_REF], authors show that under some strong assumptions, the cotangent laplacian on a triangular mesh converges to the continuous one when the mesh interpolates a smooth manifold with increasing precision (with a continuous oneto-one map between the mesh and the manifold which would not be the case on digital surfaces). The operator converges in the sense of distributions and the authors show that pointwise convergence in the l 2 sense does not generally hold. On triangular meshes, Belkin et al. [START_REF] Belkin | Discrete laplace operator on meshed surfaces[END_REF] have proposed a first Laplace-Beltrami operator that converges in the l ∞ uniform case. The digital laplacian operator we propose is an extension to digital surfaces of such operator.

In digital geometry, many estimators of differential quantity have been proposed and there exists multigrid convergent estimators for many quantities such as length, tangent and curvature in 2D (see [START_REF] Coeurjolly | Multigrid Convergence of Discrete Geometric Estimators[END_REF] for a complete survey) surface area [START_REF] Lachaud | Properties of gauss digitized shapes and digital surface integration[END_REF], normal vectors and curvature tensor in dimension 3 [START_REF] Coeurjolly | Multigrid convergent principal curvature estimators in digital geometry[END_REF]. A preliminary approach can be found in [START_REF] Mercat | Discrete complex structure on surfel surfaces[END_REF][START_REF] Cartade | Mesh parameterization with generalized discrete conformal maps[END_REF]. It focuses on the conformal map computation of a digital surface, which is a related problem involving the definition of a Laplace-Beltrami operator. However, their definition is based on the cotangent formula, which lacks pointwise convergence. When designing a discrete version of the laplacian operator, not all properties of the continuous one can be expected at the same time. In [START_REF] Wardetzky | Discrete Laplace operators: No free lunch[END_REF], entitled "Discrete Laplace operators: No free lunch.", the authors have proposed a formal evaluation of such properties. We position our new digital Laplace-Beltrami operator with respect to this analysis.

Outline After introducing mathematical definitions, we review the classical approaches to define a discrete laplacian operator and compare their properties in Section 2. We then formalize our operator in Section 3. In Section 4, we experimentally evaluate our proposal in terms of multigrid convergence and geometry processing applications.

Discretizations of the Laplace-Beltrami operator and their properties

We first describe various discretizations of the Laplace-Beltrami operator on triangular meshes. We then check several desired properties of laplacian operators using [START_REF] Wardetzky | Discrete Laplace operators: No free lunch[END_REF] as a baseline for comparisons.

Preliminaries and classical discretizations on triangular meshes

Let M be a 2-Riemannian manifold with or without boundary embedded in R 3 . More precisely, we consider a pair (M, g) where M is a smooth manifold and g is a Riemannian metric on M (i.e. with a known intrinsic notion of distances).

The intrinsic smooth Laplace-Beltrami operator [START_REF] Rosenberg | The Laplacian on a Riemannian Manifold[END_REF] is defined as:

∆ : C 2 (M ) → C 2 (M ) u → div(∇ u), (1)
be where C 2 is the set of twice differentiable functions with the second derivative continuous of M (in the literature, alternative definitions may consider "-div(∇ u)" for ∆ u).

Let Γ be a combinatorial structure (a triangular mesh for instance), V (Γ ) its set of vertices and F (Γ ) its faces. Let u : M → R be a twice differentiable function. We suppose that V (Γ ) is a sampling of M (i.e. V (Γ ) ∈ M ). In other words, u(w) is perfectly defined for w ∈ V (Γ ).

A first simple discretization only considers the combinatorial structure of Γ . Such laplacian is either called graph laplacian or combinatorial laplacian of Γ [START_REF] Zhang | Discrete combinatorial laplacian operators for digital geometry processing[END_REF]:

(L COM BI u)(w) := p∈link0(w) u(p) -deg(w)u(w) , (2)
for all w ∈ V (Γ ) where link 0 (w) is the set of points V (Γ ) adjacent to w and deg(w) is the degree of w in Γ .

A more complex approach can be defined using DEC operators [START_REF] Hirani | Discrete exterior calculus[END_REF][START_REF] Desbrun | Discrete exterior calculus[END_REF]. Using an arbitrarily embedded dual structure of Γ , the laplace operator can be shown to be a classical weighted double finite difference:

(L DEC u)(w) := 1 | w| p∈link0(w) | e wp | |e wp | (u(p) -u(w)), (3)
where is the Hodge-duality star operator acting on discrete forms (see [START_REF] Hirani | Discrete exterior calculus[END_REF]), and | • | the measure of a k-cell. As illustrated in Fig. 1, | e wp | would be the length of the segment orthogonal to e wp . If we set all measures to one, L DEC coincides with L COM BI . By fixing the dual of Γ to be the voronoi diagram of its vertices and by computing the measures as Euclidean lengths and areas of such dual complex, the DEC operator coincides exactly with the famous cotan laplacian [START_REF] Pinkall | Computing discrete minimal surfaces and their conjugates[END_REF]:

(L COT u)(w) := 1 2A w p link0(w) (cot(α wp ) + cot(β wp )) (u(p) -u(w)) , (4)
where A w is one third of the area of all incident triangles to vertex w, α wp and β wp are the angles opposing the corresponding edge e wp (see Fig. 1).

Fig. 1. Illustration of LDEC (left), and LCOT (right) on triangular meshes. For LCOT the area of integration Aw is one third the area of all triangles incident on vertex w in green. For LDEC the dual structure is in orange and the dual of the edge ewp is in blue.

Finally, we detail the definition of the mesh laplacian from [START_REF] Belkin | Discrete laplace operator on meshed surfaces[END_REF]. Let g : M × (0, T ) → R be a time-dependent function which solves the partial differential equation called the heat equation:

∆g(x, y) = ∂ ∂t g(x, t), (5)
with initial condition g 0 = g(•, 0) : M → R which is the initial temperature distribution. An exact solution called the heat kernel [START_REF] Rosenberg | The Laplacian on a Riemannian Manifold[END_REF] is:

g(x, t) = y∈M h(t, x, y)g 0 (y)dV, (6)
where h ∈ C ∞ (R + × M × M ) is called the heat kernel. Using one of the various results on the heat kernel approximation [START_REF] Varadhan | On the behavior of the fundamental solution of the heat equation with variable coefficients[END_REF][START_REF] Molchanov | Diffusion processes and riemannian geometry[END_REF] we have

g(x, t) = 1 4πt y∈M e -||x-y|| 2 4t g 0 (y)dV, (7)
which corresponds to the heat kernel in R 3 . Injecting Eq.( 7) into Eq.( 5), applying a finite time difference and knowing that the integral of e over M is one:

∆g(x, t) = lim t→0 1 t y∈M h(t, x, y)(g 0 (y) -g 0 (x))dV. ( 8
)
The previous equation can be seen as a convolution between differences of g and a time dependent Gaussian e (see Fig. 2). Then, the mesh laplace operator [START_REF] Belkin | Discrete laplace operator on meshed surfaces[END_REF] Fig. 2. Illustration of the convolution between the heat kernel e in green and differences of g evaluated at x and y. The kernel distance is in orange.

on Γ is:

(L M ESH u)(w) := 1 4πt 2 f ∈F (Γ ) A f 3 p∈V (f ) e -||p-w|| 2 4t (u(p) -u(w)), (9)
where A f is the area associated to the face f .

Desired properties of a discrete laplacian

As discussed in [START_REF] Wardetzky | Discrete Laplace operators: No free lunch[END_REF] all properties of the continuous Laplace-Beltrami operator may not be preserved when discretizing it. We consider the discrete laplacian as a linear operator acting on values u := {u p } on V (Γ ) (represented as a vector in R |V (Γ )| ). Such operator can thus be denoted as a matrix L with components l ij . Hence, v := Lu would be the resulting laplacian of u. Expected properties of the discrete Laplace-Beltrami operators are:

Symmetry (SYM) ij = ji for 0 ≤ i < |V (Γ )|
. This is very useful when it comes to solve linear systems as solvers are usually more performant when the matrix is symmetric. Locality (LOC) ij = 0 if and only if i and j shares a common edge. The locality property gives very sparse matrices decreasing drastically memory consumption. It also opens a panel of very fast linear system solvers. Linear Precision (LIN) Lu = 0 whenever u is a linear function restricted to a plane. Positive Weights (POS) ij ≥ 0 for i = j. Furthermore, for each vertex i, there exists a vertex j such that ij > 0. Positive Semi-Definiteness (PSD) The matrix is symmetric positive semidifinite regarding the standard inner product and has a one-dimensionnal kernel. (SYM) and (POS) imply (PSD), but (PSD) does not imply (POS). This property ensures that the basis generated by the eigenvectors of L is orthogonal and that the eigenvalues are real. Dirichlet Convergence (CON) When considering a sequence of meshes {Γ i } converging to M in a given sense as i → ∞, we want the laplacian sequence L i to converge to ∆ with respect to the discrete Dirichlet problem. The convergence is mendatory when we seek approximate solutions of partial differential equations.

The (CON) property requires a formal definition of the sequence {M i }. In addition to [START_REF] Wardetzky | Discrete Laplace operators: No free lunch[END_REF], we add the following property:

Pointwise Convergence (PCON) For a given sequence meshes {Γ i }, we want the associated laplacians {L i } to converge to ∆ in a pointwise l 2 or l ∞ sense. This notion of convergence is stronger than (CON) which is implied by (PCON).

For the laplacian operators on digital surfaces (e.g. L h we define in Section 3), we consider a sequence of combinatorial structures defined as the boundaries of the Gauss digitizations of M (see next section). We thus consider the (DP-CON) property as the pointwise convergence of the operator for multigrid digital surfaces. Table 1 compares classical laplacian discretizations with respect to these properties. The (DPCON) property has been evaluated experimentally in Section 4. L COM BI being purely combinatorial, the same operator can be used for both meshes and digital surfaces. For L COT and L M ESH , we have considered the Marching-cubes representation of the digital surfaces (more precisely, a continuous analog of the dual surface [START_REF] Lachaud | Continuous analogs of digital boundaries: A topological approach to iso-surfaces[END_REF]). Properties of L h are discussed in the next section.

New Laplace-Beltrami operator on digital surfaces

In the discretization schemes discussed above, the points of the combinatorial structure interpolate the underlying manifold. In order to define a discrete Laplace-Beltrami operator on digital surfaces, a challenge is to work with a combinatorial structure that only approximates the underlying manifold in a Hausdorff sense. First, we extend Eq.( 9) to digital surfaces. Proper definitions of a digital surface can be found in [START_REF] Klette | Digital geometry : geometric methods for digital picture analysis[END_REF][START_REF] Lachaud | Properties of gauss digitized shapes and digital surface integration[END_REF]. Let us recall the Gauss digitization process:

Definition 1 (Gauss digitization). Let h > 0 be the sampling grid step. The Gauss Digitization of an Euclidean shape X ⊂ R n is defined has D h (X) := X ∩ (hZ) d where d is the dimension.

For smooth object X in dimension 3 with boundary M := ∂X, the digital surface is defined as the topological boundary of D h (X), denoted ∂ h X (see [START_REF] Lachaud | Properties of gauss digitized shapes and digital surface integration[END_REF] for details). More precisely, the digital surface has a cellular representation in a cartesian cubical grid and is composed of points of dimension 0 (pointels, E 0 ), straight segments of dimension 1 (linels, E 1 ) and squares of dimension 2 (surfels, E 2 ). The Gauss digitization ∂ h X is an O(h)-Hausdorff approximation of M [START_REF] Lachaud | Properties of gauss digitized shapes and digital surface integration[END_REF]. As a consequence, we need to map the smooth function u defined on M to ∂ h X: Definition 2 (Extension of u to ∂ h X). Given a smooth function u on M , we define the extension ũ of u to ∂ h X as

ũ(s) := u(ξ(ṡ)) ,
where ṡ is the centroid of the surfel s ∈ E 2 , and ξ is the map that projects a point of ∂ h X onto the closest point of M .

We show below how to adapt the definition of [START_REF] Belkin | Discrete laplace operator on meshed surfaces[END_REF], recalled in Eq.( 9), to digital surfaces. We chose this approach because L M ESH has an interesting pointwise convergence for triangular meshes. Our Laplace-Beltrami operator is thus defined as follows:

Definition 3 (Digital Laplace-Beltrami operator). The digital Laplace-Beltrami operator is defined on ∂ h X as:

(L h ũ)(s) := 1 4πt 2 h r∈E 2 e - || ṙ-ṡ|| 2 4t h µ(r)(ũ(r) -ũ(s)), ( 10
)
where the sum is taken over all surfels of ∂ h X, ṙ is the centroid of the surfel r, µ(s) is equal to the dot product between an estimated normal and the trivial normal orthogonal to the surfel s and t h is a function of h tending to zero as h tends to zero.

The quantity µ(s) is called the measure of the surfel s: it is the area of the projected surfel s onto the tangent plane induced by the estimated normal.

Normal vectors are estimated using the estimator presented in [START_REF] Coeurjolly | Multigrid convergent principal curvature estimators in digital geometry[END_REF][START_REF] Levallois | Parameter-free and multigrid convergent digital curvature estimators[END_REF] which has the multigrid convergence property. Note that summing µ for each surfel of the surface leads to an estimation of the global area of the shape boundary, which itself has a multigrid convergence property [START_REF] Lachaud | Properties of gauss digitized shapes and digital surface integration[END_REF]. This surfel measure is a key ingredient of the digital formalization of the operator leading to an experimental multigrid convergence and isotropic properties when used for diffusion (see Sect.

4).

If we index surfels in E 2 , L h has a |E 2 |×|E 2 | matrix representation L h defined as follows:

(L h ) ij =          1 4πt 2 h e - || ṡj -ṡi || 2 4t h µ(s j ) if i = j - k =i (L h ) ik if i = j (11)
In other words, if ũ is constant and equal to 1, for a surfel s ∈ E 2 of index i, we have (L h ũ) i = (L h ũ)(ṡ i ).

For the properties listed in Table .1, we can first observe that Eq.( 11) implies that we have property (POS) but not (SYM). As L h performs a convolution on the complete surface with a Gaussian kernel, we do not have the (LOC) (similarly to L M ESH ). Although we do not provide an theoretical proof of (PSD) for L h , eigenvalues of L h , have always been positive through all experiments. The (PCON) property is not applicable to our framework. The pointwise convergence (DPCON) is observed experimentally and discussed in Section 4.1.

Experiments

Experimental Convergence

We first evaluate the multigrid convergence of our Laplace-Beltrami operator. We consider a unit sphere S 2 and three different smooth functions u : S 2 → R, namely cos(z), x 2 and e x (see. Fig. 3). Note that in our framework, u is extended in the normal direction to ũ as defined in Sect. 3. Let θ be the azimutal angle, and φ the polar angle. The spherical laplacian is then:

∆ S 2 u(θ, φ) = 1 sin 2 φ ∂ 2 u ∂θ 2 + 1 sin φ ∂ ∂φ sin φ ∂u ∂φ . (12)
We compute the Gauss digitization D h (S 2 ) of the sphere for decreasing grid step h. Since the elements of ∂ h S 2 does not interpolate the sphere, u is extended to ũ as defined in Sect. 3. We compute L COM BI and L h directly on ∂ h S 2 , but L COT and L M ESH on the associated Marching-Cubes triangulation. Since the vertices of this mesh coincide with the centroids of the surfels of ∂ h S 2 in D h (S 2 ), all these operators L COM BI , L M ESH and L h are evaluated at the same points.

For L h , we use the normal vector estimator described in [START_REF] Coeurjolly | Multigrid convergent principal curvature estimators in digital geometry[END_REF] to estimate the measure of the surfels. In addition, for both L h and L M ESH , the parameter t h is set to 0.1 × h 1 3 . As the discretization becomes finer, the standard deviation σ := √ 2t h of the Gaussian function decreases and the number of points within the standard deviation σ increases. The constant factor 0.1 is a scale term derived from the unit sphere that sets the kernel to 1/10 of the sphere.

For comparison, in order to mimic the setting of [START_REF] Belkin | Discrete laplace operator on meshed surfaces[END_REF], we have also considered the laplacian L P M ESH , which corresponds to L M ESH when the vertices of the Marching-Cubes are projected onto the sphere.

For all the above operators, we plot in Fig. 3 the l 2 and l ∞ error between the computed laplacian and the true spherical laplacian against the grid step h. First we observe that errors for L COT (in blue

) and L COM BI (in green ) are constant: clearly both operators are non-convergent. Non-convergence is also observed for L M ESH (in red ) but with lower errors. On the opposite, L P M ESH (in orange ) shows convergence behavior for both l 2 and l ∞ error, as expected in [START_REF] Belkin | Discrete laplace operator on meshed surfaces[END_REF]. Concerning L h (in purple ), experimental convergence holds for the three functions. For the periodic function cos(z), the convergence speed is slower than L P M ESH whereas for the non linear functions we can see that convergence speed is the same. Moreover, the l 2 error for L h tends toward L P M ESH , and its l ∞ error is close to the one of L P M ESH .

Shape approximation using eigenvectors decomposition

In this section, we consider the spectral analysis framework to process shapes geometry [START_REF] Lévy | Spectral Mesh Processing[END_REF]: Given a shape and its Laplace-Beltrami operator, we compute the eigenvalues and eigenvectors of the operator and project the geometry onto the eigenvector basis of the first k eigenvalues. More formally, given on operator L, we denote by e 1 , e 2 , . . . , e n its normalized eigenvectors and the matrix E whose columns are those eigenvectors. By λ 1 , λ 2 , . . . , λ n we denote the associated increasing eigenvalues where n is the number of rows of L. Given three input vectors (X, Y , Z) encoding the vertex positions in R 3 , we can approximate the input shape using a fixed number k of eigenvectors:

X (k) = E (k) (E (k) ) T X, Y (k) = E (k) (E (k) ) T Y , Z (k) = E (k) (E (k) ) T Z,
where E (k) is a matrix of size n×k containing the first k eigenvectors columnwise. We compute the eigen decomposition on D h (M ) for L COM BI , and L h in Fig. 4 on a bunny object (64 3 , 13236 eigenvectors). We illustrate the reconstruction for increasing number of eigenvectors k. For low frequencies (k ≤ 100), we observe that L h captures more geometrical details than L COM BI For k = 100, we clearly have a better approximation of both ears of the bunny shape. As k increases, both reconstructions converge to the original bunny shape (E (n) (E (n) ) T = I). In Fig. 5, we show the first 20 eigenvectors on a axis aligned cube.

Heat diffusion

In this section, we highlight an interesting isotropic property of L h compared to the combinatorial one. We compute a heat diffusion when the source is a Dirac in the center of a rotated cube face (heat diffusion is a preliminary step of [START_REF] Crane | Geodesics in heat: a new approach to computing distance based on heat flow[END_REF] to estimate geodesics on a manifold). We can derive an expression of g(x, t) from Eq.( 5): (I -tL)g(x, t) = g 0 (x), where I is the identity matrix, and L is the laplacian matrix associated to any discrete laplace operator. This diffusion only makes sense for small t and around the Dirac g 0 (x). As the computed heat diffusion decreases exponentially, we display the absolute value of its log. For L COM BI (first column in Fig. 6), both small and large values of t lead anisotropic estimations of the intrinsic metric due to the staircase effect of the face (concentric rhombi or concentric ellipses depending on t, see similar discussion in [START_REF] Crane | Geodesics in heat: a new approach to computing distance based on heat flow[END_REF]). When using Eq.( 13) with our matrix based operator L h (second column), even if numerical instabilities occur far from the Dirac depending on t, the intrinsic metric is perfectly estimated (concentric circles). Note that both L M ESH and L h already compute a diffusion without the approximation of Eq. ( 13). The third column shows the associated diffusion. Even if L M ESH does not have multigrid convergence properties, on this specific object, both methods provide good isotropic behaviors.

Conclusion and future-works

In this paper, we have investigated different discretization scheme of the Laplace-Beltrami operator on digital surfaces. The contribution is twofold: First, we have shown that classical schemes either do not asymptotically converge to the expected operator, or contain anisotropic artifacts when used for geometry processing tasks. Second, we have proposed a new Laplace-Beltrami operator that incorporates multigrid convergent surfel measures allowing us to have both an experimental multigrid convergence and isotropic properties on digital surfaces. A natural future-work consists in focusing on the multigrid convergence proof of the operator. In dimension 2, a preliminary proof has been derived using digital integration results from [START_REF] Lachaud | Properties of gauss digitized shapes and digital surface integration[END_REF] but it is still an open problem in dimension 3.

Fig. 3 .

 3 Fig. 3. Multigrid convergence graphs for various functions on S 2 the unit sphere. Both l2 error in plain line and l∞ in dashed line are displayed for LCOMBI , LMESH ,L P M ESH , and L h .

Fig. 4 .

 4 Fig. 4. Images of the reconstruction using an increasing number of eigenvectors k. (First row ) using LCOMBI , (second row ) with L h (r = 6 for [4] and t h = 3).

Fig. 5 .

 5 Fig. 5. Eigenfunctions are displayed on a simple cube with faces aligned with the grid axes (with a red to blue colormap and zero-crossing in white).

Fig. 6 .

 6 Fig.6. Heat diffusion on a cube aligned with R 3 axis. (First column) using LCOMBI , (second column) using L h and (third column, with t h = 4) using a direct diffusion. The rightmost picture shows staircases on the rotated cube.

Table 1 .

 1 Properties of various laplacians. See text for the description of each property.

		Ref. SYM LOC LIN POS PSD CON PCON DPCON
	MEAN VALUE	[8]			?	N.A.
	INTRINSINC DEL [2]		?	?	N.A.
	LCOMBI	[25]
	LCOT	[7, 9]
	LMESH	[1]	?
	L h	here	?	exp. ?	N.A.	exp.


23-24 March, 2017, London (pp. 66-68).

What do we see when we look at data? This recurrent question when confronted to Big Data, is largely answered by two different disciplinary visions dominating the debate during the last years, concluding that data are raw and Big Data is a hubris. While disciplines still treat just a part of them, we experience a wider and wider spreading of Big Data. The notion of "monuments of cyberspace" discussed here helps understand their peculiar nature and delimit related issues of method, of wealth and of experience.

In Knowledge Management, data are traditionally seen as a resource, a mine of information to be processed into knowledge. This view, dominant in the realm of the knowledge economy, conceptualised data according to the specific objective of the production of knowledge that is useful for the manager. In this context, however, data are defined in a negative way, being just "unstructured" [START_REF] Gandomi | Beyond the hype: Big data concepts, methods, and analytics[END_REF][START_REF] Bellinger | Data, information, knowledge, and wisdom[END_REF][START_REF] Earl | Knowledge Management Strategies: Toward a Taxonomy[END_REF]:

Data is raw.. It simply exists and has no significance beyond its existence (in and of itself). It can exist in any form, usable or not. It does not have meaning of itself [START_REF] Bellinger | Data, information, knowledge, and wisdom[END_REF].

More recently, along with the expansion of Big Data, researchers examining large amounts of data coming from digital platforms have been engaged in a process of identifying their meaning. Thus, we have entered a realm where data are to be understood as such and not only as a means to achieve a goal, a part of a procedure having an objective of production.

Yet, most research methodologies still describe data as an effect of the very action of the researcher. For instance, in quantitative research [START_REF] Yin | Case study research: design and methods[END_REF], data are collected on the basis of a research question, usually resulting from the identification of a gap in the literature. Hence, the structure of the data to be collected is normally designed by the researcher. Consequently, from the perspective of the research community, Big Data have been characterised as an "hubris", as they "are not the output of instruments designed to produce valid and reliable data amenable for scientific analysis" [START_REF] Lazer | The parable of Google flu: traps in big data analysis[END_REF].

Overall, it seems that the Knowledge Management approach, admitting that data exists beyond the action of the researcher appears to be more relevant to Big Data than the methodological approaches considering data as the result of a researcher's action. Still, management is not the only discipline that seeks value in them. Moreover, Big Data as a research object are yet illy defined. How to qualify their "raw" and "unstructure" nature ?

A unified view

To clarify the nature of contemporary data and, thus, the corresponding data situations that scholars and practitioners are called to address, this paper proposes the distinction of two conceptual levels: (1) resource versus result and (2) action versus research.

Thus, the problem of data can be analysed as shown in the Table 1.

Action Research

Resource Knowledge Management Unknown 1

Result

Unknown 2 Methodology

Table 1: Data situations.

As indicated by Table 1, we are today confronted with two new situations, two situations of unknown in relation to the well-known disciplinary views. The first situation concerns the condition where data are seen as resource for research, yet they do not originate from the action of the researcher and where their relevance to a research domain (research question or even related discipline) remain unknown. The second condition is when data are the result of action (of a system or a collectivity) but their usefulness remains to be discovered (we don't know what to do with them. . . ).

Treating the unknown data situations as monuments

To treat the unknown data situations we propose that an investigation of the experience of data generation is a prerequisite. We propose to treat data as monuments (e.g., similarly to the study of ancient manuscripts), positioning them in a sequence, in a sort of a chronicle, portraying this way the evolution and the ruptures in the process of their creation.

Thus, we the intrinsic mode of knowledge edification in digital settings can be described as "monuments of cyberspace" [START_REF] Chrysos | Monuments of cyberspace: Designing the Internet beyond the network framework[END_REF]. More specifically, each digital interaction is not only limited to an exchange of information, a transaction, nor is it always the effect of a mobilisation of already existent relations. A digital interaction is always inscribed on the technical substratum, the cyberspace. A common place is thus developed during an event of digital interaction. A monument of cyberspace can be visited long after its creation by unexpected visitors seeing the technical substratum as well as the traces left by the interaction in it.

Implications

The configuration of the four different data situations enables a better understanding of the unknown situations and opens the way for more relevant methods of data analysis. Moreover, the distinction of two contemporary situations of unknown helps avoid misunderstandings among disciplines enabling the identification of the specific data situation under study and mapping the field of research of this new phenomenon.

In addition, seeing data as monuments of cyberspace helps identify the dimensions to explore in the unknown situations. In the condition where data are the resource for research, the common place of their generation has to be further explored, as well as the chronicle of the interactions that occurred in it. In the condition where data are the result of an action, beyond its objective, the means of the action and specifically the design of the technical substratum has to be studied from a new perspective: the one of enabling new types of interaction, unforeseen during the initial design process of the setting used.

  Between disciplines and experience: Understanding Big Data as "monuments of cyberspace" Paris Chrysos, ISC Paris School of Management Proceedings of the International Conference 'Between Data and Senses; Architecture, Neuroscience and the Digital Worlds'. Organised by UEL, supported by RIBA, ARUP and the Museum of Architecture