The proposed foot contact prediction framework is verified on an electric motor actuated hopping robot DynJump, as shown in Figure 1. The hip joint is actuated directly by the motor while the knee joint is actuated through train driving. The parameters of the robot are listed in Table 1.

To focus on the study of foot contact prediction, the robot hopping motion is constrained in a vertical direction through a guide rail. The control board communicates to the motor driver through the CAN protocol at 250 Hz.

As we only focus on verifying the foot contact detection, a classical finite state machine is adopted as a high-level control scheme. The robot states are divided into the flight phase and contact phase based on the foot contact states as shown in Figure 2.

In the flight phase, a pre-calculated foot point trajectory is tracked through a position control scheme. The controller would bring the leg length to a predefined nominal leg length and leg angle for the next landing. After the leg landing, the leg motion is planned based on a virtual spring–mass model (Piovan and Byl, 2015). The detection algorithm for the state switching control will be discussed in the following section. The robot trunk vertical velocity is estimated at the landing moment to initialize the calculation of the robot trajectory after landing. As the robot is only subjected to the force of gravity in the flight phase, the initial velocity can be derived as v c = v 0 − g ⋅ ( t c − t 0 ) (1) where v _c and t _c are the velocity and time at the moment the foot contacts the ground and g is the gravity constant. For determining the exact moment when the foot leaves the ground, we actively shorten the leg at the later stage of the contact phase, and the corresponding time and robot vertical velocity are recorded as t ₀ and v ₀.

As the contact is an impact between the robot and the environment, variables associated with the interaction can be used to indicate the contact state of the foot. Though foot force can reflect the physical interaction between the foot and the ground, we can only estimate it based on the dynamics of the leg due to the absence of the foot force sensor. In robot dynamics, sensor noise and transmission clearance may introduce errors in foot force estimation. Furthermore, the kinetic parameters are always changing and unknown, so it is not very reliable if we only used the foot force to predict the contact. Although we can adopt a higher threshold to increase the reliability, a larger time delay in contact detection will degenerate the robot control performance.

The foot height and the gait cycle represent the kinematic interaction between the robot and the environment. For a running gait of the robot, like trotting gait, the robot will completely leave the ground during running, and the accurate foot height is difficult to obtain based on proprioceptive perception. Similarly, the gait cycle is highly dependent on the terrain and gait. So limited prediction performance would be obtained if we used kinematics information only.

Based on the analysis above, an indicator vector s  = [f_yh_foott_g]^T is defined to estimate the contact probability of the foot contact. Here f_y is the vertical foot force, h_foot is the foot height relative to the leaving ground level, and t_g is the gait cycle. Though Kalman filter is suitable for fusing the measuring data, a method which can learn from real data online would have advantages in environment adaption and be more robust to robot kinetic parameter uncertainty.

Though we expect that the robot can adapt to different terrains online, an off-line learning to obtain a group of initial parameters can speed up the online learning process. We used the classic two-step expectation maximization (EM) method for training. The two steps can be summarized as follows.

E-step: The responsibilities can be computed as γ i k = π k p ( x i | θ k t − 1 ) ∑ k ′ π k ′ p ( x i | θ k ′ t − 1 ) (10)

M-step: The parameters in the Gaussian component can be calculated as μ k = ∑ i r i k x i r k (11) Σ k = ∑ i r i k ( x i − μ k ) ( x i − μ k ) T r k (12) And for the next iteration, π k = 1 N ∑ i r i k = r k N (13)

In the data training procedure, a series of hopping experiments were conducted for data collection. The robot hopped on flat ground, upstairs, and downstairs. A high landing force threshold was adopted for reliable contact detection as the GMM had not yet been established. The threshold was determined based on the estimated foot force in a robot free-fall experiment, as shown in Figure 4. The corresponding fly phase time, estimated leg force, and foot height were simultaneously recorded at 200 Hz by the control board. We found that there was a non-negligible time delay in landing detection due to the high landing force threshold. We manually selected the contact data on the force cure after we collected all the experimental data. For each terrain, 30 contact data points were collected, and the learning process was done in the MATLAB 2018b environment.

As the terrain is unknown and the robot parameters may change over time, an online GMM adaption is desirable for robot control. As there are many online data stream clustering methods (Kokate et al., 2018), we would focus on dealing with the outlier point caused by the chatter effect, which balances the response time and model stability. Another problem in online clustering for leg contact detection is how to identify the sensory data corresponding to the foot–ground contact moment without the foot force sensor or external force sensors mounted on the ground. We cannot use the predicted contact results of GMM as the data are used for training the prediction model. A likelihood criterion of the foot contact event is used to find the exact contact data. This method would lead to a large time delay that cannot be used for contact prediction, but the data can then be used for online training of the GMM.

Here we use a modified version of SWEM (Dang et al., 2009). Three independent data sets are established corresponding to three GMM components. When the contact data were acquired, the probability density of the data in each of the three GMM components is calculated. The contact data would be put into the data set corresponding to the GMM component, which has the maximum probability density. When the number of newly arriving data in any of the three data sets reaches 5, the GMM will be updated as in the algorithm presented in this section.

We define that g(10) is the GMM of the data from the beginning 1 to N, and the newly arriving data are x_N+1, ⋅⋅⋅, x_N+M as in Song and Wang (2005). Different from many other applications, only one contact datum is available in one gait cycle, so we define a time window with five data points, and only one basic Gaussian model is used to model the five data points. The maximum likelihood method is used to estimate the parameters. ∂ ln L ( θ ) ∂ θ j = 0 (14) where j is the number of the parameters and L(θ) is the likelihood function, L ( θ ) = ∏ i = 1 n f ( x i , θ ) . This calculation is conducted for every five newly arriving data. When we obtain the mean and variance matrix of the Gaussian model, a T² statistic testing for the equality of the mean to the existing mean value in GMM is conducted (Song and Wang, 2005). Once the mean value is tested to be equal to any existing basic model, the Gaussian model will be merged to the corresponding basic Gaussian sub-model in GMM with the rule μ = N π j μ j + M k μ k N π j + M k (15) Σ = N π j Σ j + M k Σ k N π j + M k + N π j μ j μ j T + M k μ k μ k T N π j + M k − μ μ T (16) π = N π j + M k N + M k (17) where μ _j and Σ _j are the mean and covariance of the jth sub-component of GMM; π _j is the weight of the jth sub-component in GMM; μ , Σ , and π are the new mean, covariance, and weight, respectively; N is the total number of data that have been used for update; and M _k is the data number in the time window.

If the mean value is tested to be different from any basic Gaussian model in the GMM, the new Gaussian model will be labeled as a temporary outlier Gaussian model. For the outlier Gaussian model, it will also be merged to the nearest basic Gaussian sub-model in GMM, but a fading rule will be applied to it until it is determined as an outlier model or a shift of the existing basic Gaussian model. The merging procedure will be carried out in two steps.

Firstly, the temporary outlier Gaussian model will be merged to the existing outlier Gaussian model with the fading rule as μ = λ M 0 μ 0 + M k μ k M 0 + M k (18) Σ = λ M 0 Σ 0 + M k Σ k M 0 + M k + λ M 0 μ 0 μ 0 T + M k μ k μ k T M 0 + M k − μ μ T (19) where μ ₀ and Σ ₀ are the mean and covariance of the existing outlier Gaussian model, respectively; μ _k and Σ _k are the mean and covariance of the temporary outlier Gaussian model, respectively; μ and Σ are the new mean and covariance, respectively; M ₀ is the total number of data in the outlier set; and λ is the fading factor.

We should note that this procedure will be carried out in every model update no matter if the data are outliers or not. When the mean of the newly coming data is tested equal to the existing mean in GMM, the μ _k and Σ _k in (18) and (19) are set to zero; thus, the existing μ ₀ and Σ ₀ are faded down by factor λ.

Secondly, the resulting new outlier model will be merged with the nearest basic Gaussian component. But the parameters μ , Σ, and π of the present basic Gaussian component are reserved in the memory. When the next time window data come, if it is tested equal to any basic Gaussian component, the model will be merged as formulas (15)–(17), or if it is tested as an outlier model, the model will be moved to the outlier set and merged to the existing outlier model as formulas (18) and (19). Then the merged model will be merged with the nearest basic model.

When the outlier data situation emerges three times in succession, we assume that a concept shift happens, which means that the terrain or gait changes. Then the corresponding outlier data will be moved out of the outlier set and merged to the nearest basic component. The resulting GMM is considered as a new initial model for the newly arriving data. On the other hand, if the outlier data situation does not emerge in three consecutive time windows, the mean and the covariance will soon decrease to zeros due to the fading factor.

As we learn the contact model online without the force sensor mounted on the ground or foot, identifying which data are the contact data becomes a difficult task. Though the contact prediction module presented above will give an indication, we need to get the data in another way in order to update the prediction model and improve the prediction performance.

To collect the contact data, a “trace back” module is proposed. The algorithm includes two steps: firstly, we would detect the contact event using the GMM or leg force threshold, and then we would trace back from the contact event time to find the exact contact data.

Two situations are considered in the contact event detection: one is when the contact is triggered by the GMM, and the other is when the contact is triggered by the estimated foot force. It is obvious when the contact event happens in the first situation as the GMM will give an indication. However, when the robot encounters a new terrain or if the robot is unstable, the estimated foot force will firstly trigger the contact event. A typical estimated foot force trajectory during hopping is shown in Figure 5A. A fact we should note is that the estimated foot force changes to positive when the foot just leaves the ground, which is inconsistent with the physical truth as the foot force should either be negative or zero under the definition of foot force direction in this paper. This estimation bias occurs when the leg shortens rapidly to force the robot to enter the flight phase, which will converge to zero in about 20 sampling cycles before the next foot landing. A time criterion is added to filter this abnormal condition. A likelihood technique is used to determine if the contact event is triggered. The distribution function of the contact force is used to verify the contact probability at each time stamp. We test five consecutive foot force data to ensure that the contact event detection is reliable. The adopted distribution function is shown in Figure 5B, and the contact event criterion is defined in (20). 1 − F ( x p ) ⋅ F ( x p − 1 ) ⋅ F ( x p − 2 ) ⋅ F ( x p − 3 ) ⋅ F ( x p − 4 ) < 0.01 (20) where x_p is the coming data corresponding to the present sampling time and x_p−n is the data at n previous sampling time. The criterion means that if the misjudgment probability is below the threshold in a five-sample time window, we assume that the foot is in contact with the ground in the present sampling time. During robot running, a sliding window containing five data points is built, and the data in the sliding window will be tested if they satisfy the criterion in (20) in every control cycle. As the contact test is a conservative estimation, the resulting contact moment is several time steps after the leg contact, which is the reason why we cannot use this method to indicate the leg contact event for robot control.

As shown in Figure 6, as the lower point on the curve satisfies the contact criterion, a contact event happens. Once a contact event is confirmed, a “trace back” method is invoked to find the exact contact time moment. The contact probability is calculated in reverse order along the curve from the present data until the contact probability of the data is less than 0.5. The last data point in this procedure, with a probability of more than 0.5, will be adopted as the contact data. The “trace back” algorithm is shown in Table 2. In the algorithm, if the contact event is triggered, i.e., the trigger flag variable trigger = 1, a “trace back” procedure is employed to calculate the contact probability of F _p−k until the corresponding probability is less than 0.5.

Once five contact data are collected, a local Gaussian model is built using (14). To determine which sub-component of the GMM is to be updated, a testing procedure is adopted to determine if the mean value of the new arriving data is equal to any of the existing sub-models. Here, we use Hotelling’s T² test, which is applicable for multivariate normal data. The H0 hypothesis is the mean value of the contact data µ = µ _j, where µ _j the is one of the mean values of the sub-component of the GMM. When the sampling size is small, the T² statistics can be expressed as T 2 ∼ n ( x ¯ − μ j ) T Σ 0 − 1 ( x ¯ − μ j ) (21) where n is the number of the samples, Σ ₀ is the covariance matrix of the sub-model in GMM, and   x ¯ is the mean of the samples. If the H0 hypothesis is satisfied, we can derive that n − d d ( n − 1 ) T 2 ∼ F d , n − d (22) where F_d,n−d refers to an F distribution with d numerator degrees of freedom and n − d denominator degrees of freedom, n is number of the samples, and d is the dimension of the data vector. We can refer to the F distribution table to find out if the derived random variable in (22) obeys the F distribution, thus determining whether the H0 hypothesis is true. If the H0 hypothesis is true, the data will be used to update the GMM; otherwise, the data would be considered as outlier points if the mean is not equal to any means of the sub-components of the GMM.

For comparison, the robot hopping data of contact detection with only foot force were collected. The foot force threshold was 40 N in the hopping experiment, which means that the landing control program would be triggered when the estimated foot force exceeded 40 N. The threshold was selected based on experimental data, which could prevent an incorrect estimation of the contact state due to the evaluated error of the foot force and the noise in data measurement. A time criterion was added to the contact detection in this experiment, in which the contact event was triggered only 50 ms after the leg left the ground. Hopping data on flat ground are shown in Figure 7. Data of hopping upstairs and downstairs are shown in Figure 8; the t _delay is the time delay between the leg in contact with the ground and the leg force exceeding 40 N; t _u and t _d are the flight times of hopping upstairs and downstairs, respectively. The time delay is unavoidable because we have to ensure that the contact is not triggered by the noise or evaluated error. In hopping experiments, the time delay is about 16 ms.

To establish the initial GMM, we collected three groups of contact data, corresponding to hopping on flat ground, hopping upstairs, and hopping downstairs. In hopping-upstairs and hopping-downstairs experiments, we placed a wood plank under the foot when the robot was in the flight phase to mimic hopping upstairs and took away the plank after the robot took off to mimic hopping downstairs. The thickness of the wood plank was 1.5 cm, which was limited to the ability of the motor and transmission system. Though the height is relatively small compared with the size of the leg, it does verify the effectiveness of the prediction method as the higher the stair is, the farther different contact data in the phase space are away from each other, which will ease contact detection.

The contact data were processed in MATLAB 2018b and were selected manually in this stage. Then, the contact data were clustered using the fitgmdist function in MATLAB with a component parameter of 3. The classification result is shown in Figure 9.

After the learning process, the GMM could predict the contact time more accurately, as shown in Figure 10. In each control cycle, the maximum probability of coming data in the three components is calculated first. Then, the mean value of the coming data vector will be projected into a space where the three projected variables are independent as in (8). If all of the values are within 1.5σ, we believe that the contact happens. By using the actual contact data to renew the GMM, the model detects the contact event as soon as the leg lands.

The GMM predicted the contact in stepping upstairs and downstairs, which was also more precise than the prediction using contact force, as shown in Figure 11.

The main advantage of the learning method is that it can adapt to the change in terrains. In a terrain-changing experiment, the robot firstly hopped on flat ground, and then some wood strips were placed on the ground to mimic the changing of terrain roughness. The estimated vertical foot forces and estimated contact time are shown in Figure 12. Three gait cycles on flat ground are shown in the figure, and the GMM predicted the contact time precisely.

As the GMM model was updated only when five new gait data had arrived, the GMM prediction performance degenerated when the terrain just changed. The leg force prediction module was evoked to detect the foot landing, leading to a detection time delay of about six control cycles. In the first updating phase of the GMM model, including the last gait cycle on flat ground and the first four gait cycles on rough terrain, there was only one gait cycle where the GMM predicted the foot contact. In the second GMM updating phase, the GMM successfully predicted the foot contact in two gait cycles. After two updating phases, the GMM could always predict the foot contact with a higher time delay of about two control cycles as the variance of the GMM gets higher, as shown in Table 3.

Another effective contact detection approach is the data fusion technique, which combines different data to enhance accuracy. In the robot hopping experiments, the robot could detect the contact in one cycle delay, while the data fusion method applied to MIT cheetah 3 introduces a delay of four to five control cycles (Bledt et al., 2018a). This result shows that the proposed learning method can perform better if the robot runs on the same terrain with cyclic gait, which is true in most cases. Every animal or human would have preferred walking and running speeds, and the terrain types are limited in daily life. The gait contact parameters of hopping are distributed over a small area if the robot moves in the same terrain with the same gait.

On flat ground, the elements in the covariance matrix would be small, and the prediction would be very accurate, while in rough terrain, the absolute values of the matrix elements increase, and the prediction becomes a little rough. The most important point is that the elements in the covariance matrix will go back to being small if the terrain is a flat ground again, which means that the learning prediction method can predict the contact as accurately as possible. Owing to the online learning method, the prediction algorithm can adapt to different terrains that the robot has never encountered before.

A difficulty in learning the leg contact detection online is that the contact data for updating the GMM cannot be collected directly, as there is no force sensor to indicate the moment when the foot is just in contact with the ground. A “trace back” strategy is presented in this article. When the contact event is triggered by the GMM prediction module or leg force criterion, the update module will trace back the contact data in the memory stack until the data meet with the probability criterion.

A state machine control strategy is usually adopted to control a legged robot’s running gait. A contact event will trigger a shift in state, and the control module will change as well. So an accurate detection of the contact event is critical, which is even more important in high-speed running. In our experiment, the control cycle is 4 ms, and the learning detection module can predict the contact event four control cycles in advance. In the fast running gait of a quadruped, the total gait cycle could be less than 300 ms, and the contact period for each leg can be less than 40 ms (Hudson et al., 2012). And the contact period includes the leg compressing phase and leg extending phase, with each period being 20 ms. So the proposed prediction method would ease the control strategy design and enhance the control performance.