• Due Date: 19 November 2022, 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 "local" (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 crate your solution. Do not submit any executables. The report may be either HTML or PDF and must be accompagnied by all images. When an image is unreadable in your report and missing from your submission archive, you will not receive any points for that task. Make sure to mention the names and student IDs of both team members. The other team member must submit a single HTML file containing only the coordinates of both team members. You may use this template. This will allow us to put in grades for both team members in BlackBoard.
• 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: Randy Paredis.

In this assignment, you will learn how to use Causal-Block Diagrams (CBD) for modeling and simulation. You will take a closer look at specific simulation techniques (mainly numerical analysis, accuracy and stability) in order to get a better understanding about how tools like OpenModelica and MatLab Simulink work.

This assignment is to be made using this Python-based CBD framework. The ZIP-file contains a full framework to allow CBD modelling in Python. To build the documentation, use the docs.bat or docs.sh scripts. All additional information about this framework can be found in this documentation. Additionally, the documentation provides a set of examples which may help you in completing this assignment. Especially the "Sine Generator" example. Note that the internal implementation will be much more complicated copared to the algorithms seen in class. If something is not clear, or if you have found an issue/bug, or if you have a feature request, please contact the TA.

Once the simulator is installed, it works with any IDE, including Jupyter Notebooks. For creating plots, Bokeh, MatPlotLib and SeaBorn are supported.

DrawioConvert is an additional tool that allows the creation, visualization and code generation of coupled CBDs. The tool may still have small bugs, and its use is optional, however very useful. It uses DrawIO as a front-end. A small tutorial on how to use the tool is provided. Please contact the TA if you experience any problems.

• DrawIO Library
• Download Python Code
• Requires: Jinja2

A tutorial on how to use this framework and DrawioConvert is presented in this video.

1. Harmonic Oscillator
2. Integration Methods
3. Stiff Systems of Equations

A simple harmonic oscillator with no friction exhibits behaviour that can be modelled by the following second-order ODE: $$\dfrac{d^2(x)}{dt^2} = -x$$ where $x(0) = 0$ and $\dfrac{d(x)}{dt}(0) = 1$. This system can be modeled in the CT-CBD formalism using IntegratorBlocks and/or DerivatorBlocks. The exact (analytical) solution of this equation is $$x(t) = \sin(t)$$

One approach to simulate this system is to make use a small, fixed stepsize. The smaller, the more accurate the solution will be. Another is to make use of adaptive step-size (also known as variable step-size or step-size control). Instead of computing your function using a fixed $\Delta t$, the value of $\Delta t$ changes throughout its execution, where the delta between iterations $h_i$ (with $i$ the iteration number) is defined by the rate of change of the function to compute. A common family of adaptive step-size algorithms are known as the Runge-Kutta (RK) methods. The most famous algorithm is known as the Runge-Kutta-Fehlberg of the 4th and 5th order (RKF45). Luckily for you, you do not have to concern yourself with the math behind it. The CBD simulator that is provided allows a model transformation which transforms your CBD model into an equivalent, RK-based CBD. Take a look at "Continuous Time Simulation" in the documentation to read more about this.

Tasks

Your job is to find the best approach to model the above harmonic oscillator. This is done as described below. Make sure to include all your source code and everything you believe is required for your solution.
Hint: it may be possible to do multiple of these tasks in a single for-loop.

1. Build a continuous-time, coupled CBD that makes use of IntegratorBlocks. Let's call this block CBDA and its output $x_A$. Include an image of your CBD.
2. Simulate CBDA for $10$ seconds and include a plot of $x_A$ over the time.
3. Build a continuous-time, coupled CBD that makes use of DerivatorBlocks. Let's call this block CBDB and its output $x_B$ Include an image of your CBD.
4. Simulate CBDB for $10$ seconds and include a plot of $x_B$ over the time.
5. Create a coupled CBD to output the analytical solution $\sin(t)$. Include an image of your CBD.
6. Create a coupled CBD called ErrorA, which computes the error between $x_A$ and $\sin(t)$. The error metric is given as $e_A(t) = \int\vert\sin(t) - x_A(t)\vert\ dt$. Include an image of your CBD.
7. Simulate ErrorA for $50$ seconds with $\Delta t = 0.1$, $\Delta t = 0.01$ and $\Delta t = 0.001$. Include a plot of $e_A(t)$.
8. Create a coupled CBD called ErrorB, which computes the error between $x_B$ and $\sin(t)$. The error metric is given as $e_B(t) = \int\vert\sin(t) - x_B(t)\vert\ dt$. Include an image of your CBD.
9. Simulate ErrorB for $50$ seconds with $\Delta t = 0.1$, $\Delta t = 0.01$ and $\Delta t = 0.001$. Include a plot of $e_B(t)$.
10. Discuss the obtained results. Which deltas are better? What results in a more accurate system; an integrator or a derivator? Why?
11. Prove the validity of all the given formulas using the $\LaTeX$ generator. Take a look at "Generate LaTeX from CBD Models" (in CBD.converters) in the documentation for more details.
12. Transform CBDA to a new CBD using RKF45. Notice that this has fully transformed your CBD model, presumably to a point of no recognition. Include an image of the flattened version of this model. Highlight the parts that you presume to be "remnants" of your original model.
    • Note that you must use the addFixedRateClock() method to validly do a conversion. Use $\Delta t = 0.1$ (this is the initial step-size).
    • Make sure to use the CBD with the builtin IntegratorBlock for this task.
    • The flattened model may be massive and therefore unreadable in your report. Do not forget to add the figure (and preferrably the source files (e.g., *.gv, *.dot, *.puml...)) to your report.
    • Use $atol = 2\cdot 10^{-5}$, $h_{min} = 0.01$ and $safety = 0.84$.
13. Simulate the RKF45 CBD for $10$ seconds; do not set a simulator delta. Plot the new $x$ and compare this result to the other CBDs, as well as the analytical solution.
14. Add an indication of when the simulator computed the values to this plot, similar as to the examples given in the documentation. Describe what you see. How many computations have happened as compared to the other approaches?

Computing the integral of a function is at the heart of solving continuous-time models. Their computation on digital computers always requires a numerical approximation. This task will take a look at multiple integration methods and their accuracy. Table 1 shows the four integration methods we will concern ourselves with for this assignment (note: there exist many, many more).

Table 1: Four numerical integration schemas.

Backward Euler Method	Forward Euler Method	Trapezoid Rule	Simpson's 1/3 Rule
$$I_0 = f(t_0)\cdot \Delta t$$$$I_1 = f(t_1)\cdot \Delta t$$$$I = I_0 + I_1$$	$$I_0 = f(t_1)\cdot \Delta t$$$$I_1 = f(t_2)\cdot \Delta t$$$$I = I_0 + I_1$$	$$I_0 = \dfrac{f(t_0) + f(t_1)}{2}\cdot \Delta t$$ $$I_1 = \dfrac{f(t_1) + f(t_2)}{2}\cdot \Delta t$$ $$I = I_0 + I_1$$	$$I = \dfrac{\Delta t}{3}\left(f(t_0) + 4\cdot f(t_1) + f(t_2)\right)$$
computed over 1 iteration	computed over 1 iteration	computed over 1 iteration	computed over 2 iterations

The RK method that was mentioned earlier do not only allow adaptive step-size, but also provide a way to more accurately solve integrals. Notice how the forward Euler method, the trapezoid rule and Simpson's 1/3 rule are special cases of RK.

Tasks

Your job is to compare the given integration methods.

1. For each of the given integrators, implement a coupled CBD.
   • The backward Euler method is already implemented in the CBD framework in CBD.lib.std.IntegratorBlock.
   • Each integrator should have an initial condition input (which is called IC), that outputs this value at the first iteration.
   • Include images that show the created CBDs visually.
   • Use the $\LaTeX$ generator to prove the validity of the used formulas.
   • For Simpson's 1/3 rule, the trapezoid rule should be used if you have but two data points.
   • You are allowed to model the composite Simpson's rule. This is a more complicated formula, which will be reflected in your grade (given a valid implementation).
2. For the function $g(t) = \dfrac{t}{t^2 + 1}$, it is known that the integral over the range $[0, 100]$ is: $$\int_0^{100} \dfrac{t}{t^2 + 1} dt = \dfrac{1}{2}\left[\ln(t^2+1)\right]_0^{100} = 4.60522018$$ Implement $g(t)$ as a coupled CBD block.
3. Compute and compare the integral of $g(x)$ for $\Delta t = 0.1$, $\Delta t = 0.01$ and $\Delta t = 0.001$, using the four different integrators. Make sure to also compare against the analytical solution! Which one(s) is/are the most accurate?
   Hint: this will result in 13 different values (= 3 different $\Delta t$, using 4 integrators + 1 analytical solution) to compare.

A set of equations is called stiff if it is numerically unstable to solve stepwise, unless the $\Delta t$ is extremely small. There is no precise definition of stiffness, but a characteristic is that very small stepsizes are required even though the solution curve is smooth. We will consider a simple "stiff" Initial Value Problem (IVP) in this task: $$\dfrac{d(y)}{dt} = -15y;\ y(0) = 1$$ The exact solution of this equation is $y(t) = e^{-15t}$.

Tasks

Analyze the stiffness of the above-given IVP.

1. Create a coupled CBD for the IVP for each of the four integrators.
2. Simpluate each CBD for $1$ second, using $\Delta t = 0.25$, $\Delta t = 0.125$, $\Delta t = 0.1$ and $\Delta t = 0.05$.
3. Plot the output $y$ for each of these CBDs. Compare this result against the real (analytical) solution.
   Hint: You can use numpy/scipy to create a plot for the real solution.
4. Describe which CBDs produced the best result.
   Hint: you can define an error function to mathematically find the "best" approximation.
5. Analyze the usage of the forwards Euler integrator against the backwards Euler. What is the difference between the two (in terms of this problem description)?

• All parts of this assignment use Python 3.6+. Do not use features from Python 2.7 (discontinued as of January 1, 2020).
• The CBD simulator can be found here.
• The files for DrawioConvert can be found here:
  • Small Tutorial
  • DrawIO Library
  • Download Python Code
  • Requires: Jinja2
• Useful documentation:
  • The CBDsim documentation
  • A small video tutorial on the use of the CBDsim framework and drawioconvert
  • The CBD lecture slides
  • Background material: an elaborated version of the CBD story. Claudio Gomes, Joachim Denil, and Hans Vangheluwe. Causal-Block Diagrams: A family of languages for causal modelling of cyber-physical systems. In Paulo Carreira, Vasco Amaral, and Hans Vangheluwe, editors, Foundations of Multi-Paradigm Modelling for Cyber-Physical Systems, chapter 4, pages 97-125. Springer International Publishing, 2020. [pdf]