Next Article in Journal
Analysis of Employees’ Competencies in the Context of Industry 4.0
Previous Article in Journal
Non-Renewable and Renewable Energies, and COVID-19 Pandemic: Do They Matter for China’s Environmental Sustainability?
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Quadrotor Model for Energy Consumption Analysis

Institute of Aeronautics and Applied Mechanics, Warsaw University of Technology, St. Nowowiejska 24, 00-665 Warsaw, Poland
*
Author to whom correspondence should be addressed.
Energies 2022, 15(19), 7136; https://doi.org/10.3390/en15197136
Submission received: 29 August 2022 / Revised: 17 September 2022 / Accepted: 21 September 2022 / Published: 28 September 2022
(This article belongs to the Section E: Electric Vehicles)

Abstract

:
In this paper, a quadrotor dynamic model’s energy efficiency was investigated. A method for the design of the dynamic model which assures energy consumption estimation was presented. This model was developed to analyze the energy efficiency of the quadrotor during each maneuver. A medium-class quadrotor (4.689 kg) was used as a test platform. Thrust force correction factors obtained with FLIGHTLAB software were used to predict object behavior in forward flight. Model validation and long-duration flight tests in outdoor windy conditions are also presented. Monte-Carlo simulation was used to study the influence of uncertainties in model parameters on the simulation reliability. The developed model might be used for practical purposes (for example, energy-efficient coverage path planning).

1. Introduction

In the last few decades, the number of unmanned aircraft operations has increased significantly. At the beginning, the Unmanned Aerial Vehicles (UAVs) were used mainly for military purposes. Today, plenty of different types of UAVs are performing missions which are too dull, dirty, dangerous, or just not cost-effective for manned aircrafts, as well as missions in the civil airspace. Rotorcrafts are a significant group of unmanned aerial vehicles used for civil purposes. Their vertical take-off, landing, and hover capabilities give them an advantage over the fixed-wing UAVs, but their energy efficiency is significantly smaller. One of the most popular types of UAV rotorcraft in the civil market is the quadrotor. There are several reasons for this fact. Firstly, the quadrotor is configuration-easy to design and build. Additionally, the control of quadrotors is simpler compared with other rotorcraft configuration control systems, the cross-couplings between the degrees of freedom are not as strong as for a single-rotor helicopter, and the problem is easy to cope with using simple classical control algorithms [1,2,3,4]. Quadrotors have very good handling qualities and are easy to fly. Their main advantages are stable behavior, low vibrations, quiet flight, and simple maintenance. These features make them perfect for use in observation system [5].
Small, unmanned quadrotors are powered by electric motors driving fixed-pitch propellers. This multi-rotor propulsion system provides the lift, thrust, and control of the quadrotor, which results in a significant demand for electric power. Electric power is also consumed by the onboard systems (navigation system, communication system, control system, and others). The only power source aboard the quadrotor is the battery, and the amount of energy stored is limited by its weight. That is why the most significant disadvantage of a quadrotor is limited available energy restricting its flight duration and range [6,7,8]. A typical flight time of such drones on a single battery (charge) is in the order of several minutes.
Quadrotor energy effectiveness might be improved by appropriate design, flying qualities, and mission planning. The structural design aims toward weight reduction and aerodynamic efficiency improvement, as well as propulsion system [3,9] and onboard equipment optimization with respect to the energy consumption. The following actions can be used to improve quadrotor energy effectiveness: designing a structure optimized for weight reduction, increasing aerodynamic efficiency, and optimizing the propulsion system [1,2] and onboard equipment for energy consumption. Another way to increase the flight duration and range is to minimize energy expenditure during the flight by optimizing the mission plan and flight trajectory [10,11]. To achieve the mission goals the drone has to realize a set of maneuvers (e.g., take-off, hover, turn, cruise, approach, and landing). Each maneuver costs a certain amount of energy. The mission scenario and flight profile might be designed in such a manner as to minimize energy consumption [4,12,13,14,15]. To realize this approach, a reliable, mathematical quadrotor simulation model is required [16,17].
The subject of energy consumption by multirotor UAVs has been addressed by several researchers. Zhang et al. [18] discussed the comparison of such models and concluded that there exists a need to validate such models through field tests. A detailed review of existing models was also presented by Beigi et al. [19]. The existing models might be categorized into three main groups.
First, the most common approach is to use models based on physics principles. Lu et al. [2] developed a dynamic model of quadrotors and electric motors. Morbidi et al. [20] presented an energy consumption model and demonstrated how it can be used to optimize a drone’s trajectory. Yage et al. [21,22] proposed an energy consumption model and realized trajectory optimization. Jee and Cho [23] studied the energy consumption of electric motors. The polynomial expression on the consumed energy was derived by Li et al. [24]. This approach allows a detailed understanding of quadrotor dynamics. However, the main disadvantage of these models lies in the fact that they require a set of parameters that are often difficult to obtain. Often, only the energy spent on propulsion is considered [25].
Second, up to this time, a number of black box models have been proposed [26,27,28]. Alyassi et al. [29] presented extensive experiments for three different drones and proposed a model that uses nine coefficients. Abeywickrama et al. [30,31] proposed a model that is a function of on-ground power consumption, communication activities, hovering, vertical and horizontal movements, speed, payload, and wind. Recently, machine learning techniques were also used by Steup et al. to create a generic energy model [32]. These methods do not require detailed knowledge about quadrotor parameters [29]. The main limitation of such methods is that they will not allow deep insight into the drone dynamics.
Third, aeromechanics are also used to predict the energy consumption [33,34,35,36]. Momentum Theory [37,38,39,40] and Blade Element Theory [41] are used to calculate aerodynamic power. These models allow taking into account the forward speed regarding the consumed energy [42]. The disadvantage of these models is the necessity of iteratively calculating the induced velocity that degrades the computational efficiency [43].
Simple models of the loads produced by a quadrotor propulsion system were used in most of the abovementioned works. Only thrust force and torque were modeled as functions of rotor angular rate, and do not provide a full picture of energy consumption. This method is only valid in hover conditions, but during the cruise flight, where the free stream through the rotor is more horizontal and the blade tip vortices are sparse, the energy required for rotor driving is reduced. That is why the influence of the forward velocity on the rotor loads and energy consumption must be considered. Moreover, in the literature, frequently, only the energy consumed by electric motors is considered and the study of energy consumption is limited mainly to indoor tests. This simplification might lead to the overestimation of energy consumption, especially in long-duration flights. To render the simulation more realistic, the battery performance should be taken into account, as omitting battery dynamics might lead to significant errors. It is common that only the simulation results of energy consumption are presented without comparing the proposed models with real data. Eventually, the validation is performed using only data from a few flight tests [44]. Due to these facts, it is difficult to determine the reliability of existing energy consumption models. Even the common dataset various models could provide very significant differences in the amount of predicted energy. The dynamic model of quadrotor, which utilizes a detailed model of energy consumption, is needed in order to plan missions and analyze and synthesize their energy costs.
The aim of the present work was to develop and validate a quadrotor dynamics model for the analysis of energy consumption during long-duration outdoor flight. A six degrees-of-freedom nonlinear dynamic model with a reliable model of energy consumption was derived. The model was developed in the MATLAB/Simulink environment, which facilitates easy design and simulation of additional features such as atmospheric disturbance, control system, and experimental data processing. A detailed model of the rotor, which was designed with use of an advance modeling tool for rotorcraft dynamic simulation FLIGHTLAB, was used to estimate the velocity correction factors of propulsion system loads simulated in the main model. The model was validated using the experimental data that were acquired during laboratory and in-flight tests. The capability of the model to simulate energy consumption was analyzed as well.
The organization of the remaining portion of the paper is as follows. In Section 2, the quadrotor dynamics, energy consumption model, and battery model are described. In Section 3, the results of the flight tests for two example trajectories are presented. Then, flight data are compared with simulation results. Monte-Carlo simulation is used to investigate the model robustness on parameter uncertainties. The manuscript ends with a summary of the main findings and suggestions for possible further research directions.

2. Materials and Methods

2.1. Quadrotor Used as a Test Platform

The modeled UAV was a medium class, off-the-shelf M690A quadrotor (Figure 1).
The parameters of the quadrotor model were obtained by measurements and from the producer’s manual. The laboratory tests were initially aimed at estimation of the mass and geometry parameters of the quadrotor. Moments of inertia were measured experimentally using a trifilar pendulum (Figure 2). The above-mentioned mass parameters were confirmed using a quadrotor 3D CAD model.
The resulting parameters are presented in Table 1.
Four brushless direct current electric motors (BLDC) directly drive carbon fiber propellers. The drone structure is highly integrated and optimized to achieve the maximum flight duration. It should be mentioned that the drone is able to fly continuously for more than one hour (please see Table 1), which makes it suitable for inspecting large areas. The fuselage is composed of aluminum alloys to reduce its mass. The disadvantage of this design approach is that the replacement of individual components (e.g., motors) is quite difficult. The drone is equipped with an Orange Cube Flight Controller [45] containing several redundant measurement devices: ICM20948 (9 axis IMU: 3× accelerometers, 3× gyroscopes, 3× magnetometers), ICM20649 (6 axis IMU), ICM20602 (6 axis IMU), and two MS5611 barometric pressure sensors. The drone instrumentation also includes a Here 3 GPS GNSS module, RFD868x telemetry modem, and Raspberry Pi. The onboard equipment facilitates obtaining a set of data pertaining to the quadrotor state (e.g., drone velocity, position, angular rates, accelerations, Euler angles, battery parameters, etc.). The quadrotor is additionally equipped with a gimbaled digital camera.

2.2. Quadrotor Nonlinear Dynamic Model

The quadrotor is modeled as a rigid body with six degrees of freedom and constant mass. The lift, thrust, and control forces and moments are produced by four rotors driven by the electric motors. The Earth’s rotation effects were omitted due to the short flight range.
The quadrotor equations of motion are derived in the body coordinate system O b x b y b z b (Figure 3) fixed to the airplane’s fuselage. The center O b of the system was placed at the UAV gravity center. The O b x b axis lays in the plane of quadrotor symmetry and was directed forward. The O b y b axis was perpendicular to the plane of symmetry and pointed right, while the O b z b axis pointed “down”.
The translations and attitude angles were calculated in the inertial coordinate system O n x n y n z n ; the center of this system O n was placed at an arbitrary point on the earth surface. The O n z n axis was along the vector of gravity acceleration, and it pointed down. The O n x n z n plane was horizontal, tangent to the Earth’s surface, and the O n x n axis pointed to the north and O n y n the axis to the east.
The vector y = x n y n z n Φ Θ Ψ T defines the position and attitude of the quadrotor (Figure 3). It is composed of the vector of the position r n = x n y n z n T in the ground system of coordinates O n x n y n z n and roll Φ , pitch Θ , and yaw Ψ angles describe its attitude. The quadrotor state vector x = v ω T is composed of linear velocity v = U V W T and angular rate ω = P Q R T components.
The vectors of quadrotor states, position, and attitude are related by the following kinematic equation:
y ˙ = T x .
The matrix T is composed of two matrices: T V relating to velocities and T Ω relating to angular rates:
T = T V 0 0 T Ω ,
where [46,47]:
T V = cos Θ cos Ψ sin Φ sin Θ cos Ψ cos Φ sin Ψ cos Φ sin Θ cos Ψ + sin Φ sin Ψ cos Θ sin Ψ sin Φ sin Θ sin Ψ + cos Φ cos Ψ cos Φ sin Θ sin Ψ sin Φ cos Ψ sin Θ sin Φ cos Θ cos Φ cos Θ ,
and [46,47]:
T Ω = 1 sin Φ tan Θ cos Φ tan Θ 0 cos Φ sin Φ 0 sin Φ sec Θ cos Φ sec Θ .  
The quadrotor equations of motion were obtained by summing up the forces and moments from inertia, gravity f G , aerodynamic f A , and rotor f R loads [47]:
A x ˙ + B x x = f G y + f A x , y + f R x , y , Ω R ,
where Ω R is the vector of the angular velocity of the rotors.
Matrix A describes the inertia properties of the quadrotor, and matrix B x = Ω x A results from the inertia loads not depending on accelerations [47]:
A = m 0 0 0 0 0 0 m 0 0 0 0 0 0 m 0 0 0 0 0 0 I x I x y I x z 0 0 0 I x y I y I y z 0 0 0 I x z I y z I z , Ω x = 0 R Q 0 0 0 R 0 P 0 0 0 Q P 0 0 0 0 0 W V 0 R Q W 0 U R 0 P V U 0 Q P 0 ,
where m is the quadrotor mass, I x ,   I y ,   I z are the moments of inertia, and I x y ,   I y z ,   I x z are the products of inertia.
The gravity loads f G vector is composed of a gravity force vector which results from an assumption that the origin of the body frame is located in the quadrotor center of gravity:
f G y = F G 0 0 0 T ,
where:
F G = m g sin Θ sin Φ cos Θ cos Φ cos Θ T .  
and g is gravity acceleration.
The vector of aerodynamic loads f A consists of the fuselage’s aerodynamic force and moment vectors:
f A x , y = F A M A T ,
where
F A = q ¯ S C X cos α cos β C X sin β C X sin α cos β T ,
and
M A = q ¯ S d C L C M C N T ,
where q ˉ is the dynamic pressure of the free stream, S and d are the reference area and length respectively, C X is the fuselage drag force coefficient, and C L , C M , C N are aerodynamic rolling, pitching, and yawing moment coefficients respectively.
The drone airspeed V t o t is:
V t o t = U U W   2 + V V W 2 + W W W 2 .
where U W , V W , W W are wind velocity components in O b x b y b z b frame. The angle of attack α is defined as [48]:
α = atan W W W U U W ,
In the numerical simulation, the function atan 2 W W W ,   U U W was used to ensure the values of α from −180° up to 180°. The angle of sideslip β is [48]:
β = asin V V W V t o t .
The four independent rotors produce the load vector f R that consists of rotor force and moment vectors:
f R x , y , Ω R = F R M R T ,
The rotors were arranged in the X (Figure 4) configuration where rotors 2 and 4 rotate clockwise and rotors 1 and 3 counterclockwise.
The total force and moment produced by the rotors are a sum of the forces and moments produced by those of the individual rotors:
F R = i = 1 4 F R i = i = 1 4 Δ X R i 0 T i + Δ T i T ,
and
M R = i = 1 4 M R i = i = 1 4 r R i × F R i + 1 i 0 0 M i + Δ M i T ,  
where r R i is a position vector of the rotor in the O b x b y b z b frame, T i and M i are the thrust and torque produced by the i -th rotor ( i 1 , 2 , 3 , 4 ), Δ T i and Δ M i are the velocity correction function of the thrust and torque, and Δ X R i is the drag force produced by the i -th rotor during forward flight. The horizontal force and banking and pitching moments produced by the rotors were neglected due to their low values and the symmetry of the rotor pairs.
The rotors’ thrust and torque are described by the following equations [2,49,50]:
T i = ρ S p R p 2 k f Ω i 2 ,
and [2,49]:
M i = ρ S p R p 3 k m Ω i 2 ,
where ρ is the air density, R p is the radius of the rotor, S p = π R p 2 is the rotor’s disc area, k f and k m are the rotor’s thrust and torque coefficients respectively, and Ω i is the rotor’s angular rate.
The rotor-induced velocity and interaction between rotors has a significant impact on the produced aerodynamic loads, which affect the energy consumption. That is why the rotor thrust and torque correction factors and drag force were obtained using the FLIGHTLAB software for the rotor model, which included the aerodynamic effects from the quadrotor’s configuration. FLIGHTLAB is a world-class software used by leading rotorcraft manufacturers and R&D institutions to model, analyze, and simulate rotors and propellers. In the modeling of the quadrotor propeller, the T-Motor 18” ×6 twin-blade propeller geometric data were used. Due to the lack of manufacturer information on the blade twist distribution and airfoils, the data from a similar propeller were used. The aerodynamic characteristics were taken for the DAE51 airfoil and modeled as linear quasi-steady. The induced velocity was modeled with the Peters-He 3-state model [51]. During the validation, the blades’ root pitch angles were modified to match the modeled and manufacturer characteristics of thrust and torque for the hover. The resulting geometric parameters of the propeller blades are shown in Table 2. Results of the simulation compared with the manufacturer data are presented in Figure 5. The developed model was used to analyze the loads in flight states other than hover.
The analysis showed some differences between the characteristics of the front and rear rotors but the means were taken into account in this study. The approximated models of the corrections factors can be formed as polynomial functions (Figure 6):
Δ T = 2.698 10 3 V t o t 3 3.636 10 2 V t o t 2 2.499 10 2 V t o t ,
and
Δ M = 3.702 10 7 V t o t 5 1.043 10 5 V t o t 4 + 6.307 10 5 V t o t 3 + 8.792 10 5 V t o t 2 + 3.363 10 4 V t o t ,
and
Δ X = 1.207 10 4 V t o t 3 2.507 10 3 V t o t 2 + 4.553 10 2 V t o t ,
where V t o t is the drone airspeed given by (12).

2.3. Wind Model

The wind field model that was used in the presented study was composed of two submodels: uniform inflow and Dryden wind turbulence model [52]. The presence of wind gusts was neglected. Total wind velocities U W n , V W n , W W n in the O n x n y n z n coordinate system are:
U W n V W n W W n = U W n u V W n u W W n u + U W n t V W n t W W n t ,
where U W n u , V W n u , W W n u are the uniform wind velocities and U W n t , V W n t , W W n t are turbulence velocities.
The uniform wind speed expressed in the O n x n y n z n coordinate frame is V W t o t . This means wind speed describes the low frequency variations and is averaged over a specific time interval (for simplicity it was assumed that V W t o t is constant). The direction of the oncoming wind was defined by the angle Ψ W (clockwise when looking from above, for example 0°, wind from the north, 90°, from the east). The wind velocities of the uniform wind field in the O n x n y n z n frame are:
U W n u V W n u W W n u = V W t o t cos Ψ W V W t o t sin Ψ W 0 ,
Turbulence is a stochastic process that is difficult to model. The continuous representation of Dryden velocity spectra with positive vertical and lateral angular rates spectra (q+ r+) were used (details could be found in [53], though herein only a brief description is presented for brevity). This model is based on the pseudorandom white noise that is passed through band-limited filters. This approach is standard in modeling wind in the case of large aircrafts but Watkins and Vino [54] showed that such kind of model could be applied to small UAVs. The transfer functions that were used to calculate turbulence velocities U W n t , V W n t , W W n t are as follows [53]:
H u s = σ u 2 L u π V t o t 1 1 + L u V t o t s
H v s = σ u L v π V t o t 1 + 3 L v V s 1 + L u V t o t s 2
H w s = σ u 2 L w π V t o t 1 + 3 L w V s 1 + L u V t o t s 2
where L u , L v , L w are turbulence scale lengths and σ u , σ v , σ w are the turbulence intensities. For altitudes under 1000 feet, the scale lengths of the turbulence are calculated as follows:
L w = h
L u = L v = h 0.177 + 0.000823 h 1.2
where h is flight altitude (must be expressed in feets). The turbulence intensities (for low altitude only) are as follows:
σ w = 0.1 W 20
σ u σ w = σ v σ w = 1 0.177 + 0.000823 h 0.4
where W 20 is wind speed at altitude 6 m (expressed in knots; this parameter was set according to results from the flight trials, please see section “Results”). The probability of exceeding the high-altitude intensity was 0.01. The random noise seeds used to generate three turbulence velocities were 23,341, 23,342, and 23,343, respectively.
Next, the total wind velocities U W n , V W n , W W n given by (23) were transformed from O n x n y n z n to the body-fixed coordinate system O b x b y b z b :
U W V W W W = cos   Θ   cos   Ψ cos Θ sin Ψ sin Θ sin Φ sin Θ cos Ψ cos Φ sin Ψ sin Φ sin Θ sin Ψ + cos Φ cos Ψ sin Φ cos Θ cos Φ sin Θ cos Ψ + sin Φ sin Ψ cos Φ sin Θ sin Ψ sin Φ cos Ψ cos Φ cos Θ U W n V W n W W n .
Finally, these values could be inserted in (12)–(14).

2.4. Quadrotor Control System Model

2.4.1. Autopilot Structure

The quadrotor is an underactuated object (six degrees of freedom and four independent inputs). The quadrotor attitude, velocity and attitude are controlled by differentiating the rotors’ angular rates. The propagation of the control signals to individual rotors is described by the equation:
Ω c 1 Ω c 2 Ω c 3 Ω c 4 = 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 u 1 u 2 u 3 u 4 ,
where Ω c i is the i -th rotor-demanded value of angular rates, and u 1 , u 2 , u 3 , u 4 are control signals for climb rate, roll rate, pitch rate, and yaw rate control respectively. The change in quadrotor altitude was obtained by increasing or decreasing the rotational speed of all four propellers simultaneously.
The rotor and electric engine dynamics were modeled as the linear first order system with T e time constant:
T e d Ω i t d t + Ω i t = Ω c i t ,
Parameter T e was assumed to be 0.05 (s) (according to ref. [49]). A simple automatic flight control system was utilized to validate the quadrotor dynamic model. The control system consists of four independent control paths for each control signal [55]. Each control path (roll, pitch, yaw angles, and altitude) forms a double-cascade control system with PID (proportional-integral-derivative) control laws. The inputs to the system are position coordinates. Figure 7 shows the overall structure of the automatic flight control system.
The anti-windup clamping method was used to prevent the controllers from performance degradation. Simulink built-in tools were used to properly tune the PID settings.

2.4.2. Attitude Channel Autopilot

The control signal U 1 was obtained from the following equation:
U 1 = K P W n e W n e + K I W n e 0 t W n e d t + K D W n e d W n e d t
where W n e = W n c W n (the difference between commanded and actual vertical velocity) and K P W n e , K I W n e , K D W n e are PID settings (−3556.149, −538.572, and −112.917, respectively). The commanded vertical velocity is calculated as:
W n c = K P z n e z n e + K I z n e 0 t z n e d t + K D z n e d z n e d t
where z n e = z n c z n (difference between commanded and actual vertical position) K P z n e , K I z n e , K D z n e are PID controller constants (2.122, 0.035, and −0.387).

2.4.3. Roll Channel Autopilot

The output control signal U 2 from the roll channel was calculated as:
U 2 = K P P e P e + K I P e 0 t P e d t + K D P e d P e d t
where P e = P c P (commanded roll rate and actual roll rate) and K P P e , K I P e , K D P e are PID settings (−6.249, −0.904, and −0.219). The commanded roll rate P c was obtained from the following equation:
P c = K P Φ e Φ e + K I Φ e 0 t Φ e d t + K D Φ e d Φ e d t
where Φ e = Φ c Φ and K P Φ e , K I Φ e , K D Φ e are the PID settings (9.584, 0.798, and 0.192 respectively).

2.4.4. Pitch Channel Autopilot

Control signal U 3 was obtained as follows:
U 3 = K P Q e Q e + K I Q e 0 t Q e d t + K D Q e d Q e d t
where Q e = Q c Q and K P Q e , K I Q e , K D Q e are PID settings (6.304, 1.174, and 0.433 respectively). Commanded pitch rate P c is calculated as follows:
Q c = K P Θ e Θ e + K I Θ e 0 t Θ e d t + K D Θ e d Θ e d t
where Θ e = Θ c Θ and K P Θ e , K I Θ e , K D Θ e are PID settings (5.191, 0.228, and 0.127).

2.4.5. Yaw Channel Autopilot

The control signal U 4 for yaw autopilot is as follows:
U 4 = K P R e R e + K I R e 0 t R e d t + K D R e d R e d t
where R e = R c R (error between the commanded yaw rate and actual yaw rate) and K P R e , K I R e , K D R e are PID coefficients (112.662, 22.778, and 3.419). Commanded yaw angular rate R c was obtained as follows:
R c = K P Ψ e Ψ e + K I Ψ e 0 t Ψ e d t + K D Ψ e d Ψ e d t
where Ψ e = Ψ c Ψ and K P Ψ e , K I Ψ e , K D Ψ e are PID settings (5.288, 0.230, and −0.151 respectively).

2.4.6. Position Controller

The roll and pitch angles in the O n x n y n z n frame are calculated as follows:
Φ c * = K P U n e U n e + K I U n e 0 t U n e d t + K D U n e d U n e d t
Θ c * = K P V n e V n e + K I V n e 0 t V n e d t + K D V n e d V n e d t
where K P U n e , K I U n e , K D U n e , K P V n e , K I V n e , K D V n e (−118.004, −189.289, −17.433, 44.373, 2.774, and 2.398) are as follows:
U n e = U n c U n
V n e = V n c V n
where U n c and V n c commanded velocities in the O n x n y n z n frame (in O n x n and O n y n directions, respectively) are as follows:
U n c = K P x n e x n e + K I x n e 0 t x n e d t + K D x n e d x n e d t
V n c = K P y n e y n e + K I y n e 0 t y n e d t + K D y n e d y n e d t
where: K P x n e , K I x n e , K D x n e , K P y n e , K I y n e , K D y n e are controller settings (0.762, 0.005, 0.014, 0.664, 0.003, and 0.001). The position errors were defined as follows:
x n e = x n c x n
y n e = y n c y n
where x n c and y n c are demanded position coordinates in the O n x n y n z n frame. The roll Φ c * and pitch Θ c * angles were transformed from the O n x n y n z n coordinate system to O b x b y b z b :
Θ c Φ c 0 = cos Ψ sin Ψ 0 sin Ψ cos Ψ 0 0 0 1 Θ c * Φ c * 0 .
The “Waypoint Follower” block from the SIMULINK UAV Toolbox library was used to follow the quadrotor path. It requires a set of predefined waypoints, current drone position, and lookahead distance as inputs. This block calculates the coordinates of lookahead point in the O n x n y n z n coordinate frame and desired yaw angle. To increase the simulation speed significantly, the path follower was executed using generated C code.

2.5. Coverage Path Planning Method

The goal of the mission is to scan a given single region and take a series of photos. The drone is equipped with a stabilized, gimbaled camera so it is reasonable to assume that the camera principal axis is pointed vertically down during the flight. The search should be realized at a constant altitude. The path should totally cover the predefined zone. It was assumed that the region of interest is approximately rectangular and could be defined by the four vertices. Moreover, the environment is known, and there are no obstacles or no-fly zones inside the area of interest (this assumption is suggested, for example, in refs. [56,57]).
Many methods have been reported to calculate the most appropriate path [58,59]. Several search patterns might be used to inspect the predefined area: parallel track, creeping line search, expanding square, sector search, etc. [60,61,62]. The trajectory was defined by a set of waypoints. After visiting all waypoints the drone should return to its initial point.
At first, the quadrotor realized creeping line search and then the parallel track. In the creeping line search the search legs are perpendicular to the major axis of the search area. This method is preferred when location of the target is more probable at one end of the region of interest. In the parallel track, the search legs were parallel to the longer axis of the area. In this way, the area was scanned two times. Algorithms presented by Andersen [63] were used to calculate the waypoint locations. The obtained waypoints were implemented manually in the Mission Planner software. The example of the trajectory generator (MATLAB code) is included in “Supplementary Materials”.

2.6. Quadrotor Energy Consumption Model

The total amount of the energy E T consumed by the quadrotor might be divided into two parts: the energy spent on propulsion E R and the energy consumed by the onboard electronics (e.g., sensors, autopilot, and other circuits) E E [14,64]:
E T = E R + E E ,
The total energy consumed by the four electric motors could be estimated as follows [4,20,22,65]:
E R = t 0 t f i = 1 4 U i t I i t d t ,
where U i and I i are the voltage and current of the i -th motor, and t f is the time of flight. This model can also be rewritten in the form of mechanical parameters of the motor as follows [8,66]:
E R = t 0 t f i = 1 4 τ i t Ω i t d t ,
where τ i is the torque and Ω i is the angular rate of the i -th motor. The dynamics of the electric motor were modeled as follows:
I z p Ω ˙ i t = τ i t k m Ω i 2 t D v Ω i t ,
where I z p is the moment of inertia of the rotating parts (propeller + motor shaft) and D v is the viscous damping coefficient of the motor.
The efficiency f r , i τ i t , Ω i t of the i -th motor is included in the model to make it more reliable [8,66]:
E R = t 0 t f i = 1 4 I p e Ω ˙ i t + k m Ω i 2 t + D v Ω i t f r , i τ i t , Ω i t Ω i t d t .
The efficiency function was obtained using polynomial fitting and is described by the following formula:
f r , i τ i t , Ω i t = a 5 Ω i 5 t + a 4 Ω i 4 t + a 3 Ω i 3 t + a 2 Ω i 2 t + a 1 Ω i t + a 0 ,
where a 5 = −7.349·10–19, a 4 = 1.173·10–14, a 3 = −5.824·10–11, a 2 = 3.328·10–8, a 1 = 0.0004759, a 1 = −0.006304, and 0 Ω i t 6000 [RPM].
The torque generated by the electric motor is proportional to the electric current that flows though the motor [22]:
τ i t = K T I i t ,
where K T is the torque constant of the electric motor. The electric current consumed by the i -th motor was calculated as follows [2,67,68]:
I i t = τ i t K T = 1 K T I p e Ω ˙ i t + k m Ω i 2 t + D v Ω i t ,
Assuming that the electric power required to operate the onboard systems E s u b is constant during the flight, the energy consumed by the electronics (other than motors) is calculated using the following formula:
E E = t 0 t f E s u b d t ,
The value of E s u b was obtained in stationary laboratory experiments. The experimental setup is presented in Figure 8.
Initially, the battery was fully charged. During the test, the electric motors were switched off. Only the onboard electronic subsystems were powered and consumed energy from the battery. The energy meter was connected between the quadrotor and the battery. Data (time, voltage, current) were analyzed online on a laptop computer and logged on a memory card.
The model parameters are listed in Table 3.

2.7. Battery Model

The battery voltage during a flight was predicted using the modified Shepherd model [69,70,71,72,73]. This model requires only a few parameters that could be obtained from laboratory tests and the datasheet.
It was assumed that the battery’s internal resistance is constant during the discharge process. The Peukert effect was not present (meaning that the capacity of the battery did not change with the amplitude of the electric current). Battery ageing, temperature effects, and self-discharge phenomena were not included in the presented model. The battery memory effect was omitted.
The battery voltage during discharge changed with time and was calculated as follows [70,72]:
U b a t = E 0 K Q Q t 0 t k I d t I * K Q Q t 0 t k I d t I + A e B t 0 t k I d t R I ,
where E 0 is the open-circuit battery voltage [V] (when no external load is connected), K is the polarization constant (V/Ah), Q is the battery capacity (Ah), I * represents the low-frequency current dynamics (A), I is the current (A), A is the exponential voltage (V), B is the exponential capacity (Ah)−1, and R is internal battery resistance [Ω].
The internal resistance during discharge was obtained from the flight logs. To obtain the I * , the first-order transfer function was used for filtration. The time constant that is required for calculating the filtered current was difficult to estimate from the manufacturer data. Due to this reason, it was assumed to be 30 s. The detailed tests of the battery were also used to fit the model parameters to the existing data. The values of the abovementioned battery model parameters are presented in Table 4.
The state of the battery charge was estimated using formula [22,71]:
S O C t = 100 S O C 0 0 t I Q d t ,
where S O C 0 is the initial state of charge, and Q is the battery capacity (Ah). S O C = 100% for a fully charged and 0% for a discharged battery. Simscape Electrical library that is available in Simulink was used to model the battery behavior [74].

3. Results

3.1. Model Implementation

The developed mathematical model of the quadrotor was implemented in MATLAB/Simulink R2020b. Simulations were evaluated on a laptop computer with an Intel(R) Core(TM) i7-8750H CPU @ 2.20 GHz, 16 GB RAM, and Windows 10 operating system. The equations of motion were integrated numerically using the fourth-order Runge–Kutta (RK4) fixed-step solver. The time step size of the simulation was set to 0.001 s. Simulink Accelerator mode option was used to speed up the model execution. The inertial measurement model was also included (the details are beyond the scope of the current study).

3.2. Flight Test Methodology

Next, flight tests in the outdoor environment were performed to acquire appropriate data (Figure 9).
The test plans consisted of plenty of visual line-of-sight flights in various flight conditions, including forward flight, climbing and descending, and turns. Each flight scenario was realized at least twice to achieve data repeatability. The test results were used to fit the simulation model of the quadrotor. Finally, a study was conducted to investigate the suitability of the simulation model for energy consumption analysis. The tests of the flight along predefined trajectories were performed, and the results were compared with the simulation cases.
Autopilot settings limited the drone’s maneuverability. Roll and pitch angle saturations were set to ±15°. The maximum allowed yaw rate was ±90°/s. Ascent speed was limited to 2 m/s and descent speed to 2.4 m/s.
Before each flight, it was ensured that the drone was airworthy, and the measurement system was properly calibrated. The onboard batteries were new to intentionally remove the influence of battery ageing effects on the results. After each flight, two kinds of data logs were available: onboard logs and telemetry data (the structures of these logs are explained in refs. [75,76], respectively). The trajectory waypoints and the information about control system configuration were logged on the memory card. The data logging frequency was 10 Hz. After the experiments, the collected data were analyzed offline using UAV Log Viewer [77] and Mission Planner software [78]. The data were divided into two sets: the first one was used to fit the model parameters and the second one was used to validate the model. Next, the data were imported using a custom developed library to MATLAB and then compared with the results of the numerical simulation.

3.3. Case 1 (Short Duration Trajectory)

At first, the quadrotor model was validated for relatively simple short-duration trajectories. An example of such a triangular trajectory is presented in Figure 10. The route was defined manually by 10 waypoints. The waypoint radius was set to 2 m. The drone started vertically, moved between the waypoints, and returned to the initial location. The initial conditions for the simulation (linear velocity, angular rates, position, and attitude) were taken from the flight logs of the drone. The mean wind velocity estimated during the flight test with use of a ground-based, digital anemometer was 4.5 m/s and azimuth 220° (the direction the wind is coming from). These values were also used in the numerical flight simulation.
The linear velocities (expressed in the O n x n y n z n frame) are presented in Figure 11, the angular rates (measured in the body-fixed O b x b y b z b frame) in Figure 12, the quadrotor position in Figure 13, and the Euler angles in Figure 14.
The maximum flight attitude was 30 m (Figure 13c).
Small inconsistencies between the roll and pitch angles (Figure 14a,b) resulted from model simplifications (e.g., uniform wind was assumed). In Figure 14c the line that represents the model output was nearly indistinguishable from the line representing the real data.
The comparison of energetic parameters is presented in Figure 15.
The consumed energy increased linearly. The initial battery state of charge was 99.8% and during the flight it decreased to approximately 97%. In Figure 15c,d some disparities were observed. In reality, the internal battery resistance also changes with time. This effect could also be included in the presented model. Load generated a voltage drop of the battery. This voltage drop was predicted correctly.

3.4. Case 2 (Long Duration Trajectory)

Next, the robustness of the model against different input data was checked. In this way, one might ensure that the model parameters are tuned properly not only for a single trajectory, but the model is also able to produce reliable results for a wide range of flight scenarios. A more complicated example of the validation result is presented in Figure 16, Figure 17 and Figure 18. The quadrotor flew along the trajectory that was defined by the 47 waypoints (Figure 16) and it performed flight along a mixed parallel track and a creeping line search pattern. The mission was realized autonomously. The flight started on the ground. Next, the drone took off, climbed to the altitude of 30 m, and carried out the planed mission. Finally, it returned to the initial point.
The mean wind speed and direction measured during the flight test were 7 m/s and 225° and these values were adopted for the simulation test. The map of the flight trajectory used in both real flight and simulation is shown in Figure 16.
The comparison of the quadrotor velocity (with respect to the ground), angular rates, position, and attitude for the test flight and simulation are presented in Figure 17, Figure 18, Figure 19 and Figure 20.
The total flight time was 1271.55 [s]. Autopilot settings limited the maximum speed in the horizontal and vertical directions.
The drone started from the ground and the vertical coordinate changed to −30 m (Figure 19c). Most of the mission was realized at a constant altitude of 30 m.
The results showed the good fit of the quadrotor dynamic model. The small differences resulted from the quality of the onboard sensors used in the test flight, simulation model of the autopilot, and estimated constant wind conditions in the simulation [79].
The results of the comparison of the energy consumption, battery state of charge, voltage, and current are presented in Figure 21.
The energy consumption is nearly a linear function of flight time. At the end of the flight the SOC was 65%. The total amount of energy consumed in the real test was 537.063 × 103 Ws and in the simulation was 549.711 × 103 Ws, which proves the good fit of the simulation model. The battery state was not measured during the test flight, but its time plot seems to be correct according to the authors’ knowledge. The maximum current did not exceed 50 A. The battery voltage decreased slowly with time. The results of the simulation are close enough to the results from the flight tests.
The trajectory pattern consists of 22 legs (13 in E–W and 9 in N–S directions) and 20 turns. The statistical analysis of real energy consumption for each pattern portion was evaluated. The histograms and fitted probability distribution plots are shown in Figure 22. It was assumed that the distribution for the N–S leg and turn is normal, since for the E–W leg the Weibull distribution was used.
The two-sided 95% confidence intervals are presented in Table 5. The values of the consumed energy in each pattern type obtained in the simulation were 15.123 kWs for the N–S leg, 22.033 kWs for the E–W leg, and 4.843 kWs for turn. All results are included in corresponding confidence intervals; therefore, it can be assumed that the model of energy consumption works properly.
The root mean square errors (RMSE) were used as a measure of similarity between the model and the real flight (Table 6).
Additionally, Theil’s Inequality Coefficient (TIC) was applied to check the model reliability [80]. TIC is a measure of degree of conformance between two time series [81]:
T I C i = 1 N k = 1 N z i t k y i t k 2 1 N k = 1 N z i t k 2 + 1 N k = 1 N y i t k 2 ,             i = 1 n y
where y is a time history of data obtained from simulation, z is a vector of samples from the experiment, N is the number of data points, n y is the number of outputs, and t k is the k -th discrete data point. In this way, Equation (63) provides n y separate TIC coefficients. Values of TIC remained between 0 and 1. TIC = 0 means that a strong correlation between data and model existed. On the other hand, TIC = 1 when the model did not represent the reality. For practical purposes, TIC < 0.25 is preferred [80]. As mentioned before, TIC was calculated for several (12) flight parameters separately (Table 6).
The largest disparities were obtained for quadrotor roll and pitch angular rates (TIC > 0.7). The model predictions for vertical position z n were very close to the experimental results.

3.5. Monte-Carlo Validation of the Energy Consumption

The model parameters were obtained with some uncertainty. Next, to investigate the model sensitivity on the abovementioned parameter’s uncertainties, the Monte-Carlo simulation was evaluated. Quadrotor mass, moments of inertia, initial conditions, wind speed, and azimuth were taken into consideration. Each parameter was calculated as a sum of two components: nominal value μ and normally distributed pseudorandom disturbance with zero mean and standard deviation σ [82,83,84]. The input data to the simulation are presented in Table 7.
The values of standard deviations were obtained by laboratory experiments. To generate pseudorandom input data, the Marsenne–Twister algorithm [85] was used (initial seed was set to 0). For each case, 100 runs were evaluated. The average time of a single simulation run for case 1 was 14.9 s and for case 2 was approximately 179.8 s.
The obtained trajectories are presented in Figure 23.
The mean, minimum, and maximum values of consumed energy are presented in Table 8.
At first, the influence of all disturbances was tested. In Figure 24a,b the kernel density estimators of the total consumed energy are presented.
Next, the same scenarios were repeated but only with disturbances in wind conditions V W t o t and Ψ W (Figure 24c,d). Later, only the disturbances in mass m and moments of inertia I x x , I y y , I z z were tested (Figure 24e,f). Wind speed and direction are some of the major factors that could affect the drone dynamics and energy consumption. Additionally, the energy consumed in the real flight was marked by a vertical line as a reference.
The uncertainty of parameters could influence the prediction of total consumed energy. For case 1 (Figure 24a) the minimum consumed energy was between 46 × 103 Ws and 51 × 103 Ws. The model underestimated the total energy consumption. The opposite situation was obtained in case 2 (Figure 24b). The consumed energy fell into the range from 520 × 103 Ws to 570 × 103 Ws and the model overestimated the results.
Considering only wind disturbances resulted in quantitatively different results. The curve presented in Figure 24c is tall and narrow. The predicted amount of energy was smaller than in reality. The curves presented in Figure 24e,f are very similar for all disturbances (Figure 24a,b). The obtained results qualitatively agree well with the analysis presented in ref. [86].

3.6. Wind Influence on the Energy Consumption: Practical Example

Finally, a practical example of the model use is presented. To investigate the effect of wind on the amount of consumed energy, four simulation scenarios were considered: headwind and tailwind combined with two different wind velocities. Initially, the quadrotor was at hover at altitude 40 m above the ground. Next, the quadrotor moved from point (0, 0, −40) to (200, 0, −40) along a straight line at a constant altitude. In the first and second scenarios (headwind), the mean wind azimuth was 0°, but in the third and fourth scenarios (tailwind) this angle was set to 180°. The standard deviation of wind azimuth was the same for both analyzed cases and equal to 15°. Wind speed was set to 4.5 m/s and 9 m/s, respectively. All other disturbances were omitted (standard deviations for them were set to 0). For each scenario 100 Monte-Carlo runs were evaluated (total number of runs was 400). The obtained results are presented in Figure 25.
For headwind, the right-skewed histograms were obtained. When the wind azimuth was disturbed then the drone had to roll to eliminate the lateral drift effect and keep a proper, straight flight direction. The consumed energy in the first scenario was higher than in the second case. On the other hand, for tailwind the left-skewed plots were generated. For headwind and wind speed 4.5 m/s (Figure 25a) the consumed energy was lower than for 9 m/s (Figure 25c). For tailwind and wind speed 4.5 m/s (Figure 25b) the total amount of energy was higher than for 9 m/s (Figure 25d). This simple example illustrates how the model might be used in energy analysis.

4. Conclusions

The methodology of the quadrotor dynamics model design was presented. The nonlinear quadrotor dynamic model was validated using real flight data and laboratory tests. Several contributions of this paper might be mentioned.
First, the rotor thrust, torque, and drag force correction factors were used to model the quadrotor dynamics correctly in forward flight. FLIGHTLAB software was used to calculate abovementioned coefficients. In this way, the model might predict the drone behavior for a wide range of flight conditions.
Second, an analysis of energy consumption estimation during a long duration outdoor flight was performed as well. This partially fills the gap in the literature, because such validation results for several minutes of flight are not commonly available.
Third, the simulation model was tuned using the data of real quadrotor obtained in laboratory experiments and flight tests.
Fourth, the energy consumed by the onboard electronics was taken into account. In this way, the presented study also extends the energy consumption models presented by Yacef et al. [14,21,22,55] and Morbidi [4,20]. The obtained results confirm conclusions presented by Cabreira et al. [87].
Fifth, the results of Monte-Carlo simulation indicated that model parameter uncertainties might influence the predicted amount of consumed energy. This method is not commonly used in the existing literature for validation quadrotor energy models. However, it allows understanding the impact of parameters on the model output. The uncertainty in the predicted amount of energy increases with the flight duration.
Sixth, the raw data from the flight trials are included as an Appendix A to the presented study. This dataset might be useful for other researchers.
The time domain analysis of the simulation and flight tests results indicated its good capability for estimating quadrotor energy consumption and usefulness for mission plan optimization. Scientific significance of this work is that this model might be adopted for other quadrotors and used for energy estimating purposes.
Further research might concentrate on flight tests in various flight conditions to gather more data and increase model reliability. Wind tunnel tests of the fuselage and isolated propellers could be performed to obtain aerodynamic coefficients that are used in the simulation model. A detailed model of rotors might be also used. Battery ageing effects might be included into the simulation to make it even more realistic. The developed model could be used for optimal coverage path planning, taking consumed energy as a performance metric.

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/en15197136/s1, flight data that have been used to produce Figure 10, Figure 11, Figure 12, Figure 13, Figure 14, Figure 15, Figure 16, Figure 17, Figure 18, Figure 19, Figure 20 and Figure 21.

Author Contributions

Conceptualization, R.G., P.B. and M.Ż.; methodology, M.Ż.; software, M.J.; validation, M.J., M.Ż.; formal analysis, M.Ż.; investigation, M.J., P.B.; resources, R.G.; data curation, M.J.; writing—original draft preparation, M.J. and M.Ż.; writing—review and editing, M.J.; visualization, M.J. and M.Ż.; supervision, M.Ż.; project administration, R.G.; funding acquisition, R.G. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by National Center for Research and Development (NCBiR), grant number 4/1.1.1/2020 “Renewable Energy Sources (RES) in transport”.

Data Availability Statement

Flight data are available as Supplementary Materials. Simulation information will be available on request to the corresponding author’s email with appropriate justification.

Acknowledgments

Special thanks to Marcin Kasprzyk for conducting the flight trials.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Nomenclature

The following abbreviations and symbols are used in this manuscript:
Latin symbols
A exponential voltage, (V)
A aircraft inertia matrix
B exponential capacity, (Ah)–1
B gyroscopic matrix
c propeller chord, (mm)
C L , C M , C N rolling, pitching and yawing moment coefficients, (–)
C X , C Y , C Z drag, side and lift force coefficients, (–)
d reference linear dimension, (m)
D p diameter of the rotor, (m)
D v viscous damping coefficient of the motor, (Nm·s/rad)
E 0 open circuit battery voltage, (V)
E E energy consumed by the onboard electronics, (J)
E T total amount of energy consumed by the quadrotor, (J)
E R energy spent on propulsion, (J)
E s u b electric power required by the onboard subsystems, (J/s)
f r , i τ i t , Ω i t efficiency of the i -th electric motor, (–)
f A vector of aerodynamic loads,
f G vector of gravity loads,
f R vector of propulsion loads,
F G vector of gravity forces, (N)
F A vector of aerodynamic forces, (N)
F R vector of propulsion forces, (N)
g gravity acceleration, (m/s2)
h flight altitude, (m)
H u s , H v s , H w s transfer functions for Dryden wind model,
I i current of the i -th electric motor, (A)
I x , I y ,   I z quadrotor moments of inertia, (kg·m2)
I x y , I y z ,   I x z quadrotor products of inertia, (kg·m2)
I z p moment of inertia of the of the rotating parts (propeller + motor shaft), (kg·m2)
k f rotor thrust coefficient, (N/RPM2)
k m rotor torque coefficient, (Nm/RPM2)
K battery polarization constant, (V/Ah)
K P W n e , K I W n e , K D W n e vertical speed controller settings,
K P z n e , K I z n e , K D z n e altitude controller settings,
K P P e , K I P e , K D P e roll rate controller settings,
K P Φ e , K I Φ e , K D Φ e roll angle controller settings,
K P Q e , K I Q e , K D Q e pitch rate controller settings,
K P Θ e , K I Θ e , K D Θ e pitch angle controller settings,
K P R e , K I R e , K D R e yaw rate controller settings,
K P Ψ e , K I Ψ e , K D Ψ e yaw angle controller settings,
K P U n e , K I U n e , K D U n e U n speed controller settings,
K P V n e , K I V n e , K D V n e V n speed controller settings,
K P x n e , K I x n e , K D x n e x n position controller settings,
K P y n e , K I y n e , K D y n e y n position controller settings,
L u , L v , L w turbulence scale lengths, (m)
m drone mass, (kg)
M A vector of aerodynamic moments, (Nm)
M i torque produced by the i -th rotor, (Nm)
M R vector of moments generated by propellers, (Nm)
N number of discrete sample points,
P ,   Q ,   R quadrotor angular velocities in O b x b y b z b reference frame, (rad/s)
q ˉ dynamic pressure, (kg/(m·s2))
Q battery capacity, (Ah)
r radial distance from the propeller axis of rotation, (mm)
r n = x n y n z n T vector of the quadrotor position in the O n x n y n z n frame, (m)
r R i position vector of the i -th rotor in the O b x b y b z b frame, (m)
R battery internal resistance, (Ω)
R p radius of the rotor, (m)
S reference area, (m2)
S p rotor disc area, (m2)
t time, (s)
t 0 initial time, (s)
t f final time, (s)
t k k -th discrete time step, (s)
T e time constant, (s)
T i thrust produced by the i -th rotor, (N)
T matrix of linear velocities and angular rates
T V velocity transformation matrix
T Ω angle transformation matrix
u 1 , u 2 , u 3 , u 4 control signals for climb rate, roll rate, pitch rate and yaw rate control respectively, (rad/s)
U ,   V ,   W components of linear velocity in O b x b y b z b reference frame, (m/s)
U n ,   V n ,   W n components of linear velocity in O n x n y n z n reference frame, (m/s)
U W ,   V W ,   W W components of wind velocity in O b x b y b z b reference frame, (m/s)
U W n u ,   V W n u ,   W W n u components of uniform wind velocity in O n x n y n z n reference frame, (m/s)
U W n t ,   V W n t ,   W W n t components of turbulence velocity in O n x n y n z n reference frame, (m/s)
U W n ,   V W n ,   W W n components of wind velocity in O n x n y n z n reference frame, (m/s)
U b a t battery voltage, (V)
U i voltage of the i -th motor, (V)
V t o t airspeed, (m/s)
v = [ U , V , W ] T linear velocity vector in O b x b y b z b reference frame, (m/s)
V t o t drone airspeed, (m/s)
V W t o t total wind speed, (m/s)
W 20 wind speed at altitude 6 m, (kts)
x state vector
x n ,   y n , z n coordinates of aircraft position in O n x n y n z n reference frame, (m)
y model output vector
y aircraft position and attitude vector
z measurement vector
Greek symbols
α angle of attack, (rad)
β angle of sideslip, (rad)
γ rotor blade twist, (rad)
Δ c chord center shift distribution along the rotor span, (mm)
Δ T i velocity correction function of the thrust, (N)
Δ M i velocity correction function of the torque, (Nm)
μ nominal value of the parameter
ρ air density, (kg/m3)
σ standard deviation of the pseudorandom disturbance
σ u , σ v , σ w turbulence intensities,
τ i torque generated by the i -th motor, (Nm)
Φ ,   Θ ,   Ψ quadrotor attitude angles, (rad)
Ψ W wind direction (direction of oncoming flow), (rad)
ω = [ P Q R ] T vector of angular velocity, (rad/s)
Ω i angular rate of i -th motor, (rad/s)
Ω c i i -th rotor demanded value of angular rates, (rad/s)
Ω R the vector of the angular velocity of the rotors, (rad/s)
Abbreviations
BLDCBrushless Direct Current
CADComputer Aided-design
GPSGlobal Positioning System
GNSSGlobal Navigation Satellite System
IMUInertial Measurement Unit
UAVUnmanned Aerial Vehicle
PIDproportional-integral-derivative
RMSERoot Mean Square Errors
RPMrevolutions per minute
S O C battery State Of Charge, (%)
S O C 0 initial battery State Of Charge, (%)
TICTheil’s Inequality Coefficient

Appendix A. Flight Data and MATLAB Scripts

The folder “Supplementary Materials” includes the two subfolders “Case 1” and “Case 2” with flight data logs in “TLOG” and “BIN” formats. These data could be viewed using UAV Log Viewer online service that is available at:
https://plot.ardupilot.org/#/ (accessed on 13 September 2022).
and Mission Planner software:
Additionally, files with “*.log” extension includes data about quadrotor configuration. Files with “*.waypoints” extension include the coordinates of trajectory waypoints.
In “Trajectory” there are MATLAB scripts that could be used to reproduce similar trajectory as presented in “Case 2” (the user could change some parameters).

References

  1. Bibik, P.; Narkiewicz, J.; Zasuwa, M.; Żugaj, M. Modeling of Quadrotor Dynamics for Research and Training Simulator. In Proceedings of the 39th European Rotorcraft Forum, Moscow, Russia, 3–6 September 2013. [Google Scholar]
  2. Lu, H.; Chen, K.; Zhai, X.B.; Chen, B.; Zhao, Y. Tradeoff between Duration and Energy Optimization for Speed Control of Quadrotor Unmanned Aerial Vehicle. In Proceedings of the ISPCE-CN 2018—IEEE International Symposium on Product Compliance Engineering—Asia, Shenzhen, China, 5–7 December 2018. [Google Scholar]
  3. Aleksandrov, D.; Penkov, I. Energy Consumption of Mini UAV Helicopters with Different Number of Rotors. In Proceedings of the 11th International Symposium Topical Problems in the Field of Electrical and Power Engineering, Pärnu, Estonia, 16–21 January 2012; pp. 259–262. [Google Scholar]
  4. Morbidi, F.; Pisarski, D. Practical and Accurate Generation of Energy-Optimal Trajectories for a Planar Quadrotor. In Proceedings of the IEEE International Conference on Robotics and Automation, Xi’an, China, 30 May 2021–5 June2021; pp. 355–361. [Google Scholar]
  5. Bibik, P.; Narkiewicz, J.; Zasuwa, M.; Żugaj, M.; Górski, T.; Komorniczak, W. Development of an Unmanned Quadrotor: System and Simulator. Available online: https://5dok.net/document/ozlrv66z-development-of-an-unmanned-quadrotor-system-and-simulator.html (accessed on 18 August 2022).
  6. Gandolfo, D.C.; Salinas, L.R.; Brandao, A.; Toibero, J.M. Stable Path-Following Control for a Quadrotor Helicopter Considering Energy Consumption. IEEE Trans. Control Syst. Technol. 2017, 25, 1423–1430. [Google Scholar] [CrossRef]
  7. Roberts, J.F.; Zufferey, J.C.; Floreano, D. Energy Management for Indoor Hovering Robots. In Proceedings of the 2008 IEEE/RSJ International Conference on Intelligent Robots and Systems, IROS, Nice, France, 22–26 September 2008; pp. 1242–1247. [Google Scholar]
  8. Chan, C.W.; Kam, T.Y. A Procedure for Power Consumption Estimation of Multi-Rotor Unmanned Aerial Vehicle. J. Phys. Conf. Ser. 2020, 1509, 012015. [Google Scholar] [CrossRef]
  9. Penkov, I.; Aleksandrov, D. Analysis and Study of the Influence of the Geometrical Parameters of Mini Unmanned Quad-Rotor Helicopters to Optimise Energy Saving. Int. J. Automot. Mech. Eng. 2017, 14, 4730–4746. [Google Scholar] [CrossRef]
  10. Henninger, H.C.; Von Ellenrieder, K.D.; Licht, S.C. Energy-Minimal Target Retrieval for Quadrotor UAVs: Trajectory Generation and Tracking. In Proceedings of the 2020 28th Mediterranean Conference on Control and Automation, MED 2020, Saint-Raphaël, France, 15–18 September 2020; pp. 727–732. [Google Scholar]
  11. Korneyev, A.; Gorobetz, M.; Alps, I.; Ribickis, L. Adaptive Traction Drive Control Algorithm for Electrical Energy Consumption Minimisation of Autonomous Unmanned Aerial Vehicle. Electr. Control Commun. Eng. 2019, 15, 62–70. [Google Scholar] [CrossRef]
  12. Thu, A.; Lupin, S.; Oo, T.M.; Khaing, M.T. Comparing a Quadrotor Energy Consumption for Different Flight Trajectories in Windy Conditions. In Proceedings of the 2021 IEEE Conference of Russian Young Researchers in Electrical and Electronic Engineering, ElConRus 2021, St. Petersburg, Moscow, Russia, 26–29 January 2021; pp. 2064–2066. [Google Scholar]
  13. Shivgan, R.; Dong, Z. Energy-Efficient Drone Coverage Path Planning Using Genetic Algorithm. In Proceedings of the IEEE International Conference on High Performance Switching and Routing, HPSR, Newark, NJ, USA, 11–14 May 2020. [Google Scholar]
  14. Yacef, F.; Rizoug, N.; Degaa, L.; Hamerlain, M. Energy-Efficiency Path Planning for Quadrotor UAV under Wind Conditions. In Proceedings of the 7th International Conference on Control, Decision and Information Technologies, CoDIT 2020, Prague, Czech Republic, 29 June–2 July 2020; pp. 1133–1138. [Google Scholar]
  15. Dietrich, T.; Krug, S.; Zimmermann, A. An Empirical Study on Generic Multicopter Energy Consumption Profiles. In Proceedings of the 11th Annual IEEE International Systems Conference, SysCon 2017—Proceedings, Montreal, QC, Canada, 24–27 April 2017; pp. 406–411. [Google Scholar]
  16. Chen, Y.; Baek, D.; Bocca, A.; Macii, A.; Macii, E.; Poncino, M. A Case for a Battery-Aware Model of Drone Energy Consumption. In Proceedings of the INTELEC, International Telecommunications Energy Conference (Proceedings), Turin, Italy, 7–11 October 2019. [Google Scholar]
  17. Lee, C.; Son, J.J.; Lee, H.; Han, S. Energy Consumption Analysis of Downward-Tethered Quadcopter. In Proceedings of the 21st International Conference on Control, Automation and Systems (ICCAS 2021), ICROS: Ramada Plaza Hotel, Jeju, Korea, 12–15 October 2021; pp. 1212–1215. [Google Scholar]
  18. Zhang, J.; Campbell, J.F.; Sweeney, D.C.; Hupman, A.C. Energy Consumption Models for Delivery Drones: A Comparison and Assessment. Transp. Res. Part D Transp. Environ. 2021, 90, 102668. [Google Scholar] [CrossRef]
  19. Beigi, P.; Rajabi, M.S.; Aghakhani, S. An Overview of Drone Energy Consumption Factors and Models. arXiv 2022, arXiv:2206.10775. [Google Scholar]
  20. Morbidi, F.; Cano, R.; Lara, D.; Morbidi, F.; Cano, R.; Lara, D.; Generation, M.P.; Morbidi, F.; Cano, R.; Lara, D. Minimum-Energy Path Generation for a Quadrotor UAV. In Proceedings of the International Conference on Robotics and Automation, Stockholm, Sweden, 16–21 May 2016; pp. 2–8. [Google Scholar]
  21. Yacef, F.; Rizoug, N.; Degaa, L.; Bouhali, O.; Hamerlain, M. Trajectory Optimisation for a Quadrotor Helicopter Considering Energy Consumption. In Proceedings of the 2017 4th International Conference on Control, Decision and Information Technologies, CoDIT 2017, Barcelona, Spain, 5–7 April 2017; pp. 1030–1035. [Google Scholar]
  22. Fouad, Y.; Rizoug, N.; Bouhali, O.; Hamerlain, M. Optimization of Energy Consumption for Quadrotor UAV. In Proceedings of the International Micro Air Vehicle Conference and Flight Competition (IMAV) 2017, Toulouse, France, 18–21 September 2017; pp. 215–222. [Google Scholar]
  23. Jee, S.; Cho, H. Comparing Energy Consumption Following Flight Pattern for Quadrotor. J. IKEEE 2018, 22, 747–753. [Google Scholar] [CrossRef]
  24. Li, B.; Li, Q.; Zeng, Y.; Rong, Y.; Zhang, R. 3D Trajectory Optimization for Energy-Efficient UAV Communication: A Control Design Perspective. IEEE Trans. Wirel. Commun. 2022, 21, 4579–4593. [Google Scholar] [CrossRef]
  25. Wang, Y.; Wang, Y.; Ren, B. Energy Saving Quadrotor Control for Field Inspections. IEEE Trans. Syst. Man, Cybern. Syst. 2022, 52, 1768–1777. [Google Scholar] [CrossRef]
  26. Salameh, I.M.; Ammar, E.M.; Tutunji, T.A. Identification of Quadcopter Hovering Using Experimental Data. In Proceedings of the 2015 IEEE Jordan Conference on Applied Electrical Engineering and Computing Technologies, AEECT 2015, Amman, Jordan, 3–5 November 2015. [Google Scholar]
  27. Gao, N.; Zeng, Y.; Wang, J.; Wu, D.; Zhang, C.; Song, Q.; Qian, J.; Jin, S. Energy Model for UAV Communications: Experimental Validation and Model Generalization. China Commun. 2021, 18, 253–264. [Google Scholar] [CrossRef]
  28. Aguilar-Lopez, J.M.; Garcia, R.A.; Bordons, C.; Camacho, E.F. Development of the Energy Consumption Model of a Quadrotor Using Voltage Data from Experimental Flights. In Proceedings of the 2022 IEEE 17th International Conference on Control & Automation (ICCA), Naples, Italy, 27–30 June 2022; pp. 432–437. [Google Scholar] [CrossRef]
  29. Alyassi, R.; Khonji, M.; Karapetyan, A.; Chau, S.C.-K.; Elbassioni, K.; Tseng, C.-M. Autonomous Recharging and Flight Mission Planning for Battery-Operated Autonomous Drones. IEEE Trans. Autom. Sci. Eng. 2022, 1–13. [Google Scholar] [CrossRef]
  30. Abeywickrama, H.V.; Jayawickrama, B.A.; He, Y.; Dutkiewicz, E. Empirical Power Consumption Model for UAVs. In Proceedings of the IEEE Vehicular Technology Conference, Chicago, IL, USA, 27–30 August 2018; pp. 1–5. [Google Scholar]
  31. Abeywickrama, H.V.; Jayawickrama, B.A.; He, Y.; Dutkiewicz, E. Comprehensive Energy Consumption Model for Unmanned Aerial Vehicles, Based on Empirical Studies of Battery Performance. IEEE Access 2018, 6, 58383–58394. [Google Scholar] [CrossRef]
  32. Steup, C.; Parlow, S.; Mai, S.; Mostaghim, S. Generic Component-Based Mission-Centric Energy Model for Micro-Scale Unmanned Aerial Vehicles. Drones 2020, 4, 63. [Google Scholar] [CrossRef]
  33. Kreciglowa, N.; Karydis, K.; Kumar, V. Energy Efficiency of Trajectory Generation Methods for Stop-and-Go Aerial Robot Navigation. In Proceedings of the 2017 International Conference on Unmanned Aircraft Systems, ICUAS 2017, Miami, FL, USA, 13–16 June 2017; pp. 656–662. [Google Scholar]
  34. Rodrigues, T.A.; Patrikar, J.; Wagner, B.; Scherer, S.; Samaras, C. Development of an Energy Model for Quadcopter Package Delivery Drones. Available online: https://www.microstrain.com/sites/default/files/energy_model_for_quadcopter_drones_lord-compressed.pdf (accessed on 13 August 2022).
  35. Wu, F.; Yang, D.; Xiao, L.; Cuthbert, L. Energy Consumption and Completion Time Tradeoff in Rotary-Wing UAV Enabled WPCN. IEEE Access 2019, 7, 79617–79635. [Google Scholar] [CrossRef]
  36. Zeng, Y.; Xu, J.; Zhang, R. Energy Minimization for Wireless Communication with Rotary-Wing UAV. IEEE Trans. Wirel. Commun. 2019, 18, 2329–2345. [Google Scholar] [CrossRef]
  37. Pradeep, P.; Park, S.G.; Wei, P. Trajectory Optimization of Multirotor Agricultural UAVs. In Proceedings of the 2018 IEEE Aerospace Conference, Big Sky, MT, USA, 3–10 March 2018; pp. 1–7. [Google Scholar]
  38. Sekander, S.; Tabassum, H.; Hossain, E. On the Performance of Renewable Energy-Powered UAV-Assisted Wireless Communications. arXiv 2019, arXiv:1907.07158. [Google Scholar]
  39. Bangura, M.; Mahony, R. Nonlinear Dynamic Modeling for High Performance Control of a Quadrotor. In Proceedings of the Australasian Conference on Robotics and Automation 2012, Wellington, New Zealand, 3–5 December 2012; pp. 1–10. [Google Scholar]
  40. Huang, H.; Hoffmann, G.M.; Waslander, S.L.; Tomlin, C.J. Aerodynamics and Control of Autonomous Quadrotor Helicopters in Aggressive Maneuvering. In Proceedings of the IEEE International Conference on Robotics and Automation, Kobe, Japan, 12–17 May 2009; pp. 3277–3282. [Google Scholar]
  41. Liu, Z.; Sengupta, R.; Kurzhanskiy, A. A Power Consumption Model for Multi-Rotor Small Unmanned Aircraft Systems. In Proceedings of the 2017 International Conference on Unmanned Aircraft Systems, ICUAS 2017, Miami, FL, USA, 13–16 June 2017; pp. 310–315. [Google Scholar]
  42. Stolaroff, J.K.; Samaras, C.; O’Neill, E.R.; Lubers, A.; Mitchell, A.S.; Ceperley, D. Energy Use and Life Cycle Greenhouse Gas Emissions of Drones for Commercial Package Delivery. Nat. Commun. 2018, 9, 1–13. [Google Scholar] [CrossRef]
  43. Ware, J.; Roy, N. An Analysis of Wind Field Estimation and Exploitation for Quadrotor Flight in the Urban Canopy Layer. In Proceedings of the 2016 IEEE International Conference on Robotics and Automation (ICRA), Stockholm, Sweden, 16–21 May 2016; pp. 1507–1514. [Google Scholar] [CrossRef]
  44. Rodrigues, T.A.; Patrikar, J.; Choudhry, A.; Feldgoise, J.; Arcot, V.; Gahlaut, A.; Lau, S.; Moon, B.; Wagner, B.; Matthews, H.S.; et al. In-Flight Positional and Energy Use Data Set of a DJI Matrice 100 Quadcopter for Small Package Delivery. Sci. Data 2021, 8, 6–13. [Google Scholar] [CrossRef]
  45. The Cube Orange. Available online: https://docs.px4.io/master/en/flight_controller/cubepilot_cube_orange.html (accessed on 18 August 2022).
  46. Benic, Z.; Piljek, P.; Kotarski, D. Mathematical Modelling of Unmanned Aerial Vehicles with Four Rotors. Interdiscip. Descr. Complex Syst. 2016, 14, 88–100. [Google Scholar] [CrossRef]
  47. Żugaj, M.; Bibik, P.; Jacewicz, M. UAV Aircraft Model for Control System Failures Analysis. J. Theor. Appl. Mech. 2016, 54, 1405–1415. [Google Scholar] [CrossRef] [Green Version]
  48. Zipfel, P. Modeling and Simulation of Aerospace Vehicle Dynamics; American Institute of Aeronautics and Astronautics, Inc.: Reston, VA, USA, 2000. [Google Scholar]
  49. Pounds, P.; Mahony, R.; Corke, P. Modelling and Control of a Large Quadrotor Robot. Control Eng. Pract. 2010, 18, 691–699. [Google Scholar] [CrossRef]
  50. Bezzo, N.; Mohta, K.; Nowzari, C.; Lee, I.; Kumar, V.; Pappas, G. Online Planning for Energy-Efficient and Disturbance-Aware UAV Operations. In Proceedings of the IEEE International Conference on Intelligent Robots and Systems, Daejeon, Korea, 9–14 October 2016; pp. 5027–5033. [Google Scholar]
  51. Peters, D.A.; He, C.J. Finite State Induced Flow Models. II—Three-Dimensional Rotor Disk. J. Aircr. 2012, 32, 323–333. [Google Scholar] [CrossRef]
  52. MIL-F-8785C Military Specification. Flying Qualities of Piloted Airplanes. Available online: http://everyspec.com/MIL-SPECS/MIL-SPECS-MIL-F/MIL-F-8785C_5295/ (accessed on 11 September 2022).
  53. Dryden Wind Turbulence Model (Continuous). Available online: https://www.mathworks.com/help/aeroblks/drydenwindturbulencemodelcontinuous.html (accessed on 13 September 2022).
  54. Watkins, S.; Vino, G. The Turbulent Wind Environment of Birds, Insects and MAVs. In Proceedings of the 15th Australasian Fluid Mechanics Conference, Sydney, Australia, 13–17 December 2004. [Google Scholar]
  55. Yacef, F.; Bouhali, O.; Hamerlain, M.; Rizoug, N. Observer-Based Adaptive Fuzzy Backstepping Tracking Control of Quadrotor Unmanned Aerial Vehicle Powered by Li-Ion Battery. J. Intell. Robot. Syst. Theory Appl. 2016, 84, 179–197. [Google Scholar] [CrossRef]
  56. Tang, G.; Tang, C.; Zhou, H.; Claramunt, C.; Men, S. R-Dfs: A Coverage Path Planning Approach Based on Region Optimal Decomposition. Remote Sens. 2021, 13, 1525. [Google Scholar] [CrossRef]
  57. Di Franco, C.; Buttazzo, G. Coverage Path Planning for UAVs Photogrammetry with Energy and Resolution Constraints. J. Intell. Robot. Syst. Theory Appl. 2016, 83, 445–462. [Google Scholar] [CrossRef]
  58. Majeed, A.; Hwang, S.O. A Multi-Objective Coverage Path Planning Algorithm for Uavs to Cover Spatially Distributed Regions in Urban Environments. Aerospace 2021, 8, 343. [Google Scholar] [CrossRef]
  59. Fevgas, G.; Lagkas, T.; Argyriou, V.; Sarigiannidis, P. Coverage Path Planning Methods Focusing on Energy Efficient and Cooperative Strategies for Unmanned Aerial Vehicles. Sensors 2022, 22, 1235. [Google Scholar] [CrossRef] [PubMed]
  60. Otote, D.A.; Li, B.; Ai, B.; Gao, S.; Xu, J.; Chen, X.; Lv, G. A Decision-Making Algorithm for Maritime Search and Rescue Plan. Sustainability 2019, 11, 2084. [Google Scholar] [CrossRef]
  61. Cabreira, T.M.; Brisolara, L.B.; Ferreira Paulo, R. Survey on Coverage Path Planning with Unmanned Aerial Vehicles. Drones 2019, 3, 4. [Google Scholar] [CrossRef]
  62. Choutri, K.; Lagha, M.; Dala, L. A Fully Autonomous Search and Rescue System Using Quadrotor UAV. Int. J. Comput. Digit. Syst. 2021, 10, 403–414. [Google Scholar] [CrossRef]
  63. Andersen, H. Path Planning for Search and Rescue Mission Using Multicopters. Master’s Thesis, Norwegian University of Science and Technology, Trondheim, Norway, 2014. [Google Scholar]
  64. Jaafar, W.; Yanikomeroglu, H. Dynamics of Quadrotor UAVs for Aerial Networks: An Energy Perspective. arXiv 2019, arXiv:1905.06703. [Google Scholar]
  65. Pradeep, P.; Wei, P. Energy-Efficient Arrival with RTA Constraint for Multirotor EVToL in Urban Air Mobility. J. Aerosp. Inf. Syst. 2019, 16, 263–277. [Google Scholar] [CrossRef]
  66. Prasetia, A.S.; Wai, R.J.; Wen, Y.L.; Wang, Y.K. Mission-Based Energy Consumption Prediction of Multirotor Uav. IEEE Access 2019, 7, 33055–33063. [Google Scholar] [CrossRef]
  67. Li, M.; Jia, G.; Gong, S.; Guo, R. Energy Consumption Model of BLDC Quadrotor UAVs for Mobile Communication Trajectory Planning. TechRxiv 2022, arXiv:19181228.v1. [Google Scholar]
  68. Aoun, C.; Daher, N.; Shammas, E. An Energy Optimal Path-Planning Scheme for Quadcopters in Forests. In Proceedings of the IEEE Conference on Decision and Control, Nice, France, 11–13 December 2019; pp. 8323–8328. [Google Scholar]
  69. Tremblay, O.; Dessaint, L.A. Experimental Validation of a Battery Dynamic Model for EV Applications. World Electr. Veh. J. 2009, 2, 289–298. [Google Scholar] [CrossRef]
  70. Mousavi, G.S.M.; Nikdel, M. Various Battery Models for Various Simulation Studies and Applications. Renew. Sustain. Energy Rev. 2014, 32, 477–485. [Google Scholar] [CrossRef]
  71. Azam, S.M. Battery Identification, Prediction and Modelling. Master’s Thesis, Colorado State University, Fort Collins, CO, USA, 2018. [Google Scholar]
  72. Raszmann, E.; Baker, K.; Shi, Y.; Christensen, D. Modeling Stationary Lithium-Ion Batteries for Optimization and Predictive Control. In Proceedings of the 2017 IEEE Power and Energy Conference at Illinois (PECI), Champaign, IL, USA, 23–24 February 2017. [Google Scholar] [CrossRef]
  73. Hemi, H.; M’sirdi, N.K.; Naamane, A. A New Proposed Shepherd Model of a Li-Ion Open Circuit Battery Based on Data Fitting. In Proceedings of the International Conference on Integrated Modeling and Analysis in Applied Control and Automation, Lisbon, Portugal, 18–20 September 2019; pp. 83–92. [Google Scholar]
  74. Simulink Generic Battery Model. Available online: https://www.mathworks.com/help/physmod/sps/powersys/ref/battery.html (accessed on 16 August 2022).
  75. Onboard Message Log Messages. Available online: https://ardupilot.org/copter/docs/logmessages.html (accessed on 18 August 2022).
  76. MAVLink Messages. Available online: https://mavlink.io/en/messages/ardupilotmega.html#messages (accessed on 13 August 2022).
  77. UAV Log Viewer. Available online: https://plot.ardupilot.org/#/ (accessed on 13 August 2022).
  78. Mission Planner Home. Available online: https://ardupilot.org/planner/ (accessed on 16 August 2022).
  79. Jimenez, P.; Lichota, P.; Agudelo, D.; Rogowski, K. Experimental Validation of Total Energy Control System for UAVs. Energies 2020, 13, 14. [Google Scholar] [CrossRef]
  80. Dorobantu, A.; Seiler, P.J.; Balas, G.J. Validating Uncertain Aircraft Simulation Models Using Flight Test Data. In Proceedings of the AIAA Atmospheric Flight Mechanics (AFM) Conference, Boston, MA, USA, 19–22 August 2013. [Google Scholar]
  81. Jategaonkar, R.V. Flight Vehicle System Identification: A Time-Domain Methodology, 2nd ed.; American Institute of Aeronautics and Astronautics: Reston, VA, USA, 2015. [Google Scholar]
  82. Jacewicz, M.; Głębocki, R.; Ożóg, R. Monte-Carlo Based Lateral Thruster Parameters Optimization for 122 mm Rocket. In Proceedings of the Automation 2020: Towards Industry of the Future. AUTOMATION 2020, Warsaw, Poland, 18–20 March 2020; Advances in Intelligent Systems and Computing. Springer: Berlin/Heidelberg, Germany, 2020. [Google Scholar]
  83. Głębocki, R.; Jacewicz, M. Parametric Study of Guidance of a 160-Mm Projectile Steered with Lateral Thrusters. Aerospace 2020, 7, 61. [Google Scholar] [CrossRef]
  84. Szklarski, A.; Głębocki, R.; Jacewicz, M. Impact Point Prediction Guidance Parametric Study for 155 mm Rocket Assisted Artillery Projectile with Lateral Thrusters. Arch. Mech. Eng. 2020, 67, 31–56. [Google Scholar] [CrossRef]
  85. Matsumoto, M.; Nishimura, T. Mersenne Twister. ACM Trans. Model. Comput. Simul. 1998, 8, 3–30. [Google Scholar] [CrossRef]
  86. Charles, C.; Jean-Francois, G.; Abolfazl, M.; Ugo, C.; Sofiane, A. Applying Robust Design Methodology to a Quadrotor Drone. In Proceedings of the 21st International Conference on Engineering Design, ICED 17, Vancouver, BC, Canada, 21–25 August 2017. [Google Scholar]
  87. Cabreira, T.M.; Franco, C.D.; Ferreira, P.R.; Buttazzo, G.C. Energy-Aware Spiral Coverage Path Planning for UAV Photogrammetric Applications. IEEE Robot. Autom. Lett. 2018, 3, 3662–3668. [Google Scholar] [CrossRef]
Figure 1. The quadrotor used in the experiments (CAD model).
Figure 1. The quadrotor used in the experiments (CAD model).
Energies 15 07136 g001
Figure 2. Moments of inertia measurements using trifilar pendulum (a) I x x (b) I y y (c) I z z .
Figure 2. Moments of inertia measurements using trifilar pendulum (a) I x x (b) I y y (c) I z z .
Energies 15 07136 g002
Figure 3. Coordinate systems overview.
Figure 3. Coordinate systems overview.
Energies 15 07136 g003
Figure 4. Rotors configuration (top view).
Figure 4. Rotors configuration (top view).
Energies 15 07136 g004
Figure 5. Modeled and manufacturer thrust and torque comparison for hover conditions.
Figure 5. Modeled and manufacturer thrust and torque comparison for hover conditions.
Energies 15 07136 g005
Figure 6. Correction factors of (a) thrust, (b) torque, and (c) drag.
Figure 6. Correction factors of (a) thrust, (b) torque, and (c) drag.
Energies 15 07136 g006
Figure 7. Automatic flight control system structure.
Figure 7. Automatic flight control system structure.
Energies 15 07136 g007
Figure 8. Laboratory tests of the battery discharge.
Figure 8. Laboratory tests of the battery discharge.
Energies 15 07136 g008
Figure 9. Quadrotor during the flight tests: (a) drone ready to fly and (b) take-off.
Figure 9. Quadrotor during the flight tests: (a) drone ready to fly and (b) take-off.
Energies 15 07136 g009
Figure 10. The comparison between the real and simulated trajectories: (a) 3D view of the trajectory and (b) 2D plot of the drone trajectory.
Figure 10. The comparison between the real and simulated trajectories: (a) 3D view of the trajectory and (b) 2D plot of the drone trajectory.
Energies 15 07136 g010
Figure 11. Linear velocities: (a) north–south, (b) west–east, and (c) vertical.
Figure 11. Linear velocities: (a) north–south, (b) west–east, and (c) vertical.
Energies 15 07136 g011
Figure 12. Angular rates: (a) roll, (b) pitch, and (c) yaw.
Figure 12. Angular rates: (a) roll, (b) pitch, and (c) yaw.
Energies 15 07136 g012
Figure 13. Quadrotor position: (a) north–south, (b) west–east, and (c) altitude.
Figure 13. Quadrotor position: (a) north–south, (b) west–east, and (c) altitude.
Energies 15 07136 g013
Figure 14. Euler angles: (a) roll, (b) pitch, and (c) yaw.
Figure 14. Euler angles: (a) roll, (b) pitch, and (c) yaw.
Energies 15 07136 g014
Figure 15. Energetic model: (a) total energy consumed by the drone, (b) battery state of charge, (c) voltage, and (d) current.
Figure 15. Energetic model: (a) total energy consumed by the drone, (b) battery state of charge, (c) voltage, and (d) current.
Energies 15 07136 g015
Figure 16. The comparison of trajectories between reality and simulation: (a) 3D view of the trajectory and (b) 2D plot of the drone trajectory.
Figure 16. The comparison of trajectories between reality and simulation: (a) 3D view of the trajectory and (b) 2D plot of the drone trajectory.
Energies 15 07136 g016
Figure 17. Linear velocities: (a) north–south, (b) west–east, and (c) vertical.
Figure 17. Linear velocities: (a) north–south, (b) west–east, and (c) vertical.
Energies 15 07136 g017
Figure 18. Angular rates: (a) roll, (b) pitch, and (c) yaw.
Figure 18. Angular rates: (a) roll, (b) pitch, and (c) yaw.
Energies 15 07136 g018
Figure 19. Quadrotor position: (a) north–south, (b) west–east, and (c) altitude.
Figure 19. Quadrotor position: (a) north–south, (b) west–east, and (c) altitude.
Energies 15 07136 g019
Figure 20. Euler angles: (a) roll, (b) pitch, and (c) yaw.
Figure 20. Euler angles: (a) roll, (b) pitch, and (c) yaw.
Energies 15 07136 g020
Figure 21. Energetic model: (a) total energy consumed by the drone (b) battery state of charge, (c) voltage, and (d) current.
Figure 21. Energetic model: (a) total energy consumed by the drone (b) battery state of charge, (c) voltage, and (d) current.
Energies 15 07136 g021aEnergies 15 07136 g021b
Figure 22. Histogram and probability distribution fit of (a) the N–S leg (b), the E–W leg, and (c) turn (real data).
Figure 22. Histogram and probability distribution fit of (a) the N–S leg (b), the E–W leg, and (c) turn (real data).
Energies 15 07136 g022
Figure 23. Quadrotor trajectories: (a) case 1 and (b) case 2.
Figure 23. Quadrotor trajectories: (a) case 1 and (b) case 2.
Energies 15 07136 g023
Figure 24. Kernel density estimators of consumed energy: (a) case 1: all disturbances, (b) case 2: all disturbances, (c) case 1: wind disturbances, (d) case 2: wind disturbances, (e) case 1: mass disturbances, and (f) case 2: mass disturbances.
Figure 24. Kernel density estimators of consumed energy: (a) case 1: all disturbances, (b) case 2: all disturbances, (c) case 1: wind disturbances, (d) case 2: wind disturbances, (e) case 1: mass disturbances, and (f) case 2: mass disturbances.
Energies 15 07136 g024
Figure 25. Histograms of consumed energy: (a) mean wind azimuth 0°, wind speed 4.5 m/s, (b) mean wind azimuth 180°, wind speed 4.5 m/s, (c) mean wind azimuth 0°, wind speed 9 m/s, and (d) mean wind azimuth 180°, wind speed 9 m/s.
Figure 25. Histograms of consumed energy: (a) mean wind azimuth 0°, wind speed 4.5 m/s, (b) mean wind azimuth 180°, wind speed 4.5 m/s, (c) mean wind azimuth 0°, wind speed 9 m/s, and (d) mean wind azimuth 180°, wind speed 9 m/s.
Energies 15 07136 g025
Table 1. Quadrotor parameters.
Table 1. Quadrotor parameters.
Parameter [Unit]SymbolValue
Drone take-off mass (kg) m 4.689
Moments of inertia (kg·m2) I x x 0.075716
I y y 0.084124
I z z 0.126437
Products of inertia (kg·m2) I x y , I y z , I x z 0
First propeller position vector (m) r R 1 = x R 1 , y R 1 , z R 1 T [0.2475, 0.2475, −0.074]
Second propeller position vector (m) r R 2 = x R 2 , y R 2 , z R 2 T [−0.2475, 0.2475, −0.074]
Third propeller position vector (m) r R 3 = x R 3 , y R 3 , z R 3 T [−0.2475, −0.2475, −0.074]
Fourth propeller position vector (m) r R 4 = x R 4 , y R 4 , z R 4 T [0.2475, −0.2475, −0.074]
Propeller diameter (m) D p 0.4572
Propeller pitch (m)0.1549
Propeller mass (kg)0.017
Propeller moment of inertia (kg) I z p 0.0002964
Battery mass (kg)2.013
Max. propeller angular rate (RPM)6000
Propeller angular rate at hover (119 m above sea level) (RPM)3203.82
Maximum flight altitude (m)6500
Flight endurance (500 g payload) (min)≥71
Flight endurance (1000 g payload) (min)≥60
Table 2. Rotor chord, chord center shift, and blade twist distribution along the rotor span.
Table 2. Rotor chord, chord center shift, and blade twist distribution along the rotor span.
r (mm) c (mm) Δ c (mm) γ (°)
2519.75−0.121.16
37.525.130.6423.50
5032.83.2323.09
62.540.335.9620.83
7544.157.3717.00
87.543.57.0414.84
10041.696.2113.49
12537.264.4612.17
15032.872.811.30
17528.621.379.92
20023.920.18.69
21021.99−0.198.42
22016.720.997.15
22511.193.467.83
Table 3. Motor model parameters used in the numerical simulation.
Table 3. Motor model parameters used in the numerical simulation.
Parameter (Unit)SymbolValue
torque constant of the electric motor (Nm/A) K T 0.9 × 10–2
power required by the onboard subsystems (J/s) E s u b 11
rotor thrust coefficient (N/RPM2) k f 9.32 × 10–5
rotor torque coefficient (Nm/RPM2) k m 9.37 × 10–6
viscous damping coefficient of the motor (Nm·s/rad) D v 0.17 × 10–3
reference linear dimension (m) d 1
Table 4. Battery model parameters used in the numerical experiment.
Table 4. Battery model parameters used in the numerical experiment.
Parameter (Unit) *SymbolValue
constant voltage (V) E 0 16.8
polarization constant (V/Ah) K 0.038603
battery capacity (Ah)/(Wh) Q 29.7/439.6
exponential voltage (V) A 0.2468
exponential capacity (Ah)−1 B 30
internal battery resistance (Ω) R 0.025
* Maximum Continuous discharge current: 45 (A), Working and storage temperature from −10 °C up to +65 °C. 14.8 (V), maximum charge current: 15 (A).
Table 5. Energy consumed by the drone in various maneuvers (real data).
Table 5. Energy consumed by the drone in various maneuvers (real data).
Pattern PartMean (kWs)Up. CI Limit (kWs)Low. CI Limit (kWs)
N–S leg14.70913.97515.442
E–W leg22.45721.69623.218
Turn5.1375.4394.835
Table 6. Root mean square errors.
Table 6. Root mean square errors.
ParameterRMSE (Case 1)RMSE (Case 2)TIC (Case 1)TIC (Case 2)
U n 0.5482760.5795510.1040420.057415
V n 0.5932530.5229130.1215440.212000
W n 0.1761070.1137290.1012740.746461
P 0.1229860.1451710.7297820.701991
Q 0.1588960.1640720.7493230.723037
R 0.0698340.0752420.1457920.143996
x n 3.4001826.4491900.0310590.020022
y n 3.6732671.6588880.0498050.006499
z n 0.3943010.0380040.0081080.000633
Φ 3.2522053.0700440.3763300.279653
Θ 3.3980183.5906910.2647160.223927
Ψ 2.5939322.7424820.0117900.005076
Table 7. Input parameters for the Monte-Carlo simulation.
Table 7. Input parameters for the Monte-Carlo simulation.
ParameterUnit Nominal   Value   μ (Case 1) Nominal   Value   μ (Case 2) Standard   Deviation   σ
m kg4.6894.6890.05
I x x kg·m20.0757160.0757160.005
I y y kg·m20.0841240.0841240.005
I z z kg·m20.1264370.1264370.005
U n 0 m/s−0.1194184.4531110.5
V n 0 m/s−0.1360430.3052090.5
W n 0 m/s−0.377093−0.0119630.5
P 0 °/s4.846231−7.6675913
Q 0 °/s0.2170453.0239153
R 0 °/s1.8863097.6528483
x n 0 m−2.904065−48.2280103
y n 0 m0.151600−117.6519003
z n 0 m−0.095660−30.0403233
Φ 0 °−0.570000−2.6100001
Θ 0 °0.3900000.2800001
Ψ 0 °198.020000289.0500005
V W t o t m/s4.571
Ψ W °22022510
Table 8. Statistics for the obtained data.
Table 8. Statistics for the obtained data.
CaseMean (kWs)Minimum (kWs)Maximum (kWs)
1 (all disturbances)48.674747.268849.9513
1 (wind only)48.694948.633548.9322
1 (mass, moments of inertia)48.656347.250049.8289
2 (all disturbances)542.9812528.1782558.7417
2 (wind only)543.3368539.3126548.5094
2 (mass, moments of inertia)543.2884527.8764556.1354
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Jacewicz, M.; Żugaj, M.; Głębocki, R.; Bibik, P. Quadrotor Model for Energy Consumption Analysis. Energies 2022, 15, 7136. https://doi.org/10.3390/en15197136

AMA Style

Jacewicz M, Żugaj M, Głębocki R, Bibik P. Quadrotor Model for Energy Consumption Analysis. Energies. 2022; 15(19):7136. https://doi.org/10.3390/en15197136

Chicago/Turabian Style

Jacewicz, Mariusz, Marcin Żugaj, Robert Głębocki, and Przemysław Bibik. 2022. "Quadrotor Model for Energy Consumption Analysis" Energies 15, no. 19: 7136. https://doi.org/10.3390/en15197136

APA Style

Jacewicz, M., Żugaj, M., Głębocki, R., & Bibik, P. (2022). Quadrotor Model for Energy Consumption Analysis. Energies, 15(19), 7136. https://doi.org/10.3390/en15197136

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop