Google Summer of Code, Mentor Summit 2012

I probably should have written a blog post during the summer, displaying Angadh’s work . . . oh well.

SymPy sent me as one of 3 mentors to the GSoC Mentor Summit this year, at the Googleplex. I was very excited to go, meet some of the other SymPy developers (Matthew and Stefan), and meet others in the open source community. My overall experience was a little mixed though.

Matthew covered some feelings/conclusions I also share in his write-up on the event, but I have some others.

Meeting and interacting with people in the open source community was probably the highlight of the trip for me. It was interesting to talk to a few people who had taken projects that were originally academic/research code and translated them into more successful (open source) products. I’m not sure if it will be relevant to SymPy’s future, but it might be to mine. A common thread seemed to be: jumpstart a project with grants, keep improving it for long enough with that money, and then put it in a position where it can join an umbrella organization or can be used by professionals who will pay for support. I certainly don’t speak for any other SymPy developers, but I’ve never got the feeling that this was the intended trajectory of the project (at least the selling support part).

There was a talk on how to structure student/mentor/organization interactions for future GSoC projects. Some of the organizations have a much more defined structure than SymPy though, and benefit from things like daily group meetings. SymPy seems to have a more distributed organizational structure though – there are a lot of different modules, that all have some independence from each other. Despite this, SymPy’s code base is of very high quality, with credit going to the review-process/reviewers and the high standards that are enforced.

The talk on forming non-profit organizations was also interesting. One major takeaway was that the IRS would rather you put your project under an existing umbrella, rather than grant you 501(c) status, due to the potential for abuse. Also, a lot of work is involved in managing money properly once a certain amount of cash flowing through. Although there are hurdles, getting to form a board of directors sounded interesting. Scheduling board meetings could also be fun (with the money being spent responsibly, of course…).

On the less positive side of things: the unconference format. I think that it could have worked a lot better. Again, Matthew touched on this in his post, but all of the sessions (that were not just presentations) were very unfocused, with some more productive than others. The lack of moderation was a serious impediment to keeping the discussions on track. There were definitely a few times were people would hijack a session to try and talk about or show off their project. While I don’t have a problem with people showing off their work, there was limited time for each session.

The time limitations on the sessions was another issue. Each one was ~45 minutes, and almost always there was another session waiting for you to leave the room so they could start. I think in a few of the talks I went to, people would have been happy to stay in the room and continue to discuss the subject. Perhaps the expectations, of the session initiators (and myself), on what we would accomplish were too high. Perhaps the 45 minutes should have been spent networking with other people thinking about that topic, and then spawning a more detailed discussion in the future? Again, I think moderation would have helped this.

There were also some talks, with interesting sounding titles, that were just unproductive. There was too much recounting of what one organization did, with little generalization to what others could use. There was also little consideration of what decisions were made which led to important decisions; e.g. the fact that a project was split into parts A and B which were developed separately was recounted, but not what went into making that decision. Perhaps another example of expectations being too high…

I’m certainly not going to write off the unconference format – I think it could have led to some really cool things. But I don’t think I will attend another one that has sessions which are so short and under-moderated.

GSoC, Week 12

Well, this past week was the final “full” week of the 2011 Google Summer of Code. The ‘soft’ pencils down date is Monday.

This week, I mostly rewrote interfaces to some functions and classes. I should probably go back through previous blog posts and update those to show how they are different, so if anyone comes across this in the future, there won’t be example code which doesn’t work. I also tried to finish up the documentation this week. I decided a little late in the week to add another page, a “advanced/future interfaces” page, to describe how some functions/classes have different interfaces (and their advantages/disadvantages) as well as some discussed extensions to the current interface (that probably won’t be implemented right now). Some of these things were: more ways to access basis vectors, more printing options, and more ways to initialize the Kane object. Like I said, I don’t think I’ll be able to implement these now. I think I’ll try and do those three things this Fall though.

I also went through an example of bringing non-contributing forces into evidence (I discussed this topic a little in last week’s blog post). With the changes to the Kane interface it has ended up working out quite nicely. I put some code for this into the examples section, so hopefully it will be relatively easy for others to figure out. I didn’t make as much progress as I’d like with the code output function, but I think it can do an adequate job at the moment. I also spent most of a day rebasing all of my code of the latest master. I’m not sure what I did wrong, but somehow, I managed to get double commits again, so I took care of this. I think it might have been how I branched, rebased on one branch, then merged? It took some time to sort out.

I see there being two more things to do. The first is updating the images I have in the documentation; they’re all hand-drawn and of low quality. I plan on spending a little time this weekend on that, as well as Monday. The other thing is updating my pull request, getting others to look at it, and making necessary changes (and hopefully during this time, I’ll be able to check out other people’s pull requests and help them). Working with the pull request will probably take longer than updating the images…

Anyways, this is the last required blog post, but I think I’ll do at least one more as part of the Summer of Code, to wrap things up, and maybe introduce one more example.

GSoC, Week 11

Well, it looks like there is a little over 1 full week left. This week’s post probably won’t be too long. I have a list of the last few things left to work on. The last few (big) pieces are: finishing code output, making sure you can bring non-contributing forces into evidence, and documentation work.

Finishing the code output is fairly straight forward; there needs to be some thought given to: output quantities, time to write to file, correct use of cse(), and C code output. I already have code which writes to a file MATLAB or SciPy code for numerical integration, with an option to use cse() on the expressions. Using cse() makes it take longer to generate the code, but the code can execute significantly faster (especially with the bike, where there are a few hundred “common sub-expressions”). Writing to a file unfortunately seems to cost a similar amount of time as printing to the screen; in the case of the bike, a few minutes. There also needs to be some though of how to handle output quantities. When numerically integrating equations of motion, it’s common to look at things like total energy, or a body’s angular momentum, etc. A nice, clean way to work these into the integration would be ideal. There is also the issue of dealing with things like matrices as part of the code output; one might want to output a matrix for doing animations. Deciding on how to accept a matrix and format it for output also needs to be done.

One of the advantages of Kane’s Method is that “non-contributing forces” do not come into the equations of motion. Non-contributing forces are things like normal forces, contact (but not friction forces) such as pinned or rolling connections, along with some others. Unfortunately, this can also be a disadvantage of Kane’s Method; frequently one wants to know the value of these forces, in cases like: checking that normal forces do not go negative, calculating friction forces, joint constraint forces, etc. There is a way of using Kane’s Method which involves introducing “auxiliary speeds” and “fictitious forces”; the auxiliary speeds are defined as zero, additional equations are introduced, and then the forces are solved for (in terms of the “real” speeds).  I need to go through this to make sure my code is compatible with it; additionally, I’ve never done this for a system with non-holonomic constraints though, so I’ll have to check this out too.

Finally I’m finishing up the documentation. I have to go through and check docstring and doctest coverage; for some functions I don’t know if a doctest makes sense, such as generating a file from the code output, and I’m not sure what to do here, so I’ll have to look into this. Also, I have to finish the sphinx-documentation; most of it is done, but I need to go through and make sure I’ve covered some of the newer things I’ve added.

I hope by the end of next week to sort out my branches; right now they’re a little messy; I want to make them a little more organized. I’m planning on going back to having just a pydy branch, and a pydy pull branch. I’ll just work with the pydy-pull brach this time. Previously I made that branch, continued development on pydy, but also made corrections to pydy-pul in request to people’s comments. With development “stopping”, I think that will make things a bit less confusing by just working with issues on the pydy-pull branch.

GSoC, Week 10

Week 10….?

This week, there really isn’t that much to report. Week 10 was supposed to be finishing up the pretty and LaTeX printer, and I have done most of that. Tests and documentation need to be written for both though; I think I’ve learned that things haven’t really been coded successfully until you’ve done those two steps.

I also got some more work done with the documentation. It’s a little further behind than I had hoped for, but it’s coming along. I definitely need to redo the images I have throughout it so far; they are just scans of sketches I have done, and are of low-quality. There also needs to be a lot more talk about linearization in the documentation. That in fact leads to the last two big coding things.

I think the last two major things are: more work in linearization and code output for numeric integration. The linearization routine still needs some work; mainly in dealing with situations where the qdot’s are defined with coefficients that depend on the q’s. The question of what to do with user defined dynamic symbols (such as forces or specified position) is also still up in the air. Forces will probably be easy to linearize (or really, take the partial with respect to), but with a specified position, it becomes more complicated; during the process of forming the equations of motion, most likely the derivative of that specified position will be brought into the expressions. I haven’t really decided what to do in this situation, as clearly some distance (say, l), and its derivative (say, l’) are not independent. I’m not sure if adding the position as a system state is the right thing to do here. I think I’ll play around with some simple examples to see what makes the most sense, and consult some other people. The part about code output also relates here, in that I need to decide how to deal with this situation (a value and its derivative). I can imagine just putting both as empty for the user to fill in, and telling them to do it the right way. Hopefully I’ll have this figured out more by next week. I also do need to write the code output stuff. I don’t think it will be that hard, I just have to make the decisions about the formats of of the output. The size of the expressions generated can be problematic though, taking tens of minutes to print out. I’m also not sure how to deal with this situation. Hopefully I can again find out more information on this issue by next weekend.

 

GSoC, Week 9

Week 9 already….? Also, my internet connection is down, so this post is a little later than I planned….

So there was some progress this week. A lot of it was in doing the math for the linearization process. Linearization is normally easy, but when dealing with dependent quantities becomes more complicated. It requires treating the dependent quantities as functions of the independent quantities, then taking partial derivatives. It ends up being simplest (that is the current theory at least) when using the chain rule; this was done by hand in order to ease the work done by SymPy’s routines. These expressions can get quite large, so minimizing the buildup of expression size is important. Or, to put it another way, we figured out the best way to operate on smaller chunks of expressions, rather than larger ones, and then combine them afterwards, minimizing the size of some of the expressions which have to be dealt with.

In our lab at UC Davis, one of the big projects currently is studying bicycle dynamics. A paper was written in 2007 which sought to provide benchmark/reference values to validate models of the bicycle against: J.P Meijaard, Jim M Papadopoulos, Andy Ruina and A.L Schwab. Linearized dynamics equations for the balance and steer of a bicycle: a benchmark and review. Proceedings of The Royal Society (2007) 463, 1955-1982 doi: 10.1098/rspa.2007.1857 . One of the goals for PyDy this year was to be able to generate the full nonlinear equations for the bicycle, and the linear equations for comparison of eigenvalues using the reference values. Well, it looks like my code can do that successfully now:

The code at the bottom of the post calculates the “A” matrix, and writes a few eigenvalues out. You can compare them to those presented in the paper. There are still 2 issues left to investigate as part of the linearization process, there is pretty printing, latex printing, and finally code output and all of the functionality will be “done”. There is also still lots of documentation to write…

One other issue, that is proving to be somewhat serious, is the speed of the code. The execution time of the code is usually not too terrible, and if it does get bad, turning of Vector’s auto simplification (Vector.simp = False) usually helps. The linearization can also be kinda slow, but not overwhelmingly so. What seems to really be an issue is both printing (to screen or file) and substituting. The expressions get fairly big (I think the forcing term for the bicycle is on the order of 1MB when in Ascii form) so I’m not completely unhappy or surprised about the printing issue. I am more unhappy about the substitution being slow. Right now I’m using .subs(), and it seems to take a few minutes. I’ve also tried .evalf(subs=dict) and lambdify, but .evalf hits a maximum recursion depth error, and I’m getting a syntax error with lambdify (it actually looks like this is due to the use of unknown functions; I’m sure I can solve this somehow).

Use the pydy-pull-funcderiv branch; that’s the current branch. Some pushes have been forced in, so beware….sorry. Also, Brian Granger’s code got pulled in, so I’ll be defining ‘dynamic symbols’ as undefined functions of time. My branches need a little work to get them organized and unified again.

Code Output:

Calculation of Linearized Bicycle "A" Matrix, with States: Roll, Steer, Roll Rate, Steer Rate
Before Forming the List of Nonholonomic Constraints.
Before Handling of Dependent Speeds
Before Forming Generalized Active Forces, Fr
Before Forming Generalized Inertia Forces, Fr*
Base Equations of Motion Computed
Before Linearization of the "Forcing" Term
Before Substitution of Numerical Values
[               0,                                           0,                  1.0,                    0]
[               0,                                           0,                    0,                  1.0]
[9.48977444677355, -0.891197738059088*v**2 - 0.571523173729246, -0.105522449805691*v, -0.330515398992311*v]
[11.7194768719633,   -1.97171508499972*v**2 + 30.9087533932407,   3.67680523332153*v,  -3.08486552743311*v]
v = 1
{-3.13423125066578: 1, 3.52696170990069 - 0.80774027519931*I: 1, 3.52696170990069 + 0.80774027519931*I: 1, -7.11008014637441: 1}
v = 2
{-8.67387984831737: 1, -3.07158645641514: 1, 2.68234517512745 + 1.68066296590676*I: 1, 2.68234517512745 - 1.68066296590676*I: 1}
v = 3
{-2.63366137253665: 1, 1.70675605663973 + 2.31582447384324*I: 1, 1.70675605663973 - 2.31582447384324*I: 1, -10.3510146724592: 1}
v = 4
{0.41325331521124 - 3.07910818603205*I: 1, -12.1586142657644: 1, 0.41325331521124 + 3.07910818603205*I: 1, -1.42944427361326: 1}
v = 5
{-14.0783896927982: 1, -0.322866429004087: 1, -0.775341882195845 + 4.46486771378823*I: 1, -0.775341882195845 - 4.46486771378823*I: 1}

Actual Code:

from sympy import *
from sympy.physics.mechanics import *

# Code to get equations of motion for a bicycle modeled as in:
# J.P Meijaard, Jim M Papadopoulos, Andy Ruina and A.L Schwab. Linearized
# dynamics equations for the balance and steer of a bicycle: a benchmark and
# review. Proceedings of The Royal Society (2007) 463, 1955-1982
# doi: 10.1098/rspa.2007.1857

print('Calculation of Linearized Bicycle \"A\" Matrix, with States: Roll, '
      'Steer, Roll Rate, Steer Rate')

# Note that this code has been crudely ported from Autolev, which is the reason
# for some of the unusual naming conventions. It was purposefully as similar as
# possible in order to aide debugging.

# Vector's simplification routines need to be turned off, otherwise the
# expressions get too big and too slow to simplify
Vector.simp = False
mechanics_printing()

# Declare Coordinates & Speeds
q1, q2, q4, q5 = dynamicsymbols('q1 q2 q4 q5')
q1d, q2d, q4d, q5d = dynamicsymbols('q1 q2 q4 q5', 1)
u1, u2, u3, u4, u5, u6 = dynamicsymbols('u1 u2 u3 u4 u5 u6')
u1d, u2d, u3d, u4d, u5d, u6d = dynamicsymbols('u1 u2 u3 u4 u5 u6', 1)

# Declare System's Parameters
WFrad, WRrad, htangle, forkoffset = symbols('WFrad WRrad htangle forkoffset')
forklength, framlength, forkcg1 = symbols('forklength framlength forkcg1')
forkcg3, framcg1, framcg3, Iwr11 = symbols('forkcg3 framcg1 framcg3 Iwr11')
Iwr22, Iwf11, Iwf22, Ifram11 = symbols('Iwr22 Iwf11 Iwf22 Ifram11')
Ifram22, Ifram33, Ifram31, Ifork11 = symbols('Ifram22 Ifram33 Ifram31 Ifork11')
Ifork22, Ifork33, Ifork31, g = symbols('Ifork22 Ifork33 Ifork31 g')
mframe, mfork, mwf, mwr = symbols('mframe mfork mwf mwr')

# Set up reference frames for the system
# N - inertial
# Y - yaw
# R - roll
# WR - rear wheel, rotation angle is ignorable coordinate so not oriented
# FRAM - bicycle frame
# TEMPFRAM - statically rotated frame for easier reference inertia definition
# FORK - bicycle fork
# TEMPFORK - statically rotated frame for easier reference inertia definition
# WF - front wheel, again posses a ignorable coordinate
N = ReferenceFrame('N')
Y = N.orientnew('Y', 'Simple', q1, 3)
R = Y.orientnew('R', 'Simple', q2, 1)
FRAM = R.orientnew('FRAM', 'Simple', q4 + htangle, 2)
WR = ReferenceFrame('WR')
TEMPFRAM = FRAM.orientnew('TEMPFRAM', 'Simple', -htangle, 2)
FORK = FRAM.orientnew('FORK', 'Simple', q5, 1)
TEMPFORK = FORK.orientnew('TEMPFORK', 'Simple', -htangle, 2)
WF = ReferenceFrame('WF')

# Declaration of the RigidBody containers
BodyFram = RigidBody()
BodyFork = RigidBody()
BodyWR = RigidBody()
BodyWF = RigidBody()

# Setting the masses for the bodies
BodyFram.mass = mframe
BodyFork.mass = mfork
BodyWF.mass = mwf
BodyWR.mass = mwr

# Assigning the appropriate frames to each body
BodyFram.frame = FRAM
BodyFork.frame = FORK
BodyWR.frame = WR
BodyWF.frame = WF

# Kinematics of the Bicycle
# First block of code is forming the positions of the relevant points
# rear wheel contact -> rear wheel mass center -> frame mass center +
# frame/fork connection -> fork mass center + front wheel mass center -> front
# wheel contact point
WRhat = Point('WRhat')
WRmc = WRhat.newpoint('WRmc', WRrad * R.z)
STEER = WRmc.newpoint('STEER', framlength * FRAM.z)
FRAMmc = WRmc.newpoint('FRAMmc', -framcg1 * FRAM.x + framcg3 * FRAM.z)
FORKmc = STEER.newpoint('FORKmc', -forkcg1 * FORK.x + forkcg3 * FORK.z)
WFmc = STEER.newpoint('WFmc', forklength * FORK.x + forkoffset * FORK.z)
WFhat = WFmc.newpoint('WFhat', WFrad*(dot(FORK.y, Y.z)*FORK.y - Y.z).unit)

# Set the angular velocity of each frame.
# Angular accelerations end up being calculated automatically by
# differentiating the angular velocities when first needed.
# u1 is yaw rate
# u2 is roll rate
# u3 is rear wheel rate
# u4 is frame pitch rate
# u5 is fork steer rate
# u6 is front wheel rate
Y.set_ang_vel(N, u1 * Y.z)
R.set_ang_vel(Y, u2 * R.x)
WR.set_ang_vel(FRAM, u3 * FRAM.y)
FRAM.set_ang_vel(R, u4 * FRAM.y)
FORK.set_ang_vel(FRAM, u5 * FORK.x)
WF.set_ang_vel(FORK, u6 * FORK.y)

# Form the velocities of the previously defined points, using the 2 - point
# theorem (written out by hand here).
# Accelerations again are calculated automatically when first needed.
WRhat.set_vel(N, 0)
WRmc.set_vel(N, WRhat.vel(N) + (WR.ang_vel_in(N) ^ WRmc.pos_from(WRhat)))
STEER.set_vel(N, WRmc.vel(N) + (FRAM.ang_vel_in(N) ^ STEER.pos_from(WRmc)))
FRAMmc.set_vel(N, WRmc.vel(N) + (FRAM.ang_vel_in(N) ^ FRAMmc.pos_from(WRmc)))
FORKmc.set_vel(N, STEER.vel(N) + (FORK.ang_vel_in(N) ^ FORKmc.pos_from(STEER)))
WFmc.set_vel(N, STEER.vel(N) + (FORK.ang_vel_in(N) ^ WFmc.pos_from(STEER)))
WFhat.set_vel(N, WFmc.vel(N) + (WF.ang_vel_in(N) ^ WFhat.pos_from(WFmc)))

# Assign the relevant points to each body.
BodyFram.mc = FRAMmc
BodyFork.mc = FORKmc
BodyWF.mc = WFmc
BodyWR.mc = WRmc

# Sets the inertias of each body. Uses the inertia frame to construct the
# inertia dyadics. Wheel inertias are only defined by principle moments of
# inertia, and are in fact constant in the frame and fork reference frames; it
# is for this reason that the orientations of the wheels does not need to be
# defined. The frame and fork inertias are defined in the 'TEMP' frames which
# are fixed to the appropriate body frames; this is to allow easier input of
# the reference values of the benchmark paper. Note that due to slightly
# different orientations, the products of inertia need to have their signs
# flipped; this is done later when entering the numerical value.
BodyFram.inertia = (inertia(TEMPFRAM, Ifram11, Ifram22, Ifram33, 0, 0, Ifram31), FRAMmc)
BodyFork.inertia = (inertia(TEMPFORK, Ifork11, Ifork22, Ifork33, 0, 0, Ifork31), FORKmc)
BodyWR.inertia = (inertia(FRAM, Iwr11, Iwr22, Iwr11), WRmc)
BodyWF.inertia = (inertia(FORK, Iwf11, Iwf22, Iwf11), WFmc)

print 'Before Forming the List of Nonholonomic Constraints.'

# The kinematic differential equations; they are defined quite simply. Each
# entry in this list is equal to zero.
kd = [q1d - u1, q2d - u2, q4d - u4, q5d - u5]

# The nonholonomic constraints are the velocity of the front wheel contact
# point dotted into the X, Y, and Z directions; the yaw frame is used as it is
# "closer" to the front wheel (1 less DCM connecting them). These constraints
# force the velocity of the front wheel contact point to be 0 in the inertial
# frame; the X and Y direction constraints enforce a "no-slip" condition, and
# the Z direction constraint forces the front wheel contact point to not move
# away from the ground frame, essentially replicating the holonomic constraint
# which does not allow the frame pitch to change in an invalid fashion.
conlspeed = [WFhat.vel(N) & Y.x, WFhat.vel(N) & Y.y, WFhat.vel(N) & Y.z]
# The holonomic constraint is that the position from the rear wheel contact
# point to the front wheel contact point when dotted into the normal-to-ground
# plane direction must be zero; effectively that the front and rear wheel
# contact points are always touching the ground plane. This is actually not
# part of the dynamic equations, but instead is necessary for the lineraization
# process.
conlcoord = [WFhat.pos_from(WRhat) & Y.z]

# The force list; each body has the appropriate gravitational force applied
# at its mass center.
FL = [(FRAMmc, -mframe * g * Y.z), (FORKmc, -mfork * g * Y.z), (WFmc,
      -mwf * g * Y.z), (WRmc, -mwr * g * Y.z)]
BL = [BodyFram, BodyFork, BodyWR, BodyWF]

# The N frame is the inertial frame, coordinates are supplied in the order of
# independent, dependent coordinates, as are the speeds. The kinematic
# differential equation are also entered here.
KM = Kane(N)
KM.coords([q1, q2, q5, q4])
KM.speeds([u2, u3, u5, u1, u4, u6])
KM.kindiffeq(kd)

print 'Before Handling of Dependent Speeds'

# Here the dependent speeds are specified, in the same order they were provided
# in earlier, along with the non-holonomic constraints.
# The dependent coordinate is also provided, with the holonomic constraint.
# Again, this is only provided for the linearization process.
KM.dependent_speeds([u1, u4, u6], conlspeed)
KM.dependent_coords([q4], conlcoord)

print 'Before Forming Generalized Active Forces, Fr'
fr = KM.form_fr(FL)
print 'Before Forming Generalized Inertia Forces, Fr*'
frstar = KM.form_frstar(BL)
print 'Base Equations of Motion Computed'

# This is the start of entering in the numerical values from the benchmark
# paper to validate the eigen values of the linearized equations from this
# model to the reference eigen values. Look at the aforementioned paper for
# more information. Some of these are intermediate values, used to transform
# values from the paper into the coordinate systems used in this model.
PaperRadRear                    =  0.3
PaperRadFront                   =  0.35
HTA                             =  evalf.N(pi/2-pi/10)
TrailPaper                      =  0.08
rake                            =  evalf.N(-(TrailPaper*sin(HTA)-(PaperRadFront*cos(HTA))))
PaperWb                         =  1.02
PaperFramCgX                    =  0.3
PaperFramCgZ                 	=  0.9
PaperForkCgX                  	=  0.9
PaperForkCgZ                	=  0.7
FramLength                  	=  evalf.N(PaperWb*sin(HTA)-(rake-(PaperRadFront-PaperRadRear)*cos(HTA)))
FramCGNorm                     	=  evalf.N((PaperFramCgZ-PaperRadRear-(PaperFramCgX/sin(HTA))*cos(HTA))*sin(HTA))
FramCGPar                      	=  evalf.N((PaperFramCgX/sin(HTA) + (PaperFramCgZ-PaperRadRear-PaperFramCgX/sin(HTA)*cos(HTA))*cos(HTA)))
tempa                         	=  evalf.N((PaperForkCgZ - PaperRadFront))
tempb                         	=  evalf.N((PaperWb-PaperForkCgX))
tempc                         	=  evalf.N(sqrt(tempa**2+tempb**2))
PaperForkL                      =  evalf.N((PaperWb*cos(HTA)-(PaperRadFront-PaperRadRear)*sin(HTA)))
ForkCGNorm                  	=  evalf.N(rake+(tempc * sin(pi/2-HTA-acos(tempa/tempc))))
ForkCGPar              		=  evalf.N(tempc * cos((pi/2-HTA)-acos(tempa/tempc))-PaperForkL)

# Here is the final assembly of the numerical values. The symbol 'v' is the
# forward speed of the bicycle (a concept which only makes sense in the
# upright, static equilibrium case?). These are in a dictionary which will
# later be substituted in. Again the sign on the *product* of inertia values is
# flipped here, due to different orientations of coordinate systems.
v = Symbol('v')
val_dict = {WFrad: PaperRadFront,
            WRrad: PaperRadRear,
            htangle: HTA,
            forkoffset: rake,
            forklength: PaperForkL,
            framlength: FramLength,
            forkcg1: ForkCGPar,
            forkcg3: ForkCGNorm,
            framcg1: FramCGNorm,
            framcg3: FramCGPar,
            Iwr11: 0.0603,
            Iwr22: 0.12,
            Iwf11: 0.1405,
            Iwf22: 0.28,
            Ifork11: 0.05892,
            Ifork22: 0.06,
            Ifork33: 0.00708,
            Ifork31: 0.00756,
            Ifram11: 9.2,
            Ifram22: 11,
            Ifram33: 2.8,
            Ifram31: -2.4,
            mfork: 4,
            mframe: 85,
            mwf: 3,
            mwr: 2,
            g: 9.81,
            q1: 0,
            q2: 0,
            q4: 0,
            q5: 0,
            u1: 0,
            u2: 0,
            u3: v/PaperRadRear,
            u4: 0,
            u5: 0,
            u6: v/PaperRadFront}

# Here a dictionary is formed using the kinematic differential equations. The
# expression is perhaps slightly more complicated then necessary in this case,
# but should work in all cases in order the generate a dictionary in the form
# {qd: f(u)}.
sub_dict = solve_linear_system_LU(Matrix([KM._k_kqdot.T,
        -(KM._k_ku*Matrix(KM._u) + KM._f_k).T]).T, KM._qdot)

print 'Before Linearization of the \"Forcing\" Term'
# Linearizes the forcing vector; the equations are set up as MM udot = forcing,
# where MM is the mass matrix, udot is the vector representing the time
# derivatives of the generalized speeds, and forcing is a vector which contains
# both external forcing terms and internal forcing terms, such as centripital
# or coriolis forces.
# This actually returns a matrix with as many rows as *total* coordinates and
# speeds, but only as many columns as independent coordinates and speeds.
forcing_lin = KM.linearize().subs(sub_dict)

# As mentioned above, the size of the linearized forcing terms is expanded to
# include both q's and u's, so the mass matrix must have this done as well.
# This will likely be changed to be part of the linearized process, for future
# reference.
MM_full = (KM._k_kqdot).row_join(zeros((4, 6))).col_join((zeros((6,
                4))).row_join(KM.mass_matrix))

print 'Before Substitution of Numerical Values'
# I think this is pretty self explanatory. It takes a really long time though.
# I've experimented with using evalf with substitution, this failed due to
# maximum recursion depth being exceeded; I also tried lambdifying this, and am
# not sure what the error message I got there meant.
MM_full = MM_full.subs(val_dict).evalf()
forcing_lin = forcing_lin.subs(val_dict).evalf()

# Finally, we construct an "A" matrix for the form xdot = A x (x being the
# state vector), although in this case, the sizes are a little off. The
# following line extracts only the minimum entries required for eigenvalue
# analysis, which correspond to rows and columns for lean, steer, lean rate,
# and steer rate.
Amat = MM_full.inv() * forcing_lin
A = Amat.extract([1,2,4,6],[1,2,3,5])
print A
print 'v = 1'
print A.subs(v, 1).eigenvals()
print 'v = 2'
print A.subs(v, 2).eigenvals()
print 'v = 3'
print A.subs(v, 3).eigenvals()
print 'v = 4'
print A.subs(v, 4).eigenvals()
print 'v = 5'
print A.subs(v, 5).eigenvals()

GSoC, Week 8

I don’t have as much to talk about this week, as I spent two days in travel. Mostly, I worked on the sphinx documentation for mechanics. I am continuing with the previous organization of the documentation: first the mathematical descriptions, then the SymPy implementations. I’m still not convinced that this is the best way to do things, but at the moment, at least the subsubmodule’s documentation is consistent. Also, there will most likely be changes to the code, probably relating to the switch to using Brian Granger’s functional derivative branch and the different printing. I’m also trying to use the same examples in both the mathematical sections and in the code sections. I think this is the most logical approach, and if I choose to rearrange the documentation, I won’t have to rewrite much. This is all in my ‘pydy’ branch.

The other thing I am working on is forming the equations of motion for the bicycle. Right now they are being computed, but I’m not sure of the accuracy. There are already tests in place for holonomic systems, but non for non-holonomic systems (for reference: http://en.wikipedia.org/wiki/Nonholonomic_system). Unfortunately, the equations for non-holonomic systems get very, very big. In our lab, with the bicycle, we form our equations, linearize them, substitute in the system parameters, and examine the eigenvalues in order to validate our equations. I’m currently using my SymPy code to do this with the bicycle model, but am encountering some difficulties. I’ve got the mass matrix and forcing terms solved for and linearized, and the udots have a solution, in the form udot = A x, where x = [q2,q4,q5,u2,u3,u5,u1,u4,u6].T . What should happen next is that the relevant entries from this matrix are extracted, as our final state is [u2,u5,q2,q5] (this is lean rate, steer rate, lean angle, steer angle), then a smaller matrix is constructed and the eigenvalues are examined. Now, I know that the eigenvalues are off already. Examining this “A” matrix though shows that some of the elements are correct, compared to a reference “A” matrix. It would appear that the partial derivative of the udots w.r.t. the lean angle and lean rate are correct, but are not right w.r.t. the steer angle and steer rate. I’m not sure how to interpret this, as it is showing that half of the correct numbers are there (and these are correct to ~13 significant digits when evaluated numerically), and the other half are wrong. I’m hoping that I have just performed the linearization wrong, or the order in which I have arranged my matrices is incorrect; I don’t see many other options.

The full equations come out in the form: MM udot = forcing, or the mass matrix multiplied by the time derivatives of the generalized speeds equals the forcing terms. MM is 6×6, and forcing is 6×1 (as is udot). What is done next is substitute in the current state into the mass matrix, find the jacobian of the forcing vector w.r.t. the vector x (from above), and substitute in the numerical values. Next, udot is found by performing MM.inv() * forcing (where forcing is now a 6×9 matrix). Since I am getting some of the correct values, it would appear that there are elements in both the MM and the forcing matrix which are correct. I’m going to continue to play around with things, as again, I feel that I am most likely assembling these things in the incorrect order. Hopefully this will prove successful though, and will provide a good, non-trivial example of what physics.mechanics can do.

GSoC, Week 7

I submitted my pull request a little while ago for my code. Currently, it looks like I’m waiting on a few things. I think the biggest thing is the pull request Brian Granger has opened. I’ve already written a branch which uses his code, and everything works. I just have to wait for his pull request to get accepted in order to switch what is in my pull request. Then of course, that has to get approved. This delay has had some positive benefits though.

Working with other people’s branches has definitely helped my familiarity and confidence with using git. Using the fetch and merge commands with multiple remotes, along with getting merge errors and fixing them, has given some me good experience. Also, rebasing to play with commit history, and amending commits has been helpful.

I feel like I had the functionality I wanted before from my ‘DynamicSymbol’ class, but with the ability to take the derivatives of functions now, with respect to functions means that this object is no longer needed really. What I did before was have ‘DynamicSymbol’ extend ‘Symbol’. I then used the ability of the ‘symbols’ function to to supply my own class for creation in the ‘dynamicsymbols’:

def dynamicsymbols(names):
    """Wraps sympy.symbols to use DynamicSymbol. """
    return symbols(names, cls=DynamicSymbol)

This was an OK solution. Now though, one can use functions to represent these time varying quantities; instead of q1, q1d, q1dd we’ll have q1(t), Derivative(q1(t), t), etc. This is certainly more consistent with the rest of SymPy. Now my ‘dynamicsymbols’ function uses symbols with ‘cls=Function’ to create a number of undefined function objects, which then are called with ‘(t)’ as an argument (I actually do something like dynamicsymbols._t = Symbol(‘t’). This way the time value is specified in only one place, but that place is still associated with dynamic symbols. The ability of Python to set an attribute to a function is certainly interesting…). I also added the ability to specify the level of differentiation; you can call:

In [11]: dynamicsymbols('q1 q2')
(q1, q2)

In [12]: dynamicsymbols('q1 q2', 1)
(q1', q2')

This means you can’t create q1, q1′, q1” on one line, but at least all dynamic symbols of the same level of differentiation can be called together. The above snippet leads to the next topic though: printing.

When you start up ipython, and do the following:

In [1]: from sympy.physics.mechanics import *

In [2]: dynamicsymbols('q1 q2', 1)
Out[2]: (Derivative(q1(t), t), Derivative(q2(t), t))

the output isn’t actually the same as what I presented a few lines above. This is because of the default printing of Derivative (shown in the line immediately above). When dealing with multibody dynamics problems, you get a lot of derivatives of these dynamic symbols (generalized coordinates and speeds). It gets out of control pretty fast if we print out every time derivative in the above form. So, what I did was to write my own printer which shows derivatives as q’ instead of Derivative(q(t), t), under certain conditions (and q” for Derivative(q(t),t,t) and so forth). The conditions are that that the derivative has to have been taken with respect to the time symbol (stored in dynamicsymbols._t, as previously mentioned) and that its first argument (the value which is being differentiated) has to be an ‘UndefinedFunction’ and that the derivative can only be taken with respect to t.

    def _print_Derivative(self, e):
        from sympy.core.function import UndefinedFunction
        t = dynamicsymbols._t
        if (bool([i == t for i in e.variables]) &
            isinstance(type(e.args[0]), UndefinedFunction)):
            ol = str(e.args[0].func)
            for i, v in enumerate(e.variables):
                ol += '\''
            return ol
        else:
            return StrPrinter().doprint(e)

The above block of code shows how it is done; the first part of the undefined function function is added to our output string (e.args[0].func), then for each entry in e.variables we print a ‘. In this method, ‘e’ is the expression to be printed. Then our final list is returned, or just the normal printer’s output if the given expression did not meet our criteria. Now, this actually doesn’t print out everything perfectly:

In [23]: Derivative(q1,t)
q1'

In [24]: Derivative(q1*x,t)
Derivative(x*q1(t), t)

In [25]: Derivative(q1*x,t, evaluate=True)
x*q1'

We can still get examples where it prints out the long way: when evaluate is not true and the ‘Derivative’ object does not only contain an undefined function. I’m willing to live with this at the moment though, as I think that all my code uses ‘diff’ where evaluate always gets called. I also wrote a ‘_print_Function’ method, in order to print q1(t) as q1:

    def _print_Function(self, e):
        from sympy.core.function import UndefinedFunction
        t = dynamicsymbols._t
        temp = StrPrinter().doprint(e)
        if isinstance(type(e), UndefinedFunction):
            return temp.replace('(t)', '')
        return temp

I think this one is even simpler. If the Function is an undefined function, print it to a string, then replace the ‘(t)’ in the string with nothing. Checking that it is an undefined function is definitely important, otherwise you get:

In [10]: cos(t)
cos

And I don’t think anyone wants that. So the last two bits of this are interactive printing and integration of my non-SymPy objects with the printing system.

In order to get my Vector and and Dyad to print out with a SymPy printer correctly, I had to do the following: have one method (I chose __str__) which creates the string representation of the object, and then set __repr__ = __str__, _sympystr = __str__, _sympyrepr = __str__. Those last two functions are what the printer looks for, I believe. I’m actually still not completely clear on the difference between repr and str in Python, but I think one is supposed to satisfy: x == eval(repr(x)), except that I’m not sure if it is repr or str there. Additionally, in doctests, unless you do print, a returned value which you check against is called with repr I believe. SymPy does not exactly follow the above obj=eval(repr(obj)), and now, my code doesn’t always either, which I think is an acceptable tradeoff considering the readability issues. The way that the printer (which currently only has _print_Derivative and _print_Function methods overridden) is used by vector is as follows: the string value of a vector’s non-zero measure number is obtained with MechanicsStrPrinter().doprint(ar[i][0]). Here ar[i][0] is what I am supplying, but whatever is put in .doprint() will be printed by the printer.

Now for most of the work in physics.mechanics, you’ll be dealing with Vector quantities, so my printer is used automatically. But at the end of forming Kane’s equations, you are left with Fr + Fr* = 0. This is a r x 1 vector represented by a SymPy Matrix of the same shape. It will also be full of qdots and udots (generalized coordinate and speed time derivatives). This can end up being the biggest expression generated in forming the equations of motion, so being concise where it is practical is a good idea. Unfortunately my printer won’t get called here, because Matrix uses the default StrPrinter. Enter the display hook:

def mechanics_printing():
    def mydhook(ar):
        print MechanicsStrPrinter().doprint(ar)
    sys.displayhook = mydhook
    return 'displayhook set'

This function, mechanics_printing(), will be called by the user at the beginning of an interactive session in order to print out Derivatives and Functions which are not part of a mechanics Vector in the way I have chosen. You can set sys.displayhook to use a custom printing function of your creation. Here, a one line function to just use my printer is defined, then the system displayhook is set to use my displayhook. At this point I learned there is no such thing as a “void”  function in Python; everything has to return a value, and if you don’t, None is automatically returned. I just had it print that the displayhook was set, rather than it printing ‘None’ when you call this function. I believe this only sets the interactive printing though; I’m not sure what will be printed when you have a script file you have written and it prints out the screen or a file; this still needs to be worked out.

Hopefully all of this will help out someone in the future with their SymPy printing issues.

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 (http://ecommons.library.cornell.edu/handle/1813/638), 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()
    BodyD.mc = 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])
    KM.kindiffeq(kd)
    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: http://en.wikipedia.org/wiki/Dyadics, 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.kindiffeq(kd)
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.