Kinematic and Dynamic Modeling of 3DOF Variable Stiffness Links Manipulator with Experimental Validation

Abstract

Variable stiffness link (VSL) manipulators are robotic arms that can adjust their link stiffness in real time to improve their adaptability and precision. They are particularly useful in industrial environments where safe collaboration with human workers is required. However, modeling and controlling these non-linear systems is a major challenge due to their complexity. This research paper presents a mathematical model for a 3DOF VSL manipulator, which is the first step towards optimizing performance, improving safety, and reducing costs. The accuracy and reliability of the model are demonstrated through verification experiments that strengthen confidence in its validity for engineering and scientific research. This study contributes to the understanding of the dynamics of VSL manipulators and provides insights for future advances in the use of such robots. By using the proposed model, the efficiency and precision of VSL manipulators can be improved while ensuring safe human–robot interaction in various industrial applications.

1. Introduction

Soft robotics focuses on creating robots capable of working together safely with humans by adapting to their surroundings. The use of flexible link manipulators could bring several advantages: higher payload capacity, lower energy consumption, a more cost-effective design, faster movements, greater reach, and safer operation [1]. Variable stiffness (VS) is essential for this purpose, and numerous techniques have been proposed to achieve this, including soft actuators, shape-memory materials, additive manufacturing, and VS joints [2]. Researchers have developed various soft robotic systems and variable stiffness grippers, including a 3D-printed monolithic soft gripper and a fluidic actuation system [3,4,5,6,7,8,9,10,11]. Composite materials such as shape-memory alloys and magnetorheological elastomers (MREs) have also been used to modulate stiffness [12]. The modeling of flexible joints is becoming increasingly important in the development and control of robotic manipulators, as they need to move with high precision and speed. Various methods have been introduced for the manipulators’ dynamic modeling with flexible links.

The decoupling method, based on according to the decoupled natural orthogonal complement matrices and Euler–Lagrange equations, has been shown to provide computational efficiency and stability and has been extended to various systems, including serial, parallel, closed-loop, flexible, and hyper-degree-of-freedom systems [13,14]. The equivalent rigid-link system approach is proposed to model spatial lightweight industrial robots, obtaining satisfactory alignment between the numerical outputs and the measured data points [15]. The formulation based on an Equivalent Rigid-Link System has been compared to other techniques for reducing flexible multibody dynamics model order [16]. An ERLS-based 3D dynamics formulation has also been developed for flexible articulated robots that can handle large displacements and small elastic deformations [17]. The Euler–Bernoulli method has been employed to develop single-link manipulators with cables dynamic models and payloads based on the extended Hamilton principle and validated by simulations and experiments [18,19,20]. The assumed mode method has also been applied to control flexible manipulators, investigating NN controllers with output feedback and full state and trajectory tracking fuzzy NN controllers and vibration suppression, with simulations and experiments showing their feasibility [21,22,23]. Other researches have examined the consequences of the elastic deflection boundary conditions and a flexible manipulator’s torques [24]. A flexible manipulator model is formulated through Timoshenko beam theory. Divergences are quantified, and robust controllers are synthesized with l-synthesis, ensuring effective tip tracking and robustness to uncertainties and noise [25].

The popular technique finite element method (FEM) is a technique to model hyper-redundant or serpentine robots [26,27,28]. The Hamilton method is used for modeling single-link flexible beams with a tip mass [29]. The lumped parameter method is an accurate and efficient technique for spatially compliant and multi-link manipulators [30,31,32,33]. The Newton–Euler method is used for modeling highly flexible 3D manipulators featuring links of arbitrary shape [34]. Finally, the pseudo-rigid body method is used for modeling compliant mechanisms, where the generalized PRBM PPRR model has better accuracy compared to other PRBM models [35,36]. The Timoshenko beam theory method is used to investigate the frequency sensitivity and orthogonality of the modes of a flexible beam [37]. The sensitivity of flexible manipulators' vibration frequencies is analyzed to system parameters using variational methods, discussing both Euler-Bernoulli and Timoshenko models, and introduces a novel orthogonal relations derivation method [38].

On the other hand, the Lagrangian technique is applied for developing flexible manipulators’ dynamic models, including payload and hub inertia, using finite elements and assumed mode methods [39,40]. Advanced control strategies, such as fuzzy set theory and passivity-based controllers, have also been proposed to achieve precise control and energy dissipation in these systems.

In summary, variable stiffness manipulators have been designed and developed using various approaches and materials that enable soft robots to perform different applications and tasks. These manipulators have the potential to revolutionize several industries, including manufacturing and healthcare. These modeling methods provide insights into the flexible link manipulators’ dynamics and can aid in the development and control of robotic systems. This manipulator’s dynamic model helps to predict and control the robot’s motion and the forces acting on it.

An innovative idea for the development of a collaborative 3DOF robot is presented by Stilli et al. using VSL. Using a hybrid of silicone rubber construction, polypropylene plastic mesh, and fabric materials, they have created links whose stiffness is controllable, and which are pneumatically actuated. These links can be adjusted to various stiffness levels by changing the pressure value within their structure, allowing the robot’s performance to be tailored to different applications [11,41].

The pressure sensor is used to sense collisions between the robot’s body and a human. When a collision is detected, the pressure in the links is released and the links become soft, so it does not cause any harm to the human worker. The performance of the VSL is evaluated when it is integrated into a robotic manipulator. They analyzed the different loads’ effects and pressures on the manipulator’s workspace and evaluated the effectiveness detection of the collision control system and hardware [11,41].

But the mathematical model is not derived for this 3D VSL manipulator. The lack of a derived mathematical model is a significant obstacle to the implementation of model-based control systems, a fundamental requirement for robust control methods.

Dynamic mathematical models play a significant role in the manipulators with flexible links development It enables accurate prediction of the complex behavior of the system, including link deflection, vibrations, and coupling effects. Accurate models form the foundation for designing effective control strategies and optimizing performance. They facilitate parameter estimation, system identification, and the development of vibration suppression techniques. Therefore, it is essential to address this critical gap in our current understanding to unlock the full capabilities of model-based control strategies and improve overall system performance.

This research makes an important contribution to the modeling of VSL robotic manipulators presented in [11,41]. In particular, a 3DOF VSL prototype comprehensive kinematic and dynamic models are derived. The dynamic model exclusively uses the Euler–Lagrange formulation with lumped parameters for representing the system efficiently yet accurately. The drawbacks of the other modeling methods and the advantages of this method are discussed in the Dynamics section. The complete experimental validation proves the proposed techniques effectiveness in capturing system behavior. By providing a precise mathematical representation of VSL dynamics, this work lays the foundation for optimizing performance, stability, and human–robot collaboration through model-based control strategies tailored to these VSL manipulators. As VSL robots are increasingly used in safety-critical tasks, the rigorous modeling methodology presented will support efforts to analyze, predict, and regulate the behavior of similar systems.

This paper contains Seven sections covering various aspects of the research. Section 2 deals with the mathematical model derivation. Section 3 deals with the simulation of the derived mathematical model. For Section 4, it discusses the derived model experimental validation. Section 5 contains the conclusion of the study. Appendix A contains the detailed matrices of the dynamic model. The last section contains a list of references used in this paper.

2. Mathematical Modelling

VSL manipulators offer a high degree of versatility and user-friendly operation. However, their non-linear dynamics, due to the flexibility of the links, pose a major challenge for modeling and control. Existing techniques, such as the assumed mode method, cannot effectively capture complex variable stiffness effects. Few studies analyze the VSL dynamically and do not validate its performance, despite the complex interactions between the links. Few works have derived a single-link manipulator’s dynamics or planner two-link manipulator. Therefore, precise, computationally efficient, and validated mathematical dynamics models for the VSL manipulator with 3DOF are essential to optimizing speed, precision, and stability in new VSL robots through model-based control.

Comprehensive forward kinematic and inverse kinematic models are derived following the Denavit–Hartenberg conventions. The dynamic model of a 3DOF VSL manipulator, using Euler–Lagrange and lumped parameters, captures kinetic and potential energy with discrete mass–spring–damper elements, ensuring accuracy. It takes into account important factors such as link masses, inertia, joint stiffness, damping, actuator torques, and gravity. The model efficiently represents the non-linear relationships between the end effector posture and the joint variables. Rigid-body links with torsional elements characterize variable stiffness links.

2.1. Kinematics

Kinematics is of central importance in the field of variable stiffness link (VSL) manipulators, as it provides essential insights into their movement and geometry. VSL manipulators offer a high degree of flexibility and compliance with their adaptable stiffness characteristics, making them suitable for tasks in unpredictable environments. Kinematic analysis focuses on understanding the correlation between joint angles or positions and the resultant the end effector’s orientation and position. This knowledge is essential for precise control and influences the manipulator’s ability to handle objects with finesse. Kinematic models enable the prediction and control of manipulator movements without considering complex forces and are therefore invaluable for various applications, from medical procedures to industrial automation [42].

2.1.1. Forward Kinematics

Forward kinematics are critical for variable stiffness link (VSL) manipulators to enable the end effector’s precise positioning and orientation control. Since VSLs are designed to adapt to changing environments, accurate prediction of effector positions based on joint variables is essential for controlled object interactions. Forward kinematics enables the mapping of joint angles/positions to Cartesian end-effector coordinates through kinematic equations. This enables the planning and execution of precise, adaptable movements tailored to the task at hand. In applications such as medical robotics, which require delicate interventions, forward models ensure safety and effectiveness. This VSL manipulator uses pressure-mediated links with variable stiffness, with two-bar actuators providing rigidity while enabling pressure feedback for collision detection. Consequently, the flexibility of the links results in deviations from the end effector during motion. The Denavit–Hartenberg convention is used to develop joint coordinate frames and homogeneous transformations from which trigonometric relations predict the end-effector pose of the manipulator for given joint angles. The experiments presented later support the assumptions of the kinematic model in predicting motion geometry. The Sketch for the robotic arm, along with the arrangement of each joint’s frame and lumped parameters is shown Figure 1.

Table 1. DH parameters to derive forward kinematics using Denavit–Hartenberg method.

i	${\mathit{\theta }}_{\mathit{i}}$	${\mathit{\alpha }}_{\mathit{i}}$	${\mathit{a}}_{\mathit{i}}$	${\mathit{d}}_{\mathit{i}}$
1	$q_1$	90	0	l1
2	$q_2$	0	l2	0
3	$q_3$	0	l3	0

presents the Denavit-Hartenberg (DH) parameters to derive forward kinematics using Denavit–Hartenberg method.

Equation (1) determines the manipulator forward kinematics. This equation represents the transformation matrix between frame O₀ and frame O₃.

$${}^{0}T_3={}^{0}T_1{}^{1}T_2{}^{3}T_2$$

(1)

Therefore, ${}^{0}T_1, {}^{1}T_2 \mathrm{a}\mathrm{n}\mathrm{d} {}^{3}T_2$ can be obtained by Equations (2)–(4), as follows:

$${}^{0}T_1=\left[\begin{array}{cccc}cos\left(q_1\right)& 0& sin\left(q_1\right)& 0\\ sin\left(q_1\right)& 0& -cos\left(q_1\right)& 0\\ 0& 0& 0& 0.04\\ 0& 0& 0& 1\end{array}\right]$$

(2)

$${}^{1}T_2=\left[\begin{array}{cccc}cos\left(q_2\right)& -sin\left(q_2\right)& 0& 0.29cos\left(q_2\right)\\ sin\left(q_2\right)& cos\left(q_2\right)& 0& 0.29sin\left(q_2\right)\\ 0& 0& 0& 0\\ 0& 0& 0& 1\end{array}\right]$$

(3)

$${}^{2}T_3=\left[\begin{array}{cccc}cos\left(q_3\right)& -sin\left(q_3\right)& 0& 0.24cos\left(q_3\right)\\ sin\left(q_3\right)& cos\left(q_3\right)& 0& 0.24sin\left(q_3\right)\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]$$

(4)

while (5) shows the transformation matrix between O₀ and O₃, where $n_x, n_y,$ $n_z,$ $o_x,$ $o_y,$ $o_z,$ $a_x,$ $a_y$, and $a_z$ are for the orientation and $p_x, p_y, \mathrm{a}\mathrm{n}\mathrm{d} p_z$ are for the x-axis, y-axis, and z-axis position, respectively.

$${}^{0}T_3=\left[\begin{array}{cccc}{n}_x& o_x& a_x& p_x\\ n_y& o_y& a_y& p_y\\ n_z& o_z& a_z& p_z\\ 0& 0& 0& 0\end{array}\right]$$

(5)

Thus, the end effector position can be obtained through extracting the transformation matrix fourth column from base to the manipulator end effector.

$$p_x=\cos \left(q_1\right)\left(l_2\cos \right(q_2)+l_3\cos (q_2+q_3\left)\right)$$

(6)

$$p_y=\sin \left(q_1\right)\left(l_2\cos \right(q_2)+l_3\cos (q_2+q_3\left)\right)$$

(7)

$$p_z=l_1+l_2\sin \left(q_2\right)+l_3\sin (q_2+q_3))$$

(8)

2.1.2. Inverse Kinematics

Inverse kinematics is a crucial element for the functionality of manipulators with variable stiffness links. These manipulators are designed to be versatile and adaptable, and inverse kinematics plays a central role in achieving this adaptability. Variable stiffness link manipulators must respond to external forces and changing demands and dynamically adjust their joint variables to precisely position the end effector. BY applying inverse kinematics, the manipulator can calculate the joint variables needed to accomplish the desired end-effector position. This capability is important for tasks that need controlled interaction with the environment, such as grasping objects with varying levels of safety or avoiding obstacles. In applications such as advanced manufacturing and human–robot collaboration, where the manipulation environment is complex and dynamic, inverse kinematics in variable stiffness link manipulators ensures not only precise positioning but also increased safety and adaptability [43]. Having established that the VSL manipulator can be treated as a rigid link manipulator using forward kinematics, inverse kinematics can also be approached using the same framework for rigid link manipulators. Equations (9)–(12) define the inverse kinematics models for a rigid anthropomorphic manipulator with three rotational degrees of freedom, which have been solved and used extensively in previous studies [44,45,46].

$$r=\sqrt{x^2+y^2+z^2}$$

(9)

$${\theta }_1={\tan }^{-1}\left(\frac{y}{x}\right)$$

(10)

$${\theta }_2=-{\mathit{cos}}^{-1}\left(\frac{x^2+y^2+z^2-l_2^2 l_3^2}{2l_2{l}_3}\right)$$

(11)

$${\theta }_3={\mathit{sin}}^{-1}\left(\frac{z}{r}\right)+{\tan }^{-1}\left(\frac{l_3\sin \left({\theta }_2\right)}{l_2+l_3\cos \left({\theta }_2\right)}\right)$$

(12)

2.2. Dynamics

Manipulator dynamics is the study of the equations governing motion, in particular how the manipulator responds to torques exerted by its actuators or external forces. There are different methods for solving the dynamics of flexible manipulators, as discussed in the introduction.

The Euler–Lagrange method with lumped parameters is a good choice for modeling a manipulator with flexible links for several reasons. First, the Euler–Lagrange method is known for its simplicity and accuracy in describing the complex mechanical systems dynamics. It allows the lumped parameters inclusion, which simplifies the representation of distributed parameters in flexible links and facilitates accurate modeling of real-world scenarios. It provides a comprehensive framework to accurately capture the dynamic behavior of the manipulator. Unlike decoupling methods that simplify the system by treating each link independently, the Euler–Lagrange method takes into account the interactions between all links, resulting in a more realistic representation. In addition, the equivalent rigid-link method, the Euler–Bernoulli method, and the assumed mode method may not be suitable for complex manipulator systems due to their limitations in capturing higher-order modes and non-linearities. The finite element method offers high accuracy but can be computationally intensive, making it less efficient for real-time simulations. The Hamilton method, the Newton–Euler method, and the pseudo-rigid body method also have their own advantages, but the Euler–Lagrange method balance between computational efficiency and accuracy.

The Timoshenko beam theory method focuses primarily on beam deformations, which limits its applicability to manipulator systems with other types of deformations. Finally, when modeled for control purposes, the Euler–Lagrange method is compatible with various control strategies, such as PID controllers or advanced techniques like optimal control, which makes it suitable for practical applications. Therefore, the Euler–Lagrange method with lumped parameters proves to be advantageous as it provides an accurate and efficient approach to modeling manipulators with flexible links.

A model with lumped parameters for a flexible link is a mathematical model used to signify the performance of a flexible link in a robotic or mechanical system. It approximates the flexible link as a series of discrete elements, each of which is represented by a lumped parameter such as mass, damping, or stiffness coefficient. This model is useful for predicting the flexible link behavior under varying operating conditions. In this model, it is assumed that the lumped parameter for the flexibility of the link is represented by a spring and a damper at each joint and the mass at the links’ center of gravity.

To derive the 3DOF manipulator with flexible links dynamic model applying the Euler–Lagrange technique with lumped parameters, the system parameters are defined, including the mass and inertia, the length, and the distance from the joint to the center of mass of each link. These parameters correspond to those of the test rig. The stiffness coefficients of the flexible links are determined experimentally, and the damping coefficient is assumed. For each link, the center of mass and the second moment of inertia I are calculated by accurately drawing the parts in SOLIDWORKS 2020 with the defined materials used in the test rig. The products of the inertia terms are assumed to be zero as they are relatively small compared to the principal terms. The characteristics of the motors are neglected to reduce the model complexity, and the weight of the motor is added as a part of the link. Link 2 assembly using SOLIDWORK and sectional view of the VSL are shown in Figure 2.

For the 3DOF manipulator, three generalized coordinates ($q_1,q_2,q_3$) can be defined to represent the joint angles, and three generalized velocities (${\dot{q}}_1,{\dot{q}}_2,{\dot{q}}_3$) to represent the joint angular velocities.

The flexible links are modeled as a lumped model consisting of a rigid body with an additional spring and a damper at each joint, representing the links’ flexibility, to compensate for the links flexibility in the two bending directions [47]. Axial flexibility is ignored as the fabric layer prevents the flexible link from stretching, and any stretching is very small and has no significant effect on the link. The weight is assumed to be at the center of mass of each link.

The kinetic energy is the summation of each link the kinetic energy (K). While the potential energy is the summation of the potential energy (P). For the kinetic energy, it can be declared in terms of the angular velocities and the inertia matrices of each link, meanwhile the potential energy is only a function of the gravitational constant and the center of mass height of the individual links [48].

The general form of kinetic energy is shown in (10), where I_i is the second moment of inertia of the ith link at the center of mass, J_vi is the linear velocity of the ith link, J_ωi is the angular velocity of the ith link, R_i is the rotational matrix from the link i frame to the base frame, m_i is the mass of the ith link, l_i is the length of the ith link, (X_i,Y_i,Z_i) is the position of the center of mass from the i − 1 to the ith link, ki_f are the stiffness coefficients of the flexible ith link, and ci_f is the damping coefficient of the ith link. Let $S_i=\mathrm{s}\mathrm{i}\mathrm{n}\left(q_i\right)$, $C_i=\mathrm{c}\mathrm{o}\mathrm{s}\left(q_i\right)$, $S_{ij}=\mathrm{s}\mathrm{i}\mathrm{n}\left(q_i+q_j\right)$, and $C_{ij}=\mathrm{c}\mathrm{o}\mathrm{s}(q_i+q_j).$

$$\begin{array}{cc}\hfill k& =\frac{1}{2} {\dot{q}}^2 {\displaystyle \sum _{i=1}^n} [ m_i{J}_{vi}(q{)}^T{J}_{vi}\left(q\right) + J_{\omega i}{\left(q\right)}^T{R}_i\left(q\right)I_i{R}_i{\left(q\right)}^T{J}_{\omega i}\left(q\right)]\dot{q}\hfill \\ & =\frac{\dot{1}}{2}{q}^{T}D\left(q\right)\dot{q}\hfill \end{array}$$

(13)

From Equation (13), we obtain the following:

$$\begin{array}{cc}\hfill k& ={\dot{q}}_1^2\left(m_1{X}_1^2+m_1{Z}_3^2+I_{1yy}\right)/2+I_{3zz}{\left({\dot{q}}_2+{\dot{q}}_3\right)}^2/2+I_{2zz}{\dot{q}}_2^2/\hfill \\ & +(m_2({\left({\dot{q}}_2{Y}_2-{\dot{q}}_1{Z}_2{S}_2\right)}^2+{\left({\dot{q}}_2{l}_2+{\dot{q}}_2{X}_2-{\dot{q}}_1{Z}_2{S}_2\right)}^2\hfill \\ & +{\dot{q}}_1^2{{(l}_2{C}_2+X_2{C}_2)-Y_2{S}_2)}^2+m_3(X_3\left({\dot{q}}_2+{\dot{q}}_3\right)+{\dot{q}}_2{l}_3+{\dot{q}}_3{l}_3+{\dot{q}}_2{l}_2{C}_3\hfill \\ & -{\dot{q}}_1{Z}_3{S}_{23}{)}^2)/2+m_3({\dot{q}}_2{l}_2{S}_3-Y_3\left({\dot{q}}_2+{\dot{q}}_3\right)+{\dot{q}}_1{Z}_3{C}_{23}{)}^2/2\hfill \\ & +{m_{3}q}_1^2(l_3{C}_{23}+X_3{C}_{23}+l_2{C}_2{)}^2/2+{m_{3}q}_1^2(l_3{C}_{23}+X_3{C}_{23}\hfill \\ & +l_2{C}_2{)}^2/2+{\left(I_{3yy}{\dot{q}}_1^2{C}_{23}\right)}^2)/2+(I_{2yy}{\dot{q}}_1^2{C}_2^2)/2+(I_{2xx}{\dot{q}}_1^2{S}_2^2)/2\hfill \end{array}$$

(14)

The potential energy depends on the gravitational terms and the spring terms.

$$\begin{array}{cc}\hfill P& =g{(m}_2\left(l_1+Z_2+l_2{S}_2\right)+m_1\left(l_1+z_3\right)+m_3(l_1+Z_3\hfill \\ & +l_3{S}_{23}+l_2{S}_2)+\frac{1}{2}{k}_{1f}{q}_1^2+\frac{1}{2}{k}_{2f}{q}_2^2+\frac{1}{2}{k}_{3f}{q}_3^2\hfill \end{array}$$

(15)

The Lagrangian is formulated as the subtraction of the system's kinetic energy from its potential energy, and it is expressed as a function of the generalized coordinates and their corresponding velocities [48].

$$L=K-P$$

(16)

This equation is used for the equations of motion derivation for the system and involves taking the of the Lagrangian’s partial derivatives regarding the generalized coordinates and velocities [3]. D is the Rayleigh dissipation function in this model it represents the damping.

$$D=\frac{1}{2}{c}_{1f}{\dot{q}}_1^2+\frac{1}{2}{c}_{2f}{\dot{q}}_2^2+\frac{1}{2}{c}_{3f}{\dot{q}}_3^2$$

(17)

The Lagrangian’s partial derivatives become the following [49]:

$$\frac{d}{dt}\left(\frac{\partial L}{\partial \dot{q_i}}\right)-\frac{\partial L}{\partial q_i}+\frac{\partial D}{\partial {\dot{q}}_i}={\tau }_i$$

(18)

where ${\tau }_i \mathrm{i}\mathrm{s}$ all the generalized non-conservative forces.

The resulting equations of motion consist of combined non-linear differential equations that are numerically solved to obtain the motion of the manipulator over time. The manipulators dynamic equations have the following general form [48]:

$$M\left(q\right)\ddot{q}+C\left(q,\dot{q}\right)\dot{q}+G\left(q\right)=\tau +{\tau }_f$$

(19)

where M(q) is the symmetric positive definite mass inertia matrix of the system, C(q, $\dot{q}$) is the matrix of the Coriolis and centrifugal terms, G(q) is the vector of the gravity terms, τ is the input vector, and τ_f is the torque from the springs (flexibility of the links). The detailed elements of matrixes M and C are in Appendix A.

$$M=\left[\begin{array}{ccc}{M}_{11}& M_{12}& M_{13}\\ M_{21}& M_{22}& M_{23}\\ M_{31}& M_{32}& M_{33}\end{array}\right]$$

(20)

$$C=\left[\begin{array}{c}{C}_{11}\\ C_{21}\\ C_{31}\end{array}\right]$$

(21)

$$G=\left[\begin{array}{c}0\\ g{m}_3\left(l_3{C}_{23}+l_2{C}_2\right)+g{l}_2{m}_2{C}_2\\ g{l}_3{m}_3{C}_{23}\end{array}\right]$$

(22)

$$\tau =\left[\begin{array}{c}{\tau }_1\\ {\tau }_2\\ {\tau }_3\end{array}\right]$$

(23)

$${\tau }_f=\left[\begin{array}{c}{\tau }_{1f}\\ {\tau }_{2f}\\ {\tau }_{3f}\end{array}\right]$$

(24)

Equation (19) is also called the inverse dynamic equation, and the forward dynamic equation is obtained by solving the joint acceleration ($\ddot{q}$) from Equation (19) [48], we obtain:

$$\ddot{q}=M^{-1}\left(q\right)(\tau +{\tau }_f-C\left(q,\dot{q}\right)\dot{q}-G\left(q\right))$$

(25)

3. Simulation

3.1. Simulation Environment

The previously derived and discussed kinematic and dynamic mathematical models were used to create a MATLAB model for Equations (6)–(8) and (25) of the VSL robot. Since the input to the model is torque, an angle–torque transformer must be added to compare it to the test rig, where the angle is the input and the motor’s torque cannot be controlled [49]. The modeling parameters are measured from the test rig. For each link, the center of gravity and the second moment of inertia I are calculated by accurately drawing the parts in SOLIDWORKS with the defined materials used in the test rig. The mathematical model implemented in MATLAB/SIMULINK 2022a is shown in Figure 3. The test rig parameters used for the model in SIMULINK are tabulated in

Table 2. Test rig parameters used in the modeling in SIMULINK.

Modeling Parameters
m1 = 0.77189 kg	Y1 = 0.04 m	I1yy = 1.83 × 10−3
m2 = 0.215 kg	Z1 = 0.0245 m	I1zz = 1.86 × 10−3
m3 = 0.14752 kg	X2 = 0.157 m	I2xx = 3.82 × 10−5
l1 = 0.04 m	Y2 = 0	I2yy = 1.6 × 10−3
l2 = 0.29 m	Z3 = 0	I2zz = 1.59 × 10−3
l3 = 0.240 m	X3 = 0.125 m	I3xx = 4.26 × 10−5
k1f = k2f = k3f = 1.66	Y3 = 0	I3yy = 7.45 × 10−4
c1f = c2f = c3f = 1.3	Y4 = 0	I3zz = 7.25 × 10−4
X1 = 0	I1xx = 2.47 × 10−4

.

3.2. Workspace

The workspace of the 3DOF VSL manipulator shows the maximum and minimum values that the manipulator can achieve, taking into account any joint limitations, physical obstacles, or constraints in the environment. It helps in understanding the reachable space and the potential applications and limitations of the manipulator in various tasks and operational scenarios. Figure 4 shows the 3D workspace, while Figure 5 shows a slice of the 3D workspace at x = 0 plane.

4. Experimental Validation

Two experiments were carried out to validate the mathematical models. The first experiment focused on the validation of the kinematic model, while the second experiment aimed at validating the dynamic model. Both experiments followed a similar setup, as shown in Figure 6. Two types of cameras were used in the experiments. The Olympus SP-320 camera was used in the experiments for kinematic validation, while a Sony PlayStation Eye camera (30 frames per second) was used for the experimental validation of the dynamic model. Both cameras were used in the same two positions for data recording.

In the first position, the camera was positioned in front of a tripod to record the movements of the second joint (θ₂), the third joint (θ₃), and the end tip position in the X-Z plane. The second position of the camera was from the top to capture the movement of the first joint (θ₁) and the end-tip position in the X-Y plane.

4.1. Experimental Test Rig

A 3 DOF VSL manipulator prototype was developed in our laboratory from a combination of silicone rubber structures, polypropylene plastic mesh, and fabrics (see Figure 7). VSLs have the ability to control their stiffness and are pneumatically actuated. By adjusting the pressure value within their structure, the links can be tuned to different stiffness levels. Our prototype is equipped with a pressure sensor (BD SENSORS 26.600 G-11002-R—1-5-100-300-1-000) to sense the robot’s body possible collisions with a human. The base and elbow motors were serial servomotors with built-in encoders to measure the angle of rotation. The shoulder motor was a stepper motor with a gearbox and an encoder connected to the shaft. All experiments were carried out at a pressure of 2 bar in the links. This pressure was selected experimentally. The highest pressure for maximum stiffness and, at the same time, the pressure change at collision could be determined.

4.2. Kinematic Model Experimental Results

The first experiment confirmed the end effector’s position based on derived from the forward kinematics derived from the DH parameters. The VSL robot prototype was subjected to various poses, and photos were taken of each pose. Using the SOLIDWORKS 2020 program, the different angles for each pose were measured, along with the end effector’s position.

The photos of each pose were inserted into a 2D drawing and measured manually (see Figure 8 and Figure 9). These angles were then inserted into the forward kinematics equations to calculate the end effector position. That calculated values for the end-effector position were then compared with the values obtained from the photos.

The first experiment results show an average error of 2.90 mm in the x-direction, 2.29 mm in the y-direction, and 2.40 mm in the z-direction, so there is still room for improvement. To achieve better results, more precise assembly and alignment are required. It should also be noted that links 2 and 3 of the prototype were slightly curved. These adjustments can ensure that the system reaches its full potential and achieves an even higher level of accuracy and precision. The plots for the experimental manipulator‘s end tip position and the calculated position for different poses along the X-axis, Y-axis, and Z-axis are shown in Figure 10, Figure 11, and Figure 12, respectively.

Table 3. The data along the X-, Y-, and Z-axes and the error between the theoretically calculated values and experimentally measured values.

	θ1	θ2	θ3	Experimental (mm)			Theoretical (mm)			Error (mm)
				x	y	z	x	y	z	x	y	z
1	88.69	71.61	−86.82	-	326.45	262.15	-	324.66	256.96	-	1.79	5.19
2	88.69	45.48	−74.17	-	414.1	137.76	-	417.37	135.11	-	3.27	2.65
3	88.69	85.41	−2.27	-	51.9	565.35	-	52.27	572.33	-	0.37	6.97
4	88.69	5.74	−52.09	-	456.38	−101.6	-	459.18	−104.1	-	2.8	2.57
5	88.69	6.01	−6.01	-	527.89	70.54	-	533.37	70.88	-	5.48	0.33
6	88.69	6.17	−85.62	-	336.84	−160.7	-	337.23	−164.2	-	0.39	3.51
7	88.69	17.48	−6.9	-	514.11	171.89	-	517.29	172.67	-	3.18	0.78
8	88.69	17.54	−69.1	-	428.36	−54.91	-	430.49	−59.07	-	2.13	4.16
9	88.69	29.9	−5.93	-	472.4	283.21	-	475.03	284.55	-	2.63	1.34
10	88.69	29.7	−39.05	-	489.56	148.71	-	493.05	147.16	-	3.49	1.55
11	88.69	45	−4.17	-	388.62	403.97	-	390.19	405.51	-	1.57	1.54
12	88.69	45.78	−85.86	-	386.48	100.25	-	389.37	96.89	-	2.89	3.36
13	88.69	58.52	−4.14	-	292.44	484.43	-	293.82	486.67	-	1.38	2.24
14	88.69	58.35	−43.16	-	384.95	353.38	-	386.4	354	-	1.45	0.62
15	88.69	72.04	−3.98	-	180.18	541.9	-	180.63	543.24	-	0.45	1.34
16	88.69	71.44	−73.11	-	330.83	313.12	-	333.79	312.66	-	2.96	0.46
17	88.69	85.61	−24.8	-	139.88	540.29	-	139.63	543.65	-	0.25	3.36
18	88.69	85.07	−85.07	-	262.85	332.92	-	265.35	333.9	-	2.5	0.98
19	88.69	1.88	6.6	12.13	530.63	-	12.23	534.86	-	0.1	4.23	-
20	75.41	1.88	6.6	134.46	516.42	-	134.76	517.74	-	0.30	1.32	-
21	60.24	1.88	6.6	263.01	459.91	-	265.55	464.44	-	2.54	4.53	-
22	43.28	1.88	6.6	385.56	363.09	-	389.48	366.77	-	3.92	3.68	-
23	26.76	1.88	6.6	474.75	239.37	-	477.7	240.88	-	2.95	1.51	-
24	10.56	1.88	6.6	519.42	101.38	-	525.93	98.04	-	6.51	3.34	-
25	−3.09	1.88	6.6	530.18	−28.58	-	534.22	−28.83	-	4.04	0.25	-
26	88.69	1.88	6.6	-	532.18	30.27	0	533.88	29.92	0	1.7	0.1225
Average error										2.90	2.29	2.40
Average percentage error										0.81%%	0.71%	1.4%

lists the errors along the X-, Y-, and Z-axes, which indicate the theoretically calculated values deviations from the corresponding experimental measurements.

4.3. Dynamic Model Experimental Results

In this experiment, the VSL robot was given a pose, followed by the recording of its movement from its pose to the given pose at a frame rate of 30 fps. First, the second joint (θ₂) was moved from 90° to 16°, then, the third joint (θ₃) was moved from 0° to −45°, and, finally, the first joint was moved from 0° to −43°. Then, using the SOLIDWORKS program, the photos of each frame were inserted into a 2D drawing, the angles of the robot’s movement were measured manually, and the final position of the tip was tracked, frame by frame, from the beginning until it came to a complete stop. Examples of these images are shown in Figure 13, Figure 14 and Figure 15. The same pose was then applied to the model in SIMULINK and its motion was also recorded. Finally, the collected experimental and simulated data were recorded, plotted, and compared.

The results of the experiment show that, although the model’s responses were adequate, they did not have the degree of accuracy we would have liked. The average error of the angular displacement for the first joint was 0.63°, for the second joint 1.06°, and for the third joint 1.27°. In order to achieve greater precision, it is essential that we carry out further investigations into the damping and stiffness coefficients. The plots for the three angles of the joints versus time for the experimental and simulated responses of the manipulator are illustrated in Figure 16, Figure 17, and Figure 18, respectively.

In general, the plots show that the simulation model captured the whole trends of the 3DOF flexible link manipulator experimental joint angles data. However, there were some discrepancies between the experimental and simulated responses that could be attributed to various factors such as modeling assumptions (e.g., damping coefficient), measurement noise, or unmodeled dynamics (e.g., motors). The three recorded data plots for the three joint angles showed some degree of oscillation, which was due to the flexibility of the links. For θ₃ in Figure 18, there were more oscillations just before the steady state than for θ₂ in Figure 19, which was due to the second link flexibility effect on the third link, while for θ₂ in Figure 19, there was only the second link flexibility effect.

The oscillations of the links decrease or disappear with higher rotational speed when the centrifugal force exceeds the flexibility torque but increases the overshoot when it stops. The stronger the acceleration or deceleration, the more oscillations occur at the beginning and end of the movement.

The results of our research show that the end position curve obtained from the experiments with a camera and the positions from the simulated model agree well, as shown in Figure 19. The average errors in the x, y, and z directions were 6.24 mm, 6.4 mm, and 2.93 mm, respectively, demonstrating reasonable accuracy between the experimental and simulated data. However, further analysis of the model is required to increase accuracy. A sample of the experimental and simulated data and the error between them can be found in

Table 4. Sample of the data and dynamic errors along the x-, y-, and z-axes between the simulation values and those measured experimentally.

	Time (sec)	Experimental (mm)			Theoretical (mm)			Error (mm)
		X	y	z	x	y	z	x	y	z
1	1.7	278.442	0	492.48	285.39	0	486.59	6.95	0	5.88
2	1.73	281.2061	0	491.55	295.12	0	480.22	13.91	0	11.33
3	1.76	305.3148	0	474.64	312.81	0	467.84	7.49	0	6.8
4	1.8	320.0337	0	464.22	323.09	0	460.13	3.06	0	4.09
5	1.83	330.6901	0	455.15	331.74	0	453.33	1.05	0	1.82
6	1.86	344.3517	0	445.08	348.3	0	439.47	3.95	0	5.61
7	1.9	354.0388	0	437.74	358.58	0	430.27	4.54	0	7.46
8	1.93	365.542	0	425.37	371.04	0	418.45	5.5	0	6.91
9	1.96	375.4381	0	416.44	376.71	0	412.8	1.27	0	3.63
10	2	386.3407	0	405.04	388.09	0	400.94	1.75	0	4.1
11	4	510.1	0	163.63	515.59	0	159.16	5.49	0	4.47
12	4.13	519.23	0	155.47	516.45	0	153.69	2.78	0	1.77
13	4.16	513.53	0	153.83	517.52	0	145.14	3.99	0	8.68
14	4.2	517.64	0	152.18	518.09	0	139.2	0.45	0	12.97
15	4.23	515.41	0	143.87	518.25	0	137.1	2.83	0	6.76
16	4.26	514.45	0	138.89	518.36	0	135.56	3.91	0	3.32
17	4.3	520.68	0	137.01	518.56	0	132.18	2.11	0	4.82
18	4.36	518.12	0	123.38	518.82	0	125.77	0.7	0	2.39
19	4.4	510.1	0	163.63	515.59	0	159.1	5.49	0	4.47
20	7.13	484.12	0	2.81	488.66	0	2.81	4.53	0	0
21	7.16	477.51	−12	2.81	488.52	−11.7	2.81	11	0.29	0
22	7.20	483.1	−27	2.81	487.76	−29.7	2.81	4.65	2.65	0
23	7.23	491	−40.8	2.81	486.91	−41.3	2.81	4.18	0.49	0
24	7.26	485.8	−50.1	2.81	485.81	−52.7	2.81	0.05	2.59	0
25	7.30	481.83	−63.1	2.81	483.54	−70.5	2.81	1.71	7.39	0
26	7.33	485.17	−83.1	2.81	491.11	−80.5	2.81	5.93	2.65	0

.

Moreover, the system variables’ identification, especially the stiffness and damping coefficients related to flexibility, is essential for the development of a more accurate model. Taking these key factors into account leads to a more accurate representation of the system and, therefore, even more reliable results. By investing the necessary effort and resources in improving these modeling techniques, we can ensure that the system meets the highest standards and fulfills the most stringent requirements.

5. Conclusions

In conclusion, the development of VSL manipulators represents a significant advance in robotics, offering greater adaptability and safety in a range of applications. The VSL manipulator’s ability to adjust its stiffness in real time increases precision and efficiency, making it an ideal choice for tasks such as assembly line work. However, the non-linear dynamics of the VSL manipulator, resulting from the flexibility of its links, pose a challenge in modeling and control. This requires advanced techniques and leads to increased computational costs.

The mathematical model presented for this 3 DOF VSL manipulator is a first step toward enhancing its performance in order to make the manipulator more feasible, increase safety, and reduce costs. The verification experiments have proven the model reliability and accuracy and underline the importance of mathematical modeling in engineering and scientific research. The kinematic model shows an average error percentage of 0.81% along the X-axis, 0.71% along the Y-axis, and 1.4% along the Z-axis for the end-tip position. The dynamic model has an average error of 0.63° for the first joint angular displacement of 1.06° for the second joint, and 1.27° for the third joint. The average percentage error in the position of the end tip along the X-, Y-, and Z-axes is 3.36%, 3.29%, and 2.49%, respectively.

The results of the present study dynamic model of the present study exceed the results reported in [50]. In [50], a discrete Newton–Euler model is developed to analyze flexible manipulators, focusing on the fundamental dynamics of generic links. The rigid body motion is described using the Eulerian formulation, while elastic deformation is described using the Lagrangian formulation. The Jourdain principle guides the assumption of the Rayleigh–Ritz expansion for elastic variables. Manipulator dynamics are derived by assembling individual links, which leads to both articulated body (AB) and composite inertia (CI) methods. Only a single flexible link is used for validation. The maximum error reaches 44.6% at some points at the beginning of the fluctuation and 22% at some points towards the end [50], while the derived model in our study has a maximum error of about 18% at the beginning and an error of about 3% thereafter.

The results and analysis of this study offer significant insights into the kinematics and dynamics of VSL manipulators. Future work will aim to implement a control algorithm tailored to the unique characteristics of this manipulator, enabling precise and adaptive manipulation in dynamic and uncertain environments.