Thomas Rousseau
email: thomas.rousseau@ls2n.fr

Nicolò Pedemonte
email: nicolo.pedemonte@irt-jules-verne.fr

Stéphane Caro
email: stephane.caro@ls2n.fr

François Chaumette
email: francois.chaumette@inria.fr

Stability Analysis of Profile Following by a CDPR using Distance and Vision Sensors ⋆

Keywords: Stability Analysis, Sensor-based control, Cable-Driven Parallel Robot

Cable-Driven Parallel Robots (CDPRs) form a class of robots well-adapted to large workspaces since they replace rigid links by cables. However, they lack in positioning accuracy. In a previous work, a control law has been proposed to enable a CDPR to perform a profile-following task, based on the data measured by two different and redundant sensors that are fused using the Gradient Projection Method (GPM). However, its robustness had not been assessed yet. This paper analyzes the stability of such a control law in function of the systematic errors on the parameters of the system. Numerical simulations show that the characteristics of the system ensure the stability of the control law in the robot workspace even in the presence of significant errors, provided the initial angle between the surface to follow and the moving-platform is smaller than 11°.

Introduction

Cable-driven parallel robots (CDPRs) are robots where the end-effector, named the moving-platform (MP), is actuated by a set of cables instead of rigid links. The length of the cables is usually controlled by a drum-pulley system, actuated by motors fixed to the main frame. These robots are able to carry loads in large workspaces and could thus help in large industrial assembly tasks. However, their accuracy, especially for positioning with regards to a given object, needs to be improved before they can find a place on production lines. While many research works brought significant improvements to CDPR models [START_REF] Pott | Influence of pulley kinematics on cable-driven parallel robots[END_REF][START_REF] Mishra | Forward kinematics for suspended under-actuated cable-driven parallel robots with elastic cables: A neural network approach[END_REF][START_REF] Schmidt | Modeling techniques and reliable real-time implementation of kinematics for cable-driven parallel robots using polymer fiber cables[END_REF], sensor-based control, and especially vision-based control [START_REF] Dallej | Towards vision-based control of cable-driven parallel robots[END_REF][START_REF] Zake | Vision-based control and stability analysis of a cable-driven parallel robot[END_REF][START_REF] Ramadour | Grasping objects with a cable-driven parallel robot designed for transfer operation by visual servoing[END_REF], can also enable such robots to reach a better precision, without needing the tedious computations required for the Forward Kinematics. However, vision-based control has some important limitations: embedded cameras suffer from a limited field of view, while external cameras provide a global overview of the workspace but at the cost of worse accuracy. A solution to these shortcomings was proposed in [START_REF] Rousseau | Constant Distance and Orientation Following of an Unknown Surface with a Cable-Driven Parallel Robot[END_REF] where measurements from both distance sensors and an external camera are combined using the Gradient Projection Method (GPM) [START_REF] Nakamura | Task-priority based redundancy control of robot manipulators[END_REF][START_REF] Siciliano | Kinematic control of redundant robot manipulators: A tutorial[END_REF] to prioritize tasks. The robustness of this control law was not considered in [START_REF] Rousseau | Constant Distance and Orientation Following of an Unknown Surface with a Cable-Driven Parallel Robot[END_REF]. While the stability analysis for N-tasks problems was investigated in [START_REF] Antonelli | Stability analysis for prioritized closed-loop inverse kinematic algorithms for redundant robotic systems[END_REF], the proposed stability criterion is only valid in the ideal case while a specific study based on the quality of the estimations and possible sensor biases is required in our case. The robustness of CDPR control for one pure visual task has been tackled in [START_REF] Zake | Vision-based control and stability analysis of a cable-driven parallel robot[END_REF]. A novel workspace named Control Stability Workspace (CSW), enclosing the stable poses for given perturbation bounds, was also introduced in [START_REF] Zake | Control Stability Workspace for a Cable-Driven Parallel Robot Controlled by Visual Servoing[END_REF].

In this paper, the robustness of the GPM control for a profile following task is assessed. First, a stability criterion is derived from the control law of the system. Then, the impact on stability of all the identified error sources is investigated. Finally, a combined analysis is run in the robot operating conditions to confirm its robustness. Since experimental results have already been presented in [START_REF] Rousseau | Constant Distance and Orientation Following of an Unknown Surface with a Cable-Driven Parallel Robot[END_REF], only simulation results are discussed in this paper. A view of the experimental platform used in this paper, ACROBOT, a six degrees of freedom (DoF) suspended CDPR, is shown in Fig. 1.

Fig. 1: ACROBOT, a CDPR located at IRT Jules Verne

This paper is organized as follows: Section 2 presents the two sensors and their kinematics modeling. The control law is recalled from [START_REF] Rousseau | Constant Distance and Orientation Following of an Unknown Surface with a Cable-Driven Parallel Robot[END_REF] in Section 3. Section 4 focuses on the design of the stability criterion and Section 5 presents the numerical results of the stability analysis. Conclusions are drawn in Section 6.

Task definition

The CDPR considered in this paper is subject to a control law enabling its MP to perform a profile following task. This task can be divided into two subtasks: maintaining the MP parallel to the followed surface at a constant distance, denoted T 1 , and moving from a starting pose to a target pose along the profile, denoted T 2 . The control of the robot stems from sensor-based control and combines measurements provided by two types of sensors: an array of three distance sensors embedded on the MP, and an external camera located in front of the robot. The boarded sensors and the surface to follow are shown in Fig. 2.

Sensors

Distance sensors Three ultrasonic distance sensors are placed on the bottom of the MP, facing downwards along the z-axis of the MP frame F p , as shown in Fig. 2. Since all sensors are coplanar and their directions are parallel, for the MP to be locally parallel to the surface, the three sensors must return the same distance. The vector of distances measured by the sensors is denoted as s = [s 1 , s 2 , s 3 ] T , s i being the distance measured by the ith sensor S i . The time derivative of the measured distance vector, noted ṡ, is then expressed as a function of the MP twist v p :

ṡ = p L 1 v p (1)
with

v p = [ν p T , ω p T ]
T where ν p is the MP linear velocity vector and ω p is its angular velocity vector, expressed in F p . p L 1 is the interaction matrix of the distance sensors [START_REF] Samson | Robot control: the task function approach[END_REF][START_REF] Chaumette | Visual servoing[END_REF]. Using a constant value of the interaction matrix corresponding to the desired configuration is a usual choice in sensor-based control [START_REF] Chaumette | Visual servoing[END_REF]. The interaction matrix p L 1 thus only depends on the respective positions (x i , y i ) of the distance sensors expressed in F p , and can be evaluated for the particular case s = s * , where s * = [s * , s * , s * ] T and with s * = 0.2 m the common reference distance so that the MP is locally parallel to the followed surface:

p L 1|s=s * =   0 0 -1 -y 1 x 1 0 0 0 -1 -y 2 x 2 0 0 0 -1 -y 3 x 3 0   (2)
External camera The images acquired by the external camera are processed by an algorithm recognizing several AprilTags, placed on the MP as well as on the starting and target points of the trajectory. The relative pose of the reference frames attached to these April Tags is then estimated with regard to the camera frame F c . Then, this pose is expressed in the world frame F w , and, as for Pose-Based Visual Servoing (PBVS), the visual features are selected as:

s v = w t p θu (3)
where w t p is the position vector of the MP expressed in F w and θu is the axisangle representation of the rotation matrix p * R p . The target pose is given by:

s * v = w t * p 0 3 (4)
where w t * p is the position of the target AprilTag in F w . The interaction matrix w L 2 of these visual features can be found in [START_REF] Chaumette | Visual servoing[END_REF]:

w L 2 = I 3 -[ w t p ] × 0 3 L ω (5)
with [ w t p ] × the skew-symmetric matrix of w t p . L ω can be approximated by I 3 , the (3 × 3) identity matrix. To express this matrix in F p , the adjoint matrix w V p is used [START_REF] Zake | Control Stability Workspace for a Cable-Driven Parallel Robot Controlled by Visual Servoing[END_REF]:

w V p = w R p [ w t p ] × w R p 0 3 w R p (6)
The change of reference frame is then performed:

p L 2 = w L 2 w V p (7)
This interaction matrix only depends on the estimated position and orientation of the MP in F w .

Fusion using the Gradient Projection Method

The fusion between the data from the external camera and the three distance sensors is performed with the GPM. In this method, a first main task is considered as a priority, while only the part of the second task compatible with the realization of the first one is retained [START_REF] Siciliano | Kinematic control of redundant robot manipulators: A tutorial[END_REF].

Main task

The main task T 1 is selected from the distance sensor measurements, such that the task error e 1 is given by:

e 1 = s -s * (8)
The next step is to determine the Jacobian of the main task. As described above in Section 2.1, the time variation of the main task is related to the MP twist v p by [START_REF] Pott | Influence of pulley kinematics on cable-driven parallel robots[END_REF]. As the inputs of the low-level controller of the CDPR are the velocities ˙l of the m cable lengths, they are related to v p by: ˙l = A v p [START_REF] Siciliano | Kinematic control of redundant robot manipulators: A tutorial[END_REF] where the forward Jacobian matrix A of the CDPR is given by [START_REF] Zake | Vision-based control and stability analysis of a cable-driven parallel robot[END_REF]:

A = u 1 . . . u i . . . u m b 1 × u 1 . . . b i × u i . . . b m × u m T ( 10
)
where u i are the cable direction unit vectors, pointing from the exit points (pulleys) of the CDPR to the anchor points placed on the MP, and where b i are the position vectors of the anchors points, known from the platform design. All these vectors are expressed in F p . u i vectors are calculated geometrically from the MP pose expressed in F w using the straight and inelastic cable model [START_REF] Zake | Vision-based control and stability analysis of a cable-driven parallel robot[END_REF].

Equations ( 1) and ( 9) can be combined to obtain:

ė1 = J 1 ˙l ( 11
)
where the Jacobian J 1 of the main task is a (3 × m) matrix, given by

J 1 = L 1 A + (12)
with A + the Moore-Penrose pseudo-inverse of A.

Secondary task T 2 consists in minimizing the error between the current and the desired MP poses, expressed in F w , and noted as follows:

e 2 = s v -s * v (13)
From ( 7) and ( 9), its Jacobian is a (6 × m) matrix given by:

J 2 = w L 2 w V p A + (14)
Similarly to [START_REF] Zake | Control Stability Workspace for a Cable-Driven Parallel Robot Controlled by Visual Servoing[END_REF], ė2 satisfies the following relationship: ė2 = J 2 ˙l (15)

Control strategy

In practice, the true Jacobians J 1 and J 2 can be subject to approximations, which is the case for the form given in (2), noise measurements, and calibration errors. The Jacobians used in the control scheme are thus different from the true ones and are denoted J in the following. The control law, pictured in Fig. 3 and performing the fusion described above, is detailed in [START_REF] Rousseau | Constant Distance and Orientation Following of an Unknown Surface with a Cable-Driven Parallel Robot[END_REF]. It relies on the projection operator P 1 on the kernel of the primary task given by:

P 1 = I m -J + 1 J 1 (16)
The control law is expressed in terms of e 1 and e 2 , with the gain ratio

r λ = λ 2 /λ 1 ˙l = -λ 1 J + 1 e 1 -λ 1 (J 2 P 1 ) + (r λ e 2 -J 2 J + 1 e 1 ) (17)
To simplify the notations, the term ẽ2 = r λ e 2 -J 2 J + 1 e 1 and matrix J2 = J 2 P 1 are defined such that the control law has the simple form:

˙l = -λ 1 ( J + 1 e 1 + J+ 2 ẽ2 ) (18)
Fig. 3: Control scheme of the CDPR with the fusion of the distance and pose errors.

Stability criterion

The stability of the control scheme can be classicaly assessed using a Lyapunov stability analysis. Since the main task is of dimension three, the projection of the secondary task on the kernel of the main task causes the secondary task to lose three degrees of freedom. Hence, the stability of the control law cannot be assessed directly by considering ||e 1 || 2 + ||e 2 || 2 as candidate Lyapunov function.

Secondary task error redefinition

To solve the analysis problem, a three-dimensional secondary task e 2r can be defined by projecting the six-dimensional secondary task e 2 on the null space of e 1 , expressed in the MP frame. First, the projected error e ′ 2 is considered:

e ′ 2 = p P 1 e 2 , with p P 1 = I 6 -p L + 1|s=s * p L 1|s=s * (19)
Considering only the three non-zero terms of e ′ 2 , the reduced error vector e 2r can be defined as:

e 2r =   w x p -w x * p w y p -w y * p θu z   (20)
The stacked vectors e 1 and e 2r hence form a six-dimensional error. The corresponding stability problem can be assessed since the dimension of the variables matches the number of DoF of the MP. The modified interaction matrix p L 2r is then obtained by selecting the rows corresponding to the retained features and the corresponding modified matrices J 2r and J2r are computed similarly by substituting L 2 by L 2r in their definition. A modified control law follows:

˙l = -λ 1 ( J + 1 e 1 + J+ 2r ẽ2r ) (21)

Lyapunov stability criterion

By considering the reduced secondary task, the Lyapunov candidate function L, strictly positive and continuously differentiable, is defined as:

L = 1 2 ||e 1 || 2 + 1 2 ||e 2r || 2 (22)
Its time derivative is nothing but:

L = e T 1 ė1 + e T 2r ė2r = e T 1 J 1 ˙l + e T 2r J 2r ˙l (23)
Assuming the low level robot controller is able to perfectly apply the cable velocity vector computed by the control law (21), we obtain:

L = -λ 1 e T 1 J 1 J + 1 e 1 -e T 1 J 1 J+ 2r J 2r J + 1 e 1 + r λ e T 1 J
The previous equation can be expressed in a matrix form such that the Lyapunov candidate function derivative L is expressed as:

L = -λ 1 [e T 1 e T 2r ] Π e 1 e 2r (25)
which leads to the stability criterion Π> 0 with:

Π =   J 1 (I m - J+ 2r J 2r ) J + 1 r λ J 1 J+ 2r J 2r (I m - J+ 2r J 2r ) J + 1 r λ J 2r J+ 2r   (26)
In the ideal case, that is when the Jacobians used in the control scheme correspond to the real ones, i.e., J 1 = J 1 and J2 = J2 , it follows that J 1 J+ 2r = 0 3 and the expression of matrix Π becomes:

Π = I 3 0 3 J 2r J + 1 -J 2r J+ 2r J 2r J + 1 r λ J 2r J+ 2r ( 27
)
As expected, the resulting expression in the ideal case matches with the relationship already established in [START_REF] Antonelli | Stability analysis for prioritized closed-loop inverse kinematic algorithms for redundant robotic systems[END_REF].

Stability analysis

The sign of the Lyapunov function derivative L depends on the sign of the criterion Π. If the criterion is positive definite, the derivative is negative for all input vectors, thus guaranteeing the stability of the control law. This sign depends on the quality of the estimation of a series of parameters. These are the position error of the exit points of the pulleys and the anchor points of the cables on the MP δa i and δb i , the vision-based pose estimate error δ w t p , the positionning error of the distance sensors in the platform reference frame δ L s,d , the relative orientation of the surface to follow with regard to the platform frame, given by the vectors n Ti , as well as the gain ratio r λ . Since there is an infinity of combinations of poses and parameter values, the path P validated experimentally in [START_REF] Rousseau | Constant Distance and Orientation Following of an Unknown Surface with a Cable-Driven Parallel Robot[END_REF] on the ACROBOT is chosen for this study. An analysis parameter by parameter is run to determine which ones have an important impact on the stability and which ones can be neglected. Then, a combined analysis using known parameter values is conducted to provide us with the expected stability domain of the robot during its standard behaviour. The analytical sign analysis of Π is rather complex, due the multiple interactions between the parameters. A numerical analysis is thus proposed. For this analysis, an initial angle between the MP and the tangent plane to the surface followed of 15°was retained.

CDPR geometrical parameters

The precision of the CDPR geometry parameters a i and b i values has an impact on A and hence on Π. Considering only errors in the estimation of these parameters along the path P, the stability margin versus both a i and b i parameters is shown Fig. 4. The control law remains stable up to an error of 0.08 m for a i and to 0.10 m for b i . The stability isocontours, i.e., the lines where the lowest eigenvalue remains constant, closely follow the altitude profile of the path, hinting at a significant dependency on this parameter regarding the stability of the control law. For this experiment system, the accuracy of these parameters is respectively 0.01 m for δa i and 0.008 m for δb i . Hence, the stability of the robot is not threatened by these errors.

Vision-based pose estimation The pose estimation impacts the stability since it is involved in the estimation of the robot Jacobian A and p L 2 , necessary to compute the estimates of J 1 and J 2r . The impact of the position measurement bias resulting of a low image resolution is shown Fig. 5. However, this bias is not uniformly distributed along the three directions, and the estimated position error is larger along the depth of the image observed by the camera. The norm of the position error plotted in Fig. 5 was thus distributed accordingly. The impact of the orientation error about the three directions of F p was computed along P, and is represented for a central position in the robot cell in Fig. 7. In both cases, the accuracy of the pose-estimation, 70 mm for the position and 3 degrees for the orientation at worst, does not overshoot the stability limits. However, the resulting margin is lower than the one obtained for the CDPR geometrical parameters errors and pose-estimation errors have a larger stability impact. [START_REF] Mishra | Forward kinematics for suspended under-actuated cable-driven parallel robots with elastic cables: A neural network approach[END_REF]. The impact of a measurement error on the stability of the control law is shown Fig. 6. From this figure, it is clear that the quality of such measurements significantly impacts the stability since the maximum norm of the error ensuring stability for these parameters is twice lower than for the parameters studied above. However, an adequate accuracy on these parameters can easily be achieved with a robust design of the MP or following a characterization of this MP. For ACROBOT, this accuracy is estimated to 1 mm. Gain ratio Eventually, the gain ratio between the main task and the secondary task is also involved in the computation of the stability criterion Π 2 . The corresponding plot is not shown due to lack of space but the result obtained shows that this parameter cannot prevent the stability of the control law. However, the

Combination of parameters

To analyze the combination of these factors, another analysis is conducted, combining all the error sources. A series of 20.000 trajectories was run, selecting randomly the errors applied. Each of these parameters was assigned an error box whose side corresponds to its respective estimated accuracy. For each trajectory and for each parameter discussed above, a random vertex of the corresponding error box is selected. The distance measurements returned by the distance sensors are generated randomly with a Gaussian centered on µ = 0.2 m and σ = 0.025. The distance difference between the three sensors generates an angle α between the normal n S and the normal to the tangent plane n T , equivalent to the angle between the plane x p O p y p and the tangent plane.

The results of this analysis with regard to the angle α are presented Fig. 8. In blue: stable configurations and in red: unstable configurations.

When observing the minimum eigenvalues of the stability criterion versus the angle α, all configurations are stable when α ≤ 11°. However, after this threshold, the proportion of stable configurations decreases when the angle increases. No stable configurations were found when α ≥ 33°. This implies that an initial configuration of the system with α lower than 11°and that, while executing the two tasks, the tracking error with regard to the main task remains below 11°is a sufficient condition for the convergence of the system. On the tested trajectory P, the orientation of the tracking error does not exceed 8°, and the stability of the control law with the measured parameter errors is thus ensured.

Conclusion

This paper proposed a stability analysis of a sensor-based control law for profile following using a CDPR. First, a stability criterion was derived, then a numerical parameter by parameter analysis was conducted, highlighting the most critical parameters that ensure the stability of the system. The most critical parameters are the coordinates of the distance sensor in the MP frame as well as the bias on the vision-based pose estimation. A combined analysis taking into account the most relevant parameters for various configurations of the robot was also presented. This analysis shows that all these configurations for the biases of the relevant parameters of the prototype ACROBOT are stable when the angle between the MP and the tangent plane to the target surface is smaller than 11°. A sufficient condition for the stability of the system on the initial configuration of the robot and on the tracking error hence results from this analysis. This implies that the response time of the main task must remain short so that the moving MP during the trajectory following does not overshoot the stability angle threshold. While these results are specific to the ACROBOT prototype, the method is transferable. Stability domains of other CDPRs using the same control scheme can be computed similarly if the relevant parameters are given.

Fig. 2 :

 2 Fig. 2: Parametrization of the three distance sensors embedded on the MP and facing the target surface. The external camera observes the MP.

( a )

 a Stability versus the norm of the exit point error δa i . (b) Stability of the robot versus the norm of the anchor point error δb i .

Fig. 4 :Fig. 5 :Fig. 6 :

 456 Fig. 4: Stability of the robot along P versus the CDPR geometrical parameters.

Fig. 7 :

 7 Fig. 7: Stability of the robot versus the three components of the orientation estimation error δθu

Fig. 8 :

 8 Fig. 8: Stability margin of the tested configurations versus orientation angle α.In blue: stable configurations and in red: unstable configurations.

⋆ This work is supported by IRT Jules Verne in the framework of the PERFORM program.


Tahir Rasheed
email: tahir.rasheed@univ-nantes.fr

Loic Michel
email: loic.michel@ec-nantes.fr

Stéphane Caro
email: stephane.caro@ls2n.fr

Jean-Pierre Barbot
email: barbot@ensea.fr

Yannick Aoustin
email: yannick.aoustin@univ-nantes.frwww.ls2n.fr

Dynamic Parameter Identification for Cable-Driven Parallel Robots ⋆

Keywords: Cable-Driven Parallel Robots, Inverse Dynamic Model, Dynamic Parameter Identification, Least-squares estimator, Semi-Implicit Euler Homogeneous Differentiator, Experiments

This paper presents an initial work of identifying the dynamic parameters for cable-driven parallel robot (CDPR), CRAFT with a rigorous protocol. An orbital trajectory of the platform is designed in order to get a pure translation movement of the platform. This trajectory evolves in a plane and allows the identification of four dynamic parameters, the mass of the platform and the first three moments along the three main axes. The identification results obtained respectively from the data from the experimental measurements of the moving platform, and from the experimental measurements treated with the application of a semi-implicit homogeneous differentiator, show that in spite of the complexity of CRAFT an identification of all its essential dynamic parameters is possible. Moreover, in the long term, an online identification of CRAFT during handling tasks is envisaged.

Introduction

The objective is to initiate a scientific approach in order to identify the essential dynamic parameters of a cable-driven parallel robot (CDPR). A CDPR consists of a moving-platform (M P ) that is connected to a rigid frame by means of cables and actuators, the latter being generally mounted on the ground. These robots are very attractive for handling tasks [START_REF] Picard | Control Strategies for a Cable-Driven Parallel Robot with Varying Payload Information[END_REF] because of their low inertia, a higher payload to weight ratio and a large workspace compared to conventional manipulator robots with articulated rigid limbs. Their possible application fields can be industrial, or dedicated to search-and-rescue operations. To deal with various restrictions on cable tensions, cable elasticity, collisions and obstacle avoidance, over-actuation of the M P is actually a challenging scientific problem [START_REF] Gouttefarde | Geometry selection of a redundantly actuated cable-suspended parallel robot[END_REF], [START_REF] Hussein | Smallest maximum cable tension determination for cable-driven parallel robots[END_REF].

The control of a CDPR is complex because, among other things, whatever its movement, the tension of its cables must always be positive. There are still open problems in the control of CDPRs. One key example of an open problem is that control design requires a consistent dynamic model with a relatively good knowledge of the dynamic parameters of the robot such as the terms of inertia, mass, and even friction at the level of the actuators, which manage the winding of the cables. Several contributions about identification of a robot cable exist. Kraus et al. [START_REF] Kraus | System identification and cable force control for a cable-driven parallel robot with industrial servo drives[END_REF] present an identification method for the complete actuator of a cable robot. A second-order system is established with a dead time as an analogous model. An inverse dynamic method (IDM) is used to identify the model parameters of a cable-driven finger joint for surgical robot [START_REF] Yu | Dynamical model and experimental identification of a cable-driven finger joint for surgical robot[END_REF] . But To our best knowledge the identification parameters of the handling platform of a cable robot during a movement in space is not currently investigated.

A CDPR, named CRAFT and located at LS2N, Centrale Nantes campus, is equipped with eight actuators and a M P . Each motor has an encoder sensor measuring the angular velocity of its output shaft allowing to evaluate the performances of the differentiation solutions. The M P has six degrees of freedom. This M P is thus over-actuated [START_REF] Mishra | Kinematic stability based afg-rrt* path planning for cable-driven parallel robots[END_REF].

The identification of the parameters of a robot is essential to evaluate its behavior in simulation or to synthesize a control based model. The objective of the identification is not to find physical parameters of a mechanism such as a robot, whose value is the most exact possible but to provide a modeling tool that is consistent in order to make simulation, control or other. The identification methodology is usually based on an inverse dynamic model that is linear in the dynamic parameters to be identified [START_REF] Gautier | On the identification of the inertial parameters of robots[END_REF], [START_REF] Janot | A generic instrumental variable approach for industrial robot identification[END_REF], or [START_REF] Ardiani | Contribution to the parametric identification of dynamic models: Application to collaborative robotics[END_REF]. Any matrix inversion, which is generally a generator of numerical problems, is thus avoided. However the identification of dynamic parameters is not easy when the dynamic model is complex. Among the many challenges is the search for exciting trajectories that identify the largest number of dynamic parameters and overcome noise problems [START_REF] Khalil | Modeling, identification and control of robots[END_REF]. However the identification of dynamic parameters is not easy when the dynamic model is complex. The dynamic behavior of the CRAFT that moves in 3D space with rotation and translation combinations effectively requires a complex modeling with a lot of inertial parameters. In addition, one or more cables may slacken unexpectedly. As a consequence to initiate an identification work, only the model of the platform maneuvered by the cables is considered. Moreover, the defined trajectories generate only translation movements of the platform, i.e. without any rotational movement of the platform with respect to itself, in a horizontal plane. The model to describe the dynamic behavior of CRAFT becomes simpler with only four dynamic parameters to identify. To manage the noise problems semi-implicit Euler homogeneous differentiators are proposed (see e.g. [START_REF] Michel | An experimental investigation of discretized homogeneous differentiators: pneumatic actuator case[END_REF] for a comparison of the techniques) and [START_REF] Michel | Semi-implicit homogeneous euler differentiator for a second-order system: Validation on real data[END_REF]. Off-line identification was conducted in order to compare, respectively from the planned trajectories of the platform, the measured filtered signals and their usual Euler discretization, and measured signals processed with an original semi-implicit homogeneous differentiators without extra filter. The main advantage of this differentiator is to offer the perspective of a real-time identification of dynamic parameters ofCRAFT. The contribution of this work is based on three parts: Elaboration of a protocol in order to work with a dynamic model for translational movements of the platform and to have only four parameters to identify and experimental definition of this movement; Identification of the four parameters starting respectively from the reference trajectory, the tension measurement of the cable tension, the position measurement and the application of the usual Euler discretization to obtain the speed and acceleration signals. Furthermore these speed and acceleration signals are compared with those obtained an original semi-implicit Euler homogeneous differentiator. The long-term idea of this work is to develop a complete framework that brings together numerical identification and derivation tools to identify the dynamic parameters of the on-line CRAFT as it performs handling tasks. The remaining of the paper is outlined as follows. Section 3 is devoted to the presentation of CRAFT namely, its geometric, kinematic, dynamic model and its adaptation to the presented strategy. The identification methodology, which is based on the least-squares (LS) estimator is presented in Sec. 4, with the solving way. A semi-implicit homogeneous differentiator, which is considered in this study to get the velocities and the accelerations from the measured position, is presented Sec. 5. The experimental results are presented in Sec. 6. Conclusions and future work are offered in Sec. 7.

Robot CRAFT

The cable-driven parallel robot prototype, named CRAFT is located at LS2N, France. The base frame of CRAFT is 4 m long, 3.5 m wide, and 2.7 m high, see Figure 1). The three-DoF translational and the three-DoF rotational motions of its suspended M P are controlled with eight cables being respectively wound around eight actuated reels fixed to the ground. The M P is 0.28 m long, 0.28 m wide, and 0.2 m high, its overall mass being equal to 5 kg. For each of the eight electrical motors an encoder sensor measures the angular variable of its shaft. The eight motors are equipped with a gearbox reducer of ratio n = 8. The measured value is divided by n in order to obtain the angular position of the output shaft of the gearbox reducer. The robot CRAFT has no tachometer.

The main hardware of the prototype consists of a PC (equipped with © MAT-LAB and © ControlDesk software), eight © PARKER SME60 motors and TPD-M drivers, a © dSPACE DS1007-based real-time controller and eight custom made winches. Each cable can exert a tension up to 150 N to the M P . The maximum velocity of each cable is equal to 5.9 m/s. The cable tensions are measured using eight FUTEK FSH04097 sensors, one for each cable, attached to cable anchor points. Their signal is amplified using eight FSH03863 voltage amplifiers and sent to the robot controller by a coaxial cable. CRAFT is equipped with a motion capture, OptiTrack that the sampling period is equal to 0.01 s.

Robot Modeling

Figure 1 depicts the main geometric parameters of CRAFT and the i th cable where i ∈ {1, . . . , m}, m being the number of cables attached to the n Degree-offreedom (DoF) M P (for CRAFT, m = 8 and n = 6). The dynamic equilibrium wrench equation can be expressed as:

Wτ -I p p -C ṗ + w g = 0 m (1)
with W being the robot wrench matrix expressed as:

W = u 1 . . . u i . . . u m b 1 × u 1 . . . b i × u i . . . b m × u m (2)
where u i is the unit vector of the i th cable vector pointing from cable anchor points B i to exit points A i , expressed in the reference frame F b . Vector b i pointing from point P to point B i and is expressed in the platform frame F p . From Eq. ( 1), τ = [τ 1 , . . . , τ i , . . . , τ m ] is the cable tension vector, I p is the inertia tensor, C is the Coriolis matrix and w g is the gravity wrench.

w g = m I 3 b R p ŜG g (3)
m is the M P mass, I 3 is the 3 × 3 identity matrix, b R p is the rotation matrix of the M P F p frame w.r.t the base frame F b , ŜG is the skew-symmetric matrix associated to s G = [x G , y G , z G ] ⊤ which represents the coordinate vector of the center of mass of the M P expressed in F p as:

ŜG =   0 -z G y G z G 0 -x G -y G x G 0   (4)
m b R p ŜG represents the first momentum of the M P defined with respect to base frame F b . g = [0, 0, -g] is the gravity acceleration vector where g = 9.81 m s -2 . The matrix I p in Eq. ( 1) can be expressed as:

I p = m I 3 -m b R p ŜG m b R p ŜG I p ( 5
)
with I p being the M P inertia tensor expressed in F p . The term C, which represents the centrifugal and Coriolis effects in Eq. ( 1) is expressed as:

C = 0 3 -m ω ŜG 0 3 ω I p ( 6
)
Using Eqs. ( 3), ( 5), and (6) in Eq. ( 1), gives:

m I 3 -m b R p ŜG m b R p ŜG I p ẗ α + 0 3 -m ω ŜG 0 3 ω I p ṫ ω -m I 3 b R p ŜG = Wτ (7)
As a first step of the identification, in this article, we will only consider platform translations, i.e., the platform does not rotate in relation to itself. Hence, the angular velocities and accelerations of the MP in Eq. ( 7) vanish and is reduced to:

m I 3 m b R p ŜG ẗ - mI 3 m b R p ŜG g = Wτ (8)
As the platform orientation is always constant, we can consider b R p = I 3 . Eq. ( 8) can be expressed as in the form as:

Ax = b (9)
where,

A =         ẗx 0 0 0 ẗy 0 0 0 ẗz + g 0 0 0 0 0 ẗz + g -ẗy 0 -ẗz -g 0 ẗx 0 ẗy -ẗx 0         , x =     m mx G my G mz G     , b = Wτ (10)
In Eq. ( 10), we need the linear acceleration of the M P ẗ and the cable tensions τ . As the robot is already equipped with the cable tension senors so we can directly acquire τ .

Identification Methodology

The identification methodology is based on the used of the inverse dynamic model (ID) Eq. ( 9), [START_REF] Gautier | On the identification of the inertial parameters of robots[END_REF], [START_REF] Gautier | Identification of standard dynamic parameters of robots with positive definite inertia matrix[END_REF], [START_REF] Ardiani | Contribution to the parametric identification of dynamic models: Application to collaborative robotics[END_REF]. This model is linear as function of the four parameters m, mx G , my G , and mz G . An off-line identification of these parameters is then considered, given measured or estimated off-line data for the linear acceleration of the M P ẗ and the cable tensions τ , collected while the robot CRAFT is tracking a planned trajectory. The model Eq. ( 9) is sampled at a sufficient number of time samples t i , for i = 1, • • • , n e , with (6n e ) ≫ 4, in order to get an overdetermined linear system of (6n e ) equations [START_REF] Ardiani | Contribution to the parametric identification of dynamic models: Application to collaborative robotics[END_REF]:

Y ID (τ ) = W ID x + ρ ( 11
)
where

Y ID (τ ) =    b 1 (6 × 1) . . . b ne (6 × 1)    ; W ID =    A 1 (6 × 4) . . . A ne (6 × 4)    , ( 12
)
ρ is the (6n e × 1) vector of errors between the data of the measurement of the torques in Y ID (6n e ×1) and the data W ID (6n e ×4) predicted by the model. These errors are due to noise measurement and modelling error. The identification problem consists in finding x the norm squared of the error ρ:

ρ 2 = Y ID -W ID x 2 (13)
The LS estimator models a process by fitting the parameter vector x according to the minimisation of Eq. ( 13). Then model [START_REF] Mishra | Kinematic stability based afg-rrt* path planning for cable-driven parallel robots[END_REF] becomes

Y ID (τ ) = W ID x + ρ (14)
The estimated parameter x is equal to:

x = W + ID Y ID ( 15
)
where W + ID is the pseudo-inverse of W ID . Since 6n e > 4 then

W + ID = (W ⊤ ID W ID ) -1 W ⊤ ID .
x is the unique LS solution of Eq. ( 15). The standard deviation σ x is estimated assuming that W ID is a deterministic matrix and ρ is a zero-mean additive independent Gaussian noise, with a covariance matrix C ρρ

C ρρ = E(ρρ ⊤ ) = σ 2 ρ I (16)
E is the expectation operator and I(6n e × 6n e ) is the identity matrix. By replacing Y ID in Eq. ( 15) with its definition Eq. ( 14) we obtain:

x = x + W + ID ρ ( 17
)
Applying the expectation operator at both sides of Eq. ( 17) ρ being a zero-mean additive independent Gaussian noise we get:

E( x) = E(x) ( 18
)
Thus the estimated values are unbiased. An unbiased estimation of the standard deviation σ ρ is:

σ 2 ρ = Y ID -W ID x 2 (6n e -4) (19)
W ID and being respectively a non-stochastic matrix and a stochastic matrix, the covariance matrix of the LS estimation error is then given by:

C x x = E(x-x ⊤ ) = (W ⊤ ID W ID ) -1 W ⊤ ID σ 2 ρ I W ID (W ⊤ ID W ID ) -1 = σ 2 ρ (W ⊤ ID W ID ) -1 (20) σ 2 x i = C x x (i, i) is the i th diagonal coefficient of C x x .
The relative standard deviation % σ x i for each identified parameters is given by:

% σ x i = σ x i | x i | (21)
The calculation of Eq. ( 14) and the condition number of W ID can be obtained using the singular value decomposition (SVD) of W ID .

Semi-implicit Homogeneous differentiator

The purpose is to estimate the velocity of the angular variable exclusively from the measured position of the output shaft for each of the eight motors. To estimate the velocity let us introduce the continuous-time state:

Σ :    ẋ1 = x 2 ẋ2 = p(t) y = x 1 (22)
where x 1 and x 2 are respectively the angular variable and its velocity; y is the measure of x 1 with additional noise η. Here p(t) is a bounded perturbation, which is unknown and assumed to be a constant parameter or a slowly variable. This implies that for a sufficient small sampling-time h > 0, p ≡ p + , with the notation for a discretized variable:

•(t = (k + 1)h) = • + and •(t = kh) = •.
Homogeneity approach is interesting because if a local stability is obtained due to the dilatation, this framework allows extending this local property to global settings, [START_REF] Rosier | Homogeneous lyapunov function for homogeneous continuous vector field[END_REF]. This differentiator can be written such as [START_REF] Perruquetti | Finite-time observers: application to secure communication[END_REF], [START_REF] Ghanes | A new varyinggain-exponent-based differentiator/observer: An efficient balance between linear and sliding-mode algorithms[END_REF],

   ż1 = z 2 + λ 1 µ⌈ǫ 1 ⌋ α ż2 = λ 2 µ 2 ⌈ǫ 1 ⌋ 2α-1 ŷ = z 1 (23)
where α ∈ ]0.5 1[ has to be fixed [START_REF] Hong | On an output feedback finite-time stabilization problem[END_REF], ǫ 1 = yz 1 , and the notation ⌈•⌋ α = | • | α sgn(•) is adopted along the paper. The degree of homogeneity of the differentiator (23) d is equal to α -1 with respect to dilatation Λ r with r = (r 1 = 1, r 2 = 1) [START_REF] Perruquetti | Finite-time observers: application to secure communication[END_REF]. Moreover, λ i > 0, i = 1, 2 are the linear part gains, which are considered and allow to have the eigenvalues of the differentiation error ǫ 1 sufficiently stables, while the coefficient µ is chosen sufficiently large to cancel the effect of the unknown perturbation p(t).

The discretization of the differentiator ( 23) is based on the so-called implicit projection that acts as corrective terms and aim to "generalize" the sign function in sliding-based differentiator in order to reduce the chattering and preserve stability properties for high time steps). Two projectors N 1 and N 2 are used respectively to design the correction terms with the differentiator SIHD-2. In [START_REF] Michel | Experimental validation of two semi-implicit homogeneous discretized differentiators on the CRAFT cable-driven parallel robot[END_REF], it has been highlighted that the SIHD-2 algorithm offers better performances than with only one projector since the projectors N 1 and N 2 are respectively dedicated to the estimation of z 1 and z 2 . Considering the signal to differentiate y, the error ǫ 1 = yz + 1 , and the definition of the notation ⌈ǫ + 1 ⌋ ≡ |ǫ 1 |sign(ǫ + 1 ), As a result, the semi-implicit homogeneous Euler discretization based on two projectors (SIHD-2) reads as:

z + 1 = z 1 + h z + 2 + λ 1 µ|ǫ 1 | α d N 1 z + 2 = z 2 + E + 1 hλ 2 µ 2 |ǫ 1 | 2 α d -1 N 2 . ( 24
)
The projector N 1 , and N 2 are defined:

N 1 (ǫ 1 ) :=    ǫ 1 ∈ SD → N 1 = ⌈ǫ 1 ⌋ 1-α d λ 1 µh , E + 1 = 1 ǫ 1 / ∈ SD → N 1 = sign(ǫ 1 ), E + 1 = 0, (25)
with the domain of attraction

SD = {ǫ 1 / |ǫ 1 | ≤ (λ 1 µh) 1 1-α }. N 2 (ǫ 1 ) :=    ǫ 1 ∈ SD ′ → N 2 = ⌈ǫ 1 ⌋ 2(1-α d ) λ 2 h 2 µ 2 ǫ 1 / ∈ SD ′ → N 2 = sign(ǫ 1 ), (26)
where

SD ′ = {ǫ 1 ∈ SD/ |ǫ 1 | ≤ (λ 1 µ 2 h 2 ) 1 2(1-α d ) ≡ |ǫ 2 | ≤ (λ 1 µ 2 ) 1 2(1-α d ) h α d 1-α d }.
The filtering properties of the differentiator avoid using an extra filter before the differentiation, and the cascade offers the possibility of adjusting the homogeneous parameter α d with some flexibility. The measured angular positions are noisy such as y becomes y m = x 1 + η.The λ i , i = 1, 2 parameters are chosen such as the linear part is stable. The value of homogeneous exponent α d is chosen to allow better filtering properties of the estimated differentiation. The µ parameter is also chosen by numerical test trial and error in order to determine the best possible action of both projectors N 1 and N 2 . The numerical values of these five parameters are tuned as follows:

λ 1 = 210, λ 2 = 210, λ 3 = 525, λ 4 = 525 α d = 0.95, µ = 1. ( 27
)

Experiments

The platform performs four orbital movements. A video of these movements is available at3 . A sampling period of 10 ms is used to measure the position vector t and the cable tension τ with the FUTEK FSH04097 sensors. The eight measured cable tensions are plotted Fig. 2 a). The Euler angles of the platform during its motion are shown in Fig. 2 b). This motion of the platform is indeed translational as it is predicted by the trajectory planning. The identification of the parameters m, mx G , my G , and mz G needs the linear acceleration vector of the M P ẗ. For the best knowledge of this acceleration vector a comparison is made from the measured platform position t between the estimation of the velocity vector ṫ and the acceleration vector ẗ respectively thanks to the proposed SIHD-2 algorithm and the estimation of these vectors thanks the usual Euler discretization method. With the usual Euler discretization method the numerical tests shows that the estimation of ṫ and ẗ requires filtering the position vector t of the platform with a cutoff frequency of 4Hz. The filter applies a zero-phase forward and reverse digital infinite impulse response (IIR) filtering. Figure 3 shows the estimated velocity and acceleration vectors in 3D by using the SIHD-2 algorithm and the usual Euler discretization method. The curves are similar between the two algorithms. Table 1 gathers the four identified parameters. The standard deviations are not computed based on the standard deviation from the actual parameter values that are unknown. The numerical values obtained by the two methods for each identified parameter are very close to each other. However, having to filter the position measurement of the platform is a strong constraint in the perspective of doing real time identification of the dynamic parameters of CRAFT. That is why the SIDH -2 algorithm, which is iterative, is a relevant solution in the perspective of real time identification.

Discussion and Future Work

This research paper presents a novel methodology for identifying the dynamic parameters of Cable-Driven Parallel Robots (CDPRs). Initially, the focus is on

i σ x i % σ x i x i σ x i % σ x i m (kG)
9.88 0.0027 0.027 9.88 0.0027 0.027 mx G (kGm) -0.064 0.0027 4.21 -0.064 0.0027 4.16 my G (kGm) 0.019 0.0027 13.86 0.019 0.0027 13.65 mz G (kGm) -1.059 0.070 6.67 -1.210 0.078 6.41 identifying four parameters, specifically the platform mass and the location of its center of mass, using an orbital trajectory that has been specifically designed to achieve pure platform translation movement in all the three xyz-plane. The trajectory, which occurs within a single plane, allows for the identification of these four dynamic parameters. The values determined by the identification of the inertial parameters are consistent with the actual platform mass (measured) and the estimated center of mass location from the computer aided Design. The results have also been compared with other types of trajectories, such as pointto-point linear trajectories. The reason for choosing a circular trajectory is the availability of continuous platform accelerations which facilitates the use of the proposed methodology. The results obtained from experimental measurements of the platform position, as well as measurements treated with the application of a homogeneous semi-implicit differentiator, indicates that despite the complexity of the CRAFT, these results allow for the real-time identification of all the essential dynamic parameters.

In future work, the proposed approach will be used to identify the remaining six dynamic terms associated with the inertia tensor, specifically three terms associated with the platform moment of inertia and other three are products of inertia. To manage the future identification tasks, an IMU sensor will be used in order to provide acceleration measurements that are synchronized with those of the cable tensions. Additionally, we will also aim to identify the friction that acts throughout the actuators of CRAFT.

Fig. 1 :

 1 Fig. 1: (left) CRAFT prototype located at LS2N, Nantes, France and (right) CRAFT Geometric Parametrization.

  ṗ = [ ṫ, ω] ⊤ and p = [ ẗ, α] ⊤ are the vectors of the M P velocity and acceleration, where, ṫ = [ ṫx , ṫy , ṫz ] ⊤ and ẗ = [ ẗx , ẗy , ẗz ] ⊤ are the vectors of the M P linear velocity and acceleration, while ω = [ω x , ω y , ω z ] T and α = [α x , α y , α z ] ⊤ are the vectors of the M P angular velocity and angular acceleration, respectively. The gravity wrench is:

Fig. 2 :

 2 Fig. 2: a) Cable tensions τ as function of time, b) Orientation of the mobile platform.

Fig. 3 :

 3 Fig. 3: ṫ and ẗ along the orbital trajectory of the platform: velocity (a) and acceleration (b) with SIDH-2; velocity (c) and acceleration (d) with Euler discretization.

Table 1 :

 1 Numerical values of the identified parameters

	Identified	differentiation	Euler
	parameters	SIHD-2	discretization
		x

https://youtu.be/I-IOcAGha3o

⋆ Supported by the EquipEx+ TIRREX project, grant ANR-21-ESRE-0015 https://www.ls2n.fr/le-projet-tirrex-est-laureat-de-lami-equipex/


G Radtke

A Saúl

Y Klein

G Rousse

Magnetism of Ba 4 Ru 3 O 10 revealed by density functional calculations : structural trimers behaving as coupled magnetic dimers

Keywords: numbers: 75.47.Lx, 75.25.-j, 71.27.+a

Density functional calculations have been carried out to investigate the origin of the peculiar magnetic properties of Ba4Ru3O10. These calculations suggest that, under the influence of Ru-Ru covalent bonding, structural trimers behave in this compound as strong S = 1 antiferromagnetic dimers ordering below TN under the influence of weaker inter-trimer couplings.

I. INTRODUCTION

In the last 20 years, ruthenium oxides have been the topic of numerous studies because they often exhibit a strong relationship between structural, orbital and spin degrees of freedom leading to unusual physical properties 1 . For example, this interplay is responsible for the spin-singlet dimerization observed in La 4 Ru 2 O 10 2-6 and for the spontaneous formation of spin (S = 1) Haldane chains in Tl 2 Ru 2 O 7 7 . Another illustration of the complexity of Ru-oxides is the Ca 2-x Sr x RuO 4 prototype, in which one may considerably tune the transport and magnetic properties by an isovalent substitution, the size effect being the only driving force 8 . While on the Sr-rich side, unconventional superconductivity is found with the possibility of spin-triplet pairing 9 , Ca 2 RuO 4 shows an orbital selective Mott transition near room temperature 10 and was proposed to be used as a temperature sensor for bolometers 11 .

Recently, we have shown that the antiferromagnetic (AFM) phase transition at T N ∼ 105 K in Ba 4 Ru 3 O 10 is accompanied by drastic variations of the transport properties 12 . This compound is built from Ru 3 O 12 trimers of face-shared Ru(2)O 6 -Ru(1)O 6 -Ru(2)O 6 octahedra (Fig. 1(a) and (b)). Each trimer is connected to four neighbouring trimers through corner sharing of Ru(2)O 6 octahedra to form corrugated layers that are then stacked together along the b axis. This arrangement leads to the formation of zig-zag chains of Ru(2)O 6 octahedra running along the a axis in each layer (Fig. 1(c)). While the high-temperature susceptibility was interpreted in term of paramagnetic (PM) Ru 4+ in the lowspin (S = 1) configuration, low-and intermediatetemperature magnetic properties asked for more clarifications. In particular, neutron powder diffraction patterns below T N indicate that the Néel order is of a new kind since outer Ru (2) order in an AFM lattice while inner Ru (1) do not exhibit any ordered magnetic moment.

In this paper, based on Density Functional Theory (DFT) calculations, we investigate the electronic struc-ture of Ba 4 Ru 3 O 10 and identify the metal-metal covalent bonding within the trimers as a fundamental parameter influencing the magnetic behaviour of this compound. In particular, we show how structural trimers might be considered as strong AFM dimers ordering below T N under the influence of weaker inter-trimer couplings.

II. COMPUTATIONAL DETAILS

Electronic structure and total energy calculations were performed using the Wien2k code [START_REF] Blaha | Wien2k, An Augmented Plane Wave Plus Local Orbitals Program for Calculating Crystal Properties[END_REF] . This code is an implementation of the full-potential linearized augmented plane-wave method based on density-functional theory. The calculations presented hereafter were performed using the generalized gradient approximation of Perdew, Burke, and Ernzerhof (GGA) [START_REF] Perdew | [END_REF] for exchange and correlation. Muffin-tin spheres radii were set to 2.43 a.u. for Ba, 1.91 a.u. for Ru, and 1.69 a.u. for O in all of these calculations. When included, on-site Hubbard U term was taken into account using the method proposed by Anisimov et al. 15 with U eff = 2.5 eV (J = 0 eV) typically employed for Ru-4d states 4 and Spin-orbit coupling (SO) was treated using scalar-relativistic wavefunctions in a second-variational procedure 16 . All these calculations have been performed using the low-temperature (determined at 10 K) experimental crystal structure 12 .

III. RESULTS

A. Electronic structure Ba 4 Ru 3 O 10 is described in the orthorhombic Cmca space-group, with lattice parameters a = 5.767 Å, b = 13.245 Å and c = 13.064 Å and containing four formula units per cell. Ru atons are distributed among two Wyckoff sites, 4a at (0,0,0) corresponding to the Ru(1) in the middle of the trimer and 8f at (0,y,z) with y ≈ 0.87 and z ≈ 0.14 corresponding to the Ru(2) at (1) and Ru(2)-4d paramagnetic LDOS decomposed on the real t 0 2g , t + 2g , t - 2g , e + g , e - g orbitals and with the z axis along the metal-metal bond direction within the trimer.

the edges. The basic features of the electronic structure of Ba 4 Ru 3 O 10 are shown in the total and local densities of states (LDOS) obtained from a non-magnetic calculation and presented in Fig. 2(a). The valence band is essentially built from O-2p states located between -7 and -2 eV whereas Ru-4d states are found at higher energies, mainly between -2 and +4 eV. A metallic solution is found, consistent with a filling of two third of the Ru-4d t 2g bands for both Ru sites as expected for Ru 4+ . A closer inspection of the Ru(1) and Ru(2) LDOS (Fig. 2(b)) however, reveals a first interesting feature of this compound : the LDOS of Ru(2) at the Fermi level (N (E f ) ≈ 3.7 eV -1 atom -1 ) is almost twice as large as that of Ru(1) (N (E f ) ≈ 1.8 eV -1 atom -1 ). Although not strictly applicable to these LDOS, Stoner criterion might still appear here as a useful guide predicting a qualitatively different behaviour of these two ions towards ferromagnetic instability, as already noted by Streltsov and Khomskii 17 . One might indeed expect a weaker energy lowering through spin-polarization for Ru(1) than for Ru(2). This is precisely what a (GGA) ferromagnetic calculation reveals : a more stable solution is found through spin polarization (by about 115 meV/f.u.) but a larger magnetic moment is found on Ru(2) ions (∼ 0.94 µ B ) than on Ru(1) ions (∼ 0.34 µ B ). To understand the mechanism leading to this peculiar imbalance in the behaviour of the two inequivalent Ru ions, it is instructive to analyse the symmetry of the Ru-4d states close to the Fermi level. To account for the particular face-shared geometry of the coordination octahedra within the trimer, Ru-4d LDOS have been decomposed on a set of alternative real orbitals obtained for a quantization z axis taken along the threefold axis of the perfect octahedron. The expression of these orbitals in the basis of the conven- tional cubic orbitals is given by 18 2) sites is therefore not O h but we keep the t 2g and e g notations here for convenience), this choice proves particularly useful as it brings the local z axis along the metalmetal direction. t 0 2g orbitals have therefore their main lobes pointing through the face shared by Ru(1)O 6 and Ru(2)O 6 coordination octahedra, t + 2g and t - 2g lobes point toward the middle of the edges shared by these octahedra whereas e + g and e - g orbitals point directly toward the O ions. LDOS for Ru (1) and Ru(2), shown respectively in Fig. 2(c) and (d), display large crystal-field splittings between low-lying t 2g and high-lying e g states. A striking feature however occurs in the t 2g states : whereas t + 2g and t - 2g states are found around the Fermi level and exhibit a moderate energy dispersion, the t 0 2g states appear for Ru(2) in three regions respectively located below, at and above the Fermi level separated from each other by about 1 eV. The t 0 2g LDOS for Ru(1) echoes that of Ru(2) with however the absence of states at the Fermi level. This pattern is directly related to the large σ bonding covalent interaction experienced by these orbitals. Molecular orbitals arising within the trimer can be pictured using a simple three-site model accounting for hopping t < 0 and overlap S > 0 integrals (the Ru-Ru distance of ∼ 2.55 Å here is much shorter than in metallic Ru) between nearest-neighbours. Three bonding, nonbonding and anti-bonding molecular orbitals are formed, of respective energies

|t 0 2g = |d z 2 |t + 2g = 2 3 |d x 2 -y 2 - 1 3 |d xz |t - 2g = 2 3 |d xy + 1 3 |d yz and |e + g = 1 3 |d x 2 -y 2 + 2 3 |d xz |e - g = 1 3 |d xy - 2
b = √ 2t/(1 + √ 2S), nb = 0 and ab = - √ 2t/(1 - √ 2S
). The expression of these molecular orbitals in term of Ru-t 0 2g states for a trimer |φ (2) , |φ (1) and |φ (2) are respectively :

|ψ b = 1 4 1 + √ 2S |φ (2) + √ 2|φ (1) + |φ (2) |ψ nb = 1 √ 2 |φ (2) -|φ (2) |ψ ab = 1 4 1 - √ 2S |φ (2) - √ 2|φ 1 + |φ (2)
The non-bonding molecular orbital is therefore built from Ru(2)-t 0 2g states only and explains, at least partially, the reduced density of states at the Fermi level observed for Ru(1) as seen in Fig. 2(b). These results are consistent with the work published recently by Streltsov and Khomskii 17 .

B. Magnetic order and couplings

To determine the ground-state magnetic order of this compound, the total energies of 30 distinct spin configurations compatible with a single orthorhombic conventional cell of Ba 4 Ru 3 O 10 have been computed in GGA and GGA+U. In both cases, the ground-state configuration is the same and corresponds to an AFM arrangement Optimized Ising model (K) DFT relative energy (K) of the Ru(2) moments within the trimers and between adjacent trimers of the same layer, i.e. on the zig-zag chains running along a. The staking of these AFM layers along b is ferromagnetic (FM). It should be added that in this ground-state AFM structure, Ru(1) ions are nonmagnetic in both GGA and GGA+U. A summary of the magnetic moments calculated in the muffin-tin spheres for both a FM and the ground state AFM spin arrangements is presented in Tab. I. These results are in perfect agreement with the magnetic structure deduced from neutron diffraction, which indicates that Ru(1) moments are not ordered whereas Ru(2) moments order with a propagation vector k = (0, 0, 0). The eight Ru(2) ions on the 8f site exhibit magnetic moments whose signs follow the (+, +, -, -, +, +, -, -) sequence 12 . Total DOS and 4d-LDOS of Ru(1) and one of the two Ru(2) of the trimer calculated in GGA+U for this ground-state configuration are shown in Fig. 3(a),(b) and (c) respectively. By symmetry, the magnetic moment of Ru(1) cancels through a double occupation of both t + 2g and t - 2g orbitals. Ru(2) ions, on the contrary, are strongly spin-polarized and adopt an electronic configuration close to the ionic low-spin S = 1 configuration with a single occupation of both t + 2g and t - 2g orbitals and a weaker polarization of the t 0 2g states. The resulting magnetic moment on Ru(2) is 1.22 µ B , largely reduced by covalent interactions from the expected value for a low-spin configuration (2 µ B ). It should be added here that the polarization of Ru(1) ions depends on the spin configuration under study but never exceeds 0.11 µ B in GGA+U. The next important step of this work consists in the evaluation of the couplings occurring in Ba 4 Ru 3 O 10 between magnetic Ru(2) ions. In the broken symmetry formalism, the evaluation of these couplings is based on the mapping of total energies corresponding to various spin arrangements within a supercell onto an Heisenberg Hamiltonian

J chain = 86 K J trimer = 344 K J plane = 25 K +139 K +289 K +469 K +599 K +670 K +745 K +932 K Ground-state (b) (a) GGA GGA GGA+U +429 K +889 K +1411 K +1802 K +1173 K +1354 K +1802 K Ground-state (
Ĥ = Ĥ0 + i>j J ij Ŝi Ŝj ( 1
)
where Ĥ0 is the spin-independent part of the Hamiltonian, J ij are the magnetic couplings to determine and Ŝi and Ŝj are, in our case, the S = 1 spin operators localized on Ru(2) sites i and j respectively. It is straightforward to show that the expectation value of Hamiltonian (1) on a DFT state |α (obtained by preparing the initial electron density according to a particular collinear spin arrangement in the supercell and performing a selfconsistent calculation until convergence) can be simply written under the form of an Ising Hamiltonian 19

DF T α = α| Ĥ|α = 0 + S 2 i>j J ij σ i σ j (2)
with σ i = ±1 and here, S = 1. In strongly localized systems, such as 3d transition metal oxides, Eq. ( 2) can be employed to model large sets of spin configurations and a least-square minimization of the difference between DFT and Ising relative energies, can be applied to obtain a numerical evaluation of the couplings 20,21 .

In our present case however, this procedure fails to fit accurately the energies of the 30 spin configurations, i.e. no single set of magnetic couplings can be used in Eq. ( 2) to describe correctly their relative energies. This procedure is however well adapted to describe a subset of low-energy configurations found right above the ground-state AFM order. In both cases, the same subset of 8 configurations (shown in blue in Fig. 4) including the AFM ground-state has been used. Three magnetic couplings, shown in Fig. 1, have been employed to obtain these results : the intra-trimer coupling J trimer

Ru(2)

Ru( 1) Ru(2) t2g +/- FIG. 5: (color online) Schematic of the kinetic superexchange mechanism proposed to explain the AFM coupling within the trimers.

(5.10 Å), the inter-trimer coupling J chain (3.92 Å) occurring between trimers of the same corrugated planes, i.e. between Ru(2) ions forming the zig-zag chains running along a and finally J plane (5.87 Å) between nearestneighbour Ru(2) ions of adjacent planes. It should be noted here that the use of these 30 spin configurations allows for the calculation of additional couplings between the planes but that their inclusion in the fitting procedure leads to negligible values and does not modify the values of the dominant couplings nor improves the quality of the fit.

The resulting magnetic couplings are summarized in Tab. II. Several conclusions can be drawn from these results. Firstly, both approximations lead to a consistent overall picture of the magnetic interactions in this compound where a dominant AFM interaction occurs within the structural trimers (J trimer ) coupled by weaker AFM interactions either in the corrugated planes (J chain ) or between these planes (J plane ). As shown in Fig. 5, the dominant AFM coupling J trimer can be understood by the particular occupation scheme of the t + 2g and t - 2g orbitals in the Ru(2)-Ru(1)-Ru(2) trimer supporting the hypothesis of an AFM interaction between moments localized on the Ru(2) ions through a kinetic superexchange 23 mechanism mediated by the Ru(1) orbitals or, more precisely, by Ru(1)-O hybrids.

Secondly, inclusion of the Hubbard U term in the calculations, employed here to improve the description of electron correlation occurring within the Ru-4d states, leads to a reduction of the amplitude of these AFM couplings. This effect, also observed in other Ru-based systems 22 , is consistent with the expected U-dependence of the AFM contribution to the magnetic couplings derived from the one-band Hubbard model J AF M ≈ 4t 2 /U where t and U are respectively the hopping and on-site Coulomb repulsion parameters.

Thirdly, as mentioned above and as it can be clearly seen in Fig. 4(a) and (c), large discrepancies are found between DFT and Ising relative energies for most of the high-lying spin configurations. To explain this behaviour, one has to consider that the spin Hamiltonian described by Eq. (1) relies on the assumption that the spatial part of the electronic wavefunctions does not depend on the particular spin arrangement employed in the supercell. If this assumption is relatively well verified in strongly localized 3d systems, it is certainly much more questionable in the case of more extended Ru-4d orbitals and in particular in Ba 4 Ru 3 O 10 where strong covalent metal-metal bonding occurs. One might therefore expect that only a (low-energy) part of the wide range of magnetic excitations covered by our set of spin configurations verifies approximately this assumption.

C. Magnetic susceptibility

These results led us to reinterpret the magnetic susceptibility of this compound above T N as that of an ensemble of coupled S = 1 spin AFM dimers, i.e. to consider a structural trimer as a magnetic dimer where S = 1 spins are hold by the Ru(2) ions and where the Ru(1) is non-magnetic. Modelling the magnetic excitations of independent S = 1 spin dimers with an Heisenberg Hamiltonian characterized by the intra-dimer coupling J trimer , one can easily derive the expression of the molar magnetic susceptibility for this system as

χ dim = 2βg 2 µ 2 B N A [1 + 5 exp(-2βJ trimer )] [3 + exp(βJ trimer ) + 5 exp(-2βJ trimer )] (3)
where β = 1/k B T , N A is the Avogadro number, µ B is the Bohr magneton and g = 2. Treating the inter-dimer interactions through a mean-field approximation, one further obtains 24,25

χ = χ dim (1 + λχ dim ) + χ 0 (4)
where λ = J inter /N A (gµ B ) 2 , J inter ≈ 2J chain + 2J plane , χ 0 is a temperature-independent contribution and χ dim is the susceptibility of S = 1 dimers as defined in Eq. (3). Magnetic susceptibility has been measured within a Physical Property Measurement System (PPMS) from Quantum Design equipped with a vibrating sample magnetometer. Data have been collected under a magnetic field of 0.1 T in the field-cooled mode. The fit of this experimental magnetic susceptibility for the 110-350 K range is shown in Fig. 6. This fit was obtained for J trimer = 171 K, J inter = 114 K and χ 0 = 2.10 -4 emu/mol. Both intra-and inter-trimer couplings are AFM with a dominant intra-trimer coupling. The picture drawn from this analysis is therefore compatible with the one deduced from DFT of S = 1 spins hold by the Ru(2) ions only coupled by a strong AFM interaction within the structural trimers and interacting with each other through weaker AFM inter-trimer couplings. It should be mentioned here that this fit is also compatible with the intermediate-and high-temperature measurements from T N up to ∼ 612 K performed by Dussarrat et al. 26 .

D. Magnetocrystalline anisotropy

To complete our analysis and lift the ambiguity on the orientation of the magnetic moments left by neutron diffraction measurements 12 , GGA+U+SO calculations have been carried out for the ground-state AFM order within the framework of collinear magnetism. This approximation is mainly justified by the fact that only the (002) Bragg reflection presents an extra magnetic contribution below T N in neutron diffraction 12 , limiting considerably the data available to refine potential noncollinear structures. The DFT total energy dependence upon the Ru(2) magnetization orientation has been systematically calculated in the ab and bc planes, as shown in Fig. 7(a). These results clearly point out a as the magnetization direction minimizing the total energy. These results can be understood qualitatively by considering the crystal-field levels of Ru (2) shown in Fig. 3(c) and treating the SO as a local perturbation. The SO Hamiltonian Ĥso = ξ L • Ŝ can be written, for an orientation of the magnetization along n(θ, φ) (θ and φ are respectively the polar and azimuthal angles in the local coordinate system shown in Fig. 1 matrix elements, both proportional to -(iξ sin θ)/2, the SO stabilization is maximized for θ = π/2 and thus fixes the magnetization in the plane perpendicular to the local z axis, i.e. perpendicular to the Ru-Ru bond. The sole magnetization direction able to satisfy this condition in Ba 4 Ru 3 O 10 within the approximation of collinear magnetism is therefore a. This result is consistent with the neutron powder diffraction as shown in Fig. 7(b) through the Rietveld refinement of the magnetic structure using the AFM ground-state order with a magnetization laying along a and a magnetic moment of 1.11 (11) µ B on the Ru(2) ions only (to be compared to 1.22 µ B from GGA+U calculations, see Tab. I).

IV. SUMMARY

In conclusion, this work emphasizes the dominant role played by Ru-Ru covalent interactions on the magnetism of Ba 4 Ru 3 O 10 along with crystal-field effects, electron correlation and spin-orbit coupling. In particular, it provides a new picture of this compound where structural trimers behave as magnetic dimers of S = 1 spins hold by the Ru(2) ions in the low-spin electronic configuration, the Ru(1) ions being essentially non-magnetic. The phase transition observed at T N might therefore be considered as a transition from a state dominated by short-range intra-trimer AFM interactions to a Néel order mediated by weaker inter-trimer AFM interactions.

FIG. 1 :

 1 FIG. 1: (color online) (a) Atomic structure of Ba4Ru3O10. (b) Local system of coordinates employed to decompose Ru-4d LDOS. (c) Structure of a corrugated layer observed along the b axis. Dominant magnetic couplings occurring between Ru(2) ions are indicated in (a) and (c) and circles highlight the structural trimers.

FIG. 2 :

 2 FIG. 2: (color online) (a) Non-magnetic DOS of Ba4Ru3O10 calculated in GGA (b) Detail of Ru(1) and Ru(2) paramagnetic LDOS. (c) and (d) Ru(1) and Ru(2)-4d paramagnetic LDOS decomposed on the real t 0 2g , t + 2g , t - 2g , e + g , e - g orbitals and with the z axis along the metal-metal bond direction within the trimer.

3

 3 |d yz and the local systems of coordinates for both Ru(1) and Ru(2) sites are shown in Fig. 1(b). Although the particular structure of Ba 4 Ru 3 O 10 allows for local distortions of the octahedra by breaking the third order axis of the trimer (the point group symmetry of the Ru(1) and Ru(

FIG. 3 :

 3 FIG. 3: (color online) (a) DOS in Ba4Ru3O10 calculated in GGA+U for the ground-state AFM order. (b) Corresponding 4d-LDOS of Ru(1) and (c) of one the Ru(2) of a trimer.

FIG. 4 :

 4 FIG. 4: (color online) (a) Results of the least-square fitting procedure employed to calculate magnetic couplings in GGA: only the low-energy configurations shown in blue have been used, the configurations shown in green have been discarded. The energy of the ground-state AFM configuration is used as the reference energy. (b) Detailed presentation of the 8 low-energy configurations considered for the evaluation of magnetic couplings. (c) and (d) show the results of the same procedure applied to GGA+U. Red (blue) Ru(2)O6 coordination octahedra correspond to up (down) arrangements of the Ru(2) magnetic moments. According to the convention used in Eq. (1), AFM couplings are positive.

FIG. 6 :

 6 FIG. 6: (color online) Magnetic susceptibility of Ba4Ru3O10 recorded up to 350 K fitted using Eq.(3).

FIG. 7 :

 7 FIG. 7: (color online) (a) Dependence of energy on the direction of the Ru(2) magnetic moments in the ab and bc planes calculated for the ground-state AFM order in GGA+U+SO. (b) Rietveld refinement of the neutron powder diffraction pattern recorded at 1.8 K 12 . Red circles, black and green lines represent the observed, calculated and difference patterns, respectively. The expected positions of the Bragg reflections are shown as vertical bars: the upper series correspond to the nuclear phase, while the lower series marks the magnetic contribution. The inset shows an enlargement of the (020)/(002) doublet (marked by an arrow), with the blue line corresponding to magnetic contribution to the total intensity.

2 L-e iφ sin θ + 1 2 Ŝ+- 2 + L-e iφ cos 2 θ 2

 2222 (b)), as follows27 :Ĥso = ξ Ŝn Lz cos θ + 1 2 L+ e -iφ sin θ + 1 Lz sin θ -L+ e -iφ sin 2 θ Lz sin θ + L+ e -iφ cos 2 θ 2 -Le iφ sin 2 θ2Assuming, as it could be expected from Fig.3(c), that leading terms in the second-order non-degenerate perturbation expansion imply t +↑ 2g | Ĥso |t -↓ 2g and t -↑ 2g | Ĥso |t +↓ 2g

TABLE I :

 I Magnetic moments (in µB) in the Ru muffin-tin spheres calculated for a FM and ground-state AFM orders.

		GGA		GGA+U
		Ru(1)	Ru(2)	Ru(1)	Ru(2)
	FM	0.34	0.94	0.09	1.12
	AFM	0.0	0.88	0.0	1.22

TABLE II :

 II Magnetic couplings between Ru(2) ions treated as S = 1 spins and extracted from the optimization of Eq. (2) for the 8 low-energy configurations shown in Fig.4(b) and (d).According to the convention used in Eq. (2), AFM couplings are positive.

		Jtrimer (K)	J chain (K)	J plane (K)
	GGA	344	86	25
	GGA+U	265	30	10

Acknowledgements

Y.K. thanks the staff of the Low Temperature Physics plateforme (UPMC).