Nonlinear Model Predictive Control for Mobile Robot Using Varying-Parameter Convergent Differential Neural Network

Abstract

The mobile robot kinematic model is a nonlinear affine system, which is constrained by velocity and acceleration limits. Therefore, the traditional control methods may not solve the tracking problem because of the physical constraint. In this paper, we present the nonlinear model predictive control (NMPC) algorithm to track the desired trajectory based on neural-dynamic optimization. In the proposed algorithm, the NMPC scheme utilizes a new neural network named the varying-parameter convergent differential neural network (VPCDNN) which is a Hopfifield-neural network structure with respect to the differential equation theory to solve the quadratic programming (QP) problem. The new network structure converges to the global optimal solution and it is more efficient than traditional numerical methods. In the simulation, we verify that the proposed method is able to successfully track reference trajectories with a two-wheel mobile robot. The experimental validation has been conducted in simulation and the results show that the proposed method is able to precisely track the trajectory maintaining a high robustness based on the VPCDNN solver.

1. Introduction

Over past decades, the motion control of mobile robots has increasingly attracted a lot of researchers, especially for applications in industry and transportation. The two-wheel mobile robot is a special differential wheeled robot and the physical system is the nonholonomic system, so the two-wheel robot system is in asymptotic stabilization when it is controlled by time-varying methods based on Brockett’s theory [1].

There are various approaches that have been developed for mobile robot control and in some applications under certain conditions. In Reference [2], the authors presented the back-stepping control algorithm for the curved weld seam tracking of a mobile robot and achieved good tracking precision. In Reference [3], a nonlinear control algorithm based on sliding mode control (SMC) was developed for mobile robots, where the four wheel motion mode was converted to an inverted pendulum mode, which has a small turning radius and high-speed observation features. In Reference [4], the iterative learning control method was developed for path-tracking, which can reduce the errors in Cartesian space for repeat motions. In Reference [5], the fuzzy control is applied to autonomous mobile robot control. In Reference [6,7], the author presents a robust adaptive feedback controller for a tractor trailer wheeled robot under uncertain parameters.

In Reference [8], a non-iterative controller is presented for a nonlinear mobile robot system without constraint. The new method has a low worst-cost ration than nonlinear model predictive control method (NMPC), but cannot ensure safety in its running because no physical limits are considered. In Reference [9], an adaptive optimal control method is presented to solve the $H\infty$ problem under changing environmenalt conditions, which can adapt to a complex environment, but it seems that the gradient-descent based solver is used to solve the optimal problem.

For safety reasons, the velocity and acceleration constraints (equivalent to velocity incremental constraint) and other constraints need to be taken into account during operation, which is not suitable for some special applications with physical limits and is very dangerous in real time, running over the limit speed. In this paper, the nonlinear model predictive control algorithm is employed to optimize the performance of the mobile robot during the tracking of the reference trajectories. Due to the excellent performance of the NMPC control method with nonlinear input-output over a finite horizon, it is widely used to track trajectories with the kinematic constraints of a mobile robot.

There are many traditional model predictive control methods developed for the tracking of mobile robots. In References [10,11], the linear model predictive control is used for tracking problems with fast small-scale solvers. However, for strong nonlinear robot systems, linear models find it difficult to describe the complex motion process, and it is hard to predict the next step accurately. Therefore, in this paper, the nonlinear model predictive control method is proposed to control the mobile robot system. Compared with the linear MPC method, NMPC describes the complex dynamic processes of the robot system more accurately.

The NMPC method has the nonlinear input-output over a finite horizon and it is widely used for tracking in mobile robots or auto driving vehicle or mobile robots [12,13]. In Reference [14], NMPC is applied to control vehicle velocity which demonstrated the NMPC’s energy cost.

In References [15,16], a novel NMPC method is presented to control vehicles, which has complex nonlinear terms and uncertainties. In References [17,18], the stochastic MPC method is employed to penalize undesired factors. In References [19,20], the traditional numerical optimization method is used for solving the quadratic programming (QP) problem of MPC. However, it has the computational burden problem which cannot meet online computing in a real time running application.

In Reference [21], the authors present the NMPC method to follow customers and maintain a certain distance, but the solver for NMPC seems to be based on the gradient-descent method. In this case, the methods cannot track the theoretical solutions because the theoretical solutions change over time, which causes an undesired performance and cannot converge to 0 in a limited amount of time. Therefore, time-varying neural network methods have been designed to solve the optimization problem of NMPC.

There are many works on the neural network methods for solving MPC optimization [22]. In Reference [23], the nonlinear model predictive control scheme using general projection neural networks (GPNN) is employed to optimize the chained systems with nonholonomic constraints. The GPNN is designed for the QP problem and it globally converges to the optimal solution with a quick convergence rate and a low computational burden. In addition, in References [24,25,26], the primal dual neural network (PDNN) is used to optimize the QP problem of NMPC, where the PDNN has no matrix inversion and matrix-matrix multiplication, which is suitable for the online optimization in References [27,28].

However, these are all time-invariant methods based on gradient-descent for optimization. Inspired by the works mentioned above, in this paper a new neural network named the varying-parameter convergent differential neural network (VPCDNN) method is employed to solve the QP problems of MPC. Our VPCDNN has the same neurons in the network as the decision variables for NMPC optimization. The proposed VPCDNN can obtain an optimal solution for the QP problem in real-time application.

The main contributions of this paper are:

• A varying parameter convergent differential neural network method is proposed to solve the time varying QP problem of MPC and all state variables can converge quickly to the optimal value using the neural network with physical. This is the first time that be presented to optimize the MPC problem.
• The convergence analysis of VPCDNN and the simulation results demonstrate good performance on the convergence speed and robustness of VPCDNN.

The remainder of this paper is organized as follows. The kinematic model of mobile robots and the nonlinear model predictive control scheme are described in Section 2. The proposed VPCDNN method for solving the online QP problem of MPC is detailed in Section 3. The convergence and robustness analysis of VPCDNN is shown in Section 4. Simulation results are reported and discussed in Section 5. Finally, Section 6 concludes this paper.

2. Model Predictive Control Scheme

In this section, the kinematic model of a mobile robot is reformulated as a nonlinear affine system, then the model predictive control scheme is presented for the tracking problem. The varying-parameter neural network base model predictive control scheme is proposed to solve the QP problem with physical constraints.

2.1. Mobile Robot Control System

The kinematic model of two wheels mobile robot are shown in Figure 1. According to the relationship between robot velocity v and two driving wheels velocity ($v_L$, $v_R$), the robot velocity and angle velocity are expressed as: $v=(v_L+v_R)/2$, $\omega =(v_L-v_R)/l_f$, respectively. $l_f$ denotes the distance of two wheels. The kinematic model formulation of mobile robot is expressed as,

$$\dot{X}=\left[\begin{array}{c}\dot{x}\\ \dot{y}\\ \dot{\theta }\end{array}\right]=\left[\begin{array}{c}vcos\theta \\ vsin\theta \\ \omega \end{array}\right]=\left[\begin{array}{cc}cos\theta & 0\\ sin\theta & 0\\ 0& 1\end{array}\right]u$$

(1)

where $(x,y)$ is the mobile robot position in Cartesian space and $\theta$ is the orientation; $u={(v,\omega )}^T$ is the control input and X is the state vector. From the formulation in (1), we can obtain a kinematic model of the reference trajectory by state vector $X_r={(x_r,y_r,{\theta }_r)}^T$ and control input $u_r={(v_r,{\omega }_r)}^T$ and it is expressed as:

$${\dot{X}}_r=\left[\begin{array}{c}{\dot{x}}_r\\ {\dot{y}}_r\\ {\dot{\theta }}_r\end{array}\right]=\left[\begin{array}{c}{v}_{r}cos{\theta }_r\\ v_{r}sin{\theta }_r\\ {\omega }_r\end{array}\right]=\left[\begin{array}{cc}cos{\theta }_r& 0\\ sin{\theta }_r& 0\\ 0& 1\end{array}\right]u_r$$

(2)

Therefore, we can get the robotic kinematic errors,

$$X_e=\left[\begin{array}{ccc}cos\theta & sin\theta & 0\\ -sin\theta & cos\theta & 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}{x}_r-x\\ y_r-y\\ {\theta }_r-\theta \end{array}\right]$$

(3)

where $X_e={[x_e,y_e,{\theta }_e]}^T$. The derivative of ${\dot{X}}_e$ can be obtained from the error state in (3),

$$\left\{\begin{array}{c}{\dot{x}}_e=\omega y_e-v+v_{r}cos{\theta }_e\hfill \\ {\dot{y}}_e=-\omega x_e+v_{r}sin{\theta }_e\hfill \\ {\dot{\theta }}_e={\omega }_r-\omega \hfill \end{array}\right.$$

(4)

The control inputs are reformulated as:

$$u_e=\left[\begin{array}{c}{u}_1\\ u_2\end{array}\right]=\left[\begin{array}{c}{v}_{r}cos{\theta }_e-v\\ {\omega }_r-\omega \end{array}\right]$$

(5)

Substituting (5) into (4), the kinematic model of $X_e$ can be rewritten as

$${\dot{X}}_e=\left[\begin{array}{c}{\dot{x}}_e\\ {\dot{y}}_e\\ {\dot{\theta }}_e\end{array}\right]=\left[\begin{array}{ccc}0& \omega & 0\\ -\omega & 0& 0\\ 0& 0& 0\end{array}\right]\left[\begin{array}{c}{x}_e\\ y_e\\ {\theta }_e\end{array}\right]+\left[\begin{array}{c}{u}_1\\ v_{r}sin{\theta }_e\\ u_2\end{array}\right]$$

(6)

The formulation in (6) can be reformulated as the following equation:

$${\dot{\overline{X}}}_e=\left[\begin{array}{ccc}0& {\omega }_r& 0\\ -{\omega }_r& 0& v_r\\ 0& 0& 0\end{array}\right]{\overline{X}}_e+\left[\begin{array}{cc}1& 0\\ 0& 0\\ 0& 1\end{array}\right]u_e$$

(7)

Therefore, the motion control problem of the mobile robot in (1) is converted into a stabilization problem of a nonlinear affine system in (7).

2.2. Nonlinear Model Predictive Control

The NMPC can be converted into a closed-loop optimal control problem in a finite time-horizon, which is subject to the input vector and state vector constraints. The NMPC formulation can be presented as

$$X\left(j+1\right)={\psi }_1\left(X\left(j\right)\right)+{\psi }_2\left(X\left(j\right)\right)u\left(j\right)$$

(8)

the input and state vector constraints are defined as:

$$\begin{array}{c}X\left(j\right)\in R^X,j=1,2,\dots ,N\hfill \\ u\left(j\right)\in R^U,j=1,2,...N_u\hfill \end{array}$$

where X and u denote state vector and input vector, respectively. In a nonlinear affine system (8), ${\psi }_1(.)$ and ${\psi }_2(.)$ denote the nonlinear continuous functions. The initial conditions satisfy ${\psi }_1\left(0\right)=0$, and $N\ge 1$ and $1\le N\le N_u$ that denotes the prediction horizon.

In each sampling period, the optimal input vector can be obtained by a given cost function online optimization, thus the NMPC scheme is formulated by iterative solution of optimal control problems. The cost function can be described as the following form function:

$$S(X,u)=\sum _{i=j}^{j+N-1}{L}_1\left(X\left(i\right),u\left(i\right)\right)+L_2\left(X(j+N)\right)$$

(9)

where $L_1(X,u)$ denotes the immediate cost and the satisfying condition is $L_1(0,0)=0$. The cost function $S(X,u)$ can be defined as a quadratic form as follows:

$$S\left(j\right)=\sum _{i=1}^N{∥X\left(j+i|j\right)∥}_Q^2+\sum _{i=0}^{N_u-1}{∥∆u\left(j+i|j\right)∥}_R^2$$

(10)

where $X(j+i|j)$ is the predicted future horizon state; $∆u(j+i|j)=u(j+i\left|j\right)-u(j-1+i|j)$ which is the increment of the input vector; the parameters R and Q represent the constant design weight matrix; the symbol $∥·∥$ represents the Euclidean norm of the corresponding vector. The Equation (8) can be reformulated as:

$$X_e\left(j+1\right)={\psi }_1\left(X_e\left(j\right)\right)+{\psi }_2\left(X_e\left(j\right)\right)u_e\left(j\right)$$

(11)

$$\begin{array}{ccc}\hfill \mathit{subject} \mathit{to}& & u^{-}\le u_e\left(j\right)\le u^{+}\hfill \end{array}$$

(12)

$$\begin{array}{ccc}& & ∆u^{-}\le ∆u\left(j\right)\le ∆u^{+}\hfill \end{array}$$

(13)

$$\begin{array}{ccc}& & X^{-}\le X_e\left(j\right)\le X^{+}\hfill \end{array}$$

(14)

$${\psi }_1\left(X_e\right)=\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\end{array}\right]+T\left[\begin{array}{c}{\omega }_r{x}_2\\ -{\omega }_r{x}_1+v_r{x}_3\\ 0\end{array}\right]$$

$${\psi }_2\left(X_e\right)=T\left[\begin{array}{cc}1& 0\\ 0& 0\\ 0& 1\end{array}\right]$$

where $X_e={\left[x_1,x_2,x_3\right]}^T={\left[x_e,y_e,{\theta }_e\right]}^T$ is the state vector, and T is the sampling period; $u_e={\left[u_1,u_2\right]}^T$, $u_1=v_{r}cos\left({\theta }_e\right)-v$, $u_2={\omega }_r-\omega$ is the input vector; $(u^{-},u^{+})$ are the upper and lower limits of input variable, and $(∆u^{-},∆u^{+})$, $(x^{-},x^{+})$ are also the upper and lower limits of the input increment variable and state variable, respectively.

To construct the QP problem for online optimization, we introduce the vectors as:

$$\begin{array}{c}\overline{X}={\left[X_e(j+1|j),...,X_e(j+N|j)\right]}^T\end{array}$$

(15)

$$\begin{array}{c}\overline{u}\left(j\right)={\left[u_e\left(j\right|j),...,u_e(j+N_u-1|j)\right]}^T\end{array}$$

(16)

$$\begin{array}{c}∆\overline{u}\left(j\right)={\left[∆u_e\left(j\right|j),...,∆u_e(j+N_u-1|j)\right]}^T\end{array}$$

(17)

According to (11), and (15)–(17), we can get the predicted output,

$$\overline{X}\left(j\right)=H∆\overline{u}\left(j\right)+\widehat{{\nu }_1}+\widehat{{\nu }_2}$$

(18)

$$H=\left[\begin{array}{ccc}{\nu }_2\left(X_e\left(j\right|j-1)\right)& \cdots & 0\\ {\nu }_2\left(X_e(j+1|j-1)\right)& \cdots & 0\\ ⋮& \ddots & ⋮\\ {\nu }_2\left(X_e(j+N-1|j-1)\right)& \cdots & {\nu }_2\left(X_e(j+N-1|j-1)\right)\end{array}\right]$$

$$\widehat{{\nu }_1}=\left[\begin{array}{c}{\nu }_1\left(X_e\left(j\right|j-1)\right)\\ {\nu }_1\left(X_e(j+1|j-1)\right)\\ ⋮\\ {\nu }_1\left(X_e(j+N-1|j-1)\right)\end{array}\right]$$

$$\widehat{{\nu }_2}=\left[\begin{array}{c}{\nu }_2\left(X_e\left(j\right|j-1)u(k-1)\right)\\ {\nu }_2\left(X_e(j+1|j-1)u(j-1)\right)\\ ⋮\\ {\nu }_2\left(X_e(j+N-1|j-1)u(j-1)\right)\end{array}\right]$$

where $H\in R^{3N\times 2N_u}$, $\widehat{{\nu }_1}\in R^{3N}$, and $\widehat{{\nu }_2}\in R^{3N}$. Then, the optimization problem in (10) can be rewritten,

$$\begin{array}{cc}\hfill \mathit{minimum}& {∥H∆\overline{u}\left(j\right)+\widehat{{\nu }_1}+\widehat{{\nu }_2}∥}_Q^2+{∥∆u\left(j\right)∥}_R^2\end{array}$$

(19)

$$\begin{array}{cc}\hfill \mathit{subject} \mathit{to}\hspace{1em}& ∆{\overline{u}}^{-}\le ∆\overline{u}(j+1)\le ∆{\overline{u}}^{+}\end{array}$$

(20)

$$\begin{array}{c}{\overline{u}}^{-}\le \overline{u}(j-1)\le {\overline{u}}^{+}\end{array}$$

(21)

$$\begin{array}{c}{\overline{u}}^{-}\le \overline{u}(j-1)+\widehat{I}∆\overline{u}\left(j\right)\le {\overline{u}}^{+}\end{array}$$

(22)

$$\begin{array}{c}{X}^{-}\le \widehat{{\nu }_1}+\widehat{{\nu }_2}+H∆\overline{u}\left(j\right)\le X^{+}\end{array}$$

(23)

where $\widehat{I}=\left[\begin{array}{cccc}I& 0& 0& 0\\ I& I& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ I& I& \cdots & I\end{array}\right]\in R^{2N_u\times 2N_u}$

The optimization problem in (19)–(23) can be reformulated as a QP problem:

$$\begin{array}{cc}\hfill \mathit{minimum}& \frac{1}{2}∆{\overline{u}}^{T}M∆\overline{u}+c^T∆\overline{u}\end{array}$$

(24)

$$\begin{array}{cc}\hfill \mathit{subject} \mathit{to}\hspace{1em}& E∆\overline{u}\le r\end{array}$$

(25)

$$\begin{array}{c}∆{\overline{u}}^{-}\le ∆\overline{u}\le ∆{\overline{u}}^{+}\end{array}$$

(26)

and the parameter details are given,

$$M=2(H^T{\eta }_{1}H+{\eta }_2)\in R^{2N_u\times 2N_u}$$

$$c=-2H^T{\eta }_1(\widehat{{\nu }_2}+\widehat{{\nu }_1})\in R^{2N_u}$$

$$\mathrm{E}=\left[\begin{array}{c}-\widehat{I}\\ \widehat{I}\\ -H\\ H\end{array}\right]\in R^{(4N_u+6N)\times 2N_u},r=\left[\begin{array}{c}-u^{-}+u(j-1)\\ u^{+}-u(j-1)\\ -X^{-}+\widehat{{\nu }_1}+\widehat{{\nu }_2}\\ X^{+}-\widehat{{\nu }_1}-\widehat{{\nu }_2}\end{array}\right]\in R^{4N_u+6N}$$

3. Varying-Parameter Convergent Differential Neural Network (VPCDNN)

Firstly, we can rewrite the QP problem as follows:

$$\begin{array}{cc}\hfill \mathit{minimum}& \frac{1}{2}{d}^{T}Md+c^{T}d\end{array}$$

(27)

$$\begin{array}{cc}\hfill \mathit{subject} \mathit{to}\hspace{1em}& Ed\le r\end{array}$$

(28)

$$\begin{array}{c}{\zeta }^{-}\le d\le {\zeta }^{+}\end{array}$$

(29)

where $d=∆\overline{u}$, ${\zeta }^{-}=∆{\overline{u}}^{-},{\zeta }^{+}=∆{\overline{u}}^{+}$, and $d\in [{\zeta }^{-},{\zeta }^{+}]$.

In order to find the optimal solution of the time-varying convex quadratic programming problem, we design the Lagrange function for Equations (27)–(29) as follows:

$$L\left(d,y\right)=\frac{d^{T}Md}{2}+c^{T}d+y^T\left(Ed-r\right)$$

(30)

where variable y is the Lagrange-multiplier. From the Lagrangian theorem in Reference [29], if the derivatives of $\frac{\partial L}{\partial d}$, $\frac{\partial L}{\partial y}$ are continuous, we give the following algebraic equations of the Lagrange necessary condition by setting the partial derivatives of $L(d,y)$ to zero,

$$\left\{\begin{array}{c}\frac{\partial L}{\partial d}=Md+c+E^{T}y=0\hfill \\ \frac{\partial L}{\partial y}=Ed-r=0\hfill \end{array}\right.$$

(31)

We can rewrite the equation in (31) as the matrix form,

$$Wp=l$$

(32)

$$\begin{array}{c}W=\left[\begin{array}{cc}M& E^T\\ E& 0\end{array}\right]\in R^{\left(6N_u+6N\right)\times \left(6N_u+6N\right)}\hfill \\ p=\left[\begin{array}{c}d\\ y\end{array}\right]\in R^{6N_u+6N},l=\left[\begin{array}{c}-c\\ r\end{array}\right]\in R^{6N_u+6N}\hfill \end{array}$$

where W is a square matrix and invertible; p denotes the vector that needs to be solved by physical limits. It should be noted that the vectors c and r are time-varying parameters. Therefore, the time-varying solver is considered to solve the time-varying QP problem in (27)–(29).

Actually, the optimal solution is under the condition that the equation in (32) is equivalent to zero. So, the error expression is defined as

$$e\left(t\right)=Wp\left(t\right)-l\left(t\right)$$

(33)

where the error $e\left(t\right)\in R^{6Nu+6N}$ and we hope the variable e is close to zero. Therefore, the varying-parameter convergent differential neural network is proposed for the problem in (33),

$$\dot{e}\left(t\right)=-\lambda exp\left(\kappa t\right)P_{\mathsf{\Omega }}\left(e\left(t\right)\right)$$

(34)

$$P_{\mathsf{\Omega }}\left(e\left(t\right)\right)=\left\{\begin{array}{c}{e}_i^{-},\mathrm{if}{e}_i<e_i^{-}\hfill \\ \hfill \\ e_i,\mathrm{if}{e}_i^{-}\le e_i\le e_i^{+},\forall i\in \left\{1,\cdots ,6N_u+6N\right\}\hfill \\ \hfill \\ e_i^{+},\mathrm{if}{e}_i>e_i^{+}\hfill \end{array}\right.$$

where the parameters $\lambda$ and $\kappa$ are positive constants, which are used to scale the convergence rate; the symbol $P_{\mathsf{\Omega }}(·)$ is the active function. Then, substituting (33) into (34), the differential equation of the solver for the QP problem can be concluded as

$$W\dot{p}\left(t\right)=-\dot{W}p\left(t\right)-\lambda exp\left(\kappa t\right)P_{\mathsf{\Omega }}(Wp\left(t\right)-l\left(t\right))+\dot{l}\left(t\right)$$

(35)

The block diagram of the VPCDNN model (35) is shown in Figure 2. Finally, for the online solving process, the neural network consists of $N_p$ neurons and the neural network is constructed as

$$\begin{array}{c}{\dot{p}}_i=-{\displaystyle \sum _{j=1}^{N_p}}{\dot{w}}_{ij}{p}_i+{\displaystyle \sum _{j=1}^{N_p}}\left({\rho }_{ij}-w_{ij}\right){\dot{p}}_i-\lambda exp\left(\kappa t\right)P_{\mathsf{\Omega }}\left\{{\displaystyle \sum _{j=1}^{N_p}}\left(w_{ij}{p}_i-l_i\right)\right\}+{\dot{l}}_i\hfill \end{array}$$

(36)

where $p_i$, ${\dot{p}}_i$, $l_i$, ${\dot{l}}_i$ are the i-th elements of variable p, $\dot{p}$, $l_i$, ${\dot{l}}_i$, respectively; $w_{ij}$, ${\dot{w}}_{ij}$, ${\eta }_{ij}$, ${\dot{\eta }}_{ij}$ are the i-th row and j-th column elements of variable W, I, respectively. Figure 3 exhibits the detailed structure of the VPCDNN model.

4. Convergence and Robustness Analysis of VPCDNN

4.1. Convergence Analysis

It should be noted that the VPCDNN algorithm is a special type of Hopfield-neural network because the network fits the Hopfield-neural network structure with reference to the differential equation theory.

If there exists the optimal solution $d^{\ast }$ for (27)–(29), the vector p converges to equilibrium point $p^{\ast }$ from all initial states $p\left(0\right)$, which means the VPCDNN can achieve the convergence and the stable point $p^{\ast }$ of VPCDNN is the optimal solution of the QP problem. Firstly, we design the Lyapunov function of error variable $e\left(t\right)$ and then find the minimum point of the Lyapunov function. The detailed proof of convergence is given as follows.

Theorem 1.

For the time-varying QP problem in (27)–(29), assume that there is the optimal solution $d^{\ast }$ for the QP problem when the monotone nondecreasing activation function of $P_{\mathsf{\Omega }}$ acts on the error variable. From any initial state $p\left(0\right)\in R^{6N_u+6N}$, the state vector $p\left(t\right)={[d^T\left(t\right),y^T\left(t\right)]}^T$ of VPCDNN globally converges to equilibrium point $p^{\ast }\left(t\right)={\left[d^{\ast T}\left(t\right),y^{\ast T}\left(t\right)\right]}^T$, where the first n-elements of $p^{\ast }$ ($x^{\ast }$) are the optimal solution of (27)–(29).

Proof. 

Considering the candidate Lyapunov function as

$$F\left(t\right)=\frac{1}{2}{e}^T\left(t\right)e\left(t\right),\forall e\left(t\right)\ne 0$$

(37)

It is obvious that the Lyapunov function in (37) is a positive function. The time derivative of $F\left(t\right)$ can be obtained as

$$\dot{F}\left(t\right)=\frac{dF\left(t\right)}{dt}=e^T\left(t\right)\dot{e}\left(t\right)$$

(38)

Then, substituting (34) into (38), we obtain

$$\begin{array}{c}\dot{F}\left(t\right)=-\lambda exp\left(\kappa t\right)e^T\left(t\right)P_{\mathsf{\Omega }}\left\{e\left(t\right)\right\}\hfill \\ =-\lambda exp\left(\kappa t\right){\displaystyle \sum _{i=1}^{N_p}}{e}_i\left(t\right)P_{\mathsf{\Omega }}\left(e_i\left(t\right)\right)\hfill \end{array}$$

(39)

where $e_i\left(t\right)$ represents the ith element of variable $e\left(t\right)$; $P_{\mathsf{\Omega }}\left\{e\left(t\right)\right\}$ is the ith projection element of variable $P_{\mathsf{\Omega }}\left\{e\left(t\right)\right\}$. Since $P_{\mathsf{\Omega }}\left\{e\left(t\right)\right\}$ is the monotone nondecreasing projection function, it is easy to obtain

$$e_i\left(t\right)P_{\mathsf{\Omega }}\left(e_i\left(t\right)\right)=\left\{\begin{array}{c}>0,\mathrm{if}{e}_i\left(t\right)\ne 0\hfill \\ =0,\mathrm{if}{e}_i\left(t\right)=0\hfill \end{array}\right.$$

(40)

Therefore, we can obtain the following,

$$\dot{F}\left(t\right)=\left\{\begin{array}{c}>0,\mathrm{if}{e}_i\left(t\right)\ne 0\hfill \\ =0,\mathrm{if}{e}_i\left(t\right)=0\hfill \end{array}\right.$$

(41)

According to the (41), we know that if and only if $e_i\left(t\right)=0$, $\dot{F}\left(t\right)$ is equal to zero. Therefore, $p\left(t\right)-p^{\ast }\left(t\right)$ globally converges to 0 based on Lyapunov theory [30], which also means $d\left(t\right)-d^{\ast }\left(t\right)$ globally converges to zero. So, the proof of VPCDNN is finished. □

4.2. Robustness Analysis

Considering the disturbance model in (35) as follows,

$$\begin{array}{c}W\dot{p}\left(t\right)=-(\dot{W}+∆B\left(t\right))p\left(t\right)-\lambda exp\left(\kappa t\right)P(Wp\left(t\right)-l\left(t\right))+\dot{l}\left(t\right)+∆\xi \left(t\right)\hfill \end{array}$$

(42)

where $∆B\in R^{\left(6N_u+6N\right)\times \left(6N_u+6N\right)}$ is the disturbing term of W; $∆\xi \in R^{n+m}$ is the error from the VPCDNN model.

Theorem 2.

If $∥∆B\left(t\right)∥\le {\mu }_B$, $∥∆\xi \left(t\right)∥\le {\mu }_{\xi }$, $∥W^{-1}∥\le {\mu }_W$, $∥l\left(t\right)∥\le {\mu }_l$, and ${\mu }_B$, where ${\mu }_{\xi }$, ${\mu }_W$, ${\mu }_l$ are all positive parameters, the error $e\left(t\right)$ will converge to 0 under the condition of $\lambda aexp\left(\kappa t\right)-{\mu }_B{\mu }_A>0$.

Proof. 

Substituting the (33), (34) into (42), we can conclude

$$\begin{array}{c}\dot{e}\left(t\right)=-\lambda exp\left(\kappa t\right)P_{\mathsf{\Omega }}\left(e\left(t\right)\right)-∆B\left(t\right)W^{-1}e\left(t\right)+∆\xi \left(t\right)-∆B\left(t\right)W^{-1}l\left(t\right)\hfill \end{array}$$

(43)

We choose a Lyapunov function as

$$V\left(t\right)=\frac{1}{2}{e}^T\left(t\right)e\left(t\right)=\frac{1}{2}\sum _{j=1}^{6N_u+6N}{e}_i^2\left(t\right)$$

(44)

where $V\left(t\right)$ is the non-negative variable.

Then, we can obtain the time derivative of $V\left(t\right)$,

$$\begin{array}{c}V\left(t\right)=e^T\left(t\right)\dot{e}\left(t\right)\hfill \\ =e^T\left(t\right)(-\lambda exp\left(\kappa t\right)P_{\mathsf{\Omega }}\left(e\left(t\right)\right)-∆B\left(t\right)W^{-1}e\left(t\right)+∆\xi \left(t\right)-∆B\left(t\right)W^{-1}l\left(t\right))\hfill \\ =-\lambda exp\left(\kappa t\right)e^T\left(t\right)P_{\mathsf{\Omega }}\left(e\left(t\right)\right)+e^T\left(t\right)\psi \left(t\right)e\left(t\right)+e^T\left(t\right)∆\xi \left(t\right)+e^T\left(t\right)\left(-∆B\left(t\right)W^{-1}l\left(t\right)\right)\hfill \\ =-\lambda exp\left(\kappa t\right)e^T\left(t\right)P_{\mathsf{\Omega }}\left(e\left(t\right)\right)+e^T\left(t\right)\frac{\psi \left(t\right)+{\psi }^T\left(t\right)}{2}e\left(t\right)+e^T\left(t\right)∆\xi \left(t\right)+e^T\left(t\right)\left(-∆B\left(t\right)W^{-1}l\left(t\right)\right)\hfill \end{array}$$

(45)

where $\psi \left(t\right)=-∆B\left(t\right)W^{-1}$. Moreover,

$$\begin{array}{c}{e}^T\left(t\right)\frac{\psi \left(t\right)+{\psi }^T\left(t\right)}{2}e\left(t\right)\le e^T\left(t\right)e\left(t\right)\left|{\lambda }_{max}\left(\frac{\psi \left(t\right)+{\psi }^T\left(t\right)}{2}\right)\right|\hfill \\ \le e^T\left(t\right)e\left(t\right)∥∆B\left(t\right)W^{-1}\left(t\right)∥\le e^T\left(t\right)e\left(t\right){\mu }_B{\mu }_W\hfill \end{array}$$

(46)

$$e^T\left(t\right)∆\xi \left(t\right)\le \sum _{i=1}^{6N_u+6N}\left|e_i\right|{\mu }_{\xi }$$

(47)

$$\begin{array}{c}{e}^T\left(t\right)\left(-∆B\left(t\right)W^{-1}l\left(t\right)\right)\le {\displaystyle \sum _{i=1}^{6N_u+6N}}\left|e_i\right|·∥-∆B\left(t\right)W^{-1}l\left(t\right)∥\le {\displaystyle \sum _{i=1}^{6N_u+6N}}\left|e_i\right|{\mu }_B{\mu }_{\xi }{\mu }_l\hfill \end{array}$$

(48)

Therefore, substituting (46)–(47) into (45),

$$\begin{array}{c}\dot{V}\left(t\right)\le -\lambda exp\left(\kappa t\right)e^T\left(t\right)P_{\mathsf{\Omega }}\left(e\left(t\right)\right)+e^T\left(t\right)e\left(t\right){\mu }_B{\mu }_W+{\displaystyle \sum _{i=1}^{6N_u+6N}}\left|e_i\right|{\mu }_{\xi }+{\displaystyle \sum _{i=1}^{6N_u+6N}}\left|e_i\right|{\mu }_B{\mu }_{\xi }{\mu }_l\hfill \\ =-{\displaystyle \sum _{i=1}^{6N_u+6N}}\left|e_i\right|\left(\lambda exp\left(\kappa t\right)P_{\mathsf{\Omega }}\left(\left|e_i\right|\right)-{\mu }_B{\mu }_W\left|e_i\right|-{\mu }_{\xi }-{\mu }_B{\mu }_{\xi }{\mu }_l\right)\hfill \end{array}$$

(49)

We define the ${\mathsf{\Theta }}_1=\lambda exp\left(\kappa t\right)P_{\mathsf{\Omega }}\left(\left|e_i\right|\right)-{\mu }_B{\mu }_W\left|e_i\right|-{\mu }_{\xi }-{\mu }_B{\mu }_{\xi }{\mu }_l$. We know variable ${\mathsf{\Theta }}_1$ may be positive or negative.

• if ${\mathsf{\Theta }}_1\ge 0$, so $\dot{V}\le 0$. It is obvious that the error variable $e\left(t\right)$ converges to zero from the Lyapunov theorem and the state variable p converges to the optimal solution $p^{\ast }$.
• if ${\mathsf{\Theta }}_1<0$, then $\dot{V}<\delta (\delta >0)$. Therefore, $\dot{V}$ may be positive or negative.

  (a)

  if $\dot{V}\le 0$, we know the error variable $e\left(t\right)$ converges to zero, also the state variable p will converge to optimal solution $p^{\ast }$.

  (b)

  if $\dot{V}>0(0<\dot{V}<\delta )$, and consider the linear activation function ${\mathsf{\Theta }}_1\left(\left|e\left(t\right)\right|\right)=a\left|e\left(t\right)\right|(a\ge 1)$, and $\lambda aexp\left(\kappa t\right)-{\mu }_B{\mu }_W>0$, so we can obtain,

  $$\begin{array}{c}\dot{V}\le -{\displaystyle \sum _{i=1}^{6N_u+6N}}\left|e_i\right|(\lambda exp\left(\kappa t\right)a\left|e_i\right|-{\mu }_B{\mu }_W\left|e_i\right|-{\mu }_{\xi }-{\mu }_B{\mu }_W{\mu }_l\left)\right)\hfill \\ =-\left(\lambda aexp\left(\kappa t\right)-{\mu }_B{\mu }_W\right){\displaystyle \sum _{i=1}^{n+m}}\left|e_i\right|\left(\left|e_i\right|-\frac{{\mu }_{\xi }+{\mu }_B{\mu }_W{\mu }_l}{\lambda aexp\left(\kappa t\right)-{\mu }_B{\mu }_W}\right)\hfill \end{array}$$

  (50)

It is easy to obtain $\lambda >\frac{{\mu }_B{\mu }_W}{a}$. According to (50), in this case, $0<\dot{V}<\delta$, so $V\left(t\right)$ increase that means $\left|e_i\right|$ will increase, and $\dot{V}$ will decrease. Therefore, $\dot{V}$ always exists a moment $\dot{V}\le 0$, then the control system will stabilize again.

It should be noted that when $\dot{v}\left(t\right)=0$, $\left|e_i\right|={\left|e\right|}^{+}$, ${\left|e\right|}^{+}$ is the upper bound. We define ${\mathsf{\Theta }}_2=\frac{{\mu }_{\xi }+{\mu }_B{\mu }_W{\mu }_l}{\lambda aexp\left(\kappa t\right)-{\mu }_B{\mu }_W}$. If $\dot{V}\left(t\right)=0$, ${\displaystyle \sum _{i=1}^{6N_u+6N}}\left|e_i\right|\left(\left|e_i\right|-{\mathsf{\Theta }}_2\right)=0$ and the variable $\left|e_i\right|$ can be seen as the input, the function $\left|e_i\right|\left(\left|e_i\right|-{\mathsf{\Theta }}_2\right)$ will obtain the minimum output, if $\left|e_i\right|=0.5{\mathsf{\Theta }}_2$. Moreover, $\left|e_i\right|\left(\left|e_i\right|-{\mathsf{\Theta }}_2\right)>0$ when $\left|e_i\right|>{\mathsf{\Theta }}_2$

We know ${\displaystyle \sum _{i=1}^{6N_u+6N}}\left|e_i\right|\left(\left|e_i\right|-{\mathsf{\Theta }}_2\right)=0$ and the function mentioned above has a negative minimum output. We assume that $e_j(i=j)$ is the upper bound $e^{+}$, so $e_j$ will be achieved if and only if the rest $6N_u+6N-1$ terms $\left|e_i\right|\left(\left|e_i\right|-{\mathsf{\Theta }}_2\right)$ obtain the minimum point. Therefore,

$$\sum _{i=1}^{6N_u+6N}\left|e_i\right|\left(\left|e_i\right|-{\mathsf{\Theta }}_2\right)=\sum _{i=1,i\ne j}^{6N_u+6N}\left|e_i\right|\left(\left|e_i\right|-{\mathsf{\Theta }}_2\right)+\left|e_j\right|\left(\left|e_j\right|-{\mathsf{\Theta }}_2\right)$$

(51)

Then, ${\mathsf{\Theta }}_2=\frac{{\mu }_{\xi }+{\mu }_B{\mu }_W{\mu }_l}{\lambda aexp\left(\kappa t\right)-{\mu }_B{\mu }_W}$ is substituted in (51) and obtain

$$\begin{array}{c}{\displaystyle \sum _{i=1}^{6N_u+6N}}\left|e_i\right|\left(\left|e_i\right|-{\mathsf{\Theta }}_2\right)={\left|e_j\right|}^2-\left|e_j\right|\left(\frac{{\mu }_{\xi }+{\mu }_B{\mu }_W{\mu }_l}{\lambda aexp\left(\kappa t\right)-{\mu }_B{\mu }_W}\right)-\frac{6N_u+6N-1}{4}{\left(\frac{{\mu }_{\xi }+{\mu }_B{\mu }_W{\mu }_l}{\lambda aexp\left(\kappa t\right)-{\mu }_B{\mu }_W}\right)}^2=0\hfill \end{array}$$

(52)

According to analysis, we can conclude the upper bound $\left|e_j\right|$,

$$\left|e_j\right|=\frac{1}{2}\left((1+\sqrt{6N_u+6N})\eta \right)$$

where $\eta =\frac{{\mu }_{\xi }+{\mu }_B{\mu }_W{\mu }_l}{\lambda aexp\left(\kappa t\right)-{\mu }_B{\mu }_W}$, $\left|e_{end}\right|(end=6N_u+6N)$ converge to zero. □

5. Simulation

In this section, the proposed VPCDNN algorithm is tested with a two-wheel mobile robot in simulation. Two type trajectories will be tested for tracking: circle-shape trajectory (the linear velocity v and angular velocity w are time constant); ’8’-shape trajectory (the linear velocity v and angular velocity w are time varying) and the trajectory function in Cartesian space will be given in Section 5.2. Both of the two trajectories will be tested in mobile robot tracking and the detailed analysis of the simulation results will be given in Section 5.1 and Section 5.2.

5.1. Circle-Shape Tracking

In this part, we control the mobile robot to track the circle-shape trajectory, where the reference linear velocity v and angular velocities $\omega$ are set as $v=0.2$ m/s, $w=0.2$ rad/s, respectively. In order to prove the tracking performance of our method, the initial position of the mobile robot is set as $X_r\left(0\right)=[1.0;0;0]$ which does not coincide with the reference initial point $X\left(0\right)=[1.2;0;0]$.

The parameters of nonlinear MPC are listed as: $N=3$, $N_u=2$, ${\eta }_1=3I$, ${\eta }_2=I$. The neural network parameter $\lambda$ and $\kappa$ are set as: $\lambda =0.1$, $\kappa =1e-4$. The sampling period of MPC is $∆t=0.1s$.

From Figure 4 and Figure 5, it can be seen that the trajectories of the mobile robot quickly converge to the reference trajectory, regardless of the initial position of the mobile robot. From Figure 6 and Figure 7, we can see the mobile robot’s linear and angular velocity are extremely close to the reference values. We can also see that the tracking errors of the mobile robot quickly converge to zero as shown in Figure 8 and Figure 9.

5.2. ‘8’-Shape Tracking

The mobile robot’s task in this scenario is to track the ‘8’-shape trajectory for the periodic movement with changing linear and angular velocity, starting from initialized positions and orientations. The reference trajectory of ‘8’-shape is given as

$$\begin{array}{c}x=sin\left(0.1t\right)\hfill \\ y=2sin\left(0.05t\right)\hfill \end{array}$$

The parameters of nonlinear MPC are listed as: $N=3$, $N_u=2$, ${\eta }_1=3I$, ${\eta }_2=I$. The convergence parameter $\lambda$ and $\kappa$ are set as: $\lambda =0.1$, $\kappa =1e-4$. The sampling period of MPC is set $∆t=0.1$ s.

The initial position of the reference trajectory is $X_r\left(0\right)=[0;0;0.7854]$, and the initial position of the mobile robot is $X\left(0\right)=[0.3;0;0.7]$. From Figure 10 and Figure 11, we can see that the trajectories of the mobile robot quickly converge to the reference trajectory even changing the linear and angular velocity of the robot periodically. From Figure 12 and Figure 13, we can see that the mobile robot’s linear and angular velocity follow the reference velocity quite well. Figure 14 and Figure 15 show that the tracking errors of the mobile robot quickly converge to zero.

In both scenarios, the mobile robot is able to track the reference trajectory successfully. The convergence rate is very fast which is very important for running in real time. The mobile robot solves the trajectory tracking task with the optimal solution from different initialized positions and orientations. The tracking errors also converge very fast to zero. These results show the robustness and effectiveness of the proposed method. Moreover, the overall success rate, fast convergence, and the performance of the proposed VPCDNN based MPC method indicates its strong ability in real time.

In order to verify the effectiveness of VPCDNN, we have added contrast tests using LNN methods which include the time- invariant method called the recurrent neural network (gradient-descent based) in Reference [31]. The comparison results are shown in Figure 16, Figure 17 and Figure 18. We can conclude that both VPCDNN and RNN can track the given trajectory, but VPCDNN has the faster convergence speed which is better for the time-varying optimization problem.

In addition, in this paper, we connect CarSim software to a Matlab/Simulink model and test the VPCDNN-MPC method to track the sinusoidal trajectory which is shown in Figure 19. The requirements of simulator platform were Matlab2018b/Simulink, Carsim2016, Windows10.

6. Conclusions

In this paper, the MPC scheme is presented to the control problem over finite time-horizon with the input and state vector of a mobile robot. In the reformation of the nonlinear affine system of a mobile robot, the nonlinear non-convex optimization problem is converted to the nonlinear convex optimization problem. Therefore, the new network of VPCDNN is also presented to solve the online optimization QP problem of NMPC. In order to test the proposed method, the mobile robot tracks the circle-shape and ’8’-shape trajectories in simulation. The results show that the proposed method successfully tracked the desired trajectories. The method converges quickly to the solution with the VPCDNN. The tracking errors of the VPCDNN based NMPC are maintained at a very low level. All these characteristics enable our proposed method to be a powerful and efficient solution to the control problem of mobile robots.