*Article* **Modeling Parallel Robot Kinematics for 3T2R and 3T3R Tasks Using Reciprocal Sets of Euler Angles**

### **Moritz Schappler \*, Svenja Tappe and Tobias Ortmaier**

Institute of Mechatronic Systems, Leibniz University Hannover, 30167 Hannover, Germany **\*** Correspondence: moritz.schappler@imes.uni-hannover.de

Received: 27 June 2019; Accepted: 1 August 2019; Published: 6 August 2019

**Abstract:** Industrial manipulators and parallel robots are often used for tasks, such as drilling or milling, that require three translational, but only two rotational degrees of freedom ("3T2R"). While kinematic models for specific mechanisms for these tasks exist, a general kinematic model for parallel robots is still missing. This paper presents the definition of the rotational component of kinematic constraints equations for parallel robots based on two reciprocal sets of Euler angles for the end-effector orientation and the orientation residual. The method allows completely removing the redundant coordinate in 3T2R tasks and to solve the inverse kinematics for general serial and parallel robots with the gradientdescent algorithm. The functional redundancy of robots with full mobility is exploited using nullspace projection.

**Keywords:** parallel robot; five-DoF task; 3T2R task; functional redundancy; task redundancy; redundancy resolution; reciprocal Euler angles; inverse kinematics

### **1. Introduction**

Industrial tasks like welding, gluing, milling or drilling represent a major part of the applications of industrial robots, which generally have full mobility, i.e., the operational space of their end-effector has three translational and three rotational ("3T3R") degrees of freedom ("DoF"). Parallel robots like the Stewartplatform have especially been proposed for milling tasks regarding their high structural stiffness. The task space of the named applications can be defined by three translational DoF and only two rotations due to a symmetry around the tool axis ("3T2R"). This results in a functional or task redundancy, which is not exploited to full extend yet for *parallel* robots.

### *1.1. Inverse Kinematics and Resolution of Task Redundancy for Serial-Link Robots*

Various general gradient-based methods exist to solve the inverse kinematics for *serial* robots; either by augmenting the joint space [1] or by reducing the task space [2–5]. The different approaches each define a residual vector and a gradient matrix considering the properties of 3T2R tasks, e.g., by adding a virtual joint axis [1], orthogonal decomposition of the task space [2], rotation of the residual into a task frame and removing the corresponding component [3], defining the tool axis by two points for constructing a nullspace [4] or by defining the absolute orientation and the orientation residual with two reciprocal sets of Euler angles [5]. The gradient matrices corresponding to the different residuals are used for an iterative Newton-Raphson algorithm [6,7] by exploiting the functional redundancy with a null space projection of additional performance criteria [6]. Without the definition of a proper gradient, a global optimization has to be performed outside of the inverse kinematics algorithm [8,9].

### *1.2. Overview of Parallel Robots Structures for 3T2R Tasks*

Parallel robots in 3T2R tasks can be ordered in classes according to their kinematic structure into

I mechanisms with full platform mobility (3T3R) that are redundantly controlled to five DoF,

P S


The classes I and II were introduced in [10], where class IV is analyzed regarding leg symmetry and singularities. Class III is mainly influenced by the systematic synthesis of [11] and several existing prototypes and is demarcated against class II by the absence of the passive constraining leg. Class IV can be seen as a subclass of III, but is differentiated in this paper due to its characteristics. Other classifications are provided e.g., by [12], where IV and II are termed "families".

Examples for the first class are hexapods (6UPS) [13] or the Eclipse [14] machine tool (2PPRS-PPRS). (The joint structure of the parallel robots is denoted by the number of the legs and the order of universal ("U"), prismatic ("P"), spherical ("S"), helical ("H") and revolute ("R") joints in the leg chains. Different actuated legs are connected by hyphens ("-"), passive constraint legs are connected by a slash ("/"). Actuated joints are underlined.) Any other parallel robot with full mobility (see e.g., [11,15,16]) may be used as well.

The second class allows for more variety, since the six-DoF mechanism and the five-DoF constraining leg can manifest in different kinematic structures: The UPS structure is used for the six-DoF part of the mechanism by [17] with a focus on kinematic analysis of 5UPS/US, by [18] with a focus on kinetostatic modeling at the example of 5UPS/RUU (see Figure 1a), by [19] with focus on trajectory control of 5UPS/PRPU and by [20] for pose measurement with the passive leg of 5UPS/PRPU. Other possible general base structures are RUS at the 5RUS/US example in [17], PUS, which has been investigated for the control of a redundantly actuated 6PUS/UPU regarding the force control of the redundant leg with inverse dynamics compensation [21] or force optimization [22].

The most-straightforward member of the third class is the 4UPS-UPU of Figure 1b, which is investigated in [23] for a simulation and feasibility study together with a survey on possible architectures for a technical realization of this class. Other possible structures are the 4URS-URU, which is analyzed kinematically in [24] and the 4PSU-PU\*U, which has a special parallelogram structure in one leg (termed "U\*") and is presented in [25]. A subclass of III consists of mechanisms [26–28], where the last joint axis of the legs is coaxial with the tool axis and is constructed as rotating ring. It contains the Metrom machine tool (4SPRR-SPR), depicted in Figure 1c, which is analyzed regarding inverse and forward kinematics in [26] or its variants, the redundant 4SPRR-PSPR from [28] or the hybrid 4URHU-URHR with an additional linear actuator at the platform [27]. A structural synthesis based on linear transformations and evolutionary morphology [11] led e.g., to the Isoglide5 mechanisms (3PRRRRR-2PRRRR), which are analyzed and optimized regarding the isotropy of the Jacobian in [29].

The simplest member of class IV, the 5UPU is shown in [30] with the help of screw theory to only have local mobility and no global mobility, since the twist systems of the leg chains have no intersection and the resulting twist system of the platform is empty. Members of class IV have been found by systematic structural synthesis with screw theory, which has been performed in [31] for symmetric 3T2R mechanisms. The resulting 5RPUR of Figure 1d and 5PRUR are analyzed in [10,32].

In a practical application with competing requirements on workspace, stiffness, costs and precision, each of the existing systems has its legitimization. Nevertheless, each of the classes has inherent disadvantages: In tasks like milling with high process forces and requirements on stiffness and precision, robots from class II with one constraining leg have the drawback that the passive leg takes the complete reaction wrench in the blocked rotational degree of freedom, which strongly affects the mechanisms stiffness in this direction [18]. The same is argued by [27] at the example of the Metrom machine tool, but can be extended to all the members of class III, where most leg chains have six DoF and usually only one leg chain has five DoF. This leg chain also has to take the reaction moments in the blocked DoF which affects the overall stiffness. Therefore, members of the classes I and IV can be expected to reach a higher stiffness. Mechanisms of the class IV may further suffer from an increased sensitivity of manufacturing tolerances, which may cause a high pretension of the bearings or even reduce the DoF, since the five DoF of all leg chains have to coincide exactly to allow the platform to also have five DoF. Additionally, only members of class I provide redundancy which allows using the additional DoF for performance optimizations, e.g., to avoid singularities and to compensate the smaller workspace caused by the sixth leg.

Therefore, the remainder of the paper focuses on mechanisms of the first class to allow an optimization of their performance criteria using the degree of task redundancy.

### *1.3. Inverse Kinematics of Parallel Robots for 3T2R Tasks*

The parallel robots with five DoF presented above have a kinematic structure which allows for an analytic model of the inverse kinematics. All references define the end-effector orientation with two consecutive elementary rotations, i.e., define two Euler angles to represent the tool axis orientation in minimal coordinates [10,13,19–22,26,28,32,33], which is called "partial pose" in [13]. The inverse kinematics problem (IKP) is first solved for the first chain, which is called "leading chain" in this paper. Due to the geometry of the leading chain, this solution can be found algebraically. Then the IKP is solved for the other "following" chains with the given orientation from the leading chain and standard methods. For robots of class II, the constraining leg is selected as the leading leg chain and for class III the 3T2R leg is selected.

To the best knowledge of the authors, only one reference [13] for the IKP of functionally redundant parallel robots of the class I is known. The reason presumably is that a solution of the 3T3R IKP for these robots is possible with standard methods, as used in [14] for the 3T3R Eclipse. It is always possible to transfer the 3T2R IKP into 3T3R by adding an arbitrary value for the desired rotation around the tool axis. An optimization of additional performance criteria is possible by varying the redundant rotation angle [8,9]. This approach was chosen in [13] by first defining the IKP with the redundant rotation as a parameter and then performing an optimization of this parameter using analytical computation of the dexterity and interval analysis to ensure a minimum determinant of the inverse Jacobian. The drawback of this method is the need for a cascaded optimization which is more complex than the gradient-based approach presented in this paper.

### *1.4. Motivation and Summary of the State of the Art*

The overview over the literature shows that no general, machine-independent method for the resolution of functional redundancy for 3T3R PKM in 3T2R tasks exists. The works either focus on a general structural synthesis of machines, e.g., via screw theory [31] or linear transformations [29] or the description and improvement of specific, manually selected machines. To choose the best machine for given requirements, a structural synthesis is only the first step. Additionally, a dimensional synthesis should be performed for all possible structures to select the most suitable mechanism. This combined structural and dimensional synthesis [34] is sketched in Figure 2. To be able to perform the dimensional synthesis for all structures, the inverse kinematics has to be implemented in a general form to calculate the performance criteria over a given trajectory for further optimization of the structures and their comparison. For the generation of task redundant parallel robots, the inverse kinematics has to include an optimization of the performance criteria and the restrictions such as joint limits to ensure a comparability of the results.

**Figure 2.** Overview of the procedure for combined structural and dimensional synthesis.

To address this, the contributions of this paper are


The remainder of the paper is structured as follows: The description of the inverse kinematics problem and prior definitions are given in Section 2 and the concept of reciprocal Euler angles from [5] is adapted in Section 3 for parallel robots. This leads to the full kinematic constraints of parallel robots in 3T2R tasks, introduced in Section 4 and applied to the differential kinematics in Section 5. The theoretical analysis is followed by examples and simulations in Section 6. The appendix contains proofs and additional details on the mathematical formulation.

### **2. Inverse Kinematics Problem for Parallel Robots**

Before addressing the specific model for parallel robots in 3T2R tasks in the next section, the standard kinematics model of parallel kinematic machines ("PKM") is repeated in the following, corresponding to the state of the art [11,15,35] and serving as a reference to highlight its shortcomings for 3T2R tasks. The regarded parallel robot consists of *m* legs, which each have the joint coordinates *q<sup>i</sup>* . All joints are considered to be single-DoF and additionally to the active joints *qi*,a explicitly all passive joints at the base and at the platform *qi*,p are included in the coordinates *q<sup>i</sup>* of leg *i*. The coordinates

$$\mathbf{x} = \begin{bmatrix} \mathbf{x}\_{\mathbf{t}}^{\mathrm{T}} & \mathbf{x}\_{\mathbf{t}}^{\mathrm{T}} \end{bmatrix}^{\mathrm{T}} \in \mathbb{R}^{6} \tag{1}$$

of the end-effector platform describe the position and orientation of the end-effector frame F*<sup>D</sup>* with respect to the base frame F0. In the equations, this is marked with left subscript "(0)" for vectors and left superscript "0" for rotation matrices. The platform-related end-effector frame is the desired frame in the inverse kinematics problem and is therefore abbreviated with "D". The position

$$\mathbf{x}\_{\mathbf{t}} = {}\_{(0)}\mathbf{r}\_{D} \in \mathbb{R}^{3} \tag{2}$$

is defined as the origin of the platform frame and the rotation matrix

$$\mathbf{u}^{0}\mathbf{R}\_{D}(\mathbf{x}\_{\mathbf{r}}) = \begin{bmatrix} \mathbf{u}\_{D} & \mathbf{o}\_{D} & \mathbf{a}\_{D} \end{bmatrix} \in \text{SO}(3) \tag{3}$$

of the platform frame is expressed with Euler angles

$$\mathbf{x\_r} = \begin{bmatrix} \beta\_1 & \beta\_2 & \beta\_3 \end{bmatrix}^T =: \mathcal{J} \in \mathbb{R}^3 \tag{4}$$

as a minimal representation of the orientation coordinates. The symbol "*β*" will be used to denote orientations relative to the base frame throughout this paper. The *XYZ* convention

$$\mathcal{R}(\boldsymbol{\beta}) = \mathcal{R}\_x(\beta\_1)\mathcal{R}\_y(\beta\_2)\mathcal{R}\_z(\beta\_3) \in \text{SO}(3) \tag{5}$$

is used for the Euler angles without loss of generality. The relation between joint coordinates *q* and platform coordinates *x* is established with the kinematic constraint equations, for which most commonly the vector loop

$$\Phi\_{\mathbf{t},i}(q\_i, \mathbf{x}) = -\,\_{(0)}r\_{A\_iB\_i}(\mathbf{x}) + \,\_{(0)}r\_{A\_iB\_i}(q\_i) \tag{6}$$

between the position of the platform coupling point *Bi* relative to the base coupling point *Ai* is used for each leg chain *i* [15]. The second term (0)*rAiBi* (*qi* ) corresponds to the forward kinematics of the serial leg chain. The vector

$$
\hat{\mathbf{r}}\_{(0)} \mathbf{r}\_{A\_i \mathbf{B}\_i}(\mathbf{x}) = -\mathbf{r}\_{(0)} \mathbf{r}\_{A\_i} + \mathbf{x}\_{\mathbf{t}} + \,^0 \mathbf{R}\_D(\mathbf{x}\_{\mathbf{f}})\_{(D)} \mathbf{r}\_{B\_i} \tag{7}
$$

includes the term <sup>0</sup>*RD*(*x*r) that depends on the full orientation *x*<sup>r</sup> of the end-effector. For the bigger part of existing parallel robots, the passive joint coordinates can be eliminated analytically from Equation (6), e.g., by using the Euclidean distance for UPS or RPR leg chains or via trigonometry for RRRchains. This is termed "minimal kinematics set" in [15] and leads to the scalar constraint equation

$$
\Phi\_i = \Phi\_i(q\_{i, \mathbf{a}'} \mathbf{x}) \tag{8}
$$

for each leg *i*, which can be assembled to the vector of constraint equations

$$\Phi(q\_{\mathbf{u}'}\mathbf{x}) = \begin{bmatrix} \Phi\_1 & \Phi\_2 & \cdots & \Phi\_m \end{bmatrix}^\mathrm{T} \tag{9}$$

for all *m* legs of the PKM. The differential kinematics of the PKM is calculated with the time derivative

$$\frac{d}{dt}\Phi(q\_{\rm u},\mathbf{x}) = \Phi\_{\rm \partial q\_{\rm u}}\dot{q}\_{\rm u} + \Phi\_{\rm \partial x}\dot{\mathbf{x}} = \mathbf{0} \tag{10}$$

where the passive joint coordinates *q*<sup>p</sup> do not occur, since they have been eliminated in a previous step. Following [11] to avoid the term "Jacobian" for the elements of (10), the inverse-kinematics matrix

$$
\Phi\_{\partial q\_{\mathbf{a}}} = \frac{\partial \Phi}{\partial q\_{\mathbf{a}}} = \begin{bmatrix}
\Phi\_{1, \partial q\_{1, \mathbf{a}}} & 0 & 0 & 0 \\
0 & \Phi\_{2, \partial q\_{2, \mathbf{a}}} & \ddots & 0 \\
0 & \ddots & \ddots & 0 \\
0 & 0 & 0 & \Phi\_{m, \partial q\_{m, \mathbf{a}}}
\end{bmatrix} \tag{11}
$$

of this model has diagonal form and the direct-kinematics matrix

$$\boldsymbol{\Phi}\_{\partial \mathbf{x}} = \frac{\partial \boldsymbol{\Phi}}{\partial \mathbf{x}} = \begin{bmatrix} (\partial \boldsymbol{\Phi}\_1 / \partial \mathbf{x})^\mathsf{T} & (\partial \boldsymbol{\Phi}\_2 / \partial \mathbf{x})^\mathsf{T} & \cdots & (\partial \boldsymbol{\Phi}\_m / \partial \mathbf{x})^\mathsf{T} \end{bmatrix}^\mathsf{T} \tag{12}$$

is fully populated. This definition of the constraints has the following drawbacks:


An alternative kinematic model to encounter the combination of these points is presented in the next section, where the concept of reciprocal sets of Euler angles for the inverse kinematics problem for serial-link robots [5] is transferred to the leading leg of parallel robots.

### **3. Reciprocal Sets of Euler Angles for the Kinematics of a Serial Leg Chain**

To take the rotational symmetry around the tool axis in 3T2R tasks into account, a new set of task space coordinates

$$\boldsymbol{\eta} = \begin{bmatrix} \boldsymbol{\eta}\_{\mathrm{t}}^{\mathrm{T}} & \boldsymbol{\eta}\_{\mathrm{r}}^{\mathrm{T}} \end{bmatrix}^{\mathrm{T}} \in \mathbb{R}^{5} \tag{13}$$

has to be defined. The translational part

$$
\boldsymbol{\eta}\_t = \mathbf{x}\_t = {}\_{(0)}\boldsymbol{r}\_D \in \mathbb{R}^3 \tag{14}
$$

remains unchanged relative to the operational space coordinates *x*. The rotational part

$$\eta\_{\mathbf{r}} = \begin{bmatrix} \beta\_1 & \beta\_2 \end{bmatrix}^{\mathrm{T}} = \underbrace{\begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \end{bmatrix}}\_{=P\_{\eta\_{\mathbf{r}}}} \mathbf{x}\_{\mathbf{r}} \in \mathbb{R}^2 \tag{15}$$

only contains the first two rotational coordinates of *x*. The last operational space coordinate *β*3, the rotation around the *z*-axis *a<sup>D</sup>* of F*D*, is excluded from the task space by the selection matrix *Pη*<sup>r</sup> . To be able to set the rotational DoF around the tool axis in 3T2R tasks arbitrarily and use gradient-based inverse kinematics, *β*<sup>3</sup> has to be eliminated completely from the kinematics Equations (6). To simplify the following elaborations, the platform frame F*<sup>D</sup>* is still identified as the desired frame of the inverse kinematics problem and the end-effector frame that results from the joint angles of leg *i* is now termed F*E*, which corresponds to the forward kinematics of the leg chain. For a formulation without the tool axis rotation, a different constraint definition

$$\Phi\_{\rm t,i}(q\_i, \mathbf{x}) = -\,\_{(0)}r\_D + \,\_{(0)}r\_E(q\_i) = -\mathbf{x}\_{\rm t} + \,\_{(0)}r\_E(q\_i) \in \mathbb{R}^3 \tag{16}$$

containing the vector loop from the robot base frame F<sup>0</sup> to the platform F*<sup>D</sup>* and the leg chain end-effector F*<sup>E</sup>* can be used, where in contrast to (6) only the translational part *x*<sup>t</sup> of the end-effector coordinates appears and not the rotational part *x*r. The vector loop is depicted in Figure 3 for a planar robot with opened (Figure 3a,b) and closed loops (Figure 3c). The triangle represents the end-effector platform and only one leg chain is drawn in the figure.

**Figure 3.** Different cases for the kinematic constraints of the leading chain for the 3RRR example: (**a**) no constraints complied; (**b**) position and tool axis rotation complied; (**c**) all constraints complied.

As a drawback of (16), all joint angles *q<sup>i</sup>* of the leg *i* and not only the coordinates of the first joints counted from the base are now included in the vector

$$
\sigma\_{(0)} r\_E(q\_i) = {}\_{(0)}r\_{A\_i} + {}\_{(0)}r\_{A\_iB\_i}(q\_i) + {}^0\mathcal{R}\_{B\_i}(q\_i)\_{(B\_i)}r\_E.\tag{17}
$$

This implies that the platform is now part of the last link of the considered leg chain, as sketched in Figure 3 by the dashed triangle. To account for the increased number of included joints in (17), the full kinematic constraints

$$\boldsymbol{\Phi}\_{i} = \begin{bmatrix} \boldsymbol{\Phi}\_{\mathrm{t},i}^{\mathrm{T}} & \boldsymbol{\Phi}\_{\mathrm{r},i}^{\mathrm{T}} \end{bmatrix}^{\mathrm{T}} \in \mathbb{R}^{6},\tag{18}$$

have to be considered, including the rotational part

$$\boldsymbol{\Phi}\_{\mathbf{r},i}(\boldsymbol{q}\_{i},\mathbf{x}) = \begin{bmatrix} \boldsymbol{a}\_{1} & \boldsymbol{a}\_{2} & \boldsymbol{a}\_{3} \end{bmatrix}^{\mathrm{T}} = \mathfrak{a}\left(^{\mathrm{D}}\mathcal{R}\_{\mathrm{E}}(\mathbf{x}\_{\mathrm{r}},\boldsymbol{q}\_{i})\right) = \mathfrak{a}\left(^{\mathrm{0}}\mathcal{R}\_{\mathrm{D}}^{\mathrm{T}}(\mathbf{x}\_{\mathrm{r}})^{\mathrm{0}}\mathcal{R}\_{\mathrm{E}}(\boldsymbol{q}\_{i})\right),\tag{19}$$

which is needed to generate enough equations for an invertible matrix in the differential equations. The constraints again contain the deviation between the desired end-effector frame F*<sup>D</sup>* expressed with *x* and the actual robots end-effector frame F*<sup>E</sup>* expressed with *q*. Figure 3b,c show cases, where the translational constraints are met, but the rotational constraints have different values. For 3T3R tasks, only Figure 3c represents a valid solution of the inverse kinematics. For 3T2R tasks, Figure 3b,c represent valid solutions.

The goal of eliminating the tool rotation *β*<sup>3</sup> from the equations is not achieved yet, since all three components of the platform orientation *x*<sup>r</sup> affect the rotation matrix <sup>0</sup>*RD*. This can be addressed by the selection of the Euler angles: Similar to the definition of the rotational operational space coordinates *x*r in (4), the constraints *Φ*r,*<sup>i</sup>* are also expressed with a set of Euler angles *α*. In the following, "*α*" will always refer to the rotation error/residual and "*β*" to an orientation relative to the base frame. The Euler angle convention of *α* can be chosen independently of the choice for the orientation representation in *β*. The intuitive approach of choosing

$$\mathcal{R}(\mathfrak{a}^\*) := \mathcal{R}\_{\mathfrak{X}}(\mathfrak{a}\_1^\*) \mathcal{R}\_{\mathfrak{Y}}(\mathfrak{a}\_2^\*) \mathcal{R}\_{\mathfrak{z}}(\mathfrak{a}\_3^\*) \in \text{SO}(\mathfrak{z}) \tag{20}$$

the same way as *β* leads to a set of transformations depicted in Figure 4a, where the intermediate steps of the single elementary rotations are omitted since they have no technical meaning. The superscript asterisk in (20) demarcates this specific example and the following elaborations on the calculation of *α*.

To be able to remove the redundant coordinate *β*<sup>3</sup> from the rotational constraints of (19), it is necessary to change the expression of the orientation error *α* to be reciprocal to the expression of the absolute orientation *β*. By using the *ZYX* Euler angles with

$$\mathcal{R}(\mathfrak{a}) := \mathcal{R}\_z(a\_1)\mathcal{R}\_y(a\_2)\mathcal{R}\_x(a\_3) \in \text{SO}(\mathfrak{z}),\tag{21}$$

only the error component *α*<sup>1</sup> is affected by rotations around the tool axis, which is the *z*-axis of the intermediate frames F*A*1, F*A*<sup>2</sup> and the platform frame F*<sup>D</sup>* in Figure 4b, where the frame rotations with the reciprocal set of Euler angles are sketched.

**Figure 4.** Overview of the different frames (**a**) for six-DoF tasks with standard Euler angle convention and (**b**) for five-DoF tasks with reciprocal Euler angle convention; taken from [5].

The mathematical proof is given in Appendix A.1 and in [5]. The new, reduced rotational part of the kinematic constraints

$$\Psi\_{\mathbf{r},i}(q\_{i\prime},\boldsymbol{\eta}) = \begin{bmatrix} \boldsymbol{a}\_2 & \boldsymbol{a}\_3 \end{bmatrix}^{\mathrm{T}} = \overbrace{\begin{bmatrix} 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix}}^{-P\mathbf{v}\_{\mathbf{r}}} \Phi\_{\mathbf{r},i}(q\_{i\prime},\mathbf{x}) \in \mathbb{R}^2 \tag{22}$$

does not contain the tool rotation any more. The full kinematic constraints

$$\mathbf{Y}\_{i} = \begin{bmatrix} \boldsymbol{\Phi}\_{\mathbf{t},i}^{\mathrm{T}} & \mathbf{Y}\_{\mathbf{t},i}^{\mathrm{T}} \end{bmatrix}^{\mathrm{T}} \in \mathbb{R}^{5} \tag{23}$$

for the reduced coordinates can be used for the inverse kinematics of the leg chain *i* in 3T2R tasks, where *Ψ<sup>i</sup>* = **0** leads to a valid position and orientation of the tool axis. The condition *Φ<sup>i</sup>* = **0** leads to a valid configuration of leg *i* of the parallel robot in 3T3R tasks. In the following, "*Φ*" is always used for 3T3R kinematic descriptions and "*Ψ*" for 3T2R. Note that by omitting the corresponding lines in the operational space coordinates *x* and the constraint equations *Φ*, it is also possible to use the 3T3R approach for systems with reduced mobility of 2T1R, 3T0R and 3T1R platform DoF.

### **4. Full Kinematic Constraints for Parallel Robots Using Reciprocal Sets of Euler Angles**

The definition of the full kinematic constraints (16,19,18) of a single leg chain of the parallel robot from the previous chapter can be used to write the kinematic constraints in a general form. The full kinematic constraint equations can only be defined for 3T3R tasks without further adjustments as

$$\boldsymbol{\Phi} = \begin{bmatrix} \boldsymbol{\Phi}\_1^T & \boldsymbol{\Phi}\_2^T & \cdots & \boldsymbol{\Phi}\_m^T \end{bmatrix}^T. \tag{24}$$

The constraints *Ψ<sup>i</sup>* from (23) for the reduced coordinates *η* can only be defined for one leg chain: Figure 5a show an open-loop second leg chain for a given first leg chain from Figure 3.

By also closing the *3T2R* kinematic constraints *Ψ*<sup>2</sup> for the second loop, as depicted in Figure 5b, the tool axis stays arbitrary and the platform pose demanded from the two legs would be different and therefore would not be a valid solution for the complete mechanism, i.e., *Φ*<sup>2</sup> = **0**. Only if the second leg fulfills the *3T3R* kinematic constraints for all platform coordinates, as shown in Figure 5c, a valid configuration of the mechanism emerges. This approach has already been used for many specific robots systems, as introduced in Section 1. As a generalization, the first leg of the parallel robot is now termed the "leading leg chain" (index "1") and the other legs are termed as "following leg chains" (index "*j*").

The translational part of the constraints is not coupled by the platform orientation and therefore left unchanged relative to (16) with

$$\Phi\_{\mathbf{t}, \mathbf{j}}(\boldsymbol{\varrho}\_{\mathbf{j}}, \mathbf{x}) = -\mathbf{}\_{(0)} \mathbf{r}\_{D} + \, \_{(0)} \mathbf{r}\_{\mathbb{E}}(\boldsymbol{\varrho}\_{\mathbf{j}}) = -\mathbf{x}\_{\mathbb{t}} + \, \_{(0)} \mathbf{r}\_{\mathbb{E}}(\boldsymbol{\varrho}\_{\mathbf{j}}) \in \mathbb{R}^{3} \tag{25}$$

for the following legs *j*.

**Figure 5.** Different cases for the kinematic constraints of the following chain: (**a**) wrong position and orientation; (**b**) correct position and wrong orientation; (**c**) all kinematic constraints are complied.

The orientation for the platform is given with the rotation matrix

$$\prescript{0}{}{\mathcal{R}}\_L(q\_1) := \prescript{0}{}{\mathcal{R}}\_E(q\_1) \tag{26}$$

which gives the reference end-effector frame F*<sup>L</sup>* resulting from the leading leg 1. The rotational part of the kinematic constraints

$$\boldsymbol{\Phi}\_{\mathbf{r},\mathbf{j}}(\boldsymbol{q}\_{\mathbf{j}},\boldsymbol{q}\_1) = \boldsymbol{\mathfrak{a}}(\prescript{0}{}{\mathbf{R}}\_L^{\mathrm{T}}(\boldsymbol{q}\_1)^0 \mathbf{R}\_E(\boldsymbol{q}\_{\mathbf{j}})) \tag{27}$$

for the following leg is given by the Euler angle representation of the deviation between the orientation of the platform frame F*<sup>L</sup>* given by the leading ("L") leg and the frame F*<sup>E</sup>* given by the respective following leg *j*. The choice of the Euler angle convention is arbitrary. The full kinematic constraints for the complete parallel robot with *m* legs for 3T2R tasks

$$\mathbf{Y} = \begin{bmatrix} \mathbf{Y}\_1^T & \mathbf{\Phi}\_2^T & \cdots & \mathbf{\Phi}\_m^T \end{bmatrix}^T \tag{28}$$

are assembled from the 3T2R constraints *Ψ*<sup>1</sup> from (23) for the leading leg and the 3T3R constraints *Φj*, 2 ≤ *j* ≤ *m* from (25,27) for the following legs. The index "*j*" is used to distinguish the following legs of the 3T2R case and all legs "*i*" of the general case in Section 3. The full constraints of (28) lead to a 35-dimensional vector for the kinematic constraints for parallel robots with six legs in 3T2R tasks, which are aggregated as class I in Section 1.2. This formulation can be reduced by combining the mechanism-specific approach for the constraints from (6) or (8) with the principle of leading and following legs of this section. For the 6UPS structure this would result in a 10-dimensional constraint vector with five entries for the leading leg and only one entry for each of the five following legs.

### **5. Differential Kinematics for Parallel Robots**

To be able to compute the differential kinematics of the constraints *Φ* (24) and *Ψ* (28) to

$$\frac{d}{dt}\Phi(q,\mathbf{x}) = \Phi\_{\partial q}\dot{q} + \Phi\_{\partial \mathbf{x}}\dot{\mathbf{x}} = \mathbf{0} \quad \text{and} \quad \frac{d}{dt}\Psi(q,\eta) = \Psi\_{\partial q}\dot{q} + \Psi\_{\partial \eta}\dot{\eta} = \mathbf{0},\tag{29}$$

the full inverse-kinematics matrices

$$\boldsymbol{\Phi}\_{\dot{\mathbf{q}}q}(\boldsymbol{q},\boldsymbol{x}) = \begin{bmatrix} \boldsymbol{\Phi}\_{1,\dot{\boldsymbol{\lambda}}q\_{1}} & \mathbf{0} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \boldsymbol{\Phi}\_{2,\dot{\boldsymbol{\lambda}}q\_{2}} & \ddots & \mathbf{0} \\ \vdots & \ddots & \ddots & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \boldsymbol{\Phi}\_{m,\dot{\boldsymbol{\lambda}}q\_{m}} \end{bmatrix} \quad \text{and} \quad \mathbf{\varPsi}\_{\dot{\mathbf{z}}q}(\boldsymbol{q},\boldsymbol{\eta}) = \begin{bmatrix} \boldsymbol{\Psi}\_{1,\dot{\boldsymbol{\lambda}}q\_{1}} & \mathbf{0} & \mathbf{0} & \mathbf{0} \\ \boldsymbol{\Phi}\_{2,\dot{\boldsymbol{\lambda}}q\_{1}} & \boldsymbol{\Phi}\_{2,\dot{\boldsymbol{\lambda}}q\_{2}} & \ddots & \mathbf{0} \\ \vdots & \ddots & \ddots & \mathbf{0} \\ \boldsymbol{\Phi}\_{m,\dot{\boldsymbol{\lambda}}q\_{1}} & \mathbf{0} & \mathbf{0} & \boldsymbol{\Phi}\_{m,\dot{\boldsymbol{\lambda}}q\_{m}} \end{bmatrix} \tag{30}$$

and the full direct-kinematics matrices

$$\boldsymbol{\Phi}\_{\partial \mathbf{x}}(\boldsymbol{q}, \mathbf{x}) = \frac{\partial \boldsymbol{\Phi}}{\partial \mathbf{x}} = \begin{bmatrix} \boldsymbol{\Phi}\_{1, \text{idx}} \\ \boldsymbol{\Phi}\_{2, \text{idx}} \\ \vdots \\ \boldsymbol{\Phi}\_{m, \text{idx}} \end{bmatrix} \quad \text{and} \quad \boldsymbol{\Psi}\_{\partial \boldsymbol{\eta}}(\boldsymbol{q}, \boldsymbol{\eta}) = \frac{\partial \boldsymbol{\Psi}}{\partial \boldsymbol{\eta}} = \begin{bmatrix} \boldsymbol{\Psi}\_{1, \text{\hat{d}} \boldsymbol{\eta}} \\ \boldsymbol{\Phi}\_{2, \text{\hat{d}} \boldsymbol{\eta}} \\ \vdots \\ \boldsymbol{\Phi}\_{m, \text{\hat{d}} \boldsymbol{\eta}} \end{bmatrix} \tag{31}$$

have to be calculated for the 3T3R and the 3T2R case respectively. The gradient *Φ∂<sup>q</sup>* has block diagonal form, indicating that the inverse kinematics problem can be solved for each leg independently. The structures of the gradient *Ψ∂<sup>q</sup>* results from the coupling of the leading and following joints in the rotational constraints equation.

The gradient matrices *Φ∂<sup>q</sup>* and *Φ∂<sup>x</sup>* contain nested nonlinear functions related to the orientation error. They are calculated with the chain rule and a syntax for stacking matrix columns to avoid differentiating matrices or with respect to matrices, which was introduced in [5]. The product operator

Π, the stacking operator *R* and the transpose operator *P*<sup>T</sup> used for the implementation in the next section are explained in Appendix A.4.

### *5.1. Constraint Gradients for the Leading Leg of the 3T2R and All Legs of the 3T3R Case*

The constraint definition *Ψ*<sup>1</sup> for the leading leg of the 3T2R case (28) and *Φ<sup>i</sup>* with *i* = 1, ..., *m* for all legs of the 3T3R case (24) are subject to the same model of (16), (19), (18). In the following, the 3T3R constraints are displayed. The form *Ψ*<sup>1</sup> for the 3T2R case is obtained by removing the corresponding line of the rotational component according to (22) and replacing "*Φ*" by "*Ψ*" in the following equations. For the analysis, the constraint gradient matrix w. r. t. the joint coordinates has to be divided out to

$$\boldsymbol{\Phi}\_{1,\bar{\partial}q\_{1}} = \begin{bmatrix} \boldsymbol{\Phi}\_{\mathrm{t},1,\bar{\partial}q\_{1}}^{\mathrm{T}} & \boldsymbol{\Phi}\_{\mathrm{r},1,\bar{\partial}q\_{1}}^{\mathrm{T}} \end{bmatrix}^{\mathrm{T}},\tag{32}$$

where the translational component can be calculated with the geometric Jacobian of the leg chain, as derived in Appendix A.3. The rotational part is written down as a function composition of the three functions *α* (Euler angles), ∏ (matrix product) and <sup>0</sup>*R<sup>E</sup>* (rotation matrix) as

$$\boldsymbol{\Phi}\_{\mathrm{r},1;\bar{\boldsymbol{\alpha}}q\_{1}} = \frac{\partial}{\partial q\_{1}}\mathfrak{a}\left({}^{0}\boldsymbol{\mathcal{R}}\_{D}^{\mathrm{T}}(\mathbf{x}){}^{0}\boldsymbol{\mathcal{R}}\_{\mathrm{E}}(q\_{1})\right) = \frac{\partial}{\partial q\_{1}}\mathfrak{a}\left(\prod\left({}^{0}\boldsymbol{\overline{\mathcal{R}}}\_{D}^{\mathrm{T}}(\mathbf{x}),{}^{0}\boldsymbol{\overline{\mathcal{R}}}\_{\mathrm{E}}(q\_{1})\right)\right),\tag{33}$$

which is then expanded with the chain rule for differentiation and the stack operators to

$$\Phi\_{\mathbf{r},1\cdot \mathbf{\hat{z}}q\_{1}} = \underbrace{\frac{\partial \mathbf{a}}{\partial \mathbf{\overline{R}}}}\_{\mathbf{I} \in \mathbb{R}^{3 \times 9}} \underbrace{\frac{\partial \prod \left(0 \, \overline{\mathbf{R}}\_{D}^{\mathsf{T}} \, ^{0} \mathbf{\overline{R}}\_{E}\right)}{\partial ^{0} \overline{\mathbf{R}}\_{E}}}\_{\mathbf{II} \in \mathbb{R}^{9 \times 9}} \underbrace{\frac{\partial \, \mathbf{\overline{R}}\_{E}(q\_{1})}{\partial q\_{1}}}\_{\mathbf{III} \in \mathbb{R}^{9 \times \text{dim}(q\_{1})}} \in \mathbb{R}^{3 \times \text{dim}(q\_{1})}.\tag{34}$$

The two first partial derivatives from (34) are sparse matrices and can be calculated efficiently as shown in Appendix A.4. The factor "I" contains (A24) with *<sup>R</sup>* = *<sup>D</sup>RE*(*x*r, *<sup>q</sup>*1) and the factor "II" is (A27) with the contents of <sup>0</sup>*R*<sup>T</sup> *<sup>D</sup>*(*x*r). The last partial derivative "III" can be derived with computer algebra systems from the analytic expression of the rotation matrix <sup>0</sup>*RE*(*q*1) or from the geometric Jacobian, see Appendix A.2. The leading legs constraint gradient matrix w. r t. the platform coordinates can be expanded in the same manner into

$$
\boldsymbol{\Phi}\_{\mathrm{l},\bar{\mathrm{}}\mathbf{x}} = \begin{bmatrix}
\boldsymbol{\Phi}\_{\mathrm{t},\mathrm{l},\bar{\mathrm{}}\mathbf{x}\_{\mathrm{l}}} & \boldsymbol{\Phi}\_{\mathrm{t},\mathrm{l},\bar{\mathrm{}}\mathbf{x}\_{\mathrm{r}}} \\
\boldsymbol{\Phi}\_{\mathrm{r},\mathrm{l},\bar{\mathrm{}}\mathbf{x}\_{\mathrm{t}}} & \boldsymbol{\Phi}\_{\mathrm{r},\mathrm{l},\bar{\mathrm{}}\mathbf{x}\_{\mathrm{r}}}
\end{bmatrix} = \begin{bmatrix}
\mathbf{0} & \boldsymbol{\Phi}\_{\mathrm{r},\mathrm{l},\bar{\mathrm{}}\mathbf{x}\_{\mathrm{r}}}
\end{bmatrix},\tag{35}
$$

where the definitions from (16) and (19) only leave the rotational part

$$\begin{split} \boldsymbol{\Phi}\_{\mathrm{r},1,\partial\mathbf{x}\_{\mathrm{r}}} &= \frac{\partial}{\partial\mathbf{x}\_{\mathrm{r}}} \mathbf{a} \left( \left( {}^{0}\boldsymbol{\mathsf{R}}\_{E}^{\mathrm{T}}(\boldsymbol{q}\_{1}) {}^{0}\boldsymbol{\mathsf{R}}\_{D}(\mathbf{x}\_{\mathrm{r}}) \right)^{\mathrm{T}} \right) = \frac{\partial}{\partial\mathbf{x}\_{\mathrm{r}}} \overline{\mathbf{a}} \left( {}^{0}\overline{\mathbf{R}}\_{E}^{\mathrm{T}}(\boldsymbol{q}\_{1}), {}^{0}\overline{\mathbf{R}}\_{D}(\mathbf{x}\_{\mathrm{r}}) \right) \Big| \\ &= \underbrace{\frac{\partial\overline{\mathbf{a}}}{\partial\overline{\mathbf{R}}}}\_{\mathbf{I} \in \mathbb{R}^{3 \times 9}} \underbrace{\mathcal{P}\_{\mathrm{T}}}\_{\Pi \in \mathbb{R}^{9 \times 9}} \underbrace{\frac{\partial\overline{\mathbf{I}} \left( {}^{0}\overline{\mathbf{R}}\_{E}^{\mathrm{T}} {}^{0}\overline{\mathbf{R}}\_{D} \right)}\_{\Pi \in \mathbb{R}^{9 \times 9}} \underbrace{\partial^{0} \overline{\mathbf{R}}\_{D}(\mathbf{x}\_{\mathrm{r}})}\_{\Pi^{0} {}^{0}\overline{\mathbf{R}}\_{\Gamma}} \in \mathbb{R}^{3 \times 3} \right. \end{split} \tag{36}$$

where the simplicity of the single expression "I"-"IV" is demonstrated in Appendix A.4. The factors are (A24) with *<sup>R</sup>* = *<sup>D</sup>RE*(*x*r, *<sup>q</sup>*1) in "I", the permutation matrix for transposition from (A23) in "II", (A27), where the contents of <sup>0</sup>*R*<sup>T</sup> *<sup>E</sup>* are inserted in "III" and (A25) with the elements of *x*<sup>r</sup> for *β* in "IV".

### *5.2. Constraint Gradients for the Following Leg in the 3T2R Case*

As explained regarding (24), the constraints (16) and (19) and their gradients (34) and (36) are used for all legs in the 3T3R case and the leading leg in the 3T2R case. For the following legs in the 3T2R case, the gradients *Φj*,*∂q*<sup>1</sup> , *Φj*,*∂q<sup>j</sup>* and *Φj*,*∂<sup>x</sup>* with *j* = 2, ..., *m* from the right part of (30) and (31) have to be calculated in a similar way. Due to the absence of the platform orientation in (27), (35) simplifies for the following leg to 

$$
\boldsymbol{\Phi}\_{j, \partial \boldsymbol{x}} = \begin{bmatrix} -\mathbf{1} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} \end{bmatrix}. \tag{37}
$$

The gradient w. r. t. the joint coordinates of the following leg contains again the translational part of the legs Jacobian regarding the end-effector platform position in *Φ*t,*j*,*∂q<sup>j</sup>* and has the rotational part

$$\begin{split} \boldsymbol{\Phi}\_{\mathrm{r},j,\mathrm{d}q\_{j}} &= \frac{\partial}{\partial q\_{j}} \boldsymbol{\mathfrak{a}} \left( {}^{0}\mathbf{R}\_{L}^{\mathrm{T}}(q\_{1}) {}^{0}\mathbf{R}\_{E}(q\_{j}) \right) = \frac{\partial}{\partial q\_{j}} \boldsymbol{\mathfrak{a}} \left( \prod \left( {}^{0}\overline{\mathbf{R}}\_{L}^{\mathrm{T}}(q\_{1}) {}^{0}\overline{\mathbf{R}}\_{E}(q\_{j}) \right) \right) \\ &= \underbrace{\frac{\partial \boldsymbol{\mathfrak{a}}}{\partial \overline{\mathbf{R}}}}\_{\mathbf{I} \in \mathbb{R}^{3 \times 9}} \underbrace{\frac{\partial \overline{\mathbf{T}} \left( {}^{0}\overline{\mathbf{R}}\_{D}^{\mathrm{T}} {}^{0}\overline{\mathbf{R}}\_{E} \right)}\_{\mathbf{II} \in \mathbb{R}^{9 \times 9}} \underbrace{\partial^{0} \overline{\mathbf{R}}\_{E}(q\_{j})}\_{\mathbf{III} \in \mathbb{R}^{9 \times \mathrm{dim}(q\_{j})}} \in \mathbb{R}^{3 \times \mathrm{dim}(q\_{j})} \end{split} \tag{38}$$

which is similar to the expression in (34). The factors of Equation (38) are (A24) with *<sup>R</sup>* = *<sup>L</sup>RE*(*q*1, *<sup>q</sup><sup>j</sup>* ) in "I", (A27), where the elements of <sup>0</sup>*R*<sup>T</sup> *<sup>L</sup>*(*q*1) have to be inserted in "II" and the partial derivative of the platform orientation calculated from leg *j* w. r. t. the legs joint coordinates in "III", similar to term "III" from (34). The gradient w. r. t. the joint coordinates of the leading leg

$$\begin{split} \boldsymbol{\Phi}\_{\mathrm{r},j,\mathbbm{dq}\_{1}} &= \frac{\partial}{\partial q\_{1}} \boldsymbol{\mathfrak{a}} \left( \left( {}^{0}\mathbf{R}\_{\mathrm{E}}^{\mathrm{T}}(q\_{j}){}^{0}\mathbf{R}\_{\mathrm{L}}(q\_{1}) \right)^{\mathrm{T}} \right) = \frac{\partial}{\partial q\_{1}} \boldsymbol{\mathfrak{a}} \left( {}^{0}\mathbf{T}\_{\mathrm{T}} \overline{\prod} \left( {}^{0}\mathbf{R}\_{\mathrm{E}}^{\mathrm{T}}(q\_{j}){}^{0}\overline{\mathbf{R}}\_{\mathrm{L}}(q\_{1}) \right) \right) \\ &= \underbrace{\frac{\partial \overline{\mathbf{a}}}{\partial \overline{\mathbf{R}}}}\_{\mathrm{I} \in \mathbb{R}^{3 \times 9}} \underbrace{\mathrm{P}\_{\mathrm{T}}}\_{\mathrm{II} \in \mathbb{R}^{9 \times 9}} \underbrace{\frac{\partial \overline{\prod} \left( {}^{0}\overline{\mathbf{R}}\_{\mathrm{E}}^{\mathrm{T}} {}^{0}\overline{\mathbf{R}}\_{\mathrm{L}} \right)}\_{\mathrm{III} \in \mathbb{R}^{9 \times 9}}}\_{\mathrm{II} \in \mathbb{R}^{9 \times \mathrm{dim}(q\_{1})}} \underbrace{\frac{\partial {\partial} \overline{\mathbf{R}}\_{\mathrm{L}}(q\_{1})}{\partial q\_{1}}}\_{\mathrm{IV} \in \mathbb{R}^{3 \times \mathrm{dim}(q\_{1})}} \in \mathbb{R}^{3 \times \mathrm{dim}(q\_{1})}. \end{split} \tag{39}$$

is similar to the expression in (36). The order of the residual expression (27) has to be switched by exploiting the associative property of matrix transposition to avoid differentiating a transposed matrix. The factors of equation (39) are (A24) with *<sup>R</sup>* = *<sup>L</sup>RE*(*q*1, *<sup>q</sup><sup>j</sup>* ) in "I", the permutation matrix for transposition from (A23) in "II", (A27) with <sup>0</sup>*R*<sup>T</sup> *<sup>E</sup>*(*q<sup>j</sup>* ) in "III" and the term "III" from (34) in "IV".

### *5.3. Gradient-Based Solution of the Inverse Kinematics Problem with Redundancy Resolution*

The presented kinematic constraints and their gradient matrices can be used to solve the inverse kinematics problem (IKP) of single leg chains and complete parallel robots. Since all active and passive joint angles are involved for the case of parallel robots, solving their IKP results in solving the IKP for all leg chains. As first introduced in [7] for Euler angle residuals in the IKP, the Taylor series expansion of *Φ*(*q*, *x*) leads to the definition of

$$\Phi(q^{k+1}, \mathbf{x}) = \Phi(q^k, \mathbf{x}) + \frac{\partial}{\partial q} \Phi(q, \mathbf{x}) \Big|\_{q^k} (q^{k+1} - q^k) \tag{40}$$

in an iterative algorithm at the step *k* + 1, which can be used to solve the IKP using a given initial value *q*<sup>0</sup> and the condition

$$\Phi(q^{k+1}, \mathbf{x}) = \mathbf{0}.\tag{41}$$

Defining the solution of the IKP as the main task ("T"), the stepwise solution for the joint coordinates results to

$$
\Delta \boldsymbol{q}^k = \Delta \boldsymbol{q}\_{\mathbb{T}}^k = \boldsymbol{q}^{k+1} - \boldsymbol{q}^k = -\boldsymbol{\Phi}\_{\partial \boldsymbol{q}}^\dagger(\boldsymbol{q}^k, \boldsymbol{\mathsf{x}}) \boldsymbol{\Phi}(\boldsymbol{q}^k, \boldsymbol{\mathsf{x}}).\tag{42}
$$

Depending on the dimension, (·)† denotes the matrix inverse or the pseudo-inverse. Again, the Equations (40)–(42) can be written with "*Ψ*" from (28) instead of "*Φ*" from (24) for the 3T2R case.

In the latter case, the corresponding gradient matrix *Ψ∂q*(*q*, *η*) from (30) allows defining of a nullspace in the case of dim(*q*1) > dim(*η*). This redundancy can be exploited by using the nullspace ("N") projection Δ*q*<sup>N</sup> from [6] additionally to the solution Δ*q*<sup>T</sup> of the IKP in (42) with the new increment

$$
\Delta \boldsymbol{q}^{k} = \boldsymbol{q}^{k+1} - \boldsymbol{q}^{k} = \Delta \boldsymbol{q}\_{\rm T}^{k} + \Delta \boldsymbol{q}\_{\rm N}^{k} = -\boldsymbol{\Psi}\_{\partial q}^{\rm t} \boldsymbol{\Psi} + (\mathbf{1} - \mathbf{Y}\_{\partial q}^{\rm t} \boldsymbol{\Psi}\_{\partial q}) \boldsymbol{h}\_{\partial q} \tag{43}
$$

in the iterative algorithm. The optimization of additional performance criteria *h* requires their gradient *h∂<sup>q</sup>* w. r .t the joint positions. One criterion is the summed *W*1-weighted quadratic distance

$$h\_1(q) = \frac{1}{2}(q - \bar{q})^T \mathcal{W}\_1(q - \bar{q}), \quad h\_{1\bar{\beta}q} = \frac{\partial h\_1}{\partial q} = \mathcal{W}\_1(q - \bar{q}) \tag{44}$$

of the joint positions *q* from their respective reference position *q*¯, e.g., used in [2,36]. Defining *q*¯ to be in the middle of the joint limits and minimizing *h*1(*q*) reduces the risk of joints reaching their technical limits, but does not guarantee it, since exceeding the limit for one joint can be compensated by improving other joints. The *W*2-weighted hyperbolic joint limit distance

$$h\_2(q) = \frac{1}{n} \sum\_{i=1}^n w\_{2,i} \frac{(q\_{i,\text{max}} - q\_{i,\text{min}})}{8} \left( \frac{1}{(q\_i - q\_{i,\text{min}})^2} + \frac{1}{(q\_i - q\_{i,\text{max}})^2} \right) \tag{45}$$

from [8] (written element-wise for *n* = dim(*q*)) circumvents this problem by generating infinitely high values when reaching the limits. In contrast to *h*1, the criterion *h*<sup>2</sup> is only defined for joints within their limits with *qi*,min < *qi* < *qi*,max, which is ensured by setting *w*2,*<sup>i</sup>* = 0 for joints exceeding their limits and *w*2,*<sup>i</sup>* = 1 otherwise. To combine the effect of drawing joint positions to their middle with *h*<sup>1</sup> of (44) and of strongly rejecting joints directly near their limits with *h*<sup>2</sup> of (45), their weighted sum

$$h\_3(q) = \mathcal{K}\_{h\_1} h\_1(q) + \mathcal{K}\_{h\_2} h\_2(q) \tag{46}$$

is used in the simulation studies of Section 6. Other criteria not related to the joint limits are for example stiffness [9] or singularity avoidance via Frobeniusnorm condition number [8] or squared condition number [4]. The method can be used for serial-link robots as well by removing all entries for the following legs from the formulas, as presented in [5]. The platform pose *x*/*η* corresponds to the desired pose for the serial robots end-effector and the kinematic constraints *Φ*/*Ψ* correspond to the residual of the IKP.

In the practical implementation, it has proven to be useful to extend the basic principle of (43) to

$$
\Delta q^k = K\_{\rm lim}(q^k) K\_{\rm Rel}(q^k) (K\_{\rm T} \Delta q\_{\rm T}^k + K\_{\rm N} \Delta q\_{\rm N}^k), \tag{47}
$$

where the constant damping coefficients *K*<sup>T</sup> for Δ*q<sup>k</sup>* <sup>T</sup> and *<sup>K</sup>*<sup>N</sup> for <sup>Δ</sup>*q<sup>k</sup>* <sup>N</sup> were introduced to avoid overshooting of the solution for the prize of slower convergence. The damping term *K*<sup>N</sup> has to be chosen according to the optimization criterion. Further damping was introduced for the 3T2R case with task redundancy to reduce a Δ*q<sup>k</sup>* that would lead to overshoot over the joint limits with *K*Lim(*qk*). The value *K*Lim = 1 is set if no limits would be violated by the increment Δ*qk*. For the 3T3R case, *K*Lim := 1 is set permanently, since slowing down when approaching the limits does not change the direction of the increment and violating the limits is inevitable. The maximum step size for one iteration Δ*q<sup>k</sup>* was ensured with *K*Rel(*qk*) to stay below 5% of the joint limit range to prevent leaving the validity of the first-order approximation of (40). For smaller increments, *K*Rel = 1 holds. The damping terms are always applied to the full vector and not to single elements and therefore only change the norm and not the direction of Δ*qk*.

### *5.4. Differential Kinematics for the Parallel Robot and Its Applications*

The reasoning so far only considered the inverse kinematics of the parallel robot. The kinematic definitions can also be used in the differential kinematics (29) to establish the connection between joint and platform velocity. This was already presented in general form in [15] and also corresponds to the theory of linear transformation which is the base of the works of Gogu on structural synthesis [11]. The derivation in this paper is based on the position *and* orientation, but comes to the same result as the already-existing velocity-based approach. Furthermore, using *Ψ*/*η* as elaborated before allows for the first time defining differential kinematics specifically for 3T2R tasks in a general form. The differential equation of (29) is expanded to

$$\frac{d}{dt}\Phi(q\_{\rm a},q\_{\rm p},\mathbf{x}) = \Phi\_{\rm \partial q\_{\rm a}}\dot{q}\_{\rm a} + \Phi\_{\rm \partial q\_{\rm p}}\dot{q}\_{\rm p} + \Phi\_{\rm \partial x}\dot{\mathbf{x}} = \mathbf{0} \tag{48}$$

to distinguish active ("a") and passive ("p") joints. The latter also contain the coordinates of the platform-connecting joints. Reordering the equation leads to the full inverse differential kinematics

$$
\boldsymbol{\Phi}\_{\partial \mathbf{x}} \dot{\mathbf{x}} = - \begin{bmatrix} \boldsymbol{\Phi}\_{\partial q\_{\mathbf{a}}} & \boldsymbol{\Phi}\_{\partial q\_{\mathbf{p}}} \end{bmatrix} \begin{bmatrix} \dot{q}\_{\mathbf{a}} \\ \dot{q}\_{\mathbf{p}} \end{bmatrix} = \boldsymbol{\Phi}\_{\partial \mathbf{a} \mathbf{p}} \begin{bmatrix} \dot{q}\_{\mathbf{a}} \\ \dot{q}\_{\mathbf{p}} \end{bmatrix}, \qquad \begin{bmatrix} \dot{q}\_{\mathbf{a}} \\ \dot{q}\_{\mathbf{p}} \end{bmatrix} = - \boldsymbol{\Phi}\_{\partial \mathbf{a} \mathbf{p}}^{-1} \boldsymbol{\Phi}\_{\partial \mathbf{x}} \mathbf{x}, \tag{49}
$$

which has been addressed in the previous sections, and the full direct differential kinematics

$$
\boldsymbol{\Phi}\_{\partial q\_{\boldsymbol{\mathfrak{a}}}} \dot{\boldsymbol{q}}\_{\boldsymbol{\mathfrak{a}}} = - \begin{bmatrix} \boldsymbol{\Phi}\_{\partial \mathbf{x}} & \boldsymbol{\Phi}\_{\partial q\_{\boldsymbol{\mathfrak{p}}}} \end{bmatrix} \begin{bmatrix} \dot{\boldsymbol{\mathfrak{x}}} \\ \dot{\boldsymbol{q}}\_{\boldsymbol{\mathfrak{p}}} \end{bmatrix} = \boldsymbol{\Phi}\_{\partial \mathbf{x}\mathbf{p}} \begin{bmatrix} \dot{\boldsymbol{\mathfrak{x}}} \\ \dot{\boldsymbol{q}}\_{\boldsymbol{\mathfrak{p}}} \end{bmatrix}, \qquad \begin{bmatrix} \dot{\boldsymbol{\mathfrak{x}}} \\ \dot{\boldsymbol{q}}\_{\boldsymbol{\mathfrak{p}}} \end{bmatrix} = - \boldsymbol{\Phi}\_{\partial \mathbf{x}\mathbf{p}}^{-1} \boldsymbol{\Phi}\_{\partial q\_{\boldsymbol{\mathfrak{a}}}} \dot{\boldsymbol{q}}\_{\boldsymbol{\mathfrak{a}}}.\tag{50}
$$

By only selecting the first rows for *q*˙ <sup>a</sup> in (49) and for *x*˙ in (50), the well-known analytic Jacobian of the parallel robot (called "Euler angles Jacobian matrix" in [15] and "design Jacobian" in [11]), relating actuator velocities *q*˙ <sup>a</sup> and platform velocities *x*˙ , can be obtained from both equations. For the case of task or kinematic redundancy, the pseudo-inverse can be used for *Φ∂*ap in (49) together with a nullspace optimization as shown in Section 5.3. This allows exploiting the functional redundancy in 3T2R tasks also for trajectories *η*(*t*), *η*˙(*t*) and not only for single poses *η*. When solving the IKP for a trajectory, the IK on position level of Section 5.3 is only needed to correct linearization errors. The case of task redundancy does not affect (50), since the full platform velocity *x*˙ is obtained from given actuator velocities ˙*q*a.

### **6. Results**

To evaluate the inverse kinematics (IK) algorithm presented in the previous Section 5.3, first the solution of the IKP is shown for the trajectory of a serial-link industrial robot in Section 6.1 and for the trajectory of a parallel robot in Section 6.2. The results are generalized by the statistical analysis of random point-to-point movements of arbitrary serial link chains in Section 6.3.

### *6.1. Resolution of Functional Redundancy of a Serial-Link Six-DoF Robot in 3T2R tasks*

The first evaluation of the IK algorithm from Section 5.3 is performed with simulations at the basic example of a six-DoF industrial robot with a 500 mm × 800 mm rectangular trajectory centered at 1200 mm/−200 mm/200 mm with a constant orientation of the *z*-axis pointing into the ground plane. The manipulator Fanuc M-710 iC/50 was taken from the example of [36] with the tabulated kinematics parameters and a sketch of the trajectory in Figure 6. Deviations in the parameters relative to [36] result from the use of the *modified* Denavit–Hartenberg notation from [35] for the joint transformations and the use of only positive axis alignment parameters *α<sup>i</sup>* for consistency with the results of the structural synthesis from [34].


**Figure 6. Left**: Table with the kinematic parameters of the industrial manipulator Fanuc M-710 iC/50. **Right**: sketch of the robot scenario.

The IKP is solved with two settings: setting the tool axis rotation to different constant values *β*<sup>3</sup> with the 3T3R algorithm and solving the IKP only for the desired pointing direction with the 3T2R algorithm. The algorithm from (43) was used in the extended version of (47) for both cases with different settings caused by their nature. For the 3T3R case, *K*<sup>T</sup> = 0.7 and *K*<sup>N</sup> = 0.7 were set. The terms *K*<sup>N</sup> and *K*Lim have no effect, since no nullspace movement is possible. For the 3T2R case, with *Kh*<sup>1</sup> = 0 and *Kh*<sup>2</sup> = 1 only the hyperbolic limit rejection criterion from (45) was used. The first criterion was not used, since in the trajectory example the limits are not even temporarily exceeded by principle. All IKP algorithms had the same initial value from Figure 6.

The results of the inverse kinematics for different settings are given in Figure 7, where the representative joint coordinates *q*<sup>1</sup> and *q*5, the redundant coordinate of the end-effector orientation *β*3, as well as the optimization criterion (45) are depicted over time for the trajectory from Figure 6. The positions are normalized to the joint limits from −1 to 1. The first three lines in Figure 7 represent IKP solutions with a given constant end-effector orientation *β*<sup>3</sup> of −150◦, −15◦ and 45◦ and the 3T3R algorithm. The 3T2R algorithm without nullspace optimization is plotted with dotted lines for each first sample of the 3T3R cases as initial value with the same colors. Using these initial values for a 3T2R IK with optimization leads to strong nullspace movements at the beginning, quickly converging to a local minimum. Therefore, the 3T2R case with optimization, plotted as the green line with triangle markers, is shown only for the initial condition from Figure 6. It can be observed that the optimization of the criterion leads to the best solution of the IKP. The lines for the criterion for *β*<sup>3</sup> = −150◦ and *β*<sup>3</sup> = −15◦ partly exceed the limits of the plot, indicating that the limit is violated, which can also be seen at the plot for *q*5. This exposes the need for keeping the solution always within the limits by the measures described. The 3T2R IK without optimization with dotted lines tends to lower changes in the joint positions than the 3T3R IK, since this corresponds to the solution of the matrix pseudo-inverse in (43).

**Figure 7.** Results of the inverse kinematics with different settings for the trajectory of Figure 6.

### *6.2. Resolution of Functional Redundancy of a Parallel Robot in 3T2R Tasks*

As elaborated in Section 4 and 5, the solution of the IKP for 3T2R and 3T3R tasks is necessary to solve the problem for parallel robots. Therefore, the trajectory evaluation for a 6UPS parallel robot in this section is preceded by the trajectory example for a serial link chain in the previous section. This robot belongs to the first class presented in Section 1.2, which is primarily addressed in this paper.

The robot has a Goughstructure [15] with symmetric alignment of the universal joint base couplings on a circle with radius *r*0*Ai* = 1 m and the spherical joint platform couplings on a circle with radius *rBiE* <sup>=</sup> 0.4 m. The initial pose was set to a center position *<sup>x</sup>*<sup>T</sup> <sup>t</sup> = [0, 0, 0.5 m] and the initial orientation *x*r was set to zero, meaning an alignment of base and platform frame. The joint positions for each leg were defined to have the initial values *q*<sup>T</sup> *<sup>i</sup>* = [30◦, −30◦, 0.583 m, 0◦, 30◦, 60◦] for the given initial platform pose to avoid switching ±*π* within the trajectory and to avoid gimballocksingularities. The joint limits were set around the resulting zero position to ±0.5 m for the prismatic joint and ±60◦ for all single revolute joints representing the universal and spherical joints. The values are higher than typical values for real robots to emphasize the effect of the nullspace movement in a bigger simulated workspace of the robot. The settings for the IK solver are similar to those in Section 6.1, since both cases regard solving the IKP for a trajectory.

The time evolution of platform pose and optimization criteria is depicted in Figure 8. The reference trajectory can be seen at the platform position in Figure 8a and the platform orientation expressed in *XYZ* Euler angles (*β*1-*β*3) relative to the base frame in Figure 8b. The IKP is solved using two different methods: Only solving the IKP for the legs separately, called "ser. IK" in Figure 8 and solving the IKP for all legs together, called "par. IK" in Figure 8. Both methods perform an optimization with only *h*<sup>2</sup> of (45), as justified in Section 6.1. The first approach only performs this optimization according to Section 3 for the first leg using the 3T2R method and then solves the IKP for all other legs with the 3T3R method. The second approach uses the optimization for all legs together according to the 3T2R method from Section 4. This results in improved values for the performance criteria depicted for *h*<sup>1</sup> in Figure 8c and for *h*<sup>2</sup> in a logarithmic scale in Figure 8d. Since the first approach does not regard the limits of the following legs, the optimization criterion gives high values indicating many joint limit violations. The second approach only shows peaks at *t* = 1.5 s in Figure 8d that result from a joint position getting near to the limit, but not exceeding it. For the practical implementation, the computation time is only weakly influenced by the selection of the method, since calculating the (pseudo-)inverse for six 5 × 6 and 6 × 6 matrices or one 35 × 36 matrix does not present a challenge for current computing hardware. Therefore, the "par. IK"method should be preferred.

**Figure 8.** Results of the inverse kinematics of a 6UPS robot in a 3T2R task. (**a**) platform positions, (**b**) platform orientation in Euler angles, (**c**,**d**) optimization criteria *h*1(*q*) and *h*2(*q*).

### *6.3. Statistic Results for the Inverse Kinematics of Serial Link Chains*

To emphasize the generality of the presented approach, the inverse kinematics is solved for a set of 309 serial kinematic chains with six joints. This set of six-DoF kinematics is generated by permutations of their Denavit–Hartenberg parameters and is reduced with the isomorphism detection of [34] to a minimal set, representing all possible six-DoF serial kinematics with full mobility. The approach is similar to the results of the evolutionary morphology of parallel robot leg chains of [11]. In contrast to the trajectory evaluations in the previous sections focusing on nullspace movement, the inverse kinematics is solved in this section for arbitrary reachable poses of the serial chain in its individual workspace. Therefore, different settings proved to be necessary, since for point-to-point movements, intermediate steps may be outside of the joint limits. In the trajectory case, the initial value for the IKP of the continuous trajectory is always very close to the desired pose of the next trajectory sample. Preventing the algorithm completely from leaving the allowed joint positions reduces the IK success rate. Therefore, the damping term for limit violation was not used in this evaluation, resulting in a constant *K*Lim(*q*) = 1 in (47). To reach again an allowed configuration when approaching the goal pose from intermediate steps with limit violations, the combined criterion *h*3(*q*) from (46) with *Kh*<sup>1</sup> = 0.99 and *Kh*<sup>2</sup> = 0.01 was used. Further empirically determined values for all different serial chains were the damping coefficients *K*<sup>T</sup> = 0.6 and *K*<sup>N</sup> = 0.01 in (47). Since these values provide good results for all serial chains with random geometric parameters and for random configurations, they can be regarded as a good choice generally.

To create a general evaluation case, the poses for testing the IK algorithm were generated by the forward kinematics of 50 different joint configurations of the chains uniformly distributed between the joint limits of ±*π* for revolute joints and ±0.5 m for prismatic joints. Additionally, the Denavit–Hartenberg parameters were set to 50 different sets of uniformly distributed parameters between 0 and 1 meters or radians resulting in 2500 combinations for each of the 309 chains in total. The initial value *q*<sup>0</sup> for the solution of the IKP of (47) was set to random values from a uniform distribution within the joint limits. The inverse kinematics was calculated for the full pose with the 3T3R algorithm and only using the pointing direction together with the resolution of functional redundancy in the 3T2R algorithm. A maximum of 15 tries with random initial values was allowed to search for a solution of the IKP within the limits. After that, five more tries were allowed to find a solution violating the limits, but presenting a solution of the IKP to be able to distinguish the two cases, which allows further reasoning on the functionality and possible improvements. A success of the IK is defined as a solution within the joint limits.

The aggregated results are presented as histograms in Figure 9 for different settings of the algorithm. The histograms show that for the worst case in 3T2R (3T3R) tasks, the success rate is 87% (6%), marked by the position of the first bars in Figure 9a,c. These results can be vastly improved by setting the initial guess *q*<sup>0</sup> within 20% (w. r. t. the joint limit range) around the pose, from which the desired end-effector pose has been calculated. This improves the worst success rate of all kinematic chains to 98% for 3T2R (Figure 9b) and 95% for 3T3R tasks (Figure 9d).

**Figure 9.** Histograms with cumulated frequency of the IK success for all kinematic chains with different settings: 3T2R tasks (**a**,**b**) vs. 3T3R tasks (**c**,**d**) and arbitrary initial value (**a**,**c**) vs. initial near goal pose.

A detailed investigation on the success rates of all possible serial chains is performed in Figure 10.

**Figure 10.** Detailed Statistics of the success of the inverse kinematics algorithm for 3T2R tasks (**a**,**b**) and 3T3R tasks (**c**,**d**). The Success of the IK solver is shown in different shades of green for increasing numbers of required tries. Different initial values *q*<sup>0</sup> are distinguished in (**a**,**c**) and (**b**,**d**).

The 309 serial kinematics are sorted according to their number of revolute joints and are listed on the horizontal axis of the figure: They contain three revolute joints up to no. 98, four Rjoints up to no. 240, and five Rjoints up to no. 301. The eight structures from 302 to 309 with six Rjoints differ in the parallelism of their joint axes. The first 240 kinematic chains with more than one prismatic joint can be seen as a rather academic example and are listed for the sake of completeness. The most prominent chains are the UPSchain from Section 6.2 at no. 266 and the six-DoF industrial robot from Section 6.1 at no. 309. Each bar represents the stacked relative frequency of the IK result state in percent for one kinematic chain. The result state is defined as the number of tries or the success. All bars add up to 100%, which corresponds to the 2500 configurations per chain. Beginning at the bottom, the number of tries necessary for the solution of the IKP is marked with colors from green to orange. Only cases with a violation of the limits (bright red) or wrong position (dark red) correspond to a failure of the algorithm, which has been addressed in the analysis of Figure 9, representing an aggregated form of Figure 10. The subfigures a–d of Figure 10 correspond to the ones in Figure 9. It can be observed that the quality of the results is clustered according to the kinematic groups. Structures with at most one Pjoint show a considerably better performance of the algorithm with a worst success rate of 97.16% for five Rjoints and 99.36% for six Rjoints for the 3T2R case (a), which can be seen at the very small red top parts of the bars in the corresponding range of the diagram. The worse performance of the 3T3R algorithm, mostly caused by limit violations, can be explained by joints changing their configuration, i.e., from "elbow up" to "elbow down", which causes limit violations but does not affect the 3T3R IK.

### **7. Discussion**

A general kinematics model for parallel robots was introduced to solve the inverse kinematics problem for any kind of parallel robot in tasks with one redundant rotational DoF (3T2R). The prize of the generality of the approach is the increased size of the inverse-kinematics matrices, which is 10 × 11 for robots with a simple UPS structure instead of 6 × 6 for the non-redundant kinematics and grow up to 35 × 36 for general task redundant parallel robots with full mobility. This makes symbolic calculations of the kinematic matrices impossible, allowing only studies on mobility, singularities and other properties of the Jacobian matrix based on numeric calculations. The application of the proposed method can therefore be seen mainly in finding optimal trajectories for task redundant parallel robots in milling or drilling scenarios regarding stiffness, dexterity or joint limits. Due to the performance of the method demonstrated at exemplary cases, an online implementation is possible but has to be proved in future works to converge in real-time conditions for specific machines. The generality of the approach allows using it in a combined structural and dimensional synthesis sketched in Figure 2, extending the purely structural synthesis of parallel robot kinematics from [11,31] to a dimensional synthesis of all structures as shown in [34] for serial robots.

**Author Contributions:** conceptualization, M.S. and S.T.; methodology, M.S.; software, M.S.; validation, M.S.; formal analysis, M.S.; writing—original draft preparation, M.S.; writing—review and editing, M.S., S.T., T.O.; supervision, S.T. and T.O.; project administration, T.O.; funding acquisition, S.T. and T.O.

**Funding:** The financial support from the Deutsche Forschungsgemeinschaft (German Research Foundation, DFG) under grant number OR 196/33-1 is gracefully acknowledged.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **Abbreviations**

The following abbreviations are used in this manuscript:


### **Appendix A. Mathematical Symbols for Reciprocal Euler Angles in Inverse Kinematics**

The following appendix contains additional detailed information about the kinematic constraint formulation of this paper. Appendix A.1 contains a mathematical proof for the properties of reciprocal Euler angles in inverse kinematics, which is only outlined in equ. 18 of [5]. The relations of the partial derivatives to the geometric and analytic Jacobian are derived in Appendixes A.2 and A.3. The matrix operations for the partial derivatives are replicated from [5] in Appendix A.4 and the contents of the single partial derivatives are given in Appendix A.5 to facilitate the understanding and implementation by the reader.

### *Appendix A.1. Proof for the Properties of Reciprocal Euler Angles*

This section derives the effect of the reciprocity of Eulerangles at the example of the kinematics description of Section 3 and the frames of Figure 4b: An end-effector orientation *β* = *x*<sup>r</sup> gives the rotation matrix

$${}^{D}\mathcal{R}\_{E}(\mathfrak{E},\mathfrak{q}) = {}^{0}\mathcal{R}\_{D}^{\mathrm{T}}(\mathfrak{E}){}^{0}\mathcal{R}\_{E}(\mathfrak{q}),\tag{A1}$$

which rotates vectors from the actual end-effector frame F*<sup>E</sup>* to the desired or platform end-effector frame F*D*. Using the *XYZ* Euler angles and exploiting the properties of SO(3) rotation matrices yields

$$^0\mathcal{R}\_D^T(\mathcal{B}) = \mathcal{R}\_z(-\beta\_3)\mathcal{R}\_y(-\beta\_2)\mathcal{R}\_x(-\beta\_1),\tag{A2}$$

as introduced in (5). With an additional rotation −*δ* around the *z*-axis for the desired orientation, the resulting new Euler angles *β* are

$$
\beta\_1' = \beta\_1, \quad \beta\_2' = \beta\_2, \quad \beta\_3' = \beta\_3 - \delta. \tag{A3}
$$

The additional rotation corresponds to the tool axis defined in Section 3 and leads to a new residual orientation error expressed as a rotation matrix

$$\begin{split}{}^{D}\mathbf{R}\_{E}(\boldsymbol{\mathfrak{g}}^{\prime},\boldsymbol{\mathfrak{q}}) &= {}^{0}\mathbf{R}\_{D}^{\mathrm{T}}(\boldsymbol{\mathfrak{g}}^{\prime}){}^{0}\mathbf{R}\_{E}(\boldsymbol{\mathfrak{q}}) \\ &= \left({}^{0}\mathbf{R}\_{D}(\boldsymbol{\mathfrak{g}})\mathbf{R}\_{z}(-\boldsymbol{\delta})\right)^{\mathrm{T}}{}^{0}\mathbf{R}\_{E}(\boldsymbol{\mathfrak{q}}) \\ &= {}^{\mathbf{R}}\boldsymbol{\mathfrak{g}}\_{z}(\boldsymbol{\delta}){}^{0}\mathbf{R}\_{D}^{\mathrm{T}}(\boldsymbol{\mathfrak{g}}){}^{0}\mathbf{R}\_{E}(\boldsymbol{\mathfrak{q}}) \\ &= {}^{\mathbf{R}}\boldsymbol{\mathfrak{g}}\_{z}(\boldsymbol{\delta}){}^{D}\mathbf{R}\_{E}(\boldsymbol{\mathfrak{g}},\boldsymbol{\mathfrak{q}}). \end{split} \tag{A4}$$

The first residual orientation error from (A1) corresponding to *β* is defined as a rotation matrix

$$\prescript{D}{}{\mathcal{R}}\_{E}(\mathcal{J}, \boldsymbol{\mathfrak{q}}) = \begin{bmatrix} n\_{\boldsymbol{x}} & o\_{\boldsymbol{x}} & a\_{\boldsymbol{x}} \\ n\_{\boldsymbol{y}} & o\_{\boldsymbol{y}} & a\_{\boldsymbol{y}} \\ n\_{\boldsymbol{z}} & o\_{\boldsymbol{z}} & a\_{\boldsymbol{z}} \end{bmatrix} \tag{A5}$$

and as a *ZYX* Euler angle representation

$$\mathbf{a} = \begin{bmatrix} a\_1 \\ a\_2 \\ a\_3 \end{bmatrix} = \begin{bmatrix} \arctan2\left(n\_{y'}, n\_x\right) \\ \arctan2\left(-n\_{z'}\sqrt{a\_z^2 + o\_z^2}\right) \\ \arctan2\left(o\_z, a\_z\right) \end{bmatrix} \tag{A6}$$

Using the sign-aware operator arctan2(*y*, *x*) instead of arctan(*y*/*x*) allows angles to be in (−*π*, +*π*], removes ambiguities and provides global differentiability. The second residual corresponding to *β* only differs regarding the additional rotation *δ*. Combining (A4) and (A5) leads to

$$\begin{aligned} ^{D}\mathbf{R}\_{E}(\boldsymbol{\mathfrak{g}}^{\prime},\boldsymbol{\mathfrak{q}}) = \begin{bmatrix} n\_{\boldsymbol{x}}^{\prime} & o\_{\boldsymbol{x}}^{\prime} & a\_{\boldsymbol{x}}^{\prime} \\ n\_{\boldsymbol{y}}^{\prime} & o\_{\boldsymbol{y}}^{\prime} & a\_{\boldsymbol{y}}^{\prime} \\ n\_{\boldsymbol{z}}^{\prime} & o\_{\boldsymbol{z}}^{\prime} & a\_{\boldsymbol{z}}^{\prime} \end{bmatrix} = \begin{bmatrix} \mathbb{C}\_{\delta} \boldsymbol{n}\_{\boldsymbol{x}} - \operatorname{\mathbb{S}}\_{\delta} \boldsymbol{n}\_{\boldsymbol{y}} & \mathbb{C}\_{\delta} \boldsymbol{o}\_{\boldsymbol{x}} - \operatorname{\mathbb{S}}\_{\delta} \boldsymbol{o}\_{\boldsymbol{y}} & \mathbb{C}\_{\delta} \boldsymbol{a}\_{\boldsymbol{x}} - \operatorname{\mathbb{S}}\_{\delta} \boldsymbol{a}\_{\boldsymbol{y}} \\ \mathbb{C}\_{\delta} \boldsymbol{n}\_{\boldsymbol{y}} + \operatorname{\mathbb{S}}\_{\delta} \boldsymbol{n}\_{\boldsymbol{x}} & \mathbb{C}\_{\delta} \boldsymbol{o}\_{\boldsymbol{y}} + \operatorname{\mathbb{S}}\_{\delta} \boldsymbol{o}\_{\boldsymbol{x}} & \mathbb{C}\_{\delta} \boldsymbol{a}\_{\boldsymbol{y}} + \operatorname{\mathbb{S}}\_{\delta} \boldsymbol{a}\_{\boldsymbol{x}} \\ \boldsymbol{n}\_{\boldsymbol{z}} & o\_{\boldsymbol{z}} & \boldsymbol{a}\_{\boldsymbol{z}} \end{bmatrix} \end{aligned}$$

where *C<sup>δ</sup>* = cos(*δ*) and *S<sup>δ</sup>* = sin(*δ*). The *ZYX* Euler angles from this rotation matrix are

$$\mathbf{a}' = \begin{bmatrix} a'\_1 \\ a'\_2 \\ a'\_3 \end{bmatrix} = \begin{bmatrix} \arctan2\left(n'\_{y'}, n'\_x\right) \\ \arctan2\left(-n'\_{z'}\sqrt{a'^2\_z + o^2\_z}\right) \\ \arctan2\left(o'\_z, a'\_z\right) \end{bmatrix} = \begin{bmatrix} \arctan2\left((\mathbb{C}\_\delta n\_y + \mathbb{S}\_\delta n\_x), (\mathbb{C}\_\delta n\_x - \mathbb{S}\_\delta n\_y)\right) \\ \arctan2\left(-n\_{z,\prime}\sqrt{a^2\_z + o^2\_z}\right) \\ \arctan2\left(o\_z, a\_z\right) \end{bmatrix},\tag{A7}$$

where *δ* only influences the first component *α* <sup>1</sup>. This allows the conclusion that *β*<sup>3</sup> only influences *α*<sup>1</sup> and results in the dependencies

$$
\alpha\_1' = \alpha\_1'(\mathfrak{q}, \mathfrak{k}\_1, \mathfrak{k}\_2, \mathfrak{k}\_3) \tag{A8}
$$

$$
\pi\_2' = \pi\_2'(\mathfrak{q}, \beta\_1, \beta\_2) = \mathfrak{a}\_2 \tag{A9}
$$

$$
\mathfrak{a}'\_3 = \mathfrak{a}'\_3(\mathfrak{q}, \mathfrak{z}\_1, \mathfrak{z}\_2) = \mathfrak{a}\_3 \tag{A10}
$$

with the consequences for the kinematic modeling of robots in 3T2R tasks described in Section 3.

### *Appendix A.2. Relation of the Geometric Jacobian and the Partial Derivative of the Rotation Matrix*

To be able to calculate the gradient matrices of Section 5, the very uncommon partial derivative of the end-effector rotation matrix <sup>0</sup>*R<sup>E</sup>* with respect to the joint coordinates *<sup>q</sup>*<sup>1</sup> of the kinematic chain is needed as term "III" in (34) and (38). It can either be derived with computer algebra systems from the analytic expression of the rotation matrix, or by first devising the time derivative

$${}^{0}\dot{\mathbb{R}}\_{E}(q\_{1},\dot{q}\_{1}) = \frac{\mathbf{d}}{\mathbf{d}t} {}^{0}\mathbb{R}\_{E}(q\_{1}) = \frac{\partial}{\partial q\_{1}} {}^{0}\mathbb{R}\_{E}(q\_{1})\dot{q}\_{1} \tag{A11}$$

and then performing the derivative w.r.t *q*˙ <sup>1</sup>. Using the operators of Appendix A.4 on the relation between rotation matrices and angular velocities leads to the formulation

$$\frac{\partial}{\partial q\_1} \, ^0 \overline{\mathcal{R}}\_E(q\_1) = \frac{\partial}{\partial \dot{q}\_1} \, ^0 \dot{\overline{\mathcal{R}}}\_E = \frac{\partial}{\partial \dot{q}\_1} \prod \left( \overline{\mathcal{S}}(\omega\_E), \, ^0 \overline{\mathcal{R}}\_E \right) \quad \text{with} \quad \omega\_E = f\_{\omega, 1}(q\_1) \dot{q}\_{1'} \tag{A12}$$

where the end-effector angular velocity *ω<sup>E</sup>* can be expressed with the rotational part *Jω*,1 of the geometric Jacobian. It is used as the argument of the cross product matrix operator in stacked form

$$\mathbf{S}(\omega) = \begin{pmatrix} 0 & \omega\_z & -\omega\_y & -\omega\_z & 0 & \omega\_x & \omega\_y & -\omega\_x & 0 \end{pmatrix}^T \quad \text{with} \quad \omega = \begin{pmatrix} \omega\_x & \omega\_y & \omega\_z \end{pmatrix}^T. \tag{A13}$$

Applying again the chain rule for differentiation of the function composition of the three functions ∏ (matrix product), *S* (cross product matrix) and *ω<sup>E</sup>* (linear relation of velocities), as practiced in Section 5, yields

$$\frac{\partial}{\partial q\_1} \, ^0 \overline{\mathbb{R}}\_E = \underbrace{\frac{\partial}{\partial \overline{\mathbb{S}}} \prod}\_{\mathbf{I} \in \mathbb{R}^{9 \times 9}} \underbrace{\left(\overline{\mathbb{S}} \, \prescript{0}{}{\mathbf{R}}\_E\right)}\_{\mathbf{I} \in \mathbb{R}^{9 \times 9}} \underbrace{\frac{\partial}{\partial \omega} \overline{\mathbb{S}}(\omega)}\_{\mathbf{II} \in \mathbb{R}^{9 \times 3}} \underbrace{J\_{\omega, 1}(q\_1)}\_{\mathbf{III} \in \mathbb{R}^{3 \times \text{dim}(q\_1)}} \in \mathbb{R}^{9 \times \text{dim}(q\_1)}.\tag{A14}$$

The first factor is (A28), where the contents of <sup>0</sup>*R<sup>E</sup>* are inserted. The second factor is a sign matrix which can be derived directly from (A13) and only contains 0/−1/1. The third factor is the rotational part of the geometric Jacobian. This completes the calculation of the gradient matrix (34) for the Euler angles orientation residual, which was first introduced, but not pursued further in [7].

### *Appendix A.3. Relation of the Inverse- and Direct-Kinematics Matrices to the Analytic Jacobian*

The properties of the gradient matrices from Section 5 can be illustrated further by comparing them to the analytic Jacobian of one kinematic chain, which gives the platform or end-effector velocity

$$
\dot{\mathbf{x}} = J\_1 \dot{q}\_{1'} \quad \begin{bmatrix} \dot{\mathbf{x}}\_{\mathbf{t}} \\ \dot{\mathbf{x}}\_{\mathbf{r}} \end{bmatrix} = \begin{bmatrix} J\_{\mathbf{t},1}(q\_1) \\ J\_{\mathbf{r},1}(q\_1) \end{bmatrix} \dot{q}\_1. \tag{A15}
$$

Also using the differential form (29) for the first kinematic leg chain of the parallel robot gives

$$\frac{d}{dt}\Phi\_1(q\_1,\mathbf{x}) = \Phi\_{1,\bar{\partial}q\_1}(q\_1,\mathbf{x})\dot{q}\_1 + \Phi\_{1,\bar{\partial}\mathbf{x}}(q\_1,\mathbf{x})\dot{\mathbf{x}} = \mathbf{0} \tag{A16}$$

and results reorganized to the form of (A15) with (32) and (35) written component-wise to

$$
\begin{bmatrix}
\dot{\boldsymbol{x}}\_{\rm t} \\
\dot{\boldsymbol{x}}\_{\rm t}
\end{bmatrix} = -\begin{bmatrix}
\mathbf{0} & \boldsymbol{\Phi}\_{\rm t,1,\partial\mathbf{x}\_{\rm t}}^{-1}
\end{bmatrix} \begin{bmatrix}
\boldsymbol{\Phi}\_{\rm t,1,\partial q\_{1}} \\
\boldsymbol{\Phi}\_{\rm t,1,\partial q\_{1}}
\end{bmatrix} \dot{q}\_{1} = \begin{bmatrix}
\boldsymbol{\Phi}\_{\rm t,1,\partial q\_{1}} \\
\end{bmatrix} \dot{q}\_{1}.
\tag{A17}
$$

By equating coefficients of (A15) and (A17) the relations

$$
\Phi\_{1,\exists q\_1}(q\_1) = -\Phi\_{1,\exists \mathbf{x}}(q\_1, \mathbf{x}) f\_1(q\_1),\tag{A18}
$$

$$\Phi\_{\mathbf{t},1,\bar{\partial}q\_1}(q\_1) = f\_{\mathbf{t},1}(q\_1) \quad \text{and} \tag{A19}$$

$$
\Phi\_{\mathbf{r},1,\partial q\_1}(q\_{1'},\mathbf{x}) = -\Phi\_{\mathbf{r},1,\partial \mathbf{x}\_{\mathbf{r}}}(q\_{1'},\mathbf{x}) f\_{\mathbf{r},1}(q\_1) \tag{A20}
$$

can be obtained between the gradient matrices and the analytic Jacobian *J*<sup>1</sup> of the serial leg chain. The dependency on *q*<sup>1</sup> and *x* has been added to highlight the main requirement, namely the zero equality condition of (A16): (A18) and (A20) only hold if the orientation residual is zero. In the inverse kinematics procedure of Section 5.3 the residual at step *k* in (40) is unequal to zero, which means that the Equation (A20) does not hold in this case. The translational part is unaffected, as can be seen by the missing argument *x* in (A19).

### *Appendix A.4. Matrix Operations for Partial Derivatives*

To simplify the calculations of the gradient matrices of the residuals in Section 5, operators for matrices are replaced by operators for vectors, to avoid differentiating matrices or w.r.t. matrices which would require multi-dimensional tensors. The column operator *R* for rotation matrices *R* to stack the coordinate systems unit vectors *<sup>n</sup>*, *<sup>o</sup>*, *<sup>a</sup>* <sup>∈</sup> <sup>R</sup><sup>3</sup> vertically instead of horizontally is defined as

$$\overline{\mathcal{R}}(\mathcal{R}) = \begin{bmatrix} n \\ o \\ a \end{bmatrix} \in \mathbb{R}^9 \quad \text{with} \quad \mathcal{R} = \begin{bmatrix} n & o & a \end{bmatrix} = \begin{bmatrix} n\_x & o\_x & a\_x \\ n\_y & o\_y & a\_y \\ n\_z & o\_z & a\_z \end{bmatrix} \in \text{SO}(3) \tag{A21}$$

to avoid differentiating matrices or w.r.t. matrices. The special properties of the SO(3) group are not exploited and the operator can be used for R3×<sup>3</sup> as well. Matrix multiplication is expressed with the matrix product operator Π such that

$$\mathbf{R}^1 \overline{\mathbf{R}}\_3 = \overline{\prod} \left( {}^1 \overline{\mathbf{R}}\_2 \, {}^2 \overline{\mathbf{R}}\_3 \right) = \overline{\mathbf{R}} ({}^1 \mathbf{R}\_3) \quad \text{with} \quad {}^1 \mathbf{R}\_3 = {}^1 \mathbf{R}\_2 {}^2 \mathbf{R}\_3. \tag{A22}$$

The transposition operator *P*<sup>T</sup> is a 9 × 9 permutation matrix such that

$$\mathbb{P}^2 \overline{\mathbb{R}}\_1 = \mathbb{P}\_{\mathbb{T}}^1 \overline{\mathbb{R}}\_2 = \overline{\mathbb{R}}(^1 \mathbb{R}\_2^{\mathbb{T}}) = {}^1 \overline{\mathbb{R}}\_2^{\mathbb{T}} \in \mathbb{R}^9 \text{ with } \,^2 \mathbb{R}\_1 = {}^1 \mathbb{R}\_2^{\mathbb{T}} \in \text{SO}(3) \text{ and } {}^1 \overline{\mathbb{R}}\_2 = \overline{\mathbb{R}}(^1 \mathbb{R}\_2). \tag{A23}$$

Writing <sup>1</sup>*R*<sup>T</sup> <sup>2</sup> instead of *P*<sup>T</sup> <sup>1</sup>*R*<sup>2</sup> serves for the clarity of the expressions (34,36,38,39) and overloads the transposition operator for R<sup>9</sup> noted with the bar.

### *Appendix A.5. Contents of the Partial Derivatives*

The single expressions derived in Section 5 can be calculated with low computational effort from the definition of the *XYZ* and *ZYX* Euler angles from (5), (21) and (A6). With *R* = [*nx*, *ny*, *nz*, *ox*, *oy*, *oz*, *ax*, *ay*, *az*] <sup>T</sup> the gradient "I" in (34,36,38,39) for *ZYX* angles becomes

$$\frac{\partial \mathfrak{a}(\overline{\mathbf{R}})}{\partial \overline{\mathbf{R}}} = \begin{bmatrix} -\frac{n\_y}{n\_x^2 + n\_y^2} & \frac{n\_z}{n\_x^2 + n\_y^2} & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & -\sqrt{a\_z^2 + a\_z^2} & 0 & 0 & \frac{n\_z a\_z}{\sqrt{a\_z^2 + a\_z^2}} & 0 & 0 & \frac{n\_z a\_z}{\sqrt{a\_z^2 + a\_z^2}}\\ 0 & 0 & 0 & 0 & 0 & \frac{a\_z}{a\_z^2 + a\_z^2} & 0 & 0 & -\frac{a\_z}{a\_z^2 + a\_z^2} \end{bmatrix} \tag{A24}$$

and the reciprocal gradient "IV" in (36) for *XYZ* angles yields

*∂R*(*β*) *<sup>∂</sup><sup>β</sup>* <sup>=</sup> ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 0 −*S*<sup>2</sup> *C*<sup>3</sup> −*C*<sup>2</sup> *S*<sup>3</sup> *C*<sup>1</sup> *S*<sup>2</sup> *C*<sup>3</sup> − *S*<sup>1</sup> *S*<sup>3</sup> *S*<sup>1</sup> *C*<sup>2</sup> *C*<sup>3</sup> −*S*<sup>1</sup> *S*<sup>2</sup> *S*<sup>3</sup> + *C*<sup>1</sup> *C*<sup>3</sup> *S*<sup>1</sup> *S*<sup>2</sup> *C*<sup>3</sup> + *C*<sup>1</sup> *S*<sup>3</sup> −*C*<sup>1</sup> *C*<sup>2</sup> *C*<sup>3</sup> *C*<sup>1</sup> *S*<sup>2</sup> *S*<sup>3</sup> + *S*<sup>1</sup> *C*<sup>3</sup> 0 *S*<sup>2</sup> *S*<sup>3</sup> −*C*<sup>2</sup> *C*<sup>3</sup> −*C*<sup>1</sup> *S*<sup>2</sup> *S*<sup>3</sup> − *S*<sup>1</sup> *C*<sup>3</sup> −*S*<sup>1</sup> *C*<sup>2</sup> *S*<sup>3</sup> −*S*<sup>1</sup> *S*<sup>2</sup> *C*<sup>3</sup> − *C*<sup>1</sup> *S*<sup>3</sup> −*S*<sup>1</sup> *S*<sup>2</sup> *S*<sup>3</sup> + *C*<sup>1</sup> *C*<sup>3</sup> *C*<sup>1</sup> *C*<sup>2</sup> *S*<sup>3</sup> *C*<sup>1</sup> *S*<sup>2</sup> *C*<sup>3</sup> − *S*<sup>1</sup> *S*<sup>3</sup> 0 *C*<sup>2</sup> 0 −*C*<sup>1</sup> *C*<sup>2</sup> *S*<sup>1</sup> *S*<sup>2</sup> 0 −*S*<sup>1</sup> *C*<sup>2</sup> −*C*<sup>1</sup> *S*<sup>2</sup> 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (A25)

with *Ci* = cos(*βi*), *Si* = sin(*βi*). The property

$$\mathbb{P}\left(\frac{\partial \mathcal{B}}{\partial \overline{\mathcal{R}}}\right)\left(\frac{\partial \overline{\mathcal{R}}(\mathcal{B})}{\partial \mathcal{B}}\right) = \mathbf{1} \in \mathbb{R}^{3 \times 3} \tag{A26}$$

can be used to test the implementation, if (A24, A25) are defined for the same Eulerangle convention.

The gradient of the matrix product (A22) w.r.t. the second factor used in (34)/II, (36)/III, (38)/II and (39)/III is

$$\frac{\partial}{\partial \overline{\mathbf{R}}\_2} \prod \left( \overline{\mathbf{R}}\_1, \overline{\mathbf{R}}\_2 \right) = \begin{bmatrix} \mathbf{R}\_1 & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{R}\_1 & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{R}\_1 \end{bmatrix} \tag{A27}$$

and to complete the enumeration the gradient w.r.t. the first factor used in (A14)/I is

$$\frac{\partial}{\partial \overline{\mathcal{R}}\_1} \prod (\overline{\mathcal{R}}\_1, \overline{\mathcal{R}}\_2) = \begin{bmatrix} \text{diag}(n\_x) & \text{diag}(o\_x) & \text{diag}(a\_x) \\ \text{diag}(n\_y) & \text{diag}(o\_y) & \text{diag}(a\_y) \\ \text{diag}(n\_z) & \text{diag}(o\_z) & \text{diag}(a\_z) \end{bmatrix}, \tag{A28}$$

where *nx*, *ny*, ... are the entries of *R*<sup>2</sup> and the diagmatrices are 3 × 3. By transposing the elements of the matrix product (A22), only the first form (A27) had to be used in Section 5.

### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
