Najah Chridi
email: chridi@loria.fr

Mathieu Turuani
email: turuani@loria.fr

Michael Rusinowitch

Towards a Constrained-based Verification of Parameterized Cryptographic Protocols

Although many works have been dedicated to standard protocols like Needham-Schroeder very few address the more challenging class of group protocols. We present a synchronous model for group protocols, that generalizes standard protocol models by permitting unbounded lists inside messages. In this extended model we propose a correct and complete set of inference rules for checking security properties in presence of an active intruder for the class of well-tagged protocols. Our inference system generalizes the ones that are implemented in several tools for a bounded number of sessions and fixed size lists in message. In particular when applied to protocols whose specification does not contain unbounded lists our inference system provides a decision procedure for secrecy in the case of a fixed number of sessions.

Introduction

Cryptographic protocols are crucial for securing electronic transactions. They rely on cryptographic functions to ensure security properties such as secrecy or authentication. The confidence in these protocols can be increased by a formal analysis in order to verify that the security properties are met at least at the logical level, that is, even when abstracting from the cryptographic functions and considering messages as first-order terms. Verification at the logical level is nevertheless a non-trivial task since cryptographic protocols are infinite state systems and for instance the set of potential messages that can be generated by an intruder is unbounded. Recently numerous works have been dedicated to the design of automated verification tools for cryptographic protocols. Such tools are often based on model-checking, modal logics, equational reasoning, and resolution theorem-proving (see e.g., [START_REF] Weidenbach | Towards an automatic analysis of security protocols in firstorder logic[END_REF][START_REF] Armando | Automatic sat-compilation of protocol insecurity problems via reduction to planning[END_REF][START_REF] Meseguer | Symbolic reachability analysis using narrowing and its application to verification of cryptographic protocols[END_REF]). Checking whether a protocol is flawed in the abstract Dolev Yao model [START_REF] Dolev | On the security of public key protocols[END_REF] can often be reduced to a constraint solving problem in a term algebra (modulo an equational theory). This constraint-based approach has proved to be quite effective on standard benchmarks and has also permitted the discovery of new flaws in several protocols (see e.g., [START_REF] Basin | An On-The-Fly Model-Checker for Security Protocol Analysis[END_REF]).

However to our knowledge it has never been applied to the more challenging group protocols. In fact very few formal verification results are available for such protocols. The difficulty relates to the fact that group protocols may perform an arbitrary number of steps since the group of communicating agents is a priori unbounded. This allows one to encode easily undecidable problems. Related work Several works have considered protocols with unbounded number of participants and recursive steps. The formal analysis of such protocols goes up with Paulson in [START_REF] Paulson | Mechanized proofs for a recursive authentication protocol[END_REF] who studied the Recursive Authentication (RA) protocol [START_REF] Bull | The authentication protocol[END_REF] for an unbounded number of participants using Isabelle/HOL theorem prover [START_REF] Paulson | Isabelle: A generic theorem prover[END_REF]. However if the protocol is defective there is no automatic mechanism to find the attack.

The validation of group protocols has been investigated in the CLIQUES project [START_REF] Steiner | Cliques: A new approach to group key agreement[END_REF], based on group Diffie-Hellman (A-GDH) protocols. Several analysis methods have been applied in this project, from manual to automatic ones. An interesting result has been obtained by Pereira and Quisquater [START_REF] Pereira | Some attacks upon authenticated group key agreement protocols[END_REF] who found several attacks on the CLIQUES suite and have shown that it is impossible to design a correct authentication group key agreement protocol built on A-GDH for a number of participants greater than three [START_REF] Pereira | Generic insecurity of cliques-type authenticated group key agreement protocols[END_REF]. Recently Kremer and al. [START_REF] Kremer | Proving group protocols secure against eavesdroppers[END_REF] have developped an automata-based approximation technique to analyse this class of protocols and check the absence of flaw in presence of a passive intruder. In [START_REF] Meadows | Formalizing gdoi group key management requirements in npatrl[END_REF], the NRL protocol analyser, which is based on a combination of model checking and theorem-proving techniques, has been adapted to handle the GDOI's protocols. Although Diffie-Hellman exponentiation has been encoded in the tool, it was not able to rediscover Pereira-Quisquater attacks on the CLIQUES suite [START_REF] Meadows | Extending formal cryptographic protocol analysis techniques for group protocols and low-level cryptographic primitives[END_REF]. Coral system [START_REF] Steel | Attacking group multicast key management protocols using CORAL[END_REF] has analysed an improved version of the multicast group key management protocol by Tanaka and Sato [START_REF] Taghdiri | A lightweight formal analysis of a multicast key management scheme[END_REF]. Two serious attacks have been found on this protocol. Coral has also discovered other attacks on Asokan-Ginzboorg [START_REF] Asokan | Key agreement in ad hoc networks[END_REF] and Iolus [START_REF] Mittra | Iolus: A framework for scalable secure multicasting[END_REF] protocols.

Some works have focused on the modelling of recursive computations performed by some participants (such as a server) in group protocols. In [START_REF] Küsters | Automata-based analysis of recursive cryptographic protocols[END_REF] tree transducers are introduced to model recursion and to allow the protocol participants to output structured messages. This work gives a decision algorithm for secrecy in the case of atomic keys and bounded message size in the Dolev Yao setting. However messages cannot be tested for equality without losing decidability. Similarly using composed keys or adding equational theories for XOR or Diffie-Hellman exponentiation in their model leads to undecidability. In [START_REF] Truderung | Selecting theories and recursive protocols[END_REF], Truderung introduces a class of Horn clauses to model the recursive behavior of participants. In this model protocol participants may receive messages of unbounded sizes, send multiple messages in a single step, compare and store messages. He gives a decision procedure to check whether protocols in this model satisfy secrecy properties. His algorithm is in NEXPTIME and is based on the derivation of an exponential bound on the size of minimal attacks. Hence this nice result is rather of theoretical flavour and is not suitable for an implementation. Only atomic keys are allowed for encryption. Moreover, Truderung's model cannot model some computations such as list mapping or functional symbol mapping. Note that non-atomic keys can be handled by our verification procedure (to be presented in the following sections). In [START_REF] Küsters | On the automatic analysis of recursive security protocols with xor[END_REF] an extension of [START_REF] Truderung | Selecting theories and recursive protocols[END_REF] has been proposed to handle XOR operator. Security can then be decided for a class of recursive protocols where principals are forbidden to XOR several messages (depending on messages) received from the network. Another extension [START_REF] Kürtz | Selecting theories and nonce generation for recursive protocols[END_REF] has been designed to model freshness of nonces and keys more accurately. Contribution. We present a synchronous model for parameterized protocols, that can be viewed also as an extension of classical protocols model to handle uniform lists of messages whose length is a fixed parameter. Such kind of protocols can be found in different domains. As example, we can cite web services protocols where lists of messages are commonly used. We can consider for instance, a basic action of an honest participant who received a messages list whose length is not known in advance and sends a message based on the informations he got:

{m 1 } k , .., {m i } k , .., {m n } k -→ f ( m 1 , .., m i , .., m n )
In this step, the participant decrypts the different messages composing the received list using the public key of the sender. He then builds the message to be sent using the list of messages m i and a function f . Note that this kind of protocols are used in web services where messages may contain unbounded lists of encrypted XML nodes. Moreover, our model can handle in particular, group protocols that have an unbounded number of participants, and thus, admit this number as a parameter. For this model we propose a complete and correct set of inference rules that allows one to check the security of well-tagged protocols in this class in presence of an active intruder. We show that relaxing the conditions on well-tagged protocols leads immediately to undecidability. Moreover, the proposed rules can model list mapping and allow non atomic keys. The class of protocols that we study admits tagged messages. Tagging basically avoids some unifications between messages that could be exploited for attacks. Several works on protocols have considered tagging techniques on messages as ours in order to enforce decidability. But these works ([BP03,RS03]) do not consider group protocols, or protocol with unbounded lists. Moreover in our case tagging is limited to messages that contains indexed variables, that is variables to be instanciated by items of unbounded lists. The other messages do not need to be tagged. Organization of the paper. We introduce in Section 2 our protocol model. We define attacks and show that their detection is undecidable. This motivate us to introduce the class of well-tagged protocols that is a good candidate for decidable security. In Section 3, we introduce auxiliary predicates and their semantics. They are meant to express message constructibility (from intruder knowledge) and they are used to build constraint system whose satisfiability is equivalent to the existence of attacks on a protocol. In Section 4, we introduce a set of simplification rules to reduce these constraints. Finally, in Section 5, we give an algorithm for applying these rules on a constraint system modelling a protocol security problem. Full details with proofs of correctness, completeness and decidability of normalized constraints are given in a technical report [Rep].

Motivating Example: Synchronous Group Protocols

As a motivating and running example, we introduce the Asokan-Ginzboorg group protocol which is an application level protocol. Let the group be of size n + 1 for n ≥ 1. The protocol describes the establishment of a session key between a leader (a n+1 ) and a random number n of participants (a i , i = 1..n). Indeed, the leader starts the communication by sending a key (e). As a response, each participant generates a symmetric key (r i ) and a contribution to the group key (s i ) and sends them encrypted under e. The group key would be f (s 1 , . . . , s n+1 ). towards a synchronous version. In the initial version, agents receive from the leader some message R i (for the agent i) and send the message S i . Messages R i have the same pattern, as well as messages S i . We abstract all these agents into a single agent : the simulator who, while receiving the list of R i represented by the variadic list constructor mpair( , )., sends another list containing S i represented by the same operator. This protocol is called synchronous protocol.

1. a n+1 → ALL : a n+1 , {e} p 2. a i → a n+1 : a i , { r i , s i } e 1 ≤ i ≤ n 3. a n+1 → a i : { s 1 , .., s n+1 } ri 1 ≤ i ≤ n 4. a i → a n+1 : a i , { s i , h(s 1 , . . . , s n+1 ) } K some i, K = f (s 1 , . . . ,
For our example, the simulator is called S and it simulates agents a 1 , . . . , a n . Agent L simulates the leader a n+1 . This way, we obtain below a protocol with a fixed number of steps to the expense of introducing a variadic list constructor mpair( , ):

1. L -→ S : mpair(t, l, {e} p ) 2. S -→ L : mpair(i, a i , { r i , s i } e ) 3. L -→ S : mpair(i, { mpair(j, s j ), s ′ } ri ) 4. S -→ L : mpair(i, a i , { s i , h( mpair(k, s k ), s ′ ) } f ( mpair(k,s k ),s ′ ) )
Note that All parametric lists have the same length n.

We extend the protocol model [START_REF] Rusinowitch | Protocol insecurity with a finite number of sessions, composed keys is np-complete[END_REF] in order to deal with parametric lists (whose length is the parameter). They are constructed with a new operator denoted by mpair( , ). The intuition is that a mpair message is equivalent to a list of messages built with the same pattern.

Names, Operators and Messages

Let X be a set of variables represented by capital letters. Let I be a countable set of index variables. Let -→ X be a set of symbols represented by overarrowed capital letters, disjoint from X . Let X I = {Y i s.t.

-→ Y ∈ -→ X and i ∈ I} be a countable set of (indexed) variables. Similarly, let C and -→ C be (disjoint) sets of symbols, represented by (overarrowed) lower case letters, and let C I = {c i s.t. c ∈ -→ C and i ∈ I}. Elements in C and C I are called constants. In this paper, terms can be (optionally) tagged by an index. That is, let -→ e ∈ -→ C be a symbol that we reserve for tagging operations only. Then a term is an element of T in the following language:

T s = {T } p T | {T } s T | h(T ) | T , T | mpair(I, T ) | X | X I | C | C I \{e i | i ∈ I} T = [e i , T s ] | T s with i ∈ I
where T s is by definition the set of untagged terms. The operators { } p and { } s represent asymmetric and symmetric encryptions respectively, , is a pairing operator, h is a hash function and mpair(i, t) is a symbolic representation of a list (or tuple) of terms, built from the common pattern t by iterating i along integers. The translation function defined in Section 2.2 gives the semantics of this operator. We denote by signature G the set of operators in T s . To simplify the syntax, in the following we will write t i instead of [e i , t], and call it a tagged term. We also omit the tag i of a term t i , whenever the tag i is not relevant to the discussion. We denote by T g the set of ground terms, i.e. any term t ∈ T with no variable in X or X I and no mpair symbol. Ground terms will be used to describe messages that are circulated in a protocol run. Given a term t we denote by V ar(t) (resp. Cons(t)) the set of variables (resp. constants) occuring in t. We denote by Atoms(t) the set V ar(t) ∪ Cons(t).

In order to represent a list of terms we iterate the pairing operator , . For instance to represent a, b, c, d we can use the term a, b, c, d and we shall write this term in a shorthand: a, b, c, d . However we do not assume any associativity property of pairing.

A substitution σ assigns terms to variables. A ground substitution assigns ground terms to variables. The application of σ to a term t is written tσ. These notations are extended to sets of terms E in a standard way: Eσ = {tσ|t ∈ E}. The set of subterms of t is denoted by Subterm(t). It is defined recursively as follows: If t is a variable or a constant then

Subterm(t) = {t}. If t = f (t 1 , .., t n ) or t = f (t 1 , .., t n ) i with t ∈ G, then Subterm(t) = {t} ∪ ∪ n i=1 Subterm(t i ). Note that
u is not considered as a subterm of u i . We denote by ≤ the subterm relation on T . We define the relation m over T × T as the smallest reflexive and transitive relation such that if t = f (t 1 , .., t n ) or t = f (t 1 , .., t n ) j with f = mpair, then for all i = 1, . . . , m we have t i m f (t 1 , .., t m ). Note that t m u implies t ≤ u.

Finally, we define two kind of Index-operations: replacements, used in the inference rules over constraints, and substitutions, used to define the solutions of constraints (See Section 3).

Definition 1 (Index-Replacement δ and Index-Substitution τ ). An Index-Replacement δ (resp. Index-Substitution τ ) is an application from I to I (resp. to non-negative integers) that is extended to indexed variables and constants with δ(

X i ) = X δ(i) and δ(c i ) = c δ(i) (resp. τ (X i ) = X τ (i) and τ (c i ) = c τ (i)
) and extended to terms and sets of terms in the natural way.

We will use the notations δ i,j (resp. τ i,j ) to denote the replacement (resp. substitution) of i ∈ I by j ∈ I (resp. j ∈ N). We also use δ k i,j to denote the replacement of i ∈ I by j ∈ I and the other indexes apart from i by k ∈ I. We define the set of indexes occuring in a term as follows:

Definition 2 (Term Indexes). Given a term t ∈ T , we denote by V ar I (t) the set of indexes in t, recursively defined as follows:

V ar

I (mpair(i, t)) = V ar I (X) = V ar I (c) = ∅ with X ∈ X and c ∈ C V ar I (X i ) = V ar I (c i ) = {i} with X i ∈ X I and c i ∈ C I V ar I (f (t 1 , .., t n )) = V ar I (t 1 ) ∪ .. ∪ V ar I (t n ) otherwise V ar I (t i ) = V ar I (t) ∪ {i}

Protocol Specification, Intruder and attacks

A protocol is given by a set of principals and a finite list of steps for each. We associate to each principal A a partially ordered finite set (W A , < WA ) steps R ι ⇒ S ι where R ι is an expected message and S ι his reply. Init and End are fixed messages used to initiate and close a protocol session. Our notion of correct execution of a protocol session (or protocol run) follows [START_REF] Rusinowitch | Protocol insecurity with a finite number of sessions, composed keys is np-complete[END_REF].

We follow the intruder model of Dolev and Yao [START_REF] Dolev | On the security of public key protocols[END_REF]. The actions of the intruder are simulated by a sequence of rewrite rules on sets of messages. These rules are defined as follows. We note -→ * DY their reflexive and transitive closure. We define a predicate Dy. This predicate checks whether a message can be constructed by the intruder from some known messages.

Decompositition Rules L d ( a 1 , . . . , a n ) : a 1 , . . . , a n → a 1 , . . . , a n , a 1 , . . . , a n L d ({a} p K ) : {a} p K , K -1 → {a} p K , K -1 , a L d ({a} s K ) : {a} s K , b → {a} s K , b, a L d (t i ) : t i → t Compositition Rules L c ( a 1 , . . . , a n ) : a 1 , . . . , a n → a 1 , . . . , a n , a 1 , . . . , a n L c ({a} p K ) : a, K → a, K, {a} p K L c ({a} s K ) : a, b → a, b, {a} s K L c (t i ) : t → t i for any i ∈ I Definition 3 (Non-Redundant Derivation). Let D = E 0 -→ L1 . . . -→ L l E l be a derivation such that u ∈ E l and ∀i < l u / ∈ E i . We call u a goal of D and we note D as D u (E 0 ). D is a non redundant derivation if ∀i ∀t ∈ E i , if L c (t) ∈ D then ∄L d ( ) ∈ D that
Definition 4 (Dy, Dy c and Dy d ). Let E and K be sets of ground terms and t be a ground term such that there exists D ∈ N RD with goal t without decomposing any term of K. Then, we say that t is forged from E and we denote it by t ∈ Dy(E, K).

Moreover, if D = D ′ .L c (t) then t ∈ Dy c (E, K), otherwise t ∈ Dy d (E, K).
We interpret the mpair( , ) operator in the standard Dolev-Yao signature by defining a translation function that replaces any mpair by a sequence of pair applications. The number of such applications is given as a parameter e to the translation function. The integer represents the common length of lists of terms represented by any mpair(, ).

Definition 5 (Translation of terms). Let T DY be the set of terms without any mpair(, ). Given any integer e and function f = mpair, the function -e from T to T DY is defined as follows: We can now define attacks on protocols in our model, based on Dy predicate.

mpair(i, t) e = τ i
Definition 6 (Attack). Given a protocol execution P = {R i ⇒ S i |i = 1, .., k}, a secret Sec and assuming the intruder has as initial knowledge S 0 , an attack is described by a ground substitution σ, an Index-Substitution τ , and an integer e, such that ∀i = 1, . . . , k, we have: R i τ σ e ∈ Dy(S 0 e , S 1 e τ σ, . . . , S i-1 e τ σ, ∅) Sec e ∈ Dy(S 0 e , S 1 e τ σ, . . . , S k e τ σ, ∅)

Undecidability and Well-Tagged Protocols

Unfortunately, the insecurity problem (i.e. the existence of an attack) is undecidable in the general case. This can be shown by encoding Post Correspondance Problem (PCP) with two letters. Note that this encoding requires only atomic keys.

Definition 7 (PCP protocol). Let J = {(α 1 , β 2 ) , .., (α p , β p )} be an instance of PCP on the alphabet {a, b}. We define the protocol specification P (J) coding J as the following, with C = {a, b, 0, t, u},

X = {Z}, -→ X = -→ A , -→ B , -→ X , -→
Y and only one honest participant :

1. Init ⇒ a, b, 0, { 0, 0 } t 2. mpair(i, A i , B i ) ⇒ mpair(i, { A i , B i } t ) 3. mpair(i, { X i , Y i } t ) ⇒ mpair(i, { α 1 (X i ) , β 1 (Y i ) } u ), . . . mpair(i, { α p (X i ) , β p (Y i ) } u ) 4. mpair(i, { A i , B i } u ), { Z, Z } u ⇒ Sec
Theorem 1. An instance J of PCP has a solution iff P (J) has an attack on Sec.

The proof of this theorem is detailed in technical report [Rep]. We will therefore introduce the class of Well-Tagged protocols for which decidability is expected.

To do this, we first introduce the notion of autonomy: Definition 8 (Autonomy). A term mpair(i, u) is sayed to be autonomous when V ar I (u) ⊆ {i}. A term t ∈ T DY is autonomous if #V ar I (t) ≤ 1 and ∀t ′ < t, t ′ is autonomous. A protocol P = {R ı ⇒ S ı |ı ∈ J} is autonomous iff for all ı ∈ J, R ı and S ı are autonomous and V ar I (R ı ) = ∅ and V ar I (S ı ) = ∅.

For instance, the term t = mpair(i, mpair(j, {a i } cj )) is not autonomous. We remark that the autonomy property alone is not enough to guarantee decidability, since the PCP protocol of Definition 7 is autonomous.

Definition 9 (Well-Tagged protocols).

P = {R ı ⇒ S ı |ı ∈ J} is Well- Tagged iff: 1. ∀ı ∈ J, ∀X i ∈ X I ∩ Subterm(R ı ) ∩ ı ′ <ı Subterm(R ı ′ ), X i is tagged; 2. ∀ı ∈ J, ∀X i ∈ X I ∩ Subterm(S ı ), X i is tagged; 3. ∀ı ∈ J, ∀t = f (s 1 , .., s k ) ∈ Subterm(S ı ) with f = mpair, if ∃j = 1..k s.t. s j
is tagged, then t is tagged too; 4. P is autonomous.

In this definition, Conditions 1 and 2 state that any indexed variable of the protocol must be tagged, except for its first occurrence w.r.t. the partial step ordering. Moreover, Condition 3 (when combined with Condition 2) states that, for any subterm t of some acquired intruder knowledge S ı , if an indexed variable is accessible from t by decompositions without opening any mpair, then t must be tagged. Note that as a consequence of mpair autonomy, an indexed variable X i can only be tagged by its index (as in X i i ) or untagged. The idea underlying the tagging of variables is to add enough information on terms in mpair so that the protocol cannot be used to test or ensure relations between elements of the same mpair, such as ∀i = 2..n, ∃i ′ = 1..n s.t. X i = f (X i ′ ). This is precisely the kind of relations that the encoding of PCP is able to exploit. Thus, adding tags to the PCP-encoding protocol will generate a new protocol that cannot be run.

Constraints for Protocol Verification

We will use a symbolic constraint system to represent all runs of a protocol given a step ordering. This system uses (universal or existential) quantifiers on index variables and includes an (implicit) universal quantification on the number n of elements in any mpair. Before defining our constraint system, some basic notions have to be introduced. For terms s, s ′ (resp. sets of terms E, E ′ ) we note s ∼ s ′ (resp. E ∼ E ′ ) if they are equal once we erase their tags.

Definition 10 (Relation L E for accessible subterms). We consider a relation on T × 2 T × 2 T × T . We write s L E t for s, t terms in T and E and L finite subsets of T . Note that this can be used for T s too. This relation is defined as the smallest relation such that:

t ∅ ∅ t ∀t ∈ T s L E t iff s ′ L ′ E ′ t ′ where s ∼ s ′ , t ∼ t ′ , E ∼ E ′ , L ∼ L ′ If {m} p k L E t then m L ′ E,k -1 t where L ′ = L ∪ {{m} p k } If {m} s b L E t then m L ′ E,b t where L ′ = L ∪ {{m} s b } If t 1 , . . . , t n L E t then ∀i ≤ n, t i L ′ E t where L ′ = L ∪ { t 1 , . . . , t n } If m L E ′ t and E ′ ⊂ E then m L E t
We write u E t when u L E t for some L. Note that by construction, u E t implies u ∈ Subterm(t). We say that u is a subterm of t that is accessible, i.e. can be obtained from t by decompositions using keys in E. For simplicity, we note b1,..,b k instead of {b1,..,b k } . We define also the set of strict accessible terms by s< F t (resp. s< L F t) if s F t (resp. s L F t) and s = t. Given t L F u or t < L F u), we call length of t L F u or (resp t < L F u) the number of elements of L. Before defining constraint systems, we need to introduce the environment and elementary constraints.

Definition 11 (Environment). We call an environment a finite set of equalities X = u whose left-hand sides are variables (X ∈ X ∪ X I ). We usually denote it by E.

Definition 12 (Elementary Constraint). An elementary constraint is an expression

(t = t ′ ), (t ∈ F orge(E, K)), (t ∈ F orge c (E, K)), (t ∈ Sub(t ′ , E, E, K)), or (t ∈ Sub d (t ′ , E, E, K)) with an environment E, t, t ′ ∈ T and E ⊂ T .
An elementary constraint represents a basic relation on terms. Indeed, t ∈ F orge(E, K) if the term t is derivable from the knowledge E without decomposing elements of K; t ∈ F orge c (E, K) if t is derived by composition; t ∈ Sub(t ′ , E, E, K) if t is an accessible subterm from t ′ with knowledge E with none of the intermediate terms between t and t ′ in K, and this modulo replacements using equations of E; t ∈ Sub d (t ′ , E, E, K) if t is accessible by decomposition of t ′ , also modulo replacements using E; and t = t ′ if t and t ′ are equal.

Definition 13 (Negative Constraint). A negative constraint is an expression

(∀i X m = u) or (X m / ∈ F orge c (E, K)) with X m ∈ X I , u ∈ T , E, E ′ ⊂ T and i ∈ V ar I (u).

The set of solutions of a constraint S, denoted by |[ S ]| e

τ where e is a value of n and τ is an Index-Substitution is a set of ground substitutions to be defined in the following. We define GS to be the set of all ground substitutions. Definition 14 (Solutions of an Elementary Constraint).

|[ t = t ′ ]| e τ = {σ ∈ GS | t e τ σ = t ′ e τ σ} |[ t ∈ F orge(E, K) ]| e τ = {σ ∈ GS | t e τ σ ∈ Dy(E e τ σ, K e τ σ)} |[ t ∈ F orge c (E, K) ]| e τ = {σ ∈ GS | t e τ σ ∈ Dy c (E e τ σ, K e τ σ)} |[ t ∈ Sub(w, E, E, K) ]| e τ = {σ ∈ GS | ∃u ∃F, L s.t. u L F w e τ, F σ ⊆ Dy(E e τ σ, K e τ σ),
Lσ ∩ K e τ σ = ∅, and either uσ = t e τ σ or ∃v, δ, k, τ ′ s.t

              either u ∈ X , (u = v) ∈ E, δ = ∅, and τ ′ = τ or ∃ -→ Z ∈ -→ X , i, j ∈ I s.t u = Z i τ ′ , (Z j = v) ∈ E, or δ = δ k j,i , τ ⊆ τ ′ and Dom(τ ′ ) = Dom(τ ) ∪ {k, i} k, i / ∈ V ar I ({t, w, E}), uσ / ∈ Dy c (E e τ σ, K e τ σ) uσ = v e δτ ′ σ, and σ ∈ |[ t ∈ Sub(vδ, E, E, K) ]| e τ ′ ∪{k, i}} |[ t ∈ Sub d (w, E, E, K) ]| e τ is defined in a similar way as |[ t ∈ Sub(w, E, E, K) ]| e τ
with the difference that for the first case (when uσ = t e τ σ), we have u < L F w e τ .

Definition 15 (Solutions of a Negative Constraint).

|[ X m / ∈ F orge c (E, K) ]| e τ = GS \ |[ X m ∈ F orge c (E, K) ]| e τ |[ (∀i X m = u) ]| e τ = GS \ x=1...e |[ X m = u ]| e [i←x],τ
We describe our constraint system by blocks in the following way:

Definition 16 (Constraint System). First, we define a constraints block B as a conjunction of constraints together with an environment E: B = (ctr 1 ∧ . . . ∧ ctr l , E). We will sometimes handle blocks as set of elementary or negative constraints for ease of notations. For instance we write c ∈ B to express that the elementary constraint c is a conjunct of B.

We can now define the constraint system that we will use to represent protocol runs. Given two finite lists of index variables Q = i 1 , . . . , i k and R = j 1 , . . . , j l , we write the quantifier prefix ∀i 1 • • • ∀i k ∃j 1 • • • ∃j l in short: ∀Q ∃R. A constraint system, denoted by S, is a disjunction of blocks with a quantifier prefix: S = ∀Q ∃R (B 1 ∨ . . . ∨ B p ) Now, we define the set of solutions of the constraint system as follows:

Definition 17 (Solutions of the Constraint System). Consider a constraint system S, B i , for i = 1 . . . p (blocks of S) and ctr i,j , for j = 1 . . . l i (constraints of the block B i ) given in Definition 16. The set of solutions of a constraint system CS is defined inductively using the following cases:

|[ ∀i S ]| e τ = ∩ x=1,...,e |[ S ]| e [i←x],τ |[ S ]| e τ = i=1...p |[ B i ]| e τ |[ ∃i S ]| e τ = ∪ x=1,...,e |[ S ]| e [i←x],τ |[ B i ]| e τ = j=1...li |[ ctr i,j ]| e τ
The idea will be to use a constraint system based on blocks to represent all possible ways the intruder can construct a list of terms represented by an mpair. Roughly, there will be one block in the system for each way. Note also that blocks are extended to admit labeled constraints:

Notation 1 (labels of constraints) A constraint ctr may be equiped with a label (ctr) m or (ctr) sm or (ctr) f to denote respectively a master constraint or a submaster constraint or a final constraint. The two first labels allow us to keep track of the "official" formal value of some indexed or non-indexed variable. For example, we will prove that we have exactly one master constraint for every indexed variable in each block, and we will use master or sub master constraints to instanciate variables when needed; The third label will be used to prevent any further rewriting on some constraint. We introduce also the notation (ctr) to refer to labeled or non labeled constraint. The solutions of labeled constraints are the solutions of the constraints obtained by removing labels.

To simplify the use of (sub)master constraints, we group them into sets:

Definition 18 (Set of (sub)master constraints). Let S = ∀Q∃R B 1 ∨..∨B p be a constraint system, Y ∈ X , W ⊆ X and -→ X ∈ -→ X . We define M(S, -→ X ) = {ctr | ∃i (ctr) m ∈ B i and ∃j ∈ Q such that ctr = (X j ∈ F orge c (E, K)) or ctr = (X j = u)} and ∀i SM(B i , W ) = {ctr | (ctr) sm ∈ B i , ctr = (Y = u) and Y ∈ W }. Also, SM(B i , Y ) = SM(B i , {Y }).
When S is clear from the context we omit it in M(S, -→ X ) and write simply M( -→ X ).

Normalisation of a Constraint System

In this section we present the rules applied in the normalization function over constraint systems. The result of applying a rule is put in disjunctive normal form and existential quantifiers are moved up to the prefix of the system using first order logic. These rules are organized in six groups G 1 , . . . , G 6 .

G 1 aims at maintaining syntactic properties over a constraint system. Some rules handle labelling of master constraints by adding new labels or transfering existing ones. We call them respectively the labelling and the label transfer rules. Other rules format constraints in order to get preferably variables on their left hand-side, or replace indexed variables by non-indexed ones.

G 2 contains the F orge rules. Here, we have for example :

t ∈ F orge(E, K) -→ t ∈ F orge c (E, K) ∨ w∈E t ∈ Sub(w, E, E, K)
This generic rule illustrates the two possible ways for forging a term t: either by composing or by decomposing one of the knowledge. The other rules enumerate all possible ways a term can be composed by the Intruder: we have exactly one rule for decomposing each kind of operator in the signature G. In particular, there is one rule for the mpair operator where mpair autonomy is used to justify the quantification.

G 3 contains the Sub rules. These rules are similar to the F orge ones, but they decompose Intruder knowledge. In this group, there is a generic rule :

t ∈ Sub(u, E, E, K) -→ (t = u) if u ∈ K (t = u) ∨ (t ∈ Sub d (u, E, E, K)) otherwise
This rule follows precisely the intruder deduction rules: a term t is an accessible subterm of u iff t = u or there exists a direct subterm u ′ of u, derivable from u, with t being an accessible subterm of u ′ . Therefore, there exists exactly one rule for decomposing each kind of operator in G (apart from variables and constants).

G 4 encodes unification algorithm for terms in our system, and thus, the resolution of equality constraints. These rules simply consist in testing recursively the compatibility of each top operator in each term. Therefore, the only equality constraints remaining after an iteration of these rules are those assigning a value to a variable, i.e. X = u with X ∈ X ∪ X I . The goal of G 5 and G 6 is to replace variables by managing interactions between constraints.

G 5 aims at replacing variables by their value, inside one block and independently of other blocks. Consequently, these rules only consider multiple occurences of the same variable, with the same index in case of indexed variable. For instance, the interaction between two Equality constraints is managed by Rules 26 and 27 for respectively indexed and non indexed variables.

B ∧ (X i = u) ∧ (X i = v) -→ B ∧ (X i = u) ∧ (X i = v) ∧ u = v (26) (X = u) sm ∧ (X = v) -→ (X = u) sm ∧ (u = v) (27)
Besides, interaction for F orge and Sub constraints is managed as follows:

A ∈ F orge c (E, K) ∧ t ∈ Sub d (A, E ′ , E, K) -→ ⊥ where A ∈ X ∪ X I
This rule says that it is not necessary to decompose a variable. While a bit more complex than expected, our semantics of Sub d (, , ,) has been defined to prevent useless actions like this, thus ensuring the validity of G 6 . G 6 generalizes G 5 by allowing variable replacements from one block to an other one. Given a constraint containing a variable X i that must be replaced by its value, we enumerate a finite number of "candidate" terms representing all possible values of this variable according to the whole constraint system. These values are provided by master constraints. For instance, the interleaving between a Sub constraint with contraints of other blocks leads to the following rule:

t ∈ Sub d (X m , E ′ , E, K) -→ (Xi = u)∈E ∃k ′ t ∈ Sub(uδ, E ′ , E, K) ∧ (X m = uδ) f ∧X m / ∈ F orge c (E ′ , K)with δ = δ k ′ i,m
In this rule, only the case where master constraints are equality ones are taken into account since the interleaving with F orge ones leads to ⊥. This rule adds an (extra) equality representing the master constraint it used, but labeled final to prevent further reductions on it. It adds also negative constraints to eliminate the case of F orge master constraints.

Rules System with Tags We have defined our Rules system for constraints on untagged terms. Nevertheless, these rules deal also with tagged terms following the definition of our signature. For example, the equality constraint (X i ) i = u leads to ⊥ when u is untagged, since (X i ) i = [e i , X i ].

Definition 19 (Solved Constraint). A solved constraint is of type:

(X i = u) , (X = u) sm , (X i ∈ F orge c (E, K)) , X ∈ F orge c (E, K), (Y / ∈ F orge c (E, K)), or (∀j X i = u) where X ∈ X , Y ∈ X ∪ X I , X i ∈ X I , u ∈ T , j ∈ V ar I (u), E ⊂ T and K ⊂ T .
We will prove that at each step of our algorithm, the normalized constraint system contains only solved constraints.

Application to the Asokan-Ginzboorg Protocol

Let us Focus on the following step of the Asokan-Ginzboorg specification:

mpair(i, L, {E i } p ) ∈ F orge(E 1 , ∅) where E 1 = {mpair(t, l, {e} p )}. -→ (mpair(i, L, {E i } p ) ∈ F orge c (E 1 , ∅)) ∨(mpair(i, L, {E i } p ) ∈ Sub(mpair(t, l, {e} p ), E 1 , ∅, ∅)) by G 2 -→ * ∀i ∀j ((L ∈ F orge(E 1 , ∅) ∧ (E i = e) m ) ∨ ((L = l) sm ∧ (E i = e) m ) ∨((L = l) sm ∧ (E j = e) m )) by G 1 5 Well-tagged Protocol Verification
We introduce here the verification algorithm and the results that state the correctness and the completeness of the inference rules of Section 4 and decidability for protocols without mpair(, )'s and indexed variables. Given a set R of inference rules and a formula F we say that R(F ) is a closure of F by R if it is derived by a finite number of applications of rules in R and no rule can be further applied to R(F ). First of all, we define a reduction of equalities chain in an environment:

Notation 2 (⌈E⌉) We note ⌈E⌉ the closure of E by the following rules:

X = Y ∧ Y = u -→ X = u ∧ Y = u X i = Y i ∧ Y j = u -→ X i = uδ k ′ Q,i ∧ Y j = u for X, Y ∈ X ∪ X I , X i , Y j ∈ X I and u / ∈ X ∪ X I .
Second, we introduce the normalisation function denoted by S → (S)  . A normalisation of a constraint system can be defined when some closures can be computed as follows using a subset of the inference rules. This normalization operates in two main phases: Definition 20 (Normalization function). Let S be a block system. We denote by SR the whole set of inference rules except the labelling rule. We assume that we can compute: Phase 1: S 1 , a closure of S through SR except Rules 26 and 27; Labelling: S 2 , a closure of S 1 through the labelling and label transfer rules; Phase 2: S 3 , a closure of S 2 through SR. This closure is denoted by (S)

 .

The Labelling step adds labels for creating master constraints, making sure to always favour labelling of equality constraints to a forge constraints. While Phases 1 and 2 are similar by the rules they use, their behaviors differ: when used in our algorithm for a step R i ⇒ S i , Phase 1 will never use any constraint interleaving rule with a master constraint M( -→ X ) with L(X) = E i-1 . This means that during Phase 1, the variables with maximum level cannot be replaced yet from a block to an other, simply because none of them have master constraints yet. However, the second phase do not have this limitation. The verification algorithm is the following: Algorithm 1 Let P = {R i ⇒ S i |i = 1..k} be Well-Tagged, Sec ∈ T , R k+1 Sec and S 0 ⊂ T g . 1. Let CBS 0 ∀Q ∃R ⊤, with Q = R = ∅, be the initial constraint system. 2. For i from 1 to k + 1 :

(a) Assume that CBS i-1 ∀Q∃R B 1 ∨ B 2 ∨ .. ∨ B p ; (b) Let ctr i R i ∈ F orge(S 0 , S 1 , .., S i-1 , ∅); (c) Let E i = - → X M(CBS i-1 , -→ X ) and for all j = 1, ..p, X, Y ∈ X ∪ X I , E i,j = ⌈E i ∪ SM(B j , X )⌉\{(X = Y )}; (d) Let CBS i ( ∀Q∃R (B 1 ∧ ctr i , E i,1 ) ∨ .. ∨ (B p ∧ ctr i , E i,p ) )  3. Test Satisfiability of CBS k+1 (return insecure iff satisfiable).

Note that sets E i , E i,j denotes respectively the set of master constraints for vector variables and the set of submaster constraints for variables of block B j , both with variables of level strictly included in E i-1 . Notation ⊤ represents true. The algorithm chooses a "possible" protocol run represented by π, and tests if after this run Sec is derivable by the intruder for some length e of mpair(, ). We test this by increasing the initial constraint system CBS 0 with each protocol step successively, and by normalising the resulting constraint system at each step. This step-by-step normalisation is required by our inference rules which assumes that master and sub-master constraints for previous steps have been already computed. Since we will show that the normalisation preserves the semantics, |[ CBS k+1 ]| e ∅ = ∅ iff the protocol run defined by π has an attack. Thus, assuming the next lemmas, the correctness of the algorithm follows: If Algorithm 1 terminates for a subclass of Well-Tagged protocols, then the satisfiability of a normalised constraint system resulting from a protocol in this subclass is decidable.

Theorem 2. (Analysis of Well-tagged Protocols) If Algorithm 1 terminates for a subclass of Well-Tagged protocols, then, the insecurity problem is decidable for this subclass.

Full proofs can be found in Technical report [Rep]. Moreover, it is worth to notice that our algorithm always terminates for protocols without mpair(, ), thus showing that our procedure is an extension of protocol analysis in the basic case. The proof of termination for protocols without mpair(, ) and without indexed variables can be found in Technical report [Rep].

Conclusion and Further Works

We have proposed an extension of the constraint-based approach in symbolic protocol verification in order to handle a class of protocols (the well-tagged ones) which admit unbounded lists in messages. This class can be used to model in particular interesting group protocols. We conjecture that adding adequate control on our constraint simplication rules allows one to obtain termination of the constraint normalization process and therefore to derive a decision procedure for an interesting subclass of well-tagged protocols.

Fig. 1 .

 1 Fig. 1. Transformation to a Synchronous Protocol

Lemma 1 .∅

 1 (Correctness and Completeness of Normalization) Let CBS i and ctr i (i = 1..k + 1) be as in the verification algorithm, for some Well-tagged protocol P . Then for all e,|[ CBS i-1 ∧ ctr i ]| e ∅ = (CBS i-1 ∧ ctr i )  eLemma 2. (Satisfiability of normalised form)

Leader Agent 1 Agent i Agent n ... ... Agent 1 Agent i Agent n Leader Initial Version ... ... Agent n Simulator Synchronous Version

				s n+1 )
	Since here, modulo index renaming, all members of the group have identical
	actions, we are going to abstract them in one agent as it shows Figure 1.1.
	This Figure shows the transformation between an initial version of a protocol
	R1	S1
	Ri	Si	mpair(i, Ri)	mpair(i, Si)
	Rn	Sn

  generates t, and if ∃L d ( ) ∈ D generating t then L c (t) / ∈ D. We denote by N RD the set of non-redundant derivations. Remark 1. For each derivation D t (E), there exists a non redundant derivation D ′ t (E). Indeed, D ′ is obtained by elimination of L c (t) if L d ( ) ∈ D where L d ( ) generates t and by eliminating each L d ( ) ∈ D that generates t if L c (t) ∈ D.

⋆ This work was supported by the FP7-ICT-2007-1 Project no. 216471, "AVANTSSAR: Automated Validation of Trust and Security of Service-oriented Architectures"(www.avantssar.eu) and by CPER Project SeCoManet.


Vincent Fortineau

Maria Makarov

Pedro Rodriguez-Ayerbe

Isabelle A Siegler

Towards a seamless experimental protocol for human arm impedance estimation in an interactive dynamic task

The estimation of the human endpoint impedance interacting with a physical environment provides modelling insights both for the field of human movement science and for the design of innovative controllers for collaborative robotics based on physical human-robot interaction. Most of the existing human impedance estimation methods described in the literature rely on controlled environments closer to a laboratory than to an industrial setting. In this paper, a force perturbation method is proposed without any specific requirements on neither the force nor the position trajectories and without additional sensors placed on the human. The method is illustrated through an experimental study on a benchmark interaction task. The task was selected to be sufficiently variable for preventing the use of average trajectories and the perturbations used for the estimations are sufficiently low amplitude and short to not deteriorate the human performance in the task. The obtained impedance model parameters are in adequacy with the ones found in the literature, suggesting the validity of the proposed approach. The low constraints on the experimental settings make it applicable even out of the lab, for example, in humanrobot collaboration in manufacturing environments, to adjust the robot behaviour to changes in the state of the operator (fatigue, stress) or to adjust to a variety of human operators with different interaction strategies.

I. INTRODUCTION

Humans are capable of tailoring their limbs' dynamic properties [START_REF] Tsuji | Bio-mimetic impedance control of robotic manipulator for dynamic contact tasks[END_REF] to improve the performance during the exploration of a new task [START_REF] Burdet | The central nervous system stabilizes unstable dynamics by learning optimal impedance[END_REF] and achieve stable behaviour with a wide range of environments. Understanding and reproducing some properties of the humans' interactions during various tasks is thus relevant for the design of innovative robot control.

The mechanical impedance of human limbs has been studied for decades under the control prism [START_REF] Mussa-Ivaldi | Neural, mechanical, and geometric factors subserving arm posture in humans[END_REF] and is still an active research topic [START_REF] Fang | Online Model Based Estimation of Complete Joint Stiffness of Human Arm[END_REF], [START_REF] Phan | Estimating Human Wrist Stiffness during a Tooling Task[END_REF]. The notion of impedance relates kinematics to forces, allowing the modelling of human movements during physical interactions. Elementary models, such as mass (M or I for inertia) -spring (K) -damper (B, hence the acronym KBM or KBI) have been shown to approach human behaviour [START_REF] Tsuji | Bio-mimetic impedance control of robotic manipulator for dynamic contact tasks[END_REF].

The mechanical impedance cannot be directly measured and requires perturbing the limbs for estimations. Thus, the environment often has to be very well mastered to identify the impedance parameters. Two configurations can be easily distinguished: on one hand, static or isometric tasks [START_REF] Fang | Online Model Based Estimation of Complete Joint Stiffness of Human Arm[END_REF], [START_REF] Dolan | Dynamic and loaded impedance components in the maintenance of human arm posture[END_REF], where the unperturbed pose is well known and on the other vincent.fortineau@centralesupelec.fr hand, dynamic tasks, where ongoing movements make the estimation of the unperturbed trajectory more complicated. The choice of the perturbations nature, either force or position, is often dictated by kinematics knowledge. When the trajectory can be estimated online [START_REF] Burdet | A method for measuring endpoint stiffness during multijoint arm movements[END_REF], position perturbations enable stiffness identification alone. The introduction of force perturbations requires prior information on neither the movement nor the force shape. In [START_REF] Erden | Hand Impedance Measurements During Interactive Manual Welding With a Robot[END_REF], hand impedance estimations were conducted among welders, which brightly allowed the knowledge of the unperturbed trajectories. In other situations where the trajectory could not be perfectly known beforehand [START_REF] Burdet | A method for measuring endpoint stiffness during multijoint arm movements[END_REF], [START_REF] Bennett | Time-varying stiffness of human elbow joint during cyclic voluntary movement[END_REF], [START_REF] Abe | Modulation of elbow joint stiffness in a vertical plane during cyclic movement at lower or higher frequencies than natural frequency[END_REF], the movement's repeatability allowed imperfect but sufficient knowledge on the unperturbed trajectories.

Another family of methods proposed in the literature relies on measurements of muscle activity using electromyographic sensors. Impedance estimations based on such measurements are reported in an isometric task [START_REF] Wang | sEMGbased estimation of human arm force using regression model[END_REF] or right after its completion [START_REF] Castellini | sEMG-based estimation of human stiffness: Towards impedance-controlled rehabilitation[END_REF]. Dynamic tasks make it more complicated to differentiate impedance from the task muscles activation.

In previous work at the intersection of robotics, control and human movement sciences [START_REF] Avrin | The self-organization of ball bouncing[END_REF], the authors studied the human central nervous system during cyclic movements in a simulated benchmark task. They were able to reproduce human-like performances using a Central Pattern Generator based control structure. An experimental test-bed was designed for impedance measurements using interactive robotics [START_REF] Fortineau | Interactive robotics for human impedance estimation in a rhythmic task[END_REF].

The present paper aims to complete previous work by adding an impedance model of the human arm during the task. Dynamic considerations should allow finer modelling of human behaviour during a cyclic task, notably by observing possible impedance variations within a cycle or across individuals with different strategies. Such an examination provides insights for the human-robot collaboration and teleoperation as well as for control strategies design for robot interacting with various physical environments.

This paper proposes to identify the human endpoint impedance through an interactive robot [START_REF] Fortineau | Interactive robotics for human impedance estimation in a rhythmic task[END_REF], in conditions where both the restitution force and position (the differential between unperturbed and perturbed) can neither be measured, nor estimated with sufficient precision from average trajectories only. The new optimisation-based estimation of the unperturbed force trajectory proposed in this paper significantly improves the impedance parameter estimation. Experimental validation shows that the identified parameters are consistent with previously documented values (Section V-B3). Thus, the proposed estimation method allows less constrained experimental conditions and can be more easily replicated outside a laboratory for impedance estimation in a dynamic task.

The experimental test-bed is outlined in Section II, with the experimental data sets conditions. The issue of force trajectory estimation raised in [START_REF] Fortineau | Interactive robotics for human impedance estimation in a rhythmic task[END_REF] is then addressed in Section III. The impedance identification methodology is described in Section IV. Finally, the results using this implementation are presented in Section V.

II. EXPERIMENT

The experiment is briefly recalled in this section. For further details, please confer to previous work, [START_REF] Avrin | The self-organization of ball bouncing[END_REF], [START_REF] Fortineau | Interactive robotics for human impedance estimation in a rhythmic task[END_REF] for the ball bouncing task and [START_REF] Fortineau | Interactive robotics for human impedance estimation in a rhythmic task[END_REF] for the interactive robot control and a first approach of impedance estimation based on spline virtual trajectories.

A. Interactive rhythmic task

A user is manipulating an admittance controlled robot, manoeuvring its endpoint as a physical interface, to control a paddle simulated in a virtual environment (see Fig. 1). A 5-DOF KUKA TM youBot arm is used. The motion of the robot arm is restricted to the vertical plane using the control described in [START_REF] Fortineau | Interactive robotics for human impedance estimation in a rhythmic task[END_REF] and is performed using three actuated joints. To yield a rhythmic dynamic behaviour, the user has to move the paddle to bounce a virtual ball subject to gravity and moving along the vertical axis, to a target height simulated in the digital environment. When the ball hits the paddle, force feedback f i is provided by the robot to the human to simulate the impact of the ball on the paddle [START_REF] Kawazoe | Dynamics of the Ball-Racket Impact in Tennis : Contact Force, Contact Time, Coefficient of Restitution, and Deformation[END_REF]. The motion generated can be compared to industrial rhythmic tasks such as stamping, punching or extended to planar movement such as sanding.

To allow impedance estimation during the task, pulse force perturbations are introduced using the robot, in both directions of the uni-dimensional task. Both timing and direction are stochastic. These perturbations are generated at the motor torque level and are consequently filtered by the robot dynamics before reaching the endpoint. In order to have a negligible influence on the task and approach a seamless methodology, the perturbations force effects have a duration under 60 ms and a rather low amplitude, below 5 N. Perturbations are pseudo-randomly introduced at three different phases of the cyclic behaviour. A cycle is defined here as the time frame between two consecutive ball impacts on the paddle. In the considered task, the mean cycle time is of 1.18 s. According to the onset time in the cycle, the three perturbations categories are denoted as shown on Fig. 2 c 1 , c 2 , c 3 . Fig. 2: Time normalised cycles of the data set d 3 ; perturbations are indicated with circles, the dotted lines represent the average cycles while the coloured areas represent 1 and 3 standard deviations from the average cycle. The time normalisation was done using interpolations on the average cycle duration.

B. Data sets

Both the endpoint position and interaction force between the human and robot arms are acquired with external sensors to obtain unbiased data. Interaction force is measured using an ATI TM mini45 force sensor also used in the robot admittance control and the endpoint position is measured using an Optitrack TM V120:Trio motion capture system. For more information, please refer to [START_REF] Fortineau | Interactive robotics for human impedance estimation in a rhythmic task[END_REF].

A single user conducted a total of 28 physical interactions of approximately 5 min each with the experimental test-bed, without random perturbations. Half of those experiments (data set d 1 ) included force feedback for ball impacts (Section II-A). The other half (d 2 ) included no force feedback on impacts.

A third data set, d 3 , was collected by the same user. During a total of 7 manipulations of less than 2 min on average, 165 perturbed trajectories were recorded and classified according to the onset time of the force perturbation in the cycle (Section II-B). Eight of those perturbed trajectories could not be clearly classified because they occurred in unsteady cycles and were discarded. The distribution of onset time of the force perturbation is presented in Fig. 2. The three categories c 1 , c 2 and c 3 are not exactly equally represented since their generation was pseudo-random, they have respectively 60 (+28/-32), 51 (+25/-26) and 46 (+29/-17) elements each. Each of these three categories contains positive and negative perturbations (numbers indicated previously in parentheses).

III. VIRTUAL TRAJECTORIES ESTIMATION

This section depicts the new methodology proposed to estimate unperturbed trajectories, previously called virtual trajectories (VT) in [START_REF] Fortineau | Interactive robotics for human impedance estimation in a rhythmic task[END_REF].

Most impedance estimation methods rely on the fact that deviations δf (force) and δx (position) from VTs (see Section IV for the impedance model relating motion to restitution force) can be introduced either by force or position perturbations. At least one is by essence unknown and needs to be reconstructed, while the other one can be perfectly mastered when the perturbation system allows it. In the literature, the averaged previous unperturbed trajectories were used in [START_REF] Bennett | Time-varying stiffness of human elbow joint during cyclic voluntary movement[END_REF], with the assumption of sufficient reproducibility in the movements. A prediction algorithm with an unbiased error estimation, smaller than 1 cm after 200 ms, was proposed in [START_REF] Burdet | A method for measuring endpoint stiffness during multijoint arm movements[END_REF]. We could not reach better result than those with cubic splines [START_REF] Fortineau | Interactive robotics for human impedance estimation in a rhythmic task[END_REF] while replicating those methods, that were not designed for force trajectories with high variability, as is the case for the benchmark task considered in this study. The variability of position and force profiles is illustrated in Fig. 2.

A. Preliminary remarks for signal analysis

The estimation of the position and force VTs is first validated using the unperturbed data set d 2 (Section II-B). The position VT was well approximated with cubic splines in [START_REF] Fortineau | Interactive robotics for human impedance estimation in a rhythmic task[END_REF]. Accordingly, this part focuses on force trajectories.

1) Temporal analysis: As seen from Fig. 3, the force signals shape varies significantly across cycles, even in the absence of perturbations. While the positions follow roughly a sinusoidal profile (Fig. 2a), forces (Fig. 2b) are affected by higher frequency oscillations inside the cycle. These oscillations should not be considered as perturbations, and their origin can be related to the human-robot physical interaction and the admittance control tuning (while the stability of this control was experimentally verified, strict passivity was not ensured). This situation is not exceptional, all the more if several human subjects operate the same robot and employ different interaction patterns. These intra-cycle oscillations can not be simply filtered out and should be taken into account for an efficient VT estimation. 2) Frequency analysis: In the Fourier transform of the data sets d 1 and d 2 , the main frequency is observed around 0.9 Hz, that is the frequency imposed by the ball bouncing task. While looking at the spectral analysis of both case scenario (with and without force feedback), the main differences occur after 5 Hz in both position and force signals. Those differences are more pronounced in the force spectral density, since the perturbations have more impacts on the force trajectory. Therefore, it can be deduced that most of the perturbed behaviour dynamics modify the trajectories, both in force and position, within a window of 200 ms. It also is worth noticing that the force spectral density is larger, in a broader frequency band than the position, confirming that this signal is less filtered and richer in frequency content.

B. Force virtual trajectory estimation

The estimation of the virtual force using cubic spline interpolations [START_REF] Fortineau | Interactive robotics for human impedance estimation in a rhythmic task[END_REF] raised the issue of sensitivity to the landing point. Indeed, as advanced in Section III-A2, the force signal has higher power intensities, notably in higher frequencies. This makes the application of spline interpolation less reliable than for the position. In [START_REF] Abe | Modulation of elbow joint stiffness in a vertical plane during cyclic movement at lower or higher frequencies than natural frequency[END_REF], the cyclic patterns of previous unperturbed data were fitted to a sinusoidal expression to predict the virtual position. A sum of sines is proposed to approach the virtual force in an offline resolution in a similar spirit. The general form is given in (1), with n the number of sines; A i , f i and φ i the parameters of each sine. a and b the linear components, and l a boolean. Hereinafter, this general form is adapted to the considered case through optimisations, and statistical analysis is performed.

f z (t) = l(at + b) + n i=1 A i sin(2πf i + φ i ) (1)
1) Optimisation: From (1), the tested subset of functions corresponds to 2 ≤ n ≤ 5, which amounts to a total of 8 possible combinations with l = 0 or 1. To identify the parameters of these functions, an optimisation problem is defined to minimise the nonlinear quadratic error, min x ||f (x = {A i , f i , φ i , a, b})|| 2 2 , with a Levenberg-Marquardt algorithm. Thereagainst [START_REF] Burdet | A method for measuring endpoint stiffness during multijoint arm movements[END_REF], [START_REF] Abe | Modulation of elbow joint stiffness in a vertical plane during cyclic movement at lower or higher frequencies than natural frequency[END_REF], the previous unperturbed cycle is not used to solve the VT, but only data within a cycle shortly before and after the perturbation; this is only possible since the estimation is proceeded offline (see Fig. 3). A chunk of the signal is selected with a masked window, corresponding to the perturbed force. The masked interval is chosen to be 100 ms long and will not be used to minimise the quadratic errors. 40 ms before the perturbation and 110 ms after the mask are used for the optimisation definition, as shown in Fig. 3.

This optimal problem is sensitive to its initialisation. Therefore, two options are put forward using d 2 : i) for each of the 14 sessions, the optimisation for the first cycle in the session is performed with 250 starting points and the parameters obtained from this optimisation are then used as initial values for the next cycles, ii) 50 starting points are used for every cycle. 2) Evaluation: In order to evaluate the proposed method, the estimated VTs need to be compared to a ground truth, unavailable in perturbed measurements. To that purpose, the sessions d 2 conducted without perturbations nor force feedback described in Section III-A are used. A total of 3144 evaluations are performed on the rhythmic force trajectories. If the optimised parameters are capable of approximating well the 100 ms masked t mask as well as the next 100 ms t next , then we consider that this methodology is capable of estimating a sufficient portion of the virtual force trajectory, for the impedance estimation, as explained in Section IV.

A i (N) f i (Hz) φ i (rad s -1 ) a (N s -1 ) b (N) b l 0 0 -pi/2 -100 -100 x 0 1 0.9(2i -
The comparison between the different models is then done using the Akaike criterion (AIC) to favour the accuracy of the model and penalise the number of parameters and, therefore, potential overfitting. The model with three sines (n = 3), no linear components (l = 0) and systematic multi-starts (option ii) reached the lowest AIC value and was therefore retained. This model has an average adjusted coefficient of determination (R 2 adj ) of 91.0 % ± 10.0 % on t mask + t next interval and 95.1 % ± 5.3 % on t mask interval, after removing outliers outside five scaled median absolute deviations (sMAD). The distribution of the errors is Gaussianlike centred close to zero, and most of the errors occurred in the interval t mask .

IV. IMPEDANCE IDENTIFICATION METHODOLOGY

This section describes the methodology proposed for the identification of the impedance parameters.

This paper does not intend to differentiate the contribution of reflexes from the intrinsic properties of the human arm, like it was done in [START_REF] Guarin | Unbiased Estimation of Human Joint Intrinsic Mechanical Properties During Movement[END_REF], but rather observe the arm's endpoint behaviour along the uni-dimensional task. A simple viscoelastic 2 nd order model of impedance is used (2), in a comparable manner as [START_REF] Erden | Hand Impedance Measurements During Interactive Manual Welding With a Robot[END_REF], with M , B and K ∈ R nc×nc , respectively the matrices of mass, damping and stiffness in the n c Cartesians dimensions. For the studied uni-dimensional ball bouncing task, only the components along the vertical z axis are considered, n c = 1.

δf = M δ ẍ + Bδ ẋ + Kδx (2)
A. Auto Regressive model 1) Delay: Once the VTs have been computed from experimental data (see Section III), the differential estimated behaviour δf z and δx z can be exposed. Fig. 4 shows isolated perturbed behaviour of those trajectories. As expected from Section II-A, the main dynamics of the force perturbation happens in the first 60 ms, with a maximum reached around 33 ms. There is an apparent latency δ t to notice in the positioning behaviour compared to the force differential. The figure does not allow an accurate evaluation of that delay; it can solely be bounded safely between 5 and 15 ms. This delay is studied in Section V-A. A possible origin would be that the motion capture system, used for the position acquisition, has a sampling rate of 8.3 ms, while position measurements were next interpolated at 1 kHz. A delay of about 8 ms can therefore be expected. 2) Identification model: The following discrete-time difference equation is considered for the identification of the second-order system (2) sampled with a zero-order hold and taking into account the previously mentioned pure delay:

(1 + a 1 q -1 + a 0 q -2 )δx z [T k ] = (b 1 + b 0 q -1 )δf z (T k -δ t ) + e(T k )
(3) with q the delay operator, a i and b j the model coefficients related to K, B and M , δ t a pure delay and e white noise. The time window used for the identification is discussed in Section V-B1.

The eigenvalues of the discrete-time model ( 2) z 1,2 are related to the ones of the model ( 3)

z 1,2 = λ 1,2 via e λ1,2T = exp(0.5T (-B ± √ B 2 + 4M K)/M ) = 0.5(-a 1 ± a 2
1 -4a 0 ), with T the sampling time. The stiffness parameter of (2) corresponds to the static gain and thus from [START_REF] Mussa-Ivaldi | Neural, mechanical, and geometric factors subserving arm posture in humans[END_REF]

, K = (1 + a 1 + a 0 )/(b 1 + b 0 ).
Identifying the impedance parameters K, B and M using this ARX model appears to be more numerically robust than with a straight temporal least square identification technique based on velocities and accelerations. Indeed, this way, the kinematic derivatives are not required, either using sensors or using numerical derivation.

B. Model evaluation 1) Simulation:

A simulation was used to validate both the presented VT estimation and impedance parameters identification methods. In an over-simplistic spirit, the model for the arm movements (trajectory control) was assumed to be the same as the one for the impedance behaviour (disturbance rejection control). For the sake of veracity, empirical force data was injected into the arm model. Force perturbations of 5 N amplitude were generated with a second-order filter of 70 ms pseudo-period to resemble the experimental ones. Those perturbations were injected in an impedance model with known parameters and summed to the simulated arm movement. For simulated impedance behaviours given in Table II, conducted on 104 estimations, an average R 2 adj of 99 % ± 1 % was reached to reconstruct the position with the identified parameters after removing outliers (5 sMAD). The stiffness was the parameter with the highest volatility, therefore to abide by the increase in abnormal data, the threshold was reduced from five to the usual three sMAD.

It can be noticed that with perfect knowledge of the VT, the parameters were estimated with relative errors of about 0.05 %, which means that the estimation errors mainly originate from the errors of VT estimation.

K1 B1 M1 K2 B2 M2 mean 258 N m -1 12.5 N s m -1 523 g 378 N m -1
15.9 N s m -1 487 g std (1) 116 N m -1 6.81 N s m -1 110 g 166 N m -1 5.89 N s m -1 63 g re (2) 29 % 25 % 5 % 5.4 % 6.0 % 2.5 % exact (3) 200 N m -1 10.0 N s m -1 500 g 400 N m -1 15.0 N s m -1 500 g (1) standard deviation, (2) relative mean error, (3) actual simulated parameter values TABLE II: Parameters estimation in simulation.

V. RESULTS

This section presents the results obtained using the experimental protocol previously described in [START_REF] Fortineau | Interactive robotics for human impedance estimation in a rhythmic task[END_REF], with force feedback and controlled timing of force perturbations. The data set d 3 , described in Section II-B, is used. Variance analysis (ANOVA) is conducted with a threshold fixed at 5 % for significance levels. Both the degrees of freedom and the F value will be indicated in parenthesis (DOF, F) for the ANOVAs.

A. Delay

The delay δ t unveiled in Section IV-A1 was determined by running several estimations with varying delay. To evaluate the quality of estimated parameters, the measured position was compared to its reconstruction using the identified parameters and the estimated input force. A one-way ANOVA showed a significant main effect of the delay on R 2 adj (10, 14.75). The highest adjusted coefficient of determination reached a score of 95.1 % ± 5.5 % for 12 ms (5 sMAD); this delay is therefore used for the impedance estimations. Although it is important to stress that a posthoc analysis (Tukey's HSD) revealed that no significant differences were observed between the 12 ms delay and delays comprised between 8 ms to 14 ms.

B. Impedance parameters

1) Time window: As mentioned in [START_REF] Erden | Hand Impedance Measurements During Interactive Manual Welding With a Robot[END_REF], voluntary reactions need to be discarded because linear impedance models cannot correctly model them. Voluntary reactions can be observed as soon as 100 ms after a perturbation [START_REF] Libet | Unconscious cerebral initiative and the role of conscious will in voluntary action[END_REF]. Yet, at a 1 kHz rate, [START_REF] Erden | Hand Impedance Measurements During Interactive Manual Welding With a Robot[END_REF] found that using only 100 ms of data after the perturbation did not bring satisfying results. Therefore, they assumed that voluntary actions were not dominant in the first 200 ms and used that time window for the evaluation. To confirm that observation, impedance estimations were conducted with a data range varying from 50 ms to 200 ms. Fig. 5 shows the average estimated stiffness, damping and mass values as functions of the estimation time window length. Thoses average values no longer vary starting from approximately 150 ms estimation time windows. However, the standard deviation of each of the parameter seems to reach a steady-state before, for about 80 ms, without significant reductions after that milestone. At the final stage, for 200 ms, the stiffness still has a standard deviation that represents 44.8 % of its average value, the damping 34.2 % and the mass 42.9 %. This variability might be explained by the nondifferentiation of the different perturbation categories. This result tends to show that an identification on a time window of 150 ms at 1 kHz might be enough. The paper will study both the identified parameters using 150 ms (subscript i) and 200 ms (subscript ii) to carry this investigation. 2) Outliers: The impedance estimations were separated into several groups according to the cyclic phase and the perturbation direction to isolate possible sources of variability. For both identifications (i, ii), most of the outliers happened in a single class c 2 + (see Fig. 6). Among the 157 classified recordings, a total of 14 values were not considered (8.9 %). Similar results were found for the identification i, but with a total of 17 outliers, which represent 10.8 % of the data.

Mass (kg)

Fig. 6: The bar stacked chart shows in blue recordings kept for further analysis, in orange and purple respectively dismissed recordings with R 2 adj < 50 % (type 1) and stiffnesses outside 5 sMAD (type 2). The yellow indicates the outliers that belong to both to type 1 and 2. The others charts indicate the average and standard deviation for all impedance parameters estimated using the valid recordings and 200 ms (ii).

3) Estimated parameters: In [START_REF] Tsuji | Bio-mimetic impedance control of robotic manipulator for dynamic contact tasks[END_REF], [START_REF] Mussa-Ivaldi | Neural, mechanical, and geometric factors subserving arm posture in humans[END_REF], [START_REF] Dolan | Dynamic and loaded impedance components in the maintenance of human arm posture[END_REF], [START_REF] Erden | Hand Impedance Measurements During Interactive Manual Welding With a Robot[END_REF], Cartesian stiffnesses were found within a magnitude ranging in 40 N m -1 to 700 N m -1 , under static and moving configurations, with and without load for the human arm. Moreover [START_REF] Tsuji | Bio-mimetic impedance control of robotic manipulator for dynamic contact tasks[END_REF], [START_REF] Dolan | Dynamic and loaded impedance components in the maintenance of human arm posture[END_REF], [START_REF] Erden | Hand Impedance Measurements During Interactive Manual Welding With a Robot[END_REF] found damping parameters ranging between 8 N s m -1 to 44 N s m -1 and masses between 0.6 kg to 2.5 kg. The results presented here (Table III) are within those ranges but for the estimated masses that are slightly below. 4) Variance analysis: Considering the results from the identifications using 150 ms, an analysis of variance showed a significant effect of the phase on the adjusted coefficients of determination (5, 15.12). However, a posthoc analysis revealed that differences between positive and negative perturbations were only significant within c 2 . Only one significant difference could be observed between the classes c 1 and c 3 , for the specific case of c - 1 and positive c + 3 . Taking into account the impedance parameters, neither the mean variations of stiffness nor damping observed during the different phases were significant (5, 1.59), (5, 2.37).

c + 1 c - 1 c + 2 c - 2 c + 3 c - 3 + - i K 338
Similar results were found for the identifications done with 200 ms. The only noticeable discrepancy was about the damping values identified, which were found significantly different (5, 4.34), solely because of group c + 2 . Finally, comparing i and ii identifications, in first hand, if regrouping all classes, significant differences could only be found for R 2 adj (1, 12.03) with an increase of 4 % in average for the second identification and non-meaningful mean variations of the impedance parameters between 1 % and 5 %. On the other hand, if comparing separated classes in both identification, (eg. i c + 1 with ii c + 1 ), no significant differences were found.

VI. CONCLUSIONS

The identification method proposed here with the experimental test-bed design described in [START_REF] Fortineau | Interactive robotics for human impedance estimation in a rhythmic task[END_REF], offers a lowconstraint methodology for the impedance identification in the absence of knowledge on both force and position VTs during a task. Estimated endpoint impedance parameters obtained with the proposed method are within the range found in the literature. However, the identification being sensitive to errors in the VTs as well as non-linear phenomena (voluntary movements, non-linear reflexes), the method still requires numerous estimations for reliable results. The results in this paper tend to show that using a time window of only 150 ms for the identification was sufficient to bring result non significantly different from those with a 200 ms time window.

No significant stiffness variations were observed in the three cyclic phases chosen. Most of the significant variations were concentrated in class c + 2 . Compared to previous work, using an ARX methodology improved the repeatability in the estimation of both the damping and mass. Using the new force trajectory estimation method drastically reduced the errors in VTs and improved the estimation of the stiffness, for perturbation below 5 N inducing typical movements deviations δx less than 1 cm. It made the perturbations hardly distinguishable from the task force feedback without disrupting the task.

Further experiments will be conducted to evaluate the inter and intra-individual variability of end-point impedance parameters in interactive tasks. Obtained parameter values can be used for the analysis and design of feedback control strategies in interactive robotics.

Fig. 1 :

 1 Fig. 1: Experimental set-up

Fig. 3 :

 3 Fig. 3: Two examples of force trajectory estimations, with the measured unperturbed force (solid blue), the simulated perturbations (dashed orange) and estimated unperturbed force (green).Figure (a) depicts a successful estimation, while figure (b) shows a failed estimation that leads to an underestimated restitution force with oscillations.

  Fig. 3: Two examples of force trajectory estimations, with the measured unperturbed force (solid blue), the simulated perturbations (dashed orange) and estimated unperturbed force (green).Figure (a) depicts a successful estimation, while figure (b) shows a failed estimation that leads to an underestimated restitution force with oscillations.

  Lower (b l ) and upper (bu) boundaries for the random initial starting points, and first manually fed initial condition (x 0 )

Fig. 4 :

 4 Fig.4: Isolated typical effects of the perturbations on both the position and force. The average behaviours are shown with dotted lines, while the coloured areas represent the standard deviation. Positive (green) and negative (orange) perturbations are differentiated by their colour.

Fig. 5 :

 5 Fig. 5: Influence of the number of samples used for the ARX identification, the blue dotted line represents the average, and the light blue area represents one standard deviation from the average. The orange dotted lines show the 5 % boundaries of the final values.

TABLE III :

 III Results: standard deviation are indicated in as confidence intervals.

		± 158	319 ± 151	408 ± 227	294 ± 173	308 ± 133	397 ± 221	345 ± 172	343 ± 193
	i B	16.7 ± 5.2	18.5 ± 4.2	17.8 ± 12.0	13.7 ± 4.4	15.6 ± 6.6	13.2 ± 5.4	16.6 ± 7.8	14.6 ± 5.2
	i M	0.288 ± 0.061	0.369 ± 0.235	0.756 ± 0.788	0.293 ± 0.071	0.354 ± 0.123	0.484 ± 0.393	0.430 ± 0.441	0.391 ± 0.291
	i R 2 adj	83.5 ± 4.6	88.7 ± 3.9	73.6 ± 14.0	93.0 ± 2.9	81.6 ± 9.8	82.6 ± 6.9	80.3 ± 10.2	87.6 ± 6.9
	ii K	340 ± 147	310 ± 164	387 ± 177	298 ± 159	368 ± 167	411 ± 241	362 ± 160	348 ± 203
	ii B	16.3 ± 4.5	19.1 ± 4.5	20.7 ± 7.9	14.2 ± 4.7	16.2 ± 4.8	14.5 ± 6.7	17.4 ± 5.9	15.4 ± 5.9
	ii M	0.294 ± 0.057	0.378 ± 0.259	0.664 ± 0.581	0.305 ± 0.088	0.381 ± 0.147	0.471 ± 0.253	0.418 ± 0.334	0.392 ± 0.222
	ii R 2 adj	88.1 ± 3.6	91.6 ± 3.0	80.5 ± 11.4	94.0 ± 7.5	86.4 ± 7.0	84.7 ± 6.8	80.4 ± 10.7	85.6 ± 10.4

CIAMS, Université Paris-Saclay, 91405, Orsay, France, & Université d'Orléans, 45067, Orléans, France