1 Introduction↩︎

Bipedal locomotion in quadrupedal robots provides a means to extend their capabilities in environments where space is limited and mobility demands exceed the scope of traditional quadrupedal gaits. Operating in an upright posture allows quadrupeds to reach objects positioned at a higher height, negotiate obstacles that would obstruct four-legged motion, and perform manipulation tasks that require elevation of the torso. However, achieving fully dynamic bipedal balance in such constrained environments is especially challenging due to limited foothold options and the morphological design of quadrupeds. To address these challenges, we investigate wall-supported bipedal locomotion, in which the robot leverages environmental contacts with vertical surfaces for additional stability. This strategy enables robots to operate effectively in tight corridors, cluttered industrial settings, and disaster sites, where leaning against walls or bracing for balance is often essential. For example, a quadrupedal robot navigating a narrow hallway can transition into a wall-supported bipedal posture to reach and operate a control panel that would otherwise be inaccessible in a purely quadrupedal stance. By exploiting walls as supportive contacts, quadrupeds gain a robust and practical mechanism for upright locomotion in constrained environments. At the same time, this strategy introduces a new set of challenges in modeling, control, and real-time execution, which we address in this work.

Despite its potential, wall-supported bipedal locomotion introduces several challenges. First, the dynamics of quadrupedal robots in an upright posture are inherently hybrid, unstable, and high-dimensional, involving fast-varying reaction forces between the limbs and the environment together with slower variations in foot placement and contact transitions for both upper and lower limbs. This multi-rate structure complicates the design of real-time predictive controllers for force and footstep planning, since standard nonlinear model predictive control (NMPC) formulations typically assume uniform input rates. Second, operating in constrained environments requires tight coupling between locomotion and environmental support forces, demanding simultaneous regulation of the center of mass (CoM) trajectory, torso orientation, and wall interaction forces while respecting feasibility constraints such as torque limits, kinematic reachability, and friction cones. Finally, these requirements must be satisfied under stringent real-time constraints, necessitating optimization-based control frameworks that can reliably perform predictive planning and whole-body coordination at high frequency despite significant computational complexity.

The overarching goal of this paper is to develop a unified and computationally efficient multi-rate nonlinear model predictive control (MR-NMPC) framework that simultaneously plans interaction forces and footstep placements, enabling quadrupedal robots to achieve robustly stable wall-supported bipedal locomotion in constrained environments.

1.1 Related Work↩︎

Model predictive control (MPC) has emerged as a powerful framework for trajectory optimization and feedback control in legged robots [1]. A common strategy is to pair MPC with reduced-order, or template, models [2], which capture the essential features of locomotion while abstracting away the full complexity of high-dimensional robot dynamics. Among the most influential examples is the linear inverted pendulum (LIP) model [3], together with its numerous extensions, including the angular momentum LIP [4], the spring-loaded inverted pendulum (SLIP) [5], variable-height inverted pendulum model [6], the vertical LIP [7], and the hybrid LIP [8]. Beyond pendulum-based abstractions, alternative reduced-order formulations such as centroidal dynamics [9] and the single rigid body (SRB) model [10]–[14] have also been widely adopted, offering computationally tractable yet more physically representative descriptions of legged locomotion. While quadratic programming (QP)-based MPC formulations for linearized template models are highly efficient, they often struggle to capture the nonlinear effects inherent in models such as centroidal dynamics and SRB, limiting their ability to fully exploit the dynamics of legged systems. To overcome these limitations, NMPC has been employed for real-time trajectory planning of dynamic gaits and for safe motion in cluttered environments, see, e.g., [1], [15]–[21]. In addition, several efficient methods have been proposed for whole-body NMPC, such as [22]–[24], typically formulated as a nonlinear program (NLP). However, these approaches have been developed primarily for quadrupedal gaits in quadrupedal robots and have not been extended to wall-supported bipedal locomotion, where multi-rate interactions with the environment introduce fundamentally different control challenges.

Bipedal gaits of quadrupedal robots have been investigated primarily through reinforcement learning (RL)-based control approaches, which have demonstrated impressive performance in both simulation and hardware. For example, [25] proposed a multi-modal locomotion framework that combines a hand-engineered transition strategy with a learning-based controller. In [26], a transfer learning framework was introduced to enable quadrupedal robots to jump or walk on their hind legs. An RL-based bi-level framework was developed in [27] to achieve agile, human-like bipedal locomotion. [28] examined the benefits of leveraging symmetry in model-free RL to improve gait robustness, while [29] introduced an adversarial motion priors-based method for bipedal locomotion, validated in simulation. Most recently, [30] presented a learning framework with barrier-based style rewards that enabled quadrupedal, tripod, and bipedal gaits across legged robots. In addition, an RL-based algorithm was developed in [31] to endow the quadrupedal robot KLEIYN with wall-climbing capabilities.

Beyond learning-based methods, MPC-based approaches have also been explored for bipedal locomotion of quadrupedal robots. For instance, [32] proposed a distributed MPC algorithm based on the alternating direction method of multipliers (ADMM) to generate bipedal gaits. More recently, [33] introduced a contact-implicit MPC framework to enable front-leg rearing motions in quadrupedal robots.

Compared to free bipedal locomotion, wall-supported bipedal walking for quadrupedal robots introduces unique challenges, as it involves multi-contact gaits composed of four-contact, three-contact, and double-contact phases. In this setting, the rear legs enable bipedal locomotion, while the front legs exert forces against nearby walls to provide additional stability and propulsion. This renders the system inherently hybrid, with unilateral constraints acting on both the front and rear limbs. To address these challenges, we aim to design a unified MR-NMPC algorithm that simultaneously plans the optimal state trajectory, reaction forces, and foot placements for both front and rear legs.

1.2 Contributions↩︎

This paper presents a layered optimal control framework for wall-supported bipedal locomotion in quadrupedal robots (Fig. 1). The main contributions are:

• A multi-rate nonlinear model predictive controller (MR-NMPC) that jointly optimizes continuous-time CoM motion, torso orientation, and environment reaction forces (ERFs), together with discrete-time contact-point trajectories for foot placement using a single rigid-body dynamics model.

• A multi-rate formulation that incorporates foothold evolution into the prediction dynamics by propagating contact-point locations across successive gait domains using step-length control inputs.

• A nonlinear whole-body controller (WBC), formulated as a convex quadratic program with virtual constraints, that tracks the optimized state, ERF, and contact-point trajectories at 1 kHz while enforcing the full-order robot dynamics.

• Validation in simulation on a Unitree A1 robot demonstrating wall-supported bipedal locomotion over rough terrain and under external disturbances, with comparisons against a heuristic Raibert-based foot-placement controller.

2 High-Level MR-NMPC↩︎

This section establishes the mathematical foundation of the proposed high-level MR-NMPC framework (see Fig. 2). We formulate the framework for general nonlinear discrete-time dynamical systems that originate from template models of locomotion. In particular, we consider a multi-rate discrete-time system \(\Sigma\) with two sets of control inputs, defined as \[\label{eq:state95eq} \Sigma: \quad x(t+1) = f\left(x(t),u(t),v(t)\right),\tag{1}\] where \(t \in \mathbb{Z}_{\geq 0}:= \{0,1,\dots\}\) denotes the discrete-time index, \(x(t) \in \mathcal{X}\subset \mathbb{R}^{n_{x}}\) is the state vector, \(u(t) \in \mathcal{U}\subset \mathbb{R}^{n_{u}}\) is the fast control input, and \(v(t) \in \mathcal{V}\subset \mathbb{R}^{n_{v}}\) is the slow control input. The sets \(\mathcal{X}\), \(\mathcal{U}\), and \(\mathcal{V}\) represent the feasible states, admissible fast inputs, and admissible slow inputs, respectively.

The multi-rate structure arises from the fact that the fast input \(u(t)\) can be updated at every time step, whereas the slow input \(v(t)\) is allowed to change only at a predetermined subset of switching instants, denoted by \[\label{eq:switching95times} \mathcal{T}:= \{0 < t_{0} < t_{1} < t_{2} < \cdots\} \subset \mathbb{Z}_{\geq 0}.\tag{2}\] Between these switching times, the slow input remains constant, yielding a piecewise-constant profile. Formally, \[\label{eq:piecewise95constant} v(t+1) = v(t), \quad \text{if } t+1 \notin \mathcal{T}.\tag{3}\]

Returning to the general formulation in 1 , we are interested in the following problem.

Problem 1 is formulated for the original nonlinear, time-invariant system 1 with both fast and slow inputs. However, enforcing the piecewise-constant behavior of the slow input \(v(t)\) in 3 introduces significant challenges for real-time trajectory planning. To overcome this difficulty, we reformulate the problem using an indicator-function representation. This approach transforms the original steering problem into an equivalent formulation for a time-varying augmented nonlinear system with augmented inputs that are allowed to update at every time step, in contrast to the original formulation.

Specifically, we define the indicator function \(\delta_{\mathcal{T}} : \mathbb{Z}_{\geq 0}\rightarrow \{0,1\}\) as \[\label{eq:indicator95func} \delta_{\mathcal{T}}(t) := \begin{cases} 1, & t \in \mathcal{T}\\ 0, & t \notin \mathcal{T}, \end{cases}\tag{4}\] which specifies whether the current time \(t\) corresponds to a switching instant. The evolution of the slow input can then be written as \[\label{eq:v95dynamics} v(t+1) = v(t) + \delta_{\mathcal{T}}(t+1)\,\Delta v(t),\tag{5}\] where \(\Delta v(t) \in \mathbb{R}^{n_{v}}\) denotes the update applied at switching instants. Thus, \(v(t)\) remains constant between updates and changes only when \(\delta_{\mathcal{T}}(t+1) = 1\). From 5 , the trajectory of the slow input can be explicitly written as \[\label{eq:v95closed95form} v(t) = v(0) + \sum_{k=0}^{t-1} \delta_{\mathcal{T}}(k+1)\,\Delta v(k), \quad \forall t \geq 1,\tag{6}\] where \(v(0)\) denotes the initial value of the slow input.

Incorporating 1 and 5 yields the augmented state model \[\label{eq:aug95dynamics} \Sigma^{a}:\;\begin{cases} x(t+1) = f\!\left(x(t),u(t),v(t)\right), \\[6pt] v(t+1) = v(t) + \delta_{\mathcal{T}}(t+1)\,\Delta v(t), \end{cases}\tag{7}\] where the pair \(x^{a}(t):=\textrm{col}(x(t),v(t)) \in \mathcal{X}\times \mathcal{V}\) defines the augmented state and \(u^{a}(t):=\textrm{col}(u(t),\Delta v(t)) \in \mathcal{U}\times \mathbb{R}^{n_{v}}\) defines the augmented control input. Since the indicator function \(\delta_{\mathcal{T}}(t)\) vanishes outside the switching set \(\mathcal{T}\), both \((u(t),\Delta v(t))\) can be interpreted as fast inputs in this reformulated system. In particular, the evolution of \(v(t)\) is effectively gated by \(\delta_{\mathcal{T}}(t)\), ensuring that updates occur only at switching instants. For the SRB dynamics in Example 1, the augmented state is formed by stacking the CoM positions, CoM velocities, Euler angles, angular velocities, and Cartesian coordinates of the stance limbs’ end-effectors. We now present the following problem, which is equivalent to the original Problem 1.

To address Problem 2, we propose the following real-time and unified NMPC \[\begin{align} {4} &\min_{(x(\cdot),v(\cdot),u(\cdot),\Delta v(\cdot))} &&\mathcal{L}_{\textrm{terminal},x}\left(x_{t+N|t}\right) + \mathcal{L}_{\textrm{terminal},v}\left(v_{t+N|t}\right) \nonumber\\ & &&+ \sum_{k=0}^{N-1} \mathcal{L}_{\textrm{stage},x}\left(x_{t+k|t},u_{t+k|t}\right) \nonumber\\ & &&+ \sum_{k=0}^{N-1} \mathcal{L}_{\textrm{stage},v}\left(v_{t+k|t},\Delta v_{t+k|t}\right) \nonumber\\ &\textrm{s.t.} && \!\!\!\!\!x_{t+k+1|t} = f\left(x_{t+k|t},u_{t+k|t},v_{t+k|t}\right) \nonumber\\ & &&\!\!\!\!\! v_{t+k+1|t} = v_{t+k|t} + \delta_{\mathcal{T}}(t+k+1)\,\Delta v_{t+k|t} \nonumber\\ & &&\!\!\!\!\! E_{u}(t+k)\,u_{t+k|t} =0 \nonumber\\ & &&\!\!\!\!\! E_{\Delta v}(t+k)\,\Delta v_{t+k|t} =0 \nonumber\\ & &&\!\!\!\!\! x_{t+k|t} \in \mathcal{X}, \quad u_{t+k|t} \in \mathcal{U}, \quad v_{t+k|t} \in \mathcal{V},\label{eq:MR95NMPC} \end{align}\tag{8}\] where \(N\) denotes the control horizon. The variables \(x_{t+k|t}\), \(u_{t+k|t}\), \(v_{t+k|t}\), and \(\Delta v_{t+k|t}\) represent the predicted state, fast input, slow input, and slow-input update at stage \(k\), computed at time \(t\) using the augmented prediction model 7 . The initial condition is set by the measured variables at time \(t\), namely \(x_{t|t}=x(t)\) and \(v_{t|t}=v(t)\). In addition, \(x(\cdot)\), \(u(\cdot)\), \(v(\cdot)\), and \(\Delta v(\cdot)\) denote the corresponding trajectories over the control horizon. The cost function consists of the terminal and stage costs for both \(x\) and \(v\) variables, defined as \(\mathcal{L}_{\textrm{terminal},x} := \|x_{t+N|t} - x_{t+N|t}^{\textrm{ref}}\|_{P_{x}}^{2}\), \(\mathcal{L}_{\textrm{terminal},v} := \|v_{t+N|t} - v_{t+N|t}^{\textrm{ref}}\|_{P_{v}}^{2}\), and \(\mathcal{L}_{\textrm{stage},x} := \|x_{t+k|t} - x^{\textrm{ref}}_{t+k|t}\|_{Q_{x}}^{2} + \|u_{t+k|t}\|_{R_{u}}^{2}\), \(\mathcal{L}_{\textrm{stage},v}:= \|v_{t+k|t} - v^{\textrm{ref}}_{t+k|t}\|_{Q_{v}}^{2} + \|\Delta v_{t+k|t}-\Delta v^{\textrm{ref}}\|_{R_{\Delta v}}^{2}\) where \(P_{x}\), \(P_{v}\), \(Q_{x}\), \(Q_{v}\), \(R_{u}\), and \(R_{\Delta v}\) are positive-definite weighting matrices. Notably, the framework utilizes steplength obtained by Raibert heuristic as \(\Delta v^{\textrm{ref}}\) and the corresponding foot evolution as \(v^{\textrm{ref}}_{*}\). Throughout this paper, we use the notation \(\|z\|_{Q}^{2} := z^\top Q z\) for any vector \(z\).

The equality constraints of the NMPC are induced by the nonlinear, time-varying augmented model 7 , together with the assignment conditions \(E_{u}(t+k)\,u_{t+k|t}=0\) and \(E_{\Delta v}(t+k)\,\Delta v_{t+k|t}=0\), which enforce availability of the control inputs depending on the contact configuration. For example, in the SRB dynamics, the sequence of contact points with the environment changes over the control horizon, which restricts the admissible ERFs or foot placements at each stage \(t+k\). Here, \(E_{u}(t+k)\) and \(E_{\Delta v}(t+k)\) are time-varying binary matrices (composed of zeros and ones) that constrain unavailable inputs to zero. The inequality constraints originate from state feasibility and input admissibility.

Example 2 (Locomotion Pattern). For the wall-supported bipedal locomotion of the quadrupedal robot, we consider a multi-contact locomotion pattern denoted by \(\mathcal{G}\). The periodic gait cycle is modeled as a sequence of three phases: a four-contact phase, followed by a three-contact phase, and concluding with a double-contact phase (see Fig. 3). The durations of the four-contact, three-contact, and double-contact phases are set to 50 ms, 30 ms, and 120 ms, respectively. The sampling time used to discretize the SRB dynamics is \(T_{s}=10\) ms, corresponding to solving the MR-NMPC at 100 Hz. The control horizon is chosen as \(N=20\) samples, which corresponds to a prediction window of \(N T_{s}=200\) ms. Consequently, the MR-NMPC optimizes trajectories across multiple phases of the gait cycle. For example, when the NMPC is solved in the middle of the four-contact phase (for example at \(30\) ms time step), the prediction horizon encompasses the current four-contact phase, the subsequent three-contact and double-contact phases of the same cycle, as well as part of four-contact phase (3 time steps) in the next cycle.

The resulting optimal CoM and orientation trajectories, together with the ERF trajectories (fast inputs) and limb end-effector trajectories (slow inputs), are provided to the low-level nonlinear whole-body controller (WBC) for tracking. The numerical implementation details and the solver used for the high-level MR-NMPC will be discussed in Section 4.1.

3 Low-Level Nonlinear WBC↩︎

This section describes the low-level nonlinear whole-body controller (WBC) developed for the full-order dynamics of the quadrupedal robot during wall-supported bipedal locomotion. The controller is adapted from [13], [34], with modifications tailored to the proposed framework.

The full-order floating-base dynamics of the quadrupedal robot during wall-supported locomotion are expressed using the Euler–Lagrange equations and the principle of virtual work as \[\begin{align} {4} & D(q)\,\ddot{q} + H(q,\dot{q}) = B\,\tau + \sum_{\ell\in\mathcal{C}} \left(J^{\ell}(q)\right)^\top f^{\ell} \tag{9}\\ & \ddot{r}^{\textrm{foot},\ell} = 0, \quad \forall \ell \in \mathcal{C}, \tag{10} \end{align}\] where \(q \in \mathcal{Q} \subset \mathbb{R}^{n_{q}}\) denotes the generalized coordinates, \(\tau \in \mathcal{T} \subset \mathbb{R}^{n_{\tau}}\) represents the joint torques, and \(f^{\ell} \in \mathbb{R}^{3}\) is the ERF applied at the stance limb end-effector \(\ell \in \mathcal{C}\). The matrix \(D(q) \in \mathbb{R}^{n_{q}\times n_{q}}\) is the positive-definite mass–inertia matrix, \(H(q,\dot{q}) \in \mathbb{R}^{n_{q}}\) collects the Coriolis, centrifugal, and gravitational terms, \(B \in \mathbb{R}^{n_{q}\times n_{\tau}}\) is the input distribution matrix, and \(J^{\ell}(q) \in \mathbb{R}^{3\times n_{q}}\) denotes the contact Jacobian matrix associated with stance limb \(\ell \in \mathcal{C}\). The contact constraints are enforced by requiring the stance limb end effector accelerations to vanish, as specified in 10 .

To bridge the gap between reduced- and full-order locomotion models, we define a set of outputs, referred to as virtual constraints [35], which impose the full-order model to track the optimal trajectories prescribed by the high-level MR-NMPC. Specifically, we define \[\label{eq:virtual95constraints} y(t,q) := y_{a}(q) - y_{\textrm{des}}(t),\tag{11}\] where \(y_{a}(q)\) denotes the actual output variables to be controlled (controlled variables), and \(y_{\textrm{des}}(t)\) represents their desired evolution along the gait. The controlled variables include the CoM position, torso orientation (Euler angles), and the Cartesian coordinates of the swing limb end effectors. Optimal values \((x_{t+1|t}^{\star},u_{t|t}^{\star},v_{t+1|t}^{\star},\Delta v_{t|t}^{\star})\) from the SRB dynamics are used to construct the desired output function \(y_{\textrm{des}}(t)\). For the swing limb end effectors, we employ a Bézier polynomial that interpolates between the current position of the limb and the upcoming placement specified in \(v_{t+1|t}\).

We next formulate the following real-time convex quadratic program (QP) to solve for the joint-level torques: \[\begin{align} {4} & \min_{(\tau,f,\delta)}&& \frac{\gamma_{1}}{2} \|\tau\|^{2} + \frac{\gamma_{2}}{2} \|f-f_{\textrm{des}}\|^{2} && + \frac{\gamma_{3}}{2} \|\delta\|^{2} \nonumber\\ &\textrm{s.t.} && \ddot{y} + K_{D}\, \dot{y} + K_{P}\, y = \delta && \textrm{(Output dynamics)} \nonumber\\ & && \ddot{r}^{\textrm{foot},\ell} = 0,\quad \forall \ell \in \mathcal{C} && \textrm{(No slippage)} \nonumber\\ & && \tau \in \mathcal{T}, \quad\,\, f \in \mathcal{FC} && \textrm{(Feasibility)}, \label{eq:QP95WBC} \end{align}\tag{12}\] where \(\gamma_{1}\), \(\gamma_{2}\), and \(\gamma_{3}\) are positive weighting factors, \(f := \textrm{col}\{f^{\ell} \mid \ell \in \mathcal{C}\}\) denotes the stacked vector of all ERFs at stance limbs, \(f_{\textrm{des}}\) represents the desired ERFs prescribed by the high-level MR-NMPC; \(\delta\) is the slack variable that ensures feasibility of the output dynamics, and \(\mathcal{FC}\) denotes the linearized friction cone constraints. The first equality constraint enforces the desired output dynamics, \(\ddot{y} + K_{D}\,\dot{y} + K_{P}\,y = \delta\), where \(K_{P}\) and \(K_{D}\) are positive-definite gain matrices for output regulation, and \(\delta\) is used to relax the dynamics when necessary. The second equality constraint imposes zero acceleration of the stance limb end effectors. Finally, the inequality constraints ensure the feasibility of the joint-level torques and stance-limb ERFs. The cost function minimizes a weighted sum of the squared 2-norms of the joint torques, the ERF tracking error, and the slack variables, with relative importance determined by the weighting factors \(\gamma_{1}, \gamma_{2}, \gamma_{3}\). We remark that \(\ddot{y}\) and \(\ddot{r}^{\textrm{foot},\ell}\) are affine functions of \((\tau,f)\) under the Lagrangian dynamics 9 ; therefore, the optimization problem 12 is a convex QP, which will be reliably solved in real time at 1 kHz.

4 Simulation Results↩︎

To evaluate the proposed footstep planner, we conduct simulations of upright, wall-supported walking using a Unitree A1 quadruped within the RaiSim environment [36]. The simulation environment consists of a \(45\) cm narrow corridor between parallel walls with unknown terrain featuring randomly placed wooden blocks similar to the one shown in Fig. 1 to simulate a challenging, obstacle-laden path.

4.1 Setup and Controller Synthesis↩︎

This work employs the Unitree A1 quadruped for numerical validation of the proposed MR-NMPC framework. The \(13\) kg A1 (\(0.25\) m tall) features 18 DoFs, including 12 actuated leg joints (hip pitch, hip roll, and knee pitch). Numerical validation is performed using the RaiSim physics engine. The hyperparameters of the MR-NMPC for the SRB component are chosen as \(Q^{\textrm{SRB}}=\textrm{block diag}\{Q^{\textrm{SRB}}_{p},Q^{\textrm{SRB}}_{\dot{p}},Q^{\textrm{SRB}}_{\theta},Q^{\textrm{SRB}}_{\omega}\}\) with \(Q^{\textrm{SRB}}_{p}=\textrm{diag}\{1 \textrm{e}4,5\textrm{e}4,1\textrm{e}4\}\), \(Q^{\textrm{SRB}}_{\dot{p}}=\textrm{diag}\{1\textrm{e}5,1\textrm{e}4,1\textrm{e}4\}\), \(Q^{\textrm{SRB}}_{\theta}=\textrm{diag}\{8\textrm{e}4,8\textrm{e}5,3\textrm{e}4\}\), and \(Q^{\textrm{SRB}}_{\omega}=\textrm{diag}\{1\textrm{e}2,1\textrm{e}2,1\textrm{e}2\}\). The terminal cost and control penalty are set to \(P^{\textrm{SRB}}=Q^{\textrm{SRB}}\) and \(R^{\textrm{SRB}}= 0.01\,\mathbb{I}_{12\times12}\), respectively.

For the footstep placement and step length optimization of the SRB model, the weighting matrices \(Q_v\) and \(R_{\Delta v}\) from are defined as block-diagonal structures \(Q_v=\textrm{block diag}\{Q^{\textrm{FR}},Q^{\textrm{FL}},Q^{\textrm{RR}},Q^{\textrm{RL}}\}\), \(R_{\Delta v}=\textrm{block diag}\{Q_v^{\textrm{FR}},Q_v^{\textrm{FL}},Q_v^{\textrm{RR}},Q_v^{\textrm{RL}}\}\). To ensure symmetric locomotion behavior, the weights are set uniformly across all four limbs with \(Q^i=1\textrm{e}3\) and \(Q_v^i=1\textrm{e}4\) for \(i \in \{\textrm{FR},\textrm{FL},\textrm{RR},\textrm{RL}\}\). To maintain the focus on high-speed longitudinal agility, we reduce the dimensionality of the slow-input update \(\Delta v\) by considering only the footstep position and step length along the x-direction of motion. Consequently, the decision variables for footstep planning are constrained to the sagittal plane. The high-level MR-NMPC is implemented using the CasADi framework [37] with the IPOPT [38] interior-point solver on a desktop PC (Intel Core i9-12900F, 64 GB RAM). The solver achieves a mean execution time of 8.26 ms with a standard deviation of 0.87 ms, which ensures that the high-level planning loop consistently operates at 100 Hz.

4.2 Comparative Benchmark and Baseline Validation↩︎

Previous MPC formulations optimize future footholds as independent variables or static parameters within the prediction horizon [39], [40]. In contrast, the proposed MR-NMPC models foothold evolution as part of the prediction dynamics by propagating future foothold locations across successive gait domains using optimized step-length inputs.

To evaluate the benefits of this formulation, we benchmark the proposed controller against the classical Raibert heuristic [41], which computes the step adjustment as \(\Delta x=\frac{\dot{x}T_s}{2}+k_v(\dot{x}-\dot{x}_{\mathrm{des}})\), where \(k_v=\sqrt{h/g}\).

Both controllers are first evaluated under nominal conditions (flat terrain with no external disturbances). As shown in Fig. 4, the MR-NMPC and Raibert heuristic achieve nearly identical reference velocity (\(v_{\mathrm{des}}=0.8\) m/s) and base orientation tracking, establishing a common baseline. This enables the performance gains in the subsequent disturbed and uneven-terrain scenarios to be attributed to the proposed predictive foothold formulation rather than controller tuning.

4.3 Adaptive Footstep Placement and Stability↩︎

The robustness of the proposed MR-NMPC framework during obstacle negotiation is evaluated by traversing a series of piled wooden blocks. The wooden blocks are \(2\) cm tall, \(14\) cm wide, stacked up to 3 layers and cover \(70\%\) of the test track. The MR-NMPC utilizes full-state feedback and a 20-step prediction horizon to maintain the robot within the safety envelope and track the reference velocity as shown in Fig. 5 (a). To analyze the mechanism behind this stability, we examine the longitudinal foot offsets relative to the CoM, defined for the front (\(r_x^f\)) and rear (\(r_x^r\)) limbs as: \[r^f_x = p^{\textrm{FR}}_x - x_{\textrm{CoM}}, \quad r^r_x = p^{\textrm{RR}}_x - x_{\textrm{CoM}}.\]

As illustrated in Fig. 5 (b), the MR-NMPC diverges significantly from the Raibert heuristic in its treatment of the rear limbs. While the front foot placement (\(r_x^f\)) remains somewhat comparable between the two methods, the NMPC proactively modulates the rear foot placement (\(r_x^r\)) to accommodate the undulating surface. By decoupling the front and rear step lengths, the optimizer can dynamically expand or contract the support polygon. This flexibility allows the framework to generate corrective restorative torques that are unattainable under the fixed-symmetry constraints of the Raibert heuristic, effectively “steering” the robot back to the reference path despite the treacherous terrain.

4.4 Comparative and Quantitative Analysis↩︎

We evaluate the proposed MR-NMPC against the Raibert baseline across 250 randomized \(7\) m tracks featuring stochastic wooden block placement. One such terrain is shown in Fig. 1. The reference velocity is ramped up in steps of \(0.05\) m/s every two gait domains to a maximum of \(0.8\) m/s to test the algorithm across a wide speed regime. To maintain a rigorous benchmark, we define a successful trial based on a multi-modal stability envelope: nominal CoM height, lateral tracking margins, and kinematic safety limits. A trial is terminated upon the first constraint violation, and the traversal distance is recorded.

Figure 6 illustrates the survival probability across these 250 trials. While both planners demonstrate comparable reliability in the low-velocity regime (up to \(2\) m of travel), the Raibert heuristic exhibits a sharp degradation in success rate as the reference velocity increases. This divergence highlights a fundamental limitation of velocity-based heuristics in agile locomotion: at higher speeds, the increased distance covered per stance phase necessitates a predictive look-ahead to manage the coupled base-and-footstep dynamics. By co-optimizing the full robot state—the MR-NMPC proactively adapts foothold locations based on the anticipated contact configuration to allow the system to navigate irregular terrain at high speeds where reactive feedback alone is insufficient.

4.5 Disturbance Rejection and Constraint Awareness↩︎

To evaluate the robustness of the MR-NMPC against high-magnitude perturbations, we applied a \(50\) N sinusoidal force from \(t\!=\!8\) s to \(9\) s. As shown in Fig. 7 (a), the Raibert heuristic (blue) responds to the velocity spike with long step lengths that trigger joint-limit saturation, resulting in a subsequent instability as illustrated in Fig. 7 (b). In contrast, the MR-NMPC (red) leverages its prediction horizon to identify an optimal sequence of smaller, asymmetric steps. By co-optimizing state feedback with anticipated contact timings, the proposed MR-NMPC mitigates the disturbances while respecting all safety and kinematic constraints.

5 Conclusions↩︎

This paper introduces an MR-NMPC framework that unifies high-level footstep planning and low-level ERF computation within a single optimization loop. By utilizing indicator function with SRB model, the high level multi-rate NMPC dynamically plans both the slow updating discrete-time trajectories of the contact points and the fast updating continuous-time trajectory of the SRB states. By incorporating contact-point planning within the optimal control framework, this architecture provides the flexibility to employ asymmetric step lengths and user-defined contact sequences for negotiating challenging terrain. The effectiveness of the proposed framework was validated through extensive simulations using the Unitree A1 platform. The results demonstrate that the MR-NMPC framework significantly enhances the robot’s ability to navigate rough terrain and reject external disturbances in wall-supported configurations. Most notably, the proposed approach achieved a 2.9 times higher success rate in negotiating irregular terrain at high speeds compared to conventional MPC methods utilizing heuristic foot placement. This underscores the importance of unified trajectory and contact optimization for highly dynamic, non-standard locomotion tasks.

Preliminary hardware experiments indicate the real-time feasibility of the proposed framework, enabling the robot to transition from a crouched posture to upright wall-supported bipedal locomotion and maintain balance for over 60 s. The primary limitation is velocity estimation noise, which induces forward pitching and prevents the stable limit-cycle behavior observed in simulation. Future work will focus on improving state estimation through Kalman filtering and sensor fusion, increasing the front-foot contact area to enhance wall-support stability, and incorporating perception-aware planning to handle uncertainty in terrain geometry and contact conditions. We also plan to perform extensive hardware validation and comparative studies against stronger optimization-based baselines to further evaluate the benefits of the proposed multi-rate formulation.

References↩︎