Robust MPC for tracking of nonholonomic robots with
additive disturbances ⋆
Zhongqi Sun a, Li Dai a, Kun Liu a, Yuanqing Xia a, Karl Henrik Johansson b

arXiv:1703.03101v1 [cs.SY] 9 Mar 2017

a
b

School of Automation, Beijing Institute of Technology, Beijing 100081, China

ACCESS Linnaeus Centre and School of Electrical Engineering, KTH Royal Institute of Technology, SE-100 44 Stockholm,
Sweden

Abstract
In this paper, two robust model predictive control (MPC) schemes are proposed for tracking control of nonholonomic systems
with bounded disturbances: tube-MPC and nominal robust MPC (NRMPC). In tube-MPC, the control signal consists of a
control action and a nonlinear feedback law based on the deviation of the actual states from the states of a nominal system. It
renders the actual trajectory within a tube centered along the optimal trajectory of the nominal system. Recursive feasibility
and input-to-state stability are established and the constraints are ensured by tightening the input domain and the terminal
region. While in NRMPC, an optimal control sequence is obtained by solving an optimization problem based on the current
state, and the first portion of this sequence is applied to the real system in an open-loop manner during each sampling period.
The state of nominal system model is updated by the actual state at each step, which provides additional a feedback. By
introducing a robust state constraint and tightening the terminal region, recursive feasibility and input-to-state stability are
guaranteed. Simulation results demonstrate the effectiveness of both strategies proposed.
Key words: Robust control; Model predictive control (MPC); Nonholonomic systems; Bounded disturbances.

Tracking control of nonholonomic systems is a fundamental motion control problem and has broad applications in many important fields such as unmanned ground
vehicle navigation [1]; multi-vehicle cooperative control
[2]; formation control [3]; and so on. So far, many techniques has been developed for control of nonholonomic
robots [4,5,6,7,8]. However, these techniques either ignore the mechanical constraints, or require the persistent excitation of the reference trajectory, i.e., the linear and angular velocity must not converge to zero [9].
Model predictive control (MPC) is widely used in constrained systems. By solving a finite horizon open-loop
optimization problem on-line based on the current system state at each sampling instant, an optimal control

sequence is obtained. The first portion of the sequence
is applied to the system at each actuator update [10].
MPC for tracking of noholonomic systems was studied
in [2,9,11,12], where the robots were considered to be
perfectly modeled. However, when the system is uncertain or perturbed, then stability and feasibility of such
MPC may be lost. In the absence of constraints and
uncertainties, the optimal predictive control sequence
obtained by MPC is identical to that obtained by dynamic programming (DP), which provides an optimal
feedback policy or sequence of control laws [13]. Considering that feedback control is superior to open-loop
control in the aspect of robustness and that DP cannot
deal with the constrained systems, design methods for
MPC with robust guarantees is an urgent demand for
the tracking of constrained nonholonomic systems.

⋆ This paper was not presented at any IFAC meeting.
Corresponding author Y. Xia. Tel. +86 10 68914350. Fax
+86 10 68914382.
Email addresses: sunzhongqi12@gmail.com (Zhongqi
Sun), daili1887@gmail.com (Li Dai),
kunliubit@bit.edu.cn (Kun Liu),
xia yuanqing@bit.edu.cn (Yuanqing Xia), kallej@kth.se
(Karl Henrik Johansson).

There are several design methods for robust MPC. One
of the simplest approaches is to ignore the uncertainties and rely on the inherent robustness of deterministic
MPC [14,15], in which an open-loop control action solved
on-line is applied recursively to the system. However,
the open-loop control during each sampling period may
degrade the control performance even render the system unstable. Hence, feedback MPC was proposed in

1

Introduction

Preprint submitted to Automatica

10 March 2017

in this paper.

[16,17,18,19], in which a sequence of feedback control
laws is obtained by solving an optimization problem.
The determination of a feedback policy is usually prohibitively difficult. To overcome this difficulty, it is intuitive to focus on simplifying approximations by, for
instance, solving a min-max optimization problem online [17,18,19,20,21,22]. Min-max MPC provides a conservative robust solution for systems with bounded disturbances by considering all possible disturbances realizations. It is in most cases computationally intractable
to achieve such feedback laws, since the computational
complexity of min-max MPC grows exponentially with
the increase of the prediction horizon.

Robust MPC for linear systems is well studied but for
nonlinear systems is still challenging since it is usually
intractable to design a feedback law yielding a corresponding robust invariant set. Especially, the study of
robust MPC for nonholonomic systems remains open.
Motivated by the analysis above, this paper focuses on
the design of robust MPC for tracking of nonholonomic
systems with coupled input constraint and bounded additive disturbances. We discuss two robust MPC schemes
introduced above. First, a tube-MPC strategy with two
degrees of freedom is developed, in which the nominal
system is employed to generate a central trajectory and
a nonlinear feedback is designed to steer the system
trajectory of actual system within the tube for all admissible disturbances. Recursive feasibility and inputto-state stability are guaranteed by tightening the input
domain and terminal constraint via affine transformation and all the constraints are ensured. Since tube-MPC
sacrifices optimality for simplicity, an NRMPC strategy
is presented, in which the state of the nominal system is
updated by the actual one in each step. In such a way,
the control action applied to the real system is optimal
with respect to the current state. Input-to-state stability
is also established by utilizing the recursive feasibility
and the tightened terminal region.

Tube-MPC taking advantage both open-loop and feedback MPC was reported in [23,24,25,26,27,28,29]. Here
the controller consists of an optimal control action and a
feedback control law. The optimal control action steers
the state to the origin asymptotically, and the feedback
control law maintains the actual state within a “tube”
centered along the optimal state trajectory. Tube-MPC
for linear systems was advocated in [23,24,25], where the
center of the tube was provided by employing a nominal
system and the actual trajectory was restricted by an
affine feedback law. It was shown that the computational complexity is linear rather than exponential with
the increase of prediction horizon. The authors of [26]
took the initial state of the nominal system employed
in the optimization problem as a decision variable in
addition to the traditional control sequence, and proved
several potential advantages of such an approach. TubeMPC for nonlinear systems with additive disturbances
was studied in [27,28], where the controller possessed a
similar structure as in the linear case but the feedback
law was replaced by another MPC to attenuate the effect
of disturbances. Two optimization problems have to be
solved on-line, which increases the computation burden.

The remainder of this paper is organized as follows. In
Section 2, we outline the control problem and some preliminaries. Tube-MPC and NRMPC schemes are developed in Section 3 and Section 4, respectively, for tracking
of nonholonomic systems . In Section 5, Simulation results are given. Finally, we summarize the works of this
paper in Section 6.
Notation: R denotes the real space and N denotes the
collection of all positive integers. For a given matrix M ,
kM k denotes its 2-norm. diag{x1 , x2 , . . . , xn } denotes
the diagonal matrix with entries x1 , x2 , . . . , xn ∈ R.
For two vectors x = [x1 , x2 , . . . , xn ]T and y =
[y1 , y2 , . . . , yn ]T , x < y means {x1 < y1 , x2 <
y2 , . . . , xn < yn } and |x| ,√[|x1 |, |x2 |, . . . , |xn |]T denotes
its absolute value. kxk , xT x is the Euclidean
norm.
√
P -weighted norm is denoted as kxkP , xT P x, where
P is a positive define matrix with appropriate dimension.
Given two sets A and B, A ⊕ B , {a + b | a ∈ A, b ∈ B},
A ⊖ B , {a | {a} ⊕ B ⊂ A} and M A , {M a|a ∈ A},
where M is a matrix with appropriate dimensions.

In fact, tube-MPC provides a suboptimal solution because it has to tighten the input domain in the optimization problem, which may degrade the control capability.
It is natural to inquire if nominal MPC is sufficiently
robust to disturbances. A robust MPC via constraint restriction was developed in [24] for regulation of discretetime linear systems, in which asymptotic state regulation and feasibility of the optimization problem were
guaranteed. In [30], a robust MPC for discrete-time nonlinear system using nominal predictions was presented.
By tightening the state constraints and choosing a suitable terminal region, robust feasibility and input-statestability was guaranteed. In [31], the authors designed
a constraint tightened in a monotonic sequence in the
optimization problem such that the solution is feasible
for all admissible disturbances. A novel robust dualmode MPC scheme for a class of nonlinear systems was
proposed in [32], the system of which is assumed to be
linearizable. Since the procedure of this class of robust
MPC is almost the same as nominal MPC, we call this
class of robust MPC as nominal robust MPC (NRMPC)

2

Problem formulation and preliminaries

In this section, we first introduce the kinematics of the
nonholonomic robot and deduce the coupled input constraint from its mechanical model. Then, we formulate
the tracking problem as our control objective, and finally
give some preliminaries for facilitating the development
of our main results.

2

2.1

are presented as

Kinematics of the nonholonomic robot

v = (v L + v R )/2,
ω = (v R − v L )/2ρ.

Consider the nonholonomic robot described by the following unicycle-modeled dynamics:


cos θ(t) 0





˙ = f (ξ(t), u(t)) =  sin θ(t) 0  u(t),
ξ(t)


0
1

As a consequence, the control input u should satisfy the
constraint u ∈ U, where

(1)

U = {[v, ω]T :

where ξ(t) = [pT (t), θ(t)]T ∈ R2 × (−π, π] is the state,
consisting of position p(t) = [x(t), y(t)]T and orientation
θ(t), and u(t) = [v(t), ω(t)]T is the control input with
the linear velocity v(t) and the angular velocity ω(t).

w
r

r

2.2

Head position x h

q

y

x

The structure of the nonholonomic robot is shown in
Fig. 1. ρ is half of the wheelbase, v L and v R are the
velocities of the left and the right driving wheels of the
robot, respectively. Denote ξh by the head position which
is the point that lies a distance ρ along the perpendicular
bisector of the wheel axis ahead of the robot and is given
by








cos θ(t)
x(t)
xh (t)



 






ξh (t) =  yh (t)  =  y(t)  + ρ  sin θ(t) 
.
0
θ(t)
θh (t)

ξ̇f h (t) = fh (ξf h (t), uf (t)) + d(t),

cos θ(t) −ρ sin θ(t)

uf (t) ∈ U,

(7)

T
where ξf h (t) = [pT
is the state with the
f h (t), θf (t)]
head position pf h (t) = [xf h (t), yf h (t)]T , uf (t) =
T
[vf (t), ωf (t)]T is the control input, and d(t) = [dT
p (t), 0] ⊆
3
T
R , dp (t) = [dx (t), dy (t)] , is the external disturbances,
which is bounded by kd(t)k ≤ η.

(2)

The nominal system of the head position is then formulated as


(6)

The follower to be controlled is also an unicycle with
kinematics (1). Considering the existence of nonholonomic constraint, we consider its head position modeled
as (3). Furthermore, the robot is assumed to be perturbed by a disturbance caused by sideslip due to the
road ride. Therefore, we consider disturbances acting on
the linear velocity while neglecting disturbances acting
on the angular velocity. The perturbed head position
kinematics is then formulated as follows:

Fig. 1. The structure of the nonholonomic robot



Control objective

ξ˙r (t) = f (ξr (t), ur (t)).

vR



(5)

Our control objective is to track a reference trajectory
in a global frame O. The reference trajectory, which can
be viewed as a virtual leader, is described by a reference
T
2
state vector ξr (t) = [pT
r (t), θr (t)] ∈ R × (−π, π] with
pr (t) = [xr (t), yr (t)]T and a reference control signal
ur (t) = [vr (t), ωr (t)] ∈ U. The reference state vector
ξr (t) and the reference control signal ur (t) are modeled
as a nominal unicycle robot

x

O

|v| |ω|
+
≤ 1}
a
b

with b = a/ρ.

v
vL

(4)

Construct Frenet-Serret frames r O and f O for the virtual leader and the follower, respectively. They are moving frames fixed on the robots (see Fig. 2). The tracking
error prf = [xrf , yrf ]T with respect to the Frenet-Serret
frame f O is given by






ξ˙h (t) = fh (ξh (t), u(t)) = 
 sin θ(t) ρ cos θ(t)  u(t). (3)
0
1

prf (t) = R(θf (t))(pr (t) − pf (t)),
θrf (t) = θr (t) − θf (t),

It is assumed that the two wheels of the robot possess
the same mechanical properties and are bounded by
|v L | ≤ a and |v R | ≤ a, where a ∈ R is a known positive
constant. The linear and angular velocities of the robot

where R(θ) =

3

"

cos θ − sin θ
sin θ cos θ

#

(8)
(9)

is the rotation matrix.

y

y

where ũrf (t) is the input error and is given by

x
r

ũrf (t) =

O

y

xr

yrf

xrf

xrf + r

r

xf
f

"

−vf (t) + vr (t) cos θ̃rf (t)
−ρωf (t) + vr sin θ̃rf (t)

x

J(p̃rf (tk ), ũrf (tk )) =

tZ
k +T
tk

O

Taking the derivative of the tracking error yields

ṗrf (t) =

ωf (t)

−ωf (t)

0

#"

xrf (t)

#

yrf (t)
#
"
#
−vf (t) + vr (t) cos θrf (t)
dx (t)
+
+ R(θf )
. (10)
−ρωf (t) + vr (t) sin θrf (t)
dy (t)
"

2.3

Definition 1 For the nominal tracking error system
(12), the terminal region Ω and the terminal controller
ũκf (·) are such that if p̃rf (tk + T |tk ) ∈ Ω, then, for any
τ ∈ (tk + T, tk+1 + T ], by implementing the terminal
controller ũf (τ |tk+1 ) = ũκf (τ |tk+1 ), it holds that
p̃rf (τ |tk ) ∈ Ω,
ũf (τ |tk ) ∈ U,
ġ(p̃rf (τ |tk )) + L(p̃rf (τ |tk ), ũrf (τ |tk )) ≤ 0.

To distinguish the variables in the nominal system model
from the real system, we introduce ˜· as a superscript for
the variables in the nominal system. From the perturbed
system (7), the nominal dynamics can be obtained by
neglecting the disturbances as
ũf (t) ∈ U,

kprf (t)k ≤ β(kprf (t0 )k, t) + γ(η).

(11)

p̃˙ rf (t) =

0

ωf (t)

−ωf (t)

0

x̃rf (t)
ỹrf (t)

#

+ ũrf (t),

(18)

Definition 3 ([13]) A function V (·) is called an ISSLyapunov function for system (10) if there exist K∞
functions α1 (·), α2 (·), α3 (·) and a K function σ(·) such
that for all prf ∈ R2

is the state
where, similarly, ξ̃f h (t) =
of the nominal system with the position p̃f h =
[x̃f h (t), ỹf h (t)]T and orientation θ̃f (t), and ũf (t) =
[ṽf (t), ω̃f (t)]T is the control input of the nominal system.
The tracking error dynamics based on the nominal
system is then given by
#"

(15)
(16)
(17)

Definition 2 ([33]) System (10) is input-to-state stable
(ISS) if there exist a KL function β(·, ·) : R≥0 ×R≥0 → R
and a K function γ(·) such that, for t ≥ 0, it holds that

T
[p̃T
f h (t), θ̃f (t)]

"

Preliminaries

Some definitions and lemmas used in the following sections are summarized as follows.

Based on the discussion above, we will design robust
MPC strategies to drive the tracking error prf to a neighborhood of the origin. Note that the tracking system (10)
involves the disturbances but the future disturbances
cannot be predicted in advance. We will formulate the
MPC problem only involving the nominal system.

ξ̃˙ f h (t) = fh (ξ˜f h (t), ũf (t)),

(14)

in which L(p̃rf (τ |tk ), ũrf (τ |tk )) = kp̃rf (τ |tk )k2Q +
kũrf (τ |tk )k2P represents the stage cost with the positive
define matrices P = diag{p1 , p2 } and Q = diag{q1 , q2 },
g(p̃rf (τ |tk )) = 12 kp̃rf (τ |tk )k2 is the terminal penalty,
and T is the prediction horizon satisfying T = N δ,
N ∈ N.

Fig. 2. Leader-follower configuration

0

(13)

L(p̃rf (τ |tk ), ũrf (τ |tk ))dτ

+g(p̃rf (tk + T |tk )),

x

"

.

Define {tk : k ∈ N, tk+1 − tk = δ}, with δ > 0, the time
sequence at which the open-loop optimization problems
are solved. The MPC cost to be minimized is given by

x fh

O

#

α1 (kprf (tk )k) ≤ V (prf (tk )) ≤ α2 (kprf (tk )k),
(19)
V (prf (tk+1 ))−V (prf (tk )) ≤ −α3 (kprf (tk )k)+σ(η). (20)
Remark 1 It should be mentioned that both Definition 2
and Definition 3 result in input-to-state stability, which
implies that the tracking error vanishes if there is no
disturbance.

(12)

4

The following lemma provides a terminal controller and
the corresponding terminal region for the nominal error
system (12).

Lemma 2 System (3) with u ∈ U is locally Lipschitz in
ξh with Lipschitz constant a, where a is the max wheel
speed.

Lemma 1 For the nominal tracking system (12), let
ũf ∈ λf U with λf ∈ (0, 1], |vr | < √a2 λf , and λr =

Proof. Considering the function values of fh (ξh , u) at
ξh1 and ξh2 with the same u, we have

√
2
a max |vr |. Then
λr )} is a terminal

ũκf (τ |tk ) =

"

Ω = {p̃rf : k̃1 |x̃rf | + k̃2 |ỹrf | < a(λf −
region for the controller

ṽfκ (τ |tk )

ω̃fκ (τ |tk )

#

=

"

k̃1 x̃rf + vr cos θ̃rf
1
ρ (k̃2 ỹrf

+ vr sin θ̃rf )

kfh (ξh1 , u) − fh (ξh2 , u)k2


v(cos θ1 − cos θ2 ) + ρω(sin θ2 − sin θ1 )



= 
 v(sin θ1 − sin θ2 ) + ρω(cos θ1 − cos θ2 ) 
0

#

τ ∈ (tk + T, tk+1 + T ], (21)

with the parameters satisfying
pi qi <


√
√
1− 1−4pi qi 1+ 1−4pi qi
,
, i = 1, 2.
2pi
2pi

1
4

= v 2 (cos θ1 − cos θ2 )2 + ρ2 ω 2 (sin θ2 − sin θ1 )2
+v 2 (sin θ1 − sin θ2 )2 + ρ2 ω 2 (cos θ2 − cos θ1 )2
≤ 2(v 2 + ρ2 ω 2 )(θ1 − θ2 )2
≤ 2 max {v 2 + ρ2 ω 2 }(θ1 − θ2 )2

and k̃i ∈

[v,ω]T ∈U
2
= a (θ1 − θ2 )2 ,

Proof. First, consider the terminal controller
|ṽfκ | |ω̃fκ | |k̃1 x̃rf +vr cos θ̃rf )| |k̃2 ỹrf + vr sin θ̃rf |
+
=
+
a
b
a
a
1
≤ (k̃1 |x̃rf | + k̃2 |ỹrf | + |vr cos θ̃rf | + |vr sin θ̃rf |)
a
√
2
≤ λf − λr +
|vr | ≤ λf ,
a

where the mean value theorem and Lagrange multiplier
method are used in the last inequality. The maximum of
v 2 + ρ2 ω 2 , subject to |v|/a + |ω|/b ≤ 1, can be obtained
by setting v = a2 and ω = 2b . From the results above, we
conclude that
kfh (ξh1 , u) − fh (ξh2 , u)k ≤ akξh1 − ξh2 k.

which implies uκf ∈ λf U if p̃rf ∈ Ω.

(23)
✷

Next, choose g(p̃rf (τ |tk )) as Lyapunov function. The
derivative of g(p̃rf (τ |tk )) with respect to τ yields

3

2
(τ |tk )),
ġ(p̃rf (τ |tk )) = −(k̃1 x̃2rf (τ |tk ) + k̃2 ỹrf

Tube-MPC

In this section, a tube-MPC policy is developed, which
consists of an optimal control action obtained by solving
an optimization problem and a feedback law based on
the deviation of the actual state from the nominal one.
The controller forces the system state to stay within
a tube around a sensible central trajectory. The cental
trajectory is determined by the following optimization
problem.

which means that Ω is invariant by implementing the
terminal controller, i.e., p̃rf (τ |tk ) ∈ Ω holds for all τ > tk
once p̃rf (tk |tk ) ∈ Ω.
Finally, for p̃rf (τ |tk ) ∈ Ω, it follows that
ġ(p̃rf (τ |tk )) + L(p̃rf (τ |tk ), ũrf (τ |tk ))
2
2
2
+ p1 ṽrf
+ p2 ω̃rf
= x̃ef x̃˙ rf + ỹrf ỹ˙ rf + q1 x̃2rf + q2 ỹrf

Problem 1

2
2
2
2
= −k̃1 x̃2rf − k̃2 ỹrf
+ q1 x̃2rf + q2 ỹrf
+ p1 ṽrf
+ p2 ω̃rf
2
.
= (p1 k̃12 − k̃1 + q1 )x̃2rf + (p2 k̃22 − k̃2 + q2 )ỹrf

2

(22)

min

ũf (τ |tk )


√
√
i qi 1+ 1−4pi qi
,
,
Since pi qi < 14 and k̃i ∈ 1− 1−4p
2pi
2pi
i = 1, 2, the inequality ġ + L < 0 holds.

J(p̃rf (tk ), ũrf (tk )),

(24)

s.t. ξ̃f h (tk |tk ) = ξ˜f h (tk ),
ξ̃˙ f h (τ |tk ) = fh (ξ˜f h (τ |tk ), ũf (τ |tk )),
ũf (τ |tk ) ∈ Utube ,
p̃rf (tk + T |tk ) ∈ Ωtube ,



Hence, from Definition 1, Ω is a terminal region associated with the terminal controller ũκf (τ |tk ).
✷

where Utube = {[ṽf , ω̃f ]T :

The nominal system (3) is Lipschitz continuous and a
corresponding Lipschitz constant is given by the following lemma.

√
√
λtube = 22 − η a 2 ,
a(λtube − λr )}.

5

|ṽf |
a

+

|ω̃f |
b

(25)
(26)
(27)
(28)

≤ λtube } with

and Ωtube = {p̃rf : k̃1 |x̃rf |+ k̃2 |ỹrf | <

Solution of Problem 1 yields the minimizing control sequence for the nominal follower system over the interval
[tk , tk+T ]:

the feedback law are combined in a different manner compared to linear systems [23,24,25]. This increases the
difficulty of determining the tightened input constraint
set Utube such that uf ∈ U holds. The scheme is also
different from the existing works on nonlinear systems
as in [27] and [28], in which our feedback law determined
off-line is replaced with an online computation of another
MPC. Hence, two optimization problems have to be solved
in each step, which increases the computational burden.

ũ∗f (tk ) = {ũ∗f (tk |tk ), ũ∗f (tk+1 |tk ), . . . , ũ∗f (tk+N |tk )}, (29)
as well as the corresponding optimal trajectory:
ξ̃f∗h (tk ) = {[p̃∗f h (tk |tk ), θ̃f∗ (tk |tk )]T , [p̃∗f h (tk+1 |tk ),

θ̃f∗ (tk+1 |tk )]T , . . . , [p̃∗f h (tk+N |tk ), θ̃f∗ (tk+N |tk )]T }. (30)

Remark 4 From Algorithm 1, it can be observed that the
optimization problem employs only the nominal system
and thus the predictive optimal trajectory is independent
of the actual state except for the initial one. From this
point, the central trajectory of the tube can be calculated
in a parallel or even off-line way if the initial state is
known a priori. In such a way, only one feedback law is
required to be calculated on-line, which reduces the online computational burden even further.

The robust controller for the follower over the interval
[tk , tk+1 ] is designed as
uf (tk ) = M −1 (θf (tk ))M (θ̃f∗ (tk |tk ))ũ∗f (tk |tk )

+M −1 (θf (tk ))K(pf h (tk ) − p̃∗f h (tk |tk )),

where M (θ) =

"

cos θ −ρ sin θ

(31)

#

Before stating the main results of tube-MPC, the following lemma is given to show that the feedback law
renders the difference between the minimizing trajectory
and the actual trajectory bounded while guaranteeing
the satisfaction of the input constraint.

, K = diag{kx , ky },

sin θ ρ cos θ
kx < 0, ky < 0, is the feedback gain, ũ∗f (tk |tk ) is the
first control action of the optimal control sequence, and
p̃∗f h (tk |tk )) and θ̃f∗ (tk |tk ) are the first portion of the optimal position and orientation, respectively.

Lemma 3 For the tracking control system (10) with controller (31), it follows that

Based on this control strategy, the procedure of tubeMPC is summarized in Algorithm 1.

(i) the state of the real system lies in the tube T =
"
#
− kη1
∗
pf h ⊕ Pf e , where Pf e = {pf e (t) : |pf e | ≤
};
− kη2
(ii) the input constraint is satisfied, i.e., uf ∈ U.

Algorithm 1 Tube-MPC
1: At time t0 , initialize the nominal system state by the
actual state ξ̃f h (0) = ξf h (0).
2: At time tk , solve Problem 1 based on nominal system
to obtain the optimal control sequence ũ∗f (tk ) =
arg minuf (τ |tk ) Jf (tk , p̃rf , ũrf ).
3: Calculate the actual control signal for the real
system uf (tk ) = M −1 (θf )[M (θ̃f∗ )ũ∗f (tk |tk ) +
K(pf h (tk ) − p∗f h (tk )].
4: Apply the first portion of the sequence, i.e.,
ũ∗f (tk |tk ), to the nominal system, and apply uf (tk )
to the real system during the sampling interval
[tk , tk+1 ].
5: Update the state of the nominal system with
ξ̃f∗h (tk+1 ) and the state of the real system with
ξf h (tk+1 ).
6: Update the time instant tk = tk+1 and go to step 2.

Proof. Denote the deviation of the actual trajectory
from the optimal trajectory as
pf e (t) = pf h (t) − p̃∗f h (t).

(32)

Taking the derivative of (32) yields
ṗf e (t) =

"

cos θf (t) −ρ sin θf (t)

#

uf (t)
sin θf (t) ρ cos θf (t)
#
"
cos θ̃f∗ (t) −ρ sin θ̃f∗ (t)
ũ∗f (t) + dp (t). (33)
−
sin θ̃f∗ (t) ρ cos θ̃f∗ (t)

Substituting (31) into (33), we can conclude that

Remark 2 Since the optimization problem is solved online at each step and the first optimal control action is
employed to generate the control policy together with the
feedback law, the computational complexity is determined
by the nominal system. Hence, the scheme has the same
computational complexity as the deterministic MPC.

ṗf e (t) = Kpf e (t) + dp (t),

(34)

of which the solution is given by
pf e (t) = e

Remark 3 Due to the nonlinearity and nonholonomic
constraint of the system, the optimal control action and

Kt

pf e (0) +

Zt
0

6

eK(t−τ ) dp (τ )dτ.

(35)

By the initialization stage (25) and the upper-bound of
the disturbances, it follows that
|pf e (t)| ≤ η

"

1 k1 t
k1 e
1 k2 t
k2 e

1
k1
1
k2

−
−

#

.

KPf e =

Pf e (t) =

pf e (t) : |pf e (t)| ≤ η

(36)

t→∞

"

1 k1 t
k1 e
1 k2 t
k2 e

1
k1
1
k2

#)

− kη1

#)

−
−

[

θf∗ ∈(−π,π]

.

(37)

(

pf e (t) : |pf e | ≤

"

− kη2

.

pf h ∈ p∗f h ⊕ Pf e ,

(38)

For (ii), redefine the control input as

(i) Problem 1 is feasible for all t > t0 ;
(ii) the tracking control system (10) is ISS.

It can be observed that M (·) is an affine transformation,
which is equivalent to scaling ωf (ωf∗ ) by ρ and rotate
uf (u∗f ) by θ. Thus, to prove uf ∈ U if u∗f ∈ Utube is
a
equivalent to show uaf ∈ Ua if ua∗
f ∈ Utube for every
admissible θf and θf∗ . The sets Ua and Uatube are defined
as follows:
Ua = M (θf ){[v, ω]T :

(42)

Uatube

(43)

Proof. From Lemma 1, Ωtube is a terminal region by
letting λf = λtube . We assume that a feasible solution
exists and an optimal solution ũ∗f (tk ) is found at the
sampling instant tk . When applying this sequence to
the nominal system, the tracking error of the nominal
system is driven into the terminal region Ωtube , i.e.,
p̃e∗
f (tk + T |tk ) ∈ Ωtube , along the corresponding openloop trajectory ξ̃f∗ (tk ) over [tk , tk + T ]. In terms of Algorithm 1, the open-loop control ũ∗ (tk |tk ) is applied to the
nominal system, and its state measurement at time tk+1
is given by ξ̃(tk+1 ) = ξ̃ ∗ (tk+1 |tk ). Therefore, to solve
the open-loop optimal control problem at tk+1 with the
initial condition, a feasible solution can be constructed
by

Substituting (40) into (31) yields
(44)

ũf (τ |tk+1 ) =

It is obvious that

θf ∈(−π,π]

[

θf∗ ∈(−π,π]

√
2
},
U = {[v, ω] : k[v, ω] k ≤ a
2

(45)

Uatube = {[v, ω]T : k[v, ω]T k ≤ aλtube },

(46)

a

(48)

θf ∈(−π,π]

Theorem 1 For the tracking control system (10) under
Algorithm 1, if Problem 1 is feasible at time t0 , then,

(40)
(41)

uaf = ua∗
f + Kpf e .

Ua ,

The main results of tube-MPC are given in the following
theorem.

i.e., the trajectory lies in the tube T.

|v| |ω|
+
≤ 1},
a
a
|v| |ω|
= M (θf∗ ){[v, ω]T :
+
≤ λtube }.
a
a

\

Remark 5 Note that the input domain is independent
of the feedback gain K, which differs from the results
of linear systems in [23,24,25]. Meanwhile, from (i) in
Lemma 3, increasing K will reduce the difference between
the actual trajectory and the optimal one, and consequently reduce the size of the tube T. It indicates that the
steady tracking performance could be enhanced by tuning
K.

(39)

uaf = M (θf )uf ,
∗ ∗
ua∗
f = M (θf )uf .

Uatube ⊕ KPf e ⊂

which implies that uaf ∈ Ua holds for every admissible
θf and θf∗ , and uf ∈ U naturally holds.
✷

From (32) and pf e ∈ Pf e , we have

\

(47)

Thus, it can be obtained that

We further define Pf e as
Pf e = lim Pf e (t) =

Kpf e (t) : |Kpf e | ≤

" #)
η

η
n
√ o
⊂ Kpf e (t) : kKpf e k ≤ 2η .

Consequently, pf e (t) ∈ Pf e (t), where the set Pf e (t) is
defined by
(

(

T

T

(

ũ∗f (τ |tk ), τ ∈ [tk+1 , tk + T ),

ũκf (τ |tk ), τ ∈ [tk + T, tk+1 + T ),

(49)

where ũκf (τ |tk ) is the terminal controller given by (21).
Since the terminal region Ωtube is invariant with the
e
control ũκf (τ |tk ), p̃e∗
f (tk +T |tk ) ∈ Ωtube implies p̃f (tk+1 +
T |tk+1 ) ∈ Ωtube . Then, result (i) can be achieved by
induction.

7

It follows from pf r (t) = R(θf )(p̃∗f r (t) + pf e (t)) and
pf e (0) = 0 that

For (ii), we first prove that the tracking error for the
nominal system converges to the origin. Then we show
that the state of the real system converges to an invariant
set along a trajectory lying in the tube T, the center
of which is the trajectory of the nominal system. The
Lyapunov function for the nominal system is chosen as
V (tk ) = J(p̃∗rf (tk ), ũ∗rf (tk )).

kpf r (t)k ≤ β(kpf r (t0 )k, t) + γ(η).

Therefore, the solution of system (10) is asymptotically
ultimately bounded with Algorithm 1 and the closedloop system is ISS.
✷

(50)

Consider the difference of the Lypunov function at tk
and tk+1 ,

4

tk+1
Z +T

(kp̃rf (τ |tk+1 )k2Q + kũrf (τ |tk+1 )k2P )dτ

tk+1

−

tZ
k +T

(kp̃∗rf (τ |tk )k2Q + kũ∗rf (τ |tk )k2P )dτ

tk

+kp̃rf (tk+1 +
=−
+

tZk+1

tk

T |tk+1 )k2R

−

kp̃∗rf (tk

+

Problem 2

T |tk )k2R

min

ũf (τ |tk )

(kp̃∗rf (τ |tk )k2Q + kũ∗rf (τ |tk )k2P )dτ
(kp̃rf (τ |tk+1 )k2Q + kũrf (τ |tk+1 )k2P )dτ

tk +T

+kp̃rf (tk+1 + T |tk+1 )k2R − kp̃∗rf (tk + T |tk )k2R .
(51)

(52)

uf (tk ) = ũ∗f (tk |tk ).

(53)

Furthermore, because of pf e ∈ Pf e for all t > t0 , there
exists a K function γ(·) such that
∀t > t0 .

(61)

(62)

Algorithm 2 NRMPC
1: At time tk , initialize the nominal system state by the
actual one ξ˜f h (tk ) = ξf h (tk ).
2: Solve Problem 2 based on the nominal system to
obtain the minimizing control sequence ũ∗f (tk ) =
arg minũf (τ |tk ) J(p̃rf (tk ), ũrf (tk )).
3: Apply the first portion of the sequence to the real
system, i.e., uf (tk ) = ũ∗f (tk |tk )), during the interval
[tk , tk+1 ].
4: Update the state of the real system with ξf h (tk+1 ).
5: Update the time instant tk = tk+1 and go to step 1.

Due to the asymptotic stability of the nominal system,
there exists a KL function β(·, t), such that

kpf e (t)k ≤ γ(η),

(60)

The NRMPC strategy is then described in Algorithm 2.

Substituting (52) into (51), we have ∆V ≤ 0, which
implies that the tracking error for the nominal system
converges to the origin asymptotically.

∀t > t0 .

(58)
(59)

Problem 2 yields a minimizing control sequence over the
interval [tk , tk + T ] of the same form as in (29) as well
as a minimizing trajectory as in (30). The control input
over [tk , tk+1 ] is chosen as

tk +T

kp̃∗rf (t)k ≤ β(kp̃∗rf (t0 )k, t),

(57)

k̃1 +k̃2

2
kp̃rf (tk+1 + T |tk+1 )k2R − kpe∗
f (tk + T |tk )kR

L(p̃ef (τ ), ũrf (τ ))dτ ≤ 0.

(56)

a(1−λr )
, Ωε = {p̃rf : kp̃rf k ≤ ε}, and ε < r.
where r = √
2
2

By integrating (17) form tk + T to tk+1 + T , it follows
that

+

J(p̃rf (tk ), ũrf (tk )),

s.t. ξ̃f h (tk |tk ) = ξf h (tk ),
˙
ξ̃ f h (τ |tk ) = fh (ξ˜f h (τ |tk ), ũf (τ |tk )),
ũf (τ |tk ) ∈ U,
rT
kp̃rf (τ |tk )k ≤
,
τ − tk
p̃rf (tk + T |tk ) ∈ Ωε ,

tk+1
Z +T

tk+1
Z +T

NRMPC

In this section, an NRMPC strategy is developed. The
state of the nominal system is updated by the actual
state at each sampling instant. Unlike tube-MPC, the
control sequence obtained is optimal with respect to the
current actual state, and only the first control action
of the sequence is applied to the real system. The optimization problem of the NRMPC strategy is defined as
follows:

∆V = V (tk+1 ) − V (tk )
≤ J(p̃rf (tk+1 ), ũrf (tk+1 )) − J(p̃∗rf (tk ), ũ∗rf (tk ))
=

(55)

(54)

8

Remark 6 For Problem 2, the state of nominal system
is updated by the actual one at each step. As a result, the
optimization problem has to be solved on-line. However,
such an updating strategy yields an optimal control with
respect to the current state. The scheme has the same
computational burden as the deterministic MPC.

+

tZk+1

tk

kd(τ )kdτ

≤ ηδ + a

tZk+1

kξf h (τ ) − ξ˜f∗h (τ |tk )kdτ

tk

Remark 7 Note that the input domain of NRMPC is
larger than that of tube-MPC. NRMPC may therefore
have better tracking capability.

≤ ηδeaδ .

We construct a feasible control sequence at tk+1 for the
nominal system as follows.

Remark 8 As shown in Algorithm 2, an open-loop control action is applied to the real system during each sampling interval. However, the existence of disturbances
may lead to an error between the actual trajectory and
the optimal prediction. This increases the difficulty of
analyzing the recursive feasibility using the conventional
methods for MPC.

ũf (τ |tk+1 ) =

r
k̃δ ≥ ln ,
ε

−ξ˜f∗h (tk+1 |tk ) −
≤ ηδeaδ + a

tk

−ξ˜f∗h (tk |tk ) −

fh (ξ˜f∗h (s|tk ), ũ∗f (s|tk ))dsk.

tk+1

kξ˜f h (s|tk+1 ) − ξ̃f∗h (s|tk )kds.

(66)

(67)

Substituting tk + T into (67) leads to
kξ˜f h (tk + T |tk+1 ) − ξ̃f∗h (tk + T |tk )k ≤ ηδeaT .

(68)

Due to the fact kp̃rf (tk + T |tk+1 ) − p̃∗rf (tk + T |tk )k ≤
kξ˜f h (tk + T |tk+1 ) − ξ̃f∗h (tk + T |tk )k and the application
of triangle inequality, we arrive at
kp̃rf (tk + T |tk+1 )k ≤ kp̃∗rf (tk + T |tk )k + ηδeaT .


fh (ξf h (τ ), ũ∗f (·)) + d(τ ) dτ

(69)

Since p̃∗rf (tk + T |tk ) ∈ Ωε , i.e. kp̃∗rf (tk + T |tk )k ≤ ε, and
η≤

tZk+1

fh (ξ˜f∗h (τ |tk ), ũ∗f (·))dτ k

e−aT
δ

(r − ε), we obtain

kp̃rf (tk + T |tk+1 )k ≤ r,

(70)

which implies p̃rf (tk + T |tk+1 ) ∈ Ω.

tZk+1

tk

Zτ

kξ˜f h (τ |tk+1 ) − ξ˜f∗h (τ |tk )k ≤ ηδea(τ −tk+1 +δ) .

tk

≤

(65)

Applying Grönwall-Bellman inequality yields
(63)

kξf h (tk+1 ) − ξ˜f∗h (tk+1 |tk )k


Zτ

tk+1

Proof. Since Problem 2 is feasible at tk , applying the
first control of the sequence ũ∗f during the interval
(tk , tk+1 ] to the real system may lead to a difference
of the trajectory between the actual system and the
nominal system. At tk+1 , this difference is bounded by

= kξf h (tk ) +

ũκf (τ |tk ), τ ∈ (tk + T, tk+1 + T ].

tk+1

where k̃ = min{k̃1 , k̃2 }.

tZk+1

ũ∗f (τ |tk ), τ ∈ (tk+1 , tk + T ],

kξ˜f h (τ |tk+1 ) − ξ̃f∗h (τ |tk )k
Zτ
= kξf h (tk+1 ) +
fh (ξ˜f h (s|tk+1 ), ũ∗f (s|tk ))ds

Lemma 4 For the tracking control system (10), assume
that there exists an optimal control sequence ũ∗f (tk ) at
instant tk such that p̃rf (tk + T |tk ) ∈ Ωε , and apply the
first control of the sequence to the perturbed system (7).
Then, there exists a control sequence ũf (tk+1 ) at tk+1
such that p̃rf (tk+1 + T |tk+1 ) ∈ Ωε , if
e−aT
(r − ε) ,
δ

(

First, we consider the interval τ ∈ (tk+1 , tk + T ].
Since the state of the nominal system is updated by
ξ̃(tk+1 |tk+1 ) = ξ(tk+1 ), we have

The following two lemmas guarantee recursive feasibility
of Problem 2. Lemma 4 states the existence of the control
sequence that is able to drive the state of the nominal
system into Ωε in prediction horizon T , and Lemma 5
shows that the state constraint is satisfied by employing
that control sequence to the nominal system.

η≤

(64)

h
i
kfh (ξf h (τ ), ũ∗f (·)) − fh (ξ˜f∗h (τ |tk ), ũ∗f (·))k dτ

Next, consider the interval τ ∈ (tk + T, tk+1 + T ], during which the local controller ũκf (τ |tk ) is applied to the
9

kp̃∗rf (tk + T |tk )k ≤ ε. Since

nominal system

rT
τ −tk+1

Theorem 2 For the tracking control system (10), suppose that Problem 1 is feasible at time t0 and the parameters satisfy the conditions in Lemma 4 and Lemma 5.
Then,

Applying the comparison principle yields
kp̃rf (tk+1 + T |tk+1 )k ≤ kp̃rf (tk + T |tk+1 )ke−δk̃ .

(i) Problem 1 is feasible for all t > t0 ;
(ii) the tracking control system (10) is ISS if

r
that
ε

kp̃rf (tk+1 + T |tk+1 )k ≤ ε.

qε2 >

(71)

This proves the existence of a control sequence at tk+1
which is able to drive the tracking error of the nominal
system into the terminal region Ωε .
✷

r(T − δ)
.
T

(72)

rT
Proof. To prove kp̃f r (τ |tk+1 )k ≤ τ −t
, τ ∈ (tk+1 , tk+1 +
k+1
T ], we first consider the interval τ ∈ (tk+1 , tk + T ]. From
(67), it follows that

kp̃f r (τ |tk+1 )k ≤ kp̃∗rf (τ |tk )k + ηδeaT .
Due to (63) and p̃∗f r (τ |tk ) ≤
kp̃f r (τ |tk+1 )k ≤

T
τ −tk ,

(73)

δrT
δr
≤
.
T −δ
(τ − tk+1 )(τ − tk )

(74)

(78)

According to Riemann integral principle, there exists a
constants 0 < c1 ≤ δ such that
V (prf ) ≥ c1 L(tk , , urf ) , α1 (kprf k),

(79)

where α1 (·) is obviously a K∞ function.

(75)

On the other hand, from (17), we have

Substituting (75) into (74), we obtain
rT
,
kp̃f r (τ |tk+1 )k ≤
τ − tk+1

Proof. (i) Assume Problem 2 is feasible at instant tk ,
then feasibility of Problem 2 at tk+1 implies the existence
of a control sequence that is able to drive the state
of the nominal system to the terminal region Ωε while
satisfying all the constraints. In terms of Algorithm 2,
the first control of the optimal control sequence is applied to the system. From Lemma 4, at tk+1 , a feasible
control sequence in (65) renders p̃rf (tk+1 +T |tk+1 ) ∈ Ωε
while satisfying ũf (τ |tk+1 ) ∈ Ω for τ ∈ [tk+1 , tk+1 + T ].
Meanwhile, From Lemma 5, by applying the control
sequence (65) to the nominal system (11), the trajectory
rT
constraint is satisfied, i.e., kp̃rf (τ |tk+1 )k ≤ τ −t
, imk+1
plying the feasibility of Problem 2 at tk+1 . Hence, the
feasibility of Problem 2 at the initial time t0 results in
the feasibility for all t > t0 by induction.

V (prf (tk )) = J(p̃∗rf (tk ), ũ∗rf (tk )).

From (72), we have
r−ε≤

(77)

(ii) Choose a Lyapunov function as follows

we obtain

rT
+ (r − ε).
τ − tk

1 aT
q 2 η 2 δ 2aT
ηe (r + ε) +
(e
− e2aδ )
2
2a
1
2q 2 ηr T 2
1
+ √ (
− T ) 2 (e2aT − e2aδ ) 2 ,
2a δ

where q = min{q1 , q2 }.

Lemma 5 For the tracking control system (10), assume
that there exists an optimal control sequence ũ∗f (tk ) at
instant tk such that the trajectory constraint is satisfied,
i.e., p̃∗f r (τ |tk ) ≤ τ rT
−tk , and apply the first control of the
sequence to the perturbed system (7). Then, at tk+1 , by
the control sequence (65), the trajectory constraint is also
satisfied, if the parameter ε satisfies
ε≥

≥ r, kp̃f r (τ |tk+1 )k ≤

is naturally satisfied over the interval τ ∈ [tk +
T, tk+1 + T ], thereby completing the proof.
✷

d
2
kp̃rf (τ |tk+1 )k2 = −2(k̃1 x̃2rf + k̃2 ỹrf
)
dτ
≤ −2k̃kp̃rf (τ |tk+1 )k2 .

It follows from k̃δ ≥ ln

rT
τ −tk+1

V (prf (tk )) ≤ g(prf (tk )) + g(prf (tk + T |tk )), ∀prf ∈ Ωε .

(76)

Due to the decreasing property of g(·) in Ωε with respect
to time, it follows that

which proves that the state constraint is satisfied over
the interval τ ∈ [tk+1 , tk + T ].

V (prf (tk )) ≤ 2g(prf (tk )),

Next, consider the interval τ ∈ [tk + T, tk+1 + T ]. By
Lemma 4, it holds that kp̃rf (tk + T |tk+1 )k ≤ r once

Because the origin lies in the interior of Ωε and
2g(prf (t)) ≤ ε, ∀prf ∈ Ωε , it holds that 2g(prf (t)) ≥ ε
10

∀prf ∈ Ωε .

(80)

if prf ∈ R2×2 \ Ωε . Due to the feasibility of Problem
2, there exists an upper-bound c2 > ε for V (prf (t)).
c2
Thus α2 (kprf k) =
g(prf (t)) is a K∞ function such
ε
that α2 (kprf k) ≥ c2 thereby satisfying α2 (kprf (t)k) ≥
V (prf (t)).




∆V1 ≤ 

, ∆V1 + ∆V2 + ∆V3 ,

(82)

(kp̃rf (τ |tk+1 )k2Q

tk +T

+

+kp̃rf (tk+1 + T |tk+1 )k2R − kp̃∗rf (tk + T |tk )k2R

kũrf (τ |tk+1 )k2P )dτ

+kp̃rf (tk + T |tk+1 )k2R − kp̃rf (tk + T |tk+1 )k2R . (86)

+kp̃rf (tk+1 + T |tk+1 )k2R − kp̃∗rf (tk + T |tk )k2R ,

Integrating (17) from tk +T to tk+1 +T and substituting
it into (86) leads to

tZk+1

(kp̃∗rf (τ |tk )k2Q + kũ∗rf (τ |tk )k2P )dτ.

tk

∆V2 ≤ kp̃rf (tk + T |tk+1 )k2R − kp̃∗rf (tk + T |tk )k2R
1
≤ ( kp̃rf (tk + T |tk+1 ) − p̃∗rf (tk + T |tk )k)
2

× kp̃rf (tk + T |tk+1 )k + kp̃∗rf (tk + T |tk )k
1
(87)
≤ ησeaT (ε + r).
2

For ∆V1 , it holds that
tZ
k +T

tk+1

(kp̃rf (τ |tk+1 ) − p̃∗rf (τ |tk )kQ )

×(kp̃rf (τ |tk+1 )kQ + kp̃∗rf (τ |tk )kQ )dτ.

(83)

For ∆V3 , we first assume kp̃rf (tk+1 |tk )k > ε, which
implies kp̃rf (τ |tk )k > ε, τ ∈ (tk , tk+1 ], and thus we
obtain

By (67), we have

∆V1 ≤

tZ
k +T

tk+1

(kp̃rf (τ |tk+1 )k2Q + kũrf (τ |tk+1 )k2P d)τ

∆V2 =

tk +T

∆V1 ≤

(85)

tk+1
Z +T

(kp̃rf (τ |tk+1 )k2Q − kp̃∗rf (τ |tk )k2Q )dτ,

tk+1
Z +T

∆V3 = −

1
 2q 2 ηδ
kp̃∗rf (τ |tk )k2 dτ  √ (e2aT − e2aδ ) 2
2a

Rewrite ∆V2 as

tZ
k +T

tk+1

∆V2 =

 12

q 2 η 2 δ 2 2aT
(e
− e2aδ )
2a
2q 2 ηδr T 2
1
1
≤ √
(
− T ) 2 (e2aT − e2aδ ) 2
δ
2a
q 2 η 2 δ 2 2aT
+
(e
− e2aδ ).
2a
+

in which
∆V1 =

tZ
k +T

tk+1

(81)

The difference of the value Lypunov function at tk and
tk+1 satisfies
∆V = V (prf (tk+1 )) − V (prf (tk ))
≤ J(p̃rf (tk+1 ), ũrf (tk+1 )) − J(p̃∗rf (tk ), ũ∗rf (tk ))

(84)

Applying Hölder inequality to the first term of the last
inequality yields

This proves the existence of K∞ functions α1 (·) and α2 (·)
satisfying
α1 (kprf (t)k) ≤ V (prf (t)) ≤ α2 (kprf (t)k).

q 2 η 2 δ 2 2aT
(e
− e2aδ ).
2a

+

[q 2 ηδea(τ +δ−tk+1 ) (2kp̃∗rf (τ |tk )k

∆V3 < −

tZk+1

tk

kp̃∗rf (τ |tk )k2Q dτ ≤ −qδε2 .

(88)

+ηδea(τ +δ−tk+1 ) )]dτ
tZ
k +T
=
2q 2 ηδea(τ +δ−tk+1 ) kp̃∗rf (τ |tk )k

In combination with (85),(87) and (88), the inequality
(82) thus satisfies

+q 2 η 2 δ 2 e2a(τ +δ−tk+1 ) dτ
tZ
k +T
≤
2q 2 ηδea(τ +δ−tk+1 ) kp̃∗rf (τ |tk )kdτ

q 2 η 2 δ 2 2aT
1
(e
− e2aδ )
∆V < −qδε2 + ηδeaT (r + ε) +
2
2a
1
1
2q 2 ηδr T 2
+ √
(
(89)
− T ) 2 (e2aT − e2aδ ) 2 .
δ
2a

tk+1

tk+1

11

horizon and the sampling period are set to be T = 2
s and δ = 0.2 s, respectively. The positive define
matrices P and Q are chosen, according to Lemma 1, as
P = diag{0.4, 0.4} and Q = diag{0.2, 0.2}, respectively.
The feedback gains for the terminal controller are set to
be k̃1 = k̃2 = 1.2 to satisfy the requirements given by
Lemma 1.

From (77), ∆V < 0 holds. It follows from Theorem 2 of
[34] that kp̃∗rf (tk |tk )k ≤ ε for tk ≥ tf , where tf > t0 is a
finite time instant. When the tracking error enters into
the terminal region, i.e., prf (tk ) ∈ Ωε , reconsider ∆V1
and ∆V3 :
2q 2 ηδεea(τ −tk+2 ) dτ +

tk+1
2

=

q 2 η 2 δ 2 2aT
(e
− e2aδ )
2a

0.1

2 2 2

2q ηδε aT
q η δ 2aT
(e − eaδ ) +
(e
− e2aδ ).
a
2a

xrf

0

p rf

∆V1 ≤

tZ
k +T

× 10-3

yrf

1

-0.1

0
-0.2

Due to the decreasing property of kp̃∗rf (τ |tk )k2Q in Ωε , it
follows that

-1
24.8

-0.3
0
× 10

∆V3 ≤ −qδkp∗rf (tk+1 |tk )k2 .

5

10

15

24.86
20

25

30

35

40

-4

5

xfe
yfe

p fe

Since kp̃∗rf (τ |tk+1 )k ≤ kp̃∗rf (τ |tk )k + ηδeaT , we have
kprf (tk+1 )k2 ≤ kp∗rf (tk+1 |tk )k2 + η 2 δ 2 e2aδ + 2εηδeaδ .
Consequently,

0

-5
0

∆V3 ≤ −qδkprf (tk+1 )k2 + qη 2 δ 3 e2aδ + 2qεηδ 2 eaδ .

5

10

15

20

25

30

35

Time (sec)

Fig. 3. Tracking errors by using tube-MPC.

As a result, it holds that
2

∆V ≤ −qδkprf (tk+1 )k + σ(η),
2

First, let us design tube-MPC according to Lemma 3
and Theorem 1. We set the feedback gain to be K =
diag{−2.3, −2.3}. The control input constraint for Problem 1 is Utube = λtube U with λtube = 0.6636, and the
terminal region is given by Ωtube = {p̃rf : |x̃rf | + |ỹrf | ≤
0.0542}. To evaluate the tracking performance, we take
the tracking error prf and the real position deviation
from the center of the tube pf e as indexes. Applying
Algorithm 1 to the tracking system yields the tracking
performance as shown in Fig. 3. It can be found, from
Fig. 3, that the tracking error converges to a neighbourhood of the origin, and the trajectory of the follower
lies in the tube T = p∗f h ⊕ Pf e with Pf e = {|pf e | ≤
"
#
0.0017
}, which is obtained from Lemma 3. Fig. 4
0.0017
shows the control input performance of the follower. We
also take |v|/a + |ω|/b as an index to evaluate the input
constraint. The fluctuated control signal indicates the
effectiveness of the feedback part of the controller which
reduces the tracking error caused by the disturbances.
Furthermore, the input constraint Utube for the nominal
system is active over the interval t ∈ [0, 3], while the
constraint U for the real system is not active, which
indicates the weak control ability.

(90)
2 2 2

η δ
(e2aT − e2aδ ) +
where σ(η) = 2q aηδε (eaT − eaδ ) + q 2a
1
aT
(ε + r) + qη 2 δ 3 e2aδ + 2qεηδ 2 eaδ is obviously a K
2 ησe
function with respect to η. Hence, the theorem is proved.
✷

5

Simulation results

The simulation is implemented on a PC equipped with
a dual-core 3.20 GHz Intel i5 CPU, 7.88 GB RAM and
64-bit Windows 10 operating system. The optimization
problem is transcribed by Tool Box ICLOCS (Imperical
College London Optimal Control Software, see [35]), 1.2
version, and solved by NLP (Nonlinear Programming)
solver IPOPT (Interior Point OPTimizer, see [36]),
3.11.8 version.
The mechanism parameters of the two homogeneous
robots used in the simulation are taken from an
educational robot named E-puck [37], and are given by
a = 0.13 m/s, ρ = 0.0267 m and b = a/ρ = 4.8598 rad/s.
The trajectory to be tracked is a circular motion with
linear velocity vr = 0.015 m/s, angular velocity ωr =
0.04 m/s and initial configuration ξr (0) = [0, 0, π3 ]T .
The initial configuration of the follower is set to be
ξf h = [0.2, −0.2, − π2 ]T . The disturbances are bounded
by η = 0.004. For the tracking objective, the prediction

To show the effect of different choices of the feedback
gain K on the tracking performance, we set K =
diag{−1, −1} and K = diag{−4, −4}, respectively, to
observe the difference between the actual trajectory
and the optimal one. As shown in Fig. 5, increasing of

12

2

0.2
0.1
0
-0.1

Inputs

1

disturbances due to the open-loop control during each
sampling period. Fig. 7 shows the control input under
NRMPC. According to Algorithm 2, the control signal at
each time instant is optimal corresponding to its current
state, which indicates its robustness. We also note that
the input constraint is active over the interval t ∈ [0, 1.5],
which demonstrates a better tracking capability.

Linear velocity v(m/s)
Angular velocity ω(rad/s)

18 18.5 19

0
-1
-2
0

5

10

15

20

25

30

35

40
0.1
0

p rf

|v|/a+| ω |/b

1
3
2
1
0
-1

-0.1

0.5
-0.2

26

-0.3

0
0

5

10

15

20

25

30

35

0

40

Time(sec)
1

Fig. 4. Control input of tube-MPC.

5

10

xrf

× 10-3

15

yrf

26.1
20

25

30

35

× 10-3
xfe

p fe

0.5

K reduces the difference of the actual position and the
center of the tube and therefore improves the tracking
performance.

40

yfe

0
-0.5
-1
0

2

p fe

10

15

20

25

30

35

40

Time (sec)

K = diag{-1, -1}

1

xfe

Fig. 6. Tracking errors by using NRMPC.

yfe

0
3

-1

5

10

15

20

25

30

35

Inputs

0

40

Time (sec)
2

0.2

× 10-4

Linear velocity v(m/s)
Angular velocity ω(rad/s)

0.1

2

-2

0
1
18.2 18.6 19
0

K = diag{-4, -4}

xfe

1

p fe

5

× 10-3

-1

yfe

0

5

10

15

0

5

10

15

20

25

30

35

40

20

25

30

35

40

0
-1

|v|/a+| ω |/b

1

-2
0

5

10

15

20

25

30

35

40

Time (sec)

Fig. 5. Real-time position deviation from the center of the
tube pf e with different feedback gains.

0.5

0

Time(sec)

Next, design NRMPC according to Lemma 4, Lemma 5
and Theorem 2. The input constraint of NRMPC differs from the constraint of tube-MPC and is given by
ũf ∈ U according to Algorithm 2. The terminal region
is designed as Ωε = {p̃rf : kp̃rf k ≤ 0.063}, and consequently ε = 0.063, which satisfies the conditions in
Lemma 4, Lemma 5 and Theorem 2. Fig. 6 presents
the tracking performance of Algorithm 2. It can be observed that the tracking error converges to a neighbourhood of the origin. To compare the influence level by
disturbances of the two strategies proposed, we define
pf e (tk ) = pf (tk ) − p̃∗f (tk ) in NRMPC. It can be seen
that the tracking performance is directly influenced by

Fig. 7. Control input of NRMPC.

To further compare tube-MPC with NRMPC, we take
cost function J, real stage cost kprf k2Q + kurf k2P , real
state cost kprf k2Q and real input cost kurf k2P to evaluate
the converging performance. Their cost curves are plotted in Fig 8. As it can be seen, the total cost, the stage
cost and the state cost decrease faster by implementing
NRMPC than by tube-MPC. However, the input cost
of NRMPC is higher than that of tube-MPC. This is
explained by the fact that the input constraint of tubeMPC is tighter than that of NRMPC, which may degrade the control capability. This also helps explaining

13

Optimization time (sec)

Optimization time (sec)

why the tracking error decreases faster by NRMPC than
by tube-MPC. Fig. 9 provides the computation time in
solving the optimization problems. It shows that there
is no significant difference between NRMPC and tubeMPC, which implies that they have almost the same
computational complexity. However, as stated in Remark 4, the optimization problem can be solved off-line
in tube-MPC, whereas the optimization problem has to
be solved on-line in NRMPC.
Finally, we summarize the simulation study as follows:
(i) Tube-MPC presents a better steady state performance
than NRMPC. This is because the control strategy
of tube-MPC consists of two parts: optimal control
and feedback part, while NRMPC is open-loop control
during each sampling period.
(ii) NRMPC performs better in terms of dynamic property than tube-MPC due to the tighter input constraint of tube-MPC.
(iii) The computational complexities of tube-MPC
and NRMPC are almost the same. However, the
optimization problem in tube-MPC can be solved offline, which may enhance its real-time performance.
Total cost
Tube-MPC
NRMPC

J

0.15
0.1
0.05

5

15

0

3

kurf kP

Tube-MPC
NRMPC

0.03
0.02
0.01
0

× 10-3

5

10

5

10

Time(sec)

15

10

15

20

25

30

35

40

30

35

40

Sampling time (sec)
NRMPC

5
4
3
2
1
0
0

5

10

15

20

25

[4] Z. P. Jiang and H. Nijmeijer.
Tracking control of
mobile robots: a case study in backstepping. Automatica,
33(7):1393–1399, 1997.

1

0

5

10

[5] J. Yang and J. Kim. Sliding mode control for trajectory
tracking of nonholonomic wheeled mobile robots. IEEE
Transactions on Robotics and Automation, 15(3):578–587,
1999.

15

Time(sec)

Fig. 8. Costs of tube-MPC and NRMPC.

6

5

[3] G. Lafferriere, A. Williams, J. Caughman, and J.J.P.
Veerman.
Decentralized control of vehicle formations.
Systems & control letters, 54(9):899–910, 2005.

Tube-MPC
NRMPC

0
0

0

[2] P. Wang and B. Ding. Distributed RHC for tracking and
formation of nonholonomic multi-vehicle systems. IEEE
Transactions on Automatic Control, 59(6):1439–1453, 2014.

15

Time(sec)
Input cost

2

0

[1] J. Simanek, M. Reinstein, and V. Kubelka. Evaluation of
the EKF-based estimation architectures for data fusion in
mobile robots. IEEE/ASME Transactions on Mechatronics,
20(2):985–990, 2015.

0.01

Time(sec)
State cost

0.04

kprf kQ

10

1

References

0
0

2

step. It was shown that the tracking control system is ISS
under both robust MPC strategies, and their optimization problems are feasible. Simulation results illustrated
the effectiveness of the two schemes and their respective
advantages.

0.02

0

3

Fig. 9. Computation time of tube-MPC and NRMPC.

Tube-MPC
NRMPC

0.03

4

Sampling time (sec)

Stage cost

0.04

kprf kQ + kurf kP

0.2

Tube-MPC

5

[6] T. Lee, K. Song, C. Lee, and C. Teng.
Tracking
control of unicycle-modeled mobile robots using a saturation
feedback controller. IEEE Transactions on Control Systems
Technology, 9(2):305–318, 2001.

Conclusion

[7] J. A. Marshall, M. E. Broucke, and B. A. Francis. Pursuit
formations of unicycles. Automatica, 42(1):3–12, 2006.

In this paper, two robust MPC strategies have been
developed for tracking of nonholonomic systems with
coupled input constraint and bounded disturbances. We
first developed a tube-MPC strategy, where the trajectory of the real system is constrained in a tube centered
along the optimal nominal trajectory by a nonlinear
feedback law based on the deviation of the actual states
from the optimal nominal states. Tube-MPC possesses
robustness but sacrifices optimality, thus we further developed the NRMPC scheme, where the state of the
nominal system is updated by the actual state at each

[8] J. Ghommam, H. Mehrjerdi, M. Saad, and F. Mnif.
Formation path following control of unicycle-type mobile
robots. Robotics and Autonomous Systems, 58(5):727–736,
2010.
[9] D. Gu and H. Hu. Receding horizon tracking control of
wheeled mobile robots. IEEE Transactions on Control
Systems Technology, 14(4):743–749, 2006.
[10] D. Q. Mayne, J. B. Rawlings, C. V. Rao, and P. O. M.
Scokaert. Constrained model predictive control: Stability and
optimality. Automatica, 36(6):789–814, 2000.

14

[11] J. Chen, D. Sun, J. Yang, and H. Chen. Leader-follower
formation control of multiple nonholonomic mobile robots
incorporating receding-horizon scheme. The International
Journal of Robotics Research, 29:727–747, 2010.

to lipschitz nonlinear systems. Systems & Control Letters,
62(2):194–200, 2013.
[29] J. Fleming, B. Kouvaritakis, and M. Cannon. Robust tube
MPC for linear systems with multiplicative uncertainty.
IEEE Transactions on Automatic Control, 60(4):1087–1092,
2015.

[12] Z. Sun and Y. Xia. Receding horizon tracking control of
unicycle-type robots based on virtual structure. International
Journal of Robust and Nonlinear Control, 2016, DOI:
10.1002/rnc.3555.

[30] D. L. Marruedo, T. Alamo, and E. F. Camacho. Inputto-state stable MPC for constrained discrete-time nonlinear
systems with bounded additive uncertainties. In Proceedings
of IEEE Conference on Decision and Control, volume 4,
pages 4619–4624, 2002.

[13] J. B. Rawlings and D. Q. Mayne. Model predictive control:
Theory and design. Nob Hill Pub., 2009.
[14] P. O. M. Scokaert and J. B. Rawlings. Stability of model
predictive control under perturbations. In Proceedings of
the IFAC Symposium on Nonlinear Control Systems Design,
pages 1317–1322, 1995.

[31] A. Richards and J. How. Robust stable model predictive
control with constraint tightening.
In Proceedings of
American Control Conference, pages 1557–1562, 2006.
[32] H. Li and Y. Shi. Robust distributed model predictive
control of constrained continuous-time nonlinear systems:
a robustness constraint approach. IEEE Transactions on
Automatic Control, 59(6):1673–1678, 2014.

[15] D. L. Marruedo, T. Alamo, and E. F. Camacho. Stability
analysis of systems with bounded additive uncertainties
based on invariant sets: Stability and feasibility of MPC. In
Proceedings of American Control Conference, pages 364–369,
2002.

[33] E. D. Sontag and Y. Wang. On characterizations of the
input-to-state stability property. Systems & Control Letters,
24(5):351–359, 1995.

[16] M. V. Kothare, V. Balakrishnan, and M. Morari. Robust
constrained model predictive control using linear matrix
inequalities. Automatica, 32(10):1361–1379, 1996.

[34] H. Michalska and D. Q. Mayne. Robust receding horizon
control of constrained nonlinear systems. IEEE Transactions
on Automatic Control, 38(11):1623–1633, 1993.

[17] J. H. Lee and Z. Yu. Worst-case formulations of model
predictive control for systems with bounded parameters.
Automatica, 33(5):763–781, 1997.

[35] P.
Falugi,
E.
Kerrigan,
and
E. v. Wyk. Imperial college london optimal control software
(iclocs). http://www.ee.ic.ac.uk/ICLOCS/.

[18] Z. Wan and M. V. Kothare. Robust output feedback model
predictive control using off-line linear matrix inequalities.
Journal of Process Control, 12(7):763–774, 2002.

[36] A. Wächter and L. T. Biegler. On the implementation
of an interior-point filter line-search algorithm for largescale nonlinear programming. Mathematical Programming,
106(1):25–57, 2006.

[19] L. Magni, G. De Nicolao, R. Scattolini, and F. Allgöwer.
Robust model predictive control for nonlinear discrete-time
systems. International Journal of Robust and Nonlinear
Control, 13(3-4):229–246, 2003.

[37] F. Mondada, M. Bonani, X. Raemy, J. Pugh, C. Cianci,
A. Klaptocz, S. Magnenat, J. Zufferey, D. Floreano, and
A. Martinoli. The e-puck, a robot designed for education
in engineering. In Proceedings of the 9th conference on
autonomous robot systems and competitions, volume 1, pages
59–65, 2009.

[20] H. Chen, C. W. Scherer, and F. Allgöwer. A game theoretic
approach to nonlinear robust receding horizon control of
constrained systems. In Proceedings of American control
conference, volume 5, pages 3073–3077, 1997.
[21] D. Limón, T. Alamo, F. Salas, and E. F. Camacho. Input
to state stability of min–max MPC controllers for nonlinear
systems with bounded uncertainties. Automatica, 42(5):797–
803, 2006.
[22] D. M. Raimondo, D. Limon, M. Lazar, L. Magni, and
E. F. Ndez Camacho. Min-max model predictive control of
nonlinear systems: A unifying overview on stability. European
Journal of Control, 15(1):5–21, 2009.
[23] D. Q. Mayne and W. Langson. Robustifying model predictive
control of constrained linear systems. Electronics Letters,
37(23):1422–1423, 2001.
[24] L. Chisci, J. A. Rossiter, and G. Zappa. Systems with
persistent disturbances: predictive control with restricted
constraints. Automatica, 37(7):1019–1028, 2001.
[25] W. Langson, I. Chryssochoos, S. V. Raković, and D. Q.
Mayne.
Robust model predictive control using tubes.
Automatica, 40(1):125–133, 2004.
[26] D. Q. Mayne, M. M. Seron, and S. V. Raković. Robust
model predictive control of constrained linear systems with
bounded disturbances. Automatica, 41(2):219–224, 2005.
[27] D. Q. Mayne, E. C. Kerrigan, E. J. Van Wyk, and
P. Falugi. Tube-based robust nonlinear model predictive
control. International Journal of Robust and Nonlinear
Control, 21(11):1341–1353, 2011.
[28] S. Yu, C. Maier, H. Chen, and F. Allgöwer. Tube MPC
scheme based on robust control invariant set with application

15

Optimizatin time (sec)

NRMPC

5
4
3
2
1
0

Optimizatin time (sec)

0

5

10

15

20

25

30

35

40

25

30

35

40

Tube-MPC

5
4
3
2
1
0
0

5

10

15

20

Sampling time (sec)

