Neural Networks 29–30 (2012) 8–19

Contents lists available at SciVerse ScienceDirect

Neural Networks
journal homepage: www.elsevier.com/locate/neunet

The eMOSAIC model for humanoid robot control
Norikazu Sugimoto a,b,∗ , Jun Morimoto b , Sang-Ho Hyon b,c , Mitsuo Kawato baNational Institute of Communication Telecommunication, 2-2-2 Hikaridai Seika-cho, Soraku-gun, Kyoto 619-0288, Japanb

ATR Computational Neuroscience Laboratories, 2-2-2 Hikaridai Seika-cho, Soraku-gun, Kyoto 619-0288, Japanc

Department of Robotics, Ritsumeikan University, 1-1-1 Noji Higashi, Kusatsu, Shiga 525-8577, Japan

articleinfo

Article history:
Received 16 November 2010
Received in revised form 20 September
2011
Accepted 13 January 2012
Keywords:
Modular architecture
Nonlinear and non-stationary control
problem
Humanoid robot
Computational neuroscienceabstract
In this study, we propose an extension of the MOSAIC architecture to control real humanoid robots.
MOSAIC was originally proposed by neuroscientists to understand the human ability of adaptive control.
The modular architecture of the MOSAIC model can be useful for solving nonlinear and non-stationary
control problems. Both humans and humanoid robots have nonlinear body dynamics and many degrees
of freedom. Since they can interact with environments (e.g., carrying objects), control strategies need todeal with non-stationary dynamics. Therefore, MOSAIC has strong potential as a human motor-control
model and a control framework for humanoid robots. Yet application of the MOSAIC model has been
limited to simple simulated dynamics since it is susceptive to observation noise and also cannot be applied
to partially observable systems. Our approach introduces state estimators into MOSAIC architecture tocope with real environments. By using an extended MOSAIC model, we are able to successfully generatesquatting and object-carrying behaviors on a real humanoid robot.
© 2012 Elsevier Ltd. All rights reserved.

1. Introduction
Previous studies have suggested that the human central
nervous system acquires and switches internal models of outsideenvironments to adaptively perform motor control of the bodyin various environments (Imamizu, Kuroda, Yoshioka, & Kawato,
2004; Imamizu et al., 2007). Modular selection and identification
for control (MOSAIC) architecture composed of multiple linear
state predictors and controllers was originally proposed inorder to explain the motor-control strategy of the humanbrain (Haruno, Wolpert, & Kawato, 2001; Wolpert & Kawato,
1998). Humans and humanoid robots both have nonlinear body
dynamics and many degrees of freedom (DOF). Moreover, theyboth interact with objects in real environments, which naturally
requires nonlinear and non-stationary control strategies. MOSAIC
architecture has strong potential as a human motor-control modeland a control framework for humanoid robots, and its flexible
structure enables the robot to control the nonlinear and nonstationary environment. But application of the MOSAIC model
has been limited to simple simulated dynamics, which can bepartially attributed to two problems: (1) the MOSAIC model does

∗ Corresponding author at: National Institute of Communication Telecommunication, 2-2-2 Hikaridai Seiko-cho, Soraku-gun, 619-0288 Kyoto, Japan. Tel.: +81
0774 95 1064; fax: +81 0774 95 1259.
E-mail addresses: xsugi@nict.go.jp (N. Sugimoto), xmorimo@atr.jp(J. Morimoto), gen@fc.ritsumei.ac.jp (S. Hyon), kawato@atr.jp (M. Kawato).
0893-6080/$ – see front matter © 2012 Elsevier Ltd. All rights reserved.
doi:10.1016/j.neunet.2012.01.002

not explicitly consider the existence of noise input to sensory
systems and (2) it assumes full observation and thus cannot deal
with partially observable systems. In this study, we propose an
extension of the MOSAIC architecture to cope with observationnoise and partially observable systems. Each module of MOSAIC
has a forward model, and we can adopt these to construct astate estimator. Moreover, using state estimators can provide a
reasonable model of the sensorimotor function of the central
nervous system, as previously suggested in Wolpert, Ghahramani,
and Jordan (1995). A state estimation strategy using switching
linear models is also considered to be a useful approach forestimating hidden variables of complicated nonlinear dynamics(Ghahramani & Hinton, 2000). Our new method, called extended
MOSAIC with state estimators (eMOSAIC), can deal with large
observation noise and partially observable systems. Based on theseadvantages the proposed control framework can be applied to realsystems such as humanoid robots.
A popular approach to controlling nonlinear dynamics is using
the inverse dynamics model and canceling the nonlinear terms.
Then, linear error dynamics are introduced to reduce the tracking
error. This approach is called feedback linearization (FL; Slotine &Li, 1991). FL performs well if the given inverse dynamics model isaccurate. However, it is well known that estimating the inverse
model of the real system is difficult due to the existence offriction. The controller of the proposed method can be derived
as a linear optimal controller, for which it is well known thatthe sensitivity for the modeling error is low. Thus, the proposed
method is more robust against the modeling error than FL. AnotherN. Sugimoto et al. / Neural Networks 29–30 (2012) 8–19

popular nonlinear control method is gain scheduling (GS) control(Rugh, 1991; Rugh & Shamma, 2000); The GS considers nonlinear
dynamics as a linear parameter varying (LPV) system (Shamma &Athans, 1992) defined around a pre-determined state. Although ananalysis of robustness is well studied (Chesi, 2010; Lim & How,
1998), a design of the scheduling parameter requires a knowledgeof structure of task and plant, thus it is difficult to apply the GS to
a complex environment such as a humanoid robot. On the other
hand, in our proposed method appropriate models can be selectedaccording to prediction accuracy of each linear state predictor.
Adaptive control (Astrom & Wittenmark, 1989) is a standard
method for controlling a non-stationary system. In this method,
controller parameters are adaptively modified according to thetracking error. Since this method assumes that the environmentchange is gradual, it cannot cope with a sudden change ofdynamics; e.g., carrying or grabbing extra weight. To control
such an environment, the switching control method seems to besuitable (Morse, 1996). In the proposed method, again, the module
can be switched according to prediction accuracy of each linear
state predictor.
In Section 2, we introduce eMOSAIC. In this study, we adopt an
optimal control approach in the MOSAIC model (as proposed in
Doya, Samejima, Katagiri, and Kawato (2002)).
In Sections 3 and 4, we evaluate the control performance ofthe eMOSAIC model in environments with large observation noise,
partially observable setup, and modeling errors. We investigate (1)nonlinear control performance through a squatting task (Grimes,
Chalodhorn, & Rao, 2006; Nakaoka, Nakazawa, & Ikeuchi, 2004)
using either a two-link robot model or the humanoid robot, and(2) non-stationary control performance through a object-carrying
task using either a single-pendulum robot model or the humanoid
robot with a weight. We also show that eMOSAIC significantly
outperforms the original MOSAIC in these tasks. In Section 3, weconsider the squatting task by using the two-link robot model andthe lifting task by using the single-pendulum robot model in asimulated environment. In Section 4, we apply eMOSAIC to control
our humanoid robot CB-i (Fig. 1) to demonstrate that the proposedmodel can be used in a real environment. The results show thata humanoid robot can maintain its balance while performing thesquatting and object-carrying tasks using the eMOSAIC model.
2. The eMOSAIC model
In this section, we introduce the eMOSAIC model to control ahighly nonlinear system, such as our CB-i humanoid robot (seeFig. 1), in a non-stationary environment.
2.1. Optimal control problem of nonlinear and non-stationary dynamics
We consider an optimal control problem of nonlinear and nonstationary dynamics, where the dynamics are represented as:x(t + 1) = f(x(t ), u(t ), t ) + n(t ),(1)y(t ) = h(x(t ), t ) + v(t ),(2)where x ∈ ℜ , u ∈ ℜ , and y ∈ ℜ are state, action andobservation vectors, respectively, and n(t ) ∼ N (0, Σx ) and v(t ) ∼N (0, Σy ) are system and observation noises. N (0, Σ ) denotes a
Gaussian distribution with zero mean and covariance Σ .
In an optimal control framework, the learning system tries tofind the optimal controller to minimize the objective function:
N

J =E


∞

s =0

D

L


r (x(s), u(s)) ,(3)

9Fig. 1. Humanoid robot CB-i. Height: 1.6 m. Weight: 90 kg. Robot has 51 degrees
of freedom.

where r (x, u) is the cost function. To find the optimal controller forminimizing the objective function, we estimate the value function:
V (x(t )) = E


∞



r (x(s), u(s)) .

(4)

s=t

2.2. The eMOSAIC model
The eMOSAIC model has a modular architecture. Each module
is composed of a state estimator, responsibility predictor, valuefunction estimator, and controller. We approximate nonlinear
and non-stationary dynamics by switching linear models and anonlinear cost function by switching quadratic models:x(t + 1) = Ai x(t ) + Bi u(t ) + ci + n(t ),(5)y(t ) = Hi x(t ) + v(t ),(6)ri (x(t ), u(t )) =1
2

x(t )T Qi x(t ) +1
2

u(t )T Ri u(t ),(7)where Ai ∈ ℜN ×N and Bi ∈ ℜN ×D are regression parameter of the
ith linear dynamics, ci ∈ ℜN is bias parameter, and Hi ∈ ℜL×N is
an observation matrix. Qi ∈ ℜN ×N and Ri ∈ ℜD×D are parametersof the ith quadratic cost function. Therefore, each state estimatorand controller can be represented by a linear model, and the valuefunction estimator can be represented by a quadratic model.
Fig. 2 shows a schematic diagram of the eMOSAIC model. The
responsibility predictor derives the responsibility of each module
based on the state-prediction accuracy of the state estimator.
Final output from the learning system is then derived as theweighted sum of each module’s output by the responsibility signal.
Below, we explain the details of (I) learning forward models,(II) state estimator, (III) responsibility predictor, (IV) value function
estimator, and (V) controller.

10

N. Sugimoto et al. / Neural Networks 29–30 (2012) 8–19

Fig. 2. extended MOSAIC with state estimators (eMOSAIC).

2.2.1. (I) Learning forward models
The objective of learning in each linear dynamics model isto minimize the weighted prediction error. If each model isrepresented in a linear form xi (t + 1) = Wi z(t ),1 the expectedvalues for the regression parameters Wi can be derived as⟨xz⟩i (T ).
Wi =⟨zzT ⟩i (T )(8)

The notation ⟨g ⟩i (T ) annotates a weighted mean of a function gwith respect to the responsibility signal λi (t ):T⟨g ⟩i (T ) =1
T t =1g (t )λi (t ).

(9)

The responsibility signal is a probability distribution of the module
selection (described in Section 2.2.3). By iterating the responsibility
signal calculation and parameter update, the likelihood will
increase to the suboptimal.
2.2.2. (II) State estimators
The state estimator estimates the latent states of the dynamics
from an observation. We use a linear state estimator:x̂i (t + 1|t ) = Ai x̂i (t ) + Bi u(t ) + ci ,(10)x̂i (t + 1) = x̂i (t + 1|t ) + Ki (y(t ) − Hi x̂i (t + 1|t ))(11)

where x̂i is the estimated state and Ki is the parameter of thestate estimator. We derive the parameter Ki by solving the linear
optimal estimation problem (Kalman & Bucy, 1961; Lewis, 1986)(see Appendix A).
2.2.3. (III) Responsibility predictors
The contribution of each module is represented by probability
distribution λi , which we call the ‘‘responsibility signal’’. The
responsibility signal λi is given by the following Bayes’ rule:
P (i)p(x(t ) | y(1 : t ), i)λi ( t ) = i′ ∈MP (i′ )p(x(t ) | y(1 : t ), i′ )

,

(12)

where M is the set of module indices, p(x(t ) | y(1 : t ), i) is thelikelihood of the ith module, and P (i) is the prior. By assumingthat the prediction error and estimation error are Gaussian with

1 W and z represent [A B c ] and [xT uT 1]T .
i

i

i

i

covariances Σx and Σy , the likelihood of the ith module p(x(t ) |y(1 : t ), i) is given byp(x(t ) | y(t ), i) ∝ p(y(t ) | x(t |t − 1), i)p(x(t |t − 1), i),(13)p(y(t ) | x(t |t − 1), i)= 

1(2π )L |Σy |



1
T −1
exp − ei (t ) Σy ei (t ) ,2ei (t ) = y(t ) − Hi xi (t |t − 1),

(14)
(15)

p(x(t |t − 1), i) = p(xi (t |t − 1) | xi (t − 1), u(t − 1))


1
1
T −1
= 
exp − di (t ) Σx di (t ) , (16)2(2π )N |Σx |di (t ) = xi (t ) − {Ai xi (t − 1) + Bi u(t − 1) + ci } ,(17)where x̂i (t |t − 1) is the predicted state of the state estimator of theith module and ei (t ) = y(t ) − Hi x̂i (t |t − 1) is the so-called errorof innovation.
Finally, the estimated state is derived by the weighted sum ofeach module:x̂(t ) =λi (t )x̂i (t ),

(18)

i∈M

where x̂i (t ) is the estimated state of the ith module at time t.
2.2.4. (IV) Value function estimators
We derive the controller by locally solving the linear-quadratic
optimal control problem. Since we approximate nonlinear andnon-stationary dynamics by multiple linear models in Eq. (5) andthe cost function by the quadratic functions in Eq. (7), we canlocally estimate the value function by using a quadratic function:
Vi (x(t )) =12

(x̂(t ) − xvi )T Pi (x̂(t ) − xvi ),(19)

where the matrix Pi is given by solving the Riccati equation:
1 T
0 = Pi Ai + ATi Pi − Pi Bi R−
i Bi Pi + Qi .

(20)

v

The center xi of the ith value function is given by
xvi = −(Qi + Pi Ai )−1 Pi ci .

(21)

2.2.5. (V) Controllers
By using the estimated value function in Eq. (19), we can derive
linear optimal controller for the ith module as (Doya, 2000; Lewis,

N. Sugimoto et al. / Neural Networks 29–30 (2012) 8–19

a

11

b

c

Fig. 3. (a) Two-link robot model. The mass and length of each link are 5 kg and 0.5 m, respectively, and the friction coefficient of each joint is 0.1. (b) Relationship between
the phase of the periodic pattern generator φ and the posture of the two-link robot model. For details, see Appendix D. (c) Single pendulum model with payload. The masses
of weights attaches at the corner and tip are 10 kg and 5 kg, and the lengths of long and short axes are 1.5 m and 0.5 m, respectively.

1986) as:Ri BTi Pi (x̂ tui (t ) = −v( ) − xi ).

(22)

The final output of the controller is derived by the weighted sumof each module output:u(t ) =λi (t )ui (t ).

(23)

i∈M

3. Simulations
We evaluated the proposed method by using a simple robot
model (see Fig. 3(a) and (c)). The squatting and lifting tasks
evaluate the nonlinear and non-stationary control performance,
respectively.
3.1. Squatting task

apply eMOSAIC to control the two-link robot model so that it canmaintain its balance during the squatting. Note that the two-linkrobot model cannot maintain its balance by only using the outputof the periodic pattern generator.
The angle of the lower link and the center of mass (CoM) arerepresented by θ1 and θ2 (see Fig. 3(a)). The input state vector
is x = [θ1 θ2 θ˙1 θ˙2 ]T . The output of eMOSAIC is a relative
desired trajectory of the root joint. The sum of the relative desired
trajectory and the output of the periodic pattern generator is used
to derive torque at each joint based on a PD controller. The torquecommands for the actuator are generated by the PD controller
which is given byτ = −K (q − qd ),where q = [qroot qmiddle q̇root q̇middle ]T represents the positionand velocity of the joint angles, qd is the desired trajectory, τ =T[τroot τmiddle
 ] is a torque command for the actuator of the robot,
and K =

3.1.1. Simulation setup
A basic squatting behavior is provided by a periodic patterngenerator. Fig. 3(b) shows the relationship between the phase of
the periodic pattern generator (φ ) and the posture of the twolink robot model (see Appendix D). The pattern generator outputs
a simple sinusoidal trajectory to a proportional derivative (PD)controller of the robot. The frequency of the squatting movement
is 0.2 Hz. Since a robot flexes and extends its legs periodically,a complex control law is required to prevent it from falling. We

(24)

50
0

0
50

10
0

0
10

is a gain matrix of the PD controller.

The cost is given in the quadratic form (Eq. (7)). The parametersof the quadratic cost are Qi = diag{0, 1, 0, 1} and Ri = 0.1 for all
modules. One trial lasted 10 s. The simulation and observation time
step was 0.002 s.
We evaluated the prediction error, the control cost, and thecomputation time without the observation and system noise inorder to select a number of modules. Each module segmented
the squatting movement at regular phase intervals of the patterngenerator from φ = 0 (standing position) to φ = π (crouching12

N. Sugimoto et al. / Neural Networks 29–30 (2012) 8–19

a

c

b

Fig. 4. Evaluation result of (a) average prediction error, (b) average cost, and (c) computation time for different number of modules, respectively.

a

b

c

Fig. 5. Results of squatting task. (a) Average costs with different observation noise. The means and standard deviations over 100 simulation runs are plotted. (i) The dottedline indicates the result of original MOSAIC. (ii) The dashed line indicates the results of eMOSAIC without using responsibility predictor. (iii) The solid line indicates theresults of the proposed method. (b) The trajectory of CoM (θ2 ). The dotted and solid lines indicate the original MOSAIC and eMOSAIC. (c) The responsibility signal of the
proposed method (iii) is shown. Panel A indicates the phase of the periodic pattern generator, and panels B–D indicate the responsibility signals of the three modules. The
low-passed signals are plotted with a cutoff frequency of 10 Hz.

position). For example, in the case of three modules (M = 3), eachforward model predicts the dynamics around φ = 0, 12 π , π . Then,for each module, we sampled trajectory around the corresponding
posture by using random Gaussian noise as a control input andestimated the parameters of forward model (see Eq. (8)). Fig. 4
shows the relationships between the number of modules andaverage prediction error, average control cost, and computation
time for one trial. The prediction error and the control costdecreased drastically from M = 1 to M = 3 while significant
differences between these values with M = 3 and these with M =6 were not observed. On the other hand, computation time linearly
increased when the number of modules increased. Therefore,
we decided to use three modules for this task. The acquired
parameters of the linear models are presented in Appendix B.
In this simulated environment, we focused on showing threeadvantages of the eMOSAIC model: (1) it can be applied to an
environment with large observation noise (see Section 3.1.2), (2) it
can be applied to a partially observable system (see Section 3.1.3),and (3) it can cope with modeling errors (see Section 3.1.4).
3.1.2. Robustness against large observation noises
First we evaluated the robustness of eMOSAIC against large
observation noises. For this comparison, we considered two
methods: original MOSAIC and eMOSAIC.
We tested the control performance of two methods with2
2
observation noises: Σy = σobs
I where σobs
= 0–0.12 . Fig. 5(a)shows the relationship between the size of the observation noise
and the average cost in the squatting task. The mean and standard
deviation of the cost over 100 simulation runs are plotted. Dottedand solid lines represent results acquired by the two methods,
original MOSAIC and eMOSAIC, respectively. The average cost of
the original MOSAIC rapidly increased based on according to thesize of the observation noise. In contrast, the average cost ofeMOSAIC was relatively insensitive to the size of the observationnoise. These results indicate that eMOSAIC is robust against
observation noise.
Next, we evaluated the effect of the prior information for themodule selection. If the linear modules are sparsely allocated, theresponsibility signal tends to chatter around a border of the local
state spaces specified by the linear models, due to the softmax
calculation of Eq. (12). Chattering of the responsibility signals
causes damage to the hardware in the real environment. In order
to avoid the hardware damage, the temporal continuity, whichassumes that the modules switch smoothly, as prior information ofmodule selection can be introduces. We describe the details of thetemporal continuity in Appendix C. The natural frequency of thependulum represented by CoM at the standing position (φ = 0) is
0.705 Hz. We empirically found that the time constant parameterof the temporal continuity ρ around 10–20 times of the natural
frequency leads good control performance. In this simulation, we
selected the time constant parameter of the temporal continuity
as a 15 times of the natural frequency (ρ = 15×10.705 ). DashedN. Sugimoto et al. / Neural Networks 29–30 (2012) 8–19

13

a

b

c

Fig. 6. Results of squatting task in full observation and partial observation environments. (a) Trajectory of CoM (θ2 ). The solid and dashed lines represent full and partial
observation cases, respectively. (b), (c) Squatting movements of full and partial observation cases (see Movie 1).

line represents the performance of eMOSAIC with the temporalcontinuity as the prior information. The performance was slightly
better than the case of without the temporal continuity. This result
indicates that introduction of the temporal continuity dose notsignificantly deteriorate control performance. Note that we use
P (i) = 1/M where M = 3 as prior information for the case without
using the temporal continuity.
Fig. 5(b) plots the trajectory of θ2 (CoM angle) with the2observation noise with variance σobs
= 0.022 . The dotted and solidlines represent control performance of the original MOSAIC and theeMOSAIC, respectively. eMOSAIC showed better performance inmaintaining balance than the original MOSAIC. Fig. 5(c)-A show thephase of the periodic pattern generator (φ ), and Fig. 5(c)-B–D show
the responsibility signal of each module (i = 1, 2, 3) estimated byeMOSAIC. eMOSAIC was able to switch each module periodicallyfor the periodic squat movement.
Through this simulation we demonstrated that the robot
model could maintain balance and switch the modules in the
noisy environment using eMOSAIC. These results indicate thatthe adaptation capability of the MOSAIC architecture can be
implemented in a real environment.
3.1.3. Application to partially observable environments
In the real environment, there is no guarantee that all states canbe measured. Here, we apply eMOSAIC to both a fully observableand partially observable
 system (only joint angles can be observed,
i.e., Hi =

1
0

0
1

0
0

0
0

in Eq. (6)).

Fig. 6(a) plots the trajectories of CoM (θ2 ) in the fully observableand partially observable system. In the case of full observation, theangle of CoM could be maintained around zero. Table 1 shows themean squared errors, and Fig. 6(b) and (c) show the sequence of thesquatting movements for full observation and partial observation.
In the both cases, eMOSAIC was able to maintain the balance of the
two-link robot (see Movie 1). In the partially observable system,
however, the CoM angle deviated from zero, and the robot could
avoid failure even if the angular velocities of each joint could not
be observed.
Note that the two-link robot model cannot maintain its balanceby only using the output of the periodic pattern generator.

Table 1
Mean squared errors of squatting task in fully observable andpartially observable system.

MSE

Full observation

Partial observation

3.66 × 10

4.18 × 10−4

−6

3.1.4. Robustness against modeling errors
For the real system it is difficult to precisely estimate
physical parameters. Here, we evaluate the proposed method’s
robustness against modeling errors and compare the robustnessof the controller with that of a standard nonlinear tracking
controller using feedback linearization (FL) (Slotine & Li, 1991) (seeAppendix E).
FL cannot balance the CoM of the robot with the desired
trajectory used by eMOSAIC since that was a simple sinusoidal
trajectory (see Fig. 3(b)). Therefore we modified the desired
trajectory of FL (qdes ) to maintain the angle of CoM (θ2 ) zeroprecisely. (For details of the modified desired trajectory, see
Appendix D.) We optimized the parameter κ that represents theconvergence property of the tracking error in order to minimize
the tracking error: κ = 5.88 (see also Appendix D).
To simulate the modeling errors, we multiplied the parameter
α to actual mass and length of the robot links. Thus α = 1means a zero modeling error. We evaluated the robustness of bothmethods against the modeling errors (α = 0.6–2). The variance ofthe system and the observation noises ware fixed to small values,22
σsys= 0.0012 and σobs
= 0.012 . The results are plotted in Fig. 7,
in which the horizontal and vertical axes are the parameter α
and the mean of the CoM angle (θ2 ). The eMOSAIC shows better
tracking performance than FL except for the case of α ̸= 1 in which
there is no modeling error. In the implementation of eMOSAIC wederived the controller of each module as an optimal linear feedback
controller. It is well known that the optimal linear feedback control
law reduces the sensibility for the modeling error (Cruz & Perkins,
1964; Fujii & Narazaki, 1984; Perkins & Cruz, 1971). The tracking
error of FL was significantly large when the parameter α was lessthan one. This is because FL output the smaller than required
torque, and then failed to maintain the balance.

14

N. Sugimoto et al. / Neural Networks 29–30 (2012) 8–19

We again evaluated the effects of the temporal continuity.
Dashed line in Fig. 8(a) represents the performance of eMOSAIC
with the temporal continuity. The natural frequency of singlependulum without payload is 0.407 Hz, we also selected the time
constant parameter ρ as 15 times of the natural frequency (ρ =1). The performance was better than the case of without
15×0.407
the temporal continuity. This result indicates that our proposed
method can switch module fast enough even with the temporalcontinuity with the selected time constant parameter.
Fig. 8(b) and (c) represent the responsibility signals of original
MOSAIC and eMOSAIC. Original MOSAIC was not able to switch
the module according to the payload placement due to the large
observation noise. These results indicate that the state estimator
for the responsibility predictor is essential for coping with noisy
non-stationary environments.
Fig. 7. Results of robustness test of eMOSAIC and feedback linearization. The
horizontal axis is the parameter α defining a modeling error. The vertical axis is
the angle of CoM (θ2 ). The mean and standard deviation of 50 trials are plotted.

3.2. Lifting task

4. Real robot experiments
Both for the squatting task and lifting task, we implemented
eMOSAIC for the CB-i real robot (see Fig. 1).
4.1. Squatting task

3.2.1. Simulation setup
Here, we consider the lifting task as an example of a nonstationary control problem. Fig. 3(c) shows single-pendulum
models with and without payload. For the pendulum model withthe payload, the weights attached at the corner and tip represented
the mass of the robot’s body and the payload. The task was tokeep the balance in response to the existence and absence of the
payload.
One trial was 6 s and the payload was attached during the last
3 s. The state vector is x = [θ θ̇ ]T , and the output is the desired
trajectory of the joint. The torque command for the actuator isgenerated by the PD controller given byτ = −K (x − xdes ),(25)

where xdes is the desired trajectory, τ is a torque command forthe actuator, and K = [50 10] is the gain parameters of the PD
controller.
The cost is given by a quadratic formation:r (t ) =1
2

x(t )T Q x(t ) +1
2

u(t )T Ru(t ),(26)

where Q = diag{1, 1} and R = 0.1. Thus the minimum cost wasgiven at upright position θ = 0 rad.
Since there are two conditions: with and without payload,
we use two modules (i = 1, 2). The first and second modules
predict a single-pendulum model without and with the payload,
respectively. We sampled data around a standing position byusing Gaussian random input absence and presence of the payload
for the first and second modules respectively and estimatedthe parameters of forward model (see also Eq. (8)). Obtained
parameters are showed in Appendix B. We again test theperformance of the two control methods: original MOSAIC andeMOSAIC. We evaluate the robustness of two methods with2
2
observation noises (Σy = σobs
I with σobs
= 0–0.012 ).
3.2.2. Results
Fig. 8(a) shows the relationship between the observation noise
and the average cost in the lifting task. Dotted and solid lines
represent the performance of original MOSAIC and eMOSAIC. The
proposed method shows the higher performance than the original
MOSAIC. The performances of the original MOSAIC decreased
drastically with large observation noise while eMOSAIC did not.

4.1.1. Experimental setup
We sought to maintain the balance of the CB-i in the squatting
task. Fig. 9(a) shows the pitch-joint coordination of the CB-i. A
periodic pattern generator output the desired trajectories to thehip, knee, and ankle joints (see Appendix F).
We approximated the dynamics of CB-i using a two-link robot
with a foot model (Atkeson & Stephens, 2007; Hyon, Osu, &Otaka, 2009; Stephens, 2007) (Fig. 9(b)). The input state is x =Tθtorso θankle θgyro θ̇torso θ̇ankle θ̇gyro . The outputs of eMOSAIC are add add T
the relative desired torso and hip joint angles: u = θtorso
θankle .



These angles are added to the desired joint angles generated by theperiodic pattern generator. To follow the desired joint angles, thetorque output at each joint is derived by a PD controller:τ = −K (q − qd ),(27)
Twhere q = [qtorso qhip qknee qankle q̇torso q̇hip q̇knee q̇ankle ] represents the joint angles’ position and velocity, qd is the desired trajectory, τ = [τtorso τhip τknee τankle ] is torque commands for the
robot’s actuators, and K is a gain matrix of the PD controller.
Here, the purpose of each linear optimal controller is to keep
the pitch angle of the gyro sensor close to zero θgyro = 0, where thesensor is attached to the hip of CB-i (see Fig. 9(b)). The cost function
is defined by quadratic form (see Eq. (7)). The parameter for state
deviation is given as Qi = diag{0, 0, 1, 0, 0, 1} to minimize theangle of the gyro (θgyro ) and the angular velocity of the gyro(θ̇gyro ), and the parameter for output deviation is given as Ri =diag{0.01, 0.01} (i = 1, 2, 3) to minimize the control cost.
We prepared three modules: i = 1, 2, 3, same as a simulationof squatting task. The forward models of the first, second, andthird modules estimated states around three different postures.
We sampled data around the posture corresponding to the phaseof φ = 0, π2 , π by using smoothed Gaussian random input andestimated the parameters of forward models (see also Eq. (8)).
Since frequent module switching cause damage to the hardware, we introduced temporal continuity of the responsibility
signal (see also Appendix C). The height of CoM of CB-i at a standing position is about 0.98 m, the natural frequency of the pendulum model represented by CoM and the center of pressure (CoP) is0.504 Hz. By following the simulation results, we selected the time
constant parameter of the temporal continuity as a 15 times of thenatural frequency (ρ = 15×10.504 ).

N. Sugimoto et al. / Neural Networks 29–30 (2012) 8–19

a

b

15

c

Fig. 8. Results of lifting task. (a) Average costs with different observation noise are displayed in the same format as in Fig. 5. Figs. (b) and (c) show the results obtained byoriginal MOSAIC and eMOSAIC. Each figure shows the responsibility signals of the module. The low-passed signals are plotted with a cutoff frequency of 10 Hz.

a

b

c

Fig. 9. (a) Joint coordinates of humanoid robot CB-i. (b) Two-link robot with foot model of the CB-i. The two joints of the two-link model represent the torso and ankle jointof the CB-i. The gyro sensor is attached to its hip. (c) Lifting task. Payload mass is 2.9 kg. Left: without payload. Right: with payload.

4.1.2. Results
Fig. 10 shows the results of the 0.5 Hz squatting task. Fig. 10(a)A shows the phase of the periodic pattern generator (φ ). Fig. 10(a)B shows the trajectory of the pitch angle of the gyro sensor (θgyro ).
The solid line represents the eMOSAIC result (see also Movie 2). The
pitch angle of the gyro sensor (θgyro ) is maintained approximately
 2
2
π , 180
π , within which the robot does notin the range of − 180
fall over. The dashed line represents the result obtained by using
only the first module (see also Movie 3). At approximately t = 3 s,
the robot in which only the single module was used fell backward.

By using the multiple-module, we successfully balanced the robot
because the second and third modules were subsequently selected.
Fig. 10(a)-C–E show the responsibility signals of three modules.
Fig. 10(b) shows the squatting movements by using multiple
modules. The upper body of the CB-i was stable even if the lower
body was flexed and extended. In contrast, in the case of a single
module (see Fig. 10(c)), the upper body moved backward andforward. The robot could not stand without support. The robot wasable to maintain the balance with squatting frequency up to 1.5 Hz
by using the eMOSAIC model (see also Movie 4).

16

N. Sugimoto et al. / Neural Networks 29–30 (2012) 8–19

a

b

c

Fig. 10. Results of 0.5 Hz squatting task in real-robot experiments. (a) Panel A shows the phase of the periodic pattern generator (φ ). Panel B shows the pitch angle of
the gyro sensor. The solid line represents the results of the proposed method (Movie 2), and the dashed line represents the results obtained by using only the first module(Movie 3). The other panels C, E show the responsibility signals of the proposed method. (b), (c) Squatting movement of 0.5 Hz by using multiple modules and single module,
respectively. The robot was able to maintain balance during the squatting task with frequency up to 1.5 Hz (see also Movie 4).

4.2. Lifting task
4.2.1. Experimental setup
We tried to maintain the balance of the CB-i in the lifting task
(Fig. 9(c)). Again, we approximated the dynamics of our humanoid
robot by using the two-link robot with a foot model (Fig. 9(b)). We
used two modules (i = 1, 2) in eMOSAIC. The first and the second
models corresponded to the state of the CB-i absence and presenceof the payload, respectively. We again sampled data by usingsmoothed Gaussian random input and estimated the parametersof forward model (see also Eq. (8)). For the module selection, we
assumed the temporal continuity of the responsibility signal. We
selected the same time constant parameter to the squatting task
(ρ = 15×10.504 ).
In this task, each linear optimal controller tried to balance
the pitch angle of the gyro. Therefore, we used the quadratic

cost function the same as for the squatting task (Qi = diag{0,0, 1, 0, 0, 1} and Ri = diag{0.01, 0.01}, i = 1, 2).
4.2.2. Results
Fig. 11(a) shows the results of the lifting task. Fig. 11(a)-A shows
the pitch angle of the gyro sensor. We applied a payload to thetray at a time of t = 0 s. The solid and dashed lines represent the
proposed method (see also Movie 5) and the single-module method(see also Movie 6), respectively. The robot with a single module
fell forward at approximately t = 3 s. However, by the proposed
method, the pitch angle of the gyro sensor θgyro was successfully2
maintained at less than 180
π rad. Fig. 11(a)-B and C show the
responsibility signal of the modules for no-payload and payload,
respectively. First, the proposed method selected the first module,
which controls the environment without a payload. After applying
the payload, the second module, which controls the environment

N. Sugimoto et al. / Neural Networks 29–30 (2012) 8–19

17

a

b

c

Fig. 11. Results of lifting task in real-robot experiments. (a)-A The pitch angle of the gyro sensor in the multiple-module and single-module cases (see also Movie 5 and Movie
6). B and C show the responsibility signals of the proposed method. (b), (c) The movement of the lifting payload by using multiple modules and single module, respectively.
The robot was able to maintain the balance with module switching even if the payload was applied, but it could not maintain balance without module switching.

when a payload is used, was automatically selected. Fig. 11(b) and(c) show the sequence of the movements during the lifting task.
The robot was able to maintain balance even when the payload wasapplied, but could not maintain balance without module switching.
5. Conclusion
We have extended MOSAIC architecture by using state estimators. We compared the proposed method with the previous
methods with large observation noise, partially observable environment, and modeling errors. The simulation results indicate that
our extended MOSAIC architecture improves robustness against
large observation noise and partial observation. We farther investigated robustness of our proposed method against modeling errors by comparing with the feedback linearization method. We also
showed that we can roughly estimate proper number of modules
by checking the average prediction errors. The parameters of each
module can be estimated from generated data by using Gaussian
random input. We then implemented our proposed method on the
CB-i humanoid robot. We tried to stabilize the CB-i in a squatting
task and lifting task. In these tasks, the CB-i was not able to maintain balance by only using the single module. Our proposed method
however successfully stabilized the robot.
The squatting task is suitable to evaluate the control performance of nonlinear control methods. Grimes et al. (2006) proposed a imitation leaning algorithm and evaluated the algorithm

on the squatting task. They transferred the captured human motion to the humanoid robot via a latent state space composedof low-dimensional posture information, gyro sensor information,and foot pressure information. However, this imitation learning
method can only be used to generate feed-forward squatting trajectories. Therefore, the control performance can be susceptible toexternal disturbances. On the other hand, our proposed method
consists of multiple linear feedback controllers. Thus, our proposed
controller can be robust and suitable to generate fast movements.
Nakaoka et al. (2004) proposed a motion planning algorithm andthe evaluation includes squatting task. In this motion planning
method, a whole body motion sequence can be composed by combining motion primitives. Since the motion primitives need to be
designed for a particular motion sequence, generalization performance of this planning method can be limited. On the other hand,
our proposed method uses linear dynamical systems as primitive
representations. Therefor, the primitives do not directly depend on
a particular movement and the proposed method may have better
generalization performance.
Using modular architecture is a practical approach toward
dealing with a large state space. Many robotics and machinelearning studies proposed modular control methods (Doya et al.,
2002; Haruno et al., 2001; Morimoto & Doya, 2001; Samejima,Doya, & Kawato, 2003; Wiering & Schmidhuber, 1997). ‘‘Macroactions’’ (Hauskrecht, Meuleau, Kaelbling, Dean, & Boutilier, 1998;

18

N. Sugimoto et al. / Neural Networks 29–30 (2012) 8–19

Sutton, Precup, & Singh, 1999) or ‘‘sub-goals’’ (Dietterich, 2000;
Morimoto & Doya, 2001) have been used in hierarchical learning
systems. However, in the most of the previous works, each module
is represented by complex nonlinear functions. In this study, each
module is represented by a linear model. Therefore, parameterfor the controller can be solved. Thus, proposed method can beconsidered as more practical approach than the previous methods.
The cost function cannot be always modeled by a quadratic
form for all tasks. In such cases it is necessary to approximate
the complex cost function. As a future study, we will consider
a hierarchical learning algorithm. In the hierarchical system, thehigher layer can be used to properly allocate local quadratic costfunctions to approximate the complex cost function.
Acknowledgments
This research was supported by the Strategic Research Programfor Brain Sciences (SRPBS). This work was partially supported by
MEXT KAKENHI 23120004. The authors gratefully acknowledge N.
Nakano for assistance with the experimental setup.

The linear optimal state estimator (i.e., the Kalman filter)
assumes that dynamics and observation are linear (see Eqs. (5) and(6)), and system noise and observation noise (n(t ) and v(t )) are
zero-mean Gaussian. The optimal state estimator estimates latentstate by alternating two phases: prediction and update (Thrun,
Burgard, & Fox, 2005). The prediction phase is as follows:(A.1)

Σ (t + 1|t ) = AΣ (t |t )A + Σx ,(A.2)
where Σx is the covariance of system noise, x(t + 1|t ) and Σ (t +1|t ) are the mean and covariance of the predicted latent state attime t + 1 based on the information of previous step t. A, B and cT

are the parameter of linear dynamics in the same format as Eq. (5).
The update phase is as follows:
K = Σ (t + 1|t )H (H Σ (t + 1|t )H + Σy )The periodic pattern generator for the two-link robot model(see Fig. 3) outputs a simple sinusoidal trajectory at each joint:
qd1 = D1 − cos φ

qd2 = −2D

+ qrest1 ,2
1 − cos φ2+ qrest2 ,qd1 = − arctan

M2 Lg2 sin q2
M1 Lg1 + M2 L1 + M2 Lg2 cos q2(D.3)

Feedback linearization (FL) uses an inverse dynamics model tocancel nonlinear terms of original dynamics. Linear error dynamics
are then introduced. FL considers following nonlinear dynamics:
H (q)q̈ + C (q, q̇)q̇ + g (q) = τ,

(E.1)

where q = [q1 q2 · · ·]T is the position in a general coordinate, andτ = [τ1 τ2 · · ·]T is the control output. The control law to track thedesired trajectory qd is given byτ = H (q)v + C (q, q̇)q̇ + g (q),where v is the equivalent input:x(t + 1) = x(t + 1|t ) + K (y − Hx(t + 1|t )),(A.4)

Σ (t + 1) = (I − KH )Σ (t + 1|t ),where Σy is the covariance of observation noise.

(A.5)

v = q̈d − 2κ q̃˙ + κ 2 q̃,−1.

Appendix E. Feedback linearization

(A.3)

T

(D.2)

where qd1 and qd2 are the desired trajectory of the first and second
joints, qrestand qrest
are rest positions, and D is amplitude of a12periodic movement.
The modified desired trajectory of the first joint which is usedby the feedback linearization is given by

,

T

(D.1)

The two-link robot can maintain a balance of CoM (θ2 ) by tracking
this desired trajectory.

Appendix A. Linear optimal state estimatorx(t + 1|t ) = Ax(t |t ) + Bu(t ) + c

Appendix D. Periodic pattern generator for two-link robot
model

(E.2)

(E.3)

q̃ = q − q is the position tracking error and κ is positive constant(Slotine & Li, 1991).
d

Appendix B. Parameters of linear models
Table 2 shows parameters of linear models (Ai , Bi , and ci ), which
we obtained for squatting task in the real robot environment.
Table 3 also shows obtained parameters of linear models for liftingtask in the real robot environment.

Appendix F. Squatting controller
The squatting controller outputs desired trajectories at the hip,dddknee and ankle joints. Each desired trajectory (θhip, θknee, and θankle)
is given by following the central pattern generator (CPG),

Appendix C. Prior distribution of responsibility predictor
We introduce the prior probability P (i) based on ‘‘temporalcontinuity’’, for the module selection. The prior can be given by
1
P (i) ∝ √exp − 2 Ei (t − 1) .
2σ
2π σ 21(t − s)1t
Ei ( t ) =exp −∥ei (s)∥2 1t ,τ
s=0

1 − cos φ(t )
2

d
θknee
= 2D

(C.1)

Ei represents a smoothed square error (weighted sum from t = 0
to current time) of innovation ei (t ) at time t:
t


d
θhip
=D

(C.2)

d
θankle=Drest+ θhip,1 − cos φ(t )2

1 − cos φ(t )
2

(F.1)

rest
+ θknee
,

(F.2)

rest
+ θankle
,

(F.3)

rest
rest
where D is the amplitude of a squatting movement. θhip, θkneeandrestθankle define the rest posture. φ is the phase of the periodic patterngenerator.

where 0 < ρ is a parameter that controls the strength of thetemporal continuity and 1t is a time step of the observation.
Eq. (C.2) can be expanded to the following recursive form:

Appendix G. Supplementary data1t
Ei (t ) = ∥ei (t )∥ 1t + exp −Ei (t − 1).
ρSupplementary material related to this article can be foundonline at doi:10.1016/j.neunet.2012.01.002.

2



(C.3)

N. Sugimoto et al. / Neural Networks 29–30 (2012) 8–19

19

Table 2
Linear model parameters (Ai , Bi , and ci ) for squatting task in the real robot environment.
i

1

Ai

3

0.00000.00000.0000
0.0023
0.00020.0001ci
0.00000.0000
0.0000
0.0003
0.0059
−0.0001

1.0000 0.0000 0.0000−0.0801 0.0045−0.01030.00001.00000.0000−0.0043−0.2328
0.06680.00000.00001.0000−0.0054−0.0511−0.00930.00200.00000.00000.9947−0.0011−0.00030.00000.00200.00000.00120.9814−0.00180.00000.0000
0.0020
0.0001
−0.0124
0.9996

1.0000 0.0000 0.0000−0.0979
 0.0030−0.0140
1.0000 0.0000 0.0000−0.0986
 0.0003−0.00140.00001.00000.0000−0.0009−0.34220.06950.00000.00001.00000.0003−0.0731−0.0089

0.00200.00000.00000.9950−0.0027−0.0005

0.00000.00200.00000.00110.9770−0.00220.00000.0000
0.0020−0.0006
−0.0146
0.99590.00001.00000.0000−0.0044−0.2937
0.09840.00000.00001.00000.0010−0.0477−0.0057

0.00200.00000.00000.9952−0.00210.00030.00000.00200.00000.00220.98220.0004
0.00000.00000.0000 0.0000
0.0000 0.0000
0.0028 0.0002
0.0004 0.0093
0.00020.0001
0.00000.0000 0.0000 0.0000
 0.0000 0.0000
 0.0028 0.0004
 0.0006 0.0066
−0.0003 0.0005



2

Bi







0.00000.0000
0.0020
0.0001
−0.0069
0.9945
0.0000 0.0000
 0.0000−0.0175
−0.0202
0.0046
0.0000 0.0000
 0.0000−0.0216
−0.0902
0.0165
0.0000 0.0000
 0.0000−0.0237
−0.1308
0.0449

Table 3
Linear model parameters (Ai , Bi , and ci ) for lifting task in the real robot environment.
i

1Ai
1.0000 0.0000 0.0000−0.0886 0.0239−0.00380.00001.00000.0000−0.0005−0.39990.04910.00000.00001.00000.0004−0.0039−0.0191

0.00200.00000.00001.0001−0.00500.00021.0000 0.0000 0.0000−0.0867−0.0264−0.01700.00001.00000.0000−0.0029−0.3879
0.04210.00000.00001.0000−0.0001
0.0131−0.0075

0.00200.00000.00000.99940.00020.0037


2Bi

0.00000.00200.00000.00140.9873−0.01110.00000.00200.00000.00030.9872−0.0094

References
Astrom, K. J., & Wittenmark, B. (1989). Adaptive control. Addison Wesley.
Atkeson, C., & Stephens, B. (2007). Multiple balance strategies from one
optimization criterion. In The IEEE-RAS 2007 international conference onhumanoid robots (pp. 57–64).
Chesi, G. (2010). LMI techniques for optimization over polynomials in control: asurvey. IEEE Transactions on Automatic Control, 55, 2500–2510.
Cruz, J. B., & Perkins, W. R. (1964). A new approach to the sensitivity problem in
multivariable feedback system design. IEEE Transactions on Automatic Control,
9, 216–223.
Dietterich, T. G. (2000). Hierarchical reinforcement learning with the MAXQ valuefunction decomposition. Journal of Artificial Intelligence Research, 13, 227–303.
Doya, K. (2000). Reinforcement learning in continuous time and space. NeuralComputation, 12, 219–245.
Doya, K., Samejima, K., Katagiri, K., & Kawato, M. (2002). Multiple model-based
reinforcement learning. Neural Computation, 14, 1347–1369.
Fujii, T., & Narazaki, M. (1984). A complete optimality condition in the inverse
problem of optimal control. SIAM Journal on Control and Optimization, 22(2),327–341.
Ghahramani, Z., & Hinton, G. E. (2000). Variational learning for switching statespace models. Neural Computation, 12(4), 831–864.
Grimes, D. B., Chalodhorn, R., & Rao, R. P. N. (2006). Dynamic imitation in a humanoid
robot through nonparametric probabilistic inference. In Proceedings of robotics:
science and systems, RSS’06. MIT Press.
Haruno, M., Wolpert, D. M., & Kawato, M. (2001). Mosaic model for sensorimotor
learning and control. Neural Computation, 13.
Hauskrecht, M., Meuleau, N., Kaelbling, L. P., Dean, T., & Boutilier, C. (1998).
Hierarchical solution of Markov decision processes using macro-actions. In
Uncertainty in artificial intelligence (pp. 220–229).
Hyon, S., Osu, R., & Otaka, Y. (2009). Integration of multi-level postural balancing
on humanoid robots. In IEEE international conference on robotics and automation(pp. 1549–1556).
Imamizu, H., Kuroda, T., Yoshioka, T., & Kawato, M. (2004). Functional magnetic
resonance imaging examination of two modular architectures for switching
multiple internal models. Journal of Neuroscience, 24(5), 1173–1181.
Imamizu, H., Sugimoto, N., Osu, R., Tsutsui, K., Sugiyama, K., Wada, Y., et al. (2007).
Explicit contextual information selectively contributes to predictive switchingof internal models. Experimental Brain Research, 181(3), 395–408.

0.00000.00000.0000
0.0003
0.0001
0.00010.00000.0000
0.0020
0.0004
0.0147
0.9960

0.00000.0000
0.0020−0.0010
0.0098
0.9979

0.0000 0.0000 0.0000 0.0007−0.0004
0.0001ci
0.00000.0000
0.0000−0.0003
0.00690.0007
0.0000 0.0000
 0.0000−0.0382
 0.1012
−0.0125
0.0000 0.0000
 0.0000−0.0365
 0.0765−0.0169

0.00000.0000
0.0000−0.0002
0.0065
0.0008Kalman, R. E., & Bucy, R. S. (1961). New results in linear filtering and prediction
theory. Transactions of the ASME Series D-Journal of Basic Engineering, 83(1),95–108.
Lewis, F. L. (1986). Optimal estimation. John Wiley & Sons, Inc.
Lim, S., & How, J. P. (1998). Control of LPV systems using a quasi-piecewise affine
parameter-dependent Lyapunov function. In American control conference.
Morimoto, J., & Doya, K. (2001). Acquisition of stand-up behavior by a real robot
using hierarchical reinforcement learning. Robotics and Autonomous Systems,
36(1), 37–51.
Morse, A. S. (1996). Lecture notes in control and information sciences, Control using
logic-based switching. Springer-Verlag Telos.
Nakaoka, S., Nakazawa, A., & Ikeuchi, K. (2004). An efficient method for composing
while body motions of a humanoid robot. In Proceedings of international
conference on virtual systems and multimedia. VSMM2004.
Perkins, W. R., & Cruz, J. B. (1971). Feedback properties of linear regulators. IEEE
Transactions on Automatic Control, 16(6), 659–664.
Rugh, W. J. (1991). Analytical framework for gain scheduling. IEEE Control Systems
Magazine, 11, 79–84.
Rugh, W. J., & Shamma, J. S. (2000). Research on gain scheduling. Automatica, 36,
1401–1425.
Samejima, K., Doya, K., & Kawato, M. (2003). Inter-module credit assignment inmodular reinforcement learning. Neural Networks, 16, 985–994.
Shamma, J. S., & Athans, M. (1992). Gain scheduling: potential hazards and possible
remedies. IEEE Control Systems Magazine, 6, 101–107.
Slotine, J.-J. E., & Li, W. (1991). Applied nonlinear control. Prentice Hall.
Stephens, B. (2007). Humanoid push recovery. In The IEEE-RAS 2007 international
conference on humanoid robots (pp. 589–595).
Sutton, R. S., Precup, D., & Singh, S. P. (1999). Between MDPs and semi-MDPs:
a framework for temporal abstraction in reinforcement learning. Artificial
Intelligence, 112(1–2), 181–211.
Thrun, S., Burgard, W., & Fox, D. (2005). Probabilistic robotics. The MIT Press.
Wiering, M., & Schmidhuber, J. (1997). HQ-learning. Adaptive Behavior, 6(2),219–246.
Wolpert, D. M., Ghahramani, Z., & Jordan, M. I. (1995). An internal model for
sensorimotor integration. Science, 269(5232), 1880–1882.
Wolpert, D. M., & Kawato, M. (1998). Multiple paired forward and inverse models
for motor control. Neural Networks, 11, 1317–1329.