Automated Satellite Tracking for Contin- uous Internet Access Design and Implementation of a Dynamic Antenna Control System for Vehicular Internet Connectivity Master’s thesis in Systems, control and mechatronics Armin Khorsandi, Ayman Alzein DEPARTMENT OF ELECTRICAL ENGINEERING CHALMERS UNIVERSITY OF TECHNOLOGY Gothenburg, Sweden 2025 www.chalmers.se www.chalmers.se Master’s thesis 2025 Automated Satellite Tracking for Continuous Internet Access Design and Implementation of a Dynamic Antenna Control System for Vehicular Internet Connectivity Armin Khorsandi Ayman Alzein Department of Electrical Engineering Division of Systems and Control Chalmers University of Technology Gothenburg, Sweden 2025 Automated Satellite Tracking for Continuous Internet Access Design and Implementation of a Dynamic Antenna Control System for Vehicular Internet Connectivity Armin Khorsandi Ayman Alzein © Armin Khorsandi, 2025, arminkh@chalmers.se © Ayman Alzein, 2025, aymana@chalmers.se Industrial Supervisor: Badr Aldeen Alhaffar, Mechatronics SATCUBE Academic Supervisor: Knut Åkesson, Chalmers University of Technology Examiner: Knut Åkesson, Chalmers University of Technology Master’s Thesis 2025 Department of Electrical Engineering Division of Systems and Control Chalmers University of Technology SE-412 96 Gothenburg Telephone +46 31 772 1000 Cover: Control of two-axis gimbal system on moving car Typeset in LATEX, template by Kyriaki Antoniadou-Plytaria Printed by Chalmers Reproservice Gothenburg, Sweden 2025 iv Automated Satellite Tracking for Continuous Internet Access Design and Implementation of a Dynamic Antenna Control System for Vehicular Internet Connectivity Armin Khorsandi Ayman Alzein Department of Electrical Engineering Chalmers University of Technology Abstract This thesis presents the design of a 2-DoF gimbal control system for vehicle-mounted antennas to maintain stable satellite alignment for uninterrupted Internet connec- tivity. Addressing challenges posed by vehicle movement, the work encompasses developing a simulation environment incorporating mathematical model of an an- tenna system, sensor fusion algorithms, and advanced control strategies, including complementary filter, Kalman filter, PID controller, and LQR. The results revealed that while the complementary filter responded more quickly to changes in system dynamics, it suffered from drift, whereas the Kalman filter provided more stable estimations at the cost of a delayed response due to its reliance on model dynam- ics. The study adopted a hybrid control strategy, employing an LQI controller on the outer gimbal angle and a PID controller on the inner gimbal angle. Future work may integrate signal strength sensors, and adopt quaternion-based orientation to overcome gimbal lock. Conducting a robustness analysis of the controllers can further enhance system reliability. This research provides valuable insights into fil- tering and control techniques, offering a robust framework for applications requiring precise orientation estimation and efficient control in dynamic systems. Keywords: Satellite Tracking, Dynamical model, Complementary Filter, Kalman Filter, PID Controller, Linear Quadratic Regulator. v Acknowledgements The authors would like to express their heartfelt gratitude to everyone who pro- vided guidance and support throughout this thesis. We extend our sincere thanks to SATCUBE for giving us the unique opportunity to be the first students to conduct research on their innovative product. Special appreciation goes to our industrial su- pervisor, Badr Aldeen Alhaffar, for his invaluable guidance and support during our time at SATCUBE, and to our academic supervisor and examiner, Knut Åkesson, for his encouragement and expertise throughout the process. We are also deeply thankful to the academic and industrial experts at Chalmers University of Tech- nology and SATCUBE, whose insights and shared experiences enriched this project with both practical and theoretical perspectives. Armin Khorsandi, Ayman Alzein, Gothenburg, January 2025 vii List of Acronyms Below is the list of acronyms that have been used throughout this thesis listed in alphabetical order: Acc Accelerometer ARE Algebraic Riccati Equation CE Characteristic Equation DoF Degree of Freedom ESC Electronic speed controller GNSS Global Navigation Satellite System Gyro Gyroscope LHP Left Half plane LQI Linear Quadratic Integral LQR Linear Quadratic Regulator L Laplace transformation N Normal distribution RHP Right Half plane rot Rotation matrix TJV Target Joints’ Values WGN White Gaussian Noise ix Nomenclature Below is the nomenclature of indices, sets, parameters, and variables that have been used throughout this thesis. Indices ap Action point α Alpha axis β Beta axis e Extended model eff Effective est Estimated value F Friction inner Inner gimbal outer Outer gimbal ref Reference r Reaction s Sample / Satellite V Vehicle x, y, z Indices for specifying the direction Parameters a, b Friction coefficient ∆t Time discretization step (time interval) g Gravitational acceleration γ Tuning coefficient for complementary filter J Inertia xi K Controller gain Kp, Kd, Ki PID tuning parameters Q, R Weighting matrices Variables A, B, C State space matrices α Outer gimbal angle az Global azimuth angle β Inner gimbal angle d Offset vector between center of gravity and a rotation axis e Control error el Global elevation angle G, F Transfer functions H Angular momentum / Open loop transfer function / Measurement model matrix I Integral state / Indicator equation / Identity matrix L Closed loop transfer function λ Local polarization angle m Mass O Rotation matrix of the orientation angles ω Angular velocity p Poles pol Global polarization angle r Reference s Laplace Variable S ARE solution matrix T Time constant / Torque t Time θ Pitch angle φ Roll angle ψ Yaw angle u Input vector x Coordinates / State vector xii y Coordinates / Output vector z Coordinate xiii xiv Contents List of Acronyms ix Nomenclature xi List of Figures xix 1 Introduction 1 1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Problem Description . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.3 Purpose and Goals . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.4 Project Scope . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.5 Research Questions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.6 Limitations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 2 Theory 5 2.1 Dynamical modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.1.1 Equations of Motion . . . . . . . . . . . . . . . . . . . . . . . 7 2.1.1.1 Inner gimbal rotation . . . . . . . . . . . . . . . . . . 8 2.1.1.2 Outer gimbal rotation . . . . . . . . . . . . . . . . . 9 2.1.1.3 Summary . . . . . . . . . . . . . . . . . . . . . . . . 10 2.1.2 Friction modeling . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.1.3 ESC Modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.1.4 Sensor modeling . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.1.4.1 Encoders . . . . . . . . . . . . . . . . . . . . . . . . 14 2.1.4.2 Differential GNSS . . . . . . . . . . . . . . . . . . . 14 2.1.4.3 Accelerometer . . . . . . . . . . . . . . . . . . . . . . 15 2.1.4.4 Gyroscope . . . . . . . . . . . . . . . . . . . . . . . . 17 2.2 Kinematics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 2.2.1 Global Target . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 2.2.2 Target joint values . . . . . . . . . . . . . . . . . . . . . . . . 21 2.3 Processing step for sensors . . . . . . . . . . . . . . . . . . . . . . . . 22 2.3.1 Encoder Data Processing . . . . . . . . . . . . . . . . . . . . . 22 2.3.2 Gyroscope Data Processing . . . . . . . . . . . . . . . . . . . 22 2.3.3 Accelerometer Data Processing . . . . . . . . . . . . . . . . . 23 2.4 Sensor Fusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.4.1 Complementary filter . . . . . . . . . . . . . . . . . . . . . . . 24 2.4.1.1 Application in This Research . . . . . . . . . . . . . 26 xv Contents 2.4.2 Kalman filter . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 2.4.2.1 Prediction Step . . . . . . . . . . . . . . . . . . . . . 27 2.4.2.2 Update Step . . . . . . . . . . . . . . . . . . . . . . 27 2.4.2.3 Application in This Research . . . . . . . . . . . . . 28 2.4.2.3.1 Kalman based on IMU and differential GNSS 29 2.4.2.3.2 Kalman Filter Based on IMU-Only . . . . . 29 2.4.2.4 Tuning the Complementary and Kalman Filters . . . 30 2.5 Control Loop . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 2.5.1 PID Controllers . . . . . . . . . . . . . . . . . . . . . . . . . . 34 2.5.1.1 P Controller . . . . . . . . . . . . . . . . . . . . . . . 35 2.5.1.2 D Controller . . . . . . . . . . . . . . . . . . . . . . 35 2.5.1.3 PD Controller . . . . . . . . . . . . . . . . . . . . . . 36 2.5.1.4 PID Controller . . . . . . . . . . . . . . . . . . . . . 37 2.5.2 LQR Controllers . . . . . . . . . . . . . . . . . . . . . . . . . 39 2.5.2.1 LQR model . . . . . . . . . . . . . . . . . . . . . . . 40 2.5.2.2 LQI model . . . . . . . . . . . . . . . . . . . . . . . 40 2.5.2.3 Tuning LQR . . . . . . . . . . . . . . . . . . . . . . 41 3 Methodology 43 3.1 Literature Study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 3.2 Simulation Environment . . . . . . . . . . . . . . . . . . . . . . . . . 44 3.3 System Modeling Using Simscape . . . . . . . . . . . . . . . . . . . . 45 3.4 Trajectory Design for Validation . . . . . . . . . . . . . . . . . . . . . 46 3.5 3D Visualization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 4 Results 49 4.1 Modeling Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 4.1.1 The Linear Equations of Motion . . . . . . . . . . . . . . . . . 49 4.1.2 The Simplified Gimbal Transfer Function . . . . . . . . . . . . 50 4.2 Sensor Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 4.3 Sensor Processing and State Estimation . . . . . . . . . . . . . . . . . 53 4.3.1 Complementary Filter . . . . . . . . . . . . . . . . . . . . . . 54 4.3.2 Kalman Filter Based on IMU and Differential GNSS . . . . . 55 4.3.3 Kalman Filter Based Only on IMU . . . . . . . . . . . . . . . 56 4.4 Controller Design . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 5 Discussion 63 5.1 Sensors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 5.1.1 Magnetometer . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 5.1.2 GNSS and translation error . . . . . . . . . . . . . . . . . . . 64 5.1.2.1 Error in Global Target Elevation Angle . . . . . . . . 64 5.1.2.2 Error in Global Target Azimuth Angle . . . . . . . . 65 5.2 Answering Research Questions . . . . . . . . . . . . . . . . . . . . . . 65 6 Conclusion 67 7 Future Recommendations 69 xvi Contents 7.0.1 Signal Strength . . . . . . . . . . . . . . . . . . . . . . . . . . 69 7.0.2 Euler Angles to Quaternions . . . . . . . . . . . . . . . . . . . 70 7.0.3 Robustness Analysis . . . . . . . . . . . . . . . . . . . . . . . 70 Bibliography 71 xvii Contents xviii List of Figures 2.1 General block diagram of the system, showcasing the interactions between the satellite’s fixed position inputs, vehicle orientation and position, kinematics (Global target and TJV), sensor modeling, fil- tering, control, and dynamics (GGimbal). The 3D block provides a real-time visualization of the system through Simscape Multibody. . . 6 2.2 Terminal’s components and corresponding reference coordinate sys- tems. Antenna’s line of sight is along yβ . . . . . . . . . . . . . . . . 7 2.3 Different friction models: a) Coulomb’s friction; b) coulomb’s and viscous friction; c) Simplified model incorporating coulomb’s, viscous, and static friction; d) A realistic representation of the friction force [1]. 11 2.4 Comparison of the step response of the actual ESC and the first- order model. The ESC is tasked with delivering 90% of the maximum torque specified by the motor manufacturer. . . . . . . . . . . . . . . 13 2.5 Encoder model in simulation. The model takes the true angle (in de- grees) and converts it to encoder counts. Here, Max Counts denotes the total number of pulses generated by the encoder per revolution, as specified in the device datasheet. . . . . . . . . . . . . . . . . . . . 14 2.6 Differential GNSS module model in simulation. It takes the position and the heading action points at each time step and converts them to signals from the GNSS sensor. The position measurements are assumed to be ideal, whereas for the heading estimation, the WGN is characterized based on analysis of logged sensor data. . . . . . . . . 15 2.7 A top view of the terminal shown in Figure 2.2. When the car changes its heading rate, the antenna moves in the opposite direction to main- tain line-of-sight with the satellite. The difference between these an- gular accelerations results in a net angular acceleration acting on the system, introducing an induced acceleration Acci in the x direction of the accelerometer’s body frame. . . . . . . . . . . . . . . . . . . . . 16 2.8 The accelerometer model in simulation. As mentioned, external linear accelerations introduced by the car are not considered. . . . . . . . . 17 2.9 The gyroscope model in simulation. The term "R andom walk" is a standard term for WGN in gyroscope contexts. Both the bias instability and random walk are modeled based on the logged data. . 18 xix List of Figures 2.10 The function Global Target calculates the target elevation, azimuth, and polarization angles in a global reference frame, based on the po- sitions of the terminal and the satellite. Meanwhile, the function Target Joint Values computes the local joint angles β, λ and α re- quired to adjust for the terminal’s orientation, ensuring the antenna is correctly aligned with the satellite. . . . . . . . . . . . . . . . . . . 19 2.11 An illustration depicting how azimuth and elevation angles are cal- culated. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.12 Encoder processing step to estimate the angle α as a block diagram. . 22 2.13 Gyroscope processing step to estimate roll, pitch, and yaw angles as a block diagram. The estimated α is differentiated and subtracted from the gyroscope’s z-axis measurements. The data is then rotated back along the z-axis by α to isolate the raw measurements, which are subsequently integrated to estimate the orientation angles. . . . . 23 2.14 Accelerometer processing step to estimate roll and pitch angles as a block diagram. The accelerometer data is processed by removing the effects of angles α and ψest and rotating the data back along the z- axis. The roll and pitch angles are then calculated using the equations above. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.15 The simulation model representing the dynamics from requested torques to angles. The model includes two uncoupled systems, one for each axis. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 2.16 A simplified version of the block diagram shown in Figure 2.15, fo- cusing on the relationship between a single requested torque and its resulting angle. For simplicity, the loop transfer function depicted in blue is denoted as Gsub(s), while the transfer function for the satura- tion is denoted as F (s). . . . . . . . . . . . . . . . . . . . . . . . . . . 32 2.17 Closed-loop control system for GGimbal using a P-controller with gain kp. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 2.18 Closed-loop control system for GGimbal using a D-controller with gain kd. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 2.19 Closed-loop control system for GGimbal using a PD-controller. . . . . . 37 3.1 Overview of the simulation system implemented in Simulink. This diagram mirrors the conceptual block diagram shown in figure 2.1. . 44 3.2 Trajectory created for validation of the feedback control system, fea- turing varying road conditions, including pitch, roll, speed bumps, roadside inclines, and gravel roads, with 90-degree turns of different radius. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 3.3 The 3D simulation environment created using Simscape Multibody. The vehicle, equipped with the antenna system mounted on its roof is shown tracking the trajectory. The antenna remains aligned with the satellite and the line of site represented by the red line. . . . . . . 48 xx List of Figures 4.1 The graph on the right is a zoomed-in view of the left graph. A simu- lation was performed for both the nonlinear and linear models in the simulation environment. The models were subjected to torque inputs at t = 1 second, with Tα = 0.3 N and Tβ = 0.05 N. For the inner gim- bal, the β angles showed complete overlap between the two models. However, the α angles exhibited a slight deviation, reaching approx- imately 10 degrees after 10 seconds of simulation, as highlighted in the zoomed-in graph. . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 4.2 Heading angle estimations provided by the differential GNSS module recorded over an extended period with the module stationary. . . . . 52 4.3 This figure presents sub-figures A to E, illustrating the sensor signals while driving along the validation trajectory and locking the outer gimbal to ensure that the accelerometer and gyroscope data represent the vehicle frame. Sub-figure A corresponds to the trajectory segment shown in Figure 3.2. Sub-figure B depicts the action points (true values) of the vehicle orientation in the simulation. Sub-figure C shows the estimated heading angle using the differential GNSS sensor model. Sub-figure D illustrates the acceleration signal generated by the accelerometer model. Sub-figure E represents the angular velocity signal generated by the gyroscope model. . . . . . . . . . . . . . . . . 53 4.4 Orientation estimation using a complementary filter with a gain of 0.9. The figure shows the estimation After 10 laps around the vali- dation trajectory. At the end of the trajectory, the orientation action points are [roll, pitch, yaw] = [0, 0, 360] deg, while the estimation is [1.40, 1.44, 362.01] deg indicating drift. . . . . . . . . . . . . . . . . . 54 4.5 Orientation estimation using a Kalman filter based on IMU and dif- ferential GNSS. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 4.6 A comparison of the estimation performance between the complemen- tary filter and the Kalman filter. The figure shows the estimation re- sults after running 10 laps of the validation trajectory at segment 6, which simulates a roll change on a strait road. When the roll action point’s rate changes, the complementary filter responds immediately, while the Kalman filter exhibits a slight delay before adjusting to the action point. However, the complementary filter shows a noticeable drift (offset) from the action point. . . . . . . . . . . . . . . . . . . . 56 4.7 Orientation estimation using a Kalman filter based only on IMU. . . 56 4.8 A PD controller is implemented using the parameter provided in equa- tion 4.3, and the control error and the control action Treq are recorded. The errors are shown in subfigures B and C with different magnitude scaling, while subfigure C illustrates the control action where Tα and Tβ are the requested torques from the system, and Sα and Sβ are the saturation levels for each motor. The figure is obtained by running the full simulation model with the Kalman filter based on IMU and differential GNSS activated. . . . . . . . . . . . . . . . . . . . . . . . 59 xxi List of Figures 4.9 A PID controller is implemented using the parameters in equation 4.4, and the control error e and control action Treq are recorded. The errors are shown in subfigures B and C with different magnitude scaling, while subfigure C illustrates the control action, where Tα and Tβ are the requested torques from the system, and Sα and Sβ are the saturation levels for each motor. The figure is obtained by running the full simulation model with the Kalman filter based on IMU and differential GNSS activated. . . . . . . . . . . . . . . . . . . . . . . . 60 4.10 An LQR controller is implemented using the weighting matrices in equation 4.5 and the state-space model in equation 2.74. The control error e and control action Treq are recorded. The errors are shown in subfigures B and C with different magnitude scaling, while subfigure C illustrates the control action, where Tα and Tβ are the requested torques from the system, and Sα and Sβ are the saturation levels for each motor. The figure is obtained by running the full simulation model with the Kalman filter based on IMU and differential GNSS activated. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 4.11 An Extended LQR controller is implemented using the weighting ma- trices in equation 4.8 and the state-space model in equation 2.76. The control error e and control action Treq are recorded. The errors are shown in subfigures B and C with different scaling, while subfigure C illustrates the control action, where Tα and Tβ are the requested torques from the system, and Sα and Sβ are the saturation levels for each motor. The figure is obtained by running the full simulation model with the Kalman filter based on IMU and differential GNSS activated. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 5.1 The elevation error introduced by the vehicle’s height above sea level, the height of the vehicle where the terminal is mounted and the offset between the inner and outer gimbals. . . . . . . . . . . . . . . . . . . 64 5.2 Top view: GNSS estimates the vehicle to be at Pest, while it could be anywhere inside the green circle, representing the uncertainty of the circle. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 xxii 1 Introduction This chapter introduces the topic of satellite Internet, with a specific focus on satellite-based Internet on the move, and explores its significance as a vital modern solution for various organizations. It provides a concise overview of the primary problem addressed in this research. Additionally, the goals and limitations of the project are outlined to provide context for the study. 1.1 Background Internet connectivity has become a fundamental part of modern life, integrated into nearly every aspect of daily activities and professional tasks. Its availability is cru- cial in various locations and circumstances. For example, a wealthy individual on a yacht in the middle of the ocean relies on a stable Internet connection for com- munication and entertainment. A journalist reporting from a conflict zone depends on internet access to provide timely updates. Health organizations conducting mis- sions in remote or disaster-stricken areas need reliable connectivity to coordinate with headquarters. Rescue teams require Internet access to carry out lifesaving operations effectively. In all of these scenarios, satellite-based Internet emerges as a viable solution, offering uninterrupted connectivity even in the most challenging environments. 1.2 Problem Description Many satellite operators position their satellites in geostationary orbit, including those considered in this thesis. These satellites are located around the equator and remain fixed relative to the Earth. To establish a stable connection with such satellites, ground terminals equipped with directional antennas are required. This arrangement poses a significant challenge, particularly in maintaining a clear line of sight between the satellite and the ground terminal’s antenna. For stationary ground terminals, proper installation and calibration ensure reliable connectivity. However, when these terminals are mounted on moving vehicles, maintaining a sta- ble connection becomes significantly more complex. The antennas must automati- cally adjust in real-time to account for vehicle movement, continuously aligning to maintain an uninterrupted link with the satellite. 1 1. Introduction This research is conducted in collaboration with SATCUBE, a leading company spe- cializing in satellite-based Internet solutions for seamless connectivity in remote and demanding locations to address these challenges by focusing on the design and im- plementation of a simulation environment and a control system for vehicle-mounted satellite terminals. Specifically, the satellite terminal, equipped with an antenna with two degrees of freedom, is mounted on a moving vehicle. The antenna will be aligned with a known geostationary satellite to maintain continuous Internet connectivity. 1.3 Purpose and Goals The purpose of this thesis is to design and develop an automated satellite tracking system capable of controlling the antenna’s movement to maintain precise, active, and stable alignment with a satellite. 1.4 Project Scope Defining a clear and focused project scope is essential to achieving the purpose of this research. The key aspects of the project scope are outlined below: • Developing a mathematical model of the antenna system to accurately repre- sent its dynamics and behavior. • Implementing advanced sensor fusion and filtering techniques to process sensor data effectively within the simulation environment, ensuring accurate state estimation. • Designing and implementing a control system for the antenna to ensure align- ment with the known satellite. • Conducting rigorous validation of the control system through a series of ex- periments in the simulation environment, replicating real-world scenarios. This defined scope serves as the foundation for the project, guiding the develop- ment process and ensuring that all efforts remains focused on achieving the primary objective. 1.5 Research Questions This research aims to address the following key questions: • Which filtering technique offers the highest precision and reliability for pro- cessing data in this specific context? • What methods can be used to evaluate and compare different control strate- gies, and which control algorithm demonstrates superior effectiveness when implemented in MATLAB/Simulink? 2 1. Introduction 1.6 Limitations This research is conducted within a constrained timeframe, which leads to several exclusions and limitations: • Hardware Implementation: Due to time constraints, the system will not be implemented on a physical prototype and is only going to be developed and evaluated in a simulation environment. • Limited Algorithms and Methods: Only a limited range of sensors, filter- ing techniques, and controllers will be explored and implemented, reflecting the time limitations of the study. Furthermore, the sensitive nature of proprietary information related to the com- pany’s product restricts access to certain data and resources, which further shapes the scope of this research. 3 1. Introduction 4 2 Theory This chapter establishes the theoretical foundation for the research, encompassing a comprehensive analysis of system dynamics, sensor modeling, data processing, and control mechanisms. It begins with the fundamental equations of motion that govern the physical interactions within the system and explores the complexities of friction and torque dynamics. The kinematics section then provides an understanding of the spatial relationships and orientation adjustments necessary for system alignment. Subsequent sections delve into the detailed modeling and processing of sensor data, including encoders, gyroscopes, accelerometers, and Global Navigation Satellite Sys- tem (GNSS). The chapter also introduces sensor fusion techniques, comparing the effectiveness of complementary and Kalman filters in integrating multi-sensor data for optimal state estimation. Finally, the control loop design is presented, with an emphasis on feedback controllers such as PID and LQR, highlighting their roles in achieving system stability and responsiveness. This theoretical groundwork serves as the basis for the simulations and experimental validations discussed in later chap- ters. Figure 2.1 shows the overall system architecture modeled in Simulink. It highlights the interactions between key subsystems, including the satellite’s fixed position in- puts (X, Y, Z) and the vehicle’s time-varying position and orientation inputs (X, Y, Z) along with roll, pitch, and yaw. The "Global target" and "TJV" blocks are dis- cussed in Section 2.2. The "Sensors" and "Estimator" blocks are covered in Sections 2.1.4 and 2.3. The "Controller" block is presented in Section 2.5, and the "GGimbal" block, which incorporates friction, the ESC, and the antenna’s dynamics, is detailed in Section 2.1. 5 2. Theory Figure 2.1: General block diagram of the system, showcasing the interactions between the satellite’s fixed position inputs, vehicle orientation and position, kine- matics (Global target and TJV), sensor modeling, filtering, control, and dynamics (GGimbal). The 3D block provides a real-time visualization of the system through Simscape Multibody. 2.1 Dynamical modeling The mechanical system comprises two main components: the base and the an- tenna, which consists of an outer and an inner gimbal. The global coordinate frame [xg, yg, zg] serves as the reference for all angular measurements. Since the terminal is mounted on a vehicle, the orientation and position of the base are aligned with the vehicle. The outer gimbal has its own local coordinate frame [xα, yα, zα], while the inner gimbal’s frame [xβ, yβ, zβ] is shifted from the outer gimbal along the zα axis. To align the antenna with the satellite, precise calculations are needed to determine the angles α and β. A detailed explanation of these calculations is provided in section 2.2. The terminal is equipped with two motors: the beta motor, which controls the rotation of the inner gimbal around its xβ axis, corresponding to angle β, and the alpha motor, which rotates the outer gimbal around its zα axis, corresponding to angle α. A simplified representation of the terminal is illustrated in Figure 2.2. 6 2. Theory Figure 2.2: Terminal’s components and corresponding reference coordinate sys- tems. Antenna’s line of sight is along yβ 2.1.1 Equations of Motion When torque is applied to a rigid body, it causes the body to rotate about the axis aligned with the torque. The angular acceleration produced is inversely proportional to the body’s moments and products of inertia. The moments of inertia correspond to the diagonal elements of the inertia tensor, while the products of inertia are represented by the off-diagonal elements. If the products of inertia are non-zero, the applied torque can also induce rotational motion about the other two axes, a phenomenon known as the gyroscopic effect [3]. For example, when torque is applied along the x-axis, the products of inertia Jxy and Jxz cause rotation around the y and z axes, respectively. The same principle applies for torque applied along the y and z axes. According to Meriam et al. [3], the equation of motion for a rigid body can be written as follows: Ḣ = ∑ T − ω ×H (2.1) Where H = Jω is the angular momentum. If the body is constrained to rotate around only one axis (the axis where the torque is applied), the constraints are considered as a locking mechanism which applies a reaction torque, Treaction that matches the magnitude but opposes the direction of any torque that would otherwise cause rotation about the other two axes due to the products of inertia. As a result, this reaction torque cancels out any rotational motion about the other axes, ensuring the body rotates solely around the axis of applied torque. Jω̇ = Tapplied − ω × Jω − Treaction (2.2) 7 2. Theory Jxx Jxy Jxz Jxy Jyy Jyz Jxz Jyz Jzz  ω̇x ω̇y ω̇z  = Tx Ty Tz − ωx ωy ωz × Jxx Jxy Jxz Jxy Jyy Jyz Jxz Jyz Jzz  ωx ωy ωz − Trx Try Trz  (2.3) Jxxω̇x + Jxyω̇y + Jxzω̇z = Tx − ωy(Jxzωx + Jyzωy + Jzzωz) + ωz(Jxyωx + Jyyωy + Jyzωz) − Trx Jxyω̇x + Jyyω̇y + Jyzω̇z = Ty − ωz(Jxxωx + Jxyωy + Jxzωz) + ωx(Jxzωx + Jyzωy + Jzzωz) − Try Jxzω̇x + Jyzω̇y + Jzzω̇z = Tz − ωx(Jxyωx + Jyyωy + Jyzωz) + ωy(Jxxωx + Jxyωy + Jxzωz) − Trz  (2.4) In equation 2.4, Tx, Ty and Tz represent the torques applied along each axis. Mean- while, Trx, Try and Trz denote the reaction torques that constrain the body’s rotation in two of the axes. This type of motion constraint, often referred to as rotational locking, can be achieved using ball bearings. However, ball bearings introduce fric- tion torques due to their mechanical resistance. As a result, the applied torque is reduced by the friction torque, and the relationship can be expressed as Tapplied = Tmotor − Tfriction. Thus, the actual torque acting on the system is the motor torque minus the torque lost to friction. A similar interpretation applies to the rotation of each gimbal. The shapes of both the inner and outer gimbals are asymmetrical, so their inertia tensors are non- diagonal. Each gimbal has one degree of freedom (DoF) and rotates using ball bearings. 2.1.1.1 Inner gimbal rotation The entire terminal is already 3D modeled in detail, with the masses correctly positioned. As a result, the inertia tensor, J , can be easily extracted by specifying the correct coordinate system in the 3D modeling software. The inertia that the beta motor needs to overcome is that of the inner gimbal, Jβ, represented as a constant 3 × 3 matrix. Since the applied torque is only provided around the xβ axis, the applied torque vector will be Tapplied = [Tx; 0; 0]. The reaction torques must lock the rotation around the yβ and zβ axes. The angular accelerations around y and z are induced by the applied torque Tapplied on xβ, as it is the only external torque acting on the inner gimbal. This coupling between the angular accelerations in y and z and the applied torque on xβ arises from the off-diagonal elements in the inertia tensor, specifically Jxy and Jxz. To prevent angular acceleration in the y and z directions, the reaction torque Treaction must cancel out all terms in equation 2.4, including the coupling terms involving Jxy and Jxz. This ensures that the coupling between the applied torque, Tx, and the angular accelerations ω̇y and ω̇z is eliminated. Treaction = Trx Try Trz  = −Jxyω̇y − Jxzω̇z − Jxzωxωy + Jxyωxωz −Jxyω̇x − Jxyωyωz − Jxzω 2 x + Jxzω 2 x −Jxzω̇x − Jxyω 2 x + Jxyω 2 y + Jxzωyωz  (2.5) 8 2. Theory Consequently, the rotational equations describing the relationship between the ap- plied torque and the resulting rotation around all axes of the inner gimbal are as follows: Jxxω̇x = Tx − ωy(Jyzωy + Jzzωz) + ωz(Jyyωy + Jyzωz) Jyyω̇y + Jyzω̇z = −ωz(Jxxωx) + ωx(Jyzωy + Jzzωz) Jyzω̇y + Jzzω̇z = −ωx(Jyyωy + Jyzωz) + ωy(Jxxωx)  (2.6) 2.1.1.2 Outer gimbal rotation In the same way as the inner gimbal, the motion in the outer gimbal can also be expressed as shown in equation 2.4. In this case, the applied torque is Tapplied = [0; 0;Tz], as the alpha motor exerts torque along the zα-axis. The reaction torques in the outer gimbal counteract the resulting rotations around the xα and yα axes. These undesired rotations arise from the off-diagonal inertia terms, Jxz and Jyz, which couple the applied torque along zα to rotations around the other axes. Therefore, the reaction torques are as follows: Treaction = Trx Try Trz  =  −Jxzω̇z − Jxzωxωy − Jyzω 2 y + Jyzω 2 z −Jyzω̇z − Jxzω 2 z + Jxzω 2 x + Jyzωyωx −Jxzω̇x − Jyzω̇y − Jyzωxωz + Jxzωyωz  (2.7) Consequently, the rotational equations describing the relationship between the ap- plied torque and the resulting rotation around all axes of the outer gimbal are as follows:  Jxxω̇x + Jxyω̇y = −ωy(Jzzωz) + ωz(Jxyωx + Jyyωy) Jxyω̇x + Jyyω̇y = −ωz(Jxxωx + Jxyωy) + ωx(Jzzωz) Jzzω̇z = Tz − ωx(Jxyωx + Jyyωy) + ωy(Jxxωx + Jxyωy)  (2.8) The alpha motor must counteract the combined inertia of both the outer and inner gimbals, as the inner gimbal is attached to the outer gimbal. Since the inner gimbal’s moment of inertia influences the outer gimbal variably, depending on the beta angle, the J matrix in equation 2.8 should incorporate the outer gimbal’s inertia, Jα, and the transformed inertia of the inner gimbal, Jβ, which varies based on the beta angle β. Jtotal = Jα + Jβα (2.9) Jα and Jβ are obtained from the 3D modeling software, while Jβα is calculated by transforming Jβ into the outer gimbal’s coordinates and then adding the inertial offset. Jβα = rotx(β) Jβ rot −1 x (β) +mβ · (|d|2I3 − dd′) (2.10) 9 2. Theory In equation 2.10, rotx(β) is the standard rotation matrix around the x axis, where the rotation angle is β. The term rotx(β) Jβ rot −1 x (β) is derived from the similarity transformation [5]. The term mβ · (|d|2I3 − dd′) is derived from the generalized parallel axis theorem [4] and represents the effect of the inner gimbal’s center of gravity being offseted from the outer gimbal’s coordinate system. Where mβ is the mass of the inner gimbal, I3 is an identety matrix of size 3 × 3, and d is a 3 × 1 vector that includes the offsets along xα, yα, and zα. 2.1.1.3 Summary Incorporating the disturbances arising from the presence of products of inertia, and then including Treaction to counteract the rotational motion around the other two axes, restricts the rotation to the desired axis. This results in the angular velocities and accelerations of the other axes being zero. By applying this to equations 2.6 and 2.8, the outcome can be expressed as follows: Jβxx ω̇βx = Tx In the inner gimbal Jtotalzz ω̇αz = Tz In the outer gimbal (2.11) By assuming that the mass center of the inner gimbal is located on the rotational axis of the outer gimbal, and that Jβyy = Jβzz , it follows that Jtotalzz = Jαzz + Jβyy . This assumption introduces a slight modification to Equation 2.11, leading to linearized equations describing the motion of the system as follows: Jβxx ω̇βx = Tx For the inner gimbal, (Jαzz + Jβyy) ω̇αz = Tz For the outer gimbal. (2.12) 2.1.2 Friction modeling The outer gimbal is mounted to the base using a ball bearing, and the cables run from the base to the outer gimbal through the bearing using a slip ring. The combination of the slip ring and the bearing introduces friction that must be considered. The inner gimbal is mounted to the outer gimbal using another ball bearing, and there are some cables connecting the inner and outer gimbals. These cables are flexible, and their resistance is assumed to be negligible. The friction force can be described using different models, as shown in figure 2.3. According to Olsson et al, [1], the friction force in linear motion is represented by various models, each addressing certain aspects while neglecting others. The model illustrated in figure 2.3a describes the friction force as a sign function, where the fric- tion remains constant regardless of the velocity; this model only considers coulomb’s friction. The model in figure 2.3b presents a linear proportional relationship between the friction force and velocity, accounting for both coulomb’s friction and viscous 10 2. Theory friction. The model depicted in figure 2.3c incorporates the same relationships as model b but also addresses static friction when the velocity is zero. The most ac- curate representation is shown in figure 2.3d, where the model takes into account static friction and the dependencies between friction force and velocity, indicating that the friction force decreases from the static friction value at low velocities before it starts to increase again. Similarly, in rotational motion, the same models can be applied where the friction force in linear motion will be represented as friction torque in rotational motion. The model illustrated in figure 2.3b is the most effective to use in this research since it captures the relationship between velocity and friction torques, and is relatively easier to model by conducting experiments on the physical terminal. Figure 2.3: Different friction models: a) Coulomb’s friction; b) coulomb’s and viscous friction; c) Simplified model incorporating coulomb’s, viscous, and static friction; d) A realistic representation of the friction force [1]. Assuming the friction torque is linear, it can be described by the equation Tfriction = aω + b (2.13) To obtain the parameters a and b, an experiment can be conducted where the output torque of each motor is recorded while maintaining a constant angular speed. By conducting the experiment twice for each motor at two different constant angular speeds, it is possible to find the equation of the fitting line that passes through the recorded points, resulting in the values of a and b for each gimbal. 11 2. Theory 2.1.3 ESC Modeling The Electronic Speed Controller (ESC) is responsible for generating the desired torque on the motor shaft as demanded by the controller. However, there is a delay between the torque request and the actual torque reaching its desired magnitude. This delay arises due to the internal dynamics of the ESC. For modeling purposes, the ESC is treated as a black box. A step torque input is applied to the physical hardware, and the resulting torque output is recorded. The system dynamics are then identified based on this recorded torque data. The dynamics of the ESC are approximated using a first-order transfer function: 1 Ts+ 1 (2.14) In this expression, T is the time constant, which represents the time required for the motor’s response to reach 63.2% of its final value in reaction to a step input. The value of 63.2% is derived from the step response equation in the time domain. Let y(t) denote the system’s step response. The input to the system is a unit step function with a magnitude of 1, which has a Laplace transform of 1 s . The step response in the Laplace domain, denoted as Y (s), is obtained by multiplying the unit step function by the first-order transfer function. To find the time-domain step response y(t), the inverse Laplace transform of Y (s) is computed: y(t) = L −1 { Y (s) } = L −1 {1 s · 1 Ts+ 1 } = 1 − e− t T When t = T , the system reaches one time constant, and the output attains 63.2% of its final value. At this point, y(t) = 1 − e−1 ≈ 0.632. 12 2. Theory Figure 2.4: Comparison of the step response of the actual ESC and the first-order model. The ESC is tasked with delivering 90% of the maximum torque specified by the motor manufacturer. 13 2. Theory 2.1.4 Sensor modeling The whole system is designed inside a simulation environment where a simulation scenario is inputted to the system and a sensor block will take the vehicle’s position and orientation at each time step and convert them to sensors’ signals. By having the sensors’ signals inside the simulation, it is possible to validate the sensor-fusion algorithms inside the simulation environment. Sensors are characterized on the basis of both gathered data and the data sheets provided by manufacturers. 2.1.4.1 Encoders An incremental encoder is a device mounted on the motor shaft that generates pulses corresponding to the shaft’s rotation. The resolution of the encoder is determined by the number of pulses produced per revolution. To accurately monitor the shaft’s position, the encoder count must be reset to a known reference point through a process known as homing. From this reference position, the shaft angle can be tracked by accumulating the encoder pulses. Figure 2.5: Encoder model in simulation. The model takes the true angle (in degrees) and converts it to encoder counts. Here, Max Counts denotes the total number of pulses generated by the encoder per revolution, as specified in the device datasheet. 2.1.4.2 Differential GNSS A Differential GNSS sensor is a dual GNSS module that provides information about the position on Earth and the absolute angle relative to the north. This is achieved by comparing the positions of two GNSS antennas mounted on opposite sides of the terminal. In simulation, position accuracy is assumed to be ideal, as translational errors have minimal impact on the system due to the significant distance between the antenna and the satellites. However, the heading angle estimation provided by the module is noisy and requires characterization and modeling within the simulation. By logging the heading estimation from the module over several hours, it was ob- served that the only significant error source was White Gaussian Noise (WGN) with 14 2. Theory a constant standard deviation over time. The standard deviation was determined from the logged data. Figure 2.6: Differential GNSS module model in simulation. It takes the position and the heading action points at each time step and converts them to signals from the GNSS sensor. The position measurements are assumed to be ideal, whereas for the heading estimation, the WGN is characterized based on analysis of logged sensor data. 2.1.4.3 Accelerometer A three-axis accelerometer is a device that measures linear acceleration along all three cartesian axes, x, y, and z, in its own body frame. When the accelerometer is at rest, it only measures the gravitational acceleration g. The effect of gravitational acceleration on the accelerometer depends on the orientation of the accelerometer, which is mounted on the outer gimbal. By knowing the relationship between the IMU body frame and the global frame (i.e., the angle α), the impact of the gravitational acceleration is calculated as follows: Accg = Accxg Accyg Acczg  = rotz(α) · −g · cos(θ) · sin(φ) g · sin(θ) g · cos(φ) · cos(θ)  (2.15) In addition to the gravitational acceleration, the linear acceleration caused by the motion of the vehicle is assumed to be negligible, as the system is primarily used when the vehicle is stationary or moving at a relatively constant velocity. Since the accelerometer is mounted with an offset from the rotational axis of the outer gimbal, the angular acceleration of the gimbal, α̈, introduces a linear accelera- tion in the body frame of the accelerometer, as illustrated in Figure 2.7. However, as the system actively maintains alignment with the satellite, the angular acceleration (α̈) counteracts the rotation of the vehicle (ψ̈), resulting in a net angular acceleration in the world frame. This introduces an induced acceleration Acci in the x direction of the accelerometer’s body frame. 15 2. Theory Figure 2.7: A top view of the terminal shown in Figure 2.2. When the car changes its heading rate, the antenna moves in the opposite direction to maintain line-of- sight with the satellite. The difference between these angular accelerations results in a net angular acceleration acting on the system, introducing an induced acceleration Acci in the x direction of the accelerometer’s body frame. The induced acceleration is also proportional to the distance, r, between the ac- celerometer position and the rotation axis. Acci = (α̈− ψ̈) · r 0 0  (2.16) Finally, a white gaussian noise (WGN), characterized by logging and analyzing accelerometer measurements, is added to the model. The total acceleration is then expressed as: Acctotal = Accg + Acci +WGN (2.17) Expanding equation 2.17 yields the following result: Acctotal = rotz(α) · −g · cos(θ) · sin(φ) g · sin(θ) g · cos(φ) · cos(θ) + (α̈− ψ̈) · r 0 0 +WGN 16 2. Theory Figure 2.8: The accelerometer model in simulation. As mentioned, external linear accelerations introduced by the car are not considered. 2.1.4.4 Gyroscope A three-axis gyroscope is a device that measures angular velocity along all three cartesian axes, x, y, and z, in its own body frame. Since the gyroscope is mounted on the outer gimbal, its readings in the body frame must be compensated for the gimbal angle α. Consequently, the gyro readings represent the derivative of the car’s orientation over time. Gyro = rotz(α) φ̇θ̇ ψ̇  (2.18) This situation is similar to the one described in Figure 2.7. A change in the car’s heading (ψ̇) is counteracted by the antenna’s motion in the opposite direction (α̇). This results in a net angular velocity, which is the difference between these two components. Gyro = rotz(α) φ̇θ̇ ψ̇ +  0 0 −α̇  (2.19) According to the gyroscope’s datasheet, calibration is necessary to reduce drift, which occurs when gyro data is integrated to calculate angles. Calibration is per- formed by locking the motors to ensure that the gimbal remains stationary, parking 17 2. Theory the vehicle, and calculating the mean value of the readings, referred to as the Bias. This bias is then subtracted from the gyro readings. The bias is unstable and varies over time. A high-quality gyroscope exhibits good bias stability, meaning the bias can be effectively compensated. Bias instability can be modeled as a constant value added to the gyro readings. To measure bias instability, the gyroscope was operated in a stationary state while its data was recorded. From this recorded data, the bias instability was quantified, revealing the drift in degrees per second. Finally, a white gaussian noise (WGN) is added to the model, which is also char- acterized based on the recorded data. Gyro = rotz(α) φ̇θ̇ ψ̇ +  0 0 −α̇ +Bias instability +WGN (2.20) Figure 2.9: The gyroscope model in simulation. The term "R andom walk" is a standard term for WGN in gyroscope contexts. Both the bias instability and random walk are modeled based on the logged data. 18 2. Theory 2.2 Kinematics The antenna will be aligned with the satellite by first rotating around the vertical axis, starting from geographic north. This angle of rotation is called the azimuth. Then, the antenna will be tilted upward from the horizon to point towards the satellite, and this angle is known as the elevation [2]. Once the antenna is aligned with the satellite, it must rotate around the line of sight between the antenna and the satellite to match the polarization of the electromagnetic signal transmitted by the satellite. Instead of physically rotating the antenna to align with the satellite’s polarization, the signal processing hardware can compensate for the polarization mismatch through software adjustments. By knowing both the terminal’s position and the satellite’s position, the global azimuth and elevation target angles can be calculated, while the polarization target angle is provided by the satellite. Addi- tionally, by knowing the vehicle’s orientation (roll, pitch, yaw), reference values for the α and β angles can be determined. Figure 2.10: The function Global Target calculates the target elevation, azimuth, and polarization angles in a global reference frame, based on the positions of the terminal and the satellite. Meanwhile, the function Target Joint Values computes the local joint angles β, λ and α required to adjust for the terminal’s orientation, ensuring the antenna is correctly aligned with the satellite. 2.2.1 Global Target The satellite, with which the antenna needs to align, is located at a fixed position above the equator. In simulation, both the vehicle’s and the satellite’s positions are expressed in cartesian coordinates, (x, y, z), whereas in the real world, they are expressed in spherical coordinates (r, longitude, latitude). • In the real world: Given the vehicle’s position at (r0, lonv, latv) and the satel- lite’s position at (r0 + h, lons, lats), where: – lonv is the vehicle’s longitude, while latv denotes the vehicle’s latitude. The same is for the satellite, expressed as lons and lats respectively. – r0 is the Earth’s radius. – h is the satellite’s altitude (height above sea level). 19 2. Theory – The latitude of the satellite, lats , is 0 degrees since all geostationary satellites are positioned above the equator in the geostationary orbit. The global target elevation and azimuth angles are given by Maral et al, [2] , as following : Elevation = tan−1 (cos(latv) · cos(lonv − lons) − r0 r0+h√ 1 − (cos(latv) · (lonv − lons))2 ) (2.21) a = tan−1 ( tan(|lonv − lons|) sin(latv) ) Satellite EAST of antenna Satellite WEST of antenna North Hemisphere Azimuth = 180 − a Azimuth = 180 + a South Hemisphere Azimuth = a Azimuth = 360 − a • In simulation: Calculating the satellite’s orientation relative to the vehicle (or terminal’s base) can be done using the following equations, where the variables correspond to the figure 2.11. Azimuth = atan2(∆y,∆x) (2.22) The azimuth angle represents the satellite’s position in the xy-plane relative to the terminal base. Elevation = atan2(∆z, L) (2.23) The elevation angle represents the vertical angle between the terminal base and the satellite, indicating how high the satellite appears in the sky relative to the horizon from the terminal’s perspective. It is calculated based on the difference in altitude (∆z) and the horizontal distance (L) between the terminal and the satellite. Figure 2.11: An illustration depicting how azimuth and elevation angles are cal- culated. 20 2. Theory 2.2.2 Target joint values The orientation of an object can be described using three rotational matrices rep- resenting rotations about the x, y, and z axes. Given the known orientation of the satellite in the global frame and the estimated orientation of the terminal base, it is possible to calculate the reference angles β, λ, and α using the following formulas: OT arget = OBase ·OAntenna (2.24) OAntenna = O−1 Base ·OT arget (2.25) rotz(α) · rotx(β) · roty(λ) = O−1 Base ·OT arget cα cλ− sα sβ sλ −sα cβ cα sλ+ sα sβ cλ sα cλ+ cα sβ sλ cα cβ sα sλ− cα sβ cλ −cβ sλ sβ cβ cλ  = R11 R12 R13 R21 R22 R23 R31 R32 R33  (2.26) From equation 2.26, it is possible to solve: • β = sin−1(R32) • λ = atan2(−R31 R33 ) • α = atan2(−R12 R22 ) 21 2. Theory 2.3 Processing step for sensors The data collected from encoders, gyroscope, accelerometer, and differential GNSS must be processed before being fused to accurately estimate the vehicle’s orienta- tion at each time step. The processing of sensors data is essential to extract roll, pitch, and yaw angles, which represent the vehicle’s orientation. The encoders, gy- roscope, and accelerometer data require specific processing steps, while the heading measurements from the differential GNSS sensor do not need any processing since they directly provide the estimated yaw angle of the vehicle. The following sections explain the processing steps for each sensor in detail. 2.3.1 Encoder Data Processing The encoder measurements are recorded in counts, representing the rotational posi- tion of a shaft. To calculate the angle α (in degrees) from these counts, a straight- forward conversion is applied. Since one revolution equals 360 degrees, the angle is determined by dividing 360 by the maximum number of counts per revolution (as specified in the encoder’s datasheet) and then multiplying the result by the encoder’s measured counts. This process is illustrated in Figure 2.12. Figure 2.12: Encoder processing step to estimate the angle α as a block diagram. 2.3.2 Gyroscope Data Processing The gyroscope provides angular velocity measurements and is affected by rotations due to the angle α as explained in section 2.1.4.4. To accurately estimate roll, pitch, and yaw angles, this influence of α must first be removed. The corrected data is then rotated back along the z-axis to eliminate any effects from α. Afterward, the angular velocity data are integrated to compute the orientation angles. This process is visualized in Figure 2.13. 22 2. Theory Figure 2.13: Gyroscope processing step to estimate roll, pitch, and yaw angles as a block diagram. The estimated α is differentiated and subtracted from the gyroscope’s z-axis measurements. The data is then rotated back along the z-axis by α to isolate the raw measurements, which are subsequently integrated to estimate the orientation angles. 2.3.3 Accelerometer Data Processing The accelerometer measures linear accelerations along the x, y and z axes. However, these measurements are influenced by additional effects from the angles α and ψap as explained in section 2.1.4.3. These effects must be removed, and the data must be rotated back along the z-axis to provide usable measurements. After removing all the additional effects, the accelerometer data looks as in the equation 2.27: Accg = Accxg Accyg Acczg  = −g · cos(θ) · sin(φ) g · sin(θ) g · cos(φ) · cos(θ)  (2.27) From this processed data, it is possible to directly estimate the roll (φ) and pitch (θ) angles. The roll angle (φ) is estimated as in equation 2.28. Using the Pythagorean trigonometric identity, it is possible to estimate the pitch angle (θ) as well, as shown in equation 2.29. φest = − arctan (accx accz ) (2.28) θest = arctan  accy√ acc2 x + acc2 z  (2.29) The processing workflow is depicted in Figure 2.14. 23 2. Theory Figure 2.14: Accelerometer processing step to estimate roll and pitch angles as a block diagram. The accelerometer data is processed by removing the effects of angles α and ψest and rotating the data back along the z-axis. The roll and pitch angles are then calculated using the equations above. 2.4 Sensor Fusion After processing the sensor data, the next step is to fuse the measurements from all sensors. Sensor fusion aims to combine the strengths of individual sensors while compensating for their limitations. For this research, two prominent fusion algo- rithms were considered: the complementary filter and the Kalman filter. These algorithms were tested and compared to evaluate their performance in estimating the vehicle’s orientation. 2.4.1 Complementary filter The complementary filter is a simple yet effective sensor fusion technique that com- bines the outputs of two sensors using a low-pass and a high-pass filter [6]. The low-pass filter is applied to the accelerometer data because accelerometers are prone to high-frequency noise. By suppressing this noise, the low-pass filter helps retain the reliable low-frequency components of the accelerometer readings. On the other hand, a high-pass filter is applied to the gyroscope data because gyroscopes pro- vide accurate measurements of angular velocity, particularly for detecting sudden motions. The high-pass filter ensures that these rapid changes are captured while mitigating the effects of low-frequency drift, which is a common limitation of gyro- scopes. The complementary filter equation in the continues time domain is derived as follows [6]: 24 2. Theory θest = Glpf(s) · θa +Ghpf(s) · θg (2.30) Glpf(s) = 1 1 + TCs , Ghpf(s) = TCs 1 + TCs θest = 1 1 + TCs · θa + TCs 1 + TCs · 1 s · ωg (1 + TCs) · θest = θa + TC · ωg (2.31) Here, the variables and filters are defined as follows: • Glpf(s): Low-pass filter for accelerometer data. • Ghpf(s): High-pass filter for gyroscope data. • TC : Time constant of the filter. • θa: Angle estimate from accelerometer. • θg: Integrated angle from gyroscope. • ωg: Angular velocity from gyroscope. • θest: Resultant angle estimation from the complementary filter. The next step involves converting Equation 2.31 from Laplace domain to the time domain, resulting in: θest(t) + TC · dθest(t) dt = θa(t) + TC · ωg(t) (2.32) To implement the filter in discrete time, Equation 2.32 is discretized, assuming a sampling time of Ts: dθ(t) dt → θ(n) − θ(n− 1) Ts θest(n) + TC · θest(n) − θest(n− 1) Ts = θa(n) + TC · ωg(n) θest(n) = TC TC + Ts · θest(n− 1) + Ts TC + Ts · θa(n) + TC · Ts TC + Ts · ωg(n) (2.33) To simplify the notation, the complementary filter constants are defined as: γ = TC TC + Ts , 1 − γ = Ts TC + Ts 0 ≤ γ ≤ 1 Substituting these constants, the final complementary filter equation is given as: θest(n) = γ · [θest(n− 1) + Ts · ωg(n)] + (1 − γ) · θa(n) (2.34) 25 2. Theory The constant γ determines the weighting of gyroscope and accelerometer data in the complementary filter. A value of γ = 1 implies complete reliance on gyroscope measurements, while γ = 0 means only accelerometer data is used. In practice, γ is chosen between these extremes to balance the noise characteristics and drift of the sensors, resulting in a robust angle estimation. 2.4.1.1 Application in This Research One of the main limitations of the complementary filter, as applied in this research, is its inability to estimate the yaw angle effectively using accelerometer data. Unlike pitch and roll, the yaw angle cannot be derived from accelerometer measurements, as it depends primarily on the horizontal plane, which accelerometers cannot di- rectly measure. Consequently, the yaw angle must be estimated solely through the integration of gyroscope measurements. However, this method leads to drift over time, resulting in an inaccurate yaw estimation. Additionally, the complementary filter is inherently limited to combining data from only two sensors, as shown in Equation 2.34. Since the constants γ and 1 − γ are directly applied to the estimations from these two sensors, any scenario where one of the sensors becomes unavailable or provides intermittent data leads to erroneous outputs. This limitation is particularly problematic in this research for yaw angle estimation. In such cases, the filter relies solely on the integration of gyroscope data, which worsens the drift problem over time due to the lack of corrective input from a secondary sensor. 2.4.2 Kalman filter The Kalman filter is a special case of Bayesian filtering where the process and measurement models are linear and Gaussian [7]. It has long been regarded as the optimal solution for many tracking and data prediction tasks, due to its ability to minimize the mean squared error. For the state vector xk (which represents the variables we want to estimate at time k) and the measurement yk (which represents sensors observations), the Kalman filter operates under the following process and measurement models [7]: xk = Ak−1xk−1 + qk−1, qk−1 ∼ N (0, Qk−1) (2.35) yk = Hkxk + rk, rk ∼ N (0, Rk) (2.36) Here: • xk: The state vector at time k, containing the quantities being estimated. • Ak−1: The state transition matrix, which describes how the states evolve from k − 1 to k based on the system dynamics. • qk−1: The process noise, modeled as Gaussian noise with zero mean and co- variance Qk−1, representing uncertainty in the process model. 26 2. Theory • yk: The measurement vector at time k, containing sensors observations. • Hk: The measurement model matrix, which maps the state xk to the expected measurement yk. • rk: The measurement noise, Gaussian noise with zero mean and covariance Rk, representing uncertainty in the measurements. The Kalman filter consists of two main steps: Prediction and Update. 2.4.2.1 Prediction Step The prediction step estimates the state and covariance at time k using information from time k − 1 [7]: x̂k|k−1 = Ak−1x̂k−1|k−1 (Predicted state estimate) (2.37) Pk|k−1 = Ak−1Pk−1|k−1A ⊤ k−1 +Qk−1 (Predicted state covariance) (2.38) Here: • x̂k|k−1: The predicted state at time k given the state estimate from time k−1. • Pk|k−1: The predicted covariance matrix, representing the uncertainty of the state prediction. 2.4.2.2 Update Step The update step corrects the predicted state and covariance using the new measure- ment yk [7]: x̂k|k = x̂k|k−1 +KkVk (Updated state estimate: Prediction + Correction) (2.39) Pk|k = Pk|k−1 −KkSkK ⊤ k (Updated covariance matrix) (2.40) Here: • x̂k|k: The updated state estimate at time k. • Pk|k: The updated covariance matrix, representing the uncertainty after incor- porating the measurement. • Kk: The Kalman gain, which determines the weight given to the new mea- surement compared to the prediction. • Vk: The innovation or residual, which is the difference between the actual measurement yk and the predicted measurement. • Sk: The innovation covariance, which measures the uncertainty in the innova- tion. 27 2. Theory The intermediate variables are defined as: Kk = Pk|k−1H ⊤ k S −1 k (Kalman Gain) (2.41) Vk = yk −Hkx̂k|k−1 (Innovation: Measurement residual) (2.42) Sk = HkPk|k−1H ⊤ k +Rk (Innovation Covariance: Total uncertainty in the residual) (2.43) The Kalman gain Kk determines how much the measurement yk influences the up- dated state. A higher Kalman gain gives more weight to the measurement, while a lower gain prioritizes the prediction or the process. 2.4.2.3 Application in This Research In this research, the Kalman filter was designed using a constant turn rate model, which includes the orientation angles (roll, pitch, and yaw) and their angular rates as states. A constant turn rate model was chosen because orientation changes in a vehicle typically occur gradually, making it reasonable to assume that angular accelerations are negligible. This assumption implies that the angular rates of the orientation angles remain approximately constant over time, with only minor vari- ations due to process noise. With this model as the process model, the Kalman filter state vector comprises six states: the orientation angles (roll, pitch, and yaw) and their respective angular rates. The angular rates are modeled as constant values perturbed by process noise, while the orientation angles are updated by integrating the angular rates over each time step, also with the addition of process noise. This results in the following state-space representation: φk φ̇k θk θ̇k ψk ψ̇k  =  1 Ts 0 0 0 0 0 1 0 0 0 0 0 0 1 Ts 0 0 0 0 0 1 0 0 0 0 0 0 1 Ts 0 0 0 0 0 1  ︸ ︷︷ ︸ A  φk−1 φ̇k−1 θk−1 θ̇k−1 ψk−1 ψ̇k−1  + qk−1, qk−1 ∼ N (0, Qk−1) (2.44) where φ, θ, and ψ represent the roll, pitch, and yaw angles, respectively; φ̇, θ̇, and ψ̇ are their angular rates; Ts is the sampling time; and qk−1 is the process noise, which is modeled as Gaussian with zero mean and covariance Qk−1. The tuning of the process noise covariance matrix Qk−1, both in general and in this research, is discussed in Section 2.4.2.4. Using this process model, two Kalman filter variants were developed to estimate the vehicle’s orientation. These variants differ in their measurement models and are detailed in the following sections. 28 2. Theory 2.4.2.3.1 Kalman based on IMU and differential GNSS This filter integrates measurements from a 6-axis IMU (gyroscope and accelerome- ter) and a differential GNSS sensor. The differential GNSS sensor provides heading measurements (yaw angle), which is used to refine the orientation estimation and mitigate IMU drift over time. The available measurements in this variant include roll and pitch angles derived from the accelerometer (discussed in Section 2.3.3), the yaw angle from the differential GNSS sensor, and the angular rates from the gy- roscope. These measurements are combined to produce the following measurement model:  φmk φ̇mk θmk θ̇mk ψmk ψ̇mk  =  1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1  ︸ ︷︷ ︸ H  φk φ̇k θk θ̇k ψk ψ̇k  + rk, rk ∼ N (0, Rk) (2.45) Here: • φmk , θmk , and ψmk represent the measured roll, pitch, and yaw angles, respec- tively. • φ̇mk , θ̇mk , and ψ̇mk represent the measured angular rates. • φk, θk, and ψk are the true roll, pitch, and yaw angles, while φ̇k, θ̇k, and ψ̇k are the true angular rates. • rk is the measurement noise, assumed to follow a Gaussian distribution with zero mean and covariance Rk. 2.4.2.3.2 Kalman Filter Based on IMU-Only This variant of the Kalman filter relies exclusively on measurements from the 6-axis IMU, which includes a gyroscope and an accelerometer. The gyroscope provides angular velocity data, enabling the estimation of changes in orientation over time, while the accelerometer is used to estimate the roll and pitch angles. In the absence of heading (yaw angle) measurements from a GNSS sensor, this filter estimates the yaw angle using only gyroscope data. However, gyroscope measure- ments are susceptible to drift over time due to cumulative integration errors, making the yaw angle estimation prone to long-term inaccuracies. In addition, this filter re- quires an initial sky-scanning procedure to establish a reference or starting position relative to a known heading (e.g., relative to a satellite’s position). 29 2. Theory The measurement model used for this filter is expressed as follows:  φmk φ̇mk θmk θ̇mk ψ̇mk  =  1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1  ︸ ︷︷ ︸ H  φk φ̇k θk θ̇k ψk ψ̇k  + rk, rk ∼ N (0, Rk) (2.46) Here, the measurement model matrix H directly relates the measured roll (φmk ), pitch (θmk ), and their angular velocities to the state variables. Since yaw (ψk) is not directly observable without GNSS data and only the yaw rate (ψ̇mk ) is directly observable by the gyroscope, the measruement model matrix H is a 5x6 matrix and differs from the previous filter’s measurement model matrix. 2.4.2.4 Tuning the Complementary and Kalman Filters The complementary filter has only one tuning parameter, γ, with its complement 1 − γ. Tuning this filter primarily involves understanding the characteristics of the two sensors being fused. If the gyroscope estimation is highly accurate but suffers from drift over time, γ is chosen to be close to 1, placing greater trust in the gyroscope’s short-term estimates while using the accelerometer primarily to correct for drift. In contrast, tuning the Kalman filter is more complex due to its reliance on covari- ance matrices to model process and measurement uncertainties. Specifically, the process covariance matrix Q is the primary tuning parameter. Each state in the Kalman filter requires careful adjustment of its corresponding variance in Q. For instance, when the GNSS sensor is identified as a noisy measurement source, and the gyroscope measurements are observed to be more reliable, the following adjustments are made: • The process model for states influenced by the GNSS sensor, such as the yaw angle, is given higher trust compared to the sensor’s measurements. This ensures that the filter relies more on the model’s predictions than on the noisy GNSS data. • Conversely, the states measured by the gyroscope, such as roll angular rate, are weighted more heavily in the filter compared to their corresponding process models. This allows the filter to leverage the gyroscope’s reliability for these states. 30 2. Theory 2.5 Control Loop Before discussing the control design, it is important to summarize the modeling of GGimbal, which represents the complete mechatronic system to be controlled. The model GGimbal includes the following components: • ESC and Motor: The ESC saturates the maximum requested torque Treq to a value derived from the motor’s data sheet to prevent damage to the motor. The delay in the ESC is modeled as a first-order transfer function with a time constant T , as shown in Equation 2.14. • Friction: Friction is modeled as a function of angular velocity with parameters a and b, as shown in Equation 2.13. Here, a describes the Viscous friction, which relates angular velocity to friction torque, while b describes Coulomb friction. • Equations of Motion: The relationship between the effective torque Teff and the angle is described by the equations of motion. The time-domain model is shown in Equation 2.12. By applying the Laplace transform, the model becomes: α(s) = 1 Jtotal s2 Teff (s) β(s) = 1 Jβ s2 Teff (s) (2.47) Figure 2.15: The simulation model representing the dynamics from requested torques to angles. The model includes two uncoupled systems, one for each axis. The transfer function GGimbal(s), which describes the relationship between the re- 31 2. Theory quested torque and the resulting angles, can be derived as follows: Figure 2.16: A simplified version of the block diagram shown in Figure 2.15, focusing on the relationship between a single requested torque and its resulting angle. For simplicity, the loop transfer function depicted in blue is denoted as Gsub(s), while the transfer function for the saturation is denoted as F (s). The output angle can be expressed as: Output = 1 s ·Gsub(s) · Tmotor where: Gsub(s) · Tmotor = −b a+ Js + 1 a+ Js · Tmotor Substituting this expression into the output equation gives: Output = 1 s · ( −b a+ Js + 1 a+ Js · Tmotor ) The motor torque Tmotor is given by: Tmotor = 1 Ts+ 1 · F (s) · Input Substituting this into the output equation results in: Output = 1 s · { −b a+ Js + 1 a+ Js · ( 1 Ts+ 1 · F (s) · Input )} Simplifying the equation further yields: Output = −1 Js2 + as · b+ 1 TJs3 + (Ta+ J)s2 + as F (s) · Input (2.48) 32 2. Theory Experimental results indicate that the friction parameters a and b are very small, and the time constant T of the ESC delay is considered to be negligible, as shown in Figure 2.4. Therefore, for control design purposes, it is possible to make the following assumptions: • The system does not reach saturation; therefore, F (s) = 1. • The ESC response is very fast; hence, T → 0. • The impact of friction on the system is negligible; thus, a, b → 0. With these assumptions, Equation 2.48 is simplified to: Output Input = 1 Js2 . (2.49) Thus, the simplified transfer function GGimbal(s), used for control design, becomes: GGimbal(s) =  1 Jtotal s2 , Outer gimbal, 1 Jβ s2 , Inner gimbal. (2.50) The block diagrams in Figures 2.15 and 2.16 present a simplified representation of how friction is modeled. However, the actual behavior of the friction torque is nonlinear. The sign of the friction torque, depends on the sign of ω because it always opposes the direction of rotation. Additionally, the friction torque acts as a reaction torque. In other words, if the motor torque is less than the friction torque, no effective torque will be applied to the system, Tmotor − Tf = 0. Tf (ω) =  (a |ω| + b), when ω > 0, −(a |ω| + b), when ω < 0, Tmotor, when |Tmotor| < a|ω| + b. (2.51) The controller’s primary objective is to ensure system stability and following a ref- erence point. According to control theory, a feedback system is stable if its poles are located in the left half-plane (LHP) of the complex plane [12]. Moving the poles further to the left increases the system’s responsiveness and speed, but at the cost of higher control effort. This can lead to saturation due to torque requests exceeding the motor’s maximum limits. By tuning the controller, the poles of the closed loop system will be adjusted. Placing the poles at the origin of the complex plane results in marginal stability. In this case, the system may oscillate indefinitely without damping, which is undesirable for most practical applications. Therefore, careful pole placement is crucial to balance stability, responsiveness, and control signal limitations. 33 2. Theory 2.5.1 PID Controllers A PID controller is a widely used type of closed-loop controller [8]. It consists of three components: • Proportional (P): This component applies a control action proportional to the error signal. It helps reducing the error by providing a direct response to the difference between the desired and actual values. • Derivative (D): This component corresponds to the derivative of the error signal. It anticipates changes in the error and improves the system’s respon- siveness, helping to make its reaction quicker and more stable. • Integral (I): This component is proportional to the accumulation of the error over time. It addresses steady-state errors by ensuring that any persistent deviation from the desired value is eliminated. The combination of these three terms allows the PID controller to provide a balance between responsiveness, stability, and steady-state accuracy, making it effective for a wide range of control applications. 34 2. Theory 2.5.1.1 P Controller When the system is controlled by a proportional (P) controller, the open-loop trans- fer function H(s) becomes: H(s) = GController ·GGimbal = kp · 1 Js2 (2.52) Figure 2.17: Closed-loop control system for GGimbal using a P-controller with gain kp. This leads to a closed-loop transfer function L(s) as follows: L(s) = H(s) 1 +H(s) = kp · 1 Js2 1 + kp · 1 Js2 = kp Js2 + kp (2.53) The stability of the system with a proportional gain kp can be determined by ana- lyzing the poles of the closed-loop transfer function. The poles are the roots for the characteristic equation, (CE), i.e., the denominator of L(s): Poles = ± √ −kp J (2.54) From this equation, the stability depends on the value of kp: • If kp > 0, the poles are purely imaginary and lie on the imaginary axis of the complex plane. This results in a marginally stable system. • If kp < 0, one pole is in the Left Half-Plane (LHP), and the other is in the Right Half-Plane (RHP). This results in an unstable system. 2.5.1.2 D Controller When the system is controlled by a derivative (D) controller, the open-loop transfer function H(s) is: H(s) = GController ·GGimbal = s kd · 1 Js2 (2.55) 35 2. Theory Usually, when using a D-controller, a low-pass filter is added to the error signal. This prevents the derivative action from amplifying high-frequency noise, which could cause instability or erratic behavior. However, for stability analysis, the filter can be excluded. Figure 2.18: Closed-loop control system for GGimbal using a D-controller with gain kd. The resulting closed-loop transfer function L(s) is given by: L(s) = H(s) 1 +H(s) = s kd · 1 Js2 1 + s kd · 1 Js2 = kd Js+ kd (2.56) To determine the system’s stability with the derivative gain kd, the poles of the closed-loop transfer function are analyzed. The poles are the roots for the CE: Poles = −kd J (2.57) • If kd > 0, the pole lies in the Left Half-Plane (LHP), resulting in a stable system. • If kd < 0, the pole lies in the Right Half-Plane (RHP), resulting in an unstable system. 2.5.1.3 PD Controller By combining P and D components in a PD controller, the open-loop transfer func- tion H(s) becomes: H(s) = GController ·GGimbal = (kp + s kd) · 1 Js2 (2.58) 36 2. Theory Figure 2.19: Closed-loop control system for GGimbal using a PD-controller. The resulting closed-loop transfer function L(s) is: L(s) = kd s+ kp Js2 + kd s+ kp (2.59) The poles of the closed-loop system are determined by the roots of the CE: Poles =  p1 = − kd 2J + √ k2 d 4J2 − kp J p2 = − kd 2J − √ k2 d 4J2 − kp J (2.60) To simplify the tuning process and ensuring stability, a constraint can be imposed to place the poles on top of each other. This means: p1 = p2 = −kd 2J ⇔ kp = k2 d 4J (2.61) By selecting kp and kd such that they satisfy Equation 2.61, the system will be critically damped with both poles located in LHP, ensuring stability. 2.5.1.4 PID Controller Introducing only an integral (I) controller to the system results in the closed-loop transfer function: L(s) = ki s · 1 Js2 1 + ki s · 1 Js2 = ki Js3 + ki . (2.62) The poles of this system are given by the following: Poles = 3 √ −ki J . (2.63) This results in three poles: 37 2. Theory • If ki < 0, there is one real pole in the right half-plane (RHP) and two complex poles in the left half-plane (LHP). • If ki > 0, there is one real pole in the LHP and two complex poles in the RHP. Thus, the system is unstable for any value of ki. Adding both proportional and integral components (PI controller) leads to the fol- lowing closed-loop transfer function: L(s) = kp s+ ki J s3 + kp s+ ki . (2.64) Determining the poles directly is not straightforward; instead, the Routh-Hurwitz criterion is applied to identify the conditions for stability. The characteristic equa- tion (CE) is as follows: CE = J s3 + 0 s2 + kp s+ ki. (2.65) According to the Routh-Hurwitz criterion, a necessary condition for stability is that all coefficients of the characteristic equation must be strictly positive [12]. However, the coefficient of s2 is zero, violating this condition. Therefore, the system is not stable for any combination of kp and ki. Including all components in a PID controller results in the following closed-loop transfer function: L(s) = kd s 2 + kp s+ ki J s3 + kd s2 + kp s+ ki . (2.66) The characteristic equation (CE) is: CE = J s3 + kd s 2 + kp s+ ki. (2.67) Constructing the Routh-Hurwitz table for this characteristic equation: s3 J kp 0 s2 kd ki 0 s1 kdkp−Jki kd 0 s0 ki For stability, the first column of the Routh-Hurwitz table must be strictly positive. This results in the following conditions: • kd > 0 and ki > 0 38 2. Theory • kp > J ki kd These conditions ensure that all the roots of the characteristic equation have negative real parts, making the system stable. 2.5.2 LQR Controllers Linear Quadratic Regulator (LQR) is an optimal full-state feedback linear controller. An LQR controller is considered optimal because it minimizes a trade-off between the control error and the control effort by using weight matrices [9], Q1 for the states and Q2 for the inputs. These weight matrices define the relative importance of the states and the control effort in the system. The optimization problem is mathematically expressed as [9]: min ( ∥e∥2 Q1 + ∥u∥2 Q2 ) , (2.68) where u represents the control action. In this report, the control action corresponds to the requested torque, Treq. In a proportional controller (e.g., Figure 2.17), the control signal is expressed as: Treq = kp · e, where e is the difference between the reference and estimated angles. In contrast, a full-state feedback controller uses a vector of state errors. For the gimbal system, the states include angles and their derivatives, resulting in a richer feedback mechanism. To design the LQR controller, the system dynamics is usually expressed in state- space form: ẋ(t) = Ax(t) +B u(t), y(t) = C x(t), (2.69) where x(t) is the state vector, u(t) is the input vector, and y(t) is the output vector. The controller gain matrix, kLQR, that satisfies the optimization criteria in Equation 2.68, is given by [9]: kLQR = Q−1 2 BT S, (2.70) where S is a symmetric semipositive matrix obtained by solving the Algebraic Riccati Equation (ARE) [9]: ATS + SA+ CTQ1C − SBQ−1 2 BTS = 0. (2.71) Using Equation 2.70, the resulting kLQR guarantees stability with all poles in the stability region [9]. These poles can be computed as follows: Poles = eig(A−B kLQR). 39 2. Theory 2.5.2.1 LQR model The system to be controlled is described in the Laplace domain (see Equation 2.50). By applying the inverse Laplace transform, the original variables can be expressed as:  Treqα = Jtotal · ω̇α, ωα = α̇, Treqβ = Jtotal · ω̇β, ωβ = β̇. (2.72) Let x represent the state vector, which includes [α, ωα, β, ωβ]′. The input vector, u, includes [Treqα , Treqβ ]′. The state-space model can then be written as: ẋ = Ax+Bu, (2.73) where:  α̇ ω̇α β̇ ω̇β  ︸ ︷︷ ︸ ẋ =  0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0  ︸ ︷︷ ︸ A  α ωα β ωβ  ︸ ︷︷ ︸ x +  0 0 1 Jtotal 0 0 0 0 1 Jβ  ︸ ︷︷ ︸ B [ Treqα Treqβ ] ︸ ︷︷ ︸ u . (2.74) 2.5.2.2 LQI model An LQR controller with integral action is referred to as LQI. The LQI controller can eliminate the steady-state error over time. By adding additional integrating states, the state-space model is extended to include the new integral states. The integral state is a variable that accumulates the error over time: Iα = ∫ (rα − α) dt Iβ = ∫ (rβ − β) dt ⇒ İα = (rα − α) İβ = (rβ − β) (2.75) the extended state-space model will be: 40 2. Theory  α̇ ω̇α İα β̇ ω̇β İβ  ︸ ︷︷ ︸ ẋe =  0 1 0 0 0 0 0 0 0 0 0 0 −1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 −1 0 0  ︸ ︷︷ ︸ Ae  α ωα Iα β ωβ Iβ  ︸ ︷︷ ︸ xe +  0 0 1 Jtotal 0 0 0 0 0 0 1 Jβ 0 0  ︸ ︷︷ ︸ Be [ Treqα Treqβ ] ︸ ︷︷ ︸ u +  0 0 0 0 1 0 0 0 0 0 0 1  [ rα rβ ] (2.76) 2.5.2.3 Tuning LQR The tuning of an LQR controller is achieved by adjusting the elements of the weight- ing matrices, Q1 and Q2. These matrices serve as design parameters that balance the trade-off between control accuracy and control effort [9]. • Q1 is a square matrix with a size equal to the number of states in the system. • Q2 is a square matrix with a size equal to the number of inputs. A common and practical approach to tuning is to assume that cross-coupling penal- ties are zero [9] , as long as the physical system does not impose explicit constraints (for exemple, limited power availability for simultaneously driving multiple actua- tors). In most cases, there is no significant cross-coupling, so both matrices can be assumed diagonal. Each diagonal element in Q1 corresponds to a penalty for a specific state, while each diagonal element in Q2 corresponds to a penalty for a specific input. Increasing a particular diagonal element in Q1 prioritizes minimizing the error in the correspond- ing state. This results in a controller that emphasizes accuracy in that state but may require higher control effort or compromise the performance of other states. Increasing a particular diagonal element in Q2 penalizes the corresponding control input, effectively making the input more "expensive". This leads to a slower system response, as the controller reduces its reliance on that input to avoid high control effort. By systematically adjusting these diagonal elements, the LQR controller can be tuned to achieve a desired balance between performance and control effort for the specific application. 41 2. Theory 42 3 Methodology This chapter presents the methodology employed for developing and validating the two-axis gimbal system. It begins with a comprehensive literature review that es- tablishes the theoretical and technical foundation for the research. The discussion then transitions to a detailed overview of the simulation environment, focusing on the tools and techniques used to model, simulate, and validate the system effectively. The chapter further delves into the modeling of the physical system using Simscape, leveraging an acausal approach to accurately capture the nonlinear dynamics of the gimbal system. Following this, the design of a rigorous validation trajectory is described, simulating diverse real-world road conditions and challenges to thoroughly evaluate the control system’s performance. Lastly, the chapter highlights the development of a 3D visualization environment, which provides an intuitive and immersive representation of the system’s behavior. This visualization enhances the validation process by offering a clear and dynamic perspective of how the system operates in real-world scenarios. 3.1 Literature Study This research relied exclusively on published articles from reputable journals, con- ference proceedings, and authoritative books as reference materials. Relevant litera- ture was identified through a systematic search for key topics directly related to the study’s objectives. Particular emphasis was placed on selecting works that provided theoretical and technical foundations for the modeling, simulation, and validation processes. The selected references played a critical role in guiding the research methodology, offering insights into best practices, established theories, and recent advancements in the field. This comprehensive literature review ensured that the study was grounded in a strong academic and technical framework, aligning it with state-of-the-art re- search. 43 3. Methodology 3.2 Simulation Environment To model the physical system of the two-axis gimbal and to design and implement the entire control loop, a robust simulation environment was essential. MATLAB, in conjunction with Simulink, was selected as the software suite to achieve these goals. The system parameters and required inputs, such as the fixed satellite position and test scenarios defining the vehicle’s orientation and position over time, were created in MATLAB. Since both MATLAB and Simulink are products of MathWorks, the seamless integration between the two proved highly advantageous, enabling efficient modeling and simulation of the system. Simulink provided an ideal platform for designing, simulating, and validating physi- cal systems, including control loops, due to its extensive library of built-in functions and its ability to integrate MATLAB scripts. The built-in visualization tools, such as Simulink Scopes and MATLAB plotting functions, were invaluable for monitoring system behavior throughout the research. These tools facilitated iterative refine- ment and optimization of the models and controllers, ensuring that the simulation environment evolved alongside the research. Figure 3.1 illustrates the complete simulation model implemented in Simulink, show- casing the integration of system components for modeling, control design, and valida- tion. The mechatronic system, represented by GGimbal(s), incorporates the equations of motion for the rotational components, friction modeling, and the dynamics of the motor driver (ESC). In addition, the simulation includes sensor models, enabling validation of the entire system in a virtual environment. This holistic approach allowed the simulation to accurately mimic real-world conditions and provided a reliable framework for testing and analysis. Figure 3.1: Overview of the simulation system implemented in Simulink. This diagram mirrors the conceptual block diagram shown in figure 2.1. 44 3. Methodology The simulation environment takes a predefined test scenario as input, expressed in terms of action points. A test scenario specifies the satellite’s fixed position, the vehicle’s position, and its orientation over time. Using these inputs, the global target parameters, elevation (el), polarization (pol), and azimuth (az), are computed via a function called "Global Target". With the global target angles and the vehicle’s orientation as inputs, the joint angles α and β, which correspond to the actuator states, are calculated. These angles define the positioning of the system’s actuators. Meanwhile, the vehicle states, including its orientation, are estimated using sensor signals generated within the simulation environment through a dedicated block called the "Estimator". Based on the estimated states derived from these sensor signals, a feedback control configuration computes and regulates the torques of each motor to achieve precise system operation. Additionally, the simulation environment includes a comprehensive 3D visualization, making it easier to observe the system’s behavior and validate its performance. The block labeled "3D" processes signals from the simulation environment to generate this visualization, leveraging Simscape Multibody for real-time rendering. This block is further detailed in Section 3.5, where its role in the visualization of the system is thoroughly explained. As the research progressed, the system grew increasingly complex. To maintain clarity and facilitate understanding, Simulink’s tools were utilized to partition the system into manageable subsystems. This modular approach not only enhanced the organization but also streamlined the development process, making the simulation environment more efficient and user-friendly. 3.3 System Modeling Using Simscape Simscape, an acausal and equation-based modeling language, was utilized to model the nonlinear dynamics of the two-axis gimbal system. Unlike traditional signal- flow modeling approaches, Simscape allows users to define system dynamics through equations, enabling the direct modeling of physical relationships. This feature makes it particularly effective for solving complex nonlinear systems, such as the two-axis gimbal. The use of Simscape was crucial in the early stages of the research, as it facilitated an accurate representation of the nonlinear dynamics required to understand the system’s behavior. Modeling the gimbal system by directly writing its equations of motion and solving them within the Simscape environment proved to be a more effective approach compared to signal-flow-based methods. This approach ensured that the physical system was represented with greater fidelity, providing valuable insights into the gimbal’s dynamic response. After confirming that Equations 2.6 and 2.8 could be expressed as Equation 2.11, the two-axis gimbal dynamics were implemented directly in Simulink. This approach was chosen because the dynamics, 45 3. Methodology having been simplified, no longer required the use of Simscape. 3.4 Trajectory Design for Validation Validating the performance of the feedback control system was a critical aspect of this research. To achieve this, a trajectory was designed and generated using MATLAB, simulating a variety of real-world road conditions. This trajectory was used to thoroughly test the controller and filtering algorithms under diverse and challenging scenarios. The trajectory incorporated road segments with varying pitch, roll, and yaw an- gles, as well as challenging features such as speed bumps, roadside inclines, and gravel roads. The design of these road types was a collaborative effort, involving brainstorming sessions within the research team and consultations with SATCUBE supervisors. This ensured the trajectory met operational requirements while repre- senting realistic driving conditions. The trajectory was shaped like a rectangle, measuring 480 meters in length and 180 meters in width. It featured four 90-degree turns with varying radius, ranging from 30 meters to 50 meters, to evaluate the controller’s ability to handle both sharp and smooth turns. Additionally, pitch and roll angles were designed to reflect realistic worst-case scenarios. The maximum pitch angle used in the trajectory was 20 degrees, based on industry benchmarks [10], while the maximum roll angle was set to 45 degrees. Although no specific source was available for the roll angle, it was chosen to be sufficiently high to represent a challenging scenario. Special features, including speed bumps, gravel roads, and roadside inclines, were designed to simulate worst-case conditions that a vehicle might encounter during operation. These features were introduced to assess how such challenges impact the performance of the control system. Additionally, the trajectory’s key parameters were made adjustable, allowing the trajectory to be easily modified if needed. The resulting trajectory is illustrated in Figure 3.2. 46 3. Methodology Figure 3.2: Trajectory created for validation of the feedback control system, fea- turing varying road conditions, including pitch, roll, speed bumps, roadside inclines, and gravel roads, with 90-degree turns of different radius. A critical consideration during trajectory development was the vehicle’s velocity, as it significantly influences the controller’s response. To address this, the trajectory was parameterized in MATLAB, enabling dynamic adjustment of the vehicle’s speed for each simulation run. The total duration of the run was subsequently calculated based on the selected velocity. Additionally, to assess the performance of the filtering algorithm and evaluate the potential impact of gyroscope drift over extended periods, a parameter was intro- duced to define the number of rounds the vehicle should complete before finishing the trajectory. This feature provided a valuable means of validating the filtering algorithm’s robustness and its ability to mitigate drift over time. 3.5 3D Visualization To enhance the validation of the control system, more than numerical or graphical signal observations were required. A visual representation of the system’s behavior provided an additional, intuitive layer of verification. To this end, a 3D simulation 47 3. Methodology environment was developed using CAD models imported into Simscape Multibody. This visualization features a vehicle with an antenna terminal mounted on its roof, the track the vehicle follows, and a satellite. A virtual line indicates the direc- tion the antenna is pointing, providing a clear visual representation of the system’s alignment. The 3D visualization proved invaluable in verifying the controller’s performance, particularly in relation to the kinematic behavior outlined in Section 2.2. It enabled the team to observe the vehicle navigating the trajectory while the gimbal system actively maintained alignment with the satellite, providing critical insights into the system’s real-world operation. To safeguard SATCUBE’s intellectual property, the antenna system’s 3D model was intentionally simplified. Its actual dynamics and detailed geometry were excluded while still retaining the essential features needed for functional validation. Despite this simplification, the visualization served as a highly effective tool for demonstrat- ing system functionality and validating the control system’s performance in a clear and comprehensible manner, as illustrated in Figure 3.3. Figure 3.3: The 3D simulation environment created using Simscape Multibody. The vehicle, equipped with the antenna system mounted on its roof is shown tracking the trajectory. The antenna remains aligned with the satellite and the line of site represented by the red line. 48 4 Results This chapter presents the outcomes of the study, encompassing the validation of the system modeling, sensor data processing, state estimation, and controller design evaluation. Through comprehensive simulations, the performance of the proposed approaches is analyzed and illustrated using relevant graphs and figures. These results validate the accuracy of the sensor models and demonstrate the effectiveness of the processing and estimation methods, including complementary and Kalman filtering techniques. Additionally, the controller design strategies, from PD and PID controllers to LQR and extended LQR approaches, are evaluated. Furthermore, this chapter provides the tuning parameters and some specific values for critical components that correspond to the physical prototype and hardware uti- lized by SATCUBE. These details bridge the gap between simulation and real-world implementation, ensuring the presented results align closely with the capabilities and constraints of the actual system. 4.1 Modeling Results This section evaluates the impact of several assumptions made in the theoretical analysis to simplify the model and enable controller design. Specifically, the evalu- ation focuses on the linear equation of motion and the transfer function that repre- sents the complete mechatronic system, including the ESC, friction, and the equation of motion. 4.1.1 The Linear Equations of Motion Using the equation shown in 2.11, the motion of the outer gimbal is nonlinear due to the variable Jtotalzz , which includes Jβα, a nonlinear function of the state β as shown in 2.10. Therefore, an assumption was made to simplify the system and derive the linear equation in 2.12. The results of simulating both the nonlinear model and the linear model are presente