## Using SymPy’s lambdify for generating transform matrices and Jacobians

I’ve been working in VREP with some of their different robot models and testing out force control, and one of the things that becomes pretty important for efficient workflow is to have a streamlined method for setting up the transform matrices and calculating the Jacobians for different robots. You do not want to be working these all out by hand and then writing them in yourself.

A solution that’s been working well for me (and is fully implemented in Python) is to use SymPy to set up the basic transform matrices, from each joint to the next, and then use it’s derivative function to calculate the Jacobian. Once this is calculated you can then use SymPy’s `lambdify` function to parameterize this, and off you go! In this post I’m going to work through an example for controlling VREP’s UR5 arm using force control. And if you’re just looking for code examples, you can find it all up on my GitHub.

Edit: Updated the code to be much nicer, added saving of calculated functions, and a null space controller.
Edit 2: Removed the `simplify` calls from the code, not worth the generation time increase (if execution is a great concern can generate the functions using cython).

Setting up the transform matrices

This is the part of this process that is unique to each arm. The goal is to set up a system so that you can specify your transforms for each joint (and to each centre-of-mass (COM) too, of course) and then be on your way. So here’s the UR5, cute thing that it is: For the UR5, there are 6 joints, and we’re going to be specifying 9 transform matrices: 6 joints and the COM for three arm segments (the remaining arm segment COMs are centered at their respective joint’s reference frame). The joints are all rotary, with 0 and 4 rotating around the z-axis, and the rest all rotating around x.

So first, we’re going to create a class called `robot_config`. Then, to create our transform matrices in SymPy first we need to set up the variables we’re going to be using:

```# set up our joint angle symbols (6th angle doesn't affect any kinematics)
self.q = [sp.Symbol('q%i'%ii) for ii in range(6)]
# segment lengths associated with each joint
self.L = np.array([0.0935, .13453, .4251, .12, .3921, .0935, .0935, .0935])
```

where `self.q` is an array storing all our joint angle symbols, and `self.L` is an array of all of the offsets for each joint and arm segment lengths.

Using these to create the transform matrices for the joints, we get a set up that looks like this:

```# transform matrix from origin to joint 0 reference frame
self.T0org = sp.Matrix([[sp.cos(self.q), -sp.sin(self.q), 0, 0],
[sp.sin(self.q), sp.cos(self.q), 0, 0],
[0, 0, 1, self.L],
[0, 0, 0, 1]])

# transform matrix from joint 0 to joint 1 reference frame
self.T10 = sp.Matrix([[1, 0, 0, -L],
[0, sp.cos(-self.q + sp.pi/2),
-sp.sin(-self.q + sp.pi/2), 0],
[0, sp.sin(-self.q + sp.pi/2),
sp.cos(-self.q + sp.pi/2), 0],
[0, 0, 0, 1]])

# transform matrix from joint 1 to joint 2 reference frame
self.T21 = sp.Matrix([[1, 0, 0, 0],
[0, sp.cos(-self.q),
-sp.sin(-self.q), self.L],
[0, sp.sin(-self.q),
sp.cos(-self.q), 0],
[0, 0, 0, 1]])

# transform matrix from joint 2 to joint 3
self.T32 = sp.Matrix([[1, 0, 0, L],
[0, sp.cos(-self.q - sp.pi/2),
-sp.sin(-self.q - sp.pi/2), self.L],
[0, sp.sin(-self.q - sp.pi/2),
sp.cos(-self.q - sp.pi/2), 0],
[0, 0, 0, 1]])

# transform matrix from joint 3 to joint 4
self.T43 = sp.Matrix([[sp.sin(-self.q - sp.pi/2),
sp.cos(-self.q - sp.pi/2), 0, -self.L],
[sp.cos(-self.q - sp.pi/2),
-sp.sin(-self.q - sp.pi/2), 0, 0],
[0, 0, 1, 0],
[0, 0, 0, 1]])

# transform matrix from joint 4 to joint 5
self.T54 = sp.Matrix([[1, 0, 0, 0],
[0, sp.cos(self.q), -sp.sin(self.q), 0],
[0, sp.sin(self.q), sp.cos(self.q), self.L],
[0, 0, 0, 1]])

# transform matrix from joint 5 to end-effector
self.TEE5 = sp.Matrix([[1, 0, 0, self.L],
[0, 1, 0, 0],
[0, 0, 1, 0],
[0, 0, 0, 1]])
```

You can see a bunch of offsetting of the joint angles by `-sp.pi/2` and this is to account for the expected 0 angle (straight along the reference frame’s x-axis) at those joints being different than the 0 angle defined in the VREP simulation (at a 90 degrees from the x-axis). You can find these by either looking at and finding the joints’ 0 position in VREP or by trial-and-error empirical analysis.

Once you have your transforms, then you have to specify how to move from the origin to each reference frame of interest (the joints and COMs). For that, I’ve just set up a simple function with a switch statement:

```# point of interest in the reference frame (right at the origin)
self.x = sp.Matrix([0,0,0,1])

def _calc_T(self, name, lambdify=True):
""" Uses Sympy to generate the transform for a joint or link

name string: name of the joint or link, or end-effector
lambdify boolean: if True returns a function to calculate
the transform. If False returns the Sympy
matrix
"""

# check to see if we have our transformation saved in file
if os.path.isfile('%s/%s.T' % (self.config_folder, name)):
Tx = cloudpickle.load(open('%s/%s.T' % (self.config_folder, name),
'rb'))
else:
if name == 'joint0' or name == 'link0':
T = self.T0org
elif name == 'joint1' or name == 'link1':
T = self.T0org * self.T10
elif name == 'joint2':
T = self.T0org * self.T10 * self.T21
T = self.T0org * self.T10 * self.Tl21
elif name == 'joint3':
T = self.T0org * self.T10 * self.T21 * self.T32
T = self.T0org * self.T10 * self.T21 * self.Tl32
elif name == 'joint4' or name == 'link4':
T = self.T0org * self.T10 * self.T21 * self.T32 * self.T43
elif name == 'joint5' or name == 'link5':
T = self.T0org * self.T10 * self.T21 * self.T32 * self.T43 * \
self.T54
elif name == 'link6' or name == 'EE':
T = self.T0org * self.T10 * self.T21 * self.T32 * self.T43 * \
self.T54 * self.TEE5
Tx = T * self.x  # to convert from transform matrix to (x,y,z)

# save to file
cloudpickle.dump(Tx, open('%s/%s.T' % (self.config_folder, name),
'wb'))

if lambdify is False:
return Tx
return sp.lambdify(self.q, Tx)
```

So the first part is pretty straight forward, create the transform matrix, and then at the end to get the `(x,y,z)` position we just multiply by a vector we created that represents a point at the origin of the last reference frame. Some of the transform matrices (the ones to the arm segments) I didn’t specify above just to cut down on space.

The second part is where we use this awesome `lambify` function, which lets us turn the matrix we’ve defined into a function, so that we can pass in joint angles and get back out the resulting `(x,y,z)` position. There’s also the option to get the straight up SymPy matrix return, in case you need the symbolic form (which we will!).

NOTE: You can also see that there’s a check built in to look for saved files, and to just load those saved files instead of recalculating things if they’re available. This is because calculating some of these matrices and their derivatives takes a long, long time. I used the `cloudpickle` module to do this because it’s able to easily handle saving a whole bunch of weird things that makes normal `pickle` sour.

Calculating the Jacobian

So now that we’re able to quickly generate the transform matrix for each point of interest on the UR5, we simply take the derivative of the equation for each `(x,y,z)` coordinate with respect to each joint angle to generate our Jacobian.

```def _calc_J(self, name, lambdify=True):
""" Uses Sympy to generate the Jacobian for a joint or link

name string: name of the joint or link, or end-effector
lambdify boolean: if True returns a function to calculate
the Jacobian. If False returns the Sympy
matrix
"""

# check to see if we have our Jacobian saved in file
if os.path.isfile('%s/%s.J' % (self.config_folder, name)):
(self.config_folder, name), 'rb'))
else:
Tx = self._calc_T(name, lambdify=False)
J = []
# calculate derivative of (x,y,z) wrt to each joint
for ii in range(self.num_joints):
J.append([])
J[ii].append(Tx.diff(self.q[ii]))  # dx/dq[ii]
J[ii].append(Tx.diff(self.q[ii]))  # dy/dq[ii]
J[ii].append(Tx.diff(self.q[ii]))  # dz/dq[ii]
```

Here we retrieve the `Tx` vector from our `_calc_T` function, and then calculate the derivatives. When calculating the Jacobian for the end-effector, this is all we need! Huzzah!

But to calculate the Jacobian for transforming the inertia matrices of each arm segment into joint space we’re going to need the orientation information added to our Jacobian as well. This we know ahead of time, for each joint it’s a 3D vector with a 1 on the axis being rotated around. So we can predefine this:

```# orientation part of the Jacobian (compensating for orientations)
self.J_orientation = [
[0, 0, 1], # joint 0 rotates around z axis
[1, 0, 0], # joint 1 rotates around x axis
[1, 0, 0], # joint 2 rotates around x axis
[1, 0, 0], # joint 3 rotates around x axis
[0, 0, 1], # joint 4 rotates around z axis
[1, 0, 0]] # joint 5 rotates around x axis
```

And then we just fill in the Jacobians for each reference frame with `self.J_orientation` up to the last joint, and fill in the rest of the Jacobian with zeros. So e.g. when we’re calculating the Jacobian for the arm segment just past the second joint we’ll use the first two rows of `self.J_orientation` and the rest of the rows will be 0.

So this leads us to the second half of the `_calc_J` function:

```        end_point = name.strip('link').strip('joint')
if end_point != 'EE':
end_point = min(int(end_point) + 1, self.num_joints)
# add on the orientation information up to the last joint
for ii in range(end_point):
J[ii] = J[ii] + self.J_orientation[ii]
# fill in the rest of the joints orientation info with 0
for ii in range(end_point, self.num_joints):
J[ii] = J[ii] + [0, 0, 0]

# save to file
cloudpickle.dump(J, open('%s/%s.J' %
(self.config_folder, name), 'wb'))

J = sp.Matrix(J).T  # correct the orientation of J
if lambdify is False:
return J
return sp.lambdify(self.q, J)
```

The orientation information is added in, we save the result to file, and a function that takes in the joint angles and outputs the Jacobian is created (unless `lambdify == False` in which case the SymPy symbolic form is returned.)

Then finally, two wrapper functions are added in to make creating / accessing these functions easier. First, define a couple of dictionaries

```# create function dictionaries
self._T = {}  # for transform calculations
self._J = {}  # for Jacobian calculations
```

and then our wrapper functions look like this

```def T(self, name, q):
""" Calculates the transform for a joint or link
name string: name of the joint or link, or end-effector
q np.array: joint angles
"""
# check for function in dictionary
if self._T.get(name, None) is None:
print('Generating transform function for %s'%name)
self._T[name] = self.calc_T(name)
return self._T[name](*q)[:-1].flatten()

def J(self, name, q):
""" Calculates the transform for a joint or link
name string: name of the joint or link, or end-effector
q np.array: joint angles
"""
# check for function in dictionary
if self._J.get(name, None) is None:
print('Generating Jacobian function for %s'%name)
self._J[name] = self.calc_J(name)
return np.array(self._J[name](*q)).T
```

So how you use this class (all of this is in a class) is to call these `T` and `J` functions with the current joint angles. They’ll check to see if the functions have already be created or stored in file, if they haven’t then the `T` and / or `J` functions will be created, then our wrappers do a bit of formatting to get them into the proper form (i.e. transposing or cropping), and return you your `(x,y,z)` or Jacobian!

NOTE: It’s a bit of a misnomer to have the function be called `T` and actually return to you `Tx`, but hey this is a blog. Lay off.

Calculating the inertia matrix in joint-space and force of gravity
Now, since we’re here we might as well also calculate the functions for our inertia matrix in joint space and the effect of gravity. So, define a couple more placeholders in our `robot_config` class to help us:

```self._M = []  # placeholder for (x,y,z) inertia matrices
self._Mq = None  # placeholder for joint space inertia matrix function
self._Mq_g = None  # placeholder for joint space gravity term function
```

and then add in our inertia matrix information (defined in each link’s centre-of-mass (COM) reference frame)

```# create the inertia matrices for each link of the ur5
self._M.append(np.diag([1.0, 1.0, 1.0,
self._M.append(np.diag([2.5, 2.5, 2.5,
self._M.append(np.diag([5.7, 5.7, 5.7,
self._M.append(np.diag([3.9, 3.9, 3.9,
self._M.append(np.diag([0.7, 0.7, 0.7,
```

and then using our equations for calculating the system’s inertia and gravity we create our `_calc_Mq` and _calc_Mq_g functions

```def _calc_Mq(self, lambdify=True):
""" Uses Sympy to generate the inertia matrix in
joint space for the ur5

lambdify boolean: if True returns a function to calculate
the Jacobian. If False returns the Sympy
matrix
"""

# check to see if we have our inertia matrix saved in file
if os.path.isfile('%s/Mq' % self.config_folder):
Mq = cloudpickle.load(open('%s/Mq' % self.config_folder, 'rb'))
else:
# get the Jacobians for each link's COM
J = [self._calc_J('link%s' % ii, lambdify=False)

# transform each inertia matrix into joint space
# sum together the effects of arm segments' inertia on each motor
Mq = sp.zeros(self.num_joints)
Mq += J[ii].T * self._M[ii] * J[ii]

# save to file
cloudpickle.dump(Mq, open('%s/Mq' % self.config_folder, 'wb'))

if lambdify is False:
return Mq
return sp.lambdify(self.q, Mq)

def _calc_Mq_g(self, lambdify=True):
""" Uses Sympy to generate the force of gravity in
joint space for the ur5

lambdify boolean: if True returns a function to calculate
the Jacobian. If False returns the Sympy
matrix
"""

# check to see if we have our gravity term saved in file
if os.path.isfile('%s/Mq_g' % self.config_folder):
'rb'))
else:
# get the Jacobians for each link's COM
J = [self._calc_J('link%s' % ii, lambdify=False)

# transform each inertia matrix into joint space and
# sum together the effects of arm segments' inertia on each motor
Mq_g = sp.zeros(self.num_joints, 1)
for ii in range(self.num_joints):
Mq_g += J[ii].T * self._M[ii] * self.gravity

# save to file
cloudpickle.dump(Mq_g, open('%s/Mq_g' % self.config_folder,
'wb'))

if lambdify is False:
return Mq_g
return sp.lambdify(self.q, Mq_g)
```

and wrapper functions

```def Mq(self, q):
""" Calculates the joint space inertia matrix for the ur5

q np.array: joint angles
"""
# check for function in dictionary
if self._Mq is None:
print('Generating inertia matrix function')
self._Mq = self._calc_Mq()
return np.array(self._Mq(*q))

def Mq_g(self, q):
""" Calculates the force of gravity in joint space for the ur5

q np.array: joint angles
"""
# check for function in dictionary
if self._Mq_g is None:
print('Generating gravity effects function')
self._Mq_g = self._calc_Mq_g()
return np.array(self._Mq_g(*q)).flatten()
```

and we’re all set!

Putting it all together

Now we have nice clean code to generate everything we need for our controller. Using the controller developed in this post as a base, we can replace those calculations with the following nice compact code (which also includes a secondary null-space controller to keep the arm near resting joint angles):

```# calculate position of the end-effector
# derived in the ur5 calc_TnJ class
xyz = robot_config.T('EE', q)
# calculate the Jacobian for the end effector
JEE = robot_config.J('EE', q)
# calculate the inertia matrix in joint space
Mq = robot_config.Mq(q)
# calculate the effect of gravity in joint space
Mq_g = robot_config.Mq_g(q)

# convert the mass compensation into end effector space
Mx_inv = np.dot(JEE, np.dot(np.linalg.inv(Mq), JEE.T))
svd_u, svd_s, svd_v = np.linalg.svd(Mx_inv)
# cut off any singular values that could cause control problems
singularity_thresh = .00025
for i in range(len(svd_s)):
svd_s[i] = 0 if svd_s[i] < singularity_thresh else \
1./float(svd_s[i])
# numpy returns U,S,V.T, so have to transpose both here
Mx = np.dot(svd_v.T, np.dot(np.diag(svd_s), svd_u.T))

kp = 100
kv = np.sqrt(kp)
# calculate desired force in (x,y,z) space
u_xyz = np.dot(Mx, target_xyz - xyz)
# transform into joint space, add vel and gravity compensation
u = (kp * np.dot(JEE.T, u_xyz) - np.dot(Mq, kv * dq) - Mq_g)

# calculate our secondary control signal
# calculated desired joint angle acceleration
q_des = (((robot_config.rest_angles - q) + np.pi) %
(np.pi*2) - np.pi)
u_null = np.dot(Mq, (kp * q_des - kv * dq))

# calculate the null space filter
Jdyn_inv = np.dot(Mx, np.dot(JEE, np.linalg.inv(Mq)))
null_filter = (np.eye(robot_config.num_joints) -
np.dot(JEE.T, Jdyn_inv))
u_null_filtered = np.dot(null_filter, u_null)

u += u_null_filtered
```

And there you go!

You can see all of this code up on my GitHub, along a full example controlling a UR5 VREP model though reaching to a series of targets. It looks something pretty much like this (slightly better because this gif was made before adding in the null space controller): Overhead of using `lambdify` instead of hard-coding

This was a big question that I had, because when I’m running simulations the time step is on the order of a few milliseconds, with the controller code called at every time step. So I reaaaally can’t afford a crazy overhead for using `lambdify`.

To test this, I used the handy Python `timeit`, which requires a bit awkward setup, but quite nicely calls the function a whole bunch of times (1,000,000 by default) and accounts for various other things going on that could affect the execution time.

I tested two sample functions, one simpler than the other. Here’s the code for setting up and testing the first function:

```import timeit
import seaborn

# Test 1 ----------------------------------------------------------------------
print('\nTest function 1: ')
time_sympy1 = timeit.timeit(
stmt = 'f(np.random.random(), np.random.random())',
setup = 'import numpy as np;\
import sympy as sp;\
q0 = sp.Symbol("q0");\
l0 = sp.Symbol("l0");\
a = sp.cos(q0) * l0;\
f = sp.lambdify((q0, l0), a, "numpy")')
print('Sympy lambdify function 1 time: ', time_sympy1)

time_hardcoded1 = timeit.timeit(
stmt = 'np.cos(np.random.random())*np.random.random()',
setup = 'import numpy as np')
print('Hard coded function 1 time: ', time_hardcoded1)
```

Pretty simple, a bit of a pain in the `sympy` setup, but other than that not bad at all. The second function I tested was just a random collection of `cos` and `sin` calls that resemble what gets computed in a Jacobian:

```l1*np.sin(q0 - l0*np.sin(q1)*np.cos(q2) - l2*np.sin(q2) - l0*np.sin(q1) + q0*l0)*np.cos(q0) + l2*np.sin(q0)
```

And here’s the results: So it’s slower for sure, but again this is the difference in time after 1,000,000 function calls, so until some big optimization needs to be done using the SymPy `lambdify` function is definitely worth the slight gain in execution time for the insane convenience.

The full code for the timing tests here are also up on my GitHub.

## 16 thoughts on “Using SymPy’s lambdify for generating transform matrices and Jacobians”

1. Hrishik Mishra says:

Hi.

I tried to find the kinematic link lengths from the UR5 model in VREP. On the forums, I found this: http://www.forum.coppeliarobotics.com/viewtopic.php?f=9&t=207

Could you provide how you found those out? My values are a little different from the ones you used.

• travisdewolf says:

Hello! To find link lengths the way I’ve been doing it is to select the object of interest, and click the ‘object/item shift’ button. Once the dialogue pops up for that you can choose ‘parent frame’ and that will show you the offset from the centre of that item to the centre of its parent. So I just add up those distances, from a joint to the arm segment to the next joint and that’s how I get the link lengths from joint to joint!

2. […] could go about this, but the way I’m going to do it here is by building on the post where I used SymPy to automate the Jacobian derivation. The way we did that was by defining the transforms from the origin reference frame to the first […]

3. Boris Skegin says:

Hello.

Actually I use OpenRave, but I found you script nevertheless very useful.

Just one question:
# transform matrix from joint 5 to end-effector
self.TEE5 = sp.Matrix([[0, 0, 0, self.L],
[0, 0, 0, 0],
[0, 0, 0, 0],
[0, 0, 0, 1]])

Should not the rotational part of this affine transformation be a Identity Matrix?

Regards,
boris

• travisdewolf says:

gah! yes it should, I fixed it on github (i believe) but forgot to make the change here. Thanks for the catch, and glad you found the post useful!

4. Boris Skegin says:

Due to OpenRave I have to use sympy 0.7.1.
sp.simplify often yields:
ValueError: can’t raise polynomial to a negative power

sp.trigsimp does not complain.

• travisdewolf says:

Hm, Calling simplify runs all of the functions here http://docs.sympy.org/latest/tutorial/simplification.html in addition to trigsimp. I wonder which simplify function is causing the error!

• Boris Skegin says:

This might be very much model dependent as I use a different (5 DOF) roboter, and the sympy version I am confined to is 0.7.1 , not 1.0 as you probably use.

So no worries. I was just curious if you saw smth. similar.

5. Boris Skegin says:

Sorry for disturbing you again. Just one more question:

Do I guess it right that in order to calculate the rotational part of the Jacobian one has to
use self.x = sp.Matrix([1,1,1,0]) ?

• travisdewolf says:

Hi Boris, I’m not sure I follow…the code is set up to calculate the Jacobian through the partial derivative of `Tx`, where `x=[0, 0, 0, 1]` (at the origin of the last frame of reference). If you were to calculate `Tx` with `x=[1, 1, 1, 0]` that would be the Jacobian for a point offset by 1m along (x, y, z) from the last reference frame’s origin. Does that answer your question?

• Boris Skegin says:

Indeed my question was a bit confusing. Actually I wanted to know if one can derive
the EE Jacobian for the pose orientation from T = self.T0org * self.T10 * self.T21 * self.T32 * self.T43 *self.T54 * self.TEE5 as easy as it is the case for the pose position.

• travisdewolf says:

Hi Boris, yup you sure can! Once you have the rotation matrix for the end-effector you can convert it to a quaternion, and then use the methods as described in section 4 of this paper: http://bit.ly/2jvNtFR
to calculate the error between current and desired orientation and off you go!

6. wagglefoot says:

Hey there! Cracking example. I seem to have uncovered an error with this example, but I left the details on your git instead of here, as it is a bit extensive.

• travisdewolf says:

great, thanks! i’ll look it over / correct it asap!

7. wagglefoot says:

It would appear, lest I have made a mistake, that the null-controller’s gain term \$k_{p_null}\$ as not been included?

• travisdewolf says:

ah, correct! I’ve used the same `kp` for the primary and null controllers above for compactness, but in practice I’ve found `kpn = .1 * kp` to be a value that works well (and this is how it’s set in the code).