Legged humanoid robots need to deliver stable and robust legged manipulation (also referred to as loco-manipulation) behaviors while operating in unstructured environments. For this purpose, Whole-Body Control (WBC) enables the generation of control commands for tracking multiple task trajectories effectively Khatib et al. (2004). Most WBC frameworks rely on dynamic models of the robot and classical control laws such as proportional-derivative (PD) control Feng et al. (2015) or impedance control Albu-Schäffer et al. (2007). Based on these requirements, projection-based WBC approaches have been widely used due to their intuitive, concise, and computationally efficient form. However, more versatile algorithms are frequently needed to consider inequality and unilateral constraints Lee et al. (2012) such as Contact Wrench Cone (CWC) constraints Caron et al. (2015). Optimization-based approaches have been widely used to incorporate inequality constraints explicitly Escande et al. (2014); Hong et al. (2020); Feng et al. (2014); Wiedebach et al. (2016); Kim et al. (2020). However, there is a major challenge for these optimization-based WBC approaches: how can we ensure the desired task hierarchy and how can we guarantee system stability when there exist uncertain external forces?

A well-known method to enforce task priorities employs hierarchical quadratic programming (HQP) Escande et al. (2014); Hong et al. (2020). In Feng et al. (2014); Wiedebach et al. (2016), the authors employ a single quadratic program (QP) with different weighted cost terms to impose a task hierarchy. These types of controllers that employ a single QP process reduce the number of optimizations compared with HQP and enforce the task hierarchy in an implicit manner; however, the weights are heuristic, resulting in occasional violations of the priorities. Often, there is a discrepancy between the desired configuration and the result from the WBC optimization. The joint acceleration, which is one of the results from the WBC optimization process, can be numerically integrated to resolve this problem Kuindersma et al. (2016); Koolen et al. (2016); Ahn et al. (2021). Alternatively, this problem can be solved by directly re-optimizing the trajectories considering dynamic reachability Lee et al. (2020). Another approach is to have WBC incorporate both internal and external force/torque measurements via embedded sensors within the mechanical actuation/transmission system Nori et al. (2015). These sensors can be used to improve the stability of the robots via their WBC controllers.

In general, it is complex to analyze the stability of robots controlled by the WBC approaches. For instance, stability of a priority-based kinematic control approach is verified based on Lyapunov stability in Antonelli (2009). Also, it is shown that the robot controlled by the operational space control framework is asymptotically stable Dietrich et al. (2018). However, the above stability analysis studies do not consider floating-base dynamics of the legged robot and the presence of unknown external disturbances. Compared with the above stability analysis which are applied to projection-based approaches, the stability analysis of optimization-based WBC approaches is more complicated to perform when considering external perturbations. We previously proposed an optimization-based WBC method dubbed WBLC, that combines both a projection-based controller and an optimization technique to enforce task priorities while satisfying inequality constraints Kim et al. (2020). The WBLC framework combines a whole-body kinematic controller and a whole-body dynamic controller to compute joint position, velocity, acceleration, and torque references in real-time as shown in Figure 1. WBLC has successfully demonstrated dynamic locomotion of passive ankled robots. However, the robustness of WBLC regarding the external disturbances mostly relies on the experts’ parameter tuning in experiments. Also, we did not address the theoretic stability analysis in our previous work Kim et al. (2020). In this paper, we will analyze here our WBLC framework for ensuring stability and tracking performance under external forces. This is key since if WBLC is affected by unknown external forces, the robot could easily lose balance.

Compared to other methods which focus on push recovery, here we focus on adaptive control as a means to deal with external forces. Some push recovery methods employ linear MPC to obtain CoM and foot trajectories based on the Zero Moment Point (ZMP) dynamics to re-plan walking behaviors for disturbance recovery Mason et al. (2018). The desired recovery force associated with the ZMP condition is converted to CoM positions Wang et al. (2014). It has been widely investigated how to determine stepping motions for balance recovery such as Push Recovery Model Predictive Control (PR-MPC) Stephens and Atkeson (2010b) and Capture Point control Pratt et al. (2006). Although these strategies have accomplished successful results, they are based on motion generation instead of modifying feedback gains or impedances. However, the latter strategy considered in this paper can provide a layer of additional robustness. Indeed, it is not always needed to perform an additional stepping motion to recover from external forces. If the external or contact forces are not adequately compensated for or there exist unknown external forces, the controllers may be unstable regardless of the recovery strategy Nakanishi et al. (2008).

One solution to this problem is to compensate for external forces. This has been often studied as a compliant control problem for balancing despite the external forces Hyon et al. (2007); Stephens and Atkeson (2010a); Ott et al. (2011); Herzog et al. (2016). A key idea of the above researches is to obtain proper ground reaction forces keeping the balance against the external force then the forces are utilized to compute the control command torque Hyon et al. (2007); Ott et al. (2011). LQR is employed for better momentum control of a torque-controlled humanoid and they show robust performance on balancing in the face of unknown disturbances Herzog et al. (2016). A common point of the above researches is that they rely on constant gains or impedances. Although the above researches have accomplished practical balancing ability, the stability of the controlled robots are not proven at the actuation level yet when the external forces are applied to the robots. Furthermore, they also need to modify or replan the task specifications such CoM and angular momentum for keeping the robot’s balance. By contrast, our work does not require this replanning step.

Another solution to the above problem is to directly measure external forces/moments using an F/T sensor Käslin et al. (2018) or internal, tactile, and F/T sensors, simultaneously Ivaldi et al. (2011). However, it is normally challenging to attach F/T sensors everywhere in the robot’s body. An alternative method is to accurately estimate external forces and the corresponding ground reaction forces relying on the dynamic model of the robot Xin et al. (2018). For humanoid robotic applications, an external force/moment can be estimated based on an extended Kalman filter using whole-body dynamics, force sensor on the feet, and an IMU Benallegue et al. (2018). In some cases, external forces can be measured using force-sensing resistors Hawley et al. (2019) on the feet. We can also deploy adaptive control approaches for stabilizing the robots without direct estimation of the external forces. A robust and adaptive control approach was proposed without using the robot’s dynamic model to guarantee stable control performance under uncertain disturbances Lee et al. (2016). This approach shows the robustness of WBC with respect to impact disturbance at the task level.

Overall, this paper proposes a unique gain adaptation approach for WBLC applied to legged robots under both contact and external forces. More specifically, we aim to enhance the robot’s balance and robust performance by guaranteeing feedback system stability under unknown external forces. To do that, we leverage the WBLC as proposed in Kim et al. (2020). Then, we approximately decouple the joint-space dynamics and analyze their closed-loop stability. Compared to other approaches such as Lee et al. (2016), we adapt the feedback gains to make the closed-loop system critically damped at the joint level. This results in four contributions.

• Firstly, the proposed adaptation approach is very intuitive and straightforward to be implemented. Our method does not require any modification of the task specifications defined a priori. Without any significant modifications to the controller and planners, it is also possible to employ classical WBC methods Kim et al. (2020) for controlling the robot in the presence of unknown external forces.

In this paper, our primary task is to control the position of the CoM, keeping balance in the double support phase and to stably walk against the external disturbances. Many studies on legged robots such as push recovery and reactive planning have addressed keeping balance while standing and walking. Regarding reactive planning processes, they do not change low-level controller feedback gains as we do and instead modify the desired task specifications or trajectories in the presence of the perturbations. Therefore, in such methods, it is hard to guarantee the closed-loop stability of the robot while tracking the desired task trajectories. By contrast, our method stabilizes the robot against the external disturbance by adapting the feedback gains without further modification of task trajectories. For this reason, we can ensure the closed-loop stability of the robot while the external disturbance is applied to the robot. Our proposed approach is therefore simpler and more intuitive than reactive planning approaches and it can also be used in combination with them to enhance their performance.

The remainder of this paper is organized as follows. Section 2 explains the preliminaries related to a QP-based WBC approach, e.g., WBLC and the Centroidal Dynamic model. The proposed approach is explained in Section 3. In Section 4, we apply the proposed adaptive approach to a legged robot Draco2 to show the effectiveness of the proposed method in numerical simulations and hardware experiments.

We review here the structure of our previous WBC controller. Given a configuration space Q ⊆ R n and an input space U ⊆ R n u , the rigid-body dynamics of a robot is described as follows: M q q ̈ + B q , q ̇ = S ⊤ Γ + J c ⊤ q F c (1) Where q ∈ Q , M ( q ) ∈ S + + n , B ( q , q ̇ ) ∈ R n , and Γ ∈ U denote the joint position vector, mass/inertia matrix, sum of Coriolis/Centrifugal and gravity forces, and torque command, respectively. Considering six DOFs virtual joints for the floating base, n _u = n − 6, we can define the selection matrix S = 0 n u × 6 I n u × n u ∈ R n u × n . F c ∈ R 6 m is a vertically concatenated contact force vector and the corresponding contact Jacobian is represented as J c ( q ) ∈ R 6 m × n where m is the number of contacts. A standard feedback control scheme can be designed in the joint space using a PD control law: Γ cmd = S M q K p e q + K d e ̇ q + Γ ̂ d (2) Where e _q = q ^d − q, and e ̇ q = q ̇ d − q ̇ . Also, q d ∈ Q , q ̇ d ∈ R n , K p = diag ( k p 1 , … , k p n ) and K d = diag ( k d 1 , … , k d n ) denote desired joint position, desired joint velocity, proportional gain matrix and derivative gain matrix, respectively. Γ ̂ d represents a feedforward torque command generated by the Dynamic-level WBC control block shown in Figure 1. Using a hierarchical inverse kinematic control scheme, we can compute q ^d and q ̇ d given task specifications, x ^d , x ̇ d , and x ̈ d . Let us consider N hierarchical tasks in a lexicographic order. The desired joint position is obtained by q ^d = q + Δq where Δ q = ∑ k = 1 N Δ q k and Δ q k = J k | k − 1 q + e x k − J k q Δ q k − 1 , N k q = N k − 1 q − J k | k − 1 q + J k | k − 1 q (3) Where J _k|k−1(q) = J _k (q)N _k−1(q), e x k = x k d − x k and N ₀(q) = I. In addition, q ̇ d and q ̈ d can be simply computed as q ̇ d = ∑ k = 1 N q ̇ k d and q ̈ d = ∑ k = 1 N q ̈ k d where q ̇ k d = J k | k − 1 q + x ̇ k d − J k q q ̇ k − 1 d , q ̈ k d = J k | k − 1 q + x ̈ k d − J ̇ k q , q ̇ q ̇ − J k q q ̈ k − 1 d , (4) q ̇ 0 d = 0 , and q ̈ 0 = 0 . Based on derivations of q ^d , q ̇ d , and q ̈ d , an additional optimization problem is formulated to handle full-body dynamics of the robot, contact wrench cone constraints, and torque limits as described in Kim et al. (2020), i.e.: min F c , x ̈ c , δ q ̈ F c ⊤ W F F c + x ̈ c ⊤ W c x ̈ c + δ q ̈ ⊤ W q ̈ δ q ̈ s . t . M q q ̈ + B q , q ̇ = S ⊤ Γ + J c ⊤ q F c , x ̈ c = J c q q ̈ + J ̇ c q , q ̇ q ̇ , q ̈ = q ̈ d + K p h e q + K d h e ̇ q + δ q ̈ , U q F c ≥ 0 , Γ min ≤ Γ ≤ Γ max (5) Where W F ∈ S + 6 m , W c ∈ S + 6 m , and W q ̈ ∈ S + n are weighting matrices for the objective function. K p h = diag ( k p 1 h , … , k p n h ) and K d h = diag ( k d 1 h , … , k d n h ) are feedback gains in the optimization problem. U ( q ) ∈ R 17 is a matrix for expressing the contact wrench cone constraints and upper bounds of the reaction forces for smooth contact changes. Γ _min and Γ _max are the minimum and maximum torques of the actuators. Based on the optimal decision variables F c * and δ q ̈ ⋆ , we can provide the feedforward torque command as follows: Γ ̂ d = S M q q ̈ ⋆ + B q , q ̇ − J c q ⊤ F c ⋆ (6) Where q ̈ ⋆ = q ̈ d + K p h e q + K d h e ̇ q + δ q ̈ ⋆ . Using the above desired torque command Γ ̂ d , we can compute the low-level joint torque command in (2).

We analyze the closed-loop behavior of the robotic system. Substituting (6) into (2) we get: Γ cmd = S M q q ̈ ̂ + B q , q ̇ − J c q ⊤ F c ⋆ (7) Where q ̈ ̂ = q ̈ d + K ̂ p e q + K ̂ d e ̇ q + δ q ̈ ⋆ , K ̂ p = ( K p h + K p ) , and K ̂ d = ( K d h + K d ) . By substituting the command in (7) into the full-body dynamics in (1), the closed-loop dynamics become: M v v M v a M v a ⊤ M a a q ̈ v q ̈ a = d v d a , (8) With d v = S v J c q ⊤ F c − B v , d a = M v a ⊤ q ̈ ̂ v + M a a q ̈ ̂ a + S J c q ⊤ F c − F c ⋆ (9) Where subscriptions ( . ) v ∈ R 6 and ( . ) a ∈ R n u denote the properties for virtual and actuated joints, respectively. Virtual joints are those corresponding to degrees of freedom for the floating base whereas actuated joints correspond to motorized robot articulations. As such, we define a selection matrix for the virtual joints, S v = I 6 × 6 0 6 × n u ∈ R 6 × n . We also introduce a matrix called the Schur Complement of M _aa in M(q) as M a a , s = M v v − M v a M a a − 1 M v a ⊤ . Since M(q) is positive definite and M _aa is invertible, the matrix M _aa,s is also invertible. The actuated joint acceleration can be expressed as follows: q ̈ a = M a a − 1 d a + M a a − 1 M v a ⊤ M a a , s − 1 M v a M a a − 1 d a − d v = q ̈ ̂ a + Δ c (10) Where Δ c = M a a − 1 M v a ⊤ ( q ̈ ̂ v − q ̈ v ) + S J c ( q ) ⊤ ( F c − F c ⋆ ) . Assuming that the optimal decision variable satisfies F c = F c ⋆ and q ̈ ̂ v = q ̈ v with perfect compensation for the Coriolis/centrifugal and gravitational forces, we can decouple the joint space dynamics: Δc = 0. However, the optimal decision variable F c ⋆ relies on the robot model and there is no direct feedback control-loop for regulating the contact force error F c − F c ⋆ . Also, Δc becomes significant when Centroidal dynamics are not precisely controlled due to external disturbances, since both q ̈ v and F _c significantly affect the Centroidal dynamics of the robot. So the above assumptions are frequently invalid in the real world. This is the reason to introduce the error Δc in (10). Then, we can approximately decompose the closed-loop dynamics in the joint space as follows: q ̈ i = q ̈ i f + k ̂ p i q i d − q i + k ̂ d i q ̇ i d − q ̇ i + Δ c i (11) Where q ̈ i f = q ̈ i d + δ q ̈ i ⋆ , k ̂ p i = k p i + k p i h , and k ̂ d i = k d i + k d i h for all i ∈ {7, …, n}. q _i , q ̇ i , and q ̈ i represent the joint position, velocity, and acceleration values for the ith joint, respectively. Also, Δc _i is the ith component of Δc. The relationships between the natural frequency, W _n , and the PD gains are defined as W n 2 = k ̂ p i and 2 ζ W n = k ̂ d , respectively, where ζ denotes the damping ratio. To achieve critical damping, the derivative gain must be, k ̂ d i = 2 k ̂ p i . Based on this analysis, it is possible to obtain the gain matrices of the whole-body control to achieve stable closed-loop behavior.

Next, we consider an additional unknown external force applied to the robot’s body, such as a force acting on the torso or pelvis. The optimization problem (5) cannot directly incorporate the unknown perturbation. Given such perturbations, the rigid-body dynamics of the robot changes as follows: M q q ̈ + B q , q ̇ = S ⊤ Γ + J c ⊤ q F c + J t ⊤ q F t (12) Where F t ∈ R 6 denotes the unknown external force acting on the robot’s body and J t ( q ) ∈ R 6 × n is the corresponding Jacobian. Based on Proposition 1, the closed-loop behavior of the system with the control torque input becomes: q ̈ a = q ̈ ̂ a + Δ c + M a a − 1 S J t q ⊤ F t = q ̈ ̂ a + Δ c + S + E m M q − 1 S ⊤ S J t q ⊤ F t ︸ S ψ t (13) Where E m = − M a a , s − 1 M v a ⊤ M v v − 1 M v v , s 0 . When the external disturbance is unknown, it is impossible to compute Δc and F _t , separately. Therefore, we consider the effect of the external disturbance using Sψ _t .

Proposition 1. Given a matrix E = E 1 E 12 E 12 ⊤ E 2 ∈ S + + n where E 1 ∈ R n 1 × n 1 , E 2 ∈ R n 2 × n 2 , and n = n ₁ + n ₂ , E 2 − 1 can be expressed as ( E + S 2 ) G S 2 ⊤ where G = E ⁻¹ and S 2 = 0 n 2 × n 1 I n 2 × n 2 ∈ R n 2 × n .

The detailed proof of Proposition 1 is provided in Appendix. We introduce an estimated external force F ̂ t = [ f ̂ t ⊤ , τ ̂ t ⊤ ] ⊤ satisfying the following equality: flushleft S ψ t = S M q − 1 Σ a J t q ⊤ F ̂ t (14) Where f ̂ t ∈ R 3 and τ ̂ t ∈ R 3 denote the estimated force and torque of the unknown disturbance, respectively. Also, Σ a = S ⊤ S = 0 6 × 6 0 6 × n u 0 n u × 6 I n u × n u . Unfortunately, we cannot simply express Sψ _t in terms of q and q ̇ since we need an additional explicit mapping between F ̂ t and these joint state variables.

Our approach is to employ robot centroidal dynamics to convert the external forces into corresponding stiffness and damping forces in the joint space. We can approximate the CoM velocity, acceleration, and the time derivative of its angular momentum with a discrete time interval Δt Dai et al. (2014) as follows: p ̇ ≈ 2 Δ t − 1 p − p p − J CoM q p q ̇ p , p ̈ ≈ Δ t − 1 p ̇ − J CoM q p q ̇ p , h ̇ q , q ̇ ≈ Δ t − 1 h q , q ̇ − h q p , q ̇ p (15) Where p, h ( q , q ̇ ) , J _CoM(q) denote the position of CoM, its angular momentum, and its Jacobian, respectively. q ^p , q ̇ p , p ^p , and p ̇ p are the joint position, the joint velocity, the CoM position, and the CoM velocity in the previous time step, respectively. First, we estimate the linear force corresponding to the external disturbance, f _t , using the centroidal dynamics and CoM acceleration: f ̂ t = M r p ̈ − g − ∑ ℓ = 1 m f ℓ ⋆ (16) Where f ℓ ⋆ ∈ R 3 represents the linear force part for the ℓ-th contact among F c ⋆ from WBC. Also, M _r and g denote the total mass of the robot and the gravity vector, respectively. In addition, τ ̂ t can be estimated by using the centroidal momentum equation and f ̂ t . τ ̂ t = h ̇ q , q ̇ − ∑ ℓ = 1 m r ℓ × f ℓ ⋆ + τ ℓ ⋆ − r t × f ̂ t (17) Where τ ℓ ⋆ ∈ R 3 denotes the torque part for the ℓ-th contact among F c ⋆ . In turn, we can re-express the closed-loop dynamics (13) using the estimated external force F ̂ t . We cannot decouple the closed-loop dynamics due to the term ψ _t .

Our proposed method consists of decoupling the closed-loop dynamics by approximating ψ _t . Our goal is to obtain a pair of PD gains to ensure that the approximated closed-loop dynamics are critically damped, i.e. ψ t ≈ K ̃ p q + K ̃ d q ̇ + Ω t (18) Where K ̃ p = diag ( k ̃ p 1 , … , k ̃ p n ) , K ̃ d = diag ( k ̃ d 1 , … , k ̃ d n ) , and Ω _t are approximated PD gain matrices and a constant bias term, respectively. To obtain this approximation we have used the centroidal dynamics and the external force approximations specified further above. First, p is calculated using a first order approximation: p ≈ p p + J CoM q p q − q p (19) In turn the estimated external force f ̂ t can be expressed in joint space as follows using (15) and (16): f ̂ t q , q ̇ = M r J CoM q p L f q , q ̇ − M r g − ∑ ℓ = 1 m f ℓ ⋆ , L f q , q ̇ = 2 q − q p Δ t 2 − 2 q ̇ p Δ t = K p f q + K d f q ̇ + ϵ f (20) Where both matrices, K p f = diag ( k p 1 f , … , k p n f ) and K d f = diag ( k d 1 f , … , k d n f ) , represent gain matrices. Using established adaptive techniques we chose the proportional gain as k p i f = η i f ( q i d − q i p ) and plugging it in the second equation shown in (20), we get k d i f = ξ i f Δ t k p i f where η i f and ξ i f represent positive coefficients. We also get ϵ _f as: ϵ f = 2 I − Δ t K p d f q ̈ f − K p d f q ̇ p − K p f q p (21) Where K p d f = Δ t K p f + K d f . In turn Eq. 20 becomes: f ̂ t q , q ̇ = K ̃ p f q + K ̃ d f q ̇ + γ f (22) Where K ̃ p f = M r J CoM ( q p ) K p f , K ̃ d f = M r J CoM ( q p ) K d f , and γ f = M r J CoM ( q p ) ϵ f − M r g − ∑ ℓ = 1 m f ℓ ⋆ .

Secondly, we estimate (15) using the first order approximation: h q , q ̇ ≈ h q p , q ̇ p + H G q p , q ̇ p q − q p + A G q p q ̇ − q ̇ p , A G q = ∂ h ∂ q ̇ q , q ̇ , H G q , q ̇ = ∑ k = 1 n ∂ A G k q ∂ q q ̇ k (23) Where A G k ( q ) is the kth column vector of A _G (q). Based on the above approximations, the estimated external torque can be expressed as follows: τ ̂ t q , q ̇ = H G q p , q ̇ p L τ q , q ̇ + A G q p q ̈ f − ∑ ℓ = 1 m r ℓ × f ℓ ⋆ + τ ℓ ⋆ − r t × f ̂ t q , q ̇ , L τ q , q ̇ = K p τ q + K d τ q ̇ + ϵ τ (24) Where both matrices, K p τ = diag ( k p 1 τ , … , k p n τ ) and K d τ = diag ( k d 1 τ , … , k d n τ ) , are gain matrices for estimating the external torque. More specifically, each component can be computed using adaptive control techniques in Åström and Wittenmark (2013) as: k p i τ = η i τ ( q i d − q i p ) and k d i τ = ξ i τ Δ t k p i τ , where η i τ and ξ i τ represent positive coefficients. The term ϵ _τ is described as ϵ τ = Δ t I − Δ t K p d τ q ̈ f + I − K p d τ q ̇ p − K p τ q p where K p d τ = Δ t K p τ + K d τ . Suppose r _ℓ and r _t are constant then we can re-write the estimated τ ̂ t as: τ ̂ t q , q ̇ = K ̃ p τ q + K ̃ d τ q ̇ + γ τ (25) Where K ̃ p τ = H G ( q p , q ̇ p ) K p τ − s k e w ( r t ) K ̃ p f and K ̃ d τ = H G ( q p , q ̇ p ) K d τ − s k e w ( r t ) K ̃ d f . The detailed γ _τ is γ τ = H G q p , q ̇ p ϵ τ − r t × γ f − ∑ ℓ = 1 m r ℓ × f ℓ ⋆ + τ ℓ ⋆ + A G q p q ̈ f . (26) We can construct F ̂ t using selection matrices S f = I 3 × 3 , 0 3 × 3 ⊤ and S τ = 0 3 × 3 , I 3 × 3 ⊤ as follows: F ̂ t = S f f ̂ t q , q ̇ + S τ τ ̂ t q , q ̇ = K ̃ p s q + K ̃ d s q ̇ + γ s (27) Where K ̃ p s = S f K ̃ p f + S τ K ̃ p τ , K ̃ d s = S f K ̃ d f + S τ K ̃ d τ , and γ ^s = S _f γ _f + S _τ γ _τ . Using the above estimated external force F ̂ t , we approximate ψ _t as follows: ψ t ≈ M q p − 1 Σ a J t q p K ̃ p s q + K ̃ d s q ̇ + γ s ≈ Z ̂ p q + Z ̂ d q ̇ + Ω ̂ t (28) Where Z ̂ p and Z ̂ d are diagonal matrices. We define Z ̂ p = diag ( z p 1 , … , z p n ) and Z ̂ d = diag ( z d 1 , … , z d n ) where z p i and z d i are the ith diagonal components of M ( q p ) − 1 Σ a J t ( q p ) K ̃ p s and M ( q p ) − 1 Σ a J t ( q p ) K ̃ d s , respectively. The last term takes the form Ω ̂ t = M ( q s ) − 1 J t ( q p ) ⊤ γ s . Since K ̂ p , K ̂ d , Z ̂ p , and Z ̂ d are diagonal matrices, the closed-loop dynamics in actuated joint level space becomes: e ̈ a = K ̂ p a − Z ̂ p a e a + K ̂ d a − Z ̂ d a e ̇ a + Ω ̄ t , a (29) Where e q a = q a d − q a , e ̇ a = q ̇ a d − q ̇ a , and e ̈ a = S ( q ̈ d + δ q ̈ ⋆ ) − q ̈ a . Also, we define K ̂ p a = S K ̂ p = diag ( k ̂ p 7 , … , k ̂ p n ) , K ̂ d a = S K ̂ d = diag ( k ̂ d 7 , … , k ̂ d n ) , Z ̂ p a = diag ( z p 7 , … , z p n ) , and Z ̂ d a = diag ( z d 7 , … , z d n ) , and Ω ̄ t , a = S ( Ω ̂ t − Z ̂ p q d − Z ̂ d q ̇ d ) .

We define the state ν = [ e a ⊤ , e ̇ a ⊤ ] ⊤ and express the state-space model of the above dynamics equation as follows: ν ̇ = 0 n I n − K ̂ p a + Z ̂ p a − K ̂ d a + Z ̂ d a ︸ V ν + 0 n × 1 Ω ̄ t , a . (30) Properly choosing η i f , ξ i f , η i τ , and ξ i τ , we can make Ω ̄ t , a negligible. Then, the matrix V should be Hurwitz to guarantee that the closed-loop state model is exponentially stable. To compute the characteristic equation of the matrix V, we express the determinant of the matrix V − λ I _2n as follows: det − λ I n I n − K ̂ p a + Z ̂ p a − K ̂ d a + Z ̂ d a − λ I n = λ − 1 det K ̂ d − Z ̂ q ̇ + λ − 1 K ̂ p − Z ̂ q + λ I n × n = det λ − 1 K ̂ d a − Z ̂ d a + λ − 1 K ̂ p a − Z ̂ p a ︸ X + I n , Which means that   det(V − λ I ₂) =  det(λ ⁻¹ X + I _n ). It is complicated to obtain analytic forms of the eigenvalues of V. For this reason, we can compute the determinant above using the Leibniz formula described in Bhatia (2013). det λ − 1 X + I n = 1 + λ − 1 t r a c e X + O λ − 1 2 (31) Where O ( . ) represents the big-O notation. Using the above formulation, the characteristic equation can be re-written as: det V − λ I 2 n = 1 + λ − 1 t r a c e X + O λ − 1 2 = 0 . (32) All solutions of the above characteristic equation λ should be negative to guarantee that the matrix V is Hurwitz. Such analysis could be conducted via numerical computation given the gain matrices. However, in this paper, we aim at adaptive feedback gains of the controller so that we obtain K ̃ p a = Z ̂ p a and K ̃ d a = Z ̂ d a . In turn the term O ( ( λ − 1 ) 2 ) vanishes. We can obtain the analytical expression of the approximated characteristic equation as follows: λ 2 + λ t r a c e K ̂ p a − K ̃ p a + t r a c e K ̂ d a − K ̃ d a = 0 . (33) This condition will allow us to make the closed-loop system stable.

Condition 1. The following condition makes the matrix V Hurwitz: ∑ i = 7 n k ̂ p i − z p i > 0 and ∑ i = 7 n k ̂ d i − z d i > 0 . Without any adaptation, these condition should be fulfilled to make the closed-loop system asymptotically stable given desired feedback gains.

We can approximately decouple the dynamics of (29) in the joint space as follows: q ̈ i + k ̂ d i − z d i q ̇ i + k ̂ p i − z p i q i = F ̃ i ⋆ (34) Where F ̃ i ⋆ = q ̈ i f + k ̂ p i q i d + k ̂ d i q ̇ i + Ω ̂ t ( i ) . We enforce that the joint space dynamics are critically damped given z p i and z d i for all i ∈ {7, …, n}. For the above system, the following relationships are obtained in terms of the natural frequency and the damping ratio: k ̂ p i − z p i = W ̄ n 2 and k ̂ d i − z d i = 2 ζ W ̄ n . Since we already know all gains and coefficients, k ̂ p i , k ̂ d i , z p i , and z d i , ζ can be determined a priori. We can compute adaptive variables ε p i and ε d i to make W ̄ n = W n d and ζ = ζ ^d . ε p i = W n d 2 + z q i − k ̂ p i , ε d i = 2 ζ d k ̂ p i + ε p i − z q i + z d i − k ̂ d i . (35)