Bio-inspired Thrust-Assisted Trotting on a legged-Aerial Robot via Real-Time QP-MPC

Introduction

Birds are capable of using their wings to fly and also use their feet to walk on the ground, perch, and grabbing objects. Additionally, birds such as chukar’s wing-assisted incline running [1] and widowbirds’ controlled flapping descent [2] can use their wings while walking over a steep slope or during a rapid descent to stabilize their upper body orientation and velocity. We aim to take inspiration from and mimic the biological counterparts by combining both ground and aerial locomotion modes to achieve a versatile legged-aerial multi-modal robot.

There are several studies focused on multi-modal locomotion, specifically ground-aerial robots showing their versatility in traversing various terrains. Some of these studies focus more on transitioning between locomotion modes instead of using all available actuators at the same time. Some robot design requires transformation to switch mode of locomotion, such as [3, 4], which

transitions between terrestrial and aerial locomotion by morphing its appendages, [5] is a quadrupedal robot with propellers attached to its knee-axis and can perform flight after transforming into a quadrotor configuration. On the other hand, some robot design with a fixed configuration such that they can perform ground and aerial locomotion without transforming, such as a jet powered bipedal robot capable performing large vertical jumps [6], pogo-stick mechanism with quadrotor dynamics [7], a bipedal robot [8] and a quadrupedal robot [9] with quadrotor style propulsion attached to its torso. A multi-modal robot can utilize both ground and aerial actuators at the same time if the design allows it, making it possible to augment a mode of locomotion or perform a gait that can’t be done by using either ground or aerial actuator by themselves. This work explores the capability of our legged-aerial Harpy platform, shown in Figure 1, in performing a trotting gait using both the legs and thrusters at the same time. Bipedal robot locomotion requires a regulation of ground reaction force (GRF) so that the foot does not slip. The addition of thruster forces during walking will impact the resulting GRF and must be regulated together with the walking gait. Simulation studies have laid a theoretical groundwork for using both legs and thrusters at the same time to augment the legged locomotion to perform difficult maneuvers, such as narrow path walking [10, 11] and steep incline walking [9, 12, 13, 14, 15]. Building on this foundation, the Harpy platform has also been validated in simulation across various locomotion tasks, such as jumping [16, 17], walking [18, 19, 20, 21], and inclined slope walking [22]. However, unlike in simulations, the practical usage of the algorithm in a real-time controller requires an approach that can be solved quickly and efficiently.

A common solution for GRF regulation for legged robot uses an optimization method, such as quadratic programming (QP) or nonlinear model predictive control (MPC), by modeling the robot’s dynamics and implementing the foot no-slip condition (the friction cone constraint) as constraints in the optimization. The optimization problem needs to be solved quickly to be feasible for a real-time implementation, where a reduced order model (ROM) is often used to reduce the model complexity and computational time.

The reduced-order modeling tool, such as the Linear inverted pendulum (LIP) and the spring-loaded inverted pendulum (SLIP) model, captures the essential dynamics of walking or running [23, 24]. The SLIP model, first formalized by [23] and validated through extensive biomechanical studies [24, 25], has emerged as the canonical template for understanding and controlling such dynamic gaits. The SLIP model’s elegance lies in its simplicity: a point mass connected to a massless spring leg model, where the ground contact forces is applied alternatively between stance and flight phases. The reduced-order representation has proven remarkably effective in describing the locomotion of diverse running animals [26] and has guided the development of numerous robotic systems, from Raibert’s pioneering hopper [27] to modern dynamic bipeds [28, 29, 30, 31, 32, 33, 34, 35, 36].

In this paper, we formulated an MPC optimization problem into a QP, and implemented it in our Harpy robot to satisfy ground constraints to perform a thruster-assisted trotting gait. The QP is formulated using a Thruster-SLIP (TSLIP) model, where thrusters are incorporated in the SLIP reduced order model, which then is implemented on hardware and solved in real-time using the qpSWIFT [37] algorithm at 500 Hz control loop frequency. This controller is then tested on hardware during the thruster-assisted trotting gait, where the experimental results showed the controller’s feasibility and performance. The contributions of this paper is listed as follows: 1) the derivation of the dynamic Thruster-SLIP model of a legged-aerial robot, 2) implementation of the QP-SLIP model on hardware, and 3) experimental validation of the controller performance showing the satisfaction of the GRF constraints in real-time. This paper is organized as follows: the robot design overview, followed by dynamic modeling, controller design and formulation, leg trajectory controller. Then, we present the experimental results, discussions, concluding remarks, and future work.

Robot Overview

Harpy [38, 39, 40] is designed with three joints for each leg: hip frontal joint, hip sagittal joint, and knee joint. The lower leg forms a parallelogram linkage, where a shock absorber spring is attached to the rear lower linkage to reduce the impulsive force on the leg touchdown. There are four electric ducted fans (EDFs) attached to the torso: a high-powered EDF at each side of the robot, and two smaller upper EDFs with front/rear-facing orientation. Each larger and smaller EDF can generate a maximum thrust of 5kgf and 0.9kgf, respectively. The total weight of the robot is 7.1kg. Each joint is made with 3D printed material with embedded components. Joints are further connected using oval shape 2mm thickness carbon-fiber tube. The robot is equipped with STM32 F429ZI, which acts as a communication hub between ICM-20948 (IMU), OptiTrack data received over Wi-Fi using ESP32, and the xPC target machine.

Dynamic Modeling

Consider the thruster-assisted bipedal robot shown in Figure 2. The high fidelity model of the system can be developed using Lagrangian mechanics on the unpinned model. The Lagrangian is defined as $\mathcal{L} = \mathcal{K} - \mathcal{V}$, where the kinetic energy ($\mathcal{K}$) and potential energy ($\mathcal{V}$) are given by

where $\mathcal{F} = \{B, F_L, F_R, H_L, H_R, K_L, K_R, A_L, A_R\}$ denotes the set of rigid bodies and links comprising of torso ($B$), frontal joints ($F_L, F_R$), hips ($H_L, H_R$), knees ($K_L, K_R$), and ankles ($A_L, A_R$). Here, ($L$) and ($R$) denote the left and right leg, respectively. For each body $i \in\mathcal{F}$, $m_i$ denotes the mass, $I_i \in\mathbb{R}^{3 \times3}$ is the rotational inertia tensor expressed in the body frame, $p_i \in\mathbb{R}^{3 \times1}$ is the center-of-mass position for each link and body, $\omega_i \in\mathbb{R}^{3 \times1}$ is the angular velocity, and $g$ is the gravitational acceleration constant.

The dynamics can be derived by applying Hamilton’s principle to the Lagrangian as:

where the generalized coordinates are defined as:

with $n_q = 12$, $p_B \in\mathbb{R}^{3}$ denoting the torso position, $\Theta_B \in\mathbb{R}^{3}$ representing the body orientation parameterized using Euler angles in the XYZ (roll-pitch-yaw) convention, and $q_j \in\mathbb{R}^{6}$ containing the six joint angles. Further, the terms in Eq. (1) are defined where $M(q) \in\mathbb{R}^{n_q \times n_q}$ is the symmetric positive-definite mass-inertia matrix, $h(q,\dot{q}) = C(q,\dot{q})\dot{q} + G(q)$ is the vector combining Coriolis/centrifugal forces and gravitational forces, $B_m \in\mathbb{R}^{n_q \times6}$ is the joint torque selection matrix, $u_m \in\mathbb{R}^{6}$ is the vector of joint torques, $u_t \in\mathbb{R}^{n_t}$ is the vector of thruster forces with $n_t = 4$ and $\lambda_c \in\mathbb{R}^{n_c}$ is the ground reaction forces (GRF) at the contact point with $n_c \leq6$ depending on contact configuration.

The contact Jacobian $J_c(q) \in\mathbb{R}^{n_c \times n_q}$ maps generalized velocities to contact point velocities and is constructed as:

with $p_{F_L}, p_{F_R} \in\mathbb{R}^{3}$ denoting the positions of the left and right foot contact points, respectively. The thruster Jacobian $B_t(q) \in\mathbb{R}^{n_q \times n_t}$ maps thruster forces to generalized forces:

where $J_{T_i}(q) = \frac{\partial p_{T_i}}{\partial q}$ for each thruster mounting point $T_i \in\{T_L, T_R, T_F, T_B\}$, representing the left, right, front-facing, and back-facing thrusters. $B_t(q)$ is configuration-dependent since it is applied in the body-fixed frame and its effective contribution to the generalized coordinate varies with the body’s orientation.

Decoupled controller formulation

The overall control architecture is illustrated in Figure 3. The proposed framework adopts a decoupled structure separating the sagittal locomotion problem, handled by the TSLIP-based QP-MPC, from the frontal-plane stabilization problem, handled by the orientation controller.

ROM dynamics

The ROM considers the mass ($m_b$) concentrated at single point defined by $p_c = [p_{c,x}, p_{c,z}]^{\mathsf{T}}$ and ignores the orientation of the robot. Two different forces are acting on the point mass,

where, $u_t$ is the thrust force acting at the COM position, and $u_g$ is the GRF acting at the point of contact $p_f = [p_{f,x},p_{f,z}]^{\mathsf{T}}$. Under these forces, the ROM dynamics are described as follows

where, state vector is defined as $x = [p_c^{\mathsf{T},\dot{p}_c^{\mathsf{T}}}]^{\mathsf{T}} \in\mathbb{R}^{4}$, input vector is defined as $u = [u_g^{\mathsf{T}},u_t^{\mathsf{T}}]^{\mathsf{T}} \in\mathbb{R}^{4}$ and the system matrices follows directly from Newton’s second law applied to the point mass:

Applying the forward Euler method with step size $h$ gives

where $i$ is the current time step, the terms of the discrete dynamics are defined as $A_d = (I + Ah)$, $B_d = Bh$, and $d_d = [0,0,0,-hg]^{\mathsf{T}}$.

SLIP dynamics and constraint

In the stance phase, the GRF is generated by a spring-loaded virtual leg connecting the COM to the foot contact point, which is held fixed by friction. The leg vector and length are

With rest length $l_0$ and spring stiffness $k_s$, the spring force magnitude obeys a unilateral Hooke’s law:

where $(\cdot)_+=\max(\cdot,0)$ enforces that the spring can only push, not pull. The corresponding GRF vector directed along the spring axis is

Since (11) is nonlinear in the state, we linearize it about a nominal trajectory $(\bar{\mathbf{x}},\bar{\mathbf{u}})$ via a first-order Taylor expansion:

This affine approximation can be compactly written as the equality constraint

QP-MPC formulation

The SLIP equality constraint (13) is embedded within a receding-horizon MPC framework that optimizes the states and the control inputs. Given a reference trajectory $\{\mathbf{x}_{r,k}\}_{k=0}^{N}$ over prediction horizon $N$, positive-definite weight matrices $Q$ and $R$, and the initial condition $\mathbf{x}_{\mathrm{init}}$, the MPC problem is

The input bounds enforce physical limits on the GRF and thruster forces. The friction cone requires

and the saturation limits for the thruster and GRF are

These are assembled into the block-diagonal inequality matrices $A_{\mathrm{ineq}}$ and $\mathbf{b}_{\mathrm{ineq}}$ over the full horizon. To reduce the decision variables and improve solver speed, we eliminate the state from the decision variables through recursive forward simulation of (8). Stacking states and inputs column-wise as

The rollout over the horizon yields the affine map

where,

Substituting (16) into the SLIP equality constraints (13) eliminates $X$, yielding a standard QP in $U\in\mathbb{R}^{Nn_u}$:

with cost matrices

where $\bar{Q}=\operatorname{diag}([Q,\cdots,Q])\in\mathbb{R}^{(N+1)n_x\times(N+1)n_x}$, $\bar{R}=\operatorname{diag}([R,\cdots,R])\in\mathbb{R}^{Nn_u\times Nn_u}$ and $X_r=[\mathbf{x}_{r,0}^\top,\cdots,\mathbf{x}_{r,N}^\top]^\top$ is the stacked reference trajectory over the horizon.

Roll Orientation Controller

While the sagittal dynamics are handled by the QP-MPC described above, roll stabilization during single-support is decoupled and delegated to the side thrusters. During the single-support phase, since the robot doesn’t have any ankle roll joint, the side thrusters provide the necessary actuation authority in the frontal plane.

A cascaded control architecture is employed. The outer loop is a PD controller that maps the COM lateral position error and velocity to a reference roll angle:

where $e_y=y-y_{\mathrm{in}}$ is the error in the lateral position from the initial position at the start of stance phase and $\dot{e}_y$ is the velocity error, and $K_{p,y}, K_{d,y} > 0$ are tuned gains. The inner loop is a PID controller that tracks $\phi_{\mathrm{ref}}$ by commanding side thruster force:

where $i = [L,R]$ depends on leg that is currently on swing cycle, $\delta t = 0.002$ s is the control loop rate, and $K_{p,\phi}, K_{d,\phi}, K_{I,\phi} > 0$ are the gains. The feedforward term $u_{\mathrm{t,ff}}$ applies a bias force to pre-compensate for the lateral moment that is induced by the swing leg movement. The roll angle $\phi$ is given by the onboard EKF, and $\dot{\phi}$ is given by the onboard complementary filter.

Stance leg trajectory tracking controller

The stance leg controller translates the QP-MPC output into joint-space commands for the stance leg. At each control cycle, the optimal output $u^*$ from QP are used to propagate the discrete dynamics (8) by one step, yielding a desired COM state:

The desired COM position $p^{\mathrm{des}}_c$ is extracted from $x^{\mathrm{des}}_{k+1}$. The corresponding desired foot position in the body frame is computed as follows

where $p_f$ is the current stance foot position in the inertial frame and $R_b \in SO(3)$ is the rotation matrix from the inertial to the body frame. The desired foot is then mapped to joint angles via inverse kinematics of the stance leg.

Result

Experimental Details

Two hardware experiments were conducted: Experiment A with a forward velocity reference of 0 m/s over 15 steps and Experiment B with a velocity reference of 0.1 m/s over 10 steps. All experiments were conducted using the xPC target machine running MATLAB 2020a Simulink Real-Time at a 500 Hz loop frequency. The thruster PWM-to-force mapping was characterized before the experiments and is summarized in Table 1. The SLIP parameters are described as follows $l_0 = 0.55$ m, $k_s = 1380$ N/m, and $m_b = 7.1$ kg. The spring stiffness was selected such that a nominal compression of 0.05 m produces a GRF approximately equal to body weight, ensuring the SLIP model operates within a physically meaningful force limit during stance. The QP-MPC runs at a rate of 500 Hz with $N = 5$ horizon. The friction coefficient was set as $\mu= 0.5$. The lower and upper bounds for decision variables are $[-8.7,10,-75,0]^{\mathsf{T}}$ and $[8.7,10,75,150]^{\mathsf{T}}$ respectively. The cost matrices for states and inputs are defined as $Q = \operatorname{diag}([100,40,5000,1000])$ and $R = \operatorname{diag}([0.05,0.3,1,1])$. Position tracking cost is intentional relaxed to accommodate the limit forward thrust authority.

PWM (%)	Side thruster (kgf)	Upper thruster (kgf)
10	0.092	0.0875
20	0.37	0.175
30	0.809	0.2625
40	1.346	0.35
50	1.97	0.4375
60	2.674	0.5250
70	3.33	0.6125
80	3.91	0.7
90	4.49	0.78
100	5.25	0.87

Experimental Results & Discussion

Figure 4 presents the stable trotting results snapshot from Experiment B. The body position, velocity, and orientation are calculated using the onboard EKF estimator that combined IMU data at 500 Hz and OptiTrack data at 240 Hz. The thrust and stance leg was controlled using the proposed controller, while the swing leg was driven to the desired location using Raibert’s heuristic foot placement controller.

Figure 5 illustrates the computation performance of the QP solver on hardware. The total computation time for both experiments, which includes the setup-time and solve-time, consistently remain below 0.2ms, well within the 2ms control loop. Additionally, the solver converges in max 10 iteration, confirming the QP formulation remains well-conditioned throughout the locomotion.

Figure 6 shows the QP-commanded GRF distribution within friction cone bounds of $\mu= 0.5$ for both experiments. The QP controller was successfully able to modulate the thrust force while enforcing the no-slip contact constraint across both experiments. The time-colored trajectory confirms a consistent GRF regulation with no degradation over the course of the experiments.

Figure 7 presents the thruster force commanded during both experiments. The QP-commanded forward thrust frequently approaches the hardware saturation limit of $\pm8.7$N, indicating that the optimizer consistently requests forces near the physical capability of the thruster system. This is attributed to Thrust-to-weight contribution of the forward thruster, which produce a max of 8.7N compared to the robot’s 70 N body weight. Despite operating near saturation for extended periods, the QP-TSLIP controller maintains stable locomotion throughout the experiments as evidenced by the stable orientation response.

The vertical thruster force is included as an optimization variable within the QP but is bounded to a constant value of 10N which is equally distributed to both the side thrusters. This design choice is motivated by safety considerations specific to thruster-assisted biped locomotion - aggressive variation in vertical thrust risk inducing unintended floating behavior, since high $u_{t,z}$ combined with $u_{t,r}$ could cause sudden roll disturbance on the swing side. By constraining $u_{t,z}$ to a fixed gravity compensation value within the optimizer, the QP retains a consistent formulation while ensuring predictable ground contact dynamics throughout locomotion. The left and right side thruster forces combine the QP vertical force baseline of 5N per thruster with PID orientation correction. The side thruster oscillates dynamically between 5 − 25N as the orientation controller actively rejects roll disturbances throughout locomotion.

Figure 8 presents the COM body state trajectory for both experiments. The forward velocity reference is constant at 0 m/s and 0.1 m/s for Experiment A and Experiment B respectively. The instantaneous velocity oscillates due to thrust force tracking delay, and impact forces of the swing-to-stance transitions. The mean forward velocity remains close to its respective reference. The forward position reference is re-initialized at each touchdown to the current robot position and marched forward at the commanded velocity, preventing position error accumulation across steps. The pitch and roll angles remain bounded with ±9° and ±4° respectively, confirming stable sagittal and frontal-plane dynamics under the proposed decoupled control architecture. Figure 9 presents the joint tracking errors for all six joints. Joint tracking errors remain bounded within ±1° across both experiments.

Figure 10 validates the linearized SLIP model used within the QP framework. At each control cycle, the current robot state is used as an initial condition for a 5-step open-loop rollout using the QP-optimal inputs. The prediction rollout qualitatively captures the trend of the hardware trajectory. Since jerk is not included as a state variable in the ROM, sudden velocity changes arising from contact transition and impact dynamics which are not captured within the prediction model. As a

result, brief deviations between the prediction and actual trajectories are expected at the swing-to-stance phase transition, where impulsive GRF produces rapid velocity changes that lie outside the modeling assumptions of the linearized SLIP dynamics. The frequent re-initialization inherent to the receding horizon strategy effectively prevents error accumulation.

Conclusions

This paper presented the QP-TSLIP controller for thruster-assisted bipedal locomotion on the hardware platform. The proposed decoupled architecture handles sagittal locomotion via a condensed QP-MPC with SLIP constraint, while frontal plane stability is maintained by a cascaded controller. Hardware experiments across 25 steps demonstrate real-time QP solving with 0.18 ms at 500 Hz, consistent GRF cone satisfaction, and stable pitch and roll within ±9° and ±4° respectively.

Future work will focus on three directions: extending the formulation to an actuated SLIP model enabling active spring length modulation for more aggressive locomotion; incorporating explicit thruster tracking delay into the QP constraint to further reduce the velocity tracking error, and relaxing the constant vertical thrust bound to enable adaptive z thrust for more aggressive locomotion, such as 1 − 2m high jumps.