Optimization and Analysis of Liquid Anaerobic Co-Digestion of Agro-Industrial Wastes via Mixture Design

Anaerobic co-digestion (AcoD) is a widely employed technique to produce biogas from simultaneous digestion of various biomasses. However, the selection of the optimal proportions of the substrates in the mixtures presents a challenge. This research used a mixture design to investigate the interactions between the liquid fraction of piglet manure (PM), cow manure (CWM), and starch wastewater (SWW). A modified Gompertz model was used to identify the statistically significant parameters of the methane production curves. The optimal compositions of the mixtures were identified based on multi-objective optimization of the maximal methane yield (YCH4) and maximal methane specific production rate (rCH4) parameters. The study was validated using a double mixture of PM and CWM and a triple mixture. The estimated degradation rates for both mixtures were faster than the predicted ones. The absolute relative errors of rCH4 were 27.41% for the double mixture and 5.59% for the triple mixture, while the relative errors of YCH4 were 4.64% for the double mixture and 10.05% for the triple mixture. These relative errors are within the normal limits of a process with high variability like AD. Thus, mixture design supported by the tested models is suitable for the definition of practically advisable mixtures of substrates.


Introduction
Enhancement of anaerobic digestion (AD) is a commonly researched topic. The improvement can be achieved by several methods. Among the most popular and recommended methods for improving AD one can find pretreatments or enzymes (i.e., biocatalysts), reactor engineering, coupling AD with dark fermentation, genetically improving the microbial community (bio-augmentation), and anaerobic co-digestion (AcoD) [1].
AcoD involves the simultaneous AD of two or more substrates. It has proved to be a viable option for improving biogas production because it alleviates the disadvantages of mono-digestion while increasing the economic feasibility of the process [2]. Distinct advantages of AcoD include the supply of macro and micronutrients, balanced carbonnitrogen ratio, superior buffer capacity, dilution of inhibitors, and potentially enhanced biogas production [3].
There exist different criteria to assess the performance of AD processes, however, the most well-accepted and commonly used is the bio-methane potential (BMP) procedure. BMP is defined as the capacity of a substrate to be converted into methane and carbon dioxide. Determination of BMP is the first step in evaluating the digestibility or applicability of a substrate. The BMP parameter provides valuable information about general degradability, expectable energy yield, and the economic evaluation of new biogas plants. BMP is usually determined by employing a BMP test procedure in a batch anaerobic fermentation assay. This method is reliable, straightforward, and avoids inconsistencies in collected data [4,5]. The BMP test consists of adding a known quantity of an organic substrate to an active anaerobic inoculum in an air-tight serum bottle [6,7].
BMP optimization is the first step in AcoD optimization. BMP optimization for a substrate mixture is usually conducted based on statistical methods and variables. Substrate ratios and the inoculum/substrate ratio are the experimental factors. However, additional process variables may be included [4,8].
A popular statistical optimization method is response surface methodology (RSM). RSM is a group of mathematical and statistical techniques that identify process improvements based on the fit of empirical models to measured experimental data [9]. A special type of RSM is the mixture design. Mixture design is an effective method for determining the optimal proportions of ingredients in a mixture [10].
In mixture design experiments, the independent variables are the proportions of the investigated components. The total amount of material must be held constant in a mixture design experiment. This allows for analysis of the dependency of the component proportions without confounding variability due to changes in the total amount of the mixture. The mixture design analysis provides valuable information about the interactions between independent factors. It also provides a better understanding of the response variables [11].
Mixture design has been previously used to understand how substrates interact during AcoD. Pagés-Díaz et al. [12] used a four-factor mixture simplex-centroid design, which employed solid cattle slaughterhouse wastes, manure, various crops, and municipal solid wastes. Methane yield (Y CH4 ) and specific methane production rate (r CH4 ) were the response variables. Rahman et al. [13] used two sets of mixtures. The first set consisted of poultry droppings with sugarcane bagasse. The second set consisted of press mud and poultry droppings with roots, tops, and press mud from sugar beets. An augmented simplex-centroid design was applied to describe the interactions. This paper generated experimental data for the integration of high-rate anaerobic reactors in conventional agricultural biogas plants. The objective was to find the optimal mixture composition of three substrates that are usually found in great quantity in the agro-industrial sector, especially in the region of North Rhine-Westphalia (NRW), Germany. A plant's profitability and flexibility can be improved by using information about the interactions between these substrates. The substrates used in this study were piglet manure (PM), cow manure (CWM), and starch wastewater (SWW). A three-factor mixture design was employed to analyze and describe the interactions between the substrate BMP tests. The BMP tests were only conducted using the liquid fraction of the selected substrates. This was necessary because the study used expanded granular sludge bed (EGSB) reactors which operated in a continuous mode with a dry matter (DM) content of less than 8 wt%.

Raw Materials
In this study, the three selected substrates were PM, CWM, SWW. In NRW, animal manure management is a topic of interest due to the large livestock population [14]. PM and CWM were selected to address this issue. Both PM and CWM provide a strong buffer capacity that supports pH balance in the livestock management operation. SWW was selected due to its low DM content. Low DM content avoids clogging in the reactors of a long-term operation. Also, the high content of readily degradable carbohydrates in SWW may increase the biogas yield and the co-digestion degradation rate of the mixture [15]. In Germany, SWW is highly available because it is obtained from crops that have already been treated by AD, including potatoes, corn, and wheat [16].
Previously, PM and CWM were processed using screw press systems that separated the solid and liquid phases. For PM, a dual screw press (i.e., Vakusep from BETEBE GmbH) with a sieve filter size of 100 µm was applied. CWM was collected from farmers in a pre-separated form which was then processed by a second separation using o a screw press with a sieve size of 200 µm.
The liquid fractions of the substrates were used to characterize the content of DM, volatile solids (VS), macromolecules, and nutrients. The results of these analyses are shown in Table 1.

Batch-Test Setup
Batch assays were conducted on fermentation tests of the organic substances, based on the VDI 4630 guidelines. Each test was performed using a 1000 mL glass vessel. Five grams of VS were added to the vessel of each test. A defined amount of substrate and a previously calculated quantity of inoculum were weighed into the reaction vessel. When necessary, the reaction vessel was carefully filled with warm water to reach a reaction volume of 800 mL. To prevent inhibition, the substrate inoculum rate VS substrate /VS inoculum was kept constant at 0.5. Each bottle was vigorously agitated, sealed with a rubber stopper, and clamped down with a plastic screw cap connected to an eudiometer tube (as shown in Figure 1). The produced gases were transferred through PVC hoses into 1000 mL eudiometer tubes. This process enabled daily measurements of gas volume and quality. The eudiometers were sealed with a barrier fluid comprised of water and 5 wt% sulfuric acid. The acidification of the barrier liquid prevented carbon dioxide (CO 2 ) from dissolving in the gas mixture. Additionally, 7.5 wt% sodium sulfate (Na 2 SO 4 ) was added to prevent the entry of CO 2 into the barrier fluid [17]. The batch assays were developed for 42 days.
A test vessel with cellulose as the reference substrate was produced to compare and ensure adequate biological activity by the inoculum. The biogas potential of cellulose is known. Thus, cellulose can be used as a reference for evaluating the reliability of an experiment. Also, an inoculum-only batch fermentation test or zero-test was conducted. In the zero-test, the gas production value from the inoculum was subtracted from the gas production value from the substrate inoculum. Each fermentation test was performed at least three times, including the cellulose reference sample and the zero-test sample [17,18]. The pH was measured at the beginning and the end of every batch-test because changes in pH are usually correlated with other operational parameters. The accumulation of organic acids (acidification) typically lowers the pH, while increased ammonia concentrations or CO 2 removal increases the pH [19,20].

Mathematical Modeling
Prior to modeling the methane yield curves (MYC), pre-processing of the measured data was necessary. The pre-processing of the measured data involved the application of three criteria: (1) Stopping criterion: The test was concluded when the relative increase of YCH4 was less than 1% for three consecutive days. (2) Plausibility criterion: The existence of abrupt or non-monotonic trends in the curves requires individual analysis of the affected test. (3) Reproducibility/accuracy criterion: After deleting possible outliers, a coefficient of variation (CV) smaller than 5% between the curves was required.
If any of the above criteria were not met, the sample from all three batch-tests was eliminated from the study. The three criteria are based on the recommendations by Holliger et al. [18].
For validation, BMP results from the experimental tests were compared to reported

Mathematical Modeling
Prior to modeling the methane yield curves (MYC), pre-processing of the measured data was necessary. The pre-processing of the measured data involved the application of three criteria: (1) Stopping criterion: The test was concluded when the relative increase of Y CH4 was less than 1% for three consecutive days. (2) Plausibility criterion: The existence of abrupt or non-monotonic trends in the curves requires individual analysis of the affected test. (3) Reproducibility/accuracy criterion: After deleting possible outliers, a coefficient of variation (CV) smaller than 5% between the curves was required.
If any of the above criteria were not met, the sample from all three batch-tests was eliminated from the study. The three criteria are based on the recommendations by Holliger et al. [18].
For validation, BMP results from the experimental tests were compared to reported BMP values from the literature and calculated theoretical BMP values. The theoretical BMP values and the degradation fraction (f d ) values were calculated using the equations from Raposo et al.,Ebner et al. [5,21], shown in Equations (1) and (2), respectively. (1) where Xi: fraction of the macromolecule expressed in g of macromolecule per g of total so-lids BMP theol : theoretical BMP calculated using Equation (1) (L CH4 /kg VS ) BMP measured : measured BMP obtained as a result of the practical tests (L CH4 /kg VS ) The most appropriate models to simulate MYCs include quantification of important parameters like maximal Y CH4 , r CH4 , and, if it exists, lag phase time (λ). The most commonly used of these models are the first-order one-step model [6], the first-order two-step model [22] and the modified Gompertz model [12]. All three models were fitted to the measured data. Goodness-of-fit statistics were compared to identify the model with the best fit. The non-linear regression tool of the software Minitab 19 was used to fit the models to the measured data. The root mean square error (RMSE) statistic was used to evaluate goodness-of-fit. An RMSE value less than 10 L CH4 /kg VS was designated as an acceptable goodness-of-fit. The first-order one-step model, the first-order two-step model, and the modified Gompertz models are shown in Equation (3), Equation (4), and Equation (5), respectively ( Table 2). Table 2. Models employed to describe the methane yield curves.

Model Equation Reference
First-order one-step First-order two-step Modified Gompertz BMP ∞ : Extrapolated BMP at infinite retention time in L CH4 /kg VS . k: First-order reaction constant (1/d). k hyd : First-order reaction constant of the first step (hydrolysis/acidification) (1/d). k vfa : First-order reaction constant of the second step (volatile fatty acids degradation) (1/d). R max : Maximum biomethane production rate (L CH4 /kg VS /d). λ: Lag time (d). t: Time (d).

Mixtures Characterization
The synergistic and antagonistic effects of the individual mixtures on the BMP and reaction rate were characterized by two indices: (1) the co-digestion index (CI) (Equation (7)) and (2) the kinetic index (KI) (Equation (8)). In these two equations, the parameter values for the BMP and reaction rates were estimated from the three selected models that describe the MYCs. The CI (%) expresses the relative increase in YCH4 compared to the sum of the individual substrates in the mixture. Thus, it detects if the interaction between the substrates is positive or negative. The positive interactions are interpreted as synergistic effects. However, this term should be used with caution since batch-tests are inoculum-biased. Instead, the term "acute effects" is suggested to describe the positive interactions in batchtests [17,23].
The KI (%) expresses the relative increase in the degradation rate compared to the fastest substrate in the mixture. A positive KI (%) is likely related to an improvement in continuous operation mode [22].

Experimental Design
A three-factor simplex-centroid mixture design with seven design points was used to evaluate the interaction between the substrates for the response variables. Based on the procedure described in Section 2.2, the seven design points were replicated three times for a total of 21 data points. The response variables were Y CH4 and r CH4 . The mixture design's experimental points employed in this study are detailed in Table 3. The compositions of the mixtures are expressed in VS ratios (% VS) since the total VS in each batch-test was held constant at 5 g to avoid confounding variability from the response variables. Based on Equation (9), each component's proportion varied between 0 and 1, and the variable x i represents the proportion of ith constituent in the mixture.
The effect of the mixtures on the response variables was modeled by a special cubic model shown in Equation (10) [10,11,13]: The β ij and β ijk coefficient values indicate the strength of the interaction between the substrates. The sign of the β ij and β ijk coefficients indicates whether the interaction is positive or negative. An ANOVA analysis was used to identify the terms in the model. In an ANOVA analysis, the p-value associated with the statistical confidence level determines if a term should be included in the model [24]. The response optimization method was employed to identify the combination of substrate proportions that simultaneously optimize both responses. A desirability function was applied in the optimization procedure. Desirability is an objective function that ranges from zero outside of the limits to one at Processes 2021, 9, 877 7 of 16 the goal. The desirability function approach is one of the most widely used methods for optimizing multiple response processes. This method identifies the operating conditions that provide the most desirable response values given the specified assumptions. The characteristics of a goal may be altered by adjusting the weight or importance of each variable or the ranges within the optimization performance [24].

Analysis of the Curves
The MYCs were predominantly smooth with a slight leaning toward the logistic growth and signs of small lag phases. Subsequently, the dispersion analysis of the batches identified an outlier in one of the three PM replicate runs. This run was eliminated, to keep a CV of less than 5%. Furthermore, to make symmetrical the mixture experimental design the outlier curve was replaced during the mixture design analysis by the average of the two remaining curves. The rest of the CVs were approximately ≤5%, which is the validation criteria recommended by Holliger et al. [18]. The curves of the single digestions are presented in Figure 2, the numbers in the legend represent the valid number of replicates of each curve.
Processes 2021, 9, x FOR PEER REVIEW 7 of 17 of the limits to one at the goal. The desirability function approach is one of the most widely used methods for optimizing multiple response processes. This method identifies the operating conditions that provide the most desirable response values given the specified assumptions. The characteristics of a goal may be altered by adjusting the weight or importance of each variable or the ranges within the optimization performance [24].

Analysis of the Curves
The MYCs were predominantly smooth with a slight leaning toward the logistic growth and signs of small lag phases. Subsequently, the dispersion analysis of the batches identified an outlier in one of the three PM replicate runs. This run was eliminated, to keep a CV of less than 5%. Furthermore, to make symmetrical the mixture experimental design the outlier curve was replaced during the mixture design analysis by the average of the two remaining curves. The rest of the CVs were approximately ≤5%, which is the validation criteria recommended by Holliger et al. [18]. The curves of the single digestions are presented in Figure 2, the numbers in the legend represent the valid number of replicates of each curve. Both PM and CWM had reasonably small lag phases. This was attributed to the relatively large presence of nitrogen-associated compounds, which are common in these substrates [19]. Nevertheless, the ammonia concentrations were very different from the inhibition values reported by [25]. Thus, the lag phases were probably related to an adaptation phase rather than due to inhibition. SWW was the only readily degradable substrate because of the relatively large presence of nitrogen-free extracts. Consequently, no lag phase was detected for SWW. Instead, a sharp change in slope occurred after the second day, indicating slight diauxic behavior. Thus, it was assumed that a second substrate was consumed after the second day. However, no plateau phase was observed, and the MYC was considered monotonic.
The majority of the batch-tests had small lag phases. However, none of the lag phases were long enough to cause inhibition. Similarly, the YCH4 values were compared to the theoretical and literature values (shown in Table 4). For each substrate, the measured values were smaller than the calculated theoretical values but relatively close to the values reported in the literature. The YCH4 of PM was above the interval reported by the literature. However, it has been reported that different factors like age, sex, type of feeding, and separation processes of the manure can cause significant variations in YCH4 results [26][27][28][29].
The measured YCH4 of CWM was significantly different from the calculated theoretical value, though it was within the interval reported by the literature. Since the three batch-tests behaved similarly with a CV of 2.34%, this finding was attributed to cows being ruminants. Thus, the organic matter was partially degraded before the AD. Furthermore, the presence of microorganisms provided high VS values in the analysis of Both PM and CWM had reasonably small lag phases. This was attributed to the relatively large presence of nitrogen-associated compounds, which are common in these substrates [19]. Nevertheless, the ammonia concentrations were very different from the inhibition values reported by [25]. Thus, the lag phases were probably related to an adaptation phase rather than due to inhibition. SWW was the only readily degradable substrate because of the relatively large presence of nitrogen-free extracts. Consequently, no lag phase was detected for SWW. Instead, a sharp change in slope occurred after the second day, indicating slight diauxic behavior. Thus, it was assumed that a second substrate was consumed after the second day. However, no plateau phase was observed, and the MYC was considered monotonic.
The majority of the batch-tests had small lag phases. However, none of the lag phases were long enough to cause inhibition. Similarly, the Y CH4 values were compared to the theoretical and literature values (shown in Table 4). For each substrate, the measured values were smaller than the calculated theoretical values but relatively close to the values reported in the literature. The Y CH4 of PM was above the interval reported by the literature. However, it has been reported that different factors like age, sex, type of feeding, and separation processes of the manure can cause significant variations in Y CH4 results [26][27][28][29].
The measured Y CH4 of CWM was significantly different from the calculated theoretical value, though it was within the interval reported by the literature. Since the three batch-tests behaved similarly with a CV of 2.34%, this finding was attributed to cows being ruminants. Thus, the organic matter was partially degraded before the AD. Furthermore, the presence of microorganisms provided high VS values in the analysis of the macromolecules. The VS were not available for the production of biogas. However, VS are included in the calculation of the BMP theo in Equation (1). The f d was approximately 85% in PM and SWW. This result was expected because the AD occurred in the liquid phase, and no fiber content was measured in the substrates. Thus, hydrolysis was probably not a rate-limiting step. Next, the batch-test data were processed for compatibility with the model fitting procedures.

Model Fitting
The three models to describe MYCs were fitted to the measured data from the 20 batchtests, which met the three data pre-processing criteria. Table 5 shows the average RMSE results of the fit for each substrate and mixture. The modified Gompertz model had the best fit in all cases. According to Koch et al. [30], the modified Gompertz model is a better fit when a lag phase is necessary to describe a curve. However, the lag phases detected were rather small. The fit to the first-order one-step model was not acceptable, based on the RMSE value of 10 L CH4 /kg VS . This poor fit is due to small lag phases in all MYCs and hydrolysis being unlikely the rate-limiting step. A similar explanation applies to the poor fit from the first-order two-step model. However, the first-order two-step model accounts for a second degradation constant of VFAs. Thus, the first-order two-step model fit to the MYC was better than the first-order one-step model fit. However, the first-order two-step model did not meet the RMSE criteria for four of the seven substrate or mixture cases. Moreover, the values of the constants were predominantly the same, indicating that the process had only one rate-limiting step.
None of the three models adequately fit the MYCs of the SWW. Thus, a two-substrate model would provide a better fit, despite lacking a visually observed plateau phase in the SWW curve. However, for consistency in the mixture experimental design, all substrates must be described using the same model. Table 6 summarizes the average values from the fitted modified Gompertz model along with the measured pH before and after the AD. The adjusted lag times were relatively short. The lag times from the samples containing SWW were approximately one day slower than those without it. This was due to better C/N in the mixtures from the supply of carbohydrates readily degraded by the SWW. It was also due to high nitrogen content in substrates, like manures. Also, substrates with an initial pH of ≥8.10 corresponded to a lag phase of at least two days. Substrates with initial pH values < 8.10 had a maximum lag phase of 1.27 days. This drop in performance was associated with a shift in the NH4+-NH3 equilibria. This was because a higher pH moves the equilibria to NH3 production. It was reported that ammonia inhibition is commonly found in protein-rich substrates like the digestions of the manures and their double mixtures [19]. The pH was always in the optimal recommended interval of 7.4 to 7.6 [31] at the end of the AD. This indicated that the manures provided adequate buffer capacity.
Continuous surveying of the operational parameters during the batch-tests was not feasible. However, the measured initial and final pH values combined with the continuous gas production indicate that the process was kept in stable operating conditions.
The adjusted BMP ∞ of PM and CWM were smaller than the measured BMP. This was assumed to be due to the cancellation of noise from the measured data by the model, associated with the 10% measurement uncertainty of the eudiometers.

Mixture Characterization
The mixtures were further analyzed based on the CI and KI parameters, the results are shown in Table 7. No positive interactions were found for either the Y CH4 or r CH4 of the double mixtures. However Ebner et al. [21] used CI to characterize nine double mixtures of manure and a second substrate. The CIs ranged from −32% to 21% for a mixture proportion of 70:30% w/w. The mixtures of manures and carbon-rich substrates showed statistically significant positive effects. This was attributed to the buffering of the VFAs by the manure when it is digested together with carbon-rich substrates, as explained in Mata-Alvarez et al. [32]. The mixture from CWM + SWW carries this assertion, since despite having both negative CI and KI, an increase of both the methane yield and the rate were observed when compared to the individual digestion of CWM. This suggests that the augmentation of the design could find optimal double mixtures PM + SWW and CWM + SWW by providing new data and improving the interpolating capacity of the model. For the triple mixture, both parameters showed a positive effect with a CI value of 32.26% and a KI value of 7.42%. The triple mixture was the only mixture with a pH close to neutral at the beginning and a constructive effect for both CI and KI, which in batch processes benefits the acidogenic microorganisms [31]. Also, Astals et al. [33] found that "mixing a carbohydrate and/or protein source to lipids is a feasible option to reduce long-chain fatty acids (LCFA) inhibition, mainly due to the dilution". Furthermore, he concluded that AcoD leads to an enhancement of the AD kinetics, but rarely to a methane yield increase. However, in the triple mixture, both are observed. Thus, the superior performance of the triple mixture was attributed to better macro and micronutrient balance.

Mixture Design
An extra batch-test was manually added to the twenty measured batch-tests. The added test resulted from the average between the two fitted PM curves. Consequently, all substrates had three tests in the mixture design. The response variables were Y CH4 and r CH4 and the input data were the parameters resulting from the fitted curves by the modified Gompertz model. Special cubic models were fit for each variable. Model selection was made using a stepwise procedure with a 90% confidence interval for the parameters. The goodness-of-fit statistics for the models are detailed in Table 8. Table 8. Goodness-of-fit statistics for the special cubic models by response variable. The statistics indicated very good goodness-of-fit. Thus, the model was deemed acceptable for prediction purposes. The predicted model equations for the response variables were as follows:

Variable
r CH 4 = 47.03 · PM + 18.63 · CM + 32.11 · SWW + 6.16 · PM · CM + 17.39 · PM · SWW + 13.74 · CM · SWW +373.08 · PM · CM · SWW (12) where PM: Piglet Manure CWM: Cow Manure SWW: Starch Wastewater The predicted model includes all of the mixtures that participated in the design since all terms were statistically significant in the equations. The response equation for Y CH4 indicated a very negative interaction between two double mixtures: (1) PM and SWW and (2) CWM and SWW. A slight negative interaction was observed between PM and CWM. Additionally, the strength of the negative effect in the double mixtures negatively correlated with the protein ratio in the mixture. The triple mixture interaction was very positive, as previously detected in the characterization. Furthermore, Pagés-Díaz et al. [12] found qualitatively similar results having the highest positive effects among triple, while Kashi et al. [24] found the best results in a mixture of four substrates, as well that the mixture was very sensitive to changes in their composition. Moreover, the interactions between the variables in the equation were consistent with the characterization of the mixtures by the KI and CI parameters. Therefore, it served as a practical validation of the model.
The equation that described r CH4 showed positive interactions for all terms. The positive interactions did not necessarily contradict the characterization by KI. Unlike CI, KI only compares the fastest component of the mixture and not the mixture's predicted rate from the combination of individual substrates. The finding that all rate equation influences are positive is promising for further development in this research area. This finding indicates that it is more likely that the kinetics interaction transfers to the continuous stage rather than to the yield [17,23]. Therefore, the equation for r CH4 provided a better description of the interactions than the equation for KI.
Based on the r CH4 equation, the weakest positive interaction occurred in the PM + CWM mixture. Also, the magnitude of the positive interaction between double mixtures positively correlated with the percentage of readily degradable carbohydrates. It should be noted that the model was calibrated based solely on the seven mixtures but interpolated for the entire VS fraction range of 0 to 1 for each substrate. Thus, the positive interactions in the double mixtures that were outside the range of measured ratios indicated a need to improve the model by recalibration with additional runs.
The multi-objective optimization was performed after the constructed models were validated by acceptable matching with the practical values obtained from the characterization of the curves.

Optimization
The optimization was conducted with three constraints to treat a significant fraction of the piglet manure and to increase the chance of success in a continuous long-term operation. The three constraints were as follows: 1.
Maximize the specific methane production rate; 2.
Bound the BMP with a minimum value of 450 L CH4 /kg VS ; 3.
Require a minimal fraction of 0.4 of volatile solids in the piglet manure.
The optimization goal was to detect the optimal region(s) where coupled strong positive kinetics interacted with high Y CH4 quantities. The constraint that the minimum VS fraction was at least 0.4 guaranteed that most of the PM was treated due to its low VS content in terms of FM. The CWM content was indirectly restricted by giving a higher weight in the optimization to the r CH4 rather than the Y CH4 . This was due to the CWM deceleration effect on the degradation. However, it was advantageous to limit the DM content in the subsequent high-rate continuous operation.
A single common optimal region was found. Contour plots display the optimal region in Figure 3.
Processes 2021, 9, x FOR PEER REVIEW 12 of 17 be noted that the model was calibrated based solely on the seven mixtures but interpolated for the entire VS fraction range of 0 to 1 for each substrate. Thus, the positive interactions in the double mixtures that were outside the range of measured ratios indicated a need to improve the model by recalibration with additional runs. The multi-objective optimization was performed after the constructed models were validated by acceptable matching with the practical values obtained from the characterization of the curves.

Optimization
The optimization was conducted with three constraints to treat a significant fraction of the piglet manure and to increase the chance of success in a continuous long-term operation. The three constraints were as follows: 1. Maximize the specific methane production rate; 2. Bound the BMP with a minimum value of 450 LCH4/kgVS; 3. Require a minimal fraction of 0.4 of volatile solids in the piglet manure.
The optimization goal was to detect the optimal region(s) where coupled strong positive kinetics interacted with high YCH4 quantities. The constraint that the minimum VS fraction was at least 0.4 guaranteed that most of the PM was treated due to its low VS content in terms of FM. The CWM content was indirectly restricted by giving a higher weight in the optimization to the rCH4 rather than the YCH4. This was due to the CWM deceleration effect on the degradation. However, it was advantageous to limit the DM content in the subsequent high-rate continuous operation.
A single common optimal region was found. Contour plots display the optimal region in Figure 3.  The optimal region in the contour plot of Y CH4 (Figure 3a) was bigger than the optimal rate region, but both were similarly located. Consequently, the yield region contained most of the rate region. However, the third imposed constraint bounded the optimal range to the section between the PM vertex and its inferior bound on 0.4 PM. This resulted in an optimal area that did not contain most of the high Y CH4 area. Nevertheless, the relatively large size of the optimal zone represents a practical advantage since it can provide some resilience to measurement errors in VS content and still be able to operate inside the optimal region.
The three best solutions found in the optimal desirability region are shown in Table 9. Table 9. Goodness-of-fit statistics for the special cubic models by response variable. The three optimal solutions were triple mixtures with high predicted values for both CI and KI. In the first two mixtures, at least 71% of the VS in the mixtures were formed by the manures. Consequently, large quantities of manure can be treated by AcoD. Also, due to the buffer capacity of the three mixtures and their low VS values, the mixtures provide an opportunity to reach a stable operation in continuous high-rate reactors. The mixtures are composed of macro and micronutrients that provide high yields. No significant differences were found in the response variables' values between the three mixtures. This provides practical flexibility for continuous operation since the substrate's availability is sometimes a limiting factor. [19]. Additionally, other constraints can be imposed to find optimal solutions for the double mixture, but alternative approaches were not investigated in this study.

Validation of the Models
Two mixtures were selected to validate the models, namely, PM + CWM and PM + CWM + SWW. This allowed evaluation of the effect of adding SWW to a base mixture of PM + CWM. The constraints were as follows: 1.
The mixture had to remain within the optimum zone.

2.
The PM composition had to be the same in both mixtures.
The validation used fresh substrates that were collected during a different time of year than the substrates used in the modeling. The mixtures selected after optimization are detailed in Table 10. One batch-test from the co-digestion of PM + CWM had to be eliminated to keep the CV value under 5%. Furthermore, the curves showed a smooth trend, and the lag phase was barely noticeable. This is a good sign for the practical application, despite affecting the fit to the modified Gompertz model. Moreover, as predicted by the model, the addition of SWW to the mixture of both manures produced a significant increase in the performance.
The comparison of the methane yield of the double and the triple mixture is shown in Figure 4.
was barely noticeable. This is a good sign for the practical application, despite affecting the fit to the modified Gompertz model. Moreover, as predicted by the model, the addition of SWW to the mixture of both manures produced a significant increase in the performance. The comparison of the methane yield of the double and the triple mixture is shown in Figure 4. The average MYCs were fitted by the modified Gompertz model. The fitted rates were superior to the predicted rates by the resulting mixtures' design model due to the lack of a lag phase. Thus, an increase in the rate was observed, and the modified Gompertz model was not the best possible fit. However, the modified Gompertz model was employed to establish a direct comparison between the predicted and the measured results from the validation. Still, this model was a good fit with an RMSE value of less than 10 LCH4/kgVS. The results of the predicted and measured values are detailed in Table  11. The predicted values were used as references in the relative error calculations. The relative errors in YCH4 were relatively small considering the multiple sources of error, including the age of inoculum, differences in the VS content between the substrates used for modeling, and human error in the preparation of the batch-tests. The average MYCs were fitted by the modified Gompertz model. The fitted rates were superior to the predicted rates by the resulting mixtures' design model due to the lack of a lag phase. Thus, an increase in the rate was observed, and the modified Gompertz model was not the best possible fit. However, the modified Gompertz model was employed to establish a direct comparison between the predicted and the measured results from the validation. Still, this model was a good fit with an RMSE value of less than 10 L CH4 /kg VS . The results of the predicted and measured values are detailed in Table 11. The predicted values were used as references in the relative error calculations. The relative errors in Y CH4 were relatively small considering the multiple sources of error, including the age of inoculum, differences in the VS content between the substrates used for modeling, and human error in the preparation of the batch-tests.
The relative error for r CH4 in the PM + CWM mixture was quite large, but it was encouraging because the measured value was superior than the predicted one. This occurred because the deceleration of the CWM decreased compared to the CWM used in the model fitting. The average rate of the single digestion of CWM was 22.40 to a previous value of 18.66 L CH4 /kg VS /d. Thus, it was concluded that the models were valid for practical application and can be used for prediction. However, smaller values of Y CH4 and larger values of r CH4 should be expected. These two mixtures are currently being continuously tested in two separate high-rate reactors.
The developed mixture design model provided a reliable prediction and description of the interaction of the substrates at a macro scale, although it does not present profound insights into the biochemical interactions of the substrates in the mixture. This issue can be overcome by the development of a more complex mechanistic model. However, it usually requires the estimation of several non-measurable parameters, which is time-consuming and often does not assure the same level of precision as the empirical models.