Enhance. laser-induced excit. prob.]Enhancement of the
laser-induced excitation probability of the hyperfine
ground state of muonic hydrogen by a multipass cavity setup

[1,2]Rakesh Mohan<EMAIL_ADDRESS>
3]Masahiko Iwasaki

*[1]School of Applied Sciences, KIIT University, Patia, Bhubaneswar, 751024, Odisha, India

[2]Center of Excellence in High Energy and
Condensed Matter Physics, Department of Physics, Utkal University, Vani Vihar, Bhubaneswar, 751004, Odisha, India

[3]The Institute of Physical and Chemical Research, RIKEN, 2-1 Hirosawa, Wako, Saitama, 351-0198, , Japan

We study the enhancement of the magnetic dipole induced
excitation probability of the hyperfine ground state of
Doppler-broadened muonic hydrogen ($p \mu^{-}$)
by a nanosecond laser pulse in the mid-infrared range with Gaussian
temporal shape such
that the pulse bandwidth is broader than the Doppler width at
10 K. The enhancement is achieved by shrinking the
cross-section of the laser pulse and placing the muonic hydrogen
medium in a  multipass cavity, while preserving the total
irradiated target volume.
We numerically solve a set of Maxwell-Schr$\ddot{\rm o}$dinger
equations to obtain the excitation probability and the total
efficiency for various densities of the muonic hydrogen atomic
medium and at various positions in the multipass cavity. For the
typical range of densities of muonic hydrogen atoms
at major proton accelerator facilities such as the J-PARC
(density $\sim$ $10^5$cm$^{-3}$),
the laser propagation effect is insignificant. For such cases,
the total efficiency increases by an order of two for 100
reflections with a uniform polarization.
If the density exceeds the
value of $10^{17}$cm$^{-3}$ as might be in the future advances,
the laser propagation effect has to
be taken into account, and the total efficiency decreases with
the number of reflections giving rise to a pulsed polarization of
the beam. Our study can serve as a guideline for the development of
a polarized muonic beam for a precise measurement of
the ground state hyperfine splitting of muonic hydrogen, or  for
$\mu$SR experiments.

Amongst all the exotic atoms, muon ($\mu^{\pm}$) prevails
as the probing tool
as far as the theoretical studies and experiments in atomic and
nuclear physics are concerned. This is due to its relatively
easier availability and longer lifetime of 2.2 $\mu$s.
Charged pion which is often the next in line amongst the exotic atoms,
has a lifetime of 26 ns, making it very short-lived. Besides,
unlike the pion,
the muon is not influenced by the strong force of the atomic nuclei
in the lattice, and hence has higher penetration depth.
Muon has been studied extensively due to its relevance
in the Standard model of elementary particles, and its usefulness in
various applications such as $\mu$SR experiment, and
muon-catalyzed fusion, etc [1, 2, 3, 4, 5, 6, 7, 8, 9, 10].
Commonly, for the application purpose of muon (whether
positively or negatively charged), its kinetic energy is reduced
from as high as a few MeV after being ejected from the particle
accelerator, down to a few keV $\sim$ a few eV.
In this process, positive muons capture electrons to form muonium
atoms while negative muons are captured by protons to form muonic
hydrogen atoms. In both the cases of the recombination processes,
the degree of polarization, which is 100 % upon birth, is
significantly lost which needs to be spin-repolarized.

Several works have been done on obtaining slow positive muon beam
[11, 12, 13, 14, 15, 16, 17, 18].
We turn our focus on the muonic hydrogen which
has proved to be a suitable
candidate to precisely determine the size of proton.
This is because of the proximity of the negative muon to the proton
that amplifies the atomic properties of the
muonic hydrogen due to the presence of the negative muon
(200 times heavier particle than the electron) in the
atomic orbit with a radius of about 1/200 th of
that of hydrogen and muonium.
For comparison, we have used the figure from our previous work
[19] which summarizes the various atomic parameters
of hydrogen, muonium, and muonic hydrogen in Fig. <ref>.

The famous study on the Lamb shift
($2S_{1/2} - 2P_{1/2}$ energy difference) of muonic hydrogen
at PSI [2] reported that the proton radius is
about $4 \%$ smaller than the widely accepted value. This puzzle
has spearheaded many works on spectroscopy and particle physics
[20, 21, 22].
Along this line, Adamczak et al. proposed a
method of determining the proton size with accuracy by an
independent measurement of the hyperfine splitting
($1S (F=0) - 1S (F=1)$) of the ground state
of muonic hydrogen [23, 24, 25]. The $1S (F=0) - 1S (F=1)$ transition
wavelength is in the mid-infrared range (6.76 $\mu$m), is
dipole-forbidden, and the process involves magnetic dipole
transition which is several orders of magnitude weaker
than the electric dipole transition.
The experimental scheme was originally proposed by Bakalov
et al. [26], in which
negative muons are slowed down and stopped in a pure hydrogen target
between two parallel gold or
aluminium plates to form muonic hydrogen. Recently,
we also have studied the laser-induced
hyperfine ground state excitation probability
of Doppler-broadened muonic hydrogen for various pulse parameters
that can serve as a guideline for the development of the relevant
laser source, depending on whether the purpose of
exciting muonic hydrogen is measuring the ground state hyperfine
splitting of muonic hydrogen with precision, or
producing spin-polarized muonic hydrogen
beam [19].
Adamczak et al. have
proposed a scheme of placing the target inside a
multipass cavity in order to
increase the excitation probability which is otherwise merely
$1.2 \times 10^{-5}$ (too low for the proposed experiment
to be feasible) with a laser pulse of long duration at
6.76 $\mu$m with 0.25 mJ pulse energy in a single interaction
of the laser pulse with the muonic hydrogen beam at a temperature
of 300 K [23].
Undoubtedly, this experiment comes with huge technical
From a theoretical point of view, a study of the elementary
challenges such as the deterioration
of the laser field profile as it propagates through the muonic
hydrogen medium inside the multipass cavity for various ranges
of medium densities, is essential for completeness. For
the muon beam facilities around the world at present, the
achievable beam density is not high enough for any laser
propagation effect to arise. For example, at the J-PARC,
the density of the muon beam is typically $10^5$cm$^{-3}$.
Nevertheless, it is important
to study the laser propagation effect on the polarization of the
muons, as the density is increased, as might be possible
by the technological development in future.
Besides, the prospects and limitations in tailoring the
temporal profile of the muonic hydrogen beam in such a
multipass cavity setup, needs to be discussed.
We note that, pulsed muon and CW muon sources that have been
developed at various facilities across the world, have their
own advantages and disadvantages depending upon the purpose
of the muon source. For example, the PSI and
TRIUMF [27, 28] are both home
to intense continuous muon beam
which is suitable for spectroscopy experiments, whereas
KEK, RIKEN-RAL and J-PARC MUSE [29, 30] serve as the sources of pulsed
muon beams which are better suited for $\mu$SR experiments.
It is to be seen in this context,
how the multipass cavity setup might
play its role in shaping the temporal profile of the muonic hydrogen

In the present work, we study the excitation probability of
the muonic hydrogen by nanosecond laser pulses in a multipass
cavity setup and investigate the conditions under which
the excitation probability is enhanced which leads to increased
polarization of muons. The bandwidth of the laser pulse we
consider in our study is broader as compared to the Doppler width
of the atoms. The results we present in this work would be useful
in the development of experimental setup for polarized muonic beam.

(Color online) Comparison of the hydrogen and
hydrogen-like systems involving
positive and negative muons.

§ SYSTEM DESCRIPTION

(Color online)
Schematic diagram showing the muonic hydrogen atoms inside a
multipass cavity with mirrors of perfect reflectivity. An infrared
laser pulse is shrunk in its beam cross-section and goes
through a number of reflections such that the interaction volume of
the muonic hydrogen and the laser pulse is preserved.

Figure <ref> shows the schematic diagram of the
system we have considered.
We consider the hyperfine sublevels of the ground ($1s$)
state of the
muonic hydrogen. A laser beam in the mid-infrared range (6.7 $\mu$m)
is shrunk by reducing its cross section by a factor $k$ (thereby
increasing the intensity by the same factor), and projected into a
multipass cavity with perfect reflectivity. The muonic hydrogen beam
is passed through the multipass cavity providing $k$ reflections
of the
laser beam, thereby preserving the total irradiated target volume.
The assumption of a
multipass cavity with perfect reflectivity is for the sake of
simplicity given that the intensity of the laser beam is barely
attenuated by the mirrors of reflectivity as high as 99.999 % for
the number of reflections considered in this work.
In steady state, the slow beam of muonic hydrogen, after the
fast muons being captured by the hydrogen atoms inside the cavity,
can be treated to be frozen in space when compared to the laser
pulse passing through the multipass cavity.

We assume a Gaussian pulse which
excites the muonic hydrogen from the
hyperfine singlet state, $1s (F=0)$, to the triplet state,
$1s (F=1)$,
through magnetic dipole ($M1$) interaction. We note that, the
hyperfine transition is electric dipole forbidden and it is the
magnetic dipole component which has the most significant contribution
in the laser-matter interaction.
The magnetic field of the pulse with central frequency $\omega$ is
written as

\begin{align}
\label{eqn:field1}
\mathcal{B}(t) & =  B(t) \exp{\left[-i\omega t \right]} + \mathrm{c.c.},
\end{align}

where $B(t)$ is the magnetic field amplitude of the pulse. The
magnetic transition dipole moment of the muonic hydrogen atom,
$M$, calculated by us in our previous work is
$\frac{1}{\sqrt{6}} \frac{e\hbar}{2m_\mu}$ [19].
The Rabi frequency $\Omega$, corresponding to the magnetic dipole transition in S.I. unit, is

\begin{align}
\label{eqn:rabins}
\Omega ({\rm ns}^{-1}) & = 1.77 \times 10^{-6}\sqrt{I(t)}
\end{align}

where $I(t)$ is the laser intensity in units of W/cm$^2$.
We employ a set of Maxwell-Schr$\ddot{\rm o}$dinger equation,
which can be easily obtained by
coupling the probability amplitudes, $a_0$ and $a_1$, of the
hyperfine states, $1s (F=0)$ and $1s (F=1)$, respectively and
Maxwell's equations, which reads

\begin{align}
\label{eqn:amp1}
\frac{\partial}{\partial t} a_0(z,t) = i \Omega^{*}(z,t) a_1(z,t) ,
\end{align}

\begin{align}
\label{eqn:amp2}
\frac{\partial}{\partial t} a_1(z,t) = i \Omega(z,t) a_0(z,t) ,
%	- (i \delta + \gamma)a_1(z,t) ,
\end{align}

\begin{align}
\label{eqn:amp3}
\frac{\partial}{\partial z} \Omega(z,t) =  -2i\alpha a_1^{*}
\end{align}

where $\alpha = N \omega |M|^2 \mu_0 c /2$, is defined as the laser
propagation coefficient of the muonic hydrogen medium, $N$
being the density of the medium [31]. While obtaining
the above equations, we have neglected the spontaneous decay
rate of the muonic hydrogen for the range of time duration considered.
Although the slow muonic hydrogen atoms formed are unpolarized,
i.e. 50$\%$ of the atoms are in the state $1s(F=0)$ and 50$\%$ are in
state $1s(F=1)$, the initial condition that all the atoms are in the
singlet state, $1s(F=0)$, simplifies the study without the
loss of generality. We numerically calculate the Doppler-averaged
spin-flip probability, $P$, of the ground hyperfine state, from
$1S (F=0)$ to $1S (F=1)$, of the muonic hydrogen,
taking into account the laser
propagation effect of the medium on the laser pulse as it
passes through the multipass cavity.

(Color online)
(a) Normalized distribution of the number of Doppler-broadened
muonic hydrogen atoms at a temperature of 300 K (blue line) and
10 K (green line). The area under each of the lines is unity.
(b) Same as panel (a) multiplied with a
Gaussian window function of width 220 MHz, which corresponds to
the spectral width of a 2 ns Gaussian pulse. The ratio of the
areas under the green and red line is 0.917/0.387 = 2.3.

(Color online)
Temporal variation of the population of the ground hyperfine state,
$1s(F = 1)$ and the laser field, $\Omega$ with the medium propagation coefficient,
$\alpha = 4 \times 10^{-18}$ mm$^{-1}$ [(a), (b)],
$4 \times 10^{-6}$ mm$^{-1}$ [(c), (d)],
$4 \times 10^{-5}$ mm$^{-1}$ [(e), (f)],
$1 \times 10^{-4}$ mm$^{-1}$ [(g), (h)],
$2 \times 10^{-4}$ mm$^{-1}$ [(i), (j)], and
$4 \times 10^{-4}$ mm$^{-1}$ [(k), (l)], at different reflection
numbers, $k = 1$ (black solid line), 5 (red long-dashed line),
10 (green short-dashed line), 20 (blue dotted line),
50 (magenta dot-dashed line), and 100 (cyan double dot-dashed
The laser pulse intensity is $1 \times 10^6$ W/cm$^{2}$ and the
pulse duration is 2 ns for all the graphs.

§ RESULTS AND DISCUSSIONS

At the outset, it should be noted that we confine our study to a
peak intensity, $I_0 = 1 \times 10^6$ W/cm$^2$ and a duration,
$\tau = 2$ ns of the laser pulse.
At resonance, the Doppler-averaged spin-flip probability of
the muonic hydrogen atoms at 10 K, is enhanced by a factor
of 2.3 as compared to that at 300 K. The lowest temperature
the system can go is 10 K as below which molecular $p \mu p$
are formed. Hence, we also confine
our study to the muonic hydrogen medium at a temperature of
10 K.
Additionally, at 10 K, the Doppler
broadening of the medium is about 100 MHz which is less than
half of the
bandwidth of a 2 ns pulse, namely
220 MHz. As the bandwidth of the pulse covers the Doppler width
of the medium fairly enough, we can safely calculate
the spin-flip probability at resonance without taking the
Doppler average and
apply a Gaussian window function with the parameters
equal to that of the laser pulse to effectively calculate the
Doppler-averaged spin-flip probability, $P$ at 10 K, which saves
computation time, without loss of accuracy. In Fig. <ref> (a),
we show the normalized Doppler distribution of muonic hydrogen
atoms at 300 K and 10 K, with blue and green solid lines,
respectively. The area under each of the lines are unity.
In Fig. <ref> (b), we show the normalized
Doppler distribution multiplied with a Gaussian window function
of width 220 MHz. The areas, $S_{10 \rm K}$ (= 0.917) and
$S_{300 \rm K}$ (= 0.387),
of the modified distribution lines in
Fig. <ref> (b) represent the number of atoms covered by the
bandwidth of the laser pulse at temperatures 10 K and 300 K,
respectively. The ratio of $S_{10 \rm K}$ and $S_{300\rm K}$
turns out
to be 2.3, which is exactly what we have calculated by using the
amplitude equations of the states and averaging over the
Doppler distribution of the atoms. This validates our proposed
method of calculation. Henceforth, in this work,
$P$ is simply equal to the population of the ground hyperfine
state, $1s(F=1)$.

We consider a 30 cm long cavity with a width of 10 mm which
accommodates approximately 100 reflections ($k = 100$) with an
angular incidence of the laser pulse at an angle of $10$ degrees
as clear from Fig. <ref>. It is more
convenient to use $k$ instead of the otherwise commonly used,
optical depth, $\alpha z$, where $z$ is the geometrical depth
that the pulse traverses inside the medium.
Typical value of the density of the muonic hydrogen atoms, $N$,
in the cavity, after the slow down process and losses by
transfer of muons
to heavier nucleus from the proton, is about
$10^5$ cm$^{-3}$ which is equivalent to a propagation
coefficient, $\alpha = 4 \times 10^{-18}$ mm$^{-1}$. At first
glance, this is extremely
low to cause any propagation effect on the laser pulse
for the ranges of optical depths the pulse traverses inside the
cavity. We calculate the representative values at the midpoint of
each path between the mirrors, corresponding to $k$.
In Fig. <ref>, we study the time evolution of the
population of the state $1s(F=1)$, $P$ and temporal variation of the
laser field, $\Omega$ for increasing values of
propagation coefficients, $\alpha$ down the columns.
In each panel, we show the representative values at $k = $1, 5,
10, 20, 50, and 100 with black solid, red long-dashed, green
short-dashed, blue dotted, magenta dot-dashed, and cyan double
dot-dashed lines, respectively. Panels (a), (c), (e), (g), (i),
and (k) [(b), (d), (f), (h), (j), and (l)] show the
temporal variation of
$P$ ($\Omega$) for $\alpha = 4 \times 10^{-18}$, $4 \times 10^{-6}$,
$4 \times 10^{-5}$, $1 \times 10^{-4}$, $2 \times 10^{-4}$, and
$4 \times 10^{-4}$ mm$^{-1}$, respectively. Figure <ref> (a)
shows that $P$ increases by two orders of magnitude
due to shrinking of the cross-section of the laser pulse.
It remains same through out the optical path of the pulse as
the laser field, $\Omega$ is unaffected due to the low muon
By observing the Eqs. <ref>-<ref>, we note that
it is the interplay between the pulse area $\sim \Omega \tau$ and
the optical path, $\alpha z$
which dictates the laser propagation effect on the dynamics of the

The results remain the same as Fig. <ref> (a) and (b) for
the range of values of $\alpha$ until it increases to
$4 \times 10^{-6}$ mm$^{-1}$ for which temporal variation of
$P$ and $\Omega$ just barely change with $k$. This is because,
in that region of $\alpha$,
$\Omega \tau \gg \alpha z$, hence $\Omega$ effectively remains
constant. This implies that in this range  of $\alpha$,
$P$ remains uniform through out the
cavity as the muonic hydrogen atoms are sparsely populated
inside the cavity. Therefore, we skip the results for those
values of
$\alpha$ and start from $\alpha = 4 \times 10^{-6}$ mm$^{-1}$
($N = 10^{17}$ cm$^{-3}$)
from Fig. <ref> (c) and (d). Undoubtedly, for
$\alpha = 4 \times 10^{-6}$ mm$^{-1}$,
$\alpha z$ begins to approach $\Omega \tau$,
$\alpha z = 0.0004 \Omega \tau$, to be specific.
In Fig. <ref> (e) and
(f) where $\alpha = 4 \times
10^{-5}$ mm$^{-1}$ ($N = 10^{17}$ cm$^{-3}$),
$P$ starts to fall appreciably below 0.0028 and
$\Omega$ begins to modulate and exhibit negative amplitude as the
optical depth
increases. As $\alpha$ further increases by a factor of two as
shown in the subsequent rows in Fig. <ref> (h)and (j),
$\Omega$ is significantly distorted during the propagation, and the
laser pulse breaks up into sub-pulses. For further higher value
of $\alpha$ ($ = 4 \times
10^{-6}$ mm$^{-1}$), as shown in Fig. <ref> (l), $\Omega$
breaks up into sub-pulses with positive and negative amplitudes
at $k = 100$. Similarly, as shown in Fig. <ref> (g), (i) and
(k), $P$ saturates at values lower than 0.0028 and exhibits
Rabi oscillation as optical depth increases.
We note that, if the pulse area and the pulse duration is much
shorter than the
spontaneous decay time which is 2.2 $\mu$s, the propagation
of the laser pulse follows a non-Beer's law decay profile and
propagates for much greater optical depths as compared to
CW laser fields. We expect more distortion if $\alpha$ increases
further. Clearly, the polarization of the muonic hydrogen
falls off with increased propagation depth or higher value of $k$,
with a wiggly behavior due to the modulation of $\Omega$ as
$k$ increases, as clear from Fig. <ref> (i) and (k) at
$t = 30$ s.

(Color online)
(a) Variation of the population of the ground hyperfine state, $P$,
(b) variation of the total efficiency, $Q$
with the reflection number, $k$ for medium propagation
coefficient, $\alpha = 4 \times 10^{-18}$ mm$^{-1}$ (black solid line),
$4 \times 10^{-6}$ mm$^{-1}$ (red long-dashed line),
$4 \times 10^{-5}$ mm$^{-1}$ (green short-dashed line),
$1 \times 10^{-4}$ mm$^{-1}$ (blue dotted line),
$2 \times 10^{-4}$ mm$^{-1}$ (magenta dot-dashed line), and
$4 \times 10^{-4}$ mm$^{-1}$ (cyan double dot-dashed line).

To get a clear profile of the polarization of muonic hydrogen atoms
inside the cavity, in Fig. <ref> (a) and (b), we plot
the variation of $P$ and the total efficiency, $Q$ as a function
of $k$, respectively. We define the total efficiency at a
given value of $k$
as the ratio of the total number of excited atoms to the total
number of atoms in the interaction volume between $k = 0$ to
the given value of $k$, i.e.,

\begin{align}
\label{eqn:effvsk}
Q & =  \frac{\int_k P_k N a \,dz}{\int_k N a \,dz}.
\end{align}

In the above equation, $P_k$ is the value of $P$ at the optical
depth corresponding to $k$, $N$ is the density of muonic hydrogen
atoms and $a$ is the beam cross-sectional area after shrinking
of the beam cross-section. In panel (a)[(b)], we plot
the values of $P$ ($Q$) for
$\alpha = 4 \times 10^{-18}$, $4 \times 10^{-6}$,
$4 \times 10^{-5}$, $1 \times 10^{-4}$, $2 \times 10^{-4}$, and
$4 \times 10^{-4}$ by black solid, red long-dashed, green
short-dashed, blue dotted, magenta dot-dashed, and cyan
double dot-dashed lines, respectively.
As shown in Fig. <ref> (a), $P$
falls off rapidly as
$\alpha$ increases beyond $4 \times 10^{-6}$ and has a wiggly
behavior for higher values of $\alpha$ as shown by the
magenta dot-dashed line and cyan double dot-dashed line in
Fig. <ref> (a). This is due to the oscillatory behavior of
$\Omega$ at higher values of $\alpha$. Total
efficiency, $Q$ also falls off in a similar trend but, the
wiggly behavior is absent due to the cumulative averaging of $P$
over the range of $k$ starting from the entrance of the cavity,
as shown in Fig. <ref> (b). $Q$ at a certain $k$ gives the
ratio of the total number of excited atoms to the total number of
atoms in a target volume accommodating $k$ reflections with the
beam cross-section shrunk by a factor $k$. Indeed,
$Q$ as a function
of $k$ provides information on how the total efficiency varies
as we shrink the laser beam cross-section and accordingly
increase the depth of the multipass cavity keeping the target
volume fixed.

Before closing the discussion, we note that upto a fairly
high density of the muonic hydrogen atoms for the given pulse
area of the laser, we can obtain a
uniformly polarized muonic beam with increased efficiency
by employing a multipass cavity.
The so-called threshold density of the muonic hydrogen will
further rise if we employ a pulse with higher pulse area.
By tailoring the pulse area
and repetition rate, we can construct the polarized
muonic beam to be continuous or periodic depending on the
requirement of the experiment. Keeping the pulse area lower than
the total optical depth, we can construct a pulsed polarized
muonic hydrogen beam. Lastly, although the
the spin-flip probability of the muonic hydrogen atoms
decreases down the cavity for higher densities of the atoms,
the total number of atoms excited in the interaction volume is
always higher for higher density of the atoms.

§ CONCLUSION

In conclusion, we have studied the prospects of enhancement of
the laser-induced spin-flip probability between the hyperfine
ground states of Doppler-broadened muonic hydrogen
by placing a multipass cavity at a temperature of 10 K. We have
investigated the effect of propagation of the laser pulse in the
medium of muonic hydrogen atoms inside the multipass cavity for
a range of
densities of the atoms. Assuming a pulse intensity of
$1 \times 10^6$ W/cm$^2$ and pulse duration of 2 ns, we considered
the case of shrinking the beam cross-section by a factor of 100
and employing a multipass cavity which allows 100 reflections of
the pulse, preserving the total interaction volume.
This transition is dipole-forbidden, and it is the magnetic dipole
interaction that induces the transition.
For densities of muonic hydrogen atoms ranging from $10^5$ to
$\sim 10^{17}$ cm$^{-3}$, the propagation effect is totally
negligible and the spin-flip probability remains uniform through
out the cavity and the total efficiency is 0.0028, i.e.
enhanced by two orders of magnitude. This gives rise to a uniformly
polarized muon beam.
As density increases beyond $\sim 10^{17}$ cm$^{-3}$,
the propagation
effect can no longer be neglected and spin-flip probability
is non-uniform inside the cavity. The total efficiency decreases
as the pulse undergoes more reflections inside the cavity and
eventually becomes 0.00041 after 100 reflections. This, in turn
gives rise to a pulsed polarized muonic beam. The result we
have presented will serve as a guideline for arranging
muonic beams depending on the requirement of a uniform or pulsed
polarized muons in $\mu$SR experiments.

§ ACKNOWLEDGEMENT

We acknowledge the helpful comments and suggestions of
Dr. Parvendra Kumar, IIT Chennai, India.

§ COMPETING INTERESTS

The authors declare that there they have no competing interests.

§ FUNDING

Part of this work was funded by Rashtriya Uchchatar Shiksha
Abhiyan (RUSA 2.0), MHRD, Govt. of India.

[1]
T.P. Gorringe, D.W. Hertzog,
“Precision muon physics,"
Prog. Part. Nucl. Phys. 84, 73 (2015).

[2]
R. Pohl, A. Antognini, F. Nez, et al.,
“The size of the proton,"
Nature 466, 213 (2010).

[3]
A. Antognini, F. Nez, et al.,
“Proton Structure from the Measurement of 2S-2P
Transition Frequencies of Muonic Hydrogen,"
Science 339, 6118, 417-420 (2013).

[4]
P. Bakule, and E. Morenzoni,
“Generation and applications of slow polarized muons,"
Contemporary Physics 45:3, 203-225 (2004).

[5]
J. P. Miller, E. de Rafael, and B. L. Roberts,
“Muon (g-2) experiment and theory,"
Rep. Prog. Phys. 70, 795-881 (2007).

[6]
T. Lancaster, S. J. Blundell, P. J. Baker, M. L. Brooks, W. Hayes, F. L.
Pratt, J. L. Manson, M. M. Conner, and J. A. Schlueter,
“Muon-Fluorine Entangled States in Molecular Magnets,"
Phys. Rev. Lett. 99, 267601 (2007).

[7]
A. Bungau, R. Cywinski, C. Bungau, P. King, and J. S. Lord,
“Target optimization studies for surface muon production,"
Phys. Rev. ST Accel. Beams 17, 034701 (2014).

[8]
S. Nagamiya, K. Nagamine, O. Hashimoto, and T. Yamazaki,
“Negative-muon spin rotation at the oxygen site in paramagnetic
Phys. Rev. Lett. 35, 308–311 (1975).

[9]
P. Dalmas de Réotier and A. Yaouanc,
“Muon spin rotation and relaxation in magnetic materials,"
J. Phys. Condens. Matter 9, 9113–9166 (1997).

[10]
D. V. Balin, V. A. Ganzha, S. M. Kozlov, E. M. Maev, G. E. Petrov,
M. A. Soroka, G. N. Schapkin, G. G. Semenchuk, V. A. Trofimov,
A. A. Vasiliev, A. A. Vorobyov, N. I. Voropaev, C. Petitjean, B. Gartner,
B. Lauss, J. Marton, J. Zmeskal, T. Case, K. M. Crowe, P. Kammel,
F. J. Hartmann, M. P. Faifman,
“High precision study of muon catalyzed fusion in $D_2$ and $HD$ gas,"
Phys. Particles and Nuclei 42, 185–214 (2011).

[11]
K. Nagamine, Y. Miyake, K. Shimomura, P. Birrer, J. P. Marangos, M. Iwasaki,
P. Strasser, and T. Kuga
“Ultraslow Positive-Muon Generation by Laser Ionization of Thermal
Muonium from Hot Tungsten at Primary Proton Beam,"
Phys. Rev. Lett. 74, 4811 (1995).

[12]
E. Morenzoni, F. Kottmann, D. Maden, B. Matthias, M. Meyberg, Th. Prokscha,
Th. Wutzke, and U. Zimmermann,
“Generation of very slow polarized positive muons,"
Phys. Rev. Lett. 72, 2793 (1995).

[13]
P. Bakule, Y. Matsuda, Y. Miyake, P. Strasser, K. Shimomura, S.
Makimura, and K. Nagamine,
“Slow muon experiment by laser resonant ionization method at
RIKEN-RAL muon facility,"
Spectrochim. Acta Part B 58, 1019–1030 (2003).

[14]
G. A. Beer, Y. Fujiwara, S. Hirota, K. Ishida, M. Iwasaki,
S. Kanda, H.
Kawai, N. Kawamura, R. Kitamura, S. Lee, W. Lee, G. M. Marshall, T.
Mibe, Y. Miyake, S. Okada, K. Olchanski, A. Olin, H. Ohnishi,
Y. Oishi,
M. Otani, N. Saito, K. Shimomura, P. Strasser, M. Tabata,
D. Tomono,
K. Ueno, E. Won, and K. Yokoyama,
“Enhancement of muonium emission rate from silica aerogel with a
laser-ablated surface,"
Prog. Theor. Exp. Phys. 2014, 091C01 (2014).

[15]
T. Nakajima,
“A scheme to polarize nuclear-spin of atoms by a sequence of short
laser pulses: application to the muonium,"
Opt. Express 18, 27468–27480 (2010).

[16]
T. Nakajima,
“Spin polarization of Doppler-broadened atoms by the
broadband nanosecond and transform-limited picosecond laser pulses:
case study for the muonium,"
J. Opt. Soc. Am. B 29, 2420–2424 (2012).

[17]
R. M. Das, S. Chatterjee, M. Iwasaki, and T. Nakajima,
“Ionization efficiency of Doppler-broadened atoms by
transform-limited and broadband nanosecond pulses: one-photon
resonant two-photon ionization of muoniums,"
J. Opt. Soc. Am. B 32, 1237–1244 (2015).

[18]
R. M. Das, K. Ishida, M. Iwasaki, and T. Nakajima,
“Nuclear-spin polarization of atoms by chirped laser
pulses: application to the muonium,"
J. Opt. Soc. Am. B 35, 1799–1810 (2018).

[19]
R. M. Das, K. Ishida, M. Iwasaki, and T. Nakajima,
“Excitation probability of the hyperfine ground
state of Doppler-broadened muonic hydrogen through
dipole-forbidden transition,"
Optik 200, 163380 (2020).

[20]
S. Yamamoto, K. Ninomiya, N. Kawamura,
T. Yabe, and Y. Hirano,
“Three-dimensional (3D) optical imaging
of muon beam using a plastic scintillator plate in water,"
Nucl. Instrum. Meth. A 1015, 165768 (2021).

[21]
L. Serafini and I. Drebot and A. Bacci and F. Broggi and C. Curatolo and A. Marocchino and N. Panzeri and V. Petrillo and A.R. Rossi and M. Rossetti Conti,
“A muon source based on plasma accelerators,"
Nucl. Instrum. Meth. A 909, 309-313 (2018).

[22]
The future prospects of muon colliders and neutrino factories.

[23]
A. Adamczak, D. Bakalov, L. Stoychev, and A. Vacchi,
“Hyperfine spectroscopy of muonic hydrogen and the PSI
Lamb shift experiment,"
Nucl. Instrum. Meth. B 281, 72-76 (2012).

[24]
J.R. Sapirstein, J.R. Yennie, in: T. Kinoshita,
Quantum Electrodynamics, World Scientific, 560 (1990).

[25]
M.I. Eides, H. Grotch, V.A. Shelyuto,
“Theory of light hydrogenlike atoms."
Phys. Rep. 342, 62 (2001).

[26]
D. Bakalov, E. Milotti, C. Rizzo, A. Vacchi, and E. Zavattini,
“Experimental method to measure the hyperfine splitting of muonic
hydrogen ($\mu p$)1S,"
Phys. Lett. A 172, 277 (1993).

[27]
F. Foroughi, E. Morenzoni, T. Prokscha, et al.,
“Upgrading the PSI Muon Facility,"
Hyperfine Interactions 138, 483–-488 (2001).

[28]
H.K. Walter,
“Muon Physics at the Paul Scherrer Institut
(psi) and at Triumf,"
High Intensity Muon Sources - Kek International
Workshop, World Scientific Publishing Co. Pte. Ltd, 279–290

[29]
T. Matsuzaki, K. Ishida, K. Nagamine, I. Watanabe, G.H. Eaton,
and W.G. Williams,
“The RIKEN-RAL pulsed Muon Facility,"
Nucl. Instrum. Meth. A 465, 365–383 (2001).

[30]
Y. Miyake, K. Shimomura, N. Kawamura,
P. Strasser, A. Koda, S. Makimura,
H. Fujimori, Y. Ikedo, K. Nakahara,
S. Takeshita, M. Kato, K. Kojima,
Y. Kobayashi, K. Nishiyama, R. Kadono,
W. Higemoto, T. U. Ito, K. Ninomiya,
K. Kubo, and K. Nagamine,
“J-PARC Muon Facility, MUSE,"
Phys. Proc. 30, 46–49 (2012).

[31]
G. Buica, T. Nakajima,
“Propagation of two short laser pulse trains in a
$\Lambda$-type three-level medium under conditions of
electromagnetically induced transparency,"
Opt. Comm. 332, 59–69 (2014).


# Reinforcement Learning for Mixed-Integer Problems Based on MPC

Sebastien Gros  Mario Zanon  Norwegian University of Technology, NTNU  IMT
School for Advanced Studies Lucca

###### Abstract

Model Predictive Control has been recently proposed as policy approximation
for Reinforcement Learning, offering a path towards safe and explainable
Reinforcement Learning. This approach has been investigated for $Q$-learning
and actor-critic methods, both in the context of nominal Economic MPC and
Robust (N)MPC, showing very promising results. In that context, actor-critic
methods seem to be the most reliable approach. Many applications include a
mixture of continuous and integer inputs, for which the classical actor-critic
methods need to be adapted. In this paper, we present a policy approximation
based on mixed-integer MPC schemes, and propose a computationally inexpensive
technique to generate exploration in the mixed-integer input space that
ensures a satisfaction of the constraints. We then propose a simple compatible
advantage function approximation for the proposed policy, that allows one to
build the gradient of the mixed-integer MPC-based policy.

###### keywords:

Reinforcement Learning, Mixed-Integer Model Predictive Control, actor-critic
methods, stochastic and deterministic policy gradient.

##  1 Introduction

Reinforcement Learning (RL) is a powerful tool for tackling stochastic
processes without depending on a detailed model of the probability
distributions underlying the state transitions. Indeed, most RL methods rely
purely on observed data, and realizations of the stage cost assessing the
system performance. RL methods seek to increase the closed-loop performance of
the control policy deployed on the system as observations are collected. RL
has drawn an increasingly large attention thanks to its accomplishments, such
as, e.g., making it possible for robots to learn to walk or fly from
experiments (Wang et al., 2012; Abbeel et al., 2007).

Most RL methods are based on learning the optimal control policy for the real
system either directly, or indirectly. Indirect methods typically rely on
learning a good approximation of the optimal action-value function underlying
the system. The optimal policy is then indirectly obtained as the minimizer of
the value-function approximation over the inputs. Direct RL methods, if based
on the policy gradient, seek to adjust the parameters
$\boldsymbol{\mathrm{\theta}}$ of a given policy
$\pi_{\boldsymbol{\mathrm{\theta}}}$ such that it yields the best closed-loop
performance when deployed on the real system. An attractive advantage of
direct RL methods over indirect ones is that they are based on formal
necessary conditions of optimality for the closed-loop performance of
$\pi_{\boldsymbol{\mathrm{\theta}}}$, and therefore guarantee - for a large
enough data set - the (possibly local) assymptotic optimality of the
parameters $\boldsymbol{\mathrm{\theta}}$ (Sutton et al., 1999; Silver et al.,
2014).

RL methods often rely on Deep Neural Networks (DNN) to carry the policy
approximation $\pi_{\boldsymbol{\mathrm{\theta}}}$. Unfortunately, control
policies based on DNNs provide limited opportunities for formal verifications
of the resulting policy, and for imposing hard constraints on the evolution of
the state of the real system. The development of safe RL methods, which aims
at tackling this issue, is currently an open field of research (J. Garcia,
2013). A novel approach towards providing formal safety certificates in the
context of RL has been recently proposed in (Gros and Zanon, 2019, 2020; Zanon
and Gros, 2019), where the policy approximation is based on robust Model
Predictive Control (MPC) schemes rather than unstructured function
approximators like DNNs. The validity of this choice is discussed in details
in (Gros and Zanon, 2019). In (Gros and Zanon, 2020), methodologies to deploy
direct RL techniques on MPC-based policy approximations are proposed. These
methodologies are, however, restricted to continuous input spaces and
therefore exclude integer decision variables, which are central in a number of
applications.

In this paper, we propose an extension of the policy gradient techniques
proposed in (Gros and Zanon, 2020) to mixed-integer problems. A mixed-integer
MPC is used as a policy approximation, and a policy gradient method adjusts
the MPC parameters for closed-loop performance. We detail how the actor-critic
method can be deployed in this specific context. In particular, we propose an
asymptotically exact hybrid stochastic-deterministic policy approach allowing
for computing the policy gradient at a lower computational complexity than a
full stochastic approach. We then propose a hybrid compatible advantage-
function approximator tailored to our formulation. We finally detail how the
mixed-integer MPC can be differentiated at a low computational cost, using
principles from parametric Nonlinear Programming, in order to implement the
actor-critic method. The proposed method is illustrated on a simple example,
allowing for an unambiguous presentation of the results.

The paper is structured as follows. Section 2 provides background material on
MDPs and RL. Section 3 presents the construction of a mixed-integer stochastic
policy using a mixed-integer MPC scheme to support the policy approximation.
Section 4 details an actor-critic method tailored to the proposed formulation,
and how the policy gradient can be estimated. A compatible advantage function
approximation is proposed. Section 5 details how the mixed-integer MPC scheme
can be efficiently differentiated. Section 6 proposes an illustrative example,
and Section 7 provides some discussions.

##  2 Background

In the following, we will consider that the dynamics of the real system are
described as a stochastic process on (possibly) continuous state-input spaces.
We will furthermore consider (possibly) stochastic policies $\pi$, taking the
form of probability densities:

$\displaystyle{\pi}\left[\boldsymbol{\mathrm{a}}\,|\,\boldsymbol{\mathrm{s}}\right]\,:\,\mathbb{R}^{m}\times\mathbb{R}^{n}\rightarrow\mathbb{R}_{+},$
(1)

denoting the probability density of selecting a given input
$\boldsymbol{\mathrm{a}}$ when the system is in a given state
$\boldsymbol{\mathrm{s}}$. Deterministic policies delivering
$\boldsymbol{\mathrm{a}}$ as a function of $\boldsymbol{\mathrm{s}}$ will be
labelled as:

$\displaystyle\boldsymbol{\mathrm{\pi}}\left(\boldsymbol{\mathrm{s}}\right)\,:\,\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}.$
(2)

Any deterministic policy can be viewed as a stochastic one, having a Dirac
function as a probability density (or unit function for discrete inputs),
i.e.,
${\pi}\left[\,\boldsymbol{\mathrm{a}}\,|\,\boldsymbol{\mathrm{s}}\,\right]=\delta\left(\boldsymbol{\mathrm{a}}-\boldsymbol{\mathrm{\pi}}\left(\boldsymbol{\mathrm{s}}\right)\right).$

We consider a stage cost function
$L(\boldsymbol{\mathrm{s}},\boldsymbol{\mathrm{a}})\in\mathbb{R}$ and a
discount factor $\gamma\in[0,1]$, the performance of a policy $\pi$ is
assessed via the total expected cost:

$\displaystyle
J(\pi)=\mathbb{E}_{{\pi}}\left[\left.\sum_{k=0}^{\infty}\,\gamma^{k}L(\boldsymbol{\mathrm{s}}_{k},\boldsymbol{\mathrm{a}}_{k})\,\right|\,\boldsymbol{\mathrm{a}}_{k}\sim\pi\left[\,.\,|\,\boldsymbol{\mathrm{s}}_{k}\right]\,\right].$
(3)

The optimal policy associated to the state transition, the stage cost $L$ and
the discount factor $\gamma$ is deterministic and given by:

$\displaystyle\pi_{\star}=\mathrm{arg}\min_{\pi}\,J(\pi).$ (4)

The value function $V_{\boldsymbol{\mathrm{\pi}}}$, action-value function
$Q_{\boldsymbol{\mathrm{\pi}}}$ and advantage functions
$A_{\boldsymbol{\mathrm{\pi}}}$ associated to a given policy $\pi$ are given
by (Bertsekas, 1995; Bertsekas and Shreve, 1996; Bertsekas, 2007):

$\displaystyle
V_{\boldsymbol{\mathrm{\pi}}}\left(\boldsymbol{\mathrm{s}}\right)$
$\displaystyle=\mathbb{E}\left[L(\boldsymbol{\mathrm{s}},\boldsymbol{\mathrm{a}})+\gamma
V_{\boldsymbol{\mathrm{\pi}}}(\boldsymbol{\mathrm{s}}_{+})\,|\,\boldsymbol{\mathrm{s}},\,\boldsymbol{\mathrm{a}}\right],$
(5a) $\displaystyle
Q_{\boldsymbol{\mathrm{\pi}}}\left(\boldsymbol{\mathrm{s}},\boldsymbol{\mathrm{a}}\right)$
$\displaystyle=L(\boldsymbol{\mathrm{s}},\boldsymbol{\mathrm{a}})+\gamma\mathbb{E}\left[V_{\boldsymbol{\mathrm{\pi}}}(\boldsymbol{\mathrm{s}}_{+})\,|\,\boldsymbol{\mathrm{s}},\,\boldsymbol{\mathrm{a}}\right],$
(5b) $\displaystyle
A_{\boldsymbol{\mathrm{\pi}}}\left(\boldsymbol{\mathrm{s}},\boldsymbol{\mathrm{a}}\right)$
$\displaystyle=Q_{\boldsymbol{\mathrm{\pi}}}\left(\boldsymbol{\mathrm{s}},\boldsymbol{\mathrm{a}}\right)-V_{\boldsymbol{\mathrm{\pi}}}\left(\boldsymbol{\mathrm{s}}\right),$
(5c)

where the expected value in (5b) is taken over the state transition, and the
one in (5a) is taken over the state transitions and (1).

###  2.1 Stochastic policy gradient

In most cases, the optimal policy $\pi_{\star}$ cannot be computed, either
because the system is not exactly known or because solving (5) is too
expensive. It is then useful to consider approximations
$\pi_{\boldsymbol{\mathrm{\theta}}}$ of the optimal policy, parametrized by
$\boldsymbol{\mathrm{\theta}}$. The optimal parameters
$\boldsymbol{\mathrm{\theta}}_{\star}$ are then given by:

$\displaystyle\boldsymbol{\mathrm{\theta}}_{\star}=\mathrm{arg}\min_{\boldsymbol{\mathrm{\theta}}}\,J(\pi_{\boldsymbol{\mathrm{\theta}}}).$
(6)

The policy gradient
$\nabla_{\boldsymbol{\mathrm{\theta}}}\,J({\pi}_{\boldsymbol{\mathrm{\theta}}})$
associated to the stochastic policy $\pi_{\boldsymbol{\mathrm{\theta}}}$ is
then instrumental in finding $\boldsymbol{\mathrm{\theta}}_{\star}$ by taking
gradient steps in $\boldsymbol{\mathrm{\theta}}$. The policy gradient can be
obtained using various actor-critic methods (Sutton and Barto, 1998; Sutton et
al., 1999). In this paper, we will use the actor-critic formulation:

$\displaystyle\nabla_{\boldsymbol{\mathrm{\theta}}}\,J({\pi}_{\boldsymbol{\mathrm{\theta}}})$
$\displaystyle=\mathbb{E}_{{\boldsymbol{\mathrm{\pi}}_{\boldsymbol{\mathrm{\theta}}}}}\left[\nabla_{\boldsymbol{\mathrm{\theta}}}\log{\pi}_{\boldsymbol{\mathrm{\theta}}}\,A_{{\pi}_{\boldsymbol{\mathrm{\theta}}}}\right],$
(7)

for stochastic policies, and the actor-critic formulation:

$\displaystyle\nabla_{\boldsymbol{\mathrm{\theta}}}\,J({\boldsymbol{\mathrm{\pi}}}_{\boldsymbol{\mathrm{\theta}}})$
$\displaystyle=\mathbb{E}_{{\boldsymbol{\mathrm{\pi}}_{\boldsymbol{\mathrm{\theta}}}}}\left[\nabla_{\boldsymbol{\mathrm{\theta}}}{\boldsymbol{\mathrm{\pi}}}_{\boldsymbol{\mathrm{\theta}}}\,\nabla_{\boldsymbol{\mathrm{a}}}A_{{\boldsymbol{\mathrm{\pi}}}_{\boldsymbol{\mathrm{\theta}}}}\right],$
(8)

for deterministic policies.

The value functions $V_{\boldsymbol{\mathrm{\pi}}}$,
$Q_{\boldsymbol{\mathrm{\pi}}}$ and $A_{\boldsymbol{\mathrm{\pi}}}$ associated
to a given policy $\pi$ are typically evaluated via Temporal-Difference (TD)
techniques (Sutton and Barto, 1998), and require that a certain amount of
exploration is included in the deployment of the policy. For deterministic
policies, the exploration can, e.g., be generated by including stochastic
perturbations over the policy
$\boldsymbol{\mathrm{\pi}}_{\boldsymbol{\mathrm{\theta}}}$, while stochastic
policies generate exploration by construction. Note that is is fairly common
in RL to define the stochastic policy $\pi_{\boldsymbol{\mathrm{\theta}}}$ as
an arbitrary density, e.g., the normal distribution, centered at a
deterministic policy
$\boldsymbol{\mathrm{\pi}}_{\boldsymbol{\mathrm{\theta}}}$. We shall observe
here that the deterministic policy gradient (8) is not suited as such for
integer inputs, as the gradients
$\nabla_{\boldsymbol{\mathrm{\theta}}}{\boldsymbol{\mathrm{\pi}}}_{\boldsymbol{\mathrm{\theta}}},\,\nabla_{\boldsymbol{\mathrm{a}}}A_{{\boldsymbol{\mathrm{\pi}}}_{\boldsymbol{\mathrm{\theta}}}}$
do not exist on discrete input spaces. On continuous input spaces, the choice
between the deterministic approach (8) or the stochastic approach (7) is
typically motivated by computational aspects.

##  3 Mixed-integer Optimization-based policy

In this paper, we will consider parametrized deterministic policies
$\boldsymbol{\mathrm{\pi}}_{\boldsymbol{\mathrm{\theta}}}\approx\boldsymbol{\mathrm{\pi}}_{\star}$
based on parametric optimization problems. In particular, we will focus on
optimization problems resulting from a nominal mixed-integer MPC formulation.
The results proposed in this paper extend to robust MPC - enabling the
construction of safe Reinforcement Learning methods - but this case is omitted
in this paper for the sake of brevity.

###  3.1 Policy approximation based on mixed-integer MPC

The mixed-integer MPC scheme reads as:

$\displaystyle\boldsymbol{\mathrm{u}}^{\star}\left(\boldsymbol{\mathrm{s}},\boldsymbol{\mathrm{\theta}}\right),$
$\displaystyle\,\boldsymbol{\mathrm{i}}^{\star}\left(\boldsymbol{\mathrm{s}},\boldsymbol{\mathrm{\theta}}\right)=$
$\displaystyle\mathrm{arg}\min_{\boldsymbol{\mathrm{u}},\boldsymbol{\mathrm{i}}}$
$\displaystyle\quad
T(\boldsymbol{\mathrm{x}}_{N},\boldsymbol{\mathrm{\theta}})+\sum_{k=0}^{N-1}\ell(\boldsymbol{\mathrm{x}}_{k},\boldsymbol{\mathrm{u}}_{k},\boldsymbol{\mathrm{i}}_{k},\boldsymbol{\mathrm{\theta}})$
(9a) $\displaystyle\mathrm{s.t.}$
$\displaystyle\quad\boldsymbol{\mathrm{x}}_{k+1}=\boldsymbol{\mathrm{F}}\left(\boldsymbol{\mathrm{x}}_{k},\boldsymbol{\mathrm{u}}_{k},\boldsymbol{\mathrm{i}}_{k},\boldsymbol{\mathrm{\theta}}\right),\,\,\,\boldsymbol{\mathrm{x}}_{0}=\boldsymbol{\mathrm{s}},$
(9b)
$\displaystyle\quad\boldsymbol{\mathrm{h}}_{k}\left(\boldsymbol{\mathrm{x}}_{k},\boldsymbol{\mathrm{u}}_{k},\boldsymbol{\mathrm{i}}_{k},\boldsymbol{\mathrm{\theta}}\right)\leq
0,\quad k=0,...,N-1,$ (9c)
$\displaystyle\quad\boldsymbol{\mathrm{h}}_{N}\left(\boldsymbol{\mathrm{x}}_{N},\boldsymbol{\mathrm{\theta}}\right)\leq
0,$ (9d)
$\displaystyle\quad\boldsymbol{\mathrm{i}}_{k}\in\left\\{0,1\right\\}^{m_{\mathrm{i}}},$
(9e)

where $\boldsymbol{\mathrm{x}}_{k}\in\mathbb{R}^{n}$ are the predicted system
trajectories, $\boldsymbol{\mathrm{u}}_{k}\in\mathbb{R}^{m_{\mathrm{c}}}$ the
planned continuous inputs and
$\boldsymbol{\mathrm{i}}_{k}\in\left\\{0,1\right\\}^{m_{\mathrm{i}}}$ the
planned integer inputs. Without loss of generality, we consider binary integer
inputs. Functions $\ell$, $T$ are the stage and terminal costs. Functions
$\boldsymbol{\mathrm{h}}_{0,\ldots,N-1}$ are the stage constraints and
function $\boldsymbol{\mathrm{h}}_{N}$ is the terminal constraint.

For a given state $\boldsymbol{\mathrm{s}}$ and parameters
$\boldsymbol{\mathrm{\theta}}$, the MPC scheme (9) delivers the continuous and
integer input profiles

$\displaystyle\boldsymbol{\mathrm{u}}^{\star}\left(\boldsymbol{\mathrm{s}},\boldsymbol{\mathrm{\theta}}\right)$
$\displaystyle=\left\\{\boldsymbol{\mathrm{u}}_{0}^{\star}\left(\boldsymbol{\mathrm{s}},\boldsymbol{\mathrm{\theta}}\right),\ldots,\boldsymbol{\mathrm{u}}_{N-1}^{\star}\left(\boldsymbol{\mathrm{s}},\boldsymbol{\mathrm{\theta}}\right)\right\\},$
(10a)
$\displaystyle\boldsymbol{\mathrm{i}}^{\star}\left(\boldsymbol{\mathrm{s}},\boldsymbol{\mathrm{\theta}}\right)$
$\displaystyle=\left\\{\boldsymbol{\mathrm{i}}_{0}^{\star}\left(\boldsymbol{\mathrm{s}},\boldsymbol{\mathrm{\theta}}\right),\ldots,\boldsymbol{\mathrm{i}}_{N-1}^{\star}\left(\boldsymbol{\mathrm{s}},\boldsymbol{\mathrm{\theta}}\right)\right\\},$
(10b)

with
$\boldsymbol{\mathrm{u}}_{k}^{\star}\left(\boldsymbol{\mathrm{s}},\boldsymbol{\mathrm{\theta}}\right)\in\mathbb{R}^{m_{\mathrm{c}}}$
and
$\boldsymbol{\mathrm{i}}_{k}^{\star}\left(\boldsymbol{\mathrm{s}},\boldsymbol{\mathrm{\theta}}\right)\in\left\\{0,1\right\\}^{m_{\mathrm{i}}}$.
The MPC scheme (9) generates a parametrized deterministic policy

$\displaystyle\boldsymbol{\mathrm{\pi}}_{\boldsymbol{\mathrm{\theta}}}\left(\boldsymbol{\mathrm{s}}\right)=\left\\{\boldsymbol{\mathrm{\pi}}^{\mathrm{c}}_{\boldsymbol{\mathrm{\theta}}}\left(\boldsymbol{\mathrm{s}}\right),\,\boldsymbol{\mathrm{\pi}}^{\mathrm{i}}_{\boldsymbol{\mathrm{\theta}}}\left(\boldsymbol{\mathrm{s}}\right)\right\\},$
(11)

where

$\displaystyle\boldsymbol{\mathrm{\pi}}^{\mathrm{c}}_{\boldsymbol{\mathrm{\theta}}}\left(\boldsymbol{\mathrm{s}}\right)$
$\displaystyle=\boldsymbol{\mathrm{u}}_{0}^{\star}\left(\boldsymbol{\mathrm{s}},\boldsymbol{\mathrm{\theta}}\right)\,\in\,\mathbb{R}^{m_{\mathrm{c}}},$
(12a)
$\displaystyle\boldsymbol{\mathrm{\pi}}^{\mathrm{i}}_{\boldsymbol{\mathrm{\theta}}}\left(\boldsymbol{\mathrm{s}}\right)$
$\displaystyle=\boldsymbol{\mathrm{i}}_{0}^{\star}\left(\boldsymbol{\mathrm{s}},\boldsymbol{\mathrm{\theta}}\right)\,\in\,\left\\{0,1\right\\}^{m_{\mathrm{i}}},$
(12b)

are the first elements of the continuous and integer input sequences generated
by (9). In the following, it will be useful to consider the MPC scheme (9) as
a generic parametric mixed-integer NLP:

$\displaystyle\boldsymbol{\mathrm{u}}^{\star}\left(\boldsymbol{\mathrm{s}},\boldsymbol{\mathrm{\theta}}\right),\,\boldsymbol{\mathrm{i}}^{\star}\left(\boldsymbol{\mathrm{s}},\boldsymbol{\mathrm{\theta}}\right)=\mathrm{arg}\min_{\boldsymbol{\mathrm{u}},\,\boldsymbol{\mathrm{i}}}$
$\displaystyle\quad\Phi(\boldsymbol{\mathrm{x}},\boldsymbol{\mathrm{u}},\boldsymbol{\mathrm{i}},\boldsymbol{\mathrm{\theta}})$
(13a) $\displaystyle\mathrm{s.t.}$
$\displaystyle\quad\boldsymbol{\mathrm{f}}\left(\boldsymbol{\mathrm{x}},\boldsymbol{\mathrm{u}},\boldsymbol{\mathrm{i}},\boldsymbol{\mathrm{s}},\boldsymbol{\mathrm{\theta}}\right)=0,$
(13b)
$\displaystyle\quad\boldsymbol{\mathrm{h}}\left(\boldsymbol{\mathrm{x}},\boldsymbol{\mathrm{u}},\boldsymbol{\mathrm{i}},\boldsymbol{\mathrm{\theta}}\right)\leq
0,$ (13c)
$\displaystyle\quad\boldsymbol{\mathrm{i}}\in\left\\{0,1\right\\}^{m_{\mathrm{i}}\times
N-1},$ (13d)

where function $\Phi$ gathers the stage and terminal cost functions from (9a),
function $\boldsymbol{\mathrm{f}}$ gathers the dynamic constraints and initial
conditions (9b), and function $\boldsymbol{\mathrm{h}}$ gathers the stage and
terminal constraints (9c)-(9d).

##  4 actor-critic method

In order to build actor-critic methods for (11), exploration is required
(Sutton and Barto, 1998). When the input space is constrained and mixed-
integer, the exploration becomes non-trivial to setup, as 1. it must retain
the feasibility of the hard constraints (9c)-(9d) and 2. simple input
disturbances are not possible for the integer part since they are locked on an
integer grid. To address this issue, we will adopt a stochastic policy
approach, well suited for the integer part, and consider its asymptotically
equivalent deterministic counterpart on the continuous input space, well
suited for computational efficiency.

###  4.1 MPC-based exploration

In order to generate exploration, we will build a stochastic policy (1) based
on the deterministic policy (11) where $\boldsymbol{\mathrm{a}}$ will gather
the continuous inputs $\boldsymbol{\mathrm{a}}^{\mathrm{c}}$ and integer
inputs $\boldsymbol{\mathrm{a}}^{\mathrm{i}}$ actually applied to the real
system, i.e.,
$\boldsymbol{\mathrm{a}}=\left\\{\,\boldsymbol{\mathrm{a}}^{\mathrm{c}},\,\boldsymbol{\mathrm{a}}^{\mathrm{i}}\,\right\\}$.
We will build (1) such that it generates exploration that is respecting the
constraints (9c)-(9d) with unitary probability. We propose to build (1) such
that it becomes naturally separable between the integer and continuous part in
the policy gradient computation. To that end, we consider a softmax approach
to handle the integer part of the problem. More specifically, we consider the
parametric mixed-integer NLP:

$\displaystyle\Phi^{\mathrm{i}}(\boldsymbol{\mathrm{s}},\boldsymbol{\mathrm{\theta}},\boldsymbol{\mathrm{a}}^{\mathrm{i}})=\min_{\boldsymbol{\mathrm{u}},\,\boldsymbol{\mathrm{i}}}$
$\displaystyle\quad\Phi(\boldsymbol{\mathrm{x}},\boldsymbol{\mathrm{u}},\boldsymbol{\mathrm{i}},\boldsymbol{\mathrm{\theta}})$
(14a) $\displaystyle\mathrm{s.t.}$
$\displaystyle\quad\boldsymbol{\mathrm{f}}\left(\boldsymbol{\mathrm{x}},\boldsymbol{\mathrm{u}},\boldsymbol{\mathrm{i}},\boldsymbol{\mathrm{s}},\boldsymbol{\mathrm{\theta}}\right)=0,$
(14b)
$\displaystyle\quad\boldsymbol{\mathrm{h}}\left(\boldsymbol{\mathrm{x}},\boldsymbol{\mathrm{u}},\boldsymbol{\mathrm{i}},\boldsymbol{\mathrm{\theta}}\right)\leq
0,$ (14c)
$\displaystyle\quad\boldsymbol{\mathrm{i}}_{0}=\boldsymbol{\mathrm{a}}^{\mathrm{i}},$
(14d)
$\displaystyle\quad\boldsymbol{\mathrm{i}}_{1,\ldots,N-1}\in\left\\{0,1\right\\}^{m_{\mathrm{i}}},$
(14e)

derived from (13), where the first integer input is assigned to
$\boldsymbol{\mathrm{a}}^{\mathrm{i}}$ via constraint (14d). We will consider
that
$\Phi^{\mathrm{i}}(\boldsymbol{\mathrm{s}},\boldsymbol{\mathrm{\theta}},\boldsymbol{\mathrm{a}}^{\mathrm{i}})$
takes infinite value when the selected integer input
$\boldsymbol{\mathrm{a}}^{\mathrm{i}}$ is infeasible. Let us label
$\mathbb{I}(\boldsymbol{\mathrm{s}},\boldsymbol{\mathrm{\theta}})$ the
feasible set of $\boldsymbol{\mathrm{a}}^{\mathrm{i}}$ for a given state
$\boldsymbol{\mathrm{s}}$ and MPC parameter $\boldsymbol{\mathrm{\theta}}$,
and
$\tilde{\boldsymbol{\mathrm{i}}}(\boldsymbol{\mathrm{s}},\boldsymbol{\mathrm{\theta}},\boldsymbol{\mathrm{a}}^{\mathrm{i}})$
the integer profile solution of (14). By construction
$\tilde{\boldsymbol{\mathrm{i}}}_{0}(\boldsymbol{\mathrm{s}},\boldsymbol{\mathrm{\theta}},\boldsymbol{\mathrm{a}}^{\mathrm{i}})=\boldsymbol{\mathrm{a}}^{\mathrm{i}}$
when
$\boldsymbol{\mathrm{a}}^{\mathrm{i}}\in\mathbb{I}(\boldsymbol{\mathrm{s}},\boldsymbol{\mathrm{\theta}})$.
We then define the softmax stochastic integer policy distribution using

$\displaystyle\pi_{\boldsymbol{\mathrm{\theta}}}^{\mathrm{i}}\left[\,\boldsymbol{\mathrm{a}}^{\mathrm{i}}\,|\,\boldsymbol{\mathrm{s}}\,\right]\,\propto\,e^{-\sigma_{\mathrm{i}}^{-1}\Phi^{\star}_{\mathrm{i}}(\boldsymbol{\mathrm{s}},\boldsymbol{\mathrm{\theta}},\boldsymbol{\mathrm{a}}^{\mathrm{i}})}\in\mathbb{R}_{+},$
(15)

where $\sigma_{\mathrm{i}}>0$ is a parameter adjusting the variance of
$\pi_{\boldsymbol{\mathrm{\theta}}}^{\mathrm{i}}$. In order to build the
continuous part of the policy, we will consider the continuous part
$\boldsymbol{\mathrm{a}}^{\mathrm{c}}$ of the stochastic policy as conditioned
on $\tilde{\boldsymbol{\mathrm{i}}}$, and taking the form of a probability
density:

$\displaystyle\pi^{\mathrm{c}}_{\boldsymbol{\mathrm{\theta}}}\left[\,\boldsymbol{\mathrm{a}}^{\mathrm{c}}\,\right|\,\tilde{\boldsymbol{\mathrm{i}}}(\boldsymbol{\mathrm{s}},\boldsymbol{\mathrm{\theta}},\boldsymbol{\mathrm{a}}^{\mathrm{i}}),\,\boldsymbol{\mathrm{s}}\,\left.\right]\in\mathbb{R}_{+},$
(16)

which will be constructed from the parametric NLP:

$\displaystyle\tilde{\boldsymbol{\mathrm{u}}}(\boldsymbol{\mathrm{s}},\boldsymbol{\mathrm{\theta}},\boldsymbol{\mathrm{i}},\boldsymbol{\mathrm{d}})=\mathrm{arg}\min_{\boldsymbol{\mathrm{u}}}$
$\displaystyle\quad\Phi(\boldsymbol{\mathrm{x}},\boldsymbol{\mathrm{u}},\boldsymbol{\mathrm{i}},\boldsymbol{\mathrm{\theta}})+\boldsymbol{\mathrm{d}}^{\top}\boldsymbol{\mathrm{u}}_{0}$
(17a) $\displaystyle\mathrm{s.t.}$
$\displaystyle\quad\boldsymbol{\mathrm{f}}\left(\boldsymbol{\mathrm{x}},\boldsymbol{\mathrm{u}},\boldsymbol{\mathrm{i}},\boldsymbol{\mathrm{s}},\boldsymbol{\mathrm{\theta}}\right)=0,$
(17b)
$\displaystyle\quad\boldsymbol{\mathrm{h}}\left(\boldsymbol{\mathrm{x}},\boldsymbol{\mathrm{u}},\boldsymbol{\mathrm{i}},\boldsymbol{\mathrm{\theta}}\right)\leq
0,$ (17c)

derived from (13), but where the integer input profile is entirely assigned,
and where $\boldsymbol{\mathrm{d}}\in\mathrm{R}^{m_{\mathrm{c}}}$ is a random
vector chosen as
$\boldsymbol{\mathrm{d}}\sim\mathcal{N}\left(0,\sigma_{\mathrm{c}}I\right)$.
The random variable $\boldsymbol{\mathrm{a}}_{\mathrm{c}}$ in (16) will then
be selected as:

$\displaystyle\boldsymbol{\mathrm{a}}^{\mathrm{c}}=\tilde{\boldsymbol{\mathrm{u}}}_{0}(\boldsymbol{\mathrm{s}},\boldsymbol{\mathrm{\theta}},\tilde{\boldsymbol{\mathrm{i}}}(\boldsymbol{\mathrm{s}},\boldsymbol{\mathrm{\theta}},\boldsymbol{\mathrm{a}}^{\mathrm{i}}),\boldsymbol{\mathrm{d}}).$
(18)

As previously observed in (Gros and Zanon, 2020), while
$\pi^{\mathrm{c}}_{\boldsymbol{\mathrm{\theta}}}$ is easy to sample, it is in
general difficult to evaluate.

Because $\boldsymbol{\mathrm{a}}^{\mathrm{c}}$ is conditioned on
$\tilde{\boldsymbol{\mathrm{i}}}$ and, therefore,
$\boldsymbol{\mathrm{a}}^{\mathrm{i}}$, the Kolmogorov definition of
conditional probabilities entails that the overall stochastic policy (1) reads
as the distribution:

$\displaystyle\pi_{\boldsymbol{\mathrm{\theta}}}[\,\boldsymbol{\mathrm{a}}\,|\,\boldsymbol{\mathrm{s}}\,]=\pi^{\mathrm{c}}_{\boldsymbol{\mathrm{\theta}}}\left[\,\boldsymbol{\mathrm{a}}^{\mathrm{c}}\,\left|\,\tilde{\boldsymbol{\mathrm{i}}}(\boldsymbol{\mathrm{s}},\boldsymbol{\mathrm{\theta}},\boldsymbol{\mathrm{a}}^{\mathrm{i}}),\,\boldsymbol{\mathrm{s}}\,\right.\right]\,\pi_{\boldsymbol{\mathrm{\theta}}}^{\mathrm{i}}\left[\,\boldsymbol{\mathrm{a}}^{\mathrm{i}}\,|\,\boldsymbol{\mathrm{s}}\,\right].$
(19)

We establish next a straightforward but useful result concerning the
stochastic policy (19).

###### Lemma 1

The stochastic policy (19) generates input samples $\boldsymbol{\mathrm{a}}$
that are feasible for the MPC scheme (9).

Because
$\Phi^{\mathrm{i}}(\boldsymbol{\mathrm{s}},\boldsymbol{\mathrm{\theta}},\boldsymbol{\mathrm{a}}^{\mathrm{i}})=+\infty$
when
$\boldsymbol{\mathrm{a}}^{\mathrm{i}}\notin\mathbb{I}(\boldsymbol{\mathrm{s}},\boldsymbol{\mathrm{\theta}})$,
policy (15) selects feasible integer inputs
$\boldsymbol{\mathrm{a}}^{\mathrm{i}}$ with probability 1. Furthermore, NLP
(17) is feasible for all
$\boldsymbol{\mathrm{a}}^{\mathrm{i}}\in\mathbb{I}(\boldsymbol{\mathrm{s}},\boldsymbol{\mathrm{\theta}})$
and all $\boldsymbol{\mathrm{d}}$, such that its solution satisfies
constraints (13b)-(13c). As a result, the samples
$\boldsymbol{\mathrm{a}}^{\mathrm{i}},\boldsymbol{\mathrm{a}}^{\mathrm{c}}$
generated from (19) are guaranteed to be feasible. $\hfill\qed$ The policy
gradient associated to (19) can be computed using (7). Unfortunately, it has
been observed that this approach is computationally expensive for continuous
input spaces (Gros and Zanon, 2020) when the policy is restricted by non-
trivial constraints. Hence, we now turn to detailing how the policy gradient
associated to policy (19) can be efficiently computed.

###  4.2 Policy gradient

Using policy (19), the stochastic policy gradient is separable between the
continuous and integer part and reads as:

$\displaystyle\nabla_{\boldsymbol{\mathrm{\theta}}}J\left(\pi_{\boldsymbol{\mathrm{\theta}}}\right)$
$\displaystyle=\mathbb{E}_{\pi_{\boldsymbol{\mathrm{\theta}}}}\left[\nabla_{\boldsymbol{\mathrm{\theta}}}\log\pi_{\boldsymbol{\mathrm{\theta}}}A_{\boldsymbol{\mathrm{\pi}}_{\boldsymbol{\mathrm{\theta}}}}\right]$
(20)
$\displaystyle=\mathbb{E}_{\pi_{\boldsymbol{\mathrm{\theta}}}}\left[\nabla_{\boldsymbol{\mathrm{\theta}}}\log\pi^{\mathrm{c}}_{\boldsymbol{\mathrm{\theta}}}A_{\boldsymbol{\mathrm{\pi}}_{\boldsymbol{\mathrm{\theta}}}}\right]+\mathbb{E}_{\pi_{\boldsymbol{\mathrm{\theta}}}}\left[\nabla_{\boldsymbol{\mathrm{\theta}}}\log\pi^{\mathrm{i}}_{\boldsymbol{\mathrm{\theta}}}A_{\boldsymbol{\mathrm{\pi}}_{\boldsymbol{\mathrm{\theta}}}}\right],$

where $A_{\boldsymbol{\mathrm{\pi}}_{\boldsymbol{\mathrm{\theta}}}}$ is the
advantage function associated to the stochastic policy (19). Using (15), we
then observe that the score function associated to the integer part of the
policy is simply given by:

$\displaystyle\nabla_{\boldsymbol{\mathrm{\theta}}}\log\pi^{\mathrm{i}}_{\boldsymbol{\mathrm{\theta}}}[\,\boldsymbol{\mathrm{a}}^{\mathrm{i}}\,|\,\boldsymbol{\mathrm{s}}\,]=$
$\displaystyle-\frac{1}{\sigma_{\mathrm{i}}}\nabla_{\boldsymbol{\mathrm{\theta}}}\Phi^{\star}_{\mathrm{i}}(\boldsymbol{\mathrm{s}},\boldsymbol{\mathrm{\theta}},\boldsymbol{\mathrm{a}}^{\mathrm{i}})$
(21)
$\displaystyle+\frac{1}{\sigma_{\mathrm{i}}}\sum_{\boldsymbol{\mathrm{i}}_{0}\in\mathbb{I}(\boldsymbol{\mathrm{s}},\boldsymbol{\mathrm{\theta}})}\pi_{\boldsymbol{\mathrm{\theta}}}^{\mathrm{i}}\left[\,\boldsymbol{\mathrm{i}}_{0}\,|\,\boldsymbol{\mathrm{s}}\,\right]\nabla_{\boldsymbol{\mathrm{\theta}}}\Phi^{\star}_{\mathrm{i}}(\boldsymbol{\mathrm{s}},\boldsymbol{\mathrm{\theta}},\boldsymbol{\mathrm{i}}_{0}).$

The computation of the policy gradient associated to the continuous part of
the stochastic policy ought to be treated differently. Indeed, it has been
observed in (Gros and Zanon, 2020) that deterministic policy gradient methods
are computationally more effective than stochastic ones for policy
approximations on problems having continuous input and state spaces. Defining
the deterministic policy for the continuous inputs
$\boldsymbol{\mathrm{a}}^{\mathrm{c}}$ as

$\displaystyle{\boldsymbol{\mathrm{\pi}}}^{\mathrm{c}}_{\boldsymbol{\mathrm{\theta}}}\left(\boldsymbol{\mathrm{s}},\boldsymbol{\mathrm{i}}\right)=\tilde{\boldsymbol{\mathrm{u}}}_{0}(\boldsymbol{\mathrm{s}},\boldsymbol{\mathrm{\theta}},\boldsymbol{\mathrm{i}},0),$
(22)

where $\tilde{\boldsymbol{\mathrm{u}}}_{0}$ is the first element of the
solution of (17), we consider the approximation (Silver et al., 2014)

$\displaystyle\mathbb{E}_{\pi_{\boldsymbol{\mathrm{\theta}}}}\left[\nabla_{\boldsymbol{\mathrm{\theta}}}\log\pi^{\mathrm{c}}_{\boldsymbol{\mathrm{\theta}}}A_{\boldsymbol{\mathrm{\pi}}_{\boldsymbol{\mathrm{\theta}}}}\right]\approx\mathbb{E}_{\pi_{\boldsymbol{\mathrm{\theta}}}}\left[\nabla_{\boldsymbol{\mathrm{\theta}}}{\boldsymbol{\mathrm{\pi}}}^{\mathrm{c}}_{\boldsymbol{\mathrm{\theta}}}\nabla_{\boldsymbol{\mathrm{a}}^{\mathrm{c}}}A_{\boldsymbol{\mathrm{\pi}}_{\boldsymbol{\mathrm{\theta}}}}\right],$
(23)

which is asymptotically exact for $\sigma_{\mathrm{c}}\rightarrow 0$ under
some technical but fairly unrestrictive assumptions. We can then use the
asymptotically exact hybrid policy gradient

$\displaystyle\widehat{\nabla_{\boldsymbol{\mathrm{\theta}}}J\left(\pi_{\boldsymbol{\mathrm{\theta}}}\right)}$
$\displaystyle=\mathbb{E}_{\pi_{\boldsymbol{\mathrm{\theta}}}}\left[\nabla_{\boldsymbol{\mathrm{\theta}}}\boldsymbol{\mathrm{\pi}}^{\mathrm{c}}_{\boldsymbol{\mathrm{\theta}}}\nabla_{\boldsymbol{\mathrm{a}}^{\mathrm{c}}}A_{\boldsymbol{\mathrm{\pi}}_{\boldsymbol{\mathrm{\theta}}}}\right]+\mathbb{E}_{\pi_{\boldsymbol{\mathrm{\theta}}}}\left[\nabla_{\boldsymbol{\mathrm{\theta}}}\log\pi^{\mathrm{i}}_{\boldsymbol{\mathrm{\theta}}}A_{\boldsymbol{\mathrm{\pi}}_{\boldsymbol{\mathrm{\theta}}}}\right],$
(24)

as a computationally effective policy gradient evaluation. The stochastic
policy (16) is then deployed on the system and generates exploration, while
the deterministic policy (22) is used to compute the policy gradient (24). We
propose next a compatible advantage function approximator for (24), offering a
systematic approximation of the advantage function
$A_{\boldsymbol{\mathrm{\pi}}_{\boldsymbol{\mathrm{\theta}}}}$.

###  4.3 Compatible advantage function approximation

We note that the advantage function approximation

$\displaystyle\hat{A}_{\pi_{\boldsymbol{\mathrm{\theta}}}}=\boldsymbol{\mathrm{w}}^{\top}\nabla_{\boldsymbol{\mathrm{\theta}}}\log\pi_{\boldsymbol{\mathrm{\theta}}}=\boldsymbol{\mathrm{w}}^{\top}\nabla_{\boldsymbol{\mathrm{\theta}}}\log\pi^{\mathrm{i}}_{\boldsymbol{\mathrm{\theta}}}+\boldsymbol{\mathrm{w}}^{\top}\nabla_{\boldsymbol{\mathrm{\theta}}}\log\pi^{\mathrm{c}}_{\boldsymbol{\mathrm{\theta}}},$
(25)

is compatible by construction (Silver et al., 2014) for the stochastic policy
gradient (20), in the sense that

$\displaystyle\nabla_{\boldsymbol{\mathrm{\theta}}}J\left(\pi_{\boldsymbol{\mathrm{\theta}}}\right)$
$\displaystyle=\mathbb{E}_{\pi_{\boldsymbol{\mathrm{\theta}}}}\left[\nabla_{\boldsymbol{\mathrm{\theta}}}\log\pi_{\boldsymbol{\mathrm{\theta}}}\hat{A}_{\boldsymbol{\mathrm{\pi}}_{\boldsymbol{\mathrm{\theta}}}}\right]$
(26)

holds if $\boldsymbol{\mathrm{w}}$ is the solution of the Least-Squares
problem

$\displaystyle\boldsymbol{\mathrm{w}}=\mathrm{arg}\min_{\boldsymbol{\mathrm{w}}}\,\frac{1}{2}\mathbb{E}_{\pi_{\boldsymbol{\mathrm{\theta}}}}\left[\left(A_{\boldsymbol{\mathrm{\pi}}_{\boldsymbol{\mathrm{\theta}}}}-\hat{A}_{\boldsymbol{\mathrm{\pi}}_{\boldsymbol{\mathrm{\theta}}}}\right)^{2}\right].$
(27)

Similarly, we seek a compatible advantage function approximation for the
hybrid policy gradient (24). We propose the hybrid advantage function
approximation, inspired from (Gros and Zanon, 2020):

$\displaystyle\hat{A}_{\pi_{\boldsymbol{\mathrm{\theta}}}}=\boldsymbol{\mathrm{w}}^{\top}\nabla_{\boldsymbol{\mathrm{\theta}}}\log\pi^{\mathrm{i}}_{\boldsymbol{\mathrm{\theta}}}+\boldsymbol{\mathrm{w}}^{\top}{\color[rgb]{0,0,0}\frac{1}{\sigma_{\mathrm{c}}}}\nabla_{\boldsymbol{\mathrm{\theta}}}\pi^{\mathrm{c}}_{\boldsymbol{\mathrm{\theta}}}M\left(\boldsymbol{\mathrm{e}}-\boldsymbol{\mathrm{c}}\right),$
(28)

where we label
$\boldsymbol{\mathrm{e}}=\boldsymbol{\mathrm{a}}^{\mathrm{c}}-{\boldsymbol{\mathrm{\pi}}}^{\mathrm{c}}_{\boldsymbol{\mathrm{\theta}}}$
the exploration performed on the continuous part of the input space
$\mathbb{R}^{m_{\mathrm{c}}}$, and $M\in\mathbb{R}^{m_{\mathrm{c}}\times
m_{\mathrm{c}}}$ is symmetric and
$\boldsymbol{\mathrm{c}}\in\mathbb{R}^{m_{\mathrm{c}}}$. We will show in the
following proposition that for $M$ and $\boldsymbol{\mathrm{c}}$ adequately
chosen, the advantage function approximation (28) is compatible with the
policy gradient (24).

###### Proposition 1

The hybrid function approximation (28) is asymptotically compatible, i.e.,

$\displaystyle\lim_{\sigma_{\mathrm{c}}\rightarrow
0}\widehat{\nabla_{\boldsymbol{\mathrm{\theta}}}J\left(\pi_{\boldsymbol{\mathrm{\theta}}}\right)}=$
$\displaystyle\,\,\lim_{\sigma_{\mathrm{c}}\rightarrow
0}\mathbb{E}_{\pi_{\boldsymbol{\mathrm{\theta}}}}\left[\nabla_{\boldsymbol{\mathrm{\theta}}}\boldsymbol{\mathrm{\pi}}^{\mathrm{c}}_{\boldsymbol{\mathrm{\theta}}}\nabla_{\boldsymbol{\mathrm{a}}^{\mathrm{c}}}\hat{A}_{\boldsymbol{\mathrm{\pi}}_{\boldsymbol{\mathrm{\theta}}}}\right]$
(29)
$\displaystyle\qquad+\mathbb{E}_{\pi_{\boldsymbol{\mathrm{\theta}}}}\left[\nabla_{\boldsymbol{\mathrm{\theta}}}\log\pi^{\mathrm{i}}_{\boldsymbol{\mathrm{\theta}}}\hat{A}_{\boldsymbol{\mathrm{\pi}}_{\boldsymbol{\mathrm{\theta}}}}\right]$

holds for $\boldsymbol{\mathrm{w}}$ solution of (27) and for
$M,\,\boldsymbol{\mathrm{c}}$ chosen according to (Gros and Zanon, 2020):

$\displaystyle\boldsymbol{\mathrm{c}}$
$\displaystyle=\frac{\sigma_{\mathrm{c}}}{2}\sum_{i=1}^{n_{\boldsymbol{\mathrm{a}}}}\frac{\partial^{2}\tilde{\boldsymbol{\mathrm{u}}}_{0}}{\partial\boldsymbol{\mathrm{d}}_{i}^{2}},\quad
M=\left(\frac{\partial\tilde{\boldsymbol{\mathrm{u}}}_{0}}{\partial\boldsymbol{\mathrm{d}}}\frac{\partial\tilde{\boldsymbol{\mathrm{u}}}_{0}}{\partial\boldsymbol{\mathrm{d}}}^{\top}\right),$
(30)

evaluated at the solution of (17) for $\boldsymbol{\mathrm{d}}=0$, where (17)
satisfies the regularity assumptions of (Gros and Zanon, 2020, Proposition 1).
These assumptions are technical but fairly unrestrictive, see (Gros and Zanon,
2020) for a complete discussion.

The proof delivered below is a sketch that follows the lines of the proof of
Proposition 1 in Gros and Zanon (2020). We observe that the solution
$\boldsymbol{\mathrm{w}}$ of (27) using (28) is given by:

$\displaystyle\mathbb{E}_{\pi_{\boldsymbol{\mathrm{\theta}}}}\left[\left(\nabla_{\boldsymbol{\mathrm{\theta}}}\log\pi^{\mathrm{i}}_{\boldsymbol{\mathrm{\theta}}}+{\color[rgb]{0,0,0}\frac{1}{\sigma_{\mathrm{c}}}}\nabla_{\boldsymbol{\mathrm{\theta}}}\pi^{\mathrm{c}}_{\boldsymbol{\mathrm{\theta}}}M\left(\boldsymbol{\mathrm{e}}-\boldsymbol{\mathrm{c}}\right)\right)\left(A_{\boldsymbol{\mathrm{\pi}}_{\boldsymbol{\mathrm{\theta}}}}-\hat{A}_{\boldsymbol{\mathrm{\pi}}_{\boldsymbol{\mathrm{\theta}}}}\right)\right]=0.$
(31)

Using a Taylor expansion of
$A_{\boldsymbol{\mathrm{\pi}}_{\boldsymbol{\mathrm{\theta}}}}$ at
$\boldsymbol{\mathrm{e}}=0$, as proposed in (Gros and Zanon, 2020, Proposition
1), we observe that (31) becomes:

$\displaystyle\mathbb{E}_{\pi_{\boldsymbol{\mathrm{\theta}}}}\left[\nabla_{\boldsymbol{\mathrm{\theta}}}\log\pi^{\mathrm{i}}_{\boldsymbol{\mathrm{\theta}}}\left(A_{\boldsymbol{\mathrm{\pi}}_{\boldsymbol{\mathrm{\theta}}}}-\hat{A}_{\boldsymbol{\mathrm{\pi}}_{\boldsymbol{\mathrm{\theta}}}}\right)\right]+{\color[rgb]{0,0,0}\frac{1}{\sigma_{\mathrm{c}}}}\mathbb{E}_{\pi_{\boldsymbol{\mathrm{\theta}}}}\left[\nabla_{\boldsymbol{\mathrm{\theta}}}\pi^{\mathrm{c}}_{\boldsymbol{\mathrm{\theta}}}M\left(\boldsymbol{\mathrm{e}}-\boldsymbol{\mathrm{c}}\right)\xi\right]$
$\displaystyle+\mathbb{E}_{\pi_{\boldsymbol{\mathrm{\theta}}}}\left[\nabla_{\boldsymbol{\mathrm{\theta}}}\pi^{\mathrm{c}}_{\boldsymbol{\mathrm{\theta}}}M{\color[rgb]{0,0,0}\frac{\boldsymbol{\mathrm{e}}-\boldsymbol{\mathrm{c}}}{\sigma_{\mathrm{c}}}}\left(A_{\boldsymbol{\mathrm{\pi}}_{\boldsymbol{\mathrm{\theta}}}}-\hat{A}_{\boldsymbol{\mathrm{\pi}}_{\boldsymbol{\mathrm{\theta}}}}\right)\right]+$
(32)
$\displaystyle+\mathbb{E}_{\pi_{\boldsymbol{\mathrm{\theta}}}}\left[\nabla_{\boldsymbol{\mathrm{\theta}}}\pi^{\mathrm{c}}_{\boldsymbol{\mathrm{\theta}}}M{\color[rgb]{0,0,0}\frac{\left(\boldsymbol{\mathrm{e}}-\boldsymbol{\mathrm{c}}\right)\boldsymbol{\mathrm{e}}^{\top}}{\sigma_{\mathrm{c}}}}\left(\nabla_{\boldsymbol{\mathrm{a}}^{\mathrm{c}}}A_{\boldsymbol{\mathrm{\pi}}_{\boldsymbol{\mathrm{\theta}}}}-\nabla_{\boldsymbol{\mathrm{a}}^{\mathrm{c}}}\hat{A}_{\boldsymbol{\mathrm{\pi}}_{\boldsymbol{\mathrm{\theta}}}}\right)\right]=0,$

where $\xi$ is the second-order remainder of the Taylor expansion of
$A_{\boldsymbol{\mathrm{\pi}}_{\boldsymbol{\mathrm{\theta}}}}$. Unlike (31),
all terms in (4.3) are evaluated at
$\boldsymbol{\mathrm{s}},\,\boldsymbol{\mathrm{a}}^{\mathrm{c}}=\boldsymbol{\mathrm{\pi}}^{\mathrm{c}}\left(\boldsymbol{\mathrm{s}}\right)$.
Following a similar argumentation as in (Gros and Zanon, 2020, Proposition 1),
we obtain

$\displaystyle\lim_{\sigma_{\mathrm{c}}\rightarrow 0}$
$\displaystyle\mathbb{E}_{\pi_{\boldsymbol{\mathrm{\theta}}}}\left[{\color[rgb]{0,0,0}\frac{1}{\sigma_{\mathrm{c}}}}\nabla_{\boldsymbol{\mathrm{\theta}}}\pi^{\mathrm{c}}_{\boldsymbol{\mathrm{\theta}}}M\left(\boldsymbol{\mathrm{e}}-\boldsymbol{\mathrm{c}}\right)\boldsymbol{\mathrm{e}}^{\top}\left(\nabla_{\boldsymbol{\mathrm{a}}^{\mathrm{c}}}A_{\boldsymbol{\mathrm{\pi}}_{\boldsymbol{\mathrm{\theta}}}}-\nabla_{\boldsymbol{\mathrm{a}}^{\mathrm{c}}}\hat{A}_{\boldsymbol{\mathrm{\pi}}_{\boldsymbol{\mathrm{\theta}}}}\right)\right]$
$\displaystyle=\lim_{\sigma_{\mathrm{c}}\rightarrow
0}\mathbb{E}_{\pi_{\boldsymbol{\mathrm{\theta}}}}\left[\nabla_{\boldsymbol{\mathrm{\theta}}}\pi^{\mathrm{c}}_{\boldsymbol{\mathrm{\theta}}}\left(\nabla_{\boldsymbol{\mathrm{a}}^{\mathrm{c}}}A_{\boldsymbol{\mathrm{\pi}}_{\boldsymbol{\mathrm{\theta}}}}-\nabla_{\boldsymbol{\mathrm{a}}^{\mathrm{c}}}\hat{A}_{\boldsymbol{\mathrm{\pi}}_{\boldsymbol{\mathrm{\theta}}}}\right)\right],$
(33a) $\displaystyle\lim_{\sigma_{\mathrm{c}}\rightarrow 0}$
$\displaystyle\mathbb{E}_{\pi_{\boldsymbol{\mathrm{\theta}}}}\left[{\color[rgb]{0,0,0}\frac{1}{\sigma_{\mathrm{c}}}}\nabla_{\boldsymbol{\mathrm{\theta}}}\pi^{\mathrm{c}}_{\boldsymbol{\mathrm{\theta}}}M\left(\boldsymbol{\mathrm{e}}-\boldsymbol{\mathrm{c}}\right)\xi\right]=0,$
(33b) $\displaystyle\lim_{\sigma_{\mathrm{c}}\rightarrow 0}$
$\displaystyle\mathbb{E}_{\pi_{\boldsymbol{\mathrm{\theta}}}}\left[\nabla_{\boldsymbol{\mathrm{\theta}}}\pi^{\mathrm{c}}_{\boldsymbol{\mathrm{\theta}}}M{\color[rgb]{0,0,0}\frac{\boldsymbol{\mathrm{e}}-\boldsymbol{\mathrm{c}}}{\sigma_{\mathrm{c}}}}\left(A_{\boldsymbol{\mathrm{\pi}}_{\boldsymbol{\mathrm{\theta}}}}-\hat{A}_{\boldsymbol{\mathrm{\pi}}_{\boldsymbol{\mathrm{\theta}}}}\right)_{\boldsymbol{\mathrm{e}}=0}\right]=0.$
(33c)

Equality (33b) holds from the Delta method, while equalities (33a), (33c) hold
because

$\displaystyle\lim_{\sigma_{\mathrm{c}}\rightarrow
0}\mathbb{E}_{\pi_{\boldsymbol{\mathrm{\theta}}}}\left[\frac{1}{\sigma_{\mathrm{c}}}M\left(\boldsymbol{\mathrm{e}}-\boldsymbol{\mathrm{c}}\right)\boldsymbol{\mathrm{e}}^{\top}\right]=I,$
(34) $\displaystyle\lim_{\sigma_{\mathrm{c}}\rightarrow
0}\mathbb{E}_{\pi_{\boldsymbol{\mathrm{\theta}}}}\left[M{\color[rgb]{0,0,0}\frac{\boldsymbol{\mathrm{e}}-\boldsymbol{\mathrm{c}}}{\sigma_{\mathrm{c}}}}\right]=0,$
(35)

result from (30), see (Gros and Zanon, 2020). Hence

$\displaystyle\lim_{\sigma_{\mathrm{c}}\rightarrow
0}\mathbb{E}_{\pi_{\boldsymbol{\mathrm{\theta}}}}\left[\nabla_{\boldsymbol{\mathrm{\theta}}}\log\pi^{\mathrm{i}}_{\boldsymbol{\mathrm{\theta}}}\left(A_{\boldsymbol{\mathrm{\pi}}_{\boldsymbol{\mathrm{\theta}}}}-\hat{A}_{\boldsymbol{\mathrm{\pi}}_{\boldsymbol{\mathrm{\theta}}}}\right)\right]$
(36) $\displaystyle\hskip
42.67912pt+\mathbb{E}_{\pi_{\boldsymbol{\mathrm{\theta}}}}\left[\nabla_{\boldsymbol{\mathrm{\theta}}}\pi^{\mathrm{c}}_{\boldsymbol{\mathrm{\theta}}}\left(\nabla_{\boldsymbol{\mathrm{a}}^{\mathrm{c}}}A_{\boldsymbol{\mathrm{\pi}}_{\boldsymbol{\mathrm{\theta}}}}-\nabla_{\boldsymbol{\mathrm{a}}^{\mathrm{c}}}\hat{A}_{\boldsymbol{\mathrm{\pi}}_{\boldsymbol{\mathrm{\theta}}}}\right)\right]=0.$

Using (24), (29) holds from (36). $\hfill\qed$

##  5 NLP sensitivities

In order to deploy the policy gradient techniques described above, one needs
to compute the sensitivities
$\nabla_{\boldsymbol{\mathrm{\theta}}}\boldsymbol{\mathrm{\pi}}^{\mathrm{c}}_{\boldsymbol{\mathrm{\theta}}}$
and
$\nabla_{\boldsymbol{\mathrm{\theta}}}\log\pi^{\mathrm{i}}_{\boldsymbol{\mathrm{\theta}}}$.
Computing the score function (21) requires computing the sensitivity of the
cost function $\Phi^{\star}_{\mathrm{i}}$ of the NLP (14). This sensitivity
exists almost everywhere and is given by:

$\displaystyle\nabla_{\boldsymbol{\mathrm{\theta}}}\Phi^{\star}_{\mathrm{i}}(\boldsymbol{\mathrm{s}},\boldsymbol{\mathrm{\theta}},\boldsymbol{\mathrm{a}}^{\mathrm{i}})=\nabla_{\boldsymbol{\mathrm{\theta}}}\mathcal{L}(\boldsymbol{\mathrm{y}},\boldsymbol{\mathrm{\lambda}},\boldsymbol{\mathrm{\mu}},\boldsymbol{\mathrm{s}},\boldsymbol{\mathrm{\theta}},\boldsymbol{\mathrm{a}}^{\mathrm{i}},\boldsymbol{\mathrm{d}}),$
(37)

where $\boldsymbol{\mathrm{y}}$ is the primal solution of the NLP (14),
gathering the continuous inputs and states of the NLP, and
$\boldsymbol{\mathrm{\lambda}},\boldsymbol{\mathrm{\mu}}$ the dual variables
associated to constraints (13b)-(13c), respectively, and
$\mathcal{L}=\Phi+\boldsymbol{\mathrm{d}}^{\top}\boldsymbol{\mathrm{u}}_{0}+\boldsymbol{\mathrm{\lambda}}^{\top}\boldsymbol{\mathrm{f}}+\boldsymbol{\mathrm{\mu}}^{\top}\boldsymbol{\mathrm{h}}$
is the Lagrange function associated to (14). The computation of
$\nabla_{\boldsymbol{\mathrm{\theta}}}\boldsymbol{\mathrm{\pi}}^{\mathrm{c}}_{\boldsymbol{\mathrm{\theta}}}$
is more involved. Consider:

$\displaystyle\boldsymbol{\mathrm{r}}=\left[\begin{array}[]{c}\nabla_{\boldsymbol{\mathrm{y}}}\mathcal{L}(\boldsymbol{\mathrm{z}},\boldsymbol{\mathrm{s}},\boldsymbol{\mathrm{\theta}},\boldsymbol{\mathrm{a}}^{\mathrm{i}},\boldsymbol{\mathrm{d}})\\\
\boldsymbol{\mathrm{f}}\left(\boldsymbol{\mathrm{w}},\boldsymbol{\mathrm{s}},\boldsymbol{\mathrm{\theta}}\right)\\\
\mathrm{diag}(\boldsymbol{\mathrm{\mu}})\boldsymbol{\mathrm{h}}\left(\boldsymbol{\mathrm{w}},\boldsymbol{\mathrm{\theta}}\right)+\tau\end{array}\right]=0,$
(41)

i.e., the primal-dual interior-point KKT conditions associated to (14) for a
barrier parameter $\tau>0$, and $\boldsymbol{\mathrm{z}}$ gathering the
primal-dual variables of the NLP (14), i.e.,
$\boldsymbol{\mathrm{z}}=\left\\{\boldsymbol{\mathrm{y}},\boldsymbol{\mathrm{\lambda}},\boldsymbol{\mathrm{\mu}}\right\\}$.
Then, if the solution of the NLP (14) satisfies LICQ and SOSC (Nocedal and
Wright, 2006), the sensitivity of the solution of the NLP (14) exists almost
everywhere and can be computed via the Implicit Function Theorem, providing

$\displaystyle\frac{\partial\boldsymbol{\mathrm{z}}}{\partial\boldsymbol{\mathrm{\theta}}}=-\frac{\partial\boldsymbol{\mathrm{r}}}{\partial\boldsymbol{\mathrm{z}}}^{-1}\frac{\partial\boldsymbol{\mathrm{r}}}{\partial\boldsymbol{\mathrm{\theta}}},$
(42)

see (Büskens and Maurer, 2001). Using (22), the sensitivity
$\nabla_{\boldsymbol{\mathrm{\theta}}}\boldsymbol{\mathrm{\pi}}^{\mathrm{c}}_{\boldsymbol{\mathrm{\theta}}}$
then read as

$\displaystyle\nabla_{\boldsymbol{\mathrm{\theta}}}\boldsymbol{\mathrm{\pi}}^{\mathrm{c}}_{\boldsymbol{\mathrm{\theta}}}=\nabla_{\boldsymbol{\mathrm{\theta}}}\tilde{\boldsymbol{\mathrm{u}}}_{0},$
(43)

where
$\nabla_{\boldsymbol{\mathrm{\theta}}}\tilde{\boldsymbol{\mathrm{u}}}_{0}$ is
extracted from
$\frac{\partial\boldsymbol{\mathrm{z}}}{\partial\boldsymbol{\mathrm{\theta}}}$.

##  6 Simulated example

For the sake of brevity and in order to present results that are easy to
interpret and verify, we propose to use a very low dimensional example,
allowing us to bypass the evaluation of the action-value function via
Temporal-Difference techniques, and isolate the discussions of this paper from
questions regarding TD methods. We consider the linear, scalar dynamics:

$\displaystyle s_{k+1}=s_{k}+a^{\mathrm{c}}_{k}\,i_{k}+n_{k}$ (44)

where $s_{k},a^{\mathrm{c}}_{k}\in\mathbb{R}$,
$i_{k}\in\left\\{0,\,1\right\\}$ and $n_{k}$ is uniformly distributed in
$[0,\,0.05]$. We consider the baseline stage cost:

$\displaystyle L\left(s,\boldsymbol{\mathrm{a}}\right)=$
$\displaystyle\frac{1}{2}(s-s_{\mathrm{ref}})^{2}+\frac{1}{2}(a^{\mathrm{c}}-a^{\mathrm{c}}_{\mathrm{ref}})^{2}+w\,i$
(45) $\displaystyle\qquad+c\,\max\left(|s|-0.2,0\right),$

as the reference performance, where $w,\,c\in\mathbb{R}_{+}$ are scalar weight
and $s_{\mathrm{ref}},\,a_{\mathrm{ref}}$ are references for the state and
continuous input. The MPC model is deterministic, given by:

$\displaystyle x_{k+1}=x_{k}+u_{k}\,i_{k}+b$ (46)

where $b\in\mathbb{R}$ is constant, but subject to adaptation via RL. The
baseline cost imposes a high penalty for $s\notin[-0.2,0.2]$, and constitutes
an exact relaxation of the constraint $-0.2\leq s\leq 0.2$, see (Gros and
Zanon, 2019). The MPC stage cost $\ell$ has the form (45). The MPC parameters
$s_{\mathrm{ref}},\,a^{\mathrm{c}}_{\mathrm{ref}}$, $c$ and $b$ are subject to
adaptation via RL.

The policy gradient (29) was implemented, where the advantage function
estimation was computed from (27), using the approximator
$\hat{A}_{\pi_{\boldsymbol{\mathrm{\theta}}}}$ from (28). The true advantage
function $A_{\boldsymbol{\mathrm{\pi}}_{\boldsymbol{\mathrm{\theta}}}}$ was
evaluated via classic policy evaluation (Sutton and Barto, 1998) in order to
deliver unambiguous results. On more complex examples (27) would be evaluated
via Temporal-Difference techniques. The evaluations of (29) and (27) were
performed in a batch fashion, using 30 batches of 50 time steps each, all
starting from the deterministic initial condition
$\boldsymbol{\mathrm{s}}_{0}=0$. The MPC scheme had a horizon of $N=10$ time
samples, and a terminal cost based on the Riccati matrix of the control
problem with $i=1$. A discount factor of $\gamma=0.95$ was adopted. The step-
size selected for adapting the parameters from the policy gradient was
$\alpha=2\cdot 10^{-3}$. The exploration parameters were chosen as
$\sigma_{\mathrm{i}}=2\cdot 10^{-2}$, $\sigma_{\mathrm{c}}=10^{-2}$.

The parameters $s_{\mathrm{ref}}=a^{\mathrm{c}}_{\mathrm{ref}}=0$, $w=0.2$,
$c=1$ were adopted for the baseline cost. The MPC scheme parameters were
initialized using the same values, and using $b=0$. Fig. 1 reports the
trajectories of the system at the beginning and end of the learning process,
showing how performance is gained by bringing the state trajectories in the
interval $[-0.2,0.2]$. Fig. 2 reports the policy for the continuous and
integer inputs, showing how RL reshapes the MPC policy for a better closed-
loop performance. Fig. 3 reports the estimated policy gradients via the
compatible approximation (29) and directly via (24), showing a match predicted
by Prop. 1. Fig. 4 reports the closed-loop performance of the MPC controller,
calculated from
$J\left(\boldsymbol{\mathrm{\pi}}_{\boldsymbol{\mathrm{\theta}}}\right)=V_{\boldsymbol{\mathrm{\pi}}_{\boldsymbol{\mathrm{\theta}}}}(\boldsymbol{\mathrm{s}}_{0})$,
and shows the performance gain obtained via the learning. Fig. 5 shows the MPC
parameter evolution through the learning process.

Figure 1: Closed-loop trajectories before (light grey) and after (black) the
learning. The left graph shows the extreme values of the state trajectories,
the middle graph shows the extreme values of the continuous input
$a^{\mathrm{c}}_{k}$ when $i_{k}=1$, and the right graph shows the proportion
of $i_{k}=1$. Figure 2: Policy before (light grey) and after (black) the
learning. The left graph shows the Softmax policy (15) as a function of the
state $s$, giving the probability of selecting $i=1$ for a given state $s$.
The right graph shows the MPC policy (without the stochastic choice of integer
variable). Figure 3: Policy gradients throughout the learning process
(iterations of the RL method). The dots display the policy gradient as
obtained from (29), while the circles display the policy gradients obtained
from (24). Figure 4: Evolution of the closed-loop relative performance
throughout the learning process. A reduction of the cost of over 15% is
achieved here from $\boldsymbol{\mathrm{\theta}}_{0}$. Figure 5: Evolution of
the MPC parameters throughout the learning process. The references are
adjusted so that the system trajectories are better contained in the interval
$[-0.2,0.2]$. The model bias $b$ does not match the value that a classic
Prediction Error Method would deliver ($b=\mathbb{E}[n_{k}]=0.025$, dashed
line), while the cost associated to constraints is left unchanged.

##  7 Discussion & Conclusion

This paper proposed an actor-critic approach to compute the policy gradient
associated to policy approximations based on mixed-integer MPC schemes. The
methodology is generic and applicable to linear, nonlinear and robust
approaches. The paper proposes a hybrid stochastic-deterministic policy
approach to generate the exploration and evaluate the policy gradient,
avoiding the heavy computational expenses associated to using a stochastic
policy approach on problems having continuous inputs and state constraints. A
simple, compatible advantage function approximation is then proposed, tailored
to our formulation and to MPC-based policy approximations. Some implementation
details are provided, and the methods are illustrated on a simple example,
providing a clear picture of how the proposed method is performing.

Future work will consider extensions to reduce the noise in the policy
gradient estimation resulting from the choice of advantage function
approximation, and will investigate techniques to integrate the stochastic
policy and sensitivity computations with the branch-and-bound techniques used
to solve the mixed-integer MPC problem. Future work will also investigate the
potential of using the approaches detailed here to offer computationally less
expensive approaches to solve the mixed-integer problem.

## References and Notes

  * Abbeel et al. (2007) Abbeel, P., Coates, A., Quigley, M., and Ng, A.Y. (2007).  An application of reinforcement learning to aerobatic helicopter flight.  In _In Advances in Neural Information Processing Systems 19_ , 2007\. MIT Press.
  * Bertsekas (2007) Bertsekas, D. (2007).  _Dynamic Programming and Optimal Control_ , volume 2.  Athena Scientific, 3rd edition.
  * Bertsekas (1995) Bertsekas, D. (1995).  _Dynamic Programming and Optimal Control_ , volume 1 and 2\.  Athena Scientific, Belmont, MA.
  * Bertsekas and Shreve (1996) Bertsekas, D. and Shreve, S. (1996).  _Stochastic Optimal Control: The Discrete Time Case_.  Athena Scientific, Belmont, MA.
  * Büskens and Maurer (2001) Büskens, C. and Maurer, H. (2001).  _Online Optimization of Large Scale Systems_ , chapter Sensitivity Analysis and Real-Time Optimization of Parametric Nonlinear Programming Problems, 3–16.  Springer Berlin Heidelberg, Berlin, Heidelberg.
  * Gros and Zanon (2019) Gros, S. and Zanon, M. (2019).  Data-Driven Economic NMPC using Reinforcement Learning.  _IEEE Transactions on Automatic Control_.  (in press).
  * Gros and Zanon (2020) Gros, S. and Zanon, M. (2020).  Safe Reinforcement Learning Based on Robust MPC and Policy Gradient Methods.  _IEEE Transactions on Automatic Control (submitted)_.
  * J. Garcia (2013) J. Garcia, J.F. (2013).  A comprehensive survey on safe reinforcement learning.  _Journal of Machine Learning Research_ , 16, 1437–1480.
  * Nocedal and Wright (2006) Nocedal, J. and Wright, S. (2006).  _Numerical Optimization_.  Springer Series in Operations Research and Financial Engineering. Springer, 2 edition.
  * Silver et al. (2014) Silver, D., Lever, G., Heess, N., Degris, T., Wierstra, D., and Riedmiller, M. (2014).  Deterministic policy gradient algorithms.  In _Proceedings of the 31st International Conference on Machine Learning_ , ICML’14, I–387–I–395.
  * Sutton and Barto (1998) Sutton, R.S. and Barto, A.G. (1998).  _Introduction to Reinforcement Learning_.  MIT Press, Cambridge, MA, USA, 1st edition.
  * Sutton et al. (1999) Sutton, R.S., McAllester, D., Singh, S., and Mansour, Y. (1999).  Policy gradient methods for reinforcement learning with function approximation.  In _Proceedings of the 12th International Conference on Neural Information Processing Systems_ , NIPS’99, 1057–1063. MIT Press, Cambridge, MA, USA.
  * Wang et al. (2012) Wang, S., Chaovalitwongse, W., and Babuska, R. (2012).  Machine learning algorithms in bipedal robot control.  _Trans. Sys. Man Cyber Part C_ , 42(5), 728–743.
  * Zanon and Gros (2019) Zanon, M. and Gros (2019).  Safe Reinforcement Learning Using Robust MPC.  In _Transaction on Automatic Control, (submitted))_.  Https://arxiv.org/abs/1906.04005.