Re: GSOC - Potential Student for Linearization Project

[Jason Moore]

Some comments:

- In general, you can only use the Jacobian if you have an unconstrained system and the operating point is zero. If the operating point is not zero, then you need the full linear portion of the Taylor series about some point other than zero. Also, for some subset of constrained systems and particular equilibria points the chain rule evaluates to the same thing as the Taylor series (this is the case for the upright equilibrium of the Whipple bicycle model). Jacobians are a matrix operation.

- If things are complicated enough, we typically use a class to manage data and then may write a functional wrapper around the class for a specific use case. Should your linearization code be structured as classes instead of a function with subroutines? This may warrant a "Linearizer" class or something if it will be easier to handle all of the data needed. I'd think some about that, maybe speak on pros and cons of each way.

- The PDE stuff isn't explained. I'd either remove it or expand. The PDE's will be a function of space in 3 dimensions and time for flexible systems or electromagnetic systems. We don't have those things working yet, so maybe this is thinking too far ahead. But if you designed your code to work with functions of space and time, it would work for whatever we eventually have when the PDE producing methods work.

- Maybe describe the mathematical form of the type of inputs to the linearize function. In general is it computationally inefficient to symbolically form xdot = f, we typically work with M * xdot = F as M^-1 * F should be computed numerically.

- I like the option to both output M, C, K and Q and A, B.

- Another thing we've found useful, is to linearize about an arbitrary equilibrium point, then generate numerical code from that. Then you have the equations that produce the linear system for any given equilibria. We should include some code generation functionality in the linearize class or update the current sympy code gen stuff to handle matrices easily.

- The early linearization has potential to make simplified linear forms of the linear EoMs. I'm not sure how this generalizes to all methods (Lagrange, Kane, Hamilton, etc), but Kane's book has a section detailing when exactly you can start to linearize the problem.

- This linearization stuff may or may not be enough to fill the summer. You seem very familiar with Python and you know dynamics, so you may be able to make quick headway. If you want to add optional items related to speeding up the symbolic derivations or adding another method (NewtonEulerMethod, HamiltonMethod, etc), then that may be a good idea.

- Some research into efficient methods of symbolic linearization could be done. You may be able to use the common structure and common terms in linearizing m-body systems to do some simplification along the way more efficiently.

- What happens when the equations of motion have funny things in them like Piecewise functions. Some non-linearities may be trickier to deal with and we should support them. Think about a bouncing ball simulation. How do we support those from the symbolic end? I suspect that an inverted pendulum with hard impact constraints at max angles could be tricky. In general, small angles would just remove that non-linearity, but we want to make sure these things are handled.

- Be sure to include plenty of time for creating documentation, tutorials, and examples. This always gets neglected! I wouldn't just wait till the end to do this, as it rarely gets done. Do it at each step of the process.

Jason
moorepants.info
+01 530-601-9791

On Wed, Mar 12, 2014 at 3:06 AM, James Crist <cris...@umn.edu> wrote:

All,

A draft of my application has been posted here. It could definitely use some critiquing and editing.

A couple things I wasn't sure of:

1.) The interface was based purely on my thoughts, and may not be desired by everyone

2.) Is there anything in the scope of this project that I missed? Is there anything I included that is not in scope?

3.) Mathematical background? Not enough? Too much? I'm used to Mathjax in my markdown, which is why there are $ littered everywhere. I'll clean those out later. Not sure if I can even submit math to melange.

4.) My background? I included it because sympy and pydy had it in their application template. However, I noticed that neither Sachin nor Tarun included it. Not sure how much is really necessary, as I assume most things should be project specific.


I will be leaving on Friday, and will (probably) not have access to internet until the following saturday, so I'm trying to get this done in the next few days. So comments would be highly appreciated.

Thanks!

-Jim

On Tuesday, March 11, 2014 10:53:21 PM UTC-5, James Crist wrote:

Jason,

Thanks for the clarification. With the added info, I can see how to make this into a full summer project. I'm working on finishing my proposal, and I should get a first draft up hopefully tonight. Here's my thoughts on the basic functionality:

1.) A function ("linearize", or something like it). Would work much like solve does, having multiple sub methods, and paths for different system types. Should be able to handle normal unconstrained linearization (just a jacobian), as well as constraints, and possibly pdes. Ideally this could linearize systems of equations generated in places other than mechanics as well. Call signature would be something along the lines of:

linearize(exression, *args, **kwargs)

This would allow for different types (just a state vector, state and input, etc...), as well as hints if facts about the system are known (nature of constraints, etc). The hints should not be necessary, but may reduce computation time.

2.) Methods added to lagrange, and a refactor or kane, that call linearize with hints (the classes should have some knowledge of the system, which would help)

3.) A "linearized_eom" method in kane and lagrange, that starts with linearization from the get-go. Should make deriving the eom faster, and makes sense if you know you're going to want to linearize them later on.

I envision "linearize" returning the M, K, C, and Q matrices for M*ddq + K*dq + C*q = Q. Could also add a kwarg for A and B matrices, if that standard form is more desirable. Those can easily be derived from M, K, C, and Q though.

Thoughts?

On Tuesday, March 11, 2014 10:31:30 PM UTC-5, moorepants wrote:

Another note on linearization. It is possible to linearize a system partially through the symbolic derivation. You can linearize before the formation of the equations of motion and get the linear system. It could be nice to include that in the Methods classes so that if you know you want a linear model you can linearize earlier on in the process.

Jason
moorepants.info
+01 530-601-9791

On Tue, Mar 11, 2014 at 3:45 PM, Jason Moore <moore...@gmail.com> wrote:

Jim,

Getting a really nice clean and robust linearizion class for multibody systems can potentially fill the the summer. As Aaron mentions, everyone overestimates what they can get done. Multibody systems which are simply a function of time (not space) have a specific form, essentially Mxdot = F. The paper shows how you can carefully linearize systems with complex constraints. The implementation of that would likely be straightforward, but we would want a class that can accept multiple "forms" of the EoM and then linearize the equations correctly and efficiently. You'd be developing example problems with no constraints and ones with complex constraints and using them as tests to ensure the linear equations are correct. The efficiency is important too. For the differentiation steps you can look into possible ways to improve the speed by utilizing the common form of the EoM's. There maybe some ways to speed up the differentiation from an algorithmic standpoint. Also this can be sped up by using the CSymPy core. We'd really like to get that working. Also from an academic standpoint, the linear EoMs should always be shorter and simpler than the nonlinear form, but because the algorithms for simplification are not great, we could spend time trying to ensure that our linearizer not only gets the correct equations but also simplifies them to a readable form. For many super complex nonlinear EoMs, there is often a human readable form of the linear EoMs about specific operating points. For example, if our linearizer could linearize the bicycle EoMs to it's simplest form, which can be done by deriving the linear model from the get go, then we'd have something more special than other CAS based multibody dynamic engines. So, simplification is a big challenge. Another thing would be to think about how this would play with a system that is a function of time and space (linearizing PDEs for flexible systems and the electromagnetic systems).

So the linearization project can be big, especially if you want to make it better than what currently exists on the market.

But speeding up the EoM derivations and the linearizations is important too. I think you will touch on these whatever you decide to do. We all want the code to be fast and there are lots of places to improve this.

Jason
moorepants.info
+01 530-601-9791

On Mon, Mar 10, 2014 at 8:57 PM, Aaron Meurer <asme...@gmail.com> wrote:

I don't know enough about the details to comment on this specific
project, but you should note that just about everyone underestimates
how long a given project will take. Whether they don't account for bug
fixing, debugging, documenting, or something else I don't know.

But having a little more to do at the end of the summer is not all bad
because that encourages you to stick around to finish it.

Aaron Meurer

On Mon, Mar 10, 2014 at 4:11 PM, James Crist <cris...@umn.edu> wrote:
> After having read through the paper Jason posted, I'm not sure if there's
> enough work to be done here over the summer. KanesMethod already has a
> linearization method implemented. If I understand correctly, the intent is
> that this method is pulled out and put into a class of it's own that will
> handle general linearization of equations of motion? Methods for handling
> lagrange equations would need to be added, and probably a general
> linearization routine for non-constrained systems (I think this would just
> do the jacobian). I'm just not sure how to break this up into meaningful
> steps to timeline for the proposal. If someone could explain the intentions
> here, that would be helpful.
>
> On that note, some of the other projects also look interesting. I took a
> look through simbody, and the functionality provided by their library could
> be very useful in association with pydy. Also, as can be seen in my
> lagranges method example posted yesterday, solving the equations of motion
> of a complicated system takes a considerable amount of time. Without subbing
> in constants beforehand, I let the calculations run for 45 min before I gave
> up and interrupted the operation. So improvements to equations of motion
> generation could be valuable as well.
>
> Basically, I'm looking for a project of substantial size. If the
> linearization project is more in depth than I'm thinking, I'd be fine
> submitting for that. Otherwise the Simbody wrapper or the improvements to
> E.O.M. generation efficiency sound more appealing. Could someone give me a
> brief overview of these projects and the intentions behind them?
>
> - Jim
>
>
>
>
> On Monday, March 3, 2014 4:30:02 PM UTC-6, moorepants wrote:
>>
>> Awesome, I'm looking forward to see what you have to show us both previous
>> work and a proposal.
>>
>>
>> Jason
>> moorepants.info
>> +01 530-601-9791
>>
>>
>> On Mon, Mar 3, 2014 at 4:56 PM, James Crist <cris...@umn.edu> wrote:
>>>
>>> Thanks for the info. I've actually already submitted a patch to the
>>> lambdify function, as you saw. When/if that gets merged I'll fix the related
>>> issue in the pydy codegen repo. I'll try to whip some of my past work into a
>>> nicely presented ipython notebook, and send a pull request to put it in the
>>> examples. Probably related to a robotic platform and the use of the
>>> lagrangian method. I'll also start a wiki page for my proposal. This is a
>>> busy week here, so most likely none of this will happen until around
>>> Thursday-Friday or so.
>>>
>>> Thanks,
>>> -Jim
>>>
>>>
>>> On Sunday, March 2, 2014 12:12:42 PM UTC-6, moorepants wrote:
>>>>
>>>> Some answers:
>>>>
>>>> - What can I do to get involved beforehand/even if I don't get accepted?
>>>> As I said, I make heavy use of sympy mechanics and the whole python science
>>>> stack in my research, and would love to contribute back.
>>>>
>>>> We have open issues in SymPy and the PyDy projects on Github. Browse
>>>> through them and see if you'd like to try fixing one. If you want to apply
>>>> to GSoC through SymPy then you will have to submit a pull request there.
>>>> Best to get started on that ASAP.
>>>>
>>>> We also love to have examples in the pydy_examples repo. If you can
>>>> contribute soem things you already created that would be great. We will soon
>>>> be making an example gallery like matplotlib.org has but for PyDy problems
>>>> all with nice 3D visualizations embedded in the website. If you want to
>>>> initiate that, that is also an easy way to get started.
>>>>
>>>> - What would the next step be in applying? I see your "umbrellaed" under
>>>> both Python and Sympy. I assume since sympy is closer I'd need to talk to
>>>> them?
>>>>
>>>> Key things: 1) Start on a patch in SymPy (browse through issues for
>>>> ideas), or ask here. 2) Browse the ideas page and start asking questions
>>>> about it here. 3) Start developing a proposal here:
>>>> https://github.com/pydy/pydy/wiki and ask for out feedback. The sooner the
>>>> better.
>>>>
>>>>
>>>> Jason
>>>> moorepants.info
>>>> +01 530-601-9791
>>>>
>>>>
>>>> On Sun, Mar 2, 2014 at 11:44 AM, James Crist <cris...@umn.edu> wrote:
>>>>>
>>>>> @Gilbert:
>>>>>
>>>>> My research is on estimation and control of electromechanical systems.
>>>>> The basic workflow  I use is:
>>>>> 1.) Derive the equations of motion in sympy.
>>>>> 2.) Lamdify them into numpy expressions, for faster simulation
>>>>> 3.) System id using an in-shop python tool for data collection and
>>>>> analysis (cleaning this up to put on github soon)
>>>>> 4.) Develop a controller using python-control
>>>>>
>>>>> I haven't used any of the pydy specific toolset, but have used the
>>>>> Lagrange's method code in classical mechanics to derive equations of motion
>>>>> for a robot for modeling, analysis, and design optimization. I was never
>>>>> taught kane's method (although it looks super useful), so I can't comment on
>>>>> the code for that. The lagrange's method code is slick though.
>>>>>
>>>>> @Aaron:
>>>>>
>>>>> Thanks for the heads up, I'll see what I can contribute.
>>>>>
>>>>> @Jason:
>>>>>
>>>>> Thanks for the paper, I'll give it a read later today.
>>>>>
>>>>> -Jim
>>>>>
>>>>>
>>>>> On Saturday, March 1, 2014 5:50:47 PM UTC-6, moorepants wrote:
>>>>>>
>>>>>> James,
>>>>>>
>>>>>> Here is the paper. I'll respond more in a bit.
>>>>>>
>>>>>>
>>>>>> Jason
>>>>>> moorepants.info
>>>>>> +01 530-601-9791
>>>>>>
>>>>>>
>>>>>> On Fri, Feb 28, 2014 at 4:55 PM, Aaron Meurer <asme...@gmail.com>
>>>>>> wrote:
>>>>>>>
>>>>>>> "Umbrella'd" means that the PyDy GSoC proposals will be accepted
>>>>>>> under
>>>>>>> one of those two orgs, since PyDy wan't accepted as its own org. What
>>>>>>> it means for you is that you will need to satisfy the application
>>>>>>> requirements of whichever org you apply to. I don't know if PSF has
>>>>>>> any requirements, but for SymPy, you need to have a patch merged. See
>>>>>>> https://github.com/sympy/sympy/wiki/gsoc-2014-application-template.
>>>>>>> Since you've already been using SymPy, a good option for this is to
>>>>>>> try to fix some bug that you've come across, or implement some
>>>>>>> missing
>>>>>>> feature that you would have found nice.
>>>>>>>
>>>>>>> Aaron Meurer
>>>>>>>
>>>>>>> On Fri, Feb 28, 2014 at 1:42 PM, Gilbert Gede <gilbe...@gmail.com>
>>>>>>> wrote:
>>>>>>> > Hi James,
>>>>>>> > Glad to see someone else is using our tools! Do you mind sharing
>>>>>>> > some of
>>>>>>> > what you've used it for? Seeing how others are using the software
>>>>>>> > will help
>>>>>>> > us to improve it.
>>>>>>> > I'm not going to be very involved in GSoC this year, but Jason will
>>>>>>> > be
>>>>>>> > (although he's preparing for a conference that is in 4 days).
>>>>>>> > You'll want to
>>>>>>> > follow up with him what projects would be a good fit for you.
>>>>>>> > Regarding the paper, it is still under review... so building the
>>>>>>> > tex file is
>>>>>>> > your best option right now.
>>>>>>> >
>>>>>>> > -Gilbert
>>>>>>> >
>>>>>>> >
>>>>>>> >
>>>>>>> > On Thu, Feb 27, 2014 at 9:34 PM, James Crist <cris...@umn.edu>
>>>>>>> > wrote:
>>>>>>> >>
>>>>>>> >> Hello!
>>>>>>> >>
>>>>>>> >> I'm interested in your project, having used the sympy mechanics,
>>>>>>> >> numpy/scipy, python-control stack for most of my homework and
>>>>>>> >> research. I'm
>>>>>>> >> a Mechanical Engineering Master's student focusing on system
>>>>>>> >> dynamics and
>>>>>>> >> controls at the University of Minnesota.
>>>>>>> >>
>>>>>>> >> Specifically, I'm interested in the "Linearization" project, or
>>>>>>> >> possibly
>>>>>>> >> in the "Efficient Code Generation". I've worked extensively with
>>>>>>> >> Python (4+
>>>>>>> >> years experience), and have some experience in other languages as
>>>>>>> >> well
>>>>>>> >> (Julia, C, Matlab, and a little bit of Fortran). Unfortunately, I
>>>>>>> >> just got
>>>>>>> >> on Github, so I don't have much to show for my work (there's no
>>>>>>> >> python code
>>>>>>> >> up there for starters): https://github.com/jcrist.
>>>>>>> >>
>>>>>>> >> Questions:
>>>>>>> >> - Is there a way I can find the paper "A linearization procedure
>>>>>>> >> for
>>>>>>> >> constrained multibody systems"? University database didn't turn
>>>>>>> >> anything up,
>>>>>>> >> but googling brought up the tex file. I suppose I could read that.
>>>>>>> >> - What can I do to get involved beforehand/even if I don't get
>>>>>>> >> accepted?
>>>>>>> >> As I said, I make heavy use of sympy mechanics and the whole
>>>>>>> >> python science
>>>>>>> >> stack in my research, and would love to contribute back.
>>>>>>> >> - What would the next step be in applying? I see your "umbrellaed"
>>>>>>> >> under
>>>>>>> >> both Python and Sympy. I assume since sympy is closer I'd need to
>>>>>>> >> talk to
>>>>>>> >> them?
>>>>>>> >>
>>>>>>> >> Thanks,
>>>>>>> >> - Jim
>>>>>>> >>
>>>>>>> >> --
>>>>>>> >> You received this message because you are subscribed to the Google
>>>>>>> >> Groups
>>>>>>> >> "PyDy" group.
>>>>>>> >> To unsubscribe from this group and stop receiving emails from it,
>>>>>>> >> send an
>>>>>>> >> email to pydy+uns...@googlegroups.com.
>>>>>>> >> To post to this group, send email to py...@googlegroups.com.
>>>>>>>
>>>>>>> >> Visit this group at http://groups.google.com/group/pydy.
>>>>>>> >> For more options, visit https://groups.google.com/groups/opt_out.
>>>>>>> >
>>>>>>> >
>>>>>>> > --
>>>>>>> > You received this message because you are subscribed to the Google
>>>>>>> > Groups
>>>>>>> > "PyDy" group.
>>>>>>> > To unsubscribe from this group and stop receiving emails from it,
>>>>>>> > send an
>>>>>>> > email to pydy+uns...@googlegroups.com.
>>>>>>> > To post to this group, send email to py...@googlegroups.com.
>>>>>>>
>>>>>>> > Visit this group at http://groups.google.com/group/pydy.
>>>>>>> > For more options, visit https://groups.google.com/groups/opt_out.
>>>>>>>
>>>>>>> --
>>>>>>> You received this message because you are subscribed to the Google
>>>>>>> Groups "PyDy" group.
>>>>>>> To unsubscribe from this group and stop receiving emails from it,
>>>>>>> send an email to pydy+uns...@googlegroups.com.
>>>>>>> To post to this group, send email to py...@googlegroups.com.
>>>>>>>
>>>>>>> Visit this group at http://groups.google.com/group/pydy.
>>>>>>> For more options, visit https://groups.google.com/groups/opt_out.
>>>>>>
>>>>>>
>>>>> --
>>>>> You received this message because you are subscribed to the Google
>>>>> Groups "PyDy" group.
>>>>> To unsubscribe from this group and stop receiving emails from it, send
>>>>> an email to pydy+uns...@googlegroups.com.
>>>>> To post to this group, send email to py...@googlegroups.com.
>>>>> Visit this group at http://groups.google.com/group/pydy.
>>>>> For more options, visit https://groups.google.com/groups/opt_out.
>>>>
>>>>
>>> --
>>> You received this message because you are subscribed to the Google Groups
>>> "PyDy" group.
>>> To unsubscribe from this group and stop receiving emails from it, send an
>>> email to pydy+uns...@googlegroups.com.
>>> To post to this group, send email to py...@googlegroups.com.
>>> Visit this group at http://groups.google.com/group/pydy.
>>> For more options, visit https://groups.google.com/groups/opt_out.
>>
>>
> --
> You received this message because you are subscribed to the Google Groups
> "PyDy" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to pydy+uns...@googlegroups.com.
> To post to this group, send email to py...@googlegroups.com.
> Visit this group at http://groups.google.com/group/pydy.
> For more options, visit https://groups.google.com/d/optout.

--
You received this message because you are subscribed to the Google Groups "PyDy" group.
To unsubscribe from this group and stop receiving emails from it, send an email to pydy+uns...@googlegroups.com.
To post to this group, send email to py...@googlegroups.com.
Visit this group at http://groups.google.com/group/pydy.
For more options, visit https://groups.google.com/d/optout.

--
You received this message because you are subscribed to the Google Groups "PyDy" group.
To unsubscribe from this group and stop receiving emails from it, send an email to pydy+uns...@googlegroups.com.
To post to this group, send email to py...@googlegroups.com.
Visit this group at http://groups.google.com/group/pydy.
For more options, visit https://groups.google.com/d/optout.

[Aaron Meurer]

In a similar vein:

- Expect to submit a pull request as often as possible. Whenever you
have workable code, you should submit a PR. What you want to avoid is
one mega-pull request at the end of the summer that is too big for
anyone to review.

- Each PR should come with tests and documentation. It's much easier
to test and document as you go than to do it at the end. Plus, if you
do this, you will find yourself avoiding bad design much sooner than
you otherwise would.

Aaron Meurer

[James Crist]

@Jason:

Thoughts on some of your points:

- You're correct, that was a mistake in my attempt to convert math to sympy notation for the proposal. Should be fixed now.

- I'm not sure about the best way to divy out the methods into a manageable module for linearize. I want the function call to be the same regardless of system type. The printing modules use a class, dsolve uses functions. Both have their merits, but I think a class based approach in the end seems cleaner. Here are my thoughts on ways to accomplish this:

A class "Linearizer", with an init like:

linearizer = Linearizer(expression, q, *args, **kwargs)

The main method is then "linearize_about(q_eq, *args, **kwargs), which would handle the calculations. There are 3 ways I can think of using this class:

1.) Linearizer class has submethods for each type. New types could be added with new submethods, and a control path in `linearize_about` that would result in calling that submethod.

2.) Linearizer class would be a base class, and special types would subclass it, and modify `linearize_about` to add the required functionality.

3.) Linearizer class would be an intermediate form, and each system type would be transformed to fit the data structure. `linearize_about` could then be used for all types without modification.

Choice 1 would result in fewer classes, but possibly a more complicated code base. With generality in mind, choice 2 could be implemented using only a couple classes to cover the major types. Choice 3, while desirable, would likely result in a complicated codebase due to the wide variance of system types. However, some research needs to be done into the best way of grouping.

- PDE stuff has been removed. If the linearize code is designed in a modular way, this should be able to be added later when we have pde eoms being generated.

- I'm actually thinking it may be best to transform all inputs to the form "f(q, dq, t) = 0", and then linearize f. With a linear form, it should be possible to rearrange the resulting expression into the desirable output form. This may not be true for all problems though.

- After some thought, I'm not sure if there's a general way to output A & B. For example, if linearizing about a non-zero equilibrium point, the taylor series is:

dx = f(x_bar, u_bar) + df/dx(x_bar, u_bar) * (x - x_bar) + df/du(x_bar, u_bar) * (u - u_bar)

=> dx = dx_bar + A*(x - x_bar) + B*(u - u_bar)

In this case, the A and B matrices represent the dynamics about the equilibrium. Which is different than representing the full system output. Perhaps a third term could be added, so A, B, and dx_bar - (A*x_bar + B*u_bar) are returned?

- Lagranges method can be linearized in process as well. These notes explain it fairly well. I have more on this in a text book from my dynamics class. I just got Kane's book from our inter-university library system (this book is really hard to find, not sure why). Going to take some time and try to teach myself this method over the next few weeks.

- Simplification could definitely be done in route. I've also noticed, especially for lagrange's equations, that doing the rhs() method with inv_method="GE" gives a far simpler expression most of the time. Not sure why, but that could definitely be looked into.

- I like the idea of piecewise constraints, or EOM. It would be nice to be able to handle this. Without much thought, I think the best way would be to do a separate linear form for each piecewise range. Or perhaps matrices with piecewise terms. Some thought will have to be given to this.

- Documentation and tests will definitely be added each step. I just think that giving an entire week at the end to ensure this is done, and give it a final proof reading seems to be a good idea.

----------------------------

On my application on melange you mentioned some opportunities for "additions to simplification, code gen, and cse" that could be part of this project? Code gen for sure makes sense, to add the capability for nice matrix operations generated by the linearized eom. CSE could also be nice to make these more efficient in a higher level language like python. I think c compilers may be smart enough to handle this on their own?

I'm curious to hear your thoughts on this.

-Jim

[Gilbert Gede]

Hey Jim,

Here are Kane's books, from Cornell:

http://ecommons.library.cornell.edu/handle/1813/638

http://ecommons.library.cornell.edu/handle/1813/637

Regarding some of the linearization things:

When linearizing about a point, I believe the results are only valid if that points satisfies the differential equations: i.e.  x' - f(x) = 0. I don't think we currently have a "clean" way to do this. It would also be applicable to choosing initial conditions for simulation of any constrained systems.

Also, our linearization method assumes a certain form of the equations. I might be cool to have a fallback approach though. What we're effectively doing when constraining a system is changing the dependent q/u's to be functions of independent q/u's - but no one ever rewrites them as such (e.g., u_dep(q_ind, u_ind, t)). As we know the structure of Kane's equations, we can effectively manually apply the chain rule to the differentiation that needs to happen as part of the linearization process (the current approach). Another interesting approach might be (and I'm not 100% sure it would work) would be redefining the _eval_derivative function for each dependent dynamicsymbol, so that when it has it's derivative taken wrt an independent symbol, it returns the proper result (ie d u_dep(q_ind, u_ind, t) / d q_ind). It seems like it might be a more flexible approach.

-Gilbert

[James Crist]

Gilbert,

Thanks for the books!

Responses:

- Yes, linearization is only valid about equilibrium points, or trajectories (points that satisfy the initial equation). Simply plugging in and checking may be the best way of ensuring this.

- The idea of redefining __eval__derivative method for the dependent symbols to be functions of the independent is interesting. I'm not sure if there is a robust, general way of determining the dependent and independent symbols for all problems. I'll have to think on this for a bit. If possible, it would result in a more flexible approach for systems that were just odes with constraints, which is the majority of what we deal with. Definitely worth looking into.

-Jim

...

[James Crist]

@Jason (or anyone),

I hate to pester, but I only have a day or so to make changes, so if you have some time could you respond to my above comments/questions? Thanks.

-Jim

>>>>>> &gt

...

[Jason Moore]

I'll get to them this morning.

Jason
moorepants.info
+01 530-601-9791

--

You received this message because you are subscribed to the Google Groups "PyDy" group.
To unsubscribe from this group and stop receiving emails from it, send an email to pydy+uns...@googlegroups.com.
To post to this group, send email to py...@googlegroups.com.
Visit this group at http://groups.google.com/group/pydy.

For more options, visit https://groups.google.com/d/optout.

[Jason Moore]

Jason
moorepants.info
+01 530-601-9791

On Tue, Mar 18, 2014 at 2:16 PM, James Crist <cris...@umn.edu> wrote:

@Jason:

Thoughts on some of your points:

- You're correct, that was a mistake in my attempt to convert math to sympy notation for the proposal. Should be fixed now.

- I'm not sure about the best way to divy out the methods into a manageable module for linearize. I want the function call to be the same regardless of system type. The printing modules use a class, dsolve uses functions. Both have their merits, but I think a class based approach in the end seems cleaner. Here are my thoughts on ways to accomplish this:

A class "Linearizer", with an init like:

linearizer = Linearizer(expression, q, *args, **kwargs)

The main method is then "linearize_about(q_eq, *args, **kwargs), which would handle the calculations. There are 3 ways I can think of using this class:

1.) Linearizer class has submethods for each type. New types could be added with new submethods, and a control path in `linearize_about` that would result in calling that submethod.

2.) Linearizer class would be a base class, and special types would subclass it, and modify `linearize_about` to add the required functionality.

3.) Linearizer class would be an intermediate form, and each system type would be transformed to fit the data structure. `linearize_about` could then be used for all types without modification.

Choice 1 would result in fewer classes, but possibly a more complicated code base. With generality in mind, choice 2 could be implemented using only a couple classes to cover the major types. Choice 3, while desirable, would likely result in a complicated codebase due to the wide variance of system types. However, some research needs to be done into the best way of grouping.

These all sound like potentially useful methods. Maybe mention them all as you have here, their pros/cons, and maybe detail the best one with some sample/psuedo code showing how a user would use it. These details can and will be refined if you are selected.

 

- PDE stuff has been removed. If the linearize code is designed in a modular way, this should be able to be added later when we have pde eoms being generated.

- I'm actually thinking it may be best to transform all inputs to the form "f(q, dq, t) = 0", and then linearize f. With a linear form, it should be possible to rearrange the resulting expression into the desirable output form. This may not be true for all problems though.

This could be done, but there may be advantages to the form M x' = f. The linearization paper we shared shows this. Because M is only a function of q means that all the dM/dq' are zero, no need to try to compute them. But there are use cases where people want the jacobian of f(q, dq, t) = 0. It is useful in backward Euler integration, for example. I can get some specific use cases for this. From the non-linear code gen side of things, I've been trying to devise a way to pass in multiple forms of the EoMs. This will probably be similar. We basically need to have a list of forms of non-linear EoMs that we produce and use cases for the linear EoMs to decide this.

 

- After some thought, I'm not sure if there's a general way to output A & B. For example, if linearizing about a non-zero equilibrium point, the taylor series is:

dx = f(x_bar, u_bar) + df/dx(x_bar, u_bar) * (x - x_bar) + df/du(x_bar, u_bar) * (u - u_bar)

=> dx = dx_bar + A*(x - x_bar) + B*(u - u_bar)

In this case, the A and B matrices represent the dynamics about the equilibrium. Which is different than representing the full system output. Perhaps a third term could be added, so A, B, and dx_bar - (A*x_bar + B*u_bar) are returned?

I'm not following that. I'm pretty sure the A and B matrices can be formed in the same fashion regardless if you equilibrium point is zero. If the EoMs are:

M x' = f

Then the Taylor series of f about x_eq and u_eq turn this into a function

M|x_eq x_eq' = F_a x_eq + F_b u_eq

Then

A = M|x_eq^-1 F_a

B = M|x_eq^-1 F_b

I'm not sure what the dx_bar term is that you show.

 

- Lagranges method can be linearized in process as well. These notes explain it fairly well. I have more on this in a text book from my dynamics class. I just got Kane's book from our inter-university library system (this book is really hard to find, not sure why). Going to take some time and try to teach myself this method over the next few weeks.

Awesome! That will be a great addition. Kane's book has been out of print for a while. So there are few copies floating around. I just printed the PDF from Cornell originally. Now that I have a job I bought the real deal! It may be useful for one of us to Skype with you. We could explain the method in 30 minutes or so (at least enough to get you started).

 

- Simplification could definitely be done in route. I've also noticed, especially for lagrange's equations, that doing the rhs() method with inv_method="GE" gives a far simpler expression most of the time. Not sure why, but that could definitely be looked into.

Simplification in route should always be optional. At some point in the EoM derivations the simplifications become to compute intensive. It would be nice to identify the bottlenecks and either avoid them or fix them. Simplification of these really long expressions will be hard to fix. I don't think general simplification routines are that efficient. The Fu method was recently added to SymPy which does a nice job, but the speed is still rough. This could be something that would be nice to have from a CSympy core. I'm not sure what the other methods for the inverse are in SymPy. We could compare them. Note that the Cholesky factorization can be used for many dynamics problems (if M is symmetric positive definite). It is a little faster than Gaussian Elimination but I don't know what the simplification issues are. Keep in mind that forming the right hand side symbolically is really only useful for small problems. In general, we do this numerically.

 

- I like the idea of piecewise constraints, or EOM. It would be nice to be able to handle this. Without much thought, I think the best way would be to do a separate linear form for each piecewise range. Or perhaps matrices with piecewise terms. Some thought will have to be given to this.

We need to have a some example problems that have piecewise parts that we may want to linearize and then work on a solution. Concrete examples are needed first.

 

- Documentation and tests will definitely be added each step. I just think that giving an entire week at the end to ensure this is done, and give it a final proof reading seems to be a good idea.

Yes thats fine. But you must make sure to do this along the way at the primary way. Students often neglect this.

 

----------------------------

On my application on melange you mentioned some opportunities for "additions to simplification, code gen, and cse" that could be part of this project? Code gen for sure makes sense, to add the capability for nice matrix operations generated by the linearized eom. CSE could also be nice to make these more efficient in a higher level language like python. I think c compilers may be smart enough to handle this on their own?

We want to be able to generate code in a variety of languages, both high and low. SymPy does this already. Right now the code printing, code gen, and autowrap paths in SymPy do not have nice support for matrices. For the linear systems you will want to generate code that efficiently evaluates matrices (A and B mostly).

The CSE is definitely useful for speeding up high level languages. But I've found that CSE also speeds up C code (haven't delved into the depths of complier optimization though and whether pre-CSE does help). You can test this with pydy-code-gen. I did this once, will have to dig up the results. It is probably in a SymPy or PyDy list email.

[James Crist]

Thanks for the replies. Some quick responses:

> I'm not following that. I'm pretty sure the A and B matrices can be formed in the same fashion regardless if you equilibrium point is zero. If the EoMs are:

>

>M x' = f
>

>Then the Taylor series of f about x_eq and u_eq turn this into a function
>

>M|x_eq x_eq' = F_a x_eq + F_b u_eq
>

>Then
>
>A = M|x_eq^-1 F_a

>B = M|x_eq^-1 F_b
>

>I'm not sure what the dx_bar term is that you show.

That's because there isn't one, of course dx_bar is always zero for equilibrium points (I was wrong). What I meant was, if the eq. point is non-zero (i.e. x_eq != 0, u_eq != 0), then the linearized equations are of the dynamics about the equilibrium point:

dx = A(x-x_eq) + B(u - u_eq) = A*x + B*u + (-A*x_eq - B*u_eq)

What I'm wondering, is if it would make sense to return that last term as well (call it Xeq), so that the user could in all cases say:

A, B, Xeq = linearize(expression, state_space=True)
dx_linearized = A*x + B*u + Xeq

For equilibrium points at 0, Xeq will be zero, but this isn't true for all equilibrium points.

>This could be done, but there may be advantages to the form M x' = f. The linearization paper we shared shows this. Because M is only a function of q means that all the dM/dq' >are zero, no need to try to compute them. But there are use cases where people want the jacobian of f(q, dq, t) = 0. It is useful in backward Euler integration, for example. I can >get some specific use cases for this. From the non-linear code gen side of things, I've been trying to devise a way to pass in multiple forms of the EoMs. This will probably be >similar. We basically need to have a list of forms of non-linear EoMs that we produce and use cases for the linear EoMs to decide this.

Indeed, if we can clearly get them into the form M x' = f this makes sense. For both Kanes and Lagranges method, this is already implemented. I'm trying to generalize to a system of equations passed in by the user as well, so that the expression could be generated somewhere else, but the linearize code could still be useful. This may not be as simple as I thought it to be. Starting of a list of forms:

- M(q)*q'' = f(q, q', t)
- f(q, q', t) = 0
- dx = f(x, u, t)         # This is basically f(q, dq, t) rearranged, but is a common form in control theory, and may be desired

That's all I can think of for now.

>Simplification in route should always be optional. At some point in the EoM derivations the simplifications become to compute intensive. It would be nice to identify the >bottlenecks and either avoid them or fix them. Simplification of these really long expressions will be hard to fix. I don't think general simplification routines are that efficient. The >Fu method was recently added to SymPy which does a nice job, but the speed is still rough. This could be something that would be nice to have from a CSympy core. I'm not >sure what the other methods for the inverse are in SymPy. We could compare them. Note that the Cholesky factorization can be used for many dynamics problems (if M is >symmetric positive definite). It is a little faster than Gaussian Elimination but I don't know what the simplification issues are. Keep in mind that forming the right hand side >symbolically is really only useful for small problems. In general, we do this numerically.

I don't know much about simplification routines, I'll have to spend some time experimenting and reading the code. From experience with sympy, it seems that the longest computation times I've experienced have been simplification, and matrix inversion. I agree, it would be nice to get this in the CSympy core.

>We need to have a some example problems that have piecewise parts that we may want to linearize and then work on a solution. Concrete examples are needed first.

Of course. I'll probably move this to an "if I have time" goal. With most eom being formed by kanes and lagranges equations, this seems less necessary than good linearization, and code generation.

>We want to be able to generate code in a variety of languages, both high and low. SymPy does this already. Right now the code printing, code gen, and autowrap paths in >SymPy do not have nice support for matrices. For the linear systems you will want to generate code that efficiently evaluates matrices (A and B mostly).

>The CSE is definitely useful for speeding up high level languages. But I've found that CSE also speeds up C code (haven't delved into the depths of complier optimization though >and whether pre-CSE does help). You can test this with pydy-code-gen. I did this once, will have to dig up the results. It is probably in a SymPy or PyDy list email.

Should this be part of the main project, or a "if I have time" proposal? For pure c code, would numerical library support, such as eigen or lapack, be wanted? If we're solving numerically for M and f, and then finding dx = M^-1 f, this may be desired.

----

Big question now is, how should this be prioritized? My thoughts are:

1.) Linearization for general problems with constraints (use with Kane and Lagrange, as well as equations from other libraries. No PDE, or piecewise)
2.) Pre-linearization for Kane and Lagrange
3.) Fix the easy issues with simplification
4.) Code Gen for matrix operations
5.) Fix the hard issues with simplification
6.) Piecewise constraints.

I'll try to finish up my proposal by tonight and update it on Melange. Thanks for the help!

-Jim

Jason
moorepants.info
+01 530-601-9791

>>>&g

...

[James Crist]

Woops. Tried to format the things I was replying to with ">" before every line, but that clearly ended up horribly. Sorry about that, hope everything is still readable.

-Jim

...

[Jason Moore]

Jason
moorepants.info
+01 530-601-9791

On Thu, Mar 20, 2014 at 12:34 PM, James Crist <cris...@umn.edu> wrote:

Thanks for the replies. Some quick responses:

> I'm not following that. I'm pretty sure the A and B matrices can be formed in the same fashion regardless if you equilibrium point is zero. If the EoMs are:

>

>M x' = f
>

>Then the Taylor series of f about x_eq and u_eq turn this into a function
>

>M|x_eq x_eq' = F_a x_eq + F_b u_eq
>

>Then
>
>A = M|x_eq^-1 F_a

>B = M|x_eq^-1 F_b
>

>I'm not sure what the dx_bar term is that you show.

That's because there isn't one, of course dx_bar is always zero for equilibrium points (I was wrong). What I meant was, if the eq. point is non-zero (i.e. x_eq != 0, u_eq != 0), then the linearized equations are of the dynamics about the equilibrium point:

dx = A(x-x_eq) + B(u - u_eq) = A*x + B*u + (-A*x_eq - B*u_eq)

What I'm wondering, is if it would make sense to return that last term as well (call it Xeq), so that the user could in all cases say:

A, B, Xeq = linearize(expression, state_space=True)
dx_linearized = A*x + B*u + Xeq

For equilibrium points at 0, Xeq will be zero, but this isn't true for all equilibrium points.

I think it is useful to pass in and/or out the list of the linear states and inputs as symbolic values. If you actually change the variable names from x to del_x (x - x_eq) then you should spit those out for use by the user. That is the mathematically correct terminology. The current Kane linearize function is subtle about that detail.

I might be getting lost here. We really need to be able to type in rendered math to make sure we are talking about the same thing.

delx = x - xeq
delu = u - ueq

So the equations should be in this form:

delx' = A delx + B delu

xeq is all that is needed to be passed in to the linearizer (along with the nonlinear eoms).

 

>This could be done, but there may be advantages to the form M x' = f. The linearization paper we shared shows this. Because M is only a function of q means that all the dM/dq' >are zero, no need to try to compute them. But there are use cases where people want the jacobian of f(q, dq, t) = 0. It is useful in backward Euler integration, for example. I can >get some specific use cases for this. From the non-linear code gen side of things, I've been trying to devise a way to pass in multiple forms of the EoMs. This will probably be >similar. We basically need to have a list of forms of non-linear EoMs that we produce and use cases for the linear EoMs to decide this.

Indeed, if we can clearly get them into the form M x' = f this makes sense. For both Kanes and Lagranges method, this is already implemented. I'm trying to generalize to a system of equations passed in by the user as well, so that the expression could be generated somewhere else, but the linearize code could still be useful. This may not be as simple as I thought it to be. Starting of a list of forms:

- M(q)*q'' = f(q, q', t)
- f(q, q', t) = 0
- dx = f(x, u, t)         # This is basically f(q, dq, t) rearranged, but is a common form in control theory, and may be desired

That's all I can think of for now.

That's probably it. Except we can do some thinking about the constraint equations too and whether we want to explicitly deal with them. Luke and Gilberts paper may touch on this. Luke has some good ideas about all this but he hasn't piped in yet...

 

>Simplification in route should always be optional. At some point in the EoM derivations the simplifications become to compute intensive. It would be nice to identify the >bottlenecks and either avoid them or fix them. Simplification of these really long expressions will be hard to fix. I don't think general simplification routines are that efficient. The >Fu method was recently added to SymPy which does a nice job, but the speed is still rough. This could be something that would be nice to have from a CSympy core. I'm not >sure what the other methods for the inverse are in SymPy. We could compare them. Note that the Cholesky factorization can be used for many dynamics problems (if M is >symmetric positive definite). It is a little faster than Gaussian Elimination but I don't know what the simplification issues are. Keep in mind that forming the right hand side >symbolically is really only useful for small problems. In general, we do this numerically.

I don't know much about simplification routines, I'll have to spend some time experimenting and reading the code. From experience with sympy, it seems that the longest computation times I've experienced have been simplification, and matrix inversion. I agree, it would be nice to get this in the CSympy core.

We should just get everythign working without simplification for the fastest computations and then we can see if it is useful to add simplification in at certain places in the computation workflow.

 

>We need to have a some example problems that have piecewise parts that we may want to linearize and then work on a solution. Concrete examples are needed first.

Of course. I'll probably move this to an "if I have time" goal. With most eom being formed by kanes and lagranges equations, this seems less necessary than good linearization, and code generation.

Yes, these kinds of problems can be in the "if I have time".

 

>We want to be able to generate code in a variety of languages, both high and low. SymPy does this already. Right now the code printing, code gen, and autowrap paths in >SymPy do not have nice support for matrices. For the linear systems you will want to generate code that efficiently evaluates matrices (A and B mostly).

>The CSE is definitely useful for speeding up high level languages. But I've found that CSE also speeds up C code (haven't delved into the depths of complier optimization though >and whether pre-CSE does help). You can test this with pydy-code-gen. I did this once, will have to dig up the results. It is probably in a SymPy or PyDy list email.

Should this be part of the main project, or a "if I have time" proposal? For pure c code, would numerical library support, such as eigen or lapack, be wanted? If we're solving numerically for M and f, and then finding dx = M^-1 f, this may be desired.

I think generating the code should be in. You can pass in the linear equations to pydy-code-gen right now for a "dumb" code gen (i.e. no optimizations for the linear form) and it will work. So some thought can be done to improve that for linear systems.

We use lapack through scipy right now. I've got unmerged work that calls the lapack routines from cython code. But ultimately we want to generate code for any underlying lib if the user desires.

 

----

Big question now is, how should this be prioritized? My thoughts are:

1.) Linearization for general problems with constraints (use with Kane and Lagrange, as well as equations from other libraries. No PDE, or piecewise)
2.) Pre-linearization for Kane and Lagrange
3.) Fix the easy issues with simplification
4.) Code Gen for matrix operations
5.) Fix the hard issues with simplification
6.) Piecewise constraints.

I'd move 4 above 3, but other than that looks good!

 

--

[James Crist]

Made some edits. I'm going to give it a final proof read when I get home, but this is (I think) pretty close to done. Let me know what you think:

https://github.com/pydy/pydy/wiki/GSoC-2014-Application:-Jim-Crist-%28Linearization-Routines-for-Equations-of-Motion%29

-Jim

...