Operational Space Control and the Principle of Least Constraint

In this post we are going to derive the venerable operational space controller (OSC), first introduced in (Khatib, 1987), using Gauss’s Principle of Least Constraint. For further information on OSC, I would recommend looking at (Nakanishi et al., 2008) and (Peters et al., 2008) (though both are paywalled, unfortunately).

Rigid Body Systems

Consider a mechanical system of articulated rigid bodies, like a robot arm, with generalized position vector q∈Rm\bm{q}\in\mathbb{R}^m and generalized velocity vector v∈Rn\bm{v}\in\mathbb{R}^n (note that, in general, m≥nm\geq n, such as when q\bm{q} includes a 3D rotation represented in non-minimal coordinates, like a quaternion, but for the puposes of this post you are safe to imagine m=nm=n and v=q˙\bm{v}=\dot{\bm{q}}). The dynamics of the system are given by the equation

τ=M(q)v˙+h(q,v),\begin{equation}\label{1} \bm{\tau} = \bm{M}(\bm{q})\dot{\bm{v}} + \bm{h}(\bm{q},\bm{v}), \end{equation}

where τ∈Rn\bm{\tau}\in\mathbb{R}^n is the applied generalized force vector (which includes both applied motor torques as well as environmental forces like those exerted by contacts), M∈Rn×n\bm{M}\in\mathbb{R}^{n\times n} is the (positive definite) mass matrix, v˙∈Rn\dot{\bm{v}}\in\mathbb{R}^n is the generalized acceleration vector (the dot denotes a time derivative), and h∈Rn\bm{h}\in\mathbb{R}^n is the vector of gravitational, centrifugal, and Coriolis forces. Going forward, we will suppress the dependencies of M\bm{M} and h\bm{h} on q\bm{q} and v\bm{v} to simplify notation. Note that we explored the derivation of these dynamics using Lagrangian mechanics in this previous post.

The Principle of Least Constraint

The dynamics are described by the second-order ordinary differential equation (1)\eqref{1} that tells us how the system evolves through time in response to an applied force τ\bm{\tau}. When the system is subjected to additional constraints, τ\bm{\tau} acts in such a way as to enforce those constraints. Consider a constraint that depends linearly on the generalized acceleration,

Av˙+b=0,\begin{equation}\label{2} \bm{A}\dot{\bm{v}} + \bm{b} = \bm{0}, \end{equation}

where A∈Rp×n\bm{A}\in\mathbb{R}^{p\times n} and b∈Rp\bm{b}\in\mathbb{R}^p, p≤np\leq n, can both depend on q\bm{q} and v\bm{v}. In our setting—that of the control engineer—we get to choose the constraint (2)\eqref{2}. As we will see below, OSC imposes a full end effector tracking constraint on the dynamics, but we could also impose more limited constraints like requiring an upright orientation to avoid spilling a grasped container.

Let

τ0=Mv˙0+h,τc=Mv˙c+h,\begin{align} \bm{\tau}_0 &= \bm{M}\dot{\bm{v}}_0 + \bm{h},\label{3} \\ \bm{\tau}_c &= \bm{M}\dot{\bm{v}}_c + \bm{h},\label{4} \end{align}

denote the unconstrained and constrained dynamics, respectively. Both have the form (1)\eqref{1}; the difference is that τ0\bm{\tau}_0 is the external force acting on the free system (often but not necessarily τ0=0\bm{\tau}_0=\bm{0}) while τc\bm{\tau}_c is the force necessary to ensure the system’s motion satisfies (2)\eqref{2}. From a control design perspective, we get to choose τ0\bm{\tau}_0, but we must solve for the value(s) of τc\bm{\tau}_c consistent with the constraint. The choice of τ0\bm{\tau}_0 dictates how the system behaves with its remaining unconstrained degrees of freedom; we will discuss some of these choices later in the post.

Gauss’s Principle of Least Constraint tells us that the constraints interfere with the unconstrained dynamics in a minimal way. That is, the generalized acceleration v˙c\dot{\bm{v}}_c of the constrained system is a close as possible to the unconstrained value v˙0\dot{\bm{v}}_0, in the sense of M\bm{M}-weighted least squares. Mathematically, this means that the constrained dynamics can be described by the optimization problem

min⁡τc,v˙0,v˙c(1/2)∥v˙c−v˙0∥M2subject toτ0=Mv˙0+hτc=Mv˙c+h0=Av˙c+b,\begin{equation*} \begin{aligned} \min_{\bm{\tau}_c,\dot{\bm{v}}_0,\dot{\bm{v}}_c} &\quad (1/2)\|\dot{\bm{v}}_c-\dot{\bm{v}}_0\|_{\bm{M}}^2 \\ \text{subject to} &\quad \bm{\tau}_0 = \bm{M}\dot{\bm{v}}_0 + \bm{h} \\ &\quad \bm{\tau}_c = \bm{M}\dot{\bm{v}}_c + \bm{h} \\ &\quad \bm{0} = \bm{A}\dot{\bm{v}}_c + \bm{b}, \end{aligned} \end{equation*}

where

∥v˙c−v˙0∥M2:=(v˙c−v˙0)TM(v˙c−v˙0)\begin{equation*} \|\dot{\bm{v}}_c-\dot{\bm{v}}_0\|_{\bm{M}}^2 := (\dot{\bm{v}}_c-\dot{\bm{v}}_0)^T\bm{M}(\dot{\bm{v}}_c-\dot{\bm{v}}_0) \end{equation*}

is the mass-matrix-weighted least squares objective. Notice that v˙0\dot{\bm{v}}_0 is really just a dummy optimization variable, since M\bm{M} is always strictly positive definite and can therefore be inverted to uniquely determine v˙0\dot{\bm{v}}_0 from τ0\bm{\tau}_0. Indeed, by inverting M\bm{M} to rearrange the dynamics equations and substituting them into the objective, we obtain the equivalent problem

min⁡τc,v˙c(1/2)∥τc−τ0∥M−12subject toτc=Mv˙c+h0=Av˙c+b,\begin{equation*} \begin{aligned} \min_{\bm{\tau}_c,\dot{\bm{v}}_c} &\quad (1/2)\|\bm{\tau}_c-\bm{\tau}_0\|_{\bm{M}^{-1}}^2 \\ \text{subject to} &\quad \bm{\tau}_c = \bm{M}\dot{\bm{v}}_c + \bm{h} \\ &\quad \bm{0} = \bm{A}\dot{\bm{v}}_c + \bm{b}, \end{aligned} \end{equation*}

where we have eliminated the variable v˙0\dot{\bm{v}}_0 and expressed the objective in terms of the generalized forces. This form makes it clear that we are seeking a generalized force τc\bm{\tau}_c that is closest to some nominal value τ0\bm{\tau}_0 (again in the sense of least squares but this time weighted by M−1\bm{M}^{-1}) while satisfying the dynamics (4)\eqref{4} and the constraint (2)\eqref{2}.

Let’s solve this optimization problem. Rearranging the dynamics to v˙c=M−1(τc−h)\dot{\bm{v}}_c=\bm{M}^{-1}(\bm{\tau}_c-\bm{h}) and substituting it into the constraint Av˙c+b=0\bm{A}\dot{\bm{v}}_c+\bm{b}=\bm{0}, we can further simplify the problem to

min⁡τc(1/2)∥τc−τ0∥M−12subject to0=AM−1(τc−h)+b,\begin{equation}\label{5} \begin{aligned} \min_{\bm{\tau}_c} &\quad (1/2)\|\bm{\tau}_c-\bm{\tau}_0\|_{\bm{M}^{-1}}^2 \\ \text{subject to} &\quad \bm{0} = \bm{A}\bm{M}^{-1}(\bm{\tau}_c-\bm{h}) + \bm{b}, \end{aligned} \end{equation}

where we are now down to a single variable τc\bm{\tau}_c. The Lagrangian of this problem is

L(τc,λ)=(1/2)∥τc−τ0∥M−12+λT(AM−1(τc−h)+b)\begin{equation*} \mathcal{L}(\bm{\tau}_c,\bm{\lambda}) = (1/2)\|\bm{\tau}_c-\bm{\tau}_0\|_{\bm{M}^{-1}}^2 + \bm{\lambda}^T(\bm{A}\bm{M}^{-1}(\bm{\tau}_c-\bm{h}) + \bm{b}) \end{equation*}

with Lagrange multipliers λ∈Rp\bm{\lambda}\in\mathbb{R}^p. Taking the derivative of L\mathcal{L} with respect to τc\bm{\tau}_c and setting it equal to zero, we get

0=M−1(τc−τ0)+M−1ATλ.\begin{equation*} \bm{0} = \bm{M}^{-1}(\bm{\tau}_c-\bm{\tau}_0) + \bm{M}^{-1}\bm{A}^T\bm{\lambda}. \end{equation*}

Multiplying through by M\bm{M} and rearranging yields

τc=τ0−ATλ,\begin{equation}\label{6} \bm{\tau}_c = \bm{\tau}_0-\bm{A}^T\bm{\lambda}, \end{equation}

which we substitute back into the constraint of (5)\eqref{5} to obtain

0=AM−1(τ0−ATλ−h)+b.\begin{equation*} \bm{0}=\bm{A}\bm{M}^{-1}(\bm{\tau}_0-\bm{A}^T\bm{\lambda}-\bm{h}) + \bm{b}. \end{equation*}

Solving for λ\bm{\lambda}, we get

λ=(AM−1AT)+(AM−1(τ0−h)+b),\begin{equation*} \bm{\lambda}=(\bm{A}\bm{M}^{-1}\bm{A}^T)^+(\bm{A}\bm{M}^{-1}(\bm{\tau}_0-\bm{h}) + \bm{b}), \end{equation*}

where AM−1AT\bm{A}\bm{M}^{-1}\bm{A}^T is known as the Delassus matrix and (⋅)+(\cdot)^+ denotes the Moore-Penrose pseudoinverse. Substituting λ\bm{\lambda} back into (6)\eqref{6}, we obtain

τc=τ0−AT(AM−1AT)+(AM−1(τ0−h)+b)=ATΛA(AM−1h−b)+(1n−ATΛAAM−1)τ0.\begin{equation}\label{7} \begin{aligned} \bm{\tau}_c &= \bm{\tau}_0-\bm{A}^T(\bm{A}\bm{M}^{-1}\bm{A}^T)^+(\bm{A}\bm{M}^{-1}(\bm{\tau}_0-\bm{h}) + \bm{b}) \\ &= \bm{A}^T\bm{\Lambda}_{\bm{A}}(\bm{A}\bm{M}^{-1}\bm{h} - \bm{b}) + (\bm{1}_n-\bm{A}^T\bm{\Lambda}_{\bm{A}}\bm{A}\bm{M}^{-1})\bm{\tau}_0. \end{aligned} \end{equation}

where ΛA=(AM−1AT)+\bm{\Lambda}_{\bm{A}}=(\bm{A}\bm{M}^{-1}\bm{A}^T)^+ is called the effective mass matrix (in constraint space) and 1n\bm{1}_n denotes the n×nn\times n identity matrix. We can interpret τc\bm{\tau}_c as those forces that must be applied to the system (1)\eqref{1} to make its motion satisfy the constraint (2)\eqref{2} (notice that when p<np<n, there are infinite possible values of τc\bm{\tau}_c, and the one chosen depends on the value of τ0\bm{\tau}_0).

Constrained Dynamics

Substituting τc\bm{\tau}_c back into the constrained dynamics (4)\eqref{4} to solve for v˙c\dot{\bm{v}}_c while also making use of (3)\eqref{3}, we get the equations of motion for the constrained system

v˙c=AM+b+PM(A)v˙0,\begin{equation}\label{8} \dot{\bm{v}}_c = \bm{A}^+_{\bm{M}}\bm{b} + \bm{P}_{\bm{M}}(\bm{A})\dot{\bm{v}}_0, \end{equation}

where

AM+:=M−1ATΛA\begin{equation}\label{9} \bm{A}^+_{\bm{M}} := \bm{M}^{-1}\bm{A}^T\bm{\Lambda}_{\bm{A}} \end{equation}

is the M\bm{M}-weighted pseudoinverse of A\bm{A} and

PM(A):=1n−AM+A\begin{equation*} \bm{P}_{\bm{M}}(\bm{A}) := \bm{1}_n-\bm{A}^+_{\bm{M}}\bm{A} \end{equation*}

is the M\bm{M}-weighted projector onto the nullspace of A\bm{A}. These two expressions generalize the usual (Euclidean) pseudoinverse and nullspace projector that we explored in this previous post. That is, the solution to the M\bm{M}-weighted least squares problem

min⁡x(1/2)∥x−y∥M2subject toAx=b\begin{equation*} \begin{aligned} \min_{\bm{x}} &\quad (1/2)\|\bm{x}-\bm{y}\|^2_{\bm{M}} \\ \text{subject to} &\quad \bm{A}\bm{x} = \bm{b} \end{aligned} \end{equation*}

is

x⋆=AM+b+PM(A)y\begin{equation}\label{10} \bm{x}^{\star} = \bm{A}^+_{\bm{M}}\bm{b} + \bm{P}_{\bm{M}}(\bm{A})\bm{y} \end{equation}

which we can interpret as the vector x\bm{x} satisfying the constraint Ax=b\bm{A}\bm{x}=\bm{b} that is closest to y\bm{y} in the sense of the Mahalanobis distance

dM(x,y):=(x−y)TM(x−y).\begin{equation*} d_{\bm{M}}(\bm{x},\bm{y}) := \sqrt{(\bm{x}-\bm{y})^T\bm{M}(\bm{x}-\bm{y})}. \end{equation*}

The M\bm{M}-weighted pseudoinverse AM+\bm{A}^+_{\bm{M}} is a generalized inverse of A\bm{A}, which means that it satisfies AAM+A=A\bm{A}\bm{A}^+_{\bm{M}}\bm{A}=\bm{A} (which we prove in the appendix below) and therefore PM(A)\bm{P}_{\bm{M}}(\bm{A}) satisfies

APM(A)=A−AAM+A=A−A=0,\begin{equation*} \bm{A}\bm{P}_{\bm{M}}(\bm{A}) = \bm{A}-\bm{A}\bm{A}^+_{\bm{M}}\bm{A} = \bm{A} - \bm{A} = \bm{0}, \end{equation*}

so v˙0\dot{\bm{v}}_0 in (8)\eqref{8} never interferes with the constraint (2)\eqref{2}.

For a robot manipulator (or other mechanical system), we can interpret AM+\bm{A}^+_{\bm{M}} as the solution to the problem of finding the generalized velocity v\bm{v} satisfying a constraint Av=b\bm{A}\bm{v}=\bm{b} that minimizes the instantaneous kinetic energy of the system; AM+\bm{A}^+_{\bm{M}} is therefore sometimes called the dynamically-consistent pseudoinverse. When M=1n\bm{M}=\bm{1}_n, dM(x,y)d_{\bm{M}}(\bm{x},\bm{y}) becomes the usual Euclidean distance, AM+\bm{A}^+_{\bm{M}} becomes the Moore-Penrose pseudoinverse, and PM(A)\bm{P}_{\bm{M}}(\bm{A}) becomes the orthogonal projector onto the nullspace of A\bm{A}.

Operational Space Control

We are now (finally) ready to derive the operational space controller, which is simply a particular case of the constrained dynamics we just discussed. The goal of operational space control is to generate input torques that track a desired end effector trajectory (that is, a trajectory in “task” or “operational” space). The pose of the end effector is related to the generalized positions through the forward kinematics of the robot, which is a nonlinear function of q\bm{q}. To track an end effector trajectory, we want to design an input torque u∈Rn\bm{u}\in\mathbb{R}^n that ensures the dynamics satisfy the constraint

Δx(q)=0,\begin{equation}\label{11} \Delta\bm{x}(\bm{q}) = \bm{0}, \end{equation}

where Δx(q)∈R6\Delta\bm{x}(\bm{q})\in\mathbb{R}^6 is the error between the desired and actual poses. However, we need to manipulate this constraint to get it into the same form as (2)\eqref{2}; in particular, we will differentiate (11)\eqref{11} by time until v˙\dot{\bm{v}} appears. The first time derivative is ξd−ξ(q,v)=0\bm{\xi}_d-\bm{\xi}(\bm{q},\bm{v})=\bm{0}, where ξd∈R6\bm{\xi}_d\in\mathbb{R}^6 is the desired spatial velocity and ξ(q,v)∈R6\bm{\xi}(\bm{q},\bm{v})\in\mathbb{R}^6 is the actual spatial velocity. Substituting in the differential kinematics relationship ξ(q,v)=J(q)v\bm{\xi}(\bm{q},\bm{v})=\bm{J}(\bm{q})\bm{v}, where J(q)∈R6×n\bm{J}(\bm{q})\in\mathbb{R}^{6\times n} is the manipulator Jacobian matrix, we get the velocity-level constraint

ξd−J(q)v=0.\begin{equation*} \bm{\xi}_d - \bm{J}(\bm{q})\bm{v} = \bm{0}. \end{equation*}

This expression still does not contain v˙\dot{\bm{v}}, so we differentiate once more to get the acceleration-level constraint

ξ˙d−J(q)v˙−J˙(q,v)v=0,\begin{equation}\label{12} \dot{\bm{\xi}}_d - \bm{J}(\bm{q})\dot{\bm{v}} - \Jdot(\bm{q},\bm{v})\bm{v} = \bm{0}, \end{equation}

which finally contains (a linear dependence on) v˙\dot{\bm{v}} and so we can write it in the same form as (2)\eqref{2} by taking A=−J(q)\bm{A}=-\bm{J}(\bm{q}) and b=ξ˙d−J˙(q,v)v\bm{b}=\dot{\bm{\xi}}_d-\Jdot(\bm{q},\bm{v})\bm{v}. Substituting these values into the expression for the torque input (7)\eqref{7} and suppressing dependencies, we get

τc=JTΛJ(JM−1h+ξ˙d−J˙v)+PM(J)τ0.\begin{equation}\label{13} \bm{\tau}_c = \bm{J}^T\bm{\Lambda}_{\bm{J}}(\bm{J}\bm{M}^{-1}\bm{h} + \dot{\bm{\xi}}_d-\Jdot\bm{v}) + \bm{P}_{\bm{M}}(\bm{J})\bm{\tau}_0. \end{equation}

The torque (13)\eqref{13} is what the controller computes to send to the robot to achieve motion satisfying the constraint (12)\eqref{12}; to make this more clear, let us replace the symbol τc\bm{\tau}_c with the standard control input symbol u\bm{u}, yielding the control law

u=JTΛJ(JM−1h+ξ˙d−J˙v)+PM(J)τ0.\begin{equation}\label{14} \bm{u} = \bm{J}^T\bm{\Lambda}_{\bm{J}}(\bm{J}\bm{M}^{-1}\bm{h} + \dot{\bm{\xi}}_d-\Jdot\bm{v}) + \bm{P}_{\bm{M}}(\bm{J})\bm{\tau}_0. \end{equation}

The Effective Mass Matrix

The effective mass matrix ΛJ\bm{\Lambda}_{\bm{J}} (usually denoted without a subscript) can be interpreted as the mass matrix corresponding to the motion of the end effector in task space, in the sense that it can be used as an alternative to M\bm{M} to obtain the (minimum) kinetic energy of the overall system. To see this, consider the optimization problem

min⁡v(1/2)∥v∥M2subject toJv=ξ,\begin{equation*} \begin{aligned} \min_{\bm{v}} &\quad (1/2)\|\bm{v}\|_{\bm{M}}^2 \\ \text{subject to} &\quad \bm{J}\bm{v} = \bm{\xi}, \end{aligned} \end{equation*}

which seeks the generalized velocity v\bm{v} that achieves a given end effector velocity ξ\bm{\xi} while minimizing the system’s kinetic energy. As we can see from (10)\eqref{10}, the solution is v⋆=JM+ξ\bm{v}^{\star}=\bm{J}_{\bm{M}}^+\bm{\xi}, which corresponds to kinetic energy

(1/2)(v⋆)TMv⋆=(1/2)ξT(JM+)TMJM+ξ=(1/2)ξT(M−1JTΛJ)TJTΛJξ=(1/2)ξTΛJΛJ+ΛJξ=(1/2)ξTΛJξ,\begin{equation*} \begin{aligned} (1/2)(\bm{v}^{\star})^T\bm{M}\bm{v}^{\star} &= (1/2)\bm{\xi}^T(\bm{J}_{\bm{M}}^+)^T\bm{M}\bm{J}_{\bm{M}}^+\bm{\xi} \\ &= (1/2)\bm{\xi}^T(\bm{M}^{-1}\bm{J}^T\bm{\Lambda}_{\bm{J}})^T\bm{J}^T\bm{\Lambda}_{\bm{J}}\bm{\xi} \\ &= (1/2)\bm{\xi}^T\bm{\Lambda}_{\bm{J}}\bm{\Lambda}_{\bm{J}}^+\bm{\Lambda}_{\bm{J}}\bm{\xi} \\ &= (1/2)\bm{\xi}^T\bm{\Lambda}_{\bm{J}}\bm{\xi}, \end{aligned} \end{equation*}

which shows that we obtain an equivalent expression for the minimum kinetic energy using either M\bm{M} and v\bm{v} or ΛJ\bm{\Lambda}_{\bm{J}} and ξ\bm{\xi}.

This kinetic energy perspective of ΛJ\bm{\Lambda}_{\bm{J}} is particularly interesting when the robot approaches a singularity. A singularity occurs when the robot is in a configuration that instantaneously loses a degree of freedom (i.e., the end effector cannot move in a particular direction). Mathematically, a singularity occurs when J\bm{J} loses rank, implying that the eigenvalue of ΛJ\bm{\Lambda}_{\bm{J}} in the corresponding singular direction blows up to infinity, and therefore any motion along that direction requires infinite energy (and therefore cannot be done).

Feedback

Let’s return to the OSC law (14)\eqref{14}. A basic problem with it is that enforcing the constraint (11)\eqref{11} on the acceleration level (12)\eqref{12} can lead to drift and error accumulation over time. To combat this, we replace ξ˙d\dot{\bm{\xi}}_d with the reference quantity

ξ˙r=KpΔx+KdΔξ+ξ˙d,\begin{equation*} \dot{\bm{\xi}}_r = \bm{K}_p\Delta\bm{x} + \bm{K}_d\Delta\bm{\xi} + \dot{\bm{\xi}}_d, \end{equation*}

which includes proportional-derivative feedback terms on the pose and spatial velocity errors while keeping ξ˙d\dot{\bm{\xi}}_d as a feedforward term, with positive definite gain matrices Kp\bm{K}_p and Kd\bm{K}_d. Making this modification to (14)\eqref{14}, we get

u=JTΛJ(JM−1h+ξ˙r−J˙v)+PM(J)τ0,\begin{equation}\label{15} \bm{u} = \bm{J}^T\bm{\Lambda}_{\bm{J}}(\bm{J}\bm{M}^{-1}\bm{h} + \dot{\bm{\xi}}_r-\Jdot\bm{v}) + \bm{P}_{\bm{M}}(\bm{J})\bm{\tau}_0, \end{equation}

which is a standard OSC law.

Configuration Space Compensation

The standard OSC (15)\eqref{15} compensates for the gravitational, centrifugal, and Coriolis vector h\bm{h} in the operational space, but it is typically desirable to compensate for them directly in configuration space; that is, we want to compensate for them in both the constraint space and the null space. To do so, a sensible choice is to choose τ0\bm{\tau}_0 according to (3)\eqref{3}, which gives us the control law

u=JTΛJ(ξ˙r−J˙v)+h+PM(J)Mv˙0,\begin{equation*} \bm{u} = \bm{J}^T\bm{\Lambda}_{\bm{J}}(\dot{\bm{\xi}}_r-\Jdot\bm{v}) + \bm{h} + \bm{P}_{\bm{M}}(\bm{J})\bm{M}\dot{\bm{v}}_0, \end{equation*}

where h\bm{h} is now completely compensated and we are still free to pursue secondary objectives by selecting v˙0\dot{\bm{v}}_0. The simplest approach is to chooose v˙0=0\dot{\bm{v}}_0=\bm{0}, which leaves us with

u=JTΛJ(ξ˙r−J˙v)+h,\begin{equation}\label{16} \bm{u} = \bm{J}^T\bm{\Lambda}_{\bm{J}}(\dot{\bm{\xi}}_r-\Jdot\bm{v}) + \bm{h}, \end{equation}

but we could also design v˙0\dot{\bm{v}}_0 to steer toward a nominal configuration, avoid singular configurations, or even maximize distance from nearby obstacles, for example.

Conclusion

We have now completed our derivation of the OSC, one of the most popular task-space controllers for torque-controlled robots, from first principles. In particular, we showed that the OSC law is the solution to a particular optimization problem arising from Gauss’s Principle of Least Constraint in classical mechanics, and is therefore consistent with the system dynamics.

Thanks to Karime Pereida for reading a draft of this post.

Appendix

Here we provide some extra details that may be of interest but are not essential for understanding the main body of the post.

Generalized Inverse

Let’s prove that AM+\bm{A}^+_{\bm{M}}, as defined in (9)\eqref{9}, is a generalized inverse of A\bm{A}, meaning it satisfies AAM+A=A\bm{A}\bm{A}^+_{\bm{M}}\bm{A}=\bm{A}. We have

AAM+A=AM−1AT(AM−1AT)+A=ΛA+ΛAA.\begin{equation*} \begin{aligned} \bm{A}\bm{A}^+_{\bm{M}}\bm{A} &= \bm{A}\bm{M}^{-1}\bm{A}^T(\bm{A}\bm{M}^{-1}\bm{A}^T)^+\bm{A} = \bm{\Lambda}_{\bm{A}}^+\bm{\Lambda}_{\bm{A}}\bm{A}. \end{aligned} \end{equation*}

Let P(ΛA)=1n−ΛA+ΛA\bm{P}(\bm{\Lambda}_{\bm{A}})=\bm{1}_n-\bm{\Lambda}_{\bm{A}}^+\bm{\Lambda}_{\bm{A}} be the (Euclidean) projector onto the nullspace of ΛA.\bm{\Lambda}_{\bm{A}}. We have

P(ΛA)ΛA+=(1n−ΛA+ΛA)ΛA+=ΛA+−ΛA+=0,\begin{equation*} \begin{aligned} \bm{P}(\bm{\Lambda}_{\bm{A}})\bm{\Lambda}_{\bm{A}}^+ &= (\bm{1}_n-\bm{\Lambda}_{\bm{A}}^+\bm{\Lambda}_{\bm{A}})\bm{\Lambda}_{\bm{A}}^+ \\ &= \bm{\Lambda}_{\bm{A}}^+-\bm{\Lambda}_{\bm{A}}^+ \\ &= \bm{0}, \end{aligned} \end{equation*}

and therefore P(ΛA)ΛA+PT(ΛA)=0\bm{P}(\bm{\Lambda}_{\bm{A}})\bm{\Lambda}^+_{\bm{A}}\bm{P}^T(\bm{\Lambda}_{\bm{A}})=\bm{0} as well. Taking the trace, we get

0=tr(P(ΛA)ΛA+PT(ΛA))=tr(P(ΛA)AM−1ATPT(ΛA))=∑i=1nriTM−1ri\begin{equation*} \begin{aligned} 0 &= \mathrm{tr}(\bm{P}(\bm{\Lambda}_{\bm{A}})\bm{\Lambda}^+_{\bm{A}}\bm{P}^T(\bm{\Lambda}_{\bm{A}})) \\ &= \mathrm{tr}(\bm{P}(\bm{\Lambda}_{\bm{A}})\bm{A}\bm{M}^{-1}\bm{A}^T\bm{P}^T(\bm{\Lambda}_{\bm{A}})) \\ &= \sum_{i=1}^n \bm{r}_i^T\bm{M}^{-1}\bm{r}_i \end{aligned} \end{equation*}

where each ri\bm{r}_i is a row of P(ΛA)A\bm{P}(\bm{\Lambda}_{\bm{A}})\bm{A}. Since M\bm{M} is strictly positive definite, riTM−1ri≥0\bm{r}_i^T\bm{M}^{-1}\bm{r}_i\geq0 with equality if and only if ri=0\bm{r}_i=\bm{0}. This implies that each ri=0\bm{r}_i=\bm{0}, and therefore P(ΛA)A=0\bm{P}(\bm{\Lambda}_{\bm{A}})\bm{A}=\bm{0}, so

0=P(ΛA)A=(1n−ΛA+ΛA)A=A−AM−1ATΛA⏟AM+A,\begin{equation*} \begin{aligned} \bm{0} &= \bm{P}(\bm{\Lambda}_{\bm{A}})\bm{A} \\ &= (\bm{1}_n-\bm{\Lambda}_{\bm{A}}^+\bm{\Lambda}_{\bm{A}})\bm{A} \\ &= \bm{A} - \bm{A}\underbrace{\bm{M}^{-1}\bm{A}^T\bm{\Lambda}_{\bm{A}}}_{\bm{A}^+_{\bm{M}}}\bm{A}, \end{aligned} \end{equation*}

which we rearrange to give us the desired result AAM+A=A\bm{A}\bm{A}^+_{\bm{M}}\bm{A}=\bm{A}.

Resolved Acceleration Control

Another popular task space control method is called resolved acceleration control (Hsu, Mauser, and Sastry, 1989), which (ignoring nullspace terms) has control law

u=MJ+(ξ˙r−J˙v)+h.\begin{equation}\label{17} \bm{u} = \bm{M}\bm{J}^+(\dot{\bm{\xi}}_r-\Jdot\bm{v}) + \bm{h}. \end{equation}

Notice the similarity to the OSC (16)\eqref{16}. Indeed, if we ignore the inertial parameters of the robot by simply setting M=1n\bm{M}=\bm{1}_n, then (17)\eqref{17} and (16)\eqref{16} become equal. However, in general resolved acceleration control does not satisfy the Principle of Least Constraint and is therefore not dynamically consistent.