2.1. Simplified dynamic model

Designing a balance controller for quadruped robots using traditional control methods is usually challenging due to the hybrid nature of the repeated switching between contact and swing modes of the robot leg, which results in irregularity [Reference Bledt, Powell, Katz, Carlo, Wensing and Kim46]. In this regard, full dynamics modeling generally leads to high complexity in control design. We adopt a simplified virtual force control model to optimize foot-end ground reaction forces (GRFs), as shown in Figure 2. Positioning the knee actuator on the hip reduces limb inertia, allowing leg influence simplification. Thus, we reduce the complete dynamics to a simplified model with only the center of mass (COM) and contact feet, thereby controlling the COM motion by planning GRFs. The simplified model is presented as

(1) \begin{equation} \begin{bmatrix} {\mathbf{I}}_3 & \cdots & {\mathbf{I}}_3 \\ [\boldsymbol{p}_1-\boldsymbol{p}_c]\times & \cdots & [\boldsymbol{p}_4-\boldsymbol{p}_c]\times \end{bmatrix}\boldsymbol{F}=\begin{bmatrix} m(\ddot {\boldsymbol{p}}_c+\boldsymbol{g}) \\ \boldsymbol{I}_b\dot {\boldsymbol{\omega }}_b \end{bmatrix} \end{equation}

where $\boldsymbol{F} = [f_1^T, f_2^T, f_3^T, f_4^T]^T \in \mathbb{R}^{12 \times 1}$ explicitly denotes the concatenated vector of the GRFs acting on the four contact feet; $\boldsymbol{p}_{i}$ , $i\in \{1,2,3,4\}$ denote the foot positions, $\boldsymbol{p}_c$ denotes COM position, $[\boldsymbol{p}_i - \boldsymbol{p}_c]^\times$ in Eq. (1) denotes the skew-symmetric matrix corresponding to the vector $[\boldsymbol{p}_i - \boldsymbol{p}_c]$ ; $\boldsymbol{\omega }_b$ denotes body angular velocity, $m$ is the robot’s total mass, $\boldsymbol{I}_b$ denotes the centroidal rotational inertia and $\boldsymbol{g}$ denotes the gravity vector. Note that even though the leg’s influence is disregarded in the COM balance controller problem, we still utilize a complete leg-dynamic model in the single-leg controller.

Figure 2.

Graphical interpretation of the simplified dynamic model.

Figure 2 Long description

The diagram illustrates a quadruped robot with four legs, each having a contact point labeled 1 through 4. The robot's body is represented by a rectangular box with a center of mass (COM) marked inside. The diagram includes vectors and forces such as ground reaction forces (GRF) at each contact point, gravitational force (g), and other forces and vectors like f1, f4, r1, r2, r3, r4, and p1, p2, p3, p4. The coordinate systems {W} and {B} are shown, with axes labeled Xw, Yw, Zw for the world frame and Xb, Yb, Zb for the body frame. The diagram also includes angular velocity (ωb) and linear velocity (vb) vectors.

2.2. Motor-leg model

Accurate motor-leg modeling is essential for active compliance/impedance control of quadruped robots. The motor-leg model treats the leg’s swing and stance motion as a motor-load dynamic system (with torque input and angular velocity output), connected by a high-stiffness linear spring transmission, as shown in Figure 3. Since the direct measurements of joint encoders are on the output end of the speed reducer on the investigated SCIT robot, the reduction ratio does not appear in the final model. Similar to the model in [Reference Boaventura, Buchli, Semini and Caldwell47], the feedback-incorporated rotation model of the equivalent load is developed as

(2) \begin{equation} \dot {\tau _l}=K_w(\omega _m-\omega _l) \end{equation}

In Eq. (2), $\tau _l$ denotes the actuation torque exactly acted on the equivalent load, $K_w$ denotes the coefficient in the transmission mechanism, $\omega _m$ denotes the output angular velocity of the motor, and $\omega _l$ is the actual angular velocity of the foot end. Accordingly, the complete motor-actuation control scheme for single-leg motion is designed as in Figure 2, where $\tau _{cf}$ denotes the total Coulomb friction of the motor shaft and load, $\tau _{out}$ denotes the output torque provided by the controller.

Figure 3.

The equivalent motor-load control loop for single-leg control.

Figure 3 Long description

The diagram illustrates the equivalent motor-load control loop for single-leg control. It features various components and their interactions, including torque inputs, rotational speeds, and control elements. The motor side includes a moment of inertia represented by 1/JmS and a damping coefficient Bm. The load side includes a similar moment of inertia 1/JlS and a damping coefficient Bl. The diagram shows the flow of torque and rotational speed through different stages, with feedback loops and control elements like Kw and 1/s. Arrows indicate the direction of flow and interaction between these components, highlighting the dynamic control process in legged robots.

Then, the input-to-output transfer function of the total leg system can be developed as

(3) \begin{equation} \frac {\omega _l(s)}{\tau _{out}(s)}=\frac {K_w}{(K_w+J_ls^2+B_ls)(J_ms+B_m)+K_w(J_ls+B_l)} \end{equation}

where $\tau _{out}$ denotes the output torque provided by the controller, $J_m, B_m$ are the inherent inertia and viscosity constant of the motor, $J_l, B_l$ are the inertia and viscosity constant of the load.

Note that the actuation torque generated by the leg control system is proportional to the current, but the accurate input–output relation is influenced by an intricate electrical dynamics within the driving circuit. For the simplicity and implementability of the designed controller, this work adopts the current-based controller form, where the PI gains are designed with respect to the motor current feedback and torque input.

2.3. Friction torques modeling for motor-leg control

For motor-load control, the frictional torques include two components: Coulomb and viscous [Reference Tzes, Peng and Guthy48], which affect both the rotation of the motor shaft and the equivalent load. In this section, we consider the full frictional torque on the motor as

(4) \begin{equation} \tau _{fm}=\tau _{cfm}+B_m\omega _m \end{equation}

where $\tau _{fm}$ denotes the full frictional torque acted on the motor shaft, and $\tau _{cfm}$ stands for the Coulomb component. Similarly, the full frictional torque acted on the load is

(5) \begin{equation} \tau _{fl}=\tau _{cfl}+B_l\omega _l \end{equation}

where $\tau _{fl}$ denotes the full frictional torque acted on the load, with $\tau _{cfl}$ standing for the Coulomb component. Combining the above torques, the full frictional torque to the motor-load control system is described by

(6) \begin{equation} \tau _{f}=\tau _{cfm}+B_m\omega _m+\tau _{cfl}+B_l\omega _l \end{equation}

Further, taking both the frictional torque and the control torque, the total torque acted on the load is presented as

(7) \begin{equation} \tau _{l}=\tau _{out}-\tau _{fm}-J_m\dot {\omega }_m \end{equation}

Direct numerical differentiation of optical encoder signals significantly amplifies high-frequency noise and introduces phase delays, which severely degrades the friction compensation torque [Reference Lee and Song49]. Similarly, measuring acceleration for the SCIT Dog presents significant difficulties. To circumvent these issues, we have designed a joint angular acceleration input system based on the Bezier curve [Reference Carlo, Wensing, Katz, Bledt and Kim50]. The Bezier curve provides the desired angular acceleration, angular velocity, and position trajectory in the swing phase. This approach eliminates the need for direct acceleration measurements. Instead, the Bezier curve trajectories rely on intuitive parameters, such as step length and step height, rather than angular joint displacement, to obtain accurate joint acceleration data.

3.1. Inner loop of torque control

The motor system receives voltage as its raw control input, with generated torque proportionally dependent on winding current. Nevertheless, the voltage-current dynamics involve considerable complexity, and direct voltage regulation frequently induces undesirable torque fluctuations. For the sake of high control bandwidth, the locomotion control system incorporates an optimized PI-based current regulation scheme. This controller processes real-time measurements from current sensors, though these readings contain substantial Gaussian-distributed noise. Following the implementation of conventional PID control laws, it is observed that during violent ground impacts, the derivative (D) term strongly amplifies high-frequency mechanical noise, which can lead to torque saturation and destabilize the dynamic performance in leg locomotion. In contrast, the integral (I) term proves essential for minimizing steady-state deviations under sustained loads. Consequently, PI-based torque regulation presents superior practicality relative to P, PD, or PID alternatives. While PI control may introduce minor current oscillations, these perturbations remain inconsequential in practice since the current loop operates as an inner subsystem. Experimental validation confirms that internal loop oscillations negligibly affect outer loop performance. The control architecture further incorporates post-PI filtering and output saturation stages.

Figure 5 shows the designed inner loop with PI control for actuation torque control. Therein, $\tau _{out}$ serves as the motor-leg subsystem’s input torque; $K_I$ represents the motor’s torque constant; $I_{in}$ and $I_{out}$ denote the commanded current and the measured winding current, respectively; $K_p$ and $K_i$ correspond to the proportional and integral gain coefficients, respectively.

Figure 5.

PI current controller.

Figure 5 Long description

The diagram illustrates a proportional-integral (PI) current controller. It begins with a reference input labeled tau ref, which is divided by a constant KI. The result is subtracted from the input current Iin. The error signal is then processed through two parallel paths: one through a proportional gain block labeled Kp and the other through an integral gain block labeled Ki followed by an integrator. The outputs of these paths are summed and passed through a low-pass filter. The filtered signal is then fed into a saturation block to limit its amplitude, producing the output current Iout. Finally, the output current is multiplied by a constant KI to produce the output torque tau out.

3.2. Torque compensation for motor and load friction

To incorporate the previously developed friction torque model and compensation for the effect of frictional torques, both motor and load angular velocities must be measured. For swing leg operation, the compensation law takes the form of

(8) \begin{equation} \tau _\alpha =(J_m+J_l)\theta _{ref}+\tau _{cfm}+B_m\omega _m+\tau _{cfl}+B_l\omega _l+\tau _{cor}+\tau _g \end{equation}

where $\tau _{cor}$ and $\tau _g$ denote the torques generated by Coriolis force and gravity, respectively. The motor’s angular velocity $\omega _m$ is directly measurable from encoder data, with $\omega _l$ assumed to approximate $\omega _m$ . For stance legs, the compensation control law is substituted with

(9) \begin{equation} \tau _\alpha =\tau _{blc}+\tau _{cfm}+B_m\omega _m+\tau _{cfl}+B_l\omega _l+\tau _{cor}+\tau _g \end{equation}

In this case, $\tau _{blc}$ represents the output torque of the balance controller.

3.3. Nonlinear active compliance and impedance control

Foot-ground interaction presents a crucial challenge requiring careful consideration in this study. These dynamic collisions may compromise structural integrity, particularly at vulnerable connection points. Even without immediate failure, sustained operation can induce mechanical play, progressively increasing transmission backlash. Such impacts additionally generate destabilizing oscillations in the robot’s center of mass, potentially leading to gait instability or complete loss of balance. Consequently, effective ground force modulation through adjustable joint compliance/impedance becomes essential for robotic systems.

Inspired by biological quadrupeds that dynamically regulate limb stiffness through muscular adaptation, robotic implementations must emulate this adaptive compliance capability for impact absorption. Post-impact, compliance should be progressively diminished to maintain precise foot positioning. Two distinct compliance approaches exist: passive implementations utilizing mechanical springs between actuators and loads, and active solutions incorporating software-adjustable compliance parameters within torque control algorithms. The active approach enables real-time stiffness tuning, providing robotic limbs with canine-like force interaction capabilities. This active compliance mechanism is achieved through position and velocity feedback loops, mathematically expressed as

(10) \begin{equation} \boldsymbol{\tau }_{ref}=\boldsymbol{J}^T[\boldsymbol{K}_{vs,m}(\boldsymbol{x}_{des}-\boldsymbol{x})+\boldsymbol{K}_{d,m}(\dot {\boldsymbol{x}}_{des}-\dot {\boldsymbol{x}})]+\boldsymbol{\tau }_{\alpha } \end{equation}

The compliance control law utilizes the following components: $\boldsymbol{J}$ , the foot Jacobian matrix; $\boldsymbol{K}_{vs,m}$ $\in$ $R^{3\times 3}$ , position-dependent stiffness matrix; $\boldsymbol{K}_{d,m}$ $\in$ $R^{3\times 3}$ , velocity-dependent damping matrix; $\boldsymbol{x}_{des}$ and $\boldsymbol{x}$ $\in$ $R^3$ , reference and measured foot positions. Since precise steady-state tracking is not required for compliance control, the integral term is intentionally omitted. The proposed control strategy effectively creates a virtual spring-damper system along each leg axis, with the physical parameters given by elastic coefficient $K_{vs}$ and damping coefficient $K_d$ , which take place in the diagonal elements of the stiffness and damping matrices $\boldsymbol{K}_{vs,m}$ and $\boldsymbol{K}_{d,m}$ , respectively. The proportional term provides spring-like restorative forces, while the derivative term generates velocity-dependent damping. Figure 6 demonstrates this virtual mechanical element’s implementation.

Figure 6.

Leg and virtual spring-damper.

Figure 6 Long description

The image shows a close-up view of a quadruped robot leg with a virtual spring and damping mechanism. The leg consists of a metallic structure with a knee joint and a foot end. The virtual spring and damping mechanism is illustrated with a zigzag line and arrows indicating the direction of force. The diagram highlights the interaction between the robot's leg and the ground, emphasizing the use of virtual forces to control the ground reaction forces (GRFs). The simplified model focuses on the center of mass (COM) and contact feet, allowing for the control of the COM motion by planning GRFs. The knee actuator is positioned on the hip to reduce limb inertia, simplifying the leg's influence. The overall structure demonstrates a method to optimize foot-end GRFs for better balance and control of the quadruped robot.

In the validation for the virtual spring characteristics of the legs, the single leg of the investigated quadruped robot is manually compressed, and the output force is measured while the desired foot position is maintained using joint angle and leg kinematics. Repeated experiments are conducted to obtain the equivalent spring property of the legs, as shown in Figure 7.

Figure 7.

Virtual spring characteristics. Given the virtual linear spring $K_{vs}=75N/m$ and the nonlinear spring $K_{vs}=125(e^{|x_{des}^r-x^r|}-1) N/m$ . ( $x^r$ means the radial distance of the position) (a) Linear spring characteristic. (b) Nonlinear spring characteristic.

Figure 7 Long description

A scatter plot compares actual force and expected force for linear and nonlinear springs. The x-axis represents spring displacement in meters, ranging from negative to positive values, while the y-axis represents force in newtons, ranging from zero to thirty. The plot includes two subplots: (a) for the linear spring and (b) for the nonlinear spring. Each subplot shows a series of data points marked with asterisks for actual force and a red line for expected force. In both subplots, the data points closely follow the expected force line, indicating a strong correlation. The linear spring subplot shows a symmetric parabolic pattern, while the nonlinear spring subplot displays a more complex, asymmetric pattern. The data points are densely packed, suggesting a high degree of accuracy in the measurements. All values are approximated.

Figure 7 compares the force–displacement relationships of linear versus exponential nonlinear virtual springs, with radial displacement plotted along the horizontal axis and corresponding force output on the vertical axis. The experimental configuration employs deliberately reduced spring coefficients to accentuate tracking errors, demonstrating that both spring types achieve near-ideal coefficient tracking performance.

Achieving concurrent impact absorption and precise trajectory following remains essential for both spring variants. During locomotion, stance phase dynamics require effective impact absorption, while swing phase demands accurate foot positioning. Although linear springs provide basic impact mitigation, prior research [Reference Dong, Sang, Lee and Kim32, Reference Semini, Barasuol, Boaventura, Frigerio, Focchi, Caldwell and Buchli51] demonstrates that nonlinear springs offer superior performance. Our approach implements nonlinear compliance control on quadrupedal systems, modeling each leg as a spring-damper system. The post-impact dynamics resemble second-order system responses, where knee load peaks before oscillating and stabilizing at gravitational equilibrium. Effective impact mitigation requires minimizing both load peaks and center-of-mass perturbations.

Meanwhile, optimizing the virtual spring stiffness $K_{vs}$ is crucial for achieving effective impact mitigation while maintaining accurate position tracking. Lowering $K_{vs}$ can reduce the impact peak and load fluctuations, but this can also increase the position tracking error and cause greater changes in the COM position. On the other hand, higher stiffness can reduce the position tracking error but lead to greater impact peak and load fluctuations. Therefore, it is necessary to actively adjust $K_{vs}$ after the impact to reduce the impact peak and maintain the stability of the COM. In the swing phase, high stiffness should be maintained to ensure accurate position tracking. To ensure a smooth stiffness transition without inducing abrupt mechanical jerks or discontinuities in the control torque, an infinitely differentiable and bounded exponential function is adopted. For the above objectives, the stiffness adaptive law is formulated as

(11) \begin{equation} K_{vs}=K_{vs1}+\frac {2K_{vs2}}{\textrm {exp}\big(K_{cv}\big(x_{des}^r-x^r\big)\big)+\textrm {exp}\big(\!-K_{cv}\big(x_{des}^r-x^r\big)\big)} \end{equation}

where $K_{cv}$ is a human-induced fixed parameter. The stiffness coefficient $K_{vs}$ dynamically adjusts between baseline value $K_{vs1}$ and maximum value $(K_{vs1}+K_{vs2})$ . Upon impact, the radial position error $(\dot {x}_{des}^r-\dot {x^r})$ surges, triggering automatic reduction of $K_{vs}$ . As the center-of-mass recovers post-impact, $K_{vs}$ gradually increases. During the swing phase with minimal position error, the system naturally maintains elevated stiffness.

Leg motion transitions between compression and extension states post-impact. The compression phase absorbs impact energy by leg retraction, while the extension phase restores the COM position. During compression, substantial radial velocity error $(\dot {x}_{des}^r-\dot {x^r})$ generates large corrective torque $\tau _{ref}$ , necessitating reduced $K_d$ to mitigate peak forces. Conversely, the extension phase requires increased $K_d$ for COM recovery. This adaptive damping adjustment correlates with leg state, implementable through a deviation-driven relation

(12) \begin{equation} K_d=K_{d1}+K_{d2}\frac {1+\textrm {erf}\big(\!-K_{cd}\big(\dot {x}_{des}^r-\dot {x^r}\big)\big)}{2} \end{equation}

where $\textrm {erf}(x)=\frac {2}{\sqrt {\pi }}\int _{0}^{x}e^{-{\eta }^2}\textrm {d}{\eta }$ , $K_{cd}$ is a fixed scaling factor. The damping coefficient $K_d$ operates within bounds $K_{d1}$ and $(K_{d1}+K_{d2})$ , typically settling at $K_d\approx K_{d1}+K_{d2}/2$ during the leg’s swing stage.

While joint loading varies with respect to the configuration of the leg robot, the knee joint predominantly absorbs locomotion impacts. Our control strategy dynamically modulates $K_d$ according to unidirectional velocity errors $(\dot {x}_{des}^r-\dot {x^r})$ , making it suitable for swing hip joints but not span hip joints. For span hip joints experiencing bidirectional impacts, $K_d$ remains fixed. However, the stiffness parameter $K_{vs}$ adapts regardless of impact direction and applies universally to all joints.

3.4. Stability analysis

In the following subsection, the closed-loop stability of the above-developed single-leg control system is derived, considering the effect of adaptive impedance. Although the motor-leg plant represents a third-order linear system, the controller introduces nonlinear, time-varying behavior through parameters $K_{vs}$ and $K_d$ operating within defined bounds. Consequently, rigorous stability analysis becomes essential for properly selecting $K_{vs1}$ , $K_{vs2}$ , $K_{d1}$ , and $K_{d2}$ . The Popov criterion is exploited to verify the stability of this nonlinear controller-linear plant combination.

Introducing the assumption that $K_d$ is time-invariant and $K_{vs}$ is time-varying gain describing the nonlinearity, one can obtain the stability margin in the value of $K_{vs}$ , related to $K_d$ . The analysis begins by reformulating the closed-loop system, where $K_d$ and the derivative terms are incorporated into the total leg system (3), yielding the following transfer function

(13) \begin{align} P(s)&=\frac {\dfrac {K_w}{\big(K_w+J_ls^2+B_ls\big)(J_ms+B_m)+K_w(J_ls+B_l)}}{s\left [1+K_d\dfrac {K_w}{\big(K_w+J_ls^2+B_ls\big)(J_ms+B_m)+K_w(J_ls+B_l)}\right ]} \nonumber \\ &=\dfrac {K_w}{as^4+bs^3+cs^2+ds} \end{align}

where

(14) \begin{equation} \left \{ \begin{array}{lr} a=J_lJ_m \\[3pt] b=J_mB_l+J_lB_m \\[3pt] c=J_mK_w+J_lK_w+B_lB_m \\[3pt] d=K_wB_m+K_wB_l+K_dK_w \end{array} \right . \end{equation}

Stability is analyzed using the Popov plot, graphing $\omega$ . $Im\left [P(j\omega )\right ]$ and $Re\left [P(j\omega )\right ]$ denote the real and imaginary parts of $P(j\omega )$ , respectively, where $Im\left [P(j\omega )\right ]$ versus $Re\left [P(j\omega )\right ]$ , $\omega$ denotes the frequency. The Popov criterion provides a sufficient condition for global asymptotic stability: if the entire Popov plot lies strictly to the right of a line (the Popov line) passing through the point $(\!-K_{vs,max}^{-1},j0)$ on the real axis, stability is assured. Here, $K_{vs,max}$ represents the maximum permissible $K_{vs}$ , the required quantities are computed as:

(15) \begin{equation} \left \{ \begin{array}{lr} Re\left [P(j\omega )\right ]=\dfrac {K_w(a\omega ^4-c\omega ^2)}{(a\omega ^4-c\omega ^2)^2+\omega ^2(d-b\omega ^2)^2}\\[11pt] {\omega }Im\left [P(j\omega )\right ]=\dfrac {-K_w{\omega ^2}(d-b\omega ^2)}{(a\omega ^4-c\omega ^2)^2+\omega ^2(d-b\omega ^2)^2} \end{array} \right . \end{equation}

Let ${\omega }Im\left [P(j\omega _0)\right ]$ = 0, the Popov stability criterion is described by

(16) \begin{equation} \left \{ \begin{array}{lr} \omega _0=\sqrt {\dfrac {d}{b}}\\[8pt] Re\left [P(j\omega _0)\right ]=\dfrac {K_w}{\dfrac {ad^2}{b^2}-\dfrac {cd}{b}}\\[21pt] K_{vs,max}=-\dfrac {1}{ReP(j\omega _0)} \end{array} \right . \end{equation}

From Eq. (15), $K_{vs,max}$ is calculated by

(17) \begin{align} K_{vs,max}&=\frac {(J_mK_w+J_lK_w+B_lB_m)(B_m+B_l+K_d)}{J_mB_l+J_lB_m}\nonumber \\[3pt] &\quad - \frac {J_mJ_lK_w(B_m+B_l+K_d)^2}{(J_mB_l+J_lB_m)^2} \end{align}

where $K_d\gt 0,K_{vs,max}\gt K_{vs}\gt 0$ .

A similar methodology can establish the relationship between limits on $K_d$ and a fixed $K_{vs}$ (treating $K_d$ as the nonlinear gain), although the resulting expressions become considerably more complex. The functional relationship defined by Eq. (17) between $K_{vs,max}$ and $K_d$ provides the theoretical foundation for selecting feasible controller parameters. In practical implementation, factors such as the inner-loop PI controller dynamics and unmodeled system characteristics necessitate a conservative approach. Parameter values should be chosen significantly below the theoretical stability limit derived here to ensure robust performance [Reference Lee, Isaacs, Higgins, Biewener and Mcgowan52].

Figure 8.

Graphical interpretation of Popov criterion.

Figure 8 Long description

The diagram consists of two main parts labeled as stable cases and unstable cases. Each part contains two subplots. The stable cases show a Nyquist plot with a green line indicating robustness margin and another plot showing a marginally stable system. The unstable cases display two Nyquist plots where the system is unstable. The plots feature axes labeled with Im for the imaginary axis and Re for the real axis. Key points such as -1/Kvs and ω=0 are marked, and the function G*(jω) is plotted in each subplot. Arrows indicate the direction of the frequency sweep. The stable cases demonstrate how the plot avoids encircling the critical point -1/Kvs, while the unstable cases show encirclement, indicating instability.

4.1. PD control and QP problem

The controller employs PD regulation for both COM position and body orientation while respecting foot friction cone constraints. The control law is formulated as:

(18) \begin{equation} \begin{bmatrix} \ddot {\boldsymbol{p}}_{c,des} \\[3pt] \dot {\boldsymbol{\omega }}_{b,des} \end{bmatrix} = \begin{bmatrix} \boldsymbol{K}_{p,p}(\boldsymbol{p}_{c,des}-\boldsymbol{p}_c)+\boldsymbol{K}_{d,p}(\dot {\boldsymbol{p}}_{c,des}-\dot {\boldsymbol{p}}_c) \\[3pt] \boldsymbol{K}_{p,\omega }\textrm { log}(\boldsymbol{R}_{des}\boldsymbol{R}^T)+\boldsymbol{K}_{d,\omega }(\boldsymbol{\omega }_{b,des}-\boldsymbol{\omega }_b) \end{bmatrix} \end{equation}

where $\boldsymbol{R}_{des} \in SO(3)$ and $\boldsymbol{R} \in SO(3)$ represent desired and actual body orientations (expressed as Lie group matrices), and orientation error is computed via the logarithm map [Reference Barrau and Bonnabel54].

As the dynamic model is linear, the balance controller can be formulated as a quadratic programming (QP) problem. The QP problem is solved to obtain the optimized contact force that controls the COM [Reference Gehring, Coros, Hutter, Bloesch, Hoepflinger and Siegwart55]. The QP problem is described by

(19) \begin{equation} \left \{ \begin{array}{lr} \boldsymbol{A}=\begin{bmatrix} {\mathbf{I}}_3 & \cdots & {\mathbf{I}}_3 \\[8pt] [\boldsymbol{p}_1-\boldsymbol{p}_c] ^\times & \cdots & [\boldsymbol{p}_4-\boldsymbol{p}_c] ^\times \end{bmatrix}\\[13pt] \boldsymbol{b}_{des}=\begin{bmatrix} m(\ddot {\boldsymbol{p}}_{c,des}+\boldsymbol{g}) \\ \boldsymbol{I}_b\dot {\boldsymbol{\omega }}_{b,des} \end{bmatrix}\\[9pt] \boldsymbol{F}^{*}=\textrm {min} \quad (\boldsymbol{A}\boldsymbol{F}-\boldsymbol{b}_{des})^T\boldsymbol{S}(\boldsymbol{A}\boldsymbol{F}-\boldsymbol{b}_{des})\\[5pt] +\alpha \Vert \boldsymbol{F}\Vert ^2+\beta \Vert \boldsymbol{F}-\boldsymbol{F}^{*}_{pre}\Vert ^2 \quad \textrm {s.t.} \quad \boldsymbol{C}\boldsymbol{F}\leq \boldsymbol{d} \end{array} \right . \end{equation}

where $F^*$ represents the optimal contact force vector to be solved. The matrix $S \in \mathbb{R}^{6 \times 6}$ is a positive definite weight matrix balancing the importance of each control quantity, and $\alpha$ and $\beta$ are positive weighting coefficients. To satisfy the friction cone constraint, the linear inequality constraint $CF \le d$ is enforced in the optimization problem [Reference Focchi, del Prete, Havoutis, Featherstone, Caldwell and Semini53], where $C$ represents the friction cone approximation matrix, and $d$ is the corresponding boundary vector.

While the Popov criterion ensures the high-frequency local asymptotic stability of the single-leg impedance control against sudden impacts (as discussed in Section 3.4), a legged robot inherently experiences inter-leg coupling effects transmitted through the rigid torso. In our hierarchical framework, these global coupling effects are natively managed by the centralized balance controller. Specifically, the QP problem formulated in (19) utilizes the centralized rigid-body dynamics (1), which explicitly accounts for the kinematic and dynamic coupling of all stance legs. Any unmodeled disturbance or impact absorbed by one leg translates into a state deviation of the COM. The QP solver then simultaneously optimizes the 12-DoF ground reaction force vector $F$ to generate the exact compensating wrench required by the PD control law (18), thereby dynamically redistributing the coupled loads across all contact feet to maintain whole-body stability.

4.2. Nonlinear PD parameters adjustment

To account for potential external disturbances to the robot body, we have adjusted the PD parameters in the balance controller. As the COM position has minimal influence on robot balance, the $\boldsymbol{K}_p$ value should be reduced when the deviation between the COM and the expected position is significant to prevent the robot from generating an excessive response. Conversely, the velocity of the COM is crucial for robot balance. When the deviation between the velocity of the COM and the expected velocity is significant, the value of $\boldsymbol{K}_d$ should be increased to enhance the robot’s response and prevent tipping. For example, in the x-direction, the adjustment strategy can be expressed in a form similar to Eq. (12)

(20) \begin{equation} K_{px}=K_{px1}+\frac {2K_{px2}}{\textrm {exp}(K_{cpx}(p_{x,des}-p_x))+\textrm {exp}(\!-K_{cpx}(p_{x,des}-p_x))} \end{equation}

(21) \begin{equation} K_{dx}=K_{dx1}-\frac {2K_{dx2}}{\textrm {exp}(K_{cdx}(\dot {p}_{x,des}-\dot {p}_x))+\textrm { exp}(\!-K_{cdx}(\dot {p}_{x,des}-\dot {p}_x))} \end{equation}

where $K_{px1},K_{px2},K_{cpx},K_{dx1},K_{dx2},K_{cdx}$ are constants. $K_{px}$ varies between $K_{px1}$ and $(K_{px1}+K_{px2})$ . $K_{dx}$ varies between $(K_{dx1}-K_{dx2})$ and $K_{dx1}$ . Similar strategies are adopted for the control of $y$ -direction and $z$ -direction.

4.3. Footstep locations

The location of each footstep is determined by a linear combination of the Raibert heuristic [Reference Raibert and Tello10] and the capture point [Reference Pratt, Carff, Drakunov and Goswami56] formulations, which takes into account the position of the corresponding hip and the velocity of the COM. This is the classical approach described in ref. [Reference Bledt, Powell, Katz, Carlo, Wensing and Kim46]. Specifically, the position of the $i$ -th footstep on the 2D ground plane is computed as

(22) \begin{equation} \boldsymbol{p}_{step,i}=\boldsymbol{p}_{hip,i}+\frac {T_c}{2}\dot {\boldsymbol{p}}_{c,des}+\sqrt {\frac {z_0}{\Vert \boldsymbol{g}\Vert }}(\dot {\boldsymbol{p}}_c-\dot {\boldsymbol{p}}_{c,des}) \end{equation}

where $\boldsymbol{p}_{step,i}$ is the position of footstep $i$ , $\boldsymbol{p}_{hip,i}$ is the position of the corresponding hip $i$ , $T_c$ is the contact (stance) phase time, $z_0$ is the height of locomotion.

The classical approach to footstep location selection neglects the impact of body posture, which is a crucial factor in enhancing the robot’s resilience to perturbations. Therefore, it is essential to consider the influence of body posture on footstep location determination. Specifically, when the robot’s roll angle is significant, the footstep location should be shifted to the ”left” or ”right” direction to improve the robot’s potential to balance in the subsequent gait cycle. Based on these, we propose the following heuristic adjustment scheme:

(23) \begin{align} \boldsymbol{p}_{step,i}&=\boldsymbol{p}_{hip,i}+\frac {T_c}{2}\dot {\boldsymbol{p}}_{c,des}+\sqrt {\frac {z_0}{\Vert \boldsymbol{g}\Vert }}(\dot {\boldsymbol{p}}_c-\dot {\boldsymbol{p}}_{c,des}) \nonumber \\[4pt] & \quad +\begin{bmatrix} K_{p,\theta }(\theta _{des}-\theta )+K_{d,\theta }(\dot {\theta }_{des}-\dot {\theta }) \\[4pt] K_{p,\phi }(\phi _{des}-\phi )+K_{d,\phi }(\dot {\phi }_{des}-\dot {\phi }) \end{bmatrix} \end{align}

where $\theta$ and $\phi$ represent the pitch and roll angle of the robot body, respectively.

It is noteworthy that a quadruped robot that is controlled by 12 active actuators possesses 18 degrees of freedom. In this study, we propose Eqs. (11) and (12) as a parameter adjustment scheme to enhance the robot’s ground interaction capability by addressing the reaction force from the ground. Additionally, we describe Eqs. (20) and (21) as a parameter adjustment scheme to enhance the robot’s impact mitigation capability when faced with external perturbations. These equations enable the robot to attain a center of mass with stronger impact mitigation capability. Furthermore, Eq. (23) outlines the footstep position selection strategy while considering the body posture of the robot. This approach provides the robot with an impact mitigation capability of 18 degrees of freedom.

Figure 9.

The overall autonomous navigation framework.

Figure 9 Long description

The diagram illustrates the overall autonomous navigation framework for a robot. It includes components such as Desired Center of Mass, Gait Scheduler, Balance Control, Swing Leg Control, State Estimation, Contact Detection, and Active Compliance/Impedance Leg Control. The framework shows the flow of information and control signals between these components, with desired inputs leading to state estimation and contact detection, which then inform balance and swing leg control. The output from these controls is sent to the Active Compliance/Impedance Leg Control, which interacts with the SCIT Dog Robot.

5.1. PI torque controller performance

The PI controller operates at 1 kHz (5 $\times$ faster than outer-loop controllers), with dedicated microcontrollers handling motor-level execution. Figure 10 demonstrates tracking performance under sinusoidal torque reference signals. While minor output oscillations are observed, the controller maintains satisfactory reference tracking. Crucially, inner-loop oscillations show negligible propagation to outer control layers.

Figure 10.

PI torque controller performance.

Figure 10 Long description

The line graph displays two data series: actual torque in blue and expected torque in red, plotted over a time range from 0 to 10 on the x-axis. The y-axis represents torque values ranging from -0.2 to 0.4. The blue line, representing actual torque, shows slight fluctuations around the red line, which represents expected torque. A zoomed-in inset focuses on the time interval from 6.5 to 8, highlighting detailed fluctuations in the torque values. The inset shows that actual torque closely follows expected torque but with minor deviations. All values are approximated.

Despite the presence of output torque oscillations, the PI torque controller demonstrates favorable tracking performance, with minimal impact from the oscillations of the PI controller in the inner loop on the outer loop controller.

5.2. Position tracking experiment

Evaluating swing-phase performance (Section 3.3), tests suspended legs under gravity with radial sinusoidal position commands. Figure 11 compares linear ( $K_{{vs}} = 750 \, \mathrm{N/m}$ ) and nonlinear ( $K_{{vs1}} = 500\,\mathrm{N/m}$ and $K_{{vs2}} = 500 \, \mathrm{N/m}$ ), The nonlinear scheme demonstrates superior tracking, reducing error from $\pm$ 4.4 to $\pm$ 3.6 mm (Figure 11c). During swing, high stiffness is autonomously sustained.

Figure 11.

Position tracking experiment (a) Linear spring position tracking. (b) Nonlinear spring position tracking. (c) Position tracking deviation.

Figure 11 Long description

The image contains three line graphs labeled (a), (b), and (c). Graph (a) shows the position tracking of a linear spring, with the x-axis representing time in seconds and the y-axis representing position in millimeters. The blue line indicates the expected position, while the red line shows the actual position. An inset zooms in on a specific time range for detailed comparison. Graph (b) depicts the position tracking of a nonlinear spring, with similar axes and color coding as graph (a). Graph (c) illustrates the position tracking deviation over time, comparing a nonlinear spring in blue and a linear spring in red. The x-axis represents time in seconds, and the y-axis represents position deviation in millimeters. The graphs highlight the differences in performance between linear and nonlinear springs in tracking expected positions and deviations over time.

5.3. Contact impact mitigation experiments

The robot assumed a static standing posture on all four legs. A controlled vertical impulse (peak magnitudes of 500 and 800 N, duration 200 ms) was applied to the robot body via the feed-forward controller. This impulse induced an upward jump followed by ground impact. The impact event initiation was defined as the moment all four feet made contact with the ground upon landing. The COM position was estimated using leg kinematics. Results are presented in Figures 12 and 13.

Figure 12.

Impact mitigation experiments (Feed-forward vertical impulse 500 $N$ peak and 200 ms duration) (a) Average knee load. (b) Vertical direction COM deviation. (c) Time series of $K_{vs}$ . (d) Time series of $K_d$ .

Figure 12 Long description

The image contains four line graphs. The first graph (a) shows the average knee load in Newton-meters over time in seconds for different stiffness and damping coefficient settings. The second graph (b) illustrates the center of mass deviation in centimeters over time in seconds for the same settings. The third graph (c) presents the time series of the stiffness coefficient in Newton-meters over time in seconds. The fourth graph (d) displays the time series of the damping coefficient in Newton-seconds per meter over time in seconds. Each graph uses different colored lines to represent various combinations of stiffness and damping coefficients, including variable settings. The graphs collectively demonstrate the impact mitigation performance under different conditions.

This experiment evaluated the effectiveness of linear versus nonlinear virtual spring implementations. Four linear spring-damper configurations were tested: (i) High stiffness, low damping ( $K_{vs}=1000 \, \mathrm{N/m}$ , $K_d=50 \, \mathrm{Ns/m}$ , red curves); (ii) High stiffness, high damping ( $K_{vs}=1000 \, \mathrm{N/m}$ , $K_d=150 \, \mathrm{Ns/m}$ , green curves); (iii) Low stiffness, low damping ( $K_{vs}=500 \, \mathrm{N/m}$ , $K_d=50 \, \mathrm{Ns/m}$ , blue curves); (iv) Low stiffness, high damping ( $K_{vs}=500 \, \mathrm{N/m}$ , $K_d=150 \, \mathrm{Ns/m}$ , yellow curves). Besides, the nonlinear spring with $K_{vs}=variable (500\sim 1000 \, \mathrm{N/m}$ , $K_d=variable (50\sim 150 \, \mathrm{Ns/m}$ is represented by black curves. According to Figures 12(b) and 13(b), the Higher joint stiffness reduced the steady-state error in the vertical COM position caused by gravity. However, higher $K_{vs}$ also resulted in significantly larger peak knee loads during impact (Figures 12(a) and 13(a)). The nonlinear controller demonstrated superior impact mitigation performance. It achieved both smaller peak knee loads and reduced COM displacement compared to all linear configurations (Figures 12(a, b) and 13(a, b)). The statistical results of 20 repeated experiments are listed in Table III. The active compliance/impedance controller dynamically adjusted parameters $K_{vs}$ and $K_d$ in response to the impact: (i) Both $K_{vs}$ and $K_d$ were rapidly reduced to minimize the initial impact peak force; (ii) Following the peak impact force, $K_{vs}$ and $K_d$ were actively increased to stabilize the COM position and dampen oscillations. Besides, it can be obtained that higher $K_{vs}/K_d$ combinations produced excessive average knee loads. On the contrary, lower $K_{vs}/K_d$ combinations led to larger COM deviations in the vertical direction. In conclusion, the experimental results validate the theoretical analysis in Section 3. They conclusively demonstrate that the proposed adaptive compliance/impedance controller offers superior capability for mitigating the detrimental effects of ground impact forces compared to static linear approaches. The nonlinear strategy effectively balances the conflicting requirements of impact absorption and COM stabilization.

Figure 13.

Impact mitigation experiments (Feed-forward vertical impulse 800 $N$ peak and 200 ms duration) (a) Average knee load. (b) Vertical direction COM deviation. (c) Time series of $K_{vs}$ . (d) Time series of $K_d$ .

Figure 13 Long description

The image contains four line graphs. The first graph shows the average knee load over time, with different lines representing various stiffness and damping values. The second graph illustrates the vertical direction center of mass deviation over time, again with multiple lines for different stiffness and damping values. The third graph displays the time series of K vs, showing how this parameter changes over time. The fourth graph presents the time series of K d, indicating the variation of this parameter over time. Each graph has its own unique trends and values, providing insights into the impact mitigation experiments.

Table III.

Statistical results of impact mitigation.

Table III Long description

The table presents statistical results of impact mitigation experiments, comparing peak knee load and peak center of mass (COM) deviation for different controller configurations. It includes data for two scenarios: 500 newtons for 200 milliseconds and 800 newtons for 200 milliseconds. The table has four columns: Controller, Peak knee load Mean, Peak knee load Std, Peak COM deviation Mean, and Peak COM deviation Std. It lists various linear spring-damper configurations with different stiffness (K_vs) and damping (K_d) values, as well as a nonlinear spring configuration. Each row provides the mean and standard deviation of peak knee load and peak COM deviation for the respective controller settings. Notable trends include higher joint stiffness reducing steady-state error in vertical COM position but resulting in larger peak knee loads. The nonlinear controller demonstrates superior performance with smaller peak knee loads and reduced COM displacement compared to linear configurations.

Table IV.

Comparison results with other methods.

Table IV Long description

The table presents a comparison of peak knee load and peak center of mass (COM) deviation across different methods and scenarios. It consists of four rows and six columns. The columns are labeled Scenario, Method, Peak knee Load Mean, Peak knee Load Std, Peak COM deviation Mean, and Peak COM deviation Std. The scenarios include 500 newtons for 200 milliseconds and 800 newtons for 200 milliseconds. The methods compared are Proposed, VIC, and RAIC. For the 500 newton scenario, the Proposed method shows a mean peak knee load of 9.75 with a standard deviation of 0.381, and a mean peak COM deviation of 1.57 with a standard deviation of 0.131. The VIC method has a mean peak knee load of 9.65 with a standard deviation of 0.391, and a mean peak COM deviation of 1.82 with a standard deviation of 0.152. The RAIC method shows a mean peak knee load of 10.77 with a standard deviation of 0.339, and a mean peak COM deviation of 1.50 with a standard deviation of 0.124. For the 800 newton scenario, the Proposed method shows a mean peak knee load of 15.44 with a standard deviation of 0.362, and a mean peak COM deviation of 3.10 with a standard deviation of 0.149. The VIC method has a mean peak knee load of 15.31 with a standard deviation of 0.410, and a mean peak COM deviation of 3.42 with a standard deviation of 0.179. The RAIC method shows a mean peak knee load of 17.47 with a standard deviation of 0.349, and a mean peak COM deviation of 2.72 with a standard deviation of 0.165.

Then, the proposed controller is compared to two recent methods, VIC [Reference Spyrakos-Papastavridis, Childs and Dai57] and RAIC [Reference Song, Fang and Qian58], which have been validated on legged robots and achieved state-of-the-art performance. In this part, the experiment benchmark is set the same as the above impact mitigation experiment. The repeated results are listed in Table IV. As shown in Figure 14, compared with the two previous methods, the proposed method achieves a good balance between the torque load and COM accuracy, which benefits from the adaptive controller parameters in (11) and (12). Specifically, the nonlinear spring controller exhibits the same level of impact on VIC, while the tracking error is reduced by more than $10 \%$ . In addition, although RAIC exhibits better robustness, it is subject to higher conservatism in impact mitigation, so the proposed method achieves a more ideal performance.

Furthermore, to validate the robustness of the proposed method in the movement of quadruped robots, we have conducted repeated experiments under a series of operating conditions, including three kinds of gaits (trotting/galloping/bounding) and two kinds of payloads (with/without a 3 kg payload on body), with the results listed in Table V. As shown in Figure 15, the proposed controller exhibits highly consistent transient behavior and good performance across different loads and gait conditions. Among them, the bound gait experiences the most severe impacts, while the gallop gait is relatively smoother. The addition of extra unknown load increases both the control error and impact to some extent, but does not excessively degrade the overall performance.

In summary, the proposed balance controller enabled the robot to regain balance by adjusting the footstep and optimizing the desired contact force. This was accomplished through a nonlinear PD parameter adjustment strategy, which facilitated the rapid recovery of balance under external disturbance.

5.4. Impact mitigation of external disturbance

In this experiment, the walking robot was subjected to artificial external disturbance, which allowed the effectiveness of the balance controller to be tested. The controller was designed to enable the robot to recover its balance after external disturbances. The recovery capability was demonstrated through a sequence of consecutive frames captured at $t=1s$ and presented in Figure 16 (Using the method of Eqs. (18–21). We conducted multiple repeated experiments and tested the performance on the SCIT Dog using only the balance controller described in Eqs. (18) and (19), and without the balance controller. The results indicate that without the use of a balance controller, robots have almost no ability to withstand external interference. In the case of using only the balance controller described in Eqs. (18) and (19), the robot often collapses when facing external disturbance. In repeated experiments, our proposed balance control method ensures that the robot hardly collapses, indicating that this method effectively improves the robot’s ability to resist external disturbance.

Table V.

Repeated results for different operating conditions.

Table V Long description

The table presents data on peak knee load and peak center of mass (COM) deviation for three gaits (trotting, galloping, bounding) under two payload conditions (0 kilograms and 3 kilograms). It consists of 6 rows and 5 columns. The columns are labeled as Gait, Peak knee load Mean, Peak knee load Std, Peak COM deviation Mean, and Peak COM deviation Std. For 0 kilograms payload, the peak knee load mean values are 11.72 for trotting, 8.593 for galloping, and 14.86 for bounding. The peak COM deviation mean values are 1.35 for trotting, 0.91 for galloping, and 0.71 for bounding. For 3 kilograms payload, the peak knee load mean values are 12.52 for trotting, 14.46 for galloping, and 17.69 for bounding. The peak COM deviation mean values are 2.12 for trotting, 1.12 for galloping, and 0.96 for bounding. The table highlights that the bounding gait experiences the most severe impacts, while the gallop gait is relatively smoother. The addition of a 3-kilogram payload increases both the control error and impact but does not excessively degrade overall performance.

5.5. General platform test and comparison to passive compliance

To verify the applicability of the proposed method on other general quadruped robot platforms besides the SCIT robot, this subsection employs the Cheetah-software simulator [Reference Di Carlo59] for the deployment and validation of the compliance controller on the MIT Cheetah, as shown in Figure 17. Meanwhile, the control performance of the proposed nonlinear virtual components is compared with the passive spring system (PSS) and the virtual spring controller (LSC) [Reference Hyun, Seok, Lee and Kim60] of the MIT Cheetah. The parameter settings of the simulator refer to ref. [Reference Di Carlo59]. As shown in Figure 18, the proposed controller exhibited comparable performance on the MIT robot to that on our self-developed platform, validating the applicability of the proposed method on quadruped robot platforms. Meanwhile, compared with the existing virtual spring method on this platform, the proposed approach can more effectively mitigate impacts, further reducing the load torque on the knee joints. Furthermore, in contrast to passive compliance, the proposed method demonstrates better transient response and achieves higher accuracy in periodic locomotion. Moreover, to evaluate the performance difference in power consumption, the index of average COT (cost of transport, cf. [Reference Hyun, Seok, Lee and Kim60]) in each steady gait period is used, with the result: 0.521 (Proposed), 0.556 (LSC), 0.449 (PSS). By flexibly adjusting the stiffness online, the proposed method achieves higher energy efficiency. On the other hand, in comparison with passive spring systems, this method incurs a higher control cost in order to obtain improvements in accuracy and performance.

Figure 14.

Impact mitigation compared with other active compliance methods (a) Average knee load. (b) Vertical direction COM deviation.

Figure 14 Long description

The image contains two line graphs comparing the performance of different impact mitigation methods. The first graph (a) shows the average knee load over time, measured in Newton-meters, for six different methods: Proposed (S1), RAIC (S1), VIC (S1), Proposed (S2), RAIC (S2), and VIC (S2). The second graph (b) displays the center of mass (COM) deviation in the vertical direction over time, measured in centimeters, for the same six methods. Each method is represented by a distinct color line. The x-axis for both graphs represents time in seconds, ranging from 0 to 0.5 seconds. The y-axis for graph (a) ranges from 0 to 20 Newton-meters, while the y-axis for graph (b) ranges from -4 to 0 centimeters. The graphs illustrate how different methods affect knee load and COM deviation during impact mitigation, highlighting variations in performance and stability.

Figure 15.

Impact mitigation for different conditions (a) Average knee load. (b) Vertical direction COM deviation.

Figure 15 Long description

The image contains two line graphs. The first graph shows the average knee load in Newton-meters over time in seconds for different gaits: trot, gallop, and bound, both with and without payload. The second graph illustrates the center of mass deviation in centimeters over time in seconds for the same gaits and conditions. Each line represents a different gait or gait with payload, with distinct colors for each. The graphs demonstrate how the knee load and center of mass deviation vary over time for each condition.

Figure 16.

External disturbance and balance recovery (a) t=0s, (b) t=0.2s, (c) t=0.4s, (d) t=0.6s, (e) t=0.8s, (f) t=1s.

Figure 16 Long description

The image consists of six photos labeled (a) to (f), each capturing a different moment in time as a robot recovers from an external disturbance. In photo (a), at t=0 seconds, the robot is in its initial position. In photo (b), at t=0.2 seconds, the robot begins to lose balance. In photo (c), at t=0.4 seconds, the robot is midway through its recovery process. In photo (d), at t=0.6 seconds, the robot continues to adjust its posture. In photo (e), at t=0.8 seconds, the robot is almost back to its balanced position. Finally, in photo (f), at t=1 second, the robot has fully recovered and is back in a stable stance.

Figure 17.

Cheetah-software simulator.

Figure 17 Long description

The image depicts a simulation environment where two robotic figures are interacting with a checkered floor. The robots appear to be in a crawling or crouching position, with their limbs extended. There are three colored spheres near the robots: a red sphere, a blue sphere, and a green sphere. The simulation includes a graphical user interface in the top left corner displaying real-time data such as simulation time, simulation rate, and other metrics. The checkered floor consists of alternating black and white squares, providing a clear contrast for the robotic figures and spheres.

Figure 18.

Impact mitigation comparison on MIT Cheetah (a) Average knee load. (b) Vertical direction COM deviation.

Figure 18 Long description

The image contains two line graphs comparing the performance of different methods on the MIT Cheetah robot. The first graph (a) shows the average knee load in Newton-meters over time in seconds, with three lines representing the Proposed method, LSC, and PSS. The second graph (b) displays the center of mass (COM) deviation in centimeters over time in seconds, also with three lines for the Proposed method, LSC, and PSS. The Proposed method generally shows lower knee load and COM deviation compared to LSC and PSS, indicating better impact mitigation and stability. The graphs highlight the effectiveness of the Proposed method in reducing knee load and maintaining COM stability during robot-ground interaction.

5.6. Comparison with MPC and learning-based approaches

To evaluate the proposed method against state-of-the-art locomotion frameworks, we supplemented a comparative study in the MIT Cheetah simulation environment. The selected benchmarks include: (1) A convex MPC baseline [Reference Di Carlo, Wensing, Katz, Bledt and Kim38], which optimizes ground reaction forces over a predefined horizon to rigorously track a reference COM trajectory. (2) A RL baseline [Reference Shao, Jin, Liu, He, Wang and Yang4] inspired by the phase-guided control framework, which utilizes a deep RL policy combined with domain randomization to achieve robust dynamic locomotion. All three control strategies were subjected to the identical 500 N transient vertical impact condition described in Section 5.3.

As illustrated in Figure 19(a), the standard MPC baseline experiences a severe peak knee torque exceeding 24 Nm at the instant of impact. This occurs because its rigid tracking objective aggressively combats the sudden external force to maintain the planned trajectory, leading to severe mechanical rebound. The pre-trained RL policy exhibits a moderate impact peak but is accompanied by high-frequency control noise due to the lower evaluation frequency of the policy network. Conversely, our proposed nonlinear compliance method demonstrates the most favorable impact mitigation, yielding the lowest and smoothest torque peak (around 14 Nm). While Figure 19(b) indicates that our method permits a slightly larger transient COM deviation to effectively dissipate kinetic energy, it successfully prevents destructive mechanical overloading on the transmission system. Furthermore, as analyzed previously, the proposed joint-level reactive controller operates at 1 kHz with a computational time of less than 0.1 ms. This extreme computational efficiency stands in stark contrast to the millisecond-level execution times typically required by MPC solvers or deep neural network evaluations, highlighting the practical superiority of our method for high-frequency, reactive impact mitigation on resource-constrained on-board computers.

Figure 19.

Compared with the MPC and RL approaches.

Figure 19 Long description

The image contains two line graphs. The first graph shows the average knee load in Newton-meters over time in seconds for three different approaches: Proposed, MPC Baseline, and RL Baseline. The Proposed approach is represented by a blue solid line, the MPC Baseline by a red dashed line, and the RL Baseline by a green dashed line. The second graph shows the center of mass deviation in centimeters over time in seconds for the same three approaches. The Proposed approach is represented by a blue solid line, the MPC Baseline by a red dashed line, and the RL Baseline by a green dashed line. The graphs illustrate how the average knee load and vertical direction center of mass deviation change over time for each approach.

5.7. Computational load analysis

In this subsection, the computational load of the proposed controller is evaluated to validate its real-time feasibility for quadruped robots. For the nonlinear compliance controller, the computations of the nonlinear adaptive laws are concentrated on the update of time-varying gains (11, 12) and involve only elementary mathematical operations, incurring negligible computational burden. By test, the control loop runs at a frequency of 1 kHz, and the maximum computation time is less than 0.1 ms. As for the QP problem of balance control (19), an efficient on-board solver is implemented based on OSQP. The size of the QP problem is 12 decision variables and 12 constraints. The average solving time on the Nvidia Xavier on-board computer is 0.24 ms, with the maximum solving time less than 0.49 ms, and its stable operation within the 200 Hz control loop for the balance controller. This data conclusively demonstrates the real-time feasibility of the designed system.