Robust Model Selection: Flatness-Based Optimal Experimental Design for a Biocatalytic Reaction

Abstract

Considering the competitive and strongly regulated pharmaceutical industry, mathematical modeling and process systems engineering might be useful tools for implementing quality by design (QbD) and quality by control (QbC) strategies for low-cost but high-quality drugs. However, a crucial task in modeling (bio)pharmaceutical manufacturing processes is the reliable identification of model candidates from a set of various model hypotheses. To identify the best experimental design suitable for a reliable model selection and system identification is challenging for nonlinear (bio)pharmaceutical process models in general. This paper is the first to exploit differential flatness for model selection problems under uncertainty, and thus translates the model selection problem to advanced concepts of systems theory and controllability aspects, respectively. Here, the optimal controls for improved model selection trajectories are expressed analytically with low computational costs. We further demonstrate the impact of parameter uncertainties on the differential flatness-based method and provide an effective robustification strategy with the point estimate method for uncertainty quantification. In a simulation study, we consider a biocatalytic reaction step simulating the carboligation of aldehydes, where we successfully derive optimal controls for improved model selection trajectories under uncertainty.

1. Introduction

Quality by design (QbD) and quality by control (QbC) are critical aspects in pharmaceutical manufacturing, aiming for low-cost but high-quality drugs [1]. For research and development costs in the competitive and strongly regulated pharmaceutical industry have been rising over recent decades leading to shrinking profit margins [2], QbD and QbC have attracted much notice. Mathematical models and process systems engineering might be useful tools for implementing effective QbD/QbC strategies, which is particularly true for the rapidly growing biopharmaceutical market [3,4,5]. A crucial task in modeling biopharmaceutical manufacturing processes is the reliable identification of model parameters and proper model candidates from a set of various model hypotheses. Model-based design of experiment (MBDoE) approaches, where the challenge of precise parameter identification and reliable model selection is formulated as an optimization problem, have proven beneficial over the last few decades [6,7,8,9]. Please note that the literature distinguishes between statistical design of experiment and MBDoE approaches. Methods of the former as response surface methods are popular and commonly used because of their simplicity, but, for the same reason, are not able to sufficiently handle complex dynamic systems [10]. On the other hand, the MBDoE approach frequently aims for optimal control actions, where commonly used numerical methods give rise to nonlinear programming problems that need substantial amounts of computational time and efficient solvers [11]. Moreover, robustification against model parameter uncertainties further complicates the problem, as statistical quantities must be considered when solving the underlying optimization problem. In addition to the actual MBDoE optimization framework and data acquisition, recent studies have shown that a well-posed parameter identification and model selection problem is equally important [12,13,14]. In this context, the consideration of input residuals in parameter estimation beyond the classical approach of output residuals has drawn attention in the literature [13,14,15,16]. For an effective solution to the robust model selection problem, we aim at implementing an inverse modeling technique by using the differential flatness concept [17]. In systems theory, differential flatness implies that analytical expressions of the system’s states and controls exist that are functions of so-called flat outputs and their derivatives. Hence, by implementing a differential flatness strategy to solve an MBDoE problem, we can avoid solving differential and sensitivity equations numerically as required by standard control vector parameterization strategies for model selection. Please also note that while optimizing the flat output trajectory, the number of continuous optimization variables is reduced to the number of flat outputs, i.e., an additional reduction of computational costs as compared to standard nonlinear programming techniques with control and state variable discretization [18]. Methods based on differential flatness have been heavily used in research and industry—primarily in the design of open-loop controllers and the planning of state trajectories in (electro)mechanical systems [19]. In the case of the MBDoE technique, strategies exploiting differential flatness have barely been considered [15,16,20,21], and—to our knowledge—have not been applied to model selection problems under uncertainty before. Thus, this study addresses two open research challenges: (i) to apply the flatness concept to multi-model problems, and (ii) to integrate parameter uncertainties into the flatness-based optimization problem for robust model selection. In detail, we aim at discarding inappropriate model candidates from a set of competing model hypotheses. To this end, we show that the considered model candidates satisfy the differential flatness condition, and we solve an optimization problem under uncertainty that makes use of the flatness property. For the sake of illustration, we consider a biocatalytic process from a carboligation reaction system that forms an essential precursor in the synthesis of natural products and pharmaceuticals [22]. Please note that carbon-carbon bond-forming reactions are the backbone of numerous high-value molecules in industrial organic synthesis. However, biocatalytic production of the related bulk and fine chemicals and active pharmaceutical ingredients (APIs) is insufficiently explored today. Therefore, the application of enzymes in organic synthesis is currently an important research topic to exploit the significant potential of improving manufacturing processes in the chemical and pharmaceutical industry under the agenda of green chemistry [23,24]. An essential problem with all mathematical models and model-based design concepts, including MBDoE, is the uncertainty in model parameters. To ensure reliable inferences, parametric uncertainty should be taken into account in the MBDoE approach [6,14,25,26]. In general, probability-based concepts have attracted considerable attention in robust optimization. In the literature, between classical and Bayesian settings for handling model predictive errors as a potential consequence of uncertain model parameters can be distinguished. In the former, the strategies have in common that the model output differences in the objective function are weighted by an error covariance matrix of the model prediction to avoid experimental regions where model predictions are poor [7,27,28,29]. One drawback of these strategies is that the prediction error covariance matrix is calculated from local sensitivities, where designs of experiments for local sensitivities might provide misleading information compared to designs based on global sensitivities [30]. On the other hand, Bayesian model selection with consideration of parametric uncertainties requires the determination of computationally expensive integrals, in particular, if the recommended bias-free numerical evaluation path is chosen [31,32]. More broadly speaking, uncertainty propagation and quantification still pose considerable challenges in computational efficiency and numerical accuracy. In this context, we propose to diminish this computational burden by using an efficient and effective rule to determine statistical moments circumventing a local approximation, namely the point estimate method (PEM). The PEM has proven beneficial in various engineering problems [33], including robust optimization problems in complex (bio)chemical process design [3,34].

The paper is organized as follows. In Section 2, we cover the basics of the inverse modeling and differential flatness concept. In Section 3, we introduce the flatness-based MBDoE technique with the PEM for its robustification under parameter uncertainties. In Section 4, we apply the proposed MBDoE concept to a biocatalytic reaction network, and we discuss the results in Section 5. Conclusions can be found in Section 6.

2. Inverse Modeling

In process systems engineering, dynamic process models are often mathematically described with ordinary differential equations (ODEs). Their state-space formulation is

$$\frac{\mathrm{d}\mathbf{x}}{\mathrm{d}t}=\mathbf{f}(\mathbf{x}\left(t\right),\mathbf{u}\left(t\right),\mathsf{\theta }),$$

(1)

where t is the time, $\mathbf{x}\left(t\right)\in {\mathbb{R}}^n$ is the system states, $\mathbf{u}\left(t\right)\in {\mathbb{R}}^m$ is the system inputs, and $\mathsf{\theta }\in {\mathbb{R}}^p$ is the parameter vector. The nonlinear function $\mathbf{f}$: ${\mathbb{R}}^n\times {\mathbb{R}}^m\times {\mathbb{R}}^p\to {\mathbb{R}}^n$ represents the vector field of the dynamic system. The output function is defined as

$$\mathbf{y}\left(t\right)=\mathbf{h}\left(\mathbf{x}\right(t),\mathbf{u}(t),\mathsf{\theta }),$$

(2)

with the system output vector $\mathbf{y}\left(t\right)\in {\mathbb{R}}^q$ and $\mathbf{h}$: ${\mathbb{R}}^n\times {\mathbb{R}}^m\times {\mathbb{R}}^p\to {\mathbb{R}}^q$ as the (non-)linear output function. Please note that in many cases, the system outputs are a (linear) function of the system states only; that is, $\mathbf{y}\left(t\right)=\mathbf{h}\left(\mathbf{x}\right(t\left)\right)$. Equations (1) and (2) form the process model, $\mathcal{M}$. Please note that in the following, time and parameter dependencies are specified only when the context requires it.

As opposed to solving the ODE system (Equation (1)) in a forward manner to conventionally exploit the outputs $\mathbf{y}$ given the inputs $\mathbf{u}$, inverse modeling aims at reconstructing the inputs $\mathbf{u}$ for given outputs $\mathbf{y}$ [35,36,37]. Please note that the process of inverse simulation is not to be confused with inverse problems, where in the latter the goal is to estimate unknown model parameters $\widehat{\mathsf{\theta }}$ [35]. The different principles of forward simulation, inverse modeling, and the inverse problem of parameter identification are visualized in Figure 1.

In addition to numerical inversion approaches as model inversion- or inverse simulation-based techniques [38], algebraic methods have been proposed for inverse modeling [39], among which strategies based on differential flatness are important representatives.

2.1. Differential Flatness

In systems theory, within the context of algebraic model inversion techniques and aiming at feedforward control problems, differential flatness was introduced by Fliess et al. [17] in 1992. A (non-)linear process model (Equations (1) and (2)) is called differentially flat or shortly flat if there exists an output vector

$${\mathbf{y}}^{\mathrm{flat}}={\mathbf{h}}^{\mathrm{flat}}(\mathbf{x},\mathbf{u},\dot{\mathbf{u}},\cdots ,{\mathbf{u}}^{\left(s\right)},\mathsf{\theta }),$$

(3)

with a finite value $s\in \mathbb{N}$ and the smooth mapping function ${\mathbf{h}}^{\mathrm{flat}}:{\mathbb{R}}^n\times {\left({\mathbb{R}}^m\right)}^{s+1}\times {\mathbb{R}}^p\to {\mathbb{R}}^q$, referred to as a flat output which fulfills the following conditions:

• The states and the controls can be described as a function of the flat output and its derivatives:

  $$\begin{array}{cc}\hfill \mathbf{x}& ={\mathsf{\Psi }}_x({\mathbf{y}}^{\mathrm{flat}},{\dot{\mathbf{y}}}^{\mathrm{flat}},\cdots ,{{\mathbf{y}}^{\mathrm{flat}}}^{\left(r\right)},\mathsf{\theta }),\hfill \end{array}$$

  (4)

  $$\begin{array}{cc}\hfill \mathbf{u}& ={\mathsf{\Psi }}_u({\mathbf{y}}^{\mathrm{flat}},{\dot{\mathbf{y}}}^{\mathrm{flat}},\cdots ,{{\mathbf{y}}^{\mathrm{flat}}}^{(r+1)},\mathsf{\theta }),\hfill \end{array}$$

  (5)

  with the smooth mapping functions ${\mathsf{\Psi }}_x:{\left({\mathbb{R}}^m\right)}^{r+1}\times {\mathbb{R}}^p\stackrel{}{\to }{\mathbb{R}}^n$ and ${\mathsf{\Psi }}_u:{\left({\mathbb{R}}^m\right)}^{r+2}\times {\mathbb{R}}^p\stackrel{}{\to }{\mathbb{R}}^m$.
• The dimensions of the control and the flat output vector are equal:

  $$\dim {\mathbf{y}}^{\mathrm{flat}}=\dim \mathbf{u}.$$

  (6)

Here, $r\in \mathbb{N}$ specifies the number of occurring derivatives. Generally speaking, r is not known beforehand, apart from single-input single-output (SISO) systems, where $r=n-1$ [40]. There is an infinite number of flat output candidates, and often, they are a function of the states only [40]—similar to the fact that the real output function of a dynamic system is, in many cases, exclusively a function of the states.

Flat outputs might have no direct physical meaning, or they might be identical to measurable quantities; then also referred to as flat inputs [41]. Moreover, between local and global flatness must be distinguished, where for the local form, the differential flatness is valid only for a restricted domain bound by occurring singularities [42]. For certain kinds of singularities, a globally flat system can be designed following a superposition concept. To this end, individual local model inversion domains are created from different local flat outputs where the union of the local domains provides a global framework for the flatness-based model inversion [42]. Alternatively, a transformation approach exists to bypass points of a singularity [43]. In addition to differential flatness of ODE systems, flatness-based methods are also being used for discrete-time systems and partial differential equations (PDEs); see [44,45,46], and references within. Application scenarios based on differential flatness cover predominantly (electro)mechanical problems, but also process technologies like heat exchangers [47], crystallizers [48] or (bio)reactors [45,49], lithium-ion batteries [50], and more.

2.2. Flat Output Identification

To implement a flatness-based concept, two key challenges exist: First, determining if a system is differentially flat, and second, constructing a flat output alternatively. There is no general method for either of those challenges. However, exceptions exist, and include systems that are linearizable by static feedback control, which are flat by definition. Considering a linear system, controllability conditions differential flatness and vice versa [43]. Methods of determining flat outputs and the mapping function ${\mathbf{h}}^{\mathrm{flat}}$, respectively, are currently being researched, and the interested reader is referred to [19,51,52] and the references therein. Moreover, the construction of so-called flat inputs (i.e., reconstructed inputs based on the given output function (Equation (2))) was shown for SISO systems and a limited set of multi-input and multi-output (MIMO) cases [53].

Alternatively, in systems theory, valuable information about system characterization based on the interaction of inputs, system states and outputs can also be extracted from graph theory [39,54,55]. Typically, observability or controllability measures are derived for (non)linear state-space models using directed graphs (digraphs) [56,57,58]. A digraph $D=[V,E]$ with $n+m+q$ vertices $v_i\in V$ and edges $e\in E$ can be constructed from adjacency matrices $A_u$, $A_x$, and $A_y$ of the process model for the inputs, system states, and outputs, respectively. The $a_{i,j}$-th element of $A_x$ or $A_u$ is set to 1 if derivatives $\frac{\partial f_i(\mathbf{x},\mathbf{u})}{\partial x_j}$ or $\frac{\partial f_i(\mathbf{x},\mathbf{u})}{\partial u_j}$ exist, respectively, and to 0 if this is not the case. If these derivatives exist (i.e., $\partial f_i/\partial x_j\ne 0$ or $\partial f_i/\partial u_j\ne 0$), then there is an edge from vertex $v_j$ to vertex $v_i$. Analogously, the adjacency matrix for the outputs $\mathbf{y}$ and the output function $h(\mathbf{x},\mathbf{u})$ can be defined. Regarding the flatness property, in the paper by Csercsik et al. [55], an algorithm for flatness analysis of MIMO systems is introduced that uses an explicit expressibility graph. Please note that from a given digraph, the explicit expressibility graph is formed with reversed edge directions, with the outputs replaced by the flat outputs, and with the self-loop system states dependencies omitted. The following theorem then gives a necessary condition for differential flatness based on the adjacency matrices $A_u$, $A_x$, and $A_y$ [55].

Theorem 1.

For a prospective set of flat output-input pairs, m pairwise disjoint paths must exist in the explicit expressibility graph, the union of which covers each vertex of the graph.

In summary, existing construction methods are often restrictive and require cumbersome calculations. In practice, for the identification of flat outputs (Equation (3)) a sequential strategy including an expert guessing for a flat output candidate followed by a subsequent validation step in terms of Equations (4)–(6) has proven favorable. This trial-and-error approach is facilitated by two facts. First, numerous technical systems are indeed flat systems. Second, frequently, comparable to Lyapunov functions in control theory, informative flat outputs have physical signification [40]. Once the flat output configuration is identified, the flat output functions $y^{\mathrm{flat}}$ must be parameterized using for instance spline functions.

2.3. Basis Splines

The flat outputs ${\mathbf{y}}^{\mathrm{flat}}$ have to be parameterized while limited by the fact that they have to satisfy solutions of the dynamic system (Equations (1) and (2)). For flat outputs, basis functions from a large set of possibilities, ranging from simple polynomials to more complex functions, can be chosen. Basis splines, i.e., piece-wise polynomial functions also known as B-splines, offer great freedom of action, are flexible and well-studied, and libraries for their usage are implemented in many different programming languages [59]. The classical application of splines is the approximation and interpolation of data points.

A flat output based on B-splines is defined by

$$y^{\mathrm{flat}}\left(t\right)=\sum _{i=1}^l{\theta }_i^{🟉}{N}_{i,k}\left(t\right),$$

(7)

where $N_{i,k}$ are B-splines of order k on l collocation points, and ${\theta }_i^{🟉}$ is a coefficient vector, whose elements are referred to as meta-parameters. Please note that the meta-parameters ${\theta }_i^{🟉}$ constitute the decision variables in the flatness-based optimal experimental design framework for model selection, which is introduced in the following section.

3. Robust Design of Experiments for Model Selection

To maximize the information provided by an experiment for the aim of discriminating between different model hypotheses, we incorporate the inverse modeling approach into the design of experiments setting. First, a generic scheme for a standard MBDoE strategy for model selection is given in Figure 2. Starting with the model-building process, several competing model hypotheses and candidates may exist. Thus, the MBDoE technique predicts operating conditions and experimental data that ensure better model selection and falsification. Next, the optimized experimental setup is implemented, and new informative data are recorded. With these new data, the model parameters of all model candidates are refined and updated. Finally, based on the recalculated mismatch of experimental data and simulation results, the difference in the competing model candidates may become more significant, and some model candidates might be excluded in terms of falsification. As indicated in Figure 2, the MBDoE workflow is an iterative process; that is, the prediction of improved experiments, the data acquisition, the adaption of the model parameters, and the assessment of the remaining model candidates is repeated until the most suitable model candidate is identified.

Mathematically, the MBDoE approach can be stated as an optimal control problem (OCP) of the form

$$\begin{array}{c}\hfill \begin{array}{ccc}\hfill & \underset{\mathsf{\varphi }\left(t\right)\in \mathsf{\Phi }}{\max }\hfill & \mathcal{J}\left(\mathsf{\varphi }\left(t_{\mathrm{f}}\right)\right)\hfill \\ \hfill & \mathrm{s}.\mathrm{t}.\hfill & \mathrm{Process}\mathrm{model}:\mathrm{Equation}\left(1\right),\hfill \\ \hfill & \hfill & \mathbf{x}\left(t_0\right)={\mathbf{x}}_0,\hfill \\ \hfill & \hfill & {\mathbf{g}}_{\mathrm{eq}}(\mathbf{x}\left(t\right),\mathbf{u}\left(t\right),\mathsf{\theta })=0,\hfill \\ \hfill & \hfill & {\mathbf{g}}_{\mathrm{ineq}}(\mathbf{x}\left(t\right),\mathbf{u}\left(t\right),\mathsf{\theta })\le 0,\hfill \\ \hfill & \hfill & {\mathsf{\varphi }}_{\min }\le \mathsf{\varphi }\le {\mathsf{\varphi }}_{\max }.\hfill \end{array}\end{array}$$

(8)

$\mathcal{J}$ is the scalar objective function aiming to maximize the difference between the trajectories of the output functions of competing model candidates, and $\mathsf{\varphi }\left(t\right)=\mathsf{\varphi }(\mathbf{u}\left(t\right),{\mathbf{x}}_0)$ is the design vector within the design space $\mathsf{\Phi }$ specified by the control $\mathbf{u}\left(t\right)$ and the initial states ${\mathbf{x}}_0$. Please note that $t\in [t_0,t_0+t_{\mathrm{f}}]$, where $t_0=0$ is the initial time, and $t_{\mathrm{f}}$ is the time duration of the experiment. The objective function is given in Mayer form without loss of generality as the Lagrange and Bolza problem types can be converted to the former [18]. Both equality constraints ${\mathbf{g}}_{\mathrm{eq}}:{\mathbb{R}}^{n\times m\times p}\to {\mathbb{R}}^{n_{\mathrm{eq}}}$ and inequality constraints ${\mathbf{g}}_{\mathrm{ineq}}:{\mathbb{R}}^{n\times m\times p}\to {\mathbb{R}}^{n_{\mathrm{ineq}}}$ can be considered in the optimization problem. Exemplarily, they might be introduced to avoid negative values of system states or to restrict a species’ concentration to a certain limit. [${\mathsf{\varphi }}_{\min }$,${\mathsf{\varphi }}_{\max }$] are the upper and lower boundaries for the design variables.

Solving the dynamic optimization problem (Equation (8)) can be performed following two different discretization strategies of the time-dependent functions in Equation (8): Optimize then Discretize, also called the indirect or variational approach, or Discretize then Optimize, its direct counterpart. The first strategy focuses on the solution of the first-order optimality conditions for the OCP, resulting in a boundary value problem (BVP). The difficulties associated with solving BVPs led to research activities dealing with the conversion of the OCP to constrained, finite-dimensional problems to exploit state-of-the-art large-scale nonlinear program (NLP) solvers [18]. Methods applying NLP solvers can be classified further into sequential and simultaneous strategies. In the sequential approach referred to as single shooting, the control vector is parameterized, and the resulting NLP problem is solved using control vector parameterization (CVP) methods [60,61]. In contrast, the simultaneous approach multiple shooting divides the OCP into smaller subproblems requiring next to CVP initial states in the individual problems and serves as a bridge to the direct transcription approach, where all variables, i.e., states and controls, are discretized [18]. The discretization of the states and controls is usually realized using collocation and control parameterization techniques that lead to complex optimization problems for which efficient solvers and great computational power are needed [11].

Various scalar objective functions $\mathcal{J}$ for the MBDoE approach (Equation (8)) have been proposed in the literature [7,62,63,64,65,66]. The general goal of the optimal control problem for model selection (8) is to maximize some measure on the distance between the different model candidates by searching for an optimal control profile while meeting the specified constraints, including the solving of the underlying ODE system (Equation (1)).

The proposed flatness-based approach for the MBDoE performs differently. It literally operates in an inverse manner. Optimized flat output trajectories maximize the difference between the model outputs of the competing model candidates. The related controls $\mathbf{u}$ and the initial states ${\mathbf{x}}_0$ are reconstructed from the flat model ${\mathcal{M}}^{-1}$ and have to be identical for all model candidates and their optimized flat outputs ${\mathbf{y}}^{\mathrm{flat}}$. Using the Euclidean distance as a discrepancy measure between model outputs and assuming M model candidates, the original OCP given in Equation (8) becomes

$$\begin{array}{c}\hfill \begin{array}{ccc}\hfill & \underset{{\mathbf{y}}^{\mathrm{flat}}}{\max }\hfill & \underset{t_0}{\overset{t_{\mathrm{f}}}{\int }}\sum _{i=1}^{M-1}\sum _{j=i+1}^M[{\mathbf{y}}^i({\mathbf{y}}^{\mathrm{flat},i},{\dot{\mathbf{y}}}^{\mathrm{flat},i},\cdots ,{{\mathbf{y}}^{\mathrm{flat},i}}^{\left(r\right)})\hfill \\ \hfill & \hfill & -{\mathbf{y}}^j({\mathbf{y}}^{\mathrm{flat},j},{\dot{\mathbf{y}}}^{\mathrm{flat},j},\cdots ,{{\mathbf{y}}^{\mathrm{flat},j}}^{\left(r\right)}){]}^2\mathrm{d}t\hfill \\ \hfill & \mathrm{s}.\mathrm{t}.\hfill & \Delta \mathbf{u}\left({\mathbf{y}}^{\mathrm{flat}}\right)=0,\hfill \\ \hfill & \hfill & \Delta {\mathbf{x}}_0\left({\mathbf{y}}^{\mathrm{flat}}\right)=0,\hfill \\ \hfill & \hfill & {\mathbf{g}}_{\mathrm{eq}}\left({\mathbf{y}}^{\mathrm{flat}}\right)=0,\hfill \\ \hfill & \hfill & {\mathbf{g}}_{\mathrm{ineq}}\left({\mathbf{y}}^{\mathrm{flat}}\right)\le 0.\hfill \end{array}\end{array}$$

(9)

The inputs, states and outputs in the optimization formulation are expressed according to Equations (2), (4), and (5). The differences between the recalculated inputs and the initial states for all M model candidates are measured by the delta function $\Delta (·)$ that might be the Euclidean distance alike. For the flat outputs are of functional form, time-dependent empirical basis functions ${\mathbf{y}}^{\mathrm{flat},i}={\mathbf{y}}^{\mathrm{flat},i}({\mathsf{\theta }}^{🟉,i},t)\forall i\in \{1,\dots ,M\}$ with meta-parameters ${\mathsf{\theta }}^{🟉,i}$ have to be specified; see also Section 2.3. The optimal control problem stated in Equation (8) is readily transformed in an algebraic nonlinear optimization problem as the system inputs and outputs are available in closed form, where the need for integrating the underlying ODE system and the need for control and or state vector parameterization are dropped. As a last step, the resulting optimal control profiles and initial conditions are employed to conduct a new experiment that is expected to deliver additional informative data regarding model falsification.

After conducting the optimally designed experiment, gathering of corresponding data, and the model parameters update, the model assessment is performed. Criteria for model assessment and selection are subject to different dimensions: falsifiability, explanatory adequacy, interpretability, faithfulness, goodness-of-fit, complexity/simplicity and generalizability [35]. For instance, a well-known and commonly used criterion is the Akaike information criterion ($\mathrm{AIC}$) or its small-sample analog, the corrected Akaike information criterion ${\mathrm{AIC}}_{\mathrm{c}}$ [67], which evaluates the model response residuals, as well as the model complexity in terms of the dimension of the model parameter vector. In general, however, there is no master strategy for model discrimination and therefore, no universal criterion. The method applied depends on the stage of modeling, the type of uncertainty expected, and the strategy chosen [68]. For instance, the so-called overlap approach emphasizes the fact that in reality, model parameters are uncertain, which is in turn not regarded in many common discrimination criteria as the $\mathrm{AIC}$ or the Bayesian information criterion $\mathrm{BIC}$ [31]. Assuming additive measurement noise and according to the Doob–Dynkin lemma [69], the identified model parameters $\widehat{\mathsf{\theta }}$ can be considered to be random variables, where the probability space ($\Omega ,\mathcal{F},P$) contains the sample space $\Omega$, the $\sigma$-algebra $\mathcal{F}$, and the probability measure P. Thus, the overlap approach accounts for parameter uncertainties by considering the trajectory confidence intervals, which is in contrast to the $\mathrm{AIC}$, where only the point estimate at the maximum likelihood and the corresponding model responses are used [68]. Adjusted to the requirements at hand, the overlap of two models considering parametric uncertainties is

$$OVL=\sum _{k=0}^T\frac{2\sqrt{{\sigma }_{1,k}^2{\sigma }_{2,k}^2}}{{\sigma }_{1,k}^2+{\sigma }_{2,k}^2}exp\frac{-0.5{({\mu }_{1,k}-{\mu }_{2,k})}^2}{{\sigma }_{1,k}^2+{\sigma }_{2,k}^2},$$

(10)

where ${\mu }_{i,k}$ and ${\sigma }_{i,k}$ are the expected value and the variance of the output function of model i at time point k resulting from the stochastic nature of the parameters, respectively. In the next subsection, we illustrate how these statistical quantities in Equation (10) can be quantified efficiently with the point estimate method (PEM).

3.1. Point Estimate Method

The mean and variance occurring in the overlap formula (Equation (10)) are known as the first raw and second central moments of a continuous random variable and are defined as

$$\begin{array}{cc}\hfill \mu & =\mathbb{E}\left[k\left(\mathsf{\theta }\right)\right]={\int }_{\Omega }k\left(\mathsf{\theta }\right)P\left(\mathsf{\theta }\right)\mathrm{d}\mathsf{\theta },\hfill \end{array}$$

(11a)

$$\begin{array}{cc}\hfill {\sigma }^2& =\mathbb{V}\left[k\left(\mathsf{\theta }\right)\right]={\int }_{\Omega }{(k\left(\mathsf{\theta }\right)-\mu )}^{2}P\left(\mathsf{\theta }\right)\mathrm{d}\mathsf{\theta },\hfill \end{array}$$

(11b)

where $P\left(\mathsf{\theta }\right)$ is the probability density function related to the random model parameter vector $\mathsf{\theta }$ of a process model $\mathcal{M}$ (Equations (1) and (2)). $k(·)$ represents a generic nonlinear function. Usually, no closed solutions of the integrals on the right side of Equations (11a) and (11b) are available, and must be evaluated numerically. For accurate results, Monte Carlo simulations can be performed that draw many samples from the probability density functions. Thus, in particular, in the case of complex functions $k\left(\mathsf{\theta }\right)$ and many model parameters, the evaluation of the integral becomes computationally expensive or even prohibitive due to the so-called curse of dimensionality. Alternatively, the PEM, initially developed for generic multi-dimensional integration problems over symmetrical regions, is a credible and practical method for uncertainty propagation with low computational cost; see [33] and references within. In process systems engineering, the PEM has been successfully applied to robustify various optimization problems, including non-symmetrical probability density functions, correlated model parameters, and imprecise parameter uncertainties [3,70,71].

The PEM approximates the integrals in Equation (11) by summing over $n_{PEM}$ weighted sampling points

$${\int }_{\Omega }k\left(\mathsf{\theta }\right)P\left(\mathsf{\theta }\right)\mathrm{d}\mathsf{\theta }\approx \sum _{l=1}^{n_{PEM}}{w}_{l}k\left({\mathsf{\theta }}_l\right),$$

(12)

with weight factors $w_l$ and parameter vector realizations ${\mathsf{\theta }}_l$. In detail, assuming a nominal parameter vector ${\mathsf{\theta }}_0$, dedicated model parameter vector realizations, ${\mathsf{\theta }}_l$, form a parameter vector set, ${\mathsf{\theta }}_l\in \mathcal{O}:=\{{\mathsf{\theta }}_0,{\mathcal{O}}_1,-{\mathcal{O}}_1,{\mathcal{O}}_2,-{\mathcal{O}}_2,{\mathcal{O}}_3,-{\mathcal{O}}_3\}$, where

$$\begin{array}{cc}\hfill & {\mathcal{O}}_1:=\{{\mathsf{\theta }}_0\left[i\right]+\vartheta ,\forall i\in \{1,\dots ,p\}\},\hfill \\ \hfill & {\mathcal{O}}_2:=\{{\mathsf{\theta }}_0\left[(i,j)\right]+[+\vartheta ,+\vartheta ],\forall i,\underset{j>i}{j}\in \{1,\dots ,p\}\},\hfill \\ \hfill & {\mathcal{O}}_3:=\{{\mathsf{\theta }}_0\left[(i,j)\right]+[-\vartheta ,+\vartheta ],\forall i,\underset{j>i}{j}\in \{1,\dots ,p\}\}.\hfill \end{array}$$

(13)

Here, ${\mathsf{\theta }}_0\left[i\right]$ means that the ith element of the nominal parameter vector, ${\mathsf{\theta }}_0$, is permuted via the spreading parameter, $\vartheta$, and ${\mathsf{\theta }}_0\left[(i,j)\right]$ that the ith and the jth elements are modified, respectively. Note that the weight factors, $w_l$, as well as the spreading parameter, $\vartheta$, are determined via a corresponding algebraic equation system; the interested reader is referred to [33] and references therein. Based on the parameter samples, ${\mathsf{\theta }}_l$, we can approximate the statistics of a given nonlinear function. For instance, the approximations of the expected value and the variance as stated in Equations (11a) and (11b) read as

$$\begin{array}{cc}\hfill \mu & \approx w_{0}k\left({\mathsf{\theta }}_0\right)+w_1\sum _l^{|{\mathcal{O}}_{w_1}|}k\left({\mathsf{\theta }}_l\right)+w_2\sum _l^{|{\mathcal{O}}_{w_2}|}k\left({\mathsf{\theta }}_l\right),\hfill \end{array}$$

(14a)

$$\begin{array}{cc}\hfill {\sigma }^2& \approx w_0{(k\left({\mathsf{\theta }}_0\right)-\mu )}^2+w_1\sum _l^{|{\mathcal{O}}_{w_1}|}{(k\left({\mathsf{\theta }}_l\right)-\mu )}^2+w_2\sum _l^{|{\mathcal{O}}_{w_2}|}{(k\left({\mathsf{\theta }}_l\right)-\mu )}^2,\hfill \end{array}$$

(14b)

with ${\mathcal{O}}_{w_1}:=\{{\mathcal{O}}_1,-{\mathcal{O}}_1\}$ and ${\mathcal{O}}_{w_2}:=\{{\mathcal{O}}_2,-{\mathcal{O}}_2,{\mathcal{O}}_3,-{\mathcal{O}}_3\}$. Please note that the overall parameter sample number, $n_{PEM}$, used in Equations (14a) and (14b) scales quadratically to the dimension of uncertain model parameters:

$$n_{PEM}=2p^2+1.$$

(15)

At the cost of accuracy (i.e., approximation error introduced via Equation (13)), the required sample numbers in Equations (14a) and (14b) can be reduced according to:

$$\begin{array}{cc}\hfill \mu & \approx w_{0}k\left({\mathsf{\theta }}_0\right)+w_1\sum _l^{|{\mathcal{O}}_{w_1}|}k\left({\mathsf{\theta }}_l\right),\hfill \end{array}$$

(16a)

$$\begin{array}{cc}\hfill {\sigma }^2& \approx w_0{(k\left({\mathsf{\theta }}_0\right)-\mu )}^2+w_1\sum _l^{|{\mathcal{O}}_{w_1}|}{(k\left({\mathsf{\theta }}_l\right)-\mu )}^2,\hfill \end{array}$$

(16b)

where the overall parameter sample number, $n_{PEM}$, scales linearly to the dimension of uncertain model parameters:

$$n_{PEM}=2p+1.$$

(17)

Equations (14) and (16) are exact for monomials of up to degrees five and three [33]. Thus, in what follows, we use the terms PEM5 and PEM3 to distinguish between the two PEM approximation schemes. In

Table 1. Weights and spread parameter of the PEM for a standard Gaussian distribution [33].

Approximation Scheme	${\mathit{w}}_0$	${\mathit{w}}_1$	${\mathit{w}}_2$	$\mathit{\vartheta }$
PEM5 (Equations (14a) and (14b))	$1+\frac{p^2-p}{18}$	$\frac{4-p}{18}$	$\frac{1}{36}$	$\sqrt{3}$
PEM3 (Equations (16a) and (16b))	$1-\frac{p}{{\vartheta }^2}$	$\frac{1}{2{\vartheta }^2}$	-	$\vartheta$

, we summarize the values for the weights, $w_l$, and the spread parameter, $\vartheta$, assuming a standard Gaussian distribution. In the case of PEM3, the spread parameter, $\vartheta$, can be considered to be a design parameter, which is frequently set to the corresponding PEM5 value; that is, $\vartheta =\sqrt{3}$.

Please note that any parametric or non-parametric probability distribution of relevant model parameters can be considered via a (non)linear transformation step, including parameter correlations [70].

3.2. Robust Flatness-Based Objective Function

For the robust formulation of the nominal optimization problem (Equation (9)), we aim at introducing an additional term in the objective function that penalizes the propagated uncertainty that is quantified by the variances $\mathbb{V}$ of the models’ inputs $\mathbf{u}$. In this vein, we expect the calculated overlap (Equation (10)) to decrease as the confidence bands of the trajectories tighten. The objective function is balanced by a weight factor $\alpha \in [0,1]$ leading to Pareto optimal points for its different values. The optimization problem for a flatness-based robust design of experiments for model selection is then stated as

$$\begin{array}{c}\hfill \begin{array}{ccc}\hfill & \underset{{\mathsf{\theta }}^{*}}{\max }\hfill & \sum _{k=0}^T(\alpha \sum _{i=1}^{M-1}\sum _{j=i+1}^M{\left[{\mathbf{y}}^i({\mathsf{\theta }}_i^{*},t_k)-{\mathbf{y}}^j({\mathsf{\theta }}_j^{*},t_k)\right]}^2\hfill \\ \hfill & \hfill & (1-\alpha )\sum _{i=1}^M\mathbb{V}\left[{\mathbf{u}}_i({\mathsf{\theta }}_i^{*},t_k)\right])\hfill \\ \hfill & \mathrm{s}.\mathrm{t}.\hfill & \Delta \mathbf{u}({\mathsf{\theta }}^{*},t_k)=0,\hfill \\ \hfill & \hfill & \Delta {\mathbf{x}}_0({\mathsf{\theta }}^{*},t_k)=0,\hfill \\ \hfill & \hfill & {\mathbf{g}}_{\mathrm{eq}}({\mathsf{\theta }}^{*},t_k)=0,\hfill \\ \hfill & \hfill & {\mathbf{g}}_{\mathrm{ineq}}({\mathsf{\theta }}^{*},t_k)\le 0.\hfill \end{array}\end{array}$$

(18)

The decision variables of the optimization problem are the meta-parameters present in the basis functions for the flat outputs, i.e., the control vectors of the B-spline curves. As previously in the nominal optimization problem (Equation (9)), the robust formulation is a programming problem for which all occurring functions are algebraically available. Thus, the nominal version (Equation (9)) and the robust version (Equation (18)) are readily stated as an NLP without any use of parameterization methods. Furthermore, the possibility of deriving functions for exact gradients is given as opposed to the application of automatic differentiation or finite differences methods. A flow chart for the proposed robust MBDoE strategy using differential flatness is given in Figure 3.

3.3. Computational Methods

The whole program is coded in Julia [72]. The derivation of flat systems is based on the comprehensive Julia suite on differential equations [73]. For symbolic calculation, the Python package SymPy is used where its functions can be accessed from Julia via an interface. As modeling language for the optimization problems, Julia-based JuMP was chosen [74]. The nonlinear programming problems were solved with the open-source large-scale interior point solver IPOPT [75]. All calculations were run on a standard desktop machine and finished within an hour, where the bulk of the computational load was spent in symbolic calculations and the quantification of the propagated uncertainty.

4. Case Study

Biocatalytic reactions appear in numerous syntheses of natural products and active pharmaceutical ingredients. For example, chiral hydroxy ketones are important building blocks in the pharmaceutical industry and can be produced from aldehydes via enzymatic carboligation [22]. A lack of mechanistic kinetic models for biocatalytic carboligation and simultaneous inactivation of the benzaldehyde lyase was studied in [76] using MBDoE, where the model quality of the different model candidates was analyzed using the $\mathrm{AIC}$. Motivated by the biocatalytic carboligation of ketones from aldehydes, in this simulation study, we consider two models simulating the enzymatic step from benzaldehyde to benzoin [22]. In the enzymatic reaction chain, the forward reaction forming the intermediate is significantly larger than its reverse counterpart leading to the two model candidates in

Table 2. Reaction scheme of model ${\mathcal{M}}_1$.

Model ${\mathcal{M}}_1$
$2{\mathrm{S}}_1+\mathrm{E}\stackrel{k1}{\to }{\mathrm{C}}_1\stackrel{k2}{\to }\mathrm{P}+\mathrm{E}$
${\mathrm{S}}_2+\mathrm{E}\stackrel{k3}{\to }{\mathrm{C}}_2$

and

Table 3. Reaction scheme of model ${\mathcal{M}}_2$.

Model ${\mathcal{M}}_2$
$2{\mathrm{S}}_1+\mathrm{E}\stackrel{k1}{\to }{\mathrm{C}}_1\stackrel{k2}{\to }\mathrm{P}+\mathrm{E}$

.

The two model candidates differ in that way that model candidate ${\mathcal{M}}_1$ suffers from a loss reaction where a second substrate ${\mathrm{S}}_2$ is irreversibly inhibiting. Additionally, substrate ${\mathrm{S}}_2$ is expected to be constant via proper control actions over the experimental run. The corresponding ODE systems of the two model candidates are

$$\begin{array}{c}\hfill {\mathcal{M}}_1=\left\{\begin{array}{cc}{\displaystyle \frac{\mathrm{d}{x}_1}{\mathrm{d}t}}\hfill & =-2k_1{x}_1^2{x}_3+u_1\hfill \\ {\displaystyle \frac{\mathrm{d}{x}_2}{\mathrm{d}t}}\hfill & =k_1{x}_1^2{x}_3-k_2{x}_2\hfill \\ {\displaystyle \frac{\mathrm{d}{x}_3}{\mathrm{d}t}}\hfill & =-k_1{x}_1^2{x}_3+k_2{x}_2-k_3\left[{\mathrm{S}}_2\right]x_3+u_2\hfill \\ {\displaystyle \frac{\mathrm{d}{x}_4}{\mathrm{d}t}}\hfill & =k_2{x}_2\hfill \end{array}\right.,\end{array}$$

(19)

$$\begin{array}{c}\hfill {\mathcal{M}}_2=\left\{\begin{array}{cc}{\displaystyle \frac{\mathrm{d}{x}_1}{\mathrm{d}t}}\hfill & =-2k_1{x}_1^2{x}_3+u_1\hfill \\ {\displaystyle \frac{\mathrm{d}{x}_2}{\mathrm{d}t}}\hfill & =k_1{x}_1^2{x}_3-k_2{x}_2\hfill \\ {\displaystyle \frac{\mathrm{d}{x}_3}{\mathrm{d}t}}\hfill & =-k_1{x}_1^2{x}_3+k_2{x}_2+u_2\hfill \\ {\displaystyle \frac{\mathrm{d}{x}_4}{\mathrm{d}t}}\hfill & =k_2{x}_2\hfill \end{array}\right.,\phantom{xxxx}\end{array}$$

(20)

where the indices are assigned according to $\{1,2,3,4\}=\{{\mathrm{S}}_1,{\mathrm{C}}_1,\mathrm{E},\mathrm{P}\}$ and substrate ${\mathrm{S}}_1$ and enzyme E can be dosed over the course of the experiment. Please note that in model ${\mathcal{M}}_1$, the differential equations of the second substrate ${\mathrm{S}}_2$, which is assumed to be constant, and the loss product ${\mathrm{C}}_2$ are not specified, as information about their time behavior does not influence the other differential equations, and therefore, is not relevant to the problem at hand. The measured concentrations are the concentrations of the substrate and the final product; that is, $\mathbf{y}={(x_1,x_4)}^{\top }$. The kinetic constants have the following expected values: For model ${\mathcal{M}}_1$ $(k_1,k_2,k_3)=(0.215{\mathrm{L}}^2/\left({\mathrm{mmol}}^2\min \right),91.475{\min }^{-1},0.019\mathrm{L}/\left(\mathrm{mmol}\min \right))$, and for model ${\mathcal{M}}_2$ $(k_1,k_2)=(0.204{\mathrm{L}}^2/\left({\mathrm{mmol}}^2\min \right),103.344{\min }^{-1})$. The variances of the parameters equal $5\%$ of their expected values in the absence of parameter correlations.

5. Results

For the case study given in Section 4, we analyzed the differential flatness property, and then, we designed discriminating input controls based on the proposed MBDoE approach, which makes use of the differential flatness concept. The effect of model parameter uncertainties was studied as well, and the robust MBDoE approach was applied as an appropriate countermeasure.

5.1. Differential Flatness Property

We derive the flat outputs for both models following heuristic methods to (1) obtain a flat output candidate according to Equation (3) and (2), with the help of graph theory, to show that the candidate fulfills the differential flatness conditions given Equations (4)–(6). Drawing the directed graphs (digraphs) for model ${\mathcal{M}}_1$ and model ${\mathcal{M}}_2$, we observe that they look alike. The corresponding adjacency matrices are

$$A_u=\begin{array}{c}\\ x_1\\ x_2\\ x_3\\ x_4\end{array}\begin{array}{c}\begin{array}{cc}{u}_1& u_2\end{array}\\ \left(\begin{array}{cc}1& 0\\ 0& 0\\ 0& 1\\ 0& 0\end{array}\right)\end{array},A_x=\begin{array}{c}\\ x_1\\ x_2\\ x_3\\ x_4\end{array}\begin{array}{c}\begin{array}{cccc}{x}_1& x_2& x_3& x_4\end{array}\\ \left(\begin{array}{cccc}1& 0& 1& 0\\ 1& 1& 1& 0\\ 1& 1& 1& 0\\ 0& 1& 0& 0\end{array}\right)\end{array},A_y=\begin{array}{c}{y}_1\\ y_2\end{array}\begin{array}{c}\begin{array}{cccc}{x}_1& x_2& x_3& x_4\end{array}\\ \left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 0& 0& 1\end{array}\right)\end{array}$$

and the resulting digraph is shown in Figure 4a. The digraph is composed of 8 vertices, $V=\{u_1,u_2\}\cup \{x_1,x_2,x_3,x_4\}\cup \{y_1,y_2\}$, and 13 edges corresponding to the non-zero entries in the adjacency matrices $A_u$, $A_x$, and $A_y$. The self-loops of $\{x_1,x_2,x_3\}\in V$ are related to the non-zero diagonal elements of $A_x$.

The control vector dimension in both models is 2. Therefore, to comply with condition given in Equation (6), the dimension of the flat output vector must be 2 as well. Inspecting the system of differential equations in Equation (19) one might observe that $x_4$ does not appear on the right side of any of the differential equations. The corresponding node in the digraph (see Figure 4a) represents a dead end (i.e., no edge from $x_4\in V$ to $\{x_1,x_2,x_3\}\in V$), and therefore, must be part of the flat output. Obviously, the second element of the flat output vector should be chosen from the concentrations left as they have direct physical meaning for the study case at hand. We consider the flat output candidate

$$\begin{array}{c}\hfill {\mathbf{y}}^{\mathrm{flat}}=\left(\begin{array}{c}{\gamma }_1\\ {\gamma }_2\end{array}\right)=\left(\begin{array}{c}{x}_1\\ x_4\end{array}\right)\end{array}$$

(21)

that evidently satisfies Equation (3).

The explicit expressibility graph, which can be readily obtained from the digraph, is shown in Figure 4b. Please note that in comparison to the digraph (Figure 4a), the self-loops are omitted, and the edge directions are reversed. In Figure 4b, two disjoint paths are drawn that cover all vertices related to the inputs, system states, and flat outputs. Therefore, the flat output vector (Equation (21)) is supposed to form a differentially flat system which we will confirm using its definition in the following. From the ODE system, after several reformulations and substitutions, we obtain the inverse models:

$$\begin{array}{c}\hfill {\mathcal{M}}_1^{-1}=\left\{\begin{array}{cc}{\displaystyle x_1}\hfill & ={\gamma }_1\hfill \\ {\displaystyle x_2}\hfill & =\frac{{\dot{\gamma }}_2}{k_2}\hfill \\ {\displaystyle x_3}\hfill & =\frac{\frac{{\ddot{\gamma }}_2}{k_2}-{\dot{\gamma }}_2}{k_1{\gamma }_1^2}\hfill \\ {\displaystyle x_4}\hfill & ={\gamma }_2\hfill \\ {\displaystyle u_1}\hfill & ={\dot{\gamma }}_1+2{\gamma }_1^2\frac{\frac{{\ddot{\gamma }}_2}{k_2}-{\dot{\gamma }}_2}{{\gamma }_1^2}\hfill \\ {\displaystyle u_3}\hfill & ={\mathcal{F}}_1({\gamma }_1,{\dot{\gamma }}_1,{\gamma }_2,{\dot{\gamma }}_2,{\ddot{\gamma }}_2,{\stackrel{⃛}{\gamma }}_2)\hfill \end{array}\right.\end{array}$$

(22)

and

$$\begin{array}{c}\hfill {\mathcal{M}}_2^{-1}=\left\{\begin{array}{cc}{\displaystyle x_1}\hfill & ={\gamma }_1\hfill \\ {\displaystyle x_2}\hfill & =\frac{{\dot{\gamma }}_2}{k_2}\hfill \\ {\displaystyle x_3}\hfill & =\frac{\frac{{\ddot{\gamma }}_2}{k_2}-{\dot{\gamma }}_2}{k_1{\gamma }_1^2}\hfill \\ {\displaystyle x_4}\hfill & ={\gamma }_2\hfill \\ {\displaystyle u_1}\hfill & ={\dot{\gamma }}_1+2{\gamma }_1^2\frac{\frac{{\ddot{\gamma }}_2}{k_2}-{\dot{\gamma }}_2}{{\gamma }_1^2}\hfill \\ {\displaystyle u_3}\hfill & ={\mathcal{F}}_2({\gamma }_1,{\dot{\gamma }}_1,{\gamma }_2,{\dot{\gamma }}_2,{\ddot{\gamma }}_2,{\stackrel{⃛}{\gamma }}_2)\hfill \end{array}\right..\phantom{and}\end{array}$$

(23)

Please note that we omitted to write out the equations for the second element of the output vector $u_2$ for model ${\mathcal{M}}_1$ and model ${\mathcal{M}}_2$ due to readability reasons. However, keep in mind that the functions ${\mathcal{F}}_1$ and ${\mathcal{F}}_2$ in Equations (22) and (23) differ. The inverse model (Equation (22)) with its flat output (Equation (21)) satisfies the set containing all three conditions for differential flatness (Equations (4)–(6)). Accordingly, the inverse model of model ${\mathcal{M}}_2$ (Equation (23)) fulfills the conditions for differential flatness. Thus, both models are differentially flat models. For prescribed trajectories ${\gamma }_1^{*}=x_1^{*}$ and ${\gamma }_2^{*}=x_4^{*}$, the necessary experimental conditions for the controls are readily available by differentiation; for an example, see [20].

5.2. Non-Optimized Experimental Design

We first consider the non-optimized experimental design, i.e., the initial situation before any optimization iteration. Please note that a batch process is assumed, and thus, the controls are zero and not explicitly shown. Plots for the measured states and the measured states with confidence bands are shown in Figure 5 on the left side and on the right side, respectively.

For the nominal case in Figure 5a, it is impossible to distinguish the two models from each other, even if the experimental data were close to one of the models. In Figure 5b, it is well observable that the uncertainties in the measured states given as confidence intervals are low. To better compare the results with the expected results from the optimization part, we normalize the overlap measure (Equation (10)) to the overlap corresponding to Figure 5b; that is, $OV{L}_{\mathrm{N}}=1.0$.

5.3. Optimized Experimental Design without Uncertainty: The Nominal Case

Before optimizing, suitable B-spline types must be chosen. In contrast to the classical application of splines in approximation or interpolation of available data, they are used to express smooth control and state trajectories along the time axis, and are created from their definition as opposed to a fitting process. Thus, common problems experienced in spline approximation and interpolation, like cusps and loops [77], are less likely to occur. From the different methods for splines, we choose the uniform method where the data points are uniformly distributed over the domain range. Furthermore, for each element of the flat output vectors, a B-spline curve of order 6 with 12 control points and clamped end conditions was used, resulting in a 48-dimensional decision space for the optimization problem. In the first optimization step, we follow the classical deterministic optimization by setting $\alpha =1$ in the setting (Equation (18)). The results for the measured states of the nominal case are shown in Figure 6, again, without uncertainty and with uncertainty vis-à-vis. The corresponding control curves are displayed in Figure 7 on the left side.

In the plot without propagated uncertainties, it is clearly visible that the expected values of the measured states are significantly driven apart. However, a look at the plot on the right reveals that (1) the uncertainties have substantially enlarged compared to the initial situation, and (2) the confidence bands of both models overlap noticeably. The overlap values are given in

Table 4. Results for the overlap.

Design	None	Nominal	PEM3	PEM5
$OV{L}_{\mathrm{N}}$	$1.0$	$0.94$	$0.81$	$0.80$

. The overlap has decreased versus the initial situation, i.e., it is 6 percentage points lower compared to before the optimization. However, we also came across simulation runs where the overlap worsened in the sense of a larger overlap value if the nominal optimization state was compared with the non-optimized. In this context, the question arises if we can perform better, i.e., further lower the overlap by inherently taking care of the parametric uncertainties in the optimization setting. Therefore, we robustify the optimization by considering the variances in the objective function as described in Section 3.1.

5.4. Robust Experimental Design

For the robust MBDoE setting (Equation (18)), we set $\alpha =0.5$; that is, the Euclidean distance of the flat outputs of model ${\mathcal{M}}_1$ and ${\mathcal{M}}_2$ as well as the resulting uncertainties in the recalculated model inputs are equally weighted. Compared with the nominal optimization, the computational time to solve reaching satisfying convergence increases multiple-fold despite the low number of sampling points that the PEM is using. In particular, when PEM5 is employed, the computational time increases considerably due to the higher sample number; see Equations (15) and (17). The concentration curves for the resulting states of the robust flatness-based design are shown in Figure 8.

The trajectories for the expected values in the left subfigure clearly drift apart, but less strongly than in the nominal optimization case; see Figure 6. Considering 100 equally spaced points, the norms of the model output distances in the robust case are roughly half as large as in the nominal case. The presumption that the overlap has decreased can be confirmed when looking at the overlap values in Table 4.

The overlap value has decreased by 19 percentage points and by 13 percentage points when compared with the initial situation and the nominal optimization, respectively. The propagated uncertainties can be extracted from the graphs shown in Figure 9 where the variances for both outputs of model 1 are considerably lower when the nominal case on the left is compared with the robust one on the right. Furthermore, the approximation of the propagated uncertainty using the PEM is sufficiently good, as the results are in proximity to the Monte Carlo results with ${10}^4$ samples. Thus, the inclusion of the propagated uncertainty as a penalty in the objective function has delivered the desired outcome. Prospectively, a different choice of $\alpha$ or the drawing of a Pareto front might further decrease the overlap between the two model outcomes. The corresponding optimal controls for the substrate and enzyme are drawn in Figure 7 on the right side. Please note that the realization of the profiles in a laboratory setup is constrained by the available pump devices. These constraints could be included in the robust MBDoE framework (Equation (18)), but the technical realization of these input controls is beyond the scope of this study.

6. Conclusions

We have extended a novel approach to the model-based design of experiments for model selection using the differential flatness property of the systems by robustifying the optimization problem against parametric uncertainty. In comparison with commonly used methods for model discrimination, the resulting optimal control problem does not require approximation methods for solving the underlying system of differential equations as both controls and states can be derived analytically in differentially flat systems without numerical integration. Likewise, parameterization of the controls and potentially the states is obsolete for we directly arrive at a nonlinear programming problem. Moreover, the possibility to provide analytic gradients to the optimization solver, even for highly nonlinear systems, is given. The robustification strategy comprised the consideration of propagated uncertainties; that is, the variances of the reconstructed model inputs were added as a penalty term in the objective function. The corresponding uncertainty quantification was performed using the point estimate method, a computationally cheap but reasonably accurate method compared with standard Monte Carlo simulations for uncertainty quantification. We successfully applied the strategy to a nonlinear enzymatic reaction using B-splines as the parameterization technique for expressing the flat outputs, and thus, obtained optimal experimental input controls for the subsequent experiment. We showed that a nominal optimization without the consideration of parameter uncertainties might result in unreliable optimal controls, i.e., not leading to the possibility of selecting one model over another after having performed the subsequent experiment. In contrast, if parametric uncertainties were inherently considered in the optimal experimental design problem, more reliable MBDoE results were identified. We emphasize the importance of including the quantification of uncertainties over classic deterministic optimization where the differential flatness concept shows promising characteristics for advanced system identification strategies in terms of precise parameter estimates and reliable model selection.