GSoC, Week 6

During week 5, I submitted my pull request. During this past week, I have been working on getting that code to high enough quality to be accepted. It looks like Brian Granger has implemented a change in SymPy’s derivatives in a branch which allows for differentiation of arbitrary things; this includes functions and derivatives, allowing you to take the derivative of an expression with respect to another expression. So in the upcoming week, I’ll be testing how well that can be integrated with physics.mechanics.

This past week, I also started working on the sphinx documentation for the physics.mechanics subsubmodule.  Luke already started a branch for this, so I merged this into my pydy dev branch and got to work. I am currently organizing it in the form: Vector & ReferenceFrame, Kinematics, Masses/Inertias, Forces, and Kane’s Method (and References). This actually follows the order I learned dynamics in my Advanced Dynamics course and also (roughly) the order things were coded in.

So far I filled in the 1st of these, Vector & ReferenceFrame. I first covered all of the math for this, then the code to use it. I did this because the basis vectors are what form other vectors, so you need to have access to them to do any Vector operations. The basis vectors are attributes of ReferenceFrame objects though, which I believe should come after a discussion of Vector (based off of when I learned the concepts).  I also looked through the autolev user’s book; in there they have the example code create reference frames (which automatically creates basis vectors for the user to interact with), and then wait until later for a full discussion of reference frames. I’m not sure which is the better choice for order to discuss things in. Some of the other documentation in SymPy looks more like the autolev manual, with code and math together.  I feel like I might want to switch the order I did Vector & ReferenceFrame in.  Any suggestions here would be nice though.  I think I’ll wait to push the work I’ve done here until I get the chance to organize the commits better, to avoid last week’s issue of changing commit history after pushing; I’ll try and do this on Monday.

I think I’ll also have to learn how to make good figures for this, which might involve learning PSTricks, or something of that nature.  A lot of the information/knowledge that you need to understand what is going on in physics.mechanics needs to be communicated visually.  So far I made some sketches by hand and scanned them, so I know what each image should have and can discuss it appropriately, but the quality is certainly lacking. Some vector graphics here would probably be the best choice.

There’s not that much to report this week due to a day of traveling and a lot of documentation work. I’ve learned that putting a ton of  LaTeX markup in the sphinx documentation can be time consuming. I have the feeling that image creation might be the same way. Fortunately, I have the next 2 weeks to spend on documentation. Hopefully I can keep making progress on my pull request this week too.


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.

GSoC, Week 1

Well, this week was the start of Google Summer of Code 2011.  Unfortunately, the quarter is still not over for me, so I didn’t make as much progress as I would like.  But, it looks like some of the questions about rewriting the Vector class have been answered.

When doing dynamics problems by hand, basis vectors are typically written as:

e_{1,2,3}   or   e_{x,y,z}

Where e is the reference frame within which the basis vectors are “standard”.  One thing that is immediately obvious is that the numbering starts from 1, unlike most programming languages, which start indices at 0.  PyDy is going to be printing out equations to LaTeX, for use in publication; this would most certainly use 1 indexing.  This leads to some conflicting things, where e[0] in the code is e_1 in LaTeX.  This is something I definitely did not want to happen.  So e.x, e.y, and e.z will be how the basis vectors are accessed, removing the indexing question.  There will be a LaTeX printing option to determine 123, xyz, or ijk indices.  But, there will never be 0 as an index (that the end-user will see).

My decision to keep the basis vectors in the ReferenceFrame class, versus generating them upon some function call, ultimately came down to the same issue.  If you named a set of basis vectors b1,b2,b3 , but they were in frame “Foo”, printing them would give foo1>, foo2>, foo3>.  This again, would be unacceptable.  This situation is also fairly likely to come up, as one usually wants to take shortcuts (I know I have when using Autolev).  It also ties into what I guess is my last point for this topic: readability for non-python programmers.

The professor who will teach the graduate multibody dynamics course next year said that he would consider teaching the students PyDy next winter quarter.  This is going to mean that the code should be readable by people who aren’t familiar with Python.  Here, by code, I mean the scripts that people will write to generate their equations of motion.  Anyways, readability for non-Python users, while still trying to be Pythonic, will be an important consideration this summer.

Next is finishing up the Vector and ReferenceFrame classes.  I’m trying to reuse as much code as possible, and definitely not lose any functionality.  It seems to be going OK so far and hopefully it will stay that way.

PyDy and other things

This weekend is Maker Faire, so the lab is going there.  One thing I’ll have is some info on some MEMS sensors from a class project last quarter.  We collected info on the Analog Devices ADXL 345 and Invensense ITG-3200 (tried to use a Honeywell HMC5843, but we couldn’t read it over I2C).  Anyways, I got permission from my group members to post our report online, so here it is.

Report: MAE276 Final Paper

For PyDy/SymPy & GSoC there is also some news.  It looks like the previous separation between UnitVector and Vector classes is going away; in addition they will no longer extend SymPy’s Basic or Expr classes.  There will now be only one Vector class.

It will store a list of lists; the inner list will have the 3 vector measure numbers for a frame (something like [([1],[2],[4]),’b’] ), where the measure numbers will use the SymPy Matrix class (and can have symbols in them).  The outer list will hold each of these lists for each frame.  So the data stored will be something like [[([1],[2],[4]),’b’],[([3*x],[1],[0]),’c’]]… or something like that.  The current plan is to write our own operators for addition, subtraction, scalar multiplication & division, and dot and cross products.  By not using the SymPy classes, we can better control how the Vector class will behave.  One important behavior is to have every operator (except for dot product) take in and return a vector, and not have any confusion between SymPy objects and Vectors when using operators.

There is also some debate as to whether the ReferenceFrame class should store its own basis vectors, or they should be returned upon initialization (see the sympy list post: Question regarding vectors.  I understand the argument to keep them with their frame and the awkwardness that will exist when creating the basis vectors upon ReferenceFrame initialization.  But, I feel that the Vector class already stores information about each frame (in its inner lists), and question whether basis vectors are a property of ReferenceFrames, or a frame is a property of a set of basis vectors.

I guess the orthonormal properties of a set of 3 basis vectors exist independently of a particular reference frame.  I guess what defines a frame?  I think of it as a rotation (and related time derivatives) from on set of basis vectors to another.  So it is defined by both its basis vectors and a rotation.  I still feel like its is preferable to keep the basis vectors outside of the ReferenceFrame class, but not sure my position has a strong enough argument (but I also don’t think it is too weak of an argument).  Hopefully starting to discuss more of the ReferenceFrame class methods will clear this up.