Modeling and validation of a serial-manipulator tubular bending robot for complex pipe forming

3.1. Forward kinematics solution of pipe bending robot

Forward kinematics involves determining the pose [15], [16] of the robot’s end-effector by knowing the joint angles. Once the kinematic model is established, it must be linked with the coordinate systems. The transformation between adjacent coordinate systems is typically achieved through translation and rotation. The tubular bending robot is a series connection of the end tool shaft and a standard six-axis robot arm, with a proportional transmission relationship between the sixth axis and the end tool shaft. Thus, the forward kinematics solution for the tubular bending robot can be treated as that of a six-axis robot arm, with the homogeneous transformation matrix of the end tool shaft connected in series. The steps and sequence for this process, as outlined by the improved D-H parameter method, are as follows.

(1) The ${\alpha }_{i-1}$ angle of rotation of the $Z_{i-1}$ axis around the $X_{i-1}$ axis, so that the $Z_{i-1}$ axis is parallel to the $Z_i$ axis.

(2) The origin of the coordinate system $\left\{O_{i-1}\right\}$ is translated a distance $a_{i-1}$ along the $X_{i-1}$ axis, so that the $Z_{i-1}$ axis is collinear with the $Z_i$ axis.

(3) The $X_{i-1}$ axis is rotated about the $Z_i$ axis through an angle ${\theta }_i$, such that the $X_{i-1}$ axis becomes parallel to the $X_i$axis.

(4) The origin of coordinate system $\left\{O_{i-1}\right\}$ is translated along the $Z_i$ axis by a distance $d_i$, such that coordinate system $\left\{O_{i-1}\right\}$ completely coincides with coordinate system $\left\{O_i\right\}$. Therefore, the transformation matrix between coordinate systems $\left\{O_{i-1}\right\}$ and $\left\{O_i\right\}$ is given by:

By bringing the parameters in Table 1 into Eq. (2), the transformation matrix between the joints of the robot can be calculated as:

The total transformation matrix of the forward kinematics of the bending robot is the Eq. (11):

3.2. Solving the inverse kinematics of the pipe bending robot

The inverse kinematics is the basis of robot control. Different from the forward kinematics of robots, the inverse kinematics solution is to solve the angle [17], [18] of each joint when the position and pose of the end effector of the robot are known. Inverse kinematics solutions typically involve analytical, numerical, and geometric approaches. The bending robot in this paper has seven joint axes, and its configuration does not meet the Pieper principle. If the inverse solution is directly solved for the whole, only numerical method can be adopted, which is of great computational complexity, long time and high requirements on the initial value of iteration, and is not suitable for real-time control requirements [19], [20]. A solution method using joint angle parameterization is proposed in this paper to address the robot's inverse kinematics. This novel parameterization strategy effectively reduces the system's dimensionality, transforming a complex multi-axis configuration into an analytically solvable model, providing a computationally efficient alternative to conventional iterative numerical solvers.

The 7th joint of the bending robot is driven by the 6th joint, with its output angle amplified by the rotator. This results in a proportional relationship between their angles, allowing them to be treated as a single unit during decoupling calculations. To enable flexible processing, the angles of the sixth and seventh joints are customized for attitude control. The angles of the remaining joints are computed analytically to solve the robot's inverse kinematics.

Among them, the sixth joint and the seventh joint are the custom angle, and the expression of its transformation matrix can be obtained as follows, with $s\theta =\sin \theta$ (abbreviation of $\sin \theta$), $c\theta =\cos \theta$ (abbreviation of $\cos \theta$), $k$ is the gear ratio coefficient:

The end-tool matrix can be obtained according to the D-H parameter, and its expression is:

Apply the inverse matrix transformation to the Eq. (11), you can get:

Eq. (15) is the total homogeneous transformation matrix of the first to the fifth joint. The solution of this matrix should first solve the angle of the first joint, and also use the inverse matrix transformation to rewrite Eq. (15) into the following form:

Of which:

Let the (2, 4) element in the matrix product result of Eq. (16) correspond to the (2, 4) element in Eq. (17), we can get:

The solution Eq. (18) can be solved through the trigonometric identity transformation formula to get:

At this time, the ${\theta }_1$ that has been solved has two sets of solutions, and then solve ${\theta }_5$, so that (2, 1) and (2, 2) in Eq. (16) and Eq. (17) are equal, we can get:

In combination with the above Eq. (20) and Eq. (21), it can be solved:

Corresponding to different ${\theta }_1$, ${\theta }_5$ with corresponding solutions. Put the calculated ${\theta }_5$ into the corresponding homogeneous transformation matrix, and rewrite the Eq. (16) as follows by the inverse matrix transformation:

In Eq. (23), the sine function $s_{\phi }$ represents the sum of angles of joints 2, 3 and 4. The cosine function $c_{\phi }$ representing the sum of angles of joints 2, 3, and 4. While $x$ and $-z$ representing the $X$ and $Z$ axis positions in the operating space:

This is obtained by adding the above Eq. (26) to the Eq. (27) simultaneously squared:

By utilizing the sum-angle formulas of sine and cosine functions, $\sin {\theta }_3$ and $\cos {\theta }_3$ can be determined. Subsequently, ${\theta }_3$ is calculated using the two-variable arctangent formula:

$\sin {\theta }_3$ has two values (positive and negative), there exist two solutions for ${\theta }_3$. Once the values of ${\theta }_3$ are determined, Eq. (26) and Eq. (27) can be reformulated into the following forms:

In Eq. (30) and Eq. (31), the specific expressions of $k_1$ and $k_2$ are as follows:

Through variable substitution, the expression of ${\theta }_2$ can be obtained as follows:

Finally, given ${\theta }_2$ and ${\theta }_3$, ${\theta }_4$ can be determined using the following equation:

Thus, with ${\theta }_1$, ${\theta }_2$, ${\theta }_3$, ${\theta }_4$ and ${\theta }_5$ calculated and combined with the user-defined parameters ${\theta }_6$ and ${\theta }_7$, the inverse kinematics solution for the pipe-bending robot is completed.

3.3. Singularity avoidance strategy

Furthermore, it is crucial to address the kinematic singularities commonly encountered by 6-DOF and 7-DOF robots. The analytical inverse kinematics approach proposed above inherently yields multiple solution branches. During practical trajectory planning, the controller monitors the Jacobian matrix and utilizes these multiple analytical solutions to actively switch the manipulator’s configuration, thereby effectively avoiding kinematic singularities and dead zones within the robot’s workspace.

4.1. Pipe bending robot simulation model construction

The construction of the simulation model of the pipe bending robot adopts MATLAB Multibody of MATLAB Link plug-in, which exports the 3D model of Solidworks, Inventor, Creo and other modeling software into step and xml format files that can be read by MATLAB, and effectively integrates the 3D model into MATLAB for simulation and analysis [21], [22]. In the early Solidworks to create the robot model, it must be ensured that the constraints between each joint are configured for the correct coincidence and coaxial relationship, and the initial state of the assembly should meet the state requirements of the D-H parameter, so that each joint axis [23], [24] can be accurately identified when exported. The rotating head of the bending robot consists of numerous complex components, which require significant computational resources for simulation model construction. To enhance simulation efficiency, components like the bending spindle, clamping mechanism, and guide mechanism are simplified and replaced with a connecting rod model.

After building the simulation model of the robot, the next key step is to build the motor model to drive the movement of the robot joints. This paper uses the stepper motor model in Simulink, which can be quickly created by inputting the “elec_motor_stepper” command in the MATLAB command window. The motor model provides two control modes, one is position control, the other is speed control. Since the input signal received by the robot model built in the above section is angle input, the position control mode should be selected for the motor model to meet the requirements of the robot joint input.

In the input end of the motor model, in order to directly import the results of solving the inverse kinematics of the robot and the trajectory planning, the default step signal input needs to be modified to the self-defined system input signal Input1. Considering that the entire step distance of the stepper motor used is 1.8°, it is necessary to divide the angle value of the solution result by the step distance of the motor, so as to obtain the number of steps used to drive the motor, and then input the data into the motor controller in order to achieve accurate control of the robot joint [25].

At the output end of the motor model, an adder is designed to add the compensation value of the initial angle of the robot joint, so that the initial attitude of the robot can be customized. Considering that the output of the motor model is in the unit of angle, and the robot joint needs the input in radian, so the necessary conversion components are added between the two. Finally, with the help of Simulink-PS conversion element, the Simulink input signal is converted into the physical signal output, and the MATLAB simulation model is built, as shown in Fig. 5.

Fig. 5MATLAB model diagram of the bending robot

Based on the above MATLAB simulation model and motor drive model, the model can be used to verify the proposed inverse kinematics solution method based on joint angle parameterization. Specifically, the known end-pose parameters of the robot, including position $(X,Y,Z)$ and attitude $(Roll,Pitch,Yaw)$, need to be given, and converted into a one-dimensional array and input to the MATLAB Function module. The code of inverse kinematics solving algorithm has been written inside the module, which is used to calculate the corresponding joint angle. These angles are then output to a motor driven model that is then passed on to the robot. Furthermore, to validate the robustness of the proposed method under complex spatial conditions, extensive simulation tests involving continuous multi-bend trajectories were conducted. The results demonstrated that the analytical solver consistently maintained high positional accuracy and real-time computational efficiency throughout the complex bending sequences.

4.2. Pipe bending robot experiment platform construction

The experimental platform model includes multiple modules, as illustrated in Fig. 6.

(1) Instruction receiving module: This module contains two key sub-modules, namely, data receiving module and data verification module.

(2) Motion control module: In the motion control module, three different control strategies are implemented to adapt to different operational needs. They are point to point (PTP) control mode, joint space trajectory planning (PVT) mode and Cartesian space trajectory planning (PT) mode.

(3) Joint position module: This module realizes the data storage function, which is mainly used to store the joint angle data processed by the motion control module, so as to facilitate the robot model to read it.

(4) Data broadcast module: The data broadcast module provides a communication mechanism for the robot to distribute or share information among multiple systems, modules or devices.

The operation procedure for the experimental platform is conducted through a systematic workflow. The experimental platform model is opened in MATLAB, and the appropriate control mode is selected within the motion control module. For this study, the Cartesian space trajectory planning (PT) mode is adopted as the operating example.

Using the written MATLAB script, the joint angle sequence is computed through inverse kinematics for each interpolation point of the Cartesian space trajectory. The simulation is executed to verify that all robot joints can move continuously without encountering crashes or locking issues. Upon successful simulation, the computed joint angle values are stored in the joint position module of the experimental platform model.

During compilation and download of the platform model, the corresponding joint angle data are automatically extracted from the script workspace and transmitted to the robot motors for execution. The hardware components are then connected in sequence. Power is supplied to the voltage regulator, with the voltage adjusted to 48V (48.0 V~48.6 V). Subsequently, the robot controller is activated to supply power to the robot. The compile button at the top of the model interface is clicked to download the model and data to the manipulator controller. The successful program installation is verified through the View Diagnostic Viewer.

The control system software is executed following the operational sequence: Start → Reach initial attitude → Start trajectory movement → Return to zero → Stop. This systematic approach ensures proper initialization, execution, and termination of the robotic bending operations.

Fig. 6Experimental platform of bending robot. Photo by the second author in the laboratory of the School of Engineering, Huzhou Normal University in March 2025

4.3. Kinematic experiment and analysis of pipe bending robot

Based on the simulation model and experimental platform of the pipe bending robot established above, the kinematic experiment analysis is carried out. Specifically, it is necessary to give the known robot end pose parameters including position $(X,Y,Z)$ and attitude $(Roll,Pitch,Yaw)$, and convert them into a one-dimensional array and input them into MATLAB Function module. The code of inverse kinematics solving algorithm has been written inside the module, which is used to calculate the corresponding joint angle. The experiment is mainly to verify the accuracy of the proposed inverse kinematics method for solving the bending robot.

Taking a single pose point as an example, firstly, given the end-effector pose of the bending robot as (452 mm, –173 mm, 482 mm, 60°, 0°, 90°), and the initial angle of the sixth joint as –23.64°, the angles of all joints of the robot are calculated as (0°, –31.24°, –27.32°, –36.46°, 0°, –23.64°) using the inverse kinematics solution method described above. Then, verify the accuracy of the calculated results, input these joint angles into the simulation model, and observe whether the pose of the bending robot is reasonable after the completion of movement.

Finally, the given pose is tested on the real machine of the pipe bending robot, and the actual angles of each joint of the robot can be obtained as (0.01°, –31.14°, –27.52°, –36.37°, 0°, –23.73°), and the actual pose of the end of the robot is (451.68 mm, –172.86 mm, 483.56 mm, 60.26°, 0°, 89.82°), as shown in Table 2 and Table 3. The operation results are shown in Fig. 7. Fig. 7 respectively represent the initial pose, intermediate pose and end pose of the bending robot in the actual experiment.

Fig. 7Real motion pose of pipe bending robot

a) Starting pose

b) Motion process pose 1

c) Motion process pose 2

d) End pose

In the accurate test of the pipe bending robot, a specific point is chosen as the origin for pipe fitting processing, and various spatial positions are tested. Five sets of end poses are used to gather comparison data for the theoretical, simulation, and actual values of the robot's end pose, as shown in Table 4.

The kinematics error of the pipe bending robot’s end pose was analyzed by testing different processing positions to assess the accuracy of the kinematics solution. The maximum error between the theoretical and actual values is 0.29 mm, and the minimum is 0.05 mm. It is worth noting that this maximum error of 0.29 mm is partially attributed to the inherent mechanical backlash within the gear transmission system connecting the sixth and seventh joints. Although precision gears were utilized, micro-clearances during direction reversals slightly affect the absolute positioning accuracy of the end-effector.

Fig. 8Comparison data of error values of the theory, simulation and actual values of the pose of the bending robot

For simulation versus actual values, the maximum error is 0.31 mm and the minimum is 0.04 mm. The maximum rotation angle error between theory and actual values is 0.24°, and the minimum is 0.07°. Between simulation and actual values, the maximum error is 0.23° and the minimum is 0.065°. These results show that the motion error of the robot increases with the movement distance, but all errors remain within acceptable limits, satisfying the robot's accuracy requirements. The comparison of theoretical coordinates, rotation angles, and errors between simulation and actual values are presented in Fig. 8.