weak_4D.rst

Weak-constraint 4D-Var

Nonlinear Cost Function

This note describes the formulation of weak constraint 4D-Var, and its implementation in OOPS.

Let J denote the cost function and x_k denote the state at time t_k. The weak constraint 4D-Var cost function (equation cost-nonlin) is a function a set of model states x_k defined at regular intervals over a time window [t₀, t_N). We refer to this set of states as the "four-dimensional state vector". Note that there are N sub-windows and that the four-dimensional state vector contains the N states x₀, …, x_N − 1.

$$\begin{aligned} \begin{eqnarray} && J( \mathbf{x}_0 , \mathbf{x}_1 , \ldots , \mathbf{x}_{N-1} ) = \\\ &&\qquad \phantom{+} \frac{1}{2} \left( \mathbf{x}_0 - \mathbf{x}_b \right)^{\rm T} \mathbf{B}^{-1} \left( \mathbf{x}_0 - \mathbf{x}_b \right) \nonumber \\\ &&\qquad + \frac{1}{2} \sum_{k=0}^{N-1} \sum_{j=0}^{M-1} \left( \mathbf{y}_{k,j} - {\cal H}_{k,j} (\mathbf{x}_{k,j} ) \right)^{\rm T} \mathbf{R}_{k,j}^{-1} \left( \mathbf{y}_{k,j} - {\cal H}_{k,j} (\mathbf{x}_{k,j} ) \right) \nonumber \\\ &&\qquad + \frac{1}{2} \sum_{k=1}^{N-1} \left( \mathbf{q}_k - \mathbf{\bar q} \right)^{\rm T} \mathbf{Q}_k^{-1} \left( \mathbf{q}_k - \mathbf{\bar q} \right) \nonumber \end{eqnarray} \end{aligned}$$

The cost function has three terms. The first term penalizes departures of x₀ from a prior estimate (the "model background"), x_b. The matrix B is the covariance matrix of background error.

The second term of the cost function penalizes the discrepancies between observations y_k, j and their model equivalents, x_k, j. Here, the double subscript denotes a time such that (for j = 1…M − 1) t_k, j ∈ [t_k, t_k + 1). We assume that t_k, 0 = t_k and t_k, M = t_k + 1. We refer to the interval [t_k, t_k + 1) as a "sub-window", and the times t_k, j as "observation time slots". The operator ${\cal H}_{k,j}$ is the (nonlinear) observation operator, and R_k, j is the covariance matrix of observation error.

The final term in the cost function penalizes departures of model error q_k, defined at the boundaries between sub-windows, from a separately-estimated systematic model error, q̄. The matrix Q_k is the covariance matrix of model error.

The states at times t_k, j are given by an un-forced integration of the model:

$$\begin{equation} \mathbf{x}_{k,j} = {\cal M}_{k,j} (\mathbf{x}_{k,j-1} ) \end{equation}$$

for k = 0, …, N − 1 and j = 1, …, M.

The model errors are determined as the difference between the state at the start of a sub-window and the corresponding state at the end of the preceding sub-window. There are N − 1 model errors corresponding to the N − 1 boundaries between sub-windows:

q_k = x_k, 0 − x_{k − 1, M}  for k = 1, …, N − 1.

Linear (Incremental) Cost Function

The linear (incremental) cost function is defined by linearizing the operators in equation cost-nonlin around a "trajectory", to give a quadratic approximation Ĵ of the cost function. In principle, the trajectory is known at every timestep of the model. In practice, it may be necessary to assume the trajectory remains constant over small time intervals.

Let us denote by x_k^t the trajectory state at time t_k. The approximate cost function is expressed as a function of increments δx_k to these trajectory states:

$$\begin{aligned} \begin{eqnarray} && {\hat J} ( \delta\mathbf{x}_0 , \delta\mathbf{x}_1 , \ldots , \delta\mathbf{x}_{N-1} ) = \\\ &&\qquad \phantom{+} \frac{1}{2} \left( \delta\mathbf{x}_0 + \mathbf{x}^t_0 - \mathbf{x}_b \right)^{\rm T} \mathbf{B}^{-1} \left( \delta\mathbf{x}_0 + \mathbf{x}^t_0 - \mathbf{x}_b \right) \nonumber \\\ &&\qquad + \frac{1}{2} \sum_{k=0}^{N-1} \sum_{j=0}^{M-1} \left( \mathbf{d}_{k,j} - \mathbf{H}_{k,j} (\delta\mathbf{x}_{k,j} ) \right)^{\rm T} \mathbf{R}_{k,j}^{-1} \left( \mathbf{d}_{k,j} - \mathbf{H}_{k,j} (\delta\mathbf{x}_{k,j} ) \right) \nonumber \\\ &&\qquad + \frac{1}{2} \sum_{k=1}^{N-1} \left( \delta\mathbf{q}_k + \mathbf{q}^t_k - \mathbf{\bar q} \right)^{\rm T} \mathbf{Q}_k^{-1} \left( \delta\mathbf{q}_k + \mathbf{q}^t_k - \mathbf{\bar q} \right) \nonumber \end{eqnarray} \end{aligned}$$

Here, H_k, j is a linearization of ${\cal H}_{k,j}$ about the trajectory x_k, j^t. The vector d_k, j is defined as:

$$\begin{equation} \mathbf{d}_{k,j} = \mathbf{y}_{k,j} - {\cal H}_{k,j} (\mathbf{x}^t_{k,j} ) . \end{equation}$$

The increments δx_k, j satisfy a linearized version of equation nonlin-propagate-eqn:

δx_k, j = M_k, jδx_{k, j − 1}

for k = 0, …, N − 1 and j = 1, …, M.

The model error increments are given by:

δq_k = δx_k, 0 − δx_{k − 1, M}  for k = 1, …, N − 1.

Solution Algorithm: Outer Loop

The cost function (equation cost-nonlin) is minimized by successive quadratic approximations according to the following algorithm:

Given an initial four-dimensional state {x₀, x₁, …, x_N − 1}:

1. For each sub-window, integrate equation nonlin-propagate-eqn from the initial condition x_k, 0 = x_k, to determine the trajectory and the state x_k, M at the end of the sub-window.
2. Calculate the model errors from equation nonlin-model-error-eqn.
3. Minimize the linear cost function (equation cost-linear) to determine the increments δx_k and δq_k.
4. Set x_k := x_k + δx_k and q_k := q_k + δq_k.
5. Repeat from step 1.

Solution Algorithm: Inner Loop

There are several possibilities for minimizing the linear cost function. Some of these are described in the following sub-sections.

Initial State and Forcing Formulation

The initial state and forcing formulation expresses the linear cost function as a function of the initial increment δx₀ and the model error increments δq₁, …, δq_N − 1.

The control vector for the minimization comprizes a set of three-dimensional vectors χ_k for k = 0, …, N − 1, and is defined by:

$$\begin{eqnarray} \mathbf{B}^{1/2} \mathbf{\chi}_0 &=& \left( \delta\mathbf{x}_0 + \mathbf{x}^t_0 - \mathbf{x}_b \right) \end{eqnarray}$$

$$\begin{eqnarray} \mathbf{Q}_k^{1/2} \mathbf{\chi}_k &=& \left( \delta\mathbf{q}_k + \mathbf{q}^t_k - \mathbf{\bar q} \right) \qquad \mbox{for $k=1, \ldots , N-1$} \end{eqnarray}$$

The background and observation terms of the cost function can be evaluated directly from the control vector as:

$$\begin{equation} J_b = \frac{1}{2} \mathbf{\chi}_0^{\rm T} \mathbf{\chi}_0 \qquad\mbox{and}\qquad J_q = \frac{1}{2} \sum_{k=1}^{N-1} \mathbf{\chi}_k^{\rm T} \mathbf{\chi}_k \end{equation}$$

The contribution to the gradient of the cost function from these terms is simply equal to the control vector itself.

To evaluate the observation term, we must generate the four-dimensional increment: {δx_k; k = 0, …, N − 1}. This is done by first calculating δx₀ and δq_k from equations psi-to-dx0-eqn and psi-to-dq-eqn, and then generating δx_k using equations linear-propagate-eqn and linear-model-error-eqn.

The cost function is minimized, resulting in an updated control vector, corresponding to the minimum of the linear cost function. From this cost function, we must generate the four-dimensional increment required by the outer loop. This is done using equations linear-propagate-eqn and linear-model-error-eqn, and requires a integration of the tangent linear model

Four Dimensional State Formulation

In the four dimensional state formulation, the cost function is expressed as a function of the four dimensional state, δx₀, …, δx_N − 1.

The control vector for the minimization is defined by:

$$\begin{aligned} \begin{eqnarray} \mathbf{B}^{1/2} \mathbf{\chi}_0 &=& \left( \delta\mathbf{x}_0 + \mathbf{x}^t_0 - \mathbf{x}_b \right) \\\ \mathbf{Q}_k^{1/2} \mathbf{\chi}_k &=& \delta\mathbf{x}_k \qquad \mbox{for $k=1, \ldots , N-1$} \end{eqnarray} \end{aligned}$$

With this choice for χ₀, the background term of the cost function can be evaluated as $J_b = \frac{1}{2} \mathbf{\chi}_0^{\rm T} \mathbf{\chi}_0$. However, the model error term must now be evaluated explicity as:

$$\begin{equation} J_q = \frac{1}{2} \sum_{k=1}^{N-1} \left( \delta\mathbf{q}_k + \mathbf{q}^t_k - \mathbf{\bar q} \right)^{\rm T} \mathbf{Q}_k^{-1} \left( \delta\mathbf{q}_k + \mathbf{q}^t_k - \mathbf{\bar q} \right) \end{equation}$$

where δq_k is determined from equation linear-model-error-eqn.

Note that this requires that the inverse model error covariance matrix is available, and is reasonably well conditioned.