Design and analysis of a novel cable-driven parallel robot with variable stiffness and rigid–flexible coupling

This study presents a novel cable-driven parallel robot integrated with rigid–flexible coupling and a variable stiffness system, specifically designed to replicate the motion characteristics of the human cervical spine. Kinematic control of the robot is achieved by establishing a functional relationship between its configuration space and actuation space, where cable length variations are constrained within ±40 mm. Furthermore, this work explores static and dynamic modeling methodologies, highlighting their critical role in evaluating the robot's performance under operational poses. A detailed analysis is conducted to investigate the synergistic effects between the rigid–flexible coupling structure and the variable stiffness system. By incorporating adjustable stiffness parameters, a dynamic model is developed to reveal the robot's behavioral characteristics under different poses. Simulation results demonstrate that the robot achieves a maximum bending angle of 30° under cable tensions from 0 to 80 N, while its operational stiffness can be dynamically modulated from 3900 to 4700 N m⁻¹ (with a tuning bandwidth of up to 100 %). These quantitative findings confirm the validity of the proposed model and illustrate the capabilities of the robot in flexible bending motion and stiffness regulation, providing a valuable reference for the design and analysis of other cable-driven flexible mechanisms.

1 Introduction

With the rapid advancement of industrial technology, human–robot collaboration has become a crucial driver of technological innovation and productivity enhancement (Cao et al., 2025; Yin et al., 2026). The continuous progress in robotics has significantly liberated human labor, reducing the burden of monotonous and strenuous physical tasks. Furthermore, robotics has demonstrated extensive application potential and exceptional value in industrial manufacturing, healthcare, education, and training, as well as in specialized scenarios. As an important direction in robotics research, bionics provides a solid theoretical foundation and technical support for building flexible and efficient robotic systems by mimicking biological structures and functions (Liu et al., 2022; Wang et al., 2024). In the field of humanoid robotics, research on structural design and motion control is particularly critical (Ding et al., 2022). Specifically, bio-inspired designs for flexible human body parts, such as the human cervical spine, leverage their characteristics of multiple degrees of freedom (DOF) and high-precision motion control to exhibit outstanding advantages in performing complex tasks (Giannaccini et al., 2018). These designs open up new directions and technological pathways for the research and application of humanoid robots.

In recent years, bionic human neck robots have gained widespread attention and undergone in-depth research in the fields of bionics and robotics (Yang et al., 2020). As an essential component of the human body, the neck plays a crucial role in maintaining balance, performing complex tasks, and perceiving the external environment due to its unique flexibility and multi-degree-of-freedom motion characteristics (Kang et al., 2013). The cervical spine consists of seven vertebrae, and through the coordinated action of joints, ligaments, and muscles, it is capable of executing diverse and highly flexible movements, including flexion, extension, lateral bending, and rotation (Bhalkikar et al., 2024). Studies have shown that the unique structure of the cervical spine, along with the coordination between adjacent joints, enables precise and agile movements (Thorapalli Muralidharan et al., 2026). In biomechanics research, finite element modeling and experimental analysis have revealed the motion characteristics of the cervical spine under various loading conditions, while also highlighting the critical role of muscles and ligaments in maintaining stability during movement (Kuo et al., 2019). In the design of bionic neck robots, the currently developed structures are primarily classified into serial and parallel mechanisms. Serial mechanisms are widely applied due to their high flexibility and ease of control, whereas parallel mechanisms are favored for their higher rigidity and lower inertia, making them suitable for specific scenarios (Song et al., 2023). Parallel mechanisms offer more distinct advantages in handling high loads and low inertia, which are key factors in simulating human neck movements (Renda et al., 2012; Song et al., 2025).

The ASIMO robot by Honda and the HRP-2 robot both feature a serial neck structure with two degrees of freedom (pitch and roll), enabling basic head movements (Liu et al., 2023; Song et al., 2024). To achieve more complex motion capabilities, the Dav robot, Albert Hubo, and iCub robot have expanded their serial neck design by adding a yaw degree of freedom, further enhancing the movement abilities (Jones and Walker, 2006). The four-degrees-of-freedom serial neck design proposed in Falkenhahn builds upon the three existing degrees of freedom and adds a nodding freedom, resulting in more flexible movement (Falkenhahn et al., 2015). In contrast to serial structures, parallel neck structures typically follow the principles of parallel manipulators, allowing for pitch and roll movements, while yaw motion is achieved through additional mechanisms (Trivedi et al., 2008). For example, the design of the SAYA robot combines a spring spine with pneumatic artificial muscles, bringing head movements closer to the functionality of human muscles, though the sensitivity still needs improvement (Camarillo et al., 2008). The iCub robot has further explored parallel neck structures, including a combination of a spring spine and cables, as well as a parallel manipulator with three degrees of freedom. The latter design enhances stability through a central passive spherical support. The head structure of the James humanoid robot is similar to the first parallel neck design of the iCub robot. Additionally, the literature introduces a parallel neck design based on cable-driven elastic limbs and passive body parts (Falkenhahn et al., 2017; Renda et al., 2018). This innovative structural layout enhances the performance and adaptability of the bionic neck. However, the design still lacks sufficient compliance to simulate the bending motion of the cervical spine, failing to fully replicate the motion characteristics of the human neck.

In summary, while existing neck robot designs have achieved improvements in motion performance and adaptability, they are still constrained by inherent limitations. Parallel neck designs enhance motion capabilities but fail to adequately address the compliance required for bending motions, preventing them from fully replicating the flexibility and multi-degree-of-freedom movements of the human cervical spine. Although some designs have increased the number of degrees of freedom and introduced greater complexity, issues of sensitivity and biomimicry remain unresolved.

Figure 1 Design framework of the cable-driven, rigid–flexible coupled variable stiffness parallel robot for the cervical spine.

Based on an in-depth analysis of biomimetic human cervical robots and key functional requirements, this paper proposes a novel stiff–flexible coupling variable stiffness parallel robot (as shown in Fig. 1). The robot uses a cable-driven unit that combines rigid and flexible components (Sadati et al., 2021; Yang et al., 2021). It also incorporates a variable stiffness system that simulates the muscle adjustment function of the human cervical spine. This design allows the robot to more closely resemble the motion of the human cervical spine and meets human anatomical movement requirements. Compared to the cable-driven parallel robot with a spring-based spine designed by Gao and Song, the proposed robot features an evaluated workspace distribution for its rigid–flexible coupling mechanism (Grazioso et al., 2019). By applying torque balance theory based on Newton's second law and the conservation of angular momentum, it redefines the range of bending and azimuth angles (Luo et al., 2024; Yang et al., 2024). This improvement allows for a more accurate simulation of the motion characteristics of the human cervical spine. In addition, inspired by the working mechanisms of human cervical muscles and ligaments, the robot features an innovative rotational structure (Godage et al., 2015). This structure adjusts the radius of the external flexible support system through rotational motion. As a result, it dynamically regulates the cable driving force and alters the robot's overall stiffness characteristics. Consequently, the neck structure proposed in this study seeks to optimize both motion freedom and compliance, improve sensitivity, and more accurately mimic the natural movement of the human cervical spine, with the goal of overcoming the limitations inherent in current designs.

The structure of this paper is arranged as follows: Sect. 2 provides a brief introduction to the overall structure of the stiff–flexible coupling variable stiffness parallel robot and analyzes its general kinematic model; Sect. 3 delves into the static and dynamic characteristics of the robot and establishes a classical dynamics model; Sect. 4 verifies the effectiveness of the design through simulations and evaluates the stiff–flexible coupling and variable stiffness system; Sect. 5 presents a critical analysis and discussion of the results; finally, Sect. 6 summarizes the main research findings of this paper.

2 Configuration design and kinematic analysis of the parallel robot

2.1 Configuration design of the parallel robot

Figure 2 illustrates the configuration design of a cable-driven, rigid–flexible coupling parallel robot with variable stiffness. The robot comprises five main components: the moving platform, the fixed platform, the six-bar spring system, the rotating mechanism, and the external flexible support system. It is driven by three cables with negligible mass and diameter. The structural design incorporates a rigid–flexible coupling structure, where the six-bar spring system mimics the function of the human cervical spine, and the external flexible support system simulates the neck muscles, improving the robot's compliance and adaptability during motion. In operation, the moving platform is driven by cables, while the rotating mechanism adjusts the radius of the external flexible support system to simulate the force exerted by the neck muscles, thereby enabling efficient modulation of the robot's stiffness (covering kinematic analysis, dynamic analysis, and simulation analysis).

Figure 2 Configuration design of the parallel robot.

2.2 General kinematic analysis of the parallel robot

The rigid–flexible coupling parallel robot with variable stiffness has three DOF, enabling bending and stretching motions similar to the human cervical spine. Cables u (where u=1, 2, 3) pass through points C₁, C₂, and C₃ on the fixed platform and are connected to cable holes on the moving platform, as shown in Fig. 2. The robot's motion is driven by motors that control the cables. Note that although the physical prototype can theoretically execute 3-DOF motion, in this study, the distance l₀ is assumed constant to focus on the 2-DOF bending performance.

In Fig. 3, the six-bar spring system ensures structural stability between the fixed and moving platforms and provides necessary bending compliance. The system is fixed to a rotating structure, with an angle adjustment range of 0 to 60°. A motor drives the transmission components at the bottom of the rotating structure, adjusting the variable radius ^βr_Bz of the external flexible support system to control stiffness.

Figure 3 Assembly of the six-bar spring system and the external flexible support system.

This design allows the robot to exhibit different stiffness levels in various operating environments, enhancing its adaptability.

Figures 2 and 3 illustrate the global coordinate system established on the fixed platform, with its origin O located at the geometric center of the platform. The radius of the cable hole is denoted by r_c, and the radius of the installation point for the six-bar spring system is denoted by r_z. These parameters allow for the precise determination of the coordinates of each installation point within the coordinate system. This serves as the foundation for subsequent kinematic and dynamic analyses. The establishment of this coordinate system is crucial for understanding the relationship between the parallel robot platform and its connected components, as well as for analyzing the mechanical behavior under different operating conditions.

$$\begin{array}{}\text{(1)}& {\displaystyle }\left\{\begin{array}{l}{C}_{\mathrm{1}}={\left[\begin{array}{ccc}\frac{\sqrt{\mathrm{3}}}{\mathrm{2}{r}_{\mathrm{c}}}& \frac{r_{\mathrm{c}}}{\mathrm{2}}& \mathrm{0}\end{array}\right]}^{\mathrm{T}}\\ C_{\mathrm{2}}={\left[\begin{array}{ccc}-\frac{\sqrt{\mathrm{3}}}{\mathrm{2}{r}_{\mathrm{c}}}& \frac{r_{\mathrm{c}}}{\mathrm{2}}& \mathrm{0}\end{array}\right]}^{\mathrm{T}}\\ C_{\mathrm{3}}={\left[\begin{array}{ccc}\mathrm{0}& -r_{\mathrm{c}}& \mathrm{0}\end{array}\right]}^{\mathrm{T}}\end{array}\right.\text{(2)}& {\displaystyle }\left\{\begin{array}{l}{}^{\mathrm{1}}{\mathbf{D}}_{z\mathrm{1}}={\left[\begin{array}{ccc}{r}_z\sin \mathit{\alpha }& r_z\cos \mathit{\alpha }& \mathrm{0}\end{array}\right]}^{\mathrm{T}}\\ {}^{\mathrm{1}}{\mathbf{D}}_{z\mathrm{2}}={\left[\begin{array}{ccc}-r_z\sin \mathit{\alpha }& r_z\cos \mathit{\alpha }& \mathrm{0}\end{array}\right]}^{\mathrm{T}}\\ {}^{\mathrm{2}}{\mathbf{D}}_{z\mathrm{1}}={\left[\begin{array}{ccc}-\frac{r_z}{\mathrm{2}}\left(\sqrt{\mathrm{3}}\cos \mathit{\alpha }+\sin \mathit{\alpha }\right)& \frac{r_z}{\mathrm{2}}\left(\sqrt{\mathrm{3}}\sin \mathit{\alpha }-\cos \mathit{\alpha }\right)& \mathrm{0}\end{array}\right]}^{\mathrm{T}}\\ {}^{\mathrm{2}}{\mathbf{D}}_{z\mathrm{2}}={\left[\begin{array}{ccc}\frac{r_z}{\mathrm{2}}\left(\sin \mathit{\alpha }-\sqrt{\mathrm{3}}\cos \mathit{\alpha }\right)& -\frac{r_z}{\mathrm{2}}\left(\cos \mathit{\alpha }+\sqrt{\mathrm{3}}\sin \mathit{\alpha }\right)& \mathrm{0}\end{array}\right]}^{\mathrm{T}}\\ {}^{\mathrm{3}}{\mathbf{D}}_{z\mathrm{1}}={\left[\begin{array}{ccc}\frac{r_z}{\mathrm{2}}\left(\sqrt{\mathrm{3}}\cos \mathit{\alpha }+\sin \mathit{\alpha }\right)& \frac{r_z}{\mathrm{2}}\left(\sqrt{\mathrm{3}}\sin \mathit{\alpha }-\cos \mathit{\alpha }\right)& \mathrm{0}\end{array}\right]}^{\mathrm{T}}\\ {}^{\mathrm{3}}{\mathbf{D}}_{z\mathrm{2}}={\left[\begin{array}{ccc}\frac{r_z}{\mathrm{2}}\left(\sqrt{\mathrm{3}}\cos \mathit{\alpha }-\sin \mathit{\alpha }\right)& -\frac{r_z}{\mathrm{2}}\left(\cos \mathit{\alpha }+\sqrt{\mathrm{3}}\sin \mathit{\alpha }\right)& \mathrm{0}\end{array}\right]}^{\mathrm{T}}\end{array}\right.\end{array}$$

Figure 4 Schematic of the fixed platform structure.

As shown in Fig. 4, the variable stiffness system adjusts flexibility by modifying the radius δr_z of the external flexible support system. The variable radius δr_z is determined by the rotation angle β of the rotating structure through the function f_rotate, and combined with the fixed radius r_B, the actual radius r_Bz of the external flexible support system is calculated, allowing the coordinates of the three fixed points of the system to be determined.

$$\begin{array}{}\text{(3)}& \mathit{\delta }{r}_z=\left\{\begin{array}{l}\mathrm{0}=f_{\mathrm{rotate}}\left(\mathit{\beta }=\mathrm{0}\right)\\ f_{\mathrm{rotate}}\left(\mathit{\beta }\right){|}_{\mathrm{0}}^{\mathrm{\Delta }\mathit{\beta }}\end{array}\right\}=\mathit{\delta }{r}_z\left(\mathit{\beta }\right),\mathrm{\Delta }\mathit{\beta }\in \left(\mathrm{0},{\displaystyle \frac{\mathit{\pi }}{\mathrm{3}}}\right)\end{array}$$

$$\begin{array}{}\text{(4)}& \begin{array}{rl}& {}^{\mathit{\beta }}{r}_{Bz}=r_B+\mathit{\delta }{r}_z\left(\mathit{\beta }\right)\\ & \left\{\begin{array}{l}{}^k{\mathbf{B}}_{\mathrm{1}}={\left[\begin{array}{clc}\mathrm{0}& {}^{\mathit{\beta }}{r}_{Bz}& \mathrm{0}\end{array}\right]}^{\mathrm{T}}\\ {}^k{\mathbf{B}}_{\mathrm{2}}={\left[\begin{array}{ccc}-{\frac{\sqrt{\mathrm{3}}}{\mathrm{2}}}^{\mathit{\beta }}{r}_{Bz}& -{\frac{\mathrm{1}}{\mathrm{2}}}^{\mathit{\beta }}{r}_{Bz}& \mathrm{0}\end{array}\right]}^{\mathrm{T}}\\ {}^k{\mathbf{B}}_{\mathrm{3}}={\left[\begin{array}{ccc}{\frac{\sqrt{\mathrm{3}}}{\mathrm{2}}}^{\mathit{\beta }}{r}_{Bz}& -{\frac{\mathrm{1}}{\mathrm{2}}}^{\mathit{\beta }}{r}_{Bz}& \mathrm{0}\end{array}\right]}^{\mathrm{T}}\end{array}\right.\end{array}\end{array}$$

Here, δr_z is the variable radius of the slide rail and ^kB₁, ^kB₂, and ^kB₃ are the external flexible support system fixed point coordinates.

Equations (1), (2), and (4) describe the coordinates of the cable hole, the installation points of the six-bar spring system, and the fixed points of the external flexible support system, respectively. Based on the posture of the parallel robot shown in Fig. 3, setting the bending angle θ_i and azimuth angle ϕ_i allows the configuration space q_i to be determined.

$$\begin{array}{}\text{(5)}& q_i=[{\mathit{\varphi }}_i,{\mathit{\theta }}_i{]}^{\mathrm{T}}\end{array}$$

From the state of the configuration space q_i, the actuation space ^Δl of the parallel robot can then be derived.

$$\begin{array}{}\text{(6)}& {}^{\mathrm{\Delta }}l={[}^{\mathrm{\Delta }}{L}_{\mathrm{cable}\mathrm{1}}{,}^{\mathrm{\Delta }}{L}_{\mathrm{cable}\mathrm{2}}{,}^{\mathrm{\Delta }}{L}_{\mathrm{cable}\mathrm{3}}{]}^{\mathrm{T}}\end{array}$$

The configuration space q_i and the actuation space ^Δl are related through the homogeneous transformation matrix ^RT_P (Liu et al., 2020). Based on the configuration space definition, the rotation matrix R_q (ϕ_i, θ_i) and the actuator end position P_q (ϕ_i, θ_i) can be determined. When the actuation space value is zero (i.e., the length change of all three cables is zero), both the bending angle θ_i and the azimuth angle ϕ_i are zero, indicating that the parallel robot's driving cables are pre-tensioned, as illustrated in Fig. 2. This represents the initial pre-tensioned state of the robot, where the distance between the moving platform and the fixed platform is l_o. To establish a simplified yet effective kinematic model, it is assumed that the three actuation cables are controlled synchronously during main bending motions to maintain constant overall length along the central coupling axis, ensuring that the distance between P_q and point O remains constant at l_o. The axial length variation of the flexible spring elements is primarily utilized for internal stiffness modulation rather than primary kinematic positioning.

During motion, the actuator end P_q (ϕ_i, θ_i) moves along the arc trajectory Arc(tt) shown in Fig. 2, due to the combined effect of the cable-driven six-bar spring system and the external flexible support system. The arc's center is located at the origin O of the global coordinate system, ensuring that the distance between P_q and point O remains constant at l_o. From this characteristic, the mathematical expressions for the rotation matrix R_q (ϕ_i, θ_i) and the actuator end position P_q (ϕ_i, θ_i) can be derived (Ma et al., 2021).

$$\begin{array}{}\text{(7)}& {\displaystyle }\begin{array}{rl}& {\mathbf{R}}_q\left(q_i\right)={\mathbf{R}}_q({\mathit{\varphi }}_i,{\mathit{\theta }}_i)={\left[\begin{array}{ccc}{a}_{\mathrm{1}}& a_{\mathrm{2}}& a_{\mathrm{3}}\\ b_{\mathrm{1}}& b_{\mathrm{2}}& b_{\mathrm{3}}\\ c_{\mathrm{1}}& c_{\mathrm{2}}& c_{\mathrm{3}}\end{array}\right]}_{\mathrm{3}\times \mathrm{3}}\\ & \left\{\begin{array}{l}{a}_{\mathrm{1}}=\left(\cos \right({\mathit{\theta }}_i)-\mathrm{1}){\cos }^{\mathrm{2}}\left({\mathit{\varphi }}_i\right)+\mathrm{1}\\ a_{\mathrm{2}}=b_{\mathrm{1}}=\frac{\mathrm{1}}{\mathrm{2}}\sin \left(\mathrm{2}{\mathit{\varphi }}_i\right)\left(\cos \left({\mathit{\theta }}_i\right)-\mathrm{1}\right)\\ a_{\mathrm{3}}=-c_{\mathrm{1}}=\cos \left({\mathit{\varphi }}_i\right)\sin \left({\mathit{\theta }}_i\right)\\ b_{\mathrm{2}}=\left(\cos \right({\mathit{\theta }}_i)-\mathrm{1}){\sin }^{\mathrm{2}}\left({\mathit{\varphi }}_i\right)+\mathrm{1}\\ b_{\mathrm{3}}=-c_{\mathrm{2}}=\sin \left({\mathit{\varphi }}_i\right)\sin \left({\mathit{\theta }}_i\right)\\ c_{\mathrm{3}}=\cos \left({\mathit{\theta }}_i\right)\end{array}\right.\end{array}\text{(8)}& {\displaystyle }{\mathbf{P}}_q\left(q_i\right)={\mathbf{P}}_q({\mathit{\varphi }}_i,{\mathit{\theta }}_i)=\left[\begin{array}{c}\sin {\mathit{\theta }}_i\cos {\mathit{\varphi }}_i\\ \sin {\mathit{\theta }}_i\sin {\mathit{\varphi }}_i\\ \cos {\mathit{\theta }}_i\end{array}\right]l_{\mathrm{o}}\end{array}$$

Here, R_q is rotation matrix and P_q is actuator end matrix.

Based on the configuration space q_i, the rotation matrix between the moving platform and the fixed platform, as well as the expression for the actuator end position in the global coordinate system, can be obtained from Eqs. (6) and (7) (Shao, 2015). By substituting these expressions into Eqs. (7) and (8), the homogeneous transformation matrix ^RT_P, which describes the relationship between the moving platform and the fixed platform, can be derived.

$$\begin{array}{}\text{(9)}& {}^{\mathbf{R}}{\mathbf{T}}_{\mathbf{P}}=\left[\begin{array}{cc}{\mathbf{R}}_q({\mathit{\varphi }}_i,{\mathit{\theta }}_i)& {\mathbf{P}}_q({\mathit{\varphi }}_i,{\mathit{\theta }}_i)\\ \mathrm{0}& \mathrm{1}\end{array}\right]\end{array}$$

Here, ^RT_P is homogeneous transformation matrix.

When the parallel robot is in its initial state, the actuation space state is zero. At this point, based on the definition of Eq. (8), the homogeneous transformation matrix ${}^{\mathbf{R}}{\mathbf{T}}_{\mathbf{P}}^o$ takes on a specific form.

$$\begin{array}{}\text{(10)}& {}^{\mathbf{R}}{{\mathbf{T}}_{\mathbf{P}}}^o=\left[\begin{array}{cccc}\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& l_{\mathrm{o}}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right]{(}^{\mathrm{\Delta }}l=\mathrm{0})\end{array}$$

Here, ^RT_P is initial state homogeneous transition matrix.

Using the homogeneous transformation matrix ^RT_P from Eq. (9), the actuation space ^Δl of the parallel robot can be further calculated. By controlling the changes in cable length through ^Δl, the motion of the parallel robot can be driven. Since ^RT_P is a 4 × 4 matrix, when converting the point positions on the platform to the global coordinate system, the hole position coordinates in Eqs. (1), (2), and (4) need to be expressed in a 4 × 1 vector format, as shown in Eq. (11).

$$\begin{array}{}\text{(11)}& {\stackrel{\mathrm{̃}}{\mathbf{C}}}_{\mathrm{1}}=[{\mathbf{C}}_{\mathrm{1}};\mathrm{1}];{\stackrel{\mathrm{̃}}{\mathbf{C}}}_{\mathrm{2}}=[{\mathbf{C}}_{\mathrm{2}};\mathrm{1}];{\stackrel{\mathrm{̃}}{\mathbf{C}}}_{\mathrm{3}}=[{\mathbf{C}}_{\mathrm{3}};\mathrm{1}].\end{array}$$

A function G is defined to describe the mapping relationship between the configuration space q_i and the actuation space ^Δl. This function converts parameters from the configuration space into changes in the actuation space. This enables precise control of the parallel robot's driving cable lengths. By setting the parameters of the configuration space q_i and applying function G, the changes in cable length can be determined. By accurately controlling these changes, the desired motion posture of the parallel robot can be adjusted.

$$\begin{array}{}\text{(12)}& \begin{array}{rl}\mathbf{G}\left(q_i\right)& =\left[\begin{array}{c}{}^{\mathrm{\Delta }}{L}_{\mathrm{cable}\mathrm{1}}\\ {}^{\mathrm{\Delta }}{L}_{\mathrm{cable}\mathrm{2}}\\ {}^{\mathrm{\Delta }}{L}_{\mathrm{cable}\mathrm{3}}\end{array}\right]\\ & =\left[\begin{array}{c}\Vert \left({}^{\mathbf{R}}{\mathbf{T}}_{\mathbf{P}}-E\right){\stackrel{\mathrm{̃}}{\mathbf{C}}}_{\mathrm{1}}\Vert -l_{\mathrm{o}}\\ \Vert \left({}^{\mathbf{R}}{\mathbf{T}}_{\mathbf{P}}-E\right){\stackrel{\mathrm{̃}}{\mathbf{C}}}_{\mathrm{2}}\Vert -l_{\mathrm{o}}\\ \Vert \left({}^{\mathbf{R}}{\mathbf{T}}_{\mathbf{P}}-E\right){\stackrel{\mathrm{̃}}{\mathbf{C}}}_{\mathrm{3}}\Vert -l_{\mathrm{o}}\end{array}\right]\left(q_i\right)\end{array}\end{array}$$

The function G in Eq. (12) provides a new method for analyzing the kinematics of the parallel robot, especially for the kinematic modeling of bio-inspired human cervical parallel robots. Using this function G, the mapping function A from the actuation space ^Δl to the configuration space q_i can be further derived. By precisely controlling the cables, the parallel robot's motion posture within the workspace can be adjusted, ensuring that the actuator's end effector reaches the designated target point.

$$\begin{array}{}\text{(13)}& q_i({\mathit{\varphi }}_i,{\mathit{\theta }}_i)=\mathbf{A}{(}^{\mathrm{\Delta }}l)={\mathbf{G}}^{-\mathrm{1}}(q_i)\end{array}$$

3 Dynamic analysis of the parallel robot

3.1 Static analysis

Based on the kinematics of the rigid–flexible coupling variable stiffness parallel robot, the specific posture of the robot during motion can be determined through the mapping function G from the configuration space q_i to the actuation space ^Δl. When the robot reaches a specific posture, its static behavior can be analyzed to examine the forces acting on the components during motion. When the robot is in the bent state shown in Fig. 5 and remains stationary, the driving torque applied by the cables balances the resistive torque generated by the six-bar spring system and the external flexible support system. By applying the principles of torque and force equilibrium (based on Newton's laws of motion and the conservation of angular momentum) and force equilibrium theory (derived from Newton's first and second laws), a detailed analysis of the forces on the cables in the stationary state can be conducted.

Figure 5 Bending posture and torque analysis of the cable-driven parallel robot.

By combining the homogeneous transformation matrix defined in Eq. (9) with the coordinates of each point in Eqs. (1), (2), and (4), the lengths of the driving cables, the six-bar spring system, and the external flexible support system can be calculated. Additionally, as the robot moves, the length variations of these systems can be obtained for further analysis of their dynamic response characteristics.

$$\begin{array}{}\text{(14)}& \left\{\begin{array}{l}{\mathit{\delta }}_x{\mathbf{L}}^{\mathrm{a}}=\left[\begin{array}{c}∥{}^{\mathbf{R}}{\mathbf{T}}_{\mathbf{P}}^{\mathrm{3}}{\stackrel{\mathrm{̃}}{\mathbf{D}}}_{z\mathrm{2}}{-}^{\mathrm{1}}{\stackrel{\mathrm{̃}}{\mathbf{D}}}_{z\mathrm{1}}∥\\ ∥{}^{\mathbf{R}}{\mathbf{T}}_{\mathbf{P}}^{\mathrm{1}}{\stackrel{\mathrm{̃}}{\mathbf{D}}}_{z\mathrm{1}}{-}^{\mathrm{3}}{\stackrel{\mathrm{̃}}{\mathbf{D}}}_{z\mathrm{2}}∥\end{array}\right]-d_{\mathrm{o}}{\left[\begin{array}{cc}\mathrm{1}& \mathrm{1}\end{array}\right]}^{\mathrm{T}}\\ \Rightarrow \left[\begin{array}{cc}∥{\mathit{\delta }}^{\mathrm{a}}{x}_{d\mathrm{1}}∥& \mathrm{0}\\ \mathrm{0}& ∥{\mathit{\delta }}^{\mathrm{a}}{x}_{d\mathrm{2}}∥\end{array}\right]-d_{\mathrm{o}}\mathbf{E}\\ {\mathit{\delta }}_x{\mathbf{L}}^{\mathrm{b}}=\left[\begin{array}{c}∥{}^{\mathbf{R}}{\mathbf{T}}_{\mathbf{P}}^{\mathrm{1}}{\stackrel{\mathrm{̃}}{\mathbf{D}}}_{z\mathrm{2}}{-}^{\mathrm{2}}{\stackrel{\mathrm{̃}}{\mathbf{D}}}_{z\mathrm{1}}∥\\ ∥{}^{\mathbf{R}}{\mathbf{T}}_{\mathbf{P}}^{\mathrm{2}}{\stackrel{\mathrm{̃}}{\mathbf{D}}}_{z\mathrm{1}}{-}^{\mathrm{1}}{\stackrel{\mathrm{̃}}{\mathbf{D}}}_{z\mathrm{2}}∥\end{array}\right]-d_{\mathrm{o}}{\left[\begin{array}{cc}\mathrm{1}& \mathrm{1}\end{array}\right]}^{\mathrm{T}}\\ \Rightarrow \left[\begin{array}{cc}∥{\mathit{\delta }}^{\mathrm{b}}{x}_{d\mathrm{1}}∥& \mathrm{0}\\ \mathrm{0}& ∥{\mathit{\delta }}^{\mathrm{b}}{x}_{d\mathrm{2}}∥\end{array}\right]-d_{\mathrm{o}}\mathbf{E}\left(d_{\mathrm{o}}=\frac{l_{\mathrm{o}}}{\cos \mathit{\sigma }}\right)\\ {\mathit{\delta }}_x{L}^{\mathrm{c}}=\left[\begin{array}{c}∥{}^{\mathbf{R}}{\mathbf{T}}_{\mathbf{P}}^{\mathrm{2}}{\stackrel{\mathrm{̃}}{\mathbf{D}}}_{z\mathrm{2}}{-}^{\mathrm{3}}{\stackrel{\mathrm{̃}}{\mathbf{D}}}_{z\mathrm{1}}∥\\ ∥{}^{\mathbf{R}}{\mathbf{T}}_{\mathbf{P}}^{\mathrm{3}}{\stackrel{\mathrm{̃}}{\mathbf{D}}}_{z\mathrm{1}}{-}^{\mathrm{2}}{\stackrel{\mathrm{̃}}{\mathbf{D}}}_{z\mathrm{2}}∥\end{array}\right]-d_{\mathrm{o}}{\left[\begin{array}{cc}\mathrm{1}& \mathrm{1}\end{array}\right]}^{\mathrm{T}}\\ \Rightarrow \left[\begin{array}{cc}∥{\mathit{\delta }}^{\mathrm{c}}{x}_{d\mathrm{1}}∥& \mathrm{0}\\ \mathrm{0}& ∥{\mathit{\delta }}^{\mathrm{c}}{x}_{d\mathrm{2}}∥\end{array}\right]-d_{\mathrm{o}}\mathbf{E}\end{array}\right.\end{array}$$

Here, d_o is initial length of six-bar spring system.

In Eq. (14), when the parallel robot is in the zero actuation space state, d_o represents the initial length of each push rod in the six-bar spring system, and σ is the angle between the spring system's installation position and the midline, as shown in Fig. 4. The six-bar spring system consists of three intersecting structures, each comprising two spring push rods. The three sets of crossed structures are evenly distributed between the fixed platform and the moving platform, with the center axis at O–P_q, and each pair of structures forms an angle of 120° with each other. The length changes are denoted as δ_xL^a, δ_xL^b, and δ_xL^c, respectively. During bending, the real-time lengths of the three intersecting structures are represented as δ^ϰx_d1 and δ^ϰx_d2 (where ϰ = a, b, c).

$$\begin{array}{}\text{(15)}& \left\{\begin{array}{ll}{\mathit{\delta }}^Y{\mathbf{B}}_{\mathrm{\Delta }}^{\mathrm{all}}& =-l_{\mathrm{o}}+\left[\begin{array}{c}\Vert \mathit{\delta }{B}_{\mathrm{\Delta }}^{\mathrm{a}}\Vert \\ \Vert \mathit{\delta }{B}_{\mathrm{\Delta }}^{\mathrm{b}}\Vert \\ \Vert \mathit{\delta }{B}_{\mathrm{\Delta }}^{\mathrm{c}}\Vert \end{array}\right]\\ & =-l_{\mathrm{o}}\left[\begin{array}{c}\mathrm{1}\\ \mathrm{1}\\ \mathrm{1}\end{array}\right]+\left[\begin{array}{c}∥{\left({}^{\mathbf{R}}{\mathbf{T}}_{\mathbf{P}}-E\right)}^k{\stackrel{\mathrm{̃}}{B}}_{\mathrm{1}}∥\\ ∥{\left({}^{\mathbf{R}}{\mathbf{T}}_{\mathbf{P}}-E\right)}^k{\stackrel{\mathrm{̃}}{B}}_{\mathrm{2}}∥\\ ∥{\left({}^{\mathbf{R}}{\mathbf{T}}_{\mathbf{P}}-E\right)}^k{\stackrel{\mathrm{̃}}{B}}_{\mathrm{3}}∥\end{array}\right]\left(\mathit{\delta }{r}_z\right(\mathit{\beta })=\mathrm{0})\\ {\mathit{\delta }}^N{\mathbf{B}}_{\mathrm{\Delta }}^{\mathrm{all}}& =\left[\begin{array}{c}\Vert \mathit{\delta }{B}_{\mathrm{\nabla }}^{\mathrm{a}}\Vert \\ \Vert \mathit{\delta }{B}_{\mathrm{\nabla }}^{\mathrm{b}}\Vert \\ \Vert \mathit{\delta }{B}_{\mathrm{\nabla }}^{\mathrm{c}}\Vert \end{array}\right]-l_{\mathrm{o}}\\ & =\left[\begin{array}{l}∥{\left({}^{\mathbf{R}}{\mathbf{T}}_{\mathbf{P}}-E\right)}^k{\stackrel{\mathrm{̃}}{B}}_{\mathrm{1}}∥-l_{\mathrm{o}}\\ ∥{\left({}^{\mathbf{R}}{\mathbf{T}}_{\mathbf{P}}-E\right)}^k{\stackrel{\mathrm{̃}}{B}}_{\mathrm{2}}∥-l_{\mathrm{o}}\\ ∥{\left({}^{\mathbf{R}}{\mathbf{T}}_{\mathbf{P}}-E\right)}^k{\stackrel{\mathrm{̃}}{B}}_{\mathrm{3}}∥-l_{\mathrm{o}}\end{array}\right]\left(\mathit{\delta }{r}_z\right(\mathit{\beta })\ne \mathrm{0})\end{array}\right.\end{array}$$

Here, ${\mathit{\delta }}^Y{B}_{\mathrm{\Delta }}^{\mathrm{all}}$ is the length variation of the external flexible support system when δr_z(β)=0 and ${\mathrm{\Delta }}^N{B}_{\mathrm{\Delta }}^{\mathrm{all}}$ is the length variation of the external flexible support system when δr_z(β)≠0.

Equation (15) describes the deformation of the external flexible support system in the bent state. Using Eq. (4) and the transformation in Eq. (11), the coordinate vector satisfying the homogeneous transformation matrix is obtained. By incorporating the value of the variable radius δr_z(β), the calculation formula for the deformation of the external flexible support system is derived. When δr_z(β)=0, the parallel robot cannot adjust stiffness during operation. However, when δr_z(β)≠0, the stiffness of the external flexible support system can be adjusted to achieve dynamic stiffness modulation of the robot. When the parallel robot is in a static state, the forces acting on the three driving cables are denoted as $F_n^{\mathrm{C}}$ (n=1, 2, 3). The vertical components of these forces are denoted as ${}^Z{F}_n^{\mathrm{C}}$ (n=1, 2, 3). Using force and moment balance theory, the static equilibrium of the parallel robot can be analyzed. The force balance equations for the parallel robot are established, where the forces in the vertical direction include the gravity of the moving platform, the tension in the driving cables, and the deformation forces from the six-bar spring system and the external flexible support system. By setting up the vector equation of forces, the force relationships for the three driving cables can be determined.

$$\begin{array}{}\text{(16)}& \begin{array}{rl}& \left.\begin{array}{l}{}^{\mathrm{out}}{k}_{\mathrm{spri}}\sum _{k=\mathrm{1}}^{\mathrm{3}}{\mathit{\delta }}^Y{B}_{\mathrm{\Delta }}^{\mathrm{all}}\left(k\right)\\ \mathrm{OR}\\ {}^{\mathrm{out}}{k}_{\mathrm{spri}}\sum _k^{\mathrm{3}}{\mathit{\delta }}^N{B}_{\mathrm{\Delta }}^{\mathrm{all}}\left(k\right)\end{array}\right\}+k_{\mathrm{rod}}\sum _{k=\mathrm{1}}^{\mathrm{2}}\left[{\mathit{\delta }}_x{L}^{\mathrm{a}}\left(k\right)+{\mathit{\delta }}_x{L}^{\mathrm{b}}\left(k\right)+{\mathit{\delta }}_x{L}^{\mathrm{c}}\left(k\right)\right]\\ & =\left[\begin{array}{c}\mathrm{\dots }\\ \mathrm{\dots }\\ F_{\mathrm{o}/\mathrm{s}}^Z\end{array}\right]{\left[\left[\begin{array}{c}{F}_{\mathrm{1}}^{\mathrm{C}}\\ F_{\mathrm{2}}^{\mathrm{C}}\\ F_{\mathrm{3}}^{\mathrm{C}}\end{array}\right]\left[\begin{array}{c}\frac{{}^{\mathbf{R}}{\mathbf{T}}_{\mathbf{P}}^{-\mathrm{1}}\left({}^{\mathbf{R}}{\mathbf{T}}_{\mathbf{P}}-E\right){\stackrel{\mathrm{̃}}{C}}_{\mathrm{1}}}{∥\left({}^{\mathbf{R}}{\mathbf{T}}_{\mathbf{P}}-E\right){\stackrel{\mathrm{̃}}{C}}_{\mathrm{1}}∥}\\ \frac{{}^{\mathbf{R}}{\mathbf{T}}_{\mathbf{P}}^{-\mathrm{1}}\left({}^{\mathbf{R}}{\mathbf{T}}_{\mathbf{P}}-E\right){\stackrel{\mathrm{̃}}{C}}_{\mathrm{2}}}{∥\left({}^{\mathbf{R}}{\mathbf{T}}_{\mathbf{P}}-E\right){\stackrel{\mathrm{̃}}{C}}_{\mathrm{2}}∥}\\ \frac{{}^{\mathbf{R}}{\mathbf{T}}_{\mathbf{P}}^{-\mathrm{1}}\left({}^{\mathbf{R}}{\mathbf{T}}_{\mathbf{P}}-E\right){\stackrel{\mathrm{̃}}{C}}_{\mathrm{3}}}{∥\left({}^{\mathbf{R}}{\mathbf{T}}_{\mathbf{P}}-E\right){\stackrel{\mathrm{̃}}{C}}_{\mathrm{3}}∥}\end{array}\right]\right]}^{\mathrm{T}}\\ & =\left[\left[\begin{array}{l}\mathrm{\dots }\\ \mathrm{\dots }\\ {}^z{F}_{\mathrm{1}}^{\mathrm{C}}\\ \mathrm{\dots }\end{array}\right]\left[\begin{array}{l}\mathrm{\dots }\\ \mathrm{\dots }\\ {}^z{F}_{\mathrm{2}}^{\mathrm{C}}\\ \mathrm{\dots }\end{array}\right]\left[\begin{array}{l}\mathrm{\dots }\\ \mathrm{\dots }\\ {}^z{F}_{\mathrm{3}}^{\mathrm{C}}\\ \mathrm{\dots }\end{array}\right]\right]\sum _{u=\mathrm{1}}^{\mathrm{3}}{}^Z{F_u}^{\mathrm{C}}+m_i{g}_{\mathrm{true}}={F_{\mathrm{o}/\mathrm{s}}}^Z\end{array}\end{array}$$

Here, m_i is the mass of the movement platform and g_true is the acceleration of gravity.

Equation (16) represents the force balance equation of the parallel robot in a static posture. It includes the forces $F_{\mathrm{o}/\mathrm{s}}^Z$ provided by both systems in the vertical direction, as well as the vertical components of the driving cable forces ${}^Z{F}_u^{\mathrm{C}}$ along the Z axis of the global coordinate system. Depending on the adjustable radius of the external flexible support system, it can be divided into two forms. The force distribution ratios for the three cables can be calculated based on the robot's actual motion state. In the global coordinate system, taking the cross product of the position vector of the application point and the force direction vector yields the corresponding 3D moment vector. Using the principles of mechanics, the moment balance equation for the robot in a static state is established to ensure that the driving moment Q_D and the resisting moment Q_R are balanced.

$$\begin{array}{}\text{(17)}& \begin{array}{rl}\mathit{\delta }{F}^{\mathrm{a}}& =k_{\mathrm{rod}}{\mathit{\delta }}_x{L}^{\mathrm{a}}\left[\begin{array}{c}\frac{{}^{\mathbf{R}}{\mathbf{T}}_{\mathbf{P}}^{-\mathrm{1}}\left({\mathit{\delta }}^{\mathrm{a}}{x}_{d\mathrm{1}}\right)}{\Vert {\mathit{\delta }}^{\mathrm{a}}{x}_{d\mathrm{1}}\Vert }\\ \frac{{}^{\mathbf{R}}{\mathbf{T}}_{\mathbf{P}}^{-\mathrm{1}}\left({\mathit{\delta }}^{\mathrm{a}}{x}_{d\mathrm{2}}\right)}{\Vert {\mathit{\delta }}^{\mathrm{a}}{x}_{d\mathrm{2}}\Vert }\end{array}\right]\\ & =\left[\begin{array}{c}{\left[\begin{array}{cccc}\mathrm{\dots }& \mathrm{\dots }& {\mathit{\delta }}^z{F}_{\mathrm{1}}^{\mathrm{a}}& \mathrm{\dots }\end{array}\right]}^{\mathrm{T}}\\ {\left[\begin{array}{cccc}\mathrm{\dots }& \mathrm{\dots }& {\mathit{\delta }}^z{F}_{\mathrm{2}}^{\mathrm{a}}& \mathrm{\dots }\end{array}\right]}^{\mathrm{T}}\end{array}\right]\\ \mathit{\delta }{F}^{\mathrm{b}}& =k_{\mathrm{rod}}{\mathit{\delta }}_x{L}^{\mathrm{b}}\left[\begin{array}{c}\frac{{}^{\mathbf{R}}{\mathbf{T}}_{\mathbf{P}}^{-\mathrm{1}}\left({\mathit{\delta }}^{\mathrm{b}}{x}_{d\mathrm{1}}\right)}{\Vert {\mathit{\delta }}^{\mathrm{b}}{x}_{d\mathrm{1}}\Vert }\\ \frac{{}^{\mathbf{R}}{\mathbf{T}}_{\mathbf{P}}^{-\mathrm{1}}\left({\mathit{\delta }}^{\mathrm{b}}{x}_{d\mathrm{2}}\right)}{\Vert {\mathit{\delta }}^{\mathrm{b}}{x}_{d\mathrm{2}}\Vert }\end{array}\right]\\ & =\left[\begin{array}{c}{\left[\begin{array}{cccc}\mathrm{\dots }& \mathrm{\dots }& {\mathit{\delta }}^z{F}_{\mathrm{1}}^{\mathrm{b}}& \mathrm{\dots }\end{array}\right]}^{\mathrm{T}}\\ {\left[\begin{array}{cccc}\mathrm{\dots }& \mathrm{\dots }& {\mathit{\delta }}^z{F}_{\mathrm{2}}^{\mathrm{b}}& \mathrm{\dots }\end{array}\right]}^{\mathrm{T}}\end{array}\right]\\ \mathit{\delta }{F}^{\mathrm{c}}& =k_{\mathrm{rod}}{\mathit{\delta }}_x{L}^{\mathrm{c}}\left[\begin{array}{c}\frac{{}^{\mathbf{R}}{\mathbf{T}}_{\mathbf{P}}^{-\mathrm{1}}\left({\mathit{\delta }}^{\mathrm{a}}{x}_{d\mathrm{1}}\right)}{\Vert {\mathit{\delta }}^{\mathrm{c}}{x}_{d\mathrm{1}}\Vert }\\ \frac{{}^{\mathbf{R}}{\mathbf{T}}_{\mathbf{P}}^{-\mathrm{1}}\left({\mathit{\delta }}^{\mathrm{a}}{x}_{d\mathrm{2}}\right)}{\Vert {\mathit{\delta }}^{\mathrm{c}}{x}_{d\mathrm{2}}\Vert }\end{array}\right]\\ & =\left[\begin{array}{c}{\left[\begin{array}{cccc}\mathrm{\dots }& \mathrm{\dots }& {\mathit{\delta }}^z{F}_{\mathrm{1}}^{\mathrm{c}}& \mathrm{\dots }\end{array}\right]}^{\mathrm{T}}\\ {\left[\begin{array}{cccc}\mathrm{\dots }& \mathrm{\dots }& {\mathit{\delta }}^z{F}_{\mathrm{2}}^{\mathrm{c}}& \mathrm{\dots }\end{array}\right]}^{\mathrm{T}}\end{array}\right]\end{array}\end{array}$$

$$\begin{array}{}\text{(18)}& \begin{array}{rl}\mathit{\delta }{F}_{\mathrm{all}}^B& =\left.\begin{array}{l}\left({\mathit{\delta }}^Y{B}_{\mathrm{\Delta }}^{\mathrm{all}}\right)\\ \mathrm{OR}\\ \left({\mathit{\delta }}^N{B}_{\mathrm{\Delta }}^{\mathrm{all}}\right)\end{array}\right\}{\cdot }^{\mathrm{out}}{k}_{\mathrm{spri}}\cdot \left[\begin{array}{c}\frac{{}^{\mathbf{R}}{\mathbf{T}}_{\mathbf{P}}^{-\mathrm{1}}\mathit{\delta }{B}_{\mathrm{\Delta }}^{\mathrm{a}}}{\Vert \mathit{\delta }{B}_{\mathrm{\Delta }}^{\mathrm{a}}\Vert }\\ \frac{{}^{\mathbf{R}}{\mathbf{T}}_{\mathbf{P}}^{-\mathrm{1}}\mathit{\delta }{B}_{\mathrm{\Delta }}^{\mathrm{b}}}{\Vert \mathit{\delta }{B}_{\mathrm{\Delta }}^{\mathrm{b}}\Vert }\\ \frac{{}^{\mathbf{R}}{\mathbf{T}}_{\mathbf{P}}^{-\mathrm{1}}\mathit{\delta }{B}_{\mathrm{\Delta }}^{\mathrm{c}}}{\Vert \mathit{\delta }{B}_{\mathrm{\Delta }}^{\mathrm{c}}\Vert }\end{array}\right]\\ & =\left[\begin{array}{c}{\left[\begin{array}{cccc}\mathrm{\dots }& \mathrm{\dots }& {\mathit{\delta }}^z{F}_{\mathrm{1}}^B& \mathrm{\dots }\end{array}\right]}^{\mathrm{T}}\\ {\left[\begin{array}{cccc}\mathrm{\dots }& \mathrm{\dots }& {\mathit{\delta }}^z{F}_{\mathrm{2}}^B& \mathrm{\dots }\end{array}\right]}^{\mathrm{T}}\\ {\left[\begin{array}{cccc}\mathrm{\dots }& \mathrm{\dots }& {\mathit{\delta }}^z{F}_{\mathrm{3}}^B& \mathrm{\dots }\end{array}\right]}^{\mathrm{T}}\end{array}\right]\end{array}\end{array}$$

Here, δF^a, δF^b, and δF^c are the force generated by a six-bar spring system during bending and $\mathit{\delta }{F}_{\mathrm{all}}^B$ is the force generated by the external flexible support system during bending.

By substituting Eqs. (14) and (15) into Eqs. (17) and (18), the deformation forces of the six-bar spring system and the external flexible support system during the bending motion can be determined. Here, k_rod represents the stiffness of the six-bar spring system, while the external flexible support system is divided into three parts, with its stiffness denoted as ^outk_spri. In the bent posture, the resisting moment Q_R, provided by both systems, balances the driving moment Q_D exerted by the driving cables. This balance establishes the total moment τ_all equilibrium equation.

$$\begin{array}{}\text{(19)}& \begin{array}{rl}& {\mathit{\tau }}_{\mathrm{all}}=Q_{\mathrm{D}}+Q_{\mathrm{R}}=\mathrm{0}\\ & \left\{\begin{array}{ll}{Q}_{\mathrm{D}}& =\sum _{u=\mathrm{1}}^{\mathrm{3}}\left[\begin{array}{c}\mathrm{0}\\ \mathrm{0}\\ {}^z{F_u}^{\mathrm{C}}\end{array}\right]\times C_u\\ Q_{\mathrm{R}}& =\sum _{u=\mathrm{1}}^{\mathrm{3}}\left[\begin{array}{l}{\begin{array}{c}\mathrm{0}\\ \end{array}}^z\\ {}^z{F}_u^B\end{array}\right]{\times }^k{B}_u\\ & +\left(\left[\begin{array}{c}\mathrm{0}\\ \mathrm{0}\\ {\mathit{\delta }}^z{\dot{F}}_{\mathrm{1}}^{\mathrm{a}}\end{array}\right]{\times }^{\mathrm{3}}{\stackrel{\mathrm{̃}}{D}}_{z\mathrm{2}}+\left[\begin{array}{c}\mathrm{0}\\ \mathrm{0}\\ {\mathit{\delta }}^z{\dot{F}}_{\mathrm{2}}^{\mathrm{a}}\end{array}\right]{\times }^{\mathrm{1}}{\stackrel{\mathrm{̃}}{D}}_{z\mathrm{1}}\right)\\ & +\left(\left[\begin{array}{c}\mathrm{0}\\ \mathrm{0}\\ {\mathit{\delta }}^z{\dot{F}}_{\mathrm{1}}^{\mathrm{b}}\end{array}\right]{\times }^{\mathrm{1}}{\stackrel{\mathrm{̃}}{D}}_{z\mathrm{2}}+\left[\begin{array}{c}\mathrm{0}\\ \mathrm{0}\\ {\mathit{\delta }}^z{\dot{F}}_{\mathrm{1}}^{\mathrm{b}}\end{array}\right]{\times }^{\mathrm{2}}{\stackrel{\mathrm{̃}}{D}}_{z\mathrm{1}}\right)\\ & +\left(\left[\begin{array}{c}\mathrm{0}\\ \mathrm{0}\\ {\mathit{\delta }}^z{\dot{F}}_{\mathrm{1}}^{\mathrm{c}}\end{array}\right]{\times }^{\mathrm{2}}{\stackrel{\mathrm{̃}}{D}}_{z\mathrm{2}}+\left[\begin{array}{c}\mathrm{0}\\ \mathrm{0}\\ {\mathit{\delta }}^z{\dot{F}}_{\mathrm{2}}^{\mathrm{c}}\end{array}\right]{\times }^{\mathrm{3}}{\stackrel{\mathrm{̃}}{D}}_{z\mathrm{1}}\right)\end{array}\right.\end{array}\end{array}$$

By solving Eqs. (16) and (19) together, the force conditions of the driving cables under different configuration spaces q_i can be derived. This analysis provides a basis for selecting materials for the parallel robot, making it more aligned with the motion characteristics of the human cervical spine. At the same time, the static analysis lays the theoretical foundation for dynamics research, providing preliminary validation of the robot's structure and mathematical model, while effectively simplifying the complexity of dynamic analysis. Furthermore, the joint torques and forces obtained from the static dynamics calculations can serve as important references for initial conditions or specific operating scenarios in dynamic analysis, thereby improving the accuracy and stability of dynamic analysis.

3.2 Dynamic analysis

Dynamic analysis is a key aspect of robot system design, with the primary goal of describing the mechanical changes the robot experiences during motion. Unlike static analysis, which focuses on the forces acting on the robot when it is stationary, dynamic analysis considers the dynamic responses caused by factors such as inertia effects, Coriolis force, and centrifugal force. These factors not only influence the robot's motion trajectory but also determine the torque Q_D required at the driving joints and the overall performance. Therefore, accurate dynamic modeling is crucial for effective robot control, path planning, and performance optimization.

Building on static analysis, the dynamic behavior of the parallel robot is further investigated. By analyzing the mass distribution of the platform, kinematic constraints, and coupling effects between the platforms, the mechanical behavior during the robot's motion is modeled. During this process, dynamic effects such as inertial force, Coriolis force, centrifugal force, and gravity are considered and incorporated into the dynamic modeling. Using Newton–Euler equations and the Lagrangian function, corresponding dynamic equations are constructed to accurately describe the robot's dynamic response during motion.

$$\begin{array}{}\text{(20)}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{\displaystyle \frac{\partial L_{\mathrm{s}}}{\partial \stackrel{\mathrm{‾}}{q_i}}}-{\displaystyle \frac{\partial L_{\mathrm{s}}}{\partial q_i}}=Q_{\mathrm{D}}\end{array}$$

Equation (20) represents the Euler–Lagrange equation, where L_s is the Lagrangian of the parallel robot and $\stackrel{\mathrm{‾}}{q_i}$ is the derivative of the configuration space q_i. Q_D is the required torque for the robot system, which is generated by the driving force of the cables.

Equation (21) represents the result of differentiating the configuration space q_i, while the Lagrangian in Eq. (20) can be expressed by Eq. (22). The Lagrangian is obtained from the difference between the system's kinetic energy and potential energy. By substituting Eqs. (21) and (22) into Eq. (20), the new dynamic equations for the parallel robot can be derived.

$$\begin{array}{}\text{(21)}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{\displaystyle \frac{\partial T_{\mathrm{s}}}{\partial {\stackrel{\mathrm{‾}}{q}}_i}}-{\displaystyle \frac{\partial T_{\mathrm{s}}}{\partial q_i}}+{\displaystyle \frac{\partial U_{\mathrm{s}}}{\partial q_i}}=Q_{\mathrm{D}}\end{array}$$

During the motion of the parallel robot, the kinetic energy T_s in the Lagrangian L_s can be decomposed into two parts: linear velocity kinetic energy ^qT_V and angular velocity kinetic energy ^qT_ω.

$$\begin{array}{}\text{(22)}& \begin{array}{rl}& T_{\mathrm{s}}{=}^q{T}_V{+}^q{T}_{\mathit{\omega }}\\ & {}^q{T}_V={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\cdot m_i{\cdot }^J{V}_q^{\mathrm{T}}{\cdot }^J{V}_q\\ & {}^q{T}_{\mathit{\omega }}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\cdot }^J{\mathit{\omega }}_q^{\mathrm{T}}{\cdot }^{\mathit{\omega }}{\mathbf{I}}_q{\cdot }^J{\mathit{\omega }}_q\end{array}\end{array}$$

Based on Eq. (22), the linear velocity kinetic energy ^qT_V and angular velocity kinetic energy ^qT_ω of the parallel robot system can be further derived. To calculate these two types of kinetic energy, the linear velocity V_q, angular velocity ω_q, and the system's inertia matrix ^ωI_q must first be determined. The calculation of linear velocity ^JV_q and angular velocity ^Jω_q depends on the corresponding Jacobian matrices ^qJ_V and ^qJ_ω in the configuration space q_i. Additionally, the system's inertia matrix ^ωI_q is determined by the mass m_i of the moving platform and the corresponding radius lengths.

$$\begin{array}{}\text{(23)}& {\displaystyle }\begin{array}{rl}{}^q{\mathbf{J}}_V& ={\displaystyle \frac{\partial P_q}{\partial q_i}}=\left[{\displaystyle \frac{\partial P_q}{\partial {\mathit{\varphi }}_i}}{\displaystyle \frac{\partial P_q}{\partial {\mathit{\theta }}_i}}\right]\\ {}^J{V}_q& =\left[\begin{array}{c}{V}_x\\ V_y\\ V_z\end{array}\right]{=}^q{\mathbf{J}}_V\cdot {\stackrel{\mathrm{‾}}{q}}_i\\ & =\left[\begin{array}{cc}\frac{\partial P_q}{\partial {\mathit{\varphi }}_i}& \frac{\partial P_q}{\partial {\mathit{\theta }}_i}\end{array}\right]\cdot {\left[\begin{array}{cc}{\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_i& {\stackrel{\mathrm{‾}}{\mathit{\theta }}}_i\end{array}\right]}^{\mathrm{T}}\end{array}\text{(24)}& {\displaystyle }\begin{array}{rl}& {\stackrel{\mathrm{‾}}{R}}_q\left({\mathit{\varphi }}_i,{\mathit{\theta }}_i\right)={\displaystyle \frac{\partial R_q\left({\mathit{\varphi }}_i,{\mathit{\theta }}_i\right)}{\partial {\mathit{\theta }}_i}}\cdot {\stackrel{\mathrm{‾}}{\mathit{\theta }}}_i+{\displaystyle \frac{\partial R_q\left({\mathit{\varphi }}_i,{\mathit{\theta }}_i\right)}{\partial {\mathit{\varphi }}_i}}\cdot {\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_i\\ & {}^J{\mathit{\omega }}_q=\left[\begin{array}{c}{\mathit{\omega }}_x\\ {\mathit{\omega }}_y\\ {\mathit{\omega }}_z\end{array}\right]{=}^q{\mathbf{J}}_{\mathit{\omega }}\cdot {\stackrel{\mathrm{‾}}{q}}_i\leftrightarrow {\stackrel{\mathrm{‾}}{R}}_q\left({\mathit{\varphi }}_i,{\mathit{\theta }}_i\right)\cdot R_q{\left({\mathit{\varphi }}_i,{\mathit{\theta }}_i\right)}^{-\mathrm{1}}\\ & =\left[\begin{array}{ccc}\mathrm{0}& -{\mathit{\omega }}_z& {\mathit{\omega }}_y\\ {\mathit{\omega }}_z& \mathrm{0}& -{\mathit{\omega }}_x\\ -{\mathit{\omega }}_y& {\mathit{\omega }}_x& \mathrm{0}\end{array}\right]\end{array}\end{array}$$

Using Eqs. (23) and (24), the linear velocity ^JV_q and angular velocity ^Jω_q in the configuration space q_i can be obtained and expressed as three-dimensional vectors.

$$\begin{array}{}\text{(25)}& \begin{array}{rl}{}^{\mathit{\omega }}{\mathbf{I}}_q& =\left[\begin{array}{cll}{}^x{I}_{\mathit{\omega }i}& & \\ & {}^y{I}_{\mathit{\omega }i}& \\ & & {}^z{I}_{\mathit{\omega }i}\end{array}\right]\\ & =\left[\begin{array}{ccc}\frac{\mathrm{1}}{\mathrm{4}}{m}_i{r}_{\mathrm{c}}& & \\ & \frac{\mathrm{1}}{\mathrm{4}}{m}_i{r}_{\mathrm{c}}& \\ & & \frac{\mathrm{1}}{\mathrm{2}}{m}_i{r}_{\mathrm{c}}\end{array}\right]\end{array}\end{array}$$

Equation (25) provides the inertia matrix of the moving platform, whose values can be determined based on the shape of the platform to calculate the moments of inertia in the XYZ directions. By substituting Eqs. (23), (24), and (25) into Eq. (22), the kinetic energy T_s of the parallel robot system can be obtained. For Eq. (21), the potential energy of the parallel robot system still needs to be determined, including both the gravitational potential energy of the robot and the elastic potential energy from the six-bar spring system and the external flexible support system.

$$\begin{array}{}\text{(26)}& U_{\mathrm{s}}\left(q_i\right)=U_{\mathrm{gray}}+U_{\mathrm{spri}}\end{array}$$

In analyzing the gravitational potential energy U_gray in Eq. (26), the gravitational potential energy variation of the six-bar spring system and the external flexible support mechanism is neglected in the dynamic formulation. This simplification is justified because these components possess relatively lightweight structures and undergo constrained symmetric micro-deformations, resulting in vertical center-of-mass displacements whose potential energy contribution is orders of magnitude smaller than the primary rotational dynamic work and cable tension strain energy. The gravitational force acting on the moving platform can be calculated using its mass and the actual gravitational acceleration. The gravitational potential energy U_gray changes in response to the real-time Z axis position of the moving platform in the global coordinate system. During this process, the moving platform is treated as a particle located at the center of mass P_q, allowing the gravitational potential energy changes during the robot's motion to be analyzed.

$$\begin{array}{}\text{(27)}& \left\{\begin{array}{c}{g}_m={\left[\begin{array}{ccc}\mathrm{0}& \mathrm{0}& g_{\mathrm{true}}\end{array}\right]}^{\mathrm{T}}\\ U_{\mathrm{gray}}=m_i{g}_m^{\mathrm{T}}{P}_q({\mathit{\varphi }}_i,{\mathit{\theta }}_i)\end{array}\right.\end{array}$$

The elastic potential energy U_spri is provided by the six-bar spring system and the external flexible support system. During the parallel robot's bending motion, these systems generate elastic potential energy due to deformation, and the direction of deformation aligns with the axis. Equations (14) and (15) describe the deformation generated during the bending motion. Based on Hooke's law in elasticity, and using the elastic stiffness k_rod of the six-bar spring system and the stiffness ^outk_spri of each spring in the external flexible support system, the variable radius δr_z is used to determine r_Bz, which is then used to analyze the elastic potential energy U_spri of the parallel robot system.

$$\begin{array}{}\text{(28)}& {}^{\mathrm{sp}}\mathit{\delta }{L}_{\mathrm{all}}=\sum _{k=\mathrm{1}}^{\mathrm{2}}\left[{\mathit{\delta }}_x{L}^{\mathrm{a}}(k{)}^{\mathrm{2}}+{\mathit{\delta }}_x{L}^{\mathrm{b}}(k{)}^{\mathrm{2}}+{\mathit{\delta }}_x{L}^{\mathrm{c}}(k{)}^{\mathrm{2}}\right]\text{(29)}& {\displaystyle }{U}_{\mathrm{spri}}=\left\{\begin{array}{l}{}_{\mathrm{1}}^s{U}_{\mathrm{spri}}\left(\mathit{\delta }{r}_z\left(\mathit{\beta }\right)=\mathrm{0}\right)=\\ \frac{\mathrm{1}}{\mathrm{2}}\left[{}^{\mathrm{sp}}\mathit{\delta }{L}_{\mathrm{all}}{k}_{\mathrm{rod}}{+}^{\mathrm{out}}{k}_{\mathrm{spri}}\sum _{k=\mathrm{1}}^{\mathrm{3}}{\left(\mathit{\delta }{B}_{\mathrm{\Delta }}^{\mathrm{all}}\left(k\right)\right)}^{\mathrm{2}}\right]\\ {}_{\mathrm{2}}^s{U}_{\mathrm{spri}}\left(\mathit{\delta }{r}_z\left(\mathit{\beta }\right)\ne \mathrm{0}\right)=\\ \frac{\mathrm{1}}{\mathrm{2}}\left[{}^{\mathrm{sp}}\mathit{\delta }{L}_{\mathrm{all}}{k}_{\mathrm{rod}}{+}^{\mathrm{out}}{k}_{\mathrm{spri}}\sum _k^{\mathrm{3}}{\left({\mathit{\delta }}^N{B}_{\mathrm{\Delta }}^{\mathrm{all}}\left(k\right)\right)}^{\mathrm{2}}\right]\end{array}\right.\end{array}$$

When analyzing elastic potential energy, the square of the deformation is calculated using Hooke's law and the work principle in Eq. (28), providing a basis for further calculation in Eq. (29). By considering the variable radius δr_z, an analysis model of the elastic potential energy for the six-bar spring system and the external flexible support system is established using Eq. (29). By combining the gravitational potential energy of the parallel robot from Eq. (27) with the elastic potential energy from Eq. (29), the overall potential energy expression for the parallel robot during motion can be determined.

$$\begin{array}{}\text{(30)}& \begin{array}{rl}{U}_{\mathrm{s}}^{\mathrm{all}}\left(q_i\right)=& U_{\mathrm{gray}}+U_{\mathrm{spri}}\\ =& m_i{g}_m^{\mathrm{T}}{P}_q({\mathit{\varphi }}_i,{\mathit{\theta }}_i)E_{\mathrm{2}\times \mathrm{2}}+{{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}}^{\mathrm{sp}}\mathit{\delta }{L}_{\mathrm{all}}{k}_{\mathrm{rod}}{E}_{\mathrm{2}\times \mathrm{2}}\\ & {+}^{\mathrm{out}}{k}_{\mathrm{spri}}\left[\begin{array}{c}\sum _{k=\mathrm{1}}^{\mathrm{3}}{\left({\mathit{\delta }}^Y{B}_{\mathrm{\Delta }}^{\mathrm{all}}\left(k\right)\right)}^{\mathrm{2}}\\ \sum _k^{\mathrm{3}}{\left({\mathit{\delta }}^N{B}_{\mathrm{\Delta }}^{\mathrm{all}}\left(k\right)\right)}^{\mathrm{2}}\end{array}\right]E_{\mathrm{2}\times \mathrm{2}}\end{array}\end{array}$$

Furthermore, by substituting Eqs. (22) and (30) into Eq. (21), the dynamic equations of the parallel robot can be derived. Since Eqs. (22) and (30) are constructed based on the configuration space q_i, partial derivatives of q_i can be computed to more accurately describe the robot's dynamic characteristics. This dynamic equation not only helps optimize control strategies to improve trajectory tracking accuracy and motion stability but also predicts key parameters such as system inertia and torque distribution, providing theoretical support for practical applications.

4 Simulation analysis of the parallel robot

Based on the principles of bionics and robotics, this paper designs a novel cable-driven, variable-stiffness, compliant parallel robot by referring to the structural characteristics of human cervical bones and muscles. Through the study of the functional relationship between the bones and muscles during cervical flexion, the structure of traditional parallel robots is innovatively optimized. The design adopts a compliant structure that simulates the interaction between bones and muscles, incorporating a unique variable stiffness mechanism to emulate the working characteristics of human neck muscles. Building on the novel configuration and functional innovation, this paper establishes the robot's kinematic model. It derives the G function and its inverse function A, which map from the configuration space q_i to the actuation space ^Δl. Additionally, static and dynamic analyses are performed, leading to the complete dynamic equations of the robot.

Table 1 Basic parameters of the parallel robot.

To validate the posture control and dynamic characteristics of the cable-driven, variable-stiffness, compliant parallel robot, kinematic and dynamic simulations are conducted. As shown in Table 1, the basic parameters required for the simulation are provided, including the radius of each mounting hole and the corresponding angle values. Figure 6 visually presents the robot's form in different spatial poses, and the corresponding driving cable lengths are calculated to verify the accuracy and reliability of the established kinematic model.

Figure 6 Spatial posture and cable length values of the parallel robot.

According to the G function defined in Eq. (12), the length of the driving cables in the parallel robot can be calculated from the configuration space q_i, which determines the specific posture of the parallel robot shown in Fig. 6. Based on the normal human neck flexion range, and considering the reduction of neck movement by about 30 % when wearing heavy wearable protective equipment, the maximum flexion angle for the parallel robot is set to 30°, closely mimicking the true motion characteristics of the human neck (Zhang et al., 2023).

Using the kinematic equations of the parallel robot and the range of flexion and azimuth angles, the robot's task workspace can be further analyzed. This analysis helps define the reachable range of the robot's end effector, ensuring that the design of the cable-driven, stiffness-variable, compliant parallel robot effectively simulates human neck movement. During the task space analysis, as shown in Fig. 6, the robot's motion shape in 3D space is displayed in lines, and simulation results verify the robot's spatial performance during actual task execution, ensuring that the design meets the expected goals.

Figure 7 Workspace analysis of the parallel robot.

Figure 7 illustrates the workspace of the parallel robot. When the actuation space ^Δl is in the initial state, the end effector P_q is positioned at the highest point during the robot's motion. By analyzing the movement of the platform, the motion state of the driving cables can be observed. From a top-down view, the robot's movement range in the XY plane is clearly visible, and the driving cable's motion trajectory always stays within the platform's motion range. The workspace analysis provides data support for evaluating neck flexibility and the range of motion.

Through kinematic simulations of the parallel robot, its motion posture is determined, providing the foundation for further analysis of its static and dynamic characteristics. For each posture, the tension $F_d^{\mathrm{C}}$ (d=1, 2, 3) in the driving cables and the variation in cable length relative to the initial state (i.e., the change in actuation space ^Δl) are calculated. In the robot's initial state, pre-tensioning the driving cables not only stabilizes its posture but also ensures that the initial length between the platforms remains l_o. By analyzing the variations in actuation space ^Δl and the corresponding cable tensions, the robot can achieve motion control through force-position control, effectively simulating the natural motion of the human neck.

Figure 8 Variation in cable tension and actuation space of the parallel robot.

As shown in Fig. 8, by analyzing the range of values for the bending angle θ_i and azimuth angle ϕ_i, along with the dynamic equations and the G function, the force variation surface of the driving cables and the actuation space ^Δl variation surface are derived. The results show that the cable tension surface and the actuation space ^Δl variation curve exhibit similar trigonometric characteristics. On the left and right sides of Fig. 8, the tension of the same driving cable and the change in actuation space ^Δl are shown. Comparative analysis indicates that this relationship helps to better understand the parallel robot's dynamic behavior and improves the precision of motion trajectory control. As the actuation space variation of cable u increases, its tension also increases. Additionally, as the bending angle θ_i of the parallel robot increases, both the cable tension and the actuation space variation gradually increase, though the overall trend still follows a trigonometric curve. When the external flexible support system reaches the limit of its extension and contraction deformation, the tension in the flexible cable reaches its minimum and maximum values.

The parallel robot is connected through a six-bar spring system between the moving platform and the fixed platform, along with an external flexible support system, forming a rigid–flexible coupling structure. This effectively simulates the motion characteristics of the human neck. In the system, the moving platform is rigidly connected to the six-bar spring system and the flexible support system. At the boundary where the moving platform and flexible support system are coupled, force balance is achieved, and the end displacements of both the rigid and flexible bodies remain continuous. Based on this, the dynamic equations of the rigid and flexible bodies are coupled, leading to the establishment of the parallel robot's rigid–flexible coupling equation.

$$\begin{array}{}\text{(31)}& \left\{\begin{array}{l}{\mathit{\sigma }}_{\mathit{\omega }}=\frac{\mathrm{d}}{\mathrm{d}t}{\Vert }^J{\mathit{\omega }}_q\Vert \\ {\mathit{\sigma }}_V=\frac{\mathrm{d}}{\mathrm{d}t}{\Vert }^J{V}_q\Vert \end{array}\right.\end{array}$$

By substituting Eqs. (23) and (24) into Eq. (31), the angular acceleration and linear acceleration during the robot's motion can be calculated. During this process, the torque Q_D applied by the driving cables on the moving platform balances with the torque generated by the deformation of the flexible support system. Based on Eq. (31), the rigid–flexible coupling equation for the parallel robot is established.

$$\begin{array}{}\text{(32)}& \begin{array}{rl}& {\mathbf{I}}_B={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{m}_i{r}_{\mathrm{c}}\\ & {\mathbf{I}}_B{\mathit{\sigma }}_{\mathit{\omega }}={\mathit{\tau }}_{\mathrm{in}}=Q_{\mathrm{R}}=-Q_{\mathrm{D}}\end{array}\end{array}$$

Equation (32) describes the rigid–flexible coupling equation of the parallel robot, where I_B is the inertia matrix. The inertia coupling matrix can be calculated from the angular acceleration of the robot during motion. The resistance torque generated by the flexible system forms the torque coupling matrix τ_in, determined by Eq. (19), which balances with the driving torque Q_D applied by the driving cables. This leads to the rigid–flexible coupling equation of the parallel robot. The rigid–flexible coupling analysis curve is shown in Fig. 9.

Figure 9 Rigid–flexible coupling analysis of the parallel robot.

For the parallel robot mimicking the human neck muscles with a variable stiffness system, the stiffness is determined by the stiffness of the driving cables ΔK_cable and the stiffness of the flexible system K_robot. By rotating the structure to adjust the motion position of the slide rail, the radius of the external flexible support system can be changed, thereby adjusting the stiffness of the driving cables. When the parallel robot is in a fixed posture, the length of the driving cables cannot change, but adjusting the radius of the flexible support system causes a variation in the resistance torque Q_R of the flexible system. To maintain the posture of the robot, the driving torque Q_D of the cables is adjusted accordingly, achieving the function of mimicking the variable stiffness of human neck muscles.

$$\begin{array}{}\text{(33)}& {\displaystyle }\left\{\begin{array}{l}{K}_{\mathrm{robot}}=\mathrm{6}{k}_{\mathrm{rod}}+{\mathrm{3}}^{\mathrm{out}}{k}_{\mathrm{spri}}\\ \mathrm{\Delta }{K}_{\mathrm{cable}}=\sum _{d=\mathrm{1}}^{\mathrm{3}}\frac{F_d^{\mathrm{C}}}{\Vert \left({}^{\mathbf{R}}{\mathbf{T}}_{\mathbf{P}}-E\right){\stackrel{\mathrm{̃}}{C}}_d\Vert }\end{array}\right.\text{(34)}& {\displaystyle }{K}_{\mathrm{all}}=K_{\mathrm{robot}}+\mathrm{\Delta }{K}_{\mathrm{cable}}\end{array}$$

In Eq. (33), the stiffness of the driving cables ΔK_cable and the stiffness of the flexible system K_robot are key components of the overall stiffness K_all of the parallel robot. By substituting Eq. (33) into Eq. (34), the stiffness variation rule of the parallel robot can be obtained. To quantitatively evaluate the static and dynamic performance from the numerical results, specific simulation evaluation metrics are defined in Eqs. (30)–(34).

Figure 10 Stiffness variation of the variable stiffness system.

Figure 10 illustrates the variation in the stiffness of the parallel robot. By adjusting the rotational structure to modify the radius of the external flexible support system, the robot's stiffness can be controlled. The figure displays four surface variations corresponding to four different values of the variable radius δr_z(β) of the external flexible support system. As the angle β increases, the variable radius δr_z(β) also increases, resulting in a corresponding rise in the overall stiffness K_all. The value of K_all reflects the robot's ability to resist external disturbances, analogous to the stabilizing function of human neck muscles when subjected to external forces. The simulation analysis clearly demonstrates the system's ability to adjust stiffness at different postures, further validating its effectiveness and feasibility in mimicking the variable stiffness behavior of human neck muscles. This approach not only provides valuable design insights and a theoretical foundation for the development of biomimetic robots that replicate human muscle functions but also broadens the potential applications of parallel robots in complex and dynamic environments.

5 Discussion

Through the aforementioned simulation study, the biomimetic human cervical parallel robot proposed in this paper outperforms the cable-driven parallel robot with a spring-based spine designed by Gao and Song, particularly in simulating the actual range of motion of the human cervical spine (Grazioso et al., 2019). By expanding the range of the bending angle θ_i, the parallel robot increases its task space, as shown in the blue area of Fig. 11. Additionally, the robot more accurately replicates the characteristics of the cervical bones and muscles in its biomimetic structure. Detailed kinematic and dynamic analyses validate the feasibility of the proposed design and highlight its performance advantages. The kinematic simulation results show that the parallel robot can achieve high motion control precision, effectively adjusting the cervical posture. Furthermore, dynamic analysis demonstrates that the robot's stiffness adjustment mechanism can flexibly alter its overall stiffness, significantly enhancing the structural stability and resistance to external disturbances.

Figure 11 Workspace distribution and optimized region of the proposed parallel robot: (a) Y–Z plane view, (b) X–Z plane view showing the optimization area, (c) top view in X–Y plane, and (d) 3D spatial workspace plot.

Nevertheless, from the perspective of human cervical muscle movement mechanisms, the parallel robot still has room for improvement. The current variable stiffness system relies on adjusting the radius of the external flexible support system to achieve stiffness changes. However, this design may cause variations in the robot's size, affecting its biomimetic fidelity and practical application performance. Therefore, future work will focus on optimizing the structural design to better match the shape of the human cervical spine, while maintaining its functional capabilities, in order to enhance both its biomimicry and practicality.

Additionally, the friction effects were neglected in the dynamic calculations of the system in this paper. However, in practical applications, friction can have a non-negligible impact on the robot's motion performance and the accuracy of the end effector. By using lubricants, reducing the contact area along the drive path, and incorporating self-lubricating bushings at the cable fixation points, the negative impact of friction can be effectively reduced, thereby improving the overall performance of the robot.

In summary, the biomimetic human cervical parallel robot still requires further optimization in its structural design to enhance performance, stability, and biomimetic characteristics. Additionally, the frictional effects should be incorporated into the system analysis to account for their potential impact on the robot's motion. Furthermore, differences in motion time during kinematic simulations should be considered, as they may be limited by the computational capacity of the devices used. Future work should focus on optimizing and adjusting the design and simulation strategies to meet the specific needs of various application scenarios.

The robot system designed in this study has demonstrated clear advantages in structural design and stiffness adjustment, but there are still several limitations. Firstly, the current stiffness adjustment relies on the radius change of the external flexible support system, and the response speed and adjustment accuracy need to be further verified through prototype experiments. Secondly, the overall dynamic performance and control strategy of the robot have not been fully developed and improved, and complex working conditions in practical applications place high demands on system stability. In addition, this study focuses on simulating the motion range and stiffness characteristics of the bionic cervical spine and has not addressed the functional adaptability and safety assessment of the system in specific application scenarios such as humanoid robots and rehabilitation. There are still many challenges in the development of practical applications in humanoid robotics and rehabilitation training systems, such as real-time response and robustness of control systems and multi-sensor information fusion.

6 Conclusions

This paper presents an innovative cable-driven parallel robot with a variable stiffness mechanism, inspired by the motion mechanism of human cervical spine. It improves upon traditional parallel robots by incorporating rigid–flexible coupling and variable stiffness functions, allowing the robot to simulate the collaborative interaction between muscles and bones during neck motion. The robot consists of a motion platform, a fixed platform, drive cables, a six-bar spring system, and an external flexible support system. It is driven by three cables and successfully mimics the motion of the human cervical spine.

The robot adopts an innovative design combining rigid–flexible coupling with a variable stiffness system. This design provides necessary compliance during bending and allows for adjustable stiffness. It effectively simulates the resistance of neck muscles to external disturbances. This novel design offers a new approach for the structural concept of biomimetic cervical robots. Based on this, a systematic analysis of the general kinematics of the parallel robot was conducted. The G and A functions between the configuration space q_i and the actuation space ^Δl are established. The kinematic equations are optimized, leading to a significant improvement in the efficiency of kinematic calculations. In-depth static and dynamic analyses are performed. Using force-position control methods, the precise pose of the robot in space is determined. This greatly enhances the motion accuracy and trajectory stability of the end effector. The rigid–flexible coupling structure and variable stiffness system are also thoroughly studied. Adjusting the radius of the external flexible support system effectively controls the robot's stiffness. This improves its ability to resist external disturbances. These findings provide theoretical support for the rigid–flexible coupling and variable stiffness system. Finally, simulations confirm the effectiveness of the variable stiffness system within its adjustable range. This innovative design opens new research directions for developing biomimetic cervical parallel robots.

Despite these promising results, certain limitations of the present study should be acknowledged. First, the current model omits non-linear cable friction and structural hysteresis, while simplified assumptions were made regarding the mass of specific components. Second, the coupling relationship between bones and muscles has not yet been fully optimized based on detailed cervical spine anatomy. To address these limitations, future research will focus on the systematic optimization of the robotic architecture. Specifically, the dynamic mass and volume distributions of the human head will be fully integrated into the end-effector load profile. Furthermore, the bone–muscle coupling mechanism will be refined using precise anatomical data of the cervical spine to emulate human neck kinematics with higher fidelity. Additionally, practical environmental non-linearities, such as friction, will be incorporated into the governing equations, and force/pose sensors will be integrated into the end effector to establish a closed-loop feedback control scheme, thereby enhancing the operational accuracy and adaptability of the platform in physical applications. Overall, this work provides a novel structural paradigm for biomimetic human cervical parallel robots and offers fresh insights into biological muscle–bone interaction, which is expected to broaden the potential application scenarios in the field of advanced biomimetic robotics.