From Diagram to Simulation: How Tools like Wolfram System Modeler Solve Physical Models with Modelica

Ever wondered what actually happens when you simulate a Modelica model? Let’s peek under the hood of tools like Wolfram System Modeler to see how your diagrams become dynamic simulations

What Really Happens When You Click “Simulate”?

Modelica-based modeling tools like Wolfram System Modeler make it easy to drag, drop, and simulate complex physical systems. But under the hood, these tools carry out a sophisticated process: translating high-level models into solvable systems of equations, transforming them symbolically, and solving them numerically.

In this article, I’ll walk you through each step—using a simple thermal system model—to reveal how tools like System Modeler convert your diagram into real-time simulation results. Whether you’re building energy systems, vehicles, or control models, understanding this internal workflow helps you build more efficient models, troubleshoot simulation issues, and optimize your modeling architecture.

Step 0: Define the Model

Consider a simple thermal circuit consisting of:

• A heat reservoir (a capacitor that stores thermal energy),
• A thermal resistor (limiting heat flow), and
• A fixed temperature source representing the ambient environment.

Here is the corresponding Modelica code:

model ThermalCircuit
  Modelica.Thermal.HeatTransfer.Sources.FixedTemperature ambient(T = 298.15);
  Modelica.Thermal.HeatTransfer.Components.ThermalResistor resistiveElement(R = 1);
  Modelica.Thermal.HeatTransfer.Components.HeatCapacitor heatReservoir(C = 10);
equation
  connect(ambient.port, resistiveElement.port_a);
  connect(resistiveElement.port_b, heatReservoir.port);
end ThermalCircuit;        

Once the model structure is defined, the solver needs to understand how components are linked. This happens via the connect() statements, which are usually acausal and require transformation.

Step 1: Expand connect() Statements

Modelica supports acausal modeling using connect() statements. These define how components are linked through standardized ports, but not how data flows explicitly. Behind the scenes, the tool performs connection expansion, applying:

• Potential variables equality at each node (e.g., T1 = T2)
• Flow variable conservation (e.g., Q_flow1 + Q_flow2 = 0)

For instance, the connection between the ambient and the resistor expands to :

resistiveElement . port_a . T = ambient . port . T;
ambient . port . Q_flow + resistiveElement . port_a . Q_flow = 0.0;        

Each thermal port contains:

• A potential variable (T, temperature),
• A flow variable (Q_flow, heat flow rate).

These expansions form the basis of the physical interactions.

Note: For causal connections (e.g., between input and output variables), the tool replaces connect() statements with direct assignments, such as input = output.

Step 2: Flatten the Model

Modelica models are hierarchical. Each component is itself defined as a model, potentially composed of submodels. Flattening means:

• Replacing component instances with their full definitions,
• Propagating parameter values and modifications,
• Resolving inheritance and redeclarations,
• Collecting all equations into a single equation system.

For example, the HeatCapacitor expands to:

  parameter SI.HeatCapacity C;
  SI.Temperature T(start = 293.15, displayUnit = "degC");
  SI.TemperatureSlope der_T(start = 0);
  Interfaces.HeatPort_a port;
equation
  T = port.T;
  der_T = der(T);
  C * der(T) = port.Q_flow;        

After flattening, the system is no longer object-oriented but instead resembles a large set of coupled equations—algebraic, differential, and discrete.

Here is the flattened version of the thermal model.

class ThermalCircuit

  parameter Real ambient.T(quantity = "ThermodynamicTemperature", unit = "K", displayUnit = "degC", min = 0.0, nominal = 300) = 298.15;
  ...
  ...
  Real heatReservoir.port.Q_flow(quantity = "Power", unit = "W");

equation

  heatReservoir.port.Q_flow + resistiveElement.port_b.Q_flow = 0.0;
  ambient.port.Q_flow + resistiveElement.port_a.Q_flow = 0.0;
  resistiveElement.port_a.T = ambient.port.T;
  heatReservoir.port.T = resistiveElement.port_b.T;
  ambient.port.T = ambient.T;
  resistiveElement.dT = resistiveElement.port_a.T - resistiveElement.port_b.T;
  resistiveElement.port_a.Q_flow = resistiveElement.Q_flow;
  resistiveElement.port_b.Q_flow = -resistiveElement.Q_flow;
  resistiveElement.dT = resistiveElement.Q_flow * resistiveElement.R;
  heatReservoir.T = heatReservoir.port.T;
  heatReservoir.der_T = der(heatReservoir.T);
  heatReservoir.C * der(heatReservoir.T) = heatReservoir.port.Q_flow;

end ThermalCircuit;        

Step 3: Sort and Transform Equations

The flattened model is now processed to:

• Identify initial conditions
• Separate differential vs. algebraic equations
• Handle events and when conditions

Since many systems yield Differential-Algebraic Equations (DAEs), tools use symbolic techniques like:

• Pantelides Algorithm or dummy derivative method for index reduction
• Block lower-triangular decomposition for equation sorting

The goal is to convert the system into a solvable form like:

dx/dt = f (x, y, t)
0 = g (x, y, t)        

Let’s look at an example workflow:

The following are the flattened equations generated after expanding component definitions and resolving connections.

heatReservoir.port.Q_flow[t]+resistiveElement.port_b.Q_flow[t]==0
ambient.port.Q_flow[t]+resistiveElement.port_a.Q_flow[t]==0
resistiveElement.port_a.T[t]==ambient.port.T[t]
heatReservoir.port.T[t]==resistiveElement.port_b.T[t]
ambient.port.T[t]==ambient.T
resistiveElement.dT[t]==resistiveElement.port_a.T[t]-resistiveElement.port_b.T[t]
resistiveElement.port_a.Q_flow[t]==resistiveElement.Q_flow[t]
resistiveElement.port_b.Q_flow[t]==-resistiveElement.Q_flow[t]
resistiveElement.dT[t]==resistiveElement.R resistiveElement.Q_flow[t]
heatReservoir.T[t]==heatReservoir.port.T[t]
heatReservoir.der_T[t]==der(heatReservoir.T)[t]
heatReservoir.C der(heatReservoir.T)[t]==heatReservoir.port.Q_flow[t]        

Remove assignment operations like a = b, and substitute parameters with their actual values:

ambient.port.T[t]==298.15
resistiveElement.dT[t]==ambient.port.T[t]-heatReservoir.T[t]
resistiveElement.dT[t]==resistiveElement.Q_flow[t]
heatReservoir.der_T[t]==der(heatReservoir.T)[t]
10 der(heatReservoir.T)[t]==resistiveElement.Q_flow[t]        

Bipartite Graph Construction and Equation Analysis

Once the symbolic system of equations has been simplified (e.g., by substituting parameter values and eliminating trivial assignments), the tool constructs a bipartite incidence graph to assist in further analysis and simulation preparation.

In this graph:

• One set of nodes represents equations.
• The other set represents variables.
• An edge connects an equation to a variable if the variable appears in that equation.

This representation allows Modelica tools to apply standard graph algorithms for tasks such as:

• Tearing, which helps break algebraic loops to make equation systems more solvable.
• Matching, which attempts to associate equations with variables in a way that supports an efficient solution process.
• Block decomposition, used to identify independent or weakly connected subsystems.

From Bipartite to Computational Graph

To sort equations into an executable order, the bipartite graph is transformed into a directed graph of dependencies, also known as an incidence graph or computational graph. In this version:

• Variables point to the equations that define them.
• Equations point to the variables they depend on.

Using Kahn’s algorithm or other topological sorting techniques, the tool identifies:

• Independent nodes, which can be solved first.
• The correct order of evaluation based on dependency chains.

Kahn’s Algorithm

• Find nodes (variables) with no dependencies (in-degree = 1).
• Remove these nodes and their outgoing edges.
• Repeat until all nodes are ordered.

However, this process can break down if the graph contains cycles, which indicate algebraic loops or simultaneous equations that must be solved together. In these cases, symbolic methods (e.g., tearing, index reduction, state selection) are applied to minimize loop size and isolate solvable blocks.

Let’s look at the cycle in the thermal model:

To break this cycle, the variable heatReservoir.T is selected as a state variable. This promotes its derivative, der(heatReservoir.T), to be handled by the integrator—removing the need to solve it algebraically within a feedback loop.

In addition to cycle breaking and state selection, Modelica tools also handle events such as zero-crossings and discrete changes, enabling the simulation of hybrid systems that combine continuous and discrete dynamics.

Once symbolic preprocessing is complete, the numerical solver (e.g., DASSL or CVODES) performs time integration, computing the system’s behavior across the simulation interval.

Summary

When you simulate a Modelica model, the tool performs a series of symbolic transformations to prepare it for numerical solving.

• First, it expands connect() statements into physical equations based on potential and flow variables.
• Next, it flattens the hierarchical model into a unified system of equations—resolving inheritance, parameters, and component definitions.
• These equations are then sorted, simplified, and transformed using techniques like index reduction (e.g., Pantelides algorithm) and block decomposition. A dependency graph is constructed and topologically sorted to determine a solvable computation sequence.
• Finally, numerical solvers like DASSL or CVODE perform time integration, simulating the system’s behavior dynamically.

Each of these steps—connect, flatten, sort, simulate—is essential for turning your visual model into accurate, time-domain simulation results.