Kick Motions for the NAO Robot Using Dynamic Movement Primitives

Abstract

In this paper, we present the probably first application of the popular Dynamic Movement Primitives (DMP) approach to the domain of soccer-playing humanoid robots. DMPs are known for their ability to imitate previously demonstrated motions as well as to flexibly adapt to unforeseen changes to the desired trajectory with respect to speed and direction. As demonstrated in this paper, this makes them a useful approach for describing kick motions. Furthermore, we present a mathematical motor model that compensates for the NAO robot’s motor control delay as well as a novel minor extension to the DMP formulation. The motor model is used in the calculation of the Zero Moment Point (ZMP), which is needed to keep the robot in balance while kicking. All approaches have been evaluated on real NAO robots.

1 Introduction

Kick motions are an essential part of robot soccer. In recent years, the speed of the game has increased a lot with most teams now being able to stably walk at high speeds. Thus, fights for the ball are more common. A flexible kick motion that is able to adapt to different and changing ball locations as well as to different kick speeds on the fly while keeping the robot in balance during pushes from other robots is a huge advantage in such situations.

Several methods to design and execute flexible kick motions have already been developed. For instance, Müller et al. [8] model the kick foot trajectory using hand-crafted piecewise Bézier curves, which are modified on the fly to adapt to different ball positions. However, handcrafting Bézier curves is a complex and time consuming task. Wenk and Röfer [18] tackle this problem by automatically inferring trajectories based on the ball position, kick velocity, and kick direction. While this method works, it does not allow the user to influence the resulting trajectory, i. e. creating special purpose kicks like backward kicks is not possible.

In this paper, we present a middle ground between the two above-mentioned approaches: A kick motion that can be hand-crafted easily by using kinematic teach-in or be created by a multitude of optimization algorithms while retaining the ability to adapt to different ball positions and kick velocities. This is done by using a modified version of Dynamic Movement Primitives (DMPs) [4] to describe the kick trajectory. During the kick, the robot is dynamically balanced using a Linear Quadratic Regulator (LQR) with previews to keep the Zero Moment Point (ZMP) inside the support polygon. This involves a new way of estimating the ZMP based on a model of the motor behavior of the NAO.

The remainder of the paper is organized as follows: Sect. 2 introduces the motor model that is the basis of the ZMP calculation, Sect. 3 explains the ZMP estimation and introduces the balancing algorithm while Sect. 4 introduces DMPs and explains how we use them to model a kick trajectory. Sections 5 and 6 wrap up the paper with an evaluation of the kick motion and a conclusion.

2 Model-Based Motor Position Prediction

As described below, the NAO’s motor response delay is usually 30 ms. Thus, the ZMP balancer needs to take into account that the motor will not be at the currently measured position when the current command reaches the motor. Our solution to this problem is to predict the current motor position using a mathematical model and to use this prediction in our control algorithms.

2.1 Determining the NAO’s Motor Response Delay

The NAO’s motors are position-controlled using the proprietary NaoQi software. It processes the commands and relays them to an ARM-7 micro controller at 100 Hz. The controller distributes the commands over RS-485 to several dsPIC micro controllers which are responsible for controlling the actual motors. Measured motor positions travel back the same chain [3]. This chain together with the slow control rate of 100 Hz induces a delay between sending a command and being able to measure a reaction of the motor.

To determine the actual delay, a motor is moved from a resting position into a random direction and the time between sending the command and measuring a movement is recorded. A movement is registered as soon as the measured motor position deviates from the position that the motor was in when the command was issued. No threshold is used. For measuring, the internal sensor is used. This is done 100 times for each leg motor of four different NAOs, thus, we get 4400 measurements in total. As shown in Fig. 1, the vast majority of motor reactions occurs at 30 ms. The measured average distance of the reactions at 10 and 20 ms is $0.087^\circ$. This is below the maximum accuracy of the motors, which is $0.1^\circ$. Therefore, we can assume that the measurements at 10 and 20 ms are due to sensor noise. However, the measurements at 40 ms cannot be discarded as noise. Taking a closer look, it seems that some joints in some robots are more prone to responding after 40 ms than others, suggesting that hardware wear or defects might cause a delayed measurement.

Thus, for a fully repaired robot it is safe to assume that the motor response delay is 30 ms. Actually, the delay might be anywhere between 20 and 30 ms, but due to the 100 Hz duty cycle, more precise measurements are not possible.

Fig. 1.

Histograms of motor response delay. Each bar is one motor of one robot. The zoomed part shows the motors that responded after 40 ms, color-coded by robot. The experiment was conducted using the leg motors of four robots and repeating each trial 100 times. Thus, 4400 responses were measured in total.

2.2 A Model to Estimate the Motor Position

We propose to model the behavior of a motor as second order dynamical system based on a mass spring damper:

$$\begin{aligned} T^2\ddot{y}(t) + 2DT min(\dot{y}(t), V_{max}) + y(t) = u(t), \quad T, D, y, u \in \mathbb {R} \end{aligned}$$

(1)

T is the time constant, y(t) is the motor position at time t, D is a dampening constant, $V_{max}$ is the maximum motor velocity that is used to limit $\dot{y}(t)$, and u(t) is the requested motor position at time t.

The parameters $(T, D, V_{max})$ need to be set in a way that the model optimally mimics a motor. This can be achieved by minimizing the error function J:

$$\begin{aligned} J = \sum _{s \in S} \sum _{i=0}^{|s|} d(i) (m(i) - s(i))^2 \end{aligned}$$

(2)

S is a set of step responses for a given motor, |s| the number of measurements in step response s, m(i) the position of the model at the i-th step, s(i) the actual motor position at the i-th step and $d(i) = 0.85^i$ a decay function.

The decay function d emphasizes the short term model quality over the long term, i. e. we prefer parameters that provide a better short term prediction over parameters that provide an overall good prediction. This is done because in our use case the model is only used to predict a short amount of time.

Sets of step responses can be generated by applying step functions, which jump from zero to their respective values instantly, with different step heights to the motor. Figure 2(a) shows a set of step responses that has been recorded and the respective optimal model response. For fitting, the first three samples of the step response should be ignored to make up for the motor response delay.

Fig. 2.

(a) Step responses and their optimal fit. The blue lines are the step functions, red are the step responses, and green the response of the best fitting model. The model of the depicted LHipPitch joint has an error that is slightly above average (see Table 1). (b) Prediction and actual motor response of a motor while walking (Color figure online)

2.3 Model Evaluation

To show that the model can be used to predict the real world motor positions of the NAO, we compared the model response and the actual motor position while executing a five second walking motion. For the comparison, the real motor values have been shifted in time to remove the measurement delay. Figure 2(b) shows an excerpt of the experiment of the LHipRoll motor. The average absolute error over all leg joints is $0.185^{\circ }$ with a variance of 0.0573. Thus, the model seems to be able to predict the motor behavior with sufficient accuracy. Detailed results for each motor can be seen in Table 1.

Table 1. Error in leg motor predictions

Table 2. Errors when applying model parameters that have been optimized for one NAO to five different NAOs

We were also interested in the portability of the model parameters. As documented in Table 2, a repetition of the experiment on different robots (always using the same model again) was successful.

3 ZMP-Based Balancing

The robot needs to be kept in balance while kicking. A good measure for the balance of a robot is the Zero Moment Point (ZMP). A robot is said to be dynamically stable if the ZMP is inside the support polygon [17].

The center of pressure of the support polygon and the Zero Moment Point are coincident [12]. Therefore, it is possible to use the pressure sensors under the NAO’s feet to measure the ZMP, if it exists. However, the sensors are quite inaccurate and often faulty. To avoid using these sensors, we estimate the ZMP $(z_x, z_y)^T$ based on the cart-table model proposed by Kajita et al. [5]:

$$\begin{aligned} \begin{pmatrix} z_x\\ z_y\end{pmatrix} = \begin{pmatrix} c_x\\ c_y\end{pmatrix} - \frac{c_z}{g} \begin{pmatrix}\ddot{c_x}\\ \ddot{c_y}\end{pmatrix} \end{aligned}$$

(3)

$(c_x, c_y, c_z)^T$ is the center of mass (COM) and $g\approx 9.81$ the gravitational force. Due to the measurement delay and the high sensor noise, we estimate the COM using the motor model and forward kinematics $(c_x^m, c_y^m, c_z^m)^T$. The motor model is initialized using sensor readings at the beginning of the motion. While the kick motion is being executed, the model is not updated from sensor readings.

Tilting the robot over the edges of the supporting foot does not influence the estimated ZMP. To detect such situations, we calculate the scaled difference $P\varDelta \varTheta =P(\gamma - \phi )$ between the expected torso orientation $\gamma$ as provided by the motor model and the measured torso orientation $\phi$ as provided by the IMU and scale the COM accordingly [1]. The unitless constant factor P needs to be adjusted manually. We chose $P=(30, -30)^T$. Figure 3(a) shows how the ZMP behaves with and without tilt detection.

Fig. 3.

(a) The red and blue lines show the estimated ZMP position with and without tilt detection while the robot is being tilted backwards. (b) The balancing process. (Color figure online)

To be able to measure outside influences, e. g. someone pushing the robot, we replaced the COM acceleration by the acceleration of the torso $\ddot{o}$ as measured by the NAO’s IMU. This can be done because the COM is usually inside the torso and thus both accelerations are similar.

Thus, the final ZMP is calculated by:

$$\begin{aligned} \begin{pmatrix} z_x\\ z_y\end{pmatrix} = P\varDelta \varTheta \begin{pmatrix} c_x^m\\ c_y^m\end{pmatrix} - \frac{c_z^m}{g} \begin{pmatrix}\ddot{o_x}\\ \ddot{o_y}\end{pmatrix} \end{aligned}$$

(4)

To finally balance the robot, an LQR preview controller as described in [2, 16, 18] has been implemented. The inputs of the controller are the current ZMP and current COM as well as the next 50 desired ZMP positions. Figure 3(b) shows an overview of the balancing sub-system.

4 Describing Kicks Using Dynamic Movement Primitives

The kick trajectory is described by using Dynamic Movement Primitives (DMPs) [4, 13]. For the sake of simplicity, this chapter only considers one-dimensional DMPs but they can be easily scaled to n dimensions as long as the dimensions are independent of each other, i. e. for translational movements. For rotational movements, a special DMP formulation has been introduced by Ude et al. [15]. However, for kick motions it is sufficient to simply keep the foot level to the ground, thus no rotational DMP was used in this paper.

DMPs model goal-directed movements as weakly non-linear dynamical systems. They consist of the canonical system and the transformation system.

The canonical system s describes the phase of the movement:

$$\begin{aligned} \tau \dot{s}&= -\alpha _ss \end{aligned}$$

(5)

The phase s replaces the time in the transformation system. Intuitively, it drives the transformation system similar to a clock [9]. s conventionally starts at 1 and monotonically converges to zero. $\tau$ is the execution time of the movement. As long as $\tau$ remains constant, a closed solution for the canonical system exists [9]:

$$\begin{aligned} s(t) = exp(\frac{-\alpha _s}{\tau } t) \end{aligned}$$

(6)

$\alpha _s$ determines how fast s converges. The value of $\alpha _s$ has to be chosen in a way that s is sufficiently close to zero at the end of the execution. We chose $\alpha _s$ by setting $s(\tau ) = 0.01$ and solving (6) for $\alpha _s$.

The transformation system is defined by two first order differential equations:

$$\begin{aligned} \tau \dot{z}&= \alpha _z(\beta _z(g-y)-z) + sf(s)\end{aligned}$$

(7)

$$\begin{aligned} \tau \dot{y}&= z \end{aligned}$$

(8)

$\tau$ is the execution time of the movement, $\alpha _z$ and $\beta _z$ are damping constants, g is the goal position of the movement, y is the current position of the movement, s is the current phase as defined by Eq. 5 and f(s) is called the forcing term.

The dampening constants are set for critical dampening, i. e. $\beta _z = \alpha _z/4$. We chose $\beta _z=25$ and $\alpha _z=6.25$, but the exact values do not matter as long as the system is critically dampened.

The forcing term defines the movement’s shape. With $f(s) = 0$, the system is just a PD controller converging to g and reaching it at time $\tau$. One could say that f superimposes its shape onto the PD controller. To ensure that f does not keep the system from reaching the desired goal position, f is scaled by the phase, thus diminishing the influence of f towards the end of the movement. To be able to express arbitrary movements, f is typically chosen to be a radial basis function approximator:

$$\begin{aligned} f(s)&= \frac{\sum _{i=1}^N \psi _i(s)w_i}{\sum _{i=1}^N \psi _i(s)}, w_i \in \mathbb {R} \end{aligned}$$

(9)

$\psi _i$ is the i-th Gaussian radial basis function with mean $c_i$ and variance $\sigma _i^2$:

$$\begin{aligned} \psi _i(s) = exp(-\frac{1}{2\sigma _i^2}(s - c_i)^2) \end{aligned}$$

(10)

The weights $w_i$ can be chosen to create any desired function and thus define the shape of the whole movement. Different learning and optimization approaches can be used to find the weights for a certain movement, e. g. Schaal et al. [14] describe how to imitate a given trajectory.

Fig. 4.

(a) Behavior of the kick motion when the goal position is adapted. (b) Temporal scaling capability of the DMP.

Since the goal position g is part of Eq. 7, the shape also depends on g. This means that the weights can only force the system into a certain shape for one specific value of g. When g changes, the shape “bends” to reach the new goal position. In general, this is undesired behavior and is solved by scaling f if g changes. However, we found that kick motions “bend” in a natural way (Fig. 4), i. e. if the goal position is moved closer to the start position, the robot will swing further back and vice versa. This is exactly the behavior that one would expect from a kicking motion because it ensures that the distance between the inflection point and the goal position remains the same. This is important because the kicking velocity depends on this distance.

Thus, the DMP formulation allows us to define an arbitrary movement in task space and later scale its execution time as well as to move the goal position while retaining a sane shape.

A major downside of the formulation is that the final velocity is always zero making it unsuitable for kick motions because kick motions need to reach the target with a specific velocity. Solutions to this problem were proposed by Kober et al. [6] and Mülling et al. [9]. They replaced the goal g in Eq. 7 with the position, velocity and acceleration of a moving target $g_p$.

$$\begin{aligned} \tau \dot{z} = \alpha _z(\beta _z(g_p - y) + \tau \dot{g_p} - z) + \tau ^2\ddot{g_p} + sf(s) \end{aligned}$$

(11)

While Kober et al. used a target that is moving on a straight line, Mülling et al. used a fifth order polynomial.

$$\begin{aligned} g_p(t) = \sum _{i=0}^5b_it,\quad \dot{g_p}(t) = \sum _{i=i}^5ib_it^{i-1},\quad \ddot{g_p}(t) =\sum _{i=2}^5(i^2-i)b_it^{i-2} \end{aligned}$$

(12)

The coefficients $b_i$ are calculated by applying the bounding conditions:

$$\begin{aligned} g_p(t_0)&= y_0, \quad \dot{g_p}(t_0) = \dot{y_0}, \quad \ddot{g_p}(t_0) = \ddot{y_0}\end{aligned}$$

(13)

$$\begin{aligned} g_p(\tau )&= g, \quad \dot{g_p}(\tau ) = \dot{g}, \quad \ddot{g_p}(\tau ) = 0 \end{aligned}$$

(14)

Due to the time dependency, the coefficients need to be recalculated if $\tau$ changes.

Fig. 5.

Reaction of the DMP to changes in the final velocity. (a) Shows the reaction of the original DMP while (b) shows how the DMP reacts with our new scaling term A.

In this way, a new parameter $\dot{g}$ is introduced. It represents the velocity at the end of the movement. However, the weights now depend on $\dot{g}$ as well. This means that if the goal velocity is changed, the shape of the movement will change. As shown in Fig. 5(a), the trajectory reacts to changes in the goal velocity with huge changes and becomes inexecutable. We propose to fix this by scaling the forcing term with the novel factor $A = (\dot{g}_{new} - \dot{y}_0) / (\dot{g} - \dot{y}_0)$, where $\dot{y}_0$ is the starting velocity of the trajectory and $\dot{g}_{new}$ is the new goal velocity. Figure 5(b) shows that this produces much better results. If the velocity is increased, the wind up phase gets longer, if it is reduced, the wind up phase gets shorter until it completely disappears if the requested goal velocity is zero. This is exactly the behavior that one would expect from a kick motion.

Thus, the final form of the DMP used in our experiments is:

$$\begin{aligned} \tau \dot{z}&= \alpha _z(\beta _z(g_p - y) + \tau \dot{g_p} - z) + \tau ^2\ddot{g_p} + sf(s)A\end{aligned}$$

(15)

$$\begin{aligned} \tau \dot{y}&= z\end{aligned}$$

(16)

$$\begin{aligned} \tau \dot{s}&= -\alpha _ss \end{aligned}$$

(17)

This DMP responds well to changes in goal position and goal velocity. It is noteworthy that both parameters can be changed mid-execution without causing discontinuities. The implementation used in our experiments has been released as part of the B-Human Code Release 2015 [11] and is available online¹.

5 Evaluation

Several experiments have been done to evaluate the kick motion. The setup is identical for all experiments: The robot is standing at the side line of the field and kicks the ball into the field, as depicted in Fig. 7. All experiments have been done with the official RoboCup SPL ball of 2015, a Mylec street hockey ball that is 65 mm in diameter and weighs 55 g. Each experiment consists of 30 kicks.

We used the Royston H-Test [10] with a significance level of 0.05 to determine that the measured kick distances are normally distributed. Normally distributed kick results indicate that the results have only been influenced by natural noise, i. e. there is probably no systematic error in the test setup or the implementation.

To compare the performance of the kick motion to an existing one, we used imitation learning [14] to learn weights that imitate the kick motion of team B-Human of 2015 [8, 11] and executed 30 kicks with each leg and each kick motion. The results can be seen in Table 3 and Fig. 6.

Table 3. Kick distance comparision between the B-Human kick motion of 2015 and an imitation of that motion. The data for the imitated kick with the right leg contained two outliers. If those outliers are removed the result is normally distributed as well.

Fig. 6.

Kick positions of the B-Human and imitated kick motions: The red cross is the kick origin. The dots are the positions where the balls came to a halt. The blue dots originate from the B-Human kick, the red dots from the imitated kick. (a) Shows the results of kicks with the left foot while (b) shows the right foot. (Color figure online)

To test the generalization qualities of the kick motion, we conducted four experiments with different ball positions and velocities. In the first experiment, the ball is positioned 65 mm to the left. In the second experiment it is moved 80 mm forward, the third and fourth experiment reduced the kick velocity by 1/4 and 1/2 of the original kick velocity respectively. The results can be seen in Table 4. To reach the position of the left ball, the robot had to fully stretch the leg. Therefore, the knee motor could not be used to generate a forward force, thereby significantly reducing the reached kick distance. The other experiments show a reasonable scaling towards the desired kick distance. Videos showing the kick generalization can be found at https://youtu.be/g73pPCWcQvw and https://youtu.be/eANtiAiMmTg.

Fig. 7.

Ball position generalization experiments: In (a) the ball was moved 65 mm to the left, in (b) it was moved 80 mm to the front.

Table 4. Kick distance results for generalized kicks. The data of the left generalization contained one outlier. If it is removed, the result is normally distributed.

6 Conclusion

We presented a kick motion for the NAO robot that can imitate arbitrary kick trajectories and adapt them to different ball positions as well as different kick velocities. The kick motion is modeled using a slightly modified variant of DMPs. While executing the kick, the robot is kept dynamically stable using a ZMP preview controller. Additionally, we proposed a model of the NAO’s motors and used it to improve the calculation of the ZMP.

We have shown that this method of generating kick motions works but cannot kick the ball as far as a manually tuned motions. However, they are more versatile. The ball does not need to be placed perfectly to be kicked and the kick speed can be adjusted. Additionally, the underlying DMPs are easy to extend and lend themselves well to a multitude of optimization algorithms [7].

The kick motion presented in this paper has been successfully used in the corner kick challenge competition at RoboCup 2015.