Towards Exact Interaction Force Control for Underactuated Quadrupedal Systems with Orthogonal Projection and Quadratic Programming

The dynamics equation of a quadruped robot can be expressed as:

$$\bm{M}\ddot{\bm{q}}+\bm{h}=\bm{B}\bm{\tau}+\bm{J}_{c}^{T}\bm{\lambda}+\bm{J}_{x}^{T}\bm{F}_{x},$$

(1)

where $\bm{q}=\left[\bm{q}_{b}^{T},\,\bm{q}_{j}^{T}\right]^{T}$ is the generalized coordinate vector including unactuated floating base coordinates $\bm{q}_{b}\in SE(3)$ and actuated joint configuration $\bm{q}_{j}\in\mathbb{R}^{n_{j}}$, $\bm{M}\in\mathbb{R}^{(6+n_{j})\times(6+n_{j})}$ denotes the inertia matrix, $\bm{h}\in\mathbb{R}^{(6+n_{j})}$ is the non-linear effect consisting of Coriolis, centrifugal and gravitational forces, $\bm{B}=\left[\bm{0}_{n_{j}\times 6},\,\bm{I}_{n_{j}\times n_{j}}\right]^{T}$ is the selection matrix¹¹1Unlike the selection matrix in previous papers, the selection matrix here is not square. of actuated joints, $\bm{\tau}\in\mathbb{R}^{n_{j}}$ is the actuated joint torques, $\bm{J}_{c}\in\mathbb{R}^{k\times(6+n_{j})}$ represents the constraint Jacobian with $k=3n_{c}$ ($n_{c}$ denotes the number of legs in contact with the environment and the contact is assumed as a point contact), $\bm{\lambda}\in\mathbb{R}^{k}$ denotes the generalized constraint force vector used to control unactuated $\bm{q}_{b}$, $\bm{J}_{x}\in\mathbb{R}^{n_{e}\times(6+n_{j})}$ is the Jacobian at $\bm{x}$, and $\bm{F}_{x}\in\mathbb{R}^{n_{e}}$ is the external disturbances due to the interaction from human or environments.

During locomotion, the support feet should not slip, i.e., the constraint $\bm{J}_{c}\dot{\bm{q}}=\bm{0}$ must be satisfied. This constraint indicates that any admissible $\dot{\bm{q}}$ lies in the constraint null space $\mathcal{N}\left(\bm{J}_{c}\right)$. According to [10], the dynamics equation $\eqref{dynEqa}$ can be projected into two subspaces by using the orthogonal projection matrix $\bm{P}$ and $\bm{I}-\bm{P}$ respectively as:

$$\bm{PM}\ddot{\bm{q}}+\bm{Ph}=\underbrace{\bm{PB\tau}}_{\bm{\tau}_{m}}+\bm{P}\bm{J}_{x}^{T}\bm{F}_{x},$$

(2)

$$(\bm{I}-\bm{P})(\bm{M}\ddot{\bm{q}}+\bm{h})=\underbrace{(\bm{I}-\bm{P})\bm{B}\bm{\tau}}_{\bm{\tau}_{c}}+\bm{J}_{c}^{T}\bm{\lambda}+(\bm{I}-\bm{P})\bm{J}_{x}^{T}\bm{F}_{x},$$

(3)

where $\bm{P}=\bm{I}-\bm{J}_{c}^{+}\bm{J}_{c}$ implies that the orthogonal null space projection matrix is computed from the Moore-Penrose pseudoinverse of $\bm{J}_{c}$ and projects vectors into the null space of the constraint, such that $\bm{P}=\bm{P}^{2}=\bm{P}^{T}$, $\bm{P}\bm{J}_{c}^{T}=\bm{0}$, and $\bm{P}\dot{\bm{q}}=\dot{\bm{q}}$ for all $\dot{\bm{q}}\in\mathcal{N}\left(\bm{J}_{c}\right)$. Then $\bm{I}-\bm{P}$ represents the complementary projection into $\mathcal{N}^{\perp}(\bm{J}_{c})$.

To solve the equation of $\bm{\lambda}$, $\ddot{\bm{q}}$ must be computed first. However, $\ddot{\bm{q}}$ cannot be obtained directly through pre-multiplying the inverse of $\bm{PM}$ in $\eqref{dynInP}$, as $\bm{PM}$ cannot be inverted attributed to the rank deficiency of $\bm{P}$. With the help of additional equations $(\bm{I}-\bm{P})\dot{\bm{q}}=\bm{0}$ and its derivative $(\bm{I}-\bm{P})\ddot{\bm{q}}=\dot{\bm{P}}\dot{\bm{q}}$ that are derived from $\bm{P}\dot{\bm{q}}=\dot{\bm{q}}$, $\ddot{\bm{q}}$ can be solved by substituting the latter equation into $\eqref{dynInP}$ as:

$$\ddot{\bm{q}}=\bm{M}_{c}^{-1}(\bm{\tau}_{m}-\bm{Ph}+\dot{\bm{P}}\dot{\bm{q}}+\bm{P}\bm{J}_{x}^{T}\bm{F}_{x}),$$

(4)

where $\bm{M}_{c}=\bm{PM}+\bm{I}-\bm{P}$. Equations (4) shows that only the motion torques $\bm{\tau}_{m}$ contributes to the motion of the system, therefore $\mathcal{N}\left(\bm{J}_{c}\right)$ is called the motion space as described in [6]. Similarly, because in $\eqref{dynInIMinusP}$ the $\ddot{\bm{q}}$ is fixed and constraint forces $\bm{\lambda}$ are free to choose depending on the constraint torques $\bm{\tau}_{c}$, $\mathcal{N}^{\perp}(\bm{J}_{c})$ is named the constraint space. Eventually, the constraint forces are obtained by inserting $\ddot{\bm{q}}$ into $\eqref{dynInIMinusP}$ and yields:

	$\displaystyle\bm{\lambda}=(\bm{J}_{c}^{T})^{+}\Big{[}$	$\displaystyle(\bm{I}-\bm{P})[\bar{\bm{M}}(\bm{\tau}_{m}-\bm{Ph}+\dot{\bm{P}}\dot{\bm{q}})+\bm{h}]-\bm{\tau}_{c}$		(5)
		$\displaystyle+(\bm{I}-\bm{P})(\bar{\bm{M}}\bm{P}-\bm{I})\bm{J}_{x}^{T}\bm{F}_{x}\Big{]},$

with $\bar{\bm{M}}=\bm{M}\bm{M}_{c}^{-1}$, and is the constraint inertia matrix and is always invertible [10]. More details can be found in [6].

To achieve the desired locomotion behavior, a Cartesian impedance controller is applied to the robot. We start with a closed-loop system equation that reflects a mechanical impedance of the end-effector in multi-dimension to the external disturbance $\bm{F}_{x}$, as:

$$\bm{\Lambda}_{d}\ddot{\bm{e}}+\bm{K}_{d}\dot{\bm{e}}+\bm{K}_{p}\bm{e}=\bm{F}_{x},$$

(6)

where $\bm{e}=\bm{x}-\bm{x}_{d}$ denotes the pose error of the end-effector between the current and desired one, $\bm{\Lambda}_{d}$, $\bm{K}_{d}$, and $\bm{K}_{p}$ are the desired task-space inertia, damping, and stiffness matrices. Note that the torso and swing feet can be treated as the end-effector during the controller design. Again, we apply the equation $(\bm{I}-\bm{P})\ddot{\bm{q}}=\dot{\bm{P}}\dot{\bm{q}}$, as we compute the expression of $\ddot{\bm{q}}$ in Section II-A, to equation (2), then pre-multiply the resultant equation by $\bm{J}_{x}\bm{M}_{c}^{-1}$, replace $\bm{J}_{x}\ddot{\bm{q}}$ with $\ddot{\bm{x}}-\dot{\bm{J}_{x}}\dot{\bm{q}}$, and pre-multiply it by the operational space inertia matrix $\bm{\Lambda}_{c}=(\bm{J}_{x}\bm{M}_{c}^{-1}\bm{P}\bm{J}_{x}^{T})^{-1}$, the equation becomes:

$$\bm{\Lambda}_{c}\ddot{\bm{x}}+\bm{h}_{c}=\bm{\Lambda}_{c}\bm{J}_{x}\bm{M}_{c}^{-1}\bm{PB}\bm{\tau}+\bm{F}_{x},$$

(7)

where $\bm{J}_{x}$ is the task Jacobian defined by $\dot{\bm{x}}=\bm{J}_{x}\dot{\bm{q}}$, and $\bm{h}_{c}=\bm{\Lambda}_{c}\bm{J}_{x}\bm{M}_{c}^{-1}(\bm{Ph}-\dot{\bm{P}}\dot{\bm{q}})-\bm{\Lambda}_{c}\dot{\bm{J}}_{x}\dot{\bm{q}}$ is the operational space non-linear effect. At this moment, we consider the robot is fully actuated, such that:

$$\bm{PB\tau}=\bm{P}\bm{J}_{x}^{T}\bm{F}_{d,x}=\bm{\tau}_{m,d},$$

(8)

with $\bm{F}_{d,x}$ is the designed control law, and $\bm{\tau}_{m,d}$ is the desired motion torques. Then, equation $\eqref{unsimplified modified projection dynamics}$ will be simplified as:

$$\bm{\Lambda}_{c}\ddot{\bm{x}}+\bm{h}_{c}=\bm{F}_{d,x}+\bm{F}_{x}.$$

(9)

The control law $\bm{F}_{d,x}$ can be defined to achieve the impedance response of equation $\eqref{impedance_response}$ as:

$$\bm{F}_{d,x}=\bm{h}_{c}+\bm{\Lambda}_{c}\ddot{\bm{x}}_{d}-\bm{K}_{d}\dot{\bm{e}}-\bm{K}_{p}\bm{e},$$

(10)

where $\bm{\Lambda}_{d}=\bm{\Lambda}_{c}$ is assigned [11].

As described in Section II-A, once $\ddot{\bm{q}}$ is determined, the contact forces can be controlled by applying corresponding constraint torques. This has been shown in [12, 13], where the robot executed a desired motion and applied a desired force at the contact point simultaneously. Here, we extend the constraint torques control law $\bm{\tau}_{c,d}$ with the external disturbance by:

	$\displaystyle\bm{\tau}_{c,d}=(\bm{I}-\bm{P})[$	$\displaystyle\bar{\bm{M}}(\bm{\tau}_{m,d}-\bm{Ph}+\dot{\bm{P}}\dot{\bm{q}}+\bm{P}\bm{J}_{x}^{T}\bm{F}_{x})+\bm{h}]$		(11)
		$\displaystyle-(\bm{I}-\bm{P})\bm{J}_{x}^{T}\bm{F}_{x}-\bm{J}_{c}^{T}\bm{\lambda}_{d}.$

The complete torque commands that can accomplish the desired motion and contact force are the sum of $\bm{\tau}_{m,d}$ and $\bm{\tau}_{c,d}$, i.e.:

$$\bm{B}\bm{\tau}=\bm{\tau}_{d}=\bm{\tau}_{m,d}+\bm{\tau}_{c,d}$$

(12)

In Section II-B and II-C, the systems are assumed to be fully actuated when deriving the desired motion and constraint torque commands. However, $\bm{B\tau}\neq\bm{\tau}_{d}$ in general for underactuated robots. To satisfy the underactuation constraint $\bm{B\tau}=\bm{\tau}_{d}$, [7] extended the PIDC formulation for the underactuated systems to realize the contact-consistent motion and contact force control by adding an additional constraint torques $(\bm{I-P})\bm{\tau}_{u}$ to equation (12) as:

	$\displaystyle\bm{\tau}_{d}$	$\displaystyle=\bm{\tau}_{m,d}+\bm{\tau}_{c,d}+(\bm{I-P})\bm{\tau}_{u}$		(14)
		$\displaystyle=[\bm{PB}]^{+}\bm{\tau}_{m,d}+\left[\bm{I}-[(\bm{I}-\bm{B})(\bm{I}-\bm{P})]^{+}\right](\bm{I}-\bm{P})\bm{\tau}_{c,d}.$

Indeed, the additional constraint torques $(\bm{I-P})\bm{\tau}_{u}$ does not generate any motion to violate the underlying task, but it derails the force control from its desired value, as $\bm{\tau}_{c,d}$ and $(\bm{I-P})\bm{\tau}_{u}$ are both in the constraint space. In other words, this approach fails to achieve exact force control even at any of all contact points.

Exact interaction force control is an essential technique to accomplish many real tasks, e.g., grinding, polishing, and screwing. This control problem has been widely studied for fully actuated robots. Nevertheless, our work focuses on extending the PIDC framework for the underactuated systems, specifically, the floating base quadrupedal robots, towards exerting an accurate desired active force at the contact point.

Compared to the methods used in [6, 4, 7], we aim to directly find an optimal solution that solves the optimization problem $\min\left\|\bm{B\tau}-\bm{\tau}_{d}\right\|^{2}$ without inducing any additional null space components or constraint forces. However, because of the underactuation as mentioned in [5, 14], the desired motion and contact forces cannot be implemented by directly solving this optimization problem. Fortunately, [14] offers an idea of constrained hierarchical optimization that orthogonally decomposes the error norm into two subspaces at first,

$$\left\|\bm{B\tau}-\bm{\tau}_{d}\right\|^{2}=\left\|\bm{B\tau}-\bm{\tau}_{d}\right\|^{2}_{\bm{P}}+\left\|\bm{B\tau}-\bm{\tau}_{d}\right\|^{2}_{\bm{I}-\bm{P}},$$

(15)

where $\left\|\bm{B\tau}-\bm{\tau}_{d}\right\|^{2}_{\bm{P}}$ and $\left\|\bm{B\tau}-\bm{\tau}_{d}\right\|^{2}_{\bm{I}-\bm{P}}$ are the motion space and constraint space error norm, respectively. Then, the optimization problem of (15) will be solved hierarchically in the motion and constraint space. Following the constrained hierarchical optimization formulation, the rest of this section outlines our PIDC-based method for the contact force control of the underactuated robots without using FS at the contact point.

As stated in [10], the constraint forces $\bm{\lambda}$ are affected by the constraint torques $\bm{\tau}_{c}$ directly and motion torques $\bm{\tau}_{m}$ indirectly via a cross-coupling factor $\mu=(\bm{I}-\bm{P})\bar{\bm{M}}$. It is apparent that $\bm{\tau}_{m}$ must be determined first, otherwise determining $\bm{\tau}_{m}$ afterward would affect the contact force control. Similar ideas have also been seen in [12, 13, 7, 15, 8]. Therefore, the motion space error norm is minimized at first.

The selection matrix $\bm{B}$ can be decomposed into two selection matrices $\bm{B}=\bm{B}_{m}+\bm{B}_{f}$, where $\bm{B}_{m}$ selects the actuators to perform the desired underlying motion, and $\bm{B}_{f}$ selects the actuators for the desired contact force control. An intuitive instance would be a quadruped robot using its front right leg to impose a desired contact force on the button while using the other legs to hold the desired torso pose, as shown in Fig. 1. In this case, $\bm{B}_{f}$ activates the front right leg to exert the desired interaction force, whereas $\bm{B}_{m}$ activates the remaining legs for the motion task. This decomposition of the selection matrix can also be generalized to other underactuation systems, e.g., collaborative manipulation of the multi-robot team on a free-floating object [7, 15, 8]. Based on the selection matrices design, (15) can be further decomposed into:

		$\displaystyle\left\|\bm{B\tau}-\bm{\tau}_{d}\right\|^{2}_{\bm{P}}+\left\|\bm{B\tau}-\bm{\tau}_{d}\right\|^{2}_{\bm{I}-\bm{P}}$		(16)
	$\displaystyle=$	$\displaystyle\|\bm{B}_{m}\bm{\tau}+\bm{B}_{f}\bm{\tau}-\bm{\tau}_{d}\|^{2}_{\bm{P}}+\|\bm{B}_{m}\bm{\tau}+\bm{B}_{f}\bm{\tau}-\bm{\tau}_{d}\|^{2}_{\bm{I}-\bm{P}}$

where $\bm{B}_{f}\bm{\tau}$ in motion space and $\bm{B}_{m}\bm{\tau}$ in constraint space are set to zero vector, respectively, because we do not expect that the joints used for the force control contribute to the motion task, and vice versa. Thus, the optimization problem in motion space becomes:

	$$\mathrm{\textbf{QP}}_{\mathrm{\textbf{I}}}:\bm{\tau}_{\mathrm{I}}^{*}=\underset{\bm{\tau}}{\arg\min}\left\|\bm{B}_{m}\bm{\tau}-\bm{\tau}_{d}\right\|_{\bm{P}}^{2},$$		(17a)
	$$\text{s. t. \,\,\,}\bm{\tau}_{\text{min}}\leq\bm{B}^{T}\bm{P}\bm{B}_{m}\bm{\tau}\leq\bm{\tau}_{\text{max}},$$		(17b)

where (17b) is the actuation torque limits of the robot. Here, the pre-multiplication of $\bm{B}^{T}$ is to extract the actuated joints from the vertical concatenation vector of all the underactuated and actuated joints, because it is not a square matrix. Note the other physical constraints including unilateral contact and friction cones are considered within the optimization in the constraint space.

Since the motion space error norm has been minimized, we need to optimize the constraint space error norm subject to the physical constraints.