Re: dynamics: comment about MBdyn...

Herman Bruyninckx wrote:
> Some weeks ago, there were some posts about the free software dynamics
> package MBDyn <http://www.aero.polimi.it/~mbdyn/index.html>, and I
> took a closer look at it... Since there is apparently no mailinglist
> for MBDyn,

apparently, setting up a mailing list at our University requires
enough traffic to justify it, ... a sort of catch 22.  The address
mailto:[email protected] is read by all the developers, and we
may decide to set a mailing list up in the near future.

> and since the topic might be of interest to some Orocos
> people, I post my comments here, but address the MBDyn developers in
> the second person :-)
> 
> I have read with much interest the documentation that is available on
> the MBdyn webpage, and I am impressed by the functionality you have
> already. Of course, I do have lots of questions and comments... :-)

which, as usual, are quite welcome :)

> 
> First of all, I am personally quite a bit biased by the Bond Graph
> approach, both because I have many good friends that work in this
> area and so I know it better than other domains, and because I really
> think it is the best fit for an "object-oriented" software engineering
> solution to general purpose dynamics software.
> 
> Probably the most visible effect of this bias is that I am looking for
> separation of software modules between 
> - structure and dynamical properties (i.e., modelling the
>   interconnectivy separately from the dynamics attached to nodes in
>   the connectivity graph(s)).
> - what is application-dependent (e.g., the fact that MBdyn and Orocos
>   work basically in the mechanical domain) and application-independent
>   things (such as the tools to model and build connected dynamical
>   systems, the needs for symbolic pre-processing of DAE systems, the
>   GUI, etc.).
> 
> My basic remark is that you have a very tight "coupling" between
> stuff that should be decoupled. 

Yes.  When MBDyn started, I had very little programming
experience, so all I felt was some vague need for separation
and uncoupling.  Trust me: earlier versions of MBDyn
were much more coupled.  Moreover, whem it started, nobody
even vaguely realized it could become what it's now, and
even less we knew it could become an open source project.
As I mentioned in a recent paper, it is there now and it is
working, but there is a lot to do to ease its open development.

> For example, modelling, simulation,
> and file IO; topology and properties; mathematical representations and
> models; particular mathematical library dependencies; etc.
> It is my strong believe that this leads to unmaintainable projects,
> unless maintained only by the original developers.
> (Our own robot control software was in the same condition, so I know,
> and I also know that _that_ is the most important reason I tried
> to start up Orocos :-)
> 
> Here follow more detailed technical things. (The order is quite
> random.) Some of them are real questions, to which I think you might
> give an answer; some others are more open issues listed here for
> further thought and discussion :-)
> 
> About the "input" manual.
> That PDF document renders very slowly because you don't use PDF
> fonts, but bitmapped fonts. It's better to use "pdflatex" or "dvips
> -Ppdf" for production of the PDF file.)
> When trying to (pdf)latex the TeX files myself, I encountered the
> error that the file "version.tex" is not found. Another error: you use
> "time_step" instead of "time\_step".
> "sintax"-> "syntax"
> Line 615 of data.tex needs an ending "}" for the \emph.
> The reference MASARATI-LANZ-MANTEGAZZA-2001 is not included.

I think you already noted this; doc sources are there because
I wanted to distribute everything that's available; however,
they are not part of the build process because latex is not
popular enough to make everything build seamlessly.  You should
get docs from the website.  Of course I'll try to make
distributed doc sources compile and build as well as I can.

The reference is

P. Masarati, M. Lanz, and P. Mantegazza
``Multistep Integration of Ordinary, Stiff and Differential-Algebraic 
Problems for Multibody Dynamics Applications''
presented at the XVIo Congresso Nazionale AIDAA, 24-28 Settembre 2001, 
Palermo, Italy.

http://mbdyn.aero.polimi.it/~masarati/Publications/MULTISTEP_INTEGRATION-AIDAA16.ps.gz

> 
> About discussion of symbolic needs. 
> - Where exactly do you need symbolic support? I mean, what do you use
>   it for, and what do you _want_ to use it for in the future?
>   (I guess for the reduction of sets of DAEs to computationally
>    deterministic smaller sets?)

No.  At present, symbolic manipulation is in a very early
stage, and all the support that's there is in "interpreted"
drives (see the discussion about drives that occurs later),
and, recently, in a symbolic constitutive law (that requires
GiNaC; the rest is built-in) for symbolic differentiation.
In MBDyn's present state there is no low-level symbolic
manipulation; modelling occurs by using existing building
blocks in existing element libraries.  One can "easily" add
new elements with limited effort by hacking the code, or by
building dynamic modules that are linked run-time.

> - Are these symbolic needs decomposable in smaller modules, that can
>   be developed and used independently?

The short answer is: I guess they are.

This is in a very preliminary stage; we haven't decided
yet how far we want to go with symbolic manipulation.
In our approach, there is very limited need for equation
differentiation, because we directly integrate index 3
DAE problems without index reduction, and we are moving
towards matrix-free solvers which do not require accurate
jacobians.  On the other hand, using symbolic manipulation
to __reduce__ the equations before moving to numeric is
something very complicated, which, in my opinion, can be
of use when playing with mostly linear networks of lumped
elements, but it is much less effective when dealing with
deformable bodies if the ratio of algebraic constraints
to structural degrees of freedom is low.  In this case,
structural dofs cannot be eliminated because of dynamics.

> 
> About File IO.
> - Structure of models and their properties?

Can you detail your question?  What I understand is "how is
a MBDyn model structured?", but it is very broad.  If this
is what you mean, please have a coke/beer and take some spare
time to read the answer ;)

A model is: analysis data + control data + nodes + elements
(partial) modularity is obtained by using include directives,
much like in C.  Details about the above mentioned entities
will be answered later in dedicated questions.

MBDyn solves Initial Value Problems (IVP), so it needs
a model, initial values, and an integrator.  Initial
values are closely related to the model, so there is very
little separation (if any); the integrator has very little
to do with the model, so it is specified in a different
part of the input file.  A model is made of public degrees
of freedom, the nodes, which are connected by elements;
elements can have private degrees of freedom, so in some
sense they are supersets of nodes.  Nodes provide dofs,
elements write equations.

* analysis data: represents an attempt to separate what's
related to a specific simulation from what's the model.
Basically, it contains initial and final time, time step,
integration method parameters, nonlinear solver parameters,
linear solver parameters, diagnostics output requests, and more.

* control data: contains the count of entities of each type,
mostly for historical reasons, and info on how the model
consistency should be checked and repaired if needed.
That is, elements that are used in structural analysis
and require algebraic constraints to be satisfied from
the beginning, participate in a pseudo-dynamic solution
that tries to preserve IV while satisfying algebraic
constraints.  By deciding IVP priorities one in principle
can use this procedure to bring together separate parts
or change system configuration before the real analysis
starts.  There are still some problems in some cases,
and it's not worth spending too much time.

* nodes: they are in a separate block because in MBDyn's
early logic all the public dofs had to be known in
advance (and, in its early implementation, the count
of all the entities had to be known in advance because
of "optimal" memory allocation and addressing.  I know
this sounds silly, now; this explains the need to count
entities in the control data).

* elements: they're in the final block.  One of the
requirements is that dependencies must be fixed as soon
as they are read.  So, if an element depends on another,
they must be input in the appropriate sequence.  This
need was avoided with respect to nodes by moving them in
a separate block which is read first.

There are other secondary data blocks that are not
mentioned to keep this description simple.


> - Useful to transform to XML + Schema?

Wy not? We've been talking for a long time about
how to improve the input.  At some point, before
it derived towards a complete reinvetion of the
wheel (say OO language) we were also considering
the possibility to interpret (at least partially)
Modelica.

> - Reuse of shared objects and properties with Orocos?

Why not? In this sense, it would be useful to do
a sort of review of what's needed by MBDyn, what's
available and what could be rewritten in terms of
functional decoupling, and then to check what's
already available in OROCOS, and what OROCOS could
benefit from MBDyn.

> - Useful to separate files in functionally decoupled files:
>   - physical modelling.
>   - data
>   - integrator properties
>   - simulation properties (initial values, timing, logging, ...)

Partially, it's already there.  It's not so clearly
separated, but, for instance, the IVP solution
is basically unaware of what's the model, and the model
is basically unaware of what's the integration scheme
and the solution pattern; the model, to access the current
time, needs a specific drive!

> - What exactly is "control data"? 

see above

> ("The control data block mostly
>   contains information about the problem that is required to ensure
>   that a consistent model will be generated."
>   - aren't all generated models consistent by construction? Is this
>     possible/desirable?

In principle, I'm talking about two different types
of consistency (I might have been unclear with regard
to this part).  First, "consistency" in a loose way
means the model is not what the user intended it to be,
say one wants to build a model that requires three
joints, but only two are defined.  In this case, the
item count in the "control data" should help detecting
there's a problem. Second, "consistency" in terms of
algebraic constraints means that, for different reasons,
one could find easier to input a system which does
not satisfy the algebraic constraints to the desired
accuracy.  In this case, the "control data" block
contains info about how to attempt to restore constraint
consistency.


> - "Nodes": generic frameworks such as Bond Graphs and Modelica offer
>   already well-developed terminology and definitions, so it could be
>   interesting to follow those. Personally, I am very interested in
>   discussing this in detail.

I was not aware of bond graphs at the time I started
working at MBDyn, and I currently have a vague idea
of what they are.  I'm more biased towards finite
elements; in this sense we have "state" and "flow"
fields, but while in bond graphs a "force" is typically
an "effort", in finite elements it would be in the "flow"
category.  MBDyn's nodes are "efforts", and MBDyn's
elements mostly write "flow" equations in terms of
forces.  I note that in bond graphs causality seems
to play an important role, while Modelica is based
on the negation of causality (the same is true for
MBDyn), so the examples you mention, at least on this
point, are contradictory ...

> - "Elements"...? Seems to me that nodes describe the "connectivity"
>   of a dynamic system ("bonds in Bond Graph parlance) and elements
>   describe the "constitutive relationships" of dynamic objects (=
>   energy storing elements, and energy tranforming elements).

Elements describe both connectivity and constitutive
relationships.

> - "Drivers": these represent two concepts, I think:
>   - "sources" of motion or force.

They are __used__ to "represent sources of motion or
force", but they're more.  They're values that may
depend on time (remember? IVP) or some other
context-dependent scalar variable.

>   - particular "setpoint generation algorithms" (linear, cubic, ...).
>   I think this coupling is not good.

Not really; they might have been used for this purpose
somewhere, but it is not what they're intended for.

Every time a constant can be proficiently turned into
a parameter that depends on time or something else,
a drive is used.  Drives are: the magnitude of forces
and couples; the magnitude of initial strains in elastic
elements; the magnitude of trim tab deflections in simple
aerodynamic elements, and more.

A drive can be:
* a certain function of one variable (zero, one, const,
time, linear, parabolic, sine, cosine, exp, ...)
* a value extracted from a node (a dof, some elaboration
of nodal values, parameters indirecty extracted from
elements, ...);
* a string containing a math expression that is interpreted
at every call;
...

The variable usually is time, but in some cases it
can be anything else; they can be added, or, in some
cases, called recursively, ...

> - you only seem to use "constitutive relationships" for elastic
>   and viscous elements; I think also the inertial properties are
>   "constitutive relationships", so should be modelled by the same
>   software objects. (Which leads straightforwardly to Bond Graphs :-)

I agree in principle; however, in our current
implementation we decided to get rid of
accelerations, so we cannot use it the way
I guess you're thinking.

> 
> Peter Soetens is about to finish a very powerful (both in
> efficiency and in scope) library for "property management", i.e.,
> the setting and getting of properties of objects in a complex software
> syste, with interfaces to XML parsing etc. I guess it should be
> directly useable in MBdyn too.

Definitely a good suggestion.

> your file syntax is very much position dependent (see the simple
> example of the free body on p. 5), which is very error-prone; 

I agree.  I didn't try to verify if there are "dangling"
cases (e.g. possible multiple interpretations of data
sequences), but, since there's no grammar specification,
all of this is very flurry.  In my experience, every time
an input card is erroneously written, the error is caught
somewehre (usually very close to the point where it occurs).

One improvement I was considering is to group structured
data into structure symbols (e.g. vectors in {}, matrices
in {{},{},{}} and so).

 > I would definitely suggest to add the name of the property to the value,
> as would be normal in an XML-based model. 

Agree.  Currently, I have the problem of preserving
backwards compatibility with existing models, but
I understand MBDynòs input syntax is outdated and
should be improved.

> It should also be
> possible/required to add the physical units, I guess.

This is another big improvement I considered; however,
the field we're familiar with is again that of FEM, where
physical units have been introduced only at the level
of GUI, while codes require users to be consistent.
This is why I never added any support for that (and it's
true that it could remain bonded to the GUI without any
loss of generality).

> 
> About some code.
> - I browsed through it for an hour or two, but it is not so easy to
>   grasp what every file is supposed to encode. And the reason is that
>   I have no "large scale" framework (Modelica, Bond Graphs) to attach
>   the terminology and the objects to. I really think this is vital for
>   the survival of a complex piece of code such as MBdyn: using
>   terminology which is very well documented in these two big
>   frameworks...

...

> - most of the comments in the files are in Italian. Not so much of a
>   problem for me, but it might be difficult for most people :-)

As it is stated in the manual, this IS GOOD, because people
should not trust comments :); when I first realized the code
could have an international interest, I turned all the
diagnostics and started naming variables in English; I didn't
consider translating comments, however.

> - I guess that the code in which you use ASSERTs and THROWs is 
>   only used during setup or configuration? (Otherwise you will get
>   some efficiency problems in your RTAI integration, I guess.)

Yes, you need to explicitly enable it by means of the
compile switch -DDEBUG (--enable-debug in configure),
as opposed to standard assert() which must be explicitly
disabled.

> - do you have some documentation about all the math library "wrapping"
>   you have done? This would be highly relevant information for Orocos
>   (and many other projects), because we have to do it to :-)

Not yet.  I don't think this wrapping is optimal, but,
of course, you are free to use it; before you start using
it, I suggest you browse boost (we might move in that
direction, at least by adding it as an option;
I don't like to require people to download too many pieces
of software to be able to compile mine :) and, in any case,
we need to clean it up a bit: it's not been reviewed for
long time.  I can provide some documentation in a reasonable
time if you're interested in it.

> - if I am not mistaken, you use operator overloading for matrix
>   multiplications? Also not real-time safe, in principle. (But the
>   places where I have seen it used are not real-time parts of the
>   code.)

I don't see any performance limitation in operator
overloading.  Apart from the fact that only matrix
summation overloading is currently used, all the
operations are optimized as much as possible;
overloading of Mat3x3, Vec3 and similar is mandatory
to allow readable coding of fundamental physics of
structural components.  In some cases this is
creating a lot of temporaries, but this is part of
optimization (e.g. large operations unrolling, or
use of special tools that automatically do that).

> - what is a "DriveHingeJoint"?

a spherical joint where the rotation vector
is prescribed (by means of a 3D drive ...)

> - what are the differences between a "ElasticHingeJoint" and a
>   "DeformableHingeJoint"?

A DeformableHingeJoint is a deformable joint whose
internal torque depends on the relative configuration;
an ElasticHingeJoint is a DeformableHingeJoint whose
internal torque depends on the relative rotation only.
The internal torque of a ViscoElasticHingeJoint
depends also on the relative angular velocity.

>   (My BG bias tells me to decouple these things in a separare ideal
>    joint and a separare elastic element :-) )

Agree.

> - I see files that depend on Ginac in the same directory as
>   non-dependent files; I think it would be very worthwhile to separate
>   the symbolic functionalities from the numeric.

I think you overlooked the importance of GiNaC in the
current implementation.  GiNaC is used for an elastic
constitutive law that requires symbolic differentiation
(and that's experimental).  Moreover, since the
constitutive laws are C++ templates, I found it easier
to group them in the header file that must be included
by those files that require the template definition.

>   Similar remarks for dependencies on sockets and file IO.

Anyway, I agree with your remark, although I was trying
to avoid too much directory proliferation.

> - I couldn't guess what the TplDrive stuff is all about...

It is an attempt to have a drive (see above) that applies
to an arbitrary dimensionality by means of templates.
It is required by the constitutive laws because they
need prestrains and more with different dimensionality
(double for scalars, Vec3 for deformable joints, Vec6
for beams).

> - you have a "struct" directory; what does "struct" stand for?
>   (I guess a/the finite-dimensional subset of the mechanical domain?)

Yes.  I include beam elements because they're discrete
models, though.

> - in the "elec" directory I had expected to see electrical object
>   models, but I encounter stuff like "DynamicForgettingFactor" and
>   "GPCInv"...

Leftovers from my Ph.D.; in principle, all electric
stuff will go there, and currently a few more sensors
and the electric motor are there in the new release;
the discrete control stuff, on the contrary, should
move into "black box" stuff.

> 
> 
> 
> About PDE modelling and simulation.
> What we miss (well, nobody really _misses_ it at this moment :-)
> in the Orocos community are the "inifinite dimensional"
> dynamical objects, (with the Timoshenko beam as the classical
> example), which requires PDEs instead of ODEs. So, it would be very
> useful to brainstorm about how to incude both in the same software
> structure.

I'm puzzled by the use that people tend to do of PDE
acronym  when talking about structural analysis.
It is clear to me that at some point you need to move
to finite dimensional models, e.g. discretize, and I
don't find it too much useful to preserve PDE description
of the problem until the very last moment because it is
well accepted that efficient and accurate FEM models do
not allow PDE formulation and discretization separation,
unfortunately.  As a consequence I'm much more comfortable
talking directly about a beam element rather than
Timoshenko PDE beams.

>   
> About the MBdyn Tutorial.
> - Free rigid body
>   Why separate the inertial stuff from the momentum stuff? THese
>   things are both inherent properties of a mechanical body, so I would
>   model and store them together.
>   Or, do you mean that the "degrees of freedom owner" is a
>  "junction" in Bond Graph terms, and that you can attach several
>  inertial components ("I" objects in BGs) to the same junction?

Exactly.  In principle, in a node you may have momentum
and momenta moment coming from multiple lumped masses and
from multiple beams connecting to it.  Currently beams
do not have distributed inertia because of historical
reasons, however we can add it at any moment, since the
formulation is fully developed and tested (in another
software, though).

In principle, momentum is associated with the dof, while
inertia is a property of those elements that are connected
to the node.

> 
> - the "body" example on p. 5 seems to suggest a strong coupling
>   between the specification of the dynamics and the topology.

I'm not sure I understand your point here.  The body
must be attached to a node; the "natural" way is to
rely on the element's position and orientation for
the lumped inertial properties, but you can use different
reference frames for each body property.  The only
constraint is that the inertia properties are constant,
and rigidly attached to the node (e.g. the offset
and the relative orientation do not change).

Please do not hesitate in pointing out other problems,
I'm open to any kind of suggestions, and I'll try to
learn as much as possible from OROCOS' lessons.  As
I mentioned may times in the past, MBDyn's structure
and design is the result of a long and not always
linear development, which has been problem driven
for many aspects and thus may not always reflect a clear
and general view.  This explains most of the couplings
that hurt your sensibility.  We're currently working
on version 1.2 with a lot of rewrking in many parts.
There is no radical architectural change, because
we are in the middle of experimenting many new modules
(new nonlinear solvers, new time integration algorithms,
new specialized elements) and we are trying to test
this with a minimal impact on the overall architecture
(you, and other interested developers, can have
a tarball for preview; it's not being released yet
because most of this work is still unpublished, and
some features are not working, or not well integrated
yet; I expect a release this fall or by the end
of the year).  Then, of course, we're thinking about
MBDyn 2, with new architecture, to exploit some of the
features that are missing and that you highlighted
so accurately, to make the investigation of new
features, as the ones I mentioned above, much easier.

Pierangelo.

-- 
Dr. Pierangelo Masarati               | voice: +39 02 2399 8309
mailto:[email protected]  | fax:   +39 02 2399 8334
http://www.aero.polimi.it/~masarati
Dip. Ing. Aerospaziale Politecnico di Milano,
via La Masa 34, 20156 Milano, Italy