Posts Tagged ‘PyDy’

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 »

Still working on trigonometric.py

Monday, June 22nd, 2009

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

Tags: , ,
Posted in Uncategorized | No Comments »

NSF Success!

Wednesday, June 17th, 2009

Today I found out that our NSF proposal to the control systems division was accepted! This means I will be funded for two years to study and work bicycles and control! This is a dream come true for me and I am very excited about it.

Soon, I hope to use PyDy to derive the model of the bicycle for our studies. Once trigsimp() is made better, it will be fully capable of this :)

Tags: , , ,
Posted in Uncategorized | 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 »

PyDy Google group website and emailing list

Friday, May 22nd, 2009

The Google group page is here:

http://groups.google.com/group/pydy

The emailing list is here:

pydy@googlegroups.com

Tags:
Posted in Uncategorized | No Comments »