Swinging chain

[moorepants]

https://github.com/pydy/pst-notebooks/blob/main/chain-link.ipynb

A long time ago, I saw a swinging chain as an example of some simulation system, I do not recall which one. Having nothing better to do, I tried to do it with sympy mechanics.

I modelled it as a pendulum consisiting of n links, suspended at one end.
I modelled the links as thin rods.
The total length of the pendulum is l, hence each link has length l/n.

I noticed the following:

When I defined the rotation of frame A[i] with respect to A[i-1] ( A[0] is the frame attached to the suspension point ), the number of operations shot up. For n = 15, the mass matrix had 25 mio operations.
When I defined the rotation of each frame A[i] w.r.t. the inertial frame O, the operations dropped a lot: For n = 50, MM has 40,000 operations, force has 50,000 operations. Normal numbers.

I guess, this has an explanation similar to the one with orient_body_fixed, which Timo explained to me a while back.

As per my experience, it takes a long time to set up Kane's equations.
For a system with 40,000 operations in MM and 50,000 in force, this should take maybe 20 sec.
Here it took about 120 sec.
I have no explanation of this.

Also the lambdification took over 20 sec, unusually long for such a small system.
The numerical integration also seemed to be on the longer side.

Anyway, my program shows, that such a chain can easily be modelled with symy mechanics!

[tjstienstra]

Here are my quick thoughts about the remarks.

1. When computing the generalized forces, the velocity and acceleration of all masses and frames is needed with respect to the inertial frame. If you express each frame with respect to the previous in the chain, then those velocities and acceleration of that frame will also include all the previous rotations. If possible it is always better to directly express things with respect to the inertial frame. Would have to say that CSE should reduce the number of operations massively, though it will stay a little higher I think.
2. We could profile you program, though it also depends on what sympy version you are currently using. master should have a bit faster KanesMethod (though that speed up should not explain this difference, just noticed it is not even in 1.12). If I have to guess it probably has to do with this system consisting out of a lot more small parts, so there just are a lot of small computations.
3. Lambdification should scale with the number of operations, so would have expected that that would be similar. The numerical integration may indeed be slightly slower, as the number of operations after CSE may actually be relatively high compared to systems with similar size in terms of operations in the mass matrix and forcing vector. Besides that takes exponentially more time to solve a larger system matrix, so that should also slow down the integration.

[Peter230655]

Dear Timo,

This makes sense to me after what you explained to me earlier regarding setting up the kinematic equations.

On my iPad ( the only place I can use sympy until I retire and buy a computer) I use sympy 1.12

Do I understand you correctly: In this example, there are not so many "repetitions", so applying cse=True does not really speed up integration so much?

Thanks, Peter

[tjstienstra]

Do I understand you correctly: In this example, there are not so many "repetitions", so applying cse=True does not really speed up integration so much?

It still gives a speed up, rather sure of that. However I would expect that the speed up is lower when comparing it to the time gain you'll have on for example a bicycle or a rolling axle.

[Peter230655]

When 'cse' first came out, I compared the integration speeds with cse=True/False on various examples. I saw improvements by a factor 100 and more. But these were generally 3D 'monsters' with up to 1,000,000 operations in the force.
But then I stopped these comparisons.

[tjstienstra]

Yeah, it is something you'd better always turn on. Main points was that this is probably a system on which the effect is a bit lower.

[Peter230655]

I will run both ways on this pendulum as soon as I will get to my hotel and let you know the difference in this case.

[tjstienstra]

Just checked quickly myself. The cse on the mass matrix only gives an reduction of 39180 to 5293 and 49281 to 10446 for the forcing vector. This is as expected a low reduction when compared with other systems. However it might be worth profiling the sympy internals with for example pyinstrument.

[Peter230655]

Thanks!
If it is easy to explain, could you tell me how you did this? Thanks!

[tjstienstra]

I just compared sm.count_ops(MM) and sm.count_ops(sm.cse(MM)) to see the reduction for the mass matrix.

[Peter230655]

Thanks!
I will try this on a few of my examples.

[Peter230655]

I tried on my 'rolling disc on very uneven street'. Reduction by about a factor of 40 in the force vector.
I am very unclear, what exactly cse does. Where could I read up on it?
Thanks as always, Timo!

[tjstienstra]

CSE stands for Common Subexpression Elimination. Here is the wikipedia page. It is overall a relatively simple algorithm where you make sure that if you have a + b two times in you equation, then you make a temporary variable _x0 = a + b, which you then compute once and substitute in both locations of the equation.

[Peter230655]

Thanks, Timo!
You simply know everything about sympy, it seems to me!
Thanks for your patience with me!

[tjstienstra]

You're welcome, there is actually really a lot I don't know about sympy, but I have read quite a few parts of the source code and played around with it quite often.