#  Boundary metrics on soliton moduli spaces

Paul Sutcliffe
Department of Mathematical Sciences,
Durham University, Durham DH1 3LE, United Kingdom.
Email<EMAIL_ADDRESS>
(October 2021)

###### Abstract

The geodesic approximation is a powerful method for studying the dynamics of
BPS solitons. However, there are systems, such as BPS monopoles in three-
dimensional hyperbolic space, where this approach is not applicable because
the moduli space metric defined by the kinetic energy is not finite. In the
case of hyperbolic monopoles, an alternative metric has been defined using the
abelian connection on the sphere at infinity, but its relation to the dynamics
of hyperbolic monopoles is unclear. Here this metric is placed in a more
general context of boundary metrics on soliton moduli spaces. Examples are
studied in systems in one and two space dimensions, where it is much easier to
compare the results with simulations of the full nonlinear field theory
dynamics. It is found that geodesics of the boundary metric provide a
reasonable description of soliton dynamics.

##  1 Introduction

In the geodesic approximation the motion of slowly moving BPS solitons is
described by geodesic motion in the moduli space of static solitons, equipped
with a metric induced by the kinetic energy of the field theory [1]. This
approach has been applied to a wide variety of soliton systems, following its
initial application to the study of magnetic monopole dynamics in three-
dimensional Euclidean space [2]. The physical intuition underlying the success
of this approximation is that slowly moving solitons possess only a small
amount of kinetic energy and therefore energy conservation prevents the field
configuration from deforming far from the moduli space of minimal energy
static soliton solutions. Projecting the dynamics to an adiabatic evolution in
the moduli space therefore appears a reasonable simplification. This physical
argument is supported by rigorous mathematical justification in the case of
vortices and monopoles [3, 4].

There are some BPS soliton systems where the geodesic approximation is not
applicable because restricting the dynamics to motion in the moduli space
yields infinite kinetic energy, and therefore is not a prescription to obtain
a metric on the soliton moduli space. A notable example is the case of BPS
monopoles in three-dimensional hyperbolic space. The static properties of
monopoles in hyperbolic space share many of the features of monopoles in
Euclidean space [5], but this crucial difference in the divergence of the
metric means that nothing is known about the dynamics of hyperbolic monopoles.
Braam and Austin defined a metric on the moduli space of hyperbolic monopoles
using the abelian connection on the sphere at infinity [6], but its relation
to the dynamics of hyperbolic monopoles remains unclear. The moduli space of a
single hyperbolic monopole equipped with the Braam-Austin metric is three-
dimensional hyperbolic space, which is encouraging, as interpreting its
geodesics as the free motion of a single hyperbolic monopole seems reasonable.
Simulating the full field dynamics of hyperbolic monopoles would be a worthy
endeavour, although this is by no means an easy task.

In this paper an alternative approach is adopted to gain some support for the
relevance of the Braam-Austin metric to the dynamics of hyperbolic monopoles.
This metric is placed in a more general context of boundary metrics on BPS
soliton moduli spaces. The metric is obtained by a simple renormalization of
the divergent metric induced from the field theory kinetic energy. By
introducing a cutoff, the leading divergence can be identified and used to
renormalize the metric via multiplication by a factor that depends only on the
cutoff. Multiplication by this factor, of course, does not change the geodesic
equation, so the basic fact being exploited by this approach is that there is
a convergence for the geodesic equation as the cutoff is removed, even though
the natural metric diverges. The geodesic equation obtained in the limit is
the equation for the geodesics of the renormalized metric, which is termed a
boundary metric because it depends only the soliton fields at spatial
infinity. The next section provides the details of the construction of the
boundary metric for hyperbolic monopoles, and further sections consider a
radial lump in the $\mathbb{CP}^{1}$ sigma model in the plane, and kink
dynamics on a line. The advantage of low-dimensional examples is that
comparisons with full nonlinear field theory simulations are much easier to
investigate.

##  2 Hyperbolic monopoles

This section concerns $SU(2)$ magnetic monopoles in three-dimensional
hyperbolic space of curvature $-1$, with metric

$ds^{2}({\mathbb{H}}^{3})=\frac{4(dx_{1}^{2}+dx_{2}^{2}+dx_{3}^{2})}{(1-r^{2})^{2}},$
(2.1)

where $r^{2}=x_{1}^{2}+x_{2}^{2}+x_{3}^{2}<1.$ The Higgs field $\Phi$ and the
gauge potential $A_{i}$ take values in $\mathfrak{su}(2)$, and the monopole
scale relative to the curvature of hyperbolic space is fixed by the boundary
condition $|\Phi|^{2}=-\frac{1}{2}\mbox{tr}(\Phi^{2})\to\frac{1}{4}$ as $r\to
1.$

Hyperbolic monopoles are solutions of the Bogomolny equation

$F_{ij}=\frac{2}{1-r^{2}}\varepsilon_{ijk}D_{k}\Phi,$ (2.2)

where $\varepsilon_{ijk}$ is the totally anti-symmetric tensor,
$F_{ij}=\partial_{i}A_{j}-\partial_{j}A_{i}+[A_{i},A_{j}]$ is the gauge field,
and $D_{i}\Phi=\partial_{i}\Phi+[A_{i},\Phi]$ is the covariant derivative of
the Higgs field.

The Higgs field at infinity defines a map between two-spheres,
$\Phi|_{r=1}:S^{2}\to S^{2}$, and the monopole charge $N$ is the associated
positive integer element of the homotopy group $\pi_{2}(S^{2})={\mathbb{Z}}.$
The moduli space, ${\cal M}_{N}$, of charge $N$ monopole solutions of the
Bogomolny equation (2.2), up to gauge transformations, is a
$(4N-1)$-dimensional manifold [5]. In regions of this moduli space where the
monopoles are well-separated, these parameters may be interpreted as a
position for each monopole and a set of relative phases.

Monopole solutions of the Bogomolny equation are static solutions of a
dynamical $SU(2)$ Yang-Mills-Higgs theory with metric
$ds^{2}=dt^{2}-ds^{2}({\mathbb{H}}^{3})$ and kinetic energy

$T=\int\limits_{r\leq
1}-\mbox{tr}\bigg{\\{}\frac{2}{(1-r^{2})^{3}}(D_{t}\Phi)^{2}+\frac{1}{2(1-r^{2})}E_{i}^{2}\bigg{\\}}dx_{1}dx_{2}dx_{3},$
(2.3)

where $D_{t}\Phi=\dot{\Phi}+[A_{t},\Phi]$ and
$E_{i}=\dot{A}_{i}-\partial_{i}A_{t}+[A_{t},A_{i}],$ with a dot denoting
differentiation with respect to $t.$ The equation of motion for $A_{t}$ that
follows from the variation of (2.3) is Gauss’ law

$D_{i}\bigg{(}\frac{E_{i}}{1-r^{2}}\bigg{)}=\frac{4}{(1-r^{2})^{3}}[D_{t}\Phi,\Phi].$
(2.4)

In the gauge $A_{t}=0$ the kinetic energy simplifies to

$T=\int\limits_{r\leq
1}-\mbox{tr}\bigg{\\{}\frac{2}{(1-r^{2})^{3}}\dot{\Phi}^{2}+\frac{1}{2(1-r^{2})}\dot{A}_{i}^{2}\bigg{\\}}d^{3}x,$
(2.5)

subject to the constraint from Gauss’ law

$D_{i}\bigg{(}\frac{\dot{A}_{i}}{1-r^{2}}\bigg{)}=\frac{4}{(1-r^{2})^{3}}[\dot{\Phi},\Phi],$
(2.6)

to ensure that the motion is orthogonal to the gauge orbit.

In the geodesic approximation the time evolution of the fields is restricted
to time evolution in the moduli space. Taking the time derivative of (2.2)
yields

$D_{i}\dot{A}_{j}-D_{j}\dot{A}_{i}=\frac{2}{1-r^{2}}\varepsilon_{ijk}(D_{k}\dot{\Phi}+[\dot{A}_{k},\Phi]).$
(2.7)

The solutions $(\dot{\Phi},\dot{A}_{i})$ of (2.7) and (2.6) are the tangent
vectors, but each has a length-squared, given by the kinetic energy (2.5),
that is infinite. This is the reason that the geodesic approximation is not
applicable to hyperbolic monopoles. To examine the details of this divergence
let $T_{b}$ denote the regularized kinetic energy evaluated with a cutoff, so
that the region of integration is $r\leq b$, with $T_{b}$ finite for $0<b<1.$
A naive examination of (2.5) might conclude that the first term leads to a
divergence as $b\to 1,$ due to the appearance of the factor $(1-r^{2})^{-3}.$
However, the condition of orthogonality to gauge orbits yields a finite limit
for $\dot{\Phi}(1-r^{2})^{-2}$ as $r\to 1$, so there is no divergence of the
kinetic energy as $b\to 1$ from this first term. The divergence arises solely
from the contribution of $\dot{A}_{i}$ in the direction of $\Phi,$ because the
abelian component $\dot{a}_{i}=-\mbox{tr}(\dot{A}_{i}\Phi)$ is finite as $r\to
1.$ This produces a logarithmic divergence with $-T_{b}/\log(1-b)$ finite as
$b\to 1.$

The boundary metric is defined by the associated renormalized kinetic energy

$\widehat{T}=-\lim_{b\to
1}T_{b}/\log(1-b)=\int\limits_{S^{2}}\frac{1}{2}\dot{a}_{i}^{2}\sin\theta\,d\theta
d\varphi,$ (2.8)

where the integration is over the boundary sphere $r=1$, with $\theta$ and
$\varphi$ the usual angular coordinates on $S^{2}.$ On this boundary sphere
Gauss’ law (2.6) projected onto the abelian component reduces to $a_{r}=0$ and
the requirement that the angular components must satisfy

$\partial_{\theta}(\dot{a}_{\theta}\sin\theta)+\frac{\partial_{\varphi}\dot{a}_{\varphi}}{\sin\theta}=0.$
(2.9)

Finally, with this condition, the renormalized kinetic energy on the monopole
moduli space becomes that of an abelian gauge theory on the sphere

$\widehat{T}=\int\limits_{S^{2}}\frac{1}{2}\bigg{(}\dot{a}_{\theta}^{2}+\frac{\dot{a}_{\varphi}^{2}}{\sin^{2}\theta}\bigg{)}\sin\theta\,d\theta
d\varphi.$ (2.10)

The fact that this defines a metric on the monopole moduli space relies on the
result that, up to gauge transformations, the abelian connection
$a=a_{\theta}d\theta+a_{\varphi}d\varphi$ on $S^{2}$ determines the monopole
[6].

The boundary metric defined by the renormalized kinetic energy (2.10) is the
metric introduced by Braam and Austin [6]. They used a symmetry argument to
show that the single monopole moduli space, ${\cal M}_{1}$, equipped with this
metric, is three-dimensional hyperbolic space ${\mathbb{H}}^{3}$. More
explicitly, let ${\bf X}$ be the point in the interior of the unit ball that
corresponds to the position of the monopole, given by the location at which
$\Phi$ vanishes. In a suitable gauge

$\dot{a}_{\theta}=-\mbox{tr}\bigg{(}\Phi\dot{\bf X}\cdot\frac{\partial
A_{\theta}}{\partial{\bf
X}}\bigg{)}\bigg{|}_{r=1},\qquad\dot{a}_{\varphi}=-\mbox{tr}\bigg{(}\Phi\dot{\bf
X}\cdot\frac{\partial A_{\varphi}}{\partial{\bf X}}\bigg{)}\bigg{|}_{r=1},$
(2.11)

satisfy (2.9) and then (2.10) yields

$\widehat{T}=\frac{4\pi}{3}\frac{|\dot{\bf X}|^{2}}{(1-|{\bf X}|^{2})^{2}},$
(2.12)

which is indeed proportional to $ds^{2}({\mathbb{H}}^{3}).$

To even consider the relevance to hyperbolic monopole dynamics of the
geodesics of the boundary metric given by (2.10), it is necessary to clarify
the type of field theory initial conditions that might be appropriate.
Clearly, the fields $(\Phi,A_{i})$ at an initial time should be taken to be
those of a point in the moduli space ${\cal M}_{N}$, but the initial time
derivatives cannot be taken to be a tangent vector $(\dot{\Phi},\dot{A}_{i})$,
as this requires infinite kinetic energy and is therefore unphysical. The
suggestion here is to curtail the tangent vector by multiplying it by a
smoothed indicator function $\chi$, satisfying $\chi\approx 1$ in the core of
the monopoles and $\chi\approx 0$ far from the core of any monopole, with
$\chi\to 0$ as $r\to 1$ sufficiently rapidly so that the kinetic energy is
finite. A suitable choice is $\chi=1-(2|\Phi|)^{m}$, where $m$ is a large
positive integer. For any initial condition of this form, causality
considerations imply that for a sufficiently large $m$, the evolution within
any given spatial region will be insensitive to the value of $m$ for evolution
up to some time limit that can be increased by increasing $m$.

##  3 Radial lumps in the $\mathbb{CP}^{1}$ sigma model

For a second example of a boundary metric, consider the $\mathbb{CP}^{1}$
sigma model, in three-dimensional Minkowski spacetime, given by the Lagrangian
density

${\cal L}=\frac{\partial_{\mu}W\partial^{\mu}\overline{W}}{(1+|W|^{2})^{2}}.$
(3.1)

Here $W$ takes values in the Riemann sphere and satisfies the boundary
condition $W\to 0$ as $x_{1}^{2}+x_{2}^{2}\to\infty.$ This one-point
compactification of the plane implies that at any fixed time
$W:{\mathbb{R}}^{2}\cup\\{\infty\\}\mapsto\mathbb{CP}^{1}$, with an associated
topological charge $N\in{\mathbb{Z}}=\pi_{2}(\mathbb{CP}^{1})$ that counts the
number of solitons, often known as lumps in this context.

For each positive integer $N$, the general static $N$-lump solution is given
by taking $W$ to be a rational function of $x_{1}+ix_{2}$ with degree $N$.
Restricting attention to static radially symmetric charge $N$ lumps positioned
at the origin gives a two-dimensional moduli space of solutions

$W=\frac{\lambda^{N}e^{i\psi}}{(x_{1}+ix_{2})^{N}},$ (3.2)

with $\lambda>0$ a measure of the size of the lump and $\psi$ an internal
phase.

For $N>1$ the geodesic approximation is applicable to motion in this two-
dimensional moduli space as the kinetic energy is finite. Explicitly, allowing
$\lambda$ and $\psi$ to be time-dependent produces the kinetic energy

$T=\int\frac{|\dot{W}|^{2}}{(1+|W|^{2})^{2}}\,d^{2}x=(N^{2}\dot{\lambda}^{2}+\lambda^{2}\dot{\psi}^{2})\frac{\pi^{2}}{N^{2}\sin(\pi/N)},$
(3.3)

yielding a flat metric on this moduli space. However, for $N=1$ there is a
logarithmic divergence in the kinetic energy for motion in the moduli space,
making it a relevant example for the study of boundary metrics. As in the
previous section, the boundary metric is defined by the associated
renormalized kinetic energy $\widehat{T}$. Introducing polar coordinates
$r,\theta$ in the plane, let $T_{b}$ denote the regularized kinetic energy
defined by integration over the disc $r\leq b$, and set

$\widehat{T}=\lim_{b\to\infty}\frac{T_{b}}{\log
b}=\int_{0}^{2\pi}\frac{|\dot{W}|^{2}r^{2}}{(1+|W|^{2})^{2}}\bigg{|}_{r=\infty}d\theta\
=(\dot{\lambda}^{2}+\lambda^{2}\dot{\psi}^{2})2\pi,$ (3.4)

where the integration is over the boundary circle at infinity.

Note that the geodesics of this boundary metric fit naturally as the $N=1$
member of the geodesics of the general family of flat metrics

$ds^{2}=N^{2}d\lambda^{2}+\lambda^{2}d\psi^{2},$ (3.5)

proportional to those obtained from the geodesic approximation with $N>1.$ For
all $N$ the qualitative dynamics is the same. If the internal phase is
constant, $\dot{\psi}=0$, then the size $\lambda$ of the lump either increases
without limit, or shrinks to zero in a finite time, as this family of metrics
are geodesically incomplete. If there is motion in the internal phase then the
non-zero conserved quantity $\lambda^{2}\dot{\psi}$ prevents the lump from
shrinking to zero size and any initial shrinking is eventually reversed and
ultimately the lump expands without limit.

The qualitative features predicted by the geodesic approximation using the
boundary metric are in good agreement with the results from field theory
simulations [7, 8]. Numerical studies of the shrinking to zero size in finite
time for a $1$-lump confirm that a linear decrease in the size of the lump
remains an accurate description right up until times very close to the
critical time at which the size vanishes. The simulations are performed on a
numerical grid with a finite spatial extent, taken large enough so that
causality prevents any influence of the boundary on the field at the centre of
the lump during the simulation time [8]. This is equivalent to using an
indicator function with $\chi=1$ in the simulation region, with the value
outside this region being irrelevant for the evolution of the field at the
centre of the lump, within the limited time of the simulation.

The main difference between the field theory evolution and the geodesic
approximation is attributed to the fact that the latter neglects radiation,
which is essential in expelling excess energy from the core of the lump [9].
This applies equally well to all values of $N$, but for $N>1$ an initial slow
evolution tangential to the moduli space implies that the kinetic energy is
small compared to the potential energy, whereas this is not the case for
$N=1$. The geodesic approximation using the boundary metric models field
theory evolution in which the kinetic energy can be arbitrarily large,
depending upon the details of the indicator function $\chi,$ or equivalently
the spatial extent of the simulation region in numerical simulations.

A radially symmetric lump located at the origin is an appropriate system for
the application of the geodesic approximation using the boundary metric,
because all tangent vectors in the two-dimensional moduli space have infinite
length. This is complementary to previous studies of lump dynamics using the
geodesic approximation restricted to only those tangent vectors that have
finite length [10, 11]. Neither approach is applicable for motion in the full
moduli space of $N$-lumps, because of the existence of tangent vectors of both
finite and infinite length.

##  4 Kink dynamics

The simplest system to illustrate and investigate a boundary metric on a
soliton moduli space is a kink in one spatial dimension. Consider a modified
$\phi^{4}$ theory given by the Lagrangian density

${\cal L}=\frac{1}{2}\dot{\phi}^{2}\cosh(\kappa x)-\frac{1}{2}\phi^{\prime
2}-\frac{1}{2}(1-\phi^{2})^{2},$ (4.1)

where $\kappa$ is a non-negative constant and prime denotes differentiation
with respect to the spatial coordinate $x.$ The standard $\phi^{4}$ theory is
obtained for $\kappa=0.$ As the static sector is independent of $\kappa$ then
there is a one-dimensional moduli space of static kink solutions given by the
usual $\phi^{4}$ kink

$\phi(x,t)=\tanh(x-a),$ (4.2)

where the arbitrary real constant $a$ is a coordinate on the moduli space
corresponding to the position of the kink.

The geodesic approximation restricts the motion to the moduli space by
promoting $a$ to be time dependent, producing the tangent vector

$\dot{\phi}=-\dot{a}\,\mbox{sech}\,^{2}(x-a).$ (4.3)

If $0\leq\kappa<4$ then the associated kinetic energy

$T=\int_{-\infty}^{\infty}\frac{1}{2}\dot{\phi}^{2}\cosh(\kappa x)\,dx$ (4.4)

is finite and the standard geodesic approximation is applicable. For example,
if $\kappa=1$ then $T=\frac{\pi}{4}\dot{a}^{2}\cosh a$ with geodesics given by
$\ddot{a}=-\frac{1}{2}\dot{a}^{2}\tanh a.$ The upper red curve in Fig.1 shows
the geodesic with $a(0)=a_{0}=0$ and $\dot{a}(0)=v=0.1.$

Figure 1: The black curves give the position $a$ of the kink obtained from a
field theory simulation with an initial condition given by $a_{0}=0$ and
$v=0.1.$ The red curves are the corresponding geodesic approximation, using
the boundary metric in the case $\kappa=4$.

It is a simple task to compare the results of the geodesic approximation with
field theory simulations using the initial condition

$\phi(x,0)=\tanh(x-a_{0}),\quad\dot{\phi}(x,0)=-v\,\mbox{sech}\,^{2}(x-a_{0}).$
(4.5)

The results presented below use a numerical scheme with 6001 lattice points
and a lattice spacing $\Delta x=0.0025$ giving a simulation interval
$[-7.5,7.5]$. Spatial derivatives are approximated using fourth-order accurate
finite differences and time evolution is performed using a fourth-order Runge-
Kutta method with a timestep $\Delta t=\frac{1}{2}\Delta x.$ No flux boundary
conditions are imposed at the ends of the interval. The position of the kink,
$a$, is extracted using linear interpolation between the two lattice points at
which $\phi$ changes sign. In Fig.1 there is in fact an upper black curve
showing the results of the field theory simulation, but it is almost totally
obscured by the red curve, because the geodesic approximation is so accurate
in this case.

To see a slight difference between the geodesic approximation and field theory
simulations, consider the case $\kappa=2.$ The kinetic energy is
$T=\frac{4}{3}\dot{a}^{2}\cosh(2a)$ with geodesics given by
$\ddot{a}=-\dot{a}^{2}\tanh(2a).$ The middle red curve in Fig.1 again shows
the geodesic with $a_{0}=0$ and $v=0.1,$ and the black curve is the field
theory result that displays just enough difference to allow separate curves to
be distinguished.

Turning now to the case $\kappa=4$, then the kinetic energy is infinite if
motion is restricted to the moduli space, making it a suitable example for the
application of a boundary metric. Let $T_{b}$ denote the regularized kinetic
energy integrated over the interval $[-b,b].$ As $T_{b}$ is linearly divergent
as $b\to\infty,$ define the renormalized kinetic energy

$\widehat{T}=\lim_{b\to\infty}\frac{T_{b}}{b}=\lim_{x\to\infty}\bigg{[}\frac{1}{2}\dot{\phi}^{2}\cosh(4x)\bigg{]}+\lim_{x\to-\infty}\bigg{[}\frac{1}{2}\dot{\phi}^{2}\cosh(4x)\bigg{]}=8\dot{a}^{2}\cosh(4a).$
(4.6)

The geodesics of this boundary metric are given by
$\ddot{a}=-2\dot{a}^{2}\tanh(4a).$ The lower red curve in Fig.1 again shows
the geodesic with $a_{0}=0$ and $v=0.1,$ and the black curve is the field
theory result. This demonstrates that geodesics of the boundary metric provide
a reasonable description of kink dynamics. As earlier, the field theory
computation should be regarded as equivalent to using an indicator function
with $\chi=1$ in the simulation region, with the value outside this region
being irrelevant for the evolution of the kink over the timescale of the
simulation. This has been justified by performing the same computation using
double the spatial region and confirming that there is no discernible
difference in the evolution of the position of the kink. Changing the boundary
conditions from no flux to Dirichlet also yields no perceptible change,
providing another numerical check.

The expectation is that radiation is the main source of the discrepancy
between field theory evolution and the geodesic approximation. In theories
with a boundary metric, the slow motion of a soliton no longer implies that
the kinetic energy is a small fraction of the total energy. A moving soliton
may therefore generate a considerable amount of radiation and energy
conservation does not provide such a useful constraint to prevent the field
configuration from deforming far from the moduli space. To illustrate this
issue, the black curve in the left-hand-side plot in Fig.2 shows the field
$\phi$ at the end of the simulation ($t=200$) associated with the data from
Fig.1. Also shown (red curve) is the static kink with the same position. This
reveals that the two fields agree in the core of the kink, but that there are
significant deformations in the regions around $x\approx\pm 3.$ To check that
these deformations are consistent with an initial burst of radiation leaving
the core of the kink, consider the equation for the characteristics that pass
through the origin at $t=0$

$\int_{0}^{x}\sqrt{\cosh(4X)}\,dX=\pm t.$ (4.7)

Evaluating this equation at $t=200$ yields $x=\pm 3.17$, which is indeed
consistent with the location of the deformation.

Figure 2: Left: The moving kink (black curve) at $t=200$ and the static kink
(red curve) with the same position. Right: The moving kink (black curve) at
$t=200$. The field at $t=200$ from the initial condition $\phi=-1$ (blue
curve) and $\phi=1$ (orange curve) with the initial time derivative
$\dot{\phi}$ the same as for the moving kink.

Further evidence for the nature of the deformation is provided by considering
initial conditions given by the vacuum fields $\phi(x,0)=-1$ and
$\phi(x,0)=1$, with the initial time derivative $\dot{\phi}(x,0)$ the same as
that of the moving kink. The results are displayed in the right-hand-side plot
in Fig.2, where the field at time $t=200$ is shown in blue for $\phi(x,0)=-1$
and orange for $\phi(x,0)=1$. For comparison, the moving kink field is
displayed again as the black curve. This clearly shows that the moving kink
may be decomposed into a soliton component, well-approximated by a point in
the moduli space, and a substantial radiation component, which nonetheless has
little influence on the soliton core.

##  5 Conclusion

The concept of a boundary metric allows the geodesic approximation of soliton
dynamics to be applied to situations in which motion in the soliton moduli
space yields infinite kinetic energy. Numerical studies of a low-dimensional
example provide support for this approach and suggest that it provides a
reasonable description of the soliton component of the dynamics, even though
there is more radiation than in the usual application of the geodesic
approximation. This work provides motivation for a future study of the Braam-
Austin metric on the moduli space of hyperbolic monopoles and an
interpretation of the results in terms of hyperbolic monopole scattering. To
date, only the metric on the moduli space of a single hyperbolic monopole has
been calculated. Numerical field theory studies of hyperbolic monopole
dynamics would also be interesting, for comparison to predictions from the
geodesic approximation, although this is a significant numerical challenge.

It would be interesting to see if any rigorous mathematical results can be
derived for the geodesic approximation with a boundary metric, along the lines
of the results obtained for the usual geodesic approximation [3, 4].

## References

  * [1] N. S. Manton, A remark on the scattering of BPS monopoles, Phys. Lett. B110, 54 (1982).
  * [2] M. F. Atiyah and N. J. Hitchin, The Geometry and Dynamics of Magnetic Monopoles, Princeton University Press, 1988.
  * [3] D. Stuart, Dynamics of abelian Higgs vortices in the near Bogomolny regime, Commun. Math. Phys. 159, 51 (1994).
  * [4] D. Stuart, The geodesic approximation for the Yang-Mills-Higgs equations, Commun. Math. Phys. 166, 149 (1994).
  * [5] M. F. Atiyah, Magnetic monopoles in hyperbolic spaces, in M. Atiyah: Collected Works, vol. 5, Oxford, Clarendon Press, 1988.
  * [6] P. J. Braam and D. M. Austin, Boundary values of hyperbolic monopoles, Nonlinearity 3, 809 (1990).
  * [7] R. A. Leese, M. Peyrard and W. J. Zakrzewski, Soliton stability in the $O(3)$ $\sigma$-model in (2+1) dimensions, Nonlinearity 3, 387 (1990).
  * [8] B. Piette and W. J. Zakrzewski, Shrinking of solitons in the (2+1)-dimensional $S^{2}$ sigma model, Nonlinearity 9, 897 (1996).
  * [9] P. Bizoń, T. Chmaj and Z. Tabor, Formation of singularities for equivariant $(2+1)$-dimensional wave maps into the 2-sphere, Nonlinearity 14, 1041 (2001).
  * [10] R. S. Ward, Slowly-moving lumps in the $\mathbb{CP}^{1}$ model in (2+1) dimensions, Phys. Lett. B158, 424 (1985).
  * [11] R. A. Leese, Low energy scattering of solitons in the $\mathbb{CP}^{1}$ model, Nucl. Phys. B344, 33 (1990).


# Modeling, scientific computing and optimal control
for renewable energy systems with storage

Nicola Cantisani, Tobias K. S. Ritschel, Christian A. Thilker, Henrik Madsen
and John Bagterp Jørgensen N. Cantisani, T.K.S. Ritschel, C.A. Thilker, H.
Madsen, J.B. Jørgensen are with Department of Applied Mathematics and Computer
Science, Technical University of Denmark, DK-2800 Kgs. Lyngby, Denmark

###### Abstract

This paper presents models for renewable energy systems with storage, and
considers its optimal operation. We model and simulate wind and solar power
production using stochastic differential equations as well as storage of the
produced power using batteries, thermal storage, and water electrolysis. We
formulate an economic optimal control problem, with the scope of controlling
the system in the most efficient way, while satisfying the power demand from
the electric grid. Deploying multiple storage systems allows flexibility and
higher reliability of the renewable energy system.

##  I Introduction

As society has progressively become industrialized and the use of energy has
been increasing over the last century, greenhouse gas emissions have risen at
an alarming pace, leading to global warming. The biggest contribution to
greenhouse emissions is carbon dioxide, which is mainly released when burning
fossil fuels. Consequently, we need to push development towards the use of
renewable energy sources (RESs), such as wind and solar. However, solar and
wind cannot be controlled and are highly volatile/intermittent. Energy storage
systems (ESSs) can be a solution to this problem, by serving as a buffer when
the power production is low or absent, and making the energy more dispatchable
[1, 2]. Fig. 1 presents the main ESSs, providing information about the
capacity, the discharging time at rated power, and the energy conversion
efficiency. Even though ESSs with bigger capacities generally have lower
energy efficiencies, they can support energy systems for longer periods of
time.

Figure 1: Characterization of the main energy storage systems (ESS) by storage
capacity and discharging time at rated power (i.e. maximum output power). The
percentages refer to energy efficiency.

In this paper, we demonstrate mathematical modeling, simulation, and control
of such systems. Based on the scientific literature, we model the uncertainty
coming from the wind and the sun using stochastic differential equations. We
express the power extraction from wind turbines as an optimization problem and
present models for energy storage. We develop and propose simple models for
battery and thermal storage, based on energy balance equations. Electrolysis
is modeled using an electrochemical model, as in [3, 4]. We will consider a
hybrid energy system consisting of stochastic power production from wind and
solar and the combination of battery, thermal storage, and hydrogen for energy
storage.

Ultimately, we want to achieve optimal operation of such a system. Therefore,
an optimal control problem is formulated. The goal is to maximize the profit
generated from selling electricity, heat, and hydrogen, while assuring that
the grid’s power demand is always satisfied. Fig. 2 provides an overview of
the entire energy system and the controllable power flows.

Figure 2: Overview of the controllable renewable energy system.

There is growing research on how to best control a renewable energy system
with ESSs, where model predictive control (MPC) is often the preferred choice.
[5, 6, 7, 8] show how to control a system with one stochastic power source
(wind or solar) and a single storage device, often a battery. Integrating and
controlling different storage solutions together is also attempted in the
literature [9, 10]. The control strategy proposed by the mentioned resources
is based on maximizing the utilization of renewable resources. In this work,
we explore the possibility of integrating and controlling three different
ESSs, maximizing the economic revenue when such a system is interfaced with
the market.

The paper is structured as follows. Sections II and III present models for
power production from wind turbines and solar photovoltaic (PV) panels,
respectively. In Sections IV, V and VI we model energy storage via battery,
thermal storage and water electrolysis, respectively. Section VII introduces
the optimal control problem. Finally, Section VIII presents a simulation of
the system.

##  II Wind energy

In this section, we present control relevant simulation models for horizontal
axis wind turbines (HAWT). The model consists of a model for the wind speed
and one for the wind turbine.

###  II-A Wind speed model

The wind speed is modeled as the sum of two components [11]:

$v(t)=v_{m}(t)+v_{t}(t).$ (1)

$v_{m}(t)$ is the mean wind speed and ${v}_{t}(t)$ is the turbulent wind
speed. Each part is modeled by a stochastic differential equation (SDE):

$\displaystyle dv_{t}(t)$
$\displaystyle=-\frac{\pi\,v_{m}(t)}{2L}v_{t}(t)\,dt+\sqrt{\frac{\pi\,v_{m}(t)^{3}\,t_{i}^{2}}{L}}\,dw_{1}(t),$
(2a) $\displaystyle dv_{m}(t)$ $\displaystyle=\sigma_{2}\,dw_{2}(t).$ (2b)

$L$ is the turbulence length and $t_{i}$ is the turbulence intensity. The
parameter values for the wind speed model are reported in Table I. Notice that
the mean wind speed can also be set and fixed for simulation purposes. Observe
that we do not consider the wind direction here, but it is relevant for a wind
farm.

TABLE I: Parameters of the wind speed model [11]. Parameter | Value | Unit
---|---|---
$L$ | 170.1 | m
$t_{i}$ | 0.2 | -
$\sigma_{2}$ | $\sqrt{4/600}$ | m $\cdot$ s-3/2

###  II-B Power extraction

The available power in the wind is

$P_{w}(t)=\frac{1}{2}\rho\pi R^{2}v(t)^{3}.$ (3)

$\rho$ is the density of the air and $R$ is radius of the blade. The tip speed
ratio is the ratio between the tangential speed at the tip of a blade of the
turbine and the wind speed, i.e.

$\lambda(t)=\frac{R\Omega_{r}(t)}{v(t)}.$ (4)

$\Omega_{r}(t)$ is the rotational speed of the rotor. The blade pitch angle,
$\theta$, is defined as the angle between the chord of the airfoil and the
rotational plane. The extracted power by the wind turbine rotor is only a part
of the full available power, $P_{w}(t)$. The fraction of the extracted power
depends on the turbine’s power coefficient, $C_{P}(\lambda(t),\theta(t))$,
i.e.

$P_{r}(t)=P_{w}(t)C_{P}(\lambda(t),\theta(t)).$ (5)

The power coefficient values are table-based. The torque is

$Q_{r}(t)=\frac{P_{r}(t)}{\Omega_{r}(t)}.$ (6)

###  II-C Stationary optimal operation

TABLE II: Relevant NREL 5MW parameters and constraints [12]. Parameter | Value | Unit | Description
---|---|---|---
$\rho$ | 1.225 | kg/m3 | Density of air
$R$ | 62.94 | m | Length of rotor blade
$\eta_{g}$ | 94.4% | - | Power generator efficiency
$v_{cut-in}$ | 3 | m/s | Cut in velocity
$v_{rated}$ | 11.4 | m/s | Rated velocity
$v_{cut-out}$ | 25 | m/s | Cut out velocity
$\Omega_{r,min}$ | 6.9 | rpm | Minimum rotor speed
$\Omega_{r,max}$ | 12.1 | rpm | Maximum rotor speed
$\theta_{min}$ | -5 | deg | Minimum pitch angle
$\theta_{max}$ | 25 | deg | Maximum pitch angle

The maximum power extraction of the turbine is achieved as the solution of the
following optimization problem, assuming steady state conditions (i.e.
constant wind speed).

$\displaystyle\max_{\theta,\Omega_{r}}\quad$ $\displaystyle\frac{1}{2}\rho
R^{2}v^{3}C_{P}(\lambda(\Omega_{r},v),\theta)$ (7a)
$\displaystyle\mathrm{s.t.}$
$\displaystyle\theta_{\min}\leq\theta\leq\theta_{\max},$ (7b)
$\displaystyle\Omega_{r,\min}\leq\Omega_{r}\leq\Omega_{r,\max}.$ (7c)

Notice that the problem depends on the wind speed, $v$. Every value of the
wind speed will lead to a different solution. The solution for every wind
speed between the cut-in and cut-out velocity produces the power curve of the
wind turbine. The actual generated power by the generator, $P_{g}$, in the
turbine is obtained by multiplying the power produced by the rotor, $P_{r}$,
by the efficiency of the electrical generator, $\eta_{g}$. The turbine is, in
general, constrained between a range of pitch angles and rotational speed,
according to the manufacturer’s regulations, hence constraints (7b) and (7c).
We can reformulate the optimization problem, emphasizing the constraint using
$\lambda$, and discarding the constant part of the objective function.

$\displaystyle\max_{\theta,\lambda}\quad$ $\displaystyle
C_{P}(\lambda,\theta)$ (8a) $\displaystyle s.t.$
$\displaystyle\theta_{\min}\leq\theta\leq\theta_{\max},$ (8b)
$\displaystyle\lambda_{\min}(v)\leq\lambda\leq\lambda_{\max}(v),$ (8c)

where

$\lambda_{\min}(v)=\frac{R\Omega_{r,\min}}{v},\quad\lambda_{\max}(v)=\frac{R\Omega_{r,\max}}{v}.$
(9)

The solution, for each wind speed $v$, consists of $\theta^{*}=\theta^{*}(v)$,
$\lambda^{*}=\lambda^{*}(v)$, and the derived quantities
$C_{P}^{*}(v),\,P_{w}^{*}(v),\,P_{g}^{*}(v)$ and $\Omega_{r}^{*}(v).$

We consider an NREL 5MW wind turbine [12] and solve the problem, using the
parameters in Table II. The $C_{p}(\lambda,\theta)$ values and the solution of
the optimal stationary problem (8) for every $v$ is shown in Fig. 3. Fig. 4
reports the entire solution as a function of the wind speed.

Figure 3: $C_{p}(\lambda,\theta)$ values of a NREL 5MW wind turbine and the
solution of the optimal stationary problem 8, for $v\in[v_{cut-in},v_{cut-
out}]$. Figure 4: The optimal pitch angle $\theta^{*}$, tip speed ratio
$\lambda^{*}$, power coefficient $C_{p}^{*}$, rotor speed $\Omega_{r}^{*}$ and
generated power $P_{g}^{*}$ as a function of the wind speed.

##  III Solar Energy

The global radiation, $I_{G}$, is the total amount of short-wave radiation
received from the sun by a surface horizontal near the ground. This value is
of particular interest to photovoltaic installations and depends on both
direct radiation, $I_{N}$, and diffuse radiation, $I_{D}$.

In the literature, the solar radiation is commonly calculated as the intensity
that would be obtained in clear sky conditions, multiplied by a factor that
depends on the cloud cover [13, 14, 15, 16, 17]. Consequently, simulation
models for direct, diffuse, and global solar radiation have been suggested in
the literature where stochastic models for cloud cover have been used as the
driving entity [18, 14, 15]. Such methods have been used for simulation and
prediction of weather observations in optimal control of heat supply to
buildings [19].

In this paper, the cloud cover $\kappa(t)$ is modeled using the mean-reversion
SDE model from [15]

$\displaystyle\begin{split}d\kappa(t)&=\theta(\kappa(t))[\mu(\kappa(t))-\kappa(t)]\,dt\\\
&\qquad+\sigma\kappa(t)[1-\kappa(t)]\,d\omega_{\kappa}(t),\end{split}$ (10a)
$\displaystyle\mu(\kappa(t))$
$\displaystyle=\frac{\exp[P_{7}(\kappa(t))]}{1+\exp[P_{7}(\kappa(t))]},$ (10b)
$\displaystyle\theta(\kappa(t))$
$\displaystyle=\tilde{\theta}\sqrt{\kappa(t)(1-\kappa(t))}.$ (10c)

The term $P_{7}(.)=\sum_{k=1}^{7}p_{k}L_{k}(.)$ is the linear combination of
the first seven Legendre polynomials. The equation provides a cloud cover such
that $\kappa(t)\in[0,1]$. A transformation of the data is needed to map the
data to the Okta scale ($0$ = absence of cloud cover, $8$ = total coverage of
the sky). For simulation purposes, the mean value, $\mu$, can also be fixed.

We now present the dynamical models that describe the solar irradiance [18].
The model for the direct radiation, $I_{N}$, is

$\displaystyle\frac{I_{N}(\kappa(t),h(t))}{I_{0,N}(h(t))}$
$\displaystyle=f_{N}(\kappa(t))+\epsilon_{N}(t),$ (11a) $\displaystyle
I_{0,N}(h(t))$ $\displaystyle=a_{N}[1-\exp(-b_{N}h(t))],$ (11b) $\displaystyle
f_{N}(\kappa(t))$ $\displaystyle=\frac{a^{\prime}+b^{\prime}\cos(\frac{2\pi
d_{t}}{365}+c^{\prime})}{1+\exp(\kappa(t)-\alpha_{N})},$ (11c)
$\displaystyle\epsilon_{N}(t)$
$\displaystyle\sim\mathcal{N}\left(0,r_{N}^{2}\frac{\sigma_{D,\kappa(t)}^{2}}{\sin^{2}h(t)I_{0,N}(h(t))^{2}}dt\right).$
(11d)

$I_{0,N}$ is the direct radiation in clear skies (non-stochastic), $d_{t}$ is
the day of the year and $h(t)$ is the sun elevation angle at the location of
interest. The cloud cover is here provided in Okta scale
($\kappa(t)\in[0,8]$). More variability is expected when the sun elevation,
$h(t)$, is low. Madsen and Thyregod [20] suggest a stochastic model for direct
radiation in clear skies which also takes into account the auto-correlation
caused by the persistence in variations of the transparency of the atmosphere
due to ozone, dust, water, etc.

The model for the diffuse radiation, $I_{D}$, is

$\displaystyle\frac{I_{D}(\kappa(t),h(t))}{I_{0,D}(h(t))}$
$\displaystyle=f_{D}(\kappa(t))+\epsilon_{D}(t),$ (12a) $\displaystyle
I_{0,D}(h(t))$ $\displaystyle=c_{D}+a_{D}[1-\exp(-b_{D}h(t))],$ (12b)
$\displaystyle\begin{split}f_{D}(\kappa(t))&=a_{0}(t)+a_{1}(t)(1-\kappa(t)/8)\\\
&\qquad+a_{2}(\kappa(t)/8)^{\alpha_{D}}(1-\kappa(t)/8),\end{split}$ (12c)
$\displaystyle\epsilon_{D}(t)$
$\displaystyle\sim\mathcal{N}(0,\sigma_{N,\kappa(t)}^{2}r_{D}^{2}\,dt).$ (12d)

$I_{0,D}$ is the diffuse radiation in clear skies and

$\displaystyle a_{0}=r_{1}+r_{2}\cos\left(\frac{2\pi
d_{t}}{365}+r_{3}\right),$ (13a) $\displaystyle
a_{1}=k_{1}+k_{2}\cos\left(\frac{2\pi d_{t}}{365}+k_{3}\right).$ (13b)

The variances, $\sigma_{N,\kappa(t)}$ and $\sigma_{D,\kappa(t)}$, depend on
the value of the cloud cover, $\kappa(t)$ (data is available in [18]).

The global radiation is

$I_{G}(t)=I_{N}(\kappa(t),h(t))\sin h(t)+I_{D}(\kappa(t),h(t)).$ (14)

Notice that, due to the dependency of the solar elevation in the variance
structure of (11), the variance of the global radiation becomes independent of
the solar elevation.

The power generated by solar radiation is

$P_{s}(t)=\eta_{s}A_{s}I_{G}(t),$ (15)

where $A_{s}$ is the area of PVs and $\eta_{s}$ the PVs efficiency. For real-
time solar power forecasting tools, like
SOLARFOR111https://enfor.dk/services/solarfor/, the efficiency depends, e.g.,
on the dirtiness of the PV panels. Table III reports all the parameters of the
models.

TABLE III: Parameters of the solar power model. Cloud cover model [15]
---
Parameter | $\tilde{\theta}$ | $\sigma$ | $p_{1}$ | $p_{2}$ | $p_{3}$
Value | 0.187 | 0.835 | -53.1 | 14.6 | -42.3
Parameter | $p_{4}$ | $p_{5}$ | $p_{6}$ | $p_{7}$ |
Value | 8.8 | -58.1 | -30.3 | -45.7 |
Direct radiation model [18]
Parameter | $a_{N}$ | $b_{N}$ | $r_{N}$ | $\alpha_{N}$ | $a^{\prime}$
Value | 842.3 | 0.0614 | 1.0430 | 4.6368 | 1.1354
Parameter | $b^{\prime}$ | $c^{\prime}$ |  |  |
Value | 0.1965 | -0.2571 |  |  |
Diffuse radiation model [18]
Parameter | $a_{D}$ | $b_{D}$ | $c_{D}$ | $r_{1}$ | $r_{2}$
Value | 161.1 | 0.0333 | 3.68 | 0.7067 | -0.2456
Parameter | $r_{3}$ | $k_{1}$ | $k_{2}$ | $k_{3}$ | $a_{2}$
Value | 0.5625 | 0.1946 | 0.1549 | 0.6034 | 6.7033
Parameter | $\alpha_{D}$ | $r_{D}$ |  |  |
Value | 2.2993 | 1.0170 |  |  |

##  IV Battery

We model the accumulated energy in a battery considering its energy balance.
The change in the energy in the battery $E_{b}(t)$ over time is expressed by
the following differential equation

$\dot{E}_{b}(t)=-\alpha_{b}E_{b}(t)+\eta_{b,in}P_{b,in}(t)-\frac{1}{\eta_{b,out}}P_{b,out}(t),\\\
$ (16)

where $P_{b,in}(t)$ is the incoming power and $P_{b,out}(t)$ is the outgoing
power. The factors $\eta_{b,in}$ and $\eta_{b,out}$ are the charging and
discharging efficiencies. The first term of the equation accounts for the loss
of charge over time. Notice that the charge is constrained between a minimum
and a maximum value,

$E_{b,min}\leq E_{b}(t)\leq E_{b,max}.$ (17)

##  V Thermal storage

Thermal energy storage (TES) is a relatively inexpensive and efficient way to
store big amounts of energy for a long period of time. Even though the TES
systems can differ quite significantly, they all use the ability of some
material to accumulate energy through heat exchange. When energy is needed,
the material releases heat (cools down), which can be used in several ways.

TES can be either latent or sensible, according to whether a phase transition
of the medium is allowed or not. We consider the latter case. In this case, a
heat exchange with the material always implies a change in temperature. We
express the heat exchange using the first law of thermodynamics (in absence of
work):

$\Delta U=Q=m\,c_{P}\,\Delta T.$ (18)

$m$ is the mass and $c_{P}$ is the specific heat capacity. Common materials
used for TES are molten salts. [21] reports density, fusion temperature,
$T_{\min}$, decomposition temperature, $T_{\max}$, and specific heat capacity,
$c_{P}(T)$, for the most commonly used molten salts. We can compute how much
energy it is possible to store in a unit of mass of the material, considering
the specific storage capacity (J/kg), as

$\Delta u=\int_{T_{\min}}^{T_{\max}}c_{p}(T)dT.$ (19)

As for the battery, the thermal storage model is based on energy balance. The
factors contributing to an internal energy change are incoming power,
$P_{t,in}$, heat flow (loss) to the ambient, $Q_{ta}$, and outgoing heat flow,
$Q_{t,out}$.

$\dot{U}(t)=\eta_{t,in}P_{t,in}(t)-Q_{ta}(t)-Q_{t,out}(t).$ (20)

Expressing all the terms, we get the full model

$\begin{split}&m\,c_{P}(T(t))\dot{T}(t)=\\\
&\quad\eta_{t,in}P_{t,in}(t)-(UA)_{t}(T(t)-T_{a}(t))-Q_{t,out}(t),\end{split}$
(21)

subject to

$T_{\min}\leq T(t)\leq T_{\max}.$ (22)

$(UA)_{t}$ is the overall conductivity of the storage tank and depends on the
material and the exchange area. $T_{a}(t)$ is the ambient temperature.

##  VI Water electrolysis

###  VI-A Thermodynamics

The general reaction of electrolysis of water is

$\mathrm{2H_{2}O(l)+\text{Electrical power}\longrightarrow
2H_{2}(g)+O_{2}(g)}.$ (23)

The thermodynamics of the reaction needs to be considered. The energy that is
required for the reaction to happen is the enthalpy of the process

$\Delta H=\Delta G+T\Delta S.$ (24)

The reaction is endothermic because $\Delta H>0$ and it is equal to $\Delta
H^{0}=285.8$ kJ/mol at standard conditions. The change in Gibbs free energy is
equal to $\Delta G^{0}=237.2$ kJ/mol at standard conditions. A positive value
in the change of Gibbs free energy indicates that the reaction does not happen
spontaneously, but it is spontaneous in the opposite direction. To let the
reaction happen in this direction, electric work must be supplied, and it is
equal to the change in Gibbs free energy. The remaining energy requirement is
brought as heat, because of the irreversibilities of the process and it is
equal to $Q=T\Delta S=48$ kJ/mol. We can re-write the total reaction,
including the energy demands

$\begin{split}\mathrm{H_{2}O(l)+237.2\,\text{kJ}_{\text{electricity}}/\text{mol}+48.6\,\text{kJ}_{\text{heat}}/\text{mol}}\\\
\mathrm{\longrightarrow H_{2}(g)+\frac{1}{2}O_{2}(g)}.\end{split}$ (25)

The lowest voltage that is needed to initiate the process is called reversible
voltage. It is given by

$U_{rev}=\frac{\Delta G}{zF}.$ (26)

$z=2$ is the number of electrons transferred and $F=96485$ C/mol is the
Faraday constant. At standard conditions, $U_{rev}^{0}=1.229$ V.

In order for the reaction to happen, heat must also be provided. In real
electrolyzers, additional electric work is provided, instead of heat. In this
case, the minimum voltage requirement is higher, and it is called the
thermoneutral voltage

$U_{tn}=\frac{\Delta H}{zF}.$ (27)

At standard conditions, $U_{tn}^{0}=1.481$ V. Both the reversible and
thermoneutral voltages are state functions.

###  VI-B Electrochemical model

We consider the electrochemical model from [4]. This has been used extensively
in the scientific literature, because of its simplicity and its empiric basis.
We consider only alkaline electrolyzers. We model the kinetics at the
electrodes, i.e. we investigate the relationship between current intensity and
voltage in the electrolyzer ($I-U$). In general, we can express the real total
voltage in one electrolytic cell as the sum of several terms

$U_{cell}=U_{rev}+U_{ohm}+U_{act}+U_{con}.$ (28)

$U_{act}$ is the activation overvoltage, which is caused by the anode and
cathode half-reactions kinetics [22]. $U_{ohm}$ is the ohmic overvoltage. This
is due to ohmic losses in the cell, generated by the resistance to the
electron flux in the components of the electrolyzer. $U_{con}$ is
concentration overvoltage, which is negligible for alkaline electrolyzers. [4]
proposed the following electrochemical model for ohmic and activation
overvoltages

$\displaystyle U_{ohm}$ $\displaystyle=(r_{1}+r_{2}T)\frac{I}{A},$ (29a)
$\displaystyle U_{act}$
$\displaystyle=s\log\left(\left(t_{1}+\frac{t_{2}}{T}+\frac{t_{3}}{T^{2}}\right)\frac{I}{A}+1\right).$
(29b)

$A$ is the electrode area of a single cell and $T$ is the temperature of the
electrodes. The parameters $r_{1}$ and $r_{2}$ represent the ohmic resistance,
while $s$ and $t_{1},t_{2},t_{3}$ are the activation overvoltage coefficients.
Table IV reports the parameter estimates. The total power in a stack of
$n_{c}$ electrolytic cells connected in series is

$P_{el}=I\,n_{c}\,U_{cell}.$ (30)

The total production flow rate of hydrogen of an electrolyzer with $n_{c}$
cells connected in series is a function of the current, according to

$f_{H_{2}}=n_{c}\left(\eta_{F}\frac{I}{zF}\right).$ (31)

$\eta_{F}$ is called Faraday efficiency (or current efficiency). This is
experimentally defined as the ratio of the ideal hydrogen production rate to
the actual hydrogen output from the electrolyzer. This is due to the fact that
not all the electrons contribute to the water-splitting reaction. From the
stoichiometry of the reaction, we notice that

$f_{H_{2}O}=f_{H_{2}}=2f_{O_{2}}.$ (32)

The following model for the Faraday efficiency is proposed [4]

$\eta_{F}=\frac{(I/A)^{2}}{f_{1}+(I/A)^{2}}f_{2}.$ (33)

$f_{1}$ and $f_{2}$ are parameters to be estimated empirically and they depend
on the temperature. Table IV presents the estimated values for
$T=80\,{}^{\circ}\text{C}$.

Remark. The overvoltage $(U_{cell}-U_{tn})$ is responsible for generating heat
in the electrolyzer (this is equal to $Q_{gen}=n_{C}(U_{cell}-U_{tn})I$).
Therefore, a cooling system is generally deployed to keep the temperature
constant. The exchanged heat could be stored or used, e.g., for district
heating. We assume that the temperature in the electrolyzer is always kept
constant.

###  VI-C Hydrogen tank

We model the hydrogen (and oxygen) tank considering the mass balance equation

$\dot{n}_{H_{2}}^{stored}=f_{H_{2}}^{in}-f_{H_{2}}^{out}.$ (34)

The pressure in the tank can be obtained using a cubic equation of state
(EOS), e.g. the Peng-Robinson EOS.

TABLE IV: Electrochemical model parameters values. Parameter | Value | Unit
---|---|---
U-I model [3]
$r_{1}$ | 0.8 | $\Omega\,\text{cm}^{2}$
$r_{2}$ | -0.00763 | $\Omega\,\text{cm}^{2}\,{}^{\circ}\text{C}^{-1}$
$s$ | 0.1795 | V
$t_{1}$ | 20 | $\text{cm}^{2}\,\text{A}^{-1}$
$t_{2}$ | 0.1 | $\text{cm}^{2}\,{}^{\circ}\text{C}\,\text{A}^{-1}$
$t_{3}$ | 3.5 $\times 10^{5}$ | $\text{cm}^{2}\,{}^{\circ}\text{C}^{2}\,\text{A}^{-1}$
Faraday efficiency model [4]
$f_{1}(T=80\,{}^{\circ}\text{C})$ | 250 | $\text{mA}^{2}\text{cm}^{-4}$
$f_{2}(T=80\,{}^{\circ}\text{C})$ | 0.980 | -

##  VII Optimal control

We want to control the renewable energy system shown in Fig 2. The goal is to
satisfy electricity demand from the grid and to maximize the profit generated
from selling excess electricity, heat, and hydrogen. We also allow re-
conversion of heat and hydrogen to electricity. This is achieved using a fuel
cell for the hydrogen and a steam turbine for the heat coming from the thermal
storage. Electricity can also be bought from the market and stored in the
battery, when needed. The setup of the system is described in Fig. 2. We
approximate the power generated by a generic fuel cell with $\eta_{fc}$
efficiency and hydrogen flow rate $d_{H_{2}}(t)$:

$d_{fc}(t)=141800\cdot\eta_{fc}\cdot MM_{H_{2}}\cdot d_{H_{2}}(t).$ (35)

$MM_{H_{2}}$ is the molar mass of hydrogen. We introduce auxiliary variables
to simplify the formulation of the optimal control problem. Given the total
produced power, $P_{tot}(t)$, by solar PVs and wind turbines and the
electricity demand, $d_{el}(t)$, we define the excess power,
$P_{sto}(t)=P_{tot}(t)-d_{el}(t)$ (when $P_{tot}(t)>d_{el}(t)$), which is to
be stored. We define the variable $d_{sto}(t)$ as the power demand from the
storage systems, when the produced power is insufficient
($d_{el}(t)>P_{tot}(t)$). Using these variables, the optimization problem only
controls the energy surplus and deficit, while what is produced until surplus
goes directly to the grid without passing through the storage. We define the
following economic nonlinear optimal control problem

$\displaystyle\max_{u(t),x(t)}\quad$
$\displaystyle\phi=\phi([x(t),u(t)]_{t_{0}}^{t_{f}})$ (36a) s.t.
$\displaystyle x(t_{0})=\hat{x}_{0},$ (36b) $\displaystyle
P_{sto}(t)=P_{b,in}(t)+P_{t,in}(t)+P_{el}(t),$ (36c) $\displaystyle
d_{sto}(t)=d_{b}(t)+d_{t}(t)+d_{fc}(t),$ (36d)
$\displaystyle\dot{E}_{b}(t)=-\alpha_{b}E_{b}(t)+\eta_{b,in}\left(P_{b,in}(t)+P_{b,bo}(t)\right)$
$\displaystyle\quad\quad\quad\quad-\frac{1}{\eta_{b,out}}\left(P_{b,out}(t)+d_{b}(t)\right),$
(36e) $\displaystyle
m\,c_{P}(T)\dot{T}(t)=\eta_{t,in}P_{t,in}(t)-(UA)_{t}(T(t)-T_{a})$
$\displaystyle\quad\quad\quad\quad-
Q_{t,out}(t)-\frac{1}{\eta_{steam}}d_{t}(t),$ (36f)
$\displaystyle\text{Electrolyzer model \eqref{eq:el1}, \eqref{eq:el2},
\eqref{eq:el3}, \eqref{eq:el4}},$ $\displaystyle\text{Fuel cell model
\eqref{eq:fuel_cell}},$
$\displaystyle\dot{n}_{H_{2}}^{stored}(t)=f_{H_{2}}^{in}(t)-f_{H_{2}}^{out}(t)-d_{H_{2}}(t),$
(36g) $\displaystyle x_{min}\leq x(t)\leq x_{max},$ (36h) $\displaystyle
u_{min}\leq u(t)\leq u_{max}.$ (36i)

The objective function is

$\begin{split}\phi&=\int_{t_{0}}^{t_{f}}c^{e}(t)P_{b,out}(t)\,dt-\int_{t_{0}}^{t_{f}}c^{e}(t)P_{b,bo}(t)\,dt+\\\
&\int_{t_{0}}^{t_{f}}c^{h}(t)Q_{t,out}(t)\,dt+MM_{H_{2}}\int_{t_{0}}^{t_{f}}c^{H_{2}}(t)f_{H_{2}}^{out}(t)\,dt\\\
&+\text{Profit to-go}.\end{split}$ (37)

$c^{e}(t),c^{h}(t)$ and $c^{H_{2}}(t)$ are the cost functions of electricity,
heat and hydrogen over time, respectively. The subscript $out$ indicates
outgoing flows to the markets (energy sold), while $in$ indicates incoming
flows to the storage devices from wind and solar. Constraint (36c) dispatches
the energy surplus to the storage devices. Constraint (36d) ensures that power
is withdrawn from the storage devices when the production is not sufficient.
This constraint can be formulated as a soft constraint, to ensure feasibility
of the problem even when there is not enough energy in the system. Notice that
the nonlinearity is due to the dynamics in the electrolyzer.

The states consist of the dynamics in the ESSs:

$x(t)=\Big{[}E_{b}(t),T(t),n_{H_{2}}^{stored}(t)\Big{]}.$ (38)

The vector of manipulated variables is

$\begin{split}u(t)=&\Big{[}P_{b,in}(t),P_{t,in}(t),P_{el}(t),P_{b,out}(t),P_{b,bo}(t),\\\
&Q_{t,out}(t),f_{H_{2}}^{out}(t),I(t)\Big{]}.\end{split}$ (39)

The disturbances of the system are

$d(t)=\Big{[}c^{e}(t),c^{h}(t),c^{H_{2}}(t),P_{sto}(t),d_{sto}(t)\Big{]}.$
(40)

##  VIII Simulation results

We present a 3 days simulation. A constant power demand $d_{el}(t)$ of 4 MW is
assumed. The configuration of the system is the following: one NREL 5 MW
turbine, a 6 MW solar PV park, 50 MWh thermal storage, 5 MWh battery, 2.4 MW
electrolyzer, and a 1000 kg (30 m3) $\simeq$ 39 MWh hydrogen tank. The time
series for the prices are randomly generated from normal distributions. We
remark that modeling energy prices is outside of the scope of this paper. Our
main purpose is to demonstrate optimal control of renewable energy systems
with storage. Fig. 5 shows the simulated power production (solar, wind and
total).

Figure 5: Simulation of the amount of power produced from wind and solar over
three days. Only the mean wind speed is considered.

The optimal control problem (36) is discretized, using 10 min as sampling time
and 1 h as control interval. The problem is set up using CasADi [23] and
solved using IPOPT. Fig. 6 presents the states $x(t)$ of the system.

Figure 6: Solution of the optimal control problem: states $x(t)$.

Fig. 7 displays how the electricity demand is satisfied by the control
program, i.e. the power distribution. The program controls the power flows
such that the demand is always met, optimally withdrawing energy from the
storage systems.

Figure 7: Solution of the optimal control problem: power balance of the energy
system.

##  IX Conclusions

In this paper, we present models of energy production from wind and solar and
energy storage via a battery, thermal storage, and electrolysis. Such an
energy system can be operated optimally using the optimal control program that
we formulate in section VII. This program can be used as the regulator in an
MPC algorithm. The program smartly charges the storage systems before a low
power production is expected, and discharges them to satisfy the demand. The
excess power is sold when possible, and when prices are high. We demonstrate
that the power demand can be successfully satisfied, while maximizing the
economic revenue. Incorporating multiple and different ESSs is beneficial to
renewable energy systems, providing flexibility and reliability. Moreover,
they increase the economic value of the system, by eliminating or minimizing
energy curtailment, storing the surplus of energy.

## References

  * [1] D. Braga, “Integration of energy storage systems into the power system for energy transition towards 100 % renewable energy sources,” in 2021 10th International Conference on Energy and Environment (CIEM), 2021.
  * [2] O. Krishan and S. Suhag, “An updated review of energy storage systems: Classification and applications in distributed generation power systems incorporating renewable energy resources,” International Journal of Energy Research, vol. 43, no. 12, pp. 6171–6210, 2019.
  * [3] G. Sakas, A. Ibáñez-Rioja, V. Ruuskanen, A. Kosonen, J. Ahola, and O. Bergmann, “Dynamic energy and mass balance model for an industrial alkaline water electrolyzer plant process,” International Journal of Hydrogen Energy, vol. 47, no. 7, pp. 4328–4345, 2022.
  * [4] Øystein Ulleberg, “Modeling of advanced alkaline electrolyzers: a system simulation approach,” International Journal of Hydrogen Energy, vol. 28, no. 1, pp. 21–33, 2003.
  * [5] A. Belloni, L. Piroddi, and M. Prandini, “A stochastic optimal control solution to the energy management of a microgrid with storage and renewables,” in 2016 American Control Conference (ACC), pp. 2340–2345, 2016.
  * [6] A. Khatamianfar, M. Khalid, A. V. Savkin, and V. G. Agelidis, “Wind power dispatch control with battery energy storage using model predictive control,” in 2012 IEEE International Conference on Control Applications, pp. 733–738, 2012.
  * [7] R. Malkowski, P. Bućko, M. Jaskólski, W. Pawlicki, and A. Stoltmann, “Simulation of the dynamics of renewable energy sources with energy storage systems,” in 2018 15th International Conference on the European Energy Market (EEM), 2018.
  * [8] S. Teleke, M. E. Baran, A. Q. Huang, S. Bhattacharya, and L. Anderson, “Control strategies for battery energy storage for wind farm dispatching,” IEEE Transactions on Energy Conversion, vol. 24, no. 3, pp. 725–732, 2009\.
  * [9] M. Trifkovic, W. Marvin, M. Sheikhzadehy, and P. Daoutidis, “Dynamic real-time optimization and control of a hybrid energy system,” in 2013 European Control Conference (ECC), pp. 2669–2674, 2013.
  * [10] I. Aldaouab, M. Daniels, and R. Ordóñez, “Model predictive control energy dispatch to optimize renewable penetration for a microgrid with battery and thermal storage,” in 2018 IEEE Texas Power and Energy Conference (TPEC), 2018.
  * [11] M. N. Soltani, T. Knudsen, M. Svenstrup, R. Wisniewski, P. Brath, R. Ortega, and K. Johnson, “Estimation of rotor effective wind speed: A comparison,” IEEE Transactions on Control Systems Technology, vol. 21, no. 4, pp. 1155–1167, 2013.
  * [12] J. Jonkman, S. Butterfield, W. Musial, and G. Scott, “Definition of a 5-MW reference wind turbine for offshore system development,” tech. rep., 2009.
  * [13] G. W. Paltridge and C. M. R. Platt, Radiative processes in meteorology and climatology.  Elsevier, 1976.
  * [14] H. Madsen, H. Spliid, and P. Thyregod, “Markov models in discrete and continuous time for hourly observations of cloud cover,” Journal of Applied Meteorology and Climatology, vol. 24, no. 7, pp. 629–639, 1985.
  * [15] C. A. Thilker, H. Madsen, and J. B. Jørgensen, “Advanced forecasting and disturbance modelling for model predictive control of smart systems,” Applied Energy, vol. 292, p. 116889, 2021.
  * [16] R. Taesler and C. Andersson, “A method for solar radiation computations using routine meteorological observations,” Energy and Buildings, vol. 7, no. 4, pp. 341–352, 1984.
  * [17] P. Bacher, H. Madsen, and H. A. Nielsen, “Online short-term solar power forecasting,” Solar Energy, vol. 83, no. 10, pp. 1772–1783, 2009.
  * [18] H. Madsen, Statistically Determined Dynamical Models for Climate Processes.  PhD thesis, Technical University of Denmark, 1985.
  * [19] H. Madsen, J. Holst, and P. Thyregod, “A method for using hourly predictions of weather observations in optimal control of heat supply to buildings,” in International Conference on Building Energy Management, (Lausanne), pp. 380–389, 1987.
  * [20] H. Madsen and P. Thyregod, “Modelling the time correlation in hourly observations of direct radiation in clear skies,” Energy and Buildings, vol. 11, no. 11, pp. 201–211, 1988.
  * [21] A. Caraballo, S. Galán-Casado, A. Caballero, and S. Serena, “Molten salts for sensible thermal energy storage: A review and an energy performance analysis,” Energies, vol. 14, no. 4, 2021.
  * [22] J. Koponen, Energy efficient hydrogen production by water electrolysis.  PhD thesis, Lappeenranta-Lahti University of Technology LUT, 2020.
  * [23] J. A. E. Andersson, J. Gillis, G. Horn, J. B. Rawlings, and M. Diehl, “CasADi – A software framework for nonlinear optimization and optimal control,” Mathematical Programming Computation, vol. 11, no. 1, pp. 1–36, 2019.