# Introducing Periodic Parameters in a Marine Ecosystem Model Using Discrete Linear Quadratic Control

Mustapha El Jarbi, Thomas Slawig, Andreas Oschlies

► To cite this version:

Mustapha El Jarbi, Thomas Slawig, Andreas Oschlies. Introducing Periodic Parameters in a Marine Ecosystem Model Using Discrete Linear Quadratic Control. 25th System Modeling and Optimization (CSMO), Sep 2011, Berlin, Germany. pp.481-490, 10.1007/978-3-642-36062-6_48. hal-01347570

HAL Id: hal-01347570

https://inria.hal.science/hal-01347570

Submitted on 21 Jul 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Introducing Periodic Parameters in a Marine Ecosystem Model using Discrete Linear Quadratic Control

Mustapha El Jarbi, Thomas Slawig, and Andreas Oschlies

Institute for Computer Science, Christian-Albrechts Universitaet zu Kiel,
24098 Kiel, Germany
GEOMAR, Düsternbrooker Weg 20, 24105 Kiel
{mej,ts}@informatik.uni-kiel.de
{aoschlies}@ifm-geomar.de

**Abstract.** This paper presents the application of the *Discrete Linear Quadratic Control (DLQC)* method for a parameter optimization problem in a marine ecosystem model. The ecosystem model simulates the distribution of nitrogen, phytoplankton, zooplankton and detritus in a water column with temperature and turbulent diffusivity profiles taken from a three-dimensional ocean circulation model. We present the linearization method which is based on the available observations. The linearization is necessary to apply the DLQC method on the nonlinear system of state equations. We show the form of the linearized time-variant problems and the resulting two algebraic Riccati Equations. By using the DLQC method, we are able to introduce temporally varying periodic model parameters and to significantly improve – compared to the use of constant parameters – the fit of the model output to given observational data.

**Keywords:** Optimal Control, Non-linear Systems, Parameter Optimization, Biogeochemical Modelling, Discrete Linear Quadratic Regulator Problem, Periodic Parameter, Discrete Riccati Equation

## 1 Introduction

We consider nonlinear partial differential diffusion-advection systems of the form

$$ \frac{\partial x^i}{\partial t} = -w^i \frac{\partial x^i}{\partial z} + \frac{\partial}{\partial z} \left( \nu_{\rho} \frac{\partial x^i}{\partial z} \right) + q^i(\mathbf{x}, \mathbf{u}), \quad i = 1, 2, 3, 4 \quad (1) $$

$$ x^i: [0, T] \times [-H, 0] \rightarrow \mathbb{R}. $$

Here $z$ denotes the vertical spatial coordinate, $H$ the depth in the water column, $q^i$ represents the biogeochemical coupling terms for the four species and $\mathbf{u} = (u_1, \dots, u_p)$ is the vector of unknown physical and biological parameters. The circulation data are the turbulent mixing coefficient $\nu_{\rho} = \nu_{\rho}(z, t)$ and the temperature $\Theta = \Theta(z, t)$, which goes into the non-linear coupling terms $q^i$, see
(3). The vertical sinking velocity $w^i$ is a parameter of the biological model that is nonzero only for $x^4$, i.e. $w^1 = w^2 = w^3 = 0$, $w^4 = ws > 0$.

The state of the system is denoted by $\mathbf{x} = (x^1, x^2, x^3, x^4)^\top$ and the control by $\mathbf{u}$. A control problem is defined as

$$ \min_{\mathbf{u}} \mathcal{F}(\mathbf{x}, \mathbf{u}) \quad \text{subject to} \quad (1), \qquad (2) $$

where $\mathcal{F}$ is a cost functional which will be introduced later.

Our main goals are:

- to minimize a least-squares type cost functional,

- to allow the parameters to vary temporally over the year while remaining periodic over all years of the considered time interval.

The work presented in this paper is motivated by results obtained for a typical marine ecosystem model, namely the NPZD model introduced in [1], [2]. As was reported in several publications with different optimization algorithms, the quality of the model-to-dat fit was not optimal, and in some cases it was difficult to identify the parameters uniquely, see for example [3],[5],[4]. In most cases, the parameters of the marine ecosystem models are assumed to be temporally constant. This reflects the aim to obtain a model that is applicable for arbitrary time intervals. To solve this problems, we discretize and linearize the nonlinear state (1) around a reference trajectories and we interpret it as a Discrete Linear Quadratic Control (DLQC) problem. Therein, we allow the parameters to be time-dependent, apply a well-established method for optimal control, and additionally impose the constraint of annual periodicity. This avoids the process of parametrization in the sense that we do not have to know or assume how the above mentioned periodic functions look like. In contrast, the method itself will generate an optimal periodic function for each parameter. Moreover, it allows to balance the two aims that we have: By introducing weight matrices we can choose if it is more important to obtain a very good fit or nearly perfect periodicity. The method requires a reference trajectory and a reference control, i.e., the vector of model parameters. The former can be taken from the measurement data, and for the latter we use an initial guess for the parameters which can be the output of an optimization with constant parameters. The outline of this paper is as follows. In the next section we briefly described the model Equation and optimization problem (2), the DLQC problem formulation in section 3. A application of the DLQC method on the NPZD model is presented in section 4.3. Afterwards, we present our results with respect to the quality of the fit and the periodicity of the parameters and end the paper with some conclusions.

## 2 Model Equations and Optimization Problem

In this section we give the formulations of the NPZD model and of the corresponding parameter optimization problem and we formulate the optimization problem for the discrete model.
## 2.1 Model Equations

This section describes the ecosystem model. The considered system (1) is a spatially one-dimensional marine biogeochemical model, that simulates the interaction of dissolved inorganic nitrogen N, phytoplankton P, zooplankton Z and detritus D. It was developed with the aim of simultaneously reproducing observations at three North Atlantic locations by the optimization of free parameters within credible limits, see [4]. The model uses the ocean circulation and temperature field in an off-line modus, i.e. these are used only as forcing, but no feedback on them is modeled. The model simulates one water column at a given horizontal position, which is motivated by the fact that there have been special time series studies at fixed locations, one of which was used here. In the model, the concentrations (in mmol N m⁻³) of dissolved inorganic nitrogen N, phytoplankton P, zooplankton Z, and detritus D, denoted by **x** = (xᵢ)ᵢ=₁,...,₄ = (N, P, Z, D) are described by the PDE system (1).

The biogeochemical source-minus-sink terms **q** = (qᵢ)ᵢ=₁,...,₄ are explicit by given in [1]:

$$
\left.
\begin{aligned}
q^1(\mathbf{x}, \mathbf{u}) &= -\bar{J}(z, t, N)P + \gamma_2 Z + \mu_D D, \\
q^2(\mathbf{x}, \mathbf{u}) &= \bar{J}(z, t, N)P - \mu_X P - G(P)Z, \\
q^3(\mathbf{x}, \mathbf{u}) &= \gamma_1 G(P)Z - \gamma_2 Z - \mu_Z Z^2, \\
q^4(\mathbf{x}, \mathbf{u}) &= (1 - \gamma_1)G(P)Z - \mu_Z Z^2 + \mu_X P - \mu_D D - w_s \frac{\partial D}{\partial z}
\end{aligned}
\right\}
\quad (3)
$$

where $\bar{J}$ is the daily averaged phytoplankton growth rate as a function of depth $z$ and time $t$, and $G$ is the grazing function (see below). The remaining parameter in the above equations are defined in [1],

$$ G(\epsilon, g) = \frac{g\epsilon P^2}{g + \epsilon P^2} \qquad \bar{J}(z, t, N) = \min \left( L(z, t), J_{\max} \frac{N}{K_1 + N} \right), \quad (4) $$

where $L$ denotes the purely light-limited growth rate, and $J_{\max}$ is the light-saturated growth. For more details of $\bar{J}$, $L$ and the parameters see [1], [4].

## 2.2 The Optimization Problem

The aim of the optimization is to fit the aggregated model output $\mathbf{y} = C\mathbf{x}$ ($C$ is called the output matrix) to the given observational data $\mathbf{y}^{obs}$. There are five types of measurement data $\mathbf{y}^{obs} = (y_m^{obs})_{m=1,...,5}$, which correspond to aggregated values $\mathbf{y} := (y_m)_{m=1,...,5}$ of the model output see also [3]. Thus the cost function can be written as:

$$ \mathcal{F}(\mathbf{x}, \mathbf{u}) := \|C\mathbf{x} - \mathbf{y}^{\text{obs}}\|_{2,\sigma}, \quad (5) $$

where $\|\cdot\|_{2,\sigma}$ is a Euclidean norm weighted using the vector

$$ \sigma = (\sigma_l)_{l=1,...,5} = (0.1, 0.01, 0.01, 0.0357, 0.025) $$

of uncertainties corresponding to the five types of measurement data.
## 3 DLQC Problem Formulation

We use a discrete linear time-varying (LTV) system, i.e. we assume that the dynamical system is already discretized in time, namely at discrete times $t_k, k = 1, \dots, M$. In the context of the DLQC, one usually considers a discrete-time system of the form:

$$ \begin{aligned} \mathbf{x}_{k+1} &= A_k \mathbf{x}_k + B_k \mathbf{u}_k, \quad k = 1, 2, \dots, M-1 \\ x_1 &\text{(the given initial value),} \end{aligned} \quad (6) $$

where in every time step $k$

- $\mathbf{x}_k = \mathbf{x}(t_k) \in \mathbb{R}^n$ is called the state vector (here the model output),

- $\mathbf{u}_k = \mathbf{u}(t_k) \in \mathbb{R}^p$ is the control (here the model parameter) vector, with the parameter vector from the model (3).

- The matrix $A_k \in \mathbb{R}^{n \times n}$ and $B_k \in \mathbb{R}^{n \times p}$ are called the system matrix and the input matrix, respectively.

We will use the notations

$$ \begin{aligned} \mathbf{x} &= (\mathbf{x}_k)_{k=1,\dots,M} \in \mathbb{R}^{M \times n} \cong \mathbb{R}^{Mn}, \\ \mathbf{u} &= (\mathbf{u}_k)_{k=1,\dots,M-1} \in \mathbb{R}^{(M-1) \times p} \cong \mathbb{R}^{(M-1)p} \end{aligned} $$

for the whole discrete trajectories of state and control vector, respectively. The quadratic cost function of this optimal control problem is defined by:

$$ \mathcal{J}(\mathbf{u}) = \frac{1}{2} \mathbf{x}_M^\top Q_M \mathbf{x}_M + \frac{1}{2} \sum_{k=1}^{M-1} \mathbf{x}_k^\top Q_k \mathbf{x}_k + \mathbf{u}_k^\top R_k \mathbf{u}_k, \quad (7) $$

where in every time step $k$

- $Q_k$ is a positive semidefinite diagonal weighting matrix for the state vector for every model time step $k = 1, \dots, M$,

- $R_k$ is a positive definite diagonal weighting matrix for the control vector for every model time step $k = 1, \dots, M-1$.

For the solution of a discrete linear quadratic optimal control problem with LTV systems, there exists the following theorem, see [6].

**Theorem 1** If the $Q_k, k = 1, \dots, M$, are positive semi-definite and the $R_k, k = 1, \dots, M-1$, are positive definite, then there exists a unique solution of the DLQC (6), (7). The optimal control is given by the feedback law

$$ \begin{aligned} \mathbf{u}_k &= -K_k \mathbf{x}_k, \quad k = 1, \dots, M-1. \\ K_k &:= (R_k + B_k^\top \mathbf{X}_{k+1} B_k)^{-1} B_k^\top \mathbf{X}_{k+1} A_k, \quad k = 1, \dots, M-1 \\ \mathbf{x}_{k+1} &= (A_k - B_k K_k) \mathbf{x}_k, \quad k = 1, \dots, M-1. \end{aligned} $$

where the $(\mathbf{X}_k)_{k=1,\dots,M-1}$, is the unique symmetric solution of the Discrete Riccati Equation (DRE).

$$ \left. \begin{array}{l} \mathbf{X}_k = Q_k + A_k^\top \mathbf{X}_{k+1} A_k - A_k^\top \mathbf{X}_{k+1} B_k (R_k + B_k^\top \mathbf{X}_{k+1} B_k)^{-1} B_k^\top \mathbf{X}_{k+1} A_k, \\ \qquad k = 1, \dots, M-1. \end{array} \right\} (8) $$
# 4 Application of DLQC to the NPZD Model

In this section we apply the LQOC method to the discretized version of the NPZD model. We present the details of discretization, linearization and the enforcement of the periodicity of the parameters (controls).

## 4.1 Discretization Scheme

We use a discrete linear quadratic control (DLQC). For this purpose we present the original discretization scheme of the model.

The NPZD model is forced by output from the OCCAM global circulation model, namely the hourly vertical profiles of temperature t and vertical diffusivity $v_p$. The vertical grid consists of 32 layers with thickness increasing with depth. The time integration of the system (1) is performed by an operator splitting method:

- At first, the nonlinear coupling operators $q_k = (q_k^1, q_k^2, q_k^3, q_k^4)_{k=1,...,M-1}^\top$ are computed at every spatial grid point and integrated by four explicit Euler steps, each of which is described by the operator:

$$B_k(\mathbf{x}_k, \mathbf{u}_k) := (\mathbf{x}_k + \frac{\tau}{4}\mathbf{q}_k(\mathbf{x}_k, \mathbf{u}_k)). \quad (9)$$

This gives an intermediate iterate

$$\hat{\mathbf{x}}_k := B_k \circ B_k \circ B_k \circ B_k(\mathbf{x}_k, \mathbf{u}_k).$$

- Then, an explicit Euler step with full step-size $\tau$ is performed for the sinking term, which is spatially discretized by an upstream scheme. This step is summarized in a matrix $S$. Since the sinking velocity is temporally constant, this matrix does not depend on the time step $k$. Thus, at the end of this step, an intermediate tracer vector $\tilde{\mathbf{x}}_k$ is computed as

$$\tilde{\mathbf{x}}_k := S\hat{\mathbf{x}}_k, \quad (10)$$

where $S = (I_k + \tau A^{adv})$.

- Finally, an implicit Euler step is applied for the diffusion operator discretized with second order central differences. The resulting matrix $D_k$ for the diffusion depends on $k$ since the diffusion coefficient depends on time. The matrix is tridiagonal, and the system is solved directly for $\mathbf{x}_{k+1}$

$$\tilde{D}_k \mathbf{x}_{k+1} = \tilde{\mathbf{x}}_k, \quad (11)$$

where $\tilde{D}_k = (I_k - \tau D_k)\mathbf{x}_{k+1}$.

Summarizing, the discrete system can be written as

$$\begin{aligned} \mathbf{x}_{k+1} &= \tilde{D}_k^{-1} S B_k \circ B_k \circ B_k \circ B_k(\mathbf{x}_k, \mathbf{u}_k) \\ &= \tilde{D}_k^{-1} S G(\mathbf{x}_k, \mathbf{u}_k), \quad k = 1, \dots, M-1, \end{aligned} \quad (12)$$

The function $G$ is nonlinear and represents the discretized source minus sink terms.
## 4.2 Linearization of the Model

The LDQC approach is based on a linearization of (12) to obtain a linear time-varying problem. The linearization is performed around reference trajectories $(\mathbf{x}_k^r, \mathbf{u}_k^r)_{k=1,...,M-1}$. For the reference state trajectory we take available the observational data, is taken from the Bermura Atlantic Time-series Study (BATS) see also [7], the choice of the reference control trajectory is described in below. The linearized state equation now reads

$$ \tilde{\mathbf{x}}_{k+1} = A_k \tilde{\mathbf{x}}_k + B_k \mathbf{v}_k + r_k, \quad k = 1, \dots, M-1, \quad (13) $$

where

$$ \begin{aligned} A_k &= \tilde{D}_k^{-1} S \frac{\partial G}{\partial x}(\mathbf{x}_k^r, \mathbf{u}_k^r), & A_k &\in \mathbb{R}^{n \times n} \\ B_k &= \tilde{D}_k^{-1} S \frac{\partial G}{\partial u}(\mathbf{x}_k^r, \mathbf{u}_k^r), & B_k &\in \mathbb{R}^{n \times p}, \\ r_k &= \tilde{D}_k^{-1} S G(\mathbf{x}_k^r, \mathbf{u}_k^r) - \mathbf{x}_{k+1}^r, & r_k &\in \mathbb{R}^n \\ \tilde{\mathbf{x}}_k &= \mathbf{x}_k - \mathbf{x}_k^r, & \mathbf{v}_k &= \mathbf{u}_k - \mathbf{u}_k^r, & \tilde{\mathbf{x}}_k &\in \mathbb{R}^n, \mathbf{v}_k &\in \mathbb{R}^p. \end{aligned} $$

Now we write the linearized problem in the form of a (LDQC) problem, therefore we set:

$$ \hat{\mathbf{x}}_k := \begin{pmatrix} \tilde{\mathbf{x}}_k \\ 1 \end{pmatrix}, \quad \hat{A}_k := \begin{pmatrix} A_k & r_k \\ 0 & 1 \end{pmatrix}, \quad \hat{B}_k := \begin{pmatrix} B_k \\ 0 \end{pmatrix}, \quad \hat{Q}_k := \begin{pmatrix} Q_k & 0 \\ 0 & 0 \end{pmatrix} $$

where $\hat{\mathbf{x}}_k \in \mathbb{R}^{n+1}$, $\hat{A}_k \in \mathbb{R}^{(n+1) \times (n+1)}$, $\hat{B}_k \in \mathbb{R}^{(n+1) \times p}$, $\hat{Q}_k \in \mathbb{R}^{(n+1) \times (n+1)}$. The linearized state equation (13) can be written in a form similar to (6), namely as:

$$ \hat{\mathbf{x}}_{k+1} = \hat{A}_k \hat{\mathbf{x}} + \hat{B}_k \mathbf{v}_k \quad k = 1, \dots, M-1. \quad (14) $$

**Enforcing Periodicity of the Parameters.** A main objective of this work is to enforce periodicity of the parameters/controls. For this purpose, let us assume that the length of a time period – measured in time steps – is $T > 0$ and that $M \bmod T = 0$. We now chose the reference trajectory for the control $\mathbf{u}^r = (\mathbf{u}_k^r)_{k=1,...,M-1} \in \mathbb{R}^{(M-1)p}$ to be

$$ \mathbf{u}_k^r := \begin{cases} \mathbf{u}_0, & \text{if } k \le T \\ \mathbf{u}_{k-T}, & \text{if } k > T. \end{cases} \quad (15) $$

Where $\mathbf{u}_0$ is the parameter vector determined by optimization in [1], that was used as an initial guess here. we will enforce periodicity of $\mathbf{u}_k = \mathbf{u}_k^r + \mathbf{v}_k = \mathbf{u}_{k-T} + \mathbf{v}_k$ for $k \ge T+1$.
## 4.3 Application to NPZD Model

From Theorem 1 in section 3, the optimal control is given by

$$v_k = -(R_k + \hat{B}_k^\top \hat{\mathbf{X}}_{k+1} \hat{B}_k)^{-1} \hat{B}_k^\top \hat{\mathbf{X}}_{k+1} \hat{A}_k \hat{\mathbf{x}}_k, \quad k = 1, \dots, M-1.$$

According to (15), we find

$$\mathbf{u}_k = \begin{cases} \mathbf{u}_0 - (R + \hat{B}_k^\top \hat{\mathbf{X}}_{k+1} \hat{B}_k)^{-1} \hat{B}_k^\top \hat{\mathbf{X}}_{k+1} \hat{A}_k \hat{\mathbf{x}}_k, & \text{if } k \le T, \\ \mathbf{u}_{k-T} - (R_k + \tilde{B}_k^\top \hat{\mathbf{X}}_{k+1} \hat{B}_k)^{-1} \hat{B}_k^\top \hat{\mathbf{X}}_{k+1} \hat{A}_k \hat{\mathbf{x}}_k, & \text{if } k > T. \end{cases}$$

Here the $\hat{\mathbf{X}}_k$ can be computed backwards in discrete time, starting from

$$\hat{\mathbf{X}}_M = \hat{Q}_M, \quad (16)$$

as the unique symmetric solutions of the Discrete Riccati equations (8). We set

$$\hat{\mathbf{X}}_k = \begin{bmatrix} \mathbf{X}_k & h_k \\ h_k^\top & \alpha_k \end{bmatrix} \quad (17)$$

with $h_k \in \mathbb{R}^n$ and $\alpha_k \in \mathbb{R}$ for $k = 1, \dots, M-1$. We easily get

$$\mathbf{u}_k = \begin{cases} \mathbf{u}_0 + K_k \mathbf{z}_k + S_k, & \text{if } k \le T, \\ \mathbf{u}_{k-T} + K_k \mathbf{z}_k + S_k, & \text{if } k > T, \end{cases}$$

where $K_k$ and $S_k$ are given by

$$K_k = -(R_k + B_k^\top \mathbf{X}_{k+1} B_k)^{-1} B_k^\top \mathbf{X}_{k+1} A_k, \quad k = M-1, \dots, 1$$

$$S_k = -(R_k + B_k^\top \mathbf{X}_{k+1} B_k)^{-1} B_k^\top (\mathbf{X}_{k+1} r_k + h_{k+1}), \quad k = M-1, \dots, 1.$$

Now, the system (16), (17) to compute the $\mathbf{X}_k$ can be separated into

$$\mathbf{X}_M = Q_M,$$

$$\mathbf{X}_k = Q_k + A_k^\top \mathbf{X}_k A_k - A_k^\top \mathbf{X}_{k+1} B_k (R_k + B_k^\top \mathbf{X}_{k+1} B_k)^{-1} B_k^\top \mathbf{X}_{k+1} A_k, \quad k = M-1, \dots, 1.$$

To evaluate the $\mathbf{X}_k$ and an additional difference equation for the $h_k$, namely

$$h_M = 0,$$

$$h_k = A_k^\top (\mathbf{X}_{k+1} r_k + h_{k+1}) - A_k^\top \mathbf{X}_{k+1} B_k (R_k + B_k^\top \mathbf{X}_{k+1} B_k)^{-1} B_k^\top (\mathbf{X}_{k+1} r_k + h_{k+1}), \quad k = M-1, \dots, 1.$$

For the application on the NPZD Model, $Q_k$ is to be constant for all $k$, this can be written as following

$$Q_k = \mathrm{diag}(\frac{1}{\sigma_l^2})_{l=1,...,5}, \quad k = 1, \dots, M-1,$$

and the matrix $R_k$ can be written as

$$R_k = \mathrm{diag} \left\{ \begin{array}{ll} \frac{1}{|(\mathbf{u}_0)_i|^2}, & i = 1, \dots, p, k = 1, \dots, T \\ \frac{1}{|(\mathbf{u}_{k-T,i})|^2}, & i = 1, \dots, p, k = T+1, \dots, M-1 \end{array} \right.$$
# 5 Optimization Results

## 5.1 Fit of Model Output to Observational Data

This section shows a comparison between the optimized model output obtained by the DLQC method with periodic parameters and the observational data. As a reference we also compare the results to those obtained by a direct optimization of the nonlinear model using *constant parameters* with a *Sequential Quadratic Programming (SQP)* method that takes into account parameter bounds. This method was used in [4]. We performed the optimization for the years 1994 to 1998, in contrast to the years 1991 to 1996 that were used in [4]. The reason for this is that no zooplankton data are available at BATS for the years 1991 to 1993, which would be disadvantageous for the linearization procedure in the DLQC method. In [4] a minimum value of the cost function (5) of $F \approx 70$ was obtained for optimized constant parameters for the five year time interval [1991, 1996]. For the time interval [1994, 1998] a computation with the method used in [4] gave a very similar value. In contrast to these and other (as in [3]) earlier results obtained for constant model parameters, the DLQC method gives a nearly perfect fit of the data. Figure 1 shows the model results $\mathbf{y}$ obtained with the DLQC method together with the observational data $\mathbf{y}^{\text{obs}}$ for the years 1994 to 1998. The model-data fit for $\mathbf{y}_2 = P$ (chlorophyll a) is nearly perfect. Even substantial concentration changes that occur between some neighboring measurement points (e.g. for $\mathbf{y}_4 = P + Z + D$ (particulate organic nitrogen), in 1994, 1995 or 1997) can be captured by the optimized trajectory. There are only some parts of the time interval where the trajectories are slightly farther away from the data, for example in 1996 for zooplankton and in the last two years of the simulated time interval for PON.

## 5.2 Periodicity of the Parameters

we show here that the above model-to-data fit can be achieved by almost annually periodic parameters. This was possible due to an appropriate adjustment of the matrices $Q_k$ and $R_k$, $K = 1, \dots, M-1$, in the cost function (7) used in the DLQC framework, see section 4.3. Thus both the annual periodicity of the parameters. Due to the choice of the reference values for the parameters in the first year, we could also keep the parameters in their desired bounds, although these bounds need not to be imposed explicitly. Figures 2 illustrate the temporal behavior of the selected four parameters that were optimized with the DLQC method. In these figure, the temporal changes of the parameters are plotted against the actual times of the linearization points which are determined by the available measurements. Obviously, the DLQC method then leads to perfectly periodic parameters.

# 6 Conclusion

In this paper, we successfully applied the method linear quadratic optimal control to the optimization of an one-dimensional marine ecosystem model. The model
<figure><img src="image_2.png" /><figcaption>Fig. 1: Observational data $\mathbf{y}_i^{\text{obs}}$, $i = 1, \dots, 4$ and aggregated model trajectories $\mathbf{y}_i$, $i = 1, \dots, 4$, optimized with periodic parameters obtained by the DLQC method. Values are shown for the upper layer (depth less than 5 meters) and years 1994-1998.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 2: Periodicity of the selected optimal parameters $u_{n,1} = \gamma_1, u_{n,4} = \gamma_2, u_{n,6} = g, u_{n,8} = \mu_z$, obtained by the DLQC method.</figcaption></figure>

has to be linearized to fit in the LQOC frame work. The method permits perfect periodic evolution of model parameters and additionally notably improves the fit of the data in comparison with the solution with fixed model parameters. We demonstrated that the LQOC optimization is suitable for the considered prob-
lem and furthermore have shown that this method provides a very reasonable solution. Even with the available small number of observational data, which is typical to oceanographic time series sites, its quality is very high. Temporal deviations of individual parameters about the annual mean can be analyzed further to help making inferences about processes that the model cannot describe well when constant parameters are used. This analysis should contribute to a better understanding of model deficiencies and may improve marine ecosystem models. A next step could be to use only a part of the time horizon to estimate the periodic parameters and verifying the model and the parameters on the remaining part of the data.

## Acknowledgements

The authors would like to thank Johannes Rückelt, Institute of Computer Science, Christian-Albrechts-Universität zu Kiel, for a direct optimization of the nonlinear model using constant parameters with a Sequential Quadratic Programming (SQP) method over the years 1994 to 1998.

## References

1. Oschlies, A., Garon, V.: An eddy-permitting coupled physical-biological model of the North Atlantic 1. Sensitivity to advection numerics and mixed layer physics, Global Biogeochemical Cycles, vol. 13, pp. 135-160, (1999)

2. Schartau, M., Oschlies, A.: Simultaneous data-based optimization of a 1d-ecosystem model at three locations in the north Atlantic: Part I - method and parameter estimates, Journal of Marine Research vol. 61, pp. 765-793, (2003)

3. Ward, B.: Marine Ecosystem Model Analysis Using Data Assimilation, http://web.mit.edu/benw/www/Thesis.pdf, (2009)

4. Rückelt, J., Sauerland, V., Slawig, T., Srivastav, A., Ward, B., Patvardhan, C.: Parameter Optimization and Uncertainty Analysis in a Model of Oceanic CO₂-Uptake using a Hybrid Algorithm and Algorithmic Differentiation, Nonlinear Analysis B Real World Applications, vol. 10, pp. 3993-4009, (2010)

5. Ward, B.A., Anderson, M.A.M., Friedrichs, T.R., Oschlies, A.: Parameter optimization techniques and the problem of underdetermination in marine biogeochemical models, Journal of Marine Systems and Control Letters, vol. 81, pp. 34-43, (2010)

6. Rugh, W. J.: Linear System Theory, Second Edition, Prentice-Hall, Upper Saddle River, New Jersey 07458, (1996)

7. Bermuda Atlantic Time-Series Study, http://www.bios.edu/research/bats.html