Rafael Balderas Hill

Sébastien Briot

Abdelhamid Chriette

Philippe Martinet

Exploiting Natural Dynamics in order to Increase the Feasible Static-Wrench Workspace of Robots

Keywords: static-wrench workspace, natural dynamics, shooting method

The reachable workspace of a robot carrying a payload is usually limited by the maximal value of the torque that each actuator can deliver. This results in limiting the zones for the robot to operate with the payload due to a possible division of its static-wrench workspace into several disconnected aspects. In order to increase the reachable workspace areas, this paper proposes to exploit the natural oscillations in dynamics, so that the robot can carry a payload which is out of its feasible static-wrench workspace, i.e. to perform motions between two disconnected aspects, while constraining the torques of the actuators. This is done thanks to the solution of a boundary value problem, which seeks to smartly exchange the gravity potential energy and the kinetic energy in order to connect two desired payload positions, which are placed in two disconnected aspects. Simulations of the suggested approach on a 2-degree-of-freedom robot are performed and show the efficiency of the proposed approach.

Introduction

It is well-known that in several industrial robotic processes, robot manipulators are required to have as large as possible workspace. Nevertheless, when for some specific tasks, the robots are required to carry heavy payload objects, the workspace can be limited by the actuator maximum torques. The remaining workspace, limited by the torque limits, is called feasible static-wrench workspace, and it is generally defined by the set of all feasible end-effector poses for which the platform with payload can be reached statically, i.e. taking into account only the gravity effects [START_REF] Gouttefarde | Wrench-feasible workspace of parallel cable-driven mechanisms[END_REF][2] [START_REF] Hiller | Design, analysis and realization of tendon-based parallel manipulators[END_REF]. As a consequence, if the maximal values of the actuator torques are not sufficiently high for a given payload, the reachable robot workspace will be reduced. This will restrict the zones where the robot can place the payload, and it may lead to disconnect the robot workspace into several aspects.

Most of the techniques developed for increasing static-wrench workspace feasibility have been applied for cable-driven parallel robots (CDPRs), in which the cables have been used as passive elements to find dynamic trajectories that seek to match the free response of the system with the desired payload motions [START_REF] Jiang | Periodic Trajectory Planning Beyond the Static Workspace for 6-DOF Cable-Suspended Parallel Robots[END_REF] [START_REF] Gosselin | Dynamic trajectory planning of a two-dof cable-suspended parallel robot[END_REF][6] [START_REF] Jiang | Dynamic point-to-point trajectory planning of a three-DOF cable-suspended parallel robot[END_REF]. The main idea is to replace the CDPR motor-cable arrangements with linear springs, from which an undamped mass-spring system is approximated. Then, the natural frequencies from such passive mechanical system are computed, and periodic trajectories are derived. Moreover, the natural frequencies of the virtual spring-based passive system are used for exploiting the dynamics of the system in order to find natural free-motion trajectories, thus increasing the feasible static-wrench workspace. Even if such techniques have shown their effectiveness when working with CDPRs, it is not the same case for rigid-link robots, in which no elastic element can be considered unless a spring is added in parallel to the active joints, as it is done in [START_REF] Barreto | The concept of natural motion for pick and place operations[END_REF][9] [START_REF] Hill | Minimizing Input Torques of a High-Speed Five-Bar Mechanism by using Variable Stiffness Springs[END_REF]. In those works, torsional springs are placed in parallel to the actuators, and they are used as energy storage in order to exploit the natural oscillations of pick-and-place robots in order to put the system near resonance, thus considerably decreasing the input torques.

Contrary to what it is done in [START_REF] Barreto | The concept of natural motion for pick and place operations[END_REF][9] [START_REF] Hill | Minimizing Input Torques of a High-Speed Five-Bar Mechanism by using Variable Stiffness Springs[END_REF], in which the spring potential energy is used to exploit the natural dynamics of the system, in this paper we propose to keep the robot unmodified, i.e. without any additional elastic element, and instead, to exploit the natural dynamics by using the payload potential energy. This will allow to increase the robot feasible static-wrench workspace thanks to a smart control of the natural oscillations in order to perform large amplitude motions between two disconnected aspects of the static-wrench workspace, similarly to the phenomenon occurring on an oscillating pendulum, in which all the energy is stored as potential energy when the payload is at the highest position of its swing, and it is then transferred continuously to kinetic energy. This paper is divided as follows: Section 2 presents a recall on dynamic modeling of robot manipulators. Section 3 shows the algorithm formulation for performing optimal trajectories that connect two aspects of the static-wrench workspace, divided due to the maximal actuator torques. Finally, in order to show the effectiveness of the proposed approach, section 4 presents the results on a 2-DOF serial robot with different payloads and different maximal actuator torque values.

Exploiting natural dynamics for enlarging feasible static-wrench workspace

It is well-known that near natural frequencies of oscillatory systems, the actuator torques can be drastically reduced [START_REF] Barreto | The concept of natural motion for pick and place operations[END_REF][9] [START_REF] Hill | Minimizing Input Torques of a High-Speed Five-Bar Mechanism by using Variable Stiffness Springs[END_REF]. For instance, we know that, near natural frequencies of a spring-mass system, large oscillatory motions can be obtained while conserving low torques in the actuator, which leads to lowering the energy expenditure. This is possible because of a smart exchange of energies between the kinetic energy and the spring potential energy of the system. Thus, in similar vein, in this section, we propose a strategy to smartly exploit the exchange between the gravity potential energy and the kinetic energy of a robot with payload, so that it can be put in near its free response, thus decreasing the input torques. In order to do so, we will firstly recall the dynamic modeling of robot manipulators. On a second stage, we will propose the formulation of a boundary value problem (BVP), which will allow us to compute the necessary optimal velocity, associated to the kinetic energy required to cross from one disconnected aspect to another with minimum effort. Then, similarly to an oscillating pendulum, thanks to the formulation of an optimisation problem, we will exploit the natural oscillations of the robot with payload to achieve the optimal velocity that match the system to its free response, obtained from the solution of the BVP.

Physical Background

Based on [START_REF] Dombre | Modeling, identification and control of robots[END_REF], let us briefly recall the dynamic modeling of robot manipulators. We will assume a kinematic architecture composed of n degrees of freedom and driven by n active joints. The position and velocity of the robot can be described by q and q representing the n-dimensional vectors of active joint variables and of active joint velocities, respectively.

Then, by using the Lagrange formalism, the dynamic model of a robot manipulator can be written as follows:

τ = M(q)q + h(q, q) + g r (q) + f s ( 1
)
where τ is an n-dimensional vector of the robot input efforts, M is an (n × n) definite positive matrix of inertia depending on the active joints coordinates q. h is an n-dimensional vector of Coriolis and centrifugal effects and its value depends on the active joint coordinates q, and their time derivatives q. g r is an n-dimensional vector grouping the gravitational effects. f s is an n-dimensional vector grouping the active joint friction terms.

Based on the aforementioned computations, it is thus necessary to find a strategy to exploit the natural dynamics of ( 1), so that we can find free-motion trajectories that constrain the input torques grouped in τ while performing large amplitude oscillations that connect aspects of the static-wrench workspace.

Shooting method for solving Boundary Value Problems

In the field of numerical analysis, a boundary value problem (BVP) is defined when, for solving a differential equation, two conditions are given at different values of the independent variable of the differential equation [START_REF] Roberts | Two-Point Boundary Value Problems: Shooting Methods[END_REF], i.e.:

h = f (t, h, h ), a < t < b (2)
with the boundary conditions defined by:

h(a) = A, h(b) = B (3)
where ( 2) is a second order differential equation, t is the independent variable, A and B are the desired boundary conditions defined at t = a and t = b of the BVP, respectively. The classical way to solve the BVP ( 2)-( 3), is by using the so-called shooting method [START_REF] Roberts | Two-Point Boundary Value Problems: Shooting Methods[END_REF], which treats the BVP as an initial value problem (IVP) [START_REF] Roberts | Two-Point Boundary Value Problems: Shooting Methods[END_REF], with the initial conditions given by:

h(a) = A, h (a) = α (4)
where α represents the first derivative of h(t) at t = a, or in other words, the slope of the solution, and the objective is to choose α, so that the solution of (2) satisfies the remaining boundary condition h(b) = B from (3). Practically speaking, the main idea of the shooting method is to numerically integrate (2), for different values of the slope α, and then to take the solution of the differential equation that satisfies the boundary condition h(b) = B. Thus, similarly, we propose to use a BVP formulation in order to find the optimal velocity qA , or slope, associated to the necessary kinetic energy required to pass between two different disconnected aspects of the static-wrench workspace (See Fig. 1a).

Shooting method applied on robot manipulators with heavy payloads

To develop the BVP for the robot, we will denote qA , as the optimal velocity at point A in one disconnected aspect of the static-wrench workspace, that will allow to connect the aspects in point B with null final velocity (See Fig. 1a).

In order to exploit the natural dynamics of the robot with payload to join two different aspects of the robot disconnected static wrench-workspace, in this section we propose to formulate a BVP, which can be solved with the shooting method [START_REF] Roberts | Two-Point Boundary Value Problems: Shooting Methods[END_REF]. The BVP formulation seeks to find the optimal velocity grouped in qA , or slope, and associated to the optimal kinetic energy required to connect two aspects, that solve (1) for τ = 0, while ensuring to meet the desired boundary positions defined in two disconnected aspects (A and B in Fig. 1a), i.e.: q = -M -1 (h(q, q) + g r (q) + f s ) [START_REF] Gosselin | Dynamic trajectory planning of a two-dof cable-suspended parallel robot[END_REF] with the boundary conditions defined as the desired positions in two disconnected aspects of the feasible static-wrench workspace by q(t 0 ) = q * A , q(t f ) = q * B with the formulation of the following position error vector:

E q ( qA ) = q(t f , qA ) -q * B ( 6
)
where t f is the final time, and q and q are obtained from numerical single and double integration of q, respectively. The shooting method is applied to iteratively find the proper joint velocities, grouped in qA , such that expression (6) converge to zero as follows:

while

E q > n , k ≤ max k do E q ( qA ) = q(t f , qA ) -q * B qA k : E q ( qA k ) = 0 qA k = qA k + ∆ qA k end
where qA k = qA k+1 representing the optimal set of joint velocities, max k is the maximum number of iterations, and n is the error threshold.

Since we will consider that all motions of the robot will start from the equilibrium position of the system, i.e. lowest end-effector location, at q 0 = [-π, 0] T , in the following section, we develop a strategy to optimally go from the equilibrium configuration velocity of the robot, to the optimal aspect-connecting velocity qA .

Initial oscillations for achieving optimal aspect-connecting velocity

In the aforementioned formulation, we have presented a strategy based on a BVP to compute the optimal velocity qA , associated to the kinetic energy required to join two disconnected aspects. Nevertheless, since we have chosen to start the motion from the equilibrium configuration of the robot ( q0 = 0), i.e. null initial stored potential energy U 0 to attain qA , it is thus necessary to find a strategy to reach qA from rest state. One could think of planning a classical motion based on a fifth-degree polynomial which can fix the position, velocity and acceleration at the trajectory extremities, thus allowing to fix qA . Nevertheless, since we seek to constrain the torques to operate within a range of motor limits, it is thus necessary to develop an optimal motion planner that integrates the effects of the dynamics to reach qA . That is why, contrary to the aformentioned BVP formulation, which allows us to compute the kinetic energy required to connect boundary positions in different aspects, in this part we seek to propose an optimal strategy to attain that aspect-connecting velocity qA . Similar to an oscillating pendulum, this can be done by exploiting the oscillations of the robot with payload in order to store the optimal payload potential energy to reach qA . Thus, we propose to find such exciting initial oscillatory motions that accelerates the robot to reach qA , while minimizing the input torques, by defining the following optimization formulation:

minimize J = t f 0 n i=1 τ 2 i dt over x DecVar subject to |τ i | ≤ τ max , t min ≤ t ≤ t max , |q i | ≤ q max ( 7
)
where the decision variable vector x DecVar contains the following parameters: q i k , v i k , a i k and t i k . Similar to the strategy based on the via-point motion planner of [START_REF] Hill | Increasing energy efficiency of high-speed parallel robots by using variable stiffness springs and optimal motion generation[END_REF], each of these parameters represents the amplitude, velocity, acceleration Fig. 1: On the left, the graphical interpretation of the disconnection between aspects due to actuation limits on a 2R serial robot workspace. A and B represent the boundary conditions which seek to be joined thanks to the BVP. On the right, the kinematic structure of the 2-DOF robot with payload is shown and time of oscillation for the exciting trajectories parameterized by using fifthdegree order polynomials, respectively (See Fig. 4, exciting oscillations). k is an integer representing the number of required oscillations (via points) for attaining qA . Finally, t f represents the final time of the optimization, and τ i represents the torque of the actuated joints grouped in q.

3 Case Study

2-DOF serial robot

In this Section, numerical validations are made thanks to a co-simulation between Matlab and ADAMS in order to show the effectiveness of our approach. The proposed strategy for enlarging the feasible static-wrench workspace was validated on a 2-DOF serial robot with payload, as shown in Fig. 1b. The active joint coordinates of the robot are parameterized by q = [q 1 , q 2 ] T . The vector of payload pose is denoted as x = [x, y] T . The link lengths l AB and l BC are identical, and the center of mass of each link is located at a distance l AS1 and l BS2 from the rotation center of the motors q 1 and q 2 , respectively. Moreover, the gravity field is directed along y 0 and is equal to g = [0, -g, 0] T , g > 0, being equal to 9.81 m/s 2 . Each link mass is parameterized by m 1 and m 2 , and the payload mass is denoted by m. Finally, the inertia of each link is given by zz 1 and zz 2 , and it is measured about its center of mass S 1 and S 2 , respectively. The dynamic and geometric parameters for the simulationes, were defined by the following numerical values: inertia of the two links zz 1 = zz 2 = 0.133 kg/m 2 , mass of each link m 1 = m 2 = 0.1 kg, Coulomb friction terms f s1 = f s2 = 0.3 Nm, link lengths l AB = l BC = 0.25 m and the distance to center of mass for each link l AS1 = l BS2 = 0.125 m. The dynamic equations for the 2-DOF serial robot with payload are not given here due to page limitations. However, by using the Lagrange formalism shown in section 2.1, the expression for the robot can take the form (1).

Results

The algorithm proposed in section 2 was tested by defining two desired points in disconnected aspects of the static-wrench workspaces, as it is shown in the interpretation from Fig. 1a. The aim is to connect both separated aspects of the static-wrench workspace by solving the corresponding BVP. Moreover, in order to show the effectiveness of the proposed approach, four cases corresponding to different payloads, maximal actuator torques and desired boundary positions (m, τ max , A, B, respectively) have been studied: Case 1: m = 0.5 kg,

τ max = 1 Nm, A = [0, -0.43], B = [0.05, 0.4]; Case 2: m = 1 kg, τ max = 2 Nm, A = [0.15, -0.4], B = [0, 0.4]; Case 3: m = 1.25 kg, τ max = 2.25 Nm, A = [0, -0.45], B = [0, 0.45]; Case 4: m = 1.5 kg, τ max = 3 Nm, A = [0, -0.45], B = [0.1, 0.4].
Firstly, in Fig 2, the total static-wrench feasible workspace is computed for the case 2, i.e. all the configurations of the robot for which statically with the given payload, respects the actuator limits (only one case is shown due to page limitations). Moreover, a static torque map based on contours is shown in order to have an idea of the value, that in statics, the actuator torque must be able to deliver for each robot configuration. It is worth mentioning that a more accurate contour map could be computed based on [START_REF] Ruiz | ARACHNIS: Analysis of robots actuated by cables with handy and neat interface software[END_REF], however, since for our case we have a simple system, we can easily compute the static torque at each position of the workspace. Once that the total static-wrench feasible workspace is computed, the algorithm proposed in section 2 is tested. The results in input torques for the four cases, for both BVP and exciting trajectories, are shown Table 1. In addition to that, the optimal values qA from the BVP, for the different boundary conditions from the four cases, are shown as well in Table 1. It can be seen that, based on the maximum value of the input torques, the efforts from the trajectory based on the shooting method are minimum, while the input efforts from the exciting trajectories are within the maximal motor torque values, thus constraining the actuator torques. In addition to that, Fig. 3 presents the input torques computed along the motion from case 2 in order to show the effectiveness of the approach to constrain the input efforts within a range of maximal actuator values that allows to join two disconnected aspects from the static-wrench workspace (only one case is shown due to page limitations). In Fig. 4 the optimal motions from the shooting method and from the exciting trajectories are shown in the workspace of the 2-DOF serial robot for case 2. It can be seen that thanks to the solution of the BVP, two points in different aspects of the divided staticwrench workspace can be joined while constraining the torques of the actuators. Furthermore, as it is shown in Table 1 and in Fig. 3, the input torques computed all configurations for which, the static torques are within the range of maximum torques |τ i | ≤ τ max . The center of both figures is empty due to the fact that a joint limit has been set for the second joint of q 2 = 2.2 rad. Each contour in these figures represent a value for the static torque of the first and second joint, respectively, for the case 2, i.e. a payload of 1kg and a maximal actuator torque τ 1max = τ 2max = 2 Nm. Furthermore, the static torques contours are computed for all robot configurations within the workspace boundaries.

from the BVP are closed to zero, which means a correct optimal dynamics exploitation since we were seeking to solve (1) for τ = 0. For a visualization of the results, a video of the co-simulations of the 2-DOF serial robot connecting two different aspects, can be found in https://bit.ly/2QuLQdb

Conclusions and future works

This paper proposes a strategy for increasing the feasible static-wrench workspace for robots by exploiting the robot dynamics through the solution of a BVP which seeks to smartly control the free response of the system to perform motions that connect two different aspects of the static-wrench workspace. This has allowed to considerably increase the reachability to disconnected areas of the workspace of a 2-DOF serial robot with different payloads and for multiple desired boundary positions. Simulations led to a successful connection of two different aspects of the feasible static-wrench workspace for different actuator torque limits and different payload weights. Future work on this subject includes the design of a prototype to experimentally validate the approach.

  (a) Division of static-wrench workspace in different aspects (b) Parameterization of 2R serial robot

  Slice XY in which each contour represents a static torque for the first joint. (b) Slice XY in which each contour represents a static torque for the second joint.

Fig. 2 :

 2 Fig.2: In red, the total static-wrench feasible workspace for case 2 is shown, i.e. all configurations for which, the static torques are within the range of maximum torques |τ i | ≤ τ max . The center of both figures is empty due to the fact that a joint limit has been set for the second joint of q 2 = 2.2 rad. Each contour in these figures represent a value for the static torque of the first and second joint, respectively, for the case 2, i.e. a payload of 1kg and a maximal actuator torque τ 1max = τ 2max = 2 Nm. Furthermore, the static torques contours are computed for all robot configurations within the workspace boundaries.

Fig. 3 :

 3 Fig. 3: Torques computed from the exciting-plus-BVP trajectory for the case 2.

Fig. 4 :

 4 Fig. 4: Optimal trajectory for connecting two points by solving the BVP from boundaries of case 2. A and B represent the boundary conditions of case 2 and they are joined thanks to the solution of the BVP.

Table 1 :

 1 Maximum values of input torques and optimal values for qA from the solution of the BVP.

		Max. Torques	Max. Torques	Optimal Velocity
		Shooting method for BVP Exciting trajectories	BVP
		|τ |max (Nm)	|τ |max (Nm)	qA (rad/s)
	Case 1	[0.007, 0.001]	[0.597, 0.604]	[-7.257, 4.494]
	Case 2	[0.007, 0.007]	[1.153, 2]	[-10.052, 10.179]
	Case 3	[0.011, 0.013]	[2.062, 2.105]	[-10.916, 10.625]
	Case 4	[0.032, 0.022]	[2.993, 2.999]	[-11.802, 12.649]

This work was conducted with the support of the Mexican Council for Science and Technology (CONACYT).


Rafael Balderas Hill

Damien Six

Abdelhamid Chriette

Sébastien Briot

Philippe Martinet

Crossing Type 2 Singularities of Parallel Robots without Pre-planned Trajectory with a Virtual-constraint-based Controller

The presence of Type 2 singularities in parallel robots severely affects their performances, mainly because the platform motion control is partially lost. It also leads to a size reduction of the operational workspace. Moreover, the dynamic model of the parallel mechanism degenerates and locally, the robot becomes underactuated in the singularity. It has been proven that it is possible to cross Type 2 singularities by respecting a dynamic criterion. Nevertheless, the controllers designed up to now require a pre-planned optimized trajectory including this criterion, and as a result, this strategy can only be used by qualified users. In order to avoid this drawback and to cross these types of singularities even if the trajectory is not pre-planned, this paper proposes a controller based on virtual constraints. Furthermore, the controller is integrated in a multi-control architecture in order to switch between a classical computed torque control far from the singularity and the virtual-constraint-based control law near to the singularity locus. Experimental results on a five-bar mechanism validated the automatic Type 2 singularity crossing.

I. INTRODUCTION

Parallel robots have several advantages compared to serial robots, such as higher stiffness, higher payload, lowest inertia and higher speed and acceleration. However, due to their geometric and kinematic complexity, singularity issues appear, which lead to the division of their workspace into several aspects corresponding to different assembly modes [START_REF] Merlet | Parallel Robots[END_REF] [START_REF] Chablat | Working modes and aspects in fully-parallel manipulator[END_REF]. For a general overview on the classification of the singularities in parallel robots the reader can refer to [START_REF] Gosselin | Singularity analysis of closed-loop kinematic chains[END_REF] [START_REF] Conconi | A new assessment of singularities of parallel kinematic chains[END_REF]. The literature review shows that some designers have overcome the problems of workspace reachability by designing mechanisms without singularities [START_REF] Gogu | Structural synthesis of fully-isotropic parallel robots via theory of linear transformations and evolutionary morphology[END_REF], but this usually leads to robots with small workspace or parallel robots with low practicability. One can also use redundancy [START_REF] Kurtz | Multiple-goal kinematic optimization of a parallel spherical mechanism with actuator redundancy[END_REF] but the use of additional actuators is costly. Another approach was presented in [START_REF] Briot | Optimal Force Generation of Parallel Manipulators for Passing through the Singular Positions[END_REF] in which it has been shown that it is possible to pass through Type 2 singularities by planning a trajectory respecting a physical criterion defined from the analysis of the degeneracy conditions of the dynamic model. This last approach is promising since it has shown its efficiency to cross the Type 2 singularities and to increase the reachable workspace thanks to a modification of the robot trajectory. However, the trajectory planning approach should be set more robust by combining a dedicated controller with it. Thus, a recent solution [START_REF] Pagis | Enlarging parallel robot workspace through Type-2 singularity crossing[END_REF] proposed the implementation of a Multi-Model Controller in order to track the trajectory that respects the physical criterion, and by switching between two different dynamic models. Nevertheless, the main issue with the current approach is that the definition of the optimal trajectory for crossing is required to be done by a wellskilled user, having some advanced knowledge in parallel robot dynamics. This is a brake towards automatizing the singularity crossing. The aim of this paper is to propose an advanced feedback control law dedicated to ensure unplanned singularity crossing allowing for an unskilled user to design a trajectory without specific knowledge in parallel robot dynamics. In order to solve this problem, we need to find a strategy which is able to control the underactuation of the end-effector of the parallel robot which occurs locally in Type 2 singularity [START_REF] Gosselin | Singularity analysis of closed-loop kinematic chains[END_REF]. In control of underactuated systems, several approaches exist [START_REF] Shiriaev | Constructive tool for orbital stabilization of underactuated nonlinear systems: Virtual constraints approach[END_REF], however, one of the most powerful tools is the enforcement of virtual constraints through feedback control. Virtual constraints are geometrical relations defined between controlled and uncontrolled variables that can be enforced through control action. Virtual constraints strategy has been applied to several applications, e.g. humanoid walking in [START_REF] Westervelt | Feedback Control of Dynamic Bipedal Robot Locomotion[END_REF], inverted pendulums in [START_REF] Canudas-De-Wit | Orbital stabilization of underactuated mechanical systems[END_REF], aerial systems in [START_REF] Ahmed | Kite generator system periodic motion planning via virtual constraints[END_REF]. Moreover, rigorous mathematical formulations have been presented in [START_REF] Canudas-De-Wit | Orbital stabilization of underactuated mechanical systems[END_REF] [START_REF] Shiriaev | Transverse linearization for controlled mechanical systems with several passive degrees of freedom[END_REF]. Thus, the use of virtual constraints for solving the problem of Type 2 singularity crossing seems to be a promising solution for dealing with the underactuated behavior of the robot at singularity given that at this specific singular posture, controlled and uncontrolled variables can be identified. Thus, in this paper, we develop a controller based on virtual constraint approach in order to cross the singularity without pre-planned trajectory. Our main contributions are listed below:

• Separation of the controlled and free dynamics, locally, in the singularity locus; • Feedback control law based on virtual constraints that respects the dynamic criterion/free dynamics, even if the trajectory is not pre-planned, avoiding the degeneracy of the dynamic model; • Multi-control architecture in order to continuously perform a trajectory tracking far from the singularity and in the neighborhood of the singularity.

II. PHYSICAL BACKGROUND

A. Dynamic modeling of parallel robots

Let us briefly recall the dynamic modeling of parallel robots presented in [START_REF] Briot | Dynamics of Parallel Robots: From Rigid Bodies to Flexible Elements[END_REF]. We will assume a parallel robot composed of n degrees of freedom. The actuation is given by n active joints. The position and velocity of the parallel robot can be described by using:

• q and q representing the n-dimensional vectors of active joint variables and of active joint velocities, respectively. • x and ẋ representing the n-dimensional vectors of platform pose and of its time derivatives, respectively. The input-output kinematic constraint that relates ẋ with the active joint velocities q is:

A ẋ + B q = 0 (1)
where A and B are the (n × n) parallel and serial kinematic Jacobian matrices, respectively [START_REF] Merlet | Parallel Robots[END_REF]. Then, by using the Lagrange formalism, the dynamic model of the robot can be written as follows:

τ = τ ta -B T λ ( 2
)
w p = A T λ (3)
where:

• τ is the n-dimensional vector of the robot input efforts;

• λ is the n-dimensional vector of Lagrange multipliers;

• τ ta and w p are the n-dimensional vectors related to the Lagrangian L of the system:

τ ta = d dt ∂L ∂ q T - ∂L ∂q T (4)
w p = d dt ∂L ∂ ẋ T - ∂L ∂x T (5)

B. General solution for the inverse dynamic model and degeneracy conditions of the dynamic model

This section seeks to give the general solution for the inverse dynamic model and to recall its degeneracy conditions due to the loss of rank of the kinematic Jacobian matrix A. Based on the classification of the singularities in [START_REF] Gosselin | Singularity analysis of closed-loop kinematic chains[END_REF], the parallel robots can encounter three singular configurations: Type 1 singularities occur when the mechanism loses its ability to move along one given direction. In such singular configurations the kinematic matrix B becomes rank deficient. Type 2 singularities occur when the mechanism gains some uncontrollable motions. In such singularities the kinematic matrix A becomes rank deficient. Type 3 singularities are configurations when matrices A and B are simultaneously rank deficient. The general solution for the inverse dynamic model of a parallel robot assuming the matrix A can be inverted, meaning that the lagrange multipliers λ from (3) can be substituted into (2), is given by:

τ = τ ta -B T A -T w p (6)
The solution of the dynamic model under the form [START_REF] Kurtz | Multiple-goal kinematic optimization of a parallel spherical mechanism with actuator redundancy[END_REF] assumes that the kinematic matrix A is a full rank matrix. Furthermore, for a prescribed trajectory in the Cartesian space, the values of the vectors q, q and q can be calculated from the inverse kinematics. Thus, considering that the parallel robot is not in a Type 1 singularity, the terms τ ta and w p can be computed. Nevertheless, for a trajectory passing through a Type 2 singularity, the matrix A is rank deficient [START_REF] Gosselin | Singularity analysis of closed-loop kinematic chains[END_REF]. In such singular configurations, the parallel mechanism gains one or more uncontrollable degrees of freedom making the system locally underactuated. An example of a Type 2 singularity is shown in Fig. 1 on a five-bar mechanism.

The degeneracy studied in the present research work will be limited to one uncontrollable degree of freedom at the singularity locus. The case of parallel robots with more degrees of degeneracy rarely appears.

III. MODELING OF THE FREE AND CONTROLLED DYNAMICS A. Separation of the free and controlled dynamics in the singularity locus

Locally in the singularity locus it is possible to separate the controlled and free dynamics. In order to do that let us express the term w p in (3), (5) under the following form given in [START_REF] Briot | Dynamics of Parallel Robots: From Rigid Bodies to Flexible Elements[END_REF]:

w p = M r ẍ + c r = A T λ (7)
where M r is an (n × n) definite positive matrix of inertia depending on the platform coordinates [START_REF] Briot | Dynamics of Parallel Robots: From Rigid Bodies to Flexible Elements[END_REF], and c r is an ndimensional vector of Coriolis, centrifugal and gravitational effects and its value depends on the platform coordinates x and their time derivatives ẋ. Knowing that the matrix A is rank deficient in Type 2 singularity, there exists a unit vector t s in its kernel, such that:

At s = 0 (8)
We take advantage of this property in order to build an orthonormal matrix V described as follows:

V =    t T s ----- V c    (9)
whose first row is the transpose of t s and V c is an ((n -1) × n) matrix grouping the other rows of V in order to fill the orthonormal basis. Then, multiplying [START_REF] Briot | Optimal Force Generation of Parallel Manipulators for Passing through the Singular Positions[END_REF] by the orthonormal matrix V, and decomposing the resulting equations, we obtain:

t T s w p = t T s (M r ẍ + c r ) = 0 (10)
V c w p = A T s λ (11)
where A s is an ((n -1) × n) reduced kinematic Jacobian matrix, as follows:

A T s = V c A T (12)
Expression [START_REF] Westervelt | Feedback Control of Dynamic Bipedal Robot Locomotion[END_REF] represents the local free dynamics of the robot in singularity, and equation ( 11) represents the controlled dynamics, given that the control inputs appear by substituting λ, computed from expression (2), into [START_REF] Canudas-De-Wit | Orbital stabilization of underactuated mechanical systems[END_REF]:

V c w p = A T s (B -T (τ -τ ta )) (13)
Expression ( 13) is valid since matrix B is not singular in a Type 2 singularity as disclosed in the singularity conditions of Section II-B.

B. Controlled and uncontrolled variables in the singularity locus

Based on the separation of the free and controlled dynamics, let us identify the controlled and uncontrolled variables in the singularity locus by performing a change of variable as follows:

ξ = VM r x (14)
where VM r is an invertible matrix considered locally constant since the underactuated behavior occurs restrictively only in the singularity locus. Then, the platform dynamics can be expressed with the new variable ξ by computing the following expression:

ξ = VM r ẍ ( 15
)
where ξ is the acceleration of the variable change. By introducing the change of variable [START_REF] Briot | Dynamics of Parallel Robots: From Rigid Bodies to Flexible Elements[END_REF] into equations [START_REF] Westervelt | Feedback Control of Dynamic Bipedal Robot Locomotion[END_REF] and [START_REF] Canudas-De-Wit | Orbital stabilization of underactuated mechanical systems[END_REF], the free and controlled dynamics in the coordinate ξ become:

ξ1 + t T s c r = 0 (16) ξc + V c c r = A T s λ (17)
where ξ 1 is the first term of ξ and is the uncontrolled variable and ξ c = [ξ 2 , ξ 3 , ..., ξ n ] T is the (n-1)-dimensional vector of controlled variables, grouping all the remaining components of ξ.

IV. DYNAMIC MODEL IN THE SINGULARITY LOCUS

The inverse dynamic model can be expressed as a function of the acceleration of the controlled variables only. According to [START_REF] Briot | Dynamics of Parallel Robots: From Rigid Bodies to Flexible Elements[END_REF], equation ( 4) can be obtained under the form:

τ ta = M a q + c a (18)
with M a as an (n × n) matrix and c a as an n-dimensional vector. Let us now substitute (18) into (2) and compute w p in the coordinate ξ, then the dynamic model can be expressed:

τ = M a q + c a -B T λ (19)
w p = ξ + Vc r (20)
The expression (20) can be split into the free and controlled dynamics by the equations ( 16) and (17). In order to find an exact solution for the inverse dynamic model under the form (19) and (20), it is possible to multiply (19) by A T s B -T :

A T s B -T τ = A T s B -T (M a q + c a ) -A T s λ (21)
Then, A T s λ from (21) can be substituted by using the expression of the controlled dynamics (17), and equation (21) becomes:

A T s B -T τ = A T s B -T (M a q + c a ) -ξc -V c c r (22)
In order to express the dynamic model in the Cartesian space, q is mapped into the Cartesian space by using the second order inverse kinematic model, obtained from the differentiation with respect to time of equation ( 1):

q = J inv ẍ + b (23)
where

J inv = -B -1 A and b = -B -1 ( Ȧ ẋ + Ḃ q).
We can perform the change of variable from ẍ into ξ by using expression [START_REF] Liberzon | Switching in Systems and Control[END_REF], as follows:

ẍ = M -1 r V T ξ (24)
Then, to separate the controlled and uncontrolled variables from (24), the following computations are performed:

V T ξ = t s ξ1 + V T c ξc (25)
By introducing expression (25) into equation ( 24), and substituting into the second order inverse kinematic model ( 23), the model (22) becomes:

A T s B -T τ = A T s B -T (M a J inv M -1 r (t s ξ1 + V T c ξc ) +M a b + c a ) -ξc -V c c r ( 26
)
The expression of the free dynamics ( 16) can be rewritten as follows:

ξ1 = -t T s c r (27)
Therefore, by substituting ( 27) into (26) and grouping the controlled dynamics, the following expression is obtained:

J T s τ = M s ξc + c s ( 28
)
with the matrix J s = -B -1 A s , a generalized inertia matrix reduced to the controllable degrees of freedom

M s = J T s M a J inv M -1 r V T c + I (n-1)
and the term that contains the Coriolis, centrifugal, gravity and friction effects

c s = J T s (M a b + c a -M a J inv M -1 r c r ) + V c c r .
From the properties of the Moore-Penrose pseudo-inverse, the inverse dynamic model that minimizes the input torques at the singularity locus is expressed by:

τ = J T + s (M s ξc + c s ) (29)
Since the uncontrolled variable has been removed from expression (29), the matrix J T s is full row rank and the obtained inverse dynamic model does not degenerate in the singularity locus.

V. MODELING OF THE VIRTUAL CONSTRAINTS AND CONTROL LAW IN THE SINGULARITY

In sections III-B and IV, we showed that locally in a Type 2 singularity, the dynamic model of a parallel robot can be separated in a free dynamics, associated to the local uncontrolled direction, and a controlled dynamics in the other directions. Virtual-constraint-based controllers have been shown to be efficient in the control of underactuated systems. In this section, we will apply this tool to maintain the controllable variables on a trajectory defined by the virtual constraints while the uncontrolled dynamics of the system remains free. The controller developed naturally respects the controlled and uncontrolled dynamics of the robot, avoiding model degeneracy issues obtained with a classic dynamic controller as Computed Torque Control. In [START_REF] Shiriaev | Transverse linearization for controlled mechanical systems with several passive degrees of freedom[END_REF], the authors have defined the theoretical framework in the modeling of virtual constraints. They defined virtual constraints as geometrical relations among the generalized coordinates that virtually constrain the system through the application of control inputs for keeping those constraints invariant. Recalling that the controlled and uncontrolled variables at the singularity locus are

• ξ 1 : the uncontrolled variable; • ξ c = [ξ 2 , ξ 3 , ..., ξ n ] T : the controlled variables.
It is thus possible to define a relation between the controlled and uncontrolled variables in order to define the virtual constraint to be enforced around the Type 2 singularity as follows:

ξ c = h(ξ 1 ) (30)
In order to enforce the constraints, outputs are defined as

y = ξ c -h(ξ 1 ) (31)
The controller will drive the outputs y towards zero through feedback linearization. To perform this linearization, the outputs are differentiated with respect to time twice

ÿ = ξc -( ξ1 h (ξ 1 ) + ξ2 1 h (ξ 1 )) (32)
where h and h represent the first and second derivatives with respect to the passive variable ξ 1 of the function h(ξ 1 ), respectively. The uncontrolled dynamics, known by equation ( 27) are introduced in equation (32).

ÿ = ξc -(-t T s c r h (ξ 1 ) + ξ2 1 h (ξ 1 )) (33)
The objective of the controller is to obtain a convergence of the output y towards 0 as a second order ordinary differential equation. To obtain it, an auxiliary input υ corresponding to the acceleration of the controlled variables ξ c is introduced, and defined to ensure the convergence of the virtual constraints

υ = (-t T s c r h (ξ 1 ) + ξ2 1 h (ξ 1 )) -K p y -K d ẏ (34)
where K p and K d are positive definite matrices of proportional and derivative gains to be tuned. Imposing torques such that ξc = υ will lead to the closed-loop equation

ÿ + K d ẏ + K p y = 0 (35)
ensuring the convergence of y towards zero. The auxiliary input is obtained by an input-output linearization through the dynamic model of the controlled variables (29).

τ = J T + s (M s υ + c s ) (36)
The designed virtual-constraint-based control law (Fig. 2) only applies for the situation when the robot becomes underactuated in the singularity locus. Nevertheless, the goal is to continuously and automatically track a trajectory in the workspace far from the singularity and in the singularity locus. Thus, the present paper proposes the use of a multicontrol architecture that switches between two control laws.

In [START_REF] Pagis | Enlarging parallel robot workspace through Type-2 singularity crossing[END_REF] it has been shown that in order to track a given trajectory out of the singularity locus, a classical computed torque control can be implemented. The second control law is the virtual-constraint-based controller for tracking in the singularity locus.

In order to design a Multi-Controller scheme, a supervisory control architecture is implemented (for a global overview on Supervisory Control the reader can refer to [START_REF] Liberzon | Switching in Systems and Control[END_REF]). Fig. 3 shows the Multi-Controller scheme integrating the control law based on virtual constraints u 1 and the classical computed torque control u 2 .

-0.

→ B → A → B → A → C → D → A → F → E → A.
to the other. The five-bar mechanism crosses the singularity without torque discontinuities. In order to validate the virtualconstraint-based control law in the multi-control architecture, an IMU was attached in the platform of the robot. Fig. 7 shows the tracking error of the magnitude of the acceleration of the end-effector. The error is expressed as a function of the desired trajectory. Since the controller automatically corrects the trajectory when crossing the singularity with respect to the controlled and uncontrolled dynamics, the perturbation of the trajectory tracking when crossing is expected. Furthermore, model errors may have an influence over the perturbations observed. Nevertheless, the overall error in acceleration remains converging towards zero, and varying in a range of 10 -1 m/s 2 . The robustness in experimentation of the proposed control architecture was observed by testing the trajectory at least 10 times. In all the tests the robot successfully crosses the singularities automatically.

VII. CONCLUSIONS

This paper proposes a new solution to cross Type 2 singularities on a parallel robot. An analysis of the dynamic equations in the singularity locus provided the separation of the uncontrolled and controlled dynamics. Then, a specific controller using an approach based on virtual constraints was developed to deal with the local underactuation in the singularity. A multi-control law was then implemented to toggle between the local controller near the singularity locus and a classic computed torque control in the rest of the workspace. This controller allowed the crossing of singularities without the pre-planning of an optimal trajectory based on the dynamic criterion identified in previous studies. Experimental implementations led to the successful crossing of singularities in several configurations on a five-bar mechanism. Future work on this subject includes the theoretical proof of robustness of the multi-control law implemented and the design of alternative constraints to minimize the divergence from the initial trajectory when performing the singularity crossing.

Fig. 1 .

 1 Fig. 1. A five-bar mechanism gains an uncontrollable motion perpendicular to -----→ A 12 A 13 and-----→A 22 A 13 . The vector ts represents the direction of the uncontrollable motion, which makes the robot underactuated. f 1 and f 2 are two forces applied by the actuators on the platform through the passive legs.

Fig. 5 .

 5 Fig. 5. Multiple crossing trajectory with the sequence: A→ B → A → B → A → C → D → A → F → E → A.

2 Fig. 6 .

 26 Fig. 6. Input torques for the first and the second actuator.

Fig. 7 .

 7 Fig. 7. Tracking error of the magnitude of the acceleration for the multicontrol law.

Ecole Centrale de Nantes, 44321 Nantes, France

Laboratoire des Sciences du Numérique de Nantes (LS2N), UMR 6004, 44321 Nantes, France {Rafael.Balderas-Hill, Damien.Six, Abdelhamid.Chriette, Sebastien.Briot, Philippe.Martinet}@ls2n.fr

Parallel Robot

Linearized System PD Controller Controlled Outputs Virtual Constraints Fig. 2. Virtual-constraint-based Control scheme is depicted with the vector of control inputs τ representing the input torques in the singularity locus. The box in blue highlights the linearized system, and the PD controller ensures the convergence of the virtual constraints.

Virtual-constraint-based Controller

VI. CASE STUDY A. Five-bar mechanism: Benchmark in Type 2 singularity crossing

The proposed control approach for crossing Type 2 singularities was validated in a five-bar mechanism from Meca-demic, as shown in Fig. 4. The actuation is provided by q = [q 11 , q 21 ] T . The vector of passive joints is given by q d = [q 12 , q 22 , q 13 ] T . The vector of moving platform pose is given by x = [x, y] T . a is the distance between the actuated joints. All link lengths l are identical. The dynamic model of the five-bar mechanism from Mecademic with identified parameters was given in [START_REF] Briot | Controlling Parallel Robots during Singular Assembly Mode Changing[END_REF]. The dynamic model has the following form:

where A and B are computed from the input-output kinematic constraint relation [START_REF] Merlet | Parallel Robots[END_REF]. The matrix ZZ is positive diagonal resulted from the identified dynamic parameters, and f groups the active joint friction terms.

B. Virtual-constraint-based control implementation

In the case of the five-bar mechanism, the general expression of the matrix A T in a Type 2 singular configuration as in Fig. 1

is given by

with β the orientation of the alignment of the passive legs β = q 11 + q 12 + q 13 = q 21 + q 22 (40)

In order to obtain the separation between the free and the controlled dynamics, the matrix V of equation ( 9) is defined by

and the equation ( 38) is pre-multiplied by the matrix

The variable ξ = ξ 1 ξ c T is defined by

Giving the uncontrolled dynamics

and the inverse dynamic model (29) associated to the controllable variables for the five-bar mechanism at singular configuration τ = J T + s (M s ξc + c s ) (45)

with

. This inverse dynamic model will be used to perform the inputoutput linearization in the control loop (see section V). The enforced virtual constraint is simply ξ c = 0. Then the auxiliary control law υ applied on the controlled acceleration ξc is given by υ

C. Multi-Control architecture for integrating the control laws

The performance index considered in the supervisor block from Fig. 3, for the control laws transition, is the proximity to the singularity locus. We use the approach of computing the condition number of the matrix A, which is the kinematic matrix that degenerates in the singularity locus. Any other criterion characterizing the proximity to singularity could be used. Thus, the multi-control law can be synthesized by:

where σ is equal to 0 far from the singularity and equal to 1 in the singularity locus:

is a threshold to be tuned experimentally.

D. Trajectory generation

In order to validate the theoretical formulations, the controller is tested by defining sections of trajectory crossing the singularity using fifth-degree order polynomials which can fix the position at the sections extremities with null initial and final velocities and accelerations. As the uncontrolled dynamic is taken into account in the controller, it is not necessary to design those sections such that the dynamic criterion established in previous studies [START_REF] Briot | Optimal Force Generation of Parallel Manipulators for Passing through the Singular Positions[END_REF] is respected. The trajectory will be automatically modified to respect the free dynamics when crossing the singularity and the classic CTC controller catches up with the initial trajectory after the crossing. The final trajectory consists of multiple sections crossing the singularity locus switching between the two different assembly modes of the five-bar mechanism. The sections of trajectories are defined between the points given by A = [0, 0.

E. Experimental results

This section presents the results of testing the controller for the trajectory defined in section VI-D. Fig. 6 shows the input torques computed from the multi-control law along the trajectory with multiple crossings from one assembly mode


Adrien Koessler
email: adrien.koessler@sigma-clermont.fr

Alexandre Goldsztejn
email: alexandre.goldsztejn@ls2n.fr

Sébastien Briot
email: se-bastien.briot@ls2n.fr

Nicolas Bouton
email: nicolas.bouton@sigma-clermont.fr

Certified Detection of Parallel Robot Assembly Mode under Type 2 Singularity Crossing Trajectories

Increasing the size of operationnal workspace is one of the main problems parallel robots are faced with. Among all the proposed solutions to that, crossing Type 2 singularities using dedicated trajectory generation and multimodel controller has a great potential. Yet, this approach is not sufficient for the robot to operate autonomously, as assembly mode detection during the motion currently requires additional redundant information.

To tackle this problem, we propose an algorithm based on Interval Analysis (IA) that is able to track the end-effector of the robot even under assembly mode change. IA-based solvers for the forward kinematic problem of parallel robots are well known, but they cannot be used under assembly mode change. Compared to those classical approaches, the major modification introduced is the tracking of end-effector velocity in addition to its pose. Using this new information of velocity, the algorithm is capable to monitor the assembly mode change of the robot happening when the singularities are crossed. The behavior and the reliability of this algorithm are analyzed experimentally on a five-bar planar parallel mechanism.

I. INTRODUCTION

It is acknowledged that parallel manipulators have advantages in terms of rigidity, cycle time and payload-to-weight ratio over their serial counterparts. Yet, their workspace is usually split into several assembly modes separated by Type 2 (or parallel) singularities [START_REF] Gosselin | Angeles 90 -Singularity Analysis of Closed-Loop kinematic Chains[END_REF] [START_REF] Conconi | A New Assessment of Singularities of Parallel Kinematic Chains[END_REF].

These are difficult to deal with as one (or more) manipulator's degree of freedom becomes uncontrollable. Assembly modes separated by Type 2 singularities cannot be interconnected to each other using classical control algorithms, resulting in a decrease in the size of the operationnal workspace.

Therefore, approaches leading to operationnal workspace enlargement have been developed, such as:

• Optimal design of parallel manipulators to maximize the size of the operationnal workspace [START_REF] Figielski | Towards development of a 2-DOF planar parallel robot with optimal workspace use[END_REF] or to eliminate Type 2 singularities [START_REF] Carricato | Singularity-free fully-isotropic translational parallel manipulators[END_REF][5] • Redundantly actuated robots [START_REF] Merlet | Redundant Parallel Manipulators[END_REF] [START_REF] Müller | Redundant Actuation of Parallel Manipulators[END_REF], which is a costly solution or robots with variable actuation modes [START_REF] Arakelian | Mechanism and Machine Theory Increase of singularity-free zones in the workspace of parallel manipulators using mechanisms of variable structure[END_REF][9] • Assembly mode changing trajectories based on cusp points detection [START_REF] Zein | Singular Curves in the Joint Space and Cusp Points of 3-RPR Parallel Manipulators[END_REF] • Dedicated trajectory generation and control to cross directly the singularity locus [START_REF] Briot | Optimal Force Generation in Parallel Manipulators for Passing through the Singular Positions[END_REF][12] The latter approach seems more versatile to us, as it can be generalized to any type of parallel robot. It emerged with the discovery of a physical criterion allowing the non-degeneracy of the manipulator's dynamic model in Type 2 singular configurations [START_REF] Briot | Optimal Force Generation in Parallel Manipulators for Passing through the Singular Positions[END_REF]. Using this criterion, trajectory planners and dedicated multi-model Computed Torque Control schemes [START_REF] Pagis | Design of a Controller for Enlarging Parallel Robots Workspace through Type 2 Singularity Crossing[END_REF] [START_REF] Ozdemir | A switching inverse dynamics controller for parallel manipulators around drive singular configurations[END_REF] have been developed and tested, proving their efficiency in most cases. Nevertheless, unavoidable tracking errors with respect to this desired trajectory can have a strong impact in the vicinity of the singularity, up to missing the assembly mode change. Hence monitoring the assembly mode during the robot motion is needed.

Methods to detect the assembly mode of a robot have already been investigated. The introduction of measurement redundancy (using passive instrumented legs or encoders in passive joints) allows to achieve it, but the placement of encoders needs to be chosen carefully [START_REF] Tancredi | De la simplification et la résolution du modèle géométrique direct des robots parallèles[END_REF]. External measurement with vision, either to monitor directly the pose of the platform [START_REF] Martinet | Vision Based Control Law using 3D Visual Features[END_REF] or the orientation of the legs of the mechanism [START_REF] Ozgür | Dynamic Control of the Quattro Robot by the Leg Edgels[END_REF], has also been investigated. Both of those approaches come at an extra cost and might be difficult to integrate in an industrial environment.

An interesting alternative is based on interval analysis. IA is a field of mathematics that treats variables as intervals, meaning that the variable's value must lie between the lower bound and the upper bound of the interval [START_REF] Moore | [END_REF]. Thus, it can be very useful for reliable computations in boundederror context. Numerous applications of IA can be found in robotics [START_REF] Merlet | Interval Analysis and Reliability in Robotics[END_REF] and control, such as mobile robot localization [START_REF] Bars | Interval SLAM for underwater robots; A new experiment[END_REF] [START_REF] Jaulin | Range-Only SLAM with Occupancy Maps: a Set-Membership Approach[END_REF], robot kinematic calibration [START_REF] Daney | Interval method for calibration of parallel robots: Vision-based experiments[END_REF], robust control [START_REF] Vehí | A Survey of Applications of Interval Analysis to Robust Control[END_REF], positioning precision analysis and synthesis [START_REF] Pac | Interval Analysis for Robot Precision Evaluation[END_REF] [START_REF] Goldsztejn | A three-step methodology for dimensional tolerance synthesis of parallel manipulators[END_REF]. In parallel robotics in particular, methods to compute the aspects of the workspace have been discovered [START_REF] Oetomo | An Interval-Based Method for Workspace Analysis of Planar Flexure-Jointed Mechanism[END_REF] [START_REF] Caro | A branch and prune algorithm for the computation of generalized aspects of parallel robots[END_REF].

Merlet has applied IA to solve the Forward Kinematic (FK) problem of a Gough-Stewart platform [START_REF] Merlet | Solving the Forward Kinematics of a Gough-Type Parallel Manipulator with Interval Analysis[END_REF], using the equations of the inverse kinematics and a start point. Thanks to constraint propagation techniques, his algorithm showed ability to distinguish the different solutions of FK problem even when the platform is close to Type 2 singularities. Such a technique, which is time-efficient and does not require additional sensors, is highly interesting in an industrial context. Yet it has to be adapted to assembly mode changing detection, as it does not handle it: if the robot reaches a Type 2 singularity and exits it, the algorithm will compute enclosures of the pose that diverge, and therefore will not be able to state in which assembly mode it exited [START_REF] Merlet | Solving the Forward Kinematics of a Gough-Type Parallel Manipulator with Interval Analysis[END_REF].

Hence, our goal is to provide an IA-based algorithm that is able to track the assembly mode of a robot attempting to cross Type 2 singularities.

Notations: In the following, bracketed symbols are used to denote interval enclosures of the corresponding quantities. For instance, a scalar a ∈ R car be enclosed by the interval [a] = [a, a] ∈ IR such as a ≤ a ≤ a, where IR is the set of scalar real intervals. An interval vector is the cartesian product of the scalar interval of each component. For instance, a vector p = (p

1 • • • p n ) T ∈ R n admits a bracketing [p] = [p 1 , p 1 ] × • • • × [p n , p n ] ∈ IR n .
This explains why a vector of intervals is generally called a box.

II. POSE AND VELOCITY TRACKING ALGORITHM

A. Kinematic modelling of parallel mechanisms

This section aims to briefly recall the kinematic equations of a parallel manipulator with n degrees of freedom and driven by n actuators. The position and the speed of the manipulator can be fully described using:

• q and q that represent respectively the vector of active joints variables and active joints velocities, • x and ẋ that are the end-effector pose parameters and their derivatives with respect to time, respectively. In the following, the sampling time is t s with time step k, while x k = x(kt s ), ẋk = ẋ(kt s ), q k = q(kt s ) and qk = q(kt s ).

The loop-closure equation of the mechanism can then be written as [START_REF] Merlet | Parallel Robots[END_REF]:

f(x, q, ξ ξ ξ ) = 0 (1)
where ξ ξ ξ is the vector of geometrical parameters. The firstorder kinematic equation is:

J x (x, q, ξ ξ ξ ) ẋ + J q (x, q, ξ ξ ξ ) q = 0 (2)
where J x and J q are the (n × n) Jacobian matrices of the robot.

B. Outline of the algorithm

Tracking the pose given joint values information qk is usually done by solving (1) using a local solver, usually a Newton method, and the previous pose as initial condition. This is however unreliable, especially near singularities (See [START_REF] Merlet | Parallel Robots[END_REF], chapter 4). IA gives rise to a reliable variant of this tracking process [START_REF] Merlet | Solving the Forward Kinematics of a Gough-Type Parallel Manipulator with Interval Analysis[END_REF]: Given an enclosure [q] k on the joint coordinates at step k, one can compute an enclosure [x] k of the pose solving (1) with an IA-based solver. Merlet [START_REF] Merlet | Solving the Forward Kinematics of a Gough-Type Parallel Manipulator with Interval Analysis[END_REF] also proposed to use the previous pose enclosure [x] k-1 and an upper bound on the speed of the end-effector, given here under the form of an enclosure of all possible speed vectors [v max ] ∈ IR n , to obtain an initial domain [x] k for the IA-based solver:

x k = x k-1 + t s 0 v(t) dt ∈ [x] k-1 + t s [v max ] =: [x] k . (3)
This initial domain drastically simplifies and speeds up the solving process, since it is usually quite small and contains only the solution corresponding to the tracked trajectory. As mentioned in the previous section, this IA-based tracking process fails and computes divergent enclosures in the vicinity of singularities [START_REF] Merlet | Solving the Forward Kinematics of a Gough-Type Parallel Manipulator with Interval Analysis[END_REF]. We propose here to extend this IA-based tracking process by tracking the end-effector velocity in addition to its pose. Using the end-effector velocity information, we might be able to tell in which assembly mode the robot has exited the singularity. For instance, imagine that the robot's endeffector follows a path intersecting a Type 2 singularity; it may cross the singularity or "bounce" on it. If we add the information that the speed along this path did not change its sign during the whole transfer, it is obvious that the robot did cross the singularity. This is the type of information that we need to create within our algorithm.

An example is shown on Fig. 1. With the knowledge that the end-effector only had negative velocity along the Y -axis, the only feasible final configuration is f 2 .

The end-effector speed is tracked similarly to its pose. An initial domain for the end-effector velocity at step k is computed using the enclosure at step k -1 and a given maximal acceleration [a max ] ∈ IR n :

ẋk = ẋk-1 + t s 0 a(t) dt ∈ [ẋ] k-1 + t s [a max ] =: [ẋ] k . (4)
This enclosure can be safely intersected with [v max ]. Then using an enclosure [ q] k of the joint velocity, this initial domain is sharped by solving (2) using an IA-based solver. In fact, this velocity enclosure can be used to improve the initial pose enclosure (3

): Since 0 ∈ [a max ] we have t s [a max ] = [0,t s ][a max ], and therefore for all t ∈ [(k -1)t s , kt s ] ẋ(t) = ẋk-1 + t 0 v(t) dt ∈ [ẋ] k-1 + [0,t s ][a max ] = [ẋ] k , (5)
i.e., the end-effector velocity enclosure (4) at time kt s is actually an enclosure of the velocity for all t ∈ [(k -1)t s , kt s ]. As a consequence, (3) can be improved as follows:

[x] k := [x] k-1 + t s [ẋ] k . (6)
As previously, the solver will fail sharpening the initial end-effector pose and velocity domains in the vicinity of singularities, since the matrix J x in (2) is ill-conditionned. Thus, in such configurations, the maximal acceleration is the only information that can be used to enclose the end-effector velocity. It follows from (4) that the width of its enclosure will grow approximately as k t s [a max ] . Despite its constant growth, it is expected that the velocity enclosure will give enough information to remove all the parasite solutions to the FK problem, thus allowing to detect the current assembly mode. For instance, in the situation of Figure 1, the endeffector velocity enclosure can remain negative in spite of its constant growth during the singularity crossing provided that the speed is high enough.

C. Description of the algorithm

Algorithm 1 provides the pseudo-code of the end-effector pose and velocity tracking algorithm. The procedure read reads measurement of the joint coordinate and velocity (the latter being usually approximated using possibly filtered finite differences). It returns true in case of successful acquisition, false otherwise when the motion is finished. The procedure contract uses standard IA-based algorithms (numerical constraint propagation, including the interval Newton and the interval Gauss-Seidel operators [START_REF] Moore | [END_REF], the forward-backward contractor [START_REF] Benhamou | Revising hull and box consistency[END_REF], etc.) in order to contract the domains [x] k and [ẋ] k subject to the given constraints. The domains [ξ ξ ξ ], [ε ε ε q ] and [ε ε ε q] are considered as parameters of the constraint system and are not contracted.

The inputs of the algorithm are:

• The kinematic model f, from which are computed the derivatives J x and J q . These expressions are suitable to perform interval evaluations and numerical constraint programming. • The joint and joint velocity measurement error bounds [ε ε ε q ] and [ε ε ε q], given under the form of intervals which contain all possible errors. • The sampling time t s . Provided that the initial enclosure and the bounds involved in the algorithm are valid, the IA-based constraint solving process ensures that the true end-effector pose and velocity are enclosed correctly for all times. The next section discuss in detail the DexTAR parallel robot, its model and the bounds that will be used in the experiments.

III. CASE STUDY: ASSEMBLY MODE CHANGE OF A FIVE-BAR ROBOT

A. Robot kinematic analysis

For the purpose of the experiments, a five-bar planar parallel robot is used: the DexTAR, manufactured by the company Mecademic. The design of the robot is described in [START_REF] Figielski | Towards development of a 2-DOF planar parallel robot with optimal workspace use[END_REF]. A picture of the robot and a scheme are presented in Fig. 2. It exhibits the two joint inputs q 1 , q 2 (corresponding to motors located at A 1 , A 2 ) and the seven geometrical Algorithm 1: Pose and velocity tracking algorithm

Input: f : R n × R n × R m → R n ; [x] 0 , [ẋ] 0 , [v max ], [a max ] ∈ IR n ; [ξ ξ ξ ] ∈ IR m ; [ε ε ε q ], [ε ε ε q] ∈ IR n ; t s > 0; 1 k ← 0; 2 while read( qk , qk ) do 3 [q] k ← qk + [ε ε ε q ]; 4 [ q] k ← qk + [ε ε ε q ]; 5 constract [x] k and [ẋ] k subject to (1) and (2); 6 [ẋ] k+1 ← ([ẋ] k + t s [a max ]) ∩ [v max ]; 7 [x] k+1 ← [x] k + t s [x] k ; 8 k ← k + 1; 9 end A 1 q 1 +δ 1 A 2 q 2 +δ 2 B 1 B 2 E(x, y)
f i = (x ± d 2 -i1 cos(q i1 )) 2 + (y -i1 sin(q i1 )) 2 -2 i2 = 0. (7)
Enclosures [ξ ξ ξ ] of the values of the geometrical parameters ξ ξ ξ are obtained through calibration procedure. It is performed using a Leica AT901 laser tracker. The adapted DETMAX algorithm described in [START_REF] Joubair | Comparison of the efficiency of five observability indices for robot calibration[END_REF] gave the 50 optimal poses to be used for calibration. Using a least-square routine [START_REF] Khalil | Modeling, Identification and Control of Robots[END_REF], estimated values ξi were obtained for parameters. Under the assumption that the measurement noise is white and that the identification matrix is deterministic, an estimate σi of the standard deviation of each paramenter was also computed. For every parameter ξ i , the interval [ξ i ] = [ ξi -5 σi ; ξi + 5 σi ] was consequently set. Nominal values and parameter intervals are presented in Table I.

Daney remarked that such a result can be rigorously verified by using a purely geometric IA-based solver [START_REF] Daney | Interval Methods for Certification of the Kinematic Calibration of Parallel Robots[END_REF]. If the interval [ξ i ] are valid, it is ensured that the enclosure [x] computed by (1) cannot be equal to the empty set at any point of the trajectory. Sweeping the whole workspace (singularity locus included), we verified it was never the case in our application, thus validating the presented results.

Determining the bounds [v max ] and [a max ] is by far the most challenging problem in relation with our algorithm. To the best of our knowledge, no method exist to determine the fullcycle maximum values of end-effector velocity and acceleration for parallel mechanisms, as the problem degenerates in Type 2 singular configurations. For our experimentations, the real value of end-effector acceleration norm ẍ max is monitored using an inertial unit, so that an upper bound for end-effector acceleration can be deduced. The maximal acceleration measured during the three experiments presented in the next section are respectively 65.5 m • s -2 , 79.9 m • s ). The sampling time is t s = 0.001s and [ε ε ε q ] = [-π 34000 , π 34000 ] × [-π 34000 , π 34000 ] (rad) is given by the coder precision. The joint velocity is classically approximated by a finite difference qk = qk -qk-1

t s with error bound [ε ε ε q] = 2[ε ε ε q ] t s .

B. Trajectory Generation and Control

As the goal is to cross Type 2 singularities to perform our tests, an adapted controller is needed. The choice we made is the trajectory generator and the multi-model Computed Torque Controller (CTC) developed by Pagis in [START_REF] Pagis | Design of a Controller for Enlarging Parallel Robots Workspace through Type 2 Singularity Crossing[END_REF].

Adapting CTC schemes to singularity crossing is tricky, since it relies on the inverse dynamic model (IDM) which degenerates in Type 2 singularities. Yet, a physical criterion can be used to keep the model consistent. It states that the acceleration undergone by the robot's end-effector must be normal to the direction of the uncontrolled motion of the end-effector. The reader is referred to [START_REF] Briot | Optimal Force Generation in Parallel Manipulators for Passing through the Singular Positions[END_REF] for comprehensive discussion about this criterion.

This result was used by Pagis [START_REF] Pagis | Design of a Controller for Enlarging Parallel Robots Workspace through Type 2 Singularity Crossing[END_REF] to create a trajectory generator dedicated to crossing Type 2 singularities. In addition, a multi-model computed torque control scheme was created by the same author, in which a nondegenerating expression of the IDM is used in the neighbourhood of the singularity locus. Complete information can be found in [START_REF] Pagis | Enlarging parallel robot workspace through Type-2 singularity crossing[END_REF].

This control algorithm, left unmodified, is implemented in Matlab/Simulink development environment. A Quanser Q-PID FPGA is used to run the control programs and to communicate with the robot.

IV. EXPERIMENTS

A. Tracking Algorithm Implementation

The algorithm is implemented in C++ language and compiled using MinGW. The program uses the IBEX library [START_REF] Ninin | Global Optimization based on Contractor Programming : an Overview of the IBEX library[END_REF] to perform interval computations. Experiments were carried out on a quad-core Intel i7 processor at 2.40 GHz.

Contractors are also implemented by IBEX. ( 7) is solved using the composition C f of a Newton iteration [START_REF] Hansen | Nonlinear Equations and Optimization[END_REF] (CtcNewton) and a Forward-Backward contractor (CtcFwdBwd) based on HC4revise algorithm [START_REF] Benhamou | Revising hull and box consistency[END_REF]. The velocity equation is implemented through symbolic differentiation of [START_REF] Müller | Redundant Actuation of Parallel Manipulators[END_REF], and related contractor C J is obtained through an interval Gauss-Seidel method (See [START_REF] Moore | [END_REF], Chapter 7). C f and C J are composed and iterated using CtcFixPoint contractor. Experiments proved those choices to be simple and time-efficient: more advanced algorithms as Acid or shaving operators did not improve the contraction.

B. Algorithm Behaviour

The first experiment consists in tracking a straight trajectory. It gives a good insight on the way the algorithm operates.

In this experiment, the desired trajectory crossing Type 2 singularity locus is directed along negative Y -axis of the plane represented in Fig. 2. Fig. 3 shows the end-effector pose enclosures [x] k . When approaching the singularity locus, the enclosures tend to become loose. For a short time after the crossing, the enclosures [x] k contains both solutions to the FK problem. Afterwards, the width of the boxes shrink rapidly due to an efficient contraction of the IA-based solver and the new assembly mode change is detected.

The fact that the assembly mode is detected can be explained from velocity information, as said in section II. Fig. 4 explains it: It displays the enclosure of the end-effector velocity ẏ along the Y -axis, and its integral

∆y(k) = kt s k i t s ẏ(t)dt (8)
where k i = 12000 is the time step when the crossing trajectory begins. ∆y represents the displacement that is undergone by the end-effector. The figure shows that at time step k = 12103, when the pose boxes shrink, the real displacement value lies in ∆y(k) = [-205.4, -57.0] (mm). This allows to eliminate the possible FK solution corresponding to assembly mode 1, whose estimated displacement is only -31.3 mm using nominal values for ξ ξ ξ . Fig. 5 shows the sizes of each interval during the crossing phase. It can be seen that all intervals grow before and right after the crossing. Still, a contraction occurs when the singularity is crossed, and the intervals locally shrink, most notably the interval [x]. The explanation is geometrical: in a singular configuration, the distal elements 12 and 22 are fully stretched along the x-axis, thus the interval on Out of 20 tries, the average execution time of the algorithm is 5.97 s for a 6.66 second long trajectory, on an off-line implementation. This execution time is encouraging for a real-time application.

C. Algorithm operating range

When a singularity is crossed successfully by the robot, there are two outputs possible for the algorithm: either it states that the assembly mode has changed, or it cannot state on the assembly mode. Of course, the second case does not give useful information on the behavior of the manipulator. During the tests, it has been observed that the algorithm often falls into the second case when playing slower trajectories.

To prove this, several trajectories are created and played by the robot. They all follow the same path as in the previous experiment and they all succeed in crossing the Type 2 singularity. However, their respective desired velocity at singularity v sing differs. The outputs for our experiment are the maximal widths of [x] and [ẋ] along the whole trajectory, respectively denoted by m x and m ẋ. They are shown in Table II. It is also mentioned if the algorithm was able to detect the change in assembly mode.

Results show that the faster the trajectory, the tighter the intervals, both for pose and velocity. It also tends to confirm the point that was made to conclude Section II-B: Faster trajectories are more likely to get their assembly mode change detected. The explanation is that at slow crossing speeds, more time is spent around the singularity, where the algorithm is unable to contract the intervals. Thus, [ẋ] grows at the rate k t s [a max ] so that [∆y] contains both solutions to FK problem, making it impossible to conclude on the assembly mode like it was done in section IV-B.

The behaviour of the algorithm has to be explained in the case where the new assembly mode is never detected. The result for the slowest trajectory is plotted on Fig. 6. It confirms that the boxes [x] grow up to the geometrical limits of the five-bar mechanism, always enclosing the two different solutions to the direct geometric problem. Boxes [ẋ] are not shown for concision matters, but they follow a simple behaviour in that case: they grow up to infinity at the rate a max , since no contraction is occuring anymore.

D. Reliability and eventual false positives

It was already mentionned that using valid bracketings for all the parameters, the algorithm should give reliable enclosures on x and ẋ. Two cases have been observed: either the enclosure only contains one solution to the FGM, or it contains both. The goal of this experiment is to prove that the enclosures remain reliable even in a case of failed crossing trajectory.

We define as false positive the fact that the algorithm states that a crossing attempt succeeded whereas it failed in reality. This is the worst situation that can exist, as the supervision of the robot will wrongly think that the robot is operating normally. To test this, we add an overconstraint on the torque output of the control scheme, despite using a trajectory that respects the crossing criterion. It results in an unsuccessful crossing attempt where the end-effector of the robot appears to bounce on the singularity locus.

The result is presented on Fig. 7, where the robot fails to cross the singularity. It can be seen that the algorithm keeps choosing the right solution to the FK problem (that is, the one above the singularity locus), i.e. it keeps track of the end-effector pose, thus certifying the reliability of our algorithm: no false positive was detected.

V. CONCLUSION

In order to enlarge the reachable workspace of parallel robots and to improve their reconfigurability, several methods have been developed. Using dedicated motion generators and control algorithms to cross Type 2 singularities is an interesting solution to this problem, as it can be widely generalized to all robots. However, the robot is not aware of which assembly mode it is currently in, limiting its autonomy. An algorithm that is able to solve this problem is proposed. As IA-based FK solvers cannot deal with singularity crossing, end-effector velocity tracking was introduced, allowing the algorithm to successfully detect assembly mode change. The reliability of the algorithm has been proven experimentally on a five-bar planar parallel robot. Its behaviour has been analyzed.

The experiments have notably shown that the output created by the algorithm in not useful when the crossing trajectory is too slow. This seems logical, as inertia is needed too in reality to achieve a successful crossing. Nonetheless, using the algorithm as it is remains possible in an industrial context, since robots are generally put to the limit of their acceleration capabilities.

Several perspectives for further investigation have been identified. An obvious improvement would be to implement this algorithm in real-time to detect assembly mode while a trajectory is playing on the robot. Reversing the inputs and outputs of the algorithm can help programming a robust motion generator that creates only trajectories which ensure a sucessful crossing of the Type 2 singularities. Both of these perspectives are under investigation.

Fig. 1 :

 1 Fig. 1: An example of crossing attempt. i is the initial configuration, f are the possible final configurations. The dotted line represents the singularity locus, the arrows the possible trajectories.

•

  Interval enclosures [x] 0 and [ẋ] 0 of the initial endeffector pose and velocity. • Upper bounds on the end-effector velocity and acceleration given under the form of interval enclosures [v max ] and [a max ] of all possible end-effector velocity and acceleration vectors. • An enclosure [ξ ξ ξ ] of the geometric parameters.

Fig. 2 :

 2 Fig. 2: DexTAR geometry

  -2 and 50 m • s -2 . The bounds [a max ] chosen are respectively [-80, 80] × [-80, 80], [-80, 80] × [-80, 80] and [-70, 70] × [-70, 70]. We have no information of the maximal speed, and have used [v max ] = [-∞, +∞] × [-∞, +∞] in all experiments, showing that this information is not crucial for the success of the IA-based tracking. Finally, the homing process of the DexTAR gives rise to [x] 0 = [0.07, 0.110] × [0.048, 0.088] (m) and [ẋ] 0 = [-0.1, 0.1] × [-0.1, 0.1] (m • s -1

Fig. 3 :

 3 Fig. 3: Pose boxes [x] along the trajectory, reliably enclosing the real pose of the end-effector

  (a) Velocity enclosure [ ẏ] along the trajectory. (b) Displacement enclosure [∆y] along the trajectory.

Fig. 4 :

 4 Fig. 4: Velocity information. The exact singular configuration is denoted by the solid black line.

  (a) Pose intervals widths. (b) Velocity intervals widths.

Fig. 5 :

 5 Fig. 5: Interval widths during singularity crossing. The exact singular configuration is denoted by the solid black line.

Fig. 6 :Fig. 7 :

 67 Fig. 6: Pose boxes [x] when the algorithm fails to detect assembly mode change

TABLE I :

 I Geometrical parameters values

	ξ i	unit	[ξ i ]	nominal
	11	mm	[89.832,89.995]	90
	21	mm	[89.929,90.115]	90
	12	mm	[89.899,90.050]	90
	22	mm	[89.877,90.084]	90
	d	mm	[117.843,118.107]	118
	δ 1	rad	[-0.004,0.004]	0
	δ 2	rad	[-0.004,0.004]	0

TABLE II :

 II Operation range experiment result

	v sing (m/s)	1.458 1.354	1.264	1.185	1.053
	m x (mm)	47.18 49.15	53.45	178.30	178.23
	m ẋ (m/s)	4.248 4.440	4.746	517.51	455.98
	AM detected?	yes	yes	yes	no	no