7  Tutorial - Simulation of Maneuvering Motions in Python

In this tutorial we will learn to simulate the motions of a ship in python using the nonlinear equations of motion derived in the class.

After completing this tutorial, the student will be able to:

7.1 Getting Started

We will be using the same repository as the previous tutorial. Once the previous tutorial deadline is over, your instructor will update the upstream repository with the necessary files for this tutorial and raise a pull request on your team repository. You will need to pull in the changes from the upstream repository to your team repository to start this tutorial. The instructions to pull in the changes are given in Appendix A.

After you complete the steps shown in Appendix A to synchronize your repository with the upstream repository, navigate to the root directory of your repository. Pull the latest changes made to your team repository onto your local machine.

cd <your-repository-name>
git pull origin main

The relevant folder structure of the repository for this tutorial is shown below.

<your-repository-name>
├── ros2_ws/
│   ├── src/
│   │   ├── student/
│   │   │   ├── tut_03/
│   │   │   │   ├── hyd.yml
│   │   │   │   ├── input.yml
│   │   │   │   ├── simulate.py
│   │   │   │   ├── read_input.py
│   │   │   │   ├── class_vessel.py
│   │   │   │   ├── module_control.py
│   │   │   │   ├── module_kinematics.py
│   │   │   │   ├── module_dynamics.py
│   │   │   │   └── plot_results.py
│   │   │   └── ...
├── README.md
├── ros2_run_devdocker.sh
├── ros2_run_devdocker.bat
└── LICENSE

Start your own branch from the main branch.

git checkout -b <your-name>

All your group members can work independently on their own branches. However, you will have to merge your branches back to the main branch before the end of the tutorial. We learnt how to do this in Chapter 6 using pull requests.

7.2 Description of Files

The tutorial contains several Python files that work together to simulate vessel motion:

7.2.1 Configuration Files

  • input.yml: Contains vessel parameters including:
    • Physical properties (mass, dimensions, density)
    • Initial conditions for position and velocity
    • Simulation parameters (time step, duration)
    • Control settings
    • Path to hydrodynamic coefficients file
  • hyd.yml: Contains hydrodynamic derivatives including:
    • Surge coefficients (\(X_u\), \(X_{\dot{u}}\))
    • Sway coefficients (\(Y_v\), \(Y_{\dot{v}}\), \(Y_r\), \(Y_{\dot{r}}\))
    • Yaw coefficients (\(N_v\), \(N_{\dot{v}}\), \(N_r\), \(N_{\dot{r}}\))
    • Rudder coefficients (\(Y_{\delta}\), \(N_{\delta}\))

7.2.2 Python Files

  • simulate.py: Main simulation script that:
    • Reads input parameters
    • Creates vessel instance
    • Runs simulation
    • Generates plots
  • read_input.py: Handles reading and parsing of configuration files:
    • Reads YAML files
    • Converts parameters to appropriate data types
    • Returns structured dictionaries of vessel and hydrodynamic parameters
  • class_vessel.py: Contains the Vessel class implementation:
    • Manages vessel state
    • Implements vessel dynamics (vessel_ode)
    • Handles simulation stepping and history tracking
    • Calculates hydrodynamic and control forces
  • module_control.py: Implements control strategies:
    • Fixed rudder control
    • Switching rudder control
  • module_kinematics.py: Handles vessel kinematics:
    • Rotation matrices
    • Euler rate matrix calculation
    • Expressing cross product in matrix form
    • Quaternion operations
    • Conversions between NED, LLH, ECEF frames
  • module_dynamics.py: Handles vessel dynamics:
    • Coriolis matrix calculation
  • plot_results.py: Generates visualization of simulation results:
    • Trajectory plot
    • Velocity time history
    • Turn rate history
    • Rudder angle history
    • Heading angle history

The simulation uses numerical integration (RK45 method) to solve the vessel’s equations of motion, taking into account hydrodynamic forces, control inputs, and vessel kinematics. The results are stored in the vessel’s history attribute and plotted using matplotlib.

7.3 Steps to Execute the Code

7.3.1 Start the Course Docker Container

First start docker desktop and then follow the instructions below to run the course docker container. Make sure to navigate in the terminal to the root folder of your repository.

For Windows:

./ros2_run_devdocker.sh

For Ubuntu (or MacOS):

./ros2_run_devdocker.bat

You would have noticed that when entering the docker container your path is /workspaces/mavlab. The root directory of the repository is mounted to the folder /workspaces/mavlab inside the docker container. The files you edit on your system inside the repository folder are immediately accessible inside the docker container.

7.3.2 Run the Python Script to Simulate the Motion

You can now execute the following command to run the Python code that is present in your repository.

python3 ros2_ws/src/student/tut_03/simulate.py

Now we really see the power of running things inside a docker container. If we were working in our own operating systems, the path of the files will be different for different students. For running the python file, each person would have to edit the above command to make sure python is pointed towards the correct path of the file simulate.py. However, as we are running things inside a docker container, the path of the file simulate.py will not change.

Similarly, observe that you did not have to install any python packages like numpy, scipy or matplotlib. These were pre-installed for you in the docker container. Now everyone in the class is using the exact same version of Python and the same version of packages (numpy, scipy and matplotlib).

Running the above command will result in a plot named simulation_results.png being generated in the root directory of your repository. This plot contains the following subplots:

  1. The trajectory of the ship in the x-y plane
  2. The time history of the velocity of the ship
  3. The time history of the turn rate of the ship
  4. The time history of the rudder angle
  5. The time history of the yaw angle

7.4 Implementation Tasks

There are a total of 25 tasks in this tutorial, organized into three main modules:

Kinematics Module Tasks (module_kinematics.py)

Task 1: Skew-Symmetric Matrix (Smat())

Location: module_kinematics.py, Smat() function

Description: Implement the Smat() function to calculate the cross product of a vector in matrix form.

Input:

  • vec: A 3x1 numpy array representing the vector to be cross-multiplied.

Output:

  • A 3x3 numpy matrix representing the cross product matrix such that when this matrix is multiplied with another vector, it gives the cross product of the input vector with the other vector.

Task 2: Rotation Matrix (eul_to_rotm())

Location: module_kinematics.py, eul_to_rotm() function

Description: Convert Euler angles to rotation matrix.

Input:

  • eul: A 3x1 numpy array representing the Euler angles (roll, pitch, yaw).

Output:

  • A 3x3 numpy matrix representing the rotation matrix such that when this matrix is multiplied with a vector in the body coordinate system, it gives the corresponding vector in the global coordinate system.

Task 3: Euler from Rotation Matrix (rotm_to_eul())

Location: module_kinematics.py, rotm_to_eul() function

Description: Convert rotation matrix back to Euler angles.

Input:

  • rotm: A 3x3 numpy matrix representing the rotation matrix.

Output:

  • A 3x1 numpy array representing the Euler angles (roll, pitch, yaw).

Task 4: Euler Rate Matrix (eul_rate_matrix())

Location: module_kinematics.py, eul_rate_matrix() function

Description: Calculate the Euler rate matrix J2.

Input:

  • eul: A 3x1 numpy array representing the Euler angles (roll, pitch, yaw).

Output:

  • A 3x3 numpy matrix representing the Euler rate matrix J2.

Task 5: Euler rate (eul_rate())

Location: module_kinematics.py, eul_rate() function

Description: Calculate Euler angle rates from angular velocity.

Input:

  • eul: A 3x1 numpy array representing the Euler angles (roll, pitch, yaw).
  • ang: A 3x1 numpy array representing the angular velocity (roll rate, pitch rate, yaw rate).

Output:

  • A 3x1 numpy array representing the Euler angle rates (roll rate, pitch rate, yaw rate).

Task 6: Quaternion from Euler (eul_to_quat())

Location: module_kinematics.py, eul_to_quat() function

Description: Convert Euler angles to quaternion representation.

Input:

  • eul: A 3x1 numpy array representing the Euler angles (roll, pitch, yaw).

Output:

  • A 4x1 numpy array representing the quaternion (w, x, y, z).

Task 7: Euler from Quaternion (quat_to_eul())

Location: module_kinematics.py, quat_to_eul() function

Description: Convert quaternion back to Euler angles.

Input:

  • quat: A 4x1 numpy array representing the quaternion (w, x, y, z).

Output:

  • A 3x1 numpy array representing the Euler angles (roll, pitch, yaw).

Task 8: Rotation Matrix from Quaternion (quat_to_rotm())

Location: module_kinematics.py, quat_to_rotm() function

Description: Convert quaternion to rotation matrix.

Input:

  • quat: A 4x1 numpy array representing the quaternion (w, x, y, z).

Output:

  • A 3x3 numpy matrix representing the rotation matrix such that when this matrix is multiplied with a vector in the body coordinate system, it gives the corresponding vector in the global coordinate system.

Task 9: Quaternion Multiplication (quat_multiply())

Location: module_kinematics.py, quat_multiply() function

Description: Implement quaternion multiplication.

Input:

  • quat1: A 4x1 numpy array representing the first quaternion (w, x, y, z).
  • quat2: A 4x1 numpy array representing the second quaternion (w, x, y, z).

Output:

  • A 4x1 numpy array representing the product of the two quaternions.

Task 10: Quaternion Conjugate (quat_conjugate())

Location: module_kinematics.py, quat_conjugate() function

Description: Calculate quaternion conjugate.

Input:

  • quat: A 4x1 numpy array representing the quaternion (w, x, y, z).

Output:

  • A 4x1 numpy array representing the conjugate of the quaternion.

Task 11: Vector Rotation by Quaternion (rotate_vec_by_quat())

Location: module_kinematics.py, rotate_vec_by_quat() function

Description: Rotate a vector using quaternion.

Input:

  • vec: A 3x1 numpy array representing the vector to be rotated.
  • quat: A 4x1 numpy array representing the quaternion (w, x, y, z).

Output:

  • A 3x1 numpy array representing the rotated vector.

Task 12: Quaternion Rate Matrix (quat_rate_matrix())

Location: module_kinematics.py, quat_rate_matrix() function

Description: Calculate the quaternion rate matrix.

Input:

  • quat: A 4x1 numpy array representing the quaternion (w, x, y, z).

Output:

  • A 4x4 numpy matrix representing the quaternion rate matrix.

Task 13: Quaternion Rate (quat_rate())

Location: module_kinematics.py, quat_rate() function

Description: Calculate the quaternion rate from the angular velocity.

Input:

  • ang: A 3x1 numpy array representing the angular velocity (roll rate, pitch rate, yaw rate).

Output:

  • A 4x1 numpy array representing the quaternion rate (w rate, x rate, y rate, z rate).

Task 14: Smallest Signed Angle (ssa())

Location: module_kinematics.py, ssa() function

Description: Wrap angles to [-180, 180] or [-π, π].

Input:

  • ang: A scalar representing the angle to be wrapped.

Output:

  • A scalar representing the wrapped angle.

Task 15: NED to LLH (ned_to_llh())

Location: module_kinematics.py, ned_to_llh() function

Description: Convert NED coordinates to Latitude-Longitude-Height.

Input:

  • ned: A 3x1 numpy array representing the NED coordinates (north, east, down).

Output:

  • A 3x1 numpy array representing the LLH coordinates (latitude, longitude, height).

Task 16: LLH to NED (llh_to_ned())

Location: module_kinematics.py, llh_to_ned() function

Description: Convert LLH coordinates to NED frame.

Input:

  • llh: A 3x1 numpy array representing the LLH coordinates (latitude, longitude, height).

Output:

  • A 3x1 numpy array representing the NED coordinates (north, east, down).

Task 17: NED to ECEF Rotation (rotm_ned_to_ecef())

Location: module_kinematics.py, rotm_ned_to_ecef() function

Description: Calculate rotation matrix from NED to ECEF.

Input:

  • ned: A 3x1 numpy array representing the NED coordinates (north, east, down).

Output:

  • A 3x3 numpy matrix representing the rotation matrix from NED to ECEF.

Task 18: ECEF to LLH (ecef_to_llh())

Location: module_kinematics.py, ecef_to_llh() function

Description: Convert ECEF coordinates to LLH.

Input:

  • ecef: A 3x1 numpy array representing the ECEF coordinates (x, y, z).

Output:

  • A 3x1 numpy array representing the LLH coordinates (latitude, longitude, height).

Task 19: LLH to ECEF (llh_to_ecef())

Location: module_kinematics.py, llh_to_ecef() function

Description: Convert LLH coordinates to ECEF.

Input:

  • llh: A 3x1 numpy array representing the LLH coordinates (latitude, longitude, height).

Output:

  • A 3x1 numpy array representing the ECEF coordinates (x, y, z).

Control Module Tasks (module_control.py)

Task 20: Fixed Rudder Control (fixed_rudder())

Location: module_control.py, fixed_rudder() function

Description: Implement fixed angle rudder control.

Input:

  • rudder_angle: A scalar representing the rudder angle.

Output:

  • A scalar representing the rudder angle.

Task 21: Switching Rudder Control (switching_rudder())

Location: module_control.py, switching_rudder() function

Description: Implement alternating rudder control.

Input:

  • rudder_angle: A scalar representing the rudder angle.

Output:

  • A scalar representing the rudder angle.

Dynamics Module Tasks (module_dynamics.py)

Task 22: Coriolis Matrix (coriolis_matrix())

Location: module_dynamics.py, coriolis_matrix() function

Description: Calculate the Coriolis matrix using the mass matrix and velocity.

Input:

  • M: A 6x6 numpy matrix representing the mass matrix.
  • vel: A 6x1 numpy array representing the velocity of the vessel.

Output:

  • A 6x6 numpy matrix representing the Coriolis matrix.

Vessel Class Tasks (class_vessel.py)

Task 23: Mass Matrix Generation (_generate_mass_matrix())

Location: class_vessel.py, _generate_mass_matrix() function

Description: Generate the mass matrix.

Input:

  • None

Output:

  • None

Task 24: Dimensionalize Coefficients (_dimensionalize_coefficients())

Location: class_vessel.py, _dimensionalize_coefficients() function

Description: Convert non-dimensional to dimensional coefficients. This function is already implemented for you and hence this task is now redundant.

Input:

  • rho: A scalar representing the density of the fluid.
  • L: A scalar representing the length of the vessel.
  • U: A scalar representing the velocity of the vessel.

Output:

  • None

Task 25: Vessel ODE (vessel_ode())

Location: class_vessel.py, vessel_ode() function

Description: Implement vessel equations of motion.

Input:

  • t: A scalar representing the time.
  • state: A 13x1 numpy array representing the state of the vessel.

Output:

  • A 13x1 numpy array representing the derivative of the state.

Simulation Task (simulate.py)

Task 26: Simulation (simulate())

Location: simulate.py, simulate() function

Description: Implement main simulation loop. This function is already implemented for you. But this will check if your vessel class and simulation code is implemented correctly.

Input:

  • None

Output:

  • None

Each task must pass its corresponding test in the GitHub Actions workflow. The tests verify both the mathematical correctness and proper implementation of each function.

7.5 Testing the Code

To test the code, commit your changes and push them to the remote repository. This will trigger the GitHub Actions workflow to run the tests. You can view the results of the tests by navigating to the “Actions” tab in the GitHub repository. If the tests pass, the workflow will be marked as successful. If the tests fail, you can view the detailed logs to see which tests failed and why. The GitHub workflow will also show the difference between the expected output and the output generated by your code.

7.6 Evaluation

The evaluation of the tutorial will be performed in the classroom. Notice that this is a group activity. You will be working in a group of 4 students. The evaluation will be performed on the group level. In order to pass the tutorial, the group must pass all the tests in the GitHub Actions workflow and also be able to answer all the questions during the evaluation in the classroom. Those who are not present in the classroom on the day of evaluation will not be able to participate in the evaluation and will automatically fail the tutorial.

While one of the objectives of these tutorials is to give you hands on experience with coding and practical implementation, the more important goal is to help you relate what you are doing in your tutorials to the underlying concepts discussed in the class. Therefore, while you are working on the implementation tasks, make sure you are also thinking about the underlying concepts and how they are being applied in the code.

Some questions to think about while you are working on the implementation tasks:

  • Is your code able to pass all the tests in the GitHub Actions workflow?
  • Are you able to show running of the simulation live during the evaluation?
  • What is the difference between fixed and switching rudder control?
  • How have you implemented the switching rudder control?
  • Are we integrating the dynamic equations and kinematics equations separately or together? Where in the code is this happening?
  • How do you think the simulation results would change if we used a different integration method?
  • What is the difference between actual rudder angle and commanded rudder angle?
  • Did you run the simulation and check the results?
  • Do you think the simulation results are in agreement with the expected results? Why or why not?
  • Where is the commanded rudder angle being used in your simulation?

7.7 Instructor Feedback through Pull Requests

Commit your changes to the remote repository often. This will help you to track your progress and also help you to revert to a working version if you make a mistake. In addition, the instructor and TA will be able to provide you with feedback on your code through the Feeback pull request available in your repository. Please do not close this pull request.