Data assimilation for massive autonomous systems based on a second-order adjoint method

Data assimilation (DA) is a fundamental computational technique that integrates numerical simulation models and observation data on the basis of Bayesian statistics. Originally developed for meteorology, especially weather forecasting, DA is now an accepted technique in various scientific fields. One key issue that remains controversial is the implementation of DA in massive simulation models under the constraints of limited computation time and resources. In this paper, we propose an adjoint-based DA method for massive autonomous models that produces optimum estimates and their uncertainties within reasonable computation time and resource constraints. The uncertainties are given as several diagonal elements of an inverse Hessian matrix, which is the covariance matrix of a normal distribution that approximates the target posterior probability density function in the neighborhood of the optimum. Conventional algorithms for deriving the inverse Hessian matrix require $O(C{N}^2+N^3)$ computations and $O\left(N^2\right)$ memory, where $N$ is the number of degrees of freedom of a given autonomous system and $C$ is the number of computations needed to simulate time series of suitable length. The proposed method using a second-order adjoint method allows us to directly evaluate the diagonal elements of the inverse Hessian matrix without computing all of its elements. This drastically reduces the number of computations to $O\left(C\right)$ and the amount of memory to $O\left(N\right)$ for each diagonal element. The proposed method is validated through numerical tests using a massive two-dimensional Kobayashi phase-field model. We confirm that the proposed method correctly reproduces the parameter and initial state assumed in advance, and successfully evaluates the uncertainty of the parameter. Such information regarding uncertainty is valuable, as it can be used to optimize the design of experiments.

Figure 1

Time evolution of phase field $\varphi$ starting from the initial state shown in (a) for $m=0.1$. (b)–(f) $\varphi$ at $t=5.0\tau ,10.0\tau ,20.0\tau ,50.0\tau$, and $100.0\tau$, respectively. The color indicates the magnitude of $\varphi$.

Figure 2

Results of twin experiment I-(i). The length of each error bar for the optimum parameter $\widehat{m}$ in (a) corresponds to the estimated uncertainty $\delta \widehat{m}$ in (b). The uncertainty cannot be determined when $T_{\text{max}}=0.2\tau$. The dashed black line in (b) indicates a power function of order $-1.5$.

Figure 3

Results of twin experiment I-(ii). The length of each error bar for optimum parameter $\widehat{m}$ in (a) corresponds to the estimated uncertainty $\delta \widehat{m}$ in (b). The dashed black line in (b) indicates a power function of order of 0.5.

Figure 4

Results of twin experiment I-(iii). The length of each error bar for optimum parameter $\widehat{m}$ in (a) corresponds to the estimated uncertainty $\delta \widehat{m}$ in (b). The dashed black line in (b) indicates a linear function.

Figure 5

Results of twin experiment II. (a) How the LBFGS method updates $m$ when the observation data noise has a small (solid red line) and large (dashed green line) standard deviation. (b) The optimum initial state of the phase field ${\varphi }_i\left(0\right)$ in the case of small noise, and (c) that in the case of large noise.

Figure 6

Results of twin experiment III. Estimated initial states of the phase field ${\varphi }_i\left(0\right)$ (a) after the 31st step and (b) after the final step in the iteration of the LBFGS method. (c) Improvement in the cost function $J^{\prime }$.

Figure 7

Phase diagram for an axisymmetric two-dimensional Kobayashi PF model. The solid red line indicates the critical radius as a function of the parameter $m$. The region above or below the critical line corresponds to a spot growing or decaying with time, respectively.

Figure 8

Time evolutions of the phase fields starting from different initial states ${\varphi }_i\left(0\right)$: (a) the estimated initial state obtained after the 31st step [Fig. 6], (b) that after the final step [Fig. 6], and (c) the true initial state [Fig. 1].