# Privacy-preserving smart meter control strategy including energy storage
losses

Ramana R. Avula, Tobias J. Oechtering and Daniel Månsson  School of Electrical
Engineering and Computer Science
KTH Royal Institute of Technology, Stockholm, Sweden

###### Abstract

Privacy-preserving smart meter control strategies proposed in the literature
so far make some ideal assumptions such as instantaneous control without
delay, lossless energy storage systems etc. In this paper, we present a one-
step-ahead predictive control strategy using Bayesian risk to measure and
control privacy leakage with an energy storage system. The controller
estimates energy state using a three-circuit energy storage model to account
for steady-state energy losses. With numerical experiments, the controller is
evaluated with real household consumption data using a state-of-the-art
adversarial algorithm. Results show that the state estimation of the energy
storage system significantly affects the controller’s performance. The results
also show that the privacy leakage can be effectively reduced using an energy
storage system but at the expense of energy loss.

###### Index Terms:

Smart meter privacy, Bayesian hypothesis testing, partially observable Markov
decision process (PO-MDP), energy storage losses, dynamic programming

##  I Introduction

A smart grid (SG) is a next-generation energy network with capabilities to
improve grid reliability and efficiency of power generation and distribution
with smooth integration of renewable energy sources. In this automated
network, a smart meter (SM) is a crucial component which measures the energy
consumption of the user and transmits the readings to the utility provider at
regular intervals of time. This raises privacy concerns [1] since high-
resolution readings can allow anyone who has access to this data to infer
about consumer’s behavior. Since its introduction in [2], non-intrusive load
monitoring (NILM) techniques are known to be quite effective in disaggregating
the smart meter readings and thereby detecting the states of most of the
general types of household appliances [3]. A comparative study was done in
[3], which shows that the existing state of the art NILM algorithms are
capable of achieving detection accuracy up to 99% for certain appliance types,
which is quite concerning in the privacy context.

Addressing this issue, several privacy-preserving techniques have been
proposed in the literature, which are surveyed in [4, 5]. Secure communication
and cryptographic approaches [6, 7, 8] may succeed in preventing the
unauthorized third party access, but they would fail to protect the consumer
privacy from a greedy authorized or compromised utility provider. A promising
physical layer privacy approach is load signature moderation (LSM), where an
energy storage system (ESS) is used to moderate the consumer’s load profile in
order to hide appliances’ usage information. LSM using rechargeable battery
has previously been investigated in [9, 10, 11, 12, 13, 14] to obtain optimal
control strategy under different privacy settings. However, all the works so
far make some ideal assumptions such as instantaneous control without delay,
lossless ESS etc. These idealized strategies may provide theoretic performance
bounds but the feasibility of such strategies in practical situations must be
further investigated.

Figure 1: Schematic of the proposed smart metering system where the energy
management unit controls privacy leakage to an adversary by using energy
storage system with a model describing its losses and one-step-ahead
predictive control.

In this paper, we present a one-step-ahead predictive control scheme modeled
in a PO-MDP framework using an ESS. Similar to [13], we use a privacy metric
based on Bayesian risk. The overview of the proposed system is shown in Fig.
1. In this work, we restrict our analysis to the electrochemical battery as an
energy storage technology. Nonetheless, the same approach can be followed for
other storage technologies by modeling them as their equivalent electrical
circuits [15]. For a battery system, we present a model describing its losses
in power conversion, losses due to internal resistance and self-dissipation.
To the best of our knowledge, this is the first work to consider the non-
idealities in ESSs in the context of smart meter privacy.

The rest of the paper is organized as follows. In Section II, we present a
model for ESS considering the steady state energy losses. We also present the
charge and discharge bounds of ESS and also quantify the energy loss
associated with a discrete control action. In Section III, we present an
overview of the system along with the control strategy. In Section IV, we
evaluate the performance of the controller with real household data using a
state-of-the-art NILM algorithm. Lastly, we conclude the paper in Section V.

##  II Energy Storage system Model

Since batteries stores the energy as chemical potential in their electrodes,
it can only be interfaced with a DC (Direct Current) system. Hence, power
converters are needed to integrate the battery with an AC (Alternating
Current) system. The battery along with the power converters form the ESS. In
this work, we model the ESS using three simple electrical circuits as shown in
Fig. 2 to account for the steady state energy losses. Even though several
other processes of the ESS such as capacity fade, increase in internal
resistance, temperature dependence etc., can also be considered, as a first
step, we restrict our focus to the steady state energy losses. In the
following, we present our analysis of the three-circuit model in more detail.

Figure 2: Three-circuit energy storage system model

###  II-A Losses due to self-dissipation

Self-dissipation occurs even if the ESS is not connected to any load. Similar
to [16], we model this phenomenon using an RC circuit as shown in Fig. 2(a).
The capacitor C holds the charge content, Q of the battery and dissipates
through a parallel resistor R. In this circuit, the power converter and load
are together represented as a current source controlled by the current flowing
in the second circuit shown in Fig. 2(b). Given the self-discharge rate
$\gamma$ and a constant current $\text{I}_{\text{bat}}$ flowing into the
battery, the charge content of the battery is updated as

$\displaystyle{\text{Q}}_{t+\Delta
t}=(1-\gamma)\cdot{\text{Q}}_{t}+\beta\cdot{{\text{I}}_{\text{bat}}}$ (1)

where,

$\displaystyle\gamma=1-\exp\Big{(}\frac{-\Delta
T}{\text{RC}}\Big{)};\qquad\beta=\frac{-\gamma\Delta T}{\log(1-\gamma)}$ (2)

###  II-B Losses due to internal resistance

Similar to [16], the losses that occur in battery during its charging and
discharging operations are modeled using a series resistor $r$, as shown in
Fig. 2(b). The open circuit voltage of the battery, $\text{V}_{\text{OC}}$ is
represented as a voltage source controlled by the capacitor charge Q in the
first circuit shown in Fig. 2(a). For an input power of P from the power
converter, the current flowing into the battery is given as

${\text{I}}_{\text{bat}}=\frac{\sqrt{\text{V}_{\text{OC}}^{2}+4\text{rP}}-\text{V}_{\text{OC}}}{2\text{r}}$
(3)

###  II-C Losses due to power converters

As shown in Fig. 2(c), we model power converters as elements with a constant
efficiency factor within their operating region. For an input power of D from
the load, the power at the battery terminals can be written as

P $\displaystyle=\text{D}\cdot\left(\eta_{c}\mathds{1}\\{\text{D}\geq
0\\}+\eta_{d}^{-1}\mathds{1}\\{\text{D}<0\\}\right)$
$\displaystyle=\text{D}\cdot\delta(\text{D})$ (4)

where $\eta_{c},\eta_{d}$ are the efficiency factors of AC-DC and DC-AC
converters respectively and $\mathds{1}\\{A\\}$ is equal to 1 if A is true,
and 0 otherwise. $\delta(\text{D})$ is the common factor for both operations.

###  II-D Three-circuit ESS model

Integrating the three circuits by combining (1), (3) and (II-C), the
controller updates the energy state of the battery evolving over time using
the equation given as

${\text{Z}}_{t+\Delta
t}=(1-\gamma){\text{Z}}_{t}+\frac{\beta\text{V}_{\text{OC}}}{2\text{r}}\left(\sqrt{\text{V}_{\text{OC}}^{2}+4\text{r}\text{D}_{t}\delta_{t}}-\text{V}_{\text{OC}}\right)$
(5)

where $\gamma$, $r$ are time-invariant parameters, $\beta$ depends on the time
step $\Delta t$ and $\delta_{t}$ depends on the control variable
$\text{D}_{t}$. By limiting the battery current to $\text{I}_{\text{max}}$,
from (3), we have the control space limited as

$\displaystyle\text{D}_{\text{max}}$
$\displaystyle=\frac{1}{4\text{r}\eta_{d}}\left({(\text{V}_{\text{OC}}+2\text{rI}_{\text{max}})^{2}-\text{V}_{\text{OC}}^{2}}\right)$
(6) $\displaystyle\text{D}_{\text{min}}$
$\displaystyle=\frac{\eta_{c}}{4\text{r}}\left({\text{V}_{\text{OC}}^{2}-(\text{V}_{\text{OC}}-2\text{rI}_{\text{max}})^{2}}\right)$
(7)

Due to the finite energy capacity of the battery $\text{Z}_{max}$, from (5) we
have the following constraints on $\text{D}_{t}$:

$\displaystyle\text{D}_{t,max}$
$\displaystyle=\frac{\text{V}_{\text{OC}}^{2}}{4\text{r}\eta_{d}}\left(\left(\frac{2\text{r}\big{[}\text{Z}_{max}-(1-\gamma){\text{Z}}_{t}\big{]}}{\beta\text{V}_{\text{OC}}^{2}}+1\right)^{2}-1\right)$
(8) $\displaystyle\text{D}_{t,min}$
$\displaystyle=\frac{\eta_{c}\text{V}_{\text{OC}}^{2}}{4\text{r}}\left(\left(\bigg{[}\frac{-(1-\gamma)2\text{r}{\text{Z}}_{t}}{\beta\text{V}_{\text{OC}}^{2}}+1\bigg{]}^{+}\right)^{2}-1\right)$
(9)

where $[x]^{+}$ is equal to $x$ if $x\geq 0$, and 0 otherwise. In comparision
to (5), the energy state of an ideal lossless battery evolves over time as

${\text{Z}}_{t+\Delta t,ideal}={\text{Z}}_{t}+\text{D}_{t}\cdot\Delta t$ (10)

Using (5) and (10), the energy loss associated with a discrete control action
can be given as

$\mathcal{E}_{loss}({\text{Z}}_{t},\text{D}_{t})={\text{Z}}_{t}+\text{D}_{t}\cdot\Delta
t-{\text{Z}}_{t+\Delta t}$ (11)

##  III System overview and Control strategy

The proposed smart metering system uses an ESS for load signature moderation.
The ESS can be placed either in series or in parallel configurations, in
between the SM and house as shown in Fig. 3. Both these configurations have
been used in the literature for SM privacy. Under ideal assumptions, the two
configurations are equivalent. However, considering the energy losses, we have
the following proposition.

######  Proposition 1.

The average energy loss in the parallel configuration is strictly less than
that of series configuration.

###### Proof.

Let X be the average energy demand by the house and Y be the average energy
request from the grid. Assuming that the energy from the battery is not
discharged back into the grid, only $({\text{Y}}-{\text{X}})$ flows through
the ESS components in parallel configuration, however, the total energy Y from
the grid flows through the ESS components in series configuration, leading to
higher energy losses. ∎

Figure 3: Placement of ESS between smart meter and house.

In order to reduce the energy losses, we consider a system with ESS in the
parallel configuration as shown in Fig. 1. The discrete time system is
controlled for every time slot $k$ within a finite time horizon
$\\{1,2,\dots,N\\}$. Each time slot $k$ is of a fixed time duration $T$. Let
$e$ and $q$ be the resolution of energy and power measurements respectively.
In the following analysis, we use the capital letters to denote random
variables, their realizations by the lower-case letters and the range space by
calligraphic letters.

For each time slot $k$, $X_{k}$ denotes the aggregate power drawn by all the
appliances in the house and is defined on $\cal{X}$$=\\{0,q,2q,$
$\dots,x_{max}\\}$. $Z_{k}$ defined on $\cal{Z}$$=\\{0,e,2e,\dots,z_{max}\\}$
denotes the energy available in the battery. The power drawn by the ESS is
denoted as $D_{k}$ and it is the control variable which is defined on
$\mathcal{D}=\\{-d_{min},\dots,-q,0,q,\dots,d_{max}\\}$, where $d_{min}$ and
$d_{max}$ are the maximum discharge and charge power of the ESS respectively,
given in (6) and (7). $D_{k}^{*}$ defined on $\cal{D}$ denotes the desired
battery power consumption scheduled by the energy management unit (EMU). In
the presence of an ESS, the SM records the aggregate power demands of consumer
and ESS. In this work, we allow the energy from the battery to be discharged
to the grid resulting in negative values of SM measurements. It is represented
by the random variable $Y_{k}=X_{k}+D_{k}$ which is defined on
$\mathcal{Y}=\\{-d_{min},\dots,-q,0,q,\dots,x_{max}+d_{max}\\}$. $H_{k}$
defined on $\cal{H}$ denotes the $n$-ary joint hypothesis of all the
appliances in the house and $\hat{H}_{k}$ defined on $\cal{H}$ denotes the
hypothesis detected by the adversary having access to the consumer’s
statistical and real-time data as well as the control strategy employed by the
EMU.

###  III-A Bayesian risk

Similar to [13], we use detection-theoretic approach by formulating the smart
meter privacy problem into an adversarial Bayesian hypothesis testing where an
adversary having access to the consumer’s statistical and real-time data tries
to make a guess on the hypothesis state using a decision strategy. In the
Bayesian formulation, each of the hypothesis test outcomes is assigned a cost
and the decision strategy that minimizes the average decision-making cost will
be employed by the adversary [17]. The average cost or _Bayesian risk_
function, $\mathcal{R}$, is given as

$\mathcal{R}_{k}=\sum_{i,j\in\mathcal{H}^{2}}C_{i,j}\cdot P(\hat{H}_{k}=i\mid
H_{k}=j)\cdot P(H_{k}=j)$ (12)

where $C_{i,j}$ is the cost of deciding $\hat{H}_{k}=i$ when $H_{k}=j$ is
true. By setting the cost of a correct decision to zero and the cost of an
error to unity, the risk function gives the average error probability of an
adversarial detection strategy.

In this work, the _accumulated minimum Bayesian risk_ (AMBR) is chosen as a
privacy metric, which is given as

$\text{AMBR}=\sum_{k=1}^{N}\mathcal{R}_{k}^{*}$ (13)

where $\mathcal{R}_{k}^{*}=\min\\{\mathcal{R}_{k}\\}$. The AMBR is a good
choice for measuring privacy due to its operational meaning. It explicitly
characterizes the best possible detection performance achievable by any
adversary.

###  III-B Control strategy

Ideally, the controller uses all the information available until time $k$
(denoted as $\mathcal{I}_{k}$) to choose an action $d_{k+1}$. However, as
described in [18], since $\mathcal{I}_{k}$ is increasing in dimension with
$k$, its sufficient statistic given by the posterior distribution of the
Markov chain $H_{k}$ conditioned on $\mathcal{I}_{k}$ (denoted as $\pi_{k}$)
is used instead of $\mathcal{I}_{k}$ to choose the action $d_{k+1}$. For a
given initial battery state $z_{0}$, this posterior distribution forms a
_information state_ or _belief state_ at time $k$, given as

$\displaystyle\pi_{k}(i)$ $\displaystyle={P}(H_{k}=i\mid\mathcal{I}_{k})$
$\displaystyle={P}(H_{k}=i\mid\pi_{k-1},x_{k})$ (14)

where
$\mathcal{I}_{k}=\\{z_{0},\pi_{0},x_{1},y_{1},\pi_{1},\dots,\pi_{k-1},x_{k}\\}$.
The control system is modeled as a PO-MDP controlled sensor, as shown in Fig.
4, by making the following assumptions:

  * •

The hypothesis of the house $H_{k}$ evolves over time following a first-order
Markov chain with a time-invariant transition probability $P_{H_{k}\mid
H_{k-1}}$.

  * •

The controller observes the Markov chain $H_{k}$ only through a noisy
measurement $X_{k}$ made with a time-invariant observation probability
$P_{X_{k}\mid H_{k}}$.

  * •

The control signal $D_{k}^{*}$ is generated by the controller using time-
dependent control strategy $P_{Y_{k}\mid X_{k-1},Z_{k-1},\Pi_{k-1}}$.

Figure 4: EMU modelled as a PO-MDP controlled sensor.

Given the initial energy state, $z_{0}$, the controller estimates the state of
ESS at any time $k$ using the equation

$\displaystyle z_{k}$ $\displaystyle=f(z_{k-1},d_{k})$ (15)

where $d_{k}=y_{k}-x_{k}$ and $f$ is a deterministic function of ESS model
given by (5). As described in [17], the minimum Bayesian risk function based
on our assumptions is given as

$\displaystyle\mathcal{R}_{k}^{*}(\pi_{k-1},z_{k-1},\mu_{k})=\sum_{y\in\mathcal{Y}}\min_{\hat{h}\in\mathcal{H}}\Bigl{\\{}\sum_{g,h,x\in\mathcal{H}^{2}\times\mathcal{X}}C(\hat{h},h)\cdot$
$\displaystyle\qquad\qquad\qquad{P}_{Y_{k}\mid X_{k-1},Z_{k-1}}(y\mid
x,z)\cdot{P}_{X_{k-1}\mid H_{k-1}}({x\mid g})\cdot$
$\displaystyle\qquad\qquad\qquad{P}_{H_{k}\mid H_{k-1}}({h\mid
g})\cdot{P}_{H_{k-1}}({g})\Bigr{\\}}$

For the finite time horizon, the optimal control strategy is the solution to
the nonlinear optimization problem with objective function given as

$\displaystyle\mu^{*}=\operatorname*{argmax}_{\\{\mu_{1},\cdots,\mu_{N}\\}}\sum_{k=1}^{N}\mathcal{R}_{k}^{*}(\pi_{k-1},z_{k-1},\mu_{k})$
(16)

subject to the constraints given as

$\displaystyle P_{Y_{k}\mid Z_{k-1}}(y_{k}\mid z_{k-1})=0\text{ if
}\begin{cases}y_{k}<{d}_{k,min}\\\ \qquad\text{ or}\\\
y_{k}>{d}_{k,max}+x_{max}\end{cases}$

The _belief state space_ (denoted as $\Pi$) is a $|\mathcal{H}|-1$ dimensional
unit simplex. Solving this optimization problem requires discretization of
$\Pi$ in order to get a finite set. The optimization variable in (16) is of
dimension $N\times|\mathcal{Y}|\times|\mathcal{X}|$ and solving it in its
original form is computationally complex as the dimensionality of the the
problem increases with $N$. This can be formulated into a recursive dynamic
programming problem as given in the following proposition, the proof of which
follows from [18].

######  Proposition 2.

For the finite horizon PO-MDP with model given in Section III, the optimal
control strategy $\mu^{*}=\\{\mu_{1}^{*},\mu_{2}^{*},$ $\dots,\mu_{N}^{*}\\}$
is the solution to the following backward recursion: Initialize
$\mathcal{V}_{N}(\pi,z)$ and then for $k=N-1,\dots,1$ iterate

$\displaystyle
L_{k}(\pi_{k-1},z_{k-1},\mu_{k})=\mathcal{R}_{k}^{*}(\pi_{k-1},z_{k-1},\mu)~{}+$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\sum_{x_{k},y_{k}\in\mathcal{X}\times\mathcal{Y}}\mathcal{V}_{k+1}(\pi_{k},z_{k})\cdot
P(x_{k},y_{k})$
$\displaystyle\mathcal{V}_{k}(\pi_{k-1},z_{k-1})=\max_{\mu_{k}}\Bigl{\\{}L_{k}(\pi_{k-1},z_{k-1},\mu_{k})\Bigr{\\}}$
$\displaystyle\delta_{k}^{*}(\pi_{k-1},z_{k-1})=\operatorname*{argmax}_{\mu_{k}}\Bigl{\\{}L_{k}(\pi_{k-1},z_{k-1},\mu_{k})\Bigr{\\}}$

∎

With the designed optimal strategy $\mu^{*}$, a real-time PO-MDP controller is
implemented as shown in the following algorithm.

Algorithm 1 Realtime PO-MDP controller

1:$\pi_{0},z_{0}$

2:for $k=1$ to $N$ do

3:Pre-process :

4: Choose action $y_{k}^{*}=\mu_{k}^{*}(\pi_{k-1},z_{k-1})$

5:ESS control :

6: if ($y_{k}^{*}<x_{k}+d_{k,min}$) then

7: Limit $y_{k}=x_{k}+d_{k,min}$

8: else if ($y_{k}^{*}>x_{k}+d_{k,max}$) then

9: Limit $y_{k}=x_{k}+d_{k,max}$

10: else

11: Allow $y_{k}=y_{k}^{*}$

12: end if

13:Post-process :

14: Update the belief state $\pi_{k}=T(\pi_{k-1},x_{k})$

15: Update the ESS state $z_{k}=f(z_{k-1},y_{k}-x_{k})$

16:end for

##  IV Numerical Experiments

The simulation experiments to validate our control scheme are implemented in
MATLAB using real household consumption data from ECO reference dataset [19].
The control strategy is obtained by solving the optimization problem using the
nonlinear programming solver based on interior point algorithm [20]. We
simulate a scenario where the controller is tasked to protect the events of a
water kettle every day between 8 AM and 9 AM. The controller chooses an action
every minute by observing the real-time appliance consumption data. For this
objective, a 12V 100Ah lithium-ion battery is selected, which can sufficiently
satisfy the power requirements of the kettle. To simplify the problem, we
assume a fixed $\text{V}_{\text{OC}}$ equal to the nominal battery voltage.
The parameters used in the simulation are listed in Table II. The Markov chain
probabilities of the PO-MDP control model are estimated from 30 days of
labeled training data and listed in the Table II.

###  IV-A Visualization of control actions

With this setup, the control actions for different initial states of the
battery are simulated and are shown in Fig. 7. Due to the measurement
quantization, switching events are noticed as peaks in the smart meter
measurements as shown in Fig. 7. These residual peaks are informative to an
adversary operating with high precision measurements. However, the cardinality
of the state space increases with measurement precision, which increases the
dimensionality of the optimization problem by $\mathcal{O}(n^{2})$. Fig. 7
shows the evolution of the state of charge (SOC) of the battery due to control
actions. It is interesting to notice that without any design objective on the
battery state, the control scheme is steering the battery towards the full
charge state. This result in the degradation of controller’s performance which
is discussed in the following.

TABLE I: Simulation parameters

Parameter | Symbol | Value
---|---|---
Max. appliance power demand |  $x_{\text{max}}$ (W) | 1700
Time slot duration | T (s) | 60
Time horizon length | N | 60
Power measurement resolution | q (W) | 500
Energy measurement resolution | e (Wh) | 5
Battery nominal voltage |  $\text{V}_{\text{nom}}$ (V) | 12
Battery capacity |  $\text{Q}_{\text{max}}$ (Ah) | 100
Max. allowed battery charging current |  $\text{I}_{\text{max}}$ (A) | 80
Max. allowed battery discharging current |  $\text{I}_{\text{min}}$ (A) | 80
Battery internal resistance | r ($\Omega$) | 0.006
Battery self-discharge rate |  $\gamma$ (%/month) | 3
Power converter efficiency |  $\eta_{c},\eta_{d}$ (%) | 95
Cardinality of $\mathcal{X}$ | $|\mathcal{X}|$ | 4
Cardinality of $\mathcal{Y}$ | $|\mathcal{Y}|$ | 8
Cardinality of $\mathcal{Z}$ | $|\mathcal{Z}|$ | 241
Cardinality of $\mathcal{H}$ | $|\mathcal{H}|$ | 2
Cardinality of ${\Pi}$ | $|{\Pi}|$ | 11
Max. allowed battery input power | $\text{D}_{\text{max}}$ | +2q
Min. allowed battery input power | $\text{D}_{\text{min}}$ | -2q
ESS model parameter | $\beta$ | 0.017

TABLE II: PO-MDP control parameters

Parameter | Value
---|---
$C_{i,j}$ | $\begin{bmatrix}0&1\\\ 1&0\end{bmatrix}$
$\pi_{0}$ | $\begin{bmatrix}0.95&0.05\end{bmatrix}^{T}$
$P(H_{k}|H_{k-1})$ | $\begin{bmatrix}0.98&0.34\\\ 0.02&0.65\end{bmatrix}$
$P(X_{k}|H_{k})$ | $\begin{bmatrix}1&0&0&0\\\ 0&0.17&0.14&0.17\end{bmatrix}^{T}$

###  IV-B Evaluation of controller’s performance

To evaluate the performance of the controller, we simulate an adversary using
NILM algorithm. In particular, we simulated Weiss’ algorithm [21] which
extracts switching events from the aggregate smart meter data and assigns each
event to the appliance with the best match in a signature database. This
algorithm is implemented using NILM toolbox developed by [19]. We use 30 days
of labeled training data to create the signature database. Weiss’ algorithm
utilizes three-dimensional consumption data (i.e., real, reactive, and
distortion powers) in order to match an event signature. We tested its
detection performance by injecting the controlled battery current in-phase to
the supply voltage resulting in corrupted real power measurements. The
accuracy of the adversarial detection is measured using F-score, which is
given as

$\text{F-score}=\frac{1}{1+(\text{FN}+\text{FP})/(2\text{TP})}$ (17)

where FN, FP and TP denote false negative, false postive and true positive
respectively. The F-score lies between 0 and 1, where F-score = 0 indicates no
detection and F-score = 1 indicates complete detection.

Figure 5: Control actions of battery.

Figure 6: Smart meter readings vs household consumption.

Figure 7: Evolution of SOC of battery due to control actions.

The Weiss’s algorithm is simulated under different test conditions using 30
days of validation data and the obtained average F-scores, energy losses and
the AMBR are listed in Table III. It can be seen from the results that the
AMBR and the F-score are correlated. The test case without a battery resulted
in the highest F-score. While using a battery, the controller is able to
reduce the F-score significantly. However, it can be observed that the ability
of the controller to preserve privacy depends on the initial state of the
battery. For this simulation setup, the controller performed better when
initialized with a battery of 25%-50% SOC compared to full charge. This
indicates that if the battery state is steered towards 25%-50% SOC by the end
of the control time horizon, it would result in better performance for the
next control horizon. However, this improved performance is achieved at the
cost of increased energy loss.

###  IV-C ESS model comparision

For the simulated battery, Fig. 8 shows the difference between the % change in
the state of charge of the battery estimated by three-circuit model and an
ideal lossless model for different input powers. The model difference is
particularly significant at high power levels. Due to very low $\gamma$ for
electrochemical batteries, the difference in state estimation is negligible
when comparing three-circuit models with and without considering self-
dissipation. However, for energy storage systems with high self-dissipation
rate such as flywheels, the self-dissipation phenomenon cannot be neglected.

TABLE III: Evaluation of controller against NILM algorithm with different initial battery states Initial battery SOC (%) | F-score | Energy loss (Wh) | AMBR
---|---|---|---
0 | 0.1333 | 40.421 | 152.57
25 | 0.0357 | 36.217 | 153.51
50 | 0.0357 | 36.230 | 153.51
75 | 0.2667 | 26.770 | 153.50
90 | 0.4833 | 14.174 | 153.49
100 | 0.6333 | 9.779 | 148.89
Without battery | 0.7931 | 0 | 0

Figure 8: Comparision of ESS models at 50% SOC.

##  V Conclusion

In this paper, we presented a privacy-preserving control scheme based on
Bayesian risk and a three-circuit model to estimate the energy loss associated
with a control action. The controller is modeled as a PO-MDP controlled sensor
to maximize the Bayesian risk function of an adversarial hypothesis testing
and the resulting nonlinear optimization objective is solved in a backward
recursion. Extensive numerical experiments were carried out to evaluate the
performance of the controller thoroughly. Especially, we tested the
controller’s performance against a state-of-the-art NILM algorithm using real
energy consumption data. We investigated the effect of the initial state of
the energy storage system on the controller’s performance. An important
conclusion from this work is that the privacy leakage can be reduced by using
an energy storage system but at the expense of energy loss. Without an
accurate model, the error in state estimation propagates and if not corrected,
leads to suboptimal privacy control.

Future work will focus on the trade-off between the privacy and energy loss,
time dependency of the model parameters as well as control strategy and
including more energy storage technologies.

## References

  * [1] P. McDaniel and S. McLaughlin, “Security and privacy challenges in the smart grid,” _IEEE Security & Privacy_, vol. 7, no. 3, 2009.
  * [2] G. W. Hart, “Nonintrusive appliance load monitoring,” _Proceedings of the IEEE_ , vol. 80, no. 12, pp. 1870–1891, 1992.
  * [3] A. Zoha, A. Gluhak, M. A. Imran, and S. Rajasegarar, “Non-intrusive load monitoring approaches for disaggregated energy sensing: A survey,” _Sensors_ , vol. 12, no. 12, pp. 16 838–16 866, 2012.
  * [4] M. R. Asghar, G. Dán, D. Miorandi, and I. Chlamtac, “Smart meter data privacy: A survey,” _IEEE Communications Surveys & Tutorials_, vol. 19, no. 4, pp. 2820–2835, 2017.
  * [5] G. Giaconi, D. Gunduz, and H. V. Poor, “Privacy-aware smart metering: Progress and challenges,” _arXiv preprint arXiv:1802.01166_ , 2018.
  * [6] M. LeMay, G. Gross, C. A. Gunter, and S. Garg, “Unified architecture for large-scale attested metering,” in _System Sciences, 2007. HICSS 2007. 40th Annual Hawaii International Conference on_. IEEE, 2007, pp. 115–115.
  * [7] M. Jawurek, M. Johns, and K. Rieck, “Smart metering de-pseudonymization,” in _Proceedings of the 27th Annual Computer Security Applications Conference_. ACM, 2011, pp. 227–236.
  * [8] A. Rial and G. Danezis, “Privacy-preserving smart metering,” in _Proceedings of the 10th annual ACM workshop on Privacy in the electronic society_. ACM, 2011, pp. 49–60.
  * [9] G. Kalogridis, C. Efthymiou, S. Z. Denic, T. A. Lewis, and R. Cepeda, “Privacy for smart meters: Towards undetectable appliance load signatures,” in _Smart Grid Communications (SmartGridComm), 2010 First IEEE International Conference on_. IEEE, 2010, pp. 232–237.
  * [10] D. Varodayan and A. Khisti, “Smart meter privacy using a rechargeable battery: Minimizing the rate of information leakage,” in _Acoustics, Speech and Signal Processing (ICASSP), 2011 IEEE International Conference on_. IEEE, 2011, pp. 1932–1935.
  * [11] O. Tan, D. Gunduz, and H. V. Poor, “Increasing smart meter privacy through energy harvesting and storage devices,” _IEEE Journal on Selected Areas in Communications_ , vol. 31, no. 7, pp. 1331–1341, 2013.
  * [12] M. Backes and S. Meiser, “Differentially private smart metering with battery recharging,” in _Data Privacy Management and Autonomous Spontaneous Security_. Springer, 2014, pp. 194–212.
  * [13] Z. Li, T. J. Oechtering, and M. Skoglund, “Privacy-preserving energy flow control in smart grids,” in _Acoustics, Speech and Signal Processing (ICASSP), 2016 IEEE International Conference on_. IEEE, 2016, pp. 2194–2198.
  * [14] J.-X. Chin, T. T. De Rubira, and G. Hug, “Privacy-protecting energy management unit through model-distribution predictive control,” _IEEE Transactions on Smart Grid_ , 2017.
  * [15] C.-T. Pham and D. Månsson, “On the physical system modelling of energy storages as equivalent circuits with parameter description for variable load demand (part i),” _Journal of Energy Storage_ , vol. 13, pp. 73–84, 2017\.
  * [16] M. Chen and G. A. Rincon-Mora, “Accurate electrical battery model capable of predicting runtime and iv performance,” _IEEE transactions on energy conversion_ , vol. 21, no. 2, pp. 504–511, 2006.
  * [17] P. K. Varshney, _Distributed detection and data fusion_. Springer Science & Business Media, 2012.
  * [18] V. Krishnamurthy, _Partially Observed Markov Decision Processes_. Cambridge University Press, 2016.
  * [19] C. Beckel, W. Kleiminger, R. Cicchetti, T. Staake, and S. Santini, “The eco data set and the performance of non-intrusive load monitoring algorithms,” in _Proceedings of the 1st ACM International Conference on Embedded Systems for Energy-Efficient Buildings (BuildSys 2014). Memphis, TN, USA_. ACM, Nov. 2014, pp. 80–89.
  * [20] R. H. Byrd, J. C. Gilbert, and J. Nocedal, “A trust region method based on interior point techniques for nonlinear programming,” _Mathematical Programming_ , vol. 89, no. 1, pp. 149–185, 2000.
  * [21] M. Weiss, A. Helfenstein, F. Mattern, and T. Staake, “Leveraging smart meter data to recognize home appliances,” in _Pervasive Computing and Communications (PerCom), 2012 IEEE International Conference on_. IEEE, 2012, pp. 190–197.


††thanks: A. Wikner and J. Pathak contributed equally to this work.
Correspondence should be addressed to<EMAIL_ADDRESS>A. Wikner and
J. Pathak contributed equally to this work. Current address: Lawrence Berkeley
National Laboratory (LBNL), Berkeley, CA 94720

# Using Data Assimilation to Train a Hybrid Forecast System that Combines
Machine-Learning and Knowledge-Based Components

Alexander Wikner  Department of Physics, University of Maryland, College Park,
MD 20742  Jaideep Pathak  Department of Physics, University of Maryland,
College Park, MD 20742  Brian R. Hunt  Department of Mathematics, University
of Maryland, College Park, MD 20742  Institute for Physical Science and
Technology (IPST), College Park, MD 20742  Istvan Szunyogh  Department of
Atmospheric Sciences, Texas A& M University, College Station, TX 77843
Michelle Girvan  Department of Physics, University of Maryland, College Park,
MD 20742  Institute for Physical Science and Technology (IPST), College Park,
MD 20742  Institute for Research in Electronics and Applied Physics (IREAP),
College Park, MD 20742  Edward Ott  Department of Physics, University of
Maryland, College Park, MD 20742  Institute for Research in Electronics and
Applied Physics (IREAP), College Park, MD 20742  Department of Electrical and
Computer Engineering, University of Maryland, College Park, MD 20742

###### Abstract

We consider the problem of data-assisted forecasting of chaotic dynamical
systems when the available data is in the form of noisy partial measurements
of the past and present state of the dynamical system. Recently there have
been several promising data-driven approaches to forecasting of chaotic
dynamical systems using machine learning. Particularly promising among these
are hybrid approaches that combine machine learning with a knowledge-based
model, where a machine-learning technique is used to correct the imperfections
in the knowledge-based model. Such imperfections may be due to incomplete
understanding and/or limited resolution of the physical processes in the
underlying dynamical system, e.g., the atmosphere or the ocean. Previously
proposed data-driven forecasting approaches tend to require, for training,
measurements of all the variables that are intended to be forecast. We
describe a way to relax this assumption by combining data assimilation with
machine learning. We demonstrate this technique using the Ensemble Transform
Kalman Filter (ETKF) to assimilate synthetic data for the 3-variable Lorenz
system and for the Kuramoto-Sivashinsky system, simulating model error in each
case by a misspecified parameter value. We show that by using partial
measurements of the state of the dynamical system, we can train a machine
learning model to improve predictions made by an imperfect knowledge-based
model.

> A common desired function is that of cyclic forecasting of a chaotic
> system’s state, in which a set of forecasts to various time increments in
> the future is periodically produced once every cycle time $T$ (e.g., for
> weather forecasting, $T=6$ hours is common). To do this, at the beginning of
> each cycle, available measurements taken during the last cycle are combined
> with (“assimilated into”) a forecast of the current state made at the
> beginning of the previous cycle to produce an estimate of the current system
> state. This estimate then serves as the initial condition for the prediction
> model which produces the next set of forecasts, and so on. The process of
> cyclic formation of the required initial conditions is called “Data
> Assimilation.” Recently, prediction models utilizing machine learning have
> been developed and shown to have great promise. Thus, the training and use
> of machine-learning-assisted prediction models in cyclic forecasting systems
> is of considerable interest. To enable this, it is necessary to train the
> machine-learning system, using measured data, to forecast variables that are
> not measured. Here we show this can be accomplished using data assimilation,
> and demonstrate the ability of machine learning to improve the forecast of
> variables for which no directly measured training data is available.

##  I Introduction

The application of machine learning (ML) to forecasting chaotic dynamical
systems (e.g., [Lorenz, 1963]) has been a subject of much recent interest
[Jaeger and Haas, 2004; Lukoševičius, 2012; Pathak  _et al._ , 2017, 2018a,
2018b; Bocquet  _et al._ , 2019; Chattopadhyay, Hassanzadeh, and Subramanian,
2020; Brajard  _et al._ , 2020; Bocquet  _et al._ , 2020; Wikner  _et al._ ,
2020; Arcomano  _et al._ , 2020; Bocquet, Farchi, and Malartic, 2020].
Promising results of tests promote the possibility of using ML for large scale
applications such as weather forecasting. Particularly, Ref. Pathak  _et al._
, 2018a indicates that machine-learning techniques can be scaled effectively
to forecast arbitrarily high-dimensional chaotic attractors. In cases like
weather forecasting where a skillful knowledge-based model is available,
hybrid approaches that use ML to correct model deficiencies are especially
promising [Pathak  _et al._ , 2018b; Reichstein  _et al._ , 2019]. By
“knowledge-based”, we mean derived from fundamental principles such as the
laws of physics.

So far, most data-driven forecasting approaches assume that training data is
available for all the variables to be forecast. This assumption is relaxed in
[Bocquet  _et al._ , 2019; Brajard  _et al._ , 2020; Bocquet  _et al._ , 2020;
Bocquet, Farchi, and Malartic, 2020], which describe techniques to develop a
data-driven forecast model on a specified, but potentially only partially
measured, model state space. Using an approach different from those in Refs.
Bocquet  _et al._ , 2019; Brajard  _et al._ , 2020; Bocquet  _et al._ , 2020;
Bocquet, Farchi, and Malartic, 2020, in this paper, we develop a technique for
hybrid knowledge-based/ML systems in the partially measured scenario. Both our
hybrid approach and the purely data-driven approaches in [Bocquet  _et al._ ,
2019; Brajard  _et al._ , 2020; Bocquet  _et al._ , 2020; Tomizawa and Sawada,
2020] are based on the paradigm of data assimilation (DA), which has long been
used in weather forecasting and other applications.

When a forecast model is available, DA (e.g., [Law, Stuart, and Zygalakis,
2015; Asch, Bocquet, and Nodet, 2016; Fletcher, 2017]) determines initial
conditions for an ongoing sequence of forecasts, based on an ongoing sequence
of measurements. DA does not require that current measurements are sufficient
to determine the current model state; instead, it typically uses the forecast
model to propagate information from past measurements to augment current
information. This technique helps both to constrain the model state and to
filter noise in the measurements.

Ideas from ML and DA have been combined in various ways in the literature. DA
has been used with data-driven forecast models based on delay-coordinate
embedding [Hamilton, Berry, and Sauer, 2016] and analog forecasting [Lguensat
_et al._ , 2017]. On the other hand, with appropriate training data, machine-
learning techniques can be taught to perform DA, without explicitly
formulating a forecast model [Lu  _et al._ , 2017]. DA has also been employed
as a paradigm and a tool for the training of machine-learning systems
[Puskorius and Feldkamp, 1994; Abarbanel, Rozdeba, and Shirman, 2018; Bocquet
_et al._ , 2019; Brajard  _et al._ , 2020; Bocquet  _et al._ , 2020; Bocquet,
Farchi, and Malartic, 2020].

We give a brief introduction to DA in Sec. II.1. In Sec. II.2, we formulate a
potentially effective scheme combining DA (Sec. II.1) with machine-learning-
assisted prediction (as described in Ref. Pathak  _et al._ , 2018b) to create
a cyclic forecast system. In Sec. II.3, we propose an iterative procedure to
improve the DA cycle and potentially produce improved ML-assisted forecast
models. In Sec. III, we show the results of numerical experiments.

##  II Methods

We consider the problem of forecasting a chaotic, nonlinear dynamical system
represented by the equation

$\displaystyle\frac{{d\tilde{\textbf{x}}}(t)}{dt}$
$\displaystyle=\textbf{F}[\tilde{\textbf{x}}(t)].$ (1)

We assume that we are only able to obtain partial measurements of the state
$\tilde{\textbf{x}}$ at regular time intervals. We denote the measurements
made at time $k\Delta t$ by $\textbf{y}_{k}$, so that

$\displaystyle\textbf{y}_{k}$
$\displaystyle=\textbf{H}\tilde{\textbf{x}}(k\Delta t)+\eta_{k}.$ (2)

In Eq. (2), H is often called the observation operator; here, we call it the
measurement operator. For our exposition, we assume H to be linear, as will be
true in the numerical experiments performed in this paper. (Many DA methods
allow a nonlinear measurement operator, and our ML method is compatible with
this case.) The random variable $\eta_{k}$ represents measurement noise, which
is typically assumed to be Gaussian with mean 0 and covariance matrix R.

In applications, we do not have full knowledge of the dynamics of the system,
and, as such, Eq. (1) is unknown. What often is available is an imperfect
model of the dynamics, with model error due to imperfect representations of
some of the physical processes that govern the system. We assume that we have
access to such an imperfect “knowledge-based” model denoted by G and that the
dynamical equation,

$\displaystyle\frac{{d\textbf{x}}(t)}{dt}$
$\displaystyle=\textbf{G}[\textbf{x}(t)],$ (3)

can, via evolution of x, be used to forecast an approximation of the future
$\tilde{\textbf{x}}$-dependent state properties of interest. For our
exposition below, we assume that x and $\tilde{\textbf{x}}$ lie in the same
space, but in an application such as weather forecasting, x can represent the
projection of the physical state $\tilde{\textbf{x}}$ onto a model grid.

###  II.1 Data Assimilation

In DA, the initial state for Eq. (3), $\textbf{x}(0)$, is estimated from the
current measurement $\textbf{y}_{0}$ and all past measurements
$\\{\textbf{y}_{k}\\}_{-T\leq k<0}$, using the model, Eq. (3), to imperfectly
propagate information from the past to the present. Many DA methods are based
on the following nonlinear least-squares problem. Given a dynamical model, Eq.
(3), and a set of measurements $\\{\textbf{y}_{k}\\}_{-T\leq k\leq 0}$, find a
trajectory $\textbf{x}(t)$ of the model G that minimizes the cost function
given by

$\displaystyle
J_{0}(\textbf{x}(t))=\sum_{k=-T}^{0}\rho^{k}\left[\textbf{y}_{k}-\textbf{H}\textbf{x}(k\Delta
t)\right]^{\intercal}\textbf{R}^{-1}\left[\textbf{y}_{k}-\textbf{H}\textbf{x}(k\Delta
t)\right].$ (4)

Here, $\rho\geq 1$ is an adjustable constant, referred to as the covariance
inflation parameter, that discounts past measurements, in part to compensate
for model error in G. In the situation where the measurement errors
($\eta_{k}$) are Gaussian and the model G is a perfect representation of the
dynamical system, i.e., $\textbf{G}=\textbf{F}$, minimizing Eq. (4) with
$\rho=1$ will give us the trajectory of the model $\textbf{x}(t)$ that is the
most likely trajectory of the dynamical system in the Bayesian sense (here, we
assume a “flat” prior likelihood; the choice of prior distribution becomes
increasingly unimportant as $T$ increases.)

In forecast applications, the purpose of DA is to determine an estimate
$\textbf{x}(0)$ of the current system state, which is used as an initial
condition for Eq. (3) to forecast $\textbf{x}(t)$ for $t>0$. As such, the
trajectory $\textbf{x}(t)$ for $t<0$ that minimizes the cost function $J_{0}$
is often not computed explicitly; only its final state $\textbf{x}(0)$ is
needed as the algorithmic output. DA is typically done sequentially, using the
state estimate from an earlier time to compute the estimate for the current
time. To be more precise, consider the cost function $J_{j}$ that generalizes
Eq. (4) as follows:

$\displaystyle
J_{j}(\textbf{x}(t))=\sum_{k=-T}^{j}\rho^{k-j}\left[\textbf{y}_{k}-\textbf{H}\textbf{x}(k\Delta
t)\right]^{\intercal}\textbf{R}^{-1}\left[\textbf{y}_{k}-\textbf{H}\textbf{x}(k\Delta
t)\right].$ (5)

Then,

$\displaystyle
J_{j}(\textbf{x}(t))=\rho^{-1}J_{j-1}(\textbf{x}(t))+\left[\textbf{y}_{j}-\textbf{H}\textbf{x}(j\Delta
t)\right]^{\intercal}\textbf{R}^{-1}\left[\textbf{y}_{j}-\textbf{H}\textbf{x}(j\Delta
t)\right].$ (6)

Suppose for some $-T<j<=0$ that $\textbf{x}^{a}_{j-1}$ represents an estimate
of $\tilde{\textbf{x}}((j-1)\Delta t)$ based on the cost function $J_{j-1}$;
i.e., based on measurements $\textbf{y}_{k}$ up to and including $k=j-1$. The
next estimate $\textbf{x}_{j}^{a}$ is computed in two steps as follows:

  1. 1.

Forecast: The model, Eq. (3), is applied from time $(j-1)\Delta t$ to $j\Delta
t$, using $\textbf{x}^{a}_{j-1}$ as the initial condition, to obtain an
estimate $\textbf{x}^{b}_{j}$ to $\tilde{\textbf{x}}(j\Delta t)$, called the
“background” state at time $j\Delta t$.

  2. 2.

Analysis: The background $\textbf{x}^{b}_{j}$ is used along with the current
measurement vector $\textbf{y}_{j}$ to compute the “analysis” state
$\textbf{x}^{a}_{j}$ that estimates the minimizer of $J_{j}$ at time $j\Delta
t$.

These two steps are repeated to compute $\textbf{x}_{j+1}^{b}$,
$\textbf{x}_{j+1}^{a}$, and so on.

Different DA methods perform the analysis step differently. In some methods,
including those based on the Kalman filter, the cost functions are
approximated by quadratic functions of $\textbf{x}(j\Delta t)$, so that the
right side of Eq. (6) becomes

$\displaystyle\rho^{-1}\left[\mathbf{x}(j\Delta
t)-\mathbf{x}_{j}^{b}\right]^{\intercal}\left(\mathbf{P}_{j}^{b}\right)^{-1}\left[\mathbf{x}(j\Delta
t)-\mathbf{x}_{j}^{b}\right]+\left[\textbf{y}_{j}-\textbf{H}\textbf{x}(j\Delta
t)\right]^{\intercal}\textbf{R}^{-1}\left[\textbf{y}_{j}-\textbf{H}\textbf{x}(j\Delta
t)\right].$ (7)

Thus, the quantity to be minimized is a sum of (squared) weighted distances
from $\textbf{x}(j\Delta t)$ to the background state (which is forecasted from
previous observations) and to the current observations. Here
$\textbf{P}^{b}_{j}$ is a covariance matrix associated with the background
state at time $j\Delta t$. In an ensemble Kalman filter [Evensen, 2003;
Houtekamer and Mitchell, 1998; Burgers, Jan van Leeuwen, and Evensen, 1998],
$\textbf{P}^{b}_{j}$ is computed as the sample covariance of an ensemble of
forecasts. The parameter $\rho$ effectively inflates this covariance,
compensating for (in addition to model error) sampling errors due to finite
ensemble size and the effect of nonlinearity on the error dynamics [Bishop,
Etherton, and Majumdar, 2001; Whitaker and Hamill, 2002].

In the Appendix, we describe the particular DA method we use, the “Ensemble
Transform Kalman Filter” (ETKF) [Bishop, Etherton, and Majumdar, 2001; Wang,
Bishop, and Julier, 2004], a type of ensemble Kalman filter. We emphasize that
the approach we describe in Sec. II.2 can be used with any DA method. The
forecast and analysis steps we have described are applied iteratively to
produce a time series of analysis states $\textbf{x}^{a}_{j}$. We will use
this time series in Sec. II.2 to train the ML component of our method.

###  II.2 ML-DA Algorithm

In what follows, we use a reservoir computer implementation of our ML-DA
algorithm based on Refs. Jaeger and Haas, 2004; Lukoševičius, 2012, which is
similar to that of Ref. Pathak  _et al._ , 2018b. The main component of a
reservoir computer is the reservoir network. The reservoir network is a
sparse, randomly generated graph with weighted, directed edges and dynamical
nodes whose states evolve in response to input signals transported along
network edges. The network topology can be represented completely by a
weighted adjacency matrix, denoted A. We choose A to be a $D_{r}\times D_{r}$
sparse, randomly generated matrix. The network is constructed to have an
average in-degree (number of incoming edges per node) denoted by $\langle
d\rangle$, and the nonzero elements of A, representing the edge weights in the
network, are initially chosen independently from the uniform distribution over
the interval $[0,1]$. All the edge weights in the network are then uniformly
scaled via multiplication of the adjacency matrix by a constant factor to set
the largest magnitude eigenvalue of the matrix to a desired value $\omega$,
which is called the “spectral radius” of A. The state of the reservoir, given
by the vector $\textbf{r}(t)$, consists of the components $r_{n}$ for $1\leq
n\leq D_{r}$ where $r_{n}(t)$ denotes the scalar state of the $n^{th}$ node in
the network. The reservoir is coupled to the $M$ dimensional input through a
$D_{r}\times M$ dimensional matrix $\textbf{W}_{in}$. Each row of the matrix
$\textbf{W}_{in}$ has exactly one nonzero element which is independently and
randomly chosen from the uniform distribution on the interval
$[-\zeta,\zeta]$. We choose the locations of these nonzero elements such that
each node receives input from only one input variable and each input variable
is input to an equal number of nodes.

As outlined at the beginning of Sec. II, we assume that the finite-dimensional
representation of the true dynamical system is given by Eq. (1). We have
$T+T_{s}$ measurements given by $\\{\textbf{y}_{j}\\}$, in the interval
$-T-T_{s}\leq j\leq 0$, which are related to the true state of the dynamical
system by Eq. (2). We refer to $T$ and $T_{s}$ as the training and
synchronization times, where $T\gg T_{s}$. We further assume that the model G
(Eq. (3)) is imperfect. Using the model G and the DA algorithm outlined in
Sec. Appendix, we can obtain an analysis state $\textbf{x}^{a}_{j}$ at each
time step $j$, $-T-T_{s}\leq j\leq 0$. We thus create a set of “analysis”
states $\\{\textbf{x}_{j}^{a}\\}_{-T-T_{s}\leq j\leq 0}$. We are interested in
forecasting the state of the dynamical system for $j>0$. In the standard
forecast approach (without ML), one would predict the future state of the
dynamical system (Eq. (1)) using $\textbf{x}^{a}_{0}$ as the initial condition
and Eq. (3) as the dynamical model. We will call this forecast the baseline
forecast, to which we will compare the forecast made by our ML-DA algorithm.
The ML-DA algorithm is a natural extension of the algorithm proposed in Ref.
Pathak  _et al._ , 2018b for hybrid forecasting of chaotic dynamical systems
by using ML in conjunction with a knowledge-based model. Figure 1 shows a
schematic illustration of the ML-DA scheme described below.

Figure 1: (a) Training the ML-DA prediction scheme and (b) forecasting as
described in Sec. II.2. The output matrix in (a), $\mathbf{W}_{out}$, is
unknown at the time of training. The goal of training is to determine
$\mathbf{W}_{out}$ such that $\mathbf{x}^{H}_{j}$ best approximates
$\mathbf{x}^{a}_{j}$ in the regularized least-squares sense.

####  II.2.1 Training

  1. 1.

Use the model G to create a set of forecasts from the latest $T$ analysis
states $\textbf{x}^{a}_{j}$, $-T\leq j\leq 0$. We denote these time-$\Delta t$
forecasts by $\textbf{x}^{M}_{j}$, $-T+1\leq j\leq 1$.

  2. 2.

Initialize the reservoir state $\textbf{r}_{-T-T_{s}}$ to a random value.

  3. 3.

Evolve the reservoir computer using the following equation

$\displaystyle\textbf{r}_{j+1}$
$\displaystyle=\tanh[\textbf{A}\textbf{r}_{j}+\textbf{W}_{in}\textbf{x}^{a}_{j}]$
(8)

for $-T-T_{s}\leq j\leq 0$. Here, $\tanh$ is applied element-by-element to its
input vector.

  4. 4.

Find a set of output weights $\textbf{W}_{out}$ so that

$\displaystyle\mathbf{x}^{H}_{j}=\textbf{W}_{out}\left[\begin{matrix}\textbf{r}_{j}\\\
\textbf{x}^{M}_{j}\end{matrix}\right]\simeq\textbf{x}^{a}_{j}$ (9)

for $-T+1\leq j\leq 0$. The matrix $\textbf{W}_{out}$ is computed using
regularized least squares regression. Thus, we find the $\textbf{W}_{out}$
that minimizes the following cost function:

$\displaystyle\ell(\textbf{W}_{out})=\sum_{j=-T+1}^{0}\lVert\textbf{W}_{out}\left[\begin{matrix}\textbf{r}_{j}\\\
\textbf{x}^{M}_{j}\end{matrix}\right]-\textbf{x}^{a}_{j}\rVert^{2}+\beta\lVert\textbf{W}_{out}\rVert^{2}.$
(10)

Notice that we have not used the first $T_{s}$ reservoir states to find the
optimal $\textbf{W}_{out}$. This is to ensure that the reservoir states used
in the computation of $\textbf{W}_{out}$ do not depend strongly on the initial
random value of the reservoir state, $\mathbf{r}_{-T-T_{s}}$, and thus are
“synchronized” to the measurements that drive Eq. 8 through the analysis
states $\textbf{x}^{a}_{j}$.

####  II.2.2 Prediction

We now describe the hybrid forecast model that uses the trained reservoir
along with the knowledge-based model G to forecast the state
$\textbf{x}(j\Delta t)$ for $j>0$.

The hybrid reservoir runs in a closed-loop feedback mode according to the
following steps:

  1. 1.

Obtain the hybrid ML prediction at time $j\Delta t$ (denoted
$\textbf{x}^{H}_{j}$) according to:

$\displaystyle\textbf{x}^{H}_{j}=\textbf{W}_{out}\left[\begin{matrix}\textbf{r}_{j}\\\
\textbf{x}^{M}_{j}\end{matrix}\right].$ (11)

The reservoir state and knowledge-based forecasts at the first step,
$\mathbf{r}_{1}$ and $\mathbf{x}^{M}_{1}$, are computed from the final
training analysis $\textbf{x}^{a}_{0}$ according to steps 1 and 3 of the
training procedure.

  2. 2.

Compute the imperfect model forecast at time $(j+1)\Delta t$,
$\textbf{x}^{M}_{j+1}$, by evolving $\textbf{x}^{H}_{j}$ using the imperfect
model G.

  3. 3.

Evolve the reservoir according to

$\displaystyle\textbf{r}_{j+1}$
$\displaystyle=\tanh\left[\textbf{A}\textbf{r}_{j}+\textbf{W}_{in}\textbf{x}_{j}^{H}\right]$
(12)

  4. 4.

Increment $j$ and repeat steps 1 to 3.

If one wants to initialize a forecast at a time $J\Delta t$ for $J>0$, one
should continue to run the DA procedure for $J$ more steps after the training
time period, set $\textbf{x}^{H}_{J}=\textbf{x}^{a}_{J}$, then proceed as
above starting with step 2 for $j=J$. The $J$ DA cycles after time 0 can be
performed using the trained hybrid model in place of the knowledge-based
model, similarly to the approach described in the next section.

We note that the ensemble-based DA method we use could be used to initialize
an ensemble forecast. However, for simplicity, our results use only forecasts
from the ensemble mean.

###  II.3 Iterated ML-DA Algorithm

Once we have trained a hybrid model using the process described in Sec.
II.2.1, we can iterate the training procedure in the following sense. An
iterative approach, initialized differently, is described in Ref. Brajard  _et
al._ , 2020. As shown in Fig. 2, we replace the knowledge-based model used
during DA with the previously trained hybrid model. At iteration $i=1$, the
hybrid is trained using the analysis states produced using only the knowledge-
based model. We denote the output matrix to be determined by the training, the
state of the reservoir at time $j\Delta t$, the analysis state at time
$j\Delta t$, and the hybrid model predictions for time $j\Delta t$ at
iteration $i$ as $\mathbf{W}_{out,i}$, $\mathbf{r}_{j,i}$,
$\mathbf{x}^{a}_{j,i}$, and $\mathbf{x}^{H}_{j,i}$, respectively.

Training proceeds as in Sec. II.2.1 with the following modifications.

Figure 2: Training subsequent iterations ($i\geq 2$) of the ML-DA algorithm
using the synchronized hybrid model during DA as described in Sec. II.3. The
index $i$ denotes the current iteration of the algorithm, and we consider the
training procedure shown in Fig. 1(a) to occur at $i=1$. At iteration $i$,
$\mathbf{W}_{out,i-1}$ is known from the previous training iteration, while we
endeavor to determine $\mathbf{W}_{out,i}$ such that $\mathbf{x}^{H}_{j,i}$
best approximates $\mathbf{x}^{a}_{j,i}$ in the regularized least-squares
sense.

####  II.3.1 Training

  1. 1.

For $-T-T_{s}\leq j<-T$, we use the knowledge-based model to perform DA,
producing analysis states $\mathbf{x}^{a}_{j,i}$. We do not use the hybrid
model during this synchronization time period because the reservoir states
$\textbf{r}_{j,i}$ might not yet be synchronized with the measurements. These
reservoir states are generated according to step 3 of Sec. II.2.1, i.e.,

$\displaystyle\mathbf{r}_{j+1,i}=\tanh[\mathbf{A}\mathbf{r}_{j,i}+\mathbf{W}_{in}\mathbf{x}^{a}_{j,i}].$
(13)

  2. 2.

For $T\leq j\leq 0$, we use the hybrid model to perform DA in place of the
knowledge-based model, producing analysis states $\mathbf{x}^{a}_{j,i}$, as
depicted in Fig. 2. We continue to use Eq. (13) to evolve the reservoir
states.

  3. 3.

Using these new analysis states, we follow the steps contained in Sec. II.2.1
to compute a new hybrid output matrix, $\mathbf{W}_{out,i+1}$.

  4. 4.

Once the new hybrid output matrix has been computed, predictions can be
produced by following Sec. II.2.2.

In the case of an ensemble-based assimilation method (such as the ETKF
described in Appendix Appendix), we augment each ensemble member to include,
in addition to the system state $\mathbf{x}^{a,k}_{j}$, a corresponding
reservoir state $\mathbf{r}^{k}_{j}$. This reservoir state is synchronized to
the corresponding ensemble system state by evolving it as in Eq. (13):

$\displaystyle\mathbf{r}^{k}_{j+1,i}=\tanh[\mathbf{A}\mathbf{r}^{k}_{j,i}+\mathbf{W}_{in}\mathbf{x}^{a,k}_{j,i}].$
(14)

Then, $\textbf{r}^{k}_{j+1,i}$ is used in the hybrid model to produce the
$k^{\textrm{th}}$ background ensemble member for the next analysis cycle. We
note that the reservoir state(s) at $j=-T+1$ needed to begin DA using the
hybrid model must only be computed at the first iteration. Afterwards, these
initial synchronized states may be reused at the beginning of subsequent DA
cycles.

##  III Results

We demonstrate the forecast performance of the ML-DA algorithm in comparison
with the performance of baseline for forecasts prepared without using MD. We
consider two dynamical systems as our examples: the Lorenz ‘63 system [Lorenz,
1963] and the Kuramoto-Sivashinsky (KS) system [Kuramoto and Tsuzuki, 1976;
Sivashinsky, 1977]. For the latter system, we additionally demonstrate the
decrease in error in the analysis states from using the iterative method
described in Sec. II.3.

###  III.1 Experimental Design

For the purpose of evaluating forecast performance, simulated true states
$\hat{\textbf{x}}_{j}$ are generated by the perfect model dynamics in the time
interval $-T-T_{S}\leq j\leq P$, where $T$ and $T_{S}$ are described in Sec.
II.2. The data in the interval $1\leq j\leq P$ is set aside to verify the
baseline and ML-DA forecasts but is not otherwise used, since it is considered
to be in the future. We assume that we only have access to the measurements
$\\{\textbf{y}_{j}\\}_{-T-T_{S}\leq j\leq 0}$ generated from the simulated
true states via Eq. (2) and do not know the perfect model, Eq. (1). We further
assume that we have access to an imperfect model G (Eq. 3). We use this
imperfect model to perform DA, to train the ML-DA, and to prepare both the
baseline and ML-DA forecasts for $j>0$.

For DA, we use the ETKF, as described in the Appendix, with 15 ensemble
members for the Lorenz ’63 system and 30 ensemble members for the KS system.
To initialize the forecasts and the ML training, we perform DA using the
measurements $\textbf{y}_{j}$ ($-T-T_{S}\leq j\leq 0$) and the model G to
obtain the analysis states
$\textbf{x}^{a}_{-T-T_{s}},\dots,\textbf{x}^{a}_{0}$.

####  III.1.1 Baseline Forecasts:

Using $\textbf{x}^{a}_{0}$ as the initial condition, we compute a forecast
using the model G from time $j=1$ to $j=P$. We call this forecast
$\textbf{x}^{f,base}_{j}$, $1\leq j\leq P$.

####  III.1.2 ML-DA Forecasts:

We follow the ML-DA scheme detailed in Sec. II.2 using a reservoir computer
with the parameters listed in Table 1 for the Lorenz ‘63 and the Kuramoto-
Sivashinsky (KS) dynamical systems. We use the measurements $\textbf{y}_{j}$
($-T-T_{S}\leq j\leq 0$) and the imperfect model G corresponding to the Lorenz
and KS systems respectively to train the reservoir computer. The exact form of
the measurements and of the imperfect model is described in Secs. III.2 and
III.3. We then forecast using the ML-DA scheme from $j=1$ to $j=P$. We call
this forecast $\textbf{x}^{f,ml}_{j}$, $1\leq j\leq P$.

Hyperparameter | Lorenz 63 | KS
---|---|---
$D_{r}$ | $1000$ | $2000$
$\langle d\rangle$ | 3 | 3
$\omega$ | 0.9 | 0.6
$\zeta$ | 0.1 | 1.0
Table 1: Reservoir Hyperparamenters

####  III.1.3 Normalized RMS error

The normalized RMS error in the baseline forecast ($e_{j}^{base}$) and the RMS
error in ML-DA forecast ($e^{ml}_{j}$) are calculated as follows:

$\displaystyle e_{j}$
$\displaystyle=\frac{\lVert\textbf{x}^{f}_{j}-\hat{\textbf{x}}_{j}\rVert}{\sqrt{\langle\lVert\hat{\textbf{x}}_{j}\rVert^{2}\rangle}},\hskip
36.135pt1\leq j\leq P,$ (15)

where $||\dots||$ denotes the Euclidean norm, $\langle\dots\rangle$ denotes an
average over the prediction interval ($1\leq j\leq P$), and $e_{j}$,
$\textbf{x}^{f}_{j}$ denote either $e^{base}_{j}$, $\textbf{x}^{f,base}_{j}$
or $e^{ml}_{j}$, $\textbf{x}^{f,ml}_{j}$. The superscript $f$ in
$\textbf{x}^{f}_{j}$ stands for forecast.

####  III.1.4 Valid Time

We define the Valid Time (VT) as the time $j\Delta t$ at which the RMS error
($e_{j}$) first exceeds a threshold $\kappa$ chosen to be $0.9$. In order to
report the forecast quality on a meaningful natural time scale, we normalize
the valid time by the ‘Lyapunov time’ of the dynamical system. The Lyapunov
time is defined as the inverse of the Largest Lyapunov Exponent (LLE), denoted
$\Lambda_{max}$, of the dynamical system being forecast. In a chaotic system,
infinitesimal errors grow exponentially on average as $\exp(\Lambda_{max}t)$.
Thus, $\Lambda_{max}^{-1}$ is a natural time scale for evaluating forecast
quality in a chaotic system.

Since forecast Valid Time can vary considerably depending on what part of the
attractor forecasts are initialized from, we present Valid Time results as box
plots reflecting 100 different runs for the Lorenz System and 50 runs for the
KS system. Each run uses a different initial condition for the training data
(and therefore for the following forecast) and a different random realization
of the reservoir. All box plots depict the median, $25^{\textrm{th}}$ and
$75^{\textrm{th}}$ percentiles, and the $5^{\textrm{th}}$ and
$95^{\textrm{th}}$ percentiles (depicted by dashed lines).

####  III.1.5 Analysis RMS Error

For the iterative method in Sec. II.3, we evaluate the error in the analysis
state at a particular iteration $i$ by computing the normalized RMS error over
the entire analysis period,

$\displaystyle
e^{a}_{tot,i}=\sqrt{\frac{\langle\lVert\mathbf{x}^{a}_{j,i}-\mathbf{x}_{j}\rVert^{2}\rangle}{\langle\lVert\mathbf{x}_{j}\rVert^{2}\rangle}},$
(16)

where $\langle\dots\rangle$ denotes an average over the DA interval
($-T-T_{S}\leq j\Delta t\leq 0$). In addition, we evaluate the analysis error
in the measured or unmeasured system variables by including only the
particular set of variables in the norms in Eq. (16).

###  III.2 Lorenz 63

The Lorenz system is described by the equations,

$\displaystyle\frac{dX_{1}}{dt}$ $\displaystyle=-aX_{1}+aX_{2}$ (17)
$\displaystyle\frac{dX_{2}}{dt}$ $\displaystyle=bX_{1}-X_{2}-X_{1}X_{3}$ (18)
$\displaystyle\frac{dX_{3}}{dt}$ $\displaystyle=-cX_{3}+X_{1}X_{2}$ (19)

where $a=10$, $b=28$, and $c=8/3$. In this example, the true dynamical system,
Eq. (1), is given by Eqs. (17-19). We obtain simulated data by integrating
Eqs. (17-19) using a fourth order Runge-Kutta solver with time step
$\tau=0.01$. We sample the time series at intervals $\Delta t=0.01$. This data
is taken to represent the true system state $\hat{\textbf{x}}(j\Delta t)$,
$-T-T_{S}\leq j\leq P$. We use the data in the interval $-T-T_{S}\leq j\leq 1$
to create simulated measurements

$\displaystyle\textbf{y}_{j}=\textbf{H}\hat{\textbf{x}}_{j}+\boldsymbol{\eta}_{j}.$
(20)

We will study two examples with two different forms of the measurement
operator H, one in which only the $X_{1}$ variable is measured,

$\displaystyle\textbf{H}=\begin{bmatrix}1&0&0\end{bmatrix},$ (21)

and another in which the $X_{1}$ and $X_{3}$ variables are measured,

$\displaystyle\textbf{H}=\begin{bmatrix}1&0&0\\\ 0&0&1\end{bmatrix}.$ (22)

Measurement noise $\boldsymbol{\eta}_{j}$ with a Gaussian distribution that
has mean zero and covariance matrix R is added to the measurements. The
covariance matrix is diagonal, with the diagonal elements equal to
$\sigma_{noise}^{2}$, which corresponds to the situation in which the errors
in the measurements of the different state vector components are statistically
independent. In Sec. III.2.1 and Sec. III.2.2, we let $\sigma_{noise}=0.1$.
For comparison, the standard deviations
$\sigma_{X_{1}},\sigma_{X_{2}},\sigma_{X_{3}}$ of the variables
$X_{1},X_{2},X_{3}$ are $7.9$, $8.9$ and $8.6$, respectively. In Sec. III.2.3,
we consider the effect of varying $\sigma_{noise}$.

We let the imperfect model equations be given by

$\displaystyle\frac{dX_{1}}{dt}$ $\displaystyle=-aX_{1}+aX_{2},$ (23)
$\displaystyle\frac{dX_{2}}{dt}$
$\displaystyle=b(1+\epsilon)X_{1}-X_{2}-X_{1}X_{3},$ (24)
$\displaystyle\frac{dX_{3}}{dt}$ $\displaystyle=-cX_{3}+X_{1}X_{2}.$ (25)

Thus, the imperfect model G differs from the perfect model F in the parameter
$b$ by a multiplicative factor $(1+\epsilon)$. Notice also that the model
error is in the evaluation of $X_{2}$, which is unmeasured in our experiments.

####  III.2.1 Optimizing the Covariance Inflation

The quality of both the baseline and ML-DA forecasts depends on the covariance
inflation parameter $\rho$ (see Eqs. (5) and (31).)This parameter is used in
the DA cycle that initializes both forecasts, as well as forming the training
data for the ML-DA. It is thus crucial that the covariance inflation parameter
is optimized independently for both the forecast schemes. In Fig. 3, for
$\textbf{H}=\begin{bmatrix}1&0&0\end{bmatrix}$, we demonstrate the dependence
of the forecast Valid Time on the covariance inflation parameter $\rho$ when
the model error is $\epsilon=0.1$. Each Valid Time is shown for a set of 100
forecasts (see Sec. III.1.4) made with the baseline scheme (magenta box plots)
and the ML-DA scheme (light blue box plots). Figure 3 shows that the ML-DA
scheme dramatically improves forecast valid time when the model has
substantial error. In particular, the best median VT for the ML-DA (which
occurs at $\rho=1.2$) exceeds the best median VT baseline (which occurs at
$\rho=1.05$) by about a factor of 3. We speculate that the optimal value of
$\rho$ is larger for ML-DA than the baseline for the following reason. Larger
values of $\rho$ cause the measurements to be weighted more heavily relative
to the background forecast from the imperfect model in the DA procedure. This
can be advantageous for the training, which seeks to use the measurements to
correct the model error.

Figure 3: Dependence of the Forecast Valid Time (VT) for the Lorenz ’63 model
on the covariance inflation factor ($\rho$) when
$\textbf{H}=\begin{bmatrix}1&0&0\end{bmatrix}$, measurement noise has
magnitude $\sigma_{noise}=0.1$, and the model error is $\epsilon=0.1$. The ML-
DA forecasts (light blue box plots) outperform the baseline forecasts (magenta
box plots) significantly.

In the results that follow, for each choice of method, model error, and noise
magnitude, we choose the covariance inflation $\rho$ to be the value that
yields the highest median Valid Time, among the values shown in Fig. 3. For
our usual values $\epsilon=0.1$ and $\rho_{noise}=0.1$, Figure 4 shows a
typical example comparing the true and predicted time series; notice that the
ML-DA scheme improves the forecast for all three variables, including those
those that are unobserved.

####  III.2.2 Dependence on Model Error

Figure 5 demonstrates the effectiveness of the ML-DA algorithm at different
values of the model error ($\epsilon$). As expected, we find that the ML-DA
optimized over $\rho$ yields the most improvement relative to the baseline
optimized over $\rho$ when the model error is large. In the cast when the
model error is zero, the ML-DA algorithm degrades the forecasts somewhat;
since there is no model error to correct, the modification made to the model
by the ML is spurious.

Figure 4: Forecast time series plots for the Lorenz ’63 model using the ML-DA
and baseline scheme for $\textbf{H}=\begin{bmatrix}1&0&0\end{bmatrix}$ and
model error $\epsilon=0.1$. Figure 5: Dependence of the Forecast Valid Time
for the Lorenz ’63 model on the model error ($\epsilon$). The box plots
represent the result of 100 independent trials comparing the baseline and the
ML-DA forecasts, each using a different Lorenz ’63 dataset starting from a
random initial condition. The covariance inflation parameter is independently
optimized for both sets of forecasts at each value of the model error
($\epsilon$). (a) and (b) differ in the measurement operator used, which is
indicated in the figures.

####  III.2.3 Effects of measurement Noise

Figure 6 shows the effects of varying amounts of measurement noise on the
baseline (magenta box plots) and ML-DA (light blue box plots) forecast valid
time for the Lorenz system. The measurement operator is
$\textbf{H}=\begin{bmatrix}1&0&0\end{bmatrix}$ and simulated uncorrelated
Gaussian noise with standard deviation $\sigma_{noise}$ is added to the
measurements. For comparison, the standard deviation of the $X_{1}$ variable
is $\sigma_{X_{1}}\simeq 7.9$. We see that the ML-DA outperforms the baseline
forecast by a larger amount when the noise is smaller. We note that smaller
amounts of noise allow the measurements to be weighted more relative to the
background forecast in the DA procedure, yielding training data more suitable
to correcting forecast errors.

Figure 6: Effect of varying amounts of measurement noise on the baseline
(magenta box plots) and ML-DA (light blue box plots) forecast valid time for
the Lorenz ’63 system, with model error $\epsilon=0.1$. The measurement
operator is $\textbf{H}=\left[\begin{matrix}1&0&0\end{matrix}\right]$, and
simulated uncorrelated Gaussian noise with standard deviation $\sigma_{noise}$
is added to the measurements. For comparison, the standard deviation of the
$X_{1}$ variable is $\sigma_{X_{1}}\simeq 7.9$.

###  III.3 Kuramoto-Sivashinsky (KS) system

In this section we consider an example of a spatiotemporal chaotic dynamical
system called the the Kuramoto-Sivashinsky (KS) system defined by the Partial
Differential Equation (PDE),

$\displaystyle\frac{\partial u(x,t)}{\partial t}+u\frac{\partial
u(x,t)}{\partial x}+\frac{\partial^{2}u(x,t)}{\partial
x^{2}}+\frac{\partial^{4}u(x,t)}{\partial x^{4}}$ $\displaystyle=0.$ (26)

Equation (26) defines the evolution of a one-dimensional spatiotemporal scalar
field $u(x,t)$ on the spatial domain $x\in[0,L)$. We assume periodic boundary
conditions so that $u(x+L,t)=u(x,t)$. In this example, we let Eq. (26)
represent the true, unknown space and time continuous dynamics. We obtain a
time series of finitite-dimensional representations of the state that we
assume to satisfy Eq. (1) by discretizing the domain $[0,L)$ into $Q$ grid-
points and integrating Eq. (26) using a pseudo-spectral PDE solver [Kassam and
Trefethen, 2005] with time step $\tau=0.25$. We sample the time series at
intervals of $\Delta t=0.25$. Thus, our simulated true system states takes the
form of a $Q$-dimensional vector $\hat{\textbf{x}}(j\Delta t)$, $-T-T_{S}\leq
j\leq P$. As in the previous example of the Lorenz ‘63 dynamical system (Sec.
III.2), we generate simulated measurements by adding random Gaussian noise to
the simulated true states. We use the simulated measurements in the interval
$-T-T_{S}\leq j\leq 0$ as the training data for the ML-DA scheme and set aside
the data in the interval $1\leq j\leq P$ for forecast verification.

We assume that we have access to an imperfect model, G, of the KS dynamical
system given by the equations

$\displaystyle\frac{\partial u(x,t)}{\partial t}+u\frac{\partial
u(x,t)}{\partial x}+(1+\epsilon)\frac{\partial^{2}u(x,t)}{\partial
x^{2}}+\frac{\partial^{4}u(x,t)}{\partial x^{4}}$ $\displaystyle=0.$ (27)

Thus, if $\epsilon=0$, then the model represents the true dynamics perfectly,
while a nonzero value of $\epsilon$ indicates an imperfect model. We also
assume that our measurements are of the form given by Eq. (2). The measurement
operator H is a projection operator that outputs $\Theta$ uniformly spaced
measurements (of the $Q$ total variables) in $\hat{\textbf{x}}_{j}$. The
measurement noise is Gaussian with zero mean and a diagonal covariance matrix
with diagonal elements equal to $\sigma_{noise}^{2}$, where
$\sigma_{noise}=0.1$. For comparison, the standard deviation of the KS time
series variable is $\sigma_{KS}=1.3$.

####  III.3.1 Forecast Quality and Dependence on Model Error

To demonstrate the effectiveness of the ML-DA technique, we consider a KS
system with $L=35$, $Q=64$ and a measurement operator H that corresponds to
$\Theta=16$ uniformly spaced measurements. The upper three panels of Fig. 7
shows (from top to bottom) the time dependence of the true discretized
Kuramoto-Sivashinksy equation dynamics, the error in the baseline prediction,
and the error in the ML-DA prediction when the model error is $\epsilon=0.1$.
The bottom panel in Fig. 7 shows the normalized RMS error of both predictions.
Notice that the ML-DA again significantly outperforms the baseline scheme.

Figure 8 shows the forecast Valid Time for $50$ different forecasts per box
plot made with the ML-DA and baseline scheme at four different values of the
model error $\epsilon$. As for the Lorenz system, the ML-DA scheme outperforms
the baseline scheme when the error in the imperfect model is high, but cannot
improve on the baseline when the model error is zero. For this system, the ML-
DA also does not improve the forecast Valid Time for $\epsilon=0.01$, as it is
likely to be the case in any system for sufficiently small model error.

Figure 7: The top 3 panels in this figure are spatiotemporal plots with the 64
spatial grid point locations plotted vertically and the Lyapunov time plotted
horizontally. From top to bottom in these panels, the color axis shows the
true value of $\hat{\textbf{x}}_{i}(t),i=1,2,\dots,Q$ from the Kuramoto-
Sivashinky system, the difference between the true
$\hat{\textbf{x}}_{i}(t),i=1,2,\dots,Q$ and the prediction from the baseline
scheme, and the difference between the true
$\hat{\textbf{x}}_{i}(t),i=1,2,\dots,Q$ and the prediction from the ML-DA
scheme. The bottom plot shows the normalized RMS error as function of time for
the baseline (red) and ML-DA (blue) schemes, while the dotted vertical black
lines mark the valid time of prediction using each scheme. For both schemes,
the model error was $\epsilon=0.1$, and the state was measured at $\Theta=16$
equally spaced locations of the full $Q=64$ grid points with measurement noise
magnitude $\sigma_{noise}=0.1$. Figure 8: Forecast Valid Time of the ML-DA
(light blue box plots) and baseline (magenta box plots) schemes at different
values of the model error ($\epsilon$).

####  III.3.2 Analysis Error using Iterative ML-ETKF

Figure 9: Normalized RMS error during DA as we iteratively use the ML-DA to
produce analysis time series and use those time series to train the hybrid
model. A covariance inflation of $\rho=2$, which we found to produce the
lowest analysis error, was used to produce analysis time series at each step
of the iteration. In addition, all knowledge-based models used in this plot
had an error of $\epsilon=0.2$. The plotted values are the median normalized
RMS error values computed over 50 different sets of measurement data. The
variance in the error over these different data sets is very small; therefore,
we have chosen not to plot the standard deviation. At iteration 0, the error
values are those from using only the knowledge-based model as our forecast
model during DA. We also plot the error values computed over only the measured
or the unmeasured variables.

To demonstrate that the ML-DA hybrid model is also capable of improving the
analysis states in a DA cycle, we re-do the DA run that produced the training
data using the hybrid model in place of the knowledge-based model. We then re-
train the hybrid model, and iterate as described in Sec. II.3. Figure 9 shows
the error in analysis states computed using iteratively produced hybrid
models. Here, iteration $0$ denotes the original set of analysis states
computed using only the knowledge-based model. We see that using the first 2
generations of the trained ML-DA to perform DA results in a significant
reduction in error in the analyses of both the measured and unmeasured
variables. Later iterations give smaller reductions in error. We also measured
the forecast Valid Time for each hybrid model trained on the iteratively-
produced analysis time series. We found that there was only a small increase
in the median Valid Time after one iteration ($\sim 10$%), and that successive
iterations did not improve the Valid Time.

##  IV Conclusions and Discussion

We have demonstrated that DA can successfully extend the hybrid approach
introduced in [Pathak  _et al._ , 2018b], which trains a machine-learning
component based on reservoir computing to improve a knowledge-based forecast
model, to the case where a time series of directly-measured model states is
not available for training. Though we use a specific ensemble DA scheme, ETKF,
in our numerical experiments, any DA scheme could be used to process the
available measurements into estimated model states for use as training data.
Using a knowledge-based model rather than a purely data-driven approach allows
us to perform DA before training the machine-learning component. A potential
drawback to this approach is that biases present in the knowledge-based model
will affect the training data, reducing the effectiveness of the training in
correcting for these biases. Using the iterative approach described in Sec.
II.3, however, we are able to reduce the effect of these biases on the
analysis errors by using the previously trained hybrid model as the forecast
model on the DA. In applications like weather forecasting, where multiple
knowledge-based models are available, another approach would be to use a
different model in the DA procedure than in the hybrid training and prediction
procedure.

Another natural extension of the method described in this article would be to
integrate the machine-learning training with the data-assimilation cycle in
such a way that the training is updated each time new data becomes available.
The use of reservoir computing as the machine-learning component of our hybrid
model allows a large number of parameters (the elements of the output weight
matrix) to be trained efficiently using linear regression, which could be done
sequentially (like DA) through techniques such as recursive least squares. It
could be fruitful to more closely integrate the state estimation and the
estimation of the parameters of the hybrid model. The unified Bayesian
framework described in [Bocquet  _et al._ , 2020] and [Bocquet, Farchi, and
Malartic, 2020] could be applied to a hybrid model, though a successful
implementation may require limiting the number of machine-learning parameters.

## Acknowledgement

This work was supported by DARPA grant HR00111890044. The work of Istvan
Szunyogh was also supported by ONR Award N00014-18-2509.

## Appendix

The particular data assimilation algorithm used in this paper is the Ensemble
Transform Kalman Filter (ETKF) [Bishop, Etherton, and Majumdar, 2001; Wang,
Bishop, and Julier, 2004], a variant of the Ensemble Kalman Filter (EnKF)
[Evensen, 2003; Houtekamer and Mitchell, 1998; Burgers, Jan van Leeuwen, and
Evensen, 1998]. Our formulation is based on the local analysis used in the
LETKF [Ott  _et al._ , 2004; Hunt, Kostelich, and Szunyogh, 2007; Szunyogh
_et al._ , 2008] for spatially extended systems larger than the KS system used
in Sec. III.3.

Recall from Sec. II.1 that the analysis state at time $j\Delta t$ is computed
from the analysis state at time $(j-1)\Delta t$ by a forecast step followed by
an analysis step. An EnKF uses an ensemble $\\{\textbf{x}^{a,k}_{j}\\}_{1\leq
k\leq E}$ of analysis states whose mean $\bar{\textbf{x}}^{a}_{j}$ represents
the estimated state and whose covariance $\textbf{P}^{a}_{j}$ represents the
uncertainty in the estimate. We now describe the forecast and analysis steps
for the ensemble.

In an EnKF, the background ensemble members $\textbf{x}^{b,k}_{j}$ are
forecast separately using the model, Eq. (3), so that

$\displaystyle\textbf{x}^{b,k}_{j}=\textbf{x}^{k}(j\Delta
t)\quad\text{where}\quad\frac{d\textbf{x}^{k}}{dt}=\textbf{G}[\textbf{x}^{k}(t)],\quad\textbf{x}^{k}((j-1)\Delta
t)=\textbf{x}^{a,k}_{j-1}.$ (28)

The sample covariance $\textbf{P}_{j}^{b}$ of the background ensemble is

$\displaystyle\textbf{P}^{b}_{j}=(E-1)^{-1}\sum_{i=1}^{E}(\textbf{x}^{b,k}_{j}-\bar{\textbf{x}}^{b}_{j})(\textbf{x}^{b,k}_{j}-\bar{\textbf{x}}^{b}_{j})^{T},$
(29)

where $\bar{\textbf{x}}^{b}_{j}$ is the mean of the background ensemble
members. The goal of the analysis step, in which the types of EnKF differ, is
to compute an analysis ensemble whose mean $\bar{\textbf{x}}^{a}_{j}$ and
sample covariance $\textbf{P}^{a}_{j}$ are consistent with the Kalman filter
[Kalman, 1960].

We now present the computations that form the analysis step of the ETKF, as
formulated in Ref. Hunt, Kostelich, and Szunyogh, 2007 (in which the
consistency with the Kalman filter is explained). Recall that $\textbf{y}_{j}$
is a vector of measurements made at time $j\Delta t$.

  1. 1.

Create the matrix

$\displaystyle\textbf{X}^{b}_{j}$
$\displaystyle=\left[\quad\textbf{x}^{b,1}_{j}-\bar{\textbf{x}}^{b}_{j}\quad|\quad\textbf{x}^{b,2}_{j}-\bar{\textbf{x}}^{b}_{j}\quad|\quad\dots\quad|\quad\textbf{x}^{b,E}_{j}-\bar{\textbf{x}}^{b}_{j}\quad\right],$
(30)

whose columns represent deviations of the ensemble members from the ensemble
mean.

  2. 2.

Compute the matrix $\textbf{Y}^{b}_{j}=\textbf{H}\textbf{X}^{b}_{j}$ and the
vector $\bar{\textbf{y}}_{j}^{b}$.

  3. 3.

Create a matrix

$\displaystyle\tilde{\textbf{P}}^{a}_{j}=\left[(E-1)I/\rho+\textbf{C}\textbf{Y}^{b}_{j}\right]$
(31)

where $\textbf{C}=(\textbf{Y}^{b}_{j})^{T}\textbf{R}^{-1}$ and $\rho$ is the
covariance inflation parameter described in Sec. II.1.

  4. 4.

Compute
$\textbf{W}^{a}_{j}=\left[(E-1)\tilde{\textbf{P}}^{a}_{j}\right]^{1/2}$, where
the $[\cdot]^{1/2}$ denotes the symmetric square root.

  5. 5.

Compute
$\bar{\textbf{w}}^{a}_{j}=\tilde{\textbf{P}}^{a}\textbf{C}(\textbf{y}_{j}-\bar{\textbf{y}}^{b}_{j})$
and add it to the $k^{\text{th}}$ column of $\textbf{W}^{a}_{j}$ to obtain
$\textbf{w}^{a,k}_{j}$.

  6. 6.

Compute
$\textbf{x}^{a,k}_{j}=\textbf{X}^{b}_{j}\textbf{w}^{a,k}_{j}+\bar{\textbf{x}}^{b}_{j}$.
The analysis ensemble is given by $\\{\textbf{x}^{a,k}_{j}\\}_{1\leq k\leq
E}$.

For our results, the analysis ensemble is initialized arbitrarily in a
neighborhood of the model attractor. Then, the forecast and analysis steps
described above are applied iteratively to generate the analysis time series
used to train the ML system.

## References

  * Abarbanel, Rozdeba, and Shirman (2018) Abarbanel, H. D. I., Rozdeba, P. J., and Shirman, S., Neural Computation 30, 2025 (2018).
  * Arcomano  _et al._ (2020) Arcomano, T., Szunyogh, I., Pathak, J., Wikner, A., Hunt, B. R., and Ott, E., Geophysical Research Letters 47, e2020GL087776 (2020).
  * Asch, Bocquet, and Nodet (2016) Asch, M., Bocquet, M., and Nodet, M.,  _Data assimilation: Methods, Algorithms, and Applications_, Fundamentals of Algorithms (SIAM, 2016).
  * Bishop, Etherton, and Majumdar (2001) Bishop, C. H., Etherton, B. J., and Majumdar, S. J., Monthly Weather Review 129, 420 (2001).
  * Bocquet  _et al._ (2020) Bocquet, M., Brajard, J., Carassi, A., and Bertino, L., Foundations of Data Science 2, 55 (2020).
  * Bocquet  _et al._ (2019) Bocquet, M., Brajard, J., Carrassi, A., and Bertino, L., Nonlinear Processes in Geophysics 26, 143 (2019).
  * Bocquet, Farchi, and Malartic (2020) Bocquet, M. B., Farchi, A., and Malartic, Q., Foundations of Data Science 0 (2020), 10.3934/fods.2020015.
  * Brajard  _et al._ (2020) Brajard, J., Carrassi, A., Bocquet, M., and Bertino, L., Journal of Computational Science 44, 101171 (2020).
  * Burgers, Jan van Leeuwen, and Evensen (1998) Burgers, G., Jan van Leeuwen, P., and Evensen, G., Monthly Weather Review 126, 1719 (1998).
  * Chattopadhyay, Hassanzadeh, and Subramanian (2020) Chattopadhyay, A., Hassanzadeh, P., and Subramanian, D., Nonlinear Processes in Geophysics 27, 373 (2020).
  * Evensen (2003) Evensen, G., Ocean Dynamics 53, 343 (2003).
  * Fletcher (2017) Fletcher, S.,  _Data Assimilation for the Geosciences: From Theory to Applications._ (Elsevier, 2017).
  * Hamilton, Berry, and Sauer (2016) Hamilton, F., Berry, T., and Sauer, T., Physical Review X 6, 011021 (2016).
  * Houtekamer and Mitchell (1998) Houtekamer, P. L. and Mitchell, H. L., Monthly Weather Review 126, 796 (1998).
  * Hunt, Kostelich, and Szunyogh (2007) Hunt, B. R., Kostelich, E. J., and Szunyogh, I., Physica D: Nonlinear Phenomena 230, 112 (2007).
  * Jaeger and Haas (2004) Jaeger, H. and Haas, H., Science 304, 78 (2004).
  * Kalman (1960) Kalman, R. E., Journal of Basic Engineering 82, 35 (1960).
  * Kassam and Trefethen (2005) Kassam, A.-K. and Trefethen, L. N., SIAM Journal on Scientific Computing 26, 1214 (2005).
  * Kuramoto and Tsuzuki (1976) Kuramoto, Y. and Tsuzuki, T., Progress of Theoretical Physics 55, 356 (1976).
  * Law, Stuart, and Zygalakis (2015) Law, K., Stuart, A., and Zygalakis, K.,  _Data Assimilation: A Mathematical Introduction_, Vol. 62 (Springer International, 2015).
  * Lguensat  _et al._ (2017) Lguensat, R., Tandeo, P., Ailliot, P., Pulido, M., and Fablet, R., Monthly Weather Review 145, 4093 (2017).
  * Lorenz (1963) Lorenz, E. N., Journal of the Atmospheric Sciences 20, 130 (1963).
  * Lu  _et al._ (2017) Lu, Z., Pathak, J., Hunt, B., Girvan, M., Brockett, R., and Ott, E., Chaos: An Interdisciplinary Journal of Nonlinear Science 27, 041102 (2017).
  * Lukoševičius (2012) Lukoševičius, M., in  _Neural Networks: Tricks of the Trade_ (Springer, 2012) pp. 659–686.
  * Ott  _et al._ (2004) Ott, E., Hunt, B. R., Szunyogh, I., Zimin, A. V., Kostelich, E. J., Corazza, M., Kalnay, E., Patil, D., and Yorke, J. A., Tellus A 56, 415 (2004).
  * Pathak  _et al._ (2018a) Pathak, J., Hunt, B., Girvan, M., Lu, Z., and Ott, E., Physical Review Letters 120, 024102 (2018a).
  * Pathak  _et al._ (2017) Pathak, J., Lu, Z., Hunt, B. R., Girvan, M., and Ott, E., Chaos: An Interdisciplinary Journal of Nonlinear Science 27, 121102 (2017).
  * Pathak  _et al._ (2018b) Pathak, J., Wikner, A., Fussell, R., Chandra, S., Hunt, B. R., Girvan, M., and Ott, E., Chaos: An Interdisciplinary Journal of Nonlinear Science 28, 041101 (2018b).
  * Puskorius and Feldkamp (1994) Puskorius, G. V. and Feldkamp, L. A., IEEE Transactions on Neural Networks 5, 279 (1994).
  * Reichstein  _et al._ (2019) Reichstein, M., Camps-Valls, G., Stevens, B., Jung, M., Denzler, J., Carvalhais, N., and Prabhat,, Nature 566, 195 (2019).
  * Sivashinsky (1977) Sivashinsky, G., Acta Astronautica 4, 1177 (1977).
  * Szunyogh  _et al._ (2008) Szunyogh, I., Kostelich, E. J., Gyarmati, G., Kalnay, E., Hunt, B. R., Ott, E., Satterfield, E., and Yorke, J. A., Tellus A 60, 113 (2008).
  * Tomizawa and Sawada (2020) Tomizawa, F. and Sawada, Y., Geoscientific Model Development (preprint), https://doi.org/10.5194/gmd-2020-211, in review, (2020).
  * Wang, Bishop, and Julier (2004) Wang, X., Bishop, C. H., and Julier, S. J., Monthly Weather Review 132, 1590 (2004).
  * Whitaker and Hamill (2002) Whitaker, J. S. and Hamill, T. M., Monthly Weather Review 130, 1913 (2002).
  * Wikner  _et al._ (2020) Wikner, A., Pathak, J., Hunt, B., Girvan, M., Arcomano, T., Szunyogh, I., Pomerance, A., and Ott, E., Chaos: An Interdisciplinary Journal of Nonlinear Science 30, 053111 (2020).