# Human-in-the-Loop Mixed-Initiative Control under Temporal Tasks

Meng Guo, Sofie Andersson and Dimos V. Dimarogonas The authors are with the
ACCESS Linnaeus Center, School of Electrical Engineering, KTH Royal Institute
of Technology, SE-100 44, Sweden (Email: mengg, sofa, dimos@kth.se). This work
is supported in part by the H2020 ERC Starting Grant BUCOPHSYS, the EU H2020
Research and Innovation Programme under GA No.731869 (Co4Robots), SSF COIN
project and the Knut and Alice Wallenberg Foundation (KAW).

###### Abstract

This paper considers the motion control and task planning problem of mobile
robots under complex high-level tasks and human initiatives. The assigned task
is specified as Linear Temporal Logic (LTL) formulas that consist of hard and
soft constraints. The human initiative influences the robot autonomy in two
explicit ways: with additive terms in the continuous controller and with
contingent task assignments. We propose an online coordination scheme that
encapsulates (i) a mixed-initiative continuous controller that ensures all-
time safety despite of possible human errors, (ii) a plan adaptation scheme
that accommodates new features discovered in the workspace and short-term
tasks assigned by the operator during run time, and (iii) an iterative inverse
reinforcement learning (IRL) algorithm that allows the robot to asymptotically
learn the human preference on the parameters during the plan synthesis. The
results are demonstrated by both realistic human-in-the-loop simulations and
experiments.

##  I Introduction

Autonomous systems are becoming increasingly prevalent in our daily life, with
examples such as self-driving vehicles, package delivery drones and household
service robots [1]. Nevertheless, these autonomous systems often perform the
intended tasks under the supervision or collaboration with human operators
[2]. On the high level, the human operator could assign tasks for the robot to
execute or monitor the task execution progress during run time. On the low
level, the operator could directly influence or even overtake the control
commands of the robot from the on-board autonomous controller, which can be
useful to guide the robot through difficult parts of the task [2, 3, 4]. On
the other hand, the autonomous controller should take into account possibly
erroneous inputs from the operator and ensure that safety constraints are
never violated. Thus, addressing properly these online interactions between
the autonomous system and the operator during the design process is essential
for the safety and efficiency of the overall system.

In this work, we consider the interactions on both levels. Particularly, on
the high level, the operator assigns (i) offline a local task as LTL formulas
for hard and soft task constraints, and (ii) online temporary tasks with
deadlines. On the low level, the operator’s control inputs is fused directly
with the autonomous controller via a mixed-initiative controller. The proposed
motion and task planning scheme ensures that the hard task constraints
regarding safety are obeyed at all time, while the soft constraints for
performance are improved gradually as the robot is guided to explore the
workspace and more importantly, learn about the human preference over the
synthesized plan.

We rely on LTL as the formal language [5] to describe complex high-level tasks
beyond the classic point-to-point navigation. Many recent papers can be found
that combine robot motion planning with model-checking-based task planning,
e.g., a single robot under LTL motion tasks [6, 7, 8, 9], a multi-robot system
under a global task [10], or a multi-robot system under local independent [11]
or dependent tasks [12, 13]. However, none of the above addresses directly the
human initiative, neither in the continuous control nor in the discrete
planning. On the other hand, human inputs are considered in [14] via GR(1)
task formulas that require the robot to be reactive to simple sensory signals
from the human. The high-level robot-human interaction is modeled as a two-
player Markov Decision Process (MDP) game in [15, 16] where they take turns to
influence the system evolution. The goal is to design a shared control policy
to satisfy a LTL formula and minimize a cost function. Another recent work
[17] addresses the control problem of MDP under LTL formulas where the
autonomy strategy is blended into the human strategy in a minimal way that
also ensures safety and performance. However, the direct interaction on the
low level is not investigated in the aforementioned work. More importantly, an
all-time involvement of the human is required for these frameworks, while we
only assume human intervention whenever preferred by the operator itself.

Furthermore, the notion of mixed-initiative controller is firstly proposed in
[4] that combines external human inputs with the traditional navigation
controller [18], while ensuring safety of the overall system. The work in [19,
20] proposes a systematic way to compose multiple control initiatives using
barrier functions. However, high-level complex temporal tasks are not
considered in these work.

The main contribution of this work lies in a novel human-in-the-loop control
framework that allows human interaction on both high level as complex tasks
and low level as continuous inputs. We ensure all-time safety during the
interaction and accommodation of short-term contingent tasks assigned during
run time. Lastly, the proposed IRL algorithm enables the robot to
asymptotically learn and adapt to the human operator’s preference in the plan
synthesis.

The rest of the paper is organized as follows: Section II introduces some
preliminaries of LTL. Section III formulates the problem. Main algorithmic
parts are presented in Section IV, which are integrated into the complete
framework in Section V. Numerical and experiment studies are shown in Sections
VI. We conclude in Section VII.

##  II Preliminaries

###  II-A Linear Temporal Logic (LTL)

A LTL formula over a set of atomic propositions $AP$ that can be evaluated as
true or false is defined inductively according to the following syntax [5]:
$\varphi::=\top\;|\;a\;|\;\varphi_{1}\wedge\varphi_{2}\;|\;\neg\varphi\;|\;\bigcirc\varphi\;|\;\varphi_{1}\,\textsf{U}\,\varphi_{2},$
where $\top\triangleq\texttt{True}$, $a\in AP$ and $\bigcirc$ (_next_), U
(_until_). There are other useful operators like $\square$ (_always_),
$\Diamond$ (_eventually_), $\Rightarrow$ (_implication_). The full semantics
and syntax of LTL are omitted here due to limited space, see e.g., [5].
Syntactically co-safe LTL (sc-LTL) is a subclass of LTL that can be fulfilled
by finite satisfying prefix [21].

###  II-B Büchi Automaton

Given a LTL formula $\varphi$, there always exists a Nondeterministic Büchi
Automaton (NBA) $\mathcal{A}_{\varphi}$ that accepts all the languages that
satisfy $\varphi$ [5]. It is defined as
$\mathcal{A}_{\varphi}=(Q,\,2^{AP},\,\delta,\,Q_{0},\,\mathcal{F})$, where $Q$
is a finite set of states; $Q_{0}\subseteq Q$ is the set of initial states,
$2^{AP}$ is the set of input alphabets; $\delta:Q\times
2^{AP}\rightarrow{2^{Q}}$ is a transition relation and $\mathcal{F}\subseteq
Q$ is a set of accepting states. There are fast translation algorithms [22] to
obtain $\mathcal{A}_{\varphi}$. Moreover, denote by
$\chi(q_{m},\,q_{n})=\\{\ell\in 2^{AP}|\,q_{n}\in\delta(q_{m},\,\ell)\\}$ the
set of all input alphabets that enable the transition from $q_{m}$ to $q_{n}$
in $\delta$. Then the distance between $\ell\in 2^{AP}$ and $\chi\subseteq
2^{AP}$($\chi\neq\emptyset$) is defined by $\texttt{Dist}(\ell,\,\chi)=0$ if
$\ell\in\chi$ and $\min_{\ell^{\prime}\in\chi}\;|\\{a\in
AP\,|\,a\in\ell,a\notin\ell^{\prime}\\}|$, otherwise. Namely, it returns the
minimal difference between $\ell$ and any element in $\chi$.

##  III Problem Formulation

###  III-A Dynamic Workspace and Motion Abstraction

The bounded workspace where the robot is deployed is denoted by
$\mathcal{W}\subset\mathbb{R}^{2}$. It consists of $N>0$ regions of interest,
denoted by $\Pi=\\{\pi_{1},\pi_{2},\cdots,\pi_{N}\\}$, where
$\pi_{n}\subset\mathcal{W}$. Furthermore, there is a set of $M>0$ properties
(atomic propositions) associated with $\Pi$, denoted by
$AP=\\{a_{0},a_{1},\cdots,a_{M}\\}$, e.g., “this is a public area”, “this is
office room one” and “this meeting room is in use”.

The robot’s motion within the workspace is abstracted as a labeled transition
system $\mathcal{T}\triangleq(\Pi,\,\rightarrow,\,\Pi_{0},\,AP,\,L)$, where
$\Pi,\,AP$ are defined above, $\rightarrow\subseteq\Pi\times\Pi$ is the
transition relation that $(\pi_{i},\pi_{j})\in\rightarrow$ if the robot can
move from region $\pi_{i}$ to region $\pi_{j}$ without crossing other regions
in $\Pi$, $\Pi_{0}\in\Pi$ is where the robot starts initially,
$L:\Pi\rightarrow 2^{AP}$ is the labeling function where $L(\pi_{i})$ returns
the set of properties satisfied by $\pi_{i}$. Since the workspace is assumed
to be only _partially-known_ and dynamic, the labeling function and the
transition relation may change over time.

###  III-B Mixed-initiative Controller

For the simplicity of discussion, we assume that the robot satisfies the
single-integrator dynamics, i.e., $\dot{x}=u$, where $x,\,u\in\mathbb{R}^{2}$
are the robot position and control inputs. For each transition
$(\pi_{s},\pi_{g})\in\rightarrow$, the robot is controlled by the mixed-
initiative navigation controller [4] below:

$u\triangleq u_{r}(x,\pi_{s},\pi_{g})+\kappa(x,\Pi)\,u_{h}(t)$ (1)

where $u_{r}(x,\pi_{s},\pi_{g})\in\mathbb{R}^{2}$ is a given autonomous
controller that navigates the robot from region $\pi_{s}$ to $\pi_{g}$, while
staying within $\mathcal{W}$ and without crossing other regions in $\Pi$; the
function $\kappa(x,\Pi)\in[0,\;1]$ is a smooth function to be designed; and
$u_{h}(t)\in\mathbb{R}^{2}$ is the human input function, which is
_uncontrollable_ and _unknown_ by the robot.

###### Remark 1

The proposed motion and task coordination scheme can be readily extended to
robotic platforms with other dynamics and different navigation controllers,
such as potential-field-based [4] and sampling-based [7, 23]. $\blacksquare$

###  III-C Robot Task Assignment

The robot is assigned by the human operator a local task as LTL formulas over
$AP$, which has the following structure:

$\varphi\triangleq\varphi_{\text{hard}}\wedge\varphi_{\text{soft}}\wedge\varphi_{\text{temp}}$
(2)

where $\varphi_{\text{soft}}$ and $\varphi_{\text{hard}}$ are “soft” and
“hard” sub-formulas that are assigned _offline_. Particularly,
$\varphi_{\text{hard}}$ includes safety constraints such as collision
avoidance: “avoid all obstacles” or power-supply guarantees: “visit the
charging station infinitely often”; $\varphi_{\text{soft}}$ contains
additional requirements for performance such as surveillance: “surveil all
bases infinitely often”. Introducing soft and hard constraints is due to the
observation that the partially-known workspace might render parts of the task
infeasible initially and thus yielding the need for them to be relaxed, while
the safety-critical parts should not be relaxed; Lastly,
$\varphi_{\text{temp}}$ contains short-term contingent tasks that are assigned
as sc-LTL formulas _online_ and unknown beforehand. The structure in (2)
provides an effective way for the operator to handle both standard operational
tasks and contingent demands.

###  III-D Control Objective

Given the abstraction model $\mathcal{T}$ and the task formula $\varphi$, the
control objective is to design function $\kappa(\cdot)$ and control input $u$
in (1) such that: (I) the hard constraints in $\varphi_{\text{hard}}$ are
always satisfied, given all possible human inputs; (II) each time a temporary
task $\varphi_{\text{temp}}$ is assigned, it is satisfied in finite time; and
(III) the satisfaction of the soft constraints in $\varphi_{\text{soft}}$
adapts to the human inputs.

##  IV Algorithmic Components

In this section, we present the four algorithmic components of the overall
solution presented in Section V. Particularly, we start from constructing a
parameterized product automaton for the plan synthesis. Then we present a
mixed-initiative controller that guarantees safety and meaningful inputs from
the operator. Furthermore, we discuss a plan adaptation algorithms for real-
time updates of the workspace model and contingent task assignment. At last,
we describe a IRL algorithm to learn about the human preference.

###  IV-A Initial Discrete Plan Synthesis

Denote by
$\mathcal{A}_{\text{hard}}=(Q_{1},\,2^{AP},\,\delta_{1},\,Q_{1,0},\,\mathcal{F}_{1})$
and
$\mathcal{A}_{\text{soft}}=(Q_{2},\,2^{AP},\,\delta_{2},\,Q_{2,0},\,\mathcal{F}_{2})$
as the NBAs associated with $\varphi_{\text{hard}}$ and
$\varphi_{\text{soft}}$, respectively, where the notations are defined
analogously as in Section II-B. Now we propose a way to compose $\mathcal{T}$,
$\mathcal{A}_{\text{hard}}$ and $\mathcal{A}_{\text{soft}}$ into a product
automaton.

###### Definition 1

The _parameterized_ product automaton
$\mathcal{A}_{p}\triangleq(Q_{p},\,\delta_{p},\,Q_{p,0},\,\mathcal{F}_{p})$ is
defined as: $Q_{p}=\Pi\times Q_{1}\times Q_{2}\times\\{1,2\\}$ are the states
with $q_{p}=\langle\pi,\,q_{1},\,q_{2},c\rangle\in Q_{p}$, $\forall\pi\in\Pi$,
$\forall q_{1}\in Q_{1}$, $\forall q_{2}\in Q_{2}$ and $\forall
c\in\\{1,2\\}$; $\delta_{p}:Q_{p}\times Q_{p}\rightarrow(\mathbb{R}_{\geq
0}\cup\\{\infty\\})^{3}$ maps each transition to a column vector such that
$\delta_{p}(\langle\pi,\,q_{1},\,q_{2},c\rangle,\langle\check{\pi},\,\check{q}_{1},\,\check{q}_{2},\check{c}\rangle)=[\alpha_{1},\alpha_{2},\alpha_{3}]^{\intercal}$,
where

  * •

$\alpha_{1}$ is the control cost for the robot to move from $\pi$ to
$\check{\pi}$, where $\alpha_{1}>0$ if $(\pi,\,\check{\pi})\in\rightarrow$,
otherwise $\alpha_{1}\triangleq\infty$;

  * •

$\alpha_{2}$ is the indicator for whether a transition violates the hard
constraints. It satisfies that $\alpha_{2}\triangleq 0$ if the following
conditions hold: (i) $L(\pi)\in\chi_{1}(q_{1},\,\check{q}_{1})$; (ii)
$\chi_{2}(q_{2},\,\check{q}_{2})\neq\emptyset$; (iii)
$q_{1}\notin\mathcal{F}_{1}$ and $\check{c}=c=1$; or
$q_{2}\notin\mathcal{F}_{2}$ and $\check{c}=c=2$; or
$q_{1}\in\mathcal{F}_{1}$, $c=1$ and $\check{c}=2$; or
$q_{2}\in\mathcal{F}_{2}$, $c=2$ and $\check{c}=1$. Otherwise,
$\alpha_{2}\triangleq\infty$.

  * •

$\alpha_{3}$ is the measure of how much a transition violates the soft
constraints, where
$\alpha_{3}\triangleq\texttt{Dist}(L(\pi),\,\chi_{2}(q_{2},\check{q}_{2}))$,
where the functions $\texttt{Dist}(\cdot)$ and $\chi_{2}(\cdot)$ for
$\mathcal{A}_{\text{soft}}$ are defined in Section II-B.

and $Q_{p,0}=\Pi_{0}\times Q_{1,0}\times Q_{2,0}\times\\{1\\}$,
$\mathcal{F}_{p}=\Pi\times\mathcal{F}_{1}\times Q_{2}\times\\{1\\}$ are the
sets of initial and accepting states, respectively. $\blacksquare$

An accepting run of $\mathcal{A}_{p}$ is an infinite run that starts from any
initial state and intersects with the accepting states infinitely often. Note
that the component $c$ above to ensure that an accepting run intersects with
the accepting states of both $\mathcal{A}_{\text{hard}}$ and
$\mathcal{A}_{\text{soft}}$ infinitely often. More details can be found in
Chapter 4 of [5]. Furthermore, since the workspace $\mathcal{T}$ is partially-
known, we denote by $\mathcal{T}^{t}$ the workspace model at time $t\geq 0$,
and the associated product automaton by $\mathcal{A}_{p}^{t}$.

To simplify the notation, given a finite run $R=q_{p}^{0}q_{p}^{1}\cdots
q_{p}^{S}$ of $\mathcal{A}_{p}$, where $q_{p}^{s}\in Q_{p}$, $\forall
s=0,1,\cdots,S$, we denote by
$\boldsymbol{\delta}(R)=\sum_{s=0}^{S-1}\delta_{p}(q_{p}^{s},q_{p}^{s+1}),$
where $\boldsymbol{\delta}(R)\in\mathbb{R}^{3}$ is the accumulated cost vector
$\delta_{p}$ along $R$. Similar definitions hold for
$\boldsymbol{\alpha}_{k}(R)\in\mathbb{R}$ as the accumulated $\alpha_{k}$ cost
along $R$, $\forall k=1,2,3$. We consider an accepting run of
$\mathcal{A}_{p}$ with the prefix-suffix structure: $R_{p}\triangleq
q_{p}^{1}q_{p}^{2}\cdots q_{p}^{S}\big{(}q_{p}^{S+1}q_{p}^{S+2}\cdots
q_{p}^{S+F}\big{)}^{\omega}$, where $q_{p}^{j}\in Q_{p}$, $\forall
j=1,2,\cdots,S+F$, where $S,F>0$. The plan prefix
$R_{p}^{\text{pre}}\triangleq q_{p}^{1}q_{p}^{2}\cdots q_{p}^{S}$ is executed
only once while the plan suffix $R_{p}^{\text{suf}}\triangleq
q_{p}^{S+1}q_{p}^{S+2}\cdots q_{p}^{S+F}$ is repeated infinitely often. Then
the total cost of ${R}_{p}$ is defined as:

$\texttt{C}_{\beta}(R_{p})\triangleq[1,\,\;\gamma]\otimes\begin{bmatrix}1\\\
1\\\
\;\beta\end{bmatrix}^{\intercal}\cdot\begin{bmatrix}\boldsymbol{\delta}(R_{p}^{\text{pre}})\\\
\boldsymbol{\delta}(R_{p}^{\text{suf}})\end{bmatrix},$ (3)

where $\texttt{C}_{\beta}(R_{p})\geq 0$; $\otimes$ is the Kronector product;
$\gamma\geq 0$ is a weighting parameter between the cost of the plan prefix
and suffix; $\beta\geq 0$ is a weighting parameter between total control cost
of the plan and the satisfaction of the soft task $\varphi_{\text{soft}}$.
Note that $\gamma$ is normally constant [11] (set to $1$ in this work), while
$\beta$ can change according to the robot’s internal model or the operator’s
_preference_. For instance, as the robot has more accurate workspace model,
$\beta$ can be increased to penalize the violation of $\varphi_{\text{soft}}$
such that $R_{p}$ satisfies $\varphi_{\text{soft}}$ more. Or the operator
prefers that $\beta$ is decreased so that $R_{p}$ satisfies
$\varphi_{\text{soft}}$ less and the robot reserves more power.

Given the initial values of $\gamma$ and $\beta$, an initial accepting run of
$\mathcal{A}_{p}$, denoted by $R_{p}^{0}$, can be found that minimizes the
total cost in (3). The algorithms are based on the nested Dijkstra’s search,
which are omitted here and details can be found in [11]. As a result, the
robot’s initial plan, denoted by $\tau_{r}^{0}$, can be derived by projecting
$R_{p}^{0}$ onto $\Pi$, as a sequence of regions that the robot should reach:
$\tau_{r}^{0}=\pi^{1}\pi^{2}\cdots\pi^{S}\big{(}\pi^{S+1}\pi^{S+2}\cdots\pi^{S+F}\big{)}^{\omega}$,
where $\pi^{j}$ is the projection of $q_{p}^{j}$ onto $\Pi$, $\forall
j=1,2,\cdots,S+F$.

###  IV-B Mixed-initiative Controller Design

After the system starts, the robot executes the initial plan $\tau_{r}^{0}$ by
reaching the sequence of regions defined by it. However, as described in
Section III-B, the robot controller is also influenced by the human input. In
the following, we show how to construct function $\kappa(\cdot)$ in (1) such
that the hard task $\varphi_{\text{hard}}$ is respected at all times for all
human inputs.

First, we need to find the set of product states $\mathcal{O}_{t}\subset
Q_{p}$ in $\mathcal{A}_{p}^{t}$ at time $t\geq 0$, such that once the robot
belongs to any state in $\mathcal{O}_{t}$ it means that
$\varphi_{\text{hard}}$ can not be satisfied any more.

###### Lemma 1

Assume that the robot belongs to state $q_{p}\in Q_{p}$ at time $t>0$. Then
the hard task $\varphi_{\text{hard}}$ can not be satisfied in the future, if
$\mathcal{A}_{p}^{t}$ remains unchanged and the minimal cost of all paths from
$q_{p}$ to any accepting state in $\mathcal{F}_{p}$ is $\infty$.

###### Proof:

Omitted as it is a simple inference of (3). ∎

Thus denote by $Q_{t}\subset Q_{p}$ the set of _reachable_ states by the robot
at time $t>0$. For each $q_{p}\in Q_{t}$, we perform a Dijkstra search to
compute the shortest distance from $q_{p}$ to all accepting states in
$\mathcal{F}_{p}$. Lastly, $\mathcal{O}_{t}$ is given as the subset of $Q_{t}$
that have an infinite cost to _all_ accepting states, i.e.,

$\mathcal{O}_{t}=\\{q_{p}\in
Q_{t}\,|\,\texttt{C}_{\beta}(\overline{R}_{q_{p},\,q_{F}})=\infty,\forall
q_{F}\in\mathcal{F}_{p}\\},$ (4)

where $\overline{R}_{q_{p},\,q_{F}}$ is the shortest path from $q_{p}$ to
$q_{F}$.

Given $\mathcal{O}_{t}$ above, we now design the function $\kappa(x,\Pi)$ in
(1) such that $\mathcal{O}_{t}$ can be avoided. Consider the function:

$\kappa(x,\Pi)=\kappa(x,\mathcal{O}_{t})\triangleq\frac{\rho(d_{t}-d_{s})}{\rho(d_{t}-d_{s})+\rho(\varepsilon+d_{s}-d_{t})}$
(5)

where
$d_{t}\triangleq\textbf{min}_{\langle\pi,q_{1},q_{2},c\rangle\in\mathcal{O}_{t}}\;\|x-\pi\|$
is the minimum distance between the robot and any region within
$\mathcal{O}_{t}$; $\rho(s)\triangleq e^{-1/s}$ for $s>0$ and
$\rho(s)\triangleq 0$ for $s\leq 0$, and $d_{s},\,\varepsilon>0$ are design
parameters as the safety distance and a small buffer. Thus the mixed-
initiative controller is given by

$u\triangleq u_{r}(x,\pi_{s},\pi_{g})+\kappa(x,\mathcal{O}_{t})\,u_{h}(t).$
(6)

As discussed in [4], the function $\kappa(\cdot)$ above is $0$ on the boundary
of undesired regions in $\mathcal{O}_{t}$, and close to $1$ when the robot is
away from $\mathcal{O}_{t}$ to allow for meaningful inputs from the operator.
This degree of closeness is tunable via changing $d_{t}$ and $d_{s}$. However,
the original definition in [4] only considers static obstacles, instead of the
general set $\mathcal{O}_{t}$.

###### Lemma 2

Assume that the navigation control $u_{r}$ is perpendicular to the boundary of
regions in $\mathcal{O}_{t}$ and points inwards. Then the robot can avoid
$\mathcal{O}_{t}$ under the mixed-initiative controller by (6) for all human
input $u_{h}$.

###### Proof:

The proof follows from Proposition 3 of [4]. Namely, for
$x\in\partial\mathcal{O}_{t}$ on the boundary of $\mathcal{O}_{t}$,

$\dot{x}^{\intercal}u_{r}=\|u_{r}(x)\|^{2}+\kappa(x)\,u_{h}(t)^{\intercal}u_{r}(x)>0$

since $\kappa(x)=0$ for $x\in\partial\mathcal{O}_{t}$. Thus the workspace
excluding all regions in $\mathcal{O}_{t}$ is positive invariant under
controller (6). In other words, if the navigation control avoids
$\mathcal{O}_{t}$, the same property is ensured by (6) for all human inputs. ∎

###  IV-C Discrete Plan Adaptation

In this section, we describe how the discrete plan $\tau_{r}^{t}$ at $t\geq 0$
can be updated to (i) accommodate changes in the partially-known workspace
model, and (ii) fulfill contingent tasks that are assigned by the operator
during run time.

####  IV-C1 Updated Workspace Model

The robot can explore new features of the workspace while executing the
discrete plan $\tau_{r}^{t}$ or being guided by the human operator. Thus the
motion model $\mathcal{T}^{t}$ can be updated as follows: (i) the transition
relation $\rightarrow$ is modified based on the status feedback from the
navigation controller, i.e., whether it is feasible to navigate from region
$\pi_{i}$ to $\pi_{j}$ without human inputs; (ii) the labeling function
$L(\pi)$ is changed based on the feedback from the robot’s sensing module,
i.e., the properties that region $\pi$ satisfies. For example, “the corridor
connecting two rooms is blocked” or “the object of interest is no longer in
one room”. Given the updated $\mathcal{T}^{t}$, the mapping function
$\delta_{p}$ of the product automaton $\mathcal{A}^{t}_{p}$ is re-evaluated.
Consequently, the current plan $\tau_{r}^{t}$ might not be optimal anymore
regarding the cost in (3). Thus we consider the following problem.

###### Problem 1

Update $\tau_{r}^{t}$ such that it has the minimum cost in (3) given the
updated model $\mathcal{T}^{t}$. $\blacksquare$

Given the set of reachable states $Q_{t}\subset Q_{p}$ at time $t>0$, for each
state $q_{p}\in Q_{t}$, we perform a nested Dijkstra’s search [11] to find the
accepting run that starts from $q_{p}$ and has the prefix-suffix structure
with the minimum cost defined in (3). Denote by $R^{t}_{+}(q_{p})$ this
optimal run for $q_{p}\in Q_{t}$. Moreover, let
$\mathbf{R}^{t}_{+}\triangleq\\{R^{t}_{+}(q_{p}),\,\forall q_{p}\in Q_{t}\\}$
collect all such runs. Then we find among $\mathbf{R}^{t}_{+}$ the accepting
run with the minimum cost, which becomes the updated run $R_{p}^{t+}$:

$R_{p}^{t+}\triangleq\textbf{argmin}_{R_{p}\in\mathbf{R}^{t}_{+}}\texttt{C}_{\beta}(R_{p}).$
(7)

Thus the updated plan $\tau_{r}^{t}$ is given by the projection of
$R_{p}^{t+}$ onto $\Pi$. Note that the above procedure is performed whenever
the motion model $\mathcal{T}^{t}$ is updated.

####  IV-C2 Contingent Task Fulfillment

As defined in (2), the operator can assign contingent and short-term tasks
$\varphi_{\text{temp}}$ to the robot during run time. Particularly, we
consider the following “ _pick-up and deliver_ ” task with deadlines, i.e.,

$(\varphi_{\textup{temp}}^{t},\,T_{sg})\triangleq(\Diamond(\pi_{s}\wedge\Diamond\pi_{g}),\,T_{sg}),$
(8)

where
$\varphi_{\textup{temp}}^{t}\triangleq\Diamond(\pi_{s}\wedge\Diamond\pi_{g})$
is the temporary task assigned at time $t>0$, meaning that the robot needs to
pick up some objects at region $\pi_{s}$ and deliver them to $\pi_{g}$ (note
that action propositions are omitted here, see [12]), where
$\pi_{s},\pi_{g}\in\Pi$; $T_{sg}>0$ is the _preferred_ deadline that task
$\varphi_{\textup{temp}}^{t}$ is accomplished. It can be verified that
$\varphi_{\textup{temp}}^{t}$ are sc-LTL formulas and can be fulfilled in
finite time. Assume that $\varphi_{\textup{temp}}^{t}$ is satisfied at time
$t^{\prime}>t$ then the delay is defined as $\overline{t}_{sg}\triangleq
t^{\prime}-T_{sg}$. We consider to following problem to incorporate
$\varphi_{\textup{temp}}^{t}$.

###### Problem 2

Update $R_{p}^{t}$ such that $\varphi_{\textup{temp}}^{t}$ can be fulfilled
_without_ delay (if possible) and with _minimum_ extra cost, while respecting
the hard constraints $\varphi_{\textup{hard}}$. $\blacksquare$

Assume the remaining optimal run at time $t>0$ is given by
$R_{p}^{t}=q^{k_{0}}_{p}q_{p}^{k_{0}+1}\cdots(q_{p}^{S}\cdots
q_{p}^{S+F})^{\omega}$ and $q^{k_{s}}_{p},q^{k_{g}}_{p}\in Q_{p}$ are the
$k_{s}$-th, $k_{g}$-th state of $R_{p}^{t}$, where $k_{s}\geq k_{g}\geq
k_{0}$. Since $R_{p}^{t}$ is optimal for the current $\mathcal{T}^{t}$, we
search for the index $k_{s}$ where the robot can deviate from $R_{p}^{t}$ to
reach region $\pi_{s}$ and back, and another index $k_{g}$ where the robot can
deviate from $R_{p}^{t}$ to reach $\pi_{g}$ and back. Denote by $R_{p}^{t+}$
the updated run after incorporating $\pi_{s},\pi_{g}$. In this way,
$\varphi_{\textup{temp}}^{t}$ is satisfied when $\pi_{g}$ is reached after
$\pi_{s}$ is reached, where
$t^{\prime}=\sum_{j=k_{0}}^{k_{g}}\alpha_{1}(q_{p}^{j},q_{p}^{j+1})$ is the
total time. Moreover, the total cost of $R^{t}_{p}$ in (3) is changed by
$\overline{\texttt{C}}_{\beta}(R^{t}_{p})\triangleq\texttt{C}_{\beta}(R^{t+}_{p})-\texttt{C}_{\beta}(R^{t}_{p})$.
Thus we formulate a 2-stage optimization below: first, we solve

$\displaystyle\overline{d}_{sg}$
$\displaystyle=\textbf{min}_{\\{k_{g}>k_{s}\geq
0\\}}\;\\{\overline{\texttt{C}}_{\beta}(R^{t}_{p})\\},\quad\textbf{s.t.}\quad\overline{t}_{sg}\leq
0,$ (9a) in order to find out whether it is possible to avoid delay while
satisfying $\varphi_{\textup{temp}}$. If no solution is found, we solve the
relaxed optimization that allows the deadline to be missed:
$\displaystyle\overline{d}_{sg}$
$\displaystyle=\textbf{min}_{\\{k_{g}>k_{s}\geq
0\\}}\;\\{\overline{t}^{\prime}_{sg}+\overline{\texttt{C}}_{\beta}(R^{t}_{p})\\},$
(9b)

where $\overline{d}_{sg}\geq 0$; $\overline{t}^{\prime}_{sg}=0$ if
$\overline{t}_{sg}\leq 0$ and $\overline{t}^{\prime}_{sg}=\overline{t}_{sg}$,
otherwise. Note that $\overline{\texttt{C}}_{\beta}(R^{t}_{p})$ is $\infty$ if
$\varphi_{\textup{hard}}$ is violated by $R_{p}^{t+}$.

Since the suffix of $R_{p}^{t}$ is repeated infinitely often, the choice of
indices $k_{s},k_{g}$ for (9) is finite. Thus (9) can be solved as follows:
starting from $k_{s}=k_{0}$, we iterate through
$k_{g}\in\\{k_{0}+1,\cdots,\,S+F\\}$ and compute the corresponding
$\overline{t}_{sg}$ and $\overline{d}_{sg}$ for both cases in (9). Then we
increase $k_{s}$ incrementally by $k_{s}=k_{0}+1$, iterate through
$k_{g}\in[k_{0}+2,\cdots,\,S+F]$ and compute $\overline{t}_{sg}$ and
$\overline{d}_{sg}$ for both cases. This procedure repeats itself _until_
$k_{s}=S+F-1$. Then, we find among these candidates if there is a pair
$k^{\star}_{s},k^{\star}_{g}$ that solves (9a). If so, they are the optimal
choice of $k_{s},k_{g}$. Otherwise, we search for the optimal solution to
(9b), of which the solution always exists as it is unconstrained. At last,
$R_{p}^{t+}$ is derived by inserting the product states associated with
$\pi_{s}$, $\pi_{g}$ at indices $k^{\star}_{s}$, $k^{\star}_{g}$ of
$R_{p}^{t}$, respectively.

###  IV-D Human Preference Learning

As discussed in Section IV-B, the mixed-initiative controller (6) allows the
operator to interfere the robot’s trajectory such that it deviates from its
discrete plan $\tau_{r}^{t}$, while always obeying $\varphi_{\textup{hard}}$.
This is beneficial as the robot could be guided to (i) explore _unknown_
features to update its workspace model, as described in Section IV-C1; and
(ii) follow the trajectory that is _preferred_ by the operator.

Particularly, as discussed in Section IV-A, the initial run $R_{p}^{0}$ is a
balanced plan between reducing the control cost and improving the satisfaction
of $\varphi_{\textup{soft}}$, where the weighting parameter is $\beta$ in (3).
Clearly, different choices of $\beta$ may result in different $R^{0}_{p}$. The
initial plan $R_{p}^{0}$ is synthesized under the initial value $\beta_{0}\geq
0$, which however might _not_ be what the operator prefers. In the following,
we present how the robot could learn about the preferred $\beta$ from the
operator’s inputs during run time.

Consider that at time $t\geq 0$, the robot’s past trajectory is given by
$\zeta|_{0}^{t}\triangleq\pi_{0}\pi_{1}\cdots\pi_{k_{t}}$. Assume now that
during time $[t,\,t^{\prime}]$, where $t^{\prime}>t>0$, via the mixed-
initiative controller in (6), the operator guides the robot to reach a
sequence of regions that s/he prefers, which is defined by:

$\zeta_{h}|_{t}^{t^{\prime}}\triangleq\pi^{\prime}_{1}\pi^{\prime}_{2}\cdots\pi^{\prime}_{H}$
(10)

where $\pi^{\prime}_{h}\in\Pi$, $\forall h=1,2\cdots H$ and $H\geq 1$ is the
length of $\zeta_{h}$ that can vary each time the operator acts. Afterwards,
the robot continues executing its current plan $\tau_{r}^{t}$. Thus, the
actual robot trajectory until time $t^{\prime}$ is given by
$\zeta_{h}|_{0}^{t^{\prime}}\triangleq\zeta|_{0}^{t}\,\zeta_{h}|_{t}^{t^{\prime}}$,
which is the concatenation of $\zeta|_{0}^{t}$ and
$\zeta_{h}|_{t}^{t^{\prime}}$.

###### Problem 3

Given the actual robot trajectory $\zeta_{h}|_{0}^{t^{\prime}}$, design an
algorithm to estimate the preferred value of $\beta$ as $\beta_{h}^{\star}$
such that $\zeta_{h}|_{0}^{t^{\prime}}$ corresponds to the optimal plan under
$\beta_{h}^{\star}$. $\blacksquare$

The above problem is closely related to the inverse reinforcement learning
(IRL) problem [24, 25], where the robot learns about the cost functions of the
system model based on demonstration of the preferred plans. On the other hand,
in reinforcement learning [26, 27] problem, the robot learns the optimal plan
given these functions.

As mentioned in [24], most problems of IRL are ill-posed. In our case, it
means that there are more than one $\beta_{h}^{\star}$ that render
$\zeta_{h}|_{0}^{t^{\prime}}$ to be the optimal plan under
$\beta_{h}^{\star}$. In order to improve _generalization_ such that the robot
could infer the human preference based on the human’s past inputs (instead of
simply repeating them), our solution is based on the maximum margin planning
algorithm from [25]. The general idea is to iteratively update $\beta$ via a
sub-gradient descent, where the gradient is computed based on the difference
in cost between $\zeta_{h}|_{0}^{t^{\prime}}$ and the optimal plan under the
current $\beta$.

First, we compute the set of all finite runs within
$\mathcal{A}_{p}^{t^{\prime}}$, denoted by $\mathbf{R}_{h}^{t^{\prime}}$, that
are associated with $\zeta_{h}|_{0}^{t^{\prime}}$. It can be derived
iteratively via a breadth-first graph search [5]. Among
$\mathbf{R}_{h}^{t^{\prime}}$, we find the one with the minimal cost over
$\alpha_{3}$, i.e.,

$R^{\star}_{h}\triangleq\textbf{argmin}_{R\in\mathbf{R}_{h}^{t^{\prime}}}\;\boldsymbol{\alpha}_{3}(R).$
(11)

Let $R^{\star}_{h}\triangleq q_{1}q_{2}\cdots q_{H}$, where $q_{h}\in Q_{p}$,
$\forall h=1,2,\cdots,H$. Denote by $\beta_{k}$ the value of $\beta$ at the
$k$-th iteration, for $k\geq 0$. Note that $\beta_{0}\triangleq\beta_{t}$,
where $\beta_{t}$ is the value of $\beta$ at time $t>0$. For the $k$-th
iteration, we find the optimal run from $q_{1}$ to $q_{H}$ under $\beta_{k}$
with certain margins, i.e.,

$\hat{R}_{\beta_{k}}^{\star}\triangleq\textbf{argmin}_{R\in\mathbf{R}_{q_{1}q_{H}}}\;\Big{(}\texttt{C}_{\beta_{k}}(R)-M(R,R^{\star}_{h})\Big{)}$
(12)

where $\mathbf{R}_{q_{1}q_{H}}$ is the set of _all_ runs from $q_{1}$ to
$q_{H}$ in $\mathcal{A}_{p}^{t^{\prime}}$; and $M:Q^{H}\times
Q^{H}\rightarrow\mathbb{N}$ is the margin function [25]:

$M(R,\,R^{\star}_{h})=|\\{(q_{s},q_{t})\in R\,|\,(q_{s},q_{t})\notin
R_{h}^{\star}\\}|,$ (13)

which returns the number of edges within $R$ that however do not belong to
$R_{h}^{\star}$. The margin function decreases the total cost
$\texttt{C}_{\beta_{k}}(R)$ by the difference between $R$ and $R_{h}^{\star}$.
It can improve generalization and help address the ill-posed nature of Problem
3. To solve (12), we first modify $\mathcal{A}_{p}^{t^{\prime}}$ by reducing
the $\alpha_{1}$ cost of each edge $(q_{s},q_{t})\in R_{h}^{\star}$ by one.
Then a Dijkstra shortest path search can be performed over the modified
$\mathcal{A}_{p}$ to find the shortest run from $q_{1}$ to $q_{H}$ that
minimizes the cost with margins in (12). Given $\hat{R}_{\beta_{k}}^{\star}$,
we can compute the sub-gradient [28] that $\beta_{k}$ should follow:

$\nabla\beta_{k}=\lambda\cdot\beta_{k}+\big{(}\boldsymbol{\alpha}_{3}(R^{\star}_{h})-\boldsymbol{\alpha}_{3}(\hat{R}_{\beta_{k}}^{\star})\big{)},$
(14)

where $\nabla\beta_{k}\in\mathbb{R}$ and $\lambda>0$ is a design parameter.
Thus, at this iteration the value of $\beta_{k}$ is updated by

$\beta_{k+1}=\beta_{k}-\theta_{k}\cdot\nabla\beta_{k},$ (15)

where $\theta_{k}>0$ is the step size or learning rate [26]. Given the updated
$\beta_{k+1}$, the same process in (11)-(15) is repeated until the difference
$|\beta_{k+1}-\beta_{k}|$ is less than a predefined threshold $\varepsilon>0$.
At last, the value of $\beta_{t}$ is updated to $\beta_{k+1}$. The discussion
above is summarized in Alg. 1. Each time the human operator guides the robot
to reach a new sequence of regions, the estimation of the value of $\beta_{t}$
is updated by running Alg. 1. In the following, we show that Alg. 1 ensures
the convergence of $\\{\beta_{k}\\}$.

Input: $\mathcal{A}^{t}_{p}$, $\zeta_{h}|_{0}^{t^{\prime}}$, $\beta_{t}$,
$\varepsilon$

1 Initialize $\beta_{k}=\beta_{t}$ for iteration $k=0$;

2 while  _$|\beta_{k+1}-\beta_{k}| >\varepsilon$_ do // Iteration $k$

3 Compute $R_{h}^{\star}$ in (11) given $\zeta_{h}|_{0}^{t^{\prime}}$;

4 Find $\hat{R}_{\beta_{k}}^{\star}$ in (12) given $\beta_{k}$ and
$R_{h}^{\star}$;

5 Compute $\nabla\beta_{k}$ by (14) and update $\beta_{k}$ by (15);

return $\beta_{t}^{+}=\beta_{k+1}$

Algorithm 1 On-line IRL algorithm for $\beta$.

###### Lemma 3

The sequence $\\{\beta_{k}\\}$ in Alg. 1 converges to a fixed
$\beta_{l}^{\star}\geq 0$ and the optimal plan under $\beta_{l}^{\star}$ is
$\zeta_{h}|_{0}^{t^{\prime}}$.

###### Proof:

Firstly, the optimal run $R_{h}^{\star}$ associated with
$\zeta_{h}|_{0}^{t^{\prime}}$ under $\beta_{h}^{\star}$ minimizes the balanced
cost $\texttt{C}_{\beta}$ from (3), i.e.,

$\texttt{C}_{\beta}(R_{h}^{\star})\leq\texttt{C}_{\beta}(R),\quad\forall
R\in\mathbf{R}_{q_{1}q_{H}},$ (16)

where $\mathbf{R}_{q_{1}q_{H}}$ is defined in (12). Solving (16) directly can
be computationally expensive due to the large set $\mathbf{R}_{q_{1}q_{H}}$.
We introduce a slack variable $\xi\in\mathbb{R}$ to relax the constraints:

$\begin{split}&\textbf{min}_{\beta\geq
0}\quad\frac{\lambda}{2}\beta^{2}+\xi\\\
&\textbf{s.t.}\quad\texttt{C}_{\beta}(R^{\star}_{h})-\xi\leq\min_{R\in\mathbf{R}_{q_{1}q_{H}}}\Big{(}\texttt{C}_{\beta}(R)-M(R,R^{\star}_{h})\Big{)},\\\
\end{split}$ (17)

where $\lambda>0$ is the same as in (14) and the margin function $M(\cdot)$ is
from (13). Thus, by enforcing the slack variables to be tight, $\beta$ also
minimizes the combined cost function:

$\frac{\lambda}{2}\beta^{2}+\texttt{C}_{\beta}(R^{\star}_{h})-\min_{R\in\mathbf{R}_{q_{1}q_{H}}}\Big{(}\texttt{C}_{\beta}(R)-M(R,R^{\star}_{h})\Big{)},$
(18)

which is convex but non-differentiable. Instead, we compute the sub-gradient
[28] of (18):
$\nabla\beta=\lambda\beta+\big{(}\boldsymbol{\alpha}_{3}(R^{\star}_{h})-\boldsymbol{\alpha}_{3}(\hat{R}_{\beta}^{\star})\big{)}$
and
$\hat{R}_{\beta}^{\star}=\textbf{argmin}_{R\in\mathbf{R}_{q_{1}q_{H}}}{(}\texttt{C}_{\beta}(R)-M(R,R^{\star}_{h}){)}$,
which is equivalent to (14).

Lastly, by the strong convexity of (18) and Theorem 1 of [25], the estimation
$\beta_{t}$ approaches the optimal $\beta_{l}^{\star}$ with linear convergence
rate under constant stepsize $\theta_{k}=\theta$, i.e.,
$|\beta_{k}-\beta_{l}^{\star}|^{2}\leq(1-\theta\lambda)^{k+1}|\beta_{0}-\beta_{l}^{\star}|^{2}+\frac{\theta|\nabla\beta|_{\max}}{\lambda}$.
A detailed analysis on this deviation can be found in [25, 28]. ∎

###### Remark 2

It is worth noting that the convergent value $\beta_{l}^{\star}$ might be
_different_ from the preferred $\beta^{\star}_{h}$, while they both satisfy
(16) with the same optimal run $R_{h}^{\star}$. However, the margin function
in (17) ensures that $\beta_{l}^{\star}$ is _better_ than or at least
equivalent to $\beta_{h}^{\star}$ in terms of the similarity between
$R_{h}^{\star}$ and the run with the second minimum cost by (3).

Method | $|\mathcal{A}_{p}|$ | $\beta_{l}^{\star}$ | NO. of Dijkstra | Time[$\mathrm{s}$]
---|---|---|---|---
Alg.1 | 25 | 13.4 | 8 | 3.8
M1 | 25 | 10.0 | 200 | 124.4
M2 | 25 | 11.7 | 350 | 337.2
Alg.1 | 100 | 16.5 | 12 | 150.8
M1 | 100 | 14.2 | 200 | 2203.5
M2 | 100 | – | 800+ | 3000+

Table I: Comparison of computational complexity and performance of Alg. 1 and
two alternative methods in Example 1.

Now we show the computational efficiency of Alg. 1 compared with two straight-
forward solutions: (M1) choose the optimal $\beta$ among a set of guessed
values of $\beta$, denoted by $S_{\beta}$; (M2) solve (16) directly by
enumerating all runs in $\mathbf{R}_{q_{1}q_{H}}$. The first method’s accuracy
relies on $S_{\beta}$ being large, which however results in high computational
cost. Similarly, the second method relies on evaluating _every_ run in
$\mathbf{R}_{q_{1}q_{H}}$, the size of which is combinatorial to the size of
$\mathcal{A}^{t}_{p}$. The following example shows some numerical comparison.

###### Example 1

Assume that $\beta_{h}^{\star}=15$ and initially $\beta_{0}=0$. We use three
methods: Alg. 1, M1 and M2 above to estimate $\beta_{h}^{\star}$. As shown in
Table I, we compare the final convergence $\beta_{l}^{\star}$ and the
computation time under varying sizes of $\mathcal{A}_{p}$. It can be seen that
the computation time for Alg. 1 is significantly less than M1 and M2, where
for the second case M2 fails to converge within $50$\mathrm{min}$$.
$\blacksquare$

##  V The Integrated System

In this section, we describe the the real-time execution of the integrated
system given the components in Section IV. Then we discuss the computational
complexity.

###  V-A Human-in-the-loop Motion and Task Planning

The complete algorithm is shown in Alg. 2. Before the system starts, given the
initial model $\mathcal{T}^{0}$ and the task formulas in (2), the initial plan
$\tau_{r}^{0}$ is synthesized by the algorithm from Section IV-A under the
initial $\beta_{0}$. From $t=0$, the robot executes $\tau_{r}^{0}$ by
following the sequence of goal regions, see Lines 1-3. Meanwhile, the operator
can directly modify the control input $u(\cdot)$ via (6) to change the robot’s
trajectory. Thus, the robot can explore regions that are not in its initial
plan and update the model $\mathcal{T}^{t}$ as described in Section IV-C1. As
a result, the plan $\tau_{r}^{t}$ is updated by (7) accordingly, see Lines
4-5. Moreover, as described in Section IV-C2, the operator can assign
temporary tasks with deadlines as in (8), for which $\tau_{r}^{t}$ is modified
by solving (9), see Lines 6-7. Last but not least, each time the operator
guides the robot to follow a new trajectory, the parameter $\beta$ is updated
via Alg. 1 to estimate the human preference. Then, its current plan
$\tau_{r}^{t}$ is updated using the updated $\beta$, see Lines 8-12. The above
procedure repeats until the system is terminated.

###### Theorem 4

Alg. 2 above fulfills the three control objectives of Section III-D, i.e., (I)
$\varphi_{\textup{hard}}$ is satisfied for all time; (II) each
$\varphi_{\textup{temp}}$ is satisfied in finite time; and (III) the
satisfaction of $\varphi_{\textup{soft}}$ adapts to the human inputs.

###### Proof:

(I) Firstly, both the initial synthesis algorithm in Section IV-A and the plan
adaptation algorithm in Section IV-C1 ensure $\varphi_{\textup{hard}}$ is
satisfied by minimizing the total cost in (3). Then Lemma 2 ensures that
$\varphi_{\textup{hard}}$ is respected for all possible inputs from the human
operator. (II) The combined cost (9) ensures that $\varphi_{\textup{temp}}$ is
satisfied within finite time. (III) Convergence of the learning Alg. 1 is
shown in Lemma (3). Thus, the updated plan $\tau_{r}^{t}$ under the learned
value of $\beta$ adapts to the plan preferred by the operator. ∎

Input: $\mathcal{T}^{t}$, $\varphi_{\text{hard}}$, $\varphi_{\text{soft}}$,
$\beta_{0}$, $u_{h}(t)$, $(\varphi^{t}_{\text{temp}},T_{sg})$

1 Compute $\mathcal{A}_{p}^{0}$ and construct initial plan $\tau_{r}^{0}$
under $\beta_{0}$;

2 forall  _$t\geq 0$_ do

3 Compute $u_{r}(\cdot)$ in (6) to reach next $\pi^{j}\in\tau_{r}^{t}$;

4 if  _$\mathcal{T}^{t}$ updated_ then // Model update

5 Update product $\mathcal{A}_{p}^{t}$ and plan $\tau_{r}^{t}$ by (7);

6 if  _$(\varphi_{\textup{temp}}^{t},T_{sg})$ received_ then // Temp. task

7 Update plan $\tau_{r}^{t}$ by solving (9);

8 if  _$\|u_{h}(t)\| >0$_ then // Human input

9 Compute control $u(t)$ by (6);

10 Compute $\zeta_{h}|_{0}^{t}$ by (10);

11 Learn $\beta_{l}^{\star}$ by Alg. 1 and set
$\beta_{t}^{+}=\beta_{l}^{\star}$;

12 Update $\tau_{r}^{t}$ by (7) given the learned $\beta_{t}^{+}$;

13 return $u(t)$, $\tau_{r}^{t}$, $\beta_{t}^{+}$

Algorithm 2 Mixed-initiative Motion and Task Planning

###  V-B Computational Complexity

The process to synthesize $\tau_{r}^{t}$ given $\mathcal{A}^{t}_{p}$ via Alg.
2 (in Line 1) and the plan revision given $\mathcal{T}^{t}$ (in Line 5) both
have complexity $\mathcal{O}(|\mathcal{A}^{t}_{p}|^{2})$ [11]. The adaptation
algorithm for temporary tasks (in Line 7) has complexity
$\mathcal{O}(|R_{p}^{t}|^{2})$. Lastly, the learning Alg. 1 (in Line 11) has
complexity $\mathcal{O}(|R_{h}^{\star}|^{2})$, where $|R_{h}^{\star}|$ is the
length of the optimal run from (11).

##  VI Case Study

In this section, we present numerical studies both in simulation and
experiment. The Robot Operation System (ROS) is used as the simulation and
experiment platform. All algorithms are implemented in Python 2.7 and
available online [29]. All computations are carried out on a laptop (3.06GHz
Duo CPU and 8GB of RAM).

###  VI-A Simulation

Figure 1: Office environment in Gazebo with TIAGo robot, where the regions of
interest and allowed transitions are marked.

####  VI-A1 Workspace and Robot Description

Consider the simulated office environment in Gazebo as shown in Fig. 1 with
dimension $100m\times 45m$, in which there are 9 regions of interest (denoted
by $r_{0},\cdots,r_{8}$) and 4 corridors (denoted by $c_{1}$, $\cdots$,
$c_{4}$). The transition relations are determined by whether there exists a
collision-free path from the center of one region to another, without crossing
other regions.

We simulate the TIAGo robot from PAL robotics, of which the navigation control
$u_{r}(\cdot)$ with obstacle avoidance, localization and mapping are all based
on the ROS navigation stack. The human operator monitors the robot motion
through Rviz. Moreover, the control $u_{h}(\cdot)$ from the operator can be
generated from a keyboard or joystick, while the temporary task in LTL
formulas $\varphi_{\text{temp}}$ are specified via ROS messages. More details
can be found in the software implementation [29] and simulation video [30].

Figure 2: The robot’s trajectory in simulation case One, where the robot’s
initial plan is in blue, the human-guided part is in red, and the updated plan
is in green.

####  VI-A2 Case One

The hard task for _delivery_ is given by
$\varphi_{1,\text{hard}}=\big{(}\square\Diamond(r_{0}\wedge\Diamond(r_{7}\wedge\Diamond
r_{8}))\big{)}\wedge\big{(}\square\Diamond(r_{2}\wedge\Diamond(r_{3}\vee
r_{6}))\big{)}\wedge\big{(}\square\neg r_{5}\big{)}$, i.e., to transfer
objects from $r_{0}$ to $r_{7}$ (then $r_{8}$) and from $r_{2}$ to $r_{3}$ (or
$r_{6}$), while avoiding $r_{5}$ for all time. The soft task is
$\varphi_{1,\text{soft}}=(\square\neg c_{4})$, i.e., to avoid $c_{4}$ if
possible. It took $0.2s$ to compute the parameterized product automaton, which
has $312$ states and $1716$ transitions. The parameter $\beta$ is initially
set to a large value $30$, thus the initial plan satisfies both the soft and
hard tasks but with a large cost due to the long traveling distance, as shown
in Fig. 2. During $[700s,950s]$, the operator drives the robot to go through
corridor $c_{4}$ and reach $r_{8}$, which violates the soft task
$\varphi_{1,\text{soft}}$. As a result, $\beta$ is updated by Alg. 1 and the
final value is $16.35$ after $20$ iterations with $\varepsilon=0.2$, as shown
in Fig. 5. Namely, the robot has learned that the operator allows _more
violation_ of the soft task to reduce the total cost. The resulting updated
plan is shown in Fig. 2. Moreover, to demonstrate the ensured safety in Lemma
2, the human operator drives the robot towards $r_{5}$ during $[1250s,1350s]$,
which is not allowed by $\varphi_{\textup{hard}}$. The weighting function
$\kappa(\cdot)$ in the mixed controller (6) approaches $0$. Thus the robot
still follows its updated plan and avoids $r_{5}$. The mixed control inputs
during these periods are shown in Fig. 3.

Figure 3: The mixed control for linear velocity during two time periods when
the human control (in red) is active.

Figure 4: The robot’s trajectory in simulation case Two. The robot’s
trajectory while performing the temporary task is in magenta.

Figure 5: Evolution of the learned value of $\beta$ in simulation case One
(left) and case Two (right).

####  VI-A3 Case Two

The hard task for _surveillance_ is given by
$\varphi_{2,\text{hard}}=(\square\Diamond r_{2})\wedge(\square\Diamond
r_{3})\wedge(\square\Diamond r_{8})$, i.e., to surveil regions $r_{2}$,
$r_{3}$ and $r_{8}$ infinitely often. The soft task for extra performance is
$\varphi_{2,\text{soft}}=\square\Diamond\big{(}r_{4}\rightarrow(\neg
r_{5}\mathsf{U}\Diamond r_{6})\big{)}$, i.e., to collect goods from region
$r_{4}$ and drop it at $r_{6}$ (without crossing $r_{5}$ before that).
Moreover, the workspace model in this case is _different_ from the initial
model that the corridor $c_{2}$ has been blocked. By following Alg. 2, it took
$0.17s$ to compute the product automaton, which has $418$ states and $3360$
transitions. Initially, $\beta=0$ meaning that the initial plan $\tau_{r}^{0}$
only satisfies $\varphi_{2,\text{hard}}$ while $\varphi_{2,\text{soft}}$ is
fully relaxed. During $[150s,250s]$, the operator drives the robot to sense
that the corridor $c_{2}$ has been blocked. As a result, the discrete plan
$\tau_{r}^{0}$ is updated such that the robot chooses to reach $r_{8}$ from
$r_{2}$ via $c_{1}$, as shown in Fig. 4. Afterwards, during $[1100s,1200s]$,
the operator drives the robot to $r_{4}$ after reaching $r_{2}$, which
satisfies part of $\varphi_{2,\text{soft}}$. As a result, $\beta$ is increased
by Alg. 1 to $11.3$ after $12$ iterations with $\varepsilon=0.1$, as shown in
Fig. 5. Namely, the robot has learned that the soft task should be satisfied
_more_. Lastly, at time $2100s$, the operator assigns a temporary task
$\varphi_{\text{temp}}=\Diamond(r_{1}\wedge\Diamond r_{7})$ with a deadline
$2700s$, i.e., to deliver an object from $r_{1}$ to $r_{7}$. This temporary
task is incorporated into $\tau_{r}^{t}$ and is fulfilled at $2400s$, which is
shown in Fig. 4.

###  VI-B Experiment

The experiment setup involves a TurtleBot within the office environment at
Automatic Control Lab, KTH. Details are omitted here due to limited space,
which are given in the the software implementation [29] and experiment video
[30].

Figure 6: The human-in-the-loop experiment setup, where the robot is
controlled by both its autonomous controller and the human inputs.

####  VI-B1 Workspace and Task Specification

The office environment consists of three office rooms ($r_{1}$, $r_{2}$,
$r_{3}$) and one corridor $r_{0}$, as shown in Fig. 6. The robot’s task
specification is similar to case study Two above, i.e., the hard task is given
by $\varphi_{\text{hard}}=\square\Diamond r_{0}\wedge\square\Diamond r_{1}$
(to surveil regions $r_{0}$ and $r_{1}$) while the soft task is
$\varphi_{\text{soft}}=\square\Diamond r_{2}\wedge\square\Diamond r_{3}$ (to
surveil regions $r_{2}$ and $r_{3}$). The TurtleBot is controlled via ROS
navigation stack and behaves similarly to the TIAGo robot in Section VI-A.

####  VI-B2 Experiment Results

Since $\beta$ is initially set to $0$, the robot only surveils $r_{0}$ and
$r_{1}$ for the hard task, as shown in Fig. 7. From $t=59s$, the operator
starts driving the robot towards $r_{2}$ and back to $r_{0}$ until $t=137s$.
As a result, the estimated $\beta_{t}$ is updated by Alg. 1 given the robot’s
past trajectory. The final convergence value is $1.58$ with $\varepsilon=0.01$
after 15 iterations. Then updated plan is shown in Fig. 7 which intersects
with not only regions $r_{0}$ and $r_{1}$ for the hard task, but also regions
$r_{2}$ and $r_{3}$ for the soft task. Notice that the operator only needs to
interfere the robot’s motion for a small fraction of the operation time.

Figure 7: The robot’s trajectory in the experiment study, where the robot’s
initial plan is in blue (left), the human-guided segment is in red (middle),
and the updated plan is in green (right).

##  VII Summary and Future Work

In this paper, we present a human-in-the-loop task and motion planning
strategy for mobile robots with mixed-initiative control. The proposed
coordination scheme ensures the satisfaction of high-level LTL tasks given the
human initiative both through continuous control inputs and discrete task
assignments. Future work includes consideration of multi-robot systems.

## References

  * [1] M. Dunbabin and L. Marques, “Robots for environmental monitoring: Significant advancements and applications,” _Robotics & Automation Magazine, IEEE_, vol. 19, no. 1, pp. 24–39, 2012.
  * [2] T. Fong, I. Nourbakhsh, and K. Dautenhahn, “A survey of socially interactive robots,” _Robotics and autonomous systems_ , vol. 42, no. 3, pp. 143–166, 2003.
  * [3] M. A. Goodrich and A. C. Schultz, “Human-robot interaction: a survey,” _Foundations and trends in human-computer interaction_ , vol. 1, no. 3, pp. 203–275, 2007.
  * [4] S. G. Loizou and V. Kumar, “Mixed initiative control of autonomous vehicles,” in _Robotics and Automation (ICRA), IEEE International Conference on_ , 2007, pp. 1431–1436.
  * [5] C. Baier and J.-P. Katoen, _Principles of model checking_. MIT press Cambridge, 2008.
  * [6] G. E. Fainekos, A. Girard, H. Kress-Gazit, and G. J. Pappas, “Temporal logic motion planning for dynamic robots,” _Automatica_ , vol. 45, no. 2, pp. 343–352, 2009.
  * [7] A. Bhatia, L. E. Kavraki, and M. Y. Vardi, “Sampling-based motion planning with temporal goals,” in _Robotics and Automation (ICRA), IEEE International Conference on_ , 2010, pp. 2689–2696.
  * [8] X. C. Ding, S. L. Smith, C. Belta, and D. Rus, “Mdp optimal control under temporal logic constraints,” in _Decision and Control (CDC), IEEE Conference on_ , 2011, pp. 532–538.
  * [9] C. Belta, B. Yordanov, and E. A. Gol, _Formal Methods for Discrete-Time Dynamical Systems_. Springer, 2017, vol. 89.
  * [10] A. Ulusoy, S. L. Smith, X. C. Ding, C. Belta, and D. Rus, “Optimality and robustness in multi-robot path planning with temporal logic constraints,” _The International Journal of Robotics Research_ , vol. 32, no. 8, pp. 889–911, 2013.
  * [11] M. Guo and D. V. Dimarogonas, “Multi-agent plan reconfiguration under local LTL specifications,” _The International Journal of Robotics Research_ , vol. 34, no. 2, pp. 218–235, 2015.
  * [12] ——, “Task and motion coordination for heterogeneous multiagent systems with loosely coupled local tasks,” _IEEE Transactions on Automation Science and Engineering_ , vol. 14, no. 2, pp. 797–808, 2017.
  * [13] J. Tumova and D. V. Dimarogonas, “Multi-agent planning under local LTL specifications and event-based synchronization,” _Automatica_ , vol. 70, pp. 239–248, 2016.
  * [14] H. Kress-Gazit, G. E. Fainekos, and G. J. Pappas, “Temporal-logic-based reactive mission and motion planning,” _Robotics, IEEE Transactions on_ , vol. 25, no. 6, pp. 1370–1381, 2009.
  * [15] L. Feng, C. Wiltsche, L. Humphrey, and U. Topcu, “Synthesis of human-in-the-loop control protocols for autonomous systems,” _IEEE Transactions on Automation Science and Engineering_ , vol. 13, no. 2, pp. 450–462, 2016.
  * [16] J. Fu and U. Topcu, “Pareto efficiency in synthesizing shared autonomy policies with temporal logic constraints,” in _Robotics and Automation (ICRA), IEEE International Conference on_ , 2015, pp. 361–368.
  * [17] N. Jansen, M. Cubuktepe, , and U. Topcu, “Synthesis of shared control protocols with provable safety and performance guarantees,” in _In Proc. of the American Control Conference (ACC)_ , 2017.
  * [18] D. E. Koditschek and E. Rimon, “Robot navigation functions on manifolds with boundary,” _Advances in Applied Mathematics_ , vol. 11, no. 4, pp. 412–442, 1990.
  * [19] D. Panagou, D. M. Stipanovic, and P. G. Voulgaris, “Multi-objective control for multi-agent systems using lyapunov-like barrier functions,” in _Decision and Control (CDC), IEEE Conference on_. IEEE, 2013, pp. 1478–1483.
  * [20] L. Wang, A. D. Ames, and M. Egerstedt, “Multi-objective compositions for collision-free connectivity maintenance in teams of mobile robots,” in _Decision and Control (CDC), IEEE Conference on_. IEEE, 2016, pp. 2659–2664.
  * [21] O. Kupferman and M. Y. Vardi, “Model checking of safety properties,” _Formal Methods in System Design_ , vol. 19, no. 3, pp. 291–314, 2001.
  * [22] P. Gastin and D. Oddoux, “Fast LTL to büchi automata translation,” in _Computer Aided Verification_. Springer, 2001, pp. 53–65.
  * [23] S. M. LaValle, _Planning algorithms_. Cambridge university press, 2006.
  * [24] A. Y. Ng, S. J. Russell _et al._ , “Algorithms for inverse reinforcement learning.” in _Icml_ , 2000, pp. 663–670.
  * [25] N. D. Ratliff, J. A. Bagnell, and M. A. Zinkevich, “Maximum margin planning,” in _Proceedings of the 23rd international conference on Machine learning_. ACM, 2006, pp. 729–736.
  * [26] R. S. Sutton and A. G. Barto, _Reinforcement learning: An introduction_. MIT press Cambridge, 1998, vol. 1, no. 1.
  * [27] D. P. Bertsekas and J. N. Tsitsiklis, _Neuro-Dynamic Programming_ , 1st ed. Athena Scientific, 1996.
  * [28] N. Z. Shor, _Minimization methods for non-differentiable functions_. Springer Science & Business Media, 2012, vol. 3.
  * [29] Software, http://github.com/MengGuo/mix_initiative.
  * [30] Video, https://vimeo.com/230487800,232727691.


distribution function, it is convenient to rewrite Eq. (A) into the following
form

$\displaystyle f=f_{\mathrm{T=0}}+f_{\mathrm{T\neq
0}}+\tilde{f}_{\mathrm{T\neq 0}}$ (A.28)

where the function are defined as

$\displaystyle
f_{\mathrm{T=0}}=\Theta\left(\frac{\widetilde{\mu}-E}{T}\right),$ (A.29)
$\displaystyle f_{\mathrm{T\neq
0}}=\frac{1}{2}e^{-|E-\widetilde{\mu}|/T}\mathrm{sgn}\left(\frac{E-\widetilde{\mu}}{T}\right),\quad\tilde{f}_{\mathrm{T\neq
0}}=\frac{1}{2}e^{-|E-\widetilde{\mu}|/T}\tanh\left(\frac{E-\widetilde{\mu}}{2T}\right)$
(A.30)

In Fig. (A.2) we plot the zero (purple lines, $f_{T=0}$) and finite
temperature components of the Fermi distribution as a function of energy
choosing in this example the chemical potential $\widetilde{\mu}=0.461$ MeV at
temperature $T=0.02$ MeV and at $T=0.2$ MeV.

Figure A.2: The zero and finite temperature components of the decomposition
here considered for Fermi distribution as a function of energy with chemical
potential $\widetilde{\mu}=0.461$ MeV at temperature $T=0.02$ MeV and $T=0.2$
MeV. The purple line represents the zero temperature component
$f_{\mathrm{T}=0}$, and blue and red lines represent the finite temperature
components $f_{\mathrm{T\neq 0}}$ and $\tilde{f}_{\mathrm{T\neq 0}}$
respectively.

This Fig. (A.2) shows that the two finite temperature components
$f_{\mathrm{T\neq 0}}$ and $\tilde{f}_{\mathrm{T\neq 0}}$ are naturally
centered around the Fermi surface of the system $E=\tilde{\mu}$. Both finite
$T$ always have the same sign, so there is no opportunity generate numerical
noise; this can be the case when considering brute force removal of $T=0$
limit from the low temperature distribution. The dominant contribution (blue
lines, $f_{T\neq 0}$) are containing the decaying exponential function of the
Fermi distribution at large argument. The remainder (red lines,
$\tilde{f}_{T\neq 0}$) has the shape of a derivative of a finite support
representation of delta-function.

## References

  * Aaij et al. (2020a) Aaij, R. et al. (2020a).  Measurement of CP observables in B± $\rightarrow$ DK± and B± $\rightarrow$ D$\pi$± with D $\rightarrow$ ${K}_{\mathrm{S}}^{0}{K}^{\pm}{\pi}^{\mp}$ decays.  _JHEP_ , 06, p. 058.  10.1007/JHEP06(2020)058.
  * Aaij et al. (2020b) Aaij, R. et al. (2020b).  Observation of Several Sources of $CP$ Violation in $B^{+}\to\pi^{+}\pi^{+}\pi^{-}$ Decays.  _Phys. Rev. Lett._ , 124(3), p. 031801.  10.1103/PhysRevLett.124.031801.
  * Aarts et al. (2011) Aarts, G., C. Allton, S. Kim, M. P. Lombardo, M. B. Oktay, S. M. Ryan, D. K. Sinclair, and J. I. Skullerud (2011).  What happens to the $\Upsilon$ and $\eta_{b}$ in the quark-gluon plasma? Bottomonium spectral functions from lattice QCD.  _JHEP_ , 11, p. 103.  10.1007/JHEP11(2011)103.
  * Ade et al. (2014) Ade, P. A. R. et al. (2014).  Planck 2013 results. XVI. Cosmological parameters.  _Astron. Astrophys._ , 571, p. A16.  10.1051/0004-6361/201321591.
  * Ade et al. (2016) Ade, P. A. R. et al. (2016).  Planck 2015 results. XIII. Cosmological parameters.  _Astron. Astrophys._ , 594, p. A13.  10.1051/0004-6361/201525830.
  * Aghanim et al. (2020a) Aghanim, N. et al. (2020a).  Planck 2018 results. VI. Cosmological parameters.  _Astron. Astrophys._ , 641, p. A6.  10.1051/0004-6361/201833910.  [Erratum: Astron.Astrophys. 652, C4 (2021)].
  * Aghanim et al. (2020b) Aghanim, N. et al. (2020b).  Planck 2018 results. VI. Cosmological parameters.  _Astron. Astrophys._ , 641, p. A6.  10.1051/0004-6361/201833910.  [Erratum: Astron.Astrophys. 652, C4 (2021)].
  * Amhis et al. (2021) Amhis, Y. S. et al. (2021).  Averages of b-hadron, c-hadron, and $\tau$-lepton properties as of 2018\.  _Eur. Phys. J. C_ , 81(3), p. 226.  10.1140/epjc/s10052-020-8156-7.
  * Anderson and Witting (1974) Anderson, J. and H. Witting (1974).  A relativistic relaxation-time model for the Boltzmann equation.  _Physica_ , 74(3), pp. 466–488.  ISSN 0031-8914.  https://doi.org/10.1016/0031-8914(74)90355-3.
  * Arnold and McLerran (1987) Arnold, P. B. and L. D. McLerran (1987).  Sphalerons, Small Fluctuations and Baryon Number Violation in Electroweak Theory.  _Phys. Rev. D_ , 36, p. 581.  10.1103/PhysRevD.36.581.
  * Barenboim et al. (2017) Barenboim, G., W. H. Kinney, and W.-I. Park (2017).  Resurrection of large lepton number asymmetries from neutrino flavor oscillations.  _Phys. Rev. D_ , 95(4), p. 043506.  10.1103/PhysRevD.95.043506.
  * Barenboim and Park (2017) Barenboim, G. and W.-I. Park (2017).  A full picture of large lepton number asymmetries of the Universe.  _JCAP_ , 04, p. 048.  10.1088/1475-7516/2017/04/048.
  * Bazavov et al. (2018) Bazavov, A., N. Brambilla, P. Petreczky, A. Vairo, and J. H. Weber (2018).  Color screening in (2+1)-flavor QCD.  _Phys. Rev. D_ , 98(5), p. 054511.  10.1103/PhysRevD.98.054511.
  * Birrell (2014) Birrell, J. (2014).  _Non-Equilibrium Aspects of Relic Neutrinos: From Freeze-out to the Present Day_.  Ph.D. thesis, Arizona U.
  * Birrell and Rafelski (2015) Birrell, J. and J. Rafelski (2015).  Quark–gluon plasma as the possible source of cosmological dark radiation.  _Phys. Lett. B_ , 741, pp. 77–81.  10.1016/j.physletb.2014.12.033.
  * Birrell et al. (2015) Birrell, J., J. Wilkening, and J. Rafelski (2015).  Boltzmann Equation Solver Adapted to Emergent Chemical Non-equilibrium.  _J. Comput. Phys._ , 281, pp. 896–916.  10.1016/j.jcp.2014.10.056.
  * Birrell et al. (2013) Birrell, J., C.-T. Yang, P. Chen, and J. Rafelski (2013).  Fugacity and Reheating of Primordial Neutrinos.  _Mod. Phys. Lett. A_ , 28, p. 1350188.  10.1142/S0217732313501885.
  * Birrell et al. (2014a) Birrell, J., C.-T. Yang, P. Chen, and J. Rafelski (2014a).  Relic neutrinos: Physically consistent treatment of effective number of neutrinos and neutrino mass.  _Phys. Rev. D_ , 89, p. 023008.  10.1103/PhysRevD.89.023008.
  * Birrell et al. (2014b) Birrell, J., C.-T. Yang, and J. Rafelski (2014b).  Relic Neutrino Freeze-out: Dependence on Natural Constants.  _Nucl. Phys. B_ , 890, pp. 481–517.  10.1016/j.nuclphysb.2014.11.020.
  * Blennow et al. (2012) Blennow, M., E. Fernandez-Martinez, O. Mena, J. Redondo, and P. Serra (2012).  Asymmetric Dark Matter and Dark Radiation.  _JCAP_ , 07, p. 022.  10.1088/1475-7516/2012/07/022.
  * Boehm et al. (2012) Boehm, C., M. J. Dolan, and C. McCabe (2012).  Increasing Neff with particles in thermal equilibrium with neutrinos.  _JCAP_ , 12, p. 027.  10.1088/1475-7516/2012/12/027.
  * Brambilla et al. (2010) Brambilla, N., M. A. Escobedo, J. Ghiglieri, J. Soto, and A. Vairo (2010).  Heavy Quarkonium in a weakly-coupled quark-gluon plasma below the melting temperature.  _JHEP_ , 09, p. 038.  10.1007/JHEP09(2010)038.
  * Brambilla et al. (2018) Brambilla, N., M. A. Escobedo, J. Soto, and A. Vairo (2018).  Heavy quarkonium suppression in a fireball.  _Phys. Rev. D_ , 97(7), p. 074009.  10.1103/PhysRevD.97.074009.
  * Carraro et al. (1988) Carraro, C., A. Schafer, and S. E. Koonin (1988).  Dynamic Screening of Thermonuclear Reactions.  _ApJ_ , 331, p. 565.  10.1086/166582.
  * Choquet-Bruhat (2008) Choquet-Bruhat, Y. (2008).  _General relativity and the Einstein equations_.  OUP Oxford.
  * Coc et al. (2006) Coc, A., K. A. Olive, J.-P. Uzan, and E. Vangioni (2006).  Big bang nucleosynthesis constraints on scalar-tensor theories of gravity.  _Phys. Rev. D_ , 73, p. 083525.  10.1103/PhysRevD.73.083525.
  * Cugnon and Lombard (1984) Cugnon, J. and R. M. Lombard (1984).  K+ PRODUCTION IN A CASCADE MODEL FOR HIGH-ENERGY NUCLEUS NUCLEUS COLLISIONS.  _Nucl. Phys. A_ , 422, pp. 635–653.  10.1016/0375-9474(84)90369-5.
  * Czarnecki et al. (2004) Czarnecki, A., W. J. Marciano, and A. Sirlin (2004).  Precision measurements and CKM unitarity.  _Phys. Rev. D_ , 70, p. 093006.  10.1103/PhysRevD.70.093006.
  * Davidson et al. (2008) Davidson, S., E. Nardi, and Y. Nir (2008).  Leptogenesis.  _Phys. Rept._ , 466, pp. 105–177.  10.1016/j.physrep.2008.06.002.
  * Dicus et al. (1982) Dicus, D. A., E. W. Kolb, A. M. Gleeson, E. C. G. Sudarshan, V. L. Teplitz, and M. S. Turner (1982).  Primordial Nucleosynthesis Including Radiative, Coulomb, and Finite Temperature Corrections to Weak Rates.  _Phys. Rev. D_ , 26, p. 2694.  10.1103/PhysRevD.26.2694.
  * Dodelson (2003) Dodelson, S. (2003).  _Modern Cosmology_.  Academic Press, Amsterdam.  ISBN 978-0-12-219141-1.
  * Dreiner et al. (2012) Dreiner, H. K., M. Hanussek, J. S. Kim, and S. Sarkar (2012).  Gravitino cosmology with a very light neutralino.  _Phys. Rev. D_ , 85, p. 065027.  10.1103/PhysRevD.85.065027.
  * Elmfors et al. (1997) Elmfors, P., K. Enqvist, G. Raffelt, and G. Sigl (1997).  Neutrinos with magnetic moment: Depolarization rate in plasma.  _Nucl. Phys. B_ , 503, pp. 3–23.  10.1016/S0550-3213(97)00382-9.
  * Elze et al. (1980) Elze, H. T., W. Greiner, and J. Rafelski (1980).  The relativistic ideal Fermi gas revisited.  _J. Phys. G_ , 6, pp. L149–L153.  10.1088/0305-4616/6/9/003.
  * Famiano et al. (2016) Famiano, M. A., A. B. Balantekin, and T. Kajino (2016).  Low-lying Resonances and Relativistic Screening in Big Bang Nucleosynthesis.  _Phys. Rev. C_ , 93(4), p. 045804.  10.1103/PhysRevC.93.045804.
  * Fields (2011) Fields, B. D. (2011).  The primordial lithium problem.  _Ann. Rev. Nucl. Part. Sci._ , 61, pp. 47–68.  10.1146/annurev-nucl-102010-130445.
  * Formanek et al. (2018) Formanek, M., S. Evans, J. Rafelski, A. Steinmetz, and C.-T. Yang (2018).  Strong fields and neutral particle magnetic moment dynamics.  _Plasma Phys. Control. Fusion_ , 60, p. 074006.  10.1088/1361-6587/aac06a.
  * Formanek et al. (2021) Formanek, M., C. Grayson, J. Rafelski, and B. Müller (2021).  Current-conserving relativistic linear response for collisional plasmas.  _Annals Phys._ , 434, p. 168605.  10.1016/j.aop.2021.168605.
  * Fornengo et al. (1997) Fornengo, N., C. W. Kim, and J. Song (1997).  Finite temperature effects on the neutrino decoupling in the early universe.  _Phys. Rev. D_ , 56, pp. 5123–5134.  10.1103/PhysRevD.56.5123.
  * Freedman et al. (1977) Freedman, D. Z., D. N. Schramm, and D. L. Tubbs (1977).  The Weak Neutral Current and Its Effects in Stellar Collapse.  _Ann. Rev. Nucl. Part. Sci._ , 27, pp. 167–207.  10.1146/annurev.ns.27.120177.001123.
  * Fromerth et al. (2012) Fromerth, M. J., I. Kuznetsova, L. Labun, J. Letessier, and J. Rafelski (2012).  From Quark-Gluon Universe to Neutrino Decoupling: 200 $<$ T $<$ 2MeV.  _Acta Phys. Polon. B_ , 43(12), pp. 2261–2284.  10.5506/APhysPolB.43.2261.
  * Fromerth and Rafelski (2002) Fromerth, M. J. and J. Rafelski (2002).  Hadronization of the quark Universe.  _arXiv astro-ph/0211346_.
  * Fukugita and Yazaki (1987) Fukugita, M. and S. Yazaki (1987).  Reexamination of Astrophysical and Cosmological Constraints on the Magnetic Moment of Neutrinos.  _Phys. Rev. D_ , 36, p. 3817.  10.1103/PhysRevD.36.3817.
  * Giovannini (2004) Giovannini, M. (2004).  The Magnetized universe.  _Int. J. Mod. Phys. D_ , 13, pp. 391–502.  10.1142/S0218271804004530.
  * Giudice et al. (2004) Giudice, G. F., A. Notari, M. Raidal, A. Riotto, and A. Strumia (2004).  Towards a complete theory of thermal leptogenesis in the SM and MSSM.  _Nucl. Phys. B_ , 685, pp. 89–149.  10.1016/j.nuclphysb.2004.02.019.
  * Giunti and Kim (2007) Giunti, C. and C. W. Kim (2007).  _Fundamentals of Neutrino Physics and Astrophysics_.  Oxford University Press.  ISBN 978-0-19-850871-7.
  * Giunti and Studenikin (2009) Giunti, C. and A. Studenikin (2009).  Neutrino electromagnetic properties.  _Phys. Atom. Nucl._ , 72, pp. 2089–2125.  10.1134/S1063778809120126.
  * Grayson et al. (2023) Grayson, C., C. T. Yang, M. Formanek, and J. Rafelski (2023).  Electron–positron plasma in BBN: Damped-dynamic screening.  _Annals Phys._ , 458, p. 169453.  10.1016/j.aop.2023.169453.
  * Griffiths (2008) Griffiths, D. (2008).  _Introduction to elementary particles_.  ISBN 978-3-527-40601-2.
  * Gruzinov (1998) Gruzinov, A. V. (1998).  Dynamic screening and thermonuclear reaction rates.  _Astrophys. J._ , 496, p. 503.  10.1086/305349.
  * Heckler (1994) Heckler, A. F. (1994).  Astrophysical applications of quantum corrections to the equation of state of a plasma.  _Phys. Rev. D_ , 49, pp. 611–617.  10.1103/PhysRevD.49.611.
  * Hwang et al. (2021) Hwang, E., D. Jang, K. Park, M. Kusakabe, T. Kajino, A. B. Balantekin, T. Maruyama, C.-M. Ryu, and M.-K. Cheoun (2021).  Dynamical screening effects on big bang nucleosynthesis.  _JCAP_ , 11, p. 017.  10.1088/1475-7516/2021/11/017.
  * Jedamzik and Saveliev (2019) Jedamzik, K. and A. Saveliev (2019).  Stringent limit on primordial magnetic fields from the cosmic microwave background radiation.  _Physical review letters_ , 123(2), p. 021301.  10.1103/PhysRevLett.123.021301.
  * Karsch et al. (1988) Karsch, F., M. T. Mehr, and H. Satz (1988).  Color Screening and Deconfinement for Bound States of Heavy Quarks.  _Z. Phys. C_ , 37, p. 617.  10.1007/BF01549722.
  * Kislinger and Morley (1976) Kislinger, M. B. and P. D. Morley (1976).  Collective Phenomena in Gauge Theories. 1. The Plasmon Effect for Yang-Mills Fields.  _Phys. Rev. D_ , 13, p. 2765.  10.1103/PhysRevD.13.2765.
  * Koch et al. (1986) Koch, P., B. Muller, and J. Rafelski (1986).  Strangeness in Relativistic Heavy Ion Collisions.  _Phys. Rept._ , 142, pp. 167–262.  10.1016/0370-1573(86)90096-7.
  * Kolb et al. (1996) Kolb, E. W., A. D. Linde, and A. Riotto (1996).  GUT baryogenesis after preheating.  _Phys. Rev. Lett._ , 77, pp. 4290–4293.  10.1103/PhysRevLett.77.4290.
  * Kolb and Turner (1990a) Kolb, E. W. and M. S. Turner (1990a).  _The early universe_.  CRC Press.
  * Kolb and Turner (1990b) Kolb, E. W. and M. S. Turner (1990b).  _The Early Universe_ , volume 69.  ISBN 978-0-201-62674-2.  10.1201/9780429492860.
  * Kronberg (1994a) Kronberg, P. P. (1994a).  Extragalactic magnetic fields.  _Rept. Prog. Phys._ , 57, pp. 325–382.  10.1088/0034-4885/57/4/001.
  * Kronberg (1994b) Kronberg, P. P. (1994b).  Extragalactic magnetic fields.  _Reports on Progress in Physics_ , 57(4), p. 325.  10.1088/0034-4885/57/4/001.
  * Kuzmin et al. (1985) Kuzmin, V. A., V. A. Rubakov, and M. E. Shaposhnikov (1985).  On the Anomalous Electroweak Baryon Number Nonconservation in the Early Universe.  _Phys. Lett. B_ , 155, p. 36.  10.1016/0370-2693(85)91028-7.
  * Kuzmin et al. (1987) Kuzmin, V. A., V. A. Rubakov, and M. E. Shaposhnikov (1987).  Anomalous Electroweak Baryon Number Nonconservation and GUT Mechanism for Baryogenesis.  _Phys. Lett. B_ , 191, pp. 171–173.  10.1016/0370-2693(87)91340-2.
  * Kuznetsova et al. (2008) Kuznetsova, I., D. Habs, and J. Rafelski (2008).  Pion and muon production in e-, e+, gamma plasma.  _Phys. Rev. D_ , 78, p. 014027.  10.1103/PhysRevD.78.014027.
  * Kuznetsova et al. (2010) Kuznetsova, I., D. Habs, and J. Rafelski (2010).  Thermal reaction processes in a relativistic QED plasma drop.  _Phys. Rev. D_ , 81, p. 053007.  10.1103/PhysRevD.81.053007.
  * Kuznetsova and Rafelski (2010a) Kuznetsova, I. and J. Rafelski (2010a).  Unstable Hadrons in Hot Hadron Gas in Laboratory and in the Early Universe.  _Phys. Rev. C_ , 82, p. 035203.  10.1103/PhysRevC.82.035203.
  * Kuznetsova and Rafelski (2010b) Kuznetsova, I. and J. Rafelski (2010b).  Unstable hadrons in hot hadron gas: In the laboratory and in the early Universe.  _Phys. Rev. C_ , 82, p. 035203.  10.1103/PhysRevC.82.035203.
  * Kuznetsova and Rafelski (2012) Kuznetsova, I. and J. Rafelski (2012).  Electron-Positron Plasma Drop Formed by Ultra-Intense Laser Pulses.  _Phys. Rev. D_ , 85, p. 085014.  10.1103/PhysRevD.85.085014.
  * Letessier and Rafelski (2002) Letessier, J. and J. Rafelski (2002).  _Hadrons and Quark–Gluon Plasma_.  Oxford University Press.  ISBN 978-1-00-929075-3, 978-0-511-53499-7, 978-1-00-929070-8, 978-1-00-929073-9.  10.1017/9781009290753.
  * Letessier and Rafelski (2008) Letessier, J. and J. Rafelski (2008).  Hadron production and phase changes in relativistic heavy ion collisions.  _Eur. Phys. J. A_ , 35, pp. 221–242.  10.1140/epja/i2007-10546-7.
  * Lewis (1980) Lewis, R. R. (1980).  Coherent detector for low-energy neutrinos.  _Phys. Rev. D_ , 21, pp. 663–668.  10.1103/PhysRevD.21.663.
  * Mangano et al. (2002) Mangano, G., G. Miele, S. Pastor, and M. Peloso (2002).  A Precision calculation of the effective number of cosmological neutrinos.  _Phys. Lett. B_ , 534, pp. 8–16.  10.1016/S0370-2693(02)01622-2.
  * Mangano et al. (2005) Mangano, G., G. Miele, S. Pastor, T. Pinto, O. Pisanti, and P. D. Serpico (2005).  Relic neutrino decoupling including flavor oscillations.  _Nucl. Phys. B_ , 729, pp. 221–234.  10.1016/j.nuclphysb.2005.09.041.
  * Mangano et al. (2006) Mangano, G., G. Miele, S. Pastor, T. Pinto, O. Pisanti, and P. D. Serpico (2006).  Effects of non-standard neutrino-electron interactions on relic neutrino decoupling.  _Nucl. Phys. B_ , 756, pp. 100–116.  10.1016/j.nuclphysb.2006.09.002.
  * Morgan (1981) Morgan, J. A. (1981).  Cosmological upper limit to neutrino magnetic moments.  _Phys. Lett. B_ , 102, pp. 247–250.  10.1016/0370-2693(81)90868-6.
  * Morrissey and Ramsey-Musolf (2012) Morrissey, D. E. and M. J. Ramsey-Musolf (2012).  Electroweak baryogenesis.  _New J. Phys._ , 14, p. 125003.  10.1088/1367-2630/14/12/125003.
  * Mukhanov (2005) Mukhanov, V. (2005).  _Physical Foundations of Cosmology_.  Cambridge University Press, Oxford.  ISBN 978-0-521-56398-7.  10.1017/CBO9780511790553.
  * Neronov and Vovk (2010) Neronov, A. and I. Vovk (2010).  Evidence for strong extragalactic magnetic fields from Fermi observations of TeV blazars.  _Science_ , 328(5974), pp. 73–75.  10.1126/science.1184192.
  * Nicolescu (2013) Nicolescu, G. (2013).  A heuristic approach to the detection of solar neutrinos.  _Journal of Physics G: Nuclear and Particle Physics_ , 40, p. 055201.  10.1088/0954-3899/40/5/055201.
  * Nielsen and Takanishi (2001) Nielsen, H. B. and Y. Takanishi (2001).  Baryogenesis via lepton number violation in anti-GUT model.  _Phys. Lett. B_ , 507, pp. 241–251.  10.1016/S0370-2693(01)00357-4.
  * Offler et al. (2019) Offler, S., G. Aarts, C. Allton, J. Glesaaen, B. Jäger, S. Kim, M. P. Lombardo, S. M. Ryan, and J.-I. Skullerud (2019).  News from bottomonium spectral functions in thermal QCD.  _PoS_ , LATTICE2019, p. 076.  10.22323/1.363.0076.
  * Papavassiliou et al. (2006) Papavassiliou, J., J. Bernabeu, and M. Passera (2006).  Neutrino-nuclear coherent scattering and the effective neutrino charge radius.  _PoS_ , HEP2005, p. 192.  10.22323/1.021.0192.
  * Patrignani et al. (2016) Patrignani, C. et al. (2016).  Review of Particle Physics.  _Chin. Phys. C_ , 40(10), p. 100001.  10.1088/1674-1137/40/10/100001.
  * Pitrou et al. (2018) Pitrou, C., A. Coc, J.-P. Uzan, and E. Vangioni (2018).  Precision big bang nucleosynthesis with improved Helium-4 predictions.  _Phys. Rept._ , 754, pp. 1–66.  10.1016/j.physrep.2018.04.005.
  * Pshirkov et al. (2016) Pshirkov, M. S., P. G. Tinyakov, and F. R. Urban (2016).  New limits on extragalactic magnetic fields from rotation measures.  _Phys. Rev. Lett._ , 116(19), p. 191302.  10.1103/PhysRevLett.116.191302.
  * Rafelski (2020) Rafelski, J. (2020).  Discovery of Quark-Gluon-Plasma: Strangeness Diaries.  _Eur. Phys. J. ST_ , 229(1), pp. 1–140.  10.1140/epjst/e2019-900263-x.
  * Rafelski and Birrell (2014) Rafelski, J. and J. Birrell (2014).  Traveling Through the Universe: Back in Time to the Quark-Gluon Plasma Era.  _J. Phys. Conf. Ser._ , 509, p. 012014.  10.1088/1742-6596/509/1/012014.
  * Rafelski et al. (2023) Rafelski, J., J. Birrell, A. Steinmetz, and C. T. Yang (2023).  A short survey of matter-antimatter evolution in the primordial universe.
  * Rafelski et al. (2001) Rafelski, J., J. Letessier, and G. Torrieri (2001).  Strange hadrons and their resonances: A Diagnostic tool of QGP freezeout dynamics.  _Phys. Rev. C_ , 64, p. 054907.  10.1103/PhysRevC.64.054907.  [Erratum: Phys.Rev.C 65, 069902 (2002)].
  * Rafelski and Yang (2021) Rafelski, J. and C. T. Yang (2021).  The muon abundance in the primordial Universe.  _Acta Phys. Polon. B_ , 52, p. 277.  10.5506/APhysPolB.52.277.
  * Rafelski and Yang (2022) Rafelski, J. and C. T. Yang (2022).  Reactions Governing Strangeness Abundance in Primordial Universe.  _EPJ Web Conf._ , 259, p. 13001.  10.1051/epjconf/202225913001.
  * Riess et al. (2018a) Riess, A. G., S. Casertano, W. Yuan, L. Macri, J. Anderson, J. W. MacKenty, J. B. Bowers, K. I. Clubb, A. V. Filippenko, D. O. Jones, et al. (2018a).  New parallaxes of galactic cepheids from spatially scanning the hubble space telescope: Implications for the hubble constant.  _The Astrophysical Journal_ , 855(2), p. 136.
  * Riess et al. (2018b) Riess, A. G. et al. (2018b).  Milky Way Cepheid Standards for Measuring Cosmic Distances and Application to Gaia DR2: Implications for the Hubble Constant.  _Astrophys. J._ , 861(2), p. 126.  10.3847/1538-4357/aac82e.
  * Riotto and Trodden (1999) Riotto, A. and M. Trodden (1999).  Recent progress in baryogenesis.  _Ann. Rev. Nucl. Part. Sci._ , 49, pp. 35–75.  10.1146/annurev.nucl.49.1.35.
  * Sajjad Athar et al. (2023) Sajjad Athar, M., A. Fatima, and S. K. Singh (2023).  Neutrinos and their interactions with matter.  _Prog. Part. Nucl. Phys._ , 129, p. 104019.  10.1016/j.ppnp.2022.104019.
  * Sakharov (1967) Sakharov, A. D. (1967).  Violation of CP Invariance, C asymmetry, and baryon asymmetry of the universe.  _Pisma Zh. Eksp. Teor. Fiz._ , 5, pp. 32–35.  10.1070/PU1991v034n05ABEH002497.
  * Sakharov (1990) Sakharov, A. D. (1990).  Baryon asymmetry of the universe.  pp. 65–80.  10.1070/PU1991v034n05ABEH002504.
  * Salpeter (1954) Salpeter, E. E. (1954).  Electron screening and thermonuclear reactions.  _Austral. J. Phys._ , 7, pp. 373–388.  10.1071/PH540373.
  * Salpeter and van Horn (1969) Salpeter, E. E. and H. M. van Horn (1969).  Nuclear Reaction Rates at High Densities.  _ApJ_ , 155, p. 183.  10.1086/149858.
  * Schroedter et al. (2000) Schroedter, M., R. L. Thews, and J. Rafelski (2000).  $B_{c}$ meson production in nuclear collisions at RHIC.  _Phys. Rev. C_ , 62, p. 024905.  10.1103/PhysRevC.62.024905.
  * Sehgal and Wanninger (1986) Sehgal, L. M. and M. Wanninger (1986).  Atomic Effects in Coherent Neutrino Scattering.  _Phys. Lett. B_ , 171, pp. 107–112.  10.1016/0370-2693(86)91008-7.
  * Smith (1984) Smith, P. F. (1984).  Coherent neutrino scattering—Relativistic and nonrelativistic.  _II Nuovo Cimento A (1965-1970)_ , 83, pp. 263–274.  10.1007/BF02902601.
  * Steinmetz et al. (2023a) Steinmetz, A., C. T. Yang, and J. Rafelski (2023a).  Antimatter origin of cosmic magnetism: A first look.  _(in preparation)_.
  * Steinmetz et al. (2023b) Steinmetz, A., C. T. Yang, and J. Rafelski (2023b).  Matter-antimatter origin of cosmic magnetism.
  * Tanabashi et al. (2018) Tanabashi, M. et al. (2018).  Review of Particle Physics.  _Phys. Rev. D_ , 98(3), p. 030001.  10.1103/PhysRevD.98.030001.
  * Taylor et al. (2011) Taylor, A. M., I. Vovk, and A. Neronov (2011).  Extragalactic magnetic fields constraints from simultaneous GeV–TeV observations of blazars.  _Astronomy & Astrophysics_, 529, p. A144.  10.1051/0004-6361/201116441.
  * Tully (2019) Tully, A. M. (2019).  _Doubly charmed $B$ decays with the LHCb experiment_.  Ph.D. thesis, Cambridge U.  10.17863/CAM.44796.
  * Vazza et al. (2021) Vazza, F. et al. (2021).  Magnetogenesis and the Cosmic Web: A Joint Challenge for Radio Observations and Numerical Simulations.  _Galaxies_ , 9(4), p. 109.  10.3390/galaxies9040109.
  * Vernstrom et al. (2021) Vernstrom, T., G. Heald, F. Vazza, T. J. Galvin, J. L. West, N. Locatelli, N. Fornengo, and E. Pinetti (2021).  Discovery of magnetic fields along stacked cosmic filaments as revealed by radio and X-ray emission.  _Monthly Notices of the Royal Astronomical Society_ , 505(3), pp. 4178–4196.  10.1093/mnras/stab1301.
  * Vogel and Engel (1989) Vogel, P. and J. Engel (1989).  Neutrino Electromagnetic Form-Factors.  _Phys. Rev. D_ , 39, p. 3378.  10.1103/PhysRevD.39.3378.
  * Wang et al. (2011) Wang, B., C. A. Bertulani, and A. B. Balantekin (2011).  Electron screening and its effects on Big-Bang nucleosynthesis.  _Phys. Rev. C_ , 83, p. 018801.  10.1103/PhysRevC.83.018801.
  * Weber (1988) Weber, J. (1988).  Apparent observation of abnormally large coherent scattering cross sections using keV and MeV range antineutrinos, and solar neutrinos.  _Phys. Rev. D_ , 38, pp. 32–39.  10.1103/PhysRevD.38.32.
  * Widrow et al. (2012) Widrow, L. M., D. Ryu, D. R. G. Schleicher, K. Subramanian, C. G. Tsagas, and R. A. Treumann (2012).  The First Magnetic Fields.  _Space Sci. Rev._ , 166, pp. 37–70.  10.1007/s11214-011-9833-5.
  * Wolfenstein (1978) Wolfenstein, L. (1978).  Neutrino oscillations in matter.  _Phys. Rev. D_ , 17, pp. 2369–2374.  10.1103/PhysRevD.17.2369.
  * Workman et al. (2022) Workman, R. L. et al. (2022).  Review of Particle Physics.  _PTEP_ , 2022, p. 083C01.  10.1093/ptep/ptac097.
  * Yang et al. (2018a) Yang, C. T., J. Birrell, and J. Rafelski (2018a).  Lepton Number and Expansion of the Universe.
  * Yang et al. (2018b) Yang, C. T., J. Birrell, and J. Rafelski (2018b).  Temperature Dependence of the Neutron Lifespan.
  * Yang and Rafelski (2020) Yang, C. T. and J. Rafelski (2020).  Possibility of bottom-catalyzed matter genesis near to primordial QGP hadronization.
  * Yang and Rafelski (2022) Yang, C. T. and J. Rafelski (2022).  Cosmological strangeness abundance.  _Phys. Lett. B_ , 827, p. 136944.  10.1016/j.physletb.2022.136944.