Comment by abricq

Working in the field of robotics, I often have to write big vectors (usually computed as the result of 3D transformations) and then compute their Jacobian (their derivative with respect to several state-variables), which quickly becomes very nasty equations. I use sympy to (i) compute these big vectors, expressed in a very declarative way, (ii) compute the jacobian and (iii) export the results in C-code, immediately importable in my code case.

To illustrate what I mean by "expressing systems of equations in a declarative way", here is a toy-example of how to estimate the position of the a sensor with respect to a robot's center if you have access to a dataset containing robot positions and the sensor positions. The 'naive' approach (it very often works) is to solve an over-constrained system with a gradient-descent. To perform a gradient descent, you need a residual function and its jacobian. Here's how you would do to compute it with Sympy. (Note: you'd just have to define the `transform` and `invert` functions...)

    # Pose of a sensor in robot frame (to be estimated)
    xa, ya = symbols("xa, ya")
    a = Matrix([xa, ya, ta])

    # Position of the robot at time t
    rx, ry, rt = symbols("rx, ry, rt")
    rk = Matrix([rx, ry, rt])
    
    # Measure of the sensor at time t
    gx, gy = symbols("gx, gy")
    gk = Matrix([gx, gy, 0])
    
    # Estimated x (from the measures)
    estimated_a = transform(invert(rk), gk)
     
    # compute the norm of the gk, squared
    n2_mat = norm2(estimated_a - a)
    n2 = sympy.collect(sympy.expand(n2_mat[0, 0]), a).simplify()
    
    # Compute the jacobian
    J = n2_mat.jacobian([xa, ya, ta])
    
    # Print what is necessary for Guass-Markov regression
    print("\n\nres =", n2)
    print("\nJacobian = ", J)