Multibody Simulation using sympy, scipy and vpython
This is a modal window.
The media could not be loaded, either because the server or network failed or because the format is not supported.
Formal Metadata
| Title |
| |
| Title of Series | ||
| Part Number | 170 | |
| Number of Parts | 173 | |
| Author | ||
| License | CC Attribution - NonCommercial - ShareAlike 3.0 Unported: You are free to use, adapt and copy, distribute and transmit the work or content in adapted or unchanged form for any legal and non-commercial purpose as long as the work is attributed to the author in the manner specified by the author or licensor and the work or content is shared also in adapted form only under the conditions of this | |
| Identifiers | 10.5446/20177 (DOI) | |
| Publisher | ||
| Release Date | ||
| Language | ||
| Production Place | Bilbao, Euskadi, Spain |
Content Metadata
| Subject Area | ||
| Genre | ||
| Abstract |
| |
| Keywords |
EuroPython 2015170 / 173
5
6
7
9
21
27
30
32
36
37
41
43
44
45
47
51
52
54
55
58
63
66
67
68
69
72
74
75
77
79
82
89
92
93
96
97
98
99
101
104
108
111
112
119
121
122
123
131
134
137
138
139
150
160
165
167
173
00:00
BackupComputer fontComputer simulationInformationAssembly languageLecture/ConferenceXMLComputer animation
00:30
BackupIntegrated development environmentComputer simulationVarianceLecture/ConferenceComputer animation
01:02
Computer simulationBasis <Mathematik>Product (business)Nichtlineares GleichungssystemLecture/Conference
01:58
Integrated development environmentComputer simulationIndependence (probability theory)Metropolitan area networkNichtlineares GleichungssystemPhysical systemNichtnewtonsche FlüssigkeitInclusion mapProcess (computing)Software developerComputer simulationEuler's equations (rigid body dynamics)Disk read-and-write headScripting languageComputer animation
02:55
Computer simulationMetropolitan area networkCASE <Informatik>Process (computing)Computer simulationConstraint (mathematics)Mechanism designRoboticsExpressionSystem dynamicsAsynchronous Transfer ModeMultiplicationTable (information)Forcing (mathematics)Assembly languageBranch (computer science)Software developerLecture/ConferenceComputer animation
04:20
Computer simulationMetropolitan area networkBranch (computer science)Computer simulationSystem dynamicsLecture/ConferenceComputer animation
04:43
Computer simulationMetropolitan area networkAlgebraBasis <Mathematik>PiProcess (computing)AlgebraLevel (video gaming)Product (business)Symbol tableMechanism designINTEGRALMoment (mathematics)Derivation (linguistics)Right angleComputer animationLecture/Conference
05:54
Metropolitan area networkOrdinary differential equationAlgebraLinear mapComputer programmingNumeral (linguistics)Core dumpLecture/ConferenceXMLComputer animation
06:17
Computer simulationOrdinary differential equationAlgebraMetropolitan area networkOrdinary differential equationComputer animationLecture/Conference
06:41
Computer simulationCuboidSpring (hydrology)System programmingBit rateNichtlineares GleichungssystemGraph (mathematics)Operator (mathematics)Primitive (album)Software testingSpring (hydrology)Computer simulationComputer animation
07:41
Derivation (linguistics)Physical systemBlock (periodic table)Computer simulationMassMoment of inertiaMoment (mathematics)Derivation (linguistics)Physical systemBlock (periodic table)Ordinary differential equationBoiling pointLine (geometry)Degrees of freedom (physics and chemistry)Multiplication signTranslation (relic)Mechanism designTheoryBuildingDifferential (mechanical device)RotationComputer animation
08:18
Derivation (linguistics)Moment of inertiaMoment (mathematics)MassBlock (periodic table)Metropolitan area networkPoint (geometry)Type theoryCartesian coordinate systemPhysical systemComputer simulationMassMoment of inertiaFigurate numberBoiling pointDegrees of freedom (physics and chemistry)NumberPoint (geometry)Graph coloringSkeleton (computer programming)Computer animationLecture/Conference
09:07
Block (periodic table)Cartesian coordinate systemPhysical systemComputer simulationForceMetropolitan area networkData Encryption StandardLogical constantType theorySurface of revolutionFilm editingNichtlineares GleichungssystemMoment <Mathematik>Forcing (mathematics)Computer animation
09:33
ForceBlock (periodic table)Computer simulationData Encryption StandardLogical constantPhysical systemPosition operatorAngleProgrammschleifePoint (geometry)Ordinary differential equationPort scannerCASE <Informatik>NumberConstraint (mathematics)Core dumpSystem callDomain nameNichtlineares GleichungssystemAngleAdditionMultiplication signVibrationProcess (computing)Degrees of freedom (physics and chemistry)Forcing (mathematics)Algebraic equationPhysical systemDifferential algebraic equationComputer simulationMoment <Mathematik>Type theoryGravitationMechanism designDifferent (Kate Ryan album)SurfaceIndependence (probability theory)Cartesian coordinate systemAuditory maskingOffice suiteMathematicsPairwise comparisonCoordinate systemMeasurementFlow separationLoop (music)Special unitary groupChainMultibody systemSet (mathematics)Computer animation
13:21
Port scannerProgrammschleifeBlock (periodic table)Propositional formulaPhysical systemMetropolitan area networkVirtual realityData typeInternet service providerForceOrdinary differential equationFront and back endsMultiplicationGamma functionMathematical analysisJacobi methodLinear mapInclusion mapTransformation (genetics)Computer simulationParameter (computer programming)Logical constantPhysical quantityCoordinate systemEquations of motionConstraint (mathematics)Ordinary differential equationParameter (computer programming)Algebraic equationComplete metric spaceLinearizationNichtlineares GleichungssystemSet (mathematics)Physical systemObject (grammar)NumberSign (mathematics)Reflection (mathematics)ExpressionPoint (geometry)Computer simulationLoop (music)FrequencyVector spaceFormal languageVelocityFood energyMathematical analysisForcing (mathematics)MereologyProcess (computing)PiFunctional (mathematics)Computer animation
16:07
Point (geometry)Internet service providerLogical constantProgrammschleifeForceMathematical analysisLinear mapJacobi methodInclusion mapSet (mathematics)Transformation (genetics)Port scannerParameter (computer programming)Computer simulationPhysical systemMetropolitan area networkMassThetafunktionRepresentation (politics)Computer simulationSimulationExpressionMeasurementFunctional (mathematics)Forcing (mathematics)Physical systemConstraint (mathematics)XMLComputer animation
16:36
Metropolitan area networkPhysical systemData typeMoment of inertiaStorage area networkMach's principleSoftware developerNumeral (linguistics)Stability theoryLogical constantEquations of motionSoftware developerCalculationMathematical analysisINTEGRALLinearizationEigenvalues and eigenvectorsMultiplication signSampling (statistics)Process (computing)NumberPiComputer animation
17:36
Software developerEigenvalues and eigenvectorsComputer simulationFunction (mathematics)Regulärer Ausdruck <Textverarbeitung>Metropolitan area networkTabu searchStorage area networkRight angleCore dumpSign (mathematics)ExpressionMereologyFunctional (mathematics)Computer animation
18:01
Software developerRegulärer Ausdruck <Textverarbeitung>Function (mathematics)Computer simulationFunctional (mathematics)ExpressionWordOrdinary differential equationComputer simulationComputer animation
18:31
Software developerFunction (mathematics)Regulärer Ausdruck <Textverarbeitung>Metropolitan area networkOrdinary differential equationIcosahedronComputer simulationOpen sourceMoment (mathematics)Enterprise architectureComputer animation
19:08
Linear mapAnalog-to-digital converterRotationComputer simulationMetropolitan area networkMaxima and minimaRegulärer Ausdruck <Textverarbeitung>Storage area networkData acquisitionTouchscreenReflection (mathematics)Constraint (mathematics)Slide rulePhysical systemResultantRoundness (object)Spring (hydrology)Forcing (mathematics)LinearizationPosition operator
21:40
SineExecutive information systemMetropolitan area networkRadio-frequency identificationArmStatisticsState of matterNichtlineares GleichungssystemForcing (mathematics)Vector spaceExistenceMechanism design2 (number)Equations of motionComputer animation
22:53
Gamma functionMetropolitan area networkLogical constantCloud computingBoom (sailing)Read-only memoryForceMaxima and minimaComputer simulationArithmetic meanMechanism designComputer animation
23:18
Metropolitan area networkLogical constantCartesian coordinate systemComputer simulationPhysical systemSystem dynamicsDegree (graph theory)Computer animation
23:44
Sample (statistics)Gamma functionCodierung <Programmierung>Logical constantForcePiMetropolitan area networkPlane (geometry)Degrees of freedom (physics and chemistry)Real-time operating systemLevel (video gaming)CodeMoment (mathematics)BitComputer animation
24:11
Metropolitan area networkNewton's law of universal gravitationMaxima and minimaParticle systemComputer simulationBitComputer animation
24:42
Logical constantSample (statistics)Metropolitan area networkKey (cryptography)Gamma functionSet (mathematics)ForceArmMaxima and minimaAreaMultiplication signEigenvalues and eigenvectorsCalculationConstraint (mathematics)Moving averageForcing (mathematics)Computer simulationComputer animation
25:45
Uniformer RaumMaxima and minimaLogical constantMetropolitan area networkGamma functionSample (statistics)ForceKey (cryptography)Variable (mathematics)Data Encryption StandardDegrees of freedom (physics and chemistry)Degree (graph theory)Computer simulationSign (mathematics)Computer animation
26:23
Value-added networkMetropolitan area networkAlpha (investment)Uniformer RaumGraph (mathematics)Function (mathematics)
26:51
Maxima and minimaMetropolitan area networkLogical constantValue-added networkGamma functionArmRegulärer Ausdruck <Textverarbeitung>Computer iconRothe-VerfahrenVotingComputer animation
27:28
Metropolitan area networkLoop (music)Analog-to-digital converterTabu searchStorage area networkComputer simulationAssembly languageSerial portView (database)Object (grammar)Multiplication signBitComputer animation
28:09
Computer simulationMetropolitan area networkStandard deviationInterface (computing)Computer simulationMultiplication signSoftware testingValidity (statistics)Process (computing)Compilation albumJust-in-Time-CompilerComputer animation
28:50
Compilation albumMoment (mathematics)Computer simulationData conversionNumberPhysical systemRight angleMatrix (mathematics)Lecture/Conference
Transcript: English(auto-generated)
00:06
So this is the outline of my talk today. First of all an introduction about multi-body simulation. Next. Yeah, okay some background information. Then I would like to show you some assemblies. So I would fill this talk,
00:24
at least half of it, to show you some assemblies so that you get an impression what is in fact possible with this package in the end. And this package is not in public already. It will be published maybe in
00:44
September. And in the end I will give you a short note of future work and some backup. Package name is multibody-symbolic, which says
01:02
multibody-symbolic. So in fact there are two different ways to approach this problem. You can use symbolic equations and you can stick more to the numerical side. Most of the
01:21
industrial products stick to the numerical side, but we decided to approach it analytically. Our aim is to provide, on basis of existing Python packages, to provide once in a while a complete multi-body simulation tool.
01:43
So why is this important for us? I mean you can guess it. We want to be independent from market leaders. Market leaders, for example, Simpeg or also VL or also Adams, costs a whole lot of money.
02:02
The license and all of these companies are now not, they are now yeah, bought by very big companies, so it is not possible to yeah, to stick in this development process. It's much easier.
02:24
Second, scripting ability is included. Of course, it's in Python, so it's it's of course included. And third, educational purposes. So what is multi-body simulation? Multi-body simulation deals with
02:44
systems which can be described by these equations. These are the Newtonian Euler equations. They are differential equations, and so they are well known since several hundred years. What is the problem? Writing them down and integrating them.
03:01
Well, the problem here is, as you will see in my talk, this F is not only an expression which you can write down. It happens that most of multi-body simulation assemblies include constraints, and constraints are like a ball which runs on a table
03:26
or sliders or whatsoever. So these constraint forces are in fact the difficult thing in multi-body simulations. We look back at the development process of 30 years, more or less, not us, but in general.
03:41
Yeah, the community and the scientific community is working on this problem since more or less 30 years. And you can imagine there come up a whole tons of papers out of there, and yeah, and you cannot expect to climb up within half a year on top.
04:05
But we are very enthusiastic, and we think we can go ahead this way. So use cases of this multi-body simulation is mechanical engineering, ground vehicle dynamics, robotics and biomechanics. Each of these branches are more, each of these branches
04:22
is modern and relevant. So I'm working in ground vehicle dynamics since more or less seven years, and I'm doing ride and handling simulation for another company right now.
04:40
And I know how, what the difficulties are in this field. So, what package are we using? SymPy, it's more or less the basis of this one. SymPy is, these guys who produced SymPy, they were doing a really, really, really great job. And
05:05
I would like to thank you those people. It's a symbolic algebra package, so it replaces more or less these three products, and it's on a very, it's right at the moment. It's on an advanced level, so it's not that you can just do some
05:26
algebra derivatives or integrals. It's on a high advanced level, and it also includes already advanced mechanics. So why do not puzzle these together, stick these loose ends together, and
05:43
produce something with which you can do your multi-body simulation? So fine, we tried to do this. Yeah, I not, I would like not to forget these packages which are really, really helpful and maybe the core of
06:04
numerical and scientific Python programming, NumPy and SymPy. And without these, it would not happen anything about of these. For linear algebraic solvers, we use NumPy. They're very, very, very advanced and
06:24
not to forget, well tested. So these packages are so well tested that you can trust on this. And this is good. There are some ODE solvers in SymPy available. So
06:41
in the end, we need some graphics. Graphics we did with vPyson. vPyson is a, yeah, a medium advanced. It's just really, really nice, but it's, it has some primitives, called primitives, some rods, springs, and so on.
07:03
And you can put these together and make out of these your graphics simulation. In fact, you need this to make sure that your assembled equations behave well. So just a visual check of your
07:23
solution. I mean, you cannot just solve your equations and then put some some graphs and then you cannot judge if this is done really nice or really good. You have some visual
07:41
checkup. Now some background theory. I will try to make this short. Building blocks of a mechanical system. You have your bodies. Everybody has six degrees of freedom. Translation rotational. Including the time derivatives, you end up with twelve degrees of freedom. So each body should
08:03
lead to yeah, to twelve lines in your system of ordinary differential equations, if there wouldn't be some possibilities to boil it down somehow. Each of these bodies has mass and moments of inertia.
08:24
Intrinsic. So this is the figure citation. This is not mine. This is also not my figure, so I will cite it. You can boil down you can boil down the numbers of degrees of freedom by joints.
08:44
Joints you can think about as two bodies are connected. We are joints like in your skeleton. In fact, it's technically done more or less the same. Each of these joints can be produced in a technical sense and
09:01
in the end, with these joints, you reduce the numbers of degrees of freedom. Here we talk about these type of joints, cardionic, axis, revolute. For example, if you think about a tire on a car, it's just revolute.
09:20
Except that, okay, the front tire can also steer, but more or less revolute, so it's a revolute joint. In these equations appears forces and torques. Forces and torques accelerates your masses, so
09:41
forces most of the time appear pairwise, so if one body is attracting another, it's the same done for the other body. But you have different types of forces, pairwise forces. So think about Sun and Earth, for example.
10:01
The external forces, if you just want to find a system of a mechanical system on Earth, you would add the gravitational force, of course, and you wouldn't include the Earth as an independent mass, but you just would add your gravitational field as an kind of external force.
10:26
And the difficult thing is the constraint forces. For example, if you have here a surface and a ball is running on the surface, the surface would act a force on the ball to be on this surface.
10:47
So in multibody simulation, you see many times these kind of drawings, because these kind of drawings show most of the things included. Here, for example,
11:02
you have a chain of bodies and here are joints, and you can see if you switch from Cartesian coordinates to just the angle here, you can reduce it somehow, the number of degrees of freedom. Nobody would write this down in Cartesian coordinates, but in angles. But, you know,
11:24
if the generalized coordinates are minimal, it equals the number of degrees of freedom. But it is not always the case, so you can have more generalized coordinates, so you would have add on top of it some equations, which gives you the constraints.
11:45
Even the generalized coordinates are not unique. So you can, for example, measure the angle towards the set axis or towards the axis of the previous body. So it is not at all unique, the set of generalized coordinates.
12:11
One problem which is really maybe the main or one of the core problems in multibody simulation is this one. So if you have, this is called a constraint loop.
12:25
If you do not just have your joints, but if you have another constraint in the end, which builds a constraint loop, you can say, okay, these angles here on top, they are somehow connected to each other. They cannot be independent at all.
12:44
And this is called constraint loop, and this constraint loop produces another equation, most of the time algebraic equation. And this algebraic equation you can take care of in, there's several possibilities to solve these kind of problems. Most of the time there are
13:07
differential algebraic equation methods, but they are very costly, and in Python, maybe too slow. So we propose here a solution according to Lagrangian with additional drawback force.
13:21
For linearization, you can always put this algebraic equation into your set of equations. Here, this is the set of methods we can use to make, to generate our equations of motions. In this package, they produced here already Cane's method, which is also called
13:44
d'Alembert's principle of virtual work. And Lagrange's method number two, Hamilton's equations is another possibility, but not used here. Okay, so I will show you how it is used. You just need to, if you're a user, you just need to have an object which set up your world, just give you a world coordinate
14:06
system and a marker. A marker is always a coordinate system in the language of multi-body simulation. There are methods provided to add bodies and markers, so extra coordinate systems,
14:22
external forces, extra constraints, and, for example, reflective walls. An ODE solver is connected, more or less automatically. A 3D graphical back-end is connected automatically. So each body signs up in the graphical back-end and appears
14:42
somehow as you want it to be appeared. I will show you examples later on. Some physical quantities are provided, like energy, forces, velocity, and so on. You can, if you, I mean, we are inside Python, so you can calculate whatever you want here.
15:02
And this is the advantage. New and interesting. So what did we put into this work? We would like to have a completeness of joints and tools. And the Jacobian is
15:21
calculated for linearization analysis, so this is on top. Nice feature, very nice feature. The linearization tool, which is already in SymPy, is kind of completed. You can, we can detect automatically independent and dependent coordinates and so on.
15:43
Constraints loop are more or less kind of solved here. External models and parameter span can be included. This is a very important point here, the external models, which I will show you in one example. We have also some B-splines, which can be
16:03
used if you don't have an analytical expression for your force. You can also include some kind of B-splines for having a representative of your force function. So sometimes measurements don't give you analytical expressions. Okay, coding style, set up the system.
16:24
Here, this is your multi-body simulation world. You add your bodies, your markers, your forces, your force models, your external force model, your geometric constraint. Solving, assembling and solving is just this one. Here you have maybe some constant,
16:46
like the gravitational constant. You give it a number, you produce equations of motion and integrate them. Post-processing is calculation of linear analysis,
17:04
linearization, so stability analysis. You prepare and animate your result. Okay, if you would like to be a developer once in a time or if you would like to work with SymPy once in a time, you should be aware of these three things. I mean, you should be aware
17:24
of many things, but these three especially never use SymPy for numerics. Do not try to use, for example, to use SymPy to solve for eigenvectors. It will never happen. It will never be as fast as NumPy. So find the right step between SymPy to NumPy. So I got the sign
17:50
for ten minutes right now, so I will hurry a little bit up. Lambdify is one of the core functions of SymPy. I mean, if you have an algebraic expression, what can you do with it?
18:04
In a computer, you would like to use it as a function, and to use it as a function, you use this word which is called Lambdify. You Lambdify expressions into Python functions, and this is maybe the most impressive invention here. Not our from SymPy. And ODE solvers don't use your own,
18:30
so even if it looks like fun to produce one, it is not never as good as those which are around.
18:41
Use those. We use LSODR solver. Oh yeah, this is also a call for sundials. Maybe if you are interested in these topics, you know that sundials is an open source solver, which is out there but not connected to Python at the moment. It is connected, but it stopped kind of at Python 2.6 or so. So I would like, really, I would like to connect this and use this sundial solver,
19:06
which is really good in my opinion. Okay, I will show you examples, and for these examples, not only the pictures, but also the results, so the movies.
19:35
Okay, this is bad because now it appears here on my screen and not on the screen. I wanted
20:06
it to be appeared. So this is what comes out, more or less. This is a crank slide, something called crank slide, and this is a simple example for a constraint loop.
20:22
Here, this goes round and round, and this goes linear, so this is rotating linear. You see, we put in also some forces, and here an extra constraint force. So this, if we don't put this constraint force here, it would be just a pendulum.
20:52
So I will skip some of these because I have 12, and this is too much, maybe this one.
21:13
Yeah, some spring, and it's a reflective wall here at this position. There's a reflective wall.
21:25
It's interesting because it's not so fluent, but I don't know why, because it's in the reflective wall example. I will show you this one. This is a simple car,
21:48
and this is non-trivial. I mean, maybe you can see this, that assembling these equations takes around 60 seconds. Let's take a look at the assembly here. So these equations,
22:08
they just come up by this mechanism, and this mechanism, as I told you, is quite easy. So we just add bodies, markers, special forces, and so on, and in the end here, this is the steering,
22:24
more or less, and in the end here, this is a starting vector, so you need starting vectors, and in the end, we just call here Canify, and then these equations are assembled.
22:47
And now it's around 60 seconds, assembling equations of motion. And why is this non-trivial? Because it is not only the mechanics, which you can see here, it's also kind of an external model. It's a tire model, which we plugged in,
23:04
and the tire model gives you, I mean, if you think about tires, it's a non-trivial external model, and I'm working in ground vehicle dynamics, so it matters a lot to plug into
23:23
your system non-trivial external models, and the nice thing here is that it is possible to do this. Now it's integrating. We have around 27 degrees of freedom, in fact,
23:50
altogether. No, 28, 28. So it's okay, it's fast. It's not real time at the moment, but it's fast, so it's not bad. I mean, if you would like to make it real time,
24:04
you would have to export it in a C or Fortran code, then it would be real time, but at this level it is not. But it's okay, hope so. Somehow it's a little bit slower,
24:24
I don't know why, because maybe it's some resources here out. Yeah, I can tell you this, we used here a Pacheca model. Pacheca, maybe some of you may know this, Pacheca is a
24:45
professor, he's working on tire models, and he's really, really famous in this area of research, and tires, they are quite interesting to model, because it's rubber, and rubber is always
25:01
difficult. So, for example, if you think about rolling tires or so, you sometimes include a no-slip constraint, but with tires this would never work, because
25:20
tires produce the longitudinal and vertical force only with slip, so you have to include slip in your calculation. So now, here on top of this, we have calculated in the end the eigenvalues of the Jacobian, so you can see this, and now I will go into the animation.
25:49
In this step, all of these degrees of freedom has to be prepared in Cartesian coordinates, that it works. So now here, let's go. This is our car model, it looks simple, but it's okay.
26:06
And it is in fact more or less fixed on zero, so that it, this is a sign steering, it's a simple maneuver, sign steering, but it's nice that it's working well.
26:31
You can see it's breaking in the end here on top, I can show you just this. We have also some outputs calculated for the tires, special values for the tires, and put
26:46
this into these graphs. Okay, so I would skip the examples and go to some future work.
27:18
This is not so easy. We thought about making the assembly equations persistent,
27:38
this is non-trivial because the assemblies are not, they are quite complex objects,
27:48
so ZODB, which I really like, but it's based on serialization, and very complex objects cannot be serialized this way. And I would like to do this to make it persistent, because to skip
28:05
the assembly, which may take a little bit of time. So graphics, always some improvements to be done, model validation, and just-in-time compilation is another nice idea to speed up. Post-processing with Pandas, maybe you know this package. Okay, so future work, complete the basics until September,
28:27
and full vehicle simulation, so to end up with a nice full vehicle simulation, better than that one I showed to you, October to December 2015. And thank you for your attention,
28:41
and I would like to invite you to ask questions. Thank you very much.
29:08
You talked about transitioning from SymPy to NumPy to do numerical computations, is there an automated path from SymPy to NumPy, perhaps?
29:22
At the moment, I'm not aware that they are closely connected, and this is kind of a pitfall. I think you have to, yeah, except for this Lambdify, you need sometimes, if you plug your numbers in, just convert your SymPy matrices to NumPy matrices one by one.
29:46
This is how I did it, and this works fine. Thanks. Okay. Thank you very much.
Recommendations
Series of 2 media