GSoC 2019 Application Animesh Sinha: Classical Mechanics Composite Rigid Body Algorithm - sympy/sympy GitHub Wiki

Introduction

Animesh Sinha

First Year student perusing Bachelors in Computer Science and Masters in Computational Natural Sciences at International Institute of Information Technology, Hyderabad.

Contact Information

Me as a Programmer

Why did I choose this project

Classical Mechanics and Robotics are both very close to my heart. I have loved Physics forever, as a school student I had managed to get into the Indian astronomy olympiad camp and I spent a great deal of time learning Astrophysics,which essentially boiled down to solving Mechanics, orbital and stellar dynamics problems. And on the computation side I start off programming with and still continue coding on Microcontrollers to control small sensor/actuator systems, "robots".

I feel that this aligns with my current interests and is an opportunity for me to learn some new theoretical ideas on how control systems and classical mechanics gel together.

But more important than anything else, I studied Classical Mechanics only studied for it's beautiful ways of computing things and how it's ideas stretch into quantum, but this project made me see a new side of the subject, how it's relevant and beautiful today in a more modern, computational form, and I really want to explore it.

Academic experience in this domain

Currently I am pursuing a Bachelors degree in Computational Natural Science and have undergone a semester both in Linear Algebra and Classical Mechanics, so my intuition is relatively fresh and decent in both of these. I am also currently undergoing a course in these computational techniques which is based around implementing mathematical and scientific models in Python scripts.

Preferred Editor

I oscillate between Visual Studio code and JetBrains PyCharm (PyCharm being my love for all Python Projects). JetBrains IDEs manage to take a lot of routine tasks off my hands and they speed me up. They are indispensable when starting on a new project, it's very convenient even as a beginner to navigate a large code base with ease and ensure that the my code quality is decent.

Experience with Python

I started learning Python from the MIT Course 6.00.2x on edX, over 6 years back. Since then I have implemented some django websites, taken the UC Berkeley Introductory AI course and made the Pacman projects in python, and had quite a few other projects. And I have read some books, Automate the Boring Stuff with Python being one of my Favourites.

Over the past year, in college I have had courses specifically target at scientific computation using the Sympy, Scipy, Numpy stack. I have tried a bit to contribute to other Python projects as well, BeeWare's Voc being one.

Why Python?

Its interoperability with C++. C++ gives the speed and control, while Python provides the highest order of convenience of coding, this lands a perfect combo. Despite other languages having an arguably more solid C interface, Python is hands down my favourite here because of the stark differences between the priorities of Python and C++, and that both are integral to every project.

Favourite Feature of Sympy

It's hard to find out a favourite feature, I need many of them on and off. One module that I do not often use, but I find incredibly interesting is the physics.quantum, specially because this offers a very different capability set from Microsoft Q#.

The Project

Problem Statement

Implementing the utility classes required for the Composite Rigid Body Algorithm for Forward Dynamics, and then to implement the algorithm itself.

Description of the Stages and Modules

We need to Find a method for forward Dynamics such that: $$\ddot{q} = FD(model, q, \dot{q}, \tau)$$ Solvable in $O(n)$ time. Or more descriptively in $O(n \times d^2)$ time where the depth of the constraints tree d is bounded by a constant value.

To do this I shall start off by writing some auxiliary classes that help find convenient notation and some level of abstraction. The motion of a rigid body and equations relating to it can be more conveniently described as vectors embedded in 6-D vectors space that combine linear and angular aspects called the Spatial Coordinates, developing these would be the first step.

Next would be describing the system and subsystems. We need to collect equations of motion by writing down the obvious ones for bodies/sub-systems, then combining them for bigger systems and applying additional constraints of motion.

Methods for computing joint constraints would be needed. Each joint can be of several types, revolute, helical, cylindrical, planar, etc, and the constraints can be explicit or implicit. This method would give us the Joint transform $X^{s}_{p}$, which is the transform from the predecessor frame to the successor frame.

Finally I will start on the meat of the CRBA Algorithm. We are solving the following equation: $$\tau = H(q) \ddot{q} + C(q, \dot{q}, f^x)$$ This would need 3 steps:

  1. Calculating C, the vector of Centrifugal and Coriolis forces. $O(n)$
  2. Calculating H, the inertia matrix. $O(nd)$
  3. Solving the equation above to get $\ddot{q}$. $O(nd^2)$

The first two of these steps would be performed by calling the Inverse Dynamics methods (Kane's or Newton-Lagrange). One call to these methods can get the vector C, and two calls can get the product of H with any vector, allowing us to recover H column by column.

This would leave me only with final stage writing the equations in body coordinates and collecting the equations of motion. Following is the pseudocode of the algorithm (from Featherstone, Chapter-6 Page-107):

H = 0
for i in range(1, N_b):
    I_c[i] = I[i]
for i in range(N_b, 0, -1):
    if Lambda(i) != 0:
        I_c[L(i)] += X_s(L(i), i) * I_c[i] * X(i, L(i))
    F = I_c[i] * S[i]
    H[i][i] = np.transpose(S[i]) * F
    j = i
    while Lambda(j) != 0:
        F *= X(Lambda(j), j)
        j = Lambda(j)
        H[i][j] = np.transpose(F) * S[j]
        H[j][i] = np.transpose(H)[i][j]

where the _c is used as superscript c in the equations representing that it's expressed in body coordinates.

Project Timeline

  • Bonding Period - Reading up on Books and Papers Go through Rigid Body Dynamics Algorithms (Roy Featherstone). Discussing the API that would be needed in all the proposed classes and the way objects would be represented (overview of code structure).
  • Week 1 - Draft up the Implementation plans of the Helper classes: This would be for Spatial Vector Algebra (Plucker Notation). Write out the exact set of API this class will need to support. Writing tests for the same. These would include testing dot products, Cross products, addition of like vectors, differentiation, etc.
  • Week 2 - Implement classes for Spatial Vector Algebra: Computation of Spatial velocity, Spatial forces, Spatial inertia from the same physical quantity in regular 3-D space with linear and angular components. This should support Addition, Dot Products, Cross Products, etc.
  • Week 3 - Complete Implementation of the Spatial Vector Algebra module: Specifically for the Dyadic representation of the Inertia matrix which may differ from the other vectors. Completing any left-over segments.
  • Week 4 - Calling Inverse Dynamics methods to Get H and C: Starting off by digging into the current codebase and finding the calls to the Inverse Dynamics functions (Newton-Euler method implemented in sympy.physics.lagrange).
  • Week 5 - Figuring out the interface: Deciding what function will be made available to the user and exactly how the user will specify the System. Coding up this interface and collecting the data I would need.
  • Week 6 - Finishing off finding H and C: Finishing off the implementation and Writing tests. Any catch up time that may be required for the previous modules.
  • Week 7 - Representation of the Objects Planning how the elements of the kinematic tree will be represented, how the equation of motion will be stored, and how to get all of this from the input.
  • Week 8 - Joint Model Calculation routines Writing methods that computes explicit and implicit joint constraints. Method explained in chapter 4 of Featherstone.
  • Week 9 - Draft the Implementation of the Composite Rigid Body Algorithm: Map out the specifics of what methods will be needed and how the code will be structured.
  • Week 10 - Start on the CRBA: Start out on the final implementation of the algorithm. This would solely be focused on implementing the exact methods in Chapter 6 of the book, some part of the pseudocode for which is above
  • Week 11 - Finish implementing the CRBA. Additional time needed to finish off the work in the previous week. This may also involve some work in improving the user interface for the module.
  • Week 12 - Speed benchmarks, additional tests and Code cleanup: Wrapping up the project, cleanup the code and test if what the speed of the algorithm is, finishing up if anything is left.

This I believe will take up the full 12 weeks, and if at all I manage to finish this off earlier than targeted, I would like to take up the implementation of the vector algebra classes in SymEngine.

Deliverable and Pull Request Plan

  • Spatial Vector Algebra class supporting API for Addition, Dot and Cross Products, Differentiation, etc. as well as initialising a Spatial vector from two ordinary vectors. (In 4 weeks, before 1st evaluation)
  • Joint Model transform calculator methods (By week 8), and the method computing the Joint Space Inertia Matrix (By Week 10).
  • Full equations of motion module (By the end of the project)

Pull Request Plans

I would complete the code for the Vector Algebra classes before the 1st evaluation (in the first 4 weeks). Next would be code for the representing objects in the Kinematic tree with their equations of motion, and applying the different types of Joint constraints on them.

References

Major Resources

Additional Links

Portions of Current Codebase used

The implementation of the Forward Dynamics algorithm would include multiple calls to the already implemented methods in the physics.lagrange module as a GSoC project in 2016.

Time Commitments

The full 40 hours will be possible throughout the summer with no breaks in the Schedule. This would be my sole focus through the summers, so if needed, I would be more than happy to spend extra hours and work Weekends if I am unable to complete my targets.

I would also keep contributing to the Classical and Quantum Mechanics modules well after this GSoC cycle is over, and will further and finish anything left on this project.

Pull Requests

To SymPy

  • #16292: Adding Hamiltonian Method for Equations of Motion OPEN: Implemented a function that takes in a Hamiltonian and returns a list of equations of motion.
  • #16527: Breaking up the Rational Factors OPEN: Split the numerator and denominator when factorizing a term, fixing the collection of constants issue.
  • #16281: Made Product pretty-print neater MERGED: Minor change to the printing of PI for products.
  • #15668: Removed special handling of sqrt(I) MERGED: Stopped the oversimplification of expressions like i^{1/2} to -1^{1/4}.
  • #15635: Printing of Fancy Sets MERGED: Adding a few lines to print and pretty print S.Integers, S.Naturals, etc.
  • #15622: Docstrings improved in calculus/singularities.py MERGED: Updated documentation to Numpydoc format.

Other Projects

  • BeeWare/Voc #935: Making sets work with complex elements MERGED: Implementing Python hash_functions in the Java transpiler to give identical behavior in working of sets across types.
  • BeeWare/Voc #937: Hashing Functions OPEN: Updating the hashing function to let complex numbers be put in sets.
  • AimaCode/Aima-JavaScript #171: Implementing Wumpus World OPEN: Implementing a JavaScript simulation of a AI Logical agent tring to reason it's way through a game grid.
  • Processing/p5js-website #343 to #366 ES6 Standard Port MERGED: Updating JavaScript code to ES6 standards.