• Due Date: Sunday 26 November 2023, before 23:59.
• Team Size: 2 (pair design/programming)! Please contact the TA ASAP if you encounter any team issues.
  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 accompanied by all images. When an image is unreadable in your report and missing/unreadable 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. 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 behind the scenes.

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 incredibly 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

This assignment consists of three main parts:

1. Harmonic Oscillator
2. Integration Methods
3. Inline Integration (and Cosimulation)

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 CBDs using IntegratorBlocks and/or DerivatorBlocks. The real (analytical) solution of this system equals $\sin(t)$.

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 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?

Computing the integral of a function is a common practice in simulation systems. It can be used for specifying time-dependent simulations, constructing (PID) controllers and even identifying accumulative errors. Any complex Cyber-Physical process that happens over a certain period of time makes use of integrators, however their computation is always 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: Comparison of multiple integration formulas.

Backwards Euler Method	Forwards Euler Method	Trapezoid 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$$
computed over 1 iteration	computed over 1 iteration	computed over 1 iteration

Tasks

Your job is to compare the given integration methods.

1. For each of the given integrators, implement a coupled CBD.
   • The backwards 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. Make sure they are clearly readable.
2. For the function $g(t) = \dfrac{t - 1}{(t + 1)^2}$, it is known that the integral over the range $[0, 100]$ is: $$\int_0^{100} \dfrac{t - 1}{(t + 1)^2} dt = \left[\dfrac{2}{t+1}+\ln(t+1)\right]_0^{100} = 2.634922497$$ 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 10 different values (= 3 different $\Delta t$, using 3 integrators + 1 analytical solution) to compare.

Note: This task assumes you know the C programming language. A simple tutorial can be found on the W3 Schools Website.
Additionally, it is assumed you have fully completed the previous (Modelica) assignment.

Even though you have a CBD simulator, its execution in Python is slow and does a lot of things behind the scenes (feel free to take a look at the source code). Sometimes, models need to be executed on hardware that may not have the computational power (and storage requirements) to run the model you have built. Remember how Modelica solved this problem: it first flattens all models, next it converts it into executable C-code and finally, it compiles and executes the model.

The executable C-code is created by doing a topological sort for all blocks and iteratively turning each block into their corresponding equation(s). This also takes the internal connections into account. This method is called "Inline Integration" This paper (see especially table 2) somewhat explains the process for inline integration of CBDs (mainly focusing on LaTeX generation). Notice that the dependency graph at iteration 0 will differ from the graph at all other times, so the inline integration needs to happen twice. For the purposes of this assignment, we will ignore how algebraic loops can/should be solved.

Another reason for doing this is that we might want to combine multiple formalisms together. For instance, when you have a causal controller in CBDs and an acausal plant in Modelica. Combining these purely on a code level is an easy way to allow multi-formalisms. However, sometimes there may be algebraic loops between multiple formalisms, causing complicated code and a lot of headaches. Luckily, there is another way of doing this: cosimulation, in which an orchestrator manages interaction between multiple formalisms. The industrial standard for cosimulation is the Functional Mock-up Interface (FMI). A single instance is called a Functional Mock-up Unit (FMU). For this assignment, we will use FMI standard 2.0.

For this task, you will use inline integration to create your own FMU that will be used in cosimulation. Because understanding FMIs is way too out of scope for this course, we have provided you with a template FMU that you can fill with your own code. In short, an FMU is a zipped package of a specifically structured set of files. The template file can thus be uncompressed and you can browse the internal files. For your convenience, the only files you need to edit are:

• sources/defs.h: Header file that uses #defines for human-readable coding.
• sources/eq0.c: Should contain the equations for the first iteration.
• sources/eqs.c: Should contain the equations for all other iterations.
• modelDescription.xml: Contains the simulation setup and the exposed interface of the FMU.

sources/defs.h

When you start the process of inline integration, you will get M internal variables. These are the outputs and inputs of each of the blocks, after the CBD is flattened. The variables allow you to easily write the equations. For instance, an AdderBlock with two inputs will yield three variables: IN1, IN2 and OUT. For the sake of convenience, we will expose all variables to the user. This allows us to store them as an array (called modelData) in the CBD struct (see sources/model.c). To allow user-readability, you should use the sources/defs.h file to predefine some variables as macros. For the AdderBlock, an example can be:

				#define M 3  /* Number of Internal Variables */
				#define _Adder_IN1 (cbd->modelData[0])
				#define _Adder_IN2 (cbd->modelData[1])
				#define _Adder_OUT (cbd->modelData[2])
				

Notice that M is required to be set to the total number of variables. Additionally, you can assume they have access to a variable called cbd, which is a pointer to the current instance of the CBD struct. Furthermore, you may notice that the C-code has no knowledge of the past (i.e., the modelData array does not contain timing information). This is especially useful to know when dealing with DelayBlocks. For the macro names, you are encouraged to use (and encode) the individual block names. This way, you can understand, debug and verify your own generated code.

sources/eq0.c and sources/eqs.c

The normal inline integration will be yielded in these files. You have access to a variable called cbd, which is a pointer to the current instance of the CBD struct. You should only need its internal modelData array (see above). Additionally, there is a variable called delta, which is the time since the previous computation. For the first iteration, this will be set to 0. For the adder, both the sources/eq0.c and sources/eqs.c will look as follows:

				_Adder_OUT = _Adder_IN1 + _Adder_IN2;
				

modelDescription.xml

This file contains a lot of information about the general simulation setup. Please only change the described tags of this file. The other files are constructed such that they are linked correctly. Note that the description below and the provided FMU are a massive oversimplification of the total FMI 2.0 reference.

• <ModelVariables>: Contains the M internal variables, denoted as unique <ScalarVariable>s.
• <ScalarVariable>: Has 5 attributes:
  • name (the name of the variable, as exposed to the user),
  • valueReference (the index in the modelData array that refers to this variable),
  • initial (how the variable is initialized; "exact" implies that it is known, "calculated" otherwise) (this attribute may not be set for inputs!),
  • causality (kind of parameter; "input" implies model input, "output" is model output, "local" otherwise),
  • variability (how the value will change; "constant" for ConstantBlocks, "continuous" otherwise).
  Additionally, there is a <Real> tag as a child, which has a start attribute (indicating the starting value of the variable) if the <ScalarVariable>'s initial equals "exact", or when it is of causality "input" and variability "continuous".
• <ModelStructure>: Special cases of the variables. It contains an <Outputs> tag, which has a set of <Unknown> children. There is also an <InitialUnknowns> tag, which will be the exact same as the <Outputs> tag (for the purposes of this exercise). Each <Unknown> has an attribute index, referring to the corresponding (one-based) index w.r.t. the ordering of the ScalarVariables. Only the model outputs are <Unknown>s (for the purposes of this exercise).

An example for the adder is given in this file.

Tasks

In the previous exercise, you helped to keep the president safe by implementing and tuning a PID controller, by having both a plant and controller in Modelica. In reality, you would like to have a controller algorithm that you can embed on a real system. The plant is given as an FMU. It will look a bit like the following figure (notice the correspondence with the previous assignment):

1. Go to Modelica and open your previous solution for the Plant-Controller setup. Now, remove your PID components and mark $u(t)$ as the input (called u) and $e(t)$ as the output (called e). Make sure you don't remove your solution for the previous assignment if you did not submit this yet!
   Ensure that the target car's position can be accessed by using forward_car.y!
   Go to Tools > Options > FMI and make sure that Version is 2.0, Type is "Model Exchange and Co-Simulation" and Platforms is "Static". Now, right-click on your Modelica class that corresponds to the plant and environment and select Export > FMU (this might take a few seconds). The output panel will show where the FMU has been generated, so you can access it. Rename it to "Plant.fmu".
2. Create a PID controller using CBDs. Use ConstantBlocks for $K_p$, $K_i$ and $K_d$, which should (by default) have your solution from the last exercise. Note that this will hardcode these values in your solution. We will do this for the sake of simplicity. The controller's input (called IN) will be given the error value $e(t)$, as defined in the previous assignment. Its output (called OUT) will yield $u(t)$. Include a graphical representation of this CBD model in your report.
3. Create a script that generates the required code for the sources/eq0.c, sources/eqs.c, sources/defs.h files, as well as the required tags for the modelDescription.xml file. The CBD simulator comes with a lot of features to help you accomplish this task. Take a look at the depGraph.createDepGraph() and scheduling.TopologicalScheduler.schedule() methods (as well as table 2 of this paper) to help you get started. You may flatten your model (CBD.flatten() or CBD.flattened()) for a simpler solution, but you are allowed to keep hierarchy for a smarter solution. Your choice will be reflected in your grade.
   The input and output ports of the PID FMU are called PID.IN and PID.OUT, so you do not have to remove prefixes after flattening. While FMI has many ways of solving integrals, you are not allowed to use any of these.
   You may use Jinja2 or string manipulation for the code generation.
4. Generate the code for the PID controller, using your script. Copy-paste it into the correct files. Compress your FMU using the zip-algorithm and rename it "PID.fmu".
5. Install FMPy and download the compile_and_run.py orchestration script.
6. Run the orchestration script. It will first compile your created FMU and yield any errors it occurs. Next, it links both FMUs together and simulates them. The resulting plot should be the same as the ones obtained from Modelica. Include both in this report and show/prove they are the same.

Please contact the TA if you have any questions, or get stuck whilst working on this assignment. The goal of this task is to help you realize how CBD models can be converted into executable code, how it is done in industry. You are not required to understand all FMI details, or to get stuck in debugging C-code.

• 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 or here (mirror).
• The files for DrawioConvert can be found here:
  • Small Tutorial
  • DrawIO Library
  • Download Python Code
  • Requires: Jinja2
• Template FMU (empty): http://msdl.uantwerpen.be/people/hv/teaching/MoSIS/assignments/CBD/template.fmu
• FMU Co-Simulation Orchestration Script: http://msdl.uantwerpen.be/people/hv/teaching/MoSIS/assignments/CBD/compile_and_run.py
• Useful documentation:
  • The CBD Simulator's documentation
  • CBDs: http://msdl.uantwerpen.be/people/claudio/pub/Gomes2016a.pdf
  • C programming: https://www.w3schools.com/c/index.php
  • FMPy: https://pypi.org/project/FMPy/