GSoC, Week 5

At the end of last week, a lot of the functionality I was aiming for was completed.  I think the code can now compute equations of motion.  I’m still working on some more advanced tests, but for simple cases I know it works. I submitted a pull request for pulling my work into sympy:master.  It involved cleaning up a lot of ugly commit history.  I believe the last time I tried to rebase, Luke and I both rebased and pushed, leading to a lot of double commits.  This was obviously bad.  So I rebased again to remove the double commits, squashed a lot of the small commits, and added a lot of comments to commits. I think now that a lot of the desired functionality is present future commits won’t be all over the place, and I’ll do a better job with commit messages.

One thing that is important on github, and in bright letters on the rebasing help page, is that if you rebase and edit the commit history and someone else has pulled from your branch pre-commit, then they will have a bad day when they pull next…. So I would recommend avoiding that. I wish I could say I felt more comfortable with git, but every time I try and do something new, I have trouble.  I guess that means I’m learning though, even though I am scared of git still…

During the week, Luke and I discussed formulation of Fr* and Fr in a different manner when using Kane’s Method. Typically R* = -ma, where a is the translational acceleration vector of a body’s mass center. Now when you consider the final form of the equations, Fr + Fr* = 0, you have all the terms on one side. For analysis, often this will be arranged in a form: M(q) udot = f(u,q,t,…). M(q) is the mass matrix, and it is the coefficients of the accelerations. Now, accelerations are the time derivatives of velocity, and can be represented in the form: c udot + g(q,u,t,…). This c matrix is the start of contributions to the mass matrix from the accelerations (it needs to be multiplied by the body’s mass and partials). This matrix is actually in the form v1, v2, v3, … where each v is the partial velocity for the respective u dot.  This is a little confusing, so I’ll try and tie this back a step.  Velocities can be defined as follows: the sum of (V_r* u_r) + other terms. Vr is said to be the partial velocity of V with respect to generalized speed u_r. Then there is a remainder term.

Next when you take the derivative of the velocity, and use the chain rule, each element in the sum is evaluated (in a frame, N) as ^N d/dt(V_r*u_r) which is ^N d/dt(u_r) * V_r + ^N d/dt(V_r) * u_r.  The d/dt(u_r) * V_r term is then V_r * Udot_r.  Then we can see that it is the partial velocities which are the elements of the “c” matrix in the acceleration definition. This goes back the acceleration definition:  c udot + g(q,u,t,…). The next step when forming the generalized active force, Fr*, is to take multiply R* by the V_r for each generalized speed. It can then be seen that this will be of the form [v_1, v_2, v_3,…]. Next is taking the “outer product” of this with itself.  Outer product is not really correct here, but it leads to a matrix in the form -m * [[v1&v1, v1&v2, v1&v3],[v2&v1,v2&v2,v2&v3],[v3&v1,v3&v2,v3&v3]] where & represents the dot product (remember, partial velocities are still vectors). So it is this matrix which is the mass matrix.  What I’ve gone over is just the mass matrix contribution from the translational components of a body; it extends as you would expect.  Also, all those additional terms end up in the “f” term in our EoM formulation: M(q) udot = f(u,q,t,…).

This is really confusing, but I hope by adding some more text to what is in Kane’s book (, that this is more visible.  He covers everything in there, but it can be really hard to follow.  I’m going to be writing all of this up in the documentation for PyDy in the next few weeks, and I’ll be using LaTeX & mathjax with the sphinx documentation so all these things will show up correctly. I hope by that point, I’ll be able to describe these things better (and I’m sure displaying actual math will help).

Finally, a little example of how this will work.  This is a rolling disc, infinitely thin, with no slip constraints.  This is actually taken right from the test code, so you can copy, paste, and run it (assuming you have a branch of sympy with my pydy code).

    # Rolling Disc Example
    # Here the rolling disc is formed from the contact point up, removing the
    # need to introduce generalized speeds. Only 3 configuration and three
    # speed variables are need to describe this system, along with the disc's
    # mass and radius, and the local graivty (note that mass will drop out).
    q1, q2, q3, q1d, q2d, q3d = dynamicsymbols('q1 q2 q3 q1d q2d q3d')
    u1, u2, u3, u1d, u2d, u3d = dynamicsymbols('u1 u2 u3 u1d u2d u3d')
    r, m, g = symbols('r m g')

    # The kinematics are formed by a series of simple rotations. Each simple
    # rotation creates a new frame, and the next rotation is defined by the new
    # frame's basis vectors. This example uses a 3-1-2 series of rotations, or
    # Z, X, Y series of rotations. Angular velocity for this is defined using
    # the second frame's basis (the lean frame).
    N = ReferenceFrame('N')
    Y = N.orientnew('Y', 'Simple', q1, 3)
    L = Y.orientnew('L', 'Simple', q2, 1)
    R = L.orientnew('R', 'Simple', q3, 2)
    R.set_ang_vel(N, u1 * L.x + u2 * L.y + u3 * L.z)
    R.set_ang_acc(N, R.ang_vel_in(N).dt(R) + (R.ang_vel_in(N) ^ R.ang_vel_in(N)))

    # This is the translational kinematics. We create a point with no velocity
    # in N; this is the contact point between the disc and ground. Next we form
    # the position vector from the contact point to the disc mass center.
    # Finally we form the velocity and acceleration of the disc.
    C = Point('C')
    C.set_vel(N, 0)
    Dmc = C.newpoint('Dmc', r * L.z)
    Dmc.v2pt(C, N, R)
    Dmc.a2pt(C, N, R)

    # This is a simple way to form the inertia dyadic.
    I = inertia(L, m / 4 * r**2, m / 2 * r**2, m / 4 * r**2)

    # Kinematic differential equations; how the generalized coordinate time
    # derivatives relate to generalized speeds.
    kd = dict({q1d: (u3/cos(q3)),
        q2d: (u1),
        q3d: (u2 - u3 * tan(q2))})

    # Creation of the force list; it is the gravitational force at the mass
    # center of the disc. Then we create the disc by assigning a Point to the
    # mass center attribute, a ReferenceFrame to the frame attribute, and mass
    # and inertia. Then we form the body list.
    ForceList = [(Dmc, - m * g * Y.z)]
    BodyD = RigidBody() = Dmc
    BodyD.inertia = (I, Dmc)
    BodyD.frame = R
    BodyD.mass = m
    BodyList = [BodyD]

    # Finally we form the equations of motion, using the same steps we did
    # before. Specify inertial frame, supply generalized speeds, supply
    # kinematic differential equation dictionary, compute Fr from the force
    # list and Fr* fromt the body list, compute the mass matrix and forcing
    # terms, then solve for the u dots (time derivatives of the generalized
    # speeds).
    KM = Kane(N)
    KM.gen_speeds([u1, u2, u3])
    fr = KM.form_fr(ForceList)
    frstar = KM.form_frstar(BodyList)
    MM = KM.mass_matrix()
    forcing = KM.forcing()
    rhs = MM.inv() * forcing
    print rhs.expand()

GSoC, Week 4

Things went pretty well this week.  Last weekend, I implemented a Dyad class.  We decided to do this for a number of reasons; but largely it was to allow the end-user the freedom to express a body’s rotational inertia as sums of components in different frames.  More information about dyads can be found here:, or in Kane’s 1985 book: Dynamics, Theory and Applications.  A dyad (or dyadic, I guess?) represents the juxtaposition of two vectors.  My experience with dyadics has involved treating them as a matrix, basically; but implementing this class showed me some cool things they can do. I also made a pretty simple convenience function for creating an inertia dyadic by supplying 3 or 6 scalar values and a frame; this will allow users to treat the inertia dyadic like they would treat the inertia tensor; but they could still write the more complicated Dyad expressions out if they wanted.  Also, the outer product was implemented for Vector, allowing for easy creation of Dyads.

So, Dyad was the last main class which is used to describe physical parameters/quantities/etc. After that, RigidBody and Particle classes were made.  Luke and I discussed how we implement Kane’s Method when doing it by hand, and realized that basically we do all of the kinematics, then just write down tables for storing the partial velocities and forces.  We both came to the realization that the terms “particle” and “rigid body” were nothing more than associating a mass/inertia with a point/frame, and that forces are really just a Vector associated with a Point or ReferenceFrame.  So, we decided to make RigidBody and Particle container classes effectively, where they both have attributes (and getters/setters for sanitizing input) for all their relevant information.  Then, I decided that storage of forces didn’t have to be anything more than a list of tuples, in the form (point/frame, force/torque).  I think we’ll make some more convenience functions here, probably relating to gravitational forces.

Finally, Kane’s Method.  Luke has done a lot of thinking about this, so a lot of credit goes to him for the successful implementation of this.  When one does Kane’s Method by hand (or at least when I do it), I follow a series of steps; something along the lines of:

  1. Set up generalized coordinates and speeds
  2. Do kinematics
  3. List forces
  4. Form relevant partial velocities
  5. Calculate Fr, Fr*
  6. Rewrite in desired form
So, I’ve basically written a class which stores the relevant information as you go through these steps. Here’s a little example:
KM = Kane(N)
KM.gen_speeds([u1, u2, u3, u4, u5, u6])
KM.dependent_speeds([u1, u4, u6], conl)
fr = KM.form_fr(FL)
frstar = KM.form_frstar(BL)

First, you create the object specifying the inertia frame.  Then you supply the generalized speeds and kinematic differential equations, which relate all of the q’s (gen. coords.) to u’s (gen. speeds) by means of qdot = f(u).  This is stored in a dictionary, and I need to figure out a better/easier way for the user to generate these for bodies oriented by Euler Angles (perhaps another convenience function?).  Fr is the generalized active force, and is formed from the list of forces/points and their partial velocities, Fr* is the generalized inertia force, and is formed from the list of particles/bodies and their partial velocities.

Something that Luke has thought about for a while is dealing with nonholonomic systems; these are systems with velocities constraints (example: ideal skate, where it can move forward and backward and it can yaw, but never move left and right).  Kane describes it in his 1985 book, but it’s a little unclear; the way Luke has described is the same thing, but a little easier conceptually.  This will get described in a lot more detail in the documentation for PyDy.  One other cool thing which I wanted to implement (so I did) was dealing with systems with time varying mass.  Newton’s second equation (f = m a) is really f = d/dt (m v); or force = time derivative of linear momentum.  This also extends to angular momentum/torques.  I make sure to check the time derivative of the mass or inertia, and if it is zero, do things the “normal” way; if it is time varying, I calculate the momenta and take the time derivatives thereof.  So, my code should work with systems with non-constant mass/inertia.  Finally, one more cool thing.  Sometimes, a system will be defined by parameters which are time-varying, but not dynamics (ie user-specified speeds or positions or such).  These will be defined as DynamicSymbols instead of Symbols.  Then, when there are DynamicSymbols which are not generalized speeds or coordinates, the code will find and identify them.  The plan is for the SciPy output code to have empty functions for these parameters for the user to then fill in.  One thing I haven’t thought about (mainly because I have never done any examples) is flexible dynamic systems; I’d really like for my code to work with flexible systems, but I’ll need to try a few first.  Hopefully it will transfer nicely into the framework I have built so nothing new will have to be implemented.

So basically, Kane’s Method has been written.  The functionality that was supposed to be complete this week was supposed to include SciPy code generation, but I’m a little behind due to my late start.  I feel like I’ve made a good amount of progress though…… The schedule for next week has creation of equations of motion unittests, so I’ll probably work out a bunch of examples (or take them from Kane’s book) and implement them as tests, in addition to implementing SciPy code generation.  The week after is integration into SymPy, so I’ll probably start talking to Aaron this week about that.  Then it’s examples, documentation and LaTeX output.  Next week I’ll have some examples of generated EoM to show.

GSoC, Week 3

Another week of Summer of Code progress.

This week, there was a lot less new functionality added, and a lot more refinement of behaviors (and more documentation added).  The code got split into multiple files over the weekend, to match most of the other SymPy modules’ organization.  Tests were added for a lot of things; but more are still needed.  I’ve tried to write tests for the most basic functionality, and make that functionality more robust.

Something interesting comes up when testing though; frequently, outputs are generated which are equal to the “reference value” in the test file, but they don’t equate to True in the SymPy code.  I’ve been trying to make compromises between using expand()/simplify() (actually, I am using trigsimp(), due to the nature of the expressions generated), and getting more complicated expressions to work.  The further out you go when creating ReferenceFrames (A is rotated from N, B is rotated from A….), the more complicated the direction cosine matrix between the two frames are.  I’ve tried to keep the examples in the test files to what I think is a reasonable limit…. which is somewhat arbitrary, but is based on experience I’ve had with more complicated multibody system.  It’s important to note that the expressions generated are still valid; they just don’t compare correctly.

The other big thing was switching to using dictionaries to represent positions, orientations, velocities, etc. By storing a dictionary in each frame of it’s DCM (or ang_vel, or whatever) relative to another frame, AND in that other frame the inverse DCM to this frame, immediate comparisons between the two are simpler.  But, it also has the advantage of making it possible to have code of a reasonable length to find the shortest number of intermediate frames between two frames (this is a lot better, and the code is simpler and shorter too).  Basically, knowing the start and end frame, lists are made from the start frame to every other frame using the dictionary entries as the next possible steps (with no backtracking).  Then, lists that don’t end with the end frame are thrown away.  Finally, we take the shortest list.  This is nice as it should make the expression as simple as possible.  Also, it sets things up so that if the user defines, say, angular velocity between two non-adjacent frames, that this defined angular velocity will be used instead of the auto-generated angular velocity, as it will have a shorter list.  That was one important behavior; always use the user defined values over auto-generated ones.  That’s why the dictionary has been so helpful; you just update the key: value pair and everything is good. It’s also worked out well for angular velocities, as now they are generated on DCM formation and stored in an intelligent manner, and are overwritten by the user if desired.

So, it looks like the code is more robust now, and can scale to more complicated systems (I hope).  Next is better code for velocities & accelerations of points, writing the body & particle classes, then finally starting to implement Kane’s Method as an algorithm to get the equations of motion.  Oh, and more tests.  And documentation. And getting Sphinx to work. And ….

GSoC, Week 2

My blog post is a little late this week…

Good progress was made this week, with lots of code added to the kinematics file, starting a tests file, and starting the inertias file.  I think most of the kinematics classes are done…. but more tests need to be implemented, and then the code corrected if there are any errors.  It’s looking like the inertias file will be pretty easy to fill out, as those classes are extending a lot of what has already been written.  The next code that will be a challenge to figure out is the algorithm implementations (e.g. Kane’s Method, Lagrange’s Method, etc.).  In forming the equations of motion, expressions can get very big very fast, so we will have to take care here.

I think the best piece of news was from yesterday.  Throughout the week, I had been trying to figure out how to be able to take the partial derivative of a SymPy expression with respect to a generalized speed, without using any substituions (that the user would see, at least).  The problem was that SymPy’s diff can only take in symbols, and for an undefined function (say, x(t)), you get a Derivative object returned (say, D(x(t), t)).
At first I tried extending Symbol, unsuccessfully, to return a new symbol when its derivative was taken.  This didn’t work once my extended Symbol was inside a Add or Mul, as the new methods I had written weren’t called.  I posted this to the mailing list (, and after some discussion, someone suggested that I add Symbol(‘t’) to the free_symbols property of my extended Symbol class.  I also got the name for the class, DynamicSymbol, from this discussion.  Basically, now when you take the time derivative of a DynamicSymbol object, even if it is within other SymPy objects, it’s _eval_derivative() method is called, and now it will return what you want.  Here’s a little example:

In [1]: run

In [2]: x = TVS(‘x’)

In [3]: y = TVS(‘y’)

In [4]: t = Symbol(‘t’)

In [5]: diff(2 * x**2 + 4 + y,t)
Out[5]: 4*x*xdt + ydt

In [7]: diff(diff(2 * x**2 + 4 + y,t),TVS(‘xdt’))
Out[7]: 4*x

This is pretty exciting, as now things will work fairly intuitively, when it comes to derivatives and partial derivatives for generalized coordinates and speeds.  Also, if you have time-varying symbols (maybe some specified position), it will now be identified as such, and when the final equations of motion are formed with the output code, it can be treated as a variable defined by a user-specified function.

Next is filling out the inertia classes, adding more tests and ensuring correct functionality, and starting on the algorithm implementations.