Adaptive dynamic programming-based controller with admittance adaptation for robot–environment interaction

Abstract

The problem of optimal tracking control for robot–environment interaction is studied in this article. The environment is regarded as a linear system and an admittance control with iterative linear quadratic regulator method is obtained to guarantee the compliant behaviour. Meanwhile, an adaptive dynamic programming-based controller is proposed. Under adaptive dynamic programming frame, the critic network is performed with radial basis function neural network to approximate the optimal cost, and the neural network weight updating law is incorporated with an additional stabilizing term to eliminate the requirement for the initial admissible control. The stability of the system is proved by Lyapunov theorem. The simulation results demonstrate the effectiveness of the proposed control scheme.

Introduction

Robot applications are becoming more and more widespread, such as rehabilitation therapy, assembly automation and surgery. ^{1 –4} They can either work independently to accomplish tasks or cooperate with their human partners for certain tasks. In the actual application process, the robot will inevitably interact with the external environments. ^{5 –7} Consequently, in recent years, interaction control between the robot and environment has attracted great concern and is considered to be greatly important.

In existing research, two main approaches are applied to achieve compliant behaviour of the robot, that is, hybrid position/force control and impedance control. ^8,9 The first approach requires the position subspace and force subspace decomposition, task planning and control law switching in the execution process. Without considering the dynamic coupling of the environment and the robot, the accuracy of the hybrid position/force control cannot be guaranteed. ¹⁰ In contrast, the second approach aims to adjust the mechanical impedance to a target one, which will guarantee the robot to be complaint with the interaction force imposed by the external environment. Impedance control ensures the safety of the robot and the environment and it has been proved to be more feasible and has better robustness. According to the causality of the controller, impedance control has two implementation methods, one is named impedance control and the other is admittance control. In impedance control system, the interaction force can be estimated from the desired motion trajectory and impedance model, while in admittance control system, the reference trajectory is obtained from the measured environmental external force and the desired admittance model. Therefore, in this article, admittance control is adopted to solve the problem of robot–environment interaction control.

In admittance control system, force and the admittance model are two important parts. When robot–environment interaction exists, the force can be detected and measured by the sensors installed on the end-effector of the robot arm. But, how to derive optimal parameters of the admittance model is non-trivial. On the one hand, it is usually difficult to derive the desired admittance model because of the complexity of environmental dynamics; on the other hand, a fixed admittance model cannot satisfy all cases. Taking human–robot cooperation as an example, variable admittance control is necessary to ensure more efficient performance. ¹¹ To solve these problems, iterative learning has been studied in robot intelligent control area. It has been investigated to obtain admittance parameters to adapt to unknown environment. The aim of this approach is to introduce human learning skills into the robot and improve control performance by repeating a task. Cohen and Flash ¹² proposed an impedance learning control scheme using an associative search network to complete a wall-following work. Neural network (NN) is introduced into the impedance control to regulate the parameters. ¹³ However, the iterative learning method requires the robot to operate repeatedly, which brings inconvenience in practical process and is not feasible in many situations. Love and Book, ¹⁴ Uemura and Kawamura, ¹⁵ Gribovskaya et al., ¹⁶ Stanisic and Fernández, ¹⁷ Landi et al. ¹⁸ and Yao et al. ¹⁹ have proposed to utilize adaptation approaches to address the problems stated above.

Robotic motion control is a challenging task as it is difficult to obtain accurate model concerning that the robot is a non-linear and highly coupled system. Proportional–integral–derivative (PID) control, NN control, adaptive control and other control methods have been applied to the robot system. ^{20 –27} As a classical control method, PID control is employed to the robot system and can track the given reference trajectory well. ²⁸ It is acknowledged that PID control has some advantages, such as simple structure and good robustness, but it is not easy to select suitable PID parameters if the controlled plant is complex. In addition, when dynamic uncertainties exist in the system, PID control cannot satisfy the performance requirements for the magnitude of overshoot, the rising and settling time and so on. NN has the fundamental characteristics of human brain and can simulate human behaviour for information processing, therefore it is widely used in the control field for unknown system identification. NN control can model the uncertain dynamics online to improve the system performance. ²⁹ An admittance adaptation method and the NN-based controller are applied into the robot system. ³⁰

Tracking control is a significant research issue in the domain of robot intelligent control. For a controlled system, stability is just the minimum requirement. Optimal control needs to be considered, that is, it is required to design an optimal tracking controller, which could ensure system stability of the robot while minimizing the cost function. Werbos ³¹ proposed adaptive dynamic programming (ADP) strategy and it is considered to be an effective approach to resolve the optimal control problem. ³² The key of ADP method is to find a solution of Hamilton–Jacobi–Bellman (HJB) equation. However, because it is a partial differential equation, when the controlled system is non-linear but not linear, its analytical solution will be very difficult to obtain, or even impossible. To solve the above problem, policy iterative is considered as an effective method to find the approximate solution, which requires initial stability control. ³³ However, in practical process, the initial admissible control is usually very difficult to satisfy. Then, NN is introduced to derive an approximate solution of the HJB equation. The approximate solution is obtained by NN-based method, meanwhile the requirement of initial stability is eliminated with the incorporation of an additional term. ^34,35

Yang et al. ³⁰ paid attention to the robot–environment interaction control, but did not consider the optimization problem. However, for the robot, how to perform path tracking optimization and minimize the cost function is very important. Based on the above discussion, the optimal tracking control problem for robot–environment interaction is studied in this article. Moreover, the admittance control and ADP approach are adopted to improve the system performance. The contributions of this article are listed below:

The environment with unknown dynamics is modelled as a linear system. An admittance adaptation method with iterative linear–quadratic regulator (LQR) is obtained to achieve a compliant behaviour.

ADP approach is introduced into the robot system to solve the optimal tracking problem. The critic network with radial basis function (RBF) is developed to approximate the minimum cost function. In addition, to eliminate the requirement for initial admissible control, a stabilizing term is incorporated into NN weight updating law.

The rest of this article is arranged as follows. Firstly, the robot and environment systems and control objectives are described. Next, the control scheme including admittance adaptation and optimal control using ADP is developed. Then, simulation studies are given. Finally, the conclusion is drawn.

Preliminaries and problem formulation

Robot dynamics

The n-link robot manipulator dynamics is showed as the following Lagrangian form

$M(q)\ddot{q}+C(q,\dot{q})\dot{q}+G(q)=\mu$

1

where $q=[q_1,q_2,\dots ,q_n{]}^{\text{T}}\in {ℝ}^n$ , $\dot{q}=[{\dot{q}}_1,{\dot{q}}_2,\dots ,{\dot{q}}_n{]}^{\text{T}}\in {ℝ}^n$ and $\ddot{q}=[{\ddot{q}}_1,{\ddot{q}}_2,\dots ,{\ddot{q}}_n{]}^{\text{T}}\in {ℝ}^n$ represent the robot position vector, velocity vector and acceleration vector in joint space, respectively. $\mu \in {ℝ}^n$ is the joint torque, while $M(q)\in {ℝ}^{n\times n}$ , $C(q,\dot{q})\in {ℝ}^{n\times n}$ and $G(q)\in {ℝ}^n$ are known matrices and denote the inertial matrix, Coriolis/centrifugal matrix and gravity vector, respectively. For convenience, M, C and G denote the known matrices $M(q)$ , $C(q,\dot{q})$ and $G(q)$ in the following section, respectively.

Define the reference trajectory as $q_r\in {ℝ}^n$ , and the tracking error $q_e\in {ℝ}^n$ is shown as follows

$q_e=q-q_r$

2

Then, the first and second time derivative of q_e are given below

$\begin{array}{ll}{\dot{q}}_e\hfill & =\dot{q}-{\dot{q}}_r\hfill \\ {\ddot{q}}_e\hfill & =\ddot{q}-{\ddot{q}}_r\hfill \end{array}$

3

We define the sliding motion surface ξ as follows

$\xi =\Lambda q_e+{\dot{q}}_e$

4

where $\Lambda \in R^{n\times n}$ is a constant positive matrix. According to equations (2) to (4), we can get

$\begin{array}{ll}\dot{q}=\hfill & \xi -\Lambda q_e+{\dot{q}}_r\hfill \\ \ddot{q}=\hfill & \dot{\xi }-\Lambda q_e+{\ddot{q}}_r\hfill \end{array}$

5

Substituting equation (5) into equation (1), the error dynamics is obtained as follows

$\begin{array}{ll}\dot{\xi }=\hfill & -M^{-1}C(\xi -\Lambda q_e+{\dot{q}}_r)-M^{-1}G\hfill \\ \hfill & -{\ddot{q}}_r+\Lambda {\dot{q}}_e+M^{-1}\mu \hfill \end{array}$

6

Then, the following system is obtained

$\dot{\xi }=f(\xi )+g(\xi )\mu$

7

The non-linear functions $f:{ℝ}^n\to {ℝ}^n$ and $g:{ℝ}^n\to {ℝ}^{n\times n}$ in equation (7) are specified by

$\begin{array}{l}f(\xi )=-M^{-1}C(\xi -\Lambda q_e+{\dot{q}}_r)-M^{-1}G-{\ddot{q}}_r+\Lambda {\dot{q}}_e\\ g(\xi )=M^{-1}\end{array}$

8

Environment dynamics

It is assumed that the dynamics of environmental interaction force subject to the equation given below

$C_E\dot{x}+G_{E}x=-F$

9

where C_E and G_E represent the unknown damping and stiffness of the environment, respectively. F denotes the interaction force and can be detected and measured by a force sensor. x is the end-effector position in Cartesian space and the corresponding desired trajectory x_d is defined as

${\dot{x}}_d=U_d{x}_d$

10

where $U\in {ℝ}^{m\times m}$ is a known matrix. Subsequently, we define $\eta =[x,x_d{]}^{\text{T}}$ . Thus, combining equation (9) with equation (10), dynamics of the unknown environment and the desired trajectory are generated by

$\begin{array}{c}\dot{\eta }=\left[\begin{array}{cc}-C_E^{-1}{G}_E& 0\\ 0& U_d\end{array}\right]\eta +\left[\begin{array}{c}-C_E^{-1}\\ 0\end{array}\right]F\\ =A_e\eta +B_{e}F\end{array}$

11

If we take equation (11) as a linear system with F as its control input and η as its states to be controlled, this equation relates x with x_d via the optimal feedback control law $F=-K_e\eta$ whose aim is to minimize the cost function

${\text{Γ}}_1={\displaystyle {\int }_0^{\infty }}\left(x_e^{\text{T}}{Q}_{E1}{x}_e+F^{\text{T}}{R}_{E}F\right)\text{d}t$

12

This cost function also indicates that our motivation of modifying a desired trajectory x_d is to balance the contact force F with the tracking error $x_e=:x-x_d$ . And this balance can be tuned via the user-defined $Q_{E1}$ and R_E .

In this section, the robot and environment dynamics are modelled. Then, we will design a control strategy to achieve the compliant behaviour and optimal tracking control in case the robot interacts with the environment.

Control scheme

A control scheme consisting of three parts as shown in Figure 1 including an optimal trajectory modifier using admittance control, a closed-loop inverse kinematics (CLIK) solver and a trajectory tracking controller based on ADP technique is designed in this section.

Figure 1.

An illustration of the proposed control scheme.

Trajectory modification using admittance control

The solution to equation (12) is an analogy with the LQR problem. It can be rewritten as

$\begin{array}{l}\text{Γ}={\displaystyle {\int }_0^{\infty }}\left({\eta }^{\text{T}}{Q}_E\eta +F^{\text{T}}{R}_{E}F\right)\text{d}t\\ Q_E=\left[\begin{array}{cc}{Q}_{E1}& -Q_{E1}{U}_d\\ -U_d^{\text{T}}{Q}_{E1}& U_d^{\text{T}}{Q}_{E1}{U}_d\end{array}\right]\end{array}$

13

whose system counterpart is consistent with equation (11).

In this subsection, an algorithm proposed by Jiang and Jiang ³⁶ is adopted to solve the algebraic Riccati equation (ARE) in equation (14) with unknown environment parameters C_E , G_E to derive the feedback gain K_e

$\begin{array}{l}P{A}_e+A_e^{\text{T}}P+Q_E-P{B}_e{R}_E^{-1}{B}_e^{\text{T}}P=0\\ K_e=-R_E^{-1}{B}_e^{\text{T}}P\end{array}$

14

Some notations are outlined here. n, m and d are the length of η, F and the sample times integer, respectively. The sampled signal together with the historical ones comprising the matrix as follows

$\begin{array}{l}\widehat{p}={\left[p_{11}\mathrm{,2}{p}_{12},\dots \mathrm{,2}{p}_{1n},p_{22}\mathrm{,2}{p}_{23},\dots ,p_{nn}\right]}^{\text{T}}\\ \overline{\eta }={\left[{\eta }_1^2,{\eta }_1{\eta }_2,\dots ,{\eta }_1{\eta }_n,{\eta }_2^2,{\eta }_2{\eta }_3,\dots ,{\eta }_n^2\right]}^{\text{T}}\\ d_{\overline{\eta }}={\left[\overline{\eta }\left(t_1\right)-\overline{\eta }\left(t_0\right),\overline{\eta }\left(t_2\right)-\overline{\eta }\left(t_1\right),\dots ,\overline{\eta }\left(t_d\right)-\overline{\eta }\left(t_{d-1}\right)\right]}^{\text{T}}\\ I_{\eta }^{\eta }={\left[{\displaystyle {\int }_{t_0}^{t_1}}\eta \otimes \eta \hspace{0.17em}\text{d}t,{\displaystyle {\int }_{t_1}^{t_2}}\eta \otimes \eta \hspace{0.17em}\text{d}t,\dots ,{\displaystyle {\int }_{t_{d-1}}^{t_d}}\eta \otimes \eta \hspace{0.17em}\text{d}t\right]}^{\text{T}}\\ I_F^{\eta }={\left[{\displaystyle {\int }_{t_0}^{t_1}}\eta \otimes f\hspace{0.17em}\text{d}t,{\displaystyle {\int }_{t_1}^{t_2}}\eta \otimes F\hspace{0.17em}\text{d}t,\dots ,{\displaystyle {\int }_{t_{d-1}}^{t_d}}\eta \otimes F\hspace{0.17em}\text{d}t\right]}^{\text{T}}\end{array}$

15

where $p\in {ℝ}^{\frac{1}{2}n(n+1)}$ , $\overline{\eta }\in {ℝ}^{\frac{1}{2}n(n-1)}$ , $d_{\overline{\eta }}\in {ℝ}^{d\times \frac{1}{2}n(n-1)}$ , $I_{\eta }^{\eta }\in {ℝ}^{d\times n^2}$ , $I_F^{\eta }\in {ℝ}^{d\times nm}$ and ⊗ stand for the Kronecker product, and $p_{ij}$ and ${\eta }_i$ denote entries of P and η, respectively

$\text{rank}\left(\left[I_{\eta }^{\eta },I_F^{\eta }\right]\right)=\frac{n(n+1)}{2}+nm$

16

When the number of sampled data is large enough and the rank condition in equation (16) is satisfied, the algorithm can solve K_e by iteratively calculating equation (17) until ${\widehat{p}}^{(k)}$ converge to an acceptable range ε, that is, $||{\widehat{p}}^{(k)}-{\widehat{p}}^{(k-1)}||<\epsilon$ with ∥ * ∥ denoting the 2-norm of *

$\begin{array}{l}{Q}_E^{(k)}=Q_E+K_e^{(k)\text{T}}{R}_E{K}_e^{(k)}\\ {\Theta }^{(k)}=\left[d_{\overline{\eta }},-2I_{\eta }^{\eta }\left(I_n\otimes K_e^{(k)\text{T}}{R}_E\right)-2I_F^{\eta }\left(I_n\otimes R_E\right)\right]\\ {\Xi }^{(k)}=-I_{\eta }^{\eta }\text{vec}\left(Q_E^{(k)}\right)\\ \left[\begin{array}{l}{\widehat{p}}^{(k)}\hfill \\ \text{vec}\left(K_e^{(k+1)}\right)\hfill \end{array}\right]={\left({\Theta }^{(k)\text{T}}{\Theta }^{(k)}\right)}^{-1}{\Theta }^{(k)\text{T}}{\Xi }^{(k)}\end{array}$

17

where the superscript $(k)$ denotes the index of the iteration, $\text{vec}(*)$ denotes the column vectorization of * and $I_n\in {ℝ}^{n\times n}$ is an identity matrix.

Once the optimal feedback gain K_e is obtained, we can use it to modify x_d . Formulations are given as below

$F=-K_e\eta =-\left[\begin{array}{cc}{K}_{e1}& K_{e2}\end{array}\right]\left[\begin{array}{c}x\\ x_d\end{array}\right]$

18

where $K_{e1}$ and $K_{e2}$ are compatible matrix from K_e . Finally, the modified trajectory x_r to be tracked is calculated, which is equivalent to the x in equation (18)

$x_r=-K_{e1}^{-1}F-K_{e1}^{-1}{K}_{e2}{x}_d$

19

Inverse kinematics using CLIK

The CLIK algorithm is employed to resolve the Cartesian reference trajectory x_r into the one q_r in joint space. ³⁷ Let the solution error $e:=\kappa (q_r)-x_r$ where $\kappa (*)$ denotes the forward kinematics and e is given by

$\dot{e}=-K_{f}e$

20

where K_f is a positive user-defined matrix that decides the convergent rate of e. Expanding the above equations and combining with $\dot{x}=J_{\text{co}}\dot{q}$ and $J_{\text{co}}=\partial \kappa (q)/\partial q$ , the following equation holds

${\dot{q}}_r=J_{\text{co}}^{†}({\dot{x}}_r-K_f(\kappa (q_r)-x_r))$

21

integrating of which yields the CLIK method

$q_r={\displaystyle {\int }_0^t}\left(J^{†}{\dot{x}}_r-J_{\text{co}}^{†}{K}_f(\kappa (q_r)-x_r)\right)\text{d}t$

22

where $q(0)={\kappa }^{-1}(x_r(0))$ , $J_{\text{co}}^{†}=J_{\text{co}}^{\text{T}}{\left(J_{\text{co}}{J}_{\text{co}}^{\text{T}}+\sigma I_n\right)}^{-1}$ , and $\sigma \in ℝ$ is introduced to avoid the singularity problem which is recommended to be assigned small enough for improving the solution accuracy.

Optimal control using ADP

As mentioned in the Introduction section, it is very important to optimize the trajectory tracking while minimizing the design cost for robots. On the basis of optimal theory, the optimal control of the system (7) can be derived by solving the HJB equation in the frame of ADP. Consequently, in this subsection, our target is to find such an optimal control μ.

Assume that the functions $f(\xi )$ and $g(\xi )$ are Lipschitz continuous in ${ℝ}^{2n}$ and system (7) is controllable, then the optimal control ${\mu }^{*}$ should minimize the cost function which is expressed as

$J(\xi (t))={\displaystyle {\int }_t^{\infty }}\left[\text{Φ}(\xi (\tau ))+U(\xi (\tau ),\mu (\xi (\tau )))\right]\text{d}\tau$

23

where $\text{Φ}(\xi (\tau ))=\xi {(\tau )}^{T}Q\xi (\tau )$ , $U(\xi (\tau ),\mu (\xi (\tau )))=\mu {(\xi (\tau ))}^{T}R\mu (\xi (\tau ))$ , $Q\in {ℝ}^{n\times n}$ and $R\in {ℝ}^{n\times n}$ are symmetric positive definite matrices. For robot system (7), the optimal control ${\mu }^{*}$ should not only guarantee system stability but also can make the cost function finite, that is, the control law should be in the admissible control set which defined as χ. Additionally, for any admissible control law $\mu \in \chi$ , if $J(\xi )$ given in equation (23) is continuously differentiable, we will have the non-linear Lyapunov equation which is an infinitesimal version of equation (23) is shown as follows with $J(0)=0$

$\begin{array}{ll}0=\hfill & \text{Φ}(\xi (\tau ))+U(\xi (\tau ),\mu (\xi (\tau )))\hfill \\ \hfill & +{(\nabla J(\xi ))}^{\text{T}}(f(\xi )+g(\xi )\mu (\xi ))\hfill \end{array}$

24

where $J(\xi (t))$ is short for $J(\xi )$ for convenience and the notation $\nabla *\triangleq \frac{\partial *}{\partial \xi }$ denotes the partial derivative of *.

Then, the Hamiltonian function and the optimal cost function of robot system (7) are defined as below

$\begin{array}{ll}H(\xi ,\mu (\xi ),\nabla J(\xi ))=\hfill & \text{Φ}(\xi (\tau ))+U(\xi (\tau ),\mu (\xi (\tau )))\hfill \\ \hfill & +{(\nabla J(\xi ))}^{\text{T}}(f(\xi )+g(\xi )\mu (\xi ))\hfill \end{array}$

25

$J{(\xi )}^{*}=\underset{\mu \in \chi }{{\displaystyle \text{min}}}{\displaystyle {\int }_t^{\infty }}\left[\text{Φ}(\xi (\tau ))+U(\xi (\tau ),\mu (\xi (\tau )))\right]\text{d}\tau$

26

We can obtain the HJB equation shown as

$0=\underset{\mu \in \chi }{{\displaystyle \text{min}}}H(\xi ,\mu (\xi ),\nabla J^{*}(\xi ))$

27

Suppose that the minimum value on the right side of formula (27) exists and also is unique, from $\frac{\partial H(\xi ,\mu (\xi ),\nabla J^{*}(\xi ))}{\partial \mu }=0$ , then the following optimal control ${\mu }^{*}(\xi )$ can be derived as

${\mu }^{*}(\xi )=-\frac{1}{2}{R}^{-1}{g}^{\text{T}}(\xi )\nabla J^{*}(\xi )$

28

Substituting the optimal control law (28) into equation (24) yields another form of HJB equation with respect to $\nabla J^{*}(\xi )$ is obtained as

$H(\xi ,{\mu }^{*}(\xi ),\nabla J^{*}(\xi ))=0$

29

Inspired by Liu et al., ³⁴ we know that if the optimal function $J^{*}(\xi )$ is assumed to be continuously differentiable, $J^{*}(\xi )$ can be rebuilded by RBFNN which can be shown as below

$J^{*}(\xi )=w^{\text{T}}S(\xi )+\epsilon (\xi )$

30

where $w\in {ℝ}^l$ represents the ideal constant weight, $S:{ℝ}^{2n}\to {ℝ}^l$ denotes the activation function, l denotes the node number in the hidden layer and $\epsilon (\xi )$ denotes the unknown approximation error of NN. Then, the derivation of equation (30) involving ξ is derived as

$\nabla J^{*}(\xi )=(\nabla S(\xi {))}^{\text{T}}w+\nabla \epsilon (\xi )$

31

From equations (28) and (31), the following ${\mu }^{*}$ can be obtained as

${\mu }^{*}(\xi )=-\frac{1}{2}{R}^{-1}{g}^{\text{T}}(\xi )((\nabla S(\xi {))}^{\text{T}}w+\nabla \epsilon (\xi ))$

32

Then, substituting equations (31) and (32) into equation (29), we have

$\begin{array}{ll}{H}^{*}(\xi ,{\mu }^{*}(\xi ),\nabla J^{*}(\xi ))=\hfill & \text{Φ}(\xi )+w^{\text{T}}\nabla S(\xi )f(\xi )\hfill \\ \hfill & -\frac{1}{4}{w}^{\text{T}}\nabla S(\xi )D\nabla S{(\xi )}^{\text{T}}w\hfill \\ \hfill & +e_c=0\hfill \end{array}$

33

where

$e_c=(\nabla \epsilon (\xi {))}^{\text{T}}(f(\xi )+g(\xi ){\mu }^{*}(\xi ))$

34

$D=g(\xi )R^{-1}g{(\xi )}^{\text{T}}$

35

In fact, the ideal weight w and $J^{*}(\xi )$ in equation (30) are unknown, then the estimate weight and optimal cost function, respectively, denoted as $\widehat{w}$ and $J^{*}(\xi )$ can be obtained by the constructed critic NN. Therefore, the approximate optimal cost $J^{*}(\xi )$ is given as below

$\widehat{J}(\xi )={\widehat{w}}^{\text{T}}S(\xi )$

36

Then, the derivative of equation (36) is

$\nabla \widehat{J}(\xi )=(\nabla S(\xi {))}^{\text{T}}\widehat{w}$

37

Based on equations (28) and (37), the approximate optimal control is obtained as

$\widehat{\mu }(\xi )=-\frac{1}{2}{R}^{-1}{g}^{\text{T}}(\xi )(\nabla S(\xi {))}^{\text{T}}\widehat{w}$

38

Similarly, applying equations (25), (37) and (38), the approximate Hamiltonian function $\widehat{H}(\xi ,\widehat{\mu }(\xi ),\nabla \widehat{J}(\xi ))$ can be derived as

$\begin{array}{ll}\widehat{H}(\xi ,\widehat{\mu }(\xi ),\nabla \widehat{J}(\xi ))=\hfill & \text{Φ}(\xi )+{\widehat{w}}^{\text{T}}\nabla S(\xi )f(\xi )\hfill \\ \hfill & -\frac{1}{4}{\widehat{w}}^{\text{T}}\nabla S(\xi )D{(\nabla S(\xi ))}^{\text{T}}\widehat{w}\hfill \end{array}$

39

Define e_H as the error between $H^{*}$ and $\widehat{H}$ , $\stackrel{\sim }{w}$ as the approximate NN weight error, then they are shown as below

$\begin{array}{ll}{e}_H=\hfill & \widehat{H}(\xi ,\widehat{\mu }(\xi ),\nabla \widehat{J}(\xi ))\hfill \\ \hfill & -H^{*}(\xi ,{\mu }^{*}(\xi ),\nabla J^{*}(\xi ))\hfill \end{array}$

40

$\stackrel{\sim }{w}=w-\widehat{w}$

41

According to equations (33), (39) and (41), e_H in equation (40) can be described as

$\begin{array}{c}{e}_H=\widehat{H}(\xi ,\widehat{\mu }(\xi ),\nabla \widehat{J}(\xi ))\\ =-{\stackrel{\sim }{w}}^{\text{T}}\nabla S(\xi )f(\xi )+\frac{1}{2}{\stackrel{\sim }{w}}^{\text{T}}\nabla S(\xi )D(\nabla S(\xi {))}^{\text{T}}w\\ -\frac{1}{4}{\stackrel{\sim }{w}}^{\text{T}}\nabla S(\xi )D(\nabla S(\xi {))}^{\text{T}}\stackrel{\sim }{w}-e_c\end{array}$

42

To train RBFNN, an appropriate weight updating law $\widehat{w}$ should be designed to both minimize the objective function $E=\frac{1}{2}{e}_H^2$ and ensure the approximate optimal weight $\widehat{w}$ converge to the ideal weight w. To eliminate the requirement for the initial admissible control law, the weight $\widehat{w}$ is tuned according to the standard gradient descent algorithm with an additional stabilizing term. The weight updating law is given as

$\begin{array}{ll}\dot{\widehat{w}}=\hfill & -(1-h){\alpha }_H\left(\frac{\partial E}{\partial \widehat{w}}\right)\hfill \\ \hfill & +\frac{1}{2}h{\alpha }_c\left(\frac{\partial {(\nabla J_s(\xi ))}^{\text{T}}(f(\xi )+g(\xi )\widehat{\mu })}{\partial \widehat{w}}\right)\hfill \\ =\hfill & -(1-h){\alpha }_H\left(\frac{\partial E}{\partial \widehat{w}}\right)+\frac{1}{2}h{\alpha }_c\nabla S(\xi )D\nabla J_s(\xi )\hfill \end{array}$

43

$\begin{array}{ll}\frac{\partial E}{\partial \widehat{w}}=\hfill & e_H\frac{\partial e_H}{\partial \widehat{w}}\hfill \\ =\hfill & \widehat{H}(\xi ,{\widehat{\mu }}^{*}(\xi ),\nabla \widehat{J}(\xi ))\frac{\partial \widehat{H}}{\partial \widehat{w}}\hfill \\ =\hfill & [\nabla S(\xi )f(\xi )-\frac{1}{2}\nabla S(\xi )D\nabla S{(\xi )}^{\text{T}}\widehat{w}][\text{Φ}(\xi (\tau ))\hfill \\ \hfill & +{\widehat{w}}^{\text{T}}\nabla S(\xi )f(\xi )-\frac{1}{4}{\widehat{w}}^{\text{T}}\nabla S(\xi )D\nabla S{(\xi )}^{\text{T}}\widehat{w}]\hfill \end{array}$

44

where ${\alpha }_H$ and ${\alpha }_c$ are the basic learning rate of the standard gradient descent algorithm and the learning rate of the stabilizing term, respectively. h is defined as follows

$h=\left\{\begin{array}{l}0,\hspace{1em}\text{if}{(\nabla J_s(\xi ))}^{\text{T}}(f(\xi )+g(\xi )\widehat{\mu })<0\hfill \\ 1,\hspace{1em}\text{else}\hfill \end{array}\right.$

45

where $J_s(\xi )$ is selected as a Lyapunov function candidate which is continuously differentiable. And assume that a positive definite matrix N exists, then the following equation is satisfied

$\begin{array}{ll}{\dot{J}}_s(\xi )\hfill & =(\nabla J_s(\xi {))}^{\text{T}}(f(\xi )+g(\xi ){\mu }^{*})\hfill \\ \hfill & =-{(\nabla J_s(\xi ))}^{\text{T}}N\nabla J_s(\xi )<0\hfill \end{array}$

46

It should be noted that $J_s(\xi )$ is a polynomial with the state variable and can be chosen appropriately, such as the form $J_s(\xi )=\frac{1}{2}{\xi }^{\text{T}}\xi$ .

Stability analysis

In this subsection, we will analyse the stability of the system and give the detailed proof that the approximate error $\stackrel{\sim }{w}$ of the NN weight and the state ξ are convergent.

Theorem 1

Consider the robot system (7) with approximate optimal control (38) and the NN weight updating law (43), then it is concluded that the approximate error $\stackrel{\sim }{w}$ of the NN weight and the state ξ are convergent.

Proof

See the Appendix.

Numerical simulation

Simulation settings

A two-degree-of-freedom (2-DOF) planar manipulator is adopted to verify the proposed control scheme. It is constructed by the robotics toolbox with parameters shown in Table 1. ³⁸ The numerical simulation shown in Figure 2 runs on the MATLAB 2018a software where an ode3 solver is chosen with a fixed time step of 0.01 s, simulation time 20 s and other settings remain default. The initial joint position is $q_0{=[0.08211,1.897]}^{\text{T}}$ and the user-defined trajectory is $x_d(t)=[0.3\text{exp}(-t{),0.5]}^{\text{T}}$ . The environment dynamics is simulated as

$-F=C_E\dot{x}+G_E(x-x_0)$

47

where C_E , G_E and x ₀ are chosen as $\text{diag}(0.1,0.1)$ , $\text{diag}(1.0,1.0)$ , 0.2, respectively, which are unknown during the simulation. For simplicity and without losing generality, only the trajectory along the x-axis is modified and interfered with the external forces.

Table 1.

Parameters of the robot manipulator.

Parameters	Values
l 1	0.50 m
$l_{c1}$	0.25 m
l 2	0.50 m
$l_{c2}$	0.25 m
m 1	5 kg
m 2	5 kg

Figure 2.

Settings of the numerical simulation.

For the proposed control scheme, parameters are set as below: to calculate the optimal trajectory in equation (13), $Q_{E1}=1.0$ , $R_E=1.0$ and $U_d=-0.3$ ; the feedback gain in the inverse kinematics in equation (22), $K_f=30$ and $\sigma =1e-6$ ; as for the ADP controller, in equations (4), (23), (38) and (43), $\Lambda =\text{diag}(5,5)$ , $R=\text{diag}(0.02,0.02)$ , $Q=\text{diag}(2.0,2.0)$ , ${\alpha }_H=0.5$ and ${\alpha }_c=2.5$ . Besides, an RBFNN is selected to approximate the cost function in equation (23), where $\widehat{J}={\widehat{w}}^{\text{T}}S(\xi )$ , $S_i(\xi )=\text{exp}(\parallel \xi -c_{\text{rbf}}\parallel /{\sigma }_{\text{rbf}}^2)$ with $\widehat{w}\in {ℝ}^9$ , $S(\xi )\in {ℝ}^9$ , $\widehat{w}(0)=0$ , ${\sigma }_{\text{rbf}}=0.55$ , $c_{\text{rbf}}\in [-0.2,0.0,0.2]\times [-0.2,0.0,0.2]$ .

Simulation results

In this subsection, two cases will be compared to demonstrate the validity of the proposed scheme. Note that, the environment dynamics of the simulation is not totally consistent with that in equation (9), and x ₀ is unknown. Therefore, two different K_e values are considered and examined. Case 1: the feedback gain $K_e^{\text{pro}}=[-0.5367,0.22840]$ acquired from the proposed scheme, which is different from Case 2: the ideal feedback gain $K_e^{\text{opt}}=[-0.4142,0.6604]$ obtained by calculating offline with the exact values of G_E and C_E (the unknown x ₀ is ignored in this case). For fair comparison, in Case 2, the trajectory will be modified at the time as Case 1.

Simulation results are shown in Figures 3 to 6. Figure 3 shows the modification process of the user-defined trajectory along the x-axis of both cases. It is not until around 4.1 s that the rank condition in equation (16) is satisfied following that the trajectory starts being modified. During the transient process, it can be found that the modified trajectory of Case 2 has a slight oscillation, and this subsequently triggers larger tracking errors compared with Case 1. The steady state and force pair of Case 1 and Case 2 trajectories at 10.28 s are 0.13 m/−0.07 N and 0.14 m/−0.06 N, respectively, which is in line with the time series of the cost function in equation (12) of both cases as shown in Figure 4. From the figure, we can see that after the modification of trajectory, the cost function of Case 1 is smaller than that of Case 2, which implies that in this simulation settings where the actual existence of unknown x ₀ cannot be neglected, the feedback gain obtained from the proposed scheme is more appropriate. Note that, due to the unknown x ₀, the environment dynamics in equation (9) used for the designing of the trajectory modifier differentiates from that in equation (47) used for simulation. Therefore, under this situation, actually neither the K_e of Case 1 nor Case 2 is the optimal one. However, the proposed method still works and regards the dynamics in equation (47) as a linear one with an appropriate feedback gain. This has demonstrated the effectiveness of the proposed admittance control method.

Figure 3.

Simulation results of trajectory modification. (a) Case 1, $K_e=K_e^{\text{pro}}$ and (b) Case 2, $K_e=K_e^{\text{opt}}$ .

Figure 4.

The time series of cost function in equation (12), $Q_{E1}=1.0$ and $R_E=1.0$ .

Figure 5.

Modified trajectories corresponding to the different choice of $Q_{E1}$ in equation (12).

Figure 6.

Simulation results of control performance. (a) Case 1, $K_e=K_e^{\text{pro}}$ and (b) Case 2, $K_e=K_e^{\text{opt}}$ .

Figures 5 to 7 are plotted for analysing the performance of the ADP-based controller. Figure 6 shows the control torques τ and sliding mode surface z of Case 1 and Case 2. On the whole, the proposed scheme tracks the both modified trajectory well, given that only nine neurons are used in the RBFNN, and the control torques are within the physical limitation. Besides, weights convergence can be observed in Figure 7. Note that, because of the introduced additional term $\nabla J_s$ , the initial admissible policy requirement is relaxed. Thus, in the simulation we choose the weights w to be zeros, without worrying about the control stability. This can be observed from Figure 6 that despite initial errors are large, they finally converge to zeros after some oscillations. Table 2 shows the feedback gain K_e calculated online using the proposed admittance control under the choices of different $Q_{E1}$ in equation (12). Its corresponding reference trajectories are shown in Figure 5 where the dashed lines denote the reference trajectories after modification and the solid lines stand for the actual trajectories of the robot end-effector under the control of the proposed ADP controller. Obviously, as the $Q_{E1}$ is selected larger, the reference trajectories tend to get closer to the user-desired trajectory, which is consistent with the fact that more cost is applied to the modified error X_e . Furthermore, although the reference trajectory varies, the proposed ADP controller is still eventually able to track the input signals with the same set of parameters. These also reflect the effectiveness of the proposed ADP controller and admittance control.

Figure 7.

Weights of the RBFNN. (a) Case 1, $K_e=K_e^{\text{pro}}$ and (b) Case 2, $K_e=K_e^{\text{opt}}$ .

Table 2.

Feedback gains of different $Q_{E1}$ .

	$Q_{E1}=0.8$ $R_E=1.0$	$Q_{E1}=1.0$ $R_E=1.0$	$Q_{E1}=1.2$ $R_E=1.0$	$Q_{E1}=1.4$ $R_E=1.0$	$Q_{E1}=1.6$ $R_E=1.0$	$Q_{E1}=1.8$ $R_E=1.0$	$Q_{E1}=2.0$ $R_E=1.0$
$K_{e1}$	−0.3625	−0.5367	−0.7205	−0.9055	−1.0870	−1.2620	−1.4300
$K_{e2}$	0.1090	0.2284	0.3606	0.4968	0.6317	0.7627	0.8884

Conclusion

The optimal control of robots interacting between unknown environment was studied in this article. An ADP-based controller with admittance adaptation was proposed. The unknown environment was regarded as a linear system and a compliant behaviour was guaranteed by the admittance adaptation control. In addition, NN was introduced into ADP controller to ensure trajectory tracking of the robot with minimal cost. The stability of the robot system was proved and simulation studies demonstrated the effectiveness of the proposed control scheme.

Because of the complexity of the robot system, dynamic uncertainties and input constraints such as saturation and dead zone are very common in robot systems, which will not only affect system performance but also may lead to system instability. ^24,39,40 Therefore, in the frame of ADP, the optimal control problem with dynamic uncertainties and input constraints will be considered in our future work.