# Collective fast neutrino flavor conversions in a 1D box: Initial conditions
and long-term evolution

Meng-Ru Wu<EMAIL_ADDRESS>Institute of Physics, Academia Sinica,
Taipei, 11529, Taiwan  Institute of Astronomy and Astrophysics, Academia
Sinica, Taipei, 10617, Taiwan  Manu George<EMAIL_ADDRESS>Institute
of Physics, Academia Sinica, Taipei, 11529, Taiwan  Chun-Yu Lin
<EMAIL_ADDRESS>National Center for High-performance computing, National
Applied Research Laboratories, Hsinchu Science Park, Hsinchu City 30076,
Taiwan  Zewei Xiong<EMAIL_ADDRESS>GSI Helmholtzzentrum für
Schwerionenforschung, Planckstraße 1, 64291 Darmstadt, Germany

###### Abstract

We perform numerical simulations of fast collective neutrino flavor
conversions in an one-dimensional box mimicking a system with the periodic
boundary condition in one spatial direction and translation symmetry in the
other two dimensions. We evolve the system over several thousands of the
characteristic timescale (inverse of the interaction strength) with different
initial $\bar{\nu}_{e}$ to $\nu_{e}$ number density ratios and different
initial seed perturbations. We find that small scale structures are formed due
to the interaction of the flavor waves. This results in nearly flavor
depolarization in a certain neutrino phase space, when averaged over the
entire box. Specifically, systems with initially equal number of $\nu_{e}$ and
$\bar{\nu}_{e}$ can reach full flavor depolarization for the entire neutrino
electron lepton number ($\nu$ELN) angular spectra. For systems with initially
unequal $\nu_{e}$ and $\bar{\nu}_{e}$, flavor depolarization can only be
reached in one side of the $\nu$ELN spectra, dictated by the net neutrino
$e-x$ lepton number conservation. Quantitatively small differences depending
on the initial perturbations are also found when different perturbation seeds
are applied. Our numerical study here provides new insights for efforts aiming
to include impact of fast flavor conversions in astrophysical simulations
while calls for better analytical understanding accounting for the evolution
of fast flavor conversions.

##  I Introduction

The discovery of the flavor oscillations of neutrinos by terrestrial
experiments and astrophysical observations is one of the most exciting
milestones in neutrino physics. Ongoing and planned experimental projects are
expected to further pin down the yet-unknown parameters in neutrino mixing –
the mass ordering and the _CP_ -violating phase – while searching for
potential signatures of physics beyond the Standard Model Zyla  _et al._
(2020).

However, despite the success of the theory of neutrino mixing in explaining
the majority of experimental data, one aspect that is poorly understood yet is
how neutrinos oscillate in astrophysical environments dense in neutrinos.
Analytical and numerical works over the past decades have shown that in an
environment where the self-interactions between neutrinos cannot be ignored,
various collective phenomena can arise due to the nonlinear and strong
coupling nature (in the flavor space) of the system; see e.g., Duan  _et al._
(2006a, b); Hannestad  _et al._ (2006); Esteban-Pretel  _et al._ (2008);
Dasgupta  _et al._ (2009); Malkus  _et al._ (2012); Cherry  _et al._ (2012);
Raffelt  _et al._ (2013); Vlasenko  _et al._ (2014); Volpe  _et al._ (2013);
Wu  _et al._ (2016); Abbar and Duan (2015); Sawyer (2016); Izaguirre  _et al._
(2017); Capozzi  _et al._ (2019); Abbar (2021); Xiong and Qian (2021); Johns
(2021) and review articles Duan  _et al._ (2010); Mirizzi  _et al._ (2015);
Duan (2015); Tamborra and Shalgar (2020). Exploratory works also demonstrated
that such collective neutrino flavor oscillations may largely affect our
understanding of important astrophysical events such as the core-collapse
supernovae and the merger of two neutron stars Duan  _et al._ (2011); Wu  _et
al._ (2015); Sasaki  _et al._ (2017); Wu  _et al._ (2017); Wu and Tamborra
(2017); Stapleford  _et al._ (2020); Xiong  _et al._ (2020); George  _et al._
(2020); Li and Siegel (2021).

Among these efforts, an important aspect that was recently pointed out is the
potential occurrence of the “fast” neutrino flavor conversions Sawyer (2016).
Fast flavor conversions happen when the angular distribution of the neutrino
electron lepton numbers ($\nu$ELN) take both positive and negative values –
dubbed “crossings”. Under the two-flavor approximation in $\nu_{e}$–$\nu_{x}$
subspace where $\nu_{x}$ refers to a linear combination of $\nu_{\mu}$ and
$\nu_{\tau}$, a $\nu$ELN crossing means that the effective differential
neutrino $e-x$ number density $d{\tilde{n}}_{\nu}/d\Omega$, where
$\tilde{n}_{\nu}=n_{\nu_{e}}-n_{\bar{\nu}_{e}}-n_{\nu_{x}}+n_{\bar{\nu}_{x}}$
and $\Omega$ is the solid angle, transits from positive to negative (or vice
versa) in the angular phase space. This has been confirmed by means of linear
instability analysis, valid when the flavor conversion probabilities remain
small, as well as numerical studies Izaguirre  _et al._ (2017); Dasgupta  _et
al._ (2017); Capozzi  _et al._ (2017); Abbar and Volpe (2019); Yi  _et al._
(2019); Martin  _et al._ (2020); Padilla-Gay  _et al._ (2021); Johns  _et al._
(2020a); Bhattacharyya and Dasgupta (2020, 2021); Morinaga (2021); Richers
_et al._ (2021); Zaizen and Morinaga (2021); Kato  _et al._ (2021). Meanwhile,
searches for the conditions where fast conversions can develop using neutrino
angular distributions provided by hydrodynamical simulations of supernovae and
mergers were carried out extensively Morinaga  _et al._ (2020); Nagakura  _et
al._ (2019); Johns  _et al._ (2020b); Delfan Azari  _et al._ (2020); Abbar
_et al._ (2020); Glas  _et al._ (2020); Abbar (2020); Johns and Nagakura
(2021); Nagakura  _et al._ (2021).

In particular, Refs. Martin  _et al._ (2020) and Bhattacharyya and Dasgupta
(2021) examined how fast flavor conversions develop in the nonlinear regime
using numerical simulations in a one-dimensional (1D) box which possesses
translation symmetry in the $x$ and $y$ directions while periodic in the $z$
direction. Reference Martin  _et al._ (2020) showed that coherent and wavelike
patterns in flavor space can develop with a point-source like perturbation. On
the other hand, Ref. Bhattacharyya and Dasgupta (2021) found that the system
can quickly settle into a state where flavor depolarization happens when
averaged over the $z$ domain, for a major part of neutrino spectrum with
either point-source like or random perturbation seeds. These findings seem to
contradict each other and require further examinations.

In this work, we systematically investigate the dependence of the outcome of
fast neutrino flavor conversions in an 1D box on the initial condition of the
perturbation seeds. By adopting advanced numerical schemes, we evolve the
systems to late times when quasi steady states are achieved. We also explore
the dependence of the steady-state outcome on the initial $\bar{\nu}_{e}$ to
$\nu_{e}$ number density ratio.

In Sec. II, we describe our models, the initial and boundary conditions, and
the adopted numerical schemes. We also briefly review the notation of the
neutrino polarization vectors. In Sec. IV, we discuss our numerical simulation
results and the implications. We summarize our main findings in Sec. V.
Throughout the paper, we adopt natural units and take $\hbar=c=1$.

##  II Models

###  II.1 Equation of motion

We consider a simple neutrino-dense system which has translation symmetry in
the $x$\- and $y$\- directions in space, the same as in Refs. Martin  _et al._
(2020); Bhattacharyya and Dasgupta (2020, 2021). For simplicity, initially we
assume a pure system consisting of only $\nu_{e}$ and $\bar{\nu}_{e}$ (before
setting small flavor perturbations; see below). Focusing on fast flavor
conversions, we neglect the momentum changing collisions and only keep the
neutrino–neutrino forward scattering contributions Fuller  _et al._ (1987);
Pantaleone (1992); Sigl and Raffelt (1993) in the effective Hamiltonian, i.e.,
we drop the terms originating from vacuum neutrino mixing and neutrino–matter
forward scattering contributions.

Assuming the differential neutrino angular distribution
$dn_{\nu_{e}(\bar{\nu}_{e})}/d\Omega$ is uniform in $z$ and taking the two-
flavor approximation, the equation of motion, which governs the evolution of
the normalized neutrino density matrices (see below) $\varrho$ (for neutrinos)
and $\bar{\varrho}$ (for antineutrinos) is given as follows111When only
considering the neutrino–neutrino self-interaction term, one can redefine the
antineutrino density matrix by
$\bar{\rho}\rightarrow-\sigma_{y}^{\dagger}\bar{\rho}\sigma_{y}$ such that
neutrinos and antineutrinos are treated in equal footing Duan  _et al._
(2006a). However, here we choose to treat them separately for the purpose of
consistency with follow-up works that may include other terms in the
Hamiltonian and the collisions.:

$\displaystyle\frac{\partial}{\partial
t}\varrho(t,z,v_{z})+v_{z}\frac{\partial}{\partial
z}\varrho(t,z,v_{z})=-i[H(t,z,v_{z}),\varrho(t,z,v_{z})],$ (1a)
$\displaystyle\frac{\partial}{\partial
t}\bar{\varrho}(t,z,v_{z})+v_{z}\frac{\partial}{\partial
z}\bar{\varrho}(t,z,v_{z})=-i[\bar{H}(t,z,v_{z}),\bar{\varrho}(t,z,v_{z})].$
(1b)

with

$\varrho(t,z,v_{z})=\left[\begin{array}[]{cc}\varrho_{ee}&\varrho_{ex}\\\
\varrho_{ex}^{*}&\varrho_{xx}\\\
\end{array}\right],~{}\bar{\varrho}(t,z,v_{z})=\left[\begin{array}[]{cc}\bar{\varrho}_{ee}&\bar{\varrho}_{ex}\\\
\bar{\varrho}_{ex}^{*}&\bar{\varrho}_{xx}\\\ \end{array}\right],$ (2)

in the flavor basis. Since we omit the vacuum mixing contribution, the flavor
evolution of $\rho$ and $\bar{\rho}$ do not depend on the neutrino energy.

In Eq. (1), the Hamiltonian $H$ and $\bar{H}$ can be explicitly written down
as:

$\displaystyle H(t,z,v_{z})=$
$\displaystyle~{}\mu\int_{-1}^{1}dv_{z}^{\prime}(1-v_{z}v_{z}^{\prime})\times$
(3a) $\displaystyle[g_{\nu}(v_{z}^{\prime})\varrho(t,z,v_{z}^{\prime})-\alpha
g_{\bar{\nu}}(v_{z}^{\prime})\bar{\varrho}^{*}(t,z,v_{z}^{\prime})],$
$\displaystyle\bar{H}(t,z,v_{z})=$
$\displaystyle-\mu\int_{-1}^{1}dv_{z}^{\prime}(1-v_{z}v_{z}^{\prime})\times$
(3b)
$\displaystyle[g_{\nu}(v_{z}^{\prime})\varrho^{*}(t,z,v_{z}^{\prime})-\alpha
g_{\bar{\nu}}(v_{z}^{\prime})\bar{\varrho}(t,z,v_{z}^{\prime})],$

where $\mu=\sqrt{2}G_{F}n_{\nu_{e}}$ with $G_{F}$ being the Fermi coupling
constant. The asymmetry parameter $\alpha=n_{\bar{\nu}_{e}}/n_{\nu_{e}}$
quantifies the ratio between the number density of $\bar{\nu}_{e}$ and
$\nu_{e}$. The angular distribution function $g_{\nu(\bar{\nu})}(v_{z})$
relates to the physical neutrino phase-space distribution function
$f_{\nu_{e}(\bar{\nu}_{e})}$ without any flavor perturbations as

$g_{\nu(\bar{\nu})}(v_{z})=\frac{1}{4\pi^{2}n_{\nu_{e}(\bar{\nu}_{e})}}\int
dE_{\nu}E_{\nu}^{2}f_{\nu_{e}(\bar{\nu}_{e})},$ (4)

where $E_{\nu}$ is the energy of a neutrino. Note that the azimuthal symmetry
in phase space of $f_{\nu_{e}(\bar{\nu}_{e})}$ with respect to the $z$
direction is implicitly taken. In this notation, $\rho_{ee}=\bar{\rho}_{ee}=1$
for our system without any initial flavor perturbations while all other matrix
elements are zero.

Since $\mu$ is the only dimensional quantity in Eq. (1), we define $\mu\equiv
1$ and express $t$ and $z$ in dimensionless form (in the units of $\mu^{-1}$)
hereafter.

###  II.2 Initial and boundary conditions

Figure 1: The initial $\nu$ELN angular distribution $G(v_{z})$ for systems
with different $\bar{\nu}_{e}$ to $\nu_{e}$ number density ratios $\alpha$.
The $\nu$ELN crossing where $G(v_{z})=0$ shifts to smaller $v_{z}$ with larger
$\alpha$.

We simulate the flavor evolution of neutrinos inside a box of size
$L_{z}=1200$ in $z\in[-600,600]$ with the periodic boundary condition. Similar
to Ref. Martin  _et al._ (2020), we parametrize $g_{\nu(\bar{\nu})}$ as

$g_{\nu(\bar{\nu})}(v_{z})\propto\exp[-(v_{z}-1)^{2}/(2\sigma_{\nu(\bar{\nu})}^{2})],$
(5)

with normalization condition $\int_{-1}^{1}dv_{z}g_{\nu(\bar{\nu})}(v_{z})=1$.
Throughout the paper, we fix $\sigma_{\nu}=0.6$ and $\sigma_{\bar{\nu}}=0.5$,
such that $\bar{\nu}_{e}$ are more forward-peaked than $\nu_{e}$. For the
asymmetry parameter $\alpha$, we take $\alpha=0.9,1.0$, $1.1$, $1.2$, and
$1.3$. Figure 1 shows the corresponding $\nu$ELN distributions
$G_{\nu}(v_{z})\equiv g_{\nu}(v_{z})-\alpha g_{\bar{\nu}}(v_{z})$. With
increasing value of $\alpha$, the $\nu$ELN crossing where $G_{\nu}(v_{z,c})=0$
occurs at smaller $v_{z,c}$ with $|G_{\nu}(v_{z})|$ being larger (smaller) for
$v_{z}>v_{z,c}$ ($v_{z}<v_{z,c}$). The values of $v_{z,c}$ for
$\alpha=0.9,1.0,1.1,1.2$, and $1.3$ are 0.65, 0.45, 0.33, 0.23, and 0.15,
correspondingly.

For a system without the vacuum Hamiltonian, it requires a small perturbation
in $\rho_{ex}$ ($\bar{\rho}_{ex}$) to trigger the flavor instabilities. For
simplicity, we assume that the perturbations are independent of $v_{z}$.
Specifically, we use the following description for our initial condition at
$t=0$ as:

$\displaystyle\varrho_{ee}(z,v_{z})=\bar{\varrho}_{ee}(z,v_{z})=(1+\sqrt{1-\epsilon^{2}(z)})/2,$
(6a)
$\displaystyle\varrho_{xx}(z,v_{z})=\bar{\varrho}_{xx}(z,v_{z})=(1-\sqrt{1-\epsilon^{2}(z)})/2,$
(6b)
$\displaystyle\varrho_{ex}(z,v_{z})=\bar{\varrho}_{ex}(z,v_{z})=\epsilon(z)/2.$
(6c)

For $\epsilon(z)$, we explore the following two different choices:

  1. 1.

Point-source like perturbations centered at $z=0$:
$\epsilon(z)=\epsilon_{0}\exp[-z^{2}/50]$.

  2. 2.

Random perturbations: $\epsilon(z)$ is randomly assigned by a real value
between $0$ and $\epsilon_{0}$.

We take $\epsilon_{0}=10^{-6}$ throughout the paper.

###  II.3 Numerical methods

We discretize $z$ and $v_{z}$ with $N_{z}$ and $N_{v_{z}}$ grids, and evolve
$\varrho$ and $\bar{\varrho}$ defined at the grid centers. In evaluating the
advection terms $\partial\varrho/{\partial z}$ and
$\partial\bar{\varrho}/{\partial z}$, we use two different methods to check
for consistency. The first method uses the forth-order (central) finite
difference method with artificial dissipation using the third-order Kreiss-
Oliger formulation Kreiss and Oliger (1973). For the second method, we use
finite volume method plus the seventh order weighted essentially
nonoscillatory (WENO) scheme Shu (1998, 2003). For time evolution, we use
fourth-order Runge-Kutta method. The numerical implementation details will be
reported in a separate publication together with the public release of our
simulation code COSE$\nu$ George  _et al._ (2021).

Using the above initial and boundary conditions, we perform simulations with
fiducial numerical parameters of $N_{Z}=12000$ and $N_{v_{z}}=200$, with a
fixed time step size $\Delta t=C_{\rm CFL}\Delta z$, where $\Delta
z=L_{z}/N_{Z}$ and the Courant–Friedrichs–Lewy number $C_{\rm CFL}=0.4$. In
Appendixes A and B, we show the comparison of results obtained with different
resolutions and with two different numerical methods for advection. For the
rest of the paper, all results are obtained using the finite volume method
with seventh order WENO scheme, which provides better numerical accuracy given
the same resolution.

###  II.4 Polarization vectors

Before we discuss our results, let us define the neutrino and antineutrino
polarization vectors $\mathbf{P}$ and $\mathbf{\bar{P}}$ that are often used
in literature. The three components of the polarization vector $\mathbf{P}$
are defined as

$\displaystyle P_{1}$ $\displaystyle=2{\rm Re}(\varrho_{ex}),$ (7a)
$\displaystyle P_{2}$ $\displaystyle=-2{\rm Im}(\varrho_{ex}),$ (7b)
$\displaystyle P_{3}$ $\displaystyle=(\varrho_{ee}-\varrho_{xx}),$ (7c)

which satisfy $\rho\equiv(\mathbf{P}\cdot\bm{\sigma}+P_{0}\mathcal{I})/2$,
where $\sigma_{i}$ are the Pauli matrices and $\mathcal{I}$ is the identity
matrix. For antineutrinos, we have

$\displaystyle{\bar{P}}_{1}$ $\displaystyle=2{\rm Re}(\bar{\varrho}_{ex}),$
(8a) $\displaystyle{\bar{P}}_{2}$ $\displaystyle=2{\rm
Im}(\bar{\varrho}_{ex}),$ (8b) $\displaystyle{\bar{P}}_{3}$
$\displaystyle=(\bar{\varrho}_{ee}-\bar{\varrho}_{xx}).$ (8c)

With this definition, the equation of motion for $\mathbf{P}$ has the same
form as for $\mathbf{\bar{P}}$, and we can treat neutrinos and antineutrinos
in equal footing using the $\nu$ELN spectrum $G_{\nu}$ defined in Sec. II.2 to
account for the contributions from both neutrinos and antineutrinos in the
Hamiltonian. Dropping the bar for antineutrinos, we have:

$\left(\frac{\partial}{\partial t}+v_{z}\frac{\partial}{\partial
z}\right)\mathbf{P}(t,z,v_{z})=\mathbf{H}(t,z,v_{z})\times\mathbf{P}(t,z,v_{z}),$
(9)

where
$\mathbf{H}(t,z,v_{z})=\int_{-1}^{1}dv^{\prime}_{z}G_{\nu}(v^{\prime}_{z})\mathbf{P}(t,z,v^{\prime}_{z})(1-v_{z}v^{\prime}_{z})$.
Clearly, in the absence of collisions, the length of the polarization vectors
are unity for all $z$ and $v_{z}$, given that our initial condition has
$|\mathbf{P}|=1$ uniformly in $z$. Another physical quantity that is conserved
globally is the net neutrino $e-x$ lepton number in the entire simulation
domain Bhattacharyya and Dasgupta (2020, 2021):

$L_{e-x}\propto\int_{-L_{z}/2}^{L_{z}/2}dz\int_{-1}^{1}dv_{z}G_{\nu}(v_{z})P_{3}(t,z,v_{z}).$
(10)

Note that locally this net lepton number can be transferred from one location
to another due to the advection. However, globally $L_{e-x}$ must be a
conserved quantity in time for the entire system.

Figure 2: Panels (a) and (b): The dispersion relation of $\omega(k_{z})$ for
unstable solutions containing nonzero $\rm{Im}(\omega)$ for $\alpha=0.9$ and
$1.1$. Panels (c) and (d): Time evolution of the perpendicular components of
the neutrino polarization vectors $P_{\perp}(z)$ in the linear regime
($|P_{\perp}|\ll 1$) for $v_{z}=0.995$. Darker curves correspond to earlier
times. The same point-source like flavor perturbations are applied for both
cases with $\alpha=0.9$ [panel (c)] and $\alpha=1.1$ [panel (d)].

##  III Results: flavor evolution of the system

In this section, we focus on results obtained with $\alpha=0.9$ and
$\alpha=1.1$. In Sec. III.1, we examine how the point-source like
perturbations in the off diagonal elements of the density matrices grow in the
linear regime and compare their evolution with detailed analysis using the
dispersion relation. In Secs. III.2 and III.3, we further discuss how flavor
evolution proceeds in the nonlinear regime for cases with point-source like
perturbations and random perturbations, respectively.

###  III.1 Linear regime

In the regime where $|\rho_{ex}|\ll 1$ and $|\bar{\rho}_{ex}|\ll 1$, the
evolution of the system can be qualitatively understood by performing the
standard technique of linearized instability analysis Banerjee  _et al._
(2011); Izaguirre  _et al._ (2017). By inserting
$\rho_{ex}\propto\exp[-i(\omega t-k_{z}z)]$ (same for $\bar{\rho}_{ex}$) and
keeping terms in the lowest order of perturbations, Eq. (1) leads to a
dispersion relation of the collective mode of the system $\omega(k_{z})$
Izaguirre  _et al._ (2017); Capozzi  _et al._ (2017); Yi  _et al._ (2019). The
complex branch of $\omega$ for real $k_{z}$ signifies the existence of flavor
instabilities such that small off diagonal perturbations can grow
exponentially in time as the system evolves. Note that here we only include
the solution that preserves the the axial symmetry. For the symmetry-breaking
solutions, we omit them here because our simulation setup does not allow
symmetry-breaking solutions.

In Fig. 2, we show in panels (a) and (b) the dispersion relation
$\omega(k_{z})$ for systems with $\nu$ELN asymmetry parameter $\alpha=0.9$ and
$1.1$. For $\alpha=0.9$, there are two regions in $k_{z}$ where complex
$\omega$ exist. The maximal value of Im$(\omega_{\rm max})\simeq 0.044$
locates at $k_{z,{\rm max}}\simeq-0.165$. At the same $k_{z,{\rm max}}$,
$(d\omega/dk_{z})_{\rm max}\simeq 0.515$. For $\alpha=1.1$ where the $\nu$ELN
crossing is deeper in larger $v_{z}$ (see Fig. 1), Im$(\omega_{\rm max})\simeq
0.107$ at $k_{z,{\rm max}}\simeq-0.225$. The corresponding
$(d\omega/dk_{z})_{\rm max}\simeq 0.278$.

In the bottom panels [(c) and (d)] of Fig. 2, we show the evolution of
$|P_{\perp}(z)|=\sqrt{P_{1}^{2}+P_{2}^{2}}=2|\rho_{ex}|$ for $v_{z}=0.995$ in
the liner regime for the same $\nu$ELNs with point-source like perturbations.
Clearly, an initial Gaussian packet of $P_{\perp}$ grows exponentially with
time and the peak position of $P_{\perp}$ moves toward positive $z$ direction
for both the cases. The growth rate of the maximal value of $P_{\perp,{\rm
max}}$ and the velocity $dz(P_{\perp,{\rm max}})/dt$ both agree well with the
values of Im$(\omega_{\rm max})$ and $(d\omega/dk_{z})_{\rm max}$ obtained in
the dispersion relation above. The growth of perturbations using $\alpha=1.1$
is indeed much faster than that with $\alpha=0.9$. Moreover, the
$P_{\perp}(z)$ with $\alpha=1.1$ spreads over both positive and negative $z$
directions. For the case with $\alpha=0.9$, although the width of
$P_{\perp}(z)$ also increases with time, it mainly drifts to the positive $z$
direction without spreading over the negative $z$ direction.

Figure 3: Snapshots of $P_{3}(z,v_{z})$ at different simulation times for the
same systems shown in Fig. 2 when the initial perturbations have grown to the
nonlinear regime. Flavor waves with coherent structures develop and propagate.
Small scale structures form when flavor waves interact and breaks the coherent
structure.

###  III.2 Nonlinear regime: Point-source like perturbations

When the perturbations grow to the nonlinear regime, wavelike oscillatory
features develop Martin  _et al._ (2020). In Fig. 3, we show the snapshots of
$P_{3}(z,v_{z})$ at different times for $\alpha=0.9$ and $1.1$ using the same
point-source like perturbations as in Sec. III.1. For the case of
$\alpha=0.9$, the flavor evolution behavior of the system is in general
similar to those reported in Ref. Martin  _et al._ (2020). Flavor waves
develop and mainly propagate toward the positive $z$ direction. This is
consistent with the growth of perturbations in the linear regime discussed
earlier. An interesting feature shown here is that although the flavor
oscillations initially can affect all $v_{z}$ modes, illustrated by the
vertical stripes in the upper two subpanels in Fig. 3(a), this effect
diminishes when time proceeds and flavor conversions are roughly confined in
$v_{z}\gtrsim v_{z,c}\simeq 0.65$.

Another interesting feature is when the forward-traveling wave front interacts
with the slowly backward propagating part after $t\gtrsim 1200$, it pushes the
whole pattern to drift toward positive $z$ direction. Meanwhile, this
interaction breaks the coherent wavelike pattern. Substructures in smaller
scale develop such that the orientation of the neutrino polarization vectors
varies rapidly in $z$ for $v_{z}\gtrsim v_{z,c}$. Consequently, although at
each location $z$, neutrinos with different $v_{z}\gtrsim v_{z,c}$ still have
$|\mathbf{P}|=1$, the “average polarization vector” over the $z$ domain can
shrink due to their misalignment. Such a flavor state was referred as “flavor
depolarization” in Refs. Bhattacharyya and Dasgupta (2020, 2021) and we will
discuss its behavior in more detail in the next section. For $v_{z}\lesssim
v_{z,c}$, most neutrinos remain unaffected.

For $\alpha=1.1$ shown in panel (b), flavor conversions quickly develop toward
both the positive and negative $z$ directions and produces coherent and
wavelike structure, once again consistent with that indicated by the growth of
perturbations in the linear regime. Similarly, when the forward and backward
propagating modes interact after $t\gtrsim 665$, smaller structures develop
and cause a major part of $v_{z}$ reaching closer to flavor depolarization.
One important difference from the previous $\alpha=0.9$ case is now flavor
depolarization happens mostly in $v_{z}\lesssim v_{z,c}\simeq 0.45$. We will
also discuss this feature and its consequences in Sec. IV.2.

For the other $\nu$ELN spectra with $\alpha=1.0$, $1.2$, and $1.3$, the
behaviors are qualitatively similar to that with $\alpha=1.1$. Full simulation
animations are available at https://mengruwuu.github.io/COSEv1dbox/.

###  III.3 Nonlinear regime: Random perturbations

Figure 4: Snapshots of $P_{3}(z,v_{z})$ for systems with the same $\nu$ELN
shown in Fig. 3 but with different seed perturbations. With random
perturbations, interaction of flavor waves happen much earlier. No large-scale
coherent structures can be formed when compared to Fig. 3, which adopts point-
source like perturbations.

Now, let us look at how the flavor systems evolve when we adopt different seed
of random perturbations. We show again $P_{3}(z,v_{z})$ at different time
snapshots in Fig. 4 for $\alpha=0.9$ and $1.1$. With random perturbations,
some locations have larger $P_{\perp}$ initially such that flavor conversions
develop faster around those locations (see the top subpanels). Similar to what
discussed in Sec. III.2 for single point-source like perturbations, flavor
waves can initially develop independently now at different locations. The
flavor waves transported in space interact and break the coherent pattern to
form small-scale structures. Thus, with perturbations everywhere in space, the
interaction of flavor waves happen at much earlier times and no large-scale
coherent structure can be formed, differently from cases with point-source
like perturbations. The systems can then reach closer to flavor depolarization
in $v_{z}\gtrsim v_{z,c}$ ($v_{z}\lesssim v_{z,c}$) for $\alpha=0.9$ ($1.1$)
within a much shorter amount of time.

Once again, for the other $\nu$ELN spectra with $\alpha=1.0$, $1.2$, and
$1.3$, the behaviors are qualitatively similar to that with $\alpha=1.1$. Full
simulation animations are also available at
https://mengruwuu.github.io/COSEv1dbox/.

##  IV Results: evolution of averaged quantities

In the previous Sec. III, we have discussed how fast flavor conversions of
neutrinos can develop in space and time with different seed perturbations. In
this section, we further examine the time evolution of relevant quantities
averaged over $z$ and/or $v_{z}$ in Sec. IV.1 and Sec. IV.2.

###  IV.1 Overall flavor conversion probability in the box

Figure 5: Time evolution of the flavor survival probabilities [see Eq. (11)]
for $\nu_{e}$ (upper subpanels) and $\bar{\nu}_{e}$ (lower subpanels) averaged
over $z$ and $v_{z}$ for systems with point-source like perturbations [panel
(a)] and random perturbations [panel (b)] with different $\bar{\nu}_{e}$ to
$\nu_{e}$ number density ratios $\alpha$. See text in Sec. IV.1 for details.

We define the overall flavor survival probabilities of neutrinos and
antineutrinos, $\overline{\langle P_{ee}\rangle}$ and
$\overline{\langle{\bar{P}}_{ee}\rangle}$ in the simulation domain by

$\displaystyle\overline{\langle P_{ee}\rangle}=\int
dzdv_{z}\varrho_{ee}(z,v_{z})g_{\nu}(v_{z})\Bigg{/}\int
dzdv_{z}g_{\nu}(v_{z}),$ (11a)
$\displaystyle\overline{\langle{\bar{P}}_{ee}\rangle}=\int
dzdv_{z}\bar{\varrho}_{ee}(z,v_{z})g_{\bar{\nu}}(v_{z})\Bigg{/}\int
dzdv_{z}g_{\bar{\nu}}(v_{z}),$ (11b)

where the brackets and overline denote averaging over $z$ and $v_{z}$,
respectively. The integration limits over $z$ and $v_{z}$ are from $-L_{z}/2$
to $L_{z}/2$ and from $-1$ to $1$, respectively.

Fig. 5 shows the time evolution of $\overline{\langle P_{ee}\rangle}$ and
$\overline{\langle{\bar{P}}_{ee}\rangle}$ for all five different $\nu$ELN
spectra that we considered (see Fig. 1). The cases with point-source like
perturbations and random perturbations are shown in panels (a) and (b),
respectively. First, we see that cases with larger $\alpha$ reach the final
asymptotic states earlier in time for both point-source like and random
perturbations. Meanwhile, all cases with random perturbations reach the
asymptotic states much earlier than those with point-source like perturbations
as discussed in Sec. III.3.

Second, systems with $\alpha=1$ (equal number of neutrinos and antineutrinos)
have final asymptotic values of $\overline{\langle
P_{ee}\rangle}_{f}\simeq\overline{\langle{\bar{P}}_{ee}\rangle}_{f}\simeq
0.5$, i.e., full flavor depolarization is almost reached for the entire
system. For systems that are more asymmetric, the asymptotic
$\overline{\langle P_{ee}\rangle}_{f}$ and
$\overline{\langle{\bar{P}}_{ee}\rangle}_{f}$ are larger, i.e., on average
less $\nu_{e}$ and $\bar{\nu}_{e}$ are converted. Comparing $\overline{\langle
P_{ee}\rangle}_{f}$ to $\overline{\langle{\bar{P}}_{ee}\rangle}_{f}$ for a
given $\alpha$, one finds that cases with $\alpha>1$ have larger values of
asymptotic $\overline{\langle{\bar{P}}_{ee}\rangle}_{f}$ than those of
$\overline{\langle P_{ee}\rangle}_{f}$ (vice versa for $\alpha<1$). This is
related to the fact that for $\alpha>1$ ($\alpha<1$), neutrinos with
$v_{z}\lesssim v_{z,c}$ ($v_{z}\gtrsim v_{z,c}$) experience more flavor
conversions. Since for our adopted neutrino angular spectra,
$g_{\bar{\nu}}(v_{z})$ is more forward peaked in positive $v_{z}$ than
$g_{\nu}(v_{z})$, more flavor conversions in larger $v_{z}$ naturally lead to
smaller $\overline{\langle{\bar{P}}_{ee}\rangle}_{f}$ than $\overline{\langle
P_{ee}\rangle}_{f}$.

Figure 6: The final flavor survival probabilities $\langle
P_{ee}(v_{z})\rangle_{f}$ averaged over $z$ [Eq. (12)] for the same systems
shown in Fig. 5. Nearly flavor depolarization can be reached in one side of
the $\nu$ELN spectra for cases with $\alpha\neq 1$. For $\alpha=1$, the entire
$\nu$ELN spectra can achieve nearly flavor depolarization.

Comparing the final values of $\overline{\langle P_{ee}\rangle}_{f}$ and
$\overline{\langle{\bar{P}}_{ee}\rangle}_{f}$ for cases with point-source like
and random perturbations, the point-source like ones generally lead to a
slightly smaller $\overline{\langle P_{ee}\rangle}_{f}$ and
$\overline{\langle{\bar{P}}_{ee}\rangle}_{f}$ than the random ones. This is
due to the reason discussed in Sec. III.2: although the interaction of the
flavor waves lead to states close to flavor depolarization, for cases with
point-source like perturbations, some large-scale coherent structures can
remain until the end of the simulations. This can also be seen in Fig. 6,
which shows the spatially averaged flavor survival probabilities $\langle
P_{ee}(v_{z})\rangle_{f}$ as functions of $v_{z}$ at the end of our
simulations, where $\langle P_{ee}(v_{z})\rangle$ is defined by

$\displaystyle\langle P_{ee}(v_{z})\rangle=\int
dz\varrho_{ee}(z,v_{z})\Bigg{/}\int dz.$ (12a)

Note that $\langle{\bar{P}}_{ee}(v_{z})\rangle$ that can be similarly defined
are equal to $\langle P_{ee}(v_{z})\rangle$, when only neutrino-neutrino self-
interaction terms are included here.

Figure 6 clearly shows several interesting features. First, for cases with
either point-source like or random perturbations, those with $\alpha>1$
($\alpha<1$) lead to nearly full flavor depolarization ($\langle
P_{ee}(v_{z})\rangle_{f}\simeq 0.5$) for velocity modes $v_{z}<v_{z,c}$
($v_{z}>v_{z,c}$). For velocity modes in the other side of the $\nu$ELN,
nearly flavor equipartition cannot be achieved and $\langle
P_{ee}(v_{z})\rangle_{f}$ gradually increase from $\simeq 0.5$ to larger
values for larger (smaller) $v_{z}$ for $\alpha>1$ ($\alpha<1$). For the case
with $\alpha=1$, i.e., symmetric in neutrinos and antineutrinos, full flavor
depolarization for the entire $\nu$ELN can be achieved. Comparing results with
different perturbations, one sees that due to the remaining large-scale
coherent structures (see e.g., Fig. 3), systems with point-source like
perturbations give rise to smaller $\langle P_{ee}(v_{z})\rangle_{f}\simeq
0.4$ in some regions in $v_{z}<v_{z,c}$ ($v_{z}>v_{z,c}$) for $\alpha>1$
($\alpha<1$) than those with random perturbations.

Based on our simulation results, we make a conjecture as follows. For a system
with the periodic boundary condition in one spatial direction (while having
perfect translation symmetry in the other two directions) where interaction of
flavor waves can happen, the system can evolve to an asymptotic state where
nearly flavor equipartition can be reached for velocity modes in one side of
the $\nu$ELN, when averaged over space. In other words, when averaged over
space, the $\nu$ELN crossing nearly vanishes. For such a scenario to happen,
the net neutrino $e-x$ lepton number conservation $L_{e-x}$ [see Eq. (10)]
would enforce that nearly flavor depolarization can only happen for velocity
modes in the shallower side of the $\nu$ELN.

An additional remark is: comparing Fig. 5 and Fig. 6, we also point out that
although Fig. 6 seems to suggest that a narrower range of velocity modes in
$\nu$ELN distribution reaches nearly flavor depolarization for $\alpha=0.9$
than those with $\alpha>1$, on average, more neutrinos and antineutrinos are
being converted. Once again, this is because the spectra $g_{\nu}(v_{z})$ and
$g_{\bar{\nu}}(v_{z})$ are more forward peaked. This means that although a
system with $\alpha<1$ may appears to be affected in less phase space volume
in $v_{z}$, flavor conversions there can actually have a larger impact on
physical processes that are flavor dependent in realistic astrophysical
environments.

###  IV.2 Evolution of different velocity modes

Figure 7: Time evolution of different spatially averaged quantities. See text
in Sec. IV.2 for details.

Refs. Bhattacharyya and Dasgupta (2020, 2021) proposed that the behavior of
the system can be understood by examining the time evolution of
$\langle\mathbf{P}(v_{z})\rangle$. From Eq. (9), one obtains

$\frac{d}{dt}\langle\mathbf{P}(v_{z},t)\rangle=\langle\mathbf{M}_{0}\times\mathbf{P}(v_{z},t)\rangle-
v_{z}\langle\mathbf{M}_{1}(t)\times\mathbf{P}(v_{z},t)\rangle,$ (13)

where $\mathbf{M}_{i}(z,t)\equiv\int
dv_{z}L_{i}(v_{z})G_{\nu}(v_{z})\mathbf{P}(z,v_{z},t)$ with $L_{i}(v_{z})$ the
Legendre Polynomials. The authors of Refs. Bhattacharyya and Dasgupta (2020,
2021) argued that the spatial averaging of the cross products in Eq. (13) can
be approximated as the cross products of two vectors that are spatially
averaged separately:

$\frac{d}{dt}\langle\mathbf{P}(v_{z},t)\rangle\simeq[\langle\mathbf{M}_{0}\rangle-
v_{z}\langle\mathbf{M}_{1}(t)\rangle]\times\langle\mathbf{P}(v_{z},t)\rangle.$
(14)

By doing so, the time evolution of $\langle\mathbf{P}(v_{z},t)\rangle$ may be
understood as a vector precessing around an effective Hamiltonian vector
$\langle\mathbf{H}(v_{z})\rangle=\langle\mathbf{M}_{0}\rangle-
v_{z}\langle\mathbf{M}_{1}(t)\rangle$. Moreover, the depolarization of
$\langle\mathbf{P}(v_{z},t)\rangle$ happens when the $\langle
H\rangle_{\perp}=\sqrt{\langle H\rangle_{1}^{2}+\langle H\rangle_{2}^{2}}$
roughly exceeds $|\langle{\tilde{H}}\rangle_{3}|$ Bhattacharyya and Dasgupta
(2021), where $\langle{\tilde{H}}\rangle_{3}$ differs from $\langle
H\rangle_{3}$ by a factor related to $\langle\mathbf{M}_{2}\rangle$ when going
to a co-rotating frame (see Ref. Bhattacharyya and Dasgupta (2020) for
details).

We carefully examine the validity of these claims. In Fig. 7, we show in panel
(a) the time evolution of $|\langle\tilde{H}\rangle_{3}|$ and $\langle
H_{\perp}\rangle$ for $\alpha=0.9$ of $v_{z}=1$ and $\alpha=1.1$ of $v_{z}=-1$
using random perturbations as examples. Despite these modes undergo nearly
flavor depolarization, Fig. 7(a) shows that the values of $\langle
H_{\perp}\rangle$ remain much smaller than $|\langle\tilde{H}\rangle_{3}|$
nearly during the whole time, in contrast to what discussed in Ref.
Bhattacharyya and Dasgupta (2021).

In panels (b) and (c) of Fig. 7, we show the evolution of
$|\langle\mathbf{P}(v_{z})\rangle|$ and the relative angles $\theta_{r}$
between the two vectors on the right-hand sides of Eqs. (13) and (14) for ten
$v_{z}$ modes for the case with $\alpha=1.1$ with random perturbations.
Clearly, for modes with $v_{z}\lesssim v_{z,c}$ that reach flavor
depolarization (darker colors), i.e.,
$|\langle\mathbf{P}(v_{z})\rangle|\rightarrow 0$, their $\cos\theta_{r}$
oscillate rapidly between -1 and 1. This indicates that Eq. (14) is in fact
not a good approximation of Eq. (13).

In sum, our findings here suggest that the picture which explains the late-
time behavior and the depolarization of the system proposed in Refs.
Bhattacharyya and Dasgupta (2020, 2021) using the spatially averaged
polarization vectors needs to be revisited.

##  V Summary and discussions

We have performed long-term numerical simulations of fast collective neutrino
flavor conversions for systems with translation symmetry in two spatial ($x$
and $y$) directions using a newly developed code COSE$\nu$. Our numerical
calculations were conducted based on the finite volume method with seventh
order WENO scheme and the finite difference method with Kreiss-Oliger
dissipation. Both methods showed good capability of numerical error
suppression which allows simulations being carried to longer than a few
thousand characteristic timescale of the system. Assuming uniform $\nu$ELN
distributions in $z$ and taking periodic boundary conditions, we have
investigated how flavor conversions happen with different $\nu$ELN spectra,
controlled by the initial $\bar{\nu}_{e}$ to $\nu_{e}$ asymmetric ratio
parameter $\alpha$, as well as how the results depend on the chosen flavor
perturbations. Our main findings can be summarized as follows.

First, we found that for systems with point-source like initial perturbations,
flavor waves with coherent structures can develop and propagate. With periodic
boundary conditions, these flavor waves can then interact and break into
smaller scale structures. When adopting perturbation seed randomly placed in
$z$, the flavor waves originated from different locations can interact much
faster such that no coherent structure can be formed.

The interactions of the flavor waves lead the system to a final state where
part of the velocity space ($v_{z}$) is close to flavor depolarization when
averaging over the space. For any asymmetric system with $\alpha\neq 1$, only
one side of the $\nu$ELN spectrum relative to the $\nu$ELN crossing point can
reach close to an averaged flavor depolarization while neutrinos in the other
side of $\nu$ELN experience less flavor conversions, as constrained by the
conservation of the net neutrino $e-x$ number. Specifically, for $\alpha>1$
($\alpha<1$), nearly flavor depolarization can be reached for $v_{z}\lesssim
v_{z,c}$ ($v_{z}\gtrsim v_{z,c}$) when the antineutrinos angular distribution
are more forward peaked than that of neutrinos. On the other hand, for systems
with $\alpha=1$, the entire neutrino spectra can reach close to averaged
flavor depolarization. This phenomenon is qualitatively similar for systems
with either point-source like or random perturbations. Quantitatively, the
developed large-scale coherent pattern in cases with point-source like
perturbations allow part of the velocity space to have on average more flavor
conversions than perfect flavor depolarization.

Comparing our results with those reported in Refs. Martin  _et al._ (2020);
Bhattacharyya and Dasgupta (2020, 2021), our results with point-source like
perturbations agree with Martin  _et al._ (2020), which, however, only evolves
the system for a shorter period of time without allowing flavor waves to
interact. On the other hand, Refs. Bhattacharyya and Dasgupta (2020, 2021)
obtained results nearly independent of whether the initial perturbations are
being point-source like or random. In fact, behavior of random perturbations
emerged in their simulations using point-source like perturbations, different
from what we obtained here. One potential reason is that Refs. Bhattacharyya
and Dasgupta (2020, 2021) used fast Fourier transform to evaluate the
derivative terms, which might artificially generate errors of random nature in
the spatial domain. Moreover, we have tried to verify the mechanisms proposed
in Refs. Bhattacharyya and Dasgupta (2020, 2021) in explaining the occurrence
of the flavor depolarization. However, our results do not support the proposed
mechanisms and suggest that better understanding is needed.

Practically, our numerical findings may provide insights for efforts which
attempt to include the impact of fast neutrino flavor conversions in
hydrodynamical simulations; e.g., Li and Siegel (2021). For instance, if the
adopted neutrino transport scheme can provide the antineutrino-to-neutrino
asymmetry ratio and the crossing points in the $\nu$ELN spectrum, partial
flavor depolarization to one end of the $\nu$ELN together with the net $e-x$
neutrino lepton number may be applied to the neutrino distribution functions
at where $\nu$ELN crossings are identified.

Needless to say, there are still several improvements to be made in future.
For example, our simulations were performed in reduced dimensions. How the
symmetry-breaking solutions may develop and affect our conclusions need to be
further examined. In our simulations, we have completely omitted the vacuum
and the matter terms in the Hamiltonian, as well as the collision terms.
Adding these terms and incorporating full treatment with three neutrino
flavors may introduce new effects recently investigated in works without
advection Capozzi  _et al._ (2020); Shalgar and Tamborra (2021a, b); Martin
_et al._ (2021). Full inclusion of them are to be implemented in future. Last
but not least, the potential impact of the many-body nature of the problem
remain to be further elucidated Rrapaj (2020); Cervia  _et al._ (2019);
Roggero (2021a, b).

###### Acknowledgements.

We thank Geng-Yu Liu and Herlik Wibowo for useful discussions as well as
Basudeb Dasgupta, Soumya Bhattacharyya, and Irene Tamborra for valuable
feedback. M.-R. W. and M. G. acknowledge supports from the Ministry of Science
and Technology, Taiwan under Grants No. 109-2112-M-001-004, No. 110-2112-M-001
-050, and the Academia Sinica under Project No. AS-CDA-109-M11. M.-R. W. also
acknowledge supports from the Physics Division, National Center for
Theoretical Sciences, Taiwan. M.-R. W. and M. G. appreciate the computing
resources provided by the Academia Sinica Grid-computing Center. C.-Y. L.
thank the National Center for High-performance Computing (NCHC) for providing
computational and storage resources. Z. X. acknowledge supports from the
European Research Council (ERC) under the European Union’s Horizon 2020
research and innovation programme (ERC Advanced Grant KILONOVA No. 885281).

##  Appendix A Convergence tests

Figure 8: Panel (a): Difference in $\overline{\langle P_{ee}\rangle}$ between
runs using different $(N_{z},N_{v_{z}})$ and $\overline{\langle
P_{ee}\rangle}_{0}$ obtained with $(N_{z},N_{v_{z}})=(12000,200)$ with
$\alpha=0.9$ and point-source like perturbations. Panels (b) and (c):
Evolution of $\overline{\langle\delta|\mathbf{P}|\rangle}$ and
$(\delta|\mathbf{P}|)_{\rm max}$ with different simulation resolutions.

In this appendix, we quantify errors associated with the resolution of our
numerical simulations. As described in the main text, our fiducial resolution
is with $N_{z}=12000$ and $N_{v_{z}}=200$. In panel (a) of Fig. 8, we compare
the absolute difference of $\overline{\langle P_{ee}\rangle}$ between runs
with different resolutions and our fiducial case, taking $\alpha=0.9$ with
point-source like perturbations. It shows that the difference of
$\overline{\langle P_{ee}\rangle}$ amounts to $\sim 10^{-3}$ and $\sim
10^{-2}$, when we reduce $N_{v_{z}}$ by a factor of $2$ and $4$, respectively
while keeping $N_{z}$ being the same. For runs with fixed $N_{v_{z}}=200$ but
with reduced $N_{z}=6000$ and $3000$, the relative differences are smaller
than $10^{-6}$ and are thus not shown in the figure. This can also be inferred
by looking at the nearly identical black and red curves, which represent the
differences with $(N_{z},N_{v_{z}})=(12000,100)$ and $(6000,100)$,
respectively.

However, this does not mean that we can adopt a much reduced resolution in
$N_{z}$. In panels (b) and (c) of Fig. 8, we show the deviation of the length
of the polarization vector from unity, averaged spatially and spectrally over
$g_{\nu}(v_{z})$:

$\overline{\langle\delta|\mathbf{P}|}\rangle=\int
dzdv_{z}g_{\nu}(v_{z})(|\mathbf{P}(z,v_{z})|-1)\Bigg{/}\int
dzdv_{z}g_{\nu}(v_{z}),$ (15)

and the maximal deviation of the length of the polarization vectors in the
entire simulation domain $(\delta|\mathbf{P}|)_{\rm max}$. Both panels clearly
show that the errors are mostly depending on the resolution in $N_{z}$. For
our fiducial case with $(N_{z},N_{v_{z}})=(12000,200)$, the associated errors
in $\overline{\langle\delta|\mathbf{P}|\rangle}$ and
$(\delta|\mathbf{P}|)_{\rm max}$ are smaller than $10^{-8}$ and $10^{-5}$,
respectively, for $\alpha=0.9$ with point-source like perturbations. For all
of our simulations with different $\alpha$ and different perturbation seeds,
we obtain $\overline{\langle\delta|\mathbf{P}|}\rangle<10^{-4}$. For
$(\delta|\mathbf{P}|)_{\rm max}$, all cases have $(\delta|\mathbf{P}|)_{\rm
max}<10^{-2}$ except for $\alpha=1.2$ and $1.3$ with point-source like
perturbations, for which their $(\delta|\mathbf{P}|)_{\rm max}$ reach
$\mathcal{O}(1)$ by the end of the simulation. We note, however, that
demanding $(\delta|\mathbf{P}|)_{\rm max}\ll 1$ may not be a practical
convergence criterion because $(\delta|\mathbf{P}|)_{\rm max}$ is usually
associated with the grids with largest $|v_{z}|$, which has minor contribution
to the Hamiltonian and the averaged quantities. In fact, we have compared the
errors in the averaged quantities for all cases with lowered resolutions
($N_{z}=6000$) and confirmed that all the averaged quantities agree within
$5\times 10^{-3}$.

##  Appendix B Results using different numerical schemes

Figure 9: Comparison of the same quantities as those in Fig. 8 for simulations
with $\alpha=0.9$ and point-source like perturbations, using the finite volume
method with seventh order WENO scheme (FV+WENO) and the finite difference
method supplied by the third-order Kreiss-Oliger dissipation term (FD+KO3).
Note that in panels (b) and (c), we show additionally the quantities derived
using the finite difference method only without dissipation (FD only).

In this appendix, we compare our results obtained with two different numerical
methods described in Sec. II.3. Panel (a) of Fig. 9 shows the difference of
$\overline{\langle P_{ee}\rangle}$ between the run using the finite volume
method with seventh order WENO scheme (FV+WENO) and that with the finite
difference method supplied by the third-order Kreiss-Oliger dissipation term
(FD+KO3), for $\alpha=0.9$ with point-source like perturbations and our
fiducial resolution of $(N_{z},N_{v_{z}})=(12000,200)$. One sees that the
differences are smaller than $10^{-5}$ throughout the whole time. Panels (b)
and (c) of the same figure show once again
$\overline{\langle\delta|\mathbf{P}|\rangle}$ and $(\delta|\mathbf{P}|)_{\rm
max}$. Here we see that although both quantities remain much smaller than one
with either the FV+WENO or the FD+KO3 scheme, the FV+WENO scheme gives rise to
errors much smaller than those using FD+KO3 scheme. In addition, we show
values of $\overline{\langle\delta|\mathbf{P}|\rangle}$ and
$(\delta|\mathbf{P}|)_{\rm max}$ using simply the finite difference method
without applying any error suppression mechanisms (FD only) in panels (b) and
(c) of Fig. 9. It shows clearly that both our FV+WENO and FD+KO3 schemes help
suppress numerical errors in $\overline{\langle\delta|\mathbf{P}|\rangle}$ by
more than two orders of magnitudes for this particular parameter set.

## References

  * Zyla  _et al._ (2020) P. Zyla _et al._ (Particle Data Group), PTEP 2020, 083C01 (2020).
  * Duan  _et al._ (2006a) H. Duan, G. M. Fuller, and Y.-Z. Qian, Phys. Rev. D 74, 123004 (2006a), arXiv:astro-ph/0511275 .
  * Duan  _et al._ (2006b) H. Duan, G. M. Fuller, J. Carlson, and Y.-Z. Qian, Phys. Rev. D74, 105014 (2006b), arXiv:astro-ph/0606616 [astro-ph] .
  * Hannestad  _et al._ (2006) S. Hannestad, G. G. Raffelt, G. Sigl, and Y. Y. Y. Wong, Phys. Rev. D 74, 105010 (2006), [Erratum: Phys.Rev.D 76, 029901 (2007)], arXiv:astro-ph/0608695 .
  * Esteban-Pretel  _et al._ (2008) A. Esteban-Pretel, A. Mirizzi, S. Pastor, R. Tomas, G. G. Raffelt, P. D. Serpico, and G. Sigl, Phys. Rev. D 78, 085012 (2008), arXiv:0807.0659 [astro-ph] .
  * Dasgupta  _et al._ (2009) B. Dasgupta, A. Dighe, G. G. Raffelt, and A. Y. Smirnov, Phys. Rev. Lett. 103, 051105 (2009), arXiv:0904.3542 [hep-ph] .
  * Malkus  _et al._ (2012) A. Malkus, J. P. Kneller, G. C. McLaughlin, and R. Surman, Phys. Rev. D 86, 085015 (2012), arXiv:1207.6648 [hep-ph] .
  * Cherry  _et al._ (2012) J. F. Cherry, J. Carlson, A. Friedland, G. M. Fuller, and A. Vlasenko, Phys. Rev. Lett. 108, 261104 (2012), arXiv:1203.1607 [hep-ph] .
  * Raffelt  _et al._ (2013) G. Raffelt, S. Sarikas, and D. de Sousa Seixas, Phys. Rev. Lett. 111, 091101 (2013), [Erratum: Phys. Rev. Lett.113,no.23,239903(2014)], arXiv:1305.7140 [hep-ph] .
  * Vlasenko  _et al._ (2014) A. Vlasenko, G. M. Fuller, and V. Cirigliano, Phys. Rev. D89, 105004 (2014).
  * Volpe  _et al._ (2013) C. Volpe, D. Väänänen, and C. Espinoza, Phys.Rev. D87, 113010 (2013).
  * Wu  _et al._ (2016) M.-R. Wu, H. Duan, and Y.-Z. Qian, Phys. Lett. B752, 89 (2016), arXiv:1509.08975 [hep-ph] .
  * Abbar and Duan (2015) S. Abbar and H. Duan, Phys. Lett. B 751, 43 (2015), arXiv:1509.01538 [astro-ph.HE] .
  * Sawyer (2016) R. F. Sawyer, Phys. Rev. Lett. 116, 081101 (2016), arXiv:1509.03323 [astro-ph.HE] .
  * Izaguirre  _et al._ (2017) I. Izaguirre, G. Raffelt, and I. Tamborra, Phys. Rev. Lett. 118, 021101 (2017), arXiv:1610.01612 [hep-ph] .
  * Capozzi  _et al._ (2019) F. Capozzi, B. Dasgupta, A. Mirizzi, M. Sen, and G. Sigl, Phys. Rev. Lett. 122, 091101 (2019), arXiv:1808.06618 [hep-ph] .
  * Abbar (2021) S. Abbar, Phys. Rev. D 103, 045014 (2021), arXiv:2007.13655 [astro-ph.HE] .
  * Xiong and Qian (2021) Z. Xiong and Y.-Z. Qian, Phys. Lett. B 820, 136550 (2021), arXiv:2104.05618 [astro-ph.HE] .
  * Johns (2021) L. Johns, (2021), arXiv:2104.11369 [hep-ph] .
  * Duan  _et al._ (2010) H. Duan, G. M. Fuller, and Y.-Z. Qian, Ann. Rev. Nucl. Part. Sci. 60, 569 (2010).
  * Mirizzi  _et al._ (2015) A. Mirizzi, G. Mangano, and N. Saviano, Phys. Rev. D92, 021702 (2015), arXiv:1503.03485 [hep-ph] .
  * Duan (2015) H. Duan, Int. J. Mod. Phys. E24, 1541008 (2015), arXiv:1506.08629 [hep-ph] .
  * Tamborra and Shalgar (2020) I. Tamborra and S. Shalgar,  (2020), 10.1146/annurev-nucl-102920-050505, arXiv:2011.01948 [astro-ph.HE] .
  * Duan  _et al._ (2011) H. Duan, A. Friedland, G. C. McLaughlin, and R. Surman, J. Phys. G38, 035201 (2011).
  * Wu  _et al._ (2015) M.-R. Wu, Y.-Z. Qian, G. Martinez-Pinedo, T. Fischer, and L. Huther, Phys. Rev. D91, 065016 (2015).
  * Sasaki  _et al._ (2017) H. Sasaki, T. Kajino, T. Takiwaki, T. Hayakawa, A. B. Balantekin, and Y. Pehlivan, Phys. Rev. D 96, 043013 (2017), arXiv:1707.09111 [astro-ph.HE] .
  * Wu  _et al._ (2017) M.-R. Wu, I. Tamborra, O. Just, and H.-T. Janka, Phys. Rev. D96, 123015 (2017), arXiv:1711.00477 [astro-ph.HE] .
  * Wu and Tamborra (2017) M.-R. Wu and I. Tamborra, Phys. Rev. D95, 103007 (2017), arXiv:1701.06580 [astro-ph.HE] .
  * Stapleford  _et al._ (2020) C. J. Stapleford, C. Fröhlich, and J. P. Kneller, Phys. Rev. D 102, 081301 (2020), arXiv:1910.04172 [astro-ph.HE] .
  * Xiong  _et al._ (2020) Z. Xiong, A. Sieverding, M. Sen, and Y.-Z. Qian, Astrophys. J. 900, 144 (2020), arXiv:2006.11414 [astro-ph.HE] .
  * George  _et al._ (2020) M. George, M.-R. Wu, I. Tamborra, R. Ardevol-Pulpillo, and H.-T. Janka, Phys. Rev. D 102, 103015 (2020), arXiv:2009.04046 [astro-ph.HE] .
  * Li and Siegel (2021) X. Li and D. M. Siegel, Phys. Rev. Lett. 126, 251101 (2021), arXiv:2103.02616 [astro-ph.HE] .
  * Dasgupta  _et al._ (2017) B. Dasgupta, A. Mirizzi, and M. Sen, JCAP 1702, 019 (2017), arXiv:1609.00528 [hep-ph] .
  * Capozzi  _et al._ (2017) F. Capozzi, B. Dasgupta, E. Lisi, A. Marrone, and A. Mirizzi, Phys. Rev. D 96, 043016 (2017), arXiv:1706.03360 [hep-ph] .
  * Abbar and Volpe (2019) S. Abbar and M. C. Volpe, Phys. Lett. B 790, 545 (2019), arXiv:1811.04215 [astro-ph.HE] .
  * Yi  _et al._ (2019) C. Yi, L. Ma, J. D. Martin, and H. Duan, Phys. Rev. D 99, 063005 (2019), arXiv:1901.01546 [hep-ph] .
  * Martin  _et al._ (2020) J. D. Martin, C. Yi, and H. Duan, Phys. Lett. B 800, 135088 (2020), arXiv:1909.05225 [hep-ph] .
  * Padilla-Gay  _et al._ (2021) I. Padilla-Gay, S. Shalgar, and I. Tamborra, JCAP 01, 017 (2021), arXiv:2009.01843 [astro-ph.HE] .
  * Johns  _et al._ (2020a) L. Johns, H. Nagakura, G. M. Fuller, and A. Burrows, Phys. Rev. D 102, 103017 (2020a), arXiv:2009.09024 [hep-ph] .
  * Bhattacharyya and Dasgupta (2020) S. Bhattacharyya and B. Dasgupta, Phys. Rev. D 102, 063018 (2020), arXiv:2005.00459 [hep-ph] .
  * Bhattacharyya and Dasgupta (2021) S. Bhattacharyya and B. Dasgupta, Phys. Rev. Lett. 126, 061302 (2021), arXiv:2009.03337 [hep-ph] .
  * Morinaga (2021) T. Morinaga, (2021), arXiv:2103.15267 [hep-ph] .
  * Richers  _et al._ (2021) S. Richers, D. E. Willcox, N. M. Ford, and A. Myers, Phys. Rev. D 103, 083013 (2021), arXiv:2101.02745 [astro-ph.HE] .
  * Zaizen and Morinaga (2021) M. Zaizen and T. Morinaga, (2021), arXiv:2104.10532 [hep-ph] .
  * Kato  _et al._ (2021) C. Kato, H. Nagakura, and T. Morinaga, (2021), arXiv:2108.06356 [astro-ph.HE] .
  * Morinaga  _et al._ (2020) T. Morinaga, H. Nagakura, C. Kato, and S. Yamada, Phys. Rev. Res. 2, 012046 (2020), arXiv:1909.13131 [astro-ph.HE] .
  * Nagakura  _et al._ (2019) H. Nagakura, T. Morinaga, C. Kato, and S. Yamada,  (2019), 10.3847/1538-4357/ab4cf2, arXiv:1910.04288 [astro-ph.HE] .
  * Johns  _et al._ (2020b) L. Johns, H. Nagakura, G. M. Fuller, and A. Burrows, Phys. Rev. D 101, 043009 (2020b), arXiv:1910.05682 [hep-ph] .
  * Delfan Azari  _et al._ (2020) M. Delfan Azari, S. Yamada, T. Morinaga, H. Nagakura, S. Furusawa, A. Harada, H. Okawa, W. Iwakami, and K. Sumiyoshi, Phys. Rev. D 101, 023018 (2020), arXiv:1910.06176 [astro-ph.HE] .
  * Abbar  _et al._ (2020) S. Abbar, H. Duan, K. Sumiyoshi, T. Takiwaki, and M. C. Volpe, Phys. Rev. D 101, 043016 (2020), arXiv:1911.01983 [astro-ph.HE] .
  * Glas  _et al._ (2020) R. Glas, H. T. Janka, F. Capozzi, M. Sen, B. Dasgupta, A. Mirizzi, and G. Sigl, Phys. Rev. D 101, 063001 (2020), arXiv:1912.00274 [astro-ph.HE] .
  * Abbar (2020) S. Abbar, JCAP 05, 027 (2020), arXiv:2003.00969 [astro-ph.HE] .
  * Johns and Nagakura (2021) L. Johns and H. Nagakura, Phys. Rev. D 103, 123012 (2021), arXiv:2104.04106 [hep-ph] .
  * Nagakura  _et al._ (2021) H. Nagakura, L. Johns, A. Burrows, and G. M. Fuller, (2021), arXiv:2108.07281 [astro-ph.HE] .
  * Fuller  _et al._ (1987) G. M. Fuller, R. W. Mayle, J. R. Wilson, and D. N. Schramm, Astrophys. J. 322, 795 (1987).
  * Pantaleone (1992) J. T. Pantaleone, Phys.Lett. B287, 128 (1992).
  * Sigl and Raffelt (1993) G. Sigl and G. Raffelt, Nucl.Phys. B406, 423 (1993).
  * Kreiss and Oliger (1973) H. O. Kreiss and J. Oliger, GARP publication series No. 10, Geneva (1973).
  * Shu (1998) C.-W. Shu, “Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws,” in  _Advanced Numerical Approximation of Nonlinear Hyperbolic Equations: Lectures given at the 2nd Session of the Centro Internazionale Matematico Estivo (C.I.M.E.) held in Cetraro, Italy, June 23–28, 1997_, edited by A. Quarteroni (Springer Berlin Heidelberg, Berlin, Heidelberg, 1998) pp. 325–432.
  * Shu (2003) C.-W. Shu, International Journal of Computational Fluid Dynamics 17, 107 (2003), https://doi.org/10.1080/1061856031000104851 .
  * George  _et al._ (2021) M. George, C.-Y. Lin, G.-Y. Liu, M.-R. Wu, and Z. Xiong, in preparation (2021).
  * Banerjee  _et al._ (2011) A. Banerjee, A. Dighe, and G. Raffelt, Phys. Rev. D 84, 053013 (2011), arXiv:1107.2308 [hep-ph] .
  * Capozzi  _et al._ (2020) F. Capozzi, M. Chakraborty, S. Chakraborty, and M. Sen, Phys. Rev. Lett. 125, 251801 (2020), arXiv:2005.14204 [hep-ph] .
  * Shalgar and Tamborra (2021a) S. Shalgar and I. Tamborra, JCAP 01, 014 (2021a), arXiv:2007.07926 [astro-ph.HE] .
  * Shalgar and Tamborra (2021b) S. Shalgar and I. Tamborra, Phys. Rev. D 103, 063002 (2021b), arXiv:2011.00004 [astro-ph.HE] .
  * Martin  _et al._ (2021) J. D. Martin, J. Carlson, V. Cirigliano, and H. Duan, Phys. Rev. D 103, 063001 (2021), arXiv:2101.01278 [hep-ph] .
  * Rrapaj (2020) E. Rrapaj, Phys. Rev. C 101, 065805 (2020), arXiv:1905.13335 [hep-ph] .
  * Cervia  _et al._ (2019) M. J. Cervia, A. V. Patwardhan, A. B. Balantekin, t. S. N. Coppersmith, and C. W. Johnson, Phys. Rev. D 100, 083001 (2019), arXiv:1908.03511 [hep-ph] .
  * Roggero (2021a) A. Roggero, (2021a), arXiv:2102.10188 [hep-ph] .
  * Roggero (2021b) A. Roggero, (2021b), arXiv:2103.11497 [hep-ph] .


# Iterative Model Predictive Control for Piecewise Systems

Ugo Rosolia and Aaron D. Ames Ugo Rosolia and Aaron D. Ames are with the AMBER
lab at Caltech, Pasadena, USA. E-mails: {urosolia<EMAIL_ADDRESS>The
authors would like to acknowledge the support of the NSF award #1932091 .

###### Abstract

In this paper, we present an iterative Model Predictive Control (MPC) design
for piecewise nonlinear systems. We consider finite time control tasks where
the goal of the controller is to steer the system from a starting
configuration to a goal state while minimizing a cost function. First, we
present an algorithm that leverages a feasible trajectory that completes the
task to construct a control policy which guarantees that state and input
constraints are recursively satisfied and that the closed-loop system reaches
the goal state in finite time. Utilizing this construction, we present a
policy iteration scheme that iteratively generates safe trajectories which
have non-decreasing performance. Finally, we test the proposed strategy on a
discretized Spring Loaded Inverted Pendulum (SLIP) model with massless legs.
We show that our methodology is robust to changes in initial conditions and
disturbances acting on the system. Furthermore, we demonstrate the
effectiveness of our policy iteration algorithm in a minimum time control
task.

##  I Introduction

Robots performing complex tasks can be described as hybrid systems, which are
characterized by continuous dynamics and discrete events. Therefore,
controllers designed for such systems can take control actions based on
continuous and discrete decision variables. Yet the presence of discrete
variables make planning and control problems challenging, as it is required to
reason about all possible combinations of discrete events. This challenge can
be mitigated by designing hierarchical strategies, where a high-level planner
computes the discrete variables and a low-level controller optimizes the
system trajectory described by continuous variables [1, 2, 3, 4].

A popular methodology to synthesize policies, which can jointly plan over
discrete and continuous states is Model Predictive Control (MPC) [5, 6, 7, 8].
MPC is a control strategy which systematically uses forecast to compute
control actions. At each time step, an MPC plans a trajectory over a short
time window, then the first control action is applied to the system and the
process is repeated at the next time step based on new measurements. When the
system dynamics are hybrid, the MPC planning problem is a Mixed Integer
Program (MIP) that is hard to solve online with limited computational
resources. For this reason, significant work has focused on explicit MPC
strategies where the solution to the MIP is solved offline as a parametric
optimization problem [9, 7, 8]. For hybrid systems described by piecewise
affine dynamics the parametric optimization problem can be solved exactly.
Once the solution is computed offline, the MPC policy is given by a look-up
table of feedback gains that can be efficiently implemented online in real-
time [10, 11]. However, computing the explicit solution to hybrid MPC problems
is computationally demanding.

Another strategy to speed-up the computation of the MPC policy is to leverage
warm-starting strategies, where the optimization algorithm is initialized
using a candidate solution. Several strategies have been proposed for warm-
starting hybrid MPC problems [12, 13, 14, 15]. These approaches leverage the
trajectory computed at the previous time step to warm-start both the
continuous and discrete variables. As the complexity of solving MIPs is given
by the computation of the optimal integer variables, recent works have
investigated the possibility of leveraging learning algorithms to predict the
set of active discrete variables used to warm-start the MPC [16, 17, 18].

In this work, we focus on control tasks where the goal is to steer the system
from a starting configuration to a goal state in finite time, while satisfying
state and input constraints. We assume that a feasible trajectory that is able
to perform the task is available. Then, we synthesize a control policy, which
plans the system trajectory over a finite horizon that is shorter than the
control task duration and may cause the controller to take unsafe shortsighted
control actions. Thus, building upon [19, 20], we present a methodology to
construct the MPC terminal components in order to guarantee satisfaction of
the safety constraints and convergence in finite time of the closed-loop
system to the goal set.

Compared to previous works [19, 20], we show how to handle piecewise systems
by warm-starting the integer variables and we present a shrinking horizon
strategy tailored to finite time control tasks. We present an algorithm which
solves at most $M$ smooth optimization problems and is guaranteed to find a
feasible solution to the original MIP planning problem. Our approach is based
on a sub-optimal trajectory that can complete the task and affects the closed-
loop performance of the controller. Therefore, we present a policy iteration
algorithm, where simulations are used to iteratively update the controller. We
prove that our algorithm returns safe policies that have non-decreasing
performance. Finally, we demonstrate the effectiveness of our approach on the
discretized Spring Loaded Inverted Pendulum (SLIP) [21].

##  II Problem Formulation

In this section, we describe the system model and the control synthesis
objectives. We consider discrete time piecewise nonlinear systems defined over
$R$ disjoint regions $\mathcal{D}_{i}\subseteq\mathbb{R}^{n}$ for
$i\in\\{1,\ldots,R\\}$:

$x_{t+1}=f_{i}(x_{t},u_{t}),\quad\mbox{if }x_{t}\in\mathcal{D}_{i},$ (1)

where the state $x_{t}\in\mathbb{R}^{n}$ and the input
$u_{t}\in\mathbb{R}^{d}$. In the above equation
$f_{i}:\mathbb{R}^{n}\times\mathbb{R}^{d}\rightarrow\mathbb{R}^{n}$ represents
the system dynamics, which describe the evolution of the discrete time system
when the state $x_{t}$ belongs to the region
$\mathcal{D}_{i}\subseteq\mathbb{R}^{n}$. Furthermore, the system is subject
to the following state and input constraints:

$u_{t}\in\mathcal{U}\subseteq\mathbb{R}^{d}\text{ and
}x_{t}\in\mathcal{X}\subseteq\mathbb{R}^{n},~{}\forall
t\in\\{0,\ldots,T-1\\},$ (2)

where $T\in\\{0,1,\ldots\\}$ is the duration of the control task. Notice that
several robotic systems can be described by piecewise constrained nonlinear
models defined over disjoint regions, such as the SLIP model presented in
Section V-A.

Objective:  Our goal is to design a control policy
$\pi:\mathbb{R}^{n}\rightarrow\mathbb{R}^{d}$ which maps states to actions,
i.e.,

$u_{t}=\pi(x_{t}).$ (3)

The above control policy should guarantee that the state and input constraints
from (2) are satisfied and that the closed-loop system (1) and (3) converges
in finite time to a goal set $\mathcal{G}\subset\mathcal{X}$. More formally,
the control policy (3) should guarantee that for an initial condition $x_{I}$
in a neighborhood of a starting state $x_{S}\in\mathcal{X}$, the trajectory of
the closed-loop system (1) and (3) is a feasible solution to the following
Finite Time Optimal Control Problem (FTOCP):

$\displaystyle\min_{\begin{subarray}{c}u_{0},\ldots,u_{T-1}\\\
i_{0},\ldots,i_{T-1}\end{subarray}}$
$\displaystyle\sum_{t=0}^{T-1}l(x_{t},u_{t})$ (4) s.t. $\displaystyle
x_{t+1}=f_{i_{t}}(x_{t},u_{t}),$ $\displaystyle
x_{t}\in\mathcal{D}_{i_{t}}\cap\mathcal{X},~{}u_{t}\in\mathcal{U},~{}i_{t}\in\\{1,\ldots,R\\},$
$\displaystyle x_{0}=x_{I},x_{T}\in\mathcal{G},$ $\displaystyle\forall
t\in\\{0,\ldots,T-1\\},$

where the stage cost
$l:\mathbb{R}^{n}\times\mathbb{R}^{d}\rightarrow\mathbb{R}$. Notice that in
the above problem the system dynamics are a function of the integer variables
$i_{t}\in\\{1,\ldots,R\\}$. Therefore, for a feasible set of continuous inputs
$[u_{0},\ldots,u_{T-1}]$ and integer variables $[i_{0},\ldots,i_{T-1}]$, we
have that the resulting vector of states $[x_{0},\ldots,x_{T}]$ must satisfy
$x_{t}\in\mathcal{D}_{i_{t}}\cap\mathcal{X},~{}\forall
t\in\\{0,\ldots,T-1\\}$. Throughout the paper we make the following
assumptions.

###### Assumption 1.

We are given the state-input trajectories which are feasible for the FTOCP (4)
with $x_{I}=x_{S}$:

${\mathbf{x}}^{0}=[x_{0}^{0},\ldots,x_{T}^{0}]\text{ and
}{\mathbf{u}}^{0}=[u_{0}^{0},\ldots,u_{T-1}^{0}],$

where for all $t\in\\{0,\ldots,T-1\\}$ the state $x_{t}^{0}\in\mathcal{X}$,
the input $u_{t}^{0}\in\mathcal{U}$ and $x_{T}^{0}\in\mathcal{G}$.

###### Assumption 2.

For any $x\in\mathcal{G}$, the stage cost $l(x,u)=0$, for all
$u\in\mathcal{U}$. Moreover, the set $\mathcal{G}$ is control invariant, i.e,
for all $x\in\mathcal{G}$, there exists a control $u\in\mathcal{U}$ and index
$i\in\\{1,\ldots,R\\}$ such that $f_{i}(x,u)\in\mathcal{G}$ and
$x\in\mathcal{D}_{i}$.

###### Remark 1.

The proposed methodology requires only feasibility of a trajectory
${\mathbf{x}}^{0}$. However, the optimality of the trajectory
${\mathbf{x}}^{0}$ affects the performance of the proposed control synthesis
strategy. For this reason, in Section III-C we present a policy iteration
scheme that may be used to iteratively improve the closed-loop performance of
the policy.

##  III Control Design

In this section, we first introduce an FTOCP which can be recast as a non-
linear program (NLP). Afterwards, we present the proposed strategy which
leverages this NLP and the feasible trajectory from Assumption 1. Finally, we
present a policy iteration scheme which can be used to improve the performance
of the control policy.

###  III-A Model Predictive Control

The FTOCP (4) is challenging to solve as at each time $t$ the system dynamics
change as a function of the state $x_{t}$, i.e., $x_{t+1}=f_{i}(x_{t},u_{t})$
if $x_{t}\in\mathcal{D}_{i}$. However, the computational complexity may be
reduced when searching for a feasible sub-optimal solution. In particular, a
feasible solution to the FTOCP (4) may be computed fixing a priori a sequence
of regions $\\{\mathcal{D}_{i_{0}},\ldots,\mathcal{D}_{i_{T-1}}\\}$, where the
system should be constrained at each time $t$. Clearly, a trajectory which
steers the system from the starting state $x_{I}$ to the goal set
$\mathcal{G}$, while visiting the sequence of regions
$\\{\mathcal{D}_{i_{0}},\ldots,\mathcal{D}_{i_{T-1}}\\}$ is a feasible
trajectory for the original FTOCP (4).

In order to reduce the computational complexity, we introduce an FTOCP defined
over a horizon $N$ shorter than $T$ and for a set of indices
$\mathcal{I}_{t}=\\{i_{t},\ldots,i_{t+N-1}\\}$ associated with a sequence of
$N$ regions $\\{\mathcal{D}_{i_{t}},\ldots,\mathcal{D}_{i_{t+N-1}}\\}$. In
particular, given a set of indices $\mathcal{I}_{t}$, the terminal state
$x_{F}$ and an associated terminal cost $q_{F}\in\mathbb{R}$ we define the
FTOCP:

$\displaystyle J(x_{t},x_{F},q_{F},\mathcal{I}_{t},N)=\min_{\mathbf{u}_{t}}$
$\displaystyle\sum_{k=t}^{t+N-1}l(x_{k|t},u_{k|t})+q_{F}$ (5) s.t.
$\displaystyle x_{k+1|t}=f_{i_{k}}(x_{k|t},u_{k|t}),$ $\displaystyle
x_{k|t}\in\mathcal{D}_{i_{k}}\cap\mathcal{X},~{}u_{k|t}\in\mathcal{U},$
$\displaystyle x_{t|t}=x_{t},x_{t+N|t}=x_{F},$ $\displaystyle\forall
k\in\\{t,\ldots,t+N-1\\},$

where $\mathbf{u}_{t}=[u_{t|t},\ldots,u_{t+N-1|t}]$ The optimal state-input
sequences to the above FTOCP

$[x_{t|t}^{*},\ldots,x_{t+N|t}^{*}]\text{ and
}[u_{t|t}^{*},\ldots,u_{t+N-1|t}^{*}],$ (6)

steer the system from the starting state $x_{t}$ to the terminal state $x_{F}$
while satisfying state and input constraints.

In the FTOCP (5), at each predicted time $k$ the system state
$x_{k|t}\in\mathcal{D}_{i_{k}}\cap\mathcal{X}$, and therefore problem (5) can
be recast as an NLP, which is easier to solve than problem (4) where the
optimization is carried out over continuous and integer variables. Next, we
present the proposed algorithm which chooses the set
$\mathcal{I}_{t}=\\{i_{t},\ldots,i_{t+N-1}\\}$ associated with the sequence of
regions $\\{\mathcal{D}_{i_{t}},\ldots,\mathcal{D}_{i_{t+N-1}}\\}$, the
terminal state $x_{F}$, and terminal cost $q_{F}$ that are used in the FTOCP
(5).

Algorithm 1 Control Policy $\pi$

1:Init Parameters: ${\mathbf{q}}^{0}$, ${\mathbf{x}}^{0}$, ${\mathbf{i}}^{0}$,
$M$, $k_{0}=N$, $N_{0}=N$, $T$

2:Input: $x_{t}$

3:$i_{t}=\texttt{getRegion}(x_{t})$

4:for $m=[0,\ldots,M-1]$ do $\triangleright$ Solve $M$ FTOCPs

5: set $t_{F}=\min(k_{t}+m,T)$

6: set $x_{F}=x^{0}_{t_{F}}$

7: set $q_{F}=q^{0}_{t_{F}}$

8: set
$\mathcal{I}_{t}=\\{i_{t},i^{0}_{t_{F}-N_{t}+1},\ldots,i^{0}_{t_{F}-1}\\}$

9: solve the FTOCP $J(x_{t},x_{F},q_{F},\mathcal{I}_{t},N_{t})$ from (5)

10: store $c_{m}=J(x_{t},x_{F},q_{F},\mathcal{I}_{t},N_{t})$

11: store $\bar{u}_{m}=u^{*}_{t|t}$

12: store $\bar{x}_{F,m}=x^{*}_{t+N_{t}|t}$

13: if $m>0$ and $c_{m-1}<c_{m}$ then

14: $c^{*}_{t}=c_{m-1}$ $\triangleright$ Pick best cost

15: $m^{*}_{t}=m-1$ $\triangleright$ Pick best cost index

16: break

17: end if

18:end for

19:if $\bar{x}_{F,m^{*}_{t}}=x_{T}^{0}$ then $\triangleright$ Horizon and
Parameter Update

20: store $N_{t+1}=\max(1,N_{t}-1)$

21: store $k_{t+1}=T$

22:else

23: store $N_{t+1}=N_{t}$

24: store $k_{t+1}=k_{t}+m^{*}_{t}+1$

25:end if

26:$u_{t}=\bar{u}_{m^{*}_{t}}$

27:Outputs $c_{t}^{*}$, $u_{t}$

###  III-B Policy Synthesis

This section describes the proposed strategy. For each state $x_{t}^{0}$ of
the feasible trajectory $\mathbf{x}^{0}$ from Assumption 1, we define the
vector

${\mathbf{i}}^{0}=[i_{0}^{0},\ldots,i_{T}^{0}],$

where $i_{t}^{0}\in\\{1,\ldots,R\\}$ identifies the region containing the
system’s state at time $t$, i.e., $x_{t}^{0}\in\mathcal{D}_{i_{t}^{0}}$.
Moreover, for each state $x_{t}^{0}$ we introduce the cost-to-go $q_{t}^{0}$
given by the recursion:

$q_{t}^{0}=l(x_{k}^{0},u_{k}^{0})+q_{t+1}^{0},$ (7)

for $q_{T}^{0}=0$ and the stage cost
$l:\mathbb{R}^{n}\times\mathbb{R}^{d}\rightarrow\mathbb{R}$ from (4). Finally,
we define the cost vector

${\mathbf{q}}^{0}=[q_{0}^{0},\ldots,q_{T}^{0}].$

The cost vector ${\mathbf{q}}^{0}$, the feasible trajectory ${\mathbf{x}}^{0}$
and the vector of indices ${\mathbf{i}}^{0}$, together with the parameters
$M\in\\{0,1,2,\ldots\\}$ and $N\in\\{0,1,\ldots\\}$, are used to initialize
the proposed control policy, which is described by Algorithm 1. At each time
$t$, Algorithm 1 takes as input the state of the system $x_{t}$ and it returns
the best cost and the control action $u_{t}$, which is applied to system (1).
First, given the state $x_{t}$ we identify the region’s index $i_{t}$ such
that $x_{t}\in\mathcal{D}_{i_{t}}$ (line 1). Afterwards, we solve $M$ times
the FTOCP (5) with a different terminal state $x_{F}=x^{0}_{t_{F}}$, terminal
cost $q_{F}=q^{0}_{t_{F}}$ and set of indices $\mathcal{I}_{t}$ (lines 4-18).
Note that the set of indices $\mathcal{I}_{t}$ (line 8) is computed appending
to the current region’s index $i_{t}$ a sequence of $N_{t}-1$ indices selected
backward from time $t_{F}$, which is chosen independently of $N_{t}$ for
$t>0$. As a result, at each time $t$ the controller solves the FTOCP (4) using
a sequence of indices which is different from the one associated with the
first feasible trajectory111For example, at time $t=0$, for $m=2$, $k_{0}=N=4$
and $T=100$, we have that $t_{F}=k_{0}+m=6$; therefore
$\mathcal{I}_{0}=\\{i_{0},i^{0}_{3},i^{0}_{4},i^{0}_{5}\\}$ and
$x_{F}=x_{6}^{0}$.. As we will discuss in Section IV at time $t$ and for
$m=0$, the FTOCP (5) is feasible when the terminal components are defined as
in lines 5-7. Therefore if for $m>0$ we have that $c_{m-1}<c_{m}$, we stop
solving the set of $M$ NLPs and we save the best cost and the associated index
$m_{t}^{*}$. Then, we update the parameter $k_{t}$ that is used to define the
terminal MPC components and we shrink the prediction horizon, if the terminal
predicted state $\bar{x}_{F,m^{*}_{t}}$ equals the terminal state $x_{T}^{0}$.
Finally, we return the optimal control action $u_{t}=\bar{u}_{m^{*}_{t}}$.

###  III-C Policy Iteration

The control policy from Algorithm 1 leverages the feasible trajectory
$\mathbf{x}^{0}$ to compute the terminal components used in the MPC problem
(5) and, as a result, the performance of the proposed methodology is affected
by the optimality of $\mathbf{x}^{0}$. In this section, we discuss a policy
iteration strategy that may be used to improve the performance of the control
policy from Algorithm 1. In particular, we simulate the closed-loop system and
we iteratively update the control policy.

At iteration $j\geq 1$, we define the control policy
$\pi^{j}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{d}$ which is given by Algorithm
1 initialized using the feasible trajectory ${\mathbf{x}}^{j-1}$, the vector
of indices ${\mathbf{i}}^{j-1}$, and the cost vector ${\mathbf{q}}^{j-1}$. For
an initial condition $x_{0}^{j}=x_{S}$, the policy $\pi^{j}$ may be used to
simulate the trajectory of the closed-loop system:

$x^{j}_{t+1}=f_{i}(x_{t}^{j},\pi^{j}(x_{t}^{j})),\quad\mbox{if
}x_{t}^{j}\in\mathcal{D}_{i},$ (8)

which is then used to update the policy $\pi^{j+1}$ at the next iteration
$j+1$, as shown in Algorithm 2. In what follows, we show that the closed-loop
performance of the control policy $\pi^{j}$ is non-decreasing at each $j$th
policy update.

Algorithm 2 Iterative Policy Update

1:Init Parameters: ${\mathbf{q}}^{0}$, ${\mathbf{x}}^{0}$, ${\mathbf{i}}^{0}$,
$M$, $N$, $x_{S}$, $T$

2:Input: $j$

3:define $\pi^{1}$ via Algorithm 1 initialized with ${\mathbf{q}}^{0}$,
${\mathbf{x}}^{0}$, ${\mathbf{i}}^{0}$, $M$, $k_{0}=N$, $N_{0}=N$, $T$

4:for $i\in\\{1,\ldots,j\\}$ do $\triangleright$ Policy iteration loop

5: simulate the closed-loop system (8) from $x_{0}^{i}=x_{S}$

6: set $\mathbf{x}^{i}=[x_{0}^{i},\ldots,x_{T}^{i}]$

7: compute the region indices $\mathbf{i}^{i}$ from $\mathbf{x}^{i}$

8: compute $\mathbf{q}^{i}$ from (7) with $q_{T}^{i}=0$

9: define $\pi^{i+1}$ via Algorithm 1 initialized with ${\mathbf{q}}^{i}$,
${\mathbf{x}}^{i}$, ${\mathbf{i}}^{i}$, $M$, $k_{0}=N$, $N_{0}=N$

10:end for

11:Outputs
$(\pi^{0},{\mathbf{x}}^{0},{\mathbf{q}}^{0}),\cdots,(\pi^{j+1},{\mathbf{x}}^{j+1},{\mathbf{q}}^{j+1})$

##  IV Properties

###  IV-A Recursive Feasibility and Finite-Time Convergence

We show that at each time $t$ Algorithm 1 in closed-loop with system (1)
guarantees constraint satisfaction and that the closed-loop system converges
in finite time to the set $\mathcal{G}$.

###### Theorem 1.

Consider the closed-loop system (1) and (3), where the policy $\pi$ is given
by Algorithm 1. Let Assumptions 1-2 hold. If at time $t=0$ the initial
condition $x_{0}=x_{S}$, then the closed-loop system (1) and (3) satisfies
state and input constraints from (2) and it reaches the terminal set
$\mathcal{G}$ at time $T$, i.e.,
$x_{t}\in\mathcal{X},u_{t}\in\mathcal{U},\forall t\in\\{0,\ldots,T-1\\}\text{
and }x_{T}\in\mathcal{G}.$

Proof: The proof follows by standard MPC arguments [8, rawlings2009model] and
the construction of the time-varying component from Algorithm 1. Assume that
at time $t$ Algorithm 1 returns a feasible control action $u_{t}$ and let

$[x_{t|t}^{*},\ldots,x_{t+N_{t}|t}^{*}=x_{F}]\text{ and
}[u_{t|t}^{*},\ldots,u_{t+N_{t}-1|t}^{*}],$ (9)

be the optimal solution associated with the $m^{*}_{t}$th FTOCP (5) solved at
time $t$. Next, we consider three cases to show that the time-varying
components defining the FTOCPs solved at line 9 of Algorithm 1 guarantee that
Algorithm 1 returns a feasible control action $u_{t+1}$ at the next time step
$t+1$:

Case 1: If $\bar{x}_{F,m^{*}_{t}}=x_{T}^{0}$ and the horizon $N_{t}=1$, then
$k_{t+1}=T$, $N_{t+1}=1$ and $x_{t+1}=x_{T}^{0}$, therefore by the invariance
of $\mathcal{G}$ for $m=0$ the FTOCP
$J(x_{t+1},x_{F},q_{F},\mathcal{I}_{t+1},N_{t+1})$ with
$x_{F}=x_{k_{t+1}}=x_{T}^{0}$ is feasible.

Case 2: If $\bar{x}_{F,m^{*}_{t}}=x_{T}^{0}$ and the horizon $N_{t}>1$, then
$k_{t+1}=T$ and for $m=0$ the following state-input sequences
$[x_{t+1|t}^{*},\ldots,x_{t+N_{t}|t}^{*}]\text{ and
}[u_{t+1|t}^{*},\ldots,u_{t+N_{t}-1|t}^{*}]$ are feasible for the FTOCP
$J(x_{t+1},x_{F},q_{F},\mathcal{I}_{t+1},N_{t+1})$ with
$x_{F}=x_{k_{t+1}}=x_{T}^{0}$ as $x_{t+N_{t}|t}^{*}=x_{T}^{0}$ and
$N_{t+1}=N_{t}-1$.

Case 3: If $\bar{x}_{F,m^{*}_{t}}\neq x_{T}^{0}$, then
$k_{t+1}=k_{t}+m^{*}_{t}+1$ and for $m=0$ the state-input sequences
$[x_{t+1|t}^{*},\ldots,x_{t+N_{t}|t}^{*}=x_{k_{t}+m_{t}^{*}}^{0},x_{k_{t}+m_{t}^{*}+1}^{0}]$
and $[u_{t+1|t}^{*},\ldots,u_{t+N_{t}-1|t}^{*},u_{k_{t}+m_{t}^{*}}^{0}]$ are
feasible for the FTOCP $J(x_{t+1},x_{F},q_{F},\mathcal{I}_{t+1},N_{k+1})$ with
$x_{F}=x_{k_{t+1}}=x_{k_{t}+m_{t}^{*}+1}$.

From Cases $1$–$3$, we have that if at time $t$ Algorithm 1 returns a feasible
action $u_{t}$, then at time $t+1$ Algorithm 1 returns a feasible control
action $u_{t+1}$. Now, we notice that at time $t=0$ the sequence of actions
$[u_{0}^{0},\ldots,u_{N-1}^{0}]$ is feasible for the FTOCP
$J(x_{0},x_{F},q_{F},\mathcal{I}_{0},N_{0})$ with $x_{F}=x_{N_{0}}^{0}$, which
in turns implies that Algorithm 1 returns a feasible action $u_{t}$ at all
times and that state and input constraints are satisfied.

Finally, we show finite time convergence of the closed-loop system to the goal
set $\mathcal{G}$. From Algorithm 1, we have that $k_{t}$ increases at each
time step until $k_{t}=T$ after at most $T-N$ time steps and, afterwards, that
the horizon shrinks. Therefore, at time $T-1$ we have that $N_{T-1}=1$,
$t_{F}=k_{T-1}=T$ and $x_{F}=x_{T}^{0}$, thus the predicted state at time
$T-1$ of the optimal trajectory from (9) satisfies $x^{*}_{T|T-1}=x_{T}^{0}$,
which in turns implies that $x_{T}=x_{T|T-1}^{*}=x_{T}^{0}\in\mathcal{G}$.
$\mathbin{{\rule{5.38193pt}{5.38193pt}}}$

###### Corollary 1.

Consider the closed-loop system (1) and (3), where the policy $\pi$ is given
by Algorithm 1. Let Assumptions 1-2 hold. If at time $t$ Algorithm 1 returns a
feasible control action $u_{t}\in\mathcal{U}$, then the closed-loop system (1)
and (3) satisfies constraints (2) and it converges to the goal set
$\mathcal{G}$.

The above corollary highlights the advantage of computing a policy that maps
states to actions and it can be used to deal with perturbed initial conditions
and uncertainties. In the result section we perform an empirical study where
we test the robustness of the proposed methodology by changing initial
conditions and simulating disturbances acting on the system.

###  IV-B Iterative Improvement

This section discusses the properties of the iterative Algorithm 2. In
particular, we show that at each policy update the cumulative cost associated
with the closed-loop trajectories from (8) is non-increasing. Notice that our
strategy guarantees non-increasing cost at each update, but no guarantees are
given about the optimality of the trajectory at convergence.

###### Theorem 2.

For $i\in\\{0,\ldots,j\\}$ consider the closed-loop trajectories
$\mathbf{x}^{i}$ from Algorithm 2. If Assumptions 1-2 hold, then we have that
at each policy update the cost associated with the closed-loop trajectories is
non-increasing, i.e., $q_{0}^{i-1}\geq q_{0}^{i},~{}\forall
i\in\\{1,\ldots,j\\},$ where
$q_{0}^{i}=\sum_{t=0}^{T-1}l(x_{t}^{i},u_{t}^{i})$.

Proof: From Theorem 1 it follows that at each $j$th update the policy
$\pi^{j}$ from Algorithm 1 returns a feasible action $u_{t}^{j}$ and the
feasible state-input trajectories $\mathbf{x}^{j}$ and $\mathbf{u}^{j}$. At
time $t$, let $c^{*,j}_{t}$ be the optimal cost of the $m_{t}^{*}$th FTOCP and
let $[x_{t|t}^{*,j},\ldots,x_{t+N_{t}|t}^{*,j}]\text{ and
}[u_{t|t}^{*,j},\ldots,u_{t+N_{t}-1|t}^{*,j}]$ be the optimal solution. Then
we write the optimal cost as

$\displaystyle c_{t}^{*,j}$
$\displaystyle=\sum_{k=t}^{t+N_{t}-1}l(x_{k|t}^{*,j},u_{k|t}^{*,j})+q_{k_{t}+m_{t}^{*}}^{j-1}$
(10)
$\displaystyle=l(x_{t|t}^{*,j},u_{t|t}^{*,j})+\sum_{k=t+1}^{t+N_{t}-1}l(x_{k|t}^{*,j},u_{k|t}^{*,j})+q_{k_{t}+m_{t}^{*}}^{j-1}.$

Next, we consider three cases to analyze the time evolution of the optimal
cost $c_{t}^{*,j}$ associated with the $m_{t}^{*}$th FTOCP solved at line 9 of
Algorithm 1:

Case 1: If $\bar{x}_{F,m^{*}_{t}}=x_{T}^{j-1}$ and the horizon $N_{t}=1$, then
$k_{t+1}=T,N_{t+1}=1$, $x_{t+1}=x_{T}^{j-1}$ and the FTOCP for $m=0$ is
feasible (Case 1 of Theorem 1) which together with Assumption 2 imply that
$c_{t+1}^{*,j}=0$. Thus, we have that

$\displaystyle c_{t}^{*,j}$
$\displaystyle=l(x_{t|t}^{*,j},u_{t|t}^{*,j})+q_{m_{t}^{*}}^{j-1}=l(x_{t|t}^{*,j},u_{t|t}^{*,j})+c_{t+1}^{*,j},$
(11)

as the terminal cost $q_{m_{t}^{*}}^{j-1}=0=c_{t+1}^{*,j}$ for
$\bar{x}_{F,m^{*}_{t}}=x_{T}^{j-1}$.

Case 2: If $\bar{x}_{F,m^{*}_{t}}=x_{T}^{j-1}$ and the horizon $N_{t}>1$, then
we have that the last two terms in (10) represent the cost $\bar{c}_{t+1}^{j}$
associated with the state-input sequences
$[x_{t+1|t}^{*,j},\ldots,x_{t+N_{t}|t}^{*,j}]$ and
$[u_{t+1|t}^{*,j},\ldots,u_{t+N_{t}-1|t}^{*,j}]$, which are feasible at time
$t+1$ (Case 2 of Theorem 1) and therefore we have that

$c_{t}^{*,j}=l(x_{t|t}^{*,j},u_{t|t}^{*,j})+\bar{c}_{t+1}^{j}\geq
l(x_{t|t}^{*,j},u_{t|t}^{*,j})+c_{t+1}^{*,j}.$ (12)

Case 3: If $\bar{x}_{F,m^{*}_{t}}\neq x_{T}^{j-1}$, then by definition (7) we
have that the optimal cost can be written as

$\displaystyle c_{t}^{*,j}$
$\displaystyle=l(x_{t|t}^{*,j},u_{t|t}^{*,j})+\sum_{k=t+1}^{t+N_{t}-1}l(x_{k|t}^{*,j},u_{k|t}^{*,j})+q_{k_{t}+m_{t}^{*}}^{j-1}$
$\displaystyle=l(x_{t|t}^{*,j},u_{t|t}^{*,j})+\sum_{k=t+1}^{t+N_{t}-1}l(x_{k|t}^{*,j},u_{k|t}^{*,j})$
$\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad+l(x_{k_{t}+m_{t}^{*}}^{j-1},u_{k_{t}+m_{t}^{*}}^{j-1})+q_{k_{t}+m_{t}^{*}+1}^{j-1}.$

Notice that the last three terms in the above equation represent the open-loop
cost $\bar{c}_{t+1}^{j}$ associated with the state-input sequences
$[x_{t+1|t}^{*,j},\ldots,x_{k_{t}+m_{t}^{*}|t}^{j-1},x_{k_{t}+m_{t}^{*}+1|t}^{j-1}]$
and
$[u_{t+1|t}^{*,j},\ldots,u_{t+N_{t}-1|t}^{*,j},u_{k_{t}+m_{t}^{*}|t}^{j-1}]$,
which are feasible at time $t+1$ (Case 3 of Theorem 1) and therefore we have
that

$c_{t}^{*,j}=l(x_{t|t}^{*,j},u_{t|t}^{*,j})+\bar{c}_{t+1}^{j}\geq
l(x_{t|t}^{*,j},u_{t|t}^{*,j})+c_{t+1}^{*,j}.$ (13)

Finally, we notice that $c_{T}^{j,*}=0$ as from Theorem 1
$x_{T}^{j}\in\mathcal{G}$ and therefore from equations (11)–(13) we have that

$c_{0}^{*,j}\geq
l(x_{1}^{j},u_{1}^{j})+c_{2}^{*}\geq\ldots\geq\sum_{t=0}^{T-1}l(x_{t}^{j},u_{t}^{j})+c_{T}^{*,j}=q_{0}^{j},$

as $x_{t|t}^{*,j}=x_{t}^{j}$ and $u_{t|t}^{*,j}=u_{t}^{j}$ for all
$t\in\\{0,\ldots,T\\}$. Furthermore, at time $t=0$ and for $m=0$ the sequence
of open-loop actions $[u_{0}^{j},\ldots,u_{N-1}^{j}]$ is feasible for
$J(x_{0}^{j},x_{F},q_{F},\mathcal{I}_{0},N_{0})$ with $x_{F}=x_{N}^{j}$ and
the associated cost is $q_{0}^{j-1}$, which in turns implies $q_{0}^{j-1}\geq
c_{0}^{j}\geq q_{0}^{j}$. $\mathbin{{\rule{5.38193pt}{5.38193pt}}}$

Figure 1: Trajectory of the CoM and foot step locations associated with the
feasible trajectory $\mathbf{x}^{0}$ used to initialize the control policy
from Algorithm 1.

##  V Application to SLIP Walking

###  V-A Spring Loaded Inverted Pendulum Model

This section describes the Spring Loaded Inverted Pendulum (SLIP) model with
massless legs [21]. The system state is
$x=[p_{x},p_{y},v_{x},v_{y},z_{x}^{l},z_{x}^{r},z_{y}^{l},z_{y}^{r}]$, where
$(p_{x},p_{y})$ is the position of the Center of Mass (CoM), $(v_{x},v_{y})$
the velocity of the CoM, and $(z_{x}^{k},z_{y}^{k})$ the position of the $k$th
foot for $k\in\\{l,r\\}$. The input vector
$u=[\delta,v_{z},v_{y}^{l},v_{y}^{r}]$, where $\delta$ represents the change
of leg stiffness, $v_{z}$ the velocity along the $x$-axis of the foot not in
contact with the ground, and $(v_{y}^{l},v_{y}^{r})$ the velocities of the two
feet on the $y$-axis. Given the state of the system, we can compute the angle
$\theta^{k}=\text{tg}^{-1}((p_{x}-z_{x}^{k})/p_{y}),\forall k\in\\{l,r\\}$ and
the length $l^{k}=l_{0}-\sqrt{(p_{x}-z_{x}^{k})^{2}+p_{y}^{2}},\forall
k\in\\{l,r\\}$, where $l_{0}=0.55$ is the rest leg length, for more details
please refer to [21]. These quantities are used to define the dynamics:

$\displaystyle f(x,u)$ $\displaystyle=x$
$\displaystyle+\begin{bmatrix}v_{x}\\\ v_{y}\\\
(\delta+\delta_{0})[\gamma^{l}l^{l}\sin(\theta^{l})+\gamma^{r}l^{r}\cos(\theta^{r})]\\\
(\delta+\delta_{0})[\gamma^{l}l^{l}\cos(\theta^{l})+\gamma^{r}l^{r}\cos(\theta^{r})]-mg\\\
(1-\gamma^{l})v_{z}\\\ (1-\gamma^{r})v_{z}\\\ v_{y}^{l}\\\
v_{y}^{r}\end{bmatrix}dt$

where the integer variable $\gamma^{k}$ equals one if the $k$th foot is in
contact with the ground. In particular, we have that
$(\gamma^{l},\gamma^{r})=(1,1)\text{ if }x\in\mathcal{D}_{ds}$,
$(\gamma^{l},\gamma^{r})=(1,0)\text{ if }x\in\mathcal{D}_{l}$, and
$(\gamma^{l},\gamma^{r})=(0,1)\text{ if }x\in\mathcal{D}_{r}$, where the
regions are defined as follows:

$\displaystyle\mathcal{D}_{ds}$ $\displaystyle=\\{x\in\mathbb{R}^{n}|l^{r}\leq
l_{\mathrm{max}},l^{l}\leq l_{\mathrm{max}},z_{y}^{l}=0\text{ and
}z_{y}^{r}=0\\},$ $\displaystyle\mathcal{D}_{l}$
$\displaystyle=\\{x\in\mathbb{R}^{n}|l^{r}\leq
l_{\mathrm{max}},z_{y}^{l}=0\text{ and }z_{y}^{r}>0\\},$
$\displaystyle\mathcal{D}_{r}$ $\displaystyle=\\{x\in\mathbb{R}^{n}|l^{l}\leq
l_{\mathrm{max}},z_{y}^{l}>0\text{ and }z_{y}^{r}=0\\}.$

Finally, we notice that when a sequence of regions is fixed, the vertical
motion of the feet can be computed independently from the other states.

###  V-B Simulations

First, we compute a feasible trajectory that steers the system from standing
still at $(p_{x},p_{y})=(0,0.85)$ to a goal state
$x^{g}=(10,0.85,0,0,9.9,10.2,0,0)$. In order to compute a feasible trajectory
$\mathbf{x}^{0}$, we fixed a sequence of regions
$\\{\mathcal{D}_{t}\\}_{t=0}^{T}$ where the system should be at each time $t$
and we solved the resulting NLP with IPOPT [22] using CASADI [23].
Furthermore, we added a slack variable to the terminal constraint, we used

$\displaystyle l(x,u)$
$\displaystyle=||p_{x}-p_{x}^{\textrm{g}}||_{2}^{2}+10||p_{y}-p_{y}^{\textrm{g}}||_{2}^{2}+||v_{x}||_{2}^{2}+||v_{y}||_{2}^{2}$
(14) $\displaystyle+||\delta||_{2}^{2}+0.1||v_{z}||_{2}^{2},$

and the input constraints $\delta\in\\{\bar{\delta}:||\bar{\delta}||\leq
10\\}$ and $v_{z}\in\\{\bar{v}_{z}:||\bar{v}_{z}||\leq 10\\}$. Code available
at https://github.com/urosolia/SLIP.

Figure 2: Closed-loop trajectories of the CoM and foot steps for different
initial conditions.

Figure 1 shows the feasible trajectory which steers the system from the
starting configuration to the goal position while transitioning through phases
of double support (black) and single support with the left (red) and right
(blue) legs. The figure also shows the locations where the feet are in contact
with the ground. This feasible trajectory is used to initialize the control
policy from Algorithm 1. We tested the proposed strategy with $N=30$ on $10$
different initial conditions in the neighborhood of the starting state
$x_{S}$. We set $M=1$ to test the robustness to disturbances and changes in
initial conditions of Algorithm 1, when only one FTOCP is solve at each time
$t$. Figure 2 shows that for all initial conditions the controller is able to
steer the system to the goal state. Notice that in order to stabilize the
system the proposed strategy is able to plan a sequence of foot steps, which
are different from the one associated with the feasible trajectory. Finally,
it takes on average less than $0.05s$ to run Algorithm 1 with a maximum
computation time of $0.114s$ across all simulations.

Furthermore, we compared the proposed strategy with a tracking controller
which is defined removing the terminal constraint from (4) and using a
tracking cost instead of (14). Both our method and the tracking controller are
able complete the task. However, as shown in Figure 3, the tracking controller
fails to reach the goal when a disturbance hits the system. In Figure 3, it is
interesting to notice that the proposed approach initially deviates more from
the offline trajectory compared to the tracking controller. This deviation
allows the controller to compensate for the disturbance and to stabilize the
system back to a periodic gait.

Finally, we tested Algorithm 2 to iteratively update the control policy. We
set $j=40$ and we changed the stage cost to encode the objective of steering
the system from the starting state to the goal state in minimum time. In
particular, we defined the stage cost
$\bar{l}(x,u)=\mathds{1}_{\mathcal{G}}(x)+0.0001l(x,u)$, where the function
$\mathds{1}_{\mathcal{G}}(x)=0$ if $x\in x_{T}^{0}$ and
$\mathds{1}_{\mathcal{G}}(x)=1$ otherwise. Moreover, we set $M=40$ to allow
Algorithm 1 to solve multiple instances of the FTOCP (4). Note that as $M$
gets larger the controller can better explore the state space, but as a trade
off the computational cost increases. In this example it takes on average less
than $0.4s$ and at most $0.9s$ to run Algorithm 1. We initialized the proposed
policy iteration strategy with the feasible trajectory $\mathbf{x}^{0}$ from
Figure 1 that steers the system from the starting point to the goal state in
$396$ time steps, and our algorithm returned a policy which completes the task
in $193$ time steps. The closed-loop trajectory is shown in Figure 4. First
the controller accelerates and as a result the CoM oscillates more compared to
the first feasible trajectory from Figure 1. Finally, the controller slows
down and reduces the oscillation to reach the goal position with zero speed
and two feet on the ground.

##  VI Conclusions

We presented an algorithm to synthesize control policies for hybrid systems by
leveraging a feasible trajectory for the control task. We showed that the
proposed methodology guarantees constraints satisfaction and convergence in
finite time. Building upon the proposed synthesis strategy, we presented a
policy iteration algorithm which guarantees that at each policy update the
closed-loop performance is non-decreasing. Finally, we tested the proposed
strategy on a discretize SLIP model.

##  VII Acknowledgements

The authors would like to thank Andrew J. Taylor, Noel Csomay-Shanklin, and
Wenlong Ma for suggestions and proofreading the manuscript.

Figure 3: Comparison between an MPC tracking controller and the proposed
strategy, when a disturbance hits the system.

## References

  * [1] F. Jenelten, T. Miki, A. E. Vijayan, M. Bjelonic, and M. Hutter, “Perceptive locomotion in rough terrain–online foothold optimization,” _IEEE Robotics and Automation Letters_ , vol. 5, no. 4, pp. 5370–5376, 2020.
  * [2] O. Villarreal, V. Barasuol, P. M. Wensing, D. G. Caldwell, and C. Semini, “Mpc-based controller with terrain insight for dynamic legged locomotion,” in _2020 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2020, pp. 2436–2442.
  * [3] R. Grandia, A. J. Taylor, A. D. Ames, and M. Hutter, “Multi-layered safety for legged robots via control barrier functions and model predictive control,” _arXiv preprint arXiv:2011.00032_ , 2020.
  * [4] R. Grandia, A. J. Taylor, A. Singletary, M. Hutter, and A. D. Ames, “Nonlinear model predictive control of robotic systems with control lyapunov functions,” _arXiv preprint arXiv:2006.01229_ , 2020.
  * [5] A. Bemporad and M. Morari, “Control of systems integrating logic, dynamics, and constraints,” _Automatica_ , vol. 35, no. 3, pp. 407–427, 1999.
  * [6] A. Bemporad, F. Borrelli, and M. Morari, “Piecewise linear optimal controllers for hybrid systems,” in _Proceedings of the 2000 American Control Conference. ACC (IEEE Cat. No. 00CH36334)_ , vol. 2. IEEE, 2000, pp. 1190–1194.
  * [7] R. Oberdieck and E. N. Pistikopoulos, “Explicit hybrid model-predictive control: The exact solution,” _Automatica_ , vol. 58, pp. 152–159, 2015\.
  * [8] F. Borrelli, A. Bemporad, and M. Morari, _Predictive control for linear and hybrid systems_. Cambridge University Press, 2017.
  * [9] V. Dua, N. A. Bozinis, and E. N. Pistikopoulos, “A multiparametric programming approach for mixed-integer quadratic engineering problems,” _Computers & Chemical Engineering_, vol. 26, no. 4-5, pp. 715–733, 2002.
  * [10] I. Nașcu, R. Oberdieck, and E. N. Pistikopoulos, “Explicit hybrid model predictive control strategies for intravenous anaesthesia,” _Computers & Chemical Engineering_, vol. 106, pp. 814–825, 2017.
  * [11] G. Darivianakis, K. Alexis, M. Burri, and R. Siegwart, “Hybrid predictive control for aerial robotic physical interaction towards inspection operations,” in _2014 IEEE international conference on robotics and automation (ICRA)_. IEEE, 2014, pp. 53–58.
  * [12] T. Marcucci, R. Deits, M. Gabiccini, A. Bicchi, and R. Tedrake, “Approximate hybrid model predictive control for multi-contact push recovery in complex environments,” in _2017 IEEE-RAS 17th International Conference on Humanoid Robotics (Humanoids)_. IEEE, 2017, pp. 31–38.
  * [13] T. Marcucci and R. Tedrake, “Warm start of mixed-integer programs for model predictive control of hybrid systems,” _IEEE Transactions on Automatic Control_ , 2020.
  * [14] D. Frick, A. Domahidi, and M. Morari, “Embedded optimization for mixed logical dynamical systems,” _Computers & Chemical Engineering_, vol. 72, pp. 21–33, 2015.
  * [15] P. Hespanhol, R. Quirynen, and S. Di Cairano, “A structure exploiting branch-and-bound algorithm for mixed-integer model predictive control,” in _2019 18th European Control Conference (ECC)_. IEEE, 2019, pp. 2763–2768.
  * [16] J.-J. Zhu and G. Martius, “Fast non-parametric learning to accelerate mixed-integer programming for online hybrid model predictive control,” _arXiv preprint arXiv:1911.09214_ , 2019.
  * [17] D. Bertsimas and B. Stellato, “The voice of optimization,” _Machine Learning_ , pp. 1–29, 2020.
  * [18] A. Agrawal, S. Barratt, S. Boyd, and B. Stellato, “Learning convex optimization control policies,” in _Learning for Dynamics and Control_. PMLR, 2020, pp. 361–373.
  * [19] U. Rosolia and F. Borrelli, “Learning model predictive control for iterative tasks. a data-driven control framework,” _IEEE Transactions on Automatic Control_ , vol. 63, no. 7, pp. 1883–1896, 2017.
  * [20] U. Rosolia and F. Borrelli, “Minimum time learning model predictive control,” _International Journal of Robust and Nonlinear Control_ , 2020.
  * [21] M. Shahbazi, R. Babuška, and G. A. Lopes, “Unified modeling and control of walking and running on the spring-loaded inverted pendulum,” _IEEE Transactions on Robotics_ , vol. 32, no. 5, pp. 1178–1195, 2016.
  * [22] A. Wächter and L. T. Biegler, “On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming,” _Mathematical programming_ , vol. 106, no. 1, pp. 25–57, 2006.
  * [23] J. A. Andersson, J. Gillis, G. Horn, J. B. Rawlings, and M. Diehl, “Casadi: a software framework for nonlinear optimization and optimal control,” _Mathematical Programming Computation_ , vol. 11, no. 1, pp. 1–36, 2019.

Figure 4: Closed-loop trajectory computed from Algorithm 2 after $40$ updates
of the control policy.