Fabrication and in vivo evaluation of Ti6Al4V implants with controlled porous structure and complex shape

Abstract

Electron beam melting process was used to fabricate porous Ti6Al4V implants. The porous structure and surface topography of the implants were characterized by scanning electron microscopy (SEM) and digital microscopy (DM). The results showed that the pore size was around 600 and the porosity approximated to 57%. There was about±50 μm of undulation on implants surfaces. Standard implants and a custom implant coupled with porous sections were designed and fabricated to validate the versatility of the electron beam melting (EBM) technique. After coated with bone-like apatite, samples with fully porous structures were implanted into cranial defects in rabbits to investigate the in vivo performance. The animals were sacrificed at 8 and 12 weeks after implantation. Bone ingrowth into porous structure was examined by histological analysis. The histological sections indicated that a large amount of new bone formation was observed in porous structure. The newly formed bone grew from the calvarial margins toward the center of the bone defect and was in close contact with implant surfaces. The results of the study showed that the EBM produced Ti6Al4V implants with well-controlled porous structure, rough surface topography and bone-like apatite layer are beneficial for bone ingrowth and apposition.

Keywords

electron beam melting process / implant / porous structure / bone ingrowth

1 1 Introduction

Parallel manipulators (PMs), with their notable advantages of high payload, precision, and stiffness, have been widely used in industrial production, aerospace, and many other fields [1‒3]. The design, calibration, and control of PMs are crucial aspects of their production and application, and forward kinematics (FK) is often closely related to these aspects. For example, based on FK equations, we can calculate the current position and posture (collectively referred to as “pose”) of a PM by reading the encoders of the actuated joints, thereby providing precise and reliable information for spatial positioning or accuracy calibration. Over the past few decades, the issue of the FK analysis of PMs has received extensive attention. Whether it is the analysis of specific configurations or research on modeling and solving kinematic equations, substantial experience and resources have been accumulated, contributing significantly to advancements in this field. The current paper focuses on an original analysis technique for the FK analysis of PMs, known as the finite-step-integration (FSI) method. This method, introduced in a past study [3], demonstrates outstanding versatility and dependable computational efficiency. Although FSI has been initially applied to tackle a particular engineering issue without a systematic theoretical foundation, the current work aims to investigate the underlying principles and modeling theories of the FSI method, thus presenting new insights and findings.

Prior to formalizing the theoretical framework of the FSI method, reaching a consensus on PM terminology and FK analysis concepts is necessary. A PM typically consists of a fixed chassis, a movable end-effector, and motion branches (Fig.1). Actuated joints in conventional PMs, called actuators, are typically revolute or prismatic, with parameters representing joint angles or lengths. These parameters serve as input for FK analysis, resulting in the end-effector’s position and posture, which are collectively termed as pose.

Fig.1 Topology structure and three basic components of general PMs.

Full size|PPT slide

A rigid body’s pose can be represented by a 4 × 4 homogeneous matrix. Merlet [4] introduced two key concepts in the FK problem of Gough-type PMs: being given the six leg lengths, (${\mathcal{P}}_1$) find the current pose of the end-effector, or (${\mathcal{P}}_2$) find all the possible poses. While most scholars focus on problem ${\mathcal{P}}_2$, which is more theoretical, practical PM poses are typically confined within their accessible workspace [5,6]. Consequently, problem ${\mathcal{P}}_1$ holds greater research value and can be decomposed into three steps: (${\mathcal{S}}_1$) formulating kinematic equations, (${\mathcal{S}}_2$) solving them, and (${\mathcal{S}}_3$) obtaining real solutions consistent with the PM assembly configuration.

For step ${\mathcal{S}}_1$, kinematic equations that are rooted in the geometric motion laws of spatial rigid bodies can be formulated across diverse algebraic systems, such as classical geometric algebra [7,8], screw theory [9], and conformal geometric algebra [10]. Furthermore, the choice of algebraic representation directly impacts the complexity of step ${\mathcal{S}}_2$.

From a geometric perspective, problem ${\mathcal{P}}_1$ involves determining all vertices of the current spatial polytope. Classical geometric algebra-based constraint equations are adept at capturing polytope geometries. However, these equations often feature multi-variables, high dimensionality, and strong coupling, thus posing challenges for simplification and solution [11]. Existing analytical methods demonstrate their capabilities when solving equations of varying complexity. When equations involve fewer parameters or lower coupling between parameters, using elimination methods (e.g., Sylvester’s dialytic elimination method [12,13]) to decouple the equations is appropriate, as doing so can produce univariate equations while preserving the accuracy of the analytical equations. In comparison, when the equations are highly complex, the elimination method can be arduous and may not yield results; thus, its general applicability needs improvement. Related to this, the interval analysis proposed by Merlet [4] and Chablat et al. [14] remains a versatile method. Such an approach divides the solution space into multiple intervals and uses the boundary conditions of these intervals for search and filtering. Clearly, interval analysis must be combined with numerical algorithms, such as the Newton method, and incorporating Kantorovich’s theorem and inflation operators can ensure the uniqueness of solutions within an interval. The multivariate homotopy continuation method [15] can trace the homotopy path from decoupled equations to target equations, while related homotopy continuation methods can effectively find isolated real solutions of polynomial equations [11,16]. This class of methods has guided the design of polynomial tools, such as PHCpack [17] and Bertini [18]. In addition, Shen et al. [19] have introduced positional and orientation characteristic equations for FK analysis of partially decoupled PMs, providing an effective method for computing analytical solutions with symbolic parameters. Koupaei and Hosseini [20] combined the capabilities of chaotic maps with the golden section search method, transforming the solution process of non-linear equations into an optimization and search process. Liu et al. [21] used Jacobian elliptic function expansions to convert non-linear equations into non-degenerate elliptic function solutions, obtaining soliton and singular solutions. Omran and Kassem [22] used a predictive squared error cost function to seek the optimal polynomial approximation model for the FK of a Stewart platform. Overall, these methods exhibit strong mathematical logic but still have some distance from engineering applications.

The development of computer science has enabled the use of iterative and optimization-based intelligent algorithms for equation solving. These algorithms typically transform problem ${\mathcal{P}}_1$ into a task of testing, evaluating, and optimizing random samples in the solution space. Wu et al. [23] used and improved a back propagation neural network using the all class-one network strategy to establish an FK model for a hybrid manipulator. Zhu and Zhang [24] proposed a novel algorithm that integrated artificial neural networks with the global Newton–Raphson with monotonic descent algorithm, effectively reducing the training set of neural networks and preventing divergence issues. Wang et al. [25] applied the differential evolution method in this context, demonstrating its superiority over genetic algorithms and particle swarm optimization (PSO) through comparative analysis. Similarly, Mao et al. [26] devised a strategy combining the differential evolution method with enhanced PSO, thus leveraging the strengths of both approaches. Despite the availability of various methods for solving FK problems, researchers often invest significant time in kinematic analysis, particularly when dealing with novel mechanisms. This approach is often taken due to the complexity of certain scenarios, such as end-effectors coupling translation and rotation motions [27], which pose challenges for providing input conditions to the inverse kinematic equations.

In summary, existing analytical methods are feasible only for analyzing models with lower complexity; thus displaying limitations in their scope of use. Numerical methods and intelligent algorithms are commendable for their versatility and solving capabilities, but their core concepts are often constrained by trial-and-error limitations. Such a constraint necessitates further optimization in their dependency on inverse kinematics equations and the robustness of the algorithms.

Problem ${\mathcal{P}}_1$ has been rephrased as follows in Ref. [28]: (${\mathcal{P}}_3$) determining the motion of the end-effector as a function of joint motions, rather than seeking one (or more) solutions corresponding to a given set of joint variables. The underlying logic of this viewpoint can be understood as solving the process of low-dimensional variations within a higher-dimensional space by introducing a temporal dimension. Inspired by this perspective, we established the fundamental principles and feasibility of the FSI method in our initial research.

In this paper, R represents a revolute joint, P represents a prismatic joint, S represents a spherical joint, and U represents a universal joint. Following research and synthesis, Section 2 of this paper will establish a rigorous mathematical theoretical framework for the FSI method, providing standardized representations for relevant terminologies. In Section 3, we summarize four modeling methods for linear systems and elucidate their geometric interpretations. Section 4 primarily discusses the identification of singularity in mechanical configurations and associated principles during the solving process. In Section 5, we demonstrate the modeling process and results of the FSI method through two general cases, namely 6-UPS (unconstrained configuration) and 3-UPS/S (lower-mobility configuration), validating its timeliness and accuracy through numerical computations. Section 6 contrasts FSI’s strengths and weaknesses with existing algorithms, while Section 7 uses numerical examples from earlier sections to elucidate FSI’s application in tracking motion characteristic curves. The FSI method can virtually encompass all analytical aspects related to parallel mechanism kinematics, effectively addressing the long-standing challenge of efficiently modeling diverse PMs with a unified mathematical theory. Finally, we summarize our findings in Section 8.

2 2 Establishment of the modeling environment

For problem ${\mathcal{P}}_3$, the motion of a spatial rigid body occurs in six dimensions (6D: three rotations and three translations), which is challenging to intuitively express in 3D space. Fortunately, 3D space can represent points, lines, and volumes; thus, the pose of a rigid body can be entirely determined by the coordinates of its three non-collinear points. During motion, these points (denoted as mark points) generate spatiotemporal curves, as illustrated in Fig.2.

Fig.2 Sweeping process of a rigid body’s continuous motion.

Full size|PPT slide

In this paper, we employ a set of symbolic algebraic systems to describe the relevant geometric elements and motions, as shown in Tab.1.

Tab.1 Symbolic algebraic system of this paper

Interpretation	Symbol
World frame matrix fixed to the chassis, with associated element set $\left\{O,\mathit{x},\mathit{y},\mathit{z}\right\}$	$\mathit{P}$
Dynamic frame matrix fixed to the end-effector, with associated element set $\left\{O_{\text{D}},{\mathit{x}}_{\text{D}},{\mathit{y}}_{\text{D}},{\mathit{z}}_{\text{D}}\right\}$	${\mathit{P}}_{\text{D}}$
Number of actuators	$m$
Motion parameter of the ith actuator	$c_{i⟨\circ ⟩}$
Initial parameter set $\left\{c_{1⟨0⟩},c_{2⟨0⟩},\dots ,c_{m⟨0⟩}\right\}$, which is a known condition	${\mathbb{C}}_{\circ ⟨0⟩}$
Target parameter set $\left\{c_{1⟨\text{E}⟩},c_{2⟨\text{E}⟩},\dots ,c_{m⟨\text{E}⟩}\right\}$, which is a given condition	${\mathbb{C}}_{\circ ⟨\text{E}⟩}$
Total motion time	$t_{\text{E}}$
Time domain $\left[0,t_{\text{E}}\right]$ of the motion	$T$
Function of the parameter $c_{i⟨\circ ⟩}$ with respect to time t ($t\in T$); thus, there are $f_i\left(0\right)=c_{i⟨0⟩}$ and $f_i\left(t_{\text{E}}\right)=c_{i⟨\text{E}⟩}$	$f_i(t)$
Function set $\left\{f_1(t),f_2(t),\dots ,f_m(t)\right\}$	$\mathbb{F}$
Number of mark points, where $h⩾3$	$h$
Position node of the ith mark point at a certain time node; e.g., while $t=0$, there is $P_{i⟨0⟩}$	$P_{i⟨\circ ⟩}$
Initial node set $\left\{P_{1⟨0⟩},P_{2⟨0⟩},\dots ,P_{h⟨0⟩}\right\}$, which is a known condition like ${\mathbb{C}}_{\circ ⟨0⟩}$	${\mathbb{P}}_{\circ ⟨0⟩}$
Target node set $\left\{P_{1⟨\text{E}⟩},P_{2⟨\text{E}⟩},\dots ,P_{h⟨\text{E}⟩}\right\}$, which is the solution of problem ${\mathcal{P}}_1$	${\mathbb{P}}_{\circ ⟨\text{E}⟩}$
Spatiotemporal curve of all nodes $P_{i⟨\circ ⟩}$, i.e., a curve is a set of countless nodes	${\mathbb{L}}_i$
Curve set $\left\{{\mathbb{L}}_1,{\mathbb{L}}_2,\dots ,{\mathbb{L}}_h\right\}$, which is the solution of problem ${\mathcal{P}}_3$	${\mathbb{L}}_{\text{PM}}$

It is important to note that the symbol “$\circ$” within the system serves as a placeholder, facilitating standardized expression. Therefore, we can establish the environment for building mathematical models with the following steps:

1) Consider a PM system with m actuators, where the motion parameter of each actuator is denoted by $c_{i⟨\circ ⟩}$. If the actuator corresponds to a prismatic joint, $c_{i⟨\circ ⟩}$ represents the joint length; if it corresponds to a revolute joint, $c_{i⟨\circ ⟩}$ represents the joint angle. Here, the variable i represents the index of the motion branch.

2) Assume that the PM has an accurately defined initial pose (a feasible aspect in practical design), denoted as ${\mathbb{C}}_{\circ ⟨0⟩}$, representing the set of motion parameters for all actuators at this moment.

3) Based on the nature of FK problems, although the desired end-effector pose is unknown, once the PM reaches the target pose, the motion parameters of all actuators can be represented as the set ${\mathbb{C}}_{\circ ⟨\text{E}⟩}$.

4) The ith actuator is controlled by a motion function $f_i(t)$, and the set of these functions is denoted as $\mathbb{F}$, which shares the time domain $t\in T$.

5) During the motion, a sequence of position points left along the spatiotemporal curve by the marked point $P_i$ on the moving component is referred to as node $P_{i⟨\circ ⟩}$.

6) As time $t$ progresses, the marked point traces out the corresponding spatiotemporal curve before reaching its final position $P_{i⟨\text{E}⟩}$.

7) The countless nodes of the marked point $P_i$ constitute the spatiotemporal curve ${\mathbb{L}}_i$, and the collection of all ${\mathbb{L}}_i$ can be denoted as ${\mathbb{L}}_{\text{PM}}$.

Based on the aforementioned conditions and definitions, the mapping is as follows:

(1)$\mathcal{T}:\left[\begin{array}{cc}{\mathbb{P}}_{\circ ⟨0⟩}& \mathbb{F}\end{array}\right]\to {\mathbb{L}}_{\text{PM}}.$

Each ${\mathbb{L}}_i$ is continuous in 3D space, with its analytical formula being complex and non-linear. To eliminate non-linear factors, the differential theory is introduced, enabling us to analyze the geometric configuration of ${\mathbb{L}}_{\text{PM}}$ within a microscopic time interval. Therefore, the following additional definitions are provided:

1) The time interval $T$ is divided into n segments, and the kth subinterval is denoted as $T_{\circ ⟨k⟩}=[(k-1)t_{\text{E}}\text{/}n,k{t}_{\text{E}}\text{/}n]$.

2) When $t=k{t}_{\text{E}}/\phantom{k{t}_{\text{E}}n}n$, the motion parameters of all actuators are denoted as $c_{i⟨k⟩}$, and the nodes of mark point $P_i$ are denoted as $P_{i⟨k⟩}$. Hence, $c_{i⟨k⟩}=f_i(k{t}_{\text{E}}\text{/}n)$.

3) Following the definition of sets ${\mathbb{C}}_{\circ ⟨0⟩}$ and ${\mathbb{P}}_{\circ ⟨0⟩}$, we define sets ${\mathbb{C}}_{\circ ⟨k⟩}$ and ${\mathbb{P}}_{\circ ⟨k⟩}$. Hence, ${\mathbb{C}}_{\circ ⟨n⟩}={\mathbb{C}}_{\circ ⟨\text{E}⟩}$ and ${\mathbb{P}}_{\circ ⟨n⟩}={\mathbb{P}}_{\circ ⟨\text{E}⟩}$.

4) Within the interval $T_{\circ ⟨k⟩}$, the scalar $d_{i⟨k⟩}=c_{i⟨k⟩}-c_{i⟨k-1⟩}$ and the set ${\mathbb{D}}_{\circ ⟨k⟩}=\{d_{1⟨k⟩},d_{2⟨k⟩},\dots ,d_{m⟨k⟩}\}$ are defined.

5) Within the interval $T_{\circ ⟨k⟩}$, the vector ${\mathit{q}}_{i⟨k⟩}=\rightharpoonup P_{i⟨k-1⟩}{P}_{i⟨k⟩}$ and the set ${\mathbb{Q}}_{\circ ⟨k⟩}=\{{\mathit{q}}_{1⟨k⟩},{\mathit{q}}_{2⟨k⟩},\dots ,{\mathit{q}}_{h⟨k⟩}\}$ are defined.

The geometric representations and algebraic structures of the aforementioned definitions are illustrated in Fig.3, while Fig.4 present the modeling and computation flowchart of the FSI method. Clearly, while finite segmentation results in some loss of precision in the motion process, when the value of n is sufficiently large, the polyline formed by $P_{i⟨\circ ⟩}$ is able to approximate the curve ${\mathbb{L}}_i$, which is the core of the FSI method.

Fig.3 Principle explanation of the FSI method: (a) geometric tracking process and (b) algebraic structure.

Full size|PPT slide

Fig.4 Modeling and computation flowchart of the FSI method.

Full size|PPT slide

3 3 Step motion equations

Based on the algebraic environment established above, it becomes evident that the motion trajectory of this PM can be subdivided into n discrete steps. Further inspired by mathematical induction, it is discernible that if we establish a parameterized iterative algebraic system capable of resolving motion variations at each step—thereby representing all motion vectors of mark points—then the sought-after ${\mathbb{P}}_{\circ ⟨n⟩}$ can be obtained by successively superimposing $n$ motion vector sets ${\mathbb{Q}}_{\circ ⟨k⟩}$ onto ${\mathbb{P}}_{\circ ⟨0⟩}$. This solving process resembles integration computation; however, step durations cannot approach infinitesimal values akin to integration. Consequently, the number of steps for solving is finite, leading to certain errors in the final result. Nonetheless, this seemingly impractical method, following extensive computations and validations, has demonstrated its robust modeling capability and reliable computational efficacy. In the aforementioned solving process, the central challenge lies in establishing an algebraic structure that is fixed yet parameterizable for an iterative step operation system. Moreover, this operation system must have high computational efficiency to compensate for the potential drawback of excessive steps. Therefore, we propose four methods for establishing linear step motion equations within the time interval $T_{\circ ⟨k⟩}$, and these equations are composed into a linear system of equations. This linear system is referred to as the “step equation system,” with its coefficient matrix termed as the “step matrix”.

The step equations established using these four methods are related to actuator variables, the geometric motion of moving components, the stable relationship of vectors on the rigid body structures, and the specific dimensional layout of the configuration, respectively. The step equations related to actuator variables are denoted as ${\text{E}}_{\text{A}}$. Notably, when the actuator is directly connected to the chassis, changes in the actuator can be directly reflected in the motion of the next joint. In this case, ${\text{E}}_{\text{A}}$ can be interchanged with the second type of equation, which is related to the geometric motion of the moving components, denoted as ${\text{E}}_{\text{M}}$. The step equations related to vectors on the rigid body, denoted as ${\text{E}}_{\text{P}}$, use the stable projection quantity between two vectors on the rigid body to establish the equations, thus reflecting the characteristics of the rigid body. The fourth type of equations is denoted as ${\text{E}}_{\text{S}}$ and can only be used in special circumstances. Step equations essentially constrain step motions; thus, this classification is not strictly defined. The specific modeling principles and methods are outlined as follows:

● Equations arising from actuator variations (${\text{E}}_{\text{A}}$)

In step k, the scalar set ${\mathbb{D}}_{\circ ⟨k⟩}$ is the primary cause of motion occurrence; thus, there must exist a mapping relationship between the vector set ${\mathbb{Q}}_{\circ ⟨k⟩}$ and the scalar set ${\mathbb{D}}_{\circ ⟨k⟩}$. As the actuators of the PM are all on the motion branches, we will elucidate the method of establishing equation ${\text{E}}_{\text{A}}$ by discussing several common branches (Fig.5).

Fig.5 Several commonly used branches in PM design: (a) UPS branch driven by prismatic joint, (b) PUS branch driven by prismatic joint, and (c) RUS branch driven by revolute joint.

Full size|PPT slide

For Fig.5(a), ${\text{E}}_{\text{A}}$ can be established based on the dot product operation of vectors as follows:

(2)$\begin{array}{rl}{c}_{i⟨k-1⟩}{\mathit{p}}_{i⟨k-1⟩}\cdot {\mathit{q}}_{i⟨k⟩}=& c_{i⟨k-1⟩}\left|{\mathit{q}}_{i⟨k⟩}\right|\cos (\theta ({\mathit{p}}_{i⟨k-1⟩},{\mathit{q}}_{i⟨k⟩}))\\ =& c_{i⟨k-1⟩}(c_{i⟨k⟩}\cos (\theta ({\mathit{p}}_{i⟨k-1⟩},{\mathit{p}}_{i⟨k⟩}))-c_{i⟨k-1⟩}),\end{array}$

where the function $\theta ({\mathit{u}}_1,{\mathit{u}}_2)$ represents the included angle between vectors ${\mathit{u}}_1$ and ${\mathit{u}}_2$, and the vectors ${\mathit{p}}_{i⟨k⟩}$ and ${\mathit{r}}_i$ represent the unit axis vectors of the prismatic joint and the revolute joint, respectively. When the step number n is sufficiently large and the angle $\theta ({\mathit{p}}_{i⟨k-1⟩},{\mathit{p}}_{i⟨k⟩})$ tends to zero, Eq. (2) can be simplified through the Maclaurin series as follows:

(3)$c_{i⟨k-1⟩}{\mathit{p}}_{i⟨k-1⟩}\cdot {\mathit{q}}_{i⟨k⟩}=c_{i⟨k-1⟩}(c_{i⟨k⟩}-c_{i⟨k-1⟩})\Rightarrow {\mathit{p}}_{i⟨k-1⟩}\cdot {\mathit{q}}_{i⟨k⟩}=d_{i⟨k⟩}.$

Equation (3) is a linear equation with three unknowns. In fact, given that ${\text{E}}_{\text{A}}$ does not need to consider other joint constraints on the branch, this equation can be applied to most branches where the intermediate prismatic joint serves as the actuator, such as branches RPS, UPU, and RPU. Abbreviations like RPS are notations for serial joint branches. The definition of R/P/S/U can be found in Section 1.

Fig.5(b) and Fig.5(c) depict two scenarios where the actuator is fixed to the chassis. Here, ${\text{E}}_{\text{A}}$ cannot be directly established because the action of the actuator has been transferred to the next joint, resulting in the known pose of the next joint at each step containing parameter $d_{i⟨k⟩}$. Specifically, in both cases, the position of nodes $U_{i⟨k⟩}$ is always known, indicating that the influence of the actuator has been transferred to other constraint equations.

● Rigid body motion equations (${\text{E}}_{\text{M}}$)

In step k, the motion of any rigid body component can be equivalently represented as a linear combination of small translations and rotations. Thus, ${\mathbb{Q}}_{\circ ⟨k⟩}$ can be decomposed into a translation vector ${\mathit{t}}_{\circ ⟨k⟩}$ and a set of rotation vectors given by:

(4)${\mathbb{W}}_{\circ ⟨k⟩}=\left\{{\mathit{w}}_{1⟨k⟩},{\mathit{w}}_{2⟨k⟩},\dots ,{\mathit{w}}_{h⟨k⟩}\right\},$

where ${\mathit{w}}_{i⟨k⟩}={\mathit{q}}_{i⟨k⟩}-{\mathit{t}}_{\circ ⟨k⟩}$. Similarly, we provide three cases of end-effector configurations to illustrate the establishment process of ${\text{E}}_{\text{M}}$ (Fig.6).

Fig.6 Three fundamental rigid body motions of end-effectors: (a) pure translation motion, (b) pure rotation motion, and (c) composite motion.

Full size|PPT slide

The purpose of establishing ${\text{E}}_{\text{M}}$ is to describe the motion law of components and to possess good generality. In Fig.6(a), all ${\mathit{q}}_{i⟨k⟩}$ are equal to ${\mathit{t}}_{\circ ⟨k⟩}$, representing the simplest case. In Fig.6(b), only the rotational motion exists on the end-effector, and the following equation can be obtained through vector dot product operation:

(5)$\begin{array}{rl}\rightharpoonup O_{\text{D}}{P}_{i⟨k-1⟩}\cdot {\mathit{q}}_{i⟨k⟩}=& \left|\rightharpoonup O_{\text{D}}{P}_{i⟨k-1⟩}\right|\left|{\mathit{q}}_{i⟨k⟩}\right|\mathrm{c}\mathrm{o}\mathrm{s}\left(\theta \left(\rightharpoonup O_{\text{D}}{P}_{i⟨k-1⟩},{\mathit{q}}_{i⟨k⟩}\right)\right)\\ =& \left|\rightharpoonup O_{\text{D}}{P}_{i⟨k-1⟩}\right|\left(\left|\rightharpoonup O_{\text{D}}{P}_{i⟨k⟩}\right|\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{i,k}-\left|\rightharpoonup O_{\text{D}}{P}_{i⟨k-1⟩}\right|\right),\end{array}$

where ${\theta }_{i,k}=\theta \left(\rightharpoonup O_{\text{D}}{P}_{i⟨k-1⟩},\rightharpoonup O_{\text{D}}{P}_{i⟨k⟩}\right)$. Thus, we can linearize ${\text{E}}_{\text{M}}$ using the same simplification method as ${\text{E}}_{\text{A}}$, as follows:

(6)$\begin{array}{rl}\rightharpoonup O_{\text{D}}{P}_{i⟨k-1⟩}\cdot {\mathit{q}}_{i⟨k⟩}=& \left|\rightharpoonup O_{\text{D}}{P}_{i⟨k-1⟩}\right|\left(\left|\rightharpoonup O_{\text{D}}{P}_{i⟨k⟩}\right|-\left|\rightharpoonup O_{\text{D}}{P}_{i⟨k-1⟩}\right|\right)\\ \Rightarrow & \rightharpoonup O_{\text{D}}{P}_{i⟨k-1⟩}\cdot {\mathit{q}}_{i⟨k⟩}=0.\end{array}$

This type of ${\text{E}}_{\text{P}}$ can be used for components that are capable of pure rotation motion only, such as the end effector of a 3-UPS/S PM. For the more general scenario depicted in Fig.6(c), the vector ${\mathit{t}}_{\circ ⟨k⟩}$ can be defined as the translation vector of the origin $O_{\text{D}}$. This step enables the establishment of ${\text{E}}_{\text{M}}$ as follows:

(7)$\rightharpoonup O_{\text{D}⟨k-1⟩}{P}_{i⟨k-1⟩}\cdot {\mathit{w}}_{i⟨k⟩}=0\Rightarrow \rightharpoonup O_{\text{D}⟨k-1⟩}{P}_{i⟨k-1⟩}\cdot ({\mathit{q}}_{i⟨k⟩}-{\mathit{t}}_{i⟨k⟩})=0,$

where the origin $O_{\text{D}}$ also serves as a mark point, with $O_{\text{D}⟨k⟩}$ being its nodes.

Furthermore, Fig.5(b) and 5(c) subtly depict two other scenarios, although their principles still involve vector rotation motion. Taking Fig.5(b) as an example, ${\text{E}}_{\text{M}}$ can be established as follows:

(8)$\rightharpoonup U_{i⟨k-1⟩}{S}_{i⟨k-1⟩}\cdot ({\mathit{q}}_{i⟨k⟩}-d_{i⟨k⟩}{\mathit{p}}_i)=0.$

● Rigid body projection equations (${\text{E}}_{\text{P}}$)

In a rigid body, given that the distance between any two nodes in set ${\mathbb{P}}_{\circ ⟨k⟩}$ remains constant, the projection theorem can be applied to establish ${\text{E}}_{\text{P}}$ through dot products between microscopic vectors. The purpose of establishing ${\text{E}}_{\text{P}}$ is to describe the spatial relationship between two mark points on a component, thereby reflecting the characteristics of the rigid body. Therefore, ${\text{E}}_{\text{P}}$ can be established between any two component vectors, as shown in Fig.7, and in the following equations:

Fig.7 Projection relationship between two vectors on a rigid body during step motion.

Full size|PPT slide

(9)$\rightharpoonup O_{\text{D}⟨k⟩}{P}_{i⟨k⟩}\cdot \rightharpoonup O_{\text{D}⟨k⟩}{P}_{j⟨k⟩}=\rightharpoonup O_{\text{D}⟨0⟩}{P}_{i⟨0⟩}\cdot \rightharpoonup O_{\text{D}⟨0⟩}{P}_{j⟨0⟩},$

where

(10)$\rightharpoonup O_{\text{D}⟨k⟩}{P}_{i⟨k⟩}=\rightharpoonup O_{\text{D}⟨k-1⟩}{P}_{i⟨k-1⟩}+{\mathit{q}}_{i⟨k⟩}-{\mathit{t}}_{\circ ⟨k⟩}.$

Expanding Eq. (9) in accordance with Eq. (10) and directly eliminating non-linear terms based on infinitesimal theory, the linearized form of ${\text{E}}_{\text{P}}$ is obtained as follows:

(11)$\begin{array}{rl}& \rightharpoonup O_{\text{D}⟨k-1⟩}{P}_{j⟨k-1⟩}\cdot {\mathit{q}}_{i⟨k⟩}+\rightharpoonup O_{\text{D}⟨k-1⟩}{P}_{i⟨k-1⟩}\cdot {\mathit{q}}_{j⟨k⟩}\\ -& \left(\rightharpoonup O_{\text{D}⟨k-1⟩}{P}_{i⟨k-1⟩}+\rightharpoonup O_{\text{D}⟨k-1⟩}{P}_{j⟨k-1⟩}\right)\cdot {\mathit{t}}_{\circ ⟨k⟩}=0.\end{array}$

● Special equations (${\text{E}}_{\text{S}}$)

For certain PMs with specially designed features, leveraging their inherent geometric properties allows for a more expedient establishment of step equations, termed as ${\text{E}}_{\text{S}}$. Fig.7 provides two specific examples to illustrate the establishment of ${\text{E}}_{\text{S}}$.

For the scenario depicted in Fig.8(a), when the three mark points exhibit circularly distributed properties around the origin $O_{\text{D}}$, ${\text{E}}_{\text{S}}$ can be established as follows:

Fig.8 Two special cases often considered in PM design: (a) joints with circularly distributed design at the joint center, and (b) mutually perpendicular axial vectors ${\mathit{p}}_{i⟨k⟩}$ and ${\mathit{r}}_i$ in branch RPS.

Full size|PPT slide

(12)$\begin{array}{rl}\sum _{i=1}^3\rightharpoonup O_{\text{D}⟨k⟩}{P}_{i⟨k⟩}=& \sum _{i=1}^3\left(\rightharpoonup O_{\text{D}⟨k-1⟩}{P}_{i⟨k-1⟩}+{\mathit{q}}_{i⟨k⟩}-{\mathit{t}}_{\circ ⟨k⟩}\right)\\ =& \sum _{i=1}^3{\mathit{q}}_{i⟨k⟩}-3{\mathit{t}}_{\circ ⟨k⟩}={\mathbf{0}}_{3\times 1}.\end{array}$

For the branch RPS depicted in Fig.8(b), where the mark point $S_i$ can only rotate about the axis vector ${\mathit{r}}_i$, its motion space is confined to a plane. Based on this scenario, ${\text{E}}_{\text{S}}$ can thus be established as follows:

(13)${\mathit{r}}_i\cdot (c_{i⟨k⟩}{\mathit{p}}_{i⟨k⟩})={\mathit{r}}_i\cdot (c_{i⟨k-1⟩}{\mathit{p}}_{i⟨k-1⟩}+{\mathit{q}}_{i⟨k⟩})={\mathit{r}}_i\cdot {\mathit{q}}_{i⟨k⟩}=0.$

The above cases almost cover common situations. When the number of unknown coordinates equals the number of equations, a linear equation system can be established based on parameter set ${\mathbb{P}}_{\circ ⟨k-1⟩}\cup {\mathbb{D}}_{\circ ⟨k⟩}$ and independent variable set ${\mathbb{Q}}_{\circ ⟨k⟩}$, thereby obtaining a series of mappings depicted as follows:

(14)$\{\begin{array}{c}{\mathcal{L}}_{\circ ⟨k⟩}:\left[\begin{array}{cc}{\mathbb{P}}_{\circ ⟨k-1⟩}& {\mathbb{D}}_{\circ ⟨k⟩}\end{array}\right]\to {\mathbb{Q}}_{\text{C}⟨k⟩},\\ {\mathcal{L}}_{\circ ⟨k⟩}\mapsto {\mathit{L}}_{\circ ⟨k⟩}\cdot {\mathit{q}}_{\circ ⟨k⟩}={\mathit{b}}_{\circ ⟨k⟩}.\end{array}$

Notably, the number of available step equations often exceeds the number of unknowns, and some step equations impose the same constraints on the step motion, leading to linear dependence between them. Therefore, the priority of selecting the four types of step equations is different. ${\text{E}}_{\text{A}}$ is necessary, as is ${\text{E}}_{\text{M}}$, which contains actuator parameters. On this basis, ${\text{E}}_{\text{S}}$ should be prioritized because it does not introduce errors from eliminating non-linear terms, thus providing higher accuracy. Finally, ${\text{E}}_{\text{M}}$ should be used to complete the equation set, while ${\text{E}}_{\text{P}}$ is the least considered due to its higher error margin.

In the mappings mentioned above, the ${\mathcal{L}}_{\circ ⟨k⟩}$ represents a linear mapping within the interval $T_{\circ ⟨k⟩}$. The difference between ${\mathbb{Q}}_{\text{C}⟨k⟩}$ and ${\mathbb{Q}}_{\circ ⟨k⟩}$ lies in the calculation error generated by eliminating the non-linear factors. Within the algebraic structure of mapping ${\mathcal{L}}_{\circ ⟨k⟩}$, ${\mathit{L}}_{\circ ⟨k⟩}$ represents the step matrix. The vector ${\mathit{q}}_{\circ ⟨k⟩}$ consists of all coordinates of ${\mathit{q}}_{i⟨k⟩}$ in frame $\mathit{P}$, which serves as the integration variable of the system. Vector ${\mathit{b}}_{\circ ⟨k⟩}$ consists of all constant terms.

Research indicates that when the step number n is sufficiently large, the error between ${\mathbb{Q}}_{\text{C}⟨k⟩}$ and ${\mathbb{Q}}_{\circ ⟨k⟩}$ can be reduced to negligible levels while maintaining the computational efficiency of the system. Therefore, ${\mathbb{Q}}_{\text{C}⟨k⟩}$ can be used instead of ${\mathbb{Q}}_{\circ ⟨k⟩}$ to solve ${\mathbb{P}}_{\circ ⟨k⟩}$, with the mapping operation as follows:

(15)$\{\begin{array}{l}{{\mathcal{P}}_{\circ ⟨k⟩}}_{}:\left[\begin{array}{cc}{\mathbb{P}}_{\circ ⟨k-1⟩}& {\mathbb{Q}}_{\circ ⟨k⟩}\end{array}\right]\to {\mathbb{P}}_{\circ ⟨k⟩},\\ {{\mathcal{P}}_{\circ ⟨k⟩}}_{}\mapsto \rightharpoonup O{P}_{i⟨k-1⟩}+{\mathit{q}}_{i⟨k⟩}=\rightharpoonup O{P}_{i⟨k⟩}.\end{array}$

Therefore, a significant step operation within the interval $T_{\circ ⟨k⟩}$ can be denoted as step ${\mathcal{S}}_{\circ ⟨k⟩}$, as shown in Fig.9. The $det\left({\mathit{L}}_{\circ ⟨k⟩}\right)$ represents the determinant of matrix ${\mathit{L}}_{\circ ⟨k⟩}$, and the determination of singularity will be elaborated upon later. At this point, the principles and modeling logic of the FSI method have been thoroughly explained.

Fig.9 Operational procedure of the FSI method and the transfer relationship between various datasets.

Full size|PPT slide

4 4 Singularity detection and discrimination

The computational process of the FSI method resembles the motion trajectory tracking of PMs. Therefore, considering the singularities during a PM’s motion ensures the smooth tracking of its trajectory. Singularities are a special class of configurations featuring abnormal stiffness or degree of freedom (DOF) in PMs, which may have negative effect on PM control [29‒31]. Thus, avoiding singularities is an important issue in robot design and control. Cheng et al. [12] investigated the relationship between the singularities of PMs and actuator parameters and determined the algebraic form of singularity boundaries through geometric analysis. Wu and Bai [31] proposed a method based on analytic geometry that finds shape singularities in parallel platforms, thus avoiding large-scale computations, and provided three configuration examples for illustration. Merlet [32] proposed an algorithm that combines formal and numerical computations to detect singularities or near-singularities in arbitrary workspaces or trajectories. Ben-Horin and Shoham [33] and Kanaan et al. [34] analyzed the singularity conditions of PMs using Grassmann–Cayley algebra. Yang and Li [35] established a mathematical model of singularity boundaries using differential manifold theory.

For the FSI method, this section begins by discussing the issue of non-existent or multiple solutions in linear equations, which may occur during the computational process. Through analysis, we discovered that different forms of solutions correspond to different states of the PM. When the equations have a unique solution, the motion of the PM can be determined; this scenario is the most common. If unique solutions are obtained for all steps, then the spatiotemporal curve set ${\mathbb{L}}_{\text{PM}}$ possesses uniqueness. Therefore, the FSI method does not require additional filtering algorithms to select the desired solution from numerous possibilities. When the determinant of the step matrix approaches zero (which is almost impossible to occur exactly in floating-point arithmetic systems), the motion characteristics of the PM tend to exhibit singularity.

From the perspective of vector space, the following two situations occur when the determinant of the step matrix is zero (where $\left[\begin{array}{cc}{\mathit{L}}_{\circ ⟨k⟩}& {\mathit{b}}_{\circ ⟨k⟩}\end{array}\right]$ is the augmented matrix on step ${\mathcal{S}}_{\circ ⟨k⟩}$):

● $\text{det}\left({\mathit{L}}_{\circ ⟨k⟩}\right)=0\text{\&\&}\text{rank}\left({\mathit{L}}_{\circ ⟨k⟩}\right)=\text{rank}\left(\left[\begin{array}{cc}{\mathit{L}}_{\circ ⟨k⟩}& {\mathit{b}}_{\circ ⟨k⟩}\end{array}\right]\right)$

This situation indicates that in step $k,$ the solution of the equation is not unique, that is, the step motion of the PM is uncertain. Thus, at least one coordinate of vector ${\mathit{q}}_{\circ ⟨k⟩}$ cannot be determined with certainty, implying that the PM loses one or more constraint capability on the end effector.

● $\text{det}\left({\mathit{L}}_{\circ ⟨k⟩}\right)=0\text{\&\&}\text{rank}\left({\mathit{L}}_{\circ ⟨k⟩}\right)<\text{rank}\left(\left[\begin{array}{cc}{\mathit{L}}_{\circ ⟨k⟩}& {\mathit{b}}_{\circ ⟨k⟩}\end{array}\right]\right)$

This situation indicates that in step $k,$ the equation has no solution, implying that the PM lacks a direction for motion. Consequently, finding any vector ${\mathit{q}}_{\circ ⟨k⟩}$ to accomplish the step motion is impossible, meaning that the PM loses one or more DOFs on the end effector.

Considering the characteristics of floating-point arithmetic systems, we design and introduce a parameter denoted as $\epsilon$ to measure the proximity between the PM configuration and singularity. The idea of handling singularity has also been employed in Refs. [36,37], and setting a threshold is particularly effective in FSI method computations. This parameter serves as a threshold to preemptively halt the PM’s motion toward singularity (Fig.10). Before each computation of mapping ${\mathcal{L}}_{\circ ⟨k⟩}$, an additional step to calculate $\left|\text{det}\left({\mathit{L}}_{\circ ⟨k⟩}\right)\right|$ is introduced. $\left|\text{det}\left({\mathit{L}}_{\circ ⟨k⟩}\right)\right|$ values exceeding threshold $\epsilon$ indicate that the PM is still far from singularity. Conversely, $\left|\text{det}\left({\mathit{L}}_{\circ ⟨k⟩}\right)\right|$ values below threshold $\epsilon$ suggest that the PM has quasi-singularity, and the computation is terminated.

Fig.10 (a) Principle of singularity detection and (b) handling in the FSI method.

Full size|PPT slide

The magnitude of threshold $\epsilon$ determines the width of the caution zone; therefore, selecting an appropriate threshold is crucial for the FSI method. Using extensive computational data and function fitting, we provide the following empirical formula for threshold $\epsilon$ as reference:

(16)$\epsilon =0.08\text{exp}\left(-\frac{\sqrt{2}n}{2000}\right)\cdot \left|det\left({\mathit{L}}_{\circ ⟨1⟩}\right)\right|.$

5 5 Error evaluation model and case studies

In this section, we first establish an error evaluation model for assessing the final computational results of the FSI method. Using this model and the theoretical methods discussed earlier, we illustrate the specific modeling process and computational results using generalized 6-UPS and 3-UPS/S configurations as examples.

5.1 5.1 Error evaluation model

Directly comparing the errors between ${\mathbb{P}}_{\circ ⟨\text{E}⟩}$ and ${\mathbb{P}}_{\circ ⟨n⟩}$ cannot effectively illustrate the overall error of the end effector. Therefore, we first define the matrix of coordinate system ${\mathit{P}}_{\text{D}}$ and similarly define the matrix of ${\mathit{P}}_{\text{C}}$ in Eq. (17):

(17)$\{\begin{array}{c}{\mathit{P}}_{\text{D}}=\left[\begin{array}{cccc}\begin{array}{c}{\mathit{x}}_{\text{D}}\\ 0\end{array}& \begin{array}{c}{\mathit{y}}_{\text{D}}\\ 0\end{array}& \begin{array}{c}{\mathit{z}}_{\text{D}}\\ 0\end{array}& \begin{array}{c}{\mathit{o}}_{\text{D}}\\ 1\end{array}\end{array}\right]=\left[\begin{array}{cccc}{\mathit{\lambda }}_{\text{D}X}& {\mathit{\lambda }}_{\text{D}Y}& {\mathit{\lambda }}_{\text{D}Z}& {\mathit{\lambda }}_{\text{DO}}\end{array}\right],\\ {\mathit{P}}_{\text{C}}=\left[\begin{array}{cccc}\begin{array}{c}{\mathit{x}}_{\text{C}}\\ 0\end{array}& \begin{array}{c}{\mathit{y}}_{\text{C}}\\ 0\end{array}& \begin{array}{c}{\mathit{z}}_{\text{C}}\\ 0\end{array}& \begin{array}{c}{\mathit{o}}_{\text{C}}\\ 1\end{array}\end{array}\right]=\left[\begin{array}{cccc}{\mathit{\lambda }}_{\text{C}X}& {\mathit{\lambda }}_{\text{C}Y}& {\mathit{\lambda }}_{\text{C}Z}& {\mathit{\lambda }}_{\text{CO}}\end{array}\right],\end{array}$

where the matrix ${\mathit{P}}_{\text{D}}$ and coordinate system ${\mathit{P}}_{\text{D}}$ in Tab.1 represent two forms of the same element and are used for the computation and description of the coordinate system, respectively; ${\mathit{x}}_{\text{D}}$, ${\mathit{y}}_{\text{D}}$, and ${\mathit{z}}_{\text{D}}$ denote the three unit axes of ${\mathit{P}}_{\text{D}}$; matrix ${\mathit{P}}_{\text{C}}$ is the pose matrix of the end-effector obtained through computation, with ${\mathit{x}}_{\text{C}}$, ${\mathit{y}}_{\text{C}}$, and ${\mathit{z}}_{\text{C}}$ representing its three unit axes; ${\mathit{o}}_{\text{D}}=\rightharpoonup O{O}_{\text{D}}$, ${\mathit{o}}_{\text{C}}=\rightharpoonup O{O}_{\text{C}}$, and point $O_{\text{C}}$ is the computed $O_{\text{D}}$ with associated errors; ${\mathit{\lambda }}_{\text{D}X}$, ${\mathit{\lambda }}_{\text{D}Y}$, ${\mathit{\lambda }}_{\text{D}Z}$, and ${\mathit{\lambda }}_{\text{DO}}$ refer to the column vectors of matrix ${\mathit{P}}_{\text{D}}$; ${\mathit{\lambda }}_{\text{C}X}$, ${\mathit{\lambda }}_{\text{C}Y}$, ${\mathit{\lambda }}_{\text{C}Z}$, and ${\mathit{\lambda }}_{\text{CO}}$ refer to the column vectors of matrix ${\mathit{P}}_{\text{C}}$; and $X$, $Y$, and $Z$ denote the identifiers for the three coordinate directions.

Given that ${\mathit{P}}_{\text{C}}$ is determined by ${\mathbb{P}}_{\circ ⟨n⟩}$, the error between ${\mathit{P}}_{\text{C}}$ and ${\mathit{P}}_{\text{D}}$ is representative. Furthermore, we define the error evaluation metrics as follows:

(18)$\left[\begin{array}{c}{\delta }_X\\ {\delta }_Y\\ {\delta }_Z\\ {\mathbf{\Delta }}_{\text{O}}\end{array}\right]=\left[\begin{array}{c}{\delta }_X\\ {\delta }_Y\\ {\delta }_Z\\ {\mathrm{\Delta }}_X\\ {\mathrm{\Delta }}_Y\\ {\mathrm{\Delta }}_Z\end{array}\right]=\left[\begin{array}{c}\theta \left({\mathit{\lambda }}_{\text{D}X},{\mathit{\lambda }}_{\text{C}X}\right)\\ \theta \left({\mathit{\lambda }}_{\text{D}Y},{\mathit{\lambda }}_{\text{C}Y}\right)\\ \theta \left({\mathit{\lambda }}_{\text{D}Z},{\mathit{\lambda }}_{\text{C}Z}\right)\\ {\mathit{\lambda }}_{\text{CO}}-{\mathit{\lambda }}_{\text{DO}}\end{array}\right],$

where ${\delta }_X$, ${\delta }_Y$, and ${\delta }_Z$ are the posture error metrics; ${\mathrm{\Delta }}_X$, ${\mathrm{\Delta }}_Y$, and ${\mathrm{\Delta }}_Z$ are the position error metrics, and ${\mathbf{\Delta }}_{\mathrm{O}}={\left[\begin{array}{ccc}{\mathrm{\Delta }}_X& {\mathrm{\Delta }}_Y& {\mathrm{\Delta }}_Z\end{array}\right]}^{\text{T}}$. The establishment of this error model is based on geometric principles as detailed in Fig.11.

Fig.11 (a) Explanation of the geometric principles and (b) algebraic structure of the error evaluation model.

Full size|PPT slide

5.2 5.2 Example of 6-UPS PM

The 6-UPS PM (Fig.12) is a 6-DOF configuration with the prismatic joint as the actuator. We define the geometric model of this configuration as follows:

Fig.12 Generalized 6-UPS PM: (a) topology structure and (b) relationship topology diagram of equation ${\text{E}}_{\text{P}}$.

Full size|PPT slide

1) In the initial pose, frame ${\mathit{P}}_{\text{D}}$ completely coincides with frame $\mathit{P}$;

2) Joint ${\text{U}}_i$ is fixed on frame ${\mathit{P}}_{\text{D}}$, and joint ${\text{S}}_i$ is fixed on frame $\mathit{P}$;

3) The node set is defined as ${\mathbb{P}}_{\circ ⟨k⟩}=\left\{O_{\text{D}⟨k⟩},S_{1⟨k⟩},S_{2⟨k⟩},\dots ,S_{6⟨k⟩}\right\}$.

According to Eqs. (3), (7), and (11), the mapping ${\mathcal{L}}_{\circ ⟨k⟩}$ for the 6-UPS configuration can be established as shown as Eq. (19), where ${\mathit{s}}_{i⟨k⟩}=\rightharpoonup O_{\text{D}⟨k⟩}{S}_{i⟨k⟩}$.

Equation (19) represents an important mapping of approximate vectors for all motions in analytical step k, with the parameters iteratively updated as the steps progress. The parameters in matrix ${\mathit{L}}_{\circ ⟨k⟩}$ and vector ${\mathit{b}}_{\circ ⟨k⟩}$ are known, making this a standard linear algebra problem. The mathematical foundations involved in the FSI method are easily translatable into program code, such as the pseudocode Algorithm 1 (Tab.2) in this example.

Tab.2 Pseudocode of the Algorithm 1

	Algorithm 1 FK calculation of 6-UPS PM based on the FSI method
	Input Set ${\mathbb{C}}_{\circ ⟨E⟩}$.
	Output Matrix ${\mathit{P}}_{\text{C}}$.
1	Set the initial node set ${\mathbb{P}}_{\circ ⟨0⟩}$ and compute ${\mathbb{C}}_{\circ ⟨0⟩}$;
2	Set the value of n (number of steps);
3	Design the set $\mathbb{F}$, and if there are no specific requirements for the motion path, it can be set as $f_i(t)=(c_{i⟨\text{E}⟩}-c_{i⟨0⟩})t+c_{i⟨0⟩}$, where $t\in T=\left[0,1\right]$;
4	for $k=1:1:n$
5	Compute or iterate ${\mathbb{D}}_{\circ ⟨k⟩}$, ${\mathit{s}}_{i⟨k-1⟩}$, and ${\mathit{p}}_{i⟨k-1⟩}$;
6	Establish the step matrix ${\mathit{L}}_{\circ ⟨k⟩}$ and vector ${\mathit{b}}_{\circ ⟨k⟩}$;
7	Compute $\left|det\left({\mathit{L}}_{\circ ⟨k⟩}\right)\right|$;
8	if $k=1$
9	Set the threshold ε;
10	else if $\left|det\left({\mathit{L}}_{\circ ⟨k⟩}\right)\right|<\mathit{\text{ε}}$
11	return “Quasi-singular configuration”;
12	break
13	end
14	Compute ${\mathit{q}}_{\circ ⟨k⟩}={\left({\mathit{L}}_{\circ ⟨k⟩}\right)}^{-1}{\mathit{b}}_{\circ ⟨k⟩}$;
15	Compute and iterate the set ${\mathbb{P}}_{\circ ⟨k⟩}$;
16	end
17	Compute the matrix ${\mathit{P}}_{\text{C}}$ based on set ${\mathbb{P}}_{\circ ⟨n⟩}$;
18	return ${\mathit{P}}_{\text{C}}$.

(19)$\left[\begin{array}{ccccccccc}{\text{E}}_{\text{A1}}\text{:}& |& {\mathit{p}}_{1⟨k-1⟩}^{\text{T}}& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}\\ {\text{E}}_{\text{A2}}\text{:}& |& {\mathbf{0}}_{1\times 3}& {\mathit{p}}_{2⟨k-1⟩}^{\text{T}}& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}\\ {\text{E}}_{\text{A3}}\text{:}& |& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}& {\mathit{p}}_{3⟨k-1⟩}^{\text{T}}& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}\\ {\text{E}}_{\text{A4}}\text{:}& |& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}& {\mathit{p}}_{4⟨k-1⟩}^{\text{T}}& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}\\ {\text{E}}_{\text{A5}}\text{:}& |& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}& {\mathit{p}}_{5⟨k-1⟩}^{\text{T}}& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}\\ {\text{E}}_{\text{A6}}\text{:}& |& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}& {\mathit{p}}_{6⟨k-1⟩}^{\text{T}}& {\mathbf{0}}_{1\times 3}\\ \\ {\text{E}}_{\text{M1}}\text{:}& |& {\mathit{s}}_{1⟨k-1⟩}^{\text{T}}& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}& -{\mathit{s}}_{1⟨k-1⟩}^{\text{T}}\\ {\text{E}}_{\text{M2}}\text{:}& |& {\mathbf{0}}_{1\times 3}& {\mathit{s}}_{2⟨k-1⟩}^{\text{T}}& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}& -{\mathit{s}}_{2⟨k-1⟩}^{\text{T}}\\ {\text{E}}_{\text{M3}}\text{:}& |& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}& {\mathit{s}}_{3⟨k-1⟩}^{\text{T}}& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}& -{\mathit{s}}_{3⟨k-1⟩}^{\text{T}}\\ {\text{E}}_{\text{M4}}\text{:}& |& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}& {\mathit{s}}_{4⟨k-1⟩}^{\text{T}}& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}& -{\mathit{s}}_{4⟨k-1⟩}^{\text{T}}\\ {\text{E}}_{\text{M5}}\text{:}& |& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}& {\mathit{s}}_{5⟨k-1⟩}^{\text{T}}& {\mathbf{0}}_{1\times 3}& -{\mathit{s}}_{5⟨k-1⟩}^{\text{T}}\\ {\text{E}}_{\text{M6}}\text{:}& |& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}& {\mathit{s}}_{6⟨k-1⟩}^{\text{T}}& -{\mathit{s}}_{6⟨k-1⟩}^{\text{T}}\\ \\ {\text{E}}_{\text{P1}}\text{:}& |& {\mathit{s}}_{2⟨k-1⟩}^{\text{T}}& {\mathit{s}}_{1⟨k-1⟩}^{\text{T}}& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}& -({\mathit{s}}_{1⟨k-1⟩}^{\text{T}}+{\mathit{s}}_{2⟨k-1⟩}^{\text{T}})\\ {\text{E}}_{\text{P2}}\text{:}& |& {\mathit{s}}_{3⟨k-1⟩}^{\text{T}}& {\mathbf{0}}_{1\times 3}& {\mathit{s}}_{1⟨k-1⟩}^{\text{T}}& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}& -({\mathit{s}}_{1⟨k-1⟩}^{\text{T}}+{\mathit{s}}_{3⟨k-1⟩}^{\text{T}})\\ {\text{E}}_{\text{P3}}\text{:}& |& {\mathit{s}}_{4⟨k-1⟩}^{\text{T}}& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}& {\mathit{s}}_{1⟨k-1⟩}^{\text{T}}& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}& -({\mathit{s}}_{1⟨k-1⟩}^{\text{T}}+{\mathit{s}}_{4⟨k-1⟩}^{\text{T}})\\ {\text{E}}_{\text{P4}}\text{:}& |& {\mathit{s}}_{5⟨k-1⟩}^{\text{T}}& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}& {\mathit{s}}_{1⟨k-1⟩}^{\text{T}}& {\mathbf{0}}_{1\times 3}& -({\mathit{s}}_{1⟨k-1⟩}^{\text{T}}+{\mathit{s}}_{5⟨k-1⟩}^{\text{T}})\\ {\text{E}}_{\text{P5}}\text{:}& |& {\mathit{s}}_{6⟨k-1⟩}^{\text{T}}& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}& {\mathit{s}}_{1⟨k-1⟩}^{\text{T}}& -({\mathit{s}}_{1⟨k-1⟩}^{\text{T}}+{\mathit{s}}_{6⟨k-1⟩}^{\text{T}})\\ {\text{E}}_{\text{P6}}\text{:}& |& {\mathbf{0}}_{1\times 3}& {\mathit{s}}_{3⟨k-1⟩}^{\text{T}}& {\mathit{s}}_{2⟨k-1⟩}^{\text{T}}& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}& -({\mathit{s}}_{2⟨k-1⟩}^{\text{T}}+{\mathit{s}}_{3⟨k-1⟩}^{\text{T}})\\ {\text{E}}_{\text{P7}}\text{:}& |& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}& {\mathit{s}}_{4⟨k-1⟩}^{\text{T}}& {\mathit{s}}_{3⟨k-1⟩}^{\text{T}}& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}& -({\mathit{s}}_{3⟨k-1⟩}^{\text{T}}+{\mathit{s}}_{4⟨k-1⟩}^{\text{T}})\\ {\text{E}}_{\text{P8}}\text{:}& |& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}& {\mathit{s}}_{5⟨k-1⟩}^{\text{T}}& {\mathit{s}}_{4⟨k-1⟩}^{\text{T}}& {\mathbf{0}}_{1\times 3}& -({\mathit{s}}_{4⟨k-1⟩}^{\text{T}}+{\mathit{s}}_{5⟨k-1⟩}^{\text{T}})\\ {\text{E}}_{\text{P9}}\text{:}& |& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}& {\mathit{s}}_{6⟨k-1⟩}^{\text{T}}& {\mathit{s}}_{5⟨k-1⟩}^{\text{T}}& -({\mathit{s}}_{5⟨k-1⟩}^{\text{T}}+{\mathit{s}}_{6⟨k-1⟩}^{\text{T}})\end{array}\right]\left[\begin{array}{c}{\mathit{q}}_{1⟨k⟩}\\ {\mathit{q}}_{2⟨k⟩}\\ {\mathit{q}}_{3⟨k⟩}\\ {\mathit{q}}_{4⟨k⟩}\\ {\mathit{q}}_{5⟨k⟩}\\ {\mathit{q}}_{6⟨k⟩}\\ {\mathit{t}}_{\circ ⟨k⟩}\end{array}\right]=\left[\begin{array}{c}{d}_{1⟨k⟩}\\ d_{2⟨k⟩}\\ d_{3⟨k⟩}\\ d_{4⟨k⟩}\\ d_{5⟨k⟩}\\ d_{6⟨k⟩}\\ {\mathbf{0}}_{15\times 1}\end{array}\right],$

A specific numerical example is provided below to evaluate the computational performance of the FSI method.

The initial node set ${\mathbb{P}}_{\circ ⟨0⟩}$ comprises the following: $U_1=\left(130,400,125\right)$, $U_2=\left(170,395,55\right)$, $U_3=\left(40,402,-180\right)$, $U_4=\left(-35,398,-175\right)$, $U_5=\left(172,405,60\right)$, and $U_6=\left(-135,396,124\right)$; $S_{1⟨0⟩}=\left(30,2,135\right)$, $S_{2⟨0⟩}=\left(132,-3,-45\right)$, $S_{3⟨0⟩}=\left(102,1,-95\right)$, $S_{4⟨0⟩}=\left(-105,-4,-90\right)$, $S_{5⟨0⟩}=\left(-130,5,-42\right)$, and $S_{6⟨0⟩}=\left(-28,-2,134\right)$.

${\mathbb{C}}_{\circ ⟨0⟩}$ can then be computed as follows:

(20)${\mathbb{C}}_{\circ ⟨0⟩}=\left\{\left|\rightharpoonup U_1{S}_{1⟨0⟩}\right|,\left|\rightharpoonup U_2{S}_{2⟨0⟩}\right|,\dots ,\left|\rightharpoonup U_6{S}_{6⟨0⟩}\right|\right\}.$

Set $\mathbb{F}$ can be designed following the method outlined in the pseudocode. To obtain comparative data for assessing the computational performance, we directly generate a random series of ${\mathit{P}}_{\text{D}⟨\text{E}⟩}$ and then compute $c_{i⟨\text{E}⟩}$ as follows:

(21)$c_{i⟨\text{E}⟩}=\left|{\mathit{P}}_{\text{D}⟨\text{E}⟩}\left[\begin{array}{c}{\mathit{s}}_{i⟨0⟩}\\ 1\end{array}\right]-\left[\begin{array}{c}\rightharpoonup O{U}_i\\ 1\end{array}\right]\right|.$

We perform the corresponding programming and calculations using MATLAB and organize the final computational results into Fig.13.

Fig.13 Computational experiments of the FSI method on the 6-UPS configuration: (a) step–computation time curve, (b) step–posture accuracy curves, and (c) step–position accuracy curves.

Full size|PPT slide

The results of the computational experiments indicate that the FSI method aligns with the initial hypothesis, that is, its computational error decreases exponentially with an increase in the number of steps. In addition, the computational time exhibits a linear growth trend with the increment of steps. In practical applications, balancing between error and computational speed can leverage the value of the FSI method. For instance, if the requirement for computational accuracy is to achieve $0^{\circ }0\text{'}0.36"$ and 0.2 μm, we can adjust step $n$ to 4500, resulting in a computation time of approximately 0.5 s. Although the FSI method may not compete with methods such as Newton method [5,38] purely in terms of computational speed, its advantages lie primarily in its systematic modeling logic and excellent modeling efficiency. Therefore, this level of computational efficiency is acceptable.

5.3 5.3 Example of the 3-UPS/S PM

The 3-UPS/S PM (Fig.14) is a 3-DOF configuration (pure point rotation DOF) with the prismatic joint as the actuator. We define the geometric model of this configuration as follows:

Fig.14 Generalized 3-UPS/S PM.

Full size|PPT slide

1) In the initial pose, frame ${\mathit{P}}_{\text{D}}$ completely coincides with frame $\mathit{P}$;

2) Joint ${\text{U}}_i$ is fixed on frame ${\mathit{P}}_{\text{D}}$, and joint ${\text{S}}_i$ is fixed on frame $\mathit{P}$;

3) The node set is defined as ${\mathbb{P}}_{\circ ⟨k⟩}=\left\{S_{1⟨k⟩},S_{2⟨k⟩},S_{3⟨k⟩}\right\}$.

Based on Eqs. (3), (6), and (11), the mapping ${\mathcal{L}}_{\circ ⟨k⟩}$ for the 3-UPS/S configuration can be established as detailed in Eq. (22). The pseudocode and computational results of this case are similar to those of the 6-UPS configuration; therefore, the relevant content is not reiterated here.

(22)$\left[\begin{array}{ccccc}{\text{E}}_{\text{A1}}\text{:}& |& {\mathit{p}}_{1⟨k-1⟩}^{\text{T}}& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}\\ {\text{E}}_{\text{A2}}\text{:}& |& {\mathbf{0}}_{1\times 3}& {\mathit{p}}_{2⟨k-1⟩}^{\text{T}}& {\mathbf{0}}_{1\times 3}\\ {\text{E}}_{\text{A3}}\text{:}& |& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}& {\mathit{p}}_{3⟨k-1⟩}^{\text{T}}\\ \\ {\text{E}}_{\text{M1}}\text{:}& |& {\mathit{s}}_{1⟨k-1⟩}^{\text{T}}& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}\\ {\text{E}}_{\text{M2}}\text{:}& |& {\mathbf{0}}_{1\times 3}& {\mathit{s}}_{2⟨k-1⟩}^{\text{T}}& {\mathbf{0}}_{1\times 3}\\ {\text{E}}_{\text{M3}}\text{:}& |& {\mathbf{0}}_{1\times 3}& {\mathbf{0}}_{1\times 3}& {\mathit{s}}_{3⟨k-1⟩}^{\text{T}}\\ \\ {\text{E}}_{\text{P1}}\text{:}& |& {\mathit{s}}_{2⟨k-1⟩}^{\text{T}}& {\mathit{s}}_{1⟨k-1⟩}^{\text{T}}& {\mathbf{0}}_{1\times 3}\\ {\text{E}}_{\text{P2}}\text{:}& |& {\mathbf{0}}_{1\times 3}& {\mathit{s}}_{3⟨k-1⟩}^{\text{T}}& {\mathit{s}}_{2⟨k-1⟩}^{\text{T}}\\ {\text{E}}_{\text{P3}}\text{:}& |& {\mathit{s}}_{3⟨k-1⟩}^{\text{T}}& {\mathbf{0}}_{1\times 3}& {\mathit{s}}_{1⟨k-1⟩}^{\text{T}}\end{array}\right]\left[\begin{array}{c}{\mathit{q}}_{1⟨k⟩}\\ {\mathit{q}}_{2⟨k⟩}\\ {\mathit{q}}_{3⟨k⟩}\end{array}\right]=\left[\begin{array}{c}{d}_{1⟨k⟩}\\ d_{2⟨k⟩}\\ d_{3⟨k⟩}\\ {\mathbf{0}}_{6\times 1}\end{array}\right].$

6 6 Comparative computational experiments

This section further illustrates the advantages and limitations of the FSI method through comparative computational experiments. Many methods have been developed to address the FK problem of PMs, and most of them, such as interval analysis and homotopy continuation methods, require numerical algorithms to achieve the final solution. However, modeling the FK problem itself makes it difficult to compare the merits of various methods. The FSI method provides a relatively comprehensive solution based on algebraic structures and logic, thereby greatly simplifying the difficulty of modeling and programming. This inherent practicality and advantage make it comparable with intelligent optimization algorithms based on target search, such as genetic algorithms, differential evolution, simulated annealing, and particle swarm optimization (PSO), in terms of modeling generality.

Known for its computational efficiency, we select the PSO algorithm for comparison and use the numerical example of the 6-UPS PM as the research subject to compare the numerical mapping relationship between the computational accuracy and computational time of the two methods. The essence of the PSO algorithm lies in designing the fitness function and determining the size of the particle swarm. The solution space of the particle swarm corresponds to the 6D pose parameters of the end-effector; therefore, the difference in the end-effector’s state should be used as the evaluation criterion for fitness. The fitness function $f_{\text{F}}$ used to quantify the accuracy of results can be designed as follows:

(23)$f_{\text{F}}=f_{\text{RMS}}\left(\mathrm{\Delta }{c}_1,\mathrm{\Delta }{c}_2,\dots ,\mathrm{\Delta }{c}_6\right),$

where $f_{\text{RMS}}$ is the root mean square of all $\mathrm{\Delta }{c}_i$, and $\mathrm{\Delta }{c}_i$ is the motion parameter difference of the ith actuator. The size of the particle swarm may affect the robustness and efficiency of the PSO algorithm. If the swarm size is too small, then the algorithm may not converge; conversely, it may lead to inefficiency.

The comparison of the two algorithms is conducted in MATLAB, and the results are depicted in Fig.15. Through this comparison, we can discern the strengths and weaknesses of the FSI method relative to the PSO method. In terms of overall error reduction rate, the PSO method is faster but the FSI method is smoother. Within the commonly used error interval $[{10}^{-4},{10}^{-3}]$, the FSI method exhibits higher computational efficiency. However, the FSI method shows the disadvantage of not achieving high precision limits, that is, it is not suitable when high computational accuracy is required. In addition, the computation of the PSO method relies on the inverse kinematics equations of the PM, which are not applicable when dealing with configurations where inverse kinematic conditions are difficult to determine (such as the 3-RPS mechanism). By contrast, the FSI method still excels in its superior versatility.

Fig.15 Time–error curves computed by the FSI and PSO methods: (a) time–posture error curves and (b) time–position error curves.

Full size|PPT slide

Overall, the FSI method’s computational accuracy and efficiency meet the basic requirements in most engineering applications. Moreover, it is a practical and excellent FK algorithm as it ensures the avoidance of singular poses without the need to consider robustness and multi-solution filtering issues.

7 7 Application of motion tracking analysis

The computational process of the FSI method resembles the motion trajectory tracking of the PM. This method can resolve all information regarding the entire motion of the PM. Given that the motion process of the PM is divided into $n$ segments by the FSI method, ${\mathit{P}}_{\text{C}⟨k⟩}$ can be calculated in real-time at each step as follows:

(24)${\mathit{P}}_{\text{C}⟨k⟩}=\left[\begin{array}{cc}{\mathit{R}}_{\text{C}⟨k⟩}& {\mathit{o}}_{\text{C}⟨k⟩}\\ 0& 1\end{array}\right]=\left[\begin{array}{cccc}\begin{array}{c}{\mathit{x}}_{\text{C}⟨k⟩}\\ 0\end{array}& \begin{array}{c}{\mathit{y}}_{\text{C}⟨k⟩}\\ 0\end{array}& \begin{array}{c}{\mathit{z}}_{\text{C}⟨k⟩}\\ 0\end{array}& \begin{array}{c}{\mathit{o}}_{\text{C}⟨k⟩}\\ 1\end{array}\end{array}\right],$

where ${\mathit{R}}_{\text{C}⟨k⟩}=\left[\begin{array}{ccc}{\mathit{x}}_{\text{C}⟨k⟩}& {\mathit{y}}_{\text{C}⟨k⟩}& {\mathit{z}}_{\text{C}⟨k⟩}\end{array}\right]$.

Hence, we can directly define set ${\mathbb{T}}_{\text{C}}$ to represent the motion trajectory of the end-effector as

(25)${\mathbb{T}}_{\text{C}}=\left\{{\mathit{P}}_{\text{C}⟨0⟩},{\mathit{P}}_{\text{C}⟨1⟩},{\mathit{P}}_{\text{C}⟨2⟩},\dots ,{\mathit{P}}_{\text{C}⟨n⟩}\right\}.$

When the number of steps is sufficient, the step equation approaches instantaneous motion. Thus, the average velocity of the step motion can be used instead of the step velocity. Therefore, the angular velocity of ${\mathit{P}}_{\text{C}}$ can be calculated as

(26)$\begin{array}{rl}{\left({\text{ω}}_{\text{C}⟨k⟩}\right)}^{\wedge }=& {\left[\begin{array}{c}{\omega }_{\text{C}X⟨k⟩}\\ {\omega }_{\text{C}Y⟨k⟩}\\ {\omega }_{\text{C}Z⟨k⟩}\end{array}\right]}^{\wedge }=\left[\begin{array}{ccc}0& -{\omega }_{\text{C}Z⟨k⟩}& {\omega }_{\text{C}Y⟨k⟩}\\ {\omega }_{\text{C}Z⟨k⟩}& 0& -{\omega }_{\text{C}X⟨k⟩}\\ -{\omega }_{\text{C}Y⟨k⟩}& {\omega }_{\text{C}X⟨k⟩}& 0\end{array}\right]\\ =& \frac{n}{t_{\text{E}}}\left({\mathit{R}}_{\text{C}⟨k⟩}-{\mathit{R}}_{\text{C}⟨k-1⟩}\right){\left({\mathit{R}}_{\text{C}⟨k-1⟩}\right)}^{\text{T}},\end{array}$

and the velocity can be calculated as

(27)${\mathit{v}}_{\text{C}⟨k⟩}=\frac{n}{t_{\text{E}}}({\mathit{o}}_{\text{C}⟨k⟩}-{\mathit{o}}_{\text{C}⟨k-1⟩}).$

Thus, set ${\mathbb{V}}_{\text{C}}$ can be defined to represent the velocity characteristics of the motion trajectory as follows:

(28)${\mathbb{V}}_{\text{C}}=\left\{\left[\begin{array}{c}{\text{ω}}_{\text{C}⟨0⟩}\\ {\mathit{v}}_{\text{C}⟨0⟩}\end{array}\right],\left[\begin{array}{c}{\text{ω}}_{\text{C}⟨1⟩}\\ {\mathit{v}}_{\text{C}⟨1⟩}\end{array}\right],\left[\begin{array}{c}{\text{ω}}_{\text{C}⟨2⟩}\\ {\mathit{v}}_{\text{C}⟨2⟩}\end{array}\right],\dots ,\left[\begin{array}{c}{\text{ω}}_{\text{C}⟨n⟩}\\ {\mathit{v}}_{\text{C}⟨n⟩}\end{array}\right]\right\}.$

Similarly, we can approximate the instantaneous acceleration with the average acceleration in the kth step as follows:

(29)$\left[\begin{array}{c}{\text{α}}_{\text{C}⟨k⟩}\\ {\mathit{a}}_{\text{C}⟨k⟩}\end{array}\right]=\frac{n}{t_{\text{E}}}\left[\begin{array}{c}{\text{ω}}_{\text{C}⟨k⟩}-{\text{ω}}_{\text{C}⟨k-1⟩}\\ {\mathit{v}}_{\text{C}⟨k⟩}-{\mathit{v}}_{\text{C}⟨k-1⟩}\end{array}\right].$

Thus, set ${\mathbb{A}}_{\text{C}}$ can be defined to represent the acceleration characteristics of the motion trajectory as follows:

(30)${\mathbb{A}}_{\text{C}}=\left\{\left[\begin{array}{c}{\text{α}}_{\text{C}⟨0⟩}\\ {\mathit{a}}_{\text{C}⟨0⟩}\end{array}\right],\left[\begin{array}{c}{\text{α}}_{\text{C}⟨1⟩}\\ {\mathit{a}}_{\text{C}⟨1⟩}\end{array}\right],\left[\begin{array}{c}{\text{α}}_{\text{C}⟨2⟩}\\ {\mathit{a}}_{\text{C}⟨2⟩}\end{array}\right],\dots ,\left[\begin{array}{c}{\text{α}}_{\text{C}⟨n⟩}\\ {\mathit{a}}_{\text{C}⟨n⟩}\end{array}\right]\right\}.$

We predefine the initial kinematic parameters of this trajectory, including velocity and acceleration. To further demonstrate the effectiveness of the FSI method in motion tracking, we continue to base our analysis on the numerical example of the 6-UPS mechanism and design the following set of driving functions:

(31)$\mathbb{F}=\left\{\begin{array}{c}{f}_1\left(t\right)\\ f_2\left(t\right)\\ f_3\left(t\right)\\ f_4\left(t\right)\\ f_5\left(t\right)\\ f_6\left(t\right)\end{array}\right\}=\left\{\begin{array}{c}40\sin (\text{π}t)+c_{1⟨0⟩}\\ -40\sin (\text{π}t)+c_{2⟨0⟩}\\ 20\sin (\text{π}t)+c_{3⟨0⟩}\\ -20\sin (\text{π}t)+c_{4⟨0⟩}\\ 40\sin (\text{π}t)+c_{5⟨0⟩}\\ -40\sin (\text{π}t)+c_{6⟨0⟩}\end{array}\right\}.$

The initial kinematic parameters of all actuators can be obtained through derivative operations as follows:

(32)$\begin{array}{rl}& \left[\begin{array}{cccc}{\dot{f}}_1(0)& {\dot{f}}_2(0)& \cdots & {\dot{f}}_6(0)\\ {\ddot{f}}_1(0)& {\ddot{f}}_2(0)& \cdots & {\ddot{f}}_6(0)\end{array}\right]\\ =& \text{π}\left[\begin{array}{cccccc}40& -40& 20& -20& 80& -80\\ 0& 0& 0& 0& 0& 0\end{array}\right].\end{array}$

From the inverse kinematic equations, the velocity and acceleration mappings for the 6-UPS PM are as follows:

(33)$\{\begin{array}{rl}{\dot{f}}_i=& {\left({\text{ω}}_{\text{D}}\times {\mathit{s}}_i\right)}^{\text{T}}\cdot {\mathit{p}}_i+{\mathit{v}}_{\text{D}}^{\text{T}}\cdot {\mathit{p}}_i,\\ {\ddot{f}}_i=& {\left({\text{α}}_{\text{D}}\times {\mathit{s}}_i\right)}^{\text{T}}\cdot {\mathit{p}}_i+{\left({\text{ω}}_{\text{D}}\times {\dot{\mathit{s}}}_i\right)}^{\text{T}}\cdot {\mathit{p}}_i+{\left({\text{ω}}_{\text{D}}\times {\mathit{s}}_i\right)}^{\text{T}}\cdot {\dot{\mathit{p}}}_1\\ & +{\mathit{a}}_{\text{D}}^{\text{T}}\cdot {\mathit{p}}_i+{\mathit{v}}_{\text{D}}^{\text{T}}\cdot {\dot{\mathit{p}}}_i,\end{array}$

where ${\text{ω}}_{\text{D}}$ and ${\mathit{v}}_{\text{D}}$ are the angular velocity and linear velocity of ${\mathit{P}}_{\text{D}}$, respectively.

We establish a kinematic simulation model in Hexagon Adams and use it to output the data for the control group. After computation, a series of motion tracking characteristic curves is obtained for the 6-UPS PM. All motion tracking characteristic curves and their errors are shown in Fig.16, where the control curves are sourced from the Adams simulation data. The results indicate that the FSI method also exhibits high applicability in the motion tracking of the PM. Scholars often use kinematic simulation software such as Hexagon Adams to obtain FK tracking curves. The FSI method can replace simulation software to a certain accuracy level, presenting promising application prospects.

Fig.16 Motion tracking characteristic curves: (a) posture, (b) angular velocity, (c) angular acceleration, (d) posture error, (e) angular velocity error, (f) angular acceleration error, (g) position, (h) velocity, (i) acceleration, (j) position error, (k) velocity error, and (l) acceleration error.

Full size|PPT slide

8 8 Conclusions and prospects

This work proposes a systematic method for FK analysis, termed FSI, that is applicable to PMs and based on the principles of calculus and mathematical induction. The FSI method transforms static spatial geometric solving problems into linear superimposed models of quasi-transient vectors, thereby converting a non-linear problem into an iterative and cumulative multi-step linear problem. This method represents a novel concept that challenges traditional views and holds innovative and practical value. Basing on the comprehensive content of this paper, we believe that the FSI method possesses the following characteristics:

● The FSI method employs a fully linearized algebraic system, where the solving process involves repeatedly solving full-rank square matrices. This process imposes minimal computational burden on the hardware and ensures the existence of solutions.

● This systematic and modularized efficient modeling approach enables scholars to quickly establish FK equations for any PM configuration and convert them into program code.

● The FSI method has friendly mathematical foundation requirements, utilizing relatively simple mathematical tools to solve complex problems and facilitate scholars’ understanding and application.

● Its unique algebraic structure allows scholars to perform kinematic analysis on the entire motion trajectory of PMs and identify singularities during the motion process using simple criteria.

● Its targeted objective solving avoids issues related to multiple solutions and filtering. The FSI method only tracks from the initial configuration to target points within the reachable workspace and cannot cross singular boundaries.

We also need to address the current limitations of the FSI method, such as its inability to achieve ultra-high precision calculations and the need for further improvement in calculation speed. In future research, we plan to investigate the capabilities of the FSI method in computing reachable workspace and singular boundaries and to integrate it with dynamic analysis. We also intend to apply the theoretical principles of the FSI method in analyzing the dynamic problems of PMs. In conclusion, the FSI method is an innovative theoretical approach that challenges traditional concepts in PMs’ FK analysis and has achieved promising results. It is expected to effectively solve all non-linear challenges related to PM kinematics, demonstrating its significant potential for widespread application.

References

[1] Pilliar R M. Porous-surfaced metallic implants for orthopedic applications. Journal of Biomedical Materials Research, 1987, 21(A1 Suppl): 1-33

[2] Kienapfel H, Sprey C, Wilke A, Griss P. Implant fixation by bone ingrowth. The Journal of Arthroplasty, 1999, 14(3): 355-368

[3] Wen C E, Mabuchi M, Yamada Y, Shimojima K, Chino Y, Asahina T. In: Processing of biocompatible porous Ti and Mg. Scripta Materialia, 2001, 45(10): 1147-1153

[4] Karageorgiou V, Kaplan D. Porosity of 3D biomaterial scaffolds and osteogenesis. Biomaterials, 2005, 26(27): 5474-5491

[5] Hollister S J. Porous scaffold design for tissue engineering. Nature Materials, 2005, 4(7): 518-524

[6] Otsuki B, Takemoto M, Fujibayashi S, Neo M, Kokubo T, Nakamura T. Pore throat size and connectivity determine bone and tissue ingrowth into porous implants: three-dimensional micro-CT based structural analyses of porous bioactive titanium implants. Biomaterials, 2006, 27(35): 5892-5900

[7] Jones A C, Arns C H, Hutmacher D W, Milthorpe B K, Sheppard A P, Knackstedt M A. The correlation of pore morphology, interconnectivity and physical properties of 3D ceramic scaffolds with bone ingrowth. Biomaterials, 2009, 30(7): 1440-1451

[8] Hutmacher D W, Sittinger M, Risbud M V. Scaffold-based tissue engineering: rationale for computer-aided design and solid free-form fabrication systems. Trends in Biotechnology, 2004, 22(7): 354-362

[9] Williams J M, Adewunmi A, Schek R M, Flanagan C L, Krebsbach P H, Feinberg S E, Hollister S J, Das S. Bone tissue engineering using polycaprolactone scaffolds fabricated via selective laser sintering. Biomaterials, 2005, 26(23): 4817-4827

[10] Seitz H, Rieder W, Irsen S, Leukers B, Tille C. Three-dimensional printing of porous ceramic scaffolds for bone tissue engineering. Journal of Biomedical Materials Research, Part B, Applied Biomaterials, 2005, 74(2): 782-788

[11] Ryan G E, Pandit A S, Apatsidis D P. Porous titanium scaffolds fabricated using a rapid prototyping and powder metallurgy technique. Biomaterials, 2008, 29(27): 3625-3635

[12] Li J P, de Wijn J R, van Blitterswijk C A, de Groot K. Porous Ti₆Al₄V scaffold directly fabricating by rapid prototyping: preparation and in vitro experiment. Biomaterials, 2006, 27(8): 1223-1235

[13] Krishna B V, Bose S, Bandyopadhyay A. Low stiffness porous Ti structures for load-bearing implants. Acta Biomaterialia, 2007, 3(6): 997-1006

[14] Hollander D A, von Walter M, Wirtz T, Sellei R, Schmidt-Rohlfing B, Paar O, Erli H J. Structural, mechanical and in vitro characterization of individually structured Ti-6Al-4V produced by direct laser forming. Biomaterials, 2006, 27(7): 955-963

[15] Heinl P, Müller L, Körner C, Singer R F, Müller F A. Cellular Ti-6Al-4V structures with interconnected macro porosity for bone implants fabricated by selective electron beam melting. Acta Biomaterialia, 2008, 4(5): 1536-1544

[16] Kokubo T, Takadama H. How useful is SBF in predicting in vivo bone bioactivity? Biomaterials, 2006, 27(15): 2907-2915

[17] Kujala S, Ryhänen J, Danilov A, Tuukkanen J. Effect of porosity on the osteointegration and bone ingrowth of a weight-bearing nickel-titanium bone graft substitute. Biomaterials, 2003, 24(25): 4691-4697

[18] Li J P, Habibovic P, van den Doel M, Wilson C E, de Wijn J R, van Blitterswijk C A, de Groot K. Bone ingrowth in porous titanium implants produced by 3D fiber deposition. Biomaterials, 2007, 28(18): 2810-2820

RIGHTS & PERMISSIONS

2014 Higher Education Press and Springer-Verlag Berlin Heidelberg