# A Distributed Model Predictive Control for Multiple Mobile Robots with the Model Uncertainty

Authors: Wei Shang; Hanzong Zhu; Yurong Pan; Xiuhong Li; Daode Zhang
Journal: Discrete Dynamics in Nature and Society
 (2021)

---

## Abstract

In this paper, a distributed predictive control with the model uncertainty which uses the data-driven strategy and robust theory (data-driven RDMPC) is proposed for the formation control of multiple mobile robots. The robust performance objective minimization is applied to replace the quadratic performance objective minimization to establish the optimization problem, where the model uncertainty is considered in the distributed system. The control policy is derived by applying the data-driven strategy, and the future predictive value is obtained by employing the linear law in the historical data. Lyapunov theory is referred to analyze the stability of the mobile robot formation system. The effectiveness of the proposed method is proved by a set of simulation experiments.


---

## Body

## 1. Introduction

In recent years, the robot formation control has become a key technology of multiple mobile robots’ research studies which can be used in the military field, production field, and service field [1–4]. In contrast with a single mobile robot, multiple mobile robots have the advantages of coordination and cooperation, which could improve the efficiency of completing tasks. The multiple mobile robots also enhance the fault tolerance and robustness of the system and strengthen the ability of robots’ environment recognition. Comparing with the leader-follower method [5–8], behavior-based method [9–11],and virtual structure method [12, 13], the model predictive control (MPC) has attracted attention because of its ability of improving the robustness of the system and having better dynamic control performance in solving the problem of consensus protocol in distributed formation control of multiple mobile robots. Also, MPC is convenient to establish the system model.To further improve the performance of MPC, researchers investigate MPC combining with robust control [14, 15], adaptive control [16, 17], and sliding model control [18]. Among them, there are lots of works focusing on overcoming the uncertainty problem of formation control of multiple mobile robots. In the research of multiagent [19], a fault-tolerant controller was designed to compensate the influence of the system uncertainty. By proposing a static nonsmooth event-triggered protocol [20], the influence of the uncertainty is inhibited. Robust and adaptive control algorithms were used [21] for each agent to deal with the system uncertainty associated with leader input and follower agent dynamics. In the study of MAVs [22], adaptive control algorithms were employed locally for each vehicle to solve the uncertainty associated with the system and flying environment. Although these methods mentioned above have disposed of the uncertainty problem of the formation control with MPC, they are not combined with a distributed control. As for the MPC applied in the distributed control (distributed model predictive control (DMPC)), a cooperative DMPC algorithm was reported [23] where the relaxed constraints are used to find the optimization problem including the global optimization objective. This algorithm can be utilized to control the UAV tracking the ground target for real-time path planning [24] and multiple cooperative autonomous ships forming a multiship formation system [25]. To transport the plants and keep the formation at the same time [26], an objective function including trajectory tracking error and formation measurement error was adopted to meet the design requirements based on the DMPC. The DMPC is widely used in the research studies above, but the uncertainty is not considered. To the best of our knowledge, although lots of results have been available, the issue of simultaneous system uncertainty and DMPC has not been fully studied to date.In practical application, the parameters in the model may not be accurate, and it is difficult to get an accurate model of system [27, 28]. Due to the uncertainty of the model, it is almost impossible to predict the future accurately, and the mismatch factors must be noticed. The model uncertainty of the system is ignored in the DMPC algorithm. In addition, the DMPC is too optimistic in the prediction process which would cause the prediction results deviating from the actual effect and lead to many restrictions on the decoupling of the distributed system [29, 30]. For this reason, it is necessary to study the DMPC with system uncertainty. A robust model predictive control approach [31] was suggested to achieve the path tracking of an omnidirectional robot with system uncertainty. The design of robust distributed MPC controller was casted [32] into solving a linear matrix inequality (LMI) optimization problem. A global performance was achieved by using an iterative algorithm. However, the iterative algorithm would cause the jitter fluctuation of the convergence curve and a low computational efficiency [33, 34]. A novel model-free adaptive SMC approach for unknown discrete nonlinear processes was designed [35] by bringing the merits of data-driven control strategy and SMC scheme together, which made the system with model uncertainty of strong robustness. A data-driven hierarchical control (DHC) structure was introduced [36] to ameliorate the performance of systems under the effect of the system uncertainty. The above references demonstrate both robust model predictive control and data-driven strategy which are effective methods to work out the problem of system uncertainty so that it is very important to study the data-driven robust distributed model predictive control with system uncertainty.In this paper, a data-driven robust distributed model predictive control (data-driven RDMPC) with model uncertainty is proposed for multiple mobile robots. The main control method in this paper has these advantages. (1) The data-driven strategy is applied in this method using the historical data representing the future predictive value so that data processing is not limited to the model and the conservatism of system can be reduced. (2) The robust performance objective not the quadratic performance objective is chosen in this method to solve the optimization problem. As a result, the system can consider the model uncertainty, the fluctuation of the convergence curve can be decreased, and also the efficiency and accuracy of the algorithm can be promoted.This paper is arranged as follows. The control problem is described in Section2. In Section 3, a data-driven distributed model predictive controller is designed with considering the system uncertainty, and the stability of the system is analyzed. Section 4 presents a simulation example to illustrate the effectiveness of the method. Section 5 puts forward some conclusions.


## 2. Problem Formation

### 2.1. Notations and Graph Theory

Throughout the paper, the following notations will be used. The functionXk+Δk|k means the predicted value of variables Xk+Δk in time k. Xk|k=Xk represents the variable x in current time and can be measured by sensors. And a quadric form function is given as xQ2=xTQx where Q>0.This paper examines the mobile robot formation control with the distributed control method. A finite index is set to mark the mobile robot, making each robot node as 1, 2,...,n. A directed and fixed graph G is V,ε,A. V=v1,v2,…,vn is a set of nodes, where vi is referred to the i-th mobile robot. ε∈V×V is the edge set. εij=vi,vj∈ε indicates that the i-th robot can obtain the information from j-th robot. A=aij is the weighted adjacency matrix of G. If εij∈ε, we have aij=aji>0 while aij=0 and aii=0.


### 2.2. Mathematical Model

The mobile roboti has qi=xiyiαiTon their path Γi, i = 1, 2, …, n. The statement of the robot j communicating with the robot i is qj=xjyjαjT, and xi,yi,αi and xj,yj,αj are the positionx,y and orientation α of robot i and j, respectively. The desired relative statement of robot i and j is pij=xijryijrαijrTwhere xijr, yijr, and αijr represent the desired relative position and relative orientation. The kinematic models of i and j are established as follows:(1)x˙iy˙iα˙i=cosαi0sinαi001viωi,(2)x˙jy˙jα˙j=cosαj0sinαj001vjωj,where vi and ωi represent the control input of the forward and angular velocities of robot i and vj and ωj represent the control input of the forward and angular velocities of robot j, and then the relative statement between them can be expressed as follows:(3)xijyijαij=cosαisinαi0−sinαicosαi0001qi−qj−pij.By taking the derivative of (3) and combining with (1) and (2), the following system is obtained:(4)x˙ij=ωiyij−yijr−vi+vjcosαij−αijr,y˙ij=−ωixij−xijr+vjsinαij−αijr,α˙ij=ωi−ωj−ωijr,where ωijr is the desired angular speed error between robot i and robot j. System (4) between robot i and robot j can be linearized around Xij=0,Uij=0 as follows:(5)X˙ij=AliXij+BliUij,where Ali=0ωi0−ωi0vj000, Bli=100001, and Xij=xij−xijryij−yijrαij−αijrT are the statement and Uij=−vi+vjcosαij−αijrωi−ωj−ωijrT is the control input. Considering the communication between robot i and other robots, the linearized system of robot i is written as follows:(6)X˙i=AliXi+BliUi,where Xi=∑s=1maijXij and Ui=∑s=1maijUij, and the elements of Xi are defined as Xi=xBiyBiαBiT; therefore, xBi=∑s=1maijxij−xijr, yBi=∑s=1maijyij−yijr, and αBi=∑s=1maijαij−αijr. It should be noted that Xi is the sum of relative position errors between each mobile robot i and robot j and m is the number of robots communicating with robot i. After discretizing equation (6), the normative system of mobile robot i is as follows:(7)Xinormk+1=A0iXik+B0iUik,where A0i=eAliTandB0i=∫0TeAliτdτBli, T is the sampling period. Because of the uncertainty, the system will produce errors. Besides, the linearization process will also bring the errors to the system. To deal with these two errors concurrently, an uncertainty term ΔXik related to control input Uik and state variable Xik is considered to describe the uncertainty of the system as follows:(8)ΔXik=δiΔCiXik+ΔDiUik,where ΔCi and ΔDi are known matrices. δi is an unknown but bounded time-independent matrix and satisfies δi2≡σδi≤1. Then, considering the external disturbance of the mobile robot, the bounded disturbance term dik∈R3×1 is added. The uncertain but bounded system model of mobile robot i is shown as follows:(9)xik+1=xinormk+1+Δxik+dik=A0i+δiΔCiXik+B0i+δiΔDiUik+dik,and the parameter of this uncertain system model is AiBi=A0i+δiΔCiB0i+δiΔDi∈Ω.In the control process, the actuator usually bears limited control force. Therefore, it is very necessary to consider the saturation constraints of the input and the hard constraints of the state variables in the controller design [37]. Otherwise, the control will be inaccurate because of the control input and state variables exceed the limits, and this situation will inevitably reduce the tracking performance of the control system [38]. Therefore, saturation constraint (10) to limit the control input and hard constraint (11) to limit the state variable are added as follows:(10)Uikc≤Uikc,max,c=1,2,(11)FicXik≤Xikc,max,c=1,2,3,where Uikc,max and Xikc,max are the maximum of c-th element of Uik and Xik and Fic is the c-th row of the assigned matrix Fi∈R3×3.


## 2.1. Notations and Graph Theory

Throughout the paper, the following notations will be used. The functionXk+Δk|k means the predicted value of variables Xk+Δk in time k. Xk|k=Xk represents the variable x in current time and can be measured by sensors. And a quadric form function is given as xQ2=xTQx where Q>0.This paper examines the mobile robot formation control with the distributed control method. A finite index is set to mark the mobile robot, making each robot node as 1, 2,...,n. A directed and fixed graph G is V,ε,A. V=v1,v2,…,vn is a set of nodes, where vi is referred to the i-th mobile robot. ε∈V×V is the edge set. εij=vi,vj∈ε indicates that the i-th robot can obtain the information from j-th robot. A=aij is the weighted adjacency matrix of G. If εij∈ε, we have aij=aji>0 while aij=0 and aii=0.


## 2.2. Mathematical Model

The mobile roboti has qi=xiyiαiTon their path Γi, i = 1, 2, …, n. The statement of the robot j communicating with the robot i is qj=xjyjαjT, and xi,yi,αi and xj,yj,αj are the positionx,y and orientation α of robot i and j, respectively. The desired relative statement of robot i and j is pij=xijryijrαijrTwhere xijr, yijr, and αijr represent the desired relative position and relative orientation. The kinematic models of i and j are established as follows:(1)x˙iy˙iα˙i=cosαi0sinαi001viωi,(2)x˙jy˙jα˙j=cosαj0sinαj001vjωj,where vi and ωi represent the control input of the forward and angular velocities of robot i and vj and ωj represent the control input of the forward and angular velocities of robot j, and then the relative statement between them can be expressed as follows:(3)xijyijαij=cosαisinαi0−sinαicosαi0001qi−qj−pij.By taking the derivative of (3) and combining with (1) and (2), the following system is obtained:(4)x˙ij=ωiyij−yijr−vi+vjcosαij−αijr,y˙ij=−ωixij−xijr+vjsinαij−αijr,α˙ij=ωi−ωj−ωijr,where ωijr is the desired angular speed error between robot i and robot j. System (4) between robot i and robot j can be linearized around Xij=0,Uij=0 as follows:(5)X˙ij=AliXij+BliUij,where Ali=0ωi0−ωi0vj000, Bli=100001, and Xij=xij−xijryij−yijrαij−αijrT are the statement and Uij=−vi+vjcosαij−αijrωi−ωj−ωijrT is the control input. Considering the communication between robot i and other robots, the linearized system of robot i is written as follows:(6)X˙i=AliXi+BliUi,where Xi=∑s=1maijXij and Ui=∑s=1maijUij, and the elements of Xi are defined as Xi=xBiyBiαBiT; therefore, xBi=∑s=1maijxij−xijr, yBi=∑s=1maijyij−yijr, and αBi=∑s=1maijαij−αijr. It should be noted that Xi is the sum of relative position errors between each mobile robot i and robot j and m is the number of robots communicating with robot i. After discretizing equation (6), the normative system of mobile robot i is as follows:(7)Xinormk+1=A0iXik+B0iUik,where A0i=eAliTandB0i=∫0TeAliτdτBli, T is the sampling period. Because of the uncertainty, the system will produce errors. Besides, the linearization process will also bring the errors to the system. To deal with these two errors concurrently, an uncertainty term ΔXik related to control input Uik and state variable Xik is considered to describe the uncertainty of the system as follows:(8)ΔXik=δiΔCiXik+ΔDiUik,where ΔCi and ΔDi are known matrices. δi is an unknown but bounded time-independent matrix and satisfies δi2≡σδi≤1. Then, considering the external disturbance of the mobile robot, the bounded disturbance term dik∈R3×1 is added. The uncertain but bounded system model of mobile robot i is shown as follows:(9)xik+1=xinormk+1+Δxik+dik=A0i+δiΔCiXik+B0i+δiΔDiUik+dik,and the parameter of this uncertain system model is AiBi=A0i+δiΔCiB0i+δiΔDi∈Ω.In the control process, the actuator usually bears limited control force. Therefore, it is very necessary to consider the saturation constraints of the input and the hard constraints of the state variables in the controller design [37]. Otherwise, the control will be inaccurate because of the control input and state variables exceed the limits, and this situation will inevitably reduce the tracking performance of the control system [38]. Therefore, saturation constraint (10) to limit the control input and hard constraint (11) to limit the state variable are added as follows:(10)Uikc≤Uikc,max,c=1,2,(11)FicXik≤Xikc,max,c=1,2,3,where Uikc,max and Xikc,max are the maximum of c-th element of Uik and Xik and Fic is the c-th row of the assigned matrix Fi∈R3×3.


## 3. Controller Design and Stability Analysis

### 3.1. Controller Design

It is worth emphasizing that each mobile robot has its own independent MPC controller. By communicating with each other, these robots receive the information which can form the coupling conditions to solve the optimal control problem online on their respective controllers. For each robot MPC controller, in the absence of an accurate model, the performance of solving the optimization problem which uses the method of quadratic performance objective minimization will be very poor, and the computational load will be large. So, the method applying the robust performance objective minimization to replace the quadratic performance objective minimization to establish the optimization problem is considered for the controller design. A feedback gainK is set to adjust Uk and is continuously updated by the optimization problem. The optimization problem is shown as follows:(12)minUimaxAiBi∈ΩJi0∞ks.t.6−10.In order to propose an appropriate objective function, a trajectory parameter is needed to be included in the objective function of the normative system to construct the coupling relationship. A parameterζi is defined to describe the trajectory Γiζi=xΓiζiyΓiζiαΓiζi of mobile robot i, which is an iterative function related to control input, forward, and angular velocities:(13)ζik+1=ζik+MiUik+NiviωiT,where Mi and Ni are the given parameter matrix. Defining ζijr as the desired value of the trajectory parameter difference between robot i and robot j, it is an important parameter to predict the trajectory of robot j through robot i and then the trajectory of robot j can be described as Γjζi=xΓiζi+xΓijζijryΓiζi+yΓijζijrαΓiζi+αΓijζijr. The information which robots interact includes their position states and trajectory parameters. The position states and the trajectory parameters are used to calculate Xij and form coupling conditions, respectively. The objective function of the normative system coupled with ζj is as follows:(14)Ji0∞k=∑Δk=0∞Xik+Δk+1|kQ¯i2+Uik+Δk|kR¯i2+∑s=1mζik+Δk|k−ζjk+Δk|k−ζijrW¯i,l2=∑Δk=0∞AiXik+Δk|k+BiUik+Δk|kQ¯i2+Uik+Δk|kR¯i2+∑s=1mHiUik+Δk|k−ψijW¯ij2,where ψij=Θjk+Δk|k+HiUjk+Δk|k+v¯j+Θijr−Θik+Δk|k−Hiv¯i,Q¯i=diagQi,Qi,Qi, Qi>0, R¯i=diagRi,Ri, Ri>0, W¯ij=diagwij,wij,wij, wij>0, Θi=ζi,ζi,ζiT, Θj=ζj,ζj,ζjT, Θijr=ζijr,ζijr,ζijrT, v¯i=viωiT, v¯j=vjωjT, and Hi=MiMiMiT.For the uncertain system (9), compared with Xik|k, we focus on Xik+1|k to control the system and get Xik+1|k from the historical data without uncertain factors. Xik+1|k should be limited to the feasible region Z=z|zTQ˜i−1z≤1 of the optimization problem, so there is XiTk+1|kQ˜i−1Xik+1|k≤1, and using the Schur complement, we can get the following:(15)1∗Xik+1|kQ˜i≥0.If the ellipsoidZ is an invariant set, then the objective function can be divided into two parts as follows:(16)Ji0∞k=A0iXik|k+B0iUik|kQ¯i2+Uik|kR¯i2+∑s=1mHiUik|k−ψijW¯ij2+∑Δk=1∞A0iXik+Δk|k+B0iUik+Δk|kQ¯i2+Uik+Δk|kR¯i2+∑s=1mHiUik+Δk|kψijW¯ij2=Ji0k+Ji1∞k.The optimization problem is transformed into solving a set of state feedback lawsUi=KiXi. In the case, the controller objective minimizes the upper bound of the worst-case objective function. Let γ1i and γ2i be the upper bound of objective function Ji0k and Ji1∞k, respectively.Theorem 1.
For the uncertain system (9), if there exists Λ1=diagλ1,λ2,λ3>0, then the system will be stable and Z is an invariant set when following conditions are satisfied:(17)Q˜i∗∗∗∗A0iQ˜i+B0iYiQ˜i∗∗∗ΔCiQ˜i+ΔDiYiA0iTQ¯iA0i0.5Q˜iSi0.5Yi000Λ1∗∗0γ2iI∗00γ2iI≥0,(18)Q˜i∗diTkγ2i0.5∑s=1mψijTW¯ijψij≥0,(19)Gi∗YiQ˜i≥0,Gicc≤Uic2,c=1,2,(20)Li∗Q˜iFiTQ˜i≥0,Licc≤Xic2,c=1,2,3,where S=B0iTQ¯iB0i+R¯i+∑s=1mHiTW¯ijHi, Gicc and Licc represent the c-th row and c-th column of the matrix, and Uic and Xic represent the c-th element of the control input and state variable.Proof.
of Theorem 1. IfXik+1|k is restricted in the ellipsoid Z and Z is an invariant set, then the objective function in infinite time domain can be regarded as two parts of (16), where Ji01k is finite time domain. And the stability of the system can be analyzed by aiming at Ji1∞k which is in infinite time domain. Consider a Lyapunov quadratic function VXi=XiTPiXi, Pi>0. For any time, Ai and Bi satisfy AiBi∈Ω. According to the Lyapunov stability criterion, VXi should meet the following inequality:(21)VXik+j+1|k−VXik+j|k≤A0iXik+Δk|k+B0iUik+Δk|kQ¯i2+Uik+Δk|kR¯i2+∑s=1mHiUik+j|k−ψi,lW¯i,l2.
Summing (21) from i=1 to i=∞ and requiring Xi∞|k=0 or VXi∞|k=0, the result is as follows:(22)maxA0iB0i∈ΩJi1∞k≤VXik+1|k=XiTk+1|kPiXik+1|k≤γ2i.
There isXiTk+1|kQ˜i−1Xik+1|k≤1, then let Pi=γ2iQ˜i−1. Substituting (9) and state feedback law Ui=KiXi into (21), we can get the following:(23)XiTk+Δk|kAi+BiKiTPiAi+BiKi−Pi+Gi+BiKiTQ¯iGi+BiKi+KiTR¯iKi+∑s=1mHiKiTW¯ijHiKiXik+Δk|k+diTkPidik−∑s=1mψijTW¯ijψij≤0,where (23) can be regarded as two parts, and they are as follows:(23a)Ai+BiKiTPiAi+BiKi−Pi+Gi+BiKiTQ¯iGi+BiKi+KiTR¯iKi+∑s=1mHiKiTW¯ijHiKi≤0,(23b)diTkPidik−∑s=1mψijTW¯ijψij≤0.
In the above formula, the feedback gain is defined asKi=YiQ˜i−1, and the solution of Yi is given by the optimization problem. Substituting AiBi=A0i+δiΔCiB0i+δiΔDi, feedback gain Ki=YiQ˜i−1 and Pi=γ2iQ˜i−1 into (23a) and (23b), and transforming (23a) and (23b) into matrix inequality form, the results are (19) and (20), respectively. Systems (19) and (20) guarantee the positive invariance of Z. In the invariant set Z, the saturation constraint (10) and hard constraint (11) are expressed as (19) and (20) in matrix form, respectively, and the control input is limited by Gi, and the state variable is limited by Li.Remark 1.
In Theorem1, equation (17) can get the constraint that the upper bound γ2i of Ji1∞k is satisfied when the system is stable. After the upper bound γ1i of Ji0k can also be expressed as some concrete constraints, the problem of minimizing the upper bound of the objective function in the worst case can be formulated and solved. And with the update of feedback gain, Xik is always limited in the ellipsoid Z at any time in the optimization process.Theorem 2.
The upper boundγ1i of Ji0k satisfies the following relation:(24)γ1i∗∗Uik|kSi≥0,where γ1i=xiTk|kA0iTQ¯iA0ixik|k−∑s=1mψijTW¯ijψij+γ1i∗.Proof.
of Theorem 2. The objective function ofJi0k is as follows:(25)Ji0k=A0iXik|k+B0iUik|kQ¯i2+Uik|kR¯i2+∑s=1mHiuik|k−ψijW¯ij2=xiTk|kA0iTQ¯iA0ixik|k+uiTSiui−∑s=1mψijTW¯ijψij≤γ1i,where xik|k can be measured by sensors and exists in historical data. ψij is the function related to ζi,ζj,anduj. ζi can also be measured by sensors, and ζjanduj can interact with mobile robot j through wireless network or Bluetooth technology. Let uiTSiui≤γ1i∗, equation (24) is obtained after Schur complement.
Theorem2 gives the condition of upper bound γ1i of Ji0k; then, the upper bound of the objective function in the infinite time in the worst case can be expressed as follows:(26)maxA0iB0i∈ΩJi0∞k≤γ1i+γ2i=xiTk|kA0iTQ¯iA0ixik|k−∑s=1mψijTW¯ijψij+γ1i∗+γ2i.
The conditions about the upper bound of the objective function, feedback gain, and other parameters are given in the above theorems with the stable system, which can be used as the constraints of the optimization problem. Then, the optimization problem (12) can be described as LMI problem as follows:(27)minQ˜i,Yi,Gi,Li,Λ1,γ1i,γ2iγ1i∗+γ2is.t.15,17,18,19,20,24.
As mentioned in the paper above, we focus onXik+1|k not Xik|k to control the system when solving the optimization problem. At this time, the difficulty is how to use the historical data from the previous time to represent the predicted value of the future time. So, consider the following linear law of historical data:(28)Xik|k=Eik1℧i21+Eik2℧i22+ΔXik,(29)Uik|k=Eik1℧i31+Eik2℧i32+ΔUik,where Eik∈R1×2 represents the coefficient matrix when the historical data are employed to indicate the future forecast value linearly. ℧i2=Xik−1|kXik−2|kT and ℧i3=Uik−1|kUik−2|kT denote storage units containing two sets of historical data of Xik|k and Uik|k, respectively. When k=0, Xik−2|k=0 and Uik−2|k=0. ΔXik and ΔUik are the remainders of the errors. Define ℧i4=Xik|kXik−1|kT as the storage unit including two sets of historical data of Xik+1|k, but this unit has no direct connection with ℧i2 and ℧i3, and the date of them exists independently. The system does not take the bounded external interference into consideration; after substituting (28) and (29) into the system model, the result is as follows:(30)Xik+1|k=Eik1℧i41,:+Eik2℧i42,:+A0iΔXik+B0iΔUik+δiΔCiXk+ΔDiΔUik.
Substituting (30) into (15), we can get(31)1∗Eik1℧i41,:+Eik2℧i42,:+A0iΔXik+B0iΔUik+δiΔCiXk+ΔDiΔUikQ˜i≥0.
After expanding (31), set Λ2=diagλ4,λ5,λ6>0, and the result is shown as follows:(32)1∗∗Eik1℧i41,:+Eik2℧i42,:+A0iΔXik+B0iΔUikQ˜i0ΔCiΔXik+ΔDiΔUik0Λ2≥0.
Substituting (29) into (24), we can get(33)γ1i∗∗Eik1℧i31,:+Eik2℧i32,:+ΔUikSi≥0.
Inequalities (32) and (33) are obtained from historical data, and then the data-driven optimization problem (27) can be expressed as follows:(34)minQ˜i,Yi,Gi,LiΛ1Λ1,Eik,ΔXik,ΔUik,γ1i′,γ2iγ1i∗+γ2is.t17,18,19,20,32,33.


### 3.2. Stability Analysis

Lemma 1.
When the system is inspired, the stored energy decreases with time and finally reaches a minimum, and then the system is gradually stable. If there is a positive definite scalar functionVXi with negative definite V˙Xi, the system is stable.
A positive definite Lyapunov functionVXi=XiTPiXi is defined in 3.1, where Pi>0. If V˙Xi is negative definite, the system is asymptotically stable. The condition for V˙Xi to be negative definite is given in equation (21). By using the robust theory, the analysis of Xi is transformed into the analysis of state feedback law Ui=KiXi, and the conditions (17) and (18) for state feedback gain Ki are derived. According to Lemma 1, when (17) and (18) are held, the system will tend to stable.


## 3.1. Controller Design

It is worth emphasizing that each mobile robot has its own independent MPC controller. By communicating with each other, these robots receive the information which can form the coupling conditions to solve the optimal control problem online on their respective controllers. For each robot MPC controller, in the absence of an accurate model, the performance of solving the optimization problem which uses the method of quadratic performance objective minimization will be very poor, and the computational load will be large. So, the method applying the robust performance objective minimization to replace the quadratic performance objective minimization to establish the optimization problem is considered for the controller design. A feedback gainK is set to adjust Uk and is continuously updated by the optimization problem. The optimization problem is shown as follows:(12)minUimaxAiBi∈ΩJi0∞ks.t.6−10.In order to propose an appropriate objective function, a trajectory parameter is needed to be included in the objective function of the normative system to construct the coupling relationship. A parameterζi is defined to describe the trajectory Γiζi=xΓiζiyΓiζiαΓiζi of mobile robot i, which is an iterative function related to control input, forward, and angular velocities:(13)ζik+1=ζik+MiUik+NiviωiT,where Mi and Ni are the given parameter matrix. Defining ζijr as the desired value of the trajectory parameter difference between robot i and robot j, it is an important parameter to predict the trajectory of robot j through robot i and then the trajectory of robot j can be described as Γjζi=xΓiζi+xΓijζijryΓiζi+yΓijζijrαΓiζi+αΓijζijr. The information which robots interact includes their position states and trajectory parameters. The position states and the trajectory parameters are used to calculate Xij and form coupling conditions, respectively. The objective function of the normative system coupled with ζj is as follows:(14)Ji0∞k=∑Δk=0∞Xik+Δk+1|kQ¯i2+Uik+Δk|kR¯i2+∑s=1mζik+Δk|k−ζjk+Δk|k−ζijrW¯i,l2=∑Δk=0∞AiXik+Δk|k+BiUik+Δk|kQ¯i2+Uik+Δk|kR¯i2+∑s=1mHiUik+Δk|k−ψijW¯ij2,where ψij=Θjk+Δk|k+HiUjk+Δk|k+v¯j+Θijr−Θik+Δk|k−Hiv¯i,Q¯i=diagQi,Qi,Qi, Qi>0, R¯i=diagRi,Ri, Ri>0, W¯ij=diagwij,wij,wij, wij>0, Θi=ζi,ζi,ζiT, Θj=ζj,ζj,ζjT, Θijr=ζijr,ζijr,ζijrT, v¯i=viωiT, v¯j=vjωjT, and Hi=MiMiMiT.For the uncertain system (9), compared with Xik|k, we focus on Xik+1|k to control the system and get Xik+1|k from the historical data without uncertain factors. Xik+1|k should be limited to the feasible region Z=z|zTQ˜i−1z≤1 of the optimization problem, so there is XiTk+1|kQ˜i−1Xik+1|k≤1, and using the Schur complement, we can get the following:(15)1∗Xik+1|kQ˜i≥0.If the ellipsoidZ is an invariant set, then the objective function can be divided into two parts as follows:(16)Ji0∞k=A0iXik|k+B0iUik|kQ¯i2+Uik|kR¯i2+∑s=1mHiUik|k−ψijW¯ij2+∑Δk=1∞A0iXik+Δk|k+B0iUik+Δk|kQ¯i2+Uik+Δk|kR¯i2+∑s=1mHiUik+Δk|kψijW¯ij2=Ji0k+Ji1∞k.The optimization problem is transformed into solving a set of state feedback lawsUi=KiXi. In the case, the controller objective minimizes the upper bound of the worst-case objective function. Let γ1i and γ2i be the upper bound of objective function Ji0k and Ji1∞k, respectively.Theorem 1.
For the uncertain system (9), if there exists Λ1=diagλ1,λ2,λ3>0, then the system will be stable and Z is an invariant set when following conditions are satisfied:(17)Q˜i∗∗∗∗A0iQ˜i+B0iYiQ˜i∗∗∗ΔCiQ˜i+ΔDiYiA0iTQ¯iA0i0.5Q˜iSi0.5Yi000Λ1∗∗0γ2iI∗00γ2iI≥0,(18)Q˜i∗diTkγ2i0.5∑s=1mψijTW¯ijψij≥0,(19)Gi∗YiQ˜i≥0,Gicc≤Uic2,c=1,2,(20)Li∗Q˜iFiTQ˜i≥0,Licc≤Xic2,c=1,2,3,where S=B0iTQ¯iB0i+R¯i+∑s=1mHiTW¯ijHi, Gicc and Licc represent the c-th row and c-th column of the matrix, and Uic and Xic represent the c-th element of the control input and state variable.Proof.
of Theorem 1. IfXik+1|k is restricted in the ellipsoid Z and Z is an invariant set, then the objective function in infinite time domain can be regarded as two parts of (16), where Ji01k is finite time domain. And the stability of the system can be analyzed by aiming at Ji1∞k which is in infinite time domain. Consider a Lyapunov quadratic function VXi=XiTPiXi, Pi>0. For any time, Ai and Bi satisfy AiBi∈Ω. According to the Lyapunov stability criterion, VXi should meet the following inequality:(21)VXik+j+1|k−VXik+j|k≤A0iXik+Δk|k+B0iUik+Δk|kQ¯i2+Uik+Δk|kR¯i2+∑s=1mHiUik+j|k−ψi,lW¯i,l2.
Summing (21) from i=1 to i=∞ and requiring Xi∞|k=0 or VXi∞|k=0, the result is as follows:(22)maxA0iB0i∈ΩJi1∞k≤VXik+1|k=XiTk+1|kPiXik+1|k≤γ2i.
There isXiTk+1|kQ˜i−1Xik+1|k≤1, then let Pi=γ2iQ˜i−1. Substituting (9) and state feedback law Ui=KiXi into (21), we can get the following:(23)XiTk+Δk|kAi+BiKiTPiAi+BiKi−Pi+Gi+BiKiTQ¯iGi+BiKi+KiTR¯iKi+∑s=1mHiKiTW¯ijHiKiXik+Δk|k+diTkPidik−∑s=1mψijTW¯ijψij≤0,where (23) can be regarded as two parts, and they are as follows:(23a)Ai+BiKiTPiAi+BiKi−Pi+Gi+BiKiTQ¯iGi+BiKi+KiTR¯iKi+∑s=1mHiKiTW¯ijHiKi≤0,(23b)diTkPidik−∑s=1mψijTW¯ijψij≤0.
In the above formula, the feedback gain is defined asKi=YiQ˜i−1, and the solution of Yi is given by the optimization problem. Substituting AiBi=A0i+δiΔCiB0i+δiΔDi, feedback gain Ki=YiQ˜i−1 and Pi=γ2iQ˜i−1 into (23a) and (23b), and transforming (23a) and (23b) into matrix inequality form, the results are (19) and (20), respectively. Systems (19) and (20) guarantee the positive invariance of Z. In the invariant set Z, the saturation constraint (10) and hard constraint (11) are expressed as (19) and (20) in matrix form, respectively, and the control input is limited by Gi, and the state variable is limited by Li.Remark 1.
In Theorem1, equation (17) can get the constraint that the upper bound γ2i of Ji1∞k is satisfied when the system is stable. After the upper bound γ1i of Ji0k can also be expressed as some concrete constraints, the problem of minimizing the upper bound of the objective function in the worst case can be formulated and solved. And with the update of feedback gain, Xik is always limited in the ellipsoid Z at any time in the optimization process.Theorem 2.
The upper boundγ1i of Ji0k satisfies the following relation:(24)γ1i∗∗Uik|kSi≥0,where γ1i=xiTk|kA0iTQ¯iA0ixik|k−∑s=1mψijTW¯ijψij+γ1i∗.Proof.
of Theorem 2. The objective function ofJi0k is as follows:(25)Ji0k=A0iXik|k+B0iUik|kQ¯i2+Uik|kR¯i2+∑s=1mHiuik|k−ψijW¯ij2=xiTk|kA0iTQ¯iA0ixik|k+uiTSiui−∑s=1mψijTW¯ijψij≤γ1i,where xik|k can be measured by sensors and exists in historical data. ψij is the function related to ζi,ζj,anduj. ζi can also be measured by sensors, and ζjanduj can interact with mobile robot j through wireless network or Bluetooth technology. Let uiTSiui≤γ1i∗, equation (24) is obtained after Schur complement.
Theorem2 gives the condition of upper bound γ1i of Ji0k; then, the upper bound of the objective function in the infinite time in the worst case can be expressed as follows:(26)maxA0iB0i∈ΩJi0∞k≤γ1i+γ2i=xiTk|kA0iTQ¯iA0ixik|k−∑s=1mψijTW¯ijψij+γ1i∗+γ2i.
The conditions about the upper bound of the objective function, feedback gain, and other parameters are given in the above theorems with the stable system, which can be used as the constraints of the optimization problem. Then, the optimization problem (12) can be described as LMI problem as follows:(27)minQ˜i,Yi,Gi,Li,Λ1,γ1i,γ2iγ1i∗+γ2is.t.15,17,18,19,20,24.
As mentioned in the paper above, we focus onXik+1|k not Xik|k to control the system when solving the optimization problem. At this time, the difficulty is how to use the historical data from the previous time to represent the predicted value of the future time. So, consider the following linear law of historical data:(28)Xik|k=Eik1℧i21+Eik2℧i22+ΔXik,(29)Uik|k=Eik1℧i31+Eik2℧i32+ΔUik,where Eik∈R1×2 represents the coefficient matrix when the historical data are employed to indicate the future forecast value linearly. ℧i2=Xik−1|kXik−2|kT and ℧i3=Uik−1|kUik−2|kT denote storage units containing two sets of historical data of Xik|k and Uik|k, respectively. When k=0, Xik−2|k=0 and Uik−2|k=0. ΔXik and ΔUik are the remainders of the errors. Define ℧i4=Xik|kXik−1|kT as the storage unit including two sets of historical data of Xik+1|k, but this unit has no direct connection with ℧i2 and ℧i3, and the date of them exists independently. The system does not take the bounded external interference into consideration; after substituting (28) and (29) into the system model, the result is as follows:(30)Xik+1|k=Eik1℧i41,:+Eik2℧i42,:+A0iΔXik+B0iΔUik+δiΔCiXk+ΔDiΔUik.
Substituting (30) into (15), we can get(31)1∗Eik1℧i41,:+Eik2℧i42,:+A0iΔXik+B0iΔUik+δiΔCiXk+ΔDiΔUikQ˜i≥0.
After expanding (31), set Λ2=diagλ4,λ5,λ6>0, and the result is shown as follows:(32)1∗∗Eik1℧i41,:+Eik2℧i42,:+A0iΔXik+B0iΔUikQ˜i0ΔCiΔXik+ΔDiΔUik0Λ2≥0.
Substituting (29) into (24), we can get(33)γ1i∗∗Eik1℧i31,:+Eik2℧i32,:+ΔUikSi≥0.
Inequalities (32) and (33) are obtained from historical data, and then the data-driven optimization problem (27) can be expressed as follows:(34)minQ˜i,Yi,Gi,LiΛ1Λ1,Eik,ΔXik,ΔUik,γ1i′,γ2iγ1i∗+γ2is.t17,18,19,20,32,33.


## 3.2. Stability Analysis

Lemma 1.
When the system is inspired, the stored energy decreases with time and finally reaches a minimum, and then the system is gradually stable. If there is a positive definite scalar functionVXi with negative definite V˙Xi, the system is stable.
A positive definite Lyapunov functionVXi=XiTPiXi is defined in 3.1, where Pi>0. If V˙Xi is negative definite, the system is asymptotically stable. The condition for V˙Xi to be negative definite is given in equation (21). By using the robust theory, the analysis of Xi is transformed into the analysis of state feedback law Ui=KiXi, and the conditions (17) and (18) for state feedback gain Ki are derived. According to Lemma 1, when (17) and (18) are held, the system will tend to stable.


## 4. Results

In this part, two simulations with a group of mobile robots are conducted to verify the proposed approach.


### 4.1. Simulation 1

In this simulation, the lane-change trajectory is considered. The desired trajectory is considered as follows:(35)Γ=xΓr=−ζ,yΓr=−21+e−2ζ−2.Letζ12r=−10,ζ23r=−10,ζ34r=15,ζ45r=−10,andζ51r=15. The trajectory error of every two mobile robots can be expressed as xΓijζijr=−ζijr and yΓij=2. The communication topology is shown in Figure 1. The parameters of robots are given as follows: Q1=100, S1=15, w1,l=4, Q2=100, S2=6, w2,l=4, Q3=120, S3=9.6, w3,l=4, Q4=100, S4=6, w4,l=4, Q5=120, S5=9.6, and w5,l=4. Set the sampling period T = 0.1 s and coefficients ΔCi=0,ΔCi=Bli, and Mi=Ni=−0.20.1. Each robot’s angle speed and line speed are ω=1 and v=1, respectively. Taking the uncertainty set as δi=diag−1,−1,−1, the simulation is carried out for 50 s. The initial state values of each mobile robot are as follows:(36)X10=0.25−0.050.1T,ζ10=10,X20=0.10.030.15T,ζ20=20,X30=0.15−0.020.1T,ζ30=30,X40=0.15−0.020.1T,ζ40=15,X50=0.15−0.020.1T,ζ50=25.Figure 1
Communication topology: (a) Robot 1; (b) Robot 2; (c) Robot 3; (d) Robot 4; (e) Robot 5.The bounded external disturbance of each mobile robot is as follows:(37)d1k=0.22k+1+21000cos0.1k0.2k+1+21000sin0.1k0.1k+1+11000cos0.2k,d2k=0.12k+1+11000cos0.1k0.2k+1+21000sin0.1k0.12k+1+11000cos0.2k,d3k=0.12k+1+11000cos0.1k0.15k+1+11000sin0.1k0.2k+1+21000cos0.2k,d4k=0.22k+1+11000cos0.1k0.15k+1+21000sin0.1k0.1k+1+11000cos0.2k,d5k=0.12k+1+11000cos0.1k0.2k+1+21000sin0.1k0.1k+1+21000cos0.2k.The simulation is compared with the traditional DMPC. Figure2 shows the state error Xi of each mobile robot. Figure 3 illustrates the absolute value of state error Xi of each mobile robot. Figure 4 displays ζij and ζijr of each mobile robot. Figure 5 presents the different values between ζij and ζijr of each mobile robot. Figure 6 demonstrates the trajectories of five mobile robots.Figure 2
The state errorXi of each mobile robot: (a) Robot 1; (b) Robot 2; (c) Robot 3; (d) Robot 4; (e) Robot 5.
(a)(b)(c)(d)(e)Figure 3
The absolute value of state errorXi of each mobile robot: (a) Robot 1; (b) Robot 2; (c) Robot 3; (d) Robot 4; (e) Robot 5.
(a)(b)(c)(d)(e)Figure 4
The trajectory parameter differenceζij and the desired trajectory parameter difference ζijr of each mobile robot: (a) Robot 1; (b) Robot 2; (c) Robot 3; (d) Robot 4; (e) Robot 5.
(a)(b)(c)(d)(e)Figure 5
The different values betweenζij and ζijr of each mobile robot: (a) Robot 1; (b) Robot 2; (c) Robot 3; (d) Robot 4; (e) Robot 5.
(a)(b)(c)(d)(e)Figure 6
Trajectories of five mobile robots.The result in Figure2 shows that compared with the traditional DMPC, the optimization process of data-driven RDMPC uses the optimal solution instead of the iterative algorithm, which is more efficient and greatly reduces the fluctuation of the convergence curve. In Figure 3, the simulation results with the control target Xi=000T indicate that the simulation results of data-driven RDMPC are always smaller than those of the traditional DMPC in the process of gradual stabilization, which highlights the more accurate results of the data-driven RDMPC algorithm. From Figure 4, it is obvious that ζij and ζijr almost coincide. To see errors, although the simulation results fluctuate, the different values representing ζij−ζijr of each mobile robot are always about 10−13as shown in Figure 5, and such a small error indicates that the parameter differences among mobile robots converge to the expected value ζijr between them. Figure 6 displays the trajectories of five robots, where the triangle, circle, and diamond represent the initial position, position at 30 s, and the position at 50 s of each mobile robot, respectively. It can be seen that five mobile robots start from their respective starting points and keep their formation shape in the process of advancing according to each ideal trajectory.


### 4.2. Simulation 2

In this simulation, we choose the sine function as the trajectory. The desired trajectory is considered as follows:(38)Γ=xΓr=ζ,yΓr=sin0.2ζ.Letζ12r=−1,ζ23r=−1,ζ34r=−1,ζ45r=−1,andζ51r=4. The trajectory error between every two mobile robots can be expressed as xΓijζijr=−ζijr and yΓij=3. The communication topology is shown in Figure 1. The parameters and the bounded external disturbances of five mobile robots are the same as simulation 1. Taking the uncertainty set as δi=diag−0.5,−0.5,−0.5, the simulation is carried out for 50 s. The initial state values of each mobile robot are follows:(39)X10=111T,ζ10=1,X20=1−21T,ζ20=2,X30=−1.521T,ζ30=3,X40=1.5−21T,ζ40=4,X50=12−1.5T,ζ50=5.Comparing with the traditional DMPC, Figure7 reveals the state error Xi of each mobile robot. Figure 8 shows the absolute value of state error Xi of each mobile robot. The trajectory parameter difference ζij and the desired trajectory parameter difference ζijr of each mobile robot are shown in Figure 9. Figure 10 indicates the different values between ζij and ζijr of each mobile robot. Figure 11 displays the trajectories of five robots.Figure 7
The state errorXi of each mobile robot: (a) Robot 1; (b) Robot 2; (c) Robot 3; (d) Robot 4; (e) Robot 5.
(a)(b)(c)(d)(e)Figure 8
The absolute value of state errorXi of each mobile robot: (a) Robot 1; (b) Robot 2; (c) Robot 3; (d) Robot 4; (e) Robot 5.
(a)(b)(c)(d)(e)Figure 9
The trajectory parameter differenceζij and the desired trajectory parameter difference ζijr of each mobile robot: (a) Robot 1; (b) Robot 2; (c) Robot 3; (d) Robot 4; (e) Robot 5.
(a)(b)(c)(d)(e)Figure 10
The different values betweenζij and ζijr of each mobile robot: (a) Robot 1; (b) Robot 2; (c) Robot 3; (d) Robot 4; (e) Robot 5.
(a)(b)(c)(d)(e)Figure 11
Trajectories of five mobile robots.Through the simulation result, we can see that the fluctuation of the convergence curve using the data-driven RDMPC method is obviously been reduced than that using the traditional DMPC method in Figure7. Taking the absolute value of Xi, the result of data-driven RDMPC is more accurate than that of traditional DMPC as indicated in Figure 8. The simulation result of path parameter differences ζij of each mobile robot shown in Figure 9 is very close to ζijr after analyzing the different values between them of each mobile robot in Figure 10. And the difference values are always about 10−14 or 10−15, and we can think that ζij tracks ζijr well. Figure 11 displays the trajectories of five robots, where the triangle, circle, and diamond represent the initial position, position at 30 s, and the position at 50 s of each mobile robot, respectively; these five mobile robots keep their formation shapes ideally.


## 4.1. Simulation 1

In this simulation, the lane-change trajectory is considered. The desired trajectory is considered as follows:(35)Γ=xΓr=−ζ,yΓr=−21+e−2ζ−2.Letζ12r=−10,ζ23r=−10,ζ34r=15,ζ45r=−10,andζ51r=15. The trajectory error of every two mobile robots can be expressed as xΓijζijr=−ζijr and yΓij=2. The communication topology is shown in Figure 1. The parameters of robots are given as follows: Q1=100, S1=15, w1,l=4, Q2=100, S2=6, w2,l=4, Q3=120, S3=9.6, w3,l=4, Q4=100, S4=6, w4,l=4, Q5=120, S5=9.6, and w5,l=4. Set the sampling period T = 0.1 s and coefficients ΔCi=0,ΔCi=Bli, and Mi=Ni=−0.20.1. Each robot’s angle speed and line speed are ω=1 and v=1, respectively. Taking the uncertainty set as δi=diag−1,−1,−1, the simulation is carried out for 50 s. The initial state values of each mobile robot are as follows:(36)X10=0.25−0.050.1T,ζ10=10,X20=0.10.030.15T,ζ20=20,X30=0.15−0.020.1T,ζ30=30,X40=0.15−0.020.1T,ζ40=15,X50=0.15−0.020.1T,ζ50=25.Figure 1
Communication topology: (a) Robot 1; (b) Robot 2; (c) Robot 3; (d) Robot 4; (e) Robot 5.The bounded external disturbance of each mobile robot is as follows:(37)d1k=0.22k+1+21000cos0.1k0.2k+1+21000sin0.1k0.1k+1+11000cos0.2k,d2k=0.12k+1+11000cos0.1k0.2k+1+21000sin0.1k0.12k+1+11000cos0.2k,d3k=0.12k+1+11000cos0.1k0.15k+1+11000sin0.1k0.2k+1+21000cos0.2k,d4k=0.22k+1+11000cos0.1k0.15k+1+21000sin0.1k0.1k+1+11000cos0.2k,d5k=0.12k+1+11000cos0.1k0.2k+1+21000sin0.1k0.1k+1+21000cos0.2k.The simulation is compared with the traditional DMPC. Figure2 shows the state error Xi of each mobile robot. Figure 3 illustrates the absolute value of state error Xi of each mobile robot. Figure 4 displays ζij and ζijr of each mobile robot. Figure 5 presents the different values between ζij and ζijr of each mobile robot. Figure 6 demonstrates the trajectories of five mobile robots.Figure 2
The state errorXi of each mobile robot: (a) Robot 1; (b) Robot 2; (c) Robot 3; (d) Robot 4; (e) Robot 5.
(a)(b)(c)(d)(e)Figure 3
The absolute value of state errorXi of each mobile robot: (a) Robot 1; (b) Robot 2; (c) Robot 3; (d) Robot 4; (e) Robot 5.
(a)(b)(c)(d)(e)Figure 4
The trajectory parameter differenceζij and the desired trajectory parameter difference ζijr of each mobile robot: (a) Robot 1; (b) Robot 2; (c) Robot 3; (d) Robot 4; (e) Robot 5.
(a)(b)(c)(d)(e)Figure 5
The different values betweenζij and ζijr of each mobile robot: (a) Robot 1; (b) Robot 2; (c) Robot 3; (d) Robot 4; (e) Robot 5.
(a)(b)(c)(d)(e)Figure 6
Trajectories of five mobile robots.The result in Figure2 shows that compared with the traditional DMPC, the optimization process of data-driven RDMPC uses the optimal solution instead of the iterative algorithm, which is more efficient and greatly reduces the fluctuation of the convergence curve. In Figure 3, the simulation results with the control target Xi=000T indicate that the simulation results of data-driven RDMPC are always smaller than those of the traditional DMPC in the process of gradual stabilization, which highlights the more accurate results of the data-driven RDMPC algorithm. From Figure 4, it is obvious that ζij and ζijr almost coincide. To see errors, although the simulation results fluctuate, the different values representing ζij−ζijr of each mobile robot are always about 10−13as shown in Figure 5, and such a small error indicates that the parameter differences among mobile robots converge to the expected value ζijr between them. Figure 6 displays the trajectories of five robots, where the triangle, circle, and diamond represent the initial position, position at 30 s, and the position at 50 s of each mobile robot, respectively. It can be seen that five mobile robots start from their respective starting points and keep their formation shape in the process of advancing according to each ideal trajectory.


## 4.2. Simulation 2

In this simulation, we choose the sine function as the trajectory. The desired trajectory is considered as follows:(38)Γ=xΓr=ζ,yΓr=sin0.2ζ.Letζ12r=−1,ζ23r=−1,ζ34r=−1,ζ45r=−1,andζ51r=4. The trajectory error between every two mobile robots can be expressed as xΓijζijr=−ζijr and yΓij=3. The communication topology is shown in Figure 1. The parameters and the bounded external disturbances of five mobile robots are the same as simulation 1. Taking the uncertainty set as δi=diag−0.5,−0.5,−0.5, the simulation is carried out for 50 s. The initial state values of each mobile robot are follows:(39)X10=111T,ζ10=1,X20=1−21T,ζ20=2,X30=−1.521T,ζ30=3,X40=1.5−21T,ζ40=4,X50=12−1.5T,ζ50=5.Comparing with the traditional DMPC, Figure7 reveals the state error Xi of each mobile robot. Figure 8 shows the absolute value of state error Xi of each mobile robot. The trajectory parameter difference ζij and the desired trajectory parameter difference ζijr of each mobile robot are shown in Figure 9. Figure 10 indicates the different values between ζij and ζijr of each mobile robot. Figure 11 displays the trajectories of five robots.Figure 7
The state errorXi of each mobile robot: (a) Robot 1; (b) Robot 2; (c) Robot 3; (d) Robot 4; (e) Robot 5.
(a)(b)(c)(d)(e)Figure 8
The absolute value of state errorXi of each mobile robot: (a) Robot 1; (b) Robot 2; (c) Robot 3; (d) Robot 4; (e) Robot 5.
(a)(b)(c)(d)(e)Figure 9
The trajectory parameter differenceζij and the desired trajectory parameter difference ζijr of each mobile robot: (a) Robot 1; (b) Robot 2; (c) Robot 3; (d) Robot 4; (e) Robot 5.
(a)(b)(c)(d)(e)Figure 10
The different values betweenζij and ζijr of each mobile robot: (a) Robot 1; (b) Robot 2; (c) Robot 3; (d) Robot 4; (e) Robot 5.
(a)(b)(c)(d)(e)Figure 11
Trajectories of five mobile robots.Through the simulation result, we can see that the fluctuation of the convergence curve using the data-driven RDMPC method is obviously been reduced than that using the traditional DMPC method in Figure7. Taking the absolute value of Xi, the result of data-driven RDMPC is more accurate than that of traditional DMPC as indicated in Figure 8. The simulation result of path parameter differences ζij of each mobile robot shown in Figure 9 is very close to ζijr after analyzing the different values between them of each mobile robot in Figure 10. And the difference values are always about 10−14 or 10−15, and we can think that ζij tracks ζijr well. Figure 11 displays the trajectories of five robots, where the triangle, circle, and diamond represent the initial position, position at 30 s, and the position at 50 s of each mobile robot, respectively; these five mobile robots keep their formation shapes ideally.


## 5. Conclusions

The traditional distributed model predictive control has been widely used in the field of control; however, it is difficult to handle the model uncertainty of the system. This paper discusses a data-driven robust distributed model predictive control with the model uncertainty. The accuracy of the algorithm is improved by applying the robust performance objective minimization to replace the quadratic performance objective minimization when establishing the optimization problem with considering the model uncertainty of system. And on this basis, the fluctuation of convergence curve and the computational efficiency are abated due to the way of substituting for the iterative algorithm with optimal solution. The system conservatism is reduced after the data-driven strategy that the future predictive value is represented by employing the linear law of historical data which is put to use in deriving the control policy. The stability of the system is proved by Lyapunov theorem, and the effectiveness of this method is illustrated by a set of simulation examples.



---

*Source: 9923496-2021-06-21.xml*