Jacobians and velocity kinematics

Forward kinematics tells you where the tool is. This lesson is about how fast the tool moves, and in which direction, when the joints turn. The map from joint rates to tool velocity is the Jacobian. It sits underneath velocity control, singularity checks, force mapping and most numerical inverse kinematics. The lesson is for readers who are comfortable with rotation matrices and want to derive a Jacobian, test it, and use it to control an arm. All code uses numpy and runs in a few seconds.

The Jacobian as a linear map

Write the tool position as a function of the joint angles, x = f(q). Differentiate with respect to time and the chain rule gives ẋ = J(q) q̇, where entry (i, j) of J is ∂fᵢ/∂qⱼ. At a fixed configuration the map is linear in the joint rates. The matrix itself changes as the arm moves, which is where all the interesting behaviour comes from.

There is a direct physical reading. Column j of J is the tool velocity you get when joint j turns at 1 rad/s and every other joint holds still. Keep that picture in mind, because the spatial Jacobian later in the lesson is built from it.

A planar two-link arm by hand

Take an arm in the x–y plane of its base frame. We use REP 103 axes throughout: x forward, y left, z up, right-handed. Both joints rotate about z. The links are l₁ = 0.5 m and l₂ = 0.4 m. The tip position is

x = l1 cos q1 + l2 cos(q1 + q2)
y = l1 sin q1 + l2 sin(q1 + q2)

Differentiating each line with respect to q₁ and q₂ gives the Jacobian, in m/rad:

J(q) = [ -l1 sin q1 - l2 sin(q1 + q2)    -l2 sin(q1 + q2) ]
       [  l1 cos q1 + l2 cos(q1 + q2)     l2 cos(q1 + q2) ]

At q₁ = 0 rad and q₂ = π/2 rad (90 degrees) the tip sits at (0.5, 0.4) m and the matrix is

J = [ -0.4  -0.4 ]
    [  0.5   0.0 ]

Command q̇ = (0.2, −0.1) rad/s. The tip velocity is ẋ = −0.4 × 0.2 + (−0.4) × (−0.1) = −0.04 m/s and ẏ = 0.5 × 0.2 = 0.10 m/s.

Look at the columns geometrically. The first column, (−0.4, 0.5) m/rad, is perpendicular to the base-to-tip vector (0.5, 0.4) m. Its length is 0.640 m, the same as the distance from joint 1 to the tip. Joint 2 sits at (0.5, 0) m, so the tip is (0, 0.4) m away from it, and the second column is (−0.4, 0) m/rad. In both cases the column equals z × r, the cross product of the joint axis with the vector from the joint to the tip. That observation is the whole idea behind the geometric Jacobian.

The determinant is l₁ l₂ sin q₂, which is 0.2 m²/rad² at this pose. It reaches zero when the arm is fully stretched or fully folded. We come back to that below.

Checking with finite differences

A hand derivation is easy to get wrong by a sign. A central finite difference catches that. Column i is approximated by (f(q + h eᵢ) − f(q − h eᵢ)) / 2h. The truncation error shrinks with h², while rounding error grows as roughly 2.2 × 10⁻¹⁶ / h. With h = 10⁻⁶ the two balance near 10⁻¹⁰, so analytic and numeric Jacobians should agree to about ten decimal places.

import numpy as np

L1, L2 = 0.5, 0.4  # link lengths (m)


def fk_planar(q: np.ndarray) -> np.ndarray:
    """Tip position (x, y) in metres for joint angles q in radians."""
    x = L1 * np.cos(q[0]) + L2 * np.cos(q[0] + q[1])
    y = L1 * np.sin(q[0]) + L2 * np.sin(q[0] + q[1])
    return np.array([x, y])


def jacobian_planar(q: np.ndarray) -> np.ndarray:
    """Hand-derived 2x2 Jacobian, units m/rad."""
    s1, c1 = np.sin(q[0]), np.cos(q[0])
    s12, c12 = np.sin(q[0] + q[1]), np.cos(q[0] + q[1])
    return np.array([
        [-L1 * s1 - L2 * s12, -L2 * s12],
        [L1 * c1 + L2 * c12, L2 * c12],
    ])


def numeric_jacobian(f, q: np.ndarray, h: float = 1e-6) -> np.ndarray:
    """Central finite differences, one column per joint."""
    columns = []
    for i in range(len(q)):
        dq = np.zeros_like(q)
        dq[i] = h
        columns.append((f(q + dq) - f(q - dq)) / (2.0 * h))
    return np.column_stack(columns)


q = np.array([0.0, np.pi / 2])
J = jacobian_planar(q)
q_dot = np.array([0.2, -0.1])  # rad/s
tip_velocity = J @ q_dot

rng = np.random.default_rng(0)
worst = 0.0
for _ in range(100):
    q_rand = rng.uniform(-np.pi, np.pi, size=2)
    diff = jacobian_planar(q_rand) - numeric_jacobian(fk_planar, q_rand)
    worst = max(worst, float(np.max(np.abs(diff))))

print(f"det J = {np.linalg.det(J):.3f}, l1*l2*sin(q2) = {L1 * L2 * np.sin(q[1]):.3f}")
print(f"tip velocity = ({tip_velocity[0]:.3f}, {tip_velocity[1]:.3f}) m/s")
print(f"worst analytic vs numeric difference: {worst:.1e}")
print(f"agree: {worst < 1e-8}")

The first two lines print det J = 0.200, l1*l2*sin(q2) = 0.200 and tip velocity = (-0.040, 0.100) m/s, matching the hand calculation. The worst difference over 100 random configurations is of order 1e-10 (the exact digits depend on rounding), and the last line prints agree: True.

The geometric Jacobian of a spatial arm

For a spatial arm the tool has a linear velocity v (m/s) and an angular velocity ω (rad/s). Stacked together they form a twist, a 6-vector. In this lesson the twist is (v, ω) with linear velocity first, v is the velocity of the tool point, and both are expressed in the base frame. The geometric Jacobian is the 6 × n matrix that maps q̇ to that twist.

Each column follows from the planar observation. Let zᵢ be the unit axis of joint i and pᵢ any point on that axis, both in the base frame, and let p be the tool point. For a revolute joint the column is (zᵢ × (p − pᵢ), zᵢ). For a prismatic joint it is (zᵢ, 0), because sliding moves the tool along the axis without rotating it. You need only the forward kinematics to fill it in, so there is no symbolic differentiation at all.

The word “geometric” separates this matrix from the analytic Jacobian, which differentiates a minimal orientation description such as Euler angles. The two differ by a matrix that depends on the angle convention and goes singular at gimbal lock. See rotation representations for why. Controllers almost always use the geometric form.

The test arm

We use a made-up 6-DoF arm, not a commercial product, so there is no datasheet. Its reach from the shoulder to the tool point is 0.85 m and payload is not modelled. Joint 1 yaws about the base z axis. The shoulder (joint 2) sits 0.30 m above the base and pitches about y. The upper arm is 0.40 m long, followed by the elbow (joint 3, pitch). A 0.35 m forearm leads to a spherical wrist: roll about x (joint 4), pitch about y (joint 5) and roll about x (joint 6), all three axes crossing at one point. The tool point is 0.10 m beyond the wrist centre. At q = 0 the arm points straight along +x. A positive pitch tips the arm downward, by the right-hand rule about y.

To check the angular rows numerically we need the derivative of a rotation. If R(t) turns with angular velocity ω in the base frame, then Ṙ Rᵀ is the skew-symmetric matrix of ω. A central difference of R gives Ṙ, and the three off-diagonal pairs give ω.

import numpy as np


def rot(axis: str, angle: float) -> np.ndarray:
    c, s = np.cos(angle), np.sin(angle)
    if axis == "x":
        return np.array([[1.0, 0.0, 0.0], [0.0, c, -s], [0.0, s, c]])
    if axis == "y":
        return np.array([[c, 0.0, s], [0.0, 1.0, 0.0], [-s, 0.0, c]])
    return np.array([[c, -s, 0.0], [s, c, 0.0], [0.0, 0.0, 1.0]])


UNIT = {"x": np.array([1.0, 0.0, 0.0]), "y": np.array([0.0, 1.0, 0.0]), "z": np.array([0.0, 0.0, 1.0])}

# Each joint: (rotation axis in its own frame, offset from the previous joint in metres).
JOINTS = [
    ("z", np.array([0.0, 0.0, 0.0])),   # 1 base yaw
    ("y", np.array([0.0, 0.0, 0.30])),  # 2 shoulder pitch
    ("y", np.array([0.40, 0.0, 0.0])),  # 3 elbow pitch
    ("x", np.array([0.35, 0.0, 0.0])),  # 4 wrist roll
    ("y", np.array([0.0, 0.0, 0.0])),   # 5 wrist pitch
    ("x", np.array([0.0, 0.0, 0.0])),   # 6 tool roll
]
TOOL = np.array([0.10, 0.0, 0.0])  # tool point in the last frame (m)


def fk(q: np.ndarray) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
    """Tool position, tool rotation, joint axes and joint origins, all in the base frame."""
    R = np.eye(3)
    p = np.zeros(3)
    axes, origins = [], []
    for (axis, offset), q_i in zip(JOINTS, q):
        p = p + R @ offset
        axes.append(R @ UNIT[axis])
        origins.append(p.copy())
        R = R @ rot(axis, q_i)
    return p + R @ TOOL, R, np.array(axes), np.array(origins)


def jacobian(q: np.ndarray) -> np.ndarray:
    """Geometric Jacobian: rows 0-2 linear velocity (m/s), rows 3-5 angular velocity (rad/s)."""
    p_tool, _, axes, origins = fk(q)
    J = np.zeros((6, len(q)))
    for i in range(len(q)):
        J[:3, i] = np.cross(axes[i], p_tool - origins[i])
        J[3:, i] = axes[i]
    return J


def numeric_jacobian(q: np.ndarray, h: float = 1e-6) -> np.ndarray:
    _, R0, _, _ = fk(q)
    J = np.zeros((6, len(q)))
    for i in range(len(q)):
        dq = np.zeros(len(q))
        dq[i] = h
        p_plus, R_plus, _, _ = fk(q + dq)
        p_minus, R_minus, _, _ = fk(q - dq)
        J[:3, i] = (p_plus - p_minus) / (2.0 * h)
        S = ((R_plus - R_minus) / (2.0 * h)) @ R0.T  # approximately skew(omega)
        J[3:, i] = 0.5 * np.array([S[2, 1] - S[1, 2], S[0, 2] - S[2, 0], S[1, 0] - S[0, 1]])
    return J


q = np.array([0.3, -0.5, 1.0, 0.4, 0.6, -0.2])
J = jacobian(q)
J_num = numeric_jacobian(q)

q_wrist = q.copy()
q_wrist[4] = 0.0
J_wrist = jacobian(q_wrist)

print(f"max |J - J_num| = {np.max(np.abs(J - J_num)):.1e}")
print(f"agree: {np.allclose(J, J_num, atol=1e-8)}")
print(f"rank at q5 = 0.6 rad: {np.linalg.matrix_rank(J)}")
print(f"rank at q5 = 0 rad: {np.linalg.matrix_rank(J_wrist)}")
print(f"columns 4 and 6 equal at q5 = 0: {np.allclose(J_wrist[:, 3], J_wrist[:, 5])}")

Expect a maximum difference of order 1e-10 and agree: True. The rank is 6 at the test pose and drops to 5 when q₅ = 0 rad, and the last line prints True. With the wrist pitch at zero, joints 4 and 6 share the same axis through the same point. Their columns are identical, so one direction of tool rotation has no joint that can produce it.

Manipulability and singular configurations

The singular value decomposition J = U Σ Vᵀ describes what the arm can do at one instant. Unit-norm joint-rate vectors map to an ellipsoid of tool velocities. Its semi-axes point along the columns of U and have lengths equal to the singular values σ. A singularity is a configuration where at least one σ is zero, so the ellipsoid flattens and some tool direction becomes unreachable.

The manipulability measure w = √det(J Jᵀ) is the product of the singular values. Yoshikawa’s 1985 paper Manipulability of Robotic Mechanisms, in Further reading, is the classic reference. It is proportional to the volume of that ellipsoid. For the square planar Jacobian it reduces to |det J| = 0.2 |sin q₂| m², with a maximum of 0.2 m² at q₂ = 90 degrees. For the spatial arm, rows 0 to 2 are in m/rad and rows 3 to 5 are dimensionless, so a single w mixes units. Compute it separately for the linear and angular blocks, or scale the linear rows by a characteristic length.

Our 6-DoF arm has three families of singularity, and you can find each by asking which columns become dependent:

• Elbow: q₃ = 0 or π rad. The forearm lines up with the upper arm and the wrist centre cannot move along that line.
• Shoulder: the wrist centre lies on the joint 1 axis. Yawing the base then spins the wrist in place without moving it.
• Wrist: q₅ = 0 or π rad. Joints 4 and 6 align, as the code above showed.

What happens to joint speeds near a singularity

Ask the planar arm to move its tip at 0.1 m/s along x with q₁ = 0. Inverting the 2 × 2 Jacobian by hand gives q̇₁ = 0.2 cos q₂ / sin q₂ rad/s and q̇₂ = −(0.25 + 0.2 cos q₂) / sin q₂ rad/s. Both blow up as 1/sin q₂ as the arm straightens. The code below tabulates the larger of the two and adds a damped least-squares solution, explained after the chart.

import numpy as np

L1, L2 = 0.5, 0.4  # link lengths (m)


def jacobian_planar(q: np.ndarray) -> np.ndarray:
    s1, c1 = np.sin(q[0]), np.cos(q[0])
    s12, c12 = np.sin(q[0] + q[1]), np.cos(q[0] + q[1])
    return np.array([
        [-L1 * s1 - L2 * s12, -L2 * s12],
        [L1 * c1 + L2 * c12, L2 * c12],
    ])


def dls(J: np.ndarray, v: np.ndarray, damping: float) -> np.ndarray:
    """Damped least-squares joint rates for a desired tip velocity v."""
    return J.T @ np.linalg.solve(J @ J.T + damping**2 * np.eye(J.shape[0]), v)


v = np.array([0.1, 0.0])  # desired tip velocity (m/s), along x
for deg in [90, 60, 45, 30, 20, 10, 5, 2, 1]:
    J = jacobian_planar(np.array([0.0, np.radians(deg)]))
    w = np.sqrt(np.linalg.det(J @ J.T))
    q_dot = np.linalg.solve(J, v)
    print(f"q2 = {deg:3d} deg  w = {w:.4f} m^2  max|q_dot| = {np.max(np.abs(q_dot)):7.3f} rad/s")

J = jacobian_planar(np.array([0.0, np.radians(1.0)]))
q_dot = dls(J, v, damping=0.05)
print(f"DLS at q2 = 1 deg: max|q_dot| = {np.max(np.abs(q_dot)):.3f} rad/s, "
      f"achieved x speed = {(J @ q_dot)[0] * 1000:.1f} mm/s")

Expected output:

q2 =  90 deg  w = 0.2000 m^2  max|q_dot| =   0.250 rad/s
q2 =  60 deg  w = 0.1732 m^2  max|q_dot| =   0.404 rad/s
q2 =  45 deg  w = 0.1414 m^2  max|q_dot| =   0.554 rad/s
q2 =  30 deg  w = 0.1000 m^2  max|q_dot| =   0.846 rad/s
q2 =  20 deg  w = 0.0684 m^2  max|q_dot| =   1.280 rad/s
q2 =  10 deg  w = 0.0347 m^2  max|q_dot| =   2.574 rad/s
q2 =   5 deg  w = 0.0174 m^2  max|q_dot| =   5.154 rad/s
q2 =   2 deg  w = 0.0070 m^2  max|q_dot| =  12.891 rad/s
q2 =   1 deg  w = 0.0035 m^2  max|q_dot| =  25.783 rad/s
DLS at q2 = 1 deg: max|q_dot| = 0.129 rad/s, achieved x speed = 0.5 mm/s

A modest 0.1 m/s tip speed costs 0.25 rad/s at a comfortable elbow angle and 25.8 rad/s at 1 degree from straight. That is about 246 rpm at the joint output. A controller that inverts the Jacobian blindly will saturate its motors, trip a speed limit or swing the arm violently. The tool barely moves in return, because the direction you asked for has almost vanished from the velocity ellipse.

Damped least squares

For the history, see two 1986 papers: Nakamura and Hanafusa’s Inverse Kinematic Solutions With Singularity Robustness for Robot Manipulator Control (in Further reading) and Wampler’s Manipulator Inverse Kinematic Solutions Based on Vector Formulations and Damped Least-Squares Methods. The damped least-squares solution minimises ‖J q̇ − v‖² + λ² ‖q̇‖². The result is q̇ = Jᵀ (J Jᵀ + λ² I)⁻¹ v. Along a singular direction with value σ the gain becomes σ / (σ² + λ²) instead of 1/σ. That gain never exceeds 1/(2λ). With λ = 0.05 m the cap is 10 rad per metre of commanded displacement, so a 0.1 m/s request cannot cost more than 1 rad/s in that direction.

At 1 degree the damped solution peaks at 0.129 rad/s and the tip moves at only 0.5 mm/s along x. That is the trade. Damped least squares declines the motion that the arm cannot make instead of commanding it at any cost. Far from singularities σ is much larger than λ and the solution is almost identical to the exact inverse. Choosing λ is an engineering decision. You can also raise λ only when the smallest singular value drops below a threshold.

Resolved-rate control

Resolved-rate control, described by Whitney in 1969, moves the tool by commanding its velocity rather than solving for joint angles. At every control step it measures the pose error, turns it into a twist, maps that twist to joint rates through the Jacobian and integrates.

The position error is p_goal − p. The orientation error is the rotation vector of R_goal Rᵀ, the axis-angle form of the rotation that would carry the current tool frame onto the goal, expressed in the base frame. With the twist set to K e and an exact Jacobian inverse, both errors obey ė = −K e in continuous time. They decay exponentially and the tool travels close to a straight line. With K = 2 s⁻¹ the error falls by a factor of e every 0.5 s.

import numpy as np


def rot(axis: str, angle: float) -> np.ndarray:
    c, s = np.cos(angle), np.sin(angle)
    if axis == "x":
        return np.array([[1.0, 0.0, 0.0], [0.0, c, -s], [0.0, s, c]])
    if axis == "y":
        return np.array([[c, 0.0, s], [0.0, 1.0, 0.0], [-s, 0.0, c]])
    return np.array([[c, -s, 0.0], [s, c, 0.0], [0.0, 0.0, 1.0]])


UNIT = {"x": np.array([1.0, 0.0, 0.0]), "y": np.array([0.0, 1.0, 0.0]), "z": np.array([0.0, 0.0, 1.0])}
JOINTS = [
    ("z", np.array([0.0, 0.0, 0.0])),
    ("y", np.array([0.0, 0.0, 0.30])),
    ("y", np.array([0.40, 0.0, 0.0])),
    ("x", np.array([0.35, 0.0, 0.0])),
    ("y", np.array([0.0, 0.0, 0.0])),
    ("x", np.array([0.0, 0.0, 0.0])),
]
TOOL = np.array([0.10, 0.0, 0.0])  # m


def fk(q: np.ndarray) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
    R = np.eye(3)
    p = np.zeros(3)
    axes, origins = [], []
    for (axis, offset), q_i in zip(JOINTS, q):
        p = p + R @ offset
        axes.append(R @ UNIT[axis])
        origins.append(p.copy())
        R = R @ rot(axis, q_i)
    return p + R @ TOOL, R, np.array(axes), np.array(origins)


def jacobian(q: np.ndarray) -> np.ndarray:
    p_tool, _, axes, origins = fk(q)
    J = np.zeros((6, len(q)))
    for i in range(len(q)):
        J[:3, i] = np.cross(axes[i], p_tool - origins[i])
        J[3:, i] = axes[i]
    return J


def rotation_log(R: np.ndarray) -> np.ndarray:
    """Rotation vector (axis times angle, rad). Accurate for small angles, not valid near pi."""
    w = 0.5 * np.array([R[2, 1] - R[1, 2], R[0, 2] - R[2, 0], R[1, 0] - R[0, 1]])
    s = np.linalg.norm(w)  # sin(angle)
    c = (np.trace(R) - 1.0) / 2.0  # cos(angle)
    if s < 1e-12:
        return np.zeros(3)
    return np.arctan2(s, c) / s * w


def dls(J: np.ndarray, twist: np.ndarray, damping: float) -> np.ndarray:
    return J.T @ np.linalg.solve(J @ J.T + damping**2 * np.eye(J.shape[0]), twist)


q = np.array([0.0, -0.4, 1.2, 0.0, 0.6, 0.0])         # start (rad)
q_goal = np.array([0.5, -0.2, 0.9, 0.3, 0.8, -0.4])   # defines a reachable goal pose (rad)
p_goal, R_goal, _, _ = fk(q_goal)
p_start, _, _, _ = fk(q)

gain = 2.0      # 1/s
dt = 0.01       # s, a 100 Hz control loop
damping = 0.01  # m
steps = 800     # 8 s
peak_speed = 0.0

for _ in range(steps):
    p, R, _, _ = fk(q)
    error = np.concatenate([p_goal - p, rotation_log(R_goal @ R.T)])
    q_dot = dls(jacobian(q), gain * error, damping)
    peak_speed = max(peak_speed, float(np.max(np.abs(q_dot))))
    q = q + q_dot * dt

p, R, _, _ = fk(q)
pos_err = float(np.linalg.norm(p_goal - p))
rot_err = float(np.linalg.norm(rotation_log(R_goal @ R.T)))
print(f"initial position error: {np.linalg.norm(p_goal - p_start) * 1000:.0f} mm")
print(f"final position error: {pos_err:.1e} m")
print(f"final orientation error: {rot_err:.1e} rad")
print(f"peak joint speed: {peak_speed:.2f} rad/s")
print(f"converged: {pos_err < 1e-6 and rot_err < 1e-6}")

The first line prints initial position error: 347 mm. Each step multiplies the error by roughly 1 − K Δt = 0.98, and 0.98⁸⁰⁰ is about 1e-7. Both final errors therefore land between 1e-8 and 1e-7 (metres and radians respectively). The exact digits depend on rounding and on the small slowdown from damping, but both sit well below the 1e-6 threshold, so the last line prints converged: True. The peak joint speed line shows what the gain costs. On the first step the commanded joint rates are linear in K, so doubling the gain doubles them.

Three points matter before this loop drives real hardware. First, joint speed and position limits must be enforced, ideally by scaling the whole q̇ vector so the tool keeps its direction. Second, a real controller sends q̇ to a velocity interface at a fixed rate and reads back measured joint angles rather than integrating its own command. Third, the target must stay inside the workspace. Ask for a pose beyond 0.85 m from the shoulder and the loop drives the arm into the elbow singularity, where only the damping keeps joint speeds bounded. The control topic hub covers the servo loops underneath, and the two-link drawing project is a good place to try velocity-level control in simulation.

Exercise

Add a third link of 0.2 m to the planar arm, so the tip Jacobian becomes 2 × 3. Derive it by hand and confirm it with numeric_jacobian to better than 1e-8. Then compute its null space with np.linalg.svd at q = (0, 1.0, −0.5) rad and show numerically that joint rates along that null-space vector leave the tip velocity below 1e-12 m/s. Finally, in the resolved-rate script, change q_goal[4] to 0.05 rad so the goal sits near the wrist singularity. Run it with damping 0.0 m and 0.05 m and record the peak joint speed and final errors for each. Explain the difference in two sentences.

Summary

• The Jacobian maps joint rates to tool velocity. Each column is the tool twist produced by one joint moving at 1 rad/s.
• For revolute joints the geometric Jacobian column is (zᵢ × (p − pᵢ), zᵢ), built from forward kinematics alone. Central finite differences with h = 1e-6 check it to about 1e-10.
• Singular values give the velocity ellipsoid. Manipulability is their product, 0.2 |sin q₂| m² for the 0.5 m plus 0.4 m planar arm.
• Near a singularity the exact inverse demands joint speeds that grow as 1/σ, 25.8 rad/s at 1 degree from straight for a 0.1 m/s tip speed. Damped least squares caps the gain at 1/(2λ).
• Resolved-rate control turns pose error into a twist and integrates damped-inverse joint rates, giving exponential convergence at rate K away from singularities.

Further reading

• Velocity Kinematics and Statics, the chapter in Modern Robotics: Mechanics, Planning, and Control by K. M. Lynch and F. C. Park (Cambridge University Press, 2017). The authors offer a free preprint of the book.
• Robotics: Modelling, Planning and Control by B. Siciliano, L. Sciavicco, L. Villani and G. Oriolo (Springer, 2009). Chapter 3, Differential Kinematics and Statics, is the matching treatment.
• Manipulability of Robotic Mechanisms by T. Yoshikawa, The International Journal of Robotics Research, 1985.
• Resolved Motion Rate Control of Manipulators and Human Prostheses by D. E. Whitney, IEEE Transactions on Man-Machine Systems, 1969.
• Inverse Kinematic Solutions With Singularity Robustness for Robot Manipulator Control by Y. Nakamura and H. Hanafusa, Journal of Dynamic Systems, Measurement, and Control, 1986.