Biogas Prediction Enhancement of a Swine Farm Bio-Digester Using a Lag-Based Surrogate Machine Learning Model

1. Introduction

The biogas generation process is a natural phenomenon that occurs in a variety of anaerobic environments, [1]. These include marine and freshwater environments, sediments, wastewater treatment plant sludge, etc. Interest in the process stems primarily from the following reasons. On the one hand, a high degree of organic matter reduction is achieved with a small increase—compared to aerobic processes—in bacterial biomass. On the other hand, biogas production can be used to generate various forms of energy (heat and electricity) or processed as automotive fuel. Nevertheless, biogas has a lower calorific value than natural gas, and in specific applications, such as automotive fuel, treatment is necessary to improve its quality, [2].

Energy production from the anaerobic digestion has been used worldwide for over 30 years, [3]. Its viability and profitability depend not only on the amount of biogas produced, the available technology, and the efficiency of the wastewater treatment operation, but also on external parameters such as the local cost of energy production and available energy resources, [4]. Besides the economic advantages, biogas has also yielded environmental benefits, [5]. Anaerobic digestion represents a mature yet still evolving technology for renewable energy production. Predicting methane yield remains challenging due strong nonlinearities and environmental coupling. Anaerobic digestion is a biological process where organic carbon is converted through subsequent oxidation and reduction to its more oxidized state ($C{O}_2$) and its more reduced state ($C{H}_4$). A wide range of microorganisms catalyzes the process in the absence of oxygen. Nitrogen, hydrogen, ammonia, and hydrogen sulfide are also generated in smaller quantities (typically less than 1% of the total gas volume). The mixture of gaseous products is called biogas, and the anaerobic degradation process is often also referred to as the biogas process, [6]. Estimating biogas production in an anaerobic digester is a complicated problem because the system is not only chemical nor only biological, but a microbial ecosystem highly related to biochemical reactions, mass transfer and variable physicochemical conditions, [7]. Current society requirements demands high efficiency biogas production to be implemented in anaerobic digesters of wastewater treatment plants, [8]. Therefore, innovative solutions must be developed in the short term, [9].

Recent advances in machine learning, particularly deep learning architectures, have shown promising results in modeling temporal dynamics [10]. Nevertheless, some open problems are still challenging. In [11], the optimization of biogas production for treating domestic wastewater is studied using different machine learning models (XGBoost and PSO), pointing out that limitations on the volume of the training data influenced the performance of the model predictions. J. Schulz et. al. [12] investigated the features of Carbon to Nitrogen ratio of the substrates that can be used for long-term continuous Anaerobic co-digestion. In general, the study of sensors calibration and fouling, modeling and optimization approaches together with the lack of high quality standardized datasets are current trends in the literature, [13]. Some solutions using frontier technologies such as digital twins, [14] might help to shorten the implementation curve provided AI-based models are accurate and simple enough to be implemented in real-time biogas production processes.

Herein we aim to improve biogas estimation using lag-based training vectors. To avoid the lack of accurate and enough datasets, in this work, we propose a biogas estimation surrogate model based on machine learning. Unlike prior studies, our approach explicitly addresses data scarcity through synthetic data generation, temporal dependencies via lag-based modeling, as well as mechanistic consistency through ODE constraints. The model was trained using a set of data obtained from the solution of a set of differential equations with experimentally validated parameters. Our analysis indicates that biogas estimation can be achieved with more than 97% of accuracy. The rest of the work is structured as follows. Section 2 describes the proposed methodology while Section 3 depicts the corresponding results to validate our approach. Section 4 presents a brief discussion and, finally Section 5 closes the paper with conclusions and future work.

2. Methodology

2.1. Experimental Design

To systematically evaluate the predictive capacity of the proposed surrogate model for the swine farm bio-digester, in this case methane is the produced biogas, the methodological framework was structured into two progressive experiments:

• Base model: 30,000 registers were generated using different stoichiometric and kinetic parameters based on the real bio-digester by means of the ODEs (1), (2), (3), (4) and features in terms of inflow and the organic substrate concentration $S_1$ to produce biogas calculated through the equation (5). An Extreme Gradient (XGBoost) model was trained and tested using a sequential 80/20 split, respectively.
• Lag-based improved model: 10,000 registers were generated incorporating a dynamic seasonal temperature factor to perturb the microbial kinetic rates and simulate real-world environmental exposure. To capture this physical complexity, integrating thermodynamic interaction variables and historical lags of the biogas production. The inclusion of this auto-regressive feature is valid given that this information is available in real world applications, it is calculated through equation (5). An XGBoost model was trained and tested using a sequential 80/20 split to evaluate the system’s inertial memory and prevent data leakage.

2.2. Mathematical Modeling and Mass Balance

The physical and biochemical dynamics of the anaerobic digestion process (acidogenic and methanogenic phases) were simulated using a system of ODEs (1), (2), (3), (4). Given that the real rural bio-digester operates with a plug-flow hydrodynamic regime, the total reactor volume ($V_{total}$) was spatially discretized into a series of Continuous Stirred-Tank Reactors (CSTRs). For this study, the system was divided into N = 3 interconnected sub-reactors.

The general mass balance for any given state variable within the i-th sub-reactor is consistent with the fundamental principle of accumulation (Accumulation = Inflow - Outflow + Net Reaction) which is described by the following system of ODEs:

$${\displaystyle \frac{d{X}_{1,i}}{dt}}=D_{sub}(X_{1,i-1}-X_{1,i})+{\mu }_{1,i}{X}_{1,i}$$

(1)

$${\displaystyle \frac{d{X}_{2,i}}{dt}}=D_{sub}(X_{2,i-1}-X_{2,i})+{\mu }_{2,i}{X}_{2,i}$$

(2)

$${\displaystyle \frac{d{S}_{1,i}}{dt}}=D_{sub}(S_{1,i-1}-S_{1,i})-K_1{\mu }_{1,i}{X}_{1,i}$$

(3)

$${\displaystyle \frac{d{S}_{2,i}}{dt}}=D_{sub}(S_{2,i-1}-S_{2,i})+K_2{\mu }_{1,i}{X}_{1,i}-K_3{\mu }_{2,i}{X}_{2,i}$$

(4)

where $X_{1,i}$ and $X_{2,i}$ for $i=1,2,...,N$ denote the acidogenic and methanogenic biomass concentrations, respectively, for the $i-th$ biodigester. $S_1$ represents the organic substrate concentration measured as Chemical Oxygen Demand (COD), and $S_2$ defines the Volatile Fatty Acids (VFA) concentration. $D_{sub}$ is the local dilution rate for each sub-reactor, calculated as $D\times N$ (The global dilution D is defined as the flow rate $Q_{in}$ divided by the total volume $V_{total}$). This formulation implies a reduction in the effective residence time within each sub-reactor, consistent with the representation of reactors in series.

For the first sub-reactor (i=1), the feed concentrations correspond to the raw inputs ($S_{1,0}$=$S_{1in}$ and $S_{2,0}$=$S_{2in}$), assuming a biomass-free feed stream ($X_{1,0}$ and $X_{2,0}=0$).

Methane production is calculated using equation (5)

$$C{H}_4=\sum _{i=1}^N{K}_4{\mu }_{2,i}{X}_{2,i}{V}_{sub}$$

(5)

2.3. Microbial Kinetics and Thermodynamic Perturbation

Microbial growth kinetics were modeled using the Monod equation for the acidogenic phase (6) and the Haldane inhibition model for methanogenesis (7). To improve the performance of the surrogate model, two scenarios were formulated. The base model analysis considered standard kinetic rates without environmental perturbations as expressed in equation (6) and equation (7).

$${\mu }_{1,i}={\mu }_{max1}\left({\displaystyle \frac{S_{1,i}}{K_{s1}+S_{1,i}}}\right)$$

(6)

$${\mu }_{2,i}={\mu }_{max2}\left({\displaystyle \frac{S_{2,i}}{K_{s2}+S_{2,i}+\frac{S_{2,i}^2}{K_I}}}\right)$$

(7)

where ${\mu }_{max}$ represents the maximum specific growth rates under standard conditions, $K_{si}$ for $i=\{1,2\}$ denotes the half-saturation constants, and $K_I$ is the inhibition constant due to Volatile Fatty Acids (VFA) accumulation.

The lag-based improved model incorporated a seasonal temperature factor to simulate real-world environmental exposure and its direct impact on bacterial growth rates. It is important to note that, although the model is simplified, this representation captures the dominant seasonal dynamics while providing a controlled and physically consistent perturbation framework for training and evaluating the surrogate model. The modified kinetic equations are shown in equation (8) and equation (9):

$${\mu }_{1,i}=({\mu }_{max1}·T_{factor})\left({\displaystyle \frac{S_{1,i}}{K_{s1}+S_{1,i}}}\right)$$

(8)

$${\mu }_{2,i}=({\mu }_{max2}·T_{factor})\left({\displaystyle \frac{S_{2,i}}{K_{s2}+S_{2,i}+\frac{S_{2,i}^2}{K_I}}}\right)$$

(9)

To model the annual environmental exposure of rural bio-digesters, the $T_{factor}$ was defined by a cyclical model over a 365-day period as shown in equation (10):

$$T_{factor}\left(t\right)=1.0+0.15sin\left({\displaystyle \frac{2\pi t}{365}}\right)$$

(10)

Finally, the total daily methane production ($C{H}_4$), representing the system’s target energy yield, was calculated using equation (5) where $V_{sub}$ is the volume of each individual sub-reactor ($V_{total}/N$), and $K_4$ is the calibrated stoichiometric yield coefficient for methane.

2.4. Synthetic Data Generation

To train the XGBoost model, extended operational periods were simulated by numerically solving the ODE’s system using the scipy.integrate.odeint library in Python. To reflect the real operation, information from the records of the real bio-digester were considered using bounded normal distributions (empirical operation ranges), the numerical integration required defining an initial condition vector ($t=0$) for each sub-reactor, representing the starting biomass and substrate concentrations inside the digester. Across all sub-reactors, the initial state was set to: $X_1=1.8$ g/L, $X_2=0.8$ g/L, $S_1=1500.0$ mg/L, and $S_2=10.0$ mmol/L.

• Input flow rate ($Q_{in}$): Modeled as $\mathcal{N}(4.5,1.0)$ and physically constrained to the range $[1.0,8.0]$ m³/d.
• Organic substrate loading ($S_{1in}$): Modeled as $\mathcal{N}(2500,500)$ and physically constrained to the range $[1000,4000]$ mg/L.

Through this computational approach, two distinct datasets were generated and exported for the machine learning pipeline:

• Base model: A dataset with 30,000 records with the exported feature space consisted of Simulation Time (Time), Input Flow Rate (Q_in_m3_d), Organic Substrate Loading(S1_in_mg_L), and the target variable Methane Biogas Production (CH4_Prod_L_d).
• Lag-based improved model: A dataset with 10,000 records with the exported feature space included the thermodynamic perturbation, consisting of Simulation Time (Time), Input Flow Rate (Q_in_m3_d), Organic Substrate Loading (S1_in_mg_L), the Seasonal Temperature Factor (Temp_Factor), and the target variable methane biogas Production (CH4_Prod_L_d).

2.5. Feature Engineering and Machine Learning Architecture

To ensure model stability and eliminate the influence of the initial unsteady-state mathematical start-up, a transient period was discarded from the synthetic datasets prior to the feature engineering phase (the first 50 days for the base model and the first 100 days for the lag-based improved model). To capture the system’s temporal dynamics and provide the predictive models with historical memory, specific systematic lags vectors were engineered consistent with the available information in operational bio-digester. Lag intervals of $\tau =\{1,5,10,15,20,25,30\}$ days and $\tau =\{1,2,3,5,10,15,20\}$ days were incorporated into the feature space for the base model and lag-based improved model, respectively.

For the base model scenario, the input feature vector ${\mathbf{x}}_t$ constructed for any given day t resulted in a 16-dimensional array, organized as shown in equation (11):

$${\mathbf{x}}_t={\left[S_{1in}\left(t\right),Q_{in}\left(t\right),S_{1in}(t-{\tau }_1),\dots ,S_{1in}(t-{\tau }_7),Q_{in}(t-{\tau }_1),\dots ,Q_{in}(t-{\tau }_7)\right]}^T$$

(11)

where the lags set is defined as $\tau \in \{1,5,10,15,20,25,30\}$ days.

However, for the lag-based improved model, the $T_{factor}$ and its physical interactions with the mass flows were considered. Additionally, the records of historical biogas production ($C{H}_4$) were considered. The inclusion of lagged methane production variables is consistent with real-world operational scenarios, where historical biogas measurements are continuously monitored and readily available. In this context, these variables provide valuable temporal information that enhances short-term predictive performance. Therefore, this information allows the model to exploit inherent temporal dependencies of the process, enhancing predictive capability within a realistic digital twin framework.

The thermodynamic interaction variables were defined as $I_{S1}\left(t\right)=T_{factor}\left(t\right)·S_{1in}\left(t\right)$ and $I_Q\left(t\right)=T_{factor}\left(t\right)·Q_{in}\left(t\right)$. Based on the refined lags set $\tau \in \{1,2,3,5,10,15,20\}$ days, the augmented input feature vector ${\mathbf{x}}_t$ expanded into a 40-dimensional array, organized as shown in equation (12):

$$\begin{array}{cc}\hfill {\mathbf{x}}_t=[& S_{1in}\left(t\right),Q_{in}\left(t\right),T_{factor}\left(t\right),I_{S1}\left(t\right),I_Q\left(t\right),\hfill \\ \hfill & S_{1in}(t-{\tau }_1),\dots ,S_{1in}(t-{\tau }_7),\hfill \\ \hfill & Q_{in}(t-{\tau }_1),\dots ,Q_{in}(t-{\tau }_7),\hfill \\ \hfill & T_{factor}(t-{\tau }_1),\dots ,T_{factor}(t-{\tau }_7),\hfill \\ \hfill & I_{S1}(t-{\tau }_1),\dots ,I_{S1}(t-{\tau }_7),\hfill \\ \hfill & C{H}_4(t-{\tau }_1),\dots ,C{H}_4(t-{\tau }_7){]}^T\hfill \end{array}$$

(12)

The XGBoost model was selected because of its robustness against overfitting and its capacity to handle non-linear interactions without explicit scaling.

2.6. Model Training and Evaluation Metrics

To prevent data leakage and rigorously evaluate predictive performance, a sequential chronological split was implemented for both experiments. The first $80\%$ of the temporal registers were used for model training, while the remaining $20\%$ were reserved as an unseen testing set.

The hyperparameters of the XGBoost regressor were tuned for each experiment to balance learning capacity and generalization:

• Base model configuration: The model used 1000 boosting rounds (estimators), a learning rate of $0.05$, and a maximum tree depth of 5. Subsampling and column sampling by tree were both set to $0.80$. Early stopping was triggered if the validation loss did not improve for 50 consecutive rounds.
• Lag-based improved model configuration: The model used 5000 estimators with a learning rate of $0.01$ and a maximum tree depth of 12. Subsampling and column samplings were set to $0.70$, with an early stopping patience of 100 rounds.

The performance of the model was assessed using the coefficient of determination ($R^2$), the Root Mean Square Error (RMSE), and the Mean Absolute Error (MAE), which are defined in equations (13), (14), and (15), respectively:

$$R^2=1-{\displaystyle \frac{{\sum }_{i=1}^n{(y_i-{\widehat{y}}_i)}^2}{{\sum }_{i=1}^n{(y_i-\overline{y})}^2}}$$

(13)

$$RMSE=\sqrt{{\displaystyle \frac{1}{n}}\sum _{i=1}^n{(y_i-{\widehat{y}}_i)}^2}$$

(14)

$$MAE={\displaystyle \frac{1}{n}}\sum _{i=1}^n|y_i-{\widehat{y}}_i|$$

(15)

where n represents the total number of observations in the testing set, $y_i$ is the ground-truth biogas production generated by the phenomenological model, ${\widehat{y}}_i$ is the predicted biogas production by the XGBoost model, and $\overline{y}$ is the mean of the actual observed values.