Stochastic Feedback Control of Systems
with Unknown Nonlinear Dynamics∗

arXiv:1705.09761v1 [cs.SY] 27 May 2017

Dan Yu1 , Mohammadhussein Rafieisakhaei2 and Suman Chakravorty1
Abstract— This paper studies the stochastic optimal control
problem for systems with unknown dynamics. First, an openloop deterministic trajectory optimization problem is solved
without knowing the explicit form of the dynamical system.
Next, a Linear Quadratic Gaussian (LQG) controller is designed
for the nominal trajectory-dependent linearized system, such
that under a small noise assumption, the actual states remain
close to the optimal trajectory. The trajectory-dependent linearized system is identified using input-output experimental
data consisting of the impulse responses of the nominal system.
A computational example is given to illustrate the performance
of the proposed approach.

I. I NTRODUCTION
Stochastic optimal control problems, also known as
Markov Decision Processes (MDPs) have found numerous
applications in the Sciences and Engineering. In general, the
goal is to control a stochastic system so as to minimize
the expected running cost of the system. It is well known
that the global optimal solution for MDPs can be found by
solving the Hamilton-Jacobi-Bellman (HJB) equations [1].
The solution techniques can be further divided into modelbased and model-free techniques according to whether the
method uses an analytical model of the system, or it uses a
black box simulation model or real experiments.
In model-based techniques, many methods rely on a
discretization of the underlying state and action space [2],
and hence, run into the curse of dimensionality, the fact
that the computational complexity grows exponentially with
the dimension of the state space of the problem. The most
computationally efficient approach among these techniques
are trajectory-based methods, first described in [3]. These
methods linearize the nonlinear system equations about a
deterministic nominal trajectory and perform a localized
version of policy iteration to iteratively improve the policy.
For instance, the Differential Dynamic Programming (DDP)
[4], [5] linearizes the dynamics and the cost-to-go function
around a given nominal trajectory, and designs a local
feedback controller using DP. The iterative LQG (iLQG)
*This material is based upon work partially supported by NSF under
Contract Nos. CNS-1646449 and Science & Technology Center Grant CCF0939370, the U.S. Army Research Office under Contract No. W911NF15-1-0279, and NPRP grant NPRP 8-1531-2-651 from the Qatar National
Research Fund, a member of Qatar Foundation, AFOSR contract Dynamic
Data Driven Application Systems (DDDAS) contract FA9550-17-1-0068
and NSF NRI project ECCS-1637889.
1 D. Yu and S. Chakravorty are with the Department of Aerospace
Engineering, and 2 M. Rafieisakhaei Electrical and Computer
Engineering, Texas A&M University, College Station, Texas,
77840 USA. {yudan198811@hotmail.com, mrafieis,

schakrav@tamu.edu}

[6], [7], which is closely related to DDP, considers the first
order expansion of the dynamics (in DDP, a second order
expansion is considered), and designs the feedback controller
using Riccati-like equations, which is shown to be computationally more efficient. In both approaches, the control policy
is executed to compute a new nominal trajectory, and the
procedure is repeated until convergence.
Alternatively, the Trajectory-optimized LQG (T-LQG) approach [8]–[10] was recently proposed by us, which shows
that under a first order approximation of the dynamics and
cost-to-go function, a near optimal solution can be found by
first solving a deterministic trajectory optimization problem,
followed by a linear time-varying closed-loop controller
design problem. This separated approach can also be extended to the model-free situation, which is the subject
of this paper: we use a gradient descent algorithm, and a
Linear Time-Varying (LTV) system identification technique,
in conjunction with a black box simulation model of the
process in order to accomplish the “separated” design.
In the model-free case, the most popular approaches in the
community are the Adaptive Dynamic Programming (ADP)
[11], [12] and Reinforcement Learning (RL) paradigms [13],
[14]. They are essentially the same in spirit and seek to
improve the control policy for a given black box system by
repeated interactions with the environment while observing
the system’s responses. The repeated interactions, or learning
trials, allow these algorithms to construct a solution to the
DP equation, in terms of the cost-to-go function, in an online
and recursive fashion. Another variant of the RL techniques
is the so-called Q-learning method, where the basic idea
is to estimate a real-valued Q(x, a) function of states, x,
and actions, a, instead of the cost-to-go function, V (x).
For continuous state and control space problems, the costto-go functions and the Q-functions are usually represented
in a functionally parameterized form; for instance, in the
linearly parametrized form Q(x, a) = θT φ(x, a), where θ is
the unknown parameter vector, and φ is a pre-defined basis
function.
Multi-layer neural networks may also be used as nonlinearly parameterized approximators instead of the linear
architecture above. The ultimate goal of these techniques is
the estimation/learning of the parameters θ from learning
trials/repeated simulations of the underlying system. However, the size of the parameter θ grows exponentially in the
size of the state space of the problem without a compact
parametrization of the cost-to-go or Q function in terms of
the a priori chosen basis functions for the approximation, and
hence, these techniques are typically subject to the curse of

dimensionality. Albeit a compact parametrization may exist,
a priori, it is usually never known.
In the past several years, techniques based on the
DDP/iLQG approach [4], [5], [15], [16], such as the RL
techniques [17], [18] have shown the potential for RL
algorithms to scale to higher dimensional continuous state
and control space problems, in particular, high dimensional
robotic task planning and learning problems. These methods
are a localized version of the policy gradient [19]–[21]
technique that seek to directly estimate the feedback policy
via a compact parameterization. For continuous state and
control space problems, the method of choice is to wrap
an LQR feedback policy around a nominal trajectory and
then perform a recursive optimization of the feedback law,
along with the underlying trajectory, via repeated simulations/iterations. However, the parametrization can still be
very large for partially observed problems (O(d2 )) where
d is the dimension of the state space) or large motion
planning problems such as systems governed by partial
differential equations wherein the (discretized) state is very
high dimensional (millions of states). Furthermore, there are
convergence problems with these techniques that can lead to
policy chatter [17].
Fundamentally, rather than solve the derived “Dynamic
Programming” problem as in the majority of the approaches
above that requires the simultaneous optimization of the
feedback law and the underlying trajectory, our approach is
to directly solve the original stochastic optimization problem
in the “separated open loop/closed loop” fashion wherein:
1) we solve an open loop deterministic optimization problem
to obtain an optimal nominal trajectory in a model-free
fashion, and then 2) we design a closed loop controller
for the linearized time-varying system around the nominal
trajectory, again in a model-free fashion. This “divide and
conquer” strategy is nonetheless theoretically sound as shown
in the companion paper [10]. The primary contributions of
the proposed approach are:
1) Compared to other RL and ADP techniques, implementation using the proposed approach is simple. The stochastic
optimal control problem is separated into two sub-problems:
deterministic open-loop trajectory optimization problem and
a linear time-varying system identification problem, and in
each sub-problem, standard approaches can be used. The
open-loop optimization problem is solved using gradient
descent and input perturbations. The linearized system is
identified via the time-varying ERA [22], using the impulse
responses of the nominal system, and an LQG controller is
designed for the resulting linearized system. All of the above
is accomplished by only considering a sequence of open loop
impulse responses of the unknown system.
2) Unlike other ADP and RL techniques, we specify
a detailed set of experiments to accomplish the closed
loop controller design for the unknown nonlinear system.
This series of experiments consists of a sequence of input
perturbations to collect the impulse response of the system,
first to find an optimized nominal trajectory, and then to
recover the LTV system corresponding to the perturbations

of the nominal system in order to design the LQG controller
corresponding to the LTV system.
3) Albeit not covered in this paper due to a lack of space,
in general, for large scale systems with partially observed
states, the time-varying ERA constructs a reduced order
model (ROM) of the LTV system, and hence, results in
a reduced order estimator and controller. For example, the
computational complexity of designing the LQG controller
using the proposed approach is O(n3r ), while it is O(n3x ) for
the original system, where nr is the order of the reduced
model, and nr  nx . Therefore, for large scale systems
such as partially observed systems and systems governed by
PDEs, the online implementation of the LQG policy using
the proposed approach is still computationally tractable.
The rest of the paper is organized as follows. Section
II formulates the problem. In Section III, we propose a
separation-based stochastic optimal control algorithm with
discussions of implementation. Last, in Section IV, we test
the proposed approach using the inverted pendulum problem.
II. P ROBLEM F ORMULATION
Consider a discrete time nonlinear dynamical system:
xk+1 = f (xk , uk , wk ),
yk = h(xk , vk ),
nx

ny

(1)
nu

are the state,
where xk ∈ R , yk ∈ R , uk ∈ R
measurement, and the control vectors at time k, respectively,
process model f (·) and measurement model h(·) are nonlinear, the process noise wk and measurement noise vk are
assumed to be zero-mean, uncorrelated Gaussian white noise
with covariances W and V , respectively.
Assumption 1: We assume the system is fully observed:
yk = h(xk , vk ) = xk .
(2)
Remark 1: We make Assumption 1 to simplify the treatment of the problem; since otherwise, the stochastic control
problem needs to be treated as a partially observed MDP
(POMDP). This generalization may be done in a reasonably
straightforward fashion (later discussed in Remark 7).
Stochastic Control Problem: For the system with unknown nonlinear dynamics, f (·) , the optimal control problem is to find the control policies π = {π0 , π1 , · · · , πN −1 }
in a finite time horizon [0, N ], where πk is the control policy
at time k, i.e., uk = πk (xk ), to minimize the cost function
N
−1
X
Js = E(
(xTk Qk xk + uTk Rk uk ) + xTN QN xN ),
(3)
k=0

where Qk , QN  0 and Rk  0.
III. S TOCHASTIC F EEDBACK C ONTROL A LGORITHM
We compute a locally optimal solution to the stochastic
control problem in a separated open loop/closed loop (SOC)
fashion, i.e., we first solve a noiseless open-loop optimization
problem to find a nominal optimal trajectory. Next, we design
a linearized closed-loop controller around the nominal trajectory, such that, with existence of stochastic perturbations,
the state stays close to the optimal open-loop trajectory. The
separation-based approach has always been used by Control

Engineers in Aerospace Guidance and Robotics problems in
a heuristic fashion [3]. However, our recent companion work
[10], using the theory of Large Deviations, shows that this
separation, results in a near-optimal policy in the small noise
case. Moreover, experimental results confirm its validity for
moderate noise levels.
The open loop optimization problem could be solved using
a general nonlinear programming solver without knowing
the explicit form of the underlying dynamics, i.e., it only
accesses a black box simulation model of the dynamics.
Next, we perform small input perturbations about the nominal trajectory in order to obtain the impulse responses
of the LTV system governing the perturbations form the
nominal trajectory, and identify the resulting linear timevarying derivation system from these impulse responses using time-varying Eigensystem Realization Algorithm (ERA)
[22]. We consider quadratic cost functions, and design an
LQR controller which results in an optimal linear control
policy around the nominal trajectory. We discuss each of the
above steps in the following section.
A. Open Loop Optimization
Consider the noiseless nonlinear system:
xk+1 = f (xk , uk , 0), yk = xk ,

(4)

with known initial state x0 , and let the N -step cost function
N
−1
X
−1
Jd (x0 , {uk }N
)
=
(xTk Qk xk + uTk Rk uk ) + xTN QN xN .
k=0

Algorithm 1: Gradient descent Algorithm
(0)

1
2
3
4

5

6
7
8

N −1
Input : Design Parameters U (0) = {uk }k=0
, α, h, 
N −1
Output: Optimal control sequence {ūk }k=0
n ← 0, ∇U Jd (x0 , U (0) ) ← ;
while ∇U Jd (x0 , U (n) ) ≥  do
Evaluate the cost function with control U (n) ;
(n)
Perturb each control variable ui by h,
i = 0, · · · , N − 1, and calculate the gradient
vector ∇U Jd (x0 , U (n) ) using (7) and (8);
Update the control policy:
U (n+1) ← U (n) − α∇U Jd (x0 , U (n) );
n ← n + 1;
end
N −1
{ūk }k=0
← U (n) .

sequence of input-output tests. Higher order approaches
other than gradient descent are possible [24], however, for
a general system with complex cost functions, the gradient
descent approach is easy to implement and is amenable to
very large scale parallelization.
B. Linear Time-varying System Identification
We linearize the system (1) around the optimal nominal
control and its corresponding state trajectory denoted by
−1
N
{ūk }N
k=0 and {x̄k }k=0 , respectively, as:
δxk+1 = Ak δxk + Bk δuk + Gk wk ,
δyk = Ck δxk + Fk vk ,

k=0

The open loop optimization problem is to find the control
−1
sequence {ūk }N
k=0 , such that for a given initial state x0 ,
−1
N −1
{ūk }N
k=0 = arg min Jd (x0 , {uk }k=0 ),
N −1
{uk }k=0

s. t. xk+1 = f (xk , uk , 0),
yk = xk .

(5)

The problem is solved using the gradient descent approach
[23], [24], and the procedure is illustrated as follows.
Starting from an initial guess of the control sequence
(0)
U (0)={uk }N−1
k=0 , the control policy is updated iteratively via
U (n+1) = U (n) − α∇U Jd (x0 , U (n) ),

(6)
(n)

until a convergence criterion is met, where U
=
(n)
−1
th
{uk }N
denotes
the
control
sequence
in
the
n
iteration,
k=0
and α is the step size parameter. The gradient vector is
defined as:


∂Jd
∂Jd ∂Jd
,
,
·
·
·
,
|x ,{u(n) }N −1 , (7)
∇U Jd (x0 , U (n) )= ∂u
∂u
∂u
0
1
N −1
0

k

k=0

and without knowing the explicit form of the cost function,
each partial derivative with respect to the ith control variable
ui is calculated as:
∂Jd
1
(n)
(n)
(n)
(n)
|x0 ,U (n) = (Jd (x0 , u0 , u1 ,· · ·, ui +h,· · ·, uN −1 )
∂ui
h
(n)
(n)
(n)
(n)
−Jd (x0 , u0 , u1 ,· · ·, ui ,· · ·, uN −1 )), (8)
where h is a small constant perturbation.
Algorithm 1 summarizes the gradient descent approach.
Remark 2: The open loop optimization problem is solved
using a black box model of the underlying dynamics, with

(9)

where δxk = xk − x̄k describes the state deviations from
the nominal trajectory, δuk = uk − ūk describes the control
deviations, δyk = yk − h(x̄k , 0) describes the measurement
deviations, and
∂f (x, u, w)
∂f (x, u, w)
|x̄k ,ūk ,0 , Bk=
|x̄k ,ūk ,0 ,
Ak =
∂x
∂u
∂f (x,u,w)
∂h(x,v)
∂h(x,v)
Gk=
|x̄k ,ūk ,0 ,Ck=
|x̄k ,0 ,Fk=
|x̄k ,0 .
∂w
∂x
∂v
Given the LTV system {Ak , Bk , Ck , Gk , Fk } between
time [0, N − 1], an LQR controller to track the nominal
trajectory could be designed. However, since the dynamics
are unknown, we first need to identify the LTV system.
The time-varying ERA is used to construct a state space realization (Âk , B̂k , Ĉk , Ĝk , F̂k } of system (9). The state space
realization is constructed using input and output experimental
data and is shown to be minimal and balanced. The details
of the time-varying ERA can be found in [22] and is briefly
summarized next.
Define the generalized Markov parameters hk,j as:


Ck Ak−1 Ak−2 · · · Aj+1 Bj , if j < k − 1,
hk,j = Ck Bk−1 ,
if j = k − 1, (10)


0,
if j > k − 1,
and the generalized Hankel matrix as:


hk,k−1
hk,k−2
···
hk,k−q
h
hk+1,k−2 · · · hk+1,k−q 

(p,q)  k+1,k−1
Hk =
,(11)
..
..
..

| {z } 
.
.
···
.
pny ×qnu

hk+p−1,k−1 hk+p−1,k−2 · · · hk+p−1,k−q

with design parameters p and q.
The time-varying ERA starts with an estimation of the
generalized Markov parameters from input-output data using
least squares solution as follows.
Consider system (9) with zero noise and for simplicity,
assume δx0 = 0. Run M simulations and in the ith simulation, choose input sequence {δuit }kt=0 , and collect the output
δyki . The superscript (·)i denotes the experiment number.
From the input-output map, the generalized Markov parameters {hk,j }kj=0 could be recovered via solving the least
squares problem:


δyk1 δyk2 · · · δykM = 0 hk,k−1 hk,k−2 · · · hk,0

 1
δuk
δu2k · · · δuM
k

δu1k−1 δu2k−1 · · · δuM
k−1

,(12)
× .
.
.
..
.. 

 ..
δu10
δu20 · · · δuM
0
where M is a design parameter and is chosen such that the
least squares solution is possible.
After recovering the generalized Markov parameters, two
(p,q)
(p,q)
Hankel matrices Hk
and Hk+1 are constructed using
(11), and here, the design parameters p and q are chosen
such that min{pny , qnu } ≥ nx , and could be tuned for best
performance.
Then we solve the singular value decomposition problem:
(p,q)

Hk

1/2

1/2

= Uk Σk Σk VkT .
| {z } | {z }
(p)

Ok

(q)

Rk−1
(p,q)

(p)↓

(p)↑

,

nr ×nr

reduced significantly. In comparison, the computational complexity of designing an LQG controller using the identified
system is O(( nnxr )3 ) using the full order system.
Note that we cannot perturb the system (9) directly.
Instead, we identify the generalized Markov parameters as
follows.
Run M parallel simulations with the noise-free system:
xik+1 = f (xik , ūk + δuik , 0),

(p)

(14)

ny ×nr

(p)↓

δyki = yki − h(x̄k , 0).

(16)

δak+1 = Âk δak + B̂k δuk + Ĝk wk ,
δyk = Ĉk δak + F̂k vk ,

nr ×nu

Ĉk = Ok (1 : ny , :),
|{z}

(15)

where (ūk , x̄k ) is the open loop optimal trajectory. Then
solve the same least squares problem with (12).
For simplicity, we assume that the process noise is independent of the state and control variables, and Gk = Inx ×nx ,
Fk = Iny ×ny , while the proposed algorithm is extendable to
identify Gk and Fk .
In general, the identified deviation system is:

(q)

B̂k = Rk (:, 1 : nu ),
|{z}

yki = h(xik , 0),
and therefore,

(13)

Suppose the rank of the Hankel matrix Hk
is nr , where
nr ≤ nx . Then Σk ∈ Rnr ×nr is the collection of all nonzero singular values, and Uk ∈ Rpny ×nr , Vk ∈ Rqnu ×nr are
the corresponding left and right singular vectors. Similarly,
(p,q)
(p)
(q)
Hk+1 = Ok+1 Rk .
Thus, the identified system using time-varying ERA is:
Âk = (Ok+1 )+ Ok
|{z}

Algorithm 2: LTV System Identification
Input : Design Parameters M, p, q
Output: {Âk , B̂k , Ĉk }
1 k ← 0;
2 while k ≤ N − 1 do
3
Identify generalized Markov parameters with input
and output experimental data using (15), (16), (12);
(p,q)
4
Construct the generalized Hankel matrices Hk ,
(p,q)
Hk+1 using (11);
5
Solve the SVD problem, and construct
{Âk , B̂k , Ĉk } using (14);
6
k ← k + 1;
7 end

where (·)+ denotes the pseudo inverse of (·), Ok+1 contains
(p)
(p)↑
the first (p − 1)ny rows of Ok+1 , and Ok
contains the
(p)
last (p − 1)ny rows of Ok . Here, we assume that nr is
constant through the time period of interest, which could
also be relaxed.
Remark 3: The uncontrollable or unobservable eigenmodes of the dynamical system are not present in the
input-output map, and hence, the state space realization
using time-varying ERA is balanced in the sense that only
the controllable and observable eigenmodes are preserved.
Hence, for systems with high dimensions, such as systems
discretized from partial differential equations (PDEs), and
with partially observed states, we have nr  nx . Therefore,
one major contribution of this work is that we design a
reduced order estimator and controller using the identified
system, which implies that the computational complexity is

(17)

where δak ∈ <nr denotes the reduced order deviation states.
Algorithm 2 summarizes the time-varying ERA.
C. Closed-loop Controller Design
Given the identified deviation system (17), we design the
closed-loop controller to track the optimal nominal trajectory,
which is to minimize the cost function
N
−1
X
Jf =
(δâTk Qk δâk + δuTk Rk δuk ) + δâTN QN δâN , (18)
k=0

where δâk denotes the estimates of the deviation state δak .
For the linear system (17), the separation principle of control
theory is used to separate the design of an estimator and a
fully observed controller.
The feedback controller is:
δuk = −Lk δâk ,
where δâk denotes the estimates from a Kalman observer,
and the feedback gain Lk is computed by solving two
decoupled Riccati equations:
Lk = (B̂kT Sk+1 B̂k + Rk )−1 B̂kT Sk+1 Âk ,

(19)

Algorithm 3: Stochastic Feedback Control
1

2
3

4
5
6
7
8
9
10

Solve the deterministic open-loop optimization problem,
−1
N
to obtain the optimal trajectory {ūk }N
k=0 , {x̄k }k=0 ;
Identify the LTV system using Algorithm 2;
Solve the decoupled Riccati equations (19, 20) for
feedback gain {Lk }N
k=0 ;
k ← 0, δâ0 ← 0, given P0 ;
while k ≤ N − 1 do
uk ← ūk − Lk δâk , xk+1 ← f (xk , uk , wk ),
yk+1 ← h(xk+1 , vk+1 );
Update δâk , Pk using (21, 22);
k ← k + 1.;
end

where Sk is determined by running the following Riccati
equation backward in time:
Sk=ÂTk Sk+1 Âk + Qk
− Âk Sk+1 B̂k (B̂kT Sk+1 B̂k + Rk )−1 B̂kT Sk+1 Âk , (20)
with terminal condition SN = QN .
The Kalman filter observer is designed as follows:
δâk+1 =Âk δâk + B̂k δuk
+ Kk+1 (δyk+1 − Ĉk+1 (Âk δâk + B̂k δuk )), (21)
with δyk = h(xk , vk ) − h(x̄k , 0), and the covariance of the
estimation is:
Pk+1 =Âk (Pk − Pk ĈkT (Ĉk Pk ĈkT + F̂k V F̂kT )−1 Ĉk Pk )ÂTk
+ Ĝk W ĜTk ,

(22)

where the Kalman gain is:
Kk = Pk ĈkT (Ĉk Pk ĈkT + F̂k V F̂kT )−1 .

(23)

Remark 4: Albeit we have only considered the fully observed problem in this paper, any implementation will have
noisy measurements, and thus, an observer will be required to
implement a feedback controller. In our case, it is simply the
LQG controller as outlined above which can be conveniently
designed using the identified LTV system .
D. Stochastic Feedback Control Algorithm
Algorithm 3 summarizes the Stochastic Control Algorithm.
Remark 5: Extension to Data-Driven Controller Design. As mentioned in data-based LQG [25], [26] and
data-driven MPC control [27], [28], the linear system
(Âk , B̂k , Ĉk , Êk , F̂k ) need not be identified to design the
LQG controller. Once the generalized Markov parameters
are recovered, then together with other input-output data
matrices, the controller can also be directly designed. The
extension to designing such a “direct” data-driven controller
is not covered in this paper.
Remark 6: Discussion of Implementation Issues.
•

The proposed approach is implemented in an offineonline fashion. Starting from an given initial state,
the open-loop trajectory optimization problem and the
LTV system identification problem are solved offline,

using the input-output experimental data. Then the LQG
policy is implemented online.
• An initial estimate of the optimal trajectory is required
in the proposed approach, which poses the challenge,
and results in a sub-optimal path.
• The proposed approach is valid under small noise
assumption. In practice, with the presence of noise,
non-linearities, and unknown perturbations, the actual
state deviates from the nominal trajectory during the
execution, and if the actual system deviates too much
from the nominal trajectory, the linearization becomes
invalid. Therefore, once the deviation is greater than
some predefined threshold, a replanning starting from
the current location can be performed.
Remark 7: Extension to Partially Observed MDP
(POMDP). The proposed approach can also be extended to
solve the stochastic optimal control problem for a general
nonlinear system with unknown dynamics f (·) and h(·).
The only difference is that for the partially observed case,
we need to solve the open loop optimization problem in
belief space instead of state space. Assume that the belief
bk = (µk , Σk ) is approximately Gaussian, where µk and Σk
represent the mean and covariance of the belief state. The
challenge is that we do not have access to the covariance
evolution equation. An effective way is to simulate the
covariance evolution using an Ensemble Kalman Filter
[29]. Then, the gradient descent approach can be used for
solving the open loop belief space optimization problem.
We linearize the system around the optimal control sequence
{ūk }, with the associated nominal belief state (µ̄k , Σ̄k ),
and follow the same system identification and closed loop
controller design procedure. The extension to POMDPs will
be considered in future work.
IV. C OMPUTATIONAL R ESULTS
We test the method using the inverted pendulum problem
[30]. The dynamics of the inverted pendulum mounted on a
motor driven cart are given as:
ẋ1 = x2 ,
u cos x1 − (M + m)g sin x1 + ml(cos x1 sin x1 )x22
,
ẋ2 =
ml cos2 x1 − (M + m)l
ẋ3 = x4 ,
u + ml(sin x1 )x22 − mg cos x1 sin x1
(24)
ẋ4 =
M + m − m cos2 x1
with state variables
x1 = θ, x2 = θ̇, x3 = x, x4 = ẋ,

(25)

where θ is the tilt angle referenced to the vertically upward
direction and x represents the cart position, u is the xdirected external force, the cart mass is M = 2.4kg, the
pendulum point mass is m = 0.23kg, the rod length is
l = 0.36m, and the standard gravity is g = 9.81m/s2 .
We use a time step of 0.1s to discretize the system in
time:
xk+1 = f (xk , uk ) + wk ,
yk = xk + vk ,

(26)

2

1st row

Initial Guess
Optimal

1.5

2nd row

3

20

2
10

1

0.5

0

h k, j

h k, j

Control u

1

-1
-2
-3

0

0
-10
-20

-4

-0.5

-5

0

10

20

3
-1.5

-30

40

0

10

20

1

2

3

4

5

Time (s)

Fig. 1. Open Loop Optimization Problem.

where f is the unknown dynamics, the process noise and
measurement noise wk , vk are Gaussian white noise with
zero mean and covariances W = 0.01I4×4 and V =
0.01I4×4 , respectively.
The total simulation time is 5 seconds. The control ob−1
jective is to find the control sequence {uk }N
k=0 to make the
pendulum swing up within 3.5 seconds and minimize the
cost function Js in (3), and maintain the pendulum to the
upright vertical position. The initial state x0 = [π; 0; 0; 0] is
known, and the final state is xN = [0; 0; xf ; 0], where the
cart position xf is not restricted. The open loop optimization
problem is solved using Matlab nonlinear optimization solver
fminunc. The initial guess and optimal control are shown
in Fig. 1. The corresponding nominal trajectory {xk }N
k=0 are
plotted in Fig. 3.
The implementation of Algorithm 2 to identify the linearized system is performed as follows.
(p,q)
The size of the generalized Hankel matrix Hk
is pny ×
qnu , and as discussed before, the design parameters p and q
should be chosen such that min{pny , qnu } ≥ nx , which for
the current problem, nx = 4, nu = 1, and ny = 4. Design
parameters p, q are selected by trial and error. We start with
some initial guess of p, q, compare the impulse responses
of the original system and identified system, and check if
the accuracy of the identified system is acceptable. Here, we
choose p = q = 5.
We run M parallel simulations to estimate the generalized
Markov parameters {hk,j }kj=0 , k = 1, 2, · · · , N − 1. We
−1
perturb the nominal control {ūk }N
k=0 with impulse, i.e.,
i N −1
{δuk }k=0 = (0, 0, · · · , 0.01, · · · , 0) is the input perturbation sequence in the ith simulation, where only the ith
element is nonzero. Therefore, we choose design parameter M = N . In each simulation, we collect the outputs
−1
{δyki }N
k=0 in (16) corresponding to the control input {ūk +
i N −1
δuk }k=0 , and solve the least squares problem using (12).
We construct the generalized Hankel matrix, and solve
the singular value decomposition problem. The rank of the
Hankel matrix nr = 4, and hence, the identified system
Âk ∈ <4×4 . To test the accuracy of the identified system, we
calculate the identified Markov parameters using Âk , B̂k , Ĉk ,
and compare with the actual generalized Markov parameters
(calculated using impulse responses as mentioned above).
Since the generalized Markov parameter hk,j ∈ <4×1 , all
four elements are shown in Fig. 2 for k = 33.

rd

row

4

0.3

1

0.25

0.8

40

th

30

40

row

0.6

0.2
0.15
0.1

0.4
0.2
0

0.05
0

30

j

h k, j

0

h k, j

-2

30

j

-1

-0.2
0

10

20

30

40

-0.4

0

j

10

20

j

Fig. 2. Comparison of the generalized Markov parameters {hk,j }kj=0 as
a function of time index j. The solid blue lines represent the identified
Markov parameters, and the red asterisk markers represent the actual Markov
parameters.

With the identified linearized system, we design the
closed-loop controller. We run 1000 individual simulations,
where performance of the closed-loop controller is shown
in Fig. 3. There are two observations: 1) the averaged state
estimates over 1000 Monte-Carlo simulations runs (plotted
in solid green lines) are close to the open-loop optimal trajectory (plotted in red), implying that the control objective to
minimize the expected cost function could be achieved using
the proposed approach; 2) in this problem, the closed-loop
controller has better control in the angle and angular velocity,
so that the corresponding 2σ bound is tight. However, the
uncertainties of position and velocity rise especially after
the pendulum reaches the upright position. This is due to
the fact that with the presence of process noise in all states,
it is not possible to stabilize all the states simultaneously.
V. C ONCLUSION
In this paper, we have proposed a separation-based design
of the stochastic optimal control problem for systems with
unknown nonlinear dynamics and fully observed states in
a separated open loop-closed loop fashion. First, we design a deterministic open-loop optimal trajectory. Then we
identify the nominal linearized system using time-varying
ERA. The open-loop optimization and system identification
are implemented offline, using the impulse responses of
the system, and an LQG controller based on the ROM
is implemented online. The offline learning procedure is
simple, and the online implementation is fast. We tested the
proposed approach on the inverted pendulum problem, and
showed the performance of the proposed approach. Future
work will generalize the proposed approach to large-scale,
partially observed systems.
R EFERENCES
[1] D. P. Bertsekas, Dynamic Programming and Optimal Control, Two
Volume Set, 2nd ed. Athena Scientific, 1995.
[2] M. Falcone, “Recent Results in the Approximation of Nonlinear Optimal Control Problems,” in Large-Scale Scientific Computing LSSC,
2013.

4.5

5

Angular Velocity θ̇ (rad/s)

4
3.5

Angle θ (rad)

3
2.5
2
1.5
1
0.5

0

-5

0
-0.5

0

1

2

3

4

-10

5

0

1

2

3

0.4

0.4

4

0.3

3

0.3

0.2

0.2
0.15
0.1

0
-0.1
-0.2

0.05

-0.3

0

-0.4

-0.05

-0.5

0

1

2

3

Time (s)

4

5

5

2

0.1

Control u

0.25

4

Time (s)

0.35

Velocity ẋ (m/s)

position x(m)

Time (s)

1
0
-1
-2
-3

0

1

2

3

Time (s)

4

5

-4

0

1

2

3

4

Time (s)

Fig. 3. Performance of the closed-loop controller. The solid green lines represent the mean estimate, the green shaded colors indicate the uncertainty with
2σ bound, the red line with asterisk markers represent the open-loop nominal trajectory, and the blue line represents the actual states in one realization.
[3] A. Bryson and H. Y.-C., Applied Optimal Control: Optimization,
Estimation and Control. Washington: Hemisphere Pub. Corp., 1975.
[4] D. Jacobsen and D. Mayne, Differential Dynamic Programming.
Elsevier, 1970.
[5] E. Theoddorou, Y. Tassa, and E. Todorov, “Stochastic Differential
Dynamic Programming,” in Proceedings of American Control Conference, 2010.
[6] E. Todorov and W. Li, “A generalized iterative LQG method for
locally-optimal feedback control of constrained nonlinear stochastic
systems,” in Proceedings of American Control Conference, 2005, pp.
300 – 306.
[7] W. Li and E. Todorov, “Iterative linearization methods for approximately optimal control and estimation of non-linear stochastic system,”
International Journal of Control, vol. 80, no. 9, pp. 1439–1453, 2007.
[8] M. Rafieisakhaei, S. Chakravorty, and P. R. Kumar, “T-lqg: Closedloop belief space planning via trajectory-optimized lqg,” in 2017 IEEE
International Conference on Robotics and Automation (ICRA). IEEE,
2017.
[9] M. Rafieisakhaei, S. Chakravorty, and P. R. Kumar, “Belief space planning simplified: Trajectory-optimized lqg (t-lqg) (extended report),”
arXiv preprint arXiv:1608.03013, 2016.
[10] M. Rafieisakhaei, S. Chakravorty, and P. R. Kumar, “A Near-Optimal
Separation Principle for Nonlinear Stochastic Systems Arising in
Robotic Path Planning and Control,” in 56th IEEE Conference on
Decision and Control (CDC), submitted, 2017.
[11] R. P. Bithmead, V. Wertz, and M. Gerers, Adaptive Optimal Control:
The Thinking Man’s G.P.C. Prentice Hall Professional Technical
Reference, 1991.
[12] X. Zhong, H. He, H. Zhang, and Z. Wang, “Optimal Control for
Unknown Diiscrete-Time Nonlinear Markov Jump Systems Using
Adaptive Dynamic Programming,” IEEE Transactions on Neural networks and learning systems, vol. 25, no. 12, pp. 2141–2155, 2014.
[13] S. G. Khan et al., “Reinforcement learning and optimal adaptive
control: An overview and implementation examples,” Annual Reviews
in Control, vol. 36, pp. 42–59, 2012.
[14] D. Mitrovic, S. Klanke, and S. Vijayakumar, Adaptive Optimal Feedback Control with Learned Internal Dynamics Models.
Berlin,
Heidelberg: Springer Berlin Heidelberg, 2010, pp. 65–84.
[15] S. Levine and P. Abbeel, “Learning Neural Network Policies with
Guided Search under Unknown Dynamics,” in Advances in Neural
Information Processing Systems, 2014.
[16] S. Levine and K. Vladlen, “Learning Complex Neural Network Policies with Trajectory Optimization,” in Proceedings of the International
Conference on Machine Learning, 2014.

[17] R. Akrour, A. Abdolmaleki, H. Abdulsamad, and G. Neumann,
“Model Free Trajectory Optimization for Reinforcement Learning,”
in Proceedings of the International Conference on Machine Learning,
2016.
[18] E. Todorov and Y. Tassa, “Iterative Local Dynamic Programming,” in
Proc. of the IEEE Int. Symposium on ADP and RL., 2009.
[19] J. Baxter and P. Bartlett, “Infinite Horizon Policy-Gradient Estimation,” Journal of Artificial Intelligence Research, vol. 15, pp. 319–350,
2001.
[20] R. S. Sutton, D. Mcallester, S. Singh, and Y. Mansour, “Policy
Gradient Methods for Reinforcement Learning with Function Approximation,” in Proc. 1999 Neural Information Proc. Sys., 1999.
[21] P. Marbach, Simulation based Optimization of Markov Reward Processes, PhD Thesis. Boston, MA: Massachusetts Institute of Technology, 1999.
[22] M. Majji, J.-N. Juang, and J. L. Junkins, “Time-varying Eigensystem
Realization Algorithm,” Journal of Guidance, Control, and Dynamics,
vol. 33, no. 1, pp. 13–28, 2010.
[23] A.E.Bryson and W. Denham, “A steepest-ascent method for solving
optimum programming problems,” Journal of Applied Mechanics,
vol. 29, no. 2, 1962.
[24] A. Gosavi, Simulation-based optimization: Parametric optimization
techniques and reinforcement learning. Norwell, MA, USA: Kluwer
Academic Publishers, 2003.
[25] G. Shi and R. E. Skelton, “Markov Data-Based LQG Control,” Journal
of Dynamic Systems, Measurement, and Control, vol. 122, pp. 551–
559, 2000.
[26] W. Favoreel et al., “Closed-loop model-free subspace-based lqgdesign,” in Proccedings of the 7th Mediterranean Conference on
Control and Automation, 1999, pp. 1926–1939.
[27] Z.-S. Hou and Z. Wang, “From model-based control to data-driven
control: Survey, classification and perspective,” Information Sciences,
vol. 235, no. 3-35, 2013.
[28] X. Wu, J. Shen, Y. Li, and K. Y.Lee, “Data-Driven Modeling and
Predictive Control for Bioler-Turbine Unit,” IEEE Transactions on
energy conversion, vol. 228, no. 3, pp. 470–481, 2013.
[29] G. Evensen, “Advanced Data Assimilation for Strongly Nonlinear
Dynamics,” Mon. weath. Rev., vol. 125, pp. 1342–1354, 1997.
[30] L. B. Prasad, B. Tyagi, and H. O. Gupta, “Optimal Control of
Nonlinear Inverted Pendulum System using PID Controller and LQR:
Performance Analysis Without and With Disturbance Input,” International Journal of Automation and Computing, vol. 11, no. 6, pp.
661–670, 2014.

