Incremental Parameter Estimation under Rank-Deficient Measurement Conditions

Abstract

The computation and modeling of extents has been proposed to handle the complexity of large-scale model identification tasks. Unfortunately, the existing extent-based framework only applies when certain conditions apply. Most typically, it is required that a unique value for each extent can be computed. This severely limits the applicability of this approach. In this work, we propose a novel procedure for parameter estimation inspired by the existing extent-based framework. A key difference with prior work is that the proposed procedure combines structural observability labeling, matrix factorization, and graph-based system partitioning to split the original model parameter estimation problem into parameter estimation problems with the least number of parameters. The value of the proposed method is demonstrated with an extensive simulation study and a study based on a historical data set collected to characterize the isomerization of $\alpha$-pinene. Most importantly, the obtained results indicate that an important barrier to the application of extent-based frameworks for process modeling and monitoring tasks has been lifted.

1. Introduction

Despite advances in model identification theory, parameter estimation can still be very challenging in practice. Such challenges include the lack of identifiability, large computational cost, the need to formulate appropriate experimental designs, and the fact that many methods, such as those for uncertainty analysis, are still being investigated and therefore not standardized. In this work, we focus on a novel method to tackle the computational challenge associated with the identification of kinetic parameters in large dynamic models. Hence, the other challenges associated with parameter identifiability, experimental design, and uncertainty analysis, though relevant, are not investigated here.

To handle the computational challenge, it is typical to devise a protocol for model fitting and model validation. Such protocols can be divided into two broad classes. In the first class, the protocols are domain specific. For instance, several protocols for the identification of activated sludge models are discussed in [1]. Similarly, a protocol for environmental system models is proposed in [2]. These protocols incorporate significant expertise specific to the particular application domain, thereby leading to protocols that are fine-tuned for that domain. While they tend to be similar on a conceptual level, it is rather difficult to apply these protocols universally.

The second class includes protocols that are more general and thus—in principle—broadly applicable. The incremental model identification framework studied in [3,4,5,6,7] is a good example. This method is grounded on the computation of extents, which are—loosely speaking—linear combinations of the original model states and capture the progress of individual dynamic phenomena, such as chemical reactions. A more precise definition will be given below.

The applicability of the extent-based framework is limited to cases where all extents can be computed based on measurements and invariant relationships. However, since this is rarely the case for biological process models, the extent-based framework has been extended in [8] to the case where each extent is either observable or non-sensed—see below for precise definitions. While this recent work provides a meaningful improvement, extent-based incremental model identification remains inapplicable for a wide range of biological scenarios found in practice. For example, this method cannot deal with the frequent case where an extent is sensed but unobservable.

The goal of this study is to present a novel method for incremental parameter identification that is more universally applicable. This method is based on the formulation of a generalized framework for extent computation and the use of a graph-based clustering algorithm. In what follows, we present the method and demonstrate its applicability to cases that could not be handled in an extent-based incremental model identification framework before. The expected impact of the newly developed tools is discussed at the end of this study.

2. Notation and Symbols

The matrix composed of the columns $\mathit{c}$ of the matrix $\mathbf{M}$ is denoted ${\mathbf{M}}_{•,\mathit{c}}$, while the matrix composed of the rows $\mathit{r}$ of the matrix $\mathbf{M}$ is denoted ${\mathbf{M}}_{\mathit{r},•}$. All vectors are column vectors unless mentioned otherwise.

Table 1. List of symbols.

Symbol	Description	Dimensions
${\mathbf{ϵ}}_h$	Measurement error at time $t=t_h$	$M\times 1$
$\mathit{\theta }\left({\mathit{\theta }}^{\left(j\right)}\right)$	Kinetic parameters (in subsystem j)	$T\times 1$$\left(T^{\left(j\right)}\times 1\right)$
${\mathbf{\Lambda }}_{\mathrm{o}}^{\left(j\right)}$, ${\mathbf{\Lambda }}_{\theta }^{\left(j\right)}$	Selection matrix	${\rho }_{\mathrm{o}}^{\left(j\right)}\times {\rho }_{\mathrm{o}}$, $T^{\left(j\right)}\times T$
${\rho }_{\mathrm{o}}\left({\rho }_{\mathrm{o}}^{\left(j\right)},{\rho }_{\mathrm{o},\mathrm{i}}^{\left(j\right)}\right)$	Number of observable extent directions (in system j, interpolated in system j)	$1\times 1$
${\rho }_{\mathrm{u}}$	Number of unobservable extent directions	$1\times 1$
${\rho }_{\mathrm{aug}}$	Number of extents and observable extent directions	$1\times 1$
${\mathbf{\Sigma }}_{\mathbf{ϵ}}$	Measurement error variance-covariance matrix	$M\times M$
$\tau$	Time (integrand)	$1\times 1$
${\Sigma }_{\overline{\chi }}$	Estimation error variance-covariance matrix	$A\times A$
${\chi }_{\mathrm{aug}}$	Extents and observable extent directions	${\rho }_{\mathrm{aug}}\times 1$
${\chi }_{\mathrm{o}}\left({\chi }_{\mathrm{o}}^{\left(j\right)}\right)$	Observable extent directions (in subsystem j)	${\rho }_{\mathrm{o}}\times 1$$\left({\rho }_{\mathrm{o}}^{\left(j\right)}\times 1\right)$
${\overline{\chi }}_{\mathrm{o}}\left({\overline{\chi }}_{\mathrm{o}}^{\left(j\right)}\right)$	Observable extents and extent directions (in subsystem j)	$A\times 1$$\left(A^{\left(j\right)}\times 1\right)$
${\tilde{\overline{\chi }}}_{\mathrm{o},h}\left({\tilde{\overline{\chi }}}_{\mathrm{o},h}^{\left(j\right)}\right)$	Computed observable extents and extent directions (in subsystem j)	$A\times 1$$\left(A^{\left(j\right)}\times 1\right)$
${\tilde{\overline{\chi }}}_{\mathrm{o},\mathrm{i}}\left({\tilde{\overline{\chi }}}_{\mathrm{o},\mathrm{i}}^{\left(j\right)}\right)$	Interpolated observable extents and extent directions (in subsystem j)	$A\times 1$$\left(A^{\left(j\right)}\times 1\right)$
${\chi }_{\mathrm{u}}$	Unobservable extent directions	${\rho }_{\mathrm{u}}\times 1$
${\overline{\chi }}_{\mathrm{u}}$	Unobservable extents and extent directions	$\left(R-A\right)\times 1$
$A\left(A^{\left(j\right)}\right)$	Number of observable extents and extent directions (in subsystem j)	$1\times 1$
$\mathbf{a}$	Indices of ambiguous extents	$R_{\mathrm{a}}\times 1$
$\mathbf{B}$	Reduced row echelon form of $\mathbf{G}$	$M\times R$
${\mathit{d}}_h\left({\mathit{d}}_h^{\left(j\right)}\right)$	Model prediction residuals at $t_h$ (in subsystem j)	$A\times 1\left(A^{\left(j\right)}\times 1\right)$
${\mathbf{e}}^{\left(j\right)}$	Indices of computed extents and extent directions simulated by subsystem j	$A^{\left(j\right)}\times 1$
$\mathcal{F}\left({\mathcal{F}}^{\left(j\right)}\right)$	Information flow graph (for subsystem j)	$--$
$\mathit{f}\left({\mathit{f}}^{\left(j\right)}\right)$	Rate expressions (in subsystem j)	$R\times 1\left(R^{\left(j\right)}\times 1\right)$
$\mathbf{G}$	Extent-based measurement matrix	$M\times R$
${\mathcal{G}}_{\mathrm{o}},{\overline{\mathcal{G}}}_{\mathrm{o}}$	Measurement matrix for the observable extent directions/observable extents and extent directions	$M\times {\rho }_{\mathrm{o}}/A$
H	Number of measurements	$1\times 1$
h	Measurement sample index	$1\times 1$
i	Subsystem index	$1\times 1$
J	Number of subsystems	$1\times 1$
$\mathbf{j}$	Indices of extents in subsystem j	$R^{\left(j\right)}\times 1$
j	Subsystem index	$1\times 1$
$\mathbf{M}$	Species-based measurement matrix	$M\times S$
M	Number of measurement samples	$1\times 1$
$\mathbf{N}$	Stoichiometric matrix	$R\times S$
${\overline{\mathcal{N}}}_{\mathrm{o}},{\overline{\mathcal{N}}}_{\mathrm{u}}$	Extent and extent direction-based stoichiometric matrices	A/$\left(R-A\right)\times S$
$\mathbf{n}$	Indices of non-sensed extents	$R_{\mathrm{n}}\times 1$
$\mathit{n}$$\left({\mathit{n}}_0\right)$	Number of moles (at time $t=0$)	$S\times 1$
$\mathbf{o}$$\left({\mathbf{o}}^{\left(j\right)}\right)$	Indices of observable extents (in subsystem j)	$R_{\mathrm{o}}\times 1$$\left(R_{\mathrm{o}}^{\left(j\right)}\times 1\right)$
$\mathbf{P}$	Projection matrix	$A\times M$
$Q_0,Q_0^{*},Q_1^{\left(j\right)}$	Objective function	$1\times 1$
$R\left(R^{\left(j\right)},R_{\mathrm{a}},R_{\mathrm{n}}\right.$$\left.R_{\mathrm{o}},R_{\mathrm{o}}^{\left(j\right)},R_{\mathrm{o},\mathrm{i}}^{\left(j\right)},R_{\mathrm{s}}\right)$	Number of reactions (in subsystem j, ambiguous/non-sensed/ observable/observable in subsystem j/interpolated in subsystem j/sensed)	$1\times 1$
$\mathit{r}$	Reaction rates	$R\times 1$
r	Reaction index	$1\times 1$
S	Number of species	$1\times 1$
$\mathbf{s}$	Indices of sensed extents	$R_{\mathrm{s}}\times 1$
T ($T^{\left(j\right)}$)	Number of parameters (in subsystem j)	$1\times 1$
$t\left(t_h\right)$	Time (of measurement sample h)	$1\times 1$
${\mathbf{U}}^{\left(j\right)}$, ${\overline{\mathbf{U}}}^{\left(j\right)}$	Mixing matrix	${\rho }_{\mathrm{o}}^{\left(j\right)}/A^{\left(j\right)}\times R^{\left(j\right)}$
${\mathbf{V}}_{\mathrm{o}},{\mathbf{V}}_{\mathrm{u}}$	Direction matrix for observable/unobservable direction directions	$R\times {\rho }_{\mathrm{o}},R\times {\rho }_{\mathrm{u}}$
V	Reactor volume	$1\times 1$
$\mathbf{W}\left({\mathbf{W}}^{\left(j\right)}\right)$	Weight matrix (for subsystem j)	$A\times A$$\left(A^{\left(j\right)}\times A^{\left(j\right)}\right)$
$\mathit{x}\left({\mathit{x}}^{\left(j\right)},x_r\right)$	Extents (in subsystem j, of reaction r)	$R\times 1$$\left(R^{\left(j\right)}\times 1,1\times 1\right)$
$\mathit{y}$$\left({\mathit{y}}_0,\mathit{y}\left(t_h\right)\right)$	Noise-free measurements (at time $t=0$, $t=h$)	$M\times 1$
${\tilde{\mathit{y}}}_h$	Measurements in sample h	$M\times 1$

lists all symbols used in this study.

3. Methods

3.1. System Representation and Extents

3.1.1. Dynamic Model in Terms of Numbers of Moles

We study batch process systems whose dynamic behavior is described by a set of differential equations of the following form:

$$\begin{array}{cccc}\hfill \dot{\mathit{n}}\left(t\right)& ={\mathbf{N}}^{\mathrm{T}}V\mathit{r}\left(t\right),\hfill & \hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\mathit{n}\left(0\right)& ={\mathit{n}}_0.\hfill \end{array}$$

(1)

In the present work, the above equations represent a single-phase reaction system in a vessel with constant volume V. $\mathit{n}\left(t\right)$ is the S-dimensional vector of numbers of moles; ${\mathit{n}}_0$ specifies the initial conditions; $\mathit{r}\left(t\right)$ is the R-dimensional vector of reaction rates; $\mathbf{N}$ is the ($R\times S$)-dimensional stoichiometric matrix.

We assume that M noisy measurements, ${\tilde{\mathit{y}}}_h$, are available at the H sampling times $t_h$ ($h=1,\dots ,H$), according to the following equations:

$$\begin{array}{cccc}\hfill {\tilde{\mathit{y}}}_h& :=\mathit{y}\left(t_h\right)+{\mathbf{ϵ}}_h=\frac{1}{V}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\mathbf{M}\mathit{n}\left(t_h\right)+{\mathbf{ϵ}}_h,\hfill & \hfill {\mathbf{ϵ}}_h& \sim \mathcal{N}\left(\mathbf{0},{\Sigma }_{\mathbf{ϵ}}\right),\hfill \end{array}$$

(2)

where $\mathit{y}\left(t_h\right)$ is the M-dimensional vector of noise-free measurements, ${\mathbf{ϵ}}_h$ the M-dimensional vector of measurement errors, ${\Sigma }_{\mathbf{ϵ}}$ the M-dimensional measurement error variance-covariance matrix, and $\mathbf{M}$ is the $(M\times S)$-dimensional species-based measurement matrix.

We further assume that kinetic laws for the reaction rates are available. These laws are functions of the component masses, $\mathit{n}\left(t\right)$, and a T-dimensional vector of kinetic parameters, $\mathit{\theta }$:

$$\begin{array}{cc}\hfill \mathit{r}\left(t\right)& :=\mathit{f}\left(\mathit{n}\left(t\right),\mathit{\theta }\right).\hfill \end{array}$$

(3)

In the present study, we assume that the rate laws are known, except for the values of the parameters $\mathit{\theta }$ that need to be estimated.

3.1.2. Dynamic Model in Terms of Extents

We now adopt the definition of the extents of reaction, hereafter simply called extents, as the number of times each reaction has occurred since $t=0$. These extents are measured in moles and defined mathematically as:

$$\begin{array}{cc}\hfill \mathit{x}\left(t\right)& :=V{\int }_0^t\mathit{r}\left(\tau \right)d\tau .\hfill \end{array}$$

(4)

It follows that

$$\begin{array}{cc}\hfill \mathit{n}\left(t\right)& ={\mathit{n}}_0+{\mathbf{N}}^{\mathrm{T}}\mathit{x}\left(t\right).\hfill \end{array}$$

(5)

Accordingly, (1)–(3) can be represented equivalently in the following form:

$$\begin{array}{cccc}\hfill \dot{\mathit{x}}\left(t\right)& =V\mathit{f}\left(\mathit{n}\left(t\right),\mathit{\theta }\right),\hfill & \hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\mathit{x}\left(0\right)& =\mathbf{0}\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill {\tilde{\mathit{y}}}_h& =\frac{1}{V}\mathbf{M}\left({\mathit{n}}_0+{\mathbf{N}}^{\mathrm{T}}\mathit{x}\left(t_h\right)\right)+{\mathbf{ϵ}}_h\hfill \end{array}$$

$$\begin{array}{cc}\hfill & ={\mathit{y}}_0+\frac{1}{V}\mathbf{G}\mathit{x}\left(t_h\right)+{\mathbf{ϵ}}_h\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill \mathit{n}\left(t\right)& ={\mathit{n}}_0+{\mathbf{N}}^{\mathrm{T}}\mathit{x}\left(t\right),\hfill \end{array}$$

(8)

where

$$\begin{array}{cc}\hfill {\mathit{y}}_0& :=\frac{1}{V}\mathbf{M}{\mathit{n}}_0\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill \mathbf{G}& :=\mathbf{M}{\mathbf{N}}^{\mathrm{T}},\hfill \end{array}$$

(10)

with the $(M\times R)$-dimensional matrix $\mathbf{G}$ being labeled the extent-based measurement matrix.

Let A denote the rank of $\mathbf{G}$, that is, $A:=\mathrm{rank}\left(\mathbf{G}\right)\le \min \left(M,R\right)$.

3.2. Labeling Extents

3.2.1. Definitions

To label the extents, the following definitions are proposed:

• The rth extent $x_r\left(t\right)$ is labeled sensed if the measurements ${\tilde{\mathit{y}}}_h$ are affected by $x_r\left(t\right)$ through the measurement Equation (7), that is, if ${\mathbf{G}}_{•,r}$ has at least one non-zero element.
• The rth extent is labeled observable if (i) it is sensed (${\mathbf{G}}_{•,r}$ not null); and (ii) the change in the measurements ${\tilde{\mathit{y}}}_h$ caused by a change in $x_r\left(t\right)$ can be unambiguously attributed to the change in that extent, that is, ${\mathbf{G}}_{•,r}$ is independent of all other column vectors in $\mathbf{G}$.

The above labeling is based on the structure of $\mathbf{G}$ and does not depend on the temporal resolution or quality of the recorded measurements (${\tilde{\mathit{y}}}_h$). As such, this means only structural observability is considered. This terminology is similar to the one used in data reconciliation [9,10] as the labels are produced without a dynamic model.

3.2.2. Labeling Procedure

To label the extents, one first computes the $(M\times R)$-dimensional reduced row echelon form of $\mathbf{G}$ denoted $\mathbf{B}$:

$$\begin{array}{cc}\hfill \mathbf{B}& :=\mathrm{rref}\left(\mathbf{G}\right).\hfill \end{array}$$

(11)

This $(M\times R)$-dimensional matrix is composed of A non-zero rows and $M-A$ zero rows. Using $\mathbf{B}$, the extents are labeled with the following procedure:

(a)

Label the $R_{\mathrm{n}}$ extents corresponding to zero columns in $\mathbf{B}$ as non-sensed and use the vector $\mathbf{n}$ to identify their positions in $\mathit{x}$. Label the $R_{\mathrm{s}}$ remaining extents as sensed and use the vector $\mathbf{s}$ to identify their positions in $\mathit{x}$.

(b)

Find all rows in $\mathbf{B}$ with a single non-zero element and find the column positions of these non-zero elements. Label the $R_{\mathrm{o}}$ extents corresponding to these column positions as observable and use the vector $\mathbf{o}$ to identify their positions in $\mathit{x}$. These extents are observable because one can compute a unique value based on the available information (measurements, extent-based measurement matrix, and initial conditions).

(c)

Label the $R_{\mathrm{a}}$ extents that are sensed but not observable as ambiguous. These extents are ambiguous because one cannot compute a unique value based on the available information. The vector $\mathbf{a}$ is used to identify their positions in $\mathit{x}$.

Based on this labeling, the vectors ${\mathit{x}}_{\mathbf{s}}\left(t\right)$, ${\mathit{x}}_{\mathbf{n}}\left(t\right)$, ${\mathit{x}}_{\mathbf{o}}\left(t\right)$, and ${\mathit{x}}_{\mathbf{a}}\left(t\right)$ are defined to represent (i) the $R_{\mathrm{s}}$ sensed extents, (ii) the $R_{\mathrm{n}}=R-R_{\mathrm{s}}$ non-sensed extents; (iii) the $R_{\mathrm{o}}$ observable extents; and (iv) the $R_{\mathrm{a}}=R_{\mathrm{s}}-R_{\mathrm{o}}$ ambiguous extents. This is illustrated in Figure 1.

With the above definitions, (5) can be reformulated as:

$$\begin{array}{cc}\hfill \mathit{n}\left(t\right)& ={\mathit{n}}_0+{{\mathbf{N}}_{\mathbf{s},•}}^{\mathrm{T}}{\mathit{x}}_{\mathbf{s}}\left(t\right)+{{\mathbf{N}}_{\mathbf{n},•}}^{\mathrm{T}}{\mathit{x}}_{\mathbf{n}}\left(t\right)\hfill \\ \hfill & ={\mathit{n}}_0+{{\mathbf{N}}_{\mathbf{o},•}}^{\mathrm{T}}{\mathit{x}}_{\mathbf{o}}\left(t\right)+{{\mathbf{N}}_{\mathbf{a},•}}^{\mathrm{T}}{\mathit{x}}_{\mathbf{a}}\left(t\right)+{{\mathbf{N}}_{\mathbf{n},•}}^{\mathrm{T}}{\mathit{x}}_{\mathbf{n}}\left(t\right).\hfill \end{array}$$

(12)

Furthermore, it follows from ${\mathbf{G}}_{•,\mathbf{n}}={\mathbf{0}}_{M\times R_{\mathrm{n}}}$ that

$$\begin{array}{cc}\hfill {\tilde{\mathit{y}}}_h& ={\mathit{y}}_0+\frac{1}{V}{\mathbf{G}}_{•,\mathbf{s}}{\mathit{x}}_{\mathbf{s}}\left(t_h\right)+{\mathbf{ϵ}}_h\hfill \\ & ={\mathit{y}}_0+\frac{1}{V}{\mathbf{G}}_{•,\mathbf{o}}{\mathit{x}}_{\mathbf{o}}\left(t_h\right)+\frac{1}{V}{\mathbf{G}}_{•,\mathbf{a}}{\mathit{x}}_{\mathbf{a}}\left(t_h\right)+{\mathbf{ϵ}}_h.\hfill \end{array}$$

(13)

3.2.3. Practical Cases

Depending on the values of A, M, and R, one can distinguish the following cases:

• Full-rank extent-based measurement matrix ($A=R$). This occurs when there are at least as many measurements as reactions ($M\ge R$) and the matrix $\mathbf{G}$ is full column rank. As a result, $R_{\mathrm{n}}=0$, $R_{\mathrm{o}}=R$, and $R_{\mathrm{a}}=0$. This is the most frequently studied case, e.g., in [5,6,7]. This case enables computing unique values for the R extents by means of a linear transformation [5], thus allowing the estimation of kinetic parameters for each reaction rate model individually.
• Rank-deficient measurement matrix ($A<R$). If $A<R$, for example because there are fewer measurements than reactions ($M<R$), it is no longer possible to compute all R extents from M measurements without additional information such as an established kinetic model. One can distinguish two situations within this case:

  (a)

  No ambiguity ($R_{\mathrm{a}}=0$). In this case, $R_{\mathrm{o}}=A$ and $R_{\mathrm{o}}+R_{\mathrm{n}}=R$. As shown in [8], it is possible in this case to implement efficient parameter estimation by identifying subsystems of the complete model that include a subset of the kinetic rate laws and their parameters.

  (b)

  Ambiguity present ($0<R_{\mathrm{a}}\le A$). This situation results in $R_{\mathrm{o}}<A$. For this case, no generally applicable method for incremental parameter estimation is available until now.

The method proposed in this work is developed with the aim of handling all the above cases in a single framework for extent-based kinetic parameter estimation.

3.3. Observable and Unobservable Extent Directions

The ambiguous extents are now investigated in more detail to determine observable directions among them.

3.3.1. Factorization of ${\mathbf{G}}_{•,\mathbf{a}}$

Let ${\rho }_{\mathrm{o}}$ be the rank of the $(M\times R_{\mathrm{a}})$-dimensional measurement matrix ${\mathbf{G}}_{•,\mathbf{a}}$. Then, ${\mathbf{G}}_{•,\mathbf{a}}$ can be factorized into the $(M\times {\rho }_{\mathrm{o}})$-dimensional measurement matrix ${\mathcal{G}}_{\mathrm{o}}$ and the $(R_{\mathrm{a}}\times {\rho }_{\mathrm{o}})$-dimensional matrix ${\mathbf{V}}_{\mathrm{o}}$:

$$\begin{array}{cc}\hfill {\mathbf{G}}_{•,\mathbf{a}}& ={\mathcal{G}}_{\mathrm{o}}{{\mathbf{V}}_{\mathrm{o}}}^{\mathrm{T}}.\hfill \end{array}$$

(14)

Using the reduced row echelon form ${\mathbf{B}}_{•,\mathbf{a}}$, the matrices ${\mathcal{G}}_{\mathrm{o}}$ and ${\mathbf{V}}_{\mathrm{o}}$ are computed as follows:

• ${{\mathbf{V}}_{\mathrm{o}}}^{\mathrm{T}}$ is obtained by selecting the ${\rho }_{\mathrm{o}}$ non-zero rows in ${\mathbf{B}}_{•,\mathbf{a}}$,
• ${\mathcal{G}}_{\mathrm{o}}$ is the matrix composed of the columns of ${\mathbf{G}}_{•,\mathbf{a}}$ corresponding to the column positions of the first non-zero elements in the rows of ${{\mathbf{V}}_{\mathrm{o}}}^{\mathrm{T}}$.

3.3.2. Definition of Observable and Unobservable Extent Directions

The term ${\mathbf{G}}_{•,\mathbf{a}}{\mathit{x}}_{\mathbf{a}}\left(t_h\right)$ in (13) can be expressed as:

$$\begin{array}{c}\hfill {\mathbf{G}}_{•,\mathbf{a}}{\mathit{x}}_{\mathbf{a}}\left(t\right)={\mathcal{G}}_{\mathrm{o}}{{\mathbf{V}}_{\mathrm{o}}}^{\mathrm{T}}{\mathit{x}}_{\mathbf{a}}\left(t\right)={\mathcal{G}}_{\mathrm{o}}{\chi }_{\mathrm{o}}\left(t\right),\end{array}$$

(15)

with ${\chi }_{\mathrm{o}}\left(t\right)$ the ${\rho }_{\mathrm{o}}$-dimensional vector of observable extent directions among the ambiguous extents:

$$\begin{array}{cc}\hfill {\chi }_{\mathrm{o}}\left(t\right)& :={{\mathbf{V}}_{\mathrm{o}}}^{\mathrm{T}}{\mathit{x}}_{\mathbf{a}}\left(t\right)={{\mathbf{V}}_{\mathrm{o}}}^{\mathrm{T}}{\mathbf{I}}_{\mathbf{a},•}\mathit{x}\left(t\right),\hfill \end{array}$$

(16)

where the $(R_{\mathrm{a}}\times R)$-dimensional matrix ${\mathbf{I}}_{\mathbf{a},•}$ includes the rows $\mathbf{a}$ of the identity matrix ${\mathbf{I}}_R$ such that

$$\begin{array}{cc}\hfill {\mathit{x}}_{\mathbf{a}}\left(t\right)& ={\mathbf{I}}_{\mathbf{a},•}\mathit{x}\left(t\right).\hfill \end{array}$$

(17)

Equation (15) indicates that while the extents ${\mathit{x}}_{\mathbf{a}}\left(t\right)$ cannot be observed individually, their combined effects on the measurements can be observed as the linear combinations ${\chi }_{\mathrm{o}}\left(t\right)$.

Remark 1.

The unobservable extent directions span the null space of ${\mathbf{G}}_{•,\mathbf{a}}$, which is also the null space of ${{\mathbf{V}}_{\mathrm{o}}}^{\mathrm{T}}$. Denoting this null space by the $(R_{\mathrm{a}}\times {\rho }_{\mathrm{u}})$-dimensional matrix ${\mathbf{V}}_{\mathrm{u}}$,

$$\begin{array}{cc}\hfill {\mathbf{V}}_{\mathrm{u}}& =\mathrm{null}\left({{\mathbf{V}}_{\mathrm{o}}}^{\mathrm{T}}\right),\hfill \end{array}$$

(18)

with ${\rho }_{\mathrm{u}}=(R_{\mathrm{a}}-{\rho }_{\mathrm{o}})$, one can define the ${\rho }_{\mathrm{u}}$-dimensional vector of unobservable extent directions${\chi }_{\mathrm{u}}\left(t\right)$ as:

$$\begin{array}{cc}\hfill {\chi }_{\mathrm{u}}\left(t\right)& :={{\mathbf{V}}_{\mathrm{u}}}^{\mathrm{T}}{\mathit{x}}_{\mathbf{a}}\left(t\right)={{\mathbf{V}}_{\mathrm{u}}}^{\mathrm{T}}{\mathbf{I}}_{\mathbf{a},•}\mathit{x}\left(t\right).\hfill \end{array}$$

(19)

3.4. Observable Extents and Extent Directions

We further define the vector ${\overline{\chi }}_{\mathrm{o}}\left(t\right)$ consisting of the $A=R_{\mathrm{o}}+{\rho }_{\mathrm{o}}$ observable extents and extent directions as follows:

$$\begin{array}{cc}\hfill {\overline{\chi }}_{\mathrm{o}}\left(t\right)& :=\left[\begin{array}{c}{\mathit{x}}_{\mathbf{o}}\left(t\right)\\ {\chi }_{\mathrm{o}}\left(t\right)\end{array}\right]=\left[\begin{array}{c}{\mathit{x}}_{\mathbf{o}}\left(t\right)\\ {{\mathbf{V}}_{\mathrm{o}}}^{\mathrm{T}}{\mathit{x}}_{\mathbf{a}}\left(t\right)\end{array}\right].\hfill \end{array}$$

(20)

With this definition, the measurement Equation (13) can be rewritten as:

$$\begin{array}{cc}\hfill {\tilde{\mathit{y}}}_h& ={\mathit{y}}_0+\frac{1}{V}{\mathbf{G}}_{•,\mathbf{o}}{\mathit{x}}_{\mathbf{o}}\left(t_h\right)+\frac{1}{V}{\mathcal{G}}_{\mathrm{o}}{\chi }_{\mathrm{o}}\left(t\right)+{\mathbf{ϵ}}_h\hfill \\ & ={\mathit{y}}_0+\frac{1}{V}{\overline{\mathcal{G}}}_{\mathrm{o}}{\overline{\chi }}_{\mathrm{o}}\left(t_h\right)+{\mathbf{ϵ}}_h,\hfill \end{array}$$

(21)

with ${\overline{\mathcal{G}}}_{\mathrm{o}}$, an $(M\times A)$-dimensional matrix, constructed as

$$\begin{array}{cc}\hfill {\overline{\mathcal{G}}}_{\mathrm{o}}& :=\left[\begin{array}{c}{\mathbf{G}}_{•,\mathbf{o}}{\mathcal{G}}_{\mathrm{o}}\end{array}\right].\hfill \end{array}$$

(22)

Since the $(R_{\mathrm{a}}\times A)$-dimensional matrix ${\overline{\mathcal{G}}}_{\mathrm{o}}$ has full column rank, (21) can be used to compute the maximum-likelihood estimates ${\tilde{\overline{\chi }}}_{\mathrm{o},h}$ of the observable extents and extent directions directly from the measurements:

$$\begin{array}{cc}\hfill {\tilde{\overline{\chi }}}_{\mathrm{o},h}& =\mathbf{P}\left({\tilde{\mathit{y}}}_h-{\mathit{y}}_0\right),\hfill \end{array}$$

(23)

with the $(A\times M)$-dimensional matrix $\mathbf{P}$ given by:

$$\begin{array}{cc}\hfill \mathbf{P}& =V{\left({{\overline{\mathcal{G}}}_{\mathrm{o}}}^{\mathrm{T}}{\Sigma }_{\mathbf{ϵ}}^{-1}{\overline{\mathcal{G}}}_{\mathrm{o}}\right)}^{-1}{{\overline{\mathcal{G}}}_{\mathrm{o}}}^{\mathrm{T}}{\Sigma }_{\mathbf{ϵ}}^{-1}.\hfill \end{array}$$

(24)

The associated expected variance-covariance matrix of the estimation errors becomes:

$$\begin{array}{cc}\hfill {\Sigma }_{\overline{\chi }}& =\mathbf{P}{\Sigma }_{\mathbf{ϵ}}{\mathbf{P}}^{\mathrm{T}}={\left({{\overline{\mathcal{G}}}_{\mathrm{o}}}^{\mathrm{T}}{\Sigma }_{\mathbf{ϵ}}^{-1}{\overline{\mathcal{G}}}_{\mathrm{o}}\right)}^{-1}.\hfill \end{array}$$

(25)

3.5. Unobservable Extents and Extent Directions

We finally define the vector ${\overline{\chi }}_{\mathrm{u}}\left(t\right)$ consisting of the $R-A=R_{\mathrm{n}}+{\rho }_{\mathrm{u}}$ unobservable extents and non-sensed extent directions as follows:

$$\begin{array}{cc}\hfill {\overline{\chi }}_{\mathrm{u}}\left(t\right)& :=\left[\begin{array}{c}{\mathit{x}}_{\mathbf{n}}\left(t\right)\\ {\chi }_{\mathrm{u}}\left(t\right)\end{array}\right]=\left[\begin{array}{c}{\mathit{x}}_{\mathbf{n}}\left(t\right)\\ {{\mathbf{V}}_{\mathrm{u}}}^{\mathrm{T}}{\mathit{x}}_{\mathbf{a}}\left(t\right)\end{array}\right]=\left[\begin{array}{c}{\mathbf{I}}_{\mathbf{n},•}\\ {{\mathbf{V}}_{\mathrm{u}}}^{\mathrm{T}}{\mathbf{I}}_{\mathbf{a},•}\end{array}\right]\mathit{x}\left(t\right).\hfill \end{array}$$

(26)

With this definition, the expression for the number of moles (12) can be rewritten as:

$$\begin{array}{cc}\hfill \mathit{n}\left(t\right)& ={\mathit{n}}_0+{{\mathbf{N}}_{\mathbf{o},•}}^{\mathrm{T}}{\mathit{x}}_{\mathbf{o}}\left(t\right)+{{\mathbf{N}}_{\mathbf{a},•}}^{\mathrm{T}}\left({{\left({\mathbf{V}}_{\mathrm{o}}\right)}^{+}}^{\mathrm{T}}{\chi }_{\mathrm{o}}\left(t\right)+{{\left({\mathbf{V}}_{\mathrm{u}}\right)}^{+}}^{\mathrm{T}}{\chi }_{\mathrm{u}}\left(t\right)\right)+{{\mathbf{N}}_{\mathbf{n},•}}^{\mathrm{T}}{\mathit{x}}_{\mathbf{n}}\left(t\right)\hfill \\ & ={\mathit{n}}_0+{{\overline{\mathcal{N}}}_{\mathrm{o}}}^{\mathrm{T}}{\overline{\chi }}_{\mathrm{o}}\left(t\right)+{{\overline{\mathcal{N}}}_{\mathrm{u}}}^{\mathrm{T}}{\overline{\chi }}_{\mathrm{u}}\left(t\right)\hfill \end{array}$$

(27)

with:

$$\begin{array}{cc}\hfill {\overline{\mathcal{N}}}_{\mathrm{o}}& :=\left[\begin{array}{c}{\mathbf{N}}_{\mathbf{o},•}\\ {\left({\mathbf{V}}_{\mathrm{o}}\right)}^{+}{\mathbf{N}}_{\mathbf{a},•}\end{array}\right]\hfill \end{array}$$

(28)

$$\begin{array}{cc}\hfill {\overline{\mathcal{N}}}_{\mathrm{u}}& :=\left[\begin{array}{c}{\mathbf{N}}_{\mathbf{n},•}\\ {\left({\mathbf{V}}_{\mathrm{u}}\right)}^{+}{\mathbf{N}}_{\mathbf{a},•}\end{array}\right]\hfill \end{array}$$

(29)

3.6. System Partitioning

Incremental parameter estimation is based on the possibility of separating the parameter estimation problem into smaller problems. Ideally, if all extents of reaction could be computed from measurements, each reaction could be identified individually, that is, independently of the other reactions. This way, the parameter estimation problem reduces to the solution of R smaller problems. As discussed in the previous section, some extents may not be observable. For these situations, it would still be nice to be able to separate the parameter estimation problem into J smaller problems ($J\le R$). The objective is therefore to partition the reaction system effectively—with as many small groups of reactions as possible—to simplify the parameter estimation task. The first step to achieve this consists of system partitioning.

An algorithm is developed to group the kinetic parameters into J parameter subsets ($j=1,\dots ,J$), each represented as a $T^{\left(j\right)}$-dimensional vector ${\mathit{\theta }}^{\left(j\right)}$ satisfying the following properties:

• The size $T^{\left(j\right)}$ of each parameter subset should be as small as possible.
• The estimates ${\widehat{\mathit{\theta }}}^{\left(j\right)}$ in the jth parameter subset can be computed without consideration of any other parameter subset ${\mathit{\theta }}^{\left(i\right)}$, $i\ne j$.
• Each parameter in $\mathit{\theta }$ appears in at most one of the parameter subsets ${\mathit{\theta }}^{\left(j\right)}$.

This objective is achieved by means of model reformulation and graph-based system partitioning. The graph-based procedure can be interpreted as a symbolic manipulation of the process model. It does not require symbolic differentiation, however.

3.6.1. Step 1—Model Reformulation

An extended model is first defined to describe the dynamics of all extents and all observable directions among the ambiguous extents. To this end, the following procedure is applied:

(a) 

Express $\dot{\mathit{x}}\left(t\right)$ as a function of ${\overline{\chi }}_{\mathrm{o}}\left(t\right)$ and $\mathit{x}\left(t\right)$

The dynamic model (6) is modified by replacing $\mathit{n}\left(t\right)$ with the right-hand side of (27):

$$\begin{array}{ccc}\hfill \dot{\mathit{x}}\left(t\right)& =V\mathit{f}\left({\overline{\chi }}_{\mathrm{o}}\left(t\right),{\overline{\chi }}_{\mathrm{u}}\left(t\right),\mathit{\theta }\right),\hfill & \hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\mathit{x}\left(0\right)=\mathbf{0}.\end{array}$$

(30)

The vector ${\overline{\chi }}_{\mathrm{u}}\left(t\right)$ is now replaced with the right-hand side of (26). As a result, the above system becomes:

$$\begin{array}{ccc}\hfill \dot{\mathit{x}}\left(t\right)& =V\mathit{f}\left({\overline{\chi }}_{\mathrm{o}}\left(t\right),\mathit{x}\left(t\right),\mathit{\theta }\right),\hfill & \hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\mathit{x}\left(0\right)=\mathbf{0}.\end{array}$$

(31)

(b) 

State augmentation

Define the ${\rho }_{\mathrm{aug}}$-dimensional vector ${\chi }_{\mathrm{aug}}\left(t\right):=\left[\begin{array}{c}\mathit{x}\left(t\right)\\ {\chi }_{\mathrm{o}}\left(t\right)\end{array}\right]=\left[\begin{array}{c}{\mathbf{I}}_R\\ {{\mathbf{V}}_{\mathrm{o}}}^{\mathrm{T}}{\mathbf{I}}_{\mathbf{a},•}\end{array}\right]\mathit{x}\left(t\right)$ that includes all extents and extent directions, with ${\rho }_{\mathrm{aug}}=R+{\rho }_{\mathrm{o}}$. The dynamic behavior of ${\chi }_{\mathrm{aug}}\left(t\right)$ can be described by a differential-algebraic system including the R differential Equation (31) and the ${\rho }_{\mathrm{o}}$ algebraic expressions (16):

$$\begin{array}{ccc}\hfill \dot{\mathit{x}}\left(t\right)& =V\mathit{f}\left({\overline{\chi }}_{\mathrm{o}}\left(t\right),\mathit{x}\left(t\right),\mathit{\theta }\right),\hfill & \hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\mathit{x}\left(0\right)=\mathbf{0}\end{array}$$

(32)

$$\begin{array}{cc}\hfill {\chi }_{\mathrm{o}}\left(t\right)& ={{\mathbf{V}}_{\mathrm{o}}}^{\mathrm{T}}{\mathbf{I}}_{\mathbf{a},•}\mathit{x}\left(t\right).\hfill \end{array}$$

(33)

(c) 

Interpolation of the observable extents and extent directions

To increase the efficiency of system partitioning, it is useful to account for the fact that the observable extents and extent directions can be expressed in terms of measurements. However, since the observable extents and extent directions are only known at H discrete measurement points, their values always need to be obtained via interpolation. In this work, we apply piece-wise linear interpolation as follows:

$$\begin{array}{cccccc}\hfill {\tilde{\overline{\chi }}}_{\mathrm{o},\mathrm{i}}\left(t\right)& :={\tilde{\overline{\chi }}}_{\mathrm{o},h}+\frac{t-t_h}{t_{h+1}-t_h}\left({\tilde{\overline{\chi }}}_{\mathbf{o},h+1}-{\tilde{\overline{\chi }}}_{\mathrm{o},h}\right),\hfill & \hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{t}_h& \le t<t_{h+1},\hfill & \hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}h& =1,\dots ,H,\hfill \end{array}$$

(34)

with which the system (32) and (33) becomes:

$$\begin{array}{ccc}\hfill \dot{\mathit{x}}\left(t\right)& =V\mathit{f}\left({\tilde{\overline{\chi }}}_{\mathrm{o},\mathrm{i}}\left(t\right),\mathit{x}\left(t\right),\theta \right),\hfill & \hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\mathit{x}\left(0\right)=\mathbf{0}\end{array}$$

(35)

$$\begin{array}{cc}\hfill {\chi }_{\mathrm{o}}\left(t\right)& ={{\mathbf{V}}_{\mathrm{o}}}^{\mathrm{T}}{\mathbf{I}}_{\mathbf{a},•}\mathit{x}\left(t\right).\hfill \end{array}$$

(36)

3.6.2. Step 2—Graph-Based System Partitioning

The equation system (35) and (36) is now analyzed by means of a graph partitioning procedure to determine the smallest groups of kinetic parameters that can be estimated separately. To this end, the following steps are performed:

(a) 

Create a graph

One creates a directed graph $\mathcal{F}$ with a vertex for every state variable in ${\chi }_{\mathrm{aug}}\left(t\right)$ and every parameter in $\mathit{\theta }$. Hence, this graph has $R+{\rho }_{\mathrm{o}}+T$ vertices. A directed arc is added from vertex v to vertex w if the vth element of $\left[\begin{array}{c}{\chi }_{\mathrm{aug}}\left(t\right)\\ \mathit{\theta }\end{array}\right]$ appears in the right-hand side of the wth equation in (35) and (36) ($v=1,\dots ,R+{\rho }_{\mathrm{o}}+T$, $w=1,\dots ,R+{\rho }_{\mathrm{o}}$). This graph represents the information flow for simulating (35) and (36). Additional arcs and vertices may be added to describe the influence of known inputs and the links between extents and measured variables. For system partitioning, this is however unnecessary and omitted for clarity.

(b) 

Extents predicted from measurements or simulation

The simultaneous approach uses a complete model of the reaction system to predict the extents (or concentrations) via simulation. If one wants to partition the reaction system into small groups of reactions, only the extents belonging to a given group can be generated via the simulation of that group. The other extents that enter the rate laws must be provided by the user as quantities known from measurements.

That information can be included in the graph $\mathcal{F}$ by annotating the various arcs. The arcs that originate at a vertex corresponding to an observable extent or an observable extent direction are labeled observation arcs, considering that observable extent or an observable extent direction can be replaced with their measured values (34). They are visualized as dashed-line arrows. The remaining arcs are labeled simulation arcs and visualized as solid-line arrows. The observation arcs represent the idea that the elements of ${\tilde{\overline{\chi }}}_{\mathrm{o},\mathrm{i}}\left(t\right)$ can be regarded as known inputs for simulating (35) and (36).

(c) 

Subgraph selection

Identify the J subgraphs ${\mathcal{F}}^{\left(j\right)}$ consisting of arcs and vertices in $\mathcal{F}$ on directed paths that (i) lead to a vertex representing an observable extent or an observable extent direction; and (ii) consist of simulation arcs only. The selected vertices represent an $R^{\left(j\right)}$-dimensional vector of extents ${\mathit{x}}^{\left(j\right)}\left(t\right)$, a ${\rho }_{\mathrm{o}}^{\left(j\right)}$-dimensional vector of directions ${\chi }_{\mathrm{o}}^{\left(j\right)}\left(t\right)$, and a $T^{\left(j\right)}$-dimensional vector of parameters ${\mathit{\theta }}^{\left(j\right)}$. The positions of ${\mathit{x}}^{\left(j\right)}\left(t\right)$ in $\mathit{x}\left(t\right)$ are given by the vector $\mathbf{j}$ so that:

$$\begin{array}{cc}\hfill {\mathit{x}}^{\left(j\right)}\left(t\right)& :={\mathbf{I}}_{\mathbf{j},•}\mathit{x}\left(t\right)\hfill \end{array}$$

(37)

$$\begin{array}{cc}\hfill {\mathit{f}}^{\left(j\right)}\left(•\right)& :={\mathbf{I}}_{\mathbf{j},•}\mathit{f}\left(•\right),\hfill \end{array}$$

(38)

and the selection matrices ${\Lambda }_{\mathrm{o}}^{\left(j\right)}$ and ${\Lambda }_{\theta }^{\left(j\right)}$ are defined so that:

$$\begin{array}{cc}\hfill {\chi }_{\mathrm{o}}^{\left(j\right)}\left(t\right)& :={\Lambda }_{\mathrm{o}}^{\left(j\right)}{\chi }_{\mathrm{o}}\left(t\right)\hfill \end{array}$$

(39)

$$\begin{array}{cc}\hfill {\mathit{\theta }}^{\left(j\right)}& :={\Lambda }_{\theta }^{\left(j\right)}\mathit{\theta }.\hfill \end{array}$$

(40)

This means that each subgraph ${\mathcal{F}}^{\left(j\right)}$ represents a subset of Equations (35) and (36) that describes the dynamics of ${\mathit{x}}^{\left(j\right)}\left(t\right)$ and ${\chi }_{\mathrm{o}}^{\left(j\right)}\left(t\right)$ without reference to any other state variable:

$$\begin{array}{ccc}\hfill {\dot{\mathit{x}}}^{\left(j\right)}\left(t\right)& =V{\mathit{f}}^{\left(j\right)}\left({\tilde{\overline{\chi }}}_{\mathrm{o},\mathrm{i}}\left(t\right),{\mathit{x}}^{\left(j\right)}\left(t\right),{\mathit{\theta }}^{\left(j\right)}\right),\hfill & \hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{\mathit{x}}^{\left(j\right)}\left(0\right)={\mathit{x}}_{\mathbf{o}}\end{array}$$

(41)

$$\begin{array}{cc}\hfill {\chi }_{\mathrm{o}}^{\left(j\right)}\left(t\right)& ={\mathbf{U}}^{\left(j\right)}{\mathit{x}}^{\left(j\right)}\left(t\right),\hfill \end{array}$$

(42)

with ${\mathbf{U}}^{\left(j\right)}:={\Lambda }_{\mathrm{o}}^{\left(j\right)}{{\mathbf{V}}_{\mathrm{o}}}^{\mathrm{T}}{\mathbf{I}}_{\mathbf{a},•}{{\mathbf{I}}_{\mathbf{j},•}}^{\mathrm{T}}$.

(d) 

Add observation arcs and vertices

For every graph ${\mathcal{F}}^{\left(j\right)}$, add (i) the observation arcs that have a target vertex belonging to ${\mathcal{F}}^{\left(j\right)}$ and (ii) the source vertices of the added observation arcs. These added source vertices represent the minimal subset of interpolants in ${\tilde{\overline{\chi }}}_{\mathrm{o},\mathrm{i}}\left(t\right)$ (34) that are required to simulate ${\mathit{x}}^{\left(j\right)}\left(t\right)$ and ${\chi }_{\mathrm{o}}^{\left(j\right)}\left(t\right)$ and are referred to as ${\tilde{\overline{\chi }}}_{\mathrm{o},\mathrm{i}}^{\left(j\right)}\left(t\right)$. This means that the graph ${\mathcal{F}}^{\left(j\right)}$ now represents all information required to simulate the observable extents ${\mathit{x}}^{\left(j\right)}\left(t\right)$ and the observable extent directions ${\chi }_{\mathrm{o}}^{\left(j\right)}\left(t\right)$. Accordingly, one can rewrite the jth equation subsystem as:

$$\begin{array}{ccc}\hfill {\dot{\mathit{x}}}^{\left(j\right)}\left(t\right)& =V{\mathit{f}}^{\left(j\right)}\left({\tilde{\overline{\chi }}}_{\mathrm{o},\mathrm{i}}^{\left(j\right)}\left(t\right),{\chi }_{\mathrm{o}}^{\left(j\right)}\left(t\right),{\mathit{x}}^{\left(j\right)}\left(t\right),{\mathit{\theta }}^{\left(j\right)}\right),\hfill & \hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{\mathit{x}}^{\left(j\right)}\left(0\right)={\mathit{x}}_{\mathbf{o}}\end{array}$$

(43)

$$\begin{array}{cc}\hfill {\chi }_{\mathrm{o}}^{\left(j\right)}\left(t\right)& ={\mathbf{U}}^{\left(j\right)}{\mathit{x}}^{\left(j\right)}\left(t\right).\hfill \end{array}$$

(44)

At the end of this procedure, the equation system (35) and (36) is approximated by J smaller equation subsystems, each including a subset of the kinetic parameters. As intended, every kinetic parameter appears in at most one of the J subsystems. In addition, the identified subsystems do not share any of the observable extents or extent directions as state variables, that is, each observable extent and extent direction is simulated in only one of the identified subsystems.

3.7. Parameter Estimation Methods

In this work, we solve the parameter estimation problem in two distinct ways. The first way consists in a simultaneous estimation of all parameters in the maximum-likelihood sense. The second way consists in an extent-based incremental parameter estimation. The next paragraphs describe how incremental parameter estimation approximates the simultaneous estimation procedure to minimize the number of parameters that are estimated together.

3.7.1. Simultaneous Parameter Estimation

Maximum-likelihood estimation of the kinetic parameters can be obtained by solving minimization of the weighted mean squared error:

Problem $P_0$

$$\begin{array}{c}\hfill \widehat{\mathit{\theta }}=arg\underset{\mathit{\theta }}{min}{Q}_0:=\frac{1}{H·M}\sum _h^H{\left({\tilde{\mathit{y}}}_h-\mathit{y}\left(t_h\right)\right)}^{\mathrm{T}}{\Sigma }_{\mathbf{ϵ}}^{-1}\left({\tilde{\mathit{y}}}_h-\mathit{y}\left(t_h\right)\right),\end{array}$$

(45)

subject to (1)–(3). This estimation problem can be equivalently formulated in terms of the computed extents as follows:

Problem $P_0^{*}$

$$\begin{array}{c}\hfill \widehat{\mathit{\theta }}=arg\underset{\mathit{\theta }}{min}{Q}_0^{*}:=\sum _h^H{\left({\tilde{\overline{\chi }}}_{\mathrm{o},h}-{\overline{\chi }}_{\mathrm{o}}\left(t_h\right)\right)}^{\mathrm{T}}\mathbf{W}\left({\tilde{\overline{\chi }}}_{\mathrm{o},h}-{\overline{\chi }}_{\mathrm{o}}\left(t_h\right)\right),\end{array}$$

(46)

subject to (6)–(8) and (20) and with $\mathbf{W}:={\Sigma }_{\overline{\chi }}^{-1}$.

3.7.2. Incremental Parameter Estimation

The incremental parameter estimation procedure is obtained by applying two modifications, A and B, to problem $P_0^{*}$.

(a) 

Modification A: Removing correlation terms. Let $A^{\left(j\right)}=R_{\mathrm{o}}^{\left(j\right)}+{\rho }_{\mathrm{o}}^{\left(j\right)}$ and define the $A^{\left(j\right)}$-dimensional vector of observable extents and extent directions ${\overline{\chi }}_{\mathrm{o}}^{\left(j\right)}\left(t\right)$. This vector includes all observable extents ${\mathit{x}}_{\mathbf{o}}^{\left(j\right)}\left(t\right)$, whose positions in ${\mathit{x}}^{\left(j\right)}\left(t\right)$ are given by the vector ${\mathbf{o}}^{\left(j\right)}$, and all observable extent directions ${\chi }_{\mathrm{o}}^{\left(j\right)}\left(t\right)$ in Subsystem j. We further define the matrix ${\overline{\mathbf{U}}}^{\left(j\right)}$ so that:

$$\begin{array}{cc}\hfill {\mathit{x}}_{\mathbf{o}}^{\left(j\right)}\left(t\right)& :={\mathbf{I}}_{{\mathbf{o}}^{\left(j\right)},•}{\mathit{x}}^{\left(j\right)}\left(t\right)\hfill \end{array}$$

(47)

$$\begin{array}{cc}\hfill {\overline{\chi }}_{\mathrm{o}}^{\left(j\right)}\left(t\right)& =\left[\begin{array}{c}{\mathit{x}}_{\mathbf{o}}^{\left(j\right)}\left(t\right)\\ {\chi }_{\mathrm{o}}^{\left(j\right)}\left(t\right)\end{array}\right]={\overline{\mathbf{U}}}^{\left(j\right)}{\mathit{x}}^{\left(j\right)}\left(t\right)\hfill \end{array}$$

(48)

$$\begin{array}{cc}\hfill {\overline{\mathbf{U}}}^{\left(j\right)}& :=\left[\begin{array}{c}{\mathbf{I}}_{{\mathbf{o}}^{\left(j\right)},•}\\ {\mathbf{U}}^{\left(j\right)}\end{array}\right].\hfill \end{array}$$

(49)

The $A^{\left(j\right)}$-dimensional vector ${\tilde{\overline{\chi }}}_{\mathrm{o},h}^{\left(j\right)}$ is then defined by selecting the elements of ${\tilde{\overline{\chi }}}_{\mathrm{o},h}$ as:

$$\begin{array}{cc}\hfill {\tilde{\overline{\chi }}}_{\mathrm{o},h}^{\left(j\right)}& :={\mathbf{I}}_{{\mathbf{e}}^{\left(j\right)},•}{\tilde{\overline{\chi }}}_{\mathrm{o},h},\hfill \end{array}$$

(50)

where the vector ${\mathbf{e}}^{\left(j\right)}$ gives the positions of ${\tilde{\overline{\chi }}}_{\mathrm{o},h}^{\left(j\right)}$ in ${\tilde{\overline{\chi }}}_{\mathrm{o},h}$. Now define the residuals ${\mathit{d}}_h^{\left(j\right)}$ associated with subsystem j:

$$\begin{array}{cc}\hfill {\mathit{d}}_h^{\left(j\right)}& :={\tilde{\overline{\chi }}}_{\mathrm{o},h}^{\left(j\right)}-{\overline{\chi }}_{\mathrm{o}}^{\left(j\right)}\left(t_h\right)\hfill \end{array}$$

(51)

This way, the objective function $Q_0^{*}$ defined above can be reformulated as:

$$\begin{array}{cc}\hfill Q_0^{*}& =\sum _h^H\sum _i^J\sum _j^J{{\mathit{d}}_h^{\left(i\right)}}^{\mathrm{T}}{\mathbf{I}}_{{\mathbf{e}}^{\left(i\right)},•}\mathbf{W}{{\mathbf{I}}_{{\mathbf{e}}^{\left(j\right)},•}}^{\mathrm{T}}{\mathit{d}}_h^{\left(j\right)}\hfill \\ & =\sum _h^H\left(\sum _j^J{{\mathit{d}}_h^{\left(j\right)}}^{\mathrm{T}}{\mathbf{W}}^{\left(j\right)}{\mathit{d}}_h^{\left(j\right)}+\sum _{i,j:i\ne j}{{\mathit{d}}_h^{\left(i\right)}}^{\mathrm{T}}{\mathbf{I}}_{{\mathbf{e}}^{\left(i\right)},•}\mathbf{W}{{\mathbf{I}}_{{\mathbf{e}}^{\left(j\right)},•}}^{\mathrm{T}}{\mathit{d}}_h^{\left(j\right)}\right)\hfill \end{array}$$

(52)

and can subsequently approximated with $Q_1$:

$$\begin{array}{cc}\hfill Q_0^{*}& \approx Q_1:=\sum _j^J\sum _h^H{{\mathit{d}}_h^{\left(j\right)}}^{\mathrm{T}}{\mathbf{W}}^{\left(j\right)},{\mathit{d}}_h^{\left(j\right)}\hfill \end{array}$$

(53)

with ${\mathbf{W}}^{\left(j\right)}:={\mathbf{I}}_{{\mathbf{e}}^{\left(j\right)},•}\mathbf{W}{{\mathbf{I}}_{{\mathbf{e}}^{\left(j\right)},•}}^{\mathrm{T}}$. This first modification results in:

Problem $P_1$

$$\begin{array}{c}\hfill \widehat{\mathit{\theta }}=arg\underset{\mathit{\theta }}{min}{Q}_1:=\sum _j^J\sum _h^H{{\mathit{d}}_h^{\left(j\right)}}^{\mathrm{T}}{\mathbf{W}}^{\left(j\right)}{\mathit{d}}_h^{\left(j\right)},\end{array}$$

(54)

subject to (6)–(8) and (20). This error due to this approximation is referred to as Type A approximation error. This error is zero if all matrices ${\mathbf{I}}_{{\mathbf{e}}^{\left(i\right)},•}\mathbf{W}{{\mathbf{I}}_{{\mathbf{e}}^{\left(j\right)},•}}^{\mathrm{T}},i\ne j$, are zero matrices. This is true when the correlation between the estimation errors of any element in ${\tilde{\overline{\chi }}}_{\mathrm{o},h}^{\left(i\right)}$ and any element of ${\tilde{\overline{\chi }}}_{\mathrm{o},h}^{\left(j\right)}$ ($i\ne j$) is zero.

(b) 

Modification B: Separation of problem $P_1$ into J smaller problems $P_1^{\left(j\right)}$. An approximation to problem $P_1$ is now obtained by optimizing each of the J terms in (54) separately and simulating the values of ${\overline{\chi }}_{\mathrm{o}}^{\left(j\right)}\left(t_h\right)$ with (43) and (44). We refer to each problem as $P_1^{\left(j\right)}$:

Problem $P_1^{\left(j\right)}$

$$\begin{array}{c}\hfill {\widehat{\mathit{\theta }}}^{\left(j\right)}=arg\underset{\mathit{\theta }}{min}{Q}_1^{\left(j\right)}:=\sum _h^H{{\mathit{d}}_h^{\left(j\right)}}^{\mathrm{T}}{\mathbf{W}}^{\left(j\right)}{\mathit{d}}_h^{\left(j\right)}\end{array}$$

(55)

subject to

$$\begin{array}{ccc}\hfill {\dot{\mathit{x}}}^{\left(j\right)}\left(t\right)& =V{\mathit{f}}^{\left(j\right)}\left({\tilde{\overline{\chi }}}_{\mathrm{o},\mathrm{i}}^{\left(j\right)}\left(t\right),{\chi }_{\mathrm{o}}^{\left(j\right)}\left(t\right),{\mathit{x}}^{\left(j\right)}\left(t\right),{\mathit{\theta }}^{\left(j\right)}\right),\hfill & \hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{\mathit{x}}^{\left(j\right)}\left(0\right)=\mathbf{0}\end{array}$$

(56)

$$\begin{array}{cc}\hfill {\chi }_{\mathrm{o}}^{\left(j\right)}\left(t\right)& ={\mathbf{U}}^{\left(j\right)}{\mathit{x}}^{\left(j\right)}\left(t\right)\hfill \end{array}$$

(57)

$$\begin{array}{cc}\hfill {\mathit{d}}_h^{\left(j\right)}& ={\tilde{\overline{\chi }}}_{\mathrm{o},h}^{\left(j\right)}-{\overline{\chi }}_{\mathrm{o}}^{\left(j\right)}\left(t_h\right)={\tilde{\overline{\chi }}}_{\mathrm{o},h}^{\left(j\right)}-{\overline{\mathbf{U}}}^{\left(j\right)}{\mathit{x}}^{\left(j\right)}\left(t_h\right).\hfill \end{array}$$

(58)

The most important feature of problems $P_1^{\left(j\right)},j=1,\dots ,J$, is that each of them involves only a small number of parameters ${\mathit{\theta }}^{\left(j\right)}$. Please note that the approximation of problem $P_1$ by this set of problems $P_1^{\left(j\right)}$ is perfect in the special case where the right-hand sides (56) do not involve interpolated extents, that is, if the vectors ${\tilde{\overline{\chi }}}_{\mathrm{o},\mathrm{i}}^{\left(j\right)}\left(t\right)$ ($j=1,\dots ,J$) are empty. Another error, named Type B approximation error, is introduced when this is not the case.

Thanks to two modifications of problem $P_0$, one obtains J smaller problems $P_1^{\left(j\right)}$, each of which includes only a fraction of the original set of parameters. The price to pay for such a simplification is a deviation from maximum likelihood due to the introduction of approximation errors (type A & B). These errors are often marginal as demonstrated below.

As in previous work [5,6,7], the parameter estimates obtained by solving problem $P_1^{\left(j\right)}$ ($j=1,\dots ,J$) can serve as reliable initial guesses to initiate the solution to problem $P_0$. The problem $P_0$ when solved with these initial guesses is named problem $P_{1+0}$.

3.8. Implementation

All results can be reproduced with the open-source Efficient Model Identification (EMI) MATLAB package for efficient model identification [6] that includes all methods and simulations used in this study. This package is added in the Supplementary Information (Section A).

4. Results

We first explain the results obtained within an extensive simulation study to demonstrate the method. After that, results obtained with an experimental data set are used to demonstrate real-world applicability.

4.1. Simulation Study

The developed methods are demonstrated via a batch reaction system and investigating 5 different measurement scenarios (A–E). The reaction system is described first. Then, scenario A and its results are discussed in detail. The results for scenarios B to E are only summarized, with the details described in the Supplementary Information (Section B).

4.1.1. Reaction System

This reaction system has $R=5$ reactions involving $S=6$ species (A to F) with the following reaction scheme:

$$\begin{array}{cccc}\hfill R_1:& \hspace{1em}A\mathrm{+}\mathrm{B}\hfill & \hfill \to & C\hfill \\ \hfill R_2:& 2\mathrm{A}\hfill & \hfill \to & \mathrm{D}\hfill \\ \hfill R_3:& 2C\hfill & \hfill \to & \mathrm{B}+\mathrm{D}\hfill \\ \hfill R_4:& \hspace{1em}\mathrm{D}\hfill & \hfill \to & \mathrm{E}\hfill \\ \hfill R_5:& 2\mathrm{D}\hfill & \hfill \to & \mathrm{E}+\mathrm{F}\hfill \end{array}$$

with

$$\begin{array}{cc}\hfill \mathbf{N}& :=\left[\begin{array}{cccccc}-1& -1& +1& 0& 0& 0\\ -2& 0& 0& +1& 0& 0\\ 0& +1& -2& +1& 0& 0\\ 0& 0& 0& -1& +1& 0\\ 0& 0& 0& -2& +1& +1\end{array}\right].\hfill \end{array}$$

(59)

4.1.2. Dynamic Model in Terms of Numbers of Moles

The simulated kinetic rate expressions are:

$$\begin{array}{cc}\hfill \mathit{f}\left(\mathit{n}\left(t\right),\mathrm{\theta }\right)& :=\left[\begin{array}{c}{k}_1\left(\frac{n_1\left(t\right)}{V}\frac{n_2\left(t\right)}{V}-K_1\frac{n_3\left(t\right)}{V}\right)\\ k_2{\left(\frac{n_1\left(t\right)}{V}\right)}^2\\ k_3\frac{n_3\left(t\right)}{V}\\ k_4\frac{n_4\left(t\right)}{V}\\ k_5{\left(\frac{n_4\left(t\right)}{V}\right)}^2\end{array}\right]\hfill \end{array}$$

(60)

with $\mathrm{\theta }:={\left[\begin{array}{cccccc}{k}_1& k_2& k_3& k_4& k_5& K_1\end{array}\right]}^{\mathrm{T}}$. The ground truth values of the kinetic parameters are given in

Table 2. Scenario A—Parameter estimates. Model simulation and parameter estimation results: kinetic parameters and model fit. All parameter estimates are reported with their standard deviation based on the Laplacian approximation of the likelihood function (i.e., $exp\left(-Q_0\right)$, $exp\left(-Q_1^{\left(j\right)}\right)$).

	Name	Unit	True Value	${\mathit{P}}_0$	${\mathit{P}}_1^{\left(\mathit{j}\right)}$	${\mathit{P}}_{1+0}$
${\theta }_1$	$k_1$	$\mathrm{L}{\mathrm{mol}}^{-1}{\mathrm{h}}^{-1}$	$2.0$	$1.9490(\pm 0.0003)$	$2.0010(\pm 3.3625)$	$1.9857(\pm 0.0001)$
${\theta }_2$	$k_2$	$\mathrm{L}{\mathrm{mol}}^{-1}{\mathrm{h}}^{-1}$	$0.5$	$0.4929(\pm 0.0002)$	$0.5006(\pm 2.4090)$	$0.4771(\pm 0.0001)$
${\theta }_3$	$k_3$	${\mathrm{h}}^{-1}$	$1.0$	$1.0056(\pm 0.0827)$	$0.8955(\pm 0.0002)$	$0.9974(\pm 0.0805)$
${\theta }_4$	$k_4$	${\mathrm{h}}^{-1}$	$0.4$	$0.4342(\pm 0.4275)$	$0.4003(\pm 2.6400)$	$0.4606(\pm 0.1309)$
${\theta }_5$	$k_5$	$\mathrm{L}{\mathrm{mol}}^{-1}{\mathrm{h}}^{-1}$	$1.6$	$1.2285(\pm 2.5160)$	$1.6009(\pm 2.9587)$	$1.2111(\pm 0.8018)$
${\theta }_6$	$K_1$	$\mathrm{mol}{\mathrm{L}}^{-1}$	$1.4$	$1.3435(\pm 0.1016)$	$1.4001(\pm 2.9081)$	$1.3986(\pm 0.1135)$
$\sqrt{\frac{Q_0}{H·M}}$	WRMSR	−		$0.44731$	$0.47934$	$0.44771$

.

The initial conditions are ${\mathit{n}}_0:={\left[\begin{array}{cccccc}0.73& 0.42& 0& 0& 0& 0\end{array}\right]}^{\mathrm{T}}$ mol, and the reactor volume is $V=1L$. In all investigated scenarios, measurements are taken at intervals of 5 min during 10 h ($H=121$, $t_h=0,1/12,2/12,\dots ,10$).

4.1.3. Scenario A

The first scenario considers the case where the concentrations of B and C and the sum of the concentrations of E and F are measured, that is,

$$\begin{array}{c}\hfill \mathbf{M}:=\left[\begin{array}{cccccc}0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 0& 1& 1\end{array}\right].\end{array}$$

(61)

The measurement error-covariance matrix is ${\Sigma }_{\mathbf{ϵ}}=\mathrm{diag}\left({\left[\begin{array}{ccc}1& 1& 2\end{array}\right]}^{\mathrm{T}}\right)$ ${10}^{-4}{\mathrm{mol}}^2{\mathrm{L}}^{-2}$. Figure 2 shows the simulated noise-free and noisy measurements.

Extent labeling

To specify the measurement model in terms of extents (7), one needs to specify ${\mathit{y}}_0$ and $\mathbf{G}$:

$$\begin{array}{cc}\hfill {\mathit{y}}_0& =\mathbf{M}{\mathit{n}}_0={\left[\begin{array}{ccc}0.42& 0& 0\end{array}\right]}^{\mathrm{T}}\hfill \end{array}$$

(62)

$$\begin{array}{cc}\hfill \mathbf{G}& =\mathbf{M}{\mathbf{N}}^{\mathrm{T}}=\left[\begin{array}{ccccc}-1& 0& +1& 0& 0\\ +1& 0& -2& 0& 0\\ 0& 0& 0& +1& +2\end{array}\right].\hfill \end{array}$$

(63)

The reduced row echelon form of $\mathbf{G}$ is

$$\begin{array}{cc}\hfill \mathbf{B}& =\left[\begin{array}{ccccc}+1& 0& 0& 0& 0\\ 0& 0& +1& 0& 0\\ 0& 0& 0& +1& +2\end{array}\right],\hfill \end{array}$$

(64)

which indicates that there are:

• $R_{\mathrm{n}}=1$ non-sensed extent, ${\mathit{x}}_{\mathbf{n}}\left(t\right)=x_2\left(t\right)$,
• $R_{\mathrm{o}}=2$ observable extents, ${\mathit{x}}_{\mathbf{o}}\left(t\right)={\left[\begin{array}{cc}{x}_1\left(t\right)& x_3\left(t\right)\end{array}\right]}^{\mathrm{T}}$,
• $R_{\mathrm{a}}=2$ ambiguous extents, ${\mathit{x}}_{\mathbf{a}}\left(t\right)={\left[\begin{array}{cc}{x}_4\left(t\right)& x_5\left(t\right)\end{array}\right]}^{\mathrm{T}}$,

which gives:

$$\begin{array}{cc}\hfill {\mathbf{G}}_{•,\mathbf{a}}& =\left[\begin{array}{cc}0& 0\\ 0& 0\\ +1& +2\end{array}\right].\hfill \end{array}$$

(65)

Observable extents and extent directions

The factorization of ${\mathbf{G}}_{•,\mathbf{a}}$ gives:

$$\begin{array}{c}\hfill {\mathcal{G}}_{\mathrm{o}}=\left[\begin{array}{c}0\\ 0\\ +1\end{array}\right],{{\mathbf{V}}_{\mathrm{o}}}^{\mathrm{T}}=\left[\begin{array}{cc}+1& +2\end{array}\right],\end{array}$$

(66)

so that

$$\begin{array}{cc}\hfill {\chi }_{\mathrm{o}}\left(t\right)& ={{\mathbf{V}}_{\mathrm{o}}}^{\mathrm{T}}{\mathbf{I}}_{\mathbf{a},•}\mathit{x}\left(t\right)=\left[\begin{array}{ccccc}0& 0& 0& 1& 2\end{array}\right]\mathit{x}\left(t\right)=x_4\left(t\right)+2x_5\left(t\right)\hfill \end{array}$$

(67)

$$\begin{array}{cc}\hfill {\overline{\chi }}_{\mathrm{o}}\left(t\right)& =\left[\begin{array}{c}{x}_1\left(t\right)\\ x_3\left(t\right)\\ {\chi }_{\mathrm{o}}\left(t\right)\end{array}\right]=\left[\begin{array}{c}{x}_1\left(t\right)\\ x_3\left(t\right)\\ x_4\left(t\right)+2x_5\left(t\right)\end{array}\right].\hfill \end{array}$$

(68)

The values of the observable extents and extent directions ${\tilde{\overline{\chi }}}_{\mathrm{o},h}:={\left[\begin{array}{c}{\tilde{x}}_{1,h}{\tilde{x}}_{3,h}{\tilde{\chi }}_{\mathrm{o},h}\end{array}\right]}^{\mathrm{T}}$ can be computed from (23) with

$$\begin{array}{cc}\hfill \mathbf{P}& =\left[\begin{array}{ccc}-2& -1& 0\\ -1& -1& 0\\ 0& 0& +1\end{array}\right].\hfill \end{array}$$

(69)

The expected estimation error variance-covariance matrix ${\Sigma }_{\overline{\chi }}$ is:

$$\begin{array}{cc}\hfill {\Sigma }_{\overline{\chi }}& =\left[\begin{array}{ccc}+5& +3& 0\\ +3& +2& 0\\ 0& 0& +2\end{array}\right]{10}^{-4}.\hfill \end{array}$$

(70)

This demonstrates that the estimation errors in the first and second observable extents are uncorrelated with the estimation error in the computed observable extent direction. Please note that the correlation between the estimation errors of the first and second observable extents is fairly high ($\frac{3}{\sqrt{5·2}}=0.95$). Figure 3 shows the computed values of the observable extents and extent direction.

Expression (27) is then fully specified with the following vectors and matrices:

$$\begin{array}{cc}\hfill {\mathit{n}}_0& =\left[\begin{array}{c}1\\ 1\\ 0\\ 0\\ 0\\ 0\end{array}\right]\hfill \end{array}$$

(71)

$$\begin{array}{cc}\hfill {\overline{\mathcal{N}}}_{\mathrm{o}}& :=\left[\begin{array}{cccccc}-1& -1& +1& 0& 0& 0\\ 0& +1& -2& +1& 0& 0\\ 0& 0& 0& -1& \frac{3}{5}& \frac{2}{5}\end{array}\right].\hfill \end{array}$$

(72)

$$\begin{array}{cc}\hfill {\overline{\mathcal{N}}}_{\mathrm{u}}& :=\left[\begin{array}{cccccc}0& 0& 0& 0& 0& 0\\ -2& 0& 0& +1& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& +\frac{2}{5}& -\frac{2}{5}\\ 0& 0& 0& 0& -\frac{1}{5}& +\frac{1}{5}\end{array}\right]\hfill \end{array}$$

(73)

Since the 2nd and 3rd column of ${\overline{\mathcal{N}}}_{\mathrm{u}}$ are null vectors, it follows that values of ${\tilde{n}}_2\left(t\right)$ and ${\tilde{n}}_3\left(t\right)$ can be computed based on the interpolated estimates ${\tilde{\overline{\chi }}}_{\mathrm{o},\mathrm{i}}\left(t\right)$ without the need to simulate $\mathit{x}\left(t\right)$. A further exploration of this idea remains out of the scope of this paper, however.

System partitioning

For simplicity of notation, we omit the time dependence in what follows. For example, the interpolants ${\tilde{\overline{\chi }}}_{\mathrm{o},\mathrm{i}}\left(t\right)$ are given as ${\tilde{\overline{\chi }}}_{\mathrm{o},\mathrm{i}}:={\left[\begin{array}{ccc}{\tilde{x}}_1& {\tilde{x}}_3& {\tilde{\chi }}_{\mathrm{o}}\end{array}\right]}^{\mathrm{T}}$. Following Steps 1(a)–(b) in Section 3.6.1, the augmented equation system becomes:

$$\begin{array}{cccc}\hfill \dot{\mathit{x}}& =\left[\begin{array}{c}\frac{k_1}{V^2}\left[\left(n_{0,1}-{\tilde{x}}_1-2x_2\right)\left(n_{0,2}-{\tilde{x}}_1+{\tilde{x}}_3\right)-K_1\left(n_{0,3}+{\tilde{x}}_1-2{\tilde{x}}_3\right)\right]\\ \frac{k_2}{V^2}{\left(n_{0,1}-{\tilde{x}}_1-2x_2\right)}^2\\ \frac{k_3}{V}\left(n_{0,3}+{\tilde{x}}_1-2{\tilde{x}}_3\right)\\ \frac{k_4}{V}\left(n_{0,4}+{\tilde{x}}_3-{\tilde{\chi }}_{\mathrm{o}}+x_2\right)\\ \frac{k_5}{V^2}{\left(n_{0,4}+{\tilde{x}}_3-{\tilde{\chi }}_{\mathrm{o}}+x_2\right)}^2\end{array}\right],\hfill & \hfill \hspace{1em}\hspace{1em}\mathit{x}\left(0\right)& =\mathbf{0}\hfill \end{array}$$

(74)

$$\begin{array}{cc}\hfill {\chi }_{\mathrm{o}}& ={{\mathbf{V}}_{\mathrm{o}}}^{\mathrm{T}}{\mathbf{I}}_{\mathbf{a},•}\mathit{x}=\left[\begin{array}{ccccc}0& 0& 0& 1& 2\end{array}\right]\mathit{x}=x_4+2x_5.\hfill \end{array}$$

(75)

Figure 4 shows the graph corresponding to the above equation system. The vertices corresponding to the observable extents and extent direction are shaded, while the other vertices are white. The simulation arcs are shown with full-line arrows, while the observation arcs are shown as dashed-line arrows. To identify possible subsystems, one removes all the observation arcs, which results in two subgraphs. The first subgraph includes the parameters $k_1$, $k_2$, $k_4$, $k_5$, and $K_1$, which affect the observable quantities $x_1$ and ${\chi }_{\mathrm{o}}$ via a network that also involves the unobservable extents $x_2$, $x_4$, and $x_5$. The second subgraph is much smaller and includes the parameter $k_3$ that influences the observable extent $x_3$.

Accordingly, the $J=2$ subsystems are:

$$\begin{array}{cccc}\hfill {\dot{\mathit{x}}}^{\left(1\right)}& =\left[\begin{array}{c}{\dot{x}}_1\\ \\ {\dot{x}}_2\\ {\dot{x}}_4\\ {\dot{x}}_5\end{array}\right]=\left[\begin{array}{c}\frac{k_1}{V^2}[\left(n_{0,1}-x_1-2x_2\right)\left(n_{0,2}-x_1+{\tilde{x}}_3\right)\\ -K_1\left(n_{0,3}+x_1-2{\tilde{x}}_3\right)]\\ \frac{k_2}{V^2}{\left(n_{0,1}-x_1-2x_2\right)}^2\\ \frac{k_4}{V}\left(n_{0,4}+{\tilde{x}}_3-{\chi }_{\mathrm{o}}+x_2\right)\\ \frac{k_5}{V^2}{\left(n_{0,4}+{\tilde{x}}_3-{\chi }_{\mathrm{o}}+x_2\right)}^2\end{array}\right],\hfill & \hfill {\mathit{x}}^{\left(1\right)}\left(0\right)& =\left[\begin{array}{c}{x}_{0,1}\\ \\ x_{0,2}\\ x_{0,4}\\ x_{0,5}\end{array}\right]=\mathbf{0}\hfill \end{array}$$

(76)

$$\begin{array}{cc}\hfill {\chi }_{\mathrm{o}}^{\left(1\right)}& ={\chi }_{\mathrm{o}}=x_4+2x_5,\hfill \end{array}$$

(77)

and

$$\begin{array}{cccc}\hfill {\dot{x}}^{\left(2\right)}& ={\dot{x}}_3=\frac{k_3}{V}\left(n_{0,3}+{\tilde{x}}_1-2x_3\right),\hfill & \hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{x}^{\left(2\right)}\left(0\right)& =x_{0,3}=0.\hfill \end{array}$$

(78)

The information flows in these subsystems are shown in the inset of Figure 4.

Incremental parameter estimation

The resulting partitioning means that estimates of the first set of parameters ${\mathit{\theta }}^{\left(1\right)}:=\left[\begin{array}{c}{k}_1,k_2,k_4,k_5,K_1\end{array}\right]$ can be obtained by minimizing a weighted least-squares deviation between the predicted $x_1\left(t_h\right)$ and ${\chi }_{\mathrm{o}}\left(t_h\right)$ and their measured counterparts ${\tilde{x}}_{1,h}$ and ${\tilde{\chi }}_{\mathrm{o},h}$, that is, by solving problem $P_1^{\left(1\right)}$:

$$\begin{array}{c}\hfill {\widehat{\mathit{\theta }}}^{\left(1\right)}=arg\underset{{\mathit{\theta }}^{\left(1\right)}}{min}\sum _h^H{{\mathit{d}}_h^{\left(1\right)}}^{\mathrm{T}}{\mathbf{W}}^{\left(1\right)}{\mathit{d}}_h^{\left(1\right)}\end{array}$$

(79)

subject to

$$\begin{array}{cccc}\hfill \left[\begin{array}{c}{\dot{x}}_1\left(t\right)\\ \\ {\dot{x}}_2\left(t\right)\\ {\dot{x}}_4\left(t\right)\\ {\dot{x}}_5\left(t\right)\end{array}\right]& =\left[\begin{array}{c}\frac{k_1}{V^2}[\left(n_{0,1}-x_1\left(t\right)-2x_2\left(t\right)\right)\left(n_{0,2}-x_1\left(t\right)+{\tilde{x}}_3\left(t\right)\right)\\ -K_1\left(n_{0,3}+x_1\left(t\right)-2{\tilde{x}}_3\left(t\right)\right)]\\ \frac{k_2}{V^2}{\left(n_{0,1}-x_1\left(t\right)-2x_2\left(t\right)\right)}^2\\ \frac{k_4}{V}\left(n_{0,4}+{\tilde{x}}_3\left(t\right)-{\chi }_{\mathrm{o}}\left(t\right)+x_2\left(t\right)\right)\\ \frac{k_5}{V^2}{\left(n_{0,4}+{\tilde{x}}_3\left(t\right)-{\chi }_{\mathrm{o}}\left(t\right)+x_2\left(t\right)\right)}^2\end{array}\right],\hfill & \hfill \left[\begin{array}{c}{x}_1\left(0\right)\\ \\ x_2\left(0\right)\\ x_4\left(0\right)\\ x_5\left(0\right)\end{array}\right]& =\mathbf{0}\hfill \end{array}$$

(80)

$$\begin{array}{cc}\hfill {\chi }_{\mathrm{o}}\left(t\right)& =x_4\left(t\right)+2x_5\left(t\right)\hfill \end{array}$$

(81)

$$\begin{array}{cc}\hfill {\mathit{d}}_h^{\left(1\right)}& =\left[\begin{array}{c}{\tilde{x}}_{1,h}\\ {\tilde{\chi }}_{\mathrm{o},h}\end{array}\right]-\left[\begin{array}{c}{x}_1\left(t_h\right)\\ {\chi }_{\mathrm{o}}\left(t_h\right)\end{array}\right],\hfill \end{array}$$

(82)

where

$$\begin{array}{cc}\hfill {\mathbf{W}}^{\left(1\right)}& :={\mathbf{L}}_{\mathrm{e}}^{\left(1\right)}{\Sigma }_{\overline{\chi }}^{-1}{{\mathbf{L}}_{\mathrm{e}}^{\left(1\right)}}^{\mathrm{T}}=\left[\begin{array}{ccc}+1& 0& 0\\ 0& 0& +1\end{array}\right]\left(\left[\begin{array}{ccc}2& -3& 0\\ -3& 5& 0\\ 0& 0& 1/2\end{array}\right]{10}^4\right){\left[\begin{array}{ccc}1& 0& 0\\ 0& 0& 1\end{array}\right]}^{\mathrm{T}}\hfill \end{array}$$

(83)

$$\begin{array}{c}=\left[\begin{array}{cc}2& 0\\ 0& 1/2\end{array}\right]{10}^4.\hfill \end{array}$$

(84)

The second set of parameters consists of $k_3$ only, ${\theta }^{\left(2\right)}:=k_3$. Its value can be obtained by minimizing the least-squares deviation between $x_3\left(t_h\right)$ and its measured counterparts ${\tilde{x}}_{3,h}$, that is, by solving problem $P_1^{\left(2\right)}$:

$$\begin{array}{c}\hfill {\widehat{\theta }}^{\left(2\right)}=arg\underset{{\theta }^{\left(2\right)}}{min}\sum _h^H{\left(d_h^{\left(2\right)}\right)}^2\end{array}$$

(85)

subject to

$$\begin{array}{cccc}\hfill {\dot{x}}_3\left(t\right)& =\frac{k_3}{V}\left(n_{0,3}+{\tilde{x}}_1\left(t\right)-2x_3\left(t\right)\right),\hfill & \hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{x}_3\left(0\right)& =0\hfill \end{array}$$

(86)

$$\begin{array}{cc}\hfill d_h^{\left(2\right)}& ={\tilde{x}}_{3,h}-x_3\left(t_h\right).\hfill \end{array}$$

(87)

Figure 5 shows the measured concentrations and the simulated profiles obtained with (i) the true parameters; (ii) simultaneous parameter estimation ($P_0$); and (iii) incremental parameter estimation ($P_1^{\left(1\right)}$ and $P_1^{\left(2\right)}$). All simulated profiles describe the experiment well and are practically indistinguishable. The parameter values obtained by solving $P_0$, $P_1^{\left(1\right)}$, $P_1^{\left(2\right)}$, and $P_{1+0}$ are shown in Table 2. All parameter estimates are in close agreement to each other. It is worth noting that the parameter $k_2$ can be estimated even if the corresponding extent $x_2$ is labeled non-sensed, i.e., our results suggest that this parameter is both structurally and practically identifiable.

4.1.4. Scenario B

Scenario B assumes concentration measurements for the species B and C only:

$$\begin{array}{c}\hfill \mathbf{M}:=\left[\begin{array}{cccccc}0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\end{array}\right].\end{array}$$

(88)

The measurement error-covariance matrix is ${\Sigma }_{\mathbf{ϵ}}=\mathrm{diag}\left({\left[\begin{array}{cc}1& 1\end{array}\right]}^{\mathrm{T}}\right)$ ${10}^{-4}{\mathrm{mol}}^2{\mathrm{L}}^{-2}$. The reduced row echelon form of $\mathbf{G}$ is:

$$\begin{array}{cc}\hfill \mathbf{B}& =\left[\begin{array}{ccccc}+1& 0& 0& 0& 0\\ 0& 0& +1& 0& 0\end{array}\right],\hfill \end{array}$$

(89)

which leads to the following extent labeling:

• $R_{\mathrm{n}}=3$ non-sensed extents, ${\mathit{x}}_{\mathbf{n}}\left(t\right)={\left[\begin{array}{ccc}{x}_2\left(t\right)& x_4\left(t\right)& x_5\left(t\right)\end{array}\right]}^{\mathrm{T}}$,
• $R_{\mathrm{o}}=2$ observable extents, ${\mathit{x}}_{\mathbf{o}}\left(t\right)={\left[\begin{array}{cc}{x}_1\left(t\right)& x_3\left(t\right)\end{array}\right]}^{\mathrm{T}}$,
• $R_{\mathrm{a}}=0$ ambiguous extents.

Please note that the extent-based method proposed in [8] applies in this case since $R_{\mathrm{a}}=0$.

Figure 6 shows the corresponding graph $\mathcal{F}$, which can be partitioned into $J=2$ subsystems and a leftover part. The first subgraph includes the parameters $k_1$, $k_2$, and $K_1$ since they all have some effect on the observable extent $x_1$ via a network that also involves the unobservable extent $x_2$. The second subgraph includes the parameter $k_3$ that influences the observable extent $x_3$. The information flows in these subsystems are shown in the inset of Figure 6. The leftover part includes the parameters $k_4$ and $k_5$ and the extents $x_4$ and $x_5$. These parameters and variables are not part of any of the identified subsystems because there are no directed paths composed of simulation arcs from the corresponding vertices to any of the observable extents. Hence, it follows that these parameters are unidentifiable. This is consistent with the method proposed in [8] (not shown). Clearly, the graph-based partitioning procedure shows potential for parameter identifiability analysis, which is discussed as an opportunity below.

4.1.5. Scenario C

Scenario C assumes concentration measurements for the species B, C, E, and F:

$$\begin{array}{c}\hfill \mathbf{M}:=\left[\begin{array}{cccccc}0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right].\end{array}$$

(90)

The measurement error-covariance matrix is ${\Sigma }_{\mathbf{ϵ}}=\mathrm{diag}\left({\left[\begin{array}{cccc}1& 1& 1& 1\end{array}\right]}^{\mathrm{T}}\right)$ ${10}^{-4}{\mathrm{mol}}^2{\mathrm{L}}^{-2}$. The reduced row echelon form of $\mathbf{G}$ is:

$$\begin{array}{cc}\hfill \mathbf{B}& =\left[\begin{array}{ccccc}+1& 0& 0& 0& 0\\ 0& 0& +1& 0& 0\\ 0& 0& 0& +1& 0\\ 0& 0& 0& 0& +1\end{array}\right],\hfill \end{array}$$

(91)

which leads to the following extent labeling:

• $R_{\mathrm{n}}=1$ non-sensed extents, ${\mathit{x}}_{\mathbf{n}}\left(t\right)=x_2\left(t\right)$.
• $R_{\mathrm{o}}=4$ observable extents, ${\mathit{x}}_{\mathbf{o}}\left(t\right)={\left[\begin{array}{cccc}{x}_1\left(t\right)& x_3\left(t\right)& x_4\left(t\right)& x_6\left(t\right)\end{array}\right]}^{\mathrm{T}}$.
• $R_{\mathrm{a}}=0$ ambiguous extents.

Figure 7 shows the corresponding graph $\mathcal{F}$, which can be partitioned into 2 subsystems. The first subgraph includes the parameters $k_1$, $k_2$, $k_4$, $k_5$, and $K_1$ since they all have some effect on the observable extents $x_1$, $x_4$, and $x_5$ via a network that also involves the unobservable extents $x_2$. The second subgraph includes the parameter $k_3$ that influences the observable extent $x_3$. The information flows in these subsystems are shown in the inset of Figure 7. Although the assignment of the parameters is the same as in scenario A, the parameter estimation problems are different since the fit of subsystem 1 is determined with respect to ${\tilde{x}}_1$, ${\tilde{x}}_4$, and ${\tilde{x}}_5$ in scenario C, compared to only ${\tilde{x}}_1$ and ${\tilde{\chi }}_{\mathrm{o}}$ in scenario A.

4.1.6. Scenario D

Scenario D assumes concentration measurements for the species A, C, and E:

$$\begin{array}{c}\hfill \mathbf{M}:=\left[\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 0& 1& 0\end{array}\right].\end{array}$$

(92)

The measurement error-covariance matrix is ${\Sigma }_{\mathbf{ϵ}}=\mathrm{diag}\left({\left[\begin{array}{ccc}1& 1& 1\end{array}\right]}^{\mathrm{T}}\right)$ ${10}^{-4}{\mathrm{mol}}^2{\mathrm{L}}^{-2}$. The reduced row echelon form of $\mathbf{G}$ is:

$$\begin{array}{cc}\hfill \mathbf{B}& =\left[\begin{array}{ccccc}+1& 0& -2& 0& 0\\ 0& +1& +1& 0& 0\\ 0& 0& 0& +1& +1\\ 0& 0& 0& 0& 0\end{array}\right],\hfill \end{array}$$

(93)

which leads to the following extent labeling:

• $R_{\mathrm{n}}=0$ non-sensed extents,
• $R_{\mathrm{o}}=0$ observable extents,
• $R_{\mathrm{a}}=5$ ambiguous extents, ${\mathit{x}}_{\mathbf{a}}\left(t\right)={\left[\begin{array}{ccccc}{x}_1\left(t\right)& x_2\left(t\right)& x_3\left(t\right)& x_4\left(t\right)& x_6\left(t\right)\end{array}\right]}^{\mathrm{T}}$.

From $A=\mathrm{rank}\left(\mathbf{B}\right)=3$ and $R_{\mathrm{o}}=0$, one sees that there are ${\rho }_{\mathrm{o}}=A-R_{\mathrm{o}}=3$ observable extent directions among the 5 ambiguous extents. The first two are linear combinations of $x_1$, $x_2$, and $x_3$, while the third one is a linear combination of $x_4$ and $x_5$. This particular appearance of extents in subsets of the observable directions stems from the subspace clustering property of the reduced row echelon form and will be discussed in some more detail below.

Figure 8 shows the corresponding graph $\mathcal{F}$, which can be partitioned into 2 subsystems. The first subgraph includes the parameters $k_1$, $k_2$, $k_3$, and $K_1$ since they all have some effect on the observable extent directions ${\chi }_{\mathrm{o},1}$ and ${\chi }_{\mathrm{o},2}$ via a network that also involves the unobservable extents $x_1$, $x_2$, and $x_3$. The second subgraph includes the parameters $k_4$ and $k_5$ that influence the observable extent direction ${\chi }_{\mathrm{o},3}$ via a network that also involves the unobservable extents $x_4$ and $x_5$. Hence, this leads to two parameter estimation problems involving 4 and 2 parameters, respectively. The information flows to simulate each of the corresponding subsystems are shown in the inset of Figure 8.

4.1.7. Scenario E

Scenario E assumes concentration measurements of all species:

$$\begin{array}{c}\hfill \mathbf{M}:={\mathbf{I}}_6.\end{array}$$

(94)

The measurement error-covariance matrix is ${\Sigma }_{\mathbf{ϵ}}=\mathrm{diag}\left({\left[\begin{array}{cccccc}1& 1& 1& 1& 1& 1\end{array}\right]}^{\mathrm{T}}\right)$ ${10}^{-4}{\mathrm{mol}}^2{\mathrm{L}}^{-2}$. The reduced row echelon form of $\mathbf{G}$ is:

$$\begin{array}{cc}\hfill \mathbf{B}& =\left[\begin{array}{ccccc}+1& 0& 0& 0& 0\\ 0& +1& 0& 0& 0\\ 0& 0& +1& 0& 0\\ 0& 0& 0& +1& 0\\ 0& 0& 0& 0& +1\\ 0& 0& 0& 0& 0\end{array}\right],\hfill \end{array}$$

(95)

which leads to the following extent labeling:

• $R_{\mathrm{n}}=0$ non-sensed extents,
• $R_{\mathrm{o}}=5$ observable extents, ${\mathit{x}}_{\mathbf{o}}\left(t\right)={\left[\begin{array}{ccccc}{x}_1\left(t\right)& x_2\left(t\right)& x_3\left(t\right)& x_4\left(t\right)& x_6\left(t\right)\end{array}\right]}^{\mathrm{T}}$,
• $R_{\mathrm{a}}=0$ ambiguous extents.

This means that the original framework for extent computation, which assumes $A=R$ (Section 3.2.3), applies.

Figure 9 shows the graph $\mathcal{F}$, which can be partitioned into 5 subsystems that include the parameters as follows: (i) $k_1$ and $K_1$; (ii) $k_2$; (iii) $k_3$; (iv) $k_4$; and (v) $k_5$. Hence, one obtains 5 entirely decoupled parameter estimation problems. This set of optimization problems is the same as those obtained with the original extent-based model identification method as is also clear from the inset of Figure 9.

4.2. Experimental Study

To explore the applicability of the novel incremental parameter estimation method with realistic laboratory data, we execute parameter estimation for the $\alpha$-pinene batch experiment first reported in [11] and later studied in [12,13,14,15,16]. The reaction network consists of 5 species. These species are labeled $\mathrm{A}$, $\mathrm{B}$, $\mathrm{C}$, $\mathrm{D}$, and $\mathrm{E}$ and correspond to $\alpha$-pinene, dipentene, allo-ocimene, pyronene, and dimer. The reaction network is:

$$\begin{array}{cccc}\hfill R_1:& \mathrm{A}\hfill & \hfill \to & \mathrm{B}\hfill \\ \hfill R_2:& \mathrm{A}\hfill & \hfill \to & \mathrm{C}\hfill \\ \hfill R_3:& \mathrm{C}\hfill & \hfill \to & \mathrm{D}\hfill \\ \hfill R_4:& \mathrm{C}\hfill & \hfill \to & \mathrm{E}\hfill \\ \hfill R_5:& \mathrm{E}\hfill & \hfill \to & \mathrm{C}\hfill \end{array}$$

As in [12], the model is defined with the following stoichiometric matrix and process rates:

$$\begin{array}{cc}\hfill \mathbf{N}& :=\left[\begin{array}{ccccc}-1& +1& 0& 0& 0\\ -1& 0& +1& 0& 0\\ 0& 0& -1& +1& 0\\ 0& 0& -1& 0& +1\\ 0& 0& +1& 0& -1\end{array}\right]\hfill \end{array}$$

(96)

$$\begin{array}{cc}\hfill \mathit{r}& :=\left[\begin{array}{c}{k}_1\frac{n_1\left(t\right)}{V}\\ k_2\frac{n_1\left(t\right)}{V}\\ k_3\frac{n_3\left(t\right)}{V}\\ k_4\frac{n_3\left(t\right)}{V}\\ k_5\frac{n_5\left(t\right)}{V}\end{array}\right].\hfill \end{array}$$

(97)

At the start of the batch experiment, only $\mathrm{A}$ is present. As explained in [12], the concentrations of all species but $\mathrm{D}$ are measured. The concentrations of $\mathrm{D}$ are obtained by assuming that 3% of the transformed $\mathrm{A}$ is present as $\mathrm{D}$ at all times. Finally, all concentrations are normalized to 1 (100%). Figure 10 shows the experimental data reported in [12]. As the initial numbers of moles and the volume are unknown, we express both the concentrations and the extents as dimensionless fractions of the numbers of moles relative to the initial numbers of moles of $\mathrm{A}$. Hence, we can write ${\mathit{n}}_0:={\left[\begin{array}{ccccc}1& 0& 0& 0& 0\end{array}\right]}^{\mathrm{T}}$.

Extent labeling

For extent labeling, we assume that all species are measured with independent zero-mean Gaussian measurement errors, similar to the least-squares approach of [12]:

$$\begin{array}{c}\hfill \mathbf{M}:={\mathbf{I}}_5.\end{array}$$

(98)

It follows that

$$\begin{array}{cc}\hfill \mathbf{G}& =\mathbf{M}{\mathbf{N}}^{\mathrm{T}}=\left[\begin{array}{ccccc}-1& -1& 0& 0& 0\\ 1& 0& 0& 0& 0\\ 0& 1& -1& -1& 1\\ 0& 0& 1& 0& 0\\ 0& 0& 0& 1& -1\end{array}\right].\hfill \end{array}$$

(99)

The reduced row echelon form of $\mathbf{G}$ is

$$\begin{array}{cc}\hfill \mathbf{B}& =\left[\begin{array}{ccccc}+1& 0& 0& 0& 0\\ 0& +1& 0& 0& 0\\ 0& 0& +1& 0& 0\\ 0& 0& 0& +1& -1\end{array}\right],\hfill \end{array}$$

(100)

which indicates that there are:

• $R_{\mathrm{n}}=0$ non-sensed extents,
• $R_{\mathrm{o}}=3$ observable extents, ${\mathit{x}}_{\mathbf{o}}\left(t\right)={\left[\begin{array}{ccc}{x}_1\left(t\right)& x_2\left(t\right)& x_3\left(t\right)\end{array}\right]}^{\mathrm{T}}$,
• $R_{\mathrm{a}}=2$ ambiguous extents, ${\mathit{x}}_{\mathbf{a}}\left(t\right)={\left[\begin{array}{cc}{x}_4\left(t\right)& x_5\left(t\right)\end{array}\right]}^{\mathrm{T}}$,

which delivers:

$$\begin{array}{cc}\hfill {\mathbf{G}}_{•,\mathbf{a}}& =\left[\begin{array}{cc}0& 0\\ 0& 0\\ -1& +1\\ 0& 0\\ +1& -1\end{array}\right].\hfill \end{array}$$

(101)

Observable extents and extent directions

The factorization of ${\mathbf{G}}_{•,\mathbf{a}}$ gives:

$$\begin{array}{c}\hfill {\mathcal{G}}_{\mathrm{o}}=\left[\begin{array}{c}0\\ 0\\ 0\\ +1\end{array}\right],{{\mathbf{V}}_{\mathrm{o}}}^{\mathrm{T}}=\left[\begin{array}{cc}+1& -1\end{array}\right],\end{array}$$

(102)

so that

$$\begin{array}{cc}\hfill {\chi }_{\mathrm{o}}\left(t\right)& ={{\mathbf{V}}_{\mathrm{o}}}^{\mathrm{T}}{\mathbf{I}}_{\mathbf{a},•}\mathit{x}\left(t\right)=\left[\begin{array}{ccccc}0& 0& 0& 1& -1\end{array}\right]\mathit{x}\left(t\right)=x_4\left(t\right)-x_5\left(t\right)\hfill \end{array}$$

(103)

$$\begin{array}{cc}\hfill {\overline{\chi }}_{\mathrm{o}}\left(t\right)& =\left[\begin{array}{c}{x}_1\left(t\right)\\ x_2\left(t\right)\\ x_3\left(t\right)\\ {\chi }_{\mathrm{o}}\left(t\right)\end{array}\right]=\left[\begin{array}{c}{x}_1\left(t\right)\\ x_2\left(t\right)\\ x_3\left(t\right)\\ x_4\left(t\right)-x_5\left(t\right)\end{array}\right].\hfill \end{array}$$

(104)

Figure 11 shows the computed values of the observable extents and extent direction assuming these definitions. The corresponding error variance-covariance matrix is:

$$\begin{array}{cc}\hfill {\Sigma }_{\overline{\chi }}& =\left[\begin{array}{cccc}+0.8& -0.6& -0.2& -0.2\\ -0.6& +1.2& +0.4& +0.4\\ -0.2& +0.4& +0.8& -0.2\\ -0.2& +0.4& -0.2& +0.8\end{array}\right]\sigma \hfill \end{array}$$

(105)

with $\sigma$ the (unknown) measurement error variance.

System partitioning

As above, we omit the time dependence in what follows. Following Steps 1(a)–(b) in Section 3.6.1, the augmented equation system becomes:

$$\begin{array}{cccc}\hfill \dot{\mathit{x}}& =\left[\begin{array}{c}{k}_1\left(n_{0,1}-{\tilde{x}}_1-2{\tilde{x}}_2\right)\\ k_2\left(n_{0,1}-{\tilde{x}}_1-2{\tilde{x}}_2\right)\\ k_3\left(n_{0,3}+{\tilde{x}}_2-{\tilde{x}}_3-{\tilde{\chi }}_{\mathrm{o}}\right)\\ k_4\left(n_{0,4}+{\tilde{x}}_3\right)\\ k_5\left(n_{0,5}+{\tilde{\chi }}_{\mathrm{o}}\right)\end{array}\right],\hfill & \hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\mathit{x}\left(0\right)& =\mathbf{0},\hfill \end{array}$$

(106)

$$\begin{array}{cc}\hfill {\chi }_{\mathrm{o}}& =x_4-x_5.\hfill \end{array}$$

(107)

Figure 12 shows the graph corresponding to the above equation system. Removing all the observation arcs generates 4 subgraphs. The first three subgraphs represent the first three reactions and include: $k_1$ and $x_1$ $(j=1)$; $k_2$ and $x_2$ $(j=2)$; and $k_3$ and $x_3$ $(j=3)$. The fourth graph includes the parameters $k_4$ and $k_5$ via a network that also involves the observable extent direction ${\chi }_{\mathrm{o}}$ and the unobservable extents $x_4$ and $x_5$ $(j=4)$.

Accordingly, the $J=4$ subsystems are:

$$\begin{array}{cccc}\hfill {\dot{\mathit{x}}}^{\left(1\right)}& ={\dot{x}}_1=k_1\left(n_{0,1}-x_1-2{\tilde{x}}_2\right),\hfill & \hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{\mathit{x}}^{\left(1\right)}\left(0\right)& =x_{0,1}=0,\hfill \end{array}$$

(108)

$$\begin{array}{cccc}\hfill {\dot{\mathit{x}}}^{\left(2\right)}& ={\dot{x}}_2=k_2\left(n_{0,1}-{\tilde{x}}_1-2x_2\right),\hfill & \hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{\mathit{x}}^{\left(2\right)}\left(0\right)& =x_{0,2}=0,\hfill \end{array}$$

(109)

$$\begin{array}{cccc}\hfill {\dot{\mathit{x}}}^{\left(3\right)}& ={\dot{x}}_3=k_3\left(n_{0,3}+{\tilde{x}}_2-x_3-{\tilde{\chi }}_{\mathrm{o}}\right),\hfill & \hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{\mathit{x}}^{\left(3\right)}\left(0\right)& =x_{0,3}=0,\hfill \end{array}$$

(110)

$$\begin{array}{cccc}\hfill {\dot{\mathit{x}}}^{\left(4\right)}& =\left[\begin{array}{c}{\dot{x}}_4\\ {\dot{x}}_5\end{array}\right]=\left[\begin{array}{c}{k}_4\left(n_{0,3}+{\tilde{x}}_2-{\tilde{x}}_3-{\tilde{\chi }}_{\mathrm{o}}\right)\\ k_5\left(n_{0,5}+{\tilde{\chi }}_{\mathrm{o}}\right)\end{array}\right],\hfill & \hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{\mathit{x}}^{\left(4\right)}\left(0\right)& =\left[\begin{array}{c}{x}_{0,4}\\ x_{0,5}\end{array}\right]=\mathbf{0},\hfill \end{array}$$

(111)

$$\begin{array}{cc}\hfill {\chi }_{\mathrm{o}}^{\left(4\right)}& ={\chi }_{\mathrm{o}}=x_4-x_5,\hfill \end{array}$$

(112)

The information flows in these subsystems are shown in the inset of Figure 12.

Incremental parameter estimation

The resulting partitioning means that estimates of the $k_1$, $k_2$, and $k_3$ can be obtained by solving the following three single-parameter estimation problems ($P_1^{\left(1\right)}$, $P_1^{\left(2\right)}$, $P_1^{\left(3\right)}$):

$$\begin{array}{cc}\hfill P_1^{\left(1\right)}:{\widehat{k}}_1& =arg\underset{k_1}{min}\sum _h^H{\left({\tilde{x}}_{1,h}-x_1\left(t_h\right)\right)}^{\mathrm{T}}{W}^{\left(1\right)}\left({\tilde{x}}_{1,h}-x_1\left(t_h\right)\right)\hfill \end{array}$$

(113)

subject to

$$\begin{array}{cccc}\hfill {\dot{x}}^{\left(1\right)}& ={\dot{x}}_1=k_1\left(n_{0,1}-x_1-2{\tilde{x}}_2\right),\hfill & \hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{x}^{\left(1\right)}\left(0\right)& =x_{0,1}=0\hfill \end{array}$$

(114)

where

$$\begin{array}{cc}\hfill W^{\left(1\right)}& :=5/4\hfill \end{array}$$

(115)

$$\begin{array}{cc}\hfill P_1^{\left(2\right)}:{\widehat{k}}_2& =arg\underset{k_2}{min}\sum _h^H{\left({\tilde{x}}_{2,h}-x_2\left(t_h\right)\right)}^{\mathrm{T}}{W}^{\left(2\right)}\left({\tilde{x}}_{2,h}-x_2\left(t_h\right)\right)\hfill \end{array}$$

(116)

subject to

$$\begin{array}{cccc}\hfill {\dot{x}}^{\left(2\right)}& ={\dot{x}}_2=k_2\left(n_{0,1}-{\tilde{x}}_1-2x_2\right),\hfill & \hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{x}^{\left(2\right)}\left(0\right)& =x_{0,2}=0\hfill \end{array}$$

(117)

where

$$\begin{array}{cc}\hfill W^{\left(2\right)}& :=5/3\hfill \end{array}$$

(118)

$$\begin{array}{cc}\hfill P_1^{\left(3\right)}:{\widehat{k}}_3& =arg\underset{k_3}{min}\sum _h^H{\left({\tilde{x}}_{3,h}-x_3\left(t_h\right)\right)}^{\mathrm{T}}{W}^{\left(3\right)}\left({\tilde{x}}_{3,h}-x_3\left(t_h\right)\right)\hfill \end{array}$$

(119)

subject to

$$\begin{array}{cccc}\hfill {\dot{x}}^{\left(3\right)}& ={\dot{x}}_3=k_3\left(n_{0,3}+{\tilde{x}}_2-x_3-{\tilde{\chi }}_{\mathrm{o}}\right),\hfill & \hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{x}^{\left(3\right)}\left(0\right)& =x_{0,3}=0\hfill \end{array}$$

(120)

where

$$\begin{array}{cc}\hfill W^{\left(3\right)}& :=5/4\hfill \end{array}$$

(121)

The fourth set of parameters, $k_4$ and $k_5$, are estimated by solving the following problem:

$$\begin{array}{cc}\hfill P_1^{\left(4\right)}:{\widehat{k}}_3,{\widehat{k}}_4& =arg\underset{k_3,k_4}{min}\sum _h^H{\left({\tilde{\chi }}_{\mathrm{o},h}-{\chi }_{\mathrm{o}}\left(t_h\right)\right)}^{\mathrm{T}}{W}^{\left(4\right)}\left({\tilde{\chi }}_{\mathrm{o},h}-{\chi }_{\mathrm{o}}\left(t_h\right)\right)\hfill \end{array}$$

(122)

subject to

$$\begin{array}{cccc}\hfill {\dot{\mathit{x}}}^{\left(4\right)}& =\left[\begin{array}{c}{\dot{x}}_4\\ {\dot{x}}_5\end{array}\right]=\left[\begin{array}{c}{k}_4\left(n_{0,3}+{\tilde{x}}_2-{\tilde{x}}_3-{\tilde{\chi }}_{\mathrm{o}}\right)\\ k_5\left(n_{0,5}+{\tilde{\chi }}_{\mathrm{o}}\right)\end{array}\right],\hfill & \hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{\mathit{x}}^{\left(4\right)}\left(0\right)& =\left[\begin{array}{c}{x}_{0,4}\\ x_{0,5}\end{array}\right]=\mathbf{0}\hfill \end{array}$$

(123)

$$\begin{array}{cc}\hfill {\chi }_{\mathrm{o}}^{\left(4\right)}& ={\chi }_{\mathrm{o}}=x_4-x_5,\hfill \end{array}$$

(124)

where

$$\begin{array}{cc}\hfill W^{\left(4\right)}& :=5/4\hfill \end{array}$$

(125)

Please note that the (unknown) measurement error variance $\sigma$ can arbitrarily be set equal to one during parameter estimation without affecting the parameter estimates. This, however, means that parameter uncertainties cannot be quantified in this case without estimating $\sigma$ as well. This estimation and subsequent uncertainty analysis is considered out of scope for this study. Figure 11 shows the simulated extent profiles obtained after solving the optimization problems ($P_1^{\left(j\right)},j=1,\dots ,4$). There is a reasonable fit in all cases, except for the third extent of reaction ($x_3$). This, in the mind of the authors, demonstrates a tangible benefit of the extent-based model framework. Indeed, one can identify which parts of the model could be improved, e.g., by reconsidering the reaction rate expression corresponding to the computed extents that are approximated poorly. However, this is not explored further in this work.

Figure 10 shows the measured concentrations and the simulated profiles obtained with (i) incremental parameter estimation ($P_1^{\left(j\right)},j=1,\dots ,4$), followed by (ii) simultaneous parameter estimation with estimates obtained through incremental parameter estimation ($P_{1+0}$). The parameter values obtained by solving $P_0$, $P_1^{\left(1\right)}$, $P_1^{\left(2\right)}$, and $P_{1+0}$ are shown in

Table 3. $\alpha$-pinene—Parameter estimates. The parameters obtained with simultaneous estimation ($P_0$) and incremental estimation ($P_{1+0}$) are practically the same.

	Name	Unit	${\mathit{P}}_0$	${\mathit{P}}_1^{\left(\mathit{j}\right)}$	${\mathit{P}}_{1+0}$
${\theta }_1$	$k_1$	$\%{\mathrm{h}}^{-1}$	0.213	0.214	0.213
${\theta }_2$	$k_2$	$\%{\mathrm{h}}^{-1}$	0.107	0.106	0.107
${\theta }_3$	$k_3$	$\%{\mathrm{h}}^{-1}$	0.074	0.074	0.074
${\theta }_4$	$k_4$	$\%{\mathrm{h}}^{-1}$	0.989	1.037	0.989
${\theta }_5$	$k_5$	$\%{\mathrm{h}}^{-1}$	0.144	0.148	0.144
$\sqrt{\frac{Q_0}{H·M}}$	WRMSR	−	$0.66$	$0.67$	$0.66$

. Note the final values obtained with both $P_0$ and $P_{1+0}$ are equal within numerical precision.

5. Discussion

A novel procedure for model parameter estimation is proposed. The procedure combines two important features: (i) a new framework for extent computation based on the computation of observable directions among the set of sensed but ambiguous extents; and (ii) a graph-based system partitioning procedure that identifies the groups of equations and kinetic parameters that can be simulated independently from the remaining part of the dynamic model. The latter is possible by approximating the original equation system via interpolation of observed extents and extent directions. The proposed procedure enables splitting the parameter estimation problem into smaller problems in cases that could so far only be handled with simultaneous model identification methods (e.g., scenario A, C, and D in Section 4.1). It also subsumes preexisting methods under special conditions such as complete observability (e.g., scenario E in Section 4.1,5,6,7]) or absence of ambiguity (e.g., scenario B in Section 4.1,8]). Thus, this procedure is applicable to parameter estimation for any multivariate differential equation model under a large range of structural observability conditions.

This work removes an important bottleneck for the implementation of incremental model identification to biological systems. Indeed, by providing a method that can be applied to any known stoichiometric matrix $\mathbf{N}$ and measurement matrix $\mathbf{M}$, incremental model identification is now applicable to biological processes where the number of measurements M is lower than the number of modeled reactions R. Computationally speaking, the extent computation and graph-based partitioning are expected to scale well with the number of modeled reactions. The effects on the computational efficiency of the parameter estimation step will however depend primarily on the lengths of the state vectors in the identified subsystems. In turn, these lengths depend strongly on model structure and the available measured variables and less so on the number of states in the complete system.

A remaining obstacle is the fact that both $\mathbf{M}$ and $\mathbf{N}$ are assumed to be known. While several techniques exist to handle this, including methods based on extents [17], some coefficients in these matrices may be unidentifiable based on atomic and stoichiometric balances alone. For this reason, dealing with unknown stoichiometry ($\mathbf{M}$) or ill-defined measurements ($\mathbf{N}$) deserves continued attention.

Another aspect open for exploration is that no global optimization method has been tested so far to solve the smaller parameter estimation problems $P_1^{\left(j\right)}$ that result from system partitioning. Inspiration can be drawn from existing parameter estimation methods even if these methods (i) might be complex yet without publicly available implementation [18]; or (ii) remain limited to cases where all extents are observable and therefore modeled as univariate processes for parameter estimation [19,20].

5.1. Generalized Framework for Extent Computation

The addition of new concepts, such as ambiguous extents and observable extent directions expands the framework for extent computation beyond its original range of applicability. In this study, the focus has been given to the problem of parameter estimation. However, the same generalized framework is likely to be applicable and useful for challenges that have been handled with earlier extent-based methods such as data reconciliation [7], model structure identification [6], and process control [21]. Moreover, there are no obvious limitations that hinder applications involving multi-phase or distributed system models [22,23] or models including algebraic constraints [6,22].

5.2. Optimal System Partitioning

According to our understanding, the proposed procedure delivers the most efficient system partitioning given the available measurements and model information $\mathbf{M}$, $\mathit{f}\left(·\right)$, $\mathbf{N}$, and ${\mathit{n}}_0$. In other words, the procedure generates the subsystems with the smallest number of parameters that can be estimated independently from parameters in other subsystems. However, this is stated here without formal proof. The factorization of ${\mathbf{G}}_{•,\mathbf{a}}$ based on the reduced row echelon form of $\mathbf{G}$, as described in Section 3.3.1, plays a crucial role. Indeed, it follows from work on subspace clustering methods for noise-free data [24,25] that this factorization minimizes the presence of simulation arcs from the ambiguous extents to the observable directions, thereby leading to the best possible system partitioning.

Please note that the obtained generalized extent framework is optimal only for the purpose of parameter estimation. For instance, one could consider a factorization of ${\mathbf{G}}_{•,\mathbf{a}}$ which identifies directions that exhibit uncorrelated estimation errors, thereby improving the statistical quality of the computed observable directions and possibly avoiding Type A approximation error. This could be achieved via singular value decomposition of ${\mathbf{G}}_{•,\mathbf{a}}$. However, singular value decomposition cannot guarantee optimal system partitioning [26]. Consequently, it follows that the optimal factorization of ${\mathbf{G}}_{•,\mathbf{a}}$ depends on the objective of this factorization. This stands in contrast to the existing body of work concerning extents, which do not involve any such purpose-dependent factorization steps.

Note also that the proposed procedure is designed for the definition of optimality used here. For instance, if one relaxes the constraint that parameters can only appear in one subsystem (Section 3.6), then the proposed graph partitioning procedure is not optimal. For example, in scenario A, one may be able to estimate the parameters in 3 subsets rather than 2, namely, (i) $k_1$, $k_2$, and $K_1$ by simulation of $x_1\left(t\right)$ and $x_2\left(t\right)$; (ii) $k_2$, $k_3$, $k_4$, and $k_5$ by simulation of $x_2\left(t\right)$, $x_4\left(t\right)$, $x_5\left(t\right)$, and ${\chi }_{\mathrm{o},1}\left(t\right)$; and (iii) $k_3$ by simulation of $x_3\left(t\right)$. This leads to 3 parameter estimation problems with 3, 4, and 1 parameters, which may be easier to solve than the 2 parameter estimation problems with 5 and 1 parameters as obtained with the proposed procedure. It is however non-trivial to identify a procedure, including both factorization of ${\mathbf{G}}_{•,\mathbf{a}}$ and system partitioning, that guarantees the best partitioning when applying this relaxed definition of optimality, nor is it clear whether such a procedure exists.

5.3. New Opportunities

This work also hints at several new applications of the generalized framework for extent computation:

(a) 

Identifiability analysis. While not a core objective of this work, it has been shown that the graph partitioning method can help identify unidentifiable parameters. Given that this labeling is based on a model reformulation and does not depend on the temporal resolution or quality of the collected measurements, it follows that the method identifies structurally unidentifiable parameters. Unlike other methods [27], the proposed approach does not require symbolic differentiation. It remains to be explored whether this can also be used to positively identify structurally identifiable parameters. For a discussion on the evaluation and use of indicators of structural and practical parameter identifiability we refer to [28].

(b) 

Soft-sensing. The appearance of observable directions among the ambiguous extents is closely related to the concepts of observability and detectability in the context of state estimation [29]. Note however that the observability labels in this work are based on the stoichiometric balances and measurement equations alone, thus excluding the dynamic model. At the same time, it is suspected that the extents corresponding to vertices that are not on directed paths to vertices representing observable extents or extent directions can be labeled as structurally unobservable, again observing that the exact timing and quality of the measurements does not play any role in this labeling. Similarly to the identifiability analysis discussed above, such an approach would not rely on symbolic differentiation. Whether this can be used to unambiguously determine observability and detectability for all model states remains to be studied.

(c) 

Experimental design. The labeling of extents and directions as observable or unobservable suggests that experimental design may be used to optimize the selection of measured variables. A method to do so has been applied in [30] to enumerate all Pareto-optimal flow sensor layouts in wastewater treatment plants. In [31], symbolic computation enabled the identification of optimal experimental designs. Similar approaches could be applied as a measurement selection method for metabolic flux analysis and the monitoring of complex systems.

6. Conclusions

In this work, an incremental parameter identification procedure has been developed and tested based on a generalized extent-based framework. This generalized framework enables incremental parameter estimation in cases where the previous methods based on the computation of extents did not permit this. Importantly, through study of simulated and experimental medium-sized examples, the generality of the developed procedure has been demonstrated and new opportunities offered by this framework have been identified.