0018-9286 (c) 2020 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TAC.2020.3038402, IEEE
Transactions on Automatic Control
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. XX, NO. XX, 2021
1
Robust Controller Design for Stochastic
Nonlinear Systems via Convex Optimization
Hiroyasu Tsukamoto,
Graduate Student Member, IEEE
, and Soon-Jo Chung,
Senior Member, IEEE
Abstract
—
This paper presents C
onV
ex optimization-
based S
tochastic steady-state T
racking E
rror M
inimization
(CV-STEM), a new state feedback control framework for a
class of It
ˆ
o stochastic nonlinear systems and Lagrangian
systems. Its innovation lies in computing the control input
by an optimal contraction metric, which greedily minimizes
an upper bound of the steady-state mean squared tracking
error of the system trajectories. Although the problem of
minimizing the bound is non-convex, its equivalent convex
formulation is proposed utilizing state-dependent coeffi-
cient parameterizations of the nonlinear system equation.
It is shown using stochastic incremental contraction anal-
ysis that the CV-STEM provides a sufficient guarantee for
exponential boundedness of the error for all time with
L
2
-robustness properties. For the sake of its sampling-
based implementation, we present discrete-time stochastic
contraction analysis with respect to a state- and time-
dependent metric along with its explicit connection to
continuous-time cases. We validate the superiority of the
CV-STEM to PID,
H
∞
, and baseline nonlinear controllers for
spacecraft attitude control and synchronization problems.
Index Terms
—
Stochastic optimal control, Optimization
algorithms, Robust control, Nonlinear systems, LMIs.
I. I
NTRODUCTION
S
TABLE and optimal feedback control of It
ˆ
o stochastic
nonlinear systems [1] is an important, yet challenging
problem in designing autonomous robotic explorers operat-
ing with sensor noise and external disturbances. Since the
probability density function of stochastic processes governed
by It
ˆ
o stochastic differential equations exhibits non-Gaussian
behavior characterized by the Fokker-Plank equation [1],
[2], feedback control schemes developed for deterministic
nonlinear systems could fail to meet control performance
specifications in the presence of stochastic disturbances.
A. Contributions
The main purpose of this paper is to propose C
onV
ex
optimization-based S
tochastic steady-state T
racking E
rror
M
inimization (CV-STEM), a new framework to design an op-
timal contraction metric for feedback control of It
ˆ
o stochastic
nonlinear systems and stochastic Lagrangian systems as in
Fig. 1. Contrary to Lyapunov theory, which gives a sufficient
condition for exponential convergence, the existence of a
The authors are with the Graduate Aerospace Laboratories (GALCIT),
California Institute of Technology, 1200 E California Blvd, Pasadena, CA,
USA. E-mail:
{
htsukamoto, sjchung
}
@caltech.edu
, Code:
https://github.com/astrohiro/cvstem
.
Fig. 1.
Illustration of the CV-STEM control:
M
(
x,t
)
denotes
the optimal contraction metric for the differential Lyapunov function
δx
>
M
(
x,t
)
δx
;
x
(
t
)
and
x
d
(
t
)
are controlled and desired system
trajectories;
u
(
t
)
is the control input computed by
M
(
x,t
)
(see Sec. III
for details).
contraction metric leads to a necessary and sufficient char-
acterization of exponential incremental stability of nonlinear
system trajectories [3], [4]. We explore this approach further
to obtain an optimal contraction metric for controlling It
ˆ
o
stochastic nonlinear systems. This paper builds upon our prior
work [5], but provides more rigorous proofs and explanations
on how we convexify the problem of minimizing
D
in Fig. 1
in a mean squared sense. We also investigate its stochastic
incremental stability properties and the impact of sampling-
based implementation on its control performance both in
detail, thereby introducing several additional theorems and
simulation results. The construction and contributions of our
CV-STEM method are summarized as follows.
1) The CV-STEM design is based on a convex combination
of multiple State-Dependent Coefficient (SDC) forms of a
nonlinear system equation (i.e.
f
(
x,t
)
written as
A
(
x,t
)
x
[6]–
[8], where
A
(
x,t
)
is not necessarily unique). The main ad-
vantage of our control synthesis algorithm lies in solving an
optimization problem, the objective of which is to find an
optimal contraction metric that greedily minimizes an upper
bound of the steady-state mean squared tracking error of
It
ˆ
o stochastic nonlinear system trajectories, constructing an
optimal feedback control gain and Control Lyapunov Function
(CLF) [9]–[11] (see Fig. 1). Although the problem of minimiz-
ing the bound is originally non-convex, we reformulate it as a
convex optimization problem with the State-Dependent Riccati
Inequality (SDRI) constraint expressed as an LMI [12], which
can be solved by various computationally-efficient numerical
methods [12]–[15]. We also propose one way to utilize non-
unique choices of SDC forms for verifying the controllability
Authorized licensed use limited to: CALIFORNIA INSTITUTE OF TECHNOLOGY. Downloaded on November 19,2020 at 16:41:38 UTC from IEEE Xplore. Restrictions apply.
0018-9286 (c) 2020 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TAC.2020.3038402, IEEE
Transactions on Automatic Control
2
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. XX, NO. XX, 2021
of the system. This result is a significant improvement over
the observer design [16], whose optimization-cost function
uses a linear combination of observer parameters without
accounting for the contraction constraint, which we express as
an LMI [12] in this paper. This approach is further extended to
the control of stochastic Lagrangian systems with a nominal
exponentially stabilizing controller, and its superiority to the
prior work [17], [18], PID, and
H
∞
control [19]–[21] is shown
using results of numerical simulations on spacecraft attitude
control and synchronization.
2) It is proven using stochastic incremental contraction anal-
ysis that any solution trajectory under the CV-STEM feedback
control exponentially converges to the desired trajectory in a
mean squared sense with a non-vanishing error term (which
will be minimized as explained above). It is also shown that the
controller is robust against external deterministic disturbances
which often appear in parametric uncertain systems, and that
the tracking error has a finite
L
2
gain with respect to the
noise and disturbances acting on the system. We note that the
mean-square bound does not imply the asymptotic almost-sure
bounds although finite time bounds could be obtained [1],
[22], as the CV-STEM-based Lyapunov function is not a
supermartingale due to the non-vanishing steady-state error
term.
3) Discrete-time stochastic incremental contraction analysis
with respect to a state- and time-dependent metric is derived
for studying the effect of sampling-based implementation
of the CV-STEM on its control performance. It is proven
that stochastic incremental stability of discrete-time systems
reduces to that of continuous-time systems if the time interval
is sufficiently small. It is shown in the numerical simulations
that the CV-STEM sampling period
∆
t
can be relaxed to
∆
t
≤
25
(s) for spacecraft attitude control and
∆
t
≤
350
(s)
for spacecraft tracking and synchronization control without
impairing its performance.
4) Some extensions of the CV-STEM are derived to ex-
plicitly incorporate input constraints and to avoid solving the
convex optimization problem at every time instant.
B. Related Work
CLFs [9]–[11] as well as feedback linearization [11], [23],
[24] are among the most widely used tools for controlling
nonlinear systems perturbed by deterministic disturbances.
Since there is no general analytical scheme for finding a
CLF, several techniques are proposed to find them utilizing
some special structure of the systems in question [25]–[29].
The state-dependent Riccati equation method [6]–[8] can also
be viewed as one of these techniques and is applicable to
systems that are written in SDC linear structure. Building on
these ideas for deterministic systems, a stochastic counterpart
of the Lyapunov methods is proposed in [30] to design
CLF-based state and output feedback control of stochastic
nonlinear systems [31], [32]. For a class of strict-feedback and
output-feedback stochastic nonlinear systems, there exists a
more systematic way of asymptotic stabilization in probability
using a backstepping-based controller [33], [34]. However,
one drawback of these approaches is that they are primarily
directed toward stability with some implicit inverse optimality
guarantees.
Some theoretical methodologies have been developed to
explicitly incorporate optimality into their feedback control
formulation. These include
H
∞
control [20], [21], [35], which
attempts to minimize the
H
∞
norm for the sake of optimal
disturbance attenuation. Although it is originally devised for
linear systems [36]–[41], its nonlinear analogues are obtained
in [20], [21] and then expanded to stochastic nonlinear
systems [19] unifying the results on the
L
2
gain analysis
based on the Hamilton-Jacobi equations and inequalities [11].
Although we could design feedback control schemes optimally
for specific types of systems such as Hamiltonian systems
with stochastic disturbances [42] or linearized and discretized
stochastic nonlinear systems [43], finding the solution to
the stochastic nonlinear state feedback
H
∞
optimal control
problem is not trivial in general.
The CV-STEM addresses this issue by numerically sampling
an optimal contraction metric and CLF that greedily minimize
an upper bound of the steady-state mean squared tracking
error of It
ˆ
o stochastic nonlinear system trajectories. We select
this as an objective function, instead of integral objective
functions which often appear in optimal control problems,
as it gives us an exact convex optimization-based control
synthesis algorithm. Also, since the problem has the SDRI
as its constraint, the CV-STEM control is robust against both
deterministic and stochastic disturbances and ensures that the
tracking error is exponentially bounded for all time. We remark
that this approach is not intended to supersede but to be
utilized on top of existing methodologies on constructing de-
sired control inputs using stochastic nonlinear optimal control
techniques [1], [44]–[47], as this is a type of feedback control
scheme. In particular, stochastic model predictive control [48],
[49] with guaranteed stability [50], [51]
assumes
the existence
of a stochastic CLF, whilst our approach explicitly
constructs
an optimal CLF which could be used for the stochastic CLF
with some modifications on the non-vanishing error term in
our formulation.
The tool we use for analyzing incremental stability [4] in
this paper is contraction analysis [3], [52], [53], where its
stochastic version is derived in [16], [22]. Contraction analysis
for discrete-time and hybrid systems is provided in [3], [54],
[55] and its stochastic counterpart is investigated in [56] with
respect to a state-independent metric. In this paper, we describe
discrete-time incremental contraction analysis with respect to
a state- and time-dependent metric. Since the differential (vir-
tual) dynamics of
δx
used in contraction analysis is a Linear
Time-Varying (LTV) system, global exponential stability can
be studied using a quadratic Lyapunov function of
δx
,
V
=
δx
>
M
(
x,t
)
δx
[3], as opposed to the Lyapunov technique
where
V
could be any function of
x
. Therefore, designing
V
reduces to finding a positive definite metric
M
(
x,t
)
[28], [57],
[58], which enables the aforementioned convex optimization-
based control of It
ˆ
o stochastic nonlinear systems.
C. Paper Organization
The rest of this paper is organized as follows. Section II
introduces stochastic incremental contraction analysis and
Authorized licensed use limited to: CALIFORNIA INSTITUTE OF TECHNOLOGY. Downloaded on November 19,2020 at 16:41:38 UTC from IEEE Xplore. Restrictions apply.
0018-9286 (c) 2020 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TAC.2020.3038402, IEEE
Transactions on Automatic Control
H. TSUKAMOTO
et al.
: ROBUST CONTROLLER DESIGN FOR STOCHASTIC NONLINEAR SYSTEMS VIA CONVEX OPTIMIZATION
3
presents its discrete-time version with a state- and time-
dependent metric. In Sec. III, the CV-STEM control for It
ˆ
o
stochastic nonlinear systems is presented and its stability is
analyzed using contraction analysis. In Sec. IV, this approach
is extended to the control of stochastic Lagrangian systems.
Section V elucidates several extensions of the CV-STEM con-
trol synthesis. The aforementioned two simulation examples
are reported in Sec. VI. Section VII concludes the paper.
D. Notation
For a vector
x
∈
R
n
and a matrix
A
∈
R
n
×
m
, we let
‖
x
‖
,
δx
,
∂
μ
x
,
‖
A
‖
,
‖
A
‖
F
,
Im(
A
)
,
Ker(
A
)
,
A
+
, and
κ
(
A
)
denote the Euclidean norm, infinitesimal variation of
x
, partial
derivative of
x
with respect to
μ
, induced 2-norm, Frobenius
norm, image of
A
, kernel of
A
, Moore–Penrose inverse, and
condition number, respectively. For a square matrix
A
, we
use the notation
λ
min
(
A
)
and
λ
max
(
A
)
for the minimum
and maximum eigenvalues of
A
,
Tr(
A
)
for the trace of
A
,
A
0
,
A
0
,
A
≺
0
, and
A
0
for the positive definite,
positive semi-definite, negative definite, negative semi-definite
matrices, respectively, and
sym(
A
) = (
A
+
A
>
)
/
2
. For a
vector
x
∈
R
n
and a positive definite matrix
A
∈
R
n
×
n
, we
denote a norm
√
x
>
Ax
as
‖
x
‖
A
. Also,
I
∈
R
n
×
n
represents
the identity matrix,
E
[
·
]
denotes the expected value operator,
and
E
ξ
[
·
]
denotes the conditional expected value operator
when
ξ
is given. The
L
p
norm in the extended space
L
pe
,
p
∈
[1
,
∞
]
, is defined as
‖
(
y
)
τ
‖
L
p
=
(
∫
τ
0
‖
y
(
t
)
‖
p
)
1
/p
<
∞
for
p
∈
[1
,
∞
)
and
‖
(
y
)
τ
‖
L
∞
= sup
t
≥
0
‖
(
y
(
t
))
τ
‖
<
∞
for
p
=
∞
, where
(
y
(
t
))
τ
is a truncation of
y
(
t
)
, i.e.,
(
y
(
t
))
τ
= 0
for
t > τ
and
(
y
(
t
))
τ
=
y
(
t
)
for
0
≤
t
≤
τ
with
τ
∈
[0
,
∞
)
.
II. S
TOCHASTIC
I
NCREMENTAL
S
TABILITY VIA
C
ONTRACTION
A
NALYSIS
We summarize contraction analysis that will be used for
stability analysis in the subsequent sections. This allows us to
utilize approaches for LTV systems theory, yielding a convex
optimization-based framework for optimal Lyapunov function
construction in Sec. III and IV.
We also present new theorems for analyzing stochastic
incremental stability of discrete-time nonlinear systems with
respect to a state- and time-dependent Riemannian metric,
along with its explicit connection to contraction analysis of
continuous-time systems.
A. Continuous-time Dynamical Systems
Consider the following continuous-time nonlinear non-
autonomous system and its virtual dynamics:
̇
x
=
f
(
x,t
)
, δ
̇
x
=
∂f
(
x,t
)
∂x
δx
(1)
where
t
∈
R
≥
0
,
x
:
R
≥
0
→
R
n
, and
f
:
R
n
×
R
≥
0
→
R
n
.
Incremental stability [4] is defined as stability of system
trajectories with respect to each other by means of differential
(virtual) dynamics. Contraction theory is used to study incre-
mental stability with exponential convergence.
Lemma 1:
The system (1) is contracting (i.e. all the solu-
tion trajectories exponentially converge to a single trajectory
globally from any initial condition), if there exists a uni-
formly positive definite metric
M
(
x,t
) = Θ(
x,t
)
>
Θ(
x,t
)
,
M
(
x,t
)
0
,
∀
x,t
, with a smooth coordinate transformation
of the virtual displacement
δz
= Θ(
x,t
)
δx
, such that
̇
M
(
x,t
) + 2 sym
(
M
(
x,t
)
∂f
∂x
)
−
2
γ
c
M
(
x,t
)
,
∀
x,t
(2)
where
γ
c
>
0
. If the system (1) is contracting, then we have
‖
δz
(
t
)
‖
=
‖
Θ(
x,t
)
δx
(
t
)
‖≤‖
δz
(0)
‖
e
−
γ
c
t
.
Proof:
See [3].
Next, consider the nonlinear system (1) with stochastic per-
turbation given by the It
ˆ
o stochastic differential equation
dx
=
f
(
x,t
)
dt
+
G
(
x,t
)
d
W
, x
(0) =
x
0
(3)
where
G
:
R
n
×
R
≥
0
→
R
n
×
d
is a matrix-valued function,
W
(
t
)
is a
d
-dimensional Wiener process, and
x
0
is a ran-
dom variable independent of
W
(
t
)
[59]. In this paper, we
assume that
∃
L
1
>
0
s.t.
‖
f
(
x
1
,t
)
−
f
(
x
2
,t
)
‖
+
‖
G
(
x
1
,t
)
−
G
(
x
2
,t
)
‖
F
≤
L
1
‖
x
1
−
x
2
‖
,
∀
t
∈
R
≥
0
and
∀
x
1
,x
2
∈
R
n
, and
∃
L
2
>
0
s.t.
‖
f
(
x
1
,t
)
‖
2
+
‖
G
(
x
1
,t
)
‖
2
F
≤
L
2
(1 +
‖
x
1
‖
2
)
,
∀
t
∈
R
≥
0
and
∀
x
1
∈
R
n
, for the sake of existence and
uniqueness of the solution to (3). Now, consider the following
two systems with trajectories
ξ
1
(
t
)
and
ξ
2
(
t
)
driven by two
independent Wiener processes
W
1
(
t
)
and
W
2
(
t
)
:
dξ
=
[
f
(
ξ
1
,t
)
f
(
ξ
2
,t
)
]
dt
+
[
G
1
(
ξ
1
,t
)
0
0
G
2
(
ξ
2
,t
)
][
d
W
1
d
W
2
]
(4)
where
ξ
(
t
) = [
ξ
1
(
t
)
>
,ξ
2
(
t
)
>
]
>
∈
R
2
n
. The following
theorem analyzes stochastic incremental stability of the two
trajectories
ξ
1
(
t
)
and
ξ
2
(
t
)
with respect to each other in
the presence of stochastic noise. The trajectories of (3) are
parameterized as
x
(0
,t
) =
ξ
1
and
x
(1
,t
) =
ξ
2
. Also, we
define
G
(
x,t
)
as
G
(
x
(0
,t
)
,t
) =
G
1
(
ξ
1
,t
)
and
G
(
x
(1
,t
)
,t
) =
G
2
(
ξ
2
,t
)
.
Theorem 1:
Suppose that there exist bounded positive con-
stants
m
,
m
,
g
1
,
g
2
,
m
x
, and
m
x
2
s.t.
m
≤ ‖
M
(
x,t
)
‖ ≤
m
,
‖
G
1
(
x,t
)
‖
F
≤
g
1
,
‖
G
2
(
x,t
)
‖
F
≤
g
2
,
‖
∂
(
M
ij
)
/∂x
‖ ≤
m
x
, and
∥
∥
∂
2
(
M
ij
)
/∂x
2
∥
∥
≤
m
x
2
,
∀
x,t
. Suppose also that
(2) holds (i.e., the deterministic system (1) is contracting).
Consider the generalized squared length with respect to a
Riemannian metric
M
(
x
(
μ,t
)
,t
)
defined by
V
(
x,∂
μ
x,t
) =
∫
1
0
∂x
∂μ
>
M
(
x
(
μ,t
)
,t
)
∂x
∂μ
dμ
(5)
s.t.
V
(
x,∂
μ
x,t
)
≥
m
‖
ξ
1
−
ξ
2
‖
2
. Then we have
L
V
≤−
2
γ
1
V
+
m
C
c
(6)
for
γ
1
=
γ
c
−
((
g
2
1
+
g
2
2
)
/
2
m
)(
ε
c
m
x
+
m
x
2
/
2)
and
C
c
=
(
m/m
+
m
x
/
(
ε
c
m
))(
g
2
1
+
g
2
2
)
, where
L
is an infinitesimal
differential generator [16],
γ
c
is the contraction rate for the
deterministic system (1), and
ε
c
>
0
is an arbitrary constant.
Further, if we have
γ
1
>
0
, (6) implies that the mean squared
distance between the two trajectories of (4), whose initial
conditions given by a probability distribution
p
(
a
0
,b
0
)
that are
independent of
W
1
(
t
)
and
W
2
(
t
)
, is exponentially bounded as
follows:
E
[
‖
ξ
1
(
t
)
−
ξ
2
(
t
)
‖
2
]
≤
C
c
2
γ
1
+
E
[
V
(
x
(0)
,∂
μ
x
(0)
,
0)]
e
−
2
γ
1
t
m
.
(7)
Authorized licensed use limited to: CALIFORNIA INSTITUTE OF TECHNOLOGY. Downloaded on November 19,2020 at 16:41:38 UTC from IEEE Xplore. Restrictions apply.
0018-9286 (c) 2020 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TAC.2020.3038402, IEEE
Transactions on Automatic Control
4
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. XX, NO. XX, 2021
Proof:
Using the property
Tr(
AB
)
≤ ‖
A
‖
Tr(
B
)
for
A,B
0
, we have
Tr(
G
i
(
ξ
i
,t
)
>
M
(
ξ
i
,t
)
G
i
(
ξ
i
,t
))
≤
mg
2
i
.
Therefore, computing
L
V
as in the proof given in Lemma 2
of [16] yields (6). Taking expectation on both sides of (6) along
with Dynkin’s formula [1, pp. 10] completes the derivation of
(7).
Remark 1:
The contraction rate
γ
1
and uncertainty bound
C
c
depend on the choice of an arbitrary constant
ε
c
. One way
to select
ε
c
is to solve
dF/dε
c
= 0
with
F
(
ε
c
) =
C
c
/
(2
γ
1
)
,
whose solution minimizes the steady-state bound
F
(
ε
c
)
with
the constraint
γ
1
>
0
[16]. Line search algorithms could also
be used to select their optimal values [57], [58]. We will utilize
the fact that
C
c
is a function of
m/m
to facilitate the convex
optimization-based control synthesis in Sec. III and IV.
B. Main Result 1: Connection between Continuous and
Discrete Stochastic Incremental Contraction Analysis
We establish a similar result to Lemma 1 for the following
discrete-time nonlinear system and its virtual dynamics:
x
k
+1
=
f
k
(
x
k
,k
)
, δx
k
+1
=
∂f
k
(
x
k
,k
)
∂x
k
δx
k
(8)
where
x
k
∈
R
n
and
f
k
:
R
n
×
N
→
R
n
.
Lemma 2:
The system (8) is contracting if there ex-
ists a uniformly positive definite metric
M
k
(
x
k
,k
) =
Θ
k
(
x
k
,k
)
>
Θ
k
(
x
k
,k
)
,
M
k
(
x
k
,k
)
0
,
∀
x
k
,k
, with a smooth
coordinate transformation of the virtual displacement
δz
k
=
Θ
k
(
x
k
,k
)
δx
k
s.t.
∂f
k
∂x
k
>
M
k
+1
(
x
k
+1
,k
+ 1)
∂f
k
∂x
k
(1
−
γ
d
)
M
k
(
x
k
,k
)
,
∀
x
k
,k
(9)
where
γ
d
∈
(0
,
1)
. If the system (8) is contracting, then we
have
‖
δz
k
‖
=
‖
Θ
k
(
x
k
,k
)
δx
k
‖≤‖
δz
0
‖
(1
−
γ
d
)
k
2
.
Proof:
See [3], [55].
We now present a discrete-time version of Theorem 1, which
can be extensively used for proving stability of discrete-time
and hybrid stochastic nonlinear systems, along with known
results for deterministic systems [54], [55]. Consider the
discrete-time nonlinear system (8) with stochastic perturbation
modeled by the stochastic difference equation
x
k
+1
=
f
k
(
x
k
,k
) +
G
k
(
x
k
,k
)
w
k
(10)
where
G
k
:
R
n
×
N
→
R
n
×
d
is a matrix-valued function and
w
k
is a
d
-dimensional sequence of zero mean uncorrelated
normalized Gaussian random variables. Consider the following
two systems with trajectories
ξ
1
,k
and
ξ
2
,k
driven by two
independent stochastic perturbation
w
1
,k
and
w
2
,k
:
ξ
k
+1
=
[
f
k
(
ξ
1
,k
,k
)
f
k
(
ξ
2
,k
,k
)
]
+
[
G
1
,k
(
ξ
1
,k
,k
)
0
0
G
2
,k
(
ξ
2
,k
,k
)
][
w
1
,k
w
2
,k
]
(11)
where
ξ
k
= [
ξ
>
1
,k
,ξ
>
2
,k
]
>
∈
R
2
n
. The following theorem
analyzes stochastic incremental stability for discrete-time non-
linear systems, but we remark that this is different from [56],
[60] in that the stability is studied in a differential sense
and its Riemannian metric is state- and time-dependent. We
parameterize
x
k
and
G
k
in (10) as
x
k
(
μ
= 0) =
ξ
1
,k
,
x
k
(
μ
= 1) =
ξ
2
,k
,
G
k
(
x
k
(
μ
= 0)
,k
) =
G
1
,k
(
ξ
1
,k
,k
)
, and
G
k
(
x
k
(
μ
= 1)
,k
) =
G
2
,k
(
ξ
2
,k
,k
)
.
Theorem 2:
Suppose that the system (11) has the following
bounds,
m
I
M
k
(
x
k
,k
)
mI,
∀
x
k
,k
,
‖
G
1
,k
(
ξ
1
,k
,k
)
‖
F
≤
g
1
d
, and
‖
G
2
,k
(
ξ
2
,k
,k
)
‖
F
≤
g
2
d
,
∀
ξ
1
,k
,ξ
2
,k
,k
, where
m
,
g
1
d
,
and
g
2
d
are bounded positive constants. Suppose also that (9)
holds for the discrete-time deterministic system (8) and there
exists
γ
2
∈
(0
,
1)
s.t.
γ
2
≤
1
−
(
m/m
)(1
−
γ
d
)
, where
γ
d
is the
contraction rate of (8). Consider the generalized squared length
with respect to a Riemannian metric
M
k
(
x
k
(
μ
)
,k
)
defined as
V
k
(
x
k
,∂
μ
x
k
,k
) =
∫
1
0
∂x
k
∂μ
>
M
k
(
x
k
(
μ
)
,k
)
∂x
k
∂μ
dμ
(12)
s.t.
V
k
(
x
k
,∂
μ
x
k
,k
)
≥
m
‖
ξ
1
,k
−
ξ
2
,k
‖
2
2
. Then the mean squared
distance between the two trajectories of the system (11) is
bounded as follows:
E
[
‖
ξ
1
,k
−
ξ
2
,k
‖
2
]
≤
1
−
̃
γ
k
d
1
−
̃
γ
d
C
d
+
̃
γ
k
d
m
E
[
V
0
(
x
0
,∂
μ
x
0
,
0)]
.
(13)
where
C
d
= (
m/m
)(
g
2
1
d
+
g
2
2
d
)
and
̃
γ
d
= 1
−
γ
2
∈
(0
,
1)
.
Proof:
Consider a Lyapunov-like function
V
k
in (12),
where we use
V
k
=
V
k
(
x
k
,∂
μ
x
k
,k
)
and
M
k
=
M
k
(
x
k
,k
)
for notational simplicity. Using the bounds along with (9) and
(10), we have, for
`
∈
N
, that
V
`
+1
≤
m
∫
1
0
∥
∥
∥
∥
∂f
`
∂x
`
∂x
`
∂μ
+
∂G
`
∂μ
w
`
∥
∥
∥
∥
2
dμ
(14)
≤
m
m
(1
−
γ
d
)
∫
1
0
∂x
`
∂μ
>
M
`
∂x
`
∂μ
dμ
+
m
∫
1
0
(
2
∂x
`
∂μ
>
∂f
`
∂x
`
>
∂G
`
∂μ
w
`
+
w
>
`
∂G
`
∂μ
>
∂G
`
∂μ
w
`
)
dμ
where
f
`
=
f
`
(
x
`
,`
)
and
G
`
=
G
`
(
x
`
,`
)
. Taking the
conditional expected value of
(14)
when
x
`
,
∂
μ
x
`
, and
`
are
given, we have that (see also: Theorem 2 of [60])
E
ζ
`
[
V
`
+1
]
≤
γ
m
V
`
+
mE
ζ
`
[
∫
1
0
w
>
`
∂G
`
∂μ
>
∂G
`
∂μ
w
`
dμ
]
≤
γ
m
V
`
+
∑
i
=1
,
2
mE
ζ
`
[
Tr
(
w
i,`
w
>
i,`
G
>
i,`
G
i,`
)]
≤
γ
m
V
`
+
m
∑
i
=1
,
2
Tr
(
G
>
i,`
G
i,`
)
≤
̃
γ
d
V
`
+
m
C
d
.
(15)
where
γ
m
=
m/m
(1
−
γ
d
)
, and
x
`
,
∂
μ
x
`
, and
`
are denoted
as
ζ
`
. Here, we used the condition:
∃
γ
2
∈
(0
,
1)
s.t.
γ
m
≤
1
−
γ
2
= ̃
γ
d
. Taking expectation over
ζ
`
−
1
in (15) with the
tower rule
E
ζ
`
−
1
[
V
`
+1
] =
E
ζ
`
−
1
[
E
ζ
`
[
V
`
+1
]]
gives us that
E
ζ
`
−
1
[
V
`
+1
]
≤
̃
γ
2
d
V
`
−
1
+
m
C
d
+
m
C
d
̃
γ
d
(16)
where
̃
γ
d
= 1
−
γ
2
. Continuing this operation with the relation
m
E
ζ
0
[
‖
ξ
1
,`
+1
−
ξ
2
,`
+1
‖
2
]
≤
E
ζ
0
[
V
`
+1
]
yields
E
ζ
0
[
‖
ξ
1
,k
−
ξ
2
,k
‖
2
]
−
̃
γ
k
d
m
V
0
≤
C
d
k
−
1
∑
i
=0
̃
γ
i
d
=
1
−
̃
γ
k
d
1
−
̃
γ
d
C
d
Authorized licensed use limited to: CALIFORNIA INSTITUTE OF TECHNOLOGY. Downloaded on November 19,2020 at 16:41:38 UTC from IEEE Xplore. Restrictions apply.
0018-9286 (c) 2020 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TAC.2020.3038402, IEEE
Transactions on Automatic Control
H. TSUKAMOTO
et al.
: ROBUST CONTROLLER DESIGN FOR STOCHASTIC NONLINEAR SYSTEMS VIA CONVEX OPTIMIZATION
5
where
k
=
`
+ 1
. Taking expectation over
ζ
0
and rearranging
terms result in (13).
Let us now consider the case where the time interval
∆
t
=
t
k
+1
−
t
k
is sufficiently small, i.e.,
∆
t
(∆
t
)
2
. Then the
continuous-time stochastic system (3) can be discretized as
x
k
+1
=
x
k
+
∫
t
k
+1
t
k
f
(
x
(
t
)
,t
)
dt
+
G
(
x
(
t
)
,t
)
d
W
(
t
)
'
x
k
+
f
(
x
k
,t
k
)∆
t
+
G
(
x
k
,t
k
)∆
W
k
(17)
where
x
k
=
x
(
t
k
)
,
∆
W
k
=
√
∆
tw
k
, and
w
k
is a
d
-
dimensional sequence of zero mean uncorrelated normalized
Gaussian random variables. When
∆
t
(∆
t
)
2
,
f
k
(
x
k
,k
)
and
G
k
(
x
k
,k
)
in (10) can be approximated as
f
k
(
x
k
,k
) =
x
k
+
f
(
x
k
,t
k
)∆
t
and
G
k
(
x
k
,k
) =
√
∆
tG
(
x
k
,t
k
)
. In this
situation, we have the following theorem that connects stochas-
tic incremental stability of discrete-time systems with that of
continuous-time systems.
Theorem 3:
Suppose that (15) in Theorem 2 holds with
̃
γ
d
= 1
−
γ
2
∈
(0
,
1)
. Then the expected value of
V
k
+1
up to
first order in
∆
t
is given as
E
ζ
k
[
V
k
+1
] =
V
k
+ ∆
t
L
V
k
, where
L
is an infinitesimal differential generator. Furthermore, the
following inequality holds:
L
V
k
(
x
k
,∂
μ
x
k
,t
k
)
≤−
γ
2
∆
t
V
k
(
x
k
,∂
μ
x
k
,t
k
) +
m
̃
C
c
(18)
where
̃
C
c
is a positive constant given as
̃
C
c
=
C
d
∆
t
=
m
m
∆
t
(
g
2
1
d
+
g
2
2
d
) =
m
m
(
g
2
1
+
g
2
2
)
(19)
with
g
1
and
g
2
defined in Theorem 1.
Proof:
M
k
+1
up to first order in
∆
t
is written as
M
k
+1
=
∂M
k
∂t
k
∆
t
+
n
∑
i
=1
∂M
k
∂
(
x
k
)
i
(
f
c,k
∆
t
+
G
c,k
∆
W
k
)
i
(20)
+
1
2
n
∑
i
=1
n
∑
j
=1
∂
2
M
k
∂
(
x
k
)
i
∂
(
x
k
)
j
(
G
c,k
∆
W
k
)
i
(
G
c,k
∆
W
k
)
j
+
M
k
where
f
c,k
and
G
c,k
are defined as
f
c,k
=
f
(
x
k
,t
k
)
and
G
c,k
=
G
(
x
k
,t
k
)
for notational simplicity. The subscripts
i
and
j
denote the
i
th and
j
th element of the corresponding
vectors. Similarly,
∂x
k
+1
/∂μ
up to first order in
∆
t
can be
computed as
∂x
k
+1
∂μ
=
∂x
k
∂μ
+
∂f
c,k
∂x
k
∂x
k
∂μ
∆
t
+
∂G
c,k
∂μ
∆
W
k
.
(21)
Substituting (20) and (21) into
E
ζ
k
[
V
k
+1
]
yields
E
ζ
k
[
V
k
+1
] =
E
ζ
k
[
∫
1
0
∂x
k
+1
∂μ
>
M
k
+1
∂x
k
+1
∂μ
dμ
]
=
V
k
+ (
dV
d,k
+
dV
s,k
)∆
t
+
O
(∆
t
3
/
2
)
(22)
where
dV
d,k
and
dV
s,k
are given by
dV
d,k
=
∫
1
0
∂x
k
∂μ
>
(
∂f
c,k
∂x
k
>
M
k
+
̇
M
k
+
M
k
∂f
c,k
∂x
k
)
∂x
k
∂μ
dμ
(23)
with
̇
M
k
=
∂M
k
/∂t
k
+
∑
n
i
=1
(
∂M
k
/∂
(
x
k
)
i
)
f
c,k
and
dV
s,k
=
∫
1
0
n
∑
i
=1
n
∑
j
=1
(
M
k
)
ij
(
∂G
c,k
∂μ
∂G
c,k
∂μ
>
)
ij
+2
∂
(
M
k
)
i
∂
(
x
k
)
j
∂x
k
∂μ
(
G
c,k
∂G
c,k
∂μ
>
)
ij
+
1
2
∂x
k
∂μ
>
∂
2
M
k
∂
(
x
k
)
i
∂
(
x
k
)
j
∂x
k
∂μ
(
G
c,k
G
>
c,k
)
ij
]
dμ.
(24)
We note that the properties of
w
k
as a
d
-dimensional sequence
of zero mean uncorrelated normalized Gaussian random vari-
ables are used to derive these relations. Since
dV
d,k
+
dV
s,k
=
L
V
k
where
L
is the infinitesimal differential generator,
we have
E
ζ
k
[
V
k
+1
] =
V
k
+ ∆
t
L
V
k
. Thus, the condition
E
ζ
k
[
V
k
+1
]
≤
(1
−
γ
2
)
V
k
+
m
C
d
given by (15) in Theorem 2
reduces to the following inequality:
L
V
k
(
x
k
,∂
μ
x
k
,t
k
)
≤−
γ
2
∆
t
V
k
(
x
k
,∂
μ
x
k
,t
k
) +
m
C
d
∆
t
.
(25)
Finally, (25) with the relations
̃
C
c
=
C
d
/
∆
t
and
G
k
(
x
k
,k
) =
√
∆
tG
(
x
k
,t
k
)
results in (18) and (19).
Remark 2:
The positive constant
̃
C
c
is equal to the positive
constant
C
c
in Theorem 1 when
m
x
= 0
. This is due to the
fact that we used an upper bound of
‖
M
k
‖
when obtaining
the first line of (14) in Theorem 2.
In practical control applications, we use the same control
input at
t
=
t
k
for a finite time interval
t
∈
[
t
k
,t
t
+1
)
. Theo-
rems 1 and 3 indicate that if
∆
t
is sufficiently small, a discrete-
time stochastic controller can be viewed as a continuous-
time counterpart with contraction rate
2
γ
1
=
γ
2
/
∆
t
. We will
illustrate how to select the sampling period
∆
t
large enough
without deteriorating the CV-STEM control performance in
Sec. VI. Also, the steady-state mean squared tracking error
for both discrete and continuous cases can be expressed as
a function of the condition number of the metric
M
(
x,t
)
,
which is useful in designing convex optimization-based control
synthesis as shall be seen in Sec. III and IV.
III. M
AIN
R
ESULT
2: CV-STEM C
ONTROL WITH
S
TABILITY AND
O
PTIMIZATION
This section presents the CV-STEM control for general
input-affine nonlinear stochastic systems, incremental stabil-
ity of which is analyzed using contraction theory given in
Theorems 1 and 3. Since the differential dynamics of
δx
used
in contraction analysis can be viewed as an LTV system, we
can use an optimal differential Lyapunov function of the form
δx
>
M
(
x,t
)
δx
without loss of generality [3], thereby finding
M
(
x,t
)
via convex optimization. We note that this is not for
finding an optimal control trajectory and input, which can be
used as a desired trajectory in the present control design.
In Sec. III-E, we present a convex optimization problem
for finding the optimal contraction metric for the CV-STEM
control, which greedily minimizes an upper bound of the
steady-state mean squared tracking error of It
ˆ
o stochastic
nonlinear system trajectories. It is shown that this problem
Authorized licensed use limited to: CALIFORNIA INSTITUTE OF TECHNOLOGY. Downloaded on November 19,2020 at 16:41:38 UTC from IEEE Xplore. Restrictions apply.
0018-9286 (c) 2020 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TAC.2020.3038402, IEEE
Transactions on Automatic Control
6
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. XX, NO. XX, 2021
is equivalent to the original non-convex optimization problem
of minimizing the upper bound.
A. Problem Formulation
Consider the following It
ˆ
o stochastic nonlinear systems with
a control input
u
, perturbed by a
d
-dimensional Wiener process
W
(
t
)
:
dx
=
f
(
x,t
)
dt
+
B
(
x,t
)
udt
+
G
u
(
x,t
)
d
W
dx
d
=
f
(
x
d
,t
)
dt
+
B
(
x
d
,t
)
u
d
dt.
(26)
where
u
:
R
≥
0
→
R
m
,
B
:
R
n
×
R
≥
0
→
R
n
×
m
,
G
u
:
R
n
×
R
≥
0
→
R
n
×
d
, and
x
d
:
R
≥
0
→
R
n
and
u
d
:
R
≥
0
→
R
m
are the desired state and input, respectively. The dynamical
system of the desired state is deterministic as
x
d
and
u
d
are
assumed to be given.
Remark 3:
Since
̇
x
d
−
f
(
x
d
,t
)
∈
Im
B
(
x
d
,t
)
holds for
a feasible desired trajectory,
u
d
can be obtained as
u
d
=
B
(
x
d
,t
)
+
( ̇
x
d
−
f
(
x
d
,t
))
where
(
·
)
+
denotes the Moore-
Penrose inverse. This is the unique least-squares solution
(LSS) to
B
(
x
d
,t
)
u
d
= ̇
x
d
−
f
(
x
d
,t
)
when
Ker
B
(
x
d
,t
) =
{
0
}
and an LSS with the smallest Euclidean norm when
Ker
B
(
x
d
,t
)
6
=
{
0
}
. The desired input
u
d
can also be found
by solving an optimal control problem [1], [44]–[51], [61] and
a general system with
̇
x
=
f
(
x,u
)
can be transformed into an
input-affine form by treating
̇
u
as another input.
In the proceeding discussion, we assume that
f
(
x,t
) = 0
at
x
= 0
and that
f
is a continuously differentiable function.
This allows us to use the following lemma.
Lemma 3:
Let
Ω
be the state set that is a bounded open
subset of some Euclidean space s.t.
0
∈
Ω
⊆
R
n
. Under
the assumptions
f
(0) = 0
and
f
(
x
)
is a continuously dif-
ferentiable function of
x
on
Ω
, there always exists at least
one continuous nonlinear matrix-valued function
A
(
x
)
on
Ω
s.t.
f
(
x
) =
A
(
x
)
x
, where
A
: Ω
→
R
n
×
n
is found by
mathematical factorization and is non-unique when
n >
1
.
Proof:
See [8].
Using Lemma 3, (26) is expressed as
dx
=
A
(
%,x,t
)
xdt
+
B
(
x,t
)
udt
+
G
u
(
x,t
)
d
W
dx
d
=
A
(
%,x
d
,t
)
x
d
dt
+
B
(
x
d
,t
)
u
d
dt
(27)
where
%
= (
%
1
,
···
,%
s
1
)
are the coefficients of the convex
combination of SDC parameterizations
A
i
(
x,t
)
, i.e.,
A
(
%,x,t
) =
s
1
∑
i
=1
%
i
A
i
(
x,t
)
.
(28)
Writing the system dynamics (26) in SDC form provides a
design flexibility to mitigate effects of stochastic noise while
verifying that the system is controllable as shall be seen later.
B. Feedback Control Design
We consider the following feedback control scheme (to be
optimized in Sec. III-E):
u
=
−
K
(
x,t
)(
x
−
x
d
) +
u
d
=
−
R
(
x,t
)
−
1
B
(
x,t
)
>
M
(
x,t
)(
x
−
x
d
) +
u
d
(29)
where
R
(
x,t
)
0
is a weight matrix on the input
u
and
M
(
x,t
)
is a positive definite matrix which satisfies the
following matrix inequality for
γ >
0
:
̇
M
(
x,t
) + 2 sym(
M
(
x,t
)
A
(
%,x,t
)) +
γM
2
(
x,t
)
−
M
(
x,t
)
B
(
x,t
)
R
(
x,t
)
−
1
B
(
x,t
)
>
M
(
x,t
)
0
.
(30)
Define
A
cl
(
%,y,t
)
,
∆
A
(
%,y,t
)
, and
∆
B
(
y,t
)
[7] as
A
cl
(
%,y,t
) =
A
(
%,y
+
x
d
,t
)
−
B
(
y
+
x
d
,t
)
K
(
y
+
x
d
,t
)
∆
A
(
%,y,t
) =
A
(
%,y
+
x
d
,t
)
−
A
(
%,x
d
,t
)
∆
B
(
y,t
) =
B
(
y
+
x
d
,t
)
−
B
(
x
d
,t
)
.
(31)
Substituting (29) into (27) yields
de
=
f
v
(
e,t
)
dt
+
G
u
(
e
+
x
d
,t
)
d
W
(32)
where
e
=
x
−
x
d
and
f
v
(
y,t
) =
A
cl
(
%,e,t
)
y
+ ∆
A
(
%,y,t
)
x
d
+ ∆
B
(
y,t
)
u
d
.
(33)
Lemma 4:
Suppose that the deterministic system is per-
turbed as follows:
̇
x
=
f
(
x,t
) +
B
(
x,t
)(
u
+
d
)
.
(34)
If there exists a positive definite solution
M
(
x,t
)
to the
inequality (30) with
R
(
x,t
) =
S
(
x,t
)
2
0
and
S
(
x,t
)
0
,
then the system with inputs
μ
1
=
S
(
x,t
)
d
,
μ
2
= (
√
2
/γ
)∆
d
and an output
y
= (
√
γ/
2)
M
(
x,t
)(
x
−
x
d
)
, where
∆
d
=
∆
Ax
d
+ ∆
Bu
d
, is finite-gain
L
2
stable and its
L
2
gain is less
than or equal to 1 for each input
μ
1
and
μ
2
.
Proof:
See Appendix I.
C. Incremental Stability Analysis
As we discussed earlier in Sec. II, even when a control input
at
t
=
t
k
is applied during a finite time interval
t
∈
[
t
k
,t
t
+1
)
,
Theorem 3 along with Theorem 2 guarantees that the discrete-
time controller leads to an analogous result to the continuous-
time case (29) if
∆
t
k
is sufficiently small. Thus, we perform
stability analysis for continuous-time dynamical systems. Let
us define a deterministic virtual system of (27) as follows:
̇
y
=
f
v
(
y,t
) =
A
cl
(
%,e,t
)
y
+ ∆
A
(
%,y,t
)
x
d
+ ∆
B
(
y,t
)
u
d
.
(35)
where (35) has
y
=
e
and
y
= 0
as its particular solutions.
The virtual dynamics of (35) is expressed as
δ
̇
y
=
A
cl
(
%,e,t
)
δy
+
φ
(
%,y,t
)
δy
(36)
where
φ
(
%,y,t
) =
∂
(∆
Ax
d
+ ∆
Bu
d
)
/∂y
. Using
f
v
(
y,t
)
,
the virtual system of (32) with respect
y
is defined as
dy
=
f
v
(
y
(
μ,t
)
,t
)
dt
+
G
(
y
(
μ,t
)
,t
)
d
W
(37)
where
μ
∈
[0
,
1]
is introduced to parameterize the trajectories
y
=
e
and
y
= 0
, i.e.,
y
(
μ
= 0
,t
) =
e
,
y
(
μ
= 1
,t
) = 0
,
G
(
y
(0
,t
)
,t
) =
G
u
(
e
+
x
d
,t
)
, and
G
(
y
(1
,t
)
,t
) = 0
n
×
d
. It
can be seen that (37) has
y
=
e
and
y
= 0
as its particular
solutions because we have
•
(37) reduces to (32) when
y
=
e
.
Authorized licensed use limited to: CALIFORNIA INSTITUTE OF TECHNOLOGY. Downloaded on November 19,2020 at 16:41:38 UTC from IEEE Xplore. Restrictions apply.