The kinematics analysis and improvement of the 3T1R 4-PUU parallel mechanism using the whale optimization algorithm

This paper focuses on the 4-PUU parallel mechanism with Schönflies motion, investigating its kinematic analysis and utilizing intelligent algorithms to solve its forward kinematics equations. Using screw theory as the mathematical framework, the number and nature of the degrees of freedom of the 4-PUU parallel mechanism are analyzed, proving that the mechanism can achieve 3T1R motion. The forward and inverse position kinematics equations are established with link lengths as constraints, and the velocity and acceleration equations of the mechanism are derived using the vector method. The forward kinematics equations of the mechanism are transformed into an unconstrained optimization problem, which is solved using an improved beluga whale optimization (BWO). To enhance the uniformity of the initial population distribution, a chaotic-opposition-based learning initialization strategy is introduced. Additionally, an elite strategy and a golden-sine-based position update mechanism are incorporated to improve the optimization capability of BWO. The integration of these strategies results in an enhanced optimization algorithm with superior global search ability. Given the mechanism's dimensional parameters, theoretical and motion simulations are conducted using MATLAB and Adams. The results indicate that the motion simulation curves from Adams closely match the numerical simulation curves from MATLAB, validating the theoretical derivation. Building on this, the improved beluga whale optimization (IBWO) algorithm, along with comparison algorithms (DEb1, ABC, PSO), is applied to solve the forward kinematics equations for interpolation points along the simulated trajectory. Numerical experiments demonstrate that IBWO outperforms the other comparison algorithms in solving the problem efficiently and accurately.

1 Introduction

Compared to 6-degree-of-freedom (DOF) parallel mechanisms (PMs), those with fewer DOFs offer advantages such as simpler structure, fewer drives, and easier control and have gradually gained widespread attention from both academia and industry. Among them, the 3T1R (T representing translation and R representing revolute) PM is widely used in practical engineering applications such as packaging and sorting (Zhu et al., 2022), welding (Wen et al., 2024), assembly (Li et al., 2021), and rehabilitation (Shi et al., 2024).

The original SCARA – with 3T1R DOF (where T and R represent translational and rotational degrees of freedom, respectively) – is an extension of the 3-DOF Delta robot. An additional UPU chain was incorporated between the moving and fixed platforms to provide rotation about the z axis. The kinematic analysis, performance analysis, error compensation, and trajectory planning of the PM are all based on the integration of its configuration. Therefore, it is evident that the configuration synthesis of the SCARA of PM forms the foundation for subsequent theoretical analysis and structural design.

The configuration synthesis of SCARA PMs has shown significant research achievements. Company et al. (2003) designed a structure with both symmetric and asymmetric H4 PM based on the motion used to generate Schönflies displacement subgroups. Yang et al. (2017) proposed a configuration synthesis method for a 3T1R PM with a variable rotation axis based on screw theory. Salgado et al. (2008) used the theory of groups of displacements, based on the H4, I4 (Krut et al., 2003), and Par4 (Nabat et al., 2005) mechanisms. Li et al. (2014) replaced the original parallelogram chains with R(SS)2 chains, introducing the Cross IV series PMs, which were later industrialized. Staicu et al. (2018) proposed a fully symmetric X4 PM where each of the four chains is composed of R(SS)R, resulting in improved dynamic characteristics and simplified manufacturing. Wen et al. (2024) designed a 2PUU-2PSS PM with an axially symmetric structure based on screw theory. Xie et al. (2013), using Grassmann line geometry and the Atlas method, designed a class of PMs with 3T1R (1T3R, 2T2R, and 3T1R) motion.

Kinematic analysis of PMs forms the foundation for error compensation, trajectory planning, and workspace analysis. It primarily consists of forward kinematics analysis and inverse kinematics analysis. For coupled PMs, the inverse kinematics analysis typically involves decoupling the equations, while the forward kinematics usually leads to a coupled nonlinear system of equations. The solution to the forward kinematics equations of coupled PMs is mainly divided into analytical method and numerical method. Analytical methods involve eliminating unknowns and reducing higher-order equations to first-order equations. By solving the first-order equations and performing inverse substitution of the results, a high-precision solution to the forward kinematics can be obtained. However, this method is specific to certain configurations and lacks generality. Numerical methods, on the other hand, include numerical iteration, machine learning, and intelligent algorithms, which can be applied to a wider range of mechanisms and offer more flexibility in handling complex kinematic problems.

Wen et al. (2024) proposed a mean-based dual-strategy differential evolution algorithm for solving the forward kinematics equations of the 2PUU-2PSS PM. It was demonstrated that the performance of this algorithm surpasses that of comparison algorithms such as IABC and DEr1 through numerical examples. Liu et al. (2019) combined genetic algorithms and backpropagation (BP) neural networks to solve the forward kinematics equations of the Delta robot, with results showing better performance than traditional BP neural networks. Wu et al. (2024) introduced an improved backpropagation neural network for solving the forward kinematics equations of redundant PMs, and experimental results indicated a significant improvement in solution accuracy. Wen et al. (2023) proposed a time-varying differential evolution algorithm for solving the multiple solution problem of forward kinematics equations in PMs. Using PMs with 3 to 6 DOF as examples, the feasibility and versatility of the algorithm were verified.

This paper is based on the configuration synthesis of PMs and utilizes screw theory to design a 4-PUU PM with a fully symmetric structure, which the workspace is unaffected by the vertical direction. The structure and kinematic analysis of the 4-PUU PM are presented in the following chapters: Sect. 2 provides a structural description of the 4-PUU PM and analyzes its degrees of freedom using screw theory. Section 3 derives the kinematic equations of the mechanism using vector and differentiation methods. Section 4 introduces an improved whale optimization algorithm (IBWO) to solve the forward kinematics equations of the 4-PUU PM. Section 5 verifies the correctness of the kinematic formula derivations through motion simulations using Adams software and validates the performance of the IBWO algorithm through numerical examples. Section 6 concludes the paper and outlines future research directions.

2 4-PUU structural description and DOF analysis of PMs

2.1 Structural description

The 4-PUU PM is shown in Fig. 1. This mechanism mainly consists of a moving platform, a fixed platform, and four identical structural branches. Each branch is composed of a prismatic pair (P) and universal pair (U), which are rotational pairs with axes perpendicular to each other. The topology of the branch can be represented as SCO{-P₁∥R₂⊥R₃∥R₄⊥R₅}.That is, starting from the fixed platform, the axis of the moving pair 1 is parallel to the axis of the rotating pair 2, the axis of the rotating pair 2 is perpendicular to the axis of the rotating pair 3, the axis of the rotating pair 3 is parallel to the axis of the rotating pair 4, the axis of the rotating pair 4 is perpendicular to the axis of the rotating pair 5, and the axis of the rotating pair 5 is perpendicular to the moving platform. In the PUU chain, the intersection of the axis of the pair (P) with the fixed platform is denoted as $A_i(i=\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{4})$. The center of the first universal joint intersects the axis of the pair at point B_i, and the center of the second universal joint intersects the moving platform at point C_i. To avoid the initial singularity of the mechanism, the fixed platform is set as a square with the points A_i lying on a circle with radius $\sqrt{\mathrm{2}}R$. The moving platform is set as a rectangle with side lengths of 2a and 2b.

Figure 1 Three-dimensional diagram of the 4-PUU PM.

• A fixed coordinate system {f}: OXYZ is established on the fixed platform, where O is the geometric center of the fixed platform. The OX axis is perpendicular to A₁A₄ and passes through its midpoint, the OY axis is perpendicular to A₁A₂ and passes through its midpoint, and the OZ axis is perpendicular to the fixed platform and points vertically upwards.

• A moving coordinate system {m}: Pxyz is established on the moving platform, where P is the geometric center of the moving platform. The Px axis is perpendicular to C₁C₄ and passes through its midpoint, the Py axis is perpendicular to C₁C₂ and passes through its midpoint, and the Pz axis is perpendicular to the moving platform and points vertically upwards.

Figure 2 Schematic diagram of PUU branch chain structure.

2.2 Analysis of DOF

The branch coordinate system is established as shown in Fig. 2: A_i–x_iy_iz_i, where the z_i axis coincides with A_iB_i and points toward B_i, the y_i axis coincides with A_iO and points toward O, and the x_i axis is determined by the right-hand rule. The angle between the revolute pair R3 and the x_i axis is θ_I, and the angle between the link B_iC_i and the axis of the prismatic joint P is ζ_i. Thus, in the branch coordinate system, the motion screw of the PUU branch can be expressed as

$$\begin{array}{}\text{(1)}& \left\{\begin{array}{l}{\mathit{\$}}_{i\mathrm{1}}=(\mathrm{0},\mathrm{0},\mathrm{0};\mathrm{0},\mathrm{0},\mathrm{1})\\ {\mathit{\$}}_{i\mathrm{2}}=(\mathrm{0},\mathrm{0},\mathrm{1};\mathrm{0},\mathrm{0},\mathrm{0})\\ {\mathit{\$}}_{i\mathrm{3}}=(c{\mathit{\theta }}_i,s{\mathit{\theta }}_i,\mathrm{0};h_{i}s{\mathit{\theta }}_i,-h_{i}c{\mathit{\theta }}_i,\mathrm{0})\\ \begin{array}{rl}{\mathit{\$}}_{i\mathrm{4}}=& (c{\mathit{\theta }}_i,s{\mathit{\theta }}_i,\mathrm{0};(ls{\mathit{\zeta }}_i-h_i)s{\mathit{\theta }}_i,\\ & (h_i-ls{\mathit{\zeta }}_i)c{\mathit{\theta }}_i,-lc{\mathit{\zeta }}_i)\end{array}\\ {\mathit{\$}}_{i\mathrm{5}}=(\mathrm{0},\mathrm{0},\mathrm{1};-lc{\mathit{\delta }}_{i}c{\mathit{\theta }}_i,-lc{\mathit{\delta }}_{i}s{\mathit{\theta }}_i,\mathrm{0})\end{array}\right.,\end{array}$$

where cosθ_i=cθ_i, sinθ_i=sθ_i, cosξ_i=cξ_i, and sinξ_i=sξ_i.

Therefore, the constraint screw of the PUU chain can be expressed as

$$\begin{array}{}\text{(2)}& {}^r{\mathit{\$}}_i=(\mathrm{0},\mathrm{0},\mathrm{0};s{\mathit{\theta }}_i,-c{\mathit{\theta }}_i,\mathrm{0}).\end{array}$$

Equation (2) represents a reciprocal screw pair, which indicates that the branch constrains the moving platform's rotation around the direction parallel to (sinθ_i, −cosθ_i, 0). The constraints imposed by the four symmetrically distributed PUU branches on the moving platform can be expressed using the two reciprocal screws ${\mathit{\$}}_{\mathrm{11}}^r$ and ${\mathit{\$}}_{\mathrm{13}}^r$. Clearly, these two reciprocal screws are linearly dependent. Therefore, a linear combination of the following two reciprocal screws can be used to obtain

$$\begin{array}{}\text{(3)}& \left\{\begin{array}{l}{}^r{\mathit{\$}}_{\mathrm{1}}=(\mathrm{0},\mathrm{0},\mathrm{0};\mathrm{1},\mathrm{0},\mathrm{0})\\ {}^r{\mathit{\$}}_{\mathrm{3}}=(\mathrm{0},\mathrm{0},\mathrm{0};\mathrm{0},\mathrm{1},\mathrm{0})\end{array}\right..\end{array}$$

According to Eq. (3), the constraint screws applied to the moving platform by the four branches are

$$\begin{array}{}\text{(4)}& \left\{\begin{array}{l}{\mathit{\$}}_{\mathrm{CWS},\mathrm{1}}=(\mathrm{0},\mathrm{0},\mathrm{0};\mathrm{1},\mathrm{0},\mathrm{0})\\ {\mathit{\$}}_{\mathrm{CWS},\mathrm{2}}=(\mathrm{0},\mathrm{0},\mathrm{0};\mathrm{0},\mathrm{1},\mathrm{0})\end{array}\right..\end{array}$$

$_CWS,1 and $_CWS,2 represent the constraint wrenches around the X axis and Y axis, respectively, applying rotational constraints around the X and Y axes to the moving platform, thus causing the platform to lose its rotational freedom around the X and Y axes. Therefore, the 4-PUU PM proposed in this paper has 4 DOF(3TIR_Z). In general configurations (except for singular poses), the motion pairs of each branch maintain their original geometric relationships due to structural constraints. The coordinate relations shown in Fig. 2 can still be established, and the resulting motion screws and the computed reciprocal screws remain unchanged. The constraint wrenches of the four reciprocal screw systems remain in the same plane, and their linear independence and rank remain constant. That is, the virtual constraints of the mechanism remain unchanged. Thus, the DOFs correspond to full rotational freedom.

3 4-PUU PM kinematic analysis

3.1 Forward and inverse kinematics analysis

3.1.1 Forward kinematics

The forward kinematics solution of the mechanism is to calculate the output position P= (X_P, Y_P, Z_P) and the attitude angle γ of the moving platform's reference point P, given the scale parameters of the mechanism and the input link lengths h_i (for i=1, …, 4) of the prismatic joints. The structural parameters of the platform are set as ${\mathit{\sigma }}_i=(\mathrm{2}i-\mathrm{1})\mathit{\pi }/\mathrm{4}$, and the coordinates of point A_i in the frame {f} are defined as

$$\begin{array}{}\text{(5)}& {}^{\mathit{\left\{}f\mathit{\right\}}}{\mathbf{A}}_i=R(c{\mathit{\sigma }}_i,s{\mathit{\sigma }}_i,\mathrm{0}{)}^T.\end{array}$$

The coordinate vector of point C_i in the frame {m} is

$$\begin{array}{}\text{(6)}& {}^{\mathit{\left\{}m\mathit{\right\}}}{\mathbf{C}}_i=\left(\mathrm{sign}\right(c{\mathit{\sigma }}_i)a,\mathrm{sign}(s{\mathit{\sigma }}_i)b,\mathrm{0}{)}^T.\end{array}$$

The unit vector of the axis of the driving arm A_iB_i in the frame {f} is

$$\begin{array}{}\text{(7)}& {}^{\mathit{\left\{}\mathrm{f}\mathit{\right\}}}{\mathbf{u}}_i=(\mathrm{0},\mathrm{0},\mathrm{1}{)}^T.\end{array}$$

Assume the coordinate vector of the reference point P on the moving platform in the fixed frame {f} is

$$\begin{array}{}\text{(8)}& {}^{\mathit{\left\{}f\mathit{\right\}}}\mathbf{P}=(X_P,Y_P,Z_P{)}^T.\end{array}$$

The attitude matrix of the moving platform is

$$\mathbf{M}=\left(\begin{array}{ccc}c\mathit{\gamma }& -s\mathit{\gamma }& \mathrm{0}\\ s\mathit{\gamma }& c\mathit{\gamma }& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right).$$

Then, the coordinate vector of point C_i(i=1, …, 4) in the fixed coordinate system {f} is

$$\begin{array}{}\text{(9)}& {}^{\mathit{\left\{}f\mathit{\right\}}}{\mathbf{C}}_i{=}^{\mathit{\left\{}f\mathit{\right\}}}\mathbf{P}+{\mathbf{M}}^{\mathit{\left\{}m\mathit{\right\}}}{\mathbf{C}}_i{=}^{\mathit{\left\{}f\mathit{\right\}}}{\mathbf{A}}_i+h_i^{\mathit{\left\{}f\mathit{\right\}}}{\mathbf{u}}_i+l^{\mathit{\left\{}f\mathit{\right\}}}{\mathit{\omega }}_i.\end{array}$$

The unit vector of the axis of the moving arm B_iC_i in the fixed frame {f} is

$$\begin{array}{}\text{(10)}& {\mathit{\omega }}_i={\displaystyle \frac{{}^{\mathit{\left\{}f\mathit{\right\}}}\mathbf{P}+\mathbf{M}{}^{\mathit{\left\{}m\mathit{\right\}}}{\mathbf{C}}_i{-}^{\mathit{\left\{}f\mathit{\right\}}}{\mathbf{A}}_i-h_i^{\mathit{\left\{}f\mathit{\right\}}}{\mathbf{u}}_i}{l}}.\end{array}$$

The constraint equation of the link B_iC_i is

$$\begin{array}{}\text{(11)}& \begin{array}{rl}{(}^{\mathit{\left\{}f\mathit{\right\}}}\mathbf{P}& +{\mathbf{M}}^{\mathit{\left\{}m\mathit{\right\}}}{\mathbf{C}}_i{-}^{\mathit{\left\{}f\mathit{\right\}}}{\mathbf{A}}_i-h_i^{\mathit{\left\{}f\mathit{\right\}}}{\mathbf{u}}_i{)}^T{(}^{\mathit{\left\{}f\mathit{\right\}}}\mathbf{P}+{\mathbf{M}}^{\mathit{\left\{}\mathrm{m}\mathit{\right\}}}{\mathbf{C}}_i\\ & {-}^{\mathit{\left\{}f\mathit{\right\}}}{\mathbf{A}}_i-h_i^{\mathit{\left\{}f\mathit{\right\}}}{\mathbf{u}}_i)=l^{\mathrm{2}}.\end{array}\end{array}$$

According to Eq. (11), we can obtain

$$\begin{array}{}\text{(12)}& \left\{\begin{array}{l}(R-ac\mathit{\gamma }-bs\mathit{\gamma }-X_P{)}^{\mathrm{2}}+(R+ac\mathit{\gamma }-bs\mathit{\gamma }-Y_P{)}^{\mathrm{2}}+(h_{\mathrm{1}}-Z_P{)}^{\mathrm{2}}=l^{\mathrm{2}}\\ (-R+ac\mathit{\gamma }-bs\mathit{\gamma }-X_P{)}^{\mathrm{2}}+(R-ac\mathit{\gamma }-bs\mathit{\gamma }-Y_P{)}^{\mathrm{2}}+(h_{\mathrm{2}}-Z_P{)}^{\mathrm{2}}=l^{\mathrm{2}}\\ (-R+ac\mathit{\gamma }+bs\mathit{\gamma }-X_P{)}^{\mathrm{2}}+(-R-ac\mathit{\gamma }+bs\mathit{\gamma }-Y_P{)}^{\mathrm{2}}+(h_{\mathrm{3}}-Z_P{)}^{\mathrm{2}}=l^{\mathrm{2}}\\ (R-ac\mathit{\gamma }+bs\mathit{\gamma }-X_P{)}^{\mathrm{2}}+(-R+ac\mathit{\gamma }+bs\mathit{\gamma }-Y_P{)}^{\mathrm{2}}+(h_{\mathrm{4}}-Z_P{)}^{\mathrm{2}}=l^{\mathrm{2}}\end{array}\right..\end{array}$$

Equation (12) is a system of coupled nonlinear equations, which can be solved using both numerical methods and analytical methods. The analytical method involves substitution and elimination to obtain a single high-order equation. By solving this high-order equation and performing backward substitution, a high-precision forward kinematics solution for the mechanism can be obtained. However, this method has limited applicability and is only suitable for special configurations of PMs. The numerical method, on the other hand, solves the equations using intelligent algorithms or numerical iteration, offering advantages such as simplicity and efficiency. Therefore, in this paper, a numerical method is employed to solve the forward kinematics equations, as detailed in Sect. 4.

3.1.2 Inverse kinematics

The inverse kinematics analysis of the mechanism involves the process of solving for the input h_i given the output position and orientation of the mechanism (X_P, Y_P, Z_P, γ). Thus, it can be deduced that, according to Eq. (12), we obtain

$$\begin{array}{}\text{(13)}& h_i=Z_P\pm \sqrt{l^{\mathrm{2}}-{\mathit{\lambda }}_{i\mathrm{1}}^{\mathrm{2}}-{\mathit{\lambda }}_{i\mathrm{2}}^{\mathrm{2}}}.\end{array}$$

In the equation, λ_ij is a coefficient ($i=\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{4}$; $j=\mathrm{1},\mathrm{2}$), and its expression is

$$\begin{array}{}\text{(14)}& \left\{\begin{array}{l}{\mathit{\lambda }}_{i\mathrm{1}}=X_P+\mathrm{sign}\left(\cos {\mathit{\zeta }}_i\right)a\cos \mathit{\gamma }-\mathrm{sign}\left(\sin {\mathit{\zeta }}_i\right)b\sin \mathit{\gamma }-R\cos {\mathit{\zeta }}_i\\ {\mathit{\lambda }}_{i\mathrm{2}}=Y_P+\mathrm{sign}\left(\cos {\mathit{\zeta }}_i\right)a\sin \mathit{\gamma }-\mathrm{sign}\left(\sin {\mathit{\zeta }}_i\right)b\cos \mathit{\gamma }-R\sin {\mathit{\zeta }}_i\\ l^{\mathrm{2}}-{\mathit{\lambda }}_{i\mathrm{1}}^{\mathrm{2}}-{\mathit{\lambda }}_{i\mathrm{2}}^{\mathrm{2}}=\mathrm{0}\end{array}\right..\end{array}$$

3.2 Kinematic velocity analysis

Let the inputs of the four driving joints of the mechanism be h=(h₁, h₂,h₃,h₄)^T, the output quantities of the moving platform be X=(X_P, Y_P,Z_P, γ)^T, the input velocities of the four driving joints be $\mathit{v}=({\dot{h}}_{\mathrm{1}}$, ${\dot{h}}_{\mathrm{2}}$, ${\dot{h}}_{\mathrm{3}}$, ${\dot{h}}_{\mathrm{4}}{)}^T$, the output velocities of the moving platform be $\dot{\mathbf{X}}=({\dot{X}}_P$, ${\dot{X}}_P$, ${\dot{Z}}_P{)}^T$, and the angular velocity output vector of the moving platform be $\dot{\mathit{\gamma }}=(\mathrm{0},\mathrm{0},\dot{\mathit{\gamma }}{)}^T$.

Taking the derivative of Eq. (11) with respect to time t, we get

$$\begin{array}{}\text{(15)}& [\dot{\mathbf{P}}+\dot{\mathit{\gamma }}\times ({\mathbf{Mc}}_i)-{\dot{h}}_i\mathbf{u}{]}^T{\mathit{\omega }}_i=\mathrm{0}i=\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{4}.\end{array}$$

By simplifying this, we obtain

$$\begin{array}{}\text{(16)}& {\mathit{\omega }}_i^T\dot{P}+\left[\right(M{c}_i)\times {\mathit{\omega }}_i{]}^T\dot{\mathit{\gamma }}={\dot{h}}_i{\mathit{\omega }}_i^{T}u.\end{array}$$

Let $k_i={\mathit{\omega }}_i^{T}u$ and take ${\mathbf{J}}_{\mathrm{in}}=\mathrm{diag}(k_{\mathrm{1}},k_{\mathrm{2}},k_{\mathrm{3}},k_{\mathrm{4}})$:

$${\mathbf{J}}_{\mathrm{out}}=\left(\begin{array}{cc}{\mathit{\omega }}_{\mathrm{1}}^T& \left[\right({\mathbf{Mc}}_{\mathrm{1}})\times {\mathit{\omega }}_{\mathrm{1}}{]}^T\mathbf{u}\\ {\mathit{\omega }}_{\mathrm{2}}^T& \left[\right({\mathbf{Mc}}_{\mathrm{2}})\times {\mathit{\omega }}_{\mathrm{2}}{]}^T\mathbf{u}\\ {\mathit{\omega }}_{\mathrm{3}}^T& \left[\right({\mathbf{Mc}}_{\mathrm{3}})\times {\mathit{\omega }}_{\mathrm{3}}{]}^T\mathbf{u}\\ {\mathit{\omega }}_{\mathrm{4}}^T& \left[\right({\mathbf{Mc}}_{\mathrm{4}})\times {\mathit{\omega }}_{\mathrm{4}}{]}^T\mathbf{u}\end{array}\right).$$

It can be seen that Eq. (16) can be rewritten as

$$\begin{array}{}\text{(17)}& {\mathbf{J}}_{\mathrm{in}}\dot{\mathbf{h}}={\mathbf{J}}_{\mathrm{out}}\dot{\mathbf{X}}.\end{array}$$

If J_out is non-singular, then

$$\begin{array}{}\text{(18)}& \dot{\mathbf{X}}={\mathbf{J}}_{\mathrm{out}}^{-\mathrm{1}}{\mathbf{J}}_{\mathrm{in}}\dot{\mathbf{h}}.\end{array}$$

If J_in is non-singular, then

$$\begin{array}{}\text{(19)}& \dot{\mathbf{h}}={\mathbf{J}}_{\mathrm{in}}^{-\mathrm{1}}{\mathbf{J}}_{\mathrm{out}}\dot{\mathbf{X}}.\end{array}$$

Equations (22) and (23) are the forward and inverse velocity solutions of the 4-PUU PM, respectively.

3.3 Acceleration analysis

Let the input accelerations of the mechanism's driving joints be $\ddot{\mathit{h}}=({\dot{h}}_{\mathrm{1}}$, ${\dot{h}}_{\mathrm{2}}$, ${\dot{h}}_{\mathrm{3}}$, ${\dot{h}}_{\mathrm{4}}{)}^T$ and the output acceleration of the moving platform be $\ddot{\mathbf{X}}=({\ddot{X}}_P,{\ddot{Y}}_P,{\ddot{Z}}_P{)}^T$, with the output angular velocity being $\ddot{\mathit{\gamma }}=(\mathrm{0},\mathrm{0},\ddot{\mathit{\gamma }}{)}^T$. Taking the derivative of equation with respect to time t, we get

$$\begin{array}{}\text{(20)}& \begin{array}{rl}\mathit{\{}\ddot{\mathbf{P}}+& \ddot{\mathit{\gamma }}\times \left(\mathrm{M}{\mathbf{c}}_i\right)+\dot{\mathit{\gamma }}\times [\dot{\mathit{\gamma }}\times (\mathrm{M}{\mathbf{c}}_i)-{\dot{h}}_i\mathbf{u}]{\mathit{\}}}^T{\mathit{\omega }}_i\\ & +{\displaystyle \frac{\mathrm{1}}{l}}[\dot{\mathbf{P}}+\dot{\mathit{\gamma }}\times (\mathrm{M}{\mathbf{c}}_i)-{\dot{h}}_i\mathbf{u}{]}^T[\dot{\mathbf{P}}+\dot{\mathit{\gamma }}\times \left(\mathrm{M}{\mathbf{c}}_i\right)\\ & -{\dot{h}}_i\mathbf{u}{]}^T=\mathrm{0}.\end{array}\end{array}$$

If we denote ${\mathit{\kappa }}_i=({\mathit{\omega }}_i\times \dot{\mathit{\gamma }}{)}^T[\dot{\mathit{\gamma }}\times \left(\mathrm{M}{c}_i\right)]+\frac{\mathrm{1}}{l}∥\dot{P}+\dot{\mathit{\gamma }}\times \left(\mathrm{M}{c}_i\right)-{\dot{\mathbf{h}}}_{i}u∥$, then we have $\mathit{\kappa }=({\mathit{\kappa }}_{\mathrm{1}},{\mathit{\kappa }}_{\mathrm{2}},{\mathit{\kappa }}_{\mathrm{3}},{\mathit{\kappa }}_{\mathrm{4}}{)}^T$.

It can be seen that Eq. (20) can be simplified to

$$\begin{array}{}\text{(21)}& {\mathbf{J}}_{\mathrm{in}}\ddot{\mathbf{h}}={\mathbf{J}}_{\mathrm{out}}\ddot{\mathbf{X}}+\mathit{\kappa }.\end{array}$$

If J_out is non-singular, then

$$\begin{array}{}\text{(22)}& \ddot{\mathbf{X}}={\mathbf{J}}_{\mathrm{out}}^{-\mathrm{1}}{\mathbf{J}}_{\mathrm{in}}\ddot{\mathbf{h}}+{\mathbf{J}}_{\mathrm{out}}^{-\mathrm{1}}\mathit{\kappa }.\end{array}$$

If J_in is non-singular, then

$$\begin{array}{}\text{(23)}& \ddot{\mathbf{h}}={\mathbf{J}}_{\mathrm{in}}^{-\mathrm{1}}{\mathbf{J}}_{\mathrm{out}}\ddot{\mathbf{X}}+{\mathbf{J}}_{\mathrm{in}}^{-\mathrm{1}}\mathit{\kappa }.\end{array}$$

Equations (22) and (23) are the forward and inverse acceleration solutions of the 4-PUU PM, respectively.

4 Beluga whale optimization for forward kinematics analysis

4.1 Mathematical modeling and beluga whale optimization (BWO)

4.1.1 Mathematical modeling

When solving the forward kinematics (Eq. 12) of the PM using a swarm intelligence algorithm, it is necessary to first convert it into an unconstrained optimization problem and then use the swarm intelligence algorithm to iteratively optimize the unconstrained optimization equation. Thus, we define x=(x₁, x₂, x₃, x₄)^T= (X_P, Y_P, Z_P, γ)^T, and Eq. (12) can be transformed into

$$\begin{array}{}\text{(24)}& \left\{\begin{array}{l}\begin{array}{l}{f}_{\mathrm{1}}\left(\mathbf{x}\right)=(R-a\mathrm{c}\mathit{\gamma }-b\mathrm{s}\mathit{\gamma }-X_P{)}^{\mathrm{2}}+(R+a\mathrm{c}\mathit{\gamma }-b\mathrm{s}\mathit{\gamma }-Y_P{)}^{\mathrm{2}}\\ +(h_{\mathrm{1}}-Z_P{)}^{\mathrm{2}}-l^{\mathrm{2}}\end{array}\\ \begin{array}{l}{f}_{\mathrm{2}}\left(\mathbf{x}\right)=(-R+a\mathrm{c}\mathit{\gamma }-b\mathrm{s}\mathit{\gamma }-X_P{)}^{\mathrm{2}}+(R-a\mathrm{c}\mathit{\gamma }-b\mathrm{s}\mathit{\gamma }-Y_P{)}^{\mathrm{2}}\\ +(h_{\mathrm{2}}-Z_P{)}^{\mathrm{2}}-l^{\mathrm{2}}\end{array}\\ \begin{array}{l}{f}_{\mathrm{3}}\left(\mathbf{x}\right)=(-R+a\mathrm{c}\mathit{\gamma }+b\mathrm{s}\mathit{\gamma }-X_P{)}^{\mathrm{2}}+(-R-a\mathrm{c}\mathit{\gamma }+b\mathrm{s}\mathit{\gamma }-Y_P{)}^{\mathrm{2}}\\ +(h_{\mathrm{3}}-Z_P{)}^{\mathrm{2}}-l^{\mathrm{2}}\end{array}\\ \begin{array}{l}{f}_{\mathrm{4}}\left(\mathbf{x}\right)=(R-a\mathrm{c}\mathit{\gamma }+b\mathrm{s}\mathit{\gamma }-X_P{)}^{\mathrm{2}}+(-R+a\mathrm{c}\mathit{\gamma }+b\mathrm{s}\mathit{\gamma }-Y_P{)}^{\mathrm{2}}\\ +(h_{\mathrm{4}}-Z_P{)}^{\mathrm{2}}-l^{\mathrm{2}}\end{array}\end{array}\right..\end{array}$$

The nonlinear system of equations in Eq. (24) can be transformed into an unconstrained optimization problem as follows: $minF\left(\mathbf{x}\right)=\sum _{i=\mathrm{1}}^{\mathrm{4}}{f}_i\left(\mathbf{x}\right)$.

4.1.2 Beluga whale optimization (BWO)

The BWO (Zhong et al., 2022) is a bio-inspired optimization algorithm that simulates the hunting behavior of the beluga whale. The BWO algorithm has characteristics such as a simple structure, easy expandability, few control parameters, and strong robustness. The main process of this algorithm includes the following steps: population initialization, encircling prey, bubble-net attacking, and searching for prey.

1. Population initialization

   The initialization of the population in the BWO is similar to other optimization algorithms. It involves generating an initial population of N_P individuals (where N_P is the population size) randomly within the search space of the optimization variables. Specifically, the process is as follows:

   $$\begin{array}{}\text{(25)}& \begin{array}{rl}{\mathbf{x}}_{n,d}^{\left(\mathrm{0}\right)}& ={\mathbf{x}}_{min,d}+\mathrm{rand}\left(\right)({\mathbf{x}}_{max,d}-{\mathbf{x}}_{min,d})\\ & n=\mathrm{1},\mathrm{2},\mathrm{\dots },N_P.\end{array}\end{array}$$

   In the equation, x_min,d and x_max,d represent the lower and upper bounds of the dth variable, respectively. Here, d=1, 2, …, D, where D is the total number of decision variables (the dimensionality of the search space).

2. Encircling prey

   In the BWO, since the optimal solution in the search space is not known in advance, it is assumed that the best solution obtained so far is the target prey or an approximation of the global optimum. Once the optimal solution's position is determined, the other whale individuals in the population move towards this optimal position. The position update for each whale is then given by the following equation:

   $$\begin{array}{}\text{(26)}& {\displaystyle }\mathit{\tau }=\mathrm{|}C\cdot {\mathbf{x}}^{\ast \left(t\right)}-{\mathbf{x}}^{\left(t\right)}\mathrm{|}\text{(27)}& {\displaystyle }{\mathbf{x}}^{(t+\mathrm{1})}={\mathbf{x}}^{\ast \left(t\right)}-A\cdot \mathit{\tau }.\end{array}$$

   In the equation, t represents the current iteration number, x^(t) denotes the position of the whale population at the tth iteration, x^*(t) is the position of the best solution (prey) in the whale population at iteration t, A⋅τ represents the update step size for the whale population, and A and τ are the update coefficients. Their specific expressions are as follows:

   $$\begin{array}{}\text{(28)}& \left\{\begin{array}{l}A=\mathrm{2}a{r}_{\mathrm{1}}-a\\ C=\mathrm{2}{r}_{\mathrm{2}}\end{array}\right..\end{array}$$

   In the equation, r₁ and r₂ are random numbers between [0,1]; this convergence factor is used to balance the global search ability and local search ability. It is given by the following formula: $a=\mathrm{2}-\mathrm{2}t/T$, where T is the maximum number of iterations.

3. Bubble net attack

   In the BWO, the bubble net attack is a crucial phase that combines two different movement strategies: shrink encircling and spiral updating. These two strategies allow the whales to explore the solution space globally and locally, respectively. The choice of which strategy to use is determined by a probability p, which controls the switching between the two methods. The specific operations are as follows:

   $$\begin{array}{}\text{(29)}& {\mathbf{x}}^{(t+\mathrm{1})}=\left\{\begin{array}{ll}{\mathbf{x}}^{\ast \left(t\right)}-A\cdot \mathit{\tau }& p<\mathrm{0.5}\\ {\mathit{\tau }}^{\prime }\cdot e^{bl}\cos \left(\mathrm{2}\mathit{\pi }l\right)+{\mathbf{x}}^{\ast \left(t\right)}& p\ge \mathrm{0.5}\end{array}\right..\end{array}$$

   In the equation, b is the spiral shape parameter; τ′ is the Euclidean distance between the current global best solution and the current whale's position; l is a random number uniformly distributed between [−1,1], used to introduce randomness into the movement; and p is a random number uniformly distributed between [0,1], which controls the probability of selecting a specific update method (shrink encircling or spiral updating).

4. Search for prey

   The BWO updates the positions of the whales based on the distance between the individual whales, thereby performing a random search. When $\left|A\right|\ge \mathrm{1}$, the BWO algorithm randomly selects a whale position as the target whale ${\mathbf{x}}_{\mathrm{rand}}^{\left(t\right)}$, and the other whales update their positions relative to this target in search of a better prey. The updated formula is as follows:

   $$\begin{array}{}\text{(30)}& {\mathbf{x}}^{(t+\mathrm{1})}={\mathbf{x}}_{\mathrm{rand}}^{\left(t\right)}-A\cdot {\mathit{\tau }}^{\prime \prime }.\end{array}$$

   In the equation, τ” represents the step size coefficient, ${\mathit{\tau }}^{\prime \prime }=\mathrm{|}C\cdot {\mathbf{x}}_{\mathrm{rand}}^{\left(t\right)}-{\mathbf{x}}^{\left(t\right)}\mathrm{|}$.

4.2 Improved beluga whale optimization (IBWO)

4.2.1 Chaos-inspired backward learning population initialization

The initial population of the basic differential evolution algorithm is randomly generated using Eq. (25), which can lead to an uneven distribution and clustering of the generated initial population. The use of chaotic mapping can effectively address the issue of population initialization. Commonly used chaotic maps include the Kent map (Xiang et al., 2007), Circle map (Arora et al., 2019), and Logistic map (Herbadji et al., 2019). According to the reference (Sankalap et al., 2019), the Circle map generates a more evenly distributed initial population compared to other chaotic maps. However, the initial population generated by the Circle map still tends to be denser in the range [0.2, 0.6] compared to other values. Therefore, in this paper, the Circle map is improved to achieve a relatively uniform distribution across the entire solution space. The improved Circle map formula is as follows:

$$\begin{array}{}\text{(31)}& {\mathbf{k}}_{n,D}=\mathrm{mod}\left(\mathrm{2}{\mathbf{k}}_{n-\mathrm{1},D}+\mathrm{0.35}-{\displaystyle \frac{\mathrm{0.7}}{\mathrm{2}\mathit{\pi }}}\sin \left(\mathrm{2}\mathit{\pi }{\mathbf{k}}_{n-\mathrm{1},D}\right),\mathrm{1}\right).\end{array}$$

In the equation, n=1,2,…, N_P; when n=2, $k_{n-\mathrm{1},D}=$ rand (1,D). Combining with Eq. (25), the initial population is generated using the following method:

$$\begin{array}{}\text{(32)}& {}^{\prime }{\mathbf{x}}_n={\mathbf{x}}_{min}+({\mathbf{x}}_{max}-{\mathbf{x}}_{min}){\mathbf{k}}_n.\end{array}$$

As shown in Figs. 3–5 to clearly and intuitively demonstrate the advantages of the improved Circle map, the initial population distributions for random initialization, Circle map, and improved Circle map are plotted for n=1000, along with a frequency statistical analysis.

Figure 3 Distribution diagram of rand population initialization methods showing (a) value and (b) frequency.

Figure 4 Distribution diagram of circle map population initialization methods showing (a) value and (b) frequency.

Figure 5 Distribution diagram of improved circle map population initialization methods showing (a) value and (b) frequency.

As shown in Fig. 3, the improved Circle chaotic map results in a more evenly distributed chaotic value. Therefore, the improved Circle chaotic map is used to generate the initial population, which enhances the population diversity and, in turn, improves the optimization ability of the algorithm.

Backward learning can improve the quality and diversity of the population solutions. However, the distance between the solutions generated by backward learning and the current solution is fixed, which is detrimental to population diversity and lacks a certain degree of randomness. Therefore, based on the improved Circle map, a dynamic backward learning strategy is employed to further increase population diversity and the number of elite individuals, thereby expanding the search range of the algorithm. The dynamic backward learning strategy is as follows:

$$\begin{array}{}\text{(33)}& \mathrm{"}{\mathbf{x}}_n={}^{\prime }{\mathbf{x}}_n+\mathrm{rand}\left(\mathrm{rand}\right({\mathbf{x}}_{max}+{\mathbf{x}}_{min}-{}^{\prime }{\mathbf{x}}_n)-{}^{\prime }{\mathbf{x}}_n).\end{array}$$

Combining Eqs. (32) and (34), for the minimization optimization problem, the chaos-backward learning population initialization can be expressed as

$$\begin{array}{}\text{(34)}& {\mathbf{x}}_n^{\left(\mathrm{0}\right)}=\left\{\begin{array}{ll}{}^{\prime }{\mathbf{x}}_n& f{(}^{\prime }{\mathbf{x}}_n)<f(\mathrm{"}{\mathbf{x}}_n)\\ \mathrm{"}{\mathbf{x}}_n& \mathrm{else}\end{array}\right..\end{array}$$

4.2.2 Elite strategy

According to Eq. (29), the BWO algorithm updates positions based on the current optimal position of the whale. If the current optimal whale position becomes trapped in a local optimum, all the whales will converge around the local optimal individual, causing the algorithm to stagnate. In the gray wolf optimization algorithm (GWO), to increase population diversity, a hierarchical system is proposed. This involves applying weighted processing to the three individuals with the best fitness values in the current population, using them as the current optimal position. This approach helps to avoid the risk of the population getting trapped in a local optimum.

Based on the hierarchical system of the GWO algorithm (Ghalambaz et al., 2021), this paper introduces the concept of an elite pool. Specifically, the three best individuals in the population, along with the weighted individuals of high quality, are stored in the elite pool. During the position update process, an individual is randomly selected from the elite pool to serve as a guide, enhancing the algorithm's ability to escape local optima. The elite strategy is as follows:

$$\begin{array}{}\text{(35)}& {\mathrm{Elite}}^{\left(t\right)}=[{\mathbf{x}}_{\mathrm{best},\mathrm{1}}^{\left(t\right)},{\mathbf{x}}_{\mathrm{best},\mathrm{2}}^{\left(t\right)},{\mathbf{x}}_{\mathrm{best},\mathrm{3}}^{\left(t\right)},{\mathbf{x}}_{\mathrm{mean}}^{\left(t\right)}].\end{array}$$

${\mathbf{x}}_{\mathrm{best},\mathrm{1}}^{\left(t\right)},{\mathbf{x}}_{\mathrm{best},\mathrm{2}}^{\left(t\right)},and{\mathbf{x}}_{\mathrm{best},\mathrm{3}}^{\left(t\right)}$ are the top three individuals ranked in descending order based on their fitness values. ${\mathbf{x}}_{\mathrm{mean}}^{\left(t\right)}$ is considered to denote the weighted individuals, and the calculation method is as follows:

$$\begin{array}{}\text{(36)}& \left\{\begin{array}{l}{\mathbf{x}}_{\mathrm{mean}}^{\left(t\right)}=\sum _{i=\mathrm{1}}^{\mathrm{0.2}{N}_P}{\stackrel{\mathrm{‾}}{f}}_i{\mathbf{x}}_i^{\left(t\right)}\\ {\stackrel{\mathrm{‾}}{f}}_i=\frac{f_i}{\sum _{i=\mathrm{1}}^{\mathrm{0.2}{N}_P}{f}_i}\end{array}\right..\end{array}$$

Table 1 Inverse kinematics of the 4-PUU PM.

4.2.3 Incorporating the golden sine position update strategy

The BWO algorithm uses the spiral position update strategy from Eq. (30). During the later stages of the algorithm, all individuals in the population gather near the optimal individual, which accelerates the convergence speed of the algorithm but can lead to getting stuck in local optima, resulting in poor population diversity. The golden sine algorithm (Golden-SA) is a new metaheuristic algorithm (Lu and Zhang, 2023) that has strong global search capabilities. Golden-SA introduces a golden ratio coefficient in its position update process, allowing the algorithm to thoroughly search regions that can produce excellent solutions during each iteration, thereby accelerating the convergence speed and escaping local optima. Therefore, when updating positions in the BWO algorithm, the position update strategy from Golden-SA can be adopted to enhance the global search ability and accelerate its convergence speed. Thus, the new position update strategy is as follows:

$$\begin{array}{}\text{(37)}& {\mathbf{x}}_i^{(t+\mathrm{1})}={\mathbf{x}}_i^{\left(t\right)}\mathrm{|}\sin R_{\mathrm{1}}\mathrm{|}+r_{\mathrm{1}}\sin \left(\mathrm{2}\mathit{\pi }{r}_{\mathrm{2}}\right)\mathrm{|}{\mathit{\epsilon }}_{\mathrm{1}}{\mathrm{Elite}}^{\left(t\right)}-{\mathit{\epsilon }}_{\mathrm{2}}{\mathbf{x}}_i^{\left(t\right)}\mathrm{|}.\end{array}$$

In the equation, r₁ and r₂ are random numbers between [0, 1]; R₁ is a random number between [0, 2]; and ε₁ and ε₂ are golden ratio coefficients, and their specific expressions are

$$\begin{array}{}\text{(38)}& \left\{\begin{array}{c}{\mathit{\epsilon }}_{\mathrm{1}}=-\mathit{\pi }\mathrm{\Gamma }+\mathit{\pi }(\mathrm{1}-\mathrm{\Gamma })\\ {\mathit{\epsilon }}_{\mathrm{2}}=-\mathit{\pi }(\mathrm{1}-\mathrm{\Gamma })+\mathit{\pi }\mathrm{\Gamma }\end{array}\right..\end{array}$$

In the equation, Γ is the golden ratio coefficient, and its value is $\mathrm{\Gamma }=\frac{\sqrt{\mathrm{5}}-\mathrm{1}}{\mathrm{2}}$.

Figure 6 Comparison between MATLAB theoretical calculations and adams simulation analysis results showing (a) displacement, (b) velocity, (c) acceleration, (d) angular displacement, (e) angular velocity, and (f) angular acceleration.

Figure 7 Evolution curve.

Table 2 Performance statistical results of the four algorithms.

5 Numerical experiment verification

5.1 Kinematic simulation verification

In this section, MATLAB numerical examples and Adams motion simulation are used together to verify the correctness of the kinematic equations derived for the 4-PUU PM. First, a numerical example of the mechanism's position analysis is provided. The mechanism's scale parameters are set as R=38, a=18.75, b=12.5, l=43 mm. The initial position of the fixed platform is set as (X_P, Y_P, Z_P, $\mathit{\gamma }{)}^T=(-\mathrm{20}$, 20, 20 mm, −25°)^T; the final position of the moving platform is set as (X_P, Y_P, Z_P, γ)^T=(20, −10, 50, −45°)^T. The initial and final positions of the moving platform are substituted into Eq. (13) to determine the corresponding sub-configurations of the mechanism (h_i≥^{f}Z_Ci) and the input rod lengths. The change in the driving pair of each limb is calculated as Δh_i= (65.3012, 31.5674, 13.5590, 48.1959). To avoid motion shocks caused by sudden changes in velocity and acceleration, which could affect the kinematic simulation results, the four moving driving pairs follow a 3–4–5 polynomial motion law, expressed as

$$\begin{array}{}\text{(39)}& H_i=h_i+\mathrm{\Delta }{h}_i\left(\mathrm{10}{\left({\displaystyle \frac{\mathit{\tau }}{T}}\right)}^{\mathrm{3}}-\mathrm{15}{\left({\displaystyle \frac{\mathit{\tau }}{T}}\right)}^{\mathrm{4}}+\mathrm{6}{\left({\displaystyle \frac{\mathit{\tau }}{T}}\right)}^{\mathrm{5}}\right).\end{array}$$

In the equation, H_i is the function representing the variation of the driving pair rod length of the ith limb with respect to time, and h_i is the initial rod length of the ith limb. h_i=[25.0013, 49.0150, 59.4375, 36.7018]: τ is the time variable, and its value range is [0,5]; T is the maximum simulation time.

Figure 8 IBWO algorithm error showing (a) position error and (b) pose error.

The total motion time of the mechanism is set to 5 s, and the continuous motion time T is discretized with a time interval of 0.1s. Based on this, the position forward kinematics of the mechanism is solved using the method from Sect. 3.1.2. Then, the velocity and acceleration forward kinematics are obtained using Eqs. (18) and (22). Based on this, the angular displacement, angular velocity, and angular acceleration of the moving platform under the given motion law are obtained as MATLAB theoretical numerical curves, as shown in Fig. 6. At the same time, a three-dimensional model of the mechanism is established based on the mechanism's scale parameters and assembled. The assembly is then imported into Adams software, where the driving parameters from Eq. (39) are set, resulting in the output motion curve of the moving platform, as shown in Fig. 6.

As shown in Fig. 4a and f, the mechanism has translational DOF in the X, Y, and Z directions, as well as a rotational DOF around the Z axis, thereby validating the correctness of the mechanism's DOF analysis in Sect. 2.2. At the same time, the Adams simulation results are in good agreement with the MATLAB theoretical results, confirming the correctness of the position, velocity, and acceleration equations derived for the 4-PUU PM in this paper.

5.2 Numerical experiments of the BWO algorithm

5.2.1 Single-point trajectory experiment

The mechanism scale parameters are consistent with those in Sect. 5.1. The coordinates of the end-effector reference point are (X_P, Y_P, Z_P, γ)^T= (−20, 20, 20 mm, −25°)^T. By substituting these values into Eq. (13), the inverse kinematics of the mechanism can be solved, as shown in Table 1.

For the first set of parameters in Table 1, the mechanism scale parameters are consistent with those in section 5.1. The range of the variables to be determined is ${\mathit{x}}_{\min }=(-\mathrm{50}$, −50 mm, 0, −π)T, x_max=(50, 50, 100 mm, π)T. In order to evaluate the computational performance of the IBWO algorithm in the forward kinematics of the PM, the algorithm is compared with the DEb1, PSO, and GWO algorithms. The acceptable error for all five algorithms is $f_{\mathrm{ac}}={\mathrm{10}}^{-\mathrm{8}}$, the population size is N_P=40, and the maximum number of iterations is T=500. The algorithm is independently run 50 times, and the average evolution curve of the objective function is shown in Fig. 7. Table 2 presents the results of the four algorithms for the forward kinematics problem with independent random runs. f_w, f_m, f_b, f_std, and n_s represent the worst value, average value, best value, standard deviation, and the number of times the algorithm successfully found the forward kinematics solution, respectively, after 50 independent runs of the algorithm.

From Table 2, it can be observed that the values of f_b and n_s in the IBWO algorithm are the same as those in the PSO algorithm. However, the values of f_w, f_m, and f_std obtained by the IBWO algorithm are better than those of the PSO algorithm. At the same time, as shown in Fig. 7, the IBWO algorithm exhibits faster convergence speed and greater robustness. In summary, the performance of the IBWO algorithm is superior to the other three comparison algorithms.

5.2.2 Continuous trajectory error experiment

To calculate the error of the IBWO algorithm in solving the forward kinematics equation of the 4-PUU PM, the 51 end-effector output poses from Fig. 6a are substituted into Eq. (13) to obtain the corresponding input driving rod lengths h_i. These are used as the mechanism scale parameters. The IBWO algorithm is then independently run five times to solve for the output poses. The difference between the output poses obtained by the IBWO algorithm and the theoretical parameters is considered the error of the IBWO algorithm. The specific calculation expression is as follows:

$$\begin{array}{}\text{(40)}& \mathrm{\Delta }\mathrm{Er}={\displaystyle \frac{\sum _{i=\mathrm{1}}^n(F_{\mathrm{IBWO}}-F_{\mathrm{T}})}{n}}.\end{array}$$

In the equation, ΔEr is the error generated by the IBWO algorithm in solving the forward kinematics equation of the PM, n is the number of independent runs of the IBWO algorithm, F_IBWO is the output position of the PM obtained by the IBWO algorithm, and F_T is the theoretical output pose of the PM.

Based on this, the average error curve of the IBWO algorithm in solving the forward kinematics equation of the 4-PUU PM is shown in Fig. 6.

As shown in Fig. 8a, when the IBWO algorithm is used to solve the forward kinematics equation for the continuous trajectory of the 4-PUU PM, the maximum error generated is approximately $\mathrm{2}\times {\mathrm{10}}^{-\mathrm{7}}$ mm. However, in practical engineering applications, an error of 0.01 mm is sufficient to meet the requirements. In conclusion, the roots of the forward kinematics equation of the 4-PUU PM solved by the IBWO algorithm meet the practical requirements. To improve the accuracy of the IBWO algorithm, it is only necessary to increase the acceptable root accuracy f_acc in the algorithm to achieve higher-precision roots.

6 Conclusions

This paper presents a 4-PUU PM with 3T1R DOF, whose workspace is unaffected by the actuated joints. The mechanism's DOF characteristics are analyzed using screw theory. Based on the constraint of link lengths, the forward and inverse kinematics equations of the mechanism are derived. On this basis, the first and second time derivatives of the forward kinematics equations are obtained, yielding the velocity and acceleration equations of the mechanism. Given the scale parameters and actuator parameters of the mechanism, simulation experiments are conducted using both MATLAB and Adams software. The experimental results validate the correctness of the kinematic analysis of the 4-PUU PM. An IBWO is proposed to solve the forward kinematics equations of the 4-PUU PM, and comparative algorithm experiments are performed. The experimental results show that the IBWO algorithm outperforms other comparison algorithms. Finally, based on the forward kinematics solutions for the interpolation points of continuous trajectories, it is found that the IBWO algorithm meets the accuracy requirements for the mechanism's practical applications.

In the future, the IBWO algorithm will be used for the optimization of the scale parameters of the 4-PUU PM. Based on the optimization results, a prototype will be designed, and the IBWO algorithm will be applied to the prototype for kinematic trajectory planning of the mechanism.