3.1 Challenges

Online model calibration for WFSim is challenging for a number of reasons. First of all, the model is nonlinear, and thus the common linear estimation algorithms cannot be used without linearization, which limits accuracy (Boersma et al., 2018). Secondly, an estimation solution relying on WFSim is sensitive to instability when the model state sufficiently deviates from the continuity equation. Finally, the surrogate model typically has on the order of N∼10³–10⁴ states, which is extraordinarily high for control applications. However, real-time estimation is a necessity for real-time model-based control, and thus one needs to find a trade-off between accuracy while guaranteeing state updates at a low computational cost.

3.2 General formulation

This section summarizes the basics of the KF, which is the literature standard for state estimation in control. The goal of a KF is to recursively estimate the unmeasured states of a dynamical system through noisy measurements. Assumed here is a system (the wind farm) represented mathematically by a discrete-time stochastic state-space model with additive noise,

where k is the time index, x∈ℝ^N is the unobserved system state, z∈ℝ^M are the measured outputs of the system, q∈ℝ^O and w∈ℝ^N are the controllable inputs and process noise, respectively, that drive the system dynamics, and v∈ℝ^M is measurement noise. Furthermore, we assume w and v to be zero-mean white Gaussian noise with covariance matrices

with E the expectation operator. Estimates of the state x_k, denoted by ${\widehat{\mathit{x}}}_{k|k}$, are computed based on measurements from the real system. Here, ${\widehat{\mathit{x}}}_{k|\mathrm{\ell }}$ means an estimate of the state vector x at time k, using all past measurements and inputs 𝓩_ℓ as

State estimates are based on the internal model dynamics and the measurements, weighted according to their probability distributions. We aim to find an optimal state estimate, in which optimality is defined as unbiasedness, $E[{\mathit{x}}_k-{\widehat{\mathit{x}}}_k]=\mathrm{0}$, and when the variance of any linear combination of state estimation errors (e.g., the trace of $E\left[\right({\mathit{x}}_k-{\widehat{\mathit{x}}}_k\left)\right({\mathit{x}}_k-{\widehat{\mathit{x}}}_k{)}^{\mathrm{T}}]$) is minimized.

In reality, the assumed model described by f and h always has mismatches with the true system, and the assumptions in Eq. (3) often do not hold. Further, the matrices Q_k, R_k, and S_k are usually not known and are rather considered to be tuning parameters, used to shift the confidence levels between the internal model and the measured values. For R≪Q, estimations will heavily rely on the measurements, while for Q≪R, estimations will mostly rely on the internal model. Kalman filtering remains one of the most common methods of recursive state estimation. KF algorithms typically consist of two steps.

1. A state and output forecast, including their uncertainties (covariances):

   $$\begin{array}{ll}\text{(5)}& {\displaystyle }{\widehat{\mathit{x}}}_{k|k-\mathrm{1}}& {\displaystyle }=E\left[f({\mathit{x}}_{k-\mathrm{1}},{\mathit{q}}_{k-\mathrm{1}})+{\mathit{w}}_{k-\mathrm{1}}|{\mathcal{Z}}_{k-\mathrm{1}}\right],\text{(6)}& {\displaystyle }{\widehat{\mathit{z}}}_{k|k-\mathrm{1}}& {\displaystyle }=E\left[h({\mathit{x}}_k,{\mathit{q}}_k)+{\mathit{v}}_k|{\mathcal{Z}}_{k-\mathrm{1}}\right],{\displaystyle }{\mathbf{P}}_{k|k-\mathrm{1}}^{\mathit{x}}& {\displaystyle }=\mathrm{Cov}\left({\mathit{x}}_k,{\mathit{x}}_k|{\mathcal{Z}}_{k-\mathrm{1}}\right)\\ \text{(7)}& {\displaystyle }& {\displaystyle }=E\left[\right({\mathit{x}}_k-{\widehat{\mathit{x}}}_{k|k-\mathrm{1}}\left)\right({\mathit{x}}_k-{\widehat{\mathit{x}}}_{k|k-\mathrm{1}}{)}^{\mathrm{T}}],{\displaystyle }{\mathbf{P}}_{k|k-\mathrm{1}}^{\mathit{z}}& {\displaystyle }=\mathrm{Cov}\left({\mathit{z}}_k,{\mathit{z}}_k|{\mathcal{Z}}_{k-\mathrm{1}}\right)\\ \text{(8)}& {\displaystyle }& {\displaystyle }=E\left[\right({\mathit{z}}_k-{\widehat{\mathit{z}}}_{k|k-\mathrm{1}}\left)\right({\mathit{z}}_k-{\widehat{\mathit{z}}}_{k|k-\mathrm{1}}{)}^{\mathrm{T}}],{\displaystyle }{\mathbf{P}}_{k|k-\mathrm{1}}^{\mathit{x}\mathit{z}}& {\displaystyle }=\mathrm{Cov}\left({\mathit{x}}_k,{\mathit{z}}_k|{\mathcal{Z}}_{k-\mathrm{1}}\right)\\ \text{(9)}& {\displaystyle }& {\displaystyle }=E\left[\right({\mathit{x}}_k-{\widehat{\mathit{x}}}_{k|k-\mathrm{1}}\left)\right({\mathit{z}}_k-{\widehat{\mathit{z}}}_{k|k-\mathrm{1}}{)}^{\mathrm{T}}].\end{array}$$

   In Eqs. (5) and (6), ${\widehat{\mathit{x}}}_{k|\mathrm{\ell }}$ and ${\widehat{\mathit{z}}}_{k|\mathrm{\ell }}$ are the forecasted system state vector and measurement vector, respectively.

2. An analysis update of the state vector, where the measurements are fused with the internal model:

   $$\begin{array}{}\text{(10)}& {\displaystyle }& {\displaystyle }{\mathbf{L}}_k={\mathbf{P}}_{k|k-\mathrm{1}}^{\mathit{x}\mathit{z}}\cdot {\left({\mathbf{P}}_{k|k-\mathrm{1}}^{\mathit{z}}\right)}^{-\mathrm{1}}\text{(11)}& {\displaystyle }& {\displaystyle }{\widehat{\mathit{x}}}_{k|k}={\widehat{\mathit{x}}}_{k|k-\mathrm{1}}+{\mathbf{L}}_k\left({\mathit{z}}_k-{\widehat{\mathit{z}}}_{k|k-\mathrm{1}}\right),\text{(12)}& {\displaystyle }& {\displaystyle }{\mathbf{P}}_{k|k}^{\mathit{x}}=\mathrm{Cov}\left({\mathit{x}}_k,{\mathit{x}}_k|{\mathcal{Z}}_k\right)={\mathbf{P}}_{k|k-\mathrm{1}}^{\mathit{x}}-{\mathbf{L}}_k{\mathbf{P}}_{k|k-\mathrm{1}}^{\mathit{z}}{\mathbf{L}}_k^{\mathrm{T}}.\end{array}$$

   Here, ${\left({\mathbf{P}}_{k|k-\mathrm{1}}^{\mathit{z}}\right)}^{-\mathrm{1}}$ in Eq. (10) is the pseudo-inverse of ${\mathbf{P}}_{k|k-\mathrm{1}}^{\mathit{z}}$, since this matrix is not necessarily invertible.

Traditionally, state estimation for linear dynamic models is done using the linear KF (Kalman, 1960). However, this is not a viable option here, as the surrogate model is nonlinear. Rather, a number of nonlinear KF variants are looked upon.

3.3 Extended Kalman filter (ExKF)

Linearization of the surrogate model is the most popular and straight-forward solution to the issue of model nonlinearity, as done in the extended KF (ExKF). The ExKF has shown success in academia and industry (Wan and Van Der Merwe, 2000) and is perhaps the most popular nonlinear KF. However, it has a number of disadvantages. As described in Sect. 3.1, model linearization is troublesome. Furthermore, for surrogate models with many states such as WFSim, the ExKF has an additional challenge: computational complexity. The operation in Eq. (10) includes a matrix inversion with a computational complexity of 𝒪(M³), and the ExKF furthermore includes two matrix multiplications each with a complexity of 𝒪(N³). As there are significantly fewer measurements than states (M≪N) for the problem at hand, these matrix multiplications dominate the computational cost. The ExKF has a CPU time on the order of 10¹ s for a two-turbine wind farm, which may be too large for our purposes. To reduce computational cost in the ExKF, the surrogate model and/or the covariance matrix P have to be simplified. This is not further explored here. Instead, two KF approaches will be explored that directly use the nonlinear system for forecasting and analysis updates. Doing so, we circumvent the problems with linearization, and additionally better maintain the true covariance of the system state.

3.4 The unscented Kalman filter (UKF)

The unscented Kalman filter (UKF) relies on the so-called “unscented transformation” to estimate the means and covariance matrices described by Eqs. (5) to (9). The conditional state probability distribution of x_k knowing 𝓩_k is again assumed to be Gaussian. In the UKF, firstly a number of sigma points (also referred to as “particles”) are generated such that their mean is equal to ${\widehat{\mathit{x}}}_{k|k}$ and their covariance is equal to Cov(x_k,x_k). Secondly, each particle is propagated through the nonlinear system dynamics (f, h). Thirdly, the mean and covariance of the forecasted state probability distribution is again approximated by a weighted mean of these forecasted sigma points (Wan and Van Der Merwe, 2000).

Mathematically, we define the ith particle as ${\mathit{\psi }}_{k|\mathrm{\ell }}^i\in {\mathbb{R}}^N$, which is a realization of the conditional probability distribution of x_k given 𝓩_ℓ. The UKF follows a very similar forecast and analysis update approach as the traditional KF in Eqs. (5) to (12), yet applied to a finite set of particles (Wan and Van Der Merwe, 2000).

1. For the forecast step, a particle-based approach is taken.

   • i.

     A total of $Y=\mathrm{2}N+\mathrm{1}$ particles, with N equal to the state dimension, are (re)sampled to capture the mean and covariance of the conditional state probability distribution $p\left[{\mathit{x}}_{k-\mathrm{1}}|{\mathcal{Z}}_{k-\mathrm{1}}\right]$ by

     $$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathit{\psi }}_{k-\mathrm{1}|k-\mathrm{1}}^i=\\ \text{(13)}& {\displaystyle }& {\displaystyle }\left\{\begin{array}{l}{\stackrel{\mathrm{‾}}{\mathit{\psi }}}_{k-\mathrm{1}|k-\mathrm{1}}\\ \text{ for }i=\mathrm{1},\\ {\stackrel{\mathrm{‾}}{\mathit{\psi }}}_{k-\mathrm{1}|k-\mathrm{1}}+{\left(\sqrt{(N+\mathit{\lambda })\cdot {\mathbf{P}}_{k-\mathrm{1}|k-\mathrm{1}}^{\mathit{x}}}\right)}_i\\ \text{ for }i=\mathrm{2},\mathrm{\dots },N+\mathrm{1},\\ {\stackrel{\mathrm{‾}}{\mathit{\psi }}}_{k-\mathrm{1}|k-\mathrm{1}}-{\left(\sqrt{(N+\mathit{\lambda })\cdot {\mathbf{P}}_{k-\mathrm{1}|k-\mathrm{1}}^{\mathit{x}}}\right)}_{i-N-\mathrm{1}}\\ \text{ for }i=N+\mathrm{2},\mathrm{\dots },Y,\end{array}\right.\end{array}$$

     where $\mathit{\lambda }={\mathit{\alpha }}^{\mathrm{2}}\left(N+\mathit{\kappa }\right)-N$ is a scaling parameter, α determines the spread of the particles around the mean, and κ is a secondary scaling parameter typically set to 0 (Wan and Van Der Merwe, 2000). The vector ${\stackrel{\mathrm{‾}}{\mathit{\psi }}}_{k-\mathrm{1}|k-\mathrm{1}}$ is the estimated state vector calculated as ${\stackrel{\mathrm{‾}}{\mathit{\psi }}}_{k-\mathrm{1}|k-\mathrm{1}}={\sum }_{i=\mathrm{1}}^Y\left({\mathit{w}}_{\text{mean}}^i\cdot {\mathit{\psi }}_{k-|k-\mathrm{1}}^i\right)$, where the weight of each particle's mean ${\mathit{w}}_{\text{mean}}^i$ and covariance ${\mathit{w}}_{\text{cov.}}^i$ is given by

     $$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathit{w}}_{\text{mean}}^i=\left\{\begin{array}{ll}\mathit{\lambda }(N+\mathit{\lambda }{)}^{-\mathrm{1}}& \text{ for }i=\mathrm{1},\\ \frac{\mathrm{1}}{\mathrm{2}}(N+\mathit{\lambda }{)}^{-\mathrm{1}}& \text{ otherwise},\end{array}\right.\\ {\displaystyle }& {\displaystyle }{\mathit{w}}_{\text{cov.}}^i=\left\{\begin{array}{ll}\mathit{\lambda }(N+\mathit{\lambda }{)}^{-\mathrm{1}}+(\mathrm{1}-{\mathit{\alpha }}^{\mathrm{2}}+\mathit{\beta })& \text{ for }i=\mathrm{1},\\ \frac{\mathrm{1}}{\mathrm{2}}(N+\mathit{\lambda }{)}^{-\mathrm{1}}& \text{ otherwise},\end{array}\right.\end{array}$$

     and β is used to incorporate prior knowledge on the probability distribution. In this work, β=2 is assumed, which is stated to be optimal for Gaussian distributions (Wan and Van Der Merwe, 2000).

   • ii.

     Each particle is propagated forward in time using the expectation of the nonlinear model as

     $$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathit{\psi }}_{k|k-\mathrm{1}}^i=f({\mathit{\psi }}_{k-\mathrm{1}|k-\mathrm{1}}^i,{\mathit{q}}_{k-\mathrm{1}})\text{ for }i=\mathrm{1},\mathrm{\dots },Y,\\ \text{(14)}& {\displaystyle }& {\displaystyle }{\mathit{\zeta }}_{k|k-\mathrm{1}}^i=h({\mathit{\psi }}_{k|k-\mathrm{1}}^i,{\mathit{q}}_k)\text{ for }i=\mathrm{1},\mathrm{\dots },Y,\end{array}$$

     where ${\mathit{\zeta }}_{k|\mathrm{\ell }}^i$ is defined as the system output corresponding to the particle ${\mathit{\psi }}_{k|\mathrm{\ell }}^i$.

   • iii.

     The expected state $\stackrel{\mathrm{‾}}{\mathit{\psi }}$ and expected output $\stackrel{\mathrm{‾}}{\mathit{\zeta }}$ are calculated as

     $$\begin{array}{ll}{\displaystyle }{\widehat{\mathit{x}}}_{k|k-\mathrm{1}}& {\displaystyle }={\stackrel{\mathrm{‾}}{\mathit{\psi }}}_{k|k-\mathrm{1}}=\sum _{i=\mathrm{1}}^Y\left({\mathit{w}}_{\text{mean}}^i\cdot {\mathit{\psi }}_{k|k-\mathrm{1}}^i\right),\\ \text{(15)}& {\displaystyle }{\widehat{\mathit{z}}}_{k|k-\mathrm{1}}& {\displaystyle }={\stackrel{\mathrm{‾}}{\mathit{\zeta }}}_{k|k-\mathrm{1}}=\sum _{i=\mathrm{1}}^Y\left({\mathit{w}}_{\text{mean}}^i\cdot {\mathit{\zeta }}_{k|k-\mathrm{1}}^i\right),\end{array}$$

     and the covariance matrices are (re-)estimated from the forecasted ensemble by

     $$\begin{array}{ll}{\displaystyle }{\mathbf{P}}_{k|k-\mathrm{1}}^{\mathit{x}}=\sum _{i=\mathrm{1}}^Y& {\displaystyle }({\mathit{w}}_{\text{cov.}}^i\left({\mathit{\psi }}_{k|k-\mathrm{1}}^i-{\stackrel{\mathrm{‾}}{\mathit{\psi }}}_{k|k-\mathrm{1}}\right)\\ \text{(16)}& {\displaystyle }& {\displaystyle }{\left({\mathit{\psi }}_{k|k-\mathrm{1}}^i-{\stackrel{\mathrm{‾}}{\mathit{\psi }}}_{k|k-\mathrm{1}}\right)}^{\mathrm{T}})+{\mathbf{Q}}_{k-\mathrm{1}},{\displaystyle }{\mathbf{P}}_{k|k-\mathrm{1}}^{\mathit{z}}=\sum _{i=\mathrm{1}}^Y& {\displaystyle }({\mathit{w}}_{\text{cov.}}^i\left({\mathit{\zeta }}_{k|k-\mathrm{1}}^i-{\stackrel{\mathrm{‾}}{\mathit{\zeta }}}_{k|k-\mathrm{1}}\right)\\ \text{(17)}& {\displaystyle }& {\displaystyle }{\left({\mathit{\zeta }}_{k|k-\mathrm{1}}^i-{\stackrel{\mathrm{‾}}{\mathit{\zeta }}}_{k|k-\mathrm{1}}\right)}^{\mathrm{T}})+{\mathbf{R}}_k,{\displaystyle }{\mathbf{P}}_{k|k-\mathrm{1}}^{\mathit{x}\mathit{z}}=\sum _{i=\mathrm{1}}^Y& {\displaystyle }({\mathit{w}}_{\text{cov.}}^i\left({\mathit{\psi }}_{k|k-\mathrm{1}}^i-{\stackrel{\mathrm{‾}}{\mathit{\psi }}}_{k|k-\mathrm{1}}\right)\\ \text{(18)}& {\displaystyle }& {\displaystyle }{\left({\mathit{\zeta }}_{k|k-\mathrm{1}}^i-{\stackrel{\mathrm{‾}}{\mathit{\zeta }}}_{k|k-\mathrm{1}}\right)}^{\mathrm{T}})+{\mathbf{S}}_k.\end{array}$$

2. For the analysis step, one can apply the same equations as in Eqs. (10) to (12).

The UKF has been shown to consistently outperform the ExKF in terms of accuracy, since it uses the nonlinear model for forecasting and covariance propagation. However, this does come at an increased computational cost. Namely, $Y=\mathrm{2}N+\mathrm{1}$ particles are required to capture the mean and covariance of the conditional state probability distribution. This implies that 2N+1 function evaluations are required for each UKF update. Even for a two-turbine wind farm in WFSim, a computational cost of 1×10² s per iteration ($k\to k+\mathrm{1}$) would not be surprising. While Eq. (14) can easily be parallelized, computational complexity remains troublesome, especially for larger wind farms. The issue of computational complexity is tackled by the ensemble KF.

3.5 The ensemble Kalman filter (EnKF)

The ensemble Kalman filter (EnKF; Evensen, 2003) is very similar to the UKF in that it relies on a finite number of realizations (the “sigma points” or “particles” in the UKF) to approximate the mean and covariance of the conditional probability distribution of x_k knowing 𝓩_k. However, whereas the UKF relies on a systematic way of distributing the particles such that the mean and covariance of the conditional probability distribution p[x_k|𝓩_k] are equal to that of the particles, the EnKF relies on random realizations, without guarantees that the mean and covariance are captured accurately. However, the EnKF has been shown to work well in a number of applications, with typically far fewer particles than states, i.e., Y≪N (Gillijns et al., 2006; Houtekamer and Mitchell, 2005). The forecast and update step are very similar to that of the UKF.

1. In the UKF the particles are redistributed at every time step, in contrast to the EnKF. Rather, the EnKF propagates the particles forward without redistribution. We define the ith particle as ${\mathit{\psi }}_{k|\mathrm{\ell }}^i\in {\mathbb{R}}^N$, which is a realization of the conditional probability distribution p[x_k|𝓩_ℓ]. The forecast steps are

   • i.

     Each particle is propagated forward in time using the nonlinear system dynamics, and with the realizations of noise terms w and v denoted by ${\widehat{\mathit{w}}}_{k-\mathrm{1}}^i\in {\mathbb{R}}^N$ and ${\widehat{\mathit{v}}}_k^i\in {\mathbb{R}}^M$, generated using MATLABs randn(…) function.

     $$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathit{\psi }}_{k|k-\mathrm{1}}^i=f({\mathit{\psi }}_{k-\mathrm{1}|k-\mathrm{1}}^i,{\mathit{q}}_{k-\mathrm{1}})+{\widehat{\mathit{w}}}_{k-\mathrm{1}}^i\\ {\displaystyle }& {\displaystyle }\text{for}i=\mathrm{1},\mathrm{\dots },Y,\\ {\displaystyle }& {\displaystyle }{\mathit{\zeta }}_{k|k-\mathrm{1}}^i=h({\mathit{\psi }}_{k|k-\mathrm{1}}^i,{\mathit{q}}_k)+{\widehat{\mathit{v}}}_k^i\\ \text{(19)}& {\displaystyle }& {\displaystyle }\text{for}i=\mathrm{1},\mathrm{\dots },Y.\end{array}$$
   • ii.

     The expected state and output are identically calculated as in the UKF using Eq. (15) with ${\mathit{w}}_{\text{mean}}^i={\left(Y-\mathrm{1}\right)}^{-\mathrm{1}}$. The covariance matrices are (re-)estimated from the forecasted ensemble by

     $$\begin{array}{ll}{\displaystyle }{\mathbf{P}}_{k|k-\mathrm{1}}^{\mathit{z}}={\displaystyle \frac{\mathrm{1}}{Y-\mathrm{1}}}\sum _{i=\mathrm{1}}^Y& {\displaystyle }(\left({\mathit{\zeta }}_{k|k-\mathrm{1}}^i-{\stackrel{\mathrm{‾}}{\mathit{\zeta }}}_{k|k-\mathrm{1}}\right)\\ \text{(20)}& {\displaystyle }& {\displaystyle }{\left({\mathit{\zeta }}_{k|k-\mathrm{1}}^i-{\stackrel{\mathrm{‾}}{\mathit{\zeta }}}_{k|k-\mathrm{1}}\right)}^{\mathrm{T}}),{\displaystyle }{\mathbf{P}}_{k|k-\mathrm{1}}^{\mathit{x}\mathit{z}}={\displaystyle \frac{\mathrm{1}}{Y-\mathrm{1}}}\sum _{i=\mathrm{1}}^Y& {\displaystyle }(\left({\mathit{\zeta }}_{k|k-\mathrm{1}}^i-{\stackrel{\mathrm{‾}}{\mathit{\zeta }}}_{k|k-\mathrm{1}}\right)\\ \text{(21)}& {\displaystyle }& {\displaystyle }{\left({\mathit{\psi }}_{k|k-\mathrm{1}}^i-{\stackrel{\mathrm{‾}}{\mathit{\psi }}}_{k|k-\mathrm{1}}\right)}^{\mathrm{T}}).\end{array}$$

2. For the analysis step, one applies Eq. (10) to determine the Kalman gain L_k. Then, each particle is updated individually as

   $$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathit{\psi }}_{k|k}^i={\mathit{\psi }}_{k|k-\mathrm{1}}^i+{\mathbf{L}}_k\left({\mathit{z}}_k-{\mathit{\zeta }}_{k|k-\mathrm{1}}^i\right)\\ \text{(22)}& {\displaystyle }& {\displaystyle }\text{for}i=\mathrm{1},\mathrm{\dots },Y.\end{array}$$

Note that, in contrast to the ExKF and the UKF, the state covariance matrix P^x (see Eqs. 7 and 12) need not be calculated explicitly in the EnKF. This, in combination with the small number of particles Y≪N, is what makes the EnKF computationally superior to the UKF (and often also computationally superior to the ExKF). However, this reduction in computational complexity comes at a price. The disadvantages of the EnKF are discussed in the next section.

Challenges in the EnKF for small number of particles

The caveat to representing the conditional state probability distribution with fewer particles than states, Y≪N, is the formation of inbreeding and long-range spurious correlations (Petrie, 2008). The former, inbreeding, is defined as a situation where the state error covariance matrix P^x is consistently underestimated, leading to state estimates that incorrectly rely more on the internal model. One straight-forward method to address this is called “covariance inflation”, in which P^x (or rather, the ensemble from which P^x is calculated) is “inflated” to correct for the underestimated state uncertainty (Petrie, 2008). Mathematically, this is achieved by applying

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathit{\psi }}_{k|k-\mathrm{1}}^i={\stackrel{\mathrm{‾}}{\mathit{\psi }}}_{k|k-\mathrm{1}}+r\left({\mathit{\psi }}_{k|k-\mathrm{1}}^i-{\stackrel{\mathrm{‾}}{\mathit{\psi }}}_{k|k-\mathrm{1}}\right)\\ \text{(23)}& {\displaystyle }& {\displaystyle }\text{for}i=\mathrm{1},\mathrm{\dots },Y\end{array}$$

before the analysis step, with r∈ℝ the inflation factor, typically with a value of 1.01–1.25.

The latter problem, long-range spurious correlations, can be better visualized in Fig. 2.

Figure 2Long-range spurious correlations arise in the case where a covariance matrix is described by a small number of particles. Using physical knowledge of the system, these undesired correlations can be corrected. Φ^x is the localization matrix. Applying localization, the covariance of physically nearby states are multiplied with a value close to 1, and the covariance of physically distant states are multiplied with a value close to 0. In our example case, this results in the localized covariance matrix Φ^x∘P^x, where ∘ is the element-wise product.

Download

In particle-based approaches, the covariance terms cannot be captured exactly. This may lead to the formation of small yet nonzero covariance terms between states and outputs which, in reality, are uncorrelated. This can lead to the drift of unobservable states, and eventually to instability of the KF. Increasing the number of particles is the most straight-forward solution to this problem, but comes at a huge computational cost. A better alternative is “covariance localization”, where physical knowledge of the states and measurements is used to steer the sample-based covariance matrices. Recall that in the surrogate model of Sect. 2, the model states are the velocity and pressure terms inside the wind farm at a physical location. Define that the ith state entry (x_k)_i belongs to a physical location in the farm s_i. Then, looking at an arbitrary state covariance term (i,j),

$${\left({\mathbf{P}}_{k|k-\mathrm{1}}^{\mathit{x}}\right)}_{i,j}=E\left[\left(({\mathit{x}}_k{)}_i-({\widehat{\mathit{x}}}_{k|k-\mathrm{1}}{)}_i\right){\left(({\mathit{x}}_k{)}_j-({\widehat{\mathit{x}}}_{k|k-\mathrm{1}}{)}_j\right)}^{\mathrm{T}}\right],$$

we define the physical distance between these two states as $\mathrm{\Delta }{s}_{i,j}=\left|\right|{\mathit{s}}_i-{\mathit{s}}_j|{|}_{\mathrm{2}}.$ Now, we introduce a weighting factor into our covariance matrices by multiplying physically distant states with a value close to 0, and multiplying physically nearby states with a value close to 1. A popular choice for such a weighting function is Gaspari–Cohn's fifth-order discretization of a Gaussian distribution (Gaspari and Cohn, 1999), given by

$$\begin{array}{}\text{(24)}& {\displaystyle }\mathit{\varphi }\left(c_{i,j}\right)=\left\{\begin{array}{l}-\frac{\mathrm{1}}{\mathrm{4}}{c}_{i,j}^{\mathrm{5}}+\frac{\mathrm{1}}{\mathrm{2}}{c}_{i,j}^{\mathrm{4}}+\frac{\mathrm{5}}{\mathrm{8}}{c}_{i,j}^{\mathrm{3}}-\frac{\mathrm{5}}{\mathrm{3}}{c}_{i,j}^{\mathrm{2}}+\mathrm{1}\\ \text{if }\mathrm{0}\le c_{i,j}\le \mathrm{1},\\ \frac{\mathrm{1}}{\mathrm{12}}{c}_{i,j}^{\mathrm{5}}-\frac{\mathrm{1}}{\mathrm{2}}{c}_{i,j}^{\mathrm{4}}+\frac{\mathrm{5}}{\mathrm{8}}{c}_{i,j}^{\mathrm{3}}+\frac{\mathrm{5}}{\mathrm{3}}{c}_{i,j}^{\mathrm{2}}\\ -\mathrm{5}{c}_{i,j}+\mathrm{4}-\frac{\mathrm{2}}{\mathrm{3}}\frac{\mathrm{1}}{c_{i,j}}\\ \text{if }\mathrm{1}<c_{i,j}\le \mathrm{2},\\ \mathrm{0}\\ \text{otherwise,}\end{array}\right.\end{array}$$

with $c_{i,j}=\frac{\left|\right|\mathrm{\Delta }{\mathit{s}}_{i,j}|{|}_{\mathrm{2}}}{L}$ a normalized distance measure, with L the cut-off distance. Applying Eq. (24) for the covariance matrices ${\mathbf{P}}_{k|k-\mathrm{1}}^{\mathit{z}}$ and ${\mathbf{P}}_{k|k-\mathrm{1}}^{\mathit{x}\mathit{z}}$ we can define the localization matrices

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathbf{\Phi }}^z=\left[\begin{array}{ccc}\mathit{\varphi }\left(c_{\mathrm{1},\mathrm{1}}^z\right)& \mathrm{\cdots }& \mathrm{\cdots }\mathit{\varphi }\left(c_{\mathrm{1},M}^z\right)\\ \mathrm{⋮}& \mathrm{\ddots }& \\ \mathit{\varphi }\left(c_{M,\mathrm{1}}^z\right)& & \mathit{\varphi }\left(c_{M,M}^z\right)\end{array}\right],\\ {\displaystyle }& {\displaystyle }{\mathbf{\Phi }}^{xz}=\left[\begin{array}{ccc}\mathit{\varphi }\left(c_{\mathrm{1},\mathrm{1}}^{xz}\right)& \mathrm{\cdots }& \mathrm{\cdots }\mathit{\varphi }\left(c_{\mathrm{1},M}^{xz}\right)\\ \mathrm{⋮}& \mathrm{\ddots }& \\ \mathit{\varphi }\left(c_{N,\mathrm{1}}^{xz}\right)& & \mathit{\varphi }\left(c_{N,M}^{xz}\right)\end{array}\right],\end{array}$$

where $c_{i,j}^z$ is the normalized distance between two measurements i and j, and $c_{i,j}^{xz}$ is the normalized distance between state i and measurement j, respectively. Finally, localization and inflation can be incorporated into Eqs. (20) and (21) by

$$\begin{array}{ll}{\displaystyle }{\mathbf{P}}_{k|k-\mathrm{1}}^{\mathit{z}}& {\displaystyle }={\mathbf{\Phi }}^z\circ {\displaystyle \frac{\mathrm{1}}{Y-\mathrm{1}}}\sum _{i=\mathrm{1}}^Y(\left({\mathit{\zeta }}_{k|k-\mathrm{1}}^i-{\stackrel{\mathrm{‾}}{\mathit{\zeta }}}_{k|k-\mathrm{1}}\right)\\ \text{(25)}& {\displaystyle }& {\displaystyle }{\left({\mathit{\zeta }}_{k|k-\mathrm{1}}^i-{\stackrel{\mathrm{‾}}{\mathit{\zeta }}}_{k|k-\mathrm{1}}\right)}^{\mathrm{T}}),{\displaystyle }{\mathbf{P}}_{k|k-\mathrm{1}}^{\mathit{x}\mathit{z}}& {\displaystyle }=r\cdot {\mathbf{\Phi }}^{xz}\circ {\displaystyle \frac{\mathrm{1}}{Y-\mathrm{1}}}\sum _{i=\mathrm{1}}^Y(\left({\mathit{\zeta }}_{k|k-\mathrm{1}}^i-{\stackrel{\mathrm{‾}}{\mathit{\zeta }}}_{k|k-\mathrm{1}}\right)\\ \text{(26)}& {\displaystyle }& {\displaystyle }{\left({\mathit{\psi }}_{k|k-\mathrm{1}}^i-{\stackrel{\mathrm{‾}}{\mathit{\psi }}}_{k|k-\mathrm{1}}\right)}^{\mathrm{T}}),\end{array}$$

where ∘ is the element-wise product (Hadamard) of the two matrices. The improvement in terms of computational efficiency and estimation performance is displayed in Fig. 3.

Figure 3This figure shows the estimation performance and computational cost (parallelized, 8 cores) of the EnKF for a range of ensemble sizes, with and without inflation and localization. Great improvement is seen for estimation accuracy, at no additional computational cost. The simulation scenario is described in detail in Sect. 4.2, and the results presented here are rather meant as an indication.

Download

A significant increase in performance is shown, especially for smaller numbers of particles. This is in agreement with what was seen in previous work (Doekemeijer et al., 2017). Furthermore, performance is more consistent. Additionally, note that there is no increase in computational cost, as the covariance matrices are made sparse, leading to a cost reduction in the calculation of Eq. (10), which makes up for the extra operations of Eqs. (25) and (26). Also, note that the localization matrices are time-invariant and can be calculated offline.

3.6 Synthesizing an online model calibration solution

Certain model parameters such as ℓ_s are closely related to the turbulence intensity, which vary over time. Estimation of such parameters is achieved by extending the state vector with (a subset of) the model parameters. In this work, ℓ_s is concatenated to the state vector as random walk model, with a certain standard deviation (covariance). Higher values of ℓ_s lead to more wake recovery, making the calibration solution adaptable to varying turbulence levels. This adds one scalar entry to x_k, which is a negligible addition in terms of computational cost.

Furthermore, a proposal is made for the estimation of the free-stream wind speed U_∞. This is suggested to be done using the turbine's power generation measurements, following the ideas of Gebraad et al. (2016) and Shapiro et al. (2017b). Using the wind vanes and employing a simple steady-state wake model from the literature (Mittelmeier et al., 2017), the turbines operating in free-stream flow can be distinguished from the ones operating in waked flow. Next, define Γ∈ℤ^ℵ as a vector specifying the upstream turbines, with ℵ the total number of turbines operating in free stream. Then, the instantaneous rotor-averaged flow speed at each turbine's hub can be estimated by inverting the turbine power expression from WFSim (Boersma et al., 2018). One wind-farm-wide free-stream wind speed U_∞ is then calculated using actuator disk theory. Smoothing results with a low-pass filter with time-constant $c_{u_{\mathrm{\infty }}}$ on the average of $U_{{\mathrm{\infty }}_i}$ for each upstream turbine i, we obtain

where it is assumed that $U_{{\mathrm{\infty }}_i}\approx U_{r_i}\left(\mathrm{1}+\frac{\mathrm{1}}{\mathrm{4}}\cdot C_{T_i}^{\prime }\right)$ when γ_i≈0, with $U_{r_i}$ the wind speed at the rotor of turbine i. Research is currently ongoing on how to best incorporate the effects of turbine yaw (γ≠0) into the definition of $C_T^{\prime }$. Furthermore, ρ is the air density, A is the rotor swept area, and $P_{\text{turb},i}^{\text{meas.}}$ is the measured instantaneous power capture of turbine i².

Combining these elements yields an efficient, modular, and accurate model calibration solution for WFSim. The model states are estimated using SCADA and/or lidar data, of which the former is readily available, and the latter becoming more popular. State estimation paired with parameter estimation improves the accuracy of the surrogate model, potentially leading to more accurate control. Additionally, the free-stream wind speed is estimated using readily available SCADA data. This control solution is implemented in MATLAB, and leverages the numerically efficient pre-compiled solvers and parallelization for model propagation. The EnKF is orders of magnitude faster than existing estimation algorithms due to covariance localization and inflation, while competing with the UKF in terms of accuracy.

4.1 SOWFA

High-fidelity simulation data are generated using the Simulator fOr Wind Farm Applications (SOWFA), developed by the National Renewable Energy Laboratory (NREL). SOWFA provides accurate flow data at a fraction of the cost of field tests. It solves the filtered, three-dimensional, unsteady, incompressible Navier–Stokes equations over a finite temporal and spatial mesh, accounting for the Coriolis and geostrophic forcing terms. SOWFA is a LES solver, meaning that larger-scale dynamics are resolved directly, and turbulent structures smaller than the discretization are approximated using subgrid-scale models to suppress computational cost. (Churchfield et al., 2012). The turbine rotor is modeled using an actuator line representation (ALM) as derived from Sorensen and Shen (2002). SOWFA has previously been used for lower-fidelity model validation, controller testing, and to study the aerodynamics in wind farms (Fleming et al., 2016, 2017a; Gebraad et al., 2017). The interested reader is referred to Churchfield et al. (2012) for a more in-depth description of SOWFA and LES solvers in general.

Table 2Covariance settings for the KF variants, with I_• the ℝ^•×• identity matrix. The full cov. matrices are diagonal concatenations of the entries. For example, P₀ is $\mathbf{\text{diag}}\left({\mathbf{P}}_{\mathrm{0},u},{\mathbf{P}}_{\mathrm{0},v}\right)$ and $\mathbf{\text{diag}}\left({\mathbf{P}}_{\mathrm{0},u},{\mathbf{P}}_{\mathrm{0},v},{\mathbf{P}}_{\mathrm{0},{\mathrm{\ell }}_{\mathrm{s}}}\right)$ for state-only and state-parameter estimation, respectively.

Download Print Version | Download XLSX

4.2 Two-turbine simulation with turbulent inflow

In this section, a two-turbine wind farm is simulated to analyze the effect of different measurement sources, KF algorithms, and the difference between state-only and state-parameter estimation. This simple wind farm contains two NREL 5-MW baseline turbines with D=126.4 m, separated 5D in stream-wise direction. This LES simulation was described in more detail in Annoni et al. (2016). Important simulation properties are listed in Table 1 for SOWFA and WFSim. The effect of the turbulence intensity on the wake dynamics in SOWFA is captured in WFSim through its mixing-length turbulence model. In these simulations, WFSim is purposely initialized with a too low value for ℓ_s in order to represent the realistic situation of a model mismatch. The remaining tuning parameters in WFSim were chosen such that a weighted-sum cost function of the power and flow errors was minimized.

Figure 5Comparison of absolute values of the estimation errors (in long. flow fields) for state-only estimation with the EnKF for various sensor configurations: using turbine power measurements, using flow measurements with a lidar system pointing upstream, and using flow measurements with a lidar system pointing downstream of the rotor. Here, $(\mathrm{\Delta }\mathit{u}{)}_{\mathrm{•}}=|{\mathit{u}}_{\mathrm{•}}-{\mathit{u}}_{\text{SOWFA}}|$. Here, wind flows from top to bottom. The sensors are depicted by red dots (flow meas.) or red turbines (power meas.), not to be confused with estimation error.

Download

Firstly, the three KF variants will be compared in Sect. 4.2.1. Secondly, in Sect. 4.2.2, estimation using different information sources is compared. Thirdly, the potential of joint state-parameter estimation is displayed in Sect. 4.2.3.

4.2.1 A comparison of the KF variants for state estimation

In this simulation study, four estimation cases are compared: (1) the ExKF, (2) the UKF, (3) the EnKF, and (4) the open-loop (OL) simulation, i.e., without estimation. The focus here is on state-only estimation, thus excluding ℓ_s. Flow measurements downstream of each turbine are assumed (e.g., using lidar), their locations denoted as red dots in Fig. 4, which is about 2 % of the full to-be-estimated state space. These measurements are artificially disturbed by zero-mean white noise with σ=0.10 m s⁻¹. The KF settings are listed in Tables 2 and 3. The KF covariance matrices were obtained through an iterative tuning process in previous work (Doekemeijer et al., 2017) with minor adjustments to simulate performance for untrained data. Figure 4 shows state (flow field) estimation of the three KF variants for two time instants, t=300 s and t=700 s. In this figure, $(\mathrm{\Delta }\mathit{u}{)}_{\mathrm{•}}\in {\mathbb{R}}^{N_u}$ is defined as the absolute error between the estimated and true longitudinal flow velocities in the field.

Looking at Fig. 4, the OL estimations are accurate for the unwaked and single-waked flow, yet are lacking in the situation of two overlapping wakes, for which the KFs correct. There is no significant difference in accuracy between the different KF variants, yet they differ by 2 orders of magnitude in computational cost (Table 3).

Table 3Choice of tuning parameters for the KF variants, for both the two-turbine and nine-turbine simulation case. Note that the ExKF does not support power measurements nor parameter estimation due to the lack of linearization, and does not have any additional tuning parameters. In terms of computational cost: simulations were run on a single node using 8 cores in parallel.

Download Print Version | Download XLSX

Table 4Overview of several settings for the SOWFA and the WFSim nine-turbine wind farm simulation.

Download Print Version | Download XLSX

Figure 6Comparison of forecasting performance for state-only and joint state-parameter (ℓ_s) estimation with the EnKF and UKF, where measurements are available up until the vertical dashed lines, after which the estimation becomes a forecast. Here, the 2-norm of the estimation error is plotted along the y axis, with $(\mathrm{\Delta }\mathit{u}{)}_{\mathrm{•}}=|{\mathit{u}}_{\mathrm{•}}-{\mathit{u}}_{\text{SOWFA}}|$.

Download

4.2.3 Joint state-parameter estimation

Forecasting, as used in predictive control, benefits from the calibration of model parameters in addition to the states. Joint state-parameter estimation using flow measurements downstream of each turbine (as shown in the rightmost plots in Fig. 5) disturbed by zero-mean white noise with σ=0.10 m s⁻¹ for the EnKF and UKF is displayed in Fig. 6. The KF settings are shown in Tables 2 and 3. At t=0 s, both the OL and the KF simulations start with the same (wrong) value for ℓ_s. Then, every second, (noisy) measurements are fed into the KFs, and the state vector as well as the model parameter ℓ_s is estimated. However, for the OL simulation, no measurements are fed in: the state vector is simply updated with the nominal model, and the value for ℓ_s remains the same throughout the simulation. Now, after 600 s (left plot in Fig. 6) and 900 s (right plot in Fig. 6), a forecast is started, meaning no measurements are available after that time. At that moment, the OL model still has the same (poor) value for ℓ_s as at t=0 s, while the value for ℓ_s in the KFs has improved. From Fig. 6, it becomes clear that the estimates are not only improved for the 3 min forecast, but are also consistently better than the non-calibrated (open-loop) model's 10 min forecast due to the estimation of ℓ_s³. Furthermore, the EnKF performs comparably to the UKF at a lower computational cost. Note that the EnKF even outperforms the UKF in this simulation, expected to be due to randomness in the EnKF.

Figure 7Convergence of ℓ_s and U_∞ using the EnKF. Dashed lines are the grid-searched optimal constant values for the open-loop simulation. With power measurements only, the EnKF is able to estimate these parameters successfully in addition to the model states.

Download

Figure 8Comparison of absolute values of the estimation errors (in long. flow fields) for state-parameter estimation with the EnKF. Wind is coming in from the top and flows downwards. The variables U_∞ and ℓ_s are incorrectly initialized in both the OL and the EnKF. In the EnKF, U_∞ and ℓ_s are estimated in addition to the states, using only turbine power measurements. The EnKF quickly converges for the states, and more slowly for ℓ_s and U_∞. After 300 s, the EnKF has converged to a negligible estimation error.

Download

4.3 Nine-turbine simulation with turbulent inflow

In this section, we investigate the performance of the EnKF-based model calibration solution under a more realistic nine-turbine wind farm scenario. The purpose of this case study is to highlight the need for state-parameter estimation for accurate wind farm modeling. The wind farm contains nine NREL 5-MW baseline turbines, oriented in a three by three layout, separated 5D and 3D in stream- and crosswise directions, respectively. The turbines start with a 30^∘ yaw misalignment, but are then aligned with the mean wind direction within the first 30 s of simulation. The turbine layout and numbering is shown in the top-left panel of Fig. 8. This LES simulation has been used before in the literature, and is described in more detail in Boersma et al. (2018). A number of important simulation properties are listed in Table 4 for SOWFA and WFSim, respectively.

Figure 9Comparison of power forecasting using the EnKF with measurements available up until time t=600 s. After convergence U_∞ (as seen as a positive power slope for the first row of turbines), ℓ_s is also calibrated. After convergence, forecasting is better than in open-loop. Oscillatory behavior is still present due to an oscillatory input signal ($C_T^{\prime }$), turbulent flow field, and the absence of inertia in the rotor model.

Download

Table 5Turbine-averaged RMSE in power time series of Fig. 9 (compared to SOWFA). The lower the RMSE, the better the forecast.

Download Print Version | Download XLSX

Compared to the two-turbine case, N has increased by a factor of 4. In the UKF, this would result in the same factor of additional particles. Thus, not only is each particle more expensive to calculate but there are also more particles. Rather, in the EnKF, the approach is heuristic. None of the EnKF settings needed to be changed for good performance compared to Sect. 4.2, as displayed in Tables 2 and 3.

As shown in Table 3, the EnKF has a low computational cost of 1.2 s per iteration (8 cores, parallel). In this case study, both the complete model state (flow field), the turbulence model parameter ℓ_s, and the free-stream flow speed U_∞ are estimated in real-time using exclusively (readily available) power measurements from the turbines. The EnKF and one of the open-loop simulations (OL) will deliberately be initialized with poor values for ℓ_s and U_∞ to investigate convergence. The other OL simulation will be initialized with a poor value for ℓ_s but a correct value for U_∞ for comparison.

In Fig. 7, it can be seen that the EnKF is successful in estimating U_∞ and ℓ_s after 300 s using only wind turbine power measurements. Furthermore, the flow fields of SOWFA, of the OL simulation with U_∞=9.0 m s⁻¹, and of the EnKF at various time instants are displayed in Fig. 8. From this figure, it can be seen that the EnKF has large errors at the start of the simulation. However, after 10 s, the error in flow states surrounding each turbine significantly decreases through the use of turbine power measurements. This estimated flow then propagates downstream, “clearing up” the errors in the vicinity of the wind turbines. As time further propagates, the free-stream estimation improves, and finally the estimation error converges.

The power forecasting performance is shown in Fig. 9 and Table 5. As also seen in Fig. 7, the EnKF converges after 300 s, and indeed the power forecasts outperform those of the OL simulation at t=300 s. Furthermore, it is interesting to see that the filtered power estimates of the first row of turbines (i=1, 2, 3) starts low at t=1 s, but converges to the true power at t≈200 s. This can be related to the mismatch in U_∞, which takes approximately 300 s to converge to the true value of 12 m s⁻¹, as seen in Fig. 7. The oscillatory behavior in both the OL and EnKF power predictions is due to the absence of rotor inertia in the rotor model, turbulent structures in the flow, and large fluctuations on the excitation signal $C_T^{\prime }$.

Finally, the forecasts for flow at times t=300 s and t=600 s are examined in Fig. 10. The large flow estimation mismatch in the EnKF at t<250 s quickly reduces and for t≥250 s the EnKF estimation is consistently better than both the OL cases. This has to do with the convergence of the model parameters ℓ_s and U_∞, and the estimation of the states surrounding the turbines using the power measurements.

Figure 10Comparison of flow field estimation for the nine-turbine case. Measurements are available until t=300 s (a) and t=600 s (b), respectively. The EnKF converges to the true U_∞ after 300 s. After convergence, the forecasts are significantly better than in open-loop simulations.

Download

A crucial remark with the simulations presented here is that low-frequency changes in the atmosphere are neglected. In a real wind farm, atmospheric properties such as the mean wind direction and turbulence intensity change continuously, and this will impact the estimation and forecasting performance. The EnKF uses an assumption of persistence for the atmospheric properties at the time of forecasting, and thus a change in mean wind direction may invalidate the model forecast. In future work, the algorithm presented here should be tested under high-fidelity simulations with such realistic low-frequency changes. This would provide insight into the potential of the work at hand, and advance towards a practical wind farm implementation.