Enumerating the discrete half¶
A simulation engine is two machines glued together. One is continuous — a stiff integrator
with an error norm and a convergence theorem behind it. The other is discrete: latches
flipping, when branches firing, a state set being reselected, an FMU walking its calling
sequence. Only the first has a literature of error control.
The second is finite, so it can be enumerated instead of sampled. That is what
modelica.verify does, using musil as the
checker. The design note is Verification; this notebook runs it.
from pathlib import Path
from modelica.lang import parse
from modelica.sim import System, causalize, flatten, simulate
root = next(p for p in [Path.cwd(), *Path.cwd().parents] if (p / "examples" / "models").is_dir())
def system(name: str, *libraries: str) -> System:
"""Parse, flatten and causalize one of the example models."""
sources = [(root / "examples" / "library" / f"{lib}.mo").read_text() for lib in libraries]
sources.append((root / "examples" / "models" / f"{name}.mo").read_text())
return causalize(flatten([parse(text) for text in sources], name))
1. Event iteration¶
At one instant of simulated time the continuous state does not move, so the only thing that
can change is which watched relations are currently latched true. There are two-to-the-few
of those, which means the event iteration at an instant is an exactly finite machine —
no abstraction, the emitted indicators and transition are the functions being called.
Two questions get asked. Does it stop? A cycle of events that never settles is a hung simulation. Does the order matter? The engine applies the first stale event it finds, and nothing chooses that order on purpose — so if two orders reach different fixed points, that is a bug waiting for a model that notices.
from modelica.verify import check_iteration
tank = system("DrainingTank")
report = check_iteration(tank)
print(report)
DrainingTank: event iteration settles, and settles in one place, from all 2 configurations (3 discrete states)
The tank has one watched relation, so there are two configurations to start from — and a simulation only ever visits one of them. The bug this catches is the one hiding in the other.
from modelica.verify import iteration
for start in iteration(tank).initial_states():
print(start)
latches=0 fired=[] latches=1 fired=[]
Checking at t = 0 is thin for a model whose only event is a when: nothing is pending at
the start, so nothing has to settle. The instants worth asking about are the ones a run
actually reached, so check_run simulates first and then checks every event instant.
from modelica.verify import check_run
ball = system("BouncingBall")
for report in check_run(ball, simulate(ball, stop=2.0, points=400)):
print(report)
BouncingBall: event iteration settles, and settles in one place, from all 1 configurations (1 discrete states) BouncingBall: event iteration settles, and settles in one place, from all 1 configurations (1 discrete states) BouncingBall: event iteration settles, and settles in one place, from all 1 configurations (1 discrete states) BouncingBall: event iteration settles, and settles in one place, from all 1 configurations (1 discrete states)
2. Dynamic state selection¶
The pendulum is the classic index-2 model: x² + y² = L² constrains the states, so index
reduction has to differentiate it, and the selection it lands on goes singular partway
through the swing. Structural analysis cannot see that coming — it is a statement about
numbers, not about incidence.
What it can do is enumerate the alternatives:
pendulum = system("Pendulum")
print(f"chosen: {pendulum.states}")
for candidate in pendulum.candidates:
print(f" {candidate}")
chosen: ('y', 'vy')
('y', 'vy')
('x', 'vx')
(x, vx) is the one the literature names, and it is regular exactly where (y, vy) is
not. Note it is two exchanges away from the chosen set, which is why the candidates are
enumerated rather than found by perturbing the selection one variable at a time.
Switching between them is not written yet. Designing it as a checked automaton first means modelling the linear algebra as an adversary: which set is usable is not the engine's to choose, and the engine only ever learns a set has gone singular by failing to factor it.
from modelica.verify import check_selection
selection = check_selection(pendulum)
print(selection)
for assumption in selection.safety.unverified_assumptions:
print(f"\n {assumption.name}\n {assumption.description}")
Pendulum: switching between 2 state sets always recovers a regular one, assuming index-reduction-leaves-a-usable-set, jacobians-hold-still-during-a-search
index-reduction-leaves-a-usable-set
at every point of the trajectory at least one admissible state set has a non-singular constraint Jacobian. Structural analysis cannot establish this -- it is a statement about numbers, not about incidence
jacobians-hold-still-during-a-search
the Jacobians change only as the step advances, not while the engine is part-way through trying candidates. Drop it (`quiescent=False`) and liveness fails, not safety: a set can recover *after* the engine ruled it out, so the engine exhausts every candidate at a point where one of them was usable all along, and no step is ever accepted again
Those assumptions are not decoration. The second one says the Jacobians hold still while the engine is part-way through trying candidates — drop it and the check fails, with the run that breaks it:
broken = check_selection(pendulum, quiescent=False)
print(f"safety still ok: {broken.safety.ok}")
print(f"liveness ok: {broken.liveness.ok} ({broken.liveness.kind})")
print()
print(broken.liveness)
safety still ok: True
liveness ok: False (p-unreachable)
LIVENESS VIOLATED: 'a usable state set' is unreachable from a reachable state
Init: chosen=(y,vy) singular={} rejected={}
→ (y,vy) goes singular: chosen=(y,vy) singular={(y,vy)} rejected={}
→ detect: chosen=(y,vy) singular={(y,vy)} rejected={(y,vy)}
→ switch to (x,vx): chosen=(x,vx) singular={(y,vy)} rejected={(y,vy)}
→ (y,vy) recovers: chosen=(x,vx) singular={} rejected={(y,vy)}
→ (x,vx) goes singular: chosen=(x,vx) singular={(x,vx)} rejected={(y,vy)}
→ detect: chosen=(x,vx) singular={(x,vx)} rejected={(x,vx), (y,vy)}
(24 states explored)
A set recovered after the engine had already ruled it out, so the engine exhausted every candidate at a point where one of them was usable all along.
3. The FMI 3.0 calling sequence¶
FMI does not describe its interface as a list of functions. It describes a state machine,
and says which fmi3* calls are legal in each state. Most of what goes wrong between two
tools that both "support FMI" is a sequence error: every individual call is fine, and the
order is not.
The machine is written down once, before the FMI layer that has to obey it.
from modelica.fmi import Interface, Lifecycle, ModeError
from modelica.verify import check_modes
for interface in Interface:
print(check_modes(interface))
Model Exchange: all 8 modes reachable, no dead ends, and an instance can always still be freed Co-Simulation: all 9 modes reachable, no dead ends, and an instance can always still be freed Scheduled Execution: all 8 modes reachable, no dead ends, and an instance can always still be freed
"Every mode reachable" is not a formality — it caught Event Mode sitting in the Scheduled Execution table, where the standard has Clock Activation Mode instead.
The same table drives the runtime guard. This is what an exported FMU turns into
fmi3Error:
life = Lifecycle(Interface.MODEL_EXCHANGE)
life.call("fmi3EnterInitializationMode")
life.call("fmi3ExitInitializationMode")
life.call("fmi3EnterContinuousTimeMode")
print(f"mode: {life.mode}")
try:
life.call("fmi3UpdateDiscreteStates")
except ModeError as error:
print(f"\n{error}")
mode: Continuous-Time Mode
fmi3UpdateDiscreteStates is not allowed in Continuous-Time Mode (Model Exchange); after fmi3EnterInitializationMode -> fmi3ExitInitializationMode -> fmi3EnterContinuousTimeMode. Allowed here: fmi3CompletedIntegratorStep, fmi3DeserializeFMUState, fmi3EnterEventMode, fmi3FreeFMUState, fmi3FreeInstance, fmi3GetAdjointDerivative, fmi3GetContinuousStateDerivatives, fmi3GetContinuousStates, fmi3GetDirectionalDerivative, fmi3GetEventIndicators, fmi3GetFMUState, fmi3GetNominalsOfContinuousStates, fmi3GetNumberOfVariableDependencies, fmi3GetVariableDependencies, fmi3Get{VariableType}, fmi3Reset, fmi3SerializeFMUState, fmi3SerializedFMUStateSize, fmi3SetContinuousStates, fmi3SetDebugLogging, fmi3SetFMUState, fmi3SetTime, fmi3Terminate
4. The co-simulation master¶
An SSP system is several FMUs stepped together, which makes the master a distributed algorithm wearing a numerical hat. Everybody has to end each communication step at the same time, any component may refuse the step it was given, and undoing a step someone else already took requires that they kept a copy of themselves first.
FMI is explicit that the copy is optional: canGetAndSetFMUState is a capability flag, and
an FMU is allowed to say no. With every component able to roll back, the master is fine:
from modelica.ssp import Component
from modelica.verify import check_master
print(check_master((Component("a", rollback=True), Component("b", rollback=True))))
master over a, b: every component reaches the horizon together
Now make one of them a perfectly conforming FMU that declines to save its state. Rollback is a precondition of a rollback-based algorithm, so the master checks it at load time and refuses — which is the correct outcome, not a failure:
MIXED = (Component("a", rollback=True), Component("b", rollback=False))
mixed = check_master(MIXED)
print(f"ok: {mixed.ok}, declined: {mixed.declined}")
print(mixed)
ok: True, declined: True master over a, b (no rollback): declined before stepping anything -- a rollback-based algorithm needs `canGetAndSetFMUState` from every component
That precondition is not decoration either. Drop it and the trap comes back:
print(check_master(MIXED, require_rollback=False).safety)
INVARIANT VIOLATED: a-rollback-can-put-everyone-back Init: save at 0 times=(0,0) saved=(-,-) → save a: save at 0 times=(0,0) saved=(0,-) → begin: step at 0 times=(0,0) saved=(0,-) → step b: step at 0 times=(0,1) saved=(0,-) → a discards: rollback at 0 times=(0,1) saved=(0,-) (13 states explored)
Four moves from the start. b takes its step, a refuses, and there is no move that puts
b back — the master has nothing left to offer. No test that happens to use FMUs which all
support state saving will ever see it, which is why the precondition is checked rather than
assumed.
What this cannot do¶
The limits matter as much as the claims. None of this says anything about the continuous solve — accuracy, stiffness, step-size control, Jacobian quality are out of scope, and pretending otherwise would be worse than not checking at all. It is also the wrong tool for matching, BLT sorting and index reduction, which are graph algorithms over unbounded inputs rather than finite automata; those want property-based testing.
Model checking earns its place where the state is small and the orderings are the problem. That is these four and nothing else.