Solving for the Jacobians of a Robot Leg

I'm currently working on the Python code to control a simulated version of my latest bipedal robot design in PyBullet. My focus over the past few weeks was getting the operational space control to work (many thanks to Travis DeWolf's incredibly helpful blog). After finally getting it to work properly, I have decided to share my math in the hopes of providing a useful example for anyone else having trouble with this. There really aren't enough resources on the internet that explain this in a succinct manner.

An earlier mechanical design, for reference.

Kinematic diagram. (Here the body is simply assumed to be a fixed base).

Solving for the transformation matrices and centers of motion correctly is the trickiest part of the process, and deceptively so. For more information on how to set up the transformation matrices, I recommend this Youtube tutorial and this blog post by Travis DeWolf.

I find it easiest to solve for the transformation matrices by "stringing" the robot out, as shown below, such that all of the joints' coordinate systems are oriented the same way. This prevents you from having to perform additional linear algebra operations that can introduce errors.

Shown below are the transformation matrices. The first transformation represents an x-axis rotation. See the first joint in the kinematic diagram above for reference. Transformations (2) through (4) are z-axis rotations. Despite the z-axis pointing up in the world coordinate system, the default angle (90 degrees) for the first x-axis rotation points the z-axis of the second through fourth joints parallel to the y-axis of the WCS.
​

Next, the locations for the centers of mass at each link are found, where l₀, l₁, l₂, and l₃ represent the distance of the center of mass from the base joint of each link, along a line connecting the base joint to the next joint. Keep in mind that the centers of mass may not be located directly along a line from one joint to the next in any given robot link, and that in this case I am able to make this approximation for the purpose of mathematical simplicity (but that may not always be the case).

The end-effector offset (Equation 9) should also be kept in mind, but in my case I don't need an offset from the end-effector. As such, I'm expressing it as zero.

Next, I solved for the full transformation matrices from the origin to each joint. This isn't necessary for the first joint,  because it's already in the base coordinate system.

Now, the Jacobian for each COM and the end-effector can be found as a matrix of its partial derivatives with respect to each joint qᵢ. Both linear and angular velocities must be accounted for, so I have x, y, z, ωx, ωy, and ωz  to deal with.

So far I haven't explained angular velocity. This I have solved separately, and it's quite simple as long as your robot is serial and doesn't have spherical joints. The partial derivative for the angular velocity w.r.t. each joint is expressed as 1 if that joint rotates along that axis, and 0 if it doesn't. For example, in J₀ (Equation 14), joint q₀ rotates around the x axis, so the partial derivative of ωx w.r.t. q​₀ is 1. In J₁ (Equation 15), joint q₁ rotates about the z axis, so the partial derivative of ωz w.r.t. q​1 is 1.

I can only show the first two Jacobians here, as the rest are too long to fit on the screen. You'll notice that the Jacobians for the centers of motion of each link (J₀, J₁, J₂, and J₃) are calculated as a function of the transformation matrix from the link's base joint to the origin multiplied by that link's center of motion. This confused me the first time around, and I don't want you, the reader, to make the same mistake.

For example, in the case of J₁ (the Jacobian for the second link), the transformation matrix from joint 0 to the world coordinate frame (origin) is multiplied by COM₁, the center of motion of the link between joints 0 and 1 (the first and second joints, so the second link).

For  J₀ (the Jacobian for the first link), no transformation matrix is required because the first joint (q₀) is already in the base coordinate system.

I'm probably just confusing you more because my joint indexing starts at zero, but that's for programming purposes and is standard...

So there you have it. After several headaches, I was able to verify that this setup works in simulation. The chief cause of my problems was that I had solved for the transformation matrices incorrectly! Everything else is easy--though it looks intimidating, all of the symbolic math can be computationally solved by Matlab/Octave. Just remember, garbage in = garbage out.