Aerial navigation in obstructed environments with embedded nonlinear model predictive control

TL;DR

This paper introduces a real-time nonlinear model predictive control approach for autonomous aerial navigation and obstacle avoidance in complex environments, demonstrated on a micro aerial vehicle with embedded computation.

Contribution

It presents a novel NMPC methodology using PANOC for fast, embedded obstacle avoidance in non-convex environments, with a simple battery depletion compensation method.

Findings

NMPC with PANOC runs at 20Hz onboard MAV.

The MAV successfully navigates around obstacles smoothly.

Thrust compensation improves autonomy over time.

Abstract

We propose a methodology for autonomous aerial navigation and obstacle avoidance of micro aerial vehicles (MAV) using nonlinear model predictive control (NMPC) and we demonstrate its effectiveness with laboratory experiments. The proposed methodology can accommodate obstacles of arbitrary, potentially non-convex, geometry. The NMPC problem is solved using PANOC: a fast numerical optimization method which is completely matrix-free, is not sensitive to ill conditioning, involves only simple algebraic operations and is suitable for embedded NMPC. A C89 implementation of PANOC solves the NMPC problem at a rate of 20Hz on board a lab-scale MAV. The MAV performs smooth maneuvers moving around an obstacle. For increased autonomy, we propose a simple method to compensate for the reduction of thrust over time, which comes from the depletion of the MAV's battery, by estimating the thrust constant.

Figures

Tables

Table 1. TABLE I: MAV parameters — values in SI units

parameter	value	parameter	value
$A_x$	0.1	${\tau }_r$	0.5
$A_y$	0.1	$K_r$	1
$A_z$	0.2	${\tau }_p$	0.5
		$K_p$	1

Equations

$$R({\theta }_r,{\theta }_p)=R_y({\theta }_p)R_x({\theta }_r),$$

$$F=C_T{u}_T^2,$$

$$a=\frac{F}{m}=\frac{C_T}{m}{u}_T^2=C{u}_T^2,$$

$$\dot{C}={\sigma }_C^{2}w,$$

$$\tilde{C}=\frac{a_m}{u_T^2},$$

$$u_T=\sqrt{\frac{T_d}{\hat{C}}},$$

$$\dot{x}(t)=f(x(t),u(t)),$$

$$c_{\iota }(p(t))\notin {\Theta }_j(t)\text{dfn}{O}_j(t)+{\mathcal{B}}_r(r_{\mathrm{ball}}),$$

$$\ell (x,u,t)=$$

$${\ell }_f(x,t)=$$

$$\displaystyle\operatorname*{minimize}\$$

$$\displaystyle\operatorname*{subject\ to}\$$

$$\dot{\bar{x}}(t)=f(\bar{x}(t),\bar{u}(t)),$$

$$\bar{u}(t)\in U(t),$$

$$c_{\iota }(\bar{p}(t))\notin {\Theta }_j(t),j\in {\mathbb{N}}_{[1,q(t)]},$$

$$\iota \in {\mathbb{N}}_{[1,n_{\mathrm{uav}}]},t\in [0,T]$$

$$J={\ell }_f(\bar{x}(T),T)+\sum _{k=0}^{N-1}{\int }_{k{T}_s}^{(k+1)T_s}\ell (\bar{x}(\tau ),{\bar{u}}_k,\tau )\text{dd}\tau .$$

$${\bar{x}}_{k+1}\approx f_k({\bar{x}}_k,{\bar{u}}_k),$$

$$J\approx {\ell }_N({\bar{x}}_N)+\sum _{k=0}^{N-1}{\ell }_k({\bar{x}}_k,{\bar{u}}_k).$$

$${\Theta }_j(t)=\{p\in {\mathbb{R}}^3\mid h_j^i(p,t)>0,i\in {\mathbb{N}}_{[1,m_j(t)]}\},$$

$${\psi }_{\Theta (t)}(c_{\iota }(p(t)))=0,$$

$${\psi }_{\Theta (t)}(p)\text{dfn}\frac{1}{2}\prod _{i=1}^m[h^i(p,t){]}_{+}^2.$$

$$\nabla {\psi }_{\Theta (t)}(p)=1_{\Theta (t)}(p)\sum _{i=1}^m{h}^i(p,t)\prod _{j\ne i}(h^j(p,t){)}^2{\nabla }_p{h}^i(p,t),$$

$$\tilde{\ell }(\bar{x},\bar{u},t)=$$

$${\tilde{\ell }}_f(\bar{x},\bar{u},t)=$$

$$\displaystyle\operatorname*{minimize}_{\bar{u}_{0},\ldots,\bar{u}_{N-1}}{}\$$

$$\displaystyle\operatorname*{subject\ to}\$$

$${\bar{x}}_{k+1}=f_k({\bar{x}}_k,{\bar{u}}_k),k\in {\mathbb{N}}_{[0,N-1]}$$

$${\bar{u}}_k\in U_k,k\in {\mathbb{N}}_{[0,N-1]}$$

$$\underset{\bar{u}\in U}{\mathrm{minimize}}\varphi (\bar{u};\hat{x},p^{\mathrm{ref}}(\cdot )),$$

$$F_0(\bar{u})=$$

$$F_{k+1}(\bar{u})=$$

$$\varphi (\bar{u})={\tilde{\ell }}_N(F_N(\bar{u}))+\sum _{k=0}^{N-1}{\tilde{\ell }}_k(F_k(\bar{u}),{\bar{u}}_k).$$

$$R_{\gamma }({\bar{u}}^{\star })=0,$$

$$R_{\gamma }(\bar{u})=\bar{u}-T_{\gamma }(\bar{u}),$$

$$T_{\gamma }(\bar{u})={\Pi }_U(\bar{u}-\gamma \nabla \varphi (\bar{u})).$$

Full text

I Introduction

I-A Background and motivation

The need for safe aerial navigation and increased MAV autonomy nowadays poses all the more relevant and pressing research questions, as drones make their appearance in numerous application domains, such as the inspection of critical or aging infrastructure [1], surveying of underground mines [2], visual area coverage for search-and-rescue operations [3], precision agriculture [4] and many others. In the majority of these applications, MAVs have to navigate in obstructed environments, with static or moving obstacles of arbitrary geometry in known, or partially unknown surrounding environments.

Several methods have been proposed for navigation and collision avoidance, such as potential field methods [5, 6] and graph search methods [7]. Alongside these methods, NMPC is becoming popular for the navigation control of various MAVs including fixed-wing aircrafts [8, 9] and multi-rotor vehicles [10]. NMPC uses a nonlinear dynamical model of the system dynamics to predict position and attitude trajectories from its current position to a reference point, while avoiding all obstacles on its way and minimizing a certain energy/cost function. In this way, a non-convex optimization problem needs to be solved at every sampling time instant in a receding horizon fashion. Another approach to obstacle avoidance is described in [11] where a high-level path planner generates collision-free trajectories which are followed by an MPC controller.

In [12], sequential quadratic programming (SQP) is used to solve the NMPC problem for the navigation of a multi-rotor MAV with a slung load, where the authors demonstrated the effectiveness of NMPC, however, provided neither evidence of the solution quality or solver performance, nor an experimental verification. NMPC was used in [13] for solving obstacle and collision avoidance for several MAVs flying in formation, however again, only simulations were done and the computation time was addressed.

Clearly, the presence of obstacle/collision avoidance constraints makes the MPC problems particularly hard to solve. SQP is the method of choice in the literature [12, 13, 14] that has as a main disadvantage the fact that it requires the solution of a quadratic program (QP) at every iteration of the algorithm, which requires inner iterations. SQP also requires computing and storing of the Jacobian matrices of the dynamics, and sometimes the Hessians when the Hessian of the Lagrangian is used in the QPs. Furthermore, the gradient descent method has been used to solve nonlinear MPC problems for aerial navigation [14]. This method, however, is sensitive to bad conditioning and problems with long horizons tend to become ill conditioned, while the convergence is expected to be slow.

I-B Contributions

In this article we propose a control methodology for the autonomous navigation of MAVs in obstructed environments. We allow for the obstacles to have arbitrary non-convex shapes and, contrary to distance-based methods [15], we do not require that the distance function between the MAVs and each obstacle is available.

The NMPC optimization problem is solved by using PANOC [16, 17] — a recently proposed algorithm for non-convex optimization problems, which is suitable for embedded NMPC, as it requires only simple and cheap linear operations (mainly inner products of vectors) and exhibits a fast convergence. Unlike SQP, PANOC is matrix-free and only requires the computation of Jacobian-vector products, which can be computed very efficiently by backward (adjoint) automatic differentiation. PANOC has been shown in [16, 17, 18] to significantly outperform both SQP and interior-point methods. To the authors best of knowledge, this is the first time that a fast NMPC optimization problem is being demonstrated on an aerial platform, setting the base for future developments in the aerial robotics community.

Our method for modeling has the strong merit of being independent of the mass of the MAV, whereas the norm in the community is to use mass and other detailed parameters of the specific MAV used, for example [10], [11], and [13]. This allows for our method to be used without tuning the specific mass or available thrust, improving robustness, generalization and ease of use.

Evidence of the solution quality is provided by physical laboratory experiments, where a MAV is flown completely autonomously in a laboratory equipped with a VICON motion capture system. The proposed method uses a full position and attitude model of the MAVs, which is able to run onboard, using 8-15% CPU of a single core on an Intel Atom Z8350. As is shown in Section V the onboard controller is able to successfully navigate the MAV around an obstacle running at a sampling rate of $20\,\mathrm{Hz}$ and a prediction horizon of $2\,\mathrm{s}$.

II MAV dynamics

II-A MAV* kinematics*

The model of a quadrotor MAV, defined by [9], assumes that there exists a low-level controller of roll, pitch, yaw rate and thrust. This convention is common in MAV flight controllers such as PixHawk, [19] and ROSFlight, [20]. The high-level kinematics of the MAV is given by

[TABLE]

where $p(t)=(p_{x}(t),p_{y}(t),p_{z}(t))\in\R^{3}$ and $v(t)\in\R^{3}$ are the position and velocity of the MAV in the global frame of reference, and $\theta_{r}\in\R$ and $\theta_{p}\in\R$ are the roll and pitch angles, while $\theta_{r,d}\in\R$ and $\theta_{p,d}\in\R$ are the reference angles sent to the low-level controller. Furthermore, $T_{d}\in\Rplus$ is the $z$-axis thrust acceleration, while $A_{x}$, $A_{y}$, and $A_{z}$ are the linear drag coefficients. The lower layer — the attitude control system — is modeled by simple first-order dynamics with time constant $\tau_{r}$ and $\tau_{p}$ and gains $K_{r}$ and $K_{p}$ for the roll and pitch respectively. Lastly, $R(\theta_{r},\theta_{p})\in\mathrm{SO}(3)$ describes the MAV’s attitude and is defined by the classical Euler angles in rotation matrix form as

[TABLE]

with

[TABLE]

Note that yaw is absent in this rotation matrix, as this model operates in a yaw-compensated global frame, and the position control of the MAV is therefore independent of its yaw. Moreover, it is important to note that we have chosen the acceleration, $T_{d}$, to be the manipulated variable of the system, rather than the corresponding force, for the model to be mass-free. This has the strong merit of making the controller robust to changes in the mass of the MAV, the available thrust from the motors, and the loss of thrust over time due to the decline of battery voltage.

II-B Adaptive acceleration control

In order for our design to be independent of the physical characteristics that determine the available thrust acceleration, a simplified version of [21] is used to continuously estimate the MAV’s maximum available thrust. Following [21], the force, $F$, that is exercised by the propellers, is given by

[TABLE]

where $C_{T}$ is the thrust constant and $u_{T}\in[0,~{}1]$ is a unitless normalized thrust control signal. The thrust constant is time dependent based, for instance, on battery drain and how close the MAV is to the ground, which is why identifying a constant is not sufficient to track thrust references. The issue is that there is no sensor in the system, which measures the generated force, however, the IMU can provide a measurement of the linear acceleration, albeit noisy. Then, by dividing the thrust model by the mass of the MAV, $m$, the model is now based on acceleration:

[TABLE]

Equations (1ca) and (1cb) define a nonlinear dynamical system with state variable $C$, input $u_{T}$ and output $y(C,u_{T})=Cu_{T}^{2}$. We estimate $C$ by means of an extended Kalman filter (EKF). EKF is chosen because it is simple to tune, it allows to specify an initial estimate variance, and converges fast in the first few iterations.

In additional, we employ an outlier rejection scheme based on bounds of the direct estimate $\widetilde{C}\in[1g,~{}10g]$, where

[TABLE]

which is calculated for each IMU acceleration measurement $a_{m}$, which implies that each acceleration measurement is inspected to enforce that no outliers are allowed to update the filter. These bounds result from the fact that a MAV must be able to generate at least $1g$ of thrust to take off and it is assumed that it cannot generate more than $10g$ of thrust. The bounds on $\widetilde{C}$ are inherited by the estimates $\widehat{C}$ yielding a simple constrained estimation scheme.

Once the thrust constant is estimated, an acceleration reference can be converted to the thrust control signal $u_{T}$, by solving equation (1ca) for $u_{T}$, resulting in

[TABLE]

A depiction of how the thrust constant estimation is tied to the overall scheme can be found in Fig. 1.

II-C Overall system dynamics

The state of the controlled system is defined to be $x(t)=(p(t),v(t),\theta_{r}(t),\theta_{p}(t))$ and the manipulated input is $u(t)=(T_{d}(t),\theta_{r,d}(t),\theta_{p,d}(t))$. The system is observed using a VICON motion capture system, which measures the full odometry of the system and provides the corresponding estimates of the full state of the MAV as $\hat{x}=(\hat{p}(t),\hat{v}(t),\hat{\theta}_{r}(t),\hat{\theta}_{p}(t))$. Overall, the system dynamics can be concisely written as

[TABLE]

where $f$ is implicitly defined via (1).

III Nonlinear MPC for obstacle avoidance

III-A Navigation in obstructed environments

We assume that a MAV needs to navigate towards a reference position $p^{\mathrm{ref}}\in\R^{3}$, while avoiding a set of $q(t)\in\N$ moving obstacles $\{O_{j}(t)\}_{j\in\N_{[1,q(t)]}}$.

We select $n_{\mathrm{mav}}$ corner points on the MAV and position a ball with radius $r_{\mathrm{ball}}$ centered at each such point so that the whole vehicle is contained in the union of these balls. We assume that the coordinates of the corner points in the global frame of reference are given by $c_{\iota}(p(t))$, for $\iota\in\N_{[1,n_{\mathrm{mav}}]}$.

In order for the MAV to not collide with the obstacles, we shall require that

[TABLE]

for all $j\in\N_{[1,q(t)]}$, $\iota\in\N_{[1,n_{\mathrm{mav}}]}$, where $\mathcal{B}_{r}(r_{\mathrm{ball}})$ is a ball centered at the origin with radius $r_{\mathrm{ball}}$. The set $\Theta_{j}(t)$ is an enlarged version of the original obstacle $O_{j}(t)$. The concept is illustrated in Fig. 2.

We introduce the stage cost function $\ell{}:{}\R^{n_{x}}{}\times{}\R^{n_{u}}{}\times{}\R_{+}{}\to{}\R_{+}$ and the terminal cost function $\ell_{f}{}:{}\R^{n_{x}}{}\times{}\R_{+}{}\to{}\R_{+}$ which penalize the deviation of the system state from a reference state. Typical choices are

[TABLE]

where $Q\in\R^{n_{x}\times n_{x}}$, $R\in\R^{n_{u}\times n_{u}}$ and $Q_{f}\in\R^{n_{x}\times n_{x}}$ are positive semi-definite matrices and $x^{\mathrm{ref}}$ is the reference state which has the form $x^{\mathrm{ref}}=[p^{\mathrm{ref}}~{}0_{1\times n_{x}-3}]^{\top}$.

The nonlinear model predictive control problem for navigation in an obstructed environment consists in solving the following problem

[TABLE]

where $\bar{u}(t)=(\bar{T}_{d}(t),\bar{\theta}_{r,d}(t),\bar{\theta}_{p,d}(t))$ and $\bar{x}(t)=(\bar{p}(t),\bar{v}(t),$ $\bar{\theta}_{r}(t),\bar{\theta}_{p}(t))$, for $t\in[0,T]$ are the predicted input and state signals.

In this formulation we have assumed that the future trajectories of all obstacles are exactly known and independent of the trajectory of the controlled vehicle. If this is not the case, we have to formulate appropriate robust or stochastic variants of the above obstacle avoidance problem.

The control action is exercised to the system via a zero-order hold element, that is, $\bar{u}(t)=\bar{u}_{k}$ for $t\in[kT_{s},(k+1)T_{s})$, where $T_{s}$ is the sampling period. We assume that $T=NT_{s}$ for some $N\in\N$. Then, the cost function in (1ia) can be written as

[TABLE]

Since it is not possible to derive analytical solutions of the nonlinear dynamical system (1ic), the system trajectories as well as the cost function $J$ along these trajectories has to be evaluated by discretizing the system dynamics and integrals. By doing so, the system state trajectoriy $\bar{x}(t)$ is evaluated at points $\bar{x}_{k}=\bar{x}(kT_{s})$ as follows

[TABLE]

and

[TABLE]

Any explicit integration method such as the fourth-order Runge-Kutta or Forward Euler lead to high quality approximations of MAV trajectories. This way, the original continuous-time optimal control problem is approximated by a discrete-time one which is solved at every time instant in a receding horizon fashion.

III-B Penalty functions for obstacles of general shape

Each obstacle is described by a set of $m_{j}(t)$ nonlinear constraints of the form

[TABLE]

where functions $h_{j}^{i}{}:{}\R^{3}\times\R_{+}\to\R_{+}$ are $C^{1,1}$ functions. This approach allows one to describe obstacles of very general convex or nonconvex shape. For example, by choosing functions $h_{j}^{i}$ to be affine in $p$, we can model any polytopic object. Functions of the form $h_{j}(p,t)=1-(p-p_{0}(t))^{\top}M(t)(p-p_{0}(t))$ can be used to model ellipsoidal objects or elliptic cylindrical ones. Polynomial, trigonometric and other functions can be used to model more complex geometric shapes.

For simplicity, in this section we focus in the case where there is one obstacle, that is $q(t)=1$, which we denote by $\Theta(t){}={}\{p\in\R^{3}{}\mid{}h^{i}(p,t)>0,i\in\N_{[1,m]}\}$. The constraint $c_{\iota}(p(t))\notin O(t)$ is satisfied if and only if $h^{i_{0}}(c_{\iota}(p(t)),t)\leq 0$ for some $i_{0}\in\N_{[1,m]}$, or equivalently, if

[TABLE]

for all $\iota\in\N_{[1,n_{\mathrm{uav}}]}$, where $\psi_{\Theta(t)}:\R^{3}\to\R_{+}$ is the function defined as

[TABLE]

Such a function is illustrated in Fig. 3. Function $\psi_{\Theta(t)}(p)$ takes the value [math] outside the enlarged obstacle $\Theta(t)$ and increases in the interior of it as we move away from its boundary.

Function $\psi_{\Theta(t)}$ is differentiable with gradient

[TABLE]

where $1_{\Theta(t)}$ is the characteristic function of $\Theta(t)$ with $1_{\Theta(t)}(p){}={}1$ if $p\in\Theta(t)$ and $1_{\Theta(t)}(p){}={}0$ otherwise.

Functions $\psi_{\Theta(t)}$ can be used to impose the obstacle avoidance requirements as soft constraints. To this end, we eliminate the non-convex constraints $c_{\iota}(\bar{p}(t))\notin\Theta_{j}(t)$ and introduce the modified stage and terminal cost functions

[TABLE]

where $\lambda_{j,\iota}$ and $\lambda^{f}_{j,\iota}$ are positive weight coefficients. The overall model predictive control (MPC) problem becomes

[TABLE]

The optimization is carried out over finite-dimensional vectors $\bar{u}=(\bar{u}_{0},\ldots,\bar{u}_{N-1})\in\R^{n}$ with $n{}={}n_{u}(N-1)$.

III-C Single-shooting problem formulation

We shall cast optimization problem (1o) in the following compact and simple form

[TABLE]

where $U{}={}U_{0}{}\times{}U_{1}{}\times{}\ldots{}\times{}U_{N-1}$ and $\phi{}:{}\R^{n}{}\to{}\R_{+}$ is a $C^{1,1}$ function. To this end, we need to eliminate the state sequence in (1oc). Let us introduce a sequence of functions $F_{k}{}:{}\R^{n}{}\to{}\R^{n_{x}}$ for $k\in\N_{[0,N]}$ defined recursively by

[TABLE]

Then, problem (1o) is written as in (1p) with

[TABLE]

This is known as the single shooting formulation [16].

IV Fast online nonlinear MPC using PANOC

Problem (1p) is in a form that can be solved by PANOC [16]. In particular, the gradient of $\phi$ can be computed using automatic differentiation [22] which is implemented by software such as CasADi [23]. PANOC finds a $\bar{u}^{\star}\in\R^{n}$ which solves the optimality conditions

[TABLE]

where $R_{\gamma}:\R^{n}{}\to{}\R^{n}$ is the fixed-point residual operator with parameter $\gamma{}>{}0$ defined as

[TABLE]

where $T_{\gamma}:\R^{n}{}\to{}\R^{n}$ is the projected gradient operator given by

[TABLE]

PANOC combines safe projected-gradient updates $\bar{u}^{\nu+\nicefrac{{1}}{{2}}}$ with fast Newton-type directions $d^{\nu}$ computed by L-BFGS while it uses the forward backward envelope (FBE) function $\varphi_{\gamma}$ as a merit function for globalization given by

[TABLE]

The forward-backward envelope is an exact, continuous and real-valued merit function which shares the same (local/strong) minima with (1p). That said, Problem (1p) is reduced to the unconstrained minimization of $\varphi_{\gamma}$.

PANOC is shown in Algorithm 1. L-BFGS uses a buffer of length $\mu$ of vectors $s^{\nu}{}={}\bar{u}^{\nu+1}{}-{}\bar{u}^{\nu}$ and $y^{\nu}{}={}R_{\gamma}(\bar{u}^{\nu+1}){}-{}R_{\gamma}(\bar{u}^{\nu})$ to compute the update directions $d^{\nu}$ [16], [24, Sec. 7.2]. The computation of $d^{\nu}$ requires only inner products which amount to a maximum of $4\mu n$ scalar multiplications. In particular, following [25], the L-BFGS buffer is updated only if $s^{\nu\top}y^{\nu}/\|s^{\nu}\|^{2}{}>{}\epsilon_{d}\|R_{\gamma}(\bar{u})\|$.

Overall, PANOC uses exactly the same oracle as the projected gradient method, that is it only requires the invocation of $\operatorname{\Pi}_{U}$, $\phi$ and $\nabla\phi$. Lastly, owing to the FBE-based line search, PANOC converges globally, that is, from any initial guess, $\bar{u}^{0}$.

V Experimental validation

For the experimental validation of the proposed control scheme, an inverted quadrotor using the ROSFlight [20] low-level controller was used for all trials, shown in Fig. 4. The onboard computer used is an Aaeon Up Board with an Intel Atom x5-z8350 processor with four $1.44\,\mathrm{GHz}$ cores and $2\,\mathrm{GB}$ of RAM. The board runs Ubuntu Server 16.04. The field robotics lab at LuleåUniversity of technology is equipped with a Vicon motion capture system featuring 20 infrared cameras that track the odometry of the MAV; this data is used by the NMPC controller for navigation.

The NMPC module runs simple C89 code which was generated by nmpc-codegen — an LGPLv3.0-licensed open-source code generation toolkit which is available at github.com/kul-forbes/nmpc-codegen.

An upright cylindrical obstacle, $O$, is placed so that its vertical symmetry axis runs through the origin $(0.0,0.0,0.0)$ of the global coordinate frame in the flying arena at field robotics lab (FROST). The cylinder, $O$, has a radius of $r_{\mathrm{cyl}}=0.45\,\mathrm{m}$ and height $z_{\mathrm{cyl}}{}={}2\,\mathrm{m}$. The obstacle is described by the functions $h^{1}(p,t){}={}r^{2}-p_{x}^{2}-p_{y}^{2}$, $h^{2}(p,t){}={}p_{z}$ and $h^{3}(p,t){}={}z_{\mathrm{cyl}}-p_{z}$. A single corner point is used which is positioned at the center of the MAV; the enclosing ball as in Fig. 2 has a radius of $r_{\mathrm{ball}}{}={}0.24\,\mathrm{m}$. In order to account for possible small constraint violations due to the fact that obstacle avoidance constraints are modeled via penalty functions, we consider an additional enlargement of $0.06\,\mathrm{m}$. As a result, the enlarged cylinder, $\Theta(t)$, has a radius of $0.75\,\mathrm{m}$ and height $2.3\,\mathrm{m}$. The weights of the obstacle constraints, $\lambda_{j,\iota}$ and $\lambda^{f}_{j,\iota}$, in Equation (1n) were all set to $10000$, and the continuous-time system was integrate with the forward Euler method.

The flight test performed for avoiding the obstacle consisted of alternating between two position references on opposite sides of the obstacle. The two position references given alternately were, in meters, $(-2.0,0.0,1.0)$ and $(2.0,0.0,1.5)$. These references were sent when the MAV was close to its previous reference position. The exact time the references are changed can be seen in Fig. 6.

NMPC runs at $20\,\mathrm{Hz}$ with a prediction and control horizon of 40 steps, meaning the solver predicts the states of the system $2\,\mathrm{s}$ into the future. The solver occupied between $8\%$ and $15\%$ of CPU on an Intel Atom Z8350 — an indication of the solver’s computational efficiency.

Fig. 5 shows the actual path flown by the MAV during the test where the positioning data is taken from the motion tracking system and has sub-millimeter accuracy. The path is also shown in Fig. 6 where we plot the MAV’s position versus time. The MAV does not have time to settle at the reference altitude as a new reference is sent to the controller before the position completely converges.

As the MAV passes the obstacle it violates the obstacle constraint, as shown in Fig. 7, which is expected from the penalty formulation. The maximum violation is $2.86\,\mathrm{cm}$, which is below the extra enlargement of $6\,\mathrm{cm}$ of the obstacle.

Fig. 8 shows the control signals (roll, pitch, and normalized thrust references) commanded by the NMPC. The roll and pitch angles are bound between $-0.5\,\mathrm{rad}$ and $0.5\,\mathrm{rad}$; these bounds are active as shown in Fig. 8. This further motivates the use of NMPC, allowing for bounds to be directly included in the problem formulation.

The control signals could be made less aggressive by penalizing the rate of change of the input in (1o), that is, by adding a penalty of the form $\ell_{\Delta}{}={}\|\bar{u}_{k}-\bar{u}_{k-1}\|^{2}_{R_{\Delta}}$ for a symmetric positive semidefinite matrix $R_{\Delta}\in\R^{n_{u}\times n_{u}}$. Nevertheless, the maneuvering of the MAV is smooth as shown in Figs. 5, 6 and 7 and a video of the experiment which can be found at https://youtu.be/E4vCSJw97FQ.

As shown in the second subfigure of Fig. 9, once the reference changes, the solver reaches the maximum number of iterations (200 iterations) and the solution it returns is of poor quality (fourth subfigure of Fig. 9). This happens because at each time instance, the solver is initialized with the previously computed optimal trajectory. Upon a reference change, the initial guess is rather far from optimal and this necessitates more iterations for convergence. Nonetheless, this inaccuracy is eliminated at the next time instant — $0.05\,\mathrm{s}$ later — where the solver is provided a good initial estimate and converges within the prescribed tolerance, $\epsilon=10^{-3}$. This way, NMPC is executed at $20\,\mathrm{Hz}$. As shown in the third subfigure of Fig. 9, at one time instant, the solution time exceeds the maximum allowed time. This is accommodated by delaying the dispatch of the control action by few $\mathrm{ms}$ and has no practical effect.

The infinity norm of the fixed-point residual is below $\epsilon$ at all time instants with the exception of four instants from the change of reference. Lastly, the average iteration time in every time step is shown in the third subfigure of Fig. 9, and ranges from $80\,\mathrm{\mu{}s}$ to $350\,\mathrm{\mu{}s}$ where the variability is because of the different number of line search iterations.

The parameters used in the dynamics of the MAV used in the experiment are shown in Table I. These values were chosen empirically (based on accurate values for other MAVs) and are not fine-tuned via experiments; this accentuates the fact that the closed-loop and the overall obstacle avoidance scheme is robust to errors in the determination of these parameters.

The tuning parameters used by the NMPC are

[TABLE]

and the prediction horizon $T=2\,\mathrm{s}$. For the EKF for estimating the special thrust constant we have

[TABLE]

where $P_{0}$ is the initial variance for $\widehat{C}$, $Q_{T}$ is the process variance in (1cb), and $R_{T}$ is the measurement variance.

A separate experiment was carried out where the MAV was given a position reference to hold for as long as the battery could deliver power safely. This experiment was conducted to demonstrate the thrust constant estimation described in Section II-B and the results are presented in Fig. 10. As the battery drains, the special thrust constant is decreasing and the control signal is adapted to keep the MAV hovering at a constant altitude. This experiment is part of the same video mentioned in this section, found at https://youtu.be/E4vCSJw97FQ.

VI Conclusions

We presented an obstacle and collision avoidance methodology coupled with an adaptive thrust controller that leads to increased autonomy and context awareness for MAVs. Obstacle avoidance is addressed with an NMPC controller, which is solved using PANOC — a simple and fast algorithm, which involves simple algebraic operations and, unlike SQP, does not require the solution of linear systems at each step. Experiments were performed with the solver running onboard a MAV which maneuvered gently around a virtual obstacle with a smooth trajectory. The MAV passed the edge of the virtual obstacle with a minimal constraint violation, as expected from the solver.

Moreover, experiments were performed to demonstrate that our thrust estimation method successfully compensates for the reduction of thrust over time, making the control scheme applicable to any MAV platform.

Future work will focus on experiments in presence of moving obstacles with uncertain trajectories.