An Optimized Fed-Batch Culture Strategy Based on Multidimensional Time Series Aggregation

Abstract

An optimized feeding strategy for the fed-batch fermentation of 2-keto-L-gulonic acid (2-KGA) is proposed to enhance the total profit of a workshop. Instead of assigning the same quantity of substrate to each batch in the empirical scheduling, the online feeding strategy offers a practical approach that allows for the precise management of substrate feeding to maximize economic performance. The forecasting of the profit function has been achieved through the proposed similarity measurement and multidimensional time series aggregation method. This method solves the problem of predicting the total economic benefit at the end of fermentation, in the absence of fermentation substrate measurement data. Pseudo-online simulation has been conducted, utilizing data derived from ninety industrial batches. A 6% increase in total profit is observed, which indicates that the optimized strategy utilizes the substrate and other related resources more efficiently and leads to a higher overall profit.

1. Introduction

Batch operation is a common mode in chemical process industries, such as pharmaceuticals and food processing, for flexibility in process operations and adaptability to changing market scenarios. Great attention has been paid to batch process optimization [1,2,3], which is motivated by increasing competition in the manufacturing sectors.

This study focuses on the process optimization of 2-keto-L-gulonic acid (2-KGA) fed-batch fermentation, a precursor for the synthesis of L-ascorbic acid. L-ascorbic acid (vitamin C) is a vital nutrient for humans and select animal species due to its antioxidant properties and physiological functions [4]. More importantly, vitamin C may play a role against COVID-19 infection through multiple pathways of the immune system and inflammatory response [5,6].

The industrial production of L-ascorbic acid is achieved through two primary routes, each with a distinct historical origin. The Reichstein process, developed in the 1930s [7], involves a single pre-fermentation step followed by chemical synthesis. In contrast, the modern two-step fermentation process, pioneered in China in the 1970s [8], represents a more streamlined and cost-effective approach, relying on microbial biosynthesis in conjunction with chemical processing. The 2-keto-L-gulonic acid (2-KGA) cultivation under study is practiced in a Chinese pharmaceutical plant using a mixed culture of Ketogulonicigenium vulgare (K. vulgare) and Bacillus megaterium (B. megaterium), with L-sorbose as the substrate. In the pure culture of K. vulgare, cells grow very slowly [9]. However, under mixed culture conditions with B. megaterium, the growth rate of K. vulgare is significantly increased [10]. The relationship between the two strains has been a hot topic in the field of vitamin C fermentation [11]. The help from B. megaterium covers many aspects [12], such as the supply of purines [13] and pyrimidines [14], erythrose and erythritol [15], and the initiation of the reproduction of K. vulgare [16].

The industrial mixed culture is carried out in air-lift bioreactors with a volume of 200 m³. The initial medium standard recipe contains L-sorbose, corn steep liquor, carbamide, KH₂PO₄, and MgSO₄•7H₂O. During cultivation, the alkali solution is automatically fed to control the pH value at about 7.0. Key variables, including temperature, pH, aeration rate, and liquid level, are continuously monitored in real-time using the Distributed Control System (DCS). To assess the concentration of 2-KGA, offline iodometry titration is employed at regular intervals. Similarly, the L-sorbose concentration is measured using an offline anthrone colorimetric method. These analytical techniques provide valuable insights into the behavior of the system and enable informed decision making with respect to process optimization and control.

To avoid substrate inhibition, about half of the L-sorbose allocated to each batch is added to the initial media, while the rest is fed into the bioreactor under a pulse-feeding pattern during industrial production. In traditional practice, the allocation of L-sorbose resources is uniformly distributed across all batches in the empirical scheduling process. However, this approach may not be optimal as it fails to account for potential variations in specific individual batches. Despite following a standardized recipe, variations in process performance, including cultivation time and product yield, have been observed among batches. Such variations and fluctuations are unavoidable to some extent, which also implies that an optimized feeding strategy to enhance economic efficiency is needed. The strategy should allocate the feeding amount of the substrate and other related resources more efficiently, leading to a higher overall profit.

An example of industrial penicillin fed-batch production is presented to demonstrate the establishment of the profit index, or profit function [17]. By utilizing the definition of the profit function, the profit potential can be computed, and a scheduling approach to maximize the economic profit of the whole workshop is proposed [18]. The fuzzy decision-making method has also been widely used for process scheduling [19,20,21,22,23].

However, specific to the industrial production of 2-KGA, many offline variables are not periodically measured considering the cost of production, such as the substrate and biomass concentration. This makes it impossible to accurately calculate the profit function during the fermentation process. Since the feeding decision should be made in the middle stages of the fermentation process, the crux of the feeding strategy involves the assessment of the profit-making potential for a given online batch at that precise moment.

To solve this problem, one feasible solution is to dig into the potential in the time series of the multidimensional process variables. This encompasses two phases: the aggregation of individual variables into a decision-making index and the utilization of the index for an optimized feeding strategy. For the information aggregation of time series, similarity measurement is an effective approach [24,25,26]. The shape-based method of Dynamic Time Warping (DTW) has been widely applied to various fields [27,28,29].

DTW can find the minimum path by providing non-linear alignments between the two time series [30]. For instance, a method is presented for predicting the harvest time of a bacteria fermentation batch process at an early stage in the process using DTW and lasso regression [31]. Wan et al. proposed a novel approach to align uneven batch data by identifying short-window PCA and PLS models at first and then applying these identified models to extend shorter trajectories and predict future batch end-product quality [32]. Gunasekaran et al. employed the DTW algorithm to analyze diverse foot movements captured from wearable IoT devices [33], whereas Thanchanok et al. utilized a shape-based approach with the DTW method to more accurately predict consumer energy consumption patterns at the household level [34].

In this study, an online feeding strategy for 2-KGA fed-batch fermentation is proposed where enhancing the total profit of the whole process is addressed. The profit-making potential of an online batch is evaluated using time series similarity matching with all historical batches. Firstly, the time series similarity measurement is performed based on the time series of key process variables (accumulated alkali consumption and product formation at the sampled time). Then, all the measurements between the online batch and historical batches are sorted, and the top k historical batches are selected. Finally, the weighted profit functions of the chosen batches are calculated using the proposed method. This aggregation method makes the batches with a high degree of time series similarity have a greater weight in the prediction. Based on the predicted profit function, an optimized feeding strategy is proposed and high profit-making batches are fed with more resources to extend cultivation time and make more profit. The validation of the proposed approach using data from an industrial-scale cultivation process is presented and discussed.

2. Materials and Methods

2.1. The Optimized Objective

The profit function is used to evaluate the economic performance, which is defined as the gross profit of a fermentation batch over the production time:

$$J\left(t\right)=\frac{P\left(t\right)P_{sale}-C_{\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{t}}\left(t\right)}{t+T_{gap}}$$

(1)

The formulation of the profit function involves several variables, namely, the cultivation time $t$, the production downtime between batches $T_{gap}$, the product formation at time t ($P\left(t\right)$), the sale price of the product ($P_{sale}$), and the cost of the batch up to time t ($C_{\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{t}}\left(t\right)$), encompassing various expenses such as raw materials, energy, labor, and others.

In the industrial production of 2-KGA, bioreactors are operated in parallel, and a large number of batches are conducted annually. In a workshop with multiple bioreactors for 2-KGA production, the primary goal is to maximize the profit-making potential of each batch and improve the overall profit of the workshop within the scheduling time span:

$$\mathrm{Overall} \mathrm{profit} =\sum _{i=1}^N J\left(t_{i,f}\right)\left(t_{i,f}+T_{gap}\right)$$

(2)

$N$ is the number of all the batches in the scheduling time span. When a fermentation batch is finished, $t_{i,f}$ and $J\left(t_{i,f}\right)$ are the termination time of the ith batch and the profit function at $t_{i,f}$, respectively.

As per the conventional approach, the 2-KGA fed-batch fermentation process involves uniform L-sorbose feeding for all batches. However, the proposed strategy recommends that the feeding amount of the substrate should be optimized on a batch-to-batch basis, considering the individual profit potential of each batch.

2.2. Assessing Batch Profitability Potential

During the online scheduling process, $J$($t$) should be calculated and predicted at every sampling point to evaluate the profit-making potential. However, to save labor costs, the pharmaceutical plant reduced the frequency of offline measurements in recent years, which are important for the evaluation of batch profit. For example, the L-sorbose concentration is typically determined only a few times towards the end of the cultivation to establish the batch endpoint due to limited sampling assays. Therefore, it is difficult to accurately calculate the profit function during the whole fermentation process under such conditions.

The feeding operation is basically performed in the middle stages of the whole fermentation process when the substrate concentration is below the standard value. The feeding time varies between different batches according to the operating condition at that time. The feeding time distribution of historical batches is shown in Figure 1.

To evaluate the profit function of an online batch before the feeding operation, one possible approach is to find batches in the historical database with similar time series features to the current batch. Then, the profit function of these batches can be utilized to evaluate the profit function of the current batch.

The lengths of the two time series that need to be compared for similarity cannot be equal. The similarity between different batches can be calculated by extending and shortening different time series.

Assuming the two time series as $X$ and $Y$:

$$\begin{array}{c}X=\left\{x_1,x_2,\cdots x_m\right\}\end{array}i=1,2,\cdots ,m$$

(3)

$$\begin{array}{c}\begin{array}{cc}& Y=\left\{y_1,y_2,\cdots y_n\right\}\end{array}\end{array}j=1,2,\cdots ,n$$

(4)

Figure 2 illustrates a simple instance, where $x_i$ and $y_j$ are univariate samples. The two series are morphologically very similar, but these morphological features (peaks and troughs) cannot be aligned in time. If the Euclidean distance-based method is used to calculate the similarity of the two series, it will not match intuitive understanding. However, if the matching is performed with scaling deformation on the series, the matching result will be enhanced. Figure 2 shows the matching in a way that allows the data to scale and deform on the time axis, as indicated by the yellow line.

In the time series matching, the data points in the sequences no longer correspond one-to-one, but instead have varying mappings of one-to-one, one-to-many, or many-to-one. This warping process involves minimizing the overall distance between the two sequences. The dynamic programming method is utilized to find the minimum distance between the series. The distance between X and Y is defined in Equation (5).

$$\mathrm{DTW} \mathrm{Distance} \left(X,Y\right)=D\left(m,n\right)$$

(5)

Then, the iterative process is as follows:

$$d_{ij}=d\left(x_i,y_j\right)$$

(6)

$$D\left(i,j\right)=d_{ij}+\min \left\{D\left(i-1,j\right),D\left(i,j-1\right),D\left(i-1,j-1\right)\right\}$$

(7)

$D\left(m,n\right)$ can be obtained by solving this dynamic programming problem. The best path of the dynamic warping process is shown in Figure 3. The depth of the color represents the distance of the corresponding sample points in the simple instance.

The measurements between the online batch and historical batches are performed with the sampled data before the feeding time. The key process variables are collected for batch comparison at each sampled point before that. Due to the unavailability of routine data, the L-sorbose consumption and biomass concentration are not used, even though they are significant parameters. Instead, the accumulated product formation and alkali consumption are utilized for time series matching, which can be obtained after data preprocessing from the available online measured and offline assayed data.

For an online batch X, the data pairs of all sampled points before the feeding time constitute the time series of the batch. The lth data pair at a sampled point $T_l$ is defined in Equation (8).

$$x_l={\left[T_{l}P\left(T_l\right) \mathrm{Alk} \left(T_l\right)\right]}^T$$

(8)

where $P\left(T_l\right)$ and $\mathrm{Alk} \left(T_l\right)$ are the accumulated product formation and alkali consumption at $T_l$. Since it is a multivariate time series and the fluctuation range of each variable is not consistent, each variable needs to be normalized before the time series matching calculation. Furthermore, the three process variables are assumed to be independent of each other. Therefore, the similarity distance between batches X and Y is defined with independent measures [35], as shown in Equation (9).

$$ind(M,X,Y,p)={\left(\sum _{d=1}^3 {\left|M\left(X_d,Y_d\right)\right|}^p\right)}^{1/p}$$

(9)

where $X_d$ (or $Y_d)$ represents the dth dimension in the multivariate time series of batch X (or Y). $M$ is the DTW metric for univariate time series $X_d$ and $Y_d$, and parameter $p$ is set to 2. The distance between the online batch X and all historical batches is calculated and sorted. Assume the first k series in the historical batches with the smallest distance from the current batch X are as follows:

$$Y_1,Y_2,\cdots Y_k$$

(10)

Additionally, the corresponding distances are as follows:

$$D_1,D_2,\cdots D_k$$

(11)

The profit function at the termination time of these batches is as follows:

$$J_1,J_2,\cdots J_k$$

(12)

To get the predicted profit function ${\widehat{J}}_{pre}$ before the feeding operation of the online batch, the intuitive idea is to simply take the average value of the nearest k historical batches. However, when the similarity distances of these k batches differ significantly, this approach weakens the effect of the closer batches on the results. Therefore, a weighted prediction is proposed with Equations (13) and (14).

$${\widehat{J}}_{pre}=\frac{w_1{J}_1+w_2{J}_2+\cdots +w_k{J}_k}{w_1+w_2+\cdots +w_k}=\frac{\sum _{i=1}^k w_i{J}_i}{\sum _{i=1}^k w_i}$$

(13)

$$w_i=\frac{1}{D_i}$$

(14)

This weighted prediction allows the batches with small distances to have a larger weight in the prediction and reduces the influence of batches with large fluctuations in the later stages of fermentation on the prediction results. The predicted profit function ${\widehat{J}}_{pre}$ is utilized for the evaluation of batch profit-making potential.

2.3. The Optimized Feeding Strategy

When industrial batches are terminated, the variance of the profit function is as much as 15%. The profit function profiles of 90 batches are illustrated in Figure 4 after normalization.

The confidence intervals shown in Figure 4 represent the range $J_m\pm \lambda \sigma$, where $\lambda$ is the coefficient associated with the 60% confidence limits, whereas $J_m$ and $\sigma$ represent the mean and standard deviation of the statistical data.

Based on a statistical analysis of the profit function, the fermentation batches can be classified into three categories: excellent, normal, and poor. For an online batch under consideration, the predicted profit function ${\widehat{J}}_{pre}$ is obtained before the feeding time. Then, the classification criterion is based on a comparison between ${\widehat{J}}_{pre}$ and the 60% confidence limit, as shown in Figure 4. The optimized strategy is designed based on the classification result of the online batch.

If the batch is classified into the excellent category, the feeding volume should be maximized to the limits of the bioreactor. The batch termination time ${\widehat{t}}_{i,f}$ after feeding operation is estimated with Equation (15) because fermentation will be terminated when the L-sorbose is almost consumed.

$${\widehat{t}}_{i,f}=\frac{C_{i,r}\left(T_{fed}\right)\cdot V_{i,c}+{\widehat{V}}_{i,fed}\cdot C_L}{L{C}_i}$$

(15)

where $C_{i,r}\left(T_{fed}\right)$ is the residual L-sorbose concentration at the feeding time $T_{fed}$; $V_{i,c}$ is the current volume of the ith batch; ${\widehat{V}}_{i,fed}$ is the feeding volume of L-sorbose; $C_L$ is the concentration of the feeding L-sorbose. $L{C}_i$ is the L-sorbose consumption rate of the ith batch during the latter stage of fermentation, which is approximately considered to be constant after the feeding operation.

The feeding volume of the ith batch ${\widehat{V}}_{i,fed}$ should be maximized to the limits of the bioreactor, which is determined as follows:

$$V_{i,c}+{\widehat{V}}_{i,fed}+\Delta {\widehat{V}}_{i,Alk}=V_{max}\cdot \beta$$

(16)

The maximum working volume of a bioreactor ($V_{max}$) is determined by its structural parameters, and the second feeding volume is set to 0.95 for safety. $\Delta {\widehat{V}}_{i,Alk}$ is the estimated volume of alkali solution fed during the remaining cultivation time. ${\widehat{V}}_{i,fed}$ and $\Delta {\widehat{V}}_{i,Alk}$ have a good linear correlation, and the total amount of alkali feeding is assumed to be equal to its total consumption, as it is automatically fed in the pH control loop.

For the batches being classified into the normal category, the total feeding volume is the same as the empirical scheduling, which assigns the same quantity of substrate to each batch. Other resources such as the alkali solution should be adjusted accordingly after the re-allocation of L-sorbose.

If a batch is classified into the poor category, no L-sorbose will be allotted in the feeding process, and the batch will be terminated sooner than originally planned. This is because a renewed batch is expected to generate more profit.

3. Results and Discussion

The evaluation of the proposed profit function prediction method is performed with the data of industrial 2-KGA fermentation. In Equation (13), the predicted profit function ${\widehat{J}}_{pre}$ is related to the choice of parameter k. When the time series of the current online batch is compared for similarity distance to all other batches in the historical database, the nearest k batches are utilized for the calculation of the predicted profit function using the proposed methods in Equations (11)–(14).

To find out the effect of the k value on the prediction results, a total of 600 historical batches were split into two parts: the 100 batches with the closest dates were used for model validation, whereas the remaining batches were used for the retrieval of similar batches. The model’s prediction error (root mean square error, RMSE) with different values of k (from 1 to 10) are compared. Figure 5 presents the mean and standard deviation of prediction error.

When k = 6, the minimum mean of prediction error is achieved. It is noteworthy that despite the marginally elevated standard deviation observed with k = 6, this value is still deemed optimal, as the foremost objective of the feeding strategy is to optimize the predictive accuracy of the profit function. Therefore, the application of k = 6 is chosen to minimize the prediction error, and then the predicted profit function ${\widehat{J}}_{pre}$ of an online batch can be obtained before the feeding operation.

For model testing, the proposed method is applied to 90 industrial fermentation batches, and the process data are collected. The prediction error distribution with different k values for these batches is compared in Figure 6. With the choice of k = 6, the mean and standard deviation of the prediction error for these batches are 1.91% and 2.06%, respectively.

To validate the proposed approach’s efficacy, other feasible methods have also been applied to predict the profit function, with the same limitation that only the process data before the feeding operation can be utilized. The Artificial Neural Network (ANN) [36,37], Support Vector Machine (SVM) [38,39], Random Forest (RF) [40,41], and eXtreme Gradient Boosting (XGBoost) [42,43,44] have been shown to achieve state-of-the-art results in many applications with sequential data or time series.

The moving data window technique is utilized to generate input–output data pairs for the models [45]. The data pair corresponding to the time stamp of feeding operation Tv is presented with Equations (17)–(19).

$$X={\left[T_v{x}_{T_v}{x}_{T_{v-\tau }}{x}_{T_{v-2\tau }\cdots X_{T_{v-T_D}}}\right]}^T$$

(17)

$$\begin{array}{c}{x}_{T_v}=\left[P\left(T_v\right)Alk\right(T_v){]}^T\end{array}$$

(18)

$$Y=J\left(t_{i,f}\right)$$

(19)

where $\tau$ is the time interval of the variables for discretization,$t_{i,f}$ and $J\left(t_{i,f}\right)$ are the termination time of the ith batch and the profit function at $t_{i,f}$, respectively. The input variable consists of product formation (P) and alkali consumption (Alk). Since the feeding time varies between different batches, $T_D$ is set to 20, so that the width of the input data window for each batch is consistent.

The data pairs generated from the 600 historical batches are utilized for the training and validation of the other mentioned models. The testing of these models is performed with the same 90 industrial batches with the proposed multidimensional time series aggregation method (MTSA, k = 6). Specifically, a classical three-layer back propagation network is employed in the ANN model, which encompasses three neurons in the hidden layer. In the SVM model, a three-degree polynomial kernel is utilized, and the regularization parameter C is set to 1.0. Furthermore, the random forest (RF) regressor’s parameter min_samples_split is configured to a value of 2, whereas the n_estimators is set to 300. Finally, in the XGBoost model, a learning rate of 0.1 is utilized, with the number of estimators established at 300 and the maximum depth of 5. It is worth mentioning that these parameters are manually adjusted to ensure a good performance of each model.

The average prediction error with different methods is compared in

Table 1. Mean and standard deviation of prediction error using different methods.

Model	ANN	SVM	RF	XGBoost	MTSA
Mean (%)	3.61	3.55	3.59	2.78	1.91
Standard deviation (%)	2.35	2.33	2.37	2.19	2.06

.

Table 1 demonstrates that the MTSA method achieves better results for these batches compared to the other mentioned models. According to the proposed feeding strategy, the online batches being classified into the excellent category will be fed with more resources.

The effectiveness of the entire strategy is evaluated using pseudo-online simulation, which is carried out with the process data of 90 industrial batches. The inoculation sequences of these batches are assumed to be the same with empirical scheduling. For online scheduling, after the strategy is practiced, the actual process data can be obtained through measurement. However, for pseudo-online simulation, the process data after the optimized scheduling is not available. Therefore, the termination time and profit function of the corresponding batch are evaluated using Equations (13) and (15), and then the overall profit of the workshop in the scheduling time span can be obtained with Equation (2).

The scheduling results of all batches are presented in Figure 7, revealing the comparison between the optimized scheduling and the empirical ones. Arrows denote the normalized $J\left(t_{i,f}\right)$ changes from empirical scheduling results to the corresponding optimized scheduling ones. The optimized scheduling strategy prolongs the cultivation duration of batches classified as excellent with high-profit potential and shortens the cultivation duration of the poor batches.

Table 2. Comparative analysis of the scheduling strategies employed in the cultivation of 90 industrial batches: empirical versus optimized.

	Empirical Scheduling	Optimized Scheduling
Total fermentation time (h)	5589	5675
Total profit (PU)	13,970,000	15,045,051

shows the comparison of the total cultivation time and the profit of 90 industrial batches under optimized scheduling with empirical scheduling. It is assumed that these batches are being operated in series within a single bioreactor to obtain the total cultivation time for all batches. The optimized scheduling resulted in a cultivation time extension of 86 h compared to the empirical scheduling (5675–5589), with a total profit of 15,045,051 PU (Price Unit). The cumulative profit obtained from the 90 batches under empirical scheduling amounts to 13,970,000 PU.

To facilitate comparison, it is postulated that the overall cultivation time for the empirical scheduling is extended to equal that of the optimized scheduling, with an overall profit of (13,970,000 + 86 × J_ave = 14,185,000) PU. J_ave (2500 PU/h) is the average per-batch profit of the 90 industrial batches under empirical scheduling. Then, a total profit increase of 6.06% is observed through the optimized resource allocation.

4. Conclusions

In this work, an online feeding strategy for 2-KGA fed-batch fermentation is proposed, which makes the online economic performance evaluation and optimization practicable for industrial batches. Reallocation of the substrate L-sorbose, along with the alkali solution and other relevant resources, is executed to increase the overall profitability of the workshop. The decision for optimization is based on the predicted profit function, which indicates the potential for each individual batch to generate profit. However, the profit function is not always available because of the lack of offline measurements during cultivation. Consequently, the profit-making potential is evaluated using the predicted profit function at the termination time and the online batch classification result. The pseudo-online scheduling results of the 90 industrial batches indicate that the total profit under optimized scheduling increases by more than 6% compared to empirical scheduling.