• Due Date: Friday 06 November 2026, before 23:59.
• Team Size: 2 (pair design/programming).
  Note that as of the 2017-2018 Academic Year, each International student should team up with a "local" student (i.e., whose Bachelor degree was obtained at the University of Antwerp).
• Submission Information:
  • Only one member of each team submits a full solution.
  • This must be a compressed archive (ZIP, RAR, TAR.GZ...) that includes your report and all models, images, code and other sources that you have used to create your solution.
  • The report must be in PDF and must be accompanied by all images.
  • Images that are too large to render on one page of the PDF must still be rendered on the PDF (even though unclear), and the image label should mention the appropriate file included in the submission archive.
  • Make sure to mention the names and student IDs of both team members.
  • The other team member must submit a single HTML file containing the coordinates of both team members using this template.
• Submission Medium: BlackBoard Ultra. Beware that BlackBoard's clock may differ slightly from yours. If BlackBoard is not reachable due to an (unpredicted) maintenance, you submit your solution via e-mail to the TA. Make sure all group members are in CC!
• Contact / TA: Rakshit Mittal.

In this assignment, you will learn to create and use Causal-Block Diagrams (CBD) in Python. You will learn:

• How to simulate CBD models (orchestration, dependency graphs, topological sorting, and algebraic loops) using Python.
• About the numerical analysis, accuracy and stability of simulations.
• How to represent CBDs in a human-readable textual notation.
• About the optimization of simulator performance.

To get a better understanding about how tools like OpenModelica and MatLab Simulink work behind the scenes.

You are discouraged but allowed to make use of any GenAI tool in solving the assignment or writing the report (for example, to correct grammatical mistakes) as long as you adhere to the UA Guidelines for students on responsible use of GenAI.

You are required to cite any use of GenAI and describe what portions of the assignment you used it for.
You should also mention all other sources, including collaborations.

Note that a significant part of demonstrating that the learning goals have been achieved, includes being able to

1. explain the relevant concepts in the assignment,
2. explain the design choices in your implementation, and
3. critically discuss your solution.

This will be evaluated in a short (~15 minute) oral evaluation, which will be conducted in the week following the deadline of every assignment.

DrawioConvert is a tool that allows the creation, visualization and code generation of coupled CBDs. It uses DrawIO as a front-end. A small tutorial on how to use the tool is provided.

You should first install DrawIO from this link

Essentially, using this tool, you can define your CBDs in a graphical environment (much like the graphical editor of OMEdit), and then use the DrawioConvert Python code (provided in template.zip) to generate Python code that instantiates a CBD model based on the CBD framework.

This assignment consists of four main parts:

1. Implementing a CBD simulator in Python (35%)
2. Circle test and numerical accuracy (25%)
3. Automatic generation of LaTeX equations (20%)
4. Optimization of the CBD simulator and integration of FMUs (4 × 5% = 20%)

The operational semantics of a Causal Block Diagram is defined by a simulator. The simulator turns the diagram into an order of computations, then evaluates that order at every time step. Roughly, the steps are:

1. Flatten the model. CBDs can have hierarchy but only the atomic blocks actually compute a value.
2. Build a dependency graph for the current iteration. A block depends on the blocks that produce its inputs.
3. Sort the graph into strongly connected components.
   • A component of one block is normally an ordinary computation.
   • A larger component is an algebraic loop.
   • A single block that depends on itself is also an algebraic loop.
4. Compute or solve each component.
   • An ordinary block is evaluated directly.
   • A linear algebraic loop is solved as a system of equations.
   • For this assignment: A non-linear loop stops the simulation and raises a clear error.
   • A co-simulation FMU is evaluated by calling its "doStep" method.
5. Store the signals, and copy them onto the output ports of every hierarchical CBD.
6. Advance time by one step and repeat.

You will implement this simulator in Python. Use this template to get started. It contains DrawioConvert and the skeleton of the CBD framework. The CBD framework is split into four files:

• model.py defines the data structures.
• simulator.py implements the dependency graph builder, the scheduler and the simulation loop.
• linear_solver.py implements the linear algebraic-loop solver.
• library.py defines the basic blocks.

Task 1: Flatten the CBD model

A CBD model is a tree. Its leaves are atomic blocks and only they compute values. Its intermediate nodes are also instances of CBD. A CBD groups CBDs/blocks and has ports. The simulator schedules the leaves, so the first thing it needs is the list of leaves.

Complete Simulator._collect_atomic_blocks in simulator.py.

Hint: see CBD.getBlocks.
Hint: A child can itself be a CBD.

Task 2: Build the dependency graph

Before the simulator can decide what to compute first, it has to know where each input gets its value.
Look at the data structures in model.py:

• Ports:
  • Every block owns one Port per input and per output.
  • A port keeps a history of Signal(time, value) entries.
  • Port.set appends Signal objects to that port only.
• Connections:
  • A Connection joins a source port to a target port.
  • An input port has at most one incoming connection, returned by getIncoming.
  • An output port may feed many targets, returned by getOutgoing.
• Ports on CBD objects (hierarchy):
  • A CBD has ports too.
  • A CBD port computes nothing. It is only a junction along the wire.

So the value an input reads in a simulation step actually sits in the history of an output of some atomic block, possibly several CBD boundaries away.
Complete three methods in model.py that find it:

• Port.getSourcePort walks back from a port to the port that produces its value.
  Hint: Think about which kind of port ends the walk.
  Hint: What should happen when the walk reaches a port it has already visited?
• BaseBlock.getInputSignal returns the signal an input sees at an iteration. Blocks call it from compute.
• BaseBlock.getDependencies returns the atomic blocks a block reads from at an iteration. It must never return a CBD.

Then, in simulator.py:

• Complete Simulator._build_dependency_graph(curIteration).
  Hint: DependencyGraph, which is already given, defines the data structure for the dependency graph.
• The dependency graph is passed to the method, strongly_connected_components (already given), which implement's Tarjan's algorithm to group the blocks into strongly connected components.
• The groups are passed to the method, topological_sort (already given), which orders the groups so that every group comes after the groups it reads from.

Task 3: Solve linear algebraic loops

Before we move onto finishing the rest of the simulator, we need to ensure that we can deal with the components (from the previous task) that are linear algebraic loops.

No block in an algebraic loop can be computed first, because each block needs another block's output from the same iteration. When every block in the loop is linear, the loop is a system of linear equations, and you can solve the whole loop at once.

To understand how we will solve a linear algebraic loop, let's take the following CBD:

Figure 1: Exemplar CBD with a linear algebraic loop representing $x = 2 -(0.5 * x)$
The solution is $p = 2/3$, $n = -2/3$ and $a = 4/3$, so $x = 4/3$ at every iteration.

The loop has three unknowns: the product output $p$, the negator output $n$ and the adder output $a$.
The two constants are outside the loop, so they are already known (from sorting the components, they are computed before this algebraic loop).
Write one equation per block, then move all the unknowns (the variables) to the left:

p = 0.5 * a     becomes    -p        + 0.5 a  =  0
n = -p          becomes    -p  -  n           =  0
a = 2 + n       becomes        +  n  -     a  = -2
				

With the unknowns ordered $(p, n, a)$, that is $A x = b$ with

A = [[-1,  0,  0.5],        b = [ 0,
     [-1, -1,  0  ],              0,
     [ 0,  1, -1  ]]             -2]
				

Each row in the matrix $A$ is one block, and each column is one unknown.
A and b are used by numpy.linalg.solve to solve the linear system of equations.

Complete the three methods in linear_solver.py:

• is_linear decides if an algebraic loop is linear. You need to determine whether every block of the component implements a linear equation. Decide which block types qualify.
  Hint: Adders, negators and delays are linear.
  Hint: A product is linear only when fewer than two of its inputs are unknowns of this same loop.
  Hint: An inverter, a generic block, a comparison, or a product of two unknowns is not linear.
• build_system returns (A, b) for a loop. Row i is the equation of component[i].
• solve_linear_loop solves the system with numpy.linalg.solve and appends each result to the output of its block.
  Hint: numpy raises a LinAlgError when the system has no solution or has infinitely many solutions.

Task 4: Implement the simulation step

Finally, we can piece together all the components of our CBD simulator in simulator.py.

The main Simulator.run is already given.
It calls Simulator._step(curIteration, time) once per iteration and advances time.

Complete Simulator._step so that it executes the following instructions in each iteration:

• Compute the dependency graph for the current iteration.
• Group the blocks into strongly connected components.
• Order the groups so that every group comes after the groups it reads from.
• Walk the ordered groups:
  • An ordinary block computes itself with compute(curIteration, time, delta_t).
  • A group of more than one block, or a single block that depends on itself, is an algebraic loop.
  • If the loop is linear, solve it with solve_linear_loop.
    Hint: You implemented the is_linear method in Task 3.
  • Otherwise, raise a ValueError that names the blocks involved with the message that a non-linear algebraic loop was detected.
• Update the hierarchical CBD ports so that each CBD output carries the same signal as the atomic output behind it.
  After walking the groups, every atomic output holds its signals. However, the output port of a hierarchical CBD is only a junction, so its history is still empty, and model.getSignalHistory() on a nested model returns nothing.
  Complete Simulator._update_cbd_ports so that each hierarchical CBD output carries the same signal as the atomic output behind it.
• Advance time by one step.

Task 5: Complete the basic block library

You almost have a full-fledged CBD simulator! There is however, one last thing left to do (you might have already noticed it during Task 3).
You need to define the library of basic blocks that you will reuse to build your own CBDs!

Complete the definition of the basic blocks in library.py:

• Each class already declares the ports of the block, as they need to be consistent with the DrawioConvert library.
• The ConstantBlock and AdderBlock are already completed as examples.
• The methods that raise NotImplementedError are for you to complete.

Leave the IntegratorBlock and DerivatorBlock for Part 2 of the assignment.

Task 6: Test and report

Create a CBD for the following three (sets of) equations and test your simulator with them.
Simulate the models from $t = 0$ to $t = 20$, with $\Delta t = 0.1$.
For each, show the model (DrawIO diagram AND Python code generated by DrawioConvert), and plot the value of $y_i$ against $t_i$.

1. $\displaystyle d_i = \begin{cases} 17 & \text{if } i = 0 \\ \sin(t_{i-1}) & \text{if } i \ge 1 \end{cases} $
   
   $ y_i = \dfrac{5}{3} + d_i - t_i\,\Delta t $


2. $ d_i = \begin{cases} 7 & \text{if } i = 0 \\ -d_{i-1} & \text{if } i \ge 1 \end{cases} $
   
   $ y_i = - ( d_i \cdot ( 11 + t_i + y_i ) ) $


3. $\displaystyle y = \frac{-1}{2\cdot(5 + y)}$

In the report, describe your solution for Tasks 1-5. Explain, in particular (for example, by showing your code):

• how the dependency graph is built,
• how hierarchy is handled,
• how a linear loop becomes a matrix,
• how the blocks are defined in the library,
• and how the simulator works.

In the lectures, you have seen the undamped harmonic oscillator is $\ddot{x} = -x$.
With $x(0) = 0$ and $\dot{x}(0) = 1$, the solution is $x = \sin t$ and $v = \dot{x} = \cos t$.
The point $(x, v)$ stays on the unit circle, which is why this is also called the circle test.

A numerical method drifts off that circle, and the drift is a measure of its accuracy and stability.

Task 1: Integrator and derivator blocks

Complete IntegratorBlock and DerivatorBlock in library.py.
First implement them graphically with DrawIO, and then convert them with DrawioConvert to get the Python code. Copy the code to the respective methods in library.py.
Hint: Both are hierarchical CBDs built from the blocks in Part 1.
Hint: Each has inputs IN1 ,DT, and IC, and output OUT1.

• Integrator: $y_0 = IC$, and $y_i = y_{i-1} + \Delta t \cdot u_{i-1}$.
• Derivator: $y_0 = IC$, and $y_i = (u_i - u_{i-1}) / \Delta t$.

The integrator implements the first-order forward Euler method.
The derivator implements the first-order backward Finite difference method.

Report the DrawIO diagrams and the Python code generated by DrawioConvert for both blocks.

Task 2: Model the oscillator with DrawIO

Build two models of the oscillator in DrawIO, convert them with DrawioConvert, and simulate them with your simulator:

• Integral form. $x = \int -v\,dt$ and $v = 1 + \int x\,dt$.
  Use two IntegratorBlocks, a NegatorBlock and two constants.
  
  
• Derivative form. $x = d(-v)/dt$ and $v = dx/dt$.
  Use two DerivatorBlocks, a NegatorBlock and two constants.

Report the DrawIO diagrams and the Python code generated by DrawioConvert for both models.

Are there any algebraic loops in either model? Describe them in your report.

Task 3: Simulate with different time steps

Simulate both models from $t = 0$ to $t = 20$, with $\Delta t = 0.1$, $\Delta t = 0.01$ and $\Delta t = 0.001$. For each model, give six plots:

• $x(t)$ for all three step-sizes;
• the phase plot, $v$ against $x$, for all three steps.

Task 4: Numerical accuracy

Compute the analytical $x(t)$ and $v(t)$ in Python, with math or numpy, at the same sample times as the CBD trace. For each model and each $\Delta t$, report:

• the mean-squared error for the $x$ signal by comparing the CBD trace to the analytical $x(t)$;
• the mean-squared error for the $v$ signal by comparing the CBD trace to the analytical $v(t)$;

Give these numbers in one table, then answer:

• Which form (if any) spirals outward, and which spirals inward? Why?
• How does the error change when $\Delta t$ is divided by 10? What order of accuracy does that suggest?
• Which form would you trust for a long simulation, and why?

Hint: You have to report 12 mean-squared errors in total.

Mathematical models are often documented using equations. For CBD models, we can automatically generate LaTeX code that represents the model's denotational semantics - the mathematical equations describing what each signal computes.

This is useful for:

• Documentation and presentations
• Mathematical analysis and verification
• Understanding the model at a glance
• Academic papers and reports

You will create a Python module that takes a CBD model and generates LaTeX code representing its equations in three different forms:

1. CT-CBD (Continuous Time): The original continuous-time model's denotational semantics
2. DT-CBD (Discrete Time): The inlined/discretized model as a set of equations (order doesn't matter)
3. Sorted DT-CBD: The computation order with dependencies resolved (order matters, except for algebraic loops)

Task 1: Generate LaTeX for CT-CBD Model

For the original continuous-time CBD model (before discretization), generate the denotational semantics showing the continuous behavior.

For example, an integrator in continuous time is represented as:

    \[
      \frac{dx}{dt} = f(t, x)
    \]

This represents the mathematical intent of the model before discretization.

Task 2: Generate LaTeX for DT-CBD Model (Inlined/Discretized)

After discretization (e.g., using Backwards or Forwards Euler), generate the discrete-time equations. These should be formatted as a set of equations using \left\{, indicating that the order doesn't matter semantically.

Format as a set of equations:

    \[
      \left\{
        \begin{array}{rcl}
         v(i) & = & w(i) - 2\\
         w(i) & = & w(i-1) + v(i-1) + 3\\
        \end{array}
      \right.
    \]

Which renders as:

The \left\{ notation indicates these are the denotational semantics of the DT-CBD - what the equations mean, regardless of computation order.

Task 3: Generate LaTeX for Sorted DT-CBD Model

After applying topological sort from Part 1, generate the equations in computation order. Use \[ (without the set braces) to show this is a sequential computation where order matters.

    \[
      \begin{array}{rcl}
       w(i) & = & w(i-1) + v(i-1) + 3\\
       v(i) & = & w(i) - 2\\
      \end{array}
    \]

Important: If there are algebraic loops (strongly connected components), they must still be rendered as a set with \left\{ since they need to be solved simultaneously:

    \[
      \begin{array}{rcl}
       a(i) & = & c(i) + 1\\
       \left\{
         \begin{array}{rcl}
          b(i) & = & c(i) \cdot 2\\
          c(i) & = & b(i) + d(i)\\
         \end{array}
       \right. \\
       d(i) & = & a(i) - 3\\
      \end{array}
    \]

Hints:

• Use block.getPath() to get full hierarchical names (e.g., "controller.pid.sum")
• Sanitize names for LaTeX (replace underscores with escaped underscores: \_)
• Use the topological sort from Part 1 to order equations for the sorted DT-CBD
• For DelayBlocks, you want to generate two equations or use piecewise notation
• LaTeX math operators:
  • Addition: +
  • Multiplication: \cdot or \times
  • Division: \frac{a}{b} or a / b
  • Functions: \sin, \cos, \exp, \log
  • Subscripts: x_{i-1}
  • Derivatives: \frac{dx}{dt}

Report your LaTeX generator implementation

For each of the following models, show all three LaTeX representations (CT-CBD, DT-CBD, and Sorted DT-CBD) with compiled PDF screenshots:

1. All three models from Part 1 - Task 6
2. Both forms of the harmonic oscillator from Part 2 of the assignment.

For the first three tasks, .
You should copy simulator.py into optimized_simulator.py and modify it for Tasks 1 - 3.
You do not need to create a new simulator class for each optimization.
At the end of each task, show the effect of the optimization on some concrete model.
Along with other proof of optimization (such as explanation of your code, terminal output, etc.), you must also have some timing analysis.
Hint: Time your simulation atleast 100 times and take the average, for a more robust analysis.

Task 1: Constant folding (5%)

The output of some blocks never changes during a run.
Think of ConstantBlock, DeltaTBlock, etc.
Moreover, by transitivity, if a block depends only on such blocks with constant outputs, the block's output is a constant too.

We can "fold" these blocks by recording their value at the beginning of the simulation and then replaying the value at every subsequent timestep.
We save repeated calls to compute each such block in every timestep.
How would you deal with DelayBlocks in your CBD if you want to implement "constant folding"?

Modify the Simulator class in optimized_simulator.py to implement constant folding.
Hint: To find all constant blocks, traverse the dependency graph of iteration 1, not iteration 0 (to deal with DelayBlocks).
Hint: Record the signal value at time 0, and then replay the value at every subsequent timestep for all the constant blocks.

A similar optimization can also be made for the linear algebraic loop solver.
Remember that build_system in linear_solver.py returns (A, b).
The entries of A are the coefficients of the variables in the equations. They are either 1, -1, or some multiplicative input for the unknown.
When every one of those inputs comes from a block you folded, A is the same at every steady-state iteration.

Build and cache the matrix A using the dependency graph of iteration 1. At later iterations, rebuild only b.

Task 2: Cached steady-state schedules (5%)

The schedule of computation of the components does not change between iterations.
Hence, the dependency graph generation, component search and topological sort do not have to run at every iteration.
Why? Because other than DelayBlocks, the dependency graph does not change after iteration 1. Connections between blocks are static, they do not change during the simulation.

Modify the Simulator class in optimized_simulator.py to implement this optimization.

Task 3: Scope-Based Computation & Multi-Rate Scheduling (5%)

Not all signals need to be computed at every iteration.
Memory blocks (DelayBlock, IntegratorBlock, etc.) must be updated every Integration Interval (II).
However, high-level output signals may only need updating every Communication Interval (CI).

• Integration Interval: Time step for numerical integration ($\Delta t$)
• Communication Interval: Time step for output/tracing ($CI = k \times II$)
• "Must-Compute" Part: Blocks feeding into memory blocks (must run every II)
• Scope Block: Special block marking signals to trace (run every CI)

Implement multi-rate scheduling.

1. The ScopeBlock is already provided in library.py. Use it to mark signals for tracing.
2. Identify "must-compute" blocks: find all memory blocks and traverse backwards to find their dependencies
3. Identify "scope-required" blocks: find all blocks feeding into ScopeBlocks
4. Remove unused blocks: any block not in must-compute AND not in scope-required
5. Make two functions step_ii() and step_ci() to handle the two schedules respectively:
   • step_ii(): Must-compute schedule: runs every iteration
   • step_ci(): Full schedule (including scope-required blocks): runs every K iterations

Task 4: Modelica plant FMU co-simulated with a CBD PID controller (5%)

From your Modelica assignment solution, export a Functional Mockup Unit (FMU) of the quarter-car plant without the PID controller or the setpoint input.
Hint: Before exporting your Modelica model as an FMU, ensure that it has one RealInput and one RealOutput.

Figure 2: Co-simulation model. The plant FMU is exported from Modelica, while the PID controller is built from CBD library blocks.

Create a CBD model (directly in Python, as there is no graphical counterpart for the FMU block) implementing the same control loop as in the Modelica assignment (shown above).

Complete the FMUBlock class in library.py to load the FMU.

A minimal Python interface to an FMU is provided in fmi2.py. Study the contents of this module, and remember what you have learnt in the tutorial about co-simulation with FMI.
complete the FMUBlock class in library.py to load the FMU successfully in your simulator.

At iteration 0, the FMU should be instantiated, and its output value should be obtained and attached to the block's output signal.
Hint: To do this use the _open method to instantiate the FMU, and the get_real method to get the output value.
Hint: You can get the name of the FMU's output signal from the _output_reference attribute of the class.

At other iterations, you need to find the input value for the FMU emitted from the PID blocks. This value is the force command applied to the plant, and it is computed on the previous iteration.
Hint: Use the getInputSignal method to get the value of the force command from the previous iteration.
Then you need to do set_real on the FMU to set the input value, do_step to step the FMU by $\Delta t$ and then finally get_real to get the output value and attach it to the block's output signal.

Simulate your CBD model for the same time duration and time step as in the Modelica assignment, with the best PID parameters you obtained in the Modelica assignment, and the same setpoint $0.5\,\mathrm{m}$.
Report the graph and compare the output with the Modelica simulation you performed in the previous assignment.

• All parts of this assignment use Python 3.6+. Do not use features from Python 2.7 (discontinued as of January 1, 2020).
• The template for this assignment is available here.
• Install NumPy (pip install numpy). Part 1 uses numpy.linalg.solve for linear algebraic loops.
• A Small Tutorial on DrawioConvert is provided.
  Note that DrawioConvert requires Jinja2
• Useful references:
  • CBDs: http://msdl.uantwerpen.be/people/claudio/pub/Gomes2016a.pdf
  • NumPy: numpy.linalg.solve
  • Dependency graphs and topological sorting: Lecture slides