Model-free optimal decentralized sliding mode control for modular and reconfigurable robots based on adaptive dynamic programming

Abstract

In this article, a model-free decentralized sliding mode control method is proposed based on adaptive dynamic programming algorithm to solve the problem of optimal trajectory tracking control of modular and reconfigurable robots. The dynamic model of modular and reconfigurable robot is formulated by a synthesis of joint subsystems with interconnected dynamic couplings. Based on sliding mode control technique, the optimal control problem of the modular and reconfigurable robot systems is transformed into an optimal compensation issue of unknown dynamics of each joint subsystems, in which the interconnected dynamic couplings effects among the subsystems are approximated by using the developed neural network identifier. Based on policy iteration scheme and the adaptive dynamic programming algorithm, the Hamilton–Jacobi–Bellman equation can be solved by using the critic neural network, so that optimal control policy can be obtained. The closed-loop system is proved to be asymptotically stable by using the Lyapunov theory. Finally, simulation results are provided to demonstrate the effectiveness of the method.

Introduction

Modular and reconfigurable robots (MRRs) ¹ have attracted wide attention in robotics community since they are possessed of better structural adaptability and flexibility than conventional robots. Until now, MRRs have wide applications in many extreme occasions, such as space explorations, disaster assistance, hazard survey, and medical assistance. Correspondingly, appropriate control systems are required to guarantee the accuracy and efficiency of MRRs in the face of different tasks.

As a useful tool to deal with disturbances, sliding mode control (SMC) technique may effectively improve the robustness of the nonlinear systems. There are some insightful papers that address the stabilization and tracking control problems by using SMC technique. Saleh and Fairouz ² investigated a robust adaptive second-order SMC method for tracking problem of a class of uncertain linear systems with matched and unmatched disturbances. Donya and Saleh ³ proposed an adaptive super-twisting decoupled terminal SMC technique for a class of fourth-order systems. Saleh and Fairouz ⁴ offered an adaptive global second-order sliding surface for perturbed dynamical systems with matched and unmatched external disturbances. Moreover, SMC has been widely used to design the controllers of manipulators. A sliding mode robust control method was presented for the pan-tilt joint manipulator. ⁵ A sliding mode adaptive neural network (NN) control was presented for the nonholonomic mobile manipulator. ⁶ Some investigations have presented to address the problems of manipulators,^7–9 and these methods are further implemented for controlling MRR systems. A stable adaptive fuzzy SMC method was proposed for an MRR to satisfy modular software. ¹⁰ Slotine and Sastry ¹¹ applied SMC technique on 2-degree-of-freedom (DOF) rigid MRRs to deal with the problem of tracking time-varying reference trajectories. Ficola and Cava ¹² presented an SMC method with two sliding surfaces, which gave the application on controlling two-joint MRR. However, the above-mentioned methods have not considered the problem of improving the efficiency of MRRs. Indeed, MRRs are always utilized in extreme occasions without external power supply; therefore, an ideal controller for MRRs should be possessed of the properties that guarantee the robustness of robotic systems and simultaneously consider the optimality of the composite of output power and error characteristics.

Optimal control, which was developed about five decades ago by Bellman ¹³ and Pontryagin, ¹⁴ and hitherto many practical and theoretical results have been represented.^15–18 For the perspectives of mathematics, optimal control problem is addressed by minimizing the desired performance index, and then, the solution can be acquired by addressing Hamilton–Jacobi–Bellman (HJB) equation approximatively as it gives the sufficient conditions for the optimality. For addressing the tracking problem of wheeled mobile robot, Lewis and Syrmos ¹⁸ investigated a near tracking control method based on receding-horizon dual-heuristic programming algorithm. Based on reinforcement learning theory, Bhsin et al. ¹⁹ addressed the optimal coordination control problem of multiple robots when dealing with targets with desired trajectories. Nageshrao et al. ²⁰ proposed an optimal passive control for the 2-DOF manipulator by using the energy balance theory. Tang et al. ²¹ proposed a learning-based adaptive optimal control method, which is used to solve the tracking problem of n-link robots. The mentioned methods above are all belong to the centralized control scheme. Indeed, an important property of MRRs is that their modules can be replaced, removed, appended, and optionally without adjusting control parameters of others, so that there exists physical restrictions on information interchange among the joint modules of the robotic system. Unfortunately, this kind of restrictions makes it impossible to adopt a centralized control method for MRR systems. To deal with the drawbacks of the conventional centralized control schemes, Li et al. ²² presented a decentralized robust control method for MRRs based on self-tuned feedback gain. In our previous works, we also paid attention to investigate decentralized robust control, ²³ decentralized trajectory tracking control, ²⁴ decentralized fault-tolerant control, ²⁵ and decentralized force/position control ²⁶ for MRR systems. However, the mentioned control methods for MRRs are not considered the optimal implementation of the controllers, which can guarantee the stability of robotic systems and simultaneously ensure the optimality of composite of error characteristics and output energy efficiency. Some researchers are presented by combining adaptive dynamic programming (ADP)-based optimal method with decentralized control scheme. Based on ADP theory, Bian et al. ²⁷ presented a decentralized optimal controller based on unmatched uncertainties. Zhao et al. ²⁸ proposed a proportional–integral (PI) algorithm–based decentralized scheme for large systems with mismatched interconnections. Dong et al. ²⁹ addressed the optimal control problems of MRRs by combining the model-based compensation control with ADP-based learning control, and their researches are further expanded to deal with the optimal tracking control issues of MRRs with uncertain environments. ³⁰ The ADP-based decentralized methods are proposed to solve the stabilization control problem of complex robot manipulators and nonlinear systems. However, to the best of authors’ knowledge, there are few researches concentrated on dealing with model-free optimal decentralized SMC of manipulator systems, especially for MRR systems.

In this article, a model-free decentralized SMC method is presented for MRR system via PI scheme and ADP methods. First, the dynamic formulation of MRRs is composed as a synthesis of joint subsystems with interconnected dynamic coupling (IDC) effects. Then, based on SMC technique, the optimal control problem has been transformed into an optimal compensation issue of subsystem unknown dynamics. A decentralized control strategy is designed, which uses only local dynamic information of each joint module, in which subsystem dynamic model is completely unknown. Based on the ADP method and the PI algorithm, the HJB equation is solved by using a critic NN, and then, the optimal control policy can be derived. According to the Lyapunov theory, the closed-loop MRR systems are proved to be asymptotically stable. Finally, simulations are represented to verify the advantage and effectiveness of the proposed method.

The main contributions of this article can be summarized as follows:

• To the best of our knowledge, it is the first time to extend the ADP approach to address the model-free decentralized optimal control problem of MRR systems. The proposed scheme can be utilized to MRRs with different configurations and different environment without changing control parameters. Unlike the conventional ADP methods that use action NN and critic NN, in this research, the optimal control policy is obtained by using only critic NN, and the training of action NN is no longer needed. It infers that the computational burden can be reduced effectively.

• Unlike the existing methods that consider the IDC effects, a kind of system disturbance with known upper bounds, in this article, the IDC effects, which are with larger order of magnitudes than some other system dynamics, are addressed independently and based on target by using the developed NN identifier-based compensation control law.

Problem statement

The dynamic formulation of robot system with n-DOF can be formulated as follows

$$M(q)\stackrel{··}{q}+C(q,\stackrel{·}{q})\stackrel{·}{q}+G(q)=u$$

(1)

where $u\in R^n$ represents joint torque, $q\in R^n$ means vector of joint positions, $M(q)\in R^{n\times n}$ indicates inertia matrix, $C(q,\stackrel{·}{q})\in R^{n\times n}$ denotes centripetal force and Coriolis term, and $G(q)\in R^n$ represents gravity term.

In the practical application, such as space exploration or disaster rescue, MRR consists of many joint modules, which brings the problem of heavy computational burden and complex control structure. In order to address this drawback, we consider each robotic joint as a subsystem of whole MRR system, which contains IDCs among each subsystem. The ith subsystem dynamic model is expressed as

$$\begin{array}{c}{M}_i(q_i){\stackrel{··}{q}}_i+C_i(q_i,{\stackrel{·}{q}}_i){\stackrel{·}{q}}_i+G_i(q_i)+I_i(q,\stackrel{·}{q},\stackrel{··}{q})=u_i\hfill \\ I_i(q,\stackrel{·}{q},\stackrel{··}{q})=\{\sum _{j=1,j\ne i}^n[M_{\mathrm{ii}}(q)-M_i(q_i)]{\stackrel{··}{q}}_i+M_{\mathrm{ij}}(q){\stackrel{··}{q}}_j\}\hfill \\ +\{\sum _{j=1,j\ne i}^n[-C_i(q_i,{\stackrel{·}{q}}_i)+C_{\mathrm{ii}}(q,\stackrel{·}{q})]{\stackrel{·}{q}}_i+C_{\mathrm{ij}}(q,\stackrel{·}{q}){\stackrel{·}{q}}_j\}\hfill \\ +{\overline{G}}_i(q)-G_i(q_i)\hfill \end{array}$$

(2)

where $q_i,{\stackrel{·}{q}}_i,{\stackrel{··}{q}}_i,u_i,{\overline{G}}_i$ are the ith elements of the vectors $q,\stackrel{·}{q},\stackrel{··}{q},u,G$ . $M_{\mathrm{ij}}$ and $C_{\mathrm{ij}}$ are the ijth element of the matrices $M$ and $C$ , respectively. $I_i$ is the IDC term.

Let $x_i=[x_{i1},x_{i2}{]}^T=[q_i,{\stackrel{·}{q}}_i{]}^T$ , the dynamic model can be formulated as the state space (3)

$$\{\begin{array}{c}{\stackrel{·}{x}}_{i1}=x_{i2}\hfill \\ \multicolumn{1}{c}{{\stackrel{·}{x}}_{i2}=f_i(x_i)+g_i(x_i)u_i+{\psi }_i(x)}\\ \multicolumn{1}{c}{y_i=x_i}\end{array}$$

(3)

where $x_i$ is the state of the ith subsystem

$$\begin{array}{c}{f}_i(x_i)=M_i^{-1}(q_i)[-C_i(q_i,{\stackrel{·}{q}}_i){\stackrel{·}{q}}_i-G_i(q_i)]\hfill \\ g_i(x_i)=M_i^{-1}(q_i),{\psi }_i(x)=-M_i^{-1}(q_i)I_i(q,\stackrel{·}{q},\stackrel{··}{q})\hfill \end{array}$$

where ${\psi }_i(x)$ is the IDC term and $x$ is the state vector with regard to the whole MRRs.

Assumption 1

The desired position $q_{\mathrm{id}}$ , velocity ${\stackrel{·}{q}}_{\mathrm{id}}$ , and acceleration ${\stackrel{··}{q}}_{\mathrm{id}}$ are bounded.

For getting rid of the norm-boundedness assumption of IDC, desired states of coupled subsystems are used to instead of actual ones. IDC term can be written as

$${\psi }_i(x)=\mathrm{\Delta }{\psi }_i(x,x_{\mathrm{jd}})+{\psi }_i(x_i,x_{\mathrm{jd}})$$

where $x_{\mathrm{jd}}$ denotes the desired states of the coupled subsystems with $j=1,\dots ,N(j\ne i),{\psi }_i(x_i,x_{\mathrm{jd}})$ representing certain part of IDC, where $x_{\mathrm{jd}}$ is used instead of actual one $x_j$ , and $\mathrm{\Delta }{\psi }_i(x,x_{\mathrm{jd}})$ represents the substitution error. Equation (3) can be rewritten as follows

$$\{\begin{array}{c}{\stackrel{·}{x}}_{i1}=x_{i2}\hfill \\ \multicolumn{1}{c}{{\stackrel{·}{x}}_{i2}=F_i(x_i,x_{\mathrm{jd}})+g_i(x_i)u_i(x_i)+\mathrm{\Delta }{\psi }_i(x,x_{\mathrm{jd}})}\\ \multicolumn{1}{c}{y_i=x_i}\end{array}$$

(4)

where $F_i(x_i,x_{\mathrm{jd}})={\psi }_i(x_i,x_{\mathrm{jd}})+f_i(x_i)$ .

Remark 1

The system dynamic terms ${\psi }_i(x),F_i(x_i,x_{\mathrm{jd}}),g_i(x_i)$ are always unknown, since the configuration of MRRs may change according to different tasks. This is the reason why we have to design the model-free controller for MRRs. Besides, an appropriate model-free controller should be designed with the property that fit for different configurations without adjusting control parameters.

Accordingly, in the following section, a model-free decentralized optimal SMC method is presented for MRRs to ensure that the closed-loop systems are asymptotically stable.

Model-free decentralized optimal SMC based on ADP

Derivation of the optimal SMC scheme

Define the joint position tracking error as

$$e_i=x_{i1}-x_{i1d}$$

(5)

Then, the time derivative of equation (5) can be obtained as

$${\stackrel{·}{e}}_i={\stackrel{·}{x}}_{i1}-{\stackrel{·}{x}}_{i1d}$$

(6)

According to the frame of the SMC method, we can define the sliding surface as follows

$$s_i(t)={\stackrel{·}{e}}_i(t)+k_i{e}_i(t)$$

(7)

where $k_i>0$ denotes the determined constant. Then, the time derivative of equation (7) is given as

$${\stackrel{·}{s}}_i(t)={\stackrel{··}{e}}_i(t)+k_i{\stackrel{·}{e}}_i(t)=F_i+g_i{u}_i-{\stackrel{··}{x}}_{i1d}+k_i{\stackrel{·}{e}}_i+\mathrm{\Delta }{\psi }_i$$

(8)

According to equations (7) and (8), the objective of the SMC is to satisfy the relation $s_i(t)={\stackrel{·}{s}}_i(t)=0$ . Moreover, for MRR systems, an ideal sliding mode controller should not only be possessed of the properties that guarantee the convergence of error performance in the sliding surface but also consider the optimality of the composite of error characteristics and output power simultaneously.

Therefore, the control objective can be transformed to design the optimal SMC law $u_i^{*}$ . Accordingly, one can design the following performance index function

$$\underset{u_i\in {\mathrm{\Psi }}_i({\mathrm{\Phi }}_i)}{J_i(s_i)}=\underset{0}{\overset{\infty }{\int }}({\mathrm{\Lambda }}_i{(\nabla J_i^{*}(s_i))}^2+H_i^T{H}_i+U_i(s_i,u_i))\mathrm{dt}$$

(9)

where $U_i(s_i,u_i(s_i))=s_i^T{Q}_i{s}_i+u_i^T{R}_i{u}_i$ is the utility function and $U_i(0)=0,U_i(s_i,u_i(s_i))\ge 0$ . $Q_i,R_i$ are the positive defined matrices, $H_i$ indicates the known up-bound function for $\mathrm{\Delta }{\psi }_i$ and one obtains $\mathrm{\Delta }{\psi }_i^T\mathrm{\Delta }{\psi }_i\le H_i^T{H}_i$ , and ${\mathrm{\Lambda }}_i$ is a positive constant.

It is noted that the decentralized control $u_i$ should not only stabilize the robot system but also guarantee that equation (9) is finite, that is, the control policy must be admissible. Here, we introduce the definition of admissible control.

Definition 1

In equation (4), decentralized control $u_i(s_i)$ is called admissible for the cost function (9) on the compact set ${\mathrm{\Phi }}_i$ , and $J_i(s_i)$ is finite, $\forall s_i\in {\mathrm{\Phi }}_i$ . Give an admissible control policy $u_i\in {\Xi }_i({\mathrm{\Phi }}_i)$ , where $u_i\in {\Xi }_i({\mathrm{\Phi }}_i)$ denotes admissible control set, the infinite performance index function (9) is a Lyapunov equation

$$\begin{array}{c}\multicolumn{1}{c}{0=H_i^T{H}_i+{\mathrm{\Lambda }}_i{(\nabla J_i^{*}(s_i))}^2+U_i(s_i,u_i(s_i))+{(\nabla J_i(s_i))}^T}\\ \multicolumn{1}{c}{·(\mathrm{\Delta }{\psi }_i+g_i(x_i)u_i(x_i)+F_i(x_i,x_{\mathrm{jd}})-{\stackrel{··}{x}}_{i1d}+k_i{\stackrel{·}{e}}_i(t))}\end{array}$$

(10)

where $J_i^{*}(s_i)$ is the optimal performance index function and $\nabla J_i(s_i)=\partial J_i(s_i)/\partial (s_i),J_i(0)=0$ .

Define the Hamilton function and the optimal performance index function as follows

$$\begin{array}{c}{H}_i(s_i,u_i(s_i),\nabla J_i(s_i))=H_i^T{H}_i+s_i^T{Q}_i{s}_i+{u_i}^T{R}_i{u}_i\hfill \\ +{(\nabla J_i(s_i))}^T(F_i(x_i,x_{\mathrm{jd}})+g_i(x_i)u_i-{\stackrel{··}{x}}_{i1d}+k_i{\stackrel{·}{e}}_i+\mathrm{\Delta }{\psi }_i)\hfill \\ +{\mathrm{\Lambda }}_i{(\nabla J_i^{*}(s_i))}^2\hfill \end{array}$$

(11)

$$J_i^{*}(s_i)=\underset{u_i\in (s_i)}{min}\underset{0}{\overset{\infty }{\int }}\left(\begin{array}{c}{U}_i(s_i,u_i(s_i))+H_i^T{H}_i\hfill \\ \multicolumn{1}{c}{+{\mathrm{\Lambda }}_i{(\nabla J_i^{*}(s_i))}^2}\end{array}\right)\mathrm{dt}$$

(12)

Under the framework of optimal control design, one obtains that $J_i^{*}(s_i)$ satisfies the following HJB equation

$$0=\underset{u_i\in (s_i)}{min}{H}_i(s_i,u_i(s_i),\nabla {J_i}^{*}(s_i))$$

(13)

If the ${J_i}^{*}(s_i)$ is continuously differentiable, the optimal SMC can be formulated as

$$u_i^{*}=-\frac{1}{2}{R}_i^{-1}{g}_i^T(x_i)\nabla J_i^{*}(s_i)$$

(14)

Rewriting $u_i^{*}$ as the form of $u_i^{*}=u_{i1}+u_{i2}^{*}$ , $u_{i1}$ is used to compensate the IDC and $u_{i2}^{*}$ is proposed to address the optimal compensation problem of unknown dynamics of the ith subsystem.

Combining equations (13) and (10), we get

$$\begin{array}{c}{(\nabla J_i(s_i))}^T(F_i(x_i,x_{\mathrm{jd}})+\mathrm{\Delta }{\psi }_i-{\stackrel{··}{x}}_{i1d}+k_i{\stackrel{·}{e}}_i+g_i(x_i)u_i)\hfill \\ \multicolumn{1}{c}{=-s_i^T{Q}_i{s}_i-{(u_{i1}+u_{i2}^{*})}^T{R}_i(u_{i1}+u_{i2}^{*})-{\mathrm{\Lambda }}_i{(\nabla J_i^{*}(s_i))}^2}\\ \multicolumn{1}{c}{-H_i^T{H}_i}\end{array}$$

(15)

Next, an identifier-based controller is used to compensate IDCs.

Identification of the IDC

In this section, an identifier is presented to approximate the term $\mathrm{\Delta }{\psi }_i(x,x_{\mathrm{jd}})$ .

Assumption 2 ³¹

NN activation function $\sigma (·)$ and its derivative $\sigma \prime (·)$ are bounded.

Assumption 3 ³²

The NN approximation error has been upper bounded by the unknown constant.

According to the above assumptions, $\mathrm{\Delta }{\psi }_i$ can be denoted by the single-layer NN, which is given as

$$\mathrm{\Delta }{\psi }_i=w_{\mathrm{ih}}^T{\sigma }_{\mathrm{ih}}(x_{\mathrm{ih}},x_D)+{\epsilon }_{\mathrm{ih}}(x_{\mathrm{ih}})$$

(16)

where ${\sigma }_{\mathrm{ih}}(x_{\mathrm{ih}},x_D)$ represents NN activation function, $x_{\mathrm{ih}}$ indicates the determined NN state, $w_{\mathrm{ih}}$ is an unknown ideal NN weight, $x_D=[x_{1d},x_{2d},\dots ,x_{\mathrm{md}}{]}^T,m<i$ means the reference state vector which are known and bounded, and ${\epsilon }_{\mathrm{ih}}$ means the NN approximation error. Considering nonlinear dynamic system, the bounded input $u_{\mathrm{ih}}$ is represented as

$${\stackrel{·}{x}}_{\mathrm{ih}}=\mathrm{\Delta }{\psi }_i+u_{\mathrm{ih}}=w_{\mathrm{ih}}^T{\sigma }_{\mathrm{ih}}(x_{\mathrm{ih}},x_D)+{\epsilon }_{\mathrm{ih}}(x_{\mathrm{ih}})+u_{\mathrm{ih}}$$

(17)

The NN identifier, which is used to approximate equation (17), is as follows

$${\stackrel{·}{\hat{x}}}_{\mathrm{ih}}=u_{\mathrm{ih}}+\mathrm{\Delta }{\hat{\psi }}_i={\hat{w}}_{\mathrm{ih}}^T{\sigma }_{\mathrm{ih}}+{\epsilon }_{\mathrm{ih}}+u_{\mathrm{ih}}+r_{\mathrm{ih}}$$

(18)

where ${\hat{x}}_{\mathrm{ih}}$ is the identification state of $x_{\mathrm{ih}}$ , $\mathrm{\Delta }{\hat{\psi }}_i$ is the approximation of $\mathrm{\Delta }{\psi }_i$ , ${\hat{w}}_{\mathrm{ih}}$ represents the weight estimate, and $r_{\mathrm{ih}}$ represents an integral of the sign of the error feedback term

$$r_{\mathrm{ih}}=k_{\mathrm{ih}}{e}_{\mathrm{ih}}+v_{\mathrm{ih}}$$

(19)

where $e_{\mathrm{ih}}=-{\hat{x}}_{\mathrm{ih}}+x_{\mathrm{ih}}$ represents the identification error. ${\zeta }_{\mathrm{ih}}$ means the generalized solution to

$${\stackrel{·}{\zeta }}_{\mathrm{ih}}=(k_{\mathrm{ih}}{\alpha }_{\mathrm{ih}}+{\gamma }_{\mathrm{ih}})e_{\mathrm{ih}}+{\delta }_{i1}\mathrm{sgn}(e_{\mathrm{ih}})$$

(20)

where ${\alpha }_{\mathrm{ih}},k_{\mathrm{ih}},{\delta }_{i1},{\gamma }_{\mathrm{ih}}$ indicate determined positive constants and $\mathrm{sgn}(•)$ means the signal function. Combining equation (17) with equation (18), one can obtain the identification error dynamics as follows

$${\stackrel{·}{e}}_{\mathrm{ih}}=\mathrm{\Delta }{\tilde{\psi }}_i=w_{\mathrm{ih}}^T{\sigma }_{\mathrm{ih}}-{\hat{w}}_{\mathrm{ih}}^T{\sigma }_{\mathrm{ih}}+{\epsilon }_{\mathrm{ih}}-r_{\mathrm{ih}}$$

(21)

where $\mathrm{\Delta }{\tilde{\psi }}_i=\mathrm{\Delta }{\psi }_i-\mathrm{\Delta }{\hat{\psi }}_i$ . The identification error function can be defined as follows

$${\rho }_{\mathrm{ih}}={\stackrel{·}{e}}_{\mathrm{ih}}+{\alpha }_{\mathrm{ih}}{e}_{\mathrm{ih}}$$

(22)

By using equations (21) and (22), one obtains that

$$\begin{array}{cc}\multicolumn{1}{c}{{\stackrel{·}{\rho }}_{\mathrm{ih}}=}& w_{\mathrm{ih}}^T\sigma {\prime }_{\mathrm{ih}}{\stackrel{·}{x}}_{\mathrm{ih}}-{\stackrel{·}{\hat{w}}}_{\mathrm{ih}}^T{\sigma }_{\mathrm{ih}}-{\hat{w}}_{\mathrm{ih}}^T\sigma {\prime }_{\mathrm{ih}}{\stackrel{·}{\hat{x}}}_{\mathrm{ih}}+{\stackrel{·}{\epsilon }}_{\mathrm{ih}}\hfill \\ \multicolumn{1}{c}{}& -k_{\mathrm{ih}}{s}_{\mathrm{ih}}-{\gamma }_{\mathrm{ih}}{e}_{\mathrm{ih}}-{\delta }_{i1}\mathrm{sgn}(e_{\mathrm{ih}})+{\alpha }_{\mathrm{ih}}{\stackrel{·}{e}}_{\mathrm{ih}}\hfill \end{array}$$

(23)

where the weight update law of equation (24) is given as follows

$${\stackrel{·}{\hat{w}}}_{\mathrm{ih}}=\mathrm{proj}\left({\mathrm{\Gamma }}_{\mathrm{ih}}\sigma {\prime }_{\mathrm{ih}}{\stackrel{·}{\hat{x}}}_{\mathrm{ih}}{\tilde{x}}_{\mathrm{ih}}^T\right)$$

(24)

where $\mathrm{proj}(•)$ is the projection operation and ${\mathrm{\Gamma }}_{\mathrm{ih}}>0$ represents a constant matric. Then, equation (23) can be written as

$$\begin{array}{cc}\multicolumn{1}{c}{{\stackrel{·}{\rho }}_{\mathrm{ih}}=}& {\tilde{P}}_{\mathrm{ih}1}(w_{\mathrm{ih}},{\hat{w}}_{\mathrm{ih}},e_{\mathrm{ih}})+P_{\mathrm{ih}2}(w_{\mathrm{ih}},{\hat{w}}_{\mathrm{ih}},x_{\mathrm{ih}})\hfill \\ \multicolumn{1}{c}{}& +P_{\mathrm{ih}3}({\hat{x}}_{\mathrm{ih}},{\tilde{w}}_{\mathrm{ih}})-{\gamma }_{\mathrm{ih}}{e}_{\mathrm{ih}}-{\delta }_{i1}\mathrm{sgn}(e_{\mathrm{ih}})-k_{\mathrm{ih}}{s}_{\mathrm{ih}}\hfill \end{array}$$

(25)

where

$$\begin{array}{c}{\tilde{P}}_{\mathrm{ih}1}={\alpha }_{\mathrm{ih}}{\stackrel{·}{e}}_{\mathrm{ih}}-{\stackrel{·}{\hat{w}}}_{\mathrm{ih}}^T{\sigma }_{\mathrm{ih}}+\frac{1}{2}{w}_{\mathrm{ih}}^T\sigma {\prime }_{\mathrm{ih}}{\stackrel{·}{e}}_{\mathrm{ih}}+\frac{1}{2}{{\hat{w}}^T}_{\mathrm{ih}}\sigma {\prime }_{\mathrm{ih}}{\stackrel{·}{e}}_{\mathrm{ih}}\hfill \\ P_{\mathrm{ih}2}=\frac{1}{2}{w}_{\mathrm{ih}}^T\sigma {\prime }_{\mathrm{ih}}{\stackrel{·}{x}}_{\mathrm{ih}}-\frac{1}{2}{{\hat{w}}^T}_{\mathrm{ih}}\sigma {\prime }_{\mathrm{ih}}{\stackrel{·}{x}}_{\mathrm{ih}}+{\stackrel{·}{\epsilon }}_{\mathrm{ih}}\hfill \\ P_{\mathrm{ih}3}=\frac{1}{2}{\tilde{w}}_{\mathrm{ih}}^T\sigma {\prime }_{\mathrm{ih}}{\stackrel{·}{\hat{x}}}_{\mathrm{ih}}\hfill \end{array}$$

(26)

in which ${\tilde{w}}_{\mathrm{ih}}=w_{\mathrm{ih}}-{\hat{w}}_{\mathrm{ih}}$ . Then, one can obtain the auxiliary term $P_{\mathrm{ih}3}(x_{\mathrm{ih}},{\tilde{w}}_{\mathrm{ih}})$ instead of ${\stackrel{·}{\hat{x}}}_{\mathrm{ih}}$ in ${\hat{P}}_{\mathrm{ih}3}$ by ${\stackrel{·}{x}}_{\mathrm{ih}}$ to facilitate stability analysis, and we obtain ${\tilde{P}}_{\mathrm{ih}3}={\hat{P}}_{\mathrm{ih}3}-P_{\mathrm{ih}3}$ and denote $P_{\mathrm{ih}}=P_{\mathrm{ih}2}+P_{\mathrm{ih}3}$ . According to Assumptions 2 and 3 and using equations (22), (24), and (26), the following bounds can be obtained

$$\begin{array}{c}\Vert {\tilde{P}}_{\mathrm{ih}1}\Vert \le {\mu }_{i1}(\Vert b_i\Vert )\Vert b_i\Vert ,\Vert P_{\mathrm{ih}2}\Vert \le {\xi }_{i1},\Vert P_{\mathrm{ih}3}\Vert \le {\xi }_{i2}\hfill \\ \Vert P_{\mathrm{ih}}\Vert \le {\xi }_{i3}+{\xi }_{i4}{\mu }_{i2}(\Vert b_i\Vert )\Vert b_i\Vert \hfill \\ \Vert {\stackrel{·}{e}}_{\mathrm{ih}}^T{P}_{\mathrm{ih}3}\Vert \le {\xi }_{i5}{\Vert e_{\mathrm{ih}}\Vert }^2+{\xi }_{i6}{\Vert s_{\mathrm{ih}}\Vert }^2\hfill \end{array}$$

(27)

where $b_i={\left[\begin{array}{cc}{e}_{\mathrm{ih}}^T\hfill & {\rho }_{\mathrm{ih}}^T\hfill \end{array}\right]}^T,{\mu }_{i1}(·),{\mu }_{i2}(·)$ are positive, non-decreasing, invertible functions. ${\xi }_{\mathrm{im}}>0,m=1,\dots ,6$ are constants. In order to analyze stability, denote

$$Q_{\mathrm{ih}}={\left[\begin{array}{cccc}{e}_{\mathrm{ih}}^T\hfill & {\rho }_{\mathrm{ih}}^T\hfill & \sqrt{M_i}\hfill & \sqrt{N_i}\hfill \end{array}\right]}^T$$

(28)

where $M_i$ represents the generalized solution to the differential equation

$$\begin{array}{cc}\multicolumn{1}{c}{{\stackrel{·}{M}}_i=}& -{\rho }_{\mathrm{ih}}^T(P_{\mathrm{ih}2}-{\delta }_{i1}\mathrm{sgn}\left(e_{\mathrm{ih}}\right))+{\stackrel{·}{e}}_{\mathrm{ih}}^T{P}_{\mathrm{ih}3}\hfill \\ \multicolumn{1}{c}{}& -{\delta }_{i2}{\mu }_{i2}(\Vert b_i\Vert )\Vert b_i\Vert \Vert e_{\mathrm{ih}}\Vert \hfill \end{array}$$

(29)

where ${\delta }_{i1},{\delta }_{i2}$ are chosen by following conditions to guarantee that $M_i\ge 0$

$${\delta }_{i1}>max({\xi }_{i2}+{\xi }_{i1},{\xi }_{i1}+\frac{{\xi }_{i3}}{{\alpha }_{\mathrm{ih}}}),{\delta }_{i2}>{\xi }_{i4}$$

(30)

The auxiliary function $N_i$ can be defined as

$$N_i=\frac{1}{4}{\alpha }_{\mathrm{ih}}\left[\mathrm{tr}\left({\tilde{w}}_{\mathrm{ih}}^T{\mathrm{\Gamma }}_{\mathrm{ih}}^{-1}{\tilde{w}}_{\mathrm{ih}}\right)\right]$$

(31)

in which $\mathrm{tr}(•)$ indicates the trace of matrix.

Theorem 1

The IDC term is indicated in equation (16) and the dynamic system is developed in equation (17). The utilization of the NN identifier presented in equation (18) and the weight update law given in equation (24) can guarantee the asymptotic identification of the IDCs in the sense that

$$\underset{t\to \infty }{lim}\Vert \mathrm{\Delta }{\tilde{\psi }}_i\Vert =0$$

(32)

provided $k_{\mathrm{ih}},{\gamma }_{\mathrm{ih}}$ are satisfied by the following conditions

$$k_{\mathrm{ih}}>{\xi }_{i6},{\gamma }_{\mathrm{ih}}>\frac{{\xi }_{i5}}{{\alpha }_{\mathrm{ih}}}$$

(33)

$\mathrm{where}{k}_{\mathrm{ih}},{\gamma }_{\mathrm{ih}},{\alpha }_{\mathrm{ih}},{\xi }_{i5},{\xi }_{i6}$ are in equations (20) and (27).

Proof

Define the Lyapunov function candidate $V_{\mathrm{ih}}(Q_{\mathrm{ih}}):D_{\mathrm{ih}}\to R_{\mathrm{ih}}$ as follows

$$V_{\mathrm{ih}}(Q_{\mathrm{ih}})=\frac{1}{2}{\rho }_{\mathrm{ih}}^T{\rho }_{\mathrm{ih}}+\frac{1}{2}{\gamma }_{\mathrm{ih}}{e}_{\mathrm{ih}}^T{e}_{\mathrm{ih}}+M_i+N_i$$

(34)

which satisfies the following relations

$$U_1(Q_{\mathrm{ih}})\le V_{\mathrm{ih}}(Q_{\mathrm{ih}})\le U_2(Q_{\mathrm{ih}})$$

(35)

where $U_1(Q_{\mathrm{ih}}),U_2(Q_{\mathrm{ih}})$ are continuous positive definite functions

$$\begin{array}{c}{U}_1(Q_{\mathrm{ih}})=\frac{1}{2}min({\gamma }_{\mathrm{ih}},1){\Vert Q_{\mathrm{ih}}\Vert }^2\hfill \\ \multicolumn{1}{c}{U_2(Q_{\mathrm{ih}})=max({\gamma }_{\mathrm{ih}},1){\Vert Q_{\mathrm{ih}}\Vert }^2}\end{array}$$

(36)

Define the time derivative of equation (34) as follows

$$\begin{array}{cc}\multicolumn{1}{c}{{\stackrel{·}{V}}_{\mathrm{ih}}}& =\nabla V_{\mathrm{ih}}^T{\left[\begin{array}{cccc}{\stackrel{·}{\rho }}_{\mathrm{ih}}^T\hfill & {\stackrel{·}{e}}_{\mathrm{ih}}^T\hfill & \frac{1}{2}{M}_i^{-\frac{1}{2}}{\stackrel{·}{M}}_i\hfill & \frac{1}{2}{N}_i^{-\frac{1}{2}}{\stackrel{·}{N}}_i\hfill \end{array}\right]}^T\hfill \\ \multicolumn{1}{c}{}& =\left[\begin{array}{cccc}{\rho }_{\mathrm{ih}}^T\hfill & {\gamma }_{\mathrm{ih}}{e}_{\mathrm{ih}}^T\hfill & 2M_i^{\frac{1}{2}}\hfill & 2N_i^{\frac{1}{2}}\hfill \end{array}\right]\hfill \\ \multicolumn{1}{c}{}& {\left[\begin{array}{cccc}{\stackrel{·}{\rho }}_{\mathrm{ih}}^T\hfill & {\stackrel{·}{e}}_{\mathrm{ih}}^T\hfill & \frac{1}{2}{M}_i^{-\frac{1}{2}}{\stackrel{·}{M}}_i\hfill & \frac{1}{2}{N}_i^{-\frac{1}{2}}{\stackrel{·}{N}}_i\hfill \end{array}\right]}^T\hfill \\ \multicolumn{1}{c}{}& <({\tilde{P}}_{\mathrm{ih}1}+P_{\mathrm{ih}2}+{\hat{P}}_{\mathrm{ih}3}-k_{\mathrm{ih}}{\rho }_{\mathrm{ih}}-{\delta }_{i1}\mathrm{sgn}(e_{\mathrm{ih}})-{\gamma }_{\mathrm{ih}}{e}_{\mathrm{ih}})\hfill \\ \multicolumn{1}{c}{}& {\rho }_{\mathrm{ih}}^T+{\gamma }_{\mathrm{ih}}{e}_{\mathrm{ih}}^T({\rho }_{\mathrm{ih}}-{\alpha }_{\mathrm{ih}}{e}_{\mathrm{ih}})-{\rho }_{\mathrm{ih}}^T(P_{\mathrm{ih}2}-{\delta }_{i1}\mathrm{sgn}(e_{\mathrm{ih}}))-\hfill \\ \multicolumn{1}{c}{}& {\stackrel{·}{e}}_{\mathrm{ih}}^T{P}_{\mathrm{ih}3}+{\delta }_{i2}{\mu }_{i2}(\Vert b_i\Vert )\Vert b_i\Vert \Vert e_{\mathrm{ih}}\Vert -\frac{1}{2}{\alpha }_{\mathrm{ih}}[\mathrm{tr}({\tilde{w}}_{\mathrm{ih}}^T{\mathrm{\Gamma }}_{\mathrm{ih}}^{-1}{\stackrel{·}{\hat{w}}}_{\mathrm{ih}})]\hfill \\ \multicolumn{1}{c}{}& =-{\alpha }_{\mathrm{ih}}{\gamma }_{\mathrm{ih}}{e}_{\mathrm{ih}}^T{e}_{\mathrm{ih}}-k_{\mathrm{ih}}{\rho }_{\mathrm{ih}}^T{\rho }_{\mathrm{ih}}+\frac{1}{2}{\alpha }_{\mathrm{ih}}{e}_{\mathrm{ih}}^T{\tilde{w}}_{\mathrm{ih}}^T{\sigma }_{\mathrm{ih}}^{\prime }{\stackrel{·}{\hat{x}}}_{\mathrm{ih}}\hfill \\ \multicolumn{1}{c}{}& +{\delta }_{i2}{\mu }_{i2}(\Vert b_i\Vert )\Vert b_i\Vert \Vert e_{\mathrm{ih}}\Vert -\frac{1}{2}{\alpha }_{\mathrm{ih}}\left[\mathrm{tr}\left({\tilde{w}}_{\mathrm{ih}}^T{\sigma }_{\mathrm{ih}}^{\prime }{\stackrel{·}{\hat{x}}}_{\mathrm{ih}}{e}_{\mathrm{ih}}^T\right)\right]\hfill \\ \multicolumn{1}{c}{}& +{\stackrel{·}{e}}_{\mathrm{ih}}^T({\hat{P}}_{\mathrm{ih}3}-P_{\mathrm{ih}3})+\frac{1}{2}{\alpha }_{\mathrm{ih}}{e}_{\mathrm{ih}}^T{\hat{w}}_{\mathrm{ih}}^T{\sigma }_{\mathrm{ih}}^{\prime }{\stackrel{·}{\hat{x}}}_{\mathrm{ih}}+{\rho }_{\mathrm{ih}}^T{\tilde{P}}_{\mathrm{ih}1}\hfill \end{array}$$

(37)

Canceling the common terms in equation (37), denoted as $k_{\mathrm{ih}}=k_{\mathrm{ih}1}+k_{\mathrm{ih}2},{\gamma }_{\mathrm{ih}}={\gamma }_{\mathrm{ih}1}+{\gamma }_{\mathrm{ih}2}$ , using equation (27), and completing the squares, one can obtain the upper bound of equation (37) as follows

$$\begin{array}{c}\multicolumn{1}{c}{{\stackrel{·}{V}}_{\mathrm{ih}}\le -({\alpha }_{\mathrm{ih}}{\gamma }_{\mathrm{ih}1}-{\xi }_{i5}){\Vert {\rho }_{\mathrm{ih}}\Vert }^2-(k_{i1}-{\xi }_{i6}){\Vert {\rho }_{\mathrm{ih}}\Vert }^2}\\ \multicolumn{1}{c}{+\frac{{\mu }_{i1}{(\Vert b_i\Vert )}^2}{4{\lambda }_{i2}}{\Vert b_i\Vert }^2+\frac{{\delta }_{i2}^2{\mu }_{i2}{(\Vert b_i\Vert )}^2}{4{\alpha }_{\mathrm{ih}}{\gamma }_{\mathrm{ih}2}}{\Vert b_i\Vert }^2}\end{array}$$

(38)

If equation (33) is satisfied, equation (38) can be written as

$$\begin{array}{c}{\stackrel{·}{V}}_{\mathrm{ih}}\le -{\lambda }_{i1}\Vert b_i{\Vert }^2+\frac{\mu {(\Vert b_i\Vert )}^2}{4{\lambda }_{i2}}\Vert b_i{\Vert }^2\le -U\left(Q_{\mathrm{ih}}\right)\hfill \\ \forall Q_i\in D_{\mathrm{ih}}\hfill \end{array}$$

(39)

where $U\left(Q_{\mathrm{ih}}\right)={\alpha }_{\mathrm{ic}}\Vert b_i{\Vert }^2$ , ${\lambda }_{i1}=min\{k_{i1}-{\xi }_{i6},{\alpha }_{\mathrm{ih}}{\gamma }_{\mathrm{ih}1}-{\xi }_{i5}\}$ , and ${\lambda }_{i2}=min\{k_{i2},({\alpha }_{\mathrm{ih}}{\gamma }_{\mathrm{ih}2}/{\delta }_{i2}^2)\}$ , in which a constant ${\alpha }_{\mathrm{ic}}>0$ is the positive semidefinite function defined as $D_{\mathrm{ih}}=\{Q_{\mathrm{ih}}|\Vert Q_{\mathrm{ih}}\Vert \le {\mu }_i^{-1}\left(2\sqrt{{\lambda }_{i1}{\lambda }_{i2}}\right)\}$ .

Let $R_{\mathrm{ih}}\subset D_{\mathrm{ih}}$ denote a set defined as

$$R_{\mathrm{ih}}=\{Q_{\mathrm{ih}}\subset D_{\mathrm{ih}}|U_2(Q_{\mathrm{ih}})<\frac{1}{2}{\left({\mu }_i^{-1}\left(2\sqrt{{\lambda }_{i1}{\lambda }_{i2}}\right)\right)}^2\}$$

(40)

The region of attraction in equation (40) can be adjusted arbitrarily big to contain any initial conditions. Then, we can get that

$${\alpha }_{\mathrm{ic}}\Vert b_i{\Vert }^2\to 0$$

(41)

while $\forall Q_{\mathrm{ih}}(0)\in R_{\mathrm{ih}},t\to \infty$ .

As per the definition of $b_i,{\rho }_{\mathrm{ih}},\mathrm{\Delta }{\psi }_i$ , one concludes that $\Vert e_{\mathrm{ih}}\Vert ,\Vert {\rho }_{\mathrm{ih}}\Vert ,\Vert {\stackrel{·}{e}}_{\mathrm{ih}}\Vert \to 0$ , while $\forall Q_{\mathrm{ih}}(0)\in R_{\mathrm{ih}},t\to \infty$ . According to equations (21) and (32), and using the conclusion above, we can get $\underset{t\to \infty }{lim}\Vert \mathrm{\Delta }{\tilde{\psi }}_i\Vert =0$ . The proof has been completed.

According to equations (18), (19), and (20), we can design $u_{i1}$ as

$$u_{i1}=-g_i(x_i)\left(\begin{array}{c}\underset{0}{\overset{t}{\int }}\left(\begin{array}{c}(k_{\mathrm{ih}}{\alpha }_{\mathrm{ih}}+{\gamma }_{\mathrm{ih}})e_{\mathrm{ih}}\hfill \\ \multicolumn{1}{c}{+{\delta }_{i1}\mathrm{sgn}(e_{\mathrm{ih}})}\end{array}\right)\mathrm{dt}\hfill \\ \multicolumn{1}{c}{+k_{\mathrm{ih}}{e}_{\mathrm{ih}}+{\hat{w}}_{\mathrm{ih}}^T{\sigma }_{\mathrm{ih}}}\end{array}\right)$$

(42)

where the weight ${\hat{w}}_{\mathrm{ih}}$ is updated by equation (24).

Critic NN implementation

For finding the optimal control of the MRR (equation (2)), we need to address HJB equation (13) for $\nabla J_i^{*}(s_i)$ and then obtain $u_i^{*}(s_i)$ . However, equation (13) is a nonlinear partial differential equation and difficult to be solved, because of the heavy burden of computation. Therefore, one employed a critic NN to approximate $J_i^{*}(s_i)$ .

The index function is highly non-analytic and nonlinear, the critic NN can be used to approximate $J_i(s_i)$ , and the proposed critic NN can be formulated as

$$J_i(s_i)={W_{\mathrm{ic}}}^T{\delta }_{\mathrm{ic}}(s_i)+{\epsilon }_{\mathrm{ic}}$$

(43)

where $W_{\mathrm{ic}}\in R^{{\overline{N}}_i}$ is the desired weight vector, ${\overline{N}}_i$ denotes the quantity of neurons in hidden layer, ${\delta }_{\mathrm{ic}}$ indicates the active function, and ${\epsilon }_{\mathrm{ic}}$ is the critic NN approximation error. Hence, the partial derivative of $\nabla J_i(s_i)$ can be represented as

$$\nabla J_i(s_i)={(\nabla {\delta }_{\mathrm{ic}}(s_i))}^T{W}_{\mathrm{ic}}+\nabla {\epsilon }_{\mathrm{ic}}$$

(44)

where $\nabla {\epsilon }_{\mathrm{ic}}$ and $\nabla {\delta }_{\mathrm{ic}}(s_i)=(\partial ({\delta }_{\mathrm{ic}}(s_i)))/\partial s_i\in R^{{\overline{N}}_i\times n}$ are partial derivative of the approximate error and activation function.

For system (2), combining equations (36) and (44), we get

$$\begin{array}{cc}\multicolumn{1}{c}{0=}& H_i^T{H}_i+{\mathrm{\Lambda }}_i{({(\nabla {\delta }_{\mathrm{ic}}(s_i))}^T{W}_{\mathrm{ic}}+\nabla {\epsilon }_{\mathrm{ic}})}^T\hfill \\ \multicolumn{1}{c}{}& \times ({(\nabla {\delta }_{\mathrm{ic}}(s_i))}^T{W}_{\mathrm{ic}}+\nabla {\epsilon }_{\mathrm{ic}})\hfill \\ \multicolumn{1}{c}{}& +{({(\nabla {\delta }_{\mathrm{ic}}(s_i))}^T{W}_{\mathrm{ic}}+\nabla {\epsilon }_{\mathrm{ic}})}^T{\stackrel{·}{s}}_i+U_i(s_i,u_i(s_i)).\hfill \end{array}$$

(45)

Substituting equation (44) into equation (13), the Hamilton function can be reformulated as

$$\begin{array}{c}{H}_i(s_i,u_i(s_i),W_{\mathrm{ic}})=H_i^T{H}_i+{\mathrm{\Lambda }}_i\hfill \\ ·{({(\nabla {\delta }_{\mathrm{ic}}(s_i))}^T{W}_{\mathrm{ic}}+\nabla {\epsilon }_{\mathrm{ic}})}^T\times ({(\nabla {\delta }_{\mathrm{ic}}(s_i))}^T{W}_{\mathrm{ic}}+\nabla {\epsilon }_{\mathrm{ic}})\hfill \\ +{({(\nabla {\delta }_{\mathrm{ic}}(s_i))}^T{W}_{\mathrm{ic}}+\nabla {\epsilon }_{\mathrm{ic}})}^T{\stackrel{·}{s}}_i+U_i(s_i,u_i(s_i))\hfill \\ =-\nabla {{\epsilon }_{\mathrm{ic}}}^T{\stackrel{·}{s}}_i=e_{\mathrm{iJh}},\hfill \end{array}$$

(46)

where $e_{\mathrm{iJh}}$ is the NN approximation remnants error NN. Because of the desired weight vector, $W_{\mathrm{ic}}$ is unknown, and critic NN can be approximated by

$${\hat{J}}_i(s_i)={{\hat{W}}_{\mathrm{ic}}}^T{\delta }_{\mathrm{ic}}(s_i),$$

(47)

where ${\hat{W}}_{\mathrm{ic}}$ is the estimation of $W_{\mathrm{ic}}$ and ${\hat{J}}_i(s_i)$ is the estimation of $J_i(s_i)$ .

The partial derivative of ${\hat{J}}_i(s_i)$ can be formulated as

$$\nabla {\hat{J}}_i(s_i)={({\delta }_{\mathrm{ic}}(s_i))}^T{\hat{W}}_{\mathrm{ic}}$$

(48)

Then, one can obtain the approximate Hamilton function as

$$\begin{array}{c}{H}_i(s_i,u_i(s_i),{\hat{W}}_{\mathrm{ic}})=H_i^T{H}_i+{\mathrm{\Lambda }}_i\hfill \\ {\left({(\nabla {\delta }_{\mathrm{ic}}(s_i))}^T{\hat{W}}_{\mathrm{ic}}\right)}^T\left({(\nabla {\delta }_{\mathrm{ic}}(s_i))}^T{\hat{W}}_{\mathrm{ic}}\right)\hfill \\ +{\left({(\nabla {\delta }_{\mathrm{ic}}(s_i))}^T{\hat{W}}_{\mathrm{ic}}\right)}^T{\stackrel{·}{s}}_i+U_i(s_i,u_i(s_i))={e_i}_J\hfill \end{array}$$

(49)

We use the objective function ${E_i}_J=(1/2)e_{\mathrm{iJ}}^T{e}_{\mathrm{iJ}}$ , which is minimized by gradient decent algorithm, to adjust the critic NN weight vector ${\hat{W}}_{\mathrm{ic}}$ that is updated by

$${\stackrel{·}{\hat{W}}}_{\mathrm{ic}}=-{\alpha }_i{e}_{\mathrm{iJ}}\nabla {\delta }_{\mathrm{ic}}(s_i){\stackrel{·}{s}}_i$$

(50)

where ${\alpha }_i$ is the updated rate of the critic NN. Denote ${\upsilon }_i=\nabla {\delta }_{\mathrm{ic}}(s_i){\stackrel{·}{s}}_i$ and consider that there exist a positive constant ${\upsilon }_{\mathrm{iL}}$ where $\Vert {\upsilon }_i\Vert \le {\upsilon }_{\mathrm{iL}}$ .

Define the weight approximation error as

$${\tilde{W}}_{\mathrm{ic}}=W_{\mathrm{ic}}-{\hat{W}}_{\mathrm{ic}}$$

(51)

Then, according to equations (46), (49), and (50), we conclude that

$$e_{\mathrm{iJ}}=e_{\mathrm{iJh}}-{{\tilde{W}}_{\mathrm{ic}}}^T{\upsilon }_i$$

(52)

The dynamics of the weight approximation error can be given as

$${\stackrel{·}{\stackrel{~}{W}}}_{\mathrm{ic}}=-{\stackrel{·}{\hat{W}}}_{\mathrm{ic}}={\alpha }_i{e}_{\mathrm{iJ}}{\upsilon }_i={\alpha }_i(e_{\mathrm{iJh}}-{{\tilde{W}}_{\mathrm{ic}}}^T{\upsilon }_i){\upsilon }_i$$

(53)

According to equations (40) and (44), the desired optimal control policy can be formulated as

$$u_{i2}^{*}(s_i)=-\frac{1}{2}{R}_i^{-1}{g}_i^T(x_i)({(\nabla {\delta }_{\mathrm{ic}}(s_i))}^T{W}_{\mathrm{ic}}+\nabla {\epsilon }_{\mathrm{ic}})$$

(54)

and it can be approximated as

$${\hat{u}}_{i2}^{*}(s_i)=-\frac{1}{2}{R}_i^{-1}{g}_i^T(x_i)\left({(\nabla {\delta }_{\mathrm{ic}}(s_i))}^T{\hat{W}}_{\mathrm{ic}}\right)$$

(55)

Note that the expression of equation (55) is obtained using the critic NN only and the training of action NN is no longer needed. It infers that the computational burden can be reduced effectively.

Theorem 2

Consider the cost function in equation (43), which is approximated by the single-layer critic NN, and the estimated cost function in equation (47) that is built by ${\hat{W}}_{\mathrm{ic}}$ . If the weight of the critic NN is updated by equation (53), then the weight approximation error in equation (53) can be proved to be uniformly ultimately bounded (UUB).

Proof

Choose the Lyapunov function candidate as

$$V_{\mathrm{ic}}(t)=\frac{1}{2{\alpha }_i}{\tilde{W}}_{\mathrm{ic}}^T{\tilde{W}}_{\mathrm{ic}}$$

(56)

The time derivative of $V_{\mathrm{ic}}(t)$ denotes

$$\begin{array}{cc}\multicolumn{1}{c}{{\stackrel{·}{V}}_{\mathrm{ic}}(t)}& =\frac{1}{{\alpha }_i}{\tilde{W}}_{\mathrm{ic}}^T{\stackrel{·}{\stackrel{~}{W}}}_{\mathrm{ic}}={\tilde{W}}_{\mathrm{ic}}^T(e_{\mathrm{iJh}}-{\tilde{W}}_{\mathrm{ic}}^T{\upsilon }_i){\upsilon }_i\hfill \\ \multicolumn{1}{c}{}& ={\tilde{W}}_{\mathrm{ic}}^T{e}_{\mathrm{iJh}}{\upsilon }_i-\Vert {\tilde{W}}_{\mathrm{ic}}^T{\upsilon }_i{\Vert }^2\le \frac{1}{2}\Vert e_{\mathrm{iJh}}{\Vert }^2-\frac{1}{2}\Vert {\tilde{W}}_{\mathrm{ic}}^T{\upsilon }_i{\Vert }^2\hfill \end{array}$$

Hence, we can get ${\stackrel{·}{V}}_{\mathrm{ic}}(t)\le 0$ , while $s_i$ lies outside the compact set, ${\Omega }_{i1}=\{{\tilde{W}}_{\mathrm{ic}}:\Vert {\tilde{W}}_{\mathrm{ic}}\Vert \le (e_{\mathrm{iJh}}/{\upsilon }_{\mathrm{iL}})\}$ . This completes the proof.

Combining equations (42) and (55), $u_i^{*}$ is given as

$$\begin{array}{c}{u}_i^{*}=-\frac{1}{2}{R}_i^{-1}{g}_i^T(x_i)\left({(\nabla {\delta }_{\mathrm{ic}}(s_i))}^T{\hat{W}}_{\mathrm{ic}}\right)-\hfill \\ \multicolumn{1}{c}{g_i(x_i)({\hat{w}}_{\mathrm{ih}}^T{\sigma }_{\mathrm{ih}}+k_{\mathrm{ih}}{e}_{\mathrm{ih}}+\underset{0}{\overset{t}{\int }}\left(\begin{array}{c}(k_{\mathrm{ih}}{\alpha }_{\mathrm{ih}}+{\gamma }_{\mathrm{ih}})e_{\mathrm{ih}}\hfill \\ \multicolumn{1}{c}{+{\delta }_{i1}\mathrm{sgn}(e_{\mathrm{ih}})}\end{array}\right)\mathrm{dt})}\end{array}$$

(57)

Theorem 3

Consider the n-DOF MRR system, with the subsystem dynamic model as the form of equation (2), which is completely unknown while designing the optimal controller. If the decentralized optimal SMC (equation (57)) is adopted for the MRR system, then the closed-loop robotic system is asymptotically stable.

Proof

Choosing the Lyapunov candidate function

$$V(t)=\sum _{i=1}^n{V}_i(t)=\sum _{i=1}^n\frac{1}{2}{s}_i^T{s}_i+J_i^{*}(s_i)$$

(58)

According to equations (32), (42), and (53), the time derivative of $V_i(t)$ is written as

$$\begin{array}{cc}\multicolumn{1}{c}{V(t)}& =\sum _{i=1}^n{\stackrel{·}{V}}_i(t)=\sum _{i=1}^n({s_i}^T{\stackrel{·}{s}}_i+{(\nabla J_i^{*}(s_i))}^T{\stackrel{·}{s}}_i)\hfill \\ \multicolumn{1}{c}{}& =\sum _{i=1}^n({s_i}^T\left(\begin{array}{c}{F}_i(x_i,x_{\mathrm{jd}})+g_i(x_i)u_i\hfill \\ \multicolumn{1}{c}{-{\stackrel{··}{x}}_{i1d}+k_i{\stackrel{·}{e}}_i+\mathrm{\Delta }{\psi }_i}\end{array}\right)+{(\nabla J_i^{*}(s_i))}^T{\stackrel{·}{s}}_i)\hfill \\ \multicolumn{1}{c}{}& =\sum _{i=1}^n\left(\begin{array}{c}{s_i}^T\left(\begin{array}{c}{g}_i(x_i)u_i+F_i(x_i,x_{\mathrm{jd}})\hfill \\ \multicolumn{1}{c}{-{\stackrel{··}{x}}_{i1d}+k_i{\stackrel{·}{e}}_i+\mathrm{\Delta }{\psi }_i}\end{array}\right)-s_i^T{Q}_i{s}_i\hfill \\ \multicolumn{1}{c}{-u_i^T{R}_i{u}_i-{\Upsilon }_i{(\nabla J_i^{*}(s_i))}^2-H_i^T{H}_i}\end{array}\right)\hfill \end{array}$$

(59)

As $F_i(x_i,x_{\mathrm{jd}})$ and $g_i(x_i)$ are Lipchitz, so has constant $L_{\mathrm{if}}>0,L_{\mathrm{ig}}>0$ and have $\Vert F_i\Vert \le L_{\mathrm{if}}\Vert s_i\Vert ,\Vert g_i\Vert \le w_{\mathrm{ig}}$

$$\begin{array}{cc}\multicolumn{1}{c}{\stackrel{·}{V}(t)}& =\sum _{i=1}^n({\stackrel{·}{V}}_i(t))\hfill \\ \multicolumn{1}{c}{}& \le \sum _{i=1}^n\left(\begin{array}{c}{L}_{\mathrm{if}}{\Vert s_i\Vert }^2+w_{\mathrm{ig}}\Vert s_i\Vert \Vert u_i\Vert +v_i\Vert s_i\Vert +H_i\Vert s_i\Vert \hfill \\ \multicolumn{1}{c}{-s_i^T{Q}_i{s}_i-u_i^T{R}_i{u}_i-{\Upsilon }_i{(\nabla J_i^{*}(s_i))}^2-H_i^T{H}_i}\end{array}\right)\hfill \\ \multicolumn{1}{c}{}& \le \sum _{i=1}^n\left(\begin{array}{c}{L}_{\mathrm{if}}{\Vert s_i\Vert }^2+w_{\mathrm{ig}}\Vert s_i\Vert \Vert u_i\Vert +v_i\Vert s_i\Vert +H_i\Vert s_i\Vert \hfill \\ \multicolumn{1}{c}{-{\lambda }_{min}(Q_i){\Vert s_i\Vert }^2-H_i^T{H}_i}\\ \multicolumn{1}{c}{-{\lambda }_{min}(R_i){\Vert u_i\Vert }^2-{\Upsilon }_i{(\nabla J_i^{*}(s_i))}^2}\end{array}\right)\hfill \\ \multicolumn{1}{c}{}& \le \sum _{i=1}^n\left(\begin{array}{c}{L}_{\mathrm{if}}{\Vert s_i\Vert }^2+\frac{1}{2}{w_{\mathrm{ig}}}^2{\Vert s_i\Vert }^2+\frac{1}{2}{\Vert u_i\Vert }^2+v_i\Vert s_i\Vert \hfill \\ \multicolumn{1}{c}{-{\lambda }_{min}(Q_i){\Vert s_i\Vert }^2-{\lambda }_{min}(R_i){\Vert u_i\Vert }^2}\\ \multicolumn{1}{c}{+H_i\Vert s_i\Vert -H_i^T{H}_i-{\Upsilon }_i{(\nabla J_i^{*}(s_i))}^2}\end{array}\right)\hfill \\ \multicolumn{1}{c}{}& =\sum _{i=1}^n\left(\begin{array}{c}[(\frac{1}{2}{w_{\mathrm{ig}}}^2-{\lambda }_{min}(Q_i)+L_{\mathrm{if}})\Vert s_i\Vert +H_i+v_i]\hfill \\ \multicolumn{1}{c}{·\Vert s_i\Vert -(-\frac{1}{2}+{\lambda }_{min}(R_i)){\Vert u_i\Vert }^2}\\ \multicolumn{1}{c}{-H_i^T{H}_i-{\Upsilon }_i{(\nabla J_i^{*}(s_i))}^2}\end{array}\right)\hfill \end{array}$$

(60)

where $v_i=-{\stackrel{··}{x}}_{i1d}+k_i{\stackrel{·}{e}}_i(t),k_i$ is a positive constant. There has been positive constants ${\kappa }_i,{\epsilon }_i,{\varsigma }_i,\Vert {\stackrel{··}{x}}_{i1d}\Vert \le {\epsilon }_i,\Vert H_i\Vert \le {\kappa }_i,\Vert {\stackrel{·}{e}}_i\Vert \le {\varsigma }_i$ . Then, one contains

$$\begin{array}{cc}\multicolumn{1}{c}{\stackrel{·}{V}(t)}& =\sum _{i=1}^n({\stackrel{·}{V}}_i(t))\hfill \\ \multicolumn{1}{c}{}& \le \sum _{i=1}^n\left(\begin{array}{c}-(-{\delta }_i+{\varphi }_{i1}\Vert s_i\Vert )\Vert s_i\Vert -{\varphi }_{i2}{\Vert u_i\Vert }^2\hfill \\ \multicolumn{1}{c}{-H_i^T{H}_i-{\Upsilon }_i{(\nabla J_i^{*}(s_i))}^2}\end{array}\right)\hfill \end{array}$$

(61)

where ${\varphi }_{i1}={\lambda }_{min}(Q_i)-L_{\mathrm{if}}-{w_{\mathrm{ig}}}^2/2$ and ${\varphi }_{i2}={\lambda }_{min}(R_i)-1/2,{\delta }_i={\epsilon }_i+k_i{\varsigma }_i+{\kappa }_i$ . Hence, one can get ${\stackrel{·}{V}}_i(t)\le 0$ while $s_i$ lies outside the compact set ${\Omega }_{i2}=\{s_i:\Vert s_i\Vert \le ({\epsilon }_i+k_i{\varsigma }_i+{\kappa }_i)/({\lambda }_{min}(Q_i)-L_{\mathrm{if}}-(1/2)w_{\mathrm{ig}}^2)\}$ , and if the conditions hold

$${\lambda }_{min}(Q_i)\ge L_{\mathrm{if}}+w_{\mathrm{ig}}^2/2,{\lambda }_{min}(R_i)\ge 1/2$$

(62)

Moreover, equation (61) means that $\stackrel{·}{V}(t)\le 0$ for $s_i\ne 0$ , while equation (62) is satisfied. According to the Lyapunov theory, optimal controller under the SMC is asymptotically stable. This theorem has been completed.

Remark 2

To solve the difficulty of addressing HJB equation, a local PI method is introduced by referring previous literature.^33–35 The iterative procedure of PI method with equation (9) is written in Appendix 1.

Simulations

Simulation setup

In this section, two MRRs are given for simulation. For configuration A, the dynamic formulation is given as

$$\begin{array}{c}M(q)=\left[\begin{array}{cc}0.6066+0.36\cos (q_2)\hfill & 0.1233+0.18\cos (q_2)\hfill \\ \multicolumn{1}{c}{0.1233+0.18\cos (q_2)}& 0.1233\hfill \end{array}\right]\hfill \\ C(q,\stackrel{·}{q})=\left[\begin{array}{cc}-0.36\sin (q_2){\stackrel{·}{q}}_2\hfill & -0.18{\stackrel{·}{q}}_2\sin (q_2)\hfill \\ \multicolumn{1}{c}{0.18\sin (q_2)({\stackrel{·}{q}}_1-{\stackrel{·}{q}}_2)}& 0.18{\stackrel{·}{q}}_1\sin (q_2)\hfill \end{array}\right]\hfill \\ G(q)=\left[\begin{array}{c}-5.88\sin (q_1+q_2)-17.64\sin (q_1)\hfill \\ \multicolumn{1}{c}{-5.88\sin (q_1+q_2)}\end{array}\right]\hfill \end{array}$$

and the ones for configuration B is given as

$$\begin{array}{c}M(q)=\left[\begin{array}{cc}0.17-0.1166{\cos }^2(q_2)\hfill & -0.06\cos (q_2)\hfill \\ \multicolumn{1}{c}{-0.06\cos (q_2)}& 0.1233\hfill \end{array}\right]\hfill \\ C(q,\stackrel{·}{q})=\left[\begin{array}{cc}-0.1166{\stackrel{·}{q}}_2\sin (2q_2)\hfill & -0.06\sin (q_2){\stackrel{·}{q}}_2\hfill \\ \multicolumn{1}{c}{\left(\begin{array}{c}0.06\sin (q_2){\stackrel{·}{q}}_2\hfill \\ \multicolumn{1}{c}{-0.0583\sin (2q_2){\stackrel{·}{q}}_1}\end{array}\right)}& -0.06{\stackrel{·}{q}}_1\sin (q_2)\hfill \end{array}\right]\hfill \\ G(q)=\left[\begin{array}{c}0\hfill \\ \multicolumn{1}{c}{-5.88\sin (q_1+q_2)}\end{array}\right]\hfill \end{array}$$

The desired trajectories of both configurations A and B are written as follows

$$\{\begin{array}{c}{x}_{1d}=0.4\sin t+0.1\cos (2t)\hfill \\ \multicolumn{1}{c}{x_{2d}=0.2\cos (3t)-0.3\sin (2t)}\end{array}$$

The NN weight vector estimation can be written as ${\hat{W}}_{\mathrm{ic}}=[{\hat{W}}_{\mathrm{ic}1},{\hat{W}}_{\mathrm{ic}2},{\hat{W}}_{\mathrm{ic}3}{]}^T$ . For configuration A, its initial value is given as ${\hat{W}}_{1c0}=[\begin{array}{ccc}20\hfill & 15\hfill & 10\hfill \end{array}{]}^T$ and ${\hat{W}}_{2c0}=[\begin{array}{ccc}20\hfill & 25\hfill & 30\hfill \end{array}{]}^T$ . For configuration B, its initial value is represented as ${\hat{W}}_{1c0}=[\begin{array}{ccc}20\hfill & 25\hfill & 30\hfill \end{array}{]}^T$ and ${\hat{W}}_{2c0}=[\begin{array}{ccc}20\hfill & 25\hfill & 30\hfill \end{array}{]}^T$ . The activation function is chosen by ${\delta }_{\mathrm{ic}}=[s_i^2,s_i^2,s_i^2]$ . Moreover, we select the control parameters $k_{\mathrm{ih}}=800,{\alpha }_{\mathrm{ih}}=350,R_i=Q_i=0.1I,{\gamma }_{\mathrm{ih}}=5,{\mathrm{\Gamma }}_{\mathrm{ih}}=0.1I,{\delta }_{i1}=0.5$ .

Simulation results

The simulation results are represented to improve the effectiveness of control torques, joint position tracking, the convergence of NN weights, and position tracking errors. The host computer CPU used in this article is Intel Core i7-7700 with 8.00 GB RAM, 3.60 GHz, and the suite of required software is MATLAB 2016a running in Windows 10. Two different control methods are used in the simulations that contain the existing NN-based optimal control method, such as the study by Dong and colleagues,^23,24 and the proposed ADP-based decentralized optimal SMC method.

Figure 1 shows the joint position tracking curves under the existing method. One can observe that the chattering effect is obvious at the first 2 s that is caused without dynamic decomposition and optimal compensation of the IDCs. Figure 2 illustrates the joint position curves of configuration A joint 1 under the proposed method. Comparing with Figure 1, we can obtain that the desired trajectories can be tracked within a very short time period, since the proposed method can verify the effectiveness. Figures 3 and 4 are joint position tracking error curves of configuration A under the existing method and proposed method. In Figure 3, we can obtain that the position errors are obvious at few seconds and the amplitude of steady-state error can be ±0.02 rad. From Figure 4, the position errors are tracking very fast and nearly at 0.

Figure 1.

Position tracking for configuration A with existing methods: (a) Joint 1 and (b) Joint 2.

Figure 2.

Position tracking for configuration A with proposed methods: (a) Joint 1 and (b) Joint 2.

Figure 3.

Position error for configuration A with existing methods: (a) Joint 1 and (b) Joint 2.

Figure 4.

Position error for configuration A with proposed methods: (a) Joint 1 and (b) Joint 2.

Figures 5 and 6 are the velocity tracking curves under the existing method and the proposed method of configuration A, respectively. Figures 7 and 8 illustrate the velocity error curves under the existing method and the proposed method of configuration A. Because the existing methods have not considered the compensation problem of the IDC effects, the tracking error is larger than the method that we proposed.

Figure 5.

Velocity tracking for configuration A with existing methods: (a) Joint 1 and (b) Joint 2.

Figure 6.

Velocity tracking for configuration A with proposed methods: (a) Joint 1 and (b) Joint 2.

Figure 7.

Velocity error for configuration A with existing methods: (a) Joint 1 and (b) Joint 2.

Figure 8.

Velocity error for configuration A with proposed methods: (a) Joint 1 and (b) Joint 2.

Figure 9 illustrates the control torque curve with the existing method of configuration A. For joint 1, it can be observed that the initial control torque is large and may bring the burden to the motors. That is because the NNs need time to learning for supply the big torque. For joint 2, it can be seen that at some time, the control torque has a sudden vibration. That is because the method without IDCs compensation. Figure 10 illustrates the control torque curve of configuration A with the proposed scheme. The output torques have been optimized with an appropriate behavior that may match up the output power of motors.

Figure 9.

Control torque for configuration A with existing methods: (a) Joint 1 and (b) Joint 2.

Figure 10.

Control torque for configuration A with proposed methods: (a) Joint 1 and (b) Joint 2.

Figure 11 shows the critic NN curves of configuration A joints 1 and 2 under the proposed method. Because of the implementation process of PI and critic NN training, the convergence results of weights can be obtained before 1 s. The weights of the critic NN converge to $W_{c1}=[\begin{array}{ccc}7.2897\hfill & 2.2897\hfill & -2.7103\hfill \end{array}{]}^T$ and $W_{c2}=[\begin{array}{ccc}-4.9820\hfill & 0.0180\hfill & 5.0180\hfill \end{array}{]}^T$ . Unlike the conventional ADP-based optimal control methods that rely on both action and critic NNs, the optimal control is obtained using the critic NN only and the training of action NN is no longer needed. It infers that the computational burden can be reduced effectively.

Figure 11.

Critic NN for configuration A with proposed methods: (a) Joint 1 and (b) Joint 2.

Figures 12–22 represent the trajectory tracking curves, position error curves, velocity curves, velocity error curves and control torque, and convergence results of weights for configuration B. We can obtain the similar results. The conclusion can be received comparing with configuration A. It improves the proposed method without the requirements of adjusting parameters. The weights of critic NN converge to $W_{c1}=[\begin{array}{ccc}19.0669\hfill & 23.9518\hfill & 28.8368\hfill \end{array}{]}^T$ and $W_{c2}=[\begin{array}{ccc}-4.9960\hfill & 0.0040\hfill & 5.0040\hfill \end{array}{]}^T$ .

Figure 12.

Position tracking for configuration B with existing methods: (a) Joint 1 and (b) Joint 2.

Figure 13.

Position tracking for configuration B with proposed methods: (a) Joint 1 and (b) Joint 2.

Figure 14.

Position error for configuration B with existing methods: (a) Joint 1 and (b) Joint 2.

Figure 15.

Position error for configuration B with proposed methods: (a) Joint 1 and (b) Joint 2.

Figure 16.

Velocity tracking for configuration B with existing methods: (a) Joint 1 and (b) Joint 2.

Figure 17.

Velocity tracking for configuration B with proposed methods: (a) Joint 1 and (b) Joint 2.

Figure 18.

Velocity error for configuration B with existing methods: (a) Joint 1 and (b) Joint 2.

Figure 19.

Velocity error for configuration B with proposed methods: (a) Joint 1 and (b) Joint 2.

Figure 20.

Control torque for configuration B with existing methods: (a) Joint 1 and (b) Joint 2.

Figure 21.

Control torque for configuration B with proposed methods: (a) Joint 1 and (b) Joint 2.

Figure 22.

Critic NN for configuration B with proposed methods: (a) Joint 1 and (b) Joint 2.

For the simulations, one can obtain the decentralized optimal SMC method and can guarantee stability and accuracy.

Conclusion

This article proposes a model-free decentralized SMC method for MRR system via PI scheme–based ADP. First, the dynamic formulation is expressed by a synthesis of joint subsystems with IDC effects. Then, based on the SMC technique, the decentralized optimal control problem is transformed into the optimal compensation issue of unknown dynamics of each subsystem. Based on ADP and the PI theory, the HJB equation can address by using a critic NN and one can obtain the optimal control policy. According to the Lyapunov theory, the closed-loop MRR systems are guaranteed UUB. Finally, simulations are verified by the effectiveness of the method.

As is known, ADP-based control methods have been successfully used to address the optimal control problems in battery management, residential energy management, water–gas shift reaction, and coal gasification process in theory. However, the effectiveness analysis of the mentioned works all relies on numerical analyses and simulation results. Indeed, a drawback of the researches of ADP-based optimal control methods, in a common view, is the development of experimental researches and practical applications for physical systems, especially for robotic systems. How to address the problems in establishing an experimental platform of robotic systems that satisfy the real-time and accurate requirements is a key problem for implementing the proposed decentralized SMC method to actual MRR systems, and it is also our future research topic.