Joint-Space Characterization of a Medical Parallel Robot Based on a Dual Quaternion Representation of SE(3)

Abstract

The paper proposes a mathematical method for redefining motion parameterizations based on the joint-space representation of parallel robots. The study parameters of SE(3) are used to describe the robot kinematic chains, but, rather than directly analyzing the mobile platform motion, the joint-space of the mechanism is studied by eliminating the Study parameters. From the loop equations of the joint-space characterization, new parameterizations are defined, which enable the placement of a mobile frame on any mechanical element within the parallel robot. A case study is presented for a medical parallel robotic system in which the joint-space characterization is achieved and based on a new defined parameterization, the kinematics for displacement, velocities, and accelerations are studied. A numerical simulation is presented for the derived kinematic models, showing how the medical robot guides the medical tool (ultrasound probe) on an imposed trajectory.

1. Introduction

One fundamental problem in robotics (relevant for both design and control) refers to mechanism analysis and synthesis, which, in general, provides mathematical information regarding: (i) the laws of motion of the mechanical system (through forward and inverse kinematics); (ii) the singular configurations of the robot (in which the control of the robot may be lost); and, (iii) the robot workspace. There are various mathematical methods used for mechanism analysis and synthesis, with two categories being prominent, namely the vector methods and the algebraic methods. Examples for vector methods are: the Denavit-Hartenberg convention (based on homogeneous matrix representation of SE(3)), used especially for serial mechanical chains [1], with a novel algorithm for automatic identification proposed in [2]; screw theory, which may be used even for error analysis in parallel mechanisms [3] or calibration of parallel mechanism [4]; and, Euler parameters used in the study of spherical displacement [5,6]. An overview of mathematical models used in mechanism analysis for higher order kinematics is given in [7]. One specific algebraic method is an extension of the Euler parameters (which are the unit quaternion representation of three-dimensional (3D) rotations) and it is based on the study parameters (or dual quaternion representation) of SE(3) [8]. In parallel mechanisms analysis, the study parameters method has the advantage that it describes the global kinematics. Consequently, the method reveals the whole mechanism workspace, all of the working and assembly modes of the mechanism, and all of the singularities. Moreover, the study parameters of SE(3) has a close connection with screw theory (this fact is pointed in [9,10]). In the past, the study parameters method was used on medical robots (see [11,12,13]) for complete singularity analysis, which, in turn, provided valuable information for the safe operation of those medical robots (singularity avoidance is imperative in operating a medical robot to avoid injury). An algebraic formulation for collision avoidance using the dual quaternion representation of SE(3) was shown in [14], which might be useful in many areas (e.g., maximizing safety in medical robotics).

A new parameterization method was introduced in [15], which is also based on the study parameters, but rather than describing the relative motion between a mobile platform (or end effector) and a fixed one, the method eliminates the mobile platform motion information and describes the parallel robot joint-space. In this work, the term “joint-space” refers to the mathematical space of all the mechanical joints of a parallel robot (in contrast to e.g., [16], where the authors use the term for the active joints), since one has the option to choose which joints are actuated. The scope of the new parameterization introduced in [15] was to analyze medical rehabilitation robots where specific mechanical revolute joints had one-to-one correlation with corresponding anatomic joints (similarly to an exoskeleton), while the actual pose of the end-effector had no relevance for the application [17,18,19,20].

In this work, the authors further explore the usability of the method presented in [15], by proposing a novel approach to define new parameterizations that are based on the constraint equations that describe the joint-space. The main goal is to show how the relative motion between a fixed coordinate frame and a mobile one (where both coordinate frames may be attached at any element within the parallel robot) can be redefined. The method is used on a medical robot for minimally invasive liver cancer treatment where the motion of a mobile coordinate frame (attached to the medical tool that the robot guides) is required to derive the robot kinematics (displacement, velocity, and acceleration).

Following the Introduction Section, the paper is structured, as follows: Section 2 presents the mathematical background needed to describe the proposed method. The parameterization approach is generally described; the study parameters are briefly presented showing the similarities between them and the method proposed in this work; the joint-space characterization method and the technique to redefine parameterizations in SE(3) are also generally presented. Furthermore, Section 2 presents an example of joint-space characterization based on a planar closed loop mechanism. In Section 3, a case study is presented for the PROHEP-LCT parallel medical robot where the joint-space of the robot is studied, leading to the robot kinematics for displacement, velocity, and accelerations. Section 4 presents the numerical results for the PRO-HEP-LCT parallel robot, Section 5 presents a discussion, and finally the conclusions are drawn in Section 6.

2. Materials and Methods

The main goal of the paper is to show how a specific motion parameterization method (previously introduced in [15]), which describes the joint-space of a parallel robot, may be used in achieving other parameterizations (that may be defined based on technical needs, such as choosing which joints to actuate) and how the newly defined parameterizations may be used in kinematic studies. The proposed method shares similarities with the algebraic method that is based on the Study parameters of SE(3), namely the fact that the same steps are followed (see Figure 1 steps 0–2).

As illustrated in Figure 1, for mechanism analysis (step 0) using the study parameters method, the kinematic chains are decoupled (step 1) usually at the origin of the mobile platform, and the study parameters are defined for each of the n-th kinematic chain independently (steps 1.1–1.n), thus achieving the kinematic mapping (a concept that will be discussed further in the section) of the kinematic chains (step 2).

In order to describe the relative displacement between a fixed coordinate frame (attached to the parallel robot base) and a mobile one (attached to the robot mobile platform), using the study parameters method, the passive joint parameters (due to passive joints within the parallel robot) are eliminated from each kinematic chain (step 3) and, then, the constraints of the parallel robot are obtained (step 4), since all the study parameter equations of all the kinematic chains must be fulfilled (which physically represents recoupling the kinematic chains to form the parallel robot). Based on the constraint equations (which are polynomial functions of Study parameters), the forward and inverse kinematics (for displacement) may be achieved (step 5.3), as well as the singularity (using e.g., a Jacobian method [21] or a discriminant method [22]) and the workspace analysis (step 5.2 and 5.1, respectively).

If the motion of the mobile platform is not needed, the study parameters of the kinematic chains may be eliminated (step 6) by subtracting sets of study parameter equations, thus obtaining loop equations for the mechanism (step 7). Consequently, the information that describes the relative displacement between the fixed and the mobile coordinate frames is lost. However, the proposed method aims to redefine the motion of a mobile platform, based on the loop equations obtained at step 7. The unwanted motion parameters may be eliminated (although possible this step is not mandatory) from the parallel mechanism (step 8) and explicit solutions are obtained for chosen parameters, which are used to redefine other motion parameterizations and reintroduce the information that describes the relative displacement between the fixed and mobile frames (step 9). The newly defined parameterizations may be used for workspace analysis (step 10.1), kinematics (step 10.2), and singularity analysis (step 10.3).

2.1. The Study Parameters of SE(3)

As pointed out before, the study parameters of SE(3) are not a novel concept and are extensively documented in the scientific literature (see [6,8] for an introduction and [11,12,13] for their applications in kinematics and singularity analysis for medical robots). A brief description of the study parameters is presented in this section for a better understanding and connection with the proposed method.

The isomorphism between the Study parameters and the special Euclidean group SE(3) may be expressed as following:

Let k,

$$\begin{array}{c}k:SE(3)\to P\in P^7\\ D(x_i,y_i)\to [x_0:x_1:x_2:x_3:y_0:y_1:y_2:y_3]\ne [0:0:0:0:0:0:0:0],\end{array}$$

(1)

be a mapping (also called kinematic mapping [6,8]) that maps every Euclidean displacement D onto a point P of the projective space P⁷ (the kinematic image space [6,8]). The coordinates of the point P are called the study parameters of SE(3) and they are based on a dual quaternion where the following relations must be fulfilled:

$$x_0^2+x_1^2+x_2^2+x_3^2\ne 0,$$

(2)

$$x_0{y}_0+x_1{y}_1+x_2{y}_2+x_3{y}_3=0.$$

(3)

Equation (2) represents a normalizing condition, whereas Equation (3) represents the study quadric [6,8]. Assuming a homogeneous 4 × 4 transformation matrix M ∈ SE(3) with the form:

$$M:=\left[\begin{array}{cc}1& 0\\ t& R\end{array}\right],$$

(4)

were R is a 3 × 3 rotation matrix, t[t_x, t_y, t_z]^T is a translation vector, and 0[0,0,0] is a null vector, the homogeneous transformation matrix may be computed from the study parameters, as following:

$$M:=\left[\begin{array}{cccc}1& 0& 0& 0\\ 2(x_0{y}_1-x_1{y}_0+x_2{y}_3-x_3{y}_2)& x_0^2+x_1^2-x_2^2-x_3^2& 2(-x_0{x}_3+x_1{x}_2)& 2(x_0{x}_2+x_1{x}_3)\\ 2(x_0{y}_2-x_1{y}_3-x_2{y}_0+x_3{y}_1)& 2(x_0{x}_3+x_1{x}_2)& x_0^2-x_1^2+x_2^2-x_3^2& 2(-x_0{x}_1+x_2{x}_3)\\ 2(x_0{y}_3+x_1{y}_2-x_2{y}_1-x_3{y}_0)& 2(-x_0{x}_2+x_1{x}_3)& 2(x_0{x}_1+x_2{x}_3)& x_0^2-x_1^2-x_2^2+x_3^2\end{array}\right],$$

(5)

and the study parameters may be computed from a homogeneous transformation matrix as following (where a_ij represent the entries of the 4 × 4 homogeneous matrix M):

$$\begin{array}{cc}{x}_0:x_1:x_2:x_3=& 1+a_{11}+a_{22}+a_{33}:a_{32}-a_{23}:a_{13}-a_{31}:a_{21}-a_{12}\\ & a_{32}-a_{23}:1+a_{11}-a_{22}-a_{33}:a_{12}-a_{21}:a_{31}+a_{13}\\ & a_{13}-a_{31}:a_{12}+a_{21}:1-a_{11}+a_{22}-a_{33}:a_{23}+a_{32}\\ & a_{21}-a_{12}:a_{31}+a_{13}:a_{23}+a_{32}:1-a_{11}-a_{22}+a_{33}\end{array},$$

(6)

$$\begin{gathered} y_0=-\frac{1}{2}(x_1{t}_x+x_2{t}_y+x_3{t}_z),\hspace{1em}{y}_1=-\frac{1}{2}(-x_0{t}_x+x_2{t}_z-x_3{t}_y), \\ {y}_2=-\frac{1}{2}(-x_0{t}_y-x_1{t}_z+x_3{t}_x),\hspace{1em}{y}_3=-\frac{1}{2}(-x_0{t}_z+x_1{t}_y-x_2{t}_x), \end{gathered}$$

(7)

Observation 1. 

the Study parameters x₀:x₁:x₂:x₃:y₀:y₁:y₂:y₃ encode both orientation information through the x_i parameters (being the Euler parameters) and translation information through y_i parameters.

As pointed out before, one advantage of the Study parameters method is that they describe the global kinematics through constraint varieties that formed over the study parameter space. The method describes the entire workspace of a parallel mechanism, with all of the working modes and assembly modes, which, in turn, shows all of the theoretical mechanism singularities (of course not all singularities may be physically possible). Moreover, if the substitution of the tangent of half angle is used for the trigonometric functions, the constraint equations become purely polynomial, since the denominator introduced by the substitution may be factored out (the study parameters are homogeneous). This enables the use of computer algebra software for various tasks, such as: (i) eliminating unwanted parameters (usually the motion parameters describing free motion joints in parallel robots); (ii) computing the “fundamental” polynomials describing the robot workspace in the form of reduced Groebner bases. The parameterization singularity introduced by the tangent of half angle substitution (the rotation angle of π) may be avoided by defining the parameterization, such that the rotation angle of π is not in the operational workspace or it is equivalent to a known singularity of the mechanism. For an in-depth introduction into algebraic geometry (with great focus on polynomial ideals of affine and projective varieties), and computational methods to work with these polynomial ideals (e.g., Groebner bases) and their properties, the reader is referred to [23].

2.2. Joint Space Caracterization

In some robotic applications, the relative motion between the base and the end effector does not present any real interest. One example (see e.g., [19,20,24]) is the motor rehabilitation of the limbs, where some specific parallel robots are designed to guide the limbs, and there exists an almost one-to-one correlation between some mechanical revolute joints and the anatomic joints (since geometrically speaking the revolute joints share the rotational axes with the anatomic ones).

When considering the mapping defined in Equation (1), the specific mapping for each kinematic chain of the mechanism can be defined:

Let

$$D_j(x_i,y_i)\to [x_0:x_1:x_2:x_3:y_0:y_1:y_2:y_3]={\rho }_j[r_{j,1}:r_{j,2}:r_{j,3}:r_{j,4}:t_{j,1}:t_{j,2}:t_{j,3}:t_{j,4}]$$

(8a)

be the kinematic mapping of the j th kinematic chain subject to the geometrical constraints of the mechanism, were:

$$\begin{array}{cc}\hfill r_{j,k}({\alpha }_{j,1},{\alpha }_{j,2},\dots ,{\alpha }_{j,n})-& \text{the multivariate polynomial with} n \text{rotation motion parameters};\hfill \\ \hfill t_{j,k}({\alpha }_{j,1},{\alpha }_{j,2},\dots ,{\alpha }_{j,n},p_{j,1},p_{j,2},\dots ,p_{j,k})-& \text{the multivariate polynomial with} n \text{rotation parameters and k translation motion parameters};\hfill \\ \hfill {\alpha }_{j,k}=\tan \left(\frac{{\theta }_{j,k}}{2}\right),{\theta }_{j,k}-& \text{the rotation angle on the} nth \text{revolute joint of the} jth \text{kinematic chain} \hfill \\ \hfill p_{j,k}-& \text{the} kth \text{translation parameter (fixed of variable) of the} jth \text{kinematic chain};\hfill \\ \hfill {\rho }_j-& \text{the homogenizing parameter of the} jth \text{kinematic chain}.\hfill \end{array}$$

(8b)

The r_j,k and t_j,k multivariate polynomials are consistent with Observation 1. The p_j,k parameter might appear due to different components in the kinematic chain; it can be a geometric parameter (e.g., mechanical links lengths), an active translation (due to a linear actuator), or a passive translation (due to a free motion linear joint).

Choosing two kinematic chains and subtracting their kinematic images D_j eliminates the study parameters from the sets of polynomials from the mappings. Consequently, the information that describes the relative displacement between the fixed and mobile platform is lost, and the constraints are described in the joint space of the mechanism. The subtraction of two kinematic chains in general form is:

$$L:=D_{j=a}-D_{j\ne a},$$

(9)

where L contains a maximum of eight multivariate polynomials (which must simultaneously vanish). There are two special cases in this formulation: a) the pure spherical displacement which is described by 4 polynomials (y_i = 0); b) the planar motion which is also described by four polynomials (e.g., the motion in OXY plane is obtained by setting x₁ = 0, x₂ = 0, y₁ = 0, y₂ = 0). The homogenizing parameters ρ_j are not allowed to vanish to avoid the trivial solution.

Observation 2. 

the joint space of the parallel mechanism will be invariant to the coordinate frames placements. As long as the loop is closed when parameterizing the parallel mechanism, it does not matter whether the coordinate systems were placed initially.

The multivariable constraint polynomials resulted within L form an ideal of variety over the joint-space of the parallel mechanism. The model so far describes an over-constrained mechanism; consequently, some parameters must be eliminated (by considering them as free motion parameters). The elimination can be reliably achieved using algebraic methods, such as computing elimination ideals using Groebner bases.

There are two main uses of the loop equations in algebraic form obtained after the elimination of unwanted parameters:

(1)

The description of the joint space where certain dependencies are obtained: the method is suitable for applications were these dependencies are more important than the motion of a mobile platform (e.g., rehabilitation using exoskeletons based on a parallel mechanism). Furthermore, the constraints may be used in velocity and acceleration kinematics and also reveal the singularities of the mechanism.

(2)

The loop constraints may be used in defining new parameterizations e.g., assigning coordinate systems at any element in the mechanism (which is presented on a case study in Section 3).

2.3. Defining New Motion Parameterizations Based on the Joint-Space Characterization

Equation (9) describes the joint-space of two coupled kinematic chains in a general form, the unwanted motion parameters may be eliminated. When the mechanism contains more than two kinematic chains, more sets of equations are needed to describe the mechanism joint-space.

To define new parameterizations based on the loop equations (which describe the joint-space), the authors propose the following stepwise process:

• define the input motion parameters (the actuated joints denoted q_i); at this point the passive joints and the active ones should be established;
• define the fixed coordinate frame and the mobile one, respectively; the fixed coordinate frame XYZ is defined as being attached to a chosen base and the mobile coordinate frame X’Y’Z’ is attached to the chosen mobile platform;
• define a mechanical “path” (denoted P, composed of joints and links), which connects the fixed frame with the mobile one; usually the shortest path yields the simpler equations, however it is not mandatory to define the shortest path (as it is illustrated in the example shown in Section 2.3);
• determine the motion parameters (denoted a_i) within P;
• eliminate the unwanted motion parameters (denoted b_i); the parameters b_i are all the motion parameters that are not q_i or a_i;
• express the a_i parameters with respect to the input parameters q_i (f_i represents an irrational function of the active joints q_i):

  $$\begin{array}{c}{a}_1=f_1(q_1..q_i)\\ ⋮\\ a_i=f_i(q_1..q_i)\end{array}$$

  (10)
• create a geometric parameterization of P; one example is to use the Denavit-Hartenberg convention; and,
• substitute a_i into P; at this point, the relative displacement between the fixed coordinate frame XYZ and the mobile one X’Y’Z’ (with respect to the input parameters q_i) is achieved. Consequently, the constraints of the mechanism are obtained vis-à-vis the newly defined parameterization.

2.4. Joint Space Characterization of a 6-R Closed Loop Planar Parallel Mechanism

In this section, the joint space of a classic 6-R parallel (closed loop) mechanism is modeled based on the method previously discussed in Section 2. Figure 2 illustrates the 6-R planar parallel mechanism. The mechanism must have three input (revolute) joints (or active joints) to represent a stable mechanical system. With less than three inputs, the mechanism will have self-motion, whereas with more than three inputs the mechanism will be over-constrained.

When assuming that the coordinate system OXYZ is the fixed one and the link l₆ is fixed (along the direction of the OY axis) and the parallel mechanism is decoupled at the origin of the O’X’Y’Z’ coordinate system with its origin at the middle of the link l₃ and the O’Y’ axis along the link l₃, two kinematic chains are defined with the following mappings:

$$M_1:=Q_1({\theta }_1)\cdot Q_2(l_1)\cdot Q_3({\theta }_2)\cdot Q_4(l_2)\cdot Q_5({\theta }_3)\cdot Q_6(l_3/2),$$

(11a)

$$M_2=Q_7(l_6)\cdot Q_8({\theta }_6)\cdot Q_9(l_5)\cdot Q_{10}({\theta }_5)\cdot Q_{11}(l_4)\cdot Q_{12}({\theta }_4)\cdot Q_{13}(-l_3/2),$$

(11b)

where the study parameters used in Equations (11a) and (11b) are detailed in

Table 1. The Study parameters used in the motion parameterization of the 6R mechanism.

Symbol	Explicit form	Details
Q1(θ1)	[1:0:0:α1:0:0:0:0]	θ1 = π/2 → α1 = 0 1
Q2(l1)	[1:0:0:0:0:−l1/2:0:0]	translation along OX
Q3(θ2)	[1:0:0:α2:0:0:0:0]	θ2 = π → α2 = 0 1
Q4(l2)	[1:0:0:0:0:−l2/2:0:0]	translation along OX
Q5(θ3)	[1:0:0:α3:0:0:0:0]	θ3 = π/2 → α3 = 0 1
Q6(l3/2)	[1:0:0:0:0:0:−l3/4:0]	translation along OY
Q7(l6)	[1:0:0:0:0:0:−l6/2:0]	translation along OY
Q8(θ6)	[1:0:0:α6:0:0:0:0]	θ6 = π/2 → α6 = 0 1
Q9(l5)	[1:0:0:0:0:−l5/2:0:0]	translation along OX
Q10(θ5)	[1:0:0:α5:0:0:0:0]	θ5 = π → α5 = 0 1
Q11(l4)	[1:0:0:0:0:−l4/2:0:0]	translation along OX
Q12(θ4)	[1:0:0:α4:0:0:0:0]	θ4 = π/2 → α4 = 0 1
Q13(-l3/2)	[1:0:0:0:0:0:l3/4:0]	translation along OY

¹α_i are the tangents of the half angles of θ_i.

.

Numerical values are given for the links (l_i = 1 i = 1..5, l₆ = 2) in order to facilitate the computation. After multiplying the study parameters on the right hand side of Equations (11a) and (11b) M₁ and M₂ read (also considering the homogenizing parameters ρ₁, ρ₂):

$$M_1:=\left[\begin{array}{c}{x}_0\\ x_1\\ x_2\\ x_3\\ y_0\\ y_1\\ y_2\\ y_3\end{array}\right]={\rho }_1\left[\begin{array}{c}4(({\alpha }_3+{\alpha }_2){\alpha }_1+{\alpha }_2{\alpha }_3-1)\\ 0\\ 0\\ 4({\alpha }_3{\alpha }_2{\alpha }_1-{\alpha }_2-{\alpha }_3-{\alpha }_1)\\ 0\\ (({\alpha }_2+4){\alpha }_3-1){\alpha }_1-{\alpha }_2-{\alpha }_3+4\\ (-{\alpha }_2-{\alpha }_3+4){\alpha }_1+1-({\alpha }_2+4){\alpha }_3\\ 0\end{array}\right]$$

(12a)

$$M_2:=\left[\begin{array}{c}{x}_0\\ x_1\\ x_2\\ x_3\\ y_0\\ y_1\\ y_2\\ y_3\end{array}\right]={\rho }_2\left[\begin{array}{c}4(({\alpha }_4+{\alpha }_5){\alpha }_6+{\alpha }_4{\alpha }_5-1)\\ 0\\ 0\\ 4({\alpha }_6{\alpha }_5{\alpha }_4-{\alpha }_6-{\alpha }_5-{\alpha }_4)\\ 0\\ -5((({\alpha }_5+\frac{4}{5}){\alpha }_4-1){\alpha }_6-{\alpha }_5-{\alpha }_4-\frac{4}{5})\\ 3((-{\alpha }_4-{\alpha }_5-\frac{4}{3}){\alpha }_6+1-({\alpha }_5-\frac{4}{3}){\alpha }_5)\\ 0\end{array}\right]$$

(12b)

After subtracting M₂ from M₁ the study parameters are eliminated and a loop is created. It can be checked that, if a different mobile coordinate system is chosen (e.g., X’’Y’’Z’’ in Figure 2), the method yields the same polynomials that describe the joint space of the mechanism chosen in this section as an example (fact pointed out in Observation 2).

Computing L:= M₁ − M₂ results in a set of four polynomials that form an ideal of the variety describing the mechanism in the joint space. Since the 6R closed loop parallel mechanism needs three actuators (inputs) to be mechanically stable, in this example the parameters α₁, α₂, and α₆ are chosen as inputs.

An elimination ideal J₀ (independent of ρ₁, ρ₂, α₅) may be computed, as following:

J₀ = G_L [1..2]; G_L is a Groebner basis, of L with respect to $[{\rho }_1,{\rho }_1,{\alpha }_5]$≺$[{\alpha }_3,{\alpha }_4]$ where:

$$J_0:=\langle g_1^7({\alpha }_1,{\alpha }_2,{\alpha }_3,{\alpha }_4,{\alpha }_6),g_2^8({\alpha }_1,{\alpha }_2,{\alpha }_3,{\alpha }_4,{\alpha }_6)\rangle$$

(13)

J₀ contains two polynomials of total degree 7 and 8, respectively. Moreover, the two polynomials in J₀ contain 66 monomials and 54 monomials, respectively, and they are only presented in their general form. A pure lexicographic monomial ordering Groebner base can be computed to facilitate the solution solving process of the polynomials from J₀:

$G_H-$ is a Groebner basis, of J₀ with respect to α₃ ≺ α₄ (a plex ordering), where:

$$\begin{gathered} G_H:=\langle g_1^8({\alpha }_1,{\alpha }_2,{\alpha }_4,{\alpha }_6),g_1^7({\alpha }_1,{\alpha }_2,{\alpha }_3,{\alpha }_4,{\alpha }_6)\rangle ,\mathrm{with}: \\ {g}_1:=P_1{\alpha }_4^2+P_2{\alpha }_4+P_3;g_2:=P_4{\alpha }_4^{}+P_5{\alpha }_3+P_6; \\ \begin{array}{l}{P}_1:=((3{\alpha }_6^2+8{\alpha }_6+3){\alpha }_2^2+8({\alpha }_6^2+{\alpha }_6+1){\alpha }_2+3{\alpha }_6^2+8{\alpha }_6+11){\alpha }_1^2+(8(-{\alpha }_6^2+1){\alpha }_2-16({\alpha }_6^2+{\alpha }_6+1)){\alpha }_1\\ +(8{\alpha }_6^2+3{\alpha }_6+8){\alpha }_2^2-8({\alpha }_6^2+{\alpha }_6+1){\alpha }_2+11{\alpha }_6^2+8{\alpha }_6+3;\end{array} \\ {P}_2:=-4({\alpha }_1^2+1)({\alpha }_2^2+1)({\alpha }_6^2+1); \\ \begin{array}{l}{P}_3:=((3{\alpha }_6^2+8{\alpha }_6+3){\alpha }_2^2+8({\alpha }_6^2+{\alpha }_6+1){\alpha }_2+3{\alpha }_6^2+8{\alpha }_6+11){\alpha }_1^2+(8(-{\alpha }_6^2+1){\alpha }_2-16({\alpha }_6^2+{\alpha }_6+1)){\alpha }_1\\ +(3{\alpha }_6^2+8{\alpha }_6+3){\alpha }_2^2-8({\alpha }_6^2+{\alpha }_6+1){\alpha }_2+11{\alpha }_6^2+8{\alpha }_6+3;\end{array} \\ \begin{array}{l}{P}_4:=((9{\alpha }_6^2+12{\alpha }_6+9){\alpha }_2^2+8({\alpha }_6^2+{\alpha }_6+2){\alpha }_2+{\alpha }_6^2+4{\alpha }_6+9){\alpha }_1^2+(4(-{\alpha }_6^2+1){\alpha }_2-8(3{\alpha }_6^2+2{\alpha }_6+1)-\\ 12{\alpha }_6^2-16{\alpha }_6-20){\alpha }_1+({\alpha }_6^2+4{\alpha }_6+1){\alpha }_2^2-8({\alpha }_6+1){\alpha }_2+17{\alpha }_6^2+12{\alpha }_6+9;\end{array} \\ \begin{array}{l}{P}_5:=((-3{\alpha }_6^2-8{\alpha }_6-3){\alpha }_2^2-8({\alpha }_6^2+{\alpha }_6+1){\alpha }_2-3{\alpha }_6^2-8{\alpha }_6-11){\alpha }_1^2+(8({\alpha }_6^2-1){\alpha }_2+16({\alpha }_6^2+{\alpha }_6+1)){\alpha }_1\\ -(3{\alpha }_6^2+8{\alpha }_6+3){\alpha }_2^2+8({\alpha }_6^2+{\alpha }_6+1){\alpha }_2-11{\alpha }_6^2-8{\alpha }_6-3;\end{array} \\ \begin{array}{l}{P}_6:=(-8({\alpha }_6^2+{\alpha }_6+1){\alpha }_2^2+4({\alpha }_2^2+1)){\alpha }_1^2+(-8({\alpha }_6^2+{\alpha }_6+1){\alpha }_2^2+8({\alpha }_6^2-1){\alpha }_2+8({\alpha }_6^2+{\alpha }_6+1)){\alpha }_1+\\ 8({\alpha }_6^2+{\alpha }_6+1){\alpha }_2+4{\alpha }_6^2({\alpha }_2^2-1).\end{array} \end{gathered}$$

(14)

g₁ and g₂ are two input-output (I-O) constraint equations of the 6R example shown in this section. It is easy to see from g₁ and g₂ that there are two possible solutions for the rotation parameters α₃ and α₄. A numerical example is given in

Table 2. Solutions for the 6R closed loop mechanism.

I-O	Symbol	Value	Assoc. Angle (°)
I	α1	0.06	θ1 = 6.9
	α2	0.25	θ2 = 28.1
	α6	0.06	θ6 = 6.9
O	α3	0.26; −0.36	θ3 = 29.5; θ3 = −39.7
	α4	5.46; 0.18	θ4 = 159.2; θ4 = 20.8

by evaluating the constraint equations with given inputs (α₁, α₂, and α₃) and solving for the outputs (α₃ and α₄). Furthermore, the Computer Aided Design (CAD) model of the mechanism is shown (Figure 3) for the two distinct solutions (which are part of different working modes).

Based on the Equation (14), a new parameterization may be defined for the motion of the 6R mechanism. As an example, the relative displacement between the XYZ coordinate frame and the X’’Y’’Z’’ one is defined. To simplify the problem, only one solution group is chosen, namely the one that describes the “convex” working mode of the mechanism. The solution groups could be chosen based on technical needs, but, in principle, changing the working mode of the 6R mechanism will require the mechanism to pass through singularity, which in many applications must be avoided.

The solution:

$$S_1:=\{{\alpha }_2=f_1({\alpha }_1,{\alpha }_6);{\alpha }_3=f_2({\alpha }_1,{\alpha }_6)\}$$

(15)

of Equation (14) may be substituted into the kinematic mapping:

$$M_n:=Q_1({\theta }_1)\cdot Q_2(l_1)\cdot Q_3({\theta }_2)\cdot Q_4(l_2)\cdot Q_5({\theta }_3)\cdot Q_6(l_3)\cdot Q_{12}(-l_4)\cdot Q_{11}(-l_5/2),$$

(16)

which describes the kinematic path that connects the fixed frame XYZ with the mobile one X’’Y’’Z’’ in the study parameter representation. Computing the homogeneous matrix from M_n (Equation (16)) yields the Cartesian coordinates of the origin of X’’Y’’Z’’ and its orientation. A numerical example is given for the input values that are shown in Table 2, and the result shows the origin of the mobile platform located at [X = 1.7, Y = 1.9] with an orientation of 44° with respect to the fixed frame (see Figure 3b). When considering the topology of the closed loop 6R mechanism (Figure 3), the input-output parameters α_i will always fulfil Equation (14), indifferent of the placement of the mobile coordinate frame (X’’Y’’Z’’). This is actually a consequence of Observation 2 (Section 2.2).

3. PRO-HEP-LCT Kinematics (a Case Study)

The section presents the mathematical modeling of the PROHEP-LCT [25,26] parallel robotic system kinematics while using the ideas introduced in the previous section. PRO-HEP-LCT contains two identical modules, for guiding both the needles and the ultrasound probe. Special instruments are attached to the mobile platforms of the two modules for: i) inserting needles on parallel paths [27,28] and ii) guiding (inserting and rotating) an intra-operatory ultrasound (US) probe [29]. Figure 4 illustrates one module of the PRO-HEP-LCT parallel robotic system.

Referring to Figure 4, the main mechanical subsystems of the PRO-HEP-LCT parallel robotic system are [26]:

• the upper planar mechanism, which is a PRR-PP mechanism (prismatic-revolute-revolute chain coupled with a prismatic-prismatic chain where the underline represents the actuated joint). The upper planar mechanism creates two DOF pure translational motion (actuated by q₄ and q₅) in a plane parallel to XY plane;
• the lower planar chain, which is a 6R-3R mechanism, i.e., a 6R closed loop planar chain coupled with a 3R chain which creates pure translational 2 DOF motion (actuated by q₂ and q₃) in the XY plane; in the parameterization α₁ and α₆ are the tangent of half angles of the q₂ and q₃ active joint, respectively; and,
• the mobile platform, is a UPU mechanism with a prismatic joint in-between two universal joints which are connected to the upper and the lower planar mechanisms, respectively.

One note that must be made is that the actuator q₁ generates the translational motion along the Z axis of the entire frame on which the upper and lower planar mechanism are mounted. This motion is only for positioning the robot correctly at the operating table; while the robot operates, the Z motion (the parameter D_n) comes from the medical instruments (attached to the mobile platform) that have redundant DOFs.

3.1. Joint Space Caracterization of the PRO-HEP-LCT Parallel Robotic System

The joint space characterization of PRO-HEP-LCT is achieved using the ideas that are presented in Section 2 of the paper: (1) the two kinematic chains of the parallel robotic system are decoupled and studied individually; (2) new parameterizations are defined based on the loop equations derived in the previous step; and, (3) the motion of the mobile platform is studied based on the equations derived in steps 1–2.

When considering the new defined parameterizations for the robotic system, the target in this case study is to define: the constraints of the point D [X_d, Y_d, Z_d] from the lower planar mechanism; the constraints of the point U [X_u, Y_u, Z_u] from the upper planar mechanism; and finally, the constraints of the mobile coordinate system X’Y’Z’ attached to the mobile platform of the robotic system (see Figure 4).

3.1.1. The Lower Planar Mechanism

For the joint-space characterization of the lower planar mechanism, Equation (4) are referred (the 6R part of the lower planar mechanism has the same topology). To proceed further with the computation, the numerical values for the robot links are substituted in (l₁ = 175, l₂ = 175, d = 120, l₃ = 175, l₄ = 175, l₆ = 430) [mm]. The numerical substitution at this stage has two advantages: (a) it facilitates the computation (reducing the total indeterminates in the equations); and, (b) it avoids the possible Groebner bases “singularities” (that may occur due to the division algorithm used in computing the bases), which were pointed out in [26].

The study parameters for the 6R component of the lower planar mechanism are:

$$M_1:=\left[\begin{array}{c}{x}_0\\ x_1\\ x_2\\ x_3\\ y_0\\ y_1\\ y_2\\ y_3\end{array}\right]={\rho }_1\left[\begin{array}{c}4(({\alpha }_3+{\alpha }_2){\alpha }_1+{\alpha }_2{\alpha }_3-1)\\ 0\\ 0\\ 4({\alpha }_3{\alpha }_2{\alpha }_1-{\alpha }_2-{\alpha }_3-{\alpha }_1)\\ 0\\ 120(({\alpha }_2+\frac{35}{6}){\alpha }_3-1){\alpha }_1-{\alpha }_2-{\alpha }_2+\frac{35}{6}\\ 120(-{\alpha }_2-{\alpha }_3+\frac{35}{6}){\alpha }_1+1-({\alpha }_2+\frac{35}{6}){\alpha }_3\\ 0\end{array}\right],$$

(17a)

$$M_2:=\left[\begin{array}{c}{x}_0\\ x_1\\ x_2\\ x_3\\ y_0\\ y_1\\ y_2\\ y_3\end{array}\right]={\rho }_2\left[\begin{array}{c}4(({\alpha }_5+{\alpha }_4){\alpha }_6+{\alpha }_5{\alpha }_4-1)\\ 0\\ 0\\ 4({\alpha }_6{\alpha }_5{\alpha }_4-{\alpha }_6-{\alpha }_5-{\alpha }_4)\\ 0\\ -980((({\alpha }_5-\frac{5}{7}){\alpha }_4-1){\alpha }_6-{\alpha }_5-{\alpha }_4-\frac{5}{7})\\ -740(({\alpha }_4+{\alpha }_5-\frac{35}{37}){\alpha }_6-1+({\alpha }_5+\frac{35}{37}){\alpha }_4)\\ 0\end{array}\right].$$

(17b)

In order to determine the constraints that the 3R component introduces to the lower planar mechanism, one approach is to determine the relation between the angles α₆ and α₄. Using the study parameter representation for the 3R component of the mechanism yields:

$$M_3:=Q_2(l_6-d)\cdot Q_{15}({\theta }_6)\cdot Q_{16}(l_3)\cdot Q_{17}({\theta }_5)\cdot Q_{18}(l_5)\cdot Q_{19}({\beta }_1)\cdot Q_{20}(d/2).$$

(18)

Subtracting Equation M₁ from Equation M₃ creates a polynomial ideal J₀ that describes the joint space of the two concatenated parallelogram mechanisms (where one side of the mechanism is the 3R component of the lower planar mechanism), and reads:

$$\begin{array}{l}{J}_0:=\langle -4{\rho }_2(({\beta }_1-{\alpha }_5){\alpha }_6+{\alpha }_5{\beta }_1)+4{\rho }_3(({\alpha }_4+{\alpha }_5){\alpha }_6+{\alpha }_4{\alpha }_5-1);-4{\rho }_2(({\alpha }_5{\alpha }_6+1){\beta }_1+{\alpha }_5+{\alpha }_5)+\\ 4{\rho }_3(({\alpha }_4{\alpha }_5-1){\alpha }_6-{\alpha }_5-{\alpha }_4);{\rho }_2(((500{\alpha }_5-700){\beta }_1-500){\alpha }_6-500({\beta }_1+{\alpha }_5)-700)-\\ 980{\rho }_3((({\alpha }_5-\frac{5}{7}){\alpha }_4-1){\alpha }_6-{\alpha }_4-{\alpha }_5-\frac{5}{7});{\rho }_2((740({\beta }_1+{\alpha }_5)-700){\alpha }_6+740{\beta }_1{\alpha }_5+700{\beta }_1-740)-\\ 740{\rho }_3(({\alpha }_4+{\alpha }_5-\frac{35}{7}){\alpha }_6+{\alpha }_4{\alpha }_5+{\alpha }_5\frac{35}{7}-1)\rangle \end{array}$$

(19)

An elimination ideal (J₁ independent of β₁) is computed from the ideal J₀ by computing a Groebner base, as following:

J₁ = G₀[1..4]; G₀ is a Groebner basis, of J₀ with respect to $[{\rho }_3,{\rho }_2,{\beta }_1]$≺$[{\alpha }_6,{\alpha }_5,{\alpha }_4]$ where:

$$J_0:=\langle ({\alpha }_6+{\alpha }_5)({\alpha }_6{\alpha }_5{\alpha }_4-{\alpha }_6-{\alpha }_5-{\alpha }_4);({\alpha }_5^2-1)({\alpha }_6{\alpha }_5{\alpha }_4-{\alpha }_6-{\alpha }_5-{\alpha }_4)\rangle$$

(20)

showing that the constraints corresponding to the 3R component of the lower planar mechanism are given by:

$${\alpha }_6{\alpha }_5{\alpha }_4-{\alpha }_6-{\alpha }_5-{\alpha }_4=0$$

(21)

A union is defined between M₁–M₂ and Equation (21), which represents a polynomial ideal, denoted J_l, describing the joint space of the lower planar mechanism (the intersection of two algebraic varieties is defined by the union of their ideals [20]). In order to eliminate the variables α₄, and α₅ an elimination ideal J_l is computed:

J_l = G₁[1..3]; G₁ is a Groebner basis, of J₀ with respect to $[{\rho }_1,{\rho }_1,{\alpha }_4,{\alpha }_5]$≺$[{\alpha }_2,{\alpha }_3]$, where:

$$J_1:=\langle g_1^{}({\alpha }_1,{\alpha }_2,{\alpha }_3),g_2^{}({\alpha }_1,{\alpha }_2,{\alpha }_3,{\alpha }_6),g_3^{}({\alpha }_1,{\alpha }_2,{\alpha }_3,{\alpha }_6)\rangle$$

(22)

J₁ already describes the constraints of the lower planar mechanism. However, the goal is to obtain equations that describe the solution of the angles α₃ and α₂ in a simpler manner. Consequently, a pure lexicographic Groebner basis is computed from J₁, as follows:

G_s is a Groebner basis, of J₁ with respect to α₂≺α₃ where:

$$\begin{array}{c}{G}_s:=\langle P_1{\alpha }_3^4+P_2{\alpha }_3^3+P_3{\alpha }_3^2+P_2{\alpha }_3^1+P_4;Q_1{\alpha }_2+Q_2{\alpha }_3^3+Q_3{\alpha }_3^2+Q_4{\alpha }_3^1+Q_3\rangle \\ P_1=(961{\alpha }_1^2-2170{\alpha }_1+961){\alpha }_6^2+(2170{\alpha }_1^2-2450{\alpha }_1+2170){\alpha }_6+3411{\alpha }_1^2-2170{\alpha }_1+961;\\ P_2=(2170{\alpha }_1^2-2450{\alpha }_1+2170){\alpha }_6^2+(2450{\alpha }_1^2+2450){\alpha }_6+2170{\alpha }_1^2-2450{\alpha }_1+2170;\\ P_3=(1922{\alpha }_1^2-4340{\alpha }_1+4372){\alpha }_6^2+(4340{\alpha }_1^2-4900{\alpha }_1+4340){\alpha }_6+4372{\alpha }_1^2-4340{\alpha }_1+1992;\\ P_4=(961{\alpha }_1^2-2170{\alpha }_1+3411){\alpha }_6^2+(2170{\alpha }_1^2-2450{\alpha }_1+2170){\alpha }_6+961{\alpha }_1^2-2170{\alpha }_1+3411;\\ Q_1=(961{\alpha }_1^4+1922{\alpha }_1^2+961){\alpha }_6^2+(2170{\alpha }_1^4+4340{\alpha }_1^2+2170){\alpha }_6+961{\alpha }_1^4+1922{\alpha }_1^2+961;\\ Q_2=(961{\alpha }_1^4-2170{\alpha }_1^2+961){\alpha }_6^2+(2170{\alpha }_1^4-2450{\alpha }_1^2+2170){\alpha }_6+3411{\alpha }_1^4-2170{\alpha }_1^3+961{\alpha }_1^2;\\ Q_3=(2170{\alpha }_1^4-1489{\alpha }_1^3+961{\alpha }_1){\alpha }_6^2+(2450{\alpha }_1^4+2170{\alpha }_1^3+2170{\alpha }_1){\alpha }_6+2170{\alpha }_1^4+961{\alpha }_1^3+961{\alpha }_1^{};\\ \begin{array}{l}{Q}_4=(1922{\alpha }_1^4-2170{\alpha }_1^3+2883{\alpha }_1^2+961){\alpha }_6^2+(4340{\alpha }_1^4-2450{\alpha }_1^3+6510{\alpha }_1^2+2170){\alpha }_6\\ +4372{\alpha }_1^4-2170{\alpha }_1^3+2883{\alpha }_1^3+961{\alpha }_1^{}.\end{array}\\ \end{array}$$

(23)

The basis G_s contains two polynomials, the first one being univariate in α₃ of degree 4 with only two solutions being real valued, which is expected, since, for the angles α₁ and α₆, the mechanism may only pose in two ways (one pose being “concave” and the second one “convex”).

At this point, explicit solutions may be computed for the α₂ and α₃ angles (Equation (24)), which may be further used in defining the parameterization for the point D[X_d, Y_d, 0]. The solution that describes the convex working mode of the lower planar mechanism is shown in explicit form in Equation A, and it is further used in the derivation, since it describes the intended working mode of the PRO-HEP-LCT. The disadvantage of disregarding the other solution is that some singularities of the robot might be lost. However, a complete singularity analysis of PRO-HEP-LCT was achieved in previous work [26].

$$\begin{array}{c}{\alpha }_2=\frac{\sqrt{V_1{V}_2}+V_3}{961({\alpha }_1^2+1)({\alpha }_6^2+71/31{\alpha }_6+1)}\\ {\alpha }_3=\frac{\sqrt{V_1{V}_2}+V_4}{P_1}\\ V_1=(961{\alpha }_1^2-2170{\alpha }_1+2186){\alpha }_6^2+(2170{\alpha }_1^2-2450{\alpha }_1+2170){\alpha }_6+2186{\alpha }_1^2-2170{\alpha }_1+961\\ V_2=(264{\alpha }_1^2+2170{\alpha }_1-961){\alpha }_6^2-(2170{\alpha }_1^2-2450{\alpha }_1+2170){\alpha }_6-961{\alpha }_1^2+2170{\alpha }_1+264\\ V_{4,5}=(-1085{\alpha }_1^2+1225{\alpha }_1\pm 1085){\alpha }_6^2-(1225{\alpha }_1^2\mp 1225){\alpha }_6-1085{\alpha }_1^2\mp 1225{\alpha }_1\pm 1085\end{array}$$

(24)

In order to describe the motion of the point D relative to the XYZ coordinate system, Equation (17a) may be evaluated with the explicit solutions shown in Equation (24), and the resulting study parameters may be transformed back into homogeneous matrix form yielding the X_d and Y_d Cartesian coordinates of the point D in explicit form (from the (i = 2, j = 1) and (i = 3, j = 2) entries of the homogeneous matrix):

$$\begin{array}{c}{X}_d=\frac{961(-{\alpha }_1\sqrt{V_1{V}_2}+U_1+961{\alpha }_6^2-2176{\alpha }_6-961)({\alpha }_6^2+70/31{\alpha }_6+1)}{(-2170{\alpha }_6^2+2450{\alpha }_1-2450{\alpha }_6-2170)\sqrt{V_1{V}_2}+2170U_2({\alpha }_6^2+70/31{\alpha }_6+1)}\\ Y_d=\frac{(3007{\alpha }_6^2-1225{\alpha }_1+5565{\alpha }_6+3007)\sqrt{V_1{V}_2}+K_1}{\frac{1}{2}((-2170{\alpha }_6^2+2450{\alpha }_1-2450{\alpha }_6-2170)\sqrt{V_1{V}_2}+2170U_2({\alpha }_6^2+70/31{\alpha }_6+1))}\\ U_1=(1085{\alpha }_1^3-2186{\alpha }_1^2+1085{\alpha }_1){\alpha }_6^2+(1225{\alpha }_1^3-1270{\alpha }_1^2+3675{\alpha }_1){\alpha }_6+1085{\alpha }_1^3-4636{\alpha }_1^2+3255{\alpha }_1\\ U_2=(961{\alpha }_1^2-2170{\alpha }_1+2186){\alpha }_6^2+(2170{\alpha }_1^2-2450{\alpha }_1+2170){\alpha }_6+2186{\alpha }_1^2-2170{\alpha }_1+1\\ \begin{array}{l}{K}_1=961(-3007{\alpha }_1^3+4620{\alpha }_1^2+60142/31{\alpha }_1^{}+2170){\alpha }_6^4+1085(-10943{\alpha }_1^3+11410{\alpha }_1^2+\frac{230208}{31}{\alpha }_1^{}+\\ 6510){\alpha }_6^3+105(\frac{-23823529}{105}{\alpha }_1^3+152573{\alpha }_1^2+\frac{-11817479}{105}{\alpha }_1^{}+115858){\alpha }_6^2+140(-171455{\alpha }_1^3+69440{\alpha }_1^2\\ -\frac{101583}{4}){\alpha }_6^{}-69440-31(287992{\alpha }_1^3-\frac{1761970}{31}{\alpha }_1^2-68483{\alpha }_1^{}+105245)\end{array}\end{array}$$

(25)

3.1.2. The Upper Planar Mechanism

To describe the joint space of the upper planar mechanism, the following equations:

$$N_1:=Q_{21}(D_z)\cdot Q_{22}(q_5)\cdot Q_{23}({\delta }_1)\cdot Q_{24}(l_0)\cdot Q_{25}({\delta }_2)$$

(26a)

$$N_2:=Q_{26}(D_z)\cdot Q_{27}(q_4)\cdot Q_{28}(p_1)$$

(26b)

are subtracted (where the study parameters of the right hand side are detailed in

Table 3. The study parameters used in the motion parameterization of the upper planar mechanism.

Symbol	Explicit form	Details
Q21(Dz)	[1:0:0:0:0:0:0:−Dz/2]	(geometric) translation along OZ
Q22(q5)	[1:0:0:0:0:0:−q5/2:0]	(active)translation along OY
Q23(δ1)	[1:0:0:δ1:0:0:0:0]	(passive)rotation about OZ
Q24(l0)	[1:0:0:0:0:−l0/2:0:0]	(geometric) translation along OX
Q25(δ2)	[1:0:0:δ2:0:0:0:0]	(passive) rotation about OZ
Q26(l0)	[1:0:0:0:0:0:0:−Dz/2]	(geometric) translation along OZ
Q27(q4)	[1:0:0:0:0:0:−q4/2:0]	(active)translation along OY
Q28(p1)	[1:0:0:0:0:−p1/2:0:0]	(free) translation along OX

), yielding:

$$\begin{array}{l}{J}_u:=\langle 2({\rho }_1{\delta }_1{\delta }_2-{\rho }_1-2{\rho }_2);-2{\rho }_1({\delta }_1+{\delta }_2);{\rho }_1(({\delta }_1{\delta }_2+1)l_0+({\delta }_1+{\delta }_2)q_5)+2{\rho }_2{p}_1;\\ {\rho }_1(({\delta }_1-{\delta }_2)l_0+({\delta }_1+{\delta }_2)q_5)+2{\rho }_2{q}_4\rangle \end{array}$$

(27)

where J_u is a polynomial ideal over the joint space of the upper planar mechanism. Computing an elimination ideal:

J_{u_0} = G_{u_0} [1]; G_{u_0} is a Groebner basis, of J_u with respect to $[{\rho }_1,{\rho }_1,{\delta }_1,{\delta }_2]$≺$[p_1,q_4,q_5]$ where:

$$J_{u\_0}:=\langle -l_0^2+p_1^2+{(q_4-q_5)}^2\rangle$$

(28)

yields the explicit solution for p₁:

$$p_1=\pm \sqrt{-l_0^2+{(q_4-q_5)}^2}$$

(29)

Which, if substituted back into Equation (26b), yields the Cartesian coordinates of the point U in explicit form (after computing the homogeneous transformation matrix from the study parameters):

$$\begin{array}{c}{X}_u=\sqrt{(q_4-q_5+l_0)(-q_4+q_5+l_0)}\\ Y_u=q_4\\ Z_u=D_z\end{array}$$

(30)

3.1.3. The Mobile Platform Mechanism

The motion of the mobile platform may be determined easily by using the Cartesian coordinates of the points U and D, respectively (similarly to the approach used in [26]). The two study parameter equations considered are:

$$T_1:=Q_{29}(X_d)\cdot Q_{30}(Y_d)\cdot Q_{31}(c)\cdot Q_{32}({\phi }_1)\cdot Q_{33}({\phi }_2)\cdot Q_{34}(D_n)$$

(31a)

$$T_2:=Q_{35}(D_Z)\cdot Q_{36}(X_u)\cdot Q_{37}(Y_u)\cdot Q_{38}({\lambda }_1)\cdot Q_{39}({\lambda }_1)\cdot Q_{40}(p_0)\cdot Q_{41}(D_n)$$

(31b)

where the study parameters from the right hand side are detailed in

Table 4. The study parameters used in the motion parameterization of the mobile platform.

Symbol	Explicit form	Details
Q29(Xd)	[1:0:0:0:0:−Xd/2:0:0]	(active) translation along OX
Q30(Yd)	[1:0:0:0:0:0:−Yd/2:0]	(active) translation along OY
Q31(c)	[1:0:0:0:0:−c/2:0:0]	(geometric) translation along OX
Q32(φ1)	[1:φ1:0:0:0:0:0:0]	(passive) rotation about OX
Q33(φ2)	[1:0:φ2:0:0:0:0:0]	(passive) rotation about OY
Q34(Dn)	[1:0:0:0:0:0:0:Dn/2]	(active) translation along OZ
Q35(Dz)	[1:0:0:0:0:0:0:−Dz/2]	(geometric) translation along OZ
Q36(Xu)	[1:0:0:0:0:−Xu/2:0:0]	(active) translation along OX
Q37(Yd)	[1:0:0:0:0:0:−Yu/2:0]	(active) translation along OY
Q38(λ1)	[1:λ1:0:0:0:0:0:0]	(passive) rotation about OX
Q39(λ2)	[1:0:λ2:0:0:0:0:0]	(passive) rotation about OY
Q40(p0)	[1:0:0:0:0:0:0:−p0/2]	(pasive) translation along OZ
Q41(Dn)	[1:0:0:0:0:0:0:Dn/2]	(active) translation along OZ

.

Subtracting Equation (31a) from Equation (31b) yields:

$$\begin{array}{l}{J}_p:=\langle 2({\rho }_1-{\rho }_2);2{\phi }_1({\rho }_1-{\rho }_2);-2{\phi }_2({\rho }_1-{\rho }_2);-2{\phi }_2{\phi }_1({\rho }_1-{\rho }_2);\\ (((D_z-p_0-D_n){\rho }_2+{\rho }_1{D}_n){\phi }_2-X_u{\rho }_2+{\rho }_1(X_d+c)){\phi }_1-{\phi }_2(Y_d{\rho }_1-Y_u{\rho }_2);\\ ((-{\phi }_1{Y}_u+D_n+D_z+p_0){\rho }_2+({\phi }_1{Y}_d-D_n){\rho }_2){\phi }_2+X_u{\rho }_2-(X_d+c){\rho }_1;\\ (({\phi }_2{X}_u+D_z+D_n+p_0){\rho }_2-(D_n+(X_d+c){\phi }_2){\rho }_1){\phi }_1-Y_d{\rho }_1+Y_u{\rho }_2;\\ (-{\phi }_2{X}_u-{\phi }_1{Y}_u-D_n+D_z-p_0){\rho }_2+(D_n+(X_d+c){\phi }_2+Y_d{\phi }_1){\rho }_1\rangle \end{array}$$

(32)

where J_p is the polynomial ideal over the joint space of the mobile platform and J_{p_0} is an elimination ideal computed, as follows:

J_{p_0} = G_{p_0}[1..3]; G_{p_0} is a Groebner basis, of J_p with respect to $[{\rho }_1,{\rho }_1,p_0]$≺$[{\phi }_1,{\phi }_2]$, where:

$$\begin{array}{l}{J}_{p\_0}:=\langle (X_d-X_u+c){\phi }_2^2+2((Y_d-Y_u){\phi }_1+D_z){\phi }_2-c-X_d+X_u;(Y_d-Y_u){\phi }_2^2+2D_z{\phi }_1-Y_d+Y_u;\\ {(X_d-X_u+c)}^2{\phi }_2^3+(-X_d^2+2(X_u-c)X_d-{(c-X_u)}^2-4({(Y_d-Y_u)}^2-D_z^2)){\phi }_2+2(X_d-X_u+\\ c)((Y_d-Y_u){\phi }_1+D_z)\rangle \end{array}$$

(33)

Note: to facilitate the computation of the J_{p_0} ideal, the conditions [λ₁ = φ₁, λ₂ = φ₂] are also provided. These conditions are solutions of the first four equations in J_p, are easily deduced, even from the topology of the mechanism and they do not influence the computation results (they only diminish the computation time).

The constraints of the mobile platform (described in the joint space) are already given by Equation (33). However, a pure lexicographic Groebner base is computed from J_{p_0} and it yields simpler constraint equations:

G_{p_1} is a Groebner basis, of J_{p_0} with respect to φ₂≺φ₁, where:

$$G_{p\_1}:=\langle {\phi }_1^2(X_d-Y_d)+2D_z{\phi }_1-Y_d+Y_u;2(Y_d-Y_u){\phi }_1{\phi }_2+(X_d-X_u+c){\phi }_2^2+2{\phi }_2{D}_z-X_d+X_u-c\rangle$$

(34)

with the explicit solutions for the mobile platform orientation:

$$\begin{array}{c}{\phi }_1=\frac{\sqrt{D_z^2+{(Y_d-Y_u)}^2}-D_z}{Y_d-Y_u}\\ {\phi }_2=\frac{\sqrt{D_z^2+{(Y_d-Y_u)}^2}-D_z+\sqrt{D_z^2+(2(c-X_u)+X_d)X_d+{(X_u-c)}^2+{(Y_d-Y_u)}^2}}{X_d-X_u+c}\end{array}$$

(35)

The explicit solutions from Equation (35) may be substituted into Equation (31a), yielding the constraints of the mobile frame (with five DOF) attached to the tip of the medical instrument (the Cartesian coordinates are obtained from the homogeneous transformation matrix, whereas the orientations can be already obtained from Equation (34)):

$$\{\begin{array}{c}{h}_1:=X_n-\frac{(c+X_d){\phi }_2^2+2D_n{\phi }_2+c+X_d}{{\phi }_2^2+1}\\ h_2:=Y_n-\frac{Y_d({\phi }_2^2+1){\phi }_1^2-(2(D_n{\phi }_2^2-1)){\phi }_1+Y_d({\phi }_2^2+1)}{({\phi }_2^2+1)({\phi }_1^2+1)}\\ h_3:=Z_n-\frac{D_n((1-{\phi }_1^2){\phi }_2^2+{\phi }_1^2-1)}{({\phi }_2^2+1)({\phi }_1^2+1)}\\ h_4:={\phi }_1^2(X_d-Y_d)+2D_z{\phi }_1-Y_d+Y_u\\ h_5:=2(Y_d-Y_u){\phi }_1{\phi }_2+(X_d-X_u+c){\phi }_2^2+2{\phi }_2{D}_z-X_d+X_u-c\end{array}$$

(36)

The constraint equations presented in Equation (36) are useable for both forward and inverse geometric models; the coordinates of the mobile platform (for the direct geometric model) are obtained from h₁, h₂, and h₃ (which can be linearly solved for X_n, Y_n, Z_n), whereas the orientations of the mobile frame are found in Equation (35).

For the inverse geometric model, all five equations within Equation (36) are solved simultaneously for the inputs (X_d, Y_d, X_u, Y_u, D_n), yielding:

$$\{\begin{array}{c}{X}_d=\frac{({\phi }_1-1)({\phi }_1+1)(X_n-c)({\phi }_2^2-1)+2Z_n({\phi }_1^2+1){\phi }_2}{({\phi }_2^2-1){\phi }_1^2-{\phi }_2^2+1}\\ Y_d=\frac{Y_n{\phi }_1^2-2Z_n{\phi }_1-Y_n}{{\phi }_1^2-1}\\ X_u=\frac{X_n({\phi }_1^2-1){\phi }_2^2+2({\phi }_1^2+1)(Z_n-D_z){\phi }_2-X_n({\phi }_1^2-1)}{({\phi }_2^2-1){\phi }_1^2-{\phi }_2^2+1}\\ Y_u=\frac{Y_n{\phi }_1^2-2(Z_n-D_z){\phi }_1-Y_n}{{\phi }_1^2-1}\\ D_n=\frac{-Z_n({\phi }_2^2+1)({\phi }_1^2+1)}{({\phi }_2^2-1){\phi }_1^2-{\phi }_2^2+1}\end{array}$$

(37)

where the terms X_d, Y_d, X_u, Y_u, can be substituted in Equation (25) for the lower planar mechanism and in Equation (30) for the lower planar mechanism, respectively, thus achieving the inverse geometric model of the PRO-HEP-LCT parallel robot. The active joints of the robotic system in explicit form are:

$$\{\begin{array}{c}{\alpha }_1=\frac{-350(Y_d+60)+\sqrt{-(X_d^2+{(Y_d+60)}^2)(X_d^2+{(Y_d+60)}^2-{350}^2)}}{(X_d-350)X_d+{(Y_d+60)}^2}\\ {\alpha }_6=\frac{-350(Y_d+370)+\sqrt{-(X_d^2+{(Y_d+370)}^2)(X_d^2+{(Y_d+370)}^2-{350}^2)}}{(X_d-350)X_d+{(Y_d+370)}^2}\\ q_4=Y_u\\ q_5=Y_u+\sqrt{l_0^2-X_u^2}\end{array}$$

(38)

3.2. Velocity Kinematics of the PRO-HEP-LCT Parallel Robotic System

To compute the velocity kinematics of the PRO-HEP-LCT parallel robotic system the following relation is used [21]:

$$A\cdot \dot{X}+B\cdot \dot{Q}=0$$

(39)

where A is the Jacobian matrix computed from Equation (36) with respect to [X_n, Y_n, Z_n, φ₁, φ₂], B is the Jacobian matrix computed with respect to [X_d, X_u, Y_d, Y_u, D_n], X’ is the time dependent vector of mobile platform coordinates, and Q’ is the time dependent vector of the inputs of the mobile platform:

$$\begin{gathered} A=\left[\begin{array}{ccccc}1& 0& 0& 0& \frac{2D_n({\phi }_2^2-1)}{{({\phi }_2^2+1)}^2}\\ 0& 1& 0& \frac{2D_n((1-{\phi }_1^2){\phi }_2^2+{\phi }_1^2-1)}{({\phi }_2^2+1){({\phi }_1^2+1)}^2}& \frac{8D_n{\phi }_1{\phi }_2}{({\phi }_1^2+1){({\phi }_1^2+1)}^2}\\ 0& 0& 1& \frac{4D_n({\phi }_2^2-1){\phi }_1}{({\phi }_2^2+1){({\phi }_1^2+1)}^2}& \frac{4D_n({\phi }_1^2-1){\phi }_2}{({\phi }_1^2+1){({\phi }_2^2+1)}^2}\\ 0& 0& 0& 2((Y_d-Y_u){\phi }_1+D_z)& 0\\ 0& 0& 0& 2(Y_d-Y_u){\phi }_2& 2((X_d-X_u+c){\phi }_2+(Y_d-Y_u){\phi }_1+D_z)\end{array}\right] \\ B=\left[\begin{array}{ccccc}-1& 0& 0& 0& -\frac{2{\phi }_2}{{\phi }_2^2+1}\\ 0& -1& 0& 0& \frac{2{\phi }_1({\phi }_2^2-1)}{({\phi }_2^2+1)({\phi }_1^2+1)}\\ 0& 0& 0& 0& \frac{({\phi }_2^2-1){\phi }_1^2-{\phi }_2^2+1}{({\phi }_2^2+1)({\phi }_1^2+1)}\\ 0& {\phi }_1^2-1& 0& -{\phi }_1^2+1& 0\\ {\phi }_2^2-1& 2{\phi }_1{\phi }_2& -{\phi }_2^2+1& -2{\phi }_1{\phi }_2& 0\end{array}\right] \\ \dot{X}=\left[\begin{array}{c}{\dot{X}}_n\\ {\dot{Y}}_n\\ {\dot{Z}}_n\\ {\dot{\phi }}_1\\ {\dot{\phi }}_2\end{array}\right]\hspace{1em}\hspace{1em}\dot{Q}=\left[\begin{array}{c}{\dot{X}}_d\\ {\dot{Y}}_d\\ {\dot{X}}_u\\ {\dot{Y}}_u\\ {\dot{D}}_n\end{array}\right] \end{gathered}$$

(40)

were, the velocities of the inputs α₁ and α₆ are obtained from the time derivatives of the points U(X_u, Y_u) and D(X_d, Y_d) using the chain rule of the multivariate functions. To solve the velocity kinematics the following matrix relations are used: $\dot{Q}=-B^{-1}\cdot A\cdot \dot{X}$ for the inverse kinematics and $\dot{X}=-A^{-1}\cdot B\cdot \dot{Q}$ for the forward kinematics, respectively. For the inverse kinematics, the explicit solutions are (the results are only presented in a general form due to the length of the equations):

$$\{\begin{array}{c}{\dot{\alpha }}_1=f({\alpha }_1,{\alpha }_6,{\phi }_1,{\phi }_2,D_n,{\dot{\phi }}_1,{\dot{\phi }}_2,{\dot{X}}_n,{\dot{Y}}_n,{\dot{Z}}_n)\\ {\dot{\alpha }}_6=f({\alpha }_1,{\alpha }_6,{\phi }_1,{\phi }_2,D_n,{\dot{\phi }}_1,{\dot{\phi }}_2,{\dot{X}}_n,{\dot{Y}}_n,{\dot{Z}}_n)\\ {\dot{q}}_4=f(Y_d,Y_u,D_n,{\phi }_1,{\phi }_2,{\dot{Y}}_n,{\dot{Z}}_n,{\dot{\phi }}_1)\\ {\dot{q}}_5=f(X_d,X_u,Y_d,Y_u,q_4,q_5,D_n,{\phi }_1,{\phi }_2,{\dot{X}}_n,{\dot{Y}}_n,{\dot{Z}}_n,{\dot{\phi }}_1,{\dot{\phi }}_2)\\ {\dot{D}}_n=f(D_n,{\phi }_1,{\phi }_2,{\dot{\phi }}_1,{\dot{\phi }}_2,{\dot{Z}}_n)\end{array}$$

(41)

and for the forward kinematics the explicit solutions are:

$$\{\begin{array}{c}{\dot{X}}_n=f(D_n,{\alpha }_1,{\alpha }_6,q_4,q_5,{\dot{D}}_n,{\dot{\alpha }}_1,{\dot{\alpha }}_6,{\dot{q}}_4,{\dot{q}}_5,{\phi }_1,{\phi }_2)\\ {\dot{Y}}_n=f(D_n,{\alpha }_1,{\alpha }_6,q_4,q_5,{\dot{D}}_n,{\dot{\alpha }}_1,{\dot{\alpha }}_6,{\dot{q}}_4,{\dot{q}}_5,{\phi }_1,{\phi }_2)\\ {\dot{Z}}_n=f(D_n,{\alpha }_1,{\alpha }_6,q_4,q_5,{\dot{D}}_n,{\dot{\alpha }}_1,{\dot{\alpha }}_6,{\dot{q}}_4,{\dot{q}}_5,{\phi }_1,{\phi }_2)\\ {\dot{\phi }}_1=f({\alpha }_1,{\alpha }_6,q_4,{\dot{\alpha }}_1,{\dot{\alpha }}_6,{\dot{q}}_4,{\phi }_1)\\ {\dot{\phi }}_2=f({\alpha }_1,{\alpha }_6,q_4,q_5,{\dot{\alpha }}_1,{\dot{\alpha }}_6,{\dot{q}}_4,{\dot{q}}_5,{\phi }_1,{\phi }_2)\end{array}$$

(42)

where f (in Equations (41) and (42)) represents an irrational function.

3.3. Acceleration Kinematics of the PRO-HEP-LCT Parallel Robotic System

To compute the acceleration kinematics of the PRO-HEP-LCT parallel robotic system Equation (39) must be differentiated with respect to time, yielding:

$$A\cdot \ddot{X}+\dot{A}\cdot \dot{X}+\dot{B}\cdot \dot{Q}+B\cdot \ddot{Q}=0$$

(43)

where the time dependent Jacobians and the acceleration vectors are:

$$\begin{gathered} \dot{A}=\left[\begin{array}{ccccc}0& 0& 0& 0& \dot{A}[1,5]\\ 0& 0& 0& \dot{A}[2,4]& \dot{A}[2,5]\\ 0& 0& 0& \dot{A}[3,4]& \dot{A}[3,5]\\ 0& 0& 0& 2(Y_d-Y_u){\dot{\phi }}_1+2{\phi }_1({\dot{Y}}_d-{\dot{Y}}_u)& 0\\ 0& 0& 0& 2(Y_d-Y_u){\dot{\phi }}_2+2{\phi }_2({\dot{Y}}_d-{\dot{Y}}_u)& 2((X_d-X_u+c){\dot{\phi }}_2+(Y_d-Y_u){\dot{\phi }}_1+2{\phi }_1({\dot{Y}}_u-{\dot{Y}}_d))\end{array}\right] \\ \dot{A}[1,5]=-\frac{4D_n({\phi }_2^3-3{\phi }_2){\dot{\phi }}_2-2({\phi }_2^4-1){\dot{D}}_n}{{({\phi }_2^2+1)}^3} \\ \dot{A}[2,4]=\frac{4D_n(((1-{\phi }_2^4){\phi }_1^3+3(1-{\phi }_2^4){\phi }_1){\dot{\phi }}_1+2{\phi }_2({\phi }_1^4-1){\dot{\phi }}_2)+(2({\phi }_2^4+1){\phi }_1^4+2{\phi }_2^4-2){\dot{D}}_n}{{({\phi }_2^2+1)}^2{({\phi }_1^2+1)}^3} \\ \dot{A}[2,5]=\frac{8D_n((({\phi }_1^2-1){\phi }_2^3+({\phi }_1^2-1){\phi }_2){\dot{\phi }}_1+((3{\phi }_2^2-1){\phi }_1^3+(3{\phi }_2^2-1){\phi }_1){\dot{\phi }}_2+({\phi }_2^2+1)({\phi }_1^2+1){\phi }_1{\phi }_2{\dot{D}}_n)}{{({\phi }_2^2+1)}^3{({\phi }_1^2+1)}^2} \\ \dot{A}[3,4]=\frac{4D_n(3(1-{\phi }_1^2){\phi }_1^2+{\phi }_2^4-1){\dot{\phi }}_1+4{\phi }_1{\phi }_2({\phi }_1^2+1){\dot{\phi }}_2+({\phi }_1^2+1)({\phi }_2^2-1){\phi }_1{\dot{D}}_n)}{{({\phi }_2^2+1)}^2{({\phi }_1^2+1)}^3} \\ \dot{A}[3,5]=\frac{4D_n(4{\phi }_1{\phi }_2({\phi }_2^2+1){\dot{\phi }}_1+(3({\phi }_1^4-1){\phi }_2^2-{\phi }_1^4-1){\dot{\phi }}_2-1/2({\phi }_2^2-1)({\phi }_1^2-1){({\phi }_2^2+1)}^2{\dot{D}}_n)}{{({\phi }_1^2+1)}^2{({\phi }_2^2+1)}^3} \\ \dot{B}=\left[\begin{array}{ccccc}0& 0& 0& 0& \frac{2({\phi }_2^2-1){\dot{\phi }}_2}{{({\phi }_2^2+1)}^2}\\ 0& 0& 0& 0& \dot{B}[2,5]\\ 0& 0& 0& 0& \frac{4({\phi }_1({\phi }_2^4-1){\dot{\phi }}_1+{\phi }_2({\phi }_1^4-1){\dot{\phi }}_2)}{{({\phi }_2^2+1)}^2{({\phi }_1^2+1)}^2}\\ 0& 2{\phi }_1{\dot{\phi }}_1& 0& -2{\phi }_1{\dot{\phi }}_1& 0\\ 2{\phi }_2{\dot{\phi }}_2& 2({\phi }_1{\dot{\phi }}_2+{\phi }_2{\dot{\phi }}_1)& -2{\phi }_2{\dot{\phi }}_2& -2({\phi }_1{\dot{\phi }}_2+{\phi }_2{\dot{\phi }}_1)& 0\end{array}\right] \\ \dot{B}[2,5]=\frac{2(((1-{\phi }_2^2){\phi }_1^2+{\phi }_2^4-1){\dot{\phi }}_1+4({\phi }_1{\phi }_2({\phi }_1^2+1)){\dot{\phi }}_2)}{{({\phi }_2^2+1)}^2{({\phi }_1^2+1)}^2} \\ \ddot{X}=\left[\begin{array}{c}{\ddot{X}}_n\\ {\ddot{Y}}_n\\ {\ddot{Z}}_n\\ {\ddot{\phi }}_1\\ {\ddot{\phi }}_2\end{array}\right]\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\ddot{Q}=\left[\begin{array}{c}{\ddot{X}}_d\\ {\ddot{Y}}_d\\ {\ddot{X}}_u\\ {\ddot{Y}}_u\\ {\ddot{D}}_n\end{array}\right] \end{gathered}$$

(44)

were, the accelerations of the inputs α₁ and α₆ are obtained from the second time derivatives of the points U(X_u, Y_u) and D(X_d, Y_d) while using the chain rule of the multivariate functions (for second order differentiation). To solve the inverse kinematics (for the accelerations), the $\ddot{Q}=-B^{-1}\cdot (A\cdot \ddot{X}+\dot{A}\cdot \dot{X}+\dot{B}\cdot \dot{Q})$ relation is used for the inverse kinematics and $\ddot{X}=-A^{-1}\cdot (B\cdot \ddot{Q}+\dot{A}\cdot \dot{X}+\dot{B}\cdot \dot{Q})$ for the forward kinematics, respectively (where explicit solutions may be computed for both cases).

4. Numerical Results

Based on the kinematic models that are derived in Section 3, a numerical simulation was achieved for the PRO-HEP-LCT US guiding module. The simulation consists of guiding the tip of the US probe on a prescribed trajectory in the robot operational workspace. The trajectory is defined, as following:

• the initial conditions are defined by the entry point (0) (see Figure 5) that has the Cartesian coordinates [X = 225, Y = 170, Z = 0] (dimensions in mm) where the US probe has the orientation [ϕ₁ = 0°, ϕ₁ = 45°] (dimensions in angular degrees);
• the entry point is defined as a Remote Center of Motion (RCM) point;
• the probe is inserted from point (0) to point (1) until it reaches Z = −50 by actuating the redundant DOF described by the D_n parameter;
• the tip of the probe is guided on a circular path (counter-clockwise direction) until it reaches the point (2) (a rotation of 3π/2);
• the probe is further inserted on a linear trajectory until it reaches the point (3) at Z = −100;
• the tip of the probe is guided on a full circle (counter-clockwise direction), from point (3) back to point (3);
• the tip of the probe is guided on a linear trajectory until it reaches the point (4); at this point, the tip of the probe has the Cartesian coordinates [X = 225, Y = 170, Z = −100] and the orientations [ϕ₁ = 0°, ϕ₁ = 0°]; and,
• the probe is retracted from the patients’ body.

In order to compute the inverse kinematics of the PRO-HEP-LCT with respect to the trajectory previously defined, the mathematical conditions for the Cartesian coordinates of the mobile coordinate frame and the orientations were substituted in the kinematic models derived in Section 3. The trajectory between the points (1) and (2) is defined by a three quarters of a circle with 50 mm radius, whereas the trajectory starting (and ending) at point (3) is a full circle with 100 mm radius. The velocity for the linear insertion of the US probe was set to 20 mm/s and a maximum (absolute value) of 50 mm/s for the probe retraction; the velocity of the probe tip while guided on the curve from the point (1) to the point (2) was 12 mm/s, and the velocity of the probe tip while guided on the curve from the point (3) back to the point (3) was 32 mm/s. Figure 6 presents the results for the actuated joints (displacement, velocity, and acceleration). It is worth noting that no input profiles for the actuators velocities (triangle or trapezoidal profiles) were used in the computation. Although implementing these profiles is necessary for the development of the control system, the numerical results that are provided in this section are mainly used to validate the analytical models derived in previous sections.

5. Discussion

When considering different mathematical methods used for mechanism analysis and synthesis, it should not be hard to assume that certain methods are better suited for answering specific questions in mechanism analysis. The study parameters method proved to be powerful due to the fact that it describes the global kinematics (showing all the working modes and singularities). A deviation from the study parameters method was previously introduced in [15] to solve a specific problem for a “simple” hybrid rehabilitation robot, and it was later optimized (by solving the equations using algebraic techniques) to allow efficient computation for the constraints of a “more complex” spatial rehabilitation robot [24]. In both cases, the mathematical information describing the motion of a mobile platform had no practical values (the motion amplitude, velocity, etc. of the anatomic joint is more relevant for the medical personal that operates rehabilitation robots).

Referring to Observation 2 and the example based on the 6R closed loop mechanism (Section 2.4), it can be pointed out that the joint-space characterization by the proposed method has a canonical form and it is invariant to the initial coordinate system definition. As long as the motion parameters for the kinematic chains are consistently defined (e.g., the 0° angle value between two links is defined at the 180° angle between the links), it does not matter where the coordinate frames are initially placed. Furthermore, the proposed approach that reintroduces the coordinates of a mobile platform into question allows for the attachment of the mobile coordinate system on any element of the parallel robot. The main point here is that, after obtaining the joint-space equations, one can redefine other parameterizations (by placing a mobile frame on any of the mechanism components) describing the motion of the “mobile platform” (with Cartesian coordinates and orientations). To give an example where the proposed method might be useful, consider that the PRO-HEP-LCT parallel robotic system can be adapted (with slight modifications) to have different points where the mobile platform is attached to the lower planar mechanism (e.g., the universal joint may be mounted on the link l₄ after slight modifications). This reconfiguration has the benefit that redefines the operational workspace (closer the robot base) and, in some cases, this offers advantages in the medical procedure (when the target tissue is located close to the robot). The joint space is already defined for the parallel robot and the motion parameterization of the mobile platform (in the reconfigured way) may be computed with ease.

It is possible to describe the higher order kinematics while using the study parameters (see e.g., [30]), or even the dynamics of mechanical systems [31] or define trajectories for optimal control [32]. However, for the PRO-HEP-LCT parallel robotic system, the higher order kinematics were computed with simple geometric models and the Jacobian matrices representation of kinematic systems. The approach used for velocity and acceleration kinematics was the same as in [33] (where the velocity and acceleration kinematics were derived with a classic vector method). The main difference is that, in [33], the authors achieved the kinematic modeling after imposing extra geometric constraints for specific revolute joints (e.g., adding a π value for the arguments of some trigonometric functions). Although, in this work the authors also restricted the working modes of the robot by choosing specific solutions for the lower planar mechanism (see Section 3.1), the solution was chosen naturally based on an analytic result.

The kinematic models obtained with the proposed method may be implemented in the real time control of the PRO-HEP-LCT, since (even for the complexity of the robotic system) the kinematic equations were solved by a standard PC (using Maple on a Intel Core I 7 with 16 GB RAM) in less than 0.1 s. However, further research is needed to benchmark the results of the proposed method and compare it with other methods to evaluate the feasibility (referring to the control development) of the proposed method.

Both vector-based methods and algebraic ones and their applications in mechanism analysis of parallel medical robots is well established (e.g., in surgical parallel robots [34,35,36,37,38,39]), whereas the advantages and disadvantages of the proposed method and the domain in which it may be useful are not fully determined yet. Further work is intended to investigate how the method could be used for optimization problems in mechanism analysis (e.g., choosing the optimal joints to be active ones), for robotic applications where a control with high accuracy and robustness of a parallel robot is needed (e.g., micro-manipulation, aerospace), dynamic studies for rehabilitation mechanisms that are based on the dynamic characterizations of anatomic joints (e.g., the knee [40]), or other engineering applications where parallel robots are involved. It is worth pointing that the singularities of parallel mechanism are already solvable with the joint-space method.

6. Conclusions

In this paper, the authors proposed a mathematical method that is based on the loop equations that describe the joint space of a parallel robot, and it has the potential to create motion parameterizations (in Cartesian space) for any element of the robot (link, joints, and mobile platform). The method uses the study parameters of SE(3) to describe the joint-space (by eliminating the coordinates of the mobile platform) and creates loop equations that are invariant to the initial placement of the coordinate systems (fixed and mobile). The proposed method was used in the kinematic study of a parallel robotic system (PRO-HEP-LCT) that was designed for minimally invasive liver cancer treatment, and the results show that the equations are usable in the control of the robot (for displacement, velocity, and acceleration).

Further work is intended to study the feasibility of the proposed method for kinematic modeling of parallel robots, when considering the safety in their operation.