Archive for the ‘Sympy’ Category

New examples of use with PyDy

Sunday, August 23rd, 2009

In 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.

Tags: ,
Posted in PyDy, Sympy | No Comments »

Recursive methods for computing equations of motion

Tuesday, July 28th, 2009

In 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.

Tags: , , , ,
Posted in PyDy, Sympy | No Comments »

GSoC developments

Monday, July 13th, 2009

In 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.

Tags: , , ,
Posted in PyDy, Sympy | No Comments »

Point class and doctest

Monday, July 6th, 2009

When 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

Tags: , , ,
Posted in PyDy, Sympy | No Comments »

Code cleanup and updating of examples

Tuesday, June 30th, 2009

A 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.

Tags: ,
Posted in PyDy, Sympy | No Comments »

Dyad and Inertia classes implemented

Sunday, June 7th, 2009

In order to be able to specify the inertia of a rigid body, I implemented the Dyad and Inertia classes. The Dyad class works with the UnitVector class and supports left and right dot multiplication by any other UnitVector or Vector expression (assuming the relative orientations have been specified). The Inertia class is subclassed from Dyad, and allows the user to specify the reference frame, along with the 6 inertia scalars that fully define the mass distribution of a rigid body.

An example of use for the Dyad class would be:

>>> from sympy import *
>>> from pydy import *
>>> t = Symbol(‘t’)
>>> q1 = Function(‘q1′)(t)
>>> A = ReferenceFrame(‘A’)
>>> B = A.rotate(‘B’, 1, q1)
>>> var(‘a b c’)
(a, b, c)
>>> a_Dyad = Dyad(a*A[1]*B[2] + b*A[3]*B[1] + c*A[2]*B[3])
>>> print a_Dyad
Dyad(a·a₁·b₂ + b·a₃·b₁ + c·a₂·b₃)

Most of the time in classical mechanics, dyads are used to keep track of the inertia of rigid bodies, so I also implemented an inertia class:

>>> I_dyad = Inertia(A, (1, 2, 3, 4, 5, 6))
print I_dyad
Dyad(a₁**2 + 2·a₂**2 + 3·a₃**2 + 4·a₁·a₂ + 4·a₂·a₁ + 5·a₂·a₃ + 5·a₃·a₂ + 6·a₁·a₃ + 6·a₃·a₁)

I need to polish the display of the dyad and inertia classes, as they currently don’t display quite right (a₁**2 should be a₁*a₁), but the mathematics is there and behaves as it should:

>>> print dot(A[1], a_Dyad)
a·b₂
>>> print dot(a_Dyad, A[1])
b·b₁

In addition to these classes and class methods, I also implemented a couple of convenience functions: InertiaForce() and InertiaTorque(), which simply form -m*a and -dot(alpha, I) – cross(omega, dot(I, omega)), respectively.

These new functions will be incorporated into the frame work that will allow for the derivation of the equations of motion for a system of particles and rigid bodies.

My finals will be done by Tuesday, I’m looking forward to some serious work on PyDy once exams are over :)

Tags: , ,
Posted in PyDy, Sympy | No Comments »

Progress on trigsimp(), .eval() method of sin, cos, tan, and PyDy printing

Tuesday, June 2nd, 2009

A few things I’ve been working on in the last week have been the .eval() methods for the sin, cos, and tan classes inside of sympy.functions.elementary.trigonometric. Each of these class methods control how something like

>>> sin(x + pi/2)

is evaluated. Ideally, this would nicely evaluate to cos(x), by default. Similar results are desirable for cos and tan. This is the behavior that Mathematica has. So far, I have implemented this behavior into the sin and cos functions, and should have tan done in the next few days. A few snags along the way have been dealing with the case when some of the arguments are Function or Derivative types, rather than Symbol, but I think I have figured this out. This will be a nice improvement of the default Sympy behavior and will make implementing trigsimp() easier.

Once the behavior of Table 1 of the paper by Fu et al. (see last post) has been implemented into these .eval() methods of sin, cos, tan, my next task will be to implement the rule based algorithm for trigsimp(). The first step will be to be able to compare the length of two trigonometric expressions, so I will be writing a function that will act as a metric for comparing trig expressions. Once this is done, one can compare a trig expression before and after a trig identity has been applied to it, and if the new expression is smaller, according to the metric, then we can accept it, otherwise we reject it. The algorithm is very customizable because you can try any sequence of trig identities to see if one approach works better than another, although hopefully the sequences of trig identities in the paper will prove to be sufficient for most problems.

Other developments on the PyDy front include adding a subs method to the Vector class, for the purpose of replacing time derivatives of the generalized coordinates with the expressions involving the generalized speeds. This was a simple matter of looping through the expressions corresponding to each of the UnitVector keys in the dictionary that is used to represent each Vector quantity.

The PyDy printer has been customized and now it by default will display sin(q1) as s1, cos(q1) as c1, etc. This really makes things more readable. I had a lot of issues getting unicode to work properly, but now the multiplication ‘*’ symbol is being displayed with a nice small unicode dot, and the display looks a lot cleaner:

(-r1·q5′ – r1·q3′·s4 – r2·q5′·c4)·b₁ + (r1·q4′ + r2·q4′·c4)·b₂ – r2·q4′·s4·b₃

The above output is the velocity of the mass center of the rolling torus example that I am using to test many of the functions I implement in PyDy. What isn’t visible above is that the b₁, b₂, and b₃ symbols are in bold and bright white, because they are unit vectors, so they stand out make the output easy to read. Another idea I have been toying with is to make the time derivative of the generalized coordinates a different color, and similarly with the time derivatives of the generalized speeds. Additionally I have customized the printer to display derivatives with an apostrophe after them, which is a lot easier to read and takes up a lot less space than D(q1(t), t), which was the default Sympy behavior. What would be really cool was if there was a Unicode dot syntax, like \dot{} in LaTeX, then things would be really small. Creating a custom LaTeX printer is on my list of things to do this summer, so that all of the equations that are need to include in a paper can be generated in a LaTeX form that is consistent with the syntax used in the dynamics literature.

Tags: , ,
Posted in PyDy, Sympy | 3 Comments »

Body (Euler) fixed / Space fixed angles implemented in PyDy

Monday, May 25th, 2009

The rotate() method of the Vector() class in PyDy now supports the 12 Euler angle rotations and the 12 space fixed rotations. Additionally, it automatically calculates the angular velocity of the new reference frame with respect to the parent frame. Here is how it works:

>>> t = Symbol(‘t’)
>>> q1 = Function(‘q1′)(t)
>>> q2 = Function(‘q2′)(t)
>>> q3 = Function(‘q3′)(t)
>>> A = ReferenceFrame(‘A’)
>>> B = A.rotate(‘B’, ‘BODY123′, (q1, q2, q3))
>>> dot(A[1], B[3])
sin(q2(t))
>>> dot(A[1], B[2])
-cos(q2(t))*sin(q3(t))
>>> dot(A[1], B[1])
cos(q2(t))*cos(q3(t))
>>> B = A.rotate(‘B’, ‘SPACE321′, (q1, q2, q3))
>>> dot(A[1], B[3])
sin(q2(t))
>>> dot(A[1], B[2])
-cos(q2(t))*sin(q1(t))
>>> dot(A[1], B[1])
cos(q1(t))*cos(q2(t))
>>>

The A[i], B[i] (i=1,2,3) notation indicates the i’th unit vector fixed in each of the reference frames, with the usual right hand rule applying: cross(A[1], A[2]) == A[3]

All 24 angle set conventions are implemented, and tests have been written to ensure they return the correct result. Rotations about an arbitrary axis, as well as the use of Rodrigues parameters or Euler parameters remain to be implemented, but they should be done in the next week.

This code has highlighted the need for some serious work to be done on Sympy’s trigsimp(). Especially when calculating the angular velocity vectors, things get messy quickly and unless the trigonometric simplification routines are up to the task, things will get unnecessarily messy. In the next week I will be working to implement the rule based algorithm presented by Fu et al.:

Automated and readable simplification of trigonometric expressions

Rigid body mechanics are notorious for long nasty trigonometric expressions, so implementing this will really help Sympy and how it deals with trigonometric expressions.

Tags: , ,
Posted in PyDy, Sympy | No Comments »