Problem formulation

The motion control problem of robot manipulators carries several issues related to the effectiveness of the implemented controllers and the manipulator’s capability of executing some specific tasks. For instance, the efficacy of model-based controllers directly depends on the accuracy of the system model. However, in some cases, there may be discrepancies between such model and the physical system, mainly due to parameter uncertainty and unmodeled dynamics; these discrepancies may diminish the efficacy of these controllers. Furthermore, robot manipulators that operate away from controlled environments, such as the ones attached to explorer vehicles for space and deep-water missions, are usually submitted to unknown external forces, which affects the performance of the tasks assigned to the robot manipulator.

The described problems above can be treated as one if we group parameter uncertainty, unmodeled dynamics, and unknown external forces as a single perturbation term. This problem is known as disturbance rejection and can be addressed using a robust control strategy.

Therefore, the robust control strategy has relevance because it allows decreasing the effects of the mentioned features. Moreover, with this strategy, it is possible to obtain good performance for the robot manipulator tracking control problem. This work presents a solution for this issue through a sliding-mode-based control law known as ST controller. The methodology of the present work illustrated in Figure 1.

Mathematical model

Let be the 4-joint robot manipulator shown in Figure 2a (without end-effector). Consider the coordinate frames fixed to the manipulator (see Figure 2b). The dynamic model is obtained using Lagrange’s equations, where the Lagrangian L(q,q˙) of a robot manipulator of 4 DOF is the difference between its kinetic energy K and its potential energy, U, That is:

Then the kinetic energy function K(q,q˙) and potential energy U(q) may be decomposed into the sum of four parts: Kq,q˙=K1q,q˙+K2q,q˙+K3q,q˙+K4q,q˙ and Uq=U1q+U2q+U3q+U4(q). Hence, using these equations, the Lagrangian is:

In addition, Rayleigh’s dissipation function is introduced into Euler-Lagrange’s dynamic equations, obtaining:

Where: RF=Bq˙2/2. Equation (3) corresponds to the manipulator dynamic equations, including the friction forces.

Thus, four dynamic equations correspond to the torques delivered by the actuators and can be expressed as follows:

Where: q∈Rnx1 is the position vector, q˙∈Rnx1 corresponds to the velocity vector, B∈Rnx1 is the vector that contains the viscous friction coefficients, q¨∈Rnx1 indicates the acceleration vector, Mq∈Rnxn is the symmetric positive definite manipulator inertia matrix (nonsingular and invertible), Cq,q˙∈Rnxn describes the centrifugal and Coriolis forces, G(q)∈Rnx1 is the vector of gravitational forces and τ∈Rnx1 represents the torque vector. The complete dynamic model is shown in Appendix A (^{Ibarra et al., 2019}).

Manipulator parameterization

The control algorithm design requires the parameterization of the manipulator physical parameters used in the dynamic model. A manipulator design is made in the CAE-CAD SolidWorks software with the actual values of the structure dimensions and considering aluminum alloy 1060 as the primary construction material.

By using the virtual model of the manipulator, it is possible to analyze the physical properties of each of the links that make up the manipulator and thus obtain the parameters that describe it.

Considering that “i” is the ith link in the part of the manipulator structure that connects the joint (i+1) with the ith joint and separating the virtual model of the manipulator into four assemblies corresponding to the four links, each link of the manipulator is analyzed in SolidWorks as an independent assembly.

It is necessary to assign a new reference frame to each link based on the rules of Denavit Hartenberg’s method. The physical properties obtained are the total mass of each link, the entire length of each link, the center of mass of each link concerning the new reference frame, and the inertia tensor of each link to the new reference frame.

The viscous friction parameters which are present in the servo motors coupled to each joint are considered. Such friction values are obtained from a previous system identification made for the engines used in a manipulator with similar characteristics to the one studied in this work (^{Corke, 2017}). The results of the parameterization performed are grouped in Table 1.

The notation used in Table 1 is, centers of mass correspond to lci=x, y, z, which denotes the coordinate of the center of mass for the reference frame of each link. Also, the notation for the inertia tensor indicates a diagonal matrix, whose elements on the main diagonal correspond to the main moments of inertia Ixx, Iyy and Izz.

Control strategies

a. PD+ control

The PD+ Control for trajectory tracking is a robust passive computed torque controller (^{Han et al., 2020}), which requires the exact knowledge of the manipulator’s dynamic model to be implemented (^{Kelly et al., 2006}; ^{Reyes, 2011}), that is, of Mq, C(q,q˙), G(q) and B. In addition, it is necessary to know the desired trajectories q*t, q˙*(t) and q¨*(t) as well as to have the measurements q(t) and q˙(t). The PD+ is a globally stable controller that results from a Lyapunov function candidate Vt=1/2(q~˙TMqq~˙+q~Tkpq~), obtaining the following equation (see ^{Slotine & Weiping (1988)} for more information):

Where: kp and kv are symmetric positive definite matrices, q~=q*t-q(t) denotes the position error and q~˙=q˙*t-q˙(t) the velocity error in the manipulator’s joints. For the application of the control law described in (5), it is necessary that Mq and Cq,q˙ are not independent, the matrix (M˙q-2Cq,q˙) must be skew-symmetric and that the dynamic structure is linear in terms of a suitably selected set of robot and load parameters (^{Slotine & Weiping, 1988}; ^{Han et al., 2020}). Figure 3 shows the block diagram of the PD+ Control. For this case, we consider that we have access to all system information, compensating Mq,  Cq,q˙, Gq, and B.

The gain tuning of the PD+ controller does not depend on the manipulator model nor the task assigned to it. It is usually suggested to be performed manually and based on the designer’s knowledge (^{Kelly et al., 2006}; ^{Reyes, 2011}). To propose a tuning procedure for the PD+ controller, we assume that the control law compensates the nonlinear nature of the manipulator at some instant of time, then we can describe the manipulator dynamics of the form:

The linear system described in (6) is used for tuning purposes only and can be represented by the following block diagram:

The closed-loop equation is:

The selection of kp and kv must guarantee that the characteristic polynomial is Hurwitz. Therefore, we use the following reference model:

Thus, kpi=ωni2 and kvi=2ζiωni,   for i=1,2,3,4. To have a low absorption energy behavior, the damping coefficients ζi=0.7 and undamped natural frequencies ωni=5.2, 4.8, 4, 3.5 are proposed. Then we group the gains in diagonal matrices as follows:

The controller gains, kp and kv, are positive definite symmetric matrices, this being a sufficient condition of stability (^{Kelly et al. 2006}). Finally, Figure 5 shows the block diagram of the PD+ Control in closed-loop for path tracking.

b. Control Super-Twisting control

Consider the existence of a Lipschitz disturbance term, the dynamic model of the robot manipulator is rewritten as:

Where: F(q,q˙) represents the perturbations/uncertainties vector produced by the disturbance force affecting the system, parameter uncertainties, and unmodeled dynamics. This term is assumed Lipschitz. To proceed with the design of an SMC law, we define the sliding surface as follows:

The purpose of SMC is to force the sliding surface dynamics to reach the value of zero in a finite time and keep it at zero thereafter, i.e.:

Once this is accomplished, the controlled dynamics are described by:

The solution of (14) is:

Thus, it is shown that the position error and its derivative converge to zero asymptotically. Deriving the sliding surface, we obtain:

Where it is clear that: q~¨=q¨*-q¨, so (16) can be rewritten as:

Then, substituting (11) into (17), we have that:

We propose the following control law:

Where: τeq the equivalent nominal control, responsible for compensating the known system dynamics. Thus it is defined as follows:

Thus, ueq must satisfy the condition σ˙t=0, substituting (20) in (18), we deduce:

Equation (21) is the equivalent control and indicates that the second term of the control law should be equal to the negative of the disturbance term to assure the system’s robustness. However, this control action usually cannot be implemented, since it depends on the disturbance F(q,q˙) which is usually unknown (^{Shtessel et al., 2014}).

To satisfy condition (13), the use of an SMC with the equivalent nominal Control is proposed, providing robustness to the system in the presence of disturbances. The proposed controller is a Super-Twisting (ST) control. Then the control law is rewritten as:

Where: τst is the ST control action defined as follows:

With k1, k2>0 and σtn12=σt112signσt1, …,  σtn12signσtnT.

Figure 6 depicts the block diagram of the ST controller for the robot manipulator.

Substituting (22) in (18), we obtain:

If Ω=ω+M(q)-1F(q,q˙), then (25) becomes:

Thus, Ω˙=-k2signσt+Mq-1F˙q,q˙+ddtMq-1F(q,q˙), proposing φ=Mq-1F˙q,q˙+ddtMq-1Fq,q˙ and φ≤φMax=L. In this way, we obtain the following sliding surface:

At this point, the design problem is reduced to the setting of the parameters k1 and k2. In this work, we use the tuning rules for an ST controller proposed in (^{Levant, 1998}), which indicates that:

The stability proof for this tuning procedure is not developed in this paper. In exchange, we present sufficient conditions for the controller to be stable in finite time. In the interest of doing so, we use the stability proof for Super-Twisting parameter setting shown in ^{Seeber & Horn (2017)}, from where we extract the following. Let us define the below system:

Where: δt is an absolutely continuous perturbation bounded by δt≤L.

Theorem 1: The system (30)-(31) is stable in finite time if its parameters satisfy the following conditions (^{Seeber & Horn, 2017}).

Comparing (27) and (28) with (30) and (31), we can see that the sliding surface corresponds to the system form referred to in Theorem 1. Thus, such a theorem is valid to analyze the stability of the sliding surface. Observe that the parameter setting (29) proposed in (^{Levant, 1998}) satisfy condition (32), which proves the finite-time stability of the ST controller with the use of these gains. Therefore, the sliding surface designed in this work is stable in a finite time by the following control law:

To proceed with the setting of the gains k1 and k2, we estimate the bound L via the use of a Matlab-Simulink simulation environment. We propose an oscillatory perturbation with maximum perturbation torque of 1.5 Nm. The computed values were:

Finally, the design parameter λ was proposed to be λ=2.9.

Remark: Note that the proposed sliding mode controller is capable of compensating asymptotically Lipschitz uncertainties/perturbations. In the case of a discontinuity in the perturbation vector F, the controller will lose convergence, but after the discontinuity, it will reconverge if the derivative bound condition is fulfilled.

Inverse Kinematics

We use inverse kinematics to link the robot manipulator’s position and orientation in the workspace with the robot manipulator joint positions. In this paper, we decide to employ a geometric approach to obtain the inverse kinematics of the robot manipulator. Furthermore, we choose to restrict the inverse kinematic solution to the elbow-down solution.

Consider Figure 7, which represents the robot manipulator links and where the θi angles correspond to the angles formed by the robot manipulator links on each of the joints, with i=1,2,3,4. In addition, the angles formed by the second and third links are depicted in Figure 8. Furthermore, the F reference frame corresponds to the reference frame attached to the robot manipulator end effector.

The following homogenous transformation matrix denotes the position and orientation of the end effector:

Where: R∈R3x3 denotes the orientation matrix of the end effector and O∈R3x1 represents the origin point of the F reference frame.

From Figure 7 and equation (36), we define the position vector Oc as follows:

Then, θ1 is derived from Figure 7 and equation (37) as follows:

From Figure 7 and Figure 8, we define the following parameters:

Where: c^=180-θ3, and by trigonometric identities, we get:

Then, substituting equation (42) in equation (41), we obtain:

From which we derived the following:

Then, θ3 is computed as follows:

Similarly, from Figure 8, θ2 is computed as follows:

Similarly, to compute θ4 we define the following:

Then, equation (51) is rewritten as:

Finally, θ4 is derived as follows:

Therefore, the solution to the inverse kinematics is defined by equations (38), (46), (49), and (54), which are grouped in a vector that denotes the actual joint positions of the robot manipulator given a homogenous transformation matrix for the robot manipulator end effector:

Robot task

The trajectory proposed in this work is intended to be smooth in time and be reachable by the end effector within the workspace of the robot manipulator. Such trajectory in the workspace is described by the following equations:

And finally, the variable zt is defined by a Bezier curve characterized by the following parameters: t1=0, t2=19, z0=0.066 and zf=0.066. Thus, these three variables are grouped in the following position vector:

Which denotes the desired position of the robot manipulator end effector, and where the time variable t goes from 0 to 19 seconds. The total displacement of the trajectory in the workspace is depicted in Figure 9.

The orientation of the robot end-effector can be defined in terms of the zf-axis of the robot end effector (see Figure 7). We define a variable named Yaw to denote the desired rotation angle for the robot end-effector. However, the orientation of the robot end-effector is constant at all times. Thus, Yaw is fixed as -0.3 rad, and the orientation of the robot end-effector is defined by the following rotation matrix:

Finally, we can group both, the position vector, and the rotation matrix in the following homogenous transformation matrix:

Thus, for each value of Hdt in time, a value for the vector qdt=qd1,qd2,qd3,qd4T is calculated using the inverse kinematics equations. Where qdt denotes the desired robot manipulator joint positions.

In summary, the robot’s end-effector must follow the trajectory and orientation described in (58). The inverse kinematics equations are thus used to find the necessary angles in each joint. Therefore, the reference values to be entered into the controller are obtained. Figure 10 shows the block diagram to get the desired angles for the robot’s joints.

Numerical simulations

The designed control algorithms were tested in Matlab-Simulink, using the robot manipulator’s dynamic model to evaluate the performance of each controller for a trajectory tracking task. Such desired trajectory was designed to be smooth and differentiable in time.

As it was pointed out before, the perturbation term must be Lipschitz. Taking this into account, we propose Fq,q˙=-0.8cos⁡2t+0.5q+0.2q˙Nm. The simulation fixed step time was set at 0.001s.

The performance of each control algorithm in both joint space, and workspace, for the trajectory tracking task in the presence of perturbations are shown in Figure 11 and Figure 12. From these figures, the PD+ Control is evidently more affected by the disturbance force than the Super-Twisting Control.

In addition, in Figure 13, we present the required torque by each controller. Note that the demanded torque by the Super-Twisting Control is a continuous function (see Figure 13b). Furthermore, it is essential to point out that part of the tunning criteria considers that the control outputs were bounded within a range corresponding to the actual value of torque given by commercial actuators.