Rafael Balderas Hill

Sébastien Briot

Abdelhamid Chriette
email: abdelhamid.chriette@ls2n.fr

Philippe Martinet
email: philippe.martinet@inria.fr

Increasing energy efficiency of high-speed parallel robots by using variable stiffness springs and optimal motion generation

The classical approach to decrease the energy consumption of high-speed robots is by lowering the moving elements mass in order to have a lightweight structure. Even if this allows reducing the energy consumed, the lightweight architecture affects the robot stiffness, worsening the accuracy of the mechanism. Recently, variable stiffness actuators (VSAs) have been used in order to reduce the energy consumption of high-speed pick-andplace robots. The idea is to smartly tune online the stiffness of VSA springs so that the robot is put in near a resonance mode, thus considerably decreasing the energy consumption during fast pseudo-periodic pick-and-place motions. However, the serial configuration of springs and motors in the VSA leads to uncontrolled robot deflections at high-speeds and, thus, to a poor positioning accuracy of its end-effector.

In order to avoid these drawbacks and to increase the energy efficiency while ensuring the accuracy, this paper proposes the use of parallel arrangement of variable stiffness springs (VSS) and motors, combined with an energybased optimal trajectory planner. The VSS are used as energy storage for carrying out the reduction of the energy consumption and their parallel configuration with the motors ensure the load balancing at high-speed without los-ing the accuracy of the robot. Simulations of the suggested approach on a five-bar mechanism are performed and show the increase on energy efficiency.

INTRODUCTION

It is well-known that in industrial applications, such as high-speed pick-and-place operations, parallel robots are widely used [START_REF] Krut | A high-speed parallel robot for scara motions[END_REF]2]. Repeatability and accuracy are typically the most important criteria to measure their performance. Nevertheless, the design trends to operate at high speeds are shifting to the design of robots with lightweight architectures [3] in order to decrease the energy consumed by the motors, and measure as well the robot performance based on its energy efficiency [4].

For slow motions, gravity-balancing techniques [START_REF] Baradat | Design and Prototyping of a New Balancing Mechanism for Spatial Parallel Manipulators[END_REF][START_REF] Briot | A New Energyfree Gravity-compensation Adaptive System for Balancing of 4-DOF Robot Manipulators with Variable Payloads[END_REF][START_REF] Agrawal | Reactionless Space and Ground Robots: Novel Designs and Concept Studies[END_REF][START_REF] Wang | Static Balancing of Spatial Three-Degree-of-Freedom Parallel Mechanisms[END_REF] have been proposed in order to compensate the input efforts required to move the links of a pick-and-place robot, and thus to avoid consuming energy. Even if these methods have shown their effectiveness at slow speeds, it is not the case for high-speed operations in which the inertial effects are preponderant. A first solution introduced the series elastic actuators (SEAs) [START_REF] Pratt | Series elastic actuators[END_REF] to cope with the energy storage issues. The SEAs are compliant actuators composed by a motor which is linked to a spring in series that serves as energy storage, and whose stiffness is set by the spring constant. SEAs were first used to absorb contact shocks and to reduce the peak forces due to the impacts in bipedal walking robots [START_REF] Pratt | Virtual model control: an intuitive approach for bipedal locomotion[END_REF]. The limitation of the SEAs is that the stiffness is fixed and cannot be altered during motion, thus limiting the level of compliance to adapt for different tasks. Therefore, a recent second solution proposed the use of variable stiffness actuators (VSAs) [START_REF] Van Ham | Compliant actuator designs: review of actuators with passive adjustable compliance/controllable stiffness for robotic applications[END_REF][START_REF] Vanderborght | Variable impedance actuators: moving the robots of tomorrow[END_REF][START_REF] Wolf | Variable stiffness actuators: Review on design and components[END_REF] to handle with energy storage issues. VSAs consist of a motor which is connected to the output link by a spring in series and whose stiffness is variable and can be controlled. These springs are called Variable Stiffness Springs (VSS) [START_REF] Hollander | Concepts for compliant actuation in wearable robotic systems[END_REF]. Such actuators were first introduced for human-robot interaction applications [START_REF] Tsagarakis | A compact soft actuator unit for small scale human friendly robots[END_REF], for decreasing contact shocks and for efficiently actuate legged locomotion systems [START_REF] Van Ham | MACCEPA, the mechanically adjustable compliance and controllable equilibrium position actuator: design and implementation in a biped robot[END_REF][START_REF] Vanderborght | Overview of the lucy project: dynamic stabilization of a biped powered by pneumatic artificial muscles[END_REF]. Recently, they have been used for considerably decreasing the energy consumption of a pick-and-place robot in [START_REF] Uemura | Resonancebased motion control method for multi-joint robot through combining stiffness adaptation and iterative learning control[END_REF][START_REF] Goya | Realization of high-energy efficient pick-and-place tasks of scara robots by resonance[END_REF]. Even if the results seem to be impressive, their power consumption models do not consider the energetic losses from the motors that actuate the VSS, thus missing the losses of the full actuation chain. Furthermore, the serial arrangement of springs and motors would limit their operational speeds due to uncontrolled robot deflections when performing high-speed tasks.

A novel type of actuation for pick-and-place operations has been presented in [START_REF] Barreto | The concept of natural motion for pick and place operations[END_REF], where the authors proposed the use of constant stiffness springs in parallel to the motors in order to reduce the energy consumption. The main idea is to place the constant stiffness springs in parallel so that the required torque to move the output load is distributed between the motor and the parallel spring. This is possible since the constant stiffness spring in parallel supplies part of the required input effort to move the output load, thus reducing the energy consumption. Even if the results are impressive, the level of compliance to adapt to fast pseudo-periodic pick-and-place tasks is restricted due to the fixed stiffness.

The use of parallel elastic actuation seems to be a promising solution for increasing the energy efficiency. Therefore, the aim of this paper is to propose a strategy for reducing the energy consumed in high-speed pickand-place robots by studying the concept of using VSS in parallel to the motors, contrary to [START_REF] Uemura | Resonancebased motion control method for multi-joint robot through combining stiffness adaptation and iterative learning control[END_REF][START_REF] Goya | Realization of high-energy efficient pick-and-place tasks of scara robots by resonance[END_REF] in which they were mounted in series, as said previously. By placing the VSS in parallel, we will ensure direct power connection between the motor and links of the robot (and thus to have accuracy at high-speeds), and also to have tunable energy storage through the VSS. Additionally, contrary to what was done in [START_REF] Uemura | Resonancebased motion control method for multi-joint robot through combining stiffness adaptation and iterative learning control[END_REF][START_REF] Goya | Realization of high-energy efficient pick-and-place tasks of scara robots by resonance[END_REF], we deal with the minimization of the energetic losses from the full actuation chain, i.e., both the parallel robot joints motors and the variable stiffness springs motors, by finding the optimal robot trajectory thanks to an optimization problem.

This paper is divided as follows: Section 2 presents some recalls on the computation of the dynamic model of parallel robots. In addition to that, the dynamic model of parallel robots with VSS is shown. Then, the model of power losses is formulated considering the overall motor drive system. Section 3 shows the optimal problem formulation including the trajectory definition, cost function and decision variables. Section 4 presents the results of the proposed approach on a five-bar mechanism. In order to show the effectiveness of the approach presented in this paper, four cases are analyzed and compared: 1. Nominal case, which consists of the parallel robot without any compliant element, and with a classical fifth-degree polynomial trajectory;

2. Nominal case, with the proposed optimal trajectory planning algorithm;

3. Use of fixed stiffness springs in parallel with the actuated links of the robot, and with optimal motion planner;

4. Use of VSS in parallel with the robot with optimal motion planner.

Finally, the conclusions and future works on this topic are given.

PHYSICAL BACKGROUND

The paper aims at providing a strategy for increasing the energy efficiency of high-speed robots based on the concept of parallel VSS, it is thus necessary to make first some recalls on physical preliminaries of parallel robots, elastic actuation and motors modeling.

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 robot defined as the difference between the kinetic energy E and the potential energy U of the system:

q s j q i Motor Robot Link Motor VSS Robot Link
τ ta = d dt ∂L ∂ q T - ∂L ∂q T ( 4
)
w p = d dt ∂L ∂ ẋ T - ∂L ∂x T (5)
From equations (2) and (3), considering matrix A to be full rank, the dynamic model of a parallel robot is obtained:

τ = τ ta -B T A -T w p (6)
which, according to [START_REF] Briot | Dynamics of Parallel Robots: From Rigid Bodies to Flexible Elements[END_REF], can also be written under the form:

τ = Mq + h(q, q) ( 7
)
where M is an (n × n) definite positive matrix of inertia depending on the active joints coordinates q and platform coordinates x. h is an n s -dimensional vector of Coriolis, centrifugal and gravitational effects and its value depends on the active joint coordinates q, and their time derivatives q.

2.2 Modeling of VSS in parallel configuration with the actuated links

Let us consider a parallel robot with VSS in parallel to the actuated links (Fig. 1). According to Fig. 2, the actuation of the VSS is given by n s active joints. Their position and velocity are given by q s and qs representing the n s -dimensional vectors of variable stiffness joint variables and of variable stiffness joint velocities, respectively. By considering the effects of the elastic deformation of the springs with their force/displacement relations, the dynamics in [START_REF] Agrawal | Reactionless Space and Ground Robots: Novel Designs and Concept Studies[END_REF], become:

τ = Mq + h(q, q) + τ s ( 8
)
where τ s is the n s -dimensional vector of elastic torques associated to the VSS coupled to the robot in parallel:

τ s = K(q -q s ) (9)
where K is the (n s × n s ) stiffness matrix, and the dynamics of the VSS is expressed by:

τ vss = M s qs + h s (q s , qs ) -τ s ( 10
)
where M s is an (n s × n s ) definite positive matrix of inertia depending on the variable stiffness joints coordinates q s , and h s is an n s -dimensional vector of Coriolis, centrifugal and gravitational effects and its value depends on the variable stiffness joints coordinates q s and their time derivatives qs .

Power consumption model of the full actuation chain

In order to study the energy consumed by the motors during their different operation modes, several energy consumption models have been formulated, such as in [START_REF] Goya | Realization of high-energy efficient pick-and-place tasks of scara robots by resonance[END_REF].

For this work we consider direct-drive three phase brushless motors to actuate the joints of the parallel robot and the variable stiffness joints, therefore avoiding energetic losses due to gearboxes. The electromechanical schematics of such motor drive systems vary from manufacturer to manufacturer depending on the type of driver for energizing each phase of the brushless motor. However, five sections can be identified in their electromechanical sys-tem: (i) three phase rectification, (ii) braking, (iii) sixtransistor-commutable bridge, (iv) the three phase brushless motor, and (v) mechanical coupling between motor and drive shaft (for further details on electromechanical systems and power electronics the reader can refer to [START_REF] Chapman | Electric Machinery Fundamentals, Fifth Edition[END_REF]). Thus, to define the energy consumption model of such brushless motors, five types of power losses are classified:

• Resistive losses: These are heat losses that occur when current flows through the resistances of each phase of the brushless motors and the braking resistances. The mathematical expression is given by:

P motor = R phase I 2 (11)
P brake = R brake I 2 ( 12
)
where I is the line-to-line current between each phase of the motor, and R phase and R brake are electrical resistances.

• Damping losses: These are mechanical losses due to the mechanical coupling between the motor and drive shaft modeled by:

P damp = B q2 k ( 13
)
where B is the damping coefficient and qk is the velocity of the kth motor axis, where k = 1, ..., (n + n s ).

• Conduction losses: The conduction losses are electrical losses due to the current flow through the transistors and diodes occurring in the three phase rectification and the six-transistor-commutable bridge. They are modeled as follows:

P cond = u ce I cav + r c I 2 crms transistors + u D I Dav + r D I 2 Drms diodes (14)
where u ce is the collector-to-emitter voltage and I cav is the collector current, r c is the on-state resistance and I crms is the RMS value of the collector current. u D is the diode forward voltage and I Dav is the diode current, r D is the on-state resistance and I Drms is the RMS value of the diode current.

• Switching losses: These losses occur due to the continuous switching in the six-transistor-commutable bridge of the motor driver when energizing each phase of the motor. The relationship is given by:

P switch = u ce I cav f sw transistors + u D I Dav f sw diodes ( 15
)
where f sw is the switching frequency from the sixtransistor-commutable bridge.

• Rectification losses: These are electrical losses generated from the power supply of three phase rectification modeled as follows:

P rectif ier = u DR I DRav + r DR I 2 DRrms ( 16
)
where u DR is the diode rectifier forward voltage and I DRav is the diode rectifier current, r DR is an onstate resistance and I DRrms is the RMS value of the diode rectifier current.

It should be mentioned that the currents I cav , I crms , I Dav , I Drms , I DRav and I DRrms , associated to the conduction, switching and rectification losses due to the diodes and transistors, are functions of the line-to-line currents I of the motor and thus, by using the model of the brushless motor in [START_REF] Pillay | Modeling, simulation, and analysis of permanent-magnet motor drives. Part II. The brushless DC motor drive[END_REF], and functions of the motor electromagnetive torque τ e , and input torque τ . The expressions linking all these variables are not given here due to page limitations, but can be found in [START_REF] Millman | Microelectronics: Digital and Analog Circuits and Systems[END_REF]. Finally, by grouping the aforementioned power losses, the following model of losses results:

P losses = P motor + P brake + P damp

P cond + P switch + P rectif ier (17)
By investigating the mathematical formulations (modeling the dynamics of parallel robots and the power losses of their motor drive systems) it is not evident to see how to minimize the energy consumption. For this reason, it is thus necessary to find a strategy to exploit the power consumption model in order to minimize the power losses during the motion.

3 Optimal problem formulation for minimizing the energetic losses on parallel robots

In order to cope with an optimal approach to minimize the energetic losses presented in section 2, an optimal trajectory planning algorithm is designed. The main idea is to define a motion generator such that the energetic losses, computed from numerical integration of ( 17), are minimized, i.e. to find an optimal trajectory q(t) (or x(t) if expressed in Cartesian space) and optimal spring coordinates q s (t) that minimizes the energetic losses during the motion.

For a pick-and-place manipulator, there is no restriction for the robot on how to go from the initial to the final positions, except eventual collisions within the robot links, external environment or singularity loci. That is why, similar to what it was done in [START_REF] Saravana Perumaal | Automated trajectory planner of industrial robot pickand-place task[END_REF], our trajectory planning algorithm finds a path sequence of via points (intermediate points) between the initial and final positions (Fig. 3 and4), which are then interconnected by the motion profiles defined thereafter.

Our proposed energy-based optimal trajectory planner is thus an algorithm which minimizes an objective function that computes the energy consumed during the motion by finding optimal sequence of via points between initial and final robot positions. In order to do so, we first need to define a trajectory planner that fully exploits the dynamic performances of the modified high-speed parallel robot so that the cost function based on the power losses can be minimum. It is described thereafter.

Trajectory definition

In order to develop the motion planners, we decided to denote as M s the total number of via points along the variable stiffness joints trajectory and as M the total number of via points along the parallel robot joints trajectory. The total duration of the trajectory is denoted as t f and we impose the duration of the motion between two successive via-points to be equal to δt s = t f /(M s + 1) for the variable stiffness joints trajectory and to δt = t f /(M +1) for the robot joint trajectory.

Variable stiffness joints trajectories

Different to what it has been done in [START_REF] Van Ham | Compliant actuator designs: review of actuators with passive adjustable compliance/controllable stiffness for robotic applications[END_REF][START_REF] Wolf | Variable stiffness actuators: Review on design and components[END_REF], in which the authors proposed to find motion profiles for the stiffness rather than the variable stiffness joints, this paper generates directly the joint trajectories for the motors that actuate the VSS. This avoids to handle further computations for calculating the stiffness and allows the use of a classical position controller, avoiding the issue of measuring the stiffness in the feedback. In order to generate the variable stiffness joint motion planner let us express q s from (8) as follows:

q s = -K -1 (τ -Mq -h(q, q)) + q (18)
Based on expression [START_REF] Uemura | Resonancebased motion control method for multi-joint robot through combining stiffness adaptation and iterative learning control[END_REF], it is possible to observe that the variable stiffness joint dynamics in [START_REF] Pratt | Virtual model control: an intuitive approach for bipedal locomotion[END_REF], and therefore the optimal spring stiffness, is dependent, not only on q s and its time derivatives qs and qs but also on the configuration of the parallel robot joints q, and their velocities q and accelerations q.

We build the variable stiffness joint motion planner as follows. For the jth component of the vector q s denoted as q sj , let us consider a sequence of two successive viapoints q sj k = q sj (t s,k ) and q sj k+1 = q sj (t s,k+1 ) on the time interval t ∈ [t s,k , t s,k+1 ] with t s,k = t 0 + k δt s (k = 0, ..., M s ) and t s,k+1 = t s,k + δt s (one of the viapoints can be replaced by the initial or final point of the trajectory). Taking advantage of model [START_REF] Uemura | Resonancebased motion control method for multi-joint robot through combining stiffness adaptation and iterative learning control[END_REF], we connect the via-point q sj k and q sj k+1 (Fig. 4) thanks to the following expression defined by means of several polynomials, valuable on the time interval t ∈ [t s,k , t s,k+1 ]:

q sj (t) = p (k) 1j (q j (t)) + p (k) 2j ( qj (t)) + p (k) 3j (q j (t)) + p (k) 4j (t) (19) where p (k) 1j (q j (t)) is a polynomial of q j (t), p (k) 2j ( qj (t)) is a polynomial of qj (t), p (k) 3j (q j (t)
) is a polynomial of qj (t), and finally p (k) 4j (t) is a polynomial of time that is used to model the effects of τ on q sj . The motion profile used for q j (t) will be defined later in Section 3.1.2.

In order to obtain the explicit form of each polynomial, we follow a heuristic approach of experimentally testing different functions until obtaining the best polynomial form based on the energy consumed during the motion. After trials and errors, we found the best results q j s f q j0 q j k q j k+1 q j M q j f when:

Start

p (k) 1j (q j (t)) = 0, p (k) 2j ( qj (t)) = 0 (20)
p (k) 3j (q j (t)) = K jk q3 j + K j(k+1) q2 j + K j(k+2) qj + K j(k+3) (21)
p (k) 4j (t) = (q sj k+1 -q sj k )(α 0sjk + α 1sjk t + α 2sjk t 2 + α 3sjk t 3 + α 4sjk t 4 + α 5sjk t 5 ) (22)
Parameters (K jk , ..., K j(k+3) ) are set as decision vari-ables of the optimization problem, as well as the parameters q sj k , v sj k and a sj k for k = 1, ..., M s . Note that, in our case, in order to simplify the resolution of the optimization problem, we decide to use the same values of the parameters (K jk , ..., K j(k+3) ) for any segment between two consecutive via-points.

The boundary and continuity conditions between each portion of motion are defined as follows:

                  
q sj (t s,k ) = q sj k and q sj (t s,k+1 ) = q sj k+1 qsj (t s,k ) = v sj k and qsj (t s,k+1 ) = v sj k+1 qsj (t s,k ) = a sj k and qsj (t s,k+1 ) = a sj k+1 q sj (t 0 ) = q sj0 and q sj (t f ) = q sj f qsj (t 0 ) = 0 and qsj (t f ) = 0 qsj (t 0 ) = 0 and qsj (t f ) = 0 [START_REF] Chapman | Electric Machinery Fundamentals, Fifth Edition[END_REF] where v sj k and a sj k (v sj k+1 and a sj k+1 , resp.) represent the velocity and acceleration of the variable stiffness joint j at via-point k (k + 1, resp.). For a given set of parameters v sj k , v sj k+1 , a sj k , a sj k+1 , (K jk , ..., K j(k+3) ) and robot joint trajectory q j (t), the boundary and continuity conditions [START_REF] Chapman | Electric Machinery Fundamentals, Fifth Edition[END_REF] are used in order to compute the coefficients α lsjk in [START_REF] Merlet | Parallel Robots[END_REF].

Parallel robot trajectory generation

Since our motion planner for q s is a function of q and because Eq. ( 10) depends on qs , it is necessary to constrain jerks and the derivative of the jerks of q. For the jth component of the vector q denoted as q j , let us consider a sequence of two successive via-points q j k = q j (t k ) and q j k+1 = q j (t k+1 ) on the time interval t ∈ [t k , t k+1 ] with t k = t 0 + k δt (k = 0, ..., M ) and t k+1 = t k + δt (one of the via-points can be the initial or final point of the trajectory). The trajectory for q j between via-points q j k and q j k+1 is thus defined by using ninth-degree timedependent polynomials which can fix the position, velocities, accelerations, jerks and derivatives of the jerks at both via-points. The expression of the polynomial on the time interval t ∈ [t k , t k+1 ] is given by:

q j (t) = (q j k+1 -q j k )(α 0jk + α 1jk t + α 2jk t 2 + α 3jk t 3 + α 4jk t 4 + α 5jk t 5 + α 6jk t 6 + α 7jk t 7 + α 8jk t 8 + α 9jk t 9 ) (24)
The boundary and continuity conditions are defined as follows:

                  
q j (t k ) = q j k and q j (t k+1 ) = q j k+1 qj (t k ) = v j k and qj (t k+1 ) = v j k+1 qj (t k ) = a j k and qj (t k+1 ) = a j k+1 q j (t 0 ) = q j0 and q j (t f ) = q j f qj (t 0 ) = 0 and qj (t f ) = 0 qj (t 0 ) = 0 and qj (t f ) = 0 [START_REF] Millman | Microelectronics: Digital and Analog Circuits and Systems[END_REF] and q

[3]

j (t k ) = 0 and q [4]
k (t k ) = 0 at any via-point or at the initial and final points. In [START_REF] Millman | Microelectronics: Digital and Analog Circuits and Systems[END_REF], v j k and a j k (v j k+1 and a j k+1 , resp.) represent the velocity and acceleration of the robot joint joint j at via-point k (k + 1, resp.). For given values of parameters v j k , v j k+1 , a j k , a j k+1 , the boundary and continuity conditions [START_REF] Millman | Microelectronics: Digital and Analog Circuits and Systems[END_REF] are used in order to compute the coefficients α ljk in [START_REF] Pillay | Modeling, simulation, and analysis of permanent-magnet motor drives. Part II. The brushless DC motor drive[END_REF].

Parameters q j k , v j k and a j k for k = 1, ..., M are set as decision variables of the optimization problem.

Optimization Problem Formulation: Cost function and decision variables

The optimization problem aims to minimize the energetic losses of the modified parallel robot that considers the full actuation chain, taking into account physical constraints of the robot and the motors. The cost function of the optimization problem corresponds to the power losses of the motors that actuate the parallel robot joints and the variable stiffness joints since the manipulator must be as energy efficient as possible. The cost function for the optimizing process is thus defined as follows:

J = t f 0 n+ns k=1 P losses k dt ( 26
)
where P losses k represents the power losses of the kth motor of the full actuation chain from expression [START_REF] Vanderborght | Overview of the lucy project: dynamic stabilization of a biped powered by pneumatic artificial muscles[END_REF].

Based on the description of the trajectory planner in the previous sections, the decision variable vector of the optimization problem, denoted as x DecVar , contains all aforementioned parameters q j k , q sj k , v j k , v sj k , a j k , a sj k and (K jk , ..., K j(k+3) ).

The optimal trajectory through via-points problem aims at finding the decision variable vector x DecVar that minimizes the power losses of the (n + n s ) motors while re-specting physical constraints, such as motor limits and collisions. In this paper, we handle the collision avoidance of virtual obstacles, which represents singularity loci of the parallel robot. This to prevent of finding via-points in the neighborhood of the parallel robot singularities.

The optimization problem was solved by means of MATLAB f mincon function and it is formulated as follows:

minimize J over x DecVar subject to |τ e k | ≤ τ e max |κ| < (27)
where τ e k is the electromagnetive torque of the kth motor and κ represents the condition of proximity to the singularity with denoting a threshold tuned experimentally. For a detailed overview on singularities of 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].

CASE STUDY

Five-bar mechanism with variable parallel elastic actuation

The proposed optimization approach for increasing the energy efficiency was validated in a five-bar mechanism with variable stiffness torsional springs in parallel, as shown in Fig. 5. The actuation of the parallel robot 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 actuation for the variable stiffness joints is given by q s = [q s1 , q s2 ] T . The numerical values of the five-bar mechanism geometric and dynamic parameters were defined according to a real parallel robot called RobEcolo in [START_REF]Project RobEcolo[END_REF] [30], and have the following values:

• inertia of the two proximal links zz 11 = 0.133 kg/m 2 and zz 21 = 0.122 kg/m 2 ;

• Coulomb friction terms of the actuated joints f s1 = 0.659 Nm and f s1 = 0.225 Nm;

• end-effector mass m = 1.84 kg;

• inertias of the couplings between the motors and springs J 1 = J 2 = 0.0051 kg/m 2 ;

• stiffness constants of the springs k 11 = k 21 = 4 Nm/rad;

• link lengths of the two proximal links l 1 = l 2 = 0.28 m;

• link lengths of the two distal links l 3 = l 4 = 0.4 m;

• distance between the axes of the actuated links a = 0.250 m.

The motor parameters used to examine the energetic performance with the proposed approach were defined according to the motor specifications found in [START_REF] Siemens | The performance-optimized and easy-touse servo drive system[END_REF] by the following values:

• electrical resistances R phase = 0.2 Ω, R brake = 18 Ω;

• damping coefficient B = 0.005 Nm/s;

• on-state resistances r c = r D = r DR = 0.003 Ω;

• forward and collector-to-emitter voltages u D = u DR = u ce = 0.8 V;

• switching frequency f sw = 20 KHz.

The algorithm validation was performed by using the dynamic model of the real parallel robot called RobEcolo [START_REF]Project RobEcolo[END_REF] [30], which is a wooden five-bar mechanism whose dynamic model is given by the following equations:

τ = ZZq + B T λ + f s1 sign( q1 ) f s2 sign( q2 ) + τ s (28)
w p = A T λ = mẍ (29)
where A and B are computed from the input-output kinematic constraint relation (1), and τ s is calculated from [START_REF] Pratt | Series elastic actuators[END_REF]. ZZ = diag(zz 11 , zz 21 ), with zz 11 and zz 21 representing the inertia of the two proximal links, and K = diag(k 11 , k 21 ), with k 11 , k 21 being the springs constants. Finally, the VSS dynamics from ( 10) is given by:

τ vss = J 1 0 0 J 2 qs1 qs2 - k 11 0 0 k 21 (q 11 -q s1 ) (q 21 -q s2 ) τs ( 30
)
q s2
q 21 q 11 q s1 q 22 q 12 q 13 a l l l l1 2 3 4

Figure 5: FIVE-BAR MECHANISM PARAMETERIZA-TION WITH TWO ACTUATED JOINTS q 11 AND q 21 , AND THREE PASSIVE JOINTS q 12 , q 22 AND q 13 . THE VARIABLE STIFFNESS TORSIONAL SPRINGS ARE LOCATED IN PARALLEL TO THE TWO ACTUATED JOINTS DEFINED BY q s1 AND q s2 .

Results and comparative energetic analysis on three cases: nominal, constant stiffness and variable stiffness actuation

In order to validate the theoretical formulations, the optimization algorithm is tested by defining multiple desired points in the workspace of the five-bar mechanism as shown in Fig. 6. The aim is to go multiple times to these points by finding the appropriated optimal decision variable vector x DecVar that minimizes the energetic losses of the full actuation chain. The sections of trajectories are defined between the points given by: A The sequence of the trajectory is defined at Fig. 6. For showing the effectiveness of the proposed approach, the input torques and the energy consumption of four different types of actuation are analyzed: (i) Nominal actuation, in which there is no elastic element attached at each joint, and a fifth-degree polynomial trajectory is used between each set of two successive points; (ii) nominal actuation with optimized trajectory; (iii) actuation with fixed stiffness springs in parallel with each joint of the robot, with spring constants k 11 and k 21 defined according to section 4.1, and with optimized trajectory; (iv) actuation with VSS in parallel with each actuated link of the robot and with optimized trajectory. Figures 7 and8 show the input torques computed for the four types of actuation in the parallel robot joints along the multiple trajectories from the sequence in Fig. 6 with travel times of 0.2 s between each desired set of two points. It can be seen that in both parallel robot actuated joints, the input torques are reduced by using springs, with fixed stiffness and then with variable stiffness, being the latter the type of actuation with major gain in reduction of input efforts. The torques shown in Fig. 9 correspond to the input torques of the variable stiffness joints. As compared to the parallel robot joint torque magnitudes, the variable stiffness joint torques do not have a major impact on input efforts. Then, Fig. 10 shows the energetic losses from the power consumption model [START_REF] Vanderborght | Overview of the lucy project: dynamic stabilization of a biped powered by pneumatic artificial muscles[END_REF] for the four types of actuation. It can be seen that by using VSS, it is possible to reduce the energy up to 48 % with respect to the nominal case with classical fifth-degree polynomial trajectory. It should be mentioned that this result is far from the gain in energy efficiency of 90 % shown in [START_REF] Goya | Realization of high-energy efficient pick-and-place tasks of scara robots by resonance[END_REF], but this is due to the fact that, in our work, we did not neglect any losses in the actuation chain including the VSS and the motor driver, which was not the case at all in [START_REF] Goya | Realization of high-energy efficient pick-and-place tasks of scara robots by resonance[END_REF].

In addition to that, Fig. 11 shows the energy consumption reduction achieved by using VSS for different pickand-place travel times between the same points. It can be seen that our algorithm has better performance when performing the most high-speed tasks. In Fig. 12, the influence of modifying the number of via-points for the parallel robot and variable stiffness joints is shown. It is observed that the major reduction occurs when the number of via-points is set to two for both the parallel robot joints and variable stiffness joints. This is due to the fact that increasing the number of via-points and having fast pick-and-place travel times complexified the resolution of the optimization algorithm and made it converged towards local minima. Finally, Fig. 13 shows the energy reduction separately for each identified power loss presented in section 2.3. It can be seen that the major reduction occurs on the resistive losses, which is expected since they are affected by the squared of the electric currents, which are functions of the electromagnetive torques, included as a contraint in the optimization problem.

→ B → C → D → E → F → G → B → A,

CONCLUSIONS

This paper proposes a new approach to increase the energy efficiency of high-speed pick-and-place parallel robots by placing variable stiffness springs in parallel with the motors. An analysis of the energetic losses of the full actuation chain provided a cost function which is used in an optimal problem formulation in order to minimize the power losses. Then, specific trajectories for the parallel robot joints and for the variable stiffness joints were de- nates online, avoiding the pre-computation of the motion planners. [2] Angeles, J., Caro, S., Khan, W., Morozov, A., 2006. Kinetostatic design of an innovative schoenfliesmotion generator. In Proc. of IMechE Part C: Journal of Mechanical Engineering Science, 220 (7), pp. 935-943.

[3] Clavel R., 1990. Device for the movement and positioning of an element in space. US Patent 4 976 582, December 11.

[4] Brossog, M., Bornschlegl, M., Franke, J., 2015. Reducing the energy consumption of industrial robots in manufacturing systems. International Journal of Advanced Manufacturing Technology, Vol. 78 Issue 5-8, pp. 1315-1328.

Figure 1 :

 1 Figure 1: ON THE LEFT, A FIVE-BAR MECHANISM WITH VARIABLE STIFFNESS LINEAR SPRINGS IN PARALLEL ARRANGEMENT WITH THE AC-TUATED LINKS IS SHOWN. ON THE RIGHT THE SAME PARALLEL ROBOT WITH VARIABLE STIFF-NESS TORSIONAL SPRINGS IS PRESENTED (THE GREY CIRCLES DENOTE THE ACTUATED JOINTS).

Figure 2 :

 2 Figure 2: POWER TRANSMISSION SYSTEM OF VARIABLE STIFFNESS SPRINGS IN PARALLEL TO THE MOTORS. q i AND q sj REPRESENT THE PAR-ALLEL ROBOT JOINTS AND VARIABLE STIFFNESS JOINTS COORDINATES, RESPECTIVELY, AND i = 1, ..., n, j = 1, ..., n s .

Figure 4 :

 4 Figure 4: ON THE TOP A VIA-POINT TRAJECTORY FOR THE PARALLEL ROBOT JOINTS WITH M VIA-POINTS. ON THE BOTTOM A VIA-POINT TRAJEC-TORY FOR THE VARIABLE STIFFNESS JOINT MO-TION PLANNER WITH M S VIA-POINTS.

  = [0.25, 0.3], B = [0, 0.4], C = [0.2, 0.5], D = [0, 0.45], E = [0.2, 0.4], F = [0, 0.5], G = [0.2, 0.35].

Figure 6 :

 6 Figure 6: MULTIPLE OPTIMAL TRAJECTORIES FOR THE SEQUENCE: A → B → C → D → E → F → G → B → A, WITH TRAVEL TIMES OF 0.2 s BETWEEN DESIRED POINTS.

Figure 10 :

 10 Figure 10: ENERGETIC LOSSES FOR THE FOUR CASES ANALYZED: NOMINAL FIFTH-DEGREE POLYNOMIAL, NOMINAL WITH OPTIMIZED TRA-JECTORY, FIXED STIFFNESS AND VARIABLE STIFFNESS, 0.2 s AS PICK-AND-PLACE TRAVEL TIME BETWEEN DESIRED POINTS.

Figure 11 :

 11 Figure 11: REDUCTION OF ENERGY CONSUMP-TION FOR DIFFERENT PICK-AND-PLACE TRAVEL TIMES BETWEEN DESIRED POINTS WHEN USING VSA.

Figure 12 :

 12 Figure 12: REDUCTION OF ENERGY CONSUMP-TION FOR DIFFERENT NUMBER OF VIA-POINTS WITH 0.2 s AS PICK-AND-PLACE TRAVEL TIME BETWEEN DESIRED POINTS.

Figure 13 :

 13 Figure 13: ENERGY REDUCTION FOR EACH IDEN-TIFIED POWER LOSS.

  Figure 9: INPUT TORQUES FOR THE JOINTS OF VARIABLE STIFFNESS, 0.2 s AS PICK-AND-PLACE TRAVEL TIME BETWEEN DESIRED POINTS.veloped avoiding to have a nonlinear force/displacement relation for the elastic torques. Simulations led to the successful increasing of energy efficiency with multiplepoint trajectories on a five-bar mechanism for four cases: nominal with a classical fifth-degree polynomial, nominal with optimized trajectory, use of fixed stiffness and variable stiffness springs. Results show that the reduction of energy consumption can reach up to 48 % when using VSS. Next work on this subject includes experimental validation of our approach and the design of an alternative optimization strategy in order to tune the spring coordi-

					Nominal Fifth-deg						Nominal Fifth-deg
	Energetic Losses (J)	0 0 2000 4000 6000 8000	0.2	0.4	0.6 Time (seconds) 0.8 1 Nominal Optimized 1.2 Fixed Stiffness Variable Stiffness	1.4	1.6	Torque Joint 2 (Nm)	-40 -20 0 20 40 60		Nominal Optimized Fixed Stiffness Variable Stiffness
									0 -60	0.2	0.4	0.6	0.8	1	1.2	1.4	1.6
											Time (seconds)
								Figure 8: INPUT TORQUES FOR THE SECOND
								JOINT OF THE PARALLEL ROBOT FOR THE FOUR
								CASES: NOMINAL FIFTH-DEGREE POLYNOMIAL,
								NOMINAL WITH OPTIMIZED TRAJECTORY, FIXED
								STIFFNESS AND VARIABLE STIFFNESS, 0.2 s AS
								PICK-AND-PLACE TRAVEL TIME BETWEEN DE-
								SIRED POINTS.
									Elastic Torques (Nm)	0 5	Variable Stiffness Joint 1 Variable Stiffness Joint 2
									0 -5	0.2	0.4	0.6	0.8	1	1.2	1.4	1.6
											Time (seconds)

Acknowledgment

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


Lila Kaci
email: lila.kaci@ls2n.fr

Sébastien Briot
email: sebastien.briot@ls2n.fr

Clément Boudaud
email: clement.boudaud@ecoledubois.fr

Philippe Martinet
email: philippe.martinet@ls2n.fr

RobEcolo: OPTIMAL DESIGN OF A WOODEN FIVE-BAR MECHANISM

Eco-design of robots has almost never be en explored in the past. This work investigates the potential of using bio-sourced materials, which have almost no environmental impact, instead of metals for robot design. Also, wood is one of the best candidates because of its interesting mechanical properties. However, wood performance / dimensions vary with the atmospheric conditions / external solicitations. Thus, it is challenging to design a stiff and accurate wooden industrial robot.

Therefore, the objective of this paper is to describe a new design methodology leading to the design a wooden five-bar mechanism reliable in terms of accuracy and stiffness. The design optimization problem is solved in cascade. The first optimization process proposes to use a control-based design approach in order to compute the optimal primary geometric parameters of the robot (lengths of the links). This approach takes into account the sensor-based controller performance during the design phase. The second optimization process deals with the issue of the variability of the wood mechanical performance. It is based on a reliable topology optimization approach and allows for finding the shape of the robot links for which the impact of this variability in terms of deformation is minimized. Theoretical developments are described, solved and the obtained results allowed the prototyping of an industrial wooden five-bar mechanism.

INTRODUCTION

Robots considerably impact the environment, not only during their use but also during their manufacturing phase which represents almost 50 % of the impact [1]. The impact of the robot manufacturing is due to the use of polluting materials such as, for the design of the mechanical architecture, metals (iron or aluminum alloys) or even carbon composites. However, as shown during COP21, one of the priority in the world is the Climate Change Mitigation. Thus, it is very urgent to find new solutions to design machines that limit the Environmental Impact (EI) of the human industrial activities. Robotics must also participate to this effort and robot designers should investigate the field of robot eco-design.

In order to decrease the EI during the design of the robot architecture, bio-sourced materials (BSM) could be used. BSM have a little (even no) EI [START_REF] Thakur | Green Composites from Natural Resources[END_REF][START_REF] Laurent | Vers un robot en bois: Première partie[END_REF]. In the past, wood was used in machine design, (e.g. planes up to Wold War II, cars and aircraft chassis [START_REF] Sharp | Mosquito[END_REF]). Nowadays, wood is still popular in building design [START_REF] Fleming | Measuring-up in timber: a critical perspective on mid-and high-rise timber building design[END_REF] due to its high mechanical stiffness and its low cost. In Robotics, the use of this material is limited to design mock-ups and prototypes [START_REF] Roennau | Design and kinematics of a biologicallyinspired leg for a six-legged walking machine[END_REF][START_REF] Quigley | A low-cost compliant 7-DOF robotic manipulator[END_REF]. A first attempt to use wood in robot design was presented in [START_REF] Laurent | Vers un robot en bois: Première partie[END_REF]. The obtained results validated the feasibility of designing industrial wooden robot. However, the wood was not used to guarantee the robot performance in terms of stiffness and accuracy. This is a crucial issue because wood is an organic material, and it presents a variability in its performance / dimensions with respect to atmospheric conditions / external solicitations. Its performance also depends on the conditions in which trees have grown [START_REF] Kretschmann | Ch. 5: Mechanical properties of Wood[END_REF][START_REF] Wang | From hydrophilicity to hydrophobicity: a critical review[END_REF].

In this paper we introduce a design methodology well adapted in order to design a stiff and accurate wooden industrial robot. This methodology is applied on a five-bar mechanism, i.e. a parallel robot with two translational degrees of freedom. This methodology is based on the following concepts:

• In order to ensure the robot accuracy even with the variability in the wooden link dimensions, sensor-based controllers could be used [START_REF] Chaumette | Handbook of robotics[END_REF]. However, as shown in [START_REF] Kaci | Control-based design of a five-bar mechanism[END_REF], positioning accuracy of sensor-based controllers (obviously) depends of the number and types of sensors that are used for the robot control, but also, in certain cases, of the robot geometric parameters [START_REF] Briot | The hidden robot: an efficient concept contributing to the analysis of the controllability of parallel robots in advanced visual servoing techniques[END_REF]. Therefore, the robot design process must take into account performance criteria coming from usual physical properties [START_REF] Angeles | 10: Performance Evaluation and Design Criteria[END_REF] but also from the type of controller that will be implemented on the robot. This approach is quite unusual for industrial robot design.

• In order to ensure the robot stiffness performance even with the variability of the wood mechanical properties, link shapes can be optimized thanks to reliable topology optimization algorithms [START_REF] Asadpoure | Robust topology optimization of structures with uncertainties in stiffness -Application to truss structures[END_REF].

This methodology has been defined in the context of the project RobEcolo [15] whose goal is to show that the EI of robots can be considerably reduced by using BSM instead of metals or carbon composites. This paper is divided as follows. In Section 2, the robot design requirements are specified. Section 3 introduces the sensor-based control approach that will be used to control the robot motions and presents the design methodology adapted to optimize the primary geometric parameters of the wooden five-bar mechanism taking into account the controller characteristics. In Section 4, a reliable topology optimization approach able to find the robot link shapes is proposed. Section 5 will describe the CAD design and manufacture of the prototype of the wooden fivebar mechanism. Finally, in section 6, conclusions are drawn.

ROBOT ARCHITECTURE AND SPECIFICATIONS

The five-bar mechanism geometry [START_REF] Campos | Development of a five-bar parallel robot with large workspace[END_REF] is illustrated in Fig. 1. It is a two degrees-of-freedom (dof ) planar parallel robot with two translations in the plane (x 0 Oy 0 ). The end-effector is positioned at point C.

The mechanism consists of four links: proximal links are the bodies between A i and B i (i = 1,2) while the bodies between points B i and C are called distal links. The proximal and distal links are connected by three passive revolute joints, the joints located at point A 1 and A 2 are active. The requirements that the wooden five-bar mechanism should satisfy have been fixed by our partners in the project RobEcolo [15]. These requirements are given in Tab. 1. These performance should be valid wherever in the robot workspace. Thus, of course, the workspace should be free of singularity. In addition, the robot must be as compact as possible.

Deformations along z 0 have not been specified because the robot purpose is to have a precise positioning in its plane of motion. However, these vertical deformations will be checked a posteriori in the design process in order to verify that they are "acceptable", which means for the project partners, around The robot is a prototype made for validating performance in terms of stiffness and accuracy, therefore no dynamic properties have been imposed. However, the design process should tend to reduce the mass or robot link inertia so that it avoids leading to the manufacture of bulky links.

Finally, the number and types of sensor has been imposed by the project partners. With respect to the desired 500 microns of absolute positioning accuracy specified in Tab. 1, we propose to use four cameras Toshiba Bu 238 M (1920 × 1200 pixels of resolution and a focal length of 8 mm) as sensors. With lens of 8 mm, distortion will be limited and the global visual system (4 cameras) will require to be calibrated for the real implementation. Then, to control the robot, servoing approaches will be used. These approaches are relevant with respect to the accuracy to be reached [START_REF] Chaumette | Handbook of robotics[END_REF].

Based on all these requirements, it is possible to see that the problem can be solved in cascade. Indeed, the absolute positioning accuracy requirements essentially depend on the robot link length and camera performance / position, while the deformations is more related to the shape of the robot links and material properties. Therefore, the optimization problem will be separated into two independent parts: the first one will deal with the definition of the primary geometric parameters, i.e. the link lengths and camera positions/orientations, under robot positioning error specifications; the second one will deal with the definition of the optimal link shapes under elastic performance constraints.

The next Section introduces the first optimization problem taking into account the characteristics of the controller and of the perception system.

CONTROL-BASED DE-SIGN

Brief recalls on leg-based visual servoing

As mentioned above, in order to ensure the robot positioning accuracy even with the dimensional variability of the wooden links, visual servoing will be used. The usual approaches are based on the direct observation of the end-effector pose [START_REF] Traslosheros | Visual servoing of a parallel robot system[END_REF]. In our case however the end-effector of a five-bar mechanism is a point, which is difficult to observe. This is why we decided to observe the distal links, whose dimensions are more compatible with an observation by a camera, as it was proposed in [START_REF] Andreff | Image-based visual servoing of gough-stewart parallel manipulators using legs observation[END_REF][START_REF] Andreff | Visionbased kinematic modelling of some parallel manipulators for control purposes[END_REF]. The scheme of a general vision-based control based on the observation of the robot links is presented in Fig. 2. However, it is known from [START_REF] Briot | The hidden robot: an efficient concept contributing to the analysis of the controllability of parallel robots in advanced visual servoing techniques[END_REF] that these types of controllers present singularities, that affect their accuracy performance. These singularities depend not only on the choice of the sensor localizations but also on the robot dimensions. Thus, it is highly necessary to take into account all these parameters in the design process.

In order to facilitate the observation of the distal links, in the following of the paper, their shapes are imposed to be cylindrical [START_REF] Andreff | Image-based visual servoing of gough-stewart parallel manipulators using legs observation[END_REF]. Indeed, as shown in Fig. 3, cylinder projections in the camera plane are two lines whose interpretation planes (i.e. the planes passing through the camera center and the lines) have as normal vectors c ṅ1 i and c ṅ2 i (the superscript c denotes that these vectors are expressed in the camera frame). Their expressions depend on the cylinder size and location. These two vectors are indeed the measurements that can be extracted from the camera observation. By knowing them and using the expressions given in [START_REF] Andreff | Image-based visual servoing of gough-stewart parallel manipulators using legs observation[END_REF], it is then possible to reconstruct the Plücker coordinates ( c u i , c h i ) of the line L i passing trough the cylinder axis and their derivatives with respect to time as:

c uT i c ḣT i T = H T i c ṅ1 i c ṅ2 i , i = 1, 2 (1)
where H T i is a (6 × 6) matrix whose rank is equal to 4 [START_REF] Andreff | Image-based visual servoing of gough-stewart parallel manipulators using legs observation[END_REF]. In this expression, the derivatives ( c ui , c ḣi ) of the Plücker coordinates represent the velocity of displacement of the observed distal link i. Moreover, as shown in [START_REF] Andreff | Image-based visual servoing of gough-stewart parallel manipulators using legs observation[END_REF][START_REF] Vignolo | Comparative analysis of two types of leg-observation-based visual servoing approaches for the control of the five-bar mechanism[END_REF] the relationship between the fivebar end-effector twist c τ c and the velocities of the distal link i (i = 1, 2) is given by:

c uT i c ḣT i T = M T i c τ c (2)
where M T i is the interaction matrix relating the velocity of the distal link i to the velocity of the endeffector. This interaction matrix is of dimension (6 × 6) with a rank r ≤ 4 [START_REF] Andreff | Image-based visual servoing of gough-stewart parallel manipulators using legs observation[END_REF][START_REF] Vignolo | Comparative analysis of two types of leg-observation-based visual servoing approaches for the control of the five-bar mechanism[END_REF]. From equations ( 1) and ( 2), the interaction matrix define the relationship between c τ c and the vectors c ṅj i can be found:

H T c ṅk = M T c τ c (3) with H = [H 1 H 2 ], M = [M 1 M 2 ], c ṅk = [ c ṅ1 T 1 c ṅ1 T 2 c ṅ2 T 1 c ṅ2 T 2 ]
T . Finally, the using of the pseudo-inverse M T + = (MM T ) -1 M aims at finding the following interaction model for a single camera k:

c τ c = S T k c ṅk , with S T k = M T + H T = (MM T ) -1 MH T (4)
Then, the interaction model for the system made with four cameras (as in our application) is obtained by:

0 τ c = S T tot c ṅtot (5)
where 0 τ c is the end-effector velocity in the world frame, S T tot = [ 0 S 1 , 0 S 2 , ..., 0 S 4 ] T in which the matrix 0 S k is the interaction matrix associated with the camera k but expressed in the world frame and

c ṅtot = [ c ṅT 1 , c ṅT 2 , ..., c ṅT 4 ].

Positioning error modelling

Because a visual servoing will be used in order to ensure the positioning accuracy of the robot, the position error comes from the camera observation error and no more from the error of the encoder measurements (as it is the case in classic encoder-based control schemes [START_REF] Khalil | Modeling, Identification and Control of Robots[END_REF]). In a real-life scenario, the observation errors appear due to noise in the image. We modeled this noise by a random shift in the pixels where the image projection of the links edges meet the frame boundary (Fig. 4). In our design process, an error of 0.1 pixel was considered: it can indeed be obtained by using subpixellic observation approaches such as in [START_REF] Shi | Sub-pixel layout for superresolution with images in the octic group[END_REF] and can straightforwardly be related to a variation δn tot of the measurement vector c ṅtot . The error of observation being small, we approximated the error model relating the variation of the end-effector position δx to the variation δn tot of the measurement vector by using the following first order approximation based on the use of the interaction matrix given at (5):

δx = S T tot δn tot (6)
This error model will be used in order to characterize the robot accuracy during the design optimization process.

Control-based Design: problem formulation

After defining the error model, the optimization problem for finding the link lengths and camera locations can be formulated. The objective will be to minimize the size of the five-bar mechanism footprint in the plane (Fig. 1). Moreover, we take into consideration the performance of the controller in term of accuracy. Thus, we formulated the following optimization problem:

minimize A = LH over x subject to
Wx > Wx0 and Wy > Wy0 [START_REF] Quigley | A low-cost compliant 7-DOF robotic manipulator[END_REF] where:

x = [ 0 1 2 x c1 x c2 x c3 x c4 ] with 0 =
A1A2 , 1 = A1B1 = A2B2 and 2 = B1C = B2C are the robot link lengths, and x ci = [x c yc zc φ θ ψ] T (i = 1 to 4) is a vector defining the position and orientation of the cameras: x c , y c and z c defining the position of the camera respect to the word frame, φ, θ and ψ are the ZXZ Euler angles characterizing its orientation. A = LH is the robot footprint (Fig. 1).

Wx and Wy are the dimensions along x 0 and y 0 of the rectangular dexterous workspace in which all performance must be satisfied [START_REF] Merlet | Parallel Robots[END_REF]. In our case, the robot must have a rectangular dexterous workspace of rectangular shape with dimensions along x 0 bigger than Wx0 = 800 mm and along y 0 of Wy0 = 200 mm in which the following properties are certified:

1. no Type 2 singularities of the robot and singularities of the controller, 2. the end-effector is within the image frame for all cameras and thus the legs can be observed, 3. knowing the resolution of the cameras, the resolution of the end-effector position should be lower than 0.5 mm. The end effector resolution is computed based on the error model given by (6).

4. the static forces exerted into the passive joints are proportional to 1/ sin ξ, ξ being the angle between the distal links [START_REF] Briot | Investigation on the effort transmission in planar parallel manipulators[END_REF]. Consequently, it is decided that sin ξ should be higher than 0.1 to avoid excessive efforts in the joints.

It is necessary to mention that the error model also depends on the diameter of the cylindrical distal links which was not included in our optimization problem in order to accelerate the computation of the results. This diameter was fixed at 80 mm. The optimal design variables and the optimal position / orientation of the cameras can be seen in Tab. 2 and Tab. 3 respectively. It should be mentioned that the position/orientation of cameras 1 and 2 is obtained by symmetry with respect to the y 0 axis. Moreover, solving the problem [START_REF] Quigley | A low-cost compliant 7-DOF robotic manipulator[END_REF] allows to find the position of the largest dexterous regular workspace (Fig. 5). We can clearly observe that it is free of the singularities of the robot (and also of the controller, as they are identical for this type of visual servoing, as shown in [START_REF] Vignolo | Comparative analysis of two types of leg-observation-based visual servoing approaches for the control of the five-bar mechanism[END_REF]). The link lengths given in Tab. 2 are the inputs of the next optimization process.

The value of the positioning end-effector errors in the dexterous regular workspace is illustrated in Fig. 6. These errors are lower that the admissible ones that were fixed at 0.5 mm. Moreover as can be seen in Fig. 6, the error resolution is worsen near the singularities of the robot identical to those of the controller, as mentioned above.

Pos/Ori xc [m] yc [m] zc [m] φ [rad] θ [rad] ψ [rad] Camera3 0.1 0.3 0.5 0 π 0 Camera4 0.1 0.5 0.5 0 π 0

TOPOLOGY OPTIMIZA-TION

Topology optimization aims at finding the link shapes for a mechanical structure by optimally distributing the material placement in order to satisfy performance indices: in the literature, a classical problem is to lower the link mass under compliance requirements [START_REF] Sigmund | A 99 line topology optimization code written in MATLAB[END_REF]. In order to compute the link elastic performance, finite element method (FEM) is classically used [START_REF] Rossow | A finite element method for the optimal design of variable thickness sheets[END_REF]. Based on this physical model, the extraction of some valuable indices to be included into the optimization process must be made, as well as their gradients with respect to the decision variables [START_REF] Sigmund | Topology optimization approaches: A comparative review[END_REF]. These indices and their gradients are inputs of the optimization solvers among which we can mention the Method of Moving Asymptotes (MMA) [START_REF] Svanberg | The method of moving asymptotes -a new method for structural optimization[END_REF], the Optimality Criteria (OC) (for instance, see [START_REF] Yin | Optimality criteria method for topology optimization under multiple constraints[END_REF]), the Projected Gradient [START_REF] Chang | Direct gradient projection method with transformation of variables technique for structural topology optimization[END_REF], the Convex Linearization method (CONLIN) [START_REF] Fleury | Conlin: an efficient dual optimizer based on convex approximation concepts[END_REF]) and the Linearization Method (LM) [START_REF] Pshenichnyj | The Linearization Method for Constrained Optimization[END_REF][START_REF] Briot | Topology optimization of industrial robots: Application to a five-bar mechanism[END_REF].

In what follows, we first briefly introduce the way the finite element models are made. Then we define a performance index in order to take into account the variability into the wood mechanical behavior. Finally, we formulate the optimization problem and we solve it.

Modeling of the linkages elastic behavior

In what follows, we consider the finite element modelling of the five-bar mechanism, i.e. a robot made of four bodies (Fig. 1). The body i (i = 1, . . . , 4) is meshed with m i elements. The element j of the body i is denoted as the element ij (Fig. 7). To this element, we associate a variable ρ ij = 1 representing the density of the portion of material linked to the element. These density variables are the decision variables of the optimization problem: densities ρ ij equal to 1 will represent the presence of materials while null densities will represent the absence of materials. A challenge in topology optimization is to avoid having a final link design without too many intermediate values (0 < ρ ij < 1) that are difficult to manage by the designer.

x 0 y 0 z 0 O 0 x i y i z i O i M ij ρ ij Body i
In order to solve this problem, FEM is usually combined with an interpolation scheme which is used in order to define an artificial material. This method is called the Solid Isotropic Material with Penalization (SIMP, [START_REF] Bendsoe | Material interpolation schemes in topology optimization[END_REF]). This approach has proved to be efficient and it is probably the most widely used material interpolation scheme in topology optimization. This interpolation scheme is adopted in our work.

The SIMP scheme used to parameterize the Young's modulus value E ij of the element ij leads to the following expression:

E ij = E min + ρ p ij (E 0 -E min ), with ρ ij ∈ [0, 1] (8)
where p is the penalization factor (classically p = 3), E 0 is the real elasticity modulus of the material and E min is a very small value. E min is assigned to regions without materials in order to avoid singularity of the stiffness matrix.

Then, based on this definition of the elasticity modulus for the element ij, it is possible to build its stiffness matrix. Once all elementary matrices are de-fined, the computation of the body and robot stiffness matrices is the same as in the traditional methodology [START_REF] Shabana | Dynamics of Multibody Systems[END_REF].

Then, based on these stiffness models, the elastic performance of the mechanism can be defined: in what follows, we are going to consider the deformations u e of the five-bar mechanism at its end-effector which are linked to the wrench f e applied on the endeffector by:

f e = K e (π E , q)u e (π E , q) ⇒ u e (π E , q) = K -1
e (π E , q)f e (9) in which K e (π E , q) is the reduced stiffness matrix characterizing the overall stiffness of the assembled robot (when variable loads are applied at the endeffector location only) [START_REF] Briot | Topology optimization of industrial robots: Application to a five-bar mechanism[END_REF][START_REF] Craig | Coupling of substructures for dynamic analysis[END_REF] which depends of:

• the robot configuration q,

• the link Young's moduli (which can be different for any wooden link) which are stacked into the vector π E . For the five-bar mechanism, this vector has four components.

It should be mentioned that:

• due to page limitation constraints, the assembly of the body and robot stiffness matrices is not defined here, but expressions leading to the full robot model can be found in [START_REF] Briot | Topology optimization of industrial robots: Application to a five-bar mechanism[END_REF],

• in order to considerably decrease the computational time of the topology optimization algorithm that can be very long for the design of a robot, it is preferable to take advantage of model reduction techniques [START_REF] Craig | Coupling of substructures for dynamic analysis[END_REF], such as was done in [START_REF] Briot | Topology optimization of industrial robots: Application to a five-bar mechanism[END_REF].

Modelling of the variability in the robot deformations due to the use of wood material

Choosing a proper wood for the design of our industrial wooden robot is of the utmost importance. It must have a very good stiffness-to-mass ratio while having the lowest dimensional variability due to humidity. This is the reason why we decided to design our links with a novel type of wood named acetylated beech [START_REF] Hill | Wood Modification: Chemical, Thermal and Other Processes[END_REF]. This wood is known to have many advantages over raw wood species in terms of durability and dimensional stability. Moreover, the acytelation is a environmental-friendly chemical treatment [START_REF] Hill | Wood Modification: Chemical, Thermal and Other Processes[END_REF] which is extremely important with respect to our eco-design purpose.

Wood is an orthotropic material, i.e. it has not the same mechanical behavior depending if we apply loadings along its longitudinal, radial, or tangential direction, the longitudinal direction being along the wood fiber direction. According to the theory of elasticity, the wood is thus parameterized by three Young's moduli (E L , E R and E T ), six Poisson's ratios (ν LR , ν LT , ν RT , ν RL , ν T L , ν T R ) and three shear moduli (G LR , G LT , G RT ). However, the topology optimization is performed on an isotropic model for two reasons:

• the links of parallel robot presents a slenderness corresponding in the longitudinal direction of the wood and these links mostly resist to bending moment, therefor the E L modulus is used. Also, this explain why the value for Poisson's ratio is the mean value s the mean value of the two transversal direction (R and T ) ν = (ν LR +ν LL ) 2 .

• even if the wood is orthotropic, we have the possibility during the manufacture phase to properly design links so that they have a behavior close to the isotropy. This adequate design will be explained in Section 5.

Acetylated beech, as any other type of wood, presents a normal distribution of its longitudinal Young's modulus. Let us define as E(.) the expectation operator and σ(.) the standard deviation operator: the expectation of the Young's modulus for the acetylated beech is equal to E(E) = E = 12772 MPa while its standard deviation is equal to σ(E) = 2043.5 MPa [START_REF] Kaci | Elastostatic modelling of a wooden parallel robot[END_REF]. The Poisson's ratio ν is considered to be constant and equal to 0.3 [START_REF] Kretschmann | Ch. 5: Mechanical properties of Wood[END_REF].

Thanks to this knowledge, it is thus possible to compute the expectation vector E(u e ) and the covariance matrix Cov(u e ) associated with the endeffector deformations u e , by using the expressions given in [START_REF] Asadpoure | Robust topology optimization of structures with uncertainties in stiffness -Application to truss structures[END_REF]:

E(u e ) = u e0 + 1 2 i j u eij σ ij (10) Cov(u e ) = i j u ei u T ej σ ij + 1 2 i j k u ei u T ejk + u eij u T ek σ ijk + 1 4 i j k l u eij u T ekl (σ ijkl -σ ij σ kl ) (11)
where

• u e0 is the deformation of the end-effector given by the expression (9) evaluated when the vector of Young's moduli π E for the five-bar mechanism is equal to

π E0 = [E E E E] T , • u ei = ∂ue ∂Ei π E0 , u eij = ∂ 2 ue ∂Ei∂Ej π E0
, in which E i (E j , resp.) is the Young's modulus of the ith (jth, resp.) robot link. Detailed expressions of u ei and u eij can be found in [START_REF] Asadpoure | Robust topology optimization of structures with uncertainties in stiffness -Application to truss structures[END_REF].

• σ ij = E (E i -E)(E j -E) , • σ ijk = E (E i -E)(E j -E)(E k -E) , • σ ijkl = E (E i -E)(E j -E)(E k -E)(E l -E)
It should be mentioned that expressions [START_REF] Chaumette | Handbook of robotics[END_REF] and [START_REF] Kaci | Control-based design of a five-bar mechanism[END_REF] are approximated functions obtained thanks to a Taylor series expansion at the order two of the expression of the robot end-effector deformation u e with respect to the components of the vector π E [START_REF] Asadpoure | Robust topology optimization of structures with uncertainties in stiffness -Application to truss structures[END_REF]. We tested their validity and our results shown that, for the mechanical parameters of the acetylated beech, the error of approximation is lower than 10 %. Based on these expressions, it is thus possible to define the following constraints representing the possibility of the robot to have reliable stiffness performance:

u d (f e , q, k) ≤ u max (12) u r (f e , q, k) ≤ θ max (13)
where

u d (f e , q, k) = E( u de ) + kσ( u de ) (14)
Figure 8: Design of the proximal links: initial design domain

u r (f e , q, k) = E( u re ) + kσ( u re ) (15)
in which

• u de are the components of u e characterizing the end-effector translational deformations due to the loading f e , while u re are the components of u e characterizing the end-effector rotational deformations under the same loading,

• the functions σ( u de ) and σ( u re ) representing the standard deviations of the end-effector translational and rotational deformations, respectively, can be obtained by the following expressions:

σ( u de ) = i Cov ii (u e ), σ( u re ) = j Cov jj (u e ) (16
) in which Cov ii (u e ) (Cov jj (u e ), resp.) represents all diagonal terms of the covariance matrix Cov(u e ) corresponding to the translational (rotational, resp.) components of u e ,

• k is a positive real number.

The constraints [START_REF] Briot | The hidden robot: an efficient concept contributing to the analysis of the controllability of parallel robots in advanced visual servoing techniques[END_REF] and ( 13) are derivated from the well-known Bienaymé-Tchebichev theorem, which states that the probability, for a stochastic variable X, to have | X -E(X) | ≥ k, is lower than σ(X) 2 k 2 . Thus, by increasing the value of k in ( 12) and ( 13), we reduce the probability to have a robot whose deformations can overpass the requested thresholds u max or θ max . For instance, with k = 3, the probability to have a robot whose deformations are acceptable is greater than 88 %.

Optimization Problem

In order to speed up the optimization process, we imposed that the shapes of both proximal links are identical. Moreover, let us recall that, in order to have a good observation the distal links, these links must have a cylindrical shape and their diameter was fixed at 80 mm. This is indeed not a problem, because these links are only solicited in tension/compression when efforts are applied in the plane of the robot motion and have very few deformation. So shape of their cross section has very few impact on the overall robot deformations.

Based on the previous developments and the specifications of Tab. 1, we formulated the following topology optimization problem:

min ρ zz 1R = zz 1 + m 2 2 1 under g 1i = (u 2 d (f i , q * , k) -u 2 max )/u 2 max ≤ 0, i = 1, 2 g 2i = (u 2 d (f i , q * , k) -θ 2 max )/θ 2 max ≤ 0, i = 1, 2 (17)
where

• the objective function zz 1R = zz 1 + m 2 2
1 is a grouping of inertial parameters (zz 1 is the moment of inertia at point A 1 of the proximal link while m 2 is the mass of the distal link). This term is usually preponderant in the dynamic model of the five-bar mechanism [START_REF] Pagis | Enlarging parallel robot workspace through type-2 singularity crossing[END_REF]. Thus minimizing it has a positive impact on the dynamic properties of the robot (decrease of the maximal torque, reduction of the energy consumed, etc).

• ρ is the vector containing all decision variables ρ ij .

• q * are selected exciting configurations for the computation of the translational/rotational deformations when the loading f i is applied (i=1,2 -these loadings are defined in Tab. 1). Configurations q * can be found using an adequate methodology proposed in [START_REF] Briot | Topology optimization of industrial robots: Application to a five-bar mechanism[END_REF] and are the configurations where the deformations are the worst in the workspace. Thus, thanks to this approach, we can guarantee that the constraints are respected wherever in the workspace.

We solved this problem by using the optimization algorithm called the Linearization Method (LM) [START_REF] Pshenichnyj | The Linearization Method for Constrained Optimization[END_REF] which has proven, for our class of problems, to give results in a computational time equivalent to standard topology optimization algorithms [START_REF] Briot | Topology optimization of industrial robots: Application to a five-bar mechanism[END_REF].

Result of the topology optimization

The initial design domain for the proximal links is represented in Figure 8. Each link has two empty holes of diameters 6 cm at the extremity for the joints insertion (motor shaft at points A i , and passive rotary joints for all others. For the meshing of the links, 22152 QUA4 planar elements of size 1 × 1 mm are used and with 50 mm of thickness. All models and optimization algorithms have been encoded with Matlab in the Windows 7 environment. For a given value of k, computation time took around 5000 sec (for a Pentium 4 2.70 GHz, 16 GB of RAM). The value of the objective function zz 1R as a function of k is given in Fig. 9. The changes in the link design, when the value of k increases, are essentially an increase of the material thickness near the hole at the "right-hand-side" of the link (connection with the distal link) and the presence of additional material near the "left-hand-side" of the link (connection with the actuator). The total inertia between tests with k = 2 and k = 5 is almost doubled, thus increasing the robustness of the link with respect to the wood variability, but also decreasing the robot dynamic performance.

We decided to finally achieve a prototype based on the results obtained for k = 3. Key design features for this prototype are discussed in the next Section.

CAD DESIGN AND PRO-TOTYPING

Based on the previous optimization results, a prototype of a wooden five-bar mechanism was designed (Fig. 11). In this section, some key technological solution are highlighted: Wooden distal links: are cylindrical for reason of facility of the observation with camera. They are partially hollowed out to limit their mass.

Wooden proximal links: are based on topology optimization results. They are build from 7 laser cut layers (7 mm thick), glued together with polyurethane glue. In order to improve the elastic behavior of these links along all directions, the layers are assembled by differently orienting the directions of the fibers for each layer, so that the links behave as they were made of a material with properties as close as possible from the isotropy. From top to bot- tom, angles between longitudinal axis of the link and wood grain ares 0 • ;45 • ;-45 • ;0 • ;-45 • ;45 • ,0 • . CAD design of the proximal links is shown in Fig. 10(a), it should be mentioned that the difference between the CAD model and the topology optimization results is due to local edge distances needed around the holes to ensure wood resistance, but the truss-like shape is kept.

Wooden chassis: is made with acetylated pine, and supports two direct drive actuators (SIMOTTICS S-1FL6 servomotors) as shown in Fig. 10(b). Indeed, our robot should be able to fully rotate the proximal and distal links without interferences. Thus, we chose to set the actuators under the robot links.

Revolute joints: high precision angular contact ball bearings are used to design the revolute joints. This configuration delivers high level of stability and minimal backlash for the assembly.

FEM simulations of the CAD model were achieved in order to verify the overall deformation of the robot. Desired deformations under the loadings given in Tab. 1 were verified. We also tested the robot deformation under a loading of 10 N along z 0 : results about 400 microns, which is acceptable.

At this step of the project RobEcolo, the wooden five-bar mechanism is built and future works will con-

CONCLUSION

In order to try to decrease the environmental impact of robot during their manufacturing phase, we investigated in this paper the possibility to design a wooden industrial robot whose performance in terms of accuracy and stiffness are guaranteed, even if exhibits significant dimensional ans mechanical variabilities. Thus, in order to be able to ensure stiffness and the accuracy of the wooden robot, two design methodologies are introduced in this paper. These two approaches are solved in cascade and aim at finding the robot optimal design architecture for which the impact of the wood properties variability is minimal. Based on the obtained architecture of robot, a prototype of a wooden five-bar mechanism is designed and built.

After the prototype validation phase, a first perspective of our work is to implement the sensor-based controllers. Then, several experiments for characterizing the robot stiffness and accuracy will be made to validate all the theoretical results obtained in the framework of the project RobEcolo.

Acknowledgment

This work was supported by the project RobEcolo funded by the French Région Pays de la Loire (Convention No. 2015-10773).

Figure 1 :

 1 Figure 1: Five-bar mechanism

Figure 2 :Figure 3 :

 23 Figure 2: A general sensor-based control scheme

Figure 4 :

 4 Figure 4: Result of a one-pixel error on the intersection of the image boundary and the observed line

Figure 5 :

 5 Figure 5: Largest regular dexterous workspace of the optimized five-bar mechanism

Figure 6 :

 6 Figure 6: Resolution of the end-effector pose in the LRDW

Figure 7 :

 7 Figure 7: A body in space meshed with finite elements.

Figure 9 :

 9 Figure 9: Design of the proximal links: final results for different values of k

  (a) Proximal Link: CAD Design (b) Final Prototype

Figure 10 :

 10 Figure 10: CAD Design of a wooden five-bar mechanism.

Figure 11 :

 11 Figure 11: Final prototype of a wooden five-bar mechanism.

Table 3 :

 3 Optimal Position and Orientation of Cameras

The authors would like to thank Pauline Lafoux and Victoria Safyannikova for their work on the CAD design of the RobEcolo prototype.