Consider a mechanical system of articulated rigid bodies, like a robot arm, with
generalized position vector q∈Rm and generalized velocity
vector v∈Rn (note that, in general, m≥n, such as when
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=n and
v=q˙). The dynamics of the system are given by the equation
τ=M(q)v˙+h(q,v),
where τ∈Rn 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 is the (positive definite)
mass matrix, v˙∈Rn is the generalized acceleration vector (the dot denotes a time derivative), and
h∈Rn is the vector of gravitational, centrifugal, and Coriolis forces.
Going forward, we will suppress the dependencies of M and h on
q and 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) that tells us how the system evolves through time in
response to an applied force τ. When the system is subjected to
additional constraints, τ acts in such a way as to enforce those
constraints. Consider a constraint that depends linearly on the generalized
acceleration,
Av˙+b=0,
where A∈Rp×n and b∈Rp, p≤n,
can both depend on q and v. In our setting—that of the control
engineer—we get to choose the constraint (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τc=Mv˙0+h,=Mv˙c+h,
denote the unconstrained and constrained dynamics, respectively. Both have the
form (1); the difference is that τ0 is the
external force acting on the free system (often but not necessarily τ0=0)
while τc is the force necessary to ensure the system’s motion
satisfies (2). From a control design perspective, we get to
chooseτ0, but we must solve for the value(s) of τc
consistent with the constraint. The choice of τ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 of the constrained system is a close as possible
to the unconstrained value v˙0, in the sense of M-weighted
least squares. Mathematically, this means that the constrained dynamics can be
described by the optimization problem
is the mass-matrix-weighted least squares objective. Notice that
v˙0 is really just a dummy optimization variable, since M
is always strictly positive definite and can therefore be inverted to uniquely
determine v˙0 from τ0. Indeed, by inverting M to
rearrange the dynamics equations and substituting them into the objective, we
obtain the equivalent problem
where we have eliminated the variable 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 that is closest to some nominal
value τ0 (again in the sense of least squares but this time weighted
by M−1) while satisfying the dynamics
(4) and the constraint (2).
Let’s solve this optimization problem. Rearranging the dynamics to
v˙c=M−1(τc−h) and substituting it into the
constraint Av˙c+b=0, we can further simplify the problem to
where ΛA=(AM−1AT)+ is called the
effective mass matrix (in constraint space) and 1n denotes the
n×n identity matrix. We can interpret τc as those forces
that must be applied to the system (1) to make its motion
satisfy the constraint (2) (notice that when p<n, there
are infinite possible values of τc, and the one chosen depends on
the value of τ0).
Constrained Dynamics
Substituting τc back into the constrained dynamics
(4) to solve for v˙c while also making
use of (3), we get the equations of motion for
the constrained system
v˙c=AM+b+PM(A)v˙0,
where
AM+:=M−1ATΛA
is the M-weighted pseudoinverse of A and
PM(A):=1n−AM+A
is the M-weighted projector onto the nullspace of 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-weighted least squares problem
xminsubject to(1/2)∥x−y∥M2Ax=b
is
x⋆=AM+b+PM(A)y
which we can interpret as the vector x satisfying the constraint
Ax=b that is closest to y in the sense of the
Mahalanobis distance
dM(x,y):=(x−y)TM(x−y).
The M-weighted pseudoinverse AM+ is a generalized
inverse of A, which means that it satisfies
AAM+A=A (which we prove in the appendix
below) and therefore PM(A) satisfies
APM(A)=A−AAM+A=A−A=0,
so v˙0 in (8) never interferes with
the constraint (2).
For a robot manipulator (or other mechanical system), we can interpret
AM+ as the solution to the problem of finding the generalized
velocity v satisfying a constraint Av=b that minimizes
the instantaneous kinetic energy of the system; AM+ is
therefore sometimes called the dynamically-consistent pseudoinverse.
When M=1n, dM(x,y) becomes the usual Euclidean
distance, AM+ becomes the Moore-Penrose pseudoinverse, and
PM(A) becomes the orthogonal projector onto the nullspace
of 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. To track an end effector trajectory, we want to design an
input torque u∈Rn that ensures the dynamics satisfy the
constraint
Δx(q)=0,
where Δx(q)∈R6 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); in particular, we will differentiate
(11) by time until v˙ appears. The first
time derivative is ξd−ξ(q,v)=0, where
ξd∈R6 is the desired spatial velocity and
ξ(q,v)∈R6 is the actual spatial velocity.
Substituting in the differential kinematics relationship
ξ(q,v)=J(q)v, where
J(q)∈R6×n is the manipulator Jacobian matrix,
we get the velocity-level constraint
ξd−J(q)v=0.
This expression still does not contain v˙, so we differentiate once
more to get the acceleration-level constraint
ξ˙d−J(q)v˙−J˙(q,v)v=0,
which finally contains (a linear dependence on) v˙ and so we can
write it in the same form as (2) by taking
A=−J(q) and
b=ξ˙d−J˙(q,v)v. Substituting these
values into the expression for the torque input (7) and
suppressing dependencies, we get
τc=JTΛJ(JM−1h+ξ˙d−J˙v)+PM(J)τ0.
The torque (13) is what the controller computes to send to the
robot to achieve motion satisfying the constraint (12);
to make this more clear, let us replace the symbol τc with the
standard control input symbol u, yielding the control law
u=JTΛJ(JM−1h+ξ˙d−J˙v)+PM(J)τ0.
The Effective Mass Matrix
The effective mass matrix Λ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 to obtain the (minimum) kinetic energy of the overall
system. To see this, consider the optimization problem
vminsubject to(1/2)∥v∥M2Jv=ξ,
which seeks the generalized velocity v that achieves a given end effector
velocity ξ while minimizing the system’s kinetic energy. As we can see
from (10), the solution is
v⋆=JM+ξ, which corresponds to kinetic energy
which shows that we obtain an equivalent expression for the minimum kinetic
energy using either M and v or ΛJ and ξ.
This kinetic energy perspective of Λ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 loses rank, implying that the eigenvalue of
Λ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). A basic problem with it is that
enforcing the constraint (11) on the acceleration level
(12) can lead to drift and error accumulation over time.
To combat this, we replace ξ˙d with the reference quantity
ξ˙r=KpΔx+KdΔξ+ξ˙d,
which includes proportional-derivative feedback terms on the pose and spatial
velocity errors while keeping ξ˙d as a feedforward term, with
positive definite gain matrices Kp and Kd. Making this
modification to (14), we get
u=JTΛJ(JM−1h+ξ˙r−J˙v)+PM(J)τ0,
which is a standard OSC law.
Configuration Space Compensation
The standard OSC (15) compensates for the gravitational,
centrifugal, and Coriolis vector 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 according to
(3), which gives us the control law
u=JTΛJ(ξ˙r−J˙v)+h+PM(J)Mv˙0,
where h is now completely compensated and we are still free to pursue
secondary objectives by selecting v˙0. The simplest approach is to
chooose v˙0=0, which leaves us with
u=JTΛJ(ξ˙r−J˙v)+h,
but we could also design 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+, as defined in (9), is a
generalized inverse of A, meaning it satisfies
AAM+A=A. We have
AAM+A=AM−1AT(AM−1AT)+A=ΛA+ΛAA.
Let
P(ΛA)=1n−ΛA+ΛA
be the (Euclidean) projector onto the nullspace of ΛA. We have
P(ΛA)ΛA+=(1n−ΛA+ΛA)ΛA+=ΛA+−ΛA+=0,
and therefore
P(ΛA)ΛA+PT(ΛA)=0
as well. Taking the trace, we get
where each ri is a row of P(ΛA)A. Since M is strictly
positive definite, riTM−1ri≥0 with equality if and
only if ri=0. This implies that each ri=0, and
therefore P(ΛA)A=0, so
0=P(ΛA)A=(1n−ΛA+ΛA)A=A−AAM+M−1ATΛAA,
which we rearrange to give us the desired result AAM+A=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.
Notice the similarity to the OSC (16). Indeed, if we ignore the
inertial parameters of the robot by simply setting M=1n, then
(17) and (16) become equal. However, in general
resolved acceleration control does not satisfy the Principle of Least
Constraint and is therefore not dynamically consistent.