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

# 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. = FRAMmc = FORKmc = WFmc = 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])

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,

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: 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
            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)

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

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

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)

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.