We consider the problem of risk-sensitive motion planning in the presence of randomly moving obstacles. To this end, we adopt a model predictive control (MPC) scheme and pose the obstacle avoidance constraint in the MPC problem as a distributionally robust constraint with a KL divergence ambiguity set. This constraint is the dual representation of the Entropic Value-at-Risk (EVaR). Building upon this viewpoint, we propose an algorithm to follow waypoints and discuss its feasibility and completion in finite time. We compare the policies obtained using EVaR with those obtained using another common coherent risk measure, Conditional Value-at-Risk (CVaR), via numerical experiments for a 2D system. We also implement the waypoint following algorithm on a 3D quadcopter simulation.

§ INTRODUCTION

Emerging applications in robot path planning in unknown and partially known unstructured environments, such as search and rescue missions caused by natural disasters [25, 31], inspection of planetary terrains [17], and exploration of urban underground environments [21],  motivate the need for risk-sensitive path planning. In particular,  path planning in subterranean environments [30] incurs higher risks due to  lack of Global Positioning System (GPS) signals, the absence of illumination,  decentralization, and unpredictable environment topologies [27, 24] (see Fig. 1).

Motion planning risk can be quantified in multiple ways, such as chance constraints [26, 35], exponential utility functions [20], and distributional robustness [36]. However, applications in autonomy and robotics require more “nuanced assessments of risk” [23]. Artzner et. al. [7] characterized a set of natural properties that are desirable for a risk measure.  These coherent risk measures are widely used and accepted in finance and operations research, among other fields.

CVaR is an important example of a coherent risk measure that has received significant attention in decision making problems, such as Markov decision processes (MDPs) [13, 12, 28, 8]. For stochastic discrete-time dynamical systems, a MPC technique for a class of coherent risk objectives that admit polytopic representation was proposed in [32]. These authors also proposed a Lyapunov condition for risk-sensitive exponential stability in the presence of discretely quantized process noise, but did not include constraints in their formulation. Measurement noise and moving obstacles were considered in [15], wherein the authors devised an MPC-based scheme for path planning with CVaR safety constraints when a reference
trajectory is generated by RRT$^*$ [18], and extended to a Wasserstein distributionally robust formulation in [16]. Risk-sensitive obstacle avoidance has also been tackled through CVaR control barrier functions in [2] with application to bipedal robot locomotion. Moreover, a method based on stochastic reachability analysis was proposed in [11] to estimate a CVaR-safe set of initial conditions via the solution to an MDP.

Spot and Husky robots
exploring a subterranean environment in Valentine Cave, Lava Beds National Monument, California. Avoidance of obstacles and other moving robots in unstructured environments incurs higher mission risk due to lack of global positioning [9].

Despite the popularity of CVaR in risk-sensitive path planning, CVaR is hard to compute efficiently, even for the sum of arbitrary independent random variables [4, 3]. In most cases, one has to approximate CVaR through sampling methods. Furthermore, CVaR only considers the average worst case performance beyond a threshold, while ignoring the performance before reaching that threshold [6].

EVaR is the tightest upper bound on CVaR and Value-at-Risk (VaR) in the sense of the Chernoff inequality (and hence a more risk-sensitive measure) [4, 6, 5]. Nonetheless, it does not admit a polytopic representation and therefore methods such as the ones proposed in [32] cannot be applied for risk-sensitive path planning. In [33], risk-constrained and risk-averse optimal control, amenable to arbitrary coherent risk measures, is considered. The authors reformulate the optimal control optimization as a convex conic program. This formulation, however, does not consider the nonconvex, mixed-integer nature of the optimization problems that is often a result of obstacle avoidance constraints.

In this paper, we go beyond CVaR path planning and propose a framework for receding horizon path planning with risk-sensitive obstacle avoidance and guaranteed performance in terms of EVaR. We consider discrete-time systems and a class of randomly moving obstacles for which we reformulate the MPC optimization as a convex, mixed-integer program. This is done in three steps - first we write the EVaR constraint as a cone constraint, next we reformulate the obstacle avoidance constraints to obtain a mixed-integer relaxation, and lastly we add a discrete state that tracks whether the goal has been reached. This allows us to track waypoints in a way that guarantees feasibility and finite-time task completion. We elucidate the proposed method using two examples.

This paper is organized as follows. The next section briefly reviews some relevant facts on CVaR and EVaR. Section III presents the problem under study in this paper. In Section IV, we propose a reformulation based on convex mixed integer programming to solve the EVaR receding horizon path planning problem and discuss its feasibility properties. In Section V, we introduce an algorithm to follow waypoints using the aforementioned MPC optimization and prove its finite-time completion. Section VI illustrates the method via numerical experiments. Section VII concludes the paper.

Notation:  We denote by $\mathbb{R}^n$ the $n$-dimensional Euclidean space, $\mathbb{R}_{\ge 0}$ the non-negative reals, and $\mathbb{N}_{\ge0}$ the set of non-negative integers. Throughout the paper, we use bold font to denote a vector and $(\cdot)^\top$ for its transpose, e.g., $\boldsymbol{a}=(a_1,\ldots,a_n)^\top$, with $n\in \{1,2,\ldots\}$. For vector $\boldsymbol{a}$, we use $\boldsymbol{a}\succeq (\preceq) \boldsymbol{0}$ to denote element-wise non-negativity (non-positivity) and $\boldsymbol{a}\equiv \boldsymbol{0}$ to show all elements of $\boldsymbol{a}$ are zero. For two vectors $a,b \in \mathbb{R}^n$, we denote their inner product by $\langle \boldsymbol{a}, \boldsymbol{b} \rangle$, i.e., $\langle \boldsymbol{a}, \boldsymbol{b} \rangle=\boldsymbol{a}^\top \boldsymbol{b}$. In the MPC problem, we refer to $\boldsymbol{x}(t+k|t)$ as $\boldsymbol{x}_k$. For a finite set $\mathcal{A}$, we denote its power set by $2^\mathcal{A}$.

For  a probability space $(\Omega, \mathcal{F}, \mathbb{P})$ and a constant $p \in [1,\infty)$, $\mathcal{L}_p(\Omega, \mathcal{F}, \mathbb{P})$ denotes the vector space of real valued random variables $X$ for which $\mathbb{E}|X|^p < \infty$. For two probability density functions $P(X)$ and $Q(X)$, $P \ll Q$ implies that $P$ is absolutely continuous with respect to $Q$, i.e., if $Q(X)=0$, then $P(X)=0$.

§ PRELIMINARIES

This section reviews some results on CVaR and EVaR risk measures.

Conditional Value-at-Risk

For a given confidence level $\alpha \in (0,1)$, value-at-risk ($\mathrm{VaR}_{1-\alpha}$) denotes the $({1-\alpha})$-quantile value of the cost variable $X \in \mathcal{L}_p(\Omega, \mathcal{F}, \mathbb{P})$. $\mathrm{CVaR}_{1-\alpha}$ measures the expected loss in the $({1-\alpha})$-tail given that the particular threshold $\mathrm{VaR}_{1-\alpha}$ has been crossed. $\mathrm{CVaR}_{1-\alpha}$ is given by

\begin{equation}
\begin{aligned}
\mathrm{CVaR}_{1-\alpha}(X):=&\inf_{z \in \mathbb{R}}\mathbb{E}\Bigg[z + \frac{(X-z)^{+}}{1-\alpha}\Bigg], %\\
% =& \sup_{Q\in\mathcal{Q}}\mathbb{E}_Q(X)
\end{aligned}
\end{equation}

where $(\cdot)_{+}=\max\{\cdot, 0\}$. A value of $\alpha \simeq 0$ corresponds to a risk-neutral case; whereas, a value of $\alpha \to 1$ is rather a risk-averse case.

Entropic Value-at-Risk

EVaR, derived using the Chernoff inequality for VaR, is the tightest upper bound for VaR and CVaR.  It was shown in [6] that $\mathrm{EVaR}_{1-\alpha}$ and $\mathrm{CVaR}_{1-\alpha}$ are equal only if there are no losses ($X\to -\infty$) below the $\mathrm{VaR}_{1-\alpha}$ threshold. The $\mathrm{EVaR}_{1-\alpha}$ of random variable $X$ is given by

\begin{equation} \label{eq:evar}
\mathrm{EVaR}_{1-\alpha}(X):= \inf_{z > 0 }\Bigg[z^{-1}\ln\frac{\mathbb{E}[e^{Xz}]}{1-\alpha}\Bigg]
= \sup_{Q\in\mathfrak{D}}\mathbb{E}_Q(X).
\end{equation}

Similar to $\mathrm{CVaR}_{1-\alpha}$, for $\mathrm{EVaR}_{1-\alpha}$, the limit $\alpha \to 0$ corresponds to a risk-neutral case; whereas, $\alpha \to 1$ corresponds to a risk-averse case. In fact, it was demonstrated in <cit.> that $\lim_{{\alpha}\to 1} \mathrm{EVaR}_{{1-\alpha}}(X) = \esssup(X)$.

A property of coherent risk measures is that they can be written as the worst-case expectation over a convex, bounded, and closed set of probability mass (or density) functions (pdf/pmf). This is the dual representation of a risk measure and the set is referred to as the risk envelope. For EVaR, the risk envelope $\mathfrak{D}$ for a continuous random variable with the pdf $P$ is defined as the epigraph of the KL divergence, given by,

\begin{multline}
\mathfrak{D}:= \Big\{Q\ll P\, | \\ \, D_{KL}(Q||P):=\int\frac{dQ}{dP}\bigg(\ln\frac{dQ}{dP}\bigg)dP \leq -\ln(1-\alpha)\Big\}.
\end{multline}

$D_{KL}(Q||P)$ is the KL divergence between the two distributions, $Q$ and $P$. For some $x,y \in \mathbb{R}$, $D_{KL}(x||y)$ can be written in the form of the exponential cone, $K_{exp}$:

\begin{equation*}
t\geq x\ln(x/y) \iff  (y, x, -t) \in K_{exp}.
\end{equation*}

Similarly, for a discrete random variable $X \in \{x_1, x_2,\dotsc, x_J\}$ with the pmf given by $p = [p(1), p(2), \dotsc, p(J)]^T$, where $p(j) = \mathbb{P}(X = x_j), j \in \{1,\dotsc, J\}$, the KL divergence is given as

\begin{equation*}
D_{KL}(q||p) := \sum_{j=1}^{J}q(j)\ln\bigg(\frac{q(j)}{p(j)}\bigg), \quad q, p \in \Delta_J.
\end{equation*}

$\Delta_J$ is the probability simplex, $ \Delta_J:= \{q \in \mathbb{R}^J\,|\, q \geq 0, \,\sum_{j=1}^{J}q(j) = 1\}$. Hence, the epigraph of the KL divergence, and equivalently the EVaR risk envelope for a discrete distribution, is exponential cone representable [19] as

\begin{equation} \label{eq:KLcone}
\begin{aligned}
\mathfrak{D}:= \Big\{& q \in \Delta_J \,|\,\exists \delta \in \mathbb{R}^J: \sum_{j=1}^{J}\delta(j) \leq -\ln\alpha, \\
& \big(p(j), q(j), -\delta(j)\big) \in K_{exp},~ \, \forall j\in \{1, \dotsc, J\}\Big\}.
\end{aligned}
\end{equation}

§ PROBLEM STATEMENT

We consider a class of discrete-time systems given by

\begin{equation}\label{eq:sys}
\begin{aligned}
\boldsymbol{x}(t+1) &= A\boldsymbol{x}(t) + B\boldsymbol{u}(t), \\
\boldsymbol{y}(t) &= C\boldsymbol{x}(t) + D\boldsymbol{u}(t),
\end{aligned}
\end{equation}

where $\boldsymbol{x}(t) \in \mathbb{R}^{n_x}$, $\boldsymbol{y}(t) \in \mathbb{R}^{n_y}$, and $\boldsymbol{u}(t) \in \mathbb{R}^{n_u}$ are the system state, output, and controls at time $t$, respectively.
We consider obstacles with index $l \in\mathcal{L}$ that can be approximated by a convex polytope defined by $m_l$ half-spaces in $\mathbb{R}^{n_x}$

\begin{equation} \label{eq: obs_def}
\mathcal{O}_l = \{\boldsymbol{y} \in \mathbb{R}^{n_x} \,|\,\boldsymbol{c}_{i,l}^{T}\boldsymbol{y} \leq \boldsymbol{d}_{i,l}, \, i=1, \dots, m_l\}.
\end{equation}

We allow each polytopic obstacle $O_l$, $l \in \mathcal{L}$, centered at $\boldsymbol{a}_l$ at time $t$ to move randomly. That is, the point set defining obstacle $O_l$, $l \in \mathcal{L}$, at $t+k$ can be written as a random rotation $R_l(t+k)$ and random translation $w_l(t+k)$ of the $l$th obstacle  $\mathcal{O}_l$ from time $t$ to $t+k$ as described below

\begin{equation}\label{eq:obs}
\begin{aligned}
\mathcal{O}_l(t+k) &= R_l(t+k)\mathcal{O}_l(t) + \boldsymbol{w}_l(t+k)\\
&= \bigg\{\boldsymbol{y}(t+k) = R_l(t+k)(\boldsymbol{y}(t) - \boldsymbol{a}_l) + \boldsymbol{a}_l + \\ &\qquad \boldsymbol{w}_l(t+k) \, | \, \boldsymbol{c}_{i,l}^{T}\boldsymbol{y}(t) \leq \boldsymbol{d}_{i,l}, \, i=1, \dots, m_l\bigg\} \\
&= \bigg\{\boldsymbol{y}(t+k) \, | \,\boldsymbol{c}_{i,l}^{T}\Big(R_l^{-1}(t+k)\big(\boldsymbol{y}(t+k)  - \boldsymbol{a}_l \\  & \qquad  -\boldsymbol{w}_l(t+k)\big) + \boldsymbol{a}_l\Big)  \leq \boldsymbol{d}_{i,l}, \, i=1, \dots, m_l\bigg \}.
\end{aligned}
\end{equation}

The random rotations and translations are given by a joint probability distribution such that the sample space of this joint distribution has cardinality $J$, i.e., $\Omega_l = \{(R_l^1, \boldsymbol{w}_l^1), \dotsc, (R_l^J, \boldsymbol{w}_l^J)\}$. A random rotation and translation is picked from this set with pmf given by $p_l = [p(1), p(2), \dotsc, p(J)]^T$. For this distribution, we also define the index set $\mathcal{J} = \{1, \dotsc, J\}$.

The safe set is defined as the region outside of the polytopic obstacles

\begin{equation}\label{eq:safeset_def}
\begin{aligned}
\mathcal{S}_l(t) &= \mathbb{R}^{n_y} \backslash \mathcal{O}_l(t) \\
&= \big\{\mathbf{y(t)}\,|\, \exists i \in \{1, \dots, m_l \},  c_{i,l}^T\mathbf{y(t)} \geq d_{i,l}\big\}.
\end{aligned}
\end{equation}

For obstacle avoidance, we aim to minimize the distance to the safe set, which is given by

\begin{equation} \label{eq:zeta}
\zeta(\mathcal{S}_l(t)) = \text{dist}(\boldsymbol{y}(t), \mathcal{S}_l(t)):= \min_{\boldsymbol{z} \in \mathcal{S}_l(t)} ||\boldsymbol{y}(t) - \boldsymbol{z}||.
\end{equation}

every node=[font=]

no markers, domain=0:16, range=-2:8, samples=200,
axis lines*=center, xlabel=$\zeta$, ylabel=$p(\zeta)$,
every axis y label/.style=at=(current axis.above origin),anchor=south,
every axis x label/.style=at=(current axis.right of origin),anchor=west,
height=5cm, width=17cm,
xtick=0, ytick=,
enlargelimits=true, clip=false, axis on top,
grid = major
[fill=cyan!20, draw=none, domain=7:15] gauss(1.5,2) ;
[very thick,cyan!50!black] gauss(1.5,2);

[gray] (axis cs:5,0) – (axis cs:5,);
[gray] (axis cs:4.5,0) – (axis cs:4.5,);
[gray] (axis cs:10,0) – (axis cs:10,);

[gray] (axis cs:1,0)–(axis cs:5,0);
[below] at (axis cs:7.0, -0.1)  $\mathrm{VaR}_{1-\alpha}(\zeta)$;
[below] at (axis cs:5, -0.1)  $\mathbb{E}(\zeta)$;
[below] at (axis cs:10, -0.1)  $\mathrm{CVaR}_{1-\alpha}(\zeta)$; [below] at (axis cs:13, -0.1)  $\mathrm{EVaR}_{1-\alpha}(\zeta)$;

[yshift=2cm, latex-latex](axis cs:7,0) – node [fill=white] Probability $1-\alpha$ (axis cs:16,0);

Comparison of the mean, VaR, and CVaR for a given confidence $\alpha \in (0,1)$. The axes denote the values of the stochastic variable $\zeta$, i.e., the minimum distance to the safe set as defined in (<ref>), and with pdf $p(\zeta)$. The shaded area denotes the $\%(1-\alpha)$ of the area under $p(\zeta)$. If the goal is to minimize $\zeta$, using $\mathbb{E}(\zeta)$ as a performance measure is misleading because tail events with low probability of occurrence are ignored. VaR gives the value of $\zeta$ at the $(1-\alpha)$-tail of the distribution. But, it ignores the values of $\zeta$ with probability below $1-\alpha$.  CVaR is the average of the values of VaR with probability less than $1-\alpha$ (average of the worst-case values of $\zeta$ in the $1-\alpha$ tail of the distribution). Note that $\mathbb{E}(\zeta) \le \mathrm{VaR}_{1-\alpha}(\zeta) \le \mathrm{CVaR}_{1-\alpha}(\zeta) \le \mathrm{EVaR}_{1-\alpha}(\zeta)$. Hence,  $\mathrm{EVaR}_{1-\alpha}(\zeta)$ is a more risk-sensitive measure.

Our goal is to minimize the risk of collision with the randomly moving obstacles by evaluating the EVaR of the distance from the probabilistic safe set (which is just the complement of the obstacle set) and constraining it to be below a certain threshold, $\epsilon_l$, i.e.,

\begin{equation}\label{eq:risksafety}
\mathrm{EVaR}_{1-\alpha}\big[\zeta(\mathcal{S}_l(t))] \leq \epsilon_l, \quad  \forall l \in \mathcal{L}.
\end{equation}

The obstacle avoidance constraint (<ref>) is an EVaR safety  constraint with confidence level $\alpha$ (see Fig. 2 for an illustrative comparison with VaR, CVaR and statistical mean) and risk tolerance (also referred to as risk threshold in this paper) $\epsilon_l$ for each obstacle $l \in \mathcal{L}$.  Note that this implies that we allow the EVaR of the distance from the safe set to be at most $\epsilon_l$ with $1-\alpha$ worst realizations.
We are now ready to present the problem we are interested in solving in this paper.

Consider the discrete-time system given by (<ref>) and the randomly moving obstacles $O_l$, $l \in \mathcal{L}$, as defined in (<ref>) and (<ref>). Given an initial condition $x_0 \in \mathbb{R}^{n_x}$, a goal set $\mathcal{X}_f \subset \mathbb{R}^{n_x}$, state constraints $\mathcal{X}\subset \mathbb{R}^{n_x}$, control constraints $\mathcal{U}\subset \mathbb{R}^{n_u}$, an immediate convex cost function $r:\mathbb{R}^{n_x} \times \mathbb{R}^{n_u}\to \mathbb{R}_{\ge0}$, a horizon $K \in \mathbb{N}_{\ge 0}$, and risk tolerances $\epsilon_l$, $l \in \mathcal{L}$, for each obstacle, compute the receding horizon controller $\{u_k \}_{k=0}^{K-1}$ such that $x(K) \in \mathcal{X}_f$ while satisfying the risk-sensitive safety constraints (<ref>), i.e., the solution to the following optimization

\begin{align}
\begin{split}
\min_{u} \quad &J(x(t), \boldsymbol{u}) := \sum_{k=0}^{K-1}r(\boldsymbol{x}_k, \boldsymbol{u}_k) \quad
\end{split}\\
\begin{split}\label{eq:dyn1}
\textrm{s.t.} \quad &\boldsymbol{x}_{k+1} = A\boldsymbol{x}_k + B\boldsymbol{u}_k
\end{split}\\
\begin{split}\label{eq:dyn2}
&\boldsymbol{y}_k = C\boldsymbol{x}_k + D\boldsymbol{u}_k
\end{split}\\
\begin{split}\label{eq:stcon}
&\boldsymbol{x}_k \in \mathcal{X}, \, \boldsymbol{u}_k \in \mathcal{U},
\end{split}\\
\begin{split}\label{eq:evarcon}
&\mathrm{EVaR}_{1-\alpha}\big[\zeta(\mathcal{S}_l(t+k))\big] \leq \epsilon_l, \forall l \in \mathcal{L}, %\, \forall l \in \mathcal{L}
\end{split}\\
\begin{split}
x_K  \in \mathcal{X}_F.
\end{split}\\
\begin{split}\label{eq:ic}
&\boldsymbol{x}_0 = \boldsymbol{x}(t),
\end{split}
\end{align}

Note that although the obstacles $\mathcal{O}_l$ are assumed to be represented by convex polytopes (<ref>), the safe set $\mathcal{S}_l(t+k)$ given in (<ref>) is nonconvex. Hence, the minimum distance to $\mathcal{S}_l(t+k)$, $\zeta(\mathcal{S}_l(t+k))$, is also nonconvex. Therefore, the risk-sensitive safety constraint (<ref>) is a nonconvex constraint in the decision variable $u$, which renders optimization problem (<ref>) nonconvex as well.

The next section will reformulate (<ref>) as a cone constraint in order to obtain a convex mixed-integer relaxation of (<ref>), which yields locally optimal solutions to (<ref>). Nonetheless, every such locally optimal solutions satisfies the constraints of optimization (<ref>) including the risk-sensitive safety constraint (<ref>).

§ EVAR RECEDING HORIZON PLANNING

This section breaks down the MPC optimization into three parts. First, we rewrite the EVaR constraint in the more tractable form of a cone constraint. Second, we reformulate the nonconvex safe set as a set of disjunctive inequalities that can be relaxed using binary variables. Lastly, we add a discrete state $\psi$ that signals task completion and allows us to prove feasibility of the MPC optimization. The resulting optimization is a convex mixed-integer program.

§.§ EVaR Constraint Reformulation

We reformulate the EVaR safety constraint to a cone constraint.

Let Assumption 1 hold, then the L.H.S. of constraint (<ref>) is equivalent to

\begin{equation} \label{eq:evar_reformulation}
\begin{aligned}
\min_{s_l, v_l, z_l, h_{l,k}} \, &\eta_l - \beta_l\ln\alpha + \sum_{j=1}^{J}p_l(j)s_l(j) \qquad \\
\textrm{s.t.} \quad &  \eta_l \in \mathbb{R}, \, \beta_l \in \mathbb{R}_{\ge0}, \\
& \eta_l-v_l(j) \geq h_{l,k},& \, \forall j \in \mathcal{J},\\
&  \beta_l + z_l(j) = 0, & \, \forall j \in \mathcal{J},\\
&  \big(s_l(j), v_l(j), z_l(j)\big) \in (K_{exp})_*, &  \forall j \in \mathcal{J},\\
& \boldsymbol{y}_k +\frac{\boldsymbol{c}_{i,l}}{||\boldsymbol{c}_{i,l}||}h_{l,k} \in \mathcal{S}_l^j(t+k), & \forall l \in \mathcal{L}, j \in \mathcal{J},
\end{aligned}
\end{equation}

where $(K_{exp})_*$ is the dual of the exponential cone.

We begin by finding the distance of $\boldsymbol{y}_k$ from the safe set, given by $\zeta(\mathcal{S}_l(t+k))$. To this end, we define a set of variables $h_{l,k}\ge 0$, $l \in \mathcal{L}$ and $k=0,\ldots,K-1$ satisfying

\begin{equation}\label{eq:safeset1}
\boldsymbol{y}_k +\frac{\boldsymbol{c}_{i,l}}{||\boldsymbol{c}_{i,l}||}h_{l,k} = \boldsymbol{z}
%\in \mathcal{S}_l^j(t+k), \,~~
\end{equation}

$\forall j \in \mathcal{J}, \forall k \in \{0, \dotsc, K-1\}$ with $i \in \{1, \dotsc, m_l\}$, which is the distance from every $\boldsymbol{y}_k$ to a point $\boldsymbol{z} \in \mathcal{X}$. If $\boldsymbol{z} \in \mathcal{S}_l^j(t+k)$, then  minimizing $h_{l,k} $ gives us the line segment connecting $ \boldsymbol{y}_k$ and $\boldsymbol{z}$, which is the minimum distance to the set $ \mathcal{S}_l^j(t+k)$. Therefore, we obtain

\begin{multline}\label{eq:safeset1}
\zeta(\mathcal{S}_l(t+k)) = \min_{\boldsymbol{z} \in \mathcal{S}_l(t+k)} ||\boldsymbol{y}(t+k) - \boldsymbol{z}|| \\ = \left \lbrace\begin{matrix} \min_{h_{l,k}} & h_{l,k} \qquad \\
\textrm{s.t.}&\boldsymbol{y}_k +\frac{\boldsymbol{c}_{i,l}}{||\boldsymbol{c}_{i,l}||}h_{l,k} \in \mathcal{S}_l^j(t+k), \,~~ \forall j \in \mathcal{J}, \end{matrix} \right.
\end{multline}

and define $h_{l,k}^*$ as the solution to (<ref>).

Next, substitute the dual form of EVaR from (<ref>) into the L.H.S. of (<ref>). Then replace the risk envelope $\mathfrak{D}$ with the exponential cone representation to yield a discrete probability distribution given by (<ref>). That is, $\mathrm{EVaR}_{1-\alpha}(\zeta(\mathcal{S}_l(t+k))) = \max_{Q\in\mathfrak{D}}\mathbb{E}_Q\big[\zeta(\mathcal{S}_l(t+k))\big]=\max_{Q\in\mathfrak{D}}\mathbb{E}_Q\big[h_{l,k}^*\big]$, where in the last equality we used (<ref>).   Thus, we have the following exponential cone program for computing EVaR

\begin{equation}\label{eq:dual}
\begin{aligned}
\max\limits_{q, \delta} \quad &\sum\limits_{j=1}^{J} h_{l,k}^* q(j) \\
\textrm{s.t.} \quad & \sum_{j=1}^J q(j) = 1,\\
& \sum_{j=1}^{J} \delta(j) \leq -\ln\alpha,\\
& \big((p_l(j), q(j), -\delta(j)\big) \in K_{exp},~~ \forall j \in \mathcal{J},  \\
& q(j) \in \mathbb{R}_{\ge 0}, \, \delta(j) \in \mathbb{R}, \phantom{--iiii}\forall j \in \mathcal{J}.
\end{aligned}
\end{equation}

The dual of the above maximization problem is given by [19, 33]:

\begin{equation} \label{eq:dualofdual}
\begin{aligned}
\min_{s_l, v_l, z_l} \, &\eta_l - \beta_l\ln\alpha + \sum_{j=1}^{J}p_l(j)s_l(j) \qquad \\
\textrm{s.t.} \quad & \eta_l-v_l(j) \geq h_{l,k}^*,& \, \forall j \in \mathcal{J},\\
&  \beta_l + z_l(j) = 0, & \, \forall j \in \mathcal{J},\\
&  \eta_l \in \mathbb{R}, \, \beta_l \in \mathbb{R}_{\ge0}, \\
&  \big(s_l(j), v_l(j), z_l(j)\big) \in (K_{exp})_*, &  \forall j \in \mathcal{J},
\end{aligned}
\end{equation}

where $\boldsymbol{s, v, z}$ are the dual variables.

We conclude that (<ref>) and (<ref>) are equivalent because strong duality holds by Slater's condition [10]. Slater's condition is satisfied by showing strict feasibility for (<ref>), i.e., there exists a feasible solution to (<ref>) such that the inequality constraints hold with strict inequalities. One such solution is when $q(j) = p(j), \,\delta(j) < 0, \, \forall j\in \mathcal{J}$.

Finally, substituting minimization problem (<ref>) for $h^*_{l,k}$ in  optimization (<ref>) gives  (<ref>).

with additional variables in the optimization: $\boldsymbol{s_l, v_l, z_l},h_{l,k}$

Utilizing the fact that $\mathrm{EVaR}_{1-\alpha}(\mathcal{S}_l(t+k))$ can be written as a minimization over the variables $\boldsymbol{s_l, v_l, z_l}, h_{l,k}$, we can return to our original MPC problem (<ref>) and simplify it as a one-layer optimization.

Consider the MPC optimization given by (<ref>) with confidence level $\alpha$ and risk tolerances $\epsilon_l$, $l \in \mathcal{L}$. If Assumption 1 holds, then (<ref>) is equivalent to a minimization over $\mathcal{V} = \{\boldsymbol{u,s_l,v_l,z_l},h_{l,k}\}$ given by

\begin{align}
\begin{split}
\min_{\mathcal{V}} \quad &J(\boldsymbol{x}(t), \boldsymbol{u}) := \sum_{k=0}^{K-1}r(\boldsymbol{x}_k, \boldsymbol{u}_k)
\end{split}\\
\begin{split}
\textrm{s.t.} \quad & \eta_l - \beta_l\ln\alpha + \sum_{j=1}^{J}p(j)s(j) \leq \epsilon_l \, \forall l \in \mathcal{L},
\end{split}\\
\begin{split}
& \eta_l-v_l(j) \geq  h_{l, k} \, \phantom{-----ii}\forall j \in \mathcal{J},l \in \mathcal{L},
\end{split}\\
\begin{split} \label{eq:safeset}
& \boldsymbol{y}_k +\frac{\boldsymbol{c}_{i,l}}{||\boldsymbol{c}_{i,l}||}h_{l,k} \in \mathcal{S}_l^j(t+k) \phantom{-} \forall l \in \mathcal{L},
\end{split}\\
\begin{split}
& \beta_l + z_l(j) = 0  \,\phantom{------ii}  \forall j\in \mathcal{J},l \in \mathcal{L},
\end{split}\\
\begin{split}
& \eta_l \in \mathbb{R}, \, \beta_l, h_{l,k} \in \mathbb{R}_{\geq 0} \,\phantom{---ii} \forall l \in \mathcal{L},
\end{split}\\
\begin{split}
& \big(s_l(j), v_l(j), z_l(j)\big) \in (K_{exp})_* \,\,\forall j\in \mathcal{J},l \in \mathcal{L},
\end{split}\\
\begin{split}
(\ref{eq:dyn1}),(\ref{eq:dyn2}),(\ref{eq:ic}), (\ref{eq:stcon}).
\end{split}
\end{align}

We can substitute the result from Lemma <ref> in (<ref>) to get

\begin{align}
\begin{split}
\min_{\boldsymbol{u}} \quad &J(\boldsymbol{x}(t), \boldsymbol{u}) := \sum_{k=0}^{K-1}r(\boldsymbol{x}_k, \boldsymbol{u}_k) \quad \end{split}\\
\begin{split}
\textrm{s.t.} \quad & (\ref{eq:dyn1}),(\ref{eq:dyn2}),(\ref{eq:ic}), (\ref{eq:stcon}),
\end{split}\\
\begin{split}
\eqref{eq:evar_reformulation} \leq \epsilon_l, \quad l \in \mathcal{L}.
\end{split}
\end{align}

Suppose we have an optimal solution to (<ref>) given by $\boldsymbol{u}^*$. As (<ref>) is feasible, its constraints must be satisfied; this implies the inner minimization (<ref>) must also be feasible (with solution $(\boldsymbol{s_l^*, v_l^*, z_l^*}, h_{l,k}^*)$). Hence, $(\boldsymbol{u^*,s_l^*, v_l^*, z_l^*}, h_{l,k}^*)$ must also be a feasible solution to (<ref>) and give the same objective value. Conversely, consider the optimal solution to (<ref>) to be given by $(\boldsymbol{u^*,s_l^*, v_l^*, z_l^*}, h_{l,k}^*)$. This $\boldsymbol{u}^*$ must be feasible for (<ref>) and gives the same objective value. Hence, the above optimization (<ref>) is equivalent to the one-layer optimization (<ref>).

§.§ Mixed-Integer Reformulation of the MPC optimization

This subsection frames the nonconvex safe set as a set of disjunctive inequalities. These inequalities are incorporated in our optimization by introducing a set of binary variables and using the Big-M relaxation [34].

The safe set (<ref>) is defined as the region outside the obstacle $l$. Given that an obstacle has rotated and translated by $R_l(t+k)$ and $\boldsymbol{w}_l(t+k)$ between times $t$ and $t+k$, we can write the safe set at $t+k$ as the region outside $\mathcal{O}_l(t+k)$ described in (<ref>). It can equivalently be expressed as a result of the rotation and translation of the safe set itself from $t$ to $t+k$

\begin{equation}
\begin{aligned}
\mathcal{S}_l^j(t+k) & = \mathbb{R}^{n_y} \backslash \mathcal{O}_l^j(t+k) \\
& = R_l^j(t+k)\mathcal{S}_l(t) + \boldsymbol{w}_l^j(t+k).
\end{aligned}
\end{equation}

In (<ref>), $\mathcal{S}_l^j(t+k)$ is a nonconvex set. For some obstacle $l \in \mathcal{L}$, we can be rewrite (<ref>) as

\begin{equation*}
R_l^j(t+k)^{-1}\bigg(\boldsymbol{y}_k +\frac{\boldsymbol{c}_{i,l}}{||\boldsymbol{c}_{i,l}||}h_{l,k} - \boldsymbol{w}_l^j(t+k)\bigg) \in \mathcal{S}_l^j(t).
\end{equation*}

Given that the obstacles are convex polygons of the form (<ref>), we write the safe region as the union of regions outside of the halfspaces that define an obstacle as follows

\begin{multline}\label{eq:disjunction}
\bigvee_{i=1}^{m_l} \boldsymbol{c}_{i,l}^T \Bigg[ R_l^j(t+k)^{-1}\bigg(\boldsymbol{y}_k +\frac{\boldsymbol{c}_{i,l}}{||\boldsymbol{c}_{i,l}||}h_{l,k} \\- \boldsymbol{w}_l^j(t+k) - \boldsymbol{a}_l\bigg) + \boldsymbol{a}_l \Bigg] \geq d_{i,l}.
\end{multline}

The above disjunctive inequalities, however, are hard to enforce. To overcome this difficulty, we relax the constraint using a Big-M reformulation. The reformulation converts the disjunctive inequalities into a set of constraints described using binary variables, $\gamma_i \in \{0,1\}$ and a large positive constant $M$. The value of $M$ depends on the bounds on $h_{l,k}$ (determined from the size of obstacle $l$) and $\boldsymbol{y}_k$ (dependent on the state and control inputs). It can be computed using linear programming. The Big-M relaxation of (<ref>) is as follows

\begin{align}
\begin{split}\label{eq:bigM1}
\boldsymbol{c}_{i,l}^T \Bigg[ R_l^j(t+k)^{-1}&\bigg(\boldsymbol{y}_k +\frac{\boldsymbol{c}_{i,l}}{||\boldsymbol{c}_{i,l}||}h_{l,k} - \boldsymbol{w}_l^j(t+k) - \boldsymbol{a}_l\bigg)  + \boldsymbol{a}_l \Bigg] \\& \geq d_{i,l} - M\gamma_i, \qquad  \, \forall i \in \{1, \dotsc m_l \},
\end{split}\\
\begin{split}
\sum_{i = 1}^{m_l} \gamma_i \leq m_l - 1.
\end{split}
\end{align}

§.§ Task Completion

In order to steer the system to the target region in finite time, we follow the footsteps of [29] and define a new discrete state $\psi \in \{ 0, 1\}$, such that $\psi = 0$ implies that the task has been completed at an earlier step or at the current step and $\psi = 1$ means that the task has not yet been completed. The update equation of $\psi$ is then given by

\begin{equation}\label{eq:taskCompletion}
\psi_{k+1} = \psi_{k} - \mu_{k},
\end{equation}

where ${\mu}_{k} \in \{ 0, 1\}$ is a discrete input.
The goal to drive the system to $(\boldsymbol{x}_{des}, \boldsymbol{u}_{des})$ (this desired position can be replaced by a region as well), is incorporated in the following additional constraints

\begin{equation}\label{eq:terminalState}
\begin{aligned}
\boldsymbol{x}_{k+1} - \boldsymbol{x}_{des} &\leq \mathds{1}{M}(1 - \mu_{k}), \\
-(\boldsymbol{x}_{k+1} - \boldsymbol{x}_{des}) &\leq -\mathds{1}{M}(1 - \mu_{k}).
\end{aligned}
\end{equation}

Here ${\mu}_{k} = 1$ if the task of reaching the goal is completed at time step $t+k+1$. Equation (<ref>) implies that $\psi$ jumps from $1 \rightarrow 0$, signaling completion of the task. After the task completion, all other MPC problem constraints can be relaxed by adding the term $M(1 - \psi_{k})$ to them, i.e., any constraints of the form $C_1\boldsymbol{\nu}_k + C_2\gamma_i + C_3 \geq 0$ are modified to $C_1\boldsymbol{\nu}_k + C_2\gamma_i  + C_3 + \mathds{1}M(1 - \psi_{k})\geq 0, \, \forall i, k$ where $\boldsymbol{\nu}_k = [\boldsymbol{u_k,x_k,y_k,s_k,v_k,z_k,h_k} ,\eta_k,\beta_k] $. We also add the following terminal constraint at time $t + K$ to ensure that the task is completed

\begin{equation} \label{eq:terminalConst}
\psi_{K} = 0.
\end{equation}

Note that the discrete state $\psi$ need not be a binary variable as long as we enforce the constraint,

\begin{equation}\label{eq:discreteStateConst}
0 \leq \psi_{k} \leq 1, \quad k=1,2,\ldots,K.
\end{equation}

The MPC objective function is then modified as

\begin{equation}
\min_{\mathcal{V}} \quad J(t) := \sum_{k=0}^{K-1}\big(r(\boldsymbol{u}_k) + \psi_k\big),
\end{equation}

where $\mathcal{V} = \{\boldsymbol{u,v,h,s,z,\mu,\gamma}\}$ and $r(\boldsymbol{u}_k)$ is a convex function of the control input with $r(0) = 0$.

The MPC optimization (<ref>) has the following convex mixed integer relaxation,

\begin{align}
\begin{split}
\min_{\mathcal{V}} \quad & J(t) := \sum_{k=0}^{K-1}\big(r(\boldsymbol{u}_k) + \psi_k\big)
\end{split}\\
\textrm{s.t.} \quad \begin{split}
& \eta_{l,k} - \beta_{l,k}\ln\alpha + \sum_{j=1}^{J}p_{l}(j)s_{l,k}(j) \leq \epsilon_l  + M_k,
\end{split}\\
\begin{split}
& \eta_{l,k}-v_{l,k}(j)+ M_k \geq  h_{l, k},
\end{split}\\
\begin{split}
& \text{L.H.S.}\eqref{eq:bigM1} + M_k\geq d_{i,l} - M\gamma_{i,l,k}(j),
\end{split}\\
\begin{split}
\sum_{i = 1}^{m_l} \gamma_{i,l,k}(j) \leq m_l - 1 + M_k,
\end{split}\\
\begin{split}
& -M_k \leq \beta_{l,k} + z_{l,k}(j) \leq M_k,
\end{split}\\
\begin{split}
& \beta_{l,k} +M_k,  \,h_{l,k}+ M_k \in \mathbb{R}_{\ge 0},
\end{split}\\
\begin{split}
& \big(s_{l,k}(j) + M_k, v_{l,k}(j), z_{l,k}(j)\big) \in (K_{exp})_*,
\end{split}\\
\begin{split}
(\boldsymbol{x}_k, \boldsymbol{u}_k, \psi_k) \in (\mathcal{X}, \mathcal{U}, 1) \cup (\mathbb{R}^n, \mathbb{R}^n, 0),
\end{split}\\
\begin{split}
& (\ref{eq:taskCompletion}), (\ref{eq:terminalState}), (\ref{eq:terminalConst}), (\ref{eq:discreteStateConst}), (\ref{eq:dyn1}),(\ref{eq:dyn2}),(\ref{eq:ic}),
\end{split}
\end{align}

where $M_k = M(1 - \psi_{k})$ and $\mathcal{X}, \mathcal{U}$ are assumed to be convex sets. The constraints must hold $\forall k\in \{1, \dotsc, K-1\},\, l\in \mathcal{L},\, j \in \mathcal{J},$ and $ i \in \{1, \dotsc, m_l\}$.

The above convex, mixed-integer relaxation of a nonconvex optimization problem will give us locally optimal solutions.

If the optimization (<ref>) is feasible at time $t = 0$, it is feasible for future time steps.

Assume that the feasible solution to (<ref>) at time $t$ is given by the input sequence $\{\boldsymbol{u}_0^*, \boldsymbol{u}_1^*, \dotsc , \boldsymbol{u}_{K-1}^*\}$ and the state sequence $\{(\boldsymbol{x}_0^*, \psi_0^*), (\boldsymbol{x}_1^*, \psi_1^*), \dotsc , (\boldsymbol{x}_K^*, \psi_K^*)\}$. Recall from the notation section that $\boldsymbol{x}_k := \boldsymbol{x}(t+k|t)$. Applying the first control input leads the system to the next state in the sequence $(\boldsymbol{x}_1^*, \psi_1^*)$, provided that the model of the actual system matches the one in (<ref>). The optimization is feasible at time $t+1$ if there exists a feasible input at time $t+K$ that does not violate constraints. Since $\psi_K^* = 0$ by virtue of the terminal constraint, all the constraints in the optimization are relaxed thereafter. Note that the state $\psi_K = 0$ is invariant due to (<ref>) and (<ref>) and that $\mu_k = 0$ at all time after the task has been completed. Therefore, once the optimization constraints are relaxed, they will remain this way.

A control input $u_K = 0$ ensures that the optimization is feasible. At time $t+1$, a feasible solution to (<ref>) is given by the control sequence $\{\boldsymbol{u}_1^*, \boldsymbol{u}_2^*, \dotsc , \boldsymbol{u}_{K-1}^*, 0\}$ and the state sequence $\{(\boldsymbol{x}_1^*, \psi_1^*), \dotsc ,(\boldsymbol{x}_K^*, 0), (A\boldsymbol{x}_{K}^*, 0)\}$. Hence, if the optimization (<ref>) is feasible at time $t$, then it is feasible at time $t+1$. By extension, if the optimization (<ref>) is feasible at time $t = 0$, it is feasible for all future time steps.

§ WAYPOINT FOLLOWING ALGORITHM

MPC is often used as a tool to plan trajectories locally and it is given a reference trajectory or a set of waypoints from a higher-level global planner like A* or RRT  [22, 14]. Let  $\{\boldsymbol{w}_1, \boldsymbol{w}_2, \dotsc, \boldsymbol{w}_N\}$ be a given a sequence of waypoints.

We call a waypoint $\boldsymbol{w}_{j+1}$ K-step reachable from $\boldsymbol{w}_j$, if there exists a feasible solution to (<ref>) with $\boldsymbol{x}_0 = \boldsymbol{w}_j$ and $\boldsymbol{x}_{K} = \boldsymbol{w}_{j+1}$.

Follow waypoints

Number of waypoints visited, $W = 0$
$W < N$
Initialize $(\boldsymbol{x}_0, \psi_0) = (\boldsymbol{w}_{W}, 1)$
Set desired goal $\boldsymbol{x}_{des} = \boldsymbol{w}_{W+1}$
$\psi_0 \neq 0$
Solve (<ref>) to obtain policy $\{\boldsymbol{u}_0^*, \boldsymbol{u}_1^*, \dotsc , \boldsymbol{u}_{K-1}^*\}$
Update $\boldsymbol{x}_0 = A\boldsymbol{x}_0 + B\boldsymbol{u}_0^*$
Update $\psi_0 = \psi_0 - \mu_0$
$\boldsymbol{x}_0 =\boldsymbol{x}_{des}$
$W = W + 1$

Assuming that the waypoint $\boldsymbol{w}_{j+1}$ is K-step reachable from $\boldsymbol{w}_j, \, \forall j \in \{1, \dotsc, N-1\}$, Algorithm <ref> gives a sequence of control inputs to move from $\boldsymbol{w}_0$ to $\boldsymbol{w}_N$ in finite time.

Consider the simple case of starting from $\boldsymbol{w}_{0}$ and reaching $\boldsymbol{w}_{1}$, i.e., when we have exactly two waypoints. We implement Algorithm <ref> till $\psi_0 = 0$. Let $J_t^*$ be the cost of the MPC optimization (<ref>) at time $t$. The feasible solution to (<ref>) at $t$ is given by the input sequence $\{\boldsymbol{u}_{0}^*, \boldsymbol{u}_{1}^*, \dotsc , \boldsymbol{u}_{K-1}^*\}$ and the state sequence $\{(\boldsymbol{x}_{0}^*, \psi_{0}^*), (\boldsymbol{x}_{1}^*, \psi_{1}^*), \dotsc , (\boldsymbol{x}_{K}^*, \psi_{K}^*)\}$. At time $t+1$, the cost of the MPC optimization is $J_{t+1}^* \leq J_{t}^* - r(\boldsymbol{u}_{0}^*) - \psi_{0}^*$. This is true because we know from Proposition <ref> that at time $t+ 1$, $\{\boldsymbol{u}_{1}^*, \dotsc , \boldsymbol{u}_{K-1}^*, 0\}$ is a feasible control input with $\psi(t+K|t+1) = 0$, i.e., $J_{t}^*$ will incur no additional cost from $u(t+K|t+1) = 0$ and $\psi(t+K+1|t+1) = 0$. Since $J_{t+1}^* - J_{t}^* \leq - r(\boldsymbol{u}_{0}^*) - \psi_{0}^*$, the cost decreases by at least $1$ at each time step till the task is completed. Considering that $J_t^*$ is always positive and finite, it will take a finite number of steps to reach $J_k^* = 0, \, k\geq t$. Hence, the policy to take the system from $\boldsymbol{w}_{0}$ to $\boldsymbol{w}_{1}$ will be implemented in finite time.

By induction, the system will reach $\boldsymbol{w}_N$ from $\boldsymbol{w}_0$ in finite time.

§ NUMERICAL RESULTS

This section shows the efficacy of the proposed EVaR-based risk-sensitive planning method via two numerical examples.

§.§ Monte-Carlo Simulations

To illustrate the effectiveness of the proposed method, we compare it to another risk measure, CVaR for different confidence levels, $\alpha$.
We look at the two-dimensional discrete system $ x_{k+1} = Ax_k + Bu_k$, with

\begin{equation*}{\small
A = \begin{bmatrix} 1.0475 & -0.0463 \\ 0.0463 & 0.9690\end{bmatrix}, \, B = \begin{bmatrix} 0.028 \\ -0.0195\end{bmatrix}.}
\end{equation*}

The control constraints are
$$-100 \leq u_k \leq 100.$$
One randomly moving obstacle interferes with the original MPC solution path that would be found in the absence of obstacles. We ran $100$ Monte-Carlo simulations for the two risk measures for different values of the confidence level $\alpha$.

The initial system state, $\boldsymbol{x}_0$, lies somewhere between $(3.1, 0.5)^T$ and $(4.1, 1.5)^T$. For each Monte-Carlo simulation, we randomly chose an initial condition in this range. The paths resulting from this set of initial conditions are most affected by the randomly moving obstacle present at $(-1, 4.5)^T$ with probability $0.75$ and at $(2.5, 3.5)^T$ with probability $0.25$. The risk tolerance is set to $\epsilon = 0.04$. The results are summarized in Table <ref>. Note that the percentage of collisions is not exact as we run 100 random simulations. A few such trajectories are shown in Fig. <ref> when $\alpha = 0.5$. The two rectangles show both possible obstacle configurations. The darker rectangle has a higher probability of occurrence, $0.75$ and the lighter rectangle has a lower probability of occurrence, $0.25$. The $20$ trajectories seen in the plots are randomly initialized as discussed above. We can see that  more CVaR trajectories intersect the obstacle.

$\alpha$      $0.9$      $0.7$     $0.5$     $0.3$     $0.1$
EVaR Collisions     $\%0$      $\%0$     $\%6$     $\%3$    $\%66$

CVaR Collisions     $\%7$      $\%17$     $\%17$     $\%14$     $\%74$

Summary of results from Monte-Carlo simulations

-0.1 true in

-0.2 true in

-0.15 true in
Comparison between EVaR (top) and CVaR (bottom) trajectories. The two uncertain obstacle locations are shown by gray rectangles.  Note that a path may cut through an obstacle between discrete time-steps (which are denoted by a '*').  A collision is not
defined when the system trajectory lies  outside the obstacle at the discrete time-steps.

§.§ Quadcopter

In this subsection, we implement the waypoint following algorithm - Algorithm <ref>. To this end, we consider a quadcopter that must follow given waypoints while avoiding randomly moving obstacles and adhering to state and control constraints. The quadcopter is described by the states $(x, y, z, \phi, \theta, \varphi, \dot{x}, \dot{y}, \dot{z}, \dot{\phi}, \dot{\theta}, \dot{\varphi})^T$. The position of the quadcopter in 3D space is given by the coordinates $x, y, z$ and the roll, pitch, and yaw are given by $\phi, \theta, \varphi$ respectively. The model of the quadcopter is given by

\[ \ddot{x} = -g\theta, \,
\ddot{y} = g\theta,\, \ddot{z}=-\frac{u_1}{m} - g, \, \]

\[ \ddot{\phi} = \frac{u_2}{I_{xx}}, \,
\ddot{\theta} = \frac{u_3}{I_{yy}}, \,
\ddot{\varphi} = \frac{u_4}{I_{zz}},\]

where $m$ is the quadcopter's mass, $g$ is the acceleration due to gravity, and $I_{xx}, I_{yy}, I_{zz}$ are the quadcopter moments of inertia about the $x,y,z$-axes of the system. The control inputs are given by $u_1, u_2, u_3, u_4$. We used the following parameters: $m = 0.65$kg, $I_{xx} = 0.0075$kg.m$^2$, $I_{yy} = 0.0075$kg.m$^2$, $I_{zz} = 0.0013$kg.m$^2$, $g = 9.81$m.s$^{-2}$ [14].

The risk constraint has two parameters: the confidence level, $\alpha$, and the risk-threshold, $\epsilon$. We chose $\alpha = 0.5,\, \epsilon = 0.04$. The waypoints are given by regions of size $[-0.5,0.5]^3$ around the waypoint center (denoted by o in Fig. <ref>). We chose a horizon length of $K = 15$ for the MPC optimization. We considered the case of having one randomly translating and rotating obstacle. The obstacle is a rectangular box of size $2$x$1$x$4$ m$^3$; it can rotate by $\frac{\pi}{2}$ and translate by $4$m along the y-axis and $1$m along the z-axis. Fig. <ref> shows all the different configurations of this obstacle at different times.

-0.1 true in

Snapshots of the quadcopter trajectory followed using Algorithm <ref> (x) given a sequence of waypoints (denoted by o). We can see the obstacle (gray boxes) changes position and orientation with time. The more likely obstacle configuration are shaded darker than the ones less likely to occur. We assume $z = 0$ is the ground plane and gravity pulls the quadrotor in the $-z$ direction.

§ CONCLUSIONS

In this paper, we proposed a risk-constrained motion planning framework for obstacle avoidance. We presented an MPC reformulation in the form of a convex mixed integer program. We showed the recursive feasibility of this optimization and introduced an algorithm to follow waypoints. As shown in this paper, through a comparison with CVaR, the framework is amenable to other risk measures. All coherent risk measures have a convex, bounded, and closed risk envelope. This framework allows any coherent risk measure constrained motion planning problem to be expressed as a convex mixed integer relaxation.

There are many paths of future research for this problem. We could extend this framework to include robust (and by extension risk-sensitive) feasibility to disturbances via constraint tightening [29]. It is also possible to extend this framework to continuous probability distributions using the relaxation technique that involves sample average approximation [14]. Future work also considers risk-sensitive robot planning under imperfect information [1].

[1]
M. Ahmadi, M. Ono, M. D. Ingham, R. M. Murray, and A. D. Ames.
Risk-averse planning under uncertainty.
In 2020 American Control Conference (ACC), pages 3305–3312.
IEEE, 2020.

[2]
M. Ahmadi, X. Xiong, and A. D Ames.
Risk-sensitive path planning via CVaR barrier functions:
Application to bipedal locomotion.
arXiv preprint arXiv:2011.01578, 2020.

[3]
A. Ahmadi-Javid.
Addendum to: Entropic value-at-risk: A new coherent risk measure.
Journal of Optimization Theory and Applications,
155(3):1124–1128, 2012.

[4]
A. Ahmadi-Javid.
Entropic value-at-risk: A new coherent risk measure.
Journal of Optimization Theory and Applications,
155(3):1105–1123, 2012.

[5]
A. Ahmadi-Javid and M. Fallah-Tafti.
Portfolio optimization with entropic value-at-risk.
European Journal of Operational Research, 279(1):225–241,

[6]
A. Ahmadi-Javid and A. Pichler.
An analytical study of norms and banach spaces induced by the
entropic value-at-risk.
Mathematics and Financial Economics, 11(4):527–550, 2017.

[7]
P. Artzner, F. Delbaen, J. Eber, and D. Heath.
Coherent measures of risk.
Mathematical finance, 9(3):203–228, 1999.

[8]
N. Bäuerle and J. Ott.
Markov decision processes with average-value-at-risk criteria.
Mathematical Methods of Operations Research, 74(3):361–379,

[9]
A. Bouman, M. F. Ginting, N. Alatur, M. Palieri, D. D. Fan,
T. Touma, T. Pailevanian, S. K. Kim, K. Otsu, J. Burdick, and A. A.
Autonomous spot: Long-range autonomous exploration of extreme
environments with legged locomotion.
In 2020 IEEE/RSJ International Conference on Intelligent Robots
and Systems (IROS), pages 2518–2525, 2020.

[10]
Stephen Boyd and Lieven Vandenberghe.
Convex Optimization.
Cambridge University Press, USA, 2004.

[11]
M. P. Chapman, J. Lacotte, A. Tamar, D. Lee, K. M. Smith, V. Cheng, J. F.
Fisac, S. Jha, M. Pavone, and C. J. Tomlin.
A risk-sensitive finite-time reachability approach for safety of
stochastic dynamic systems.
In 2019 American Control Conference (ACC), pages 2958–2963.
IEEE, 2019.

[12]
Y. Chow and M. Ghavamzadeh.
Algorithms for CVaR optimization in MDPs.
In Advances in neural information processing systems, pages
3509–3517, 2014.

[13]
Y. Chow, A. Tamar, S. Mannor, and M. Pavone.
Risk-sensitive and robust decision-making: a CVaR optimization
In Advances in Neural Information Processing Systems, pages
1522–1530, 2015.

[14]
A. Hakobyan, G. C. Kim, and I. Yang.
Risk-aware motion planning and control using cvar-constrained
IEEE Robotics and Automation Letters, 4(4):3924–3931, 2019.

[15]
A. Hakobyan, Gyeong C. Kim, and I. Yang.
Risk-aware motion planning and control using CVaR-constrained
IEEE Robotics and Automation Letters, 4(4):3924–3931, 2019.

[16]
A. Hakobyan and I. Yang.
Wasserstein distributionally robust motion control for collision
avoidance using conditional value-at-risk.
arXiv preprint arXiv:2001.04727, 2020.

[17]
A. Husain, H. Jones, B. Kannan, U. Wong, T. Pimentel, S. Tang, S. Daftry,
S. Huber, and W. L. Whittaker.
Mapping planetary caves with an autonomous, heterogeneous robot team.
In 2013 IEEE Aerospace Conference, pages 1–13. IEEE, 2013.

[18]
S. Karaman and E. Frazzoli.
Sampling-based algorithms for optimal motion planning.
The international journal of robotics research, 30(7):846–894,

[19]
B. Kocuk.
Conic reformulations for Kullback-Leibler divergence constrained
distributionally robust optimization and applications.
arXiv e-prints, 2020.

[20]
S. Koenig and R. G. Simmons.
Risk-sensitive planning with probabilistic decision graphs.
In Principles of Knowledge Representation and Reasoning, pages
363–373. Elsevier, 1994.

[21]
H. Kolvenbach, D. Wisth, R. Buchanan, G. Valsecchi, R. Grandia, M. Fallon, and
M. Hutter.
Towards autonomous inspection of concrete deterioration in sewers
with legged robots.
Journal of Field Robotics, 2020.

[22]
B. T. Lopez and J. P. How.
Aggressive collision avoidance with limited field-of-view sensing.
In 2017 IEEE/RSJ International Conference on Intelligent Robots
and Systems (IROS), pages 1358–1365, 2017.

[23]
A. Majumdar and M. Pavone.
How should a robot assess risk? towards an axiomatic theory of risk
in robotics.
In Robotics Research, pages 75–84. Springer, 2020.

[24]
S. S. Mansouri, C. Kanellakis, B. Lindqvist, F. Pourkamali-Anaraki, A. A.
Agha-Mohammadi, J. Burdick, and G. Nikolakopoulos.
A unified NMPC scheme for MAVs navigation with 3D collision
avoidance under position uncertainty.
IEEE Robotics and Automation Letters, 5(4):5740–5747, 2020.

[25]
K. Nagatani, S. Kiribayashi, Y. Okada, K. Otake, K. Yoshida, S. Tadokoro,
T. Nishimura, T. Yoshida, E. Koyanagi, M. Fukushima, et al.
Emergency response to the nuclear accident at the fukushima daiichi
nuclear power plants using mobile rescue robots.
Journal of Field Robotics, 30(1):44–63, 2013.

[26]
M. Ono, M. Pavone, Y. Kuwata, and J. Balaram.
Chance-constrained dynamic programming with application to risk-aware
robotic space exploration.
Autonomous Robots, 39(4):555–571, 2015.

[27]
C. Papachristos, S. Khattak, F. Mascarich, and K. Alexis.
Autonomous navigation and mapping in underground mines using aerial
In 2019 IEEE Aerospace Conference, pages 1–8. IEEE, 2019.

[28]
L. Prashanth.
Policy gradients for CVaR-constrained MDPs.
In International Conference on Algorithmic Learning Theory,
pages 155–169. Springer, 2014.

[29]
A. Richards and J. P. How.
Model predictive control of vehicle maneuvers with guaranteed
completion time and robust feasibility.
In Proceedings of the 2003 American Control Conference, 2003.,
volume 5, pages 4034–4040 vol.5, 2003.

[30]
T. Rouček, M. Pecka, P. Čížek,
T. Petříček, J. Bayer, et al.
DARPA subterranean challenge: Multi-robotic exploration of
underground environments.
In International Conference on Modelling and Simulation for
Autonomous Systesm, pages 274–290. Springer, 2019.

[31]
E. Seraj and M. Gombolay.
Coordinated control of UAVs for human-centered active sensing of
In 2020 American Control Conference (ACC), pages 1845–1852.
IEEE, 2020.

[32]
S. Singh, Y. Chow, A. Majumdar, and M. Pavone.
A framework for time-consistent, risk-sensitive model predictive
control: Theory and algorithms.
IEEE Transactions on Automatic Control, 2018.

[33]
P. Sopasakis, M. Schuurmans, and P. Patrinos.
Risk-averse risk-constrained optimal control.
In 2019 18th European Control Conference (ECC), pages 375–380,

[34]
A. Vecchietti, S. Lee, and I.E. Grossman.
Modeling of Discrete/Continuous Optimization Problems:
Characterization and Formulation of Disjunctions and Their
Computers and Chemical Engineering, 27(3):433–448, 2003.

[35]
A. Wang, A. M Jasour, and B. Williams.
Non-Gaussian chance-constrained trajectory planning for autonomous
vehicles under agent uncertainty.
IEEE Robotics and Automation Letters, 2020.

[36]
H. Xu and S. Mannor.
Distributionally robust Markov decision processes.
In Advances in Neural Information Processing Systems, pages
2505–2513, 2010.