# Isospin effect on baryon and charge fluctuations from the pNJL model

He Liu  Shanghai Institute of Applied Physics, Chinese Academy of Sciences,
Shanghai 201800, China  University of Chinese Academy of Sciences, Beijing
100049, China  Jun Xu111correspond author<EMAIL_ADDRESS>Shanghai
Institute of Applied Physics, Chinese Academy of Sciences, Shanghai 201800,
China

###### Abstract

We have studied the possible isospin corrections on the skewness and kurtosis
of net-baryon and net-charge fluctuations in the isospin asymmetric matter
formed in Au+Au collisions at RHIC-BES energies, based on a 3-flavor polyakov-
looped Nambu-Jona-Lasinio model. With typical scalar-isovector and vector-
isovector couplings leading to the splitting of $u$ and $d$ quark chiral phase
transition boundaries and critical points, we have observed dramatic isospin
effects on the susceptibilities, especially those of net-charge fluctuations.
Reliable experimental measurements at even lower collision energies are
encouraged to confirm the observed isospin effects.

###### pacs:

25.75.Nq, 12.39.-x, 25.75.Gz

Studying the hadron-quark phase transition and exploring the phase structure
of quantum chromodynamics (QCD) matter are the fundamental goals of
relativistic heavy-ion collision experiments. Lattice QCD simulations predict
that the transition between the hadronic phase and the partonic phase is a
smooth crossover at nearly zero baryon chemical potential ($\mu_{B}\sim 0$)
Ber05 ; Aok06 ; Baz12a . At larger $\mu_{B}$, the transition can be a first-
order one based on investigations from theoretical models, see, e.g., studies
from the Nambu-Jona-Lasinio (NJL) model and its extensions Asa89 ; Fuk08 ;
Car10 ; Bra13 . The knowledge of the QCD critical point (CP) in-between the
smooth crossover and the first-order phase transition boundaries is important
in mapping out the whole QCD phase diagram Ste06 . In order to find the
signature of the QCD CP at finite baryon chemical potentials, the Beam Energy
Scan (BES) program at the Relativistic Heavy-Ion Collider (RHIC) MMA10  has
been carrying out $``$low-energy$"$ relativistic heavy-ion collisions and
obtained many interesting results. However, heavy-ion collisions with neutron-
rich beams produce isospin asymmetric quark matter consisting of different net
numbers of $u$ and $d$ quarks, and the isospin degree of freedom is expected
to be increasingly important at lower collision energies with larger net
baryon densities, related to the QCD phase structure at finite isospin
chemical potentials $\mu_{I}$ Son01 . Although lattice QCD can explore the
phase transition at finite $\mu_{I}$ Kog02 ; Det12 ; Bra16 , it suffers from
the fermion sign problem at finite $\mu_{B}$ Bar86 ; Kar02 ; Mur03 . Based on
the NJL model studies, the isovector couplings may lead to the isospin
splittings of chiral phase transition boundaries and critical points for $u$
and $d$ quarks in isospin asymmetric quark matter Fra03 ; Tou03 ; Zha14 ;
Liu16 . This may influence the search of QCD CP signals in heavy-ion
collisions at RHIC-BES energies.

It was proposed Ste09  that the non-Gaussian fluctuations of observables in
relativistic heavy-ion collisions can be a signal of the phase transition
critical point, and the moment of these fluctuations can be used to measure
the magnitude of the correlation length. Typically, the susceptibilities of
conserved quantities carry information of the QCD phase boundary as well as
the position of the critical point Asa09 . The above findings have stimulated
measurements of the kurtosis and skewness of net-baryon, net-charge, and net-
strangeness fluctuations from $\sqrt{s_{NN}}=7.7$ to $200$ GeV at RHIC Agg10 ;
Ada14 ; Ada14c ; Ada16 ; Ada17 , as well as lattice QCD calculations Baz12 ;
Bor13 . For recent reviews on this topic, we refer the reader to Refs. Asa16 ;
Luo17 . However, to the best of our knowledge, the isospin effect on these
fluctuations was assumed to be small and has never been addressed. Actually,
based on studies from the statistic model And09 ; Sta14 ; Hat16 , the isospin
chemical potential can be as larger as 10 MeV at the chemical freeze-out stage
of heavy-ion collisions at RHIC-BES energies. Taking this isospin effect into
consideration, it is expected that the isospin splitting of the chiral phase
transition boundaries and the critical points may affect the moments of
fluctuations, especially for the net-charge fluctuations. In the present
study, we investigate the isospin effects based on the 3-flavor polyakov-
looped NJL (pNJL) model with scalar-isovector and vector-isovector couplings
incorporated in Ref. Liu16 .

The thermodynamic potential of the 3-flavor pNJL model with scalar-isovector
and vector-isovector couplings at temperature $T$ can be expressed as Liu16

$\displaystyle\Omega_{\textrm{pNJL}}$ $\displaystyle=$
$\displaystyle\mathcal{U}(\Phi,\bar{\Phi},T)-2N_{c}\sum_{i=u,d,s}\int_{0}^{\Lambda}\frac{d^{3}p}{(2\pi)^{3}}E_{i}$
(1) $\displaystyle-$ $\displaystyle
2T\sum_{i=u,d,s}\int\frac{d^{3}p}{(2\pi)^{3}}[\ln(1+e^{-3\beta(E_{i}-\tilde{\mu}_{i})}$
$\displaystyle+$ $\displaystyle 3\Phi
e^{-\beta(E_{i}-\tilde{\mu}_{i})}+3\bar{\Phi}e^{-2\beta(E_{i}-\tilde{\mu}_{i})})$
$\displaystyle+$
$\displaystyle\ln(1+e^{-3\beta(E_{i}+\tilde{\mu}_{i})}+3\bar{\Phi}e^{-\beta(E_{i}+\tilde{\mu}_{i})}$
$\displaystyle+$ $\displaystyle 3\Phi
e^{-2\beta(E_{i}+\tilde{\mu}_{i})})]+G_{S}(\sigma_{u}^{2}+\sigma_{d}^{2}+\sigma_{s}^{2})$
$\displaystyle-$ $\displaystyle
4K\sigma_{u}\sigma_{d}\sigma_{s}+G_{V}(\rho_{u}^{2}+\rho_{d}^{2}+\rho_{s}^{3})$
$\displaystyle+$ $\displaystyle
G_{IS}(\sigma_{u}-\sigma_{d})^{2}+G_{IV}(\rho_{u}-\rho_{d})^{2}.$

In the above, the temperature-dependent effective potential
$\mathcal{U}(\Phi,\bar{\Phi},T)$ as a function of the polyakov loop $\Phi$ and
$\bar{\Phi}$ is expressed as Fuk08

$\displaystyle\mathcal{U}(\Phi,\bar{\Phi},T)$ $\displaystyle=$
$\displaystyle-b\cdot T\\{54e^{-a/T\Phi\bar{\Phi}}+\ln[1-6\Phi\bar{\Phi}$ (2)
$\displaystyle-$ $\displaystyle
3(\Phi\bar{\Phi})^{2}+4(\Phi^{3}+\bar{\Phi}^{3})]\\},$

with $a=664$ MeV and $b$ = 0.015$\Lambda^{3}$ Fuk08 , where $\Lambda=750$ MeV
is the cutoff value in the momentum integral of the second term in Eg. (1).
The factor $2N_{c}$ with $N_{c}=3$ represents the spin and color degeneracy,
and $\beta=1/T$ represents the temperature. $G_{S}$ and $G_{V}$ are
respectively the scalar and vector coupling constants, $K$ is the coupling
constant of the six-point Kobayashi-Maskawa-t Hooft (KMT) interaction, and
$G_{IS}$ and $G_{IV}$ are respectively the strength of the scalar-isovector
and vector-isovector couplings that break the SU(3) symmetry while keeping the
isospin symmetry. For the ease of discussions, we define $R_{IS}=G_{IS}/G_{S}$
and $R_{IV}=G_{IV}/G_{S}$ as the reduced strength of the scalar-isovector and
vector-isovector coupling. The energy $E_{i}$ of quarks with flavor $i$ is
expressed as $E_{i}(p)=\sqrt{p^{2}+M_{i}^{2}}$, where $M_{i}$ is the
constituent quark mass. In the mean-field approximation, quarks can be
considered as quasiparticles with constituent masses $M_{i}$ determined by the
gap equation

$\displaystyle M_{i}$ $\displaystyle=$ $\displaystyle
m_{i}-2G_{S}\sigma_{i}+2K\sigma_{j}\sigma_{k}-2G_{IS}\tau_{3i}(\sigma_{u}-\sigma_{d}),$
(3)

where $m_{i}$ is the current quark mass, $\sigma_{i}$ stands for quark
condensate, ($i$, $j$, $k$) is any permutation of ($u$, $d$, $s$), and
$\tau_{3i}$ is the isospin quantum number of quark flavor $i$, i.e.,
$\tau_{3u}=1$, $\tau_{3d}=-1$, and $\tau_{3s}=0$. As shown in Eq. (3),
$\sigma_{d}$ and $\sigma_{s}$ contribute to the constituent quark mass of $u$
quarks as a result of the six-point interaction and the scalar-isovector
coupling. Similarly, the effective chemical potential expressed as

$\displaystyle\tilde{\mu}_{i}$ $\displaystyle=$
$\displaystyle\mu_{i}+2G_{V}\rho_{i}+2G_{IV}\tau_{3i}(\rho_{u}-\rho_{d})$ (4)

has also the contribution from quarks of other isospin states through the
vector-isovector coupling. The net quark number density of flavor $i$ can be
calculated from

$\displaystyle\rho_{i}=2N_{c}\int(f_{i}-\bar{f_{i}})\frac{d^{3}p}{(2\pi)^{3}}$
(5)

where

$\displaystyle
f_{i}=\frac{1+2\bar{\Phi}\xi_{i}+\Phi\xi_{i}^{2}}{1+3\bar{\Phi}\xi_{i}+3\Phi\xi_{i}^{2}+\xi_{i}^{3}}$
(6)

and

$\displaystyle\bar{f_{i}}=\frac{1+2\Phi{\xi}^{\prime}_{i}+\bar{\Phi}{{\xi}^{\prime}_{i}}^{2}}{1+3\Phi{\xi}^{\prime}_{i}+3\bar{\Phi}{{\xi}^{\prime}_{i}}^{2}+{{\xi}^{\prime}_{i}}^{3}}$
(7)

are the effective phase-space distribution functions for quarks and antiquarks
in the pNJL model with $\xi_{i}=e^{(E_{i}-\tilde{\mu_{i}})/T}$ and
${\xi}^{\prime}_{i}=e^{(E_{i}+\tilde{\mu_{i}})/T}$.

In the present study, we adopt the values of parameters Bra13 ; Lut92  as
$m_{u}=m_{d}=3.6$ MeV, $m_{s}=87$ MeV, $G_{S}\Lambda^{2}=3.6$, and
$K\Lambda^{5}=8.9$. Since the purpose is not to study the effect of $G_{V}$ on
the structure of phase diagram Asa89 ; Fuk08 ; Car10 ; Bra13 , it is set to 0
in the present study. The following equations

$\displaystyle\frac{\partial\Omega_{\textrm{pNJL}}}{\partial\sigma_{u}}=\frac{\partial\Omega_{\textrm{pNJL}}}{\partial\sigma_{d}}=\frac{\partial\Omega_{\textrm{pNJL}}}{\partial\sigma_{s}}=\frac{\partial\Omega_{\textrm{pNJL}}}{\partial\Phi}=\frac{\partial\Omega_{\textrm{pNJL}}}{\partial\bar{\Phi}}=0$

are used to obtain the values of $\sigma_{u}$, $\sigma_{d}$, $\sigma_{s}$,
$\Phi$, and $\bar{\Phi}$ at the minimum thermodynamic potential in the pNJL
model.

In the following, we consider fluctuation moments of conserved quantities,
such as net-baryon and net-charge fluctuations, from the above 3-flavor pNJL
model. The $n$-order susceptibility representing the cumulant of a given
conserved quantity in the grand ensemble can be expressed as the derivative of
the thermodynamic potential as

$\displaystyle\chi^{(n)}_{X}=\frac{\partial^{n}(-\Omega/T)}{\partial(\mu_{X}/T)^{n}},$
(9)

where $\mu_{X}$ represents the baryon ($\mu_{B}$) or the charge ($\mu_{Q}$)
chemical potential. Numerically, the isospin chemical potential $\mu_{I}$ and
the charge chemical potential $\mu_{Q}$ are equal to each other Baz12 ; Che15
. In the present study, we use the empirical relation between the isospin
chemical potential and the baryon chemical potential, i.e.,
$\mu_{I}=-0.293-0.0264\mu_{B}$ (MeV), determined from the statistical model
fits of the experimental data from Au+Au collisions from center-of-mass energy
$\sqrt{s_{NN}}$ = 7.7 GeV to 200 GeV And09 ; Sta14 ; Hat16 . For the
strangeness chemical potential, we take the empirical relation as
$\mu_{S}=1.032+0.232\mu_{B}$ (MeV) And09 ; Sta14 ; Hat16 . Note that the
chemical potentials of $u$, $d$, and $s$ quarks in the 3-flavor pNJL model can
be expressed in terms of $\mu_{B}$, $\mu_{I}$, and $\mu_{S}$. The higher-order
susceptibilities are related to the skewness $S$ and kurtosis $\kappa$
measured experimentally in relativistic heavy-ion collisions through the
relations

$\displaystyle
S\sigma=\frac{\chi^{(3)}}{\chi^{(2)}},\quad\kappa\sigma^{2}=\frac{\chi^{(4)}}{\chi^{(2)}},$
(10)

where $\sigma$ is the variance of the conserved quantity. The subscript of the
net baryon (B) or the net charge (Q) is omitted in the above equation.

Figure 1: (color online) Skewness and kurtosis of net-baryon (B) (left) and
net-charge (Q) (right) fluctuations in the $\mu_{B}-T$ plane with finite
scalar-isovector (upper panels) and vector-isovector (lower panels) coupling
constants and the empirical relation $\mu_{I}=-0.293-0.0264\mu_{B}$ (MeV). The
chiral phase transition boundaries and the corresponding critical points (CP)
of $u$ (black) and $d$ (white) quarks are plotted in all panels for reference.

We begin our discussion with the higher-order fluctuations of net baryon and
net charge from the 3-flavor pNJL model in the $\mu_{B}-T$ plane with various
isovector coupling constants, as shown in Fig. 1. As shown in Fig. 3 of Ref.
Liu16 , $R_{IS}=0.14$ and $R_{IV}=0.5$ lead to the splittings of $u$ and $d$
quark chiral phase transition boundaries as well as their critical points,
which are plotted in all panels of Fig. 1 for references. We note that the
empirical relation $\mu_{I}=-0.293-0.0264\mu_{B}$ (MeV) is only well valid
near the phase boundary, and this is exactly the region where we are
interested. As seen in Fig. 1, chiral phase transition boundaries separate the
red and blue areas for skewness results, representing respectively the
positive and negative values of $S\sigma$. For kurtosis results, however, the
chiral phase transition boundaries go through the blue areas, and the critical
points stand at the ends of the blue areas. It is also interesting to see that
the skewness of net-charge fluctuations gives the different orders of the red
and blue areas compared to that of net-baryon fluctuations. Using the
empirical relation between $\mu_{I}$ and $\mu_{B}$, the splitting of the
chiral phase transition boundaries from $R_{IV}=0.5$ is much weaker compared
to that observed in Ref. Liu16 .

Figure 2: (color online) Density plot of $S\sigma$ (left) and
$\kappa\sigma^{2}$ (right) in the $\mu_{B}-T$ plane as well as those along the
different hypothetical chemical freeze-out lines with the scalar-isovector
coupling constant $R_{IS}=0.14$ and the empirical relation
$\mu_{I}=-0.293-0.0264\mu_{B}$ (MeV). The solid line represents the averaged
chiral phase transition boundary for $u$ and $d$ quarks, while the other lines
are hypothetical chemical freeze-out curves by rescaling $\mu_{B}$ of the
solid line with factors of $0.98$, $0.95$, and $0.90$. The upper panels are
the net-baryon susceptibilities while the lower panels are the net-charge
susceptibilities.

In relativistic heavy-ion collision experiments, the net-baryon and net-charge
fluctuations are measured at chemical freeze-out. It is well-known that the
phase transition temperature from lattice QCD calculations at zero baryon
density is the same as that extracted from the statistical model at top RHIC
energy or LHC energy, leading to the conclusion that the chemical freeze-out
happens right after hadron-quark phase transition in extreme relativistic
heavy-ion collisions. In collisions at RHIC-BES energies, when the chemical
freeze-out happens is not known. There are several empirical criteria for
chemical freeze-out in relativistic heavy-ion collisions, such as fixed energy
per particle, baryon+antibaryon density, normalized entropy density, as well
as percolation model and so on (see Ref. Cle06  and references therein). In
order to compare qualitatively the higher-order susceptibilities from the pNJL
model with experimental results, we obtain the hypothetical chemical freeze-
out lines by rescaling $\mu_{B}$ of the averaged chiral phase transition
boundaries of $u$ and $d$ quarks with factors of $0.98$, $0.95$, and $0.90$,
corresponding respectively to the dash-dotted, dashed, and dotted curves in
Fig. 2. We note that a similar assumption of the chemical freeze-out lines was
made in Ref. Che17 , and the present hypothetical chemical freeze-out lines
are always below the chiral phase transition boundaries of both $u$ and $d$
quarks. For the net-baryon susceptibility shown in the upper panels of Fig. 2,
it is seen that $S\sigma(B)$ has one positive peak while $\kappa\sigma^{2}(B)$
has a positive and a negative peak along the chemical freeze-out lines, and
the peaks are sharper if the hypothetical chemical freeze-out lines are closer
to the chiral phase transition boundary. The critical point of the $d$ quark
chiral phase transition is always at the low-temperature or low-energy side of
the positive peaks for both $S\sigma(B)$ and $\kappa\sigma^{2}(B)$, and the
distance between the positive peaks and the critical point is related to that
between the chemical freeze-out line and the chiral phase transition boundary.
For the net-charge susceptibility shown in the lower panels of Fig. 2, we have
also observed sharper peaks if the hypothetical chemical freeze-out line is
closer to the chiral phase transition boundary. It is seen that $S\sigma(Q)$
has a negative peak at lower temperatures/energies and a positive peak at
higher temperatures/energies, and the broad positive peak of
$\kappa\sigma^{2}(Q)$ turns to two positive peaks and a negative one if the
hypothetical chemical freeze-out line is very close to the chiral phase
transition boundary. Again, the critical point of the $d$ quark chiral phase
transition is at the low-temperature or low-energy side of the negative peak
for $S\sigma(Q)$ or the positive peak for $\kappa\sigma^{2}(Q)$.

Figure 3: (color online) $S\sigma$ and $\kappa\sigma^{2}$ for net-baryon
(left) and net-charge (right) fluctuations along the closest hypothetical
chemical freeze-out line to the phase boundary with four different scenarios
of isospin chemical potentials and isovector couplings.

In order to compare the higher-order susceptibilities with and without the
isospin effect, we display in Fig. 3 the skewness and kurtosis of net-baryon
and net-charge fluctuations with four different scenarios of isospin chemical
potentials and isovector couplings: $\mu_{I}=0$, $R_{IS}=0$, $R_{IV}=0$;
$\mu_{I}=-0.293-0.0264\mu_{B}$ (MeV), $R_{IS}=0$, $R_{IV}=0$;
$\mu_{I}=-0.293-0.0264\mu_{B}$ (MeV), $R_{IS}=0.14$, $R_{IV}=0$;
$\mu_{I}=-0.293-0.0264\mu_{B}$ (MeV), $R_{IS}=0$, $R_{IV}=0.5$. In order to
illustrate the largest possible isospin effect, the susceptibilities are
calculated along the closest hypothetical chemical freeze-out line to the
phase boundary, i.e., using the rescaling factor of 0.98. It is seen that even
without isovector couplings, the skewness and kurtosis can be slightly
different for $\mu_{I}=0$ and $\mu_{I}\neq 0$, especially for net-charge
susceptibilities. With finite isovector coupling constants, the isospin effect
is largely enhanced. The peaks of the net-baryon susceptibilities move to the
low-temperature or low-energy side, especially for $R_{IS}=0.14$. The general
shape of the net-baryon susceptibility is qualitatively consistent with the
experimental data Ada14 ; Luo . The net-baryon susceptibility is not a unique
probe of the isospin effect, since it is largely affected by other effects as
well, such as the vector coupling. On the other hand, the isospin couplings
affect dramatically the net-charge susceptibilities. It is seen that the
isovector couplings lead to a negative peak at lower temperatures/energies for
$S\sigma(Q)$. For $\kappa\sigma^{2}(Q)$, the isovector couplings lead to two
peaks and $R_{IS}=0.14$ is the only scenario that leads to negative values.
The experimental results from STAR and PHENIX Collaborations for net-charge
susceptibilities are not consistent with each other yet Ada14c ; Ada16 . So
far the experimental results for $S\sigma(Q)$ seem to be positive above
$\sqrt{s_{NN}}=7.7$ GeV from both STAR and PHENIX measurements, and it is of
great interest to distinguish different scenarios if reliable measurements are
done at even lower collision energies. For the $\kappa\sigma^{2}(Q)$ results,
the STAR results lead to negative values at lower $\sqrt{s_{NN}}$ Ada14c
while the PHENIX results remain positive at all energies Ada16 . It is again
of great interest to check experimentally whether another peak appears at even
lower collision energies, as seen from the $\kappa\sigma^{2}(Q)$ results with
$R_{IS}=0.14$.

To summarize, based on the 3-flavor polyakov-looped Nambu-Jona-Lasinio model
with the scalar-isovector and vector-isovector couplings, we have studied the
higher-order susceptibilities of net-baryon and net-charge fluctuations in
isospin asymmetric matter formed in Au+Au collisions at RHIC-BES energies.
Although the general features remain the same, the isovector couplings move
the peaks of the net-baryon susceptibilities to the low-temperature or low-
energy side. On the other hand, the shape of the net-charge susceptibilities
is largely changed with the isovector couplings, i.e., an additional negative
peak appear in the skewness results and two positive peaks peaks appear in the
kurtosis results along the hypothetical chemical freeze-out line, if it is
very close to the phase boundary. It is of great interest to confirm our
findings by measuring the net-charge susceptibility at even lower collision
energies, such as that in the future FAIR-CBM program. Such analysis may also
be helpful in extracting the strength of the isovector couplings or even the
information of the isospin dependence of the QCD phase diagram.

We acknowledge helpful communications with Xiao-Feng Luo. This work was
supported by the Major State Basic Research Development Program (973 Program)
of China under Contract Nos. 2015CB856904 and 2014CB845401, the National
Natural Science Foundation of China under Grant Nos. 11475243 and 11421505,
the ”100-talent plan” of Shanghai Institute of Applied Physics under Grant
Nos. Y290061011 and Y526011011 from the Chinese Academy of Sciences, and the
Shanghai Key Laboratory of Particle Physics and Cosmology under Grant No.
15DZ2272100.

## References

  * (1) C. Bernard et al., Phys. Rev. D 71, 034504 (2005).
  * (2) Y. Aoki et al., Nature 443, 675 (2006).
  * (3) A. Bazavov et al., Phys. Rev. D 85, 054503 (2012).
  * (4) M. Asakawa and K. Yazaki, Nucl. Phys. A 504, 668 (1989).
  * (5) K. Fukushima, Phys. Rev. D 77, 114028 (2008) [Erratum-ibid. D 78, 039902 (2008)].
  * (6) S. Carignano, D. Nickel, and M. Buballa, Phys. Rev. D 82, 054009 (2010).
  * (7) N. M. Bratovic, T. Hatsuda, and W. Weise, Phys. Lett. B 719, 131 (2013).
  * (8) M. A. Stephanov, PoS LAT 2006, 024 (2006).
  * (9) M. M. Aggarwal et al. (STAR Collaboration), arXiv: 1007.2613 [nucl-ex].
  * (10) D. T. Son and M. A. Stephanov, Phys. Rev. Lett. 86, 592 (2001).
  * (11) J. B. Kogut and D. K. Sinclair, Phys. Rev. D 66, 034505 (2002).
  * (12) W. Detmold, K. Orginos, and Z. Shi, Phys. Rev. D 86, 054507 (2012).
  * (13) B. B. Brandt and G. Endrodi, arXiv: 1611.06758 [hep-lat]
  * (14) I. Barbour et al., Nucl. Phys. B 275, 296 (1986).
  * (15) F. Karsch, Lect. Notes Phys. 583, 209 (2002).
  * (16) S. Muroya, A. Nakamura, C. Nonaka, and T. Takaishi, Prog. Theor. Phys. 110, 615 (2003).
  * (17) M. Frank, M. Buballa, and M. Oertel, Phys. Lett. B 562, 221 (2003).
  * (18) Z. Zhang and H. P. Su, Phys. Rev. D 80, 054020 (2014).
  * (19) D. Toublan and J. B. Kugot, Phys. Lett. B 564, 211 (2003).
  * (20) H. Liu, J. Xu, L. W. Chen, and K. J. Sun, Phys. Rev. D 94, 065032 (2016).
  * (21) M. A. Stephanov, Phys. Rev. Lett. 102, 032301 (2009).
  * (22) M. Asakawa et al., Phys. Rev. Lett. 103, 262301 (2009).
  * (23) M. M. Aggarwal et al. (STAR Collaboration), Phys. Rev. Lett. 105, 022302 (2010).
  * (24) L. Adamczyk et al. (STAR Collaboration), Phys. Rev. Lett. 112, 032302 (2014).
  * (25) L. Adamczyk et al. (STAR Collaboration), Phys. Rev. Lett. 113, 092301 (2014).
  * (26) A. Adare et al. (PHENIX Collaboration), Phys. Rev. C 93, 011901(R) (2016).
  * (27) L. Adamczyk et al. (STAR Collaboration), arXiv: 1709.00773 [nucl-ex]
  * (28) A. Bazavov et al., Phys. Rev. Lett. 109, 192302 (2012).
  * (29) S. Borsanyi et al., Phys. Rev. Lett. 111, 062005 (2013).
  * (30) M. Asakawa and M. Kitazawa, Prog. Part. Nucl. Phys. 90, 299 (2016).
  * (31) X. F. Luo and N. Xu, Nucl. Sci. Tech. 28, 112 (2017).
  * (32) A. Andronic, P. Braun-Munzinger, and J. Stachel, Phys. Lett. B 673, 142 (2009).
  * (33) J. Stachel, A. Andronic, P. Braun-Munzinger, and K. Redlich, J. Phys. Conf. Ser. 509, 012019 (2014).
  * (34) Y. Hatta, A. Monnai, and B. W. Xiao, Nucl. Phys. A 947, 155 (2016).
  * (35) M. F. M. Lutz, S. Klimt, and W. Weise, Nucl. Phys. A 542, 521 (1992).
  * (36) J. W. Chen et al., Phys. Rev. D 92, 054019 (2015).
  * (37) J. Cleymans, H. Oeschler, K. Redlich, and S. Wheaton, Phys. Rev. C 73, 034905 (2006).
  * (38) J. W. Chen et al., Phys. Rev. D 95, 014038 (2017).
  * (39) X. F. Luo, private communications.


# Model predictive control with dynamic move blocking

Valentina Breschi, Simone Formentin and Alberto Leva This project was
partially supported by the Italian Ministry of University and Research under
the PRIN’17 project “Data-driven learning of constrained control systems ”,
contract no. 2017J89ARP.V. Breschi is with Department of Electrical
Engineering, Eindhoven University of Technology, Eindhoven, Netherlands. S.
Formentin and A. Leva are with Dipartimento di Elettronica, Informatica e
Bioingegneria (DEIB), Politecnico di Milano, Via G. Ponzio 34/5, 20133 Milano,
Italy. Corresponding author<EMAIL_ADDRESS>
###### Abstract

Model Predictive Control (MPC) has proven to be a powerful tool for the
control of systems with constraints. Nonetheless, in many applications, a
major challenge arises, that is finding the optimal solution within a single
sampling instant to apply a receding-horizon policy. In such cases, many
suboptimal solutions have been proposed, among which the possibility of
“blocking” some moves a-priori. In this paper, we propose a dynamic approach
to move blocking, to exploit the solution already available at the previous
iteration, and we show not only that such an approach preserves asymptotic
stability, but also that the decrease of performance with respect to the ideal
solution can be theoretically bounded.

###### Index Terms:

model predictive control, move blocking, computational load

##  I Introduction

Model Predictive Control (MPC) is fundamental to handle constrained control
problems providing formal guarantees [2, 4]. MPC most frequently relies on the
Receding Horizon (RH) approach: an Optimisation Problem (OP) is solved at each
step to yield a vector of future controls, only the first one is applied, and
the whole procedure is repeated for the next step.

Sometimes, however, solving the OP within one control timestep can be
computationally infeasible, and simplifying for brevity, one can choose among
the following alternatives:

  * •

increase the timestep — but this can be prohibited by the physics of the
control problem;

  * •

resort to _explicit_ MPC [1] — but this can entail too demanding memory
requirements and would not be feasible for large scale problems;

  * •

simplify the OP, for example via local linearisations — but this makes global
properties hard to enforce [3];

  * •

reduce the OP dimensionality by shrinking the control horizon or by having
some future controls depend on others, which is called “Move Blocking” (MB)
[10, 5].

Notice that to perform less optimisations one could also apply not just the
first computed control sample as in the RH case, but up to all those that
cover the (future) control horizon $N$. This is called OL-MPC for Open-Loop
MPC, as the loop actually closes just at each $N$-th step, but does not remove
the constraint of solving the OP within one timestep.

In our research, we propose a novel combination of OL-MPC and MB, accepting
that the OP cannot be solved in one timestep and just requiring an overbound
$N_{OP}$ of the number of steps required. In detail

  * •

after the $k$-th optimisation we apply $N-N_{B}$ control samples, where
$N_{B}\geq N_{OP}$ is the _blocking horizon_ ,

  * •

we start the $(k+1)$-th optimisation $N_{OP}$ steps before the above samples
are exhausted, hence overlapping the control horizons,

  * •

and we constrain the first $N_{OP}$ controls from the $(k+1)$-th optimisation
to equal those coming from the $k$-th one, applied – as said – while the
$(k+1)$-th is running.

The main novelty is that in doing so we make the MB mechanism stateful, which
is why the name our technique DMB for Dynamic MB.

In this paper, analogously to what has been done for other (static) MB
strategies (see, e.g., [6, 9]), we analyse the DMB strategy so as to assess
the level of sub-optimality with respect to a purely RH realisation of the
same MPC scheme, that we take as reference but is infeasible owing to
computational limits. We further show that asymptotic stability is preserved
under some mild assumptions, and illustrate the limited decrease of
performance on a simple numerical example. A preliminary version of this idea
was proposed in [8], where however no theoretical analysis of the proposed
scheme was provided.

The remainder of the paper is as follows. In Section II, we formally state the
problem of interest and introduce the main terminology and notation. Section
III discusses the DMB-MPC approach and illustrates its theoretical properties.
A simple numerical example is proposed in Section IV. The paper is ended by
some concluding remarks.

##  II Problem statement

Consider the _nonlinear_ discrete-time dynamical system

$x(t+1)=f(x(t),u(t)),~{}~{}x(0)=x_{0},$ (1)

where $x(t)\in\mathcal{X}\subseteq\mathbb{R}^{n}$ is the state of the system
at time $t\in\mathbb{R}^{n}$, laying in a _finite dimensional_ space $X$,
while $u(t)\in U\subseteq\mathbb{R}^{m}$ is a controllable input and
$f:X\times U\rightarrow X$. For a given state $x$, our goal would be to find a
feedback control action $u=\mu_{\infty}(x)$ that minimizes the infinite
horizon cost

$J_{\infty}(x,u)=\sum_{k=0}^{\infty}\ell(x_{u}(k),u(k)),~{}~{}x_{u}(0)=x,$ (2)

where $x_{u}(k)$ denotes the state realized for some control
sequence111$\mathcal{U}$ indicates the space of feasible control sequences.
$u\in\mathcal{U}$, and $\ell:X\times U\rightarrow\mathbb{R}_{0}^{+}$ is the
running cost of the problem. The latter is here assumed to be _positive
definite_ and _proper_ with respect to some compact set $A$.

However, minimizing (2) is often infeasible in practice, motivating the
introduction of _model predictive control_ (MPC) strategies that focus on
solving the following finite-horizon alternative:

$\displaystyle\underset{u}{\mbox{minimize}}~{}~{}J_{N}(x,u)$ (3a)
$\displaystyle\qquad~{}\mbox{s.t.
}~{}x_{u}(k\\!+\\!1)\\!=\\!f(x_{u}(k),u(k)),~{}~{}\forall
k\in\mathcal{I}_{N},$ (3b) $\displaystyle\qquad\qquad~{}x_{u}(0)\\!=\\!x,$
(3c) $\displaystyle\qquad\qquad~{}x_{u}(k)\in X,~{}u(k)\in U,~{}~{}\forall
k\in\mathcal{I}_{N},$ (3d) where $\mathcal{I}_{N=\\{0,\ldots,N-1\\}}$ and the
cost is given by $J_{N}(x,u)=\sum_{k=0}^{N-1}\ell(x_{u}(k),u(k)).$ (3e)

Once the problem is solved, only the first control action is fed to the
system, whereas the rest of the optimized input sequence is discarded. The
problem in (3) needs then to be solved again, with updated initial state, at
the following sampling instant. By introducing the optimal value function
$V_{N}(x)$ associated to this problem, namely,

$V_{N}(x)=\inf_{u\in\mathcal{U}}J_{N}(x,u),$ (4)

a direct consequence of the Bellman’s optimality principle is that the control
law is finally selected as

$\mu_{N}(x)=\arg\min_{u\in U}~{}\\{V_{N-1}(f(x,u))+\ell(x,u)\\}.$ (5)

Nonetheless, especially when $N$ is large, also this simplified strategy can
eventually become unfeasible in practice, since it requires the solution of an
optimization problem in real-time at each time step. This is indeed the case
of fast sampled systems, where a single sampling instant might not be
sufficient to solve the entire optimization problem. Under the assumption that
a number of steps $N_{OP}(t)>1$ 222We stress here that $N_{OP}(t)$ might
potentially be time-varying. are needed to solve (3) at each time instant $t$,
our objective in this paper is to propose a move-blocking strategy that allows
us to $(i)$ approximate the solution of (3) with a certified level of sub-
optimality, while $(ii)$ being compatible with the computational limits of the
application at hand.

##  III Dynamic move blocking

In this work, we focus on two elements to reduce the computational burden
characterizing standard receding horizon strategies: $(i)$ the degrees of
freedom of (3) and $(ii)$ the input actually fed to the system.

To limit the degrees of freedom and, thus, the computational complexity of the
optimal control problem at time $t$, we leverage the idea of move-blocking
strategies and consider the following “variation” of the standard MPC problem:

$\displaystyle\underset{u}{\mbox{minimize}}~{}~{}J_{N}(x(t),u)$ (6a)
$\displaystyle~{}~{}\mbox{s.t.
}~{}x_{u}(k\\!+\\!1)\\!=\\!f(x_{u}(n),u(n)),~{}~{}\forall
k\in\mathcal{I}_{N},$ (6b) $\displaystyle\qquad~{}~{}x_{u}(0)\\!=\\!x(t),$
(6c)
$\displaystyle\qquad~{}~{}u(k)\\!=\\!\tilde{u}(N\\!\\!-\\!\\!N_{B}\\!+\\!k),~{}k\\!=\\!0,\ldots,N_{B}\\!-\\!1,$
(6d) $\displaystyle\qquad~{}~{}x_{u}(k)\in X,~{}u(n)\in U,~{}~{}\forall
k\in\mathcal{I}_{N},$ (6e) where $x(t)$ is the measured (or estimated) state
of the system at the time when the optimization is carried out, and $N_{B}\geq
N_{OP}$ dictates the number of moves that are _overlapped_ , since they are
here blocked to the optimal control actions obtained by solving this problem
at the previous optimization step (here denoted as $\tilde{u}\in\mathcal{U}$).
As the values of the “blocked” inputs depends on the optimization at the
previous time instant, we will call such a strategy MPC with dynamic move
blocking (MPC-DMB) hereafter.

######  Remark 1 ($N_{B}$ and $N$)

Clearly, the number $N_{B}$ of blocked actions must also satisfy $N_{B}\leq
N$. The special case of $N_{B}=N$ would imply that the past sequence is fully
used to control the system and only the action $u(t+N-1)$ is optimally
computed. As we will show next, this choice is far from being optimal, since
the last control action is computed without “looking into the predicted
future” of the state, ultimately jeopardizing closed-loop stability.

###  III-A Properties of MPC-DMB

We now discuss some of the properties of the control law originated by the
MPC-DMB approach, namely

$\mu_{N}^{\mathrm{DMB}}(x)=\begin{bmatrix}\tilde{u}(N-N_{B})\\\ \vdots\\\
\tilde{u}(N-1)\\\ \mu_{N_{B}}(x_{N_{B}})\end{bmatrix},$ (7a) where $x_{N_{B}}$
compactly denotes the $N_{B}$-step ahead prediction of the system’s state
starting from the initial condition $x$ and using the blocked inputs
$\tilde{u}(N-N_{B}+n)$, for $n=1,\ldots,N_{B}-1$, and
$\mu_{N_{B}}(x_{N_{B}})\\!=\\!\arg\min_{u\in
U}\\{V_{N\\!-\\!N_{B}\\!-\\!1}(f(x_{N_{B}},u))\\!+\\!\ell(x_{N_{B}},u)\\}.$
(7b)

In the next Lemma, we will show that blocking the first $N_{B}$ actions in (6)
is equivalent to reformulate the problem as a standard MPC with reduced
horizon

$\displaystyle\underset{u}{\mbox{minimize}}~{}~{}J_{N-N_{B}}(x(t+N_{B}),u)$
(8a) $\displaystyle~{}~{}\mbox{s.t.
}~{}x_{u}(k\\!+\\!1)\\!=\\!f(x_{u}(k),u(k)),~{}~{}\forall
k\in\mathcal{I}_{N-N_{B}},$ (8b)
$\displaystyle\qquad~{}~{}x_{u}(0)\\!=\\!x(t+N_{B}),$ (8c)
$\displaystyle\qquad~{}~{}x_{u}(k)\in X,~{}u(k)\in U,~{}~{}\forall
k\in\mathcal{I}_{N-N_{B}},$ (8d) where
$\mathcal{I}_{N-N_{B}}=\\{0,\ldots,N-N_{B}\\}$, the cost is defined as
$J_{N-N_{B}}(x(t+N_{B}),u)=\\!\\!\\!\\!\\!\\!\sum_{k=0}^{N-N_{B}-1}\\!\\!\\!\ell(x_{u}(k),u(k)),$
(8e)

and the initial state is replaced by a suitable _prediction_ of the state at
time $t+N_{B}$ using the model equation:

$x(t+N_{B})=f(x(t+N_{B}-1),\tilde{u}(N-1)).$ (9)

######  Lemma 1 (Equivalence with reduced-horizon MPC)

Solving (6) is equivalent to tackling the reduced horizon problem (8), with
initial condition given by the predicted state $x(t+N_{B})$ defined in (9).

###### Proof:

Since the first $N_{B}$ actions are blocked, the cost in (6a) can be rewritten
as

$J_{N}(x(t),u)=J_{N_{B}}(x(t),\tilde{u})+\\!\\!\sum_{k=N_{B}}^{N-1}\\!\\!\ell(x_{u}(k),u(k)),$
(10a) where
$J_{N_{B}}(x(t),\tilde{u})=\sum_{k=0}^{N_{B}-1}\ell(x_{u}(k),\tilde{u}(N\\!-\\!N_{B}+k)),$
(10b)

and $x_{u}(k)$, for $k=1,\ldots,N_{B}$ can be computed based on the predictive
model in (1), namely333Note that both the inputs used for prediction and the
associated state satisfy the constraints by design.

$\displaystyle x_{u}(1)=x(t+1)\\!=\\!f(x(t),\tilde{u}(N-N_{B})),$
$\displaystyle
x_{u}(k)\\!=\\!f(x_{u}(k-1),\tilde{u}(N\\!-\\!N_{B}\\!+\\!k\\!-\\!1)),~{}k\\!=\\!2,\ldots,N_{B},$

starting from $x_{u}(0)=x(t)$. Accordingly, $J_{N_{B}}(x(t),\tilde{u})$ is
fixed based on this pre-computable state sequence and, thus, it does not
impact on the solution of (6). At the same time, due to the dynamics of the
controlled system, the only state that directly impacts the computation of the
optimal input sequence is $x_{u}(N_{B})$. Based on these considerations, the
optimization problem in (6) can be equivalently recast as

$\displaystyle\underset{u}{\mbox{minimize}}~{}~{}\sum_{k=N_{B}}^{N-1}\ell(x_{u}(k),u(k))$
(11a) $\displaystyle\mbox{s.t.
}x_{u}(k\\!+\\!1)\\!=\\!f(x_{u}(k),u(k)),~{}k\\!=\\!N_{B},\ldots,N\\!-\\!1,$
(11b) $\displaystyle\quad~{}x_{u}(N_{B})\\!=\\!x(t+N_{B}),$ (11c)
$\displaystyle\quad~{}x_{u}(k)\in X,~{}u(k)\in U,~{}~{}k=N_{B},\ldots,N-1.$
(11d)

The reduced-horizon MPC problem in (8) directly stems from this result based
on a simple change of indexes, _i.e.,_ considering a new index $k=n-N_{B}$ and
shifting the limits of the sum in the loss and the constraints accordingly. ∎

Algorithm 1 MPC-DMB at time $t$

Input: initial condition $x(t)$; previous optimal sequence $t-1$,
$\\{\tilde{u}(k)\\}_{k=0}^{N-1}$.

  1. 1.

for $k=0,\ldots,N_{B}-1$

    1. 1.

predict $x(t+k+1)\\!=\\!f(x(t+k),\tilde{u}(N\\!-\\!N_{B}+k))$;

  2. 2.

solve the blocked MPC problem in (8);

  3. 3.

extract the first optimal action $\mu_{N_{B}}(x(t+N_{B}))\\!=\\!u(0)$;

  4. 4.

retrieve $x(t+N_{B})$;

  5. 5.

apply

$u(t+k)=\begin{cases}\tilde{u}(N\\!-\\!N_{B}\\!+\\!k),&\mbox{if }0\leq k\leq
N_{B}\\!-\\!1,\\\ \mu_{N_{B}}(x(t+N_{B})),&\mbox{if }k=N_{B};\end{cases}$

  6. 6.

update $\tilde{u}=\\{u(t+k)\\}_{k=0}^{N_{B}}$;

  7. 7.

set $x(0)=x(t+N_{B})$;

  8. 8.

shift to $t+N_{B}$;

Once the solution of (6) is computed, the first $N_{B}$ samples of the control
sequence (namely, those fixed based on the results of the previous
optimization step and the first optimal action over the reduced horizon
$N-N_{B}$) are kept, while the rest of the optimal sequence is discarded and
the optimization window shifts of $N_{B}$ steps. The resulting MPC-DMB is
summarized in Algorithm 1.

Let us define the optimal value function associated to (7a) as

$V_{N}^{\mathrm{DMB}}(x)=\delta(x)+V_{N-N_{B}}(x_{N_{B}}),$ (12a) where
$\displaystyle V_{N\\!-\\!N_{B}}(x_{N_{B}})\\!=\\!\min_{u\in
U}\\{V_{N\\!-\\!N_{B}\\!-\\!1}(f(x_{N_{B}},u))\\!+\\!\ell(x_{N_{B}},u)\\},$
(12b)
$\displaystyle\delta(x)=\sum_{k=0}^{N_{B}-1}\ell(x_{\tilde{u}}(k),\tilde{u}(N-N_{B}+k)),$
(12c) with $x_{\tilde{u}}(0)=x$.

By comparing the value functions of the receding horizon strategy (4) and the
dynamic move blocking approach, we can formalize the following result.

######  Lemma 2 (Sub-optimality of MPC-DMB)

Consider a receding horizon $N\in\mathbb{N}$ and a blocking horizon
$N_{B}\in\mathbb{N}$, with $1\leq N_{B}\leq N$. Then the following holds:

$V_{N}(x)\leq V_{N}^{\mathrm{DMB}}(x),\ \forall x\in X.$ (13)

###### Proof:

According to Bellman’s optimality principle, the value function $V_{N}(x)$ can
be rewritten as

$V_{N}(x)=\min_{u\in U}\\{V_{N-1}(f(x,u))+\ell(x,u)\\}.$ (14)

Therefore, by picking any $\mu(x)$ different from the receding-horizon
controller $\mu_{N}(x)$ in (5), the following holds by definition:

$V_{N}(x)\\!\leq\\!\ell(x,\mu(x))+\min_{u\in
U}\\{V_{N\\!-\\!2}(f(x_{\mu},u))+\ell(x_{\mu},u)\\},$ (15)

where $x_{\mu}=f(x,\mu(x))$. As a particular case, the previous inequality is
thus verified when setting $\mu(x)=\tilde{u}(N-N_{B})$, _i.e.,_ considering
the first control action resulting from the dynamic move blocking strategy.
Analogously, by considering again $\mu(x_{\mu})$ different from
$\mu_{N-1}(x_{\mu})$, we further have that

$\displaystyle V_{N}(x)\leq$
$\displaystyle\sum_{k=0}^{1}\ell(x_{\mu}(k),\mu(x_{\mu}(k)))+$
$\displaystyle~{}~{}+\min_{u\in
U}\\{V_{N\\!-\\!3}(f(x_{\mu}(2),u))\\!+\\!\ell(\tilde{x}_{\mu}(2),u)\\},$ (16)

where

$\displaystyle x_{\mu}(0)=x,$ $\displaystyle
x_{\mu}(k)=f(x_{\mu}(k-1),\mu(x_{\mu}(k-1))),~{}~{}k=1,\ldots,2.$

Once again, specifically $\mu(x_{\mu})=\tilde{u}(N-N_{B}+1)$ is only a
particular choice for the sub-optimal law and, thus, the previous inequality
holds when using the dynamic move blocking strategy. Exploiting the same
reasoning iteratively, it is straightforward to prove that

$\displaystyle
V_{N}(x)\leq\\!\\!\\!\sum_{k=0}^{N-N_{B}}\ell(x_{\mu}(k),\mu(x_{\mu}(k)))+$
$\displaystyle~{}~{}+\min_{u\in
U}\\{V_{N\\!-\\!N_{B}\\!-\\!1}(f(x_{\mu}(N_{B}),u))\\!+\\!\ell(\tilde{x}_{\mu}(N_{B}),u)\\},$
(17)

with

$\displaystyle x_{\mu}(0)=x,$ $\displaystyle
x_{\mu}(k)=f(x_{\mu}(k-1),\mu(x_{\mu}(k-1))),~{}~{}k\\!=\\!1,\ldots,N-N_{B}.$

The proof is thus concluded by simply replacing $\mu(x_{\mu}(k))$ with the
$k$-th component $[\mu_{N}^{\mathrm{DMB}}(x)]_{k}$ of
$\mu_{N}^{\mathrm{DMB}}(x)$ in (7a), for $k=0,\ldots,N_{B}-1$. ∎

######  Remark 2 (The role of $N_{B}$ for suboptimality)

Let us consider two possible choices for $N_{B}$, denoted as $N_{B}^{(1)}$ and
$N_{B}^{(2)}$ respectively. Let $N_{B}^{(1)}<N_{B}^{(2)}<N$. As a
straightforward consequence of Bellman’s optimality principle, it holds that

$V_{N}^{\mathrm{DMB},(1)}(x)\leq V_{N}^{\mathrm{DMB},(2)}(x)$ (18)

where $V_{N}^{\mathrm{DMB},(i)}(x)$ corresponds to $V_{N}^{\mathrm{DMB}}(x)$
with $N_{B}=N_{B}^{(i)}$, for $i=1,2$, This can be easily proven by following
the same reasoning exploited in the proof of Lemma 2.

Let us now introduce the “ideal” value function

$V_{\infty}(x_{0})=\inf_{u\in\mathcal{U}}J_{\infty}(x(0),u),$ (19)

which one would attain by minimizing the (practically unfeasible) infinite-
horizon loss in (2), and let us denote with $V_{\infty}^{\mu_{N}}(x)$ and
$V_{\infty}^{\mu_{N}^{\mathrm{DMB}}}(x)$ the optimal value functions obtained
when using the receding-horizon controller and the dynamic move blocking
strategy over an infinite horizon, respectively, namely,

$\displaystyle
V_{\infty}^{\mu_{N}}(x)=\sum_{k=0}^{\infty}\ell(x_{\mu_{N}}(k),\mu_{N}(x_{\mu_{N}}(k))),$
(20a) $\displaystyle
V_{\infty}^{\mu_{N}^{\mathrm{DMB}}}(x)=\delta(x)+\sum_{k=N_{B}}^{\infty}\ell(x_{\mu_{N}^{\mathrm{DMB}}}(k),\mu_{N}(x_{\mu_{N}^{\mathrm{DMB}}}(k)))=$
$\displaystyle~{}\qquad\qquad=\delta(x)+V_{\infty}^{\mu_{N_{B}}}(x_{\mu_{N}^{\mathrm{DMB}}}(N_{B})).$
(20b)

Due to the suboptimality of the receding-horizon controller and the MPC-DMB
law, these three functions satisfy the relationship

$V_{\infty}(x)\leq V_{\infty}^{\mu_{N}}(x)\leq
V_{\infty}^{\mu_{N}^{\mathrm{DMB}}}(x),$ (21)

as formalized in the following Lemma.

######  Lemma 3 (Ranking of infinite-horizon value functions)

Let $N,N_{B}\in\mathbb{N}$ be the prediction and blocking horizons
characterizing the problem in (6), with $N_{B}\geq 1$. Then the relationship
in (21) holds for all $x\in X$.

###### Proof:

Based on the definition of $\mu_{N}^{\mathrm{DMB}}(x)$ and its last component
$\mu_{N_{B}}(x)$, it holds that

$\displaystyle V_{N}(x)=V_{N-1}(f(x,\mu_{N}(x)))+\ell(x,\mu_{N}(x)),$ (22)
$\displaystyle
V_{N}^{\mathrm{DMB}}(x)\\!=\\!\delta(x)\\!+\\!V_{N\\!-\\!N_{B}\\!-\\!1}(f(x_{\mu^{\mathrm{DMB}}},\mu_{N_{B}}(x_{\mu^{\mathrm{DMB}}})))+$
$\displaystyle\qquad\qquad\qquad\qquad\qquad+\ell(x_{\mu^{\mathrm{DMB}}},\mu_{N_{B}}(x_{\mu^{\mathrm{DMB}}})).$
(23)

Since $N_{B}\neq 0$, $\delta(x)>0$, the MPC-DMB strategy is suboptimal with
respect to a standard receding-horizon and the inequality in (13) is
satisfied. Based on this result and on the previous definition, (21) follows
when considering $N\rightarrow\infty$. ∎

######  Remark 3 (Corollary of Lemma 3)

The suboptimality of the dynamic move blocking strategy is lower-bounded by
that of a standard MPC approach, namely

$\frac{V_{\infty}^{\mu_{N}}(x)-V_{\infty}}{V_{\infty}}\leq\frac{V_{\infty}^{\mu_{N}^{\mathrm{DMB}}}(x)-V_{\infty}}{V_{\infty}},\vspace{0.1cm}$
(24)

Starting from these results, we now follow [7] toward detecting an upper-bound
on the blocked horizon $N_{B}$ preserving the properties characterizing the
standard receding horizon strategy and to characterize the suboptimality of
the proposed approach. To this end, let us recall the following proposition.

######  Proposition 1 ([7])

Consider a generic feedback law $\tilde{\mu:X\rightarrow U}$ and a function
$\tilde{V}:X\rightarrow\mathbb{R}_{0}^{+}$ satisfying the inequality

$\tilde{V}(x)\geq\tilde{V}(f(x,\tilde{\mu}(x)))+\alpha\ell(x,\tilde{\mu}(x)),$
(25)

for some $\alpha\in[0,1]$ and $x\in X$. Then, for all $x\in X$, the following
estimate holds:

$\alpha V_{\infty}\leq\alpha V_{\infty}^{\tilde{\mu}}(x)\leq\tilde{V}(x).$
(26)

By focusing on the reduced-horizon problem in (8), whose associated optimal
value function is $V_{N-N_{B}}(x_{N_{B}})$ in (12b), we introduce the
following result.

######  Lemma 4 (Infinite vs reduced horizon value functions)

Consider $N,N_{B}\in\mathbb{N}$, with $N_{B}\in[1,N]$, and the feedback law
$\mu_{N_{B}}(x_{N_{B}})$ defined as in (7b). Assume that

$\displaystyle V_{N\\!-\\!N_{B}}(x_{N_{B}}^{+})$ $\displaystyle-
V_{N\\!-\\!N_{B}\\!-\\!1}(x_{N_{B}}^{+})\leq$
$\displaystyle\qquad\qquad(1-\alpha)\ell(x_{N_{B}},\mu_{N_{B}}(x_{N_{B}}))),$
(27)

holds for some $\alpha\in[0,1]$ and all $x\in X$, with
$x_{N_{B}}^{+}=f(x_{N_{B}},\mu_{N_{B}}(x_{N_{B}})$. Then,
$V_{N\\!-\\!N_{B}}(x_{N_{B}})$ satisfies (25) and, thus

$\alpha V_{\infty}^{\mu_{N_{B}}}(x_{N_{B}})\leq V_{N\\!-\\!N_{B}}(x_{N_{B}}),$
(28)

for all $x\in X$.

###### Proof:

By exploiting the definition of $V_{N\\!-\\!N_{B}}(x_{N_{B}})$ in (12b) and
$\mu_{N_{B}}(x_{N_{B}})$ in (7b) it follows that

$\displaystyle V_{N\\!-\\!N_{B}}$
$\displaystyle(x_{N_{B}})=V_{N\\!-\\!N_{B}\\!-\\!1}(x_{N_{B}}^{+})+\ell(x_{N_{B}},\mu_{N_{B}}(x_{N_{B}}))$
$\displaystyle\geq
V_{N\\!-\\!N_{B}}(x_{N_{B}}^{+})\\!-\\!(1\\!-\\!\alpha)\ell(x_{N_{B}},\mu_{N_{B}}(x_{N_{B}}))+$
$\displaystyle\qquad\qquad\qquad\qquad\qquad+\\!\ell(x_{N_{B}},\mu_{N_{B}}(x_{N_{B}}))$
$\displaystyle=V_{N\\!-\\!N_{B}}(x_{N_{B}}^{+})+\alpha)\ell(x_{N_{B}},\mu_{N_{B}}(x_{N_{B}})),$

where the second inequality holds thanks to (4). Based on the definition of
$x_{N_{B}}^{+}$, from the previous inequality it follows that (25) holds,
directly yielding (28). ∎

Let us now introduce the following assumptios.

######  Assumption 1 ( [7])

For $N\in\mathbb{N}$, there exists $\gamma>0$ such that

$\displaystyle V_{2}(x)\leq(\gamma+1)V_{1}(x),$ (29) $\displaystyle
V_{k}(x)\leq(\gamma+1)\ell(x,\mu_{k}(x)),k=2,\ldots,N,$ (30)

holds for all $x\in X$.

Then, the relationship between $V_{N\\!-\\!N_{B}}(x_{N_{B}})$ and
$V_{N\\!-\\!N_{B}\\!-\\!1}(x_{N_{B}})$ can be characterized through $\gamma$
based on the following proposition.

######  Proposition 2 (Shrinking of $V_{N\\!-\\!N_{B}}(x_{N_{B}})$)

Let $N-N_{B}\geq 2$ and Assumption 1 be satisfied for the MPC problem (8) with
horizon $N-N_{B}$. Then, the inequality

$\eta(\gamma,N-N_{B})V_{N\\!-\\!N_{B}}(x_{N_{B}})\\!\leq\\!V_{N\\!-\\!N_{B}\\!-\\!1}(x_{N_{B}}),$
(31)

holds for all $x\in X$, where

$\eta(\gamma,N-N_{B})=\left(\\!\frac{(\gamma\\!+\\!1)^{N\\!-\\!N_{B}\\!-\\!2}}{(\gamma\\!+\\!1)^{N\\!-\\!N_{B}\\!-\\!2}+\gamma^{N\\!-\\!N_{B}}}\\!\right).\vspace{.2cm}$
(32)

The proof of this proposition straightforwardly follows from that of [7,
Proposition 4.4], by replacing the full horizon $N$ with the difference
$N-N_{B}$ and it is thus omitted.

######  Remark 4 (Constraining $N_{B}$)

Through Proposition 2, we obtain a first upper bound on $N_{B}$. Indeed, the
blocking horizon $N_{B}$ must satisfy

$N_{OP}\leq N_{B}\leq N-2,$ (33)

for the result in Proposition 2 to hold. In turn, this enforces $N$ to be at
least equal to $N_{OP}+2$.

######  Remark 5 (A note about Assumption 1)

The satisfaction of Assumption 1 can be checked by analyzing the running cost
$\ell(\cdot)$. In particular, as detailed in [7, Proposition 4.7], the loss
should satisfy

$\ell(x,u)\geq\alpha W(x),$

with $\alpha>0$ and $W:X\rightarrow\mathbb{R}_{0}^{+}$, while also verifying
an exponential controllability condition.

By leveraging this result, we can now characterize the suboptimality of
$\mu_{N_{B}}(x_{N_{B}})$ as follows.

######  Theorem 1 (Asymptotic stability)

Let $N,N_{B}\in\mathbb{N}$ be the prediction and blocking horizon of the MPC-
DMB problem in (6), respectively, and let $N_{B}$ satisfy (33). Let
$\gamma>0$, assume that

$(\gamma+1)^{N-N_{B}-2}>\gamma^{N-N_{B}},$ (34)

and that Assumption 1 holds for these $\gamma$ and for the reduced horizon
$N-N_{B}$. Then, for all $x\in X$, it holds that

$V_{\infty}^{\mu_{N_{B}}}(x_{N_{B}})\leq\frac{(\gamma+1)^{N\\!-\\!N_{B}\\!-\\!2}}{(\gamma+1)^{N\\!-\\!N_{B}\\!-\\!2}\\!-\\!\gamma^{N-N_{B}}}V_{\infty}(x_{N_{B}}),$
(35)

and, as a consequence,

$\frac{V_{\infty}^{\mu_{N_{B}}}(x_{N_{B}})\\!-\\!V_{\infty}(x_{N_{B}})}{V_{\infty}(x_{N_{B}})}\\!\leq\\!\frac{\gamma^{N\\!-\\!N_{B}}}{(\gamma+1)^{N\\!-\\!N_{B}\\!-\\!2}\\!-\\!\gamma^{N\\!-\\!N_{B}}}.\vspace{.1cm}$
(36)

###### Proof:

The proof follows the same steps of [7, Proof of Theorem 4.5]. In particular,
from proposition 2, it follows that

$\displaystyle\Delta V_{N\\!-\\!N_{B}}(x_{N_{B}})$
$\displaystyle\leq(\eta^{\\!-1}(\gamma,n-n_{b})-1)V_{N\\!-\\!N_{B}\\!-\\!1}(x_{N_{B}})$
$\displaystyle=\frac{\gamma^{N-N_{B}}}{(\gamma+1)^{N-N_{B}-2}}V_{N\\!-\\!N_{B}\\!-\\!1}(x_{N_{B}}).$

where $\Delta
V_{N\\!-\\!N_{B}}(x_{N_{B}})=V_{N\\!-\\!N_{B}}(x_{N_{B}})-V_{N\\!-\\!N_{B}\\!-\\!1}(x_{N_{B}})$.
Since Assumption 1 implies that444The reader is referred to [7] for a detail
proof of this inequality.

$\displaystyle\Delta
V_{N\\!-\\!N_{B}}(f(x_{N_{B}},\mu_{N_{B}}(x_{N_{B}}))\leq$
$\displaystyle\qquad\qquad\qquad\qquad\frac{\gamma^{N-N_{B}}}{(\gamma\\!+\\!1)^{N\\!-\\!N_{B}\\!-\\!2}}\ell(x_{N_{B}},\mu_{N_{B}}(x_{N_{B}})),$

in turn, yielding the inequality in (4). Therefore, Lemma 2 can be applied by
setting

$\alpha=1-\frac{\gamma^{N\\!-\\!N_{B}}}{(\gamma+1)^{N\\!-\\!N_{B}\\!-\\!2}}=\frac{(\gamma+1)^{N\\!-\\!N_{B}\\!-\\!2}-\gamma^{N\\!-\\!N_{B}}}{(\gamma+1)^{N\\!-\\!N_{B}\\!-\\!2}}.$
(37)

Hence, (35) results directly from (26) that, in turn, yields the bound on the
relative difference between infinite value functions in (36), thus concluding
the proof. ∎

Thanks to our assumptions about the features of the running cost and the state
space $X$, Theorem 1 implies asymptotic stability of the compact set $A$ with
respect to which $\ell(\cdot)$ is positive definite and proper, starting from
the initial state $x_{N_{B}}$ and solving the reduced-horizon problem in (8)
(as discussed in [7]). As a direct consequence of this Theorem 1, we can
further formalize our bound (33) on the blocking horizon as follows.

######  Lemma 5 (Further constraining $N_{B}$)

Given $N,N_{B}\in\mathbb{N}$, let Assumption 1 be satisfied by
$V_{N\\!-\\!N_{B}}(x_{N_{B}})$. Then, for the performance bounds (35)-(36) to
hold, the blocking horizon should be upper-bounded as follows:

$N_{B}<N-2\frac{\log(\gamma+1)}{\log(\gamma+1)-\log(\gamma)}.\vspace{.2cm}$
(38)

###### Proof:

The proof is easily obtained by simply manipulating the assumption on $\delta$
and the reduced horizon $N-N_{B}$ required by Theorem 1 (see (34)). ∎

By relying on these results for the reduced horizon problem (8), we can now
state the following sub-optimality bound for the MPC-DMB scheme.

######  Theorem 2 (Suboptimality bound for MPC-DMB)

Let $N,N_{B}\in\mathbb{N}$ be the prediction and blocking horizon of the MPC-
DMB problem in (6), respectively, and let $N_{B}$ satisfy (38). Suppose that
there exists a $\gamma>0$ such that (34) and Assumption 1 hold, with the
latter being verified for the reduced horizon $N-N_{B}$. Then, for all $x\in
X$, it holds that

$\frac{V_{\infty}^{\mu_{N}^{\mathrm{DMB}}}\\!\\!\\!\\!(x)\\!-\\!V_{\infty}(x)}{V_{\infty}(x)}\\!\leq\\!\frac{\gamma^{N\\!-\\!N_{B}}}{(\gamma+1)^{N\\!-\\!N_{B}\\!-\\!2}\\!-\\!\gamma^{N\\!-\\!N_{B}}}\\!+\\!\\!\frac{\delta(x)}{V_{\infty}(x)}.\vspace{.1cm}$
(39)

###### Proof:

Thanks to our assumptions, the bound in (35) holds. Therefore,

$V_{\infty}^{\mu_{N_{B}}}(x_{N_{B}})+\delta(x)-\delta(x)\leq\upsilon(\gamma,N-N_{B})V_{\infty}(x_{N_{B}}),$

with
$\upsilon(\gamma,N-N_{B})=\frac{(\gamma+1)^{N\\!-\\!N_{B}\\!-\\!2}}{(\gamma+1)^{N\\!-\\!N_{B}\\!-\\!2}\\!-\\!\gamma^{N-N_{B}}}$
which further implies that

$V_{\infty}^{\mu_{N}^{\mathrm{DMB}}}(x)-\delta(x)\leq\upsilon(\gamma,N-N_{B})V_{\infty}(x_{N_{B}}),$
(40)

according to (20b). Since $\ell(x,u)\in\mathbb{R}^{+}_{0}$ for all $(x,u)\in
X\times U$ and $x_{N_{B}}$ is obtained as the $N_{B}$ steps ahead prediction
of the state of (1) starting from $x$, the following further holds

$\displaystyle V_{\infty}(x_{N_{B}})$
$\displaystyle=\inf_{u\in\mathcal{U}}\sum_{k=N_{B}}^{\infty}\ell(x_{u}(k),u(k))$
$\displaystyle\leq\inf_{u\in\mathcal{U}}\sum_{k=0}^{\infty}\ell(x_{u}(k),u(k))=V_{\infty}(x),$

where $x_{u}(N_{B})=X_{N_{B}}$ while $x_{u}(0)=x$. Hence, from (40) we obtain
that

$V_{\infty}^{\mu_{N}^{\mathrm{DMB}}}(x)-\delta(x)\leq\upsilon(\gamma,N-N_{B})V_{\infty}(x).$
(41)

Straightforward manipulations of this inequality, leads to

$V_{\infty}^{\mu_{N}^{\mathrm{DMB}}}(x)-V_{\infty}(x)\leq\upsilon(\gamma,N-N_{B})+\delta(x),$
(42)

that, by dividing both terms for $V_{\infty}(x)$ leads to (39) and, thus,
concludes the proof. ∎

This Theorem 2 formalizes a rather intuitive result. Indeed, (39) implies that
the level of suboptimality of the controller resulting from the MPC-DMB scheme
is linked to both the properties of control action at time $t+N_{b}$ (that is
optimized) and the value of the cumulative loss associated with those steps
where the input is blocked.

##  IV Numerical example

The scope of this brief numerical section is to show that, at least in a
simple simulation case study, the DMB approach is feasible way to implement
MPC under computational constraints, without leading to detrimental
performance. Let the system to be controlled be the linear time-invariant
plant

$\left\\{\begin{array}[]{rcl}x_{1}(t+1)&=&0.9\cdot x_{1}(t)+0.1\cdot u(t)\\\
x_{2}(t+1)&=&0.6\cdot x_{1}(t)+0.4\cdot x_{2}(t)\\\
y(t)&=&x_{2}(t),\end{array}\right.$ (43)

while the control problem to solve is as indicated in (2), where the quadratic
cost

$\ell(x_{u}(t),u(t))=x(t)^{T}Qx(t)+u(t)^{T}Ru(t),$

is employed, $R=1$, and

$Q=\begin{bmatrix}10&0\\\ 0&100\end{bmatrix}.$

Let the control variable $u$ be constrained so that it must lie within the
[-1,1] range. No constraints on the state $x$ are instead imposed, namely
$X\equiv\mathbb{R}^{2}$. We consider the achieved closed-loop properties in
terms of regulation to zero, starting from an initial condition with $x_{1}=1$
and $x_{2}=-1$. The prediction horizon $N$ is set equal to $6$ steps.

Let us assume that the computation time needed for the solution of the control
problem amounts to three steps, thus $N_{B}=N_{OP}=3$. For a fair comparison,
we assess DMB MPC as compared to a traditional MPC control with reduced
horizon $N-N_{B}=3$.

The time histories of the state trajectories are illustrated in Figure 1,
where it can be clearly seen that the error between the (unfeasible) receding-
horizon solution and that of the (feasible) DMB-MPC approach remains limited,
even if MPC is obviously faster at the end of the transient, due to a more
rapid update of the control action. This is also confirmed by the small
difference between the Root Mean Square Errors (RMSE) reported in Table I for
the two strategies.

Figure 1: State trajectories in the example of Section IV for a traditional MPC strategy and the proposed DMB rationale. MPC | 15.9134
---|---
DMB MPC | 16.1334
TABLE I: RMSE for MPC and DMB MPC.

##  V Conclusions and future work

In this paper, we introduced a dynamic move blocking MPC approach to comply
with the case in which the optimisation cannot be solved in one sampling
period, without changing the formulation of the problem at hand. The issue
occurs whenever computational resource limitations are relevant, but the
sampling frequency cannot be decreased due to physical constraints.

The proposed strategy is simple yet effective, as it exploits the optimal
solution already available (the one computed at the previous round) while the
new optimization is ongoing. The resulting problem is shown to be equivalent
to a reduced horizon MPC. Asymptotic stability as well as error bounds have
been proven under some mild assumptions.

Future work will be directed toward an extensive experimental assessment of
the proposed strategy, as well as a comparison with the other empirical
alternatives to deal with the problem of limited resources.

## References

  * [1] Alessandro Alessio and Alberto Bemporad.  A survey on explicit model predictive control.  Nonlinear Model Predictive Control: Towards New Challenging Applications, pages 345–369, 2009.
  * [2] F. Allgöwer and A. Zheng.  Nonlinear model predictive control.  Birkhäuser, 2012.
  * [3] Julian Berberich, Johannes Köhler, Matthias A Müller, and Frank Allgöwer.  Linear tracking MPC for nonlinear systems—part I: The model-based case.  IEEE Transactions on Automatic Control, 67(9):4390–4405, 2022.
  * [4] Francesco Borrelli, Alberto Bemporad, and Manfred Morari.  Predictive control for linear and hybrid systems.  Cambridge University Press, 2017.
  * [5] Raphael Cagienard, Pascal Grieder, Eric C Kerrigan, and Manfred Morari.  Move blocking strategies in receding horizon control.  Journal of Process Control, 17(6):563–570, 2007.
  * [6] Ravi Gondhalekar and Jun-ichi Imura.  Recursive feasibility guarantees in move-blocking mpc.  In 2007 46th IEEE Conference on Decision and Control, pages 1374–1379. IEEE, 2007.
  * [7] L. Grune and A. Rantzer.  On the infinite horizon performance of receding horizon controllers.  IEEE Transactions on Automatic Control, 53(9):2100–2111, 2008.
  * [8] Alberto Leva, Simone Formentin, and Silvano Seva.  Overlapping-horizon MPC: A novel approach to computational constraints in real-time predictive control.  In Third Workshop on Next Generation Real-Time Embedded Systems (NG-RES 2022). Schloss Dagstuhl-Leibniz-Zentrum für Informatik, 2022.
  * [9] Tim Schwickart, Holger Voos, Mohamed Darouach, and Souad Bezzaoucha.  A flexible move blocking strategy to speed up model-predictive control while retaining a high tracking performance.  In 2016 European Control Conference (ECC), pages 764–769. IEEE, 2016.
  * [10] Rohan C Shekhar and Chris Manzie.  Optimal move blocking strategies for model predictive control.  Automatica, 61:27–34, 2015.


# DQN-based Beamforming for Uplink mmWave Cellular-Connected UAVs

Praneeth Susarla1, Bikshapathi Gouda1, Yansha Deng2, Markku Juntti1, Olli
Sílven1, Antti Tölli1 email<EMAIL_ADDRESS><EMAIL_ADDRESS>1University of Oulu, Finland  2King’s College
London, United Kingdom

###### Abstract

Unmanned aerial vehicles (UAVs) are the emerging vital components of
millimeter wave (mmWave) wireless systems. Accurate beam alignment is
essential for efficient beam based mmWave communications of UAVs with base
stations (BSs). Conventional beam sweeping approaches often have large
overhead due to the high mobility and autonomous operation of UAVs. Learning-
based approaches greatly reduce the overhead by leveraging UAV data, like
position to identify optimal beam directions. In this paper, we propose a
reinforcement learning (RL)-based framework for UAV-BS beam alignment using
deep Q-Network (DQN) in a mmWave setting. We consider uplink communications
where the UAV hovers around 5G new radio (NR) BS coverage area, with varying
channel conditions. The proposed learning framework uses the location
information to maximize data rate through the optimal beam-pairs efficiently,
upon every communication request from UAV inside the multi-location
environment. We compare our proposed framework against Multi-Armed Bandit
(MAB) learning-based approach and the traditional exhaustive approach,
respectively and also analyse the training performance of DQN-based beam
alignment over different coverage area requirements and channel conditions.
Our results show that the proposed DQN-based beam alignment converge faster
and generic for different environmental conditions. The framework can also
learn optimal beam alignment comparable to the exhaustive approach in an
online manner under real-time conditions.

###### Index Terms:

5G, mmWave, Beam alignment, Deep Q-Network

##  I Introduction

Unmanned aerial vehicles (UAVs) are envisioned as the vital ingredients for
future wireless systems using millimeter wave (mmWave). Especially, the
deployment of cellular-enabled UAV- user equipment (UE)s (hereafter addressed
as UAVs) adds unique features pertaining to high mobility and autonomous
operations to a myriad of civil applications, such as traffic surveillance,
mineral exploration, internet drone delivery systems, etc.[1]. The mmWave
frequencies with multiple input multiple output (MIMO) beamforming and line-
of-sight (LoS) dominant connectivity enable high-speed data access for
unmanned aerial vehicle (UAV)s and also possess challenges such as reliability
and low-latency communication besides interference during aerial-ground
communications. Solving these challenges is also essential for efficient
control of UAVs in fifth generation (5G) and beyond communications. It is
noted that more flexible three-dimensional (3D) beamforming will be deployed
in the forthcoming 5G systems, to enhance both beamforming gain and
interference mitigation by exploiting the angle resolutions in both azimuth
and elevation dimensions of UAVs in the sky[2]. UAVs for aerial-terrestrial
interference management has been extensively studied in base station (BS)-UAV
communication scenario over recent years[3, 4]. Hence, mmWave beam alignment
is critical in supporting the flying UAVs with high mobility and semi-
autonomous operations.

Fast mmWave beam alignment can enhance the data throughput for both UAV-UAV
and BS-UAV communications under 5G and beyond wireless systems. Especially,
the availability of UAV position information at lower frequencies (following
the 3rd generation partnership project (3GPP) standard [5]) may also provide
scope for reliable communication in addition to increasing throughput.
Position information for fast beam alignment has been recently studied under
vehicular context in mmWave systems [6, 7]. On the other hand, high mobility
and autonomous operation of UAVs will require frequent beam realignment as
well. Therefore, a faster and reliable beam alignment using UAV position
information is crucial in enabling high data rate for mmWave UAVs.

Existing works [8, 9] proposed beam tracking schemes using Kalman filters with
high processing complexity. However, these schemes are vulnerable to abrupt
changes in environment when UAVs are moving at high speeds and tracking error
accumulates over time. On the other hand, conventional beam sweeping solutions
[10] often have large overhead, which is also unacceptable with respect to
high speed and autonomous movement of UAVs as it requires frequent beam re-
alignment. An alternate approach to this is to perform fast beam alignment for
every change in UAV position along the BS coverage area. Existing works in
vehicular environment proposed different beam alignment approaches for
terrestrial systems based on black box function optimization[11, 12] and the
use of contextual information [7, 13].

Contextual information generally involves data from the sensors such as
position information, antenna configurations, received signal strength (RSS)
power etc. using low frequency control signals from 3GPP[14] whenever needed,
as this information is used abundantly to reduce the beam training overhead.
The authors in [13] assumed position information as context and perform beam
training over selected subset of beam-pairs by updating a database of channel
strength information in an online manner. The proposed method suits well to
the vehicular communication context with rapid channel variations in the
environment. BS-UAV environment generally involves large coverage areas and
high speed mobility of UAVs resulting in frequent change in UAV position. As a
result, pre-processing of beam pairs and maintaining a database could be
complex, non-generic and also cause significant overhead. In addition,
learning a subset of beam-pairs based on the past beam measurements might not
be essential as there are relatively less channel variations under this
environment.

Black box optimization frameworks in [15, 16], investigate application of
supervised machine learning (ML) techniques to beam alignment problem. These
frameworks are computationally efficient but assumed a separate training data
collection phase for proposed supervised learning environments.

In our work, we propose a generic deep Q-network (DQN)-based framework using
both black box and user-context optimization techniques, in an attempt to
progress from the previous works for UAV beam alignment problem. We
approximate the beam alignment optimization problem for a multi-location
environment using neural networks (NN) techniques and learn best beam-pairs
for any UAV in these locations from their past beam measurements. This
approximation also helps in using prior knowledge of the propagation
environment to significantly reduce the beam-alignment overhead. Our algorithm
takes into account UAV location and finds an optimal beam pair to maximize the
beamforming gain between BS and UAV. To identify the significant change in UAV
location, we consider a grid environment with UAV position information of the
grid element as the user-context information. The main contribution of this
paper is that we propose the reinforcement learning (RL)-based DQN framework
to optimize beam alignment between BS and UAV for any grid position inside the
BS coverage area, during uplink mmWave communication. We benchmark our results
against multi-armed bandit (MAB), exhaustive based approaches, and analyse
against different coverage area requirements under ideal urban macro-cellular
(UMa)- non-line of sight (nLoS) conditions. We also simulate the better
training performance of DQN framework generic to any UAV position in an online
manner under varying channel conditions.

The rest of the paper is organized as follows. Section II presents the problem
formulation and communication modelling, considered in this problem. Section
III discuss in detail about the proposed DQN based RL approach for beam
alignment. Section IV presents the comparison of the proposed DQN-based RL
approach against MAB-based learning and different coverage area requirements
under ideal channel conditions. The section also discuss about the online
learning of DQN and its analysis under varying channel conditions. Section VI
summarizes the conclusion and future work.

##  II System and Communication Model

In this section, we describe the system model, communication model for the
learning framework using 3GPP protocol standards and also formulate the
parameters to be used in the proposed learning framework. The objective of
this problem is to maximize the beamforming between BS and the UAV to provide
efficient communication under the defined environment and channel conditions.

Figure 1: System Model

###  II-A System Model

We consider a cellular mmWave MIMO uplink communications between a single UAV
and a single BS under UMa environments. The BS is fixed at $\mathcal{O}(0,0)$
and communicates with the moving UAV (hereafter used as UE) using a multi-path
mmWave beamforming as shown in Figure 1. The multi-antenna UE hovers randomly
and communicates with the multi-antenna BS in the urban environment following
5G new radio (NR) standard protocol[14]. We consider an analog beamforming
equipped with one radio frequency (RF) chain and uniform linear array (ULA)
structures of $N_{t}$ and $N_{r}$ antennas for both BS and UE, respectively.
The UE transmits (TX) while the BS receives (RX) a radio signal in multiple
beam directions following $\mathcal{B_{\text{TX}}}$ and
$\mathcal{B_{\text{RX}}}$ codebook, respectively with angles defined as

$b_{i}=(i-1)\frac{\pi}{N},i\in\\{1,2,...,N\\},$ (1)

where $b_{i}$ represents a RF radio beam direction with a fixed narrow beam
width ($\frac{\pi}{N}$), $N$ represents $N_{t}$, $N_{r}$ antennas for
$\mathcal{B_{\text{TX}}}$ and $\mathcal{B_{\text{RX}}}$ codebook,
respectively. The codebook values are defined using the beamforming vectors
$\mathbf{w}_{\text{TX}}$ and $\mathbf{w}_{\text{RX}}$ for UE and BS,
respectively, given by
$\mathbf{w}(b_{i})]_{n=0}^{N-1}=\frac{1}{\sqrt{N}}\exp(j\frac{2\pi
nd}{\lambda}\sin(b_{i}))$, where $b_{i}\in\mathcal{B}$, $N=N_{t}$,
$\mathcal{B}=\mathcal{B_{\text{TX}}}$ and $N=N_{r}$,
$\mathcal{B}=\mathcal{B_{\text{RX}}}$ for $\mathbf{w}=\mathbf{w}_{\text{TX}}$
and $\mathbf{w}=\mathbf{w}_{\text{RX}}$, respectively. Here,
$d=\frac{\lambda}{2}$ is the antenna spacing, $\lambda$ is the wavelength and
$b_{i}$ is the $i^{th}$ codebook direction (1). The UE moves randomly along
the 3D coverage area which is divided into multiple grids of equal size, the
set enclosing them is denoted as $\mathbb{U}$. The received signal measurement
can be observed at the BS for different TX-RX beam direction pairs and their
timing information can be estimated using subframe structure (SS) block
symbols, during beam selection procedure of the 3GPP beam-alignment
protocol[14]. For simplicity, we define the time consumed by each such symbol
as the travel time unit (TTU), to measure the communication overhead of the
learning procedure.

###  II-B Communication Model

For the communication model, we consider a multi-path link (LoS or nLoS) radio
channel between UE at time $t$ and BS location $\mathcal{O}\in\mathbb{R}^{3}$.
$\mathrm{UE}_{l}(t)$ represents the UE position on grid index
$l\in\\{0,1,....|\mathbb{U}|\\}$ in the BS coverage area $\mathbb{U}$ at any
time instant $t$, given by

$\mathrm{UE}_{l}(t)=(x_{t},y_{t},z_{t}),$ (2)

where $\\{x_{t},y_{t},z_{t}\\}\in\mathbb{U}$. We assume $\mathrm{UE}_{l}$
(with respect to BS locations) is known during each $P_{1}$ procedure of 3GPP
beam access protocol. Let $\theta_{\mathrm{tx,m}}$, $\theta_{\mathrm{rx,m}}$
be the angle of departure (AoD) and angle of arrival (AoA) of
$m^{\mathrm{th}}$ communication link between BS and UE, respectively. The UE
transmits radio signals in a codebook direction from $\mathcal{B}_{\text{TX}}$
while the BS receives the signal through one of its multiple beam directions
following $\mathcal{B}_{\text{RX}}$ (1). Baseband equivalent of the received
signal is given by

$\begin{multlined}y[k]=\underbrace{\sum_{m=0}^{M}\sqrt{P_{tx}}\beta_{m}\
\mathbf{w}_{\text{RX}}^{H}\mathbf{a}_{R}(\theta_{\mathrm{rx,m}})\mathbf{a}_{T}^{H}(\theta_{\mathrm{tx,m}})\mathbf{w}_{\text{TX}}x[k]}_{\text{$r[k]$}}+\nu[k],\end{multlined}y[k]=\underbrace{\sum_{m=0}^{M}\sqrt{P_{tx}}\beta_{m}\
\mathbf{w}_{\text{RX}}^{H}\mathbf{a}_{R}(\theta_{\mathrm{rx,m}})\mathbf{a}_{T}^{H}(\theta_{\mathrm{tx,m}})\mathbf{w}_{\text{TX}}x[k]}_{\text{$r[k]$}}+\nu[k],$
(3)

where $P_{tx}$ is transmission power, $M$ is the number of multi-paths or
reflection points in the UMa environment [17], $\beta_{m}$,
$\mathbf{a}_{R}({\theta}_{\mathrm{rx,m}})\in\mathbb{C}^{N_{r}}$,
$\mathbf{a}_{T}({\theta}_{\mathrm{tx,m}})\in\mathbb{C}^{N_{t}}$ are the
antenna channel gain and array response vectors for $\theta_{\mathrm{rx,m}}$
and $\theta_{\mathrm{tx,m}}$ along the $m^{th}$ communication link,
respectively. Here,
$\mathbf{a}(\theta)]_{l=0}^{N-1}=\frac{1}{\sqrt{N}}\exp(j\frac{2\pi
ld}{\lambda}\sin(\theta))$, where $\theta=\theta_{\mathrm{rx,m}},N=N_{r}$ and
$\theta=\theta_{\mathrm{tx,m}},N=N_{t}$ for
$\mathbf{a}_{R}({\theta}_{\mathrm{rx,m}})$ and
$\mathbf{a}_{T}({\theta}_{\mathrm{tx,m}})$, respectively.
$\mathbf{w}_{\text{RX}}\in\mathbb{C}^{N_{r}}$,
$\mathbf{w}_{\text{TX}}\in\mathbb{C}^{N_{t}}$ are the transmit and receive
unit-norm beamforming vectors, $\nu[k]\sim\mathcal{CN}(0,WN_{0})$ is the
effective noise with zero mean and two-sided power spectral density
$\frac{N_{0}}{2}$, $x[k]$ represents one orthogonal frequency division
multiple access (OFDM) symbol of the time-domain transmitted signal with
bandwidth $W$ and TTU time period with
$\frac{1}{K}\sum_{k=0}^{\text{K}}\|x[k]\|^{2}=1$. Here, $k=0,1,...K$ is the
number of samples spanned over TTU time. In this work, we assume the channel
measurements $\mathbf{H}_{m}$ of the $m^{th}$ multi-path link
$\mathbf{H}_{m}(\theta_{\text{tx},m},\theta_{\text{rx},m})\triangleq\beta_{m}\mathbf{a}_{R}({\theta}_{\mathrm{rx,m}})\mathbf{a}_{T}^{H}({\theta}_{\mathrm{tx,m}})$
to follow 3GPP UMa channel conditions[17]. We define
$r[k]=\sum_{m=0}^{M}\sqrt{P_{tx}}\mathbf{w}_{\text{RX}}^{H}\mathbf{H}_{m}(\theta_{\text{tx},m},\theta_{\text{rx},m})\mathbf{w}_{\text{TX}}x[k]$.
Then, the signal-to-noise ratio (SNR) is given as
$\mathrm{SNR}=\frac{\frac{1}{K}\sum_{k=0}^{K}\mathinner{\\!\left\lVert
r[k]\right\rVert}^{2}}{N_{0}W}$ and overall rate measurement $R$ in bits per
channel use is given by

$R=\log(1+\mathrm{SNR}).$ (4)

Thus, the beam alignment for the above formulation is realized through data
rate measurements. The optimal beam-pair for UE-BS is selected based on their
data rates under the scenario mentioned in Section II-A.

###  II-C Problem Formulation

We consider an uplink communication between BS and UE following 3GPP beam
access protocol [18, 14]. The beamforming protocol in general involves mainly
three procedures, Initial communication (used as $P_{1}$ procedure), beam
selection (used as $P_{2}$ procedure) and beam refinement (used as $P_{3}$
procedure) [18]. Here, we assume UE and BS share position information and
sequence of beam-pairs with one another at low frequencies during $P_{1}$,
respectively and sweep the selected sequence of beam pairs during $P_{2}$
following 5G communication [14]. We formulate the 3GPP based beam-pair
alignment learning through $P_{1}$ and $P_{2}$ procedures, and maximize the
beamforming gain for any UE position around the BS coverage area $\mathbb{U}$.
Timing information of these procedures assumes 3GPP frame structure and time
for each OFDM slot is defined in our work as the TTU[18]. We consider RSS of
radio signal and radio beam pair directions (both TX and RX) as the known and
unknown parameters of this multi-location environment, respectively. We define
the state and action spaces for learning based beam-pair alignment as follows:

$\displaystyle(\mathcal{E})\mathrel{\mathop{\mathchar 58\relax}}\ \mathcal{S}$
$\displaystyle=\\{\mathrm{UE}_{l},1\leq l\leq|\mathbb{U}|\\}$ (5)
$\displaystyle\mathcal{A}$ $\displaystyle=\\{(b_{p},b_{q}),1\leq
p\leq|\mathcal{B}_{\mathrm{TX}}|,1\leq q\leq|\mathcal{B}_{\mathrm{RX}}|\\},$

where $\mathrm{UE}_{l}$ is the location of UE within coverage area
$\mathbb{U}$ (2) while $\mathcal{B}_{\mathrm{TX}}$,
$\mathcal{B}_{\mathrm{RX}}$ is the beam codebook sets at UE and BS side
respectively (1).

RL is an interactive learning problem following a markov decision process
(MDP). In this work, the RL-based beam alignment problem is modelled as a
partially observable Markov decision process (POMDP). At any time instant $t$,
we define the parameters $s_{t}=\\{(s^{\prime},a^{\prime})\
s^{\prime}\in\mathcal{S},a^{\prime}\in\mathcal{A}\\}$, $a_{t}\in\mathcal{A}$
and $r_{t}\in\mathcal{R}$ where $s_{t}$, $a_{t}$, $r_{t}$ are the state,
action and reward at time instant $t$. $a^{\prime}$ corresponds to the set of
previous actions applied for state transitions until the time instant $t$.
Data rate measurements computed on applying each action are considered as the
rewards for the problem. We denote
$o_{t}=\\{a_{t-1},s_{t-1},a_{t-2},s_{t-2},....,a_{1},s_{1}\\}$ as the observed
history of all such state information and past actions. After the 3GPP initial
communication procedure with UE, BS starts with a random receiving beam
direction and then proceeds towards the maximum beamforming gain by applying
actions and undergoing state transitions, accordingly. The current applied
action becomes part of the next state, undergoing state transition. The
objective of this problem can be formulated as

$\displaystyle(P1)$ $\displaystyle\mathrel{\mathop{\mathchar
58\relax}}\max_{\\{\pi(a_{t}|o_{t})\\}}\sum_{i=t}^{\infty}\gamma^{i-t}\mathbb{E}_{\pi}[r_{a_{i}}(i)],$
(6) $\displaystyle s.t.$ $\displaystyle r_{a_{t}}(t)$
$\displaystyle=\begin{cases}1&\text{if }R(a_{t})\geq R_{max}(o_{t})\\\
-1&\text{otherwise}\end{cases},$ $\displaystyle\gamma\in(0,1],$

where $R_{max}(o_{t})$ is the optimal data rate measurement observed among the
information history $o_{t}$ until time instant $t$, $r_{a_{t}}(t)$ and
$R(a_{t})$ are the rewards and data rate measurement observed on applying
action beam-pair $a_{t}$, respectively. We follow DQN approach to solve this
RL objective problem.

##  III Implementation

DQN is a value-based approach used generally in the context of RL [19]. The
approach learns an optimal approximated policy of states mapping to actions
$\pi(s)=a$ by parameterizing and estimating state-action value function
$Q(s,a;\theta)$ using deep neural networks (DNN). We denote the primary DNN
network weight matrix and target DNN network weight matrix as $\theta$ and
$\overline{\theta}$, respectively [19]. We consider a fully connected DNN for
both the networks where $\overline{\theta}$ is updated with primary network
parameters $\theta$, after every $K$ iterations. The input of DNN is given by
the variables in $s_{t}$. The intermediate layers are fully connected linear
units with rectifier linear units (ReLU) by using the function
$f(x)=\max{(0,x)}$ and the output layer is composed of linear units, which are
in one-one corresponding relationship with the action space $\mathcal{A}$. We
consider initialization of bias and weights of these layers with zeros and
Kaiming normalization, respectively.

Input: The set of UAV x,y,z location coordinates and training episodes M

1 Algorithm hyper-parameters: learning rate $\xi\in(0,1]$, discount rate
$\gamma\in[0,1)$, $\epsilon$-greedy rate $\epsilon\in(0,1]$, update steps $K$;

2 Initialization of replay memory $M$ to capacity $C$, the primary Q-network
with parameters $\theta_{1}$, the target Q-network with parameters
$\theta_{2}$

3 $\mathcal{S},\mathcal{A}$: State and Action space of DQN agent

for  _episode $\leftarrow 1\ to\ M$ // for each episode _

4 do

5 Initialization of $s_{1}$ by executing a random action $a_{0}$

6 n=0, $N\rightarrow$Episode limit

7 while  _$True$_ do

8 if  _$p_{\epsilon} <\epsilon$_ then

9 select a random action $a_{t}\in A$

10 else

11 select $a_{t}={argmax}_{a\in\mathcal{A}}Q(s_{t},a,\theta)$

12 BS applies $a_{t}$ over the channel, receive signal for $(t+1)^{th}$
episode during uplink communication

13 UE observes $S^{t+1}$ and calculate the reward

14 Store transition $e=(s_{t},a_{t},r_{t+1},s_{t+1})$ in replay memory $D$

15 Sample random minibatch of transitions $U(D)$

16 Compute Loss and Perform gradient descent for $Q(s,a;\theta)$

17 Update the target network parameters $\theta_{2}=\theta_{1}$ after every
$K$ steps

$n=n+1$ // Increment episode time

18 if  _done or $(n=N)$_ then

break // End episode

19

20

21

Algorithm 1 DQN approach

At a time instant $t$, DQN selects action $a_{t}\in\mathcal{A}$ and perform
forward propagation of $Q(s_{t},a_{t};\theta)$ following $\epsilon$-greedy
policy. A memory buffer of experiences
$D_{t}=\\{e_{1},e_{2},e_{3},...,e_{t}\\}$,
$e_{i}=(s_{i},a_{i},r_{i+1},s_{i+1})$ are collected, where a mini batch of
them $U(D)$ are randomly sampled and sent into DQN [19]. During back
propagation, a mean squared error (MSE) loss function is computed between
primary, target networks and $\theta$ is updated using stocastic gradient
descent (SGD) and Adam Optimizer as

$\theta_{t+1}=\theta_{t}-\xi_{\mathrm{Adam}}\nabla\mathcal{L}^{\mathrm{DQN}}(\theta_{t}),$
(7)

where $\xi_{\mathrm{Adam}}$ is the learning rate,
$\nabla\mathcal{L}(\theta_{t})$ is the gradient of the DQN loss function.
Complete steps followed by DQN for RL based beam alignment problem are shown
in Algorithm 1. Here, we define episode as the consecutive set of actions
applied on the starting state until it reaches the terminal state with maximum
beamforming gain for that location. In order to prevent episodes with infinite
set of actions during training, we confine maximum episode length to
exhaustive set of beam pairs possible under the chosen antenna configuration.

As the RL learning objective formulation involves both current data rate
$R_{t}$ and best observed data rate $R_{\mathrm{max}}(o_{t})$ measurements
(shown in (6)), we consider the overall online training procedure of DQN
framework under two phases namely, Warmup phase and Training phase as shown in
Fig. 2. During the Warmup phase, the exploration is set to maximum, in order
to observe the best possible data rate for the given UE location by applying
maximum episode length of actions. During the Training phase, the algorithm
continues to reduce its exploration and move towards exploitation following
$\epsilon$-greedy policy. The episode starts with initial action corresponding
to $o_{t}$ and applies next actions to reach the terminal state as quickly as
possible. The Warmup phase results in extra training time at the start but
favours quick convergence of DQN during training phase resulting in faster
beam-alignment training for the multi-location environment.

Figure 2: DQN training procedure consisting of two phases namely, Warmup phase
and Training phase

##  IV Simulation results

As described in Section III, we implement the RL-based beam alignment using
DQN, following $(P1)$ objective (6) and steps mentioned in Algorithm 1.
Similarly, we implement the MAB-based approach for our problem (following
existing works [13, 20]) using greedy upper confidence bound (gUCB) algorithm.
Traditional approach mainly involves exhaustive search over entire action
space $\mathcal{A}$, to find the best beam pair with maximum beamforming gain
between UE and BS. In our work, the exhaustive approach in general causes
significant communication overhead ($\mathcal{O}(|\mathcal{A}|)$) with higher
antenna elements as frequent beam scanning is required for every change in
grid element unit of UE inside $\mathbb{U}$. On the other hand, RL and MAB
learning-based methods once converged, can significantly reduce the
communication overhead during $P_{2}$ procedure and maximize beamforming gain
in $\mathcal{O}(1)$ time.

In this section, we first investigate the training performance of our proposed
RL-based approach against MAB-based learning and different coverage area
requirements in order to maximize beam alignment under ideal UMa-nLoS channel
conditions. Later, we combine these observations and perform online beam
alignment in the presence of channel variation conditions.

Figure 3: DQN reward-time plots for various BS coverage area requirements for
ideal UMa-nLoS channel conditions

###  IV-A DQN training performance

As shown in Fig. 3, brown and blue plots represent the accumulated rewards
over time for DQN and gUCB, respectively, under single UE grid location
($\mathbb{U}=1$) and ideal UMa-nLoS channel conditions with 6 reflection
points. Our results show that the DQN-based approach quickly accumulate the
rewards over TTU time for learning the optimal beam pair. The gUCB approach
involves iterating over exhaustive action space $\mathcal{A}$, thus consuming
more time at every episode of its training procedure. On the other hand, DQN
uses its warmup phase to quickly determine the best possible data rate
measurement for the current state and then learn the determined best rate
during training phase, by applying actions in a iterative manner. With
increase in BS coverage area, the gUCB-based approach is observed to be
unreliable for multiple grid locations without frequent re-training. This is
due to substantial change in RSS measurements w.r.t (TX,RX) beam pairs for
every UE grid location.

On the other hand, DQN agent with the same architecture is very much reliable
to learn beam alignment across multiple UE locations for different coverage
area requirements under the same ideal channel conditions as shown in Figure
3. After the convergence is obtained, BS can quickly align using learnt
optimal beams for any UE position within the coverage area without any re-
training. Also, the learning is observed to be relatively quicker in
convergence with increase in coverage area of BS. With increase in coverage
area, neighbouring grid elements with similar optimal beam pairs converge DQN
faster as part of its MDP process, resulting in average less number of
training iterations per grid element alongside convergence. Thus, the DQN
agent can achieve faster and reliable beam alignment due to its training
procedure and the neighbourhood grid element convergence.

###  IV-B online DQN performance with channel variation and shadow fading

In this subsection, we plot DQN training performance in real-time conditions
in an online manner by considering change in channel conditions, thermal
noise, slow fading and slow channel variation at UE grid locations as shown in
Fig. 4.

Table I: Simulation Parameters cyan Parameters | Value
---|---
mmWave Channel | UMa
mmWave freq |  $30$ GHz
carrier spacing freq $df$ |  $60$ kHz
antenna element spacing $d$ | $0.5$
Transmit power $P_{\mathrm{tx}}$ |  $0$ dBm
Transmit antenna elements $N_{tx}$ | 8
Receiving antenna elements $N_{rx}$ | 8
Noise Level $N_{0}$ | -174 dBm
BS location | $[0,0,25]$
coverage xloc $\mathbb{U}_{\mathrm{xloc}}$ |  $[-60,60,20]$ m
coverage yloc $\mathbb{U}_{\mathrm{yloc}}$ |  $[-60,60,20]$ m
coverage zloc $\mathbb{U}_{\mathrm{zloc}}$ |  $[41.5,81.5]$ m
Cardinality of State Space $|\mathcal{S}|$ | $72$
Cardinality of Action Space $|\mathcal{A}|$ | $64$

Fig. 4(a) and Fig. 4(b) plot the rewards and RSS error of the DQN learning
agent, respectively. We define RSS error as the expectation of the difference
between DQN-agent and exhaustive approach RSS errors over BS coverage area.
For these simulations, we consider a ($\mathbb{U}=32$) coverage area
environment between BS and UE. The environment is equipped with thermal noise,
shadow fading, UMa-nLoS to UMa-LoS channel conditions along with a slow rician
channel variation. The parameters used for this simulation are disclosed under
Table. I. Depending on the defined coverage area, different path loss models
including the aerial view of UE for both UMa-nLoS and UMa-LoS are included
following 3GPP[17]. It is noted here that same set of hyper-parameters are
used throughout the DQN simulation. We observe that DQN agent performs similar
to previous results, converging well under both varying channel conditions.
The disturbance in the smoothness of the reward plots observed in Fig. 4(a)
could be due to the impact of channel variation under nLoS conditions. From
Fig. 4(b), we observe that the RSS error increases at first during online
training phase and then decreases to zero on model convergence under both
varying LoS and nLoS conditions. This shows that DQN-framework on convergence,
achieves optimal data rate measurements similar to that of exhaustive
approach.

(a) Reward plot for DQN based beamlearning with channel variation including
shadow fading, thermal noise and rician fading

(b) RSS Error plot for DQN based beam learning with channel variation
including shadow fading, thermal noise and rician fading

Figure 4: DQN-based beam learning performance under UMa-nLoS to UMa-LoS
varying channel conditions

##  V Acknowledgements

The research was supported by 6G Flagship programme, Finland and the
Engineering and Physical Research Council (EPSRC), U.K., under Grant
EP/W004348/1.

##  VI Conclusion and Future Work

In this paper, we developed a learning-based beam alignment framework for
mmWave MIMO uplink BS-UAV communication. We proposed a RL-based framework
using DQN to maximize the beam alignment for any UAV position within the BS
coverage area following 3GPP standard. We also analyze the same DQN
architecture over different coverage requirements under ideal channel
conditions to demonstrate the generalization of the proposed RL-based
framework for beam learning. Our results show that, the proposed approach
significantly outperforms the training performance of the MAB-based method and
also learns optimal beam pairs comparable to that of the heuristic method in
an online manner under varying channel conditions. Having shown some promising
results, we will address the full capabilities of these generic learning
architectures under higher MIMO antenna configurations, large number of beam-
directional pairs, interference mitigation etc. in future.

## References

  * [1] L. Zhang, H. Zhao, S. Hou, Z. Zhao, H. Xu, X. Wu, Q. Wu, and R. Zhang, “A Survey on 5G Millimeter Wave Communications for UAV-Assisted Wireless Networks,” _IEEE Access_ , vol. 7, pp. 117 460–117 504, July 2019.
  * [2] Y. Zeng, J. Lyu, and R. Zhang, “Cellular-connected uav: Potential, challenges, and promising technologies,” _IEEE Wireless Communications_ , vol. 26, no. 1, pp. 120–127, September 2018.
  * [3] X. Lin, V. Yajnanarayana, S. D. Muruganathan, S. Gao, H. Asplund, H.-L. Maattanen, M. Bergstrom, S. Euler, and Y.-P. E. Wang, “The sky is not the limit: Lte for unmanned aerial vehicles,” _IEEE Communications Magazine_ , vol. 56, no. 4, pp. 204–210, April 2018.
  * [4] M. M. Azari, F. Rosas, A. Chiumento, and S. Pollin, “Coexistence of terrestrial and aerial users in cellular networks,” in _2017 IEEE Globecom Workshops (GC Wkshps)_. IEEE, Dec 2017, pp. 1–6.
  * [5] E. T. S. Institute, “LTE; Evolved Universal Terrestrial Radio Access (E-UTRA); Radio Resource Control (RRC); Protcol Specification, Technical Report 36.331, version 15.10.0, Release 15,” _European Telecommunications Standards Institute_ , July, 2020.
  * [6] F. Devoti, I. Filippini, and A. Capone, “Facing the millimeter-wave cell discovery challenge in 5g networks with context-awareness,” _IEEE Access_ , vol. 4, pp. 8019–8034, Nov 2016.
  * [7] J. C. Aviles and A. Kouki, “Position-aided mm-wave beam training under nlos conditions,” _IEEE Access_ , vol. 4, pp. 8703–8714, Nov 2016.
  * [8] J. Zhao, F. Gao, W. Jia, S. Zhang, S. Jin, and H. Lin, “Angle domain hybrid precoding and channel tracking for millimeter wave massive mimo systems,” _IEEE Transactions on Wireless Communications_ , vol. 16, no. 10, pp. 6868–6880, Aug 2017.
  * [9] J. Zhao, F. Gao, L. Kuang, Q. Wu, and W. Jia, “Channel tracking with flight control system for uav mmwave mimo communications,” _IEEE Communications Letters_ , vol. 22, no. 6, pp. 1224–1227, April 2018.
  * [10] J. Saloranta, G. Destino, and H. Wymeersch, “Comparison of different beamtraining strategies from a rate-positioning trade-off perspective,” in _2017 European Conference on Networks and Communications (EuCNC)_. IEEE, June 2017, pp. 1–5.
  * [11] T. Kadur, H.-L. Chiang, and G. Fettweis, “Effective beam alignment algorithm for low cost millimeter wave communication,” in _2016 IEEE 84th Vehicular Technology Conference (VTC-Fall)_. IEEE, Sept 2016, pp. 1–5.
  * [12] P. Susarla, Y. Deng, G. Destino, J. Saloranta, T. Mahmoodi, M. Juntti, and O. Sílven, “Learning-based trajectory optimization for 5g mmwave uplink uavs,” in _2020 IEEE International Conference on Communications Workshops (ICC Workshops)_. IEEE, June, 2020, pp. 1–7.
  * [13] V. Va, T. Shimizu, G. Bansal, and R. W. Heath, “Online Learning for Position-Aided Millimeter Wave Beam Training,” _IEEE Access_ , vol. 7, pp. 30 507–30 526, Mar 2019.
  * [14] E. Dahlman, S. Parkvall, and J. Skold, _5G NR: The next generation wireless access technology_. Academic Press, 2020.
  * [15] Y. Wang, N. J. Myers, N. González-Prelcic, and R. W. Heath, “Site-specific online compressive beam codebook learning in mmwave vehicular communication,” _IEEE Transactions on Wireless Communications_ , vol. 20, no. 5, pp. 3122–3136, 2021.
  * [16] C. Jiang, H. Zhang, Y. Ren, Z. Han, K.-C. Chen, and L. Hanzo, “Machine learning paradigms for next-generation wireless networks,” _IEEE Wireless Communications_ , vol. 24, no. 2, pp. 98–105, Dec 2016.
  * [17] E. T. S. Institute, “Study on Enhanced LTE support for Aerial Vehicles 36.777 version 1.0.0 Release 15,” _European Telecommunications Standards Institute_ , Jan, 2018.
  * [18] ——, “Study on 5G New Radio(NR) access Technology, Technical Report 38.912, version 15.0.0, Release 15,” _European Telecommunications Standards Institute_ , May, 2017.
  * [19] V. Mnih, K. Kavukcuoglu, D. Silver, A. A. Rusu, J. Veness, M. G. Bellemare, A. Graves, M. Riedmiller, A. K. Fidjeland, G. Ostrovski _et al._ , “Human-level control through deep reinforcement learning,” _Nature_ , vol. 518, no. 7540, p. 529, Feb, 2015.
  * [20] M. Hashemi, A. Sabharwal, C. E. Koksal, and N. B. Shroff, “Efficient beam alignment in millimeter wave systems using contextual bandits,” in _IEEE INFOCOM 2018-IEEE Conference on Computer Communications_. IEEE, April, 2018, pp. 2393–2401.

  *[UAVs]: Unmanned aerial vehicles
  *[mmWave]: millimeter wave
  *[UE]: user equipment
  *[MIMO]: multiple input multiple output
  *[LoS]: line-of-sight
  *[UAV]: unmanned aerial vehicle
  *[5G]: fifth generation
  *[3D]: three-dimensional
  *[BS]: base station
  *[3GPP]: 3rd generation partnership project
  *[RSS]: received signal strength
  *[ML]: machine learning
  *[DQN]: deep Q-network
  *[NN]: neural networks
  *[RL]: reinforcement learning
  *[MAB]: multi-armed bandit
  *[UMa]: urban macro-cellular
  *[nLoS]: non-line of sight
  *[NR]: new radio
  *[RF]: radio frequency
  *[ULA]: uniform linear array
  *[SS]: subframe structure
  *[TTU]: travel time unit
  *[AoD]: angle of departure
  *[AoA]: angle of arrival
  *[OFDM]: orthogonal frequency division multiple access
  *[SNR]: signal-to-noise ratio
  *[MDP]: markov decision process
  *[POMDP]: partially observable Markov decision process
  *[DNN]: deep neural networks
  *[ReLU]: rectifier linear units
  *[MSE]: mean squared error
  *[SGD]: stocastic gradient descent
  *[gUCB]: greedy upper confidence bound