Autonomous drone delivery is an emerging trend in modern aviation, with many companies already testing the commercial service in cities. Even though safety precautions are required, an in-flight failure is not fully preventable and there always will be some probability of a crash after an in-flight failure.

The potential risk can be minimized by risk-aware trajectory planning, in other words, by planning the flight path such that if an in-flight failure occurs, the consequences will be minimal. Such risk-aware cost function can be employed into a trajectory planner such as RRT*, which you have implemented in the last assignment. The aim of this assignment is to implement a risk cost function assessing a risk associated with an in-flight loss of thrust on a small quad-copter.

The risk is defined as the number of possible on-ground casualties, inspired by [1]. The risk $M(x, E, \gamma)$ associated with an impact location $x$, impact energy $E$, and impact angle $\gamma$ is given as $$ M(x, E, \gamma) = p_{hit}(x, \gamma) p_{casualty}(x, E)\,.$$ The probability of hitting a person is given as $$p_{hit}(x, \gamma) = \varrho(x) A_{exp}(\gamma)\,,$$ where $\varrho(x)$ denotes the population density at $x$ and $A_{exp}$ is the area exposed to the crash expressed as $$A_{exp}(\gamma) = 2(r_p + r_{uav}) \frac{h_p}{\tan |\gamma|} + \pi (r_p + r_{uav})^2\,.$$ The $r_{uav}$ denotes the circumsribed circle of the UAV, $r_p = 0.2$ m and $h_p = 1.8$ m denotes the radius and height of an average person.

The probability of causing a casualty if person is hit, is formally defined as $$p_{casualty}(x, E) = \frac{1-k}{1-2k+\sqrt{\frac{\alpha}{\beta}} \left(\frac{\beta}{E}\right)^\frac{3}{S(x)}}\,,$$ where $k = \mathrm{min} \left(1, \left(\frac{\beta}{E}\right)^\frac{3}{S(x)} \right)$, $S(x)$ denotes a sheltering factor at impact location $x$ (how much is the environment protecting people), $\alpha$ is the impact energy needed to achieve $p_{casualty} = 50\%$ when $S = 6$, and $\beta$ is the impact energy causing a casualty for $S \rightarrow 0$; $\beta = 34$ J is used based on [2].

The ground map data consist of three layers – height, sheltering factor, and population density. The map data are taken from OpenStreetMaps [3], terrain profile is taken from [4], and population density is from [5]. The sheltering factors are arbitrarily assigned based on the object type. The data are stored in form of PNG image. A scaled version for visualization is also provided in the resource pack.

The herein considered possible in-flight failure is a loss of thrust, which leads to a ballistic fall on a small quad-copter. The fall is described by an ODE (source) \begin{equation} \frac{\mathrm{d}}{\mathrm{d}t} \left(\begin{matrix} x \\ y \\v_x \\ v_y \end{matrix}\right) = \left(\begin{matrix} v_x \\ v_y \\ -\mu v_x \sqrt{v_x^2 + v_y^2} \\ -g-\mu v_y \sqrt{v_x^2 + v_y^2}\end{matrix} \right)\,,\end{equation} where $\mu = \frac{C\rho S}{2 m}$, $C$ is the drag coefficient, $\rho$ is the air density, $S$ is the cross-section area of the UAV, and $m$ is its mass. This equation does not have a closed-form solution, and has to be solved numerically.

The parameters of a ballistic fall may be uncertain, especially the drag coefficient and the initial heading. Thus, the predicted impact location is not an exact place, but a probability map. The uncertain parameters are assumed to be of Gaussian distribution. Random samples are taken and a solution to the ODE is found; the impact location is stored. The process is repeated until the given amount of falls have been solved, and the impact locations are fitted with a multivariate Gaussian (in heading x lateral distance coordinates). Finally, the Gaussian can be sampled to obtain a grid-based impact probability map.

The risk $\mathcal{R}(q)$ associated with a failure occurring at configuration $q$ is simply given as \begin{equation}\mathcal{R}(q) = \sum\limits_{\mathbb{R}^2} p_{imp}(x | \Gamma_f) M(x, E, \gamma)\,,\end{equation} where $\Gamma_f$ is the most likely trajectory of the drone after the failure occurred. The impact probability map is independent on the configuration where the failure occurred; it depends only on the fallen altitude (and the parameters, but those are considered fixed for a single scenario). To increase the computational efficiency, the impact probability maps can be precomputed for various fallen altitudes and initial headings.

To estimate a risk $\mathcal{R}(q)$ associated with a failure occurring at $q$, the most likely fall is used to determine the predicted impact location on the terrain, and the corresponding fallen altitude. Based on the fallen altitude and initial heading, the most suitable impact probability map from the precomputed ones is selected, and the risk is calculated according to (2). The risk associated with a trajectory $\Gamma: 0 \rightarrow T$ is then simply given as \begin{equation}\mathcal{R}(\Gamma) = \int_0^T \mathcal{R}(\Gamma(t)) \mathrm{d}t\,,\end{equation} where $\Gamma(t)$ stands for a configuration $q$ at specified distance along the trajectory.

[1] S. Primatesta, A. Rizzo, A. la Cour-Harbo, Ground risk map for unmanned aircraft in urban environments, Journal of Intelligent & Robotic Systems 97 (3) (2020) 489–509.
[2] R. Standard, 321-07, common risk criteria for national test ranges. range safety group risk committee, range commanders council, US Army White Sands Missile Range, New Mexico. [3] OpenStreetMap contributors, Planet dump retrieved from https://planet.osm.org, https://www.openstreetmap.org, accessed on: 28 Feb 2023 (2023).
[4] NASA, Nasa shuttle radar topography mission global 1 arc second, https://opentopography.org/, accessed on: 28 Feb 2023 (2013).
[5] Facebook Connectivity Lab and Center for International Earth Science Information Network - CIESIN - Columbia University, High Resolution Settlement Layer (HRSL). Source imagery for HRSL © 2016 DigitalGlobe, accessed on: 28 Feb 2023 (2016).

Implement evaluating a risk associated with potential in-flight loss of thrust on a small quad-copter according to the definitions above.

Implement creation of ground map layers in class GroundMap. The map layers are provided in 16-bit PNG format, where the values of each pixel represent the particular value. Note that the loaded image is scaled to [0, 1].

In function BallisticFallODE, implement the ODE describing the ballistic fall motion according to (1). Also implement proper calls of the ODE implementation by the numerical solver in GenerateBallisticFall function.

Implement the trajectory risk estimation in get_risk function. The risk is given based to (3) as a sum of risk associated of each sample of the given trajectory path.

1. For each sample, properly align the precomputed generic ballistic fall dummy_fall with 'q', and determine the crash location and fallen altitude.
2. Based on the initial heading at q and the fallen altitude, select the most suitable impact probability map.
3. Properly offset the selected impact probability map w.r.t. q.
4. Based on the selected impact probability map and the given ground map, determine the risk associated with the failure occurring at q according to (2).

Two maneuvers are tested in the provided evaluation script t3c-risk.py. The expected results are as follows: Maneuver risk: 0.40167636334498374
Maneuver risk: 0.02036248426234383

Note that the impact risk maps are calculated on a randomized data, and so a slight risk deviation is possible. However, the risk should not deviate more that 3%.

1) Download prepared codes and configuration files.

2) Install dubins package to python3.

On Ubuntu 18.04 or Ubuntu 20.04 (or Python < 3.9), install the package using

pip3 install dubins # or our favorite way to install the package

./install_dubins.sh