Félix Rutard

Olivier Sigaud

Mohamed Chetouani

TIRL: Enriching Actor-Critic RL with non-expert human teachers and a Trust Model

Reinforcement learning (RL) algorithms have been demonstrated to be very attractive tools to train agents to achieve sequential tasks. However, these algorithms require too many training data to converge to be efficiently applied to physical robots. By using a human teacher, the learning process can be made faster and more robust, but the overall performance heavily depends on the quality and availability of teacher demonstrations or instructions. In particular, when these teaching signals are inadequate, the agent may fail to learn an optimal policy. In this paper, we introduce a trustbased interactive task learning approach. We propose an RL architecture able to learn both from environment rewards and from various sparse teaching signals provided by non-expert teachers, using an actor-critic agent, a human model and a trust model. We evaluate the performance of this architecture on 4 different setups using a maze environment with different simulated teachers and show that the benefits of the trust model.

I. INTRODUCTION

Various studies have shown that Reinforcement Learning (RL) can be used to train a robot to perform various tasks in autonomy [START_REF] Kalashnikov | Qt-Opt: Scalable deep reinforcement learning for vision-based robotic manipulation[END_REF], [START_REF] Andrychowicz | Learning dexterous in-hand manipulation[END_REF], [START_REF] Akkaya | Solving rubik's cube with a robot hand[END_REF], [START_REF] Kober | Reinforcement learning in robotics: A survey[END_REF] in labs, but the corresponding algorithms generally suffer from a very low sample efficiency which render them impractical for deployment in real-world applications. This limitation can be addressed by a human partner in charge of accelerating the learning process by providing adequate demonstrations, instructions or feedback. In such situations where humans and robots share the same social and task spaces, identifying factors for efficient, safe and successful interactions is challenging. Trust is one of such factors and has mainly be considered as the human willingness to cooperate with the robot [START_REF] Lee | Computationally modeling interpersonal trust[END_REF]. It is usually considered as a critical factor for establishing and maintaining effective relationships with robots [START_REF] Hancock | A meta-analysis of factors affecting trust in human-robot interaction[END_REF]. In [START_REF] Gaudiello | Trust as indicator of robot functional and social acceptance. an experimental study on user conformation to icub answers[END_REF], trust in the robot has been considered as a main indicator of acceptance and measured by the participants' conformation to the robot's answers to questions on functional and social tasks.

The main focus of such works was the human willingness to cooperate and interact with a robot. In this paper, we propose to reverse the question by considering situations where human and robot are engaged in an interactive task learning scenario, and in case of non-optimal teaching, we investigate mechanisms allowing a robot to efficiently learn the task by considering a human task trust assessment. We propose to estimate this trust by assessing how the human teacher guides the robot during interactive task learning.

Sorbonne Université, CNRS, Institut des Systèmes Intelligents et de Robotique, ISIR, F-75005 Paris, France II. RELATED WORK As described in [START_REF] Broekens | Towards transparent robot learning through tdrl-based emotional expressions[END_REF], Interactive Robot Learning deals with models and methodologies allowing a human teacher to guide the learning process of the robot by providing it teaching signals [START_REF] Thomaz | Teachable characters: User studies, design principles, and learning performance[END_REF]. Usual teaching signals include instructions [START_REF] Grizou | Robot learning simultaneously a task and how to interpret human instructions[END_REF], [START_REF] Paléologue | Semantic-based interaction for teaching robot behavior compositions[END_REF] advice [START_REF] Griffith | Policy shaping: Integrating human feedback with reinforcement learning[END_REF], demonstrations [START_REF] Argall | A survey of robot learning from demonstration[END_REF], guidance [START_REF] Suay | Effect of human guidance and state space size on interactive reinforcement learning[END_REF], [START_REF] Najar | Training a robot with evaluative feedback and unlabeled guidance signals[END_REF] and evaluative feedback [START_REF] Knox | Training a Robot via Human Feedback: A Case Study[END_REF], [START_REF] Najar | Training a robot with evaluative feedback and unlabeled guidance signals[END_REF]. A standard assumption in Interactive Robot Learning is that teaching signals are provided by experts. A distinguishing feature of our work is that our system can handle erroneous teaching signals. In that respect, the most closely related works are [START_REF] Najar | Training a robot with evaluative feedback and unlabeled guidance signals[END_REF], [START_REF] Celemin | A fast hybrid reinforcement learning framework with human corrective feedback[END_REF], [START_REF] Zhu | Learning sparse rewarded tasks from sub-optimal demonstrations[END_REF] and more recently [START_REF] Kwon | When humans aren't optimal: Robots that collaborate with risk-aware humans[END_REF].

To deal with this more difficult concern, our key contribution in this paper is to endow the agent with a trust model towards its teachers.

A taxonomy of trust related errors and mitigation strategies is proposed in [START_REF] Tolmeijer | Taxonomy of trust-relevant failures and mitigation strategies[END_REF]. In the literature, most of the works about trust in human-robot interaction focus on the human towards robot trust point of view, in sharp contrast with our approach. For instance, in [START_REF] Salem | Would you trust a (faulty) robot? effects of error, task type and personality on human-robot cooperation and trust[END_REF], the authors investigated the effects of robot's error, task type and personality on cooperation and trust. Similarly, in [START_REF] Robinette | Effect of robot performance on human-robot trust in time-critical situations[END_REF], the authors also investigate the effects of robot error on human trust with a specific focus on time-critical situations. Few different works studied trust from a different perspective, such as [START_REF] Surendran | Your robot is watching: Using surface cues to evaluate the trustworthiness of human actions[END_REF] in which the robot uses surface cues to estimate trustworthiness of the human it interacts with. Other approaches develop methods in which enable the agent to use its own confidence to decide when to ask a human to help it by providing teaching signals [START_REF] Chernova | Interactive policy learning through confidence-based autonomy[END_REF].

In our study, the definition of the trust of the agent towards the teacher is closely correlated to the expertise level estimation of the teacher by the agent. A similar work studied the estimation of user's expertise level by a robot [START_REF] Carreno-Medrano | Incremental estimation of users' expertise level[END_REF]. The authors proposed a method to estimate the user's expertise level relying on a Bayesian inference framework updating incrementally the human level of expertise model processed using a Boltzmann policy. Our approach requires fewer computation and can be incrementally performed within the interaction loop.

Finally, like us the authors of [START_REF] Nguyen | Active choice of teachers, learning strategies and goals for a socially guided intrinsic motivation learner[END_REF] propose a mechanism where the agent can choose which human to listen to when they are not equally trustworthy.

III. MODEL

In this section we describe TIRL (Trust-based Interactive Reinforcement Learning), our approach to endow a reinforcement learning agent with the capability to make profit of non expert teaching signals. The architecture depicted in Figure 1 is composed of two blocks. First, an actor-critic agent learns from rewards a stochastic policy represented as a distribution. Second, the agent learns one or multiple models of non-expert teacher(s) who provide sparse teaching signals as instructions and evaluative feedback. In this paper, we address the problem of selecting an action using either the suggestion of the agent's actor or the teacher guidance signal. For this purpose, we introduce a trust value estimation in human guidance.

A. Overview

Agent

A key design choice in our architecture is that evaluative feedback is always used to update the trust model through the teacher state-action value, whereas instructions are only followed if the trust of the agent towards the teacher is high enough. The trust model estimates a value of trust of the agent towards the teacher based on the difference between the actor's probability distribution over actions and the probability distribution over actions computed by the teacher state-action value. The teacher state-action value estimates how good the action is in the state according to the teacher, based on the provided evaluative feedback. A specific value propagation mechanism is necessary here to deal with sparse teaching signals.

In this architecture we used a fixed size experience replay buffer as described in [START_REF] Lin | Self-improving reactive agents based on reinforcement learning, planning and teaching[END_REF] and a mini-batch training mode for the update of all architecture estimators in order to improve training efficiency. The indice j used in the following pseudo-codes describing the architecture refers to the j th transition of the mini-batch transitions sampled randomly from the replay buffer.

See Algorithm 1 for the description of the complete algorithm.

B. Reinforcement Learning agent

We consider a Markov Decision Process (MDP) (ρ, S, A, T , R, γ) where ρ is the initial state distribution, S the state space, A the action space, T : S × A → S the transition function, R : S → R the reward function and γ ≤ 1 is a discount factor. The RL agent uses an Actor-Critic algorithm as described in [START_REF] Sutton | Reinforcement Learning: An Introduction[END_REF]. The critic is defined as an estimated state value function V (s t ) which estimates

s prev ← s, a prev ← a, f prev ← f , s ← s
Sample mini-batch of size N in replay buffer Update critic V and actor π from mini-batch for each teacher do Update teacher state-action value estimator H from mini-batch Update trust estimator T from mini-batch end for end for the true value function V (s t ) = E πt (J t |S t = s) defined as the expected long-term reward return while following policy π t . The long-term reward return is defined by J t = End k=0 γ agent k r t+k+1 with k = End corresponding to the end of episode and r t+k+1 corresponding to the reward from the environment received at step t + k + 1.

The critic is updated using:

V (s t ) ← V (s t ) + α critic δ agentt where δ agentt is the tempo- ral difference (TD) error δ agentt = r t +γ agent V (s t )-V (s t ).
The actor π(s, a)

t = P r(a t = a|s t = s) = e p(s,a) τ k e p(s,k) τ
is defined as a probability distribution over action obtained from a softmax function with temperature τ on state-action preferences p(s, a). The temperature parameter τ controls the exploration behavior of the agent: with low values the agent is more greedy and with high values of τ it explores a lot. State-action preferences p(s, a) are initialized at 0 and are updated using the previously processed TD error δ agentt as p(s t , a t ) ← p(s t , a t ) + α actor δ agentt with α actor the update parameter of the actor. We call π * the optimal policy defined such that V π * (s) = max π V π (s) | ∀s and we define the optimal path as the set of states crossed by the optimal policy starting from the initial state. See Algorithm 2 for the complete actorcritic update.

Algorithm 2 Actor π and critic V update Inputs: mini-batch of transitions of size N , actor π, critic V , actor and critic learning rates α actor and α critic , agent discount factor γ agent for each transition j in mini-batch do Compute TD error with critic V :

δ agentj = r j + γ agent V (s j ) -V (s j ) Update Critic V (s j ) ← V (s j ) + α critic δ agentj Update p(s j , a j ) ← p(s j , a j ) + α actor δ agentj Process Actor π(s, •) = e p(s,•) τ k e p(s,k) τ end for

C. Non-expert human teacher

We consider a non-expert human teacher who can provide two types of teaching signals: evaluative feedback and instructions. Evaluative feedback is a scalar provided after one or multiple agent transitions, which can take values in {-1, 0, 1}. An instruction is a human guidance signal toward an action associated with the current agent's state. These signals come with a sparsity rate: a sparsity of 0 corresponds to a silent teaching signal and a sparsity of 1 corresponds to the teacher providing it at each time step.

An evaluative feedback is erroneous if it is positive when the agent is not following the optimal policy and negative when the agent behaves optimally, whereas an instruction is erroneous if it provides actions different from the optimal policy. A human teacher who never provides erroneous teaching signals is said expert.

The agent is not aware of the sparsity rates and the details of the erroneous zone (if any) of the teacher.

D. Teacher state-action value estimator

A teacher can provide evaluative feedback for a given transition or for a set of transitions and the agent does not know which transition(s) the evaluative feedback is about. This problem is known as the temporal credit assignment problem and has already been tackled by the RL community using temporal difference (TD) learning.

The teacher state-action value estimates the expected value corresponding to the evaluative feedback return starting from state s t = s and playing action a t = a while following the agent's policy π t . To perform this estimation, we define a specific discount factor γ teacher ∈ [0, 1] which accounts for how much the agent should back-propagate the evaluative feedback signal in time. Put differently, the parameter γ teacher controls the credit decay time horizon of evaluative feedback for teacher state-action value, thus we need to set this γ teacher to fit the assumed credit decay time horizon when the teacher provides an evaluative feedback.

The teacher state-action value Ĥ(s, a) is defined as the estimate of the expected evaluative feedback return starting from state s and playing action a: Ĥ(s, a) = E(F t |S t = s, A t = a). It uses the teacher evaluative feedback return

F t = End k=0
γ teacher k f t+k+1 with f t+k+1 corresponding to the evaluative feedback provided by the teacher at step t + k + 1. Ĥ(s, a) is updated using Bellman equation with the SAFSA algorithm. The SAFSA algorithm comes from the SARSA algorithm with evaluative feedback (F) replacing the reward (R). The SARSA algorithm commonly used in RL is defined as

Q(s, a) ← Q(s, a) + α (r t + γ Q(s , a ) -Q(s, a))
with Q corresponding to the state-action value, s corresponding to the current state, a corresponding to the current action, s corresponding to the future state, a corresponding to the future action, γ corresponding to the discount factor and α corresponding to a learning rate.

We modified the SARSA algorithm using Ĥ(s, a) instead of Q(s, a) and the teacher evaluative feedback f t instead of environment reward r t and obtained Ĥ(s, a) ← Ĥ(s, a) + α teacher (f t + γ teacher Ĥ(s , a ) -Ĥ(s, a)).

The SARSA algorithm tends to update the state-action value relatively to the current policy while the update rule in Q-Learning tends to update the state-action value relatively to the greedy policy. Since Ĥ(s, a) represents the teacher value according to its own policy, we chose to use the SARSA update rule rather than the Q-learning update rule.

In our implementation, instead of using the current state, action and teacher evaluative feedback (s, a, f t ) and future state and action (s , a ), we used previous state, action and feedback (s prev , a prev , f prev ) and current state and action (s, a). This is due to SAFSA algorithm which needs two transitions in order to obtain every component of the update rule (s prev , a prev , f prev , s, a). See Algorithm 3 for the update of Ĥ.

Algorithm 3 Teacher state-action value estimator Ĥ update Inputs: mini-batch of size N, actor π, critic V , Ĥ update coefficient α teacher For each transition of the mini-batch update Ĥ: for steps j in {1, 2, ...N } do Compute TD error with SAFSA: δ Hj = f prevj + γ teacher Ĥ(s j , a j ) -Ĥ(s prevj , a prevj ) Update Ĥ(s prevj , a prevj ): Ĥ(s prevj , a prevj ) ← Ĥ(s prevj , a prevj ) + α teacher δ Hj end for As a result, the replay buffer stores any previous evaluative feedback teaching signal provided by the teacher f prev as well as previous (s prev , a prev , s) and current (s, a, s , r) agent transitions in order to update both actor-critic and the teacher state-action value estimator.

E. Trust estimator

With a trust model, the agent can filter teaching signals by evaluating the teacher considering the agent's state. By listening to teacher's instructions when the trust model is high and ignoring teacher's instructions when the trust is low, we leverage teacher's knowledge to train faster while we avoid too much wastes of time or dead-ends when the teacher is wrong. Our trust value is obtained by comparing two discrete distributions over the action space: the stochastic actor's distribution π(s, •) and teacher state-action value distribution H dist (s, •), with

H dist (s, •) = e Ĥ(s,•) k e Ĥ(s,k) (1)
for a given state s. This comparison is performed using a Bhattacharyya distance D B between π(s, •) and H dist (s, •):

D B (H dist (s, •), π(s, •)).
We want the trust value to be 1 if both distributions are identical and 0 if they are completely different. Since D B (H dist (s, •), π(s, •)) is defined on [0, +∞[ with 0 in the case of identical distributions, we compute the trust model using a exponentiated negative Bhattacharyya distance:

T (s) = e -D B (H dist (s,•),π(s,•)) = k H dist (s, k) π(s, k).
We want to update this trust model only if H(s, •) and π(s, •) are accurate. Otherwise, we may process trust based on the Bhattacharyya distance between two random distributions. Consequently, the trust model T (s) is initialized at 1 (agent starts fully confident in the teacher) and updated if both the agent and teacher TD errors associated to the state are below thresholds δ agent thr and δ teacher thr . See Algorithm 4 for details.

Algorithm 4 Trust estimator T update

Inputs: mini-batch of size N For each transition of the mini-batch update T : for steps j in {1, 2, ...N } do Compute H dist (s j , •) using (1) if δ agent (s j ) ≤ δ agent thr and δ teacher (s j ) ≤ δ teacher thr then

T (s j ) = k H dist (s j , k) π(s j , k)
end if end for

F. Action selection

During training, we want the agent to filter teacher's instructions depending on the trust value. At each time step, if an instruction was provided by the teacher and if trust is above a threshold T thr , the agent follows the teacher's instruction. Otherwise, the agent plays the action processed by its actor.

This strategy could also be employed in case of multiple teachers, where we need to choose which teacher to learn from. If only one teacher provides an instruction and the trust T (s) associated to the current state for this teacher is greater than the trust threshold T thr , we select this instruction. Otherwise, if multiple teachers provide instructions, we apply the following heuristic: we select a teacher based on the trust attributed to each teacher who provided instructions and whose trust value associated to the current state is greater than the trust threshold (T (s) > T thr ). We call this subset of teachers validated teachers V T . Hence we apply a softmax function to the trust vector T (s) of each validated teacher and we sample a teacher from this probability vector K, the selected instruction is the instruction provided by the selected teacher. See Algorithm 5 for multiple teacher instruction selection.

Algorithm 5 Multiple teacher instruction selection (MTIS)

Inputs: teacher vector T each, instructions vector I[teacher], trust vector T , Filter valid teachers:

V T = T each(T > T thr ); if V T = ∅ then return valid = F alse ; instruction = N one else
Apply softmax to trust:

K = sof tmax(V ) Sample selectedT eacher in V T with probability K return valid = T rue ; instruction = I[selectedT eacher] end if
To control the balance between agent exploration and following teacher's instructions, we use a stochastic instruction following policy parameterized by if r, the instruction following ratio. At each time step, if an instruction was provided by the teacher, we sample a value from a uniform distribution U(0, 1) and if this value is below the instruction following ratio and the trust value is greater than the trust threshold, the agent plays the instruction. Otherwise, the agent plays the action suggested by its actor. See Algorithm 6 for action selection. In our experiments, we bind to the teacher an erroneous zone in which he or she provides erroneous teaching signals.

In the maze environment the erroneous zone is defined by a geometry, which is either a square (associated with its position in the maze) or a diagonal zone, and a size, which is the edge size for the square zone, and a strip width for the diagonal zone.

In order to test our architecture we study the training performance (time to complete episode versus the episode index) in four different experimental setups:

• Exp 1: Agent trained with an expert teacher (no erroneous zone) with different sparsity rates from 0.0 (no teaching signals provided by the teacher) to 1.0 (teaching signals provided at each step);

• Exp 2: Agent trained with an expert teacher (no erroneous zone) then trained with a non-expert teacher (square erroneous zone of edge size 10 and including some optimal path's states) with different sparsity rates (from 0 to 1, with the same sparsity rates for both teachers);

• Exp 3: Agent trained with a non-expert teacher (square erroneous zone of edge size 10 and including some optimal path's states, figure 3) with different sparsity rates (from 0 to 1);

• Exp 4: Agent trained with two non-expert teachers sharing the same erroneous zone shape and size but not the same position (square erroneous zone of edge size 10 and including some states from the optimal path, figure 3) with different sparsity (from 0 to 1, with the same sparsity rates for both teachers). Fig. 3: Erroneous zones and optimal policy path in magenta

C. Hyper-parameters and settings

The mini-batch size N and the experience replay buffer size R have been arbitrarily set to values 10 and 5000 respectively. We begin by setting the actor-critic parameters τ , α actor and α critic which is done by training the agent without any teaching signals. The environment discount factor γ has been set to 0.99, a standard value in RL. The agent exploration temperature τ , actor and critic learning rates α actor and α critic have been set using a grid search method trying to maximize the performance. We then set teacher related parameters α teacher and γ teacher . To set α teacher , we start with a α teacher and increase it until learning stops converging. We apply this method for different sparsity rates and we keep the lowest acceptable value of α teacher . Setting γ teacher , δ teacher and δ actor is performed in a similar way using again different sparsity values. Finally we set T thr , using a grid search for different configurations of teachers. The settings used in this study are described in Table I

V. RESULTS

All experiments have been conducted using 20 different seeds for each set of sparsity rates. Figure 4 shows the time to complete an episode in log scale versus the episode index in the different settings. The lower this time, the better the performance. Results are presented for each sparsity rate as a mean represented by a line and trust interval of 95% represented by a shadow band area around the line. In Exp. 1, the teacher is expert, providing no erroneous teaching signal. Figure 4a shows that every sparsity rate condition converged to a near-identical value. However, the higher the sparsity the faster the training and the more robust the training. Indeed the higher the sparsity rate, the narrower the confidence interval band.

B. Second experiment

In Exp. 2, the teacher is expert for the first 120 episodes, but not for the 120 last ones. Figure 4b shows that for the 120 first episodes, results are identical to the first experiment, since we use the same training conditions. However we identify a performance drop at the beginning of the 120 last episodes of the experiment which vanishes in a few dozen episodes. The higher the sparsity rates, the bigger is the performance drop and the slower it vanishes. This can be explained by the trust value being very high in every state for the first 120 episodes since the first teacher was expert, and the agent taking more time to lose trust after episode 120 when the teaching signal is more sparse. When the signal is not sparse, the agent successfully attributes in a dozen episodes a low trust value to the teacher erroneous zone states and improves its performance using its own policy in the erroneous states which leads to a performance increase.

Figure 5a shows the trust map of the teacher just before swapping between expert and non-expert teachers (step = 119). As expected, the trust is high. Figure 5b shows the trust map just after the swap and we observe the trust in the teacher erroneous zone has not been reduced since the agent still had not enough time to discover the teacher erroneous zone. Figure 5c shows the last episode (step = 239) trust map on which we can identify very low trust attributed to states crossed by the agent's policy which overlap the erroneous zone of the teacher.

C. Third experiment

In Exp. 3, the teacher is not expert. Figure 4c shows that the higher the sparsity rate, the lower the performance at the beginning of the training (expected for sparsity rate 0.0, that is no teaching signal). Besides, the higher the sparsity rate, the higher the convergence speed except for sparsity 1.0 which is a little slower than sparsity rate 0.8. For high sparsity rate, the teacher provides way more erroneous signal than for lower sparsity rates, but the agent faster identifies the teacher's erroneous zone and so to learns faster. However, a high sparsity rate induces wastes of time or dead-ends at the beginning of training, because the agent ignores the teacher's erroneous zone. We identify an optimum at sparsity rate = 0.6 in Figure 4c.

Figure 6 shows the last episode trust map, we notice very low trust has only been attributed to states crossed by the agent's policy which overlap the erroneous zone of the teacher. We also observe the agent policy is not optimal. Indeed, during training, the agent has been misled by the teacher in the teacher erroneous space. Thus the agent spent a lot of time trapped in this zone without being rewarded. This phenomenon led the agent to estimate expected reward in states overlapping teacher's erroneous zone to be very low. As a consequence, the agent's policy avoided the teacher's erroneous zone.

D. Fourth experiment

In Exp. 4, there are two non-expert teachers. Figure 4d shows that the higher the sparsity, the faster and more robust the training (the higher sparsity rate, the narrower the trust band). The trade-off between the low performance at the beginning of training and the high convergence speed we encountered during Exp. 3 is not present anymore. This can be attributed to the complementary nature of the erroneous zones of both teachers. While one teacher is providing erroneous teaching signal, the other one is providing correct teaching signal. The trust mechanism helps the agent to ignore the erroneous teacher depending on its state. This demonstrates how useful the trust mechanism is when using multiple teachers.

Figure 7a shows the last episode trust map of the first teacher. We notice very low trust has only been attributed to states crossed by the agent's policy which overlap the erroneous zone of the first teacher. Figure 7b shows the last episode trust map of the second teacher. We also notice very low trust in states crossed by the agent's policy which overlap the erroneous zone of the second teacher. This shows that the agent correctly estimates and uses the trust mechanism when training with multiple teachers. An immediate extension for future work consists in adding more types of teaching signals such as demonstrations. Besides, in this work, the interpretation of evaluative feedback and instructions was given, but we could incorporate an interpretation learning process by using the methodology of the TICS architecture [START_REF] Najar | Training a robot with evaluative feedback and unlabeled guidance signals[END_REF]. Finally, we implemented and tested TIRL using a discrete environment to facilitate the interpretation of our results focused on trust, but our work could easily be extended to continuous state and action environments using state-of-the-art deep reinforcement learning algorithms.

Fig. 1 :

 1 Fig. 1: Architecture overview: the agent can estimate trust in human teaching signals during interactive task learning

Algorithm 6

 6 Action selectionInputs: Actor π, state s, trust threshold T thr , instructionfollowing ratio if r, instruction or (teachers vector T each, instructions vectorI[T each], trust vector T ) if multiple teachers if multiple teachers then instruction, valid = M T IS(T each, I[T each], T ) else valid = T (s) ≥ T thr end if Sample x from uniform distribution U(0, 1) if instruction is available and valid and x ≤ if r then a ← instruction else a ∼ π(s, •) end if IV. EXPERIMENTS A. EnvironmentWe use a maze with 60 × 60 cells, each action leads to a reward r = -0.1 (mazesize) 2 . The agent receives a reward r = 1 for reaching the terminal state. An episode is considered finished with success if the agent reaches terminal state within (maze size ) 2 and is considered finished with failure if the agent reaches terminal after more than (maze size ) 2 or if the number of steps of the episode exceeds 100 (maze size ) 2 steps. Training is considered finished if the agent has succeeded 100 episodes in a row. The initial state is the top left corner of the maze and the rewarding state is the bottom right corner of the maze.

Fig. 2 :

 2 Fig. 2: Maze with optimal policy path in magenta

  (a) Err. zone in Exp. 2&3 (b) Erroneous zones in Exp. 4

4 Fig. 4 :

 44 Fig. 4: Time to complete episode in log scale versus the episodes index

Fig. 5 :

 5 Fig. 5: Trust maps in Exp. 2. The color scale shows trust towards the teacher along the trajectory.

Fig. 6 :

 6 Fig. 6: Trust maps of last episode in Exp. 3

Fig. 7 :

 7 Fig. 7: Trust maps in Exp. 4

  prev , a prev , f prev , s, a, r, s ) Store current episode for future transition storage:

	Algorithm 1 Trust-based Interactive Reinforcement Learning
	Inputs: actor π(s), critic V (s),
	teacher state-action value estimator H(s, a)
	initialized trust estimator T (s),
	for steps t in {1, 2, 3, ...k} do
	for each teacher do
	if instruction is available then
	i ← instruction
	end if
	end for
	Select action:
	a ← action selection(i, π, if r, T thr )
	Play action a and observe (s , r)
	for each teacher do
	if evaluative f eedback is available then
	f ← evaluative f eedback
	end if
	end for
	Add transition to replay buffer:
	(s

  below.

TABLE I :

 I Hyper-parameters of the study

	Hyper-parameter name	value
	Maze size	60
	Mini-batch size N	10
	Experience replay buffer size R	5000
	Instruction-following ratio if r	0.6
	Trust threshold T thr	0.6
	Environment discount factor γ	0.99
	Teacher discount factor γ teacher	0.75
	Actor learning rate αactor	6
	Critic learning rate α critic	0.6
	Teacher learning rate α teacher Agent TD error threshold δagent thr Teacher TD error threshold δ teacher thr Exploration temperature τ	0.7 10 -4 10 -3 2.5


Benjamin Fasquelle
email: benjamin.fasquelle@ls2n.fr

Cable lengths calculation, modeling and identification of parameters on a manipulator inspired by bird necks

Introduction

Here we will detail the modifications to be made to the cable length calculation when taking into account the pulleys.

The application of forces on the manipulator is linked to the cables and the taking into account of the pulleys therefore modifies it.

Models of friction and elasticity of the cables will then be presented. The associated parameters will be identified in order to obtain a dynamic model allowing a good estimation of the behavior of our prototype.

Manipulator

We first introduce the different notations that will be useful to us later on:

• N: the number of mechanisms composing the manipulator. The mechanisms are numbered from bottom to top: the mechanism at the base of the manipulator is the first mechanism, the one just above the second mechanism, and so on until the highest mechanism which is the mechanism N.

• α: the vector of orientations. The orientation of the mechanism i is noted α i .

• φ and ψ: the vectors of the orientations of the diagonal bars with respect to the base bar of each mechanism. The i elements of these vectors are noted φ i and ψ i .

• r: the radius of the pulleys. It is assumed that all pulleys have the same radius.

• L: length of the diagonal bars (identical for each mechanism)

• b: length of base bars (identical for each mechanism)

• l( j) is the length of the cable associated with the motor j. The length of the cable is decomposed as follows:

l( j) = l 0 ( j) + N ∑ i=1 l i ( j)
, where l i ( j) is the length of the part of the cable j that depends only on the mechanism i (of its geometry and orientation), and l 0 ( j) is a constant corresponding to the rest of the cable.

In this report, we are interested in the manipulator represented in Figure [START_REF] Benjamin | Dynamic modeling and control of a tensegrity manipulator mimicking a bird neck[END_REF]. This one is composed of N = 3 mechanisms and is operated using N m = 4 cables.

Application of forces on the manipulator

The equation of motion of the manipulator is [2]:

M(α) α + c(α, α) + g(α) = Q (1)
where M is the inertia matrix, c the vector of Coriolis and centrifugal forces, g the vector of forces deriving from potential energy (gravity, spring forces) and Q the vector of generalized forces. The virtual work of the force f j applied on a wire cable j is the product of the force and the displacement of the wire cable:

δW = f j δ l( j) (2)
Thus, the vector Q can be written as Q = Z(α)f where Z is a matrix of size N × N m whose elements are:

Z(i, j) = - ∂ l i ( j) ∂ α i (3) α 1 α 2 α 3 c 4 c 3 c 2 c 1

Cable length 4.1 Cable extensions without pulleys

If we consider that there are no pulleys, the cables are arranged:

• either along a strut-routed bar, and the corresponding partial length is therefore that of the bar • either on the side of a mechanism, like springs.

It is therefore necessary to calculate the lengths of the springs to the left and right of the mechanisms. Let d(α i ) be the distance between the middle of the bottom bar and the middle of the top bar of mechanism i. We have:

d(α i ) = L 2 -b 2 cos 2 ( α i 2 ) (4)
Let l l (α i ) and l r (α i ) be the lengths of the springs respectively on the left and on the right.We have:

l l (α i ) = d(α i ) -b sin( α i 2 ) (5) l r (α i ) = d(α i ) + b sin( α i 2 ) (6)
The unwound lengths of the cables of our manipulator are therefore (starting from the bottom of the manipulator):

l c (1) = l r (α 1 ) (7) l c (2) = L + b + l r (α 2 ) (8) l c (3) = 2L + 2b + l r (α 3 ) (9) l c (4) = l l (α 1 ) + l l (α 2 ) + l l (α 3 ) (10)
Thus, the matrix Z can be written in the following form: where Z 1 is a diagonal matrix:

Z = (Z 1 Z 2 ) (11)
Z 1 =     ∂ l r (α 1 ) ∂ α 1 0 0 0 ∂ l r (α 2 ) ∂ α 2 0 0 0 ∂ l r (α 3 ) ∂ α 3     (12)
and Z 2 a vector:

Z 2 =     ∂ l l (α 1 ) ∂ α 1 ∂ l l (α 2 ) ∂ α 2 ∂ l l (α 3 ) ∂ α 3     (13)

Taking into account the pulleys

Decomposition

We will calculate the partial lengths of the cables passing through the mechanism i. The total length will therefore be the sum of these partial lengths.

We will study four different cases:

• The one where the cable passes on the left side of the mechanism i and is attached to the middle of the top pulley.

• The one where the cable passes on the right side of the mechanism i and is attached to the middle of the top pulley.

• The one where the cable passes on the left side of the mechanism and is wound on the top pulley to be redirected above.

• The one where the cable is strut-routed to operate upper mechanisms.

The four cases are summarized in Figure (2). In order to delimit the portions of cables corresponding to each pulley, the lowest point of the pulley that crosses the vertical when the manipulator is in a straight position (α = 0deg for each mechanism) will be taken as a reference point on the pulleys. The portion of cable corresponding to a mechanism goes from this reference point of the pulley through which this cable arrives to the attachment point if the cable stops at this mechanism, or to the reference point of the exit pulley if the cable also passes through higher mechanisms. Note that when the mechanism is oriented at an angle α, the reference point on the exit pulley is moved at the same angle with respect to the vertical.

First case

The side actuated part is summarized in Figure 2 (left mechanism). The red cable pulls on the left side, is attached in the middle of the pulley 3 and wraps around the pulley 1. Symmetrically, the blue cable pulls on the right side, is attached to the middle of the pulley 4 and wraps around the pulley 2.

The partial length of this cable that depends only on this mechanism i can be written as follows:

l c1 (α i ) = l t l (α i ) + rβ 1 1 (α i ) (14) r ll ll t R A1 A2
Figure 3: Actuation of a mechanism on the side: lengths and angles useful for the calculation. The red line represents the vertical axis. The blue line corresponds to the portion of cable between the middle of the upper pulley and its point of detachment from the lower pulley. The R point is taken as the reference point as the cable end point, and depends only on the lower mechanisms. The winding angle of the cable on the pulley is therefore

-A 1 + A 2 + π.
where l t l (α i ) is the length of the cable between the center of the pulley 3 and the place where the cable contacts the pulley 1, and β 1 1 (α i ) is the angle of winding of the cable on the pulley 1 from the reference point. The Figure (3) shows the lengths and angles useful for this calculation.

The length l t l (α i ) is obtained with the Pythagorean theorem (from the right angle shown in the Figure (3)):

l t l (α i ) = l l (α i ) 2 -r 2 (15)
This right angle is also used to calculate the angle A 1 :

A 1 = atan( l l (α i ) 2 -r 2 r ) (16)
The angle A 2 is the opposite of half of the orientation α i of the mechanism:

A 2 = - α i 2 (17)
The angle β 1 1 (α i ) is therefore:

β 1 1 (α i ) = π -atan( l l (α i ) 2 -r 2 r ) - α i 2 (18)
The equation ( 14) can therefore be rewritten:

l c1 (α i ) = l l (α i ) 2 -r 2 + r(π -atan( l l (α i ) 2 -r 2 r ) (19)
(-π 2 -α + π -A 2 ).

Second case

The second case is the symmetrical of the first one (on the right side instead of the left side). With the same reasoning, the partial length of the cable is therefore:

l c2 (α i ) = l t r (α i ) + r(π + α i 2 -atan( l r (α i ) 2 -r 2 r ) (20) with l t r (α i ) = l r (α i ) 2 -r 2 (21)

Third case

The wiring studied is that of the mechanism in the middle of Figure 2. The cable is wrapped around pulley 1 and then around pulley 3, linking the two from the outside.

The portion of cable that connects the 1 and 3 pulleys is l l (α i ) (see Figure ( 4)). So we have:

l c3 (α i ) = l l (α i ) + r(β 3 1 (α i ) + β 3 3 (α i )) (22)
where β 3 k (α i ) is the winding angle of the cable around the pulley k which depends only on the orientation of the mechanism i.

According to the Figure ( 4), we have:

β 3 1 (i) = π 2 - α i 2 (23) r r L/2 L/2 R A3 A4 A3 A4 T - r T P A5 A5 r
Figure 5: Actuation of a mechanism (fourth case): lengths and angles useful for the calculation. The notations are spread over two drawings to avoid overloading and to facilitate comprehension. The red line represents the horizontal axis. The blue line corresponds to the portion of cable between the detachment point on one pulley and the attachment point on the next pulley (either on the left or on the right). Whether the cable length along the diagonal bar or along the horizontal bar, the same method as in the previous cases is used to calculate their length, which is constant. The point R is taken as the reference point as the end point of the cable. The point T is a point through which the cable passes and which allows to make the junction between the two designs. The point P is taken as reference for the point where the cable exits. The winding angle of the cable on the lower pulley is

( π 2 + A 4 -A 3 ).
The angle of the cable on the upper left pulley to the T point is

(-α + A 4 -A 3 ) (left figure)
. The angle of winding on the same pulley from the point T is (-A 5 + π) (right figure). Finally, the winding angle of the cable on the upper right pulley is

(-A 5 + π 2 ).
According to the same figure, the winding angle on the upper pulley is:

β 3 3 (i) = 3π 2 - α i 2 (24)
We deduce that:

l c3 = l l (α i ) + r(2π -α i ) (25)
Fourth case The strut-routed part is summarized in Figure 5. The cable passes through the right side of pulley 2, winds from the left side on pulley 3 and then winds from the bottom of the pulley 4.

First of all, it can be geometrically demonstrated that the portions of cables that link the pulleys have a constant length (calculable with the Pythagorean theorem). The length linking the lower right and upper left pulleys is thus 2 ( L 2 ) 2r 2 , and the length linking the upper pulleys is 2 ( b 2 ) 2r 2 . So we can write:

l c4 (α i ) = 2 ( L 2 ) 2 -r 2 + 2 ( b 2 ) 2 -r 2 + r(β 2 2 (α i ) + β 2 3 (α i ) + β 2 4 (α i )) (26)
where β 2 k (α i ) is the winding angle of the cable around the pulley k which depends only on the orientation of the mechanism i.

The angles represented on the Figure [START_REF] Kraus | Pulley friction compensation for winch-integrated cable force measurement and verification on a cable-driven parallel robot[END_REF] are

A 3 = atan( L 2r ), A 4 = ψ i , A 5 = atan( b 2r ).
So we have:

β 2 2 (i) = π 2 + ψ i -atan( L 2r ) (27) β 2 3 (i) = -α i + ψ i -atan( L 2r ) -atan( b 2r ) + π (28) β 2 4 (i) = -atan( b 2r ) + π 2 (29)
We can therefore deduce that:

l c4 (α i ) =2 ( L 2 ) 2 -r 2 + 2 ( b 2 ) 2 -r 2 +r π 2 + ψ i -atan( L 2r ) -α i + ψ i -atan( L 2r ) -atan( b 2r ) + π -atan( b 2r ) + π 2 (30)
This result happens to be symmetrical in α. This implies that the symmetric strut-routed case has the same relative cable length.

Synthesis

The unwound lengths of the cables of our manipulator are therefore (starting from the bottom of the manipulator):

l c (1) = l c2 (α 1 ) (31) l c (2) = l c4 (α 1 ) + l c2 (α 2 ) (32) l c (3) = l c4 (α 1 ) + l c4 (α 2 ) + l c2 (α 3 ) (33) l c (4) = l c3 (α 1 ) + l c3 (α 2 ) + l c1 (α 3 ) (34)
Thus, the matrix Z can be written in the following form:

Z = (Z 1 Z 2 ) (35)
where Z 1 is a triangular matrix:

Z 1 =     ∂ l c2 (α 1 ) ∂ α 1 ∂ l c4 (α 1 ) ∂ α 1 ∂ l c4 (α 1 ) ∂ α 1 0 ∂ l c2 (α 2 ) ∂ α 2 ∂ l c4 (α 2 ) ∂ α 2 0 0 ∂ l c2 (α 3 ) ∂ α 3     (36)
and Z 2 a vector:

Z 2 =     ∂ l c3 (α 1 ) ∂ α 1 ∂ l c3 (α 2 ) ∂ α 2 ∂ l c1 (α 3 ) ∂ α 3     (37)

Dynamic model

We now present the dynamic model of our system. First, we will present a basic dynamic model, assuming that the cables are rigid. Then, we will present the dynamic model including the elasticity of the cables and the motor friction.

Basic dynamic model

Assuming that the cables are rigid and neglecting friction, a first dynamic model can be written as follows [2]:

M(α) α + c(α, α) + g(α) = Z(α)f ( 38
)
where g is the vector of forces deriving from a potential energy (gravity, spring forces

), f = [ f 1 , f 2 , f 3 , f 4 ]
T is the vector of the forces applied by the motors and transmitted by the cables, M is the matrix of intertia taking into account the inertia of the manipulator M m [START_REF] Benjamin | Dynamic modeling and control of a tensegrity manipulator mimicking a bird neck[END_REF] and the inertia of the motors:

M(α) = M m (α) + I a ( R r d ) 2 Z(α) T (39)
where I a is the inertia of the motors. c is the vector of Coriolis effects, defined by c m and the effects of the motors:

c = c m (α) + I a ( R r d ) 2 Ż(α, α)Z(α) T (40)
The control law tested in [2] was based on this dynamic model, but large oscillations in the applied forces were observable and the values of α calculated for the control were different from the observable values. In order to have a better control, an improved dynamic model is proposed below, including the motor friction and the elasticity of the ropes.

No attempt will be made to identify the parameters of the above dynamic model. Indeed, all data is available from the CAD model of our prototype and the manufacturer's motor data sheets. In order to have high available forces f, a high reduction ratio is used in the motor gearboxes (R = 25). Moreover, since tensegrity systems are inherently light, the contribution of the manipulator inertia is negligible compared to that of the gearmotors. On the other hand, we note that the latter is taken into account via the Z matrix, which depends on the cable routing and α. The inertia matrix M is therefore variable and it is important to take this into account in our models and in the order.

Dynamic model with friction and cable elasticity

All cables are now considered elastic. Assuming that θ i = 0, i = 1, ..., 4 when the manipulator is in α i = 0, i = 1, ...3 at rest (i.e. with zero applied forces and, therefore, with zero cable elongations), the vector of elastic cable elongations x c = [x c1 , x c2 , x c3 , x c4 ] T is defined as:

x c = l c -l c0 - r d R θ (41)
In other words, the elongation of each cable is defined as the difference between its theoretical length as calculated in section 4.2 and the length of cable wrapped around the associated winch.

Differentiating Eq. ( 41) with respect to time gives :

Z(α) T α = r d R θ + ẋc (42)
Let t c = [t c1 ,t c2 ,t c3 ,t c4 ] T be the vector of cable tensions and f f riction be the vector of friction forces. Most friction occurs in motor gearboxes. Modeling friction in pulleys can be interesting to improve the accuracy of cable force transducer measurements [START_REF] Kraus | Pulley friction compensation for winch-integrated cable force measurement and verification on a cable-driven parallel robot[END_REF]. In this article, no cable force transducers are used. Therefore, we focus on the friction coming from motor gearboxes. The dynamic model written for the system composed of motor gearboxes and cables provides the equation of motion of θ :

I a θ + r d R t c + r d R f f riction = r d R f (43)
The dynamic model written for the manipulator and the cables gives the equation of motion of α :

M m (α) α + c m (α, α) + g(α) = Z(α)t c (44)
By differentiating the equation ( 42) with respect to time and substituting θ for the equation ( 43), the tension in the cables can be expressed as follows :

t c = f -f f riction -I a ( R r d ) 2 (Z(α) T α + Ż(α, α) T α -ẍc ) (45)
Substituting the above expression in Eq (44) gives the dynamic model of the complete system with the friction and elasticity of the cable :

M(α) α + c(α, α) + g(α) = Z(α) f -f f riction + I a ( R r d ) 2 ẍc ( 46
)
Now it is necessary to identify the friction and elasticity of the cables.

6 Identification

Friction formulas

As we have written previously, it is assumed that friction occurs mainly in geared motors. A modeled dry friction can be written as follows:

f dry (f s , θ ) = diag(sign( θ ))f s (47)
where f s = [ f s (1), f s (2), f s (3), f s (4)] T is a constant vector that must be identified, diag(v) is a diagonal matrix of size (4 × 4) constructed with the components of v and sign is the sign function. A model ensuring continuity can be derived as follows:

f dry (f s , θ ) = diag( 2 π atan(c m θ ))f s (48)
where c m is a setting parameter that adjusts the slope of the arc tangent around 0. This parameter is arbitrarily set to c m = 0.5. The factor 2 π ensures that f f riction tends towards f s and -f s when θ has positive and negative values, respectively.

Viscous friction can be modeled as follows:

f visc (f v , θ ) = diag( θ )f v ( 49
)
where

f v = [ f v (1), f v (2), f v (3), f v (4)
] T is a constant vector that must be identified. Thus, the friction of the engines are:

f f riction (f s , f v , θ ) = f dry (f s , θ ) + f visc (f v , θ ) (50)
Substituting f f riction for Eq (46), we get :

Z(α)(diag( 2 π atan(c m θ ))f s + diag( θ )f v ) = Z(α)(f + I a ( R r d ) 2 ẍc ) -M(α) α -c(α, α) -g(α) (51)
We thus obtain a system of three equations that are linear in f s and f v . Once α, θ , their derivatives and f are known, we can identify f s and f v .

Cable elasticity formulas

During the identification experiments, α is obtained from the image processing and θ is measured by the engine encoders. The cable elongation x c can be calculated with Eq. (41), namely

x c = l c (α) -l c0 - r d R θ (52)
This calculation assumes that there is no cable elongation in the initial configuration of the right-hand manipulator. In practice, the cables must be slightly stretched to prevent slackening, but the resulting elongation is small and can be neglected. The ẋc derivative of x c is calculated numerically.

We will study three models of cable stiffness. The first is a general cable model that allows us to identify a set of parameters for each cable:

t ci = k(i)x ci + c a (i) ẋci (53)
or, as a vector :

t c = diag(x c )k + diag(ẋ c )c a ( 54
)
where k = [k(1), k(2), k(3), k(4)] T is the constant vector of cable stiffnesses and c a = [c a (1), c a (2), c a (3), c a (4)] T is the constant vector of damping coefficients. Using the equation ( 43) and the friction model, the tension in the cable can be written :

t c = f -f f riction - R r d I a θ (55)
Combination of eqs. ( 54) and (55), we obtain the following vector equation :

diag(x c )k + diag(ẋ c )c a = f -f f riction - R r d I a θ (56)
This vector equation provides four equations that are linear in the coefficients of friction, stiffness and damping of each cable.

The second and third elastic models are based on the Young's modulus [START_REF] Mohammad | Dynamic analysis and control of fully-constrained cable robots with elastic cables: variable stiffness formulation[END_REF]:

t c = diag(x c ) ES l + diag(ẋ c ) cS l ( 57
)
where S is the section of the cable. E is Young's modulus and c is a damping coefficient, both to be identified. l is the length of the cable: for the second model we will take l = l 0 the total length of the cable, while in the third model we will take l = l c the unwound length of the cable. Using the same method as in the first model, we obtain:

diag(x c ) ES l + diag(ẋ c ) cS l = f -f f riction - R r d I a θ (58)
This vector equation also provides four equations that are linear in the coefficients of friction, stiffness and damping for each cable.

Full identification model

The equations ( 51) and (56), which describe the contribution of friction in the dynamic model and the stiffness model of the cables, respectively, are now assembled to obtain a system of linear equations with the parameters to be identified :

Dv = y ( 59
)
where v is the vector of n p = 12 parameters to be identified v = k c a f s T . Now, the contribution of friction in Eqs. ( 51) and ( 56) are not of the same order because the coefficients of Z are of the order of 10 -2 . To avoid any discrepancy, Eq (56) is multiplied by a weight W = 100. The vector y is given by:

y =   f -R r d I a θ W (Z(f + I a ( R r d ) 2 ẍc ) -M α -g -c)   (60)
and D is written:

D = diag(x c ) diag(ẋ c ) diag( 2 π atan(c m θ )) 0 3×4 0 3×4 W Zdiag( 2 π atan(c m θ ) (61)
Friction parameters are also independently identified for comparison with the overall identification. In this case, there are only parameters n p = 8 to be identified, v = [f s f v ] T and :

y = f - R r d I a θ (62) D = [diag( 2 π atan(c m θ )) diag( θ )]. (63)

Method of identification

Motive forces f and positions θ and orientation angles α are sampled along a specific trajectory to obtain r sampled values. Then, v can be estimated as the least square solution v of the linear system :

Wv + ρ = Y ( 64
)
where W is an observation matrix (r × n p ) defined by sampling the regressor D defined in (61) or (63), Y is a sampling of y defined in (60) or (62), ρ is the vector of residuals or the vector of errors. The number r of samples must be greater than the number of parameters to be identified: r > n p [3].

The vector parameter v minimizing the standard 2 of ρ is given by : The quality of the solution depends on the W conditioning, which depends on the choice of the trajectory made during the identification process. If the measurement noises are independent in W and Y , the quality of the identification can be evaluated with the covariance matrix C X :

v = W + Y (65)
C X = σ 2 ρ (W T W) -1 (66)
where σ 2 ρ = ||Y-Wv|| 2 r-N p . The standard deviation of parameter i is σ i = C X (i, i) and the confidence interval is given by 2σ i . The relative standard deviation 100 σ i v i evaluates the quality of the identification of the parameter i.

Filtering

The data used for identification are the measurements of the angular position of the motor θ , the input actuating forces f (assuming perfect electrical response of the motors) and the manipulator configuration α evaluated from the image processing. Cable elongations are deduced from (41), after initialization in straight position. The derivatives of the measured quantities are calculated by centered numerical differentiation after filtering.

It is essential to filter the measured data to avoid biases resulting from discrete and noisy measurements (41). A Butterworth low-pass filter is used to filter forward and backward measurements with the Matlab function filtfilt. The derivatives of θ are obtained with a central difference algorithm to avoid phase shift. The cutoff frequency of the Butterworth filter is calculated with the method explained in [START_REF] Tu Pham | Identification of joint stiffness with bandpass filtering[END_REF].

The video frequency used to measure the orientation of the mechanisms is 30Hz, while the frequency of the control law is 500Hz. To synchronize the data, the motor encoder measurements are resampled with the Matlab decimate procedure at 31.25Hz, the value closest to 30Hz. The video data is also resampled to 31.25Hz with linear interpolation.

Results

Identification of dry and viscous friction

At first, only friction is identified, and in two different ways:

• On the one hand, dry and viscous friction is identified.

• On the other hand, only dry friction is identified, neglecting viscous friction.

The identified parameters are in the Table 1. Dry friction is well identified in both cases, but viscous friction is not: its relative standard deviation is well above 10%.

The Figure [START_REF] Tu Pham | Identification of joint stiffness with bandpass filtering[END_REF] shows the plots of the Z(α)f couples from the forces and their reconstruction with the equation (51), with and without the friction model, and with and without the viscous friction. The objective is to highlight the contribution of friction terms and to show that these terms must be taken into account. However, it shows that viscous friction, in addition to being poorly identified, does not allow for a better estimation of couples. For this reason, we will only consider dry friction in the rest of this report.

Identification of cable elasticity

Here, the dry friction and elasticity of the cables are identified, with the three elastic models offered.

The identified parameters are in the Table 2. The dry friction is well identified and similar whatever the elastic model. Stiffness parameters are well identified in all three models, but damping parameters are poorly identified in all three cases. This is due to the low values of the elongation rate.

There are two ways to calculate the tension in cables: by the applied forces and friction on the one hand, and by the elastic model and elongation on the other hand. The Figure (7) shows the tension in the cables. calculated by the forces applied and the friction identified (taking those identified in the previous section, all values identified here being strongly similar), and the tension calculated by the elastic models identified. It is shown that the three elastic models have different results but none of them can faithfully reproduce the tension calculated by the forces. The mean error of reconstruction is 6.88N for Model 1, 9.84N for Model 2 and 9.64N for Model 3.

The Figure (8) shows the reconstruction error of each model as a function of the tension calculated by the applied forces. It shows that for each model and each cable, the error increases as the tension increases. Nonlinear cable models should therefore be more suitable.

Conclusion

In this report, we have first detailed the calculations of the cable lengths, which are used in the dynamic model. These calculations are detailed without and with consideration of the pulleys, which have an impact on the shape of the matrix Z. Indeed, in comparison with the calculation without pulleys, the pulleys modify the non-zero values of Z and add some, corresponding to the strut-routed passage of the cables.

We then presented a dynamic model and added a friction model and an elasticity model of the cables.

In order to exploit these new models, the friction parameters were identified. The dry frictions are well identified and allow a better estimation of the forces by the dynamic model. Viscous friction is poorly identified and does not improve the model, so it is neglected.

Three elastic cables models were then compared, identifying the parameters specific to each one. For each model, the stiffness coefficients are well identified while the damping coefficients are not. All three models show interesting results, but they are less good when the tension in the cabes is high. Non-linear elastic models could therefore be more efficient. [2] Benjamin Fasquelle et al. "A bio-inspired 3-DOF light-weight manipulator with tensegrity X-joints". In: ICRA'2020. 2020.

[3] Wisama Khalil and Etienne Dombre. Modeling, identification and control of robots. Butterworth-Heinemann, 2004.

Figure 1 :

 1 Figure 1: Wiring of the manipulator being studied.

Figure 2 :

 2 Figure 2: The different cabling studied.

Figure 4 :

 4 Figure4: Actuation of a mechanism (third case): lengths and angles useful for the calculation. The red line represents the vertical axis. The blue line corresponds to the portion of cable between the detachment point on the lower pulley and the attachment point on the upper pulley. It is the same length as the black line connecting the centers of the two pulleys, because these two segments and the right angles formed with the tangent radii of the pulleys form a rectangle. The R point is taken as the reference point as the point of arrival of the cable, and depends only on the lower mechanisms. The winding angle of the cable on the lower pulley is therefore (-π 2 + A 2 + π). The winding angle on the upper pulley is (-π 2α + π -A 2 ).

Figure 6 :

 6 Figure 6: Plots of Z(α)f with the dynamic model, including friction or not, and from measurement data. Torques evaluated with friction model are clearly closer to the measured ones. The viscous friction does not modify a lot the curve in comparision with the sole dry friction.

Figure 7 :

 7 Figure 7: Tension in the cables, calculated by forces and elastic models.

Figure 8 :

 8 Figure 8: Error in reconstructing the tension in the cables, according to the tension calculated by the applied forces.

Table 1 :

 1 Identified friction parameters. The viscous parameters are badly identified.

		Dry friction Viscous friction		Dry friction only
	f s (1)(N)	12.91	4.68 × 10 -2	f s (1)(N)	14.13
	relative σ i	5.62%	56.85%	relative σ i	2.52%
	f s (2)(N)	16.72	4.63 × 10 -2	f s (2)(N)	17.70
	relative σ i	4.55%	72.11%	relative σ i	2.02%
	f s (3)(N)	18.11	4.64 × 10 -2	f s (3)(N)	18.97
	relative σ i	4.25%	83.42%	relative σ i	1.88%
	f s (4)(N)	16.09	0.72 × 10 -2	f s (4)(N)	16.29
	relative σ i	1.85%	105.96%	relative σ i	1.25%

Table 2 :

 2 Identified stiffness and damping parameters for three mechanisms.

		Model 1 Model 2 Model 3
	f f f s (1)(N)	14.42	14.38	14.32
	relative σ	1.68%	1.73%	1.72%
	f f f s (2)(N)	17.65	17.49	17.48
	relative σ	1.38%	1.44%	1.42%
	f f f s (3)(N)	18.77	18.41	18.48
	relative σ	1.32%	1.38%	1.35%
	f f f s (4)(N)	16.18	16.32	16.29
	relative σ	0.88%	0.90%	0.89%
	k k k(1) (N.m -1 )	10454	-	-
	relative σ	1.37%	-	-
	k k k(2) (N.m -1 )	9256	-	-
	relative σ	1.60%	-	-
	k k k(3) (N.m -1 )	8319	-	-
	relative σ	1.85%	-	-
	k k k(4) (N.m -1 )	5336	-	-
	relative σ	1.28%	-	-
	E (GPa)	-	4.04	1.23
	relative σ	-	0.83%	0.81%
	c c c a (1) (N.s.m -1 )	805	-	-
	relative σ	30.5%	-	-
	c c c a (2) (N.s.m -1 )	333	-	-
	relative σ	68.7%	-	-
	c c c a (3) (N.s.m -1 )	-218	-	-
	relative σ	94.3%	-	-
	c c c a (4) (N.s.m -1 )	370	-	-
	relative σ	64.3%	-	-
	c (MPa.s -1 )	-	139.8	57.4
	relative σ	-	45.5%	32.88%


Elsa Cattelain-Lopez

Delphine Chu

Miow Lin

André Happe

Philippe Goupille

Emmanuel Oger

Denis Mulleman
email: denis.mulleman@univ-tours.fr

Poor assessment of bone mineral density after a forearm fracture in women aged 50 years or older: Data from a French health insurance database

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Letter to the Editor,

Poor assessment of bone mineral density after a forearm fracture in women aged 50 years or older: data from a French health insurance database Forearm fractures should be considered as an early sign of osteoporosis. [START_REF] Klotzbuecher | Patients with prior fractures have an increased risk of future fractures: a summary of the literature and statistical synthesis[END_REF] In 2010, 47,131 forearm fractures occurred in France, a country that comprises 22.5% women ≥ 50 years old and 2.8 million osteoporotic patients. [START_REF] Hernlund | Osteoporosis in the European Union: medical management, epidemiology and economic burden. A report prepared in collaboration with the International Osteoporosis Foundation (IOF) and the European Federation of Pharmaceutical Industry Associations (EFPIA)[END_REF] In a recent meta-analysis, the risk of a subsequent hip fracture was reported to be increased in women who presented a forearm fracture. [START_REF] Johnson | Risk of hip fracture following a wrist fracture-A meta-analysis[END_REF] Thus, bone mineral density (BMD) testing is recommended after a forearm fracture to assist clinicians in decision making. [START_REF] Briot | update of French guidelines for the pharmacological treatment of postmenopausal osteoporosis[END_REF] We previously reported that in our district, only 10% of women ≥ 50 years old had BMD testing after a forearm fracture. [START_REF] Erny | Management of osteoporosis in women after forearm fracture: data from a French health insurance database[END_REF] However, the size of the population sample in that study was relatively small and it was restricted to the private sector. We used the French national health insurance information system (Système National d'Information Inter Régimes de l'Assurance Maladie [SNIIRAM]), [START_REF] Martin-Latry | Pharmacoepidemiological research using French reimbursement databases: yes we can![END_REF] which can be linked to the French medical program for information systems (Programme de Médicalisation des Systèmes d'Information [PMSI]), a medico-economic tool used for describing public or private hospital activity. The study included 4,120 women. Among the 4,120 women ≥ 50 years old who had a forearm fracture between 2011 and 2012 in the Centre-Val de Loire area, 546 (13.3%) had a BMD assessment in the year after the forearm fracture (Table 1). The median delay between the fracture and BMD testing was 117.5 days (range 4-365). Among the 3574 women without BMD testing, 231 (6.5%) initiated or switched an anti-osteoporotic drug in the year after the fracture versus 168/546 women (30.8%) who had BMD testing (OR=6. This finding confirms previous works and stresses the need to implement these guidelines in clinical practice. [START_REF] Erny | Management of osteoporosis in women after forearm fracture: data from a French health insurance database[END_REF][START_REF] Viprey | Low Osteoporosis Treatment Initiation Rate in Women after Distal Forearm or Proximal Humerus Fracture: A Healthcare Database Nested Cohort Study[END_REF] BMD testing was less frequent for women with than without comorbidities. This finding may be explained by a lesser attention given to the prevention of osteoporosis as compared with other diseases. In our study, aromatase inhibitors was associated with increased frequency of BMD assessment which is in good agreement with the recommendations for patients who receive such drugs. We found no association between accessibility to BMD testing and proportion of BMD tests performed. Since there is no apparent bias in identification of BMD assessment, the availability of BMD machines does not explain these inter-district differences.

To prevent secondary fracture, several successful initiatives such as Fractures Liaisons Service

(

  FLS) were promoted.[8, 9] FLS have significantly decreased the fracture rate and increased anti-osteoporotic drug prescription in one study.[10] FLS, along with local communication campaigns, may help prevent osteoporosis after a forearm fracture in women. In 2017 we launched a FLS after forearm fracture in our hospital. Based on the present work, we will investigate the further trend of BMD assessment in our setting.

A c c e p t e d m a n u s c r i p t [ 8 ]

 8 Eisman JA, Bogoch ER, Dell R, Harrington JT, McKinney RE, Jr., McLellan A, et al. Making the first fracture the last fracture: ASBMR task force report on secondary fracture prevention. J Bone Miner Res 2012;27(10):2039-46. [9] Ganda K, Puech M, Chen JS, Speerin R, Bleasel J, Center JR, et al. Models of care for the secondary prevention of osteoporotic fractures: a systematic review and meta-analysis. Osteoporos Int 2013;24(2):393-406. [10] Lih A, Nandapalan H, Kim M, Yap C, Lee P, Ganda K, et al. Targeted intervention reduces refracture rates in patients with incident non-vertebral osteoporotic fractures: a 4-year prospective controlled study. Osteoporos Int 2011;22(3):849-58.

A c c e p t e d m a n u s c r i p t

Table 1 )

 1 . BMD testing was not associated with corticosteroids exposure on univariable analysis [35/297 (11.8%) versus 511/3823 (13.4%), OR=0.9, 95% CI 0.6-1.2, p=0.439], but was associated with aromatase inhibitors [21/79 (26.6%) versus 525/4041 (13.0%), OR=3.2, 95% CI 2.0-5.3, p<0.001]. BMD testing was performed most frequently in the Loir-et-Cher district (109 BMD tests, 19.2%) and the least in Eure-et-Loir (55 BMD tests, 9.8%). By contrast, accessibility to BMD testing was lower in Loir-et-Cher (3.97 densitometers per 100,000 women ≥ 50 years old) and higher in Eure-et-Loir and Loiret districts (14.23 and 13.12

	4, 95% CI

The comorbidities in patients with or without BMD are presented in the supplemental table.

Table 1 :

 1 Characteristics of women ≥ 50 years old with a forearm fracture in 2011 and 2012 in Centre-Val de Loire area in France by bone mineral density (BMD) testing in the year after the Supplemental material: percentage of bone mineral density assessment according to the presence or absence of comorbidities in women ≥ 50 years old in Centre-Val de Loire area in France

	fracture

Supplementary Material

Click here to access/download;Supplementary Material;Cattelain-Lopez_supplemental_table.docx

A c c e p t e d m a n u s c r i p t

Acknowledgements: We thank Caroline Rault (CHU de Rennes, UPRES-EA 7449 REPERES) ; Florence Erny, François Jacquot and Philippe Tauveron (Centre de l'Ostéoporose, Tours) ; Ken Haguenoer and Emmanuel Rusch (Service d'Information Médicale, d'Epidémiologie et d'Economie de la Santé, Laboratoire de santé publique, CHRU de Tours). This work would not have been possible without the support of Mehdi Gabbas, CNAM-TS (Caisse Nationale de l'Assurance Maladie des Travailleurs Salariés) (Paris), Valérie Edel, INDS (Institut Nationale des données de Santé). This article is the result of the "Prevention of osteoporosis after distal forearm fracture (Prévention de l'Ostéoporose après fracture du Poignet [POP])" study group.

A c c e p t e d m a n u s c r i p t

Funding: This work was supported by grants from the Osteoporosis Research and Information Group (GRIO -Groupe de Recherche et d'information sur les ostéoporoses) in 2013 and the French League Against Rheumatism (AFLAR -Association Française de Lutte Anti-Rhumatismale) in 2014.

Ethics approval and consent to participate: Not applicable

Consent for publication: Not applicable

Availability of data and materials:

The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.