Fourier series-based approximation of time-varying parameters in ordinary differential equations

The approach presented in this work focuses on recovering unknown and unmeasured dynamical system parameters with time-dependent behavior, updating the parameterization of the physical model through use of a sinusoidal representation of the underlying TVP. In essence, this method utilizes a dictionary of sinusoidal functions inspired by the Fourier series terms and weighted with coefficients estimated online by ensemble Kalman filtering to sequentially update the TVP approximation. The main contributions of this work include:

• A new method for estimating TVPs in dynamical systems, wherein sinusoidal functions inspired by Fourier series terms are used in a dictionary to generate an approximation model for the TVP, thereby transforming the TVP estimation into a constant coefficient estimation inverse problem;
• Two implementation options for sequentially updating the TVP approximations, both using an ensemble Kalman filtering framework to track the system states and estimate the unknown TVP approximation model coefficients; and
• Strategies for using this method to estimate periodic parameters both when the period is known and unknown, as well as extensions to estimate TVPs that are not periodic over the time interval of available data.

In a set of numerical examples considering a forced mass-spring system, we first establish the ability of the proposed Fourier series-based approximation approach to estimate a periodic forcing parameter with known period, then extend the estimation to include the period of the TVP as an additional unknown. We further demonstrate the capability of the proposed method in estimating TVPs that are not periodic (and potentially not continuous) over the time interval of observed data and show that the resulting parameter approximation models can be used to make reasonably accurate predictions of the system dynamics for different initial conditions. In the appendix, we include additional numerical results comparing the proposed approach with a standard augmented ensemble Kalman filtering algorithm and parameter tracking techniques, as well as an example considering a periodically-forced, partially observed Lorenz system with chaotic behavior.

We remark that the focus of this work does not immediately lie on application to large-scale problems using small ensemble sizes or on systems with chaotic dynamics (as is common in the data assimilation literature) but rather on establishing a method to well estimate the functional form (shape and magnitude) of the underlying TVP. Even for small systems (e.g. where d = 2 or 3), the number of unknown coefficients grows with the number of terms used in the TVP approximation model, thereby requiring use of an ensemble size larger than the number of system states. We discuss considerations including the choice of the number of approximation model terms needed to well represent the underlying TVP in this work.

The remainder of the paper is organized as follows: section 2 briefly reviews ensemble Kalman filtering for constant parameter estimation. Section 3 details the proposed Fourier series-based TVP approximation models, providing two implementation approaches that can be used together with ensemble Kalman filtering to estimate the unknown model coefficients. Section 4 provides the main results of the numerical experiments (with some additional results provided in the appendix), and section 5 gives conclusions and future work.

Recall that the Qth order Fourier series of a univariate function h(t) is given by

$$\begin{align} h\left(t\right) = \frac{a_0}{2} + \sum_{q = 1}^{Q} \left( a_q \cos\left(\frac{2\pi qt}{P}\right) + b_q\sin\left(\frac{2\pi qt}{P}\right) \right) \end{align} \tag{ 4 }$$

where a_q and b_q are the expansion coefficients and P is the period. If h(t) is known, the coefficients a_q and b_q can be computed explicitly using the formulas

$$\begin{align} a_q = \displaystyle\frac{2}{P} \displaystyle\int_0^P h\left(t\right) \cos\left(\frac{2\pi qt}{P}\right) \textrm{d}t \end{align} \tag{ 5 }$$

and

$$\begin{align} b_q = \displaystyle\frac{2}{P} \displaystyle\int_0^P h\left(t\right) \sin\left(\frac{2\pi qt}{P}\right) \textrm{d}t, \end{align} \tag{ 6 }$$

respectively. For example, figure 2 shows the Fourier series approximations of the periodic function $h(t) = 2\sin(t) - 0.5\cos(2t/3)$ with known period $P = 6\pi$ for different choices of Q. However, when using Fourier series for function approximation, the function h(t) is generally unknown and the coefficients a_q and b_q must be estimated [15, 18, 20].

Zoom In Zoom Out Reset image size

Inspired by use of the Fourier series for approximating unknown functions, we propose to represent the unknown, univariate TVPs $\theta(t)$ in this work as linear combinations of sine and cosine pairs, such that the estimate of each $\theta(t)$ is determined by the approximation model

$$\begin{align} \theta_M\left(t\right) = \alpha_0 + \sum_{i = 1}^M \left( \alpha_i \sin\left(\omega_i t\right) + \beta_i \cos\left(\omega_i t\right) \right), \end{align} \tag{ 7 }$$

with coefficients α_i , β_i and fixed values of ω_i for each $i = 1,\dots,M$. Moving forward, we denote the $2M+1$ coefficients by c_k , $k = 0,\dots,2M$, such that the model in (7) becomes

$$\begin{align} \theta_M\left(t\right) = c_0 + c_1 \sin\left(\omega_1 t\right) + c_2 \cos\left(\omega_1 t\right) + \cdots + c_{2M-1} \sin\left(\omega_M t\right) + c_{2M} \cos\left(\omega_M t\right), \end{align} \tag{ 8 }$$

for some fixed M. Our goal in approximating $\theta(t)$ is therefore to estimate the $2M+1$ coefficients c_k that provide the best fit between the model in (8) and the true time-varying parameter given the observed system data. This essentially transforms the inverse problem at hand into a constant parameter estimation problem, and we can utilize the EnKF algorithm described in section 2 to estimate these coefficients. We remark that this problem is similar in spirit to dictionary learning, where writing the TVP approximation model in (8) equivalently as

$$\begin{align} \theta_M\left(t\right) = \left[ \begin{array}{cccccc} 1 & \sin\left(\omega_1 t\right) & \cos\left(\omega_1 t\right) & \cdots & \sin\left(\omega_M t\right) & \cos\left(\omega_M t\right) \end{array}\right] c \end{align} \tag{ 9 }$$

highlights the dictionary of sinusoidal functions, inspired by the Fourier series terms, and $c\in\mathbb{R}^{2M+1}$ is a vector containing the unknown coefficients c_k . The coefficient estimation thereby determines which of these functions make meaningful contributions in representing the underlying shape and magnitude of $\theta(t)$ given the available system observations.

To apply the TVP approximation model in (8), one must specify the value of M, which sets the number of terms in the approximation, and the values of ω_i , $i = 1,\dots,M$, which act as the angular frequencies of the sinusoidal functions. When $\theta(t)$ is periodic, we consider two different approaches in assigning the ω_i values: If we know the period of the underlying $\theta(t)$ in advance of the estimation process, we explicitly define ω_i using the formula

$$\begin{align} \omega_i = \frac{2\pi i}{P}, \quad i = 1,\dots, M \end{align} \tag{ 10 }$$

where P is the known period of $\theta(t)$; this form of ω_i follows from the terms in the Fourier series in (4) for approximating a periodic function with known period. If the period of $\theta(t)$ is unknown or uncertain, we set ω_i as in (10) and treat the period P as an additional unknown constant parameter to be jointly estimated along with the coefficients of the TVP approximation model. In more general cases, i.e. when $\theta(t)$ is not known to be periodic over the time interval of available data, we instead choose a fixed increment ω and define $\omega_i = \omega i$ for each $i = 1,\dots, M$, systematically incorporating sine and cosine pairs with different frequencies into the approximation. Following this approach, appropriate choice of the increment ω becomes an important factor in the estimation process.

In the previous section, we propose an approximation model for $\theta(t)$ that we can use directly within the ODE model in (1) in place of the TVP. As an alternative means of implementation, much in the spirit of state augmentation, we can use the approximation model in (8) to prescribe a model for $\textrm{d}\theta/\textrm{d}t$ and define a coupled ODE system of the form

$$\begin{eqnarray} \frac{\textrm{d}x}{\textrm{d}t} & = & f\left(t,x,\theta\right), \quad x\left(0\right) = x_0 \end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray} \frac{\textrm{d}\theta}{\textrm{d}t} & = & \theta^{^{\prime}}_M\left(t\right), \quad \theta\left(0\right) = \theta_0, \end{eqnarray} \tag{ 12 }$$

which we can write equivalently as the augmented system

$$\begin{align} \displaystyle\frac{dz}{\textrm{d}t} = \left[ \begin{array}{c} \displaystyle\frac{\textrm{d}x}{\textrm{d}t} \\[.4cm] \displaystyle\frac{\textrm{d}\theta}{\textrm{d}t} \end{array} \right] = \left[ \begin{array}{c} f\left(t,x,\theta\right) \\[.2cm] \theta^{^{\prime}}_M\left(t\right) \end{array} \right] , \quad z\left(0\right) = z_0 = \left[ \begin{array}{c} x_0 \\[.2cm] \theta_0 \end{array} \right] \end{align} \tag{ 13 }$$

where $z = z(t) \in \mathbb{R}^{d+1}$. In this representation, $\theta = \theta(t)$ is treated as an unobserved state of the system in (13). The equation for $\textrm{d}\theta/\textrm{d}t$ in (13) follows from taking the derivative of the approximation model $\theta_M(t)$ in (8), which gives

$$\begin{align} \theta^{^{\prime}}_M\left(t\right) = c_1 \omega_1 \cos\left(\omega_1 t\right) - c_2 \omega_1 \sin\left(\omega_1 t\right) + \cdots + c_{2M-1} \omega_M \cos\left(\omega_M t\right) - c_{2M} \omega_M \sin\left(\omega_M t\right), \end{align} \tag{ 14 }$$

with 2M unknown coefficients, $c_1, \dots, c_{2M}$. Note that the additive constant c₀ no longer explicitly appears in the system equations; instead, we estimate the initial condition θ₀ as an additional constant parameter playing a similar role.

In this example, we consider a sinusoidal forcing parameter of the form $\theta(t) = 2\sin(t) - 0.5\cos(2t/3)$ with period $P = 6\pi$ as the underlying truth. Utilizing MATLAB's ode45 to solve the ODE system in (18), we record observations every 0.5 time units over the interval [0,60], spanning just over three periods of $\theta(t)$, and corrupt the observations using Gaussian noise with zero mean and standard deviation taken to be 20% of the standard deviation of the true system states. Figure 3 shows the simulated data. In the experiments that follow, we consider two cases for the estimation procedure: (i) when the period of $\theta(t)$ is known, and (ii) when the period of $\theta(t)$ is unknown. Results focus on use of the TVP approximation model approach described in section 3.1. Appendix A.1 shows an example of the corresponding numerical results when using the derivative-based augmentation approach described in section 3.2 for the known period case.

Zoom In Zoom Out Reset image size

4.1.1. Estimation with Known Period.

Assuming a known period of $P = 6\pi$ for $\theta(t)$, we apply the formula in (10) to set $\omega_i = {i}/{3}$, $i = 1,\dots, M$, and employ the TVP approximation model in (8) with increasing values of M, ranging from M = 1 (which involves three coefficients) through M = 5 (11 coefficients). In each case, the EnKF tracks the model states and estimates the unknown coefficients comprising the Fourier series-based model approximation of $\theta(t)$. As an example, figure 4 shows the resulting EnKF time series estimates of position, velocity, and the seven TVP model coefficients when M = 3. We then use the posterior sample mean of each coefficient, $\bar{c}_k$, $k = 0,\dots,2M$, to construct an approximation of $\theta(t)$, such that

$$\begin{align} \bar{\theta}_{M}\left(t\right) = \bar{c}_0 + \bar{c}_1 \sin\left(\omega_1 t\right) + \bar{c}_2 \cos\left(\omega_1 t\right) + \cdots + \bar{c}_{2M-1} \sin\left(\omega_M t\right) + \bar{c}_{2M} \cos\left(\omega_M t\right) \end{align} \tag{ 19 }$$

for each M. Table 1 lists the posterior sample mean of each coefficient for each M value, and figure 5 shows the resulting $\bar{\theta}_{M}(t)$ approximations in each case compared with the true $\theta(t)$. To further compare the approximations with the true $\theta(t)$, we use a scaled version of the root mean square error (RMSE), where

$$\begin{align} \text{Scaled RMSE} = \frac{\text{RMSE}}{\sigma_\theta}, \qquad \text{RMSE} = \sqrt{\frac{1}{T}\sum_{j = 1}^T \left( \theta\left(t_j\right) - \bar{\theta}_{M}\left(t_j\right) \right)^2}, \end{align} \tag{ 20 }$$

and σ_θ is the standard deviation of $\theta(t)$. Note that the plots in figure 5 and corresponding scaled RMSE values in table 1 were computed at t_j values taken every 0.1 time units over [0,60], a finer time discretization than used during the filtering process.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1. EnKF posterior sample means of the TVP model coefficients for different values of M when approximating the sinusoidal forcing parameter $\theta(t) = 2\sin(t) - 0.5\cos(2t/3)$ in (18), using the data in figure 3 and assuming a known period for $\theta(t)$. Numbers are rounded to four decimal places. Note that the Fourier series coefficients obtained using the formulas in (5) and (6) to approximate the known function for orders 3 and above correspond to $c_4 = -0.5$, $c_5 = 2$, and $c_k = 0$ for all other coefficients.

Coefficient	M = 1	M = 2	M = 3	M = 4	M = 5
c0	−0.1328	−0.0306	0.0303	−0.0083	−0.0552
c1	−0.2106	−0.3373	−0.0013	−0.0000	0.0255
c2	0.3387	−0.2012	−0.0927	−0.0624	−0.0784
c3	—	0.1689	0.0545	0.0195	−0.0035
c4	—	−0.2770	−0.5236	−0.5032	−0.4689
c5	—	—	2.0299	2.0047	2.0538
c6	—	—	0.0594	−0.0032	−0.0533
c7	—	—	—	0.1099	0.1303
c8	—	—	—	0.0431	0.0053
c9	—	—	—	—	−0.1584
c10	—	—	—	—	0.1012
Scaled RMSE	1.0211	1.0087	0.0645	0.0657	0.1282

As illustrated in figure 4, the filter well tracks the model states for both p(t) and v(t), and the coefficient estimates converge to constant values after assimilating approximately one period of data (here, one period occurs at time $6\pi \approx 18.85$), with the ±2 standard deviation curves representing uncertainty around the mean estimate shrinking significantly. In comparing the results for different values of M, we note that the lowest RMSE occurs when M = 3 and is similar when M = 4, with a small increase in error when M = 5; however, all three of these M values result in reasonably close model approximations to the true $\theta(t)$, as shown in figure 5. Further, the EnKF posterior sample means for the TVP approximation model coefficients when M = 3, 4, and 5 are quite similar to the Fourier series coefficients obtained when using the formulas in (5) and (6) to approximate the function $h(t) = 2\sin(t) - 0.5\cos(2t/3)$ when Q = 3, 4, and 5, assuming that both h(t) and $P = 6\pi$ are known (i.e. $c_4 = -0.5$, $c_5 = 2$, and $c_k = 0$ for all other coefficients).

Figure 6 shows the mass-spring system model predictions using $\bar{\theta}_{M}(t)$ as the forcing parameter in (18) for each M with the initial condition $x(0) = [1; -1]$, a different initial condition than used in generating the simulated data, along with the corresponding scaled RMSE values for each model state. As might be expected given the TVP approximation results, the model predictions when using M = 1 and 2 are less accurate than when using M = 3, 4, and 5, which all provide similarly accurate predictions of both position and velocity. While not shown, similar results hold for other initial conditions in this range. For predictions using initial conditions larger in magnitude with this system, e.g. $x(0) = [100; -100]$, a similar pattern holds where the approximation models with M = 3, 4, and 5 provide the best results (i.e. lowest RMSE values), but the overall error is lower for predictions using any of the M values considered.

Zoom In Zoom Out Reset image size

4.1.2. Estimation with unknown period.

Given the same simulated data, we now assume that the period P of $\theta(t)$ is also unknown and estimate P along with the unknown approximation model coefficients for different choices of M. We draw an initial sample of P values from a uniform distribution over [15, 20], supposing that we have some reasonable prior information on a range of likely values (the true period being $6\pi\approx 18.8496$ in this case). As an example for comparison, figure 7 shows the resulting EnKF time series estimates of the seven TVP model coefficients and period P when M = 3. Table 2 lists the posterior sample mean of the estimated coefficients and P values for each M, and figure 8 shows the resulting $\bar{\theta}_{M}(t)$ approximations and time series estimates of P compared with the true $\theta(t)$ and P, respectively, when M = 1, 3, and 5.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 2. EnKF posterior sample means of the TVP model coefficients and period P for different values of M when approximating the sinusoidal forcing parameter $\theta(t) = 2\sin(t) - 0.5\cos(2t/3)$ in (18), using the data in figure 3 and treating the period of $\theta(t)$ as an additional unknown. Numbers are rounded to four decimal places. Note that the Fourier series coefficients obtained using the formulas in (5) and (6) to approximate the known function for orders 3 and above correspond to $c_4 = -0.5$, $c_5 = 2$, and $c_k = 0$ for all other coefficients.

Coefficient	M = 1	M = 2	M = 3	M = 4	M = 5
c0	−0.3468	0.1705	−0.0053	0.0962	0.0014
c1	0.5621	−0.0454	−0.0043	−0.0709	−0.1448
c2	0.2776	0.0085	−0.0760	0.0274	−0.1032
c3	—	1.7142	−0.0479	−0.0343	0.0747
c4	—	−0.0627	−0.5039	−0.5476	−1.3471
c5	—	—	2.0026	1.9663	2.4455
c6	—	—	0.0926	−0.4897	−0.8694
c7	—	—	—	0.0487	0.3266
c8	—	—	—	−0.2534	0.4760
c9	—	—	—	—	−0.0178
c10	—	—	—	—	0.3467
P	10.0386	12.5959	18.8658	18.8024	18.9363
Scaled RMSE	1.0373	0.3170	0.0554	0.2288	0.8050

Similar to the results obtained using M = 3 with a known period (shown in figure 4), we see in figure 7 that most of the coefficient estimates converge after assimilating about one period of data; the uncertainty encoded in the ±2 standard deviations around the mean is a bit wider but continues to decrease as more data are sequentially incorporated. The estimate of P takes more time to converge but, after assimilating about two periods of data, converges closely to the true underlying period (with a relative error of approximately $8.6124\times10^{-4}$). While not shown, the resulting time series estimates for both position and velocity are similar to those obtained in the known period case.

The results in table 2 and plots in figure 8 emphasize that, again for this case, M = 3 gives the best overall approximation to $\theta(t)$ (with smallest scaled RMSE) as well as the best estimate of P, with reasonably small uncertainty in this estimate by the end of the filtering process. When M = 1 and M = 2, the period estimates somewhat diverge from the truth and result in under-approximations. The results when M = 4, while not shown graphically, are similar to M = 3 but with more uncertainty in the posterior estimate of P and an increase in the scaled RMSE of the TVP approximation. When M = 5, the increasing uncertainty in the posterior estimate of P becomes more clear (as seen in figure 8), along with increased error in the TVP approximation. Following from these results, the corresponding mass-spring model predictions using different initial conditions with these TVP approximation models for $\theta(t)$ are most accurate when M = 3. However, if directly using the models when M = 4 or M = 5 in this case, the increased uncertainty in the period and corresponding increased error in the TVP approximations lead to increased error in the model predictions. This can be addressed by using the posterior mean estimate of P, fixing the period to this known value, and re-running the filtering process (now assuming a known period) to obtain improved TVP approximation model coefficient estimates.

In the previous simulations, we demonstrate the effectiveness of the proposed estimation method when approximating a periodic TVP, addressing situations when we know (and fix) and when we do not know (and estimate) the underlying period. Here we extend this approach to approximate more general TVP, in particular, considering cases when $\theta(t)$ is not periodic over the time interval of available data. As discussed in section 3.1, this prevents use of a direct formula for setting the ω_i values in our sinusoidal approximation model terms. Instead, here we choose a fixed increment ω and let $\omega_i = \omega i$, $i = 1,\dots, M$, in order to systematically include sine and cosine pairs with different frequencies in the approximation.

Using the same procedure as before, we simulate data from (18) over the time interval [0,60] with three different non-periodic forcing parameters:

• (i)  
  a linear polynomial, where

  $$\begin{align} \theta\left(t\right) = -0.07t + 2; \end{align} \tag{ 21 }$$
• (ii)  
  a cubic polynomial, where

  $$\begin{align} \theta\left(t\right) = 0.0001\left(t - 25\right)^3 - 0.001t^2 + 3; \end{align} \tag{ 22 }$$

  and
• (iii)  
  a step function, where

  $$\begin{align} \theta\left(t\right) = \begin{cases} -2 & t \unicode{x2A7D} 30 \\ \ 2 & t>30 \end{cases} . \end{align} \tag{ 23 }$$

Figure 9 shows the simulated data in each case.

Zoom In Zoom Out Reset image size

In building the TVP approximation models as in (8), we set an increment of ω = 0.01 and define $\omega_i = 0.01 i$ for each $i = 1,\dots, M$. Since each of the underlying $\theta(t)$ have different levels of complexity, we test a variety of M values in each case. Figure 10 shows the best resulting approximations (i.e. those with the smallest corresponding scaled RMSEs) for each non-periodic forcing parameter. For the linear forcing parameter, the best fit occurs when M = 1 (with scaled RMSE ≈ 0.0239), needing only three terms in the TVP approximation model to obtain an accurate approximation. Here, increasing M to somewhat larger values (e.g. M = 4 or 5) increases the approximation error and number of unknowns but still permits reasonable approximations with convergent coefficients.

Zoom In Zoom Out Reset image size

More terms are needed when approximating the cubic forcing parameter, and it is important to note that, for this example, using too small an M value results in coefficient estimates that do not converge over the time interval of available data. More specifically, the coefficients do not converge when M = 1, 2 or 3, but convergence improves beginning with M = 4. The best fit for the cubic forcing parameter occurs when M = 6 (with scaled RMSE ≈ 0.1596), but with a similarly good fit when M = 5 (scaled RMSE ≈ 0.1614) requiring the estimation of two less coefficients.

Estimating the step forcing parameter presents the most difficult challenge of the three, given the jump discontinuity in the function halfway through the time interval of available data. For this example, the TVP approximation model coefficients do not converge for smaller M values, but a significant increase in M leads to reasonable fits with convergent coefficients. The TVP approximations yield similar scaled RMSE values for M between 15 and 22, with the lowest occurring when M = 21 (scaled RMSE ≈ 0.2737). While not capturing the exact shape, the TVP approximations in this case are able to capture the jump point between constant parameter values in the underlying step function. The mass-spring model predictions for different initial conditions follow as expected using the best TVP approximation results over the time interval [0,60], but it becomes more difficult to accurately predict the behavior of the system after this time frame for non-periodic TVPs, since the parameters themselves will continue to dynamically change over time without updating their approximation models.

In this section, we apply the derivative-based augmentation approach described in section 3.2 to estimate the sinusoidal forcing parameter $\theta(t) = 2\sin(t) - 0.5\cos(2t/3)$ in (18) using the same data as in figure 3 and assuming a known period for $\theta(t)$. Compared to the approach used in section 4.1.1, the difference with this implementation is that we now treat $\theta(t)$ as unobserved system state, where

$$\begin{align} \frac{\textrm{d}\theta}{\textrm{d}t} = \theta^{^{\prime}}_M\left(t\right), \quad \theta\left(0\right) = \theta_0 \end{align} \tag{ A.1 }$$

with $\theta^{^{\prime}}_M(t)$ defined as in (14), and we track it along with p(t) and v(t), estimating the 2M coefficients $c_1,\dots,c_{2M}$ plus the TVP initial value θ₀ for a total of $2M+1$ unknown constant parameters. Figure A.1 shows the resulting EnKF time series estimates of position, velocity, and $\theta(t)$, along with the six unknown coefficients in (14) and the initial condition θ₀, when M = 3. Table A.1 lists the posterior sample mean of each coefficient and θ₀ for each M value, and figure A.2 shows the corresponding EnKF time series estimates of $\theta(t)$ in each case compared with the true $\theta(t)$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The results in figure A.1 highlight the differences between this derivative-based augmented systems implementation and the direct approximation model (with corresponding results shown in figure 4), where here $\theta(t)$ is treated as an unobserved system state and tracked along with position and velocity. The six coefficients converge after assimilating about one period of data, which is also reflected in the tracking of $\theta(t)$, while the estimate of θ₀ continues to improve over the remaining time interval of observations. As the plots in figure A.2 illustrate, the filter is not able to well track $\theta(t)$ when M = 1 and 2, noticeably under-estimating the dynamics with a damping effect beginning close to the one period mark (i.e. around time t = 19). The estimation significantly improves when M = 3, where the EnKF mean estimate very well captures the behavior of the true underlying parameter after assimilating one period of data, with fairly tight uncertainty bounds around the mean estimate. Similar behavior occurs when M = 4 and 5, although there is more initial error in the approximation over the first period and slightly wider uncertainty bounds after the coefficients converge. The results in table A.1 paint a similar picture for the posterior mean estimates of θ₀, where the estimate when M = 3 yields the smallest relative error (≈ 0.1440) when compared to the true initial value $\theta_0 = -0.5$.

In this section, we demonstrate the performance when applying the standard augmented EnKF algorithm described in section 2 to estimate the sinusoidal forcing parameter $\theta(t) = 2\sin(t) - 0.5\cos(2t/3)$ in (18) using the same data as in figure 3. Since the standard augmented EnKF is meant for estimating constant parameters, we anticipate that applying this algorithm to estimate the time-varying forcing parameter will not capture its dynamic behavior. As expected, figure A.3 shows that the standard augmented EnKF indeed fails to capture the dynamics of $\theta(t)$ and instead converges to a constant that represents somewhat of an average value. Note that the filter still reasonably well tracks both states p(t) and v(t); however, some of the behavior towards the end of the time series is not captured as the estimate is somewhat damped down (corresponding to the underestimated forcing).

Zoom In Zoom Out Reset image size

To account for time-varying behavior in the parameter estimate, we can modify the augmented EnKF algorithm to include parameter tracking (see, e.g. [10]) by incorporating a random walk model for the parameter in the filter prediction step, such that

$$\begin{align} \theta_{j+1}^{\left(n\right)} = \theta_j^{\left(n\right)} + \xi_{j+1}^{\left(n\right)}, \quad \xi_{j+1}^{\left(n\right)}\sim\mathcal{N}\left(0,\sigma^2\right) \end{align} \tag{ A.2 }$$

for each ensemble member, $n = 1,\dots, N$, where $\xi_{j+1}^{(n)}$ defines the parameter perturbation from time j to j + 1. The success of parameter tracking with the EnKF relies significantly on the choice of the parameter noise model variance σ² [10], which must be manually selected before running the filter. Figure A.4 shows the resulting estimates of $\theta(t)$ when using parameter tracking with four different choices of σ (namely, when σ = 0.1, 0.5, 1, and 5) and keeping all other filter settings (ensemble size, etc) the same as in the previous numerical experiments. When σ = 0.1, we see that the underlying form and magnitude of $\theta(t)$ are not well captured. When σ = 0.5, the approximation improves but the form is still not well captured and the dynamics of the estimate lag behind the truth. For the larger values when σ = 1 and σ = 5, note that the overall dynamics of $\theta(t)$ are generally captured (with more error as σ increases) but the functional form is not well approximated in any case and the uncertainty bounds grow increasingly with σ. While some of these parameter estimates may be reasonable overall, note that this algorithm fails to well approximate the underlying functional form of $\theta(t)$. Further, since σ is fixed, we do not see a decrease in the EnKF uncertainty bounds over time. Comparing these results with those using the Fourier series-based approaches presented in this work (perhaps most visually comparable with the derivative-based augmentation results shown in appendix A.1), we see that the Fourier series-based approximations have improved performance in capturing the functional form of the underlying TVP and by design give an approximation formula for $\theta(t)$ (or $\textrm{d}\theta/\textrm{d}t$) that can be used in making forward predictions with the system. We also see a decrease in uncertainty in the TVP approximation model coefficient estimates over time, as well as in the $\theta(t)$ approximations shown in figure A.2 with the best fits resulting when M = 3, 4, and 5.

Zoom In Zoom Out Reset image size

As an additional test system, we consider a modified version of the Lorenz-63 equations [51], commonly studied in data assimilation as a simplified model for atmospheric convection that exhibits chaotic dynamics under certain parameterizations. In particular, we consider a periodically-forced Lorenz system of the type in [52] with governing equations given by

$$\begin{eqnarray} \frac{\textrm{d}x_1}{\textrm{d}t} & = & \sigma\left(x_2 - x_1\right) \end{eqnarray} \tag{ A.3 }$$

$$\begin{eqnarray} \frac{\textrm{d}x_2}{\textrm{d}t} & = & \rho\theta\left(t\right) x_1 - x_1 x_3 - x_2 \end{eqnarray} \tag{ A.4 }$$

$$\begin{eqnarray} \frac{\textrm{d}x_3}{\textrm{d}t} & = & x_1 x_2 - \beta x_3 \end{eqnarray} \tag{ A.5 }$$

where the constant parameters σ = 10, ρ = 28, and $\beta = 8/3$ are chosen so that the system exhibits chaotic behavior and $\theta(t)$ is a periodic forcing function. Assuming that σ, ρ, and β are known, our goal is to estimate $\theta(t)$ using the Fourier series-based TVP approximation technique. While in the previous examples the TVP of interest was additive with respect to the system states, note that here $\theta(t)$ is a multiplicative TVP. Further, here we assume that the system is partially observed, such that we are able to obtain observations of $x_2(t)$ and $x_3(t)$ (relating to the temperature variations) at some discrete times, but we cannot observe $x_1(t)$ (relating to convection) and thereby must track the behavior of $x_1(t)$ along with estimating the unknown TVP. Letting $\theta(t) = 1 + \varepsilon\cos(\omega t)$ (as studied in [52]) with ɛ = 0.5 and ω = 0.5 (corresponding to period $P = 4\pi \approx 12.57$), we simulate data from the forced Lorenz system in (A.3)–(A.5) by taking observations of the states x₂ and x₃ every 0.02 time units over the interval [0,50] (covering about 4 periods of $\theta(t)$) and corrupting with Gaussian noise (zero mean and variance taken to be 20% of the standard deviation of the true states, as before). Figure A.5 shows the Lorenz butterfly of the reference solutions (uncorrupted) and the forcing parameter $\theta(t)$.

Zoom In Zoom Out Reset image size

Assuming a known period for $\theta(t)$, we apply the method in section 3.1 using the formula in (10) to set $\omega_i = 0.5{i}$, $i = 1,\dots, M$, and employ the TVP approximation model in (8) with increasing values of M, ranging from M = 1 (which involves three coefficients) through M = 5 (11 coefficients). Note that the Fourier series coefficients obtained when using the formulas in (5) and (6) to approximate the function $h(t) = 1 + 0.5\cos(0.5 t)$, assuming that both h(t) and $P = 4\pi$ are known, gives a close approximation when Q = 1 with coefficients $c_0 = 1$, $c_1 = 0$, and $c_2 = 0.5$ (and using higher Q values yields similar approximations with $c_k = 0$ for the additional coefficients), so we expect comparable results using the Fourier series-based TVP approximation with a known period to estimate $\theta(t)$. We apply the EnKF similarly as before with an ensemble size of N = 100, with the aim of tracking the model states and estimating the unknown coefficients comprising the Fourier series-based model approximation of $\theta(t)$ for each value of M. Table A.2 lists the posterior sample mean of each coefficient for each M value, and figure A.6 shows the resulting $\bar{\theta}_{M}(t)$ approximations in each case compared with the true $\theta(t)$.

Zoom In Zoom Out Reset image size

As demonstrated in table A.2, the filter well estimates the TVP approximation model coefficients, resulting in similar values to those expected from Fourier series expansion (with small approximation errors) for each M value. While not shown, the filter also well tracks all three model states (including the unobserved state x₁) for each M value. As shown in figure A.6, while the lowest RMSE occurs when M = 4 for this example, all of these M values result in reasonably close model approximations to the true $\theta(t)$ in this case. While we used an ensemble of size N = 100 to compute the results as shown, we remark that similar results can be obtained using smaller ensembles (e.g. N = 20 and above when M = 1) but the performance degrades for too small an ensemble relative to the number of unknowns. Though these initial results are promising, in future work we will perform additional studies to verify the feasibility of the proposed technique for TVP estimation on larger-scale forced chaotic systems with more unknowns.