2.1. Whole-body planner (based on MPC)

We formulate the optimal control problem as follows:

(1) \begin{align}\left\{\begin{array}{c@{\quad}c} {\min\limits_{u}} & \Phi (x(T))+\int _{0}^{T}L(x(t),u(t),t)dt\\[6pt]s.t. & {\dot{x}}(t)=f(x(t),u(t),t)\\[6pt] & g(x(t),u(t),t)=0\\[6pt] & h(x(t),u(t),t)\geq 0\\[6pt] & x(0)=x_{0}, \end{array}\right.\end{align}

where $x(t)$ and $u(t)$ represent the state vector and input vector at time step $t$ , while $L(x(t),u(t),t)$ and ${\Phi} (x(T))$ denote the stage cost and terminal cost functions. Our objective is to find a control sequence that minimizes the total cost under the constraints of initial state $x_{0}$ , system dynamics ${\dot{x}}(t)=f(x(t),u(t),t)$ , equality constraints $g(\!\cdot\! )$ , and inequality constraints $h(\!\cdot\! )$ .

2.1.1. System dynamics

This research addresses the dynamic control problem of quadrupedal mobile manipulation platforms using a whole-body centroidal dynamics modeling approach. Following the conclusions in ref. [Reference Di Carlo, Wensing, Katz, Bledt and Kim8], we emphasize that for quadrupedal systems equipped with multi-degree-of-freedom manipulators, the influence of leg mass distribution on system dynamics must be properly accounted for, rather than adopting the simplified massless limb assumption. The derived centroidal dynamics equation is expressed as:

(2a) \begin{equation}\dot{p}_{G}=\sum _{i=1}^{n_{c}}\lambda _{i}+mg\end{equation}

(2b) \begin{equation}\dot{l}_{G}=\sum _{i=1}^{n_{c}}r_{Gi}\times \lambda _{i}\end{equation}

The robot state $x=(h_{G},q_{b},q_{j})\in \mathrm{\mathbb{R}}^{12+{n_{a}}}$ comprises the centroidal momentum $h_{G}$ , base pose $q_{b}$ , and joint positions $q_{j}$ . The centroidal momentum $h_{G}=(p_{G},l_{G})\in \mathrm{\mathbb{R}}^{6}$ consists of linear momentum $p_{G}$ and angular momentum $l_{G}$ . The input vector $u=(\lambda _{i},\cdots ,\lambda _{{n_{c}}},v_{j})\in \mathrm{\mathbb{R}}^{3{n_{c}}+{n_{a}}}$ includes contact forces $\lambda _{i}$ at $n_{c}$ contact points and joint velocities $v_{j}$ , where $r_{Gi}$ denotes the position vector from the i-th contact point to the centroid.

2.1.2. Cost function

We encode all robot-centric and end-effector-centric tasks into a unified cost function formulated as:

(3) \begin{equation} L\left(x\left(t\right) \! ,u\left(t\right) \! ,t\right)=\left\| x_{r}-x_{r}^{ref}\right\| _{Q_{r}}^{2}+\left\| u-u^{ref}\right\| _{R}^{2}+\left\| r_{IE}-r_{IE}^{ref}\right\| _{Q_{eep}}^{2}+\left\| \xi _{IE}\right\| _{Q_{eeo}}^{2} \end{equation}

where $R$ are positive definite matrices; $Q_{r} , Q_{eep} , Q_{eeo}$ are positive semi-definite weighting matrices; $r_{IE},\xi _{IE}\in \mathrm{\mathbb{R}}^{3}$ represent the end-effector’s position and orientation errors respectively.

2.1.3. Constraint functions

These constraints capture the distinct modes of each leg at any given time. Assuming the mode sequence is a predefined temporal function, we denote by $N_{\textit{contact}}$ the set of all closed contacts during time period. The resulting constraints are expressed as:

(4) \begin{equation} \left\{\begin{array}{ll} v_{i}\left(x,u\right)=0,f_{c}{\left(\mu _{c}\right)}\cdot \lambda _{i}\geq 0, & \begin{array}{ll} \textit{if} & i\in N_{\textit{contact}}, \end{array}\\[4pt] v_{i}{\left(x,u\right)}\cdot \hat{n}=s{\left(t\right)},\lambda _{i}=0, & \begin{array}{ll} \textit{if} & i\notin N_{\textit{contact}}, \end{array} \end{array}\right. \end{equation}

where $v_{i}$ denotes the foot-end velocity in the world coordinate frame, enforcing that: (1) stance legs maintain ground contact through $v_{i}=0$ ; (2) swing legs follow the predefined trajectory $s(t)$ along the local surface normal direction $\hat{n}\in \mathrm{\mathbb{R}}^{3}$ . Additionally, zero contact force is enforced for swing legs, while stance legs are subject to friction cone constraints $f_{c}(\mu _{c})$ defined by the surface normal $\hat{n}$ and friction coefficient $\mu _{c}=0.7$ .

2.2. Whole-body controller (based on WBC)

The whole-body controller employs a hierarchical optimization architecture to track reference trajectories for both the base and limbs. This controller solves for optimal generalized accelerations and contact forces through QP formulation with prioritized tasks [Reference Bellicoso, Gehring, Hwangbo, Fankhauser and Hutter23, Reference Xin, Wolfslag, Lin, Tiseo and Mistry24], subsequently computing joint torques via inverse dynamics. The specific task prioritization scheme is illustrated in Figure 1.

2.2.1. Equations of motion

The mechanical system’s equations of motion establish the relationship between generalized acceleration $\dot{u}_{WBC}$ , contact forces $\lambda _{c}$ , and joint torques $\tau$ . By defining $x_{WBC}=[\begin{array}{cc} \dot{u}_{WBC} & \lambda _{c} \end{array}]^{T}$ , the robot’s centroidal dynamics model takes the following form:

(5) \begin{equation}\left[\begin{array}{cc} M_{b}(q) & -J_{b}^{T} \end{array}\right]\cdot x_{WBC}=-h_{b}(q,\dot{q})\end{equation}

where $(\!\cdot\! )_{b}$ denotes base-related components in parameter matrices.

2.2.2. No contact motion

To ensure contact consistency for stance legs (which are constrained from moving), the following tasks must be enforced:

(6) \begin{equation}\left[\begin{array}{cc} J_{i} & 0_{3\times {n_{a}}} \end{array}\right]\cdot x_{WBC}=-\dot{J}_{i}\cdot u_{WBC},~ i\in N_{\textit{contact}}\end{equation}

2.2.3. Torque limits

To enforce joint torque limitations within hardware-specified bounds, we formulate the following joint torque constraint task:

(7) \begin{equation}\left[\!\!\!\begin{array}{cccccc} \begin{array}{cc} & M_{j}(q)\\[3pt] & -M_{j}(q) \end{array} & \begin{array}{cc} & -J_{j}^{T}\\[3pt] & J_{j}^{T} \end{array} \end{array}\right]\cdot x_{WBC}\leq \left(\begin{array}{c} \tau _{\max}\cdot 1_{{n_{a}}\times 1}-h_{j}(q,\dot{q})\\[3pt] -\tau _{\min}\cdot 1_{{n_{a}}\times 1}+h_{j}(q,\dot{q}) \end{array}\right)\end{equation}

where $1_{{n_{a}}\times 1}=[\begin{array}{c} \begin{array}{ccc} 1 & \cdots & 1 \end{array} \end{array}]^{T}\in \mathbb{R}^{{n_{a}}}$ , $(\!\cdot\! )_{j}$ denotes joint-related components in parameter matrices.

2.2.4. Friction pyramid constraint

To prevent slippage, the contact force $\lambda _{c}\in \mathbb{R}^{3\times {n_{c}}}$ generated by the i-th foot must be constrained within the friction cone. To obtain a linear formulation, we approximate the friction cone using a pyramid. Thus, given the friction coefficient $\mu _{c}$ , the friction constraint for the i-th leg’s contact force can be expressed as:

(8) \begin{equation} f_{c}{\left(\mu _{c}\right)}\cdot \lambda _{i}=\left[\begin{array}{l@{\quad}l@{\quad}l} \begin{array}{l} 1\\ 1\\ 0\\ 0\\ 0 \end{array} & \begin{array}{l} 0\\ 0\\ 1\\ 1\\ 0 \end{array} & \begin{array}{l} \mu _{c}\\ -\mu _{c}\\ \mu _{c}\\ -\mu _{c}\\ 1 \end{array} \end{array}\right]\cdot \lambda _{i}\geq 0,\begin{array}{ll} & i\in N_{\textit{contact}} \end{array} \end{equation}

2.2.5. Trajectory tracking

After ensuring dynamic feasibility, hardware torque limits, and non-slipping contact constraints, we perform trajectory tracking for the center of mass (CoM), swing legs, and end-effector [Reference Barasuol, Buchli, Semini, Frigerio, De Pieri and Caldwell25, Reference Hoffman, Laurenzi, Muratore, Tsagarakis and Caldwell26]. Given reference CoM trajectories ${C}^{}{r}{}_{b}^{ref}$ , $\dot{{C}^{}{r}{}}_{b}^{ref}$ and $\ddot{{C}^{}{r}{}}_{b}^{ref}$ in the operational space, the tracking is implemented as follows [Reference Di Carlo, Wensing, Katz, Bledt and Kim8]:

(9a) \begin{equation}\left[\begin{array}{cc} _{C}^{}{J}{}_{b} & 0_{6\times {n_{a}}} \end{array}\right]\cdot x_{WBC}={_{C}^{}{\ddot r}{}}_{b}^{ref}+K_{p,b}\cdot {_{C}^{}{\widetilde r}{}}_{b}+K_{d,b}\cdot {_{C}^{}{\dot{\widetilde r}{}}}_{b}-{_{C}^{}{\dot J}{}}_{b}\cdot u_{WBC}\end{equation}

Similarly, trajectory tracking for both the swing legs and the robotic arm’s end-effector can also be achieved using the above formulation; therefore, we define:

(9b) \begin{equation}\left[\begin{array}{cc} _{C}^{}{J}{}_{i} & 0_{3\times {n_{a}}} \end{array}\right]\cdot x_{WBC}={_{C}^{}{\ddot r}{}}_{i}^{ref}+K_{p,i}\cdot {_{C}^{}{\widetilde r}{}}_{i}+K_{d,i}\cdot {_{C}^{}{\dot{\widetilde r}{}}}_{i}-{_{C}^{}{\dot J}{}}_{i}\cdot u_{WBC},\begin{array}{cc} & i\notin N_{\textit{contact}} \end{array}\end{equation}

(9c) \begin{equation}\left[\begin{array}{cc} _{C}^{}{J}{}_{ee} & 0_{6\times {n_{a}}} \end{array}\right]\cdot x_{WBC}={_{C}^{}{\ddot r}{}}_{ee}^{ref}+K_{p,ee}\cdot {_{C}^{}{\widetilde r}{}}_{ee}+K_{d,ee}\cdot {_{C}^{}{\dot{\widetilde r}{}}}_{ee}-{_{C}^{}{\dot J}{}}_{ee}\cdot u_{WBC}\end{equation}

2.2.6. Contact force tracking

We execute contact force tracking with the lowest priority, which guides the contact force solution while ensuring the primary task objectives proceed unimpeded, thereby accelerating the overall solver convergence.

(10) \begin{equation} \left[\begin{array}{ll} 0_{{n_{a}}\times \left(12+{n_{a}}\right)} & \mathrm{I }_{{n_{a}}\times {n_{a}}} \end{array}\right]\cdot x_{WBC}=\lambda ^{ref} \end{equation}

where $\lambda ^{ref}$ represents the ideal contact force computed by the MPC solver.

3.1. External force estimation model

For floating-base robotic systems, we first establish the standard dynamic equation:

(11) \begin{equation}M(q)\ddot{q}+h(q,\dot{q})=\tau +J^{T}(q)\lambda\end{equation}

By defining the nominal state dynamics (external-force-independent terms) on the left-hand side as: $\tau _{d}=M(q)\ddot{q}+h(q,\dot{q})$ , the dynamic equation simplifies to an external-force-explicit form:

(12) \begin{equation} \lambda =J^{-T}{\left(q\right)}{\left(\tau _{d}-\tau \right) }\end{equation}

In practical systems, accounting for noise interference in joint torque measurements and state estimations, we introduce joint torque residual $\delta \tau =\tau -\hat{\tau }$ and state residual $\delta \tau _{d}=\tau _{d}-\hat{\tau }_{d}$ , where $\hat{\tau }$ and $\hat{\tau }_{d}$ represent the measured values of joint torques and dynamic terms, respectively, yielding the disturbance equation:

(13) \begin{equation} \hat{\tau }_{d}+\delta \tau _{d}=\hat{\tau }+\delta \tau +J^{T}{\left(q\right)}\lambda \end{equation}

Assuming the system noise follows Gaussian distributions: $\delta \tau _{d}\sim \mathrm{N }(0,R_{d}) , \delta \tau \sim \mathrm{N }(0,R)$ , where $\{R_{d},R\}$ are their respective covariance matrices. By merging the measurement noise terms, we ultimately construct an observation model incorporating noise characteristics:

(14) \begin{equation} e=J^{T}\left(q\right)\lambda -\tau ^{\circ } \end{equation}

where $\tau ^{\circ }=\hat{\tau }_{d}-\hat{\tau }$ represents the observed value of the nominal state, and $e=\delta \tau _{d}\oplus \delta \tau$ denotes the system noise. The “ $\oplus$ ” symbol indicates probabilistic superposition between $\delta \tau _{d}$ and $\delta \tau$ . This model probabilistically combines the nominal state observation with system noise, establishing a theoretical foundation for subsequent parameter estimation.

3.2. External force estimation model

The external force estimation problem is transformed into a noise energy minimization process for dynamic systems. Based on the weighted least squares criterion, an optimization objective function incorporating a forgetting factor is constructed:

(15) \begin{equation} \lambda {\left(t_{N}\right)}=\arg \min _{\lambda }\left(\sum _{k=1}^{N}\alpha ^{N-k}\left\| e {\left(t_{k}\right)}\right\| _{W}^{2}\right) \end{equation}

where $\alpha \in (0,1]$ is the forgetting factor, which exponentially weights and decays the influence of historical data; $N$ represents the observation window length, reflecting the time constant of the system’s dynamic characteristics. This weighted least squares problem is solved in real-time using the RLS method, with its iterative form given by:

(16) \begin{equation} \left\{\begin{array}{l} K_{ef}\left(t_{k}\right)=\dfrac{P_{ef}\left(t_{k-1}\right)\cdot \left(J^{T}\left(t_{k}\right)\right)^{T}}{\alpha \cdot E_{\left(6+{n_{a}}\right)\cdot \left(6+{n_{a}}\right)}+J^{T}\left(t_{k}\right)\cdot P_{ef}\left(t_{k-1}\right)\cdot \left(J^{T}\left(t_{k}\right)\right)^{T}}\\[18pt] \lambda {\left(t_{k}\right)}=\lambda {\left(t_{k-1}\right)}+K_{ef}\left(t_{k}\right)\cdot \left(\tau ^{\circ }\left(t_{k}\right)-J^{T}\left(t_{k}\right)\cdot \lambda \left(t_{k-1}\right)\right)\\[15pt] P_{ef}{\left(t_{k}\right)}=\dfrac{P_{ef}{\left(t_{k-1}\right)}-K_{ef}\left(t_{k}\right)\cdot J^{T}{\left(t_{k}\right)}\cdot P_{ef}\left(t_{k-1}\right)}{\alpha } \end{array}\right. \end{equation}

Here, $P_{ef}(t_{k})$ is the error covariance matrix, $K_{ef}(t_{k})$ is the Kalman gain matrix, and $E_{(6+{n_{a}})\cdot (6+{n_{a}})}$ is the identity matrix. Finally, the estimated generalized interaction force $\tilde{\lambda }_{ee}=[\begin{array}{cc} \tilde{f}_{ee} & \tilde{\tau }_{ee} \end{array}]^{T}\in \mathbb{R}^{6}$ at the robot’s end-effector can be directly extracted from the last six elements of $\lambda (t_{k})$ , corresponding to the contact force and torque components.

3.3. The α-adaptive method based on control architecture

The RLS-FF algorithm serves as an efficient online parameter estimation method, particularly suitable for time-varying robotic dynamic parameters. By recursively updating parameters, this method introduces a forgetting factor $\alpha$ to exponentially decay the weight of historical data, thereby assigning higher confidence to new data and enabling effective tracking of time-varying parameters. The value of $\alpha$ critically impacts estimation performance: as $\alpha$ approaches 0, the algorithm exhibits strong tracking capability but increased noise sensitivity, leading to amplified estimation variance and parameter oscillation; as $\alpha$ approaches 1, it demonstrates excellent noise resistance but noticeable parameter update lag.

Conventional force estimation methods typically employ a fixed forgetting factor, which suffers from two major limitations: (1) difficulty in adapting to rapid parameter variations in dynamic interaction scenarios, and (2) failure to fully leverage the hierarchical nature of the control architecture. In contrast, the deep coupling framework proposed in this study adopts a dual-time-scale estimation strategy, enabling high-frequency parameter fine-tuning at the WBC layer (1 kHz) and low-frequency parameter reconstruction at the MPC layer (50–100 Hz), thereby significantly improving system adaptability. Within the MPC-WBC hierarchical control architecture, the MPC module optimizes and generates reference joint trajectories over a finite time horizon, while the WBC module performs local optimization at a higher frequency to achieve accurate trajectory tracking based on the commands from MPC. This division of labor implies that each time MPC updates its planned trajectory, the reference for the WBC module changes, thereby reducing the relevance of previously accumulated historical data for force estimation at the current moment. Accordingly, our method proactively reduces the forgetting factor in the RLS algorithm after each MPC update. This deliberate reduction diminishes the weight of older data, allowing the estimator to adapt more rapidly to dynamic variations introduced by the new trajectory commands. Compared to strategies with a fixed forgetting factor, this approach aligns more naturally with the decision-making cycle of the controller, enhancing the response speed to load variations while maintaining estimation stability.

As shown in the control framework (Figure 2), within each control cycle (1 kHz): (1) the WBC first solves for joint torque commands and sends them to the low-level controller; (2) the force estimator then updates dynamic parameters using the RLS-FF algorithm (computational complexity O(n ²)) based on discrepancies between actual states and the reference model; (3) the refined parameters are applied to the next cycle’s WBC computation.

Figure 2.

Block diagram of the proposed adaptive force-based controller.

Under the synchronous framework, the forgetting factor is set as $\alpha =\alpha _{WBC}$ . Furthermore, since the WBC algorithm operates at a fixed control cycle, Equation (15) is reformulated as $\lambda (N)=\arg \min _{\lambda }(\sum _{k=1}^{N}\alpha _{WBC}^{N-k}\| e(k)\| _{W}^{2})$ , where N denotes the number of cycles. Assuming $J(N)=\sum _{k=1}^{N}\alpha _{WBC}^{N-k}\| e(k)\| _{W}^{2}$ , we derive:

(17) \begin{equation} J {\left(N\right)}=\left[\begin{array}{llll} \alpha _{WBC}^{N-1}{\cdot }{}e^{2}{\left(1\right)} & \cdots & \alpha _{WBC}^{1}\cdot e^{2}{\left(N-1\right)} & e^{2}{\left(N\right) }\end{array}\right]\cdot 1_{N\times 1} \end{equation}

To address the lower update frequency of MPC (100 Hz), a dynamic forgetting factor adjustment strategy is designed: During each MPC policy update, the forgetting factor is automatically reduced to diminish historical data weighting, thereby enhancing the algorithm’s rapid response capability to new policies. This design ensures stable estimation within MPC cycles while maintaining adaptability to parameter changes during policy transitions. The modified forgetting factor is introduced, and the weighted sum of squared prediction errors for each MPC cycle becomes:

(18) \begin{equation} J_{MPC}=\alpha _{MPC}^{1}\cdot \left[\begin{array}{lll} \alpha _{WBC}^{N^{*}-2}{\cdot }{}e^{2}\left(1\right) & \cdots & \alpha _{WBC}^{1}\cdot e^{2}\left(N^{*}-1\right) \end{array}\right]\cdot 1_{\left({N^{*}}-1\right)\times 1}+e^{2}\left(N^{*}\right) \end{equation}

Combining Eq. (17) and (20), we obtain:

(19) \begin{equation} J{\left(N\right)}=\left[\begin{array}{c@{\quad}c@{\quad}c} \alpha_{MPC}^{n^{*}-1}{\cdot }{}J_{MPC}\left(1\right) & \cdots & \alpha _{MPC}^{1}\cdot J_{MPC}{\left(n^{*}\right)} \\[6pt] \alpha _{WBC}^{N-n^{*}N^{*}-1}\cdot e^{2}{\left(n^{*}N^{*}+1\right)} & \cdots & e^{2}{\left(N\right)} \end{array}\right]\cdot 1_{N\times 1} \end{equation}

where $n^{*}$ represents the number of MPC executions. Through this design, whenever the robot’s MPC policy updates, the system dynamically decreases the value of the forgetting factor $\alpha$ , thereby achieving exponential decay weighting of historical data while correspondingly increasing the weight of new data. This dynamic adjustment mechanism enables the force estimation algorithm to rapidly track dynamic parameter changes induced by new MPC strategies, significantly improving response capability to external disturbances while ensuring system stability.

4.1. Whole-body control extension

While traditional WBC can handle basic dynamic constraints, their performance often degrades when facing unknown external disturbances. This study enhances WBC in three key aspects:

First, by treating the external force $\lambda _{ee}=(f_{ee},\tau _{ee})\in \mathrm{\mathbb{R}}^{6}$ acting on the end-effector as a disturbance term in the centroidal momentum equation, the centroidal dynamics undergo the following modification:

(20) \begin{equation}\left[\begin{array}{cc} M_{b}(q) & -J_{b}^{T} \end{array}\right]\cdot x_{WBC}=-h_{b}(q,\dot{q})+J_{ee}^{T}\cdot \lambda _{ee}\end{equation}

where $J_{ee}$ is the Jacobian matrix of the robotic arm’s end-effector. By treating the external force $\lambda _{ee}$ acting on the end-effector as a disturbance term and incorporating it into the Task: Equations of Motion, this enhancement enables the system to accurately reflect actual force conditions. Secondly, regarding constraint handling, we have reformulated the Joint Torque Saturation Task:

(21) \begin{equation}\left[ \!\!\! \begin{array}{cccccc} \begin{array}{cc} &M_{j}(q)\\[6pt] &-M_{j}(q) \end{array} & \begin{array}{cc} & -J_{j}^{T}\\[6pt] &J_{j}^{T} \end{array} \end{array}\right]\cdot x_{WBC}\leq \left(\begin{array}{c} \tau _{\max}\cdot 1_{{n_{a}}\times 1}-h_{j}\left(q,\dot{q}\right)+J_{ee}^{\tau }\cdot F_{ee}\\[6pt] -\tau _{\min}\cdot 1_{{n_{a}}\times 1}+h_{j}{\left(q,\dot{q}\right)}-J_{ee}^{\tau }\cdot F_{ee} \end{array}\right)\end{equation}

This formulation ensures joint torque commands satisfy hardware limitations even under external force disturbances. Finally, a compensation term is introduced in the torque computation step, yielding the ideal joint torque $\tau ^{des}$ as:

(22) \begin{equation}\tau ^{des}=\left[\begin{array}{cc} M(q) & -J_{j}^{T} \end{array}\right]\cdot {x_{WBC}}^{des}+h_{j}(q,\dot{q})-J_{ee}^{T}\cdot F_{ee}\end{equation}

This enhancement enables the WBC to perform real-time compensation for measured external disturbances at 1 kHz control frequency.

4.2. Model predictive control extension

To address the instability of conventional WBC in highly dynamic scenarios, we integrate the extended dynamic model into predictive optimization, combining MPC’s predictive capability with WBC’s real-time control. The MPC generates reference trajectories accounting for future states, while the WBC performs refined optimization based on real-time dynamic constraints. This framework significantly enhances the dynamic performance of quadruped mobile manipulation platforms when handling unknown end-effector loads, enabling rapid stabilization under disturbances.

During robot-environment interaction, the end-effector experiences both linear contact forces and contact torques $\tau _{ee}$ . To accurately model this complete interactive dynamics, we systematically extend the momentum balance equations within the original single rigid-body dynamics framework:

(23) \begin{equation}\begin{cases} & \dot{p}_{G}=\sum\limits _{i=1}^{n_{c}}f_{i}+\lambda _{ee}+mg, i\in N_{\textit{contact}}\\[12pt] & \dot{l}_{G}=\sum\limits _{i=1}^{n_{c}}r_{Gi}\times f_{i}+r_{Gee}\times \lambda _{ee}+\tau _{ee} \end{cases} \end{equation}

where $\lambda _{ee}=[\begin{array}{ll} f_{ee} & \tau _{ee} \end{array}]^{T}\in \mathrm{\mathbb{R}}^{6}$ denotes the end-effector force and $r_{Gee}$ denotes the position of the end-effector relative to the CoM. The input vector is extended to include $u=(\lambda _{{c_{1}}},\cdots ,\lambda _{{c_{{n_{c}}}}},\lambda _{ee},v_{j})\in \mathrm{\mathbb{R}}^{3{n_{c}}+6+{n_{a}}}$ .

In formulating the optimization problem, equality constraints must be applied to the generalized contact forces (including linear force $f_{ee}$ and torque $\tau _{ee}$ ) at the end-effector to ensure physically plausible interactions:

(24) \begin{equation} \left\{\begin{array}{l} u_{ee,\lambda }=f_{ee}\\[6pt] u_{ee,\tau }=\tau _{ee} \end{array}\right. \end{equation}

Since the cost function already includes an optimization objective for the input terms, we employ a relaxed barrier function method to convert these hard constraints into soft constraints integrated into the cost function. Specifically, the composite cost function is constructed as:

(25) \begin{equation} L\left(x\left(t\right),u\left(t\right),t\right)=\left\| x_{r}-x_{r}^{ref}\right\| _{Q_{r}}^{2}+\left\| u-u^{ref}\right\| _{R}^{2}+\left\| r_{IE}-r_{IE}^{ref}\right\| _{Q_{eep}}^{2}+\left\| \xi _{IE}\right\| _{Q_{eeo}}^{2} \end{equation}

where $u^{ref}=[\begin{array}{lll} 0_{3{n_{c}}} & \lambda _{ee} & 0_{{n_{a}}} \end{array}]^{T}\in \mathrm{\mathbb{R}}^{3{n_{c}}+6+{n_{a}}}$ is the extended input vector.

5.1. Force estimation performance validation

As shown in Figure 3, during this experiment, the end-effector grasped a 1 kg counterweight while the robot maintained stationary standing. Based on local gravitational acceleration (9.8 m/s²), the nominal external force acting on the robot was 9.8 N.

Comparative analysis of the z-axis force estimation results (Figure 3) reveals that both the traditional method and the α-adaptive method converge near the theoretical value, confirming algorithmic correctness. Notably, the α-adaptive method demonstrates faster convergence compared to traditional methods. These improvements primarily stem from its innovative dynamic forgetting factor adjustment mechanism and deep coupling design with the WBC controller. The performance enhancement significantly boosts the robot’s adaptability in unstructured environments, providing reliable technical support for subsequent complex interaction tasks. Experimental data fully validate the dual advantages of the proposed method in both response speed and estimation accuracy, with steady-state estimation mean values deviating $\lt$ 5% from theoretical values.

Furthermore, as detailed in Table I, we conducted a comparative analysis of the computational time consumption for the WBC control loop with the external force estimation module. Experimental results demonstrate that after integrating the force estimation algorithm, the average computation time per control cycle increased by only 2.5 × 10^-5 ms. This marginal increment remains entirely within the permissible range of the 1 ms control period, ensuring the system’s capability to stably maintain 1000 Hz high-frequency control.

Figure 3.

A performance comparison between the conventional method and the α-adaptive method.

5.2. Algorithm comparison

To systematically evaluate the performance of the proposed adaptive control algorithm, comparative experiments were conducted with the non-adaptive method from ref. [Reference Sleiman, Farshidian, Minniti and Hutter21]. The experimental setup involved a mobile manipulation platform performing in-place stepping motions while monitoring end-effector positional changes under time-varying step loads (20N, 25N, 30N, 35N) applied at the end-effector, with all tests completed in Gazebo simulation environment for safety considerations.

Figure 4.

Compare the postural changes of the robot under adaptive and non-adaptive control.

Figure 5.

Comparing performance of adaptive and nonadaptive controllers.

Experimental results (Figures 4, 5) demonstrate the proposed algorithm’s superior performance: in terms of robustness, it maintains excellent trajectory tracking performance under full load conditions with end-effector position deviations controlled within ±1cm; regarding stability, it successfully handles sudden 35N loads while maintaining stable system operation; comparative algorithm performance shows the non-adaptive method exhibited significant Z-axis deviation ( $\gt$ 5 cm) at just 20N load, ultimately leading to system instability.

Figure 6.

Successful transportation to target location (0.8 m,1.0 m).

5.3. Trajectory tracking under time-varying forces

This experiment validated the effectiveness of the proposed adaptive control algorithm in handling force uncertainties through a challenging non-homogeneous object manipulation task. As shown in Figure 6, the quadruped mobile manipulation platform successfully completed grasping and transporting a 2.5 kg non-homogeneous object (with internal rolling balls simulating mass distribution variations).

The results demonstrate that the proposed algorithm maintains excellent control performance under dynamic disturbances, achieving fast and precise positioning at the target location. Notably, the system exhibited stable transportation capability even under harsh oscillatory disturbances caused by suspended weights.

Experimental results demonstrate that the proposed algorithm effectively handles uncertainties in unstructured environments, as evidenced by comprehensive experimental records, thereby providing a reliable technical solution for complex manipulation tasks. Our work builds upon the advanced framework introduced by Jean-Pierre Sleiman et al., inheriting its advantage of a large manipulator workspace. The key contribution of this study lies in overcoming the stability limitation of that framework when dealing with unknown payloads. Our method achieves a trajectory tracking accuracy within ±2 cm, demonstrating its potential for practical applications in complex interactive scenarios such as warehouse logistics or search-and-rescue operations.