I changed the code to the rolling torus example pretty significantly. Little fixes in how the arrows representing the body fixed unit vectors are now fixed. Additionally, I added some code that calculates the kinetic, potential and total energy of the system, and plots it. The scipy solver seems to do a pretty good (not perfect) job of maintaining constant energy. I think tightening up the error tolerances would improve this, but the simulation is well behaved at this point so I didn’t bother.
New rolling torus code
October 15th, 2009Gyrostat
October 9th, 2009I just added an example of use for the derivation of the equations of a gyrostat using PyDy. It is located in the examples/gyrostat/ subdirectory.
In case you aren’t familiar, here is the definition of a gyrostat:
A gyrostat is a mechanical system which is comprised of more than one body and yet has the rigid body property that its inertia components are time independent constants.
The reason I became interested in gyrostats is because in modeling the bicycle, the rear frame and the rear wheel can be treated as a gyrostat, and in doing so, the number of parameters that appear in the equations of motion will be reduced by two. The same can be done for the front fork and handlebar assembly and front wheel.
In addition to the Python script that generates the equations, I took the time to do a complete write up of the model in LaTeX and generate a nice pdf of it.
New examples of use with PyDy
August 23rd, 2009In the last several days I’ve made a lot of updates to tools for using PyDy and have come very close to settling on the way that PyDy is used to derive the equations of motion. I have added and updated several examples, including the double pendulum and a rigid body with two reaction wheels used for attitude control. Animations for each example have been added using visual python, so visualization of the dynamics is very tangible.
There are two main things that I am still trying to iron out with PyDy. The first is how to handle ignorable coordinates and how to allow for the user to control whether or not they are included in the output equations of motion. For example, for the symmetric rolling disc, there are 4 ignorable coordinates: two for the location in the ground plane, one for the heading and one for the spin. For purposes of stability analysis, the kinematic differential equations associated with these coordinates are not needed. However, for purposes of animation, these equations are needed. My goal is to make it easy to clearly specify whether or not these equations are desired in the output equations.
The second major thing that still needs work is handling dependent generalized speeds in nonholonomic systems. When the motion equations are generated, they will involve these dependent generalized speeds, and their time derivatives, but these quantities can be computed from the constraint equations and therefore can be left implicit in the final equations of motion. The derivatives of the dependent generalized can be intelligently computed by some careful formations of the gradients of the components of the matrix that relates the dependent speeds to the independent ones.
I will be working on both of these tasks this week and will be writing examples for the Whipple bicycle model in addition to a spinning top and the rattleback.
Developments in PyDy
August 12th, 2009I haven’t blogged in about two weeks, and there has been a lot of progress in PyDy. For starters, PyDy now ‘automatically’ derives the equations of motion for a rolling disc and a rigid body in space. I haven’t blogged in about two weeks, and there has been a lot of progress in PyDy. For starters, PyDy now ‘automatically’ derives the equations of motion for a rolling disc and a rigid body in space. One of the more useful methods I implemented is the output_eoms() method in the NewtonianReferenceFrame class. Once the equations of motion have been derived, this method creates python functions for them (in the form that scipy.odeint wants them in) and writes them to a file that you specify. Once you have this, you can import the function and use it for numerical integration. Additionally, this method also outputs a function that calculates the generalized speeds given the coordinates and their time derivatives, which is useful if you know the initial conditions of the coordinates time derivatives but need to determine the initial conditions of the generalized speeds when you want to perform a numerical integration. Finally, this method also optionally outputs an animate function, which takes a variable number of Points and (Axis, Angle) tuples and outputs a function which given the coordinates and parameters, will compute the inertial position and orientation of Points and/or Reference Frames. This is useful when you want to animate your simulation.
Both the rolling_disc.py and rigidbody.py files in the pydy/examples folder demonstrate the usage of these new methods.
One issue that has been repeatedly rearing its ugly head is the use of Function instances instead of Symbol instances. For everything except time differentiation, PyDy expressions could work fine with Symbols instead of Functions. I have become to realize that everything just works faster and more efficiently in Sympy with Symbols, so I will be working in the next week or so to replace the use of Functions in PyDy with Symbol instances instead. This should allow lots of things to work much better without having to call a bunch of helper methods to get expressions in their simplest form, or having to repeatedly substitute and back substitute with Symbols and Functions. For time differentiation, I will keep track of a list of symbols that represent the generalized coordinates, and then any expression that needs time differentiation will be examined to see if it has any occurrences of those coordinates. For expressions involving those coordinates, the partial derivative will be taken with respect to each coordinate, and then multiplied by the symbol that will be used to represent the time derivative of that coordinate, in essence a manual (as opposed to letting Sympy take care of it by using sympy Function) application of the chain rule. The user will be required to eclare these coordinate symbols in a special way, specifically using the .setcoords() method. Finally, I will implement method that will replace the symbols with Function instances (depending on time), so that if one wants to play with an expression interactively, all that they will have to do is call this function so that things like q1 = Symbol(’q1′) get replaced with q1 = Function(’q1′)(t), and then time differentiation will work just fine.
Recursive methods for computing equations of motion
July 28th, 2009In the last week, I added methods to the NewtonianReferenceFrame class of PyDy which start at the top of the reference frame tree and the point tree and do several things. The first method, recursive_subs(expr), is designed to take kinematic differential equations (which arise through a combination of linear nonholonomic constraints and/or holonomic constriants, and definitions of generalized speeds) and substitute them into all expressions for relative angular velocity between ReferenceFrames and relative velocity between Points. In doing so, the angular velocity and velocity expressions will involve only the independent time derivatives of coordinates. The second method, recursive_acc(), parses the Point/ReferenceFrame trees and forms the relative angular acceleration and acceleration. The third method, apply_gravitational_forces(), allows one to apply a gravitational forces to every particle and rigid body in the system with one command.
Currently, the NewtonianReferenceFrame class is acting like a container class for information about the system which is somewhat global in nature. For example, definitions of generalized coordinates, their time derivatives, and generalized speeds must all be assigned as attributes of this class in order for it to work correctly.
On the todo list is to write fr() and frstar(), which will also be recursive functions which will compute the equations of motion in one command. More work on the functions dealing with constraints is also needed, and there are at least two approaches that I could potentially implement for this. In particular, the approach taken with nonholonomic systems versus holonomic systems may necessitate different approaches. Currently, the work flow is as follows:
1) Form angular velocity and velocity expressions (linear in time derivatives of coordinates)
2) Form constraint equations arising from differentiated holonomic constraints and/or simple motion (rolling) constraints (m equations, linear in time derivatives of coordinates).
3) Define n-m generalized speeds as linear combinations of time derivatives of coordinates.
4) Solve generalized speed equations (item 3) for n-m ‘independent’ time derivatives of coordinates, in terms of the n-m generalized speeds and the m ‘dependent’ time derivatives of coordinates.
5) Solve constraint equations for m ‘dependent’ time derivatives of coordinates in terms of the ‘independent’ time derivatives of coordinates.
6) The complete set of kinematic differential equations will be the equations arising from item 5 and item 4.
7) Substitute these expressions into the expressions for angular velocity and linear velocity so that the constraints and definitions for generalized speeds are imposed, and so that the constrained partial velocities and partial angular velocities can be formed.
In taking this approach, the minimal set of n-m dynamic equations of motion will be formed, and they will obey all the constraints. These n-m equations, combined with the n kinematic differential equations comprise the complete set of motion equations that describe the system.
PyDy’s examples are being updated to reflect this new functionality, especially the rolling torus example. Additionally I have added a bicycle model which should soon have the nonlinear equations of motion for a very standard bicycle model.
New trig functions
July 20th, 2009Last week I spent most of my time writing and testing new trigonometric functions for Sympy. Sympy already had sin, cos, tan, cot, asin, acos, atan, acot, but did not have sec, csc, asec, or acsc. So I added those four new functions and am in the process of writing very extensive tests for them.
The reason for writing these new functions is so that the trig functions have the same behavior as the ones in Mathematica. There are a few things that Mathematica does that is very nice, one of them being that any arguments that have a pi shift (as in sin(a + b*pi)), will always return another trig function whose pi shift is within the interval between 0 and pi/4. For example,
>>> sin(11*pi/8)
-cos(pi/8)
>>> x = Symbol(‘x’)
>>> sin(x + 11*pi/8)
-cos(-x + pi/8)
>>> sin(acot(x))
1/(x*(1 + x**(-2))**(1/2))
>>> sin(pi/10)
-1/4 + 5**(1/2)/4
-cos(pi/8)
>>> x = Symbol(‘x’)
>>> sin(x + 11*pi/8)
-cos(-x + pi/8)
>>> sin(acot(x))
1/(x*(1 + x**(-2))**(1/2))
>>> sin(pi/10)
-1/4 + 5**(1/2)/4
Getting all the behavior to work with the inverse trig functions is also implemented, as well as the special values of the trig functions for arguments which are integer multiples of pi/10. The next task is going to be to get all of this code working with the rest of Sympy so it can be put into place.
This behavior is essential to making the trigonometric simplification routine in Sympy to be more capable, and will greatly benefit both Sympy and PyDy.
GSoC developments
July 13th, 2009In the last week, I finalized the locate method for the point class, as well as the .rel() and .vel() methods. Just as relative angular velocity between adjacent ReferenceFrame instances is stored, so is the relative position and relative velocity of points. By doing it in this fashion, the partial velocities can be formed incrementally, and there are some nice efficiencies to be realized when it comes time to form the equations of motion by avoiding the need to take redundant dot products. This will also translate into natural definitions for intermediate variables in the final form of the equations of motion.
One subtlety that is encountered in the kinematic chain of points is the possibility for closed kinematic chains. Typical examples are four and six bar linkages, or a bicycle with both wheels on the ground. These closed kinematic loops introduce holonomic constraint equations which are generally not solvable directly, but can be dealt with in differentiated form. PyDy now has a function called kinematic_chain() that deals with such loops and automatically forms the constraint equations and the differentiated holonomic constraint equations, and its use is illustrated in examples/fourbarlinkage.py in the git repository.
I have neared the completion of my updated sin,cos,tan,cot,sec,csc functions that handle arbitrary rational shifts by pi, and am in the process of testing and making things work with the hyperbolic functions that already exist in Sympy. This has been very challenging to get all the behavior that Mathematica does, but I think I have finally converged upon a solution that will be functional. Once this behavior is done, my next task will be to implement the Fu et al. paper so that trigsimp works. I will be putting PyDy development on hold until this is complete because PyDy needs to have some good trig simplification routines for many things to simplify as they should.
Point class and doctest
July 6th, 2009When I implemented the ReferenceFrame class with Ondrej, things were fairly straightfoward because the angular velocities of frames are 1) always described relative to another frame, and 2) because they add up in according to the angular velocity addition theorem. Points, on the other hand have proven a bit more complex to implement. When the user instantiates a new point using the locate method of the Point class, they provide the position Vector from the parent point to the new point. Optionally they may provide a ReferenceFrame which will determine how the velocity of the new point is determined (more on this later). Finding the position of one point relative to another is straightforward, and is implemented in the .rel() method of the Point class. What isn’t so straightforward is how to implement a velocity method to the point class.
When forming the velocity of a point, say point P1, one has to take the time derivative of a position vector of P1. This begs two questions:
1) Which relative position vector should be differentiated?
2) What frame should this position vector be differentiated in?
For forming the equations of motion, the answer to the first question is: the position from the inertial origin (or any other inertially fixed point) to the point in question. Similarly, the answer to the second question is: the inertial reference frame. This may not be what you want to do *ALL* of the time, but in every example I could think of, this is what is needed to form the equations of motion. The next question is how to do this efficiently.
The structure of the ReferenceFrame class is such that upon instantiation, ReferenceFrame objects have both an orientation and angular velocity that is relative to their parent frame. For the Point class, I have implemented the same thing — each new point has a position relative to it’s parent point. In both cases, this generates a tree structure with an inertial (Newtonian) frame at the top of the tree (in the case of the ReferenceFrame) and an inertial origin at the top of the tree (in the case of the Point class).
When the angular velocity of a frame is desired, one can use the .ang_vel() method of the ReferenceFrame class, and one needs to provide another ReferenceFrame whcih that angular velocity is computed relatve to (as in, the angular velocity of frame A relative to frame B). When the velocity of a point is desired on the other hand, what information should the user specify? My initial answer to this question is to make a method call .vel() with 2 default arguments that are None, and compute the velocity of the point with respect to the inertial origin unless otherwise specified. The issue with this is whether to pre-compute the time derivative (in the inertial frame) of each relative position vector, so that all one has to do is add the relative velocities of every point between the point in question and the inertial origin, so that the inertial velocity of that point is obtain. Additionally, when one uses the locate method, the option ReferenceFrame argument that can be specified will determine how this relative velocity is computed. Without specifying the optional argument direct time differentiation (in the inertial frame) of the position vector is performed, while specifying the optional argument results in the velocity being formed by omega (of the frame provided, relative to the inertial frame) crossed with the position vector. Ensuring that this subtlety is taken into account when implementing the .vel() method is a little tricky and that is what I’ve been working on for the last few days.
The Sympy documentation day last week was a great oppurtunity to write some documentation for how to subclass the StrPrinter, and it also exposed me to Sphinx and doctest. I started putting doctests in all of my docstrings, and I also started playing with Sphinx to figure out how to generate some nice documentation. I will be putting this up on docs.pydy.org in the next week.
I also spent a bunch of time learning how to manage a VPS. There is an excellent video blog series on www.guvnr.com where he goes through step by step (from zero to hero, as he puts it) how to setup a VPS to run the webserver nginx, a wordpress blog, etc. I’ve gotten the webserver up and running (now pydy.org is being redirected to the Google group page by my server instead of Ondrej’s), and I also got a nice secure ftp server up and running (using vsftpd in inetd mode with SSL). All in all very useful, but very time consuming!!! I’m exploring options for a wiki, probably either mediawiki or plone, but not sure yet. I managed to get plone up and running as well, but it looks like it might be more than is needed, plus I think it is very resource intensive.
Finally, over the last couple months I’ve been formally studying the book Programming in Python 3: A Complete Introduction to the Python Language by Mark Summerfield. I’ve learned about bunch of things I didn’t know, including named tuples, default dictionaries, and dictionary comprehensions (very handy for my Vector class), and I’m just about to dig into the modules chapter. Reading this book has really given me a better understanding of things like decorators (although they are still a little confusing), and also how to structure a module/package.
~Luke
Code cleanup and updating of examples
June 30th, 2009A few days ago, I discovered the .get() method of the dictionary data type. It allows you to request the value corresponding to a key which a dictionary may not have, in which case you can specify what should be returned instead of rasing a KeyError. In PyDy, I had numerous cases where I was using if/else statements to check if a dictionary had a certain value, but by using the .get() method, I was able to eliminate about 50 or 60 lines of code and make the code a lot more readable and easier to maintain.
The progress on the trigonometric functions has been slow. Ondrej and I worked on some code together that handles most of the cases in the table by Fu et al., but not all. One of the things that has been a little challenging is to identify pi shifts in the argument to the trig functions and map the shift into the first quadrant (0, pi/2) and return the appropriate result, e.g.: sin(x+17*pi/18) –> sin(x + pi/18), etc. The approach has been to match x and r in the following generic argument:
x + r*pi
As I see it, there are three cases to consider: 1) x==0, r!=0, 2) r==0, x!=0, and 3) x!=0, and r!=0. I am splitting the logic up in this fashion and am trying to deal with each case in the appropriate manner. In most cases, r will be of the Rational class, and therefore the modulo operator is needed to bring it to the interval (0, 2*pi), so I implemented the necessary code and sent a patch to the Sympy mailing.
Although the trig functions aren’t quite as they need to be yet, I was still able to manually get the simplifications that are presented by Fu et al. to work using the builtin Sympy commands, which hopefully will mean that simplifying the implementing the algorithm will be feasible with the existing Sympy tools and not too many extra functions will need to be created.
Still working on trigonometric.py
June 22nd, 2009I am still working on modifying trigonometric.py to make all of the trig methods behave as they would in Mathematica/Matlab/Maple. This involves getting the .eval method of each of them correct, and then from there I can work on implmenting a new trigsimp.