Sébastien Briot
email: sebastien.briot@irccyn.ec-nantes.fr

Alexandre Goldsztejn
email: alexandre.goldsztejn@irccyn.ec-nantes.fr

Topology Optimization of Industrial Robots: Application to a Five-bar Mechanism

Keywords: Robot design, deformation, natural frequency, topological optimization, linearization method

Recent works introduced topology optimization in the design of robots, but the proposed methodologies led to a local optimization of the performance. Moreover, most of performance indices used are not in strong relation with easy-to-understand technological requirements.

We propose a methodology that is able to perform a topology optimization for robots, valid globally in the workspace or for a set of given trajectories, and which is based on the use of technology-oriented performance criteria. In order to enforce the chosen performance indices to be valid globally, optimal robot configurations or trajectories for which extreme performance will be attained are computed, and iteratively updated.

In order to decrease the computational time associated with these performance indices, we exploit the structure of the elastic models in order to reduce their computational complexity.

Finally, we use an optimization algorithm called the Linearization Method which gives results in a computational time equivalent to standard topology optimization algorithms, but its implementation is less complex and makes it quite easy to perform modification or improvement.

The methodology is applied for the design of a five-bar mechanism. We show that our approach leaded to a robust optimization of the robot performance over the whole workspace.

Introduction

Robots are expected to perform a large variety of tasks. However, it is not wise to believe that a single robot will be able to achieve all conceivable tasks. Inherent robot limitations arise from its own physical performance (accuracy, deformation, vibrations, etc.), which are a combination of the performance of the mechanical architecture and of the controller.

Good performance of the mechanical architecture can be obtained via optimal design [START_REF] Merlet | Optimal design of robots[END_REF]. The usual design methodology proposed by French [START_REF] French | Conceptual Design for Engineers[END_REF] is illustrated in Fig. 1. The first step is to analyze the need in order to formulate the design problem. The second phase focuses on the preliminary design and aims to synthesize design concepts (for instance, new types of robot architectures [START_REF] Gogu | Structural Synthesis of Parallel Robots[END_REF][START_REF] Kong | Type Synthesis of Parallel Mechanisms[END_REF]) and to select the best design alternatives with respect to given criteria (e.g. complexity [START_REF] Caro | The rule-based conceptual design of the architecture of serial Schönflies-motion generators[END_REF], singularity [START_REF] Ben-Horin | Singularity analysis of a class of parallel robots based on grassmanncayley algebra[END_REF]). The third phase, denoted as the advanced design phase or embodiment of schemes, deals with the computation of the dimensions and shapes of the product element in order to fulfill performance criteria in terms of:

• geometric performance: e.g. workspace size and shape under joint limitations, link collisions [START_REF] Merlet | Designing a parallel manipulator for a specific workspace[END_REF][START_REF] Merlet | Workspaces of planar parallel manipulators[END_REF], accuracy under input errors [START_REF] Briot | Accuracy analysis of 3-dof planar parallel robots[END_REF], under clearance [START_REF] Binaud | The kinematic sensitivity of robotic manipulators to joint clearances[END_REF][START_REF] Venanzi | A new technique for clearance influence analysis in spatial mechanisms[END_REF] or link manufacturing errors [START_REF] Caro | Sensitivity analysis of the orthoglide: A three-dof translational parallel kinematic machine[END_REF],

• kinematics / kinetostatics: e.g. velocity transmission ratio [START_REF] Liu | Kinematics, singularity and workspace of planar 5R symmetrical parallel mechanisms[END_REF][START_REF] Merlet | Jacobian, manipulability, condition number, and accuracy of parallel robots[END_REF][START_REF] Gosselin | A global performance index for the kinematic optimization of robotic manipulators[END_REF], effort transmission [START_REF] Arakelian | Increase of singularity-free zones in the workspace of parallel manipulators using mechanisms of variable structure[END_REF][START_REF] Balli | Transmission angle in mechanisms[END_REF][START_REF] Chen | Generalized transmission index and transmission quality for spatial linkages[END_REF][START_REF] Takeda | Motion transmissibility of in-parallel actuated manipulators[END_REF][START_REF] Hubert | Static of parallel manipulators and closeness to singularity[END_REF][START_REF] Chen | Evaluation of force/torque transmission quality for parallel manipulators[END_REF],

Embodiment of schemes Detailling Selected schemes Working drawings etc.

Feedback

Figure 1: French design process [START_REF] French | Conceptual Design for Engineers[END_REF] • dynamics: e.g. moving mass reduction [START_REF] Germain | Optimal design of the IRSBot-2 based on an optimized test trajectory[END_REF], maximal input torques [START_REF] Briot | Technology-oriented optimization of the secondary design parameters of robots for high-speed machining applications[END_REF], static [START_REF] Barents | Spring-to-spring balancing as energy-free adjustment method in gravity equilibrators[END_REF] or dynamic balancing [START_REF] Chaudhary | Dynamics and balancing of multibody systems[END_REF][START_REF] V. Arakelian | Shaking force and shaking moment balancing of mechanisms: a historical review with new examples[END_REF] conditions, decoupled structure of the dynamics equations [START_REF] Asada | Decoupling of manipulator inertia tensor by mass redistribution[END_REF][START_REF] Hess-Coelho | Decoupling of dynamic equations by means of adaptive balancing of 2-dof open-loop mechanisms[END_REF],

• elasticity: static or dynamic deformations, natural frequencies and vibrations [START_REF] Briot | On the optimal design of parallel robots taking into account their deformations and natural frequencies[END_REF][START_REF] Bouzgarrou | Rigidity analysis of t3r1 parallel robot uncoupled kinematics[END_REF][START_REF] Cammarata | Elastodynamic optimization of a 3T1R parallel manipulator[END_REF].

The fourth phase is the detailed design stage and consists in obtaining the working drawings of the product elements, in synthesizing their dimensional and geometric tolerances [START_REF] Caro | Tolerance synthesis of mechanisms: a robust design approach[END_REF], and in manufacturing the prototypes.

The design optimization problem treated in the advanced design phase is usually formulated as a multicriteria optimization problem and it is most of the time solved in cascade in order to reduce its complexity [START_REF] Germain | Optimal design of the IRSBot-2 based on an optimized test trajectory[END_REF].

In a first step, the multicriteria optimization problem takes only into account geometric, kinematic and kinetostatic constraints and objectives and allows for fixing the primary geometric parameters of the robot (lengths of links, angles between the joint axes etc.) [START_REF] Briot | Optimal technology-oriented design of parallel robots for highspeed machining applications[END_REF]. In a second step, the secondary geometric parameters are found (size of the link cross-sections, link mass distribution, or more generally link shapes) taken into account dynamic, elastostatic and elastodynamic aspects [START_REF] Briot | Technology-oriented optimization of the secondary design parameters of robots for high-speed machining applications[END_REF].

The link shape optimization of robots is probably the most time-consuming step of the optimal design process. This is due to the complexity of the model involved, especially the elastic models, which must be computed thousands of times (and even more) in order to calculate the robot elastic performance in many robot configurations for a given set of design variables [START_REF] Germain | Optimal design of the IRSBot-2 based on an optimized test trajectory[END_REF]. This is necessary in order to ensure that the performance can be guaranteed in a wide range of robot configurations [START_REF] Gosselin | A global performance index for the kinematic optimization of robotic manipulators[END_REF]. As a result, in order to decrease the time of computation, a common approach is to reduce the number of design variables. It can be easily reduced by doing a parametric optimization [START_REF] Germain | Optimal design of the IRSBot-2 based on an optimized test trajectory[END_REF][START_REF] Cammarata | Elastodynamic optimization of a 3T1R parallel manipulator[END_REF], i.e. by modeling links using beam theory [START_REF] Shabana | Dynamics of Multibody Systems[END_REF] and by considering that the geometry of the beam cross-sections is fixed (for instance, circle, square, rectangle, I-shape) but parameterized by a limited number of variables (radius for circles, edge lengths for squares, rectangles, I-shapes).

This approach is known not to be the more accurate for finding the optimal design of links, contrary to topology optimization [START_REF] Sigmund | Topology optimization approaches: A comparative review[END_REF]. Topology optimization was for instance used for the design of compliant mechanisms [START_REF] Chen | The local optimum in topology optimization of compliant mechanisms[END_REF][START_REF] Venkiteswaran | Compliant mechanism design through topology optimization using pseudo-rigid-body models[END_REF][START_REF] Saxena | Topology optimization of compliant mechanisms with strength consider[END_REF]. This latter technique aims at optimizing the material distribution in a link in order to satisfy performance criteria: a classical problem met in the literature is to minimize the link mass under compliance constraints [START_REF] Sigmund | A 99 line topology optimization code written in MATLAB[END_REF]. The link shape is meshed, and deformation and vibration models are computed using Finite Element Methods (FEM) [START_REF] Rossow | A finite element method for the optimal design of variable thickness sheets[END_REF]. The presence of one element of material is parameterized by a design variable varying from 0 to 1, 1 means that there is material while 0 represents a void. As a result, in order to have a refined prediction of the link behavior and a refined visualization of the link shape, this method usually leads to a vector of design variables in the optimization process containing dozens of thousands of components. Topology optimization is thus most of the time computationally expensive due to both the complexity of the models involved and the high number of design variables. Therefore, it is few used in robot design.

However, recent works introduced this technique in the design of robots. First attempts optimized the robot topology for a single loading case. For instance, [START_REF] Albers | Upper body of a new humanoid robotthe design of ARMAR III[END_REF] optimized the torso of a humanoid robot with an objective of minimal mass under compliance constraints while [START_REF] Lohmeier | Leg design for a humanoid walking robot[END_REF][START_REF] Lohmeier | Humanoid robot LOLA[END_REF] optimized the pelvis of a walking robot. [START_REF] Kwon | Biped humanoid robot Mahru III[END_REF] employed topology optimization method to develop stiff and light frames of the lower body for stable walking of their humanoid robot. [START_REF] Yunfei | Structural Topology Optimization for a Robot Upper Arm Based on SIMP Method[END_REF][START_REF] Huang | The topology optimization for l-shape arm of Motorman-HP20 robot[END_REF] both optimized the shape of the upper-arm for a 6-degreesof-freedom (dof) industrial robot. [START_REF] Oliveri | Structural and topological optimization in robot design[END_REF] improves the shape of a chassis for a mobile robot.

Optimizing the robot link shapes for a single loading case does not take into account the intrinsic nature of a robot whose performance varies with the configuration, the loading and with the time. Therefore other works proposed alternative approaches. For instance, [START_REF] Albers | Methods for lightweight design of mechanical components in humanoid robots[END_REF][START_REF] Albers | System based topology optimization as development tools for lightweight components in humanoid robots[END_REF] considered a given set of reference control signals for the motion of the robot arm and look at some robot performance, such as overshoot, controller settling time, number of oscillations, final deviation or also actuator energy consumption. [START_REF] Kim | Topology optimization of industrial robots for system-level stiffness maximization by using part-level metamodels[END_REF] optimized the link shapes of 3-dof robots under varying configurations while addressing the problem of the reduction of the computational time by dividing the optimization problem into subproblems with lower computational complexity. They minimized the strain energy while constraining the robot mass. [START_REF] Hong | Optimization of flexible components of multibody systems via equivalent static loads[END_REF] optimized the shape of the pelvis of a humanoid robot by using equivalent static loads with the objective to minimize the strain energy.

All these works introduced the variation of the robot configuration but they constrained the robot to move on a limited set of trajectories. This is because the computation of their performance indices is extremely timeconsuming. As a result, the optimization is still local, and it is not possible to ensure that the performance will be kept along another set of trajectories. Moreover, most of authors used criteria whose physical sense is not straightforward: For instance, it is difficult to say what is a good value for the strain energy or for the compliance while it seems simpler to understand what are good performance in terms of (static or dynamic) deformations, natural frequencies, input torques or actuator energy consumption. These indices are said to be "technology oriented" as we can relate their value to some technological constraints, as it will be detailed later in the paper.

In the present paper, we propose a methodology in order to perform a topology optimization for robots, valid globally in the workspace or for a set of given trajectories, and which is based on the use of technologyoriented performance criteria. In order to enforce the chosen performance indices to be globally valid, optimal robot configurations or trajectories for which extreme performance will be attained are computed, and iteratively updated.

In order to decrease the computational time associated with the computation of performance indices in numerous configurations, we exploit the structure of the elastic models (which have the highest computational costs) in order to reduce their computational complexity. Indeed, we show that it is possible to use configuration-independent model reduction techniques in order to considerably decrease the size of the stiffness and mass matrices of each link, expressed in their local frame. Then, these reduced matrices are used in order to build the configuration-dependent robot elastic models which are shown to be computationally efficient.

Finally, we use an 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 (e.g. Method of Moving Asymptotes (MMA) [START_REF] Svanberg | The method of moving asymptotes -a new method for structural optimization[END_REF], Optimality Criteria (OC) (for instance, see [START_REF] Yin | Optimality criteria method for topology optimization under multiple constraints[END_REF]), Projected Gradient [START_REF] Chang | Direct gradient projection method with transformation of variables technique for structural topology optimization[END_REF], Convex Linearization method (CONLIN) [START_REF] Fleury | Conlin: an efficient dual optimizer based on convex approximation concepts[END_REF]) but its implementation is less complex and makes it quite easy to perform modification or improvement.

The paper is written as follows. First, in Section 2, we introduce the LM algorithm. Then, in Section 3, we discuss the several technology-oriented performance indices that could be used for our optimization problem and explain our choice. We also provide a computationally efficient method for their calculation. Section 4 introduces the procedure to choose the optimal configurations and trajectories that will be used in order to compute the robot performance. In Section 5, the topology optimization of a five-bar mechanism is performed.

Finally, in Section 6, our conclusions are drawn.

Topology optimization formulation and resolution

Problem formulation

The topology optimization problems have the standard form of a nonlinear program:

min ρ∈[0,1] n g(ρ)≤0 h(ρ)=0 f (ρ) (1)
Each component of the vector of decision variables ρ ∈ [0, 1] n ⊆ R n represents the density of a given element, defined in Section 3. The objective function and the constraints usually model some performance indices (described in the next sections) or constraints on the structure (e.g., some constant mass for easing the solving process).

The main characteristics of these optimization problems is the high number of variables, usually several thousands, and the availability of the symbolic expression of the gradient. Therefore, the aim of the solving method is basically to use as well as possible the linear models of the constraints to prevent too many function and gradient evaluations.

The most used method for solving this class of problems is MMA [START_REF] Svanberg | The method of moving asymptotes -a new method for structural optimization[END_REF]. It turns out to be more efficient than other methods, and is available under the form of a Matlab code, which is the only one that handles multiple inequality and equality constraints. However, it presents several drawbacks that lead us to consider alternative methods. First, MMA shows some unexplained surprising behaviors. In particular, it fails on badly conditioned problems, e.g., minimizing ρ 1 + ρ 2 subject to ρ 2 1 + 1000ρ 2 2 ≤ 11 . Second, MMA and its corresponding code are very complex, making very difficult any modification or improvement. For these reasons, we considered using a more standard optimization method, namely the Linearization Method.

The Linearization Method

The projected gradient method is the most used method for large scale convex optimization. It consists in choosing the search direction by projecting the cost function opposite gradient on the constraints. This projection is well defined only for convex constraints. A straightforward extension to non-convex problems consists in projecting the cost function opposite gradient on the constraint linearization instead of the constraint, and is called the Linearization Method. Up to our knowledge, this method seems to be first introduced in [START_REF] Pshenichnyi | Algorithms for general mathematical programming problems[END_REF] (see [START_REF] Pshenichnyj | The Linearization Method for Constrained Optimization[END_REF][START_REF] Poljak | Introduction to optimization[END_REF] for more comprehensive studies, including convergence properties of the method), but has not been applied to topology optimization so far. The basic idea of the method is a standard iteration

ρ k+1 = ρ k + α k d k ,
where α k is a step size and d k is chosen as illustrated on Figure 2. The feasible set of the constraint g(ρ) = ρ 2 1 + ρ 2 2 -1 ≤ 0 is the disk, and the linearization of the constraint at ρ k is the half-plane to the minimizer of the problem for sufficiently small step sizes. Formally, the Linearization Method is defined as follows:

g(ρ k ) + ∇g(ρ k )(ρ -ρ k ) ≤ 0. The vector u k = ρ k -s k ∇ f (ρ k ) is
v k = proj L k ρ k -s k ∇ f (ρ k ) (2) ρ k+1 = ρ k + α k (v k -ρ k ). (3)
with

L k = {ρ ∈ R n : g(ρ k ) + ∇g(ρ k ) T (ρ -ρ k ) ≤ 0, h(ρ k ) + ∇h(ρ k ) T (ρ -ρ k ) = 0}, (4)
where the inequality stands component-wise. The projection operation proj L k (u) projects the vector u onto the set L k , which is here defined as linear inequality and equality constraints. It is therefore convex, and the projection is well defined.

Implementation details. The projection onto the convex set L k , although well defined, is difficult in practice due to the high number of variables. We used an approximation computed by successive projections [START_REF] Combettes | Method of successive projections for finding a common point of sets in metric spaces[END_REF].

The high number of variables although prevents any line search, which requires too many function evaluations. Therefore, we implemented a constant step size strategy α k = 1 and s k = s 0 [START_REF] Bertsekas | Nonlinear Programming[END_REF]. However, several studies have shown that enforcing a maximal norm for the step increases the overall efficiency. Therefore, we implemented the following heuristic step size rule:

α k = 1 if v k -ρ k ∞ ≤ ᾱ ᾱ v k -ρ k ∞ otherwise. ( 5
)
Figure 2: Left: Projected gradient. Right: Linearization method (vectors are scaled for the sake of lisibility).

The constants s 0 and ᾱ need to be fixed for each problem with the usual drawbacks of constant step size methods. Too large initial steps cause divergence, while too small initial steps slow down the convergence.

However, this behavior is easily identified on the first iterations of the algorithm, and a correct step size is easy to find. Another heuristic dedicated to topology optimization has been implemented. Preliminary experiments showed that the Linearization Method did not remove enough material during the first iterations. We fixed this issue by considering a descent direction with a larger infinite norm, having the effect of removing more material. Equation ( 2) is replaced by

v k = proj L k ρ k -s k d k with d ki =          ∇ f (ρ k ) ∞ if ∇ f (ρ k ) i > 0 -∇ f (ρ k ) ∞ if ∇ f (ρ k ) i < 0 0 otherwise (6)
This is indeed a descent direction since ∇ f (ρ k ) T (-d k ) < 0. Finally, a simplified heavy-ball process [START_REF] Poljak | Introduction to optimization[END_REF] is added to the method to speedup its convergence. An implementation of this Linearization Method is available on the website of the second author.

Performance indices

Selection of the performance indices

In our work, we focus on the selection on technology-oriented performance indices whose computational cost can be acceptable for an optimization process.

The selection of the performance indices for the optimization of a mechanical system is not a simple task and is highly dependent on the sensitivity of the man-or-the-art. For instance, as mentioned in the introduction, past works [START_REF] Kim | Topology optimization of industrial robots for system-level stiffness maximization by using part-level metamodels[END_REF][START_REF] Hong | Optimization of flexible components of multibody systems via equivalent static loads[END_REF] used as criteria the strain energy or the compliance of the robot, but the physical sense of these indices is not straightforward. Indeed, engineers would wonder what is an acceptable value for the strain energy or for the compliance of the robot. Other works [START_REF] Albers | Methods for lightweight design of mechanical components in humanoid robots[END_REF][START_REF] Albers | System based topology optimization as development tools for lightweight components in humanoid robots[END_REF] considered performance such as overshoot, controller settling time, number of oscillations, final deviation or also actuator energy consumption. All these indices are closer from technological constraints, however the authors constrain the robot to move on a limited set of trajectories in order to decrease the computational cost of the problem. Therefore, the optimization of the robot stays local.

In order to solve these issues, we prefer to use for the optimization process the following indices (the list is not necessary exhaustive), which can be strongly related to the technological constraints:

• robot inertial parameters: the authors of [START_REF] Kim | Topology optimization of industrial robots for system-level stiffness maximization by using part-level metamodels[END_REF][START_REF] Hong | Optimization of flexible components of multibody systems via equivalent static loads[END_REF] used the robot mass as a criterion performance, which is one of the most obvious performance index, as the mass plays an important role on the value of the input efforts, but also can be linked to the robot cost (lower mass leads to the use of less material for the robot design). However, the robot mass is not necessarily the only inertial parameter to be taken into account. For instance, as shown in [START_REF] Pagis | Enlarging parallel robot workspace through type-2 singularity crossing[END_REF], the most influential term in the dynamic model of the five-bar mechanism (Fig. 3) is a grouped inertial parameter equal to zz R = zz i1 + 2 i1 m i2 , where zz i1 is the moment of inertia around z 0 of the link between the points A i and B i computed at A i , A i B i is the distance between points A i and B i and m i2 is the mass of the link between the points B i and C. Therefore, minimizing zz R is more likely to minimize the robot input torques than minimizing the robot mass.

B 1 C (x, y) x 0 y 0 O q 1 θ 1 q 2 θ 2 A 1 B 2 A 2
• robot input efforts or energy consumption: the cost of an actuator is related to its power or to its continuous or peak torques. We know that, for a usual robot, more than 90 % of the dynamic effects come from the "rigid behavior" of the mechanical architecture (remaining 10 % being due to elasticity, impacts due to clearance in joints, and unmodelled friction effects) [START_REF] Gautier | Dynamic identification of robots with power model[END_REF]. Therefore, criteria based on the input efforts, energy or power consumption can be computed with a very good accuracy neglecting the elastic behavior of the links.

• robot static deformations: static deformations under external (potentially variable) loading are classical performance indices when designing a mechanical system.

• robot natural frequencies: the first natural frequencies are associated with the highest level of energy due to vibrations, and the first modes represent the displacements with the highest amplitude. Moreover, the first natural frequencies are also used in order to set the bandwidth of controllers: it is usually said that the cutoff frequency for a robot controller must be set around the half of the first natural frequency of the mechanical architecture [START_REF] Khalil | Modeling, Identification and Control of Robots[END_REF]. If not, the controller can become unstable due to unmodeled dynamics.

We believe that natural frequencies are more important indices than the amplitude of the robot dynamic deformations used [START_REF] Albers | Methods for lightweight design of mechanical components in humanoid robots[END_REF][START_REF] Albers | System based topology optimization as development tools for lightweight components in humanoid robots[END_REF][START_REF] Kim | Topology optimization of industrial robots for system-level stiffness maximization by using part-level metamodels[END_REF][START_REF] Hong | Optimization of flexible components of multibody systems via equivalent static loads[END_REF]: The lack of accuracy due to oscillation can be compensated thanks to advanced controller (e.g. [START_REF] Singhose | Design and implementation of time-optimal negative input shapers[END_REF][START_REF] Douat | Hinfiny control applied to the vibration minimization of the parallel robot Par2[END_REF]) while low robot natural frequencies cannot be increased through control.

The computation of the robot inertial parameters and input efforts presents no difficulty, and this is the reason why they are given in Appendix A and Appendix B. We show in these appendices that the expression of all rigid inertial parameters and of the rigid dynamic model is linear with respect to the decision variables.

Below, we prefer to stress the issue of the calculation of the robot deformations and natural frequencies, whose cost of computation is huge, and which must be computed for numerous robot configurations. As mentioned in introduction, in order to considerably decrease the computation cost of these performance indices, we exploit the structure of the elastic models (which have the highest computational costs) in order to reduce their computational complexity. Indeed, we show below that it is possible to use configurationindependent model reduction techniques in order to considerably decrease the size of the stiffness and mass matrices of each link, expressed in their local frame. Then, these reduced matrices are used in order to build the configuration-dependent robot elastic models which are shown to be computationally efficient.

x 0 y 0 z 0 O 0 x i y i z i O i M ij m ij

Body i

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

Static deformations

In what follows, we consider a robot made of n bodies. The body i is meshed with m elements. The element j of the body i is denoted as the element i j (Fig. 4).

Material interpolation scheme and stiffness matrix of a link

Topology optimization uses the same physical model as in the FEM for modelling of a link, except that we use an interpolation scheme 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]) and is known to be the most effective and the most widely used material interpolation scheme. This material interpolation scheme is adopted in order to avoid having a result of optimization with too much intermediate material density, i.e. by denoting as ρ i j the density associated with the element i j, in order to have a black (ρ i j = 1) and white (ρ i j = 0) solution without too many grey elements (0 < ρ i j < 1).

The SIMP scheme is defined as follows:

E i j = E min + ρ p i j (E 0 -E min ), with ρ i j ∈ [0, 1] (7)
where E 0 is the Young's modulus of the material, E min is a very small stiffness value assigned to void regions in order to prevent the stiffness matrix from becoming singular, p (typically p = 3) is the penalization factor, and E i j is the Young's modulus of element j of the body i corresponding to the density variable ρ i j .

Then the stiffness matrix associated with the element i j is thus given by:

K i j = E i j K (0) i j = E min + ρ p i j (E 0 -E min ) K (0) i j (8)
where K (0) i j is the stiffness matrix of a single element computed for a Young's modulus equal to 1, and the potential elastic energy of the element is:

U e i j = 1 2 u T i j K i j u i j (9)
in which u i j is the vector of the element i j nodal displacements. The total elastic energy of the body i is thus equal to

U e i = m j=1 U e i j = 1 2 m j=1 u T i j K i j u i j = 1 2 u T i tot K i tot u i tot (10)
where

• u i tot = [u T i1 . . . u T im ]
T is the vector stacking all nodal displacements for all m elements of the body i

• K i tot is a block-diagonal matrix stacking on its diagonal all elementary stiffness matrices as follows:

K i tot =             K i1 0 . . . 0 K im             (11)
Finally, the link stiffness matrix K i of the body i can be obtained by taking into account the fact that the nodal displacements of the element i j are equal to the nodal displacements of its adjacent elements. As a result, the expression of the vector u i tot can be obtained from a reduced set of independent coordinates u i [START_REF] Imbert | Analyse des structures par lments finis[END_REF] as follows:

u i tot = J i u i (12)
where, for a single body, J i is constant. Introducing ( 12) into (10), we get

U e i = 1 2 u T i K i u i ( 13
)
where K i = J T i K i tot J i is the body stiffness matrix which relates the nodal displacements u i to the forces f i exerted on the nodes by the relation [START_REF] Shabana | Dynamics of Multibody Systems[END_REF]:

f i = ∂U e i ∂u i = K i u i (14)
In the following the paper, we consider that the components of the vector u i are ordered and split into two parts such as

u i = u i f u i l (15)
where u i f is the vector of displacements of the nodes on which no force or displacements are imposed, while u i l is the vector of displacements for the nodes which are constrained with forces and/or displacements. The nodes whose displacements are represented by the vector u i l are named the interface nodes.

Efficient computation of the elastostatic model

The usual computation of the robot elastostatic model is shown in Appendix C. If this approach is used for the computation of the robot stiffness matrix K, its final dimension is huge (typically, in our problems, the dimension of K is greater than (10 5 × 10 5 )), thus leading to a huge computational cost for solving the deformation problem given at equation (C.5), problem which must be solved for any tested robot configurations.

It is possible to considerably reduce the computational cost for solving the elastostatic model by using a model reduction technique for the robot links as follows.

In order to perform the reduction of the model size, we use the partitioning of the vector u i shown in [START_REF] Gosselin | A global performance index for the kinematic optimization of robotic manipulators[END_REF].

We denote as a the size of the vector u i f and as b the size of the vector u i l . From the definitions of u i f and u i l , we also partition the vector f i in ( 14) such that

f i = 0 (a×1) f i l ( 16
)
in which 0 (a×1) is a zero vector of size a which represents the absence of interaction between the environment and the nodes whose displacements are parameterized by u i f . Taking into account the separation of components in the vectors u i and f i , the equation ( 14) can be rewritten as:

0 (a×1) f i l =       K (11) i K (12) i K (21) i K (22) i       u i f u i l ( 17
)
where K (11) i is a (a × a) matrix, K (12) i = K (21) T i is a (a × b) matrix and K (22) i is a (b × b) matrix.

(17) can be split into the two following equations:

0 (a×1) = K (11) i u i f + K (12) i u i l (18)
f i l = K (21) i u i f + K (22) i u i l (19)
By using [START_REF] Chen | Generalized transmission index and transmission quality for spatial linkages[END_REF], it is possible to express u i f as a function of the displacements of the interface nodes u i l as:

u i f = Φ s i u i l , in which Φ s i = -K (11) i -1 K (12) i (20)
Φ s i is the matrix of the static modes [START_REF] Craig | Coupling of substructures for dynamic analysis[END_REF].

Introducing [START_REF] Hubert | Static of parallel manipulators and closeness to singularity[END_REF] into [START_REF] Takeda | Motion transmissibility of in-parallel actuated manipulators[END_REF], we get

f i l = K red i u i l , in which K red i = K (21) i Φ s i + K (22) i ( 21
)
while the potential elastic energy ( 13) can be rewritten as

U e i = 1 2 u T i l K red i u i l (22)
K red i
is the reduced stiffness matrix associated with the displacements of the interface nodes of the body i.

This matrix is usually of small dimension (typically, in our 2D problems, of dimension (6 × 6), while in 3D it is of dimension (12 × 12)) and it does not depend on the robot configuration, which is a great advantage because it must be computed only once at each step of the optimization algorithm, whatever the number of tested configurations for the robot. Indeed, most of the computational cost of its computation is due to the computation of the matrix Φ s i in [START_REF] Hubert | Static of parallel manipulators and closeness to singularity[END_REF]. However, for two bodies i and k, the computation of the matrices Φ s i and Φ s k is independent. Thus the computation can be made in parallel on multi-core computers in order to save computational time.

Considering now the robot made of n bodies, the full potential elastic of the system is given by:

U e = n i=1 U e i = 1 2 n i=1 u T i l K red i u i l = 1 2 u red T tot K red tot u red tot ( 23
)
where

• u red tot = [u T 1 l . . . u T n l ]
T is the vector stacking all interface node displacements for all n robot bodies

• K l is a block-diagonal matrix stacking on its diagonal all bodies stiffness matrices as follows:

K red tot =             K red 1 0 . . . 0 K red n             (24)
The reduced robot stiffness matrix K r can be obtained by taking into account the fact the robot bodies are connected altogether through the interface nodes. As a result, the expression of the vector u red tot can be obtained from a reduced set of independent coordinates u r [START_REF] Shabana | Dynamics of Multibody Systems[END_REF][START_REF] Imbert | Analyse des structures par lments finis[END_REF] as follows:

u red tot = J r (q)u r (25)
where J r (q) depends on the robot configuration q but not on the nodal displacements u r (assumption of small perturbations).

Introducing [START_REF] Chaudhary | Dynamics and balancing of multibody systems[END_REF] into [START_REF] Briot | Technology-oriented optimization of the secondary design parameters of robots for high-speed machining applications[END_REF], we get

U e = 1 2 u T r K r u r (26)
where K r = J T r (q)K red tot J r (q) is the robot stiffness matrix which relates the nodal displacements u r to external forces f r exerted on the considered nodes by the relation:

f r = ∂U e ∂u r = K r u r (27)
For a 2D problem, the typical size of the matrix K r is lower than (6n × 6n), n being the number of robot bodies ((12n × 12n) for 3D problems). Therefore, even if the problems (C.5) and ( 27) will give exactly the same results of computation, the resolution of the equation ( 27) is much more efficient than the resolution of the equation (C.5) due to the considerably reduced size of the problem.

Stiffness performance index

Usually, it is not necessary to constrain all deformations u r , but a set of them. Let us define a row vector e able to extract c components u c of u r such as u c = e u r . It is then quite usual to constrain the two-norm of this vector u c or its square (in order to avoid the use of square roots in the definition of the performance index)

given by:

c 1 = u T c u c = u T r e T e u r (28)
The derivative of this criterion with respect to the decision variable ρ i j is given in Appendix E.

Natural frequencies

Mass matrix of a link

The mass matrix associated with the element i j is given by:

M i j = ρ i j M (0) i j (29)
where M (0) i j is the mass matrix of a single element computed for a density equal to 1. Accordingly, the kinetic energy of the element due to elastic oscillations is:

T e i j = 1 2 uT i j M i j ui j ( 30
)
in which ui j is the vector of the element i j nodal velocities. The total kinetic energy of the body i is thus equal to

T e i = m j=1 T e i j = 1 2 m j=1 uT i j M i j ui j = 1 2 uT i tot M i tot ui tot (31)
where M i tot is a block-diagonal matrix stacking on its diagonal all elementary mass matrices as follows:

M i tot =             M i1 0 . . . 0 M im             (32)
Finally, differentiating [START_REF] Caro | Sensitivity analysis of the orthoglide: A three-dof translational parallel kinematic machine[END_REF] with respect to time and recalling that the matrix J i is constant, we obtain:

ui tot = J i ui (33)
Introducing ( 33) into (31), we get

T e i = 1 2 uT i M i ui ( 34
)
where M i = J T i M i tot J i is the body mass matrix.

Efficient computation of the elastodynamic model

The usual computation of the robot elastodynamic model is shown in Appendix D. If this approach is used for the computation of the natural frequencies, as previously met for the computation of the elastostatic model, the cost of their evaluation will be prohibitive due to the huge dimension of the matrices to evaluate for any tested robot configurations.

It is possible to considerably reduce the computational cost for computing the robot natural frequencies by using a Craig-Bampton model reduction technique [START_REF] Craig | Coupling of substructures for dynamic analysis[END_REF] applied to each robot body as follows.

The Craig-Bampton model reduction technique for a body is based on the assumption that the body nodal coordinates u i f (recall that u i f is the vector of displacements of the nodes on which no force or displacements are imposed) can be expressed as a function of the nodal coordinates u i l corresponding to the interface nodes plus a term characterizing their vibratory free behavior, as follows:

u i f = Φ s i u i l + Φ d i q i ( 35
)
where the term Φ s i u i l comes from the equation ( 20) and characterizes the node static displacements and the term Φ d i q i is an additional term characterizing the body oscillatory behavior. Classically, the matrix Φ d i has the following form

Φ d i = u (1) i f . . . u (s) i f (36)
where the vector u (k) i f is the kth eigenmode associated with the equation:

M (11) i üi f + K (11) i u i f = 0 ( 37
)
in which the matrix K (11) i is defined in [START_REF] Balli | Transmission angle in mechanisms[END_REF] and the matrix M (11) i is the part of the matrix M i corresponding to the variables üi f . Thus, the dimension of the matrix M (11) i is (a × a). In Φ d i , s vectors are chosen among the a eigenmodes associated with (37). Usually, s << a (typically, in our examples, s is lower than 10). Now, using [START_REF] Sigmund | Topology optimization approaches: A comparative review[END_REF], the nodal displacement vector u i is thus given by

u i = u i f u i l = B i u i r , where B i = Φ s i Φ d i 1 b 0 (b×s) and u i r = u i l q i (38)
in which 1 b is the identity matrix of dimension b. Taking into account once again that the matrix B does not depend on the robot configuration or of the link deformation (hypothesis of small perturbations), the derivative of [START_REF] Saxena | Topology optimization of compliant mechanisms with strength consider[END_REF] with respect to time leads to

ui = B i ui r (39)
By using these transformations of coordinates, the body potential and kinetic elastic energies given at [START_REF] Liu | Kinematics, singularity and workspace of planar 5R symmetrical parallel mechanisms[END_REF] and [START_REF] Cammarata | Elastodynamic optimization of a 3T1R parallel manipulator[END_REF] becomes

U e i = 1 2 u T i r K i r u i r (40)
T e i = 1 2 uT i r M i r ui r (41)
where

K i r = B T i K i tot B i and M i r = B T i M i tot B i
are the reduced body stiffness and mass matrices. These matrices are usually of small dimension (typically, in our 2D problems, of dimension (12 × 12), while in 3D it is of dimension (18 × 18)) and they not depend on the robot configuration, which is a great advantage because they must be computed only once at each step of the optimization algorithm, whatever the number of tested configurations for the robot. Indeed, similarly as for the elastostatic model, most of the computational cost is due to the computation of the matrix Φ d i in [START_REF] Chen | The local optimum in topology optimization of compliant mechanisms[END_REF]. However, for two bodies i and k, the computation of the matrices Φ d i and Φ d k is independent. Thus the computation can be made in parallel on multi-core computers in order to save computational time.

Considering now the robot made of n bodies, the full potential and kinetic elastic energies of the system are given by:

U e = n i=1 U e i = 1 2 n i=1 u T i r K red i u i r = 1 2 u dyn T tot K dyn tot u dyn tot (42)
T e = n i=1 T e i = 1 2 n i=1 uT i r M red i ui r = 1 2 udyn T tot M dyn tot udyn tot (43)
where

• u dyn tot = [u T 1 r . . . u T n r ]
T is the vector stacking all vectors u i r for all n robot bodies

• K dyn tot and M dyn tot are block-diagonal matrices stacking on their diagonal all bodies stiffness and mass matrices as follows:

K dyn tot =             K 1 r 0 . . . 0 K n r             , M dyn tot =             M 1 r 0 . . . 0 M n r             (44)
The reduced robot elastodynamic model can be obtained by taking into account the fact the robot bodies are connected altogether through the interface nodes. As a result, the expression of the vector u dyn tot can be obtained from a reduced set of independent coordinates u d [START_REF] Shabana | Dynamics of Multibody Systems[END_REF][START_REF] Imbert | Analyse des structures par lments finis[END_REF] as follows:

u dyn tot = J d (q)u d (45)
where J d (q) depends on the robot configuration q but not on the nodal displacements u d (assumption of small perturbations). Moreover, for the analysis of the oscillatory free behavior, the configuration q is considered constant, thus

udyn tot = J d (q) ud (46)
Introducing ( 45) and ( 46) into ( 42) and ( 43), we get

U e = 1 2 u T d K d u d (47
)

T e = 1 2 uT d M d ud (48)
where

K d = J T d (q)K dyn tot J d (q) and M d = J T d (q)M dyn tot J d (q)
are the reduced robot stiffness and mass matrices.

Using the Lagrange equations, which state that, in absence of external efforts, we have

d dt ∂L ∂ u - ∂L ∂u = 0 (49)
where L = T e -U e , the reduced dynamic equation characterizing the robot free oscillations is

M d üd + K d u d = 0 (50)
A solution u d k of this equation satisfies:

ω 2 k M d -K d u d k = 0 (51)
where ω k = 2π f k , f k is the natural frequency associated with the kth natural mode of vibrations and u d k is its associated eigenvector.

For a 2D problem, the typical size of matrices K d and M d is lower than (12n × 12n), n being the number of robot bodies ((18n × 18n) for 3D problems). Therefore, even if the problems (F.1) and (51) will give exactly the same results of computation for the first frequencies, the resolution of the equation ( 51) is much more efficient than the resolution of the equation (F.1) due to the considerably reduced size of the problem.

Vibration performance index

The first robot natural frequency is used as a performance index:

c 2 = f 1 (52)
Its derivative with respect to the decision variables is given in Appendix F.

Selection of the robot configurations and trajectories for the computation of the performance criteria

All indices presented in Section 3 are configuration dependent (except for the case of the inertial parameters) and are thus local by nature. In the present section, we propose a methodology that can be used in order to enforce the chosen performance indices to be valid globally in the workspace or for a set of given trajectories. This methodology is based on the choice of optimal robot configurations or trajectories for which extreme performance will be attained. This procedure was inspired from [START_REF] Torii | Robust compliance topology optimization based on the topological derivative concept[END_REF][START_REF] Blankenship | Infinitely constrained optimization problems[END_REF].

Selection of the robot trajectories for the computation of the input efforts

As presented in the introduction, the main drawback of analyzing the dynamic performances with the manipulator input efforts is that they depend on the trajectory. In order to avoid this problem, one idea could be to make the manipulator moves on selecting "exciting" trajectories that make it exhibiting the highest values of the input efforts. To generate these exciting trajectories, it is necessary to define an optimization procedure.

Several methods for exciting trajectory generation exist [START_REF] Hollerbach | Handbook of Robotics -Model Identification[END_REF], the large majority of them being developed for dynamic parameter identification purpose. However, these methods have been developed in order to minimize the uncertainty of estimation of the parameters to be identified, which is a different goal from what we target, i.e. finding trajectories on which the designed robot will be near to attain its maximal dynamic capacities.

Therefore, we propose to adopt the following strategy, which is based on the assumption that we have an a priori estimation of the near-optimal design for the robot link that will be denoted as ρ. The choice of the variables ρ is discussed in Section 4.3. From Appendix B, we know that the input efforts are linear with respect to the decision parameters ρ (see Eq. (B.4)) and that, for our estimated optimal design ρ, they will be given by:

τ(t) = Γ(q(t), q(t), q(t)) ρ ( 53
)
where t is the current time instant, τ is the vector of robot input efforts, q, q, q are the robot active joint position, velocity and acceleration, respectively, and the matrix Γ is the Jacobian matrix of τ with respect to ρ.

Let us now denote as T a trajectory described by N samples of time, leading to N samples of the joint positions, velocities and accelerations (q, q, q), (q k , qk , qk ) being the value of (q(t k ), q(t k ), q(t k )) at the time instant t k . Considering all samples, we rewrite [START_REF] Svanberg | The method of moving asymptotes -a new method for structural optimization[END_REF] as

y = W ρ, where W stacks all samples of Γ such as W =             Γ(q 1 , q1 , q1 ) . . . Γ(q N , qN , qN )             (54)
and y stacks all samples of τ such as y = [τ(t 1 ) T . . . τ ( t N ) T ] T . Knowing a priori values for the actuator performance (such as the maximal velocity and motion ranges), the "exciting" trajectory will be found as the solution of the following optimization problem:

max T c j ( ρ) (for j = 4 or j = 5) under | qi (t k )| ≤ qi max q ∈ [q] and x ∈ [x] (55)
where the functions c 4 and c 5 are defined in (B.10) and (B.11) and represent the squared values of the infinite norms of the input efforts and of their root-mean-square along the trajectory T . Usually, we know that the infinite norm of the effort should be lower than the actuator peak effort specification, while the root-meansquare should be lower than the actuator continuous effort. Moreover,

• qi (t k ) is the value of the actuator i velocity at the sample time t k while qi max is the maximal admissible velocity for the actuator i.

• [q] represents the range of admissible joint configurations while [x] is the range of admission endeffector configurations.

In order to reduce the complexity of the design problem, we parameterize the trajectory T as a succession of portions of motions, each portion being parameterized by a typical motion profile (e.g. 5th order polynomial laws, bang-bang profiles, etc). Therefore, the decision variables of the problem (55) are restricted to

• the initial and final robot configurations for each portion of the trajectory T

• plus the duration of the motion profile corresponding to each portion of the trajectory T .

Selection of the robot configurations for the computation of elastic performance

Similarly as in the previous section, it is considered here that we have an a priori estimation of the nearoptimal design ρ. Therefore, the selection of the poses for the computation of the elastic performance will be found by solving the following optimization problems

• in the case of the static deformations, for a given loading f:

max q c 1 ( ρ, f) under q ∈ [q] and x ∈ [x] (56)
Note that several loadings f can be considered, resulting in the resolution of the problem (56) each time that the loading is changed.

• in the case of the natural frequencies:

min q c 2 ( ρ) under q ∈ [q] and x ∈ [x] ( 57
)
where the functions c 1 and c 2 are defined in ( 28) and ( 52), and represent the deformations under a loading f and the first robot natural frequency, respectively.

Discussion on the selection of near-optimal design parameters ρ

The selection of the optimal robot configurations and trajectories defined above depends on an a priori value of some near-optimal design parameters ρ. This is an issue as we should not have an idea of them before solving the design optimization problem.

A first approach would be to set ρ as the starting point of the optimal design process, i.e. to set all components in ρ at 1. This is of course not the best. The approach can be refined by using an iterative procedure:

1. First, set all components in ρ at 1. Let us denote this vector as ρ0 .

2. Then, solve the problems (55), ( 56) and ( 57) in order to find the exciting configurations q 0 and trajectories T 0 corresponding to the design ρ0 .

3. Run the design optimization problem (1) during v iterations (in our case study, 100 to 150 iterations where enough in our case). Save the final value for the design parameters into the variable ρnew .

4. Start again the steps 2 and 3 until the change in the exciting configurations q new and trajectories T new is "too big" (based on a criterion to be defined by the user).

We used this iterative procedure in what follows. In practice, we ran the previous algorithm during no longer than three iterations and it was enough in order to obtain a robust solution.

Case study: topology optimization of a five-bar mechanism

Definition of the optimal design problem in 2D

In this section, we propose as a case study the topology optimization of a five-bar mechanism (Fig. 3) which is a robot composed of five revolute joints, two of them being active (at points A 1 and A 2 ) while the other are passive. This mechanism with two dof is able to position in a plane its end-effector located at point C with coordinates (x y). The bodies between points A i and B i are called the proximal links and the bodies between points B i and C are called the distal links. We denote as "leg i" (i = 1, 2) the leg made with the joints located at A i , B i and C.

In what follows, we chose as link lengths the design parameters of the DexTAR robot [START_REF] Bourbonais | Minimum-time trajectory planning and control of a reconfigurable pick-and-place parallel robot[END_REF], which is a five-bar mechanism designed for high-speed pick-and-place operations. The link lengths are thus:

OA 1 = OA 2 = 0.1375 m, A 1 B 1 = A 2 B 2 = 0.23 m and B 1 C = B 2 C = 0.23 m,
where PQ is the distance between any point P and Q (Fig. 3).

The workspace of the DexTAR robot under consideration is plotted in Fig. 5 with the Type 2 singularity loci [START_REF] Gosselin | Singularity analysis of closed-loop kinematic chains[END_REF] corresponding to the robot working mode represented on the picture (branch index -1 for leg 1, +1 for leg 2 [START_REF] Bonev | Geometric analysis of parallel mechanisms[END_REF]).

In this workspace, we define a rectangular region, associated with the aforementioned branch indices for the robot legs, in which we decide to guarantee the robot performance performance as was done in [START_REF] Germain | Optimal Design of the IRSBot-2 Based on an Optimized Test Trajectory[END_REF]. This region is called the regular dextrous workspace [START_REF] Merlet | Parallel Robots[END_REF] and is defined here as a rectangle centred in (0 0.295) m of width equal to 0.25 m and of height equal to 0.15 m (Fig. 5). This region fixes the boundaries of the space search for finding the optimal robot configurations and trajectories for the computation of the robot performance.

We impose that both proximal (distal, resp.) links have the same shape for two main reasons:

• this will allow a symmetrical distribution of the mechanism performance in the dextrous workspace with respect to the y-axis

• the number of design variables is divided by two in this case, which decreases the computational time for the robot elastic performance.

Finally, we decide to solve the following optimization problem:

• the objective is to minimize the RMS of the actuator torques for any trajectory defined as follows:

the trajectories are parameterized by fifth-order polynomial motion profiles [START_REF] Khalil | Modeling, Identification and Control of Robots[END_REF],

the regular dextrous workspace is discretized with a grid made of 21 × 11 = 231 points, as shown in Fig. 5, and the initial and final points for the trajectories are selected among these 231 points the duration t f of the trajectory between the initial and final points is parameterized by the function t f = 0.75 p for p an integer, p ∈ {0, 1, 2, . . . , 8} (thus t f ∈ [0.1001, 1]) the exciting trajectory T * selected for being used in the optimal design problem are found by solving the optimization problem (55) when considering near-optimal design parameters ρ in which all components are equal to 1 for the initial optimization round, or equal to the optimized design parameters at the previous design optimization round.

Due to the symmetry imposed in the leg design and in the location and size of the dextrous operational workspace, we consider only the RMS of the actuator of leg 1.

• the constraints are:

to ensure that the elastic translational displacement δ C of the end-effector under both external loadings f 1 = [0 N 100 N 1 Nm] T and f 2 = [100 N 0 N -1 Nm] T (where the components of f i represent the force applied at C along the x and y-axes and the moment around z, respectively) is lower than δ max = 0.2 mm (in terms of norm) and that the elastic rotational displacement θ C of the distal link of leg 1 at C is lower than θ max = 0.25 mrad to ensure that the first natural frequency f 1 is bigger than f min = 180 Hz wherever in the dextrous regular workspace. The selection of the exciting configurations q δ i , q θ i (i = 1, 2) and q f for the computation of these performance indices is based on the methodology explained in Section 4.2 (where q δ i is the selected exciting configuration for the computation of the translational deformation when the loading f i is applied, while q θ i is the selected configuration for the computation of the rotational deformation for the same loading; q f is the configuration for the computation of the first natural frequency).

The design optimization problem can thus be formalized as

min ρ F = τ 2 1 (T * ) subject to g 1 = δ 2 C (f 1 , q δ 1 ) -δ 2 max ≤ 0 g 2 = θ 2 C (f 1 , q θ 1 ) -θ 2 max ≤ 0 g 3 = δ 2 C (f 2 , q δ 2 ) -δ 2 max ≤ 0 g 4 = θ 2 C (f 2 , q θ 2 ) -θ 2 max ≤ 0 g 5 = f min -f 1 (q f ) ≤ 0 (58)
In the next section, results obtained when applying this methodology are presented.

Results and discussion for the 2D model

Results

The initial design domain for the proximal and distal links is represented in Fig. 6. Each link is made of two holes (voids) at its extremity in which the joints will be inserted (axis of the motor at points A i , and passive revolute joints for all other holes). For the meshing of the links, QUA4 finite elements (i.e. four-nodes rectangular planar elements) of size 1 × 1 mm and thickness of 2 cm are used. The planar stress assumption is used. This assumption is valid because all links are considered coplanar and the efforts are applied in the links plane. If these hypotheses were not verified, other assumptions should have been taken into account [START_REF] Hao | Understanding coupled factors that affect the modelling accuracy of typical planar compliant mechanisms[END_REF].

Links are considered to be made of steel with Young's modulus E 0 = 210 GPa, Poisson's ratio ν = 0. For all optimization rounds, in order to obtain a smoother layout without checkerboards problem, a filtering procedure is used that consists in the modification of the density variables assigned to the elements with the information of its neighborhoods as was proposed in [START_REF] Bourdin | Filters in topology optimization[END_REF].

Then, we applied our approach in order to obtain the optimal design of our mechanism. All models and optimization algorithms have been encoded with Matlab in the Windows 7 environment.

We started with the selection of the configurations and trajectories for the computation of the robot performance. The selected configurations for the elastic performance computation are given in Tabs. We then ran the optimization algorithm during 100 iterations. The results in terms of design for the two links is shown in Figs. 7(c) and 7(d).

We then selected new configurations and trajectories based on the use of the near optimal design of the links obtained at the end of the first round. The selected configurations for the elastic performance computation are given in Tabs. 1 and 2 (row corresponding to round 2) while the trajectory of the computation of the torques is the same as the one presented in Fig. 8. We can see that the configuration q θ 2 for the computation of the rotational deformations when the wrench f 2 is applied changed. This is of critical importance for the global robustness of the final solution.

We then ran the optimization algorithm during 150 iterations and we stopped the computation. The results in terms of design for the two links is shown in Figs. 7(e) and 7(f).

Once again, we selected new configurations and trajectories based on the use of the near optimal design of the links obtained at the end of the first round. The new selected configurations for the elastic performance computation are given in Tabs. 1 and 2 (row corresponding to round 3) while the trajectory of the computation of the torques is still the same as the one presented in Fig. 8. We can see now that all configurations and trajectories are unchanged. Therefore we went for the final round of optimization.

The final results in terms of design for the two links is shown in Figs. 7(g) and 7(h). For this last round, there was no constraint on the number of maximal iteration, we waited for the convergence of the algorithm.

We consider that the algorithm converged when the maximal change between two sequential iterations for any component of the density vector ρ is lower than 0.01. The algorithm stopped after 597 iterations, with a maximal constraint violation of 8 • 10 -3 % (the maximal deformation in rotation in the whole workspace is equal to 0.250008 mrad, instead of 0.25 mrad which is negligible). In totality, the full procedure of optimization (including the selections of the trajectories and configurations) took 98 minutes with a Pentium 4 2.70 GHz, For this optimal design, the RMS of the torques along the trajectory shown in Fig. 8 is of 98 Nm while it was of 27 Nm for the final design, i.e. the torque RMS was divided by almost 4.

δ 1 q δ 2 q θ 1 q θ 2 x [m] y [m] x [m] y [m] x [m] y [m] x [m] y [m]
The values of the deformations and natural frequencies within the dextrous workspace are shown in Fig. 9.

We can see that the deformations (first natural frequency, resp) are lower (bigger, resp.) than the defined acceptable values, i.e. that our approach leaded to a robust optimization of the robot performance over the whole workspace.

In order to show the importance of the update of the value of q θ 2 between Round 1 and Round 2, we optimized the robot for constraints calculated for the configurations and trajectories selected at the initial round only (row corresponding to round 1 in Tabs. 1 and 2 for the configurations, Fig. 8 for the trajectory).

The algorithm stopped after 824 iterations, and 94 minutes of computation on the same computer as previously.

The final optimized links are shown in Fig. 10 and the values of the deformations and natural frequencies are depicted in Fig. 11. It can be shown that the desired boundaries on the deformations are not respected (the maximal deformation in rotation in the whole workspace is equal to 0.35 mrad, which is far beyond the limit of 0.25 mrad -see Fig. 11(c)). This shows that it is critical to update the selected worst case configuration / trajectory during the solving process.

Discussions

Comparison of LM and MMA in terms of convergence time. In order to show the efficiency of the Linearization Method (LM) for this class of problem, we decided to solve the optimization problem [START_REF] Pshenichnyj | The Linearization Method for Constrained Optimization[END_REF] with another optimization method known for its computational efficiency: the Method of Moving Asymptotes (MMA) [START_REF] Svanberg | The method of moving asymptotes -a new method for structural optimization[END_REF].

We used the Matlab code kindly provided by Prof. K. Svanberg and ran the optimization problem ( 58) by defining the same stopping conditions as previously. The MMA algorithm stopped after 1785 iterations. In totality, the full procedure of optimization (including the selections of the trajectories and configurations) took 143 minutes on the mentioned computer, so almost 45 % of additional time with respect to LM, however, the final torque RMS was of 25.9 Nm, i.e. MMA proposes a robot whose torque RMS on the exciting trajectory is lower than with LM, but the difference is very low (1 Nm, i.e. almost 1 % of the final objective). The results in terms of robot link design with MMA are shown in Fig. 12. Surprisingly, they are sensibly different, which seems to indicate that several local minima exist.

The evolution of the torque RMS as a function of the computational time in the last optimization round is shown in Fig. 13. It can be seen than both functions reach the same values of the objective in almost the same time, even if MMA goes a bit lower (but in a more longer time) than LM. Thus, LM is a valuable and computationally-efficient method for topology optimization.

Computational efficiency of the reduction techniques. We compared the time requested for the computation of the constraints when: • they are calculated with the model reduction techniques introduced in Section 3

• they are calculated without model reduction techniques, as shown in Appendix C and Appendix D.

On the aforementioned computer, the computation of the constraints when using model reduction techniques took 4.19 sec. When using the full models, it took 8.22 sec. In order to better understand the benefits of the model reduction technique, we plotted in Fig. 14 the evolution of computational cost for the computation of the first natural frequency as a function of the number of tested configurations. For the computation of a single configuration, the computational costs are equivalent (4.41 sec for the model reduction technique, 4.19 sec without it). The benefits appear when we compute the natural frequency for two configurations: in this case, the total computation time with the model reduction technique is of 4.52 sec, to be compared with 8.37 sec in order to compute the frequency without reduction technique. For twenty configurations, the computation with model reduction techniques took 6.44 sec, instead of 85 sec without it. Therefore, using model reduction techniques allows the considerable decrease of the computational cost in robot design optimization.

Solving the problem in 3D

The optimization problem can also be solved in 3D, even if in such a case the computational time is increase. In Fig. 15, we show the proximal and distal links of a five-bar mechanism which have been obtained when solving the following optimization problem:

• the objective was to minimize the robot mass

• the constraints are to ensure that the elastic translational displacement δ C of the end-effector under the external loading f = [0 N 50 N 10 N 0 Nm 0 Nm 1 Nm] T (where the components of f represent the force applied at C along the x, y and z-axes and the moment around the x, y and z-axes, respectively) is lower than δ max = 0.025 mm (in terms of norm) and that the elastic rotational displacement θ C of the distal link of leg 1 at C is lower than θ max = 0.25 mrad wherever in the dextrous regular workspace. The selection of the exciting configurations q δ and q θ for the computation of these performance indices was the based in the methodology explained in Section 4.2.

The design optimization problem can thus be formalized as

min ρ F = m robot under g 1 = δ 2 C (f, q δ ) -δ 2 max ≤ 0 g 2 = θ 2 C (f, q θ ) -θ 2 max ≤ 0 (59)
Here, we chose as link lengths the design parameters of the DexTAR robot [START_REF] Bourbonais | Minimum-time trajectory planning and control of a reconfigurable pick-and-place parallel robot[END_REF] from the company Mecademic. The link lengths are thus: (c) Norm of the translational deformations under the application of the loading f 2 (Fig. 3). The selected dextrous regular workspace is a rectangle of dimension 12 × 4 cm centred in x = 0 m and y = 0.11 m.

OA 1 = OA 2 = 0.056 m, A 1 B 1 = A 2 B 2 =
The initial design domain for the proximal and distal links is represented in Fig. 16 (planar view, the initial width of both links is of 20 mm). For the meshing of the links, HEXA8 finite elements (i.e. eight-nodes cubic planar elements) of size 1.3 × 1.3 × 1.3 mm are used. Links are considered to be made of steel with Young's modulus E 0 = 210 GPa, Poisson's ratio ν = 0.3 and density of 7800 kg/m 3 . As a result, 75200 elements are used for meshing the proximal links while the distal links are made of 21760 elements.

We consider that the algorithm converged when the maximal change between two sequential iterations for any component of the density vector ρ is lower than 0.05. The algorithm stopped after 4 hours and 42 minutes of computation with a Pentium 4 2.70 GHz, 24 GB of RAM.

In order to accelerate the computation of the elastostatic properties by avoiding the inversion of the matrices K (11) i in [START_REF] Hubert | Static of parallel manipulators and closeness to singularity[END_REF] for both links, we computed the matrices Φ s i using a "pcg" (Preconditioned Conjugate Gradient) algorithm [START_REF] Barrett | Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods[END_REF]. Indeed, we computed each column of the matrix Φ s i by solving the system of equations with the "pcg" algorithm: K (11) i Φ s i j = -K (12) i j [START_REF] Combettes | Method of successive projections for finding a common point of sets in metric spaces[END_REF] in which Φ s i j (K (12) i j , resp.) is the jth column of Φ s i (K (12) i , resp.). The computation can be accelerated by doing a computation in parallel of each column of the matrix Φ s i (indeed, we used 6 CPUs in parallel on our computer to do so).

We also introduced in the optimization problem some constraints on the limits of the first natural frequency. The results obtained were similar to the ones shown in Fig. 15, but with much longer computational time (around 6 days and 8 hours). This considerable increase of the computational cost was due to the calculation of the matrices Φ d i in [START_REF] Chen | The local optimum in topology optimization of compliant mechanisms[END_REF] for both links which is obtained by solving an eigenvalue problem of large dimension (matrices of size (256581 × 256581) for the proximal links). Furthermore, it should be mentioned that this eigenvalue problem cannot be solved in Matlab with the classical "eig" or "eigs" functions which requests too much memory allocation (our computer went out of memory), but with a "lobpcg" (Locally Optimal Block Preconditioned Conjugate Gradient) algorithm [START_REF] Knyazev | Toward the optimal preconditioned eigensolver: Locally optimal block preconditioned conjugate gradient method[END_REF] which is known to be efficient for computing eigenvalues problems of very large dimensions.

Conclusion

Topology optimization aims at optimizing the material distribution in a link or a set of links in order to satisfy performance criteria related to the link elastic behavior. Recent works introduced topology optimization in the design of robots, but the proposed methodology led to a local optimization of the robot performance (for a few trajectories or configurations). Moreover, most of performance indices used are not in strong relation with easy-to-understand technological requirements.

In the present paper, our primary contribution was the proposition of a methodology that was able to perform a topology optimization for robots, valid globally in the workspace or for a set of given trajectories, which was based on the use of technology-oriented performance criteria. In order to enforce the chosen performance indices to be globally valid, optimal robot configurations or trajectories for which extreme performance will be attained are computed, and iteratively updated.

In order to decrease the computational time associated with the computation of performance indices in numerous configurations, we exploited the structure of the elastic models in order to reduce their computational complexity. Indeed, we showed that it was possible to use configuration-independent model reduction techniques in order to considerably decrease the size of the stiffness and mass matrices of each link, expressed in their local frame. Then, these reduced matrices are used in order to build the configuration-dependent robot elastic models which are shown to be computationally efficient.

We used an optimization algorithm called the Linearization Method which has proven, for our class of problems, to give results in a computational time equivalent to standard topology optimization algorithms but its implementation was less complex and made it quite easy to perform modification or improvement.

The methodology was successfully applied for the design of a five-bar mechanism and we showed that our approach leaded to a robust optimization of the robot performance over the whole workspace. We also showed that the Linearization Method algorithm was as competitive as the standard algorithms for topology optimization in terms of convergence time, such as the Method of Moving Asymptotes.

Finally, we would like to mention that the proposed approach is not software dependent: we used Matlab for finite element modeling for reasons of simplicity, but co-simulation with any FEA software able to communicate with Matlab (or another software allowing the encoding of the optimization procedure) could be used.

Appendix A. Inertial parameters

Inertial parameters are performance indices which are global by nature, as they do not depend on the robot configuration. As known from [START_REF] Khalil | Modeling, Identification and Control of Robots[END_REF], each robot body can be characterized by ten (rigid) inertial parameters: its mass, three first moments of inertia (also called static moments), three moments of inertia and three products of inertia.

Let us consider the body i meshed with m elements. For the element j of the body i, denoted as the element i j, we define as:

• x i j , y i j and z i j the position of the origin M i j of the element i j in the local frame (O i , x i , y i , z i ) attached to the body i (Fig. 4).

• ρ i j the density associated with the element i j,

• m i j the mass of this element, As a result, and from [START_REF] Khalil | Modeling, Identification and Control of Robots[END_REF], we have the following formulas for the computation of the inertial parameters:

• body's mass: m i = m j=1 m i j ρ i j ,

• body's first moments of inertia: mx i = m j=1 m i j x i j ρ i j , my i = m j=1 m i j y i j ρ i j and mz i = m j=1 m i j z i j ρ i j ,

• body's moments of inertia:

xx i = m j=1 m i j (y 2 i j + z 2 i j )ρ i j , yy i = m j=1 m i j (x 2 i j + z 2 i j )ρ i j , zz i = m j=1 m i j (x 2 i j + y 2 i j )ρ i j ,
• body's products of inertia: xy i = m j=1 m i j x i j y i j ρ i j , xz i = m j=1 m i j x i j z i j ρ i j , yz i = m j=1 m i j y i j z i j ρ i j .

Grouping all inertial terms for the body i in the vector χ i = [m i mx i my i mz i xx i yy i zz i xy i xz i yz i ] T , we finally have the following expression:

χ i = Ψ i ρ i , where Ψ i = m i j                                                 1 1 . . . 1 x i1 x i2 . . . x im y i1 y i2 . . . y im z i1 z i2 . . . z im y 2 i1 + z 2 i1 y 2 i2 + z 2 i2 . . . y 2 im + z 2 im x 2 i1 + z 2 i1 x 2 i2 + z 2 i2 . . . x 2 im + z 2 im x 2 i1 + y 2 i1 x 2 i2 + y 2 i2 . . . x 2 im + y 2 im x i1 y i1 x i2 y i2 . . . x im y im x i1 z i1 x i2 z i2 . . . x im z im y i1 z i1 y i2 z i2 . . . y im z im                                                 and ρ i =                  ρ i1 ρ i2 . . . ρ im                  (A.1)
in which m i j is a constant mass as we consider here that the body is meshed with equal size elements. The matrix Ψ i is the Jacobian of χ i with respect to ρ i .

Grouping all inertial terms for all bodies in the vector χ = [χ T 1 . . . χ T n ] T (n being the number of robot bodies), we finally get the expression:

χ = Ψ ρ, where Ψ =             Ψ 1 . . . 0 . . . . . . . . . 0 . . . Ψ n             and ρ =             ρ 1 . . . ρ n             (A.2)
Any linear combination of the inertial parameters can be then obtained by multiplying the vector χ by a row vector a leading to a performance criterion c 1 equal to

c 3 = aχ = aΨ ρ (A.3)
The performance criterion c 3 is linear with the decision variables ρ and the computation of its derivative with respect to ρ is straightforward.

Note that the linear combinations allowing the computation of the grouped inertial parameters with the most effects on the dynamics of the robot can be found using the methods [START_REF] Gautier | Direct calculation of minimum set of inertial parameters of serial robots[END_REF][START_REF] Gautier | Numerical calculation of the base inertial parameters[END_REF].

• the infinite norm τ i ∞ of the input effort for the motor i (τ i being the ith component of the vector τ)

along the trajectory (q, q, q) that should be lower than the actuator peak effort specification,

• the root-mean-square τ i of the input effort for the motor i along the trajectory (q, q, q) than should be lower than the actuator continuous effort specification.

In order to compute these indices, let us sample the trajectory (q, q, q) at a given frequency, leading to N samples for (q, q, q). We denote as y i the vector containing all values of τ i computed along the sampled trajectory (q, q, q): the kth component of y i is equal to τ (k) i = Γ i (q k , qk , qk ) ρ, where (q k , qk , qk ) is the value of (q, q, q) at the sample k and Γ i is the ith row of Γ. As a result, y i takes the form:

y i = W i ρ, where W i =             Γ i (q 1 , q1 , q1 ) . . . Γ i (q N , qN , qN )             (B.7)
By definition, we have

τ i ∞ = max √ y i * y i (B.8)
in which the operator " * " represents the term-by-term product and

τ i = y T i y i N (B.9)
In order to avoid conserving square roots in the expressions of the performance criteria, it is better to use as performance index the square of the infinite norm τ i ∞ and of the root-mean-square τ i . As a result, by using (B.7), the selected performance criteria based on the input efforts are:

c 4 = τ i 2 ∞ = max(y i * y i ) = max (W i ρ * W i ρ) (B.10) c 5 = τ 2 i = 1 N y T i y i = 1 N ρ T W T i W i ρ (B.11)
The derivatives of these quadratic forms with respect to the decision variables is still straightforward. Obviously, all these criteria will differ depending on the chosen trajectory (q, q, q). Discussion on the selection of an optimal trajectory is the aim of the Section 4.

Appendix C. Usual computation of the robot elastostatic model

Considering the robot made of n bodies, and by using the expression of the body potential energy (13), the full potential elastic of the system is given by:

U e = n i=1 U e i = 1 2 n i=1 u T i K i u i = 1 2 u T tot K tot u tot (C.1)
where

• u tot = [u T 1 . . . u T n ]
T is the vector stacking all nodal displacements for all n robot bodies

• K tot is a block-diagonal matrix stacking on its diagonal all bodies stiffness matrices as follows:

K tot =             K 1 0 . . . 0 K n             (C.2)
The robot stiffness matrix K can be obtained by taking into account the fact the robot bodies are connected altogether through the interface nodes. As a result, the expression of the vector u tot can be obtained from a reduced set of independent coordinates u as follows [START_REF] Shabana | Dynamics of Multibody Systems[END_REF][START_REF] Imbert | Analyse des structures par lments finis[END_REF]:

u tot = J(q)u (C.3)
where J(q) depends on the robot configuration q but not on the nodal displacements u (assumption of small perturbations).

Introducing (C.3) into (C.1), we get

U e = 1 2 u T Ku (C.4)
where K = J T (q)K tot J(q) is the robot stiffness matrix which relates the nodal displacements u to external forces f exerted on the nodes by the relation [START_REF] Imbert | Analyse des structures par lments finis[END_REF]:

f = ∂U e ∂u = Ku (C.5)

Appendix D. Usual computation of the robot natural frequencies

Considering the robot made of n bodies, and by using the expression of the body kinetic energy (31) the full kinetic elastic of the system is given by:

T e = n i=1 T e i = 1 2 n i=1 uT i M i ui = 1 2 uT tot M tot utot (D.1)
where M tot is a block-diagonal matrix stacking on its diagonal all bodies mass matrices as follows:

M tot =             M 1 0 . . . 0 M n             (D.2)
Differentiating (C.3) with respect to time, and keeping in mind that for the robot frequency analysis, the configuration q is nominal is considered as constant: utot = J(q) u (D.3) Introducing (D.3) into (D.1), we get

T e = 1 2 uT M u (D.4)
where M = J T (q)M tot J(q) is the robot mass matrix.

Using the Lagrange equations, which state that, in absence of external efforts, we have

d dt ∂L ∂ u - ∂L ∂u = 0 (D.5)
where L = T e -U e in which the expression of U e comes from (C.1), we get the dynamic equation characterizing the robot free oscillations M ü + Ku = 0 (D.6)

A solution u t k of this equation satisfies:

ω 2 k M -K u t k = 0 (D.7)
where ω k = 2π f k , f k is the natural frequency associated with the kth natural mode of vibrations and u t k is its associated eigenvector.

Appendix E. Derivative of the elastostatic performance index with respect to the decision variables

The derivative of the criterion [START_REF] Hess-Coelho | Decoupling of dynamic equations by means of adaptive balancing of 2-dof open-loop mechanisms[END_REF] with respect to the decision variable ρ i j is equal to:

∂c 1 ∂ρ i j = 2u T c e ∂u r ∂ρ i j (E.1)
The term ∂u r ∂ρ i j can be obtained by the differentiation of the equation ( 27) with respect to the decision variable ρ i j : 0 = ∂ ∂ρ i j (K r u r ) = ∂K r ∂ρ i j u r + K r ∂u r ∂ρ i j (E.2) leading thus to ∂u r ∂ρ i j = -K -1 r ∂K r ∂ρ i j u r (E.

3)

The computation of the derivative of the matrix K r with respect to the decision variable is not computationally efficient (problem of data storage due to the large number of decision variables), but, by analyzing the expression (E.3), we see that it is not necessary to compute ∂K r ∂ρ i j but the product ∂K r ∂ρ i j u r which can be obtained in a more efficient way. Indeed, by using [START_REF] Chaudhary | Dynamics and balancing of multibody systems[END_REF], this product becomes equal to ∂K r ∂ρ i j u r = J T r (q)

                            0 . . . ∂K red i ∂ρ i j . . . 0                             J r (q)u r = J T r (q)                             0 . . . ∂K red i ∂ρ i j . . . 0                             u red tot (E.4)
or also ∂K r ∂ρ i j u r = J T r (q)

                            0 . . . ∂K red i ∂ρ i j u i l . . . 0                             (E.5)
where ∂K red i ∂ρ i j is the ith block term on the diagonal, all other terms being null.

Then, by differentiating the equation ( 21), it is possible to compute the term ∂K red i ∂ρ i j u i l in (E.5) as follows:

∂K red i ∂ρ i j u i l = ∂ ∂ρ i j K (21) i Φ s i + K (22) i u i l = ∂K (21) i ∂ρ i j Φ s i u i l + K (21) i ∂Φ s i ∂ρ i j u i l + ∂K (22) i ∂ρ i j u i l (E.6)

The derivatives of K (21) i and K (22) i , which are nothing else than the derivatives of the terms of the matrix K i given in [START_REF] Balli | Transmission angle in mechanisms[END_REF], can be found from the differentiation of the matrix K i given by:

∂K i ∂ρ i j = p ρ p-1 i j E 0 J T i                           0 . . . K (0) i j . . . 0                           J i (E.7)
where K (0) i j is the jth block term on the diagonal, all other terms being null.

The other terms K (21) i ∂Φ s i ∂ρ i j in (E.6) can be found from the differentiation of the expression K (11) i Φ s i = -K (12) i deduced from (20)

∂K (11) i ∂ρ i j Φ s i + K (11) i ∂Φ s i ∂ρ i j = -∂K (12) i ∂ρ i j ⇒ ∂Φ s i ∂ρ i j = -K (11) i

-1       
∂K (11) i ∂ρ i j Φ s i + ∂K (12) i

∂ρ i j        (E.8)
in which the derivatives of K (11) i and K (12) i can be found from the partitioning of the matrix ∂K i /∂ρ i j in (E.7), using the same partitioning as for the matrix K i of [START_REF] Balli | Transmission angle in mechanisms[END_REF].

Finally, by multiplying this equation by K (21) i (which is equal to K (12) T i ) and by using the expression of Φ s i given in [START_REF] Hubert | Static of parallel manipulators and closeness to singularity[END_REF], we found

K (21) i ∂Φ s i ∂ρ i j = Φ T s i       
∂K (11) i ∂ρ i j Φ s i + ∂K (12) i ∂ρ i j

       (E.9)
Appendix F. Derivative of the elastodynamic performance index with respect to the decision variables

From [START_REF] Huang | Evolutionary topological optimization of vibrating continuum structures for natural frequencies[END_REF], we know that derivative of the criterion [START_REF] Pshenichnyj | The Linearization Method for Constrained Optimization[END_REF] with respect to the decision variable ρ i j can be obtained as follows. First, let us consider the Rayleigh quotient which is given by:

ω 2 k = u T d k Ku d k u T d k Mu d k (F.1)
Expanding (F.1) for the first mode of vibration, we have

u T d 1 K d u d 1 -ω 2 1 u T d 1 M d u d 1 = 0 (F.2)
By introducing ( 12), ( 36) and ( 46) into (F.2), we know that

u all = J d BJ tot u d 1 = J all u d 1 , where B =             B 1 0 . . . 0 B n             and J tot =             J 1 0 . . . 0 J n             (F.3)
and that K d = J T all K all J all , M d = J T all M all J all (F.4) in which J all = J d BJ tot and u all = [u T 1 tot . . . u T n tot ] T groups all nodes variables (u i tot being defined in (10)), and where

K all =             K 1 tot 0 . . . 0 K n tot             and M all =             M 1 tot 0 . . . 0 M n tot             (F.5)
Thus, introducing all these expressions in (F.2), we get u T all K all u all -ω 2 1 u T all M all u all = 0 (F.6)

Then, differentiating (F.6) with respect to ρ i j , we get the derivative of ω 1 with respect to ρ i j as: . . .

∂ω 1 ∂ρ i j = 1
0                           (F.8)
in which the expressions of ∂K i tot ∂ρ i j and ∂M i tot ∂ρ i j are given by

∂K i tot ∂ρ i j = p ρ p-1 i j E 0                           0 . . . K (0) i j . . . 0                           (F.9)
and

∂M i tot ∂ρ i j =                           0 . . . M (0) i j . . . 0                           (F.10)
where K (0) i j and M (0) i j are defined at (8) and ( 29), respectively, and they are the jth block terms on the diagonal, all other terms being null.

As, at the natural frequency f 1 , K all -ω 2 1 M all u all = 0, and taking into account that f 1 = 2πω 1 , we finally have:

∂ f 1 ∂ρ i j = 1 4πω 1 u T all M d u all u T
all ∂K all ∂ρ i j u all -ω 2 1 u T all ∂K all ∂ρ i j u all (F.11)

  projected on the linearization of the constraint to obtain v k . Then, d k is the vector v k -ρ k , depicted as a thick arrow in Figure 2. We see that the projection process modifies the cost function gradient (the thin arrow) so that the new direction (the thick arrow) accounts the constraint. The scalar s k is used as a second step size for the method. The right hand side diagram of Figure 2 shows the vectors d k obtained for different vectors ρ k , and clearly illustrates that the method converges

Figure 3 :

 3 Figure 3: A five-bar mechanism (the gray pairs denote the actuated joints).

Figure 5 :

 5 Figure 5: The five-bar mechanism under consideration and its workspace (to scale).

  3 and density of 7800 kg/m 3 . As a result, 20932 elements are used for meshing the proximal links while the distal links are made of 14208 elements (Figs. 7(a) and 7(b)).

Figure 6 :

 6 Figure 6: Initial design domain for the design of the robot link (to scale). All dimensions are in millimeters.

  1 and 2 (row corresponding to round 1) while the trajectory of the computation of the torques is given in Fig. 8. For this first round, we use links made of the full initial design domain for the computation of the elastic and dynamic performance (Figs. 7(a) and 7(b)).

Figure 7 :Figure 8 :

 78 Figure 7: Evolution of the design of the five-bar links at the end of each optimization round: the links are shown in gray-scale (black elements correspond to ρ i j = 1, white elements to ρ i j = 0, and gray elements to 0 < ρ i j < 1)

  0.09 m and B 1 C = B 2 C = 0.09 m Norm of the translational deformations under the application of the loading f 1 Absolute value of the rotational deformations under the application of the loading f 1

Figure 9 :Figure 10 :yFigure 11 :Figure 12 :Figure 13 :

 910111213 Figure 9: Evolution of the deformations and of the first natural frequency in the workspace for the optimized design

  configurations Computational time [sec] with model reduction technique without model reduction technique

Figure 14 :

 14 Figure 14: Evolution of computational cost for the computation of the first natural frequency as a function of the number of tested configurations.

Figure 15 :

 15 Figure 15: Optimal designs of the five-bar links using 3D models

Figure 16 :

 16 Figure 16: Initial design domain for the design of the robot link (planar view, to scale). The initial width of both links is of 20 mm. All dimensions are in millimeters.

Table 1 :

 1 Configurations (positions of the end-effector) used for the computation of the elastostatic performance Const. for translational displacements Const. for rotational displacements q

Table 2 :

 2 Configurations (positions of the end-effector) used for the computation of the elastodynamic performance

		q f
		x [m]	y [m]
	Round 1	0	0.25
	Round 2	0	0.25
	Round 3	0	0.25
	16 GB of RAM.

Note that providing the diagonal of the Hessian to MMA allows solving this problem, but this second order information is not available in the problems solved in the present paper.

Acknowledgements

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

The assistance on the MMA algorithm from Professor K. Svanberg is gratefully acknowledged.

Appendix B. Input efforts and energy consumption

From [START_REF] Khalil | Modeling, Identification and Control of Robots[END_REF], we know that the input efforts, power and energy functions of any robot (assuming that the robot is rigid) are linear with respect to the inertial parameters, i.e. we have: τ = Φ(q, q, q)χ (B.1)

where τ is the vector of robot input efforts, P is the power delivered by the actuators, H is the total energy of the system, q, q, q are the robot active joint position, velocity and acceleration, respectively, and the matrix Φ(q, q, q) is the Jacobian matrix of τ with respect to χ (and is also called the observation matrix) while h(q, q)

is the Jacobian matrix of H with respect to χ. By introducing (A.2) into (B.1) and (B.3), we obtain τ = Γ(q, q, q) ρ, where Γ(q, q, q) = Φ(q, q, q)Ψ (B.4) P = α(q, q, q) ρ, where α = qT Γ(q, q, q) = qT Φ(q, q, q)Ψ (B.5) H = β(q, q) ρ, where β(q, q) = h(q, q)Ψ (B.6)

The power P and energy H could be used directly as performance criteria. However, two different (but complementary) performance indices could be defined from the expression of the input efforts: