Eric Marchand
email: eric.marchand@irisa.fr.

Subspace-based Direct Visual Servoing

Keywords:

To date most of visual servoing approaches have relied on the geometric features that have to be tracked and matched in the image. Recent works have highlighted the importance of taking into account the photometric information of the entire images. This leads to direct visual servoing (DVS) approaches. The main disadvantage of DVS is its small convergence domain compared to conventional techniques, which is due to the high non-linearities of the cost function to be minimized. In this paper we propose to project the image on an orthogonal basis (PCA) and then servo on either images reconstructed from this new compact set of coordinates or directly on these coordinates used as visual features. In both cases we exhibit the analytical formulation of the interaction matrix. We show that these approaches feature a better behavior that the classical photometric visual servoing scheme allowing larger displacements and a satisfactory decrease of the error norm thanks to a well modelled interaction matrix.

I. INTRODUCTION

V ISUAL servoing uses the information provided by a vision sensor to control the movements of a dynamic system [START_REF] Chaumette | Visual servo control, Part I: Basic approaches[END_REF]. This approach requires the extraction of visual information (usually geometric features) from the image in order to design the control law. Robust extraction and realtime spatio-temporal tracking of these visual cues is to date a non trivial task.

While there has been progress in extracting and tracking relevant features, a new approach called direct visual servoing (DVS) has been emerging for almost 10 years now [START_REF] Kallem | Kernel-based visual servoing[END_REF], [START_REF] Collewet | Visual servoing set free from image processing[END_REF], [START_REF] Collewet | Photometric visual servoing[END_REF], [START_REF] Dame | Mutual information-based visual servoing[END_REF]. It has been demonstrated that only the luminous intensities of the images can be taken into account to control the robot and that conventional tracking processes can be avoided. The main drawback of DVS is its convergence domain compared to classical VS techniques, which is due to the high non-linearities of the cost function to be minimized. Various schemes have been proposed in order to improve the robustness of DVS by considering various descriptors (image intensity, gradient, color, etc.) or cost functions (mutual information [START_REF] Dame | Mutual information-based visual servoing[END_REF], histogram distances [START_REF] Bateux | Histograms-based visual servoing[END_REF], mixture of Gaussians [START_REF] Crombez | Photometric Gaussian mixtures based visual servoing[END_REF]). Recently, it has been proposed to consider convolutional neural network to bypass the modelling step [START_REF] Bateux | Training deep neural networks for visual servoing[END_REF].

Another solution to increase the convergence domain would be extract from the image a set of coefficients that could then be used as control input. The idea is not to extract geometric features from the image (as it is usually done in visual servoing) but to "compress" the original image information in order to get a compact representation. This is what has been done with photometric moments which allows retain geometric information [START_REF] Bakthavatchalam | A Direct Dense Visual Servoing Approach using Photometric Moments[END_REF]. It was shown that it provides a better behavior that a classical control based on points [START_REF] Chaumette | Visual servo control, Part I: Basic approaches[END_REF] and extends significantly the convergence of the photometric visual servoing approach [START_REF] Collewet | Photometric visual servoing[END_REF]. Another interesting approach was proposed in order to consider Wavelet [START_REF] Ourak | Direct visual servoing using wavelet coefficients[END_REF], [START_REF] Duflot | Wavelet and Shearlet-based Image Representations for Visual Servoing[END_REF] and Shearlet-based [START_REF] Duflot | Wavelet and Shearlet-based Image Representations for Visual Servoing[END_REF] image representations and thus to consider in the control law the wavelet/shearlet coefficients. The authors [START_REF] Ourak | Direct visual servoing using wavelet coefficients[END_REF], [START_REF] Duflot | Wavelet and Shearlet-based Image Representations for Visual Servoing[END_REF] then proposed an analytical formulation of the interaction matrix that links the variation of these coefficients to the camera motion.

Actually, one of the very first attempt to achieved DVS has been proposed in [START_REF] Nayar | Subspace methods for robot vision[END_REF], [START_REF] Deguchi | Visual servoing using eigenspace method and dynamic calculation of interaction matrices[END_REF], [START_REF] Deguchi | A direct interpretation of dynamic images with camera and object motions for vision guided robot control[END_REF]. In these papers, image intensity is not directly considered but an eigenspace decomposition is performed to reduce the dimensionality of image data. This is done thanks to a Principal Component Analysis (PCA) process (also known as Karhunen-Loève expansion). The control is then performed on the image coordinates in the eigenspace (an orthogonal basis). This process requires the off-line computation of this eigenspace and then, for each new frame, the projection of the image on this subspace in order to compute the set of coordinates (coefficients) in the new basis that will be used in the control law. An interest of such approach is that, when projecting an incoming image on this basis, the greatest variance comes to lie on the first coefficient, the second greatest variance lies on the second coefficient, etc. Only a few coefficients allow to grasp most of the variance of the image. Nevertheless, in [START_REF] Nayar | Subspace methods for robot vision[END_REF], [START_REF] Deguchi | Visual servoing using eigenspace method and dynamic calculation of interaction matrices[END_REF], [START_REF] Deguchi | A direct interpretation of dynamic images with camera and object motions for vision guided robot control[END_REF], the interaction matrix related to the eigenspace is not computed analytically but is estimated on-line from the estimation of the plane coordinates tangent to the cost function surface leading to unsatisfactory behavior. In [START_REF] Deguchi | A direct interpretation of dynamic images with camera and object motions for vision guided robot control[END_REF] the author points out the importance of using proper interaction matrices for visual servoing.

In this article, we clearly draw inspiration from Deguchi's previous work [START_REF] Deguchi | A direct interpretation of dynamic images with camera and object motions for vision guided robot control[END_REF]. Based on our own previous work on DVS [START_REF] Collewet | Photometric visual servoing[END_REF], we propose two new control laws based on a principal analysis decomposition of the main component of the image. In both case we first project the image on the new orthogonal basis and obtain a new and compact set of coordinates (coefficients). The former approach is a photometric visual servoing technique that considers, as input, an image reconstructed from a small number of coefficients. Indeed, an approximation of the image can be obtained as a linear combination of a small subset of coefficients and the eigenspace. Within this context, we remain in the photometric visual servoing scheme and provide an analytical formulation of the interaction matrix. Dealing with the latter approach, we consider the coordinates in the new basis as the visual features. Again, within this context we propose an explicit and analytical formulation of the interaction matrix. We show on various experiments including real 6 DoF positioning tasks, that these approaches feature a better behavior that the classical photometric visual servoing scheme allowing larger displacements and a satisfactory decrease of the error norm thanks to a well modelled interaction matrix.

In the reminder of this paper, Section II gives an overview of the DVS scheme. Section III recalls the principal component analysis concept and the way to build the eigenspace and obtain the new image coordinates. Section IV gives the details of the control laws including the derivation of the related interaction matrix. Finally, Section V illustrates the effectiveness of the approach: first with simulated results but also with experiments carried out on a 6 DoF Viper 850 robot.

II. DIRECT VISUAL SERVOING A. Positioning Task by visual servoing

The aim of a positioning task is to reach a desired pose of the camera r * , starting from an arbitrary initial pose. To achieve that goal, one needs to define a cost function that reflects, in the image space, this error. Most of the time this cost function is an error measure which needs to be minimized. Considering the actual pose of the camera r the problem can therefore be written as an optimization process:

r = arg min r e(r). (1)
For example, considering a set of geometrical features s, the task will typically have to minimize the error e(r) that is the difference between the current s(r) and the desired configuration s * , classically e(r) = s(r) -s * . This visual servoing task is achieved by iteratively applying a velocity to the camera. This requires the knowledge of the interaction matrix L s of s(r) that links the variation of ṡ to the camera velocity and which is defined as [START_REF] Espiau | A new approach to visual servoing in robotics[END_REF], [START_REF] Chaumette | Visual servo control, Part I: Basic approaches[END_REF]:

ṡ(r) = L s v ( 2
)
where v is the camera velocity. This equation leads to the expression of the velocity that needs to be applied to the robot. The control law is classically given by [START_REF] Chaumette | Visual servo control, Part I: Basic approaches[END_REF]:

v = -λL + s e(r) (3)
where λ is a positive scalar and L + s is the pseudo inverse of L s .

This approach requires to extract and track visual information (usually geometric features) from the image in order to design the control law. This difficult tracking process is one of the bottlenecks in the development of visual servoing techniques.

B. Photometric visual servoing

Recent works have tried to circumvent these problems by using directly the information provided by the entire image [START_REF] Collewet | Visual servoing set free from image processing[END_REF]. Features are no longer extracted from the image. In [START_REF] Collewet | Photometric visual servoing[END_REF], a control law was proposed that minimizes the error between the current image and the desired one. In that case the vector of visual feature in nothing but the image itself and the error to be regulated is the sum of squared differences (the SSD).

In that case, the feature s becomes the image itself (s(r) = I(r)). This means that the optimization process becomes [START_REF] Collewet | Photometric visual servoing[END_REF]:

r = arg min r (I(r) -I * ) (4)
where I(r) and I * are respectively the image seen at the position r and the template image (both of N pixels). The control law is inspired by the Levenberg-Marquardt (LM) optimization approach. It is given by:

v = -λL + I (I(r) -I * ) ( 5
)
where λ is a positive scalars and L I is the interaction matrix related to the luminance [START_REF] Collewet | Visual servoing set free from image processing[END_REF]. If we introduce the interaction matrices L x and L y related to the coordinates x and y of x, we obtain

L I = -(∇I x L x + ∇I y L y ) (6)
where ∇I x and ∇I y are the components along x and y of the image gradient ∇I. Note that it is actually the only image processing step necessary to implement the photometric visual servoing method. This approach features many advantages: it does not require any matching or tracking process. Furthermore, since the image measurements are nothing but the pixel intensity, there are no error in the feature extraction process leading to a very precise realization of the task.

III. PRINCIPAL COMPONENT ANALYSIS (PCA) It is possible to reduce the image high dimensions, by considering Principal Component Analysis (PCA) techniques. PCA is a standard technique which enables to linearly-project high-dimensional samples onto a low-dimensional feature space, that is called eigenspace. Such method has been shown to be very well-suited for object recognition [START_REF] Turk | Face recognition using eigenfaces[END_REF], appearancebased tracking [START_REF] Black | Eigentracking: Robust matching and tracking of articulated objects using a view-based representation[END_REF], keypoint tracking [START_REF] Ke | PCA-SIFT a more distinctive representation for local image descriptor[END_REF]. As stated, in visual servoing, this approach has been initially proposed by [START_REF] Nayar | Subspace methods for robot vision[END_REF] and [START_REF] Deguchi | Visual servoing using eigenspace method and dynamic calculation of interaction matrices[END_REF] [START_REF] Deguchi | A direct interpretation of dynamic images with camera and object motions for vision guided robot control[END_REF].

A. Building the eigenspace

Considering PCA for visual servoing requires two steps. The former is an off-line step which consist in building the eigenspace (a set of eigen images) from a large number of arbitrary images of the considered scene. This can be considered as a learning step. The building of the eigenspace consists in following steps:

• acquire M images (with N pixels) and build 1D column vectors that contain all the pixels of the image. We thus have M vector I 1 , . . . , I M . We assume that M N . • we then build a N × M matrix A whose columns are the normalized image vector

A •i : ∀i = 1 . . . M, A •i = I i -Ī, with Ī = 1 M M i=1 I i (7)
• Compute the eigenvalues σ i and their corresponding eigenvectors U i of the covariance matrix C given by

C = 1 M M i=1 A •i A •i = AA (8)
Note that C is then a N × N matrix. Considering a Singular Value Decomposition (SVD) of A allows to efficiently compute the eigenvectors U •i . Indeed, when decomposing A such that A = UΣV , U •i is the eigenvector associated to the eigen value σ i in Σ. For simplicity issue, we now denote U i the ith eigenvector. This a vector of size N . • Finally, we can keep only K eigenvectors corresponding to K largest eigenvalues. These vectors create a new orthogonal basis named eigenspace of dimensions K (K < M ).

B. Image decomposition and reconstruction

Once an eigenspace is built, we have a new orthogonal basis U 1 , U 2 , . . . , U M on which we can project each new image and obtain a new and more compact representation of this image. The k-th coordinates of an incoming image (with k < K) in this new basis is given by:

w k = U k (I -Ī) (9)
When projecting an incoming image on this basis, the greatest variance comes to lie on the first coordinate w 1 , the second greatest variance lies on w 2 , etc. When k increases the component w k becomes less and less significant. It can be noted that the proportion of the variance that each eigenvector represents can be calculated by dividing the corresponding eigenvalue by the sum of all eigenvalues. Following this principle, an image can be represented as linear combination of the eigenimage U i and described only by a vector w = (w 1 , . . . , w K ) with (K M ) which is a very compact representation. If one wants to reconstruct the image from the vector w, it can be computed as

I R = Ī + [Uw] K with [Uw] K = K k=1 w k U k (10)
with K ≤ M .

[Uw] K is the projection of the image I onto the subspace defined by the first K basis vectors. I R is then approximation of I.

IV. PCA-BASED VISUAL SERVOING

We now present how the PCA can be used within a visual servoing control law.

A. Photometric visual servoing on the reconstructed images

As stated in the previous section the current image I(r) can be decomposed using equation ( 9) and then reconstruct using equation [START_REF] Dame | Mutual information-based visual servoing[END_REF]. The image [Uw(r)] is an approximation of I(r) that correspond to the least-square estimate of w [START_REF] Black | Eigentracking: Robust matching and tracking of articulated objects using a view-based representation[END_REF]. The coefficients w are those that minimize [Uw(r)] -I(r) 2 . Note that this reconstruction can be achieved using a limited number of eigenvectors K. In that case, and following the methodology presented in [START_REF] Collewet | Visual servoing set free from image processing[END_REF], [START_REF] Collewet | Photometric visual servoing[END_REF], the error to be minimized is given by:

e(r) = I R (r) -I * R = [Uw(r)] K -[Uw * ] K (11)
Within equation [START_REF] Deguchi | A direct interpretation of dynamic images with camera and object motions for vision guided robot control[END_REF], U is constant and only the coefficient w(r) depends of the camera pose. [Uw(r)] being an image (we discard subscript K without loss of generality), it can be use as is within a photometric visual servoing process, as defined in section II-B, that minimizes e(r). The interaction matrix related to [Uw(r)] is given by L [Uw(r)] defined by:

L [Uw(r)] = -(∇[Uw(r)] x L x + ∇[Uw(r)] y L y ) (12)
with

L x = -1/Z 0 x/Z xy -(1 + x 2 ) y L y = 0 -1/Z y/Z 1 + y 2 -xy -x (13)
The control law is then given by (see Section II-B):

v = -λ L [Uw(r)] + ([Uw(r)] K -[Uw * ] K ). (14)
Such an approach based on reconstructed images allows to remain in a known framework introduced in [START_REF] Collewet | Visual servoing set free from image processing[END_REF]. Nevertheless, the considered images are built from the K eigenimages that contains the greatest variance. It still contains rich information and mainly discards the details.

B. Visual servoing from the compact eigenspace representation

An alternative to the image reconstruction process presented in the previous paragraph is to directly use the vector w as the visual feature. This what first was proposed in [START_REF] Deguchi | Visual servoing using eigenspace method and dynamic calculation of interaction matrices[END_REF] [START_REF] Deguchi | A direct interpretation of dynamic images with camera and object motions for vision guided robot control[END_REF]. In these works, the interaction matrix was estimated on-line leading to important approximations and sub-optimal robot motion. In this paper we proposed an analytical formulation of the interaction matrix.

The error to be minimized is then given by:

e(r) = w(r) -w * (15)
with, for k = 1..K:

w k (r) = U k (I(r) -Ī) (16)
Having defined the cost function to be minimized, one has to compute the interaction matrix L w k = ∂w k ∂r that links the variation of w k to the camera motion. Within the former equation U k and Ī are constant thus

∂w k ∂r = ∂U k (I(r) -Ī) ∂r = ∂U k I(r) ∂r (17)
since ∂U k Ī/∂r = 0 ( Ī being constant), pursuing the derivation of equation ( 17), we have:

∂w k ∂r = ∂U k ∂r I(r) + U k ∂I(r) ∂r (18)
Since U k is a constant not depending of r (it has been learnt off-line), this leads to:

L w k = U k L I . ( 19
)
where L I = ∂I(r)/∂r is given by:

L I = -(∇I x L x + ∇I y L y ) (20)
The complete control law is then given by:

v = -λL w + w(r) -w * . ( 21
)
with L w = U K L I where U K are the K first columns of U.

From a practical point of view we considered a Levenberg-Marquardt-like control law given by:

v = -λ H + µ diag(H) -1 L w w(r) -w * . ( 22
)
with H = L w L w is an approximation of the Hessian. More precisely, each component of the gradient is scaled according to the diagonal of the Hessian, which leads to larger displacements along the direction where the gradient is low. Such a control law has proven its effectiveness in a context of DVS [START_REF] Collewet | Photometric visual servoing[END_REF], [START_REF] Dame | Mutual information-based visual servoing[END_REF]. Note that beside gains λ and µ in equation ( 22) and, obviously K, no parameters are involved in these experiment. In all the experiment described below, we set λ = 1 and µ = 0.01. µ decreases by a factor 0.99 at each iteration. Therefore, the control law tends to the classical visual servoing control law that is similar to a Gauss-Newton minimization process (see equation ( 3)). With respect to [START_REF] Deguchi | A direct interpretation of dynamic images with camera and object motions for vision guided robot control[END_REF] we have there an analytic formulation of the interaction matrix (or image Jacobian) and it is not necessary to estimate it on-line from the image sequence from the estimation of the plane coordinates tangent to the cost function surface leading to a more precise calculation of L w .

V. EXPERIMENTAL RESULTS

Experiments have been carried out both in simulation and on a 6-DOF anthropomorphic robotic arm (a Viper 850 from Adept Company) equipped with a camera mounted on the end-effector. The camera calibration as well as the end-eye calibration have been done in an off-line step. The image processing and the control law computation are performed on a PC equipped with a Dual-core 2.4 Ghz Intel Pentium. The code has been written in C++ using the ViSP library [START_REF] Marchand | ViSP for visual servoing: a generic software platform with a wide class of robot control skills[END_REF]. The time required for an iteration of the VS closed loop is linear wrt. the number of selected coefficients K. In our experiment an iteration corresponds to 31ms for K = 20 and 73ms for K = 50 (including image acquisition). We mainly considered two scenes. One is a planar scene (hollywood) that also corresponds to the simulation experiments. The second is a 3D castle (castle3d) that is 0.2m high (see Figure 1). When dealing with simulation we considered the same calibration parameters as the real system leading to realistic simulations. Simulations have been done thanks to the ViSP simulator [START_REF] Marchand | ViSP for visual servoing: a generic software platform with a wide class of robot control skills[END_REF]. A video of the results can been at https://youtu.be/HncORj8Hpjk.

A. Construction of the eigenspace

The construction of the eigenspace is done thanks to the acquisition of M image I i , i = 1...M of the scene. It is done by just moving the camera handled around the scene. As far as simulation is concerned we generate random camera positions and get the corresponding images. Typically between 2500 and 8000 images are acquired (5000 for the simulations). A SVD decomposition of A (see equation ( 7)) is achieved to compute the eigenvectors U k . Since A is a large matrix (number of pixel × number M of images), it takes time (typically 522s when A is of size 66000 × 5000). This is done off-line and U is saved and used for all the experiments. Figure 2 shows the firsts and the last eigenimages for the two scenes hollywood and castle3d. As expected, the last eigenimage accounts mainly for noise.

The determination of the number K of coefficients is an important problem. As stated in section III-B, each coefficient explains a given percentage of the variance of the data which is directly related to the normalize value of the corresponding eigenvalue (see Figure 2). If we compute the cumulative value of the normalized eigenvalue, one can choose K such that it explains, eg, 75% of the variance which correspond to K = 20 on Figure 2). If the M images used to construct the eigenspace are well chosen, less coefficients are necessary to explain the same variance. In practice, we never used more than 50 coefficients.

B. Simulation results

Simulation results have first been carried out in order to validate the proposed control laws while allowing a fair comparison of different direct visual servoing approaches. We first compare the photometric visual servoing control law [START_REF] Collewet | Visual servoing set free from image processing[END_REF], [START_REF] Collewet | Photometric visual servoing[END_REF] with the proposed one. The error between the initial and desired pose is ∆r = (-0.05m, -0.1m, -0.03m, -5 o , 5 o , -15. metric visual servoing results, where the cost function to be minimized is given by e(r) = I(r)-I * , are shown on Figure 3 (we called it exp-photo). The results for the approach presented in section IV-A, where the cost function to be minimized is given by e(r) = I R (r) - 4 (this is exp-photo-rec). Finally, the proposed subspace approach, where the cost function to be minimized is given by e(r) = w(r) -w * is presented in Figure 6 for K = 6 (this is exp-w6) and in Figure 7 for K = 20 (this is exp-w20). For all the experiment we show the norm of the cost function, the camera velocity (in m/s and rad/s) and the positioning error (in m and rad).

I * R = [Uw(r)] K -[Uw * ] K is presented in Figure
When considering pure photometric visual servoing (Experiment exp-photo, Figure 3), the photometric cost function is highly non-linear which explain the perturbation in the velocity plots. Nevertheless, the visual servoing converges thanks to the redundancy of the considered information. Such experiment are in line with previous experiments (eg [START_REF] Collewet | Visual servoing set free from image processing[END_REF]). Dealing with photometric visual servoing from the reconstructed images (see Section IV-A and Experiment exp-photo-rec, Figure 4) the camera behavior is far better and the camera trajectory is closer from a geodesic. The camera velocities ( conclusion was reached by [START_REF] Crombez | Photometric Gaussian mixtures based visual servoing[END_REF] using photometric Gaussian Mixtures as visual features). Comparison between camera trajectories for exp-photo and exp-photo-rec is given in Figure 5.a.

We then consider directly the coefficients of the compact representation. Let us recall that this is a size K vector. Let us point out that K = 6, is the minimum number of coefficients necessary to have a full rank interaction matrix and thus to control the 6 DoF of the camera. We first consider two experiments the former with K = 6 and the latter with K = 20. In both cases, the visual servoing control law converges. Nevertheless, our experiments show that considering at least 15 coefficients leads to better behavior and a better monotonous decrease of the cost function norm (see 7.a vs 6.a). The convergence is also much faster than with the two previous photometric-based methods (140 vs 600 iterations). The small overshoot that can be observed in the camera velocities between iteration 35 and 80 (Figure 7) is due to decay of the parameters µ in equation ( 22). When µ is high (typically 0.01), the control law produces fast motion along the direction where the gradient is low. It then reach very quickly a position where the cost function is relatively low and where the error in the image space (I -I * ) is very small although the error in the cartesian space is still important. Then when µ decreases, the control law tends to the classic one which is more precise when the cost function decreases (thanks to a better estimation of the interaction matrix). This explains We also show the results obtained for two specific motions: a pure 0.10 m translation along the optical axis (see Figure 8) and a pure 45 o rotation around the optical axis (see Figure 9). In both cases, the control law generates a motion along the desired axis. This tends to show that, at least in these cases, there are almost no coupling between translational and rotational camera motions leading to very satisfactory camera trajectories [START_REF] Chaumette | Potential problems of stability and convergence in image-based and position-based visual servoing[END_REF]. Dealing with the rotation which is quite large, a pure photometric approach fails and a classical IBVS method based on point coordinates would features a significant motion along the z axis (camera retreat). This motion exists here (see Figure 9.d) but is almost insignificant.

C. Experimental results on a 6 DoF robot

We report here three experiments carried out on the Viper 850 robot. The first two ones consider the hollywood scene and the later the castle3d scene. a) 6 DoF positioning task: This experiment reports a 6 DoF positioning task with respect to a planar scene with K = 50. The displacement to be achieved is ∆r = (0.04m, 0.27m, 0.04m, 22.3 o , 8 o , 26.3 o ). Let us point out that the transformation between the initial and desired poses (and particularly the rotation around the x and z axes) is very large and makes this experiment very challenging. This is also illustrated by the initial and desired images depicted in Figure 10(a-b). The norm of the cost function w(r) -w * decreases monotonously (Figure 10.e). The decrease in errors (Figure 10.d) is also highly satisfactory considering the fact that only the interaction matrix at the desired position and an approximated depth were employed. One can see on Figure 10.h that the control law tends to minimize rapidly the errors accounting for the z translation and z rotation (at iteration 200, the positioning error for these axis are almost null). The final error is ∆r = (0.0012m, 0.0023m, 0.00001m, -0.32 o , 0.19 o , 0.06 o ) which show the accuracy of the proposed approach. Although we consider here K = 50, the camera converges for smaller value of K. Figure 11 show the camera trajectories for K = 20 (blue) and K = 50 (green). The camera trajectories are close to a geodesic. b) Dealing with partial occlusions: In the next experiment we add large occlusions by two (non planar) objects (see Figure 12). We still consider the eigenspace learnt without occlusions (the same as for the previous experiment). We consider K = 50. We still consider large displacement in ∆r = (-0.18m, -0.07m, 0.04m, -18 o , 14.3 o , 28.6 o ). The final error is ∆r = (-0.0024m, -0.0017m, 0.0003m, 0.25 o , 0.25 o , 0.05 o ) which shows that the occlusions and the nonplanar scene does not affect the precision of the positioning task. Nevertheless, it can be noted that positioning errors decrease less monotonously due to larger modelling errors in the calculation of the interaction matrix. This experiment demonstrates the robustness of the control law to large displacements, partial occlusions, and modelling errors.

c) Dealing with a 3D scene: In order to demonstrate further the robustness of the approach, we propose to consider in this experiment a large 3D object (Figure 13). For this last experiment, we consider ∆r = (0.15m, 0.21m, -0.13m, 16.6 o , 13.7 o , 30.9 o ). We only consider K = 20 coefficients in vectors w and w * . Although, in this case, we trained a new eigenspace, it is important to note that depth is still assumed to be constant (Z = 0.6m), whereas the object high is almost 0.2m, leading to modelling errors in the interaction matrix. These modelling errors, that could be overcome with a RGB-D camera, affect the behavior of the control law as can be seen on the error and velocities plots. Nevertheless, the system converges and the final error is ∆r = (-0.0007m, -0.0008m, 0.0007m, 0.01 o , 0.01 o , 0.01 o ) (we increase the number of iteration leading to a better precision that for the two previous experiments).

VI. CONCLUSION

In this paper we demonstrated that direct visual servoing techniques can take advantage from the projection of the

Figure 1 .

 1 Figure 1. Experimental setup: camera mounted on a Viper 850 from ADEPT and the two considered scenes

Figure 2 .

 2 Figure 2. The first four eigenimages U 1 ...U 4 along with U 50 and the last one U 8358 for the hollywood scene. As expected, the last eigenimage accounts mainly for noise. Plot on left show the cumulative eigenvalues along with the eigenvalues themselves. One can see that with K = 20, we grasp almost 75% of the variance.

Figure 3 .

 3 Figure 3. Experiment exp-photo: Photometric visual servoing [7]. (a) Initial image I (b) desired image I * (c) error norm I(r) -I * (d) camera velocity (in m/s and rad/s) (e) positioning error (in m and rad).

Figure 4 .Figure 5 .

 45 Figure 4. Experiment exp-photo-rec: Photometric visual servoing from the reconstructed image [Uw(r)] K with K = 20. (a) initial reconstructed image [Uw(r)] K (b) desired reconstructed image [Uw * ] K (c) error norm [Uw(r)] K -[Uw * ] K (d) camera velocity (in m/s and rad/s) (e) positioning error (in m and rad).

Figure 6 .Figure 7 .

 67 Figure 6. Experiment exp-w6: Positioning task using compact eigenspace representation with K = 6 (a) error norm w(r) -w * (b) errors w i (r)w * i , i = 1..6 (c) camera velocity (in m/s and rad/s) (d) positioning error (in m and rad)

Figure 8 .Figure 9 .

 89 Figure 8. Positioning task using compact eigenspace representation with K = 20 for a pure translation motion ∆r = (0, 0, 0.1m, 0, 0, 0, 0) (a) error norm w(r) -w * (b) errors w i (r) -w * i , i = 1..20 (c) camera velocity (in m/s and rad/s) (d) positioning error (in m and rad)

Figure 10 .

 10 Figure 10. Experiment with real planar scene (a) initial image acquired by the camera I(r), (b) image I * acquired from the desired position, (c,d) error [Uw(r)] 50 -[Uw * ] 50 between reconstructed image for initial and desired position (a,b,c,d) are used for visualization but are not used in the algorithm. Only the error w(r) -w * plotted in (f) is considered, (e) w(r) -w * (g) camera velocity (in m/s and rad/s) (h) positioning error (in m and rad).

  a

Figure 11 .

 11 Figure 11. 3D camera trajectories for experiment (blue: K = 20 ; green: K = 50).

Figure 12 .Figure 13 .

 1213 Figure 12. Experiment real scene with partial occlusions (a) initial image acquired by the camera I(r), (b) image I * acquired from the desired position, (c-d) error [Uw(r)] 50 -[Uw * ] 50 between reconstructed image for initial and desired position used for visualization but not used in the algorithm. (e) errors w(r) -w * (f) positioning error (in m and rad).

The following notation has been used: ∆r = (t, θu), where t describes the translation part of the homogeneous matrix related to the transformation from the current to the desired frame, while its rotation part is expressed under the form θu, where u represents the unit rotation-axis vector and θ the rotation angle around this axis. This representation is also considered in the plot reporting the positioning error.


Sana Baklouti
email: sana.baklouti@insa-rennes.fr

Eric Courteille
email: eric.courteille@insa-rennes.fr

Philippe Lemoine
email: philippe.lemoine@ls2n.fr

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

Vibration Reduction of Cable-Driven Parallel Robots through Elasto-dynamic Model-Based Control

Keywords: Cable-Driven Parallel Robot, Model-based control, Elasto-dynamic pre-compensation

This paper deals with the elasto-dynamic model-based control of Cable-Driven Parallel Robots (CDPRs), which manifests in the coupling of a PID feedback controller with a model-based feed-forward control scheme. The feed-forward controller is derived from an inverse elasto-dynamic model of CDPR, which compensates the end-effector dynamics and specifically its vibrations due to cable elasticity. The integration of cable tension calculation into this control strategy guarantees positive cable tensions along the trajectory. Simulations and experimentations while using a suspended and non-redundant CDPR show that tracking errors and vibrations can be reduced by the proposed strategy compared to conventional rigid-body model-based control.

Introduction

Cable-Driven Parallel Robots (CDPRs) are a special kind of parallel robots, where rigid links are replaced by cables. A CDPR is decomposed into an endeffector (EE) connected to a fixed base frame through cables. CDPRs may be used in some application fields where industrial robots cannot be used due to limitation of their workspace, payload and the required cycle time. These specifications have attracted the interest of many researchers [START_REF] Dallej | Visionbased modeling and control of large-dimension cable-driven parallel robots[END_REF][START_REF] Nguyen | On the analysis of largedimension reconfigurable suspended cable-driven parallel robots[END_REF][START_REF] Izard | On the improvements of a cable-driven parallel robot for achieving additive manufacturing for construction[END_REF].

Thanks to their low inertia, CDPRs can reach high velocities and accelerations in large workspaces [START_REF] Lamaury | Dual-space adaptive control of redundantly actuated cable-driven parallel robots[END_REF], which may lead to vibrations due to cable elasticity. Two main issues related to CDPRs are discussed: pose stabilization and trajectory 10 tracking. While the pose stabilization aims to stabilize the robot in a reference pose, the trajectory tracking aims to have the EE following a reference trajectory.

The improvement of the robot performance can be done through the modification of the robot architecture either by optimizing the design [START_REF] Agrawal | Optimization and design of a cable driven upper arm exoskeleton[END_REF][START_REF] Yao | Dimensional optimization design of the four-cable-driven parallel manipulator in fast[END_REF][START_REF] Gagliardini | Discrete reconfiguration planning for cable-driven parallel robots[END_REF][START_REF] Gagliardini | A reconfiguration strategy for reconfigurable cable-driven parallel robots[END_REF] or by adding other components such as in [START_REF] Viet | A cable-passive damper system for sway and skew motion control of a crane spreader[END_REF][START_REF] Weber | Active vibration canceling of a cabledriven parallel robot using reaction wheels[END_REF]. Improving accuracy is still possible once the robot is in operation through a suitable control scheme.

Several controllers have been proposed in the literature to improve CDPR accuracy locally or on trajectory tracking [START_REF] Jamshidifar | Adaptive vibration control of a flexible cable driven parallel robot[END_REF][START_REF] Zi | Dynamic modeling and active control of a cable-suspended parallel robot[END_REF][START_REF] Pott | Ipanema: a family of cable-driven parallel robots for industrial applications[END_REF]. In [START_REF] Cuevas | Assumed-mode-based dynamic model for cable robots with non-straight cables[END_REF], the control of CDPR in the operational space is presented while considering non elastic but sagging cables through the Assumed Mode Method. In [START_REF] Merlet | Simulation of discrete-time controlled cable-driven parallel robots on a trajectory[END_REF], a discrete-time control strategy is proposed to estimate the position accuracy of the EE by taking into account the actuator model, the kinematic and static behavior of the CDPR.

Many papers deal with CDPR control while considering cable elongations and their effect on the dynamic behavior. A robust H ∞ control scheme for CDPR is described in [START_REF] Laroche | A preliminary study for h ∞ control of parallel cable-driven manipulators[END_REF] while considering the cable elongations into the dynamic model of the EE and cable tension limits. Assuming flexibility in the longitudinal direction of cable, a control strategy is proposed for CDPRs in [START_REF] Khosravi | Dynamic modeling and control of parallel robots with elastic cables: singular perturbation approach[END_REF], [START_REF] Khosravi | Stability analysis and robust pid control of cable driven robots considering elasticity in cables[END_REF], [START_REF] Khosravi | Dynamic analysis and control of cable driven robots with elastic cables[END_REF]. It consists in adding an elongation compensation term in the control law of a CDPR with rigid cables, using singular perturbation theory to get rid of undesirable vibrations. Here, cables are modeled by linear axial springs, but with constant stiffness. This control method is improved in [START_REF] Khosravi | Dynamic analysis and control of fullyconstrained cable robots with elastic cables: variable stiffness formulation[END_REF] by integrating variable cable stiffness. In this context, a robust adaptive controller is presented in [START_REF] Babaghasabha | Adaptive robust control of fully-constrained cable driven parallel robots[END_REF] to attenuate vibrations in presence of kinematic and dynamic uncertainties. This control method requires the measurement of cable lengths and the knowledge of the real-time end-effector EE pose through exteroceptive measurements. However, the external-measurement-based control methods add complexity to the cabledriven manipulators and restrict the application conditions due to the need of additional devices [START_REF] Cui | Closed-loop control for a cable-driven parallel manipulator with joint angle feedback[END_REF].

The importance of the feed-forward effect on non-linear systems control is highlighted in [START_REF] Slotine | Applied nonlinear control[END_REF]. It leads to stable systems with enhanced trajectory tracking performance. Feed-forward model-based controllers are used to fulfill accuracy improvement by using a CDPR reference model [START_REF] Zhang | Dynamic control with tension compensation of a 3-dof cable-driven parallel manipulator[END_REF]. This latter predicts the mechanical behavior of the robot; and then generates an adequate reference signal to be followed by the CDPR. This control type leads to the compensation of the undesirable effects without exteroceptive measurements. A model-based control scheme for CDPR used as a high rack storage is presented in [START_REF] Bruckmann | Design and realization of a high rack storage and retrieval machine based on wire robot technology[END_REF]. This research work takes into account cable elasticity. This strategy, integrating the mechanical behavior of cables in the reference signal, enhances the CDPR performance. However, the compensation for the EE pose errors due to the cable interaction with the whole system and elasto-dynamic behavior was not considered.

Consequently, we proposed an elasto-dynamic model-based control in [START_REF] Baklouti | Elasto-dynamic model-based control of non-redundant cable-driven parallel robots[END_REF] for only non-redundant CDPRs. This control method combines a model-based feed-forward control scheme for CDPRs with a PID feedback controller. The feed-forward controller is derived from an inverse elasto-dynamic model of CDPR, which compensates the EE deflection due to cable elasticity and vibrations due to dynamic interactions. The relevance of the proposed control method was checked only numerically. Therefore, in addition of the work reported in [START_REF] Baklouti | Elasto-dynamic model-based control of non-redundant cable-driven parallel robots[END_REF], this paper gives more details on the improved control strategy and contains experimental validations of the obtained theoretical results. Moreover, this paper deals with the improvement of the elasto-dynamic model-based control strategy to make it valid for any CDPR configuration through the integration of cable tension distribution into this control strategy, guaranteeing positive cable tensions along the trajectory. This paper is organized as follows: Section 2 presents the dynamic equations of motion of CDPR required to establish the control laws. Section 3 introduces the design of the elasto-dynamic model-based control and the definition of the corresponding control laws. Then, experimental validations performed on the CREATOR prototype located at LS2N, Nantes, France, a CDPR with three cables and three Degree-Of-Freedom (DOF), are discussed in detail in Section 4. Finally, the effectiveness of the proposed elasto-dynamic control strategy is studied numerically (Section 5) and experimentally (Section 6) to check its performance with regard to the CDPR positioning accuracy and trajectory tracking compared to the conventional ones.

Dynamic modeling of CDPR

The joint coordinate vector is denoted as q = [q 1 , ..., q n ] T ∈ R n , n being the number of cables. As shown in Fig. 1, the Cartesian coordinate vectors of anchor points A i and exit points B i , i = [1..n], are denoted as a i and b i . These vectors are expressed in the EE frame F p = {P, x p , y p , z p } and in the base frame F b = {O, x b , y b , z b }, respectively. The ith cable connecting exit point B i and anchor point A i , i = [1.

B 1 B n A 1 A n B i A i x p y p z p P End-effector Fixed-base x b z b y b O p a i b i l i d i
.n], exerts a force from A i to B i of magnitude τ i >0. The ith cable length l i is expressed as follows:

l i = p -b i + b R p a i 2 . ( 1
)
The relationship between the cable length vector l = [l 1 , ..., l n ] T ∈ R n and the motor angular displacement vector q is supposed to be linear:

l = χ q, (2)
where χ = diag[χ 1 , ..., χ n ] ∈ R n×n is a diagonal matrix containing the winding ratio of the winches. The wrench matrix W ∈ R m×n is a function of x and maps the EE velocities to the cable velocity vector. It is expressed as follows:

W = d 1 ... d n b R p a 1 × d 1 ... b R p a n × d n , d i = b i -p -b R p a i b i -p -b R p a i 2 . (3)
where d i ∈ R u is unit vector of ith cable pointing from A i to B i . b R p is the rotation matrix of the mobile frame F p with respect to the fixed frame F b .

The equations of motions in this paper are derived from Newton-Euler equations while considering the center of mass of the EE coincident with its geometric center, namely,

M(x) ṫ + C(x, t) t = Wτ + w g + w e . (4)
M ∈ R m×m is the mass matrix of the moving platform. C ∈ R m×m is the matrix of Coriolis and centrifugal forces. t ∈ R m is the twist vector of the EE. τ ∈ R n is the cable tension vector. w g ∈ R m is the wrench due to gravity acceleration. w e ∈ R m amounts to the other wrenches applied to the mobile platform, except w g .

Determination of cable tensions

The determination of τ is a function of the EE pose and the wrench that the latter should support. It is about solving the dynamic equilibrium equations for a given pose x of the EE, which can be described as follows:

Wτ + w ex = 0, w ex = w g + w e -M ṫ -C t, (5)
where w ex ∈ R m denotes the external wrench. The cables can only pull and not push the EE. If the number of cables n is equal to the DOF m, the inversion of Eq. ( 5) is possible, as long as the wrench matrix W is not singular.

However, this tension set can present negative values. That means that one or more cable(s) may have to push the EE, which is not possible. When m<n, Eq. ( 5) may have an infinite number of solutions. Therefore, the redundancy allows to select a solution amongst the infinite set cable tension vectors satisfying some criteria. The problem of force distribution presents one important design issue for redundant actuated CDPRs, which is the determination of feasible cable force distribution 1 [START_REF] Pott | An improved force distribution algorithm for over-constrained cable-driven parallel robots[END_REF][START_REF] Gosselin | On the determination of the force distribution in overconstrained cable-driven parallel mechanisms[END_REF][START_REF] Borgstrom | Rapid computation of optimally safe tension distributions for parallel cabledriven robots[END_REF][START_REF] Mikelsons | A real-time capable force calculation algorithm for redundant tendon-based parallel manipulators[END_REF].

In what remains, the pre-compensation of the end-effector deflection will be based on the estimated cable elongations due to the applied tension set. Depending on the cable stiffness, different reference models for control will be defined: (i) rigid model, (ii) elasto-static model [START_REF] Bruckmann | Design and realization of a high rack storage and retrieval machine based on wire robot technology[END_REF][START_REF] Arsenault | Workspace and stiffness analysis of a three-degree-of-freedom spatial cable-suspended parallel mechanism while considering cable mass[END_REF] and (iii) elasto-dynamic model [START_REF] Baklouti | Elasto-dynamic model-based control of non-redundant cable-driven parallel robots[END_REF].

Feed-forward model-based control

The feed-forward model-based control uses a priori knowledge of the CDPR dynamics to improve the EE trajectory tracking by generating a pre-compensated control input. The pre-compensation unit is seen as an off-line compensation of the actuated joint displacement in such a way that after execution the EE follows the desired motion more accurately.

The feed-forward model-based control scheme is shown in Fig. 2. It is composed of a feed-forward block in which the inverse kinematic model is determined based on a CDPR reference model (Red block in Fig. 2). This latter is a predictive model of the dynamic behavior of the mechanism. Its input is the torque set-point ζ rg = χτ rg ∈ R n , which is induced by the cable tensions τ rg ∈ R n . It depends on the EE desired motion (x rg , t rg and ṫrg ). It is determined based on the used tension distribution algorithm.

The output of the CDPR reference model is q ref ∈ R n , which represents the reference angular displacement of the motors. It depends on the chosen reference model. The joint angle vector q rg is calculated based on x rg with respect to Eq. ( 2). The joint velocity vector qrg and joint acceleration vector qrg are the time derivatives of q rg . The vector q represents the measured joint angle vector of motors.

In this paper, it is assumed that all components of an actuator (motor, reducer and winch) are coaxial and all actuators are identical. Therefore, the motor torques vector ζ m ∈ R n is expressed as follows:

ζ m = ζ corr + ζ f ( qrg ) + ζ rg , (6)
The vector ζ corr = Ih(t) corresponds to the torque of correction. The ith diagonal term of I = I 1 n×n is the total actuator inertia matrix seen at their 1 A cable force distribution will be said to be feasible in a particular configuration and for a specified set of wrenches if the cable tensions can counteract any external wrench of the specified set applied to the EE [START_REF] Ebert-Uphoff | On the connections between cable-driven robots, parallel manipulators and grasping[END_REF].

) 7
I motor is the inertia of the motor with brakes. I reducer is the inertia of the reducer. r is the gear-head (reducer) ratio. h(t) is defined by:

h(t) = qrg + K p e + K d ė + K i t + i ti e dt, (8)
where e = q ref -q and ė are the angular displacement error vector and its time derivative, respectively. K p ∈ R n×n is the proportional gain matrix. K d ∈ R n×n is the derivative gain matrix. K i ∈ R n×n is the integrator gain matrix. ζ f ∈ R n is the friction torque vector, which is calculated with respect to the static friction model [START_REF] Khalil | Modeling, identification and control of robots[END_REF] as a function of the desired motion of the EE:

ζ f ( qrg ) = ζ d tanh( qrg ) + ζ v qrg . (9)
ζ d ∈ R n×n and ζ v ∈ R n×n are two diagonal matrices containing the dry friction coefficients and the viscous friction coefficients of the actuators, respectively. By calculating off-line an appropriate trajectory pre-compensation, the effects of the oscillatory dynamics are compensated. However, this may lead to a non robust control with respect to model uncertainties and a delay between the reference and the real CDPR behavior may occur. To handle this issue, the off-line feed-forward part should respect the findings in [START_REF] Baklouti | Sensitivity analysis of the elastogeometrical model of cable-driven parallel robots[END_REF], dealing with a sensitivity analysis allowing to highlight the crucial parameters to be upstream identified. Accordingly, the control scheme of the CDPR is based on a robust model to uncertainties. This leads to the attenuation of errors coming from modeling and parameters uncertainties as shown for several industrial robots [START_REF] Olabi | Improving the accuracy of industrial robots by offline compensation of joints errors[END_REF][START_REF] Verdonck | Improving the dynamic accuracy of industrial robots by trajectory pre-compensation[END_REF].

Discrete-time control of CDPR

The correspondence between a CDPR and its control unit is described in Fig. 3. This organization of the equivalent control scheme is based on the findings in [START_REF] Merlet | Simulation of discrete-time controlled cable-driven parallel robots on a trajectory[END_REF], dealing with a simulation tool of a discrete-time controlled CDPR. It is composed of three units: (i) Actuation systems, (ii) Feedback controller and (iii) Feed-forward controller. Each unit of the equivalent control scheme is characterized by its own discretization time. Let ∆t a be the sampling time of the actuation system, which is a function of the feed-back controller frequency. ∆t b is the sampling time of the feed-forward unit. ∆t b is a function of the actuation system response. It is chosen smaller or equal to ∆t a . The feed-forward unit should be fast enough to avoid delays and divergence.

Actuation system Feedback control Feed-forward control x rg , x rg , x rg . .. rg , q ref q m { t = ( +j ) t a } { t = t a } Offline Online { t = α Δt b }

. Actuation system

This model presents the actuators of the CDPR. Its input is the motor torque vector ζ m and its output is the measured angular coordinates q and velocities q of the actuators. This measurement is obtained at each time t = β∆t a .

Feedback control loop

This control loop permits the application of the control law mentioned in Eq. ( 6). The feedback control loop has a sampling time ∆t a . At each time t = β∆t a , β ∈ N + ,the loop calculates a new value for ζ m , which is sent to the actuation model at time t = (β + j β )∆t a , at which a new measurement of q and q is made. j β is the delay between the generation of the control set-point by the feed-back controller and the reception of the encoder information. The value of

ζ m is kept constant in the time interval [(β + j β )∆t a (β + j β + 1)∆t a ].

Feed-forward control loop

This control loop aims to generate the modified trajectory. The input of this loop is the desired motion of the EE, which manifests in the desired pose x rg , twist t rg and acceleration ṫrg of the moving platform. The output is the reference angular displacement q ref , the torque set-point ζ rg , the friction torque ζ f and the angular acceleration qrg generated from the desired motion of the moving-platform. The sampling time of the feed-forward control loop is ∆t b .

At each time α∆t b , a torque set-point ζ rg is calculated, coming from a cable tension distribution algorithm, with respect to the desired motion of the EE. The pre-compensation of the trajectory is obtained by generating the joint angular displacement q ref with respect to the chosen reference model of the CDPR.

Pre-compensation

The off-line pre-compensation strategy is applied to eliminate the trajectory tracking errors due to the robot compliance. This strategy is based on the anticipation of the joint errors and correcting them by modifying the reference actuated joint position q ref as follows:

q j ref = q rg -δq j mod , δq j mod = q j mod -q rg , (10)
where j refers to the reference model control. It can be either "rg", "es"

or "ed" to mention the rigid, elasto-static and elasto-dynamic models, resp. δq j mod ∈ R n is a vector of angular displacements errors, leading to trajectory tracking errors. This vector is calculated by simulating a CDPR model depending on the considered error sources. It anticipates the correction of the actuated joint angles q j mod with respect to the displacement vector q rg , which is calculated through the rigid model of CDPR and corresponds to the unstrained cable length l rg . The subtraction of δq j mod from q rg generates a new reference input signal q j ref allowing the EE to track more accurately the desired trajectory.

It is noteworthy that δq j mod depends on the CDPR reference model used to calculate the vector q j ref .

To the best of our knowledge, two CDPR models have been used in the literature for the feed-forward model-based control of CDPRs with non-sagging cables: (i) rigid model and (ii) elasto-static model [START_REF] Bruckmann | Design and realization of a high rack storage and retrieval machine based on wire robot technology[END_REF][START_REF] Arsenault | Workspace and stiffness analysis of a three-degree-of-freedom spatial cable-suspended parallel mechanism while considering cable mass[END_REF]. As a consequence, one contribution of this paper deals with the use of the (iii) elasto-dynamic model of CDPRs for feed-forward control [START_REF] Baklouti | Elasto-dynamic model-based control of non-redundant cable-driven parallel robots[END_REF].

Rigid pre-compensation

The CDPR rigid model considers cables as rigid links. Here, the reference signal q rg ref anticipates neither the cable elongation nor the oscillatory motions of the EE. It is expressed as:

q rg ref = q rg . (11)
The PID feedback controller uses the motor encoders response q, which is related to l rg . It should be noted that the cable elongations and EE oscillatory motions are not considered here, and as a consequence, cannot be rejected.

Elasto-static pre-compensation

The CDPR elasto-static model integrates a feed-forward cable elongation compensation [START_REF] Bruckmann | Design and realization of a high rack storage and retrieval machine based on wire robot technology[END_REF]. At each sampling time α∆t b , each cable of the CDPR is isolated and its elongation is determined with respect to the discretized EE motion (x rg , t rg and ṫrg ). Knowing the force applied on each cable, we estimate its elongation based on a cable tension model, which can be either linear or non-linear [START_REF] Baklouti | Dynamic and oscillatory motions of cable-driven parallel robots based on a non-linear cable tension model[END_REF]. In this paper, we consider the linear cable tension model, where the cable tension is linearly proportional to cable stiffness. When q rg ref is used as a reference signal in the feedback control scheme, the EE displacement δx es is obtained from cable elongation vector δl es . To compensate for the cable elongation effects, δl es is converted into winch displacement δq es mod , which estimates the joint error with respect to the rigid angular position q rg . Thus, the elasto-static reference angular displacement q es ref becomes:

q es ref = q rg -δq es mod . (12)
The reference signal q es ref corresponds to a fake position of the EE for the cable elongation compensation. Here, under the effect of cable elongations, the reference EE pose is estimated to achieve the desired pose. Although the elasto-static reference model takes into account the cable elongations, the noncompensation for the EE pose errors due to the cable interaction with the whole system and elasto-dynamic behavior is not considered.

Elasto-dynamic pre-compensation

The CDPR elasto-dynamic model takes into account the oscillatory and dynamic behavior of the EE due to cable elongations. Here, the cables are no-longer isolated and are affected by the EE dynamic behavior. The real EE pose is expressed as: x ed = x rg + δx ed . The EE displacement δx ed leads to variations in both cable lengths and cable tensions. Indeed, the ith cable tension τ i ed obtained from the elasto-dynamic model differs from τ i rg calculated based on a tension distribution algorithm:

τ i ed = τ i rg + δτ i ed = ES δl i ed δl i ed + l i rg , ( 13
)
where E is cable modulus of elasticity and S is its cross sectional area. δl i ed is the ith cable elongation assessed by considering cable elasticity and oscillations. The EE pose x ed is calculated through the direct elasto-dynamic model thanks to MATLAB ® routine ode45 by coupling Eq. ( 4) and Eq. ( 13). The cable elongation vector δl ed is converted into δq ed mod , which corrects the angular position vector q rg . Therefore, the elasto-dynamic reference angular displacement q ed ref becomes:

q ed ref = q rg -δq ed mod . (14)
The proposed control strategy based on the elasto-dynamic model leads to a feed-forward controller for EE oscillatory motion compensation in addition to the conventional rigid body feedback while considering the motor encoder measurements.

Controller tuning

No matter the reference model used in the control scheme (see Fig. 2), the controller tuning methodology is the same. Applying the Laplace transform to the tracking error differential Eq. ( 8), and assuming zero initial conditions, we obtain for each drive:

e(s)(s 3 + K p s + K d s 2 + K i ) = 0. ( 15
)
Equation. ( 15) can be rewritten as follows:

e(s)(s + Υ n )(s 2 + 2υΥ n s + Υ n 2 ) = 0, (16)
with

K p = (2υ + 1)Υ 2 n , K d = (2υ + 1)Υ n and K i = υ 3 n .
This differential equation expresses the dynamic behavior of the closed-loop tracking error in response to a perturbation. This corresponds to a first-order system followed by a second-order system. Since the parameters υ and Υ n are strictly positive, the system is always stable.

To obtain the fastest answer without oscillation, υ = 1 is chosen. Thus, a third order system with a triple real pole is obtained. In this case, the PID gains become:

K p = 3Υ 2 n , K d = 3Υ n , K i = Υ 3 n . ( 17
)
Note that K p , K d and K i are all functions of Υ n . A simple and effective method is to select Υ n is twice the breaking pulse Υ c = k t k e I of the motors where is the resistance of the motor armature, k e is the counter-electromotive force coefficient and k t represents the electromagnetic torque coefficient of the motor.

Control of a 3-DOF spatial CDPR with a point-mass end-effector

This section aims to describe the CDPR prototype to be used to check the relevance of the elasto-dynamic feed-forward control with respect to the classical rigid and elasto-static ones. The variations in the EE position error along a trajectory are numerically and experimentally studied on the CREATOR prototype located at LS2N. The CDPR configuration used for the tests is composed of three cables and a point-mass EE. The EE is requested to track the linear path defined in Section 4.2.

Test bench

Figure 5 describes the equivalent architecture of the CREATOR prototype. The test bench is prepared under the assumption that cables are elastic and keep linear shape along a trajectory. The EE and pulleys are considered as points. To provide this configuration, some precautions are taken into account on the experimental test bench to provide the particular configuration of the CREATOR prototype:

• The CREATOR prototype is actuated by three Parker ™ motors 2 with gear-heads 3 connected to 3D printed winches. The motors characteristics The friction torques of the actuators are identified with respect to the static friction model [START_REF] Khalil | Modeling, identification and control of robots[END_REF] by incrementing the joint angular velocity. The do not have an accurate measurement of the motor torques, we suppose that the viscous friction torque is zero. Here, the dry friction value is

ζ d =0.14 N.m.
• The winch is 3D printed and it is mounted to the reducer axis. Its pitch is equal to 2 mm and its radius is equal to 50 mm.

• 3D printed pieces are used as pulleys of the CREATOR protype (See Fig. 6). They were designed in such a way that their exit amount to a point. The entrance of the pulley is conical and their exit is a hole of 3 mm of diameter. Grease is added inside the pulley to reduce friction between the cable and the pulley. • The robustness analysis in [START_REF] Baklouti | Sensitivity analysis of the elastogeometrical model of cable-driven parallel robots[END_REF] highlighted the high effect of the cables elasticity onto the dynamic and oscillatory motions of the EE. The CRE-ATOR cables are made up of eight threads of polyethylene fiber with a diameter of 0.5 mm. These cables were experimentally upstream identified.

The identification method is described in [START_REF] Baklouti | Dynamic and oscillatory motions of cable-driven parallel robots based on a non-linear cable tension model[END_REF]. The absolute uncertainties in the applied force and resulting elongation measurements from the test bench outputs are estimated to be ± 1 N and ± 0.03 mm, respectively [START_REF] Baklouti | Dynamic and oscillatory motions of cable-driven parallel robots based on a non-linear cable tension model[END_REF].

The resulting modulus of elasticity of the cable is equal to 70 ± 1.51 GPa.

The following trajectory tracking experiments will confirm a good identification of the modulus of elasticity. This is shown through the measurement of a natural frequency ≈ 3.5 Hz, which corresponds to the calculated natural frequency (3.67 Hz) while using the identified modulus of elasticity.

The reference [START_REF] Kozak | Static analysis of cable-driven manipulators with non-negligible cable mass[END_REF] was used for the calculation of the natural frequencies of the CDPR. It is based on the global stiffness matrix of the CDPR. The natural frequencies of the system are obtained by solving the eigenvalue problem of the linearized free vibration equation driven by the generalized stiffness matrix of the CDPR and the mass matrix of the end-effector. • The EE was designed such that it is close to a point-mass. It is a steel machined cylinder with a hole that crosses vertically and a conical hole on its top surface. This makes the cables to intersect at the same point.

The magnetic support of the laser tracker target is connected to the EE at its bottom part. The laser tracker target is then connected to this magnetic support, Fig. 7. The acquired information from the laser tracker presents measurements with respect to the target center. An offset with respect to the EE gravity center is considered to get the appropriate trajectory tracking measurements.

The total mass of the EE used for the control gathers its own mass plus the mass of the laser tracker target and its magnetic support. The total mass of the EE is equal to 0.780 kg.

• The command of the CREATOR prototype is implemented in a host PC through a software interface generated by ControlDesk ®4 . This latter enables the real-time implementation of the control schemes, created with Matlab-Simulink ® , in the dSpace ™ control unit. The control sampling times ∆t a and ∆t b are considered equal: ∆t a = ∆t b = 1 ms.

• An API Radian laser tracker was used to define the base frame F b . A circle plane passing through points B 1 , B 2 and B 3 is created and is projected on the floor to create the base frame, whose origin is the circle center. All the measurements are expressed in that frame. Moreover, the laser tracker is used to measure the pose of the EE along the tracked trajectory in the global base frame F b . The static measurement accuracy of the used laser tracker is ±10 µm. The measurement start signal is synchronized with the control signal through the dSpace controller.

Trajectory generation

Before running trajectory tracking, a pre-compensation of the EE static 305 deflection is performed while adopting the CDPR control based on its elastostatic and the elasto-dynamic models, respectively. Thus, the static deflection due to the EE weight is compensated by correcting the joint angular displacement leading to trajectory tracking error.

A linear path along the vertical axis with a 5-th degree polynomial motion profile is considered. It goes from the initial position p 1 to the final position p 2 , 1 m high along z-axis, during t f =3 s. A pause equal to t f at the final pose is considered. The effects of the EE residual oscillations are easily detected during the steady-state phase. Therefore, the EE trajectory is parametrized as follows: The fifth degree polynomial for the linear path is expressed as:

p(t) = p 1 + α(t) (p 2 -p 1 ); t ∈ [0 t f ], ( 18a
)
p(t) = p 2 ; t > t f . ( 18b
α(t) = β 5 t 5 + β 4 t 4 + β 3 t 3 + β 2 t 2 + β 1 t + β 0 , (19)
where

β 5 = 6 t 5 f , β 4 = -15 t 4 f , β 3 = 10 t 3 f
, β 2 = 0, β 1 = 0 and β 0 = 0. These values are chosen such that the velocity and acceleration of the EE are null at the beginning and the end of the trajectory as shown in Figs. 8b and8c.

Numerical simulations

The three control schemes under study were simulated through Matlab-Simulink ® to analyze their relevance. To assess the performances of different control laws, the CDPR elasto-dynamic model is used to estimate the EE real position as it is the closest one to the real CDPR with non sagging cables. This is to be confirmed in Section 6. Therefore, it is used to predict the real behavior of the CDPR. The input of this model is ζ m , which leads to the Cartesian coordinate vector p m ed of the EE position. The trajectory error is defined as: 9b shows the norm δp of the end-effector trajectory tracking error when the proposed model-based feed-forward control law is applied. The three CDPR models are successively used to generate the reference signal. Figure 9a illustrates the end-effector trajectory tracking error δz along the z-axis, which is the main one as the CDPR under study is assembled in a suspended configuration and its end-effector tracks a vertical trajectory. The red (green, blue, resp.) curve depicts the end-effector trajectory tracking error when the elasto-dynamic (elasto-static, rigid, resp.) model is used as a reference model. The position error norm δp while using the rigid model as a reference for the feed-forward at time t = 3 s is equal to 13.14 mm. This value is equal to 2.18 mm when the elasto-static model is used as a reference for the feed-forward, which represents a relative difference of 83 % with respect to the rigid control 5 . The position error norm δp while using the elasto-dynamic model as a reference for the feed-forward at time t = 3 s is equal to 0.58 mm, which represents a relative difference of 95 % with respect to the rigid control.

δp(t) = p m ed (t) -p rg (t). Figure
This shows a good improvement of the moving-platform positioning accuracy at the end of the trajectory through the elasto-dynamic compensation comparing to classical feed-forwards. The trajectory tracking error is also reduced significantly with the CDPR elasto-dynamic control6 compared to the classical feed-forward model-based control schemes. This is due to the fact that the elastodynamic control compensates not only the static errors due to cable elongations but also the oscillatory behavior of the moving-platform. Next, experimental results are performed to verify numerical ones.

Experimental results

Experimental setup

The experimental procedure used to verify the relevance of the elasto-dynamic model-based control is the following: • A support for the EE, shown in Fig. 10, is placed on its desired initial position. The base plate of this support is 3D printed and is connected to a steel block through steel square tubes. This base is a negative form of the EE. The support is useful to define the relationship between actuated angles and the corresponding unstrained cable lengths. Besides, it provides the same initial pose for the EE for the different experiments.

As shown in Fig. 10, the EE is first located on the support. Thus, it is in the non-suspended initial pose and it undergoes neither the external effects nor its stiffness. The absolute position of the EE is measured in that non-suspended configuration.

• Modify cable lengths such that their tension is a minimum while keeping cables without sag. Then, lock the motors.

• The support is removed and the EE becomes suspended under the effect of its own weight subjected to cable elasticity. Then, the suspended position of the EE is measured.

• A command is generated through ControlDesk ® in combination with Matlab-Simulink ® . It makes the EE move along a straight line.

• When the elasto-static or elasto-dynamic models are used as references for the control, a compensation of the static deflection is performed before starting the trajectory. This is done through smooth angular displacements of joints corresponding to the correction of the static deflection.

Static analysis

The ideal initial position of the EE corresponds to the initial rigid position p rg (t = 0 -s). This latter is measured when the EE is located on the support in the absence of the stiffness effect; p rg (t = 0 -s) = [0.299, -0.047, 0.623] m.

Three tests are performed to check the positioning repeatability. This is done by suspending the EE and measuring its corresponding static deflection. The suspended position p es (t = 0 -s) of the EE is measured when the support is removed and the EE is suspended under the effect of its own weight. Table 3 displays the static deflection dp with respect to the same nonsuspended rigid position p rg (t = 0 -s) for different tests. The EE presents a different static deflection for each test. This depends on the initial tensions on the cables when the EE is in its ideal initial position. The absence of tension sensors makes it hard to start the experiments with zero cable tensions without sag, leading to different static deflections. The position repeatability is not good enough, so for each test we obtain a different position for the EE. Indeed, the trajectory does not start with the same position for the different tests. Comparing the control laws by comparing the absolute trajectory errors will not project the true improvement of the accuracy as the different trajectories do not start with the same static error. For this fact, we propose to define a new reference p * to be used for the determination of the trajectory error. p * is defined as described in Sec. 4.2 starting from the compensated position, where p * (t = 0 s) = p es (t = 0 + s). This allows to compensate the static error for the different tests at the start time of the trajectory. This static error should be seen at the end of the trajectory.

Trajectory tracking

The trajectory error δp is defined as the difference between the measured trajectory p and the reference one p * : δp(t) = p(t) -p * (t). Figure 12d shows the norm of the EE trajectory tracking error δp when the proposed feedforward control law is applied while using successively the three CDPR models to generate the reference signal. Figure 12c illustrates the EE trajectory tracking error along the z-axis δz, which is the main one as the CDPR under study is assembled in a suspended configuration and the trajectory is vertical. The red (gree, blue, resp.) curve depicts the EE trajectory tracking error when the elasto-dynamic (elasto-static, rigid, resp.) model is used as a reference. As one can see, the set-point of current is not constant during the steady-state phase. This reflects the integrator effect of the PID controller as its input is never zero. In addition, the torques ζ rg and ζ f do not compensate totally the gravity and friction effects as they are obtained by identification, which is influenced by uncertainties.

While applying i m , the resulting trajectory error of the EE is plotted in Fig. 12. It is worth noting in these plots that the frequency of EE vibrations is equal to 3.5 Hz, which is close to the calculated first natural frequency f 1 =3.67 Hz.

A good identification of the modulus of elasticity, respecting the uncertainty analysis conclusions, is confirmed. From Figs. 11 and12, the delay j β can be determined. When the rigid model of the CDPR is used, a current discontinuity at time 2.53 s is depicted. This leads to a discontinuity in δp at time 2.68 s, which corresponds to j rg β =0.15 s. When the elaso-static model of the CDPR is used, a current discontinuity at time 2.73 s is depicted. This leads to a discontinuity in δp at time 2.87 s, which corresponds to j es β =0.14 s. When the elasto-dynamic model of the CDPR is used, a current discontinuity at time 2.57 s is depicted. This leads to a discontinuity in δp at time 2.70 s, which corresponds to j ed β =0.13 s. Note that the end-effector position errors along x-axis and z-axis are mainly due to the end-effector swaying, the end-effector being not a point-mass. Besides, the experimental validations were performed with a CDPR suspended configuration, which leads to a low horizontal stiffness of the end-effector. This may lead to end-effector swaying and large horizontal oscillations under the effect of external disturbances, especially when using a light end-effector.

As shown in Fig. 12, the trajectory tracking error norm δp while using the rigid model as a reference for the feed-forward at time t=3 s is equal to 12.05 mm. This value is equal to 3.45 mm when the elasto-static model is used as a reference for the feed-forward, which represents a relative difference of 71 % with respect to the rigid control. The position error norm δp while using the elasto-dynamic model as a reference for the feed-forward at time t=3 s is equal to 2.01 mm, which represents a relative improvement of 83 % with respect to the rigid control.

This confirms a better improvement of the EE positioning errors through the elasto-dynamic compensation comparing to classical feed-forwards. This is due to the correction of static errors through the compensation of cable elongations and the EE oscillatory effects. As shown in Fig. 12, at time t=3 s, the EE achieves its final position and then oscillates. Fig. 13 depicts that when the rigid control is used, the first Peak-to-Peak amplitude of residual vibrations is equal to 2.64 mm. This value is equal to 2.48 mm when the elasto-static model is used as a reference for the feed-forward, which represents a relative difference of 6 % with respect to the rigid control. When the elasto-dynamic model is used as a reference for the feed-forward, the first Peak-to-Peak amplitude of residual vibrations is equal to 1.32 mm, which represents a relative improvement of 50 % with respect to the rigid control.

Steady-state phase

The fifth Peak-to-Peak amplitude of residual vibrations is equal to 0.31 mm when the rigid model is used as a reference for the feed-forward. This value is equal to 0.82 mm when the elasto-static model is used as a reference for the feed-forward. The fifth Peak-to-Peak amplitude of residual vibrations is equal to 0.21 mm when the elasto-dynamic control is used, which represents a relative improvement of 32 % with respect to the rigid control.

Note that the experimental first Peak-to-Peak amplitude of residual vibrations is ≈10 times higher than the numerical one as shown in Fig. 13. This is due to the oscillatory behavior of the EE, which is a non-perfect point-mass that may swing out along the trajectory. The decrement of residual vibrations is used to define the performance indexes of control laws with respect to residual vibration reduction. It is defined as the ratio between the first and the fifth Peak-to-Peak amplitudes. Figure 14 shows the decrement of the Peak-to-Peak amplitudes of δp during the steady-state phase. This decrement is equal to 6.28 (2.94, 8.51, resp.) when the elasto-dynamic (elasto-static, rigid, resp.) model is used as a reference for the feed-forward. The decrement of the Peak-to-Peak amplitudes of δp obtained through experiments presents the same tendency than the numerical one. However, this experimental decrement is higher than the simulated ones. This means that some damping effects were not taken into consideration in simulations. A part of the damping may occur in the pulley system as the different control laws take into account neither the pulley friction nor the cable sliding into the pulley grooves. Furthermore, the winch friction torques are not well identified as the drums are 3D printed, leading to uncertainties. However, there is a clear improvement of trajectory tracking and vibration attenuation.

Conclusion

This paper proposed a model-based feed-forward control strategy for CDPRs.

The elasto-dynamic model of CDPRs was proposed to anticipate the full dynamic behavior of the mechanism including the vibratory effects, cable elongations and their interaction with the whole system to compensate them. The integration of tension distribution makes the proposed control scheme valid for every configuration of CDPRs. The comparison between the trajectory tracking errors obtained when using the proposed elasto-dynamic model and the trajectory tracking errors obtained while using the classical rigid and elasto-static ones as control references shows meaningful differences. The different experimental results confirm the numerical ones stating that the elasto-dynamic control leads to better trajectory tracking and better attenuation of residual vibrations of the EE than the conventional elasto-static and rigid controls. Numerical results mention that elasto-dynamic control leads to an improvement of 35 % for the first Peak-to-Peak amplitude and an improvement of 26 % for the fifth one with respect to the rigid control. Accordingly, experimentations confirm this improvement as shown in Fig. 15. Indeed, the elasto-dynamic control leads to an improvement of 50 % for the first Peak-to-Peak amplitude and an improvement of 32 % for the fifth one with respect to the rigid control. As a conclusion, the elasto-dynamic control scheme introduced in this paper is a good way to attenuate CDPR residual vibrations. The proposed control strategy was experimentally for a CDPR prototype where cable elasticity arises. Future research work should focus on the cable dynamic sagging modeling, which may be more relevant for other CDPR configurations. Such a case will occur when cable sagging becomes more influential than cable elasticity. F p = {P, x p , y p , z p } : End-effector frame.

Figure 1 :

 1 Figure 1: The ith closed-loop of a CDPR

Figure 2 :

 2 Figure 2: Feed-forward model-based PID control

Figure 3 :

 3 Figure 3: Equivalent control scheme: Discrete-time control

  q ref , ζ rg , ζ f and qrg are sent to the feedback control loop at each time α∆t b , α ∈ N + . These vectors are kept constant in the time interval [α∆t b (α + 1)∆t b ].

Figure 4 :

 4 Figure 4: (a) CREATOR prototype CAD diagram (b) End-effector desired path in F b

Figure 5 :

 5 Figure 5: Equivalent architecture of the CREATOR protype

Input

Figure 6 :

 6 Figure 6: 3D-printed pulleys used for the CREATOR protype

Figure 7 :

 7 Figure 7: End-effector and laser tracker target

Figure 8 :

 8 Figure 8: Desired (a) Position (b) Velocity and (c) Acceleration profiles of the end-effector

Figure 9 :

 9 Figure 9: (a) Position error of the end-effector along the z-axis (b) 2-norm of the end-effector Cartesian trajectory tracking error vector

Figure 10 :

 10 Figure 10: Experimental setup: Non-suspended initial pose of the end-effector

Figure 11 : 3 m

 113 Figure 11: Experimental results: Current set-points (a) i 1 m , (b) i 2 m and (c) i 3 m Figure 11 shows the set-point of current i m = [i 1 m , i 2 m , i 3 m ] T sent to the actuators, which is an image of torque set-point ζ m . It is noteworthy that the elasto-dynamic current set-point takes into account the oscillatory behavior on the contrary to classical feed-forwards.As one can see, the set-point of current is not constant during the steady-state phase. This reflects the integrator effect of the PID controller as its input is never

Figure 12 :

 12 Figure 12: Experimental results: Position error along (a) x-axis, (b) y-axis and (c) z-axis of the end-effector (d) Position error norm

Figure 13 :

 13 Figure 13: (a) Experimental, (b) Numerical First and fifth Peak-to-Peak amplitudes of δp

Figure 14 :

 14 Figure 14: Decrement of the Peak-to-Peak amplitudes of δp during the steady-state phase

Figure 15 :

 15 Figure 15: Vibration improvement of the elasto-dynamic control with respect to rigid control

F

  b = {O, x b , y b , z b } : Base frame.

2

  Parker ™ SMB/SME Brushless servo motor series. Ref: SMEA 60601-489IZ64S54.

3

Parker ™ Economical Planetary Gearheads -PE. Ref: PE3-008-10M040/063/09/20.

Table 1 :

 1 Characteristics of the CREATOR actuators and the gains of the PID controller

	Parameter (Unit)	Imotor (kg.m 2 )	I reducer (kg.m 2 )	I (kg.m 2 )	ke (Vrms.s/rad) (Ohm)	r
	Value	42.5e-6	6.5e-6	0.0031	0.48	12.8	8

Table 2 :

 2 Kp, K d and K i gains of the PID controller

	Gain	Kp	K d	Ki
	Value 1125.8 58.12 7269.60

Table 3 :

 3 Experimental data: Static deflection of the end-effector before compensation

	Experience	Test 1 Test 2 Test 3
	Static deflection dp (mm)	9.00	6.32	8.30

first move of the actuator corresponds to the dry friction torque ζ d . As we

ControlDesk is the dSPACE experiment software for seamless electronic control unit development. It performs all the necessary tasks in a single working environment, from the start of experimentation right to the end.

"rigid control" refers to "rigid model-based control of CDPR"

"elasto-dynamic control" refers to "elasto-dynamic model-based control of CDPR"

τ : Cable tension vector. q : Angular displacement vector. l : Cable length vector.

χ : Diagonal matrix containing the winding ratio of the winches.