• Due Date: 23 October 2026, before 23:59.
• Team Size: 2 (pair design/programming)
  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 accompagnied 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 should 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 use the Modelica Language to first create a model of a Cyber-Physical System (CPS) and then create a controller for optimal control of the CPS. You will learn:

1. How Ordinary Differential Equations (ODE) can be used for the modeling and simulation of a CPS.
2. Furthermore, you will learn how a generic model must be calibrated to fit real-world data.
3. You will then design and model a controller to control the CPS.
4. You will finally tune the control parameters to optimize some goal function based on the (simulated) behavior of the CPS.

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) individual oral evaluation, which will typically be conducted in the week following the deadline of every assignment.

You will use OpenModelica to perform all the tasks. You can install the latest stable version of OpenModelica on your systems by following the directions by going to the OpenModelica website > Download > Windows/Linux/Mac (based on your OS). Choose the 'Official Release' package.

Installing OpenModelica installs a number of applications in your system:

OMEdit

The OpenModelica Connection Editor (OMEdit) is the standard graphical user interface that is used for graphical Modelica model creation and editing. We will use OMEdit for all the tasks. Launch instructions and other features of OMEdit are described on the OMEdit website. OMEdit has 4 different views but we are only interested in the Modelling and Plotting views. We will primarily use OMEdit to design the model. Information about the other OpenModelica applications are just given for information.

Modelica Standard Library

A major advantage of using Modelica is the availability of a wide variety of Modelica libraries, making it easier to re-use blocks and classes from these libraries. For the assignment, you should ensure, once you have installed OpenModelica, that the Modelica Standard Library (MSL) package is also correctly downloaded, and automatically loaded when you start OMEdit. To do so, go to File > Manage Libraries > Install Library > select Modelica with Version 4.0.0 > OK. The Modelica Standard Library should be visible when you restart OMEdit.

OMNotebook

OMNotebook is a notebook-like tool to learn the usage of OpenModelica and construct more-readable documentation which incorporates code-windows in between text. Think of Jupyter Notebooks.

OMShell

OMShell is an interactive session handler that parses and interprets commands and Modelica expressions for evaluation, simulation, plotting, etc. The session handler also contains simple history facilities, and completion of file names and certain identifiers in commands. We will not use OMShell because we want to learn and visualize graphically the models that we are building. Hence we will only use OMEdit.

We will use Python because we can then directly use the data obtained from simulations for analysis through a number of data analysis packages available in Python.

One thing to note is that a Modelica model needs to be compiled only once. Once compiled, the parameters can be changed in subsequent executions/simulation runs without having to re-compile the model. This is advantageous since compilation can prove to be a significant overhead if we are optimising for a number of different parameters and ranges (read: the number of simulation runs could go into the tens of thousands or more, when using simulation for the purpose of optimization).

Basics

While the Modelica language and its capabilities of object-oriented acausal modeling of the multi-physics have already been described in the lectures, what follows is a 'how-to' to get started building your own Modelica models for the assignment using OMEdit.

On the left of the OMEdit window is the Libraries pane which displays all the libraries that have been loaded in the current Modelica session. Modelica models are organised in Packages which can contain other packages and Modelica classes.

Create a Modelica Package

Figure 1: Example package creation.

Go to File > New > New Modelica Class > in the Create New Modelica Class dialog box > enter a valid package name > for Specialization, select Package, this is also shown in Fig.1 > OK. You should now see your Modelica package in the Libraries pane. You should save it before proceeding further by clicking on the corresponding icon in the menu bar or going to File > Save.

If you have already created and saved a package, you can load it in the current session by going to File > Open Model/Library Files, and choosing the corresponding file in the file browser.

Create your first Modelica Model

Right click on your Modelica package in the Libraries window > in the drop-down menu > select New Modelica Class > in the create New Modelica Class dialog box > enter a valid class name > for Specialization, select Model > for Insert in class ensure the textbox specifies the Modelica package within which you want to create your Modelica model > OK. You have created your first Modelica model! You should save it before proceeding further by clicking on the corresponding icon in the menu bar or going to File > Save.

Modelica classes can extend other Modelica classes, i.e. inherit the declarations from the extended class. They may also simply re-use components defined in other Modelica classes.

Build/Edit your Modelica Model

Double click on the Modelica model you want to edit. This takes you to the Modeling view. In the top right corner of the Modeling view are 3 buttons to toggle between different editing interfaces. The 2nd button is to open the graphical editor and the 3rd button is to open the textual editor.

Figure 2: Buttons to change editor views. The Icon view is only used to arrange connector points for a block.

The diagram view can only be used in case the model consists of blocks with causal/acausal connections between them.

Look at the Modelica By Example website for simple introductory examples of Modelica models based on simple Ordinary Differential Equations of physics.

If you are making a component-based model. You can add components to the model by dragging them from the Libraries pane onto the graphical editor. You will be prompted to give the component a unique name. You can draw connections between the interfaces of components by holding the left mouse button while simultaneously moving it from one interface ot the other. This highly intuitive. You can modify parameters of the components by double-clicking on them, which will open the Element Parameters box.

Simulate your Modelica Model in OMEdit

Figure 3: 'Check Model' and 'Simulation' buttons in OMEdit.

Once you have defined a complete and correct Modelica model, you can now simulate it. To do so, you should click on the S icon in the menu bar. This opens the Simulation Setup dialog box. You can set simulation parameters in this dialog box.

Make sure that you have verified the simulation start time, stop time, step-size/number of intervals, solver method (DASSL is preferred), and tolerance (the default 1e-6 is good enough, as low as possible is preferred ideally).

You can also set the output format to output the data to a file. However, we do not need to do so, since we will later perform simulations through Python where we will handle the simulation output directly. This is just for visualization and initial validation of the model.

Click OK to begin simulation.

Note: You can also 'check' the model by clicking on the tick-mark icon in the menu bar. The output from this can be observed in the message browser. The message browser will tell you if something is missing in the model, or if there are inconsistencies in the model.

Visualize simulation trace

Clicking OK in the Simulation Setup ideally takes you to the Plotting view. If not, you can go to the plotting view through the corresponding tab in the bottom-right corner.

On the right you will see the variable browser. You can select any variable/s to plot them in the pane and observe their values over the simulation. The values indicated in the variables pane is the value of that variable at the end of the simulation.

Parametrizing simulation quantities in OMEdit

To define a parameter you simply have to use the 'parameter' keyword followed by its unit name and default value in this syntax:

parameter <unit> <name> = <default value> ;

If you want to parametrize a quantity which is a parameter of a block in your model. Simply double-click on the block in the graphical editor and change the parameter value to the parameter name. This value need not be a numeric quantity, but it can be a reference to the name of the parameter within the model, or a function of the parameter as well. Once you have changed the value in the block, go to the textual editor of your model (not of the block!), and use the same syntax as above to define a parameter. Note that the name of this parameter should correspond to the name of the parameter you have defined in the value of the block. Look at the following image for an example:

Figure 4: Parametrizing the values of the block-parameters of the blocks in your model.

Multiple simulations and analyses of your Modelica model

The example Python file contained in the zip file describes how to use Python to simulate and visualize data from an example Modelica model which will be introduced in the tutorial as well.

We will perform multiple simulations by generating the executable through OMEdit. You should make sure that the executable is generated and not deleted in the right directory by:

1. go to Tools > Options > General > Adjust the value of Working Directory as you like. The executable will be stored in selected_working_directory/model_name/ . So, in the example, if I select example as the working directory, the executable will be example/NewtonCooling/NewtonCooling
2. go to Tools > Options > General > disable 'Delete entire simulation directory of the model when OMEdit is closed'

To generate these files, once the settings mentioned above are correctly configured, you should simulate the model in OMEdit at least once.

To re-simulate your model, but not from OMEdit, you can simply use the terminal command to execute the executable which is in a subdirectory of the working directory as described above. This subdirectory is named 'package_name.model_name' and the executable in this subdirectory is executed with the command './model_name' (in Linux). For the example, it will be './NewtonCooling'. Note that you should change the current directory of the terminal session to the example/NewtonCooling/ directory before doing so. This cannot be done from within the example/ directory.

The simulation results will be stored in the same directory as the executable. By default, it is a MAT-file, and you should try to use the same because it is the most efficient (compared to CSV, etc.). At the end of each simulation you should read the MAT-file, and record the variables that you are interested in before performing the next simulation, because this file will be overwritten during re-simulation. You can read the MAT-files generated by OpenModelica using the read_mat_file() function in the example code!

To re-simulate the model, but with different parameters, you can use command line as well. You simply have to structure your command using the -override flag, as the following

'./model_name -override parameter1=value, parameter2=value ..'.

This is also described and demonstrated in the example code.

To implement a loop over the parameter values, you will have to write a Python script. You can execute shell commands using the os package as demonstrated in the example. Take special note that os.changeDir function is used to change the working directory of the virtual shell to the results directory.

This assignment consists of five main parts:

1. Actuator Model Creation (textual)
2. Actuator Model Calibration
3. Plant Model Creation (graphical)
4. Controller Model Creation
5. PID Controller Tuning

Each part is aligned with one of the learning goals

Vehicle suspension systems are critical components that ensure passenger comfort and vehicle stability by isolating the cabin from road irregularities. A well-designed suspension system must balance multiple objectives: maintaining ride height, minimizing vibrations transmitted to passengers, and ensuring that accelerations stay within comfort limits defined by automotive standards.

Figure 5: Simplified illustration of a car suspension system, taken from this paper.

Modern active suspension systems use actuators that can apply controlled forces to adapt the suspension characteristics in real-time. However, these actuators often exhibit nonlinear behavior due to the physics of their spring and damper elements. Accurate modeling of these nonlinearities is essential for designing effective controllers and predicting real-world performance.

In this assignment, you will work with a quarter-car model that includes:

1. a simplified human body model with head, torso, arms, and lower body masses connected by spring-damper elements;
2. an active suspension actuator with both linear and nonlinear (cubic) spring and damper terms that must be calibrated from measured data;
3. and a PID controller that maintains the car body at a desired height (0.5 m) while the vehicle traverses a road profile.

The goal is to design and tune a controller for the active suspension system such that:

1. the car body height closely tracks the 0.5 m setpoint (minimizing the cumulative absolute error);
2. the peak head acceleration remains below 0.315 m/s² (the ISO 2631-1 "not uncomfortable" threshold for seated passengers);
3. and the car height stays within a reasonable range (0.45 m to 0.55 m) throughout the journey.

You are the control engineer and the designer/scientist/client has given you the information below. Your task is to model the whole system using the information given below, in Modelica. Special care should be taken to parametrize the model correctly.

The first step will always be to make an appropriately named Modelica package. This will serve as the package for your entire Modelica assignment. You SHOULD create NEW models ONLY within this package, for each part of the entire assignment.

The actuator sits between the wheel assembly and the car body. It consists of a spring and a damper in parallel, together with an electro-mechanical actuator that can apply an additional force $F_{act}(t)$ on command of the controller. Both the spring and the damper of this actuator behave nonlinearly: the spring stiffens when it is compressed or extended far from its natural length, and the damper resists fast motion more strongly than slow motion. We describe both effects with a cubic term added to the usual linear term.

The manufacturer of the actuator has provided us with the equations that describe its behavior. The actuator is a one-dimensional element with two ends, and its behavior only depends on the relative displacement and relative velocity between these two ends. They are: $$ s_{rel}(t) = s_b(t) - s_a(t) $$ $$ v_{rel}(t) = \frac{d}{dt} s_{rel}(t) $$ $$ f_{spring}(t) = c_1 \cdot \left(s_{rel}(t) - s_{rel,0}\right) + c_2 \cdot \left(s_{rel}(t) - s_{rel,0}\right)^3 $$ $$ f_{damper}(t) = d_1 \cdot v_{rel}(t) + d_2 \cdot \left(v_{rel}(t)\right)^3 $$ $$ f(t) = f_{spring}(t) + f_{damper}(t) + F_{act}(t) $$

Where:

• $s_a(t)$: the position of the lower end of the actuator (attached to the wheel)
• $s_b(t)$: the position of the upper end of the actuator (attached to the car body)
• $s_{rel}(t)$: the length of the actuator, i.e. the relative displacement between its two ends
• $v_{rel}(t)$: the rate of change of the length of the actuator, i.e. the relative velocity between its two ends
• $f_{spring}(t)$: the force due to the spring
• $f_{damper}(t)$: the force due to the damper
• $F_{act}(t)$: the force applied by the actuator, on command of the controller
• $f(t)$: the total force in the actuator
• $s_{rel,0}$: the natural (unloaded) length of the spring
• $c_1$: linear stiffness of the spring
• $c_2$: nonlinear (cubic) stiffness of the spring
• $d_1$: linear damping coefficient of the damper
• $d_2$: nonlinear (cubic) damping coefficient of the damper

Note that the actuator is not a stand-alone system. It must be connected to other physical components (the wheel and the car body) and to a control signal (from the controller). In Modelica, this is done through connectors. Our actuator has three of them:

• flange_a: the lower end of the actuator, of type Modelica.Mechanics.Translational.Interfaces.Flange_a. This is an acausal connector. It carries the position flange_a.s (which corresponds to $s_a$ above) and the force flange_a.f which is applied at that end.
• flange_b: the upper end of the actuator, of type Modelica.Mechanics.Translational.Interfaces.Flange_b. It carries the position flange_b.s (which corresponds to $s_b$ above) and the force flange_b.f.
• control_in: the actuator force command, of type Modelica.Blocks.Interfaces.RealInput. This is a causal connector; its value is decided by whichever block is connected to it (eventually, the controller).

Note that the designers of the mechanical domain part of the Modelica Standard Library have chosen to use position and force in a mechanical interface (known as pin in the electrical domain and flange in the mechanical domain), rather than velocity and force (with product = power), which would be more elegant and consistent with the interfaces in other physical domains. In the mechanical domain, historically, position and force were used.

The force f(t) inside the actuator acts on both of its ends, equally and in opposite directions (Newton's third law). Following the sign convention of the Modelica.Mechanics.Translational library of the MSL, this is expressed with the following two equations: $$ f_a(t) = -f(t) $$ $$ f_b(t) = f(t) $$ where $f_a$ and $f_b$ are the forces of flange_a and flange_b respectively. This means that a positive $f$ (for example, when the spring is stretched beyond its natural length) pulls the two ends of the actuator towards each other.

You are referred to the documentation of Modelica.Mechanics.Translational.Interfaces to learn more about translational connectors. Have a look in particular at PartialCompliant, which is the base class that the MSL uses for its own springs and dampers.

Tasks

You are the engineer and the manufacturer has given you the information above. Your task is to model the actuator as a re-usable Modelica block, using the information given above. Special care should be taken to parametrize the block correctly.

`
1. Create a new Modelica class in your assignment package, but this time of the type 'Block'.
2. First add the three connectors to your block. The easiest way to do this is to drag them from the Libraries pane onto the graphical editor of your block.
   Hint: you could also declare the connectors in the textual editor, but then OMEdit does not know where to draw them, and you will not be able to connect your block graphically in the next parts of the assignment. Dragging them in the graphical editor generates the necessary annotation(Placement(...)) for you.
3. Parametrize the terms $c_1$, $c_2$, $d_1$, $d_2$, and $s_{rel,0}$ in the list above.
4. Model the relative displacement, relative velocity, spring force, damper force, and total force as variables.
5. Implement the given actuator equations in the equation section of your block. Remember that the connector variables can be accessed using the dot-notation. The value of the input connector F_act can be used directly as a variable in your equations.
   Hint: you can use the der() function to represent the time-based derivatives in Modelica, and the ^ operator for powers.
   Hint: use the 'Check Model' button. It reports the number of equations and variables of your block in the message browser, and complains if your block is not balanced (i.e., the number of equations must be equal to the number of variables. If not, the mathematical problem is either over- or under-determined which means it will not have a (unique) solution).
6. Use the following values for the parameters, which are provided by the manufacturer:
   • $c_1 = 20 kN/m$
   • $d_1 = 1.5 kN \cdot s/m$
   • $s_{rel,0} = 35 cm$
   The values of the nonlinear coefficients $c_2$ (in $N/m^3$) and $d_2$ (in $N \cdot s^3/m^3$) are not provided by the manufacturer. You will find them through parameter estimation (model calibration) using measured behaviour traces in the next part of the assignment. For now, set both of them to $0$.
7. Include your code in your report.

Note: A block on its own cannot be simulated, since nothing is connected to its connectors. You will test your block in the next part of the assignment, when you connect it to a calibration experiment.

BEWARE! Often customer requirements will be given in multiple units as in the above (for example metres and centimetres). It is your task to unify and express them all in a single unitary system so that your solution is numerically correct. The recommended way is to use SI units in all cases.

BONUS: for explicitly specifying the units of all the variables and parameters in your Modelica model by reusing definitions from the MSL.

BONUS: for explicitly describing each parameter and variable with comments in your Modelica code.

Note: There are many assumptions in our model, for example, that the actuator is massless, that its force is applied instantaneously, and that it has no limits on its stroke or on the force it can apply.

You have successfully created an abstract model of the actuator, but it is not complete. Remember the nonlinear coefficients $c_2$ and $d_2$ that the manufacturer did not provide? These are parameters that cannot be estimated either by measurement (for example, weighing a component to determine its the mass) or by computation from the design (for example, if weighing is not possible, by using a model of the geometry and the mass density of the component). How will you determine their values so that your model truly represents reality? This is done through parameter estimation, also known as model calibration (though the two are often mixed and described as parameter calibration). The idea is that we try (combinations of) different parameter values, and compare the output of the simulation of the Modelica model with data obtained from measurements on the real system. The parameters, which give the least error (the smallest "distance") between the two traces are said to be the best fit, and are chosen accordingly. Note that this implicitly assumes that the model structure is correct.

An exhaustive way of finding optimal parameter values is to loop over possible values for the parameters to find the combination that gives the closest match mentioned above. This is known as a parameter sweep.

From the manufacturer, we already know that:

• $c_1 = 20 kN/m$
• $d_1 = 1.5 kN \cdot s/m$
• $s_{rel,0} = 35 cm$

Hence, we need to find the nonlinear coefficients $c_2$ and $d_2$. We cannot perform two separate experiments in which the effect of one of the coefficients is nullified: as soon as the actuator moves, both the spring and the damper produce a force. Hence, both coefficients must be estimated together, from a single experiment, by sweeping over all combinations of $c_2$ and $d_2$.

Experimental setup

In the real-world experiment, the actuator was removed from the car and mounted on a test bench:

• The lower end of the actuator (flange_a) is clamped to the (fixed) frame of the test bench.
• The upper end of the actuator (flange_b) is moved by a hydraulic press, which imposes the length of the actuator exactly, regardless of the force that is needed to do so.
• The actuator is switched off during the experiment, i.e. $F_{act} = 0$.
• The press starts at the natural length of the actuator ($35 cm$), compresses it at a constant speed to a length of $5 cm$ in $1$ second, and then extends it at a constant speed back to its natural length in $1$ second. The length of the actuator during the experiment is hence piecewise-linear through the points:
  • $t = 0 s$: $s_{rel} = 35 cm$
  • $t = 1 s$: $s_{rel} = 5 cm$
  • $t = 2 s$: $s_{rel} = 35 cm$
• The force $f$ in the actuator was measured by a load cell, during the whole experiment.
• The above set-up essentially compresses the spring far away from its natural length (where its nonlinear term is significant), while moving it at a constant speed. Notice that during the compression, the force of the damper adds to that of the spring, whereas during the extension, it opposes the spring. The jump in the measured force at $t = 1 s$, when the direction of motion reverses, is the effect of the damper.

The real life-data from the experiment is provided in this csv file. This data represents the measured force $f$ in the actuator (in Newtons), sampled every $2$ milliseconds, for $2$ seconds. As with any real measurement, the data contains some noise.

Figure 6: Imposed length of the actuator (top), and measured force in the actuator (bottom), in the real-life calibration experiment.

Tasks

The main purpose of this part of the assignment is finding the nonlinear coefficients $c_2$ and $d_2$ of the actuator, based on experimental data (i.e., observations/measurements on the real system).

1. Create a new Modelica model in your assignment package, that models the test bench described above. This time, you will use the graphical editor only. Add your actuator block to the model, and connect it to the following blocks from the MSL:
   • Modelica.Mechanics.Translational.Components.Fixed, to model the frame of the test bench.
   • Modelica.Mechanics.Translational.Sources.Position, to model the hydraulic press. Set its parameter exact = true, so that the reference position is followed exactly, without any filtering.
   • Modelica.Blocks.Sources.CombiTimeTable, to provide the reference position to the press (have a look at its documentation here). Use linear interpolation between the points described above.
   • Modelica.Blocks.Sources.Constant, to provide the (zero) actuator force.
2. Parametrize $c_2$ and $d_2$ in this model, and pass them on to the parameters of your actuator block (as described in the introduction to Modelica, at the top of this page), so that they can be changed with the -override flag without re-compiling the model.
3. Simulate the model for $2$ seconds with a time step size of $0.002$ seconds.
4. Create a Python script that executes this Modelica model.
5. This Python script should have a loop that iterates over all possible combinations of the parameters, in the following intervals:
   • $c_2$ in the interval $[0, 400000] N/m^3$, in steps of $10000 N/m^3$
   • $d_2$ in the interval $[0, 15000] N \cdot s^3/m^3$, in steps of $250 N \cdot s^3/m^3$
   That makes $41 \times 61 = 2501$ combinations, hence $2501$ simulations. This is why we compile the model only once!
6. The Python script should be able to record the trace of the force in the actuator, from all of these different simulations. The Python script may store the trace (for each simulation) in memory or save the data as files. You can choose as you prefer.
7. Compute the sum of squared errors (SSE) between each of your simulations and the given trace from real-life data. The simulation that has the smallest total error is the simulation that has the most accurate values for the parameters.
   • You are free to chose how you compute this; whether you use Python, Matlab, Excel...
8. Report the values of the parameters and the SSE in your report, with appropriate graphs comparing the simulation output versus real life data, and (possibly) a heatmap of the SSE as a function of the parameters.
9. Finally, set the default values of $c_2$ and $d_2$ in your actuator block to the values you obtained by calibration. You will use these in the rest of the assignment.

BEWARE! The csv data has a few more data-points than you would expect from a time step of $0.002$ seconds for $2$ seconds. This is because OpenModelica adds extra data-points at events, i.e. at $t = 1 s$ and $t = 2 s$, where the interpolation in the table changes. If your model is set up as described above, your simulation trace will have exactly the same time-stamps as the csv data.

Now that you have a calibrated model of the actuator, you will model the rest of the car and its passenger, i.e. the plant. In contrast to the actuator, you will build the plant using the graphical editor in OMEdit, by connecting blocks from the Modelica.Mechanics.Translational package of the Modelica Standard Library, together with your own actuator block.

TIP: Have a look at the documentation and the examples section in the Modelica.Mechanics.Translational package to learn how to use the different components. This is important to understand how to set the parameters for the different components.

The general assumptions are as follows:

• We only model the vertical motion of one corner of the car (a so-called quarter-car), as shown in Figure 5. All the bodies move along the same vertical axis, and the positive direction is upwards.
• All the bodies are modelled as point masses, i.e. they have no length.
• All the bodies are subject to gravity, with the standard gravitational acceleration $g_n$, which is available in the MSL as Modelica.Constants.g_n.
• Every element that connects two bodies (spring, spring-damper, or actuator) has its flange_a at the lower end, and its flange_b at the upper end.

The Modelica.Mechanics.Translational.Components.Mass block does not include gravitational force on the mass. Since we have six bodies that are all subject to gravity, we provide you with a re-usable block that combines a mass with its weight, in this Modelica file.

Copy the gravity_mass block into your assignment package, for example by pasting its code into the textual view of your package, just before the end of the package. You can access the internal mass of a gravity_mass block using the dot-notation, for example head.mass.s, head.mass.v, and head.mass.a for the position, velocity, and acceleration of a block named head.

The passenger is modelled with a simplified biomechanical model, which divides the human body into four masses (lower body, upper torso, arms, and head), connected to each other by spring-dampers that represent the muscles and joints. The passenger sits on a seat on the car body, which is connected to the wheel by the actuator that you created earlier. The wheel touches the road through the tire.

Figure 7: Lumped parameter quarter-car model (plus occupant) with 6 Degrees of Freedom (6DOF).

The design parameters are set as follows:

• Road: The road surface is modelled by a Modelica.Mechanics.Translational.Sources.Position block, whose reference position is the height of the road. The road profile is smoothed (exact = false) with a critical frequency of 0.8 Hz, which approximates the filtering effect of the contact patch of the tire on the road.
• Tire: The tire is modelled as a linear spring (Modelica.Mechanics.Translational.Components.Spring) with a stiffness of 200 kN/m and a natural length of 30 cm. The damping of the tire is negligible.
• Wheel: The wheel, together with the hub assembly, has a mass of 40 kg.
• Actuator: The calibrated actuator from the previous parts of the assignment connects the wheel to the car body.
• Car body: The mass of the quarter of the car body that is supported by this wheel is 300 kg.
• Seat: The seat is modelled as a spring-damper (Modelica.Mechanics.Translational.Components.SpringDamper) between the car body and the lower body of the passenger, with a stiffness of 40 kN/m, a damping coefficient of 1 kN·s/m, and a natural length of 10 cm.
• Passenger: The masses and connections of the passenger are as follows:
  • Lower body: mass of 24 kg, sitting on the seat.
  • Spine: spring-damper between the lower body and the upper torso, with a stiffness of 70 kN/m, a damping coefficient of 1.5 kN·s/m, and a natural length of 45 cm.
  • Upper torso: mass of 14 kg.
  • Neck: spring-damper between the upper torso and the head, with a stiffness of 34 kN/m, a damping coefficient of 230 N·s/m, and a natural length of 25 cm.
  • Head: mass of 5.1 kg.
  • Shoulders: spring-damper between the upper torso and the arms, with a stiffness of 6.4 kN/m, a damping coefficient of 61 N·s/m, and a natural length of 45 cm.
  • Arms: mass of 7.2 kg.

Implement the system described above, as a Modelica model, using the graphical editor in OMEdit, the provided gravity_mass block for each of the six bodies, and your actuator block. Then do the following:

1. Notice that the description above does not specify the initial positions of the bodies. We assume that the car starts at rest, in static equilibrium. In the textual editor, add an initial equation section to your model, in which you set the velocity and the acceleration of each of the six masses to zero. Modelica will then compute the initial positions of all the bodies for you, such that all the forces are in balance. You don't need to calculate the initial positions yourself!
2. For this part of the assignment, like the previous part, the actuator is switched off. Connect a Modelica.Blocks.Sources.Constant block with value 0 to its F_act input. You will replace it with the PID controller in the next part of the assignment.
3. To the s_ref input of the road, connect a Modelica.Blocks.Sources.Constant block with value 0 (a flat road). Simulate the model for 40 seconds with a step-size of 0.002 seconds.
   Does anything in your model move? It should not, since the car starts in equilibrium and nothing disturbs it. If it does, you have probably made a mistake in the parameters, the signs, or the orientations of your components. Do the positions of the masses relative to each other make sense? For example, the car body should be higher than the wheel. These are excellent self-tests for your implementation! The car body should rest at a height of just under 0.5 m.
   Note: You do not need to provide a plot for this simulation in your report. Do feel free to report IF this helped diagnose a problem in your model.
4. Now replace the constant road with a Modelica.Blocks.Sources.CombiTimeTable block (have a look at its documentation here), that describes the road profile below. The road height is linearly interpolated between the given points, and held constant after the last point. The profile describes a flat road, a long speed bump, two small bumps in quick succession, a dip, and a flat road again:

   Time (s)	0	5	9	12	16	18	19.5	20	21.5	22	24	28	30	35	40
   Road height (cm)	0	0	8	8	0	0	5	0	5	0	0	-5	-5	0	0


Figure 8: Trace of the road height, as described by the above table.

5. Simulate the model again, for 40 seconds with a step-size of 0.002 seconds.

Report your model, and graph the traces of the road height, the height of the wheel, the height of the car body, and the acceleration of the head of the passenger.

Recall the requirements from the problem statement. Does the car body stay within 0.45 m and 0.55 m? Does the acceleration of the head stay below 0.315 m/s²? What do you think is the solution (that we will eventually implement)?

BONUS: for giving meaningful names to all blocks and elements in your model/s.

NOTE: This is a highly simplified model of a car and its passenger. In practice, the car body also pitches and rolls, all four wheels interact with each other, the tire can lose contact with the road, and the human body is modelled with many more degrees of freedom.

You will now create a PID controller for the car suspension system. You will do this in the graphical editor.

A PID (Proportional-Integral-Derivative) controller is widely used in practice to drive the output of the plant towards a set-point. This set-point could change in time. We will however, consider the set-point is constant. In our example, the set-point represents the desired height of the car body. The control output from the PID controller, i.e. $u(t)$ drives the plant / process. The plant produces an output i.e. $h(t)$. The difference between the measured output $h(t)$ and the desired output $r(t)$ (a.k.a the set-point) is called the error $e(t)$. This error value is processed in the PID controller in three different ways as shown in the figure:

Figure 9: PID control loop (adapted from Wikipedia).

• Proportional. The error is multiplied by a constant gain $K_p$.
• Integral. The integral of the error (over time) is multiplied by a constant gain $K_i$.
• Derivative. The derivative of the error (over time) is multiplied by a constant gain $K_d$

Hence, the output from the PID controller $u(t)$ is given by:

$$u(t) = (K_p \cdot e(t)) + (K_i \cdot \int_0^t e(t) dt) + (K_d \cdot \dfrac{d e(t)}{dt}) $$

The gains $K_p$, $K_i$, and $K_d$ modulate the behavior of the PID controller. You may read more about a PID controller online. The essential point is that three signals are summed, which forms the PID control output $u(t)$.

So, to sum up, in our case:

• The plant is the complete quarter-car model with its passenger.
• The measured output from the plant is the height of the car body.
• The set-point is the desired height of the car body, i.e. 0.5 m.
• The control output is the force command $F_{act}$ to the actuator.
• The error signal is the difference between the set-point and the current height of the car body.

Hence, the equation for implementing the control is:

$$ e(t) = r - h(t) $$ coupled with the equations above.

Note: In our simple case, the set-point is constant

Task 1: PID controller block

1. Create another (and I promise, the final!) block in your Modelica assignment package.
2. Even though you have already been given the PID control equations in textual format, you will use the graphical modelling capabilities of Modelica to create the PID controller block.
3. You should first add the connectors to your block, by dragging them onto the graphical editor, the same way you did for the actuator: two connectors of the type Modelica.Blocks.Interfaces.RealInput (one for the set-point $r$, and one for the measured height $h$), and one connector of the type Modelica.Blocks.Interfaces.RealOutput (for the control output $u$).
4. Then, you can use the graphical editor, with blocks from the Modelica.Blocks.Math and Modelica.Blocks.Continuous libraries to create your PID controller, while graphically connecting blocks to each other.
   Hint: Read the description in the introduction to Modelica at the top of this web-page to understand how to parametrize quantities for blocks. You should parametrize the PID control parameters.
   Hint: use a Modelica.Blocks.Continuous.Integrator block with an initial output of 0 for the integral part, and a Modelica.Blocks.Continuous.Derivative block with its default time constant and a steady-state initialization for the derivative part. The latter is an approximate derivative, which avoids differentiating the measured signal exactly.

Note: Do not model any part of the PID controller using textual equations. In this task of the assignment you should use the graphical view to model the PID equations. Only use the textual syntax to define the parameters and interfaces of your block.

Task 2: PID control loop

You can finally use all the blocks that you have defined before, in your own PID control loop!

1. In your plant model, remove the constant input of 0 to the F_act input of the actuator (remember that you had set it earlier, to test the plant model). Retain the CombiTimeTable block for the road profile.
2. Modelica allows re-use of models by means of the 'extend' clause. Create a new Modelica model (lets call it the "controlled car") that extends your plant model.
3. To the controlled car model, add:
   • a Modelica.Mechanics.Translational.Sensors.PositionSensor, connected to the car body, to measure the height $h(t)$ of the car body;
   • a Modelica.Blocks.Sources.Constant block with value 0.5 (in metres), as the set-point $r$;
   • your PID controller block.
4. Connect the set-point and the measured height to the inputs of your PID controller block, and connect the output of the PID controller block to the F_act input of the actuator, as shown in Figure 9.
5. BEWARE of the sign of the control output! Recall the sign convention of the actuator that you modelled: a positive force $f$ in the actuator pulls its two ends towards each other. This means that a positive $F_{act}$ pulls the car body down towards the wheel, and a negative $F_{act}$ pushes the car body up. When the car body is too low, the error $e(t) = r - h(t)$ is positive, and so is the output $u(t)$ of a PID controller with positive gains. If you connect $u(t)$ directly to $F_{act}$, the controller will pull the car body down even further! Hence, the force applied by the actuator must be the negative of the PID output, i.e. $F_{act}(t) = -u(t)$, so that a car body that is too low is raised, and a car body that is too high is lowered. You can implement this either with a Modelica.Blocks.Math.Gain block with gain $-1$ between the controller and the actuator, or by swapping the two inputs of the error computation in your PID block (i.e. computing $h(t) - r$ instead of $r - h(t)$) as it is shown in the figure. Clearly state which option you chose in your report.
   Hint: if the height of the car body runs away from the set-point instead of towards it, you have probably made the sign mistake explained above.
6. Parametrize $K_p$, $K_i$, and $K_d$ in the controlled car model, and pass them on to the parameters of your PID controller block, so that they can be changed with the -override flag without re-compiling the model.

Simulate the controlled car with $K_p = 10000$, $K_i = 1$, and $K_d = 100$ for 40 seconds with a step-size of 0.002 seconds. Report the traces of the height of the car body, the acceleration of the head of the passenger, and the force applied by the actuator. Compare them with the traces of the uncontrolled plant model from the previous section. Do you think this is an optimal controller? Why or why not?

Vary the PID controller parameters and observe the effect on the control outcome. Report your observations.

Note: Do not model any part of the PID control loop using textual equations. Only use the textual syntax to define the parameters of your model.

Similar to the calibration of the actuator parameters (recall that we calibrated $c_2$ and $d_2$), the controller also needs to be tuned for the most ideal control outcome. However, we do not directly have reference data to compare the simulation output to like we had in parameter estimation. Instead, the designer (you) should come up with a cost function to ascertain if the plant behaves as required.

Recall the three requirements from the problem statement:

1. the height of the car body should track the set-point of 0.5 m as closely as possible;
2. the peak acceleration of the head of the passenger should remain below 0.315 m/s²;
3. the height of the car body should stay between 0.45 m and 0.55 m.

The first requirement is an objective: we want to minimize the error between the height of the car body and the set-point. The other two are constraints: a controller that violates them is unacceptable, no matter how well it tracks the set-point.

Task 1: Design of the cost function

1. The cost of a simulation is the cumulative absolute error between the height of the car body and the set-point, over the entire simulation: $$ J = \int_{0}^{t_{end}} | h(t) - r | \, dt $$ Where,
   • $h(t)$ is the height of the car body at time $t$.
   • $r$ is the set-point, i.e. 0.5 m.
   • $t_{end}$ is the duration of the simulation, i.e. 40 seconds.
   The unit of $J$ is metre-seconds (m·s).
2. Unlike in the calibration, you will compute the cost inside your Modelica model, rather than in Python. In the textual editor of your controlled car model, add a variable for the cost, with an initial value of 0, and the equation as described above to compute its value.
   The nice thing about doing this in Modelica is that the value of cost at the end of the simulation is then exactly the cost $J$, computed by the solver, independently of the step-size of the output trace.
3. Report your implementation of this cost function.

Task 2: Early termination of simulations

The tuning will require hundreds of simulations. Many of them are hopeless long before they end:

1. the passenger's head has already been shaken too hard,
2. the car body has already left the allowed band,
3. or the cost has already grown larger than that of the best controller found so far. (Since the cost can only increase over time, such a simulation can never become the best one, so there is no point in simulating it until the end.)

Recall the modelling of discontinuous behavior in Modelica using the 'when' statement, that you learnt in the theory lectures. Here is a reference to the syntax. Inside a when statement, you can call the built-in function terminate("<message>"), which ends the simulation successfully at that moment and prints the message to console.

1. In your controlled car model, parametrize the three limits: the maximum head acceleration (0.315 m/s²), the minimum and maximum heights of the car body (0.45 m and 0.55 m), and a cost limit (with a very large default value, for example $10^{30}$).
2. Add one Boolean variable for each of the three conditions, which indicates whether the corresponding limit was exceeded. For example:
   Boolean headLimitExceeded(start = false, fixed = true);
3. For each of the three conditions, add a when statement that sets the corresponding Boolean variable to true, and terminates the simulation. For example:
   when abs(head.mass.a) > headAccelerationLimit then
     headLimitExceeded = true;
     terminate("Head acceleration exceeded the limit.");
   end when;
   The three conditions are: the absolute acceleration of the head exceeds its limit, the height of the car body leaves the allowed band, and the cost exceeds the cost limit.
4. Simulate the controlled car with the gains from the previous section ($K_p = 10000$, $K_i = 1$, and $K_d = 100$). Does the simulation end early? Which condition caused it, and at what time?
5. Report your implementation of the early termination conditions.

Note: The Boolean variables are stored in the simulation results like any other variable (as 0 or 1), so your Python script can read their final values to find out if/why a simulation ended early.

Task 3: Tune PID controller

1. Create a Python script to simulate the final controlled car model (with the above early termination conditions) with all combinations of the following possible values for the control parameters:
   • $100000 \leq K_p \leq 1000000$ (only multiples of 100000)
   • $K_i \in \{0, 1, 10, 100, 1000\}$
   • $K_d \in \{0, 500, 1000, 2000, 5000, 10000, 20000, 50000\}$
   That makes $10 \times 5 \times 8 = 400$ combinations.
2. The simulation duration should be 40 seconds and the timestep 0.002 seconds.
3. For each simulation, pass the lowest cost of all the acceptable simulations so far as the cost limit to the model, using the -override flag. For the first simulation, there is no such cost yet, so keep the default (very large) cost limit.
4. After each simulation, read the final values of the cost and of the three Boolean variables:
   • If the head acceleration or height limit was exceeded, the controller is unacceptable, and it is discarded.
   • If the cost limit was exceeded, the controller is acceptable, but it is worse than the best controller so far.
   • Otherwise, the simulation ran until the end, and its cost is the new lowest cost.
5. Report the best parameters, its cost, and the traces of the height of the car body, the acceleration of the head of the passenger, and the force applied by the actuator in the best simulation. Compare them with the traces of the uncontrolled car and of the controller from the previous section.
6. Report how many simulations were terminated early because of each of the three conditions.
7. Discuss the results. Does the simulation trace make sense?

BEWARE! The best parameters do not depend on the order in which you simulate the combinations. However, the number of simulations that are terminated early because of the cost limit, and hence the total duration of the tuning, does depend on this order.

BONUS: for measuring and reporting how much time the early termination saves, by also running the tuning with the early termination disabled (for example, by overriding the limits with very large values).

Congratulations! You have successfully modeled a physical system and an optimal controller for it in Modelica!

Note: This is a very naive controller, and the gains that you obtain are only optimal for this particular road profile and passenger. In practice, active suspension controllers are tuned for many different roads and loads, and they must also respect the limits on the force and stroke of the actuator.

• While OpenModelica allows the creation of icons for your models, this is not required for the assignment.
• Student groups with access to only Macs may find that the installation of OpenModelica is not straightforward. Try the Docker image instead. If it doesn't work, please contact the teaching assistant ASAP.