The Need for Formulae Manipulation in the Object-oriented Modeling of Physical Systems

1
The Need for Formulae Manipulationin the
Object-oriented Modelingof Physical Systems

• François E. Cellier, Ph.D.
• Institute of Computational Science
• ETH Zurich
• Switzerland
• July 26, 2006

2
Object-oriented Modeling

• Todays modeling needs extend to physical systems
  that are described by many thousands of
  mathematical equations.
• Such models are being created and maintained by
  employing the principles of object-oriented
  graphical modeling.
• Each object contains a number of different layers
  (aspects, realities) that are to be described
  shortly.
• Symbolic formula manipulation plays a central
  role in converting these object-oriented
  descriptions to a form that can be simulated
  effectively and efficiently.

3
Object Descriptions

• Objects may contain an iconic representation that
  is being used to embed the object at the next
  higher hierarchical level.
• Objects may contain a geometric representation
  that is being used for animation.
• Objects may contain a topological description
  that shows its decomposition structure.
• Objects may contain a mathematical description
  that shows the equations that the model is formed
  of.
• Objects may contain a parametric description that
  shows parameters and their default values.
• Objects may contain a documentation layer that
  describes the object verbally.

4
Heterogeneous Modeling Formalisms
5
Heterogeneous Modeling Formalisms II
6
A First Example
7
A First Example II
Component equations U0 f(t) U0 v1
v0 uR1 R1 iR1 uR1 v1 v2 uR2 R2 iR2
uR2 v2 v0 iC C1 duC/dt uC v2 v0
uL L1 diL/dt uL v1 v0 v0 0 Node
equations i0 iR1 iL iR1 iR2 iC
The circuit contains 5 components and 3 nodes. ?
We require 13 equations in 13 unknowns.
8
The Structure Incidence Matrix
?
9
The Structure Digraph
?
10
The Tarjan Algorithm

• The Tarjan Algorithm is based on the structure
  digraph.
• It is a graph coloring algorithm.

? equations with only a single black line leading
to them, color that line in red and color all
black lines that emanate from the indicated
unknown in blue. Renumber the equations anew
starting with 1. ? unknowns with only a single
black line leading to them, color that line in
red and color all black lines that emanate from
the indicated equation in blue. Renumber the
equations anew starting with n, where n denotes
the number of equations.
11
The Tarjan Algorithm An Example
?
12
The Tarjan Algorithm An Example II
?
13
The Tarjan Algorithm An Example III
?
14
The Tarjan Algorithm An Example IV
?
15
The Structure Incidence Matrix II
?
16
Algebraic Loops
Component equations U0 f(t) u3 R3 i3 u1
R1 i1 uL L diL/dt u2 R2 i2 Node
equations i0 i1 iL i1 i2 i3 Mesh
equations U0 u1 u3 uL u1 u2 u3 u2
The circuit contains 5 components ?
We require 10 equations in 10 unknowns
17
Algebraic Loops An Example
?
18
Algebraic Loops An Example II
?
19
Algebraic Loops An Example III
?
20
The Tearing of Algebraic Loops

• The following heuristic may be used to find
  suitable tearing variables

Determine those equations in the digraph that
contain the largest number of black lines. For
each of these equations, find the unknown with
the largest number of black lines emanating from
it. Determine for each of these variables, how
many additional equations can be causalized, when
this variable is assumed known. Choose that
variable as the next tearing variable that
causalizes the largest number of additional
equations.
21
Algebraic Loops An Example IV
?
22
Algebraic Loops An Example V
?
23
Algebraic Loops An Example VI
?
06
07
24
Algebraic Loops An Example VII
?
06
04
07
25
Algebraic Loops An Example VIII
?
26
Algebraic Loops An Example IX
?
27
The Structure Incidence Matrix III
?
28
The Solving of Algebraic Loops I

• The Tarjan algorithm identifies and isolates
  algebraic loops.
• It transforms the structure incidence matrix to
  BLT form, whereby the diagonal blocks are made as
  small as possible.
• The selection of the tearing variables is not
  done in a truly optimal fashion. This is not
  meaningful, because the optimal selection of
  tearing variables has been shown to be an
  np-complete problem. Instead, a set of
  heuristics is being used, which usually comes up
  with a small number of tearing variables,
  although the number may not be truly minimal.
• The Tarjan algorithm does not concern itself with
  how the resulting algebraic loops are being
  solved.

29
The Solving of Algebraic Loops II
Structure diagram
30
Solution of Algebraic Loops III
?
i1 i2 i3 u2 / R2 u3 / R3 u3 /
R2 u3 / R3 ((R2 R3 ) / (R2 R3 ))
u3 ((R2 R3 ) / (R2 R3 )) (U0 - u1 )
((R2 R3 ) / (R2 R3 )) (U0 - R1 i1 )
?
31
The Solving of Algebraic Loops IV

• The algebraic loops can be solved either
  analytically or numerically.
• If the loop equations are non-linear, a Newton
  iteration on the tearing variables may be
  optimal.
• If the loop equations are linear and if the set
  is fairly large, Newton iteration may still be
  the method of choice.
• If the loop equations are linear and if the set
  is of modest size, the equations can either be
  solved by matrix techniques or by means of
  explicit formulae manipulation.

32
Structural Singularities An Example I
We compose a model using the currents, the
Voltages and the potentials. Consequently, the
mesh equations are being ignored. We have 7
circuit components plus the ground, therefore 2?7
1 15 equations. In addition, there are 4
nodes giving rise to another 3
equations. Therefore, we expect 18 equations in
18 unknowns.
For passive components, it is customary to
normalize the Voltages in the same direction as
the currents. For active components (sources),
the reverse is true.
33
Structural Singularities An Example II
?
34
Structural Singularities An Example III
?
35
Structural Singularities An Example IV
?
36
The Algorithm by Pantelides

• As soon as a constraint equation has been found,
  this equation must be differentiated.
• In the algorithm of Pantelides, the
  differentiated constraint equation is added to
  the set of equations.
• Consequently, the set of equations has now one
  equation too many.
• In order to re-equalize the number of equations
  and unknowns, one of the integrators that is
  associated with the constraint equation is being
  eliminated.

37
The Algorithm by Pantelides II
An additional unknown was created by the
elimination of the integrator. x and dx are now
algebraic variables, for which there must be
found equations.
38
The Algorithm by Pantelides III

• When differentiating constraint equations, it can
  happen that additional new variables are being
  created, e.g. v ? dv, where v is an algebraic
  variable.
• Since v was already blue (otherwise, this would
  not have been a constraint equation), there
  already exists another equation to compute v.
• This equation must also be differentiated.
• The differentiation of additional equations
  continues until no additional variables are being
  created.

39
The Algorithm by Pantelides An Example
1 I1 f1(t) 2 I2 f2(t) 3 I3 f3(t) 4
uR R iR 5 uL1 L1 diL1 /dt 6 uL2 L2
diL2 /dt 7 iC C duC /dt 8 v0 0
9 u1 v0 v1 10 u2 v3 v2 11 u3 v0
v1 12 uR v3 v0 13 uL1 v2 v0 14 uL2
v1 v3 15 uC v1 v2
16 iC iL1 I2 17 iR iL2 I2 18 I1 iC
iL2 I3 0
40
The Algorithm by Pantelides An Example II
1 I1 f1(t) 2 I2 f2(t) 3 I3 f3(t) 4
uR R iR 5 uL1 L1 diL1 /dt 6 uL2 L2
diL2 7 iC C duC /dt 8 v0 0
9 u1 v0 v1 10 u2 v3 v2 11 u3 v0
v1 12 uR v3 v0 13 uL1 v2 v0 14 uL2
v1 v3 15 uC v1 v2
?
16 iC iL1 I2 17 iR iL2 I2 18 I1 iC
iL2 I3 0 19 dI1 diC diL2 dI3 0
41
The Algorithm by Pantelides An Example III
1 I1 f1(t) 2 I2 f2(t) 3 I3 f3(t) 4
uR R iR 5 uL1 L1 diL1 /dt 6 uL2 L2
diL2 7 iC C duC /dt 8 v0 0
9 u1 v0 v1 10 u2 v3 v2 11 u3 v0
v1 12 uR v3 v0 13 uL1 v2 v0 14 uL2
v1 v3 15 uC v1 v2
16 iC iL1 I2 17 iR iL2 I2 18 I1 iC
iL2 I3 0 19 dI1 diC diL2 dI3 0
42
The Algorithm by Pantelides An Example IV
43
Inline Integration

• Higher-order models are almost invariably stiff.
• Stiff models require an implicit ODE solver for
  their simulation.
• Implicit ODE solvers call for an additional
  iteration during each integration step.
• Inline integration may be able to reduce the
  number of iteration variables by merging the
  solver equations with the model equations.

44
Inline Integration An Example
?
iL
11
45
Inline Integration An Example II
01
11
08
09
?
iL
10
11
46
Inline Integration An Example III
01
11
08
09
?
02
iL
10
11
Choice
47
Inline Integration An Example IV
01
11
03
08
09
?
05
02
04
06
07
iL
10
11
Choice
48
Inline Integration II

• DASSL would have required to use all state
  variables (1) plus all tearing variables (1) as
  iteration variables.
• With inline integration, we got away with a
  single iteration variable in this example.
• Determining a good (small) set of tearing
  variables is thus essential, as it ultimately
  determines the efficiency of the simulation.
• Dymola offers a highly sophisticated selection of
  heuristics to determine small sets of tearing
  variables. These heuristics have not been
  published, as the company views them as their
  technological edge.

49
Discontinuity Handling

• Discontinuity handling requires special
  considerations during the (symbolic) translation.
  (Switching equations change their computational
  causality as a function of the switch position.)
• Discontinuity handling also requires special
  considerations during the (numerical) simulation.
  (Extrapolation polynomials dont exhibit
  discontinuities.)

50
Mixed Symbolic and Numerical Processing