Stripe skyrmions and skyrmion crystals

X. R. Wang1,2,∗, X. C. Hu1,2, and H. T. Wu1,2

  1. 1.

Physics Department, The Hong Kong University of Science and Technology, Clear
Water Bay, Kowloon, Hong Kong.

  2. 2.

HKUST Shenzhen Research Institute, Shenzhen 518057, China.

Correspondence and requests for materials should be addressed to X.R.W.
(email<EMAIL_ADDRESS>
Skyrmions are important in topological quantum field theory for being soliton
solutions of a nonlinear sigma model and in information technology for their
attractive applications. Skyrmions are believed to be circular and stripy spin
textures appeared in the vicinity of skyrmion crystals are termed spiral,
helical, and cycloid spin orders, but not skyrmions. Here we present
convincing evidences showing that those stripy spin textures are skyrmions,
“siblings” of circular skyrmions in skyrmion crystals and “cousins” of
isolated circular skyrmions. Specifically, isolated skyrmions are excitations
when skyrmion formation energy is positive. The skyrmion morphologies are
various stripy structures when the ground states of chiral magnetic films are
skyrmions. The density of skyrmion number determines the morphology of
condensed skyrmion states. At the extreme of one skyrmion in the whole sample,
the skyrmion is a ramified stripe. As the skyrmion number density increases,
individual skyrmion shapes gradually change from ramified stripes to
rectangular stripes, and eventually to disk-like objects. At a low skyrmion
number density, the natural width of stripes is proportional to the ratio
between the exchange stiffness constant and Dzyaloshinskii-Moriya interaction
coefficient. At a high skyrmion number density, skyrmion crystals are the
preferred states. Our findings reveal the nature and properties of stripy spin
texture, and open a new avenue for manipulating skyrmions, especially
condensed skyrmions such as skyrmion crystals.

Skyrmions, originally used to describe resonance states of baryons [1], were
first unambiguously observed in the form of skyrmion crystals in various
chiral magnets and by various experimental techniques [2, 3, 4, 5, 6, 7].
Magnetic skyrmions are topological non-trivial spin textures of magnetic films
characterized by a skyrmion number of
$Q=\frac{1}{4\pi}\int\vec{m}\cdot(\partial_{x}\vec{m}\times\partial_{y}\vec{m}){\rm
d}x{\rm d}y$, here $\vec{m}$ is the unit vector of magnetization. $Q$ must be
an integer for an infinite magnetic film, and a non-zero $Q$ spin texture is
called a skyrmion of skyrmion number $Q$. It was isolated circular skyrmions,
not a skyrmion crystal, that were predicted in early theories [8, 9]. Isolated
skyrmions were indeed observed later in confined structures and films [10, 11,
12, 13, 14, 15]. It is an experimental fact that skyrmion crystals form in
very narrow magnetic-field-temperature windows. Outside of the windows, stripy
phases appear. The stripy phases, which can even coexist with skyrmion
crystals, are in fact easier to form than a skyrmion crystal does. In
contrast, those stripy phases do not appear together with isolated skyrmions.
Interestingly, the stripy phase was observed many years before the
observations of skyrmion crystals [16]. These stripy phases are called
helical, spiral, and cycloid spin orders. A one-dimensional model [9, 17, 18,
19, 20, 21] was used to describe the rotation of spins perpendicular to
stripes. To date, a holistic description of various stripy structures,
especially ramified stripes and stripy maze, [22, 23, 24, 25, 26, 27, 28, 29]
is lacking. The general belief is that those stripes are not skyrmions and
have zero skyrmion numbers [30] although some race-track stripes are called
merons with 1/2 skyrmion number or bimerons[31, 32, 27].

In this letter we show that stripy magnetic textures appear in a chiral
magnetic film with Dzyaloshinskii-Moriya interaction (DMI) are actually
irregular skyrmions. Each stripy texture has exactly skyrmion number 1. The
irregular shape is due to the negative skyrmion formation energy (relative to
ferromagnetic state) when the ferromagnetic state is not the ground state. For
a given system, the morphology of these skyrmions are random when the skyrmion
number density is low. At extremely low density, magnetic textures are
ramified stripes. The exact appearance of each pattern is very sensitive to
the initial spin configuration and its dynamical path. The basic building
blocks of irregular random skyrmions are stripes of well-defined width. The
optimal width comes from the competition between the Heisenberg’s exchange
energy and the DMI energy that respectively prefer a larger and smaller width
in order to minimize exchange energy cost and maximize the negative skyrmion
formation energy gain. Unexpectedly, this exchange energy and DMI energy
dependence of width is opposite to the skyrmion size of an isolated skyrmion
that increases with DMI interaction and decreases with exchange energy [34].

We consider an ultra-thin ferromagnetic film of thickness $d$ in the $xy$
plane. The film has an exchange energy $E_{\mathrm{ex}}$ with exchange
stiffness constant $A$, an interfacial DMI energy $E_{\mathrm{DM}}$ with DMI
coefficient $D$, an anisotropy energy $E_{\mathrm{an}}$ with a perpendicular
easy-axis anisotropy $K$, and the Zeeman energy $E_{\mathrm{Ze}}$ in a
perpendicular magnetic field $H$. The total energy $E_{total}$ reads

$E_{total}=E_{\mathrm{ex}}+E_{\mathrm{DM}}+E_{\mathrm{an}}+E_{\mathrm{Ze}},$
(1)

where $E_{\mathrm{ex}}=Ad\iint|\nabla\vec{m}|^{2}\mathrm{d}S$,
$E_{\mathrm{DM}}=Dd\iint[m_{z}\nabla\cdot\vec{m}-(\vec{m}\cdot\nabla)m_{z}]\mathrm{d}S$,
$E_{\mathrm{an}}=Kd\iint(1-m_{z}^{2})\mathrm{d}S$, and
$E_{\mathrm{Ze}}=\mu_{0}HM_{\text{s}}d\iint(1-m_{z})\mathrm{d}S$.
$M_{\text{s}}$ is the saturation magnetization and $\mu_{0}$ is the vacuum
permeability. The integration is over the whole film. The energy is set to
zero, $E_{total}=0$, for the ferromagnetic state of $m_{z}=1$. We have assumed
$\vec{m}$ is uniform in thickness direction. The demagnetization effect is
included in $E_{\mathrm{an}}$ through the effective anisotropy
$K=K_{u}-\mu_{0}M_{\text{s}}^{2}/2$ corrected by the shape anisotropy, here
$K_{u}$ is the perpendicular magnetocrystalline anisotropy. This is a good
approximation when the film thickness $d$ is much smaller than the exchange
length [34]. It is known that isolated circular skyrmions are metastable state
of energy $8\pi Ad\sqrt{1-\kappa}$ when $\kappa=\pi^{2}D^{2}/(16AK)<1$ [34].
Here we use MuMax3 simulator [35] to numerically solve the Landau-Lifshitz-
Gilbert (LLG) equation for the stable states in the opposite regime of
$\kappa>1$ where circular skyrmions are not stable states [34], see Method.

Figure 1: Numerical solutions of LLG equation under the periodical boundary
conditions and with $A=10{\rm pJ/m},D=6{\rm mJ/m^{2}},K=0.49{\rm
MJ/m^{3}},M_{s}=0.58{\rm MA/m}$ for a sample of $300\rm{nm}\times
300\rm{nm}\times 0.4\rm{nm}$. (a1,b1,c1) are different initial configurations
of a disk of diameter $20\rm{nm}$ (a1), a hexagon of side length $10\rm{nm}$
(b1) and a square of length $20\rm{nm}$ (c1). (a2,b2,c2) are intermediate
states at 0.3ns with irregular shapes due to the negative formation energy.
(a3,b3,c3) are the final stable pattern with irregular ramified stripes.
Skyrmion charge density $\rho$ is encoded by colours (the blue for positive
and the red for negative) while the gray-scale encodes $m_{z}$. The dark black
lines denote $m_{z}=0$. The positive and negative charges exist respectively
only around convex and concave areas. The spin profile across the stripes at
the green ⓝ is considered. (a4, b4, c4) are the evolution of total energy
$E_{total}$ and topological skyrmion number $Q$ ($t$ is in the logarithmic
scale). $Q$ reaches the skyrmion number 1 within $1ps$. Clearly, the skyrmion
number is a constant and the total energy is negative and approaches a
constant almost independent from the initial configurations.

For a sample of $300\rm{nm}\times 300\rm{nm}\times 0.4\rm{nm}$ under the
periodical boundary conditions and with $A=10{\rm
pJ/m},D=6\rm{mJ/m^{2}},K=0.49\rm{MJ/m^{3}}$, and $M_{s}=0.58\rm{MA/m}$ that
are typical values for those chiral magnets supporting skyrmion crystals [36,
37, 30], $\kappa>1$, so that single domain and isolated circular skyrmion are
not stable any more [34]. We will start from a small nucleation magnetic
domain with sharp domain wall to speed up skyrmion formation dynamics although
stripe skyrmions will appear spontaneously due to thermal or other
fluctuations in reality. Figure 1 shows how various initial states in the top
panel (a1, b1, c1) evolve according to the LLG equation with $\alpha=0.25$.
$\alpha$ will not change physics described below because we are interested in
the spin textures of the LLG equation that do not vary with time. However,
$\alpha$ can change the evolution path and energy dissipation rate so that it
can change the intermediate states shown in (a2-c2) and influence which fixed
point to pick when a system has many fixed points (textures) like the current
situation. This is similar to the sensitiveness of attractor basins to the
damping in a macro spin system [38]. The initial configurations are obtained
by reversing spins in the white regions in (a1-c1) from $m_{z}=-1$ to
$m_{z}=1$ such that the configurations have very sharp domain walls and zero
skyrmion number $Q=0$. After a short time of an order of picosecond, the
initial states transform into irregular structures of skyrmion number $Q=1$,
no matter whether the initial shape is circular or non-circular as shown in
(a4-c4) where the time is in logarithmic scale. As time goes on, $Q$ (the blue
lines in a4-c4) stays at 1, and system energy $E_{total}$ (the red curves) is
negative and keeps decreasing until it reaches a stable ramified stripy spin
texture. Clearly, this irregular ramified spin texture is a non-circular
skyrmion whose formation energy is negative. The negative formation energy
explains why the skyrmion prefers to stretch out to occupy the whole space to
lower its energy. This process is clearly demonstrated in the movie in the
Supplementary Material, showing how the system evolves from (a1) to (a3) and
how $Q$ grows from 0 to 1. The simulations show also that the exact pattern is
very sensitive to the initial configuration and dynamical/kinetic paths as
well as how energy dissipates. In a real system, the process shall proceed
spontaneously from any randomly generated nucleation centre no matter by a
thermal fluctuation or an intentional agitation.

We use color to encode the skyrmion charge density $\rho$ defined as
$\rho=\vec{m}\cdot(\partial_{x}\vec{m}\times\partial_{y}\vec{m})/(4\pi)$ as
indicated by the color bar in Fig. 1. Interestingly, $\rho$ is non-zero only
around convex (positive) and concave (negative) areas. $\rho$ is almost zero
along the long straight stripe. This may explain why the skyrmion number of
stripy textures were thought to be zero in the literature [30]. As shown in
Fig. 1(a4-c4), the sum of all positive and negative skyrmion charge is always
quantized to 1 protected by the topology. The continuous decrease of energy
$E_{total}$ also indicate the final morphology of the ramified stripe
skyrmions is not unique and depends on dynamical path of the system evolution
that in turn relates to the initial configuration.

One striking feature about the stripes shown in Fig. 1 is a well-defined
stripe width that is usually referred to a given wave vector in experiments.
The competition between the exchange energy and DMI energy can easily lead to
a $A/D$ dependence [9, 17, 18, 19, 20]. However, it is highly non-trivial to
understand the effect of the magnetic anisotropy on stripe width. In order to
understand the underlying physics of this well-defined width, one needs a good
spin texture profile along the direction perpendicular ($x$) to the stripes.
We found that black ($m_{z}\leq 0$) and white ($m_{z}\geq 0$) stripes can be
approximated by
$\Theta(x)=2\arctan\left[\frac{\sinh(L/2w)}{\sinh(|x|/w)}\right]$ and
$\Theta(x)=2\arctan\left[\frac{\sinh(|x|/w)}{\sinh(L/2w)}\right]$ ($|x|\leq
L/2$), respectively. $\Theta$ is the polar angle of the magnetization at
position $x$ and $x=0$ is the centre of a stripe. $L$ and $w$ measure
respectively the stripe width and skyrmion wall thickness as schematically
illustrated in Fig. 2(a). Figure 2(b) demonstrates the excellence of this
approximate profile for a set of model parameters of
$A=10\rm{pJ/m},D=6\rm{mJ/m^{2}},K=0.49\rm{MJ/m^{3}}$, and
$M_{s}=0.58\rm{MA/m}$. The $y-$axis is $m_{z}$ and $x=0$ is the stripe centre
where $m_{z}=1$. Different symbols are numerical data across different stripes
labelled by the green ⓝ in Fig. 1(a3-c3). The solid curve is the fit of
profile $\cos\Theta(x)$ with $L=10.78\rm{nm}$ and $w=3.05\rm{nm}$. All data
from different stripes falling onto the same curve demonstrates that stripes,
building blocks of pattern, are identical.

Using the excellent spin profile, one can obtain magnetic energy of a film
filled with the stripe skyrmions as a function of $L$ and $w$. Minimizing the
energy against $L$ and $w$ allows us to obtain $A$, $D$, and $K$ dependence of
stripe width $L$ and skyrmion wall thickness $w$. In terms of $L$,
$\epsilon=L/(2w)$, $\xi=A/D$, and $\kappa$, the total system energy density of
a film filled by such stripes is

$E=\frac{4D}{\xi}\left[\frac{\xi^{2}}{L^{2}}g_{1}(\epsilon)-\frac{\pi\xi}{4L}+\frac{\pi^{2}}{64\kappa}g_{2}(\epsilon)\right],$
(2)

where
$g_{1}(\epsilon)=\int_{0}^{1}\frac{[2\epsilon\sinh(\epsilon)\cosh(\epsilon
x)]^{2}}{[\sinh^{2}(\epsilon)+\cosh^{2}(\epsilon x)]^{2}}{\mathrm{d}}x$, and
$g_{2}(\epsilon)=\int_{0}^{1}\frac{[2\sinh(\epsilon)\sinh(\epsilon
x)]^{2}}{[\sinh^{2}(\epsilon)+\cosh^{2}(\epsilon x)]^{2}}{\mathrm{d}}x$. The
optimal stripe width obtained from minimizing energy $E$ is $L=a\frac{4A}{\pi
D}$, where $a$ depends weakly on $\kappa=\pi^{2}D^{2}/(16AK)>1$. The physics
of this result is clear: DMI energy is negative, and one can add more stripes
by reducing $L$ such that the total energy will be lowered. On the other hand,
the exchange energy will increase with the decrease of $L$. As a result, $L$
is proportional to $A/D$ that is opposite to the behaviour of size of an
isolated skyrmion whose size increases with $D$ and decreases with $A$ [34,
30]. These theoretical results agree very well with micro-magnetic simulations
as shown in Fig. 2(c) with $a=7.61$ for $\kappa=1.3$; $a=6.82$ for
$\kappa=1.8$; $a=6.47$ for $\kappa=2.9$. The dependence of $L$ on $\kappa$ for
a fixed $A/D=1\rm{nm}$ is shown in the inset. One can see that $L$ depends
weakly on large $\kappa\gg 1$. Unexpectedly, $L$ are the same no matter
whether we have one, two, or more ramified and non-ramified stripe skyrmions.
This is reflected in the skyrmion number density independence of $L$ as shown
in Fig. 2(d).

Figure 2: (a) Schematic diagram of spin texture of parallel stripes. (b)
Symbols, from different stripes labelled by textcircledn in Fig. 1(a3-c3), are
$m_{z}=\cos\Theta$ from MuMax3. The green solid line is the fit of
$\Theta(x)=2\arctan\left[\frac{\sinh(|x|/w)}{\sinh(L/2w)}\right]$.
($\Theta\leq\pi/2$) and
$\Theta(x)=2\arctan\left[\frac{\sinh(L/2w)}{\sinh((|x-L|))/w)}\right]$
($\Theta\geq\pi/2$) with $L=10.78\rm{nm}$ and $w=3.05\rm{nm}$. $x=0$ is the
centres of white stripes. All data from different stripes fall onto the same
curve means that stripes of the same width are basic building blocks of stripe
skyrmions. (c) $A/D$-dependence of stripe width $L$ for various
$\kappa=\pi^{2}D^{2}/(16AK)=1.3$ (the green stars); 1.8 (the blue squares);
and 2.9 (the red circles). The solid lines are the fits of $L=aA/D$ with
$a=7.61$ for $\kappa=1.3$; $a=6.82$ for $\kappa=1.8$; $a=6.47$ for
$\kappa=2.9$. Inset: the dependence of $L$ on $\kappa$ for $A/D=1\rm{nm}$.
Symbols are the numerical data and the dashed line is the theoretical
prediction without any fitting parameter. (d) The dependences of stripe width
$L$ on skyrmion number density for
$A=10\rm{pJ/m},D=6\rm{mJ/m^{2}},K=0.49\rm{MJ/m^{3}}$, and
$M_{s}=0.58\rm{MA/m}$. The number of skyrmions varies from 1 to 113 in our
sample of $300\rm{nm}\times 300\rm{nm}\times 0.4\rm{nm}$ under the periodical
boundary conditions. The dash line is $L=10.52\rm{nm}$.

To understand why condensed stripe skyrmions and skyrmion crystals can appear
together in a given chiral magnetic film when the material parameters are
fixed, we try to increase the number of skyrmions in our film of
$300\rm{nm}\times 300\rm{nm}\times 0.4\rm{nm}$. Encouraged by the results in
Fig. 1 that each nucleation domain creates one irregular stripe skyrmion, we
place 2, 100 and 169 small disk domains of diameter $10\rm{nm}$ as shown in
Fig. 3(a1-c1) and let them to evolve according to LLG equation. Fig. 3(a2-c2)
are the final steady states. As expected, we indeed obtained two irregular
ramified stripe skyrmions (a2). For the case of 100 skyrmions, some of them
have rectangular shape and are arranged in a nematic phase while the rest of
skyrmions look like disks and are in a lattice structure (b2). In the case of
169 skyrmions, skyrmions are disk-like and are in a triangular lattice (c2).
The skyrmion nature of the spin textures can be clearly confirmed by the
change of $Q$ (blue) with time as shown in Fig. 3(d) (the solid, dashed, and
dotted lines for $Q=2,\ 100,$ and 169 respectively). (d) shows also how system
energy (the red lines) changes with time. One interesting feature is that
total energy is not sensitive to skyrmion number density before skyrmion-
skyrmion distance is comparable to the optimal stripe width and skyrmions take
stripe shape. The system energy starts to increase with the skyrmion number
density, and condensed skyrmions transform from rectangular stripe skyrmions
in nematic phase into circular skyrmion crystal. This feature does not favour
skyrmion crystal formation, and may explain why skyrmion crystals were
observed in the presence of an external magnetic field: With only one skyrmion
in the whole system as those shown in Fig. 1, the net magnetic moment is zero
because spin up (white area) and down spin (black area) are almost equal so
that magnetic moment cancel each other. As skyrmion number increases, the net
magnetic moment increases and becomes non-zero.

Figure 3: Evolution of systems with different number of single domains $N=2$
(a), 100 (b) and 169 (c) spread in the sample of $300\rm{nm}\times
300\rm{nm}\times 0.4\rm{nm}$ with the periodical boundary conditions. The
material parameters are the same as those in Fig. 1. To speed up the system
reaching its final stable states, we use a large Gilbert damping constant of
$\alpha=1$. (a1,b1,c1) are the initial configurations and (a2,b2,c2) are the
final stable pattern. skyrmion charge density $\rho$ is encoded in colours
(the red and the blue for positive and negative $\rho$ respectively) while
$m_{z}$ is encoded in the grayscale. The dark black lines denote $m_{z}=0$.
The positive and negative charges exist respectively only around canvex and
concave areas. (d) is the evolution of total energy $E_{total}$ (the left
y-axis) and topological skyrmion number $Q$ (the right y-axis with three
ranges). Clearly, $Q$ reaches a constant integer in a very short time and
$Q=N=2,100$ and $169$ respectively. The total energy is negative and
approaches a constant, indicating a stable state.

We have showed that stripe and ramified stripe skyrmions are essentially the
same as the circular skyrmions in skyrmion crystals. The difference in
skyrmion shapes at different skyrmion number density come from the skyrmion-
skyrmion interaction. When the average distance between two nearby skyrmions
is order of the stripe width, the skyrmion-skyrmion repulsion compress
skyrmions into circular objects. This understanding permits a skyrmion crystal
in the absence of an external magnetic field as long as one can use other
means to add more skyrmions into a film such as a scanning tunnelling tip [6].

We have presented results for the interfacial DMI so far, similar results are
also true for bulk DMI. As shown in Fig. 4, one ramified stripe skyrmion has
very similar structures for interfacial (a) and bulk (b) DMIs when we start
from the same initial configuration with the same interaction strength. The
only difference is the change of Neel-type of stripe wall to the Bloch-type
stripe wall as shown by the color-coding in the figures.

A perpendicular magnetic field can modify the morphology and width of a stripe
without changing its skyrmion number. Stripe width increases (decreases) with
the field when out-of-plane spin component is anti-parallel (parallel) to the
field.

Figure 4: Structures of one stripe skyrmion for interfacial (a) and bulk (b)
DMIs starting from the same initial configuration. The sample size is
$300\rm{nm}\times 300\rm{nm}\times 0.4\rm{nm}$ with the same model parameters
as those for Fig. 1 except $D=4{\rm mJ/m^{2}}$. The periodical boundary
conditions are used in the simulations. The in-plane spin component is encoded
by the colors.

It requires certain amount energy to destroy a stable or metastable state by
definition. Thus, all skyrmions discussed in this paper are stable against
thermal noise as long as the thermal energy is smaller than their potential
barriers. In the Supplemental Material, we provide a video to show that the
state in Fig. 3(b2) is stable at 50K.

Skyrmions provide a fertile ground for studying fundamental physics. For
example, the topological Hall effect is a phenomenon about how non-collinear
magnetization in skyrmion crystals affect electron transport. Knowing stripy
phases are also condensed irregular skyrmions, it expands surely arena of
topological Hall effect. We can not only study how electron transport be
affected by the skyrmion crystals, but also skyrmions in other condensed
phases such as nematic phases, or how the elongation and orientation of stripe
skyrmion affect Hall transport. With the new discovery of stripe skyrmions, it
will also allow us to investigate the interplay of topology, shape, spin and
charge. One can investigate how the topology, the local and global geometries
affect spin excitations separately.

The assertion that ramified stripes and other stripes appeared together with
skyrmion crystals are irregular skyrmions are firmly confirmed by the nano-
magnetic simulations. These stripes have a well-defined width that is from the
competition between exchange interaction energy and DMI energy. The inverse of
the width is the well-known wave-vector used to describe the spiral spin order
in the literature [16]. In contrast to isolated skyrmions whose size increases
with DMI constant and decreases with the exchange stiffness constant, the
stripe width increases with the exchange stiffness constant and decreases with
the DMI constant. Counter-intuitively, skyrmion crystals are highly
compressible like a gas, not like an atomic crystal. This detail property need
a careful further study. We believe that our findings should have profound
implications in skyrmion-based applications.

Methods
Numerical simulations. Spin dynamics is governed by the Landau-Lifshitz-
Gilbert (LLG) equation,

$\frac{\partial\vec{m}}{\partial t}=-\gamma\vec{m}\times\vec{H}_{\rm
eff}+\alpha\vec{m}\times\frac{\partial\vec{m}}{\partial t},$ (3)

where $\vec{m}$, $\gamma$, $\alpha$ are respectively the unit vector of the
magnetization, gyromagnetic ratio, and the Gilbert damping. $\vec{H}_{\rm
eff}=2A\nabla^{2}\vec{m}+2Km_{z}\hat{z}+\vec{H}_{d}+\vec{H}_{\rm DM}$ is the
effective field including the exchange field characterized by the exchange
stiffness $A$, crystalline anisotropy field, demagnetizing field
$\vec{H}_{d}$, and DMI field $\vec{H}_{\rm DM}$.

In the absence of energy source like an electric current, LLG equation
describe a dissipative system whose energy can only decrease [39, 40]. Thus,
solving LLG equation is an efficient way to find the stable spin textures. In
the present case, we apply periodic boundary conditions to eliminate the edge
effects. We use the Mumax3 package [35] to numerically solve the LLG equation
with mesh size of $1\rm{nm}\times 1\rm{nm}\times 0.4\rm{nm}$ for those
skyrmions of $L>5\rm{nm}$. For skyrmions of $L<5\rm{nm}$ , the mesh size is
$0.1\rm{nm}\times 0.1\rm{nm}\times 0.4\rm{nm}$. The number of stable states
and their structures should not depend on the Gilbert damping constant, but
spin dynamics is very sensitive to $\alpha$. To speed up our simulations, we
use large $\alpha$ of 0.25 and 1. We consider only material parameters that
supports condensed skyrmion states. The skyrmion size $L$ is obtained directly
from numerical data or by fitting the simulated spin profile to
$\Theta(x)=2\arctan\left[\frac{\sinh(L/2w)}{\sinh(|x|/w)}\right]$ (the black
stripes for $m_{z}\leq 0$) and
$\Theta(x)=2\arctan\left[\frac{\sinh(|x|/w)}{\sinh(L/2w)}\right]$ (white
stripes for $m_{z}\geq 0$) with stripe width $L$ and wall width $w$. Here
$-L/2\leq x\leq L/2$ and $x=0$ is the centre of a stripe. It should be pointed
out that the physics discussed here does not depend on the boundary
conditions. However, different boundary conditions have different confinement
potentials (or potential well depth) for stripe skyrmions, and can affect the
maximal skyrmion number density below which a condensed skyrmion state is
metastable.

Energy density of a film filled by stripes.

Following similar assumptions as those in Ref. [19], we can computer the
energy density of a system filled with stripes of width $L$ that parallel to
the $y-$axis and are periodically arranged along the $x-$axis. If $\Theta$ is
the polar angle of spins and spin profile is
$\Theta(x)=2\arctan\left[\frac{\sinh((|x-nL|)/w)}{\sinh(L/2w)}\right]+n\pi$
with integer $n$ and $x\in((n-0.5)L,(n+0.5)L)$ as shown in Fig. 2(a), the
total energy of the film in the range of $y_{1}<y<y_{2}$ and $x_{1}<x<x_{2}$
is

$\displaystyle
E_{total}=d\int_{y_{1}}^{y_{2}}\int_{x_{1}}^{x_{2}}[A(\partial_{x}\Theta)^{2}-D\partial_{x}\Theta+K\sin^{2}\Theta]{\rm
d}x{\rm d}y.$

$x_{2}-x_{1}$ and $y_{2}-y_{1}$ are much bigger than $L$, we need only to
minimize energy density $E=E_{total}/[d(x_{2}-x_{1})(y_{2}-y_{1})]$. Then one
has

$E=\frac{1}{L}\int_{-L/2}^{L/2}\left[A\left(\frac{\partial\Theta}{\partial
x}\right)^{2}-D\frac{\partial\Theta}{\partial x}+K\sin^{2}\Theta\right]{\rm
d}x.$ (4)

In terms of $\epsilon=L/(2w)$, terms in $E$ are,

$\displaystyle E_{ex}$
$\displaystyle=\int_{-L/2}^{L/2}A\left[\frac{\partial\Theta(x)}{\partial
x}\right]^{2}{\mathrm{d}}x$
$\displaystyle=\frac{A}{w^{2}}\int_{-L/2}^{L/2}\left[\frac{2\sinh(L/w)\cosh(x/w)}{\sinh^{2}(L/w)+\sinh^{2}(x/w)}\right]^{2}{\mathrm{d}}x$
$\displaystyle=\frac{2A}{L}\int_{-1}^{1}\left[\frac{2\epsilon\sinh(\epsilon)\cosh(\epsilon
x)}{\sinh^{2}(\epsilon)+\sinh^{2}(\epsilon x)}\right]^{2}{\mathrm{d}}x,$
$\displaystyle=\frac{4A}{L}g_{1}(\epsilon)$

$\displaystyle E_{DM}$
$\displaystyle=-D\int_{-L/2}^{L/2}\frac{\partial\Theta(x)}{\partial
x}{\mathrm{d}}x$ $\displaystyle=-D\left[\Theta(L/2)-\Theta(-L/2)\right]=-\pi
D,$ $\displaystyle E_{an}$
$\displaystyle=\int_{-L/2}^{L/2}K\sin^{2}\Theta(x){\mathrm{d}}x$
$\displaystyle=K\int_{-L/2}^{L/2}\left[\frac{2\sinh(L/2w)\sinh(x/w)}{\sinh^{2}(L/2w)+\sinh^{2}(x/w)}\right]^{2}{\mathrm{d}}x$
$\displaystyle=\frac{K}{2}\int_{-1}^{1}\left[\frac{2\sinh(\epsilon)\sinh(\epsilon
x)}{\sinh^{2}(\epsilon)+\sinh^{2}(\epsilon x)}\right]^{2}L{\mathrm{d}}x$
$\displaystyle=KLg_{2}(\epsilon).$

Add the three terms up, the energy density $E$ is

$E=\frac{4A}{L^{2}}g_{1}(\epsilon)-\frac{\pi D}{L}+Kg_{2}(\epsilon).$

In terms of $L$, $\epsilon$, $\xi=A/D$ and $\kappa=\pi^{2}D^{2}/(16AK)$, we
obtain Eq. (2),

$E=\frac{4D}{\xi}\left[\frac{\xi^{2}}{L^{2}}g_{1}(\epsilon)-\frac{\pi\xi}{4L}+\frac{\pi^{2}}{64\kappa}g_{2}(\epsilon)\right].$

The first and the third terms on the right hand side are positive, thus

$E\geq\frac{4D}{\xi}\left[\sqrt{\frac{4\pi^{2}\xi^{2}}{64L^{2}\kappa}g_{1}(\epsilon)g_{2}(\epsilon)}-\frac{\pi\xi}{4L}\right]=\frac{\pi
D}{L}\left[\sqrt{\frac{g_{1}g_{2}}{\kappa}}-1\right].$

To have negative $E$, $\kappa$ must be larger than $\sqrt{g_{1}g_{2}}\geq 1$.

Data availability. The data that support the plots within this paper and other
findings of this study are available from the corresponding author on
reasonable request.

References

## References

  * [1] Skyrme, T. H. R. A unified field theory of mesons and baryons. Nucl. Phys. 31, 556-569 (1962).
  * [2] Mühlbauer, S. et al. Skyrmion lattice in a chiral magnet. Science 323, 915-919 (2009).
  * [3] Yu, X. Z. et al. Real-space observation of a two-dimensional skyrmion crystal. Nature 465, 901-904 (2010).
  * [4] Yu, X. Z. et al. Near room temperature formation of a skyrmion crystal in thin-films of the helimagnet FeGe. Nat. Mater. 10, 106-109 (2011).
  * [5] Heinze, S. et al. Spontaneous atomic-scale Magnetic skyrmion lattice in two dimensions. Nat. Phys. 7, 713-718 (2011).
  * [6] Romming, N. et al. Writing and deleting single magnetic skyrmions. Science 341, 636-639 (2013).
  * [7] Onose, Y., Okamura, Y., Seki, S., Ishiwata, S. & Tokura, Y. Observation of magnetic excitations of skyrmion crystal in a helimagnetic insulator $\mathrm{Cu}_{2}\mathrm{OSeO}_{3}$. Phys. Rev. Lett. 109, 037603 (2012).
  * [8] Bogdanov, A. N. & Rößler, U. K. Chiral symmetry breaking in magnetic thin films and multilayers. Phys. Rev. Lett. 87, 037203 (2001).
  * [9] Rößler, U. K., Bogdanov, A. N. & Pfleiderer, C. Spontaneous skyrmion ground states in magnetic metals. Nature 442, 797-801 (2006).
  * [10] Sampaio, J., Cros, V., Rohart, S., Thiaville, A. & Fert, A. Nucleation, stability and current-induced motion of isolated magnetic skyrmions in nanostructures. Nat. Nanotechnol. 8, 839–844 (2013).
  * [11] Li J. et al. Tailoring the topology of an artificial magnetic skyrmion. Nat. Commun. 5, 4704 (2014).
  * [12] Back, C. et al. The 2020 skyrmionics roadmap. J. Phys. D: Appl. Phys. 53, 363001 (2020).
  * [13] Jiang, W. et al. Blowing magnetic skyrmion bubbles. Science 349, 283-286 (2015).
  * [14] Du, H. et al. Edge-mediated skyrmion chain and its collective dynamics in a confined geometry. Nat. Commun. 6, 8504 (2015).
  * [15] Yuan, H. Y. & Wang, X. R. Creation and Manipulation by Nano-Second Current Pulses. Sci. Rep. 6, 22638 (2016).
  * [16] Uchida, M., Onose, Y., Matsui, Y. & Tokura, Y. Real-Space Observation of Helical Spin Order. Science 311, 359-361 (2006).
  * [17] Han, J. H., Zang, J. D., Yang, Z. H., Park, J. H. & Nagaosa, N. Skyrmion lattice in a two-dimensional chiral magnet. Phys. Rev. B 82, 094429 (2010).
  * [18] Butenko, A. B., Leonov, A. A. & Rößler, U. K. Stabilization of skyrmion textures by uniaxial distortions in noncentrosymmetric cubic helimagnets. Phys. Rev. B 82, 052403 (2010).
  * [19] Rohart, S. & Thiaville, A. Skyrmion confinement in ultrathin film nanostructures in the presence of Dzyaloshinskii-Moriya interaction. Phys. Rev. B 88, 184422 (2013).
  * [20] Leonov, A. O. et al. The properties of isolated chiral skyrmions in thin magnetic films. New J. Phys. 18, 065003 (2016).
  * [21] Bogdanov, A. & Hubert, A. Thermodynamically stable magnetic vortex states in magnetic crystals. J. Magn. Magn. Mater 138, 255-269 (1994).
  * [22] Jiang, W. et al. Nanoscale magnetic skyrmions and target states in confined geometries. Science 349, 283-286 (2015).
  * [23] Yu, G. et al. Room-temperature skyrmions in an antiferromagnet-based heterostructure. Nano. Lett. 18, 980-986 (2018).
  * [24] Birch, M. T. et al. Real-space imaging of confined magnetic skyrmion tubes. Nat. Commun. 11, 1-8 (2020).
  * [25] He, M. et al. Realization of zero-field skyrmions with high-density via electromagnetic manipulation in Pt/Co/Ta multilayers. Appl. Phys. Lett. 111, 202403 (2017).
  * [26] Jena, J. et al. Evolution and competition between chiral spin textures in nanostripes with D2d symmetry. Sci. Adv. 6, eabc0723 (2020).
  * [27] Müller, J. et al. Magnetic Skyrmions and Skyrmion Clusters in the Helical Phase of $\mathrm{Cu}_{2}\mathrm{OSeO}_{3}$ Phys. Rev. Lett. 119, 137201 (2017).
  * [28] Raju, M. et al. The evolution of skyrmions in Ir/Fe/Co/Pt multilayers and their topological Hall signature. Nat. Commun. 10, 1-7 (2019).
  * [29] Schoenherr, P. et al. Topological domain walls in helimagnets. Nat. Phys. 14, 465-468 (2018).
  * [30] Cortés-Ortun̈o, D. et al. Nanoscale magnetic skyrmions and target states in confined geometries. Phys. Rev. B 99, 214408 (2019).
  * [31] Ezawa, M. Compact merons and skyrmions in thin chiral magnetic films. Phys. Rev. B 83, 100408 (2011).
  * [32] Silva, R. L., Secchin, L. D., Moura-Melo, W. A., Pereira, A. R., & Stamps, R. L. Emergence of skyrmion lattices and bimerons in chiral magnetic thin films with nonmagnetic impurities. Phys. Rev. B 89, 054434 (2014).
  * [33] N. Nagaosa, and Y. Tokura, Topological properties and dynamics of magnetic skyrmions, Nat. Nanotech. 8, 899 (2013).
  * [34] Wang, X. S., Yuan, H. Y. & Wang, X. R. A theory on skyrmion size. Commun. Phys. 1, 31 (2018).
  * [35] Vansteenkiste, A. et al. The design and verification of MuMax3. AIP Adv. 4, 107133 (2014).
  * [36] Fert, A., Reyren, N. & Cros, V. Magnetic skyrmions: advances in physics and potential applications. Nat. Rev. Mater. 2, 17031 (2017).
  * [37] Zhang, X. et al. Skyrmion-electronics: writing, deleting, reading and processing magnetic skyrmions toward spintronic applications. J. Phys.: Condens. Matter 32, 143001 (2020).
  * [38] Sun, Z. Z. & Wang, X. R. Fast magnetization switching of Stoner particles: A nonlinear dynamics picture. Phys. Rev. B 71, 174430 (2005).
  * [39] Wang, X. R., Yan, P., Lu, J. & He, C. Magnetic field driven domain-wall propagation in magnetic nanowires. Ann. Phys. (NY) 324, 1815-1820 (2009).
  * [40] Wang, X. R., Yan, P. & Lu, J. High-field domain wall propagation velocity in magnetic nanowires. Europhys. Lett. 86, 67001 (2009).

Acknowledgement
This work is supported by Ministry of Science and Technology through grant
MOST20SC04-A, the NSFC Grant (No. 11974296 and 11774296) and Hong Kong RGC
Grants (No. 16301518 and 16301619). Partial support by the National Key
Research and Development Program of China (Grant No. 2018YFB0407600) is also
acknowledged.

Author contributions
X. R. Wang planned the project and wrote the manuscript. X.C.H. and H.T.W.
performed theoretical and numerical simulations, and prepared the figures. All
authors discussed the results and commented on the manuscript.

Additional information

Supplementary Information accompanies this paper at http://

Competing interests: The authors declare no competing interests.


# Making Human-Like Trade-offs in Constrained Environments by Learning from
Demonstrations

Arie Glazier1  Andrea Loreggia2  Nicholas Mattei1
Taher Rahgooy3  Francesca Rossi3  K. Brent Venable3,5

###### Abstract

Many real-life scenarios require humans to make difficult trade-offs: do we
always follow all the traffic rules or do we violate the speed limit in an
emergency? These scenarios force us to evaluate the trade-off between
collective norms and our own personal objectives. To create effective AI-human
teams, we must equip AI agents with a model of how humans make trade-offs in
complex, constrained environments. These agents will be able to mirror human
behavior or to draw human attention to situations where decision making could
be improved. To this end, we propose a novel inverse reinforcement learning
(IRL) method for learning implicit hard and soft constraints from
demonstrations, enabling agents to quickly adapt to new settings. In addition,
learning soft constraints over states, actions, and state features allows
agents to transfer this knowledge to new domains that share similar aspects.
We then use the constraint learning method to implement a novel system
architecture that leverages a cognitive model of human decision making, multi-
alternative decision field theory (MDFT), to orchestrate competing objectives.
We evaluate the resulting agent on trajectory length, number of violated
constraints, and total reward, demonstrating that our agent architecture is
both general and achieves strong performance. Thus we are able to capture and
replicate human-like trade-offs from demonstrations in environments when
constraints are not explicit.

##  1 Introduction

Implicit and explicit constraints are present in many decision making
scenarios, and they force us to make difficult decisions: do we always satisfy
all constraints, or do we violate some of them in exceptional circumstances?
Many techniques can be used to combine constraints and goals so that the agent
rationally minimizes constraint violations while achieving the given goal
(Noothigattu et al. 2019). However, it is well known that humans are not
rational. When we need to make a decision in a constrained environment, we
often reason by employing heuristics and approximations which are subject to
bias and noise (Busemeyer and Diederich 2002; Booch et al. 2021). This means
that optimal techniques may not be suitable if the aim is to design autonomous
artificial agents that act like humans, or decision support systems that
simulate human behavior to anticipate it and possibly alert humans by making
them aware of their reasoning and inference deficiencies.

Moreover, these constraints are often not explicitly given, but need to be
inferred from observations of how other agents act in the constrained world.
Learning constraints from demonstrations is an important topic in the domains
of inverse reinforcement learning (Scobee and Sastry 2020; Abbeel and Ng
2004), which is used to implement AI safety goals including value alignment
(Russell, Dewey, and Tegmark 2015) and to circumvent reward hacking (Amodei et
al. 2016; Ray, Achiam, and Amodei 2019). Recent work has focused on building
ethically bounded agents (Svegliato, Nashed, and Zilberstein 2021; Rossi and
Mattei 2019) that comply with ethical or moral theories of action. Following
the work of Scobee and Sastry (2020), we propose an architecture that, given
access to a model of the environment and to demonstrations of constrained
behavior, is able to learn constraints associated with states, actions, or
state features. Our method, MESC-IRL, performs comparably with the state of
the art and is more general, as it can handle both hard and soft constraints
in both deterministic and non-deterministic environments. It is also
decomposable into features of the environment, supporting the transfer of
learned constraints between environments.

Once the constraints are learned, we turn our attention to _making human like
trade-offs_. Enabling agents to make trade-offs like humans allows them to
mirror human behavior or to draw human attention to situations where decision
making could be improved (Booch et al. 2021). Additionally, making these
trade-offs explicit enables decision support tools that are able to mirror the
goals of human decision makers (Balakrishnan et al. 2018, 2019). We propose a
novel orchestration technique leveraging Multi-Alternative Decision Field
Theory (MDFT) (Busemeyer and Diederich 2002), a decision making framework
which is based on a psychological theory of how humans make decisions that is
able to capture deviations from rationality observed in humans, making trade-
offs between competing objectives in a more human-like way. We compare this
MDFT-based orchestrator with other methods both theoretically and empirically,
showing that our architecture is theoretically more expressive and obtains
better empirical performance across a range of metrics when acting in
constrained environments. The goal here is to use a cognitive model to capture
the sometimes irrational decisions made by humans. Building machines that
_act_ more like humans is a step to create effective human-machine teams or
decision support systems.

##  2 Preliminaries and Related Work

Figure 1: Grid with constraints of varying costs over actions (top right),
state occupancy (grid, left) and state features (bottom right). Note that
$\pi_{\mathcal{N}}^{*}$ for $\mathcal{M}^{\mathcal{N}}$ is the red trajectory,
obtained by an agent that does not know the constraints while
$\pi_{\mathcal{C}^{*}}^{*}$ for $\mathcal{M}^{\mathcal{C}^{*}}$ is the blue
trajectory.

We begin this section by providing the preliminary notions on the context of
our work, that is, constrained Markov Decision Processes and Reinforcement
Learning (Sutton and Barto 2018). We then review fundamental concepts and
methods on Inverse Reinforcement Learning (Ng and Russell 2000; Abbeel and Ng
2004) and background on Constrained Markov Decision Processes (Altman 1999)
including related work on learning constraints (Ziebart et al. 2008; Malik et
al. 2021; Scobee and Sastry 2020), which we will leverage to develop our novel
method for learning soft constraints (Rossi, Van Beek, and Walsh 2006) from
demonstrations (Chou, Berenson, and Ozay 2018). We conclude this section
presenting a short review of the Multi-Alternative Decision Field Theory
(Busemeyer and Diederich 2002), the cognitive model of decision making which
will be at the core of our novel approach to orchestrating competing
objectives.

###  2.1 Markov Decision Processes and Reinforcement Learning

A finite-horizon Markov Decision Process (MDP) $\mathcal{M}$ is a model for
sequential decision making over a number of time steps $t\in T$ defined by a
tuple $(\mathcal{S},\mathcal{A},P,D_{0},\phi,\gamma,R)$ (Sutton and Barto
2018). $\mathcal{S}$ is a finite set of discrete states;
$\\{\mathcal{A}_{s}\\}\subseteq\mathcal{A}$ is a set of actions available at
state $s$; $P:\mathcal{S}\times\mathcal{A}\times\mathcal{S}\rightarrow[0,1]$
is a model of the environment given as transition probabilities where
$P(s_{t+1}|s_{t},a_{t})$ is the probability of transitioning to state
$s_{t+1}$ from state $s_{t}$ after taking action
$a_{t}\in\\{\mathcal{A}_{s_{t}}\\}$ at time $t$.
$D_{0}:\mathcal{S}\rightarrow[0,1]$ is a distribution over start states;
$\phi:\mathcal{S}\times\mathcal{A}\times\mathcal{S}\rightarrow\mathbb{R}^{k}$
is a mapping from the transitions to a $k$-dimensional space of features;
$\gamma\in[0,1)$ is a discount factor; and
$R:\mathcal{S}\times\mathcal{A}\times\mathcal{S}\rightarrow\mathbb{R}$ is a
scalar reward received by the agent for being in one state and transitioning
to another state at time $t$, written as $R(s_{t},a_{t},s_{t+1})$.

An agent acts within the environment defined by the MDP, generating a sequence
of actions called a _trajectory_ of length $t$. Let
$\tau=((s_{1},a_{1},s_{2}),...,(s_{t-1},a_{t-1},s_{t}))\in(\mathcal{S}\times\mathcal{A}\times\mathcal{S})^{t}$.
We evaluate the quality of a particular trajectory in terms of the amount of
reward accrued over the trajectory, subject to discounting. Formally,
$R(\tau)=\sum^{t}_{i=1}\gamma^{i}R(s_{i},a_{i},s_{i+1})$. A policy,
$\pi:\mathcal{S}\rightarrow\mathcal{P}(\mathcal{A})$ is a map of probability
distribution to actions for every state such that $\pi(s,a)$ is the
probability of taking action $a$ in state $s$. We can also write the
probability of a trajectory $\tau$ under a policy as $\pi(\tau)$. The feature
vector associated with trajectory $\tau$ is defined as the summation over all
transition feature vectors in $\tau$,
$\phi(\tau)=\sum_{(s_{t},a_{t},s_{t+1})\in\tau}\phi(s_{t},a_{t},s_{t+1})$

The goal within an MDP is to find a policy $\pi^{*}$ that maximizes the
expected reward, $J(\pi)=\mathbb{E}_{\tau\sim\pi}[R(\tau)]$ (Malik et al.
2021). In the MDP literature, classical tabular methods are used to find
$\pi^{*}$ including value iteration (VI). Such method finds an optimal policy
by estimating the expected reward for taking an action $a$ in a given state
$s$, i.e., the $Q$-value of pair $(s,a)$, written $q(s,a)$. (Sutton and Barto
2018).

###  2.2 Constrained MDPs and Inverse Reinforcement Learning

We are interested in learning constraints from demonstrations. Our goal is to
create agents that are able to be trained to follow constraints that are not
explicitly prohibited in the MDP, but should be avoided (Rossi and Mattei
2019). (Scobee and Sastry 2020) discusses the importance of such constraints:
an MDP $\mathcal{M}$ may encode everything necessary about driving a car, e.g.
the dynamics of steering and movements, but often one wants to add additional
general constraints such as _avoid obstacles on the way to the goal_. These
constraints are often non-Markovian and engineering a reward function that
encodes these constraints may be a difficulty or impossible task (Vazquez-
Chanlatte et al. 2018).

One approach for learning constraints from demonstrations is to use techniques
from inverse reinforcement learning (IRL): given a set of demonstrated
trajectories $\mathcal{D}$ of an agent in an environment $\mathcal{M}$ with an
unknown reward function $\mathcal{M}\setminus R$, IRL provides a set of
techniques for learning a reward function $\hat{R}$ that explains the agent’s
demonstrated behavior (Abbeel and Ng 2004; Ng and Russell 2000). However, this
technique has many drawbacks: often there are many reward functions that lead
to the same behavior (Scobee and Sastry 2020), the reward functions may not be
interpretable (Vazquez-Chanlatte et al. 2018), and there may be issues such as
reward hacking – wherein the agent learns to behave in ways that create reward
but are not intended by the designer – an important topic in the field of AI
safety (Amodei et al. 2016; Ray, Achiam, and Amodei 2019) and value alignment
(Rossi and Mattei 2019; Russell, Dewey, and Tegmark 2015).

We follow the framework of Altman (1999) and Malik et al. (2021) and define a
Constrained MDP $\mathcal{M}^{\mathcal{C}}$ which is a nominal MDP
$\mathcal{M}^{\mathcal{N}}$ with an additional cost function
$\mathcal{C}:\mathcal{S}\times\mathcal{A}\times\mathcal{S}\rightarrow\mathbb{R}$
and a budget $\alpha\geq 0$. We can then define the cost of a trajectory to be
$c(\tau)=\sum^{t}_{i=1}c(s_{i},a_{i},s_{i+1})$. Setting $\alpha=0$ is
enforcing _hard constraints_ , i.e., we must never trigger constrained
transitions. In this work, unlike the work of both Scobee and Sastry (2020)
and Malik et al. (2021), we are interested in learning _soft constraints_
(Rossi, Van Beek, and Walsh 2006). Under a soft constraints paradigm, each
constraint comes with a real-valued penalty/cost and the goal is to minimize
the sum of penalties incurred by the agent.

Following Scobee and Sastry (2020), the task of constraint inference in IRL is
defined as follows. Given a nominal MDP $\mathcal{M}^{\mathcal{N}}$ and a set
of demonstrations $\mathcal{D}$ in ground-truth constrained world
$\mathcal{M}^{\mathcal{C}^{*}}$, we wish to find the most likely set of
constraints $\mathcal{C}$ that could modify $\mathcal{M}^{\mathcal{N}}$ to
explain the demonstrations. We are concerned with three types of constraints:

Action Constraints.

We may not want an agent to ever perform some (set of) action $a_{i}$.

Occupancy Constraints.

We may not want an agent to occupy a (set of) states $s_{i}$.

Feature Constraints.

Given a feature mapping of transitions $\phi$, we may not want an agent to
perform an (set of) action in presence of specific state features.

Without loss of generality, we add the state and actions to the features.
Hence, action and occupancy become specific cases of feature constraints. Note
that the set of constraints is defined as a cost function over the set of
transitions
$\mathcal{C}:\mathcal{S}\times\mathcal{A}\times\mathcal{S}\rightarrow\mathbb{R}$.
In Scobee and Sastry (2020), this definition is limited to a set of state-
actions $\mathcal{S}\times\mathcal{A}$ as they are assuming a deterministic
setting and hence are able to define $\mathcal{M}^{\mathcal{C}}$ by
substituting $\mathcal{A}=\\{A_{s}\\}$ with
$\mathcal{A}^{\mathcal{C}}=\\{A^{\mathcal{C}}_{s}\\}$ in
$\mathcal{M}^{\mathcal{N}}$. Finally, both Scobee and Sastry (2020) and Malik
et al. (2021) propose a greedy approach to infer a set of constraint
$\mathcal{C}$ that explains the demonstrations $D$ on
$\mathcal{M}^{\mathcal{C}^{*}}$. In both Scobee and Sastry (2020) and Malik et
al. (2021) the domain is restricted to deterministic MDPs, which we strictly
generalize in this work; additionally Scobee and Sastry (2020), like our
model, only works with discrete actions, while Malik et al. (2021) works for
both discrete and continuous action sets. We also generalize to the non-
deterministic setting; we additionally generalize to the setting of soft
constraints, hence our task is to learn the cost function $\mathcal{C}$.

To test our methods, we use the same grid world setup as Scobee and Sastry
(2020). Within our grid world example, shown in Figure 1, we have an action
penalty of $-4$ for the cardinal directions, $-4\times\sqrt{2}$ for taking the
diagonal actions, and reaching the goal state has a reward of $10$. In Figure
1 we set the constraint costs to various values but in all our experiments we
fix the constraint costs on the generated grids for states, actions, and
features to be $-50$. Throughout we assume a non-deterministic world with a
$10\%$ chance of action failure, resulting in a random action.

###  2.3 Multi-Alternative Decision Field Theory

Multi-alternative Decision Field Theory (MDFT) is a dynamic cognitive approach
that models human decision making based on psychological principles (Busemeyer
and Diederich 2002; Roe, Busemeyer, and Townsend 2001). MDFT models
preferential choice as an iterative cumulative process in which at each time
instant the decision maker attends to a specific attribute to derive
comparisons among options and update their preferences accordingly. Ultimately
the accumulation of those preferences informs the decision maker’s choice. In
MDFT an agent is confronted with multiple options and equipped with an initial
personal evaluation for them along different criteria, called attributes. For
example, a student who needs to choose a main course among those offered by
the cafeteria will have in mind an initial evaluation of the options in terms
of how tasty and healthy they look. More formally, an MDFT model is composed
of the following (Roe, Busemeyer, and Townsend 2001):

Personal Evaluation: Given set of options $O=\\{o_{1},\dots,o_{k}\\}$ and set
of attributes $A=\\{A_{1},\dots,A_{J}\\}$, the subjective value of option
$o_{i}$ on attribute $A_{j}$ is denoted by $m_{ij}$ and stored in matrix M. In
our example, let us assume that the cafeteria options are Salad (S), Burrito
(B) and Vegetable pasta (V). Matrix $\mathbf{M}$, containing the student’s
preferences, could be defined as shown in Figure 2 (left), where rows
correspond to the options $(S,B,V)$ and the columns to the attributes $Taste$
and $Health$.

Figure 2: Evaluation (M), Contrast (C), and Feedback (S) matrix.

Attention Weights: Attention weights are used to express the attention
allocated to each attribute at a particular time $t$ during the deliberation.
We denote them by vector $\textbf{W}(t)$ where $W_{j}(t)$ represents the
attention to attribute $j$ at time $t$. We adopt the common simplifying
assumption that, at each point in time, the decision maker attends to only one
attribute (Roe, Busemeyer, and Townsend 2001). Thus, $W_{j}(t)\in\\{0,1\\}$
and $\sum_{j}W_{j}(t)=1$, $\forall t,j$. In our example, we have two
attributes, so at any point in time $t$ we will have $\textbf{W}(t)=[1,0]$, or
$\textbf{W}(t)=[0,1]$, representing that the student is attending to,
respectively, $Taste$ or $Health$. The attention weights change across time
according to a stationary stochastic process with probability distribution w,
where $w_{j}$ is the probability of attending to attribute $A_{j}$. In our
example, defining $w_{1}=0.55$ and $w_{2}=0.45$ would mean that at each point
in time, the student will be attending $Taste$ with probability $0.55$ and
$Health$ with probability $0.45$. In other words, $Taste$ matters slightly
more than $Health$.

Contrast Matrix: Contrast matrix C is used to compute the advantage (or
disadvantage) of an option with respect to the other options. In the MDFT
literature (Busemeyer and Townsend 1993; Roe, Busemeyer, and Townsend 2001;
Hotaling, Busemeyer, and Li 2010), C is defined by contrasting the initial
evaluation of one alternative against the average of the evaluations of the
others, as shown for the case with three options in Figure 2 (center).

At any moment in time, each alternative in the choice set is associated with a
valence value. The valence for option $o_{i}$ at time $t$, denoted $v_{i}(t)$,
represents its momentary advantage (or disadvantage) when compared with other
options on some attribute under consideration. The valence vector for $k$
options $o_{1},\dots,o_{k}$ at time $t$, denoted by column vector
$\mathbf{V}(t)=[v_{1}(t),\dots,v_{k}(t)]^{T}$, is formed by
$\textbf{V}(t)=\textbf{C}\times\textbf{M}\times\textbf{W}(t)$. In our example,
the valence vector at any time point in which $\textbf{W}(t)=[1,0]$, is
$\textbf{V}(t)=[1-\nicefrac{{7}}{{2}},5-\nicefrac{{3}}{{2}},2-\nicefrac{{6}}{{2}}]^{T}$.

In MDFT, preferences for each option are accumulated across the iterations of
the deliberation process until a decision is made. This is done by using
Feedback Matrix $\mathbf{S}$, which defines how the accumulated preferences
affect the preferences computed at the next iteration. This interaction
depends on how similar the options are in terms of their initial evaluation
expressed in $\mathbf{M}$. Intuitively, the new preference of an option is
affected positively and strongly by the preference it had accumulated so far,
while it is inhibited by the preference of similar options. This lateral
inhibition decreases as the dissimilarity between options increases. Figure 2
(right) shows S for our example following the MDFT method in (Hotaling,
Busemeyer, and Li 2010).

At any moment in time, the preference of each alternative is calculated by
$\textbf{P}(t+1)=\textbf{S}\times\textbf{P}(t)+\textbf{V}(t+1)$ where
$\textbf{S}\times\textbf{P}(t)$ is the contribution of the past preferences
and $\textbf{V}(t+1)$ is the valence computed at that iteration. Starting with
$\textbf{P}(0)=0$, preferences are then accumulated for either a fixed number
of iterations (and the option with the highest preference is selected) or
until the preference of an option reaches a given threshold. In the first
case, MDFT models decision making with a specified deliberation time, while,
in the latter, it models cases where deliberation time is unspecified and
choice is dictated by the accumulated preference magnitude. In general,
different runs of the same MDFT model may return different choices due to the
attention weights’ distribution. In this way MDFT induces choice distributions
over set of options and is capable of capturing well know behavioral effects
such as the compromise, similarity, and attraction effects that have been
observed in humans and that violate rationality principles (Busemeyer and
Townsend 1993).

Figure 3: Comparision of MESC-IRL to the best performing method of Scobee and
Sastry (2020) at recovering hard constraints in a deterministic setting
according to false positive rate (left) and KL-Divergence from the
demonstrations $\mathcal{D}$ (right) as we vary the number of demonstrations.
Each point is the mean of 10 independent draws.

##  3 Learning Soft Constraints From Demonstrations

We now describe our method for learning a set soft constraints from a set of
demonstrations $\mathcal{D}$ and a nominal MDP $\mathcal{M}^{\mathcal{N}}$.
This is the first step in our goal of providing a flexible and expressive
method for orchestrating trade-offs. The method described here generalizes the
work of both Scobee and Sastry (2020) and Malik et al. (2021) to the setting
of non-deterministic MDPs and soft constraints.

###  3.1 MESC-IRL: Max Entropy Inverse Soft-Constraint Reinforcement Learning

Following Ziebart et al. (2008), our goal is to optimize a function that
linearly maps the features of each transition to the reward associated with
that transition, $R(s_{t},a_{t},s_{t+1})=\omega\phi(s_{t},a_{t},s_{t+1})$,
where $\omega$ is the reward weight vector. Ziebart et al. (2008) propose a
maximum entropy model for finding a unique solution ($\omega$) for this
problem. Based on this model, the probability of finite-length trajectory
$\tau$ being executed by an agent traversing an MDP $\mathcal{M}$ is
exponentially proportional to the reward earned by that trajectory and can be
approximated by:

$P(\tau|\omega)\approx\frac{e^{\omega^{T}\phi(\tau)}}{Z(\omega)}\prod_{(s_{t},a_{t},s_{t+1})\in\tau}P(s_{t+1}|s_{t},a_{t})$.

The optimal solution is obtained by finding the maximum likelihood of the
demonstrations $\mathcal{D}$ using this probability distribution:
$\omega^{*}=\underset{\omega}{\text{argmax}}\sum_{\tau\in\mathcal{D}}\log
P(\tau|\omega)$.

We extend the problem defined in Scobee and Sastry (2020) to that of learning
a set of _soft_ constraints which best explain a set of observed
demonstrations. This allows us to move from the notion of a constraint
forbidding an action or a state to that of a soft constraint imposing a
penalty proportional to the gravity of its violation. In other words, given
access to $\mathcal{M}^{\mathcal{N}}$ and a set of demonstrations
$\mathcal{D}$ in ground-truth constrained MDP $\mathcal{M}^{\mathcal{C}^{*}}$
we want to find the costs $\mathcal{C}$. More formally, we define the residual
reward function
$R^{\mathcal{R}}:\mathcal{S}\times\mathcal{A}\times\mathcal{S}\rightarrow\mathbb{R_{+}}$
as a mapping from the transitions to the penalties. We can now formally define
our soft-constrained MDP $\mathcal{M}^{\mathcal{C}}$ as follows:

###### Definition 1

Given
$\mathcal{M}^{\mathcal{N}}=\langle\mathcal{S},\mathcal{A},P,\mu,\phi,R^{\mathcal{N}}\rangle$
we define soft-constrained MDP
$\mathcal{M}^{\mathcal{C}}=\langle\mathcal{S},\mathcal{A},P,\mu,\phi,R^{\mathcal{C}}\rangle$
where $R^{\mathcal{C}}=R^{\mathcal{N}}-R^{\mathcal{R}}$.

Thus, the goal of our task is to find a residual reward function
$R^{\mathcal{R}}$ that maximizes the likelihood of the demonstrations
$\mathcal{D}$ given the nominal MDP $\mathcal{M}^{\mathcal{N}}$.

Our solution is based on adapting Maximum Causal Entropy Inverse Reinforcement
learning (Ziebart et al. 2008; Ziebart, Bagnell, and Dey 2010) to soft-
constrained MDPs. Following the setting of Ziebart et al. (2008) we can write
the reward function $R^{\mathcal{N}}$ (resp. $R^{\mathcal{C}}$) of
$\mathcal{M}^{\mathcal{N}}$ (resp. $\mathcal{M}^{\mathcal{C}}$) as a linear
combination of the transitions:
$R^{\mathcal{N}}(s_{t},a_{t},s_{t+1})=\omega^{\mathcal{N}}\phi(s_{t},a_{t},s_{t+1})$
and
$R^{\mathcal{C}}(s_{t},a_{t},s_{t+1})=\omega^{\mathcal{C}}\phi(s_{t},a_{t},s_{t+1})$.
As, both reward functions $R^{\mathcal{N}}$ and $R^{\mathcal{C}}$ are linear,
$R^{\mathcal{R}}$ should be linear as well
$R^{\mathcal{R}}=\omega^{\mathcal{R}}\phi(s_{t},a_{t},s_{t+1})$. From this
formulation of $R^{\mathcal{R}}$ we can infer that the reward vectors follow
$\omega^{\mathcal{C}}=\omega^{\mathcal{N}}-\omega^{\mathcal{R}}$.

At this point we can use Max Entropy IRL for learning a reward function
compatible with the trajectories in $\mathcal{D}$. The gradient for maximizing
the likelihood in this setting is defined as in Ziebart et al. (2008):

$\nabla_{\omega^{\mathcal{C}}}\mathcal{L}(\mathcal{D})=\mathbb{E}_{\mathcal{D}}[\phi(\tau)]-\sum_{(s_{t},a_{t},s_{t+1})}D_{s_{t},a_{t},s_{t+1}}\phi(s_{t},a_{t},s_{t+1})$

(1)

Where $D_{s_{t},a_{t},s_{t+1}}$ is the expected feature frequencies for
transition $(s_{t},a_{t},s_{t+1})$ using the current $\omega^{\mathcal{C}}$
weights. Given that the reward vectors follow
$\omega^{\mathcal{C}}=\omega^{\mathcal{N}}-\omega^{\mathcal{R}}$, we can write
$\nabla_{\omega^{\mathcal{C}}}=-\nabla_{\omega^{\mathcal{R}}}$. Finally, by
substituting this in Eq. 1 we obtain the gradient of likelihood of the
constrained trajectories with respect to $\omega^{\mathcal{R}}$:

$\nabla_{\omega^{\mathcal{R}}}\mathcal{L}(\mathcal{D})=\sum_{(s_{t},a_{t},s_{t+1})}D_{s_{t},a_{t},s_{t+1}}\phi(s_{t},a_{t},s_{t+1})-\mathbb{E}_{\mathcal{D}}[\phi(\tau)]$

(2)

As we estimate the residual rewards with respect to the nominal rewards, these
rewards are automatically scaled to be compatible with the nominal rewards.

###  3.2 Generalizing From Penalties to Probabilities

The estimated penalties from the previous section can effectively guide an
agent to navigate the environment optimally as well as provide estimates of
the cost of the constraints scaled to the value of the original reward signal.
However, there may be instances, such as when comparing with hard constraints,
where we desire _probabilities_ that a particular action is constrained.
Having probabilities allows us to compare constraints across environments with
possibly different scales, allows us to use this information to guide our
policies, and allows us to evaluate the confidence we have in a particular
constraint. In this section we describe a method to transition from penalties
to probabilities, as well as a generalized method to extract these
probabilities based on a subset of the features of the environment, which can
facilitate transfer learning between domains.

Intuitively, a transition where the residual reward, i.e., the penalty, is
significantly larger than zero is more likely to be a constraint. We estimate
the significance of a penalty by scaling it to the standard deviation of the
mean learned reward. Therefore, we assume that a transition penalty is a
random variable, denoted by $\mathbb{C}\sim
logistic(\sigma_{pooled},\sigma_{pooled})$, following a logistic distribution
with standard deviation $\sigma_{pooled}$, where
$\sigma_{pooled}=\sqrt{\nicefrac{{(\sigma_{\mathcal{N}}^{2}+\sigma_{\mathcal{C}}^{2})}}{{2}}}$
and $\sigma_{\mathcal{N}}$ and $\sigma_{\mathcal{C}}$ are the standard
deviations of the rewards in the nominal and learned constrained worlds,
respectively. Informally, when penalties are close to zero, we want their
probabilities to be small. To do this we set the mean of the distribution to
be $\mu=\sigma_{pooled}$.

We now want to reason about a random variable $\zeta$ that indicates our
belief that the transition $(s_{t},a_{t},s_{t+1})$ is forbidden. Hence using
the above probability distribution we can define the probability of constraint
given a transition as:

$\zeta\equiv P\left(\mathbb{C}\leq
R^{R}(s_{t},a_{t},s_{t+1})\right)=sigmoid\left(\frac{R^{R}(s_{t},a_{t},s_{t+1})-\sigma_{pooled}}{\sigma_{pooled}}\right)$

In our formulation, the residual rewards only depend on the features
associated with them. Hence, we can use this fact to reason about constraints
over only a subset of features $\mathbf{f}$, e.g., only color or state
position. Let $\phi_{f}\subseteq\phi$ be the subset of features we are
concerned with. In our grids we represent $\phi$ with a vector of length 92.
The first 81 elements represent the states, the next 8 represent the actions,
and the last 3 represent the colors. So if we are interested in only learning
about constraints over the colors, $\phi_{f}$ will be a vector equal to the
last three elements of $\phi$ that is $\phi_{color}\equiv\phi_{90,91,92}$.

Let $\phi_{\mathbf{f}}$ and $\omega^{R}_{\mathbf{f}}$ be the feature function
and residual feature weight vector for $\mathbf{f}$. We can now define the
probability of a feature value to be constrained as:

$\zeta_{\mathbf{f}}\equiv
P(\mathbb{C}\leq\omega^{R}_{\mathbf{f}}\phi_{\mathbf{f}})=sigmoid\left(\frac{\omega_{\mathbf{f}}\phi_{\mathbf{f}}-std_{pooled}}{std_{pooled}}\right)$.

###  3.3 Experimental Evaluation of MESC-IRL

In this section we empirically validate our method for soft constraint
learning against both the method of Scobee and Sastry (2020) for learning hard
constraints in deterministic settings as well as on learning soft constraints
in non-deterministic settings. Figure 3 shows the performance of MESC-IRL
compared to the method proposed by Scobee and Sastry (2020) on the same
metrics from their paper: false positives, i.e., predicting a constraint when
one does not exist, and KL-Divergence from the demonstrations set
$\mathcal{D}$. For this test we use the same single grid, hard constraints,
and a deterministic setting to allow for a direct comparison. We generate 10
independent sets of 100 demonstrations and report the mean. In order to decide
if the values returned by MESC-IRL represent a hard constraint, we threshold
the value of $\zeta$ at various levels and plot the comparison to the best
result from Scobee and Sastry (2020). MESC-IRL with $\zeta\geq 0.6$ performs
better than existing methods when the number of demonstrations is low, about
the same when there are more demonstrations, and is able to also work for soft
constraints and non-deterministic settings.

In order to evaluate MESC-IRL on soft constraints we need to adapt the notion
of false positives and false negatives. Let a false positive $fp$ be:

$fp=\frac{\Big{|}\big{\\{}x\in\mathcal{C}\mid
c(x)=0\wedge(\zeta_{\mathcal{C}}(x)-\zeta_{\mathcal{C}^{*}}(x)>\chi)\big{\\}}\Big{|}}{\text{Num.
Constraints}}$.

Where $\zeta_{\mathcal{C}}(x)$ and $\zeta_{\mathcal{C}^{*}}(x)$ are the
predicted and true probability of transition $x$ being constrained as
described in Section 3.2, and $\chi$ is a value in [0,1]. Intuitively, we
count a constraint as a false positive whenever there is no constraint in
$\mathcal{M}^{\mathcal{C}^{*}}$ and the predicted probability exceeds the true
probability by more than the threshold $\chi$. We can adapt the notion of
false negatives, $fn$ in the same way by taking $c(x)\neq 0$.

Figure 4: Performance of MESC-IRL on recovering soft constraints in non-
deterministic settings according to false negatives (left), false positives
(center), and JS-Divergence to $\mathcal{D}$ (right). We see that across all
these settings we are able to accurately recover constraints and generate
behavior similar to the $\mathcal{D}$ even with few demonstrations.

Figure 4 shows the results of our tests on recovering soft constraints in non-
deterministic settings with random grids, results for deterministic settings
can be found in the Appendix. For these tests we choose a start and a goal
state randomly at least 8 moves apart, set 6 states for blue, 6 for green
randomly, and select 6 randomly constrained states; all penalties are set to
$-50$. Again we take 10 sets of 100 demonstrations. We see a strong decrease
in both false positives and false negatives as the number of demonstrations
grows. We see that in general, and even more so when the optimal threshold
$\chi=0.2$ is selected, our method almost never adds constraints that are not
present in the ground truth and rarely underestimates the probability of
existing ones, even for relatively small demonstration sets. Likewise our
method is able to generate trajectories very close to $\mathcal{D}$, showing
that we are able to recover both constraints and behavior even with soft
constraints in non-deterministic setting. Hence MESC-IRL is able to work
across a variety of settings and accurately capture demonstrated constraints.

##  4 Orchestrating Goals and Constraints

Often humans are confronted with decisions that require making trade-offs
between collective norms and personal objectives (Rossi and Mattei 2019;
Noothigattu et al. 2019). In this section we investigate different ways to
model this orchestration in constrained grid environments. We consider
different methods for combining policies $\pi_{n}$ for the nominal
$\mathcal{M}^{\mathcal{N}}$ and $\pi_{c}$ for the learned constrained
$\mathcal{M}^{\mathcal{C}}$. For every state action pair $(s,a)$ we consider
vectors $\langle sq_{n}(s,a_{i})\rangle$ and $\langle sq_{c}(s,a_{i})\rangle$
with $i\in\\{1,\ldots,k\\}$ where $sq_{n}(s,a_{i})$ (resp. $sq_{c}(s,a_{i})$)
represents the probability of choosing action $a_{i}$ in state $s$ according
to policy $\pi_{n}$ (resp. $\pi_{c}$). They are obtained by taking the softmax
of the Q-values for each policy:

Greedy ($\pi_{G}$):

Let $\pi_{G}(s)=a$, where each $a$ is the one with highest Q-value,
$a$=$\underset{a\in\mathcal{A}_{s}}{argmax}$ $max\\{q_{c}(s,a),q_{n}(s,a)\\}$.

Weighted Average ($\pi_{WA}$):

Given weight vector $(w_{n},w_{c})$ with $w_{n}$,$w_{c}\in[0,1]$ and
$w_{c}+w_{n}=1$, action $a=\pi_{WA}(s)$ is chosen according to probability
distribution $p_{WA}(a_{i})=w_{n}sq_{n}(s,a_{i})+w_{c}sq_{c}(s,a_{i})$.

MDFT ($\pi_{MDFT}$):

Action $a=\pi_{MDFT}(s)$ is chosen via an MDFT model where: M is a $k\times 2$
matrix where rows (i.e., options) correspond to actions and columns (i.e.,
attributes) correspond to $\mathcal{M}^{\mathcal{N}}$ and
$\mathcal{M}^{\mathcal{C}}$. The $i$-th element of the respective world column
is $sq_{n}(s,a_{i})$ (resp., $sq_{c}(s,a_{i})$), i.e., we are using the
probability of choosing an action as a proxy of its preference. The weight
vector $(w_{n},w_{c})$ is defined as for $\pi_{WA}$, and serves as probability
distribution $\mathbf{w}$ defining how attention shifts between attributes
during deliberation. Matrices $\mathbf{C}$ and $\mathbf{S}$ are defined in the
standard way as described in Section 2.3. When reaching state $s$, an MDFT
deliberation process is launched to decide which action should be chosen. At
each step the focus is shifted to $\mathcal{M}^{\mathcal{N}}$ or
$\mathcal{M}^{\mathcal{C}}$ according to probability distribution
$(w_{n},w_{c})$, and the preferences of the actions according to the selected
attribute are accumulated as per Section 2.3.

Informally, Greedy is a deterministic approach that takes the most promising
action, WA allows the agent to prioritize the pursuit of the goal state and
satisfying constraints via a new policy obtained by considering the weighted
average of the nominal and constrained distributions, and the MDFT-based
orchestrator uses MDFT to chose at each step an action as suggested by the
MDFT machinery.

Figure 5: Comparison of the three orchestrators (Greedy, WA, and MDFT) on
metrics including average path length (left), normalized penalty (lower is
better, center), and average number of violated constraints (lower is better,
right). We see that across all the metrics the MDFT orchestration is able to
generate shorter paths that incur less penalty and violate fewer constraints.

###  4.1 Comparison of Orchestration Methods

We first compare theoretically the expressive power of the three
orchestrators. We focus on a single state $s$ and consider how the policies
compare in terms of being able to model a given distribution over the actions
available in $s$. We start by considering the Greedy orchestrator that is
deterministic and will pick a fixed action $a$ in state $s$. Both WA and MDFT
can model the Greedy policy by shifting all the weight to the environment
where the maximum value is obtained and zeroing all preferences except for
that of action $a$. A formal description is provided in the Appendix. This
observation, along with the fact that MDFT and WA are non-determistic, allows
us to conclude that Greedy is strictly less expressive than the other two
orchestrators.

Turning to the comparison between MDFT and WA, we can prove the following
statement.

###### Theorem 1

Given any state $s$, there exist choice probability distributions over the
actions available in $s$ that can be modeled by MDFT but not by WA.

We use an instance of the well known compromise effect (Busemeyer and
Diederich 2002) according to which a compromising alternative tends to be
chosen more often by humans than options with complementary preferences with
respect to the attributes. Consider the case of state $s$ with three actions
$a_{1}$, $a_{2}$ and $a_{3}$. Let us assume that, for example,
$sq_{n}(s,a_{1})=\nicefrac{{1}}{{6}}$, $sq_{n}(s,a_{2})=\nicefrac{{1}}{{3}}$
$sq_{n}(s,a_{3})=\nicefrac{{1}}{{2}}$ and
$sq_{c}(s,a_{1})=\nicefrac{{1}}{{2}}$, $sq_{n}(s,a_{2})=\nicefrac{{1}}{{3}}$
$sq_{n}(s,a_{3})=\nicefrac{{1}}{{6}}$. According to the compromise effect
humans will tend to choose $a_{2}$ more often than $a_{1}$ and $a_{3}$. Such a
choice distribution over the actions can be modeled by an MDFT defined over
option set $\\{a_{1},a_{2},a_{3}\\}$, with two attributes and weights
$w_{n}=0.55$ and $w_{c}=0.45$ (Busemeyer and Diederich 2002). However, if we
now consider WA, we can see that there is no way to define weights
$(w_{n},w_{c})$ such that the corresponding weighted average probability
satisfies
$w_{n}sq_{n}(s,a_{2})+w_{c}sq_{c}(s,a_{2})>max\\{w_{n}sq_{n}(s,a_{1})+w_{c}sq_{c}(s,a_{1}),w_{n}sq_{n}(s,a_{3})+w_{c}sq_{c}(s,a_{3})\\}$.
Thus, this distribution over actions cannot be modeled by the WA.

On the other hand, if we consider MDFTs in general, i.e. without the
restriction of having two attributes, we can model any distribution.
Intuitively this is achieved by defining an MDFT model over $k$ actions and
with $k$ attributes where the weight of the $i$-th attribute corresponds to
the probability of the $i$-th action. Matrix $\mathbf{M}$ is set to the
identity matrix and deliberation is halted after one iteration; see the
Appendix for details.

###### Theorem 2

Given a $s$ and the set $\mathcal{A}_{s}$ of actions available in $s$,
consider a probability distribution $p$ defined over $\mathcal{A}_{s}$. We can
define an MDFT model where the set of options corresponds to $\mathcal{A}_{s}$
and the induced choice probability distribution coincides with $p$.

As a consequence, MDFT is general enough to express the probability
distributions induced over the actions by WA. Whether this is true also in the
case of MDFT with only two attributes, as used in $\pi_{MDFT}$, remains an
open theoretical question. However, we can see this experimentally: in Rahgooy
and Venable (2019) the authors propose an RNN-based approach that starts from
samples of a choice distribution and recovers parameters of an MDFT model,
minimizing the divergence between the original and MDFT-induced choice
distributions. We adapt their code111Available at
https://github.com/Rahgooy/MDFT and generate 100 instances of WA distributions
starting from random $sq_{n}$ and $sq_{c}$ distributions and $(w_{n},w_{c})$
weights. For each of these instances we generate 100 samples (i.e, chosen
actions). We fix the $sq_{n}$ and $sq_{c}$ values as parameters for the
$\mathbf{M}$ matrix and learn the attention weight distribution $\mathbf{w}$
using 300 learning iterations. We use the learned MDFT model to generate a
choice distribution over the actions with a stopping criteria of 25
deliberations steps. The observed average JS divergence between the original
WA distributions and the ones induced by learned MDFT is 0.024 with standard
error 0.0013; showing experimentally we can learn weights for an MDFT model to
replicate any choice distribution of WA.

###  4.2 Experimentally Evaluating Orchestrators

We compare the orchestrators empirically with the goal of testing if the
combination of MESC-IRL with the orchestration techniques can be leveraged to
create agents that trade-off between conflicting objectives like humans.

We start by generating 100 different non-deterministic nominal worlds,
$\mathcal{M}^{\mathcal{N}}$, as described in Sections 2 and 3. We learn, via
VI on $\mathcal{M}^{\mathcal{C}^{*}}$ the optimal policy in the (ground truth)
constrained world, denoted $\pi_{\mathcal{C}^{*}}^{*}$ and, similarly to
Scobee and Sastry (2020), we use it to generate sets of 200 demonstrations,
$\mathcal{D}$. We then pass $\mathcal{D}$ to MESC-IRL, and the learned
constraints are added to $\mathcal{M}^{\mathcal{N}}$ yielding the learned
constrained MPD, $\mathcal{M}^{\mathcal{C}}$. We use VI on
$\mathcal{M}^{\mathcal{N}}$and $\mathcal{M}^{\mathcal{C}}$to obtain $\pi_{n}$
and $\pi_{c}$, and we consider different ways to prioritize them by sweeping
the weight values $(w_{n},w_{c})$ from $(0,1)$ to $(1,0)$ in steps of $0.1$.
Note that at $(1,0)$ (resp. $(0,1)$) WA is equivalent to $\pi_{n}$ (resp.
$\pi_{c}$), and that in both cases MDFT becomes deterministic, picking the
action with highest Q-value.

For all our results we generated 200 trajectories for each step and method
(including $\pi_{n}$ and $\pi_{c}$, denoted as Nominal and Constrained in Fig.
5), and for each of the 100 random worlds. We first perform the Kolmogorov-
Smirnov test to see if the trajectories generated by WA and MDFT induce the
same distribution ($\text{H}_{0}$). We reject $\text{H}_{0}$ at every weight
step with $p\\_value\leq 0.01$, thus the two techniques induce statistically
significantly different choice distributions.

Figure 6: Average JS Divergence between policies generated with $\pi_{c}$ and
orchestrators MDFT and WA.

Figure 5 (left) shows the average length of trajectories produced by the
orchestrators normalized so that $1.0$ is the shortest path between the start
and goal state; (center) we scale the penalty by the average penalty for
trajectories in $\mathcal{M}^{\mathcal{C}}$, lower is better; (right) we show
the average number of violated constraints. Across all these metrics, the MDFT
agent is performing better than WA by always reaching the goal in a smaller
number of steps no matter the configuration of the orchestrator. We can also
see that the MDFT agent violates fewer constraints and accumulates lower
penalties.

Finally, in Figure 6 we show the JS Divergence between the trajectories
generated by $\pi_{\mathcal{C}^{*}}^{*}$ and the trajectories generated by
MDFT and WA, as the weight vector varies. For both agents, the divergence is
small on the left and grows moving to the right, as constraints become less
important. This is not surprising, since the reference trajectories are
generated using $\pi_{\mathcal{C}^{*}}^{*}$. Furthermore, we note that the
MDFT advantage is more significant when $w_{c}$ is larger, that is when
constraints matter more. An explanation for this is that a large value $w_{c}$
results in more MDFT deliberation steps to be focused (exclusively) on
preferences relative to the constrained world. In WA, the averaging of the
values underlying the policies, although weighted, is not able to maintain the
importance of the constraints.

##  5 Conclusions and Future Work

We proposed a novel and general constraint learning method combined with a
unique agent architecture aimed at learning constraints from demonstrations
and exhibiting human-like trade-offs in environments with competing
objectives. Our theoretical and experimental results show that the MDFT-based
method exhibits superior expressive power and performance both in terms of the
quality of the produced trajectories as well as capability of capturing
initial demonstrations. Another important features of this cognitive approach
is its effectiveness in capturing behavioral traits of humans decision making,
a key factor for real life applications involving human-generated
demonstrations and decisions. Our MESC-IRL approach is a general approach for
learning soft constraints over actions, states, and features in non-
deterministic decision making environments.

We plan to run experiments with human decision makers and employ also
learning-based orchestrators Noothigattu et al. (2019). We are also working on
a novel multi-agent architecture with several orchestrators acting as either
reactive or deliberate agents, e.g., the system 1 / system 2 model of Kahneman
(2011) with a meta-cognitive agent to arbitrate, with the goal of further
advancing the performance and generality of decision making agents.

## References

  * Abbeel and Ng (2004) Abbeel, P.; and Ng, A. Y. 2004.  Apprenticeship learning via inverse reinforcement learning.  In _Proceedings of the 21st International Conference on Machine Learning (ICML)_.
  * Altman (1999) Altman, E. 1999.  _Constrained Markov Decision Processes_ , volume 7.  CRC Press.
  * Amodei et al. (2016) Amodei, D.; Olah, C.; Steinhardt, J.; Christiano, P.; Schulman, J.; and Mané, D. 2016.  Concrete problems in AI safety.  _arXiv preprint arXiv:1606.06565_.
  * Balakrishnan et al. (2018) Balakrishnan, A.; Bouneffouf, D.; Mattei, N.; and Rossi, F. 2018.  Using Contextual Bandits with Behavioral Constraints for Constrained Online Movie Recommendation.  In _Proc. of the 27th Intl. Joint Conference on AI (IJCAI)_.
  * Balakrishnan et al. (2019) Balakrishnan, A.; Bouneffouf, D.; Mattei, N.; and Rossi, F. 2019.  Incorporating Behavioral Constraints in Online AI Systems.  In _Proc. of the 33rd AAAI Conference on Artificial Intelligence (AAAI)_.
  * Booch et al. (2021) Booch, G.; Fabiano, F.; Horesh, L.; Kate, K.; Lenchner, J.; Linck, N.; Loreggia, A.; Murugesan, K.; Mattei, N.; Rossi, F.; and Srivastava, B. 2021.  Thinking Fast and Slow in AI.  In _Thirty-Fifth AAAI Conference on Artificial Intelligence, AAAI 2021, Thirty-Third Conference on Innovative Applications of Artificial Intelligence, IAAI 2021, The Eleventh Symposium on Educational Advances in Artificial Intelligence, EAAI 2021, Virtual Event, February 2-9, 2021_ , 15042–15046. AAAI Press.
  * Busemeyer and Diederich (2002) Busemeyer, J. R.; and Diederich, A. 2002.  Survey of decision field theory.  _Mathematical Social Sciences_ , 43(3): 345–370.
  * Busemeyer and Townsend (1993) Busemeyer, J. R.; and Townsend, J. T. 1993.  Decision field theory: a dynamic-cognitive approach to decision making in an uncertain environment.  _Psychological review_ , 100(3): 432.
  * Chou, Berenson, and Ozay (2018) Chou, G.; Berenson, D.; and Ozay, N. 2018.  Learning constraints from demonstrations.  _arXiv preprint arXiv:1812.07084_.
  * Hotaling, Busemeyer, and Li (2010) Hotaling, J. M.; Busemeyer, J. R.; and Li, J. 2010.  Theoretical developments in decision field theory: Comment on Tsetsos, Usher, and Chater (2010).  _Psychological Review_ , 117(4).
  * Kahneman (2011) Kahneman, D. 2011.  _Thinking, Fast and Slow_.  Macmillan.
  * Malik et al. (2021) Malik, S.; Anwar, U.; Aghasi, A.; and Ahmed, A. 2021.  Inverse Constrained Reinforcement Learning.  In Meila, M.; and Zhang, T., eds., _Proceedings of the 38th International Conference on Machine Learning_ , volume 139 of _Proceedings of Machine Learning Research_ , 7390–7399. PMLR.
  * Ng and Russell (2000) Ng, A. Y.; and Russell, S. J. 2000.  Algorithms for Inverse Reinforcement Learning.  In _Proceedings of the Seventeenth International Conference on Machine Learning_ , ICML ’00, 663–670. San Francisco, CA, USA: Morgan Kaufmann Publishers Inc.  ISBN 1-55860-707-2.
  * Noothigattu et al. (2019) Noothigattu, R.; Bouneffouf, D.; Mattei, N.; Chandra, R.; Madan, P.; Varshney, K. R.; Campbell, M.; Singh, M.; and Rossi, F. 2019.  Teaching AI agents ethical values using reinforcement learning and policy orchestration.  _IBM J. Res. Dev._ , 63(4/5): 2:1–2:9.
  * Rahgooy and Venable (2019) Rahgooy, T.; and Venable, K. B. 2019.  Learning Preferences in a Cognitive Decision Model.  In Zeng, A.; Pan, D.; Hao, T.; Zhang, D.; Shi, Y.; and Song, X., eds., _Human Brain and Artificial Intelligence_ , 181–194. Singapore: Springer Singapore.  ISBN 978-981-15-1398-5.
  * Ray, Achiam, and Amodei (2019) Ray, A.; Achiam, J.; and Amodei, D. 2019.  Benchmarking safe exploration in deep reinforcement learning.  _arXiv preprint arXiv:1910.01708_ , 7.
  * Roe, Busemeyer, and Townsend (2001) Roe, R. M.; Busemeyer, J. R.; and Townsend, J. T. 2001.  Multialternative decision field theory: A dynamic connectionst model of decision making.  _Psychological review_ , 108(2): 370.
  * Rossi and Mattei (2019) Rossi, F.; and Mattei, N. 2019.  Building Ethically Bounded AI.  In _Proc. of the 33rd AAAI Conference on Artificial Intelligence (AAAI)_.
  * Rossi, Van Beek, and Walsh (2006) Rossi, F.; Van Beek, P.; and Walsh, T. 2006.  _Handbook of Constraint Programming_.  Elsevier.
  * Russell, Dewey, and Tegmark (2015) Russell, S.; Dewey, D.; and Tegmark, M. 2015.  Research priorities for robust and beneficial artificial intelligence.  _AI Magazine_ , 36(4): 105–114.
  * Scobee and Sastry (2020) Scobee, D. R. R.; and Sastry, S. S. 2020.  Maximum Likelihood Constraint Inference for Inverse Reinforcement Learning.  In _8th International Conference on Learning Representations, ICLR 2020, Addis Ababa, Ethiopia, April 26-30, 2020_. OpenReview.net.
  * Sutton and Barto (2018) Sutton, R. S.; and Barto, A. G. 2018.  _Reinforcement Learning: An Introduction, 2nd Edition_.  Cambridge, MA, USA: A Bradford Book.
  * Svegliato, Nashed, and Zilberstein (2021) Svegliato, J.; Nashed, S. B.; and Zilberstein, S. 2021.  Ethically compliant sequential decision making.  In _Proceedings of the 35th AAAI International Conference on Artificial Intelligence (AAAI)_.
  * Vazquez-Chanlatte et al. (2018) Vazquez-Chanlatte, M.; Jha, S.; Tiwari, A.; Ho, M. K.; and Seshia, S. A. 2018\.  Learning Task Specifications from Demonstrations.  In Bengio, S.; Wallach, H. M.; Larochelle, H.; Grauman, K.; Cesa-Bianchi, N.; and Garnett, R., eds., _Advances in Neural Information Processing Systems 31: Annual Conference on Neural Information Processing Systems 2018, NeurIPS 2018, December 3-8, 2018, Montréal, Canada_ , 5372–5382.
  * Ziebart, Bagnell, and Dey (2010) Ziebart, B. D.; Bagnell, J. A.; and Dey, A. K. 2010.  Modeling Interaction via the Principle of Maximum Causal Entropy.  In Fürnkranz, J.; and Joachims, T., eds., _Proceedings of the 27th International Conference on Machine Learning (ICML-10), June 21-24, 2010, Haifa, Israel_ , 1255–1262. Omnipress.
  * Ziebart et al. (2008) Ziebart, B. D.; Maas, A. L.; Bagnell, J. A.; and Dey, A. K. 2008.  Maximum Entropy Inverse Reinforcement Learning.  In Fox, D.; and Gomes, C. P., eds., _Proceedings of the Twenty-Third AAAI Conference on Artificial Intelligence, AAAI 2008, Chicago, Illinois, USA, July 13-17, 2008_ , 1433–1438. AAAI Press.

## Additional Material for
Making Human-Like Trade-offs in Constrained Environments by Learning from
Demonstrations

##  Appendix A Additional Graphs For Learning Constraints

Additional graphs and results for comparision with the methods of Scobee and
Sastry (2020) and MESC-IRL. Figure 7 shows the stepping results for MESC-IRL.

Figure 7: Performance of MESC-IRL for various settings of $\zeta$ at
recovering hard constraints in a deterministic setting according to false
positive rate (left) and KL-Divergence from the demonstrations $\mathcal{D}$
(right) as we vary the number of demonstrations. Each point is the mean of 10
independent draws.

Figure 8 shows the performance of our best cutoff with the best method from
Scobee and Sastry (2020).

Figure 8: Comparision of the best performing version of MESC-IRL to the best
performing method of Scobee and Sastry (2020) at recovering hard constraints
in a deterministic setting according to false positive rate (left) and KL-
Divergence from the demonstrations $\mathcal{D}$ (right) as we vary the number
of demonstrations. Each point is the mean of 10 independent draws.

Figure 9 shows the performance of MESC-IRL on recovering soft constraints in
the deterministic setting.

Figure 9: Performance of MESC-IRL on recovering soft constraints in
deterministic settings according to false negatives (left), false positives
(center), and JS-Divergence to $\mathcal{D}$ (right). We see that across all
these settings we are able to accurately recover constraints and generate
behavior similar to the $\mathcal{D}$ even with few demonstrations.

Figure 10 shows the performance of MESC-IRL on recovering soft constraints in
the non-deterministic setting.

Figure 10: Performance of MESC-IRL on recovering soft constraints in
deterministic settings according to false negatives (left), false positives
(center), and JS-Divergence to $\mathcal{D}$ (right). We see that across all
these settings we are able to accurately recover constraints and generate
behavior similar to the $\mathcal{D}$ even with few demonstrations.

##  Appendix B Proof Details for Comparison of Orchestration Methods

We provide the full proofs for the theoretical comparison of the
orchestrators.

###### Theorem 3

Consider state $s$. Any choice probability distribution over the actions
available in $s$ that can be modeled by the Greedy approach can be modeled via
the MDFT or WA approaches.

Proof. We can model the (degenerate) probability distribution induced by
Greedy via an MDFT with as many options as the actions available in $s$, two
attributes with weights set to any random pair of values, and preferences in
the $\mathbf{M}$ matrix all equally to 0 except for those in the row
associated with $a$ which are set to 1. Matrices $\mathbf{C}$ and $\mathbf{S}$
can be defined in the standard way described in Section 2.3 and deliberation
can be halted after one deliberation step. In fact, when deliberation is
launched, an attribute will be selected. Regardless of which one is selected,
action $a$ will be chosen given that it is the only one with non-zero
preference.

Similarly, we can model the Greedy distribution using a weighted average where
$w_{n}=w_{c}=1/2$, and $sq_{n}(s,a)=sq_{c}(s,a)=1$ and
$sq_{n}(s,a^{\prime})=sq_{c}(s,a^{\prime})=0$, $\forall a^{\prime}\neq a$.
$\Box$

###### Theorem 4

Given state $s$ and the set $\mathcal{A}_{s}$ of actions available in $s$,
consider a probability distribution $p$ defined over $\mathcal{A}_{s}$. We can
define an MDFT model where the set of options corresponds to $\mathcal{A}_{s}$
and the induced choice probability distribution coincides with $p$.

Prrof. Consider the MDFT model defined as follows:

  * •

Matrix $\mathbf{M}$ is the $k\times k$ identity matrix;

  * •

Weight vectors $\mathbf{W}$ are defined as in Section 2.3 and select a single
attribute at each iteration. Probability distribution over attributes
$\mathbf{w}$ is defined in a way such that the probability of selecting the
$j$-th attribute, is $p(a_{j})$.

  * •

Matrices $\mathbf{C}$ and $\mathbf{S}$ are defined in the standard way as
described in Section 2.3.

  * •

The deliberation time for the model is fixed at one iteration.

It is easy to see that running the model induces a choice probability over the
actions which corresponds to $p$. In fact, in every run, which consists of a
single iteration, an attribute $A_{h}$ will be sampled according to
probability $p$. Given how $\mathbf{M}$ is defined and the fact that the
initial value of the accumulated preference $\mathbf{P}(0)=0$, action $a_{h}$
will be chosen. Thus, the probability of action $a_{h}$ being selected, given
the MDFT model, coincides with $p(a_{h})$. $\Box$