# The Ancient Globular Clusters of NGC 1291

Kyle Hixenbaugh  Rupali Chandar  Ritter Astrophysical Research Center, The
University of Toledo, Toledo, OH 43606, USA  Angus Mok  OCAD University,
Toronto, Ontario, M5T 1W1, Canada

###### Abstract

We present a new catalog of 81 ancient globular clusters (GCs) in the early-
type spiral (SB0/a) galaxy NGC 1291. Candidates have been selected from
$B$,$V$, and $I$ band images taken with the Hubble Space Telescope, which also
reveal 17 younger ($\tau{\>\hbox to0.0pt{\raise 2.0pt\hbox{$<$}\hss}\lower
3.0pt\hbox{$\sim$}\>}\mbox{few}\times 100$ Myr) clusters. The luminosity
function shows a peaked shape similar to that found for GC systems in other
spiral and elliptical galaxies. The ancient clusters have a bimodal color
distribution, with approximately 65% (35%) of the population having blue (red)
colors. The red, presumably metal-rich GCs are more centrally concentrated, as
expected for a bulge population; while the blue, presumably metal-poor GCs,
are more broadly distributed, consistent with expectations of a halo
population. The specific frequency of GCs in NGC 1291 is higher than found
previously in most spiral galaxies. However, if we consider just the blue
subpopulation, we find $S_{\rm N,blue}=0.50\pm 0.06$, quite similar to that
found for other spirals. This result supports the hypothesis of a universal
population of halo GCs in spirals. The fraction of red GCs in NGC 1291, when
compared with those found in other galaxies, suggests that these correlate
with host galaxy type rather than with host galaxy luminosity.

††software: astropy (The Astropy Collaboration 2013, 2018), Baolab(Larsen
1999)

##  1 Introduction

Globular clusters (GCs) are some of the oldest known objects in the universe,
and as such give clues to the early stages of galaxy formation. These ancient
($\tau{\>\hbox to0.0pt{\raise 2.0pt\hbox{$>$}\hss}\lower
3.0pt\hbox{$\sim$}\>}10$ Gyr), massive ($M\approx 10^{4}-10^{6}~{}M_{\odot}$)
stellar systems (Brodie & Strader, 2006) are found in the halos of nearly all
galaxies, at least those with stellar masses greater than $\sim
10^{9}~{}M_{\odot}$ (Harris et al., 2017). In earlier-type galaxies, some
globular clusters are believed to be associated with bulges (Zinn, 1985), and
in some lower mass spiral and irregular galaxies they have also been shown to
have disk kinematics (Olsen et al., 2004; Schommer et al., 1992).

The Hubble Space Telescope (HST) has enabled a number of major breakthroughs
in our understanding of globular cluster systems, primarily in elliptical
galaxies, which have fairly simple morphologies and only contain ancient
clusters (Peng et al., 2006, 2008; Jordán et al., 2007). Spiral galaxies, by
contrast, have more complicated structures due to on-going star-formation in
their disks, and can have rich populations of young clusters, (Portegies Zwart
et al., 2010; Christian & Schommer, 1988; Chandar et al., 1999, 2001), making
it more challenging (but not impossible) to identify ancient globular clusters
(Chandar et al., 2004). Outside of the Milky Way, M31 (Huxor et al., 2011,
2014) and M81 (Nantais & Huchra, 2010; Nantais et al., 2011; Santiago-Cortés
et al., 2010) are the spiral galaxies with the best-studied globular cluster
systems. Other spirals which have more recent catalogs of globular clusters
from the Hubble Space Telescope include: M101, M51, NGC 4258, and NGC 6278
(Lomelí-Núñez et al., 2022).

Lenticular and early-type spirals provide an important transition between non-
star forming ellipticals and later-type galaxies with significant current star
formation. Lenticular and early type galaxies have bimodal GC populations: one
blue and metal-poor, and the other red and metal-rich, whose colors correlate
with the host galaxy luminosity and color (Peng et al., 2006). This has been
studied in galaxies such as the lenticular galaxy NGC 1172 (Ennis et al.,
2020) and the early-type NGC 6876 (Ennis et al., 2019). These GC populations
have distinct spatial distributions, with redder, metal rich GC’s more
centrally concentrated than their metal-poor counterparts. The distinct
subpopulations, along with with high specific frequencies (the number of GCs
normalized to a galaxy with $M_{V}=-15$), are believed to result from earlier
merger events (Ennis et al., 2020, 2019).

In this work, we present the first study of the globular cluster system of NGC
1291, an early-type barred spiral (SB0/a) galaxy that is forming stars slowly
in its disk (Thilker et al., 2007), and has a prominent bulge (de Vaucouleurs,
1975). The galaxy presents nearly face-on (Bosma et al., 2010). In this work,
we adopt a distance of $8.9$ Mpc ($m-M=29.75)$ (de Vaucouleurs, 1975), based
on distance modulus estimations from bright stars, which is quite similar to
the estimate of $9.08\pm 0.29$ Mpc made by McQuinn et al. (2017) based on the
tip of the red giant branch method.

We use archival, optical images of four pointings within NGC 1291 taken with
$HST$ to identify and study it’s globular cluster system. It should be noted
the globular clusters identified in this work are candidates and have no
velocity confirmation. For simplicity, we refer to these candidates as
globular clusters throughout this paper. The rest of this paper is organized
as follows. In Section 2 we present how we select our globular cluster sample,
including photometry and an estimate of completeness. Section 3 presents the
color-magnitude and color-color distributions of the globular clusters, along
with estimates of their masses. Section 4 focuses on the luminosity/mass
functions, color/metallicity distributions, and spatial distributions of the
ancient clusters. We present and discuss the specific frequency of the GCs in
Section 5, and summarize our main results in Section 6.

##  2 Observations, Cluster Selection, and Photometry

###  2.1 Observations and Photometry

We use 4 pointings in NGC 1291 taken by the HST ACS/WFC camera in 2005 April
and 2013 September. All pointings have images in three filters: F435W ($B$
band), F555W or F606W ($V$ band), and F814 ($I$ band). Each field covers an
area of approximately 40,000 $\arcsec^{2}$ or $74~{}\mbox{kpc}^{2}$. At the
adopted distance of 8.9 Mpc to NGC 1291, the fields used cover a total area of
around $296~{}\mbox{kpc}^{2}$, and the pixel scale is 43 pc arcsec-1. Figure 1
shows the locations of the 4 pointings, and Table 1 compiles basic information
for each one, including position.

Figure 1: Ground-based optical image of NGC 1291 showing the locations of the
four $HST$ ACS/WFC pointings (labeled) used in this work. A $5\arcmin$ scale
bar is shown for reference.

Reduced images were downloaded from the Hubble Legacy Archive
(HLA).111https://hla.stsci.edu/hlaview.html These images have been reduced
through the HLA pipeline, which includes bias subtraction, dark subtraction,
flat-fielding, cosmic ray rejection, and drizzling.

We detect sources using the DAOfind routine in the Astropy Python package. The
detection algorithm is run on each $V$ band image, set to approximately a
$3\sigma$ threshold. A visual inspection confirmed that essentially all
obvious, point-like sources are detected. The list of detected objects
includes star clusters, individual stars, and background galaxies. Masks are
applied to bright, saturated foreground stars to prevent their detection.
Aperture photometry was performed in each available filter using the photutils
python package in Astropy. The brightness of these sources is measured using a
circular aperture with radius of 3 pixels, with the background level estimated
in an annulus with radii between 5 and 8 pixels. We apply the following
zeropoints (obtained from the STScI webpage) to our aperture photometry to put
it on the VEGAMAG system: $F435W=25.788$, $F555W=25.732$, $F606W=26.406$, and
$F814W=25.528$. Table 2 compiles the magnitudes, positions, and size
information (described in Section 3.2) for five clusters in our catalog. The
full table has released with this publication. The VEGAMAG photometric system
is similar to the Johnson-Cousins system for the ACS/WFC filters used in this
work, with only small differences in the zeropoints on the order of $\sim
0.03-0.05$ mag (Sirianni et al., 2005).

In order to determine aperture corrections, we measure the curve of growth out
to 20 pixels for ten isolated clusters in the $V$ band images. We find an
average aperture correction (between 3 and 20 pixels) of $-0.47$ mag, and
include an additional $-0.1$ mag to correct from 20 pixels to infinity, giving
a total $V$-band aperture correction of $-0.57$ mag. This average correction
was applied to each cluster (selection described below). In addition, PT3 in
Figure 1 is observed in the F606W rather than the F555W filter. To account for
this, we find 15 bright and isolated sources in the region of overlap between
PT3 and PT2. We calculate the magnitude difference for the sources in common
between the two filters, and find an average difference of $\sim 0.258$ mag.
We add this to offset to the F606W photometry to convert to the F555W band.

In the rest of this work, we use magnitudes that have been corrected for the
following foreground extinction values: $m_{\rm F435W}=0.047$ mag, $m_{\rm
F555W}=0.036$ mag, and $m_{\rm F814W}=0.02$ mag, which gives small shifts in
the colors of $m_{\rm F435W}-m_{\rm F555W}=0.011$ mag and $m_{\rm
F555W}-m_{\rm F814W}=0.016$ mag Schlafly & Finkbeiner (2011).222Extinction
values can be found with NED (https://ned.ipac.caltech.edu), along with the
paper producing the values, Schlafly & Finkbeiner (2011).

###  2.2 Cluster Selection

The initial source list contains many objects, including, clusters and bright,
individual stars in NGC 1291, background galaxies, and foreground stars. The
goal of this work is to select and study the ancient globular clusters in NGC
1291. To remove contaminants and select candidate ancient clusters, we use
automated selection criteria followed by a final visual inspection. Given the
relative proximity of NGC 1291, essentially all globular clusters are broader
than the PSF (see Section 3.2 for size measurements. We show some example
clusters from our dataset in Figure 2.

Figure 2: Three color images of typical candidate clusters, showing the range
of colors from reddest (top row), to intermediate (middle row), to bluest
(bottom row) of the detected clusters. Each panel covers an approximate area
of $~{}18"^{2}$, or $78~{}\mbox{pc}^{2}$ , and gives the corresponding ID
number from our catalog.

We build training sets of $\sim 15$ hand-selected stars, clusters, and
galaxies to help guide the automated cuts. Following Whitmore et al. (2021),
we identify our training set using a variety of tools, including color images
and plotting radial profiles using the imexam python package, where stars all
have similar sharp radial profiles but clusters are clearly more extended. Our
primary quantitative tool is the Concentration Index (CI), which is the
difference in the V-magnitude measured in apertures of 0.5 and 3 pixel radii.
Clusters have a larger measured CI than stars. In Figure 3, we plot the CI vs.
V magnitude for our training set of stars (star symbols) and clusters (red
circles). The clusters clearly have larger values of CI at similar magnitudes.
Based on the measurements of the sources in our training set, we assume that
sources with CI$>2.65$ are candidate clusters.

Figure 3: Plot of concentration index (CI) vs apparent magnitude (both
measured in F555W) for the training set of stars and clusters. There is a
clear separation between the two objects. The black lines represent the
CI$=2.65$ used to automatically select cluster candidates.

We eliminate all sources with $CI<2.65$, then visually inspect all remaining
sources. Background galaxies and close stellar pairs or triples, which also
have larger values of CI, were eliminated by visual inspection. This procedure
resulted in 98 candidate clusters; we show representative color images for
these sources in Figure 2. Three clusters representing the range of cluster
colors (which we describe in the next section), were selected to show the
range of cluster appearances. While we believe that our catalog of GCs is
fairly contamination free, almost all cluster catalogs selected from imaging
alone have at least some low-level of contamination from spurious sources, and
future radial velocity measurements could be used to eliminate any
contaminants.

###  2.3 Completeness

Figure 4: Completeness fraction for artificial sources recovered in the F555W
band. Details given in text.

We assess the completeness of our cluster catalogs by creating artificial
clusters and adding them to the F555W band images. We use the Baolab task
MKCMPPSF to create 4000 artificial, extended sources with an assumed King30
profile (Larsen, 1999). These artificial clusters are added randomly to the
images using the MKSYNTH task, and have a range of F555W band magnitudes
between 18 and 25 mag. For the artificial clusters, we adopt the approximate
median size of 0.5 pixels for the detected cluster population, which we found
by running Baolab on our cluser sample.

After the artificial clusters are created and added to the images, we rerun
our detection and selection algorithms. We then apply the same CI cuts that
were used on the actual data. The fraction of artificial clusters of different
magnitudes that are recovered by this process are plotted in Figure 4. The
recovery fraction is quite high, above 90%, down to a magnitude of $\approx
23$. The completeness fraction is $\approx 50$% near $m_{F555W}\sim 24$ and
drops rapidly below this value.

##  3 Cluster Properties

###  3.1 Luminosities and Colors

Figure 5: F435W-F555W vs F555W color-magnitude diagram of all candidate
clusters (red circles). Milky Way globular clusters (grey squares) from the
Harris (1996a) catalog are plotted for reference with the uncertainty due to
the conversion between the HST VEGAMAG and Johnson-Cousins BVI photometric
systems. Magnitudes reported are absolute V band magnitudes. The blue and red
lines represent the typical peak B-V color for blue, metal-poor and red,
metal-rich globular clusters, respectively, for the Milky way and early-type
galaxies.

Figure 5 shows the ($F435W-F555W$) versus F555W band color magnitude diagrams
(CMD) of our detected cluster sample as the red circles. One of the most
striking features in this figure is the contrast between faint, bluer clusters
(with F435W-F555W${\>\hbox to0.0pt{\raise 2.0pt\hbox{$<$}\hss}\lower
3.0pt\hbox{$\sim$}\>}0.5$) versus the much brighter, redder clusters
(F435W-F555W$>0.5$). These bluer clusters are almost certainly young
(${\>\hbox to0.0pt{\raise 2.0pt\hbox{$<$}\hss}\lower
3.0pt\hbox{$\sim$}\>}\mbox{few}\times 100$ Myr) clusters in the disk of NGC
1291, with fairly low masses. A small population of young, low mass clusters
is consistent with what we would expect from a galaxy that is currently
forming stars at a low rate. The majority of the clusters are redder, with
$F435W-F555W$ colors between $\sim 0.5$ and 1.0, and have magnitudes that
reach up to $M_{F555W}\sim-10.5$ mag; these are candidate globular clusters in
NGC 1291. The colors and magnitudes of globular clusters in the Milky Way
(gray circles; from Harris (1996b)), shown for comparison, nicely overlap the
colors and magnitudes of the redder clusters in NGC 1291.

In Figure 5 we also show the mean ($B-V$) colors of the two color peaks found
for GC systems in the Milky Way based on an updated version of the Harris
(1996b) catalog and a number of elliptical and lenticular galaxies, with
typical values of 0.65 (blue, metal-poor) and 0.94 (red, metal-rich; (Harris,
1996b; Kissler-Patig et al., 1997)). The globular clusters in NGC 1291 appear
to follow this bimodal trend in color. Because colors are sensitive to
metallicity in single stellar populations older than a few Gyr, the color
distribution suggests that the GCs in NGC 1291 have a bimodal distribution.
The color distribution is explored further in Section 4.2.

Figure 6: F435W-F555W vs F555W-F814W color-color plot. Predictions from two
different cluster evolution models from Bruzual & Charlot (2003) are shown,
Z=0.05 (solar metallicity) as the black solid line and Z=0.004
($\frac{1}{4}\times$ solar) as the dashed purple line. The predicted location
of clusters with ages of $10^{7}$, $10^{8}$, $10^{9}$, and $10^{10}$ years are
indicated. Clusters in NGC 1291 follow the models closely, and are color coded
by age: young clusters (blue), metal-poor globular clusters (green), metal-
rich globular clusters (red). The colors of Milky Way GCs are also shown for
reference (gray), with the uncertainty due to the conversion between the HST
VEGAMAG and Johnson-Cousins BVI photometric systems. The VEGAMAG colors of
bright cluster in M51 (shown in yellow) show that the blue clusters in NGC
1291 have similar colors as the young clusters formed in the disk of M51.

In Figure 6 we plot the $F435W-F555W$ vs. $F555W-F814W$ color-color diagram of
clusters in NGC 1291. These are compared with the predictions from stellar
evolution models from Bruzual & Charlot (2003) for two different
metallicities, Z=0.05 (solar), shown as a solid line, and Z=0.004
($\frac{1}{4}\times$ solar), shown as a purple dashed line. These models start
at 1 Myr (upper left) and go through ages older than 10 Gyr (lower-right), as
indicated. Predictions for the evolution at the different metallicities is
fairly similar, except at older ages where metal-poor clusters have bluer
colors than metal-rich ones at the same age. This is known as the age-
metallicity degeneracy, and affects clusters with ages $\tau{\>\hbox
to0.0pt{\raise 2.0pt\hbox{$>$}\hss}\lower 3.0pt\hbox{$\sim$}\>}1$ Gyr.

The color-color diagram shows that there are a handful of young clusters in
NGC 1291 (blue points), with ages ranging from a few to a few hundred million
years. These have similar colors to bright young clusters in M51 taken from
Chandar et al. (2016). However, as already seen in the color-magnitude
diagram, the majority of clusters in NGC 1291 are quite red, and overlap with
the colors of Milky Way globular clusters (shown as the gray squares in Figure
6). We assume that these are ancient clusters, with ages ${\>\hbox
to0.0pt{\raise 2.0pt\hbox{$>$}\hss}\lower 3.0pt\hbox{$\sim$}\>}10$ Gyr. The
distribution of colors for these ancient GCs in NGC 1291 appears to be
bimodal. To approximately define the 2 different regimes, we use two
perpendicular bisectors to the stellar models, and indicate the clusters that
fall into the bluer, presumably metal-poor clump (shown as green circles), and
those that fall into the red, presumably metal-rich clump (red circles).

###  3.2 Size Estimates

We estimate the size of each cluster using 2 different methods. The first
method is the concentration index (CI), which is the difference in F555W-band
magnitude measured in a 0.5 pixel and a 3 pixel aperture. The CI is a very
simple measure of size that works well for most clusters. CI is generally
preferred to estimate the sizes of clusters in crowded regions because the
largest aperture is only 3 pixels. Based on the CI values measured for stars
and clusters in our training sets, we adopted a minimum size of $CI=2.65$ to
select candidate clusters. We find that the maximum size for a cluster in our
sample is $CI=3.38$, the mean value is 2.86 and the median is 2.82.

The effective radius ($R_{eff}$) of each cluster is estimated from the Ishape
program in the Baolab software Larsen (1999). This is the radius that contains
half the light. Ishape convolves analytic profiles for different $R_{eff}$,
which model the profile of the surface brightness of a cluster, with the PSF
to determine the best-fit size for each cluster. A subsampled PSF (as required
by Baolab) was generated from $\sim 30$ bright, fairly isolated point sources.
These were selected from different locations within an ACS $F555W$-band image,
since this filter is the most sensitive and hence stars tend to have the
highest signal-to-noise in this filter. To run Ishape, a King 30 model is
assumed for each cluster King (1962), although the resulting $R_{eff}$ is not
sensitive to the exact profile shape that we select. Ishape returns a best-fit
FWHM value broader than the PSF in pixels, which we convert to $R_{eff}$ in
parsecs using the pixel scale of $0.05\arcsec\mbox{pix}^{-1}$, assuming that
$1\arcsec$ subtends 43 pc at the adopted distance of 8.9 Mpc, and including a
factor of 1.48 as stated in the Ishape manual. The Ishape measurements confirm
that all of our globular cluster candidates are broader than the PSF. We find
a mean FWHM ($R_{eff}$) of 1.03 pix (3.28 pc) and a median of 0.77 pix (2.46
pc), with a range between 0.31 pix (0.99 pc) and 5.0 pix (15.96 pc). These are
similar to values found for the Milky Way GC system, which has a median
$R_{eff}$ around 2.98 pc, and spiral galaxy M101, whose bright GCs have a
median of 2.41 pc Simanton et al. (2015). The full distribution of sizes is
shown in Figure 7.

Figure 7: The distribution of measured values for $R_{eff}$ in pc for our GC
sample. The median size for Milky Way GC’s is overplotted as a black line.

We compare the relationship between CI and FWHM in Figure 8 and find a strong
correlation between the 2 different size measurements, although the
relationship does not appear to be strictly linear (i.e. the CI values appear
to ”saturate”). There are four obvious outliers in this figure (circled). We
find these are all in crowded regions, and therefore their FWHM values are
more likely to be suspect than CI.

Figure 8: The measured values for CI vs FWHM from Ishape show good general
agreement for the sizes of the clusters. There are four obvious outliers
(circled) that all come from crowded regions.

##  4 Results

In the rest of this paper we will focus on the ancient globular clusters
(shown in green and red in Figure 6).

###  4.1 Luminosity Function

The luminosity (mass) function of the globular cluster systems in the Milky
Way, M31, and elliptical galaxies have a distinctive shape, where they
initially rise towards lower luminosities and masses333The luminosities of
ancient globular clusters are a good proxy for their masses, although there
may be some metallicity-dependence to the mass-to-light ratio., then flatten
near $M_{V}\sim-7.4$ mag ($\sim 10^{5}~{}M_{\odot}$) (Kundu & Whitmore, 2001);
this peak magnitude is similar in the Johnson-Cousins and VEGAMAG systems.
This shape is different from the one observed for young cluster populations
($<\mbox{few}\times 100$ Myr), which continue to rise in power-law fashion at
lower luminosities and masses (Zhang & Fall, 1999). More traditionally, when
plotted as a histogram in magnitudes, GC lumniosity functions can be modeled
by a Gaussian distribution that peaks near $M_{V}\approx-7.4$ mag. This
distribution corresponds to a lognormal distibrution of their masses (Jordán
et al., 2007). The peaked distribution for old globular clusters is believed
to result from the dynamical evolution of the population, where lower mass
clusters are disrupted faster than their more massive counterparts, at a given
stellar density (McLaughlin, 2000; Fall & Zhang, 2001; Zhang & Fall, 1999).

The globular cluster luminosity function in NGC 1291 is shown in Figure 9. We
see a peak in our CG luminosity function around a magnitude off 22.5, which
corresponds to $M_{V}\approx-7.4$, similar to the peak found in the Milky Way
and many other galaxies. This peak occurs above the 90% completeness level of
our sample, so is robust against incompleteness.

Figure 9: The luminosity function for globular cluster candidates in NGC1291
is shown in blue. The red dashed lines shows the completeness corrected
lumniosity function derived by boosting counts based on the recovery fraction
detailed in Section 2.3.

###  4.2 Color/Metallicity Distributions of Globular Clusters

The metallicity distribution of GC systems sheds light on the formation
history of the parent galaxy. For the GC systems in elliptical galaxies,
widespread bimodality in the color (and by extension metallicity)
distributions suggests that these galaxies experienced multiple epochs of
cluster formation early in their histories, possibly due to mergers (Kundu &
Whitmore, 2001). However, much less is known concerning the metallicity
distributions of GC systems in spirals. The two best-studied spirals, the
Milky Way Galaxy and M31, both have bimodal GC metallicity distributions
(Côté, 1999; Perrett et al., 2002; Huxor et al., 2014).

The colors of globular clusters are a proxy for their metallicity. The
intrinsic metallicity distribution of GCs in the Milky Way is bimodal (Côté,
1999; Perrett et al., 2002). This manifests as an extended color distribution
in Figure 6 (gray dots). A similar extent in color is seen in the ancient
cluster population of NGC 1291, as highlighted by our color coding bluer
(green) vs. redder globular clusters (red). Our clusters have mean and median
F435W-F555W colors of 0.79 and 0.73, and 1.04 and 1.00 respectively, in
F555W-F814W. As mentioned in Section 3, these colors align closely with colors
of GC’s in other early-type spiral galaxies and in the Milky Way.

We develop a ”two color index” for globular clusters by measuring their
distance in color-color space (shown in Figure 6) along the $Z=0.02$ model
track, starting from the perpendicular model bisector at younger ages (the
results are nearly identical if we use the lower metallicity model instead).
The distance of each globular cluster along the model as measured from this
zero-point constitutes the ”two-color index”. The top panel of Figure 10 plots
the distribution of this index for the candidate globular clusters in NGC
1291, with consistent colors between this panel and Figure 6. For comparison,
we also show the distributions of single colors, with (F555W-F814W) in the
middle panel and (F435W-F555W) in the bottom. The two-color index in the top
panel shows the bimodality in globular cluster colors more clearly than either
of the single color indices.

Figure 10: Histograms showing the distribution of $F555W-F814W$ (middle) and
$F435W-F555W$ (bottom) colors. The top panel shows clear bimodality in the
”two-color” index, where we measure the distance of each cluster along the
model in $F435W-F555W$ vs $F555W-F814W$ color-color space (starting from the
arbitrary zero point shown as the higher black line in Figure 6). Two
Gaussians were fit to the two-color distribution.

The two-color index for globular clusters in NGC 1291 is well fit by a bimodal
distribution, but not a unimodal one. We therefore fit two Gaussians to the
distribution shown in the top panel of Figure 10, using the curve_fit function
in the python package Scipy. The best fit values of each peak, height, and
width $\sigma$ were returned. The blue, metal-poor (green) clusters have a
mean color index located at 0.29, a peak of 13.14, and a FWHM of 0.19. The
red, metal-rich (red) clusters have a mean color index of 0.63, a lower peak
of 5.98, and a FWHM of 0.19.

###  4.3 Spatial Distributions

Figure 12 displays the location of each cluster in the galaxy on a ground
based image obtained from the DSS survey 444(https://archive.stsci.edu). The
young, blue clusters appear to be located in the disk (see more below), while
the red, presumably metal-rich GCs appear somewhat more centrally concentrated
than the green, metal-poor GCs.

We show histogram distributions of the radial locations of the young (blue),
metal-poor GC (green), and metal-rich GC (red) populations in Figure 11. This
figure shows that the young, blue clusters avoid the central region,
consistent with the expectations that they reside in the disk. The green and
red (metal-poor and metal-rich ancient GCs) are more centrally concentrated
than the young clusters, consistent with expectations of spheroidal
populations relative to those in galaxy disks. While there are some gaps in
the numbers because of incomplete coverage, the red clusters tend to be more
evenly distributed across galactocentric distances than the green (and the
coverage is identical for both populations). This is consistent with metal-
rich (red) GCs being somewhat more centrally concentrated than the metal-poor
(green) counterparts. This follows the expectations of a population of metal-
rich GCs associated with the bulge and another of metal-poor GCs associated
with the halo in NGC 1291.

Figure 11: A histogram of the distances in kpc for young cluster (blue with
horizontal hatches), metal poor GCs (green with diagonal hatches), and metal
rich GCs (red with vertical hatches). Distances are calculated by projecting
the face-on images onto the xy plane and finding the distance from the center
of the galaxy.  Figure 12: The locations of all selected clusters in NGC 1291
are shown on a ground-based image with the following color-coding: blue
represents young, ${\>\hbox to0.0pt{\raise 2.0pt\hbox{$<$}\hss}\lower
3.0pt\hbox{$\sim$}\>}\mbox{few}\times 100$ Myr clusters; green represents
ancient metal-poor globular clusters, and red represents ancient metal-rich
globular clusters. The orange boxes represent our four pointings, and the
inner (outer) ellipse shown in red (blue) represent our defined bulge (disk)
regions, respectively. The red circle represents the bulge region and reaches
out to $\sim 5.5$ kpc from the center of the galaxy, and the disk (shown by
the blue circle) out to 13.5kpc.

##  5 The Globular Cluster Specific Frequency

The specific frequency ($S_{N}$) allows us to compare the number of globular
clusters between galaxies by normalizing the population to a galaxy of
luminosity $M_{V}=-15$. $S_{N}$ is defined as:

$S_{N}=N_{GC}*10^{0.4*(M_{V}+15)}$ (1)

where $N_{GC}$ is the total number of globular clusters and $M_{V}$ is the
absolute $V$-band magnitude of the galaxy (Harris & van den Bergh, 1981). An
additional parameter, $T_{L}$ (specific mass), is sometimes used because
spiral galaxies can be dominated by different age stellar populations, so
their luminosities do not perfectly track their masses, which is the more
fundamental parameter. $T_{L}$ is defined as the number of GC’s per unit
stellar mass ($10^{9}M_{\odot}$) (Goudfrooij et al., 2003). Because there are
fairly large uncertainties in the appropriate mass-to-light ratio to adopt for
spiral galaxies, in this work we focus on specific frequency.

To calculate specific frequency and specific mass, we need the total number of
clusters in NGC 1291. Along with the 81 detected globular clusters, we must
account for incompleteness in our number from not including faint objects and
due to incomplete coverage of the galaxy, in order to estimate the total
number of globular clusters in NGC 1291.

As discussed in Section 2.3, we have a high recovery fraction, above 90%, down
to $m_{F555W}\approx 23$, and in Section 4, we found that the shape of the
luminosity function is similar to that of GCs in the Milky Way and other
galaxies. The turnover magnitude for the Milky Way, and other early type
spirals is $\approx-7.4$ (Chandar et al., 2004). At the distance modulus of
29.75 for NGC 1291, this corresponds to an apparent magnitude of 22.35. Since
we have such high completeness brighter than the turnover magnitude, we do not
apply any completeness correction for potentially missing clusters at
magnitudes down to the turnover.

To correct for incomplete spatial coverage of NGC 1291, we make corrections
for the bulge, disk, and outer halo regions (shown in Figure 12) separately.

  * •

Bulge region: we find 22 globular clusters brighter than the turnover
magnitude of $M_{V}\approx-7.4$ mag in the bulge, and our images cover
$\approx 73\%$ of this region. We multiply by a factor of 1.27 to account for
the incomplete coverage, resulting in 28 clusters brighter than the turnover.
The bulge area is defined by the red ellipse in Figure 12.

  * •

Disk region: we repeat this procedure for the disk region (area between the
blue and red ellipses in Figure 12), using a factor of 1.52 to correct for
incomplete spatial coverage of the disk, and estimate a total of 64 globular
clusters brighter than the turnover in this region.

  * •

Outer halo region: our images cover out to a galactocentric distance of
$R_{gc}=13.5$ kpc. In the Milky Way, we find only 20 globular clusters
brighter than the turnover outside of this distance. We assume that the outer
halo globular clusters in NGC 1291 have a radial distribution that is similar
to that in the Milky Way.

  * •

We estimate a total of 112 globular clusters brighter than $M_{V}=-7.4$ by
summing up the contributions from the coverage-corrected bulge, disk, and halo
regions. We then multiply this number by a factor of two to account for the
clusters fainter than the turnover magnitude, resulting in an estimated total
number of 224. We assume Poisson statistics, for an estimated total globular
cluster population of $224\pm 15$ in NGC 1291.

For NGC 1291 we define intrinsic properties to aid in calculations of specific
frequency and mass. We adopt an absolute V band magnitude of $M_{V}=-21.05$
(estimated from the total apparent $V$ band magnitude given in HyperLeda and
corrected for the distance modulus of 29.75 assumed here). We adopt a mass of
$5.8\pm{2.9}\times 10^{10}M_{\odot}$ (Bittner et al., 2021) for NGC 1291,
which includes an uncertainty of 50% on the assumed $M/L_{V}$. By using Eq 1
and the stated parameters for NGC 1291, we find a specific frequency of
$S_{N}=0.84\pm 0.06$ and a $T_{L}=3.79\pm 1.90$. If we consider the two
different color populations separately, we find that the blue (red) GCs
comprise 65% (35%) of the observed sample. If we assume they comprise a
similar fraction of the total population, we estimate $S_{\rm N,blue}=0.55\pm
0.06$. If instead we estimate the total number of blue GCs by correcting for
incomplete spatial coverage as described above, we find $S_{N}=0.51\pm 0.06$.
Therefore, regardless of the specific method used to estimate the number of
blue GCs, we find values of $S_{\rm N,blue}$ close to 0.5.

Previous studies of the GC systems in a handful of spiral galaxies have
correlated their properties with galaxy type, luminosity, and mass. One key
finding was that later-type spirals (Sb and later) appear to have a fairly
constant $S_{N}=0.5\pm 0.2$ (and $T=1.3\pm 0.2$)(Goudfrooij et al., 2003;
Chandar et al., 2004). Because these galaxies are dominated by blue, metal-
poor GCs, this suggests that spirals may form a fairly ”universal” halo GC
population. Another trend that has been suggested but not established
definitively because of the poorer data quality of previous studies is that
$S_{N}$ and $T_{N}$ increase from late to early-type spirals, presumably
because earlier-type galaxies have formed an additional population of metal-
rich GCs. This trend appears to be unaffected by galaxy luminosity or mass.

In elliptical galaxies, the globular cluster specific frequency shows somewhat
different trends, although they are also dominated by the blue GCs.
Midluminosity ellipticals, which have brightnesses similar to NGC 1291, have
$S_{N,blue}\approx 0.5$, very similar to results for spirals of all
luminosities. However, brighter (giant) and fainter (dwarf) ellipticals have
significantly higher blue GC specific frequencies, reaching values as high as
5 for some dwarf ellipticals (Peng et al., 2008). The fraction of red GCs
increases with luminosity, ranging from $\approx 0.1$ to $\approx 0.5$, until
a luminosity of $\approx-22$. At this value, the fraction of red GC’s appears
to turn over (Peng et al., 2008).

As a very early-type spiral (SB0/a), the GC system in NGC 1291 presents an
important test of previous suggestions that spiral galaxies form universal
halo GC populations. When we only consider the blue GCs which are typically
associated with halos, we found $S_{N,blue}\approx 0.5$, right in line with
the values found for later-type spirals that only form halo populations. Our
results for NGC 1291 therefore supports the idea of a ”universal” relationship
between blue, metal-poor GCs and their host galaxy luminosities.

Red GCs provide clues to whether NGC 1291 better follows the trends found for
spiral or elliptical galaxies. There are only a few spiral galaxies with
estimate fractions of red GCs, which we compile below along with the galaxy
type and V-band magnitude.

NGC 1291: galaxy type$=$SB0/a, $M_{V}=-21.05$, $\sim 35$% red GCs (this work)

M81: galaxy type$=$Sab, $M_{V}=-21.63$ (Chandar et al., 2004), $\approx 40$%
red GCs (Chandar et al., 2004)

M31: galaxy type$=$Sb, $M_{V}=-21.8$ (Chandar et al., 2004), $\sim 34$% red
GCs (Barmby et al., 2000)

Milky Way: galaxy type$=$Sbc, $M_{V}\approx-21.3$ (Chandar et al., 2004),
$\sim 30$% red GCs (Harris, 1996a; Côté, 1999; Barmby et al., 2000; Rhode &
Zepf, 2003).

Spirals that are later-type than the Milky Way (Sc and later) appear to have
mostly blue, metal-poor GCs (e.g., Chandar et al. (2004)), and hence a
significantly lower fraction of red GCs than the galaxies listed above. The
available information for NGC 1291, M81, M31, and the Milky Way do not show
any obvious trend with galaxy luminosity, but when later-type spirals are
included there is a general trend that the fraction of red GCs increases
towards earlier spiral-type/bulge contribution. These trends (or lack thereof)
suggests that NGC 1291 more closely follows spiral GC systems than those
formed in ellipticals.

##  6 Summary and Conclusion

  * •

We have produced a new catalog of compact star clusters in the nearby ($D=8.9$
Mpc), early-type (SB0/a) spiral galaxy NGC 1291 based on $BVI$ images taken
with the Hubble Space Telescope. 81 of the clusters have colors expected of
ancient globular clusters, while 17 have colors expected of young
($\sim$few$\times$100 Myr) clusters. This is the first published cluster
catalog in NGC 1291.

  * •

The ancient globular clusters have a bimodal color distribution, with $\approx
65$% ($\approx 35$%) being blue (red). The red, presumably metal-rich clusters
are more centrally concentrated than the blue, presumably metal-poor, ones.

  * •

The correlation between color/metallicity and spatial distribution is
consistent with the expectations that the red GCs are associated with the
bulge, while the blue GCs are associated with the halo. The small number of
young clusters are clearly located in the disk of NGC 1291.

  * •

The total globular cluster specific frequency is higher than that found in
typical, later-type spirals. The value when only the blue globular clusters
are considered, $S_{N,\rm blue}=0.50\pm 0.06$, is quite similar to that for
later-type spirals. This is consistent with the expectations of a universal
halo population that scales with host galaxy luminosity.

  * •

The fraction of red GCs for NGC 1291, when compared with published results for
other spiral galaxies, shows no overall trend with galaxy luminosity but is
consistent with an increasing fraction from late- to early-type spirals. These
trends are dissimilar to those found for globular cluster systems in
elliptical galaxies.

We thank the anonymous referee for their comments that improved our paper.

## References

  * Barmby et al. (2000) Barmby, P., Huchra, J. P., Brodie, J. P., et al. 2000, AJ, 119, 727, doi: 10.1086/301213
  * Bittner et al. (2021) Bittner, A., de Lorenzo-Cáceres, A., Gadotti, D. A., et al. 2021, A&A, 646, A42, doi: 10.1051/0004-6361/202039505
  * Bosma et al. (2010) Bosma, A., Gadotti, D. A., de Blok, W. J. G., & Athanassoula, E. 2010, in Astronomical Society of the Pacific Conference Series, Vol. 421, Galaxies in Isolation: Exploring Nature Versus Nurture, ed. L. Verdes-Montenegro, A. Del Olmo, & J. Sulentic, 53
  * Brodie & Strader (2006) Brodie, J. P., & Strader, J. 2006, ARA&A, 44, 193, doi: 10.1146/annurev.astro.44.051905.092441
  * Bruzual & Charlot (2003) Bruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000, doi: 10.1046/j.1365-8711.2003.06897.x
  * Chandar et al. (1999) Chandar, R., Bianchi, L., & Ford, H. C. 1999, ApJS, 122, 431, doi: 10.1086/313225
  * Chandar et al. (2001) —. 2001, A&A, 366, 498, doi: 10.1051/0004-6361:20000263
  * Chandar et al. (2004) Chandar, R., Whitmore, B., & Lee, M. G. 2004, ApJ, 611, 220, doi: 10.1086/421934
  * Chandar et al. (2016) Chandar, R., Whitmore, B. C., Dinino, D., et al. 2016, ApJ, 824, 71, doi: 10.3847/0004-637X/824/2/71
  * Christian & Schommer (1988) Christian, C. A., & Schommer, R. A. 1988, AJ, 95, 704, doi: 10.1086/114668
  * Côté (1999) Côté, P. 1999, AJ, 118, 406, doi: 10.1086/300930
  * de Vaucouleurs (1975) de Vaucouleurs, G. 1975, ApJS, 29, 193, doi: 10.1086/190341
  * Ennis et al. (2019) Ennis, A. I., Bassino, L. P., Caso, J. P., & De Bórtoli, B. J. 2019, MNRAS, 488, 770, doi: 10.1093/mnras/stz1660
  * Ennis et al. (2020) Ennis, A. I., Caso, J. P., Bassino, L. P., Salinas, R., & Gómez, M. 2020, MNRAS, 499, 2554, doi: 10.1093/mnras/staa2967
  * Fall & Zhang (2001) Fall, S. M., & Zhang, Q. 2001, ApJ, 561, 751, doi: 10.1086/323358
  * Goudfrooij et al. (2003) Goudfrooij, P., Strader, J., Brenneman, L., et al. 2003, MNRAS, 343, 665, doi: 10.1046/j.1365-8711.2003.06706.x
  * Harris (1996a) Harris, W. E. 1996a, AJ, 112, 1487, doi: 10.1086/118116
  * Harris (1996b) —. 1996b, AJ, 112, 1487, doi: 10.1086/118116
  * Harris et al. (2017) Harris, W. E., Blakeslee, J. P., & Harris, G. L. H. 2017, ApJ, 836, 67, doi: 10.3847/1538-4357/836/1/67
  * Harris & van den Bergh (1981) Harris, W. E., & van den Bergh, S. 1981, AJ, 86, 1627, doi: 10.1086/113047
  * Huxor et al. (2011) Huxor, A. P., Ferguson, A. M. N., Tanvir, N. R., et al. 2011, MNRAS, 414, 770, doi: 10.1111/j.1365-2966.2011.18450.x
  * Huxor et al. (2014) Huxor, A. P., Mackey, A. D., Ferguson, A. M. N., et al. 2014, MNRAS, 442, 2165, doi: 10.1093/mnras/stu771
  * Jordán et al. (2007) Jordán, A., McLaughlin, D. E., Côté, P., et al. 2007, ApJS, 171, 101, doi: 10.1086/516840
  * King (1962) King, I. 1962, AJ, 67, 471, doi: 10.1086/108756
  * Kissler-Patig et al. (1997) Kissler-Patig, M., Richtler, T., Storm, J., & della Valle, M. 1997, A&A, 327, 503.  https://arxiv.org/abs/astro-ph/9706303
  * Kundu & Whitmore (2001) Kundu, A., & Whitmore, B. C. 2001, AJ, 121, 2950, doi: 10.1086/321073
  * Larsen (1999) Larsen, S. S. 1999, A&AS, 139, 393, doi: 10.1051/aas:1999509
  * Lomelí-Núñez et al. (2022) Lomelí-Núñez, L., Mayya, Y. D., Rodríguez-Merino, L. H., Ovando, P. A., & Rosa-González, D. 2022, MNRAS, 509, 180, doi: 10.1093/mnras/stab2890
  * McLaughlin (2000) McLaughlin, D. E. 2000, ApJ, 539, 618, doi: 10.1086/309247
  * McQuinn et al. (2017) McQuinn, K. B. W., Skillman, E. D., Dolphin, A. E., Berg, D., & Kennicutt, R. 2017, AJ, 154, 51, doi: 10.3847/1538-3881/aa7aad
  * Nantais & Huchra (2010) Nantais, J. B., & Huchra, J. P. 2010, AJ, 139, 2620, doi: 10.1088/0004-6256/139/6/2620
  * Nantais et al. (2011) Nantais, J. B., Huchra, J. P., Zezas, A., Gazeas, K., & Strader, J. 2011, AJ, 142, 183, doi: 10.1088/0004-6256/142/6/183
  * Olsen et al. (2004) Olsen, K. A. G., Miller, B. W., Suntzeff, N. B., Schommer, R. A., & Bright, J. 2004, AJ, 127, 2674, doi: 10.1086/383297
  * Peng et al. (2006) Peng, E. W., Jordán, A., Côté, P., et al. 2006, ApJ, 639, 95, doi: 10.1086/498210
  * Peng et al. (2008) —. 2008, ApJ, 681, 197, doi: 10.1086/587951
  * Perrett et al. (2002) Perrett, K. M., Bridges, T. J., Hanes, D. A., et al. 2002, AJ, 123, 2490, doi: 10.1086/340186
  * Portegies Zwart et al. (2010) Portegies Zwart, S. F., McMillan, S. L. W., & Gieles, M. 2010, ARA&A, 48, 431, doi: 10.1146/annurev-astro-081309-130834
  * Rhode & Zepf (2003) Rhode, K. L., & Zepf, S. E. 2003, AJ, 126, 2307, doi: 10.1086/378711
  * Santiago-Cortés et al. (2010) Santiago-Cortés, M., Mayya, Y. D., & Rosa-González, D. 2010, MNRAS, 405, 1293, doi: 10.1111/j.1365-2966.2010.16531.x
  * Schlafly & Finkbeiner (2011) Schlafly, E. F., & Finkbeiner, D. P. 2011, ApJ, 737, 103, doi: 10.1088/0004-637X/737/2/103
  * Schommer et al. (1992) Schommer, R. A., Suntzeff, N. B., Olszewski, E. W., & Harris, H. C. 1992, AJ, 103, 447, doi: 10.1086/116074
  * Simanton et al. (2015) Simanton, L. A., Chandar, R., & Whitmore, B. C. 2015, ApJ, 805, 160, doi: 10.1088/0004-637X/805/2/160
  * Sirianni et al. (2005) Sirianni, M., Jee, M. J., Benítez, N., et al. 2005, PASP, 117, 1049, doi: 10.1086/444553
  * Thilker et al. (2007) Thilker, D. A., Bianchi, L., Meurer, G., et al. 2007, ApJS, 173, 538, doi: 10.1086/523853
  * Whitmore et al. (2021) Whitmore, B. C., Lee, J. C., Chandar, R., et al. 2021, MNRAS, 506, 5294, doi: 10.1093/mnras/stab2087
  * Zhang & Fall (1999) Zhang, Q., & Fall, S. M. 1999, ApJ, 527, L81, doi: 10.1086/312412
  * Zinn (1985) Zinn, R. 1985, ApJ, 293, 424, doi: 10.1086/163249

## Tables

Table 1: $HST$ Pointings of NGC 1291

Pointing | HLA Dataset | RA | DEC | Start Time | Filters
---|---|---|---|---|---
|  | J2000 | J2000 |  |
PT1 | hst_10402_01_acs_wfc_ | 3:17:16.46 | -41:04:27.8 | 2005-04-07 15:17:33 | F435W,F555W,F814W
PT2 | hst_10402_04_acs_wfc_ | 3:17:28.12 | -41:07:04.1 | 2005-08-04 17:10:13 | F435W,F555W,F814W
PT3 | hst_13364_87_acs_wfc_ | 3:17:37.87 | -41:08:34.0 | 2013-09-05 19:33:44 | F435W,F606W,F814W
PT4 | hst_10402_02_acs_wfc_ | 3:17:19.50 | -41:10:17.7 | 2005-06-12 13:01:13 | F435W,F555W,F814W

Note. — Start times are different for each filter. The given Hubble data set
name is not complete. The actual listing on HLA is column 2 in the table plus
the corresponding filter.

Table 2: Cluster Catalog

ID | RA | DEC | F435W | F435W Error | F555W | F555W Error | F814W | F814W Error | CI | $R_{eff}$
---|---|---|---|---|---|---|---|---|---|---
| (J2000) | (J2000) | (mag) | (mag) | (mag) | (mag) | (mag) | (mag) | (mag) | (pc)
0 | 49.312229 | -41.103973 | 21.463 | 0.009 | 20.749 | 0.008 | 19.724 | 0.008 | 2.677 | 1.533
1 | 49.306077 | -41.097649 | 24.145 | 0.101 | 23.828 | 0.131 | 23.220 | 0.182 | 2.971 | 3.065
2 | 49.335306 | -41.097234 | 22.873 | 0.033 | 22.008 | 0.026 | 20.905 | 0.022 | 3.041 | 3.768
3 | 49.303996 | -41.096582 | 23.949 | 0.084 | 22.989 | 0.060 | 21.684 | 0.044 | 3.226 | 6.833
4 | 49.318185 | -41.093568 | 23.657 | 0.067 | 22.869 | 0.056 | 21.742 | 0.048 | 2.735 | 1.245

Note. — Sample data of the first 5 clusters in the catalog.


[2]Philipp Sauerteig

# Surrogate Models in Bidirectional Optimization of Coupled Microgrids

Manuel Baumann  Max Planck Institute for Dynamics of Complex Technical
Systems<EMAIL_ADDRESS>Sara Grundel  Max Planck
Institute for Dynamics of Complex Technical Systems, {baumann,grundel}@mpi-
magdeburg.mpg.de Karl Worthmann  Technische Universität Ilmenau,
<EMAIL_ADDRESS>
###### Abstract

The energy transition entails a rapid uptake of renewable energy sources.
Besides physical changes within the grid infrastructure, energy storage
devices and their smart operation are key measures to master the resulting
challenges like, e.g., a highly fluctuating power generation. For the latter,
optimization based control has demonstrated its potential on a microgrid
level. However, if a network of coupled microgrids is considered, iterative
optimization schemes including several communication rounds are typically
used. Here, we propose to replace the optimization on the microgrid level by
using surrogate models either derived from radial basis functions or neural
networks to avoid this iterative procedure. We prove well-posedness of our
approach and demonstrate its efficiency by numerical simulations based on real
data provided by an Australian grid operator.

Zusammenfassung: Die Energiewende bringt einen raschen Zuwachs eneuerbarer
Energiequellen mit sich. Neben den physikalischen Veränderungen der
Netzinfrastruktur spielen Energiespeichereinheiten und deren intelligente
Nutzung eine entscheidende Rolle, um die sich ergebenden Probleme wie z.B. die
stark schwankende Energieerzeugung zu bewältigen. In Bezug auf Letztere haben
optimierungsbasierte Steuerungstechniken ihr Potential auf Microgrid-Ebene
unter Beweis gestellt. Betrachtet man jedoch ein Netzwerk gekoppelter
Microgrids, werden üblicherweise iterative Optimierungsansätze gewählt, welche
mit mehreren Kommunikationsrunden einhergehen. Um derartigen
Kommunikationsschleifen vorzubeugen, schlagen wir vor, den Optimierungsschritt
auf Microgrid-Ebene durch den Einsatz geeigneter Ersatzmodelle zu vermeiden.
Den hier verwendeten Ersatzmodellen liegen zum einen radiale Basisfunktionen
und zum anderen neuronale Netze zugrunde. Wir zeigen, dass unser Ansatz
wohlgestellt ist und demonstrieren die Effizienz anhand numerischer
Simulationen basierend auf realen Daten eines australischen
Verteilnetzbetreibers.

Schlagwörter: Smart Grids, Modellprädiktive Regelung, Verteilte Optimierung,
Ersatzmodelle, Bidirektionale Optimierung, Neuronale Netze, Radiale
Basisfunktionen

##  1 Introduction

The share of renewable energy sources rapidly increases; also due to more and
more installed devices like e.g., solar panels at household-level. Hence,
households become prosumers, i.e., power is not only consumed but also
produced. Therefore, energy generation and distribution takes place in a
distributed way. In particular, energy can be transmitted bidirectionally
between the grid and the prosumers, which results in a paradigm shift in the
grid organization. In addition, prosumers may also possess some kind of energy
storage device in order to manipulate their power demand profiles by either
charging or discharging. From the grid operator's perspective it might be
beneficial that charging decisions are not made based on local information
only. Instead taking into account information on the entire grid may improve
the system-wide operation, e.g., to flatten the overall power demand within
the grid in order to facilitate the power supply. [1]. In order to achieve
this goal, communication is needed. In the future, each household shall be
equipped with a smart meter which yields so-called _smart homes_. Smart meters
collect data and communicate with the grid operator automatically.

A straight-forward way to optimally operate the overall system is to formulate
one large-scale optimization problem and to solve it in a centralized way,
see, e.g. [2]. This approach, however, is hard to realize in practice. One of
the disadvantages is that some central node needs the complete information
about the grid, which is, e.g. due to data privacy, not desirable.
Alternatives are decentralized or distributed optimization algorithms. In [3]
the authors propose a decentralized approach to steer energy storage systems
in order to avoid over-capacity of pole transformers while maintaining a high
charging amount of energy storage systems in low-voltage distribution systems.
The other option mentioned above are distributed optimization methods such as
distributed dual ascent [4], Alternating Direction Method of Multipliers
(ADMM) [5] or Augmented Lagrangian based Alternating Direction Inexact Newton
(ALADIN) [6]. These algorithms use a star-shaped communication topology, i.e.
each smart home communicates only with the grid operator and does not share
any information with its neighbours. Nevertheless, in every iteration each
household has to transmit specific (personal) data to the grid operator, see
also [7] and [8] for an application of ADMM and ALADIN to electrical networked
systems, respectively. In order to exploit the potential of these algorithms
they are typically embedded within a Model Predictive Control (MPC) framework.
MPC is a state-of-the-art technique to tackle optimal control problems by
solving finite-dimensional optimization problems successively, see e.g. [9]
for an introduction to MPC and [10, 11] for MPC approaches in electrical
networks.

An alternate option to steer the power demand of local agents besides battery
control is to schedule so-called controllable loads. Controllable loads can be
shifted in time to avoid bottlenecks in the energy supply, see e.g. [12, 13].
There is also a large potential in the context of stochastic optimization of
smart grids. For weather forecasting methods we refer to [14]. How to
integrate electrical vehicles into the electricity network under uncertainties
is described in [15].

Considering the power networks described so far, it is assumed that exchange
of energy within the grid is possible at any time and does not cause any
losses or additional costs, which might (approximately) hold for domestic
nets, e.g. a town. In this paper, we refer to these grids as microgrids (MGs).
In [16, 17], the concept of coupled MGs is used to tackle large-scale problems
incorporating several MGs. In the latter, the authors show that even if each
single MG is optimally operated, there is still room for improvement if energy
can be exchanged among MGs. Therefore, a second optimization problem is solved
on a higher grid level in order to optimally exchange energy resulting in a
bilevel optimization problem [18].

In [19], the authors propose to replace the distributed optimization routine
on the lower grid level by a surrogate model in order to speed-up the
calculation and further reduce communication effort. Here, Radial Basis
Functions (RBFs) [20] are used to approximate the input-output behaviour of
ADMM within the framework of coupled MGs established in [17]. Besides RBFs
there are various methods to learn the behaviour of a complex function.
Artificial Neural Networks (NNs) are one of the most popular representatives
of modern artificial intelligence techniques and are often used in practice
due to their success in various application fields, see e.g. the survey
article [21]. In [22] the authors forecast loads in a power grid using NNs,
whereas in [23] NNs are used in an optimal power flow framework. The main
advantage of using surrogates is that communication effort can be reduced.

In this paper, we extend the idea of coupled microgrids established in [17] by
proposing an iterative _bi-directional_ optimization routine in order to
improve the overall performance. Due to its iterative structure, however, our
method comes along with a strong need for communication between smart homes
and grid operator. As a remedy we present two approaches to reduce the
communication effort by substituting the optimization on microgrid level via
surrogate models. A main difference compared to [19] lies in the different
input-output map that is replaced by the surrogate models, for which we can
show that each input uniquely determines an (optimal) output. Furthermore, we
also take NNs as potential surrogate models into account and study the
performance of the resulting approximations numerically in an MPC framework.
Our simulations show that the proposed method approximately recovers the
performance based on using ADMM but significantly reduces the communication
burden. The effect of applying surrogate models within MPC extends our
previous work [19] where a surrogate model based on RBFs was only applied in a
static optimization problem.

The paper is structured as follows: In Section 2 we formulate a mathematical
model for coupled microgrids that consists of two hierarchy levels, and
introduce optimization problems corresponding to each of them. In the
consecutive section, we propose an iterative scheme that requires the solution
of a distributed optimization problem on the lower level which is solved using
ADMM. In Section 4, we investigate the impact of disturbances w.r.t. the
lower-level solution on the performance measured in terms of the upper-level
objective function. Based on the results, we propose to replace ADMM by
surrogates in order to reduce communication effort and computation time. The
performance of the optimization scheme incorporating surrogates is analysed in
an MPC framework in Section 5.

##  2 A model for coupled microgrids

We consider a system of coupled microgrids (MGs) and call it a smart grid.
Each MG consists of several residential energy systems (agents) coupled
through the grid operator, which can be seen as Central Entity (CE). The
coupling of the microgrids is done through a network, where some MGs are
connected by a transmission line and others are not connected, cf. Figure 1.

Fig. 1: Upper-level model: Schematic representation of four coupled MGs.
Energy exchange is possible only via transmission lines.

###  2.1 Upper-level model: Energy exchange

We assume that we have $\Xi\in\mathbb{N}$ many MGs which are partially coupled
via transmission lines and can be interpreted as nodes of a non-complete
graph, see Figure 1 where $\Xi=4$. Each MG $\kappa$, for $\kappa\in[1:\Xi]$,
consists of $\mathcal{I}_{\kappa}\in\mathbb{N}$ agents modeled in detail in
Subsection 2.2. We assume that each MG $\kappa$ has an average power demand
$\bar{z}_{\kappa}(n)$ at time $n$. Given this, we can compute the total power
demand $\mathcal{I}_{\kappa}\bar{z}_{\kappa}$ of a MG. The control goal is to
exchange power among the MGs in a way such that a desired quantity
$\bar{\zeta}(n)$ is targeted by controlling the residential storage units. We
specify $\bar{\zeta}(n)$ in Subsection 2.2, and assume for the moment that
this quantity is known to the grid and has advantages for the grid operation.
We assume that this desired quantity is independent of $\kappa$, but this is
not necessary for the rest of the discussion.

Let $\delta_{\nu\kappa}$ describe the percentage of power
$\mathcal{I}_{\nu}\bar{z}_{\nu}(n)$ that is transferred from MG $\nu$ to MG
$\kappa$. We enforce $\delta_{\nu\kappa}$ equals zero if there is no
transmission line between the two MGs. Otherwise, the power demand of a MG
$\kappa$ is given by its own total power demand
$\delta_{\kappa\kappa}\mathcal{I}_{\kappa}\bar{z}_{\kappa}$, where
$\delta_{\kappa\kappa}$ is what remains at the MG, and the sum over the power
received from connected MGs,
$\sum_{\nu\neq\kappa}\delta_{\nu\kappa}\mathcal{I}_{\nu}\bar{z}_{\nu}$. For
each time step in our prediction horizon of length $N\in\mathbb{N}_{\geq 2}$,
we want to match this to the desired power demand starting at time $k$ for $N$
timesteps of each MG in a least-squares sense. The objective function is thus
given by $\mathcal{J}:\mathbb{R}^{\Xi N}\times\mathbb{R}^{\Xi\times\Xi\times
N}\to\mathbb{R}$,

$\displaystyle(\bar{z},\delta)\mapsto\sum_{n=k}^{k+N-1}\sum_{\kappa=1}^{\Xi}\left(\bar{\zeta}(n)\mathcal{I}_{\kappa}-\sum_{\nu=1}^{\Xi}\delta_{\nu\kappa}(n)\eta_{\nu\kappa}\mathcal{I}_{\nu}\bar{z}_{\nu}(n)\right)^{2}.$
(1)

Here, the vector $\bar{z}=(\bar{z}(k),\ldots,\bar{z}(k+N-1)^{\top})^{\top}$
with $\bar{z}(\cdot)\in\mathbb{R}^{\Xi}$ stacks the average power demand per
MG and time step while the matrix
$\eta=(\eta_{\nu\kappa})_{\kappa,\nu=1}^{\Xi}\in[0,1]^{\Xi\times\Xi}$
incorporates efficiencies along the transmission lines.

We are interested in minimizing (1) under the following constraints: All
exchange rates $\delta_{\nu\kappa}$ are within the interval $[0,1]$, sum up to
$1$, and only transfer power in one direction, meaning that either
$\delta_{\nu\kappa}$ or $\delta_{\kappa\nu}$ is zero. Moreover, note that at
this grid level, the average power demands per MG are known. Following [17],
the optimization problem of the upper-level is, thus, formulated over the
exchange rates $\delta$,

$\displaystyle\min_{\delta\in\Delta}\quad$
$\displaystyle\mathcal{J}(\bar{z},\delta)$ (2a)
$\displaystyle\mathrm{s.t.}\quad$
$\displaystyle\sum_{\kappa=1}^{\Xi}\delta_{\nu\kappa}(n)=1$ (2b)
$\displaystyle\delta_{\nu\kappa}(n)\cdot\delta_{\kappa\nu}(n)\leq 0,\
\kappa\neq\nu$ (2c)
$\displaystyle\forall\,\nu,\kappa\in[1:\Xi],\,n\in[k:k+N-1],$

where $\delta_{\nu\kappa}(n)$ denotes the power exchange rate from MG $\nu$ to
MG $\kappa$ at time instance $n$, and
$\Delta\\!:=\\!\left\\{\delta\\!\in\\![0,1]^{\Xi\times\Xi\times
N}\\!\mid\\!\delta_{\nu\kappa}(n)\\!=\\!0,\text{ if no line from }\nu\text{ to
}\kappa\right\\}$. Constraints (2b) and (2c) ensure that the whole energy of
each MG is _scheduled_ and that exchanges via transmission lines can only
occur in one direction during one time step. The definition of the set
$\Delta$ encodes the grid topology and, hence, avoids scheduling energy
exchange in the absence of a transmission line, cf. Figure 1. We denote the
feasible set of (2) by

$\displaystyle\mathbb{D}^{\delta}\;=\;\Set{\delta\in\Delta\subseteq[0,1]^{\Xi\times\Xi\times
N}}{\text{\eqref{Constraint:conservation} - \eqref{Constraint:oneway} hold}}.$

The efficiency of a transmission line does not depend on the direction of the
transfer, i.e. the matrix $\eta$ is symmetric. Furthermore, we assume no loss
without transport, i.e. $\eta_{\kappa\kappa}\equiv 1$ for all
$\kappa\in[1:\Xi]$, in the rest of the paper.

###  2.2 Lower-level model: Single microgrid

As we have seen in the previous section, we consider an average power demand
at each MG as well as some desired quantity $\bar{\zeta}(n)$. In order to
understand these quantities better, we explain the modeling of the MG in all
detail. The basis of our considerations forms the model presented in [24] and
its extension [7]. Therefore, let each subsystem be equipped with an energy
generation device (e.g. roof top photo-voltaic panels) and some storage device
(e.g. a battery). Then, the $i$-th system, $i\in[1:\mathcal{I}_{\kappa}]$ in
MG $\kappa$, can be described by the discrete time system dynamics,

$\displaystyle x_{\kappa_{i}}(n+1)\;$
$\displaystyle=\;\alpha_{\kappa_{i}}x_{\kappa_{i}}(n)+T(\beta_{\kappa_{i}}u_{\kappa_{i}}^{+}(n)+u_{\kappa_{i}}^{-}(n))$
(3a) $\displaystyle z_{\kappa_{i}}(n)\;$
$\displaystyle=\;w_{\kappa_{i}}(n)+u_{\kappa_{i}}^{+}(n)+\gamma_{\kappa_{i}}u_{\kappa_{i}}^{-}(n),$
(3b)

where $x_{\kappa_{i}}(n)$ and $z_{\kappa_{i}}(n)$ denote the State of Charge
(SoC) of the battery and the power demand at time instance
$n\in\mathbb{N}_{0}$, respectively. The latter incorporates the net
consumption, $w_{\kappa_{i}}(n)=\ell_{\kappa_{i}}(n)-g_{\kappa_{i}}(n)$, as
the difference of load and power generation, cf. Figure 2.

Fig. 2: Lower-level model: Star-shaped network of residential energy systems
and central entity (CE). The quantity $w=\ell-g$ in (3b) is obtained from the
data set [25].

The system can be controlled by charging $u_{\kappa_{i}}^{+}$ and discharging
$u_{\kappa_{i}}^{-}$ the battery at each time step. The length of a time step
in hours is denoted by $T>0$, e.g. $T=0.5$ corresponds to a half-hour time
window. The constants
$\alpha_{\kappa_{i}},\beta_{\kappa_{i}},\gamma_{\kappa_{i}}\in(0,1]$ represent
efficiencies w.r.t. self discharge and energy conversion. Furthermore, the
initial SoC $x_{\kappa_{i}}(k)=\hat{x}_{\kappa_{i}}$, where
$k\in\mathbb{N}_{0}$ denotes the current time instance, is assumed to be
known. State and input are subject to the inequality constraints,

$\displaystyle 0\quad\leq\quad$ $\displaystyle x_{\kappa_{i}}(n)$
$\displaystyle\quad\leq\quad C_{\kappa_{i}}$ (4a)
$\displaystyle\underaccent{\bar}{u}_{\kappa_{i}}\quad\leq\quad$ $\displaystyle
u_{\kappa_{i}}^{-}(n)$ $\displaystyle\quad\leq\quad 0$ (4b) $\displaystyle
0\quad\leq\quad$ $\displaystyle u_{\kappa_{i}}^{+}(n)$
$\displaystyle\quad\leq\quad\bar{u}_{\kappa_{i}}$ (4c) $\displaystyle
0\quad\leq\quad$
$\displaystyle\frac{u_{\kappa_{i}}^{-}(n)}{\underaccent{\bar}{u}_{\kappa_{i}}}+\frac{u_{\kappa_{i}}^{+}(n)}{\bar{u}_{\kappa_{i}}}$
$\displaystyle\quad\leq\quad 1.$ (4d)

Here, $C_{\kappa_{i}}\geq 0$ denotes the battery capacity. The last constraint
ensures that the bounds on discharging (4b) and charging (4c) hold even if the
battery is both discharged and charged during one time step. Note that the
case $C_{\kappa_{i}}=0$ covers the case, where not all systems have a storage
device. Since the future net consumption is not known in advance, it is
assumed to be reliably predictable on a short time horizon of size $N$,
$N\in\mathbb{N}_{\geq 2}$, time steps.

For a concise notation we introduce the set
$\mathbb{X}_{\kappa_{i}}=[0,C_{\kappa_{i}}]$ of feasible states, the set
$\mathbb{U}_{\kappa_{i}}=\Set{(u_{\kappa_{i}}^{+},u_{\kappa_{i}}^{-})^{\top}\in\mathbb{R}^{2}}{\eqref{Constraint:discharge}-\eqref{Constraint:chargeNdischarge}\text{
hold}}$ of feasible control pairs and the set

$\displaystyle\mathbb{D}_{\kappa_{i}}\\!=\\!\Set{z_{\kappa_{i}}\\!\in\\!\mathbb{R}^{N}}{\begin{matrix}z_{\kappa_{i}}=\left(z_{\kappa_{i}}(k),\ldots,z_{\kappa_{i}}(k\\!+\\!N\\!-\\!1)\right)^{\top}\\\
\exists\,u_{\kappa_{i}}\in\mathbb{U}^{N}\text{such that }\\\
x_{\kappa_{i}}(k)=\hat{x}_{\kappa_{i}},\text{\eqref{SysDyn} and
\eqref{Constraints} hold}\end{matrix}}$

of feasible outputs over the next $N$ time steps,
$i\in[1:\mathcal{I}_{\kappa}]$, $\kappa\in[1:\Xi]$. Referring to the feasible
sets of a MG $\kappa$ we use the Cartesian product, e.g.
$\mathbb{D}^{(\kappa)}=\mathbb{D}_{\kappa_{1}}\times\ldots\times\mathbb{D}_{\kappa_{\mathcal{I}_{\kappa}}}$
and $z^{(\kappa)}\in\mathbb{D}^{(\kappa)}$. Note that the sets
$\mathbb{D}_{\kappa_{i}}$, $i\in[1:\mathcal{I}_{\kappa}]$, and hence
$\mathbb{D}^{(\kappa)}$ are non-empty, compact, and convex.

The output quantity in (3b) is the power demand $z_{\kappa_{i}}$ of an
individual agent in MG $\kappa$. The average power demand
$\bar{z}_{\kappa}=\frac{1}{\mathcal{I}_{\kappa}}\sum_{i=1}^{\mathcal{I}_{\kappa}}z_{\kappa_{i}}$
in each MG can then be computed from the individual power demands, and is used
as an input to (1). We define $\bar{\zeta}$ in (1) as a stable reference
trajectory by averaging over a past time horizon and over all individual
residential units of all microgrids, the so-called overall average net
consumption as proposed in [7],

$\displaystyle\bar{\zeta}(n)\,=\,\frac{1}{\mathcal{I}\cdot\min\\{N,n+1\\}}\;\sum_{j=n-\min\\{n,N-1\\}}^{n}\;\sum_{i=1}^{\mathcal{I}}\;w_{i}(j)$
(5)

where $\mathcal{I}=\sum_{\kappa=1}^{\Xi}\mathcal{I}_{\kappa}$ denotes the
total number of agents within the entire smart grid. Due to averaging, the
trajectory $\bar{\zeta}$ has little fluctuations and yields advantages for the
grid operation.

Let us for the moment ignore the coupling described in Subsection 2.1. Then,
$\delta(n)$ equals the identity for all $n\in[k:k+N-1]$ and $\mathcal{J}$
becomes

$\displaystyle\bar{z}\mapsto\sum_{n=k}^{k+N-1}\sum_{\kappa=1}^{\Xi}\left(\bar{\zeta}(n)\mathcal{I}_{\kappa}-\mathcal{I}_{\kappa}\bar{z}_{\kappa}(n)\right)^{2}.$

Therefore, the overall objective $\mathcal{J}$ can be decoupled yielding the
local optimization problem

$\displaystyle\min_{z^{(\kappa)}\in\mathbb{D}^{(\kappa)}}\quad$ $\displaystyle
g(\bar{z}_{\kappa})$ (6)

per MG with local objective function $g:\mathbb{R}^{N}\to\mathbb{R}_{\geq 0}$,

$\displaystyle
g(\bar{z}_{\kappa})=\left\|\bar{\zeta}-\bar{z}_{\kappa}\right\|_{2}^{2}.$

###  2.3 Fully coupled optimization problem

We are interested in optimizing the function (1). This function, in general,
depends on $\delta$ as well as on $\bar{z}$. As seen in the previous section,
the average power demand $\bar{z}=\bar{z}(u)$ depends on the control $u$,
which we have to find in such a way that $\mathcal{J}$ is optimal. The overall
optimization problem can be written as

$\min_{\delta\in\mathbb{D}^{\delta},z^{(\kappa)}\in\mathbb{D}^{(\kappa)}}\quad\mathcal{J}(\bar{z},\delta).$
(7)

Note that due to constraint (2c) the optimization of $\mathcal{J}$ w.r.t.
$\delta$ is non-convex. Furthermore, the large scale of the optimization
w.r.t. $\bar{z}$ causes the use of a centralized solver to be expensive. In
addition, using a centralized solver assumes the existence of a global entity
gathering the information of the whole grid, in particular the personal data
of each household, which is undesirable in practice. Hence, solving (7)
centralized is impractical. In the subsequent section we present an approach
to tackle (7) by solving the upper and lower-level problem iteratively. Doing
so we avoid a node with full knowledge in the grid and only communicate
specific aggregated information in each iteration among the agents.

##  3 Bidirectional optimization

We propose to tackle the optimization problem (7) in a bidirectional way, i.e.
we first find an optimal $\bar{z}$ for $\delta$ being the identity and then
optimize (7) w.r.t. $\delta$ for fixed $\bar{z}$ in order to find the optimal
exchange strategy. This typically already gives a considerable improvement,
and has been also done e.g. in [17]. To solve (2) for a fixed $\bar{z}$ is
straight forward, and we use a standard Sequential Quadratic Programming (SQP)
solver. We refer to [26] for an introduction to SQP methods. In this paper we
show how to incorporate the computed energy exchange from the upper level into
the lower-level optimization problem in order to improve the overall
performance.

###  3.1 Bidirectional optimization scheme

Assume that each MG $\kappa$, $\kappa\in[1:\Xi]$, within the smart grid has
already solved its inherent optimization problem (6) and based on the
corresponding solutions $\bar{z}_{\kappa}$ an energy exchange policy
$\delta^{\star}$ has been computed according to (2). This exchange yields an
updated power demand

$\displaystyle\bar{z}_{\kappa}^{+}(n)=\frac{1}{\mathcal{I}_{\kappa}}\sum_{\nu=1}^{\Xi}\delta^{\star}_{\nu\kappa}(n)\,\eta_{\nu\kappa}\,\mathcal{I}_{\nu}\,\bar{z}_{\nu}(n),\
n\in[k:k+N-1],$

and hence, the difference
$\Delta\bar{z}_{\kappa}(n)=\bar{z}_{\kappa}(n)-\bar{z}_{\kappa}^{+}(n)$ in
power demand for all MGs. We are interested in updating $\bar{z}$ in such a
way that our cost function (1),

$\displaystyle\mathcal{J}(\bar{z},\delta^{\star})$
$\displaystyle=\sum_{n=k}^{k+N-1}\sum_{\kappa=1}^{\Xi}\\!\left(\\!\bar{\zeta}(n)\,\mathcal{I}_{\kappa}\\!-\\!\sum_{\nu=1}^{\Xi}\delta^{\star}_{\nu\kappa}(n)\,\eta_{\nu\kappa}\,\mathcal{I}_{\nu}\,\bar{z}_{\nu}(n)\\!\right)^{\\!2}$
$\displaystyle=\sum_{\kappa=1}^{\Xi}\mathcal{I}_{\kappa}^{2}\left\|\bar{\zeta}-\bar{z}_{\kappa}^{+}\right\|_{2}^{2}=\sum_{\kappa=1}^{\Xi}\mathcal{I}_{\kappa}^{2}\left\|(\bar{\zeta}+\Delta\bar{z}_{\kappa})-\bar{z}_{\kappa}\right\|_{2}^{2}$

is minimized further. One could think about fixing $\delta$ and finding an
optimal $\bar{z}$. This, however, leads to an optimization problem not
avoiding communication and coupling all microgrids. The trick here is now to
fix not only the $\delta$ but also the $\bar{z}$-components from all the
microgrids but one. This leads to optimizing
$\mathcal{I}_{\kappa}^{2}\left\|(\bar{\zeta}+\Delta\bar{z}_{\kappa})-\bar{z}_{\kappa}\right\|_{2}^{2}$
locally in each MG, where $\Delta\bar{z}_{\kappa}$ is computed by using the
$\bar{z}_{\kappa}$'s and $\delta$'s from the previous optimization step.
Intuitively, the difference $\Delta\bar{z}_{\kappa}$, $\kappa\in[1:\Xi]$, can
be interpreted as an additional load or generation, and, therefore, as a
change of the desired power demand profile for MG $\kappa$. This yields the
modified lower-level optimization problem

$\displaystyle\min_{z^{(\kappa)}\in\mathbb{D}^{(\kappa)}}\quad
g_{\kappa}(\bar{z}_{\kappa})=\left\|\bar{\zeta}_{\kappa}^{+}-\bar{z}_{\kappa}\right\|_{2}^{2},$
(8)

where $\bar{\zeta}^{+}_{\kappa}=\bar{\zeta}+\Delta\bar{z}_{\kappa}$ In this
formulation the updated reference trajectories $\bar{\zeta}_{\kappa}^{+}$,
$\kappa\in[1:\Xi]$, differ among the single MGs and depend on a given $\delta$
and a given previous $\bar{z}_{\kappa}$. The solution of the newly derived
lower-level optimization problem can be solved with ADMM for all microgrids
independently and in a parallel way.

Based on the updated reference value we solve (8) and (2) to improve the
battery usage and the energy exchange and repeat the optimization until some
terminal condition is satisfied, e.g. performance improvement less than a pre-
defined tolerance or maximal number of iterations exceeded. This procedure is
summarized in Algorithm 1. Note that we only update the reference
$\bar{\zeta}$ on the lower level, since the upper-level optimization problem
(1) does not change.

Algorithm 1 Iterative bidirectional optimization scheme

Input: Current time instance $k\in\mathbb{N}_{0}$, current SoC
$x_{\kappa_{i}}(k)\in\mathbb{X}_{\kappa_{i}}$, prediction horizon
$N\in\mathbb{N}_{\geq 2}$, predicted net consumption
$(w_{i_{\kappa}}(k),\ldots,w_{i_{\kappa}}(k+N-1))^{\top}\in\mathbb{R}^{N}$,
${i_{\kappa}}\in[1\\!:\\!\mathcal{I}_{\kappa}]$, $\kappa\in[1\\!:\\!\Xi]$,
reference trajectory
$(\bar{\zeta}(k),\ldots,\bar{\zeta}(k+N-1))^{\top}\in\mathbb{R}^{N}$, maximal
number $j_{\max}\in\mathbb{N}$ of iterations, and tolerance $\varepsilon>0$.

  1. Initialization:
  2. 1.

Set $j=0$ and $\delta^{0}(n)=I_{\Xi}$ for all $n\in[k:k+N-1]$.

  3. 2.

_Lower level (parallel in $\kappa$)._ Compute $\bar{z}_{\kappa}^{j}$ as the
solution of (6) using ADMM.

  4. 3.

_Upper level._ Given $\bar{z}_{\kappa}^{j}$, solve (2) for $\delta^{j}$ using
SQP.

  5. While

$\quad j<j_{\max}$ and
$\mathcal{J}(\bar{z}^{j-1},\delta^{j-1})-\mathcal{J}(\bar{z}^{j},\delta^{j})>\varepsilon$

  6. Do:
  7. 4.

_Lower level (parallel in $\kappa$)._

    1. (a)

Compute $\bar{\zeta}_{\kappa}^{+}$ from $\bar{z}_{\kappa}^{j}$ and
$\delta^{j}$.

    2. (b)

Solve (8) using ADMM and send $\bar{z}_{\kappa}^{j+1}$ to the upper level.

  8. 5.

_Upper level._ Given $\bar{z}_{\kappa}^{j+1}$, solve (2) for $\delta^{j+1}$
using SQP.

  9. 6.

$j\to j+1$

Neither convergence nor the interpretation of a potential limit of Algorithm 1
is clear a priori. Figure 3,

Fig. 3: Evolution of the costs before and after energy exchange computed
according to the bidirectional optimization scheme described in Algorithm 1,
i.e. $\mathcal{J}(\bar{z}^{j},\delta^{j-1})$ and
$\mathcal{J}(\bar{z}^{j},\delta^{j})$, resp. Note that
$\mathcal{J}(\bar{z}^{1},\delta^{0})$ yields the costs without microgrid
coupling.

however, experimentally shows convergence of the proposed scheme and a
continuous improvement of the upper-level performance index. Here, we ran 10
iterations of the optimization scheme and plotted both the objective function
values before and after the energy exchange within each iteration. The values
stagnate after four iterations indicating that additional iterations do not
further improve the overall performance. The next subsection elaborates on how
to solve (8) in a fully distributed way using ADMM.

######  Remark 1.

The results displayed in Figure 3 and our numerical investigations indicate
convergence to the global minimum of (7). However, a formal (and rigorous)
derivation of this conjecture is left for future research.

###  3.2 Distributed optimization via ADMM

In this section we briefly discuss how to solve the lower-level optimization
problem (6) or (8) using an Alternating Direction Method of Multipliers (ADMM)
approach. We consider a single MG and therefore omit the index $\kappa$. Since
the averaged output quantity appears in the objective function (6) or (8), we
need to introduce an auxiliary variable $a$ in order to decouple the lower-
level optimization in the following way,

$\displaystyle\min_{z,a}\quad$ $\displaystyle
g(\bar{a})=\|\bar{a}-\bar{\zeta}\|_{2}^{2}$ (9a)
$\displaystyle\mathrm{s.t.}\quad$
$\displaystyle\frac{1}{\mathcal{I}}\sum_{i=1}^{\mathcal{I}}a_{i}-\bar{a}=0,\quad
z_{i}-a_{i}=0$ (9b) $\displaystyle
z_{i}\in\mathbb{D}_{i}\quad\forall\,i\in[1:\mathcal{I}].$ (9c)

Note that (9c) is a short-hand notation for the battery dynamics (3)-(4), and
yields a fully decoupled constraint in the variable $z$. ADMM is an
optimization scheme to solve (9) based on the augmented Lagrangian
$\mathcal{L}_{\rho}:\mathbb{R}^{\mathcal{I}N}\times\mathbb{R}^{\mathcal{I}N}\times\mathbb{R}^{\mathcal{I}N}\to\mathbb{R}$,
for $\rho>0$,

$\displaystyle\mathcal{L}_{\rho}(z,a,\lambda)=g(\bar{a})+\sum_{i=1}^{\mathcal{I}}\left(\lambda_{i}^{\top}(z_{i}-a_{i})+\frac{\rho}{2}\left\|z_{i}-a_{i}\right\|_{2}^{2}\right),$

in a distributed way, cf. [5]. Following [7], the ADMM algorithm for (9)
yields the three-step iteration $\ell\mapsto\ell+1$,

$\displaystyle
z_{i}^{\ell+1}=\operatorname*{arg\,min}_{z_{i}\in\mathbb{D}_{i}}z_{i}^{\top}\lambda_{i}^{\ell}+\frac{\rho}{2}\left\|z_{i}-a_{i}^{\ell}\right\|_{2}^{2}$
(10a) $\displaystyle
a^{\ell+1}=\operatorname*{arg\,min}_{a\in\mathbb{R}^{\mathcal{I}N}}g(\bar{a})-\sum_{i=1}^{\mathcal{I}}a_{i}^{\top}\lambda_{i}^{\ell}+\frac{\rho}{2}\|z_{i}^{\ell+1}-a_{i}\|_{2}^{2}$
(10b)
$\displaystyle\lambda_{i}^{\ell+1}=\lambda_{i}^{\ell}+\rho(z_{i}^{\ell+1}-a_{i}^{\ell+1})$
(10c)

until some termination condition is satisfied. Note that (10b) is an
unconstrained optimization problem and can be solved explicitly. The problem
(10a) can be solved in parallel by each battery in the MG introduced for our
model in Section 2.2. Note that scheme (10) assumes communication within the
MG, more precisely, each system $i$, $i\in[1:\mathcal{I}]$, sends its optimal
solution $z_{i}$ to the CE and receives both the updated auxiliary $a_{i}$ and
dual variable $\lambda_{i}$ in return. The variant discussed in [7] avoids
unnecessary communication overhead by returning a broadcast variable instead,
which only incorporates information on the aggregated values.

According to Theorem 3.1 in [7] the optimization scheme (10) converges in the
following sense.

######  Theorem 1.

Consider Problem (9) with $g$ being strictly convex, closed and proper and let
the iterates $(z^{\ell},a^{\ell},\lambda^{\ell})$ be computed according to
(10). Then the following following statements hold true:

  1. 1.

$(z^{\ell}-a^{\ell})_{\ell\in\mathbb{N}_{0}}$ converges to zero for
$\ell\to\infty$,

  2. 2.

$(g(\bar{a}^{\ell}))_{\ell\in\mathbb{N}_{0}}$ converges to the optimal value
$g^{\star}$ of (9),

  3. 3.

$(\lambda^{\ell})_{\ell\in\mathbb{N}_{0}}$ converges to the dual optimal
$\lambda^{\star}$ of (9).

According to [7] and the references [5, Section 3] and [4, Appendix C]
therein, problem (6) fulfils the assumptions of Theorem 1.

##  4 Surrogate models for ADMM

This section is dedicated to surrogate models for the optimization routine
(10) within a single MG. For simplicity of notation we omit the index
$\kappa$.

Due to the distributive structure of ADMM, the residential energy systems do
not need to share information with their neighbours but only with the CE, see
also the star-structure in Figure 2. In each iteration $\ell$ of ADMM,
subsystem $i$ has to transmit its solution $z_{i}^{\ell}$ of (10a) to the CE.
The optimization scheme presented in Algorithm 1, however, suggests to run
ADMM more than once in order to improve the performance w.r.t. to the
objective function (1). In order to avoid unnecessary communication, we
propose to use surrogate models to approximate the optimization routine (10).
More precisely, we are interested in a function which approximates
$\varphi:\mathbb{R}^{N}\times\mathbb{R}^{\mathcal{I}}\times\mathbb{R}^{N}\to\mathbb{R}^{N}$:

$\displaystyle\varphi(\bar{w},x(k),\bar{\zeta})\quad=\quad\bar{z},$ (11)

for all feasible
$(\bar{w},x(k),\bar{\zeta})\in\mathbb{R}^{N}\times\mathbb{X}\times\mathbb{R}^{N}$.
Note that we do neither assume knowledge on the local net consumption $w_{i}$,
$i\in[1:\mathcal{I}]$, nor on the future SoC $x(n)$, $n>k$.

Figure 4 (top) shows that if the approximation (11) is sufficiently accurate,
the impact on the performance of the optimization scheme is negligible.

Fig. 4: Effect of mapping error in (11) (top) and of the approximation via
radial basis functions (RBFs) and a neural network (NN) on the open-loop
performance $\mathcal{J}(\bar{z},\delta)$ within 48 consecutive time steps
(bottom). We use $T=0.5$h in (3a).

Here, the costs $\mathcal{J}(\bar{z},\delta)$ after optimization are
visualized for 48 consecutive time steps (equals $24$-hours simulation time).
In the experiment, we disturbed the ADMM solution in Algorithm 1 by uniformly
distributed additive noise, i. e., $\bar{z}+10^{-p}\cdot d$, where the vector
$d\sim\mathcal{U}(-1,1)$, and $p\in\mathbb{N}_{0}$ denotes the intensity of
the disturbance.

Note that (11) might yield approximations to the solution $\bar{z}$ that
either aggravate the performance w.r.t. (6) compared to the solution $\bar{w}$
associated with $u\equiv 0$ or solutions that correspond to an infeasible
control $\hat{u}\notin\mathbb{U}$. As a remedy, we propose to apply ADMM once
after replacing it by a surrogate in the optimization scheme. More precisely,
first we run Algorithm 1 using a surrogate in Step 4(b) until the while loop
terminates and then we additionally repeat Steps 4 and 5 using ADMM.

###  4.1 Well-posedness

The following proposition states that for equality in (11), a proper mapping
is defined. For a concise notation we replace the index $\kappa_{i}$ by $i$
here.

######  Proposition 1.

Consider $\varphi$ given by (11), where $\bar{z}$ describes the optimal
solution of (6) computed via ADMM, i.e. $\bar{z}=\bar{z}(u^{\star})$. We
assume all hyper-parameter to be fixed meaning that
$\\{T,\alpha_{i},\beta_{i},\gamma_{i},C_{i},\bar{u}_{i},\underaccent{\bar}{u}_{i}\\}$
in (3)-(4) are constant over time for all $i\in[1:\mathcal{I}]$. Then
$\varphi$ is a mapping, i.e. for all
$(\bar{w},x(k),\bar{\zeta})\in\mathbb{R}^{N}\times\mathbb{X}\times\mathbb{R}^{N}$,
there exists a uniquely determined $\bar{z}\in\mathbb{R}^{N}$ such that
$\bar{z}$ is the solution to the optimization problem (6).

  * Proof.

First note that ADMM yields the unique solution of (6), see e.g. [5].
Furthermore, there are no constraints on $z_{i}$, $i\in[1:\mathcal{I}]$, and
the future SoC can be interpreted as an affine function of the current SoC and
the future (dis-) charging rate. Hence, expansion of (3a) and averaging of
(3b) yield,

$\displaystyle\min_{u}\quad\left\|\bar{z}(u)-\bar{\zeta}\right\|_{2}^{2},\quad\text{
subject to}$ $\displaystyle
x_{i}(k\\!+\\!1\\!+\\!n)=\alpha_{i}^{n+1}x_{i}(k)\\!+\\!T\\!\sum_{\ell=k}^{k+n}\\!\alpha_{i}^{n+k-\ell}\big{(}\beta_{i}u_{i}^{+}(\ell)\\!+\\!u_{i}^{-}(\ell)\\!\big{)},$
$\displaystyle x_{i}(k)=\hat{x}_{i},\ \text{ and
constraints\leavevmode\nobreak\ \eqref{Constraints}},$
$\displaystyle\bar{z}(k+n)=\bar{w}(k+n)+\bar{u}^{+}(k+n)+\bar{\gamma}\bar{u}^{-}(k+n),$
$\displaystyle n\in[0:N-1],$

where $\bar{\cdot}$ denotes the corresponding average value w.r.t. all
subsystems, in particular
$\bar{z}(n)=\frac{1}{\mathcal{I}}\sum_{i=1}^{\mathcal{I}}z_{i}(n)$. This
representation of (6) illustrates that the (predicted) average values
$\bar{\zeta}=(\bar{\zeta}(k),\ldots,\bar{\zeta}(k+N-1))^{\top}$ and
$\bar{w}=(\bar{w}(k),\ldots,\bar{w}(k+N-1))$ and the current SoC
$\left\\{x_{i}(k)\right\\}_{i=1}^{\mathcal{I}}$, uniquely determine the
optimal solution $\bar{z}(u^{\star})$ obtained by ADMM. ∎

######  Remark 2.

Note that $\bar{\zeta}\mapsto\bar{z}$ as introduced in [19] does not define a
mapping since the solution $\bar{z}$ of (6) not only depends on the reference
value $\bar{\zeta}$ but also on the future net consumption $\bar{w}$ and the
current SoC $x_{i}(k)$.

###  4.2 Radial basis functions approximation

Radial Basis Functions (RBFs) are used to interpolate functions based on a set
of sampling data. We briefly recap some basics on RBFs. For a detailed
introduction to theory and application see e.g. [20], for a similar approach
where RBFs are used to replace ADMM we refer to [19].

Let $M\in\mathbb{N}$ denote the number of samples. Then, the interpolation
function of (11) is given as the sum of basis functions
$\psi_{m}:\mathbb{R}^{N}\times\mathbb{X}\times\mathbb{R}^{N}\to\mathbb{R}$,
$m\in[1:M]$, and a regularization term
$q:\mathbb{R}^{N}\times\mathbb{X}\times\mathbb{R}^{N}\to\mathbb{R}^{N}$. More
precisely,

$\displaystyle\bar{z}\approx\varphi^{\mathrm{RBF}}(\chi)=\sum_{m=1}^{M}\psi_{m}(\chi)\alpha_{m}+q(\chi),$
(12)

where $\chi=(\bar{w},x(k),\bar{\zeta})$ is the joint inputs of Proposition 1,
and $\alpha_{m}\in\mathbb{R}^{N}$, $m\in[1:M]$. The basis functions are so-
called radial basis functions of the form,
$\psi_{m}(\chi)=\psi(\left\|\chi-\chi_{m}\right\|)$, where the kernel $\psi$
yields support close to the sampling data $\chi_{m}$, $m\in[1:M]$. We choose
an affine linear regularization $q(\chi)=\beta_{0}+B\chi$. Note that different
choices are possible. The missing parameters $\alpha_{m}$, $\beta_{0}$ and $B$
are determined by interpolation conditions, cf. [19, 20].

In Figure 5, a possible fit via RBFs is visualized. Here, we interpolated
given data from two-weeks of optimization ($4540$ data points) based on
sampling data picking each $25$-th data point to train (12). Then, we tested
$\varphi^{\mathrm{RBF}}$ on the following day, and plotted the fitting.

Fig. 5: RBF and neural net fitting of the first component $\bar{z}(k)$ of
$\bar{z}=(\bar{z}(k),\ldots,\bar{z}(k+N-1))^{\top}$ within 48 consecutive time
steps indicating the quality of the approximation (11).

Our implementation is based on the Matlab toolbox DACE [27]. Note that the
evaluation time of the RBF approximation grows with the number of data points
used. Already with 180 data points to train (12) with $N=6$ causes the
function evaluation of $\varphi^{\mathrm{RBF}}$ to be expensive, see Table 1.
Using more data points would no longer yield an advantage over using ADMM
w.r.t. computation time.

###  4.3 Neural networks approximation

Neural Networks (NNs) are a state-of-the-art method in artificial intelligence
frameworks. Based on huge amounts of data $M\gg 1$ they are able to learn and
recognize patterns in complex systems. We consider a NN of $l$-layers as an
approximation to the mapping (11), i.e.,

$\displaystyle\bar{z}\approx\varphi^{\mathrm{NN}}(\chi)=\sigma\left(W^{[l]}\
...\ \sigma(W^{[2]}\chi+b^{[2]})\ ...\ +b^{[l]}\right),$ (13)

where $\sigma$ denotes the sigmoid function, and the weights $W^{[l]}$ and
biases $b^{[l]}$ are determined during the training phase. Here, the number of
neurons at layer $l-1$ and at layer $l$ determine the number of rows and
columns of $W^{[l]}$, respectively. Note that a separate neural network is
trained for each MG. For an introduction to deep learning and neural networks
we refer the reader to [28, 29]. To train the NNs, we used Matlab's built-in
toolbox nftool.

The overall goal of the approximation (13) is to be sufficient in the sense of
the MPC performance shown in Figure 6. Our experiments in Figures 4 (bottom)
and Figure 5 show that with one hidden layer of ten neurons only, a satisfying
approximation on a 24-hours time window can be achieved if the training data
is large enough. Note that NNs benefit from big data. In our case study, we
trained the NN only on data corresponding to two weeks.

##  5 Numerical proof-of-concept

Model Predictive Control (MPC) is a method to tackle optimal control problems
on an infinite time horizon by solving a series of finite dimensional
optimization problems instead, see e.g. [9] for an introduction to non-linear
MPC.

###  5.1 Model predictive control (MPC)

Consider the optimal control problem (6). In order to provide an optimal
control sequence over an arbitrary long time horizon we use MPC. To this end,
at current time instance $k\in\mathbb{N}_{0}$ we assume the future net
consumption
$(w_{i}(k),w_{i}(k+1),\ldots,w_{i}(k+N-1))^{\top}\in\mathbb{R}^{N}$ to be
predicted for all subsystems $i\in[1:\mathcal{I}]$. Based on the prediction,
Algorithm 1 is executed (per MG) to determine control sequences $u_{i}$,
$i\in[1:\mathcal{I}]$, and an exchange strategy $\delta$. Then, only the first
instances $u_{i}(k)$ and $\delta(k)$ are implemented and the time instance is
incremented. Algorithm 2 outlines this MPC scheme.

Algorithm 2 MPC for coupled MGs

Input: Current time instance $k\in\mathbb{N}_{0}$, current SoC
$x_{i}(k)\in\mathbb{X}_{\kappa_{i}}$, prediction horizon $N\in\mathbb{N}_{\geq
2}$, and reference trajectory
$(\bar{\zeta}(k),\ldots,\bar{\zeta}(k+N-1))^{\top}\in\mathbb{R}^{N}$.

  1. Repeat:
  2. 1.

Measure current state $x_{i}(k)$ and update the forecast
$(w_{i}(k),w_{i}(k+1),\ldots,w_{i}(k+N-1))^{\top}$, $i\in[1:\mathcal{I}]$.

  3. 2.

Run Algorithm 1 for all MG to get optimal control sequences
$u_{i}^{\star}=(u_{i}^{\star}(k),\ldots,u_{i}^{\star}(k+N-1))^{\top}$ for all
subsystems $i\in[1:\mathcal{I}]$, and an optimal exchange strategy
$\delta^{\star}=(\delta^{\star}(k),\ldots,\delta^{\star}(k+N-1))$.

  4. 3.

Implement $u_{i}^{\star}(k)$, $i\in[1:\mathcal{I}]$, and $\delta^{\star}(k)$
and shift the time instance $k\to k+1$.

Note that Problems (6) or (8) and (2) have to be solved in order to determine
$\bar{z}^{\star}$ and $\delta^{\star}$ in each MPC iteration. Therefore, the
open-loop costs $\mathcal{J}(\bar{z}^{\star},\delta^{\star})$ can be computed
in each iteration as well (cf. Figure 4). However, since Step 3 in Algorithm 2
suggests to only implement the first instance of the controls computed in Step
2, these costs are not attained. Instead the stage costs

$\displaystyle\sum_{\kappa=1}^{\Xi}\left(\bar{\zeta}(k)\mathcal{I}_{\kappa}-\sum_{\nu=1}^{\Xi}\delta_{\nu\kappa}^{\star}(k)\eta_{\nu\kappa}\mathcal{I}_{\nu}\bar{z}_{\nu}^{\star}(k)\right)^{2}$
(14)

are realized at each time step $k\in\mathbb{N}_{0}$ (cf. Figure 6).

###  5.2 Usage of surrogate models in MPC

We compare the performances using ADMM, RBFs, and NNs on the lower-level, i.e.
in Step 4(b) of Algorithnm 1. In all numerical simulations we set $T=0.5$,
$N=6$, $\Xi=4$, $\mathcal{I}_{1}=50$, and
$\mathcal{I}_{2}=\mathcal{I}_{3}=\mathcal{I}_{4}=10$.111Note that this setting
yields the global optimization problem (7) with more than 1000 variables. The
battery parameters were randomly chosen with mean values $C=0.98$,
$\underaccent{\bar}{u}=-0.24$, and $\bar{u}=0.25$. Based on the battery
capacities we set $\hat{x}_{i}=0.5C_{i}$. In order to incorporate losses along
the transmission lines, we used the efficiency matrix,

$\displaystyle\eta=\begin{bmatrix}1.0&0.9&0.9&0.85\\\ 0.9&1.0&0.0&0.85\\\
0.9&0.0&1.0&0.0\\\ 0.85&0.85&0.0&1.0\end{bmatrix}\text{ in
}\eqref{OP_upper_cost_function}.$

For simplicity of the numerical computation, we only replaced the lower-level
optimization routine for MG 1 and, thus, avoid training a separate surrogate
model for each MG. We used Matlab for implementation.

Results on the MPC closed loop can be found in Figure 6 and Table 1.

Fig. 6: Impact of mapping error (top) and approximation via RBF and NN (bottom) on the stage costs (14) within 48 consecutive time steps. | closed-loop cost | runtime [ms]
---|---|---
no control | 12,228 | —
ADMM | 4,416 | 2.5
RBFs | 4,529 | 1.2
NNs | 5,598 | 0.05
Tab. 1: Comparison of the summed MPC closed-loop performance
$\sum_{k=0}^{47}\sum_{\kappa=1}^{\Xi}\left(\bar{\zeta}(k)\mathcal{I}_{\kappa}-\sum_{\nu=1}^{\Xi}\delta_{\nu\kappa}^{\star}(k)\eta_{\nu\kappa}\mathcal{I}_{\nu}\bar{z}_{\nu}^{\star}(k)\right)^{2}$
and runtime (per call): ADMM vs. RBFs vs. NNs.

In Figure 6 the closed-loop performances of ADMM (black line) compared to
perturbed ADMM, and ADMM (black line) compared to the two surrogate models are
visualized. Similar to the open-loop case, small disturbances in ADMM have
little impact and RBFs outperform the NN. The first column of Table 1 compares
the sum of all MPC closed-loop performances using ADMM, RBFs and a NN while in
column 2 the average runtimes of these approaches are reported. Note that when
using a surrogate, we call ADMM once per MPC iteration. As elaborated in [7]
in each ADMM iteration an $N$-dimensional vector has to be transmitted twice.
Hence, both surrogates reduce the need for communication. Two great advantages
of ADMM are that the local optimization (10a) can be parallelized and the
global optimization is independent of the size of the MG. However, a single
function evaluation such as (12) or (13) is faster than running the entire
ADMM optimization routine.

Note that in column 2 of Table 1 we ignored the communication between smart
homes and CE which is needed to apply ADMM in practice. However, the runtime
of ADMM impairs when executed in an actual smart grid while surrogates do not
require additional communication.

In order to improve the performance of the NN, more sampling data has to be
generated to increase the training set significantly. To avoid large offline
computation times, we chose $N=6$, i.e. a prediction horizon of three hours,
which is rather short compared to [7, 19].

######  Remark 3.

We point out two implementation details to solve (2) efficiently. First, note
that the optimization (2) can be parallelized in $n$, since there is no
coupling. Furthermore, we replace (2c) by

$\displaystyle\delta_{\kappa\nu}(n)\cdot\delta_{\nu\kappa}(n)\;\leq\;\varepsilon$

for some tolerance $\varepsilon>0$ to smooth the feasible set.

##  6 Conclusions

In this paper we recalled an optimization problem arising in large-scale
electrical networks. We proposed an iterative _bidirectional_ optimization
scheme to tackle this problem in a distributed way, and showed numerically
that a small error on the lower level does not have noticeable impact on the
performance. Based on this observation, we replaced the lower-level
optimization by surrogate models using radial basis functions and artificial
neural networks. The numerical results show the potential of using these
surrogates to reduce communication effort and computational time in MPC while
preserving the overall performance.

Funding: The authors gratefully acknowledge funding by the Federal Ministry
for Education and Research (BMBF; grants 05M18EVA and 05M18SIA). Karl
Worthmann is indebted to the German Research Foundation (DFG; grant WO
2056/6-1).

## References

  * [1] S. Parhizi, H. Lotfi, A. Khodaei, and S. Bahramirad.  State of the Art in Research on Microgrids: A Review.  IEEE Access, 3(1):890–925, 2015.
  * [2] D. E. Olivares, C. A. Ca$\tilde{\text{n}}$izares, and M. Kazerani.  A centralized optimal energy management system for microgrids.  In 2011 IEEE Power and Energy Society General Meeting, pages 1–6, 2011.
  * [3] R. Okubo, S. Yoshizawa, Y. Hayashi, S. Kawano, T. Takano, and N. Itaya.  Decentralized Charging Control of Battery Energy Storage Systems for Distribution System Asset Management.  In 2019 IEEE Milan PowerTech, pages 1–6, 2019.
  * [4] D. P. Bertsekas and J. N. Tsitsiklis.  Parallel and Distributed Computation: Numerical Methods.  Belmont, MA, USA: Athena Scientific, 1989.
  * [5] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein.  Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers.  Foundations and Trends in Machine Learning, 3(1):1–122, 2011.
  * [6] B. Houska, J. Frasch, and M. Diehl.  An Augmented Lagrangian Based Algorithm For Distributed Nonconvex Optimization.  SIAM J. Optim., 26(2):1101–1127, 2016.
  * [7] P. Braun, T. Faulwasser, L. Grüne, C. M. Kellett, S. R. Weller, and K. Worthmann.  Hierarchical distributed ADMM for predictive control with applications in power networks.  IFAC J. Syst. Control, 3:10 – 22, 2018.
  * [8] A. Engelmann, Y. Jiang, T. Mühlpfordt, B. Houska, and T. Faulwasser.  Towards Distributed OPF using ALADIN.  IEEE Trans. Power Syst., 34(1):584–594, 2019.
  * [9] L. Grüne and J. Pannek.  Nonlinear Model Predictive Control. Theory and Algorithms.  Springer, London, 2 edition, 2017.
  * [10] M. Khalid and A. V. Savkin.  A model predictive control approach to the problem of wind power smoothing with controlled battery storage.  Renewable Energy, 35(7):1520–1526, 2010.
  * [11] A. Parisio, E. Rikos, and L. Glielmo.  A model predictive control approach to microgrid operation optimization.  IEEE Trans. Control Syst. Technol., 22(5):1813–1827, 2014.
  * [12] G. Graditi, M. L. Di Silvestre, R. Gallea, and E. R. Sanseverino.  Heuristic-Based Shiftable Loads Optimal Management in Smart Micro-Grids.  IEEE Trans. Ind. Informat., 11(1):271–280, 2015.
  * [13] P. Braun, L. Grüne, C. M. Kellett, S. R. Weller, and K. Worthmann.  Model Predictive Control of Residential Energy Systems Using Energy Storage & Controllable Loads.  Progress in Industrial Mathematics at ECMI 2014. Mathematics in Industry, 22:617–623, 2016.
  * [14] R. R. Appino, J. Á. G. Ordiano, R. Mikut, T. Faulwasser, and V. Hagenmeyer.  On the use of probabilistic forecasts in scheduling of renewable energy sources coupled to storages.  Applied Energy, 210:1207–1218, 2018.
  * [15] R. R. Appino, M. Mu$\tilde{\text{n}}$oz-Ortiz, J. Á. G. Ordiano, R. Mikut, V. Hagenmeyer, and T. Faulwasser.  Reliable Dispatch of Renewable Generation via Charging of Time-varying PEV Populations.  IEEE Trans. Power Syst., 34(2):1558–1568, 2018.
  * [16] R. H. Lasseter.  Smart distribution: Coupled microgrids.  Proceedings of the IEEE, 99(6):1074–1082, 2011.
  * [17] P. Braun, P. Sauerteig, and K. Worthmann.  Distributed optimization based control on the example of microgrids.  In M. J. Blondin, P. M. Pardalos, and J. S. Sáez, editors, Computational Intelligence and Optimization Methods for Control Engineering, volume 150 of Springer Optimization and Its Applications, pages 173–200. Springer International Publishing, 2019.
  * [18] A. Sinha, P. Malo, and K. Deb.  A Review on Bilevel Optimization: From Classical to Evolutionary Approaches and Applications.  arXiv preprint arXiv: 1705.06270, 2017.
  * [19] S. Grundel, P. Sauerteig, and K. Worthmann.  Surrogate Models For Coupled Microgrids.  In I. Faragó, F. Izsák, and P. Simon, editors, Progress in Industrial Mathematics at ECMI 2018, volume 30. Springer International Publishing, 2019.
  * [20] M. D. Buhmann.  Radial basis functions: Theory and implementations, volume 12.  Cambridge university press, 2003.
  * [21] O. I. Abiodun, A. Jantan, A. E. Omolara, K. V. Dada, N. A. Mohamed, and H. Arshad.  State-of-the-art in artificial neural network applications: A survey.  Heliyon, 4(11), 2018.
  * [22] H. Zhang, F. Xu, and L. Zhou.  Artificial neural network for load forecasting in smart grid.  In 2010 International Conference on Machine Learning and Cybernetics, volume 6, pages 3200–3205, 2010.
  * [23] P. Siano, C. Cecati, H. Yu, and J. Kolbusz.  Real Time Operation of Smart Grids via FCN Networks and Optimal Power Flow.  IEEE Trans. Ind. Informat., 8(4):944–952, 2012.
  * [24] K. Worthmann, C. M. Kellett, P. Braun, L. Grüne, and S. R. Weller.  Distributed and Decentralized Control of Residential Energy Systems Incorporating Battery Storage.  IEEE Trans. Smart Grid, 6(4):1914–1923, 2015.
  * [25] E. L. Ratnam, S. R. Weller, C. M. Kellett, and A. T. Murray.  Residential load and rooftop PV generation: an Australian distribution network dataset.  Internat. J. Sustain. Energy, 2015.
  * [26] J. Nocedal and S. J. Wright.  Numerical Optimization.  Springer, 2006.
  * [27] S. N. Lophaven, H. B. Nielsen, and J. Søndergaard.  DACE-A Matlab Kriging toolbox, version 2.0.  Technical report, 2002.
  * [28] I. Goodfellow, Y. Bengio, and A. Courville.  Deep Learning.  MIT Press, 2016.
  * [29] C. F. Higham and D. J. Higham.  Deep Learning: An Introduction for Applied Mathematicians.  arXiv preprint arXiv:1801.05894, 2018.