RE: Re: dynamics: comment about MBdyn...

> From: Herman Bruyninckx [mailto:[email protected]]
> Subject: RE: [Orocos] Re: dynamics: comment about MBdyn...
> 
> 
> On Mon, 14 Jul 2003, Tiller, Michael (M.M.) wrote:
> 
> [...]
> > > Well, parsing is not really the big software hurdle... 
> It's defining
> > > the objects that correspond to the models, as well as the tools to
> > > build systems and prepare them for simulation. (And the numeric
> > > solvers too, of course.)
> > 
> > I want to just clarify my comments.  I said "parsing *and* semantic
> > processing tools".  In this context, semantic processing means
> > elaboration of the underlying Modelica code (i.e. things like type
> > checking, name lookup, declaration processing, etc).  I would argue
> > this *is* a big step.  My tools are not complete, but they go well
> > beyond parsing.

> Ok, sounds interesting! But I must confess I do not completely
> understand what you mean exactly. Isn't type checking etc performed
> during parsing?

No not really.  You can think of the process as involving at least three steps:

1) Lexical analysis (i.e. finding the "words")
2) Grammatical analysis (i.e. recognizing the "sentences")
3) Tree walking

#1 and #2 are what I call "parsing" because they involve working with the original source file and converting the textual representation into data structures.  This is true for nearly any parsing tool and any language.

#3 is where much of the actual semantic processing occurs.  During #1 and #2, you really don't know what is a type and what isn't.  That is what is getting sorted out during #1 and #2.  Granted, some people do try and squeeze some basic semantic processing into #2 (e.g. building symbol tables and/or determining scoping), but generally speaking, you need at least two passes for most languages (i.e. #2 and #3, lexical analysis doesn't really count as a "pass").

> I mean, I expect a parser to fail when the types of
> the processed language constructs do not match...
> Probably I am overlooking something here :-)

Syntax errors are detected during phase #2, but all the other stuff (semantic errors) generally doesn't get figured out until #3.

> > Perhaps I'm mistaken, but there seems to be quite a distinct
> > difference in perspectives about what is required here.  You
> > mention "objects that correspond to the models" and "build
> > systems", etc.  The "objects that correspond to the models" are
> > represented in the Modelica code itself.  You seem to be (correct me
> > if I'm wrong) thinking more along the lines of building a C++
> > framework in which you would implement class definitions for nodes,
> > elements, etc and then inherit from these framework components to
> > build specific models.  

> Indeed, that sounds like what I had in mind.

Good.  I wasn't very comfortable about putting words in your mouth.  I just wanted to lay that out (even if I was only speculating about what was in your mind) so that I could contrast the two approaches.

> > I used to be quite deeply in the "framework"
> > camp when I was in grad school (In fact, I built a PDE framework
> > around these same ideas).  Eventually, I abandoned that way of
> > thinking in favor of Modelica.  I'm sure we could argue endlessly
> > about frameworks vs. special purpose languages but let's not.  I'm
> > perfectly happy to accept that for reasons related to your specific
> > goals (e.g. writing embedded control code), you prefer to work with a
> > framework. !

> Well, I do not see the relevance of mentioning a "framework" :-)

When I say framework, I mean that the modeling is done at the C++ (Java, Python) level and that you take advantage of interface definitions, abstract types, etc. at that level to make the work of building the models as minimal as possible.  In other words, you have a framework of supporting code in a certain language that allows you to focus just on the modeling and have the support code take responsibility for the processing stages.

> What I do care about is that one needs runnable code, for simulation
> and control purposes, and having the model description correctly
> parsed is just the small first step, isn't it?

Yes and no.  It is a first step, but it gets you pretty far.  For example, once you have done the complete elaboration that I talked about in the previous message, you have access to the equations and algorithms as data structures.  Then you need to decide what to do with them.  Generating runnable C code is one possibility.  Another possibility is writing an interpreter.  Yet another possibility would be targeting a virtual machine (e.g. for Java or Python).  I'm not saying that elaboration gets you everything, I just wouldn't call it small.

> > I'm not sure what you mean by "build systems".  I can think of two
> > interpretations.  First, you could mean the software build (e.g.
> > autoconf, make, etc).  Or, you could mean the formulation of the
> > system of equations.

> The latter. Including the processing into a form that the DAE solver
> can work with efficiently.

Elaboration of the model gets you a system of hybrid DAEs.  That is it.  If you want to do additional symbolic processing, that is entirely implementation specific.  In other words, the semantics defined in the Modelica specification define rules for transforming a given Modelica model into a set of hybrid DAEs.  That is where the Modelica specification stops and the implementations begin.  The nice thing about having a tool to do the elaboration is not just that it gives you the hybrid DAEs, but it gives them to you in symbolic form so that it is easy to pick up where the semantic processing left off and continue on with any additional processing.  

Believe it or not, getting at the equations isn't of much interest to me (at least at the moment).  My tools are more like "code metric tools".  As such they don't do the complete semantic processing to get at the underlying equations (yet).

> > After semantic processing of a Modelica model is complete, the
> > result is a system of hybrid DAEs.  There are no objects at this
> > point, only variables and equations.  That is it.  The object
> > hierarchy has "washed away" and only a "flattened" representation of
> > the problem remains.

> Ok, I follow you until here. But in order to have a consistent and
> correct flattened representation, quite some model checking and
> connection construction had to take place. And now I begin to see
> where we probably misunderstood each other: what you seem to call
> "parsing", I called "symbolic processing" :-) So, as far as I am
> concerned, we are synchronized again!

Just to make sure we are talking about the same thing, let me give you an example.  A Modelica model for a simple electrical circuit could be described in Modelica as follows:

  model Circuit 
    Modelica.Electrical.Analog.Basic.Ground g;
    Modelica.Electrical.Analog.Basic.Inductor L(L=100e-3);
    Modelica.Electrical.Analog.Basic.Resistor R1(R=15);
    Modelica.Electrical.Analog.Basic.Resistor R2(R=5000);
    Modelica.Electrical.Analog.Basic.Capacitor C(C=100e-6);
    Modelica.Electrical.Analog.Sources.StepVoltage vs(startTime=1);
  equation 
    connect(vs.n, g.p);
    connect(vs.p, L.p);
    connect(R2.p, g.p);
    connect(R2.n, L.n);
    connect(R1.p, R2.n);
    connect(C.n, R1.n);
    connect(C.p, g.p);
  end Circuit;

I think the syntax is obvious enough that I don't need to explain what this circuit is composed of.  Individual components, like ta Resistor include additional detail, e.g.:

  model Resistor "Ideal linear electrical resistor" 
    extends Interfaces.OnePort;
    parameter SI.Resistance R "Resistance";
  equation 
    R*i = v;
  end Resistor;

Note the use of inheritance.

The point is that once you elaborate the model (in Modelica, we call this instantiation but I'm using the more general term "elaboration" because I think it is more intuitive), you get a flattened model.  This is normally available to a given tool as a set of internal data structures organized as a tree.  If you were to formulate the model as "equivalent Modelica code" that consisted only of variable declarations and equations, you would get:

model Circuit
parameter Real L.L = 0.1 "Inductance";
parameter Real R1.R = 15 "Resistance";
parameter Real R2.R = 5000 "Resistance";
parameter Real C.C = 0.0001 "Capacitance";
parameter Real vs.offset = 0 "Voltage offset";
parameter Real vs.startTime = 1 "Time offset";
parameter Integer vs.signalSource.nout = max([size(vs.signalSource.height, 1); 
  size(vs.signalSource.offset, 1); size(vs.signalSource.startTime, 1)]) 
  "Number of outputs";
parameter Integer vs.signalSource.outPort.n = vs.signalSource.nout 
  "Dimension of signal vector";
parameter Real vs.signalSource.offset[:] = {vs.offset} "offset of output signal";
parameter Real vs.signalSource.startTime[:] = {vs.startTime} "output = offset for time < startTime";
parameter Real vs.signalSource.height[:] = {vs.V} "Heights of steps";
parameter Real vs.signalSource.p_height[vs.signalSource.nout] = if size(
  vs.signalSource.height, 1) == 1 then vs.signalSource.height[1]*ones(
  vs.signalSource.nout) else vs.signalSource.height;
parameter Real vs.signalSource.p_offset[vs.signalSource.nout] = if size(
  vs.signalSource.offset, 1) == 1 then vs.signalSource.offset[1]*ones(
  vs.signalSource.nout) else vs.signalSource.offset;
parameter Real vs.signalSource.p_startTime[vs.signalSource.nout] = if size(
  vs.signalSource.startTime, 1) == 1 then vs.signalSource.startTime[1]*ones(
  vs.signalSource.nout) else vs.signalSource.startTime;
parameter Real vs.V = 1 "Height of step";

Real g.p.v "Potential at the pin";
Real g.p.i "Current flowing into the pin";
Real L.v "Voltage drop between the two pins (= p.v - n.v)";
Real L.i "Current flowing from pin p to pin n";
Real L.p.v "Potential at the pin";
Real L.p.i "Current flowing into the pin";
Real L.n.v "Potential at the pin";
Real L.n.i "Current flowing into the pin";
Real R1.v "Voltage drop between the two pins (= p.v - n.v)";
Real R1.i "Current flowing from pin p to pin n";
Real R1.p.v "Potential at the pin";
Real R1.p.i "Current flowing into the pin";
Real R1.n.v "Potential at the pin";
Real R1.n.i "Current flowing into the pin";
Real R2.v "Voltage drop between the two pins (= p.v - n.v)";
Real R2.i "Current flowing from pin p to pin n";
Real R2.p.v "Potential at the pin";
Real R2.p.i "Current flowing into the pin";
Real R2.n.v "Potential at the pin";
Real R2.n.i "Current flowing into the pin";
Real C.v "Voltage drop between the two pins (= p.v - n.v)";
Real C.i "Current flowing from pin p to pin n";
Real C.p.v "Potential at the pin";
Real C.p.i "Current flowing into the pin";
Real C.n.v "Potential at the pin";
Real C.n.i "Current flowing into the pin";
Real vs.v "Voltage drop between the two pins (= p.v - n.v)";
Real vs.i "Current flowing from pin p to pin n";
Real vs.p.v "Potential at the pin";
Real vs.p.i "Current flowing into the pin";
Real vs.n.v "Potential at the pin";
Real vs.n.i "Current flowing into the pin";
Real vs.signalSource.outPort.signal[vs.signalSource.outPort.n] "Real output signals";
Real vs.signalSource.y[vs.signalSource.nout];

equation
g.p.v = 0;
L.v = L.p.v-L.n.v;
0 = L.p.i+L.n.i;
L.i = L.p.i;
L.L*der(L.i) = L.v;
R1.v = R1.p.v-R1.n.v;
0 = R1.p.i+R1.n.i;
R1.i = R1.p.i;
R1.R*R1.i = R1.v;
R2.v = R2.p.v-R2.n.v;
0 = R2.p.i+R2.n.i;
R2.i = R2.p.i;
R2.R*R2.i = R2.v;
C.v = C.p.v-C.n.v;
0 = C.p.i+C.n.i;
C.i = C.p.i;
C.i = C.C*der(C.v);
vs.v = vs.p.v-vs.n.v;
0 = vs.p.i+vs.n.i;
vs.i = vs.p.i;
vs.signalSource.y = vs.signalSource.outPort.signal;
for i in (1:vs.signalSource.nout) loop
  vs.signalSource.outPort.signal[i] = vs.signalSource.p_offset[i]+(if time < 
    vs.signalSource.p_startTime[i] then 0 else vs.signalSource.p_height[i]);
end for;
vs.v = vs.signalSource.outPort.signal[1];
C.n.i+R1.n.i = 0;
R1.n.v = C.n.v;
C.p.i+R2.p.i+g.p.i+vs.n.i = 0;
R2.p.v = C.p.v;
g.p.v = C.p.v;
vs.n.v = C.p.v;
L.n.i+R1.p.i+R2.n.i = 0;
R1.p.v = L.n.v;
R2.n.v = L.n.v;
L.p.i+vs.p.i = 0;
vs.p.v = L.p.v;

end Circuit;

The above was generated automatically by Dymola.  Remember when I said that once you had the data structures you could do several things and C code generation was only one of them.  Well tranformed Modelica output is another (as you can see from the above).  Keep in mind that this listing may seem overly complex for a simple RLC circuit, but much of this comes from the "object-oriented" nature of the formulation.  As a consequence of building these models, each component may have its own local name for something that is identical across several components.  Almost trivial symbolic processing will be able to eliminate these variables.  For example, in the case above, you get the following information from Dymola about the symbolic processing:

STATISTICS
  
Original Model
  Number of components: 8
  Variables: 49
    Constants: 0
    Parameters: 15 (15 scalars)
    Unknowns: 34 (34 scalars)
    Differentiated variables: 2 scalars
  Equations: 30
    Nontrivial : 23
   
Translated Model
  Constants: 6 scalars
  Free parameters: 7 scalars
  Parameter depending: 6 scalars
  Inputs: 0
  Outputs: 0
  Continuous time states: 2 scalars
  Time-varying variables: 6 scalars
  Alias variables: 19 scalars
  Number of mixed real/discrete systems of equations: 0
  Sizes of linear systems of equations: {4}
  Sizes after manipulation of the linear systems: {0}
  Sizes of nonlinear systems of equations: { }
  Sizes after manipulation of the nonlinear systems: { }
  Number of numerical Jacobians: 0

Here are a few interesting observations about the symbolic processing.  First, Dymola recognized that despite the fact that 49 individual variables were declared, there were only 6 unique time varying signals.  Of those 6, two are states (the others are evaluated algebraicly).  Of the 49 variables, 19 were aliases for other variables.  Finally, note that there was a 4x4 linear system of equations present after the BLT transformation but Dymola eliminated it (my guess is that they replaced it with an analytical solution).

An interesting thing to note about the symbolic processing is that it doesn't include any domain specific knowledge about resistors or robots or anything.  It just looks at the Modelica code, applies the semantic rules for constructing the system of hybrid DAEs.

Now before you say "but that doesn't help open source users", keep in mind that I'm not trying to sell you a copy of Dymola.  I'm just trying to provide you with some detailed information about how Modelica tools function.  In fact, I think downloading the demo copy of Dymola would be a great way to at least develop some intuition and understanding about what kinds of symbolic processing is possible.  I realize that you don't want to rely on closed-source tools for your project but *that doesn't mean you cannot learn something from them*.

By the way, some people might be curious to know what the circuit schematic for the example above looks like.  If so, you can download the demo of Dymola and paste the complete model (provided at the bottom of this message) into the "Modelica Source" window or you could look at Figure 3.3 in "Introduction to Physical Modeling with Modelica".

> > > BTW, do you have any idea about the quality of the DAE solvers in
> > > projects like OCtave and Scilab? (They use some form of 
> DASSL, if I am
> > > not mistaken...)
> > 
> > I have not used the DAE solvers in Octave and Scilab.  The DAE
> > solvers I am familiar with are also based on DASSL but since they
> > exist in closed source tools, I'm sure they are of no interest to
> > you.
> Unless you have an idea about what they changed or added to the DASSL
> code...

I'm not sure that they changed much in DASSL (I really don't know), but I can say that they probably benefit enormously from the symbolic pre-processing (especially index reduction) that can be done *BEFORE* you pass the problem off to a numerical solver.

> Herman

--
Mike

P.S. - Here is the complete Modelica source for my example including graphical annotations:

model Circuit 
  Modelica.Electrical.Analog.Basic.Ground g
    annotation (extent=[-10, -40; 10, -20]);
  Modelica.Electrical.Analog.Basic.Inductor L(L=100e-3)
    annotation (extent=[-40, 50; -20, 70]);
  Modelica.Electrical.Analog.Basic.Resistor R1(R=15)
    annotation (extent=[20, 50; 40, 70]);
  Modelica.Electrical.Analog.Basic.Resistor R2(R=5000)
    annotation (extent=[-10, 10; 10, 30], rotation=90);
  Modelica.Electrical.Analog.Basic.Capacitor C(C=100e-6)
    annotation (extent=[50, 10; 70, 30], rotation=90);
  Modelica.Electrical.Analog.Sources.StepVoltage vs(startTime=1)
    annotation (extent=[-40, -30; -20, -10]);
  annotation (Diagram);
equation 
  connect(vs.n, g.p) annotation (points=[-20, -20; 0, -20], style(color=3));
  connect(vs.p, L.p)
    annotation (points=[-40, -20; -60, -20; -60, 60; -40, 60], style(color=3));
  connect(R2.p, g.p)
    annotation (points=[-6.12303e-016, 10; 0, 10; 0, -20], style(color=3));
  connect(R2.n, L.n) annotation (points=[6.12303e-016, 30; 0, 30; 0, 60; -20, 60], style(color=3));
  connect(R1.p, R2.n) annotation (points=[20, 60; 6.12303e-016, 60; 6.12303e-016, 30], style(color=3));
  connect(C.n, R1.n)
    annotation (points=[60, 30; 60, 60; 40, 60], style(color=3));
  connect(C.p, g.p)
    annotation (points=[60, 10; 60, -20; 0, -20], style(color=3));
end Circuit;