Development of Path Planning Tool for Unmanned System Considering Energy Consumption

Abstract

Increasing the flight endurance of unmanned aerial vehicles (UAVs) has received attention recently. To solve this problem, two research topics have generally appeared: Shortest-path planning (SPP) and remaining-flying-range estimation. In this work, energy-efficient path planning by considering the distance between waypoint nodes, the minimum and maximum speed of the UAV, the weight of the UAV, and the angle between two intersecting edges is proposed. The performances of energy-efficient path planning (EEPP) and generic shortest-path planning are compared using extended-Kalman-filter-based state-of-charge and state-of-power estimation. Using this path-planning tool and considering energy consumption during flight operation, two different path plans can be obtained and compared in advance so that the operator can decide which path to choose by consulting a comparison chart. According to the experimental results, the EEPP algorithm results in 0.96% of improved SOC leftover and 11.03 ($\mathrm{W}$) of lowered SOP compared to the SPP algorithm.

1. Introduction

As unmanned aerial vehicle (UAV) technology matures, the longtime flight endurance capability of UAVs is gaining attention [1,2]. Because the energy density ($\mathrm{W}/\mathrm{kg}$) of the battery pack (in particular, Li-ion) is closely correlated with the overall flight time of the UAV, it is recommended to choose cells possessing as high an energy density as possible [3,4]. However, in contrast to electric vehicles (EVs), the operation times of electric aircraft, including the UAV, drastically decrease as the system weights increase. Hence, technology development for increasing UAV flight time has been delayed for a long time.

Lately, the EV industry has been growing quickly, and automotive industry researchers have focused their attention on accurately predicting the energy leftover so that the driving range can be extended. Accurate energy leftover monitoring (state of charge (SOC) estimation in the EV industry), aims to extend the driving range, and it is related to shortest-path planning (SPP), which is essential to saving energy. Both accurate SOC estimation and SPP should ultimately enhance the fuel economy of EVs, and the same philosophy applies to UAVs.

There have been numerous studies regarding path planning. Chang et al. proposed a method to achieve an effective delivery route for trucks carrying delivery drones with moving shift-weights-based K-means clustering and traveling-salesman-problem modeling [5]. Rubio et al. presented adaptive path planning for multiple UAVs using evolutionary computation combined with market-based cooperation strategies [6].

Further, there has been research regarding the accurate estimation of the remaining driving or flying range. Oliva et al. suggested a model-based approach by combining unscented filtering and Markov chains to predict the remaining driving range (RDR) [7]. Ondruska et al. proposed the efficient computation of RDR confidence level using a feature-based linear regression framework with a varying probabilistic attainability map in real-time [8]. Ferreira et al. presented an interaction method between a driver and EV to estimate EV range using a regression-model-based data mining approach [9].

Unlike the previously listed works, this study contributes to the energy-efficient path planning (EEPP) of the UAV to ultimately increase the flight endurance compared with the general SPP. To verify the improved fuel efficiency, the concept of the extended Kalman filter (EKF)-based SOC estimation method is adopted. To mimic cell behavior, a physics-based equivalent circuit model (ECM) is used [10,11,12,13].

The overall mission hierarchy to achieve an energy-efficient path proposed in this paper is as follows.

First, with the given satellite map, the EEPP algorithm generates the most energy-efficient path with step-by-step processes; (1) obstacle detection, (2) extract obstacle, (3) obstacle separation, (4) add buffer zone, (5) Voronoi, (6) filter out useless candidates, (7) 1D optimal control, and (8) SPP or EEPP selection. Here, the format of the output path has both Euclidean and GPS coordinates. In contrast to the SPP algorithm, which uses Dijkstra algorithm, the EEPP algorithm uses three step-by-step processes; (1) find all trajectory candidates, (2) calculate total unit energy consumption (UEC) of all trajectory candidates, and (3) find a trajectory having the minimum UEC.

Second, either flight simulation with equivalent circuit model (ECM) of the Li-Ion battery pack or the flight experiment with a real Li-Ion battery pack is performed using the previously obtained energy-efficient path and gives out the voltage and current profiles of each cell.

Third, a battery state estimator uses the voltage and current profiles of cells to calculate the SOC. In particular, the battery state estimator uses the EKF to calculate the SOC.

The flow of this paper is as follows: In Section 2, the overall mission hierarchy, including four subparts (path planning, flight simulation, flight experiment, and battery state estimation), is described. In Section 3, the path planning, including the SPP and EEPP algorithms, is explained. Section 4 covers the mathematical formulation of the UAV and battery pack. Section 5 presents the simulation and experimental setups, while Section 6 shows the simulation and experiment results. Section 7 contains the conclusion of this paper and describes future works.

2. Mission Hierarchy

Figure 1 depicts the overall mission hierarchy to achieve the most energy-efficient path and corresponding estimated SOC value.

First, using a Google Earth satellite photograph, and the takeoff and landing locations, the most energy-efficient path is determined by the path-planning part. Here, the path is saved in both the Cartesian system unit, $Trajector{y}_{xyz}$, and the geographic coordinate system unit, $Trajector{y}_{GPS}$.

Second, the achieved path is utilized for the indoor flight simulation and outdoor flight experiment to calculate the cell open-circuit voltage (OCV), $v_{ocv}$, and load current, $i_o$, data.

Third, the achieved cell voltage and current data are utilized for the battery state estimation to calculate the SOC. Both ${\mathrm{SOC}}_{Sim}$ and ${\mathrm{SOC}}_{Exp}$ are calculated to compare results.

With the achieved SOC value, one can finally determine whether the EEPP algorithm results in an energy-efficient path.

3. Path Planning

Figure 2 displays a magnified version of the path-planning part of Figure 1. This part results in either the shortest path or energy-efficient path. The detailed description of each subpart is given elsewhere [14].

3.1. SPP

The Dijkstra algorithm is used for the SPP algorithm, which utilizes the weights of each edge between nodes to calculate the minimum sum of the total distance from the start node to the goal node [15].

The overall flowchart of the SPP algorithm is shown in Figure 3.

3.2. EEPP

In contrast to the general SPP algorithm, which only accounts for the distance between nodes, the EEPP algorithm calculates the energy-efficient path by considering factors including the distance between nodes, the minimum and maximum speed of the UAV, and the weight of the UAV as shown in Figure 4.

The overall amount of battery consumption of the UAV, $B_o$, would be

$$B_o=B_t+B_h+B_a+B_n+B_d+B_c+B_l,$$

(1)

where $B_t$ is the battery consumption during takeoff ($\mathrm{W}$), $B_h$ is the battery consumption during hovering ($\mathrm{W}$), $B_a$ is the battery consumption during acceleration ($\mathrm{W}$), $B_n$ is the battery consumption during normal flight ($\mathrm{W}$), $B_d$ is the battery consumption during deceleration ($\mathrm{W}$), $B_c$ is the battery consumption during the circular flight ($\mathrm{W}$), and $B_l$ is the battery consumption during landing ($\mathrm{W}$).

Because the takeoff and landing occur in a short period and the amount of battery consumption is insignificant compared with the total consumption, $B_t$ and $B_l$ are assumed to be zero, and Equation (1) becomes

$$B_o\approx B_h+B_a+B_n+B_d+B_c.$$

(2)

If there are $n_e$ number of edges in the given path,

$$B_o\approx {\displaystyle \sum _{i=1}^{n_e}\left(B_h+B_a+B_n+B_d+B_c\right)}.$$

(3)

From the equation of a single rotor helicopter in hover [16,17], the required hovering power of the UAV can be derived as

$$B_h=m^{\frac{3}{2}}\sqrt{\frac{g^3}{2\rho \zeta n_r}},$$

(4)

where $m$ is the mass of the UAV ($\mathrm{kg}$), $g$ is the gravitational acceleration ($\mathrm{m}/{\mathrm{s}}^2$), $\rho$ is the fluid density of air ($\mathrm{kg}/{\mathrm{m}}^3$), $\zeta$ is the total area of the spinning rotor disk (${\mathrm{m}}^2$), and $n_r$ is the number of rotors (no unit).

The required power for the acceleration can be derived as

$$B_a=\frac{1}{2t_a}{\left({\displaystyle \underset{0}{\overset{t_a}{\int }}{a}_{\max }dt}\right)}^2=\frac{1}{2}m{a}_{\max }^2{t}_a,$$

(5)

where $a_{max}$ is the maximum acceleration ($\mathrm{m}/{\mathrm{s}}^2$), and $t_a$ is the time taken for the acceleration ($\mathrm{s}$).

The required power for the deceleration can be derived as

$$B_d=\frac{1}{2t_d}{\left({\displaystyle \underset{0}{\overset{t_d}{\int }}{a}_{\min }dt}\right)}^2,$$

(6)

where $a_{min}$ is the minimum acceleration ($\mathrm{m}/{\mathrm{s}}^2$), and $t_d$ is the time taken for the deceleration ($\mathrm{s}$). Here, with assumptions that $t_a=t_d$ and $a_{max}=-a_{min}$, one can conclude $B_a=B_d$.

The required power for the normal (constant speed) flight can be derived as

$$B_n=\frac{1}{2t_n}m{v}_{\max }^2.$$

(7)

Here, $v_{max}$ is the maximum velocity ($\mathrm{m}/{\mathrm{s}}^2$), and $t_n$ is the time taken for the normal flight ($\mathrm{s}$).

The required power for the circular flight can be derived as

$$B_c=\frac{1}{2t_c}m{v}_{\min }^2$$

(8)

where $v_{min}$ is the minimum velocity ($\mathrm{m}/{\mathrm{s}}^2$), and $t_c$ is the time taken for the circular flight ($\mathrm{s}$).

The unknown value, $t_a$, in Equation (5) can be calculated as

$$t_a=t_d=\frac{v_{\max }-v_{\min }}{a_{\max }}.$$

(9)

The unknown value, $t_n$, in Equation (7) can be calculated as

$$\frac{1}{2}{t}_c=\frac{1}{v_{\min }}\left\{\frac{v_{\min }^2}{a_{\max }}\tan \left(\frac{\pi }{2}-\frac{\alpha }{2}\right)\right\},$$

(10)

therefore,

$$t_c=2\frac{v_{\min }}{a_{\max }}\tan \left(\frac{\pi }{2}-\frac{\alpha }{2}\right),$$

(11)

where $\alpha$ is the intersection angle of two edges, which satisfies $0\le \alpha \le \pi$.

Because

$$\begin{array}{l}{b}_1=t_a{v}_{\min }+\frac{1}{2}{t}_a\left(v_{\max }-v_{\min }\right),\\ b_2=\frac{v_{\min }^2}{a_{\max }}\tan \left(\frac{\pi }{2}-\frac{\alpha }{2}\right),\end{array}$$

(12)

the unknown value, $t_c$, in Equation (8) can be calculated as

$$t_n=\frac{1}{v_{\max }}\left\{\Vert W_{i-1},W_i\Vert -2\left(b_1+b_2\right)\right\},$$

(13)

where $W_i$ is the waypoint (no unit), $b_1$ is the distance in which the UAV accelerates or decelerates ($\mathrm{m}$), and $b_2$ is the distance at which the UAV flies in a circular motion ($\mathrm{m}$).

If Equations (4)–(13) are combined into Equation (3), then

$$B_0={\displaystyle \sum _{i=1}^{n_e}\left[m^{\frac{3}{2}}\sqrt{\frac{g^3}{2\rho \zeta n_r}}+m{a}_{\max }\left(v_{\max }-v_{\min }\right)+\frac{m{v}_{\max }^3}{2\left\{\Vert W_{i-1},W_i\Vert -2\left(b_1+b_2\right)\right\}}+\frac{m{v}_{\min }{a}_{\max }}{4\tan \left(\frac{\pi }{2}-\frac{\alpha }{2}\right)}\right]}.$$

(14)

To calculate the total unit energy consumption (UEC) ($\mathrm{Wh}$), $U_b$, Equation (14) is converted as

$$B_0=\frac{t_t}{3600}{\displaystyle \sum _{i=1}^{n_e}\left[m^{\frac{3}{2}}\sqrt{\frac{g^3}{2\rho \zeta n_r}}+m{a}_{\max }\left(v_{\max }-v_{\min }\right)+\frac{m{v}_{\max }^3}{2\left\{\Vert W_{i-1},W_i\Vert -2\left(b_1+b_2\right)\right\}}+\frac{m{v}_{\min }{a}_{\max }}{4\tan \left(\frac{\pi }{2}-\frac{\alpha }{2}\right)}\right]}$$

(15)

where $t_t=t_a+t_n+t_d+t_c$.

The overall flowchart of the EEPP algorithm is shown in Figure 5.

4. Mathematical Modeling of the UAV and Battery Pack

4.1. UAV

The configuration and governing control inputs of the quadrotor are described in Figure 6, where $X_E$, $Y_E$, and $Z_E$ represent earth coordinates, $x_q$, $y_q$, and $z_q$ represent quadrotor coordinates, and $F_i$ represents a thrust force generated by the i-th rotor in the indicated direction.

To extract current and voltage profiles along the predetermined path obtained from the previous section, an enhanced quadrotor dynamic model can be represented as

$$\begin{array}{l}{w}_1^2=b{\left(T+{\tau }_x+{\tau }_y+{\tau }_z\right)}^{-1}\\ w_2^2=lb\left(-{\tau }_x+{\tau }_z\right)\\ w_3^2=lb\left(-T+{\tau }_y\right)\\ w_4^2=d\left(-T+{\tau }_x-{\tau }_y+{\tau }_z\right)\end{array}$$

(16)

where $T$ is the total upward thrust of the UAV ($\mathrm{N}$), b is the thrust factor, and $w_i$ is the i-th rotor speed ($\mathrm{rad}/\mathrm{s}$), $l$ is the lever length from the center of gravity (CG) to rotors ($\mathrm{m}$), d is the drag factor, and ${\tau }_x$, ${\tau }_y$, and ${\tau }_z$ are torques applied on the UAV’s body along x-, y-, and z-axes ($\mathrm{N}·\mathrm{m}$) [21,22].

Here, $T$, ${\tau }_x$, ${\tau }_y$, and ${\tau }_z$ can be calculated using proportional–integral–derivative (PID) controllers as

$$\begin{array}{l}T=k_{p_T}\left(z^{*}-z\right)+k_{i_T}{\displaystyle \underset{t_1}{\overset{t_2}{\int }}\left(z^{*}-z\right)dt}+k_{d_T}\left({\dot{z}}^{*}-\dot{z}\right)+\sqrt{W/4b},\\ {\tau }_x=k_{p_x}\left({\varphi }^{*}-\varphi \right)+k_{i_x}{\displaystyle \underset{t_1}{\overset{t_2}{\int }}\left({\varphi }^{*}-\varphi \right)dt}+k_{d_x}\left({\dot{\varphi }}^{*}-\dot{\varphi }\right),\\ {\tau }_y=k_{p_y}\left({\theta }^{*}-\theta \right)+k_{i_y}{\displaystyle \underset{t_1}{\overset{t_2}{\int }}\left({\theta }^{*}-\theta \right)dt}+k_{d_y}\left({\dot{\theta }}^{*}-\dot{\theta }\right),\\ {\tau }_z=k_{p_z}\left({\psi }^{*}-\psi \right)+k_{i_z}{\displaystyle \underset{t_1}{\overset{t_2}{\int }}\left({\psi }^{*}-\psi \right)dt}+k_{d_z}\left({\dot{\psi }}^{*}-\dot{\psi }\right),\end{array}$$

(17)

where $k_p$, $k_i$, and $k_d$ are PID controller gains, $\varphi$, $\theta$, and $\psi$ are roll, pitch, and yaw angles ($\mathrm{rad}$), uppercase letter * represents reference data, and $W$ is the UAV weight ($\mathrm{kg}$).

The full quadrotor UAV dynamic model can be written as

$$\begin{array}{l}\ddot{x}=\left(\cos \varphi \sin \theta \cos \psi +\sin \varphi \sin \psi \right)\frac{1}{m}b\left(w_1^2+w_2^2+w_3^2+w_4^2\right),\\ \ddot{y}=\left(\cos \varphi \sin \theta \sin \psi -\sin \varphi \cos \psi \right)\frac{1}{m}b\left(w_1^2+w_2^2+w_3^2+w_4^2\right),\\ \ddot{z}=\left(\cos \varphi \cos \theta \right)\frac{1}{m}b\left(w_1^2+w_2^2+w_3^2+w_4^2\right)-g,\\ \ddot{\varphi }=\dot{\theta }\dot{\psi }\left(\frac{I_y-I_z}{I_x}\right)-\frac{J_r}{I_x}\dot{\theta }\left(w_1+w_3-w_2-w_4\right)+\frac{l}{I_x}b\left(w_4^2-w_2^2\right),\\ \ddot{\theta }=\dot{\varphi }\dot{\psi }\left(\frac{I_z-I_x}{I_y}\right)-\frac{J_r}{I_y}\dot{\varphi }\left(w_1+w_3-w_2-w_4\right)+\frac{l}{I_y}b\left(w_3^2-w_1^2\right),\\ \ddot{\psi }=\dot{\varphi }\dot{\theta }\left(\frac{I_x-I_y}{I_z}\right)+\frac{1}{I_z}d\left(w_2^2+w_4^2-w_1^2-w_3^2\right).\end{array}$$

(18)

Here, $J_r$ is the rotor inertia ($\mathrm{kg}·{\mathrm{m}}^2$).

Using the fact that power can be calculated by multiplying the thrust and velocity, one can calculate the required propulsive thrust power from the total upward thrust $T$ in Equation (15) as

$$P=T\cdot \Vert v_0\Vert .$$

(19)

Then, the load current can be calculated using

$$\frac{P}{v_t}=i_0$$

(20)

where $v_t$ is the terminal voltage ($\mathrm{V}$), and $i_0$ is the load current ($\mathrm{A}$).

The above equations are used for the operation of flight simulation, as shown in Figure 7.

4.2. Battery Pack

4.2.1. ECM

By using the two-ladder battery ECM and impedance voltage response, as shown in Figure 8, the $v_t$ can be determined as

$$\begin{array}{ll}{v}_t& =v_{ocv}-i_0{R}_0-v_{c1}-v_{c2}\\ & =v_{ocv}-i_0{R}_0-i_0{R}_1\left(1-e^{-\Delta t/{\tau }_1}\right)-i_0{R}_2\left(1-e^{-\Delta t/{\tau }_2}\right)\end{array}$$

(21)

where $v_{ocv}$ is the OCV voltage ($\mathrm{V}$), $R_0$ is the resistor ($\mathsf{\Omega }$), $v_{ci}$ is the voltage across the i-th RC-network ($\mathrm{V}$), ${\tau }_i$ is the time constant ($\mathrm{s}$), which is calculated as ${\tau }_i=R_i{C}_i$ in which $R_i$ is the resistor ($\mathsf{\Omega }$) and $C_i$ is the capacitor ($\mathrm{F}$) of the i-th ladder, and $t$ is the time ($\mathrm{s}$).

4.2.2. TBLI Parameter Identification

The ECM parameters, including $R_0$, $R_1$, $C_1$, $R_2$, and $C_2$, are extracted at SOC 10% intervals using the high-performance pulse characterization current pattern, as shown in Appendix A [12]. Here, the table-based linear interpolation (TBLI) method is used to determine ECM parameters corresponding to varying SOC values. Although $R_0$, $R_1$, $C_1$, $R_2$, and $C_2$ are the functions of the SOC, C-rate, and temperature, only the SOC variance effect is considered in this research.

4.2.3. SOC and SOP State Estimation

The SOC estimation is performed using the EKF for each cell. For the EKF algorithm application, Equation (21) and SOC calculation based on current integration method are fused and converted into the discrete-time state space model as

$$\begin{array}{l}{\left[\begin{array}{c}{v}_{c1}\\ v_{c2}\\ z\end{array}\right]}_{k+1}=\left[\begin{array}{ccc}{e}^{-\Delta t/{\tau }_1}& 0& 0\\ 0& e^{-\Delta t/{\tau }_2}& 0\\ 0& 0& 1\end{array}\right]{\left[\begin{array}{c}{v}_{c1}\\ v_{c2}\\ z\end{array}\right]}_k+{\left[\begin{array}{c}{R}_1\left(1-e^{-\Delta t/{\tau }_1}\right)\\ R_2\left(1-e^{-\Delta t/{\tau }_2}\right)\\ \eta \Delta t/\left(3600Q_{norm}\right)\end{array}\right]}_k{i}_{0,k}+{\left[\begin{array}{c}{n}_{w1}\\ n_{w2}\\ n_{w3}\end{array}\right]}_k,\\ {\left[\begin{array}{c}{v}_{c1}\\ v_{c2}\\ z\end{array}\right]}_k=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right]{\left[\begin{array}{c}{v}_{c1}\\ v_{c2}\\ z\end{array}\right]}_k+{\left[\begin{array}{c}0\\ 0\\ 0\end{array}\right]}_k{i}_{0,k}+{\left[\begin{array}{c}{n}_{v1}\\ n_{v2}\\ n_{v3}\end{array}\right]}_k,\end{array}$$

(22)

where $z$ is the SOC ($\%/100$), $\Delta t$ is 0.1 s, and $n_{wi}$ and $n_{vi}$ are system and measurement noises ($\mathrm{V}$), respectively. Here, the term $z$ can be calculated as $z\left(t\right)=Q_{\mathrm{init}}+\eta i_0\left(t\right)/Q_{\mathrm{norm}}\left(t\right)$ where $Q_{\mathrm{init}}$ is the initial battery capacity ($\mathrm{Ah}$), $\eta$ is the battery efficiency (assumed as 1), and $Q_{\mathrm{norm}}$ is the current battery capacity ($\mathrm{Ah}$).

With an initial estimation vector, ${\widehat{x}}_0$, and covariance vector, $P_0$, the estimation vector and covariance vector can be propagated using

$$\begin{array}{l}{\widehat{x}}_{k+1}^{-}=f\left({\widehat{x}}_k\right),\\ P_{k+1}^{-}=A_{k+1}{P}_k{A}_{k+1}^T+Q_k.\end{array}$$

(23)

The Kalman gain, $K_k$, can be calculated as

$$K_{k+1}=P_{k+1}^{-}{H}_{k+1}^T{\left(H_{k+1}{P}_{k+1}^{-}{H}_{k+1}^T+R_{k+1}\right)}^{-1}.$$

(24)

With the updated Kalman gain, $K_{k+1}$, and updated measurement, $z_{k+1}$, the updated estimation vector, ${\widehat{x}}_{k+1}$, can be calculated as

$${\widehat{x}}_{k+1}={\widehat{x}}_{k+1}^{-}+K_{k+1}\left(z_{k+1}-h_{k+1}\left({\widehat{x}}_{k+1}^{-}\right)\right),$$

(25)

and the updated covariance vector, $P_{k+1}$, can be calculated as

$$P_{k+1}=P_{k+1}^{-}-K_{k+1}{H}_{k+1}{P}_{k+1}^{-}.$$

(26)

The final SOC value of the battery pack is conservatively chosen as the minimum value of each cell as

$$SO{C}_{pack,k}=\min \left(z_k^1,\dots ,z_k^n\right)\times 100,$$

(27)

where $SO{C}_{pack,k}$ is the SOC (%) of a battery pack at the k-th time step.

Regarding the state of power (SOP) estimation, the maximum discharging current, $i_{max}^{dis}$, can be calculated as

$$i_{\max }^{dis}=C_{rate}{Q}_{norm}$$

(28)

where $C_{rate}$ is the C-rate of the cell (no unit).

The maximum discharging current considering the SOC, $i_{max,k}^{dis,SOC\left(i\right)}$, can be calculated as

$$i_{\max ,k}^{dis,SOC(i)}=\left\{\left(z_k^{(i)}-3{\sigma }_{z,k}^{(i)}\right)-z_{\min }\right\}{C}_{rate}{Q}_{norm}$$

(29)

where $\mathsf{\sigma }$ is the standard deviation of the cell SOCs ($\%/100$), $z_{min}$ is the minimum SOC ($\%/100$), $\eta$ is the battery efficiency (assumed as one), and $T_k-T_o$ is the time flown from the beginning ($\mathrm{s}$). Here, $3\mathsf{\sigma }$ was used for the 99.7% confidence interval of the cell SOC estimation values.

The maximum discharging current considering the cell voltage, $i_{max,k}^{dis,volt\left(i\right)}$, can be calculated as

$$i_{\max ,k}^{dis,volt(i)}=\left|i_{0,k}\right|\left\{v_{ocv}\left(z_k^{(i)}\right)-v_{\min }^{(i)}\right\}\left\{\frac{1}{\left|v_{ocv}\left(z_k^{(i)}-v_{\min }^{(i)}\right)\right|}+\frac{1}{v_{\max }^{(i)}-\left\{v_{ocv}\left(z_k^{(i)}\right)-v_{\min }^{\left(i\right)}\right\}}\right\},$$

(30)

where $v_{ocv}$ is the open-circuit voltage (OCV) (V), and $v_{min}$ is the minimum cell voltage (V). Here, $v_{min}^{\left(1\right)}=3.679\mathrm{V},$ $v_{min}^{\left(2\right)}=3.660\mathrm{V},$ $v_{min}^{\left(3\right)}=3.663\mathrm{V},$ $v_{min}^{\left(4\right)}=3.626\mathrm{V},$ $v_{max}^{\left(1\right)}=4.171\mathrm{V},$ $v_{max}^{\left(2\right)}=4.176\mathrm{V},$ $v_{max}^{\left(3\right)}=4.174\mathrm{V},$ $v_{max}^{\left(4\right)}=4.174\mathrm{V},$ $z_{min}=0,$ $C_{rate}=18 \mathrm{C},$ $N_p=1,$ $N_s=4,$ and $Q_{norm}=3.3 \mathrm{Ah}$.

With the previously found $i_{max}^{dis}$, $i_{max,k}^{dis,SOC\left(i\right)}$, and $i_{max,k}^{dis,volt\left(i\right)}$, one can calculate $i_{max,k}^{dis}$ as

$$i_{\max ,k}^{dis}=\min \left(i_{\max }^{dis},i_{\max ,k}^{dis,SOC(i)},i_{\max ,k}^{dis,volt(i)}\right)$$

(31)

and can calculate $P_{max,k}^{dis\left(i\right)}$ as

$$P_{\max ,k}^{dis(i)}=i_{\max ,k}^{dis}\left\{v_{ocv}\left(z_k^{(i)}\right)-i_{\max ,k}^{dis}\left|\frac{v_{ocv}\left(z_k^{(i)}\right)-v_{\min }}{i_{0,k}}\right|\right\}.$$

(32)

Then, the maximum power that each cell could perform, $P_{max}^{dis\left(i\right)}$, is calculated as

$$P_{\max }^{dis(i)}=v_{\max }^{(i)}{i}_{\max }^{dis}$$

(33)

and compared with the maximum power that each cell performs during the discharging process, $P_{max,k}^{dis\left(i\right)}$, as

$$P_{\max ,k}^{dis(i)}=N_p\min \left(N_s{P}_{\max }^{dis(i)},{\displaystyle \sum _{i=1}^{N_s}{i}_{\max ,k}^{dis}\left(v_{ocv}\left(z_k^{(i)}\right)-i_{\max ,k}^{dis}{R}_{dis}^{(i)}\right)}\right)$$

(34)

where the total initial resistance of each cell, $R_{dis}^{\left(i\right)}$, can be calculated as

$$R_{dis}^{(i)}=\frac{v_{ocv}\left(z_k^{(i)}\right)-v_{\min }}{i_{0,k}}.$$

(35)

Figure 9 presents the flowchart of the battery state estimation [23].

5. Simulation and Experimental Setups

5.1. Simulation Setup

All simulations were performed with a quadrotor vehicle with initial system properties $P_o=\left[0,0,0\right] \mathrm{m}$, ${\theta }_o=\left[0,0,0\right] \mathrm{rad}$, $w_o=\left[0, 0, 0\right] \mathrm{rad}/\mathrm{s}$, $v_o=\left[0, 0, 0\right] \mathrm{m}/\mathrm{s}$, $a_o=\left[0, 0, 0\right] \mathrm{m}/{\mathrm{s}}^2$, $\Vert v_{min}\Vert {}_2=0 \mathrm{m}/\mathrm{s}$, $\Vert v_{max}\Vert {}_2=5 \mathrm{m}/\mathrm{s}$, $\Vert a_{max}\Vert {}_2=0.5 \mathrm{m}/{\mathrm{s}}^2$, $Al{t}_{min}=3 \mathrm{m}$, and $Al{t}_n=30 \mathrm{m}$, where $P_o$, ${\theta }_o$, $w_o$, $v_o$, $a_o$, $Al{t}_{min}$, and $Al{t}_n$ are initial position, attitude, angular velocity, velocity, acceleration, minimum flight altitude, and normal flight altitude, respectively. Also, the time step used in MATLAB/Simulink is 0.01 s, installed in a desktop with the following system specifications; Intel Core i5-4590 CPU 3.40-$\mathrm{GHz}$ processor, 64-bit operating system, and 4.00-$\mathrm{GB}$ RAM. Every simulation flight was performed with a fully charged 4s 3300-$\mathrm{mAh}$ battery pack.

5.2. Experiment Setup

All experiments were performed with the parts shown in Figure 1—a PX4 Pro powered F450 quadrotor, shown in Figure 10a; Pixhawk 2.1 flight controller, as in Figure 10b; QGroundControl (QGC), shown in Figure 10c; Pixhawk 2.1 Here + RTK, as in Figure 10d; eLogger V4 voltage and current data logging system, as shown in Figure 10e; and a 4s 3300-$\mathrm{mAh}$ battery pack, shown in Figure 10f—for the outdoor autopilot experiment with a previously achieved energy-efficient path. In particular, the Pixhawk 2.1 Here + RTK is used for the centimeter-level autopilot path following accuracy.

The overall experiment setup is presented in Figure 11. On the day of the experiment, the wind speed was approximately 2 $\mathrm{m}/\mathrm{s}$.

For the experiment, the SPP and EEPP trajectories represented as the Cartesian system unit, $Trajector{y}_{xyz}$, obtained from the path-planning algorithm shown in Figure 2 are converted into the geographic coordinate system unit, $Trajector{y}_{GPS}$. The detailed description can be found in previous work [24].

In particular, the waypoint list of $Trajector{y}_{GPS}$ is saved as a text file following the Pixhawk waypoint file format, as depicted in

Table 1. Waypoint file format.

QGC WPL < VERSION>
<INDEX>	<CURRENT WP>	<COORD FRAME>	<COMMAND>	<PARAM1>	<PARAM2>
<PARAM3>	<PARAM4>	<PARAM5/X/	< PARAM6/Y/	< PARAM5/Z/	<AUTOCONTINUE>
		LONGITUDE>	LATITUDE >	ALTITUDE >

[25]. For the waypoints lists of SPP and EEPP trajectories, see Appendix B.

6. Simulation and Experimental Results

6.1. Simulation Result

From the flight simulation, the obtained current and voltage profiles are shown in Figure 12. Here, negative current represents discharging.

Energy consumption comparison results of the SPP and EEPP algorithms are listed in

Table 2. Comparison result (simulation).

	SPP		EEPP
	Prediction Value	Actual Value	Prediction Value	Actual Value
Total Trajectory Distance ($\mathrm{m}$)	204.63	193.30	208.76	201.08
Min Speed of UAV ($\mathrm{m}/\mathrm{s}$)	0	0	0	0
Max Speed of UAV ($\mathrm{m}/\mathrm{s}$)	5	0.34	5	0.33
Weight of UAV ($\mathrm{kg}$)	1.41	1.41	1.41	1.41
Total UEC ($\mathrm{Wh}$)	4.47	11.54	5.71	11.99
Algorithm Run Time ($\mathrm{s}$)	11.51		20.34
SOC Leftover ($\%$)	72.00		70.80
SOP Peak ($\mathrm{W}$)	566.20		543.47

and Figure 13.

The percentage difference between the actual values of SPP and EEPP is shown in

Table 3. Comparison results in percentage (simulation).

	$\frac{\mathit{S}\mathit{P}\mathit{P}-\mathit{E}\mathit{E}\mathit{P}\mathit{P}}{\mathit{S}\mathit{P}\mathit{P}}\times 100$$\left(\mathit{\%}\right)$
Total Trajectory Distance ($\mathrm{m}$)	−4.02
Min. Speed of UAV ($\mathrm{m}/\mathrm{s}$)	0
Max. Speed of UAV ($\mathrm{m}/\mathrm{s}$)	2.94
Weight of UAV ($\mathrm{kg}$)	0
Total UEC ($\mathrm{Wh}$)	−3.90
Algorithm Run Time ($\mathrm{s}$)	−76.72
SOC Leftover ($\%$)	1.67
SOP Peak ($\mathrm{W}$)	4.01

and Figure 14. According to the result, the SOC leftover is larger and the SOP peak is lower when using SPP.

6.2. Experimental Result

From the flight experiment, the obtained current and voltage profiles are shown in Figure 15.

Energy consumption comparison results of the SPP and EEPP algorithms are shown in

Table 4. Energy consumption comparison (experiment).

	SPP		EEPP
	Prediction Value	Actual Value	Prediction Value	Actual Value
Total Trajectory Distance ($\mathrm{m}$)	205.75	174.11	212.17	183.64
Min. Speed of UAV ($\mathrm{m}/\mathrm{s}$)	0	0	0	0
Max. Speed of UAV ($\mathrm{m}/\mathrm{s}$)	5	4.20	5	4.61
Weight of UAV ($\mathrm{kg}$)	1.41	1.41	1.41	1.41
Total UEC ($\mathrm{Wh}$)	4.47	7.94	5.71	7.41
Algorithm Run Time ($\mathrm{s}$)	11.51		20.34
SOC Leftover ($\%$)	83.73		84.53
SOP Peak ($\mathrm{W}$)	448.85		498.37

and Figure 16.

The percentage difference between the actual values of SPP and EEPP is shown in

Table 5. Comparison result in percentage (experiment).

	$\frac{\mathit{S}\mathit{P}\mathit{P}-\mathit{E}\mathit{E}\mathit{P}\mathit{P}}{\mathit{S}\mathit{P}\mathit{P}}\times 100$$\left(\mathit{\%}\right)$
Total Trajectory Distance ($\mathrm{m}$)	−5.47
Min. Speed of UAV ($\mathrm{m}/\mathrm{s}$)	0
Max. Speed of UAV ($\mathrm{m}/\mathrm{s}$)	−9.76
Weight of UAV ($\mathrm{kg}$)	0
Total UEC ($\mathrm{Wh}$)	6.68
Algorithm Run Time ($\mathrm{s}$)	−76.72
SOC Leftover ($\%$)	−0.96
SOP Peak ($\mathrm{W}$)	−11.03

and Figure 17. According to the result, the SOC leftover is larger and the SOP peak is lower when using EEPP.

The overall reference, simulation, and experiment flight trajectories of the SPP and EEPP algorithms are drawn on Google Earth, as shown in Figure 18. Here, yellow polygons represent imaginary obstacles that the UAV should avoid. The blue line represents the reference trajectory that the UAV should follow. The red line represents the simulation flight trajectory. The green line represents the experiment flight trajectory.

7. Conclusions

In this work, the EEPP algorithm is proposed and its performance is compared to the generic SPP algorithm using EKF-filter-based SOC estimation. Using the proposed path planning tool, two different path plans can be obtained and compared in advance. According to the experimental results, the EEPP algorithm results in an energy-efficient path.

Regarding the SOC, 0.96% is saved when using the EEPP algorithm, and this implies that the fuel economy is improved. Regarding the SOP, 11.03 ($\mathrm{W})$ is lowered, and this implies that the internal shock of the battery pack is reduced. Although the rate of improved fuel economy and internal shock of the battery pack is not as good in this study, these two key factors can be enhanced in mission planning with longer trajectories. Here, the difference between simulation and experiment results appears to be caused by factors such as environmental effects and a less mature quadrotor simulation model.

In the future, instead of 2D path planning, 3D path planning with an enhanced EEPP algorithm are planned to greatly improve the fuel economy of the UAV. Also, rather than the offline method presented here, the development of online real-time energy-efficient path planning based on the probabilistic method is planned.