A Path Planning Framework for a Flying Robot in Close Proximity of Humans

TL;DR

This paper introduces a path planning framework for flying robots that considers human safety perception by estimating model parameters through a Hidden Markov Model and generating safe trajectories using Bernstein polynomials.

Contribution

It proposes a novel HMM-based approach for learning human safety perception and integrates it into an optimal control framework for safe robot navigation near humans.

Findings

HMM approach improves model likelihood over standard regression.

Generated paths maintain safe distance from humans.

HMM-based planning shows robustness to unknown factors.

Abstract

We present a path planning framework that takes into account the human's safety perception in the presence of a flying robot. The framework addresses two objectives: (i) estimation of the uncertain parameters of the proposed safety perception model based on test data collected using Virtual Reality (VR) testbed, and (ii) offline optimal control computation using the estimated safety perception model. Due to the unknown factors in the human tests data, it is not suitable to use standard regression techniques that minimize the mean squared error (MSE). We propose to use a Hidden Markov model (HMM) approach where human's attention is considered as a hidden state to infer whether the data samples are relevant to learn the safety perception model. The HMM approach improved log-likelihood over the standard least squares solution. For path planning, we use Bernstein polynomials for discretization, as the resulting path remains within the convex hull of the control points, providing guarantees for deconfliction with obstacles at low computational cost. An example of optimal trajectory generation using the learned human model is presented. The optimal trajectory generated using the proposed model results in reasonable safety distance from the human. In contrast, the paths generated using the standard regression model have undesirable shapes due to overfitting. The example demonstrates that the HMM approach has robustness to the unknown factors compared to the standard MSE model.

Equations

$$z_{n}:=\left\{\begin{array}[]{@{}ll@{}}1,&\text{if the human is attentive to the robot,}\\ 2,&\text{otherwise}.\end{array}\right.$$

$${\pi }_{n+1}={\pi }_n\mathbf{A}.$$

$$y_n=\text{mathds}{1}_{\{z_n=1\}}(f_{\beta }(x_n)+ϵ)+\text{mathds}{1}_{\{z_n=2\}}\delta ,$$

$$p(\delta |w_n=k)=\mathcal{N}(\delta |{\mu }_k,{\sigma }_k^2),p(w_n=k)=c_k,$$

$$\underset{\theta }{\mathrm{argmax}}p(\mathbf{y}|\mathbf{x},\theta )=\underset{\theta }{\mathrm{argmax}}\sum _{\mathbf{z}}\sum _{\mathbf{w}}p(\mathbf{y},\mathbf{z},\mathbf{w}|\mathbf{x},\theta ),$$

$$p(\mathbf{y},\mathbf{z},\mathbf{w}|\mathbf{x},\theta )=$$

$$p(z_1|{\pi }_1)\left[\prod _{n=2}^{N}p(z_n|z_{n-1},\mathbf{A})\right]\prod _{n=1}^{N}p(w_n|\theta )p(y_n|z_n,w_n,x_n,\theta ),$$

$${\beta }^{\ast }:=\underset{\beta }{\mathrm{argmin}}\sum _{n=1}^{N}P(z_{n,1}|\mathbf{x},\mathbf{y},{\theta }^{old})(y_n-f_{\beta }(x_n){)}^2,$$

$$Q(\theta ,{\theta }^{old}):=\sum _{\mathbf{z},\mathbf{w}}p(\mathbf{z},\mathbf{w}|\mathbf{x},\mathbf{y},{\theta }^{old})\log p(\mathbf{y},\mathbf{z},\mathbf{w}|\mathbf{x},\theta ),$$

$${\theta }^{\ast }=\underset{\theta }{\mathrm{argmax}}Q(\theta ,{\theta }^{old}),$$

$$H_0:\phi \in {\Theta }_0\text{versus}{H}_1:\phi \in {\Theta }_1,$$

$$\lambda =2\log \left(\frac{{\sup }_{\theta \in {\Theta }_1}\mathcal{L}(\theta )}{{\sup }_{\theta \in {\Theta }_0}\mathcal{L}(\theta )}\right),$$

$$\mathbf{p}(t)=\sum _{k=0}^n{\bar{\mathbf{p}}}_k{b}_k^n(\zeta (t)),\zeta :[0,t_f]\to [0,1],\zeta (t):=\frac{t}{t_f},$$

$$\mathbf{p}(t)=\sum _{k=0}^n{\mathbf{p}}_k{\ell }_k(\eta (t)),\eta :[0,t_f]\to [-1,1],\eta (t):=\frac{2t}{t_f}-1,$$

$$\underset{{\bar{\mathbf{p}}}_0,\dots ,{\bar{\mathbf{p}}}_n,t_f}{\mathrm{argmin}}J({\bar{\mathbf{p}}}_0,\dots ,{\bar{\mathbf{p}}}_n,t_f),$$

$$\text{collision avoidance constraint},$$

$$\text{velocity and acceleration constraint},$$

$$L(\mathbf{p}(t),\dot{\mathbf{p}}(t)):=1+\gamma \max (0,f_{\beta }(\mathbf{x}(t))-b_a{)}^2,$$

$$J({\bar{\mathbf{p}}}_0,\dots ,{\bar{\mathbf{p}}}_n,t_f)=t_f+\gamma {\int }_0^{tf}\max (0,f_{\beta }(\mathbf{x}(t))-b_a{)}^{2}dt.$$

$$p(\mathbf{z},\mathbf{w}|\mathbf{x},\mathbf{y},{\theta }^{old})=p(\mathbf{z}|\mathbf{x},\mathbf{y},{\theta }^{old})p(\mathbf{w}|\mathbf{x},\mathbf{y},{\theta }^{old}).$$

$$a(z_{n,i})$$

$$b(z_{n,i})$$

$$a(z_{n,i})$$

$$b(z_{n,i})$$

$$p(y_n|z_{n,i},\mathbf{x},\theta )$$

$$\displaystyle=\left\{\begin{array}[]{@{}ll@{}}\mathcal{N}(y_{n}-f_{\beta}(x_{n})|0,\,\sigma^{2}),&\text{if}\ i=1\\ \sum_{k=1}^{K}\phi_{k}\mathcal{N}(y_{n}|\mu_{k},\,\sigma_{k}^{2}),&\text{if}\ i=2.\end{array}\right.$$

$$p(z_n|\mathbf{x},\mathbf{y},{\theta }^{old})=\frac{a(z_n)b(z_n)}{p(\mathbf{y}|\mathbf{x},{\theta }^{old})},$$

$$p(z_{n-1,j},z_{n,k}|\mathbf{x},\mathbf{y},{\theta }^{old})=\frac{a(z_{n-1,j})p(y_n|z_{n,k},x_n,{\theta }^{old})A_{jk}b(z_{n,k})}{p(\mathbf{y}|\mathbf{x},{\theta }^{old})},$$

$$p(\mathbf{y}|\mathbf{x},{\theta }^{old})=\sum _{k=1}^{2}a(z_{N,k}).$$

$$p(w_{n,i}|x_n,y_n,{\theta }^{old})=\frac{{\varphi }_i\mathcal{N}(y_n|{\mu }_i,{\sigma }_i^2)}{{\sum }_{k=1}^K{\varphi }_k\mathcal{N}(y_n|{\mu }_k,{\sigma }_k^2)},$$

$$Q(\theta ,{\theta }^{old})$$

$$:=\sum _{z,w}p(\mathbf{z},\mathbf{w}|\mathbf{x},\mathbf{y},{\theta }^{old})\log p(\mathbf{y},\mathbf{z},\mathbf{w}|\mathbf{x},\theta )$$

$$=\sum _{i=1}^{2}p(z_{1,i}|\mathbf{x},\mathbf{y},{\theta }^{old})\log {\pi }_{1,i}$$

$$+$$

Full text

I INTRODUCTION

In the last decade, multi-rotor copters have seen immense growth in popularity, not only as a research platform, but also as a commercial and industrial device. By 2020, the market for these devices is expected to attain a value of $11.2 billion with an annual growth of over 30% [1]. The mechanical simplicity, the ability to hover and the maneuverability of these flying robots justify their use in civilian applications such as media production, inspection, and precision agriculture. The inclusion of these micro unmanned aerial vehicles (UAVs) in our day-to-day lives brings immediate benefits to society. As an example, by using fast and cheap UAVs, delivery from major retailers like Amazon and Walmart, can keep a reduced inventory resulting in cost-effective warehouse management. Additionally, by exploring their small and lightweight form factor and substantially leveraging their agility and reliability, applications in elderly care, medicine, transportation and mobile surveillance are being developed. In all these examples, it is important to fly safely near people and navigate in densely populated areas. Unlike current mobile robots that autonomously operate without considering humans, these flying collocated and cooperative robots (co-robots) are intended to interact and cooperate with people in a shared and constrained environment.

It is a long tradition in robot control and motion planning to focus on the robot’s actual safety, i.e., the ability to generate safe paths that avoid collisions with obstacles. However, this is insufficient for robots operating in human congested areas. Studies of human perception have shown that there is a sharp distinction between human perceived safety and the actual safety. This paper presents a path planning framework that takes into account the human’s safety perception in the presence of a flying robot. Human’s safety perception is predicted based on data collected from physiological experiments in a virtual reality (VR) environment. In the VR experiment, the participants experience a robot flying in the proximity, and the physiological signals and the position coordinates of the flying robot are recorded simultaneously. The target variable is the physiological arousal signal, and the feature variables are the position and the velocity coordinates of the flying robot.

There are a number of unknown factors present in the data collected from these experiments. Naively assuming the unknown factor to be an independent identically distributed (i.i.d.) Gaussian noise model would not be suitable for the data where only partial observation of the state of the system is provided. However, the i.i.d. Gaussian noise assumption is popularly used for regression tasks, since maximizing the likelihood for estimation task conveniently reduces to the mean squared error (MSE) minimization problem. In [2], a recurrent neural network (RNN) is applied to predict music mood (valence), where the coefficient of determination111The coefficient of determination $R^{2}$ is a performance metric used in regression task. $R^{2}=1-\frac{\sum_{i}(y_{i}-\hat{y}_{i})^{2}}{\sum_{i}(y_{i}-\bar{y})^{2}}$, where $\bar{y}$ denotes empirical mean and $\hat{y}_{i}$ denotes the prediction of the target variable given a feature variable $x_{i}$. with a recurrent neural network (RNN) is at most 50% [2]. The undesired goodness of fit despite the highly complex models used in the previous papers suggests that the other factors not contained in the data may influence the outcome of the human test. To overcome the issue, we propose to use a Hidden Markov Model [3] approach, which divides the data samples into two clusters: (i) relevant samples where the target variable (physiological arousal) can be predicted as a function of the feature (robot’s position and velocity); (ii) irrelevant samples where it is better predicted by a random source than a function of the feature. The prediction model estimated with the help of the relevant samples is incorporated in the optimal path planning framework to take into account the human’s safety perception. In the optimal path planning, the flight path is parameterized using Bernstein polynomials, which ensures that the path remains inside the convex hull of its control points [4]. This feature helps to have collision avoidance guarantees at low computational cost.

I-A BACKGROUND

Human’s perception of a flying robot dependent upon its spatial and temporal variables (e.g. distance and speed) has been studied using virtual reality testbed, and as well using a real flying robot. A comfortable distance for a flying robot to approach a human was studied in [5], where the authors tested the effects of the size of the robot on the comfort levels of the human subjects using behavioral, physiological, and survey measures. Distancing with the robot and interaction preference between two differently behaving robots were investigated for speed and repeating behavior (cyclicity) using VR experiments, [6].

On the other hand, various design approaches to the human-aerial-robot system have been explored again for the purpose of ensuring comfort for humans. Laban effort222A method to interpret human motion used in choreography[7] is employed to design affective locomotion for a flying robot in [8], and the effect on arousal and valence due to the design parameter of the locomotion is tested. Emotional encoding in a flight path of a robot was investigated in [9], where the encoding is derived from characterizing stereotypes of personality and motion parameters using interaction vocabulary. In [10], the authors propose a flight path design approach, which improves the ease of human’s perception of the robot’s motion, and the proposed design is tested using survey measures. A signaling device that resembles the turn signals of a car is proposed to communicate the robot’s intent to humans, [11]. Using gestures to communicate the user’s intent to the robot is investigated in [12].

In the papers cited above, the focus is on either discovering a general model in human-aerial-robot interaction based on the empirical data or devising a heuristic method to improve the acceptability of the robot for humans. However, the sparse and qualitative model stated in the form of null hypothesis testing is not straightforward to apply for engineers, who intend to use it with an optimal design technique. In the proposed framework, we estimate the uncertain parameters of a human’s physiological signal model, then the estimated parametric model is considered as a cost to minimize in the optimal path planning task.

The remainder of the paper is organized as follows. In Section II, the VR experiment set-up is described, and the preliminary finding that relates the physiological arousal to the perceived safety is introduced. In Section III, we propose a model to address the influence of unknown factors and validate it against the VR experimental data. In Section IV, the optimal path planning that takes into account the human arousal model is presented. Section V summarizes and discusses future directions.

II VR EXPERIMENT AND DATASET GENERATION

Virtual Reality offers a safe, low-cost, and time efficient method to collect data [6]. For example, the precise coordinates of the human and robot can easily be recorded in real-time, which is useful for studying spatial-temporal variables in human behavior. To this end, we have developed a VR test environment to explore human-aerial-robot interactions in a variety of experimental scenarios [13, 14]. Concurrent psychophysiological reactions of participants are recorded in terms of head motion kinematics and electrodermal activity (EDA), and time-aligned with attributes of the robot’s flight path, e.g. velocity, altitude, and audio profile.

During the experiment, participants were introduced to the virtual environment (VE) and told that they would experience a simulation of an urban scene lasting approximately ten minutes. Participants were seated at the junction of a three-way intersection with unoccluded paths in the forward, left, and right direction. Three arbitrary trajectories conforming to the shape of the intersection were chosen and reversed, for a total of six unique trajectories (1.6 m altitude) (see Figure 2(a)). The simulation started with a 90 seconds baseline period – without any flying robot – allowing time for the EDA signal to plateau. The first robot then appeared and completed its trajectory. After a pause of 30 to 40 seconds, the next robot appeared, and the process repeated itself for the duration of the entire experiment. Figure 2(b) shows the position and velocity profiles of all these flying robots. We collected the data from 56 participants (20 males / 36 females) recruited from our university.

The skin conductance signal is preprocessed by EDA analysis package, Ledalab, to generate the phasic activation signal [15]. The EDA toolbox decomposes the skin conductance signal into phasic and tonic signal as shown in Figure 3. The phasic signal is then deconvolved to determine phasic activation; phasic activation represents an instantaneous arousal response.

To our knowledge, there has been no standard index of perceived safety in the literature. Although physiological measurements of arousal (e.g., EDA) alone are not necessarily equivalent to people’s perceived safety, several pieces of evidence suggest that the EDA measure of physiological arousal in our study is closely related to people’s anxiety induced by the approaching drone. For example, in a follow-up experiment examining the effects of path height on people’s EDA responses, we found that the EDA phasic response was significantly stronger for a drone approaching at eye-height, where a potential collision was possible, than when the drone was flying at a height beyond the observer’s head, where there is no danger of collision, even though all other characteristics of the drone movement were the same. These results suggest that such arousal was most likely due to people’s anxiety in response to approaching danger rather than general excitement caused by watching flying robots. Moreover, analysis of the simultaneous head motion showed that as the EDA signal increased with the approaching drone, people made characteristic collision-avoidance movements by jerking their heads away from the drone. These findings again suggest that the physiological arousal signals observed in our study are most likely a result of people’s anxiety to avoid impeding danger, rather than general excitement. Thus, in the following sections we consider the EDA signal as an operational approximation of the human’s perceived safety for the optimal path generation algorithm.

III THE PROPOSED MODEL

We aim to develop a data-driven model that predicts the phasic activation (arousal), given the robot’s position and velocity. Let $y_{n}\in\mathbb{R}$ denote the phasic activation, where $n$ is the time index. The input (feature) variable, denoted by $x_{n}\in\mathbb{R}^{8}$, is the vector that contains the distance to the robot, the rate of change of the distance, the Cartesian position coordinates, and the velocity coordinates. Despite the high-fidelity test environment, it is impossible to measure every stimulus on the subject, i.e. there are unknown factors in the data. As an example, one of our collected datasets, shown in Figure 4, illustrates the unknown factors’ presence in the data. One can notice an increase in the phasic activation in the shaded area, although the flying robot is far away and virtually invisible to the subject.

To account for the unknown factors in the data, we hypothesize that the unexpected spike of the phasic activation in Figure 4 is due to the change of the participant’s focus of attention, i.e. the participant is distracted by some other stimulus. Inspired by the work in [16], we model the sequential dependence of the (hidden) human’s focus of attention using a Hidden Markov Model (HMM). The HMM has two states represented by the latent variable that models the human’s attention state, which we denote by

[TABLE]

Then $z_{n}$ is modeled by a homogeneous Markov chain with the following probability transition equation:

[TABLE]

The vector $\pi_{k}:=[p(z_{n}=1),p(z_{n}=2)]$ is the stochastic row vector for the distribution over the state $z_{n}$, and $\mathbf{A}\in\mathbb{R}^{2\times 2}$ denotes the transition probability matrix of the Markov chain333We used $p(\cdot)$ for both probability and probability density; its distinction easily follows from the context.. The initial distributions $\pi_{0}$ and $\mathbf{A}$ are the parameters of the Markov chain.

The attention state variable $z_{n}$ assigns one of the two output emission models $f_{\beta}(x_{n})+\epsilon$ or an independent random source $\delta$ as follows

[TABLE]

where $\mathds{1}_{A}$ denotes the indicator function, and $f_{\beta}:\mathbb{R}^{8}\rightarrow\mathbb{R}$ is a function $f_{\beta}(x):=\beta^{\top}\phi(x)$, which is linearly parameterized with $\beta$ and basis $\phi(x)$4443rd order polynomial basis functions were chosen for the $\phi(\cdot)$., $\epsilon\sim\mathcal{N}(0,\sigma^{2})$, and $\delta$ denotes the random source. As seen in (2), $y_{n}$ depends on $x_{n}$ when $z_{n}=1$; however, $y_{n}=\delta$ when $z_{n}=2$, i.e. it is modeled as an independent random signal. In (2), it can be seen that one of the two regression functions of the models $y_{n}=f_{\beta}(x_{n})$ and $y_{n}=\delta$ is chosen based on the likelihood given the observation $(x_{n},y_{n})$.

In addition to the hypothetical binary HMM that models the change of human’s attention, we further hypothesize that the unexpected spike of the activation signal follows a multimodal distribution. We employ a mixture of Gaussians to model $\delta$ for the multimodal distribution. The Gaussian mixture model (GMM) allows a multi-modal and skewed distribution in contrast to the Gaussian distribution. The density of the mixture model for $\delta$ is defined by another latent variable $w_{n}\in\{1,\dots,K\}$ as follows:

[TABLE]

where $c_{k}$ are the mixing coefficients such that $\sum_{k=1}^{K}c_{k}=1$, $\mathcal{N}(\delta|\mu,\,\sigma)$ denotes a Gaussian density function of $\delta$ with the mean $\mu$ and the variance $\sigma^{2}$. We assume that $w_{n}$ is independent and identically distributed. Also, it is assumed that $w_{n}$ and $z_{n}$ are independent, and furthermore, $w_{n}$ and $z_{n}$ are conditionally independent, given the observation $(x_{n},y_{n})$.

III-A Model Parameter Estimation

The model defined by (1) - (3) has a set of parameters denoted by $\theta:=\{\beta,\mu,\sigma,\pi_{1},\mathbf{A},\{c_{i},\mu_{i},\sigma_{i}\}^{K}_{i=1}\}$. Given the dataset $\mathbf{x}:=\{x_{1},\dots,x_{N}\}$ and $\mathbf{y}:=\{y_{1},\dots,y_{N}\}$, the parameter of the model is estimated by the maximum likelihood estimation (MLE) through the following conditional likelihood equation:

[TABLE]

where the summation takes place over all possible sequences $\mathbf{z}:=\{z_{1},\dots,z_{N}\}$ and $\mathbf{w}:=\{w_{1},\dots,w_{N}\}$. The number of terms for the summation is $2^{N}K^{N}$, which makes the optimization intractable for large number of samples. Due to this challenge, Expectation-Maximization (EM) algorithm [17] is widely used to obtain the MLE for HMM.

EM Algorithm: Assume the complete data $(\mathbf{x},\mathbf{y},\mathbf{z},\mathbf{w})$ is available. Using the independence assumption on $w_{n}$, $z_{n}$ and the Markov chain property as defined in (1), the conditional likelihood is calculated as follows:

[TABLE]

where $\mathbf{z}:=\{z_{1},\dots,z_{N}\}$ and $\mathbf{w}:=\{w_{1},\dots,w_{N}\}$. The EM algorithm iteratively maximizes the likelihood using the posterior $p(\mathbf{z},\mathbf{w}|\mathbf{x},\mathbf{y},\theta^{old})$ and the likelihood for the supposed complete data $p(\mathbf{y},\mathbf{z},\mathbf{w}|\mathbf{x},\theta)$, as summarized in Algorithm 1. The detailed calculation of EM algorithm for the proposed model is in Appendix.

In contrast to standard regressions that minimize the mean squared error (MSE), the latent variable model (HMM) determines the parameter of the prediction model $f_{\beta}(x)$ as the weighted least squared error solution with the weight of the posterior $P(z_{n,1}|\mathbf{x},\mathbf{y},\theta^{old})$ as follows:

[TABLE]

where $P(z_{n,1}|\mathbf{x},\mathbf{y},\theta^{old})$ denotes $P(z_{n}=1|\mathbf{x},\mathbf{y},\theta^{old})$. The proposed method puts greater weight on the samples, which are more relevant to the input based on the posterior of the attention state.

III-B RESULT

We choose the following initial parameter $\theta^{0}$:

$\beta^{0}:=\operatorname*{argmin}_{\beta}\sum_{n=1}^{N}(y_{n}-f_{\beta}(x_{n}))^{2}$, and $\sigma^{0}=0.5$, 2. 2.

$\mathbf{A}^{0}:=\begin{bmatrix}1/2&1/2\\ 1/2&1/2\end{bmatrix}$, and $\pi_{1}^{0}:=[1/2,1/2]$, 3. 3.

$c_{k}^{0}=1/K$, $\sigma_{k}^{0}=1$, and $\mu_{k}$ are randomly chosen from the interval,$[-1,1]$,

where $K$ is the number of the Gaussian basis of the GMM. For the linear function $f_{\beta}(x):=\beta^{\top}\phi(x)$, we choose the basis functions $\phi(x)$ as polynomials with degree 3.

The i.i.d. Gaussian noise model (MSE minimization approach) is contained in the proposed model structure as a special case, $p(z_{n}=1)=1$, i.e. the arousal is always explained by $f_{\beta}(x)$ than the random source $\delta$. We would like to know whether the proposed model is better than the MSE minimization approach. Likelihood ratio test is typically used to decide if the additional complexity in the modeling is desired compared to a simple model (e.g. i.i.d. Gaussian noise model), [18]. Since HMM is non-identifiable in general [19], the likelihood ratio model comparison test with the training dataset does not apply. To determine which model is more suitable, we calculate the likelihood with test dataset by employing the approach from [20]. We randomly partition the data from 56 subjects into a training set with 38 subjects and a test set with the other 18 subjects. Figure 6 shows the log-likelihood with test data set for the models: (i) the Gaussian i.i.d. noise model, (ii) the proposed model with a different basis $K$ of GMM. We see that there is a significant increase in the likelihood using the proposed model as compared to the Gaussian i.i.d. noise model (or the MSE minimization model).

Remark 1

The improvement of the log-likelihood by using the proposed model suggests that the null hypothesis555The null hypothesis assumes that the true parameter $\varphi$ is in the set of parameters $\Theta_{0}$, which corresponds to the Gaussian i.i.d. noise model. is not true. By rejecting the null hypothesis, we show that the proposed model is more suitable than the Gaussian noise model. Consider the following log-likelihood ratio test (see 11.7.4 in [18]):

[TABLE]

where $\varphi$ denotes the true parameter, $\Theta_{0}$ denotes the set of parameters for the Gaussian noise model, and $\Theta_{1}$ denotes the set of parameters for the proposed model. Note that the Gaussian noise model is a special case of the proposed model, i.e. $\Theta_{0}\subseteq\Theta_{1}$. The likelihood ratio statistics $\lambda$ is calculated as

[TABLE]

where $\mathcal{L}(\theta):=p(\mathbf{y},\mathbf{z},\mathbf{w}|\mathbf{x},\theta)$ denotes the likelihood for $\theta$ in (5). Figure 6 shows that $\lambda>4000$. The relative degree is $r=10$, as the proposed model has 10 more parameters than the Gaussian noise model. The likelihood ratio test is: reject $H_{0}$, when $\lambda>\chi^{2}_{r,\alpha}$. We reject $H_{0}:\varphi\in\Theta_{0}$ with p-value at $0.01$.

We fix the proposed model with $K=2$, since a greater number of basis does not result in significant improvement in likelihood, as shown in Figure 6. The function $f_{\beta}$ with the fixed model is used to predict the phasic activation (arousal) as shown in Figure 5(a). Figure 5(b) shows that the MSE minimization approach has signs of over-fitting (oscillation and spiky shape). In the following section, optimal path planning with the prediction model is presented.

IV OPTIMAL PATH PLANNING

The ability to generate a safe path while considering human’s safety perception (using the prediction function $f_{\beta}(\cdot)$) is necessary for our optimal path planning task. To meet the requirement, we employ the trajectory generation method from [4, 21]. To ensure spatial separation of the robot’s path from obstacles, we use Bernstein polynomials to discretize trajectories. Bernstein polynomials are useful for checking collision avoidance, as the convex hull of the vertices determined from the coefficients of the Bernstein polynomial contains the flight path, Figure 7. The convex hull is used to check for collision between the vehicle and the obstacles, as demonstrated in [4].

Polynomial interpolation of state trajectories has been used to determine numerical approximate solutions to optimal path planning. For example, an interpolation polynomial with the Legendre-Gauss-Lobatto (LGL) time nodes has been used to solve an optimal control problem in [23, 24]. In [21] it is shown that using Bernstein polynomials the optimal control solution can be approximated sufficiently closely as the number of time-nodes increases, while at the same time ensuring spatial separation from obstacles. In our path planning framework, we use the LGL quadrature [23] to approximately calculate the cost to minimize, and we use Bernstein polynomials to check for collision avoidance [4].

IV-A Finite Dimensional Optimization

Both Bernstein polynomial curves and interpolation polynomial curves with LGL nodes are two equivalent parametrizations for a polynomial trajectory. Consider a 2D trajectory time function $\boldsymbol{\mathbf{p}}(t):=[x(t),y(t)]^{\top}$ as an $n^{th}$ order polynomial in $x(t)$ and $y(t)$. We present a brief overview of both equivalent representations.

A degree $n$ Bernstein polynomial is given by:

[TABLE]

where $b^{n}_{k}(\zeta):=\binom{n}{k}(1-\zeta)^{n-k}\zeta^{k},\quad\zeta\in[0,1]$, represents the polynomial basis, and the coefficients $\bar{\boldsymbol{\mathbf{p}}}_{k}$ are called control points666Note that the polynomial is a vector equation, so the coefficient $\bar{\boldsymbol{\mathbf{p}}}_{k}$ is also a vector. of the Bernstein polynomial. 2. 2.

The interpolation curve is represented as:

[TABLE]

where $\boldsymbol{\mathbf{p}}_{k}$ are interpolation points at time nodes $t_{k}$, and $\ell_{k}(\eta(t)):=\prod_{0\leq i\leq n,i\neq k}\left(\frac{\eta(t)-\eta(t_{i})}{\eta(t_{k})-\eta(t_{i})}\right)$ are the Lagrange polynomial basis.

The $n^{th}$ order polynomial trajectories can be parameterized by either $n+1$ control points $\bar{\boldsymbol{\mathbf{p}}}_{k}$, or $n+1$ interpolation points $\boldsymbol{\mathbf{p}}_{k}$, and the transformation between the control points and the interpolation points can be done using matrix multiplication.

The optimal path planning is formulated as the following finite dimensional optimization:

[TABLE]

where $J(\bar{\boldsymbol{\mathbf{p}}}_{0},\dots,\bar{\boldsymbol{\mathbf{p}}}_{n},t_{f})$ is the LGL quadrature of $\int_{0}^{t_{f}}L(\boldsymbol{\mathbf{p}}(t),\dot{\boldsymbol{\mathbf{p}}}(t))dt$ calculated by the method in [23], and $L(\boldsymbol{\mathbf{p}}(t),\dot{\boldsymbol{\mathbf{p}}}(t))$ is the running cost to be minimized in time average. Constraint equations for collision avoidance and velocity acceleration can be written as functions of control points $\bar{\boldsymbol{\mathbf{p}}}_{0},\dots,\bar{\boldsymbol{\mathbf{p}}}_{n}$ by following the methodology in [4].

IV-B Optimal Path Planning in the Presence of Humans

Define $\boldsymbol{\mathbf{x}}(t)$ in the same way as we defined $x_{n}$ in (2), where $\boldsymbol{\mathbf{x}}(t)\in\mathbb{R}^{8}$ contains the distance to the robot, the rate of change of the distance, the position coordinates and the velocity coordinates at time $t$. Notice that with the polynomial path $\boldsymbol{\mathbf{p}}(t)$ and $\dot{\boldsymbol{\mathbf{p}}}(t)$ one can directly construct $\boldsymbol{\mathbf{x}}(t)$. For this reason, to simplify the notation we can use $\boldsymbol{\mathbf{x}}(t)$ and $(\boldsymbol{\mathbf{p}}(t),\dot{\boldsymbol{\mathbf{p}}}(t))$ interchangeability as arguments of the functions $f_{\beta}(\cdot),\>J(\cdot)$ and $L(\cdot)$.

In the optimal path planning, we only consider values of $f_{\beta}$ larger than a threshold $b_{a}$, where $b_{a}~{}\geq~{}0$ is essentially a tuning parameter. Intuitively, we ignore arousal levels below the threshold. To make the optimization problem tractable, instead of adding a strict constraint to the minimization problem, the constraint is incorporated in the running cost as a penalty function [25]:

[TABLE]

where $\gamma$ is the penalty coefficient. The corresponding cost function $J(\bar{\boldsymbol{\mathbf{p}}}_{0},\dots,\bar{\boldsymbol{\mathbf{p}}}_{n},t_{f})$ becomes

[TABLE]

The two arousal prediction functions are used in the optimal flight trajectory generation, as shown in Figure 8 and Figure 9. The smaller value of $b_{a}$ results in a path that is more safety conscious, as intended by the running cost function in (6). Flight paths generated with the proposed model show the desirable behavior, as shown in Figure 5(a) (decreasing $b_{a}$ results in greater distance from the human). However, the paths with the MSE minimization model have unconvincing shapes. This undesirable behaviour of the MSE minimization model is due to over-fitting, as we have seen in Figure 5(b). It shows that the arousal prediction model which only minimizes MSE does not generalize in the optimal path generation task.

V CONCLUSION

We present a path planning framework that takes into account human’s safety perception. Psychophysiological reactions for different paths of flying robots in VR were collected to estimate the arousal prediction model. To consider the unknown factor in the data, we proposed a hidden Markov model approach. Compared to the mean squared error (MSE) minimization approach (due to i.i.d. Gaussian noise model), the proposed model has improved the likelihood significantly. When the arousal prediction functions are implemented in the optimal path planning, the flight paths with the proposed model show desirable behaviors, in contrast to the unconvincing flight paths with the MSE minimization approach.

APPENDIX

The EM algorithm for the proposed model is calculated by the following three subsequent steps:

V-1 Determine the Posterior $p(\mathbf{z},\mathbf{w}|\mathbf{x},\mathbf{y},\theta^{old})$

Using the conditional independence assumption of the latent variables, the posterior is factorized as

[TABLE]

Now, we can calculate $p(\mathbf{z}|\mathbf{x},\mathbf{y},\theta^{old})$ and $p(\mathbf{w}|\mathbf{x},\mathbf{y},\theta^{old})$ separately.

Forward-Backward Algorithm [26]:

Define $a(z_{n,i})$ and $b(z_{n,i})$ as follows:

[TABLE]

where $z_{n,i}$ denotes the event $\{z_{n}=i\}$. To calculate the posterior, the following recursive equations are used:

[TABLE]

where $p(y_{n}|z_{n,i},\mathbf{x},\theta^{old})$ is calculated using (2) as

[TABLE]

The boundary values $a(z_{1})$ and $b(z_{N})$ are determined as $a(z_{1})=p(z_{1}|\pi_{1})p(y_{1}|\mathbf{x},\theta^{old})$ and $b(z_{N})=1$. After calculating $a(z_{n})$ recursively and $b(z_{n})$, the posterior is determined as follows:

[TABLE]

[TABLE]

where the likelihood is calculated as

[TABLE]

Posterior for the GMM: Due to the conditional independence assumption of $w_{n}$ and $z_{n}$, the posterior $p(w_{n,i}|x_{n},y_{n},\theta^{old})$ is calculated independently from $p(z_{n}|\mathbf{x},\mathbf{y},\theta^{old})$, by directly using the result on GMM from [26] as follows

[TABLE]

where $w_{n,i}$ denotes the event $\{w_{n}=i\}$.

V-2 Calculate $Q(\theta,\theta^{old})$

Using the poseterior calculated in (9) and (12) and the likelihood in (5), $Q(\theta,\theta^{old})$ in Algorithm 1 is calculated by expanding the log term:

[TABLE]

where $z_{n,i}$ denotes the event $\{z_{n}=i\}$, $w_{n,k}$ denotes the event $\{w_{n}=k\}$, $A_{i,j}$ is the $(i,j)$ element of the matrix $\mathbf{A}$, and $p(y_{n}|z_{n,i},w_{n,k},x_{n},\theta)$ is calculated using the model equation (2) as follows

[TABLE]

V-3 Find the maximizer $\theta^{*}$

As seen in (13) and (14), each term has a distinct set of parameters. Hence, we can determine the maximizer for each term independently from the other terms. From [26] we have the maximizers as follows:

[TABLE]

[TABLE]

The maximizer for the last term in (13) is calculated using the model equation (2). Let $L$ denote the last term in (13). Using (14), $L$ is written as follows:

[TABLE]

The variable $\beta^{*}:=\operatorname*{argmax}_{\beta}Q(\beta,\theta^{old})$ is calculated as

[TABLE]

The $\phi_{i}^{*}$, $\mu_{i}^{*}$, $\sigma_{i}^{*}$, and $\sigma^{new}$ are determined using KKT (Karush-Kuhn-Tucker) condition as follows:

[TABLE]

[TABLE]

[TABLE]

[TABLE]