Modeling and control strategy of a haptic interactive robot based on a cable-driven parallel mechanism

This study proposes a haptic interactive robot (HIR) configuration and a control strategy based on a cable-driven parallel mechanism. The ball screw drives the cable to improve the motion control accuracy. The robot system control strategy improves the accuracy and stability of haptic interaction. Through configuration optimization design and analysis, eight cables are used to ensure that the robot end effector exerts force and enables motion. Moreover, a forward and inverse kinematics model of the robot is developed. According to the configuration of the HIR, an improved cable tension distribution algorithm can facilely determine the cable tension. Hence, each cable is consistently in a tight state, and the change in tension is not sudden. Drive unit and robot system control strategies are proposed to render the haptic interaction accurate and stable. A simulation experiment of a complex space motion track is implemented through the robot end effector, thus verifying the accuracy of the established forward and inverse kinematics model. The accuracy of the tension distribution algorithm, control strategy, and robot stability are verified through simulation experiments, considering different forces and motion tracks of robot end effectors.

1 Introduction

Compared with rigid parallel robots, a cable-driven parallel robot (CDPR) is characterized by its light weight, large bearing capacity, fast speed, and rapid reconfiguration. Because of these advantages, many scholars have recognized the importance of CDPRs and have applied them to rehabilitation training robots (Hamida et al., 2021; Ghrairi et al., 2023; Chen et al., 2019; Duan et al., 2015), radio telescopes (Yao et al., 2010; Duan et al., 2011), force interactive robots (i.e., haptic interactive robots, HIRs) (Fortin-Côté et al., 2014; Meziane et al., 2018), camera robots (Pannequin et al., 2020), handling robots (Gosselin, 2014; Lu et al., 2017), and three-dimensional (3D) printing mechanisms (Amare et al., 2019; Qian et al., 2020; Zi et al., 2019). The application of the above CDPR adopts the method of winding the drum with the cable.

Since the cables can only bear tension, they must be in tension when the CDPR operates. Therefore, in general, a CDPR with n degrees of freedom (DOF) requires at least n+1 cables to drive, or one cable is replaced by the action of gravity on the end effector (Hussein et al., 2021; Chawla et al., 2023; Aflakian et al., 2019; Jamshidifar et al., 2017; Cho et al., 2019; Ouyang and Shang, 2016). The CDPR can be classified into two categories: suspended cable parallel robot and constrained or fully constrained cable parallel robot. Scholars in the field of CDPR have focused their research on cables on various aspects, such as trajectory planning, kinematics, workspace, control, and collision (Gouttefarde et al., 2015; Zake et al., 2019; Shang et al., 2020; Zarafshan et al., 2022).

They have also proposed a numerical solution approach for the forward and inverse kinematics of the CDPR through an iterative operation. However, the convergence of the solution is difficult to determine using an iterative algorithm, and the stability of the analytical solution is difficult to ensure (Pott and Schmidt, 2015; João et al., 2022; Hao et al., 2021). For the tension distribution problem in CDPRs, scholars have proposed obtaining the numerical solution by setting the objective function of the iterative operation or solving the tension feasible region by computational geometry (Fang et al., 2004; Sun et al., 2021; Cui et al., 2019). Consequently, accounting for the stability of the change in cable tension becomes less necessary. To improve the accuracy of CDPR motion, scholars have analyzed the variation in cable length and mass caused by winch winding. However, the calculation process is complex and computation is time-consuming (Ryan and James, 2017; Rushton et al., 2019; Rijk et al., 2018); therefore, many compensation variation factors (such as cable mass and deformation) must be considered. Moreover, the controller is complex, and stability is uncertain (Picard et al., 2021; Korayem et al., 2020; Nguyen and Caverly, 2021).

To improve the HIR motion flexibility, stability of cable tension distribution, control accuracy, and haptic interaction accuracy, an optimal HIR configuration based on parallel cables is proposed herein. The control strategy of driving the cable and drive unit through the ball screw improves the motion control accuracy of the robot. By optimizing the configuration of the end effector, the isotropy and flexibility of the end effector are guaranteed. A forward and inverse kinematics model is also devised, and the cable tension distribution algorithm is improved. Furthermore, a control strategy for the drive unit and robot system is proposed to improve the accuracy and stability of the robot interaction force.

The rest of this paper is organized as follows. Section 2 discusses the optimization of the HIR configuration. Section 3 elaborates on the formulation of the forward and inverse kinematics model of the new robot. The correctness of the kinematics is verified by simulation experiments. Section 4 presents the development of the dynamic model of the robot, and a cable tension distribution algorithm is proposed. In Sect. 5, the driving unit and robot system control strategies are proposed. Simulation experiments verify the accuracy and stability of robot interaction force. Finally, the conclusions are summarized in Sect. 6.

2 HIR configuration

The HIR based on a cable-driven parallel mechanism is composed of a fixed platform, a ball screw drive unit, an end effector, and a transmission device. The ball screw drives one end of the cable to move in a straight line, and the other end is connected to the end actuator. The cable transfers movement and force to the end actuator to realize its free, stable, and flexible movement in multiple spaces.

2.1 General robot configuration

The general robot model is shown in Fig. 1. Assume that the global coordinate system is (O, x, y, z); the local coordinate system is (O₁, x, y, z); A_i is the output point of the cable on the fixed platform; B_i is the connection point between the cable and end actuator; l_i is the length of the cable; a_i is the vector of the cable output point on the fixed platform in the global coordinate system; b_i is the vector from the end actuator center to the cable connection point; p is the vector in the local coordinate system, O₁, at the end actuator center in the global coordinate system, O; and the attitude angle is given by (α, β, γ).

Figure 1 Schematic diagram of an n-DOF robot with m cables.

2.2 Optimization of configuration topology

Since trainers must stand inside the fixed platform for operation training, the arms must be able to move conveniently in the space without interference from the cables. This improves the mobility of the end effector, reduces the interference between the cables and the end effector, and decreases the difficulty in implementing control.

The dimension synthesis optimization strategy reported in Hamida et al. (2021) is used to optimize the configuration of the end effector. The size of the mechanism is determined according to the working requirements of the HIR. To maximize the spatial rotation angle of the end effector and minimize the interference among cables, the objective function and constraints are set as follows:

$$\begin{array}{}\text{(1)}& \left\{\begin{array}{l}\mathrm{maximize}\mathit{\varpi }\left(\mathbf{G}\right)\\ \mathrm{subject}\mathrm{to}{H}_k\left(\mathbf{G}\right)\le \mathrm{0}\\ G=[x_{\mathrm{1}},\mathrm{\cdots },x_m]\\ x_i^{min}<x_i<x_i^{max}\end{array}\right.,\end{array}$$

where ϖ(G) is the spatial rotation angle of the end actuator, H_k(G) is the distance between the cable and boundary of the end actuator, $G=[x_{\mathrm{1}},\mathrm{\cdots },x_m]$ is the change in the value of the hinge point of the end actuator, and ($x_i^{min}$, $x_i^{max}$) is the variable search range.

The side length of the square fixed platform is set considering the operation and sensor installation size requirement, Γ = 4000 mm. As shown in Fig. 2, the initial shape of the end actuator is a cylinder (radius (r) = 300 mm; height (h) = 200 mm). The optimization is implemented according to the scale synthesis optimization algorithm. The hinge point varies along the direction indicated by the arrow in Fig. 2.

Figure 2 Initial shape of end actuator.

By approximating the configuration and rounding the side length, the configuration of the end effector is obtained, as shown in Fig. 3. The upper and lower end faces are regular pentagons with a side length, Ω, of 500 mm. The top view in Fig. 3a shows that centroids coincide. The 3D configuration is shown in Fig. 3b.

Figure 3 Optimized structure of end effector.

Figure 4 Cable layout scheme.

Based on the optimized configuration of the end effector, two cable layout schemes are shown in Fig. 4. To allow the end effector to rotate more freely, the scheme shown in Fig. 4a is selected as the final configuration scheme for the HIR.

Table 1 Local positions of eight anchor points (unit: millimeters).

The positions of the hinge point between the cable and end actuator of the force feedback robot are summarized in Table 1.

3 Kinematics model

The HIR requires high control accuracy. It must solve the forward and inverse kinematics of the robot and cable tension. The robot kinematics are inverse and forward. Inverse kinematics determines the length of the cable by determining the position and attitude of the center of mass of the end effector. Forward kinematics solves the position and attitude of the center of mass of the robot end effector with the known cable length.

3.1 Inverse kinematics model

The movement of the cable is driven by the ball screw; hence, the cable is not wound around the drum. Moreover, the mass and elasticity of the cable slightly change; thus, the influence of these two factors is not considered. The robot configuration is shown in Fig. 5. The length of the cable is measured by the encoder on the drive unit, and the interaction force exerted by the human body on the end effector is detected by the six-axis force sensor.

Figure 5 Space configuration of HIR.

The coordinate system, O–xyz, is established at point A₁, and O₁–xyz is established at the center of the mass of the end actuator; this centroid is at position χ₀. According to the layout of this new configuration scheme, the inverse kinematics model of the robot is established as

$$\begin{array}{}\text{(2)}& {\displaystyle }{\mathbf{L}}_i={\mathbf{a}}_i-{\mathbf{Rb}}_i-\mathbf{p},\text{(3)}& {\displaystyle }{\mathbf{u}}_i={\displaystyle \frac{{\mathbf{L}}_i}{l_i}},\text{(4)}& {\displaystyle }\begin{array}{rl}\mathbf{R}& ={\mathbf{R}}_x\left(\mathit{\alpha }\right){\mathbf{R}}_y\left(\mathit{\beta }\right){\mathbf{R}}_z\left(\mathit{\gamma }\right)\\ & =\left[\begin{array}{ccc}\cos \mathit{\beta }\cos \mathit{\gamma }& \sin \mathit{\alpha }\sin \mathit{\beta }\cos \mathit{\gamma }-\cos \mathit{\alpha }\sin \mathit{\gamma }& \cos \mathit{\alpha }\sin \mathit{\beta }\cos \mathit{\gamma }+\sin \mathit{\alpha }\sin \mathit{\gamma }\\ \cos \mathit{\beta }\sin \mathit{\gamma }& \sin \mathit{\alpha }\sin \mathit{\beta }\sin \mathit{\gamma }+\cos \mathit{\alpha }\cos \mathit{\gamma }& \cos \mathit{\alpha }\sin \mathit{\beta }\cos \mathit{\gamma }-\sin \mathit{\alpha }\sin \mathit{\gamma }\\ -\sin \mathit{\beta }& \sin \mathit{\alpha }\cos \mathit{\gamma }& \cos \mathit{\alpha }\cos \mathit{\beta }\end{array}\right].\end{array}\end{array}$$

The derivative of Eq. (2) with respect to time is

$$\begin{array}{}\text{(5)}& {\dot{l}}_i{u}_i+l_i{\dot{u}}_i={\dot{\mathbf{a}}}_i-\dot{\mathbf{R}}{\mathbf{b}}_i-\mathbf{R}{\dot{\mathbf{b}}}_i-\dot{\mathbf{p}},\end{array}$$

where $\dot{\mathbf{R}}{\mathbf{b}}_i=\mathit{\omega }\times \mathbf{R}$, and ω is the angular velocity of the end actuator around the coordinate axis.

Because a_i and b_i are constants in the global and local coordinate systems, further rearrangement yields

$$\begin{array}{}\text{(6)}& {\dot{l}}_i{u}_i+l_i{\dot{u}}_i=-\dot{\mathbf{R}}{\mathbf{b}}_i-\dot{\mathbf{p}}.\end{array}$$

Multiply Eq. (6) by $u_i^T$ and simplify

$$\begin{array}{}\text{(7)}& \begin{array}{rl}{\dot{l}}_i=& -{\left[\begin{array}{cccc}{u}_{\mathrm{1}}& u_{\mathrm{2}}& \mathrm{\cdots }& u_m\\ \left(R{b}_{\mathrm{1}}\right)\times u_{\mathrm{1}}& \left(R{b}_{\mathrm{2}}\right)\times u_{\mathrm{2}}& \mathrm{\cdots }& \left(R{b}_m\right)\times u_m\end{array}\right]}_{m\times n}^T\\ & \left[\begin{array}{c}\dot{p}\\ \mathit{\omega }\end{array}\right].\end{array}\end{array}$$

The Jacobi matrix, A^T, and end actuator speed, $\dot{\mathbf{X}}$, are as follows:

$$\begin{array}{}\text{(8)}& {\displaystyle }{\mathbf{A}}^T={\left[\begin{array}{cccc}{u}_{\mathrm{1}}& u_{\mathrm{2}}& \mathrm{\cdots }& u_m\\ \left(R{b}_{\mathrm{1}}\right)\times u_{\mathrm{1}}& \left(R{b}_{\mathrm{2}}\right)\times u_{\mathrm{2}}& \mathrm{\cdots }& \left(R{b}_m\right)\times u_m\end{array}\right]}_{m\times n}^T,\text{(9)}& {\displaystyle }\dot{X}=\left[\begin{array}{c}\dot{p}\\ \mathit{\omega }\end{array}\right].\end{array}$$

The cable speed changes to

$$\begin{array}{}\text{(10)}& {\dot{\mathbf{L}}}_i=-{\mathbf{A}}^T\dot{\mathbf{X}}.\end{array}$$

The change in cable acceleration is

$$\begin{array}{}\text{(11)}& {\ddot{\mathbf{L}}}_i=-{\mathbf{A}}^T\ddot{\mathbf{X}}-{\dot{\mathbf{A}}}^T\dot{\mathbf{X}}.\end{array}$$

3.2 Forward kinematics model

The forward kinematics of HIR based on a cable-driven parallel mechanism calculates the location and pose of the end effector corresponding to a given cable length; the Newton–Raphson iterative algorithm is typically used (Pott and Schmidt, 2015; Lau et al., 2016). Set the initial value, solve the nonlinear equation by iteration, and halt the iteration when the solution encounters an error.

According to Eq. (1), the objective function is constructed as

$$\begin{array}{}\text{(12)}& \mathfrak{R}\left(X\right)={\mathbf{L}}_i-\left({\mathbf{a}}_i-{\mathbf{Rb}}_i-\mathbf{p}\right).\end{array}$$

According to the Newton–Raphson iterative algorithm,

$$\begin{array}{}\text{(13)}& {\displaystyle }{X}_{k+\mathrm{1}}=X_k-{\displaystyle \frac{\mathfrak{R}\left(X\right)}{\dot{\mathfrak{R}}\left(X\right)}},\text{(14)}& {\displaystyle \frac{\mathfrak{R}\left(X\right)}{\dot{\mathfrak{R}}\left(X\right)}}=-{\left(J^{T}J\right)}^{-\mathrm{1}}{J}^T\mathfrak{R}\left(X\right).\end{array}$$

The aforementioned algorithm locally converges, and the selection of an appropriate initial value can ensure convergence. To verify the correctness of the forward and inverse kinematics model, a joint simulation of the forward and inverse kinematics is implemented. Moreover, set the motion track of the end actuator center of mass point, χ₀, in 3D space as

$$\begin{array}{}\text{(15)}& \left\{\begin{array}{l}x=\mathrm{100}\sin \left(t\right)\\ y=\mathrm{100}\cos \left(t\right)\\ z=\mathrm{100}\sin \left(t\right)\cos \left(\mathrm{2}t\right)\end{array}\right..\end{array}$$

The space trajectory of Eq. (15) is shown in Fig. 6. The trajectory motion is relatively complex, ensuring that the robot can complete complex tasks.

Figure 6 Space motion track.

Calculate the length of each cable according to Eq. (2), substitute the cable length value into Eq. (12), and calculate the motion track of the end actuator center of mass according to Eq. (12). The forward and inverse kinematics joint simulation algorithm is described as follows.

Through this joint simulation algorithm, the expected and actual trajectories of different planes are obtained, as shown in Fig. 7a–c; the trajectory error is shown in Fig. 7d.

Figure 7 Projection and error of motion track on each plane.

As shown in Fig. 7a–c, the expected trajectory projected on each plane virtually coincides with the actual trajectory. The order of magnitude of the trajectory error is 10⁻¹³, as shown in Fig. 7d. This proves that the expected trajectory is fundamentally consistent with the actual trajectory. This proves the correctness of the forward and inverse kinematics model of the HIR formulated in this study.

4 Cable tension distribution

The force of the HIR is complex, and the cable tension distribution to obtain the difference in cable tension is critical to the stability of the robot.

4.1 Dynamic model

The sampling Lagrangian modeling method for the dynamic model of HIR (Shang et al., 2020; Babaghasabha et al., 2014; Khosravi and Taghirad, 2014) yields the kinetic energy, E_D, of the end actuator as

$$\begin{array}{}\text{(16)}& E_D={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{v}^T{m}_e{\mathbf{I}}^{\mathrm{3}}v+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{w}^T\left({}_O^{O_{\mathrm{1}}}R\cdot {\mathbf{{\rm Y}}}_{O_{\mathrm{1}}}\cdot {}_O^{O_{\mathrm{1}}}{R}^T\right)w,\end{array}$$

where m_e is the mass of the end actuator; ${\mathbf{I}}^{\mathrm{3}}\in {\mathbf{R}}^{\mathrm{3}\times \mathrm{3}}$ is the identity matrix; v∈R³ and w∈R³ are the translation and angular rotation speeds of the end effector, respectively; and ${\mathbf{{\rm Y}}}_{O_{\mathrm{1}}}\in R^{\mathrm{3}\times \mathrm{3}}$ represents the inertia matrix of the end effector in the local coordinate system.

Similarly, the potential energy, S, of the end actuator is

$$\begin{array}{}\text{(17)}& S=m_{e}gx,\end{array}$$

where g is the gravitation acceleration, and x is the displacement of the end actuator center of mass in the direction of gravity.

The combination of Eqs. (16) and (17) yields the Lagrange equation of the end actuator as

$$\begin{array}{}\text{(18)}& {\displaystyle }\begin{array}{rl}{E}_C=& E_D-S\\ =& {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{v}^T{m}_e{\mathbf{I}}^{\mathrm{3}}v+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{w}^T\left({}_O^{O_{\mathrm{1}}}R\cdot {\mathbf{{\rm Y}}}_{O_{\mathrm{1}}}\cdot {}_O^{O_{\mathrm{1}}}{R}^T\right)w\\ & -{\mathbf{m}}_{e}gx,\end{array}\text{(19)}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left({\displaystyle \frac{\partial E_D}{\partial \dot{\mathbf{X}}}}\right)-{\displaystyle \frac{\partial E_D}{\partial \mathbf{X}}}=\mathbf{AF}+\mathbf{W},\end{array}$$

where $\mathbf{A}\in {\mathbf{R}}^{n\times m}$ is a Jacobian matrix; F∈R^m is the cable tension matrix; and W∈Rⁿ represents the external force, F_h, and moment, M_h.

Substitute Eq. (18) into Eq. (19),

$$\begin{array}{}\text{(20)}& \mathbf{M}\left(\mathbf{X}\right)\ddot{\mathbf{X}}+\mathbf{C}(\mathbf{X},\dot{\mathbf{X}})\dot{\mathbf{X}}+\mathit{G}=\mathbf{AF}+\mathbf{W},\end{array}$$

where M(X) is the inertia matrix, $\mathbf{C}(\mathbf{X},\dot{\mathbf{X}})$ is the centrifugal and Coriolis force matrices, and G is the gravity vector.

Set a group of determined dynamic parameters. Re-express Eq. (20) as (An et al., 1985; Khosla and Kanade, 1985)

$$\begin{array}{}\text{(21)}& \mathbf{AF}=\mathbf{M}\left(\mathbf{X}\right)\ddot{\mathbf{X}}+\mathbf{C}(\mathbf{X},\dot{\mathbf{X}})\dot{\mathbf{X}}+\mathit{G}-\mathbf{W}=\mathbf{N},\end{array}$$

where N∈Rⁿ represents the resultant force and moment of the end actuator.

4.2 Tension distribution

Based on the dynamic model, the tension in each cable is

$$\begin{array}{}\text{(22)}& \mathbf{AF}=\mathbf{N}.\end{array}$$

According to the cable characteristics, the solution of Eq. (22) is (Ouyang and Shang, 2016)

$$\begin{array}{}\text{(23)}& {\mathbf{F}}_{\mathrm{d}}={\mathbf{A}}^{+}\left(\mathbf{N}\right)+\mathbf{Null}\left(\mathbf{A}\right)\mathit{\lambda }.\end{array}$$

In Eq. (23), ${\mathbf{F}}_{\mathrm{s}}={\mathbf{A}}^{+}\left(\mathbf{N}\right)$ is the special solution, F_g=Null(A)λ is the general solution, ${\mathbf{A}}^{+}={\mathbf{A}}^T({\mathbf{AA}}^T{)}^{-\mathrm{1}}$ is the Moore–Penrose inverse of A, $\mathbf{Null}\left(\mathbf{A}\right)\in R^{\mathrm{8}\times \mathrm{2}}$ is the null space vector of A, and $\mathit{\lambda }=[{\mathit{\lambda }}_{\mathrm{1}},{\mathit{\lambda }}_{\mathrm{2}}]$ is an arbitrary vector.

Set the tension condition of the cable according to the physical characteristics of the drive unit and cable, as follows:

$$\begin{array}{}\text{(24)}& {\mathbf{F}}_{max}\ge {\mathbf{F}}_i\ge {\mathbf{F}}_{min}>\mathrm{0},\left\{i=\mathrm{1},\mathrm{2},\mathrm{\dots },\mathrm{8}\right\}.\end{array}$$

In Eq. (14), F_max and F_min are the maximum and minimum tension values in the cable, respectively.

To generate a smooth initial cable tension change and stabilize the system, reduce the system vibrations due to sudden changes in cable tension. The change in cable tension must also satisfy the following condition:

$$\begin{array}{}\text{(25)}& \mathrm{6}\ge {\displaystyle \frac{{\mathbf{F}}_{imax}}{{\mathbf{F}}_{jmin}}}\ge \mathrm{1},\left\{i,j=\mathrm{1},\mathrm{2},\mathrm{\dots },\mathrm{8}\right\},\end{array}$$

where F_imax and F_jmin are the maximum and minimum tension values of the eight cables. The smaller the ratio, the more gradual the change in tension.

Considering the real-time requirements of the HIR, a quadratic programming (QP) solution is implemented to minimize the tension, λ. The optimization algorithm is

$$\begin{array}{}\text{(26)}& \begin{array}{ll}\mathrm{find}& \mathit{\lambda }\\ min& F_{\mathrm{g}}\\ \mathrm{s}.\mathrm{t}.& F_{max}-F_{\mathrm{s}}\ge \mathbf{Null}\left(\mathbf{A}\right)\mathit{\lambda }\ge F_{min}-F_{\mathrm{s}}.\end{array}\end{array}$$

The optimization process in the foregoing algorithm is simple. The feasible solution consists of two-dimensional convex sets of λ with extreme values (Shang et al., 2020). The tension distribution process is shown in Fig. 8.

Figure 8 Tension distribution process.

The space trajectory of the HIR is within the feasible workspace of the force screw. The mass of the end effector is 18.62 kg. The minimum and maximum tensile forces in the cable are 30 and 1000 N, respectively. The position of the hinge point of the end actuator relative to the local coordinate system is listed in Table 1. The tension distribution and change cable in length are shown in Fig. 9a and b, respectively.

Figure 9 Tension distribution and change in cable length.

As shown in Fig. 9, the cable tension is between the minimum and maximum values, and the ratio of the maximum to the minimum tensile forces in the eight cables is greater than 1 but less than 6. The change in the cable length is continuous. The foregoing proves that the tensile force yielded by the proposed tension distribution algorithm is considerably stable. Moreover, the change in cable length is stable and continuous.

5 Control strategy

A double closed-loop scheme is adopted for the design of the driver controller. The inner loop speed suppresses interference, and the outer loop position protection ensures driver position accuracy. The design requires that the speed loop bandwidth is not less than 5 times the bandwidth of the position loop. When the controller is designed with a known transfer function, a stable and reliable classical control theory is adopted. The HIR adopts the same drive unit servo motor, ball screw module, cable, and drive chain structure; thus, only one of the drive unit controllers is designed.

5.1 Mathematical model of drive unit

The drive unit is a direct curve (DC) electric servo system. The rotor of the DC motor generates the driving torque under the action of the armature current. Then, the motor transfers the torque to the cable through the ball screw mechanism. The cable drives the load movement. The drive unit can be simplified as the model shown in Fig. 10 (u_a is the armature voltage, L_a is the armature inductance, R_a is the armature circuit resistance, C_m is the torque constant, e is the back electromotive force (EMF), and C_e is the back EMF constant). The drive unit model is mainly composed of the electric circuit model of the motor and the force balance equation of the ball screw.

Figure 10 Mechanism of drive unit model.

The electric circuit equation of the motor is

$$\begin{array}{}\text{(27)}& L_{\mathrm{a}}{\displaystyle \frac{\mathrm{d}{i}_{\mathrm{a}}}{\mathrm{d}t}}+R_{\mathrm{a}}{i}_{\mathrm{a}}+e=u_{\mathrm{a}}.\end{array}$$

The equation of the back EMF is

$$\begin{array}{}\text{(28)}& e=C_e\cdot {\dot{\mathit{\theta }}}_m.\end{array}$$

The moment equation is

$$\begin{array}{}\text{(29)}& T_m\left(t\right)=C_m\cdot i_{\mathrm{a}}\left(t\right),\end{array}$$

where T_m is the driving torque, and i_a is the armature current.

The force balance equation is

$$\begin{array}{}\text{(30)}& T_m-K_{\mathrm{s}}{F}_{\mathrm{1}}=J_{\mathrm{m}}{\ddot{\mathit{\theta }}}_m+B_{\mathrm{m}}\cdot {\dot{\mathit{\theta }}}_m,\end{array}$$

where J_m is the comprehensive converted moment of inertia of the motor rotor and load; B_m is the comprehensive converted damping of the motor rotor and load; K_s is the transmission ratio of the motor rotor to the sliding block translation; and F₁ is the cable tension.

Equations (28) and (29) are obtained by simultaneous Laplace transformation. The following equation is also transformed:

$$\begin{array}{}\text{(31)}& I_{\mathrm{a}}\left(s\right)={\displaystyle \frac{U_{\mathrm{a}}\left(s\right)}{L_{\mathrm{a}}s+R_{\mathrm{a}}}}-{\displaystyle \frac{C_e{\dot{\mathit{\theta }}}_m\left(s\right)}{L_{\mathrm{a}}s+R_{\mathrm{a}}}}.\end{array}$$

After the Laplace transformation of Eqs. (29), (30), and (31), the mathematical model of the drive unit is obtained as follows:

$$\begin{array}{}\text{(32)}& V\left(s\right)=M_{\mathrm{1}}\left(s\right)U_{\mathrm{a}}\left(s\right)-M_{\mathrm{2}}\left(s\right)F_{\mathrm{1}},\end{array}$$

where M₁(s) and M₂(s) are

$$\begin{array}{}\text{(33)}& {\displaystyle }{M}_{\mathrm{1}}\left(s\right)={\displaystyle \frac{C_m{K}_{\mathrm{s}}}{(J_{\mathrm{m}}s+B_{\mathrm{m}})(L_{\mathrm{a}}s+R_{\mathrm{a}})+C_e{C}_m}},\text{(34)}& {\displaystyle }{M}_{\mathrm{2}}\left(s\right)={\displaystyle \frac{-K_{\mathrm{s}}^{\mathrm{2}}(L_{\mathrm{a}}s+R_{\mathrm{a}})}{(J_{\mathrm{m}}s+B_{\mathrm{m}})(L_{\mathrm{a}}s+R_{\mathrm{a}})+C_e{C}_m}}.\end{array}$$

The model block diagram of the drive unit obtained from Eqs. (31) and (32) is shown in Fig. 11.

Figure 11 Block diagram of drive unit model.

As shown in Fig. 11, the output speed of the drive unit is determined by the armature voltage and slider load torque. When the load force is fixed, the slider speed can be controlled by changing the armature voltage of the motor. When the armature voltage of the motor is fixed, the change in the load force can cause the motor speed to change. Therefore, to eliminate the fluctuation of the slider speed caused by the change in the load torque, the slider speed must be accurately controlled through the closed-loop control of the motor speed feedback.

Figure 12 Block diagram of improved drive unit model.

The controlled drive unit must also include the mechanical and control parts. By ignoring the influence of drive unit stiffness and considering the electric circuit equation (Eq. 27) of the motor, the drive unit model improves, as shown in Fig. 12. The drive unit limits current by limiting voltage. Resolving the problem of limiting the current in the control strategy also protects the drive unit because the voltage and current of the unit are regulated at the same time. Conventional hardware current chopping, current inner loop, and other methods typically allow current overload and adopt shutdown to protect equipment. The proposed drive unit model avoids the foregoing situations and protects equipment without halting equipment operation.

Table 2 Base parameters of drive unit.

The drive unit parameters are summarized in Table 2. Equipment parameters are provided by component manufacturers and measured by experiments. All controller design and simulation experiments adopt the listed parameters.

5.2 Controller design

5.2.1 Speed controller

The motor speed is fed back by the drive unit encoder. The control block diagram of the closed-loop system of the motor speed is shown in Fig. 13.

Figure 13 Block diagram of speed control system.

Figure 14 Zero-pole distribution of transfer function and Bode diagram.

Without considering the change in the cable tension of the drive unit, the driving function of the drive unit is expressed as

$$\begin{array}{}\text{(35)}& G_v\left(s\right)={\displaystyle \frac{V\left(s\right)}{U_{\mathrm{a}}\left(s\right)}}={\displaystyle \frac{C_m{K}_{\mathrm{s}}}{(J_{\mathrm{m}}s+B_{\mathrm{m}})(L_{\mathrm{a}}s+R_{\mathrm{a}})+C_e{C}_m}}.\end{array}$$

Substitute the drive unit parameters listed in Table 2 into Eq. (35). This derives the transfer function of the drive unit:

$$\begin{array}{}\text{(36)}& G_v\left(s\right)={\displaystyle \frac{\mathrm{59.51}}{s^{\mathrm{2}}+\mathrm{100.6}s+\mathrm{4004}}}.\end{array}$$

According to Eq. (36), the two poles of the transfer function are $p_{\mathrm{1}}=-\mathrm{50.3}+\mathrm{38.3915}i$ and $p_{\mathrm{2}}=-\mathrm{50.3}-\mathrm{38.3915}i$. The poles are on the left-half plane of the S domain, and the system is in the minimum phase. The system's natural frequency, damping ratio, and open-loop gain are w_n=63.2772 rad s⁻¹, ξ_n=0.7949, and K=59.51, respectively, as shown in Fig. 14a. When the drive unit is within the operating frequency range, the amplitude–frequency characteristic of the system's forward channel transfer function is stable. Moreover, the phase frequency characteristic is approximately linearly transformed, and the correction system forward channel is relatively simple, as shown in Fig. 14b.

Figure 15 Open-loop Bode diagram and sine tracking curve of system after speed loop correction.

Because the forward channel of the drive unit is a type-0 system, the loading of the high-order static error-free system cannot be satisfied, and the hysteretic characteristics lead to the inadequate dynamic performance of the system. The proportional–integral–derivative (PID) controller is used to correct the forward channel, improve the steady-state accuracy of the drive unit, and eliminate the steady-state error. The transfer function, G_c2, of the PID controller is

$$\begin{array}{}\text{(37)}& G_{\mathrm{c}\mathrm{2}}\left(s\right)=k_{pi}\left(\mathrm{1}+{\displaystyle \frac{k_{vi}}{s}}\right).\end{array}$$

The designed controller parameters are k_pi=0.0047 and k_vi=33.5319. The system's open-loop Bode diagram after calibrating the speed loop controller is shown in Fig. 15a. The shear frequency and phase margin are w_c=57.6 rad s⁻¹ and 66.6^∘, respectively. The amplitude–frequency characteristic crosses the 0 dB line with a slope of −20 dB per decade. After correcting the speed loop, the system stability and rapidity are ensured, and the steady-state error of the system is eliminated. The situation of the system when the sine signal with an amplitude of 1000 is inputted is shown in Fig. 15b. The system is stable with a maximum error of 5^∘. The system phase lag time is 0.0081 s.

5.2.2 Double closed-loop controller

The block diagram of the position and speed double closed-loop control system is shown in Fig. 16.

Figure 16 Control block diagram of position and speed double closed-loop control system.

Figure 17 Proposed controller compared with traditional PID Bode diagram and sinusoidal tracking curve.

Figure 18 Control strategy of HIR system.

Figure 19 Subsequent condition of each cable in space track completed by HIR.

When the PID controller is started or the amplitude of the input signal significantly varies, the system output deviates considerably. Moreover, the integral of the controller operation accumulates, thus overshooting the system and oscillation. To stabilize the system and improve its dynamic tracking control accuracy, the position loop adopts a variable speed integral separation PID controller. This PID controller can adjust integration and modify integration speed according to the deviation of the system. A pure differential link of the actual system is physically impossible, and high-frequency noise easily infiltrates. Given this, the differential link is connected in series using a low-pass filter with a cutoff frequency of N rad s⁻¹. Coefficient f(e), which is related to deviation change, is introduced to modify the integral strength. The variable speed integral separation PID transfer function is

$$\begin{array}{}\text{(38)}& G_{\mathrm{c}\mathrm{1}}\left(s\right)=k_{pj}\left(\mathrm{1}+{\displaystyle \frac{k_{vj}f\left(e\right)}{s}}+{\displaystyle \frac{k_{dj}s}{\frac{s}{N}+\mathrm{1}}}\right).\end{array}$$

Set the relationship between deviation e and f(e) as

$$\begin{array}{}\text{(39)}& f\left(e\right)=\left\{\begin{array}{ll}\mathrm{1}& \left|e\right|\le {\mathit{\epsilon }}_{\mathrm{1}}\\ \frac{{\mathit{\epsilon }}_{\mathrm{2}}-\left|e\right|+{\mathit{\epsilon }}_{\mathrm{1}}}{{\mathit{\epsilon }}_{\mathrm{2}}}& {\mathit{\epsilon }}_{\mathrm{1}}<\left|e\right|\le {\mathit{\epsilon }}_{\mathrm{1}}+{\mathit{\epsilon }}_{\mathrm{2}}\\ \mathrm{0}& \left|e\right|>{\mathit{\epsilon }}_{\mathrm{1}}+{\mathit{\epsilon }}_{\mathrm{2}}\end{array}\right..\end{array}$$

For this system, ε₁=0.2 mm, ε₂=0.8 mm, k_pj = 72 000, k_vj=100, k_dj=0.12, and N=200 are selected for simulation debugging. When the influence of nonlinear factors is ignored, the open-loop Bode diagram of the system after correcting the proposed position and PID controllers is shown in Fig. 17a. The turning frequency of the proposed position controller exceeds that of the PID controller and is considerably stable. The input is the typical sine signal, S=100sin (t). The resulting situation of the system is shown in Fig. 17b. The performance of the proposed control system is better than that of the PID controller, and the system's dynamic performance is satisfactory. When t=1.71 s, system integration is complete.

5.3 Control strategy of robot system

The control strategy of the HIR system is shown in Fig. 18.

Figure 20 Projection and error of space track on each plane.

To verify the correctness of the proposed control strategy for the HIR system, the strategy is used to complete the space trajectory of Eq. (15). The combination of the force exerted by the operator and the force generated by the virtual environment is given by $F_h=F_{h\mathrm{1}}+F_e$. The results are shown in Fig. 19.

The dotted line in Fig. 19 represents the theoretical value of the corresponding cable, generated by the admittance model, and the solid line represents the drive unit following the theoretical value. In 6 s, the drive unit fully enters the following state with satisfactory tracking performance; tracking is stable for 100 s.

To demonstrate the stability of the system control strategy further, a 10 N force is applied to the end actuator, causing it to move along a straight line; the simulation results are shown in Fig. 20. The dotted line is the theoretical value of the corresponding cable generated by the admittance model, and the solid line is the theoretical value of the following drive unit. In 5 s, the drive unit completely enters the following state with satisfactory tracking performance; tracking is stable for 100 s. Therefore, the proposed drive unit and robot system control strategies can stabilize the system.

6 Conclusion

In this study, a new drive unit configuration and an end effector configuration are proposed for HIRs. The drive unit configuration can improve the motion accuracy of the cable, and the end effector configuration can render the cable motion flexible. The algorithm for distributing cable tension is improved to render the change in cable tension gradual. The mathematical model, control strategy, over-current protection, and motion control accuracy of the drive unit are enhanced. Moreover, a control strategy for stabilizing the control of the HIR system is proposed.

In future studies, the HIR will be optimized, considering human–robot interaction safety and other conditions to reduce volume and build a compact robot. The control strategy of the robot system will be improved, such that the drive unit can rapidly and accurately follow the theoretical trajectory.