# Coherence Control of Directional Nonlinear Photocurrent
in Spatially Symmetric Systems

Yongliang Shi  Center for Spintronics and Quantum System, State Key Laboratory
for Mechanical Behavior of Materials, Xi’an Jiaotong University, Xi’an,
710049, China  Jian Zhou  Email<EMAIL_ADDRESS>Center for Alloy
Innovation and Design, State Key Laboratory for Mechanical Behavior of
Materials, Xi’an Jiaotong University, Xi’an, 710049, China

###### Abstract

The interplay between crystal symmetry and its optical responses is at the
heart of tremendous recent advances in light-matter interactions and
applications. Nonlinear optical processes that produce electric currents, for
example bulk photovoltaic (BPV) effect, require inversion symmetry broken
materials, such as ferroelectrics. In the current work, we demonstrate that
such BPV current could be generated in centrosymmetric materials with
excitation of out-of-equilibrium coherent phonons. This is much different from
the generally studied static or thermally excited states. We show that
depending on the oscillating phase factor of the coherent phonon, uni-
directional static electric current can be observed, in addition to some
terahertz alternating currents. We also generalize the conventional injection
charge current into angular momentum (spin and orbital) degrees of freedom,
and demonstrate spin and orbital BPV photocurrents under coherent phonons. Our
findings open the pathway to exploring the exotic phonon-photon-electron
coherent interactions in quantum materials.

##  .1 Introduction

Light-matter interaction is one of the key routes to manipulating electrons
and converting light energy into electrical energy. Conventional light energy
harvesting and conversion into electric current relies on semiconductor
heterostructures (such as p-n junctions), in which electrons and holes can be
separated spatially under a built-in electric field. This requirement can be
eliminated when one uses a single crystal with non-centrosymmetric structure
($\mathcal{P}$-broken), and light illumination could generate direct electric
current, known as bulk photovoltaic (BPV) effect Fridkin (2001); Aversa and
Sipe (1995). The inversion asymmetry yields anharmonic electric field
potential, leading to a collective motion of electrons under light irradiation
(when photon energy $\hbar\omega$ is typically larger than semiconductor
bandgap $E_{g}$). Here, $\mathcal{P}$-broken is an acknowledged condition to
generate bias-free photocurrent.

The requirement of $\mathcal{P}$-broken limits the candidate material set to
produce photocurrents (including shift current and injection current,
depending on the crystal symmetry and light polarization). Both shift current
and injection current processes have been widely observed in various
materials, including ferroelectrics, quantum wells, organic crystals, and two-
dimensional (2D) interfaces. One has to note that most previous works focus on
ferroelectrics in their static or thermally equilibrium state, which are
inherently $\mathcal{P}$-broken Gong  _et al._ (2018); Tan  _et al._ (2016).
In this work, we suggest that injection photocurrent can be triggered in
spatially symmetric materials with coherent phonon oscillation excited using a
THz pulse irradiation. We perform first-principles density functional theory
calculations Hohenberg and Kohn (1964); Kohn and Sham (1965); Kresse and
Furthmüller (1996); Perdew  _et al._ (1996); Blöchl (1994) and discuss a
spatially centrosymmetric system, 2D phosphorene ($\alpha$-P) monolayer Liu
_et al._ (2014) to illustrate our theory. The 2D materials are advantageous
for their easy optical access and manipulation Zhou  _et al._ (2018); Guo  _et
al._ (2021); Xu  _et al._ (2021a), owing to their ultrathin nature. In
addition, the direct optical absorption is much reduced in 2D thin films,
hence the hot carrier mobility drift current can be marginal, yielding better
experimental observations. We show that the zone-center ($q=0$) infrared (IR)
active phonon can be excited under an THz laser pulse with intermediate
fluence, which brings coherent odd-parity mode oscillations. Under such
dynamic state, circularly polarized light (CPL) irradiation could produce a
nonzero unidirectional propagating current, according to injection
photocurrent mechanism. One has to note that even though the average
displacements of such coherent phonon vibration averages to be zero, which
corresponds to the static centrosymmetric geometry, the photon irradiation
induced injection current does not vanish and contains a static current
component. Such counter-intuitive phenomenon arises due to the injection
current microscopic mechanism, which reflects the velocity mismatch between
the valence and conduction bands. The direction of such photocurrent can be
determined by the phonon vibration phase when the CPL irradiates, thus a
coherence control of charge current can be realized Zhao  _et al._ (2006).
Furthermore, we propose that under linearly polarized light (LPL) irradiation,
one could generate angular momentum (AM) related photocurrent (both spin BPV
current and orbital BPV current), which has similar injection current nature.
This generalizes the photo-induced nanoelectronics into spin and orbital AM
degrees of freedom, which carries magnetic moment information. Therefore,
dynamic but nonzero photocurrent with tunable transport direction (rather than
a pure THz alternating current which averages to be zero over time) and with
charge, spin, and orbital information can be produced via a synergistic
cooperation of coherent phonon and second order nonlinear injection BPV
approaches. Such mechanism can be extended to other spatially (vertical
mirror) symmetric materials, such as H-MoS2 monolayer Wang  _et al._ (2019).

Figure 1: (a) THz optics irradiation inducing coherent IR-active phonon
vibration in $\alpha$-P monolayer. The mode displacement amplitude is enlarged
for clarity reason. (b) Nonlinear injection photocurrent move along the
$x$-(zigzag) direction. (c) Quasi-static band dispersion variation under
$\lambda$.

According to nonlinear optics theory, once a CPL irradiates onto a nonmagnetic
($\mathcal{T}$-symmetric) material, its helicity phase generates injection
current that is proportional to the interband Berry curvature (usually known
as circular photogalvanic effect Ji  _et al._ (2019)). In this process, light
excites carriers into velocity-asymmetrical states in the Brillouin zone (BZ)
and produces a net charge current Sipe and Shkrebtii (2000); Fregoso (2019);
Wang and Qian (2019); Holder  _et al._ (2020); Taghizadeh  _et al._ (2017).
The photocurrent conductivity $\sigma(0;\omega,-\omega)$ grows linearly with
time, and its rate is

$\begin{split}\eta_{bc}^{a}(\omega)&=\frac{d}{dt}\sigma_{bc}^{a}(0;\omega,-\omega)=-\frac{\pi
e^{3}}{2\hbar^{2}}\int_{BZ}\frac{d^{2}\bm{k}}{(2\pi)^{2}}\\\
&\times\mathrm{Im}\sum_{n,m}f_{nm}\Delta_{nm}^{a}[r_{mn}^{b},r_{nm}^{c}]\delta(\omega_{mn}-\omega)\end{split}$
(1)

Here $a,b,c$ indicate 2D Cartesian indices $x$ and $y$. $f_{nm}=f_{n}-f_{m}$
and $\Delta_{nm}^{a}=v_{nn}^{a}-v_{mm}^{a}$ are the occupation and velocity
difference between band-$n$ and $m$, respectively. Interband transition dipole
is $r_{mn}^{b}=\langle m|r^{b}|n\rangle$. Note that all explicit
$\bm{k}$-dependence of quantities are omitted. The commutator
$[r_{mn}^{b},r_{nm}^{c}]=r_{mn}^{b}r_{nm}^{c}-r_{mn}^{c}r_{nm}^{b}=-i\epsilon_{abc}\Omega_{mn}^{a}$
corresponds to interband Berry curvature. We use Lorentzian function to
represent the Dirac $\delta$-function, with a broadening factor of
$0.02\,\mathrm{eV}$. It could phenomenologically incorporate the temperature,
disorder, and scattering effects. The total current evolves as
$j_{\mathrm{inj}}^{a}=2\sigma_{bc}^{a}(0;\omega,-\omega)E_{b}(\omega)E_{c}(-\omega)$,
where $E_{i}$ indicates optical electric field component. According to Eq. 1,
the different velocities at the valence and conduction bands produces the net
injection photocurrent. For the 2D BZ integration which lacks the out-of-plane
$k_{z}$ component, the calculated $\eta$ (and $\sigma$) includes a surplus
length unit (nm or Å) compared with that in conventional 3D crystals. In order
to eliminate such difference and adopt the 3D unit systems, we rescale them by
utilizing an effective thickness of the system,
$\eta^{\mathrm{eff}}=\eta^{\mathrm{2D}}/d$. Here $d$ refers to the monolayer
effective thickness (taken to be $6\,\mathrm{\AA}$ in the current work), and
the superscripts “2D” and “eff” indicate 2D values and the rescaled effective
results. This is based on a parallel capacitor model and is widely used
previously Zhou  _et al._ (2018); Laturia  _et al._ (2018); Jiang  _et al._
(2017); Mu  _et al._ (2021). In this work, we report these rescaled values.

##  .2 Results

Coherent vibration of $\alpha$-P monolayer and injection current rate under
CPL. According to symmetry arguments, we limit our discussion on the zone-
centered phonon modes (unit cell wavelength coherent oscillation), which can
be effectively excited via THz optics (Supplemental Material, SM). In
$\alpha$-P monolayer, the only one is $B_{1u}$ (at
$4.2\,\mathrm{THz}=140\,\mathrm{cm}^{-1}$), which breaks $\mathcal{P}$ and
introduces dynamic in-plane polarization ($P_{y}(t)\neq 0$, see SM). This mode
mainly displaces P atoms out-of-plane, as shown in Fig. 1(a). We find that the
anharmonicity of $B_{1u}$ is very weak (see SM), suggesting its small damping
effect and long lifetime during operation. Now we displace atoms according to
their phonon eigenvectors, and use $\lambda$ to measure the total
displacements of all four atoms in a unit cell. When $\lambda=\pm
0.1\,\mathrm{\AA}$, the energy gain is $10\,\mathrm{meV}$ per atom. This
corresponds to a THz laser pulse with pump fluence of
$0.25\,\mathrm{mJ/cm^{2}}$ centered at $4.2\,\mathrm{THz}$ (see SM). The
positive and negative $\lambda$ correspond to spatially opposite images, i.e.,
atomic coordinates $\mathbf{x}_{i}^{\lambda}=-\mathbf{x}_{i}^{-\lambda}$ for
each atom $i$. Hence, they show same electronic energy band structure [Fig.
1(c)].

Figure 2: (a) Injection photocurrent rate $\eta_{xy}^{x}$ of $\alpha$-P in
various $\lambda$ states. (b) $k$-resolved contribution of injection
photocurrent rate at $\omega_{0}=2.7\,\mathrm{eV}$. (c) Linear fitting of
$\eta_{xy}^{x}(\omega_{0})$ with $\lambda$. Figure 3: (a) Spin and (b) orbital
photocurrent rate under different $\lambda$. (c) and (d) show their
$\lambda$-dependence at $\omega_{0}=2.7\,\mathrm{eV}$. Insets show the
$k$-resolved contributions at $\lambda=-0.1\,\mathrm{\AA}$.

The coherent $B_{1u}$ reduces the space group from static $Pmna$ to
$Pmn2_{1}$, which breaks the $\mathcal{M}_{y}$ mirror symmetry but preserves
$\mathcal{M}_{x}$. According to our symmetry analysis (see SM), dynamic
injection current then emerge under CPL and would flow along the
$x$-direction. Here we will only discuss $\eta_{xy}^{x}$ since
$\eta_{yx}^{x}=-\eta_{xy}^{x}$. Since the ionic vibration periodic is on a few
hundreds of picosecond, one can safely assume that electrons effectively
experience a quasi-static potential field under different $\lambda$. Figure
2(a) shows the injection current rate as a function of different displacement.
One clearly sees that positive and negative $\lambda$ give opposite
photoconductivity rates ($\eta_{xy}^{x}(\lambda)=-\eta_{xy}^{x}(-\lambda)$),
even though their band dispersions are the same. We also plot the $k$-resolved
integrand of Eq. 1 at a selected optical incident energy
$\omega_{0}=2.7\,\mathrm{eV}$, which shows large contribution near the
$\Gamma-Y$ paths [Fig. 2(b)]. This distribution is a bit different from the
joint density of states (see SM), suggesting anisotropic interband transition
contribution $\langle m|\bm{r}|n\rangle$ at different $\bm{k}$. We then fit
$\eta_{xy}^{x}(\omega_{0})$ with respect to $\lambda$ by linear regression
$\eta_{xy}^{x}(\omega_{0})=C_{1}\lambda+o(\lambda^{2})$ [Fig. 2(c)], which
gives $C_{1}=6.2\times
10^{3}\,\frac{\mu\mathrm{A}}{\mathrm{V}^{2}\mathrm{\AA}\cdot\mathrm{ps}}$.

Angular momentum current under LPL irradiation. Photocurrent generated by CPL
also includes LPL shift current feature, while LPL irradiation does not have
such problem Ahn  _et al._ (2020). Hence, next we consider LPL induced
injection current. For time-reversal symmetry systems, we have recently
demonstrated that LPL could generate AM photocurrents, showing injection
current feature Mu  _et al._ (2021); Xu  _et al._ (2021b). Photo-induced AM
currents include both spin and orbital currents. They indicate that electrons
carrying opposite spin and orbital moments are flowing oppositely. Note that
orbital angular momentum is a largely overlooked degree of freedom that has
attracted increasingly attention recently Bernevig  _et al._ (2005). We
evaluate the spin (and orbital) photocurrent conductivity rate Bhat  _et al._
(2005); Mu  _et al._ (2021); Fei  _et al._ ; Fei  _et al._ (2021)

$\begin{split}\eta_{bc}^{a;\hat{O}_{z}}(\omega)&=\frac{d}{dt}\sigma_{bc}^{a;\hat{O}_{z}}(0;\omega,-\omega)=-\frac{\pi
e^{3}}{2\hbar^{2}}\int_{BZ}\frac{d^{2}\bm{k}}{(2\pi)^{2}}\\\
&\times\mathrm{Re}\sum_{n,m}f_{nm}j_{nm}^{a;\hat{O}_{z}}\\{r_{mn}^{b},r_{nm}^{c}\\}\delta(\omega_{mn}-\omega)\end{split}$
(2)

Here $j_{nm}^{a;O_{z}}$ denotes the difference of angular momentum current for
band-$n$ and $m$, and
$j^{a;\hat{O}_{z}}=\\{v_{a},\hat{O}_{z}\\}/2e,\,(\hat{O}=\hat{S},\hat{L})$ Shi
_et al._ (2006); Cysne  _et al._ (2021). The $\\{r_{mn}^{b},r_{nm}^{c}\\}$
refers to anticommutation between $r_{mn}^{b}$ and $r_{nm}^{c}$. Here we focus
on the spin-$z$ and orbital-$z$ components, since they are uni-directionally
pointing along the out-of-plane directions in most $k$-points (singly
degenerate) Pan and Zhou (2020); Mu  _et al._ (2021). This is because $\bm{S}$
and $\bm{L}$ are pseudovectors and the $Pmn2_{1}$ contains glide mirror
operation $\\{\mathcal{M}_{z}|\frac{1}{2}\frac{1}{2}0\\}$.

Figure 3 shows our calculated nonlinear spin and orbital photocurrent rates
under $y$-LPL for different $\lambda$. The $x$-LPL responses are plotted in
SM. Similar as in the $\eta_{xy}^{x}$ case, the static $\mathcal{P}$-symmetric
system constraints that $\eta_{yy}^{x;S_{z}}(\lambda=0)=0$ and
$\eta_{yy}^{x;L_{z}}(\lambda=0)=0$. Positive and negative $\lambda$ gives
opposite values for the AM photocurrent rates. One has to note that under LPL,
there is no shift current along $x$, due to presence of $\mathcal{M}_{x}$.
Thus, there is no charge current mixture in the current case (pure spin and
orbital currents). We observe that the orbital photocurrent has much larger
values than the spin photocurrent, owing to small spin-orbit coupling in P.
This is in contrast with the common picture that orbital moments are usually
quenched in magnetic systems under strong and isotropic crystal field effect,
where magnetisms are mainly from spin degree of freedom. We also fit the these
photocurrent rates at $\omega_{0}=2.7\,\mathrm{eV}$ by odd power polynomial
function using $\lambda$,
$\eta_{yy}^{x;S_{z}}(\omega_{0})=C_{1}^{S}\lambda+C_{3}^{S}\lambda^{3}+o(\lambda^{4})$
and
$\eta_{yy}^{x;L_{z}}(\omega_{0})=C_{1}^{L}\lambda+C_{3}^{L}\lambda^{3}+o(\lambda^{5})$.
Our results are $C_{1}^{S}=-10.8$, $C_{3}^{S}=-2.3\times 10^{3}$,
$C_{1}^{L}=-4.6\times 10^{3}$, and $C_{3}^{L}=-6.7\times 10^{4}$. The units of
$C_{1}^{S,L}$ and $C_{3}^{S,L}$ are
$\frac{\mu\mathrm{A}}{\mathrm{V}^{2}\mathrm{\AA}\cdot\mathrm{ps}}\frac{\hbar}{2e}$
and
$\frac{\mu\mathrm{A}}{\mathrm{V}^{2}\mathrm{\AA}^{3}\cdot\mathrm{ps}}\frac{\hbar}{2e}$,
respectively. Both of them show strong nonlinear $\lambda^{3}$ relationship
due to the interplay of bandgap and wavefunction variations under $B_{1u}$.

Time-dependent injection charge, spin, and orbital currents. We now integrate
$\eta$ in time to explore these second order photoconductivities under
coherent $B_{1u}$ oscillation

$\frac{d\sigma_{\mathrm{inj}}(t)}{dt}=\eta(\omega_{0},\lambda(t))$ (3)

where $\lambda(t)=\lambda_{0}\sin(2\pi\omega_{\mathrm{IR}}t+\phi)$ is the
$B_{1u}$ mode vibration equation of motion with amplitude of
$\lambda_{0}=0.1\,\mathrm{\AA}$ and frequency
$\omega_{\mathrm{IR}}=4.2\,\mathrm{THz}$. $\phi\in(-\pi,\pi]$ is the vibration
phase at the initial state ($t=0$). The initial condition is
$\sigma_{\mathrm{inj}}(t=0)=0$.

As shown in Fig. 4(a)$-$(c), the photocurrent oscillates around an average
value, which is determined by $\phi$. In the coherent vibration case, the
$\phi$ is the same in different unit cells. When $\phi=\pm\pi/2$, one always
obtains an oscillating (charge, spin, and orbital) current. The time average
of such current is zero $-$ indicating a pure alternating current (with
propagating directions opposite for $\phi=\pm\pi/2$). When we perform Fourier
transformation of such currents [dashed curves in Fig. 4(d)$-$(f)], it shows
sharp peak at $\omega_{\mathrm{IR}}=4.2\,\mathrm{THz}$. For the spin and
orbital photocurrents, addtional peaks at
$3\omega_{\mathrm{IR}}=12.6\,\mathrm{THz}$ can be seen, consistent with the
$\lambda^{3}$-dependence of $\eta_{yy}^{x;O_{z}}(\omega_{0})$ in these cases.
Depending on the incident photon energy, such non-linearity may be tuned.
Remarkably, when $\phi=0$ and $\phi=\pi$, one can see that in additional to
the oscillatory nature of these currents ($\omega_{\mathrm{IR}}$ and
$3\omega_{\mathrm{IR}}$ for AM photocurrent peaks), they also contain large
static current (or AM currents) features with $0\,\mathrm{THz}$ contributions.
This indicates that one could observe a uni-directional injection
photocurrent, the magnitude of which wiggles but it always flows toward one
direction. Depending on the sign of $C$ (and $C^{S}$, $C^{L}$), the
photocurrent can be positive or negative. The $\phi=0$ and $\phi=\pi$ give
opposite propagating currents. Note that the time integral cannot be taken for
arbitrarily long time, as electron relaxation time $\tau$ is always finiteXu
_et al._ (2021b); Barik and Sau (2020) (but usually longer than a few
picoseconds).

Figure 4: Time evolution of (a) injection photoconductivity under CPL, (b)
spin and (c) orbital photoconductivity under LPL. The current propogation
strongly depends on phase $\phi$. In (d)$-$(f) we show frequency distributions
of (a)$-$(c) photocurrents by Fouier transformations. Incident photon energy
of $\omega_{0}=2.7\,\mathrm{eV}$ is used.

Transition metal dichalcogenide monolayer case. So far we have demonstrated
that coherent IR-active oscillation could induce finite injection current in a
static $\mathcal{P}$-symmetric system, $\alpha$-P monolayer. In order to show
that it is a universal effect, we now briefly discuss such oscillating
injection current in another well-studied monolayer, transition metal
dichalcogenide MoTe2 Wang  _et al._ (2017). This system does not possess
$\mathcal{P}$ at its static structure, but it includes a $\mathcal{M}_{x}$
vertical mirror plane [inset of Fig. 5(a)]. It symmetrically forbids CPL
injection current and LPL AM photocurrent along $y$. In the unit-cell coherent
vibrations, there is a doubly degenerate $E^{\prime}$ mode (at
$\omega_{\mathrm{IR}}^{\prime}=7.1\,\mathrm{THz}$) being IR-active. It
contains two vertical vibration eigenvectors, along $x$ and $y$ directions,
respectively. They can be individually excited using linearly polarized THz
pulses. Hence, we only focus on the $x$-mode which breaks $\mathcal{M}_{x}$
dynamically.

Figure 5: (a) IR-active mode (along $x$, zigzag direction) driven variation of
total energy and electric polarization ($P_{x}$) in MoTe2. Inset shows atomic
geometry and vibration mode. (b) shows the variation of injection photocurrent
rate under CPL (unit of
$\frac{\mu\mathrm{A}}{\mathrm{V}^{2}\cdot\mathrm{ps}}$) and injection AM
photocurrent rate under LPL (unit of
$\frac{\mu\mathrm{A}}{\mathrm{V}^{2}\cdot\mathrm{ps}}\frac{\hbar}{2e}$). All
currents propagate along $y$ (armchair).

We again consider vibrations with $\lambda$ between $-0.1\,\mathrm{\AA}$ and
$0.1\,\mathrm{\AA}$. The total energy increases quadratically up to 43
meV/f.u. The system shows a finite electric polarization along zigzag
direction, $P_{x}$, which linearly increases with $\lambda$. Then we consider
the injection photocurrents. When the ionic system is coherently excited
($E_{x}^{\prime}$), the $\mathcal{M}_{x}$ symmetry constraints lift. From Fig.
5 we plot their values at a chosen incident photon energy
($1.7\,\mathrm{eV}$). Nonlinear $\lambda$-dependence can be clearly observed
in this case, especially for orbital and spin photocurrents. If we integrate
the photocurrent rate as in Eq. 3, we obtain similar uni-directional injection
photocurrent ($0\,\mathrm{THz}$ static current component), and the nonlinear
contributions ($3\omega_{\mathrm{IR}}^{\prime}$ and
$5\omega_{\mathrm{IR}}^{\prime}$ contributions) can be obtained (see SM).

##  .3 Discussion

Before concluding, we would like to note that the coherent IR-active vibration
at $\Gamma$ is different from other ionic vibrations or motions, such as
thermal excitations, Raman active modes, or BZ boundary modes. Due to equi-
partition theorem, the thermally excited phonons are random in phase and do
not possess a definite wavelength. Thus, there is no collective and coherent
motion of ions under thermal excitation. The phase factor $\phi$ does not have
a long range correlation in different unit cells. This will result in
qualitatively similar behaviors as the static configuration, and both
alternative (at $\omega_{\mathrm{IR}}$ and static injection currents vanish.
Collective $\Gamma$-centered Raman modes and BZ boundary modes usually do not
break the spatial symmetry $\mathcal{P}$ or $\mathcal{M}$. Hence, even though
they would affect the electronic structure (band dispersion, etc) and can
change the otherwise existing current values, the spatial symmetry arguments
do not break. Then the oscillating effect will be different and may not be as
strong as that in our current case. We will study their effects in a separate
work. We also want to note that recently Schlipf and Giustino proposed dynamic
Rashba-Dresselhaus (dRD) effect under coherent IR-active phonon oscillations
Schlipf and Giustino , where they showed that the dRD can be observed at a
given time $t$. But in their case the dRD would diminish if time average is
performed. In our proposal, on the contrary, we can obtain static (zero
frequency) photocurrent component, which does not require an ultrafast
dynamics observation setup. Of course, the recently developed ultrafast
dynamic approaches could provide high resolution observation of the injection
photocurrent Luo  _et al._ (2021), and would be useful for experimental
verification and applications of our prediction.

In our current cases, the $\alpha$-P and MoTe2 monolayers possess out-of-plane
glide and mirror symmetries, respectively. They forbid in-plane angular
momentum components, since $\bm{L}$ and $\bm{S}$ are both axial vectors.
However, one has to note that $L_{z}$ does not commute with the Hamiltonian as
it is determined by crystal field effect. Hence, the orbital angular momentum
is generally not a good quantum number. This would yield additional orbital
torque effectGo  _et al._ (2020) in the continuum equation of motion.
Actually, similar problem occurs in quantum spin Hall insulator when Rashba
spin-orbit coupling is present. Various approaches have been suggested to
avoid such problem, for example, one can project $L_{z}$ onto eigenstate space
via
$\mathbb{L}_{z}=\sum_{n}|\psi_{n\bm{k}}\rangle\langle\psi_{n\bm{k}}|L_{z}|\psi_{n\bm{k}}\rangle\langle\psi_{n\bm{k}}|$Cysne
_et al._ (2021); Phong  _et al._ (2019); Go  _et al._ (2018), which
corresponds to the diagonal terms in density matrix. Then the
$[\mathbb{L}_{z},\hat{H}]=0$. Our test calculations show that such a scheme
gives qualitatively the same results as discussed here with the symmetry
arguments still hold. But a well-defined orbital operator still requires
theoretical efforts.

The definition of angular momentum current is another open problem. When the
out-of-plane spin-$z$ component is conserved, one may define the spin current
as $j_{\mathrm{sp}}=j_{\uparrow}-j_{\downarrow}$, where $\uparrow$ and
$\downarrow$ indicate spin-up and spin-down channels, respectively. However,
in a general case, there are several spin current expressions. A widely-
adopted approach (also in the current study) is
$j_{\mathrm{sp}}^{a;S_{z}}=\\{v_{a},S_{z}\\}$, which has been used to evaluate
spin Hall effect and spin Seebeck effect, etc. On the other hand, it has also
been proposed that a proper definition takes the form of
$j_{\mathrm{sp}}^{a;S_{z}}=\frac{d(r_{a}S_{z})}{dt}$. We have compared both
these approaches and showed that they yield similar results with small
differencesXu  _et al._ (2021b). In the current study, our focus is on the
prediction of static injection current components induced by coherent phonon
vibration, while such differences would not affect the main results.

In this study, we use the quadratic recursive response theory to evaluate the
injection current rate, which essentially corresponds to
$\mathrm{Tr}(j^{(0)}\rho^{(2)})$, where $j^{(0)}$ is electric field
$E$-independent current operator and $\rho^{(2)}$ is the second order
perturbative density matrix operator. They can be calculated according to Kubo
theory. Note that there are other contributions, such as self energy and
vertex correctionsAvdoshkin  _et al._ (2020), which cannot be simply included
in the current density functional theory calculations and is not computed
here. In addition, the $\mathrm{Tr}(j^{(1)}\rho^{(1)})$ contribution, where
$j^{(1)}$ originates from anomalous velocity and is linearly dependent on
electric field $E$, would vanish in such topologically trivial and insulating
systems. The band edge dephasing time is calculated to be a few femtoseconds,
increases with temperature reduces (see SM). This is consistent with the
previous work Zereshki  _et al._ (2018), and is sufficiently long for
nonlinear optical process (absorbing and releasing photons during the two-
photon-one-electron process).

As discussed previously, bulk photovoltaic effect includes shift current and
injection current, which belong to different microscopic mechanismsTan  _et
al._ (2016). Taking charge current as an example (similar arguments hold for
AM current), the shift current arises from the mismatch of wavefunction center
position in the valence and conduction bands. On the other hand, injection
current originates from the mismatch of velocity operator in the valence and
conduction bands. Hence, the injection current rate is more fundamental. Most
previous studies evaluate the injection current for static configurations via
$\sigma_{\mathrm{inj}}=\eta\tau$ where $\tau$ is electron relaxation lifetime,
determined by disorder, temperature, impurities, and electron-phonon
scattering, etcWang and Qian (2019); Xu  _et al._ (2021b); Fei  _et al._ . In
this study we show that coherent vibration gives a time dependent rate
$\eta(t)$, hence the ultimate injection current can be dynamic, nonzero, and
even contains static components. The shift current, on the contrary, does not
have such property. In realistic experimental setups, above bandgap photon
absorption produces hot carriers (electrons and holes), and their drifts
contribute to photoconductivities as well (e.g., ballistic current). This is
not the situation discussed in the current work. One also has to note that
damping effect would gradually suppress the static current component (see SM).

##  .4 Conclusion

In summary, we apply first-principles calculations to show that
$\Gamma$-centered coherent IR-active phonon oscillation that breaks spatial
symmetry could trigger nonlinear injection photocurrent in otherwise silent
systems. This photocurrent could be pure alternating or uni-directionally
propagating, depending on the phase of the coherent phonon. Our work shows
that unlike the widely studied static or thermally equilibrium conditions, the
synergistic interplay between the coherent phonon, electron, and photon could
lead to intriguing and exotic nonlinear optical responses in ultrafast
dynamics. The system setup we propose here may serve as a future THz
electromagnetic wave detector for the next generation communication.

###### Acknowledgements.

This work is supported by the National Natural Science Foundation of China
(NSFC) under Grant Nos. 11974270, 21903063, and 11904353. J.Z. acknowledges
valuable discussions with H. Xu (MIT) and Y. Gao (USTC), and the computational
resources provided by HPCC platform of Xi’an Jiaotong University.

## References

  * Fridkin (2001) V.M. Fridkin, “Bulk photovoltaic effect in noncentrosymmetric crystals,” Crystallogr. Rep. 46, 654–658 (2001).
  * Aversa and Sipe (1995) C. Aversa and J. E. Sipe, “Nonlinear optical susceptibilities of semiconductors: Results with a length-gauge analysis,” Phys. Rev. B 52, 14636–14645 (1995).
  * Gong  _et al._ (2018) S.-J. Gong, F. Zheng, and A. M. Rappe, “Phonon influence on bulk photovoltaic effect in the ferroelectric semiconductor gete,” Phys. Rev. Lett. 121, 017402 (2018).
  * Tan  _et al._ (2016) L. Z. Tan, F. Zheng, S. M. Young, F. Wang, S. Liu, and A. M. Rappe, “Shift current bulk photovoltaic effect in polar materials—hybrid and oxide perovskites and beyond,” npj Comput. Mater. 2, 16026 (2016).
  * Hohenberg and Kohn (1964) P. Hohenberg and W. Kohn, “Inhomogeneous electron gas,” Phys. Rev. 136, B864–B871 (1964).
  * Kohn and Sham (1965) W. Kohn and L. J. Sham, “Self-consistent equations including exchange and correlation effects,” Phys. Rev. 140, A1133–A1138 (1965).
  * Kresse and Furthmüller (1996) G. Kresse and J. Furthmüller, “Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set,” Phys. Rev. B 54, 11169–11186 (1996).
  * Perdew  _et al._ (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, “Generalized gradient approximation made simple,” Phys. Rev. Lett. 77, 3865–3868 (1996).
  * Blöchl (1994) P. E. Blöchl, Phys. Rev. B 50, 17953–17979 (1994).
  * Liu  _et al._ (2014) H. Liu, A. T. Neal, Z. Zhu, Z. Luo, X. Xu, D. Tománek, and P. D. Ye, “Phosphorene: An unexplored 2d semiconductor with a high hole mobility,” ACS Nano 8, 4033–4041 (2014).
  * Zhou  _et al._ (2018) J. Zhou, H. Xu, Y. Li, R. Jaramillo, and J. Li, “Opto-mechanics driven fast martensitic transition in two-dimensional materials,” Nano Lett. 18, 7794–7800 (2018).
  * Guo  _et al._ (2021) H. Guo, W. Chu, O. V. Prezhdo, Q. Zheng, and J. Zhao, “Strong modulation of band gap, carrier mobility and lifetime in two-dimensional black phosphorene through acoustic phonon excitation,” J. Phys. Chem. Lett. 12, 3960–3967 (2021).
  * Xu  _et al._ (2021a) H. Xu, J. Zhou, H. Wang, and J. Li, “Light-induced static magnetization: Nonlinear edelstein effect,” Phys. Rev. B 103, 205417 (2021a).
  * Zhao  _et al._ (2006) H. Zhao, E. J. Loren, H. M. van Driel, and A. L. Smirl, “Coherence control of hall charge and spin currents,” Phys. Rev. Lett. 96, 246601 (2006).
  * Wang  _et al._ (2019) C. Wang, B. Lian, X. Guo, J. Mao, Z. Zhang, D. Zhang, B.-L. Gu, Y. Xu, and W. Duan, “Type-ii ising superconductivity in two-dimensional materials with spin-orbit coupling,” Phys. Rev. Lett. 123, 126402 (2019).
  * Ji  _et al._ (2019) Z. Ji, G. Liu, Z. Addison, W. Liu, P. Yu, H. Gao, Z. Liu, A.M. Rappe, C.L. Kane, E.J. Mele, and R Agarwal, “Spatially dispersive circular photogalvanic effect in a weyl semimetal,” Nat. Mater. 18, 955 (2019).
  * Sipe and Shkrebtii (2000) J. E. Sipe and A. I. Shkrebtii, “Second-order optical response in semiconductors,” Phys. Rev. B 61, 5337–5352 (2000).
  * Fregoso (2019) B. M. Fregoso, “Bulk photovoltaic effects in the presence of a static electric field,” Phys. Rev. B 100, 064301 (2019).
  * Wang and Qian (2019) H. Wang and X. Qian, Sci. Adv. 5, eaav9743 (2019).
  * Holder  _et al._ (2020) T. Holder, D. Kaplan, and B. Yan, “Consequences of time-reversal-symmetry breaking in the light-matter interaction: Berry curvature, quantum metric, and diabatic motion,” Phys. Rev. Research 2, 033100 (2020).
  * Taghizadeh  _et al._ (2017) A. Taghizadeh, F. Hipolito, and T. G. Pedersen, “Linear and nonlinear optical response of crystals using length and velocity gauges: Effect of basis truncation,” Phys. Rev. B 96, 195413 (2017).
  * Laturia  _et al._ (2018) A. Laturia, M.L. Van de Put, and M.G. Vandenberghe, npj 2D Mater. Appl. 2, 6 (2018).
  * Jiang  _et al._ (2017) Z. Jiang, Z. Liu, Y. Li, and W. Duan, “Scaling universality between band gap and exciton binding energy of two-dimensional semiconductors,” Phys. Rev. Lett. 118, 266401 (2017).
  * Mu  _et al._ (2021) X. Mu, Y. Pan, and J. Zhou, “Pure bulk orbital and spin photocurrent in two-dimensional ferroelectric materials,” npj Comput. Mater. 7, 61 (2021).
  * Ahn  _et al._ (2020) J. Ahn, G.-Y. Guo, and N. Nagaosa, “Low-frequency divergence and quantum geometry of the bulk photovoltaic effect in topological semimetals,” Phys. Rev. X 10, 041041 (2020).
  * Xu  _et al._ (2021b) H. Xu, H. Wang, J. Zhou, and J. Li, “Pure spin photocurrent in non-centrosymmetric crystals: bulk spin photovoltaic effect,” Nat. Commun. 12, 4330 (2021b).
  * Bernevig  _et al._ (2005) B. A. Bernevig, T. L. Hughes, and S.-C. Zhang, “Orbitronics: The intrinsic orbital current in $p$-doped silicon,” Phys. Rev. Lett. 95, 066601 (2005).
  * Bhat  _et al._ (2005) R. D. R Bhat, F. Nastos, A. Najmaie, and J. E. Sipe, “Pure spin current from one-photon absorption of linearly polarized light in noncentrosymmetric semiconductors,” Phys. Rev. Lett. 94, 096603 (2005).
  * (29) R. Fei, X. Lu, and L. Yang, “Intrinsic spin photogalvanic effect in nonmagnetic insulator,” arXiv , 2006.10690.
  * Fei  _et al._ (2021) R. Fei, S. Yu, Y. Lu, L. Zhu, and L. Yang, “Switchable enhanced spin photocurrent in rashba and cubic dresselhaus ferroelectric semiconductors,” Nano Lett. 21, 2265–2271 (2021).
  * Shi  _et al._ (2006) J. Shi, P. Zhang, D. Xiao, and Q. Niu, “Proper definition of spin current in spin-orbit coupled systems,” Phys. Rev. Lett. 96, 076604 (2006).
  * Cysne  _et al._ (2021) T. P. Cysne, M. Costa, L. M. Canonico, M. B. Nardelli, R. B. Muniz, and T. G. Rappoport, “Disentangling orbital and valley hall effects in bilayers of transition metal dichalcogenides,” Phys. Rev. Lett. 126, 056601 (2021).
  * Pan and Zhou (2020) Y. Pan and J. Zhou, “Toggling valley-spin locking and nonlinear optical properties of single-element multiferroic monolayers via light,” Phys. Rev. Appl. 14, 014024 (2020).
  * Barik and Sau (2020) Tamoghna Barik and Jay D. Sau, “Nonequilibrium nature of nonlinear optical response: Application to the bulk photovoltaic effect,” Phys. Rev. B 101, 045201 (2020).
  * Wang  _et al._ (2017) C. Wang, X. Liu, L. Kang, B.-L. Gu, Y. Xu, and W. Duan, “First-principles calculation of nonlinear optical responses by wannier interpolation,” Phys. Rev. B 96, 115147 (2017).
  * (36) M. Schlipf and F. Giustino, “Dynamic rashba-dresselhaus effect,” arXiv , 2004.10477.
  * Luo  _et al._ (2021) L Luo, D. Cheng, B. Song, L.-L. Wang, C. Vaswani, P. M. Lozano, G. Gu, C. Huang, R. H. J. Kim, Z. Liu, J.-M. Park, Y. Yao, K. Ho, I. E. Perakis, Q. Li, and J. Wang, Nat. Mater. 20, 329–334 (2021).
  * Go  _et al._ (2020) Dongwook Go, Frank Freimuth, Jan-Philipp Hanke, Fei Xue, Olena Gomonay, Kyung-Jin Lee, Stefan Blügel, Paul M. Haney, Hyun-Woo Lee, and Yuriy Mokrousov, “Theory of current-induced angular momentum transfer dynamics in spin-orbit coupled systems,” Phys. Rev. Research 2, 033401 (2020).
  * Phong  _et al._ (2019) Võ Tien Phong, Zachariah Addison, Seongjin Ahn, Hongki Min, Ritesh Agarwal, and E. J. Mele, “Optically controlled orbitronics on a triangular lattice,” Phys. Rev. Lett. 123, 236403 (2019).
  * Go  _et al._ (2018) Dongwook Go, Daegeun Jo, Changyoung Kim, and Hyun-Woo Lee, “Intrinsic spin and orbital hall effects from orbital texture,” Phys. Rev. Lett. 121, 086602 (2018).
  * Avdoshkin  _et al._ (2020) Alexander Avdoshkin, Vladyslav Kozii, and Joel E. Moore, “Interactions remove the quantization of the chiral photocurrent at weyl points,” Phys. Rev. Lett. 124, 196603 (2020).
  * Zereshki  _et al._ (2018) P. Zereshki, Y. Wei, F. Ceballos, M. Z. Bellus, S. D. Lane, S. Pan, R. Long, and H. Zhao, “Photocarrier dynamics in monolayer phosphorene and bulk black phosphorus,” Nanoscale 10, 11307–11313 (2018).


# Safe Reinforcement Learning Using Robust Action Governor

Yutong Li<EMAIL_ADDRESS>
University of Michigan  Ann Arbor  MI  USA  Nan Li<EMAIL_ADDRESS>
University of Michigan  Ann Arbor  MI  USA  H. Eric Tseng<EMAIL_ADDRESS>
Ford Motor Company  Dearborn  MI  USA  Anouck Girard<EMAIL_ADDRESS>
University of Michigan  Ann Arbor  MI  USA  Dimitar Filev dfilev@@ford.com
Ford Motor Company  Dearborn  MI  USA  Ilya Kolmanovsky<EMAIL_ADDRESS>
University of Michigan  Ann Arbor  MI  USA

###### Abstract

Reinforcement Learning (RL) is essentially a trial-and-error learning
procedure which may cause unsafe behavior during the exploration-and-
exploitation process. This hinders the application of RL to real-world control
problems, especially to those for safety-critical systems. In this paper, we
introduce a framework for safe RL that is based on integration of a RL
algorithm with an add-on safety supervision module, called the Robust Action
Governor (RAG), which exploits set-theoretic techniques and online
optimization to manage safety-related requirements during learning. We
illustrate this proposed safe RL framework through an application to
automotive adaptive cruise control.

###### keywords:

safety-critical systems, reinforcement learning, action governor, automotive
applications

##  1 Introduction

In Reinforcement Learning (RL), the agent interacts with the environment by
perceiving environment’s states and selecting the action that maximizes the
long-term return based on a real-valued reward signal (Sutton and Barto,
2018). The success of RL has been apparent in a broad spectrum of applications
(Mnih et al., 2013; Kober et al., 2013; Abbeel et al., 2007). However, RL is
essentially a trial-and-error learning process which may cause unsafe behavior
during the learning process. This hinders the RL real-world applications,
especially to safety-critical systems.

One approach to addressing such safety issues is the safe RL approach (Garcıa
and Fernández, 2015). In particular, risk-sensitive safe RL aims to promote
the constraint satisfaction via balancing the long-term return and the risk of
reaching the unsafe region (Geibel and Wysotzki, 2005). Methods based on
policy optimization with constraints were also proposed, where a constraint on
the probability of the system being able to return to the safe region is
enforced (Wachi et al., 2018; Moldovan and Abbeel, 2012). However, safety
cannot be guaranteed via these model-free approaches, as the agent needs to
learn to operate safely via its interactions with the environment, which may
lead to constraint violations.

To guarantee the safety constraint satisfaction during the learning process,
an effective approach is to incorporate the system model information. In
general, the safe region can be determined by the system model, and a control
policy that keeps system state staying inside this safe region can then be
computed (Aswani et al., 2013; Fisac et al., 2018; Larsen et al., 2017; Sloth
et al., 2012). The main advantage of these model-based methods is that within
the interior of the safe region, the RL agent can explore safely to improve
the performance. However, the optimality of the trained policy is highly
dependent on the size of the safe region computed based on the system model.

Following the idea of this model-based direction, in this paper, we propose a
novel safe RL framework that also exploits the system model to design a safe
set that regulates the RL explorations to guarantee system safety. In
particular, the proposed framework exploits an add-on module, called the
Robust Action Governor (RAG), to manage safety. With the RAG, an arbitrary RL
algorithm can be integrated into the framework and lead to safe RL. The RAG
enforces safety constraints by monitoring, and minimally modifying when
necessary, the control signal produced by the nominal RL policy to a
constraint-admissible one.

The proposed safe RL framework based on RAG operates on the basis of set-
theoretic techniques and online optimization. Similar approaches were proposed
by integrating the Reference Governor (RG) with RL to enforce constraint
satisfaction (Li et al., 2018, 2019). A distinguishing feature of our safe RL
framework based on RAG is that, unlike the RG which modifies the reference
input to the controller, the RAG modifies controller output signal, i.e. the
RAG can be placed closer to the environment/plant. A direct consequence is
that our safe RL framework based on RAG can be used to train lower-level
controls, while the approach based on RG in Li et al. (2018, 2019) can only be
used to train higher-level planners where the plant must have already been
stabilized by some controller that is fixed. Another advantage of RAG compared
to RG is that the use of RAG yields a larger safe set, and thereby, can
potentially achieve better control performance (Li et al., 2021a). On the
other hand, unlike safe RL approaches based on the use of control barrier
functions (Cheng et al., 2019; Sloth et al., 2012), our RAG is formulated
based on discrete-time models from the start and thereby can be more directly
applied in a digital setting (e.g., one does not need to modify the algorithm
to account for sampling time and discrete updates).

In summary, the contributions of this paper include: 1) establishing a safe RL
framework based on RAG, 2) extending the theoretical results and computational
methods in Li et al. (2021a) for the Action Governor for discrete-time linear
systems to their robust versions for linear systems with additive bounded
disturbances and non-convex constraints, and 3) illustrating effectiveness of
the proposed safe RL framework using an example relevant to automated driving
and adaptive cruise control.

##  2 Conventional RL and Safe RL

Conventional RL optimizes the agent’s action via trial-and-error to maximize
its accumulated reward (or minimize accumulated cost) through continuously
interacting with the environment, as illustrated in Figure 1a. Specifically,
at each time instant $k\in\mathbb{Z}_{\geq 0}$, the agent takes a measurement
of the state $x(k)\in\mathbb{R}^{n}$, executes a control
$u_{\phi}(k)\in\mathbb{R}^{m}$ and collects a reward $r(k)\in\mathbb{R}$. A
control policy, $\pi:\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}$, which is a
mapping from the state space $\mathbb{R}^{n}$ to the action space
$\mathbb{R}^{m}$, describes the agent’s behavior. The control policy is
learned from the experience $\\{x(k),u_{\phi}(k),r(k)\\}$ to maximize the
long-term reward $R=\sum_{k=0}^{\infty}\gamma^{k}r(k)$, where $\gamma\in(0,1)$
is a discount factor. However, to maximize the long-term reward, RL agent
needs to explore within the action space, which may cause unsafe behaviors
(i.e., violation of certain safety constraints) during the learning process.
This feature hinders the application of RL to actual engineering systems, and
motivates us to propose a safe RL framework enabling the agent to learn a
policy safely.

The proposed safe RL framework is illustrated in Figure 1b. An add-on module,
termed Robust Action Governor (RAG), is introduced between the RL agent and
the environment (which represents the system that the control policy is acting
on, e.g., the plant in a conventional control setting). The RAG monitors the
control signal $u_{\phi}(k)$ generated by the RL agent and corrects the ones
that may cause unsafe behavior. We will describe the properties and design
procedure of the RAG in the subsequent section.

Figure 1: Schematic diagram of conventional RL (a) and safe RL (b).

##  3 Robust Action Governor

###  3.1 Problem Formulation

In this section, we review the main results in Li et al. (2021a) on the Action
Governor and extend them from disturbance-free systems to systems subject to
bounded disturbances, enabling the Action Governor to be applicable to more
general cases. Consider a discrete-time linear model with additive
disturbances as follows:

$x(k+1)=Ax(k)+Bu(k)+Ew(k),$ (1)

where $x(k)\in\mathbb{R}^{n}$ is the state at discrete time instant
$k\in\mathbb{Z}_{\geq 0}$ and $u(k)\in\mathbb{R}^{m}$ is the control input.
The vector $w(k)$ represents an unmeasured disturbance input and we assume
that it is bounded in a known set, i.e., $w(k)\in\mathcal{W}=\\{w:Mw\leq
m\\}$, where $\mathcal{W}$ is a known polytope in $\mathbb{R}^{n_{w}}$. In
applications to RL, $w(k)$ can account for model mismatch, e.g., differences
between the actual lead vehicle acceleration and assumed zero lead vehicle
acceleration in an adaptive cruise control (ACC) setting. We further assume
that state, input and disturbance constraints are not coupled, i.e.,
$\mathcal{T}=\\{(x,u,w)|x\in\mathcal{X},u\in\mathcal{U},w\in\mathcal{W}\\}=\mathcal{X}\times\mathcal{U}\times\mathcal{W}$,
where $\mathcal{X}$ represents the feasible set of states and $\mathcal{U}$ is
the control input set. We assume that a nominal control policy $\phi$ has been
designed for the system (1),

$u_{\phi}=\phi(x(k),x_{r}(k),w(k),k),$ (2)

where $x_{r}(k)\in\mathbb{R}^{r}$ denotes a reference signal for the system
(1), which defines the control objective. Note that there is no further
assumption on the nominal control $\phi$, and it can take any form, e.g.
nonlinear and time-varying. In particular, in this paper, $\phi$ represents
the RL policy.

We assume that the system is subject to a safety requirement of the form

$x(k)\in\mathcal{X}=\mathbb{R}^{n}\setminus X_{0},\ \ \ \forall
k\in\mathbb{Z}_{\geq 0},$ (3)

where the unsafe set $X_{0}$ can be expressed as a finite union of polytopes,
i.e.,

$X_{0}=\bigcup\limits_{i=1}^{n_{g}}\\{x\in\mathbb{R}^{n}:G_{i}x<g_{i}\\},$ (4)

where $G_{i}\in\mathbb{R}^{n_{g}\times n}$ and $g_{i}\in\mathbb{R}^{n_{g}}$.
Note that the expressions (3) and (4) can represent a broad range of safety
requirements, including box constraints, obstacle avoidance constraints, etc.
In general, the feasible set $\mathcal{X}$ can be non-convex.

As the safety requirement (3) may not be strictly handled by the nominal
control policy (2) (especially when (2) represents an RL policy), this
motivates us to design an add-on scheme to enforce (3). In particular, we
exploit a supervisory solution, the RAG, to satisfy this requirement, as
illustrated in Figure 1b. The RAG monitors the nominal control input
$u_{\phi}$ and, if necessary, minimally modifies $u_{\phi}$ to guarantee that
the system state can stay outside $X_{0}$ for the present and all future time
instants even in the presence of disturbances.

In particular, at each time instant, the RAG solves the following constrained
optimization problem to enforce the safety requirement (3),

$\displaystyle u(k)=$
$\displaystyle\operatorname*{arg\,min}_{u\in\mathcal{U}}\
{\|u-u_{\phi}(k)\|}_{S}^{2}$ (5a) $\displaystyle\text{ subject to
}\,Ax(k)+Bu+Ew\in X_{\text{safe}},\forall w\in\mathcal{W}$ (5b)

where $X_{\text{safe}}\subset\mathcal{X}$ is a “safe set” which will be
introduced in the next section. The function
$\|\cdot\|_{S}=\sqrt{(\cdot)^{T}S(\cdot)}$ is used to penalize the difference
between the nominal control $u_{\phi}$ and the modified control $u$, where the
matrix $S\in\mathbb{R}^{m\times m}$ is positive-definite.

###  3.2 Safe Set and Unrecoverable Sets

To enforce both present and future safety, the safe set $X_{\text{safe}}$ is
characterized by the following requirements: There exists a state-feedback
control $u(x)$ such that

  * •

$u(x)\in\mathcal{U}$ for all $x\in X_{\text{safe}}$;

  * •

Given $x(0)\in X_{\text{safe}}$, this control can keep all future states
$\\{x(1),x(2),...\\}$ within $\mathcal{X}$ in the presence of disturbances
$\\{w(0),w(1),...\\}\subset\mathcal{W}$.

The derivation of $X_{\text{safe}}$ relies on its complementary sets, termed
unrecoverable sets, which are recursively defined as follows,

$\displaystyle\begin{split}X_{k}&=X_{0}\cup\\{x\in\mathbb{R}^{n}:\forall
u\in\mathcal{U},\exists w\in\mathcal{W},\text{s.t.}\ Ax+Bu+Ew\in X_{j}\\\
&\quad\text{for some}\,j=0,...,k-1\\}\\\
&=X_{0}\cup\left\\{x\in\mathbb{R}^{n}:Ax\in\bigcup\limits_{j=0}^{k-1}X_{j}\oplus((-E)\circ\mathcal{W})\sim(B\circ\mathcal{U})\right\\},\end{split}$
(6)

where $\oplus$, $\sim$, and $\circ$ represent the Minkowski sum, Pontryagin
difference, and affine mapping operations of sets. Note that in the derivation
of (6) we have used the fact that $S=\\{v:\forall u\in U,\exists z\in
Z,\text{s.t.}\ v=z-u\\}=\\{v:\forall u\in U,\exists z\in Z,\text{s.t.}\
z=v+u\\}=\\{v:v+u\in Z,\forall u\in U\\}=Z\sim U$.

The unrecoverable set $X_{k}$ in (6) has the following properties elaborated
in Propositions 1 to 3.

Proposition 1: 1) If $x_{0}\in X_{k}$, then for any state-feedback control
sequence $\\{u_{0}(x_{0})...,u_{k-1}(x_{k-1})\\\
\\}\in\mathcal{U}\times...\times\mathcal{U}$, there exists a disturbance
sequence $\\{w_{0},...,w_{k-1}\\}\in\mathcal{W}\times...\times\mathcal{W}$
such that $x_{j}\in X_{0}$ for some $0\leq j\leq k$ .
2) Let $x_{0}$ be given. If for any state-feedback control sequence
$\\{u_{0}(x_{0}),...,u_{k-1}(x_{k-1})\\}\in\mathcal{U}\times...\times\mathcal{U}$,
there exists a sequence
$\\{w_{0},...,w_{k-1}\\}\in\mathcal{W}\times...\times\mathcal{W}$ such that
$x_{j}\in X_{0}$ for some $0\leq j\leq k$, then $x_{0}\in X_{k}$.

The proof of Proposition 1 can be constructed following similar steps as those
in the proofs of Propositions 1 and 2 in Li et al. (2021a). It is omitted here
but will be included in our subsequent publication Li et al. (2021b).
Proposition 1 implies that for any state-feedback control sequence
$\\{u_{0}(x),...,u_{k-1}(x)\\}\in\mathcal{U}\times...\times\mathcal{U}$, there
exists a disturbance sequence
$\\{w_{0},...,w_{k-1}\\}\in\mathcal{W}\times...\times\mathcal{W}$ that causes
the state trajectory to enter $X_{0}$ during the steps $0,...,k$ if and only
if $x_{0}\in X_{k}$. Or equivalently, for any admissible disturbances
$\\{w_{0},...,w_{k-1}\\}\in\mathcal{W}\times...\times\mathcal{W}$, there
exists a state-feedback control sequence
$\\{u_{0}(x),...,u_{k-1}(x)\\}\in\mathcal{U}\times...\times\mathcal{U}$ that
can prevent the state trajectory from entering $X_{0}$ over steps $0,...k$ if
and only if $x_{0}\in\mathcal{X}\setminus X_{k}$.

Proposition 2: For each $k=0,1,...,$ we have $X_{k}\subset X_{k+1}$, i.e.
$X_{k}$ is an increasing sequence of sets. In turn,
$X_{\infty}=\lim_{k\to\infty}X_{k}$ exists and satisfies $X_{k}\subset
X_{\infty}$ for all $k$.

On the basis of the unrecoverable sets $X_{k}$, we define the safe set as
$X_{\text{safe}}=\mathcal{X}\setminus
X_{\infty}=\lim_{k\to\infty}(\mathcal{X}\setminus X_{k})$.

Proposition 3: For any $x\in X_{\text{safe}}$, it holds that (1)
$x\in\mathcal{X}=\mathbb{R}^{n}\setminus X_{0}$ and (2) there exists
$u\in\mathcal{U}$ such that $Ax+Bu+Ew\in X_{\text{safe}}$ for all
$w\in\mathcal{W}$.

The proofs of the above Propositions 2 and 3 can also be constructed following
proofs of the Propositions 3 and 4 in Li et al. (2021a) and will be included
in Li et al. (2021b). Note that one needs to use the formula
$\bigcup\limits_{k=0}^{k}X_{k}\oplus W=\bigcup\limits_{k=0}^{k}(X_{k}\oplus
W)$ to complete these proofs, which can be easily shown to hold using the
definitions of Minkowski sum and union of sets.

Proposition 3 ensures that if the RAG operates based on (5), then a feasible
solution exists to (5) for all $k$, and the safety requirement (3) is
satisfied for all $k$ in spite of disturbance $w\in\mathcal{W}$. Note that the
exact determination of $X_{\text{safe}}$ relies on the set $X_{k}$ iteratively
computed according to (6) with $k\rightarrow\infty$. In practice, we
approximate $X_{\text{safe}}$ by $X_{\text{safe},k}=\mathcal{X}\setminus
X_{k}$ with $k$ being sufficiently large. We also note that according to
Proposition 1, the safe set $X_{\text{safe}}$ used by RAG approach is maximal.

###  3.3 Offline and Online Computations

Proposition 4: Suppose $A$ is invertible. Then, for each $k=1,2,...,$ we have
(1) $X_{k}$ can be represented as the union of a finite number of polytopic
sets, i.e., $X_{k}=\bigcup\limits_{j=1}^{r_{k}}X_{k,j}$ where $X_{k,j}$ is a
polytopic set for each $j=1,...,r_{k}$; and (2) $X_{k}$ can be numerically
computed using Algorithm 1.

The proof of Proposition 4 is similar to that of Proposition 6 in Li et al.
(2021a). RAG operates by solving the optimization problem in (5) at each time
step. Using the fact that
$\mathcal{X}\backslash(X_{k}\oplus(-E\circ\mathcal{W}))=\mathcal{X}\backslash
X_{k}\sim E\circ\mathcal{W}$, and
$X_{k}\oplus(-E\circ\mathcal{W})=\bigcup\limits_{j=1}^{r_{k}}\bigcap\limits_{i=1}^{s_{j}}\\{x\in\mathcal{X}:G_{i,j}x<g_{i,j}\\}$,
(5) can be transformed into a Mixed-Integer Quadratic Programming (MIQP)
problem with the constraints as following

$\displaystyle G_{i,j}(Ax(k)+Bu)\geq g_{i,j}-M(1-\delta_{i,j}),$ (7a)
$\displaystyle\delta_{i,j}\in\\{0,1\\},\forall i=1,...,s_{j},\forall
j=1,...,r_{k^{{}^{\prime}}},$ (7b)
$\displaystyle\sum\limits_{i=1}^{s_{j}}\delta_{i,j}=1,\forall
j=1,...,r_{k^{{}^{\prime}}},$ (7c)

where $M>0$ is a sufficiently large number.

Algorithm 1 Offline computation of $X_{k}$ for systems in (1)

Input: $A,B,E,X_{0},X_{k-1},\mathcal{U},\mathcal{W}$
Output:$\ X_{k}$

1:$\mathcal{H}\leftarrow\text{convhull}(X_{k-1}\oplus((-E)\circ\mathcal{W}))$

2:$\mathcal{D}\leftarrow\mathcal{H}\sim((B)\circ\mathcal{U}))$

3:$\mathcal{E}\leftarrow\mathcal{H}\setminus(X_{k-1}\oplus((-E)\circ\mathcal{W})))$

4:$\mathcal{F}\leftarrow\mathcal{E}\oplus((-B)\circ\mathcal{U})$

5:$\mathcal{G}\leftarrow\mathcal{D}\setminus\mathcal{F}$

6:$X_{k}\leftarrow X_{0}\bigcup A^{-1}\mathcal{G}$

##  4 Safe RL With RAG

Given the constraint enforcement property of RAG, in this section, we
integrate RL with RAG to achieve safe learning without violating system
constraints. Indeed, the proposed safe RL framework can be combined with any
RL algorithms. In this paper, we consider the Neural-Fitted Q-learning (NFQ)
as the base RL module due to its ability to deal with continuous state and
control spaces (Riedmiller, 2005).

In NFQ, we use a neural network to approximate the Q-function, and we update
this Q-function approximation using the following formula,

$Q(x(t),u(t))\leftarrow\lambda\tilde{Q}(x(t),u(t))+(1-\lambda)\big{(}R(x(t),u(t))+\gamma\tilde{V}(x(t+1))\big{)},$
(8)

where $\lambda\in(0,1)$ is the learning rate, $\gamma\in(0,1)$ is a factor to
discount future rewards and avoid the Q-values increasing to infinity, and

$\tilde{V}(x)=\max_{u\in\mathcal{U}}\tilde{Q}(x,u),$ (9)

where $\tilde{Q}(x,u)$ represents the approximated Q-values extracted from the
neural network. Note that the data points used to update Q-function in (8) can
be obtained through black-box simulations or hardware experiments, which makes
the NFQ a model-free method.

The proposed safe RL framework based on RAG is formally presented as Algorithm
2. We first initialize a neural network to approximate Q-function for the
continuous state and action spaces. The action space $\mathcal{U}$ is
discretized with step $d_{s}$ to reduce the computational cost of selecting
the optimal control action $u(t)$ based on currently approximated Q-function
in (9). The $n_{t}$ is the collected trajectory number in one episode, and $T$
is the length of one trajectory. We use a replay buffer $B$ to store the
experience collected in the most recent episode.

The RL agent balances exploration and exploitation using an $\epsilon$-greedy
action selection rule in Lines 5-9. Line 10 is the control modification step.
By solving (5), any action that may lead to constraint violation will be
modified to $u^{\text{safe}}$. Note that we only discretize the action space
$\mathcal{U}$ in Line 8, (5) is solved on the continuous action space.

Lines 11-16 collect experience and store it into replay buffer $B$ for later
use in training the Q-function NN $\tilde{Q}(x,u)$. In particular, we update
the Q-value of the current state $x(t)$ and nominal control $u(t)$,
$Q(x(t),u(t))$, with the reward $R(x(t),u^{\text{safe}}(t))$ brought by safe
action $u^{\text{safe}}(t)$. Furthermore, we store the tuple of current state,
nominal control and updated Q-value, $(x(t),u(t),Q(x(t),\\\ u(t)))$, to the
replay buffer. This way, the action modification by RAG is not perceived by
the agent, and the agent will explore all actions. This may be beneficial as
this will potentially improve convergence to the optimal control policy.

Algorithm 2 Safe RL algorithm

Input Initialized Q-value NN $\tilde{Q}(x,u)$, empty replay buffer $B$,
discretized action space $\mathcal{U}_{d}$, and the maximum trajectory number
within a training episode $N$.

1:for each training episode do

2: while $n_{t}<N$ do

3: Pick a safe initial state $x(0)\in X_{\text{safe}}$

4: while $t<T$ do

5: if $rand()<\epsilon$ then

6: $u(t)$ takes a random value within $\mathcal{U}_{d}$$\triangleright$
Exploration

7: else

8:
$u(t)\in\operatorname*{arg\,max}_{u\in\mathcal{U}_{d}}\tilde{Q}(x(t),u)$$\triangleright$
Exploitation

9: end if

10: Solve (5) to modify $u(t)$ to $u^{\text{safe}}(t)$$\triangleright$ RAG

11: Apply $u^{\text{safe}}(t)$ to system

12: Observe next state $x(t+1)$ and reward $R(x(t),u^{\text{safe}}(t))$

13: Update Q value:

14: $\tilde{V}(x(t+1))=\max_{u\in\mathcal{U}_{d}}\tilde{Q}(x(t+1),u)$

15:
$Q(x(t),u(t))\leftarrow\lambda\tilde{Q}(x(t),u(t))+(1-\lambda)\big{(}R(x(t),u^{\text{safe}}(t))+\gamma\tilde{V}(x(t+1))\big{)}$

16: Store $(x(t),u(t),Q(x(t),u(t)))$ to replay buffer $B$

17: end while

18: end while

19: Train Q-value NN $\tilde{Q}(x,u)$ using sampled data from replay buffer
$B$

20: $n_{t}\leftarrow 0$

21:end for

Note that the proposed RAG approach can be combined with an arbitrary RL
algorithm to achieve safe RL. The discretization step of action space
$\mathcal{U}$ is only for reducing the computational complexity and is not a
necessary step of the proposed safe RL framework. Furthermore, although in
Algorithm 2 the Q-learning algorithm produces a nominal control input that
takes values in the discretized space $\mathcal{U}_{d}$, the modified safe
control input $u^{\text{safe}}$ after RAG takes values in the original
continuous space.

##  5 Safe RL for Adaptive Cruise Control

In this section, we apply the proposed safe RL framework to the Adaptive
Cruise Control (ACC) problem for an automated vehicle, to demonstrate the
major advantage of the proposed safe RL approach over conventional RL in terms
of guaranteed constraint satisfaction both during and after training. The
dynamics of relative motion between the lead vehicle and the following ego
vehicle are represented as

$\displaystyle\begin{split}\begin{bmatrix}\Delta s(k+1)\\\ \Delta v(k+1)\\\
v^{\text{ego}}(k+1)\end{bmatrix}=&\begin{bmatrix}1&T_{s}&0\\\ 0&1&0\\\
0&0&1\end{bmatrix}\begin{bmatrix}\Delta s(k)\\\ \Delta v(k)\\\
v^{\text{ego}}(k)\end{bmatrix}+\begin{bmatrix}-\frac{T_{s}^{2}}{2}\\\
-T_{s}\\\ T_{s}\end{bmatrix}u(k)+\begin{bmatrix}\frac{T_{s}^{2}}{2}\\\
T_{s}\\\ 0\end{bmatrix}w(k),\end{split}$ (10)

where $\Delta s$[m] and $\Delta v$[m/s] are the longitudinal distance and
relative speed between the lead and ego vehicles, respectively. The state
$v^{\text{ego}}$[m/s] represents the ego vehicle’s velocity. The control input
$u$[m/s2] represents the ego vehicle’s acceleration, and to improve the
driving comfort, we impose the following input constraints
$u\in\mathcal{U}=\\{u:-3\leq u\leq 3\\}$. $T_{s}=0.5$s is the sampling period.
We treat the lead vehicle’s acceleration $w\in\mathcal{W}=\\{w:-1.5\leq w\leq
1.5\\}$ as a disturbance. The upper and lower bounds of the lead vehicle’s
acceleration are derived based on the FTP75 driving cycle, which represents
the typical city driving behavior.

Our goal is to achieve a safe and efficient car-following performance. We use
headway time as the metric and set its target value to 1.5s. The reward
function of the RL agent is design as following:

$R=\begin{cases}-(\frac{\Delta s}{v^{\text{ego}}}-1.5)^{2}&\text{if
$v^{\text{ego}}\geq 5$},\\\ -(\Delta s-7.5)^{2}&\text{if $0\leq
v^{\text{ego}}\leq 5$}.\\\ \end{cases}$ (11)

Also, we consider the following constraints on the headway time to enforce
safety:

$1\leq\frac{\Delta s}{\text{max}(v^{\text{ego}},5)}\leq 2.$ (12)

As the set defined by (12) is the region where the system state is supposed to
be in, $X_{0}$ is its complement and can be expressed in the form of (4).

The nominal control $u_{\phi}$ (before training) is a state-feedback control
policy tracking the desired headway time of 2.5s, which is different from the
target, i.e. 1.5s. Two RL schemes, namely, conventional RL and safe RL, as
shown in Figure 1, are employed to train the nominal control policy
$u_{\phi}$. During training, we randomly sample segments (with length of
$T=30$s) within the FTP75 driving cycle as the lead vehicle’s speed
trajectory. The unrecoverable set $X_{k}$ and safe set $X_{\text{safe},k}$ are
computed with the MPT3 toolbox offline (Herceg et al., 2013). The online MIQP
optimization problem in (5) is solved by OPTI with SCIP (Currie and Wilson,
2012; Achterberg, 2009). We compute $X_{k}$ according to Algorithm 1 with
$k=10$.

The training histories of conventional and safe RL algorithms are illustrated
in Figure 2. As the nominal feedback control policy tracks the headway time of
2.5s before training, which is not within the range of constraints in (12),
the constraint violation rate of conventional RL is high at the beginning of
training. In contrast, no constraint violation is exhibited for safe RL during
the entire training process. With RAG, the RL agent can learn the control
policy safely without any constraint violation. Moreover, with RAG the headway
time of ego vehicle is always within the range of $[1,2]$s, leading to a
smaller reward variation and faster learning compared with conventional RL, as
shown in Figure 2b.

Figure 2: Training histories of conventional RL and safe RL. (a) Constraint
violation rates of each episode. (b) Average reward values of each episode.
Solid lines represent average values and shaded areas represent the standard
deviation values over 20 experiments.

The validation results of trained policies are shown in Figure 3. The policy
trained with safe RL is implemented with RAG. As shown in Figure 3a, 3b and
3c, the velocity and headway time tracking performances of conventional and
safe RL policies are both satisfactory. However, under the control of
conventional RL policy, there are still occasional constraint violations in
$\Delta s$, as shown in Figure 3d. In comparison, there is no constraint
violation with the control of safe RL policy, which is attributed to the fact
that the RAG monitors and modifies the control input to guarantee the
constraint satisfaction as illustrated in Figure 3f.

We remark that the fact that our safe RL algorithm based on the use of RAG
guarantees no safety constraint violation during both the training and the
operating phases yields that it can be used for onboard applications. For
instance, it can be used to continuously improve the performance of a
controller during its onboard operation.

Figure 3: Validation results of conventional RL and safe RL. (a, b) Speed
tracking performances. (c) headway time tracking performances (target value:
1.5s). (d) Relative distance between lead vehicle and ego vehicle. (e, f)
Control inputs.

The computational time performances are shown in Figure 4. The simulations are
performed on the Matlab R2019b platform using an Intel Xeon E5-1650 3.50 GHz
PC with Windows 10 and 16.0 GB of RAM. As the optimization problem (5) needs
to be solved at each time instant, safe RL takes around 130 s on average to
complete an episode, which is longer compared with the one with conventional
RL. Furthermore, we can observe that the average online solving time for (5)
is around 30 ms, which is feasible for real-time control scenarios.

Figure 4: Comparison of computational times. (a) Training time of each
episode. (b) Computational time of each step during validation. Solid lines
represent average values and shaded areas represent the standard deviation
values over 20 experiments.

##  6 Conclusions and Future Work

In this paper, we developed a safe RL framework based on RAG that integrates
model-based safety supervision and model-free learning. We exploited the
underlying dynamics and exclusion-zone requirement to construct a safety set
for constraining learning exploration using set-theoretic techniques and
online optimization. We applied the proposed safe RL framework to an Adaptive
Cruise Control system and showed that we could conduct the online learning
with no safety constraint violations. Future work will include the
investigation into the computational complexity, scalability, and approaches
to improving the exploration efficiency of the proposed safe RL framework.

This work is supported by Ford Motor Company. We also thank Mr. Xintao Yan for
helpful discussions about implementing reinforcement learning algorithm in the
ACC example.

## References

  * Abbeel et al. (2007) Pieter Abbeel, Adam Coates, Morgan Quigley, and Andrew Y Ng.  An application of reinforcement learning to aerobatic helicopter flight.  In _Advances in neural information processing systems_ , pages 1–8, 2007.
  * Achterberg (2009) Tobias Achterberg.  SCIP: solving constraint integer programs.  _Mathematical Programming Computation_ , 1(1):1–41, Jul 2009.  ISSN 1867-2957.  10.1007/s12532-008-0001-1.  URL https://doi.org/10.1007/s12532-008-0001-1.
  * Aswani et al. (2013) Anil Aswani, Humberto Gonzalez, S Shankar Sastry, and Claire Tomlin.  Provably safe and robust learning-based model predictive control.  _Automatica_ , 49(5):1216–1226, 2013.
  * Cheng et al. (2019) Richard Cheng, Gábor Orosz, Richard M Murray, and Joel W Burdick.  End-to-end safe reinforcement learning through barrier functions for safety-critical continuous control tasks.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 33, pages 3387–3395, 2019.
  * Currie and Wilson (2012) Jonathan Currie and David I. Wilson.  OPTI: Lowering the Barrier Between Open Source Optimizers and the Industrial MATLAB User.  In Nick Sahinidis and Jose Pinto, editors, _Foundations of Computer-Aided Process Operations_ , Savannah, Georgia, USA, 8–11 January 2012\.
  * Fisac et al. (2018) Jaime F Fisac, Anayo K Akametalu, Melanie N Zeilinger, Shahab Kaynama, Jeremy Gillula, and Claire J Tomlin.  A general safety framework for learning-based control in uncertain robotic systems.  _IEEE Transactions on Automatic Control_ , 64(7):2737–2752, 2018.
  * Garcıa and Fernández (2015) Javier Garcıa and Fernando Fernández.  A comprehensive survey on safe reinforcement learning.  _Journal of Machine Learning Research_ , 16(1):1437–1480, 2015.
  * Geibel and Wysotzki (2005) Peter Geibel and Fritz Wysotzki.  Risk-sensitive reinforcement learning applied to control under constraints.  _Journal of Artificial Intelligence Research_ , 24:81–108, 2005.
  * Herceg et al. (2013) M. Herceg, M. Kvasnica, C.N. Jones, and M. Morari.  Multi-Parametric Toolbox 3.0.  In _Proc. of the European Control Conference_ , pages 502–510, Zürich, Switzerland, July 17–19 2013.  http://control.ee.ethz.ch/ mpt.
  * Kober et al. (2013) Jens Kober, J Andrew Bagnell, and Jan Peters.  Reinforcement learning in robotics: A survey.  _The International Journal of Robotics Research_ , 32(11):1238–1274, 2013.
  * Larsen et al. (2017) Rie B Larsen, Andrea Carron, and Melanie N Zeilinger.  Safe learning for distributed systems with bounded uncertainties.  _IFAC-PapersOnLine_ , 50(1):2536–2542, 2017.
  * Li et al. (2021a) N. Li, K. Han, A. Girard, H. E. Tseng, D. Filev, and I. Kolmanovsky.  Action governor for discrete-time linear systems with non-convex constraints.  _IEEE Control Systems Letters_ , 5(1):121–126, 2021a.  10.1109/LCSYS.2020.3000198.
  * Li et al. (2021b) Y. Li, N. Li, H. E. Tseng, A. Girard, D. Filev, and I. Kolmanovsky.  Robust action governor for discrete-time piecewise affine systems with additive disturbances.  _IEEE Control Systems Letters_ , 2021b.
  * Li et al. (2018) Zhaojian Li, Uroš Kalabić, and Tianshu Chu.  Safe reinforcement learning: Learning with supervision using a constraint-admissible set.  In _2018 Annual American Control Conference (ACC)_ , pages 6390–6395. IEEE, 2018.
  * Li et al. (2019) Zhaojian Li, Tianshu Chu, and Uroš Kalabić.  Dynamics-enabled safe deep reinforcement learning: Case study on active suspension control.  In _2019 IEEE Conference on Control Technology and Applications (CCTA)_ , pages 585–591. IEEE, 2019.
  * Mnih et al. (2013) Volodymyr Mnih, Koray Kavukcuoglu, David Silver, Alex Graves, Ioannis Antonoglou, Daan Wierstra, and Martin Riedmiller.  Playing atari with deep reinforcement learning.  _arXiv preprint arXiv:1312.5602_ , 2013.
  * Moldovan and Abbeel (2012) Teodor Mihai Moldovan and Pieter Abbeel.  Safe exploration in markov decision processes.  _arXiv preprint arXiv:1205.4810_ , 2012.
  * Riedmiller (2005) Martin Riedmiller.  Neural fitted Q iteration–first experiences with a data efficient neural reinforcement learning method.  In _European Conference on Machine Learning_ , pages 317–328. Springer, 2005.
  * Sloth et al. (2012) Christoffer Sloth, George J Pappas, and Rafael Wisniewski.  Compositional safety analysis using barrier certificates.  In _Proceedings of the 15th ACM international conference on Hybrid Systems: Computation and Control_ , pages 15–24, 2012.
  * Sutton and Barto (2018) Richard S Sutton and Andrew G Barto.  _Reinforcement learning: An introduction_.  MIT press, 2018.
  * Wachi et al. (2018) Akifumi Wachi, Yanan Sui, Yisong Yue, and Masahiro Ono.  Safe exploration and optimization of constrained mdps using gaussian processes.  In _AAAI_ , pages 6548–6556, 2018.