Differential Dynamic Programming with Temporally Decomposed
Dynamics
Akihiko Yamaguchi1 and Christopher G. Atkeson1
Abstract— We explore a temporal decomposition of dynamics
in order to enhance policy learning with unknown dynamics.
There are model-free methods and model-based methods forpolicy learning with unknown dynamics, but both approaches
have problems: in general, model-free methods have less generalization ability, while model-based methods are often limitedby the assumed model structure or need to gather many samples
to make models. We consider a temporal decomposition of
dynamics to make learning models easier. To obtain a policy,
we apply differential dynamic programming (DDP). A featureof our method is that we consider decomposed dynamics evenwhen there is no action to be taken, which allows us to decompose dynamics more flexibly. Consequently learned dynamics
become more accurate. Our DDP is a first-order gradientdescent algorithm with a stochastic evaluation function. In DDP
with learned models, typically there are many local maxima. In
order to avoid them, we consider multiple criteria evaluation
functions. In addition to the stochastic evaluation function, we
use a reference value function. This method was verified withpouring simulation experiments where we created complicateddynamics. The results show that we can optimize actions with
DDP while learning dynamics models.

I. I NTRODUCTION
Manipulation of non-rigid materials is an important challenge to make robots useful in the home. Pouring a range ofmaterials, such as liquids, and granular particles, is a good
example. Folding towels is another example [1]. A difficultyis that the behavior of such material is too complicated tomake a precise model. Leaning methods are used in thosecases, such as reinforcement learning (e.g. [2]). Kormushev
et al. applied a reinforcement learning method to obtain apolicy of flipping a pancake in a frying pan [3].
Yamaguchi et al. created a pouring behavior that cangeneralize well [4]. They used simple learning methods
combined with planning methods to achieve generalization.
They did a careful primitive skill design which was further
improved by learning from practice methods. Pouring is stilla difficult task. An example is shown in Fig. 1 where a robot
is pouring granular material (dried peas) by shaking. Somepeas spilled out of the container. The problem is due to thecomplicated dynamics of flow. Designing behaviors is very
difficult even for professionals.
Therefore we consider a policy learning method for unknown dynamics. There are model-free methods (e.g. [2])and model-based methods (e.g. [5]) for policy learning withunknown dynamics, but both approaches have problems: ingeneral, model-free methods have less generalization ability,
1 A. Yamaguchi and C. G. Atkeson are with the Robotics Institute,Carnegie Mellon University, 5000 Forbes Avenue, Pittsburgh PA 15213,
United States info@akihikoy.net

Fig. 1. PR2 pouring materials (dried peas) by shaking. Dotted circles inright image show flow; a little amount spilled out. See the accompanying
video or: https://youtu.be/OrjTHw0CHew

while model-based methods need to gather many samples tomake dynamic models.
We explore a temporal decomposition of dynamics tomake learning dynamics easier. For example in pouring, thewhole dynamics is complicated: the state is positions ofcontainers, types of containers, and material type, the action
is grasping parameters and flow-control parameters, and theoutput is material flow. But a small part of the dynamics,
for example the relationship between the mouth positionof the source container and the flow trajectory, is easier tomodel. We believe that humans have a model of ketchup flow
since humans can pour ketchup correctly on a target surface.
Thus we explore a method to learn temporally decomposed
dynamics and plan a behavior with dynamic programming.
Since we are working with continuous actions, we use
differential dynamic programming (DDP) which is a gradient based optimization algorithm. Since (1) learned models
typically have modeling (prediction) error, and (2) flow is
a probabilistic process, we consider probability distributionsof states and an expectation of the evaluation function (i.e.
stochastic DDP). For simplicity, we use a first-order gradientdescent algorithm. Especially when using learned dynamics,
there are many local maxima in the evaluation function, a
DDP often converges to poor local maxima. In order to avoid
them, we consider multiple criteria evaluation functions.
In addition to the stochastic evaluation function, we use
reference value functions made from reference states. The
states of the system are constrained by the dynamics, but
the reference states are optimized while ignoring those constraints. A reference value function is defined as a quadratic
centered on a reference state. DDP with reference valuefunctions performs rough optimization. A quadratic valuefunction is easier to optimize, but the solution quality is
not good due to dynamics not being taken into account toobtain the reference states. On the other hand, DDP withthe original stochastic evaluation function is useful for finetuning. Combining them gives better solutions. In addition
to this, we test multiple gradient candidates in each iterationof DDP to avoid getting stuck on saddle points.
A remarkable feature of our method is that we consider
a dynamics stage even when an action is not taken, whichallows us to decompose dynamics more flexibly. Reference
states are also considered in such stages. For example, we
consider where the flow center should be in the receiving
container, which corresponds to a reference state in our
method.
We use locally weighted regression (LWR) to learn the
dynamics (process) models [6]. In order to obtain an adequate expectation, we provide a simple extension to locallyweighted regression.
We applied our method to a pouring task on a dynamics
simulator. We used a large bouncing parameter to make theflow less predictable. Through comparing several conditions,
we found: (1) temporal decomposition of a process helps
to learn dynamics, (2) value functions with reference states
help to find a good solution especially when the processmodels are not learned, (3) our numerical method to computean expectation with LWR is useful to treat uncertainty,(4) especially when the process models are not learned,
a brute force dynamic programming using CMA-ES [7]outputs poor local maxima even with decomposed processmodels, and (5) multiple criteria optimization is helpful toavoid poor local maxima. Therefore we contribute practicalmethods for DDP with learned models, and show empirical
benefits of dynamics decomposition.
Related WorkRegarding planning behaviors under totally or partiallyunknown dynamics, there are two approaches: a model-free
approach and a model-based approach. In a model-based
approach, we learn the dynamics of the system, and then weapply dynamic programming, such as differential dynamic
programming [8]. Schaal et al. [5] used locally weighted
regression to learn a dynamics model, and applied linear
quadratic control, which was a successful approach to robot
juggling. Morimoto et al. [9] used a similar approach tolearn a dynamics model, then applied robust DDP where
the system disturbances were considered. A problem withmodel-based approaches is that we need to gather many
samples to construct a dynamics model which is especiallyproblematic in higher dimensional systems and/or systems
with complicated dynamics such as pouring. In contrast ina model-free approach, we do not learn dynamics models,
but learn policies directly. Many reinforcement learningalgorithms are proposed in this category [10], [11], [12],[3], [13]. These were successful in learning specific tasks
such as flipping a pancake, a ball-in-cup task, and a crawling
motion. However their issue is poor generalization ability
compared to that of model-based methods, although several
attempts are taken to increase generalization ability (e.g.
[2]). Another interesting approach is a hybrid of model-

Fig. 2. Temporal decomposition considered in this paper. The dotted box
denotes the whole process.

based and model-free methods. A well-known architecture is
Dyna [14]. Originally it was developed for a discrete stateaction domain, and later a new version with a linear function
approximator was developed for a continuous domain [15].
Recently Levine et al. [16] proposed a more practical approach where a trajectory optimization method for unknown
dynamics is combined with a local linear model learned fromsamples.
We introduce another insight into this model-based v.s.
model-free argument: a temporal decomposition of dynamics. Learning dynamics becomes easier when the relationship between the input and output are simple. Levine’s
approach [16] makes learning dynamics easier by extracting
a series of local dynamics (in both time and state space).
We focus on a temporal decomposition. Typically in discrete
time DDP, we divide dynamics where an action is taken.
In this case dynamics models are considered in a form of
xn+1 = Fn (xn , an ), a mapping from state xn and action
an to a next state xn+1 . We consider dynamics even whenan action is not taken: i.e. we use xn+1 = Fn (xn ) as wellas xn+1 = Fn (xn , an ) in order to make the decomposition
more flexible.
Among many DDP methods, our DDP is stochastic; thatis, we consider probability distributions of states and expectations of rewards. This design of evaluation functions issimilar to recent DDP methods (e.g. [17]). On the other hand,
we use a simple gradient descent method to optimize actions,while traditional DDP [8] and recent methods (e.g. [18], [17])use a second-order algorithm like Newton’s method. Our
DDP is inferior to second-order DDP algorithms regarding
convergence speed. The quality of the solution would be
the same as far as we use the same evaluation function.
In addition, we introduce practical methods to increase thesolution quality: multiple criteria evaluation functions withreference states and a gradient descent method with multiple
gradient candidates.
II. D IFFERENTIAL DYNAMIC P ROGRAMMING
We introduce our stochastic DDP with a first order gradient descent algorithm. We assume the dynamics are alreadydecomposed and learned. We also introduce multiple criteria
evaluation functions using reference states, and a gradientdescent method with multiple gradient candidates.
A. Problem Formulation
We consider a system with N decomposed processes
illustrated in Fig. 2. The input of n-th process is a state
xn and an action an , and its output is the next state xn+1 .
Note that some actions may be omitted. A reward is givento the outcome state. Let Fn denote the process model:
xn+1 = Fn (xn , an ), and Rn denote the reward model:
rn = Rn (xn ). Typically the most important reward is the
final reward RN (xN ) and Rn (xn ) = 0 for n < N , but weconsider a general form Rn (xn ) so that we can give rewards(or penalties) for intermediate states. We assume that every
state is observable.
We use a regression model for each Fn learned fromsamples. We assume that the regression model gives the∂Fnngradients ∂FXn = ∂F
∂xn , ∂FAn = ∂an for given xn , an .
In addition, for an input distribution xn ∼ N (x̄n , Σn ) and adeterministic an , the regression model gives the output distribution xn+1 ∼ N (x̄n+1 , Σn+1 ) as x̄n+1 = F̃n (x̄n , Σn , an )and Σn+1 = Sn (x̄n , Σn , an ) (we use F̃n (xn , an ) andSn (xn , an ) to denote them shortly), and the gradients of∂ F̃nnF̃n w.r.t. x̄n and an as ∂ F̃Xn = ∂∂F̃x̄n , ∂ F̃An = ∂an .
The purpose of dynamic programming for the n-th process
is to choose actions {an , . . . , aN −1 } so that an evaluation
function Jn (xn , {an , . . . , aN −1 }) is maximized. Jn is an
expected sum of future rewards, defined as follows:∑N
Jn (xn , {an , . . . , aN −1 }) = E[ n′ =n+1 Rn′ (xn′ )]∑N∑N= n′ =n+1 E[Rn′ (xn′ )] = n′ =n+1 R̃n′ ,(1)
where R̃n′ (xn′ ) = E[Rn′ (xn′ )]. The expectation is takenover {xn , . . . , xN } where xn′ ∼ N (x̄n′ , Σn′ ). We also
denote R̃n′ (xn′ ) = R̃n′ (x̄n′ , Σn′ ). Note that we can use
a different evaluation criteria by changing R̃n′ (xn′ ). For
example, an upper confidence bound (UCB) is given by:√
R̃n′ (xn′ ) = E[Rn′ (xn′ )] + fUCB var[Rn′ (xn′ )], (2)where fUCB is a coefficient (such as 1).
B. Simulating the Processes
At the n-th process, we estimate future states and rewards
from the state distribution xn ∼ N (x̄n , Σn ) and the actions {an , . . . , aN −1 }. The process models give xn+1 ∼N (x̄n+1 , Σn+1 ) = N (F̃n (xn , an ), Sn (xn , an )), so we obtain xn′ ∼ N (x̄n′ , Σn′ ) for n′ = n + 1, . . . , N .
In order to evaluate the expected reward, we consider
quadratic models of reward functions. Consider a Taylor
series expansion of Rn (xn ) around x̄n : Rn (x̄n + δxn ) ≈⊤δx⊤
n An δxn + bn δxn + Rn (x̄n ). Then the expectation andthe variance of Rn are given by:
E[Rn (xn )] = Rn (x̄n ) + Tr(An Σn ),(3)

var[Rn (xn )] = 2Tr(An Σn An Σn ) + b⊤n Σn bn .

(4)

The Taylor series expansion is done in each iteration of DDP.
Thus we can evaluate Jn of Eq. (1).
C. Gradient of Jn
We apply the gradient chain rule to obtain the gradient of
Jn w.r.t. {an , . . . , aN −1 }. First we consider a case of N = 2.
The evaluation function at n = 0 is given by:
J0 (x0 , {a0 , a1 }) = R̃2 (x̄2 , Σ2 ) + R̃1 (x̄1 , Σ1 )(5)

= R̃2 (F̃1 (F̃0 (x̄0 , Σ0 , a0 ), S0 (· · · ), a1 ), S1 (· · · ))+ R̃1 (F̃0 (x̄0 , Σ0 , a0 ), S0 (· · · )),(6)

where the omitted arguments of Sn (· · · ) are the same as
those of corresponding F̃n . We approximate the derivatives∂Snnby assuming ∂S∂an = ∂an−1 = · · · = 0. Then the derivativesof J0 w.r.t. a0 and a1 are obtained by the chain rule:∂JA0 = ∂ F̃A0 ∂ F̃X1 ∂ R̃X2 + ∂ F̃A0 ∂ R̃X1 ,(7)∂JA1 = ∂ F̃A1 ∂ R̃X2 ,

(8)

R̃n
n
where ∂ R̃Xn = ∂∂x, and ∂JAn = ∂J∂an . Note that since an
n
∂J0
∂Jn
′
affects only Jn′ (n ≥ n), ∂an = ∂an = ∂JAn . In a generalcase N = N , we start the calculation from an = aN −1 :∂JAN −1 = ∂ F̃AN −1 ∂ R̃XN ,∂JAN −2 = ∂ F̃AN −2 ∂ F̃XN −1 ∂ R̃XN + ∂ F̃AN −2 ∂ R̃XN −1 ,∂JAN −3 = ∂ F̃AN −3 ∂ F̃XN −2 ∂ F̃XN −1 ∂ R̃XN
+ ∂ F̃AN −3 ∂ F̃XN −2 ∂ R̃XN −1 + ∂ F̃AN −3 ∂ R̃XN −2 ,···∂JAn = ∂ F̃An Ωn ,where

Ωn = ∂ F̃Xn+1 Ωn+1 + ∂ R̃Xn+1 ,ΩN = 0.

(9)(10)

Thus we can calculate the gradient of Jn w.r.t.
{an , . . . , aN −1 } by computing backwards from n = N − 1.
The derivatives ∂ F̃Xn , ∂ F̃An are given by the processmodels. The derivative ∂ R̃Xn is obtained from a local
quadratic reward model (Taylor series expansion aroundx̄n ); ∂ R̃Xn = bn .
D. Value Functions with Reference States
We assume a reference state at state n is given as x∗n . Then
another optimization criterion is choosing an action an sothat the next state xn+1 is close to x∗n+1 . In this case, we do
not need to consider the gradient chain. For a deterministic
xn , we define a value function with a reference state as:Vn (xn ) = −(x∗n − xn )⊤ Wrs (x∗n − xn ),(11)

where Wrs is a weight matrix. Note that since x∗n is
optimized by considering Rn (x∗n ) and Vn+1 (xn+1 ), Vn (xn )doest not need to consider Rn (xn ). The expectation overa state distribution xn ∼ N (x̄n , Σn ) is E[Vn (xn )] =Vn (x̄n ) − Tr(Wrs Σn ). With the value function criteria, we
optimize an so that E[Vn+1 (xn+1 )] = E[Vn+1 (F̃n (xn , an ))]is maximized. The gradient of E[Vn+1 (xn+1 )] w.r.t. an is
obtained by:∂E[Vn+1 ]n+1 ]∗= ∂ F̃An ∂E[V∂an∂xn+1 = ∂ F̃An (2Wrs (xn+1 − x̄n+1 )),n+1where we also assumed ∂S= 0. Thus, an action is∂anupdated only considering making the next state close to thereference state; i.e. further future states are not considered.
We expect the numerical stability with such a gradient might
be more than a gradient obtained by the chain rule. For
multiple criteria optimization, we consider E[Vn+1 (xn+1 )]n+1 ]as an alternative of Jn , and ∂E[V∂xn+1 as a corresponding Ωn .

E. Gradient Descent Algorithm for Dynamic Programming
The core procedure of a gradient descent algorithm for
dynamic programming is StepGD listed in Algorithm 1.
StepGD updates current actions with multiple gradients.
Since the gradients are computed backwards from n = N −1,
StepGD calls itself recursively until n reaches N −1. Multiple
gradients computed in different ways are considered. All
gradients are tested and the one which most improves theevaluation function is actually used. We use a gradientof the evaluation function and two gradients of the valuefunctions (best two reference states are used independently).
In line A of Algorithm 1, {Ωn } ← {bn+1 , 2Wrs (x∗1
n+1 −x̄n+1 ), 2Wrs (x∗2
−
x̄
)}.
In
line
B
of
Algorithm
1,
n+1
n+1
∗2
{Ωn−1 } ← {∂ F̃Xn Ωn + bn , 2Wrs (x∗1
−
x̄
),
2W
(x
−
n
rs n
n
x̄n )}. The gradients are updated for each iteration. Note thatregardless of whether we use the evaluation function or thevalue function as criteria, we use these gradients.
PlanGD is the higher level routine listed in Algorithm 2,
which starts with an initial guess to generate a initial
sequence of actions. Then it uses StepGD repeatedly untilconvergence. In the initial guess, we search for good action
samples from a database storing every state-action-reward
from past executions. Some additional samples are generated
randomly if not many samples are found in the database. The
best action sample is chosen from these samples as an initialvalue of actions.
PlanGD has two loops. The inside loop is to run StepGD
repeatedly until it converges. The outside loop is to apply
several different criteria. In early stage of planning when theactions are far from optimal, optimizing the actions w.r.t.
value functions might be efficient because of the quadratic
objective function. On the other hand, using the evaluation
function as the optimization criterion may improve thequality of actions since this criterion is the true purpose of
dynamic programming. Thus using several criteria must be
better numerically. We use the value function criterion twice
then use the evaluation function criterion several times1 . In
PlanGD, in order to avoid a poor local maxima, Gaussian
noise is added to an action (line C of Algorithm 2).
F. Planning Reference States
We also apply a gradient descent algorithm to update
the reference states and actions. Though the states areconstrained in the actual processes, we optimize the states aswell as the actions. Here we consider deterministic variables
only. The objective criterion for optimizing a state x∗n andan action a∗n is maximizing the reward Rn (x∗n ) and making
the next state close to a reference value. This objective is
described as:
Ln (x∗n , a∗n ) = Rn (x∗n )− (x∗n+1 − Fn (x∗n , a∗n ))⊤ Wrs (x∗n+1 − Fn (x∗n , a∗n )).

(12)

1 The combination and repeat number of criteria are decided from preliminary experiments. Increasing the repeat number will increase the solution
quality but also increase the computation time. The evaluation function
criterion should be used last.

Algorithm 1: StepGD
Input:n, xn ∼ N (x̄n , Σn ), an:N −1 = {an , . . . , aN −1 }1: if n = N − 1 thenA:
Calculate {Ωn }3: else4:
Get xn+1 : x̄n+1 ← F̃n (xn , an ), Σn+1 ← Sn (xn , an )5:
a′n+1:N −1 , {Ωn } ← StepGD(n + 1, xn+1 , an+1:N −1 )6: ebest ← −∞7: for each Ωn in {Ωn } do
8:a ← an + α∂ F̃An Ωn
9:if Jn (xn , {a, a′n+1:N −1 }) > ebest then10:
a′n ← a11:
ebest ← Jn (xn , a′n:N −1 )12: if n > 0 thenB:
Calculate {Ωn−1 }14:return a′n:N −1 , {Ωn−1 }15: else
16:
return a′n:N −1 , {}Algorithm 2: PlanGD
Input:n, xn ∼ N (x̄n , Σn ), {criteria}
1: Initial guess: an:N −1
2: abest
n:N −1 ← an:N −1
3: for each criteria j in {criteria} do
4:
loop
5:
a′n:N −1 ← StepGD(n, xn , an:N −1 )6:
if Jn (xn , a′n:N −1 ) ≤ Jn (xn , an:N −1 ) then
7:
break
8:
an:N −1 ← a′n:N −1
9:
if Jn (xn , an:N −1 ) > Jn (xn , abest
n:N −1 ) then
10:
abest
n:N −1 ← an:N −1
11:
if criteria j = criteria j+1 thenC:
an:N −1 ← abest
n:N −1 + Gaussian noise
13:
else
14:
an:N −1 ← abest
n:N −1
15: return abest
n:N −1
Algorithm 3: StepGDRef∗i′
Input:n, {x∗in′ , an′ } (n = n, . . . , N , i = 1, . . . )1: if n = N then2:for each i do
∗i
3:
x∗i
n ← xn + α∂RXn
4:return5: Call StepGDRef for n + 1
6: for each i do
∗i
∗i
∗i
∗i
7:
x∗i
n ←xn +α(∂RXn +2∂FXn Wrs (xn+1 −Fn (xn , an )))
∗i
∗i
∗i
∗i
8:
a∗i
←
a
+
2α∂F
W
(x
−
F
(x
,
a
))
An
rs
n
n
n
n+1
n
n

The second term becomes zero for n = N . The partial
derivatives of Ln w.r.t. x∗n and a∗n are given by:∂Ln= ∂RXn + 2∂FXn Wrs (x∗n+1 − Fn (x∗n , a∗n )),
∂x∗
n
∂Ln
= 2∂FAn Wrs (x∗n+1 − Fn (x∗n , a∗n )).
∂a∗
n

(13)
(14)

The process models give ∂FXn , ∂FAn , and a Taylor series
expansion around the current reference state gives ∂RXn .
In order to avoid poor local maxima, we use a multi
point search. Each point is updated with a gradient descentmethod. Some search points that have smaller scores (Ln )are replaced by newly generated points, since these points
are probably getting stuck in poor local maxima. The newsearch points are generated from the database and randomly.
Algorithm 3 describes the one step update routine StepGDRef. Since the objective function includes the next state,
StepGDRef calls itself recursively until n reaches n = N ,then updates states and actions backwards.
G. Complete Action Selection
The complete action selection procedure consists of updating reference states and PlanGD. First we apply the
replacement of search points of reference states and actions.
In our experiments, we used 6 points for each state andaction, three of them are replaced by database samples,and one of them is replaced by a random value2 . Then weapply StepGDRef to update the reference states and actions.
StepGDRef is repeated several times (e.g. 3). Finally weuse PlanGD to obtain an action sequence. For exploration
purpose in an on-line learning setting, we add a zeromean Gaussian noise whose∑ variance is proportional to theuncertainty of the rewardvar[Rn (xn )].
III. L EARNING P ROCESS M ODELSIn order to model each process, we use locally weighted
regression (LWR) [6]. LWR is a memory-based method: fora query point, weights between every sample are calculated
by a kernel, then a local model is obtained by a weighted
regression. The kernel width is decided for each sample
point: assigning a small value for dense points, and a large
value for coarse points.
A. Prediction and Gradient
Let x1 , . . . , xM denote input sample vectors andy1 , . . . , yM corresponding output sample vectors. For aquery point x, an output y is obtained by:
y = β ⊤ x,

(15)

⊤

−1

β = (X WX + λI)

⊤

X WY,

⊤

(16)
⊤

where X = [x1 , . . . , xM ] , Y = [y1 , . . . , yM ] , W =diag(ϕ(x, x1 , c1 ), . . . , ϕ(x, xM , cM )), λ is a regularization
parameter, I is an identity matrix, and ϕ(x, xk , ck ) is a kernel
function with a width ck . Typically we append 1 to each input
vector of both samples and a query in order to represent
a constant term. We assume that the prediction error of
LWR is a zero-mean normal distribution N (0, Q(x)) whose
covariance matrix is given by:
Q(x) = (Xβ − Y)⊤ W(Xβ − Y)/(Tr(W)(1 − D/M )),(17)

where D is a number of dimensions of x.
In our problem setting, an input x of LWR is a concatenated vector of a state xn , an action an , and 1: x⊤ =⊤[x⊤n , an , 1], and an output y is the next state xn+1 . We use
an LWR for each process. Training an LWR is simple: adda sample of an input and an output to a memory, and update
the kernel widths.
2 These values are chosen from preliminary experiments.

We ignore the dependency of x on β to obtain the gradient;i.e. the gradient around a query point is given by β. β isa matrix consisting of gradients w.r.t. a state xn and anaction an , and a constant term: β ⊤ = [FXn , FAn , F0n ]. Thus⊤
∂FXn = F⊤
Xn and ∂FAn = FAn . For an input distribution
xn ∼ N (x̄n , Σn ) and a deterministic an , we compute theoutput distribution xn+1 ∼ N (x̄n+1 , Σn+1 ) with⊤⊤x̄n+1 = β ⊤ [x̄⊤n , an , 1] ,(18)

Σn+1 = FXn Σn F⊤
Xn + Qn (x̄n , an ),(19)

where we considered a local linear model, and we used
the MAP estimate Qn (x̄n , an ) although the covariance
Qn (xn , an ) depends on xn .
B. Improvement of Expectation
Computing the expectation using a local linear model (e.g.
Eq. (18)) or a local quadratic model (e.g. Eq. (3)) sometimes
does not work, for example when the original function is a
step function: f (x) = 1 if |x| < 0.5 otherwise 0. For x ∼N (x̄, 1), the expectation E[f (x)] computed with Eq. (18) or(3) will take 1 if |x| < 0.5, 0 if |x| > 0.5, and a strange
value if |x| ≈ 0.5. The true expectation smoothly changes
from 0 to 1 around |x| = 0.5. Applying a numerical Taylor
series expansion to f (x) with a wider window will reduce
this issue, but in order to obtain a good result, we need toobtain many samples to compute derivatives. A numericalsampling-based computation of the expectation is also time
consuming.
We solve this issue with LWR. Our simple approach isthat we use a Gaussian kernel with a diagonal covariancematrix diag(c2k , . . . , c2k ) (ck ∈ R is the width of the original
kernel), and modify it with diagonal elements of cov[x]where x is a multi dimensional query point. That is, we use adiagonal covariance matrix diag(c2k , . . . , c2k ) + diag(cov[x])to compute a weight of a query point w.r.t. each sample k.
Finally we use Eq. (18) and (19) with the new β.
This method has two advantages. One is that the computational complexity is just a single evaluation of LWR.
The other is that the gradient β obtained with the modified
covariance matrix also takes into account the uncertainty ofa query point, i.e. cov[x], which contributes a lot in our
gradient descent dynamic programming.
Fig. 3 shows an example where f (x) = 1 if |x| <0.5 otherwise 0. LWR is trained with 50 samples fromf (x) without noise. The other curves entitled with “E[f]”are expectations taken with variance 0.02 computed withthree different methods: (1) E[f (x)] obtained numerically(sampling based), (2) E[f (x)] computed similarly with (3)where the Taylor
series expansion is used numerically with√a window 0.02, and (3) E[f (x)] computed with LWR
considering the query uncertainty. (3) is closest to (1) (theyare almost overlapping), which is showing that the above
expectation improvement method works.
IV. E XPERIMENTS
We verify our method in simulated experiments. Although
our final goal is to improve actual robot performance, we use

Fig. 3.

Fig. 4.

Comparison of three calculations of E[f (x)].

Simulation environment. Right figure is after pouring material.

Fig. 5.

Fig. 6.

Typical failures of pouring.

Process model of pouring simulation experiments.

a simulator in order to compare various conditions.
We built a simulator using Open Dynamics Engine [19]where we simulate source and receiving containers, poured
material, and a robot gripper grasping the source container(Fig. 4). We simulate poured material with many (100)spheres. The gripper is modeled as fixed blocks on the source
container; we can change the grasping position, but it does
not affect the grasp quality. This gripper possibly pushes
the receiving container during pouring. In order to simulate
the complicated behavior of material such as tomato sauce
during shaking, we modified some contact model parameters;such as increasing the bouncing parameters. As a result, although the spheres are rigid objects, their trajectories during
flow are complicated and sometimes unpredictable. Typicalfailure cases are shown in Fig. 5. See also the accompanying
video or: https://youtu.be/OrjTHw0CHew
We design state machines for pouring, which is a simplified version of [4]. Those state machines have someparameters: a grasping height and pouring position. These
parameters are planned.
The process for planning is illustrated in Fig. 6. There are
five decomposed processes: x0 = (prcvx0 , prcvy0 ), a0 =(hg ), x1 = (prcvx1 , prcvy1 , hg1 ), a1 = (ppourx , ppourz ),x2 = (⟨prcvx2 − prcvx1 ⟩, ⟨prcvy2 − prcvy1 ⟩, ⟨vrcv2 ⟩, ppourx2 −prcvx2 , ppourz2 − prcvz2 ), a2 = (), x3 = (⟨pflowx3 −prcvx3 ⟩, ⟨pflowy3 − prcvy3 ⟩, wflow3 , ppourz3 − prcvz3 ), a3 =(), and x4 = (arcv4 , −aspill4 ), where (prcvx , prcvy , prcvz )is the receiving container position, hg is the graspingheight, (ppourx , ppoury , ppourz ) is the pouring position thatthe source container is rotated around during flow control,
vrcv is the speed of receiving container, (pflowx , pflowy , 0) is
the filtered flow position, wflow is the filtered flow variance,
arcv is the amount in the receiving container, and aspill isthe amount spilled out of the receiving container. The reason
we are using relative values in x2 and x3 is to reducethe modeling complexity. We apply the arctangent function
to some elements enclosed with ⟨⟩ since sometimes they
become large, which causes undesirable learning results.
The reward consists of the amount arcv4 , the spilled
penalty −aspill4 , and the penalty for the movement of the receiving container −(⟨prcvx2 −prcvx1 ⟩2 +⟨prcvy2 −prcvy1 ⟩2 +⟨vrcv2 ⟩2 ).
A. Dynamic Programming ComparisonFirst we compare the performance of DDP with parametric
optimization. Here we use process models learned beforehand; a set of process models is well learned with 90 samples(Model-1) and another set is not well learned with 30 samples (Model-2). We apply our DDP with several conditions,and also apply CMA-ES [7] directly to optimize the evaluation function w.r.t. an action sequence. The conditions are:Grad/Full: our method with five criteria at the minimum(two value function based, and three evaluation function
based), Grad/Light: our method with two criteria at theminimum (one value function based, and one evaluation
function based), Grad/NoLWREx: Grad/Full without
the expectation modification with LWR, Grad/NoDBInit:
Grad/Full without the initial guess using a database,Grad/NoRef: Grad/Full without the reference states(no value function based update, and no gradient candidates using reference states), Grad/Single: our methodwith a single criterion (evaluation function based only),
CMA/DBInit: brute force CMA-ES with the initial guessusing the database, and CMA/NoDBInit: brute force CMAES without the initial guess using the database.
Fig. 7(a) and 7(b) show the result, acquired reward andplanning time respectively. For each condition, we executed

Fig. 8. Learning curves of on-line learning. Moving average filter with 10
episode window is applied.
(a) Comparison of reward (planning quality).

Fig. 9.

Learning curves in different decompositions.

(b) Comparison of computation time.
Fig. 7.
Results of dynamic programming with well-learned models(Model-1) and not well-learned models (Model-2). Average, maximum, andminimum values are plotted.

20 trials. Briefly, using well-learned process models (Model1) is better than the other (Model-2). Using the database for
an initial guess is useful. By comparing Grad/Full withother methods, we can find that: (1) value functions withreference states helps to find a better solution especially in
the Model-2 case, (2) the expectation modification with LWR
is useful to increase the stability to obtain the solutions,(3) CMA-ES [7] sometimes outputs poor local maxima
especially in the Model-2 case.
Regarding the computation time, the gradient based methods with a small number of criteria are fast (around orless than one second). However, Grad/Full takes longer
than CMA-ES. Since by reducing the number of criteria
significantly decreases the computation time, we think we
can optimize the criteria combination to make it faster than
CMA-ES. The reason why Model-1 cases take a longer timethan Model-2 cases is that the prediction with LWR depends
on the number of samples.
B. On-line Learning
Next we explore an on-line setting. We start with nosamples, and apply the action selection and updating alternately. We compare the following conditions and methods: Sum/Grad: our method with the sum of expected

reward as the evaluation function, Sum/Grad/NoNoise:
Sum/Grad without the exploration noise in the action
selection, Sum/CMA: brute force CMA-ES with the sum of
expected reward, and UCB/Grad: our method with the UCB
as the evaluation function.
Fig. 8 shows the average learning curves of 10 trials.
There are no significant difference, but in the early stagesof learning (around 30 episodes), Sum/Grad/NoNoise isbetter than Sum/CMA. This result is the same as the previous
section: especially when the process models are not learned
well, a brute force CMA-ES goes to poor local maxima.
C. Comparison of Dynamics Decomposition
We investigate the effect of temporal decomposition of
dynamics by comparing different patterns. We conduct an
on-line setting using the same learning method, and compare
the following decomposition patterns: St5/Sum/Grad: thesame setup mentioned so far, St4/Sum/Grad: x3 is replaced by (arcv4 − aspill4 ) which is equivalent to the reward
for the final state, and x4 is removed, St3/Sum/Grad: x2
is replaced by (arcv4 −aspill4 −(⟨prcvx2 −prcvx1 ⟩2 +⟨prcvy2 −prcvy1 ⟩2 +⟨vrcv2 ⟩2 )) which is equivalent to the total rewards,
and x3 and x4 are removed.
Fig. 9 shows the learning curves averaging 10 trials. Thereis no significant difference between St4 and St5 since
the process complexity does not change by removing thestate x3 . However St3 has a big difference. The states
x2 , x3 , and x4 of the original setup are not necessary tochoose the actions a0 and a1 since only the final reward

is that we consider a dynamics decomposition even when an
action is not taken, which allows us to decompose dynamics
more flexibly. In order to learn process models, we used
locally weighted regression with an improved expectation
computation. To verify the method, we conducted simulation
experiments where we used many spheres with high bouncing parameters to simulate a complicated flow. Despite thecomplexity of the dynamics, our method worked and we got
good results.
ACKNOWLEDGMENT(a) Using the result of St5/Sum/Grad.

We thank Professor Maxim Likhachev’s Search-based
Planning Lab in Carnegie Mellon University for making their
PR2 robot available for our research.
R EFERENCES(b) Using the result of St3/Sum/Grad.
Fig. 10. Plot of estimated evaluation at x1 . Though actually the evaluationis computed for x1 = (prcvx1 , prcvy1 , hg1 ) and a1 = (ppourx , ppourz ),we changed only prcvx1 (Rcv x) and ppourx (Pour x), and used a
median value of samples for the others.

is important to optimize. These three conditions have thesame dynamics regarding x0 , a0 , x1 , a1 , and rewards, thus
if the learning methods were perfect, the results would
be the same. Fig. 10 shows an estimated evaluation at x1
plotted for varying prcvx1 (an element of the state x1 )and ppourx (an element of the action a1 ). ppourx shouldtake a proportional value to prcvx1 ; St5/Sum/Grad could
learn this, but St3/Sum/Grad could not. Thus, the bad
performance of St3 is caused by the inaccuracy of learned
process models. This result empirically supports the idea of
temporal decomposition of dynamics.
V. C ONCLUSION
We explored a temporal decomposition of dynamics inorder to enhance policy learning with unknown dynamics.
We expected decomposed processes to be easier to learn,
especially when we considered complicated dynamics, forexample pouring a range of materials like water, ketchup,and sugar. To obtain a policy, we applied differential dynamic programming using a gradient decent method withmultiple criteria. In order to take advantage of model-freemethods, we used a reference value function, which is aquadratic centered on a reference state. Reference states
are obtained through an optimization without considering
dynamics constraints. A remarkable feature of our method[1] J. Maitin-Shepard, M. Cusumano-Towner, J. Lei, and P. Abbeel,“Cloth grasp point detection based on multiple-view geometric cueswith application to robotic towel folding,” in the IEEE International
Conference on Robotics and Automation (ICRA’10), 2010, pp. 2308–
2315.
[2] J. Kober, A. Wilhelm, E. Oztop, and J. Peters, “Reinforcement
learning to adjust parametrized motor primitives to new situations,”
Autonomous Robots, vol. 33, pp. 361–379, 2012.
[3] P. Kormushev, S. Calinon, and D. G. Caldwell, “Robot motor skill
coordination with EM-based reinforcement learning,” in the IEEE/RSJ
International Conference on Intelligent Robots and Systems (IROS’10),
2010, pp. 3232–3237.
[4] A. Yamaguchi, C. G. Atkeson, and T. Ogasawara, “Pouring skills withplanning and learning modeled from human demonstrations,” International Journal of Humanoid Robotics, vol. 12, no. 3, p. 1550030, 2015.
[5] S. Schaal and C. Atkeson, “Robot juggling: implementation ofmemory-based learning,” in the IEEE International Conference on
Robotics and Automation (ICRA’94), vol. 14, no. 1, 1994, pp. 57–
71.
[6] C. G. Atkeson, A. W. Moore, and S. Schaal, “Locally weighted
learning,” Artificial Intelligence Review, vol. 11, pp. 11–73, 1997.
[7] N. Hansen, “The CMA evolution strategy: a comparing review,” in
Towards a new evolutionary computation, J. Lozano, P. Larrañaga,
I. Inza, and E. Bengoetxea, Eds. Springer, 2006, vol. 192, pp. 75–
102.
[8] D. Mayne, “A second-order gradient method for determining optimal
trajectories of non-linear discrete-time systems,” International Journalof Control, vol. 3, no. 1, pp. 85–95, 1966.
[9] J. Morimoto, G. Zeglin, and C. Atkeson, “Minimax differential dynamic programming: Application to a biped walking robot,” in theIEEE/RSJ International Conference on Intelligent Robots and Systems(IROS’03), vol. 2, 2003, pp. 1927–1932.
[10] J. Kober and J. Peters, “Policy search for motor primitives in robotics,”
Machine Learning, vol. 84, no. 1-2, pp. 171–203, 2011.
[11] E. Theodorou, J. Buchli, and S. Schaal, “Reinforcement learning of
motor skills in high dimensions: A path integral approach,” in the IEEE
International Conference on Robotics and Automation (ICRA’10), may
2010, pp. 2397–2403.
[12] D. Ernst, P. Geurts, and L. Wehenkel, “Tree-based batch mode reinforcement learning,” Journal of Machine Learning Research, vol. 6,
pp. 503–556, 2005.
[13] A. Yamaguchi, J. Takamatsu, and T. Ogasawara, “DCOB: Action spacefor reinforcement learning of high dof robots,” Autonomous Robots,vol. 34, no. 4, pp. 327–346, 2013.
[14] R. S. Sutton, “Integrated architectures for learning, planning, and reacting based on approximating dynamic programming,” in the SeventhInternational Conference on Machine Learning. Morgan Kaufmann,
1990, pp. 216–224.
[15] R. S. Sutton, C. Szepesvári, A. Geramifard, and M. Bowling, “Dynastyle planning with linear function approximation and prioritizedsweeping,” in Proceedings of the 24th Conference on Uncertainty in
Artificial Intelligence, 2008, pp. 528–536.
[16] S. Levine, N. Wagener, and P. Abbeel, “Learning contact-rich manipulation skills with guided policy search,” in the IEEE International
Conference on Robotics and Automation (ICRA’15), 2015.
[17] Y. Pan and E. Theodorou, “Probabilistic differential dynamic programming,” in Advances in Neural Information Processing Systems 27,
Z. Ghahramani, M. Welling, C. Cortes, N. Lawrence, and K. Weinberger, Eds. Curran Associates, Inc., 2014, pp. 1907–1915.
[18] S. Levine and V. Koltun, “Variational policy search via trajectory optimization,” in Advances in Neural Information Processing Systems 26,
C. Burges, L. Bottou, M. Welling, Z. Ghahramani, and K. Weinberger,
Eds. Curran Associates, Inc., 2013, pp. 207–215.
[19] R. Smith, Open dynamics engine (ODE), http://www.ode.org/.