Spatiotemporal synchronization of biped walking patterns with multiple external inputs by style–phase adaptation

Abstract

In this paper, we propose a framework for generating coordinated periodic movements of robotic systems with multiple external inputs. We developed an adaptive pattern generator model that is composed of a two-factor observation model with a style parameter and phase dynamics with a phase variable. The style parameter controls the spatial patterns of the generated trajectories, and the phase variable manages its temporal profiles. By exploiting the style–phase separation in the pattern generation, we can independently design adaptation schemes for the spatial and temporal profiles of the pattern generator to multiple external inputs. To validate the effectiveness of our proposed method, we applied it to a user–exoskeleton model to achieve user-adaptive walking assistance for which the exoskeleton robot’s movements need to be coordinated with the user walking patterns and environment. As a result, the exoskeleton robot successfully performed stable biped walking behaviors for walking assistance even when the style of the observed walking pattern and the period were suddenly changed.

Appendices

Appendix 1: Learning procedure of pattern generator model

Data alignment by phase information Focusing on the periodicity of the target motions, we utilize autocorrelative and cross-correlative coefficients for the alignment process. First, we maximize the autocorrelative coefficient for identifying the period T of each sequence as: $T^s \leftarrow \arg \max _j A^{s}(j)$, where $A^{s}(j) =\sum _n^N { \mathbf{y} ^s_n}^T \mathbf{y} ^s_{n+j}$ is the autocorrelative coefficient with the self-index shift j in phase with a style indexed by s. Next, we maximize the cross-correlative coefficient to find the optimal cross-index shift h in phase as: $h^s \leftarrow \arg \max _j C^{s}(j)$, where $C^{s}(j) = \sum _n^N { \mathbf{y} ^{b}_n}^T \mathbf{y} ^s_{j + \mathrm{rd}\left( \frac{n T^s}{T^b}\right) }$ is the cross-correlative coefficient and index b is the style index corresponding to the sequence that has the shortest period. $T^b$ is the period of the shortest period indexed by b. The function $\mathrm{rd(\cdot )}$ is a round-off function. The above procedures yield an aligned data matrix:

$$\begin{aligned} \mathbf{Y} ^\mathrm{all}_a = \left[ \begin{array}{ccc} \mathbf{y} ^1_{\mathrm{rd}\left( h^1+\frac{T^1}{T^b}\right) } &{} \cdots &{} \mathbf{y} ^1_{\mathrm{rd}\left( h^1+T^1\right) } \\ \vdots &{} \ddots &{} \vdots \\ \mathbf{y} ^S_{\mathrm{rd}\left( h^S+\frac{T^S}{T^b}\right) } &{} \cdots &{} \mathbf{y} ^S_{\mathrm{rd}\left( h^S+T^S\right) } \\ \end{array} \right] \end{aligned}$$

(20)

where $\mathbf{Y} ^\mathrm{all}_a = [ \mathbf{Y} _a^1 \cdots \mathbf{Y} _a^S]^T \in \mathbb {R}^{DS \times C}$, $\mathbf{Y} _a^s \in \mathbb {R}^{C \times D}$, and C is the number of aligned data samples. Note that each row of the matrix $\mathbf{Y} ^\mathrm{all}_a$ is the observations corresponding to the same value of the phase $\phi$. Each column indicates a corresponding motion sequence indexed by s.

Extraction of observation bases Since the aligned data matrix $\mathbf{Y} ^\mathrm{all}_a$ is a rectangular matrix, singular value decomposition (SVD)-based matrix factorization can be applied to extract the observation bases. By applying the factorization, we can form a style–content factorial model [referred to as the asymmetric bilinear model in Tenenbaum and Freeman (2000)].

Let ${ \mathbf{Y} ^\mathrm{all}_a}^\mathrm{VT}$ be an $S \times DC$ matrix stacked from $DS \times C$ Matrix $\mathbf{Y} ^\mathrm{all}_a$. Then, SVD for this matrix leads to the following factorial representation as

$$\begin{aligned} {{ \mathbf{Y} }^\mathrm{all}_a}^\mathrm{VT} = \mathbf{U} \mathbf{S} \mathbf{V} ^T \approx \mathbf{W} \tilde{ \mathbf{Y} }. \end{aligned}$$

(21)

We define the style parameter matrix $\mathbf{W} = [ \mathbf{w} ^1 \cdots \mathbf{w} ^S]^T \in \mathbb {R}^{S \times J}$ to be the first $J (\le S)$ rows of $\mathbf{U}$, and the content parameter matrix $\tilde{ \mathbf{Y} } = ( {[\bar{ \mathbf{Y} }^1 \cdots {\bar{ \mathbf{Y} }^J}]}^T )^\mathrm{VT} \in \mathbb {R}^{J \times DC}$ to be the first J columns of $\mathbf{S} \mathbf{V} ^T$. As a result, we can obtain an approximated form as $\mathbf{Y} ^s_a \approx \sum _{j=1}^{J} w_j^s \bar{ \mathbf{Y} }^j$. The obtained matrix $\bar{ \mathbf{Y} }^j \in \mathbb {R}^{C \times D}$ is named the jth observation basis, and the vector $\mathbf{w} ^s \in \mathbb {R}^J$ is the sth style parameter.

Learning a two-factor observation model with bases With the extracted bases $\bar{ \mathbf{Y} }^j$ for all j, we can learn the observation model. Since each point of a basis corresponds to a value of phase, we learn a mapping between $\bar{ \mathbf{y} }^j$ and $\mathbf{z}$ using all data of each basis. Here we utilize a Gaussian process regression (Rasmussen and Williams 2006) because it allows us to derive an analytically tractable predictive distribution and to learn hyper-parameters by maximization of the marginalized likelihood.

Each basis $\bar{ \mathbf{Y} }^j$ is independently modeled as a Gaussian process as $p(\bar{ \mathbf{Y} }^j| \mathbf{Z} ,\varvec{\beta }^j) \propto \exp \left( -\frac{1}{2} \mathrm{Tr} \left( ( \mathbf{K} _y^j)^{-1} \bar{ \mathbf{Y} }^j (\bar{ \mathbf{Y} }^j)^T \right) \right)$ where $\mathbf{Z}$ is the phase-aligned matrix corresponding to $\bar{ \mathbf{Y} }^j$ and $\mathbf{z} = \mathbf{h} (\phi ) = [\cos (\phi );\sin (\phi )]^T$ is a transformed variable of the phase $\phi$ to approximately measure the geodesic distance between points on $\mathbb {S}$ as the Euclidean distance in $\mathbb {R}^2$ and $\mathbf{K} ^j_y \in \mathbb {R}^{C \times C}$ is the gram matrix in which (p, q) entry is $k^j_y( \mathbf{z} _p, \mathbf{z} _q)=\beta ^j_1 \exp \left( -\frac{\beta ^j_2}{2}|| \mathbf{z} _p- \mathbf{z} _q||^2\right) + \delta _{ \mathbf{z} _p, \mathbf{z} _q}/\beta ^j_3$, and the hyper-parameter is represented as $\varvec{\beta }^j=\{\beta ^j_1, \beta ^j_2, \beta ^j_3\}$. With the above settings, the predictive distribution of the jth basis $\bar{ \mathbf{y} }^{j*}$ given a novel input $\mathbf{z} ^*$ can be easily derived as $p(\bar{ \mathbf{y} }^{j*}| \mathbf{z} ^*,\bar{ \mathbf{Y} }^j, \mathbf{Z} ) = \mathcal {N}(\bar{\mu }^j( \mathbf{z} ^*),\bar{\varSigma }^j( \mathbf{z} ^*))$ where $\bar{\mu }^j( \mathbf{z} ^*) = (\bar{ \mathbf{Y} }^j)^T ( \mathbf{K} ^j_y)^{-1} \mathbf{k} ^j_y( \mathbf{z} ^*)$, $\bar{\varSigma }^j( \mathbf{z} ^*) = \left( k^j_y( \mathbf{z} ^*, \mathbf{z} ^*) - \mathbf{k} ^j_y( \mathbf{z} ^*)^T ( \mathbf{K} ^j_y)^{-1} \mathbf{k} ^j_y( \mathbf{z} ^*) \right) \mathbf{I}$ and $\mathbf{k} ^j_y( \mathbf{z} ^*) = [k^j_y( \mathbf{z} ^*, \mathbf{z} _1) \cdots k^j_y( \mathbf{z} ^*, \mathbf{z} _C)]^T$ (Rasmussen and Williams 2006).

As the result, we extract the mean of the prediction distribution as the basis $\mathbf{f} _j(\cdot ) = \bar{\mu }^j(\cdot )$ and obtain the observation model of the adaptive pattern generator in Eq. (2) as

$$\begin{aligned} \mathbf{y} _c= & {} g(\phi _c, \mathbf{w} _c) = \sum _{j=1}^{J} w_{c,j} \mathbf{f} _j( \mathbf{h} (\phi _c)). \end{aligned}$$

(22)

Appendix 2: Two-dimensional coupled oscillator model

Here we consider the behavior of two-dimensional coupled-phase oscillator model of $\{\phi _1, \phi _2\}$ as:

$$\begin{aligned} \dot{\phi }_1= & {} \omega _1 + K_1 \sin (\phi _2 - \phi _1), \end{aligned}$$

(23)

$$\begin{aligned} \dot{\phi }_2= & {} \omega _2 + K_2 \sin (\phi _1 - \phi _2) \end{aligned}$$

(24)

where $\omega _1 > 0$ and $\omega _2 > 0$ are natural frequencies of each oscillator, and $K_1$ and $K_2$ are positive coupling gains. The dynamics of the phase difference $\varPhi = \phi _1 - \phi _2$ is given as $\dot{\varPhi } = \dot{\phi }_1 - \dot{\phi }_2 = (\omega _1 - \omega _2) - (K_1 + K_2)\sin (\varPhi )$.

Because of its simple structure, we can easily find two fixed points if $|\omega _1 - \omega _2|< K_1 + K_2$. There is no fixed point if $|\omega _1 - \omega _2|> K_1 + K_2$. A saddle-node bifurcation occurs when $|\omega _1 - \omega _2|= K_1 + K_2$. If $|\omega _1 - \omega _2|< K_1 + K_2$, the coupled oscillator runs with the phase difference $\varPhi ^* = |\phi _2 - \phi _1| = \sin ^{-1} \left( \left( \omega _2 - \omega _1 \right) /\left( K_1 + K_2\right) \right)$ and with the compromise frequency $\omega ^* = (K_2 \omega _1 + K_1 \omega _2)/(K_2 + K_1)$ (Strogtz 1994).

Appendix 3: Style–phase estimation procedures

E-step The goal of this step is to compute the estimate for the current state from the estimate for the previous time step and current observation. To this end, we first compute the predictive distribution:

$$\begin{aligned} \mathbf{x} _{t|t-1}= & {} \mathbf{A} \hat{ \mathbf{x} }_{t-1}, \end{aligned}$$

(25)

$$\begin{aligned} \varSigma _{x, t|t-1}= & {} \mathbf{A} \hat{\varSigma }_{x,t-1} \mathbf{A} ^T + \mathbf{Q} . \end{aligned}$$

(26)

Then, by using the current observation $\mathbf{y} _t$, the estimate is updated for the current time step:

$$\begin{aligned} \hat{ \mathbf{x} }_{t}= & {} \mathbf{x} _{t|t-1} + \mathbf{K} _t \left( \mathbf{y} _t - \mu ( \mathbf{x} _{t|t-1};\hat{ \mathbf{w} }_{t-1}) \right) , \end{aligned}$$

(27)

$$\begin{aligned} \hat{\varSigma }_{x,t}= & {} \left( \mathbf{I} - \mathbf{K} _t \mathbf{H} _{t|t-1} \right) \varSigma _{x,t|t-1} \end{aligned}$$

(28)

where

$$\begin{aligned} \mathbf{H} _{t|t-1}= & {} \frac{\partial \mu ( \mathbf{x} ;{\hat{ \mathbf{w} }}_{t-1}) }{\partial \mathbf{x} }\Bigg |_{ \mathbf{x} = \mathbf{x} _{t|t-1}}, \end{aligned}$$

(29)

$$\begin{aligned} \mathbf{K} _t= & {} \mathbf{A} \varSigma _{x,t|t-1} \mathbf{H} _{t|t-1}^T \left( \mathbf{H} _{t|t-1} \varSigma _{x,t|t-1} \mathbf{H} _{t|t-1}^T + \mathbf{R} \right) ^{-1}.\nonumber \\ \end{aligned}$$

(30)

M-step The goal of this step is to update the new style parameter $\hat{ \mathbf{w} }$ with the estimate of the state given in E-step and the observation at the current time step. To this end, we define the time-forgetting weighted mean as follows:

$$\begin{aligned} \ll \cdot \gg _T= & {} \eta _T \sum _{t=1}^{T} \left\{ \prod _{s=t+1}^{T} \lambda \right\} \left( \cdot _t \right) , \end{aligned}$$

(31)

$$\begin{aligned} \eta _T= & {} \left\{ \sum _{t=1}^T \left\{ \prod _{s=t+1}^T \lambda \right\} \right\} ^{-1} \end{aligned}$$

(32)

where $\lambda (0 \le \lambda \le 1)$ is a time-dependent forgetting factor which is introduced for forgetting the effect of the past observations. By using the mean, we can update the style parameter recursively:

$$\begin{aligned} \hat{ \mathbf{w} }_t = \ll {\varvec{\mu }}^T {\varvec{\mu }} \gg ^{-1}_t \ll {\varvec{\mu }}^T \mathbf{y} \gg _t \end{aligned}$$

(33)

where

$$\begin{aligned} \ll {\varvec{\mu }}^T {\varvec{\mu }} \gg _t= & {} (1-\eta _t) \ll {\varvec{\mu }}^T {\varvec{\mu }} \gg _{t-1} + \eta _t {\varvec{\mu }}(\hat{ \mathbf{x} }_t)^T {\varvec{\mu }}(\hat{ \mathbf{x} }_t),\nonumber \\ \end{aligned}$$

(34)

$$\begin{aligned} \ll {\varvec{\mu }}^T \mathbf{y} \gg _t= & {} (1-\eta _t) \ll {\varvec{\mu }}^T \mathbf{y} \gg _{t-1} + \eta _t {\varvec{\mu }}(\hat{ \mathbf{x} }_t)^T \mathbf{y} _t, \end{aligned}$$

(35)

$$\begin{aligned} \eta _t= & {} \left\{ 1 + \frac{\lambda _t}{\eta _{t-1}} \right\} ^{-1}, \end{aligned}$$

(36)

$$\begin{aligned} {\varvec{\mu }}( \mathbf{x} _{t|t-1})= & {} [\mu ^1_y( \mathbf{x} _{t|t-1}) \cdots \mu _y^J( \mathbf{x} _{t|t-1})], \end{aligned}$$

(37)

$$\begin{aligned} {\varvec{\mu }}(\hat{ \mathbf{x} }_{t})= & {} [\mu ^1_y(\hat{ \mathbf{x} }_{t}) \cdots \mu _y^J(\hat{ \mathbf{x} }_{t})]. \end{aligned}$$

(38)