# Reduced order modelling using parameterized non-uniform boundary conditions
in room acoustic simulations

Hermes Sampedro Llopis<EMAIL_ADDRESS>Acoustic Technology Group,
Department of Electrical and Photonics Engineering, Technical University of
Denmark, Kongens Lyngby, Denmark  Rambøll Denmark, Copenhagen, Denmark  Cheol-
Ho Jeong  Acoustic Technology Group, Department of Electrical and Photonics
Engineering, Technical University of Denmark, Kongens Lyngby, Denmark  Allan
P. Engsig-Karup  Scientific Computing Section, Department of Applied
Mathematics and Computer Science, Technical University of Denmark, Kongens
Lyngby, Denmark

###### Abstract

Quick simulations for iterative evaluations of multi-design variables and
boundary conditions are essential to find the optimal acoustic conditions in
building design. We propose to use the reduced basis method (RBM) for
realistic room acoustic scenarios where the surfaces have inhomogeneous
acoustic properties, which enables quick evaluations of changing absorption
materials for different surfaces in room acoustic simulations. The RBM has
shown its benefit to speed up room acoustic simulations by three orders of
magnitude for uniform boundary conditions. This study investigates the RBM
with two main focuses, 1) various source positions in diverse geometries,
e.g., square, rectangular, L-shaped, and disproportionate room. 2)
Inhomogeneous surface absorption in 2D and 3D by parameterizing numerous
acoustic parameters of surfaces, e.g., the thickness of a porous material,
cavity depth, switching between a frequency independent (e.g., hard surface)
and frequency dependent boundary condition. Results of numerical experiments
show speedups of more than two orders of magnitude compared to a high fidelity
numerical solver in a 3D case where reverberation time varies within one just
noticeable difference in all the frequency octave bands.

##  I  Introduction

Room acoustic simulations are typically used during the design stage of a
building to find the optimal amount, choice and position of materials
according to the use of the room. Poor acoustic conditions can produce a
negative effect, e.g., a decrease in the work productivity [1, 2], a decrease
in the learning quality [3, 4], and an increase in the stress level [5, 6].
Performing room acoustic simulations by solving the acoustic wave equation
using numerical methods, such as Finite Difference Time-Domain methods
(FDTD)[7], Finite Element Methods (FEM) [8] or Spectral Element Methods
(SEM)[9] is accurate as all important wave phenomena can be accounted for.
However, it is computationally expensive compared to geometric acoustics
methods [10]. As a result, active research seeks to find ways to speed up the
numerical methods to be used in building design. Common strategies include
using model order reduction (MOR) [11], parallel computing on multi-core
processing units [12], and more recently also machine learning techniques
[13].

A review of the state-of-the-art in MOR across disciplines can be found in the
literature [14]. The reduced basis method (RBM) [11] is a method belonging to
the class of MOR techniques. It exploits the parametric dependence in the
solution of a partial differential equation (PDE) by combining different
solutions given by the variation of a set of parameter values. Problems of
large systems of PDEs are effectively reduced to a low dimensional subspace to
achieve computational acceleration and transformed back to the original
problem size [15, 16, 17, 18, 19, 20, 21, 22], however, the speedup is
dependent on the problem of interest. MOR have been recently presented in some
acoustic applications [23, 24, 25] as well as in many different fields, e.g.,
electromagnetics [26, 27], computational fluid dynamics [28, 29], heat
transfer [30], vibroacoustics [31, 32, 33] and many others [34, 35, 36, 37,
38], demonstrating, in general, an efficient reduction in the computational
burden. Using RBM in room acoustic simulations with parameterized boundary
conditions during the design stage of an indoor space may allow exploring many
acoustic conditions by solving the reduced system and exploiting the benefit
of a reduced computational cost under the variation of several parameters,
e.g., the thickness of a porous material, air gap distance, etc. Today, the
application of MOR to room acoustic simulations with boundary parameterization
is scarce despite the significant potential for acceleration. A recent study
[39] presents a model order reduction strategy using a Krylov subspace
algorithm in the time domain with a FEM solver where speedups of 11–36 were
demonstrated. The study is performed for a simple domain where all the
surfaces are assumed to be rigid boundaries, and only the floor is modelled
with a surface impedance. Moreover, another study demonstrates the potential
of RBM in room acoustic applications, however, only for simplified homogeneous
boundary conditions [40]. Two to three orders of speedup factors were achieved
in the online stage.

The study of RBM with realistic inhomogeneous boundary conditions, where the
absorption materials are distributed inhomogeneously among the surfaces, is
scarce in the literature. For example, classrooms typically have highly
absorbing ceilings and scattering objects near the rear surface and window on
one side wall. In this study, we investigate how to effectively apply RBM in a
realistic setting with inhomogeneous boundary conditions, how the RBM
performance varies with the room geometry and source/receiver location, and
what acceleration is expected in such conditions.

The novelty of this investigation is to construct and analyze the performance
of a RBM strategy for realistic scenarios with two focuses, one being the RBM
performance in various room shapes and source locations and two being the
performance when including numerous parameters of acoustic materials
distributed inhomogeneously. First, this study presents a conceptual proof-of-
concept for a complex 2D case, as it simplifies the computational burden.
Second, two 3D rooms are analyzed. This study is essential for future work,
e.g., scaling the method to extend the RBM for large building projects.

In Sec. II the governing equations, the boundary conditions, the reduced basis
method and the error measures are described together with the different
domains and simulation parameters. Section III handles the numerical
experiments and results, which are analyzed and discussed in Sec. IV.

##  II Methods

This section presents an overview of the methods and simulation conditions
used in this study. A more detailed description of the full order model (FOM)
solver used as a reference model and the ROM is deeply described in previous
work [40].

###  II.1 Governing equations and boundary conditions

We consider the acoustic wave equation in the Laplace domain

$s^{2}p-c^{2}\Delta p=0,$ (1)

where $p(\bm{x},t)$ is the sound pressure, $\bm{x}\in\Omega$ the position in
the domain $\Omega\subset\mathbb{R}^{d}$ with $d=\\{2,3\\}$, $t$ is the time
in the interval $(0,T]$ s and $c$ is the speed of sound ($c=343$ m/s).
Equation (1) is discretized using the SEM formulation, which is well known and
an overview can be found [41, 9, 42]. The final formulation written in the
Laplace domain is given by

$\displaystyle\big{(}s^{2}\mathcal{M}+c^{2}\mathcal{S}+sc^{2}\frac{\rho}{Z_{s}}\mathcal{M}_{\Gamma}\big{)}\mathbf{p}=0,$
(2)

where $\mathcal{M}$ refers to the mass matrix, $\mathcal{S}$ is the stiffness
matrix, $\rho$ is the density of the medium ($\rho=1.2$ kg/m3) and $Z_{s}$ is
the surface impedance. Note that the impedance boundaries are considered
denoting the boundary domain as $\Gamma$.

The frequency dependent boundary conditions are implemented via the method of
auxiliary differential equations (ADE) [43, 44, 9]. The surface impedance of a
porous absorber is modelled using Miki’s model [45] in conjunction with a
transfer matrix method [46], and mapped to a six pole rational function by
using a vector fitting algorithm [47] so that the surface admittance
$Y_{s}=1/Z_{s}$ can be written as a rational function and expressed using
partial fraction decomposition [44]. Then, the system (2) can be stated in the
form of a linear system of equations and solved using a sparse direct solver
as presented in [9, 40]

$\displaystyle\mathbf{K}\mathbf{p}=0,\quad\mathbf{K}\in\mathbb{R}^{N\times
N},$ (3)

where, $\mathbf{K}$ refers to the operators shown in (2) and $N$ corresponds
to the degree of freedom (DOF). Note that if the system is split into real and
imaginary parts for implementation purposes, the size of the operator is
$\mathbf{K}\in\mathbb{R}^{2N\times 2N}$ [48, 40]. The system (2) is
initialized using a Gaussian pulse with a spatial distribution $\sigma_{g}$
that determines the frequencies to span, by adding the right hand side term
$s\mathcal{M}p_{0}$, where $p_{0}$ is the initial sound pressure state in the
time domain.

The solution in the Laplace domain is finally transformed to the time domain
by means of the Weeks method [49, 40].

###  II.2 The reduced order model

The purpose of using RBM is to substantially reduce the size of the problem
while ensuring a certain level of accuracy. Specifically, the DOF are reduced,
and the techniques succeed when RDOF $\ll$ DOF, RDOF being the corresponding
degrees of freedom in the numerical scheme after applying RBM. The RBM
consists of two stages. First, one or more variables present in the partial
differential equation or its discretized form (2) are chosen as parameters,
e.g., $Z_{s}$, spanning a discrete range of values. Then, in the first stage,
referred to as the offline stage, the parameter space is explored to generate
a problem-dependent basis by collecting FOM solutions for different parameter
values within the range of interest. A Galerkin projection takes place to
reduce the dimensionality of the problem by utilizing the generated basis. In
the second stage, referred to as the online stage, the reduced problem is
solved for a new parameter value that was not explored in the offline stage at
a much lower computational cost. The offline stage is typically
computationally costly as it requires multiple FOM solutions that capture
relevant information on the parameters variations and allow generating of
representative basis functions of that variation. The reduction of the size of
the computational problem comes with the truncation of the basis. The solver
is stated in the Laplace domain to ensure the stability of the ROM solution
[48]. The basis generation can be seen as a data-driven technique based on
proper orthogonal decomposition (POD). It relies on a proper symplectic
decomposition (PSD) with a symplectic Galerkin projection. Specifically, the
cotangent-lift method introduced in [50] is applied, which preserves the
structure of the operators when the problem is split and solved into real and
imaginary parts. The ROM solution is expressed as an expansion of the basis
functions $\phi_{i}$ and coefficients $a_{i}$, which is represented as

$\displaystyle p_{rom}=\mathbf{\Phi}\mathbf{a},$ (4)

where $\Phi_{ij}\equiv\phi_{i}(x_{j})$. Inserting (4) into (3) yields to a
similar problem, where now the system is solved for the coefficients
$\mathbf{a}$

$\displaystyle\mathbf{K}_{rom}\mathbf{a}=0,$ (5)

where $\mathbf{K}_{rom}=\bm{\Phi}_{cl}^{T}\mathbf{K}\bm{\Phi}_{cl}$ and
$\bm{\Phi}_{cl}$ defines the symplectic basis constructed as

$\bm{\Phi}_{cl}=\begin{bmatrix}\bm{\Phi}&\bm{0}\\\
\bm{0}&\bm{\Phi}\end{bmatrix}.$ (6)

Moreover, the reduced operator can be written also as

$\displaystyle\mathbf{K}_{rom}=s^{2}\bm{\Phi}^{T}\mathcal{M}\bm{\Phi}+c^{2}\bm{\Phi}^{T}\mathcal{S}\bm{\Phi}+sc^{2}\frac{\rho}{Z_{s}}\bm{\Phi}^{T}\mathcal{M}_{\Gamma}\bm{\Phi},$
(7)

where a new parameter value can be chosen during the online stage (for
example, $Z_{s}$). The generation of a basis comes by first collecting all the
FOM solutions obtained during the offline stage in the snapshot matrix[40]
$\bm{S}_{N}\in\mathbb{R}^{N\times 2N_{s}}$, where $N$ is already described
above, and $N_{s}$ is the number of evaluated complex frequencies determined
by the Weeks method [48]. $\bm{S}_{N}$ can be decomposed using the proper
orthogonal decomposition (POD) technique based on a singular value
decomposition (SVD) to get the corresponding basis functions, defined as
$\bm{\Phi}=[U_{1},...,U_{N_{rb}}]\in\mathbb{C}^{N\times N_{rb}}$. The reduced
basis is chosen through truncation of this basis relying on the rate of decay
of the singular values.

The singular value decay shows the energy distribution among the basis,
providing information about the reduction of the problem. It is defined as

$E/E_{0}=\frac{\textrm{diag}(\sigma_{1},..,\sigma_{N})}{\sum_{i=1}^{N}\sigma_{i}},$
(8)

where $\Sigma=\textrm{diag}(\sigma,..,\sigma_{N})$ is obtained by SVD where
$S=U\Sigma V^{T}$. The number of basis can be chosen so that the projection
error is smaller than a given tolerance $\epsilon_{POD}$

$I(N_{rb})=\frac{\sum_{i=1}^{N_{rb}}\sigma^{2}_{i}}{\sum_{i=1}^{N}\sigma^{2}_{i}}\geq
1-\epsilon_{POD},$ (9)

where $N_{rb}$ denotes the number of basis functions, $I(N_{rb})$ represents
the percentage of the energy of the collection of FOM solutions captured by
the first. For a given $\epsilon_{POD}$, the faster the energy decays, the
smaller number of basis needed, and thus, a better reduction of the problem is
expected. The reduction comes with a truncation of the basis, which determine
the size of (7) as $N_{rb}$.

###  II.3 Error measures

In this study, two type of errors are considered to compare the ROM against
the FOM. First, the relative error using the root mean square (rms) pressure
is introduced

$\epsilon_{rel}=\frac{p_{rms_{ROM}}-p_{rms_{FOM}}}{p_{rms_{FOM}}}\times
100\quad(\%).$ (10)

Second, the error in the frequency domain expressed in dBs is considered

$\Delta L(f)=20\log_{10}\Big{|}\frac{p_{FOM}(f)}{p_{ROM}(f)}\Big{|},$ (11)

where, $p_{FOM}$ and $p_{ROM}$ are the sound pressure of the FOM and ROM
respectively along the frequency spectrum.

The performance of the ROM is measured in terms of speedups ($sp$) defined as

$sp=\frac{CPU_{FOM}}{CPU_{ROM}},$ (12)

where $CPU_{FOM}$ and $CPU_{ROM}$ corresponds to the computational time of the
FOM and ROM respectively.

###  II.4 Test rooms and simulation conditions

First, Section III.1, deals with ROMs with several 2D geometries with
different source locations illustrated in Figure 1. Four different geometries
and two source positions are considered, one at the corner and one at the
centre. First a $4\times 4$ m2 square domain with the source placed at
$(s_{x_{1}},s_{y_{1}})_{SQ}=(0.2,0.2)$ m (SQ1), and
$(s_{x_{2}},s_{y_{2}})_{SQ}=(2,2)$ m (SQ2) is introduced (Figure 1a). Second,
a $4\times 2.5$ m2 rectangular domain where
$(s_{x_{1}},s_{y_{1}})_{RC}=(0.2,0.2)$ m (RC1), and
$(s_{x_{2}},s_{y_{2}})_{SQ}=(2,1.25)$ m (RC2) is considered (Figure 1b).
Third, an L-shaped room where the long side is $4$ m and the short side is $2$
m is considered with only one source position at the corner
$(s_{x},s_{y})_{LS}=(0.2,0.2)$ m (LS1) (Figure 1c). Finally a corridor shape
of size $10\times 1$ m2 is presented where $(s_{x},s_{y})_{CO}=(0.2,0.2)$ m
(CO1) (Figure 1d). The maximum element size is selected considering triangular
high-order elements ($P=4$) and using $4$ points per wavelength ($PPW$)
leading into an upper frequency $f_{u}=2.8$ kHz, which is approximately the
upper cutoff frequency of the $2$ kHz octave band. The model is excited with a
Gaussian pulse as initial condition with $\sigma_{g}=0.1$ m2 [51]. The ROMs
for each room type and source position are built by sampling the surface
impedance of all the sides at the following values $Z_{s}=[500,5250,10000]$
Nsm-3.

Second, section III.2 deals with several 2D ROMs for an inhomogeneous
distribution of the acoustic material, which in this study is named
inhomogeneous boundary conditions. A 2D rectangular room ($4$ m $\times 2.5$
m) with inhomogeneous boundary conditions is considered with the same
simulation parameters as before. The source is placed at
$(s_{x},s_{y})=(3,1.2)$ m and the receiver is at $(r_{x},r_{y})=(1,1.2)$ m.
Moreover, an additional ROM is constructed for this section with an upper
frequency of $f_{u}=4$ kHz.

Third, for section III.3, two 3D models are considered. A $1$ m cube (CB)
enclosure is presented, where three of the six surfaces are parameterized.
Simulations are carried out using a polynomial order of $P=4$ with $N=35937$
elements. Assuming $PPW=4$, the upper frequency is given by $f_{u}=2.8$ kHz.
The model is excited with a Gaussian pulse as an initial condition with
$\sigma_{g}=0.1$ m2 placed at ($s_{x},s_{y},s_{z})_{CB}=(0.7,0.5,0.5)$ m. The
receiver position is placed at ($r_{x},r_{y},r_{z})=(0.25,0.25,0.50)$ m. A
second 3D room is considered to follow a good ratio (GR) of 1.9:1.4:1, which
assures an even distribution of the room modes[52]. The room size is
$(L_{x},L_{y},L_{z})=(1.615,1.190,0.850)$ m, sound source is placed at
($s_{x},s_{y},s_{z})_{GR}=(1.200,0.600,0.425)$ m and the receiver
($r_{x},r_{y},r_{z})=(0.500,0.200,0.425)$ m. Again, a polynomial order of
$P=4$ with $N=35937$ is considered. Assuming a spatial resolution
corresponding to about 4 $PPW$ for the highest frequencies ($f_{u}=1.7$ kHz).

fig1.eps8cm

Figure 1: The geometries and source positions.

###  II.5 Inhomogeneous boundary parameterization

First, the 2D domain is considered in section III.2. The ceiling (CE) is
modelled with a porous absorber. The key parameters affecting the absorption
characteristic of a porous absorber are the flow resistivity, thickness, and
air cavity depth [53]. To parameterize the porous ceiling, FOMs with
$\sigma_{mat}=[10,30,50]$ kNsm-4, $d_{mat}=[0.02,0.12,0.22]$ m, and
$d_{0}=[0.02,0.12,0.22]$ m are simulated in the offline stage. The floor (FL)
is designed with two different options: as a hard surface modelled as a
frequency independent boundary and covered with a carpet modelled as a porous
layer. Frequency independent boundary conditions are ranges $Z_{s}=[10,50,90]$
kNsm-3, while frequency dependent boundary conditions to model the carpet are
computed with a fixed thickness of $0.02$ m but varying
$\sigma_{mat}=[10,30,50]$ kNsm-4. Finally, the left wall (WL) and right (WR)
walls are modelled as porous panels where $d_{mat}=0.03$ m, $d_{0}=0$ m and
$\sigma_{mat}=[5,12,19]$ kNsm-4. Figure 2 shows the absorption coefficients
for the most and least absorptive cases of each surface, whose corresponding
parameter values are given in Table 1.

fig2a.eps0.5(a) fig2b.eps0.5(b) fig2c.eps0.5(c)

Figure 2: Absorption coefficients. a) Ceiling (CE), b) Floor (FL), c) Walls
(W). Table 1: Parameter values of the presented absorption coefficients of
Figure 2.

|  $\sigma_{mat}$ [kNsm-4] |  $d_{mat}$ [$m$] |  $d_{0}$ [$m$] |  $Z_{s}$ [kNsm-3]
---|---|---|---|---
CEA | 10 | 0.02 | 0.02 | -
CEB | 30 | 0.12 | 0.12 | -
CEC | 10 | 0.12 | 0.22 | -
CED | 30 | 0.02 | 0.22 | -
FLA | 10 | 0.02 | 0 | -
FLB | 30 | 0.02 | 0 | -
FLC | 50 | 0.02 | 0 | -
FLD | - | - | - | 10
FLE | - | - | - | 90
WA | 5 | 0.03 | 0 | -
WB | 12 | 0.03 | 0 | -
WC | 19 | 0.03 | 0 | -

A way to construct the ROM is by performing FOM simulations for each
combination of the parameter values. To cover the inhomogeneous boundary
variation, a total of $3^{7}$ $(2187)$ FOM simulations are possible. This is
way too many for practical runtime constraints. Instead, we have chosen only
$3\times 7$ $(=21)$ FOM simulations, which is shown in the Appendix (Table A1)
indicating for each FOM simulation which parameter values are chosen and which
are fixed. The table rows correspond to the 21 simulations, and the columns
correspond to the different parameter values shown in Table 2. The chosen
parameters are marked with an X. To the author’s knowledge, this study is the
first to report MOR performances with such complicated room acoustic setups
and is, therefore, useful to establish the feasibility of the approach.

Table 2: Parameter values chosen to construct the ROM in 2D. The parameters
under variation that are included in the ROM are marked with *.

| CE | FL | WL | WR
---|---|---|---|---
$\sigma_{mat}$ ${[}kNsm^{-4}$] | [10, 30, 50]* | [10, 30, 50]* | [5, 12, 19]* | [5, 12, 19]*
$d_{mat}$ [$m$] | [0.02, 0.12, 0.22]* | 0.02 | 0.03 | 0.03
$d_{0}$ [$m$] | [0.02, 0.12, 0.22]* | 0 | 0 | 0
$Z_{s}$ [kNsm-3] | - | [10, 50, 90]* | - | -

Second, section III.3 considers the 3D domain for both cube and good ratio
shapes. In both cases, the ceiling (CE) is modelled as a porous acoustic
material of a fixed thickness $d_{mat}=0.05$ m, where
$\sigma_{mat}=[10,30,50]$ kNsm-4 and $d_{0}=[0.02,0.12,0.22]$ m are
parameterized. The floor (FL) is modelled as a frequency independent boundary
with $Z_{s}=50$ kNsm-3. The east wall (WE) is modelled as a frequency
independent brick surface with $Z_{s}=50$ kNsm-3. The south wall (WS) is not
parameterized, and it is covered with a porous material where $\sigma_{mat}=7$
kNsm-4, $d_{mat}=0.02$ m and $d_{0}=0$ m. The west wall (WW) is modelled with
a porous material of $d_{mat}=0.05$ m and $d_{0}=0$ m whose flow resistivity
is parameterized with values of $\sigma_{mat}=[5,12,19]$ kNsm-4. The north
wall (WN) is modelled as a hard surface with $Z_{s}=50$ kNsm-3. In addition, a
$0.5\times 0.5$ m2 square acoustic panel made of porous material is placed at
the centre of the wall. The panel has a fixed thickness and flow resistivity
of $d_{mat}=0.1$ m and $\sigma_{mat}=30$ kNsm-4 respectively. The air gap
between the panel and the wall is parameterized $d_{0}=[0.2,0.12,0.22]$ m.
Table 3 summarizes the parameter values of each surface.

The ROM is constructed for each 3D room constructed with four different
parameters (2 for CE, 1 for WW and 1 for WN) with three different values each.
Thus, a total number of $3^{4}=81$ FOM simulations are possible. Instead,
$3\times 4=12$ FOM simulations were carried out to construct the ROM in the
same way described for the 2D case. The chosen parameter values for each FOM
simulation are shown in the Appendix (Table A2).

Table 3: Parameter values chosen to construct the ROM in 3D. The parameters
under variation that are included in the ROM are marked with *.

| CE | FL | WE | WS | WW | WN
---|---|---|---|---|---|---
$\sigma_{mat}$ [$kNsm^{-4}$] | [10, 30, 50]* | - | - | 70 | [5, 12, 19]* | 30
$d_{mat}$ [$m$] | 0.05 | - | - | 0.02 | 0.05 | 0.1
$d_{0}$ [$m$] | [0.02, 0.12, 0.22]* | - | - | 0 | 0 | [0.02, 0.12, 0.22]*
$Z_{s}$ [$kNsm^{-3}$] | - | 50 | 50 | - | - | -

##  III  Results

This section presents results mostly with the singular energy decay, E/Eo,
relative error shown in (10), speedups (12), sound pressure level (SPL)
spectrum, and the reverberation time (RT).

###  III.1 2D - Influence of the source position and geometry

For the geometries and source positions tested, the singular value decay is
shown in Figure 3. A faster decay means that a larger portion of the energy is
concentrated in the first singular values, effectively indicating that a
smaller number of basis $N_{rb}$ is needed for a given error tolerance. Slower
decays would need a larger number of $N_{rb}$ to provide the same error. Note
that the smaller $N_{rb}$ is, the higher speedups are achieved. In Figure 3,
one can see that the more symmetric the problem is, both in terms of geometry
and source location, the faster the decay is, as measured in the singular
values translating into more efficiency (speedup). For example, $SQ$ shows a
faster decay than $RC$ and $LS$. However, the difference is not significant
until the energy reaches the value of $10^{-10}$, and it can be concluded that
the room geometry does not significantly change the energy distribution among
the basis. On the other hand, the centred sound source locations lead to a
faster decay of the singular values compared to the corner source. $SQ2$ shows
a faster decay than $SQ1$, and the same for $RC2$ in comparison to $RC1$. This
is because placing a source at the centre of the room will fail to excite some
room modes, of which the nodal lines/points coincide with the source location.
This leads to a smaller number of basis needed to describe the physical
dynamics of the wave propagation in the room accurately.

fig3.eps7cm

Figure 3: Singular value decay for the first 2500 modes of the basis and
energy distribution $(E/E_{0})$ among the basis with frequency independent
boundaries.

###  III.2 2D - Parameterization of different absorption properties

In Figure 4, the singular value decay when a different number of parameters
are included in the model, i.e., only the ceiling parameters are included in
the ROM (CE); the ceiling and the floor parameters (CE+FL); the ceiling, the
floor and the left wall (CE+FL+WL) and all the sides (CE+FL+WL+WR). The more
parameters added, the slower the decay curve, so adding more parameters to the
ROM has a clear effect on the singular value decay and, thus, the choice of
the basis for ROM. However, the decay remains similar in the first basis
functions, which are the ones used to construct the ROM.

fig4.eps7cm

Figure 4: Singular value decay for different number of parameters for case 1
and case 2. Table 4: Parameter values for the online stage of the 2D ROM.
Values marked with * denotes the parameters which are parameterized.

|  $\sigma_{mat}$ [$kNsm^{-4}$] |  $d_{mat}$ [$m$] |  $d_{0}$ [$m$] |  $Z_{s}$ [$kNsm^{-3}$]
---|---|---|---|---
CE1 | 2* | 0.1* | 0.1* | -
FL1 | 12* | 0.02 | 0 | -
WL1 | 10* | 0.03* | 0 | -
WR1 | 15* | 0.03 | 0 | -
CE2 | 45* | 0.05* | 0.2* | -
FL2 | - | - | - | 7*
WL2 | 10* | 0.2* | 0 | -
WR2 | 6* | - | 0 | -

In the online stage, two ROM cases are simulated and compared with their
corresponding FOM. Case 1 deals with frequency dependent boundaries, while
case 2 allows to change the floor to a rigid surface, modelled as a frequency
independent boundary. The boundary parameter values are shown in Table 4, and
the IR and SPL are shown in Figure 5. Several different online ROM simulations
are compared for $N_{rb}=69,158,274,410,571,752$ for the ROM with an upper
frequency of 2 kHz corresponding to values of
$\epsilon_{POD}=10^{-2},10^{-3},...,10^{-7}$; and $N_{rb}=132,310,531,810$ for
the ROM with an upper frequency of 4 kHz corresponding to values
$\epsilon_{POD}=10^{-2},...,10^{-5}$. Figure 5 shows the results up to 4 kHz
with $N_{rb}=531$. Figure 5 shows the impulse response, and Figure 5 presents
the sound pressure level (SPL) from 20 Hz to 4 kHz showing that $\Delta L$ is
nearly zero below $600$ Hz and increasing with the frequency. Note that around
$4$ kHz, a roll-off is presented due to the source excitation. This case
enables a computational speedup by a factor of 37 for an error of
$\epsilon_{rel}=0.03\%$. For the ROM constructed with $N_{rb}=274$ and an
upper frequency of 2 kHz, the error is $\epsilon_{rel}=0.2\%$ for case 1 with
a speedup of 70. Moreover, case 2 presents an error of $\epsilon_{rel}=0.3\%$
and a speedup of 50. Note that the accuracy depends on $N_{rb}$. The speedup
against $\epsilon_{rel}$ given in (10) is shown in Figure 6. These results
show speedups around two orders of magnitude for 2D.

fig5a.eps.4(a) fig5b.eps.4(b)

Figure 5: (a) Impulse response and (b) spectrum of 2D FOM and ROM of case 1
for $N_{rb}=531$ and $f_{u}=4$ kHz (using the parameter in Table 4).

fig6.eps7cm

Figure 6: Speedup against the relative rms error for the 2D domain (case 1)
with $f_{u}=2$ kHz, and 3D CB.

###  III.3 3D - Parameterization of different absorption properties

This section presents a similar analysis with 3D test cases. By varying
$N_{rb}$, simulations are compared. The chosen parameters for the online
simulation are shown in Table 5. For the CB the following number of basis are
considered $N_{rb}=30,81,175,275,543,837$ corresponding to values of
$\epsilon_{POD}=10^{-2},...,10^{-9}$ to be compared against the corresponding
FOM solution for verification purposes. Moreover, for the GR room, the ROM was
constructed with $N_{rb}=62,303,478,649,1500$ corresponding to values of
$\epsilon_{POD}=10^{-2},...,10^{-9}$. Table 5 shows the chosen parameter
values for all the surfaces for both the cubic and good ratio domains. Figure
7 and Figure 7 present the impulse response and frequency response,
respectively, for the CB. Moreover, Figure 7 and Figure 7 show the impulse
response and frequency response, respectively, for the GR room. In both cases,
the $\Delta L$ is included, which confirms a good agreement between FOM and
ROM for the given $N_{rb}$. The CB is constructed with $N_{rb}=273$ where the
error is $\epsilon_{rel}=0.04\%$ with a speedup of 143. On the other hand, the
GR room is constructed with $N_{rb}=478$, where the error is
$\epsilon_{rel}=0.66\%$ with a speedup of 90. Note $\Delta L$ in Eq. (11) and
$\epsilon_{rel}$ are also presented for the different number of basis $N_{rb}$
in Figure 8 for the CB (Figure 8) and the good ratio room (Figure 8). Again,
the $\Delta L$ is nearly zero at lower frequencies and increases with
frequency, showing more differences at the anti-resonances. The speedup
against the error for CB case 1 is presented in Figure 6, showing values
around three orders of magnitude.

In order to understand the behaviour of different room ratios, the number of
basis functions $N_{rb}$ for a given tolerance $\epsilon_{POD}$ described in
(9) is compared. A new ROM of the cubic room is computed with $f_{u}=1.7$ kHz.
Figure 9 shows the comparison for GR1.7kHz, CB1.7kHz and CB2.8kHz. Results
show that for a fixed upper frequency, a less symmetric geometry GR1.7kHz
needs more basis functions for a given $\epsilon_{POD}$ compared to a
symmetric geometry CB1.7kHz. This would not necessarily lead to lower speedups
considering that GR has a larger number of DOF. Moreover, comparing the curves
corresponding to GR1.7kHz and CB2.8kHz for a fixed number of DOF, GR1.7kHz
results in a larger number of $N_{rb}$ for a given $\epsilon_{POD}$. Thus, it
will lead to lower speedups as the numerator of equation (12) remains the same
for GR1.7kHz and CB2.8kHz. At the same time, the denominator becomes larger
for GR1.7kHz at a given $\epsilon_{POD}$ compared to CB2.8kHz.

The reverberation time is one of the most widely used acoustic parameters
defined in ISO 3382-1 [54] as the time needed for the energy to decrease by 60
dB. A way to evaluate the error of the reverberation time between the FOM and
ROM is by means of the JND, which is the minimum change in the RT that can be
perceptually perceived. The JND of the RT is 5% as defined in the standard
[54]. Note that if the difference is larger than 1 JND, the IR can be
potentially differently heard. $T_{20}$ is calculated from IRs via FOM and ROM
to quantify how ROM degrades the accuracy of RT. Figure 10 shows the RT for
different frequency octave bands. The CB ROM is performed with $N_{rb}=185$
while the GR ROM with $N_{rb}=649$. The RT difference is below one JND in all
the frequency octave bands in both cases, which indicates that the present ROM
would not be perceptually different compared to the FOM. Moreover, Figure 10
and Figure 10 show the numbers of JND for $T_{20}$ for various $N_{rb}$.
Decreasing $N_{rb}$ increases the RT difference in the higher frequency bands.
Note that in this case the ROM is 365 times faster than the FOM and the CB
needs fewer basis functions than GR, which supports the finding that symmetric
conditions are more favourable for higher reductions.

Table 5: Online stage boundary parameters for the 3D rooms. Values marked with
* denotes parameterization.

CB domain
---
| CE | FL | WE | WS | WW | WN
$\sigma_{mat}$ [$kNsm^{-4}$] | 12* | - | - | 7 | 10* | 30
$d_{mat}$ [m] | 0.05 | - | - | 0.02 | 0.05 | 0.1
$d_{0}$ [m] | 0.06* | - | - | 0 | 0 | 0.1*
$Z_{s}$ [$kNsm^{-3}$] | - | 50 | 50 | - | - | 50
GR domain
| CE | FL | WE | WS | WW | WN
$\sigma_{mat}$ [$kNsm^{-4}$] | 10.5* | - | - | 7 | 5.5* | 30
$d_{mat}$ [m] | 0.05 | - | - | 0.02 | 0.05 | 0.1
$d_{0}$ [m] | 0.025* | - | - | 0 | 0 | 0.03*
$Z_{s}$ [$kNsm^{-3}$] | - | 50 | 50 | - | - | 50

fig7a.eps.5(a) fig7b.eps.5(b) fig7c.eps.5(c) fig7d.eps.5(d)

Figure 7: Simulated pressure using the 3D FOM and ROM. a) CB domain sound
pressure with $N_{rb}=275$, b) CB domain frequency response with $N_{rb}=275$,
c) GR domain sound pressure with $N_{rb}=478$, d) GR domain frequency response
with $N_{rb}=478$.

fig8a.eps.5(a) fig8b.eps.5(b)

Figure 8: $\Delta L$ error. a) CB case 1 for $N_{rb}=81,275,543,837$. b) GR
for $N_{rb}=303,649,1500$.

fig9.eps7cm

Figure 9: Number of basis functions $N_{rb}$ against the tolerance
$\epsilon_{POD}$ introduced in (9) for GR and CB at two different upper
frequencies $1.7$ kHz and $2.8$ kHz

fig10a.eps9cm(a) fig10b.eps9cm(b) fig10c.eps9cm(c)

Figure 10: Comparison of reverberation time between FOM and ROM and their
difference in terms of numbers of JND. a) Cube domain with $N_{rb}=185$ (CB)
and good ratio domain with $N_{rb}=649$ (GR) including the number of JNDs per
octave band, b) Number of JNDs for CB, c) Number of JNDs for GR.

##  IV Discussion

The performance behaviour of the ROM in terms of speedups is case-dependent
and can be challenging to predict, especially when including a large number of
parameters in a complex scene. This study analyzes the behaviour of realistic
room acoustic scenarios by varying the geometry, source location, and
inhomogeneous boundary in 2D and 3D. According to the singular energy decay,
the geometry of the room does not have a significant impact on choosing the
reduced basis and, therefore, the speedup. Moreover, the performance of the
ROM when adding new parameters can also be estimated based on previous
calculations, as it has been shown in Figure 4 that the singular value decay
is practically the same in the first basis functions when adding new
parameters to the ROM, which is a clear indication of the potential of ROM for
the type of applications studied.

When increasing the dimension of a space (from 2D to 3D), a higher speedup is
resulted for the same error values (Figure 6). Speedups of up to two orders of
magnitude are found for 2D and up to three orders for 3D simulations when
compared to FOM. These results agree with previous studies with homogeneous
boundary conditions [40].

A recent MOR in time domain room acoustic simulations using an automated
Krylov subspace algorithm reported a speedup of 11–36 without introducing
audible differences for a simple scenario [39]. The present study shows higher
speedups for a larger domain (considering the higher frequency). For the 2D
domain (cases 1) with $f_{u}=2$ kHz, a reduction of the degree of freedom from
DOF$=12039$ to $N_{rb}=158$ and $N_{rb}=752$ resulted in a speedup of $142$
($\epsilon_{rel}=0.35\%$) and $9$ ($\epsilon_{rel}=3.2\times 10^{-3}\%$)
respectively. Moreover, for the 3D cubic case, a reduction of the degree of
freedom from DOF$=35937$ to $N_{rb}=81$ and $N_{rb}=837$ resulted in a speedup
of $800$ ($\epsilon_{rel}=1.7\%$) and $47$ ($\epsilon_{rel}=5.8\times
10^{-3}\%$) respectively.

The difference in the reverberation between FOM and ROM has been quantified in
relation to $5\%$ JND of RT defined in [54]. Note that different studies show
that the perception of the reverberation varies depending on the sound decay
[55] and the nature of the stimuli [56, 57, 58, 59], where the JND range from
$3\%$-$20\%$.

##  V Conclusion

This study is concerned with the ability of ROM for use with different
boundary conditions under the variation of a considerable number of parameters
and an inhomogeneous distribution of absorption across the different surfaces
of a room. First, our results confirm that the RBM is more favourable in terms
of computational reduction for symmetric problems, e.g., source positioned at
the centre than in a corner. Second, results show that the singular value
decay becomes more gentle when including more parameters into the ROM.
Thirdly, speedups of one-two orders of magnitude are found for 2D, while two-
three orders of magnitude are found for 3D. Fourthly, an analysis of
reverberation time confirms that ROM produces IRs of which the RTs are less
than $5\%$ from FOM. A smaller number of basis modes is needed in the
truncated basis for the symmetric case in 3D, which shows a performance 365
times faster than the FOM.

It can be concluded that complex ROMs with a large number of parameters and
with acoustic materials distributed inhomogeneously behave similarly to simple
and homogeneous models and can achieve similar speedup performance. However,
non-symmetric source positions and geometries and a large number of parameters
can lead to a slower singular value decay, which may decrease the reduction if
a large number of basis functions are included in the ROM. Although the FOM
simulations are computationally costly, the price paid to create the ROM using
FOM simulations is worthy for multiple design evaluations, where different
parameter configurations are to be explored or optimized for room acoustics.

###### Acknowledgements.

This research was done at the Technical University of Denmark and partly
supported by Innovationsfonden, Denmark (Grant ID 9065-00115B), Rambøll
Danmark A/S and Saint-Gobain Ecophon A/S, Sweden.

## References

## References

  * [1] D. Abbot, “Calming the office cacophony,” Safety Health Pract. 22(1), 34–36 (2004).
  * [2] N. A. Oseland and A. Burton, “Quantifying the impact of environmental conditions on worker performance for inputting to a business case to justify enhanced workplace design features,” J. Build. Survey, Appr. Val. 1(2), 151–164 (2012).
  * [3] B. M. Shield and J. E. Dockrell, “The effects of noise on children at school: A review,” Build. Acoust. 10(2), 97–116 (2003).
  * [4] M. Klatte, J. Hellbrück, J. Seidel, and P. Leistner, “Effects of classroom acoustics on performance and well-being in elementary school children: A field study,” Environ. Behav. 42(5), 659–692 (2010).
  * [5] A. Seddigh, E. Berntson, F. Jönsson, C. B. Danielson, and H. Westerlund, “The effect of noise absorption variation in open-plan offices: A field study with a cross-over design,” J. Environ. Psych. 44, 34–44 (2015).
  * [6] G. W. Evans and D. Johnsson, “Stress and open-office noise,” J. Appl. Psych. 85(5), 779–783 (2000).
  * [7] D. Botteldooren, “Finite-difference time-domain simulation of low-frequency room acoustic problems,” J. Acoust. Soc. Am 98(6), 3302–3308 (1995).
  * [8] A. Craggs, “A finite element method for free vibration of air in ducts and rooms with absorbing walls,” J. Sound Vib 73(4), 568–576 (1994).
  * [9] F. Pind, A. P. Engsig-Karup, C. H. Jeong, J. S. Hesthaven, M. S. Mejling, and J. S. Andersen, “Time domain room acoustic simulations using the spectral element method,” J. Acoust. Soc. Am 145(6), 3299–3310 (2019).
  * [10] L. Savioja and U. P. Svensson, “Overview of geometrical room acoustic modeling techniques,” JASA 138(2), 708–730 (2015).
  * [11] J. S. Hesthaven, G. Rozza, and B. Stamm, _Certified Reduced Basis Methods for Parametrized Partial Differential Equations_ (Springer, 2016).
  * [12] A. Melander, E. S. . F. Pind, A. Engsig-Karup, C. Jeong, T. Warburton, N. Chalmers, and J. S. Hesthaven, “Massive parallel nodal discontinuous galerkin finite element method simulator for room acoustics,” Int. J. of High Performance Computing Applications (under review in 2022) (2020).
  * [13] N. Borrel-Jensen, A. Engsig-Karup, and C.-H. Jeong, “Physics-informed neural networks for one-dimensional sound field predictions with parameterized sources and impedance boundaries,” JASA Express Letters. 1(12), 122402 (2021).
  * [14] P. Benner, S. Gugercin, and K. Willcox, “A survey of projection-based model reduction methods for parametric dynamical systems,” SIAM 57(4), 10.1137/130932715 (2015).
  * [15] A. Quarteroni, G. Rozza, and A. Manzoni, “Certified reduced basis approximation for parametrized partial differential equations and applications,” J. Math Ind 1(3), . (2011).
  * [16] N. N. Cuong, K. Veroy, and A. T. Patera, _Certified real-time solution of parametrized partial differential equations_ (Handbook of Materials Modeling, 2005, pp. 1529-1564).
  * [17] G. Rozza, D. B. P. Huynh, and A. T. Patera, “Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial differential equations,” Archives of Computational Methods in Engineering 15(3), 229–275 (2008).
  * [18] M. Drohmann and K. Carlberg, “The ROMES method for statistical modeling of reduced-order-model error,” SIAM/ASA J. on Uncer. Quant. 3(1), 116–145 (2015).
  * [19] C. Prud’homme and A. Patera, “Reduced-basis output bounds for approximately parametrized elliptic coercive partial differential equations,” Computing and Visualization in Science 6(2-3), 147–162 (2004).
  * [20] A. Quarteroni, A. Manzoni, and F. Negri, _Reduced Basis Methods for Partial Differential Equations: An Introduction_ (Springer, 2015).
  * [21] G. Berkooz, P. Holmes, J. L. Lumleyand, and C. W. Rowley, _Turbulence, Coherent Structures, Dynamical Systems and Symmetry_ (Cambridge University Press, 1996).
  * [22] M. Ohlberger and S. Rave, “Reduced basis methods: Success, limitations and future challenges” , Proceedings Of The Conference Algoritmy, 1-12 (2016).
  * [23] Q. Aumanna, E. Deckersb, S. Jonckheerec, W. Desmetc, and G. Muller, “Automatic model order reduction for systems with frequency-dependent material properties,” Comput. Meth. Appl. Mech. Eng. 397, 115076 (2022).
  * [24] Q. A. G. Muller, “Predicting near optimal interpolation points for parametric model order reduction using regression models,” App. Math and Mech. 20(1), 10 (2021).
  * [25] X. Xiea, W. Wanga, K. Hea, and G. Lia, “Fast model order reduction boundary element method for large-scale acoustic systems involving surface impedance,” Comput. Meth. Appl. Mech. Eng. 400, 115618 (2022).
  * [26] M. Ganesh, J. Hesthaven, and B. Stamm, “A reduced basis method for electromagnetic scattering by multiple particles in three dimensions,” J. Comput. Physics 231(23), 7756–7779 (2012).
  * [27] Y. Chen and J. S. Hesthaven, “Certified reduced basis methods and output bounds for the harmonic maxwell’s equations,” SIAM J. Sci. Comput., 32(2), 970–996 (2010).
  * [28] D. Amsallem, J. Cortial, and C. Farhat, “Towards realtime computational-fluid-dynamics-based aeroelastic computations using a database of reduced-order information,” AIAA J., 48(9), 2029–2037 (2010).
  * [29] S. Giere, “Numerical and analytical aspects of pod-based reduced-order modeling in computational fluid dynamics,” PhD thesis, Free University of Berlin, Germany, 2016.
  * [30] M. A. Grepl, “Reduced-basis approximation and a posteriori error estimation for parabolic partial differential equations,” PhD thesis, Massachusetts Institute of Technology, United States, 2005\.
  * [31] R. S. Puri, D. Morrey, A. J. Bell, J. F. Durodola, E. B. Rudnyi, and J. G. Korvink, “Reduced order fully coupled structural– acoustic analysis via implicit moment matching,” Appl. Math. Mod 33(11), 4097–4119 (2009).
  * [32] U. Hetmaniuk, R. Tezaur, and C. Farhat, “Review and assessment of interpolatory model order reduction methods for frequency response structural dynamics and acoustics problems,” Int. J. Num. Meth. Eng. 90(13), 1636–1662 (2012).
  * [33] J. Herrmann, M. Maess, and L. Gaul, “Substructuring including interface reduction for the efficient vibro-acoustic simulation of fluid-filled piping systems,” Mech. Sys. Sig. Proc 24(1), 153–163 (2010).
  * [34] H. Basir, A. Javaherian, Z. Shomali, R. Firouz-Abadi, and S. Gholamy, “Acoustic wave propagation simulation by reduced order modelling,” Exploration Geophysics 49(3), 386–397 (2017).
  * [35] L. Borcea, J. Garnier, A. Mamonov, and J. Zimmerling, “Reduced order model approach for imaging with waves,” Inverse Problems 38(2), 0025004 (2022).
  * [36] J. Willhite, K. Frampton, and D. W. Grantham, “Reduced order modeling of head related transfer functions for virtual acoustic displays,” JASA 113(4), 2270–2270 (2003).
  * [37] S. Bukka, Y. Law, H. Santo, and R. Chan, “Reduced order model for nonlinear multi-directional ocean wave propagation,” Physics of Fluids 33(11), 117115 (2021).
  * [38] A. Tello and R. Codina, “Field-to-field coupled fluid structure interaction: A reduced order model study,” International Journal for Numerical Methods in Engineering 122(1), 53–81 (2021).
  * [39] M. Miller, S. Ophem, E. Deckers, and W. Desmet, “Time-domain impedance boundary conditions for acoustic reduced order finite element simulations,” Computer Methods in Applied Mechanics and Engineering. 387(15), 114173 (2021).
  * [40] H. S. Llopis, A. P. Engsig-Karup, C.-H. Jeong, F. Pind, and J. S. Hesthaven, “Reduced basis methods for numerical room acoustic simulations with parametrized boundaries,” JASA 152(2), 851–865 (2022).
  * [41] J. S. Hesthaven and T. Warburton, _Nodal Discontinuous Galerkin Methods—Algorithms_ (Springer, New York, 2008).
  * [42] H. Xu, C. Cantwell, C. Monteserin, C. Eskilsson, A. P. Engsig-Karup, and S. Sherwin, “Spectral/hp element methods: Recent developments, applications, and perspectives,” Journal of Hydrodynamics 30(1), 1–22 (2018).
  * [43] B. Cotté, P. Blanc-Benon, C. Bogey, and F. Poisson, “Time-domain impedance boundary conditions for simulations of outdoor sound propagation,” AIAA Journal 47, 10 (2009).
  * [44] R. Troian, D. Dragna, C. Bailly, and M.-A. Galland, “Broadband liner impedance eduction for multimodal acoustic propagation in the presence of a mean flow,” J. Sound Vib. 392, 200–216 (2017).
  * [45] Y. Miki, “Acoustical properties of porous materials—modifications of delany-bazley models,” J. Acoust. Soc. Jpn. 11(1), 19.24 (1990).
  * [46] J. F. Allard and N. Atalla, _Propagation of Sound in Porous Media: Modelling sound absorbing materials_ , 3rd ed. (Wiley, West Sussex, 2009).
  * [47] B. Gustavsen and A. Semlyen, “Rational approximation of frequency domain responses by vector fitting,” IEEE Trans. Power Delivery 14(3), 1052–1061 (1999).
  * [48] C. Bigoni and J. S. Hesthaven, “Simulation-based anomaly detection and damage localization: an application to structural health monitoring,” Comput. Meth. Appl. Mech. Eng 63, 12896 (2020).
  * [49] W. T. Weeks, “Numerical inversion of laplace transform using Laguerre functions,” J. Assc. Comp. Mach. 13(3), 419–429 (1966).
  * [50] L. Peng and K. Mohseni, “Symplectic model reduction of hamiltonian systems,” SIAM J. Sci. Comput. 38(1), A1–A27 (2016).
  * [51] S. Sakamoto, “Phase-error analysis of high-order finite difference time domain scheme and its influence on calculation results of impulse response in closed sound field,” Acoust. Sci.Tech 28(5), 295–309 (2007).
  * [52] C. Jeong and J. lh, “Effects of source and receiver locations in predicting room transfer functions by a phased beam tracing method,” J. Acoust. Soc. Am 131(5), 3864–3875 (2012).
  * [53] T. Cox and P. DAntonio, _Acoustic absorbers and diffusers : theory, design and application : Theory, design and application_ (Spon, 2004).
  * [54] “EN ISO 3382-1:2009, Measurement of room acoustic parameters, part 1: Performance spaces.,” Standard, International Organization for Standardization (2009).
  * [55] P. Luizard, B. F. G. Katz, and C. Guastavino, “Perceptual thresholds for realistic double-slope decay reverberation in large coupled spaces,” J. Acoust. Soc. Am 137(1), 75–84 (2015).
  * [56] M. G. Blevins, A. T. Buck, Z. Peng, and L. M. Wang, “Quantifying the just noticeable difference of reverberation time with band-limited noise centered around 1000 hz using a transformed up-down adaptive method” , in Proc. Int. Symp. Room Acoustics (ISRA) Toronto, Canada (2013).
  * [57] M. Karjalainen and H. Järveläinen, “More about this reverberation science: Perceptually good late reverberation” , In Proc. 111th Audio Eng. Soc. Conv., New York, USA (2001).
  * [58] Z. Meng, F. Zhao, and M. He, “The just noticeable difference of noise length and reverberation perception” , In Proc. Int. Symp. Com. Inf. Tech., Bangkok, Thailand (2006).
  * [59] K. Prawda, S. J. Schlecht, and V. Välimäki, “Improved reverberation time control for feedback delay networks” , In Proc. Int. Conf. Digital Audio Effects, Birmingham, UK (2019).

## Appendix

Table A1: Parameters selected to build the 2D ROM. Each row corresponds to a
FOM simulation with the corresponding parameters marked for each column whose
values are presented in Table 2.

| CEσ1 | CEσ2 | CEσ3 | CEd1 | CEd2 | CEd3 | CE${}_{d_{0}1}$ | CE${}_{d_{0}2}$ | CE${}_{d_{0}3}$ | FLZ1 | FLZ2 | FLZ3 | FLσ1 | FLσ2 | FLσ3 | WLσ1 | WLσ2 | WLσ3 | WRσ1 | WRσ2 | WRσ3
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
FOM1 | X |  |  |  | X |  |  | X |  |  |  |  |  | X |  |  | X |  |  | X |
FOM2 |  | X |  |  | X |  |  | X |  |  |  |  |  | X |  |  | X |  |  | X |
FOM3 |  |  | X |  | X |  |  | X |  |  |  |  |  | X |  |  | X |  |  | X |
FOM4 |  | X |  | X |  |  |  | X |  |  |  |  |  | X |  |  | X |  |  | X |
FOM5 |  | X |  |  | X |  |  | X |  |  |  |  |  | X |  |  | X |  |  | X |
FOM6 |  | X |  |  |  | X |  | X |  |  |  |  |  | X |  |  | X |  |  | X |
FOM7 |  | X |  |  | X |  | X |  |  |  |  |  |  | X |  |  | X |  |  | X |
FOM8 |  | X |  |  | X |  |  | X |  |  |  |  |  | X |  |  | X |  |  | X |
FOM9 |  | X |  |  | X |  |  |  | X |  |  |  |  | X |  |  | X |  |  | X |
FOM10 |  | X |  |  | X |  |  | X |  | X |  |  |  |  |  |  | X |  |  | X |
FOM11 |  | X |  |  | X |  |  | X |  |  | X |  |  |  |  |  | X |  |  | X |
FOM12 |  | X |  |  | X |  |  | X |  |  |  | X |  |  |  |  | X |  |  | X |
FOM13 |  | X |  |  | X |  |  | X |  |  |  |  | X |  |  |  | X |  |  | X |
FOM14 |  | X |  |  | X |  |  | X |  |  |  |  |  | X |  |  | X |  |  | X |
FOM15 |  | X |  |  | X |  |  | X |  |  |  |  |  |  | X |  | X |  |  | X |
FOM16 |  | X |  |  | X |  |  | X |  |  |  |  |  | X |  | X |  |  |  | X |
FOM17 |  | X |  |  | X |  |  | X |  |  |  |  |  | X |  |  | X |  |  | X |
FOM18 |  | X |  |  | X |  |  | X |  |  |  |  |  | X |  |  | X | X |  |  |
FOM19 |  | X |  |  | X |  |  | X |  |  |  |  |  | X |  |  | X |  | X |  |
FOM20 |  | X |  |  | X |  |  | X |  |  |  |  |  | X |  |  | X |  |  | X |
FOM21 |  | X |  |  | X |  |  | X |  |  |  |  |  | X |  |  | X |  |  |  | X

Table A2: Parameters selected to build the 3D ROM. Each row corresponds to a
FOM simulation with the corresponding parameters marked for each column whose
values are presented in Table 3.

| CEσ1 | CEσ2 | CEσ3 | CE${}_{d_{0}1}$ | CE${}_{d_{0}2}$ | CE${}_{d_{0}3}$ | WWσ1 | WWσ2 | WWσ3 | WN${}_{d_{0}1}$ | WN${}_{d_{0}2}$ | WN${}_{d_{0}3}$
---|---|---|---|---|---|---|---|---|---|---|---|---
FOM1 | X |  |  |  | X |  |  | X |  |  | X |
FOM2 |  | X |  |  | X |  |  | X |  |  | X |
FOM3 |  |  | X |  | X |  |  | X |  |  | X |
FOM4 |  | X |  | X |  |  |  | X |  |  | X |
FOM5 |  | X |  |  | X |  |  | X |  |  | X |
FOM6 |  | X |  |  |  | X |  | X |  |  | X |
FOM7 |  | X |  |  | X |  | X |  |  |  | X |
FOM8 |  | X |  |  | X |  |  | X |  |  | X |
FOM9 |  | X |  |  | X |  |  |  | X |  | X |
FOM10 |  | X |  |  | X |  |  | X |  | X |  |
FOM11 |  | X |  |  | X |  |  | X |  |  | X |
FOM12 |  | X |  |  | X |  |  | X |  |  |  | X


# State-of-Charge Aware EV Charging

Yize Chen and Baosen Zhang Y. Chen is with the Lawrence Berkeley National
Laboratory, email<EMAIL_ADDRESS>B. Zhang is with the University of
Washington, email<EMAIL_ADDRESS>
###### Abstract

Recent proliferation in electric vehicles (EVs) are posing profound impacts
over the operation of electrical grids. In particular, due to the physical
constraints on charging stations’ capacity and uncertainty in charging demand,
it becomes an emerging challenge to design high performance scheduling
algorithms to better serve charging sessions. In this paper, we design a
predictive charging controller by actively incorporating each EV’s state-of-
charge (SOC) information, which has strong effects on the utilization of
dispatchable power during peak hours. Simulation results on both synthetic and
real-world EV session and charging demand data demonstrate the proposed
algorithm’s benefits on maximizing charging throughput and achieving higher
rate of feasible charging sessions while satisfying battery and station
physical constraints at the same time.

##  I Introduction

EV proliferation is just unfolding worldwide. In the United States, a recent
report predicts under medium scenario, EV consumption will increase 30-fold
from 2020 to 2050 [1]. However, such rapid electrification of transportation
sector presents unprecedented challenge to the charging infrastructure. For
uncontrolled EV charging system, the unplanned charging sessions may lead to
transformer overloading, reduced equipment lifetime, low charging capacity
utilization ratio, and even bring distribution grid operational challenges
[2]. Compared to unmanaged charging, recent report found specifically designed
smart charging EVs could not only maximize the usage of charging facilities,
but also further reduce greenhouse gas emissions by additional $32\%$ [3]. By
properly controlling the charging rate for individual EV, system operator can
achieve the operational goals such as alleviating over-capacity demand during
peak hours, fulfilling more charging requests, or participating in demand-side
management and other utility services.

In order to realize highly flexible and reliable charging support for EVs,
researchers have looked into the scheduling and operation problem of a
charging station. [4] focused on earliest-deadline-first charging scheme, and
proved the average-cost optimality of such policy. If the full information
regarding each charging session is beforehand revealed to the charging station
operator, the offline optimal charging strategy can be derived by solving an
optimization problem [5, 6]. Online formulation and hardware interface were
introduced in [7]. [8, 9] designed pricing schemes to either optimize social
welfare or minimize system costs. There are also learning-based algorithms to
tackle the problem of charging session uncertainty and associated resource
allocation under congestion [10]. For a survey of EV charging algorithms, we
refer readers to [11].

Figure 1: The schematic of proposed algorithm on station-level EV charging
control. By explicitly modeling the relationship between peak charging rate
($\bar{u}_{i}(t)$) and the battery SOC, proposed algorithm is able to better
allocate available charging power to individual EVs while avoiding power loss.

However, in most of existing charging station models, the physical
characteristic of battery peak charging rate with respect to state-of-charge
(SOC) is not explicitly taken into consideration. For the typical battery in
EV, once the SOC is greater than a relatively small ratio (e.g., $15\%$), the
maximum allowed charging power decreases gradually with respect to SOC level
due to the electrochemical characteristics [11]. Overlooking such peak rate -
SOC curve may cause system operator allocating redundant and over-limit power
for individual EVs with high SOC. It is thus becoming critical to quantify how
battery SOC impact the peak charging rate, which in turn affects the
utilization of available power during peak hours.

In this paper, we take the perspective from operating a charging station with
unknown future charging demand. We explicitly take EV battery SOC into
account, and examine how to incorporate time-variant peak charging rate into
charging session scheduling. By focusing on modeling both the public charging
stations and EV charging sessions, we propose a fast predictive controller and
a tractable solving procedure. Our algorithm is illustrated in Fig. 3. At each
timestep, the controller determines the optimal charging rate for each EV in
order to fulfill operator’s objective such as maximizing total power delivered
to the vehicles while satisfying station and battery constraints. The key
insight comes from the flexibility provided by the heterogeneous deadlines of
charging tasks, while we can successfully defer charging load during peak
hours by modeling the variable peak charging rate as a function of battery
SOC. In the numerical experiments, we demonstrate that for different setups on
peak rate - SOC curve and charging demand, the proposed SOC-aware EV charging
scheduler can always accommodate more feasible charging tasks when charging
capacity is limited.

##  II EV Charging Model

In this section we describe charging system model, and illustrates how to
model both the EV level and station level constraints. Along with the
operator’s objective, we formulate the optimization problem for scheduling
charging sessions, and illustrate the importance of procuring charging
flexibility during peak hours.

###  II-A Charging Station Model

We consider a single charging station which can provide simultaneous charging
services for multiple EVs. Let $\mathscr{V}:=\\{1,2,3,...\\}$ be the set of
all EVs over an optimization horizon $\mathscr{T}:={1,...,T}.$ For each EV
charging session $i$, it is represented by a 5-element tuple
$(x_{i,initial},x_{i,final},t_{i,arrival},t_{i,depart},\bar{u}_{i})$, where
$x_{i,initial}$ is the initial state of the battery when arriving at the
station, and $x_{i,final}$ is user-specified required SOC before departure.
$t_{i,arrival}$ and $t_{i,depart}$ denote the arrival and departure time for
vehicle $i$ respectively. For each EV, there is also a hard constraint
$\bar{u}_{i}$ on the EV’s peak charging rate, which we will detail the
modeling of such variable.

The system operator of the EV charging station is interested in finding a
feasible schedule of charging rate $u_{i}(t),i\in\mathscr{V},t\in\mathscr{T}$
which can satisfy each customer’s charging needs. There are a set of
operational constraints which have to be satisfied considering both the
charging station and charging sessions physical limits:

####  II-A1 Charging session temporal constraints

For each EV $i$, the charging power is only available when the charging
session is live:

$u_{i}(t)=0,\quad t\notin[t_{i,arrival},t_{i,depart}),\;i\in\mathscr{V}.$ (1)

####  II-A2 Available charging power

At each timestep, there is a limit on the total power drawn from the charging
station:

$\sum_{i\in\mathscr{V}}u_{i}(t)\leq P(t),\quad t\in\mathscr{T}.$ (2)

We note that $P(t)$ can be a limit imposed by the parking garage
infrastructure such as transformer capacity or available renewable
generations. In this work, we assume there are adequate chargers for onsite
EVs.

####  II-A3 Valid charging rate

For each EV $i$, the charging power is nonnegative and always limited by the
maximum acceptable charging rate $\bar{u}_{i}(t)$:

$0\leq u_{i}(t)\leq\bar{u}_{i}(t),\;i\in\mathscr{V}.$ (3)

####  II-A4 EV battery state

For each EV with a live charging session, the battery state is initiated with
EV’s arrival state, and is updated based upon the current charging rate
$u_{i}(t)$:

$\displaystyle x_{i}(0)=$ $\displaystyle x_{i,initial},$ (4a) $\displaystyle
x_{i}(t)=$ $\displaystyle x_{i}(t-1)+\delta u_{i}(t),\;t\geq 1$ (4b)

where $\delta$ is a constant with unit of hours of sampling time intervals.

###  II-B Scheduling Objective

There are several different objectives possibly achieved by the charging
facility operator. In this work, we encourage maximizing the sum of each
customer’s utility throughout $\mathscr{T}$. Other possible choices include
following pre-defined load profile [7], reducing load variance with demand
response goals [5, 12], or minimizing system cost considering variable price
[8], while such objectives can be flexibly integrated into our scheduling
framework. The utility maximization problem (MPC) can be formulated as
follows:

$\displaystyle\max_{\mathbf{u}}\quad$
$\displaystyle\sum_{t\in\mathscr{T}}\sum_{i\in\mathscr{V}}\log
u_{i}(t)-\lambda\sum_{i\in\mathscr{V}}||x_{i}(T)-x_{i,depart}||_{2}^{2}$ (5a)
$\displaystyle s.t.\quad$
$\displaystyle\eqref{equ:time},\eqref{equ:sum_p},\eqref{equ:u},\text{and}\eqref{equ:state};$
(5b)

where $\mathbf{u}$ denotes the collection of
$u_{i}(t),\;i\in\mathscr{V},\;t\in\mathscr{T}$ and $\lambda$ is a weighting
parameter.

The objective (5a) is a strictly concave function, and it will incur a
proportionally fair vector for each timestep, which maximizes the sum of all
vechicle’s logarithmic utility functions. We choose the objective function to
motivate fair sharing of available $P(t)$ and meeting the terminal energy
demand at the same time. By solving (5) based on current availability of
charging power and connected EVs, we can find optimal charging action which
satisfy all operating constraints.

##  III SOC Aware Charging Scheduling

In this section, we describe how to integrate each vehicle’s SOC-dependent
peak charging rate into the design of charging sessions scheduling. We would
like to emphasize that by leveraging the SOC and peak charging rate
information, the proposed algorithm can manage to serve more EVs with higher
throughput.

One intrinsic limitation of previous work is treating the peak charging rate
$\bar{u}_{i}(t)$ as a time-invariant variable. Under such assumption, solving
the ordinary MPC problem (5) leads to inefficient use of charging capacity
especially when battery is approaching fully charged states, as the dispatched
$u_{i}(t)$ may exceed the actual peak charging rate. These unutilized energy
during congested hours can alternatively meet greater charging demands. To
overcome such limitation, we propose to model $\bar{u}_{i}(t)$’s dependency on
SOC using the following function:

$\bar{u}_{i}(t)=\bar{u}_{i}^{*}-\alpha_{i}\cdot x_{i}(t-1),\quad
i\in\mathscr{V}$ (6)

where we term $\alpha_{i}$ as a “decaying factor” to model the peak charging
rate decrease with respect to battery SOC. And $\bar{u}_{i}^{*}$ is the
nominal value for the maximum charging rate. Adding such evolution of
$\bar{u}_{i}(t)$ results to the following optimization problem (SOC_MPC):

$\displaystyle\max_{\mathbf{u}}\quad$
$\displaystyle\sum_{t\in\mathscr{T}}\sum_{i\in\mathscr{V}}\log
u_{i}(t)-\lambda\sum_{i\in\mathscr{V}}||x_{i}(T)-x_{i,depart}||_{2}^{2}$ (7a)
$\displaystyle s.t.\quad$
$\displaystyle\eqref{equ:time},\eqref{equ:sum_p},\eqref{equ:u},\eqref{equ:state},\text{and}\eqref{equ:bar_u}.$
(7b)

We also note that such modeling assumption on SOC is practical and can be
easily integrated into real-world EV charging stations. Current charging
station design realizes the possibility of collecting and utilizing vehicle’s
SOC information with both hardware and software infrastructure [13]. When an
EV plugs in, the vehicle sends information about its current SOC, so that
$\bar{u}_{i}(t)$ can be adjusted in real time for both operator’s and EV
owner’s benefits.

However, as $\bar{u}_{i}(t)$ is a function of $x_{i}(t)$, while $u_{i}(t)$ is
further constrained by $\bar{u}_{i}(t)$, optimization problem (7) can not be
solved directly using off-the-shelf solvers due to the nested nature of
constraints and time-varying bounds in inequality constraints. In [14], the
authors proposed a heuristic based method to solve the optimization problem
involving time-changing $\bar{u}(t)$, yet the optimality and solution
feasibility are not guaranteed. We will illustrate our solution procedure to
find the exact solution in the next subsection.

###  III-A Modeling of Peak Charging Rate under Variable SOC

Figure 2: The arrival time statistics for our synthetic data test case and ACN
real data. The charging peaks arrive in the morning.

###  III-B SOC Aware Adaptive Charging

For each vehicle $i$, to resolve the dependence between optimization variables
and time-varying $\bar{u}_{i}(t)$, we explicitly write out the relationship
between $\bar{u}_{i}(t)$ and $x_{i}(t)$ for all timestep:

$\begin{bmatrix}\bar{u}_{i}(1)\\\ \bar{u}_{i}(2)\\\ \vdots\\\
\bar{u}_{i}(T)\end{bmatrix}=\bar{u}_{i}^{*}\cdot\mathbf{1}-\alpha_{i}\begin{bmatrix}1&&&\\\
&1&&\\\ &&\ddots&\\\ &&&1\end{bmatrix}\begin{bmatrix}x_{i}(0)\\\ x_{i}(1)\\\
\vdots\\\ x_{i}(T-1)\end{bmatrix}$

Then we can further replace (3) in the vector form
$\mathbf{0}\leq\mathbf{u}_{i}\leq\mathbf{\bar{u}}_{i}$ where
$\mathbf{u}_{i},\mathbf{\bar{u}}_{i}\in\mathbb{R}^{T}$. The state evolution
function (4) can be also rolling out in the matrix form:

$\begin{bmatrix}x_{i}(1)\\\ x_{i}(2)\\\ \vdots\\\
x_{i}(T)\end{bmatrix}=\begin{bmatrix}x_{i}(0)\\\ x_{i}(1)\\\ \vdots\\\
x_{i}(T-1)\end{bmatrix}+\delta\begin{bmatrix}u_{i}(1)\\\ u_{i}(2)\\\ \vdots\\\
u_{i}(T)\end{bmatrix}$

In this way, we can explicitly model all the constraints on EV peak charging
rate, EV charging capacity and charging states, and thus we can find the
optimal solution for the SOC_MPC problem with low computation complexity.

Input: Charging session information

Output: Charging actions $\mathbf{u}(1),...,\mathbf{u}(T)$

Parameters: Charging station capacity $P(t)$

1 for  _t=0,…, T_ do

2 $\hat{\mathscr{V}_{t}}:=\\{i\in\mathscr{V}|t_{i,arrival}\leq
t\;\textbf{AND}\;t_{i,depart}>t\\}$;

3 $\\{u_{i}(t)\\}$=SOC_MPC$(x_{i}(t-1),P(t),t)$

4 Validate feasibility of $u_{i}(t)$;

5 $x_{i}(t)=x_{i}(t-1)+\delta u_{i}(t)$;

6 Check terminal constraint based on $x_{i}(t)$ and $x_{i,depart}$;

7Check EV departure and arrival information.

Algorithm 1 SOC Aware EV Charging

The detailed algorithm is listed in Algorithm 1. In our implementation, the
proposed SOC-aware charging scheme resolves the problem every time when there
is a new EV coming/leaving, or the terminal states for required energy are met
so that the list of active sessions
$\hat{\mathscr{V}}_{t}:=\\{i\in\mathscr{V}|t_{i,arrival}\leq
t\;\textbf{AND}\;t_{i,depart}>t\\}$ needs to be updated. At time $t$, the
predictive controller receive last step’s EV states and final required energy
$x_{i,depart}$, and implements the action $\mathbf{u}(t)$. We note that the
EV’s state evolution is always based on actual $\bar{u}_{i}(t)$, so that over-
the-limit charging power will not be used. We will show in the next section
such SOC information is valuable especially when there is a congestion in the
charging station (e.g., the total available charging power $P(t)$ can not meet
all EVs to charge at their full capacity), and can allocate limited power more
efficiently.

##  IV Simulations

In this section, we will evaluate the performance of proposed SOC-aware EV
charging mechanism (SOC_MPC), and compare it against other algorithms such as
equal share (ES) [15], earliest-deadline-first (EDF) [16, 4], and model
predictive control (MPC) [7] (without considering SOC information). We make
the code for EV station modeling and proposed algorithm publicly
available111github.com/chennnnnyize/State-Demand_Aware_EV_Charging.

Figure 3: Comparison on the total delivered energy of four methods on two
datasets and varying battery $\alpha$ setup. SOC_MPC can maximize the overall
throughput for all test cases.

###  IV-A Test Cases and Simulation Setup

We use both synthetic dataset and real-world EV charging station data from
ACN-Data [17]. In the synthetic data, we model the commuting pattern for a
workspace charging station using a Poisson process, and randomly sample the
initial SOC and required energy for each individual EV. For the ACN-Data, we
process the EV data in the week of September 5th to September 11th in 2021 for
the Caltech site. The charging session temporal distribution is illustrated in
Fig. 2. We choose charging station power capacity $P(t)$ such that offline
feasible solution can be found given the EV charging demand. We notice the
synthetic data has a sharper morning peak with regard to charging sessions.

We include a brief introduction of other implemented EV charging schemes for
method comparisons:

  * •

The ES algorithm allocates available charging power supply equally to all
active EV charging sessions at each timestep while satisfying each EV’s
minimum and maximum charging rate.

  * •

The EDF algorithm assigns $\bar{u}_{i}$ to EVs with earliest deadlines while
keeping the sum of charging power less or equal to $P(t)$.

  * •

The ordinary MPC algorithm iteratively solves (5) with a fixed $\bar{u}_{i}$,
and other settings are kept the same as SOC_MPC.

For all tested cases and implemented algorithms, we update the charging status
when either a EV connects, departs or the charging session terminal
constraints have been met. Without loss of generality, we choose simulation
interval as $15$ minutes, and set $\delta=1$ in all our simulation cases. We
look into three different settings on the peak charging rate-SOC parameter
$\alpha$222The choice of $\alpha$ is based on the battery capacity, and thus
we choose different set of values for the two datasets., and illustrate the
value of including $SOC$ information into efficient EV charging algorithm
design. We use CVXPY [18] to solve the optimization problem and find MPC
solutions.

Figure 4: The comparison of stepwise delivered energy between MPC and SOC_MPC
on the ACN dataset with $\alpha=0.1$.

###  IV-B Performance Evaluation

We first compare the amount of delivered energy using all four algorithms. In
all settings with different $\alpha$ for both datasets, we find proposed
algorithm can always deliver the most energy to EVs, and the average
performance gain with respect to ordinary MPC is $10.27\%$. ES and EDF
algorithms are not performing stably when faced with different settings of EV
arrival process and peak charging rate decay. The EDF algorithm are performing
well when the charging demand is low, such that charging sessions with
earliest deadline can be always satisfied. Once the charging demand are high
while EDF is dispatching too much power based on schedule priority, a portion
of the power is not fully utilized on the EVs with high SOC. The ES algorithm
is performing more evenly under different setting of $\alpha$, but can not
produce many feasible charging sessions with regard to the final energy
requests $x_{i,depart}$, since many sessions requesting higher quantity of
energy during peak hours are only allocated small amount of power. This
illustrates that when there are limited charging power capacity, EDF and ES
are not ideal policies in terms of maximizing the usage of charging power or
satisfying individual EV charging needs.

Performance gain achieved by SOC_MPC majorly comes from better utilization of
total available power $P(t)$ during peak hours. The awareness of SOC lets the
algorithm distributes the charging power to more EVs rather than allocating
large amount of power to single EV. This is further illustrated in Fig. 4,
where during peak hours with multiple simultaneous charging sessions, SOC_MPC
is able to make use of almost all of the available charging power. On the
other hand, MPC finds charging actions that are over the charging limit
$\bar{u}_{i}(t)$. Similar deficiencies are observed for ES and EDF algorithm,
as EDF and ES often allocate more-than-allowed power to batteries which are
almost fully charged.

Figure 5: The final battery states on the synthetic dataset with
$\alpha=0.05$.

We also compare the final states of batteries in Fig. 5 between the MPC and
SOC_MPC algorithm. We can observe that for almost all vehicles, the final
battery states by using the proposed algorithm are higher than using standard
MPC controller. This also shows the promise that SOC information can help
design algorithm that brings benefits for EV customers.

##  V Conclusion and Discussion

In this work, we show that state-of-charge information for EV batteries are
valuable for designing charging station algorithms in terms of maximizing
energy usage and satisfying more charging session requests. By explicitly
considering the time-varying constraints on EV peak charging rate with respect
to SOC, we design EV charging controller to better allocate limited charging
power particularly during peak hours. This calls for a more accurate modeling
of the battery characteristics especially the peak charging rate given
different SOC values. In the future work, we will investigate the predictive
model for the future demand, and analyze how to harness flexibility as grid
resources via smart charging interfaces.

## References

  * [1] T. T. Mai, P. Jadun, J. S. Logan, C. A. McMillan, M. Muratori, D. C. Steinberg, L. J. Vimmerstedt, B. Haley, R. Jones, and B. Nelson, “Electrification futures study: Scenarios of electric technology adoption and power consumption for the united states,” National Renewable Energy Lab.(NREL), Golden, CO (United States), Tech. Rep., 2018.
  * [2] J. A. P. Lopes, F. J. Soares, and P. M. R. Almeida, “Integration of electric vehicles in the electric power system,” _Proceedings of the IEEE_ , vol. 99, no. 1, pp. 168–183, 2010.
  * [3] T. Lipman, D. Callaway, T. Peffer, and A. von Meier, “Open-source, open-architecture softwareplatform for plug-inelectric vehicle smartcharging in california,” 2020.
  * [4] E. Bitar and Y. Xu, “Deadline differentiated pricing of deferrable electric loads,” _IEEE Transactions on Smart Grid_ , vol. 8, no. 1, pp. 13–25, 2016\.
  * [5] L. Gan, U. Topcu, and S. H. Low, “Optimal decentralized protocol for electric vehicle charging,” _IEEE Transactions on Power Systems_ , vol. 28, no. 2, pp. 940–951, 2012.
  * [6] Z. Ma, D. S. Callaway, and I. A. Hiskens, “Decentralized charging control of large populations of plug-in electric vehicles,” _IEEE Transactions on control systems technology_ , vol. 21, no. 1, pp. 67–78, 2011.
  * [7] Z. J. Lee, D. Chang, C. Jin, G. S. Lee, R. Lee, T. Lee, and S. H. Low, “Large-scale adaptive electric vehicle charging,” in _2018 IEEE International Conference on Communications, Control, and Computing Technologies for Smart Grids (SmartGridComm)_. IEEE, 2018, pp. 1–7.
  * [8] M. Alizadeh, H.-T. Wai, M. Chowdhury, A. Goldsmith, A. Scaglione, and T. Javidi, “Optimal pricing to manage electric vehicles in coupled power and transportation networks,” _IEEE Transactions on control of network systems_ , vol. 4, no. 4, pp. 863–875, 2016.
  * [9] Z. J. Lee, J. Z. Pang, and S. H. Low, “Pricing ev charging service with demand charge,” _Electric Power Systems Research_ , vol. 189, p. 106694, 2020.
  * [10] A. Al Zishan, M. M. Haji, and O. Ardakanian, “Adaptive control of plug-in electric vehicle charging with reinforcement learning,” in _Proceedings of the Eleventh ACM International Conference on Future Energy Systems_ , 2020, pp. 116–120.
  * [11] J. C. Mukherjee and A. Gupta, “A review of charge scheduling of electric vehicles in smart grid,” _IEEE Systems Journal_ , vol. 9, no. 4, pp. 1541–1553, 2014.
  * [12] L. Zhang, V. Kekatos, and G. B. Giannakis, “Scalable electric vehicle charging protocols,” _IEEE Transactions on Power Systems_ , vol. 32, no. 2, pp. 1451–1462, 2016.
  * [13] “BMW chargeforward,” https://www.bmwchargeforward.com, accessed: 2021-08-12.
  * [14] Y. Cao, S. Tang, C. Li, P. Zhang, Y. Tan, Z. Zhang, and J. Li, “An optimized ev charging model considering tou price and soc curve,” _IEEE Transactions on Smart Grid_ , vol. 3, no. 1, pp. 388–393, 2011.
  * [15] J. Kay and P. Lauder, “A fair share scheduler,” _Communications of the ACM_ , vol. 31, no. 1, pp. 44–55, 1988.
  * [16] J. A. Stankovic, M. Spuri, K. Ramamritham, and G. C. Buttazzo, _Deadline scheduling for real-time systems: EDF and related algorithms_. Springer Science & Business Media, 2012, vol. 460.
  * [17] Z. J. Lee, T. Li, and S. H. Low, “Acn-data: Analysis and applications of an open ev charging dataset,” in _Proceedings of the Tenth ACM International Conference on Future Energy Systems_ , 2019, pp. 139–149.
  * [18] S. Diamond and S. Boyd, “Cvxpy: A python-embedded modeling language for convex optimization,” _The Journal of Machine Learning Research_ , vol. 17, no. 1, pp. 2909–2913, 2016.