The shape of the equations¶

Matching and BLT sorting are described in the literature as graph algorithms, and implemented here as ones — but what they do is permute a matrix until it is lower triangular. That is a picture, so this notebook draws it.

Three stages get looked at: the incidence matrix a flat model comes out with, what alias elimination deletes from it, and the staircase that causalization leaves behind. The last cell is a dropdown over every example model in the repository, and it needs a live kernel. The file itself is attached to this page — structure.ipynb — but it reads those models from examples/, so it wants a checkout rather than a bare download:

git clone https://gitlab.com/jorgeecardona/pymodelica && cd pymodelica
pip install "modelica[sim]" matplotlib ipywidgets
jupyter lab docs/notebooks/structure.ipynb

In [1]:

from pathlib import Path

import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.colors import ListedColormap

# The models come from `examples/` in the repository rather than from strings in here: a copy
# pasted into a notebook is a copy that goes stale, and the reason this page is executed when the
# site is built is that it then cannot show output from a model the project does not have.
root = next(
    (p for p in [Path.cwd(), *Path.cwd().parents] if (p / "examples" / "models").is_dir()), None
)
if root is None:
    message = (
        "This notebook reads its models from `examples/` in the pymodelica repository. "
        "Run it from a checkout:\n"
        "    git clone https://gitlab.com/jorgeecardona/pymodelica\n"
        "    cd pymodelica && jupyter lab docs/notebooks/"
    )
    raise RuntimeError(message)

SURFACE, INK, MUTED, LINE = "#fcfcfb", "#0b0b0b", "#52514e", "#d4d3cd"
BLUE, ORANGE, AQUA = "#2a78d6", "#eb6834", "#1baf7a"

mpl.rcParams.update(
    {
        "figure.figsize": (8.4, 3.6),
        "figure.dpi": 110,
        "figure.facecolor": SURFACE,
        "axes.facecolor": SURFACE,
        "savefig.facecolor": SURFACE,
        "axes.edgecolor": LINE,
        "axes.labelcolor": MUTED,
        "axes.spines.top": False,
        "axes.spines.right": False,
        "axes.grid": False,
        "axes.titlelocation": "left",
        "axes.titlepad": 10,
        "axes.titlesize": 11,
        "lines.linewidth": 2.0,
        "legend.frameon": False,
        "text.color": INK,
        "xtick.color": MUTED,
        "ytick.color": MUTED,
        "font.size": 10,
    }
)
%config InlineBackend.figure_formats = ["svg"]

from pathlib import Path import matplotlib as mpl import matplotlib.pyplot as plt import numpy as np from matplotlib.colors import ListedColormap # The models come from `examples/` in the repository rather than from strings in here: a copy # pasted into a notebook is a copy that goes stale, and the reason this page is executed when the # site is built is that it then cannot show output from a model the project does not have. root = next( (p for p in [Path.cwd(), *Path.cwd().parents] if (p / "examples" / "models").is_dir()), None ) if root is None: message = ( "This notebook reads its models from `examples/` in the pymodelica repository. " "Run it from a checkout:\n" " git clone https://gitlab.com/jorgeecardona/pymodelica\n" " cd pymodelica && jupyter lab docs/notebooks/" ) raise RuntimeError(message) SURFACE, INK, MUTED, LINE = "#fcfcfb", "#0b0b0b", "#52514e", "#d4d3cd" BLUE, ORANGE, AQUA = "#2a78d6", "#eb6834", "#1baf7a" mpl.rcParams.update( { "figure.figsize": (8.4, 3.6), "figure.dpi": 110, "figure.facecolor": SURFACE, "axes.facecolor": SURFACE, "savefig.facecolor": SURFACE, "axes.edgecolor": LINE, "axes.labelcolor": MUTED, "axes.spines.top": False, "axes.spines.right": False, "axes.grid": False, "axes.titlelocation": "left", "axes.titlepad": 10, "axes.titlesize": 11, "lines.linewidth": 2.0, "legend.frameon": False, "text.color": INK, "xtick.color": MUTED, "ytick.color": MUTED, "font.size": 10, } ) %config InlineBackend.figure_formats = ["svg"]

A model big enough to have a shape¶

The example models are small on purpose, and a 5×5 matrix has no texture to look at. So this one is generated: an RC ladder of n stages, written as Modelica text and parsed like any other source. The component library it draws on is the same examples/library/Electrical.mo the hand-written examples use.

In [2]:

from modelica.lang import parse, unparse_equation

electrical = parse((root / "examples" / "library" / "Electrical.mo").read_text())


def ladder(stages: int) -> str:
    """Modelica source for a resistor–capacitor ladder of `stages` stages."""
    head = [f'model Ladder "An RC ladder of {stages} stages"']
    head += ["  Electrical.ConstantVoltage source(V = 1);", "  Electrical.Ground ground;"]
    for i in range(stages):
        head += [
            f"  Electrical.Resistor r{i}(R = 1);",
            f"  Electrical.Capacitor c{i}(C = 1, v(start = 0, fixed = true));",
        ]
    body = ["equation", "  connect(source.p, r0.p);"]
    for i in range(stages):
        body += [f"  connect(r{i}.n, c{i}.p);", f"  connect(c{i}.n, source.n);"]
        if i + 1 < stages:
            body += [f"  connect(c{i}.p, r{i + 1}.p);"]
    body += ["  connect(source.n, ground.p);", "end Ladder;"]
    return "\n".join([*head, *body]) + "\n"


print(ladder(2))

from modelica.lang import parse, unparse_equation electrical = parse((root / "examples" / "library" / "Electrical.mo").read_text()) def ladder(stages: int) -> str: """Modelica source for a resistor–capacitor ladder of `stages` stages.""" head = [f'model Ladder "An RC ladder of {stages} stages"'] head += [" Electrical.ConstantVoltage source(V = 1);", " Electrical.Ground ground;"] for i in range(stages): head += [ f" Electrical.Resistor r{i}(R = 1);", f" Electrical.Capacitor c{i}(C = 1, v(start = 0, fixed = true));", ] body = ["equation", " connect(source.p, r0.p);"] for i in range(stages): body += [f" connect(r{i}.n, c{i}.p);", f" connect(c{i}.n, source.n);"] if i + 1 < stages: body += [f" connect(c{i}.p, r{i + 1}.p);"] body += [" connect(source.n, ground.p);", "end Ladder;"] return "\n".join([*head, *body]) + "\n" print(ladder(2))

model Ladder "An RC ladder of 2 stages"
  Electrical.ConstantVoltage source(V = 1);
  Electrical.Ground ground;
  Electrical.Resistor r0(R = 1);
  Electrical.Capacitor c0(C = 1, v(start = 0, fixed = true));
  Electrical.Resistor r1(R = 1);
  Electrical.Capacitor c1(C = 1, v(start = 0, fixed = true));
equation
  connect(source.p, r0.p);
  connect(r0.n, c0.p);
  connect(c0.n, source.n);
  connect(c0.p, r1.p);
  connect(r1.n, c1.p);
  connect(c1.n, source.n);
  connect(source.n, ground.p);
end Ladder;

1. What flattening hands over¶

The incidence matrix has one row per equation and one column per unknown, and a mark where the equation mentions the unknown. Nothing in it says which equation solves for which unknown — Modelica's = is acausal, and this matrix is exactly that ignorance written down.

In [3]:

from modelica.ir.rewrite import names
from modelica.sim import System, causalize, flatten


def flat_of(stages: int):
    return flatten([electrical, parse(ladder(stages))], "Ladder")


flat = flat_of(6)
unknowns = [v.name for v in flat.unknowns]
column = {name: i for i, name in enumerate(unknowns)}

incidence = np.zeros((len(flat.equations), len(unknowns)), dtype=int)
for row, equation in enumerate(flat.equations):
    for name in {*names(equation.lhs), *names(equation.rhs)}:
        if name in column:
            incidence[row, column[name]] = 1

density = incidence.sum() / incidence.size
print(f"{incidence.shape[0]} equations × {incidence.shape[1]} unknowns, {density:.1%} filled")

from modelica.ir.rewrite import names from modelica.sim import System, causalize, flatten def flat_of(stages: int): return flatten([electrical, parse(ladder(stages))], "Ladder") flat = flat_of(6) unknowns = [v.name for v in flat.unknowns] column = {name: i for i, name in enumerate(unknowns)} incidence = np.zeros((len(flat.equations), len(unknowns)), dtype=int) for row, equation in enumerate(flat.equations): for name in {*names(equation.lhs), *names(equation.rhs)}: if name in column: incidence[row, column[name]] = 1 density = incidence.sum() / incidence.size print(f"{incidence.shape[0]} equations × {incidence.shape[1]} unknowns, {density:.1%} filled")

80 equations × 80 unknowns, 2.8% filled

Drawn, with the marks in blue and the empty entries left as surface:

In [4]:

CELLS = ListedColormap([SURFACE, BLUE])


def spy(ax, matrix, title: str, cmap=CELLS) -> None:
    """One binary matrix, equations down, unknowns across."""
    ax.imshow(matrix, cmap=cmap, interpolation="nearest", aspect="equal", vmin=0, vmax=cmap.N - 1)
    ax.set_title(title)
    ax.set_xlabel("unknowns →")
    ax.set_ylabel("← equations")
    ax.set_xticks([])
    ax.set_yticks([])
    for side in ("top", "right"):
        ax.spines[side].set_visible(True)
    for spine in ax.spines.values():
        spine.set_color(LINE)


fig, ax = plt.subplots(figsize=(4.6, 4.6))
spy(ax, incidence, f"Ladder(6), as flattened — {incidence.sum()} marks")
plt.show()

CELLS = ListedColormap([SURFACE, BLUE]) def spy(ax, matrix, title: str, cmap=CELLS) -> None: """One binary matrix, equations down, unknowns across.""" ax.imshow(matrix, cmap=cmap, interpolation="nearest", aspect="equal", vmin=0, vmax=cmap.N - 1) ax.set_title(title) ax.set_xlabel("unknowns →") ax.set_ylabel("← equations") ax.set_xticks([]) ax.set_yticks([]) for side in ("top", "right"): ax.spines[side].set_visible(True) for spine in ax.spines.values(): spine.set_color(LINE) fig, ax = plt.subplots(figsize=(4.6, 4.6)) spy(ax, incidence, f"Ladder(6), as flattened — {incidence.sum()} marks") plt.show()

2. Most of it is a = b¶

The band down the diagonal and the scatter around it are two different things. A connect(a, b) between potential variables becomes a.v = b.v, and a large fraction of what a flattener emits is exactly that: equations whose only content is that two names are the same variable. Alias elimination spends them, and the system keeps a map from every name it removed to what replaced it — so a user can still ask for r3.p.v afterwards.

In orange, the columns whose variable did not survive:

In [5]:

from modelica.ir import nodes as ir

system = causalize(flat)
survives = np.array([name not in system.aliases for name in unknowns])

kept = np.zeros_like(incidence)
kept[(incidence == 1) & survives] = 1
kept[(incidence == 1) & ~survives] = 2

fig, ax = plt.subplots(figsize=(4.6, 4.6))
spy(
    ax,
    kept,
    "Blue survives causalization, orange is aliased away",
    ListedColormap([SURFACE, BLUE, ORANGE]),
)
plt.show()

trivial = sum(1 for e in flat.equations if isinstance(e.lhs, ir.Ref) and isinstance(e.rhs, ir.Ref))
computed = len(unknowns) - len(system.aliases) - len(system.states)
print(f"{trivial} of {len(flat.equations)} equations are literally `a = b`\n")
print(f"{len(unknowns)} unknowns in the flat model:")
print(f"  {len(system.aliases):3} eliminated as aliases")
print(f"  {len(system.states):3} became states — integrated, not computed")
print(f"  {computed:3} are computed by the {len(system.blocks)} blocks")
print(
    f"\n(the blocks also compute {len(system.derivatives)} derivative variables the flat model never had)"
)

from modelica.ir import nodes as ir system = causalize(flat) survives = np.array([name not in system.aliases for name in unknowns]) kept = np.zeros_like(incidence) kept[(incidence == 1) & survives] = 1 kept[(incidence == 1) & ~survives] = 2 fig, ax = plt.subplots(figsize=(4.6, 4.6)) spy( ax, kept, "Blue survives causalization, orange is aliased away", ListedColormap([SURFACE, BLUE, ORANGE]), ) plt.show() trivial = sum(1 for e in flat.equations if isinstance(e.lhs, ir.Ref) and isinstance(e.rhs, ir.Ref)) computed = len(unknowns) - len(system.aliases) - len(system.states) print(f"{trivial} of {len(flat.equations)} equations are literally `a = b`\n") print(f"{len(unknowns)} unknowns in the flat model:") print(f" {len(system.aliases):3} eliminated as aliases") print(f" {len(system.states):3} became states — integrated, not computed") print(f" {computed:3} are computed by the {len(system.blocks)} blocks") print( f"\n(the blocks also compute {len(system.derivatives)} derivative variables the flat model never had)" )

33 of 80 equations are literally `a = b`

80 unknowns in the flat model:
   57 eliminated as aliases
    6 became states — integrated, not computed
   17 are computed by the 23 blocks

(the blocks also compute 6 derivative variables the flat model never had)

3. The staircase¶

What causalization returns is a sequence of blocks: each one is a set of equations and the same number of unknowns, sorted so that a block only depends on blocks before it. Drawn with the blocks in evaluation order, that ordering is the lower-triangular shape — every mark above the diagonal would be a dependency on something not computed yet.

Marks on the diagonal are the unknown a block solves for; marks below it are the values it reads. The diagonal is orange to separate the two.

In [6]:

def staircase(system: System) -> np.ndarray:
    """Block-ordered incidence: 2 on the diagonal a block solves, 1 where it reads."""
    order = [u for block in system.blocks for u in block.unknowns]
    at = {name: i for i, name in enumerate(order)}
    matrix = np.zeros((len(order), len(order)), dtype=int)
    row = 0
    for block in system.blocks:
        for equation in block.equations:
            for name in {*names(equation.lhs), *names(equation.rhs)}:
                if name in at:
                    matrix[row, at[name]] = 1
            for name in block.unknowns:
                matrix[row, at[name]] = 2
            row += 1
    return matrix


sorted_blocks = staircase(system)

fig, ax = plt.subplots(figsize=(4.6, 4.6))
spy(
    ax,
    sorted_blocks,
    f"Ladder(6) causalized — {len(system.blocks)} blocks",
    ListedColormap([SURFACE, BLUE, ORANGE]),
)
plt.show()

print(
    f"largest block: {system.largest_block}   explicit: {system.is_explicit}   states: {system.states}"
)

def staircase(system: System) -> np.ndarray: """Block-ordered incidence: 2 on the diagonal a block solves, 1 where it reads.""" order = [u for block in system.blocks for u in block.unknowns] at = {name: i for i, name in enumerate(order)} matrix = np.zeros((len(order), len(order)), dtype=int) row = 0 for block in system.blocks: for equation in block.equations: for name in {*names(equation.lhs), *names(equation.rhs)}: if name in at: matrix[row, at[name]] = 1 for name in block.unknowns: matrix[row, at[name]] = 2 row += 1 return matrix sorted_blocks = staircase(system) fig, ax = plt.subplots(figsize=(4.6, 4.6)) spy( ax, sorted_blocks, f"Ladder(6) causalized — {len(system.blocks)} blocks", ListedColormap([SURFACE, BLUE, ORANGE]), ) plt.show() print( f"largest block: {system.largest_block} explicit: {system.is_explicit} states: {system.states}" )

largest block: 1   explicit: True   states: ('c0.v', 'c1.v', 'c2.v', 'c3.v', 'c4.v', 'c5.v')

is_explicit is the claim that matters: every block is 1×1 and carries an assignment, so evaluating the right-hand side is straight-line arithmetic with no inner solve anywhere. The generated code for this model has no while in it.

How the three counts scale¶

The interesting number is not how many equations flattening produces — it is how few survive. For the ladder, the equations grow at twelve per stage and the blocks at four, because the rest were connect statements:

In [7]:

rows = []
for n in range(1, 11):
    model = flat_of(n)
    causalized = causalize(model)
    rows.append((n, len(model.equations), len(causalized.blocks), len(causalized.states)))

stages, equations, blocks, states = (np.array(values) for values in zip(*rows, strict=True))

fig, ax = plt.subplots(figsize=(6.2, 3.4))
for series, colour, label in (
    (equations, BLUE, "flat equations"),
    (blocks, ORANGE, "blocks after sorting"),
    (states, AQUA, "states"),
):
    ax.plot(stages, series, color=colour, marker="o", markersize=5)
    ax.annotate(
        label,
        xy=(stages[-1], series[-1]),
        xytext=(6, 0),
        textcoords="offset points",
        color=colour,
        fontsize=9,
        va="center",
    )

ax.set_xlim(0.6, 12.8)
ax.set_xticks(range(1, 11))
ax.set_title("An RC ladder, stage by stage")
ax.set_xlabel("stages")
ax.set_ylabel("count")
ax.grid(visible=True, color="#e8e7e2", linewidth=0.8)
ax.set_axisbelow(True)
plt.show()

rows = [] for n in range(1, 11): model = flat_of(n) causalized = causalize(model) rows.append((n, len(model.equations), len(causalized.blocks), len(causalized.states))) stages, equations, blocks, states = (np.array(values) for values in zip(*rows, strict=True)) fig, ax = plt.subplots(figsize=(6.2, 3.4)) for series, colour, label in ( (equations, BLUE, "flat equations"), (blocks, ORANGE, "blocks after sorting"), (states, AQUA, "states"), ): ax.plot(stages, series, color=colour, marker="o", markersize=5) ax.annotate( label, xy=(stages[-1], series[-1]), xytext=(6, 0), textcoords="offset points", color=colour, fontsize=9, va="center", ) ax.set_xlim(0.6, 12.8) ax.set_xticks(range(1, 11)) ax.set_title("An RC ladder, stage by stage") ax.set_xlabel("stages") ax.set_ylabel("count") ax.grid(visible=True, color="#e8e7e2", linewidth=0.8) ax.set_axisbelow(True) plt.show()

When the staircase cannot be built: the pendulum¶

Everything above sorts completely because every block came out 1×1. The pendulum does not. In Cartesian coordinates its constraint x² + y² = L² involves no derivative at all, so no amount of permuting produces a triangular form — the matrix is structurally singular until index reduction differentiates the constraint and adds the result to the system.

In [8]:

pendulum = causalize(
    flatten([parse((root / "examples" / "models" / "Pendulum.mo").read_text())], "Pendulum")
)

print(
    f"index: {pendulum.index}   blocks: {len(pendulum.blocks)}   largest: {pendulum.largest_block}   explicit: {pendulum.is_explicit}"
)
for block in pendulum.blocks:
    kind = "assignment" if block.is_linear_scalar else f"{len(block)}×{len(block)} solve"
    print(f"\n  {kind}: {', '.join(block.unknowns)}")
    for equation in block.equations:
        print(f"    {unparse_equation(equation)}")

pendulum = causalize( flatten([parse((root / "examples" / "models" / "Pendulum.mo").read_text())], "Pendulum") ) print( f"index: {pendulum.index} blocks: {len(pendulum.blocks)} largest: {pendulum.largest_block} explicit: {pendulum.is_explicit}" ) for block in pendulum.blocks: kind = "assignment" if block.is_linear_scalar else f"{len(block)}×{len(block)} solve" print(f"\n {kind}: {', '.join(block.unknowns)}") for equation in block.equations: print(f" {unparse_equation(equation)}")

index: 3   blocks: 3   largest: 3   explicit: False

  1×1 solve: x
    ['x ^ 2 + y ^ 2 = L ^ 2;']

  assignment: vx
    ['2.0 * x * vx + 2.0 * y * vy = 0.0;']

  3×3 solve: F, der_vy, der_vx
    ['m * der_vx = -F * x / L;']
    ['m * der_vy = -F * y / L - m * g;']
    ['2.0 * vx * vx + 2.0 * x * der_vx + (2.0 * vy * vy + 2.0 * y * der_vy) = 0.0;']

Read that from the top. The first block solves x out of x² + y² = L² — 1×1, and still no assignment, because the equation is quadratic in x and there is no expression to write down; that block is a scalar Newton iteration at run time. The second gets vx from the once-differentiated constraint, which is linear in vx, so it is an assignment.

The last one is the interesting one: F and both accelerations, three unknowns determined together by three equations, one of them the twice-differentiated constraint. That is a linear solve inside every evaluation of the right-hand side, and no permutation removes it — it is a property of the model, which is what largest_block reports.

pendulum.index is 3, because the constraint had to be differentiated twice before the accelerations appeared at all. That is also why the model has two admissible state sets rather than one, which is the next notebook.

In [9]:

fig, ax = plt.subplots(figsize=(3.4, 3.4))
spy(
    ax,
    staircase(pendulum),
    "Pendulum — the 3×3 block, bottom right",
    ListedColormap([SURFACE, BLUE, ORANGE]),
)
plt.show()

fig, ax = plt.subplots(figsize=(3.4, 3.4)) spy( ax, staircase(pendulum), "Pendulum — the 3×3 block, bottom right", ListedColormap([SURFACE, BLUE, ORANGE]), ) plt.show()

Every example, on a dropdown¶

The same two pictures for anything in examples/models. DrainingTank and BouncingBall are tiny; RLCCircuit and MassSpringDamper show the alias story at hand-written scale; Pendulum is the one that does not sort.

In [10]:

LIBRARIES = [
    parse(path.read_text()) for path in sorted((root / "examples" / "library").glob("*.mo"))
]
MODELS = sorted(path.stem for path in (root / "examples" / "models").glob("*.mo"))


def show(model: str) -> None:
    """Raw incidence beside the causalized staircase, for one example model."""
    flat = flatten(
        [*LIBRARIES, parse((root / "examples" / "models" / f"{model}.mo").read_text())], model
    )
    system = causalize(flat)
    unknowns = [v.name for v in flat.unknowns]
    column = {name: i for i, name in enumerate(unknowns)}
    raw = np.zeros((len(flat.equations), len(unknowns)), dtype=int)
    for row, equation in enumerate(flat.equations):
        for name in {*names(equation.lhs), *names(equation.rhs)}:
            if name in column:
                raw[row, column[name]] = 1

    fig, (left, right) = plt.subplots(1, 2, figsize=(8.4, 4.2))
    spy(left, raw, f"{model}: {raw.shape[0]}×{raw.shape[1]} as flattened")
    spy(
        right,
        staircase(system),
        f"{len(system.blocks)} blocks, largest {system.largest_block}",
        ListedColormap([SURFACE, BLUE, ORANGE]),
    )
    plt.tight_layout()
    plt.show()

    print(f"index {system.index}   explicit {system.is_explicit}   states {system.states}")
    print(f"state sets {system.candidates}")


show("RLCCircuit")

LIBRARIES = [ parse(path.read_text()) for path in sorted((root / "examples" / "library").glob("*.mo")) ] MODELS = sorted(path.stem for path in (root / "examples" / "models").glob("*.mo")) def show(model: str) -> None: """Raw incidence beside the causalized staircase, for one example model.""" flat = flatten( [*LIBRARIES, parse((root / "examples" / "models" / f"{model}.mo").read_text())], model ) system = causalize(flat) unknowns = [v.name for v in flat.unknowns] column = {name: i for i, name in enumerate(unknowns)} raw = np.zeros((len(flat.equations), len(unknowns)), dtype=int) for row, equation in enumerate(flat.equations): for name in {*names(equation.lhs), *names(equation.rhs)}: if name in column: raw[row, column[name]] = 1 fig, (left, right) = plt.subplots(1, 2, figsize=(8.4, 4.2)) spy(left, raw, f"{model}: {raw.shape[0]}×{raw.shape[1]} as flattened") spy( right, staircase(system), f"{len(system.blocks)} blocks, largest {system.largest_block}", ListedColormap([SURFACE, BLUE, ORANGE]), ) plt.tight_layout() plt.show() print(f"index {system.index} explicit {system.is_explicit} states {system.states}") print(f"state sets {system.candidates}") show("RLCCircuit")

index 1   explicit True   states ('l.i', 'c.v')
state sets (('l.i', 'c.v'),)

In [11]:

import ipywidgets as widgets

_ = widgets.interact(
    show, model=widgets.Dropdown(options=MODELS, value="RLCCircuit", description="model")
)

import ipywidgets as widgets _ = widgets.interact( show, model=widgets.Dropdown(options=MODELS, value="RLCCircuit", description="model") )