Optimized Kalman Filter based State Estimation and Height Control in Hopping Robots

Abstract

Rotor-based hopping locomotion significantly improves efficiency and operation time as compared to purely flying systems; where most hopping robots use the liftoff states and an assumed ballistic trajectory to determine the hopping height. However, significant aerial phase force (e.g., thrust and drag) can invalidate this assumption and lead to poor estimation performance. To combat this issue, a group has implemented multiple sensors (active and passive optical, inertial, and contact) and significant computational power to achieve full state estimation. This, however, poses a significant challenge to the development of light-weight, high-performance, low observable, jamming and electronic interference resistant hopping systems; especially in perceptually degraded environments (e.g., dust, smoke). Here we show a training procedure for a coupled hopping phase and Kalman filter-based vertical state estimator, requiring only inertial measurements, which is able to learn the characteristics of the target system, sensors, locomotion behaviors, environment, and acceleration measurement aliasing conditions. The resulting estimator, given hop heights up to 4 m and velocities up to $\pm 7$ m/s, achieves a mean absolute percent error in the hop apex height of 12.5% with an aerial trajectory average normalized mean absolute error in position and velocity of 19% and 16.5%, respectively; while operating at 840 Hz, on a dual-core 240 MHz processor, with a total robot mass of 672 g. Due to the low mass and computational power, the presented estimator could also be used as a degraded operational mode in cases of sensor damage, malfunction, or occlusion in more complex robots.

I Introduction

Hopping robots combine the agility of aerial systems with the efficiency of terrestrial systems [1, 2, 3]; where the focus has been placed on predicting the vertical states for control, from the liftoff characteristics [4, 5]. However, this method is not applicable to the recent generation of rotor-based hopping robots that may not follow strictly ballistic trajectories in the aerial phase [1, 6, 3, 7, 2]. Therefore, a new method for continuous state estimation is necessary. This is particularly challenging in light-weight hopping robots, without the payload for significant increases in sensors, processing, and power; furthermore high-performance hopping robots, where aliasing of sensors measurements can be significant, especially during the impacts associated with hopping. The challenge is further compounded in perceptually degraded environments (e.g., dust, smoke, darkness), in contested environments in which active jamming may be present, and where sensing modalities that radiate energy into the environment may lead to detection and destruction. Currently, robot rotational state estimation can be achieved through inertial measurement units (IMUs), and translational state estimation can be achieved through GPS and combined optical flow and range finder senors; however, range finder sensors radiate energy into the environment for sensing, and can have difficultly operating in perceptually degraded environments. As IMUs are ubiquitous in robots and small in size, this work seeks to develop an IMU-based vertical state estimation methodology for hopping height control that could be extended to horizontal states in the future; fully replacing the optical flow sensor for operation in highly perceptually degraded environments.

To date numerous non-continuous jumping and hopping robots have been developed including miniature jumping robots [8, 9, 10, 11, 12, 13, 14], high specific energy robots [15, 16], soft jumping robots [17, 18], multimodal robots [19], with wings for gliding [20, 21, 22, 23, 24, 25, 26, 27, 28], and inertial tails or reaction wheels for orientation control [29, 30, 31, 32, 33]. These robots have also used a wide variety of hopping leg designs including: linkage [34, 35, 36, 37], and linear [38, 39, 40, 41] legs. These robots represent a highly diverse range of characteristics and constraints. However, given the minimal impact of an IMU integration, the addition of the proposed method is highly feasible for any system.

Currently, five untethered continuous hopping robots exist: MultiMo-MHR (our system) [1], PogoDrone [6], Hopcopter [2], Salto/Salto-1P [42, 43, 36, 44, 30, 31, 32, 33], and PogoX [3, 7] and one tethered continuous insect-scale hopping robot [45]. While PogoDrone and the insect-scale hopping robot shows no form of state estimation, both the Hopcopter and Salto robots use the estimated or measured liftoff states and an assumed ballistic aerial trajectory to estimate the hopping height. However, the ballistic trajectory assumption invalids this technique for robots with applied aerial phase forces. PogoX (robot mass = 3.65 kg) displays the most thorough state estimation to date, using an extended Kalman filter, moving horizon estimation, and visual inertial odometry, requiring multiple sensors, advanced computational hardware, large weight budget, and direct measurement of both vertical position and acceleration to estimate vertical velocity; thus posing a challenge to broad applicability. This work seeks to estimate both position and velocity using only direct acceleration measurements allowing for a broader integration across platforms while also being used in a degraded operational mode when primary sensors are damaged or malfunction in more complex systems; such as PogoX.

This paper is organized as follows. Section 1 introduces the topic, and Section 2 briefly covers the robot design with focus on the changes from previous work [1]. Section 3 discusses the hopping phase estimator, and Section 4 presents the hopping state estimator with associated training procedure and results. Section 5 explores the effects of the sensing frequency on measurement aliasing and ultimately state estimation. Section 6 presents the experiential procedure and results of the deployed state estimation technique on the MultiMo-MHR, and Section 7 summarizes the work.

Figure 1: Hopping robot. a) Operational phases (drop, stance down, stance up, rebound) and transitions (touchdown (TD), max squat (MS), liftoff (LO), hop apex (HA)) a.1) Photo of the MultiMo-MHR platform with components labeled. a.2) CAD model with overall dimensions labeled. a.3) CAD model section view with spring extension (stance phases) showing the combined body and leg structures. b) Dynamic simulation showing the locally, highly non-linear, behavior of the body at TD and LO for a 3 m hop.

II Robot Design

The MultiMo-MHR, developed in previous work [1], is a monopedal hopping robot that utilizes rotors for orientation control and energy input (Fig. 1a, Appendix Robot Parameters). The robot consists of two bodies where a single unactuated prismatic joint, with the motion-axis aligned to the leg-axis, connects the body to the leg. At touchdown the body’s momentum causes it to move down the leg, extending the main power springs. The springs consist of multiple rubber bands that are able to store a maximum usable locomotion energy of approximately 34 J resulting in a maximum hopping height of over 5 meters; where additional bands can increase the height.

The MultiMo-MHR uses an on-board microcontroller (ESP32-WROOM-32, two Tensilica Xtensa LX6 cores at 240MHz) for real-time control (Fig. 1a). Given the significant difference in the accelerations experienced during ground contact and the aerial phases, the platform has two inertial measurement units (IMUs), including a low-g (MPU-6050: $\pm 16$g, $\pm$2000 deg/s, 1000 Hz) and high-g (H3LIS331DL: $\pm 100$g, 1000Hz), ensuring accurate measurements tailored to the varying magnitudes of forces encountered; where the static measurement noises are $\pm 0.013g$ (MPU-6050) and $\pm 0.32g$ (H3LIS331DL). A data logger (Sparkfun OpenLog) is used to store both the sensor data and calculated values during operation. The robot is powered by a 14.7v battery pack (two 7.4v, 2s, 2200 mAh, 35C batteries in series) to provide sufficient power to the four electronic speed controllers (ESC, 20 Amp) and motors (SpeedyBee 1404 4500 kv). Finally, a radio frequency transceiver (nRF24L01) receives data from an external PC which is used to simulate additional sensors, via a motion capture system (Vicon vantage V5), while focusing on hopping height control (e.g. PMW3901 Optical Flow Sensor [horizontal positions] and WT901 [orientation]).

The previously developed controller, composed of the horizontal [X,Y], orientation [roll,pitch,yaw], and locomotion energy (LE) input [Z or height] control, is augmented with a vertical state estimator and a new LE input controller for hopping height control (Fig. 2). Due to the highly damped nature of the vertical states, a proportional controller is used for the LE input control. The P-controller output is mapped to a thrust-to-weight ratio (TWR) of between 0% and 83.7% (i.e., 100% equates to the same thrust required for hover); where the minimum motor command is set to 5% duty-cycle to increase the motor’s responsiveness while producing almost no thrust. The controller operates as a sliding mode control at ranges greater than 0.33 m from the desired height.

Figure 2: Hopping robot control block diagram.

Figure 3: State machine diagram of the Hopping Phase Estimator (HPE), Hopping Vertical State Estimator (HVSE), and controllers showing both the prediction updates and Inferred Measurement Updates (IMUPTs) of the HVSE. Shown are phases (green), phase transition criteria (red), and transitions (blue)

III Hopping Phase Estimator (HPE)

A hopping phase estimator (HPE) is developed (Fig. 3) where the operational phases (green), phase transition criteria (red), and transitions (blue) are shown. The HPE determines the drop, stance down, stance up, and rebound phases through identification of the four transitions; touchdown (TD), maximum squat (MS), liftoff (LO), and hop apex (HA). Previous work in hopping robots has indicated the TD and LO transitions can be determined from the states [4], with two specific methods including leg motion through the use of a leg encoder [42, 30, 31, 32, 33] and ground contact through a force sensor on the foot [3, 7]. However, both methods require an additional sensor with complex mounting requirements beyond the typical 9-axis inertial measurement unit (IMU) required for aerial phase control. Instead, inspired by work in inertial pedestrian navigation or pedestrian dead-reckoning [46, 47, 48, 49], multiple acceleration measurements (IMU) are used to calculate the slope (jerk) and a set positive threshold determines TD. This same measurement is then used to determine MS as the point at which the jerk becomes negative. LO is determined by a single filtered acceleration measurement that is less than zero; indicating gravity is the dominant acceleration. Finally, previous works have predicted the HA from measured liftoff states and an assumed ballistic trajectory. However, as hopping height increases so to will drag, resulting in increased deviation from gravity-based ballistic motion. Furthermore, rotor-based hopping robots can significantly alter their aerial trajectory from that of ballistic motion; due to thrust produced in the aerial phase. While range finding sensors can directly measure the height and thus the HA, they significantly increase the observability of the system, provide and opening for external interference, and increased height can degrade the measurement. As proprioceptive sensors are unable to directly or indirectly detect it, the transition criteria for the HA requires an accurate estimation of the vertical velocity.

TABLE I: Kalman Filters

Filter	Filter Equations	Filter Parameters	IMUPTs	Errors
Type			(Phase/Transition)	(998 hops)
	KF - Prediction	$z_{k}=H\hat{x}_{k}+n_{a},\,u_{k}=[0,0,1]\,a_{k}^{w}$
KF1	$x_{k}^{-}=Fx_{k-1}+Gu_{k}$ $P_{k}^{-}=FP_{k-1}F^{T}+Q$	$F=\begin{bmatrix}1&\Delta t\\ 0&1\end{bmatrix},\,G=\begin{bmatrix}0.5\Delta t^{2}\\ \Delta t\end{bmatrix},\,x_{k}=\begin{bmatrix}z\\ \dot{z}\end{bmatrix}$	Position (TD, LO) Velocity (MS, LO)	$\gamma_{1}=0.108$ $\gamma_{2}=0.329$ $\gamma_{3}=0.594$
	KF - Meas. Update	$z_{k}=H\hat{x}_{k}+n_{a},\,u_{k}=[0,0,1]\,a_{k}^{w}$
KF2	$K=P_{k}^{-}H^{T}(HP_{k}^{-}H^{T}+R)^{-1}$ $x_{k}=x_{k}^{-}+K(z_{k}-Hx_{k}^{-})$ $P_{k}=(I-KH)P_{k}^{-}$	$F=\begin{bmatrix}1&\Delta t&-0.5\Delta t^{2}\\ 0&1&-\Delta\\ 0&0&1\end{bmatrix},\,G=\begin{bmatrix}0.5\Delta t^{2}\\ \Delta t\\ 0\end{bmatrix}$	Position (TD, LO) Velocity (MS, LO) Accel. Bias (Aerial Phases)	$\gamma_{1}=0.111$ $\gamma_{2}=0.497$ $\gamma_{3}=0.903$
		$x_{k}=[z\,\dot{z}\,\ddot{z}_{bias}]^{T}$
	ESKF - Nominal State Prediction	$z_{k}=H\hat{x}_{k}+n_{a}-Hx_{k}^{-},\,u_{k}=[0,0,1]\,a_{k}^{w}$
ESKF1	$x_{k}^{-}=f(x_{k-1})+g(x_{k-1})u_{k}$ ESKF - Prediction $P_{k}^{-}=FP_{k-1}F^{T}+Q$	$F=\begin{bmatrix}1&\Delta t&-0.5\Delta t^{2}\\ 0&1&-\Delta t\\ 0&0&1\end{bmatrix},\,G=\begin{bmatrix}0.5\Delta t^{2}\\ \Delta t\\ 0\end{bmatrix}$	Position (TD, LO) Velocity (MS)	$\gamma_{1}=0.127$ $\gamma_{2}=0.417$ $\gamma_{3}=0.724$
	ESKF - Meas. Update
ESKF2	$K=P_{k}^{-}H^{T}(HP_{k}^{-}H^{T}+R)^{-1}$ $\delta x_{k}=\delta x_{k}^{-}+K(z_{k}-H\delta x_{k}^{-})$ $P_{k}=(I-KH)P_{k}^{-}$ $x_{k}=x_{k}^{-}+\delta x_{k}$	$x_{k}=[z\,\dot{z}\,\ddot{z}_{bias}]^{T}$ $\delta x_{k}=[\delta z\,\delta\dot{z}\,\delta\ddot{z}_{bias}]^{T}$	Position (TD, LO) Velocity (MS) Accel. Bias (Aerial Phases)	$\gamma_{1}=0.134$ $\gamma_{2}=0.423$ $\gamma_{3}=0.725$

• Parameters: state transition matrix, $F$, control matrix, $G$, estimated covariance, $P$, process noise covariance, $Q$, observation matrix, $H$, measurement noise covariance, $R$, Kalman gain, $K$, input, $u_{k}$, and measurements, $z_{k}$. Abbreviations and descriptions: Measurement (Meas.), Acceleration (Accel.), Kalman Filter (KF), Error State Kalman Filter (ESKF), true state ($\hat{x}_{k}$).

IV State Estimation Methodology

Two common approaches involve a Kalman filter (KF) with states ($x_{k}=[p_{k},v_{k},\theta_{k},b_{ak},b_{gk}]$) or an error state Kalman filter (ESKF) with error states ($\delta x_{k}=[\delta p_{k},\delta v_{k},\delta\theta_{k},\delta b_{ak},\delta b_{gk}]$) and nominal states ($x_{k}$); where, the subscript $k$ is the time step, $p_{k}$ is the position, $v_{k}$ is the velocity, $\theta_{k}$ is attitude, $b_{ak}$ is the accelerometer bias, $b_{gk}$ is the gyroscope bias, and $\delta$ indicates the error between the true and estimated states. The inertial measurement unit (IMU) and accelerometer provide both the linear acceleration and rotational velocity as,

$$a^{IMU} =\hat{a}+b_{a}+n_{a}+\alpha_{a}$$ (1)

$$\omega^{IMU} =\hat{\omega}+b_{g}+n_{g}$$ (2)

where, $\hat{a}$ is the true acceleration, $n_{a}$ is the accelerometer noise, $\alpha_{a}$ is the accelerometer aliasing, $\hat{\omega}$ is the true rotational rate, and $n_{g}$ is the gyroscope noise. As with [50], neglecting the Coriolis forces and curvature of the earth, the state transition equations are as follows,

$$p_{k+1}^{w}=p_{k}^{w}+v_{k}^{w}\Delta t+\frac{1}{2}((R_{k}^{bw}(a_{k}^{IMU}-b_{ak})+g^{w})\Delta t^{2}$$ (3)

$$v_{k+1}=v_{k}+(R_{k}^{bw}(a_{k}^{IMU}-b_{ak})+g^{w})\Delta t$$ (4)

$$R_{k+1}^{bw}=R_{k}^{bw}(I+[\omega^{IMU}_{k}-b_{gk}]_{\times}\Delta t)$$ (5)

$$b_{ak+1}=b_{ak}+n_{ba}$$ (6)

$$b_{gk+1}=b_{gk}+n_{bg}$$ (7)

where, the superscript $w$ indicates the world frame, $R^{bw}$ is the rotation matrix from the body to the world frame, $I+[\omega^{IMU}_{k}-b_{gk}]_{\times}\Delta t$ is the first-order Taylor expansion of Rodrigues’ rotation formula, $a^{IMU}$ is the measured acceleration in the IMU frame, $\omega^{IMU}$ is the measured rotation rate in the IMU frame, the subscript $\times$ indicates the cross-product operator matrix for the bias compensated angular rate measurement, and $g^{w}$ is the gravity vector.

Developing a KF from the state transition equations 3-7, shows that measurement updates impact only the measured state and its time integrals and derivatives; thereby allowing for extraction and independent development of the KF for the vertical states. Extending this to an ESKF, it does show a small amount of coupling between non-integral and non-derivative states, due to the attitude error impacting the velocity error [47]. However, assuming small attitude errors and therefore limited coupling, this isolation and development can be applied to both KFs and ESKFs. This allows for a focused study on managing significant and variable aliasing in acceleration measurements (MultiMo-MHR: Nyquist frequency $\sim$420 Hz). The challenge is further compounded by the fact that the most significant aliasing occurs at the regions of maximum energy transfer during TD and LO; as seen in the previously developed model (Fig. 1b) and experimental testing [1].

IV-A Hopping Vertical State Estimator (HVSE)

Differing from [6, 2], with thrust-to-weight ratios (TWRs) $\geq$ 1, the MultiMo-MHR operates with a maximum TWR limited to 83.7%. Operating at TWR $<$ 1 poses an additional challenge as the true drop phase transition is determined by the dynamic behavior and not the controller; where an early estimated drop phase transition reduces the maximum performance while a late drop phase transition reduces hopping efficiency. Vertical state estimation at a TWR of $\sim$90% has been shown, however, the work required both position and acceleration to estimate vertical velocity [3, 7].

Previous studies in pedestrian tracking employing Kalman Filters (KF) to estimate both position and velocity, have shown that assuming zero velocity of the foot during ground contact can allow for inferred velocity measurements, termed zero-velocity updates (ZUPTs), which can significantly reduce drift in the state estimation [49, 46]. Further studies have added zero-altitude updates [47] to reduce the drift in the vertical axis as well. Expanding on the concept for vertical state estimation in high-performance hopping robots, three of the four phase transitions including touchdown (TD), max squat (MS), and liftoff (LO) as well as the motor drive commands allow for inferred-measurements updates (IMUPTs) including:

IMUPT-position (TD,LO)

$=[0,0,1]\,R^{bw}[0,0,-L_{f}]^{T}$ at TD and LO where the vertical position is inferred from the known distance between the robot’s center-of-mass and the foot ($L_{f}=(m_{B}L_{3}+m_{L}L_{1})/(m_{B}+m_{L})-(L_{1}+L_{2})$) without spring extension.

IMUPT-velocity (MS)

$=0$ at MS where the vertical velocity is inferred from the body’s change in direction.

IMUPT-velocity (LO)

$=v_{LO}\delta_{vLO}$ at LO where the vertical velocity is inferred from an understanding of the system losses, momentum transfer, and inherent acceleration aliasing at liftoff. The scaling factor $\delta_{vLO}=(c_{vel2}v_{LO}^{2}+c_{vel1}v_{LO}+c_{vel0})(c_{ch1}h_{ch}+c_{ch0})$; where, liftoff velocity estimate, $v_{LO}$, and preset commanded height for the subsequent hop, $h_{ch}$, determine the velocity adjustment (Fig. 7b.1 $\sim$2.06% mark). The $c_{vel2,vel1,vel0}$, and $c_{ch1,ch0}$ are the $v_{LO}$ and $h_{ch}$ fit coefficients respectively, and account for estimation errors during the stance and subsequent rebound phases.

IMUPT-acceleration (Aerial Phases)

$=[0,0,1](R^{bw}F_{U1}/m-g^{w})$ during the aerial phase where the acceleration is inferred from the total commanded motor thrust, $F_{U1}=c_{m1}(\sum_{i=1}^{4}d_{mi})+c_{m0}$; where $c_{m1,m0}$ are the fit coefficients to convert from duty-cycle to thrust, $d_{mi}$ are the commanded motor duty-cycles, $m$ is the overall robot mass.

where, Appendix Robot Parameters shows the variables.

Based on the identified IMUPTs, two Kalman filters (KF1, KF2) and two error state Kalman filters (ESKF1, ESKF2) are developed (Table I). The vertical acceleration input $u_{k}=[0,0,1]\,a_{k}^{w}$ is composed of the world frame acceleration measurement $a_{k}^{w}=R_{k}^{bw}(a_{k}^{IMU}-b_{ak})+g^{w}$. These filters required both a process and measure noise covariance matrix containing up to four additional parameters including the standard deviation of the acceleration ($\sigma_{az}$), position ($\sigma_{pz}$), velocity ($\sigma_{vz}$), and bias ($\sigma_{bz}$) measurements. Additional adaptability is included in the state estimation, by first low pass filtering the accelerometer measurements resulting in two cutoff frequency parameters, $f_{HVSE}$ and $f_{HPE}$, for the HVSE and HPE, respectively. Finally, as two accelerometers are used, the switching g-force parameter, $g_{s}$, is included to determine at which point to switch from the low-g to high-g accelerometer measurements. The total Kalman filter parameters are presented in Table II; highlighting the challenge to hand tuning these filters.

IV-B HPE and HVSE Training

Initial experimental trials (n=998 hops using the true state for control) were conducted at hopping heights between 1 and 4 meters; where the acceleration readings from both the low-g and high-g accelerometers were captured (robot approx. 840 Hz) along with the true state data from a motion capture system (Vicon Vantage, 600 Hz). Implementing in simulation (MATLAB) a particular HPE and HVSE on the captured experimental data, allows for an exact recreation of what the robot would experience and the associated state estimation, without the risks and time associated with testing on the real robot. Therefore, assuming the initial experiments span the complete range of desired robot operational characteristics, this allows for a complete understanding of the HVSE error including all environmental (e.g. wind, substrate properties), locomotion (e.g., hopping heights, traveling speeds), control (e.g., controller driven behaviors), and sensor (e.g., noise, positioning) characteristics.

The ability to capture large amounts of experimental data and exactly simulate the HVSE errors, suggests that the Kalman filter parameters can be learned from the data. Therefore, due to the challenging parameter space with significant local minima, a genetic algorithm was developed to optimize the HPE and HVSE parameters listed in Table II. To control the hopping height of the robot, the most important parameter is the hop apex height, $z_{HA}$. Therefore, the cost function ($\mathcal{L}_{c}$) is set as the mean-absolute-percentage-error (MAPE, $\gamma_{1}$) in the apex height as,

$$\text{min}\,\mathcal{L}_{c}=\begin{cases}w_{1}\gamma_{1}&\text{if }n_{HA}=\hat{n}_{HA}\\ w_{2}\gamma_{2}+w_{3}\gamma_{3}&\text{otherwise}\end{cases}$$ (8)

$$\gamma_{1}=\frac{1}{n_{HA}}\sum_{i}^{n_{HA}}\bigg{|}\frac{z_{HA_{i}}-\hat{z}_{HA_{i}}}{\hat{z}_{HA_{i}}}\bigg{|}$$ (9)

$$\gamma_{2}=\sqrt{\frac{1}{n_{p}}\sum_{i}^{n_{p}}(z_{i}-\hat{z}_{i})^{2}},\,\gamma_{3}=\sqrt{\frac{1}{n_{v}}\sum_{i}^{n_{v}}(\dot{z}_{i}-\hat{\dot{z}}_{i})^{2}}$$ (10)

where, the ground true values [$\hat{z}_{HA},\hat{z},\hat{\dot{z}}$] of the hop apex height, $z_{HA}$, vertical position, $z$, and vertical velocity, $\dot{z}$, are as shown respectively. The RMSE in both vertical position, $\gamma_{2}$, and velocity, $\gamma_{3}$, are used to maintain a gradient in the event that the estimation does not capture the correct number of hops. The number of estimated hop apexes, $n_{HA}$, the true hop apexes, $\hat{n}_{HA}$, vertical heights, $n_{p}$, and vertical velocities, $n_{v}$, are represented accordingly. Finally, the weights are as follows, $[w_{1}=100,\,w_{2}=w_{3}=10]$.

The genetic algorithm is initialized as follows: individuals = 15 parameters (Table II), initial population = 1000 randomly generated individuals, maximum generations $n_{ga}$ = 20, elites carried over to the subsequent generation = 5%, crossover fraction = 80%, and the mutation fraction = 15%. Stochastic universal sampling is used to select two parents and the children are created through a uniform scattering crossover; where each optimization parameter is given a random one (parent-1) or zero (parent-2), determining which parent provides the parameter for the child. Finally, a constraint-aware adaptive mutation algorithm is used that creates the child as follows, $x_{ga}^{\prime}=x_{ga}+\alpha_{0}\,\exp(n_{ga_{i}}/n_{ga})\,d/||d||$; where the direction vector, $d_{(15\times 1)}$, is a random vector in which each component is drawn from a normal distribution $\sim\mathcal{N}(0,1)$. The $n_{ga_{i}}$ is the current generation number, and $\alpha_{0}$ is the initial step size that can be reduced if the constrains are violated. Finally, as stated previously, learning the parameters requires that the training data covers all desired robot operational characteristics. However, initial testing has shown that a much more limited set of data than that collected sufficiently well cover the operating range allowing for more efficient optimization on a reduced set containing 111 hops spread across the [1,2,3,4] m hop trials. This results in a reduction from multiple days per optimization cycle, to approximately 5 hours of parallel processing across eight cores; where code optimization could improve this further.

TABLE II: Optimization Parameters

Parm.	Description	Opt. Value	Bounds	Units
$f_{HVSE}$	Cutoff Freq. (HVSE)	7	[400,5]	Hz
$f_{HPE}$	Cutoff Freq. (HPE)	8	[400,5]	Hz
$g_{s}$	Accel. Switching g-force	14.24	[14.5,12]	g-force
$\sigma_{az}$	Accel. Meas. STD	9.9857	[10,0.0001]	m/s2
$\sigma_{bz}$	Accel. Bias Meas. STD	NA	[10,0.0001]	m/s2
$\sigma_{vz}$	Velocity Meas. STD	9.5722	[10,0.0001]	m/s
$\sigma_{pz}$	Position Meas. STD	0.0091	[10,0.0001]	m
$c_{vel2}$	Fit Constant Coef. vel2	-1.1246	[10,-10]
$c_{vel1}$	Fit Constant Coef. vel1	5.9203	[10,-10]
$c_{vel0}$	Fit Constant Coef. vel0	6.7054	[10,-10]
$c_{ch1}$	Fit Constant Coef. ch1	6.2138	[10,-10]
$c_{ch0}$	Fit Constant Coef. ch0	8.4355	[10,-10]
$c_{m0}$	Fit Constant Coef. $d_{mi}^{1}$	NA	[10,-10]
$c_{m1}$	Fit Constant Coef. $d_{mi}^{0}$	NA	[10,-10]

• Abbreviations: Measurement (Meas.), Acceleration (Accel.), standard deviation (STD), coefficient (Coef.).

IV-C Training Results

Optimizing across the four filters (KF1, KF2, ESKF1, ESKF2) for the MAPE in the apex height ($\gamma_{1}$) results in KF1 with the lowest error parameters (Table I Error column). The addition of the accelerometer bias measurement (IMUPT-Accel. Bias (Aerial Phase)) to KF2 results in increased errors. This suggests that the measurement adds more uncertainty than useful information; likely due to the delay between motor command and thrust production. Shifting to ESKF1 and ESKF2, IMUPT-Velocity(LO) was removed because of an observed instability suggesting that the short-term significant change in the error, cause by the aliasing at TD and LO, is not well managed by the ESKF.

Selecting KF1, the optimized parameter values can be seen in Table II, which result in a MAPE in the apex height of 10.77% (n=998 hops) over the complete set of training data. Analyzing the optimized values in more detail, we can further understand the system’s characteristics. First, the cutoff frequencies for both the HVSE and HPE show a trend towards low values, causing significant phase lag and smoothing of the acceleration measurements. This suggests a preference towards removing the inherently arbitrary aliasing of the acceleration measurements over phase accuracy for both the HVSE and HPE; where the colored dots in Fig. 8 show the HPE’s transitions during the subsequent experimental trials. Second, the liftoff velocity and commanded height show a correlation with the HVSE error as the exponential coefficients ($c_{vel2}$, $c_{vel1}$, $c_{ch1}$) are not equal to zero. Third, as expected, the smoothing does show an increase in the average error as evident by the non-unity constant coefficient values ($c_{vel0}$, $c_{ch0}$). Finally, arbitrary initialization of the estimated covariance matrix, $P$, yielded poor initial state estimations at the first TD. Therefore, to improve the initial state estimation and accelerate convergence of the KF, the $P$ is initialized to an average $P$ after convergence; yielding a much closer estimation of the real value for the first TD (MultiMo-MHR: initial $P=[0.0582,0.0774;0.0774,0.1441].*10^{-4}$).

Figure 4: Sensitivity to changes in the optimized parameters; $\mathcal{L}_{c}=100\gamma_{1}+10\gamma_{2}+10\gamma_{3}$ from equations 9, 10.

Figure 5: Change in the estimated ground height over multiple hops; where STD is the standard deviation, and the true ground height is zero.

Figure 4 shows the sensitivity of the optimized parameters to variation, where the cutoff frequencies ($f_{HVSE}$, $f_{HPE}$) and the switching g-force parameter ($g_{s}$) are the most sensitive; resulting in larger percent changes to the cost function than the percent change to the parameter. Furthermore, the significant change caused by increasing the $g_{s}$, suggests that the low-g accelerometer’s noise characteristics change abruptly as its limits are approached, and that the optimization methodology can learn the sensor characteristic from the data.

As can be seen in [47], zeroing altitude (IMUPT-position) removes drift but also the ability to track absolute height in environments where the ground is not flat; resulting in an error relative to the previous contact surface of approximately 0.05 m per step. Absolute height can be recovered through summing the change in height of each stepping or hopping cycle prior to the zeroing. The change in the ground height between cycles can be calculated in two ways including,

$$\Delta h_{1n} =(h_{HAn}-h_{TDn})-(h_{HAn-1}-h_{TDn-1})$$ (11)

$$\Delta h_{2n} =h_{TDn}$$ (12)

where the estimated height at the hop apex ($h_{HA}$) and TD ($h_{TD}$) are represented accordingly, and $n$ indicates the hop number. The difference between them ($h_{HA}-h_{TD}$) gives a measure of the drop height, and the difference between subsequent drop heights, provides a measure of the change in the ground height ($\Delta h_{1}$). Second, the estimated TD height ($h_{TD}$) provides a direct measure of the change in the ground height ($\Delta h_{2}$). Taking the average of these two measures,

$$h_{gn} =h_{gn-1}+(\Delta h_{1n}+\Delta h_{2n})/2$$ (13)

yields an estimate of the change in the ground height for the subsequent hop or step ($h_{g}$). Comparing the estimated ground height error of KF1 (Figure 5) to that of an untrained KF [47], shows that not only can drift be significantly reduced but that the absolute ground height can be accurately tracked. With an overall average ground height error of 0.018 m after 30 hops of up to 4 meters, the trained KF1 can achieve a 64% improvement on the single step error and significantly better than the sum over the same 30 steps. This is due to the optimization tending towards zero mean error with minimum standard deviation; resulting in a minimum error with an equal chance to be positive or negative for each hop.

Figure 6: Simulated State Estimation Error. a) Unbound and large bound position error, given state estimation at 10 to 400 Hz. b) Unbound and large bound velocity error, given state estimation at 10 to 400 Hz. c) Bound position error, given state estimation at 600 to 3360 Hz. b) Bound velocity error, given state estimation at 600 to 3360 Hz.

V Sensing Frequency

The previously developed dynamic model in [1] can be used to understand the impact of the sensing frequency on the HVSE; which fundamentally controls the aliasing of the accelerometer measurements. Specifically, the state estimation error boundedness and growth as a function of the limiting frequency; that being sensor read speeds or state estimation frequency. To explore the estimation error behavior, the measurement and estimation frequencies are varied over 12 values including [3360, 1680, 840, 700, 600, 500, 400, 300, 200, 100, 50, 10] Hz; where the MultiMo-MHR’s HVSE is optimized to, and operates at, 840 Hz. The simulation (MATLAB) results are shown in Fig. 6 with a control frequency of 400Hz, commanded height, $h_{ch}$, of 3 m, and the robot parameters in Appendix Robot Parameters. The results are shown for simulations both with and without sensor noise, and the controller is using the true state data, at the measurement frequency, as reference for control. The simulation without sensor noise (mean (w/o Noise)) can highlight effects of aliasing alone while the simulation with noise (mean and STD. (w/ Noise)) can provide an understanding of the true error growth characteristics.

As can be seen in Fig. 6a,b, frequencies from 10 to 100 Hz show an unbounded state estimation error with significant growth rate in the standard deviation as indicated by the large monotonic increase over time (Fig. 6a,b STD. (w/ Noise)); at $\sim$7 seconds both the 10 and 50 Hz trials are unable to sustain hopping and simply remain stationary on the ground. This unboundedness is similar to the drift associated with direct accelerometer integration (DAI) and therefore, provides no improvement in state estimation over DAI.

Frequencies from 200 to 400 Hz show a bound error characteristic in which the standard deviation in the error is not increasing monotonically (Fig. 6a,b STD. (w/ Noise)). However, both mean and standard deviation remain sufficiently large so as to cause variation in the subsequent hops (Fig. 6a,b mean and STD. (w/ Noise)).

Frequencies from 600 to 3360 Hz show both a bound mean error and a significant reduction in the standard deviation of the error, in both position and velocity, within each hop cycle (Fig. 6c,d mean and STD. (w/ Noise)); creating a deterministic starting point for the subsequent hop. Additionally, from Fig. 6c,d, we see both a low sensitivity to minor changes in frequency and diminishing returns as the frequency is increased. This suggests that the impact of accelerometer aliasing is saturating, and that the HVSE, trained at 840 Hz, is robust to frequency changes in the range of 600 to 3360 Hz. Interestingly, as seen in Fig. 6d, it is not required that the mean error in velocity be reduced; only that the standard deviation be reduced as the velocity scaling factor ($\delta_{vLO}$) accounts for mean errors in velocity.

VI Hopping Height Control Experiments

Hopping experiments were conducted with the MultiMo-MHR (Fig. 1a) at target heights of 1, 2, 3, and 4 m. Trials had two possible measurement inputs for vertical control; true state (GT, motion capture system, n = 70 hops) or the HVSE position and velocity (n = 121 hops). The HVSE is recorded in all scenarios to identify if the HVSE and vertical controller are decoupled and to evaluate the estimator performance.

The experimental procedure begins with zeroing of the accelerometers and motion capture system with the robot positioned vertically on the ground. The robot is then flown to the set initial height with a separate controller (only used in the startup procedure) before entering the first drop state, beginning the hop cycle and activating the hopping height controller. For a complete breakdown of the hop cycle and corresponding HPE and HVSE see Fig. 3; where, Fig. 1 and 7 provide a visual reference of the hop cycle and HVSE.

Figure 7: Example hopping trials using the hopping vertical state estimator (HVSE) based control a-b) Displays the vertical position and velocity error of the first 3 and last 2 hops for hopping trials at heights of 1,2,3&4 m. Each has a zoomed in region (a.1, b.1) showing the stance phase behavior. c) Snapshots from a high-speed video of the MultiMo-MHR hopping with the HVSE at a desired height of 3.25 m (Supplemental Video: Video1).

VI-A Results

Figure 8 shows the experimental results of the vertical position and velocity estimation using the HVSE. Table III compares the HVSE to other vertical position and velocity estimation techniques used in hopping robots across five error metrics (Appendix: Error Metrics) including: the mean of the normalized MAE in the vertical position (M1), the mean of the normalized MAE in the vertical velocity (M2), mean absolute percent error (MAPE) in the apex height (M3), mean absolute error (MAE) in the apex time (M4), and MAE in the desired apex height (M5). The table compares three vertical state estimation techniques to the HVSE. These include a Newtonian approach that assumes a ballistic aerial trajectory where the LO state is determined by a leg encoder (BA1) [4, 33, 2]. However, as rotor based hopping robots may not follow a ballistic trajectory, a natural extension is to dead-reckoning the accelerometer measurements (LO to TD) to improve the overall estimates (DR1). Finally, as shown in [49, 46, 47], a Kalman filter with a zero altitude update at TD can be employed to further improve the state estimation (KF3 = KF1 with only IMUPT-position(TD)). Four variations of the HVSE are shown including HVSE1, using the true state (GT) for control, and HVSE2, using the estimated state (SE) for control; where HVSE3 and HVSE4 only include the aerial phase errors. The controller for each technique can use either the true state (GT), focusing on the estimation characteristic, or the estimated state (SE), focusing on the overall hop height control performance.

Validating the simulation (Sim.) data, the experimental (Exp.) results show a similar reduction in both the error mean and standard deviation in position estimation at TD; as a result of the IMUPT-Position(TD) ((Sim.) Fig. 6c, (Exp.) Fig. 8a, a.1). The IMUPT-Velocity(MS) in both cases reduces the velocity mean error growth rate but with limited impact on the standard deviation ((Sim.) Fig. 6d, (Exp.) Fig. 8f). Finally, in both cases at LO (IMUPT-Position(LO), Velocity(LO)), a significant reduction in the velocity error’s mean and standard deviation is observed ((Sim.) Fig. 6d, (Exp.) Fig. 8g.1). The significant reductions in both the mean and standard deviations in the position and velocity errors creates a deterministic starting point for the subsequent hop cycle estimation; minimizing the error cause by the significant aliasing of the accelerometer measurements.

Table III highlights a number of key concepts. First, comparing the common hopping robot vertical state estimation technique (BA1) to its natural extension for rotor based hopping (DR1), suggests that the ballistic trajectory assumption may not be valid. Second, comparing DR1 to the Kalman filter based pedestrian navigation technique (KF3), suggests that a leg encoder can achieve similar results to that of a zeroing altitude condition at TD in a Kalman filter based approach. Third, comparing KF3 to the HVSE, suggests that a trained Kalman filter with the specified IMUPTs can improve the state estimation under significant aliasing conditions. Fourth, comparing HVSE1 to HVSE2 only shows a significant different in M5, suggesting that the controller and HVSE are independent; however, as expected M5 does show a significant difference, suggesting that the HVSE does introduce error beyond that caused by the controller itself. Fifth, comparing HVSE1 and HVSE2 to HVSE3 and HVSE4, suggests that the estimation of the stance phase states is more challenging; however, as the vast majority of the motion of the robot occurs in the aerial phase this will have a minimal overall impact. Sixth, from basic optics, the M1 metric is shown to be good approximation of the error in the horizontal state estimation, given that an optical flow sensor’s range finder (height measurement) is replaced with the HVSE (Appendix: Optical Flow Error). Seventh, given the MultiMo-MHR travels up to $\sim 8$ meters at velocities that range between $\sim\pm 7$ m/s in $\sim$2.6 seconds coupled with the significant and stochastic nature of aliasing, a certain amount of error is inherent without further direct measurements of the states. Finally, using experimental data, a KF or ESKF can be trained to extract as much information as possible from a highly aliased signal, and be incorporated back into a full state KF or ESKF.

Figure 8: Mean and standard deviation of the position (a-d) and velocity (e-h) error for the experimental trials (n=70 hops) utilizing the true state (motion capture data) for height control. a&e) Stance down phase, TD to MS. a&e.1) Zoomed in region of (a&e) to show the behavior after TD). b&f) Stance up phase, MS to LO. c&g) Rebound phase, LO to HA. c&g.1) Zoomed in region from 0-15% of Rebound phase c&g.2) Zoomed in region from 15-100% Rebound phase. d&h) Drop phase, HA to next TD. Note) The phases are segmented by the true state data, while the dots indicate the points at which the robot identified the specified transition.

TABLE III: HVSE Performance Metrics

Metric	Est.	Meas.	1-4m	Units	Statistics
(M1) Mean	BA1	GT	196.8	%	t(138)=30.97, p=0.00
nMAE Traj	DR1	GT	44.05	%	t(138)=5.62, p=0.00
Pos.	KF3	GT	60.45	%	t(138)=7.49, p=0.00
	HVSE1	GT	21.24	%	t(189)=0.99, p=0.32
	HVSE2	SE	19.46	%	Comparison
	HVSE3†	GT	20.21	%
	HVSE4†	SE	18.97	%
(M2) Mean	BA1	GT	179.9	%	t(138)=24.15, p=0.00
nMAE Traj	DR1	GT	28.15	%	t(138)=3.17, p=0.00
Velo.	KF3	GT	37.84	%	t(138)=5.24, p=0.00
	HVSE1	GT	19.58	%	t(189)=0.31, p=0.76
	HVSE2	SE	19.11	%	Comparison
	HVSE3†	GT	17.08	%
	HVSE4†	SE	16.50	%
(M3)	BA1	GT	53.85	%	t(138)=32.01, p=0.00
MAPE	DR1	GT	43.36	%	t(138)=5.90, p=0.00
Apex	KF3	GT	50.11	%	t(138)=7.38, p=0.00
Height	HVSE1	GT	13.77	%	t(189)=0.83, p=0.41
	HVSE2	SE	12.49	%	Comparison
(M4)	BA1	GT	0.7140	s	t(138)=16.83, p=0.00
MAE	DR1	GT	0.1323	s	t(138)=2.66, p=0.00
Apex Time	KF3	GT	0.1388	s	t(138)=4.18, p=0.00
	HVSE1	GT	0.0641	s	t(189)=0.96, p=0.34
	HVSE2	SE	0.0588	s	Comparison
(M5) MAE	HVSE1	GT	0.0912	m	t(170)=5.87, p=0.00
Apex Height	HVSE2	SE	0.2597	m	Comparison
(desired)

• Error metrics M1-M5 (n=191 hops, Appendix Error Metrics); where Meas. column indicates the measurement used for control (true state (GT), estimated state (SE)). ^* Values indicate the mean of the RMSE value for a complete hop cycle; TD to subsequent TD. ^† Data uses only the aerial phase as detected by the robot for analysis. Comparing the 1-4 m data to the HVSE2 using a two-tailed t-test given the null hypothesis with a 95% confidence interval. Related works, BA1 [4, 33, 2], KF3 [49, 46, 47].

VI-B Agility Metric

The commonly used hopping performance metric, vertical jumping agility ($\nu_{VJA}$), defined as an average climb speed in a gravitational environment, allows for quantifying and comparing jumping and hopping robot’s approximate average locomotion speed [42]. The metric is defined as follows, $\nu_{VJA}=\frac{h_{1}}{t_{apogee}}=\frac{h_{1}}{t_{s}+t_{r}}$ where, the apex height, $h_{1}$, stance through apex time, $t_{apogee}$, stance time, $t_{s}$, and rebound time, $t_{r}$, define the metric. However, being focused on climb rate, the $\nu_{VJA}$ is unable to account for the drop phase time of continuous hopping robots. A subsequent group amended this metric to the hopping agility metric ($\nu_{HA}$) adding the drop phase time, $t_{d}$, to quantify the overall average locomotion speed of continuous hoppers. The metric is defined as follows [2], $\nu_{HA}=\frac{2h_{1}}{t_{cycle}}=\frac{2h_{1}}{t_{s}+t_{r}+t_{d}}$ where, the total cycle time is $t_{cycle}=t_{s}+t_{r}+t_{d}$. However, the base metric, $\nu_{VJA}$, was developed to quantify existing systems’ performance from videos, and not as a robot design tool. Therefore, while increasing height and decreasing cycle time can increase both metrics, little can be said about control schemes, system losses, and non-constant hopping heights.

The similarity between the metrics allows for their unification, and a mathematical re-derivation allows for a deeper understanding of the effect of system design and operation parameters (Appendix: Hopping Agility). The unified hopping agility (UHA) metric is as follows:

$$\nu_{UHA} =\frac{h_{1}+\beta{h_{0}}}{t_{s}+t_{r}+\beta t_{d}}$$ (14)

$$=\frac{h_{1}+\beta{h_{0}}}{t_{s}+\sqrt{\frac{2h_{1}}{g(1-\gamma_{r}+\gamma_{lr})}}+\beta\sqrt{\frac{2h_{0}}{g(1-\gamma_{d}-\gamma_{ld})}}}$$ (15)

$$h_{0} =\zeta_{s}^{2}\frac{h_{1}(1-\gamma_{r}+\gamma_{lr})}{(1-\gamma_{d}-\gamma_{ld})}$$ (16)

where, the previous hop height, $h_{0}$, and gravitational acceleration, $g$, are represented accordingly. The average rebound phase acceleration, $a_{rebound}=-g(1-\gamma_{r}+\gamma_{lr})$, and the average drop phase acceleration, $a_{drop}=-g(1-\gamma_{d}-\gamma_{ld})$, are composed of both energy input, $[\gamma_{r},\gamma_{d}]<1$, representing average thrust forces in rebound and drop, and energy loss, $[\gamma_{lr},\gamma_{ld}]\geq 0$, representing average drag forces in rebound and drop; where positive $[\gamma_{r},\gamma_{d}]$ indicates upward thrust and $[\gamma_{lr},\gamma_{ld}]$ are always opposite the velocity vector. The change in stance phase velocity between TD and LO, $\zeta_{s}>0$, can capture energy lost, $0<\zeta_{s}<1$, no energy change, $\zeta_{s}=1$, and energy gained, $\zeta_{s}>1$ in the stance phase. The previous vertical jumping agility metric ($\nu_{VJA}$) can be derived from $\beta=0$ whereas, the hopping agility metric ($\nu_{HA}$) can be derived from $\beta=1$, thus unifying the two metrics within a mathematical framework to inform future hopping robot design. Comparing the MultiMo-MHR to existing hopping robots (Table IV) shows that it is able to outperform not only all existing robots but the Galago animal as well.

TABLE IV: Hopping Platforms Vertical Agility

Platform	Hop	$t_{apogee}$	$t_{cycle}$	$\nu_{VJA}$	$\nu_{HA}$
	Height(m)	(s)	(s)	(m/s) [42]	(m/s) [2]
Salto[42]	1.008	0.58	1.03	1.75	1.96
Salto-1P[30]	1.252	0.68	1.15	1.83	2.18
Galago[42]	1.742	0.78	1.55	2.24	2.24
PogoDrone[6]	0.7	4.18	8.25	0.17	0.17
PogoX[3]	0.6	0.44	0.88	1.37	1.37
HopCopter[2]	1.63	0.71	1.37	2.30	2.37
MultiMo-MHR	3.92	1.58	2.60	2.47	3.01
MultiMo-MHR	4.02	1.85	2.86	2.18	2.81

• Note: Vertical Jumping Agility ($\nu_{VJA}$) [42]. Hopping Agility ($\nu_{HA}$) [2]. : assumes drop time equals rise time given ballistic motion,: estimated by paper data,: assumes rise and drop time are half the total air time,: MultiMo-MHR max agility,: MultiMo-MHR max achieved height.

VII Summary

Here we have shown a light-weight, low computational load, IMU only, Kalman filter-based hopping vertical state estimator (position, velocity) coupled to a hopping phase estimator (drop, stance down, stance up, rebound) and the training procedure to learn the specific characteristics of the target system, including significant aliasing of the acceleration measurements. As discussed in the paper, the resulting filter can be reincorporated into a full state Kalman filter and error state Kalman filter. Finally, unlike previous works that trade absolute for relative height estimation to remove drift, the characteristics of the learned filter result in both an absolute height estimation and a significantly reduced drift (average ground height drift of 0.018 m after 30 hops of up to 4 m). The filter was experimentally validated on the MultiMo-MHR with hops up to 4.02 m, 2.47x greater than the next highest continuous hopping robot, at velocities approaching $\pm 7$ m/s. This results in a mean absolute percent error in hop apex height (M3) of 12.5% with an aerial trajectory average normalized mean absolute error in position (M1) and velocity (M2) of 19% and 16.5%, respectively.

Appendix

Hopping Agility

Expanding upon the Vertical Jumping Agility ($\nu_{VJA}$)[42] and Hopping Agility ($\nu_{HA}$)[2] metrics, we developed a unified agility metric (UHA: Unified Hopping Agility). Using plots from [42], the hopping cycle for the agility metric proceeds from the prior HA to current HA. While this is irrelevant for constant hopping height cycles, it can create a critical difference in agility for non-constant hopping height cycles.

Using fundamental Newtonian equations and starting with the sequence of prior HA to TD, the fundamental equations are,

$$a_{drop} =(-g+a_{d}+a_{ld})=-g(1-\gamma_{d}-\gamma_{ld})$$ (S1)

$$v_{TD} =v_{0}+a_{drop}t_{d}$$ (S2)

$$h_{TD} =h_{0}+v_{TD}t_{d}+\frac{1}{2}a_{drop}t_{d}^{2}$$ (S3)

where we have a control input in the drop phase of $a_{d}=g\gamma_{d}$ and accompanying loss shown as $a_{ld}=g\gamma_{ld}$. The system starts at an initial height and velocity ($h_{0},v_{h0}$) proceeds to a touchdown position and velocity of ($h_{TD},v_{TD}$) over a drop phase time window ($t_{d}$). In the stance phase, between TD and the subsequent LO, the total energy may change do to loss or addition,

$$\zeta_{s}=\frac{v_{LO}}{-v_{TD}}$$ (S4)

over a time period ($t_{s}$). The last segment is the rebound phase, between LO and the current HA, in which the equations are as follows,

$$a_{rebound} =(-g+a_{r}-a_{lr})=-g(1-\gamma_{r}+\gamma_{lr})$$ (S5)

$$v_{1} =v_{LO}+a_{rebound}t_{r}$$ (S6)

$$h_{1} =h_{LO}+v_{LO}t_{r}+\frac{1}{2}a_{rebound}t_{r}^{2}$$ (S7)

where we have a control input in the rebound phase of $a_{r}=g\gamma_{r}$ with accompanying loss shown as $a_{lr}=g\gamma_{lr}$. The system starts at LO with an initial height and velocity ($h_{LO},v_{LO}$) and continues to an ultimate height and velocity of ($h_{1},v_{1}$) over a rebound time window ($t_{r}$).

Solving for the drop phase time, $t_{d}$, from Eq. S3 and S2 and the rebound phase time, $t_{r}$, from Eq. S6 and S7, assuming $[v_{0},v_{1},h_{TD},h_{LO}]=0$ results in,

$$t_{d} =\sqrt{\frac{2h_{0}}{g(1-\gamma_{d}-\gamma_{ld})}}$$ (S8)

$$t_{r} =\sqrt{\frac{2h_{1}}{g(1-\gamma_{r}+\gamma_{lr})}}$$ (S9)

Using a switching parameters, $\beta$, the two agility metrics, $\nu_{VJA}$ and $\nu_{HA}$, can then be combined as,

$$\nu =\frac{h_{1}+\beta{h_{0}}}{t_{s}+t_{r}+\beta t_{d}}$$

$$=\frac{h_{1}+\beta{h_{0}}}{t_{s}+\sqrt{\frac{2h_{1}}{g(1-\gamma_{r}+\gamma_{lr})}}+\beta\sqrt{\frac{2h_{0}}{g(1-\gamma_{d}-\gamma_{ld})}}}$$ (S10)

where, $\beta=0$ gives the vertical jumping agility metric and $\beta=1$ gives the hopping agility metric, thus unifying the two metrics within a mathematical framework to inform future hopping robot design. To include changes in energy during the stance phase ($\zeta_{s}$), first solve eq. S2 and S6 for their time and substitute them into eq. S3 and S7, respectively. Solve each of these for their respective velocity at TD ($v_{TD}$) or LO ($v_{LO}$), substitute those into eq. S4, and solve for $h_{0}$ resulting in,

$$h_{0} =\zeta_{s}^{2}\frac{h_{1}(1-\gamma_{r}+\gamma_{lr})}{(1-\gamma_{d}-\gamma_{ld})}$$ (S11)

where, the previous hop height, $h_{0}$, and gravitational acceleration, $g$, are represented accordingly. The average rebound phase acceleration, $a_{rebound}=-g(1-\gamma_{r}+\gamma_{lr})$, and the average drop phase acceleration, $a_{drop}=-g(1-\gamma_{d}-\gamma_{ld})$, are composed of both energy input, $[\gamma_{r},\gamma_{d}]<1$, representing average thrust forces in rebound and drop, and energy loss, $[\gamma_{lr},\gamma_{ld}]\geq 0$, representing average drag forces in rebound and drop; where positive $[\gamma_{r},\gamma_{d}]$ indicates upward thrust and $[\gamma_{lr},\gamma_{ld}]$ are always opposite the velocity vector. The change in stance phase velocity between TD and LO, $\zeta_{s}>0$, can capture energy lost, $0<\zeta_{s}<1$, no energy change, $\zeta_{s}=1$, and energy gained, $\zeta_{s}>1$ in the stance phase. Comparing the MultiMo-MHR at differing hopping heights highlights how the hopping height affects the metrics (Table S1); where an increase is observed at lower heights and a saturation appears to occur at higher heights.

Analysis of the metric can highlight key characteristics for enhancing robot agility. First, given Eq. 15, assuming ballistic hopping at a constant height with no inputs or losses, and zero stance time ($[\gamma_{lr},\gamma_{ld}]=0,\zeta_{s}=1,[\gamma_{r},\gamma_{d}]=0,h_{1}=h_{0},t_{s}=0$), the resulting equation, $\nu=\sqrt{2gh_{1}}/2$, shows the metric will increase with increasing gravity or hopping height. Moreover, as the metric switching parameter, $\beta$, cancels out, the two metrics are equivalent under these conditions.

Second, given Eq. 15, it is clear that increasing the metric requires decreasing the denominator time parameters for the specified heights. The $t_{s}$ is inversely related to the robot’s natural frequency ($\sim\sqrt{m_{B}/K_{s}}$); where, the spring constant, $K_{s}$, and body mass, $m_{B}$ are represented accordingly. Therefore, the stance time ($t_{s}$) can be minimized by increasing $K_{s}$ and decrease $m_{B}$. The drop and rebound phase times ($t_{d},t_{r}$) can be minimized by maximizing the acceleration magnitudes ($a_{drop},a_{rebound}$); thus requiring maximum aerial forces. As the MultiMo-MHR can only generate thrust in the positive vertical direction, $0\leq[\gamma_{r},\gamma_{d}]<1$, the metric can be maximized by minimizing $[\gamma_{r},\gamma_{d},\gamma_{ld}]=0$, and maximizing the rebound loss $\gamma_{lr}$. This will inherently result in non-constant hopping height unless energy is added in stance; the amount of which can be determined by solving Eq. 16 for the velocity ratio between TD and LO, $\zeta_{s}$. This does highlight an additional important consideration of the metric. The metric can be increased with a higher previous hop, $h_{0}$, followed by a lower hop, $h_{1}$. However, this would reduce the hopping height with each cycle, and therefore could be used to increase agility for individual hops at the expense of long-term agility.

TABLE S1: MultiMo-MHR Vertical Agility

Meas.	Hop	$t_{s}$	$t_{apogee}$	$t_{cycle}$	$\nu_{VJA}$	$\nu_{HA}$
(n=hop num.)	Height(m)	(s)	(s)	(s)	(m/s) [42]	(m/s) [2]
GT (n=16)	1.03 $\pm$ 0.05	0.059 $\pm$ 0.001	0.68 $\pm$ 0.04	1.18 $\pm$ 0.04	1.52	1.74
HVSE (n=18)	1.00 $\pm$ 0.00	0.058 $\pm$ 0.002	0.69 $\pm$ 0.06	1.19 $\pm$ 0.08	1.45	1.68
GT (n=19)	2.02 $\pm$ 0.02	0.067 $\pm$ 0.001	0.98 $\pm$ 0.02	1.68 $\pm$ 0.01	2.05	2.41
HVSE (n=28)	1.97 $\pm$ 0.11	0.07 $\pm$ 0.001	0.96 $\pm$ 0.11	1.65 $\pm$ 0.13	2.04	2.39
GT (n=16)	3.03 $\pm$ 0.06	0.069 $\pm$ 0.001	1.26 $\pm$ 0.03	2.14 $\pm$ 0.04	2.40	2.84
HVSE (n=39)	2.85 $\pm$ 0.15	0.07 $\pm$ 0.001	1.29 $\pm$ 0.11	2.13 $\pm$ 0.13	2.20	2.68
GT (n=19)	3.78 $\pm$ 0.13	0.071 $\pm$ 0.001	1.71 $\pm$ 0.07	2.72 $\pm$ 0.08	2.21	2.78
HVSE (n=36)	3.41 $\pm$ 0.29	0.071 $\pm$ 0.001	1.64 $\pm$ 0.19	2.59 $\pm$ 0.24	2.08	2.63
GT (max)	3.92	0.071	1.58	2.60	2.47	3.01
HVSE (max)	4.02	0.072	1.85	2.86	2.18	2.81

• Note: Vertical Jumping Agility ($\nu_{VJA}$)[42]. Hopping Agility ($\nu_{HA}$) [2].

TABLE S2: Robot Parameters

Parm.	Description	Value	Units
$m_{B}$	Body Mass	0.5619	kg
$m_{L}$	Leg Mass	0.0981	kg
$I_{cmB}$	Body Rotational Inertia	0.0012	kg m2
$I_{cmL}$	Leg Rotational Inertia	0.0019	kg m2
$K_{s}$	Main Power Spring Const.	704	N/m
$K_{lb}$	Leg-Body Spring Const.	$400K_{s}$	N/m
$b_{lb}$	Leg-Body Damping Coef.	100	N s/m
$L_{1}$	Leg Top to CML	0.1053	m
$L_{2}$	Leg Bottom to CML	0.2821	m
$L_{3}$	Body CMB from Top Leg	0.1191	m
$d_{mtrx}$	Motor Offset from Body CM	0.0820	m
$d_{mtry}$	Motor Offset from Body CM	0	m

• Full schematic diagram shown in [1].

Robot Parameters

Table S2 shows the robot parameters used in the dynamic simulation from pervious work [1].

Error Metrics

Five error metrics are shown in Table S3 including:

$$M1 =1/n_{d}\sum_{j}^{n_{d}}\Bigg{(}\frac{\frac{1}{n_{p}}\sum_{i}^{n_{p}}\big{|}z_{i}-\hat{z}_{i}\big{|}}{\frac{1}{n_{p}}\sum_{i}^{n_{p}}\hat{z}_{i}}\Bigg{)}_{j}$$ (S12)

$$M2 =1/n_{d}\sum_{j}^{n_{d}}\Bigg{(}\frac{\frac{1}{n_{p}}\sum_{i}^{n_{p}}\big{|}\dot{z}_{i}-\hat{\dot{z}}_{i}\big{|}}{\frac{1}{n_{p}}\sum_{i}^{n_{p}}|\hat{\dot{z}}_{i}|}\Bigg{)}_{j}$$ (S13)

$$M3 =\frac{1}{n_{HA}}\sum_{i}^{n_{HA}}\bigg{|}\frac{z_{HA_{i}}-\hat{z}_{HA_{i}}}{\hat{z}_{HA_{i}}}\bigg{|}$$ (S14)

$$M4 =\frac{1}{n_{HA}}\sum_{i}^{n_{HA}}\big{|}t_{HA_{i}}-\hat{t}_{HA_{i}}\big{|}$$ (S15)

$$M5 =\frac{1}{n_{HA}}\sum_{i}^{n_{HA}}\big{|}z_{HA_{i}}-\bar{z}_{HA_{i}}\big{|}$$ (S16)

where, $n_{d}$ is the number of hops, $n_{p}$ is the number of points per hop (TD to subsequent TD), $n_{HA}$ is the number of hop apexes, and $[z,\hat{z}],[\dot{z},\hat{\dot{z}}],[z_{HA},\hat{z}_{HA}],[t_{HA},\hat{t}_{HA}]$ represent the estimated and true vertical position, velocity, hop apex height, and hop apex time, accordingly. Finally, $\bar{z}_{HA}$ represents the desired hop apex height.

TABLE S3: HVSE Performance Metrics

			Set Height (m)
Metric	Est.	Meas.	1	2	3	4	1-4m	Units	Comparison to
Description	Type								HVSE
(M1) Mean	BA1	GT	171.9	157.7	193.5	248.8	196.8	%	t(138) = 30.97, p = 0.00
nMAE Traj Pos.	DR1	GT	39.10	30.06	29.10	74.01	44.05	%	t(138) = 5.62, p = 0.00
	KF3	GT	50.02	44.12	51.64	89.64	60.45	%	t(138) = 7.49, p = 0.00
	HVSE1	GT	24.64	20.29	13.81	25.60	21.24	%	t(189) = 0.99, p = 0.32
	HVSE2	SE	25.18	22.48	13.88	19.79	19.46	%	Statistical Comparison
	HVSE3†	GT	21.36	19.71	13.51	25.38	20.21	%
	HVSE4†	SE	25.17	21.91	13.48	19.52	18.97	%
(M2) Mean	BA1	GT	143.1	136.5	174.7	245.0	179.9	%	t(138) = 24.15, p = 0.00
nMAE Traj Velo.	DR1	GT	24.24	19.14	17.97	48.46	28.15	%	t(138) = 3.17, p = 0.00
	KF3	GT	30.91	27.19	31.70	57.35	37.84	%	t(138) = 5.24, p = 0.00
	HVSE1	GT	28.52	21.52	12.44	16.13	19.58	%	t(189) = 0.31, p = 0.76
	HVSE2	SE	25.53	21.95	15.73	17.36	19.11	%	Statistical Comparison
	HVSE3†	GT	24.88	18.78	11.37	13.63	17.08	%
	HVSE4†	SE	21.98	19.59	13.35	14.78	16.50	%
(M3)	BA1	GT	57.34	51.72	54.53	53.23	53.85	%	t(138) = 32.01, p = 0.00
MAPE	DR1	GT	29.04	22.13	22.47	73.30	43.36	%	t(138) = 5.90, p = 0.00
Apex Height	KF3	GT	37.28	33.24	40.75	81.82	50.11	%	t(138) = 7.38, p = 0.00
	HVSE1	GT	13.50	12.25	8.81	19.68	13.77	%	t(189) = 0.83, p = 0.41
	HVSE2	SE	14.79	14.36	8.11	14.62	12.49	%	Statistical Comparison
(M4)	BA1	GT	0.3305	0.4653	0.6987	0.9950	0.7140	s	t(138) = 16.83, p = 0.00
MAE	DR1	GT	0.0554	0.0583	0.0659	0.2430	0.1323	s	t(138) = 2.66, p = 0.00
Apex Time	KF3	GT	0.0726	0.0830	0.1178	0.2483	0.1388	s	t(138) = 4.17, p = 0.00
	HVSE1	GT	0.0611	0.0581	0.0369	0.0956	0.0641	s	t(189) = 0.96, p = 0.34
	HVSE2	SE	0.0524	0.0611	0.0505	0.0695	0.0588	s	Statistical Comparison
(M5) MAE Apex	HVSE1	GT	0.0477	0.0294	0.0562	0.2190	0.0912	m	t(170) = 5.87, p = 0.00
Height (desired)	HVSE2	SE	0.0353	0.0959	0.1773	0.5885	0.2597	m	Statistical Comparison

• Error metrics M1-M5 (n=191 hops, Appendix Error Metrics); where Meas. column indicates the measurement used for control (true state (GT), estimated state (SE)). ^* Values indicate the mean of the RMSE value for a complete hop cycle; TD to subsequent TD. ^† Data uses only the aerial phase as detected by the robot for analysis. Comparing the 1-4 m data to the HVSE2 using a two-tailed t-test given the null hypothesis with a 95% confidence interval. Related works, BA1 [4, 33, 2], KF3 [49, 46, 47].

Optical Flow Error

To determine the potential error in the horizontal distance traveled, measured by optical flow. From basic optics, assuming the sensor is traveling parallel to a surface in a single direction, we know the the ratio between the distance an object moves on a camera’s image plane, $\hat{d}_{c}$, divided by the focal distance, $\hat{h}_{c}$, is equal to the distance the object moves in the world, $d_{w}$, divided by the distance to the surface, $h_{w}$, as follows $d_{w}/h_{w}=\hat{d}_{c}/\hat{h}_{c}$. If the measurement of the distance to the object includes error, $h_{w}=\hat{h}_{w}+\delta h_{w}$, then, so too will the distance traveled, $d_{w}=\hat{d}_{w}+\delta d_{w}$. Substituting this into the previous equation results in,

$$\hat{d}_{w}+\delta d_{w}=(\hat{d}_{c}/\hat{h}_{c})\hat{h}_{w}(1+\alpha_{h})$$ (S17)

where, the $\hat{}$ represents the true value, $\delta$ represents the error, and the height error is converted to a percent of the true height, $\delta h_{w}=\alpha_{h}\hat{h}_{w}$. The $\alpha_{h}$ can be approximated as the mean of the normalized mean-absolute-error (nMAE) in the height measurement (Table S3, M1); defined as the mean of the mean-absolute-error in the height divided by the mean of the true height per hop cycle. The error as a percent of the distance traveled can therefore be determined from the ratio of the horizontal distance both with and without the error as,

$$\frac{\hat{d}_{w}+\delta d_{w}}{\hat{d}_{w}}=1+\alpha_{h}\rightarrow\alpha_{h}=\frac{\delta d}{d_{w}}$$ (S18)

where, it can be seen that $\alpha_{h}$ also represents the horizontal error, as a percent of the distance traveled, that would be incurred by replacing an optical flow sensor’s range finder with the HVSE.