# Matching small $\beta$ functions using centroid jitter and two beam position
monitors

C. A. Lindstrøm<EMAIL_ADDRESS>R. D’Arcy  M. J. Garland  Deutsches
Elektronen-Synchrotron DESY, Notkestraße 85, 22607 Hamburg, Germany  P.
Gonzalez  Deutsches Elektronen-Synchrotron DESY, Notkestraße 85, 22607
Hamburg, Germany  Universität Hamburg, Luruper Chaussee 149, 22761 Hamburg,
Germany  B. Schmidt  Deutsches Elektronen-Synchrotron DESY, Notkestraße 85,
22607 Hamburg, Germany  S. Schröder  Deutsches Elektronen-Synchrotron DESY,
Notkestraße 85, 22607 Hamburg, Germany  Universität Hamburg, Luruper Chaussee
149, 22761 Hamburg, Germany  S. Wesch  J. Osterhoff  Deutsches Elektronen-
Synchrotron DESY, Notkestraße 85, 22607 Hamburg, Germany

###### Abstract

Matching to small beta functions is required to preserve emittance in plasma
accelerators. The plasma wake provides strong focusing fields, which typically
require beta functions on the mm-scale, comparable to those found in the final
focusing of a linear collider. Such beams can be time-consuming to
experimentally produce and diagnose. We present a simple, fast, and
noninvasive method to measure Twiss parameters in a linac using two beam
position monitors only, relying on the similarity of the beam phase space and
the jitter phase space. By benchmarking against conventional quadrupole scans,
the viability of this technique was experimentally demonstrated at the
FLASHForward plasma-accelerator facility.

##  I Introduction

Plasma-wakefield accelerators TajimaDawsonPRL1979 ; ChenPRL1985 ; RuthPA1985
can provide accelerating gradients in the GV/m-range BlumenfeldNature2007 ;
LitosNature2014 , promising smaller and cheaper accelerators JoshiPT2003 ;
LeemansPT2009 . Reaching high energies, needed for X-ray free-electron lasers
MadeyJAP1971 ; CouprieJPB2014  and linear colliders RosenzweigNIMA1998 ;
SchroederPhysRevSTAB2010 ; AdliSnowmass2013 ; LindstromThesis2019  in
particular, will require multiple accelerator stages SteinkeNature2016 ;
LindstromNIMA2016  and hence some form of external beam injection into the
plasma wake.

Since the focusing field from an exposed ion column in a plasma accelerator is
typically very strong, beams must be tightly focused for the beam size not to
oscillate, as this would lead to significant and unacceptable emittance growth
MehrlingPRAB2012 . In terms of Twiss or Courant-Snyder parameters
CourantSnyderAoP1957 , the beta function needs to be matched to

$\beta_{m}=\sqrt{\frac{2E\epsilon_{0}}{ne^{2}}}$ (1)

where $E$ is the beam energy, $n$ is the plasma density, $\epsilon_{0}$ is the
vacuum permittivity, and $e$ is the electron charge. Injecting a GeV-level
beam into a typical plasma accelerator requires beta functions on the mm-
scale. While plasma density ramps MarshPAC2005 ; DornmairPRAB2015 ;
ZhaoPRAB2020  can relax the matching condition by increasing $\beta_{m}$ at
the entry and exit of the accelerator stage, it will nevertheless be
challenging and time-consuming to experimentally produce and diagnose the
required tightly-focused beams.

Conventional beam-focus diagnostics include wire scanners and high-resolution
screens around the focal point, or downstream quadrupoles that point-to-point
image the beam onto a screen—all of which require nontrivial experimental
setups and careful data analysis. This can be inconvenient when matching beams
into a plasma accelerator—a slow multi-parameter optimization process where
fast feedback will be crucial.

In this paper, we present an alternative method for simple, fast, and
noninvasive measurement of small beta functions by using two beam position
monitors (BPMs) to measure the centroid jitter. The technique is based on the
observation that the phase space of the jitter often has similar Twiss
parameters to that of the beam, and can therefore be used as a proxy. While
the method is approximate in nature, it allows online monitoring and iterative
adjustment of the waist location and beta function. This technique was
successfully implemented and experimentally demonstrated at the FLASHForward
AschikhinNIMA2016 ; DArcyRSTA2019  plasma-accelerator facility at DESY.

##  II Beam and jitter phase spaces

Figure 1: Evolution of a matched beam phase space and various possible
mismatched jitter phase spaces through a strong-focusing lattice. Starting out
moderately mismatched (a), the jitter beta functions evolve (b) and appear to
diverge from that of the matched beam beta function. Nevertheless, in the
focus region (c) the mismatched jitter phase spaces are all focused to a
similar waist beta function and waist location as the beam.

The phase space of a beam consists of its particle distribution in
$x$–$x^{\prime}$ space (in one transverse plane). Similarly, the phase space
of the beam centroid jitter—the jitter phase space—is the distribution of beam
centroid offsets in $x$–$x^{\prime}$ space when integrated over a large number
of shots. Therefore, the jitter has its own Twiss parameters and emittance.

The central assumption underpinning this technique is that the Twiss
parameters of the jitter are similar to those of the beam. A significant
consequence of this connection is that it is possible to simply and
noninvasively measure the phase space of the jitter using BPMs, which then
acts as an approximate measurement of Twiss parameters of the beam. It should
be noted that this technique is generally not suitable for measuring the beam
emittance, but this is also not required for matching (see Eq. 1).

While the similarity of the beam and the jitter is not guaranteed, it is
motivated by both experimental observation and theoretical considerations.
Linear accelerators usually have long FODO-like lattices with beta functions
on the 1–10 m scale. This means that magnets and accelerating cavities—sources
of jitter—are typically distributed across a range of phase advances. As a
consequence, the jitter-phase-space ellipse gradually expands while it rotates
to acquire a similar shape to the beam-phase-space ellipse. Conversely, if
there were only a few dominant jitter sources—such as the gun or a single
vibrating quadrupole—the jitter-phase-space ellipse would be
disproportionately stretched in the $x^{\prime}$-dimension at these phases.
Similar beam and jitter-phase-space ellipses can therefore be expected in any
well-commissioned machine where such dominant jitter sources have been
removed. Even if the beam and jitter phase spaces are moderately mismatched,
both will evolve and be focused similarly in a linear-optics lattice—also in
the case of strong focusing, as demonstrated by the example in Fig. 1. It
should be noted that in the presence of strong focusing, the final-focusing
quadrupoles can contribute disproportionately to the jitter phase space due to
large beta functions—special care therefore needs to be taken to ensure their
stability BalikIPAC2017 ; BettNIMA2018 .

The most interesting quantities in the context of matching are the location
and beta function of the focus waist. How inaccurate should we expect the
jitter-based measurement to be? Consider a lattice that focuses the beam to a
small waist, where the beam size is demagnified by a factor $B$. Starting from
matched Twiss parameters $\beta_{0}$ and $\alpha_{0}$, the resulting waist
beta function would be $\beta_{0}/B^{2}$. The transfer matrix of such a
lattice CourantSnyderAoP1957  can be expressed as

$R=\begin{bmatrix}\frac{\cos\psi+\alpha_{0}\sin\psi}{B}&\frac{\beta_{0}}{B}\sin\psi\\\
\frac{B}{\beta_{0}}\left(\alpha_{0}\cos\psi-\sin\psi\right)&B\cos\psi\end{bmatrix},$
(2)

where the phase advance $\psi$ is a free parameter. Consider then a mismatched
jitter with an initial betatron amplitude matrix

$\Sigma_{0}=\begin{bmatrix}\beta&-\alpha\\\ -\alpha&\gamma\end{bmatrix},$ (3)

where $\gamma=(1+\alpha^{2})/\beta$ is the Twiss gamma function. The overall
mismatch can be quantified by the mismatch parameter SandsSLAC1991

$\mathcal{M}=\frac{1}{2}\left(\tilde{\beta}_{e}+\tilde{\gamma}_{e}+\sqrt{(\tilde{\beta}_{e}+\tilde{\gamma}_{e})^{2}-4}\right),$
(4)

where $\tilde{\beta}_{e}=\beta/\beta_{0}$,
$\tilde{\alpha}_{e}=\alpha-\alpha_{0}\beta/\beta_{0}$ and
$\tilde{\gamma}_{e}=(1+\tilde{\alpha}_{e}^{2})/\tilde{\beta}_{e}$ quantify the
normalized error of each Twiss parameter. The mismatch parameter $\mathcal{M}$
is invariant in a linear-optics lattice, whereas the individual Twiss errors
are not.

We can transport the mismatched jitter to the beam waist location (i.e., the
end of the lattice) using

$\Sigma=R\Sigma_{0}R^{T}.$ (5)

The $\Sigma_{11}$ element corresponds to the jitter beta function at the beam
waist location. However, the beam waist does not generally coincide with the
jitter waist, and therefore $\Sigma_{11}$ does not correspond to the waist
beta function of the jitter. Instead, assuming that the focus region consists
only of a drift, the waist beta function equates to the inverse gamma function
($1/\Sigma_{22}$), which can be expressed as

$\beta_{w}=\frac{\beta_{0}}{B^{2}}\left(\frac{1}{\tilde{\alpha}_{e}\sin
2\psi+\tilde{\beta}_{e}\sin^{2}\psi+\tilde{\gamma}_{e}\cos^{2}\psi}\right).$
(6)

Similarly, the shift of the jitter waist location is given by the ratio of the
alpha and the gamma function ($-\Sigma_{12}/\Sigma_{22}$), which is derived to
be

$\delta s_{w}=\frac{\beta_{0}}{B^{2}}\left(\frac{\tilde{\alpha}_{e}\cos
2\psi+\frac{1}{2}(\tilde{\beta}_{e}-\tilde{\gamma}_{e})\sin
2\psi}{\tilde{\alpha}_{e}\sin
2\psi+\tilde{\beta}_{e}\sin^{2}\psi+\tilde{\gamma}_{e}\cos^{2}\psi}\right).$
(7)

Although lengthy, it is easy to see that if the Twiss errors
$\tilde{\alpha}_{e}$, $\tilde{\beta}_{e}$, and $\tilde{\gamma}_{e}$ are all of
order one (i.e., moderately mismatched), the brackets in both Eqs. 6 and 7
become numerical factors also of order one, regardless of the free parameter
$\psi$. This means that the waist beta function of the mismatched jitter
remains similar to the waist beta function $\beta_{0}/B^{2}$ of the matched
beam. Moreover, it implies that the offset of the waist location is also
approximately $\beta_{0}/B^{2}$—of the order of the waist beta function
itself.

In a plasma accelerator, this mismatch leads to an emittance growth for beams
of finite energy spread, as the phase-space ellipse of each energy slice
rotates at a different rate. Fully decohered, the relative emittance growth
saturates at MehrlingPRAB2012

$\frac{\epsilon_{\textmd{sat}}}{\epsilon_{0}}=\frac{1}{2}\left(\mathcal{M}+\frac{1}{\mathcal{M}}\right),$
(8)

which also agrees with simulations. This implies that for a moderate mismatch
($\mathcal{M}$ of order one) the emittance growth is relatively small—e.g., a
mismatch of $\mathcal{M}=2$ leads to an emittance growth of only 25%. Using
the jitter as a proxy is therefore appropriate for a quick first-pass matching
to the plasma, before a final in-situ optimization using the beam.

##  III Two-BPM measurement method

Having connected the phase space of the beam to that of the jitter, the
problem has been reduced to measuring the jitter phase space. This can be done
quickly and noninvasively with a multi-shot measurement using two BPMs—see
Fig. 2 for a conceptual setup. Correlated offset data is required to measure
the position and angle of each shot, which for a ballistic orbit (i.e., no
magnets between the BPMs) is given by

$x^{\prime}=\frac{x_{2}-x_{1}}{\Delta s},$ (9)

where $x_{1}$ and $x_{2}$ are the upstream and downstream centroid offsets,
respectively, and $\Delta s$ is the separation of the two BPMs. Both
transverse planes can be measured simultaneously. As the number of shots
increases, the jitter phase space will gradually build up, assuming that the
optics remains unchanged. Whenever the optics does change, the measurement
must be restarted.

Given that no scan is performed, the data can be analyzed immediately from the
start of the measurement, then re-analyzed with every additional shot,
gradually increasing the precision. As the number of shots $N$ increases, the
relative measurement error of Twiss parameters and jitter emittance will be
approximately $1/\sqrt{N}$. Since the connection between the beam and the
jitter phase space is only approximate, it will rarely be necessary to require
better than about 10% precision (corresponding to 50–500 shots). In a typical
accelerator with a 1–10 Hz repetition rate, this allows quasi-online
monitoring on a few-tens-of-seconds timescale.

Figure 2: Basic experimental setup, with two BPMs surrounding the focus
region, separated by a ballistic orbit. The measurement can also be
generalized to nonballistic orbits, where there are magnets also between the
BPMs.

The measurement can also be generalized to nonballistic orbits (i.e., with
magnets between the BPMs), as is relevant to for instance plasma accelerators
with strong permanent quadrupoles close to the plasma entrance YaminIPAC2019 .
In this case, the angle at the upstream BPM can be calculated using

$x_{1}^{\prime}=\frac{x_{2}-M_{11}x_{1}}{M_{12}},$ (10)

where $M$ is the transfer matrix between the two BPMs. However, predicting the
evolution of beta functions with mm-level accuracy then requires very accurate
(per-mille-level) measurements of monitor locations, quadrupole locations,
field strengths and beam energy—just like for a quadrupole scan. A ballistic
measurement is comparatively simple, and hence always preferable if possible,
as only an accurate measurement of BPM locations (and relative-offset
calibrations) is required. For the remainder of this paper we will, therefore,
assume that the BPMs are separated by only a drift space in order to
facilitate ballistic measurements.

##  IV Resolution limits

The main limitation of this technique stems from the finite resolution of
BPMs. In measuring the jitter-phase-space ellipse, the width of each angle-
slice (i.e., the position jitter at the waist) must be well resolved, which
limits how small a waist beta function can be measured.

###  IV.1 Analytic model

To calculate this resolution limit, we consider the apparent covariance matrix
of the jitter at the upstream BPM location

$\textmd{cov}(x,x^{\prime})=\begin{bmatrix}\langle
x^{2}\rangle+\sigma^{2}&\langle xx^{\prime}\rangle-\frac{\sigma^{2}}{\Delta
s}\\\ \langle xx^{\prime}\rangle-\frac{\sigma^{2}}{\Delta s}&\langle x^{\prime
2}\rangle+2\frac{\sigma^{2}}{\Delta s^{2}}\end{bmatrix},$ (11)

where $\sigma$ is the BPM resolution. The true covariances of the jitter can
be expressed in terms of its waist parameters as $\langle
x^{2}\rangle=\epsilon(\beta_{w}+s_{w}^{2}/\beta_{w})$, $\langle
xx^{\prime}\rangle=-\epsilon s_{w}/\beta_{w}$, and $\langle x^{\prime
2}\rangle=\epsilon/\beta_{w}$, where $\epsilon$ is the geometric jitter
emittance and $s_{w}$ is the distance from the upstream BPM to the jitter
waist.

The measured jitter emittance for this finite BPM resolution is given by the
determinant of Eq. 11

$\hat{\epsilon}^{2}=\epsilon^{2}+\frac{\sigma^{2}}{\Delta
s^{2}}\frac{\epsilon}{\beta_{w}}\left(s_{w}^{2}+(\Delta
s-s_{w})^{2}+2\beta_{w}^{2}\right)+\frac{\sigma^{4}}{\Delta s^{2}}.$ (12)

Employing the same logic as in Sec. II (for Eqs. 6 and 7), we can find the
measured waist beta function from the inverse of the measured gamma function

$\hat{\beta}_{w}=\frac{\hat{\epsilon}\beta_{w}}{\epsilon+2\frac{\beta_{w}\sigma^{2}}{\Delta
s^{2}}},$ (13)

where $\hat{\epsilon}$ can be substituted from Eq. 12, as well as the measured
waist location from the ratio of the measured alpha and gamma functions

$\hat{s}_{w}=\frac{s_{w}+\frac{\beta_{w}\sigma^{2}}{\epsilon\Delta
s}}{1+2\frac{\beta_{w}\sigma^{2}}{\epsilon\Delta s^{2}}}.$ (14)

Equations 12–14 establish three resolution regimes: (1) well-resolved, (2)
distorted, and (3) fully saturated. To avoid any distortion whatsoever, the
BPM resolution must be better than

$\sigma\ll\sqrt{\frac{\epsilon\beta_{w}\Delta s^{2}}{s_{w}^{2}+(\Delta
s-s_{w})^{2}+2\beta_{w}^{2}}},$ (15)

found by requiring the quadratic $\sigma^{2}$-term in Eq. 12 to be smaller
than the constant $\epsilon^{2}$-term. To avoid saturation (i.e., noise
dominating the signal), the resolution should be better than

$\sigma\ll\sqrt{\frac{\epsilon\Delta s^{2}}{2\beta_{w}}},$ (16)

found by demanding the $\sigma^{2}$-term in the denominator of Eq. 13 be
smaller than the $\epsilon$-term. Encouragingly, the measurement of the waist
location is not affected by the distortion limit, and instead only by the
significantly larger saturation limit. This is because the waist location is
only related to the phase-space correlation and not its area. These regimes
are demonstrated by the example in Fig. 3, which also shows exact agreement
with Monte Carlo simulations of two finite-resolution BPMs.

Figure 3: Simulated two-BPM measurements over a large range of finite BPM
resolutions. The BPMs are spaced 1 m apart and the jitter is focused 0.3 m
from the upstream BPM with a 10 mm beta function and 0.1 mm mrad normalized
emittance. Monte Carlo simulations (average of 107 shots per resolution)
demonstrate that this analytic model is exact.

In a typical case where the waist beta function is small compared to the BPM
separation ($\beta_{w}\ll\Delta s$) and the waist is approximately half way
between the BPMs ($s_{w}\approx\Delta s/2$), Eq. 15 simplifies to
$\sigma\ll\sqrt{2\epsilon\beta_{w}}$—therefore the BPM resolution should be
smaller than the position jitter at the waist. This limit informs the choice
of BPM technology required for the application in question.

###  IV.2 Overcoming the resolution limit

For matching into a plasma accelerator with mm-scale beta functions and sub-µm
jitter emittances, a very high BPM resolution is required. Ideally, this is
achieved using state-of-the-art cavity BPMs, which can provide sub-100 nm
resolution (depending on the charge distribution) KimPRAB2012 ; BettIPAC2018 .
However, if such BPMs are not available, measuring the waist beta function may
require going beyond the resolution limit. This is in principle possible to
do, if the emittance of the jitter is already known.

Just like the emittance of the beam, the jitter emittance is preserved in a
linear-optics lattice (assuming it contains no significant jitter sources).
Therefore one of two alternative measurement techniques can be utilized: (1)
Relax the strength of the focusing until the jitter waist is well
resolved—giving different Twiss parameters, but the same emittance. (2)
Simultaneously perform a similar measurement with two other BPMs just upstream
or just downstream, where the focusing is relaxed compared to the focus
region. The first method requires the jitter emittance to persist in time,
whereas the second requires it to persist in space. Both methods assume
negligible chromaticity or that energy slices are measured separately (see
Sec. V.2).

When the jitter emittance is known, the analysis simplifies greatly. The waist
beta function can be calculated using

$\beta_{w}=\frac{\epsilon\Delta s^{2}}{\langle(x_{2}-x_{1})^{2}\rangle},$ (17)

based on the variance of the angle jitter (Eq. 9), and the waist location is
simply

$s_{w}=\frac{\Delta s}{1-\frac{\partial x_{2}}{\partial x_{1}}},$ (18)

where $\frac{\partial x_{2}}{\partial x_{1}}$ is the slope of the correlation
between the two BPM readings.

##  V Measurements at FLASHForward

Experimental demonstration of the two-BPM method was performed at the
FLASHForward facility at DESY, which uses a 1 GeV electron beam from the FLASH
free-electron-laser facility AckermannNatPhot2007 . FLASH provides high-charge
(up to 1 nC), low-emittance (1 mm mrad) bunches with relatively small centroid
jitter. After an approximately 150 m long linac, the bunches are diverted into
the FLASHForward beamline. Here, a dispersive section allows for advanced
energetic collimation SchroederIOP2020 , then a final-focusing section
[depicted in Fig. 1(b)] tightly focuses the beam into a plasma accelerator.
Downstream of the plasma is a suite of beam diagnostics, in particular a
dipole spectrometer with quadrupoles for point-to-point beam imaging.

###  V.1 Comparison to quadrupole scans

Figure 4: Measured jitter phase spaces at the location of the upstream BPM for
two different optics settings, indicating cm-scale waist beta functions
focused at two waist locations 60 mm apart. A small distortion from a finite
BPM resolution was taken into account when calculating the jitter parameters.
These measurements should be compared to the corresponding quadrupole scans in
Fig. 5. Each dataset consists of 210 shots, giving an estimated relative error
of 7%. Figure 5: Quadrupole scans on a downstream spectrometer screen
performed during the two-BPM measurements in Fig. 4, (a) imaging the beam from
a range of object planes around the beam waist. (b) The variation of
horizontally projected beam size for each object plane indicates that the beam
was focused to a small waist beta function (27–29 mm) close to the center of
the focus region—only moderately mismatched from measured jitter. The screen
resolution was accounted for in the calculation of beam parameters. Error bars
represent the standard error of the mean.

To test the assumptions in Sec. II and the applicability of the method, a
detailed comparison of the measured jitter and beam phase spaces was
performed. A strong-focusing optic was set up to focus bunches with an energy
of 678 MeV and charge 290 pC down to a cm-scale beta function at the location
of the plasma accelerator module (which had been removed from the beam path).
Surrounding this focus region were two cavity BPMs LipkaIBIC2015 ;
NolleIBIC2019  with a resolution of 0.9 µm, separated by 1.073 m, and
approximately equidistant from the nominal focus point.

Two datasets were collected, using slightly different final-focusing optics
with the beam focused at two locations 60 mm apart. Figure 4 shows the
measured jitter phase space for each of these two settings using the two-BPM
method. The presence of outliers (as seen in Fig. 4) can significantly skew
the calculation of phase-space parameters, and thus an outlier-cleaning method
was applied: (i) translate the jitter to the waist location from the BPM
correlation (Eq. 18), (ii) perform Gaussian fits of both the $x$ and
$x^{\prime}$ distributions, (iii) remove all shots beyond $\pm 5\sigma$, and
then (iv) undo the translation from (i). Finally, a small distortive effect
from the finite BPM resolution was removed by numerically solving Eqs. 12–14
for the true jitter-phase-space parameters.

At the same time, an object-plane scan was performed with the downstream
quadrupoles (after the second BPM), imaging the beam onto a LANEX screen with
a resolution of 57 µm. Figure 5 shows the corresponding measurement of the
beam waist. Note that the spectrometer limits the measurement to the
horizontal plane, as the dipole disperses vertically. No chromaticity was
observed on the screen.

The two-BPM measurement agrees with the quadrupole scan measurement to an
acceptable level. The waist beta function of the jitter (18–19 mm) differs
from that of the beam (27–29 mm) by about 35%, and the jitter waist location
is offset from the beam waist location by 15 mm—on the same scale as the waist
beta function, as expected. Based on these numbers, the beam–jitter mismatch
parameter was calculated to be $\mathcal{M}=2.1$–$2.2$, implying that the
phase space of the jitter was indeed closely matched to that of the beam. If
used to match into a plasma accelerator (where the jitter phase space would be
matched), the expected emittance growth of the beam from mismatching (Eq. 8)
would be 28–33%.

Figure 6: Comparison of the centroid jitter measured on the spectrometer
screen and with the two-BPM method (artificially transported through the
quadrupoles to the spectrometer location). Error bars represent the standard
error of the mean. A statistically significant agreement is observed,
verifying the accuracy of both measurement methods.

As an additional cross-check of the jitter measurements, the centroid jitter
was also measured directly on the spectrometer (see Fig. 6). This was used to
verify the accuracy of the distances and quadrupole field strength
calibrations used for the quadrupole scans, as well as to fine-tune the value
of the BPM resolution and calibrations.

###  V.2 Slice-by-slice measurements

Chromaticity, where Twiss parameters change with energy ZyngierLAL1977 ;
MontagueLEP1979 , can be a concern when tightly focusing beams of finite
energy spread LindstromPRAB2016 . This is especially important in energy-
chirp-based two-bunch experiments where a trailing bunch needs to be exactly
matched into the plasma wake behind a different-energy driver bunch.

Figure 7: Chromaticity measurement using an energy-slice scan. (a) Jitter
phase space for a range of energy slices. Significant dispersion (energy-
dependent offsets) can also be observed. Each step consists of 500 shots. (b)
The corresponding evolution of the beta function for each slice in the focus
region, indicating a highly chromatic focus.

Measuring chromaticity with the two-BPM technique requires it to be combined
with an energy filter—each energy slice sufficiently narrow to have an
achromatic focus. At FLASHForward this is accomplished using an energetic
collimator SchroederIOP2020 . Moving both the high- and low-energy collimators
together, thin slices with 0.1% root-mean-square (rms) energy spread could be
made. Figure 7 shows the result of such an energy-slice scan around a mean
energy of 1120 MeV, indicating a highly chromatic focus in the vertical plane.
The waist beta functions are relatively consistent (10–20 mm) across all
slices, whereas the waist location shifted significantly between the highest
and lowest energy slice (by 80 mm). In the horizontal plane (not shown in Fig.
7), the waist location spanned only 10 mm—considerably less chromatic. This
asymmetric chromaticity is expected in a quadrupole-based final-focus system,
where the beam is more strongly defocused in one plane (typically the vertical
plane) before being focused to a waist.

Taking into account all the information gathered in an energy-slice scan, we
can extract a partial 5D beam tomography. As seen in Fig. 7(a), the average
position and angle of each individual energy slice is also measured, and
therefore both the beam centroid and the (emittance-normalized) beam size of
each slice is known in both planes. This in-situ tomography allows not only
slice-specific matching, but also measurement and removal of any bunch
dispersion. For a linearized longitudinal phase space, dispersion corresponds
to a bunch tilt or curvature, which in a plasma wake leads to emittance growth
AssmannNIMA1998 ; LindstromIPAC2016  and potentially a hosing instability
WhittumPRL1991 ; HuangPRL2007 ; MehrlingPRL2017 . Finally, combining such a
two-BPM tomography with longitudinal-phase-space data from a transverse-
deflecting cavity allows a 6D phase space to be reconstructed—important for
realistic simulations and detailed optimization of the external injection
process.

##  VI Conclusions

We have shown that centroid jitter measured by two BPMs can be used to quickly
estimate Twiss parameters in a region of strong focusing. While being an
approximate measurement, it can significantly speed up the complex and
delicate beam setup procedure needed to properly match into a plasma
accelerator, and allows noninvasive online monitoring of the beam focus.
Experiments were successfully performed at FLASHForward to verify this
technique, by comparing the two-BPM measurement to a conventional quadrupole
scan. Already in routine use for plasma-wakefield experiments at FLASHForward,
it is clear that the power of this method lies in its simplicity.

###### Acknowledgements.

The authors would like to thank K. Flöttmann, S. Karstensen, K. Ludwig, F.
Marutzky, A. Rahali, A. Schleiermacher and S. Thiele, as well as the FLASH
management for their scientific and technical support. This work was supported
by the Helmholtz ARD program.

## References

  * (1) T. Tajima and J. M. Dawson, Laser electron accelerator, Phys. Rev. Lett. 43, 267 (1979).
  * (2) P. Chen, J. M. Dawson, R. W. Huff and T. Katsouleas, Acceleration of electrons by the interaction of a bunched electron beam with a plasma, Phys. Rev. Lett. 54, 693 (1985).
  * (3) R. D. Ruth, A. W. Chao, P. L. Morton and P. W. Wilson, A plasma wake field accelerator, Part. Accel. 17, 171 (1985).
  * (4) I. Blumenfeld, C. E. Clayton, F.-J. Decker, M. J. Hogan, C. Huang, R. Ischebeck et al., Energy doubling of 42 GeV electrons in a metre-scale plasma wakefield accelerator, Nature (London) 445, 741 (2007).
  * (5) M. Litos, E. Adli, W. An, C. I. Clarke, C. E. Clayton, S. Corde et al., High-efficiency acceleration of an electron beam in a plasma wakefield accelerator, Nature (London) 515, 92 (2014).
  * (6) C. Joshi and T. Katsouleas, Plasma accelerators at the energy frontier and on tabletops, Phys. Today 56, 6, 47 (2003).
  * (7) W. P. Leemans and E. Esarey, Laser-driven plasma-wave electron accelerators, Phys. Today 62, No. 3, 44 (2009).
  * (8) J. M. J. Madey, Stimulated emission of bremsstrahlung in a periodic magnetic field, J. Appl. Phys. 42, 1906 (1971).
  * (9) M. E. Couprie, A. Loulergue, M. Labat, R. Lehe and V. Malka, Towards a free electron laser based on laser plasma accelerators, J. Phys. B 47, 234001 (2014).
  * (10) J. Rosenzweig, N. Barov, A. Murokh, E. Colby and P. Colestock, Towards a plasma wake-field acceleration-based linear collider, Nucl. Instrum. Methods Phys. Res. A 410, 532 (1998).
  * (11) C. B. Schroeder, E. Esarey, C. G. R. Geddes, C. Benedetti and W. P. Leemans, Physics considerations for laser-plasma linear colliders, Phys. Rev. ST Accel. Beams 13, 101301 (2010).
  * (12) E. Adli, J.-P. Delahaye, S. J. Gessner, M. J. Hogan, T. O. Raubenheimer, W. An, C. Joshi and W. B. Mori, A beam driven plasma-wakefield linear collider: from Higgs factory to multi-TeV, SLAC Report No. SLAC-PUB-15426, 2013.
  * (13) C. A. Lindstrøm, Emittance growth and preservation in a plasma-based linear collider, Ph.D. thesis, University of Oslo, 2019.
  * (14) S. Steinke, J. van tilborg, C. Benedetti, C. G. R. Geddes, C. B. Schroeder, J. Daniels et al., Multistage coupling of independent laser-plasma accelerators, Nature (London) 530, 190 (2016).
  * (15) C. A. Lindstrøm, E. Adli, J. M. Allen, J. P. Delahaye, M. J. Hogan, C. Joshi, P. Muggli, T. O. Raubenheimer and V. Yakimenko, Staging optics considerations for a plasma wakefield acceleration linear collider, Nucl. Instrum. Methods Phys. Res. A 829, 224 (2016).
  * (16) T. Mehrling, J. Grebenyuk, F. S. Tsung, K. Floettmann and J. Osterhoff, Transverse emittance growth in staged laser-wakefield acceleration, Phys. Rev. ST Accel. Beams 15, 111303 (2012).
  * (17) E. D. Courant and H. S. Snyder, Theory of the alternating-gradient synchrotron, Ann. Phys. 3, 1 (1958).
  * (18) K. A. Marsh, C. E. Clayton, D. K. Johnson, C. Huang, C. Joshi, W. Lu et al., Beam matching to a plasma wake field accelerator using a ramped density profile at the plasma boundary, Proceedings of PAC2005, Knoxville, TN, USA (JACoW, Geneva, 2005), p. 2702.
  * (19) I. Dornmair, K. Floettmann and A. R. Maier, Emittance conservation by tailored focusing profiles in a plasma accelerator, Phys. Rev. ST Accel. Beams 18, 041302 (2015).
  * (20) Y. Zhao, W. An, X. Xu, F. Li, L. Hildebrand, M. J. Hogan, V. Yakimenko, C. Joshi and W. B. Mori, Emittance preservation through density ramp matching sections in a plasma wakefield accelerator, Phys. Rev. Accel. Beams 23, 011302 (2020).
  * (21) A. Aschikhin, C. Behrens, S. Bohlen, J. Dale, N. Delbos, L. di Lucchio et al., The FLASHForward facility at DESY, Nucl. Instrum. Methods Phys. Res. A 806, 175 (2016).
  * (22) R. D’Arcy, A. Aschikhin, S. Bohlen, G. Boyle, T. Brümmer, J. Chappell et al., FLASHForward: plasma wakefield accelerator science for high-average-power applications, Philos. Trans. Royal Soc. A 377, 20180392 (2019).
  * (23) G. Balik, B. Aimard, L. Brunetti and B. Caron, Proof of concept of CLIC final focus quadrupoles stabilization, Proceedings of IPAC2017, Copenhagen, Denmark (JACoW, Geneva, 2017), p. 1290.
  * (24) D. R. Bett, C. Charrondière, M. Pateckia, J. Pfingstner, D. Schulte, R. Tomás et al., Compensation of orbit distortion due to quadrupole motion using feed-forward control at KEK ATF, Nucl. Instrum. Methods Phys. Res. A 895, 10 (2018).
  * (25) M. Sands, A beta mismatch parameter, SLAC Report No. SLAC-AP-85, 1991.
  * (26) S. Yamin, R. W. Assmann, U. Dorda, F. Lemery, B. Marchetti, E. Panofski and P. A. Walker, Design considerations for permanent magnetic quadrupole triplet for matching into laser driver wake field acceleration experiment at SINBAD, Proceedings of IPAC2019, Melbourne, Australia (JACoW, Geneva, 2019), p. 143.
  * (27) Y. I. Kim, R. Ainsworth, A. Aryshev, S. T. Boogert, G. Boorman, J. Frisch et al., Cavity beam position monitor system for the Accelerator Test Facility 2, Phys. Rev. ST Accel. Beams 15, 042801 (2012).
  * (28) D. R. Bett, N. Blaskovic Kraljevic, R. M. Bodenstein, T. Bromwich, P. N. Burrows, G. B. Christian et al., Performance of nanometre-level resolution cavity beam position monitors at ATF2, Proceedings of IPAC2018, Vancouver, BC, Canada (JACoW, Geneva, 2018), p. 1212.
  * (29) W. Ackermann, G. Asova, V. Ayvazyan, A. Azima, N. Baboi, J. Bähr et al., Operation of a free-electron laser from the extreme ultraviolet to the water window, Nat. Photonics 1, 336 (2007).
  * (30) S. Schröder, K. Ludwig, A. Aschikhin, R. D’Arcy, M. Dinter, P. Gonzalez et al., Tunable and precise two-bunch generation at FLASHForward, arXiv:2005.12071 [J. Phys.: Conf. Ser. (to be published)].
  * (31) D. Lipka, N. Baboi, D. Noelle, G. Petrosyan, S. Vilcins, R. Baldinger et al., FLASH undulator BPM commissioning and beam characterization results, Proceedings of IBIC2014, Monterey, CA, USA (JACoW, Geneva, 2015), p. 315.
  * (32) D. Nölle et al., The diagnostic system at the European XFEL; commissioning and first user operation, Proceedings of IBIC2018, Shanghai, China (JACoW, Geneva, 2019), p. 162.
  * (33) H. Zyngier, Strategy for correcting for chromaticity, LAL-77/35 (Laboratoire de l’Accélérateur Linéaire, Orsay, 1977).
  * (34) B. W. Montague, Linear optics for improved chromaticity correction, LEP Note 165 (CERN, Geneva, 1979).
  * (35) C. A. Lindstrøm and E. Adli, Design of general apochromatic drift-quadrupole beam lines, Phys. Rev. Accel. Beams 19, 071002 (2016).
  * (36) R. Assmann and K. Yokoya, Transverse beam dynamics in plasma-based linacs, Nucl. Instrum. Methods Phys. Res. A 410, 544 (1998).
  * (37) C. A. Lindstrøm, E. Adli, J. Pfingstner, E. Marín and D. Schulte, Transverse tolerances of a multi-stage plasma wakefield accelerator, Proceedings of IPAC2016, Busan, Korea (JACoW, Geneva, 2016), p. 2561.
  * (38) D. H. Whittum, W. M. Sharp, S. S. Yu, M. Lampe and G. Joyce, Electron-hose instability in the ion-focused regime, Phys. Rev. Lett. 67, 991 (1991).
  * (39) C. Huang, W. Lu, M. Zhou, C. E. Clayton, C. Joshi, W. B. Mori et al., Hosing instability in the blow-out regime for plasma-wakefield acceleration, Phys. Rev. Lett. 99, 255001 (2007).
  * (40) T. J. Mehrling, R. A. Fonseca, A. Martinez de la Ossa and J. Vieira, Mitigation of the hose instability in plasma-wakefield accelerators, Phys. Rev. Lett. 118, 174801 (2017).


# Load Shifting Versus Manual Frequency Reserve:
Which One is More Appealing to Flexible Loads?

Peter A. V. Gade, Trygve Skjøtskift, Charalampos Ziras, Henrik W. Bindner, and
Jalal Kazempour P. A. V. Gade, C. Ziras, H. W. Bindner and J. Kazempour are
with the Technical University of Denmark , 2800 Kgs. Lyngby, Denmark (e-mail:
{pega, chazi, hwbi, jalal}@dtu.dk).P. A. V. Gade and T. Skjøtskift are with
IBM Client Innovation Center, Copenhagen, Denmark (e-mail:
Trygve.Skjotskift@ibm.com)This work was supported in part by Innovation Fund
Denmark, Grant Agreement No. 0153-00205B.

###### Abstract

This paper investigates how a thermostatically controlled load can deliver
flexibility either in form of manual frequency restoration reserves (mFRR) or
load shifting, and which one is financially more appealing to such a load. A
supermarket freezer is considered as a representative flexible load, and a
grey-box model describing its temperature dynamics is developed using real
data from a supermarket in Denmark. Taking into account price and activation
uncertainties, a two-stage stochastic mixed-integer linear program is
formulated to maximize the flexibility value from the freezer. For practical
reasons, we propose a linear policy to determine regulating power bids, and
then linearize the mFRR activation conditions through the McCormick relaxation
approach. For computational ease, we develop a decomposition technique,
splitting the problem to a set of smaller sub-problems, one per scenario.
Examined on an out-of-sample simulation based on real Danish spot and
balancing market prices in 2022, load shifting shows to be more profitable
than mFRR provision, but is also more consequential for temperature deviations
in the freezer.

###### Index Terms:

Demand-side flexibility, thermostatically controlled loads, supermarket
freezer, manual frequency restoration reserves, load shifting.

##  I Introduction

Power consumers are looking to reduce costs by any means due to sky-rocketing
energy and gas prices. In Denmark, supermarkets are especially exposed due to
their high energy consumption, and they are actively looking into initiatives
to reduce their electricity bill. One of these initiatives is the application
of demand-side flexibility, whereby a supermarket can shift consumption in
time when spot (day-ahead) prices are comparatively lower or when the power
system needs up-regulation services. There are thousands of supermarkets in
Denmark willing to provide flexibility, and their freezers are especially
suited for this as they constitute a big share of their energy consumption.

IBM is currently developing a software platform, called Flex Platform, aiming
to harness flexibility from industrial and commercial consumers for bidding in
ancillary service markets via a Balance Responsible Party (BRP). One large
supermarket chain in Denmark has already signed up for the Flex Platform, and
another supermarket chain is currently exploring the possibility of joining.
While these supermarket chains are keen on the provision of flexibility, the
estimation of their monetary benefit still remains an open question. We aim to
address it by considering two potential ways for supermarkets to exploit their
consumption flexibility: (i) load shifting, and (ii) the provision of manual
Frequency Restoration Reserve (mFRR) services. We develop a stochastic
optimization tool capturing the temperature dynamics and therefore consumption
flexibility of supermarkets, and compute the monetary benefit of flexibility
implemented via load shifting or mFRR service provision. We use two scenario
generation strategies: one using a five-day lookback on spot and balancing
market prices in Denmark, and the other one using those price data in 2021. We
then systematically compare these two strategies ex-post via an out-of-sample
simulation based on price data in 2022.

###  I-A Dilemma: Load Shifting Versus mFRR

For power consumers who are actively looking into cost saving solutions, it is
important to estimate the monetary benefit of flexibility. In Nord Pool, spot
prices are announced at 2pm for the next day, so supermarkets who purchase
power through a retailer at spot prices, have a chance to shift their
consumption plan for the next day, accordingly. From their perspective, it is
natural to first look into load shifting as the most obvious and
straightforward way to utilize flexibility. Here, load shifting simply refers
to shifting consumption to another point in time, typically when spot prices
are lower.

Nevertheless, load shifting, although might be attractive individually, is not
necessarily in favor of the power system. If a significant number of consumers
start to act similarly, it could be detrimental to the power system by moving
the peak consumption to a different time of the day. In this case, spot prices
are no longer reflecting supply-demand equilibrium, and the demand
distribution over the day. Intriguingly, many industrial and residential
consumers have already started shifting demand to some degree, as a response
to the energy crisis. As a result, load shifting might provide short-term
monetary benefit, but it is not necessarily helping the power system in the
long run.

As an alternative option which is certainly in favor of the power system,
consumers could choose to provide flexibility by offering ancillary services
through an aggregator, for which they would be compensated. Here, we focus on
frequency-supporting services, and in particular mFRR, and left a potential
extension by considering more services for the future work. This service is a
slow-responding reserve used to stabilize frequency in the power system after
fast frequency reserves are depleted. The market for mFRR is usually operated
by the national Transmission System Operator (TSO), which is Energinet in
Denmark. Note that mFRR resembles load shifting in the sense that it is the
largest energy reserve operated by the TSO, and can be provided by consumers
via shifting their consumption in time. Likewise, generators provide mFRR by
increasing their power generation.

In order for a supermarket freezer to deliver mFRR, it must be part of a
larger portfolio with potentially many small assets to meet the minimum bid
size requirement for entering the mFRR market. In Denmark, this minimum bid
size is currently 5 MW, but it is expected to reduce to 1 MW in 2024.
Furthermore, there might be additional synergy effects such as increased
operational flexibility in the control response, rebound, and temperature
deviations [1].

For simplicity, we assume that one freezer delivering mFRR can be scaled to
many freezers in the same way. Therefore, demand response for consumers
participating in the mFRR market versus load shifting will look similar. From
the perspective of consumers, it is of great interest to explore the monetary
benefit of mFRR versus load shifting. Similarly, it is important for TSOs to
understand the incentive of consumers to provide mFRR services opposed to
simply shifting load.

###  I-B Research questions and literature review

This work aims to investigate and answer the following research questions: (i)
how can the flexibility of a supermarket freezer be characterized? (ii) how
much monetary benefit can be earned by the mFRR provision opposed to load
shifting? and finally (iii) what are the pros and cons of participating in the
mFRR market versus load shifting? Some of these questions have been looked at
in the literature, so we provide a literature review to pinpoint where new
work is required, and how our work differs from the current literature.

The approach taken in each study in the literature is very much dependent on
the perspective from which flexibility is utilized. Often, full knowledge and
control of all assets is assumed. A prevalent implicit assumption is that the
aggregator and the trading entity are the same, with an exclusive business
relationship to flexible consumers [2]. The incentives and investments
required for delivering flexibility for the flexible consumer are often
overlooked as well. Many also assume an idealized market mechanism or simply
propose a new mechanism for trading flexibility.

In [3] and [4], a complete white-box model of a supermarket refrigeration
system is presented and validated against real-life data. It is also shown how
such a system can provide demand response. This work serves as a benchmark for
any grey-box model of a supermarket refrigeration system. However, such an
approach is hard to scale and requires complete knowledge of every
refrigeration system. In [5], a second-order model is used to model the food
and air temperature in a freezer in an experimental setting. It is shown how
food temperature has much slower dynamics than air temperature. In [6], a
simple first-order virtual battery model for a Thermostatically Controlled
Load (TCL) is presented. This model also constitutes the starting point for
modeling temperature dynamics in this work, as it has an intuitive
interpretation. A similar bucket model is introduced in [7].

In [8], a Mixed-Integer Linear Programming (MILP) problem is developed to
address up- and down-regulation hours from a baseline consumption. It is
assumed that the energy for down-regulation is equal to the energy not
consumed when up-regulating. A similar approach is taken in our work, except
we will define the use of down-regulation until the temperature state is
(approximately) back to its setpoint.

In [9], [10], [11], and [11], residential air condition units are modeled
using up- and down-regulation blocks characterizing flexibility. The blocks
are obtained from grey-box models of households [12]. The authors show how
such a block formulation can be solved using exact and stochastic dynamic
programming in the context of peak shaving demand for a utility in the US. The
demand-side flexibility then functions as a hedge towards extreme electricity
prices. Although such a use case has not been relevant in Denmark yet, the
estimation of flexibility blocks is a novel idea as it avoids having to
explicitly integrate a physical model into an optimization problem. However, a
tradeoff is that the curse of dimensionality quickly makes dynamic programming
computationally intractable for a large portfolio of heterogeneous demand-side
assets. Furthermore, it might be difficult to assure the Markov property [13]
when reformulating the problem even slightly. Flexibility blocks are also used
in [14] for an offering strategy. Here, the flexibility blocks are derived
from measurements of residential appliances in the Ecogrid 2.0 project [15].

In [16] and [17], it is described how pre-defined asset characteristics allow
aggregators to utilize and contract their flexibility through an interface.
This assumes that such characteristics are fully known beforehand. However,
this type of information flexibility interface is still useful for aggregators
when engaging with consumers with well-known assets at the point where
aggregation and penetration of demand-side flexibility is a mature and
prevalent business.

In [18], a portfolio of residential heat pumps is modeled using a linear
first-order model. Individual units are lumped together as essentially one
big, aggregated heat pump as described in [19]. The first-order model is a
grey-box model incorporating the physics, i.e., the temperature dynamics in a
household. Such an approach will be used in this work as well. The authors
also consider how the portfolio of heat pumps can be used for load shifting
and real-time bidding in the balancing market, but do not consider mFRR
service provision. However, they assume that all the assets can be controlled
continuously in a model predictive control setting. Temperature deviations are
minimized directly in the objective function using the integrated error of
total deviations. This introduces a trade-off parameter in the objective
function that must be tuned to weigh the economic value versus the temperature
deviation compared to the setpoint.

###  I-C Our contribution and outline

To the best of our knowledge, no work has investigated the monetary value and
incentives for providing mFRR services versus load shifting, although both
have been studied extensively. This work aims to fill that gap by taking a
holistic view on the incentives from the flexible consumer’s perspective.
Furthermore, when stating the objective functions for mFRR, most studies have
assumed a simplified market structure, and neglect the fact that the
aggregator and the consumer are not necessarily the same entity. We discuss
the consequences of this assumption in relation to the market structure in
Denmark [2], and in the context of mFRR and load shifting.

We provide a realistic model formulation of the mFRR bidding, taking into
account the sequence of decisions needed to be made. We propose two solution
strategies for mFRR bidding: (i) A simple one learned using the Alternating
Direction Method of Multipliers (ADMM) on 2021 price data which can be applied
for all of 2022, and (ii) A dynamic one which computes a new policy every day
by looking at spot prices for the past five days. Real data of a supermarket
freezer in Denmark is used.

The rest of the paper is organized as follows. Section II describes modeling a
supermarket freezer as a TCL. Thereafter, mFRR and load shifting are described
in detail and their objective functions are stated. Section III presents the
proposed stochastic optimization problem. Section IV provides results for a
case study. Section V concludes the paper. Finally, an Appendix provides the
full model formulation.

##  II Describing flexibility of a TCL

There are several ways to monetize flexibility from TCLs. In this section, we
focus on mFRR and load shifting. First, we describe how to mathematically
model a TCL as a flexible resource. Second, we describe how to monetize the
flexibility from TCLs for mFRR and load shifting, and provide the objective
functions in both cases. For mFRR, the objective function includes all costs
and revenues for the BRP, while for load shifting, the objective function only
includes the flexible demand’s perspective. This approach explicitly shows the
situation where flexible demand is activated without including the BRP, as
this is more realistic.

###  II-A Modeling TCL as a flexible resource

TCLs are characterized by being controlled such that the temperature is kept
at a specified setpoint. Examples include heat pumps, freezers, air condition
units, etc. They constitute an important part of demand-side flexibility due
to inherent thermal inertia of such temperature-driven systems [6].

We focus on freezers only, which are a very common type of TCLs. Specifically,
we focus on a single freezer display in a Danish supermarket. Freezers are
characterized by a large thermal inertia due to the frozen food, which makes
them suitable for flexibility. On the other hand, there is a risk of food
degradation when utilizing flexibility. Therefore, it is important to model
the temperature dynamics in the freezer for a realistic and risk-aware
estimation of its flexibility.

The rest of the section is organized as follows. First, we visualize the
measurements from a real supermarket freezer. Second, we introduce a second-
order grey-box model that characterizes the supermarket freezer. Third, we
validate the second-order model and show how it can be used to simulate demand
response from a freezer.

####  II-A1 Supermarket freezer description

We use real data from a single freezer operating in a large Danish supermarket
as a case study. In Fig. 1, the top plot shows the 15-minute average air
temperature of the freezer, whereas the middle plot depicts the opening degree
of the valve. Temperature fluctuates around its setpoint at -18 ∘C with the
exception of hours 7 and 8, where defrosting is scheduled. While defrosting, a
heating element is briefly turned on, and the expansion valve is closed, such
that the flow of refrigerant stops. Afterwards, while recovering the
temperature, the expansion valve is fully opened. The electric power of the
variable-speed compressor rack, scaled to one freezer is shown in the bottom
plot. Since temperature dynamics are similar for all freezers, homogeneity is
assumed. Hence, an equal consumption level is assumed for all freezers.

Consumption is highest during opening hours, and it is lowest during closing
hours. During opening hours, food is being replaced and customers open the
display case constantly. Furthermore, most supermarkets put additional
insulation on the display cases during closing hours which reduces thermal
losses. For these reasons, there are effectively two regimes for a supermarket
freezer plus a short defrosting regime.

Figure 1: Top: temperature of a single freezer in a supermarket. Middle:
opening degree (OD) of the freezer expansion valve. Bottom: electric power of
the compressor rack feeding a single freezer.

####  II-A2 Thermal modeling of freezer

In [6], it is described how a simple TCL model can be developed. We extend it
to a second-order state-space model that accounts for the thermal mass of the
food, which essentially provides the flexibility in freezers:

$\displaystyle T^{\rm{f}}_{t+1}$ $\displaystyle=T^{\rm{f}}_{t}+dt\
\frac{1}{C^{\rm{f}}}\left(\frac{1}{R^{\rm{cf}}}(T^{\rm{c}}_{t}-T^{\rm{f}}_{t})\right)$
(1) $\displaystyle T^{\rm{c}}_{t+1}$ $\displaystyle=T^{\rm{c}}_{t}+dt\
\frac{1}{C^{\rm{c}}}\Bigl{(}\frac{1}{R^{\rm{cf}}}(T^{\rm{f}}_{t}-T^{\rm{c}}_{t})+\frac{1}{R^{\rm{ci}}}(T^{\rm{i}}_{t}-T^{\rm{c}}_{t})$
$\displaystyle\mspace{50.0mu}-\eta\ OD_{t}\
P_{t}\Bigr{)}+\epsilon\mathbbm{1}^{\rm{df}}_{t},$ (1)

where $T^{\rm{c}}_{t}$ is the measured air temperature in the freezer at time
$t\in\\{1,\ldots,J\\}$, and $T^{\rm{f}}_{t}$ is the food temperature which is
a latent unobserved state. Parameters $C^{\rm{f}}$ and $C^{\rm{c}}$ are the
thermal capacitance of the food and air in the freezer, respectively. In
addition, $R^{\rm{cf}}$ and $R^{\rm{ci}}$ are the thermal resistance between
food and air in the freezer, and air and indoor temperature, respectively.
There is essentially a low-pass filter of the air temperature in the freezer
with the time constant $C^{\rm{f}}$ times $R^{\rm{cf}}$. Parameter $\epsilon$
represents the temperature change when defrosting, whereas
$\mathbbm{1}^{\rm{df}}_{t}$ is the indicator function for when defrosting
happens. Parameter $R^{\rm{ci}}$ is either one of two values,
$R^{\rm{ci},\text{day}}$ and $R^{\rm{ci},\text{night}}$, to capture the
differences between opening hours (6am to 10pm) and closing hours (10pm to
6am). The opening degree $OD_{t}$, indoor air temperature $T^{\rm{i}}_{t}$,
and power $P_{t}$ are exogenous inputs as the opening degree is assumed to be
fixed during a demand-response event, while only $P_{t}$ is controllable, and
the indoor air temperature is unaffected by the freezers. Parameter $\eta$ is
the compressor efficiency. The model is discretized with a time step of 15
minutes, therefore $dt=0.25$ and $J=96$.

####  II-A3 Model validation

Using the R library of CTSM-R [20], all parameters in (1) have been estimated
as given in Table I. Notice that the thermal capacitance of the air in the
freezer is significantly smaller than the thermal capacitance of the food,
indicating that the food temperature changes comparatively slower. The thermal
resistance between the food and air inside the freezer, $R^{\rm{cf}}$, is also
significantly smaller than the thermal resistance between the air in the
freezer and the indoor temperature in the supermarket, $R^{\rm{ci}}$, both
during the day and night. This makes sense as the lid acts as a physical
barrier insulating the freezer. Furthermore, the thermal resistance to the
indoor air temperature is higher during the night, which means that less power
is needed, as observed in Fig. 1.

TABLE I: Parameter Estimates of (1). Parameter | Value | Unit | Parameter | Value | Unit
---|---|---|---|---|---
$C^{\rm{f}}$ | 6.552 | kWh/∘C | $R^{\rm{ci},\text{night}}$ | 61.25 |  ∘C/kW
$C^{\rm{c}}$ | 0.077 | kWh/∘C | $\eta$ | 1.561 |
$R^{\rm{cf}}$ | 5.010 |  ∘C/kW | $\epsilon$ | 3.372 |  ∘C/h
$R^{\rm{ci},\text{day}}$ | 41.05 |  ∘C/kW |  |  |

The one-step residuals for the air temperature should ideally resemble white
noise in order for a model to capture all dynamics observed in the data [21].
Fig. 2 shows the auto-correlation and cumulative periodogram of the residuals.
The autocorrelation shows two significant lags for lag two and seven, but
looks satisfactory otherwise. Likewise, it seems the model is able to capture
most dynamics at all frequencies. Since there are only $J=96$ time steps, the
defrosting period can result in relatively large residuals as it is difficult
to capture such fast transient dynamics. Since up-regulation is not allowed
during defrosting, the residuals are not too important during that period. To
decrease their effect in the parameter estimation procedure, the term
$\epsilon\mathbbm{1}^{\rm{df}}_{t}$ was added to (1).

Figure 2:  Validation of the state-space model in (1). Left: auto-correlation
function (acf) of the residuals. Right: cumulated periodogram of the
residuals.

Furthermore, Fig. 3 (left) shows a 24-hour simulation of (1). It is observed
that the simulation is reasonable and closely follows the measured air
temperature, although a slight bias is observed as the simulated temperature
is a bit higher. Such a visual validation is important because the model will
be embedded later in an optimization model in Section III.

Ideally, the validation of (1) would also include real measurements from the
air and food temperature in a freezer during demand-response events. However,
by adhering to the fundamental physics governing the temperature dynamics as
shown, the model is trusted to be accurate during demand-response events too.
It brings an intuitive interpretation which can be used to understand the
impact on temperature during a demand-response event. In Fig. 3 (right), such
a simulated example of a demand-response event is shown. It can be observed
how the air temperature increases when the power is turned off, and how it
decreases when the power is turned back on. The food temperature is much more
stable and only changes slightly, as expected. The rebound occurs until the
food temperature is back to its normal value.

Figure 3:  Left: Simulation of food and air temperatures using (1) with
parameters in Table I, and its comparison to the measured air temperature.
Right: Simulation where power is turned off for two hours with a subsequent
rebound at the nominal power until the food temperature is back to its normal
value.

Notational remark: Hereafter, $t\in\\{1,\ldots,J=96\\}$ represents the index
for 15-minutes time steps, whereas $h\in\\{1,\ldots,24\\}$ denotes the index
for hours.

###  II-B mFRR

Fig. 4 shows the timeline of the mFRR market in Denmark. There is a market for
up-regulation only. One should note that BRPs can choose not to bid in the
reserve market and only bid in the real time for up- and down-regulation. We
consider a case wherein the flexible demand delivers reservation through a
BRP, since the payment is received for both reservation and activation. All
prices considered are for the bidding zone DK2 (eastern Denmark).

[width=]tex/figures/timeline_mfrr_tikz

Figure 4: Timeline of the Danish mFRR market. Symbols $\bm{p}$ denote power
(in MW), whereas symbols $\bm{\lambda}$ refer to price (in DKK/MW).

According to Fig. 4, at 9:30am of day $\rm{D}$-1, the BRP submits reserve
capacity bids $p_{h}^{\rm{r},\uparrow}$ $\forall{h}$ to the mFRR market for
the next day $\rm{D}$. If accepted, she is paid at the reservation price
$\lambda_{h}^{\rm{r}}$. Therefore, she earns
$\lambda_{h}^{\rm{r}}p_{h}^{\rm{r},\uparrow}$. This happens before the day-
ahead (spot) market clearing at noon for which the BRP buys energy for her
expected demand $P_{h}^{\rm{Base}}$ at the spot price $\lambda_{h}^{\rm{s}}$.
After that, at 5pm, a regulating power bid, $\lambda_{h}^{\rm{bid}}$, must be
submitted by the BRP for each hour of day $\rm{D}$, where
$p_{h}^{\rm{r},\uparrow}>0$ [22]. In hour $h$ of day $\rm{D}$, the reserves
are activated if the following three conditions hold: (i) there is a
reservation accepted, i.e., $p_{h}^{\rm{r},\uparrow}>0$, (ii) the balancing
price is higher than the spot price, i.e.,
$\lambda_{h}^{\rm{b}}>\lambda_{h}^{\rm{s}}$, and the regulating power bid is
lower than or equal to the balancing price minus the spot price, i.e.,
$\lambda_{h}^{\rm{bid}}\leq\lambda_{h}^{\rm{b}}-\lambda_{h}^{\rm{s}}$. If all
these three conditions are met, the BRP is paid at the balancing price
$\lambda_{h}^{\rm{b}}$ times her actual up-regulation
$p_{h}^{\rm{b},\uparrow}$. The BRP may incur a cost due to any subsequent
rebound $p_{h}^{\rm{b},\downarrow}$. Furthermore, the BRP should pay at the
penalty price $\lambda^{\rm{p}}$ for
$s_{h}=\text{max}\\{0,p_{h}^{\rm{r},\uparrow}-p_{h}^{\rm{b},\uparrow}$}, i.e.,
if she fails delivering her promised reserve. In reality,
$p_{h}^{\rm{b},\uparrow}$ is determined by the TSO, but for simplicity, we
assume a bid is always activated in its entirety. Accordingly, the objective
function of the BRP maximizing the flexibility value over a day (via mFRR
provision) is

$\displaystyle-\underbrace{\sum_{h=1}^{24}\lambda^{\rm{s}}_{h}P^{\rm{Base}}_{h}}_{\textrm{Energy
cost}}+\underbrace{\sum_{h=1}^{24}\lambda_{h}^{\rm{r}}p^{\rm{r},\uparrow}_{h}}_{\textrm{Reservation
payment}}$
$\displaystyle\quad\quad+\underbrace{\sum_{h=1}^{24}\lambda_{h}^{\rm{b}}p^{\rm{b},\uparrow}_{h}}_{\textrm{Activation
payment}}-\underbrace{\sum_{h=1}^{24}\lambda_{h}^{\rm{b}}p^{\rm{b},\downarrow}_{h}}_{\textrm{Rebound
cost}}-\underbrace{\sum_{h=1}^{24}\lambda^{\rm{p}}s_{h}}_{\textrm{Penalty
cost}}.$ (2)

Note that the first term of (II-B), i.e., the energy cost, contains
$\lambda^{\rm{s}}_{h}$ and $P^{\rm{Base}}_{h}$, which are both parameters, and
therefore, can be omitted from the objective function. However, we include it
here for the completeness.

###  II-C Load shifting

Another option for utilizing flexibility is to shift the load to a different
time according to the spot market prices which are known already 12-36 hours
in advance. Then, it is simply a matter of consuming in low-price hours. For a
TCL, there are additional constraints to how energy can be shifted due to the
rebound effect. First, there can be temperature constraints which will result
in less energy being shifted. Second, the rebound must happen immediately
after reducing power consumption. Otherwise, the temperature deviation becomes
too big for too long. The cost saving from load shifting is directly
proportional to the volume and price difference between the baseline load and
shifted load as given by

$\sum_{h=1}^{24}\lambda^{\rm{s}}_{h}\big{(}p^{\rm{Base}}_{h}-p_{h}\big{)},$
(3)

where $p_{h}$ is the power profile when the load is shifted, whereas
$p^{\rm{Base}}_{h}$ is the baseline in the day-ahead stage.

Since the load shifting action only occurs after the day-ahead market clearing
(cf. Fig. 4), the BRP has already bought
$\lambda^{\rm{s}}_{h}p^{\rm{Base}}_{h}$ and any deviation results in an
imbalance for the BRP. In this work, we look at the case where the flexible
demand acts selfishly and excludes the BRP from her load shifting action.
Therefore, the objective function for the flexible demand simply reduces to
$\sum_{h=1}^{24}\lambda_{h}^{\rm{s}}p_{h}$.

##  III Optimization model and solution strategy

This section presents the optimization problem. The rest of this section is
organized as follows. First, the time sequence of decisions for participation
in the mFRR market is presented. Second, we explain how scenarios for price
data are generated. Third, we present the model formulation. Fourth, we
discuss how the bidding policy is implemented. Lastly, we show how a scenario
decomposition method with an ADMM strategy is used to solve the optimization
model.

###  III-A Time sequence for decision making

Fig. 5 shows the stages for making decisions in the mFRR market. In the first
stage in day $\rm{D}$-1, the BRP makes a reservation bid decision
$p_{h}^{\rm{r},\uparrow}$ for every hour $h$ of day $\rm{D}$, while being
uncertain about input parameters in the next stages, including day-ahead
market prices $\lambda_{h}^{\rm{s}}$ and balancing market prices
$\lambda_{h}^{\rm{b}}$. We assume the BRP is not uncertain about mFRR market
prices $\lambda_{h}^{\rm{r}}$.

For simplicity and avoid the need for developing a multi-stage stochastic
program, we merge the second, third, and fourth stages in Fig. 5 as one, and
call it the second stage. By this, $p_{h}^{\rm{r},\uparrow}$ is the first-
stage variable, whereas the second stage variables, indexed by scenario
$\omega$, include the regulating power bid $\lambda_{h,\omega}^{\rm{bid}}$ and
the set of real-time variables
$\Gamma_{h,\omega}=\\{p_{h,\omega},p_{h,\omega}^{\rm{b},\uparrow}$,
$p_{h,\omega}^{\rm{b},\downarrow}$, $s_{h,\omega}$, $T_{h,\omega}^{\rm{c}}$,
$T_{h,\omega}^{\rm{f}}$, $T_{h,\omega}^{\rm{c,\text{Base}}}$,
$T_{h,\omega}^{\rm{f,\text{Base}}}$, $\phi_{h,\omega}$,
$g_{h,\omega},u^{\uparrow}_{h,\omega},u^{\downarrow}_{h,\omega},y^{\uparrow}_{h,\omega},y^{\downarrow}_{h,\omega},z^{\uparrow}_{h,\omega},z^{\downarrow}_{h,\omega}\\}$.
This set contains the real-time power activated, as well as auxiliary
variables for identifying up- and down-regulation, temperature dynamics, and
when to deliver up-regulation according to the bid and prices. See the
Appendix for a detailed description. Note that down-regulation refers to the
rebound effect of the freezer.

[width=]tex/figures/timeline_mfrr_variables_tikz

Figure 5: Variables related to the mFRR bidding decisions.

###  III-B Scenario generation

To make decisions on reservation capacity $p_{h}^{\rm{r},\uparrow}$ and
regulating power bid $\lambda_{h,\omega}^{\rm{bid}}$, we generate a set of
scenarios $\omega\in\Omega$ for spot and balancing price data., i.e.,
$\lambda_{h,\omega}^{\rm{s}}$ and $\lambda_{h,\omega}^{\rm{b}}$. Each scenario
contains spot and balancing price data for the entire day in question. The
number of scenarios is $|\Omega|$, and we assume all scenarios are
equiprobable. We refer to them as in-sample scenarios. We use two different
strategies for in-sample scenario generation: (i) considering historical spot
and balancing prices in DK2 in 2021. (ii) considering prices of the most
recent five days (lookback strategy). In both scenario generation strategies,
balancing prices $\lambda_{h,\omega}^{\rm{b}}$ are sampled in the following
way: First, an integer $v$ is sampled uniformly from $\\{0,\ldots,24\\}$ which
represents the total number of up-regulation hours in a day. We then sample
spot and balancing price differentials from a single day within the set of all
days where up-regulation happened $v$ times. This constitutes one scenario and
is repeated $|\Omega|$ times. In this way, days where up-regulation happened
are essentially up-sampled, and the model learns more when up-regulation
happens than otherwise. The in-sample results are systematically compared
against the same set of unseen Out-of-Sample (OOS) scenarios, which are DK2
prices for 2022.

For load shifting, the solution approach is simply to solve a deterministic
optimization problem for the next day, since the day-ahead (spot) market
clearing happens in advance, and thereby the spot market prices are known.

###  III-C Model formulation

Recall that (II-B) gives the objective function of a BRP participating in the
mFRR market, but in a deterministic setup. The optimization problem of such a
BRP in a two-stage stochastic programming format is

where bold symbols represent vectors. For example, the vector
$\bm{p}^{\rm{r},\uparrow}$ includes $p_{h}^{\rm{r},\uparrow}\ \forall{h}$.
Parameter $\pi_{\omega}$ is the probability assigned to in-sample scenario
$\omega$. Here, we present optimization (4) in a compact form using functions
$f(.)$, $g(.)$, and $h(.)$. The detailed formulation is given in the Appendix.
Constraint (LABEL:P1:eq2) includes all power and activation related limits,
whereas (LABEL:P1:eq3) models temperature dynamics in the freezer. Constraints
(LABEL:P1:eq4)-(LABEL:P1:eq5) declare continuous and binary variables. Note
that $\mathbb{R}$ and $\mathbb{R}^{+}$ denote free and non-negative real
numbers, respectively. The optimization model (4) is a MILP problem.

For optimal decision making for load shifting, (4) is simplified by removing
scenarios as well as reservation and bid constraints, and replacing the
objective function (LABEL:P1:eq1) by minimizing the total power purchase cost
$\bm{\lambda}^{\rm{s}}\bm{p}$.

###  III-D Regulating power bidding implementation

Recall from Fig. 5 that we have merged three stages, by which the regulating
power bidding decisions $\lambda_{h,\omega}^{\rm{bid}}$ and real-time
operational decisions $\Gamma_{h,\omega}$ become second-stage variables. This
is the reason both set of variables $\lambda_{h,\omega}^{\rm{bid}}$ and
$\Gamma_{h,\omega}$ are similarly indexed by $\omega$, while in reality, the
decision $\lambda_{h,\omega}^{\rm{bid}}$ should be made before
$\Gamma_{h,\omega}$. This also challenges the ex-post OOS simulation. To
resolve it, we use a learning policy, such that we replace the scenario-
indexed variable $\lambda_{h,\omega}^{\rm{bid}}$ in (4) by $\alpha
q(\lambda_{h,\omega}^{\rm{s}})+\beta$, where $q(.)$ is an arbitrarily selected
function. In addition, $\alpha$ and $\beta$ are non-negative first-stage
variables (they are not indexed by $\omega$). This replacement shrinks the
degree of freedom for the BRP compared to (4), but makes it more practical to
be used. The reason is that, by this trick, the regulating power bidding
decision to be made at 5pm of day $\rm{D}$-1 becomes a first-stage decision,
dependent not only on policies $\alpha$ and $\beta$, but also on uncertain
spot prices $\lambda_{h,\omega}^{\rm{s}}$. In other words, by using the in-
sample scenarios $\omega$, the BRP obtains optimal values for $\alpha$ and
$\beta$ at 9:30am of day $\rm{D}$-1. Then, she waits to see the spot prices at
noon of day $\rm{D}$-1, and submits her regulating power bids at 5pm of day
$\rm{D}$-1. In our simulations, we found out that the selection of
$q(\lambda_{h,\omega}^{\rm{s}})$ as the difference of spot prices in
subsequent hours works comparatively more satisfactory ex-post, as it suits
better to accommodate the rebound effect. By this, the BRP prefers to be
activated when the spot price during rebound (down-regulation) is
comparatively lower. Therefore, we replace $\lambda_{h,\omega}^{\rm{bid}}$ in
(4) by

$\displaystyle\alpha(\lambda_{h+1,\omega}^{\rm{s}}-\lambda_{h,\omega}^{\rm{s}})+\lambda_{h,\omega}^{\rm{s}}+\beta,\
\forall{\omega},\forall{h}\in\\{1,\ldots,23\\}.$ (4)

Another implementation challenge is to enforce price conditions under which
the mFRR reservation is activated. Recall from Section II-B, for the
activation at hour $h$ under scenario $\omega$, it is necessary to hold
$\lambda_{h,\omega}^{\rm{bid}}\leq\lambda_{h,\omega}^{\rm{b}}-\lambda_{h,\omega}^{\rm{s}}$
and $\lambda_{h,\omega}^{\rm{b}}>\lambda_{h,\omega}^{\rm{s}}$. These
conditions can be equivalently enforced as

$p^{\rm{b},\uparrow}_{h,\omega}+s_{h,\omega}\geq
p^{\rm{r},\uparrow}_{h}\mathbbm{1}^{\big{(}\lambda_{h,\omega}^{\rm{bid}}\leq\lambda_{h,\omega}^{\rm{b}}-\lambda_{h,\omega}^{\rm{s}}\
\text{and}\ \lambda^{\rm{b}}_{h,\omega}>\lambda^{\rm{s}}_{h,\omega}\big{)}},$
(4)

where $\mathbbm{1}^{(.)}$ is 1 if conditions (.) are met, otherwise it is
zero. If the non-negative slack variable $s_{h,\omega}$ takes a non-zero
value, it shows that the BRP fails in the activation stage, and therefore will
be penalized. The challenge is that (4) makes a condition on variable
$\lambda_{h,\omega}^{\rm{bid}}$, or equivalently on $\alpha$ and $\beta$ as
defined in (4). This makes (4) non-linear. To linearize it, we use the
McCormick relaxation technique [23] and define auxiliary variables
$\phi_{h,\omega}\in\mathbb{R}^{+}$ and $g_{h,\omega}\in\\{0,1\\}$. By this, we
replace (4) for every hour $h$ and scenario $\omega$ by a set of mixed-integer
linear constraints as

where $M$ is a large enough positive constant. Constraint
(LABEL:con_bid:subeq1) ensures that $g_{h,\omega}=1$ when
$\lambda_{h,\omega}^{\rm{b}}-\lambda_{h,\omega}^{\rm{s}}\geq\lambda^{\rm{bid}}_{h,\omega}$.
Otherwise, it sets $g_{h,\omega}=0$. Constraints
(LABEL:con_bid:subeq3)-(LABEL:con_bid:subeq4) set the up-regulation equal to
$\phi_{h,\omega}$ (or incurs a penalty through $s_{h,\omega}$) if there is an
up-regulation event in the system, i.e., if
$\mathbbm{1}^{\lambda^{\rm{b}}_{h,\omega}>\lambda^{\rm{s}}_{h,\omega}}=1$.
Constraint (LABEL:con_bid:subeq5) enforces $\phi_{h,\omega}=0$ when
$g_{h,\omega}=0$, implying the balancing price minus the spot price is smaller
than the power regulating bid. Constraint (LABEL:con_bid:subeq7) ensures that
$\phi_{h,\omega}$ is equal to the reservation capacity
$p^{\rm{r},\uparrow}_{h}$ whenever $g_{h,\omega}=1$, i.e., if
$\lambda_{h,\omega}^{\rm{b}}-\lambda_{h,\omega}^{\rm{s}}\geq\lambda^{\rm{bid}}_{h,\omega}$.
Note that in the final model, $\lambda^{\rm{bid}}_{h,\omega}$ in (4) should be
replaced as defined in (4). The resulting model formulation, which is a MILP,
is given in the Appendix.

###  III-E Scenario decomposition with ADMM

By increasing the number of scenarios, the proposed stochastic program quickly
becomes computationally intractable due to the number of binary variables and
intertemporal constraints. To resolve it, we apply a scenario decomposition
approach built upon the ADMM algorithm [24]. To do so, we relax non-
anticipativity constraints by adding a scenario index to the first-stage
variables, i.e.,

$p_{h}^{\rm{r},\uparrow}\rightarrow p^{\rm{r},\uparrow}_{h,\omega},\
\alpha\rightarrow\alpha_{\omega},\ \beta\rightarrow\beta_{\omega}.$ (4)

This decomposes the original problem to a set of deterministic sub-problems,
one per scenario, to be solved in parallel. A quadratic regularizer is added
to the objective function of every subproblem, making it a mixed-integer
quadratic program. The ADMM algorithm is iterative. The convergence happens
when we achieve a consensus over sub-problems on first-stage variables. Due to
having binary variables in the original problem, this ADMM algorithm is
eventually a heuristic [25], i.e., it may not converge to optimality. However,
we have observed in a case with a limited number of scenarios for which we can
also solve the original MILP problem directly, the proposed ADMM exhibits a
satisfactory performance.

##  IV Numerical Results and Discussion

We consider five models to calculate the operational (energy) cost of a single
freezer, including three models where the freezer provides mFRR services, one
model where the freezer shifts the load in response to spot prices, and the
last model, the so-called base cost, where the freezer neither shifts the load
nor provides mFRR services. The data availability for these five models is
summarized in Table II. All source codes and relevant data are publicly shared
in [26].

###  IV-A Load shifting vs mFRR: Which one is more appealing?

TABLE II: Data availability for each model. Model | Curve(s) in Fig. 6 | $\bm{\lambda}^{\rm{r},\uparrow}$ | $\bm{\lambda}^{\rm{s}}$ | $\bm{\lambda}^{\rm{b}}$
---|---|---|---|---
mFRR oracle | Dashed black | $\checkmark$ | $\checkmark$ | $\checkmark$
mFRR cases | Blue and red | $\checkmark$ | $\ast$ | $\ast$
Load shifting | Yellow | N/A | $\checkmark$ | N/A
Base cost | Solid black | N/A | $\checkmark$ | N/A
$\checkmark$: True (realized) data are available.
$\ast$: Forecast data (in form of scenarios) are available.
N/A: Not applicable.

For all five models, Fig. 6 shows the out-of-sample cumulative operational
cost of the single freezer during the first nine months of 2022. Recall that
all these models have been compared fairly by an out-of-sample simulation
against identical scenarios, i.e., real spot and balancing market prices from
2022. The base cost (solid black) has only access to spot price data, and
leads to the highest cost in Fig. 6. Blue and red curves (mFRR cases), both
below the base cost curve, correspond to the cases where the freezer provides
mFRR services, and uses scenarios to model spot and balancing market price
uncertainties. Their difference comes from in-sample scenarios used: while the
red curve uses five equiprobable scenarios coming from the most recent
historical data (lookback strategy), the blue curve uses 50 equiprobable
scenarios coming from 2021. The yellow curve shows the reduced cost of the
freezer due to load shifting. Finally, the dotted black curve (lowest in Fig.
6) provides an oracle (ideal benchmark), where the freezer perfectly knows the
spot and balancing market prices.

An interesting observation is that, in comparison to the base cost, load
shifting (yellow curve) bring a 13.9% cost reduction, whereas the mFRR
provision reduces the cost by 11.6% (red curve) and 10.1% (blue curve). Note
that the cost savings reported for red and blue curves are not far from the
mFRR oracle, showing that the scenarios used are sufficiently adequate. While
this observation might be appealing to the freezer as load shifting requires a
comparatively simpler decision-making process than the mFRR provision, it is
not necessarily a desirable outcome for the power system. The mFRR provision
helps the system keep the supply-demand balance, while load shifting as a
response to spot prices is not necessarily helping the system, and even in the
worse case, the spot prices fixed before load shifting may no longer represent
the supply-demand equilibrium. This calls system operators and regulators for
potential changes, e.g., market redesign, to make the mFRR provision more
attractable.

Figure 6: Out-of-sample cumulative operational cost for the freezer during the
first nine months of 2022.

In the case of load shifting, it is worth mentioning that the settlement cost
of the BRP is ignored — this reflects the self-interested flexible loads.
While the flexible load reduces her cost by load shifting, the corresponding
BRP may incur a settlement cost due to the resulting imbalance between the
true consumption and the day-ahead schedule. In practice, the BRP would have
to pay for such an imbalance or at least buy/sell energy in the intra-day
market to adjust her day-ahead schedule and consider potential load shifting
actions by flexible loads. However, loads are not necessarily aware of those
costs for the BRP. Hence, they may have a strong incentive to exploit their
flexibility for load shifting. Some larger flexible consumers, such as
industrial and commercial loads, might not be exposed fully to spot prices,
and for those consumers, load shifting could be less profitable. However, they
may still have an incentive to change their deal with the BRP to get a full
exposure.

In the case of mFRR provision, it is assumed that the entire revenue from
reservation and activation go to the flexible load. This neglects the fact
that, in practice, the BRP requires a share of the revenue. Furthermore, there
might be an aggregator or a technology provider who facilitates the
aggregation and communication of the flexibility. In such a case, they would
also request a share of the revenue. This can potentially reduce the revenue
of flexible loads by the mFRR provision.

Figure 7: Example of load shifting in a representative day (an in-sample
scenario). Top: Baseline consumption of the freezer and the power profile
after load shifting. Middle: Air and food temperature dynamics. Bottom: Spot
market prices. Figure 8: Example of mFRR provision in a representative day (an
in-sample scenario). Top: Reservation capacities and baseline power of the
freezer. Upper middle: Air and food temperature dynamics. Lower middle: mFRR
activation in this specific scenario, i.e., when
$\lambda_{h}^{\rm{bid}}\leq\lambda_{h}^{\rm{b}}-\lambda_{h}^{\rm{s}}$,
$\lambda_{h}^{\rm{b}}>\lambda_{h}^{\rm{s}}$, and $p^{\rm{r},\uparrow}_{h}>0$.
Bottom: Spot and balancing market price as well as regulating power bids in
this specific scenario.

###  IV-B Technical results

Fig. 7 shows technical results for the case of load shifting in a
representative day. For given spot prices (bottom plot), it is evident that
the freezer shifts the load to low-price hours, especially in the later hours
of the day (top plot). However, this has a remarkable effect on the air and
food temperatures with large deviations from the normal set-point (middle
plot).

Fig. 8 presents technical results for the case of mFRR provision in a
representative day, whose spot and balancing market prices are given in the
lower plot. The top plot shows the baseline power profile
$P^{\text{Base}}_{h}$ and the mFRR reservation $p^{\rm{r},\uparrow}_{h}$ sold.
For a better illustration, we have plotted
$P^{\text{Base}}_{h}-p^{\rm{r},\uparrow}_{h}$, i.e., the power profile of the
freezer if the full activation happens during the entire day. The third plot
of Fig. 8 shows, for this specific day (as an in-sample scenario), the freezer
is activated three times, i.e., in hours 6, 9, and 10. In these three hours,
all conditions
$\lambda_{h,\omega}^{\rm{bid}}\leq\lambda_{h,\omega}^{\rm{b}}-\lambda_{h,\omega}^{\rm{s}}$,
$\lambda_{h,\omega}^{\rm{b}}>\lambda_{h,\omega}^{\rm{s}}$, and
$p^{\rm{r},\uparrow}_{h}>0$ hold (see the bottom plot). The rebound happens in
hours 11 to 19. The majority of rebound occurs in hour 19, as the balancing
price is comparatively low. Finally, the second plot shows that the food and
air temperature deviations in the freezer due to mFRR provision are smoother
in comparison to load shifting.

Our second main observation is as follows: Although the cost saving is higher
for load shifting in comparison to mFRR provision, it is directly proportional
to the energy shifted, and therefore, the temperature deviation in the
freezer. This is not the case for mFRR, as the reservation is not always
activated.

##  V Conclusion

We investigated how a supermarket freezer can provide flexibility for mFRR and
load shifting in Denmark, and explored which one provides a greater monetary
incentive for a flexible load. To this end, we used actual data from a Danish
supermarket. This was done by developing a second-order grey-box model of the
temperature dynamics in the freezer with the food temperature as a latent
state. In the state-space form, the model was directly incorporated as
constraints into a two-stage stochastic MILP problem, whose objective is to
maximize the monetary value from the freezer’s flexibility. Two scenario
generation strategies were implemented: one with a five-day lookback strategy
on the spot and balancing market prices in DK2, and the other one based on
price data for those markets in 2021. For mFRR, we used a linear policy, and
then linearized the conditions for activation via the McCormick relaxation
method. For computational ease, we used an ADMM-based scenario decomposition
technique. An out-of-sample evaluation was done on unseen 2022 price data. We
observed that load shifting is more profitable, but has a greater impact on
the air and food temperatures in the freezer as opposed to mFRR that depends
on the system state and bid price for activation.

We made a set of simplifications and assumptions, whose impacts need to be
explored for the future work. The revenue share of BRP and aggregator was not
considered. This may change our finding that load shifting is more financially
appealing to a flexible load. As mentioned earlier, a single freezer or
supermarket must be part of a larger portfolio through an aggregator in order
to participate in the mFRR market. Such an aggregated portfolio has some
issues that are neglected here, such as the baseline estimation for
verification of the demand response, allocation of profits within the
portfolio, and an accurate capacity estimation of the whole portfolio that
bids in the mFRR market. Furthermore, the European mFRR markets will change
from a 60-minute resolution to a 15-minute market in the next few years [27].
This makes it more feasible for TCLs to participate in the mFRR market, given
their sensitivity to large temperature deviations, and the fact that TCLs can
get a passive income when there is no up-regulation need in the power grid.

## Appendix: MILP problem formulation

The two-stage stochastic MILP problem for the freezer to optimally bid in the
mFRR market reads as

The objective function (LABEL:P2:1) maximizes the expected flexibility value
of the freezer. Constraint (LABEL:P2:2) contains (4), representing the
McCormick relaxation of activation conditions. Recall that
$\lambda^{\rm{bid}}_{h,\omega}$ in (4) should be replaced as defined in (4).
Constraint (LABEL:P2:3) declares continuous and binary variables. Aligned with
(1), constraints (LABEL:P2:state-space-1)-(LABEL:P2:state-space-2) are the
state-space model for the food and air temperature dynamics. Similarly,
(LABEL:P2:state-space-3)-(LABEL:P2:state-space-4) include the baseline air
temperature $T^{\rm{c},\text{Base}}_{t}$ and the baseline food temperature
$T^{\rm{f},\text{Base}}_{t}$, and model temperature dynamics for the baseline
power. Recall in case the hour index $h$ runs from 1 to 24, index $t$ runs
from 1 to $J=96$. Constraint (LABEL:power:6) sets the real-time power
consumption $p_{h,\omega}$ equal to the baseline power $P^{\rm{Base}}_{h}$
unless there is up-regulation $p^{\rm{b},\uparrow}_{h,\omega}$ or down-
regulation $p^{\rm{b},\downarrow}_{h,\omega}$. Constraint (LABEL:power:7)
binds the mFRR reservation $p_{h}^{\rm{r},\uparrow}$ to the baseline power.
Constraint (LABEL:power:8) ensures that up-regulation is zero when there is no
need for up-regulation, and at the same time binds it to the reservation
power. Constraint (LABEL:power:9) includes the binary variable
$u^{\uparrow}_{h,\omega}$, indicating whether the freezer is up-regulated in
hour $h$ under scenario $\omega$. This constraint ensures that up-regulation
is zero whenever $u^{\uparrow}_{h,\omega}=0$, and otherwise restricted to the
maximum up-regulation service $P^{\rm{Base}}_{h}-P^{\rm{Min}}$ that can be
provided. Note that $P^{\rm{Min}}$ is the minimum consumption level of the
freezer. Constraint (LABEL:power:10) works similarly for down-regulation. Note
that the binary variable $u^{\downarrow}_{h,\omega}$ indicates whether down-
regulation happens, whereas $P^{\rm{Nom}}$ is the nominal (maximum)
consumption level of the freezer. Constraint (LABEL:power:11) restricts the
power consumption to lie within the minimum and nominal rates. Constraint
(LABEL:power:12) binds the slack variable $s_{h,\omega}$, representing the
service not delivered as promised. Constraint (LABEL:power:14) ensures that
down-regulation is equal to at least 10% of the down-regulation capacity.
Constraint (LABEL:power:15) prohibits any up-regulation when defrosting
occurs. Constraints (LABEL:aux:1)-(LABEL:aux:7) define auxiliary binary
variables $y^{\uparrow}_{h,\omega}$, $y^{\downarrow}_{h,\omega}$,
$z^{\uparrow}_{h,\omega}$, and $z^{\downarrow}_{h,\omega}$, identifying
transitions from/to up-regulation and down-regulation. During all hours with
up-regulation, $y^{\uparrow}_{h,\omega}=1$. In the hour that up-regulation is
stopped, $z^{\uparrow}_{h,\omega}$ is 1. There is a similar definition for
$y^{\downarrow}_{h,\omega}$ and $z^{\downarrow}_{h,\omega}$ related to down-
regulation. See Chapter 5 of [28] for complete details. Constraint
(LABEL:temp:1) restricts the food temperature for the last time period $J$.
Constraints (LABEL:rebound:1)-(LABEL:rebound:4) control the rebound behavior
such that the rebound finishes when the temperature is below the baseline
temperature. Note that $M$ is a sufficiently big positive constant such that
the food temperature is allowed to deviate from the baseline. Also, they
ensure that the rebound happens right after up-regulation. Lastly,
(LABEL:up_reg_first) ensures that up-regulation happens first. This makes
sense since it impossible (or at least difficult) to anticipate potential up-
regulation events in the power system. As such, it does not make sense to pre-
cool (or pre-heat) a TCL in the context of mFRR.

## Acknowledgement

The authors would also like to thank Coop and AK-Centralen for providing the
freezer and power data from a supermarket.

## References

  * [1] S. Koch, J. L. Mathieu, and D. S. Callaway, “Modeling and control of aggregated heterogeneous thermostatically controlled loads for ancillary services,” in _17th Power Systems Computation Conference (PSCC)_ , Stockholm, Sweden, 2011, pp. 1–8.
  * [2] P. A. Gade, T. Skjøtskift, H. W. Bindner, and J. Kazempour, “Ecosystem for demand-side flexibility revisited: The Danish solution,” _The Electricity Journal_ , vol. 35, no. 9, p. 107206, 2022.
  * [3] L. N. Petersen, H. Madsen, and C. Heerup, “ESO2 optimization of supermarket refrigeration systems: Mixed integer MPC and system performance,” 2012. [Online]. Available: https://eudp.dk/files/slutrapporter/3.pdf
  * [4] R. Pedersen, J. Schwensen, S. Sivabalan, C. Corazzol, S. E. Shafiei, K. Vinther, and J. Stoustrup, “Direct control implementation of a refrigeration system in smart grid,” in _2013 American Control Conference_ , Washington, DC, USA, 2013, pp. 3954–3959.
  * [5] R. Pedersen, J. Schwensen, B. Biegel, T. Green, and J. Stoustrup, “Improving demand response potential of a supermarket refrigeration system: A food temperature estimation approach,” _IEEE Transactions on Control Systems Technology_ , vol. 25, no. 3, pp. 855–863, 2017.
  * [6] H. Hao, B. M. Sanandaji, K. Poolla, and T. L. Vincent, “Aggregate flexibility of thermostatically controlled loads,” _IEEE Transactions on Power Systems_ , vol. 30, no. 1, pp. 189–198, 2014.
  * [7] M. K. Petersen, K. Edlund, L. H. Hansen, J. Bendtsen, and J. Stoustrup, “A taxonomy for modeling flexibility and a computationally efficient algorithm for dispatch in smart grids,” in _2013 American Control Conference_ , Washington, DC, USA, 2013, pp. 1150–1156.
  * [8] G. De Zotti, S. A. Pourmousavi, J. M. Morales, H. Madsen, and N. Kjølstad Poulsen, “Consumers’ flexibility estimation at the TSO level for balancing services,” _IEEE Transactions on Power Systems_ , vol. 34, no. 3, pp. 1918–1930, 2019.
  * [9] J. R. Schaperow, S. Gabriel, M. Siemann, and J. Crawford, “A simulation-based model for optimal demand response load shifting: A case study for the Texas power market,” _Journal of Energy Markets_ , vol. 12, no. 4, pp. 53–80, 2019\.
  * [10] P. Chanpiwat, S. A. Gabriel, R. L. Moglen, and M. J. Siemann, “Using cluster analysis and dynamic programming for demand response applied to electricity load in residential homes,” _ASME Journal of Engineering for Sustainable Buildings and Cities_ , vol. 1, no. 1, p. 011006, 2020.
  * [11] R. L. Moglen, P. Chanpiwat, S. A. Gabriel, and A. Blohm, “Optimal thermostatically-controlled residential demand response for retail electric providers,” _Energy Systems_ , vol. 14, pp. 641–661, 2023.
  * [12] M. Siemann, “Performance and applications of residential building energy grey-box models,” 2013, PhD Dissertation, University of Maryland, College Park.
  * [13] I. I. Gihman and A. V. Skorohod, _The theory of stochastic processes II_. Springer, 2004.
  * [14] L. Bobo, S. Delikaraoglou, N. Vespermann, J. Kazempour, and P. Pinson, “Offering strategy of a flexibility aggregator in a balancing market using asymmetric block offers,” in _Power Systems Computation Conference (PSCC)_ , Dublin, Ireland, 2018, pp. 1–7.
  * [15] Ecogrid, “Ecogrid 2.0.” [Online]. Available: http://ecogrid.dk/
  * [16] B. Biegel, P. Andersen, J. Stoustrup, L. H. Hansen, and D. V. Tackie, “Information modeling for direct control of distributed energy resources,” in _American Control Conference_ , Washington, DC, USA, 2013, pp. 3498–3504.
  * [17] S. Harbo and B. Biegel, “Contracting flexibility services,” in _IEEE PES ISGT Europe_ , Lyngby, Denmark, 2013, pp. 1–5.
  * [18] B. Biegel, P. Andersen, T. S. Pedersen, K. M. Nielsen, J. Stoustrup, and L. H. Hansen, “Electricity market optimization of heat pump portfolio,” in _IEEE International Conference on Control Applications (CCA)_ , Hyderabad, India, 2013, pp. 294–301.
  * [19] B. Biegel, P. Andersen, J. Stoustrup, M. B. Madsen, and L. H. Hansen, “Lumped thermal household model,” in _IEEE PES ISGT Europe 2013_ , Lyngby, Denmark, 2013, pp. 1–5.
  * [20] R. Juhl, J. K. Møller, and H. Madsen, “ctsmr–continuous time stochastic modeling in R,” 2016. [Online]. Available: https://arxiv.org/abs/1606.00242
  * [21] H. Madsen, _Time series analysis_. Chapman and Hall/CRC, 2007.
  * [22] “Introduktion til systemydelser,” _Energinet_ , 2021, accessed: 2021-09-06. [Online]. Available: https://energinet.dk/-/media/8F2B403B3A724197A179FF6DFD9A926C.pdf
  * [23] G. P. McCormick, “Computability of global solutions to factorable nonconvex programs: Part I–Convex underestimating problems,” _Mathematical Programming_ , vol. 10, pp. 147–175, 1976.
  * [24] 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_ , vol. 3, no. 1, pp. 1–122, 2011.
  * [25] M. Hong, Z.-Q. Luo, and M. Razaviyayn, “Convergence analysis of alternating direction method of multipliers for a family of nonconvex problems,” _SIAM Journal on Optimization_ , vol. 26, no. 1, pp. 337–364, 2016.
  * [26] “Github repository,” 2023. [Online]. Available: https://github.com/PeterAVG/paper_flex
  * [27] “MARI,” _ENTSO-e_ , 2021. [Online]. Available: https://www.entsoe.eu/network_codes/eb/mari/
  * [28] J. M. Morales, A. J. Conejo, H. Madsen, P. Pinson, and M. Zugno, _Integrating renewables in electricity markets: Operational problems_. Springer Science & Business Media, 2013.