1 Introduction

A typical nonlinear identification problem consists of selecting an appropriate identification model and adjusting its parameters according to some adaptive law, such that the response of the identification model to input signals approximates the response of the real system to the same inputs. The identification of nonlinear systems can be posed as a nonlinear functional approximation problem.

Since neural networks (NNs) have good approximation capabilities, even with one hidden layer, and inherent adaptivity features, they provide a powerful tool for identification of nonlinear systems [1, 2]. In [3,4,5] the approximation and learning properties of one class of recurrent neural networks known as recurrent high–order NNs (RHONNs) have been studied. This architecture was applied for identification of a two degrees of freedom (DOF) mathematical planar robot manipulator model, on which gravitational forces affecting the robot dynamics and terms related to friction were not considered. Stability and convergence properties of RHONNs as models of nonlinear dynamical systems were also studied, showing that if a sufficiently large number of high–order connections between neurons are allowed then the RHONN model is capable of approximating a large class of dynamical systems. Using Lyapunov stability arguments, it was shown for these algorithms that the identification error converges exponentially fast to zero.

Neural control systems are defined as control systems, in which at least one NN is directly involved in generating the control command. In [6], the most popular types of neural control systems are briefly introduced and their main features are reviewed. In the literature, few works on which continuous–time HONN structures for identification and control of robot manipulators have been reported. A continuous–time decentralized control strategy using a recursive method and based on a recurrent NN (RNN) identifier with a block control structure for both identification and control was proposed in [7]. The RNN learning law was performed using the filtered error (FE) technique [8]. This approach was tested via simulation on a two DOF robot manipulator model, on which terms due to friction were not considered. An adaptive tracking controller algorithm in continuous-time for nonlinear multiple input-multiple output (MIMO) systems was presented in [9]. This algorithm includes the combination of a RHONN with high–order sliding modes techniques. The RHONN is used to identify the dynamics for the plant where the FE algorithm is used to train the neural identifier. A decentralized high-order sliding modes in twisting modality was used to design reduced-chattering independent controllers to develop through simulation trajectory tracking for a three DOF robot arm. It has recently been recognized that the accuracy of the conventional control approaches in high-speed applications is greatly affected by the presence of unmodeled dynamics caused by friction forces, vibration, backlash and uncertainties in the parameters describing the dynamic properties of the grasp load (e.g., unknown moments of inertia). To compensate the existing uncertainties and unmodeled dynamics, many research works have been proposed adaptive and learning strategies for controlling rigid robot manipulators that ensure asymptotic trajectory tracking. The development of effective adaptive controllers represents an important step toward high speed-precision robotic applications. In [10], results about trajectory tracking control for a two DOF Rotational-Rotational planar robot manipulator using feedback linearization and an adaptive interval type-2 fuzzy neural compensator were presented, showing that the joint positions were controlled under wide variation of operation conditions and existing uncertainties. In [11], a NN control was presented for a rehabilitation robot with unknown system dynamics. To deal with the system uncertainties and improve the system robustness, adaptive NNs were used to approximate the unknown model of the robot and to adapt interactions between the robot and the patient. Both full state feedback control and output feedback control were considered. Uniform Ultimate Boundedness (UUB) of the closed–loop system was achieved in the context of Lyapunov stability theory and its associated techniques. It was proven that the state of the system converges to a small neighborhood of zero by an appropriate choice of the design parameters. In [12], an adaptive neural tracking control was proposed for a robotic manipulator with uncertainties in both manipulator dynamics and joint actuator dynamics. The manipulator joints were subject to inequality constraints, i.e., the joint angles are required to remain in some compact sets. Integral Barrier Lyapunov Functionals were employed to address the joint space constraints directly without performing an additional mapping to the error space. NNs were used to compensate for the unknown robot dynamics and external force. Adapting parameters were developed to estimate the unknown bounds on NN approximations. By the Lyapunov synthesis, Semi-Global Uniform Ultimate Boundedness of the closed-loop system was guaranteed, and the practical tracking of joint reference trajectory was achieved without the violation of predefined joint space constraints. In [13], the authors focused on adaptive neural control of a robot manipulator with unknown system dynamics under the limitation of prescribed performance. A performance function was introduced to express the prescribed constraints of tracking errors. Subsequently, a performance transformation method was proposed in order to solve the problem of the prescribed performance. The unknown dynamics of the robot were approximated by radial basis function NNs (RBFNNs). An adaptive neural control law was proposed based on Lyapunov stability. The proposed control method guaranteed that all the signals in the closed-loop system were UUB and simultaneously the tracking errors converge to small preassigned regions of zero with satisfying prescribed performance. In [14], a RBFNN based adaptive control was designed for nonlinear robot manipulators. Barrier Lyapunov Function (BLF) technique and terminal sliding mode (TSM) technique were seamlessly integrated to achieve both tracking and NN learning convergence for a finite time. In [15], the authors have studied the tracking control problem for an uncertain m-link robot with full-state constraints. The rigid robotic manipulator was described as a MIMO system. An adaptive NN control for the robotic system with full-state constraints was designed. In the control design, the adaptive NNs are adopted to handle system uncertainties and disturbances. The Moore–Penrose inverse term was employed in order to prevent the violation of the full-state constraints. A BLF was used to guarantee UUB of the closed-loop system. The control performance of the closed-loop system was guaranteed by an appropriate choice of the design parameters. In general, the robotic behavior was divided into free and constrained motions. When a robot moves in space, it is unconstrained and its end-effector is free of external forces. In this case, compensation of friction in joints and transmissions of robotic actuators is important to achieve accurate tracking performance of the robot manipulator. The joint friction compensation of robots has been reviewed in [16]. Likewise [17] deals with a robust non-smooth non-linearity compensation scheme for a direct-drive robot manipulator with an asymmetric, deadzone, uncertainty, dynamic friction in joints, and between environmental contact space and end-effector. A model-free recurrent fuzzy NN control system to approximate the ideal backstepping control law was designed to replace the traditional model-based adaptive controller. Other applications of neural networks in robotics can be reviewed in [18,19,20].

In this work, we have been implementing in real-time the proposal of a RHONN structure in a block strict-feedback form in order to identify the dynamic behavior of a two DOF robotic arm of unknown parameters, of our own design and powered by industrial servomotors, evolving in the vertical plane on which friction is inherently present. The training of the RHONN is performed online, in a series-parallel configuration using the FE technique. The stabilizing controller is derived from the block strict-feedback form proposed for the RHONN using the integrator backstepping procedure. The effectiveness of the proposed neural control scheme is verified via experimental results. As it is well known, nonlinear physical systems are continuous in nature. The transformation of a continuous-time nonlinear system to a discrete-time equivalent model is not straightforward and often does not maintain certain useful properties of the continuous-time model, e.g., affine in the control [21]. Since digital control systems may be treated as continuous-time systems in analysis and design if high sampling rates are used [22], hence, the nonlinear system analysis and neural control scheme design are performed in continuous–time form.

This paper is organized as follows: In Sect. 1 an introduction about state of art in this topic is presented. The RHONN model description, associated theorems and assumptions are given in Sect. 2; Sect. 3 describes the design of the neural control scheme proposed via integrator backstepping approach; Sect. 4 gives a description of the physical plant and real-time results are exhibited in order to evaluate the effectiveness of the neural controller scheme proposed; and finally, conclusions are given in Sect. 5.

2 Recurrent High–Order Neural Networks Model

Let consider a RHONN consisting of i neurons and j inputs, where the state space variable of each neuron is governed by a differential equation of the form [8]

$$\begin{aligned} {\dot{x}}^{i}_{j}=-a^{i}_{j} x^{i}_{j}+b^{i}_{j}\sum _{k=1}^L w^{i}_{jk}\prod _{j \in I_k} {y_{j}}^{d_j(k)}, \end{aligned}$$
(1)

where \(\{I_1, I_2,\ldots , I_L\}\) is a collection of L non-ordered subsets of \(\{1, 2, \ldots ,j+i\}\), \(d_j(k)\) are non-negative integers, and

$$\begin{aligned} {{{\varvec{y}}}}=\begin{pmatrix} y_1 \\ \vdots \\ y_{i} \\ y_{i+1} \\ \vdots \\ y_{i+j} \\ \end{pmatrix} = \begin{pmatrix} s(x^1) \\ \vdots \\ s(x^{i}) \\ u^1 \\ \vdots \\ u^{j} \\ \end{pmatrix} \end{aligned}$$
(2)

is the input vector to each neuron, where \({{\varvec{u}}}=[u^1 \; u^2 \; \cdots \; u^{j}]^\top \) is the external control input vector to the network and the superscript “\(\top \)” denotes transpose vector. The function \(s(\cdot )\) is a differentiable sigmoidal function, monotonic increasing, of the form:

$$\begin{aligned} s(x)=\frac{\alpha }{1+e^{-\beta x}}-\varsigma , \end{aligned}$$
(3)

where the parameters \(\alpha \) and \(\beta \) represent the bound and slope of the curvature of the sigmoid function and \(\varsigma \) is a bias constant. In order to avoid singularities, from the above equation, the term \({e^{-\beta x}}\) must be different from \(-1\) or the term \({1+e^{-\beta x}}\) different from 0.

Introducing a L-dimensional vector defined as

$$\begin{aligned} {{{\varvec{z}}}}=\begin{pmatrix} z_1 \\ z_2 \\ \vdots \\ {z_L} \\ \end{pmatrix} = \begin{pmatrix} {\prod }_{j \in I_1} y_j^{d_j(1)}\\ {\prod }_{j \in I_2} y_j^{d_j(2)}\\ \vdots \\ {\prod }_{j \in I_L} y_j^{d_j(L)}\end{pmatrix}, \end{aligned}$$
(4)

the RHONN model (1) can be rewritten as

$$\begin{aligned} {\dot{x}}^{i}_{j}=-a^{i}_{j}x^{i}_{j}+b^{i}_{j} \sum _{k=1}^L w^{i}_{j k} z_{k}. \end{aligned}$$
(5)

The adjustable parameter vector is defined by \({{\varvec{w}}}^{i}_{j}=b^{i}_{j}[w^{i}_{j1} \; w^{i}_{j2} \; \cdots \; w^{i}_{jL}]^\top \), so (5) becomes

$$\begin{aligned} {\dot{x}}^{i}_{j} = -a^{i}_{j} x^{i}_{j} + ({{{{\varvec{w}}}}^{i}_{j}})^\top {{{\varvec{z}}}}, \end{aligned}$$
(6)

where the weights \(w_{jk}^{i}\), for \(k=1,2,\ldots ,L\), represent the adjustable weights of the network, while the coefficients \(a_{j}^{i}\) for \(i=1,2,\ldots ,n\) are part of the underlying network architecture and are fixed during training. Expressing Eq. (6) in vector notation, the dynamic behavior for the overall network is described by

$$\begin{aligned} \dot{{{\varvec{x}}}}_{j}^{i} = {{\varvec{A}}}{{\varvec{x}}}_{j}^{i} + {{\varvec{W}}}^\top {{\varvec{z}}}_{j}^{i}, \end{aligned}$$
(7)

where \({{\varvec{x}}}_{j}^{i} = [x_{j}^1 \; x_{j}^2 \; \cdots \; x_{j}^n]^\top \in {\mathfrak {R}}^{n},\ {{\varvec{W}}}=[{{\varvec{w}}}_{j}^{1} \; {{\varvec{w}}}_{j}^{2} \; \cdots \; {{\varvec{w}}}_{j}^{n}]^\top \in {\mathfrak {R}}^{L{\times }n}\), and \({{\varvec{A}}} = \text {diag} \{-a_{j}^{1}, -a_{j}^{2}, \ldots , -a_{j}^{n}\} \in \mathfrak {R}^{n \times n}\) is a diagonal matrix. Since each \(a_{j}^i\) is positive, \({{\varvec{A}}}\) is a stability matrix. It must be noticed that the vector \({{\varvec{z}}}_{j}^{i}\) is a function of both the neural network state \({{\varvec{x}}}_{j}^{i}\) and the external input \({{\varvec{u}}}\).

2.1 Approximation Properties of the RHONN

In the following, the problem of approximating a general nonlinear dynamical system by a RHONN is described. The input–output behavior of the system to be approximated is given by

$$\begin{aligned} \dot{\varvec{\chi }}_{j}^{i}={{\varvec{F}}}(\varvec{\chi }_{j}^{i},{{\varvec{u}}}^{i}), \end{aligned}$$
(8)

where \({{\varvec{u}}}\in {\mathfrak {R}}^{i}\) is the input control vector to the system, \(\varvec{\chi }_{j} \in \mathfrak {R}^j\) is the state vector of the system and \({{\varvec{F}}}: {\mathfrak {R}}^{i+j}\mapsto {\mathfrak {R}}^{j}\) is a smooth vector field defined on a compact set \({\mathcal {Y}}\subset {\mathfrak {R}}^{i+j}\), where i and j are constants.

The approximation problem consists of determining, by allowing enough high–order connections, if there exists a matrix of synaptic weights \({{\varvec{W}}}\) such that system (7) approximates the input–output behavior of an arbitrary dynamical system of the form (8). Assume that \({{\varvec{F}}}(\cdot )\) is continuous and satisfies a local Lipschitz condition such that (8) has a unique solution and \((\varvec{\chi }_{j}^{i}(t),{{\varvec{u}}}^{i}(t)) \in {{\mathcal {Y}}}\) for all t in some time interval \(J_T = \{ t: 0 \le t \le T \}\) where \(J_T\) represents the time period over which the approximation is performed.

Theorem 1

[4] Suppose that system (8) and the RHONN model (7) are initially at the same state \(\varvec{\chi }_{j}^{i}(0)={{\varvec{x}}}_{j}^{i}(0)\). Then, for any \(\varepsilon > 0\) and any finite \(T>0\) there exists an integer L and a matrix \({{\varvec{W}}}^{*} \in \mathfrak {R}^{L \times j}\) such that the state \({{\varvec{x}}}_{j}^{i}(t)\) of the RHONN model (7), with L high–order connections and weight values \({{\varvec{W}}} = {{\varvec{W}}}^{*}\), satisfies

$$\begin{aligned} {\text {sup}}_{0 \le t \le T}\ |{{\varvec{x}}}_{j}^{i}(t)-\varvec{\chi }_{j}^{i}(t)| \le \varepsilon . \end{aligned}$$

The above theorem, which is strictly an existence result and does not provide any constructive method in order to obtain the correct synaptic weight matrix \({{\varvec{W}}}^{*}\), proves that if a sufficiently large number of high–order connections are allowed in (7), then it is possible to approximate any dynamical system to any degree of accuracy.

2.2 Filtered Error Training Algorithm

Under the assumption that the unknown system to be controlled is exactly modeled by a RHONN architecture of the form (6), the weight adjustment law and the FE training algorithm for this RHONN are next summarized. Based on the assumptions of no modeling error, there exist unknown weight vectors \({{{\varvec{w}}}^{*}}_{j}^{i}\) such that each state \(\varvec{\chi }_{j}^{i}\) of the unknown dynamical system (8) satisfies

$$\begin{aligned} \dot{\varvec{\chi }}_{j}^{i} = -a_{j}^{i} \varvec{\chi }_{j}^{i} + ({{\varvec{w}}}_{j}^{i*})^{\top } {{\varvec{z}}}_{j}^{i}(\varvec{\chi },{{\varvec{u}}}),\quad \varvec{\chi }_{j}^{i}(0)=\varvec{\chi }_{j0}^{i} \, \end{aligned}$$
(9)

where \(\varvec{\chi }_{j0}^{i}\) is the initial state of the system. As it is standard in systems identification procedures, here it is assumed that the input \({{\varvec{u}}}^{i}(t)\) and the state \(\varvec{\chi }_{j}^{i}(t)\) remain bounded for all \(t \ge 0\). Based on the definition for \({{\varvec{z}}}(\varvec{\chi }_{j}^{i},{\varvec{u}}^{i})\) given by (4), this implies that \({{\varvec{z}}}(\varvec{\chi }_{j}^{i},{{\varvec{u}}}^{i})\) is also bounded. In the sequel, unless there exists confusion, the arguments of the vector field \({{\varvec{z}}}\) will be omitted. Next, the approach for estimating the unknown parameters \({{\varvec{w}}}_{j}^{i*}\) of the RHONN model (9) is described.

Considering (9) as the differential equation describing the dynamics of the unknown system to be controlled, the identifier structure is chosen with the same form as in (6), where \({{\varvec{w}}}_{j}^{i}\) is the estimate of the unknown weight vector \({{\varvec{w}}}_{j}^{i*}\). The identifier structure constructs the regressor \({{\varvec{z}}}\) and other signals related by the identification error equation [23]

$$\begin{aligned} \varvec{\xi }_{j}^{i}= {{\varvec{x}}}_{j}^{i} - \varvec{\chi }_{j}^{i} \, \end{aligned}$$
(10)

which satisfies the differential equation

$$\begin{aligned} \dot{\varvec{\xi }}_{j}^{i} = \dot{{{\varvec{x}}}}_{j}^{i} - \dot{\varvec{\chi }}_{j}^{i} \, \end{aligned}$$
(11)

when considering the parameterization of the plant in a series–parallel configuration for parameter estimation [24], i.e., from substituting (6) and (9) in (11) this latter yields

$$\begin{aligned} \dot{\varvec{\xi }}_{j}^{i} = -a_{j}^{i}{{\varvec{x}}}_{j}^{i}+({{{\varvec{w}}}}_{j}^{i})^{\top }{{\varvec{z}}}_{j}^{i}-\left( -a_{j}^{i}\varvec{\chi }_{j}^{i}+\left( {{{\varvec{w}}}^{i*}_{j}}\right) ^{\top } {{\varvec{z}}}_{j}^{i}\right) \ . \end{aligned}$$
(12)

Thus, expanding terms it results

$$\begin{aligned} \dot{\varvec{\xi }}_{j}^{i} = -a_{j}^{i}{{\varvec{x}}}_{j}^{i}+({{\varvec{w}}}_{j}^{i})^{\top }{{\varvec{z}}}_{j}^{i}+a_{j}^{i}\varvec{\chi }_{j}^{i}-({{\varvec{w}}}^{i*}_{j})^{\top }{{\varvec{z}}}_{j}^{i} \ . \end{aligned}$$
(13)

Then, rearranging terms it can be rewritten as

$$\begin{aligned} \dot{\varvec{\xi }}_{j}^{i} = -a_{j}^{i}({{\varvec{x}}}_{j}^{i}-\varvec{\chi }_{j}^{i})+(({{\varvec{w}}}_{j}^{i})^{\top }-({{\varvec{w}}}_{j}^{i*})^{\top }){{\varvec{z}}}_{j}^{i} \, \end{aligned}$$
(14)

yielding

$$\begin{aligned} \dot{\varvec{\xi }}_{j}^{i} = -a_{j}^{i}\varvec{\xi }_{j}^{i}+({\tilde{{{\varvec{w}}}}}_{j}^{i})^{\top } {{\varvec{z}}}_{j}^{i} \, \end{aligned}$$
(15)

where \((\tilde{{{\varvec{w}}}}_{j}^{i})^{\top }=({{\varvec{w}}}_{j}^{i})^{\top }-({{\varvec{w}}}_{j}^{i*})^{\top }\) denotes the parametric error.

In the series-parallel configuration, the neural network and the plant receive the same external inputs, but the output of the plant is part of the inputs to the neural network. Thus, the dynamic behavior of the neural network is affected by the plant [25]. The series-parallel configuration admits the generation of stable adaptive laws for estimating the unknown parameters of the plant. It uses the actual plant output rather than the estimate of the plant output in the generation of the estimates for some signals of the regressor vector. Since a part of the observer related to the output is in series with the plant and the part related to the input is in parallel with it, the observer is referred to as series-parallel observer [26]. The estimation method based on a series-parallel observer is referred to as the equation error method [27]. The equation error method is based on the principle of least squares which minimizes a quadratic cost function of the error in the state equations to estimate the parameters when assuming that the states, their derivatives, and control inputs are available or accurately measured. It is applicable to linear as well as linear–in–parameters systems. In addition, since it does not involve any numerical integration of the dynamic system, that would otherwise cause divergence, it can be applied to unstable systems [28].

The weights \({{\varvec{w}}}_{j}^{i}\) are adjusted online according to the learning law

$$\begin{aligned} \dot{{{\varvec{w}}}}_{j}^{i}=-\varvec{{\varGamma }}_{j}^{i}{{\varvec{z}}}_{j}^{i}\varvec{\xi }_{j}^{i} \, \end{aligned}$$
(16)

where the adaptive gain, or learning rate, \(\varvec{{\varGamma }}_{j}^{i} \in \mathfrak {R}^{L \times L}\) is a fixed positive definite matrix chosen appropriately to slow down or speed up the rate of adaptation of the parameters. In discrete–time the corresponding adaptive gain, sometimes referred to as step size, needs to be sufficiently small in order to guarantee convergence; however, in continuous–time \(\varvec{{\varGamma }}_{j}^{i}\) simply needs to be positive definite, due to the infinitesimal change of the derivative \(\dot{{{\varvec{w}}}}_{j}^{i}\) [21].

The following theorem establishes that the identification scheme described above has convergence properties with the gradient method for adjusting the weights.

Theorem 2

[4] Consider the filtered error RHONN model given by (9) whose weights are adjusted according to (16). Then,

  1. 1.

    \(\varvec{\xi }_{j}^{i}\), \({\tilde{{{\varvec{w}}}}}_{j}^{i} \in \mathcal {L_{\infty }}\) (i.e., \(\varvec{\xi }_{j}^{i}\) and \({\tilde{{{\varvec{w}}}}}_{j}^{i}\) are uniformly bounded).

  2. 2.

    \(lim_{t \rightarrow \infty } \varvec{\xi }_{j}^{i}(t)=0\).

Stability and convergence properties of the weight adjustment law given by (16) are well known in the adaptive control literature [29, 30]. By proposing a candidate Lyapunov function in terms of the identification error and parametric error, then using (15) and (16), its time derivative results negative semidefinite where only the quadratic term for the identification error appears and can be showed that it is square integrable. Therefore, by applying Barbalat’s lemma, \(lim_{t \rightarrow \infty } \varvec{\xi }_{j}^{i}(t)=0\).

We use a neural model due to the fact that the parameters and external disturbances of the plant to be identified are unknown; then, to approximate the plant complete model, we propose a neural identifier which reproduces the behavior of the plant, without the identification of the plant parameters [31].

2.3 RHONN Structure

As it is well known, the model for a robot manipulator of n DOF is given in the following form

$$\begin{aligned} \varvec{\tau }={{\varvec{M}}}({{\varvec{q}}})\ddot{{{\varvec{q}}}}+ {{\varvec{C}}}({{\varvec{q}}},\dot{{{\varvec{q}}}})\dot{{{\varvec{q}}}}+ {{\varvec{g}}}({{\varvec{q}}})+ {{\varvec{f}}}(\dot{{{\varvec{q}}}}), \end{aligned}$$
(17)

where \(\varvec{\tau }\in \mathfrak {R}^{n}\) is an input vector that represents the torques applied to each joint; \({{\varvec{q}}},\dot{{{\varvec{q}}}}\), \(\ddot{{{\varvec{q}}}} \in \mathfrak {R}^n\) are the states of the system corresponding to position, velocity and acceleration, respectively; \({{\varvec{M}}}({{\varvec{q}}})\ddot{{{\varvec{q}}}}\)\(\in \mathfrak {R}^{n{\times }n}\) represents the contribution of the inertial forces to the dynamical equation; hence the matrix \({{\varvec{M}}}({{\varvec{q}}})\) is the inertia matrix of the robot manipulator; \({{\varvec{C}}}({{\varvec{q}}},\dot{{{\varvec{q}}}})\dot{{{\varvec{q}}}}\)\(\in \mathfrak {R}^n\) represents the Coriolis forces, \({{\varvec{g}}}({{\varvec{q}}}) \in \mathfrak {R}^n\) represents the gravitational forces, and \({{\varvec{f}}}(\dot{{{\varvec{q}}}}) \in \mathfrak {R}^n\) is a vector that combines both viscous and Coulomb friction terms, i.e., the so-called friction vector [29, 30].

The efficiency of an identification procedure depends mainly on the following: error and speed of convergence, stability in cases of abrupt input changes, and its performance during the online training. In the literature, RHONN schemes in conjunction with either EKF [32] or FE algorithm [9, 31], for discrete–time and continuous–time respectively, have shown to be an advantageous approach for identification and control of nonlinear systems. The RHONN identification structure for a 2 DOF robotic arm is proposed as

$$\begin{aligned} \begin{aligned} \dot{{{\varvec{x}}}}_1&=-{{\varvec{A}}}_1 {{\varvec{x}}}_1+{{\varvec{W}}}_1 {{\varvec{z}}}_1+{{\varvec{x}}}_{2} \, \\ \dot{{{\varvec{x}}}}_{2}&=-{{\varvec{A}}}_{2}{{\varvec{x}}}_{2} +{{\varvec{W}}}_{2}{{\varvec{z}}}_{2}+{{\varvec{B}}}{{\varvec{u}}} \, \end{aligned} \end{aligned}$$
(18)

where \({{\varvec{A}}}_1\), \({{\varvec{A}}}_{2}\), and \({{\varvec{B}}} \in {\mathfrak {R}}^{n \times n}\) are diagonal matrices which contain the fixed parameters of the neural network, \({{\varvec{x}}}_1=[x_{1}^1 \; x_{1}^2]^\top \) and \({{\varvec{x}}}_{2}=[x_{2}^1 \; x_{2}^2]^\top \) are neural network state vectors, where \({{\varvec{x}}}_1\) approximates the measurable position vector for the robot manipulator, i.e., \({{\varvec{q}}}=[q^1 \; q^2]^\top \) and \({{\varvec{x}}}_{2}\) approximates the calculated velocity vector \(\dot{{{\varvec{q}}}}=[{\dot{q}}^1 \; {\dot{q}}^2]^\top \), and \({{\varvec{u}}}=[u^1\ u^2]^\top \) is the vector of control inputs. It should be noticed that the structure (18) is similar to that in (7) but given in a block strict-feedback form (BSFF) [33] or regular form [34]. The weight matrices, \({{\varvec{W}}}_1\) and \({{\varvec{W}}}_{2} \in {\mathfrak {R}}^{n \times 2L}\), are defined as follows

$$\begin{aligned} {{\varvec{W}}}_1= & {} \begin{bmatrix} {{\varvec{w}}}_{1}^{1\top }&0_{1\times 5} \\ 0_{1\times 5}&{{\varvec{w}}}_{1}^{2\top } \\ \end{bmatrix}, \end{aligned}$$
(19)
$$\begin{aligned} {{\varvec{W}}}_2= & {} \begin{bmatrix} {{\varvec{w}}}_{2}^{1\top }&0_{1\times 5} \\ 0_{1\times 5}&{{\varvec{w}}}_{2}^{2\top } \\ \end{bmatrix}, \end{aligned}$$
(20)

where \({{\varvec{w}}}_{1}^1\) and \({{\varvec{w}}}_{1}^2 \in {\mathfrak {R}}^{L}\) represent the weight vectors for the neurons that approximate position in both links one and two, respectively. \({{\varvec{w}}}_{2}^1\) and \({{\varvec{w}}}_{2}^2\) correspond to the weight vectors for the neurons that approximate velocity in both links one and two, respectively. The subscript \(1\times 5\) denotes the dimension for a zero vector.

Higher-order interconnections for the NN are established through vectors \({{\varvec{z}}}_1\) and \({{\varvec{z}}}_2 \in {\mathfrak {R}}^{2L}\), which are defined as follows

$$\begin{aligned} \begin{aligned} {{\varvec{z}}}_1&=[{{\varvec{z}}}_{1}^1 \, {{\varvec{z}}}_{1}^2]^\top ,\\ {{\varvec{z}}}_2&=[{{\varvec{z}}}_{2}^1 \; {{\varvec{z}}}_{2}^2]^\top , \end{aligned} \end{aligned}$$
(21)

with

$$\begin{aligned} {{\varvec{z}}}_{1}^1= & {} [s(q^1)s(q^1)\quad s(q^1)\quad s(q^2)\quad s(q^1)s(q^2)\quad 0]^{\top }, \\ {{\varvec{z}}}_{1}^2= & {} [s(q^2)s(q^2)\quad s(q^2)\quad s(q^1)\quad s(q^1)s(q^2)\quad 0]^{\top }, \\ {{\varvec{z}}}_{2}^1= & {} [s({\dot{q}}^1)\quad s({\dot{q}}^1)s({\dot{q}}^1)\quad s({\dot{q}}^1)s(q^1)\quad s({\dot{q}}^1)s(q^2)\quad s({\dot{q}}^1)s({\dot{q}}^2)]^{\top }, \\ {{\varvec{z}}}_{2}^2= & {} [s({\dot{q}}^2)\quad s({\dot{q}}^2)s({\dot{q}}^2)\quad s({\dot{q}}^2)s(q^1)\quad s({\dot{q}}^2)s(q^2)\quad s({\dot{q}}^1)s({\dot{q}}^2)]^{\top }. \end{aligned}$$

The identification error between the neural identifier and the joint variables for each i-th joint is defined as

$$\begin{aligned} \begin{aligned} \varvec{\xi }_{1}&={{\varvec{x}}}_1-{{\varvec{q}}} \, \\ \varvec{\xi }_{2}&={{\varvec{x}}}_2-{\dot{{{\varvec{q}}}}} \, \end{aligned} \end{aligned}$$
(22)

where \(\varvec{\xi }_{1}\) and \(\varvec{\xi }_{2}\) are the identification error for position and velocity, respectively. Then, from (16) the adaptive learning laws are established as

$$\begin{aligned} \begin{aligned} \dot{{{\varvec{w}}}}_{1}^i&=-\varvec{{\varGamma }}_{1}{{\varvec{Z}}}_{1}^{i}\varvec{\xi }_{1} \, \\ \dot{{{\varvec{w}}}}_{2}^i&=-\varvec{{\varGamma }}_{2}{{\varvec{Z}}}_{2}^{i}\varvec{\xi }_{2} \, \end{aligned} \end{aligned}$$
(23)

with \({{\varvec{Z}}}_{1}^{i}=[{{\varvec{z}}}_{1}^{1} \, {{\varvec{z}}}_{1}^{2}]\) and \({{\varvec{Z}}}_{2}^{i}=[{{\varvec{z}}}_{2}^{1} \, {{\varvec{z}}}_{2}^{2}]\).

3 Neural Controller Design

Denoting \(\varvec{\xi }_{1}\) as the identification error vector and the trajectory tracking error between the neural network’s states for position and the desired trajectory as \(\varvec{\epsilon }_{1}= |{{\varvec{x}}}_{1}-{{\varvec{q}}}_{d}|\), where \({{\varvec{q}}}_d=[q_d^1\ q_d^2]^\top \) is the desired trajectory tracking for each joint, the output tracking error is rewritten as

$$\begin{aligned} \tilde{{{\varvec{q}}}}= \varvec{\xi }_{1}+\varvec{\epsilon }_{1} \, \end{aligned}$$
(24)

where \(\tilde{{{\varvec{q}}}}=|{{\varvec{q}}}-{{\varvec{q}}}_{d}|\). Consequently, the error dynamics is given as

$$\begin{aligned} \dot{\tilde{{{\varvec{q}}}}}= \dot{\varvec{\xi }}_{1}+\dot{\varvec{\epsilon }}_{1}. \end{aligned}$$
(25)

Considering (15), (18) and

$$\begin{aligned} \dot{\varvec{\epsilon }}_{1}= \dot{{{\varvec{x}}}}_{1}-\dot{{{\varvec{q}}}}_{d} \, \end{aligned}$$
(26)

the error dynamics is then described by

$$\begin{aligned} \dot{\tilde{{{\varvec{q}}}}}= & {} -{{\varvec{A}}}_{1}\varvec{\xi }_1 +\tilde{{{\varvec{w}}}}_{1}^\top {{\varvec{z}}}_{1}-A_{1}{{\varvec{x}}}_{1} +{{\varvec{W}}}_{1}{{\varvec{z}}}_{1}+{{\varvec{x}}}_{2}-\dot{{{\varvec{q}}}}_{d} \ . \end{aligned}$$

Considering the BSFF system defined by (18), the integrator backstepping procedure is used to construct a stabilizing controller. Assuming that the system is Lipschitz, we proceed by proposing the augmented Lyapunov-like function candidate

$$\begin{aligned} {V}_{1}(\varvec{\epsilon }_{1},\varvec{\xi }_{1},\tilde{{{\varvec{w}}}}_{1}) ={V}(\varvec{\xi }_{1},\tilde{{{\varvec{w}}}}_{1}) +\frac{1}{2}\varvec{\epsilon }_{1}^\top \varvec{\epsilon }_{1} \, \end{aligned}$$
(27)

with

$$\begin{aligned} V(\varvec{\xi }_{1},\tilde{{{\varvec{w}}}}_{1})= \frac{1}{2}\sum _{i=1}^n\big (\varvec{\xi }_{1}^\top \varvec{\xi }_{1} +{\tilde{{{\varvec{w}}}}_{1}}^\top {\varvec{{\varGamma }}_{1}}^{-1}\tilde{{{\varvec{w}}}}_{1}\big )\ , \end{aligned}$$
(28)

taken from [4].

Considering (15), (16), the time derivative of (28) is then given by [4]

$$\begin{aligned} {\dot{V}}(\cdot )=-\sum _{i=1}^{n}{{\varvec{A}}}_1\varvec{\xi }_{1}^\top \varvec{\xi }_{1}\le 0. \end{aligned}$$
(29)

The time derivative of \(V_{1}(\cdot )\) along the solution of (26) is

$$\begin{aligned} {\dot{V}}_1(\cdot )={\dot{V}}(\cdot )+ \varvec{\epsilon }_{1}^\top \big [-{{\varvec{A}}}_{1}{{\varvec{x}}}_{1}+ {{\varvec{W}}}_{1}{{\varvec{z}}}_{1}+{{\varvec{x}}}_{2}-\dot{{{\varvec{q}}}}_{d}\big ] . \end{aligned}$$
(30)

In order to guarantee that (30) to be negative definite, the desired dynamics for \({{\varvec{x}}}_{2}\) is proposed as

$$\begin{aligned} {{\varvec{x}}}_{2d}=-\varvec{\zeta }_{1}\varvec{\epsilon }_{1} +{{\varvec{A}}}_{1}{{\varvec{x}}}_{1} -{{\varvec{W}}}_{1}{{\varvec{z}}}_{1}+\dot{{{\varvec{q}}}}_{d} . \end{aligned}$$
(31)

Now, proposing a new augmented Lyapunov–like function candidate

$$\begin{aligned} V_{2}(\varvec{\epsilon }_1,\varvec{\xi }_1,\tilde{{{\varvec{w}}}}_1,\dot{\tilde{{{\varvec{q}}}}}) =V(\cdot )+V_1(\cdot )+ \frac{1}{2}{\dot{\tilde{{{\varvec{q}}}}}}^\top {\dot{\tilde{{{\varvec{q}}}}}} \, \end{aligned}$$
(32)

and introducing the error

$$\begin{aligned} \dot{\tilde{{{\varvec{q}}}}}={{\varvec{x}}}_{2}-{{\varvec{x}}}_{2d} \, \end{aligned}$$
(33)

the time derivative of (32) is given by

$$\begin{aligned} {\dot{V}}_2(\cdot )= & {} {\dot{V}}(\cdot ) +{\dot{V}}_1(\cdot )+{\dot{\tilde{{{\varvec{q}}}}}}^\top \ddot{\tilde{{{\varvec{q}}}}} . \end{aligned}$$
(34)

From (33), we have

$$\begin{aligned} \ddot{\tilde{{{\varvec{q}}}}}= & {} \dot{{{\varvec{x}}}}_{2}-\dot{{{\varvec{x}}}}_{2d} \, \end{aligned}$$
(35)

which, considering (18), is rewritten as

$$\begin{aligned} \ddot{\tilde{{{\varvec{q}}}}}= & {} -{{\varvec{A}}}_{2}{{\varvec{x}}}_{2} +{{\varvec{W}}}_{2}{{\varvec{z}}}_{2}+{{\varvec{B}}}{{\varvec{u}}}- \dot{{{\varvec{x}}}}_{2d} . \end{aligned}$$

From all of the above, (34) takes the form

$$\begin{aligned} {\dot{V}}_2(\cdot )= & {} {\dot{V}}(\cdot )+ {\varvec{\epsilon }_{1}}^\top \big [-\varvec{\zeta }_{1}\varvec{\epsilon }_{1}+ \dot{\tilde{{{\varvec{q}}}}}\big ]+ {\dot{\tilde{{{\varvec{q}}}}}}^\top \big [-{{\varvec{A}}}_{2}{{\varvec{x}}}_{2} +{{\varvec{W}}}_{2}{{\varvec{z}}}_{2} +{{\varvec{B}}}{{\varvec{u}}}-\dot{{{\varvec{x}}}}_{2d}\big ], \end{aligned}$$

or equivalently

$$\begin{aligned} {\dot{V}}_2(\cdot )= & {} {\dot{V}}(\cdot )-\varvec{\epsilon }_{1}^\top \varvec{\zeta }_{1}\varvec{\epsilon }_{1} +{\dot{\tilde{{{\varvec{q}}}}}}^\top \big [-{{\varvec{A}}}_{2} {{\varvec{x}}}_{2} +{{\varvec{W}}}_{2}{{\varvec{z}}}_{2} +{{\varvec{B}}}{{\varvec{u}}}-\dot{{{\varvec{x}}}}_{2d} +\varvec{\epsilon }_{1}\big ]. \end{aligned}$$
(36)

Then, from (31) we have

$$\begin{aligned} \dot{{{\varvec{x}}}}_{2d}= & {} -\varvec{\zeta }_{1}\dot{\varvec{\epsilon }}_{1} +{{\varvec{A}}}_{1}\dot{{{\varvec{x}}}}_{1}-(\dot{{{\varvec{W}}}}_{1}{{\varvec{z}}}_{1} +{{\varvec{W}}}_{1}\dot{{{\varvec{z}}}}_{1})+\ddot{{{\varvec{q}}}}_{d} \, \end{aligned}$$

which, considering (26) and (18), it can be rewritten as

$$\begin{aligned} \dot{{{\varvec{x}}}}_{2d}= & {} ({{\varvec{A}}}_{1}-\varvec{\zeta }_{1})(-{{\varvec{A}}}_{1}{{\varvec{x}}}_{1}+{{\varvec{W}}}_{1}{{\varvec{z}}}_{1}+{{\varvec{x}}}_{2})\\&-\dot{{{\varvec{W}}}}_{1}{{\varvec{z}}}_{1}-{{\varvec{W}}}_{1}\frac{\partial ({{\varvec{z}}}_{1})}{\partial {{{\varvec{q}}}}}\dot{{{\varvec{q}}}} +\ddot{{{\varvec{q}}}}_{d}+\varvec{\zeta }_{1}\dot{{{\varvec{q}}}}_{d} . \end{aligned}$$

Finally, for (36) be negative definite, the control law is then proposed as

$$\begin{aligned} {{\varvec{u}}}= & {} {{\varvec{B}}}^{-1}[{{\varvec{A}}}_{2}{{\varvec{x}}}_{2} -{{\varvec{W}}}_{2}{{\varvec{z}}}_{2}+\dot{{{\varvec{x}}}}_{2d}-\varvec{\epsilon }_{1} -\varvec{\zeta }_{2}\dot{\tilde{{{\varvec{q}}}}}] . \end{aligned}$$
(37)

4 Results

This section gives a description of the physical plant and real-time results that are exhibited in order to evaluate the effectiveness of the neural controller scheme proposed. The platform of the 2 DOF robotic arm used for this work is depicted in Fig. 1.

Fig. 1
figure 1

Platform of the 2 DOF robotic arm

4.1 Robot Description

In order to evaluate in real–time the performance of the neural block control scheme proposed, it is implemented on a 2 DOF robotic arm, of our own design and unknown parameters shown in Fig. 2, whose displacements evolve on the vertical plane. The robotic arm consists of two rigid links, with a brushless direct-drive servo for each link; the first one with a 16:1 gear reduction for link 1 and direct connection of the last one for link 2. The robotic arm is constituted by the MTR-70-3S-AA and MTR-55-3S-AA servomotors, manufactured by FESTO Pneumatic AG, for link 1 and link 2 respectively. Incremental encoders embedded on the servomotors deliver information, via RS-422 protocol, related with the angular displacements. Both servomotors exhibit a resolution of 4096 (pulses/rev), i.e. an accuracy of 0.00153398 (rad/pulse). The angular velocities are computed via numerical differentiation from the angular position measurement.

According to the actuator’s manufacturer, the direct-drive servomotors are able to supply torques within the following bound

$$\begin{aligned} |u| \le \tau _{\text {max}}= & {} 2\ \text{ Nm }. \end{aligned}$$
(38)

Considering that servo-drives are controlled by means of a data acquisition board (DAB) using two analog channels of \(\pm \,10\) (VDC), with a sampling time of 0.0005 (s), the voltage–torque relationship has been established as it is shown in Table 1.

Fig. 2
figure 2

Direct drive vertical robotic arm test bed

Table 1 Admissible signals by the servomotors

4.2 Tracking Neural Control

For validating the experimental results, we define the following bounded desired trajectories, which possess special characteristics (as it is done in [30])

$$\begin{aligned} q_{d}^1= & {} r_1(1-e^{d_1t^3})+c_1(1-e^{d_1t^3})\text {sin}(\omega _1t) \, \\ q_{d}^2= & {} r_2(1-e^{d_2t^3})+c_2(1-e^{d_2t^3})\text {sin}(\omega _2t) \, \end{aligned}$$

where \(r_1=\pi /2\), \(c_1=2\pi /18\), \(d_{1}=-2\), and \(\omega _1=2\) (rad/s) are parameters of the desired position trajectory for the first joint; whereas \(r_2=\pi /2\), \(c_2=25\pi /36\), \(d_{2}=-1.8\) and \(\omega _2=1\) (rad/s) are parameters of the desired position trajectory for the second joint.

These trajectories have the following characteristics: (a) they include a step-like end of small magnitude such that it is possible to show the transient response of the controller; (b) they incorporate a sinusoidal term to evaluate the provision to relatively rapid periodic signals, where the nonlinearities of the robot dynamics are really important and (c) they present a term that gently increases to keep the robot in a state of operation without saturating the actuators.

4.3 Real-Time Results

To evaluate the efficiency and speed of the identification RHONN scheme and the performance and robustness of the neural controller, we propose the next experiment under the following arbitrary criterion:

For joints 1 and 2: zero initial conditions for both plant and RHONN; with a delay of 2 s for the controller of the joint 1 and a delay of 4 s for the controller of the joint 2.

Figures 3 and 4 display the neural identification and control results for trajectory tracking for both joints when considering delays to verify the robustness of the neural stabilizing controller for any point of the trajectories. Figures 5 and 6 show the identification and tracking errors for each joint respectively, and the applied torques are shown in Fig. 7. It can be noticed that the control signals are maintained inside of the prescribed limits given by the actuators manufacturer. The use of an identification RHONN scheme reduces the time and considerably minimizes the identification error. This is due to the consideration of high–order interconnections which makes the difference with respect to conventional RNN. Figures 8 and 9 show the trajectories of the weights according to positions 1 and 2, respectively.

Fig. 3
figure 3

Neural identification and control of the plant behavior for joint 1

Fig. 4
figure 4

Neural identification and control of the plant behavior for joint 2

Fig. 5
figure 5

Identification and trajectory tracking errors for link 1

Fig. 6
figure 6

Identification and trajectory tracking errors for link 2

To quantify the control performance we use the root–mean square (RMS) average of tracking error, based on the \({\mathcal {L}}_{2}\) norm of the tracking error \({\tilde{q}}\), given by

$$\begin{aligned} {\mathcal {L}}_{2}[{\tilde{q}}]=\sqrt{\frac{1}{T-t_{0}}\int _{t_0}^{T}{\tilde{q}}^{T} {\tilde{q}}\ dt} \end{aligned}$$
(39)

where T is the total experimentation time and \(t_{0}\) is for the initial time.

A comparative analysis between the control scheme RHONNBC (Recurrent High–Order Neural Network Block Controller) here proposed and RNNBC (Recurrent Neural Network Block Controller) designed from [7] but in a centralized way is included in Table 2.

Table 2 Comparison of the RMS for the joint positions of the robotic arm
Fig. 7
figure 7

Applied torques

Fig. 8
figure 8

Trajectories of the weights for position 1

Fig. 9
figure 9

Trajectories of the weights for position 2

5 Conclusion

The real–time results validate the effectiveness of the continuous–time neural block control scheme proposed for trajectory tracking of a robotic arm. Moreover, it is shown that effects due to friction and gravitational forces, which can be seen as physical interconnections between the joints of the whole system, are absorbed by the neural identifier. Furthermore, the neural stabilizing controller exhibits a good performance in spite of forced delays. As future work, we will deal with the implementation of a neural sliding mode control scheme from the RHONN model proposed as well as extending our results to controlling a 3 DOF robotic arm.