Dynamic optimization of a two-stage fractional system in microbial batch process

1.   Introduction

1,3-Propanediol (1,3-PD) is a valuable chemical raw material widely used in various industries, including food, cosmetics, and pharmaceuticals ^[1]. Two commonly used methods for producing 1,3-PD are chemical synthesis and microbial fermentation ^[2]. Compared with chemical synthesis, microbial fermentation offers distinct advantages, such as operational simplicity and minimal byproduct formation. These benefits have garnered increasing attention from both academia and industry. Given the impracticality of conducting numerous laboratory experiments to achieve high 1,3-PD conversion, it is essential to study dynamic optimization models for microbial fermentation.

There are three modes of glycerol bioconversion to 1,3-PD: batch culture, continuous culture, and fed-batch culture. In batch culture, the microorganism (Klebsiella pneumoniae) and the substrate (glycerol) are injected into the bioreactor only once at the beginning, with no further additions or removals during the culture period. In continuous culture, the substrate is continuously injected into the bioreactor at a specific rate, while the culture fluid is simultaneously removed at the same rate. In fed-batch culture, the operation alternates between batch and feeding phases. In this paper, we will focus on batch culture because of its ease of operation and high yield of 1,3-PD, while also laying a theoretical foundation for continuous and fed-batch culture processes ^[3]. To enhance the understanding of the microbial batch process, a nonlinear kinetic system for substrate consumption and product formation is proposed in ^[4]. A parameter identification problem with unknown time-delay and system parameters is introduced in ^[5]. Based on these mathematical models, a robust dynamic optimization problem with respect to uncertain system parameters is studied in ^[6]. A dynamic optimization problem with a cost sensitivity constraint is developed in ^[7]. A bi-objective dynamic optimization problem is discussed in ^[8], and a stochastic dynamic optimization problem is presented in ^[9]. More recently, a robust dynamic optimization problem governed by a nonlinear switched time-delay system involving an unknown time-varying function is formulated in ^[10]. A multi-objective dynamic optimization problem subject to a nonlinear time-delay system aimed at balancing system cost and system sensitivity is considered in ^[11]. Furthermore, artificial intelligence methods, such as deep learning and reinforcement learning, have been explored in various fields; see, for example, ^[12,13,14]. In ^[12], a multifactor prediction model incorporating a combined normalization layer and a codec is proposed to effectively address data differences and complex nonlinearity in the prediction of industrial wastewater pollutants. An antimicrobial peptide screening model based on long short-term memory neural networks with an attention mechanism is presented in ^[13]. In ^[14], multi-objective reinforcement learning is employed to obtain Pareto optimal solution sets for each objective in the control of fed-batch fermentation processes. Although the aforementioned results are of interest, they are restricted to dynamical systems with integer-order derivatives.

Fractional derivatives extend integer-order derivatives to non-integer orders, providing a powerful tool for modeling and analyzing complex systems ^[15]. Unlike integer-order derivatives, fractional derivatives are regarded as nonlocal operators with memory and hereditary because they take into account a broader range of historical influences. As a result, fractional dynamical systems are well-suited for representing complex phenomena characterized by memory effects. In recent years, many successful fractional models have been studied in bioengineering research. For example, in ^[16], fractional differential equations are used to model biological reactive systems such as tequila production, bioethanol production, and the thermal hydrolysis process, demonstrating the feasibility and effectiveness of fractional calculus in modeling biological process systems. A fractional mathematical model for erythritol and mannitol synthesis is established in ^[17], which proves to be useful for both prediction and process optimization. In ^[18], a novel cascaded control strategy based on fractional-order fuzzy proportional-derivative/proportional-integral control is proposed for temperature regulation in the fermentation process. A nonlinear fractional Michaelis-Menten enzyme kinetics model is introduced in ^[19], and a homotopy perturbation method is developed to effectively solve this biochemical reaction process. It is worth noting that fractional modeling has also been gaining popularity in the bioconversion of glycerol to 1,3-PD. In ^[20], a fractional parameter identification problem is considered in the continuous culture, where the dynamical system is solved by applying the trapezoidal method and the predictor-corrector method. A single-stage fractional dynamical system with unknown kinetic parameters and fractional orders is proposed to describe the batch culture in ^[21]. Given that the kinetic behavior of the stationary phase described in ^[21] does not align well with experimental results, a new parameter identification problem involving a two-stage fractional dynamical system with different fractional orders and kinetic parameters is introduced in ^[22]. However, to the best of our knowledge, no study on the dynamic optimization of a two-stage fractional system in 1,3-PD batch production has been reported in the literature.

Motivated by this, in this paper, we propose a dynamic optimization problem involving the two-stage fractional system in ^[22] to optimize the microbial batch process, aiming to maximize the productivity of 1,3-PD at the terminal time. The decision variables of this process include the initial concentrations of biomass and glycerol, along with the terminal time of the batch process. By applying a proposed time-scaling transformation, we equivalently transform the problem with free terminal time into one with fixed terminal time in a new time horizon. By using an exact penalty function method, we convert the equivalent problem, which involves both a terminal state inequality constraint and continuous state inequality constraints, into an optimization problem with only box constraints. We then present a novel third-order numerical scheme to solve the two-stage fractional system. An improved particle swarm optimization IPSO algorithm is developed to determine the optimal decision variables. Numerical results demonstrate that the productivity of 1,3-PD at the terminal time is higher than the previous results, verifying the effectiveness of the proposed optimization strategy.

The rest of the paper is structured as follows. A two-stage fractional system is presented in Section 2. In Section 3, a dynamic optimization problem is proposed. Section 4 develops a third-order numerical scheme and an IPSO algorithm. In Section 5, numerical results are provided. Finally, concluding remarks are given in Section 6.

2.   Two-stage fractional system in microbial batch process

Based on the work presented in ^[22], mass balance equations between biomass, substrate, and products in the microbial batch process can be modeled by the following two-stage fractional system:

where $\mathit{\boldsymbol{\alpha}}^{1} = (\alpha_{1}^{1}, \alpha_{2}^{1}, \ldots, \alpha_{5}^{1})^{\top}\in (0, 1]^{5}$ and $\mathit{\boldsymbol{\alpha}}^{2} = (\alpha_{1}^{2}, \alpha_{2}^{2}, \ldots, \alpha_{5}^{2})^{\top} \in (0, 1]^{5}$ are given vectors of fractional orders; $t$ is the process time (in hours); $\mathit{\boldsymbol{x}}(t) = (x_{1}(t), x_{2}(t), \ldots, x_{5}(t))^{\top}\in \mathbb{R}^{5}$ is the state vector representing the concentrations of biomass, glycerol, 1,3-PD, ethanol, and acetate, respectively, with $x_{1}(t)$ measured in gL$^{-1}$ and $x_{i}(t), i = 2, 3, 4, 5$, measured in mmolL$^{-1}$; $t_{1}$ is a given switching time; $t_{f}$ is a free terminal time; $\mathit{\boldsymbol{x}}^{0}\in \mathbb{R}^{5}$ is the initial state vector; $\mathit{\boldsymbol{x}}(t_{1}^{+})$ denotes the righthand limit of $\mathit{\boldsymbol{x}}(t_{1})$; $_{0}^{C}D_{t}^{\mathit{\boldsymbol{\alpha}}^{1}}\mathit{\boldsymbol{x}}(t) = \Big({_{0}^{C}}D_{t}^{\alpha_{1}^{1}}x_{1}(t), \ldots, {_{0}^{C}}D_{t}^{\alpha_{5}^{1}}x_{5}(t)\Big)^{\top}$ and $_{t_{1}}^{C}D_{t}^{\mathit{\boldsymbol{\alpha}}^{2}}\mathit{\boldsymbol{x}}(t) = \Big({_{t_{1}}^{C}}D_{t}^{\alpha_{1}^{2}}x_{1}(t), \ldots, {_{t_{1}}^{C}}D_{t}^{\alpha_{5}^{2}}x_{5}(t)\Big)^{\top}$ with $_{0}^{C}D_{t}^{\alpha_{i}^{1}}x_{i}(t)$ and $_{t_{1}}^{C}D_{t}^{\alpha_{i}^{2}}x_{i}(t), i = 1, 2, \ldots, 5,$ denoting the $\alpha_{i}^{1}$th and $\alpha_{i}^{2}$th Caputo fractional derivatives of $x_{i}(t)$ defined by

with $\Gamma(\cdot)$ representing the Gamma function and $\dot{x}_{i}(\tau)$ denoting the first-order derivative of $x_{i}(\tau)$; and $\mathit{\boldsymbol{f}}^{1}: \mathbb{R}^{5}\to \mathbb{R}^{5}$ and $\mathit{\boldsymbol{f}}^{2}: \mathbb{R}^{5}\to \mathbb{R}^{5}$ describe a two-stage dynamical system in the microbial batch process, as detailed below:

Here, $c_{j}$ and $c_{j+5}, j = 1, 2, \ldots, 5$, are given kinetic parameters that represent biomass growth, glycerol consumption, and the formation of 1,3-PD, ethanol, and acetate, respectively. Similarly, $d_{j}$ and $d_{j+5}, j = 1, 2, \ldots, 5$, are given kinetic parameters that represent the inhibitory effects of cell death on these same processes. Under anaerobic conditions at $37\; ^\circ \mathrm{C}$ and pH 7.0, the values of these fractional orders and kinetic parameters are listed in Table 1. Furthermore, Figure 1 depicts substrate consumption and product formation in the microbial batch process.

In the two-stage fractional system (1) and (2), the initial concentrations of the products (i.e., 1,3-PD, ethanol, and acetate) are set to $x_{i}^{0} = 0$ for $i = 3, 4, 5$, since no products are presented at the initial time point. In contrast, the initial concentrations of biomass and glycerol are treated as decision variables to be optimized. Let $\mathit{\boldsymbol{\theta}} = (x_{1}^{0}, x_{2}^{0})^{\top}$ and define

where $a_{m}^{\min}$ and $a_{m}^{\max}, m = 1, 2$, are given real numbers such that $a_{m}^{\min}\leq a_{m}^{\max}$. In addition, the terminal time $t_{f}$ is also considered as a decision variable. Define

where $b^{\min}$ and $b^{\max}$ are given real numbers such that $t_{1} < b^{\min}\leq b^{\max}$. Any pair $(\mathit{\boldsymbol{\theta}}, t_{f})\in \Theta \times \mathcal{T}$ is called an admissible pair for the two-stage fractional system (1) and (2).

For the two-stage fractional system (1) and (2), there exists a unique absolutely continuous solution denoted by $\mathit{\boldsymbol{x}}(\cdot|\mathit{\boldsymbol{\theta}}, t_{f})$ for each $(\mathit{\boldsymbol{\theta}}, t_{f})\in \Theta \times \mathcal{T}$ ^[23]. Moreover, in consideration of the practical production process, it is biologically meaningful to restrict the concentrations of biomass, glycerol, and products within the specified set $W$:

where $x_{*1} = 0.01$ and $x_{*i} = 0, i = 2, 3, 4, 5$, represent the lower thresholds for the concentrations of biomass, glycerol, 1,3-PD, ethanol, and acetate, respectively, which are required for cell growth; and the corresponding upper concentration thresholds are $x_{1}^{*} = 10, x_{2}^{*} = 2039, x_{3}^{*} = 939.5, x_{4}^{*} = 360.9$, and $x_{5}^{*} = 1026$ ^[22].

3.   Dynamic optimization problem

In the microbial batch process, glycerol is the substrate and 1,3-PD is the target product. Thus, we aim to maximize the productivity of 1,3-PD at the terminal time while minimizing the glycerol consumption rate. To achieve this, we consider the following productivity objective function:

and the glycerol consumption rate constraint:

where $\theta_{2}$ denotes the initial concentration of glycerol, and $\Omega > 0$ is a predefined positive real number representing the maximum allowable glycerol consumption rate. Therefore, by combining (6)–(8), the dynamic optimization problem involving the two-stage fractional system (1) and (2) can be stated as follows:

Note that Problem (P) is a nonstandard dynamic optimization problem with two notable characteristics: (ⅰ) the terminal time $t_{f}$ is free rather than fixed, and (ⅱ) constraints (6) are continuous state inequality constraints.

To effectively handle the nonstandard characteristic (ⅰ), we will transform Problem (P), which is of a free terminal time, into an equivalent problem with a fixed terminal time in a new time horizon by employing the following time-scaling transformations. This procedure is implemented in three steps.

Step 1. $t\in [0, t_{1}]\to s\in[0, 1]$

Let

Clearly, $s = 0$ and $s = 1$ correspond to $t = 0$ and $t = t_{1}$, respectively. Furthermore, let $\tau = t_{1}\eta$ and $y_{i}(\eta) = x_{i}(t_{1}\eta)$. Then, we rewrite (3) as

for $i = 1, 2, \ldots, 5$. Let $\mathit{\boldsymbol{y}}(s) = \mathit{\boldsymbol{x}}(t_{1}s)$ and ${_{0}^{C}}D_{s}^{\mathit{\boldsymbol{\alpha}}^{1}}\mathit{\boldsymbol{y}}(s) = \Big({_{0}^{C}}D_{s}^{\alpha_{1}^{1}}y_{1}(s), \ldots, {_{0}^{C}}D_{s}^{\alpha_{5}^{1}}y_{5}(s)\Big)^{\top}$. Then, fractional system (1) becomes

where $t_{1}^{\mathit{\boldsymbol{\alpha}}^{1}} = \Big(t_{1}^{\alpha_{1}^{1}}, t_{1}^{\alpha_{2}^{1}}, \ldots, t_{1}^{\alpha_{5}^{1}}\Big)^{\top}$; and $\circ$ represents the Hadamard product, which denotes element-wise multiplication between vectors or matrices of the same dimensions.

Step 2. $t\in (t_{1}, t_{f}]\to s\in (1, 2]$

Define

where $b^{\min}$ and $b^{\max}$ are as defined in (5).

Let

Clearly, $s = 2$ corresponds to $t = t_{1}+\sigma = t_{f}$. Furthermore, let $\tau = t_{1}+\sigma(\eta-1)$ and $y_{i}(\eta) = x_{i}(t_{1}+\sigma(\eta-1))$. Then, we rewrite (4) as

for $i = 1, 2, \ldots, 5$. Let $\mathit{\boldsymbol{y}}(s) = \mathit{\boldsymbol{x}}(t_{1}+\sigma(s-1))$ and ${_{1}^{C}}D_{s}^{\mathit{\boldsymbol{\alpha}}^{2}}\mathit{\boldsymbol{y}}(s) = \Big({_{1}^{C}}D_{s}^{\alpha_{1}^{2}}y_{1}(s), \ldots, {_{1}^{C}}D_{s}^{\alpha_{5}^{2}}y_{5}(s)\Big)^{\top}$. Then, fractional system (2) becomes

where $\sigma^{\mathit{\boldsymbol{\alpha}}^{2}} = \Big(\sigma^{\alpha_{1}^{2}}, \sigma^{\alpha_{2}^{2}}, \ldots, \sigma^{\alpha_{5}^{2}}\Big)^{\top}$.

Step 3. The restatement of Problem (P)

Let $\mathit{\boldsymbol{y}}(\cdot|\mathit{\boldsymbol{\theta}}, \sigma)$ denote the solution of the two-stage fractional system (9) and (10) for each $(\mathit{\boldsymbol{\theta}}, \sigma)\in \Theta \times \Sigma$. Under the time-scaling transformations (Steps 1 and 2), constraints (6) can be rewritten as

the productivity objective function (7) is transformed into

and the glycerol consumption rate constraint (8) is reformulated as

Therefore, Problem (P) can be restated as the following equivalent dynamic optimization problem with fixed terminal time:

It is worth noting that Problem (EP) is a dynamic optimization problem subject to the constraints (11) and (13). Handling constraint (11) is particularly challenging numerically because it must be satisfied at an infinite number of points over the entire time horizon. The penalty function method is a commonly used technique for handing this type of constrained optimization problems. It transforms a constrained problem into a sequence of unconstrained ones by incorporating a penalty term into the objective function. In particular, the exact penalty function method only requires a sufficiently large but finite penalty parameter to obtain a solution that satisfies the constraints and achieves the optimal objective. Therefore, to surmount the nonstandard characteristic (ⅱ), we will employ an exact penalty function method to effectively handle these constraints.

Let

Then, constraint (11) is equivalent to

Furthermore, we rewrite the constraint (13) as

Combining (14) and (15) yields the following constraint violation function:

Obviously, $\hat{H}(\mathit{\boldsymbol{\theta}}, \sigma) = 0$ ensures that the constraints (11) and (13) are satisfied. By incorporating $\hat{H}(\mathit{\boldsymbol{\theta}}, \sigma)$ as a penalty term into the objective function (12), we propose the following penalty problem:

where $\beta > 0$ is a positive penalty parameter. Furthermore, by using a similar derivation as presented in ^[26], it can be proved that $\hat{J}_{\beta}(\mathit{\boldsymbol{\theta}}, \sigma)$ is an exact penalty function.

4.   Computational approach

In this section, we present a third-order numerical scheme for discretizing the two-stage fractional system (9) and (10). Based on this scheme, we introduce an IPSO algorithm to determine the optimal decision variables for Problem $({\rm EP}_{\beta})$.

4.1.   Numerical scheme

For given positive integers $N^{l}, l = 1, 2$, we divide the intervals $(l-1, l]$ into $N^{l}$ subintervals $(s_{q-1}^{l}, s_{q}^{l}]$ with partition points $s_{q}^{l} = (l-1)+qh^{l}$, where $q = 1, 2, \ldots, N^{l}$, and $h^{l} = 1/N^{l}$. Specifically, $s_{0}^{1} = 0$ and $s_{0}^{2} = 1$. Then, the two-stage fractional system (9) and (10) can be transformed into the following Volterra integral equations ^[15]:

for $i = 1, 2, \ldots, 5$, where

Now, we subdivide the intervals $(s_{0}^{l}, s_{1}^{l}]$ for $l = 1, 2$, into $\aleph^{l}$ subintervals $(\varpi_{\kappa-1}^{l}, \varpi_{\kappa}^{l}]$ with partition points $\varpi_{\kappa}^{l} = s_{0}^{l}+\kappa \hbar^{l}, \kappa = 1, 2, \ldots, \aleph^{l}$, where $\aleph^{l} = \big\lceil 1/(h^{l})^{\frac{1}{2}}\big\rceil$ with $\lceil \cdot \rceil$ representing the ceiling function, and $\hbar^{l} = h^{l}/\aleph^{l}$. Specifically, $\varpi_{0}^{l} = s_{0}^{l}$. For any $l = 1, 2$, and $i = 1, 2, \ldots, 5$, we approximate $f_{i}^{l}(\mathit{\boldsymbol{y}}(\eta))$ on the righthand side of (18) in $(\varpi_{\kappa-1}^{l}, \varpi_{\kappa}^{l}]$ for $\kappa = 1, 2, \ldots, \aleph^{l}$, by the following second-order Taylor expansion:

where $\Lambda^{0}_{\kappa}f_{i}^{l}$ and $\Lambda^{1}_{\kappa}f_{i}^{l}$ are coefficients to be determined, and $\Lambda^{2}_{\kappa}f_{i}^{l}$ is the remainder term. Omitting this remainder term, we can use the values of $f_{i}^{l}(\mathit{\boldsymbol{y}}(\varpi_{\kappa-1}^{l}))$ and $f_{i}^{l}(\mathit{\boldsymbol{y}}(\varpi_{\kappa}^{l}))$ to determine $\Lambda^{0}_{\kappa}f_{i}^{l}$ and $\Lambda^{1}_{\kappa} f_{i}^{l}$, yielding the following divided differences:

Substituting (19) into the first term of the sums of (18), we obtain

where

Then, for any $l = 1, 2$, and $i = 1, 2, \ldots, 5$, we approximate $f_{i}^{l}(\mathit{\boldsymbol{y}}(\eta))$ on the righthand side of (18) in $(s_{j}^{l}, s_{j+1}^{l}]$ for $j = 1, 2, \ldots, q-1$, and $q = 2, 3, \ldots, N^{l}$, by the following third-order Taylor expansion:

where $\Delta^{0}_{j+1}f_{i}^{l}, \Delta^{1}_{j+1}f_{i}^{l}$, and $\Delta^{2}_{j+1} f_{i}^{l}$ are coefficients to be determined; and $\Delta^{3}_{j+1}f_{i}^{l}$ is the remainder term. Omitting this remainder term, we can use the values of $f_{i}^{l}(\mathit{\boldsymbol{y}}(s_{j-1}^{l})), f_{i}^{l}(\mathit{\boldsymbol{y}}(s_{j}^{l}))$, and $f_{i}^{l}(\mathit{\boldsymbol{y}}(s_{j+1}^{l}))$ to determine $\Delta^{0}_{j+1}f_{i}^{l}, \Delta^{1}_{j+1}f_{i}^{l}$, and $\Delta^{2}_{j+1} f_{i}^{l}$, yielding the following divided differences:

Substituting (21) into the $(j+1)$th term of the sums of (18), we obtain

where

Combining (18), (20), and (22), and then omitting the truncation errors from all remainder terms, we obtain the following numerical scheme:

for $l = 1, 2, i = 1, 2, \ldots, 5$, and $q = 1, 2, \ldots, N^{l}$, where $\mathit{\boldsymbol{y}}^{l, \kappa, 0}$ and $\mathit{\boldsymbol{y}}^{l, q}$ denote the approximations of $\mathit{\boldsymbol{y}}(\varpi_{\kappa}^{l})$ and $\mathit{\boldsymbol{y}}(s_{q}^{l})$ for each feasible $l, \kappa$, and $q$.

Based on (23), it is evident that the aforementioned numerical scheme is implicit and can be efficiently solved using Newton's method ^[24], which is a well-established technique for solving nonlinear equations. The following theorem presents the convergence rate of this numerical scheme.

Theorem 4.1. Let $\hat{\mathit{\boldsymbol{y}}}^{l, q}$ be the numerical solution of (23) obtained using Newton's method for $l = 1, 2$, and $q = 1, 2, \ldots, N^{l}$. Then, there exists a positive constant $\hat{h} > 0$ such that $h < \hat{h}$, and the following inequality is satisfied:

where $\|\cdot\|_{\infty}$ denotes the infinity norm; $C > 0$ is a positive constant independent of $\mathit{\boldsymbol{\alpha}}^{l}$; and $h: = \max\limits_{l\in \{1, 2\}}\{h^{l}\}$.

Proof. The proof is similar to those given for Theorems 3.1–3.3 in ^[25]. □

4.2.   Optimization algorithm

Problem $({\rm EP}_{\beta})$ is a parameter optimization problem with only box constraints. While gradient-based optimization methods ^[27,28] can be used to solve this problem, they are prone to getting trapped in local optima. Therefore, in this subsection, we will develop an IPSO algorithm to solve Problem $({\rm EP}_{\beta})$.

PSO, a global heuristic algorithm inspired by the social behavior of bird flocking and fish schooling during food searches, is widely regarded as one of the most successful nature-inspired optimization algorithms. Compared to other heuristic algorithms, such as the genetic algorithm, bat algorithm, ant colony optimization, and ecosystem based optimization, PSO is highly efficient and adaptable to a variety of dynamic environments. Due to its simple structure and minimal algorithmic parameters, PSO has gained widespread popularity across diverse fields such as power systems, control systems, networking, image segmentation, and more ^[29]. However, the standard PSO (SPSO) algorithm ^[30] is susceptible to falling into local optima. To enhance both its local and global search capabilities, we introduce a new strategy ^[31,32] for updating particle velocity and position by adjusting algorithm parameters, including inertia weight, cognitive factor, and social factor. For the box constraints in Problem $({\rm EP}_{\beta})$, we also propose a new strategy to handle constraint violations. For brevity, let $\mathit{\boldsymbol{\varrho}} = (\theta_{1}, \theta_{2}, \sigma)^{\top}\in \mathbb{R}^{3}$ represent the decision vector of Problem $({\rm EP}_{\beta})$ in the IPSO algorithm. The main steps of the IPSO algorithm are given in Algorithm 1.

5.   Numerical results

In the numerical simulation, Algorithm 1 is implemented in a Matlab program to solve Problem $({\rm EP}_{\beta})$, with all computations performed on a personal computer equipped with a 2.80 GHz CPU and 16.0 GB of RAM. To solve Problem $({\rm EP}_{\beta})$, the initial decision vector, switching time, maximum allowable glycerol consumption rate, and number of subintervals for systems (9) and (10) are set as $(0.2245\, {\rm gL^{-1}}, 509.8913\, {\rm mmolL^{-1}}, 2.0\, {\rm h})^{\top}$, $t_{1} = 5.75\, {\rm h}$, $\Omega = 100$ mmolL$^{-1}$h$^{-1}$, $N^{1} = 500$, and $N^{2} = 100$, respectively. In Algorithm 1, the parameters $\bar{N}$, $I_{\rm \max}$, $\xi$, $\omega_{\min}$, $\omega_{\max}$, $\mathit{\boldsymbol{\varrho}}_{\rm low}$, $\mathit{\boldsymbol{\varrho}}_{\rm upp}$, and $\beta$ are set as 50,100, 10$^{-4}$, 0.4, 0.9, $(0.05,200, 0.01)^{\top}$, $(1, 1700, 5)^{\top}$, and 100, respectively. It should be noted that the parameters in Algorithm 1 are empirically determined through numerous numerical experiments. By running Algorithm 1, we obtain the optimal decision vector $(0.05\, {\rm gL^{-1}}, 765.2177\, {\rm mmolL^{-1}}, 0.8576\, {\rm h})^{\top}$. As a result, the terminal time of the batch process is reduced to 6.6076 h, representing a 14.7406% decrease compared with the experimental data of 7.75 h in ^[4].

Based on the obtained optimal decision vector, the productivity of 1,3-PD at the terminal time is determined to be 61.8505 mmolL$^{-1}$h$^{-1}$, showing a 79.3694% improvement compared with the 34.4822 mmolL$^{-1}$h$^{-1}$ observed in the experimental data from ^[4], and a 19.1495% improvement compared with the 51.91 mmolL$^{-1}$h$^{-1}$ achieved by the single-stage integer-order system in ^[6]. The glycerol consumption rate in our optimization model is 100 mmolL$^{-1}$h$^{-1}$, which is comparable to the 100.8 mmolL$^{-1}$h$^{-1}$ reported in ^[6]. Figure 2 illustrates the corresponding changes in the concentrations of biomass, glycerol, 1,3-PD, ethanol, and acetate under the optimal decision vector. Figure 3 depicts the curve of 1,3-PD productivity obtained by using the optimal decision vector. For comparison, changes in 1,3-PD productivity based on the experimental data ^[4] and the single-stage integer-order system ^[6] are also plotted in Figure 3, from which we observe that our optimization model achieves the highest productivity of 1,3-PD at the terminal time.

To further evaluate the performance of the IPSO algorithm, the SPSO algorithm is also developed to solve Problem $({\rm EP}_{\beta})$, with the inertia weight, cognitive factor, and social factor set to be 0.5, 2, and 2, respectively. We conduct 50 independent tests for both the IPSO and SPSO algorithms. The obtained optimal cost, worst cost, average cost, and average iteration time are listed in Table 2, from which we see that the IPSO algorithm outperforms the SPSO algorithm in terms of optimal cost, worst cost, and average cost. The average iteration time of the IPSO algorithm is comparable to that of the SPSO algorithm. The convergence curves of the objective function for both the IPSO and SPSO algorithms are depicted in Figure 4. From Figure 4, we observe that the IPSO algorithm converges significantly faster than the SPSO algorithm.

6.   Conclusions

In this paper, we have proposed the dynamic optimization problem involving the two-stage fractional system subject to the terminal state inequality constraint and continuous state inequality constraints. The objective function is the productivity of 1,3-PD at the terminal time, while the decision variables are the initial concentrations of biomass and glycerol, as well as the terminal time of the batch process. The main contributions of our work include: ⅰ) the dynamic optimization problem under consideration is equivalently transformed into the one with fixed terminal time and only box constraints by applying the time-scaling transformation and the exact penalty function method; ⅱ) we present the novel third-order numerical scheme to solve the two-stage fractional system; and ⅲ) we develop the highly effective IPSO algorithm to determine the optimal decision variables. The numerical results show a significant increase in the productivity of 1,3-PD at the terminal time compared to previously reported outcomes, providing robust support for practical 1,3-PD batch production. In future research, we will validate the findings of this work in the real-world fermentation process. It is worth noting that, in this paper, the switching time between the two-stage fractional system is fixed. In the future, we will extend our developed method to solve the dynamic optimization of switched systems in the microbial batch process. Additionally, integrating the third-order numerical scheme with deep reinforcement learning could be a promising approach for solving large-scale fractional dynamic optimization problems.

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 is supported by the Ministry of Higher Education of Malaysia, Fundamental Research Grant Scheme (FRGS/1/2021/STG06/SYUC/03/1), the National Natural Science Foundation of China (No. 12271307), and the Natural Science Foundation of Shandong Province, China (No. ZR2023MA054).

Conflict of interest

The authors declare that they have no competing interests.