Mohamed Amine Rahmouni
email: mohamedamine.rahmouni@cea.fr

Richard Bearee
email: richard.bearee@ensam.eu

Mathieu Grossard
email: mathieu.grossard@cea.fr

Eric Lucet
email: eric.lucet@cea.fr

Vibration reduction control for redundant flexible robot manipulators

Keywords: Vibration Suppression, Redundancy Resolution, Flexible Robot Manipulators, Feedback Control Design

In this paper, both vibration reduction and task accuracy improvement for flexible redundant manipulators, using the resolution of redundancy, is addressed. The exciting force of flexural motion, which is induced by the motion of the manipulator, causes undesired deflection at the end-effector. The null space inherent in redundancy feature is exploited to damp out the vibration while maintaining the end-effector posture. In order to exploit all the degrees of freedom in the null space, a feedback control design based on torque optimization has been introduced using analytical method for the redundancy resolution. The simulation analysis is presented to test the effectiveness of the proposed approach.

INTRODUCTION

The robotic vision of the advanced industry pays more attention to flexible manipulators, which have high load to mass ratio. However, the presence of flexible substructures, joints and/or links, may induce oscillator behavior which deteriorates the accuracy of the end-effector position. In addition, the industrial application including robot manipulation process are in need of robot manipulators which are able to operate in various constrained environments. In this context, the kinematic redundancy of redundant robots is highly investigated. The fact that it is offering more degrees of freedom (DOFs) regarding to those which are required to execute a given task, leads to several possible configurations in joints space to hold the same Cartesian posture of the end-effector. In literature, this internal motion is referred to as self-motion or null-space. The redundancy is used to deal with different kinematic and dynamic issues related to robot manipulators, in order to enhance their performances. For instance, in [START_REF] Xinyu | Obstacle avoidance for kinematically redundant robot[END_REF] and [START_REF] Khatib | A unified approach for motion and force control of robot manipulators: The operational space formulation[END_REF], the redundancy is used to deal with obstacle avoidance and joints constraints.

The control of redundant flexible robot manipulator is an interesting research topic, since a vast investigating is in progress to date. In order to achieve a given optimal performance, in [START_REF] Benosman | Control of flexible manipulators: A survey[END_REF], the dynamic of flexible elements is included in the generalised dynamic model of the flexible manipulator, and used to design various controllers algorithms. The input shaping technique, which is based on pre-estimated frequencies of the vibration modes, generates a smooth trajectory so that the flexural dynamic will be reduced, see [START_REF] Sahinkaya | Input shaping for vibration-free positioning of flexible systems[END_REF]. Since, different shaping filters are developed, such as exponential filter in [START_REF] Biagiotti | Optimal trajectories for vibration reduction based on exponential filters[END_REF], and damped jerk-limited profile in [START_REF] Béarée | New Damped-Jerk trajectory for vibration reduction[END_REF]. On the other hand, as cited in [START_REF] Benosman | Control of flexible manipulators: A survey[END_REF], the active control schemes are established similarly to those found in nonlinear and linear control theory for dynamic system, such as PID control, singular perturbation control, and so on. The mentioned control designs are developed for both redundant and non-redundant manipulators. For example in [START_REF] Nakanishi | Operational space control: A theoretical and empirical comparison[END_REF], the operational space control design is proposed for the redundant manipulators, in velocity and acceleration levels. Although, a little attention has been paid to take advantage of the redundancy in the controller design and a few works have introduced this ability, for instance in [START_REF] Sadeghian | Dynamic multi-priority control in redundant robotic systems[END_REF], a control design in dynamic level for multi-priority control in redundant space is presented.

The joint motion is capable to affect the vibrations since the dynamic of the rigid body and the flexible parts of the manipulator are highly coupled. This ability may be used to achieve efficiently smooth motion without exciting the vibration dominant modes. The redundancy ability to affect the manipulator flexibility has been introduced and examined in some works. In [START_REF] Yue | Point-to-point trajectory planning of flexible redundant robot manipulators using genetic algorithms[END_REF] and [START_REF] Gao | Research on Joint Torque Optimization Method of Redundant Space Manipulators with Vibration Suppression[END_REF], the null-space of the joint acceleration has been used to achieve an optimal joints motion so that the vibration are reduced. Nevertheless, the optimization algorithm consumes a considerable time, and even in some cases the optimal solution for joints motion is irregular and difficult to execute on-line. In [START_REF] Nguyen | Dynamic control of flexible, kinematically redundant robot manipulators[END_REF] and [START_REF] Bian | Study on vibration reduction and mobility improvement for the flexible manipulator via redundancy resolution[END_REF], authors have established a control design based on the Jacobian null-space, which is used explicitly to damp out the vibration modes. It has been proved that the self-motion gets effectively rid of the manipulator flexibility effects only when the number of redundancy degrees is the same as the number of selected modes to be controlled. Similarly, in [START_REF] Hirano | Simultaneous control for endpoint motion and vibration suppression of a space robot based on simple dynamic model[END_REF], a flexural state estimator using force/torque sensor is added this time to the control design to reduced the induced vibration of flexible appendage in space robot. In [START_REF] Sasaki | Vibration suppression control of redundant manipulator with flexible structure by considering nullspace motion[END_REF], the resonance ratio control is designed to suppress the vibrations. The work space observer is used to estimate the flexural reactive force, which is injected afterward in the Jacobian null-space in addition to the joint local control, by assuming that the joints velocities are low. In [START_REF] Hishinuma | Singularity-consistent vibration suppression control with a redundant manipulator mounted on a flexible base[END_REF], the vibration reduction is guaranteed by resolving the kinematic redundancy with the help of an additional constraint which is obtained from vibration dynamics, this approach is known as the task-space augmentation approach which is based on the augmented Jacobian. In this case, the joints motion solution may become unfeasible around algorithm singularity.

However, in the methods mentioned above, all the redundant DOFs were devoted to damp out the vibration modes in aim to ameliorate the accuracy of the end-effector. Nevertheless, realizing more than one performance in same time via redundancy resolution, would be highly appreciated, specially in case of hyper-redundant flexible manipulator. Thus, in the proposed control algorithm, we exploit all redundancy DOFs to compensate for the nonlinear dynamics in motion, and to damp out the induced vibration as well. These two stages control approach which is investigated in the following analysis, has not yet presented in literature.

In this paper, we analyze the ability of using the redundancy to optimize and carry out multiple performance criteria for flexible redundant manipulators. A control algorithm including the local redundancy resolution scheme for flexible serial manipulators is proposed. The motion planning of the manipulator is generated so that the nonlinear dynamics are compensated, while the dominating vibratory modes are damped out. Numerical simulations are conducted on a 5 degree of freedom planar manipulator to prove the effectiveness of the proposed approach.

GENERALIZED DYNAMIC MODEL

Dynamic analysis

The dynamic structure of the flexible manipulator is the combination of both the joints motion and the flexural motion, i.e the elastic deformation in the flexible elements, which is depicted by n r rigid DOFs and n f flexural DOFs respectively. Based on Lagrange's formulation and assumed mode method, the dynamic model of the flexible manipulator with flexible links is generally described by n dynamics equations of motion (for more details see [START_REF] Gao | Research on Joint Torque Optimization Method of Redundant Space Manipulators with Vibration Suppression[END_REF] and [START_REF] Bian | Impact vibration attenuation for a flexible robotic manipulator through transfer and dissipation of energy[END_REF])

M(Θ) Θ + D( Θ, Θ) Θ + K(Θ) Θ = Γ + Γ (1)
where n = n r + n f is the total DOFs of the system, Θ = [q T r q T f ] T ∈ IR n is the generalized coordinates vector, q r ∈ IR n r and q f ∈ IR n f are respectively the joint coordinates vector and the elastic displacement coordinates vector. M(Θ) ∈ IR n×n is the global mass matrix, D( Θ, Θ) ∈ IR n×n denotes the global damping matrix, K(Θ) ∈ IR n×n is the global stiffness matrix. Γ ∈ IR n is the generalized torques vector, which depends on the control torques. Γ ∈ IR n is the vector which aggregates the Coriolis, gravitational and centripetal forces in addition to other nonlinear dynamics including the externally applied forces and the non-conservative forces such as friction forces. One may note that, the flexible links of the flexible manipulator are modeled as clamped base beams, see [START_REF] Bian | Impact vibration attenuation for a flexible robotic manipulator through transfer and dissipation of energy[END_REF]. The dynamic model of the flexible manipulator is supposed to be estimated by an observer or practically identified, see for instance [START_REF] Sasaki | Vibration suppression control of redundant manipulator with flexible structure by considering nullspace motion[END_REF].

The dynamics of both rigid and flexural motion may be described separately from (1) as follows M rr qr + Γr = τ

(2)

M f r qr + M f f q f + D f q f + K f q f = Γ f (3)
where M rr ∈ IR n r ×n r , M f r ∈ IR n f ×n r and M f f ∈ IR n f ×n f are block matrices of the global mass matrix,τ ∈ IR n r is the control torques vector applied to the manipulator joints, Γr ∈ IR n r is the non-linear dynamics in the rigid motion. Γ f ∈ IR n f is nonlinear force which is applied to the flexural dynamics. D f ∈ IR n f ×n f and K f ∈ IR n f ×n f are respectively the damping and the stiffness matrix of the flexural motion.

Redundancy resolution

The end-effector position of the flexible manipulator, which is operating in n t dimensional work space, is function of both joints angles and elastic displacement coordinates as follows

x = f (q r , q f ) (4)
where x ∈ IR n t .

The kinematic model can be derived from ( 4) by time differentiation ẋ = J r qr

+ J f q f (5) ẍ = J r qr + Jr qr + J f q f + Jf q f (6)
where J r ∈ IR n t ×n r and J f ∈ IR n t ×n f are the rigid and the flexible Jacobian matrix respectively.

From (3), one may notice that the flexural displacement can not be controlled directly by the joints torques. Usually, in the motion planning process only the nominal joints motion is considered to generate the desired trajectory for the flexible robot manipulators, see [START_REF] Sasaki | Vibration suppression control of redundant manipulator with flexible structure by considering nullspace motion[END_REF] and [START_REF] Gao | Research on Joint Torque Optimization Method of Redundant Space Manipulators with Vibration Suppression[END_REF]. In this case, the term (J f q f + Jf q f ) is not considered in redundancy resolution. The dynamical resolution of redundancy which is presented in (6) becomes ẍ = J r qr + Jr qr (7)

The flexible manipulator is kinematically redundant, i.e. the number of actuated (rigid) joints n r is greater than the dimension of the operational work-space n t . The degree of redundancy m is defined as m = n rn t .

The resolution of redundancy based on the pseudo-inverse of the Jacobian matrix is usually used to determinate the joint acceleration. In [START_REF] Gao | Research on Joint Torque Optimization Method of Redundant Space Manipulators with Vibration Suppression[END_REF], the projected joints dynamic in the Jacobian null-space is only devoted to damp out the vibration. In this case, having extra degrees of redundancy will be used to damp out effectively the vibration but not to fulfil additional tasks. In this study, we investigate the use of the extra number of redundancy DOFs in aim to perform more than one second task simultaneously. Therefore, an other method for resolution of redundancy, which is based on the manipulator dynamic analysis, has been considered in the development of the proposed strategy to exploit all the available degrees of redundancy, see [START_REF] Ghosal | Dynamical Resolution of Redundancy for Robot Manipulators[END_REF]. The joint torques are used in redundancy resolution, which is advantageous since they are the main input variables to the manipulator controller.

Substituting (2) into ( 7), the end-effector acceleration can be rewritten as ẍ = J r M -1 rr (τ -Γr ) + Jr qr (8)

The following terms are defined for the convenience of the further study development,

ϑ = J r M -1 rr (9) χ = Jr qr -J r M -1 rr Γr (10)
The equation ( 8) then becomes ẍ = ϑ τ t τ t + ϑ τ r τ r + χ (11) where the joint torques are decomposed to two sets, τ t ∈ IR n t the task torques to perform a given task, and τ r ∈ IR m the redundant torques which are extra torques regarding the task. ϑ τ t ∈ IR n t ×n t and ϑ τ r ∈ IR n t ×m are matrices that form the matrix ϑ , where the ith task torque τ t,i (respectively the jth redundant torque τ r, j ) is associated with the ith column of the matrix ϑ τ t (respectively jth column of the matrix ϑ τ r ), (i = 1, 2, ..., n t ) and ( j = 1, 2, ..., m). The term χ ∈ IR n t aggregates all non-linear dynamics.

At time instant t * , the joint coordinates ( qr (t * ), q r (t * )) are generally known, as a result the vector χ can be determined, as a consequence, the mapping h given by h : τ -→ ẍ ( 12) is affine and the admissible acceleration in this case are ẍac ∈ D τ , where

D τ = { ẍ = h(τ), s.c τ min,i ≤ τ i ≤ τ max,i , i = 1, 2, .
., n r } where τ min,i and τ max,i are, respectively, the lower and the upper limit of the ith joint torque.

The analytical redundancy resolution at acceleration-level is used to set up the proposed approach. The algorithm combines both an alternative control design to feedback linearization and a method for reducing the vibrations, in aim to improve the endeffector accuracy.

CONTROL DESIGN VIA REDUNDANCY RESOLUTION

The feedback linearization control approach is the combination of two imbricated control loops. First, the inner-loop is constructed to ensure the system linearization. Then, the outerloop is designed to guarantee the control of the resulted linear system. This control strategy is one of the robust schema for the manipulator control. The free m redundant joint torques τ r,i , (i = 1, 2, ...m) in the mapping (11) are used to cancel the nonlinear dynamics, unlike the traditional feedback linearization control which exploits all the manipulator torques for the same process. Analytically, that can be reached by setting

ϑ τ r τ r + χ = 0 (13)
Since the vector χ is in IR n t , only n t torques from m available redundant torques are sufficient to cancel completely the nonlinear term χ. But if the redundant DOFs number is less than the task dimension, i.e (m < n t ), we can cancel χ in only m preferential directions. Furthermore, in the cases where there are excess of redundant degrees of freedom, i.e (n t < m), an infinity of the redundant torques combination can cancel completely the nonlinear term χ. In this paper, the last case is deeply analyzed. The forward analysis have demonstrated that the additional (m-n t ) redundant degrees are useful to affect the flexural motion and eventually damp out the induced vibrations.

Considering the last case, the general solution for redundant torques that satisfies the cancellation of the non-linear dynamics χ is given by

τ r = -ϑ + τ r χ + (I m -ϑ + τ r ϑ τ r )τ ε ( 14
)
where ϑ + τ r ∈ IR m×n t is the generalized inverse matrix of ϑ τ r , the resulting vector (I mϑ + τ r ϑ τ r )τ ε ∈ N(ϑ τ r ), where N(ϑ τ r ) is the null-space of ϑ τ r , τ ε ∈ IR m is an arbitrary vector, I m ∈ IR m×m is the unit matrix.

After the linearization process has been achieved, the joint torques τ t and τ r , which guarantee the system linearization, are determined by solving the equations in ( 11) and ( 14) for a given end-effector acceleration ẍd , providing that the matrix ϑ τ t is not singular. One may note that, the choice of torque sets τ t and τ r depends on various factors, such as torques magnitude limits, for more details see [START_REF] Ghosal | Dynamical Resolution of Redundancy for Robot Manipulators[END_REF]. After that, a conventional linear controller can be used to guarantee the end-effector tracking. In the following analysis, we investigate the possibility of the additional self-motion to reduce the deformations due to the flexural motion. The diagram of the proposed control method is presented in the Fig. 1.

METHOD FOR VIBRATIONS REDUCTION

Analysis of the flexural motion

The presence of the arbitrary vector τ ε in ( 14) may generate an infinite number of solutions for the redundant torques τ r . This self-motion ability, which doesn't influence the end-effector motion, is used to adjust the joint motions and the flexural motion.

Substituting the joints acceleration qr value from (2) into (3) yields to

M f f q f + D f q f + K f q f = F f (15) F f = Γ f + M f r M -1 rr Γr -µ τ t τ t -µ τ r τ r (16
) where F f ∈ IR n f is the exciting force of the flexural motion, µ τ t ∈ IR n f ×n t and µ τ r ∈ IR n f ×m are blocks matrices of the matrix µ ∈ IR n f ×n r , where µ = M f r M -1 rr . The extended expression of the exciting force F f , that includes the arbitrary vector τ ε , is obtained by combining the two equations ( 14) and ( 16 Where τ(τ t , τ r ) is a torque vector which is composed by components combination of τ t and τ r . τ c is the control torque for the linearized system.

F f = -µ τ r (I -ϑ + τ r ϑ τ r )τ ε + Γ f + M f r M -1 rr Γr -µ τ t τ t + µ τ r ϑ + τ r χ (17)
Let's consider that the damping matrix D f is diagonalizable, then the modal analysis of (15) will provide n f independent equations as follows

φ f + D φ φ f + K φ φ f = F M ( 18
)
where

D φ = diag(2ξ i ω i ) ∈ IR n f ×n f , K φ = diag(ω 2 i ) ∈ IR n f
×n f , ξ i and ω i are the ith modal damping term and the ith natural frequency for (i = 1, 2, ...n f ). φ f = M -1 m q f ∈ IR n f is a vector that contains the mode coordinates, where M m ∈ IR n f ×n f is the mode shape matrix.

F M = M -1 φ M T m F f , F M ∈ IR n f is the modal generalized force, where M φ ∈ IR n f ×n f is the modal mass matrix, M φ = M T m M f f M m .
The arbitrary vector τ ε is considered as a control optimization term that can be used to damp out the flexural motion. The modal generalized force is a controllable force via self-motion and it contains two terms: nonlinear term F M,n and control term F M,c which depends on the arbitrary vector τ ε

F M,n = M -1 φ M T m • [ Γ f + M f r M -1 rr Γr -µ τ t τ t + µ τ r ϑ + τ r χ ] (19) F M,c = -M -1 φ M T m µ τ r (I -ϑ + τ r ϑ τ r )τ ε (20)

Flexural motion damping

To damp out the vibrations of oscillatory system, the exciting force should be eliminated, and more than this, the system damping should be increased as well. Since the self-motion can directly affect the vibration responses of the flexible redundant manipulator as shown in (20), the arbitrary vector τ ε could be chosen properly to increase the damping of the manipulator vibration modes and suppress the exciting force, while the desired end-effector position is tracked. The value of the control term F M,c in the modal generalized force which satisfies the vibrations damping is given by

F M,c = -D d φ φ f + D φ φ f -F M,n (21
) substituting this chosen term into the equation ( 18), the flexural motion system becomes φ

f + D d φ φ f + K φ φ f = 0 n f (22)
where, D d φ = diag(2ξ d i ω i ) ∈ IR n f ×n f is the desired modal damping matrix, ξ i is the ith desired modal damping, (i = 1, 2, ...n f ). In practice, the damping ratio D φ is negligible, a numerical estimator using torque/force sensors may be used to obtain the flexural reactive force. Furthermore, the numerical differential value of φ f , φ f , may be determined as well, see [START_REF] Hirano | Simultaneous control for endpoint motion and vibration suppression of a space robot based on simple dynamic model[END_REF].

The value of the arbitrary vector τ ε which satisfies the vibration reduction is the solution of the following equation and that by substituting ( 19) and ( 20) into ( 21)

Aτ ε = B (23) where, A = -M -1 φ M T m µ τ r (I -ϑ + τ r ϑ τ r ) (24) B = -M -1 φ M T m • [ Γ f + M f r M -1 rr Γr -µ τ t τ t + µ τ r ϑ + τ r χ ] + (D φ -D d φ ) φ f (25)

Analysis of self-motion solution

In the relation ( 23), there are m unknown components of τ ε with n f equations. The solution of this system of equations depends on the rank (ρ) of the matrix A and the augmented matrix [A B]. When the manipulator is not in singular configuration, ρ((Iϑ + τ r ϑ τ r )) = mn t . In general, the matrix µ τ r has the full rank, ρ(µ τ r ) = min(n f , m). Since the matrix M T m has the full rank, ρ(M T m µ τ r ) = ρ(µ τ r ) = min(n f , m). The rank of matrix A depends on the number of both the flexural degrees n f and the difference between the redundant degrees and the task dimension mn t . The first case considers (mn t ) ≥ n f , the redundant degrees are more than the flexural DOFs and ρ([A B]) = n f , there are infinity of solution τ ε that can damp out the vibration effectively. The second case considers (m -

n t ) < n f , either m < n f or m > n f , we have ρ(A) = m - n t , if ρ([A B]) > ρ(A)
there are no solutions to the system of equations. But, on the other hand, there are an infinity of solutions that minimize

||Aτ ε -B|| 2 when ρ([A B]) = ρ(A), the general solution τ * ε is given by τ * ε = A + B + (I -A + A)τ ε,0 (26)
where τ ε,0 ∈ IR m is arbitrary vector, A + ∈ IR m×n f is the generalized inverse matrix of A. The self-motion solution τ * ε can damp out effectively the flexural motion. The equation ( 26) shows that the solution of τ * ε still depends on another arbitrary vector τ ε,0 , which means that there are many solutions that satisfy the vibration reduction. Thus, the arbitrary vector τ ε,0 may be exploited to design secondary task for ensuring additional performance criteria. In this case, generally, the additional ability can be used to avoid singularities or obstacles, on promise of accuracy improvement for the flexible redundant manipulator.

SIMULATIONS AND DISCUSSION

To test and demonstrate the effectiveness of control strategy presented in this paper, 5R planar flexible manipulator is used in numerical simulations, see Fig. 2. The first 4 links are supposed to be rigid and the last link is flexible. Since the planar manipulator has five rigid DOFs n r = 5 and n t = 2 (since the end-effector rotation is considered to be free), three degrees of redundancy are available, m = 3. The motion of the end-effector is expressed in the origin reference frame (XY). In practice, the first mode has the most contribution to the flexural motion than the other vibration modes. In the following analysis, only the two first vibration modes of the flexural motion are retained, n f = 2. The parameters of 5R planar manipulator are given as follows: the length of each 5 links is 1m. The rigid links have X Y r,1 q r,2 q r,3 q r,4 q r,5 q f Fig. 2. Redundant 5R planar flexible manipulator. square cross-section of 0.00045m 2 and density of 8000Kg/m 3 . The flexible link has square cross-section of 0.000225m 2 with elastic modulus of 69GPa and density of 2700Kg/m 3 . All terms due to gravity in the dynamical equations are null as long as the manipulator moves in a horizontal plane (XY). Without considering any external forces, only Coriolis force is retained in the non-linear term, Γr = C( Θ, Θ). In this example, we set the two first actuators as task torques τ t,1 = τ 1 and τ t,2 = τ 2 .

The last 3 torques are considered to be the redundant torques ( τ r,1 = τ 3 , τ r,2 = τ 4 , τ r,3 = τ 5 ). In addition, no limits on the joint torques have been taken into account. The initial joint angles are q r (0) = [80 -145 120 -60 -30] T degrees, with null initial velocities qr (0) = [0 0 0 0 0] T rad/s. The desired velocity of the end-effector is set to be null at the beginning and the end of the motion. The end-effector was controlled to follow a desired accelerations which has Bang-bang profile as follows

ẍd (t) =    [0.1 0.45] T m/s 2 0s ≤ t ≤ 1.5s [0 0] T m/s 2 1.5s < t ≤ 3s [-0.1 -0.45] T m/s 2 3s < t ≤ 4.5s (27)
The results have shown the responses of the system to the proposed control approach for the two cases, i.e. with and without the utilization of torques self-motion to damp out the induced vibration.

Case 1 (τ ε is null) : In this case, the self-motion is only calculated to cancel the system non-linearity in (13). It was possible to apply the proposed control algorithm because the degree of the redundancy is greater than the work-space dimension. Besides, a conventional PD controller at joint level is designed and then applied to the resulting linear system for maintaining the end-effector position tracking. On the other hand, it is shown in Fig. 3 that the end-effector deformations are non-negligible, since the end-effector vibrations have not been reduced. The end-effector follows the desired trajectory with vibrations of maximum amplitude equals to 6.6mm.

Case 2 (τ ε is not null) : In this case, τ ε is calculated from (26). Taking into account the relation ( 26), the self-motion is recalculated to get rid of the non-linearity, and to damp out the vibration modes of the resulting flexural motion simultaneously. One may note that, the torques magnitude of τ r may become very large as the end-effector approaches singularity.

As shown in Fig. 3, the deformations of the end-effector are well reduced via self-motion. Compared the Fig. 4 with Fig. 5, it can be clearly seen that the self-motion has been changed to take into account the additional task of reducing the vibration deformations, see Fig. 6. The end-effector follows the desired trajectory with small deformations (around 0.02mm of maximum amplitude).

The torque profiles are shown in Fig. 4 and Fig. 5. The actuator torques of the two first actuators (τ 1 and τ 2 ) are significantly larger than the other actuator torques. This depends on which actuator torques are chosen as τ t or τ r . Therefore, it is preferably to have a large limit for the first two torques to develop a sufficient energy during the manipulator motion. The resulted torques in actuators τ 3 , τ 4 and τ 5 in Fig. 4 are slightly different to those in Fig. 5. This is due to the changing in the self-motion solution in purpose to damp out the vibrations.

The real time implementation of the proposed algorithm is feasible. Since the robot state such as joint coordinates are generally measurable via implemented sensors. The main challenge is to measure the flexural motion states. Nevertheless, that may be resolved using information derived from torque/force sensors to estimate both the flexural reactive force and the modal coordinates following the estimation approach in Hirano (2014). Thus, the discretization and the synchronization problem of data is avoided. Afterwards, the robot model can be simultaneously linearized since the calculation can be factorized.

In the all previous methods presented in [START_REF] Gao | Research on Joint Torque Optimization Method of Redundant Space Manipulators with Vibration Suppression[END_REF]; [START_REF] Nguyen | Dynamic control of flexible, kinematically redundant robot manipulators[END_REF]; [START_REF] Bian | Study on vibration reduction and mobility improvement for the flexible manipulator via redundancy resolution[END_REF]; [START_REF] Hirano | Simultaneous control for endpoint motion and vibration suppression of a space robot based on simple dynamic model[END_REF]; [START_REF] Sasaki | Vibration suppression control of redundant manipulator with flexible structure by considering nullspace motion[END_REF], the redundancy resolution is used for vibration reduction only and none see the advantage of having more redundant degrees specially in case of the hyper-redundant flexible manipulator. Although, we proposed algorithm that takes advantage of all redundant DOFs, first for controlling the manipulator (feedback linearisation via redundancy resolution), and secondly to suppress the induced vibrations of the flexible links. Moreover, no model simplifications are needed for control design which ensures the reliability of the proposed algorithm. The vibration behaviour is well treated compared to the previous articles, through the analysis, one can see clearly how the vibration are damped. The simulation result shows the effectiveness of the proposed approach.

CONCLUSION

In this paper, the self-motions of redundant flexible manipulators are exploited for both designing a control algorithm and establishing a novel method for vibrations reduction of its endeffector while tracking a given trajectory. The effectiveness of this approach has been validated by the results of numerical simulations. It has been shown that the self-motions have cancelled the undesirable dynamics, besides, through the dynamic coupling between the rigid and the flexural motions the vibrations of the end-effector have been reduced as well. The analysis of the self-motions satisfying the vibration reduction revealed an additional ability which may be used on promise of improving the kinematic and/or dynamic performances of flexible redundant manipulators.

Fig. 1 .

 1 Fig. 1. Diagram of the proposed feedback control (blue zone)with vibration reduction action (gray zone) for flexible redundant manipulator (red zone). Where τ(τ t , τ r ) is a torque vector which is composed by components combination of τ t and τ r . τ c is the control torque for the linearized system.

Fig. 4 .

 4 Fig. 3. End-effector deformations for both cases.

Fig. 6 .

 6 Fig. 5. Control torque profiles-case 2.

ACKNOWLEDGEMENTS

The author would like to thank the "European Union" and the french region "Hauts-de-France" for the financial support.


Mohamed Amine Rahmouni
email: mohamedamine.rahmouni@cea.fr

Richard Bearee
email: richard.bearee@ensam.eu

Adel Olabi
email: adel.olabi@ensam.eu

Mathieu Grossard
email: mathieu.grossard@cea.fr

Eric Lucet
email: eric.lucet@cea.fr

Robot trajectory generation for three-dimensional flexible load transfer

Keywords: Vibration Suppression, Elastic Deformation, Robot Manipulators, Trajectory Shaping

This paper addresses the problem of reducing the elastic deformations and the residual vibrations of flexible loads when they are handled by a robot manipulator. During the manipulation of the low-stiffness load, such as bumper or exhaust system in automotive industry, large motion-induced deformations and vibrations may be induced. These deformations will have detrimental effects on the settling time, on the accuracy and on the integrity of the operational process in a constraint environment. The trajectory shaping approaches, i.e smoothing filter or input shaping method are well-known solutions for the suppression of the residual vibration at the end of a rest-to-rest motion. However, using trajectory shaping technique alone may not be sufficient to suppress the static elastic deformations during the transfer phase of the object. Thus, the main contribution of this paper is to propose a two stages feedforward based approach that combines trajectory shaping technique for vibration reduction, with a deformation compensation trajectory. The latter exploits the rotation space of the manipulator to mitigate the flexural motion of the flexible load. The effectiveness of the proposed method is demonstrated by experimental validations on an industrial robot Kuka iiwa.

I. INTRODUCTION

The robotic manipulation is highly concerned with several industrial tasks including assembly operations, loading and unloading parts from an assembly line or a machine and the safe manipulation more specifically in physical Human-Robot Interaction. When handling flexible load at the robot tip, the previous robotics tasks accuracy, sometimes even their feasibility, may be significantly limited by undesired motioninduced vibrations and deformations of the manipulated object.

Many works on flexible objects manipulation have been done in different contexts, such as the clothes manipulation [START_REF] Pan | Motion planning for fluid manipulation using simplied dynamics[END_REF]- [START_REF] Li | Model-driven feedforward prediction for manipulation of deformable objects[END_REF], the ropes manipulation [START_REF] Huang | Leveraging appearance priors in non-rigid registration, with application to manipulation of deformable objects[END_REF] or the transfer of liquid containers with minimal sloshing [START_REF] Aribowo | Input shaping control to suppress sloshing on liquid container transfer using multi-joint robot arm[END_REF], [START_REF] Aribowo | Integrated trajectory planning and sloshing suppression for three-dimensional motion of liquid container transfer robot arm[END_REF]. The existing methods can be coarsely divided into the state feedback approach and the open-loop approach.

The state feedback approach is generally based on the use of extra force or torque sensors to measure the vibration of the flexible object. In [START_REF] Jiang | Robotic manipulation of exible objects based on vibration control using force sensors[END_REF], a vibration control scheme is synthesized based on linear feedback of vibration measured by the force sensor. The control is established by adding the vibration control scheme into the controller of the robot wrist joint. In [START_REF] Tanner | Mobile manipulation of flexible objects under deformation constraints[END_REF], authors propose a velocity field trackingcontrol scheme using proportional integral force feedback to transport and manipulate deformable material with a mobile manipulator. The efficiency of the state feedback approach may be primarily limited by the technical and economic viability of adding extra sensors.

On the open-loop approach side, input shaping technique is a well-established method to suppress motion-induced residual vibrations [START_REF] Singer | Preshaping Command Inputs to Reduce System Vibration[END_REF], [START_REF] Singhose | Command shaping for flexible systems: A review of the first 50 years[END_REF]. A set of finite impulses, called shaper, are convolved with the desired trajectory to suppress or reduce the residual vibration at the end of point-to-point motion. This technique turned to be a robust and effective algorithm in practical applications, where only the natural frequency and the damping ratio of the contributed flexible modes are required [START_REF] Singhose | A comparison of input shaping and timeoptimal flexible-body control[END_REF]. Some of the vibration suppression properties of input shaping technique can be directly found in some trajectory patterns. For instance, in [START_REF] Biagiotti | FIR filters for online trajectory planning with time-and frequency-domain specifications[END_REF], authors propose a finite impulse filter-based trajectory generation taking the vibration reduction into account. Such approach is specifically applied to the jerk-limited trajectory generation in [START_REF] Bearee | Dissociated jerk-limited trajectory applied to timevarying vibration reduction[END_REF], [START_REF] Bearee | New Damped-Jerk trajectory for vibration reduction[END_REF]. The jerk time (i.e. the acceleration slope time) can be explicitly linked to the vibration caused by a dominating flexible mode. Hence, the jerk-limited profile can be seen as a special combination of the smoothness property of polynomial profiles and the vibration reduction property of basic input shaping. The main limitation of the open-loop approach is that the vibration can be suppressed at the end of the motion, but the deformation of the flexible object during the motion still exist.

In [START_REF] Moriello | Manipulating liquids with robots: A sloshing-free solution[END_REF]- [START_REF] Feddema | Control for slosh-free motion of an open container[END_REF], the oscillatory behavior of a liquid container during motion is interpreted as a spherical pendulum. Shaping filter are used to suppress the oscillatory behavior of the liquid and the trajectory of the robot is adapted by introducing a tilting compensation value to compensate for the deformation during the motion.

Based on the previous combined open-loop approaches, this paper proposes a trajectory generation method for industrial robot manipulator, so that the deformations during threedimensional motion and the residual vibrations of the flexible load are suppressed. The sensitivity analysis of the proposed method regarding the estimation of the dominating modal frequencies is presented. This promising feedforward approach has been experimentally tested on an industrial Kuka LBR iiwa.

II. LOAD DEFORMATION AND VIBRATION COMPENSATION

A. Modeling of the load deformations

In this work, the translational displacement task of a flexible object with a n DOFs serial manipulator (n ≥ 6) is investigated. The robot manipulator is kinematically redundant for the previous 3d task and among the remaining extra DOFs, the orientation in 3d space will be used in the proposed methodology for the elastic deformation compensation. The flexible object to be manipulated is assumed to be submitted to dominating modal deformations. In practice the flexible motion of a nearlinear object, i.e. with a main axis, is dominated by the first bending natural frequency in the two transverse directions. For the sake of simplicity, Fig. 1 shows the planar definition of the problem for one bending deformation. One may note that the flexible motion in the other plane can be handled using the same approach. The flexible load is rigidly grasped by the robot gripper at the robot Tool Center Point (TCP). Its position is noted P tcp .

In the following, the elastic deformation of the object is considered. The pure flexural motion of the object tip, noted w(t), is assumed to occur in Y tcp direction only, i.e. no axial deformation is considered. Hence, the coordinates of the object tip, noted P ob , can be expressed in the TCP frame R tcp with the origin point P tcp as follows

P ob = l w R tcp (1)
with l and w respectively the length and the flexural deformation of the object. The dynamics of the points P tcp and P ob in the fixed absolute frame are submitted to the same rigid motion imposed by the robot TCP trajectory. The evolution of the position P ob results from the sum of the previous rigid mode with the elastic deformation w(t). Noting y tcp (t) and y ob (t) respectively the position of the points P tcp and P ob in Y tcp direction, the deformation w(t) can be written

w(t) = y ob (t) -y tcp (t) (2)
Considering now small variation of the deflection angle, noted θ in Fig1, the approximated deflection angle of the object can be expressed as

θ (t) = w(t) l (3)
The dynamic vibration of the object tip, which occurs in Y tcp direction only, can be modeled by the following general second order equation in continuous Laplace domain

y ob (s) y tcp (s) = ω 2 ω 2 + 2 • ξ • ω • s + s 2 (4)
with s is the Laplace operator, ω and ζ are respectively the natural frequency and the damping ratio of the considered bending mode. Substituting (2) and ( 3) into (4), the dynamics of the deflection angle θ can be written

θ (s) = - 1 l • 2 • ξ • ω • s + s 2 ω 2 + 2 • ξ • ω • s + s 2 • y tcp (s) (5)
which leads in time domain to the differential equation

θ (t) + 2 • ξ ω • θ (t) + 1 ω 2 • θ (t) = - 2 • ξ l • ω • ẏtcp (t) - 1 l • ω 2 • ÿtcp (t) (6)
Hence, the deflection angle for the considered bending mode will be driven by the acceleration and the velocity of the robot TCP in (Y tcp ) direction. In the same way, the TCP acceleration in Z ob direction may induce vibrations of the first bending mode around Y tcp . Assuming that the load vibration can be canceled or reduced, i.e. that the oscillatory behavior given by (6) vanishes and θ = θ = 0, the remaining deflection angle would be dominated by the static contribution of (6) given by

θ (t) = - 2 • ξ l • ω • ẏtcp (t) - 1 l • ω 2 • ÿtcp (t) (7)
Equation [START_REF] Jiang | Robotic manipulation of exible objects based on vibration control using force sensors[END_REF], which establishes the relation between the static angular deformation during the motion and the instantaneous acceleration and velocity of the robot, is crucial for the proposed methodology.

B. Methodology for vibration and deformation compensation

The proposed methodology, depicted in Figure 2, is a two stages feedforward-based approach. The first stage consists in compensating for the motion-induced vibration via trajectory shaping. The second stage aims at suppressing the remaining static component of the tracking error during the motion, i.e. the deformation of the load tip, using the rotational DOFs of the robot TCP for a proper compensation.

1) First stage -Vibration reduction via trajectory shaping: The dynamic model presented in ( 6) is a second order system and its excitation force is a combination of the TCP kinematics (acceleration and velocity) along (Y ob ) axis. Trajectory shaping for vibration compensation of such an oscillatory behavior is a well established method. The vibration reduction may be obtained by several types of filters. In this study, the jerk-limited profile, which is available in most robot controllers, is used for the joints trajectory. Fig. 3 describes an acceleration profile for a rest-to-rest motion using Jerk-limited pattern compared to the motion with simple acceleration-limited pattern. The extra motion duration T J compared to the acceleration-limited solution is equal to the jerk time. As depicted in Fig. 4, a Jerklimited trajectory can be easily and efficiently synthesized by applying a moving averaging Finite Impulse Response (FIR) filter, with a filter time equal to T J , to an acceleration limited profile. The FIR filter F JL can be described in continuous time domain by the function (8), with s is the Laplace domain operator and A i are the impulse values for the equivalent shaper.

F JL = 1 s (A 1 + A 2 e -s T J ) (8)
As described in [START_REF] Bearee | New Damped-Jerk trajectory for vibration reduction[END_REF], the filter F JL has some common properties with classical shapers. Considering an undamped vibration with a given natural frequency ω n , by setting the filter time T J equals to the natural period time of the vibration, the residual vibration are suppressed. This simple tuning methodology of the FIR filter F JL can be expressed as follows

T J = k 2π ω n A 1 = -A 2 = 1 T J ( 9
)
where k is a positive integer, classically set to one for time saving. A detailed analysis and the adaptation of such a method for the damped vibrations case can be found in [START_REF] Bearee | New Damped-Jerk trajectory for vibration reduction[END_REF], [START_REF] Biagiotti | Optimal Trajectories for Vibration Reduction Based on Exponential Filters[END_REF]. Hence, using the previous feedforward approach, the residual vibrations and most of the vibrations occurring during the motion are suppressed. The next stage consists in suppressing the deformation of the load tip during the motion.

2) Second stage -Residual deformation compensation: Assuming that the load vibration is suppressed by the trajectory shaping method, the residual deflection angle θ will be given by [START_REF] Jiang | Robotic manipulation of exible objects based on vibration control using force sensors[END_REF]. This remaining motion-induced deformation during the motion can be mitigated by a proper compensation scheme of the deflection angle using the rotational space of the robot TCP, i.e. by adjusting the orientation of the load tip during the motion. This action is added to the design of the feedforward action and can be written

Rot new tcp (θ (t)) = Rot tcp • Rot comp (-θ (t)) (10)
with Rot tcp and Rot comp are, respectively, the predefined TCP rotation matrix in base frame and the added rotation matrix which counteracts for the estimated deflection angle. From [START_REF] Jiang | Robotic manipulation of exible objects based on vibration control using force sensors[END_REF] it is noticed that the angle θ depends eventually on the TCP acceleration and, to a lesser extent, on the TCP velocity (if the damping is sensitive). Thus, the vibration reduction properties introduced into the TCP trajectory via the shaping method will be extended to the added TCP rotational motion, which will not introduce new vibrations to the load.

C. Sensitivity analysis

The proposed feedforward method is inherently sensitive to parameters estimation and variations. Naturally, the main parameter for the vibration and deformation compensation is the value of each modal frequency to be canceled. In input shaping framework, it is customary to define the sensitivity plot, i.e. the percent of residual vibration amplitudes as function of the normalised variation (ω/ω n ), with ω and ω n respectively the estimated and the real value of the natural frequency. Such a sensitivity plot is depicted in Fig. 5 for the jerk-limited profile, used in this study, and for two wellknown shapers, the Zero-Vibration (ZV) shaper and the Zero-Vibration-Derivative (ZVD) shaper (see [START_REF] Singhose | Command shaping for flexible systems: A review of the first 50 years[END_REF]). Here, 100% residual vibration corresponds to the vibration induced by an acceleration step command, and robustness or sensitivity is classically defined as the frequency range over which a command induces less than 5% vibration. These sensitivity curves reveal how much residual vibration will exist when there is an error in the estimation of the system frequency. Hence, the jerk-limited trajectory tune to cancel the residual vibration of one flexible mode will be around 40% more robust than a ZV shaper. Moreover, the low-pass filtering property which is introduced by the integrator is effective to reduce significantly the impact of the neglected high frequency modes without any required frequency identification, which is practical for the industrial applications.

The same methodology can be followed to express the sensitivity plot for the deformation compensation. The deformation is affected by the variation of the angle compensation θ given by [START_REF] Jiang | Robotic manipulation of exible objects based on vibration control using force sensors[END_REF] and can be expressed as

w n = lθ = - 2ξ n ω n ẏtcp - 1 ω 2 n ÿtcp (11)
Where ω n and ξ n are respectively the frequency and the damping values of the dominant mode in the load vibration. Considering a compensation action with the estimated frequency and damping value ω and ξ , the remaining residual deformation can be expressed as

∆w w n = w n -w w n = 1 - ω n ω 2 • a • ( ω ω n ) + b a + b (12)
Where

a = -2 • ξ • ẏob , b = -ÿob /ω n .
Assuming now that the damping can be neglected for the dominating deformations, the sensitivity of the deformation compensation is given by

∆w w n = 1 - ω n ω 2 (13)
Fig. 6 represents the sensitivity of the deformation compensation according to the normalized frequency. 100% residual deformation corresponds to the initial deformation without compensation. One can note that the sensitivity to keep the residual deformation below 5% is almost half the sensitivity of the vibration compensation action. However, in a more practical manner, an error of 10% for the estimated vibration frequency still allow to reduce the deformation by 80%.

III. EXPERIMENTAL VALIDATION A. Experimental setup

To validate the proposed methodology and evaluate its effectiveness, an experimental validation has been done with a 7 DOFs robot manipulator kuka iiwa 14, which manipulates a flexible beam (see Fig. 7). A laser tracker from API is used to measure the beam tip displacement along a predefined linear trajectory of the robot TCP. Only the first vibration mode of the flexural motion is considered in the motion direction. The system parameters have been identified as follows: the natural frequency is ω = 2.2Hz (0.038rad/s); the beam length is l = 0.25m ; the value of the damping ratio is close to ξ = 0.001 and therefore the damping can be reasonably neglected for the considered beam.

Four different cases are presented in the following. The first case corresponds to the initial response of the beam with no vibration and deformation compensation. The second experiment consists in adding the deformation compensation action (static compensation) only. In the third case, only the vibration compensation is used and the last case combines the vibration and deformation compensation actions. The different robot TCP trajectories have been generated off-line and then send to the robot joint controller with a sampling time equal to 4ms.

B. Discussions

Fig. 8 shows the evolution of the measured load tip position and Fig. 9 shows the measured position tracking error of the load tip for the different trajectories. Table I gives the main quantitative results for the different trajectories. For the initial measured response and for the case with the deformation compensation (static compensation) only, the reference trajectory is an acceleration-limited trajectory. For the two other cases, the reference trajectory is the jerk-limited trajectory resulting from the FIR filtering of the acceleration-limited trajectory.

Case 1 -Initial response: The initial response of the load tip is clearly characterized by two different phases. During the motion, the load tip exhibits significant vibrations superimposed to the deformation error pattern. After the motion, the undamped residual vibrations are clearly visible.

Case 2 -Deformation (static) compensation only: In this case, only the angle compensation technique is used. As described in Fig. 9, the vibration amplitude is increased due to the fact that the compensation is done on the unshaped acceleration profile. Hence, vibrations induces by the translation motion of the TCP are naturally not suppressed (see [START_REF] Aribowo | Integrated trajectory planning and sloshing suppression for three-dimensional motion of liquid container transfer robot arm[END_REF]), but worse, the added rotation trajectory of the TCP induces new vibrations. For the test trajectory, the contribution of each motion on the vibration level is clearly additive. This case reflects the fact that the proposed deformation compensation cannot be used without a vibration suppression action.

Case 3 -Vibration compensation only: The jerk-limited profile is used in this case to supress the vibrations. Fig. 9 shows that the vibrations are efficiently damped out while the deformation still exists during the motion. The residual vibration amplitude is less than 5% of the initial one. The remaining small vibration amplitude can be related to different factors such as the estimation bias of the natural frequency and the neglect of the damping effect.

Case 4 -Vibration and deformation compensation: This case shows that the proposed two-stages feedforward method is efficient to significantly reduce the vibrations and the deformation. The residual tracking error of the load tip is less than 18% and the residual vibration value is equivalent to the previous case.

IV. CONCLUSION

In this study, we proposed a feedforward strategy to compensate for the deformation and vibrations during the translational motion of a flexible load with a robot manipulator. The method is based on the combination of a trajectory shaping stage for the vibration suppression, with a trajectory compensation scheme for the elastic deformation occurring during the motion. The latter action makes use of the rotation space of the robot manipulator, which is not used for the task. The sensitivity of the two actions with regard to the estimated vibration frequencies has been derived. The experimental tests, which has been conduct using a Kuka iiwa robot and a pure flexible beam, demonstrated the effectiveness of the proposed approach. Future works will focus on the generalization of the methodology for translational and rotational motion of the flexible load.

Fig. 1 .

 1 Fig. 1. Parameters definition: robotic manipulation of a cantilevered object with one dominating bending motion.

Fig. 2 .Fig. 3 .

 23 Fig. 2. Deformation w: (a) initial deformation and vibration. (b) deformation with shaping trajectory for vibration compensation. (c) residual deformation with shaping trajectory and deformation compensation.

Fig. 4 .

 4 Fig. 4. FIR filtering scheme equivalent to jerk-limited profile. T J is the filter time (equivalent to the jerk time) and T e is the sampling time of the controller.

Fig. 5 .Fig. 6 .

 56 Fig.5. Sensitivity of the residual vibration suppression action[START_REF] Bearee | New Damped-Jerk trajectory for vibration reduction[END_REF].

Fig. 7 .

 7 Fig. 7. The experimental setup.

Fig. 8 .Fig. 9 .

 89 Fig. 8. The load tip position for different cases. The after-motion state of the no-shaped (acceleration-limited) and shaped (jerk-limited) trajectory cases is depicted by the blue and the red zones respectively.

ACKNOWLEDGMENT

The author would like to thank the "European Union" and the french region "Hauts de France" for the financial support.


José M Mendes Filho

Eric Lucet
email: eric.lucet@cea.fr

MULTI-ROBOT MOTION PLANNING: A MODIFIED RECEDING HORIZON APPROACH FOR REACHING GOAL STATES

Keywords: multi-robot motion planning, nonholonomic mobile robot, distributed planning, receding horizon

This paper proposes the real-time implementation of an algorithm for collision-free motion planning based on a receding horizon approach, for the navigation of a team of mobile robots in presence of obstacles of different shapes. The method is simulated with three robots. Impact of parameters is studied with regard to computation time, obstacle avoidance and travel time.

Introduction

The control of mobile robots is a long-standing subject of research in the robotics domain. A trending application of mobile robotic systems is their use in industrial supply-chains for processing orders and optimizing products' storage and distribution. Companies such as Amazon and the logistic provider IDEA Groupe employ mobile multi-robot systems (Kiva systems, Scallog respectively) for autonomously processing client orders [START_REF] Robarts | Autonomous robots are helping to pack your Amazon orders[END_REF][START_REF]Idea Groupe met en place Scallog pour sa prÃľparation de commandes[END_REF]. Such logistics tasks became increasingly complex as sources of uncertainty, such as human presence, are admitted in the work environment.

One basic requirement for such mobile multirobot systems is the capacity of motion planning, that is, generating admissible configuration and input trajectories connecting two arbitrary states. For solving the motion planning problem, different constraints must be taken into account, in particular, the robot's kinematic and geometric constraints.

The first constraints derive directly from the mobile robot architecture implying, in particular, in nonholonomic constraints. Geometric constraints result from the need of preventing the robot to assume specific configurations in order to avoid collisions, communication lost, etc.

We are particularly interested in solving the problem of planning a trajectory for a team of nonholonomic mobile robots, in a partially known environment occupied by static obstacles, being efficient with respect to the travel time (amount of time to go from initial to goal configuration).

A great amount of work towards collision-free motion planning for cooperative multi-robot systems has been proposed. That work can be split into centralized and distributed approaches. Centralized approaches are usually formulated as an optimal control problem that takes all robots in the team into account at once. This produces solutions closer to the optimal one than distributed approaches. However, the computation time, security vulnerability and communication requirements can make it impracticable, specially for a great number of robots [START_REF] Borrelli | MILP and NLP Techniques for centralized trajectory planning of multiple unmanned air vehicles[END_REF].

Distributed methods based in probabilistic [START_REF] Sanchez | On delaying collision checking in PRM planning: Application to multi-robot coordination[END_REF] and artificial potential fields [START_REF] Khatib | Real-time obstacle avoidance for manipulators and mobile robots[END_REF] approaches, for instance, are computationally fast. However, they deal with collision avoidance as a cost function to be minimized. But rather than having a cost that increases as paths leading to collision are considered, collision avoidance should to be considered as hard constraints of the problem.

Other distributed algorithms are based on receding horizon approaches. In [START_REF] Defoort | Motion planning for cooperative unicycle-type mobile robots with limited sensing ranges: A distributed receding horizon approach[END_REF] a brief comparison of the main distributed receding methods is made as well as the presentation of the base approach extended in our work. In this approach each robot optimizes only its own trajectory at each computation/update horizon. In order to avoid robot-to-robot collisions and lost of communication, neighbor robots exchange information about their intended trajectories before performing the update. Intended trajectories are computed by each robot ignoring constraints that take the other robots into account.

Identified drawbacks of this approach are the dependence on several parameters for achieving realtime performance and good solution optimality, the difficulty to adapt it for handling dynamic obstacles, the impossibility of bringing the robots to a precise goal state and the limited geometric representation of obstacles.

Therefore, in this paper, we propose a motion planning algorithm that extends the approach presented in [START_REF] Defoort | Motion planning for cooperative unicycle-type mobile robots with limited sensing ranges: A distributed receding horizon approach[END_REF]. In this modified algorithm goal states can be precisely reached and more complexes forms of obstacles handled. Furthermore, we perform an investigation about how the method's parameters José M. Mendes Filho, Eric Lucet http://robotics.fel.cvut.cz/demur15/ DEMUR'15 impact a set of performance criteria. Thus, this distributed algorithm is able to find collision-free trajectories and computes the corresponding angular and longitudinal velocities for a multi-robot system in presence of static obstacles perceived by the robots as they evolve in their environment.

This paper is structured as follows. The second section states the problem to be resolved, pointing out the cost function for motion planning and all constraints that need to be respected by the computed solution. The third section explains the algorithm for resolving the motion planning problem and gives some remarks on how to resolve the constrained optimization problems associated with the method. The forth section is dedicated to the results found using this method and the analysis of the specific performance criteria and how they are impacted by the algorithm parameters. Finally, in last section we present our conclusions and perspectives.

Problem Statement

Assumptions

In the development of the approach presented in this paper, the following assumptions are made: (1.) The motion of the multi-robot system begins at the instant t init and goes until the instant t f inal .

( .) The team of robots consists of a set R of B nonholonomic mobile robots. ( .) A robot has precise knowledge of the position and geometric representation of a detected obstacle, i.e., obstacles perception issues are neglected. (8.) A robot can access information about any robot in the team by using a wireless communication link. Latency, communication outages and other problems associated to the that communication link are neglected. (9.) The dynamic model of the multi-robot systems is neglected. (10.) The input (or control) vector of a mobile robot R b is bounded.

Constraints and cost functions

After giving the assumptions in the previous Subsection, we can define the constraints and the cost function for the multi-robot navigation.

( .) The solution of the motion planning problem for the robot R b represented by the pair (q * b (t), u * b (t)) q * b (t) ∈ R n being the solution trajectory for the robot's configuration and u * b (t) ∈ R p the solution trajectory for the robot's input -must satisfy the robots kinematic model equation:

q * b (t) = f (q * b (t), u * b (t)), ∀t ∈ [t init , t f inal ].
(1) where f : R n × R p → R n the is the vector-valued function modeling the robot kinematics. (2.) The planned initial configuration and initial input for the robot R b must be equal to the initial configuration and initial input of R b : .) The planned final configuration and final input for the robot R b must be equal to the goal configuration and goal input for R b :

q * b (t init ) = q b,init , ( 2
)
u * b (t init ) = u b,init . ( 3
) ( 3
q * b (t f inal ) = q b,goal , (4)
u * b (t f inal ) = u b,goal . ( 5
) (4.
) Practical limitations of the input impose the following constraint:

∀t ∈ [t init , t f inal ], ∀i ∈ [1, 2, • • • , p], |u * b,i (t)| ≤ u b,i,max . ( 6
)
( 5.) The cost for the multi-robot system navigation is defined as:

L(q(t), u(t)) = B-1 b=0 L b (q b (t), u b (t), q b,goal , u b,goal ) (7) where L b (q b (t), u b (t), q b,goal , u b,goal
) is the integrated cost for one robot motion planning (see [START_REF] Defoort | Motion planning for cooperative unicycle-type mobile robots with limited sensing ranges: A distributed receding horizon approach[END_REF]). (6.) To ensure collision avoidance with obstacles, the euclidean distance between a robot and an obstacle

(denoted d(R b , O m ) | O m ∈ O b , R b ∈ B) has to satisfy: d(R b , O m ) ≥ 0. ( 8
)
For the circle representation of an obstacle the distance d(R b , O m ) is defined as:

(x b -x Om ) 2 + (y b -y Om ) 2 -ρ b -r Om .
For the convex polygon representation, the distance was calculated using three different definitions, according to the Voronoi region [START_REF] Ericson | Real-Time Collision Detection[END_REF] R b is located.

DEMUR'15

http://robotics.fel.cvut.cz/demur15/ Trajectory Generation Approach (7.) In order to prevent inter-robot collisions, the following constraint must be satisfied:

∀ (R b , R c ) ∈ R × R, b = c, c ∈ C b , d(R b , R c ) -ρ b -ρ c ≥ 0 (9) where d(R b , R c ) = (x b -x c ) 2 + (y b -y c ) 2 and C b
is the set of robots that present a collision risk with R b . (8.) Finally, the need of a communication link between two robots (R b , R c ) yields to the following constraint:

d(R b , R c ) -min(d b,com , d c,com ) ≤ 0 ( 10

Distributed motion planning

Receding horizon approach

Since the environment is progressively perceived by the robots and new obstacles may appear as time passes, planning the whole motion from initial to goal configurations before the beginning of the motion is not a satisfying approach. Planning locally and replanning is more suitable for taking new information, as they come, into account. Besides, the computation cost of finding a motion plan using the first approach may be prohibited high if the planning complexity depends on the distance between start and goal configurations.

Therefore, an on-line motion planner is proposed. In order to do so a receding horizon control approach [START_REF] Keviczky | Decentralized receding horizon control for large scale dynamically decoupled systems[END_REF] is used.

Two fundamental concepts of this approach are the planning horizon T p and update/computation horizon T c . T p is the timespan for which a solution will be computed and T c is the time horizon during which a plan is executed while the next plan, for the next timespan T p , is being computed. The problem of producing a motion plan during a T c time interval is called here a receding horizon planning problem.

For each receding horizon planning problem, the following steps are performed:

Step 1. All robots in the team compute an intended solution trajectory (denoted (q b (t), ûb (t))) by solving a constrained optimization problem.

Coupling constraints ( 9) and ( 10) that involve other robots in the team are ignored.

Step 2. Robots involved in a potential conflict (that is, risk of collision or lost of communication) update their trajectories computed during Step 1 by solving another constrained optimization problem that additionally takes into account coupling constraints [START_REF] Defoort | Contributions à la planification et à la commande pour les robots mobiles coopératifs[END_REF] and [START_REF] Milam | Real-time optimal trajectory generation for constrained dynamical systems[END_REF]. This is done by using the other robots' intended trajectories computed in the previous step as an estimate of those robots' final trajectories. If a robot is not involved in any conflict, Step 2 is not executed and its final solution trajectory is identical to the one estimated in Step 1.

All robots in the team use the same T p and T c for assuring synchronization when exchanging information about their positions and intended trajectories.

For each of these steps and for each robot in the team, one constrained optimization problem is resolved. The cost function to be minimized in those optimization problems is the geodesic distance of a robot's current configuration to its goal configuration. This assures that the robots are driven towards their goal.

This two step scheme is explained in details in [START_REF] Defoort | Motion planning for cooperative unicycle-type mobile robots with limited sensing ranges: A distributed receding horizon approach[END_REF][START_REF] Defoort | Contributions à la planification et à la commande pour les robots mobiles coopératifs[END_REF] where constrained optimization problems associated to the receding horizon optimization problem are formulated.

However, constraints related to the goal configuration and goal input of the motion planning problem are neglected in their method. Constraints (4) and ( 5) are left out of the planning. For taking them into account, a termination procedure is proposed in the following that enables the robots to reach their goal state.

Motion planning termination

After stopping the receding horizon planning algorithm, we propose a termination planning that considers those constraints related to the goal state. This enables the robots to reach their goal states.

The criterion used to pass from the receding horizon planning to the termination planning is based on the distance between goal and current position of the robots. It is defined by the equation 11: [START_REF]SciPy -Scientific Computing Tools for Python[END_REF] This condition ensures that the termination plan will be planned for at least a d min distance from the robot's goal position. This minimal distance is assumed to be sufficient for the robot to reach the goal configuration. Before solving the termination planning problem new parameters for the solution representation and computation are calculated by taking into account the estimate remaining distance and the typical distance traveled for a T p planning horizon. This is done in order to rescale the configuration intended for a previous planning horizon not necessarily equal to the new one. Potentially, this rescaling will decrease the computation time for the termination planning.

d rem ≥ d min + T c • v max
The following pseudo code 1 summarizes the planning algorithm and the Figure 1 illustrates how plans would be generated through time by the algorithm.

In the pseudo code, we see the call of a PlanSec procedure. It corresponds to the resolution of the José M. Mendes Filho, Eric Lucet http://robotics.fel.cvut.cz/demur15/ DEMUR'15 receding horizon planning problem as defined in subsection 3.1.

PlanLastSec is the procedure solving the termination planning problem. This problem is similar to the receding horizon planning problems.

It also has the two steps presented before for computing an intended plan and for updating it, if need be, so conflicts are avoided.

The difference consists in how the optimization problems associated to it are defined.

The optimization problem defined in ( 12) and ( 13) is the problem solved at the first step. The optimization problem associated with the second step is defined ( 14) and (15). Besides, in both new constrained optimal problems, the planning horizon is not a fixed constant as before, instead it is a part of the solution to be found.

Then, for generating the intended plan the following is resolved:

min qb (t),û b (t),T f L b,f (q b (t), ûb (t), q b,goal , u b,goal ) ( 12
)
under the following constraints for τ k = kT c with k the number of receding horizon problems solved before the termination problem:

                   qb (t) = f (q b (t), ûb (t)), ∀t ∈ [τ k , τ k + T f ] qb (τ k ) = q * b (τ k-1 + T c ) ûb (τ k ) = u * b (τ k-1 + T c ) qb (τ k + T f ) = q b,goal ûb (τ k + T f ) = u b,goal |û b,i (t)| ≤ u b,i,max , ∀i ∈ [1, p], ∀t ∈ (τ k , τ k + T f ) d(R b , O m ) ≥ 0, ∀O m ∈ O b , t ∈ (τ k , τ k + T f ) ( 13
) And for generating the final solution:

min q * b (t),u * b (t),T f L b,f (q * b (t), u * b (t), q b,goal , u b,goal ) (14)
under the following constraints:

                                   q * b (t) = f (q * b (t), u * b (t)), ∀t ∈ [τ k , τ k + T f ] q * b (τ k ) = q * b (τ k-1 + T c ) u * b (τ k ) = u * b (τ k-1 + T c ) q * b (τ k + T f ) = q b,goal u * b (τ k + T f ) = u b,goal |u * b,i (t)| ≤ u b,i,max , ∀i ∈ [1, p], ∀t ∈ (τ k , τ k + T f ) d(R b , O m ) ≥ 0, ∀O m ∈ O b , ∀t ∈ (τ k , τ k + T f ) d(R b , R c ) -ρ b -ρ c ≥ 0, ∀R c ∈ C b , ∀t ∈ (τ k , τ k + T f ) d(R b , R d ) -min(d b,com , d d,com ) ≥ 0, ∀R d ∈ D b , ∀t ∈ (τ k , τ k + T f ) d(q * b (t), qb (t)) ≤ ξ, ∀t ∈ (τ k , τ k + T f ) (

Strategies for solving the constrained optimization problems

Flatness property

As explained in [START_REF] Defoort | Contributions à la planification et à la commande pour les robots mobiles coopératifs[END_REF], all mobile robots consisting of a solid block in motion can be modeled as a flat system.

Algorithm 1 Motion planning algorithm

1: procedure Plan 2:

q latest ← q initial 3:

d rem ← |Pos(q f inal )-Pos(q latest )| 4: while d rem ≥ d min + T c • v max do 5: InitSolRepresentation(• • • ) 6: q latest ←PlanSec(• • • ) 7:
d rem ← |Pos(q f inal )-Pos(q latest )| This means that a change of variables is possible in a way that states and inputs of the kinematic model of the mobile robot can be written in terms of a new variable, called flat output (z), and its lth first derivatives. The value of l | l ≤ n depends on the kinematic model of the mobile robot. Therefore, the flat output can completely determine behavior of the system.

Searching for a solution to our problem in the flat space rather than in the actual configuration space of the system presents advantages. It prevents the need for integrating the differential equations of system (constraint 1) and reduces the dimension of the problem of finding an optimal admissible trajectory. After finding (optimal) trajectories in the flat space, it is possible to retrieve back the original configuration and input trajectories.

Parametrization of the flat output by

B-splines Another important aspect of this approach is the parametrization of the flat output trajectory. As done in [START_REF] Milam | Real-time optimal trajectory generation for constrained dynamical systems[END_REF], the use of B-spline functions present interesting properties.

• It is possible to specify a level of continuity C k when using B-splines without additional constraints.

• B-spline presents a local support -changes in parameters values have a local impact on the resulting curve.

The first property is very well suited for parametrizing DEMUR'15 http://robotics.fel.cvut.cz/demur15/ Trajectory Generation Approach the flat output since its lth first derivatives will be needed when computing the system actual state and input trajectories. The second property is important when searching for an admissible solution in the flat space; such parametrization is more efficient and well-conditioned than, for instance, a polynomial parametrization [START_REF] Milam | Real-time optimal trajectory generation for constrained dynamical systems[END_REF]. This choice for parameterizing the flat output introduces a new parameter to be set in the motion planning algorithm which is the number of nonnull knots intervals (denoted simply N knots ). This parameter plus the l value determines how many control points will be used for generating the B-splines.

Optimization solver

The optimization problems associated with finding the solution q * (t), u * (t) are solved using a numerical optimization solver.

For all time dependent constraints time sampling is used. This introduces a new parameter in the algorithm: the time sampling for optimization N s . Each constraint that must be satisfied ∀t ∈ (τ k , τ k + T f ) implies in N s equations.

The need of a solver that supports nonlinear equality and inequality constrains restricts the number of numerical optimization solvers to be considered.

For our initial implementation of the motion planning algorithm, the SLSQP optimizer stood out as a good option. Besides being able to handle nonlinear equality and inequality constrains, its availability in the minimization module of the open-source scientific package Scipy [START_REF]SciPy -Scientific Computing Tools for Python[END_REF] helps to facilitate the motion planner implementation.

However, an error was experienced using this optimizer which uses the SLSQP Optimization subroutine originally implemented by Dieter Kraft [START_REF] Kraft | A software package for sequential quadratic programming[END_REF]. As the cost function value becomes too high (typically for values greater than 10 3 ), the optimization algorithm finishes with the "Positive directional derivative for linesearch" error message. This appears to be a numerical stability problem experienced by other users as discussed in [START_REF]Runtime errors for large gradients[END_REF].

For working around this problem, we proposed a change in the objective functions of the receding horizon optimization problems. This change aims to keep the evaluated cost of the objective function around a known value when close to the optimal solution instead of having a cost depending on the goal configuration (which can be arbitrarily distant from the current position).

We simply exchanged the goal position point in the cost function by a new point computed as follows:

p b,new = p b,goal -p b (τ s-1 + T c ) norm(p b,goal -p b (τ s-1 + T c )) αT p v b,max
Where p b,goal and p b (τ s-1 + T c ) are the positions associated with configurations q b,goal and q b (τ

s-1 +T c ) respectively, α | α ≥ 1, α ∈ R

Simulation results

Results and their analysis for the motion planner presented in the previous sections are presented here.

The trajectory and velocities shown in Figures 2 and3 illustrate a motion planning solution found for a team of three robots. They plan their motion in an environment where three static obstacles are present. Each point along the trajectory line of a robot represents the beginning of a T c update/computation horizon.

It is possible to see on those figures how the planner generates configuration and input trajectories satisfying the constraints associated with the goal states.

In particular, in Figure 2, the resulting plan is computed ignoring coupling constraints (Step 2 is never performed) and consequently two points of collision occur.

A collision-free solution is presented in Figure 3. Specially near the regions were collisions occurred a change in the trajectory is present from Figure 2 to Figure 3 to avoid collision. Complementary, changes in the robots (linear) velocities across charts in both figures can be noticed. Finally, the bottom charts show that the collisions were indeed avoided: inter-robot distances in Figure 3 are greater than or equal to zero all along the simulation.

For performing these two previous simulations, a reasonable number of parameters have to be set. These parameters can be categorized into two groups. Algorithm related parameters and the optimization solver related ones. Among the former group, the most important ones are:

• The number of sample for time discretization (N s );

• The number of internal knots for the B-splines curves (N knots );

• The planning horizon for the sliding window (T p );

• The computation horizon (T c ).

• The detection radius of the robot (d sen ).

The latter kind depends on the numeric optimization solver adopted. However, since most of them are iterative methods, it is common to have at least a maximum number of iterations and a stop condition parameters.

This considerable number of parameters makes the search for a satisfactory set of parameters' values a laborious task.

Therefore, it is important to have a better understanding of how some performance criteria are impacted by the changes in algorithm parameters. Inter-robot distances throughout simulation

d(R 0 ,R 1 )-ρ 0 -ρ 1 d(R 0 ,R 2 )-ρ 0 -ρ 2 d(R 1 ,R 2 )-ρ 1 -ρ 2

Parameters' impact

Three criteria considered important for the validation of this method were studied: Maximum computation time during the planning over the computation horizon (M CT /T c ratio); Obstacle penetration area (P ); Travel time (T tot ). Different parameters configuration and scenarios where tested in order to highlight how they influence those criteria.

Maximum computation time over computation horizon M CT /T c

The significance of this criterion lays in the need of assuring the real-time property of this algorithm. In a real implementation of this approach the computation horizon would have always to be superior than the maximum time took for computing a plan.

Table 1 summarizes one of the scenarios studied for a single robot. Results obtained from simulations in that scenario are presented in Figure 4, for different parameters set.

Each dot along the curves corresponds to the Inter-robot distances throughout simulation Rather than observing the absolute values, it is interesting to analyze the impact of changes in the parameters values. In particular, an increasing number of N s increases M CT /T c for a given T c /T p . Similarly, an increasing of M CT /T c as the number of internal knots N knots increases from charts 4a to 4c is noticed.

d(R 0 ,R 1 )-ρ 0 -ρ 1 d(R 0 ,R 2 )-ρ 0 -ρ 2 d(R 1 ,R 2 )-ρ 1 -ρ 2
Further analyses of those data show that finding the solution using the SLSPQ method requires O(N Another parameter having direct impact on the M CT /T c ratio is the detection radius of the robot's sensors. As the detection radius of the robot increases, more obstacles are seen at once which, in turn, increases the number of constraints in the optimization problems. The impact of increasing the detection radius d sen in the M CT /T c ratio can be seen in Figure 5 for a scenario with seven obstacles. The computation time stops increasing as soon as the A solution where the planned trajectory does not pass through an object at any instant of time gives P = 0. The solution quality decrease with increasing P. However, since time sampling is performed during the optimization, P is usually greater than zero. A way of assuring P = 0 would be to increase the obstacles radius computed by the robot's perception system by the maximum distance that the robot can run within the time spam T p /N s . However simple, this approach represents a loss of optimality and is not considered in this work.

N s =10 N s =11 N s =12 N s =13 N s =14 N s =15 (
N s =10 N s =11 N s =12 N s =13 N s =14 N s =15 (
N s =10 N s =11 N s =12 N s =13 N s =14 N s =15 (c) . Six internal knots
It is relevant then to observe the impact of the algorithm parameters in the obstacle penetration area. T c /T p ratio, N knots and d sen impact on this criteria is only significant for degraded cases, meaning that around typical values those parameters do not change P significantly. However, time sampling N s is a relevant parameter. Figure 6 shows the penetration area decreasing as the number of samples increases. Obstacle penetration decreasing as sampling increases

Travel time T tot

Another complementary metric for characterizing solution quality is the travel time T tot . Analyses of data from several simulations show a tendency that for a given value of N knots , N s and T c the travel time decreases as the planning horizon T p decreases. This can be explained by the simple fact that for a given T c , a more optimal solution (in terms of travel time) can be found if the planning horizon T p is smaller. José M. Mendes Filho, Eric Lucet http://robotics.fel.cvut.cz/demur15/

DEMUR'15

Another relevant observation is that the overall travel time is shorter for smaller N s 's. This misleading improvement does take into account the fact that the fewer the samples the greater will be the obstacle penetration area as shown previously in Figure 6.

Furthermore, Figure 7 shows travel time invariance for changes in the detection radius far from degraded values that are too small. This points out that a local knowledge of the environment provides enough information for finding good solutions. Total execution time and detection radius relationship

Conclusions

A distributed motion planner based on a receding horizon approach, modified for taking into account termination constraints, was proposed. Near the goal configuration neighborhood, the receding horizon approach is finished and a termination planning problem is solved for bringing the robots to their precise final state. The problem is stated as a constrained optimization problem. It minimizes the time for reaching a goal configuration through a collision-free trajectory securing communication between robots. Circle and convex polygon representation of obstacles are supported. Key techniques for implementing the motion planner are: system flatness property, B-spline parameterization of the flat output and SLSQP optimizer. Finally, solutions using this planner for different scenarios were generated in order to validate the method. Impact of different parameters on computation time and quality of the solution was analyzed. Future work will be performed in physics simulation environment where dynamics is taken into account as well as sensors models and communication latency.

( 3 .

 3 ) A robot (denoted R b , R b ∈ R, b ∈ {0, . . . , B-1}) is geometrically represented by a circle of radius ρ b centered at (x b , y b ).

( 4 .

 4 ) All obstacles in the environment are considered static. They can be represented by a set O of M static obstacles.

(

  ) with d b,com , d c,com the communication link reach of each robot and D b is the set of robots that present a communication lost risk with R b .

  15) A possible definition for the L b,f cost function present in the equations above can be simply T f . The sets O b , C b and D b are functions of τ k .

Figure 1 .

 1 Figure 1. Receding horizon scheme with termination plan. The timespan T f represents the duration of the plan for reaching the goal configuration.

  is a constant for controlling how far from the current position the new point is placed, the product T p v b,max the maximum possible distance covered by R b during a planning horizon and s | s ∈ [0, k), s ∈ N the current receding horizon problem index.

  http://robotics.fel.cvut.cz/demur15/

Figure 2 .

 2 Figure 2. Motion planning solution without collision handling

Figure 3 .

 3 Figure 3. Motion planning solution with collision handling

Figure 4 .

 4 Figure 4. Three obstacles scenario simulations

Figure 5 . 4 . 1 . 2 .

 5412 Figure 5. Increasing of detection radius and impact on a M T C/Tc ratio

P (cm 2 )

 2 Time sampling and obstacle penetration relationship P(N_s) Fitted Curve (9395.01 exp(-0.48 N s ) +3.13)

Figure 6 .

 6 Figure 6.Obstacle penetration decreasing as sampling increases

Figure 7 .

 7 Figure 7. Increasing of detection radius and impact on Ttot

  Om , y Om ).

	.) An obstacle (denoted O m , O m ∈ O,
	m ∈ {0, . . . , M -1}) is geometrically represented
	either as a circle or as a convex polygon. In the
	case of a circle its radius is denoted r Om centered
	at (x (6.) For a given instant t k ∈ [t init , t f inal ], any
	obstacle O m is considered detected by the robot
	R b whenever the distance between their geometric
	centers is less than or equal to the detection radius
	d b,sen of the robot R b . Therefore, this obstacle O m
	is part of the set O b (O b ⊂ O) of detected obstacles.

Table 1 .

 1 3 knots ) and O(N s ) time. Although augmenting N knots can yield to an impractical computation time, typical N knots values did not need to exceed 10 in our simulations, which is a sufficiently small value. Values for scenario definition

	DEMUR'15
	http://robotics.fel.cvut.cz/demur15/	Trajectory Generation Approach
			v max	1.00 m/s
			ω max	5.00 rad/s
			q inital [-0.05 0.00 π/2] T
			q f inal	[0.10 7.00 π/2] T
			u initial	[0.00 0.00] T
			u goal	[0.00 0.00] T
			O 0	[0.55 1.91 0.31]
			O 1	[-0.08 3.65 0.32]
			O 2	[0.38 4.65 0.16]
	MCT/T	0.2 0.2 0.4 0.6 0.8 1.0 1.8 1.6 1.4 1.2	0.3	0.4	T c /T p	0.5	0.6

c

Computation cost behavior

  a) . Four internal knots

	MCT/T	0.2 0.2 0.4 0.6 0.8 1.0 1.8 1.6 1.4 1.2	0.3	0.4	T c /T p	0.5	0.6

c

Computation cost behavior