Simulating a model, one stage at a time¶
The design rule for this project is that every stage is a value you can hold. That is hard to believe from a prose description and obvious from a notebook, so this one walks a model from source text to trajectory and stops at each stage to look at what came out.
The model is examples/models/RLCCircuit.mo — a resistor, an inductor and a capacitor,
built from the component library in examples/library/Electrical.mo.
from pathlib import Path
root = next(p for p in [Path.cwd(), *Path.cwd().parents] if (p / "examples" / "models").is_dir())
source = (root / "examples" / "models" / "RLCCircuit.mo").read_text()
library = (root / "examples" / "library" / "Electrical.mo").read_text()
print(source)
model RLCCircuit "A series RLC circuit stepped with 1 V -- damped oscillation, two states"
Electrical.ConstantVoltage source(V = 1) annotation(Placement(transformation(extent = {{-18, -64}, {18, -28}})));
Electrical.Resistor r(R = 1) annotation(Placement(transformation(extent = {{-80, 16}, {-44, 52}})));
Electrical.Inductor l(L = 1, i(start = 0, fixed = true)) annotation(Placement(transformation(extent = {{-18, 16}, {18, 52}})));
Electrical.Capacitor c(C = 1, v(start = 0, fixed = true)) annotation(Placement(transformation(extent = {{44, 16}, {80, 52}})));
Electrical.Ground ground annotation(Placement(transformation(extent = {{3, -114}, {39, -78}})));
equation
connect(source.p, r.p);
connect(r.n, l.p);
connect(l.n, c.p);
connect(c.n, source.n);
connect(source.n, ground.p);
end RLCCircuit;
1. Source → AST¶
The parser is hand-written recursive descent. It returns modelica.ir values, and it
refuses constructs it does not implement by name rather than misreading them.
from modelica.lang import parse
program = parse(source)
[c.name for c in program.classes]
['RLCCircuit']
2. AST → flat model¶
Flattening instantiates the components, resolves the connect equations into the Kirchhoff
laws they stand for, and produces a pile of scalar equations that happen to be true. No
direction of computation yet — Modelica's = is acausal, so nothing here says which
variable an equation solves for.
from modelica.sim import flatten
flat = flatten([parse(library), program], "RLCCircuit")
print(f"{len(flat.variables)} variables, {len(flat.equations)} equations")
print(f"balanced: {flat.is_balanced()}")
30 variables, 26 equations balanced: True
3. Flat model → causalized system¶
This is the stage that does the work: alias elimination, index reduction where the model needs it, matching each equation to the variable it computes, and sorting the result into blocks that can be evaluated in order. Whatever is left unmatched becomes a state.
from modelica.sim import causalize
system = causalize(flat)
print(f"states: {system.states}")
print(f"candidates: {system.candidates}")
print(f"index: {system.index}")
print(f"blocks: {len(system.blocks)}, largest {system.largest_block}")
print(f"explicit: {system.is_explicit}")
states: ('l.i', 'c.v')
candidates: (('l.i', 'c.v'),)
index: 1
blocks: 4, largest 1
explicit: True
candidates is every state set that would also have worked. Here there is only one, so
there is nothing to choose. The pendulum is the interesting case — see the
verification notebook.
4. Causalized system → generated code¶
The IR is compilable by construction, which is what lets a Python engine stand next to the
Fortran solvers. emit writes the model out as straight-line arithmetic; nothing walks a
tree of closures at evaluation time.
from modelica.sim import emit
code = emit(system)
print(code.source[: code.source.index("def rhs")])
import numpy as np
_seed = np.eye(2)
_null = np.zeros(2)
def parameters(p):
source_V = 1
r_R = 1
l_L = 1
c_C = 1
p[0] = source_V
p[1] = r_R
p[2] = l_L
p[3] = c_C
return p
def initial(p, y, w):
source_V = p[0]
r_R = p[1]
l_L = p[2]
c_C = p[3]
y[0] = 0
y[1] = 0
w[0] = source_V
w[1] = 0.0
w[2] = r_R
w[3] = 0.0
w[4] = 0
w[5] = l_L
w[6] = 0
w[7] = c_C
w[8] = 0.0
w[9] = 0.0
w[10] = 0.0
w[11] = 0.0
w[12] = 0.0
w[13] = 0.0
w[14] = 0.0
w[15] = 0.0
w[16] = 0.0
w[17] = 0.0
w[18] = 0.0
w[19] = 0.0
w[20] = 0.0
w[21] = 0.0
w[22] = 0.0
w[23] = 0.0
w[24] = 0.0
w[25] = 0.0
w[26] = 0.0
w[27] = 0.0
w[28] = 0.0
w[29] = 0.0
w[30] = 0.0
w[31] = 0.0
return y
def latch(t, y, p, w):
return w
Here is the right-hand side itself — the function the integrator calls a few times per step, and the reason the shape of the IR matters:
body = code.source[code.source.index("def rhs") :]
print(body[: body.index("\n\n\n")])
def rhs(t, y, p, w, dy):
source_V = p[0]
r_R = p[1]
l_L = p[2]
c_C = p[3]
l_i = y[0]
c_v = y[1]
r_v = r_R * l_i
l_v = -r_v + source_V - (c_v - (source_V - source_V) / -1.0)
der_l_i = -(l_v / -l_L)
der_c_v = -(l_i / -c_C)
dy[0] = der_l_i
dy[1] = der_c_v
return dy
5. Simulate¶
Radau by default: a model that came out of a Modelica flattener is stiff more often than
not, and an implicit method that is merely slower on an easy problem beats an explicit one
that fails on a hard one.
from modelica.sim import simulate
result = simulate(system, stop=0.05, points=500)
result
The chart is a self-contained HTML page with no CDN behind it: every variable the model
has is in the legend, the moving ones are ticked on to begin with, and events are drawn as
vertical marks. It writes itself to a file with result.plot().save("rlc.html").
An event, and why the loop restarts¶
BouncingBall is the one to look at for discontinuities. The integrator carries an
implicit assumption that what it is integrating is smooth over the step it is taking; a
bounce violates that exactly once per bounce. So the solver is stopped at the crossing, the
discontinuity is applied, and integration starts again from the other side.
ball = causalize(
flatten([parse((root / "examples" / "models" / "BouncingBall.mo").read_text())], "BouncingBall")
)
bounces = simulate(ball, stop=3.0, points=800)
for event in bounces.events:
print(f"t = {event.time:.6f} {event.cause}")
t = 0.451524 when h <= 0 t = 1.173961 when h <= 0 t = 1.751912 when h <= 0 t = 2.214272 when h <= 0 t = 2.584160 when h <= 0 t = 2.880071 when h <= 0
bounces.plot()