• Due Date: 29 October before 23:55.
• 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 must be called index.html. Make sure to mention the names and student IDs of all the team members. The other team member must submit a single HTML file containing only the coordinates of both team members. 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 use the Modelica Language to automate and optimize a driver-less PRT controller. You will learn how Ordinary Differential Equations (ODE) can be used to facilitate the modeling and simulation of Cyber-Physical Systems (CPS). This assignment expands on the PRT case from previous assignment. This solution will be made with OpenModelica. Make sure to clearly indicate how your solution can be executed and which files implement which models in your report.

On MAC, you can install OMEdit by installing the most recent binary.

To make you familiar with Modelica and OpenModelica's interface, let's take a look at a simple example: Predator-Prey models.

Lotka-Volterra equations are first-order nonlinear, differential equations that can be used to describe the dynamics of biological systems, specifically the relationship between predators and preys in a specific biotope. The equations are as follows: $$\begin{cases} \dfrac{dx}{dt} &=& \alpha x - \beta xy\\ \dfrac{dy}{dt} &=& \delta xy - \gamma y \end{cases}$$ Where $x$ identifies the number of prey and $y$ the number of predators. Their derivatives are therefore the instantaneous growth rates of both populations. $\alpha$ and $\beta$ represent the growth and death rates of the preys, whereas $\delta$ and $\gamma$ identify the growth and death rates of the predators.

You will need to perform the following tasks step by step:

1. Use OMEdit's Text View to create a Predator-Prey model for a fictious fox (predator) and hare (prey) population. Assume the biotope starts with 10 foxes and 10 hares. Additionally, $\alpha = 1.5$, $\beta = 0.7$, $\gamma = 0.2$ and $\delta = 0.2$.
2. Simulate your model for 100 seconds (in reality, these "seconds" corrspond to "years") and describe what happens to the system.
   • In the "Simulation Setup", change the "Method" to "rungekutta" and keep the "Tolerance" at 1e-6 to get the best results.
   • Include a plot of the population dynamics w.r.t. the time.
   • Include a phase-space plot (predator count vs prey count) by using OMEdit's "Parametric Plot".
   • Does this system stabilize? Why (not)? How can you tell?
   • What happens when you change the "Integration" parameters of the "Simulation Setup"?
3. Make $\alpha$, $\beta$, $\gamma$ and $\delta$ variable parameters for the user to manipulate. Also, add an output port for both the predator count and the prey count.
4. Recreate the same Predator-Prey model, using OMEdit's Diagram View and basic blocks. You don't have to add parameters for this model. Include a plot and compare your results to the previous results.

In a PRT system, it is ideal to have as much automation as possible. Hence, we wish to create an automated controle system for a driver-less trolley. This controller should replace a human driver.

Because of the limited person capacity and the short distances that need to be traveled, we will assume there are no seats in the rail cart. Instead, there are multiple poles on which commuters can secure themselves while standing. In a typical execution, the trolley receives the desired velocity from a centralized coordinator, based on the trolley's current GPS signal. Where a human driver would have to manually adapt the trolley's velocity based on experience, the automated controller will use this received data to accelerate, decelerate or brake the rail cart.

For the purposes of this assignment, you can assume the PRT moves in one dimension (i.e., the rail is straight and flat) and only has the capacity to transport a single passenger. The centralized coordinator can thus be represented by a pre-defined ideal velocity profile. The velocity is given in meters per second ($m/s$), always positive (or zero) and cannot vary faster than $10\ m/s^2$. For instance, the controller will never be asked to accelerate from $10\ m/s$ to $30\ m/s$ in one second. The plot below shows an example ideal velocity profile.

The control system must ensure a correct and optimal execution. This implies:

correct

The trolley will attain the ideal velocity, given enough time. On top of that, the passenger will not fall due to acceleration/deceleration.

optimal

The trolley will attain the ideal velocity as fast as possible.

Notice the trolley may overshoot/undershoot the desired velocity. A drawing of this system is given in the figure below.

As you can guess, the control system should output $F_{traction}$, based on its input: the difference between the ideal velocity $v_{ideal}$ and the actual velocity $v$. If $v_{ideal} - v$ is small, the trolley is close to the desired velocity, so $F_{traction}$ should be small as well. If $v_{ideal} - v$ is large, $F_{traction}$ should be large as well, to allow the trolley to quickly accelerate/decelerte towards its ideal velocity. A simple implementation of such a system would be a so-called Bang-Bang controller. Two thresholds are set: $d_{min}$ and $d_{max}$. If $v_{ideal} - v$ exceeds $d_{max}$, then $F_{traction}$ is set to some positive value $g$. However, if $v_{ideal} - v < d_{min}$, then $F_{traction}$ is set to $0$ (which allows the car to decelerate due to friction and air resistance). In the real world, the passenger would feel the sudden bursts of acceleration, hence the name: "Bang-Bang".

A slightly better controller is a Proportional-Integral-Derivative (PID) Controller. Instead of setting $d_{min}$ or $d_{max}$, $F_{traction}$ will be related to $v_{ideal} - v$. This relationship is obtained by computing the sum of three different controllers:

Proportional Controller:

For some constant value $K_p$, this controller sets $F_{traction} = K_p\cdot (v_{ideal} - v)$.

Integral Controller:

For some constant value $K_i$, this controller sets $F_{traction} = K_i\cdot \int (v_{ideal} - v) dt$.

Derivative Controller:

For some constant value $K_d$, this controller sets $F_{traction} = K_d\cdot \dfrac{d(v_{ideal} - v)}{dt}$.

In total, the PID-controller sets: $$F_{traction} = K_p\cdot (v_{ideal} - v) + K_i\cdot \int (v_{ideal} - v) dt + K_d\cdot \dfrac{d(v_{ideal} - v)}{dt}$$ Intuitively, the derivative component analyzes the speed at which $v_{ideal} - v$ evolves. It tries to predict the future state of the system and smoothens $F_{traction}$ accordingly. Yet, using a derivative causes an unstable and inconsistent signal. The integral component counteracts this instability by accumulating $v_{ideal} - v$ over time. Finally, the proportional component identifies the allowed rate of change.

While we can model the controller as a PID-controller, we still need access to a trolley in order to validly analyze the system's behaviour. Because we do not have access to an actual trolley, we need to create a model for this as well. Such a model is called a plant in control systems' jargon. To do this, we will use a set of Orinary Differential Equations (ODE) that represents the behaviour of the trolley. A mathematical derivation for the following ODE can be found here. $$\begin{cases} \dfrac{dv_{psgr}}{dt} &=& \dfrac{k\cdot (-x_{psgr}) + c\cdot (-v_{psgr}) - m_{psgr}\cdot\dfrac{F_{traction} - \dfrac{1}{2}\cdot p\cdot v_{trolley}^2\cdot C_D\cdot A}{m_{trolley} + m_{psgr}}}{m_{psgr}}\\ \dfrac{dv_{trolley}}{dt} &=& \dfrac{F_{traction} - \dfrac{1}{2}\cdot p\cdot v_{trolley}^2\cdot C_D\cdot A}{m_{trolley} + m_{psgr}}\\ \dfrac{dx_{psgr}}{dt} &=& v_{psgr}\\ \dfrac{dx_{trolley}}{dt} &=& v_{trolley} \end{cases}$$

You will need to perform the following tasks step by step:

1. Implement a look-up model in Modelica. This model will continuously output the ideal velocity $v_{ideal}$, given the current time (hint: Modelica stores the time in the time variable). The look-up table is given below. Simulate this model for $300$ seconds and create a plot of the result. Note that this plot will differ from the plot given above. You are free to choose the interval size, but make sure this is small enough.
   
   

   time (s)	output (m/s)
   < 10	0
   < 170	10
   < 200	8
   < 260	18
   ≥ 260	12

   
   
2. Implement a Modelica model for the given plant model. The model should have 1 input ($F_{traction}$) and 2 outputs ($v_{trolley}$ and $x_{psgr}$). Make sure all other constants ($m_{psgr}$, $m_{trolley}$, $k$, $c$, $C_D$, $p$ and $A$) can be set by the user (i.e., define them as model parameters).
3. Implement a simple Bang-Bang controller in Modelica. This has a single input $u$ and a single output $F$. Allow users to set $d_{min}$, $d_{max}$ and $g$. Next, implement the full PRT system as shown below:
   
   
   To ensure somewhat realistic results, the following set of constants can be used: $$\begin{cases} m_{psgr} &=& 77\ kg\\ m_{trolley} &=& 2376\ kg\\ k &=& 300\\ c &=& 150\\ C_D &=& 0.6\\ p &=& 1.2\\ A &=& 9.12\\ d_{min} &=& -1\\ d_{max} &=& 1\\ g &=& 2000 \end{cases}$$
4. Simulate the combined model for $300$ seconds and plot the results. You are free to choose the interval size, but make sure this is small enough. Compare people displacement ($x_{psgr}$) against the trolley's acceleration. Does the passenger fall or not? You can assume that they fall if the displacement exceeds $0.35\ m$. Additionally, plot the $F_{traction}$ profile for this trajectory.
5. Run a number of experiments for multiple values of $d_{min}$, $d_{max}$ and $g$. What happens to your results? Can you find the "meaning" of these parameters? Make sure to include your plots for all your experiments in your report. What happens if you change the if-statements to when-statements, or vice-versa? To ensure the required behaviour, you may need to change more than a single keyword!
6. Use a Modelica PID-controller (Modelica > Blocks > Continuous > PID) instead of the Bang-Bang controller. To ensure realistic results, the following set of constants can be used: $$\begin{cases} m_{psgr} &=& 77\ kg\\ m_{trolley} &=& 2376\ kg\\ k &=& 300\\ c &=& 150\\ C_D &=& 0.6\\ p &=& 1.2\\ A &=& 9.12\\ K_p &=& 400\\ K_i &=& 1\\ K_d &=& 0 \end{cases}$$ Note that OpenModelica's PID controller uses parameters $k$, $T_i$ and $T_d$ instead of $K_p$, $K_i$ and $K_d$. Their relationships are defined as follows: $$\begin{cases} K_p &=& k\\ K_i &=& K_p / T_i\\ K_d &=& K_p \cdot T_d \end{cases}$$ Now, simulate the newly combined model for $300$ seconds and plot and discuss your results.
7. Recall that the controller should be correct and optimal. Additionally, we would like to ensure the passengers do not fall. This can be done by tuning the PID-controller. In other words: find the ideal $K_p$, $K_i$ and $K_d$ values, either manually or automatically. This can be done as follows:
   1. Set some initial values for $K_p$, $K_i$ and $K_d$.
   2. Simulate the model.
   3. Evaluate the results. This can be done by sending $x_{psgr}$, $v_{ideal}$ and $v_{trolley}$ to a Cost Function Block, as shown in the figure below. The higher the cost, the worse the selection for $K_p$, $K_i$ and $K_d$. Define a simple cost function and plot it for the initial experiment and your found solution. Make sure to discard any result that causes the passenger to fall.
      
   4. Alter the values for $K_p$, $K_i$ and $K_d$.
   5. Repeat steps 2-4 until a result has been found.
   Hint: $K_p$, $K_i$ and $K_d$ respectively lie in the ranges $[300, 350]$, $[0.5, 1.5]$ and $[-20, 70]$.
   Note: You are free to chose manual and/or automated tuning for this task, but know that automated tuning is much more sophisticated. This will be reflected in your grade.
8. Write a report that explains your solution for this assigment. Include your plots and discuss them. Also clearly indicate which models can be found in which files. Make sure to mention all your hypotheses, assumptions and conclusions. See also the "submission information" at the top of this page.

Automated PID tuning in OpenModelica is, unfortunately, not straightforward. You will have to create a script that executes your model multiple times, changing some parameters for every run. Such a collection of simulations is called a simulation sweep or a batch simulation. As this is not a standard feature in OpenModelica, you will have to obtain the results in a somewhat roundabout way. What follows is a description to do a simulation sweep of an OpenModelica model from a Python script.

1. For simplicity, it is encouraged to set $K_p$, $K_i$ and $K_d$ as your model parameters. Even if you have built your model in the Diagram View, you can still alter the Text View for this. Don't forget to link the parameters correctly to your code!
2. Open the Simulation Setup and navigate to the Output tab. Change "Output Format" to "csv".
3. Compile/simulate your model in OpenModelica. This will generate an executable (and a set of config files) from which the results can be obtained. Luckily, this executable also provides a way to alter the model parameters without needing to recompile the model every time.
   • Open Tools > Options and open the General menu. The "Working Directory" field identifies the location on your system where these files will be generated.
   • If these directories are empty, reopen OMEdit and re-compile/re-simulate your model once more. To disable the removal of these files, follow the steps described in "Practical Issues" at the bottom of this page.
4. Open the *_init.xml file and make sure the outputFormat attribute is set to "csv".
5. You can now call your model from Python using the os module. The example below assumes the working directory is stored in the path_to_exec_dir variable and the model is called "MyModel". Obviously, non-Windows users should remove the ".exe" from the following code.

                   import os
                   os.chdir(path_to_exec_dir + "/MyModel")
                   os.system("MyModel.exe")
                   results_file = "MyModel_res.csv"
                 

   The results can now be read using Python's Pandas module or any other data processing module you desire.
6. The executable has a lot of model parameters. For the sweep, you can use the -override flag as follows:

                   os.system("MyModel.exe -override param1=value1,param2=value2,param3=value3")
                 

Take pride in your work. Make sure your report and all models and images are clearly readable. Additionally, produce clean and (preferrably) well documented code. If a plot is difficult to read, also include the source files from which these can be (re-)generated (*.csv, *.mo). Do not submit your executables! You will not lose points for doing too much, but you will lose points if you do too little.

• This assignment is to be created with OpenModelica and Python 3.6+. Do not use features from Python 3.10 (current alpha version), nor from Python 2.7 (discontinued as of January 1, 2020).
• $\LaTeX$ rendering in HTML can be done via MathJax. Alternatively, it is allowed to refer to *.tex and/or *.pdf files in your report.
• It is allowed to have unresolved "warnings" in your solution, but be aware they might yield different results on other machines.
• While OpenModelica allows the creation of icons for your models, this is not required for the given assignment.
• Because the relationship between ($K_p$, $K_i$, $K_d$) and ($k$, $T_i$, $T_d$); you are free to chose which triplet you tune. Keep in mind that $T_i$ will be inverted and therefore cannot be $0$.
• It is possible that OpenModelica removes the generated files and executables when closed. Open Tools > Options and navigate to the Simulation menu. At the bottom, disable "Delete entire simulation directory of the model when OMEdit is closed" to prevent this from happening.
• On MAC, you can install OMEdit by installing the most recent binary.
• Useful links:
  • Modelica: https://modelica.org/
  • OpenModelica: https://www.openmodelica.org/.
  • PID-Controller: https://en.wikipedia.org/wiki/PID_controller.
  • Modelica executable command line arguments: https://openmodelica.org/doc/OpenModelicaUsersGuide/latest/simulationflags.html.