Modelica Assignment 

Practical Information

  • Due Date: 2 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 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 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.

Goals

In this assignment, you will use the Modelica Language to create a simple motor yacht autopilot. You will learn how Ordinary Differential Equations (ODE) can be used to facilitate the modeling and simulation of Cyber-Physical Systems (CPS). Furthermore, you will learn how a generic model must be calibrated to fit real-world data and how it can be tuned to optimize some goal function.

This solution will be made with OpenModelica. MacOS users can install OpenModelica from the binaries (make sure to use the latest version). Also ensure that you clearly indicate how your solution can be executed and which files implement which models in your report!

Assignment Overview

This assignment consists of four main parts:
  1. Plant Model Creation
  2. Plant Model Calibration
  3. Controller Model Creation
  4. Controller Model Tuning

Plant Model Creation

In this part of the assignment, we will create a simple ODE for modelling the yacht movement, based on nautical physics. This ODE will be controlled such that the yacht may use cruise control. For the sake of simplicity, we shall represent this movement in one dimension (i.e., no steering shall be done). In figure 1, you can see a free body diagram of this system. The yacht moves towards the right (positive direction).



Figure 1: Free body diagram of the yacht.

The yacht moves forward with a throttle force of $\vec{F}_T$. Additionally, it receives some resistance $\vec{F}_R$ from the water. Hence, the total force can be computed as $F = F_T - F_R$. Newton's second Law of Motion states $F = m\cdot a$, with $m$ the mass of the yacht and $a$ its acceleration. This implies (if the mass does not change): $$\begin{align} F = m\cdot a = m\cdot \dot{v} = m\cdot \dfrac{d v}{dt} &= F_T - F_R\\ \dfrac{d v}{dt} &= \dfrac{F_T - F_R}{m}\\ \end{align}$$ This results in an Ordinary Differential Equation (ODE) Initial Value Problem (IVP). An IVP typically has a time-based derivative of a variable $x$ on the left hand side (LHS) and a function in terms of that $x$ on the right hand side (RHS). Given an initial value $x(0)$, it is sometimes possible to analytically compute a solution for such an equation. Luckily, in the case no analytical solution may exist, a Modelica-based numerical simulation tool such as OpenModelica is able to provide a good approximation of the actual function. The above equation can be represented in Modelica and solved numerically. However, we do need more information abot $F_R$ and $F_T$.

From nautical physics, it is known that $$\begin{cases} F_R &= \dfrac{1}{2}\cdot\rho\cdot v^2\cdot S\cdot C_f\\ C_f &= \dfrac{0.075}{\left(\log_{10}(Re) - 2\right)^2}\\ Re &= \dfrac{v\cdot L}{k} \end{cases}$$

Where $\rho$ identifies the density of the water, $v$ the yacht's velocity, $S$ the submerged area of (the full hull of) the yacht, $C_f$ the friction value, $Re$ Reynolds number, $L$ the length of the yacht and $k$ the kinematic viscosity of the water. Most of these parameters are known. The waters in which the yacht is sailing has salt water of about $15^{\circ}C$. This implies $\rho = 1025\ kg/m^3$ and $k = 1.188\cdot 10^{-6}\ m^2/s$. The yacht's length $L$ is $21.54\ m$ and its dry mass $m$ is $32,000\ kg$.

Tasks

Your task is to model the yacht's plant (i.e., the equation(s)) in Modelica.

  1. Implement the given plant equation(s).
    Hint: you can use the der(v) function to represent the time-based derivative of the velocity.
  2. Provide an output for the yacht's velcocity $v$ and its position $x$. Note that $\dfrac{d x}{dt} = v$. Assume the starting value for the velocity is $0\ m/s$. The starting position is $0\ m$.
  3. Make $\rho$, $L$, $S$, $m$, $k$ and $F_T$ parameters of the plant model and set these by default to the values given above. Set (for now) $S$ to $100\ m^2$ and $F_T$ to $400\ N$.
  4. Note that your simulation will crash when the velocity becomes $0\ m/s$ (because $\log_{10}(0) = -\infty$). Add a check in your plant that sets $C_f = 0$ if your velocity is $0\ m/s$.
  5. Simulate this model for $400$ seconds and provide three plots: one of the velocity, one of the acceleration and one of the position. Discuss the results in your report.
  6. Do not forget to discuss this task in your report.

Plant Model Calibration

While this plant model may appear to be realistic, we have made an assumption on the value of $S$. Truth be told, we have no idea about what is going on underneath the surface of the water. Even if we would know the shape of the hull, there is no guarantee on how deep the yacht lies in the water (this is also known as the draft of the yacht). All we do know for certain is that $S$ lies between $150\ m^2$ and $300\ m^2$. We do have a physical yacht and thus we are able to run some experiments, as shown in figure 2.



Figure 2: Two experiments done with the yacht.

The first experiment is an acceleration from approximately $0\ m/s$ to just over $1\ m/s$, using a throttle force of $400\ N$. Every $5$ seconds, we looked at the GPS-based speedometer on board and we have written down the value. The experiment results can be found in this CSV file. Similarly, we have released the throttle ($F_T = 0\ N$) at $10\ m/s$ and wrote down the values for the next $400$ seconds. These results can be found in this CSV file. Note that the speedometer is not fully accurate, resulting in quite some noise on the data.

Calibration is the process of estimating parameters of a model. Those parameter values are chosen which result in the closest match between simulated and measured values of pertinent, measurable variables, over time, over an entire experiment run. Such a match is often quantified in terms of a sum of squares of differences between measured and simulated values. In this case, we would like to estimate the submerged area of the yacht, the parameter $S$. An exhaustive way of finding optimal parameter values (note that in this case, there is only one parameter) is to loop over possible values for the parameters to find the one that gives the closest match mentioned above. This is known as a parameter sweep.

Tasks

The main purpose of this part is finding the submerged area of the yacht $S$, based on experimental data (i.e., observations/measurements on the real system).

  1. 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 (upon closing OMEdit), follow the steps described in "Practical Issues" at the bottom of this page.
    • You can re-simulate your model (without recompiling) from a terminal by running the executable in the "Working Directory". You can even change the Simulation Setup by using some command line parameters. For instance, to change the parameter $S$ to $789\ m^2$ and $\rho$ to $988\ kg/m^3$, you can use the following command (non-windows users should obviously forgo the "exe" extension): 
      	YachtPlant.exe -override S=789,density=988
  2. Create a script that is able to do a parameter sweep of the plant, changing the value of $S$ for each simulation run and storing the results for future analysis. You are free to choose which language to use for this script. We recommend using Python or MOS.
    • If you decide to use MOS, OpenModelica's Scripting Language, you create a *.mos file in the "Working Directory" that looks like this. Next, in OpenModelica, you navigate to Tools > OpenModelica Compiler CLI and call runScript("file.mos").
    • Make sure to store all experiment results, as they are needed in the next step.
    • By default, the generated results will be *.mat files. You can change this to *.csv (note: this will be slower) by changing the outputFormat attribute to "csv" in the *_init.xml file. Alternatively, if you use MOS, you can use the outputFormat parameter of the buildModel function. The *.mat files can be read in Python using SciPy: 
      	from scipy import io
	resDict = scipy.io.loadmat("YachtPlant_res.mat")
    • You can either brute-force over all possible values of $S$, but you may also use a smarter algorithm to find the best value faster (see also the next task).
    • Make sure to create a plot that shows some of these experiments (e.g., velocity as function of time) results!
  3. Compute the sum of squared errors between each of your experiments/simualtions and each of the given traces. The simulation that has the smallest total error is the simulation that has the most accurate value for $S$.
    • You may compute this directly from the parameter sweep script, but this will be generally slower.
    • You are free to chose how you compute this; whether you use Python, Matlab, Excel...
    • $S$ will be an integer value. There is no need to concern yourself with decimals.
    • Note that the given experiments only have datapoints every $5$ seconds!
    • Do not forget to find an ideal $S$ such that both the acceleration and deceleration profiles match!
  4. Do not forget to discuss this task in your report. Include plots ( sum of squared errors as a function of parameter values) and experiment results!

Controller Model Creation

In this part of the assignment, we will create the yacht's cruise control system. In a typical execution, the yacht receives the desired velocity from a centralized coordinator, based on its current GPS signal. Where a human skipper would have to manually adapt the yacht's velocity based on experience, the automated controller will use this received data to accelerate or decelerate.

The centralized coordinator is represented by a pre-defined ideal velocity profile. The velocity $v$ is given in meters per second ($m/s$) and always positive (or zero). Additionally, in the bay where the yacht is sailing, it is illegal to go faster than $20\ knots$ (which is approximately $10.28\ m/s$). From practice, we know that the yacht's engine cannot provide a speed variation faster than $1.5\ m/s^2$. For instance, the controller can never be asked to accelerate from $1\ m/s$ to $3\ m/s$ in one second. While this is mathematically possible in our system, it will be impossible in reality. Figure 3 shows an example velocity profile.



Figure 3: Example ideal velocity profile of the yacht.

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

correct
The yacht will attain the ideal velocity, given enough time. On top of that, the velocity must not exceed $10.28\ m/s$ and the absolute acceleration must be at most $1.5\ m/s^2$.
optimal
The yacht will attain the ideal velocity as fast as possible.
Notice the yacht may overshoot/undershoot the desired velocity. A plot of an example controller is shown in figure 4.



Figure 4: Velocity of the yacht using an automated (Bang-Bang) controller.

As you can guess, the controler has a single output $F_T$; and a single input $u$ (which is the difference between the ideal velocity $v_{ideal}$ and the actual velocity $v$). Ideally, we would like $v_{ideal} - v$ to be small, as this implies the yacht is close to the desired velocity. A simple implementation of such a system would be a so-called Bang-Bang controller. Here, an upper bound $d_{max}$ and a lower bound $d_{min}$ is set around the $v_{ideal}$. When $v$ exceeds (from below) $v_{ideal}+d_{max}$, the acceleration stops, such that the ship can decelerate due to the water resistance. Alternatively, some implementations help this by providing a "braking" force (for a ship this happens by making the engine propeller(s) turn backwards). But when $v$ goes below (from above) $v_{ideal}-d_{min}$, the ship needs to accelerate again, with a given (acceleration) force $g$. In the real world, the passenger could feel the sudden changes in acceleration, hence the name: "Bang-Bang".
Notice the "when" in the previous sentences. This implies that an action should happen exactly when a condition becomes true. The usage of "if" statements will only look the next computation if the condition is true.

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

Proportional Controller:
For some constant value $K_p$, this controller sets $F_T = K_p\cdot (v_{ideal} - v)$.
Integral Controller:
For some constant value $K_i$, this controller sets $F_T = K_i\cdot \int (v_{ideal} - v) dt$.
Derivative Controller:
For some constant value $K_d$, this controller sets $F_T = K_d\cdot \dfrac{d(v_{ideal} - v)}{dt}$.
In total, the PID-controller sets: $$F_T = K_p\cdot (v_{ideal} - v) + K_i\cdot \int (v_{ideal} - v) dt + K_d\cdot \dfrac{d(v_{ideal} - v)}{dt}$$

Tasks

Model the yacht's controller in Modelica.

  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 $600$ seconds and create a plot of the result. Note that this plot will be the same as the plot given above. You are free to choose the interval size, but make sure this is small enough.

    time (s)output (m/s)
    < 100
    < 2105
    < 30010
    < 4404
    < 5208
    ≥ 5203

  2. Implement a simple Bang-Bang controller in Modelica. This has a single input $u$ and a single output $F_T$. Allow users to set $d_{min}$, $d_{max}$ and $g$. Assume $d_{min} = -1$, $d_{max} = 1$ and $g = 12000\ N$.
  3. Implement the full plant-controller system as shown in figure 5, by combining all previous models.


    Figure 5: Plant and controller system for the yacht.
  4. Simulate the full model for $600$ seconds and plot the results ($x$, $v$ and $a$). You are free to choose the integration method and interval size. Check that the acceleration and velocity stay within the given ranges. Additionally, plot the $F_T$ 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.
  6. Use a Modelica PID-controller (Modelica > Blocks > Continuous > PID) instead of the Bang-Bang controller. You can use $K_p = 4000$, $K_i = 10$ and $K_d = 0$.
    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 $600$ seconds and plot and discuss your results (compared to the previously obtained results).
  7. Do not forget to discuss this task in your report.

Controller Model Tuning

Similar to the calibration of the plant, the controller also needs to be tuned for the most ideal solution. Instead of fitting a curve to existing experiment data, we shall try to find the controller parameters such that the given velocity profile is followed as closely as possible. Figure 6 shows an adaptation of the previous plant and controller system that also allows the computation of a cost.



Figure 6: Plant and controller system for the yacht with cost computation.

Tasks

Tune the yacht's controller.

  1. Implement a cost model in Modelica. You are free to choose your own cost function, but clearly describe what you chose and why.
  2. Implement the full plant-controller system as shown in figure 6.
  3. Create a tuning script using the same logic as defined in the Plant Model Calibration section. This script should tune the PID-controller. 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! In order to provide some guidance, you are given the following ranges: $K_p \in [6000, 7000]$, $K_i \in [300, 350]$ and $K_d \in [-1, 1]$. Furthermore, $K_p$, $T_i$ and $T_d$ will be integers.
  4. Find the ideal $K_p$, $K_i$ and $K_d$ values and show a plot for $v$, $x$, $F_T$ and $a$. Make sure that the speed limit is respected and that the absolute acceleration is at most $1.5\ m/s^2$.
  5. Do not forget to discuss this task in your report! Clearly indicate which values you have obtained for $K_p$, $K_i$ and $K_d$.

Practical Issues
Maintained by Hans Vangheluwe. Last Modified: 2023/10/24 15:05:43.