End-effector trajectory tracking of flexible link parallel robots using servo constraints

We apply the concept of servo constraints to end-effector trajectory tracking control of parallel robots with structural link flexibilities. Such servo constraints deliver the inverse robot model where solution approaches via projections are proposed, which transform the resulting differential-algebraic equations to ordinary differential equations. The applicable solution process depends on the existence and stability of the internal dynamics. When using the exact end-effector of flexible link robots as output, this internal dynamics is usually unstable. Then a two-point boundary value problem is considered in the framework of stable inversion to obtain the noncausal solution offline. This solution is used as a feedforward control, which is initially combined only with actuator feedback control. To also account for errors within the link flexibility, the well-known linear–quadratic regulator is adapted to end-effector trajectory tracking based on differential-algebraic equations. Finally, we propose a systematic input–output feedback linearization approach, which uses servo constraints for flexible link parallel robots. Here a minimum phase system is obtained by tracking a redefined end-effector output, which is an approximation of the exact end-effector position. All control concepts are validated experimentally with a parallel robot having a highly flexible link. The results allow us to compare different control approaches and show the superior performance of controllers that rely on a flexible multibody model in contrast to classical rigid multibody modeling.


Introduction
In the last years the importance of energy and material efficiency has increased substantially. In the field of robotic manipulators, these challenges can be overcome by lightweight designs. Such designs also provide more safety for human-robot interaction, which is a R. Seifried robert.seifried@tuhh.de 1 considerably growing application area. However, these advantages come at the expense of oscillations due to structural flexibilities. If these flexibilities are significant, then they need to be considered within the control design. This is often a complicated task as such systems are in general underactuated, i.e., they have less control inputs than degrees of freedom. Consequently, many control approaches from rigid-body robotics cannot be directly applied. Additionally, such flexible robots are often nonminimum phase systems when using the exact end-effector as output, which makes the control even more challenging. It should be noted that here the description of the end-effector as exact differentiates it from later introduced end-effector approximations.
In this paper, we consider such structural flexibilities in the form of flexible link robots, with focus on the end-effector trajectory tracking control. This topic has been extensively studied in the last decades [2,11,22]. Still, a significant part of the literature reports only theoretical and numerical results or considers only single-link flexible manipulators [29] or serial robots. Especially, the literature on flexible link parallel robots is rare [16,35]. Thus this research aims at contributing to this field by presenting a complete process from modeling over end-effector trajectory tracking control design to the experimental validation with a flexible two-link parallel robot.
Usually, parallel robots are described by differential-algebraic equations (DAEs), and serial robots by ordinary differential equations (ODEs). Thus, mathematically, serial robots can be regarded as a simpler subgroup of parallel robots. Therefore the developed techniques can be also applied to purely serial robots with flexible links.

Modeling
The occurring link flexibilities are modeled with the linear finite element method (FEM), and a modal reduction [26] is used to reduce the number of elastic degrees of freedom. Based on the floating frame of reference approach [32], the motion of a flexible body is separated into a large nonlinearly described reference frame motion and linear deformations described with respect to this reference frame. This yields a computationally efficient model for state estimation and control.
Still, geometric nonlinearities such as the foreshortening effect [21] can have a large influence on the actual end-effector kinematics. In this research, we efficiently describe geometric nonlinear effects by solely adapting the kinematics of the considered output.

End-effector trajectory tracking
The discussed modeling approaches are the basis for the main focus of the paper. This is the end-effector trajectory tracking control for flexible link parallel robots based on an inverse model using servo constraints. These servo constraints [3] are additional algebraic equations, which restrict the output to a desired trajectory. To simplify the solution for the considered model inversion-based controllers, we transform the obtained DAEs to ODEs based on projections. As mentioned before, when taking the exact end-effector of flexible link robots as output, the zero dynamics is usually unstable, i.e., the system is nonminimum phase; see, e.g., [30]. To still obtain bounded results when tracking the exact end-effector, we use the framework of stable inversion [9], where a two-point boundary value problem (BVP) needs to be solved offline. Alternatively to stable inversion, we use the concept of output redefinition to obtain a stable internal dynamics. This allows us to get bounded results by solving an initial value problem (IVP) via forward integration instead of a usually much more complicated BVP. The IVP may even allow a model inversion in real-time. Here we realize the redefinition by directly weighting the elastic deformation in the nonlinear end-effector output function based on [24].

Linear-quadratic regulator
Within initial experiments the feedforward control approaches based on model inversion are not combined with feedback of the flexibilities but only with feedback control of actuatorrelated quantities.
Still, the flexible model inversion computes the complete system state, including elastic coordinates, which will be used within a linear-quadratic regulator (LQR) applied to endeffector trajectory tracking. In [1,8,27], LQRs are applied experimentally to end-effector trajectory tracking of serial robots with two flexible links. The desired state trajectories come here from a rigid-body inverse kinematics. In contrast to this, the presented approach calculates the desired state trajectories with a dynamic flexible model inversion. Also, instead of serial robots, parallel robots are considered, which complicates the control design. In this regard, the DAEs are transformed to ODEs via projection, which allows us to use the standard LQR algorithm after a linearization of the minimal form.

Feedback linearization
Since the LQR does not directly reduce the end-effector trajectory tracking errors, the concept of feedback linearization [14] is considered.
Feedback linearization has been applied to flexible link robots before. Nevertheless, due to the typically nonminimum phase system property when tracking the exact end-effector, input-output feedback linearization has been mainly used to track joint trajectories for flexible link robots [2], such as in [20], where end-effector trajectory tracking of a serial robot with flexible and rigid links is considered. Hereby, a collocated output, being the joint angles, is used for the input-output feedback linearization. Also, in [15], end-effector trajectory tracking of a serial robot with two flexible links is performed via computed torque, being a special form of full state feedback linearization, but based only on a rigid-body model. In [4] the theory of a computed torque method for serial flexible link robots is discussed, which only controls the rigid coordinates onto a desired trajectory. An additional stabilization of the elastic oscillations is then superimposed. Still, in [11] a point along a one-link flexible arm is tracked via input-output feedback linearization within simulations. The idea of output redefinition is also employed in [23], where an output close to the exact end-effector is used for trajectory tracking of a serial two-link robot with one flexible link. There the tracking is realized via input-output linearization validated within experiments.
The presented research also experimentally applies the concept of input-output feedback linearization for trajectory tracking of a redefined output close to the exact end-effector position. However, in contrast to the cited literature, we consider flexible link robots that also have a parallel part. Here we propose a systematic approach to realize the input-output feedback linearization via servo constraints.

Actuator cascade control
All control concepts are designed to convert the user-defined desired end-effector trajectories to actuator position, velocity and current, i.e., force or torque, trajectories. These actuator trajectories are then sent to actuator cascade controllers, which feed back the corresponding actuator measurements. This ensures that disturbances such as friction within the actuators can be effectively compensated.

Main contribution
The main contribution of this research is the proposal of an input-output feedback linearization formulation for flexible link parallel robots, i.e., flexible multibody systems with kinematic loops. Here geometric and servo constraints are combined resulting in a compact control law. Furthermore, we present an approach that uses a dynamic flexible model inversion together with an LQR. Here we propose a design that enables the application of the standard LQR algorithm to flexible multibody systems with algebraic constraint equations. For these model-based controllers, we examine different approaches that transform the equations of motion from DAEs to ODEs. We experimentally validate all concepts within end-effector trajectory tracking scenarios with a flexible link parallel robot.

Structure
In Sect. 2, we introduce the system that is used as an application example and for validation purposes throughout the paper. In Sect. 3, we discuss the utilized modeling, inversion, and solution approaches. In Sect. 4, we experimentally apply model-based end-effector trajectory tracking controllers.

Considered system
To clarify and validate the discussed modeling and control techniques, throughout this paper, we use a real flexible link parallel robot as an application example. Its components and notation can be seen in Fig. 1. Here the motor currents are the control inputs u of two ironless linear motors from KML, which actuate the system. Three 2 mm thin spring steel sheets are used to create a short link 1 and a long highly elastic link 2. They are connected by revolute joints and mounted on the motor sliders, which results in a parallel robot. This parallel part is often also denoted as a kinematic loop. The utilized sensors comprise optical encoders to measure the f a = 2 actuator positions s = [s 1 , s 2 ] T , i.e., the slider positions, as well as f κ = 2 strain gauge measurement points on the long link 2 to obtain information on the curvature κ ∈ R fκ . For validation purposes of the end-effector trajectory tracking performance, a camera is mounted over the robot. Two reflective markers are used to reconstruct the xyposition of the exact end-effector r. The advantage of this camera setup is that it enables a direct measurement without relying on a robot model, but the disadvantage is that the image processing needs to be done offline in a post-processing step. As a result, all later presented end-effector measurements are instead obtained with an estimator running in real-time. The sample time of the real-time target is 0.5 ms, which is also the sample time of the estimators and controllers within experiments.

Modeling
Since the application focus of this paper is model-based control for flexible link parallel robots, an accurate and computationally efficient model is required. The main steps used to obtain a flexible multibody model including kinematic loops are outlined in Fig. 2.
This illustration also shows how the subsequent model inversion via servo constraints is related to the later applied end-effector trajectory tracking control approaches. The details of this overview are discussed throughout this paper.

Forward model
In this research, the elements which show significant structural oscillations are modeled as flexible bodies via the linear FEM. This enables a straightforward description of arbitrary geometries, and with a large number of finite elements (FEs), a high model accuracy is ensured. Subsequently, a modal model order reduction (MOR) [26] is applied to reduce the number of elastic degrees of freedom by retaining only the significant eigenmodes. After a conversion to the standard input data (SID) format [34], the flexible parts are combined with the rigid bodies within the framework of the floating frame of reference (FFOR) [32]. Overall, this modeling process yields a computationally efficient flexible multibody system (FMBS), which is especially important for real-time estimation and control purposes.
Nevertheless, within this model occurring kinematic loops are virtually cut open, which leads to a flexible multibody system in serial, i.e., chain structure; see Fig. 2. These loops are closed again via algebraic geometric constraints c = 0, which, together with the system dynamics, describe the forward model in its original loop, i.e., parallel structure.
These steps are applied to the considered exemplary system of Fig. 1. Here the long flexible link 2 consists of three rigid parts, which are connected by two spring steel sheets. Each sheet is described with 100 Timoshenko BEAM188 elements via ANSYS. The floating frame of reference for link 2 is chosen as a chord frame, which is fixed to the joint on the first slider, whereas the x -axis always points to the revolute joint at the point L. The short link 1 is modeled as a rigid body since it exhibits only negligible oscillations. The arising flexible multibody system in chain structure, i.e., the two serial subsystems, is described with the chain coordinates q ∈ R f . The introduced notion of chain coordinates represents a minimal set of coordinates needed to describe a system in chain structure. They consist of the joint angles α 1 and α 2 , the slider positions s 1 and s 2 , and the elastic coordinates q e ∈ R fe of link 2. For the considered system, it turns out that only the first bending eigenmode is significant in normal operation. Therefore the vector of the elastic coordinates reduces to a scalar q e = q e,1 , i.e., f e = 1, which is used to describe the elastic deformation within all application scenarios. This results in f = 5 chain coordinates The constrained flexible multibody system follows by introducing loop-closing constraints, which are formulated as n c geometric constraint equations c(q) = 0. Similar to [31], the general formulation of the forward dynamics is then given by which is a set of DAEs with differentiation index 3. Here M ∈ R f ×f is the symmetric and positive definite mass matrix. The input matrix B ∈ R f ×fa of the input-affine system is multiplied with the control inputs u ∈ R fa . The vector f ∈ R f contains the contributions of the Coriolis, centrifugal, gyroscopic, and elastic inner forces, as well as of the applied forces that are not related to the control inputs u. Moreover, C = ∂c/∂q ∈ R nc×f is the Jacobian matrix of the geometric constraints, and λ ∈ R nc are the corresponding Lagrange multipliers, which represent the reactions enforcing these constraints. For the constraints, the position of the point L, where the cut loop shall be closed, is described with respect to the two serial subsystems via the 2D position vectors r L,1 (q) and r L,2 (q). Thus here the n c = 2 geometric constraints can be written as For further details on the modeling of flexible multibody systems, we refer to the books [30,32].

Model inversion and projections
For the later discussed end-effector trajectory tracking controllers, we need an inverse model. As indicated in Fig. 2, the inverse dynamics is simply obtained by adding n ς algebraic servo constraints ς = 0 to Eqs. (2a), (2b). This gives Here the output r o ∈ R fo is restricted to a desired trajectory r d ∈ R fo through the servo constraints. For the system of Fig. 1, the 2D end-effector position or an approximation of it is the output of interest, i.e., f o = 2, and accordingly n ς = 2 servo constraints are used. Servo constraints [3] are similar to geometric constraints but are not enforced by the Lagrange multipliers λ but by the control inputs u. For the considered system type, where the geometric constraints lead to DAEs with differentiation index 3, the inversion via servo constraints does not change the problem structure. The solution of Eqs. (4a)-(4c) is a major focus of this research. It further represents the considered problem class, input-affine flexible multibody systems with holonomic geometric and servo constraints with a vector relative degree of r deg = {2, . . . , 2}. For details on the relative degree, we refer to [14,28], and for further information on the relation between the relative degree and the differentiation index, we refer to [7].
To simplify the solution of the inverse model of Eqs. (4a)-(4c), the idea is to transform the set of DAEs to ODEs, i.e., to reduce the differentiation index from 3 to 0. This will be realized with a projection that cancels the unknown Lagrange multipliers λ and the inputs u from Eq. (4a). This method is also called the null-space method.
Initially, the constraint equations (4b) and (4c) need to be differentiated twice with respect to time:ċ Even though the terms c and h vanish for the considered constraints, they are kept for generality. It is worth mentioning that and are not mathematical operators but are used for a compact notation throughout this paper. The Jacobian matrices of the constraints and the input matrix can be summarized as Here both matrices T ∈ R f ×(nc+nς ) and G T ∈ R f ×(nc+nς ) have the same size, which is important for the following computations. Thus the number of servo constraints needs to match the number of utilized inputs. If more inputs are available, then a submatrix of the input matrix B should be used within G T . Next, we discuss systematic approaches to arrive at the needed projection and the ODEs.

Manual selection via coordinate partitioning
Firstly, we extend the coordinate partitioning approach [30] to the case with servo constraints. We manually split up the chain coordinates q ∈ R f into independent coordi- Note that this coordinate partitioning is needed for the subsequent steps, but it does not imply that the first f i coordinates within q need to be selected as independent. Instead, the independent coordinates can be chosen freely with a subsequent internal reordering within q. For systems where the underactuation is due to the flexibility, the elastic coordinates are a possible choice of the independent coordinates.
Based on the partitioning of the coordinates q and Eq. (6a), Eqs. (5a) and (5c) can be rewritten as where are Jacobian matrices with γ = [c T , ς T ] T ∈ R nc+nς . Solving Eq. (7) forq d giveṡ Similarly, by Eqs. (5b) and (5d) on acceleration level it follows thaẗ To transform the DAEs (4a)-(4c) to ODEs, a projection matrix J ∈ R f ×f i needs to be found to cancel G, i.e., the Lagrange multipliers λ and the inputs u from Eq. (4a). A possible choice is Here G i ∈ R (nc+nς )×f i and G d ∈ R (nc+nς )×f d can be obtained by a partitioning of G analogously to in Eq. (7). The transposed projection matrix J T is now multiplied from the left to the system dynamics (4a), where also the accelerationsq are replaced by Eq. (10), yielding Then the equation of motion (12), which is an ODE for the independent accelerations, is rearranged toq Now, solving Eq. (9) forq i via a left multiplication with the transposed nonsquare matrix J r yields the independent velocitiesq Plugging Eqs. (13) and (14) back into Eqs. (10) and (9) gives the equations of motion in chain coordinatesq These ODEs in state-space representation for the state x = [q T ,q T ] T are in a form that can be directly used within an integrator. Here the right-hand side is evaluated leading to the left-hand side needed by such an integrator. It is worth noting thatq appears on both sides of Eq. (15a). This, however, does not indicate that it needs to be solved forq, but by inserting and projecting it with J r it is ensured that the obtainedq on the left side complies with the constraints on the velocity level. Consequently, the drift through integration is reduced to being only linear within the constraint equations (4b) and (4c) on the position level. For systems where the underactuation is due to the flexibilities, neglecting these flexibilities yields an equivalent fully actuated rigid system. Then the projected parts of Eqs. (15a), (15b) vanish as no independent coordinates exist, which givesq = θ andq = θ . Thus, the inverse is purely algebraic, and instead of time integration, it can also be solved by finding the roots of Eqs. (4b) and (4c), as well as with Eq. (7), to obtain the state trajectories. This corresponds to inverse kinematics.

Automatic selection via QR decomposition
A manual selection of the independent and dependent coordinates as in the preceding coordinate partitioning approach can cause singularities within d , and also G d can be singular, which will cause the algorithm to fail. Therefore an alternative approach is now presented, which uses, amongst others, an automatic selection of the independent and dependent coordinates. Based on [17,31], this approach relies on the QR decompositions which yield the matrices J r,q and J ,q . They have the same dimension and purpose as J r and J from the coordinate partitioning approach but contain different values. With orthogonal matrices Q ∈ R f ×f and Q G ∈ R f ×f , it follows that Thus, analogously to the coordinate partitioning approach, the columns of J r,q and J ,q span the null spaces of and G, respectively. Based on the QR decomposition, we introduce a new set of independent and dependent coordinates according to with q d,q ∈ R f d and q i,q ∈ R f i of the same dimension as in the coordinate partitioning approach. This relates the independent coordinates to an orthonormal basis of the constraint tangent plane, and the dependent coordinates are related to an orthonormal basis of the row space of the constraint Jacobian matrix . These coordinates are in general not single elements of q, i.e., they have no direct physical meaning anymore. Also, the selection changes automatically in each time step depending on q. Its time derivatives can be written aṡ with q d,q and q d,q to be determined. Inserting Eqs. (16a) and (19a) into the constraint equation (7) on the velocity level gives This yields is a square matrix. Analogously, on acceleration level, it follows that q d,q = −R −T r γ . Based on Eq. (17b), the Lagrange multipliers λ and the control inputs u are canceled by a left multiplication with J T ,q ; compare Eq. (12). With analogous steps to arrive at Eqs. (15a), (15b), we have the following equations of motion in chain coordinates: These are again ODEs in state-space representation with the same structure as for the coordinate partitioning. The only difference is that for the QR decomposition, the product J T r,q J r,q yields the identity matrix and has therefore been neglected within Eq. (21a).

Mixed coordinate partitioning
It is worth noting that it is possible to replace J with J ,q within the coordinate partitioning approach to prevent singularities, which might appear in G d , while still keeping the manually selected independent coordinates. This means that a mix of both approaches from Sects. 3.2.1 and 3.2.2 is possible to benefit from the advantages of both methods. Here we denote this mix as a mixed coordinate partitioning. Based on Eq. (13), the equation of motion in minimal form then follows as

Lagrange multipliers and inputs
Next, we can compute the previously canceled Lagrange multipliers λ and the control inputs u independently of the coordinate partitioning or the QR decomposition. For underactuated systems, the matrix G T cannot be inverted, and Eq. (4a) cannot be directly solved for λ and u. Nevertheless, they can be obtained in a form independent ofq. These accelerations are eliminated by combining the system dynamics (4a) with the second time derivatives of the constraint equations (5b) and (5d). This gives the intermediate result Analogously to [5], solving for λ and u leads to with v =r d here. This is an algebraic equation to be evaluated with the system state q, q coming, e.g., from a time integration of the equations of motion of the inverse model. The matrix = M −1 G T needs to be regular and is here denoted as an extended decoupling matrix. The reason is that it can be regarded as an extension of the usual decoupling matrix¯ = H M −1 B; compare [30]. Note that for the considered exemplary system of Fig. 1, there are constellations where becomes singular. This is the case when α 1 is 90 • , or α 2 is −90 • , or when they are close to such an angle depending on the elastic deformation, since then the actuators can move the end-effector only in the x-direction. However, as this only happens close to the actuator limits, it does not affect typical operations. Thus we further assume that is regular, leading to a differentiation index of 3 and to a vector relative degree of r deg = {2, . . . , 2} as mentioned before.

Solving the inverse model
For underactuated systems such as flexible multibody systems, the inverse model from Eqs. (4a)-(4c) has a dynamics corresponding to the so-called (driven) internal dynamics, which is commonly known from feedback linearization [14]. The transformation of the inverse model from DAEs to ODEs, as discussed in the previous sections, will be used to simplify the solution process of this dynamic inverse. In this regard, two solution approaches are  Fig. 3. One approach uses a formulation in independent coordinates with the independent state x i , and one uses a formulation in chain coordinates with the state x. These state variables are defined as with x d being the dependent state. For the minimal form, i.e., in independent coordinates, we consider only a manual selection of the independent coordinates q i . This ensures that the coordinate selection does not change throughout operation, which simplifies the solution process. The proposed process via the mixed coordinate partitioning approach is shown in Fig. 3(a). Here we use a QR decomposition for the projection to prevent singularities within G d . The superscript "iv" in the function symbol f iv i denotes that an inverse model is used. As depicted in step 1), it is required to internally solve for the dependent coordinates q d . Thus there is no drift of the constraint equations (4b) and (4c) on the position level. For the formulation in chain coordinates, illustrated in Fig. 3(b), we propose to purely use the QR decomposition since it also prevents singularities within the coordinate selection.

Boundary value problem
The stability of the internal dynamics determines how the solution process can be performed. If it is unstable, then the solution process is much more complicated. This is usually the case when tracking the exact end-effector of flexible link robots, i.e., when we use r as output r o , such as for the exemplary system of Fig. 1. Then the classical causal inversion, which originates from [12,13], via forward time integration leads to unbounded states and inputs. This can be prevented with the concept of stable inversion [9], applied in this research to obtain bounded results. Within this framework a two-point BVP needs to be solved instead of an IVP, typically much easier to handle. This means that we cannot predefine the initial state, but for the considered system type, a preactuation is necessary to drive the system to the required state while keeping the output constant. Within the preactuation phase an actuator force or torque, which can be regarded as the output of the inverse model, acts before the end-effector moves, which corresponds to the input of the inverse model. Consequently, the inverse has anticausal characteristics and needs to be calculated offline, since the desired end-effector trajectory needs to be known in advance and cannot be applied online. Also, a postactuation is needed at the trajectory end, where the actuators and the internal dynamics come to rest after the end-effector already stopped moving. This is a causal characteristics besides the explained anticausal behavior, meaning that the inverse system incorporates a combination of both characteristics denoted as noncausal [19].
For the stable inversion approach, it is required that the equilibrium points of the internal dynamics are hyperbolic, i.e., the linearization has no eigenvalues on the imaginary axis. Also, the boundary conditions of the BVP require that the internal dynamics starts on an unstable manifold and ends on a stable manifold. These stable and unstable invariant manifolds are the nonlinear analogs of the corresponding stable and unstable eigenspaces, which are tangent to these manifolds at the considered equilibrium [28]. Since these manifolds are typically difficult to handle, they are now locally approximated by the eigenspaces, which has also been done in [36]. The eigenspaces are obtained by a linearization of the internal dynamics around the corresponding hyperbolic equilibrium points at the trajectory start and end, where the output is held constant via servo constraints. This linearization of the so-called zero dynamics, i.e., the internal dynamics with constant output r o , can be performed, e.g., by applying finite differences to the state-space representationẋ Fig. 3. This allows us to obtain the partial derivative with respect to x i or x giving the corresponding constant system matrix A η of the linearized zero dynamics. It can be consequently written as the autonomous systeṁ Then the eigenvectors and eigenvalues of A η are calculated. Subsequently, the possibly complex diagonal form with the diagonal matrix of the complex eigenvalues is transformed to the real block-diagonal form, e.g., via MATLAB's cdf2rdf command. This yields the real block-diagonal matrix D η , which has the same eigenvalues as A η . The corresponding real transformation matrix T fulfills The η-coordinates of the linearized zero dynamics can then be transformed to new coordinates ϕ = T −1 η, leading toφ Due to the block-diagonal form of D η , the stable and unstable dynamics are decoupled in this equation. With T s being the rows of T −1 related to the stable eigenvalues and T u being the rows related to the unstable eigenvalues, we can construct the boundary conditions. At the trajectory start time t = t 0 the boundary condition keeps the internal dynamics on the approximated unstable manifold, and is used for the approximated stable manifold at the trajectory end time t = t F [28,30]. As a result of the real block-diagonal form, these boundary conditions are expressed with real values.
As mentioned before, the eigenvalues of the linearized zero dynamics at the equilibrium points need to be off the imaginary axis. If we apply the discussed approach to a formulation in chain coordinates, then the rows of T −1 related to the additional 2f d zero eigenvalues, introduced through the constraint equations, need to be removed. To make up for the now missing boundary conditions, e.g., c = 0 and ς = 0 may be used at one boundary andċ = 0 andς = 0 at the other.
Finally, the noncausal solution for the control inputs and state trajectories can be obtained by solving the two-point BVP, e.g., with MATLAB's bvp5c solver. Further details on stable inversion can be found, e.g., in the book [30].
Initial value problem Even if a bounded noncausal solution can be computed, often a stable internal dynamics is desired, which can be computed online by forward integration of an IVP to obtain a bounded causal solution. In the case of sufficient computational power, this enables even a real-time capable model inversion, and for the later applied feedback linearization, a minimum phase system is needed.
Since for flexible link robots, the unstable internal dynamics is usually introduced through the structural flexibility, such as for the now considered system of Fig. 1, the idea is to weight the corresponding deformation. Thus a redefined end-effector output close to the exact end-effector position will be identified, which leads to a small tracking error and a stable internal dynamics. The 2D position of the exact end-effector can be described by Here u GN x and u GN y account for geometric nonlinearity and are functions of the linear deformation φ( 2 )q e of the flexible link in the y -direction, and φ is the corresponding row vector of the shape functions. Also, 2 is the undeformed length of the long link 2, and α 1,r is the rotation of the floating frame of reference. As opposed to this, α 1 describes the complete rotation of the joint on the first slider including the part coming from the elastic rotation.
The insignificant deformations in the x -direction are neglected within r. As mentioned before, only the first bending eigenmode is essential. Therefore the elastic deformation φ( 2 )q e is replaced by φ 1 ( 2 )q e,1 , where φ 1 is the displacement shape function related to the elastic coordinate q e,1 of the first bending eigenmode. If we additionally weight this deformation by the output redefinition weight w according to wφ 1 ( 2 )q e,1 , then we obtain the redefined end-effector output .  Fig. 3(a)) stable inversion of a two-point BVP unstable For the weight w = 1, it corresponds to the exact end-effector position, and for w = 0, it corresponds to the end-effector position of an equivalent rigid system. Consequently, a weight w close to 1 needs to be identified, which results in a stable internal dynamics. Then for such w, tracking of the redefined output allows us to use standard ODE integrators to solve the dynamic inverse model via forward integration.

End-effector trajectory tracking control
In this section, we present the design and experimental application of different end-effector trajectory tracking controllers, which are the main result of this research. All control approaches are based on an inverse flexible multibody model obtained via servo constraints, as discussed in Sect. 3.2. Here the stability of the resulting internal dynamics decides which control approach is used; compare the overview in Fig. 2.

Feedforward-based control
In the following, we apply the model inversion discussed in Sect. 3.2 experimentally as feedforward control for the exemplary system of Fig. 1. We implement it via three different cases, which are summarized in Table 1. Here the rigid case corresponds to inverse kinematics for an equivalent rigid system, where the system state, without elastic terms, is obtained via root finding. The input u is computed in a subsequent step. For the redefined output case, we select a weight w that yields a stable internal dynamics. Then we solve an IVP for a formulation in chain coordinates via QR decomposition. This prevents possible singularities, which is especially important if the model inversion will be performed in real-time.
When tracking the exact end-effector position, a BVP is solved due to the unstable internal dynamics. Here the equations of motion are formulated in independent coordinates, which makes the problem more compact. For the considered system, the only retained elastic coordinate q e,1 is selected as independent. It turns out that a rigid initial guess, i.e., with zero deformation, can be used to solve the BVP for the considered trajectory. Since the BVP needs to be solved offline, possible singularities within the underlying mixed coordinate partitioning approach can be caught. These, however, do not occur here.

Feedforward control only with actuator feedback control
A straightforward implementation of the model-based feedforward control within experiments is to only apply the corresponding actuator-related trajectories. This means that only the feedforward trajectories for the actuator motions s,ṡ and for the input currents u are used and feedback controlled. The resulting control structure is shown in Fig. 4. Here "ff" denotes the corresponding feedforward quantities, and the hat symbol denotes estimated quantities. The necessary time derivatives of the position measurements s can be efficiently obtained by a standard Kalman filter, considering these signals as linear single-degree-of-freedom (DOF) systems where the neglected dynamics is included as a process noise [25,33]. This method is here denoted as single DOF filtering. Then the obtained quantities are sent to an actuator cascade control, which compensates disturbances such as friction to a large extent. It consists of three loops, with a P control for the position loop, a PI control for the velocity loop, and another PI control for the loop of the current i. With these loops, we obtain the voltage u v .
Errors within the elastic deformation are not feedback controlled. Nevertheless, for cases with an accurate model-based feedforward control, these errors will be small. Moreover, not using feedback of the elastic link deformations ensures a very stable control behavior also in the presence of significant disturbances. Now experiments are conducted for the exemplary system of Fig. 1. The results for the three model inversion cases from Table 1 are shown in Fig. 5. Here e rms is the root-meansquare error, and e max is the maximum absolute error of the 2D end-effector position. The bars here and of all following measurements represent the mean values of three different experiments. The smooth desired end-effector trajectory is a line of length 0.2 √ 2 m to be tracked in 0.5 s, which will also be used in all following experiments. This trajectory, i.e., the desired end-effector motion starts after 1 s, which allows a preactuation phase for the BVP. For this BVP, the model inversion control tracks the exact end-effector position and therefore introduces only numerical tracking errors in simulations, in contrast to the other two cases, which track end-effector approximations. Still, Fig. 5 shows that the redefined output case in experiments performs very similarly to the exact end-effector case with very small absolute error differences. On the one hand, this is caused by experimental errors, which introduce additional deviations for all cases. On the other hand, this is due to the high value of w = 0.9625, which gives a redefined end-effector output r re close to the exact end-effector r. This value also leads to a stable internal dynamics.
Note that an unscented Kalman filter (UKF) is used for the end-effector position estimation, which is based on a very similar model as the flexible controllers. The implementation details of the UKF can be found in [25]. Nevertheless, both flexible cases clearly outperform the classical rigid case, since the considered motion excites significant bending. Snapshots of the experiments in Fig. 6 confirm a huge improvement introduced through using a flexible model. Due to the poor performance of the rigid case, it will not be further considered in this paper.

Feedforward control with an LQR
In the previous Sect. 4.1.1, only using measurement feedback for the actuators performs very effectively since negligible disturbances on the link deformations exist. Still, scenarios are possible where significant deviations occur in the elastic deformations. Then the state trajectories obtained by a model inversion can be systematically tracked, for example, using an LQR, which uses state feedback. This enables flexible link robots to also react to deviations within the trajectories of the elastic coordinates to ensure a close end-effector trajectory tracking. The LQR formalism (see, e.g., [18]) is designed for linear forward models in state-space representation. Therefore, in the first step the DAEs (2a), (2b), which describe a general forward model of flexible link parallel robots, will be transformed to state space via projections. Here a representation in independent coordinates, i.e., in minimal form, is necessary as the dependent coordinates introduce zero eigenvalues, which render the exemplary system of Fig. 1 to be not stabilizable. The forward model needed for the LQR simply corresponds to the inverse model without servo constraints, i.e., n ς = 0. Therefore the equations introduced in Sect. 3.2 can be analogously applied when omitting the terms coming from servo constraints. The matrices of dimension f × f i are now identical and are denoted J c . Here only the Jacobian matrices of the geometric constraints are needed, which have the dimensions n c × f i and n c × f d , respectively. With the projection matrix J c the unknown Lagrange multipliers λ are canceled from the system dynamics (2a). With a derivation analogous to that of Eq. (13), the independent accelerations follow as Here the control input u still occurs since a forward model is considered in contrast to Eq. (13). Also, the terms related to servo constraints within θ vanish, which gives With Eq. (35), the function for the equations of motion in independent coordinates can be established analogously to Fig. 3(a) but without servo constraints. Thus the system is now time-invariant and depends on the input u. The involved steps of the resulting functioṅ Fig. 7. Here, "fw" denotes the forward model. In the second step, this function needs to be linearized to enable the application of the LQR formalism. The linearization is done with respect to x i and u to obtain a linear system with independent state variables. The linearization is performed at steady-state points, i.e., the state x i = [q T i ,q T i ] T is chosen as x i,s = [q T i,s , 0 T ] T , and assuming a planar scenario, the input u is chosen as u s = 0. For the considered exemplary system, q i is now of dimension Here the actuator positions and the elastic coordinates, containing only the first bending eigenmode, are chosen as independent, i.e., q i = [q e,1 , s 1 , s 2 ] T . This is different from the previous Sect  linearization point with zero link deformation. Using such steady-state points represents an effective trade-off, since flexible links oscillate around these points, and typically the nonlinearity comes mainly from the underlying rigid body poses. This also ensures that the LQR design is applicable to a variety of scenarios, such as tracking of initially unknown trajectories, i.e., without linearization around the considered trajectory. Analogously to f fw i within Fig. 7, the output y = r o (q) ∈ R fo can also be written as a function y = r o (q i ), which relies on root finding for the dependent coordinates q d = q d (q i ). Thus this output function can also be linearized with respect to the independent state x i , where the derivatives with respect toq i are zero. The linear state-space representation in independent coordinates follows aṡ s representing variations around the considered steady-state working point x i,s . HereÃ i ∈ R 2f i ×2f i is the constant system matrix,B i ∈ R 2f i ×fa is the constant input matrix, andC i ∈ R fo×2f i is the constant output matrix. The linearization to obtain numerical values for these matrices can be done, e.g., via finite differences. For the LQR, the cost function being minimized for infinite final time reads SelectingQ i =C T iC i ∈ R 2f i ×2f i minimizesỹ Tỹ ≥ 0, which guarantees thatQ i is positive semidefinite. Since the output r o does not depend on the velocitiesq i , theQ i -matrix contains nonzero elements only in the upper left block. Additionally using the exact endeffector r as output, i.e., y = r o = r, ensures that its error is kept small, which is crucial for the considered end-effector trajectory tracking. Note that since the LQR is still no output controller, no problems occur when using the exact end-effector. The matrixR i ∈ R fa×fa is chosen as diagonal with positive elements to obtain a positive definite matrix. The control gains K lqr ∈ R fa×2f i for the state feedback control law can then be calculated, e.g., via MATLAB's lqr command, which solves the arising algebraic Riccati equation. For each linearization point, an individual gain matrix K lqr is computed. To account for nonlinearities nonetheless, these gain matrices are interpolated, i.e., a gain scheduling is performed for K lqr . This is possible since the same fixed set of independent coordinates q i is used for all linearization points. For the considered exemplary robot, K lqr only depends on the sum of the two actuator positions. The reason is that linearization points with the same sum lead to an identical kinematic setup only with a possible shift in the x-direction of the complete robot. The currents obtained via the LQR are then added to the currents obtained via the feedforward control for the exact end-effector as described in Table 1. These currents are denoted as u lqr and u ff , respectively; see Fig. 8. Also, in this control structure, we can see that the standard LQR control law from Eq. (40) is slightly adapted by feeding back a trajectory tracking error. This is the error between the independent statex i , estimated by a UKF, and the precalculated independent state x i,ff of the feedforward control. This ensures that the state trajectories x i,ff are closely tracked. Since these trajectories are obtained via a model inversion for the desired end-effector trajectory r d , the LQR can even improve the endeffector output trajectory tracking performance.
Directly applying the input currents u lqr + u ff is usually not meaningful, e.g., because of friction effects within the actuators. Therefore the currents are converted to actuator motions via the projected system dynamics and the estimated state to realize these motions via the actuator cascade control; compare Fig. 8. Here Eq. (35) can be used, which is based on the coordinate partitioning approach. However, since the incorporated projection needs to be applied in real-time, it is advisable to obtain the projection matrix J c,q ∈ R f ×f i based on a QR decomposition. This prevents possible singularities within C d , which can occur through a manual selection of independent coordinates. With the matrices from Eqs. (6a), (6b), which reduce to of dimension n c × f for the forward model, the QR decomposition follows as This also gives Q c ∈ R f ×nc and the square matrix R c ∈ R nc×nc . Also, the unphysical independent accelerationsq i,q need to be transformed back to the accelerationsq in chain coordinates. Since here with θ c,q = −Q c R −T c c . Fromq, resulting from evaluating the right side of Eq. (44), the actuator accelerationss can be selected, which are then integrated by the simple forward Euler method. Finally, the obtained desired actuator motions s d andṡ d can be accurately reproduced via the actuator cascade control. Now experiments are conducted for the exemplary system of Fig. 1, where the LQR is used together with the feedforward control for the desired line trajectory. Here 101 steadystate points are used for the linearization for the LQR design. The linearization is performed in a preprocessing step. Also, both diagonal elements of the diagonal matrixR i are chosen as 10 −4 , which is found via tuning. Note that the SI units forR i and for further control gains are omitted for easier readability.
It turns out that the experimental end-effector trajectory tracking performance with an additional LQR is very similar when only using actuator feedback control such as in the previous Sect. 4.1.1. Nevertheless, the LQR can be highly advantageous for cases with larger disturbances or modeling errors. In the experiments presented in Fig. 9 an unmodeled mass of 0.4 kg is attached to the end-effector, which reduces the first bending eigenfrequency by around 15%. Here the feedforward control, the LQR, and the state estimator, on which the LQR relies, all lack knowledge of the added mass. Therefore, to still evaluate the control performance of the exact end-effector, a second informed UKF is used, which knows the added mass. The results in Fig. 9 are completely based on this informed UKF. They clearly show the improved end-effector trajectory tracking performance with an additional LQR in the presence of significant modeling errors.

Feedback linearization
Up to now the actuator motions and the system state have been controlled in a feedback loop to provide end-effector trajectory tracking. Still, the error of the actual tracking output r o has not been directly feedback controlled yet. Therefore the concept of feedback linearization [14] is now considered to obtain a linear input-output relation that is easy to control. Here an approach is proposed for feedback linearization of flexible multibody systems with Fig. 10 Control structure of the input-output feedback linearization for flexible link robots with geometric and servo constraints geometric and servo constraints. Due to the existence of an internal dynamics when using a flexible model, full state feedback linearization is not possible. Therefore in this section, we focus on input-output feedback linearization using a flexible multibody model.
Tracking of a redefined output, i.e., r o = r re , ensures that a minimum phase system is obtained. For the considered system of Fig. 1, this results in two output components, which will be controlled by two inputs, being the currents of the linear motors. Also, as explained in Sect. 3.2.4, the extended decoupling matrix is regular for the investigated scenarios. This leads to the fulfilled sufficient conditions for input-output feedback linearization with geometric and servo constraints via a static state feedback law: (1) The system is minimum phase, due to tracking of an appropriate redefined output r re .
(2) The numbers of inputs and outputs, introduced via servo constraints, are equal. The implemented control structure is illustrated in Fig. 10 explained in the following. Feedback linearization is usually applied to systems without geometric constraint equations. Still, for the considered system type, geometric constraints due to the parallel part exist, but also servo constraints are used. To apply the concept of input-output feedback linearization, nonetheless, these constraints are incorporated via their second time derivatives. Together with the system dynamics, Eq. (24) has been obtained. This equation already represents the static state feedback law of the input-output feedback linearization to obtain the required control inputs u and the Lagrange multipliers λ. To evaluate Eq. (24), the state x = [q T ,q T ] T is needed, which is estimated via a UKF based on the measured actuator positions s and the link curvatures κ. This means that also the current link deformations are considered in the controller. Note that no transformation to the (Byrnes-Isidori) input-output normal form [6] is needed to obtain the input-output linearizing controller. Similarly, in Sect. 3.2 the internal dynamics is directly obtained via projections (compare Fig. 3) instead of a transformation to such an input-output normal form. The internal dynamics is also not needed to calculate the input-output linearizing controller, but it is required for stability analysis. Here the internal dynamics has been used to find a redefined output r re that renders this internal dynamics stable.
Applying the control input u, obtained via Eq. (24), then theoretically cancels all nonlinear terms. Due to the incorporated servo constraints, this yieldsr re = v, i.e., a linear differential input-output relation between the new control input v and the output r re . Selecting v =r d givesr re =r d in an ideal scenario. Due to the input-output feedback linearization, occurring tracking errors within realistic scenarios can thus be compensated by a simple linear controller of the form for the considered system. This leads to a stable output error dynamics for P , D > 0. Also, via Eq. (45) the desired end-effector trajectory is now applied in real-time. Regarding the Lagrange multipliers λ as additional control inputs, we can analogously stabilize the geometric constraints c = 0 [10]. This is, however, not necessary when using a real system where these constraints cannot drift. As a next step, the control inputs u, i.e., the motor currents, are again transformed to actuator, i.e., slider accelerationss, to guarantee, amongst others, an effective friction compensation via the utilized actuator cascade control. In contrast to the LQR of Sect. 4.1.2, where the unknown Lagrange multipliers are canceled by a projection leading to Eqs. (35) and (44), the λ-vector is automatically calculated via Eq. (24) of the feedback linearization. This vector can then be simply plugged into Eq. (4a) together with the available control input u and the state x, which gives Hereq contains the needed actuator accelerationss. Evaluating Eq. (46) is computationally very efficient since the inverse of M is already known from Eq. (24) and the forward Euler method is accurate enough for the subsequent time integration. The resulting desired actuator motions s d andṡ d are then closely reproduced via the actuator cascade control. Alternatively, if the currents u are not needed, e.g., to apply the control only on position and velocity level, then the actuator accelerations can be obtained directly via Eq. (21b). Then a replacement ofr d by v within θ q is necessary.
It is worth noting that a major difference from the previously discussed flexible feedforward controllers, which are dynamic controllers, is that the input-output feedback linearization only uses an algebraic, i.e., static control law, which is computationally much more efficient. Nevertheless, feedback linearization requires a state estimator, which can be computationally expensive. This is especially the case for the utilized UKF, which internally integrates the system dynamics. Now the input-output feedback linearization is applied experimentally for the desired line trajectory using the flexible link parallel robot of Fig. 1. Within the resulting plots of Fig. 11, we can observe that a larger weight w, i.e., choosing a redefined output r re closer to the exact end-effector position output r, yields a more effective trajectory tracking behavior of r. Also, the redefined output, which is in fact tracked by the controller, is kept very close to the desired end-effector trajectory that validates the control approach.
The corresponding errors are shown in Fig. 12 for different weights w and PD control gains. Here the small errors for the redefined output r re are comparable to the flexible feedforward results from Fig. 5. Therefore only a minor improvement is visible when using the nonzero PD gains. Nevertheless, the PD gains help to reduce the errors in all cases and also have a positive effect on the exact end-effector r. Obviously, the main error within the exact end-effector comes from output redefinition needed to obtain a minimum phase system. This error can only be safely decreased if, e.g., disturbances are reduced to allow a higher value than the robust value of w = 0.75. The experimental setup under consideration turned out to become unstable roughly at w = 0.85. This is still somewhat away from the weight of w = 0.9625, which can be used within a flexible feedforward control in real-time. Thus  the redefinition error for the feedback linearization is significantly higher than for the feedforward control. As a result, the feedback linearization based on output redefinition is mainly useful for scenarios where using a state feedback is so advantageous that it can outweigh this error discrepancy. Nevertheless, for a scenario in which the desired end-effector trajectory is known in advance, such as for repetitive tasks, the error introduced through output redefinition can be precompensated. With the offline calculated states from a model inversion for the exact end-effector case formulated as BVP, the corresponding redefined output motion can be evaluated. Applying this motion, denoted as adapted "line" trajectory, a new desired trajectory will ensure that a close tracking of r re of this adapted line results in a close tracking of r of the original line. The highly reduced tracking errors of the original line for the exact end-effector via this trajectory adaption method can be seen in Fig. 13. The performance is also close to the flexible feedforward control experiments from Fig. 5 with the possibility to compensate the disturbances on the flexible link. Furthermore, the PD control gains now improve the exact end-effector trajectory tracking error by a significant percentage. Interestingly, for such an adapted trajectory, the case w = 0.3 performs better than the case w = 0.75. However, as this cannot be confirmed by simulations, it is most likely an experimental effect occurring through the different weighting of the elastic coordinate. Finally, the additional end-effector mass from Sect. 4.1.2 is used to face the input-output feedback linearization controller with significant modeling errors. Both the controller and the state estimator, which it relies on, do not know this mass of 0.4 kg. Therefore, for quantification purposes of the actual control performance, the results in Figs. 14 and 15 are estimates of a second informed UKF, which knows the added mass. Also, the output redefinition is not compensated, i.e., no adapted line trajectory is used. The significant end-effector trajectory tracking error caused by the unmodeled mass for the feedforward control is taken from Fig. 9. The input-output feedback linearization without PD control, i.e., without feedback of the output error, leads to a reduced maximum error but to a larger root-mean-square error. In contrast, the input-output feedback linearization with active PD control improves both errors. Especially the reduction of the maximum error compared to the feedforward control is significant, caused by the compensation of the overshoot at the trajectory end; see Fig. 14. Also, the errors with active PD control are comparable to the results with an additional LQR shown in Fig. 9.

Conclusions
This research presents an overview of computationally efficient modeling and end-effector trajectory tracking control approaches for flexible multibody systems with geometric constraints, e.g., due to kinematic loops, and servo constraints. Initially, the flexible multibody model is inverted with the concept of servo constraints. To simplify the solution process, the resulting set of differential-algebraic equations is transformed into ordinary differential equations by projections. Typically, for flexible link robots, such as the considered exemplary system, an unstable internal dynamics occurs when tracking the exact end-effector. Then a two-point boundary value problem is solved offline within the framework of stable inversion. Initially, the resulting feedforward control based on flexible model inversion is applied only with actuator feedback. Here tracking of a redefined end-effector output, which leads to a stable internal dynamics that can be solved as initial value problem, performs similarly to tracking of the exact end-effector position within experiments. Also, both flexible model inversions show a significantly improved end-effector trajectory tracking performance compared to feedforward control based on classical rigid body modeling.
Subsequently, to also compensate errors within the elastic deformation, not only the actuator feedback is used, but also the state feedback. Here the well-known linear-quadratic regulator approach is adapted to trajectory tracking control for flexible multibody systems with geometric constraints. This leads to a superior end-effector trajectory tracking performance within experiments where an unmodeled mass is attached to the end-effector of the considered flexible link parallel robot. Finally, for a redefined output rendering the system minimum phase, the well-established input-output feedback linearization is applied. Here a systematic approach is proposed for flexible multibody systems with geometric and servo constraints. Experimental results validate the control design. Assuming a repetitive task allows us to adapt the desired trajectory offline to precompensate the error introduced by output redefinition. This further reduces the experimental end-effector trajectory tracking errors. The input-output feedback linearization also proves to be valuable for the case of an unmodeled mass at the end-effector.

Conflict of Interest
The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.