Distributionally robust parameter estimation for nonlinear fed-batch switched time-delay system with moment constraints of uncertain measured output data

1.   Introduction

Distributionally robust optimization (DRO) is an optimization method designed to address uncertainty, particularly in situations where the underlying probability distribution is not fully known. Unlike traditional robust optimization, DRO assumes that we only have partial information about the uncertain probability distribution, rather than complete certainty. By accounting for the worst-case scenarios, DRO can reduce decision-making risks and enhance the robustness of systems across various real-world applications. It has aroused much interest among researchers to develop algorithms to solve the DRO problems in kernel density estimation ^[1], data-driven problem ^[2,3], portfolio optimization problem ^[4], decision bounding problem ^[5], and matrix moment constraint ^[6]. However, no dynamical systems are involved in the aforementioned DRO problems. Recently, the DRO problems involving linear discrete-time systems are discussed in multi-access space-time block coded MIMO system ^[7], stochastic model predictive control problem ^[8], modeling of optimal control problem ^[9], constrained stochastic system ^[10], partially observable linear stochastic system ^[11], and output feedback ^[12]. However, in these literatures, the dynamical systems involved in the DRO problems are discrete-time dynamical (DTD) systems. Continuous-time dynamical (CTD) systems are not involved. There are major difficulties in solving DRO problems governed by CTD systems such as (i) a DRO problem constrained by a CTD system is fundamentally a bi-level optimization problem, which, unlike a DTD system, cannot be simplified to a single-level optimization problem with a finite number of constraints; and (ii) the presence of CTD systems in DRO problems causes the inner-level objective functional to be expressed as the expectation of a non-convex function. In this paper, we examine a DRO problem influenced by a nonlinear continuous-time switched time-delay (NCTSTD) system encountered in fed-batch production of 1, 3-propanediol (1, 3-PD). That is to say, we directly address the DRO problem governed by a CTD system, rather than discretizing the CTD system.

1, 3-PD is a key chemical raw material for polymer production because of its unique symmetric structure ^[13]. The methods for the production of 1, 3-PD can be categorized into two types: microbial conversion and chemical synthesis ^[14]. This paper focuses on the first category due to its eco-friendliness ^[15]. 1, 3-PD microbial conversion process of glycerol can be divided into three categories ^[16]: (i) batch culture ^[17], (ii) continuous culture ^[18,19], and (iii) fed-batch culture (see Figure 1 in ^[20]). In fed-batch culture, substrate inhibition is dramatically decreased, allowing for greater consumption of glycerol and alkali, which in turn leads to increased biomass production and higher concentrations of 1, 3-PD. For fed-batch culture, many interesting papers appeared in the existing literature such as those involving state-dependent impulsive system ^[21], state-dependent switched system ^[22], multistage optimal control ^[23], optimal control ^[24], a novel downstream process ^[25], and Koopman modeling ^[26,27]. However, the references mentioned do not consider the presence of time delays.

Time delays widely exist in many practical problems, including optimal control problem ^[28], predator-prey model ^[29], perimeter control problem ^[30], the survey ^[31], sequential time scaling transformation ^[32], non-autonomous time-delayed SIR model ^[33], multiregion urban traffic network ^[34], and sparse optimal control ^[35]. At the beginning of the fed-batch process, all substrates are placed into the fermentator. However, the chemical reaction does not trigger immediately ^[36], but needs to go through the process in which all substrates need to be well-mixed. Clearly, completing this process will require a certain amount of time. Thus, this process should be characterized by time-delay systems. In the context of fed-batch culture with time-delay, numerous intriguing papers have been published in the existing literature, such as robust parameter identification ^[37], multi-objective optimization ^[38], and optimal control ^[39]. However, the aforementioned papers assumed that the probability distribution of the output measurement data is fully available. This assumption is obviously not realistic. In practical applications, typically only first-order moment information is obtainable. This is the main motivation for this paper.

This paper examines an NCTSTD system with unknown time-delay and system parameters to characterize the bioconversion process of 1, 3-PD in fed-batch culture of glycerol caused by K. pneumoniae. Furthermore, the uncertain measurement output data is regarded as a stochastic variable and the only first-order moment information of its distribution is obtainable. Our aim is to identify these unknown quantities in a setting with uncertain output data measurement. Considering these factors, we propose a distributionally robust parameter estimation (DRPE) problem governed by the NCTSTD system, where the objective function contains two levels: (i) the inner-level is to maximize the expectation of the discrepancy between the output of the NCTSTD system and the uncertain output data measurement, and (ii) the outer-level is to minimize with respect to the unknown quantities. By the duality theory, the DRPE problem is equivalently translated into a single-level parameter estimation (SLPE) problem with non-smooth term smoothed by a smoothing function. Its convergence analysis is also given. Moreover, the gradients of the objective function are derived. A hybrid algorithm is proposed for solving the SLPE problem. Finally, numerical results show the applicability of the NCTSTD system to characterize the process of fed-batch culture and the availability of the proposed hybrid algorithm for solving the SLPE problem.

The contributions of the paper are threefold:

1) Different from the above literature, the optimal solution of the DRPE problem (OSDRPEP) is between the optimal solution of the robust parameter estimation problem (OSRPEP) and the optimal solution of the stochastic parameter estimation problem (OSSPEP). OSDRPEP builds a bridge between OSRPEP and OSSPEP, neither as conservative as OSRPEP nor as strict as OSSPEP when dealing with random probability distribution.

2) The DRPE problem subject to the CTD system is actually a bi-level parameter estimation problem, which cannot be described as an SLPE problem with finitely many equality constraints as the example of the DTD system. In this paper, a method on the basis of the duality theory is proposed for solving the DRPE problem governed by the CTD system.

3) The non-smooth term in the objective function is approximated by a smoothing technique. Some noteworthy theorems are required to prove that the smoothing function can overcome the problem that the constraint qualification is not satisfied.

The remainder of the paper is structured as follows. A DRPE problem is proposed in Section 2. A computational method is proposed for solving the DRPE problem in Section 3. Numerical results are given in Section 4. In Section 5, we draw some concluding remarks and imply some research directions in the future.

2.   Problem formulation

2.1.   NCTSTD system

In fed-batch culture, we assume that the conditions given below are satisfied.

$\bullet$ The concentrations of reactants are uniformly distributed in the reactor. Nonuniform space distribution is ignored.

$\bullet$ In fed-batch culture, the reactor is fed exclusively with glycerol and alkali.

Based on ^[43], the NCTSTD system governing the fed-batch culture can be described as follows:

For $t\in(\tau^{\ell-1}, \tau^{\ell}], \ell = 2i+1, i\in\Lambda_2, {f}^\ell(x(t), x(t-h), p)$ is of the form given by

where its components are

On the other hand, for $t\in(\tau^{\ell-1}, \tau^{\ell}], \ell = 2i+2, i\in\Lambda_1, {f}^\ell(x(t), x(t-h), p)$ is of the form given by

where its components are

with

The concentrations of substrate and the product of fermentation fluid are restricted to within a set $W$ defined by

Based on sensitivity analysis in ^[43], the time-delay $h$ and the kinetic parameter vector $p$, the meanings of which are listed in Table 1, are required to be identified. Such analysis allows us to focus our optimization efforts on these parameters and reduce unnecessary complexity. Thus, they are treated as decision variables to be chosen from the sets defined by

This prior information helps narrow down the parameter space, making the optimization process more stable.

2.2.   Distributionally robust parameter estimation problem

The aim of this section is to present a DRPE problem subject to the NCTSTD system with the unknown time-delay $h$ and system parameter vector $p$ to be optimized.

Let $x(\cdot|h, {p}): = (x_1(\cdot|h, {p}), ..., x_5(\cdot|h, {p}))^\mathrm{T}$ be the solution of system (2.1) (i.e., the system output) corresponding to each pair $(h, p)\in [0, \bar{h}]\times \mathcal{P}$. In ^{[20,37,38,39]}, assume that the distributions of output data are known precisely is an idealistic assumption. In fact, the distributions of the measured output data are often not available. Thus, the measured output data vector $\tilde{d}$ is regarded as a random vector. Its first-order moment is defined by

where $\mathbf{E}_{\mathbb{P}}$ represents the expected value according to the probability distribution $\mathbb{P}$; $\tilde{d}: = ((\tilde{d}^1)^\mathrm{T}, \cdots, (\tilde{d}^\sigma)^\mathrm{T})^\mathrm{T} = (\tilde{d}^1_1, \cdots, \tilde{d}^1_3, \cdots, \tilde{d}^\sigma_1, \cdots, \tilde{d}^\sigma_3)^\mathrm{T}\in \mathbb{R}^{3\sigma}; \; \tilde{d}^\upsilon$ refers to the output vector recorded at specific sample times $t_\upsilon, \upsilon\in I_\sigma$ subject to $0\leqslant t_1 < \cdots < t_\sigma\leqslant t_f;$ and $\tilde{z}: = (\tilde{z}^1_1, \cdots, \tilde{z}^1_3, \cdots, \tilde{z}^\sigma_1, \cdots, \tilde{z}^\sigma_3)^\mathrm{T}\in \mathbb{R}^{3\sigma}.$

The support set $\widetilde{\mathcal{D}}$ of $\tilde{d}$ is defined as

where $\underline{d}_\iota^\upsilon$ and $\overline{d}_\iota^\upsilon$ are, respectively, the lower and the upper limits of $\tilde{d}_\iota^\upsilon$ with $\underline{d}_\iota^\upsilon < \overline{d}_\iota^\upsilon$.

The ambiguity set $\widetilde{\mathfrak{P}}$ of probability measures of $\tilde{d}$ is defined as

Definition 1. For a given $(h, p)\in [0, \bar{h}]\times \mathcal{P}$, the discrepancy between the system output $x_\iota(t_\upsilon| h, p), t_\upsilon\in \mathcal{T}, \iota\in I_3, \upsilon\in I_\sigma,$ and the measured output data $\tilde{d}_\iota^\upsilon, \iota\in I_3, \upsilon\in I_\sigma,$ is defined by

where $\mathcal{T}$ represents the set of specified time points in fed-batch culture, and $|\mathcal{T}| = \sigma$ represents the number of elements in the set $\mathcal{T}.$

It is assumed that $\mathbb{P}\in\widetilde{\mathfrak{P}}$. Therefore, the worst-case scenario for $R_{e}(h, p, \tilde{d})$ is defined as $\max_{\mathbb{P}\in \widetilde{\mathfrak{P}}}\mathbf{E}_{\mathbb{P}}[R_{e}(h, {p}, \tilde{d})].$ Considering this, our DRPE problem is stated below:

where $g_j(x(t| h, p)) = -x_j(t| h, p), g_{j+5}(x(t| h, p)) = x_j(t| h, p)-\bar{x}_{j}, j\in I_5.$

3.   Solution approaches

This section presents a numerical approach for solving Problem A. Figure 1 summarizes the approximation and transformation process of solving Problem A. A thorough discussion is given in Section 3.1. Based on Theorems 3 and 4, the aim of Theorem 3 is to calculate the gradients of ${J}_{\rho, \epsilon}^{\mathbf{H}}(h, p, \psi, \zeta)$ in Problem H. Theorems 1 and 2 are applied to carry out the error analysis. Algorithm 1 is employed to solve Problem H.

3.1.   Single-level parameter estimation problem

Problem A is a bi-level parameter estimation problem subject to the NCTSTD system. In this subsection, Problem A is converted into a deterministic single-level parameter estimation problem.

The inner-level problem of Problem A can be expressed as given below.

For ease of use, the expectation of the measured output data vector $\tilde{d}$ is expressed using the Lebesgue integral. Then, Problem B is reformulated as:

Problem C is difficult to solve since it is a parameter estimation problem with respect to $\mathbb{P}\in \widetilde{\mathfrak{P}}$ on $\widetilde{\mathcal{D}}$. Considering this, in view of the duality theory ^[42], Problem C can be converted into its dual problem (i.e., Problem D), which can be solved easier than Problem C. Through the application of the duality theory ^[42], we obtain the duality of Problem $\mathbf{C}$ as follows:

where $\zeta: = (\zeta^1_1, \cdots, \zeta^1_3, \cdots, \zeta^\sigma_1, \cdots, \zeta^\sigma_3)^\mathrm{T}\in \mathbb{R}^{3\sigma}, \; \psi$ is the dual variable, and $R_e(h, {p}, \tilde{d})$ is as defined in (2.15).

Based on the definition of $R_e(h, {p}, \tilde{d})$ in (2.15) and $\max$, we have

Based on the separability of the variables $d_\iota^\upsilon, \iota\in I_3, \upsilon\in I_\sigma,$ constraint (3.2) is equivalent to constraint (3.3).

From the properties of quadratic functions, it follows that constraint (3.3) $\Leftrightarrow$ constraint (3.4).

Since $\max\Big\{a_1, a_2\Big\} = \frac{a_1+a_2}{2}+\frac{|a_1-a_2|}{2},$ constraint (3.4) is equivalent to constraint (3.5).

Letting $\alpha_\iota^\upsilon(x(t_\upsilon|h, p), \zeta): = (\frac{\underline{d}_\iota^\upsilon}{\overline{d}_\iota^\upsilon})^2+1-(2\frac{x_\iota(t_\upsilon| h, p)}{\overline{d}_\iota^\upsilon}+\zeta_\iota^\upsilon)(\frac{\underline{d}_\iota^\upsilon}{\overline{d}_\iota^\upsilon}+1)+2(\frac{x_\iota(t_\upsilon| h, p)}{\overline{d}_\iota^\upsilon})^2; \beta_\iota^\upsilon(x(t_\upsilon|h, p), \zeta): = (\frac{\underline{d}_\iota^\upsilon}{\overline{d}_\iota^\upsilon}-1)(\frac{\underline{d}_\iota^\upsilon}{\overline{d}_\iota^\upsilon}+1-2\frac{x_\iota(t_\upsilon| h, p)}{\overline{d}_\iota^\upsilon}-\zeta_\iota^\upsilon).$ Based on $\max\Big\{a_1, a_2\Big\} = \frac{a_1+a_2}{2}+\frac{|a_1-a_2|}{2},$ we have

where $\Lambda(h, p, \zeta, \psi): = \frac{1}{2}\sum_{\iota = 1}^{3}\sum_{\upsilon = 1 }^{\sigma}\Bigg\{\alpha_\iota^\upsilon(x(t_\upsilon|h, p), \zeta)+2\max\Big\{\beta_\iota^\upsilon(x(t_\upsilon|h, p), \zeta), 0\Big\}-\beta_\iota^\upsilon(x(t_\upsilon|h, p), \zeta)\Bigg\}-\psi.$

Problem A is equivalently transformed into the following SLPE problem:

3.2.   Constraint transcription technique

Problem E involves an infinite number of constraints (i.e., constraint (3.7) that restricts the state variables at each point in $[0, t_f]$). Therefore, Problem E can be regarded as a semi-infinite programming problem ^[46]. Based on the constraint transformation method involving the definition of $\max\{\cdot, \cdot\}$ ^[46], Problem E is equivalently transformed into Problem $\mathbf{E_c}.$

In view of the penalty function method ^[45,46], Problem $\mathbf{E_c}$ can be approximated by Problem $\mathbf{F}.$

Remark 1. It can be demonstrated that if the penalty parameter $\rho$ is sufficiently increased, then any minimizer of Problem $\mathbf{F}$ within the region $\zeta\in \mathbb{R}^{3\sigma}, \psi\in \mathbb{R}, h\in [0, \bar{h}], p\in \mathcal{P}$ will be feasible for Problem $\mathbf{E_c}$. Therefore, a solution to Problem $\mathbf{E_c}$ can be achieved by minimizing Problem $\mathbf{F}$ with appropriately chosen values for the penalty parameter $\rho$.

3.3.   Smoothing technique and error analysis

3.3.1.   Smoothing technique

One of the main hardships to solve Problem $\mathbf{F}$ numerically is that the second and the third terms in (3.9) are non-smooth, with respect to $h$ and $p$, due to the existence of the function $\max\{\cdot, 0\}$. To surmount this hardship, the function $\max\{\cdot, 0\}$ is approximated by the smoothing technique ^[44] as follows.

where $\eta$ is a continuous inequality constraint; $\rho$ is a penalty factor, $\kappa$ is the number of continuous inequality constraints, and $\epsilon > 0$ is called the smoothing parameter.

Approximate the non-smooth terms in Problem $\mathbf{F}$ using the function (3.10), and we get

where $\Lambda_\epsilon(h, p, \zeta, \psi): = \frac{1}{2} \sum\limits_{\iota=1}^3 \sum\limits_{\upsilon=1}^\sigma\Bigg\{\alpha_\iota^\upsilon(x(t_\upsilon|h, p), \zeta)+2\gamma_{\epsilon, \rho, 3\sigma}[\beta_\iota^\upsilon(x(t_\upsilon|h, p), \zeta)]-\beta_\iota^\upsilon(x(t_\upsilon|h, p), \zeta)\Bigg\}-\psi.$

3.3.2.   Error analysis

Since the functions involved with Problem H are smooth, it can be addressed using the gradient-based nonlinear programming (NLP) algorithm. It is evident that the difference between the solution of Problem H and the solution of Problem E is of importance. This section conducts an error analysis comparing Problem H, Problem F, and Problem E. With this in mind, Lemmas 1 and 2 are needed in the process of error analysis.

Lemma 1. For any $\epsilon > 0$, the function $\gamma_{\epsilon, \rho, \kappa}[\eta]$ has the following properties:

Proof. In view of Theorem 3.2.1 in ^[44], it is easy to prove the validity of (3.11) and (3.12) and

Now, we will prove (3.13). Based on (3.12), we have

which implies

According to (3.14) and (3.16), it follows that

which completes the proof. ■□

Lemma 2. For $\rho > 0$ and $\epsilon > 0,$ it holds that

Proof. Based on Lemma 1, we have

The proof is complete. ■□

Using Lemmas 1 and 2 as a foundation, we now proceed with the error analysis comparing Problem H, Problem F, and Problem E.

Theorem 1. Let $(h^*, p^*, \psi^*, \zeta^*)$ and $(h^{**}, p^{**}, \zeta^{**}, \psi^{**})$ be the solutions of Problems $\mathbf{F}$ and $\mathbf{H},$ respectively. Then,

Proof. In view of Lemma 2,

Because $(h^*, p^*, \psi^*, \zeta^*)$ is the solution of Problem $\mathbf{F},$ we have

which is equivalent to

As $(\tau^{**}, p^{**}, \theta^{**}, \eta^{**})$ is the solution of Problem $\mathbf{SLI_2},$ it is clear that

which is equivalent to

Thus, it follows from (3.19)–(3.22) that

The proof is complete. ■□

Definition 2. If $g_j(x(t|h^{\epsilon}, p^{\epsilon}))\leqslant \epsilon, j\in I_{10}; \Lambda(h^{\epsilon}, p^{\epsilon}, \psi^{\epsilon}, \zeta^{\epsilon})\leqslant \epsilon,$ then the set of decision variables $(h^{\epsilon}, p^{\epsilon}, \psi^{\epsilon}, \zeta^{\epsilon})$ is called $\epsilon-$feasible to Problem $\mathbf{E}.$

Theorem 2. Let $(h^*, p^*, \psi^*, \zeta^*)$ and $(h^{**}, p^{**}, \zeta^{**}, \psi^{**})$ be the solution of Problems $\mathbf{F}$ and $\mathbf{H},$ respectively. Furthermore, suppose that $(h^*, p^*, \psi^*, \zeta^*)$ is feasible to Problem $\mathbf{E}$ and that $(h^{**}, p^{**}, \zeta^{**}, \psi^{**})$ is $\epsilon-$feasible to Problem $\mathbf{E}.$ Then,

Proof. Due to the feasibility of $(h^{*}, p^{*}, \zeta^{*}, \psi^{*})$ in Problem $\mathbf{E}$, we get

From the $\epsilon-$feasibility of $(h^{**}, p^{**}, \zeta^{**}, \psi^{**})$ in Problem $\mathbf{E}$, it follows that

Based on (3.15) and (3.25), we obtain

From (3.24), (3.26), and the monotonicity of the $\gamma_{\epsilon, \rho, \kappa}[\eta]$ in (3.10), it follows that

which implies

This can be rearranged to yield

It is clear from Theorem 1 that

Combining (3.23) and (3.27) gives

The proof is complete. ■□

3.4.   State variation

In the process of the numerical optimization, the gradient formulae of ${J}_{\rho, \epsilon}^{\mathbf{H}}(h, p, \psi, \zeta)$ of Problem $\mathbf{H}$ with respect to $h, p, \psi,$ and $\zeta$ are required. With these gradient formulae, a gradient-based NLP algorithm can be employed for solving Problem $\mathbf{H}$. Therefore, this section is dedicated to deriving the gradient formulas.

3.4.1.   State variation with respect to time-delay $h$

The partial derivative of $x(t|\cdot, p)$ with respect to $h$ is known as the state variation with respect to $h.$

Lemma 3. For each pair $(h, p)\in [0, \bar{h}]\times \mathcal{P},$ it holds that

where $\varpi(t|h, p)\in \mathbb{R}^{5}$ is the solution of the following variational system:

Proof. This theorem can be proven in a manner similar to that used for Theorem 5 in ^[43]. ■□

3.4.2.   The state variation with respect to $p_\imath, \imath\in I_{9}$

To solve Problem $\mathbf{H}$, we also require the information of the state variation with respect to $p_\imath, \imath\in I_{9}$.

Lemma 4. For each pair $(h, p)\in [0, \bar{h}]\times \mathcal{P},$ it holds that

where $\pi_\imath(t|h, p)\in \mathbb{R}^{5}, \imath\in I_{9},$ is the solution to the following variational system:

Proof. This theorem can be proven in a manner similar to that used for Theorem 5 in ^[43]. ■□

3.5.   Gradient formulas

The state variation with respect to $h$ and $p_\imath, \imath\in I_{9}$ can be computed through solving system (2.1) together with systems (3.30) and (3.32) forward in time. Then, the gradient formulae of ${J}_{\rho, \epsilon}^{\mathbf{H}}(h, p, \psi, \zeta)$ of Problem $\mathbf{H}$ with respect to $h, p, \psi,$ and $\zeta$ can be easily obtained. By Lemmas 3 and 4 and the chain rule, we have

Theorem 3. The gradients of the cost function ${J}_{\rho, \epsilon}^{\mathbf{H}}(h, p, \psi, \zeta)$ with respect to $\lambda\in \Big\{h; p_\imath, \imath\in I_{9}; \psi; \zeta_\iota^\upsilon, \iota\in I_3, \upsilon\in I_\sigma\Big\}$ are given by

where

Remark 2. There exist two steps to calculate the gradients of the cost function ${J}_{\rho, \epsilon}^{\mathbf{H}}(h, p, \psi, \zeta)$ with respect to $\lambda\in \Big\{h; p_\imath, \imath\in I_{9}; \psi; \zeta_\iota^\upsilon, \iota\in I_3, \upsilon\in I_\sigma\Big\}$: (i) solve systems (2.1), (3.30), and (3.32) using a numerical differential scheme, progressing forward in time to obtain their respective solutions; (ii) compute the gradients based on Theorem 3. When these gradients are derived, a gradient-based NLP method can be applied to solve Problem H.

3.6.   Hybrid optimization algorithm

The numerical techniques for solving Problem $\mathbf{H}$ can be divided into two main categories: deterministic methods and intelligent methods ^[46]. Deterministic methods rely on known information, such as gradients, to solve Problem $\mathbf{H}$, and are recognized for their fast convergence. In contrast, intelligent methods often employ probabilistic strategies that are more effective for global optimization. While stochastic methods excel at exploring the global search space, they depend on random search processes, which can be time-consuming.

Typically, gradient-based methods ^[46] are applied to local optimization problems, as they are well-suited for refining solutions near the initial starting point. To address issues like getting stuck in local optima and premature convergence, a particle swarm optimization (PSO) algorithm, an example of an intelligent method, has been introduced ^[41].

To achieve a balance between local refinement and global exploration, a hybrid optimization algorithm, referred to as Algorithm 1, has been developed to effectively address Problem $\mathbf{H}$. In this algorithm, the modified PSO method ^[41] is employed for the exploration phase, while gradient-based techniques are used for local exploitation, as noted in Line 3 of Algorithm 1. This hybrid approach combines the advantages of the PSO algorithm with the efficiency of gradient-based methods, which leverage the gradient information of the objective function in Problem $\mathbf{H}$.

Let $M_{sup}$ be a sufficiently large number. Let $\rho_{sup}$ be a predefined parameter used to ensure the termination of the penalty approach. Let $\rho^0$ and $\epsilon^0$ be the initial values of $\rho$ and $\epsilon$, respectively. Let $\epsilon$ be a predefined parameter used to ensure the termination of Algorithm 1. Let $E_J$ be the convergence tolerance of the cost function. Let $\Psi_{\rho, \epsilon}: = (h, p^\mathrm{T}, \psi, \zeta^\mathrm{T})^\mathrm{T}.$ Let $\Psi_{\rho, \epsilon}^0$ be the initial value of $\Psi_{\rho, \epsilon}$.

Remark 3. Essentially, the goals of Algorithm 1 are for decreasing the value of $\epsilon$ to approximate the $\max$ function and increase the value of $\rho$ to ensure that $|J_{present}-J_{before}| < E_J$ is satisfied.

The convergence of Algorithm 1 is assured as stated in the following. (i) Theorem 1 guarantees that if the smoothing parameter $\epsilon$ is small enough, the solution of Problem $\mathbf{H}$ is an approximate solution of Problem $\mathbf{F};$ (ii) by virtue of the fundamental principles of the penalty function, when $\rho$ is large enough, any solution to Problem $\mathbf{F}$ is a feasible solution to Problem $\mathbf{E}$ ^[40,45]. Note that in Section 3.1, Problem A can be converted into Problem $\mathbf{E}$.

The selection of hyperparameters is often based on a combination of experience, experimentation, and evaluation metrics. For example, inertia weight (w): typically chosen between 0.4 and 0.9; acceleration factors ($c_1, c_2$): generally set between 1.5 and 2; population size: usually ranges from 20 to 50.

4.   Numerical results

4.1.   Experiment design

All computations are performed using MATLAB R2016a on a computer equipped with a 3.70 GHz Intel Core i9-10900K CPU and 32.0 GB of RAM. In the numerical computation, the cultivation conditions have been studied in ^[16]. During the fed-batch culture, alkali, which is intermittently put in the fermentor to maintain its pH value at approximately 7.0, can cause a chemical reaction with acetic acid. This operation significantly affects the correctness of the extracellular concentrations of acetic acid and ethanol. Therefore, in this paper, the focus is on the discrepancy between the measured output data and the system output for the first three substances. The fed-batch culture evolves starting from the initial states $x_{0}.$

The feeding process starts at $\tau^1 = 5.33 \; (\text{h}).$ The feeding time starts at $\tau^{2i+1}, i\in \Lambda_2$ and stops at $\tau^{2i+2}, i\in \Lambda_1$. They were determined through the experiment reported in ^[43]. To obtain the value of the cost function (3.9) in Problem F, system (2.1) needs to be solved numerically. Thus, the computational load is rather large. Therefore, for saving the calculation time, the process of the fed-batch culture reported in ^[47] is grouped into ten stages (i.e., Batch Stage and Stages II–X). In each of Stages II–X, the time duration for each batch stage is 100s minus the duration of the feeding process. The feeding processes in Stages II–X last for 5, 7, 8, 7, 6, 4, 3, 2, 1, 2, and 1 seconds out of every 100 seconds, leaving 95, 93, 92, 93, 96, 97, 98, 99, 98, and 99 seconds for batch processes, respectively.

4.2.   Parameters setting

Based on ^[43], the experiment parameters used in the calculation of the solution of system (2.1) are, respectively, $x_0 = (0.1115\; \text{gL}^{-1}, \ 495 \; \text{mmolL}^{-1}, \ 0, \ 0, \ 0)^{\mathrm{T}}, \underline{x} = (0.0001, 0.1, 0, 0, 0)^{\mathrm{T}}, \overline{x} = (15, 2039,939.5, 1026,360.9)^{\mathrm{T}}, \ r = 0.75, \ C_{s0} = 10,762\; \text{mmolL}^{-1}, \ V_0 = 5\; \text{L}, \ t_f = 26.8\; \text{(h)}, N = 783, \ \overline{h} = 12\; \text{(hour)}, \ \underline{p} = (0.438, 0.5, -7.35435, 0.0039, 33.845, 5.945445, 8.8648, 2.59, 10.225)^{\mathrm{T}}, \overline{p} = (1.314, 3.3, -2.45145, 0.0117,101.535, 17.836335, 26.5944, 7.77, 90.675)^{\mathrm{T}}, \sigma = |\mathcal{T}| = 16, \mathcal{T} = \{0, 2, 4, 5.83, 7.83, 9.83, 11.83, 13.83, 15.83, 17.83, 19.83, 21.83, 23.83, 24.83, 25.83, 26.83\}.$

Euler method is used to solve system (2.1) with a step size of 1/3600. In Algorithm 1, the value of $\rho^0, \epsilon^0, M_{sup}, E_J, \rho_{sup}, \epsilon^{\min}$ are, respectively, $1, 10, 10^7, 10^{-2}, 10^6, 10^{-6}.$ In addition, the value of $\Psi_{\rho, \epsilon}^0$ is selected in its feasible region. The initial function $\phi(t)$ is obtained through the use of the cubic spline method to interpolate the experimental data ^[40].

Remark 4. Practitioners of numerical approximation are most concerned with truncation error, which in a numerical method is an error that is caused by using simple approximations to represent exact mathematical formulas. We seek information about errors on both a local and global scale. Local truncation error is the amount of truncation error that occurs in one step of a numerical approximation. Global truncation error is the amount of truncation error that occurs in the use of a numerical approximation to solve a problem.

We obtain the global truncation error at step n in any numerical approximation of (2.1) by

The local truncation error for Euler's method is $-\frac{(\Delta t)^2}{2}\ddot{x}(\tau)$. The local truncation error has two factors of $\Delta t$, and we say that it is $\mathcal{O}((\Delta t)^2).$ Based on ^[48], we have

where

We can see from the above description that the error vanishes as $\Delta t\rightarrow 0$. Therefore, it is appropriate to use the Euler method for solving differential equations with $\Delta t = {1}/{3600} \; (hour),$ which is small enough for the requirement of microbial fermentation simulation.

4.3.   Results and analysis

By virtue of Algorithm 1, the obtained optimal values of time-delay and system parameters are, respectively, $h^* = 0.824237042558447$ and $p^* = (1.142666910517339, 0.500218711380598, -3.629735870518227, 0.005853192577210, 99.392558636711343, 6.427113527307207, 8.869263769688454, 2.659555099094625, 80.975777396140828)^{\mathrm{T}}$.

The optimal time delay in literature ^[43] was obtained using precise experimental data. In contrast, the optimal time delay presented in this paper is derived under the worst-case distribution of experimental data. Therefore, the optimal time delays obtained in the two cases differ.

The relative errors between the computed values $x_\iota(t_\upsilon|h^*, p^*), t_\upsilon\in \mathcal{T}, \iota\in I_3, \upsilon\in I_\sigma,$ and the experimental data $y_{\iota}(t_\upsilon), t_\upsilon\in \mathcal{T}, \iota\in I_3, \upsilon\in I_\sigma,$ obtained with and without (denoted by "C" and "B", respectively) and taking into consideration of the worst-case probability distribution of the experimental data (WCPDED) are listed in Table 2, where the relative errors are defined by

Here, $x(t_\upsilon| h^*, p^*): = (x_1(t_\upsilon| h^*, p^*), \cdots, x_5(t_\upsilon| h^*, p^*))^\mathrm{T}\in \mathbb{R}_+^5$ denotes the solution of system (2.1) corresponding to the obtained optimal time-delay $h^*$ and system parameters $p^*$ with the initial concentration $x_0,$ and $y(t_\upsilon) = (y_1(t_\upsilon), \cdots, y_3(t_\upsilon))^\mathrm{T}\in \mathbb{R}_+^3$ is the experiment data at the time point $t_\upsilon\in \mathcal{T}.$

Table 2 reveals that the relative error of 'C' exceeds that of 'B'. This discrepancy can be attributed to the fact that the relative error for 'C' is calculated based on the WCPDED. This WCPDED reflects greater variability and uncertainty of the distribution of the experimental data in the measurements for 'C', which naturally leads to a higher relative error. In contrast, surrounding 'B' without taking into consideration of the WCPDED contributes to its comparatively lower relative error. Understanding these distinctions is crucial for accurately interpreting the reliability of our findings.

Under the obtained $h^*$ and $p^*$, the curve of the concentration of the first three substances are drawn in Figures 2–4 (the horizontal axis stands for time; the vertical axis stands for the concentrations of biomass, glycerol, and 1, 3-PD, respectively; A: experimental data; B: the concentration change curve of biomass, glycerol, and 1, 3-PD without considering the WCPDED, respectively; C: the concentration change curve of biomass, glycerol, and 1, 3-PD that takes into account the WCPDED, respectively). The analysis presented in Figures 2–4 illustrates a reasonable alignment between the computed values and the experimental data under the WCPDED framework. To facilitate comparison, we determined the relative errors between the computed concentrations of biomass, glycerol, and 1, 3-PD against the experimental results, utilizing the optimal time-delay and system parameters outlined in ^[43]. It is important to note that, unlike the optimization problem discussed in ^[43], our approach to parameter estimation relies solely on the first-order moment of the experimental data, thereby simplifying the process. Moreover, Algorithm 1 is designed to operate effectively without requiring precise knowledge of the distribution of the experimental data. The entire calculation process took approximately 15,653 seconds.

Figures 2–4 demonstrate that the NCTSTD system, utilizing the parameters $h^*$ and $p^*$, effectively captures the dynamics of the fed-batch culture process, even within the constraints of the WCPDED. This suggests that our model not only aligns well with observed behaviors but also provides a robust framework for understanding the intricacies of the cultivation process. The ability of the NCTSTD system to maintain accuracy under these conditions highlights its potential for applications in optimizing bioprocesses.

Note that our method is under the assumption that the first-order moment information of the experimental data can be accurately obtained. Obviously, it will be practically important to extend this method to cases where moment information is also uncertain. This will be an interesting further research direction.

5.   Discussion

This paper studies an NCTSTD with unknown quantities, in which the measured output data is uncertain. The focus is to identify these unknown quantities. For this, we formulate the problem as a DRPE problem subject to the NCTSTD system. A hybrid algorithm is presented to solve the DRPE problem. Finally, numerical results are obtained, which demonstrate the validity of the designed algorithm.

On the other hand, it is important to note that the designed algorithm is based on the assumption that the first-order moment information of the measured output data with respect to its probability distribution is available, highlighting the need to extend the algorithm to cases where this moment information is uncertain. This would be an interesting further research area.

In summary, the designed algorithm provides a promising method to optimally identify time-delay and system parameters in the fed-batch culture process of glycerol bioconversion to 1, 3-PD under the environment of uncertain measured output data. The potential extension of this algorithm to more general cases would enhance its scope of applications.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This work was supported in part by the National Key Research and Development Program of China under Grant 2022YFB3304600; in part by the Nature Science Foundation of Liaoning Province of China under Grant 2024-MS-015; in part by the Fundamental Research Funds for the Central Universities under Grant 3132024196, and Grant DUT22LAB305; in part by the National Natural Science Foundation of China under Grant 12271307, and Grant 12161076; in part by the China Postdoctoral Science Foundation under Grant 2019M661073; in part by the Xinghai Project of Dalian Maritime University.

Conflict of interest

The authors declare there is no conflicts of interest.