Haptic-based Shared-Control Methods for a Dual-Arm System Mario Selvaggio 1 , Firas Abi-Farraj 2 , Claudio Pacchierotti 2 , Paolo Robuffo Giordano 2 and Bruno Siciliano 1 Abstract-We propose novel haptic guidance methods for a dual-arm telerobotic manipulation system, which are able to deal with several different constraints, such as collisions, joint limits, and singularities. We combine the haptic guidance with shared-control algorithms for autonomous orientation control and collision avoidance meant to further simplify the execution of grasping tasks. The stability of the overall system in various control modalities is presented and analyzed via passivity arguments. In addition, a human subject study is carried out to assess the effectiveness and applicability of the proposed control approaches both in simulated and real scenarios. Results show that the proposed haptic-enabled shared-control methods significantly improve the performance of grasping tasks with respect to the use of classic teleoperation with neither haptic guidance nor shared control.

Index Terms-Haptics and Haptic Interfaces; Telerobotics and Teleoperation; Industrial Robots.

I. INTRODUCTION

O VER the last few years remote telemanipulation has shown significant advancements and promising results have been achieved in several fields such as minimally-invasive robotic surgery [START_REF] Freschi | Technical review of the da vinci surgical telemanipulator[END_REF], sort and segregation of nuclear waste [START_REF] Abi-Farraj | A visual-based shared control architecture for remote telemanipulation[END_REF], [START_REF] Esfahani | Human-in-the-loop optimisation: mixed initiative grasping for optimally facilitating post-grasp manipulative actions[END_REF], telemaintenance [START_REF] Aschenbrenner | Teleoperation of an Industrial Robot in an Active Production Line[END_REF], and micromanipulation [START_REF] Pacchierotti | Steering and control of miniaturized untethered soft magnetic grippers with haptic assistance[END_REF], [START_REF] Pacchierotti | Haptic feedback for microrobotics applications: a review[END_REF]. In this respect, the nuclear industry is one of the most exciting fields of application for robotic telemanipulation, where several tasks, such as sort and segregation of nuclear waste, can clearly benefit from advanced telemanipulation techniques.

Current telerobotic systems designed for such tasks provide teleoperation capabilities through extremely primitive master consoles (e.g., passive joystick or teach pendants), making these operations prohibitively slow to process large amounts of material in a reasonable time. Besides being time demanding, these tasks usually require highly-skilled human operators. Indeed, steering a remote manipulator toward a desired grasping pose is a quite complex task for an operator directly controlling the 6-degrees-of-freedom (DoF) pose of a manipulator endeffector. This is due to (i) the complexity of regulating both the position and the orientation at the same time and (ii) the presence of several constraints (e.g., collisions, joint limits, singularities) that further limit the operator's maneuvering dexterity (but of which the operator has no direct or intuitive awareness). To overcome these limitations, the user needs to constantly pay close attention to the status of the robotic system, which can be sometimes difficult and cognitively demanding (e.g., for singularity or joints limit avoidance). A possible way to reduce the operator cognitive and physical workload is to exploit the sensory information collected at the slave side to design novel haptic-guided and shared control-enabled teleoperation systems.

Kinematic limitations avoidance has been recently applied to collaborative robots to improve the performance and the intuitiveness of physical human-robot interaction [START_REF] Lebel | An anticipative kinematic limitation avoidance algorithm for collaborative robots: Threedimensional case[END_REF]. In robotic teleoperation this approach can be implemented through haptic guidance methods: instead of providing a high-fidelity haptic feedback, reflecting the actual physical contacts between a slave manipulator and the environment, the human operator is provided with an interface which informs her/him about proximity to constraints. Proposed solutions to this problem are based on potential fields [START_REF] Khatib | Real-time obstacle avoidance for manipulators and mobile robots[END_REF], virtual spring-damper systems [START_REF] Santis | The skeleton algorithm for self-collision avoidance of a humanoid manipulator[END_REF], and/or virtual fixtures [START_REF] Selvaggio | Enhancing bilateral teleoperation using camera-based online virtual fixtures generation[END_REF], [START_REF] Selvaggio | Vision based virtual fixture generation for teleoperated robotic manipulation[END_REF]. When the slave side is a complex robotic system (e.g., dual-arm, highly redundant) or environmental obstacles are present, (self-)collision avoidance methods need to be considered. In the past, several collision avoidance algorithms have been developed to implement reactive control strategies or to plan collision-free paths for redundant robots [START_REF] Dietrich | Integration of reactive, torque-based self-collision avoidance into a task hierarchy[END_REF]- [START_REF] Seto | Self-collision avoidance motion control for human robot cooperation system using robe[END_REF]. In this regard, self-collisions and joint limits have been used as criteria to find optimal inverse kinematic solutions for robotic manipulators [START_REF] Dariush | Constrained resolved acceleration control for humanoids[END_REF], [START_REF] Chan | A weighted least-norm solution based scheme for avoiding joint limits for redundant manipulators[END_REF]. The idea of combining kinematic constraints and self-collision avoidance to create safe and intuitive haptic-guided telerobotic system has been preliminarily introduced in [START_REF] Selvaggio | Towards a self-collision aware teleoperation framework for compound robots[END_REF].

On the other hand, several shared-control techniques for telerobotic systems have been proposed in the past [START_REF] Saeidi | A mixed-initiative haptic teleoperation strategy for mobile robotic systems based on bidirectional computational trust analysis[END_REF]- [START_REF] Khademian | Shared control architectures for haptic training: Performance and coupled stability analysis[END_REF]. For instance, in [START_REF] Abi-Farraj | A visual-based shared control architecture for remote telemanipulation[END_REF], a sensory-based shared-control architecture for remote manipulation is presented: a visualbased autonomous algorithm regulates a subset of the gripper's DoF to ease the approach toward an object to grasp. At the same time, the operator has control over the remaining null-space motions w.r.t. the primary (autonomous) task and is informed about possible constraints of the single-slave robotic system through haptic cues. However, all these previous works consider one slave arm operating in a static obstacle-free environment.

In this paper, we consider the case of a dual-arm robotic system: one robot is (partially) teleoperated by a human, while the other autonomously executes a task in the same workspace. Similar scenarios are considered in [START_REF] Nicolis | Occlusion-free visual servoing for the shared autonomy teleoperation of dual-arm robots[END_REF], [START_REF] Pedemonte | Visual-based shared control for remote telemanipulation with integral haptic feedback[END_REF], where an operator teleoperates a robotic arm while a camera, mounted on a second manipulator, keeps the scene visible. In [START_REF] Nicolis | Occlusion-free visual servoing for the shared autonomy teleoperation of dual-arm robots[END_REF], an occlusionavoidance algorithm controls the manipulator with the camera and keeps a target object in the field of view. No haptic guidance is considered. In [START_REF] Pedemonte | Visual-based shared control for remote telemanipulation with integral haptic feedback[END_REF], users are provided with haptic feedback about proximity to kinematic constraints of the manipulator they are controlling. However, the two manipulators are far enough to prevent any interference between them.

In this work, two manipulators (one of which is controlled by the operator) have independent and sometimes conflicting tasks. Moreover, they share a limited workspace, risking dangerous collisions. The primary goal of our paper is to define and assess a suitable haptic shared-control interface for this system which facilitates the operator's task and decrease her/his mental workload in avoiding system constraints. To this end, we adapt and improve the approaches presented in [START_REF] Abi-Farraj | A visual-based shared control architecture for remote telemanipulation[END_REF] and in [START_REF] Selvaggio | Towards a self-collision aware teleoperation framework for compound robots[END_REF], combining their features, extending their capabilities, and providing a novel, extensive human subject evaluation and results analysis. Specifically,

• we improved the shared-control approach [START_REF] Abi-Farraj | A visual-based shared control architecture for remote telemanipulation[END_REF] by (i) dividing the approach-to-grasp phase into multiple parts, according to the gripper distance from the target object, and (ii) devising different, optimized strategies for each of these parts; • we improved the self-collision-aware approach presented in [START_REF] Selvaggio | Towards a self-collision aware teleoperation framework for compound robots[END_REF] by (i) also considering collisions with the surrounding environment and (ii) replacing the convex meshes collision model of the robot with discretized sphere volumes; • we combined the proposed, new above mentioned techniques to create a novel shared-control framework able to manage a dual-arm system, where one robot is (partially) controlled by the operator and the other one autonomously performs a task in the same environment; • we presented a port-Hamiltonian model of the overall system subject to various control modalities, proving the passivity of the system w.r.t. the operator actions; • we carried out an extensive human subject evaluation in simulated and real environments, enrolling a total of 20 participants.

II. SYSTEM ARCHITECTURE Our telemanipulation system consists of two serial manipulators sharing a workspace in an industrial setting. One of the manipulators is performing an independent fully-autonomous task (in our case, a visual inspection with an onboard camera), while the second manipulator is (partially) teleoperated by a human operator whose aim is to grasp a target object (see Fig. 1). The operator faces a number of constraints while controlling such a system: kinematic/workspace constraints for the manipulator she/he is directly controlling, as well as possible interferences with the task performed by the fullyautonomous arm. As explained, the purpose of this work is to detail and assess a suitable shared-control interface to facilitate the operator's task and decrease her/his mental workload.

A. System Modeling

With reference to Fig. 1, the slave system consists of two 6-DoF velocity-controlled manipulator arms, A and B, with joint configurations q A ∈ R 6 and q B ∈ R 6 , respectively. Manipulator A is (partially) controlled by a human operator, while manipulator B autonomously performs a visual task with a potential overlap with the workspace of manipulator A. We also let

uA = vA ωA = JA qA, uB = vB ωB = JB qB (1)
be the Cartesian linear/angular velocities of the gripper (for manipulator A) and of the camera (for manipulator B) in their respective frames, and (J A , J B ) the associated geometric Jacobians. In the following, we will equivalently consider joint velocities ( qA , qB ) or Cartesian velocities (u A , u B ) as available control inputs, depending on the context. The master device is instead modeled as a generic, gravity pre-compensated, mechanical system,

Mm(xm) um + Cm(xm, um)um + Bmum = fm + f h , (2)
where x m ∈ R 6 represents the device pose in Cartesian space and u m ∈ R 6 the associated linear/angular velocity, M m (x m ) ∈ R 6×6 is the positive-definite and symmetric inertia matrix, C m (x m , u m ) ∈ R 6×6 consists of the Coriolis/centrifugal terms, B m accounts for the friction term, and f m , f h ∈ R 6 are the control and human forces, respectively. The control forces f m = [. . . f m,i . . .] T ∈ R 6 will be exploited to provide force cues to the user meant to inform of the feasibility of her/his commands against the constraints of the slave system (see Sec. II-D).

B. Constraints on the Slave Side

Constraints at the slave side are encoded in a suitable cost function H(q A , q B ), whose gradient w.r.t. the joint configuration vectors will be exploited to generate force cues f m provided to the operator and to implement the reactive behavior of manipulator B (see Sec. II-C).

1) Joint Limits: As customary, we consider that both manipulator arms are subject to joint limits. Among many possibilities, we adopt the following cost function to encode proximity to joint limits

hJ (q) = ρJ 6 i=1 e -α J (qi,max-qi) (qi,max -qi) -β J + e -α J (qi-qi,min) (qi -qi,min) -β J , ( 3
)
where q is the configuration vector of one of the slave arms, (q i,min , q i,max ) the i-th min/max joint limits, and (ρ J , α J , β J ) are scalar positive constants. The function h J (q) has a minimum at the center of the joint range and grows to infinity as one of the joint approaches its limits. 2) Singularities: As second constraint, we consider singularities. This constraint is encoded by the cost function [START_REF] Dariush | Constrained resolved acceleration control for humanoids[END_REF],

hS(q) = ρSe -α S det(ΓJ(q)) (det(ΓJ (q))) -β S , (4)
where, again, (ρ S , α S , β S ) are scalar positive constants, and q and J stand for the joint vector and Jacobian of the manipulator A or B, respectively, and Γ is a diagonal and positive-definite weight matrix that accounts for non-homogeneous physical dimensions of the Jacobian entries, namely, rotations and translations. The function h S (q) grows to infinity as det(J (q)) → 0 and it vanishes exponentially as det(J (q)) → ∞.

3) Collisions: As final constraint, we consider collisions between the two slave manipulators.

Any cost function encoding the collision avoidance constraint needs to exploit some measure of the minimum distance between the two arms. In order to facilitate the computation of the minimum distance between two bodies, we take inspiration from the line sphere-swept model [START_REF] Cascio | Smooth proximity computation for collision-free optimal control of multiple robotic manipulators[END_REF] (see Fig. 2a). However, as well known, the minimum distance between two (also convex) bounding volumes can have a discontinuous gradient when the pair of closest points on the two shapes suddenly jumps due to small perturbations of the shape locations. We therefore opted for discretizing the bounding volumes with a finite number of spheres, and exploiting all the possible intersphere distances to build the cost function associated to the collision avoidance (see Fig. 2b). An illustrative example of the resulting collision model is given in Fig. 2e, where one can appreciate how the adopted discretized-sphere model represents a fairly good approximation of the manipulator original and convex mesh models (Fig. 2c-d).

Let then d ij (q A , q B ) represent the minimum distance between the i-th sphere on manipulator A and the j-th sphere on manipulator B. The collision avoidance cost function is defined as

hC (qA, qB) = ρC S A i=1 S B j=1 e -α C d ij d -β C ij , ( 5
)
where S A and S B are the number of spheres used for discretizing the bodies of manipulators A and B, and (ρ C , α C , β C ) are scalar positive constants. One can easily verify that h C (q A , q B ) → ∞ as any of the inter-sphere distances vanishes while h C (q A , q B ) → 0 as all the inter-sphere distances become large enough. 4) Total Cost Function: Finally, the total cost function accounting for all the above mentioned constraints is H(qA, qB) = hJ (qA)+hJ (qB)+hS(qA)+hS(qB)+hC (qA, qB).

(

) 6
As explained, the gradient of H(q A , q B ) w.r.t. the joint vector q A can be used for cueing the operator about the feasibility of her/his commands against the constraints of the slave side, while the gradient of H(q A , q B ) w.r.t. the joint vector q B can be used for implementing a reactive behavior in manipulator B for avoiding possible collisions with manipulator A (see Sec. II-C). Note that these two actions (force cues and reactive behavior) are potentially coupled because of the mixed term h C (q A , q B ).

C. Slave Control

We start detailing the control architecture of manipulator B, which performs an autonomous visual task to keep an object of interest in visibility during the operation. Let s ∈ R 2 represent the image plane location of a representative point on the target object (see Fig. 1), s d a desired value for s, and L s ∈ R 2×6 the associated point feature interaction matrix [START_REF] Chaumette | Visual servo control, part i: Basic approaches[END_REF], such that ṡ = L s J B qB . Manipulator B is then controlled by employing the usual projected gradient control [START_REF] Siciliano | A general framework for managing multiple tasks in highly redundant robotic systems[END_REF] (hereafter we use H to indicate H(q A , q B ))

qB = kB1 (LsJB) † (s d -s) primary task (visual servoing) -kB2 (I -(LsJB) † (LsJB))∇q B H secondary task (minimize H) , (7)
with k B1 > 0 and k B2 ≥ 0. The primary task consists in the regulation of s toward s d . Its null-space is then exploited for the secondary task of minimizing the constraint cost function H in [START_REF] Pacchierotti | Haptic feedback for microrobotics applications: a review[END_REF]. This null-space action will keep manipulator B away from its singularities and joint limits, as well as maintain a safe distance from manipulator A. This action can be deactivated by setting k B2 = 0, asking manipulator B to carry out only the primary visual servoing task.

For manipulator A, we instead considered two possible control modalities: Full Teleoperation and Shared Control.

1) Full Teleoperation: In this modality, the user is given full control over the 6-DoF pose of the gripper on manipulator A. This is achieved by simply setting

uA = kAum, kA > 0 (8)
thereby implementing a classic velocity-velocity coupling between master and slave gripper linear/angular velocities.

2) Shared Control: The shared-control modality is meant to facilitate the user's operation in approaching the object to be grasped by letting an autonomous algorithm control a suitable subset of the gripper pose. This modality is split into two phases depending on the distance between the gripper and the object. When d < d th (threshold value), we adopt the shared-control strategy described in [START_REF] Abi-Farraj | A visual-based shared control architecture for remote telemanipulation[END_REF], whose aim is to autonomously keep the pointing direction toward the object

p = G PO G PO ∈ S 2 (9)
always aligned with the gripper approaching direction p d = Z G , while the operator is in control of the remaining DoF in the null-space of this primary task. By doing so, the gripper is always oriented toward the object, thus helping the user in completing the task. This behavior is obtained as follows: defining P p = I-pp T , [•] × the usual skew-symmetric operator, and

Lp = - 1 d Pp [p]× ∈ R 3×6 (10)
the interaction matrix associated to the task (9) (such that ṗ = L p u A ), the following controller is implemented on manipulator A,

uA = kA1L † p (p d -p) + kA2N Λum (11)
with k A1 ≥ 0 and k A2 > 0. Here, N = [n 1 . . . n 4 ] ∈ R 6×4 collects the four null-space directions

n1 = p 0 , n2 = 0 p , n3 = -[p]×ey -Ppey/d , n4 = [p]×ex Ppex/d , (12)
with e x = [1 0 0] T and e y = [0 1 0] T , which represent a basis of the four-dimensional null-space of the primary task, and Λ ∈ R 4×6 is a selection matrix for extracting four components1 of u m used by the operator for commanding a 'velocity' along the four motion directions spanned by N . Figure 3 provides a visual illustration of the four null-space motion directions associated to [START_REF] Dietrich | Integration of reactive, torque-based self-collision avoidance into a task hierarchy[END_REF], while further details can be found in [START_REF] Abi-Farraj | A visual-based shared control architecture for remote telemanipulation[END_REF].

On the other hand, when d ≥ d th , we do not constrain the gripper to be oriented toward the object, since, when far from the target object, it would unnecessarily limit the slave maneuverability. Therefore, when d ≥ d th , the shared-control algorithm switches to a modality in which the operator has control over the translational motion of the slave gripper, while the gripper orientation is autonomously controlled so as to minimize the cost function H and, thus, stay away as much as possible from the system constraints. This is obtained by replacing [START_REF] Selvaggio | Vision based virtual fixture generation for teleoperated robotic manipulation[END_REF] with

uA = kA1Sum -kA2ZJ -T A ∇q A H, (13)
where

Z = 0 0 0 I3 , S = I3 0 0 0 ∈ R 6×6 (14)
are selection matrices. As shown in the Appendix, controller ( 13) is indeed able to minimize H when Su m = 0 (i.e., in free translational motion) by acting on the angular velocity of manipulator A.

D. Haptic guidance

As mentioned before, the control forces f m on the master side (2) are exploited to inform the operator about the feasibility of her/his commands against the system constraints, encoded in the cost function H. This is obtained as follows, depending on the control modality considered for manipulator A:

1) Full Teleoperation: In this case, manipulator A is controlled by [START_REF] Khatib | Real-time obstacle avoidance for manipulators and mobile robots[END_REF], and the force cues are generated as

fm = -kM J -T A ∇q A H, kM > 0, (15)
the rationale being that (15) provides the force that, when applied to all the master device DoF, would make the master move so as to minimize H at the slave side. Therefore, cues [START_REF] Seto | Self-collision avoidance motion control for human robot cooperation system using robe[END_REF] provide the human operator with information about where to move in order to stay away from the constraints.

2) Shared Control: When manipulator A is close enough to the object (d < d th ), its behavior is controlled by [START_REF] Selvaggio | Vision based virtual fixture generation for teleoperated robotic manipulation[END_REF]. In this case, the force cues are generated as

fm = -kM Λ T N T J -T A ∇q A H, kM > 0. (16)
Similarly to the previous case, cues [START_REF] Dariush | Constrained resolved acceleration control for humanoids[END_REF] represent the forces that should be applied to the free directions of the master device for letting the manipulator A minimize H along the null-space directions spanned by N (where the operator can act). When manipulator A is instead far from the object (d ≥ d th ), its behavior is regulated by [START_REF] Dube | A model of the humanoid body for self collision detection based on elliptical capsules[END_REF], and the force cues are

fm = -kM SJ -T A ∇q A H, kM > 0, (17)
with again the idea of providing a force feedback that, when applied to the free master DoF (i.e., the translational ones in this case), would make manipulator A move so as to minimize H.

III. PASSIVITY ANALYSIS

Enforcing passivity is very important when dealing with teleoperation systems in order to guarantee a stable closed-loop behavior. Passivity has been extensively studied and applied to haptic feedback teleoperation systems, especially when time delays are present [START_REF] Lee | Passive bilateral teleoperation with constant time delay[END_REF]. In this section, our aim is to analyze the passivity of the considered teleoperation system subject to the control modalities introduced in Secs. II-C and II-D.

By considering the two velocity-controlled robots as simple integrators, the total energy of the system can be written as

V (pm, qA, qB) = 1 2 p T m M -1 m pm + H(q A , q B ) ( 18
)
where p m = M m u m is the haptic device momentum, q A , q B ∈ R 6 are the generalized coordinates of the two slave manipulators, respectively, and H(q A , q B ), introduced in ( 6), is the potential energy associated with the system constraints. We can thus analyze passivity by showing that the closedloop system equations can be put in port-Hamiltonian form. Formally, a port-Hamiltonian system (PHS) can be represented by the following set of equations

     ẋ = [J (x) -R (x)] ∂V ∂x + g (x) u y = g T (x) ∂V ∂x (19)
where x ∈ R n represents the system state, V (x) : R n → R is the Hamiltonian function, namely the sum of system energies, J (x) = -J (x) T represents the internal interconnection, R (x) ≥ 0 the internal dissipation, g (x) the input matrix, u the system input, and y the system output. It is easy to show that for PHS the following inequality holds

u T y = V (x) - ∂ T V ∂x R (x) ∂V ∂x ≤ V (x) . (20)
Equation ( 20) establishes the inherent passivity condition of a PHS with respect to the input-output pair (u, y) with storage function V(x). In the following, we then show that the three control modalities described in Sec. II lead to a closed-loop PHS formulation. For the analysis, we consider the primary task in ( 7) and [START_REF] Selvaggio | Vision based virtual fixture generation for teleoperated robotic manipulation[END_REF] to have reached a steady state (s → s d and p → p d ).

1) Full teleoperation: setting k M = k A = k, the closedloop system can be written as follows

ṗm qA qB =     0 -kJ -T A 0 kJ -1 A 0 0 0 0 0   + - Bm 0 0 0 0 0 0 0 P    ∇p m V ∇q A V ∇q B V    + 1 0 0 f h . ( 21
)
2) Shared control (d < d th ): setting k M = k A2 = k, the closed-loop system can be written as follows

ṗm qA qB =     0 -kΛ T N T J -T A 0 kJ -1 A N Λ 0 0 0 0 0   + - Bm 0 0 0 0 0 0 0 P    ∇p m V ∇q A V ∇q B V    + 1 0 0 f h . ( 22
)
3) Shared control (d ≥ d th ): setting k M = k A1 = k, the closed-loop system can be written as follows

ṗm qA qB =     0 -kSJ -T A 0 kJ -1 A S T 0 0 0 0 0   + -   Bm 0 0 0 kA2J -1 A ZJ -T A 0 0 0 P        ∇p m V ∇q A V ∇q B V    + 1 0 0 f h , (23)
where

P = I -(L s J B ) † (L s J B
) is a null-space projector such that P = P T ≥ 0. Therefore, since in all cases the closed-loop systems can be put in a PHS form (being the resulting interconnection and dissipation matrixes always skew-symmetric and positive semi-definite, respectively), one can conclude passivity of the three modalities w.r.t. the pair (u m , f h ) with energy function V(p m , q A , q B ) as sought 2 .

IV. EXPERIMENTS IN SIMULATION

A. Experimental setup, task, and participants

The experimental setup is shown in Fig. 4. The slave side is simulated using V-REP, and it is composed of our two velocity-controlled manipulators: an Adept Viper 850, controlled by the human operator (manipulator A), and an Adept Viper 650, controlled by an autonomous algorithm (manipulator B). The master side is composed of an Haption Virtuose 6D Desktop interface, which is used to control the Viper 850 robot and provides kinesthetic haptic feedback. Manipulator A is endowed with a ROBOTIQ 2-finger gripper, while manipulator B is equipped with a vision sensor. The remote environment is composed of two objects, a cube and a rectangular parallelepiped, placed on a conveyor belt. As detailed in Sec. II, the autonomous manipulator B robot is in charge of tracking the parallelepiped using standard visual servoing techniques [START_REF] Marchand | Visp for visual servoing: a generic software platform with a wide class of robot control skills[END_REF]. On the other hand, participants are required to control the motion of manipulator A to grasp the cube and lift it from the ground, avoiding collisions with the fully-autonomous robot. The task starts when the manipulator moves for the very first time, and it is considered successfully completed when the object is lifted from the ground. All the control policies described are implemented in MATLAB/SIMULINK and interfaced with V-REP using ROS through the matlab ros bridge (https://github.com/lagadic/ matlab ros bridge). The control loop runs at 100 Hz. A video of the experiment is available as supplemental material and at https://youtu.be/uWzVEGW2i5c.

Fifteen right-handed subjects (average age 25.9, 11 males, 4 females) participated in the study. Four of them had previous experience with haptic interfaces. Each subject spent about ten minutes practicing the control of the telemanipulation system before starting the experiment. Participants were briefed about all the tasks and afterwards signed an informed consent, including the declaration of having no conflict of interest. All of them were able to give the consent autonomously. The study was done in accordance with the 1964 Helsinki declaration and its later amendments.

B. Experimental conditions

We considered two different levels of human involvement in the control of manipulator A (teleoperation vs. shared control, T vs. S), two haptic feedback modalities (haptic feedback about the distance from the second robot and workspace/joints constraints vs. no haptic feedback, H vs. H), and two control policies for manipulator B (reactive vs. non-reactive, R vs. ✓ R).

1) Human involvement (T vs. S):

In condition T, subjects are able to control all the 6 DoF of manipulator A through the grounded master interface (see Sec. II-C1). Conversely, in condition S, the orientation of the gripper is controlled by an autonomous algorithm. When the gripper is close to the object to grasp, it is automatically oriented towards it; otherwise, the gripper is oriented to stay as much as possible away from the system constraints (see Sec. II-C2).

2) Haptic feedback (H vs. H): When haptic feedback is activated (H), subjects receive haptic stimuli about the feasibility of their commands against system constraints, such as joint limits, singular configurations, and collisions with manipulator B (see Sec. II-D). Conversely, in condition H, subjects do not receive any haptic feedback.

3) Control of the fully-autonomous manipulator (R vs. ✓ R): Manipulator B always performs an autonomous visual task to keep an object of interest visible. In condition R, it also uses the null-space of the above-mentioned primary visual task to minimize the constraint cost function H(q A , q B ) (see Sec. II-C). Notably, this secondary action will enable manipulator B to keep a safe distance from manipulator A, re-actively moving away when the latter comes too close. Conversely, in ✓ R, manipulator B only focuses on the primary visual task, indifferent to what the other manipulator does.

Considering all the possible combinations, we ended up with eight different experimental conditions: THR, TH ✓ R, T HR, T H ✓ R, SHR, SH ✓ R, S HR, S H ✓ R. For brevity, from now on we will omit the H and ✓ R variables (e.g., S H ✓ R is called S).

The cost functions introduced in Sec. II include several parameters, which let us control the curvature, rate of increase, and proximity to the limits of these functions. The choice of these parameters is challenging, system-dependent, and very delicate, as it must ensure a smooth force feedback. To choose the right parameters for our system and target application, we asked 2 expert operators to repeatedly carry out the task, changing the parameters at runtime (i.e., ρ X , α X , β X , k A1 , k A2 , k B1 , k B2 , k M ) to make the teleoperation as intuitive, safe, and comfortable as possible. Finally, we asked them to find a consensus on the parameters' values and we used those in our implementation.

C. Results

Average task success rate across conditions was 92.4 ± 6.3%. A Friedman test showed no statistically significant difference between the means of the eight feedback conditions. Figure 5a shows the total number of collisions occurred during the experiment between the two manipulators.

To compare other metrics, we ran three-way repeatedmeasures ANOVA tests on the data. Human involvement in the control (shared control vs. teleoperation, S vs. T), presence of haptic feedback (haptic feedback vs. no haptic feedback, H vs. H), and behavior of the fully-autonomous robot (reactive vs. non-reactive, R vs. ✓ R) were treated as within-subject factors. All data passed the Shapiro-Wilk normality test. Sphericity was assumed for all variables, since they all have two levels of repeated measures. Interaction effects between the factors were not statistically significant. Figure 5b shows the completion time. The three-way ANOVA test revealed a statistically significant change in the metric for the human involvement in the control (F(1, 14) = 52.165, p < 0.001, shared control was better) and the control behavior for manipulator B (F(1, 14) = 6.400, p = 0.024, reactive mode was better) variables. Figure 5c shows the linear motion covered by the robotic gripper during the task. The three-way ANOVA test revealed a statistically significant change in the metric for the human involvement in the control (F(1, 14) = 13.599, p = 0.002, shared control was better) and the control behavior for the second robot (F(1, 14) = 6.567, p = 0.023, reactive mode was better) variables. Immediately after the experiment, subjects were also asked to report the effectiveness of each feedback condition in completing the given task using bipolar Likerttype twenty-two-point scales. Figure 5d shows the perceived effectiveness for the eight experimental conditions. The threeway ANOVA test revealed a statistically significant change in the metric for the human involvement in the control (F(1, 14) = 34.700, p < 0.001), the presence of haptic feedback (F(1, 14) = 33.217, p < 0.001, shared control was better), and the control behavior for the second robot (F(1, 14) = 25.305, p < 0.001, reactive mode was better) variables. Finally, all fifteen subjects found conditions using the shared-control approach to be the most effective at completing the grasping task. Ten subjects chose SHR as the most effective, three SH, and two SR.

V. EXPERIMENT IN A REAL ENVIRONMENT

We also carried out an experiment in the real environment. The setup and task are similar to the simulated scenario of Sec. IV. The setup is shown in Fig. 1. Five right-handed subjects participated in the study. Three of them had previous experience with haptic interfaces. In this real-world case, we implemented control modalities TR, THR, SR, and SHR, which led to no collisions in Sec. IV. Of course, in this real-world case, we could not consider conditions which may lead to a collision between the two manipulators. All subjects successfully completed the task in all conditions and no collisions occurred. All subjects chose SHR to be the most intuitive and effective condition. A video of the experiment in the SHR condition (shared control, haptic feedback, and reactive control of manipulator B) is available as supplemental material and at https://youtu.be/uWzVEGW2i5c. A sequence of frames taken from the video is shown in Fig. 6, where it is possible to recognize all the relevant phases of the grasping task. In snapshots (a)-(b), it is possible to see the manipulator B reactively avoiding collisions while keeping track of the desired object. In snapshots (c)-(d), the shared-control algorithm is active and automatically orients the gripper of manipulator A toward the object to grasp. From the bottom graphs, we can see that the first half of the experiment is characterized by persistent haptic forces and by a significant reactive velocity of manipulator B. In the second half, the risk of collisions is lower and the operator can safely approach the object, aided by the action of the shared-control algorithm.

VI. DISCUSSION AND CONCLUSION

We presented a shared-control framework with haptic feedback, which enables the intuitive and effective teleoperation of a dual-arm robotic system. We tested the proposed architecture in simulated and real environments composed of two 6-DoF robotic manipulators. The human operator controls one arm (manipulator A), equipped with a gripper, through a 6-DoF grounded haptic interface. Haptic feedback provides the operator with information about joint and workspace limits as well as about the presence of singularity configurations and imminent collisions. The shared-control algorithm autonomously controls 2-DoF of the robotic manipulator, orienting the gripper toward the object to grasp. The other robotic arm (manipulator B) is equipped with a camera, and it moves autonomously to track a second object, placed near the one to grasp.

To demonstrate the feasibility and effectiveness of our methods, we conducted a human subject experiment in the simulated scenario. Fifteen subjects controlled manipulator A to grasp an object placed on a conveyor belt. We tested eight experimental conditions, considering two levels of human involvement in the control (shared control vs. classic teleoperation, S vs. T), two feedback modalities (haptic feedback about imminent collisions and workspace/joints constraints vs. no haptic feedback, H vs. H), and two control policies for manipulator B (reactively moving away from the other manipulator vs. non-reactive/still, R vs. ✓ R). Results proved the effectiveness and viability of our haptic-enabled sharedcontrol approaches. Using shared control (S) on manipulator A and the reactive mode (R) on manipulator B significantly improved the performance in most metrics (completion time, linear motion, perceived effectiveness). Conditions employing shared control were also the most preferred, confirming the all-round viability of such approach. Moreover, as expected, in conditions R, manipulator B was always able to prevent collisions with the other robot by moving away when the latter was approaching. Nonetheless, even in conditions ✓ R, repulsive haptic feedback (H) provided when the robots were too close showed good results (only two collisions happened in conditions TH, SH). This result is very promising, as haptic feedback acts only at the master side, leaving the action of manipulator B unaffected. In this respect, it is interesting to notice that, although the applied force f m should go to infinite as the distance between the two robots goes to zero (see eq. ( 5)), we still experienced two collisions in conditions H. This is due to the limited actuation capabilities of our haptic interface, which is obviously not able to provide arbitrarily high forces.

Finally, the experiment in the real scenario confirmed the results obtained with the simulated setup. Haptic guidance effectively steered the user toward the safe zones of the workspace, the reactive behavior enabled a safe interaction between manipulators, and the shared control made the task fast, easy, and intuitive to complete.

The proposed shared-control architecture is independent from the distribution of DoF between the autonomous controller and the human operator. While having the autonomous controller keep the gripper oriented toward the object seemed a natural choice for this work, the underlying shared-control architecture is independent from the task, the primary task variables, and the level of human involvement. For the future, we are planning to study how the number of DoF controlled autonomously affects the task performance, e.g., a system could use a highly-autonomous shared-control approach (i.e., many DoF managed by the autonomous controller) when it is operated by novices, while it could implement a lowlyautonomous shared-control approach (i.e., few DoF managed by the autonomous controller) when it is operated by experts. This flexible approach could be also useful when teaching new operators, employing different levels of autonomy according to the operator's experience. As future work, we will also investigate the effects of employing different cost functions and/or control gains on task performance. Finally, we plan to carry out more human subject experiments in real scenarios.

APPENDIX

We wish to prove that, in the shared-control modality when far from the object to be grasped, the slave angular velocity ω A in ( 13) is able to minimize the constraint cost function H(q A , q B ) in free-motion (v m = 0). This can be easily shown as follows: by considering variations w.r.t. q A (q B is taken care of by [START_REF] Lebel | An anticipative kinematic limitation avoidance algorithm for collaborative robots: Threedimensional case[END_REF]), one has Ḣ = ∂ T H ∂q A qA = ∂ T H ∂q A J † A u A . By applying [START_REF] Dube | A model of the humanoid body for self collision detection based on elliptical capsules[END_REF] when v m = 0, one has

Ḣ = -k A2 ∂ T H ∂q A J † A Z(J † A ) T ∂H ∂q A ,
which is a negative semi-definite quadratic form, therefore resulting in Ḣ ≤ 0 as claimed.

Fig. 1 .

 1 Fig.1. System architecture and main reference frames attached to the gripper and the object to be grasped.

Fig. 2 .

 2 Fig. 2. Minimum distance between two line sphere-swept bounding volumes (a) and distances among discrete sphere-swept bounding volumes (b). Solid (c), convex (d), and discrete-sphere (e) collision models of the robot.

Fig. 3 .

 3 Fig. 3. A visualization of the four null-space directions defined in (12).

Fig. 4 .

 4 Fig. 4. Experimental setup used for user study: slave side composed of simulated robots and vision; master side comprising the haptic interface.

  THR S SH SR SHR (b) Completion time (s).

  THR S SH SR SHR (c) Linear trajectory (m).

20 T

 20 TH TR THR S SH SR SHR (d) Perceived effectiveness.

Fig. 5 .

 5 Fig. 5. Experimental evaluation. Mean and 95% confidence interval of (a) total number of collisions, (b) completion time, (c) linear motion, and (d) perceived effectiveness of the eight feedback conditions are plotted.

Fig. 6 .

 6 Fig. 6. Data time history recorded during the real experiments. Upper graph: minimum distance d min and collision cost value H; bottom graph: haptic guidance fm and escaping velocity v B norms.

In this case, we assume that the motion of the master device is blocked (via software) along the two non-selected components.

We note that the presented analysis does not account for the controller switching in the shared-control modality when close/far from the target object. If this switch needs to be taken into account because of non-negligible effects on the total energy, one could employ the energy tank machinery for passifying potential instabilities due to the switching mechanism[START_REF] Secchi | Bilateral teleoperation of a group of uavs with communication delays and switching topology[END_REF]. However, we empirically found the switching to have a negligible effect on the system stability during our experiments.

supported by the EU H2020 RoMaNS project 645582 and FP7 RoDyMan project 320992.


L R Douat

I Queinnec

G Garcia

S Tarbouriech

Identification and Vibration Attenuation for the Parallel Robot Par2

Luiz R. Douat, Isabelle Queinnec, Germain Garcia, Micaël Michelin, Franc ¸ois Pierrot and Sophie Tarbouriech

Abstract-Par2 is a two-degree-of-freedom parallel robot designed for high-speed and high-accuracy industrial pick-andplace operation tasks. As a result of the high acceleration trajectories, the end-effector undergoes some undesirable vibrations after reaching the stop positions, compromising its precision and leading therefore to an increase in the operation cycle time. Accelerometer sensors placed on the end-effector and piezoelectric patch actuators wrapped around the robot arms are employed in order to actively reduce these vibrations in a non-collocated closed-loop setting. After submitting the robot to an identification procedure, the obtained nominal model is used to synthesize a reduced order controller with the H ∞ Loop Shaping technique. Performance analysis as well as simulation and experimental results exhibit that vibration reduction is achieved around the nominal operating point, but fails for some extreme operating points, due to high control efforts. An antiwindup strategy is then employed to deal with the saturation of the actuator, which allows to achieve vibration attenuation on the whole operation domain, for a given configuration of the robot at the stop point.

Index Terms-High-speed parallel robotics, robust control, vibration attenuation, piezoelectric actuators, saturation, antiwindup.

I. INTRODUCTION

In the context of industrial handling and assembly tasks, the increasing demand for productivity imposes more and more strict requirements in terms of production cycle times. Lighter than their serial counterparts, parallel robots [START_REF] Merlet | Parallel Robots[END_REF] are able to accomplish the desired trajectories faster and with a better precision. For pick-and-place operations, these higher trajectory accelerations demand correspondingly higher decelerations each time the end-effector reaches some stop position, giving rise to undesirable vibrations and therefore impairing the positioning accuracy. The gain obtained in reducing the trajectory time by using higher accelerations is thus countered by an increase in the residual vibration levels after reaching the stop positions.

For minimizing these vibrations, different measures could be undertaken:

1) Improvement of the robot design (light but stiff robot links) [START_REF] Company | Modelling and preliminary design issues of a 3-axis parallel machine-tool[END_REF], [START_REF] Marquet | Enhancing parallel robots accuracy with redundant sensors[END_REF];

2) Optimization of the pick-and-place trajectory, avoiding inflection points in the higher successive position derivatives [START_REF] Constantinescu | Smooth and time optimal trajectory planing for industrial manipulators along specified paths[END_REF], [START_REF] Gasparetto | Experimental validation and comparative analysis of optimal time-jerk algorithms for trajectory planning[END_REF]; 3) Active vibration control during the trajectory by means of the motor magnetic moments [START_REF] Da Silva | Integrating structural and input design of a 2-DOF high-speed parallel manipulator: a flexible model-based approach[END_REF], [START_REF] De Luca | Handbook of Robotics[END_REF]; 4) Active vibration control during the trajectory by means of piezoelectric actuators [START_REF] Gurses | Vibration control of a singlelink flexible manipulator using an array of fiber optic curvature sensors and pzt actuators[END_REF], [START_REF] Hassan | Active vibration control of a flexible one-link manipulator using a multivariable predictive controller[END_REF]. This case is of particular interest when the vibrations may induce hazardous consequences to the robot during the trajectory; 5) Active vibration control at the stop positions by means of piezoelectric actuators [START_REF] Jalili | Piezoelectric-based vibration control: from macro to micro/nano scale systems[END_REF], [START_REF] Moheimani | Piezoelectric transducers for vibration control and damping[END_REF], [START_REF] Preumont | Vibration control of active structures -an introduction[END_REF]. According to the type and placement of the sensors and actuators in the robot, two different architectures can be characterized: collocated and non-collocated. In the non-collocated case [START_REF] Qiu | Active vibration control of a flexible beam using a non-collocated acceleration sensor and piezoelectric patch actuator[END_REF], [START_REF] Tliba | Active vibration damping of a smart flexible structure using piezoelectric transducers: H-infinity design and experimental results[END_REF], addressed in this paper, different kinds of sensors can also be used, as for instance accelerometers [START_REF] Algermissen | Applied robust control for vibration suppression in parallel robots[END_REF], generally allowing better performances [START_REF] Verhoeven | Robust control of flexible motion systems: A literature study[END_REF]. Par2, depicted in Figure 1, is one such fast pick-andplace two-degree-of-freedom parallel robot [START_REF] Baradat | Par2: a spatial mechanism for fast planar, 2-DOF, pick-and-place applications[END_REF]. Despite its very good stiffness properties, there are still some remaining vibrations showing up at the end of each trajectory, which are measured by means of three orthogonally oriented accelerometers (in X, Y and Z directions) placed on the end-effector. Given that these vibrations hinder the accurate positioning of the end-effector at the stop positions, it is at these precise points that active control is required. Following the strategy adopted in [START_REF] Algermissen | Applied robust control for vibration suppression in parallel robots[END_REF], piezo-actuator patches are glued all along and around the arms of Par2 in order to attenuate the residual vibrations by means of the controlled deformation of the arms (section II). The control objective is then to reject disturbances, seen as the residual vibrations coming up just after the arrival of the end-effector at a stop position, for various trajectories and load conditions.

In order to synthesize a control law, it is essential to have first a representative model of the plant. Piezoelectric devices generally exhibit nonlinear response behavior involving hysteresis, drift, creep, and/or time delay [START_REF] Leang | Feedback-linearized inverse feedforward for creep, hysteresis, and vibration compensation in AFM piezoactuators[END_REF], which may be critical in certain applications of high-precision positioning of nano-devices. In the current work, however, we are mainly concerned with the vibration response of the piezo-elements further used to compensate for the vibrations induced by the motion of the robot. Set as a black box connecting the piezoelectric voltage (system input) to the measured endeffector acceleration (system output), a system identification procedure is performed using state subspace techniques [START_REF] Ljung | System identification -Theory for the user[END_REF], leading to the production of a nominal model of the system [START_REF] Douat | Flexible model identification of the parallel robot Par2[END_REF] (section III).

Aiming at synthesizing a robust controller, capable of attenuating the vibrations generated at different operating points (trajectories and load conditions), the H ∞ Loop Shaping strategy is employed (section IV). It differs from previously published papers from the authors on the subject, where a mixed sensitivity H ∞ design strategy was considered [START_REF] Douat | H ∞ control applied to the vibration minimization of the parallel robot Par2[END_REF], [START_REF] Douat | Piezoactuated vibration attenuation of the parallel robot Par2[END_REF]. The current strategy is adopted here to benefit from its easy implementation with robust and easily tunable filters.

For some operating points kept away from the nominal point, saturation of the actuator may however induce a nasty effect, from loss of performance to loss of stability. Then, an anti-windup compensator strategy may be addressed as an additional control loop on the closed-loop system, active only when saturation occurs. Many different anti-windup design strategies have been addressed in the literature [START_REF] Tarbouriech | Stability and stabilization of linear systems with saturating actuators[END_REF], [START_REF] Zaccarian | Modern anti-windup synthesis[END_REF]. We consider here a direct linear anti-windup compensator design, and we restrict to the case of a simple static anti-windup gain which acts on the output of the robust controller previously designed (section V).

At each step of the whole control design procedure, simulations and experimental results are proposed to exhibit the potential interest of the approach proposed in this paper.

II. PROBLEM STATEMENT

The manipulator end-effector trajectories (dash-dotted representation in Figure 2), constrained to the xz plane between stop positions Pos1 (x = 0.35 m, z = -0.925 m) and Pos2 (x = -0.35 m, z = -0.925m), are commanded by the so-called active arms (those directly connected to the motors, Figure 1) while the passive ones hinder, as much as possible, the perpendicular motions.

Pos1 Pos2 X Z Fig. 2. Stop positions
The piezo-actuator patches are glued on those arms (Figure 3), and are labelled with letters A and P, given accordingly to the active or passive characteristic of the arm. They are PZT elements (PI, DuraAct P-876.A12) whose deformation can be considered linear with respect to an input voltage signal ranging from 0 to 400 V . Each piezo-actuator input voltage signal is supplied by means of a high-voltage amplifier (Walcher Electronik, HVA-400-100-ISO) with bandwidth of 8 kHz, that amplifies the -10 V to 10 V control signal sent by the control unit (Adept, sMI6). Despite presenting vibrations during the whole pick-andplace trajectory, their effects are problematic only when amplified by the high decelerations on the arrival at a stop position, overly exciting some modes of the robot arms and translating into loss of end-effector precision. These amplified vibrations at the stop positions are called the disturbances to our system. The accelerometers (BrüelKjaer, DeltaTron 4506) employed to measure these vibrations are able to provide linear measures up to 700 m/s 2 (70 g) with a sensitivity of 100 mV /g and frequencies between 0.3 Hz and 3 kHz. Before sending them to the control unit, the accelerometer signals are filtered by a charge amplifier (Endevco, Model 133) of 100 kHz bandwidth and adjustable cut-off frequency.

The pick-and-place motions of the robot (from Pos1 to Pos2 and vice versa) can be accomplished for many different operating points. Each operating point is characterized by a load on the end-effector (varying from empty to 3.5 kg) and a particular trajectory peak acceleration (varying from 10 g to 30 g). Considering that the robot motor currents are limited to 30 A, the higher is the load on the end-effector, the lower is the maximum acceleration allowed for a trajectory. The operating points taken into account in this work, for which the residual vibrations have to be reduced, are depicted in Figure 4, where the upper circles represent the maximal pick acceleration for each load condition. We consider in this work the case where the piezo control is activated immediately after arriving at the stop position Pos2 (the residual vibrations observed at the stop positions Pos1 and Pos2 are symmetrical and, thus, only Pos2 is considered here). For the control law synthesis, we first need to identify a nominal open-loop dynamical system which connects the piezo-actuator to the accelerometer. Such a model is associated to a particular nominal operating point. We then want to synthesize a robust control law, that is, a controller able to reduce the vibrations not only for the nominal operating point, but also for the largest set of uncertain model representations around this point. We choose a 2.5 kg load with a 18 g trajectory to be our nominal operating point and want to verify the efficacy of its associated robust controller over the whole space of selected operating points (Figure 4).

It was experimentally verified that, due to a much lower damping factor, the residual vibrations in all desired operating points are much more substantial in the Y direction. This phenomenon is illustrated for the nominal operating point in Figure 5, in the frequency domain, and in Figure 6, in the time domain. Two experiments are plotted, dist 1 and dist 2, to illustrate that they are reproducible. Figure 5 exhibits a low damped flexible mode around 21.7 Hz in the Y direction. In Figure 6, the disturbances in the Y direction between 0.01 s and 1 s are far less damped than the disturbances in the other two directions. To control the acceleration in the Y direction, it is rec-ommended to manipulate the piezo-element patches on the passive arms (labelled with P), which are exactly provided to limit those off-line vibrations. Then in order to remain as simple as possible, rather than to consider a multiple input system with one input associated with each passive piezoactuated arm, we prefer to select a pair of passive piezoactuated arms (namely P1P2 in Figure 1), jointly and equally actuating and from now on considered as a single element, to be our system input. We then consider a SISO control problem and, moreover, we keep free the other piezo-actuators for a potential future extension of the approach to the full MIMO control problem. According to the system representation of Figure 7, this work can be initially divided into 3 steps. The first step, briefly described above and presented in detail in [START_REF] Douat | Flexible model identification of the parallel robot Par2[END_REF], is the characterization of the frequency spectra of the disturbance signal d for the nominal operating point. The second step concerns the open-loop identification of the nominal model G from the piezo-actuator input voltage u to the accelerometer in the Y direction (section III). The third set addresses the design problem of a robust controller K capable of attenuating the vibrations for the nominal operating point as well as for all operating points represented in Figure 4 (section IV). It is initially done without taking explicitly into account the saturation of the voltage delivered to the piezo-actuator. Then, the influence of this nonlinear element in the closed-loop may be further studied (section V).

III. IDENTIFICATION Dynamical model production can be either analytic or blackbox based. An analytic model is obtained based on physical principles and allows a better comprehension of what is

happening inside the model. Its main disadvantages are the often strenuous conceptual phase and the need of further experimental parameter validation (data updating). If the dynamical model is to be used only for control purposes, a fast-to-obtain and easy-to-validate alternative is to employ black-box system identification techniques. These identification techniques deal with the problem of building mathematical models of dynamical systems based only on its input-output data [START_REF] Ljung | System identification -Theory for the user[END_REF]. The validation of these models is often done by the comparison in the frequency domain of the estimated model with a nonparametric model obtained from the same input-output data. A typically used non-parametric model in the frequency domain is called ETFE (Empirical Transfer Function Estimate), given by the division of the discrete Fourier transforms of the output and input vector signals.

We can classify the identification techniques into classical methods and methods based on subspaces. The works of [START_REF] Ljung | System identification -Theory for the user[END_REF] and [START_REF] Verhaegen | Filtering and system identification: An introduction[END_REF] show that subspace identification methods allow an easy extraction of good quality models. For that reason, the open-loop models relating the piezo-actuator input voltage to the accelerometer output of Par2 are identified based on one of these subspace methods. The identification procedure is performed twice, giving rise to two models: one used for control synthesis and another for validation of the obtained controller in computer simulations.

A. Subspace identification methods

Subspace identification emerged from the association of the multivariable statistical analyses methods (stochastic realization) from the 30's and modern numerical linear algebra tools like QR and SVD decompositions [START_REF] Katayama | Subspace methods for system identification[END_REF], [START_REF] Larimore | System identification, reduced-order filtering and modeling via canonical variate analysis[END_REF]. Contrarily to classical parameter identification methods as ARMAX and Box-Jenkins [START_REF] Ljung | System identification -Theory for the user[END_REF], subspace methods allow to obtain a statespace system representation based on the estimation of certain subspaces from matrices formed with the input and output experimental data vectors. Some advantages from subspace methods over the classical ones are:

• Only the desired (expected) system order for the statespace model must be provided as an entry parameter. This unique-parameter feature makes subspace methods easier to employ than their counterpart methods (classical parametric methods of equivalent complexity, where the polynomial orders have to be specified, usually require from 3 to 5 tuning parameters); • There are no nonlinear optimization problems. Subspace methods are based on subspace matrices projections without any parametric optimization procedure involved, avoiding problems like local-minima and non-convexity; • The method is directly expandable to multivariable systems and without the drawbacks of classical methods applied to MIMO systems which often give ill-conditioned models and are usually very time demanding [START_REF] Chiuso | Some algorithmic aspects of subspace identification with inputs[END_REF]. There exist many different subspace identification algorithms: MOESP [START_REF] Verhaegen | Identification of the deterministic part of MIMO state space models given in innovations form from input-output data[END_REF], IVM [START_REF] Viberg | On subspace-based methods for the identification of linear time-invariant systems[END_REF], CVA [START_REF] Larimore | Canonical variate analysis in identification, filtering and adaptive control[END_REF], N4SID [START_REF] Van Overschee | N4SID: Subspace algorithms for the identification of combined deterministic-stochastic systems[END_REF], etc. In this work the N4SID (Numerical algorithm For Subspace IDentification) algorithm was employed because of its proven applicability to the identification of vibration system models [START_REF] Sumiya | Active damping using N4SID model referenced predictive anti-phase pitch driving method[END_REF], [START_REF] Nestorovic-Trajkov | Active control of a piezoelectric funnel-shaped structure based on subspace identification[END_REF], [START_REF] Douat | Flexible model identification of the parallel robot Par2[END_REF].

In the N4SID algorithm, an extended observability matrix of a system with state-space representation in some base T is estimated through the input-output data vectors. This observability matrix allows the extraction of a Kalman filter state sequence which is then used to obtain, by means of the solution of a least square problem, an estimate of the statespace matrices of the original system in the same base T .

B. Open-loop piezo-actuator identification

For the identification procedure, a chosen input voltage signal excited the resonant modes of the piezo-actuator P1P2. The corresponding time response vibration was then stored for the accelerometer measurements in the Y direction. These input and output data vectors were pre-processed and used to build a mathematical description of the piezo-actuator dynamical behaviour, by means of the N4SID subspace identification algorithm. The identified model of the piezo-actuator did not depend on the peak acceleration as the identification experiments were performed with the end-effector initially at rest at the position Pos2. On the other hand, it depended on the load at the end-effector. We considered here the end-effector with a nominal load of 2.5 kg. The experimental identification set-up is summarized in table I. During the identification steps, a sampling rate of 1 kHz was chosen after verification that the initially adopted 2 kHz sampling rate led to stable but ill conditioned models (poles very close to each other near the real value 1 in the unit circle). According to the identification experiments performed in [START_REF] Algermissen | Applied robust control for vibration suppression in parallel robots[END_REF] and [START_REF] Douat | Flexible model identification of the parallel robot Par2[END_REF], a chirp signal was used as input excitation. It was experimentally verified that the highest main frequency modes excited in the Y direction for all considered operating points (Figure 4) did not surpass by far the limit of 30 Hz. For that reason the chirp signal used was empirically chosen to sweep linearly a five times larger spectrum (from 0 up to 150 Hz) during 100 s (sufficient to minimize the influence of the sweeping transitory effects [START_REF] Ljung | System identification -Theory for the user[END_REF]). The output acceleration signals were analogically low-pass filtered by the charge amplifier with a cut-off frequency set to 500 Hz (Nyquist frequency) before being stored. Despite the analogical filtering process, the noise levels verified in the output data were still considerable. Therefore, a digital 5th order low-pass Butterworth filter with 150 Hz cut-off frequency was applied to the input-output stored data, shown in Figure 8. The procedure was applied twice to obtain two input-output data vectors. It has been checked that associate non-parametric frequency models show a good reproducibility between the two experiments. They were then used to build two models, initially named mod ini ctrl 2.5kg and mod ini sim 2.5kg. These models were further reduced and became mod ctrl 2.5kg, used for control synthesis purposes, and mod sim 2.5kg, employed for validating the obtained controller in numerical simulations.

Let us consider the first set of input-output data. By using the N4SID algorithm with different desired state-space model orders, the obtained candidate mod ini ctrl 2.5kg models can be seen in Figure 9 with also the non-parametric model directly issued from the experimental data (NP model). Both amplitude and phase are represented in this figure. Nevertheless, the control technique adopted in this work is based on a worst case amplitude criterion without any consideration to the phase. This reason motivated our choice to present, from now on, only the harmonic gain amplitudes. Figure 9 shows that the 40 order N4SID model was close to the non-parametric model (used for validation) in the frequencies of interest, being therefore the selected one. Considering the high order of the model, a reduction procedure was then performed in two steps, through techniques that preserve the main dynamical properties of the original model. Note that it is recommended to select a high-order model, then to reduce it, as an original low-order model would not focus enough on the low frequency part of the data, and would put too much emphasis to the high frequency part related to noise activity. In the first step of reduction, the high frequency modes were suppressed. Under the assumption of linearity, all the poles and zeros beyond the input signal excitation range were considered to be noisy measurements and did not reveal any information about the system. These singularities from 150 Hz to 500 Hz were thus eliminated and the model was reduced to a 22 order model. The second step consisted in rewriting the resulting system in a new base, called balanced, where the controllability and observability gramians were equal and diagonal. This realization allowed to isolate the states that less contributed to the output system dynamics. Once isolated, these states could be eliminated with, for instance, the singular perturbation model reduction method [START_REF] Gu | Robust Control Design with MAT LAB R[END_REF].

After writing the 22 order model in a balanced realization, a Hankel singular value analysis shows, in Figure 10, that the first 14 Hankel singular values correspond to 96% of the model energy. Using the singular perturbation reduction method it is therefore reduced to a 14 order model, being now called mod ctrl 2.5kg. The same procedure of identification and reduction was applied also to mod ini sim 2.5kg (respectively mod sim 2.5kg). Considering that the model for simulation shall be more informative than the model for control synthesis, the N4SID method was applied to mod ini sim 2.5kg using a conservative order of 80. The high frequency modes were suppressed and the model order was reduced to 43. After a balanced realization followed by a singular perturbation reduction method the final mod sim 2.5kg model obtained was of order 30 (corresponding to 98% of the energy of the partially reduced 43 order mod ini sim 2.5kg model).

The non-parametric frequency model, compared to both mod ctrl 2.5kg (for control synthesis) and mod sim 2.5kg (for controller validation in numerical simulations) models is shown in Figure 11. The models mod ctrl 2.5kg and mod sim 2.5kg were identified in the discrete time domain. However, the MAT LAB R function (hinfsyn from the Robust Control Toolbox) used further on to synthesize the H ∞ control law is available only for the continuous time case. This reason together with the fact that the visual analysis of poles and zeros is often easier in the continuous time domain prompted us to convert the discrete time models into the continuous time. These conversions were done by using a Tustin bilinear transformation [START_REF] Oppenheim | Signals and Systems[END_REF]. Figure 12, comparing the model mod ctrl 2.5kg in the discrete and continuous time, shows that the approximation is valid. [START_REF] Mcfarlane | A loop shaping design procedure using H ∞ synthesis[END_REF], [START_REF] Zhou | Essentials of Robust Control[END_REF]. The control performance goals are usually to reduce the magnitude of the sensitivity function S in some desired (normally low) frequency ranges, while preventing high control signals in high frequencies, which could lead to noise amplification, saturation and, in the particular case of resonant systems, the spillover phenomenon (excitation of non-modeled high frequency resonant modes).

Based on this objectives and taking into account Figure 13, the H ∞ LSDP can be described in four steps: ). The filters are chosen accordingly to the previously mentioned principles for conforming the desired S and KS closed-loop functions; 2) After a normalized coprime factorization of the obtained model G S , find the robust controller K S that stabilizes the closed-loop system for the largest possible set of uncertainties (associated to a stability margin ε) acting on the coprime factors [START_REF] Mcfarlane | A loop shaping design procedure using H ∞ synthesis[END_REF]. This is done by solving the H ∞ problem:

G W 1 W 2 K S G S
min K W -1 1 W 2 G (I -KG) -1 W 1 GW -1 2 ∞
3) If the resulting stability margin ε is greater than 0.2, the frequency response of K S W 2 GW 1 is close to that of W 2 GW 1 [START_REF] Mcfarlane | A loop shaping design procedure using H ∞ synthesis[END_REF], and we can conclude that our shaped function G S is close to the open-loop gain K S W 2 GW 1 ; 4) In this case, the controller K of the original system can be recovered as

K = W 1 K S W 2 .

B. Controller design and experimental validation

As we are dealing with a SISO system, only the W 1 filter is considered for the G S shaping in the LSDP (the W 2 filter can be set to 1). The G S function is shaped to favour both high gains at low frequencies, implying a reduced S at these frequencies, and a decreasing slope at the higher frequencies, to prevent high amplitude signals and the spillover effect. The W 1 filter is chosen as

W 1 = 0.055(s + 300)(s + 1600) 2 (s + 100) 3 (1)
and its effect on the nominal mod ctrl 2.5kg model is shown in Figure 14. Applying the LSDP, a controller K A of order 19 is obtained. This controller is further reduced to a 4 order controller, named K Ar , by using a balanced realization followed by a singular perturbation reduction method. It is defined by

K Ar = -0.44(s 2 + 9s + 1.32e4)(s 2 -232s + 4.2e5) (s 2 + 159s + 6987)(s 2 + 29.6s + 1.48e4) (2)
The comparison between K A and K Ar is illustrated in Figure 15, showing a good correspondence in the desired frequency range.

The sensitivity function S Ar for the nominal simulation model (mod sim 2.5kg) with the reduced order controller K Ar is given in Figure 16, showing that frequencies around 21.7 Hz (corresponding to the main Y direction disturbance mode shown in Figure 5) should be well attenuated. We consider now the experimental validation of the controller K Ar applied to all operating points of Figure 4. The performance criterium adopted to measure attenuation is based on the RMS acceleration value [START_REF] Balas | Robust control of flexible structures: theory and experiments[END_REF]. The percentage reduction From Figure 20 we can observe that the closed-loop responses with the controller K Ar are better attenuated than those in open-loop for all but two cases. These exceptional limit cases correspond, first, to the system with a 3.5 kg load and a trajectory of 15 g (representing a 3.6% closed-loop vibration augmentation), and second, to all situations without load. This second limit case is not even represented because the controller brings the system to instability for all trajectories. The two limit cases have their time responses respectively represented in Figures 21 and22. We can observe in both figures the important saturation levels leading, in the first limit case, to performance degradation and, in the second limit case, to system instabilities (illustrated for 10 g, 20 g and 30 g). The experimental results obtained with the controller K Ar have brought out that, despite attenuating the disturbances for almost all chosen operating points, there are still some remaining limit cases where this controller is not efficient. The stability and performance problems for these limit cases are associated with the high control signals delivered by the controller, which induce saturation of the actuator. This can also be exhibited thanks to a stability analysis of the closed-loop saturated system in presence of exogenous energy-bounded disturbance (see [START_REF] Tarbouriech | Stability and stabilization of linear systems with saturating actuators[END_REF], Chapter 3, Proposition 3.14), considering explicitly that the input of the piezo-actuator model is subjected to amplitude saturation -u 0 ≤ u(t) ≤ u 0 , u 0 = 10V . For an a priori given bound on the disturbance, i.e., for a disturbance belonging to the set 1

W = w : [0, ∞) → ℜ q ; ∞ 0 w(τ) w(τ) dτ ≤ δ -1 1
x denotes the transpose of vector x.

one can evaluate the L 2 performance gain γ between the disturbance and the output acceleration of the end-effector. Results are summarized in Table II, both for the saturationfree case (second column) and for the saturated input case (third column).

We can check that, as expected, the controller K Ar , designed without taking saturations into account explicitly, ensures excellent results for the nominal model in the unconstrained case (γ=0.08). As the actuator is actually subject to saturations, the saturated case (third column) exhibits that, for the nominal case, vibrations should not be amplified (γ=1.00) but the expected performance is much reduced. When the load is changed (reduced to 2kg or increased to 3.5kg), the expected performance is degraded as well for the unconstrained case and for the saturated case. This is observed experimentally, i.e. the disturbance is not attenuated but amplified in the higher load case (Figure 21). Moreover, the controller K Ar is definitively not appropriate when the robot does not carry any load (γ=165.27 without limitation and the problem is even infeasible in presence of limitations of the actuator). This explains the loss of stability experimentally observed (Figure 22) in that case.

At that point, one option may be to design another controller with the LSDP such as to reduce the control signals and avoid, as much as possible, saturation of the piezoelectric actuator. This is for example the case of controller K Br , designed with the W 1 weighting filter selected such as to reduce the gain from 0.2 to 0.1. Such a solution however also reduces the disturbance attenuation performance around the nominal operation point (see Table II, column 5, where the nominal attenuation at 2.5kg is increased from 0.08 with K Ar to 0.23 with K Br ). This is also illustrated in Figure 23, which exhibits that the vibration attenuation with the second controller is increased with respect to the open-loop behavior for all the operating points of Figure 4, even for the limit cases, but is decreased at (and around) the nominal operating point for which have been designed the controllers. As an alternative, we propose to build an additional antiwindup compensator, which does not modify the performance of the nominal controller when no saturation occurs and only acts when the saturation becomes active. A similar anti-windup strategy has been already considered for active vibration isolation systems in a slightly different problem [START_REF] Teel | An anti-windup strategy for active vibration isolation systems[END_REF], and proved its efficiency.

The anti-windup strategy employed in this paper is fully described in [START_REF] Tarbouriech | Stability and stabilization of linear systems with saturating actuators[END_REF] (Chapter 7, Proposition 7.4) and briefly commented here. Let us go back to Figure 7 and consider now an anti-windup compensator block acting on the controller K Ar (s). The state-space dynamics controller K Ar (s), between the acceleration y of the end-effector (input of the controller) and the controller output y K is then modified as follows:

ẋK (t) = A K x K (t) + B K y(t) + v x (t) y K (t) = C K x K (t) + D K y(t) (3)
where the additional input v x (t) allows to introduce the controller's anti-windup correction factor in the dynamics of the controller. In a static strategy, it is defined as

v x (t) = D x aw (u(t) -y K (t)) (4)
i.e., it uses the error between the real input u(t) = sat(y K (t)) and the output of the controller y K (t) to compute the antiwindup action. Considering the piezo-actuator model identified with a 3.5 kg load, the nominal controller K Ar , the disturbance bound δ = 0.2 and the bound u 0 = 10, one computes the static anti-windup compensator

D x aw = [-0.3032 -0.1134 -0.0368 -0.1694] (5)
Roughly speaking, the associated L 2 performance gain γ is improved (Table II, column 4) with respect to the saturated case without anti-windup (Table II To go further, one can check that, with the same antiwindup gain, vibrations are also attenuated in the case of the robot without load (see Figure 25). Moreover, using this solution where the anti-windup action is only active when the input saturates, the vibration attenuation is unchanged for all the other operating points where the actuator almost never saturated.

VI. CONCLUSIONS AND FUTURE WORKS

This work presented the results obtained for the piezoactuated robust vibration attenuation of a high-speed and high-accuracy pick-and-place parallel robot in a set operating pointsdescribed by the peak accelerations and loads and for a given configuration of the stop position.

First, based on well reproducible experimental data, the N4SID subspace identification algorithm was used to estimate models for the piezo-actuator in a nominal operating point. A controller was synthesized by employing the H ∞ Loop Shaping Design Procedure to the nominal model. This controller was capable of attenuating the vibrations for almost all operating points. In the other operating points (limit cases), the high control signals saturated the system input, leading, as a consequence, to loss of performance and, in some cases, even to loss of stability. Then, rather than designing another less active controller able to prevent, as much as possible, saturating signals, we have added a static anti-windup compensator to alleviate the malicious effect of saturation. Fair performance level has been observed in such a case in experimental evaluations on the whole domain of operation.

The results obtained underline the feasibility of such strategy involving piezoelectric patches on the arms of the robot for vibration attenuation during high-speed pick-and-place movements of the parallel robot. They pave the way for future evaluation of the strategy from an industrial point-of-view, in a MIMO context and for different configurations of the stop point. Actually, if the robot configuration at the stop point is changed, then one can expect that the model will not be any more valid. Either one can design a controller for each possible configuration Pos2 (if the robot is used within a limited and a priori given set of configurations) or the control design has to be robust to the set of configurations.

Fig. 1 .

 1 Fig. 1. Parallel robot -Par2

Fig. 3 .

 3 Fig. 3. Par2s arms equipped with the piezoelectric patches

Fig. 4 .

 4 Fig. 4. Operating points -load conditions and trajectories

Fig. 5 .Fig. 6 .

 56 Fig. 5. Disturbances dist 1 and dist 2 -Nominal operating point -Frequency domain

1 Fig. 7 .

 17 Fig. 7. Function block representation of the closed-loop system

2 )Fig. 8 .

 28 Fig. 8. Data for mod ctrl 2.5kg -Input and output signals

Fig. 9 .

 9 Fig. 9. Frequency domain models obtained for mod ini ctrl 2.5kg with the N4SID algorithm (orders 10, 20, 30 and 40)

Fig. 10 .

 10 Fig. 10. Hankel singular values for the partially reduced model mod ini ctrl 2.5kg of order 22

Fig. 11 .

 11 Fig. 11. Models mod ctrl 2.5kg (order 14) and mod sim 2.5kg (order 30)

Fig. 12 .

 12 Fig. 12. Model mod ctrl 2.5kg in discrete and continuous time

Fig. 13 .

 13 Fig. 13. Loop Shaping Design Procedure (LSDP) scheme

mod sim 2 .Fig. 14 .

 214 Fig. 14. LSDP open-loop ponderation

Fig. 15 .Fig. 16 .

 1516 Fig. 15. Complete and reduced order controllers

  .

Fig. 17 .

 17 Fig.[START_REF] Larimore | System identification, reduced-order filtering and modeling via canonical variate analysis[END_REF]. Experimental results -Robot with 2.5 kg and reduced controller K Ar -18 g trajectory -Power spectra

Fig. 18 .Fig. 19 .

 1819 Fig. 18. Experimental results -Robot with 2.5 kg and reduced controller K Ar -18 g trajectory -Simulated and experimental responses

2 )Fig. 20 .

 220 Fig. 20. Experimental results -Percentual RMS acceleration reduction during the first 2 seconds -Closed-loop system with controller K Ar in relation to the open-loop system

Fig. 22 .

 22 Fig. 22. Experimental results -Robot without load and controller K Ar -10 g, 20 g and 30 g trajectories -Time response and control signal

2 )Fig. 23 .

 223 Fig. 23. Experimental results -Percentual RMS acceleration reduction during the first 2 seconds -Closed-loop system with the controller K Br in relation to the open-loop system

Fig. 24 .

 24 Fig. 24. Experimental results -Robot with 3.5 kg and controller K Ar + anti-windup -15 g trajectory -Time response and control signal

Fig. 25 .

 25 Fig. 25. Experimental results -Robot without load and controller K Ar + anti-windup -10 g, 20 g and 30 g trajectories -Time response and control signal

  Fig. 21. Experimental results -Robot with 3.5 kg and controller K Ar -15 g trajectory -Time response and control signal

	) Acc Y (m/s 2	-5 0 5								OL experiment CL experiment
				0	0.2	0.4	0.6	0.8	1	1.2	1.4	1.6	1.8	2
	Control signal (V )	0 -10 0 10	0.2	0.4	0.6	0.8	1	1.2	1.4	1.6	1.8	2
									Time (s)
	) Acc Y (m/s 2	-10 0 10 -20							10 g trajectory 20 g trajectory 30 g trajectory
				0	0.2	0.4	0.6	0.8	1	1.2	1.4	1.6	1.8	2
	Control signal (V )	0 -10 0 10	0.2	0.4	0.6	0.8	1	1.2	1.4	1.6	1.8	2
									Time (s)

TABLE II PERFORMANCE

 II ANALYSIS OF THE CLOSED-LOOP SYSTEM (FOR δ = 0.2)

			controller K Ar		controller K Br
	model	γ unconstrained	γ with saturation	γ with saturation and AW	γ unconstrained
	empty	165.27	unfeasible	unfeasible	77.72
	2kg	0.52	1.98	1.09	0.77
	2.5kg	0.08	1.00	0.75	0.23
	3.5kg	0.80	1.97	1.45	1.02

VII. ACKNOWLEDGMENTS

This work was supported in part by the French National Research Agency (ANR), in the framework of the Objectif 100G project. The authors also gratefully acknowledge the reviewers for their fruitful comments.