Physics-informed CoKriging model of a redox flow battery

We begin by describing the GPR method for modeling cell potential $E$ of the RFB as a function of SOC and parameters describing the battery operation conditions and properties of the battery. Defining the input space of the GPR model is an important step. The number of parameters describing the battery and its operating condition is very large, see the Supplementary Material. Ideally, the input space of the GPR model would include all these parameters and SOC such that the GPR model can predict the cell voltage $E$ at any SOC for any type and size of the RFB. However, the training of such GPR model would require vast amounts of computational time and data, with both training time of the GPR model and required data exponentially increasing with the number of input parameters. Since we only have 16 experiments (15 experiments for training to account for a test data) we must reduce the number of input parameters as much as possible. Therefore, we only include the varying operation conditions in the considered experiments and SOC in the input space and use the 0d model to describe the dependence of $E$ on the rest of parameters. Since all 16 experiments are performed with fixed velocity $v=4.17$ m/s and concentrations $C_{V(II)}^{0}=0$ mol/m³, and $C_{V(V)}^{0}=0$ mol/m³, we exclude $v$, $C_{V(II)}^{0}$, and $C_{V(V)}^{0}$ from the input space. The resulting input space of the GPR model is

with $d=10$. The cell voltage, $E$, is the output (or state of interest) of the GPR model, $E\;:\;\mathbb{D}\rightarrow\mathbb{R}$. In other words, the GPR model aims to predict $E$ for any $s$, $j$, $V_{r}$, $w_{m}$, $C_{V(III)}^{0}$, $C_{V(IV)}^{0}$, $C_{H^{+},pos}^{0}$, $C_{H_{2}O,pos}^{0}$, $C_{H^{+},neg}^{0}$, and $C_{H_{2}O,neg}^{0}$ with fixed $v=4.17$ m/s and $C_{V(II)}^{0}=C_{V(V)}^{0}=0$ mol/m³.

Next, we organize the measurements of $E$ at locations $\mathbf{X}=\{{\mathbf{x}^{(i)}}\}_{i=1}^{N}$ in the observation vector $\mathbf{E}=(E^{1},E^{2},\ldots,E^{N})^{T}$, where $N=15\cdot N_{s}$ is the total number of observations for fifteen experiments, $N_{s}$ is the size of the $\mathbf{s}^{*}$ vector of SOC measurements, and $E^{i}$ is the $E$ measurement at the location $\mathbf{x}^{(i)}$.

In the GPR approach we treat $E(\mathbf{x})$ as a Gaussian process as described in the Methods section. Estimating the mean $\mu({\mathbf{x}})$ and covariance $k({\mathbf{x}},{\mathbf{x}}^{\prime})$ models is the main challenge in GPR-based methods (see Methods.) In the standard data-driven GPR, $\mu({\mathbf{x}})$ and $k({\mathbf{x}},{\mathbf{x}}^{\prime})$ would be computed from data by minimizing the so-called marginal log-likelihood function  [15]. It should be noted that the prior statistics (i.e., $\mu({\mathbf{x}})$ and $k({\mathbf{x}},{\mathbf{x}}^{\prime})$) computed in this approach do not satisfy the physical laws governing the system of interest unless the phase space is very densely sampled, which can lead to significant errors and nonphysical behavior [11, 12]. To address this limitation of the data-driven GPR method, in the following Section II.2 we propose to estimate the prior statistics from the 0d model where we treat unknown parameters as random variables. A similar approach was used in [11, 12] in GPR models of subsurface transport and power grid dynamics.

The motivation behind this approach is that the 0d model, which is described in Methods, has the unknown parameters $k_{1},k_{2},\sigma_{e}$, and $S$ denoting the reference rate constants, ionic concentration of the electrolyte, and specific surface area, which depend on operating conditions and must be calibrated for each experiment. In this work we treat these parameters as Gaussian random variables with prescribed (prior) mean and variance that turns the 0d model in a stochastic model. Next, we use the Monte Carlo (MC) method to compute the prior GP model of $E$ based on the stochastic 0d model. Specifically, we generate $N_{MC}$ realizations of the $k_{1},k_{2},\sigma_{e}$, and $S$ parameters and compute $E$ for each realization of the parameters from the 0d model as described in Methods. Then, the prior mean and covariance of $E$ are calculated from this set of $\{E^{(i)}_{L}\}_{i=1}^{N_{MC}}$ realizations as

Here, the superscript $L$ signifies that the corresponding quantities are obtained from a numerical model that is based on approximation and has a lower fidelity than the experimental data. The estimates of $E$ in terms of the low-fidelity prior statistics are given as:

where the components of ${\mathbf{c}}_{L}({\mathbf{x}^{*}})$ and ${\mathbf{C}}_{L}$ are given by $k_{L}$.

When the physics-based model of a considered system is accurate the prior statistics computed from Eqs (5) and (6) were shown to provide an accurate model of the system given the system’s measurements [11, 12]. However, if the physics-based model cannot capture certain features of the system (as is the case with the 0d model), then the prior statistics computed would produce a poor GPR model of the system [10]. In the following section we demonstrate how the discrepancy between the observation data and the 0d model can be accounted for in the GPR prediction.

The systematic errors in the 0d model (also called model bias or misspecification) can be corrected by establishing the correlation between the 0d model predictions and the experimental data using the CoKriging method [10, 16]. This idea was originally formulated for systems where the data are sparse, so measurements of other states or other parameters are additionally used to model the data [17, 18]. CoKriging has recently been used with success for data-driven multifidelity modeling to predict system responses using data with two levels of fidelity [19, 20] adapted from Kennedy and O’Hagan’s CoKriging framework [16].

In the multifidelity framework, we treat $E$ observed in the experiments as high-fidelity data. We then use the 0d model to generate the low-fidelity data $\mathbf{E}_{L}=\left(E_{L}^{(1)},\ldots,E_{L}^{(N_{L})}\right)^{T}$ at locations $\mathbf{X}_{L}=\{{\mathbf{x}^{(i)}_{L}}\}_{i=1}^{N_{L}}$, where $E_{L}^{(i)}\in{\mathbb{R}}$, and ${\mathbf{x}}_{L}^{(i)}\in{\mathbb{D}}\subseteq{\mathbb{R}}^{d}$. Note that in this work we generate $E_{L}^{(i)}$ using MC simulations (see Section V.2) and $N_{L}=N_{MC}$, the number of realizations in the MC simulations.

Following [10, 16], we define an auto-regressive model for the high-fidelity data, $E({\mathbf{x}})=\rho E_{L}({\mathbf{x}})+E_{D}({\mathbf{x}})$, where $E_{L}({\mathbf{x}})$ is the GP model of the low-fidelity data, $E_{D}({\mathbf{x}})$ is the GP model describing the mismatch between the high and low fidelity models, and $\rho$ is the regression parameter. The observation covariance matrix for this model is

where the elements of the ${\mathbf{C}}_{L}$ matrices are computed as $k_{L}(\mathbf{x},\mathbf{x}^{\prime})$, where $k_{L}$ is the covariance function of $E_{L}$ and the elements of the ${\mathbf{C}}_{D}$ matrix are given by $k_{D}(\mathbf{x},\mathbf{x}^{\prime})$, where $k_{D}$ is the covariance function of $E_{D}$.

The posterior mean and variance of $E$ at any point ${\mathbf{x}}^{*}\in{\mathbb{D}}$ conditioned on the measurements of $\mathbf{E}$ and $\mathbf{E}_{L}$ are given by

where $\tilde{\mathbf{E}}=(\mathbf{E}_{L}^{T},\mathbf{E}^{T})^{T}$, $\mu({\mathbf{x}}^{*})=\rho\mu_{L}({\mathbf{x}}^{*})+\mu_{D}({\mathbf{x}}^{*})$, $\sigma_{L}^{2}({\mathbf{x}}^{*})=k_{L}({\mathbf{x}}^{*},{\mathbf{x}}^{*})$, $\sigma_{D}^{2}({\mathbf{x}}^{*})=k_{D}({\mathbf{x}}^{*},{\mathbf{x}}^{*})$, $\tilde{\bm{\mu}}=(\bm{\mu}_{L}^{T},\bm{\mu}^{T})^{T}$, $\tilde{\mathbf{c}}({\mathbf{x}}^{*})=(\rho{\mathbf{c}}_{L}({\mathbf{x}}^{*})^{T},{\mathbf{c}}({\mathbf{x}}^{*})^{T})^{T}$, and $k({\mathbf{x}},{\mathbf{x}}^{\prime})=\rho^{2}k_{L}({\mathbf{x}},{\mathbf{x}}^{\prime})+k_{D}({\mathbf{x}},{\mathbf{x}}^{\prime}).$

In the data-driven approaches of [16, 17, 18, 19, 20], the prior statistics (i.e., $\mu_{L}({\mathbf{x}}^{*})$, $\mu_{D}({\mathbf{x}}^{*})$, $k_{L}({\mathbf{x}},{\mathbf{x}}^{\prime})$, and $k_{D}({\mathbf{x}},{\mathbf{x}}^{\prime})$) are learned by maximizing the likelihood function of $Y$. Here, we compute $\mu_{L}({\mathbf{x}}^{*})$ and $k_{L}({\mathbf{x}},{\mathbf{x}}^{\prime})$ from the Monte Carlo simulations based on the 0d model with random parameters $k_{1},k_{2},\sigma_{e}$, and $S$ as described in Section II.2. Then, we compute $\mu_{L}({\mathbf{x}}^{*})$ and $k_{L}({\mathbf{x}},{\mathbf{x}}^{\prime})$ using the data set $\mathbf{E_{D}}=\mathbf{E}-\rho\mathbf{E}_{L}^{X}$ with the corresponding locations $\mathbf{X}$ (where $\mathbf{E}_{L}^{X}=[E_{L,1},E_{L,1},...,E_{L,N}]^{T}$ and $E_{L,i}=\mu_{L}(\mathbf{x}_{i})$ ($\mathbf{x}_{i}\in\mathbf{X}$)) by minimizing the log likelihood:

Here, we set $\rho=1.0$ and assume that $E_{D}(\mathbf{x})$ has the prior exponential covariance function

and constant $\mu_{D}$. In Eq (10), the $(i,j)$ component of the covariance matrix $\mathbf{C}_{D}$ is given by $k_{D}(\mathbf{x}_{i},\mathbf{x}_{j})$ ($(\mathbf{x}_{i},\mathbf{x}_{j})\in\mathbf{X}_{D}$).

A key step in the data preparation for training the CoKriging model is to rescale SOC in different experiments according to Eqs. (1) and (2). For data from an experiment $i$, the scaling parameters in Eqs. (1) and (2) (i.e., $SOC_{i}(0)$ and $\max(SOC_{i})$) are given by the experimental data. Furthermore, for a given set of input parameters, the CoKriging model predicts voltage as a function of the scaled SOC. However, to covert the scaled SOC to time, one needs to evaluate scaling parameters as functions of the input parameters. In the following, we introduce two scaling parameters denoted by $\theta_{C}=1/(\max(SOC_{C})-SOC_{C}(0))$ and $\theta_{D}=1/(\max(SOC_{D})-SOC_{D}(0))$. When $\max(SOC_{C})$ and $\max(SOC_{D})$ are known these scaling parameters are trivial to find. However, we may not know $SOC_{C}(0)$ and $SOC_{D}(0)$ for arbitrary conditions, so we first need to estimate them or, equivalently, the scaling parameters. We compute $\theta_{C}$ using the standard GPR model trained on a data set consisting of the scaling parameters $\theta_{C}$ and corresponding initial voltages $E_{0}=E(t=0)$ from the known experiments. We note that in the considered lab experiments, $-\theta_{D}=0.05+\theta_{C}$, which allows us to estimate $\theta_{D}$ given $\theta_{C}$ without needing the voltage at any SOC during the discharge cycle. We assumed that there is uncertainty in the $\theta_{C}$ values obtained from the experiments with the variance $\sigma^{2}_{\theta}=0.5$ The non-zero $\sigma^{2}_{\theta}$ reflects the fact that the different values of $\theta_{C}$ were obtained from the experiments with very similar values of $E_{0}$. The positive value of $\sigma^{2}_{\theta}$ is also needed to make the measurement covariance matrix well-conditioned. Figure 3 shows the GPR model prediction of $\theta_{C}$ as a function of $E_{0}$ and the prediction uncertainty. The GPR model uncertainty increase outside of the range $1.3-1.5$V of $E_{0}$, where the data are available. This uncertainty can be reduced by conducting experiments with the expended range of $E_{0}$. Within $1.3-1.5$V range, the uncertainty in GPR model prediction can be reduced by considering additional variables that affect $\theta_{c}$.

In this section we conduct CoPhIK simulations as described above, except the input parameters to the Gaussian random variables in Eqs. (30)–(33) are taken from literature [8] instead of [21]. The rest of the data remains the same. Fig. 8 demonstrates that $E$ computed from the 0d model with these parameters overestimates the experimental values during charging and underestimates the discharging section by approximately 10$\%$. As a consequence, the mean of the low-fidelity dataset also over- and underestimates the experimental data. However, the results with CoPhIK can compensate for the discrepancies between the experimental and low-fidelity data. In Fig. 6, we show the 0d and CoPhIK predictions of $E$ for the same experimental data set as in Fig. 4. The mean of the low-fidelity dataset, $\langle E\rangle$, does not agree well with the experimental data. However, the CoPhIK prediction agrees very well. In some cases, the uncertainty is higher than in Fig. 4, but the absolute error is of the same order of magnitude with both sets of input parameters to the 0d model. The average $\ell_{2}$ errors are $4.09\times 10^{-4}$, $1.64\times 10^{-4}$, and $9.44\times 10^{-5}$ with one, two, and three data points, respectively, and the average $\ell_{\infty}$ errors are $8.85\times 10^{-2}$, $9.66\times 10^{-2}$, and $9.99\times 10^{-2}$ with one, two, and three data points. Comparing with Table 1, there is no significant difference between using the input parameters in this work and the input parameters from [8] for the $\ell_{2}$ error, and results are similar for the $\ell_{\infty}$ error. Very similar results are seen when using the parameters from [9].

These results show that the CoPhIK framework removes the potential bias because of the difference in the prior and actual statistics of the parameters in the 0d model, i.e., the proposed CoPhIK framework is non-sensitive to the choice of the prior statistics of the parameters in the 0d model.

In GPR we wish to find the cell voltage $E\in{\mathbb{R}}$ at input variables $\mathbf{x}\in{\mathbb{D}}$, $E(\mathbf{x})$, as a Gaussian process:

where $\mu(\cdot):{\mathbb{D}}\rightarrow{\mathbb{R}}$ and and $k(\cdot,\cdot):{\mathbb{D}}\times{\mathbb{D}}\rightarrow{\mathbb{R}}$ are the mean and covariance function, given by

where the operator $\mathbb{E}\{\cdot\}$ is the ensemble mean. We denote the covariance matrix to be the $N\times N$ matrix given by

and define the covariance vector as

The GPR prediction of $E$ at ${\mathbf{x}}^{*}$ is then given by the posterior distribution $\hat{E}_{p}({\mathbf{x}}^{*})\sim\mathcal{N}(\mu_{p}({\mathbf{x}}^{*}),\sigma^{2}_{p}({\mathbf{x}}^{*}))$, where the posterior mean and variance are given by

where $\bm{\mu}=(\mu(\mathbf{x}^{(1)}),\mu(\mathbf{x}^{(2)}),\ldots,\mu(\mathbf{x}^{(N)}))^{T}$ and $\sigma^{2}({\mathbf{x}^{*}})=k(\mathbf{x},\mathbf{x})$.

To ensure that ${\mathbf{C}}$ is invertible we substitute ${\mathbf{C}}={\mathbf{C}}+\epsilon\mathbf{I}$, where $\epsilon\ll 1$. This is equivalent to assuming there is i.i.d. observation noise with variance $\epsilon$ [10].

The 0d model of an RFB was originally proposed in [22, 8] and modified by [9] to give a fast estimate of redox-flow battery performance. The 0d model gives the cell voltage, $E$, due to an applied current density, $j_{app}$, in the form

where $E_{cell}^{rev}$ is the reversible Open Current Voltage (OCV), $\eta_{act}$ is the activation overpotential, and $\eta_{ohm}$ is the ohmic losses. A short explanation of the 0d model is given below.

The two half reactions in the RFB are given by

From Nernst’s equation [9],

where $C_{i}$ is the molar concentration of species $i$, $E_{1}^{0}$ and $E_{2}^{0}$ are the formal potentials for the reactions at the negative and positives electrodes, $R$ is the molar gas constant, $T$ is the cell temperature, and $F$ is the Faraday constant.

The activation overpotential is given by

where

and

The coefficients $k_{1}$ and $k_{2}$ are the reference rate constants for reactions (20) and (21) at $T_{ref}=293K$, respectively, and are typically determined from experiments.

The ohmic losses are given by

with the current collector conductivity $\sigma_{c}$, the membrane conductivity

where $\lambda=22$ is the membrane water content for a fully saturated membrane, and the electrolyte ionic conductivity $\sigma_{e}$ to be determined from experimental results. The widths of the membrane, electrolyte, and current collector are denoted by $w_{m}$, $w_{e}$, and $w_{c}$, respectively. The specific area is given by $A_{s}=SV_{e}$, where $S$ is the specific surface area that is determined from experiments and $V_{e}$ is the electrode volume. Mass balance for each species is used to determine $C_{i}$ for $i=V(II),$ $V(III),$ $V(IV),$ $V(V),$ $H^{+},$ and $H_{2}O$. For details see [22, 8].

While most of the parameters in the 0d model can be directly measured or estimated based on the battery design and material properties, the parameters $k_{1},k_{2},\sigma_{e}$, and $S$ denoting the reference rates for reactions (20) and (21), the electrode conductivity, and the specific surface area, must be calibrated [8, 9, 21]. Table 2 presents values of these parameters from the calibration studies in [8, 9, 21]. The fixed parameters are given in Table 3.

A comparison of the 0d model predictions with these three sets of parameters for two experiments is given in Fig. 7. For all considered experiments, the parameters in [9] underestimate $E$, while the parameters in [8] overestimate $E$. The parameters from [21] were obtained from an experiment in a cell similar to the one used in the experiments considered in this work, and so they provide the best agreement in Fig. 7 and serve as the starting point for the MC simulations.

To model the $i$th experiment ($i=1,...,16$), in the MC simulation method, we solve the 0d equations $N_{MC}=1000$ times with the input parameters

taken from Table 5 for experiment $i$. The parameters $\mathbf{v}=[k_{1},k_{2},\sigma_{e},S]^{T}$ are taken as Gaussian random variables with the mean values $\langle v_{i}\rangle$ given by [21] and standard deviation $\sigma_{v_{i}}=0.25\langle v_{i}\rangle$:

where $\xi^{m}(\mu,\sigma)$ denotes a Gaussian random variable with mean $\mu$ and standard deviation $\sigma.$ Fig. 8 shows 50 realizations of $E(t)$ computed from the 0d model for one of the experiments as well as $E(t)$ observed in the experiment and computed from the 0d model with $\mathbf{v}=\langle\mathbf{v}\rangle$. Because the 0d model involves logarithm functions, for some values of $\mathbf{v}$ the output may be undefined. Therefore, we only consider MC realizations where the value is defined for all time. This figure shows that the mean cell voltage $\mu_{L}$ obtained from the MC method cannot accurately predict the tails of the charge and discharge curves. Similar results are obtained for the other experiments.

The experiments use Nafion 115 and 212 membranes. The positive and negative electrodes are made of 3 mm thick graphite felt with active areas of 10 cm². The electrolytes are prepared by dissolving 1.5 M $VOSO_{4}$ (Aldrich, 99%) in 3.5 M $H_{2}SO_{4}$ solution (Aldrich, 96–98%). The electrolytes are placed in two glass reservoirs and circulated at a flow rate of 20 mL/min by a peristaltic pump. The flow cell electrochemical performance are measured and recorded using a potentiostat/galvanostat (Arbin Instrument, USA) with the fixed voltage range of 0.8 to 1.6 V. Table 5 lists parameters that were varied between the experiments while the parameters that were kept constant in all experiments are given in Table 4. The current density in each experiment is kept constant, although it varied between experiences.