Virtual Sensor for Kinematic Estimation of Flexible Links in Parallel Robots

The control of flexible link parallel manipulators is still an open area of research, endpoint trajectory tracking being one of the main challenges in this type of robot. The flexibility and deformations of the limbs make the estimation of the Tool Centre Point (TCP) position a challenging one. Authors have proposed different approaches to estimate this deformation and deduce the location of the TCP. However, most of these approaches require expensive measurement systems or the use of high computational cost integration methods. This work presents a novel approach based on a virtual sensor which can not only precisely estimate the deformation of the flexible links in control applications (less than 2% error), but also its derivatives (less than 6% error in velocity and 13% error in acceleration) according to simulation results. The validity of the proposed Virtual Sensor is tested in a Delta Robot, where the position of the TCP is estimated based on the Virtual Sensor measurements with less than a 0.03% of error in comparison with the flexible approach developed in ADAMS Multibody Software.


Introduction
Since their introduction in the industry back in the early 1960s, robots have been considered as a cornerstone of the mass-production system due to their capacity to combine minimal cost with adaptability, quality and high productivity. Although automotive industry has been the predominant user, in the last decade there has been a growing interest on advanced robot technology in other areas such as food, pharmaceutics or manufacturing industries.
In the current global market, smaller production time and higher quality products are required to be competitive. Robotic applications in industry need to be fast and accurate enough to satisfy these requirements. However, in serial robots, which have been traditionally used in industry, the increase in speed usually implies a loss in accuracy. Hence, when both requirements need to be met, the use of Parallel Kinematic Robots (PKR) [1] has been proposed. These robots are composed by two platforms connected by multiple kinematic chains, being one of them fixed (base platform), and the other, where the Tool Centre Point (TCP) is located, mobile.
Parallel robots have a series of advantages in comparison with their serial counterparts, such as higher load/weight ratio, and higher stiffness and precision. Furthermore, since the motors can be fixed to the base platform, the moving mass is reduced and higher acceleration and speed can be achieved. These characteristics improve the performance of the robot in terms of working efficiency, precision and energy consumption [2].
To address the need for achieving higher productivity, manufacturers have tried to reduce the moving mass of the robots by decreasing the cross-section of the limb. This, in addition to the high accelerations required for low operation cycles, results in a certain degree of elastic deformation in the limbs [2]. These deformations can cause substantial errors in motion control due to their influence in both dynamics and kinematics of these robots. In fact, one of the most difficult problems caused by elastic deformation is TCP trajectory tracking, as deformations need to be compensated to maintain the required accuracy.
In stiff or undeformable robots, the location of the TCP can be determined using the active joint position sensors, this is, the sensors attached to the actuators, and the kinematic model derived from its geometrical structure. However, in the case of flexible robots, the TCP location depends not only on joint positions but also on the deformation of its flexible links. Therefore, the stiff element assumption does not provide enough accuracy, as the deformation of the different flexible elements must be considered [3].
Control of flexible link manipulators, even single link manipulators, is a challenging task as these robots are of non-minimum phase [4]. When multiple links are considered, this problem increases in complexity due to the distributed flexibility. In addition, feedback controllers require the measurement or estimation of both joint motions and the flexible link deformation to guarantee accuracy and trajectory tracking while reducing the effect of vibrations caused by flexibility [4]. For this purpose, both rigid and flexible variables must be measured.
The measurement of deformations due to link flexibility has been carried out considering different types of sensors [5]. One of the approaches that provides better deformation data is the use of artificial vision, which has been used to determine the location of the TCP of flexible manipulators [6][7][8] and provide their calibration [9,10]. Vision systems can be used to implement visual servoing [11], that controls directly the TCP's position of the robot without the use of estimators. Most of visual servoing approaches require the calculation of the image Jacobian and robot Jacobian to map end-effector velocity to image feature velocity [12]. In flexible manipulators, however, these Jacobians require the information of the flexible variables [13], which are still needed to be measured [14]. Hence, the measurement of the deformation of each link by the use of artificial vision systems has also been proposed [6,7], though this approach requires a camera per flexible link, increasing the cost of the robot. Hence, although precise and global data can be derived from vision systems, their use is limited by their view range, the effects of visual obstruction, link interference, low sampling rate and the computational cost required to process the images.
To overcome the limitations of vision systems, conventional sensors such as accelerometers or strain gauges are widely used. The accelerometers are usually placed at the tip of each flexible link, and their measurements are traditionally used to implement vibration control approaches in order to stabilize the oscillation of the TCP due to flexibility [15]. Hence, as tip accelerations contain both information of rigid-body motions and flexible-link vibrations, it is possible to estimate the TCP location of the manipulator if joint motions are sensorized [16]. This approach, however, provides noisy measurements, contains biases, and requires to integrate the acceleration signal twice for velocity and position estimations, resulting in high accumulation of errors [6].
On the other hand, strain gauges have been widely used to measure the local deformation of flexible links [16]. For instance, in [17] strain gauges were used in a three degree-of-freedom flexible robot links for control proposes. In [18] an interpolation algorithm to determine the tip position and orientation of a flexible beam is defined from a finite set of strain measurements. Finally, in [14] the bending and torsional deformations kinematic relations of a 3D flexible-beam were verified using strain measurements.
Strain gauges require only the knowledge of geometrical parameters of the link to estimate the deformation, without considering dynamic parameters such as link masses or inertias. Hence, they have become an attractive approach. However, strain gauges are prone to temperature variations and they inherently suffer from measurement noise [19] and biasing due to electromagnetic interferences [6,7].
Thus, the measurement of the deformation of flexible links is not a trivial task even with the use of additional sensors. Based on this fact, the use of sensor fusion techniques have also been proposed to improve the measurement: In [20] the fusion of accelerometer and encoder signals using a disturbance observer to compensate the nonlinearities of the deformation was proposed, while in [21], an extended Kalman filter to process accelerometer and encoder data was developed to estimate the forward kinematics of a 6 DOF robot. However, although the use of sensor fusion techniques improves the accuracy of the measurement, the proposed techniques usually require high computational load [22].
This work presents a novel Virtual Sensor that allows to estimate the deformation of flexible links based on the use of a single high resolution optical encoder and a mathematical model. This deformation is critical to accurately estimate the position and orientation of the TCP of flexible parallel robots. The proposed approach allows to calculate the deformation of the flexible links with minimal computational cost, providing significant advantages for Real-Time implementation over the aforementioned approaches.
The rest of this paper is organized as follows: Section 2 outlines concept and theoretical development of the presented virtual sensor. In Section 3, the theoretical development is applied into a Delta robot. The simulation results are discussed in Section 4. Finally, the most important ideas are summarized in the final section.

Theoretical Development of the Kinematic Virtual Sensor
The Direct Kinematic Model of parallel robots allows to calculate the trajectory in the operational space x(t) for any given joint space trajectory q(t), this is, f(q(t)) = x(t). Therefore, it is easy to see that for any flexible link manipulator, the direct kinematic equation admits multiple solutions, as the deformation of the link provides more than one end-effector position/orientation for the same configuration of joints q [23] (Figure 1). Figure 1. A flexible manipulator with the same configuration of rigid joints q r (q 1 , q 2 ) but multiple end-effector position/orientation due to the flexibility of the links, q f (q f d , q f s ). (q f s link tip transverse flexural deflection and q f s link tip flexural slope). Therefore, the Direct Kinematic Model, which is mandatory for robot control applications, has to be defined in terms of the rigid joint motion variables q r and the flexible variables q f that model the deformation of the flexible links, where the flexible variables q f (q f d , q f s ) are composed by the link tip transverse flexural deflection q f s and the link tip flexural slope q f s . However, the number of degrees of freedom (DOF) of the system increases, and in order to estimate correctly the TCP of the robot, measurements of both rigid and flexible variables are required.
Measurement of rigid joints can be easily carried out using extra sensors attached to rotary or linear joints. However, as stated in the introduction, the measurement of the flexible variables can be a challenging tasks. In this section, a novel Virtual Sensor based on a normal mode analysis approach is proposed that allows to estimate the deformation with enough accuracy and minimal computational cost.

Fundamentals of the Virtual Sensor
The normal mode analysis is a mathematical tool which represents a pattern of motion in which all parts of a system move with the same frequency and fixed phase relation. Hence, this mathematical tool provides the kinematic relation between the different DOF of the system based on its dynamic properties, such as their material, structure and boundary conditions. Furthermore, the general motion of a system is composed as a superposition of its normal modes, being all normal modes orthogonal to each others, since the excitation of one mode will not affect to a different mode. This way, if a particular mode is considered and the modal motion of a single degree of freedom measured, the complete set of DOF q DOF modal motions can be estimated, where X k is the kth eigenvector or mode shape and k is the kth modal motion [24]. Several researchers [25][26][27] have suggested considering only the first few modes in the model neglecting high frequency since the amplitude terms related to them are much smaller. In addition, due to the damping properties of the materials, high frequency vibrations are softened more quickly than those of lower frequency [28]. On the other hand, in control applications of industrial manipulators, the bandwidth of the working frequencies is limited by the actuator system and the application itself. This, in addition to the use of rigid materials in the construction of robots, usually ensure that the second resonant frequency of the manipulator is out of the working bandwidth in all robot workspace. Hence, the deflection of the limbs can be estimated considering the mode associated to the first natural frequency of the manipulator. This way, Equation (2) can be simplified to, This constitutes the basis of the proposed Virtual Sensor. Let us assume that a robot presents a series of flexible links i = 1, . . . , n connected with several stiff links. If the modal analysis of each link is carried out and the previous facts considered, the final deflection of the ith link can be estimated with a properly placed sensor that measures a single flexible DOF and the use of the model obtained after carrying out the modal analysis. Furthermore, if only the link is considered, the relation that is obtained in the modal analysis is constant for a given mode, resulting in a low computational cost approach for the estimation of the deflection of the link.
Consider, for instance, the example of Figure 2, where the flexible link b is connected with two stiff links (a, c) using rotary joints. If Finite Element Method approach and Euler-Bernouilli Beam Theory is used to model the flexible link, one of the flexible variables, q f s i , would be the deformation slope at the tip of the flexible link. If a single high resolution optical encoder is used to measure this angular deformation and its data is introduced in Equation (3), by substituting 1 = q f s i , the total deflection at the tip of the link could be estimated, which is required to compute the TCP estimation of the robot. This approach presents several advantages over previous approaches. First, as only the link is considered, the relation that is obtained in the modal analysis is constant for a given mode, resulting in a low computational cost approach for the estimation of the deflection of the link. Second, if small deflections are assumed, the approach provides the required accuracy for control purposes. Finally, since vector X k defines the relationship between the kth mode of the motion of the flexible DOFs of the ith link, the relationship holds not only for the position problem, but also for velocity and acceleration problems.
This way, full kinematic relations between the flexible variables, Note, however, that the whole manipulator is not considered in the modal analysis, as it is highly dependant on the particular configuration of the robot. The consideration of the whole robot requires recalculation the modal analysis, and thus, high computational cost. Hence, the proposed approach presents a simplification of the procedure by considering only the flexible joints, q DOF = q f , allowing high computational efficiency while maintaining the required accuracy. In this work, the normal modes of each link are independent of the manipulator configuration, and only depend on the mechanical properties of the limbs and the discretization carried out based on the Finite Element Method approach. This allows to calculate the required matrices off-line, reducing the time to compute link deformation.
In the following section, the mathematical development of the presented virtual sensor is widely detailed.

Modelling of Flexible Links
In order to carry out the modal analysis of each flexible link and calculate their eigenvalues, the generalised inertia matrix and the stiffness matrix of each link has to be calculated. In this section, the procedure to calculate these matrices is detailed.
Finite Element Method (FEM) [29][30][31] and Assumed Modes Method (AMM) [32][33][34] approaches are the most used ones when modelling flexible links robots for control purposes. As the requirements and procedures of both approaches are different, several studies have been carried out to compare them [28,35,36]. This way, when computational cost is to be analysed, the aforementioned studies determine that FEM is a better approach due to its fewer computation requirements.
The Finite Element Method considers each flexible link i as an assemblage of a finite number, n i , of small elements of length l i which are interconnected at certain points called nodes. Each element is referenced as ij, where subscript j denotes the number of the element. As it is well known, the more number of elements per link (hence, smaller elements), the more precise the overall solution of the system will be, making it converge to the exact solution as precisely as desired at the cost of higher computational cost.
In order to solve the dynamic problem, the Euler-Lagrange method is used, in which the energy of each element is to be analysed. However, before calculating the kinematic and potential energy, the position vector r 0 i , in inertial coordinates, of each element must be defined ( Figure 3). For simplicity, this vector will be expressed in terms of a local coordinate vector r i [3]. This is accomplished by using the transformation matrix T i 0 , which relates the location of the reference system attached to each link 0 i X i Y i Z i and the inertial system OXYZ, where According to Euler-Bernoulli beam theory, all elements are considered to possess two degrees of freedom (DOF) at each end of the element: a transverse flexural deflection (u 2j−1 and u 2j+1 ) and a flexural slope (u 2j and u 2j+2 ) ( Figure 3). These flexible DOF are related by the use of the shape functions φ k (x), which describe the flexural displacement z, as where, for element j of the ith link the shape functions are defined by the following Hermitian polynomials [37,38], where l ij is the length of the element and the displacement variables x and z(x) satisfy the following boundary conditions, Once the position vector is defined, the kinetic energy T ij and potential energy V ij of each element are computed in terms of the generalised variables of the system q = (q 1 , q 2 , ..., q n ) T and their time derivativesq [39]. These energies are then summed to obtain the total kinetic T and potential V energies for the entire system. That is, and where m is the number of links in the system and the kinetic energy T ij is obtained by integrating over the ij element's length, l ij , the corresponding energy function. Thus, which can be rewritten as, where matrix M ij is the generalised inertia matrix of the element ij. This matrix is symmetrical and each (k, o) element of it is defined as, where ρ i is the mass density, A i is the cross-section area of the element, z jk is the kth element of z j = [q r , u 2j−1 , u 2j , u 2j+1 , u 2j+2 ] T , n q is the number of the variables of the system and q r = [q 1 , q 2 , ..., q n r ] is the vector of variables associated to the stiff model. The resulting generalized inertia matrix M ij with respect to a single rotation, q r = [q 1 ], is always defined as [3]: and, with, where ρ ij is the mass density, A ij the cross-section area and l ij the length of the j th element of the link i.
In the same way, by computing the potential energy for each element V ij of the link and adding their contributions, the overall potential energy of the system is obtained, Equation (11). However, unlike the kinetic energy, the potential energy is divided into two components: the potential energy due to the gravity (V g ij ) and the potential energy due to the elasticity (V e ij ). This is, By algebraically manipulating the elasticity term of Equation (19), the following expression is obtained, where K ij is the stiffness matrix of the j element of the i th link which is calculated as, and 0 n r ×n r is a n r × n r dimension null matrix, n r is the number of rigid variables of the system, E ij is the Young's modulus of the element and I i is the moment of inertia of the link i. In order to determine the deflection of each link, the generalised inertia and the stiffness submatrices concerning the flexible variables of each link (M i f and K i f respectively) have to be defined. As seen before, the generalised inertia and stiffness of each link is obtained by adding the inertia M ij and the stiffness K ij of all the elements that compose each link. From these matrices, the submatrices associated to the flexible variables q f are selected.
As length l ij and the Young's modulus E ij of each element are the same for all elements of a given link, the structure of the obtained matrices is the similar, as it can be seen in the following example for a three element link, and where Once the generalised inertia submatrix M i f and the stiffness submatrix K i f of each limb are obtained, the boundary conditions of each flexible link have to be defined, this is, the deflection, the slope, the shear force and/or the bending moment in both ends of each link (x = 0 and x = l T i ). Depending on the link's configuration, different Boundary Conditions (BC) are defined. As an example, when the link is working as a simply supported beam, both displacement and bending moment are set to zero at both ends. In a clamped configuration, on the contrary, only displacement and the slope are set to zero at the clamped node.

Modal Analysis of the Flexible Links
Once calculated the inertia and stiffness submatrices, the modal analysis procedure can be applied. For that purpose, a reduced set of equations of motion is required. In these equations the damping and the applied load are not taken into account, To solve Equation (24), a harmonic solution as shown in Equation (25) is proposed.
where ω is the natural frequency vector. The harmonic solution is not only required for the numerical solution of the problem, but also for the physical interpretation of the equation. This way, the harmonic solution defines the way in which all the degrees of freedom of the link deflect and their relationship. Therefore, for a given mode, the structural configuration does not change during motion, only its amplitude, being possible to relate the effects of flexibility on each DOF if one of them is measured. This concept is the key to develop the proposed virtual sensor, as stated in the beginning of Section 2.
By substituting Equation (25) into Equation (24), the following equation is obtained.
which after simplifying becomes Equation (27) is called the eigenequation, which is a set of homogeneous algebraic equations for the components of the eigenvector, which forms the basis for the eigenvalue problem. The basic structure of the eigenvalue problem is defined as, where A is a square matrix, λ are the eigenvalues, I is an identity matrix of the same size as A and X is the eigenvector.
In the structural analysis of the limb of a robot, as in any other mechanical structure, the use of the generalised inertia matrix M i f BC and the stiffness matrix K i f BC gives a physical meaning to the eigenequation solution, being the obtained eigenvalues the natural frequencies of the system (λ = ω 2 ) and the obtained eigenvector X the mode shapes of the system. Solving Equation (27), two mathematically possible solutions can be obtained: This is a trivial solution, which does not provide any valuable information from a physical point of view, since it represents the case of no motion. 2. If X = 0, the result of the eigenproblem is reduced to solve Equation (30).
The determinant defined in Equation (30) can only be zero at a set of discrete eigenvalues ω 2 i . Furthermore, there is an eigenvector X which satisfies the equation Equation (27). Therefore, Equation (27) can be rewritten as: where the number of eigenvalues and eigenvectors n DOF is equal to the number of DOF of the discretised limb has. Moreover, the ith eigenvalue ω 2 i is related to the ith natural frequency, f i , by Hence, by solving Equation (31), the eigenvector X i can be calculated, which, in fact, represents the relationship between the flexible DOF of the link as it has been shown in Equation (4). As stated in the introduction of Section 2, in the proposed procedure, this calculation can be carried out off-line, as the inertia and stiffness submatrices do not depend on the particular position of the robot. This way, once calculated the eigenvalues for a given mode, if one of the flexible DOF is measured, the rest can be estimated.
Note that this approach is based on a simplified flexible model based on the discretization of the flexible links. However, its simple structure makes it computationally efficient, and, as it will be proved in the next section, it allows to estimate the deformation of flexible links with accuracy.

Case of Study: Delta Robot
In the following section, the developed virtual sensor will be implemented into a flexible-link parallel robot in order to validate it. For that purpose, the procedure detailed in the previous section will be followed.

The Delta Robot
The Delta robot is one of the most popular three DOF parallel robot, being widely used in industry. Its lightweight structure provides dynamic capabilities for quick motions, being mainly used for pick and place applications, in which several models can achieve up to 300 picks per minute.
The robot is composed of 3 arms (i = 1, 2, 3) distributed uniformly (β 1 = 0 rad, β 2 = 2 π/3 rad, β 3 = 4 π/3 rad) that connect the fixed based with the mobile platform. Each of these arms is composed by two limbs: the upper one is connected to the fixed base by actuated rotational joints, this is, they are connected to the actuators (q a i for i = 1, 2, 3); the lower ones, on the other hand, are composed by an articulated parallelogram which allows passive rotations in two directions (q na i and α i for i = 1, 2, 3), as it can be seen in Figure 4. This configuration limits the motion of the mobile platform, and therefore the motion of the TCP, to a pure three dimensional translational motion with no rotations.  For this study case, the main parameters have been obtained from a commercial Omron Mini Delta CR-UGD4MINI-NR robot (Omron, Kioto, Japan), whose model has also been implemented in ADAMS Multibody Software (2014.0.1, MSC Software, Newport Beach, CA, USA). These parameters are summarised in Table 1.  I L iyy = 8.7001 × 10 −5 I l iyy = 0.0010 I p yy = 4.6225 × 10 −5 I L izz = 8.6422 × 10 −5 I l izz = 0.0010 I p zz = 9.1472 × 10 −5 In order to increase the effect of deformations, the links which compose the lower limbs have been replaced by AW5083/H111 aluminium platens (IMH, Elgoibar, Spain). These platens have 0.003 m thickness and a width of 0.015 m, and the used aluminium presents a Young's modulus of E i = 71 GPa and a mass density of ρ i = 2740 kg/m 3 . The geometry has been selected in order to limit the deformation to the direction of the z i axis of each link, and analyse the validity of the proposed approach. All these proprieties are summarised in Table 2.

Numeric Implementation of the Kinematic Virtual Sensor Applied to the Delta Robot
Based on the procedure detailed in Section 2 the first step to define the virtual sensor for the estimation of the links deflection is to define the transformation matrices between the inertial system OX 0 Y 0 Z 0 , and the body-fixed system O i4 X i4 Y i4 Z i4 , i = 1, 2, 3 ( Figure 5). In order to place these reference systems, a systematic procedure is used. The most used procedure in parallel manipulators is the Equivalent Rigid Link System (ERLS) description based on the notation of Khalil and Kleinfinger [40] which is an adaptation of the well-known Denavit-Hartenberg notation [41] for closed-loop robots. However, this procedure has been adapted to provide a more systematic approach in flexible robots. This way, the x axis is always set along the length of the link, the z axis along the deformation plane of the link and the y axis perpendicular to the xz plane. Hence, the transformation matrix T i 0 for each kinematic loop of the Delta Robot is defined as, In order to simplify the mathematical model, the articulated parallelogram has been reduced to a single equivalent link. This assumption has proven to be valid after analysing the dynamic behaviour of the articulated parallelogram link in simulations, in which it was observed that both elements of the parallelogram present the same dynamic behaviour at lower natural frequencies.
The Finite Element Method is then applied to model the deformation of the links. This approach requires to define the number of elements considered per link. Tsujisawa suggested in [42] that for robot manipulators a relatively small number of modes (two or three) are usually enough to represent the dynamics of flexible links. Furthermore, according to Przemieniecki's work [43], when an m-element FEM model is used, the system's first m vibration modes can be obtained with acceptable accuracy. Hence, for the developed virtual sensor a three element FEM model per flexible link will be used.
Once the number of elements per link is defined, the generalised inertia matrix and the stiffness matrix of each flexible link has to be calculated. As explained in Section 2.2, for the developed approach, just the submatrix concerning the flexible variables, M i f and K i f respectively, are required. These submatrices have the main advantage that are constant, since they are calculated in terms of the mechanical parameters, which are invariable in all the elements of a given link.
Hence, for a three-element FEM model platens, and the stiffness submatrix as where l = l i = l T i 3 . Finally, the boundary conditions (BC) to each flexible link must be defined. It is a general mathematical principle that the number of BC necessary to calculate the solution to a differential equation matches the order of the differential equation. The used Euler-Bernoulli beam model is a fourth-order equation, so each flexible link requires four boundary conditions to be solved.
According to the literature, different types of boundary conditions can be used. The BC associated to a Pinned-Pinned beam model, for instance, has been used for trajectory control of the TCP of the robot [44] due to its simplicity when implementing the position problem of the TCP. The boundary conditions associated to Free-Free beams have also been used with rigid motion variables [45]. Finally, other authors, have considered the flexible link as a Clamped Beam, demonstrating that their use simplifies the measurement of joint variables and the calculation of required torques [46,47].
According to the experimental verifications reported in [26,48], if the beam to hub inertia ratio is very small, (at least ten times smaller) the Clamped Beam boundary condition set yields better results compared to Pinned-Pinned boundary condition set. Hence, for the developed virtual sensor, Clamped Beam boundary conditions will be used:

•
The base of the link does not experience any deflection: z i (0) = 0.

•
The link at the base has no deformation, so that the derivative of the deflection function is zero at that point: There is no bending moment at the end of the link: • There is no shearing force acting at the end of the link: If these boundary conditions are introduced in the developed model, Therefore, the first two lines and columns of the M i f and K i f can be erased, defining the new generalized inertia submatrix as, and the stiffness submatrix as Once inertia and stiffness submatrices have been defined, based on Equation (31) the natural frequencies ω i and the eigenvector X i can be obtained, As explained in Section 2.3, the obtained eigenvector X i defines the relationship between the displacement and the slope of the nodes of the link. The same obtained eigenvector, not only holds for the position problem relationship, but also for the velocity and acceleration problems relationship, which makes the base of the proposed virtual sensor.
Hence, Equation (39) allows to calculate the modal displacement or slope of all the nodes of the ith link using the measurement of a single displacement or slope. Moreover, the measurement of the slope of the flexible link's tip (u i8 ) is quite straightforward in parallel robots ( Figure 6) by the use of precision encoders. Hence, using a single sensor, the deformation (q i f ), velocities (q i f ) and acceleration (q i f ) of the flexible link's DOF can be estimated easily using the virtual sensor Equation (39).

Simulation Setup
The virtual sensor is validated next. For that purpose, the Delta robot detailed in Section 3.1 has been modelled in ADAMS Multibody Software (2014.0.1, MSC Software, Newport Beach, CA, USA), (Figure 7). Flexibility parameters ( Table 2) have also been introduced into the software. In order to analyse the flexible behaviour of the links of the Delta robot and validate the proposed approach, the robot is excited in a wide range of frequencies within the operational bandwidth of the robot. Hence, a sinusoidal motion of different frequencies has been applied in the actuated joints q a i , where ω τ 1 = 3 rad/s, ω τ 2 = 0.75 · 3 rad/s, ω τ 3 = 1.2 · 3 rad/s and t is the time, which is defined between 0 and 5 s. With it, the obtained TCP movement if shown in Figure 8. Three additional sensors will be considered in the robot. As stated in the previous section, the virtual sensor can determine the deformation of the link by measuring a single deformation slope, which in the case of rotary joints can be easily measured by the use of three encoders installed in the passive rotary joints of the mobile platform ( Figure 6). The validation of the proposed kinematic virtual sensor has been carried out considering two factors. First, the evaluation of the virtual sensor capability to estimate the deformation of the flexible links based on a single encoder measurement will be analysed. According to FEM definition, each of the flexible links is composed by a certain number of nodes, each of which have two DOF, the displacement and the slope deflection. Hence, in order to analyse the accuracy of the virtual sensor, the estimation of each DOF is compared with the measurements obtained from ADAMS Multibody Software. As explained in previous sections, the developed virtual sensor not only estimates the position and the orientation of the DOF, but also the speed and the acceleration of them. Due to it, the accuracy of the virtual sensor will be analysed in position, speed and acceleration.
On the other hand, the performance of the virtual sensor will be evaluated for its application to the calculation of the Direct Kinematic Problem of flexible parallel robots. This model is used to estimate the motion of the TCP in terms of the joint sensors of the robot and the virtual sensors introduced to measure the flexibility. The kinematic model is derived from the loop closure equations of the mechanism, which relate the kinematic variables of the robot, such as the angle of the joints, the position of the TCP and the deformations due to flexibility [1]. For the particular case of the analysed Delta robot, these equations are (Figure 9), where i = 1, 2, 3 is the link identification and δ i is the deflection distance in the flexible link end which respect to its analogous rigid link. This last variable is the one to be estimated by the virtual sensor. Hence, the position of the TCP of the robot, p x is calculated as, Next, the results obtained from both analysis will be discussed.

Results and Discussion
First, the accuracy of the virtual sensor to estimate the deformation of each flexible link will be analysed. For that purpose, the deflection estimation of a single flexible link tip provided by the virtual sensor has been compared with the data provided by ADAMS Multibody software. Results are shown in Figure 10, where it has been assumed that the first natural frequency amplitude was dominant in the estimator. As it can be seen in Figure 10, the estimation is accurate. Data shows that the virtual sensor is able to estimate the deflection of the tip with a maximum error of less than 0.1 mm, which is less than the 2% of the maximum deflection amplitude.
Furthermore, the proposed virtual sensor does not only establish the relationship between the displacements and the orientations of each link's nodes, but also their speed and acceleration. Hence, following the same procedure, the speed error of the flexible link tip, Figure 11a, and its acceleration error, Figure 11b, can be obtained.  Figure 11a shows that the maximum error is 3 × 10 −3 m/s, which means a speed error of less than the 6% of the amplitude of the signal. In the same way, in Figure 11b the acceleration error is plotted, in which the maximum value is 0.2 m/s 2 , a relative error of less than the 13% of its amplitude.
Next, the deflection data calculated by the virtual sensor is used to compute the Direct Kinematic Problem (Equation (41)). For that purpose, the estimation error of the TCP considering the virtual sensor is compared with the TCP trajectory shown in Figure 8, which was obtained from the flexible ADAMS Multibody Software model, as it is shown in Figure 12.   This way, in Figure 13 the x, y and z component of the TCP positioning error of the developed estimator is shown. As it can be seen, the estimation error of the developed approach rises up to 10 −4 m for the x axis, 5 × 10 −5 m for the y axis and 1.5 × 10 −5 m for the z axis, which means an error of less than the 0.03%, 0.02%, 0.007% of the applied movement, respectively.  Hence, the aforementioned results demonstrate that the developed approach presents great potential for the estimation of deflections in flexible robots. However, it is important to note that in order to obtain accurate estimation results the proposed virtual sensor needs to consider the bandwidth of the motors and ensure that the range of the estimator is inside it. Furthermore, an analysis of the resolution of the encoder needs to be carried out so that the precision requirements of the estimation are satisfied. Finally, although implementation issues are not the focus of this work, it is important to analyze the computational cost requirements of the approach when applied to a particular robot, and select an appropriate hardware to ensure that the estimation is carried out within the time constraints of the control loop.

Conclusions
Parallel robots are widely used in industry for high dynamic requirement tasks that imply high speed. In order to achieve these capabilities, the structure of the robot is usually lightweight, so that small deflections can arise due to high speed motions. Modern controllers can compensate the deformations caused by the elastic deformation of the links. However, the measurement of small deformations is a challenging tasks that usually require extra sensors and complex models.
In this work, a novel virtual sensor for the estimation of the deformation of flexible links in parallel robots is presented. The developed approach is based on the modal analysis of the flexible link, which is modelled based on Finite Element Method approach and the Euler Bernouilli Beam theory. This approach is focused on control applications, providing a simple yet computationally effective approach that provides accurate estimations in comparison with other approaches.
The proposed virtual sensor has been validated in a flexible Delta Robot which has been modelled on ADAMS Multibody Software. The deformation data obtained in ADAMS was compared with the one estimated by the virtual sensor, resulting in a maximum error of less than 2% for the deformation estimation and 6% and 13% for the speed and acceleration estimation, respectively.
In addition, using a similar procedure, the virtual sensor estimation was used to estimate the Direct Kinematic Problem of the flexible Delta Robot. This data was compared with the flexible Delta model developed in ADAMS Multibody Software, resulting in a TCP error estimation of less than 0.03% (0.03%% for the x axis, 0.02% for the y axis and 0.007% for the z axis). Author Contributions: Pablo Bengoa is in charge of the modelling and control of flexible parallel robots. During his work, he has developed the kinematic virtual sensor approach presented in this paper. Asier Zubizarreta and Itziar Cabanes have supervised and coordinated all the research work. They also have provided technical support in the modelling of the robot whereas Charles Pinto has added technical advice on vibration analysis. Aitziber Mancisidor and Sara Mata have participated in the simulations done in Matlab and ADAMS Multibody Software, respectively. All authors contributed to the organization of the paper and approval of the final manuscript.

Conflicts of Interest:
The authors declare no conflicts of interest.

Abbreviations
The following abbreviations are used in this manuscript: