Enhanced Potential Field-Based Collision Avoidance in Cluttered Three-Dimensional Urban Environments

Abstract

With the various applications of unmanned aerial vehicles (UAVs), the number of UAVs will increase in limited airspace, leading to an increased risk collision. To reduce such potential risk, this work proposes a collision avoidance strategy for UAVs using an enhanced potential field (EPF) approach in cluttered three-dimensional urban environments. Using the EPF formulated in a two-dimensional environment, the avoidance maneuvers for both horizontal and vertical planes are generated by introducing rotation matrices, and these maneuvers are combined by applying a weighting factor. The numerical simulations with various meaningful scenarios are conducted to validate the performance of the proposed approach. To mimic practical situations, UAV dynamics and sensor limitations were considered. The simulation results show that the proposed approach provides an efficient, reliable, and collision-free path without local minima and unreachable goal issues.

1. Introduction

Unmanned aerial vehicles (UAVs) have recently been used in several fields to perform various missions, such as search and rescue [1], surveillance [2], monitoring [3], delivery [4], and inspection [5]. UAVs are being increasingly used due to their cost efficiency, high reliability, and good performance. These needs lead to an increased demand for UAV airspace and an increased risk of collisions in that limited airspace. To perform successful missions and reduce the risk of collisions, the safe operation of UAVs must be ensured. That is, UAVs must equip a certain sensing system to detect threats and have a collision avoidance (CA) system to avoid threats without collision. This is called sense and avoid for UAVs [6].

Among several CA techniques, one of the most widely used algorithms is the artificial potential field (APF) method [7]. With its simple principle, the APF efficiently provides a smooth and safe trajectory. However, the APF faces some problems when one or multiple obstacles block an agent’s path toward the goal position with the symmetric pattern or obstacles are located near the goal position. The former is called local minima, and this generally occurs when obstacles are exactly located between the agent and the goal location. The latter is called the goal non-reachable obstacles nearby (GNRON) problem. Indeed, when the goal location is within the influence of the repulsive potential field, the agent faces a GNRON situation.

Several researchers have attempted to resolve these two issues by applying other techniques or modifying formulations. Some researchers have combined the APF with other existing algorithms. As an illustration, Amiryan et al. [8] used the rapidly-exploring random tree method as a global planner in the presence of static obstacles to generate a global path that does not encounter the local minima and GNRON problems and then applied the APF as a local planner to the global path generated beforehand. In Ref. [9], a bug algorithm was introduced to resolve the issues when the agent faces local minima or GNRON situation. Moreover, to deal with these issues, other scholars have suggested making modifications of the potential functions or fields of the APF. Triharminto et al. [10] modified the repulsive potential field formulation by adding a fractional term and another term with the signum function to solve the local minima and GNRON problems, respectively. In Ref. [11], the relative distance between an agent and a goal location was multiplied to the repulsive potential function to resolve the GNRON problem. Furthermore, the local minima problem was resolved by multiplying a term with exponential functions into the repulsive potential function. In Ref. [12], a fuzzy inference system based on the APF was considered to deal with the local minima problem. Furthermore, Zhang et al. [13] revised the repulsive potential field formulation by multiplying an additional term and combining it with the rolling window approach to resolve the APF’s two issues, and the intermediate target position was used to avoid obstacles. In addition, exponential functions [14] and Gaussian-like functions [15] have been utilized to formulate attractive and repulsive potential fields. Choi et al. [16] introduced the concept of the curl-free vector field to resolve the local minima issue. Although several researchers have proposed various techniques for resolving the APF’s issues, the aforementioned approaches have been applied in a two-dimensional (2D) environment and are not suitable in a situation where the 3D motion of UAVs has to be considered.

This research proposes a CA approach based on an enhanced potential field (EPF) building on the authors’ previous work developed in 2D environments to expand to become applicable to 3D environments by combining the horizontal and vertical avoidance maneuvers. Although the aforementioned studies undertaken by several researchers assumed that the obstacles’ position information is always given to an agent, this study used a 3D scanner to measure the obstacles’ point cloud information for practical sensing when the obstacles are within the sensing range and field of view (FOV). In addition, the dynamics of UAV are taken into account to describe the practical behaviors of the UAV. The proposed CA approach is validated at realistic cluttered urban simulation environments generated using real urban datasets.

This paper is organized as follows: Section 2 describes the dynamics and kinematic equations, the control logic, and the formulation of the proposed CA approach. Section 3 presents the analysis results for the key parameters of the proposed CA algorithm and the numerical simulation results to validate the performance of the proposed CA algorithm. Finally, the conclusion is described in Section 4.

2. Methodologies

2.1. Dynamic and Kinematic Equations

This work considers a quadcopter as an UAV. To describe the motion of a quadcopter, the reference frames and the vectors with the quadcopter model are displayed in Figure 1. The inertial frame is defined by ${\widehat{\mathbf{n}}}_k$ for k = 1, 2, and 3, and the body frame fixed on the center of mass of the quadcopter is defined by ${\widehat{\mathbf{b}}}_k$. The quadcopter’s position and velocity vectors expressed in the body frame are defined as $\mathbf{r}\in {\mathcal{R}}^3$ and $\mathbf{v}\in {\mathcal{R}}^3$, respectively. The quadcopter’s attitude is defined using the 3–2–1 set of the Euler angles, and the yaw, pitch, and roll angles of the quadcopter with respect to the inertial frame are defined as $\psi$, $\theta$, and $\varphi$, respectively. In addition, the quadcopter’s angular velocity vector with respect to the inertial frame is defined by $\mathbf{\omega }\in {\mathcal{R}}^3$. The direction cosine matrix from the inertial frame to the body frame of the quadcopter is defined as [17]

$$C_{BN}=\left[\begin{array}{ccc}cos\theta cos\psi & cos\theta sin\psi & -sin\theta \\ sin\varphi sin\theta cos\psi -cos\varphi sin\psi & sin\varphi sin\theta sin\psi +cos\varphi sin\psi & sin\varphi cos\theta \\ cos\varphi sin\theta cos\psi +sin\varphi sin\psi & cos\varphi sin\theta sin\psi -sin\varphi cos\psi & cos\varphi cos\theta \end{array}\right].$$

(1)

Using Equation (1), the position and velocity vector expressed in the inertial frame are defined as

$$\begin{array}{cc}\hfill {}^N\mathbf{r}& =C_{BN}^{\mathrm{T}}\mathbf{r},\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill {}^N\dot{\mathbf{r}}& =C_{BN}^{\mathrm{T}}\mathbf{v},\hfill \end{array}$$

(3)

where a vector with the superscript N denotes the vector expressed in the inertial frame. Note that the superscript B is omitted in the vectors with the body frame expression.

The motion of the quadcopter is composed of rotational and translational motions. The rotational motion of the quadcopter with respect to the inertial frame is defined as [17]

$$J\dot{\mathbf{\omega }}+\left[{\mathbf{\omega }}^{\times }\right]J\mathbf{\omega }={\mathbf{\tau }}_b-{\mathbf{\tau }}_g,$$

(4)

where $J\in {\mathcal{R}}^{3\times 3}$ is the inertia matrix, $\left[{\mathbf{\omega }}^{\times }\right]$ is the skew-symmetric matrix composed of the elements of the angular velocity vector (${\omega }_1$, ${\omega }_2$ and ${\omega }_3$), ${\mathbf{\tau }}_b\in {\mathcal{R}}^3$ is the external torque, and ${\mathbf{\tau }}_g=\left[{\mathbf{\omega }}^{\times }\right]{[0,0,I_r{\Omega }_r]}^{\mathrm{T}}$ is the gyroscopic torque generated by the quadcopter’s rotors. Here, $I_r$ is the inertia of the rotor, and ${\Omega }_r$ is the sum of the rotors’ speed, which is defined as ${\Omega }_r={\Omega }_1-{\Omega }_2+{\Omega }_3-{\Omega }_4$. Here, ${\Omega }_i$ for $i=$ 1, 2, 3, and 4 represents the rotation speed of the i-th rotor.

The quadcopter’s translational motion is defined as [18]

$$m\dot{\mathbf{v}}+\left[{\mathbf{\omega }}^{\times }\right]m\mathbf{v}={\mathbf{f}}_g+{\mathbf{f}}_b,$$

(5)

where m is the mass of the quadcopter and ${\mathbf{f}}_g\in {\mathcal{R}}^3$ is the gravitational force defined as

$${\mathbf{f}}_g=C_{BN}{}^N\mathbf{g},$$

(6)

where ${}^N\mathbf{g}={[0,0,-mg]}^{\mathrm{T}}$ is the gravitational force vector expressed in the inertial frame with the gravitational constant g. Furthermore, the thrust force vector is defined by ${\mathbf{f}}_b={[0,0,f_{th}]}^{\mathrm{T}}$, and the thrust force along with the ${\widehat{\mathbf{b}}}_3$ axis generated by the quadcopter’s rotors is expressed as

$$f_{th}=k_f({\Omega }_1^2+{\Omega }_2^2+{\Omega }_3^2+{\Omega }_4^2),$$

(7)

where $k_f$ is the aerodynamic force constant.

In addition, the control input vector for controlling the quadcopter is defined as

$$\mathbf{u}=\left[\begin{array}{c}{f}_{th}\\ {\tau }_1\\ {\tau }_2\\ {\tau }_3\end{array}\right]=\left[\begin{array}{cccc}{k}_f& k_f& k_f& k_f\\ -l{k}_f& -l{k}_f& l{k}_f& l{k}_f\\ -l{k}_f& l{k}_f& l{k}_f& -l{k}_f\\ -k_{\tau }& k_{\tau }& -k_{\tau }& k_{\tau }\end{array}\right]\left[\begin{array}{c}{\Omega }_1^2\\ {\Omega }_2^2\\ {\Omega }_3^2\\ {\Omega }_4^2\end{array}\right],$$

(8)

where ${\tau }_k$ for $k=$ 1, 2, and 3 is the body axes’ torque generated by the rotors’ speed difference, $k_{\tau }$ is the aerodynamic torque constant, and l is the distance between the center of the quadcopter and the center position of propellers 1 and 2, as shown in Figure 1.

The kinematic differential equation composed of the Euler angles and the angular velocity is defined as [17]

$$\left[\begin{array}{c}\dot{\psi }\\ \dot{\theta }\\ \dot{\varphi }\end{array}\right]={\displaystyle \frac{1}{cos\theta }}\left[\begin{array}{ccc}0& sin\varphi & cos\varphi \\ 0& cos\varphi cos\theta & -sin\varphi cos\theta \\ cos\theta & sin\varphi sin\theta & cos\varphi sin\theta \end{array}\right]\left[\begin{array}{c}{\omega }_1\\ {\omega }_2\\ {\omega }_3\end{array}\right].$$

(9)

2.2. Quadcopter Control Logic

To control the position of the quadcopter with respect to the inertial frame, the translational motion of the quadcopter expressed in the inertial frame is derived as [19]

$$m{}^N\ddot{\mathbf{r}}={}^N\mathbf{g}+C_{BN}^{\mathrm{T}}{\mathbf{f}}_b,$$

(10)

and ${}^N\ddot{\mathbf{r}}$ in Equation (10) is rewritten as

$${}^N\ddot{\mathbf{r}}=\left[\begin{array}{c}\ddot{x}\\ \ddot{y}\\ \ddot{z}\end{array}\right]=\left[\begin{array}{c}{\displaystyle \frac{f_{th}}{m}}(sin\varphi sin\psi +cos\varphi cos\psi sin\theta )\\ {\displaystyle \frac{f_{th}}{m}}(-sin\varphi cos\psi +cos\varphi sin\psi sin\theta )\\ -g+{\displaystyle \frac{f_{th}}{m}}(cos\varphi cos\theta )\end{array}\right].$$

(11)

Since positions x and y are coupled with the attitude of the quadcopter, they cannot be directly controlled using thrust force. Therefore, Equation (10) is rewritten as

$$\begin{array}{cc}\hfill {\ddot{x}}_d& ={\displaystyle \frac{f_{th}}{m}}(sin{\varphi }_{d}sin\psi +cos{\varphi }_{d}cos\psi sin{\theta }_d),\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill {\ddot{y}}_d& ={\displaystyle \frac{f_{th}}{m}}(-sin{\varphi }_{d}cos\psi +cos{\varphi }_{d}sin\psi sin{\theta }_d),\hfill \end{array}$$

(13)

where ${\ddot{x}}_d$ and ${\ddot{y}}_d$ are the desired accelerations with regard to x and y, and ${\varphi }_d$, ${\theta }_d$, and ${\psi }_d$ are the desired roll, pitch, and yaw angles, respectively. To compute the desired attitude angles, small angle assumptions were considered in this work ($sin{\varphi }_d\approx {\varphi }_d$, $sin{\theta }_d\approx {\theta }_d$, $cos{\varphi }_d\approx 1,cos{\theta }_d\approx 1$). This assumption is generally valid except for the large changes of the roll and pitch angles duirng a short time period. Hence, Equations (12) and (13) are rewritten as

$$\begin{array}{cc}\hfill {\ddot{x}}_d& ={\displaystyle \frac{f_{th}}{m}}({\varphi }_{d}sin\psi +{\theta }_{d}cos\psi ),\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill {\ddot{y}}_d& ={\displaystyle \frac{f_{th}}{m}}(-{\varphi }_{d}cos\psi +{\theta }_{d}sin\psi ).\hfill \end{array}$$

(15)

The desired x and y accelerations are obtained by the PID control as follows:

$$\begin{array}{cc}\hfill {\ddot{x}}_d& =k_{P,pos}(x_d-x)+k_{I,pos}\int (x_d-x)\mathrm{d}t+k_{D,pos}({\dot{x}}_d-\dot{x}),\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill {\ddot{y}}_d& =k_{P,pos}(y_d-y)+k_{I,pos}\int \left(y_d-y\right)\mathrm{d}t+k_{D,pos}({\dot{y}}_d-\dot{y}),\hfill \end{array}$$

(17)

where $k_{P,pos}$, $k_{I,pos}$, and $k_{D,pos}$ are the proportional, integral, and derivative control gains for position control, respectively. Moreover, for the altitude control of the quadcopter, the thrust force is obtained by the PID control as follows:

$$f_{th}=k_{P,pos}(z_d-z)+k_{I,pos}\int (z_d-z)\mathrm{d}t+k_{D,pos}({\dot{z}}_d-\dot{z}).$$

(18)

After obtaining the desired accelerations, the desired attitude is obtained using Equations (14) and (15). Then, a PD controller is developed to control the attitude of the quadcopter. The control torque for attitude control is derived as

$$\begin{array}{cc}\hfill {\tau }_1& =k_{P,att}({\varphi }_d-\varphi )+k_{D,att}({\dot{\varphi }}_d-\dot{\varphi }),\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill {\tau }_2& =k_{P,att}({\theta }_d-\theta )+k_{D,att}({\dot{\theta }}_d-\dot{\theta }),\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill {\tau }_3& =k_{P,att}({\psi }_d-\psi )+k_{D,att}({\dot{\psi }}_d-\dot{\psi }),\hfill \end{array}$$

(21)

where $k_{P,att}$ and $k_{D,att}$ are the proportional and derivative control gains for attitude control, respectively.

2.3. Collision Avoidance Algorithm

The widely used CA algorithm is known as the APF due to its simple principle. The APF consists of attractive and repulsive potential fields. For example, a goal position attracts a quadcopter in the attractive potential field while obstacles repel the quadcopter in the repulsive potential field. As mentioned in Section 1, however, the APF has GNRON and local minima issues. The GNRON issue is resolved by introducing an additional term in the repulsive potential function. The total potential function is defined as [20,21]

$$U_t{(}^N\mathbf{r})=U_a{(}^N\mathbf{r})+U_r{(}^N\mathbf{r}),$$

(22)

where

$$\begin{array}{cc}\hfill U_a{(}^N\mathbf{r})& =\frac{1}{2}{k}_{a}d{{(}^N\mathbf{r},{}^N{\mathbf{r}}_g)}^2,\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill U_r{(}^N\mathbf{r})& =\left\{\begin{array}{cc}{\displaystyle \frac{1}{2}{k}_r{\left(\frac{1}{d{(}^N\mathbf{r},{}^N{\mathbf{r}}_o)}-\frac{1}{d_o}\right)}^{2}d{{(}^N\mathbf{r},{}^N{\mathbf{r}}_g)}^{n_g}}\hfill & \mathrm{if}d{(}^N\mathbf{r},{}^N{\mathbf{r}}_o)\le d_o,\hfill \\ 0\hfill & \mathrm{if}d{(}^N\mathbf{r},{}^N{\mathbf{r}}_o)>d_o.\hfill \end{array}\right.\hfill \end{array}$$

(24)

Here, $U_a{(}^N\mathbf{r})$ is the attractive potential function, $U_r{(}^N\mathbf{r})$ is the repulsive potential function, ${}^N\mathbf{r}\in {\mathcal{R}}^3$ is the quadcopter position in the inertial frame, ${}^N{\mathbf{r}}_g\in {\mathcal{R}}^3$ is the goal position, ${}^N{\mathbf{r}}_o\in {\mathcal{R}}^3$ is the obstacle’s position, $d(\mathbf{a},\mathbf{b})$ is the relative distance between generic vectors $\mathbf{a}\in {\mathcal{R}}^3$ and $\mathbf{b}\in {\mathcal{R}}^3$, $d_o$ is the distance influenced by the repulsive potential function, and $n_g$ is the order of the relative distance term. Furthermore, $k_a$ and $k_r$ are the gain coefficients of attractive and repulsive potential functions, respectively. Due to $d{{(}^N\mathbf{r},{}^N{\mathbf{r}}_g)}^{n_g}$, the total potential function solely becomes zero when the quadcopter arrives at the goal position. This is how the GNRON issue is resolved. Since the potential functions with $n_g=0$ are the same as the original one presented by Khatib [7], $n_g$ should be greater than 0. By taking the negative gradient to the potential functions, the potential fields, which are the vector form, are derived as

$$\begin{array}{cc}\hfill {\mathbf{f}}_a{(}^N\mathbf{r})& =-\nabla U_a{(}^N\mathbf{r})=k_{a}d{(}^N\mathbf{r},{}^N{\mathbf{r}}_g)\frac{\partial d{(}^N\mathbf{r},{}^N{\mathbf{r}}_g)}{\partial {}^N\mathbf{r}},\hfill \end{array}$$

(25)

$$\begin{array}{cc}\hfill {\mathbf{f}}_r{(}^N\mathbf{r})& =-\nabla U_r{(}^N\mathbf{r})=\left\{\begin{array}{cc}{\mathbf{f}}_{r,o}{(}^N\mathbf{r})+{\mathbf{f}}_{r,g}{(}^N\mathbf{r})\hfill & \mathrm{if}d{(}^N\mathbf{r},{}^N{\mathbf{r}}_o)\le d_o,\hfill \\ 0\hfill & \mathrm{if}d{(}^N\mathbf{r},{}^N{\mathbf{r}}_o)>d_o,\hfill \end{array}\right.\hfill \end{array}$$

(26)

where

$$\begin{array}{cc}\hfill {\mathbf{f}}_{r,o}{(}^N\mathbf{r})& =-k_r\frac{d{{(}^N\mathbf{r},{}^N{\mathbf{r}}_g)}^{n_g}}{d{{(}^N\mathbf{r},{}^N{\mathbf{r}}_o)}^2}\left(\frac{1}{d{(}^N\mathbf{r},{}^N{\mathbf{r}}_o)}-\frac{1}{d_o}\right)\frac{\partial d{(}^N\mathbf{r},{}^N{\mathbf{r}}_o)}{\partial {}^N\mathbf{r}},\hfill \end{array}$$

(27)

$$\begin{array}{cc}\hfill {\mathbf{f}}_{r,g}{(}^N\mathbf{r})& =\frac{1}{2}{n}_g{k}_{r}d{{(}^N\mathbf{r},{}^N{\mathbf{r}}_g)}^{n_g-1}{\left(\frac{1}{d{(}^N\mathbf{r},{}^N{\mathbf{r}}_o)}-\frac{1}{d_o}\right)}^2\frac{\partial d{(}^N\mathbf{r},{}^N{\mathbf{r}}_g)}{\partial {}^N\mathbf{r}}.\hfill \end{array}$$

(28)

Note that ${\mathbf{f}}_a{(}^N\mathbf{r})\in {\mathcal{R}}^3$ is the attractive potential field and ${\mathbf{f}}_r{(}^N\mathbf{r})\in {\mathcal{R}}^3$ is the repulsive potential field. Also, $\partial d{(}^N\mathbf{r},{}^N{\mathbf{r}}_o)/\partial {}^N\mathbf{r}$ is the direction from the quadcopter to the obstacle, and $\partial d{(}^N\mathbf{r},{}^N{\mathbf{r}}_g)/\partial {}^N\mathbf{r}$ is the direction from the quadcopter to the goal position. When the goal position is under the influence of the repelling force, the quadcopter can arrive at the goal position because the repulsive potential field also contains ${\mathbf{f}}_{r,g}$ which produces the potential field toward the goal position. That is, this term contributes to resolving the GNRON issue. By summing up two potential fields, the total potential field (${\mathbf{f}}_t$) is computed.

Since Equations (25)–(28) only resolve the GNRON issue, the local minima issue still remains. Hence, this study focuses on resolving the local minima issue by redefining the first-term ${\mathbf{f}}_{r,o}{(}^N\mathbf{r})$ of the repulsive potential field in Equation (27) as follows:

$${\mathbf{f}}_{r,e}{(}^N\mathbf{r})=-k_r\frac{d{{(}^N\mathbf{r},{}^N{\mathbf{r}}_g)}^{n_g}}{d{{(}^N\mathbf{r},{}^N{\mathbf{r}}_o)}^2}\left(\frac{1}{d{(}^N\mathbf{r},{}^N{\mathbf{r}}_o)}-\frac{1}{d_o}\right)\widehat{\mathbf{q}}.$$

(29)

It is important to note that the direction of ${\mathbf{f}}_{r,o}{(}^N\mathbf{r})$ is replaced with $\widehat{\mathbf{q}}\in {\mathcal{R}}^3$ because the local minima issue is highly related to the determination of the direction in front of the obstacles, as shown in Figure 2. This shows the difference between the generated total potential field of the conventional and proposed approaches in the local minima situation. While the conventional one has a total potential field of zero (i.e., ${\mathbf{f}}_a+{\mathbf{f}}_{r,o}+{\mathbf{f}}_{r,g}=\mathbf{0}$), the proposed approach results in the non-zero total potential field by utilizing ${\mathbf{f}}_{r,e}$ instead of ${\mathbf{f}}_{r,o}$.

The new direction vector $\widehat{\mathbf{q}}$ is defined by combining the horizontal avoidance direction with the vertical one. The horizontal and vertical avoidance direction vectors (i.e., ${\mathbf{q}}_h={[q_{h,1},q_{h,2},q_{h,3}]}^{\mathrm{T}}\in {\mathcal{R}}^3$ and ${\mathbf{q}}_v={[q_{v,1},q_{v,2},q_{v,3}]}^{\mathrm{T}}\in {\mathcal{R}}^3$) are defined as

$$\begin{array}{cc}\hfill {\mathbf{q}}_h& =R_h({{}^N\mathbf{r}}_o-{}^N\mathbf{r}),\hfill \end{array}$$

(30)

$$\begin{array}{cc}\hfill {\mathbf{q}}_v& =R_v({{}^N\mathbf{r}}_o-{}^N\mathbf{r}),\hfill \end{array}$$

(31)

where $R_h\in {\mathcal{R}}^{3\times 3}$ and $R_v\in {\mathcal{R}}^{3\times 3}$ are the rotation matrices that rotate the direction of the conventional repulsive potential field into a horizontal new direction and a vertical new direction, respectively. These new directions seek to avoid the situation in which the quadcopter becomes stuck in front of the obstacle. Once avoidance directions are determined for each plane, $\widehat{\mathbf{q}}$ is calculated as follows:

$$\widehat{\mathbf{q}}=\frac{{\mathbf{q}}^{\prime }}{\left|\right|{\mathbf{q}}^{\prime }\left|\right|},$$

(32)

where

$${\mathbf{q}}^{\prime }=\left[\begin{array}{c}\alpha q_{h,1}\\ \alpha q_{h,2}\\ (1-\alpha )q_{v,3}\end{array}\right].$$

(33)

Here, $\alpha \in [0,1]$ is a weighting factor for consolidating two avoidance directions. Note that $q_{h,1}$ and $q_{h,2}$ for the horizontal direction and $q_{v,3}$ for the vertical direction are only used because these are the dominant terms for generating the maneuvers in each direction.

To determine the proper rotation matrices for the collision-free repelling direction of the repulsive potential field, this work considers the relative angles between the quadcopter’s velocity vectors and the relative position vector between the quadcopter and the obstacle in the horizontal and vertical planes. To obtain the relative angles, an additional reference frame (i.e., E frame) is defined using ${\widehat{\mathbf{e}}}_k$ for $k=$ 1, 2, and 3, as shown in Figure 3a. Note that ${\widehat{\mathbf{e}}}_1$ is a unit vector that has the same direction along with the quadcopter’s velocity direction, and ${\widehat{\mathbf{e}}}_2$ is defined as ${\widehat{\mathbf{n}}}_3\times {\widehat{\mathbf{e}}}_1$. Then, ${\widehat{\mathbf{e}}}_3$ is determined by ${\widehat{\mathbf{e}}}_1\times {\widehat{\mathbf{e}}}_2$. Here, the E frame is not defined during take-off and landing because the CA algorithm is not activated during these phases. Then, the relative angle ${\rho }_h$ defined in the horizontal plane is shown in Figure 3b, and ${\rho }_v$ is defined in the vertical plane, as displayed in Figure 3c. Note that ${\rho }_h$ and ${\rho }_v$ provide the relative location of the obstacle from the quadcopter with respect to the horizontal and vertical planes, respectively. Based on this information, the rotation matrices to determine the avoidance direction are obtained by the criteria listed in

Table 1. Criteria for rotation matrices to determine the avoidance direction.

	Relative Angle	Rotation Matrix	Remark
Horizontal plane	${\rho }_h\ge 0$	$R_h=R_3$	$R_3\equiv \left[\begin{array}{ccc}cos\gamma & -sin\gamma & 0\\ sin\gamma & cos\gamma & 0\\ 0& 0& 1\end{array}\right]$
	${\rho }_h<0$	$R_h=R_3^{\mathrm{T}}$
Vertical plane	${\rho }_v\ge 0$	$R_v=R_2$	$R_2\equiv \left[\begin{array}{ccc}cos\gamma & 0& sin\gamma \\ 0& 1& 0\\ -sin\gamma & 0& cos\gamma \end{array}\right]$
	${\rho }_v<0$	$R_v=R_2^{\mathrm{T}}$

.

In the formulation of the proposed algorithm, there are several free parameters that need to be determined by a user. In particular, $\alpha$ and $\gamma$, which are the newly introduced parameters in this paper, are critical because they highly affect the performance of the proposed approach. Here, $\alpha$ plays a role in determining the portion of each avoidance direction. For example, when $\alpha =1$, the quadcopter only performs a horizontal avoidance maneuver, and when $\alpha =0$, it only flies on a vertical plane. When the quadcopter’s operation altitude is less than the height of buildings (i.e., obstacles to the quadcopter) in urban areas, the horizontal avoidance maneuver is sufficient. However, when the quadcopter is operating near building tops, a combination of the horizontal and vertical maneuvers is preferable since the horizontal maneuver only may require a longer trajectory than the combined maneuvers. Therefore, $\alpha$ needs to be properly determined between 0 and 1. In addition to this, $\gamma$ directly affects the direction of ${\mathbf{f}}_{r,e}$, as shown in Figure 2. When $\gamma =0$ deg, the proposed algorithm is performed as the same with the conventional approach. When $\gamma$ is greater than 90 deg, ${\mathbf{f}}_{r,e}$ (i.e., the repelling direction) is in the direction of the obstacle. This means that ${\mathbf{f}}_{r,e}$ does not play a role in repelling the quadcopter anymore. Hence, $\gamma$ is needed to be selected between 0 and 90 deg. These two parameters will be analyzed in the next section.

3. Simulation Studies

3.1. Analysis and Selection of the Key Parameters

Among several parameters of the proposed CA approach in this paper, the effects of $\alpha$ and $\gamma$ were analyzed through numerical simulations to select proper values. To plainly evaluate the effect of each parameter independently, the quadcopter’s dynamics are not considered in this subsection. In other words, the inclusion of dynamics makes it difficult to analyze the CA trajectories purely affected by $\alpha$ and $\gamma$ because the quadcopter’s motion is being affected by the controller and the quadcopter’s parameters. Simulation studies are performed with several values of $\alpha$ and $\gamma$ within the ranges mentioned in the previous section. The other simulation parameters for the proposed algorithm are listed in

Table 2. Simulation parameters for the proposed CA approach.

Parameter	Value
Time interval, $\mathrm{d}t$	0.1 (s)
Attractive gain coefficient, $k_a$	0.01
Repulsive gain coefficient, $k_r$	1
Order of relative distance term, $n_g$	2
Obstacle threshold, $d_o$	10 (m)
Physical collision	2 (m)

.

To analyze the effect of $\gamma$’s changes, only the horizontal avoidance maneuver (i.e., $\alpha =0$) was considered. The local minima situation was selected to test, and in this situation, the conventional APF approach cannot find the proper path, so the $\gamma =0$ case does not need to be tested. Figure 4 shows the trajectories of the quadcopter with various $\gamma$ values, and Figure 5 displays the total path lengths of the quadcopter and the minimum relative distances between the quadcopter and the nearest obstacle in accordance with $\gamma$. As $\gamma$ increases, the quadcopter’s trajectory tends to approach the obstacles, and this result can be confirmed in Figure 5. The shorter the total path length, the lesser the energy consumption will be, by assuming the constant velocity-based operation for example. However, since too close a relative distance may cause a collision, a safety perspective must also be considered when selecting a $\gamma$ value. This study assumed that the quadcopter is safe when the relative distance to the nearest obstacle is greater than 8 m. The minimum relative distance (i.e., safety margin) is highlighted with the red-dotted line in Figure 5. Here, the quadcopter may approach the obstacle closer than the safety margin if the dynamics of the quadcopter is considered. Furthermore, since the safety margin is 8 m, the quadcopter may collide with the nearest obstacle in 3 s when the quadcopter flies at a nominal velocity of 3 m/s. Based on this setting and information, $\gamma =45$ deg is the best selection in this study, and the selection of $\gamma =45$ deg results in the minimum relative distance of 8.3 m.

With the fixed $\gamma =45$ deg, the next analysis is that of the effect of $\alpha$’s changes. To analyze the vertical avoidance maneuver, the operational altitude of the quadcopter in this simulation is selected around the top height of the building. Similar to the analysis process for $\gamma$, Figure 6 describes the trajectories of the quadcopter with various $\alpha$ values, and Figure 7 shows the total path lengths of the quadcopter and the minimum relative distances between the quadcopter and the nearest obstacle in accordance with $\alpha$. As $\alpha$ changes, the degree of the portion of horizontal and vertical maneuvers on the planned trajectory also changes. Since the quadcopter operates near the top of the obstacle, the vertical avoidance maneuver (i.e., $\alpha =0$) provides better performance than the others in terms of path length. However, when $\alpha =0.5$, the result is best, as shown in Figure 7, because the minimum relative distance it is required to maintain is greater than the safety margin (8 m).

3.2. Numerical Simulation Results in an Urban Environment

To test the proposed CA approach in realistic environments, a high-resolution urban area was simulated using real datasets measured by an airborne LiDAR. Since the open source urban data are available, this work utilized data obtained from Open-Topography (https://opentopography.org/, accessed on 17 October 2021). Note that LiDAR data (i.e., “*.las” file here) of San Diego downtown were selected as an example of 3D urban environments. The data are classified into x, y, and z components of a point cloud, and one filters out the data so that all of z is less than 200 m. The filtered point cloud data are visualized in Figure 8.

To validate the CA performance of the proposed approach, the numerical simulations of three scenarios in the selected urban environment were performed. The local minima and GNRON cases were chosen as the first and second scenarios, respectively, to highlight the fact that the proposed approach resolves the conventional APF’s well-known issues. The third case is a complex urban environment in the presence of multiple static obstacles (e.g., buildings). In addition, a quadcopter in this study is mounted on a 3D scanner (e.g., depth camera) to simulate practical sensing. By the 3D scanner, the point cloud of objects (i.e., obstacles here) existing within the sensing range and FOV of the sensor are measured. The simulation parameters for the quadcopter, sensor’s specifications, and control gains are listed in Table 3. Here, the quadcopter’s dynamics and control logic are considered, unlike in Section 3.1. Furthermore, the quadcopter’s diameter size was chosen based on information about the commercial drone’s specifications [22]. Moreover, this work assumed that the quadcopter maintains a constant nominal speed during operation. In addition, the parameters for the proposed CA algorithm listed in Table 2, and $\alpha$ and $\gamma$ values selected in Section 3.1 were used. For each case, the start and goal positions are tabulated in Table 4.

Table 3. Simulation parameters for the quadcopter system.

Parameter		Value
Quacopter’s diameter	$d_q$	0.6 (m)
Quacopter’s nominal speed	$v_n$	3 (m/s)
Quacopter’s mass	m	0.65 (kg)
Quadcopter inertia	J	diag($[7.5,7.5,13.0]$) $\times {10}^{-3}$ (kg m${}^2$)
Gravity constant	g	9.81 (m/s${}^2$)
Rotor inertia	I	$6\times {10}^{-5}$ (kg m${}^2$)
Distance from center to rotors 1 and 2	l	0.23 (m)
Aerodynamic force constant	$k_f$	$3.23\times {10}^{-5}$ (N s${}^2$)
Aerodynamic torque constant	$k_{\tau }$	$7.5\times {10}^{-5}$ (N ms${}^2$)
Sensor’s sensing range	$d_s$	20 (m)
Sensor’s horizontal FOV	${\kappa }_h$	210 (deg)
Sensor’s vertical FOV	${\kappa }_v$	70 (deg)
Position control gains	$k_{P,pos}$	2
	$k_{I,pos}$	0.001
	$k_{D,pos}$	1
Attitude control gains	$k_{P,att}$	0.28
	$k_{D,att}$	0.05

Table 3. Simulation parameters for the quadcopter system.

Parameter		Value
Quacopter’s diameter	$d_q$	0.6 (m)
Quacopter’s nominal speed	$v_n$	3 (m/s)
Quacopter’s mass	m	0.65 (kg)
Quadcopter inertia	J	diag($[7.5,7.5,13.0]$) $\times {10}^{-3}$ (kg m${}^2$)
Gravity constant	g	9.81 (m/s${}^2$)
Rotor inertia	I	$6\times {10}^{-5}$ (kg m${}^2$)
Distance from center to rotors 1 and 2	l	0.23 (m)
Aerodynamic force constant	$k_f$	$3.23\times {10}^{-5}$ (N s${}^2$)
Aerodynamic torque constant	$k_{\tau }$	$7.5\times {10}^{-5}$ (N ms${}^2$)
Sensor’s sensing range	$d_s$	20 (m)
Sensor’s horizontal FOV	${\kappa }_h$	210 (deg)
Sensor’s vertical FOV	${\kappa }_v$	70 (deg)
Position control gains	$k_{P,pos}$	2
	$k_{I,pos}$	0.001
	$k_{D,pos}$	1
Attitude control gains	$k_{P,att}$	0.28
	$k_{D,att}$	0.05

Table 4. Start and goal positions of the quadcopter for each case.

	Description	Start Position (m)	Goal Position (m)
Case 1	Local minima	${[360,130,60]}^{\mathrm{T}}$	${[250,130,60]}^{\mathrm{T}}$
Case 2	GNRON	${[380,110,70]}^{\mathrm{T}}$	${[350,208,60]}^{\mathrm{T}}$
Case 3	Complex environment	${[370,190,78]}^{\mathrm{T}}$	${[200,50,115]}^{\mathrm{T}}$

Table 4. Start and goal positions of the quadcopter for each case.

	Description	Start Position (m)	Goal Position (m)
Case 1	Local minima	${[360,130,60]}^{\mathrm{T}}$	${[250,130,60]}^{\mathrm{T}}$
Case 2	GNRON	${[380,110,70]}^{\mathrm{T}}$	${[350,208,60]}^{\mathrm{T}}$
Case 3	Complex environment	${[370,190,78]}^{\mathrm{T}}$	${[200,50,115]}^{\mathrm{T}}$

In Case 1, an obstacle point is located exactly between the start and goal positions. Figure 9 shows the quadcopter’s avoidance trajectory and the nearest obstacle point sensed by the 3D scanner at each simulation discrete-time. Note that the blue trajectory represents the desired trajectory obtained from the proposed CA algorithm, while the black one indicates the quadcopter’s actual trajectory. At every time step, the CA algorithm computes the desired position for avoiding obstacles, and the control logic mentioned in Section 2.2 uses the desired position calculated. Furthermore, Figure 10 depicts the relative distance between the quadcopter and the nearest obstacle point. In the trajectory figure, it is observed that there are slight changes in the quadcopter’s altitude near the obstacle, and this is because of the portion of the vertical avoidance direction inclusion. In addition, the relative distance during the maneuver remains greater than the collision threshold 2 m that represents physical collision. Note that the obstacle threshold of 10 m denotes the relative distance that activates the repulsive potential field, as defined in Table 2. Since the smallest relative distance is observed at 8.42 m, one can say that the quadcopter completes the mission safely.

In Case 2, to construct the GNRON case, the goal position is placed near the obstacle. The quadcopter’s CA trajectory by the proposed approach is described in Figure 11. The quadcopter completes the mission without any collision, although the goal position is within the range (10 m) in which the repulsive potential field is activated. That is, although the repelling force is produced near the goal position by the repulsive potential field, the quadcopter smoothly arrives at the goal position. Unlike in Case 1, a few altitude changes were observed due to the smooth direction changes required. Figure 12 shows the relative distance to the nearest obstacle. The relative distance continuously decreases within the obstacle threshold until the quadcopter reaches the goal position. Once the quadcopter arrives at the goal position, the relative distance is not measured to perform the landing. Thus, there is no time history after 70 s. The closest distance 7.62 m was measured at the goal position, and the quadcopter complied with the collision threshold of 2 m.

The last case is that of a complex environment with several static obstacles. In this scenario, the quadcopter must avoid three buildings to achieve the mission. Figure 13 shows the quadcopter’s 3D and 2D trajectories. The quadcopter mainly performs the horizontal avoidance maneuver near the first and third buildings and the vertical avoidance maneuver for the second building. Since the proposed approach combines the horizontal and vertical maneuvers, the quadcopter can successfully avoid several obstacles without collision. Figure 14 illustrates the relative distance from the quadcopter to the nearest obstacle between the sensing range of 20 m and the collision threshold of 2 m. The smallest relative distance can be observed at 8.5 m so that the quadcopter safely reaches the goal position without collision—even in the presence of multiple static obstacles within the limited 3D urban environment.

4. Conclusions

This paper proposes a novel collision avoidance (CA) approach that reformulates the repulsive potential function of the artificial potential field. The new formulation of the proposed CA approach was mainly developed in two-dimensional (2D) planes, i.e., the horizontal and vertical planes, by introducing the rotation matrices that resolve the local minima issue. Then, the avoidance maneuvers for each 2D plane were combined using a weighting factor to avoid obstacles in cluttered 3D urban environments. Numerical simulations were executed to validate the performance of the proposed CA approach. Urban LiDAR data were utilized for creating an urban environment, and 3D scanner’s characteristics and quadcopter’s dynamics were considered for simulating the practical sensing and motion. The resulting quadcopter’s trajectory and relative distance between the quadcopter and the nearest obstacle evidently show that the quadcopter reaches the goal position without collision for several of the meaningful scenarios. In the future, more realistic and practical environments, such as those involving the presence of dynamic obstacles and the 3D simulator that utilizes high-performance physics engines (e.g., Gazebo), will be considered to further improve the proposed approach.